Accueil | Cast3M

* fichier : jetkei.dgibi
************************************************************************
************************************************************************
************************************************************************
*  jetkei.dgibi
*
*                 jet 2D axi monophasique incompressible
*
*                  pour fiche de validation du K-epsilon
*
*  Pierre Cornet , sept 97
*  jpm           , mars 06 : adaptation pour K-epsilon implicite (kei)
*
************************************************************************
*
 graph = faux;
*graph = vrai;
*opti trace 'PSC' ;
COMPLET = FAUX;
 
*--------------------------- Numérique ---------------------------------
Si COMPLET;
 
KMAIL=1;
zr=1./3.;
zr=1.;
NBIT=2000;
NBIT=10000;
NBIT=3000;
NBIT=10000;
CFL=5.;
CFL= 1;
NBIT=400;
NBIT=3000;
*NBIT=800;
*NBIT=3   ;
CFL=5.;
DISCR = MACRO;
KPRESS= CENTREP1;
DISCR = LINE;
KPRESS= CENTRE;
sinon;
 
KMAIL=1;
zr=2.;
NBIT=80;
CFL=15.;
DISCR = LINE;
KPRESS= CENTRE;
finsi;
 
*--------------------------- procédure ---------------------------------
*$$$$ OUTFLOW
'DEBPROC' OUTFLOW  ;
ARGU RX*TABLE ;
************************************************************************
* ZONE $Mt OPER 'OUTFLOW' $out Ro Un Muf INCO 'UN'
*
*         UN
*         MUF       viscosité dynamique effective
*                             (tourbillonnaire+moléculaire)
************************************************************************
rv=rx.'EQEX' ;
iarg=rx.'IARG' ;
*mess 'iarg='iarg;
$mt=rx.'DOMZ' ;
 
* Lecture du 1er Argument objet mmodel
Si(ega ('TYPE' rx.'ARG1') 'MMODEL  ');
$mfront=rx.'ARG1';
Dg=doma $mfront 'XXDIAGSI' ;
Sinon ;
erreur 2;
quitter outflow;
Finsi ;
 
* Lecture du 2ème Argument la densité
Si(ega ('TYPE' rx.'ARG2') 'MOT     ');
Ro = rv.inco.(rx.'ARG2');
Sinon ;
Ro = rx.'ARG2';
Finsi ;
 
* Lecture du 3ème Argument la vitesse
Si(ega ('TYPE' rx.'ARG3') 'MOT     ');
un = rv.inco.(rx.'ARG3');
Sinon ;
un = rx.'ARG3';
Finsi ;
 
* Lecture du 4ème Argument la viscosité
Si(ega ('TYPE' rx.'ARG4') 'MOT     ');
Muf = rv.inco.(rx.'ARG4');
Sinon ;
Muf = rx.'ARG4';
Finsi ;
 
**************** INITIALISATIONs ***********************
Si (Exist rx 'rxtm');
rxtm=rx.'rxtm';
lmx=rxtm.'lmx';
lm0=rxtm.'lm0';
lm=rxtm.'lm';
nc=rxtm.'nc';                                 
MSV = chai rxtm.'msv' ;
MI1 = rxtm.'MI1' ;
KPRES=rxtm.'KPRES';
Sinon;
 
rxtm=table 'KIZX';
rx.'rxtm'=rxtm;
Si (NON (EGA (dime rx.'LISTINCO') 1));
  mess ' Erreur dans la procédure Outflow !';
  mess ' Il doit y avoir une inconnue !';
  quitter outflow ;
Finsi;
MI1=EXTR rx.'LISTINCO' 1 ;
lmx=(extr (rv.inco.MI1) 'COMP');
lm0=(extr un 'COMP');
nc=dime lm0;
 
Si(nc > 1);
*mess 'Outflow Cas vectoriel ' MI1;
MSV=chai 'VECT';
lm=mots (chai 1 MI1) (chai 2 MI1);
Si (EGA (vale dime) 3);
lm=lm et (mots (chai 3 MI1));
Finsi;
kpr=rx.'KOPT'.'KPOIN';
Si(EGA KPR 2);KPRES='CENTRE';Finsi;
Si(EGA KPR 4);KPRES='CENTREP1';Finsi;
Si(EGA KPR 5);KPRES='MSOMMET';Finsi;
Sinon;
*mess 'Outflow Cas scalaire  ' MI1;
MSV=chai 'SCAL';
lm=mots MI1;
Finsi;
rxtm.'lmx'=lmx;
rxtm.'lm0'=lm0;
rxtm.'lm'=lm;
rxtm.'nc'=nc;
rxtm.'msv'=chai MSV;
rxtm.'MI1'=MI1;
rxtm.'KPRES'=KPRES;
 
rxtm.'LISTINCO'=MOTS MI1 ;
rxtm.'DOMZ'=$mfront;                  
rxtm.'KOPT'=rx.'KOPT';
rxtm.'IARG'=1;
rxtm.'EQEX'=rx.'EQEX';
rxtm.'NOMZONE'=' ';
rxtm.'TDOMZ'=0;
rxtm.'NOMOPER'=mot 'MDIA_T';
Finsi;
**************** Fin Initialisations *******************
 
 mfront=doma $mfront maillage ;
 nj=doma $mt 'NORMALEV';
*unj= vect nj 0.1 ux uy jaune ;
*trace unj mt;
 njf = redu nj mfront ;
*unj= vect njf 0.1 ux uy rouge ;
*trace unj mt TITR ' FRONTIERE';
 xn = rv.inco.MI1;
 ujf = redu un mfront ;
 us=psca ujf lm0 njf lm0;
 ius=masq us 'INFERIEUR' 0.;
 rv.inco.'us'=us;
 rv.inco.'ius'=ius;
 Si(Ega nc 1);
 us=2.*(redu xn mfront);
 Finsi;
 rxtm.'ARG1'=kcht $mfront scal sommet ((-1.)*Ro*us*ius);
 
 st1 mat1 = MDIA rxtm ;
*st1 mat1 = KOPS 'MATRIK' ;
 
 Dgi=Dg*ius;
 Si(nc > 1);
 puj=0.;
 Si(EXIST (rv.inco) 'PRESSION');
 pn=redu (elno $mt rv.inco.'PRESSION' KPRES) mfront;
 puj=Ro*pn*us;
 Finsi;
 gru = redu (Muf*(kops 'GRADS' (exco un 'UX') $mt)) mfront;
 grv = redu (Muf*(kops 'GRADS' (exco un 'UY') $mt)) mfront;
 mgunj = (nomc (Dgi*((psca gru lm0 njf lm0)-puj)) (extr lm 1))
       + (nomc (Dgi*((psca grv lm0 njf lm0)-puj)) (extr lm 2)) ;
 Si(Ega (vale dime) 3);
 mgunj = mgunj + (nomc (Dgi*((psca grv lm0 njf lm0)-puj)) (extr lm 3));
 Finsi;
 Sinon;
 grt = redu (Muf*(kops 'GRADS' xn $mt)) mfront ;
 mgunj = (nomc (Dgi*(psca grt lm0 njf lm0)) (extr lm 1)) ;
 Finsi;
 St1 = kcht $mfront MSV sommet comp lm ((-1.)*mgunj);
 
RESPRO St1 mat1 ;
FINPROC ;
*----------------------- fin procédure ---------------------------------
 
 
 
 
*--------------------------- maillage ----------------------------------
 
TITRE 'JET' ;
OPTI MODE AXIS ;
OPTI DIME 2 ELEM QUA8 ;
 
DJ = 2.e-2 ;   RJ = DJ/2. ;  RM = 50.*RJ ;                     
epsi=1.e-7;
 
* points :
DJJ=1.*DJ;
P00=0. 0.;       PJ0=RJ 0.;       PR0=RM 0.;       PJ5=RJ (50.*DJJ);
P02=0. (20.*DJJ); P03=0. (30.*DJJ); P04=0. (40.*DJJ); P05=0. (50.*DJJ);
PR2=RM (20.*DJJ); PR3=RM (30.*DJJ); PR4=RM (40.*DJJ); PR5=RM (50.*DJJ);
 
* segments verticaux :
Si (EGA KMAIL 1);
A02 = DROI P00 P02 dini (0.5*zr*DJ) dfin (1.*zr*DJ) ;
A23 = DROI P02 P03 dini (1.*zr*DJ) dfin (1.2*zr*DJ) ;
A34 = DROI P03 P04 dini (1.2*zr*DJ) dfin (1.3*zr*DJ) ;
A45 = DROI P04 P05 dini (1.3*zr*DJ) dfin (1.4*zr*DJ) ;
Finsi;
 
Si (EGA KMAIL 2);
A02 = DROI -40 P00 P02 dini (0.25*DJ) dfin (0.6*DJ) ;
A23 = DROI -12 P02 P03 dini (0.45*DJ) dfin (0.6*DJ) ;
A34 = DROI 8   P03 P04 ;
A45 = DROI -6  P04 P05 dini (0.6*DJ) dfin (1.*DJ) ;
Finsi;
 
Si (EGA KMAIL 3);
A02 = DROI -80 P00 P02 dini (0.1*DJ) dfin (0.6*DJ) ;
A23 = DROI -12 P02 P03 dini (0.45*DJ) dfin (0.6*DJ) ;
A34 = DROI 8   P03 P04 ;
A45 = DROI -6  P04 P05 dini (0.6*DJ) dfin (1.*DJ) ;
Finsi;
 
AXE = A02 ET A23 ET A34 ET A45 ;
BORD= AXE plus PR0 ;
 
* segments horizontaux :
Si (EGA KMAIL 1);
JET  = DROI 2  P00 PJ0 ;
BAS2 = DROI    PJ0 PR0 dini (zr*RJ/2.) dfin (5.*zr*RJ) ;
 BAS = JET ET BAS2 ;
Finsi;
 
Si (EGA KMAIL 2);
JET  = DROI 5 P00 PJ0 ;
BAS2 = DROI    PJ0 PR0 dini (RJ/5.) dfin (10.*RJ) ;
 BAS = JET ET BAS2 ;
Finsi;
 
Si (EGA KMAIL 3);
JET  = DROI 10 P00 PJ0 ;
BAS2 = DROI    PJ0 PR0 dini (RJ/10.) dfin (10.*RJ) ;
 BAS = JET ET BAS2 ;
Finsi;
 
HAUT = BAS plus P05;
elim (haut et axe et bord et bas) epsi;
 
* domaine total :
MT = DALL  BAS BORD (inve HAUT) (inve AXE) 'PLAN' ;
CMT = cont MT ;
 
Mmt  = chan QUAF MT  ;
Mjet = chan QUAF jet ;
Mbord= chan QUAF bord;
Mbas2= chan QUAF bas2;
Maxe = chan QUAF axe ;
Mhaut= chan QUAF haut;
Mout = chan QUAF (bord et haut);
 
elim (Mmt et Mjet et Mbord et Mbas2 et Maxe et Mhaut) epsi;
elim (Mmt et Mout) epsi;
 
$MT   = mode Mmt   'NAVIER_STOKES' DISCR;
$JET  = mode Mjet  'NAVIER_STOKES' DISCR;
$BORD = mode Mbord 'NAVIER_STOKES' DISCR;
$BAS2 = mode Mbas2 'NAVIER_STOKES' DISCR;
$AXE  = mode Maxe  'NAVIER_STOKES' DISCR;
$HAUT = mode Mhaut 'NAVIER_STOKES' DISCR;
$out  = mode Mout  'NAVIER_STOKES' DISCR;
 
Ml20d = Mhaut moins (P05 moins P02) coul verte;
Ml30d = Mhaut moins (P05 moins P03) coul verte;
Ml40d = Mhaut moins (P05 moins P04) coul verte;
 
elim (Mmt et Ml20d et Ml30d et Ml40d) epsi;
 
MT  = doma $MT  maillage ;
axe = doma $axe maillage ;
jet = doma $jet maillage ;
bas2= doma $bas2 maillage ;
bord= doma $bord maillage ;
out = doma $out maillage ;
haut= doma $haut maillage;
 
doma $mt 'IMPR';
 
Si Graph;
* On prend soin de vérifier que l'allongement n'est pas supérieur à 20
diamin=doma $mt 'DIAMIN';
diamax=doma $mt 'DIAMAX';
ksi=diamax*(inve diamin);
ksi=elno $mt ksi centre;
trace ksi mt (cont mt) TITR ' Facteur d allongemet du maillage';
trace MT;
Finsi;
 
*--------------------------- maillage ----------------------------------
 
 
 
*------------------------ donnees physiques ----------------------------
*
NUF = 1.5E-5 ;
REJ = 2.e5   ;
REJ = 100.   ;
*REJ = 30.    ;
REJ = 2.e4   ;
REJ = 3000.  ;
UJ0 = REJ*NUF/DJ ; mess 'vitesse d injection (m/s) =' UJ0;
DT=CFL*RJ/UJ0 ;
mess ' CFL inj = ' CFL ' DT=' DT;
tasser jet;
xjet=coor 1 jet;
Dgj=doma $jet 'XXDIAGSI';                        
a=(SOMT (Dgj*((RJ - xjet)**(1./7.))))*(1./PI/RJ/RJ);
UJ=UJ0/a*((RJ - xjet)**(1./7.));
*dess (evol chpo uj jet);
 
*opti donn 5;
 KJ  = (0.02*UJ0)**2.  ;     mess 'KJ =' KJ ;
 EJ  =(KJ**1.5)/DJ ; mess 'EJ =' EJ ;
NUTj = 0.09*KJ*KJ/EJ ; mess 'NUTJ =' nutj ;
KA  = KJ*1.E-5 ;
NUTAsnu=1.e-4;
cnu=0.09;
EA=cnu*ka*ka/nuf;
 
*-------------------------- equations ----------------------------------
*$out = $haut;
RV = EQEX  $MT 'DUMP' 'ITMA' NBIT
    'OPTI' 'EF' 'SUPG' 'IMPL' KPRESS
*   'ZONE' $MT OPER OUTFLOW $out 1. 'UN' 'MUF'  'INCO' 'KN'
*   'ZONE' $MT OPER OUTFLOW $out 1. 'UN' 'MUF'  'INCO' 'EN'
    'ZONE' $MT OPER OUTFLOW $out 1. 'UN' 'MUF'  'INCO' 'UN'
    'ZONE' $MT 'OPER' 'KEPSILON' 1. 'UN' NUF 'DT' 'INCO' 'KN' 'EN'
    'ZONE' $MT 'OPER' 'NS' 1. 'UN' 'MUF'  'INCO' 'UN';
RV = EQEX RV
    'OPTI' 'EFM1' 'CENTREE'
    'ZONE' $MT 'OPER' DFDT 1. 'UN'  'DT'  'INCO' 'UN'
    ;
 
    RV = EQEX RV
    'CLIM' 'UN' UIMP JET    0.      'UN' VIMP JET    UJ
           'UN' VIMP BAS2   0.      'UN' UIMP AXE    0.
           'KN' TIMP JET    KJ      'EN' TIMP JET    EJ
*                                   'UN' VIMP BORD   0.
*          'KN' TIMP Bord   0.      'EN' TIMP Bord   1.e-10
           'KN' TIMP BAS2   0.      'EN' TIMP BAS2   1.e-10;
     ;
*RV.'ALGO_KEPSILON'= MOTS 'RNG';
*RV.'ALGO_KEPSILON'= MOTS 'Ret';
 
RVP = EQEX 'OPTI' 'EF' KPRESS
 'ZONE' $MT 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES' ;
 
RV.'PROJ'= RVP ;
 
*------------------------ initialisations ------------------------------
 
RV.INCO = TABLE 'INCO' ;
RV.'INCO'.'UN'   = KCHT $MT VECT SOMMET (1.E-7 1.E-7) ;
RV.'INCO'.'PRES' = 'KCHT' $MT 'SCAL' KPRESS 0.;
RV.'INCO'.'KN'   = KCHT $MT SCAL SOMMET 1.E-7 ;
RV.'INCO'.'EN'   = KCHT $MT SCAL SOMMET 1.E-5 ;
RV.'INCO'.'NUT'  = KCHT $MT SCAL CENTRE 1.E-6 ;
RV.'INCO'.'MUF'  = KCHT $MT SCAL sommet nuf   ;
RV.'INCO'.'DT'   = DT ;
 
*------------------------ historiques ----------------------------------
 
P11 = MT POIN 'PROC' ((25.*RJ) (10.*DJ));
P12 = MT POIN 'PROC' ((25.*RJ) (20.*DJ));
P14 = MT POIN 'PROC' ((25.*RJ) (40.*DJ));
P15 = MT POIN 'PROC' ((25.*RJ) (50.*DJ));
LH = P02 et P03 et P04 et P05 et P11 et P12 et P14 et P15 ;
HIS = KHIS 'UN' 1 LH  'UN' 2 LH  'KN' LH  'EN' LH ;
RV.'HIST' = HIS ;
 
*------------------------ resolution -----------------------------------
 
 lh=   (POIN MT PROC (0.05 0.01)) et (POIN MT PROC (0.05 0.05))
    et (POIN MT PROC (0.05 0.1))  et (POIN MT PROC (0.05 0.15))
    et (POIN MT PROC (0.05 0.2)) et (POIN MT PROC (0.05 0.7))
    et (POIN MT PROC (0.05 0.3)) et (POIN MT PROC (0.05 0.6))
    et (POIN MT PROC (0.05 0.4)) et (POIN MT PROC (0.05 0.5));
 
  his=khis 'UN' 1 lh
           'UN' 2 lh
           'KN'   lh
           'EN'   lh;
      his.'KFIH'=1;
  rv.'HIST'=his;
 
 EXEC rv;
 
*opti sauv 'JETKEI1';
*sauv;
*------------------------ post-traitement ------------------------------
UN = (RV.'INCO'.'UN');
uaxe = ((exco UN 'UY') redu axe)*(1./UJ0);
evuaxe= evol chpo uaxe axe ;
 
Si graph
 trace ((cont mt) et Ml20d et Ml30d et Ml40d);
 
 dessin his.'TABD' his.'KN' TITR ' historiques:  k';
 dessin his.'TABD' his.'EN' TITR ' historiques:  epsilon';
 dessin his.'TABD' his.'1UN' TITR ' historiques:  ux';
 dessin his.'TABD' his.'2UN' TITR ' historiques:  uy';
 
UNV = VECT UN (0.1/UJ0) UX UY JAUNE;
TRAC UNV CMT TITRE 'VITESSES AIR'  ;
 
dess evuaxe TITR 'Vitesse sur l axe';
 
z= (coor 2 axe) + 1.e-3;
vaxe= 5.8 * Dj * (inve z);
evax = evol chpo vaxe axe ;                                             
dess (evax et evuaxe) ybor 0.1 1.2 TITR ' Vitesses sur l axe';
 
 u20d  = (exco un 'UY') redu Ml20d;
 um20d = maxi u20d;
 u20d = evol chpo (u20d*(1./um20d)) ml20d ;
 
 u30d  = (exco un 'UY') redu Ml30d;
 um30d = maxi u30d;
 u30d = evol chpo (u30d*(1./um30d)) ml30d ;
 
 u40d  = (exco un 'UY') redu Ml40d;
 um40d = maxi u40d;
 u40d = evol chpo (u40d*(1./um40d)) ml40d ;
 
 z20d = maxi (coor 2 ml20d);
 z30d = maxi (coor 2 ml30d);
 z40d = maxi (coor 2 ml40d);
 r20d = (extr u20d 'ABSC')*(1./z20d);
 r30d = (extr u30d 'ABSC')*(1./z30d);
 r40d = (extr u40d 'ABSC')*(1./z40d);
 
 u20d = evol manu r20d (extr u20d 'ORDO');
 u30d = evol manu r30d (extr u30d 'ORDO');
 u40d = evol manu r40d (extr u40d 'ORDO');
 
 dess (u20d et u30d et u40d) XBOR 0. 0.4
 TITR ' Profil radiaux de vitesse';
 
 
nutsnu = rv.inco.'MUF' * (1./NUF) ;
trace nutsnu mt (cont mt) TITR 'Muf / Mu' ;
trace (log nutsnu) mt (cont mt) TITR ' Log (Muf / Mu)' ;
 
trace rv.inco.'KN' mt (cont mt) TITR ' Kn ';
trace (log rv.inco.'KN') mt (cont mt) TITR ' Log Kn ';
trace rv.inco.'EN' mt (cont mt) TITR ' En ';
trace (log rv.inco.'EN') mt (cont mt) TITR ' Log En ';
Finsi;
*-------------------- fin post-traitement ------------------------------
 
 Si (NON COMPLET);
 uaxref=prog                                            
   1.8650       1.8896       1.9826       1.8759       1.4946
   1.1368      0.91643      0.74427      0.60911      0.51210
  0.43502      0.37462      0.32726      0.28785      0.25592
  0.23163      0.21113      0.19329      0.17719      0.16177
  0.14722      0.13685      0.13161      0.12988      0.12893
  0.12679      0.12180      0.11305;
zaxe = extr evuaxe 'ABSC';
list (extr evuaxe 'ORDO');
evuaxer=evol manu zaxe uaxref ;
m=(INTG evuaxer 'ABSO') ;
delt2= (INTG (evuaxer - evuaxe) 'ABSO')/m;
mess ' deltaL2 : m=' m ' delt2=' delt2;
 Si (delt2 > 1.e-3);ERREUR 5;Finsi;
 FINSI;
 
FIN ;