Accueil | Cast3M

* fichier :  jetaxi.dgibi
************************************************************************
************************************************************************
GRAPH = FAUX ;
COMPLET = FAUX ;
 
DISCR  = 'MACRO';
KPRESS = 'CENTRE';
BETA   =   1.;
 
************************************************************************
*  JETAXI.DGIBI : jet turbulent monophasique axisymétrique             *
************************************************************************
* FICHE DE VALIDATION DU K-EPSILON INCOMPRESSIBLE + FILTRE             *
* FORMULATION EF QUA8                                                  *
* P. CORNET SEMT/TTMF DECEMBRE 1998                                    *
** modifie le 15/06/2014 passage EQPR -> EQEX                          *
************************************************************************
 
SI ( COMPLET) ;
   NITER = 40000 ;
SINON ;
   NITER = 20 ;
FINSI ;
 
*--------------------------- maillage ----------------------------------
 
TITRE 'JET' ;
OPTI MODE AXIS ;
OPTI DIME 2 ELEM QUA8 ;
*
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' mt   'QUAF';
Mjet  = 'CHAN' jet  'QUAF';
Mbord = 'CHAN' bord 'QUAF';
Mbas2 = 'CHAN' bas2 'QUAF';
Mbas  = 'CHAN' bas  'QUAF';
Mhaut = 'CHAN' haut 'QUAF';
Maxe  = 'CHAN' axe  'QUAF';
 
ELIM (Mmt et Mjet et Mbord et Mbas2 et Mbas et Mhaut et Maxe) 1.e-5;
 
 
* tables domaine :
$MT   = 'MODE' MMT    'NAVIER_STOKES' DISCR;
$JET  = 'MODE' MJET   'NAVIER_STOKES' DISCR;
$BORD = 'MODE' MBORD  'NAVIER_STOKES' DISCR;
$BAS2 = 'MODE' MBAS2  'NAVIER_STOKES' DISCR;
$BAS  = 'MODE' MBAS   'NAVIER_STOKES' DISCR;
$HAUT = 'MODE' MHAUT  'NAVIER_STOKES' DISCR;
$AXE  = 'MODE' MAXE   'NAVIER_STOKES' DISCR;
DOMA $MT 'IMPR'     ;
MT   = 'DOMA' $MT   'MAILLAGE' ;
jet  = 'DOMA' $jet  'MAILLAGE' ;
bord = 'DOMA' $bord 'MAILLAGE' ;
bas2 = 'DOMA' $bas2 'MAILLAGE' ;
bas  = 'DOMA' $bas  'MAILLAGE' ;
haut = 'DOMA' $haut 'MAILLAGE' ;
axe  = 'DOMA' $axe  'MAILLAGE' ;
 
 
*------------------------ donnees physiques ----------------------------
*
NUF = 1.5E-5 ;
REJ = 1.e4   ;
UJ  = REJ*NUF/DJ ;
 
KJ = 1.E-3 ;
EJ = 6.E-3 ;
NUTj = 0.09*KJ*KJ/EJ ;
KA  = 1.E-7  ;
EA  = 1.E-5  ;
L0  = 25.*DJ  ;
 
*-------------------------- equations ----------------------------------
 
RV = EQEX  $MT 'DUMP' 'ITMA' NITER  'ALFA' 0.9
    'ZONE' $MT 'OPER' 'FILTREKE'  UJ  L0     'INCO' 'KN' 'EN'
    'ZONE' $MT 'OPER' 'NSKE'      NUF 'NUT'  'INCO' 'UN' 'KN' 'EN'
    '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
           'KN' TIMP BORD   KA      'EN' TIMP BORD   EA  ;
 
  RV=EQEX RV OPTI EFM1 'CENTREE'
  'ZONE' $MT    'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN'
  'ZONE' $MT    'OPER' 'DFDT' 1. 'KN' 'DELTAT' 'INCO' 'KN'
  'ZONE' $MT    'OPER' 'DFDT' 1. 'EN' 'DELTAT' 'INCO' 'EN'
  ;
 
RVP = EQEX 'OPTI' 'EF' KPRESS
 'ZONE' $MT    OPER  KBBT  -1. beta  INCO 'UN' 'PRES'
;
 
    rvp.'METHINV'.TYPINV=1 ;
    rvp.'METHINV'.IMPINV=0 ;
    rvp.'METHINV'.NITMAX=300;
    rvp.'METHINV'.PRECOND=3 ;
    rvp.'METHINV'.RESID  =1.e-8 ;
    rvp.'METHINV' . 'FCPRECT'=100 ;
    rvp.'METHINV' . 'FCPRECI'=100 ;
 
  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 ;
 
*------------------------ 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 ;
 
*------------------------ résolution -----------------------------------
 
EXEC RV ;
 
*------------------------ post-traitement ------------------------------
SI GRAPH ;
 
opti 'ISOV' 'LIGNE' ;
un = rv.INCO.'UN' ;
unn = un/uj ;
kn = rv.INCO.'KN' ;
en = rv.INCO.'EN' ;
nut= rv.INCO.'NUT' ;
nut= elno $mt nut  ;
nur= (1./nuf)*nut ;
pn = rv.pression.'PRESSION' ;
pe = elno $mt pn ;
 
*----------------------------- isovaleurs ------------------------------
 
UNV = VECT (RV.'INCO'.'UN') 5.E-3 UX UY vert ;
  TRAC UNV CMT TITRE 'VITESSES '  ;
  TRAC NUR  MT CMT TITRE 'NUT/NU'   ;
 
*------------------------- U sur l'axe ---------------------------------
 
uz  = exco uy un  uy ;
uzn = exco uy unn uy ;
uzne = evol chpo uzn uy axe ;
uzno = extr uzne ordo ;
zz   = extr uzne absc ; zn = zz/DJ ;
uznn = evol manu 'Z/Dj' zn 'Vaxe/Vj' uzno ;
 
*------------------------ valeurs experimentales de U ------------------
 
ue = 5.8 * (zn**(-1)) ;
uzne = evol manu 'Z/Dj' zn 'Vaxe/Vj' ue ;
tabu = table ;
tabu . 1 = 'REGU MARQ CARR ' ;
tabu . 2 = 'REGU MARQ ETOI ' ;
tabu.'TITRE' = table ;
tabu.'TITRE' . 1 = MOT 'CASTEM' ;
tabu.'TITRE' . 2 = MOT 'THEORIE' ;
 
titr 'VITESSE SUR L AXE' ;
dess (uznn et uzne)  lege tabu xbord 0. 50. ybord 0.0 1.1  ;
 
*------------------------- u' sur l'axe --------------------------------
 
kn = abs kn ;
up = kn*2./3. ; up = up**0.5 ;
upaxe = evol chpo up scal axe ;
uz1 = nomc scal uz ;
upu = up/uz1 ;
upue = evol chpo upu scal axe ;
upuo = extr upue ordo ;
zz   = extr upue absc ; zn = zz/DJ ;
upun = evol manu 'Z/Dj' zn 'u/Vaxe' upuo ;
 
*------------------------ valeurs experimentales de u' -----------------
*  Panchapakesan & Lumley 1993, et Wygnanski & Fiedler 1969
 
lzz = prog 30.   35.    40.    45.    50.   ;
lup = prog 0.222 0.2245 0.2265 0.2283 0.230 ;
lvp = prog 0.177 0.178  0.179  0.1796 0.180 ;
upe = evol manu 'Z/D' lzz 'Uprime/Umoyen' lup ;
vpe = evol manu 'Z/D' lzz 'Uprime/Umoyen' lvp ;
lz1 = prog 40. ; up1 = prog 0.268 ;
lz2 = prog 30. 40. 50. ;
vp2 = prog 0.183 0.1975 0.211  ;
wp2 = prog 0.188 0.2035 0.2075 ;
u1e = evol manu 'Z/D' lz1 'Uprime/Umoyen' up1 ;
v1e = evol manu 'Z/D' lz2 'Uprime/Umoyen' vp2 ;
w1e = evol manu 'Z/D' lz2 'Uprime/Umoyen' wp2 ;
 
tabp = table ;
tabp.1 = '          REGU ' ;
tabp.2 = 'MARQ PLUS NOLI ' ;
tabp.3 = 'MARQ CROI NOLI ' ;
tabp.4 = 'MARQ TRIA NOLI ' ;
tabp.5 = 'MARQ LOSA NOLI ' ;
tabp.6 = 'MARQ CARR NOLI ' ;
tabp.'TITRE' = table ;
tabp.'TITRE' . 1 = MOT 'CASTEM' ;
tabp.'TITRE' . 2 = MOT 'u/V P&L' ;
tabp.'TITRE' . 3 = MOT 'v/V P&L' ;
tabp.'TITRE' . 4 = MOT 'u/V W&F' ;
tabp.'TITRE' . 5 = MOT 'v/V W&F' ;
tabp.'TITRE' . 6 = MOT 'w/V W&F' ;
 
titre 'FLUCTUATIONS DE VITESSE SUR L AXE' ;
dess (upun et upe et vpe et u1e et v1e et w1e) lege tabp ;
 
*------------------------------- les rayons ---------------------------
 
P01 = MT poin proc (0. (10.*DJ));
PJ1 = MT poin proc (RJ (10.*DJ));
PR1 = MT poin proc ((50.*DJ) (10.*DJ)) ;
* ordonnee de p01 = 0.2035 si RJ = 0.01 m
l101 = droi 2 P01 PJ1 ;
l102 = droi   PJ1 PR1 dini (RJ/2.) dfin (10.*RJ) ;
L10 = l101 et l102 ;  elim mt l10 1.e-5 ;
 
PJ2 = RJ (20.*DJ) ;
l201 = droi 2 P02 PJ2 ;
l202 = droi   PJ2 PR2 dini (RJ/2.) dfin (10.*RJ) ;
L20 = l201 et l202 ;  elim mt l20 1.e-5 ;
 
PJ3 = RJ (30.*DJ) ;
l301 = droi 2 P03 PJ3 ;
l302 = droi   PJ3 PR3 dini (RJ/2.) dfin (10.*RJ) ;
L30 = l301 et l302 ;  elim mt l30 1.e-5 ;
 
PJ4 = RJ (40.*DJ) ;
l401 = droi 2 P04 PJ4 ;
l402 = droi   PJ4 PR4 dini (RJ/2.) dfin (10.*RJ) ;
L40 = l401 et l402 ;  elim mt l40 1.e-5 ;
 
PJ5 = RJ (50.*DJ) ;
l501 = droi 2 P05 PJ5 ;
l502 = droi   PJ5 PR5 dini (RJ/2.) dfin (10.*RJ) ;
L50 = l501 et l502 ;  elim mt l50 1.e-5 ;
 
*---------------------------- U sur les rayons -------------------------
 
uz20 = evol chpo uz l20 ;
u20 = extr uz20 ordo ; uzn20 = u20/(maxi u20) ;
r20 = extr uz20 absc ; rn20 = r20/(20.*DJ) ;
un20 = evol manu 'R/Z' rn20 'V/Vaxe' uzn20 ;
 
uz30 = evol chpo uz l30 ;
u30 = extr uz30 ordo ; uzn30 = u30/(maxi u30) ;
r30 = extr uz30 absc ; rn30 = r30/(30.*DJ) ;
un30 = evol manu 'R/Z' rn30 'Uz/Uzmax' uzn30 ;
 
uz40 = evol chpo uz l40 ;
u40 = extr uz40 ordo ; uzn40 = u40/(maxi u40) ;
r40 = extr uz40 absc ; rn40 = r40/(40.*DJ) ;
un40 = evol manu 'R/Z' rn40 'Uz/Uzmax' uzn40 ;
 
*------------------- valeurs experimentales ----------------------------
*     ut=formule de Tollmien, ue=valeurs exp de Trupel
 
lz=prog 0.000 0.009 0.018 0.027 0.036 0.045 0.054 0.063 0.072 0.081
        0.090 0.099 0.108 0.117 0.126 0.135 0.144 0.153 0.162 0.171
        0.180 0.189 0.198 0.207 0.216 0.225 0.234 0.243 0.252 0.261
        0.270 0.279 0.288 0.297 0.306                               ;
ut=prog 1.000 0.984 0.958 0.922 0.884 0.843 0.795 0.748 0.700 0.653
        0.605 0.555 0.510 0.470 0.425 0.378 0.340 0.300 0.265 0.230
        0.198 0.169 0.140 0.117 0.094 0.075 0.059 0.046 0.034 0.024
        0.017 0.011 0.007 0.003 0.                                  ;
ue=prog 1.000 0.991 0.965 0.936 0.898 0.864 0.809 0.762 0.707 0.700
        0.612 0.562 0.517 0.475 0.425 0.374 0.342 0.298 0.268 0.228
        0.196 0.171 0.138 0.119 0.096 0.075 0.056 0.046 0.035 0.025
        0.017 0.012 0.009 0.004 0.                                  ;
uze = evol MANU 'R/Z' lz 'U/Umax' ue ;
uzt = evol MANU 'R/Z' lz 'U/Umax' ut ;
 
tabv = table ;
tabv.1 = 'MARQ TRIA REGU ' ;
tabv.2 = 'MARQ LOSA REGU ' ;
tabv.3 = 'MARQ CARR REGU ' ;
tabv.4 = 'MARQ CROI NOLI ' ;
tabv.5 = 'MARQ PLUS NOLI ' ;
tabv.'TITRE' = table ;
tabv.'TITRE' . 1 = MOT 'CASTEM Z=20Dj' ;
tabv.'TITRE' . 2 = MOT 'CASTEM Z=30Dj' ;
tabv.'TITRE' . 3 = MOT 'CASTEM Z=40Dj' ;
tabv.'TITRE' . 4 = MOT 'THEORIE'   ;
tabv.'TITRE' . 5 = MOT 'EXPERIENCE' ;
 
titr 'PROFILS RADIAUX DE VITESSE MOYENNE AXIALE' ;
dess (un20 et un30 et un40 et uzt et uze)
                       lege tabv ybord 0. 1.  xbor 0. 0.3 ;
 
*---------------------------- historiques ------------------------------
 
titre 'HISTORIQUE DE LA VITESSE RADIALE' ;
  dess rv.hist.'1UN' ;
titre 'HISTORIQUE DE LA VITESSE AXIALE' ;
  dess rv.hist.'2UN' ;
 
FINSI ;
*----------------test sur le débit total -------------------------------
 
debent= DBIT $bas  RV.INCO.'UN' ;
debsor= DBIT $haut RV.INCO.'UN' ;
debbor= DBIT $bord RV.INCO.'UN' ;
debtot= debent + debbor - debsor ;
 
SI ( (ABS debtot) > 1.e-10 ) ;
    ERREUR 5 ;
FINSI ;
 
FIN ;