* fichier : panach1.dgibi
************************************************
* JET/PANACHE 2D - SEMI-INFINI *
* Comparaison K - Epsilon / RNG K - Epsilon *
* il suffit pour cela d'enlever ou remettre*
* l'option RNG devant NSKE *
* Mode axisymetrique *
* La procedure FILTREKE est incorporee *
* Le Richardson est parametrable Ri = -5.e-2 *
* Gregory Turbelin 29/12/1998 *
************************************************
OPTI DIME 2 ELEM QUA8;
OPTI MODE AXIS; OPTI ISOV SULI;
GRAPH='N' ;
COMPLET=FAUX ;
** PROCEDURE DE FILTRE DE LA CONCENTRATION
DEBP 4FILTRE ; ARGU RX*TABLE ;
rv=rx.'EQEX' ; rvp=rv.'PRESSION' ;
cn=rv.inco.'cn';
cn=KOPS cn '|<' 0.;
cn=KOPS cn '>|' 1.;
rv.inco.'cn'=cn ;
FINPROC ;
DEBP FILTREKE ;
ARGU RX*TABLE ;
* Filtre sur K et Epsilon
* - Echelle de vitesse (K**0.5) inférieure à une fraction (alfk)
* de Uref (vitesse caractéristique) (alfk=1 pour l'instant)
* Uref=max(UN,U0)
* - K > K0
* - Epsilon tel que l'echelle de longueur reste inférieure
* à (L0/a) où L0 = diamètre enceinte et a=f(Re)
* => Nut < Uref*L0/a
rv=rx.'EQEX' ;
rvp=rv.'PRESSION' ;
iarg=rx.'IARG' ;
si( non ( ega iarg 4)) ;
mess 'Procedure FILTREKE : nombre d arguments incorrect ' iarg ;
quitter FILTREKE ;
finsi ;
si ( ega ('TYPE' rx.'ARG1') 'MOT ') ;
U1=rv.'INCO'.(rx.'ARG1') ;
sinon ;
si ( ega ('TYPE' (rx.'ARG1')) 'FLOTTANT') ;
U1=rx.'ARG1' ;
sinon ;
mess 'Procedure FILTREKE : type argument 1 invalide ' ;
quitter FILTREKE ;
finsi ;
finsi ;
si ( ega ('TYPE' rx.'ARG2') 'MOT ') ;
L0=rv.'INCO'.(rx.'ARG2') ;
sinon ;
si ( ega ('TYPE' (rx.'ARG2')) 'FLOTTANT') ;
L0=rx.'ARG2' ;
sinon ;
mess 'Procedure FILTREKE : type argument 2 invalide ' ;
quitter FILTREKE ;
finsi ;
finsi ;
si ( ega ('TYPE' rx.'ARG3') 'MOT ') ;
NU=rv.'INCO'.(rx.'ARG3') ;
sinon ;
si ( ega ('TYPE' (rx.'ARG3')) 'FLOTTANT') ;
NU=rx.'ARG3' ;
sinon ;
mess 'Procedure FILTREKE : type argument 3 invalide ' ;
quitter FILTREKE ;
finsi ;
finsi ;
si ( ega ('TYPE' rx.'ARG4') 'MOT ') ;
UN=rv.'INCO'.(rx.'ARG4') ;
sinon ;
si ( ega ('TYPE' (rx.'ARG4')) 'CHPOINT') ;
UN=rx.'ARG4' ;
sinon ;
mess 'Procedure FILTREKE : type argument 4 invalide ' ;
quitter FILTREKE ;
finsi ;
finsi ;
nic=dime (rx.'LISTINCO') ;
si( non ( ega nic 2)) ;
mess 'Procedure FILTREKE : nombre d inconnues incorrect ' nic ;
quitter FILTREKE ;
finsi ;
nomi1=extr 1 (rx.'LISTINCO');
nomi2=extr 2 (rx.'LISTINCO');
nom1= mot (text (chai nomi1)) ;
nom2= mot (text (chai nomi2)) ;
en=rv.'INCO'.nom2 ;
kn=rv.'INCO'.nom1 ;
Rec=100.;
k0 = 1.e-10 ;
cnu=0.09;
mdu=kops un 'PSCA' un ;
mdu=mdu ** 0.5 ;
Re=kops (kops (kops mdu '*' L0) '/' nu) '+' (Rec / 10.) ;
a= exp (kops Rec '/' Re ) ;
mdu = kops mdu '|<' u1 ;
mdu2= kops mdu '*' mdu ;
kn=kops kn '|<' k0 ;
kn=kops kn '>|' mdu2 ;
E0= kops (kops kn '**' 1.5) '*' (a / L0) ;
en=kops en '|<' E0 ;
rv.'INCO'.nom2=en ;
rv.'INCO'.nom1=kn ;
FINPROC ;
****************************
* CONSTRUCTION DU MAILLAGE *
****************************
DJ=1. ;RJ=DJ/2. ;
** POINTS
P0=0. 0.; P1=RJ 0.;
P2=8. 0.; P3=8. 15.;
P4=RJ 15.; P5=0. 15.;
P1b=4. 0. ; P2b=8. 8.;
P3b=4. 15.; P5b=0. 8.;
P6=8. 5.; P7=8. 10.;
P9=0. 10.; P10=0. 5.;
Pe1=20. 0. ; Pe2=20. 8.;
Pe3=20. 15.; Pe4=20. 50.;
Pe5=8. 50. ; Pe6=4. 50.;
Pe7=RJ 50. ; Pe8=0. 50. ;
***** SEGMENTS *****
Si complet ;
N1=5 ; N2=-18 ; N3=-10 ; N4=-12 ; N5=-30 ; N6=-38 ;
sinon ;
N1=1 ; N2=-4 ; N3=-2 ; N4=-3 ; N5=-7 ; N6=-9 ;
finsi ;
ENT=P0 DROI N1 P1;
BA1=P1 DROI N2 P1b DINI (0.008*DJ) DFIN (0.1*DJ);
BA2=P1b DROI N3 P2 DINI (1.1*DJ) DFIN (1.4*DJ);
BA3=P2 DROI N4 Pe1 DINI (1.6*DJ) DFIN (3.*DJ);
BAS=BA1 ET BA2 ET BA3;
BAT=ENT ET BAS;
CO1=Pe1 DROI N5 Pe2 DINI (0.008*DJ) DFIN (0.1*DJ);
CO2=Pe2 DROI N2 Pe3 DINI (1.35*DJ) DFIN (1.5*DJ);
CO3=Pe3 DROI N6 Pe4 DINI (1.8*DJ) DFIN (3.*DJ);
COT=CO1 ET CO2 ET CO3;
HT1=Pe4 DROI N4 Pe5 DINI (3.*DJ) DFIN (1.6*DJ);
HT2=Pe5 DROI N3 Pe6 DINI (1.4*DJ) DFIN (1.1*DJ);
HT3=Pe6 DROI N2 Pe7 DINI (0.1*DJ) DFIN (0.008*DJ);
HT4=Pe7 DROI N1 Pe8;
HTT=HT1 ET HT2 ET HT3 ET HT4;
AX1=Pe8 DROI N6 P5 DINI (3.*DJ) DFIN (1.8*DJ);
AX2=P5 DROI N2 P5b DINI (1.5*DJ) DFIN (1.35*DJ);
AX3=P5b DROI N5 P0 DINI (0.1*DJ) DFIN (0.008*DJ);
AXE=AX1 ET AX2 ET AX3;
ELIM (BAT ET COT ET HTT ET AXE) 1.e-5;
MT=DALL BAT COT HTT AXE PLAN;
ORIENTER MT;
***** TABLES DOMAINE *****
$MT=DOMA MT 1.e-5 MACRO;
$COT=DOMA COT INCL $MT 1.e-5 MACRO;
$AXE=DOMA AXE INCL $MT 1.e-5 MACRO;
$BAT=DOMA BAT INCL $MT 1.e-5 MACRO;
$HTT=DOMA HTT INCL $MT 1.e-5 MACRO;
MT=$MT.maillage;
* TRAC MT;
**********************
* DONNEES PHYSIQUES *
**********************
rhoa=2*1.19; rhof=2*0.58;
mua=1.8e-5 ; nua=mua/rhoa;
cref=0. ; Ri=0. -5.e-2;
Re=3000. ; iRe=1./Re;
db=1.e-5 ; iSc=db/nua;
RS=iRe*iSc ; Sct=0.7;
***DONNEES POUR FILTREKE
U0=1.; L0=10.;
***************
* EQUATIONS *
***************
rv=EQEX $MT ITMA 80 ALFA 0.5
OPTI 'SUPG' 'RNG'
ZONE $MT OPER FILTREKE U0 L0 nua 'un' INCO 'kn' 'en'
ZONE $MT OPER NSKE iRe 'nut' Ri 'cn' cref INCO 'un' 'kn' 'en'
ZONE $MT OPER 4FILTRE
ZONE $MT OPER TSCA RS 'un' 0. 'nut' Sct INCO 'cn';
rv.nomvi='un';
rv=eqex rv
CLIM 'un' VIMP ENT 1.
'un' VIMP BAS 0.
'un' UIMP AXE 0.
'un' VIMP COT 0.
'kn' TIMP ENT 2.e-5
'en' TIMP ENT 1.e-6;
rv=eqex rv
CLIM 'cn' TIMP ENT 1.
'cn' TIMP COT 0.
'cn' TIMP BAS 0.
'kn' TIMP COT 1.e-5
'en' TIMP COT 1.e-6;
rvp=EQPR $MT KTYPI 1
ZONE $MT OPER pression 0.;
** INITIALISATIONS
rv.inco=TABLE INCO;
rv.inco.'un'=KCHT $MT VECT SOMMET (1.e-7 1.e-7);
rv.inco.'kn'=KCHT $MT SCAL SOMMET 1.e-3;
rv.inco.'en'=KCHT $MT SCAL SOMMET 1.e-4;
rv.inco.'nut'=KCHT $MT SCAL CENTRE nua;
rv.inco.'cn'=KCHT $MT SCAL SOMMET 0.;
rv.pression=rvp;
lh= (POIN MT PROC (0.75 1.)) et (POIN MT PROC (0.75 5.))
et (POIN MT PROC (0.75 10.)) et (POIN MT PROC (0.75 15.))
et (POIN MT PROC (0.75 20.)) et (POIN MT PROC (0.75 25.))
et (POIN MT PROC (0.75 33.)) et (POIN MT PROC (0.75 40.))
et (POIN MT PROC (0.75 45.)) et (POIN MT PROC (0.75 50.));
his=khis 'un' 1 lh
'un' 2 lh
'kn' lh
'cn' lh
'en' lh;
rv.hist=his;
EXEC rv;
** POST-TRAITEMENT
**----------------
** HISTORIQUES
**------------
TAB1=TABLE;
TAB1.'TITRE'= table ;
TAB1.1 ='TIRR MARQ CROI';
TAB1.'TITRE' . 1 ='P1';
TAB1.2 ='TIRR MARQ TRIA';
TAB1.'TITRE' . 2 ='P2';
TAB1.3 ='TIRR MARQ CARR';
TAB1.'TITRE' . 3 ='P3';
TAB1.4 ='TIRR MARQ ETOI';
TAB1.'TITRE' . 4 ='P4';
TAB1.5 ='TIRR MARQ PLUS';
TAB1.'TITRE' . 5 ='P5';
TAB1.6 ='TIRR MARQ LOSA';
TAB1.'TITRE' . 6 ='P6';
TAB1.7 ='TIRR MARQ TRIB';
TAB1.'TITRE' . 7 ='P7';
si ('EGA' graph 'O' );
DESS (his.'1un') LEGE TAB1 ;
DESS (his.'2un') LEGE TAB1 ;
DESS (his.'kn') LEGE TAB1;
DESS (his.'cn') LEGE TAB1;
DESS (his.'en') LEGE TAB1;
finsi ;
** PROFILS SUR LE DOMAINE TOTAL
**------------------------------
vv = vect rv.inco.'un' 5. ux uy jaune ;
kk = rv.inco.'kn' ;
ee = rv.inco.'en';
cc = rv.inco.'cn' ;
nn = elno (rv.inco.'nut') $MT;
si ('EGA' graph 'O' );
titre 'Vecteur Vitesse';
trac vv MT (CONT MT); ;
titre 'Energie Cinetique Turbulente';
trac kk MT (CONT MT);
titre 'NUT';
trac nn MT (CONT MT);
titre 'dissipation';
trac ee MT (CONT MT);
finsi ;
** PROFILS A DIFFERENTES ALTITUDES
**---------------------------------
LMT=CHAN LIGNE MT;
** ALTITUDES
alt1 = 1. ; alt2 = 5.;
alt3 = 10.; alt4 = 25.;
** LIGNES HORIZONTALES
lign1=POIN MT DROI (POIN MT PROC (0. alt1))
(POIN MT PROC (20. alt1)) 1.e-3;
lign2=POIN MT DROI (POIN MT PROC (0. alt2))
(POIN MT PROC (20. alt2)) 1.e-3;
lign3=POIN MT DROI (POIN MT PROC (0. alt3))
(POIN MT PROC (20. alt3)) 1.e-3;
lign4=POIN MT DROI (POIN MT PROC (0. alt4))
(POIN MT PROC (20. alt4)) 1.e-3;
L1=ELEM LMT APPUYE lign1;
L2=ELEM LMT APPUYE lign2;
L3=ELEM LMT APPUYE lign3;
L4=ELEM LMT APPUYE lign4;
nut=elno (rv.inco.'nut') (rv.domaine);
prvy1=EVOL CHPO (rv.inco.'un') uy L1;
prvx1=EVOL CHPO (rv.inco.'un') ux L1;
prc1=EVOL CHPO (rv.inco.'cn') scal L1;
prk1=EVOL CHPO (rv.inco.'kn') scal L1;
prnut1=EVOL CHPO nut scal L1;
prvy2=EVOL CHPO (rv.inco.'un') uy L2;
prvx2=EVOL CHPO (rv.inco.'un') ux L2;
prc2=EVOL CHPO (rv.inco.'cn') scal L2;
prk2=EVOL CHPO (rv.inco.'kn') scal L2;
prnut2=EVOL CHPO nut scal L2;
prvy3=EVOL CHPO (rv.inco.'un') uy L3;
prvx3=EVOL CHPO (rv.inco.'un') ux L3;
prc3=EVOL CHPO (rv.inco.'cn') scal L3;
prk3=EVOL CHPO (rv.inco.'kn') scal L3;
prnut3=EVOL CHPO nut scal L3;
prvy4=EVOL CHPO (rv.inco.'un') uy L4;
prvx4=EVOL CHPO (rv.inco.'un') ux L4;
prc4=EVOL CHPO (rv.inco.'cn') scal L4;
prk4=EVOL CHPO (rv.inco.'kn') scal L4;
prnut4=EVOL CHPO nut scal L4;
TAB2=TABLE;
TAB2.'TITRE'= TABLE ;
TAB2.1='TIRR ';
TAB2.'TITRE' . 1 = MOT 'Z=1';
TAB2.2='TIRC';
TAB2.'TITRE' . 2 = MOT 'Z=5';
TAB2.3='TIRL';
TAB2.'TITRE' . 3 = MOT 'Z=10';
TAB2.4='TIRM';
TAB2.'TITRE' . 4 = MOT 'Z=25';
** ZOOM SUR L'AXE X
xbinf = 0.;
xbsup = 6.;
xbsup2 =1.2;
si ('EGA' graph 'O' );
dess (prvy1 et prvy2 et prvy3 et prvy4) titr 'PROFIL VITESSE
(Uz)' Xbor xbinf xbsup LEGE TAB2 titx 'R/DJ' tity 'Uz';;
dess (prc1 et prc2 et prc3 et prc4) titr 'PROFIL
CONCENTRATION' Xbor xbinf xbsup LEGE TAB2 titx 'R/DJ' tity 'Cf';
dess (prk1 et prk2 et prk3 et prk4) titr 'PROFIL EN CINET
TURB' Xbor xbinf xbsup LEGE TAB2 titx 'R/DJ' tity 'K';
dess (prnut1 et prnut2 et prnut3 et prnut4) titr 'PROFIL VISCO
TURB' Xbor xbinf xbsup LEGE TAB2 titx 'R' tity 'Nut';
finsi ;
** DIAMETRES CARACTERISTIQUES
**----------------------------
DEBP ADIM VISU*EVOLUTIO ;
LV = EXTR VISU 'ORDO' ;
LX = EXTR VISU 'ABSC' ;
** Recherche des grandeurs d'adimensionnalisation
Vaxe = EXTR LV 1 ;
V12 = Vaxe / 2. ;
i = 1 ;
REPETER BOUCLE (dime LV) ;
Vcour = EXTR LV i ;
SI (Vcour < V12) ; quitter boucle ;
SINON ; i = i + 1 ;
FINSI ;
FIN BOUCLE ;
** Determination de X12
Xi = EXTR LX i ;
Xii = EXTR LX (i - 1) ;
Vi = EXTR LV i ;
Vii = EXTR LV (i - 1) ;
X12 = (((Vii - V12) * Xi) + ((V12 -Vi) * Xii)) / (Vii - Vi) ;
** Adimensionnalisation
LVA = LV / (Vaxe + 1e-20) ;
LXA = LX / (X12 + 1e -20) ;
VISUA = EVOL 'MANU' 'X / X1/2' LXA 'V / Vmax' LVA ;
FINP X12 ;
yaxe=COOR 2 $AXE.maillage;
lyaxe=PROG;
i=0;
repe bloc (nbno $AXE.maillage);
i=i+1;
vyaxe=EXTR yaxe SCAL ($AXE.maillage POIN i);
lyaxe=lyaxe et (prog vyaxe);
fin bloc;
lyaxe = ORDO lyaxe;
LMT=CHAN LIGNE MT;
diam2=prog;
i=0;
repeter boucl3 (dime lyaxe);
i=i+1;
*mess 'i vaut' i;
yi=extr lyaxe i;
lign6=POIN MT DROI (POIN MT PROC (0. yi))
(POIN MT PROC (10. yi)) 1.e-3;
L6=elem lmt appuye lign6;
visu1=evol chpo rv.inco.'un' uy L6;
XD=adim visu1;
diam2=diam2 et (prog XD);
menage;
fin boucl3;
dcar1=evol manu diam2 lyaxe;
si ('EGA' graph 'O' );
dess (dcar1) titx 'R' tity 'Z';
finsi ;
** GRANDEURS SUR L'AXE ADIMENSIONNEES
**------------------------------------
** Vitesse
si complet ;
evv4=EVOL CHPO (rv.inco.'un') uy (INVE AXE);
ev4=EXTR evv4 ORDO ; ev4=ENLE ev4 173;
evn4=LOG ev4;
z4=EXTR evv4 ABSC ; z4=ENLE z4 1; zn4=LOG z4;
nevv4=EVOL MANU zn4 evn4;
dess evv4 titr 'UY SUR L AXE (ADIM)';
dess nevv4 titr 'UY SUR L AXE (ADIM) ECHELLE LOG';
finsi ;
** Concentration
** Attention, on ne peut calculer la concentration
** qu'à partir d'un grand nombre de pas de temps
*
* Xa=287.; Xhe=2079.;
* cn1=KOPS (KOPS Xa '*' (rv.inco.'cn')) '/'
* (KOPS (KOPS (KOPS 1 '-' (rv.inco.'cn')) '*' Xhe)
* '+' (KOPS (rv.inco.'cn') '*' Xa));
* evc1=EVOL CHPO cn1 SCAL (INVE AXE);
* ec1=extr evc1 ORDO; ec1=ENLE ec1 1; ecn1=LOG ec1;
* zc1=extr evc1 ABSC ;zc1=ENLE zc1 1; zcn1=LOG zc1;
* nevc1=EVOL MANU zcn1 ecn1;
*dess evc1 titr 'CONCENTRATION SUR L AXE (ADIM)';
*dess nevc1 titr 'CONCENTRATION SUR L AXE (ADIM) ECHELLE LOG';
* TAB3=TABLE;
* TAB3.'TITRE'= TABLE ;
* TAB3.'TITRE' . 1 = MOT 'VITESSE AXIALE';
* TAB3.2='TIRC';
* TAB3.'TITRE' . 2 = MOT 'CONCENTRATION';
* dess (evv4 et evc1) LEGE TAB3;
* dess (nevv4 et nevc1) LEGE TAB3;
** Energie
evk1=EVOL CHPO (rv.inco.'kn') SCAL (INVE AXE);
ek1=EXTR evk1 ORDO ;ek1=ENLE ek1 1;ek1=ABS ek1; ekn1=LOG ek1;
zk1=EXTR evk1 ABSC ;zk1=ENLE zk1 1;zk1=ABS zk1; zkn1=LOG zk1;
nevk1=EVOL MANU 'Z' zkn1 'K' ekn1;
si ('EGA' graph 'O' );
dess evk1 titr 'EN. CINET. TURB. SUR L AXE (ADIM)';
dess nevk1 titr 'EN. CINET. TURB. SUR L AXE (ADIM) ECHELLE LOG';
finsi ;
* fin;
**************************
*RESULTATS COMPLEMENTAIRES*
***************************
** LIGNES DE COURANT
**-----------------
corx=COOR 1 MT;
un=rv.inco.'un';
un1=EXCO un ux; un1=-1*corx*un1;
un2=EXCO un uy; un3=corx*un2;
un1=KCHT $MT SCAL SOMMET un1;
un2=KCHT $MT SCAL SOMMET un2;
un3=KCHT $MT SCAL SOMMET un3;
uncz=NOEL $MT un2;
uncz=KCHT $MT SCAL CENTRE uncz;
rt2d=KOPS un 'ROT' $MT;
corx=KCHT $MT SCAL SOMMET corx;
corm=NOEL $MT corx;
unn1=KCHT $BAT SCAL SOMMET un1;
unn2=KCHT $COT SCAL SOMMET un3;
unn3=KCHT $HTT SCAL SOMMET (-1*un1);
unn4=KCHT $AXE SCAL SOMMET (-1*un3);
unn1=NOEL $BAT unn1;
unn1=KCHT $BAT SCAL CENTRE unn1;
unn2=NOEL $COT unn2;
unn2=KCHT $COT SCAL CENTRE unn2;
unn3=NOEL $HTT unn3;
unn3=KCHT $HTT SCAL CENTRE unn3;
unn4=NOEL $AXE unn4;
unn4=KCHT $AXE SCAL CENTRE unn4;
sw=2*uncz - (corm*rt2d);
rk=EQEX $MT 'OPTI' 'EF' 'IMPL'
ZONE $MT OPER LAPN 1. inco 'psi'
ZONE $MT OPER FIMP sw inco 'psi'
ZONE $BAT OPER FIMP unn1 inco 'psi'
ZONE $COT OPER FIMP unn2 inco 'psi'
ZONE $HTT OPER FIMP unn3 inco 'psi'
ZONE $AXE OPER FIMP unn4 inco 'psi'
CLIM 'psi' TIMP BAS 0.;
rk.inco.'psi'=KCHT $MT SCAL SOMMET 0.;
exec rk;
si ('EGA' graph 'O' );
psi=rk.inco.'psi';
OPTI ISOV LIGNE;
TRAC psi MT (CONT (MT)) titr 'LIGNES DE COURANT';
psi1=EVOL CHPO psi SCAL COT;
db1=EXTR psi1 ORDO;db2=EXTR psi1 ABSC;
npsi1=EVOL MANU 'Z' db2 'DEBIT' db1 ;
DESS npsi1;
titr 'DEBIT D ENTRAINEMENT';
FINSI ;
FIN;
© Cast3M 2003 - All rights reserved.
Disclaimer