* fichier : gridturb.dgibi
**********************************************************
* GRID TURBULENCE : convection of homogeneous turbulence *
* Analysis of the K-Epsilon TURBULENCE MODEL *
* Mohammadi/Pironneau p. 74 (Wiley) *
* H. PAILLERE/TTMF/AVRIL 1997 (à l'origine 1/2 explicite)*
* Version implicite MARS 2000 *
* Teste K-Epsilon en implicite méthode de projection *
**********************************************************
* WARNING : Il faut que Nt/Nu > 1
*opti trace 'PSC';
GRAPH = VRAI ;
GRAPH = FAUX ;
KPRES='MSOMMET';
DISCR='LINE';
KSUPG='SUPGDC';
*
D = 1 ;
NITMA =200 ;
LCFL = PROG 0.1 0.5 1. 10.;
LCFL = PROG 1. 10. 50. 100.;
lerkm= prog -1.5 -2. -4. -4. ;
lerem= prog -1.5 -2. -4. -4. ;
NX = D*40 ;
NY = D ;
opti dime 2 ;
opti elem qua8 ;
opti isov suli ;
******************************************************************
* PROCEDURE POUR ESTIMER LA CONVERGENCE VERS L'ETAT STATIONNAIRE *
******************************************************************
DEBPROC CALCUL ;
ARGU RVX*'TABLE' ;
RV = RVX.'EQEX' ;
DD = RV.PASDETPS.'NUPASDT' ;
NN = DD/10 ;
LO = (DD-(10*NN)) EGA 0 ;
SI ( LO ) ;
EN = RV.INCO.'EN' ;
ENM1 = RV.INCO.'ENM1' ;
KN = RV.INCO.'KN' ;
KNM1 = RV.INCO.'KNM1' ;
ERRE = KOPS EN '-' ENM1 ;
ELIE = MAXI ERRE 'ABS' ;
ELIE = (LOG (ELIE + 1.0E-20))/(LOG 20.) ;
ERRK = KOPS KN '-' KNM1 ;
ELIK = MAXI ERRK 'ABS' ;
ELIK = (LOG (ELIK + 1.0E-20))/(LOG 20.) ;
MESSAGE 'ITER ' RV.PASDETPS.'NUPASDT' ' ERREUR LINF ' ELIK ELIE ;
RV.INCO.'ENM1' = KCHT $DOMTOT 'SCAL' 'SOMMET' (RV.INCO.'EN') ;
RV.INCO.'KNM1' = KCHT $DOMTOT 'SCAL' 'SOMMET' (RV.INCO.'KN') ;
IT = PROG RV.PASDETPS.'NUPASDT' ;
ERK = PROG ELIK ;
ERE = PROG ELIE ;
RV.INCO.'IT' = (RV.INCO.'IT') ET IT ;
RV.INCO.'ERK' = (RV.INCO.'ERK') ET ERK ;
RV.INCO.'ERE' = (RV.INCO.'ERE') ET ERE ;
FINSI ;
as2 ama1 = 'KOPS' 'MATRIK' ;
FINPROC as2 ama1 ;
***********************
* GEOMETRY OF PROBLEM *
***********************
A1 = 0.0 0.0 ;
A2 = 1.0 0.0 ;
A3 = 1.0 0.02;
A4 = 0.0 0.02;
A1A2 = A1 'DROI' NX A2 ;
A2A3 = A2 'DROI' NY A3 ;
A3A4 = A3 'DROI' NX A4 ;
A4A1 = A4 'DROI' NY A1 ;
DOMTOT = 'DALLER' A1A2 A2A3 A3A4 A4A1 'PLAN' ;
DOMTOT = CHAN DOMTOT QUAF ;
A4A1 = CHAN A4A1 QUAF ;
A2A3 = CHAN A2A3 QUAF ;
A1A2 = CHAN A1A2 QUAF ;
A3A4 = CHAN A3A4 QUAF ;
ELIM (DOMTOT et A4A1 et A2A3 et A3A4 et A1A2) 1.E-9 ;
MCOTE = A1A2 et A3A4;
$DOMTOT = 'MODE' DOMTOT 'NAVIER_STOKES' DISCR ;
DOMA $DOMTOT IMPR ;
$INLET = 'MODE' A4A1 'NAVIER_STOKES' DISCR ;
INLET = 'DOMA' $INLET 'MAILLAGE' ;
$OUTLET = 'MODE' A2A3 'NAVIER_STOKES' DISCR ;
OUTLET = 'DOMA' $OUTLET 'MAILLAGE' ;
poutlet= 'DOMA' $OUTLET 'MSOMMET' ;
$COTE = 'MODE' MCOTE 'NAVIER_STOKES' DISCR ;
COTE = DOMA $COTE MAILLAGE ;
DOMTOT = DOMA $DOMTOT MAILLAGE ;
CNT = CONT DOMTOT ;
*****************
* PHYSICAL DATA *
*****************
Uinlet = 1.71;
Kinlet = 0.013;
Einlet = 0.029;
deltax = (coor 1 A2)/nx ;
Nu = 1.E-3 ;
*****************************************
* ANALYTICAL SOLUTION FOR K(x) AND E(x) *
*****************************************
C2 = 1.92 ;
X1 = COOR 1 (DOMA $DOMTOT SOMMET) ;
TINLET = KINLET/EINLET;
TAN=(((C2 - 1.)*X1) + (UINLET*TINLET))/UINLET;
M=1. - C2 ;
KAN=KINLET * (( 1. - ((M*X1)/UINLET/TINLET))**(1./M));
EAN=KAN*(INVE TAN) ;
********************
* SET OF EQUATIONS *
********************
RV = EQEX $DOMTOT 'ITMA' NITMA 'OMEGA' 1. 'NITER' 1 'FIDT' 20
'ZONE' $DOMTOT 'OPER' CALCUL
'OPTI' KSUPG 'IMPL' EF
'ZONE' $DOMTOT 'OPER' 'NS' 'Ro' 'UN' 'MUF' 'INCO' 'UN'
'ZONE' $DOMTOT 'OPER' 'KEPSILON' 'Ro' 'UN' 'Mu' 'DT' 'INCO' 'KN' 'EN'
OPTI EFM1 'CENTREE'
'ZONE' $DOMTOT 'OPER' 'DFDT' 'Ro' 'UN' 'DT' 'INCO' 'UN'
;
RV = EQEX RV
'CLIM' 'UN' 'UIMP' INLET Uinlet
'UN' 'VIMP' INLET 0.
'UN' 'VIMP' COTE 0.
'CLIM' 'KN' 'TIMP' INLET Kinlet
'EN' 'TIMP' INLET Einlet
;
RVP = EQEX 'OPTI' 'EF' KPRES
ZONE $DOMTOT 'OPER' KBBT (-1.) INCO 'UN' 'PRES'
CLIM 'PRES' TIMP poutlet 0.
;
rvp.'METHINV' . 'FCPRECT'=300 ;
RV.INCO = TABLE INCO ;
RV.INCO.'UN' = KCHT $DOMTOT VECT SOMMET (uinlet 0.) ;
RV.INCO.'KN' = KCHT $DOMTOT SCAL SOMMET Kinlet;
RV.INCO.'EN' = KCHT $DOMTOT SCAL SOMMET 1.e-3 ;
RV.INCO.'ENM1' = KCHT $DOMTOT SCAL SOMMET Einlet ;
RV.INCO.'KNM1' = KCHT $DOMTOT SCAL SOMMET Kinlet ;
RV.INCO.'MUF' = KCHT $DOMTOT SCAL SOMMET 1.E-10 ;
RV.INCO.'Ro' = KCHT $DOMTOT SCAL SOMMET 1. ;
RV.INCO.'Mu' = KCHT $DOMTOT SCAL SOMMET Nu ;
RV.INCO.'PRES' = KCHT $DOMTOT SCAL KPRES 0. ;
RV.INCO.'IT' = PROG 1 ;
RV.INCO.'ERK' = PROG 0.;
RV.INCO.'ERE' = PROG 0.;
RV.INCO.'EAN'=EAN;
RV.INCO.'KAN'=KAN;
RV.'PROJ' = RVP ;
EVOLKA = EVOL (BLEU) 'CHPO' KAN 'SCAL' A1A2 ;
EVOLEA = EVOL (BLEU) 'CHPO' EAN 'SCAL' A1A2 ;
NCFL = DIME LCFL ;
REPETER BCFL NCFL;
CFL = extr LCFL &BCFL ;
dt = CFL*deltax/Uinlet ;
mess ' CFL=' CFL ' deltax ' deltax ' dt = ' dt ;
RV.INCO.'DT' = dt;
RV.INCO.'IT' = PROG 1 ;
RV.INCO.'ERK' = PROG 0.;
RV.INCO.'ERE' = PROG 0.;
RV.'PASDETPS'.'NUPASDT'= 1;
EXEC RV ;
EVOLKA = EVOLKA et (EVOL (ROUG) 'CHPO' RV.INCO.'KN' 'SCAL' A1A2) ;
EVOLEA = EVOLEA et (EVOL (ROUG) 'CHPO' RV.INCO.'EN' 'SCAL' A1A2) ;
EVOLK = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
(RV.INCO.'ERK') ;
EVOLE = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
(RV.INCO.'ERE') ;
Si (EGA &BCFL 1);
EVOLER = EVOLK et EVOLE;
Sinon ;
EVOLER = EVOLK et EVOLE et EVOLER ;
Finsi ;
* TESTS D'ERREUR *
erkm= extr &bcfl lerkm;
erem= extr &bcfl lerem;
list rv.inco.'ERK' ;
mess ' Sur K' (MINI RV.INCO.'ERK') erkm;
list rv.inco.'ERE' ;
mess ' Sur E' (MINI RV.INCO.'ERE') erem;
SI ( (MINI RV.INCO.'ERE') > erem ) ;
* ERREUR 5 ;
FINSI ;
SI ( (MINI RV.INCO.'ERK') > erkm ) ;
* ERREUR 5 ;
FINSI ;
FIN BCFL ;
SI ( GRAPH ) ;
TRACE DOMTOT ;
TAB1 = TABLE ;
TAB1.1 = 'TIRC MARQ TRIA REGU' ;
TAB1.2 = 'TIRC MARQ TRIA REGU' ;
TAB1.3 = 'TIRL MARQ LOSA REGU' ;
TAB1.4 = 'TIRL MARQ LOSA REGU' ;
TAB1.5 = 'TIRM MARQ TRIB REGU' ;
TAB1.6 = 'TIRM MARQ TRIB REGU' ;
TAB1.7 = 'TIRR MARQ ETOI REGU' ;
TAB1.8 = 'TIRR MARQ ETOI REGU' ;
TAB1.'TITRE' = TABLE ;
REPETER BCL NCFL ;
TAB1.'TITRE'. (2*&bcl-1) = (CHAI 'Cv_ k : CFL=' (EXTR LCFL &bcl));
TAB1.'TITRE'. (2*&bcl) = (CHAI 'Cv_ e : CFL=' (EXTR LCFL &bcl));
FIN BCL ;
titl = chai 'Convergence Linf : Cfl = ' cfl ;
DESS EVOLER 'XBOR' 0. (2.*NITMA) 'YBOR' -15.0 0.0 'GRIL'
'TITR' titl LEGE TAB1 ;
TRACE RV.INCO.'KN' DOMTOT CNT 'TITR' 'KN' 20 ;
TRACE RV.INCO.'EN' DOMTOT CNT 'TITR' 'EN' 20 ;
TAB1 = TABLE ;
TAB1.1 = 'MARQ CROI ' ;
TAB1.2 = 'TIRC MARQ TRIA REGU' ;
TAB1.3 = 'TIRL MARQ CARR REGU' ;
TAB1.4 = 'TIRM MARQ LOSA REGU' ;
TAB1.5 = 'TIRM MARQ ETOI REGU' ;
TAB1.'TITRE' = TABLE ;
TAB1.'TITRE'. 1 = 'MOT' 'k : ANALYTIQUE' ;
REPETER BCL NCFL ;
TAB1.'TITRE'. (1+&bcl) = (CHAI 'k : CFL=' (EXTR LCFL &bcl));
FIN BCL ;
mess ' CFL=' CFL ;
titk = chai 'Grid Turbulence K : Cfl = ' cfl ;
tite = chai 'Grid Turbulence E : Cfl = ' cfl ;
DESS EVOLKA LEGE TAB1 'TITR' titk 'GRIL';
TAB1.'TITRE'. 1 = 'MOT' 'epsilon ANALYTIQUE' ;
REPETER BCL NCFL ;
TAB1.'TITRE'. (1+&bcl) = (CHAI 'e : CFL=' (EXTR LCFL &bcl));
FIN BCL ;
DESS EVOLEA LEGE TAB1 'TITR' tite 'GRIL';
FINSI ;
ERRK = KOPS RV.INCO.'KN' '-' KAN ;
ERRK = KOPS ERRK '*' ERRK ;
ERR2 = 0. ;
REPETER BLOC1 (NBNO (DOMA $DOMTOT MAILLAGE)) ;
P1 = (doma $DOMTOT SOMMET) POIN &BLOC1 ;
ERR1 = 'EXTR' ERRK 'SCAL' P1 ;
ERR2 = ERR2 + ERR1 ;
FIN BLOC1 ;
ERR2 = ERR2/(nbno (doma $DOMTOT MAILLAGE)) ;
ERR2 = ERR2 ** 0.5 ;
ERR2K = ERR2/Kinlet ;
MESSAGE 'ERREUR SUR K EN NORME L2 = ' ERR2K ;
ERRE = KOPS RV.INCO.'EN' '-' EAN ;
ERRE = KOPS ERRE '*' ERRE ;
ERR2 = 0. ;
REPETER BLOC1 (NBNO (DOMA $DOMTOT MAILLAGE)) ;
P1 = (DOMA $DOMTOT SOMMET) POIN &BLOC1 ;
ERR1 = 'EXTR' ERRE 'SCAL' P1 ;
ERR2 = ERR2 + ERR1 ;
FIN BLOC1 ;
ERR2 = ERR2/(NBNO (DOMA $DOMTOT MAILLAGE)) ;
ERR2 = ERR2 ** 0.5 ;
ERR2E = ERR2/Einlet ;
MESSAGE 'ERREUR SUR E EN NORME L2 = ' ERR2E ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer