Accueil | Cast3M

* fichier : jetplankei.dgibi
************************************************************************
************************************************************************
************************************************************************
*  jetplankei.dgibi
*
*                 jet 2D Plan monophasique incompressible
*
*                  pour fiche de validation du K-epsilon
*
*  Pierre Cornet , sept 97
*  jpm           , mars 06 : adaptation pour K-epsilon implicite (kei)
*
************************************************************************
*
COMPLET = FAUX;
*--------------------------- Numérique ---------------------------------
Si COMPLET;
NBIT=400;
DT = 0.002 ;
DISCR = MACRO;
KPRESS= CENTREP1;
sinon;
NBIT=80;
DT = 0.01  ;
DISCR = LINE;
KPRESS= CENTRE;
finsi;
 
graph = faux;
*graph = vrai;
*opti trace 'PSC' ;
*----------------------- fin Numérique ---------------------------------
 
*--------------------------- 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 DIME 2 ELEM QUA8 ;
*
DJ = 2.e-2 ;   RJ = DJ/2. ;  RM = 150.*RJ ;
DJ = 2.e-2 ;   RJ = DJ/2. ;  RM =  50.*RJ ;
 
* points :
P00=0. 0.;       PJ0=RJ 0.;       PR0=RM 0.;       PJ5=RJ (50.*DJ);
P02=0. (20.*DJ); P03=0. (30.*DJ); P04=0. (40.*DJ); P05=0. (50.*DJ);
PR2=RM (20.*DJ); PR3=RM (30.*DJ); PR4=RM (40.*DJ); PR5=RM (50.*DJ);
 
* segments verticaux :
A02 = DROI -20 P00 P02 dini (0.5*DJ) dfin (1.25*DJ) ;
A23 = DROI -6  P02 P03 dini (0.9*DJ) dfin (1.25*DJ) ;
A34 = DROI 4   P03 P04 ;
A45 = DROI -3  P04 P05 dini (1.25*DJ) dfin (2.*DJ) ;
 
B02 = DROI -20 PR0 PR2 dini (0.5*DJ) dfin (1.25*DJ) ;
B23 = DROI -6  PR2 PR3 dini (0.9*DJ) dfin (1.25*DJ) ;
B34 = DROI 4   PR3 PR4 ;
B45 = DROI -3  PR4 PR5 dini (1.25*DJ) dfin (2.*DJ) ;
 
AXE = A02 ET A23 ET A34 ET A45 ;
BORD= B02 ET B23 ET B34 ET B45 ;
 
* segments horizontaux :
JET  = DROI 2  P00 PJ0 ;
BAS2 = DROI    PJ0 PR0 dini (RJ/2.) dfin (10.*RJ) ; BAS = JET ET BAS2 ;
 
HAU1 = DROI 2  P05 PJ5 ;
HAU2 = DROI    PJ5 PR5 dini (RJ/2.) dfin (10.*RJ) ; HAUT=HAU1 ET HAU2 ;
 
* 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;
 
elim (Mmt et Mjet et Mbord et Mbas2 et Maxe et Mhaut) 1.e-5;
 
$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;
 
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) 1.e-5;
 
MT  = doma $MT  maillage ;
axe = doma $axe maillage ;
jet = doma $jet maillage ;
bord= doma $bord maillage ;                                             
bas2= doma $bas2 maillage ;
haut= doma $haut maillage ;
 
Si Graph;trace MT;Finsi;
*------------------------ donnees physiques ----------------------------
*
NUF = 1.5E-5 ;
REJ = 2.e4   ;
UJ  = REJ*NUF/DJ ; mess 'vitesse d injection (m/s) =' UJ ;
 
* KJ  = 0.05*UJ*UJ   ;     mess 'KJ =' KJ ;
* EJ  = 0.02*(UJ**3.)/DJ ; mess 'EJ =' EJ ;
KJ = 1.E-3 ;
EJ = 6.E-3 ;
NUTj = 0.09*KJ*KJ/EJ ; mess 'NUTJ =' nutj ;
KA  = 1.E-7  ;
EA  = 1.E-5  ;
L0  = 25.*DJ  ;
 
*-------------------------- equations ----------------------------------
 
RV = EQEX  $MT 'DUMP' 'ITMA' NBIT
    'OPTI' 'EF' 'SUPG' 'IMPL'
*   'ZONE' $MT OPER OUTFLOW $haut 1. 'UN' 'MUF'  'INCO' 'KN'
*   'ZONE' $MT OPER OUTFLOW $haut 1. 'UN' 'MUF'  'INCO' 'EN'
    'ZONE' $MT OPER OUTFLOW $haut 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'
    '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.
                                    'UN' VIMP BORD   0.
           'KN' TIMP JET    KJ      'EN' TIMP JET    EJ
           'KN' TIMP BORD   0.      'EN' TIMP BORD   EA  ;
 
*RV.'ALGO_KEPSILON'= MOTS 'RNG';
 
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'.'MUF'  = KCHT $MT SCAL SOMMET 1.E-6 ;
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;
 
UN = (RV.'INCO'.'UN');
uaxe = ((exco UN 'UY') redu axe)*(1./UJ);
evuaxe= evol chpo uaxe axe ;
 
*------------------------ post-traitement ------------------------------
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 5.E-3 UX UY  ;
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 / Nu' ;
 
trace rv.inco.'KN' mt (cont mt) TITR ' Kn ';
trace rv.inco.'EN' mt (cont mt) TITR ' En ';
Finsi;
*-------------------- fin post-traitement ------------------------------
 
 Si (NON COMPLET);
 uaxref=prog                                                   
  1.0000       1.0250      0.97229      1.0009      0.96387  
  0.99073      0.97449     1.0046       1.0016      1.0228    
  1.0149       1.0024      0.96638      0.93515     0.89487  
  0.85870      0.81994     0.78738      0.75198     0.72043    
  0.69096      0.66619     0.64104      0.61681     0.59295
  0.57226      0.55194     0.52699      0.50559     0.48760    
  0.47019      0.44757     0.42364      0.40301  ;
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 > 8.e-3);ERREUR 5;Finsi;
 FINSI;
 
 
 
 
 
FIN ;