* fichier : 15wedge.dgibi
******************************************************************
* CALCUL DE L'ECOULEMENT SUPERSONIQUE STATIONNAIRE DANS UN CANAL *
* AVEC RAMPE INCLINEE A 15deg. *
* POWELL AND DEZEEUW, AIAA-91-1542-CP *
* FORMULATION VF COMPRESSIBLE EXPLICITE/IMPLICITE *
* *
* H. PAILLERE/P. GALON TTMF AOUT 1997 *
* *
* MODIF BECCANTINI MARS 97 (reconstruction lineaure exacte) *
* MODIF BECCANTINI MARS 01 (implicite) *
* MODIF, 10/07/01, syntaxe de PENT changée *
******************************************************************
'MESSAGE' 'A mettre a jours' ;
'FIN' ;
TYEL = 'QUA4' ;
'OPTION' 'DIME' 2 'ELEM' TYEL 'ISOV' 'SULI'
'ECHO' 1 'TRAC' 'X' ;
*****************************************************
*****************************************************
***** PROCEDURE PROLIM *****
*****************************************************
*****************************************************
*
*
**** Cas Euler mono-espèce
*
'DEBPROC' PROLIM ;
'ARGUMENT' RV * TABLE ;
UN = RV . 'UN' ;
UNLIM = RV . 'UNLIM' ;
MAILLIM = 'EXTRAIRE' UNLIM 'MAILLAGE' ;
UNLEV = 'REDU' UN MAILLIM '*' (-1.0) ;
UN = UN '+' UNLEV '+' UNLIM ;
'FINPROC' UN ;
*********************************************************************
**** Procedure CALC *************************************************
*********************************************************************
*
* Cette procédure, utilisée comme un opérateur, calcule la différence
* entre deux pas de temps toutes les (RVX . 'FCALC') itérations.
* L'evolution de cette différence (erreur absolue) au cours du temps
* est conservée en vue d'un post-traitement.
*
*
'DEBP' CALC ;
'ARGU' RVX*'TABLE' ;
*
RV = RVX . 'EQEX' ;
UN = RV . 'UN' ;
NOMDEN = 'EXTRAIRE' (RV . 'LISTCONS') 1 ;
MOMMOM = 'EXTRAIRE' (RV . 'LISTCONS') ('LECT' 2 3) ;
RN = 'EXCO' UN NOMDEN ;
GN = 'EXCO' UN NOMMOM ('MOTS' 'UX' 'UY') ;
*
'SI' ( 'EXIS' RVX 'COMPT') ;
RVX . 'COMPT' = RVX . 'COMPT' + 1 ;
'SINON' ;
RVX . 'COMPT' = 1 ;
RVX . 'IT' = 'PROG' ;
RVX . 'TPS' = 'PROG' ;
RVX . 'ER1' = 'PROG' ;
RVX . 'ERINF' = 'PROG' ;
RVX . 'GN0' = 'COPI' GN ;
RVX . 'CFLX' = 'PROG' ;
RVX . 'CFLY' = 'PROG' ;
RVX . 'RESLINY' = 'PROG' ;
'FINSI' ;
*
DD = RVX . 'COMPT' ;
COMBIEN = RVX . 'FCALC' ;
NN = DD '/' COMBIEN ;
LO = (DD '-' (COMBIEN '*' NN)) 'EGA' 0 ;
'SI' LO ;
ERR = GN '-' (RVX . 'GN0') ;
ERRINF = 'MAXIMUM' ERR 'ABS' ;
ERR = ('XTX' ERR) '**' 0.5 ;
ERR = ERR '/' ('NBNO' ('EXTRAIRE' GN 'MAILLAGE')) ;
RVX . 'GN0' = 'COPI' GN ;
EL1 = ('LOG' (ERR + 1.0E-50))/('LOG' 10.) ;
ELINF = ('LOG' (ERRINF + 1.0E-50))/('LOG' 10.) ;
'MESSAGE' ;
'MESSAGE' 'Erreur de convergence' ;
'MESS' 'ITERATION ' (RVX . 'COMPT') ' LOG10 ERR_1(GN)' EL1 ;
'MESS' 'ITERATION ' (RVX . 'COMPT') ' LOG10 ERR_I(GN)' ELINF ;
'MESSAGE' ;
IT = 'PROG' RVX . 'COMPT' ;
ER1 = 'PROG' EL1 ;
ERINF = 'PROG' ELINF ;
TPS = 'PROG' (RV . 'PASDETPS' . 'TPS') ;
RVX . 'IT' = (RVX . 'IT') 'ET' IT ;
RVX . 'ER1' = (RVX . 'ER1') 'ET' ER1 ;
RVX . 'ERINF' = (RVX . 'ERINF') 'ET' ERINF ;
RVX . 'TPS' = (RVX . 'TPS') 'ET' TPS ;
'FINS' ;
*
* On peut augmunter le CFL (et MATTAB . 'RESID')
*
'SI' VRAI ;
COMBIEN = RVX . 'FCFL' ;
NN = DD '/' COMBIEN ;
LO = (DD '-' (COMBIEN '*' NN)) 'EGA' 0 ;
'SI' LO ;
CFL = PKONTAB . 'ALPHA' ;
CFL = CFL '*' 2. ;
CFL = 'MINIMUM' ('PROG' CFL (RVX . 'ALPMAX')) ;
PKONTAB . 'ALPHA' = CFL ;
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'CFL =' ((PKONTAB . 'ALPHA') '/' 2)) ;
* We decrease the linear iterations error
MATTAB . 'RESID' = 'MAXIMUM' ('PROG'
(MATTAB . 'RESID' '/' 3.1623) 1.0D-14) ;
'MESSAGE' ('CHAINE' 'MATTAB . RESID = ' MATTAB . 'RESID') ;
'FINSI' ;
RVX . 'CFLX' = RVX . 'CFLX' 'ET' ('PROG' DD) ;
RVX . 'CFLY' = RVX . 'CFLY' 'ET' ('PROG'
((PKONTAB . 'ALPHA') '/' 2)) ;
RVX . 'RESLINY' = (RVX . 'RESLINY') 'ET'
('PROG' MATTAB . 'RESID') ;
'FINSI' ;
IRESU IJACO ='KOPS' 'MATRIK' ;
IALPDT = 'PROG' ;
'MENAGE' ;
'RESPRO' IRESU IALPDT ;
'FINP' ;
************************************************************************
************************************************************************
***************** FIN PARTIE PROCEDURES ********************************
************************************************************************
************************************************************************
GRAPH = VRAI ;
* GRAPH = FAUX ;
***********************************************************
*********************** MAILLAGE **************************
***********************************************************
NRAFF = 8 ;
NX1 = 1 '*' NRAFF ;
NX2 = 1 '*' NRAFF ;
NX3 = 4 '*' NRAFF ;
NY = 2 '*' NRAFF ;
DY = 2.0 '/' NY ;
H1 = 0.268 ;
A1 = (1. '-' DY) 0.0 ;
A2 = 1.0 0.0 ;
A3 = 2.0 0.0 ;
A4 = 3.0 H1 ;
A5 = 7.0 H1 ;
A6 = (7.0 '+' DY) H1 ;
B1 = (1.0 '-' DY) 2.0 ;
B2 = 1.0 2.0 ;
B3 = 2.0 2.0 ;
B4 = 3.0 2.0 ;
B5 = 7.0 2.0 ;
B6 = (7.0 '+' DY) 2.0 ;
BAS1 = A2 'DROI' NX1 A3 ;
BAS2 = A3 'DROI' NX2 A4 ;
BAS3 = A4 'DROI' NX3 A5 ;
HAU3 = B5 'DROI' NX3 B4 ;
HAU2 = B4 'DROI' NX2 B3 ;
HAU1 = B3 'DROI' NX1 B2 ;
COTEG = B2 'DROI' NY A2 ;
COTED = A5 'DROI' NY B5 ;
FRONTG = 'DALLER' (A1 'DROI' 1 A2) ('INVE' COTEG)
(B2 'DROI' 1 B1) (B1 'DROI' NY A1) ;
FRONTG = FRONTG 'COUL' ROUG ;
FRONTD = 'DALLER' (A5 'DROI' 1 A6) (A6 'DROI' NY B6)
(B6 'DROI' 1 B5) ('INVE' COTED) ;
FRONTD = FRONTD 'COUL' VERT ;
FRONT = FRONTG 'ET' FRONTD ;
'ELIMINATION' FRONT 1D-6;
DOMINT = 'DALLER' (BAS1 ET BAS2 ET BAS3) COTED
(HAU3 ET HAU2 ET HAU1) COTEG 'PLAN' ;
DOMTOT = DOMINT ET FRONT ;
'ELIMINATION' DOMTOT (1D-6 '/' NRAFF) ;
$DOMTOT = 'DOMA' DOMTOT ;
$DOMINT = 'DOMA' DOMINT 'INCL' $DOMTOT (1D-6 '/' NRAFF) ;
$FRONT = 'DOMA' FRONT 'INCL' $DOMTOT (1D-6 '/' NRAFF) ;
'DOMA' $DOMINT 'XXNORMAF' ;
CHPONOR = $DOMINT . 'XXNORMAF' ;
*****************************************************************
* CONDITIONS INITIALES ET LIMITES : ECOULEMENT SUPERSONIQUE M=2 *
* Adimensionalisation cf. These HP. p. 152 *
*****************************************************************
gamscal = 1.4 ;
GAMN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gamscal;
Minf = 2.0D0 ;
alpha = 0.D0 ;
'SI' ('NEG' ALPHA 0.0) ;
'MESSAGE' 'Attention aux C.L.' ;
'ERREUR' 21 ;
'FINSI' ;
rhoinf = 1.0D0 ;
rhouinf = 'COS' alpha ;
rhovinf = 'SIN' alpha ;
rhoEtinf = 2.0D0 '+' (gamscal '*' Minf '*' Minf
'*' (gamscal '-' 1.0D0)) ;
rhoEtinf = rhoEtinf '/' (2.D0 '*' gamscal '*' Minf
'*' Minf '*' (gamscal '-' 1.0D0)) ;
RN0 = 'MANUEL' 'CHPO' ($DOMTOT . 'CENTRE') 1 'SCAL' rhoinf ;
GN0 = 'MANUEL' 'CHPO' ($DOMTOT . 'CENTRE') 2 'UX' rhouinf
'UY' rhovinf ;
RETN0 = 'MANUEL' 'CHPO' ($DOMTOT . 'CENTRE') 1 'SCAL' rhoEtinf ;
VN0 PN0 = 'PRIM' 'PERFMONO' RN0 GN0 RETN0 GAMN ;
VN2 = 'PSCAL' VN0 VN0 ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
C2 = GAMN '*' ( PN0 '/' RN0) ;
MACH2 = VN2 '/' C2;
MACHN0 = MACH2 '**' 0.5;
MOD1 = 'MODELISER' ($DOMTOT . 'MAILLAGE' ) 'THERMIQUE';
'SI' GRAPH ;
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN0 ;
CHM_PN = 'KCHA' $DOMTOT 'CHAM' PN0 ;
CHM_VN = 'KCHA' $DOMTOT 'CHAM' VN0 ;
CHM_MACH = 'KCHA' $DOMTOT 'CHAM' MACHN0 ;
'TRACER' CHM_RN MOD1
'TITR' ('CHAINE' 'RN at t=' 0.0);
'TRACER' CHM_PN MOD1
'TITR' ('CHAINE' 'PN at t=' 0.0);
'TRACER' CHM_MACH MOD1
'TITR' ('CHAINE' 'MACHN at t=' 0.0);
'TRACER' CHM_VN MOD1
'TITR' ('CHAINE' 'VN at t=' 0.0);
'FINSI' ;
NOMDEN = 'RN' ;
NOMMOM = 'MOTS' 'RUX' 'RUY' ;
NOMENE = 'RETN' ;
LISTCONS = ('MOTS' NOMDEN) 'ET' NOMMOM 'ET' ('MOTS' NOMENE) ;
UN =
('NOMC' NOMDEN RN0 'NATURE' 'DISCRET') 'ET'
('EXCO' ('MOTS' 'UX' 'UY') GN0 NOMMOM
'NATURE' 'DISCRET') 'ET'
('NOMC' NOMENE RETN0 'NATURE' 'DISCRET') ;
*****************************************************
*****************************************************
*****************************************************
*****************************************************
**************** La table RV **********************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
RV = 'TABLE' ;
* Constant time step
* RV . 'DTIMP' = 1.191 ;
*
RV . 'UN' = UN ;
* Primitive formulation
RV . 'FORMPRIM' = FAUX ;
RV . 'LISTCONS' = LISTCONS ;
RV . 'LISTPRIM' = LISTCONS ;
**** Le domaine interne
RV . 'DOMINT' = $DOMINT . 'CENTRE' ;
*
RV . 'CLIM' = VRAI ;
RV . 'MAIFAN' = $FRONT . 'CENTRE' ;
*
RV . 'INST' = VRAI ;
RV . 'CONS' = FAUX ;
*
'SI' (RV . 'INST') ;
RV . 'TFINAL' = 10.0D5 ;
RV . 'NITMAEX' = 100 ;
'FINSI' ;
*
RV . 'NITMAIN' = 1 ;
RV . 'NITMIIN' = 1 ;
RV . 'EPSINT' = 3.0D-10 ;
*
RV . 'FEXT' = 20 ;
RV . 'FINT' = 1 ;
*
RV . 'LISTOPER' = 'MOTS' 'CALC' 'PKON' ;
RV . '1CALC' = 'TABLE' ;
CALCTAB = RV . '1CALC' ;
RV . '2PKON' = 'TABLE' ;
PKONTAB = RV . '2PKON' ;
*
RV . 'MATIDE' = 'KOPS' 'MATIDE' LISTCONS ($DOMTOT . 'CENTRE') ;
*
**** Parametres de la procedure PROLIM
*
RV . 'UNLIM' = 'REDU' (RV . 'UN') (RV . 'MAIFAN') ;
* Gas model (Procedure EXEXPRIM, EXEXCONV)
RV . 'GAMMA' = GAMN ;
* Bas Mach (Procedures EXEXDT, PKON)
RV . 'ORDTPS' = 1 ;
RV . 'PRECTEMP' = FAUX ;
* RV . DUALTEXP = VRAI ;
* RV . DUALTIMP = VRAI ;
RV . CFLDTAU = 1.0 ;
RV . 'DIAM' = ('DOMA' $DOMTOT 'DIAMIN') '*' 1.0 ;
* RV . 'CO1' = 'MANUEL' 'CHPO' ($DOMTOT . 'CENTRE')
* 1 'SCAL' (0.1 * U_G) '+'
* ('MANUEL' 'CHPO' ($FRONT . 'CENTRE')
* 1 'SCAL' (0. * U_G)) ;
* RV . 'CO2' = 'MANUEL' 'CHPO' ($DOMTOT . 'CENTRE')
* 1 'SCAL' (0. * U_G) ;
*****************
** CALCUL *******
*****************
CALCTAB . 'EQEX' = RV ;
CALCTAB . 'FCALC' = 5 ;
CALCTAB . 'IMPL' = FAUX ;
* Doubling frequence for CFL
CALCTAB . 'FCFL' = 5 ;
* Maximum of CFL
CALCTAB . 'ALPMAX' = 10000.0 ;
*****************
* PKON **********
*****************
PKONTAB . 'EQEX' = RV ;
PKONTAB . 'GAZ' = 'PERFMONO' ;
PKONTAB . 'GAMN' = GAMN ;
PKONTAB . 'DOMA' = $DOMTOT ;
PKONTAB . 'LISTCONS' = LISTCONS ;
* PKONTAB . 'LISTPRIM' = LISTPRIM ;
PKONTAB . 'METHODE' = 'VLH' ;
*
*
PKONTAB . 'IMPL' = VRAI ;
'SI' (PKONTAB . 'IMPL') ;
PKONTAB . 'TYPEJACO' = 'VLH' ;
'FINSI' ;
*
PKONTAB . 'ORDREESP' = 2 ;
'SI' ((PKONTAB . 'ORDREESP') 'EGA' 2) ;
PKONTAB . 'LIMITEUR' = 'LIMITEUR' ;
CHPVID MATVID = 'KOPS' 'MATRIK' ;
PKONTAB . 'VLIM' = CHPVID ;
CACCA PIPI GRADSCAL = 'PENTE' $DOMTOT 'CENTRE' 'EULESCAL'
'NOLIMITE' ('EXCO' NOMDEN (RV . 'UN')) ;
CACCA PIPI GRADVEC = 'PENTE' $DOMTOT 'CENTRE' 'EULEVECT'
'NOLIMITE' ('EXCO' NOMMOM (RV . 'UN')
('MOTS' 'UX' 'UY')) 'CLIM' (PKONTAB . 'VLIM') ;
PKONTAB . 'GRADRN' = GRADSCAL ;
PKONTAB . 'GRADPN' = GRADSCAL ;
PKONTAB . 'GRADVN' = GRADVEC ;
'FINSI' ;
*
ELP1 = 'MANUEL' 'POI1' (123.56 124.87) ;
PKONTAB . 'FACELIM' = 'DIFF' ELP1 ELP1 ;
PKONTAB . 'ALPHA' = 10. ;
***********************************
* Inversion de la matrice *********
***********************************
* RV . 'LINEHIST' = historique de convergence des iterations internes
* RV . 'MATUPDAT' updating
* RV . 'MUPINI' we update de matrix if the external iterations
* are less than (RV . 'MUPINI')
* RV . 'MUPEXT' external iterations updating
* RV . 'MUPINT' internal iterations updating
* RV . 'MATACC' Broyden method?
* RV . 'MUPLIN' We update the matrix if in the previous
* internal iteration the number of linear
* iterations to solve the system were bigger
* than (RV . 'MUPLIN')
RV . 'LINEHIST' = 'PROG' ;
RV . 'MATUPDAT' = VRAI ;
'SI' ('NON' ( RV . 'MATUPDAT')) ;
RV . 'MUPEXT' = 4000 ;
RV . 'MUPINT' = 200 ;
RV . 'MATACC' = FAUX ;
RV . 'MUPINI'= 0 ;
RV . 'MUPLIN' = 80 ;
'FINSI' ;
RV . 'MATINV' = 'TABLE' 'METHINV' ;
MATTAB = RV . 'MATINV' ;
MATTAB . 'TYPINV' = 5 ;
MATTAB . 'IMPINV' = 0 ;
*
* Matrice pour assurer que la matrice à inverser est correctement
* assemblé
*
* MATTAB . 'MATASS' definie en EXEXIM
* MATTAB . 'MAPREC' "
* MATTAB . 'XINIT' "
*
* Methode de numerotation de DDL
MATTAB . 'TYRENU' = 'RIEN' ;
MATTAB . 'PCMLAG' = 'APR2' ;
MATTAB . 'GMRESTRT' = 500 ;
* ILU 3
* ILUT (dual truncation) 5
* Dans le cas ILUT il faut definir
* ILUTLFIL
* ILUTDTOL
MATTAB . 'PRECOND' = 3 ;
* MATTAB . 'ILUTLFIL' = 32. ;
* MATTAB . 'ILUTDTOL' = 0. ;
MATTAB . 'NITMAX' = 1000 ;
MATTAB . 'RESID' = 1.D-3 ;
* bcgstab
MATTAB . 'BCGSBTOL' = 1.0d-40 ;
'OPTION' DONN 5 ;
'TEMPS' 'ZERO' ;
EXEXIM RV ;
'TEMPS' ;
*
**** Les variables conservatives
*
'SI' (RV . 'INST') ;
TFIN = RV . 'PASDETPS' . 'TPS' ;
'SINO' ;
TFIN = 'INF' ;
'FINSI' ;
RN = 'EXCO' 'RN' (RV . 'UN') ;
GN = 'EXCO' ('MOTS' 'RUX' 'RUY') (RV . 'UN') ('MOTS' 'UX' 'UY') ;
RETN = 'EXCO' 'RETN' (RV . 'UN') ;
*
**** Les variables primitives
*
VN PN = 'PRIM' 'PERFMONO'
RN GN RETN GAMN ;
CSON2 = (GAMN '*' PN) '/' RN ;
VN2 = 'PSCAL' VN VN ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
MACHN2 = VN2 '/' CSON2 ;
MACHN = MACHN2 '**' 0.5 ;
*
*** GRAPHIQUE DES SOLUTIONS
*
MOD1 = 'MODELISER' $DOMINT 'THERMIQUE';
'SI' GRAPH ;
CHM_RN ='KCHA' $DOMINT 'CHAM' ('REDU' RN ($DOMINT . 'CENTRE')) ;
CHM_PN ='KCHA' $DOMINT 'CHAM' ('REDU' PN ($DOMINT . 'CENTRE')) ;
CHM_VN = 'KCHA' $DOMINT 'CHAM' ('REDU' VN ($DOMINT . 'CENTRE')) ;
CHM_MACH= 'KCHA' $DOMINT 'CHAM' ('REDU' MACHN ($DOMINT . 'CENTRE'));
'TRACER' CHM_RN MOD1 ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'METO =' METO ' NBEL =' NTOTEL
' RN at t=' TFIN ' OE =' 2);
'TRACER' CHM_PN MOD1 ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'METO =' METO ' NBEL =' NTOTEL
' PN at t=' TFIN ' OE =' 2);
'TRACER' CHM_MACH MOD1 ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'METO =' METO ' NBEL =' NTOTEL
' MACHN at t=' TFIN ' OE =' 2);
'TRACER' CHM_VN MOD1 ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'METO =' METO ' NBEL =' NTOTEL
' VN at t=' TFIN ' OE =' 2);
RNC = 'ELNO' $DOMINT ('REDU' RN ($DOMINT . 'CENTRE'));
'TRACER' DOMINT RNC ('CONTOUR' DOMINT)
'TITR' ('CHAINE' 'METO =' METO ' NBEL =' NTOTEL
' RN at t=' TFIN ' OE =' 2) 15 ;
'FINSI' ;
*
*** Test
*
PCEN1 = ($DOMTOT . 'CENTRE') 'POIN' 'PROC' (2.5 (H1 '/' 2.0)) ;
'SI' GRAPH ;
'TRACER' (DOMTOT 'ET' PCEN1) 'TITRE'
'Test 1: solution derrier le premier choc' ;
'FINSI' ;
ro1 = 'EXTRAIRE' RN 'SCAL' PCEN1 ;
erro = ((ro1 '-' 1.7289) 'ABS') '/' 1.7289 ;
'SI' (erro > 1.0D-1) ;
'ERREUR' 5 ;
'FINSI' ;
p1 = 'EXTRAIRE' PN 'SCAL' PCEN1 ;
erro = ((p1 '-' 0.3919) 'ABS') '/' 0.3919 ;
'SI' (erro > 1.0D-1) ;
'ERREUR' 5 ;
'FINSI' ;
m1 = 'EXTRAIRE' MACHN 'SCAL' PCEN1 ;
erro = ((m1 '-' 1.44572) 'ABS') '/' 1.44572 ;
'SI' (erro > 1.0D-1) ;
'ERREUR' 5 ;
'FINSI' ;
* Derrier le mach stem on a le mach min
* MMIN = 0.577350
* Il faut avoir un maillage assez fin pour l'avoir
MMIN = 'MINIMUM' MACHN ;
'MESSAGE' ('CHAINE' 'MMIN_ex =' 0.577350) ;
'MESSAGE' ('CHAINE' 'MMIN =' MMIN) ;
'SI' (MMIN > 1.5) ;
'ERREUR' 5 ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer