Download tube1D_deto_C2H2.dgibi
* fichier : tube1D_deto_C2H2.dgibi
****************************************************************
**** APPROCHE VF "Cell-Centred Formulation" pour la ****
**** combustion. ****
**** MODELE RDEM ****
**** ****
**** OPERATEURS 'PRIM', PRET, KONV ****
**** ****
**** PROPAGATION D'UNE CJDT DANS UN TUBE ****
**** Cas de l'acetylene ****
**** A. BECCANTINI, SFME/LTMF, 30.07.2010 ****
****************************************************************
*
* The file is divided into 6 parts
*
* 0) exact solution and procedures
* 1) mesh
* 2) initial conditions and gas properties
* 3) parameters used in the computation
* 4) the main part (where the computation is realised)
* 5) the post-treatment
*
'OPTION' 'ECHO' 1 'DIME' 2
'TRAC' 'PSC' ;
GRAPH = VRAI ;
* GRAPH = FAUX ;
* Upwind scheme
METO = 'VLH' ;
* METO = 'SS' ;
* METO = 'AUSMPUP' ;
*
****************************************************************
****************************************************************
************* EXACT SOLUTION AND PROCEDURES ********************
****************************************************************
****************************************************************
*
***** BEGIN EXACT SOLUTION
*
LEXACT = 'PROG'
0.0 .9300690945091593 3.979039320256561E-12 697580.6326340148
2563.02106222187 1.0 288819.5707231797 268622.1566190319
48.64683385104058 .9300690945091593 3.979039320256561E-12
697580.6326340148 2563.02106222187 1.0 288819.5707231797
268622.1566190319 97.29366770208105 .9300690945091593
3.979039320256561E-12 697580.6326340148 2563.02106222187 1.0
288819.5707231797 268622.1566190319
;
LEXACT = LEXACT 'ET' ('PROG'
145.9405015531216 .9300690945091593 3.979039320256561E-12
697580.6326340148 2563.02106222187 1.0 288819.5707231797
268622.1566190319 194.5873354041621 .9300690945091593
3.979039320256561E-12 697580.6326340148 2563.02106222187
1.0 288819.5707231797 268622.1566190319 ) ;
LEXACT = LEXACT 'ET' ('PROG'
243.2341692552026
.9300690945091593 3.979039320256561E-12 697580.6326340148
2563.02106222187 1.0 288819.5707231797 268622.1566190319
291.8810031062432 .9300690945091593 3.979039320256561E-12
697580.6326340148 2563.02106222187 1.0 288819.5707231797
268622.1566190319
) ;
LEXACT = LEXACT 'ET' ('PROG'
340.5278369572836 .9300690945091593 3.979039320256561E-12
697580.6326340148 2563.02106222187 1.0 288819.5707231797
268622.1566190319 389.1746708083242 .9300690945091593
3.979039320256561E-12 697580.6326340148 2563.02106222187 1.0
288819.5707231797 268622.1566190319) ;
LEXACT = LEXACT 'ET' ('PROG'
437.8215046593647 .9300690945091593 3.979039320256561E-12
697580.6326340148 2563.02106222187 1.0 288819.5707231797
268622.1566190319 486.4683385104051 .9300690945091593
3.979039320256561E-12 697580.6326340148 2563.02106222187
1.0 288819.5707231797 268622.1566190319) ;
LEXACT = LEXACT 'ET' ('PROG'
535.1151723614457 .9300690945091593 3.979039320256561E-12
697580.6326340148 2563.02106222187 1.0 288819.5707231797
268622.1566190319 583.7620062124862 .9300690945091593
3.979039320256561E-12 697580.6326340148 2563.02106222187 1.0
288819.5707231797 268622.1566190319) ;
LEXACT = LEXACT 'ET' ('PROG'
632.4088400635267 .9300690945091593 3.979039320256561E-12
697580.6326340148 2563.02106222187 1.0 288819.5707231797
268622.1566190319 681.0556739145673 .9300690945091593
3.979039320256561E-12 697580.6326340148 2563.02106222187 1.0
288819.5707231797 268622.1566190319) ;
LEXACT = LEXACT 'ET' ('PROG'
729.7025077656077 .9300690945091593 3.979039320256561E-12
697580.6326340148 2563.02106222187 1.0 288819.5707231797
268622.1566190319 778.3493416166482 .9300690945091593
3.979039320256561E-12 697580.6326340148 2563.02106222187
1.0 288819.5707231797 268622.1566190319) ;
LEXACT = LEXACT 'ET' ('PROG'
826.9961754676888 .9300690945091593 3.979039320256561E-12
697580.6326340148 2563.02106222187 1.0 288819.5707231797
268622.1566190319 875.6430093187292 .9300690945091593
3.979039320256561E-12 697580.6326340148 2563.02106222187 1.0
288819.5707231797 268622.1566190319) ;
LEXACT = LEXACT 'ET' ('PROG'
924.2898431697698 .9300690945091593 3.979039320256561E-12
697580.6326340148 2563.02106222187 1.0 288819.5707231797
268622.1566190319 972.9366770208103 .9300690945091593
3.979039320256561E-12 697580.6326340148 2563.02106222187 1.0
288819.5707231797 268622.1566190319) ;
LEXACT = LEXACT 'ET' ('PROG' 972.9366770208103 .9300690945091593
3.979039320256561E-12 697580.6326340148 2563.02106222187 1.0
288819.5707231797 268622.1566190319 1020.563696234262 .9711992885935007
42.21849439643984 736732.9074055142 2592.236892185716 1.0
321460.0704866028 313067.3251900754 1067.969151870608 1.01370750996012
84.2438857700123 777645.6009989232 2621.452722149562 1.0
354138.8712603799 362590.3908306613) ;
LEXACT = LEXACT 'ET' ('PROG' 1115.155623979032 1.057624905695079
126.0783005920503 820378.2182800922 2650.668552113408 1.0
386854.8055661046 417553.1413927785 1162.125645554223 1.102982969882543
167.7238330059224 864991.6014201408 2679.884382077254 1.0
419606.7429636912 478333.2552105426 1208.881704582817 1.149813544020885
209.1825462777236 911547.9450587952 2709.1002120411 1.0
452393.5895046745 545324.6661907386) ;
LEXACT = LEXACT 'ET' ('PROG' 1255.426245966713 1.198148817580572
250.4564741516483 960110.8116026787 2738.316042004945 1.0
485214.2871855269 618937.9308335221 1301.761673331164 1.248021328707566
291.5476221162355 1010745.146670078 2767.531871968791 1.0
518067.8134009624 699600.5972715921 1347.890350724929 1.299463965075698
332.4579685872177 1063517.29469371 2796.747701932637 1.0
550953.1803972456 787757.5764255146) ;
LEXACT = LEXACT 'ET' ('PROG' 1393.814604219309 1.352509964891132
373.1894660122936 1118495.014692988 2825.963531896483 1.0
583869.4347255011 883871.5153818108 1439.536723412327 1.407192918051686
413.7440419027419 1175747.496227202 2855.179361860329 1.0
616815.6566950218 988423.173109015 1485.058962843933 1.463546767463509
454.1235997964707 1235345.375540995 2884.395191824175 1.0
649790.9598265779 1101911.798635182) ;
LEXACT = LEXACT 'ET' ('PROG' 1530.383543327666 1.521605810517206
494.3300201567479 1297360.751913322 2913.611021788021 1.0
682794.4903057241 1224855.511818089 1575.512653203853 1.581404700725235
534.3651612105829 1361867.204221023 2942.826851751867 1.0
715825.426436109 1357791.686846843 1620.448449519048 1.642978449522024
574.2308597304395 1428939.807727903 2972.042681715713 1.0
748882.9780927855 1501277.33862052) ;
LEXACT = LEXACT 'ET' ('PROG' 1665.193059136131 1.70636242822789
613.9289317627123 1498655.151110031 3001.258511679559 1.0
781966.3861755147 1655889.512155949 1709.748579779141 1.771592370177581
653.4611733061729 1571091.353727757 3030.474341643404 1.0
815074.9220620818 1822225.67518284 1754.117081016696 1.838704373013818
692.8293609433672 1646328.083154666 3059.69017160725 1.0
848207.8870615955 2000904.114089844) ;
LEXACT = LEXACT 'ET' ('PROG' 1798.300605187562 1.907734901145937
732.0352524277539 1724446.572973421 3088.906001571096 1.0
881364.6118678055 2192564.333390289 1842.301168271707 1.978720788373339
771.0805872291835 1805529.640848111 3118.121831534942 1.0
914544.4560124031 2397867.458880638 1886.120760709965 2.05169924067314
809.9670870401542 1889661.706882407 3147.337661498788 1.0
947746.8073183382 2617496.644668788) ;
LEXACT = LEXACT 'ET' ('PROG' 1886.120760709965 1.170626887265085 0.0
101300.0 298 0.0 582562.3435327417 681963.1428475864 1904.981968317065
1.170626887265085 0.0 101300.0 298 0.0 582562.3435327417
681963.1428475864 1923.843175924165 1.170626887265085 0.0 101300.0 298
0.0 582562.3435327417 681963.1428475864 1942.704383531264
1.170626887265085 0.0 101300.0 298 0.0 582562.3435327417
681963.1428475864) ;
LEXACT = LEXACT 'ET' ('PROG' 1961.565591138364 1.170626887265085 0.0
101300.0 298 0.0 582562.3435327417 681963.1428475864 1980.426798745464
1.170626887265085 0.0 101300.0 298 0.0 582562.3435327417
681963.1428475864 1999.288006352563 1.170626887265085 0.0 101300.0 298
0.0 582562.3435327417 681963.1428475864 2018.149213959663
1.170626887265085 0.0 101300.0 298 0.0 582562.3435327417
681963.1428475864) ;
LEXACT = LEXACT 'ET' ('PROG' 2037.010421566762 1.170626887265085 0.0
101300.0 298 0.0 582562.3435327417 681963.1428475864 2055.871629173862
1.170626887265085 0.0 101300.0 298 0.0 582562.3435327417
681963.1428475864 2074.732836780962 1.170626887265085 0.0 101300.0 298
0.0 582562.3435327417 681963.1428475864 2093.594044388061
1.170626887265085 0.0 101300.0 298 0.0 582562.3435327417
681963.1428475864 ) ;
LEXACT = LEXACT 'ET' ('PROG'
2112.455251995161 1.170626887265085 0.0 101300.0 298 0.0
582562.3435327417 681963.1428475864 2131.316459602261 1.170626887265085
0.0 101300.0 298 0.0 582562.3435327417 681963.1428475864
2150.17766720936 1.170626887265085 0.0 101300.0 298 0.0
582562.3435327417 681963.1428475864 2169.03887481646 1.170626887265085
0.0 101300.0 298 0.0 582562.3435327417 681963.1428475864) ;
LEXACT = LEXACT 'ET' ('PROG'
2187.90008242356 1.170626887265085 0.0 101300.0 298 0.0
582562.3435327417 681963.1428475864 2206.76129003066 1.170626887265085
0.0 101300.0 298 0.0 582562.3435327417 681963.1428475864
2225.622497637759 1.170626887265085 0.0 101300.0 298 0.0
582562.3435327417 681963.1428475864 ) ;
LEXACT = LEXACT 'ET' ('PROG'
2244.483705244859 1.170626887265085 0.0 101300.0 298 0.0
582562.3435327417 681963.1428475864 2263.344912851958 1.170626887265085
0.0 101300.0 298 0.0 582562.3435327417 681963.1428475864) ;
LEXACT = LEXACT 'ET' ('PROG'
10000.0 1.170626887265085 0.0 101300.0 298 0.0
582562.3435327417 681963.1428475864) ;
aa = 'DIME' LEXACT ;
LEX_XOT = 'LECT' 1 PAS 8 (aa '-' 7) ;
LEX_XOT = 'EXTRAIRE' LEXACT LEX_XOT ;
LEX_RHO = 'LECT' 2 PAS 8 (aa '-' 6) ;
LEX_RHO = 'EXTRAIRE' LEXACT LEX_RHO ;
LEX_U = 'LECT' 3 PAS 8 (aa '-' 5) ;
LEX_U = 'EXTRAIRE' LEXACT LEX_U ;
LEX_P = 'LECT' 4 PAS 8 (aa '-' 4) ;
LEX_P = 'EXTRAIRE' LEXACT LEX_P ;
LEX_T = 'LECT' 5 PAS 8 (aa '-' 3) ;
LEX_T = 'EXTRAIRE' LEXACT LEX_T ;
LEX_CSI = 'LECT' 6 PAS 8 (aa '-' 2) ;
LEX_CSI = 'EXTRAIRE' LEXACT LEX_CSI ;
LEX_EIN = 'LECT' 7 PAS 8 (aa '-' 1) ;
LEX_EIN = 'EXTRAIRE' LEXACT LEX_EIN ;
LEX_RET = 'LECT' 8 PAS 8 (aa '-' 0) ;
LEX_RET = 'EXTRAIRE' LEXACT LEX_RET ;
*
****END EXACT SOLUTION
*
******************
**** PRIMCONS ****
******************
*
* From primitive to conservative variables
*
'DEBPROC' PRIMCONS ;
'ARGUMENT' PGAS*TABLE TN1*'CHPOINT ' TN2*'CHPOINT '
PN1*'CHPOINT ' PN2*'CHPOINT ' VN1*'CHPOINT ' VN2*'CHPOINT ' ;
*
* ETHER = int_0^T cv(T') dT' T < TMAX
* = int_0^TMAX cv(T') dT' '+'
* cv(TMAX) T >= TMAX
ESP1 = 'EXTRAIRE' (PGAS . 'SPECIES') 1 ;
* DY1 = y_i - y_f for the species 1
DY1 = (('EXTRAIRE' (PGAS . 'MASSFRA') 1) '-'
('EXTRAIRE' (PGAS . 'MASSFRA') 2)) ;
COEF1 = ('EXTRAIRE' (PGAS . 'CHEMCOEF') 1) '*'
(PGAS . ESP1 . 'W') ;
YFINPH1 = 1.0 ;
YFINPH2 = 1.0 ;
'SI' (COEF1 > 0) ;
YPH2 = 'EXTRAIRE' (PGAS . 'MASSFRA') 2 ;
YPH1 = YPH2 '+' DY1 ;
'SINON' ;
YPH1 = 'EXTRAIRE' (PGAS . 'MASSFRA') 2 ;
YPH2 = YPH1 '-' DY1 ;
'FINSI' ;
YFINPH1 = YFINPH1 '-' YPH1 ;
YFINPH2 = YFINPH2 '-' YPH2 ;
PRYPH1 = 'PROG' YPH1 ;
PRYPH2 = 'PROG' YPH2 ;
'REPETER' BLESP (('DIME' (PGAS . 'SPECIES')) '-' 2) ;
ESP = 'EXTRAIRE' (PGAS . 'SPECIES') (&BLESP '+' 1) ;
COEF = ('EXTRAIRE' (PGAS . 'CHEMCOEF') (&BLESP '+' 1))
'*' (PGAS . ESP . 'W') ;
DY = (DY1 * (COEF '/' COEF1)) ;
'SI' (COEF > 0) ;
YPH2 = 'EXTRAIRE' (PGAS . 'MASSFRA') (&BLESP '+' 2) ;
YPH1 = YPH2 '+' DY ;
'SINON' ;
YPH1 = 'EXTRAIRE' (PGAS . 'MASSFRA') (&BLESP '+' 2) ;
YPH2 = YPH1 '-' DY ;
'FINSI' ;
PRYPH1 = PRYPH1 'ET' ('PROG' YPH1) ;
PRYPH2 = PRYPH2 'ET' ('PROG' YPH2) ;
YFINPH1 = YFINPH1 '-' YPH1 ;
YFINPH2 = YFINPH2 '-' YPH2 ;
'FIN' BLESP ;
PRYPH1 = PRYPH1 'ET' ('PROG' YFINPH1) ;
PRYPH2 = PRYPH2 'ET' ('PROG' YFINPH2) ;
'LISTE' PRYPH1 ;
'LISTE' PRYPH2 ;
*
TMAX = (PGAS . 'TMAX') ;
* TCAL1 = MIN TN1, TMAX
TCAL1 = 0.5D0 '*' ((TMAX '+' TN1) '-' ('ABS' (TN1 '-' TMAX))) ;
DTN1 = TN1 '-' TCAL1 ;
* TCAL1 = MIN TN1, TMAX
TCAL2 = 0.5D0 '*' ((TMAX '+' TN2) '-' ('ABS' (TN2 '-' TMAX))) ;
DTN2 = TN2 '-' TCAL2 ;
*
* Internal energy (J/kg in SI)
*
ETHER1 = 0.0 ;
CV1 = 0.0 ;
ETHER2 = 0.0 ;
CV2 = 0.0 ;
FUNTN1 = 1.0 ;
FUNTN2 = 1.0 ;
'REPETER' BLPO ((PGAS . 'NORD') '+' 1) ;
'REPETER' BLESP ('DIME' (PGAS . 'SPECIES')) ;
ESP = 'EXTRAIRE' (PGAS . 'SPECIES') &BLESP ;
YCEL1 = 'EXTRAIRE' PRYPH1 &BLESP ;
YCEL2 = 'EXTRAIRE' PRYPH2 &BLESP ;
AA = 'EXTRAIRE' (PGAS . ESP . 'A') &BLPO ;
DCV1 = (AA * YCEL1 * FUNTN1) ;
DCV2 = (AA * YCEL2 * FUNTN2) ;
CV1 = CV1 '+' DCV1 ;
CV2 = CV2 '+' DCV2 ;
ETHER1 = ETHER1 '+' (DCV1 * TCAL1 '/' (&BLPO)) ;
ETHER2 = ETHER2 '+' (DCV2 * TCAL2 '/' (&BLPO)) ;
'FIN' BLESP ;
FUNTN1 = FUNTN1 '*' TCAL1 ;
FUNTN2 = FUNTN2 '*' TCAL2 ;
'FIN' BLPO ;
ETHER1 = ETHER1 '+' (CV1 '*' DTN1) ;
ETHER2 = ETHER2 '+' (CV2 '*' DTN2) ;
*
* Formation energy/enthalpy (J/kg in SI) and gas constant (J/kg/K)
*
EFORM1 = 0.0 ;
EFORM2 = 0.0 ;
RGAS1 = 0.0 ;
RGAS2 = 0.0 ;
'REPETER' BLESP ('DIME' (PGAS . 'SPECIES')) ;
ESP = 'EXTRAIRE' (PGAS . 'SPECIES') &BLESP ;
YCEL1 = 'EXTRAIRE' PRYPH1 &BLESP ;
YCEL2 = 'EXTRAIRE' PRYPH2 &BLESP ;
EFORM1 = EFORM1 '+' (YCEL1 * (PGAS . ESP . 'H0K')) ;
EFORM2 = EFORM2 '+' (YCEL2 * (PGAS . ESP . 'H0K')) ;
RGAS1 = RGAS1 '+' (YCEL1 * (PGAS . 'RUNIV') '/'
(PGAS . ESP . 'W')) ;
RGAS2 = RGAS2 '+' (YCEL2 * (PGAS . 'RUNIV') '/'
(PGAS . ESP . 'W')) ;
'FIN' BLESP ;
'LISTE' RGAS1 ;
'LISTE' RGAS2 ;
*
* Computation of the conservative variables
*
RN1 = PN1 '/' (RGAS1 '*' TN1) ;
RN2 = PN2 '/' (RGAS2 '*' TN2) ;
GN1 = RN1 * VN1 ;
GN2 = RN2 * VN2 ;
LVEL = 'MOTS' 'UX' 'UY' 'UZ' ;
ECIN1 = 0.5D0 '*' ('PSCAL' VN1 VN1 LVEL LVEL) ;
ECIN2 = 0.5D0 '*' ('PSCAL' VN2 VN2 LVEL LVEL) ;
RETN1 = RN1 '*' (ETHER1 '+' ECIN1 '+' EFORM1) ;
RETN2 = RN2 '*' (ETHER2 '+' ECIN2 '+' EFORM2) ;
*
'RESPRO' RN1 RN2 GN1 GN2 RETN1 RETN2 ;
'FINPROC' ;
*
****************************************************************
****************************************************************
*********END EXACT SOLUTION AND PROCEDURES ********************
****************************************************************
****************************************************************
*
*
****************************************************************
****************************************************************
***** 1D mesh *****
****************************************************************
****************************************************************
*
*
* A4 ------------------- A5 ----------------------- A6
* | | |
* | ZONE1 | ZONE2 |
* | | |
* A1 ------------------- A2 ----------------------- A3
* | | |
* | L1 | L2 |
* >|<--------------------->|<----------------------->|
*
* AR = activation region
* ZONE1
* ZONE2
*
RAF = 200 ;
L1 = 1. ;
DX = L1 '/' RAF ;
L2 = DX ;
N1 = RAF ;
N2 = 1 ;
* For the sake of simplicity, we will only consider ny mesh in y-direction
NY = 3 ;
H1 = NY '*' DX ;
A1 = 0.0 0.0 ;
A2 = L1 0.0 ;
A3 = (L1 '+' L2) 0.0 ;
A4 = 0.0 H1 ;
A5 = L1 H1 ;
A6 = (L1 '+' L2) H1 ;
*
'OPTION' 'ELEM' 'SEG2' ;
A1A2 = A1 'DROIT' N1 A2 ;
A2A3 = A2 'DROIT' N2 A3 ;
A4A5 = A4 'DROIT' N1 A5 ;
A5A6 = A5 'DROIT' N2 A6 ;
*
A1A4 = A1 'DROIT' NY A4 ;
A2A5 = A2 'DROIT' NY A5 ;
A3A6 = A3 'DROIT' NY A6 ;
*
*
**** 2D mesh
*
* DOM1 is zone 1
* DOM3 is zone 2
*
*
'OPTION' 'ELEM' 'QUA4' ;
DOM1 = 'DALLER' A1A2 A2A5 ('INVERSE' A4A5) ('INVERSE' A1A4) ;
DOM2 = 'DALLER' A2A3 A3A6 ('INVERSE' A5A6) ('INVERSE' A2A5) ;
LIMG = A1A4 ;
LIMD = A3A6 ;
*
* DOMTOT = total region
*
DOMTOT = DOM1 'ET' DOM2 ;
'ELIMINATION' (1.0D-3 '*' DX) DOMTOT ;
DOMLIM = 'CONTOUR' DOMTOT ;
*
**** Line for graphics
*
LIMCON = 'MANUEL' 'QUA4'
(A1 'PLUS' ((0.25 * DX) (0.25 * DX)))
(A3 'PLUS' ((-0.25 * DX) (0.25 * DX)))
(A3 'PLUS' ((-0.25 * DX) (0.75 * DX)))
(A1 'PLUS' ((0.25 * DX) (0.75 * DX))) ;
LIMCON = LIMCON 'COULEUR' ROUG ;
*
ANGLE = 45.0 ;
DOMTOT2 = DOMTOT 'TOURNER' ANGLE (-10. -10.) ;
CHAMP = DOMTOT2 'MOIN' DOMTOT ;
FORM CHAMP ;
LIMCON = LIMCON 'TOURNER' ANGLE (-10. -10.) ;
*
**** The tables 'DOMAINE'
*
$DOMTOT = 'MODELISER' DOMTOT 'EULER' ;
$DOMLIM = 'MODELISER' DOMLIM 'EULER' ;
$DOM1 = 'MODELISER' DOM1 'EULER' ;
$DOM2 = 'MODELISER' DOM2 'EULER' ;
$LIMG = 'MODELISER' LIMG 'EULER' ;
$LIMD = 'MODELISER' LIMD 'EULER' ;
*
TDOMTOT = 'DOMA' $DOMTOT 'VF' ;
TDOMLIM = 'DOMA' $DOMLIM 'VF' ;
TDOM1 = 'DOMA' $DOM1 'VF' ;
TDOM2 = 'DOMA' $DOM2 'VF' ;
TLIMG = 'DOMA' $LIMG 'VF' ;
TLIMD = 'DOMA' $LIMD 'VF' ;
*
QDOMTOT = TDOMTOT . 'QUAF' ;
QDOMLIM = TDOMLIM . 'QUAF' ;
QDOM1 = TDOM1 . 'QUAF' ;
QDOM2 = TDOM2 . 'QUAF' ;
QLIMG = TLIMG . 'QUAF' ;
QLIMD = TLIMD . 'QUAF' ;
*
'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOMLIM ;
'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOM1 ;
'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOM2 ;
'ELIMINATION' QDOMTOT (0.001 '*' DX) QLIMG ;
'ELIMINATION' QDOMTOT (0.001 '*' DX) QLIMD ;
*
**** MOD1 = model (created just to display the CHAMELEMs)
*
MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
*
**** Line for graphics
*
*
LIGEVO = (TDOMTOT . 'CENTRE') 'INCL' LIMCON ;
P1 = 'POIN' (COOR 1 LIGEVO) 'MINIMUM' ;
'ORDONNER' LIGEVO P1 ;
LIGEVO = 'QUELCONQUE' 'SEG2' LIGEVO 'COULEUR' 'VERT';
'SI' GRAPH ;
'TRACER' (LIMCON 'ET' QDOMTOT 'ET' LIGEVO)
'TITR' ('CHAINE' 'Maillage: ' ('NBEL' DOMTOT) ' elts') ;
'FINSI' ;
* 'OPTION' 'SAUVER' ('CHAINE' 'mesh.sauv_raf' RAF) ;
* 'SAUVER' RAF DOMTOT $DOMTOT $DOM1 $DOM2 $DOM3 LIGEVO MOD1 ;
****************************************************************
****************************************************************
***** Initial conditions and gas properties. *****
****************************************************************
****************************************************************
*
*************************************************
**** The table for the properties of the gas ****
*************************************************
*
FCOMB = 0.0 ;
FDISS = 0.6 ;
PGAS = 'TABLE' ;
*
**** Order of the polynomial order for cv = cv(T)
* For T > TMAX, cv(T) = cv(Tmax)
*
PGAS . 'TMAX' = 6000.0 ;
PGAS . 'NORD' = 4 ;
*
**** Species involved in the mixture (before or after
* the chemical reaction)
*
PGAS . 'SPECIES' = 'MOTS' 'C2H2' 'O2 ' 'H2O ' 'CO2 ' 'CO ' 'O '
'N2 ';
*
*
**** Coefficient of the chemical reaction.
* Note that for the first species this coefficient should be positive
* Normal, we take it equal to 1.
*
* C2H2 '+' 2.5 O2 ---> H2O + (1-F) CO2 + F CO + F O
*
PGAS . 'CHEMCOEF' = 'PROG' 1.0 2.5 -1.0 (-2.0 '+' FDISS)
(-1.0 '*' FDISS) (-1.0 '*' FDISS) 0.0 ;
*
**** Coef with the gas properties
*
PGAS . 'C2H2' = 'TABLE' ;
PGAS . 'H2O ' = 'TABLE' ;
PGAS . 'CO2 ' = 'TABLE' ;
PGAS . 'O2 ' = 'TABLE' ;
PGAS . 'N2 ' = 'TABLE' ;
PGAS . 'CO ' = 'TABLE' ;
PGAS . 'O ' = 'TABLE' ;
*
**** Runiv (J/mole/K)
*
PGAS . 'RUNIV' = 8.31441 ;
*
**** W (kg/mole). Gas constant (J/kg/K = Runiv/W)
*
PGAS . 'C2H2' . 'W' = 26.03788D-3 ;
PGAS . 'CO2 ' . 'W' = (2. * 15.9994E-3) '+' 12.0107E-3 ;
PGAS . 'O2 ' . 'W' = 2. * 15.9994E-3 ;
PGAS . 'N2 ' . 'W' = 2. * 14.0067E-3 ;
PGAS . 'H2O ' . 'W' = (2. * 1.00797E-3) '+'
(15.9994E-3) ;
PGAS . 'CO ' . 'W' = 15.9994E-3 '+' 12.0107E-3 ;
PGAS . 'O ' . 'W' = 15.9994E-3 ;
*
**** Mass fraction of the first species before and after the combustion
* Final mass fractions of the species with positive coefficients.
* Final mass fractions of the species with non-positive coefficient.
* The mass fraction of the last species is not given.
*
*
NC2H2 = 1.;
NO2 = 2.5 * NC2H2 ;
NN2 = NO2 '*' 0.79 '/' 0.21 ;
NH2O = 1.0D-12 ;
NCO2 = 1.0D-12 ;
NCO = 1.0D-12 ;
NO = 1.0D-12 ;
MC2H2 = NC2H2 * (PGAS . 'C2H2' . 'W') ;
MO2 = NO2 * (PGAS . 'O2 ' . 'W') ;
MN2 = NN2 * (PGAS . 'N2 ' . 'W') ;
MH2O = NH2O * (PGAS . 'H2O ' . 'W') ;
MCO2 = NCO2 * (PGAS . 'CO2 ' . 'W') ;
MCO = NCO * (PGAS . 'CO ' . 'W') ;
MO = NO * (PGAS . 'O ' . 'W') ;
MTOT = MC2H2 '+' MO2 '+' MH2O '+' MN2 '+' MCO2 '+' MCO '+' MO;
YC2H2 = MC2H2 '/' MTOT ;
YO2 = MO2 '/' MTOT ;
YH2O = MH2O '/' MTOT ;
YN2 = MN2 '/' MTOT ;
YCO2 = MCO2 '/' MTOT ;
YCO = MCO '/' MTOT ;
YO = MO '/' MTOT ;
SI FAUX ;
*
**** Complete combustion
*
SI (NC2H2 > (0.4 * NO2)) ;
* C2H2 en excès
YO2F = 1.0D-12 ;
YC2H2F = YC2H2 '-' ((YO2 '-' YO2F) * ((PGAS . 'C2H2' . 'W') '/'
(2.5 * (PGAS . 'O2 ' . 'W')))) ;
'SINON' ;
* O2 en excès
YC2H2F = 1.0D-12 ;
YO2F = YO2 '-'
(2.5 * (YC2H2 '-' YC2H2F) * ((PGAS . 'O2 ' . 'W')
'/'(PGAS . 'C2H2' . 'W'))) ;
'FINSI' ;
'SINON' ;
YC2H2F = FCOMB * YC2H2 ;
YO2F = YO2 '-'
(2.5 * (YC2H2 '-' YC2H2F) * ((PGAS . 'O2 ' . 'W') '/'
(PGAS . 'C2H2' . 'W'))) ;
'FINSI' ;
YCO2F = YCO2 '+' ((YC2H2 '-' YC2H2F) '*' (2.0 '-' FDISS)
'*' ((PGAS . 'CO2 ' . 'W') '/'
(PGAS . 'C2H2' . 'W'))) ;
YCOF = YCO '+' ((YC2H2 '-' YC2H2F) '*' (FDISS)
'*' ((PGAS . 'CO ' . 'W') '/'
(PGAS . 'C2H2' . 'W'))) ;
YOF = YO '+' ((YC2H2 '-' YC2H2F) '*' (FDISS)
'*' ((PGAS . 'O ' . 'W') '/'
(PGAS . 'C2H2' . 'W'))) ;
YH2OF = YH2O '+' (YC2H2 '-' YC2H2F) '+' (YO2 '-' YO2F) '+'
(YCO2 '-' YCO2F) '+' (YO '-' YOF) '+' (YCO '-' YCOF) ;
YN2F = YN2 ;
ERRO = ((YC2H2F '+' YO2F '+' YH2OF '+' YN2F '+' YCO2F '+' YCOF
'+' YOF) '-' 1.0) ;
SI (ERRO > 1.0D-12) ;
'ERREUR' 21 ;
'FINSI' ;
*
* PGAS . 'CHEMCOEF' = 'PROG' 1.0 2.5 -1.0 (-2.0 '+' FDISS)
* (-1.0 '*' FDISS) (-1.0 '*' FDISS) ;
* PGAS . 'SPECIES' = 'MOTS' 'C2H2' 'O2 ' 'H2O ' 'CO2 ' 'CO ' 'O ';
PGAS . 'MASSFRA' = 'PROG' YC2H2 YC2H2F YO2F YH2O YCO2 YCO YO ;
*
**** Polynomial coefficients
*
PGAS . 'C2H2' . 'A' = 'PROG' 0.88323615E+03 0.20414991E+01
-0.76189061E-03 0.13137105E-06
-0.83517351E-11 ;
PGAS . 'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
-1.82753232E-08 2.44485692E-12 ;
PGAS . 'O2 ' . 'A' = 'PROG' 575.012333 0.350522002 -0.000128294865
2.33636971E-08 -1.53304905E-12;
PGAS . 'N2 ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
8.78233606E-09 -3.05514485E-13 ;
PGAS . 'CO2 ' . 'A' = 'PROG' 4.0835592237504562e+02
9.8512056453381769e-01 -4.4170720111493721e-04
8.5343398699551964e-08 -5.8566495097123519e-12 ;
PGAS . 'CO ' . 'A' = 'PROG' 0.66108887E+03 0.30437019E+00
-0.90998496E-04 0.12134406E-07 -0.58302668E-12 ;
PGAS . 'O ' . 'A' = 'PROG' 0.90791174E+03 -0.19422711E+00
0.94136383E-04 -0.17707293E-07 0.12172246E-11 ;
*
**** Formation enthalpies (energies) at 0K (J/Kg)
*
PGAS . 'C2H2' . 'H0K' = 9.05431D6 ;
PGAS . 'H2O ' . 'H0K' = -1.395D7 ;
PGAS . 'O2 ' . 'H0K' = -2.634D5 ;
PGAS . 'N2 ' . 'H0K' = -2.953D5 ;
PGAS . 'CO2 ' . 'H0K' = -93.965*((4.187E3)'/'(PGAS . 'CO2 ' . 'W'));
PGAS . 'CO ' . 'H0K' = -4.06313D6 ;
PGAS. 'O ' . 'H0K' = 1.54250E+07 ;
*************************************************
**** The initial conditions *********************
*************************************************
*
eps = 1.0D-64 ;
K0 = 3000. ;
*
T1 = 298. ;
alph11 = 1.0 '-' 1.0D-4 ;
alph12 = 1.0D-4 ;
alph21 = 1.0 '-' alph11 ;
alph22 = 1.0 - alph12 ;
T2 = 4000. ;
un1 = 0.0 ;
un2 = 0.0 ;
ut1 = 0.0 ;
ut2 = 0.0 ;
pre1 = 1.013D5 ;
pre2 = 1.013D5 ;
*
CHVN1 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 2 'UX' un1
'UY' ut1) ;
CHVN2 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 2 'UX' un2
'UY' ut2) ;
*
CHTN1 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' T1) ;
CHTN2 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' T2) ;
*
CHPN1 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' pre1) ;
CHPN2 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' pre2) ;
*
CHAL1 = ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL'
alph11) '+'
('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 1 'SCAL'
alph12) ;
CHAL2 = ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL'
alph21) '+'
('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 1 'SCAL'
alph22) ;
*
CHRN1 CHRN2 CHGN1 CHGN2 CHRET1 CHRET2 = PRIMCONS
PGAS CHTN1 CHTN2 CHPN1 CHPN2 CHVN1 CHVN2 ;
*
RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS CHAL1 CHAL2
(CHAL1 * CHRN1) (CHAL2 * CHRN2) (CHAL1 * CHGN1)
(CHAL2 * CHGN2) (CHAL1 * CHRET1) (CHAL2 * CHRET2)
CHTN1 CHTN2 EPS ;
TN1M = TN1 ;
TN2M = TN2 ;
*
**** Parameter for the time loop
*
* Names of the residuum components
*
LISTINCO = ('MOTS' 'ALF1' 'RN1' 'RNX1' 'RNY1' 'RET1'
'ALF2' 'RN2' 'RNX2' 'RNY2' 'RET2') ;
*
**** BC
*
KONLIM = 'DIFF' (TLIMG . 'CENTRE') (TLIMG . 'CENTRE') ;
CHPVID CHMVID = 'KOPS' 'MATRIK' ;
RESLIM = 'COPIER' CHPVID ;
*
****************************************************************
****************************************************************
***** Parameters for the computations ******
****************************************************************
****************************************************************
*
* Iterations
* Final time
* Safety Factor fot the time step
* Second order reconstruction?
NITER = 10000 ;
TFINAL = 0.35D-3 ;
SAFFAC = 0.75 ;
LOGSO = VRAI ;
* 'OPTION' 'SAUVER' 'parameters.sauv' ;
* 'SAUVER' METO NITER TFINAL SAFFAC LOGSO ;
****************************************************************
****************************************************************
***** The computation ******
****************************************************************
****************************************************************
AL1 = CHAL1 ;
AL2 = CHAL2 ;
ARN1 = CHAL1 * CHRN1 ;
ARN2 = CHAL2 * CHRN2 ;
AGN1 = CHAL1 * CHGN1 ;
AGN2 = CHAL2 * CHGN2 ;
AREN1 = CHAL1 * CHRET1 ;
AREN2 = CHAL2 * CHRET2 ;
*
AL10 = 'COPIER' AL1 ;
AL20 = 'COPIER' AL2 ;
ARN10 = 'COPIER' ARN1 ;
ARN20 = 'COPIER' ARN2 ;
AGN10 = 'COPIER' AGN1 ;
AGN20 = 'COPIER' AGN2 ;
AREN10 = 'COPIER' AREN1 ;
AREN20 = 'COPIER' AREN2 ;
TN10 = 'COPIER' TN1 ;
TN20 = 'COPIER' TN2 ;
*
**** Geometrical coefficient to compute grad(alpha)/|grad(alpha)|
*
EPSGR = 1.0D-6 ;
CHPL1 = CHPVID ;
CHPL2 = 'MANUEL' 'CHPO' (TDOMLIM . 'CENTRE') 2
'P1DX' 0.0 'P1DY' 0.0 ;
GRALP1 GRAL = 'PENT' $DOMTOT 'FACE'
'DIAMAN2' ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY')
AL1 CHPL1 CHPL2 ;
* 'SI' VRAI ;
* V1 = 'VECTEUR'
* ('NOMC' ('MOTS' 'P1DX' 'P1DY') GRALP1 ('MOTS' 'UX' 'UY')) ;
* 'TRACER' DOMTOT V1 DOMLIM ;
* 'FINSI' ;
GRALP1 = GRALP1 '+' EPSGR ;
GRALP1 = GRALP1 '/' (('PSCAL' GRALP1 GRALP1 ('MOTS' 'P1DX' 'P1DY')
('MOTS' 'P1DX' 'P1DY')) '**' 0.5) ;
* 'SI' VRAI ;
* V1 = 'VECTEUR'
* ('NOMC' ('MOTS' 'P1DX' 'P1DY') GRALP1 ('MOTS' 'UX' 'UY')) ;
* 'TRACER' DOMTOT V1 DOMLIM ;
* 'FINSI' ;
*
**** Geometrical coefficient to compute gradients
*
GRADR0 ALR0 COEFSCAL = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
('MOTS' 'SCAL') AL10 ;
GRADR0 = GRADR0 * 0.0 ;
ALR0 = ALR0 * 0.0 ;
GRADV0 ALV0 COEFVECT = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
('MOTS' 'UX' 'UY') AGN10 ;
GRADV0 = GRADV0 * 0.0 ;
ALV0 = ALV0 * 0.0 ;
VINF1 = 'MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' 10. ;
VINF2 = 'MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' 100. ;
*
**** Temporal loop
*
TPS = 0.0 ;
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
'MESSAGE' ('CHAINE' 'SAFFAC = ' SAFFAC) ;
'MESSAGE' ;
'TEMPS' 'ZERO' ;
'REPETER' BLITER NITER ;
*
**** Primitive variables
*
RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1 AL2
ARN1 ARN2 AGN1 AGN2 AREN1 AREN2 TN1M TN2M EPS ;
TN1M = 'COPIER' TN1 ;
TN2M = 'COPIER' TN2 ;
**** Gradient of alpha
GRALP1 = 'PENT' $DOMTOT 'FACE'
'DIAMAN2' ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY')
AL1 CHPL1 CHPL2 'GRADGEO' GRAL ;
* 'SI' VRAI ;
* V1 = 'VECTEUR'
* ('NOMC' ('MOTS' 'P1DX' 'P1DY') GRALP1 ('MOTS' 'UX' 'UY')) ;
* 'TRACER' DOMTOT V1 DOMLIM ;
* 'FINSI' ;
GRALP1 = GRALP1 '+' EPSGR ;
GRALP1 = GRALP1 '/' (('PSCAL' GRALP1 GRALP1 ('MOTS' 'P1DX' 'P1DY')
('MOTS' 'P1DX' 'P1DY')) '**' 0.5) ;
'SI' LOGSO ;
GRADA1 ALA1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') AL1 'GRADGEO' COEFSCAL ;
GRADR1 ALR1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') RN1 'GRADGEO' COEFSCAL ;
GRADP1 ALP1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PN1 'GRADGEO' COEFSCAL ;
GRADV1 ALV1 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY') VN1 'GRADGEO' COEFVECT ;
GRADA2 ALA2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') AL2 'GRADGEO' COEFSCAL ;
GRADR2 ALR2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') RN2 'GRADGEO' COEFSCAL ;
GRADP2 ALP2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PN2 'GRADGEO' COEFSCAL ;
GRADV2 ALV2 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY') VN2 'GRADGEO' COEFVECT ;
CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
'PRET' 'DEM' $DOMTOT
AL1 GRADA1 ALA1
AL2 GRADA2 ALA2
* AL1 GRADA1 ALR0
* AL2 GRADA2 ALR0
RN1 GRADR1 ALR1
RN2 GRADR2 ALR2
VN1 GRADV1 ALV1
VN2 GRADV2 ALV2
PN1 GRADP1 ALP1
PN2 GRADP2 ALP2 ;
'SINON' ;
CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
'PRET' 'DEM' $DOMTOT
AL1 GRADR0 ALR0
AL2 GRADR0 ALR0
RN1 GRADR0 ALR0
RN2 GRADR0 ALR0
VN1 GRADV0 ALV0
VN2 GRADV0 ALV0
PN1 GRADR0 ALR0
PN2 GRADR0 ALR0 ;
'FINSI' ;
SI ('EGA' METO 'AUSMPUP') ;
RESIDU DELTAT SURF = 'KONV' 'VF' 'DEM' 'RESI' METO 'CONS'
$DOMTOT PGAS LISTINCO AL1 AL2 CHFAL1 CHFAL2 CHFRN1 CHFRN2
CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS KONLIM VINF1 VINF2 ;
'SINON' ;
RESIDU DELTAT SURF = 'KONV' 'VF' 'DEM' 'RESI' METO 'CONS'
$DOMTOT PGAS LISTINCO AL1 AL2 CHFAL1 CHFAL2 CHFRN1 CHFRN2
CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS KONLIM ;
'FINSI' ;
RESIDU = RESIDU '+' RESLIM ;
*
* 'REPETER' BL1 ('DIME' LISTINCO) ;
* mot1 = 'EXTRAIRE' LISTINCO &BL1 ;
* valre1 = 'MAXIMUM' ('EXCO' RESLIM mot1) 'ABS' ;
* tit1 = 'CHAINE' 'Component ' mot1 ', reference value ' valre1 ;
* evaa = 'EVOL' 'CHPO' RESIDU mot1 LIGEVO ;
* 'DESSIN' evaa 'TITRE' tit1 ;
* 'FIN' BL1 ;
* 'OPTION' DONN 5 ;
* Problem with RNY1, RNY2
DT_CON = SAFFAC '*' DELTAT ;
*
**** The time step linked to tps
*
DTTPS = (TFINAL '-' TPS) * (1. '+' 1.0D-6) ;
*
**** Total time step
*
DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ;
*
**** Increment of the variables (predictor)
*
* RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1 AL2
* ARN1 ARN2 AGN1 AGN2 AREN1 AREN2 CHTN1 CHTN2 EPS ;
RESIDU = (0.5 * DTMIN) '*' RESIDU ;
DALP1 = 'EXCO' 'ALF1' RESIDU 'SCAL' ;
DRN1 = 'EXCO' 'RN1' RESIDU 'SCAL' ;
DGN1 = 'EXCO' ('MOTS' 'RNX1' 'RNY1') RESIDU ('MOTS' 'UX' 'UY') ;
DRET1 = 'EXCO' 'RET1' RESIDU 'SCAL' ;
DALP2 = 'EXCO' 'ALF2' RESIDU 'SCAL' ;
DRN2 = 'EXCO' 'RN2' RESIDU 'SCAL' ;
DGN2 = 'EXCO' ('MOTS' 'RNX2' 'RNY2') RESIDU ('MOTS' 'UX' 'UY') ;
DRET2 = 'EXCO' 'RET2' RESIDU 'SCAL' ;
AL1B = AL1 '+' DALP1 ;
ARN1B = ARN1 '+' DRN1 ;
AGN1B = AGN1 '+' DGN1 ;
AREN1B = AREN1 '+' DRET1 ;
AL2B = AL2 '+' DALP2 ;
ARN2B = ARN2 '+' DRN2 ;
AGN2B = AGN2 '+' DGN2 ;
AREN2B = AREN2 '+' DRET2 ;
*
**** Corrector
*
RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1B AL2B
ARN1B ARN2B AGN1B AGN2B AREN1B AREN2B TN1M TN2M EPS ;
**** Gradient of alpha
GRALP1 = 'PENT' $DOMTOT 'FACE'
'DIAMAN2' ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY')
AL1B CHPL1 CHPL2 'GRADGEO' GRAL ;
* 'SI' VRAI ;
* V1 = 'VECTEUR'
* ('NOMC' ('MOTS' 'P1DX' 'P1DY') GRALP1 ('MOTS' 'UX' 'UY')) ;
* 'TRACER' DOMTOT V1 DOMLIM ;
* 'FINSI' ;
GRALP1 = GRALP1 '+' EPSGR ;
GRALP1 = GRALP1 '/' (('PSCAL' GRALP1 GRALP1 ('MOTS' 'P1DX' 'P1DY')
('MOTS' 'P1DX' 'P1DY')) '**' 0.5) ;
'SI' LOGSO ;
GRADA1 ALA1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') AL1B 'GRADGEO' COEFSCAL ;
GRADR1 ALR1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') RN1 'GRADGEO' COEFSCAL ;
GRADP1 ALP1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PN1 'GRADGEO' COEFSCAL ;
GRADV1 ALV1 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY') VN1 'GRADGEO' COEFVECT ;
GRADA2 ALA2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') AL2B 'GRADGEO' COEFSCAL ;
GRADR2 ALR2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') RN2 'GRADGEO' COEFSCAL ;
GRADP2 ALP2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PN2 'GRADGEO' COEFSCAL ;
GRADV2 ALV2 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY') VN2 'GRADGEO' COEFVECT ;
CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
'PRET' 'DEM' $DOMTOT
AL1B GRADA1 ALA1
AL2B GRADA2 ALA2
* AL1 GRADA1 ALR0
* AL2 GRADA2 ALR0
RN1 GRADR1 ALR1
RN2 GRADR2 ALR2
VN1 GRADV1 ALV1
VN2 GRADV2 ALV2
PN1 GRADP1 ALP1
PN2 GRADP2 ALP2 ;
'SINON' ;
CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
'PRET' 'DEM' $DOMTOT
AL1B GRADR0 ALR0
AL2B GRADR0 ALR0
RN1 GRADR0 ALR0
RN2 GRADR0 ALR0
VN1 GRADV0 ALV0
VN2 GRADV0 ALV0
PN1 GRADR0 ALR0
PN2 GRADR0 ALR0 ;
'FINSI' ;
SI ('EGA' METO 'AUSMPUP') ;
RESIDU DELTAT SURF = 'KONV' 'VF' 'DEM' 'RESI' METO 'CONS'
$DOMTOT PGAS LISTINCO AL1B AL2B CHFAL1 CHFAL2 CHFRN1 CHFRN2
CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS KONLIM VINF1 VINF2 ;
'SINON' ;
RESIDU DELTAT SURF = 'KONV' 'VF' 'DEM' 'RESI' METO 'CONS'
$DOMTOT PGAS LISTINCO AL1B AL2B CHFAL1 CHFAL2 CHFRN1 CHFRN2
CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS KONLIM ;
'FINSI' ;
RESIDU = RESIDU '+' RESLIM ;
*
* 'REPETER' BL1 ('DIME' LISTINCO) ;
* mot1 = 'EXTRAIRE' LISTINCO &BL1 ;
* valre1 = 'MAXIMUM' ('EXCO' RESLIM mot1) 'ABS' ;
* tit1 = 'CHAINE' 'Component ' mot1 ', reference value ' valre1 ;
* evaa = 'EVOL' 'CHPO' RESIDU mot1 LIGEVO ;
* 'DESSIN' evaa 'TITRE' tit1 ;
* 'FIN' BL1 ;
* 'OPTION' DONN 5 ;
* Problem with RNY1, RNY2
RESIDU = DTMIN '*' RESIDU ;
DALP1 = 'EXCO' 'ALF1' RESIDU 'SCAL' ;
DRN1 = 'EXCO' 'RN1' RESIDU 'SCAL' ;
DGN1 = 'EXCO' ('MOTS' 'RNX1' 'RNY1') RESIDU ('MOTS' 'UX' 'UY') ;
DRET1 = 'EXCO' 'RET1' RESIDU 'SCAL' ;
DALP2 = 'EXCO' 'ALF2' RESIDU 'SCAL' ;
DRN2 = 'EXCO' 'RN2' RESIDU 'SCAL' ;
DGN2 = 'EXCO' ('MOTS' 'RNX2' 'RNY2') RESIDU ('MOTS' 'UX' 'UY') ;
DRET2 = 'EXCO' 'RET2' RESIDU 'SCAL' ;
TPS = TPS '+' DTMIN ;
AL1 = AL1 '+' DALP1 ;
ARN1 = ARN1 '+' DRN1 ;
AGN1 = AGN1 '+' DGN1 ;
AREN1 = AREN1 '+' DRET1 ;
AL2 = AL2 '+' DALP2 ;
ARN2 = ARN2 '+' DRN2 ;
AGN2 = AGN2 '+' DGN2 ;
AREN2 = AREN2 '+' DRET2 ;
'SI' (((&BLITER '/' 20) '*' 20) 'EGA' &BLITER) ;
'MESSAGE' ('CHAINE' 'ITER =' &BLITER ' TPS =' TPS) ;
'FINSI' ;
'SI' (TPS > TFINAL) ;
'QUITTER' BLITER ;
'FINSI' ;
'FIN' BLITER ;
* 'TEMPS' 'IMPR' ;
'TEMPS' ;
* 'OPTION' 'SAUVER'
* ('CHAINE' '/RES/alberto/Dgros/result_tube1D.sauv_' RAF
* 'tps_' TPS ) ;
* 'SAUVER' ;
* 'FIN' ;
****************************************************************
****************************************************************
***** The post-treatment ******
****************************************************************
****************************************************************
* 'OPTI' 'REST' 'result.sauv_1tps_5' ;
* 'REST' ;
****************************************************************
****************************************************************
***** The post-treatment ******
****************************************************************
****************************************************************
* 'OPTI' 'REST' 'result.sauv_1tps_5' ;
* 'REST' ;
*
**** The mesh
*
* 'SI' GRAPH ;
* 'TRACER' ('DOMA' $DOMTOT 'MAILLAGE')
* 'TITR' ('CHAINE' 'Maillage: ' ('NBEL' DOMTOT) ' elts');
* 'FINSI' ;
*
**** Initial conditions
*
MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
'SI' FAUX ;
RN10 RN20 VN10 VN20 PN10 PN20 TN10 TN20 = 'PRIM' 'DEM' PGAS
AL10 AL20 ARN10 ARN20 AGN10 AGN20 AREN10 AREN20 TN1M TN2M EPS ;
CHM_AL10 = 'KCHA' $DOMTOT 'CHAM' AL10 ;
CHM_RN10 = 'KCHA' $DOMTOT 'CHAM' RN10 ;
CHM_VN10 = 'KCHA' $DOMTOT 'CHAM' VN10 ;
CHM_PN10 = 'KCHA' $DOMTOT 'CHAM' PN10 ;
CHM_TN10 = 'KCHA' $DOMTOT 'CHAM' TN10 ;
CHM_AL20 = 'KCHA' $DOMTOT 'CHAM' AL20 ;
CHM_RN20 = 'KCHA' $DOMTOT 'CHAM' RN20 ;
CHM_VN20 = 'KCHA' $DOMTOT 'CHAM' VN20 ;
CHM_PN20 = 'KCHA' $DOMTOT 'CHAM' PN20 ;
CHM_TN20 = 'KCHA' $DOMTOT 'CHAM' TN20 ;
'TRAC' CHM_AL10 MOD1
'TITR' ('CHAINE' 'alp_1 at t=' 0.0) ;
'TRAC' CHM_RN10 MOD1
'TITR' ('CHAINE' 'rho_1 at t=' 0.0) ;
'TRAC' CHM_VN10 MOD1
'TITR' ('CHAINE' 'v_1 at t= ' 0.0) ;
'TRAC' CHM_PN10 MOD1
'TITR' ('CHAINE' 'p_1 at t= ' 0.0) ;
'TRAC' CHM_TN10 MOD1
'TITR' ('CHAINE' 'T_1 at t= ' 0.0) ;
'TRAC' CHM_AL20 MOD1
'TITR' ('CHAINE' 'alp_2 at t=' 0.0) ;
'TRAC' CHM_RN20 MOD1
'TITR' ('CHAINE' 'rho_2 at t=' 0.0) ;
'TRAC' CHM_VN20 MOD1
'TITR' ('CHAINE' 'v_2 at t= ' 0.0) ;
'TRAC' CHM_PN20 MOD1
'TITR' ('CHAINE' 'p_2 at t= ' 0.0) ;
'TRAC' CHM_TN20 MOD1
'TITR' ('CHAINE' 'T_2 at t= ' 0.0) ;
'FINSI' ;
*
**** Results
*
RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1 AL2
ARN1 ARN2 AGN1 AGN2 AREN1 AREN2 TN1M TN2M EPS ;
'SI' FAUX ;
CHM_AL1 = 'KCHA' $DOMTOT 'CHAM' AL1 ;
CHM_RN1 = 'KCHA' $DOMTOT 'CHAM' RN1 ;
CHM_VN1 = 'KCHA' $DOMTOT 'CHAM' VN1 ;
CHM_PN1 = 'KCHA' $DOMTOT 'CHAM' PN1 ;
CHM_TN1 = 'KCHA' $DOMTOT 'CHAM' TN1 ;
CHM_AL2 = 'KCHA' $DOMTOT 'CHAM' AL2 ;
CHM_RN2 = 'KCHA' $DOMTOT 'CHAM' RN2 ;
CHM_VN2 = 'KCHA' $DOMTOT 'CHAM' VN2 ;
CHM_PN2 = 'KCHA' $DOMTOT 'CHAM' PN2 ;
CHM_TN2 = 'KCHA' $DOMTOT 'CHAM' TN2 ;
'TRAC' CHM_AL1 MOD1
'TITR' ('CHAINE' 'alp_1 at t=' TPS) ;
'TRAC' CHM_RN1 MOD1
'TITR' ('CHAINE' 'rho_1 at t=' TPS) ;
'TRAC' CHM_VN1 MOD1
'TITR' ('CHAINE' 'v_1 at t= ' TPS) ;
'TRAC' CHM_PN1 MOD1
'TITR' ('CHAINE' 'p_1 at t= ' TPS) ;
'TRAC' CHM_TN1 MOD1
'TITR' ('CHAINE' 'T_1 at t= ' TPS) ;
'TRAC' CHM_AL2 MOD1
'TITR' ('CHAINE' 'alp_2 at t=' TPS) ;
'TRAC' CHM_RN2 MOD1
'TITR' ('CHAINE' 'rho_2 at t=' TPS) ;
'TRAC' CHM_VN2 MOD1
'TITR' ('CHAINE' 'v_2 at t= ' TPS) ;
'TRAC' CHM_PN2 MOD1
'TITR' ('CHAINE' 'p_2 at t= ' TPS) ;
'TRAC' CHM_TN2 MOD1
'TITR' ('CHAINE' 'T_2 at t= ' TPS) ;
'FINSI' ;
*
*** Evolution Objects. We put the speed in a frame
* parallel to the mesh.
*
VN1old = 'COPIER' VN1 ;
VN2old = 'COPIER' VN2 ;
UX1N = (('EXCO' 'UX' VN1) * ('COS' ANGLE)) '+'
(('EXCO' 'UY' VN1) * ('SIN' ANGLE)) ;
UY1N = (('EXCO' 'UX' VN1) * ('SIN' ANGLE)) '-'
(('EXCO' 'UY' VN1) * ('COS' ANGLE)) ;
UX2N = (('EXCO' 'UX' VN2) * ('COS' ANGLE)) '+'
(('EXCO' 'UY' VN2) * ('SIN' ANGLE)) ;
UY2N = (('EXCO' 'UX' VN2) * ('SIN' ANGLE)) '-'
(('EXCO' 'UY' VN2) * ('COS' ANGLE)) ;
VN1 = ('NOMC' UX1N 'UX') '+' ('NOMC' UY1N 'UY') ;
VN2 = ('NOMC' UX2N 'UX') '+' ('NOMC' UY2N 'UY') ;
LXMED = 'PROG' ('DIME' LEX_XOT) '*' 1.0 ;
TABLEG = 'TABLE' ;
TABLEG . 1 = 'MOT' 'TIRC MARQ PLUS' ;
TABLEG . 2 = 'MOT' 'TIRC MARQ CROI' ;
TABLEG . 3 = 'MOT' '' ;
TABLEG . 'TITRE' = 'TABLE' ;
TABLEG . 'TITRE' . 1 = 'MOT' 'phase 1' ;
TABLEG . 'TITRE' . 2 = 'MOT' 'phase 2' ;
TABLEG . 'TITRE' . 3 = 'MOT' 'alp_1*ph_1 + alp_2 ph_2' ;
TABLE2 = 'TABLE' ;
TABLE2 . 1 = 'MOT' 'TIRC MARQ PLUS' ;
TABLE2 . 2 = 'MOT' '' ;
TABLE2 . 'TITRE' = 'TABLE' ;
TABLE2 . 'TITRE' . 1 = 'MOT' 'Num' ;
TABLE2 . 'TITRE' . 2 = 'MOT' 'Exa' ;
eval1 = ('EVOL' 'CHPO' al1 LIGEVO) 'COULEUR' 'ROUG' ;
eval2 = ('EVOL' 'CHPO' al2 LIGEVO) 'COULEUR' 'BLEU' ;
evone = (eval1 '+' eval2) 'COULEUR' 'VERT' ;
xmin = 'MINIMUM' ('EXTRAIRE' eval1 'ABSC') ;
xmax = 'MAXIMUM' ('EXTRAIRE' eval1 'ABSC') ;
evrn1 = ('EVOL' 'CHPO' rn1 LIGEVO) 'COULEUR' 'ROUG' ;
evrn2 = ('EVOL' 'CHPO' rn2 LIGEVO) 'COULEUR' 'BLEU' ;
evrn = ((eval1 * evrn1) '+' (eval2 * evrn2)) 'COULEUR' 'VERT' ;
evrna = 'EVOL' 'MANU' 'x' (LXMED '+' (LEX_XOT * TFINAL * -1))
'SCAL' LEX_RHO ;
evpn1 = ('EVOL' 'CHPO' pn1 LIGEVO) 'COULEUR' 'ROUG' ;
evpn2 = ('EVOL' 'CHPO' pn2 LIGEVO) 'COULEUR' 'BLEU' ;
evpn = ((eval1 * evpn1) '+' (eval2 * evpn2)) 'COULEUR' 'VERT' ;
evpna = 'EVOL' 'MANU' 'x' (LXMED '+' (LEX_XOT * TFINAL * -1))
'SCAL' LEX_P ;
evux1 = ('EVOL' 'CHPO' vn1 LIGEVO 'UX') 'COULEUR' 'ROUG' ;
evux2 = ('EVOL' 'CHPO' vn2 LIGEVO 'UX') 'COULEUR' 'BLEU' ;
evux = ((eval1 * evux1) '+' (eval2 * evux2)) 'COULEUR' 'VERT' ;
evuna = 'EVOL' 'MANU' 'x' (LXMED '+' (LEX_XOT * TFINAL * -1))
'SCAL' (-1 * LEX_U) ;
evuy1 = ('EVOL' 'CHPO' vn1 LIGEVO 'UY') 'COULEUR' 'ROUG' ;
evuy2 = ('EVOL' 'CHPO' vn2 LIGEVO 'UY') 'COULEUR' 'BLEU' ;
evuy = ((eval1 * evuy1) '+' (eval2 * evuy2)) 'COULEUR' 'VERT' ;
evtn1 = ('EVOL' 'CHPO' tn1 LIGEVO) 'COULEUR' 'ROUG' ;
evtn2 = ('EVOL' 'CHPO' tn2 LIGEVO) 'COULEUR' 'BLEU' ;
evtn = ((eval1 * evtn1) '+' (eval2 * evtn2)) 'COULEUR' 'VERT' ;
evtna = 'EVOL' 'MANU' 'x' (LXMED '+' (LEX_XOT * TFINAL * -1))
'SCAL' LEX_T ;
RET = AREN1 '+' AREN2 ;
evret = ('EVOL' 'CHPO' RET LIGEVO) 'COULEUR' 'VERT' ;
evreta = 'EVOL' 'MANU' 'x' (LXMED '+' (LEX_XOT * TFINAL * -1))
'SCAL' LEX_RET ;
'SI' GRAPH ;
'DESSIN' (eval1 'ET' eval2 'ET' evone) 'TITRE'
'alpha (alp_1 + alp_2)'
'LEGE' TABLEG 'XBOR' XMIN XMAX ;
'DESSIN' ((evone '-' (evone '/' evone))) 'TITRE'
'(alp_1 + alp_2), error' ;
'DESSIN' (evrn1 'ET' evrn2 'ET' evrn) 'TITRE' 'Density'
'LEGE' TABLEG 'XBOR' XMIN XMAX ;
'DESSIN' (evrn et evrna) 'TITRE' 'Density' 'LEGE' TABLE2
'XBOR' XMIN XMAX ;
'DESSIN' (evpn1 'ET' evpn2 'ET' evpn) 'TITRE' 'Pressure'
'LEGE' TABLEG 'XBOR' XMIN XMAX ;
'DESSIN' (evpn et evpna) 'TITRE' 'Pressure' 'LEGE' TABLE2
'XBOR' XMIN XMAX ;
'DESSIN' (evtn1 'ET' evtn2 'ET' evtn) 'TITRE' 'Temperature'
'LEGE' TABLEG 'XBOR' XMIN XMAX ;
'DESSIN' (evtn et evtna) 'TITRE' 'Temperature' 'LEGE' TABLE2
'XBOR' XMIN XMAX ;
'DESSIN' (evux1 'ET' evux2 'ET' evux) 'TITRE' 'ux'
'LEGE' TABLEG 'XBOR' XMIN XMAX ;
'DESSIN' (evux et evuna) 'TITRE' 'ux' 'LEGE' TABLE2
'XBOR' XMIN XMAX ;
'DESSIN' (evuy1 'ET' evuy2 'ET' evuy) 'TITRE' 'uy'
'LEGE' TABLEG 'XBOR' XMIN XMAX ;
'DESSIN' (evuy) 'TITRE' 'uy' 'XBOR' XMIN XMAX ;
'DESSIN' (evret et evreta) 'TITRE' 'ret' 'LEGE' TABLE2
'XBOR' XMIN XMAX ;
'FINSI' ;
*
**** Test
*
lxa = 'EXTRAIRE' EVRNA 'ABSC' ;
lrhoa = 'EXTRAIRE' EVRNA 'ORDO' ;
lx = 'EXTRAIRE' EVRN 'ABSC' ;
lrho = 'EXTRAIRE' EVRN 'ORDO' ;
lrhoa1 = 'IPOL' lx lxa lrhoa ;
erro = 'ABS' ((lrho '-' lrhoa1) '/' lrhoa1) ;
everro = 'EVOL' 'MANU' 'x' lx 'y' erro ;
erroint = ('PRIM' everro) ;
'SI' GRAPH ;
'DESSIN' everro 'TITRE' 'Erreur' ;
'DESSIN' erroint 'TITRE' 'Integral de l Erreur' ;
'FINSI' ;
errL1 = ('MAXIMUM' ('EXTRAIRE' erroint 'ORDO')) '/'
('MAXIMUM' lx) ;
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'ERR_L1 = ' errL1) ;
'MESSAGE' ;
SI (errL1 > 5.0D-2) ;
'ERREUR' 'errL1 too big !!! ' ;
'ERREUR' 5 ;
'FINSI' ;
*
* Conservation of the total energy
*
CHVO = TDOMTOT . 'XXVOLUM' ;
RETN0 = AREN10 '+' AREN20 ;
TOTE0 = 'XTY' RETN0 CHVO ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
RETN = AREN1 '+' AREN2 ;
TOTE = 'XTY' RETN CHVO ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
ERRO = 'ABS' (TOTE '-' TOTE0) '/' TOTE ;
SI (ERRO > 1.0D-6) ;
'MESSAGE' ('CHAINE' 'Energy error ' erro) ;
'MESSAGE' 'Conservation of the total energy' ;
'ERREUR' 21 ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer