Télécharger bc30.dgibi

Retour à la liste

* fichier :  bc30.dgibi
************************************************************************
************************************************************************
*$$$ BC30
*X   BC30 (Jeux de donnees)
 
**
**         BC30
**
**   CANAL CHAUFFE  INCLINE  ANGL = 30.  LONGUEUR H0=10 LARGEUR 1.
**   K-E  PLAN --- RE=2.E+5                                                     
**
**   teste NSKE(SUPG) TSCAL(PSI) LAPN FPU ECHI et VNIMP
**   Formulation EFM1 algorithme semi implicite
**   teste aussi les méthodes d'inversion de la pression :
**   Choleski et gradient conjugué.
**
**
**   --- 26 JANVIER  1987 ---  modifie  --- 26 MAI 1993      ---
**                             modifie  --- 17 DECEMBRE 1997 ---
**
**    A noter :
**
**   1/ la manière dont sont donnés les différents coefficients
**      du fait que df/dt a comme coefficient 1. pour toutes les
**      équations.
**   2/ La manière dont est réalisé le couplage : construction
**      du maillage de raccord (uniquement valable en 2D)
**      et procédure CALCUL .
**   3/ Le bilan thermique dans ce cas est assez mal vérifié .
**      Les causes : maillage grossier, flux thermique tres faible
**      par rapport au debit enthalpique, terme d'accumulation
**      d teta / dt non pris en compte dans le bilan.
**
**
 
GRAPH = 'N'                   ;
COMPLET = FAUX ;
 
NITMA=200;
 
SI (NON COMPLET) ;
err1=5.E-4;
NITMA = 40 ;
FINSI ;
 
option dime 2 elem tri6 ;
option dime 2 elem qua8 ;
p1=0  0 ;
p2=1. 0.;
 
nbpe=7  ; nbpp=10 ;
*nbpe=40 ; nbpp=80 ;
*nbpe=50 ; nbpp=150 ;
tole=1.e-5 ;
 H0=10.;
 
entree= p1 d nbpe p2 ;
*entree= p1 d dini 0.05 dfin 0.2  p2 ;
c1= 5. 0. ;
 
 q1=p1 tour c1 -90. ;
 q2=p2 tour c1 -90. ;
 pp1=p1 c c1 q1 nbpp ;
 pp2=p2 c c1 q2 nbpp ;
 
 q1=p1 plus (0 H0) ;
 q2=p2 plus (0 H0) ;
 pp1=p1 d    q1 nbpp ;
 eps=0.1  ;
 pri=pp1 plus ((-1.*eps) 0) ;
 pp2=p2 d    q2 nbpp ;
 
 sortie=entree tour c1 -90 ;
 sortie=entree plus (0 H0) ;
 pp1=inve pp1 ;
 sortie=inve sortie;
 elim tole   (sortie et pp1 et pp2 );
 cont=entree et pp2 et sortie et pp1 ;
*bell=surf cont ;
 bell=daller entree pp2 sortie pp1 ;
 paroi=pri trans 4 (-0.1 0) ;
 pr= cote 3 paroi ;
 entref=entree elem seg3 (lect 2 pas 1 6) ;
 
 bell=bell tour c1 -20. ;
 entree=entree tour c1 -20. ;
 entref=entref tour c1 -20. ;
 sortie=sortie tour c1 -20. ;
 pp1=pp1 tour c1 -20. ;
 pp2=pp2 tour c1 -20. ;
 
 paroi=paroi tour c1 -20. ;
 paroi=paroi coul rouge ;
 pr= pr tour c1 -20. ;
 pri= pri tour c1 -20. ;
 MT=bell et paroi ;
pp1=inve pp1 ;
 
elim (pp1 et pp2 et entree et entref et sortie et bell) tole   ;
elim ( pri et pr et paroi ) tole   ;
 
 titred= 'PLAQUES // COUDEES : TEST  TRIO-EF  NSKE FPU VNIMP' ;
 titre= 'PLAQUES // COUDEES : TEST  TRIO-EF  NSKE FPU VNIMP' ;
  Mmt    = CHAN mt     QUAF ;
  Mbell  = CHAN bell   QUAF ;
  Mparoi = CHAN paroi  QUAF ;
  Mentree= CHAN entree QUAF ;
  Mentref= CHAN entref QUAF ;
  Msortie= CHAN sortie QUAF ;
  Mpri   = CHAN pri    QUAF ;
  Mpp1   = CHAN pp1    QUAF ;
  Mpp2   = CHAN pp2    QUAF ;
 
  elim (Mmt et Mbell et Mparoi et Mentree et Mentref et Msortie
     et Mpri et Mpp1 et Mpp2) tole  ;
 
 $mt   = MODE MMT 'NAVIER_STOKES' 'MACRO' ;
 $bell = MODE Mbell 'NAVIER_STOKES' 'MACRO' ;
 $paroi= MODE Mparoi 'NAVIER_STOKES' 'MACRO' ;
  paroi = doma $paroi 'MAILLAGE' ;
 DOMA $mt 'IMPR' ; DOMA $bell 'IMPR' ; DOMA $paroi 'IMPR' ;
 
 $pri  = MODE Mpri 'NAVIER_STOKES' 'MACRO' ;
 $pp1  = MODE Mpp1 'NAVIER_STOKES' 'MACRO' ;
 $pp2  = MODE Mpp2 'NAVIER_STOKES' 'MACRO' ;
 $entree = MODE Mentree 'NAVIER_STOKES' 'MACRO' ;
 $entref = MODE Mentref 'NAVIER_STOKES' 'MACRO' ;
 $sortie = MODE Msortie 'NAVIER_STOKES' 'MACRO' ;
 
 pp1c=pp1 elem (lect 2 pas 1 (nbel pp1));
 pric=pri elem (lect 2 pas 1 (nbel pri));
 
 Mpp1c   = CHAN pp1c   QUAF ;
 Mpric   = CHAN pric   QUAF ;
 elim (Mpric et Mpp1c et Mmt ) tole ;
 
 $pp1c=MODE Mpp1c 'NAVIER_STOKES' 'MACRO' ;
 $pric=MODE Mpric 'NAVIER_STOKES' 'MACRO' ;
 pp1l=doma $pp1c MAILLAGE ;
 pril=doma $pric MAILLAGE ;
 
 rfp=racc eps pp1l pril ;
 rfpl= rfp chan ligne ;
 rfp= diff rfpl (pp1l et pril) ;
 rfpi= inve rfp ;
*************************  PROCEDURE CALCUL *****************
 DEBP CALCUL ;
 ARGU RX*TABLE ;
 rv=rx.'EQEX' ;
 
 tn=rv.'INCO'.'TN' ;
 un=rv.'INCO'.'UN' ;
 
 TF=redu tn pp1c ;
 rv.'INCO'.'TF'= kcht $pric scal sommet  (tf kpro rfp) ;
 TP=redu tn pric ;
 rv.'INCO'.'TP'= kcht $pp1c  scal sommet (tp kpro rfpi) ;
 rv.'INCO'.'UF'= kcht $bell vect sommet (redu un bell) ;
 
 as2 ama1  = 'KOPS' 'MATRIK' ;
 finproc as2 ama1 ;
 
*************************  PROCEDURE CALCUL *****************
 
*option donn 5 ;
NU =7.3E-7 ;
LAMBDAP=15.2 ;
ROCPP=3.7E6 ;
ALFP=LAMBDAP/ROCPP ;
*ALFP=4.06E-6 ;
LAMBDAE=0.62 ;
ROCPE=4.15E6 ;
ALFE=LAMBDAE/ROCPE ;
*ALFE=1.5E-7 ;
mess ' ALFE=' ALFE ' ALFP=' ALFP ;
YP =5.E-2 ;   
R0 =0.5  ;
TPSC =1500.;
TREF=0. ;
ZERO =1.E-20 ;
gb = 0. 0. ;
H=1.E4 ;
HF=H/ROCPE ;
HP=H/ROCPP ;
HP=H/ROCPE ;
CMU  = 0.09 ;
UREF = 4. ;
KE = (0.05 * UREF)*(0.05 * UREF) ;
EE = (KE**1.5) / R0 ;
nut = 0.09 * KE * KE / EE ;
  mess ' NUT entree ' nut ;
  UET1 = KCHT $PP1 SCAL CENTRE  5.E-2 ;
  UET2 = KCHT $PP2 SCAL CENTRE  5.E-2 ;
  NUT  = KCHT $BELL SCAL CENTRE nut ;
 
  RV= EQEX        'DUMP' 'ITMA' NITMA 'ALFA' 0.9 OPTI 'SUPG'
 'ZONE'  $MT      'OPER'  'CALCUL'
 'ZONE'  $BELL    'OPER'  'NSKE' NU NUT         'INCO' 'UN' 'KN' 'EN'
 'ZONE'  $PP1     'OPER'  'FPU'  NU UET1 YP     'INCO' 'UN' 'KN' 'EN'
 'ZONE'  $PP2     'OPER'  'FPU'  NU UET2 YP     'INCO' 'UN' 'KN' 'EN'
  ;
  RV=EQEX RV OPTI EFM1 'PSI'
 'ZONE'  $BELL    'OPER'  'TSCA' ALFE 'UF' 0. NUT 0.7 'INCO' 'TN'
 'ZONE'  $PAROI   'OPER'  'LAPN' ALFP           'INCO' 'TN'
 'ZONE'  $PRIC    'OPER'  'ECHI' HP   'TF'      'INCO' 'TN'
 'ZONE'  $PP1C    'OPER'  'ECHI' HF   'TP'      'INCO' 'TN'
  ;
  RV=EQEX RV OPTI EFM1 'CENTREE'
  'ZONE' $BELL    'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN'
  'ZONE' $MT      'OPER' 'DFDT' 1. 'TN' 'DELTAT' 'INCO' 'TN'
  'ZONE' $BELL    'OPER' 'DFDT' 1. 'KN' 'DELTAT' 'INCO' 'KN'
  'ZONE' $BELL    'OPER' 'DFDT' 1. 'EN' 'DELTAT' 'INCO' 'EN'
  ;
  RV=EQEX RV CLIM
  'KN' TIMP ENTREE  KE   'EN' TIMP ENTREE EE
  'TN' TIMP PR 100.   'TN' TIMP ENTREE  20.
  ;
  RV.INCO=TABLE 'INCO' ;
  UF = kcht $BELL VECT SOMMET (UREF*(0.5 0.9));
  RV.INCO.'UN'=kcht $BELL   VECT SOMMET (0. 0.) UF ;
  RV.INCO.'KN'=kcht $BELL   SCAL SOMMET   KE   ;
  RV.INCO.'EN'=kcht $BELL   SCAL SOMMET   EE   ;
  tp=kcht $paroi scal sommet 100. ;
  RV.INCO.'TN'=kcht $MT   SCAL SOMMET   20. tp  ;
 
* Suivi temporel de quelques points
lh= elem (doma $bell sommet) 'POI1' (lect 1 pas 20 101) ;
lc1= elem (doma $pp1 sommet) 'POI1' (lect 4 pas 4 20) ;
lc2= elem (doma $pri sommet) 'POI1' (lect 4 pas 4 20) ;
 lc=lc1 et lc2 ;
       his = khis 'UN' 1 lh
                  'UN' 2 lc1
                  'KN'   lh
                  'EN'   lh
                  'TN'   lc1 ;
   rv.'HIST'=his ;
 
 
 
 
DEBPROC TEST ityp*entier ksto*entier rv*table ;
 
*  1 : Choleski  (valeur par defaut)
 
**** Le reste est non maintenu et ne fonctionne pas sur machine 64bytes**
*  2 : gradient conjugue sans preconditionnement
*  3 : gradient conjugue avec preconditionnement diagonal
*  4 : gradient conjugue avec preconditionnement par Choleski incomplet
*  7 : gradient conjugue avec preconditionnement diagonal (variante)
*
*Pour 2,3,4, la matrice de pression est calculee et stockee une fois
*pour toutes. Pour 7, elle n'est jamais stockee explicitement et permet
*de ce fait d'atteindre les plus grosses tailles de problemes.
*
*ityp=1 ;
* obtenir ityp*entier ;
*
*Pour 2,3,4, on propose deux types de stockage : morse et compresse. Le
*second type n'est interessant que pour les machines vectorielles.
* ksto=0 morse    ksto=1 compresse
* ksto=0 Valeur par defaut
*
*ksto=0 ;
*obtenir ksto*entier ;
 
  RVP = EQPR $BELL  KTYPI ITYP 'ZONE' $BELL  'OPER' 'PRESSION' 0.
 'ZONE' $PP1    'OPER' 'VNIMP' 0.  'ZONE' $PP2    'OPER' 'VNIMP' 0.
 'ZONE' $ENTREE 'OPER' 'VNIMP' UREF 'ZONE' $ENTREF 'OPER' 'VTIMP' 0.;
 
 si (non (ega ityp 1));
 rvp.METHODE.KSTOCK=ksto ;
 rvp.METHODE.JTEPS=5.E-4;
 rvp.METHODE.'IPAT'=1;
 finsi ;
 
  RV.'PRESSION' =RVP ;
 
  EXEC RV ;
 
  mess ' Fin PROC ' ;
FINPROC RV ;
 
DEBPROC POST RV*TABLE GRAPH*MOT ;
 
 un=rv.'INCO'.'UN' ;
 tn=rv.'INCO'.'TN' ;
 
 qe=dbit un $entree ;
 qs=dbit un $sortie ;
 dq=(abs qs) - (abs qe) ;
 dqr=dq / (abs qe) ;
 mess ' Erreur sur la masse ' dqr ;
 ER1= abs (abs(dqr) - 1.E-10) ;
 
 srti=inve (doma $sortie MAILLAGE) ;
 lcu = extr un 'COMP';
 mu2=un lcu 'PSCA' un lcu;
 mu= mu2 ** 0.5 ;
 evolV = EVOL 'CHPO' mu SCAL srti  ;
 evy=extr evolV 'ORDO' ;
 list evy ;
 ley=prog 3.4573 3.8455 4.0470 4.0806 4.0848 4.0855 4.0853 4.0852
          4.0853 4.0855 4.0848 4.0806 4.0470 3.8455 3.4573 ;
 
 
 ER2=SOMM( (abs (ley - evy)) / 4. ) ;
 mess ' Ecart sur profil de vitesse en sortie ' er2 ;
 
 evolT = EVOL 'CHPO' tn SCAL srti  ;
 evt=extr evolT 'ORDO' ;
 list evt ;
 let=prog 20.985 20.501 20.179 20.078 20.045 20.026 20.015 20.008
          20.005 20.002 20.001 20.001 20.000 20.000 20.000  ;
 
 
 ER3=SOMM( (abs (let - evt)) / 20. ) ;
 mess ' Ecart sur profil de température en sortie ' er3 ;
 
 D=doma $pp1c 'XXDIAGSI' ;
 TF= kcht $pp1c scal sommet  (redu tn pp1c) ;
 TP=rv.inco.'TP' ;
 DTFP=kops tp '-' tf ;
 FP=somt ( hf* (kops d '*' (kops tp '-' tf) ) ) ;
 mess 'FP=' fp ;
 uteta= kops ('KCHT' $bell 'SCAL' 'SOMMET' tn) '*' un ;
 ute=dbit uteta $entree ;
 uts=dbit uteta $sortie ;
 dut=(abs uts) - (abs ute) - fp ;
 dutr=dut/FP ;
 mess (' Bilan : erreur relative = ') dutr ;
 ER4= abs (abs(dutr) - .58735) ;
 
 si (non COMPLET ) ;
 si ( ER1 > 1.E-9 ) ; erreur 5 ; finsi ;
 si ( ER2 > err1 ) ; erreur 5 ; finsi ;
 si ( ER3 > err1 ) ; erreur 5 ; finsi ;
 si ( ER4 > err1 ) ; erreur 5 ; finsi ;
 finsi ;
 
 SI ('EGA' GRAPH 'O' );
 
 TAB1=TABLE;
 TAB1.'TITRE'=TABLE ;
 TAB1 . 1 ='MARQ      REGU '     ;
 TAB1.'TITRE' . 1 = mot 'Module de U '    ;
 DESS evolV 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE  TAB1 ;
 
 TAB1.'TITRE' . 1 = mot 'Temperature sortie ' ;
 DESS evolT 'TITX' 'R (m)' 'TITY' 'T °C' LEGE  TAB1 ;
 
 
 ung1=vect un 0.1 ux uy jaune;
 pn=elno (kcht $bell scal centre (rv.'PRESSION'.pression)) $bell ;
 
 his=rv.'HIST' ;
 dessin his.'TABD' his.'1UN' ;
 dessin his.'TABD' his.'2UN' ;
 dessin his.'TABD' his.'TN' ;
 dessin his.'TABD' his.'KN' ;
 dessin his.'TABD' his.'EN' ;
 
 trace ung1 bell ;
 trace pn bell ;
 trace tn bell ;
 trace tn paroi;
 kn=rv.inco.'KN' ;
 en=rv.inco.'EN' ;
 nut=0.09*kn*kn/en ;
 trace nut bell ;
 trace kn bell ;
 trace en bell ;
 FINSI ;
 
FINPROC ;
 
******************** TESTS **************************
 
 ityp=1 ;
 ksto = 0 ;
 rv = test ityp ksto rv ;
 
* ityp=4 ;non maintenu, ne marche pas sur machine 64bytes
 ityp=1 ;
 ksto = 0 ;
 rv = test ityp ksto rv ;
 
* ityp=4 ;non maintenu, ne marche pas sur machine 64bytes
 ityp=1 ;
 ksto = 1 ;
 rv = test ityp ksto rv ;
 
* ityp=3 ;non maintenu, ne marche pas sur machine 64bytes
 ityp=1 ;
 ksto = 1 ;
 rv = test ityp ksto rv ;
 
* ityp=2 ;non maintenu, ne marche pas sur machine 64bytes
 ityp=1 ;
 ksto = 1 ;
 rv = test ityp ksto rv ;
 
 post rv graph ;
 FIN ;
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales