* fichier : crebe12.dgibi
*--------------------------------------------
* GIBIANE file for the HDR test E12.3.2 *
* in 2D *
* We make 100 iterations on the coarse mesh *
* and compare the overpressure in the room *
* R1904 with earlier computed values *
* CREBCOM combustion model is used *
*--------------------------------------------
*--------------------------------------------
* Here we construct the mesh for the *
* HDR test E12.3.2 *
*--------------------------------------------
'OPTION' 'ECHO' 1 'DIME' 2 'ELEM' 'QUA4'
'TRAC' 'X' ;
*-------------------------------------------
* Given a line and a mesh, this procedure
* computes a new line which containes the
* points of the mesh closest to the points
* of the given line
*-------------------------------------------
'DEBPROC' LIGPROC ;
******* inputs
'ARGUMENT' LIGG*'MAILLAGE' MAI1*'MAILLAGE' ;
NPOI = 'NBNO' LIGG ;
MAI2 = 'CHANGER' MAI1 'POI1' ;
* First element of the line
P1 = LIGG 'POIN' 1 ;
P1P = MAI2 'POIN' 'PROC' P1 ;
P2 = LIGG 'POIN' 2 ;
P2P = MAI2 'POIN' 'PROC' P2 ;
'REPETER' BL1 (NPOI '-' 1) ;
'SI' ('NEG' P1P P2P) ;
NPOI = (NPOI '-' &BL1) '+' 1 ;
'QUITTER' BL1 ;
'FINSI' ;
P2 = LIGG 'POIN' (&BL1 '+' 2) ;
P2P = LIGG 'POIN' 'PROC' P2 ;
'FIN' BL1 ;
LIGRES = 'MANUEL' 'SEG2' P1P P2P ;
* The other element of the line
'REPETER' BL1 (NPOI '-' 2) ;
P1P = P2P ;
P2 = LIGG 'POIN' (&BL1 '+' 2) ;
P2P = MAI2 'POIN' 'PROC' P2 ;
'SI' ('NEG' P1P P2P) ;
LIGRES = LIGRES 'ET' ('MANUEL' 'SEG2' P1P P2P) ;
'FINSI' ;
'FIN' BL1 ;
'RESPRO' LIGRES ;
'FINPROC' ;
*--------------------------------------------------------
*---refinement of the mesh (10 or 16)
RAF = 10 ;
*-- Length of the thirst and second room
L1 = 7.32D0 ;
*-- Height of the first window
H1 = 0.5D0 ;
*-- Height of the first and second room
H2 = 4.0D0 ;
*---------------------------------
* Number of cells in each region
*---------------------------------
N1 = RAF ;
*---- number of points in window1
** It is BETTER to have NY1 even!
NY1 = 'ENTIER' (RAF '/' 4) ;
*---------------------------------
DX = L1 '/' N1 ;
*---- cell size in window
DX1 = H1 '/' NY1 ;
*---- number of points in the outer window
NY2 = 'ENTIER' ((H2 '-' H1) '/' DX) ;
***---- wall thickness
DW = H1 ;
*-----------------------------------------
* First and Second Room (R1904 and R1905)
*-----------------------------------------
LB = L1 '-' (7.5 '*' H1) ;
A1 = 0.0 0.0 ;
A12 = LB 0.0 ;
A13 = LB (H1 '/' 2) ;
A14 = 0.0 (H1 '/' 2) ;
*------------------------------------
* "Lower-left" domain
*-------------------------------------
NY1M = NY1 '/' 2 ;
A1A12 = A1 'DROIT' A12 'DINI' DX1 'DFIN' (DW '/' NY1) ;
A12A13 = A12 'DROIT' NY1M A13 ;
A13A14 = A13 'DROIT' A14 'DINI' (DW '/' NY1) 'DFIN' DX1 ;
A14A1 = A14 'DROIT' NY1M A1 ;
DOM1L = 'DALLER' A1A12 A12A13 A13A14 A14A1 'PLAN' ;
*--------------------------------------
* "Lower-right" domain
*--------------------------------------
LB1 = L1 '-' (2.5 '*' H1) ;
AI1 = LB1 0.0 ;
A2 = L1 0.0 ;
AI3 = L1 (H1 '/' 2) ;
AI4 = LB1 (H1 '/' 2) ;
*------------------------------------
DENF = (DW '/' NY1) ;
AI1A2 = AI1 'DROIT' A2 'DINI' DENF 'DFIN' DENF ;
A2AI3 = A2 'DROIT' NY1M AI3 ;
AI3AI4 = AI3 'DROIT' AI4 'DINI' DENF 'DFIN' DENF ;
AI4AI1 = AI4 'DROIT' NY1M AI1 ;
DOM1R = 'DALLER' AI1A2 A2AI3 AI3AI4 AI4AI1 'PLAN' ;
DOM1 = DOM1L 'ET' DOM1R ;
'ELIMINATION' DOM1 1.0D-6 ;
*----------------------------------------------------
* Line for postreatment
*----------------------------------------------------
*** NNN is the number of points in the disc
NNN = NY1 '*' 4 ;
xcel = LB '+' (H1 '/' 2) ;
ycel = H1 '/' 2 ;
ACEL = xcel ycel ;
LIGB = A13 'DROIT' NY1M ACEL ;
*------------------
xmm = LB '+' (2.5 '*' H1) ;
'REPETER' BL1 NNN ;
ACEL0 = ACEL ;
xcel = xcel '+' DX1 ;
dop = 18.0625 '*' H1 '*' H1 ;
dop = dop '-' ((xcel-xmm) '*' (xcel-xmm)) ;
dop = dop '**' 0.5 ;
ycel = dop '-' (3.25 '*' H1) ;
ACEL = xcel ycel ;
LIGB = LIGB 'ET' (ACEL0 'DROIT' 1 ACEL) ;
'FIN' BL1 ;
LIGB = LIGB 'ET' (ACEL 'DROIT' AI4 'DINI' DX1 'DFIN' DX1 ) ;
*--------------------------------------------
LIGBG = ('INVERSE' A13A14) 'ET' LIGB 'ET'
('INVERSE' AI3AI4) ;
*-------------------------------------
* "HIGH" domain
*-------------------------------------
A3 = L1 H1 ;
A4 = 0.0 H1 ;
A5 = L1 H2 ;
A6 = 0.0 H2 ;
AI5 = LB H2 ;
*---------------------------------
AI3A5 = AI3 'DROIT' A5 'DINI' DX1 'DFIN' DX1 ;
A5AI5 = A5 'DROIT' AI5 'DINI' DX1 'DFIN' DX1 ;
AI5A6 = AI5 'DROIT' A6 'DINI' DX1 'DFIN' DX1 ;
LIGH = A5AI5 'ET' AI5A6 ;
A6A14 = A6 'DROIT' A14 'DINI' DX1 'DFIN' DX1 ;
*---------------------------
DOM2 = 'DALLER' LIGBG AI3A5 LIGH A6A14 'PLAN' ;
*---------------------------------------------
* The Channel (junction between rooms)
*---------------------------------------------
A21 = A2 'PLUS' (DW 0.0) ;
A31 = A3 'PLUS' (DW 0.0) ;
A2A21 = A2 'DROIT' NY1 A21 ;
A21A31 = A21 'DROIT' NY1 A31 ;
A31A3 = A31 'DROIT' NY1 A3 ;
A3A2 = A3 'DROIT' NY1 A2 ;
*-------------------------
DOMC = 'DALLER' A2A21 A21A31 A31A3 A3A2 'PLAN' ;
*--------------------------------------------
* Third room (R1801)
*--------------------------------------------
*---- Length of the third room
L3 = 5.0D0 ;
*--- length before the window
LB = 1.8D0 ;
*--- width of window
LW = 0.668D0 ;
*--- length after the window
LA = L3 '-' LB '-' LW ;
*----------------------------------------
* Low Domain
*----------------------------------------
H3 = (9.5 '-' H2) '-' H1 ;
MH3 = (-1.0D0) '*' H3 ;
A70 = (L1 '+' DW) MH3 ;
A71 = (L1 '+' LB '+' DW) MH3 ;
A72 = (L1 '+' LB '+' DW) 0.0 ;
*----------------------------
DW1 = LW '/' NY1 ;
A70A71 = A70 'DROIT' A71 'DINI' (DW '/' NY1) 'DFIN' DW1 ;
A71A72 = A71 'DROIT' A72 'DINI' DX1 'DFIN' DX1 ;
A72A21 = A72 'DROIT' A21 'DINI' DW1 'DFIN' (DW '/' NY1) ;
A21A70 = A21 'DROIT' A70 'DINI' DX1 'DFIN' DX1 ;
*------------------------------
DOM3A = 'DALLER' A70A71 A71A72 A72A21 A21A70 'PLAN' ;
*------------------------------
A7W = A71 'PLUS' (LW 0.0) ;
A7U = A72 'PLUS' (LW 0.0) ;
*--------------------------
A71A7W = A71 'DROIT' NY1 A7W ;
A7WA7U = A7W 'DROIT' A7U 'DINI' DX1 'DFIN' DX1 ;
A7UA72 = A7U 'DROIT' NY1 A72 ;
A72A71 = 'INVERSE' (A71A72) ;
*-------------------------
DOM3B = 'DALLER' A71A7W A7WA7U A7UA72 A72A71 'PLAN' ;
*-------------------------
A81 = A7W 'PLUS' (LA 0.0) ;
A82 = A7U 'PLUS' (LA 0.0) ;
*-------------------------
A7WA81 = A7W 'DROIT' A81 'DINI' DW1 'DFIN' DX1 ;
A81A82 = A81 'DROIT' A82 'DINI' DX1 'DFIN' DX1 ;
A82A7U = A82 'DROIT' A7U 'DINI' DX1 'DFIN' DW1 ;
A7UA7W = 'INVERSE' (A7WA7U) ;
*-------------------------
DOM3C = 'DALLER' A7WA81 A81A82 A82A7U A7UA7W 'PLAN' ;
*-------------------------
DOM3 = DOM3A 'ET' DOM3B 'ET' DOM3C ;
'ELIMINATION' DOM3 1.0D-6 ;
*-----------------------------------------
* Middle domain
*-----------------------------------------
A73 = A72 'PLUS' (0.0 H1) ;
***--- first layer
A21A72 = A21 'DROIT' A72 'DINI' (DW '/' NY1) 'DFIN' DW1 ;
A72A73 = A72 'DROIT' NY1 A73 ;
A73A31 = A73 'DROIT' A31 'DINI' DW1 'DFIN' (DW '/' NY1) ;
A31A21 = A31 'DROIT' NY1 A21 ;
*-----------------------------
DOM4A = 'DALLER' A21A72 A72A73 A73A31 A31A21 'PLAN' ;
*-------------------------------------------
A7M = A7U 'PLUS' (0.0 H1) ;
*---------------
A72A7U = A72 'DROIT' NY1 A7U ;
A7UA7M = A7U 'DROIT' NY1 A7M ;
A7MA73 = A7M 'DROIT' NY1 A73 ;
A73A72 = 'INVERSE' (A72A73) ;
*---------------------
DOM4B = 'DALLER' A72A7U A7UA7M A7MA73 A73A72 'PLAN' ;
*---------------------------
A83 = A82 'PLUS' (0.0 H1) ;
*---------------------------
A7UA82 = A7U 'DROIT' A82 'DINI' DW1 'DFIN' DX1 ;
A82A83 = A82 'DROIT' NY1 A83 ;
A83A7M = A83 'DROIT' A7M 'DINI' DX1 'DFIN' DW1 ;
A7MA7U = 'INVERSE' (A7UA7M) ;
*---------------------------
DOM4C = 'DALLER' A7UA82 A82A83 A83A7M A7MA7U 'PLAN' ;
*---------------------------
DOM4 = DOM4A 'ET' DOM4B 'ET' DOM4C ;
'ELIMINATION' DOM4 1.0D-6 ;
*-----------------------------------------
* High domain
*-----------------------------------------
A32 = A31 'PLUS' (0.0 H2) ;
A74 = A73 'PLUS' (0.0 H2) ;
*---------------------------
A31A73 = 'INVERSE' (A73A31) ;
A73A74 = A73 'DROIT' A74 'DINI' DX1 'DFIN' DW1 ;
A74A32 = A74 'DROIT' A32 'DINI' DW1 'DFIN' (DW '/' NY1) ;
A32A31 = A32 'DROIT' A31 'DINI' DW1 'DFIN' DX1 ;
*---------------------------
DOM5A = 'DALLER' A31A73 A73A74 A74A32 A32A31 'PLAN' ;
*---------------------------
A7H = A7M 'PLUS' (0.0 H2) ;
*---------------------------
A73A7M = 'INVERSE' (A7MA73) ;
A7MA7H = A7M 'DROIT' A7H 'DINI' DX1 'DFIN' DW1 ;
A7HA74 = A7H 'DROIT' NY1 A74 ;
A74A73 = 'INVERSE' (A73A74) ;
*---------------------------
DOM5B = 'DALLER' A73A7M A7MA7H A7HA74 A74A73 'PLAN' ;
*---------------------------
A84 = A83 'PLUS' (0.0 H2) ;
*---------------------------
A7MA83 = A7M 'DROIT' A83 'DINI' DW1 'DFIN' DX1 ;
A83A84 = A83 'DROIT' A84 'DINI' DX1 'DFIN' DW1 ;
A84A7H = A84 'DROIT' A7H 'DINI' DX1 'DFIN' DW1 ;
A7HA7M = 'INVERSE' (A7MA7H) ;
*-----------------------------
DOM5C = 'DALLER' A7MA83 A83A84 A84A7H A7HA7M 'PLAN' ;
*-------------------------------------
DOM5 = DOM5A 'ET' DOM5B 'ET' DOM5C ;
'ELIMINATION' DOM5 1.0D-6 ;
****----------*****----------------------******-------
* Upper vent
*---**********-----**********************------*******
A75 = A74 'PLUS' (0.0 LW) ;
A7F = A7H 'PLUS' (0.0 LW) ;
*---------------------------
A74A7H = A74 'DROIT' NY1 A7H ;
A7HA7F = A7H 'DROIT' NY1 A7F ;
A7FA75 = A7F 'DROIT' NY1 A75 ;
A75A74 = A75 'DROIT' NY1 A74 ;
*---------------------------
DOMW = 'DALLER' A74A7H A7HA7F A7FA75 A75A74 'PLAN' ;
FRONTH = A7FA75 ;
*----------------------------------*
* Ignition Zone *
*----------------------------------*
** height at which we ignite
HIGN = H2 '/' 2.0D0 ;
*-------------------------
PP1 = 0.0 HIGN ;
PP2 = PP1 'PLUS' (0.0 (DX1)) ;
PP3 = PP2 'PLUS' (0.0 (DX1)) ;
gg = (-1.0D0) '*' DX '*' 0.6 ;
PP5 = PP1 'PLUS' (0.0 gg) ;
PP51 = PP5 'PLUS' (0.0 gg) ;
PO1 = A6A14 'POIN' 'PROCHE' PP1 ;
PO2 = A6A14 'POIN' 'PROCHE' PP2 ;
PO31 = A6A14 'POIN' 'PROCHE' PP3 ;
PO5 = A6A14 'POIN' 'PROCHE' PP5 ;
PO51 = A6A14 'POIN' 'PROCHE' PP51 ;
**--------------------------------
PPP = DX1 0.0 ;
PA1 = A1A12 'POIN' 'PROCHE' PPP ;
PPP = (2.0 '*' DX) 0.0 ;
PA2 = A1A12 'POIN' 'PROCHE' PPP ;
XXX1 = 'COORDONNEE' 1 PA1 ;
XXX2 = 'COORDONNEE' 1 PA2 ;
YYY2 = 'COORDONNEE' 2 PO2 ;
YYY3 = 'COORDONNEE' 2 PO31 ;
PO4 = XXX1 YYY2 ;
PO7 = XXX1 YYY3 ;
*-----------------------------
* Just eight cells
*-----------------------------
PO1PO2 = PO1 'DROIT' 1 PO2 ;
PO2PO31 = PO2 'DROIT' 1 PO31 ;
PO4PO7 = PO4 'DROIT' 1 PO7 ;
PO5PO51 = PO5 'DROIT' 1 PO51 ;
**------------------------------
'SI' (RAF 'EGA' 10) ;
NZD = 1 ;
KEF = -1.0D0 ;
'SINON' ;
KEF = -2.0D0 ;
NZD = 2 ;
'FINSI' ;
DIG1 = 'TRANSLATION' PO1PO2 NZD ((KEF '*' XXX1) 0.0) ;
DIG5 = 'TRANSLATION' PO2PO31 NZD ((KEF '*' XXX1) 0.0) ;
DIG6 = 'TRANSLATION' PO4PO7 1 ((XXX2-XXX1) 0.0) ;
DIG7 = 'TRANSLATION' PO5PO51 1 (XXX1 0.0) ;
'SI' (RAF '>' 10) ;
YO31 = 'COORDONNEE' 2 PO31 ;
YO2 = 'COORDONNEE' 2 PO1 ;
DIFY = YO31 '-' YO2 ;
VECUP = (0.0 DIFY) ;
DIG6 = DIG1 'PLUS' VECUP ;
DIG7 = DIG5 'PLUS' VECUP ;
'FINSI' ;
*-----------------------------
'SI' (RAF '>' 10) ;
DOMI = DIG1 'ET' DIG5 'ET' DIG6 'ET' DIG7 ;
'SINON' ;
DOMI = DIG1 'ET' DIG5 ;
'FINSI' ;
'ELIMINATION' DOMI 1.0D-6;
* 'TRACER' DOMI ;
*----------------------------------------------------
* Total Domain
*----------------------------------------------------
DOMTOT = DOM1 'ET' DOM2 'ET' DOM3 'ET' DOM4 'ET' DOM5 'ET'
DOMC 'ET' DOMW 'ET' DOMI ;
'ELIMINATION' DOMTOT 1D-6 ;
*-------------------------------------
* Total Domain without ignition region
*-------------------------------------
DREST = 'DIFF' DOMTOT DOMI ;
*----------------------------------
* The fird room (R1801)
*----------------------------------
DTHIRD = DOM3 'ET' DOM4 'ET' DOM5 'ET' DOMC 'ET' DOMW ;
'ELIMINATION' DTHIRD 1.0D-6 ;
*----------------------------------------
* First and second room (R1904 and R1905)
*----------------------------------------
DSEC = DOM1 'ET' DOM2 ;
*-------------------------------
* Lignes de postreatment
*-------------------------------
LIGA = LIGBG 'PLUS' (0.0 (DX1 '/' 2.0)) ;
A3A31 = 'INVERSE' A31A3 ;
A31A73 = 'INVERSE' A73A31 ;
A73A7M = 'INVERSE' A7MA73 ;
A7MA83 = 'INVERSE' A83A7M ;
LIGT = A3A31 'ET' A31A73 'ET' A73A7M 'ET' A7MA83 ;
'ELIMINATION' LIGT 1.0D-6 ;
dis1 = (DX1 '/' 2) '-' (DW '/' 2.0) ;
LIGB2 = LIGT 'PLUS' (0.0 dis1) ;
LIGC = LIGA 'ET' LIGB2 ;
'ELIMINATION' LIGC 1.0D-6;
LIG1 = LIGC ;
* 'TRACER' (DOMTOT 'ET' LIG1) ;
*-------------------------------
* Points for capturing variables
*-------------------------------
*--------- inside the first two rooms
dish = (H2 '/' 2) '-' (DX1 '/' 2) ;
disl = (H1 '/' 2) '+' (DX1 '/' 2) ;
PCEL1 = (L1 '/' 6) dish ;
PCEL2 = (L1 '/' 6) disl ;
PCEL3 = (L1 '/' 2) dish ;
PCEL4 = (L1 '/' 2) disl ;
PCEL5 = (L1 '-' (L1 '/' 6)) dish ;
PCEL6 = (L1 '-' (L1 '/' 6)) disl ;
*-------- inside the third room
dis1 = L1 '+' (2.0 '*' DW) ;
PCEL7 = dis1 (2.0 '*' DW) ;
dis2 = (LB '+' (LW '/' 2) '+' (DX1 '/' 2.0)) '+' L1 '+' DW ;
PCEL8 = dis2 (DW '-' 5.0) ;
PCEL9 = dis2 (4.5 '-' DW) ;
* 'TRACER' (DOMTOT 'ET' PCEL1 'ET' PCEL2 'ET' PCEL3 'ET' PCEL4
* 'ET' PCEL5 'ET' PCEL6 'ET' PCEL7 'ET' PCEL8 'ET' PCEL9) ;
*
* 'TRACER' (('CONTOUR' DREST)
* 'ET' PCEL1 'ET' PCEL2 'ET' PCEL3 'ET' PCEL4
* 'ET' PCEL5 'ET' PCEL6 'ET' PCEL7 'ET' PCEL8 'ET' PCEL9) ;
*----------------------------------
* Creating the models
*----------------------------------
MDNS = 'EULER' ;
*---------------------------------
MDOMI = 'MODELISER' DOMI MDNS ;
MDREST = 'MODELISER' DREST MDNS ;
MFRONTH = 'MODELISER' FRONTH MDNS ;
MDOMTOT = 'MODELISER' DOMTOT MDNS ;
MDTHIRD = 'MODELISER' DTHIRD MDNS ;
MDSEC = 'MODELISER' DSEC MDNS ;
*----------------------------------
$DOMI = 'DOMA' MDOMI 'VF' ;
$DREST = 'DOMA' MDREST 'VF' ;
$FRONTH = 'DOMA' MFRONTH 'VF' ;
$DOMTOT = 'DOMA' MDOMTOT 'VF' ;
$DTHIRD = 'DOMA' MDTHIRD 'VF' ;
$DSEC = 'DOMA' MDSEC 'VF' ;
*----------------------------------
QDOMI = $DOMI . 'QUAF' ;
QDREST = $DREST . 'QUAF' ;
QFRONTH = $FRONTH . 'QUAF' ;
QDOMTOT = $DOMTOT . 'QUAF' ;
QDTHIRD = $DTHIRD . 'QUAF' ;
QDSEC = $DSEC . 'QUAF' ;
*----------------------------------
'ELIMINATION' QDOMTOT (1.0D-6 '/' RAF) QDOMI ;
'ELIMINATION' QDOMTOT (1.0D-6 '/' RAF) QDREST ;
'ELIMINATION' QDOMTOT (1.0D-6 '/' RAF) QFRONTH ;
'ELIMINATION' QDOMTOT (1.0D-6 '/' RAF) QDTHIRD ;
'ELIMINATION' QDOMTOT (1.0D-6 '/' RAF) QDSEC ;
*-------------------------------
LIN1 = LIGPROC LIG1 ('DOMA' MDOMTOT 'CENTRE') ;
* 'TRACER' (DOMTOT 'ET' LIN1) ;
*----------------------------------
* Ligne at the "coridor"
*----------------------------------
VECL = ((-0.5 '*' DX1) (0.5 '*' DX1)) ;
PL2 = A21 'PLUS' VECL ;
PL2 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PL2 ;
VECL = ((-0.5 '*' DX1) (-0.5 '*' DX1)) ;
PL3 = A31 'PLUS' VECL ;
PL3 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PL3 ;
LIG2 = PL2 'DROIT' (NY1 '-' 1) PL3 ;
'SI' (RAF 'EGA' 10) ;
LIN2 = LIG2 ;
'SINON' ;
LIN2 = LIGPROC LIG2 ('DOMA' MDOMTOT 'CENTRE') ;
'FINSI' ;
* 'TRACER' (DOMTOT 'ET' LIG2) ;
* 'TRACER' (DOMTOT 'ET' LIN2) ;
*--------------------------------
MDLIN2 = 'MODELISER' LIN2 MDNS ;
$DLIN2 = 'DOMA' MDLIN2 'VF' ;
QDLIN2 = $DLIN2 . 'QUAF' ;
'ELIMINATION' QDOMTOT (1.0D-6 '/' RAF) QDLIN2 ;
*----------------------------------
* "Centralising the capteurs"
*----------------------------------
PCE1 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL1 ;
PCE2 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL2 ;
PCE3 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL3 ;
PCE4 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL4 ;
PCE5 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL5 ;
PCE6 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL6 ;
PCE7 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL7 ;
PCE8 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL8 ;
PCE9 = ('DOMA' MDOMTOT 'CENTRE') 'POIN' 'PROCHE' PCEL9 ;
* 'TRACER' (DOMTOT 'ET' PCE1 'ET' PCE2 'ET' PCE3 'ET' PCE4
* 'ET' PCE5 'ET' PCE6 'ET' PCE7 'ET' PCE8 'ET' PCE9) ;
* 'TRACER' (DOMTOT 'ET' ('DOMA' MDOMI 'CENTRE')) ;
numbel = 'NBEL' DOMTOT ;
*-------------------------------------------
* Cells sizes: min, max and averaged
*-------------------------------------------
VOLC = 'DOMA' MDOMTOT 'VOLUME' ;
SIZE = VOLC '**' 0.5 ;
*------- minimal
SMIN = 'MINIMUM' SIZE ;
*------- maximal
SMAX = 'MAXIMUM' SIZE ;
*------- averaged
CHUN = 'MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE')
1 'SCAL' 1.0 'NATU' 'DISCRET' ;
SIZEAV = 'XTY' SIZE CHUN ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
SIZEAV = SIZEAV '/' numbel ;
*-------------------------------
* Saving the results
*-------------------------------
*'OPTION' 'SAUVER' ('CHAINE' '/serg/HDR12.3.2_dir/'
* 'mail' RAF '.sauv') ;
* 'SAUV' ;
*-----------------------------------------------
* Initial conditions '+' gas properties
*-----------------------------------------------
* 'OPTION' 'ECHO' 1 'DIME' 2 'ELEM' 'QUA4'
* 'TRAC' 'X' ;
*
* RAF = 10 ;
*
* 'OPTION' 'REST' ('CHAINE' '/serg/HDR12.3.2_dir/mail' RAF '.sauv') ;
* 'RESTITUER' ;
zero = 1.0D-9 ;
GRAPH = FAUX ;
***************************************************************
***** PROCEDURE POUR CALCULER LES 'CP' ET 'CV' *****
***************************************************************
'DEBPROC' CPCV TN*'FLOTTANT' PGAZ*'TABLE' ;
*
CPP = 'TABLE' ;
CVV = 'TABLE' ;
TT = TN ;
T2 = TT '*' TT ;
T3 = T2 '*' TT ;
T4 = T3 '*' TT ;
'REPETER' BL1 ('DIME' (PGAZ . 'ESPEULE')) ;
NOMCOM = 'EXTRAIRE' &BL1 (PGAZ . 'ESPEULE');
A0 = ('EXTRAIRE' 1 (PGAZ . NOMCOM . 'A')) ;
A1 = ('EXTRAIRE' 2 (PGAZ . NOMCOM . 'A')) '/' 2.0 ;
A2 = ('EXTRAIRE' 3 (PGAZ . NOMCOM . 'A')) '/' 3.0 ;
A3 = ('EXTRAIRE' 4 (PGAZ . NOMCOM . 'A')) '/' 4.0 ;
A4 = ('EXTRAIRE' 5 (PGAZ . NOMCOM . 'A')) '/' 5.0 ;
CVV . NOMCOM = A0 '+' (A1*TT) '+' (A2 '*' T2) '+'
(A3 '*' T3) '+' (A4 '*' T4) ;
CPP . NOMCOM = (CVV . NOMCOM) '+' (PGAZ . NOMCOM . 'R') ;
'FIN' BL1 ;
NOMCOM = PGAZ . 'ESPNEULE' ;
A0 = ('EXTRAIRE' 1 (PGAZ . NOMCOM . 'A')) ;
A1 = ('EXTRAIRE' 2 (PGAZ . NOMCOM . 'A')) '/' 2.0 ;
A2 = ('EXTRAIRE' 3 (PGAZ . NOMCOM . 'A')) '/' 3.0 ;
A3 = ('EXTRAIRE' 4 (PGAZ . NOMCOM . 'A')) '/' 4.0 ;
A4 = ('EXTRAIRE' 5 (PGAZ . NOMCOM . 'A')) '/' 5.0 ;
CVV . NOMCOM = A0 '+' (A1*TT) '+' (A2 '*' T2) '+'
(A3 '*' T3) '+' (A4 '*' T4) ;
CPP . NOMCOM = (CVV . NOMCOM) '+' (PGAZ . NOMCOM . 'R') ;
'FINPROC' CPP CVV ;
***************************************************************
***** PROCEDURE POUR CALCULER LES VARIABLES CONSERVATIVES *****
***************************************************************
'DEBPROC' CONS RN*'CHPOINT' VN*'CHPOINT' PN*'CHPOINT'
YN*'CHPOINT' $DM*'MMODEL' PGAZ*'TABLE' ;
*
A0 = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
A1 = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
A2 = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
A3 = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
A4 = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
RTOT = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
YTOT = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
'REPETER' BL1 ('DIME' (PGAZ . 'ESPEULE')) ;
NOMCOM = 'EXTRAIRE' &BL1 (PGAZ . 'ESPEULE');
YCEL = 'EXCO' NOMCOM YN ;
RYCEL = 'NOMC' (YCEL * RN) NOMCOM 'NATU' 'DISCRET' ;
YTOT = YTOT '+' YCEL ;
RTOT = RTOT '+' ((PGAZ . NOMCOM . 'R') * YCEL) ;
A0 = A0 '+' (('EXTRAIRE' 1 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A1 = A1 '+' (('EXTRAIRE' 2 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A2 = A2 '+' (('EXTRAIRE' 3 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A3 = A3 '+' (('EXTRAIRE' 4 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A4 = A4 '+' (('EXTRAIRE' 5 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
'SI' ('EGA' &BL1 1) ;
RYN = 'KCHT' $DM 'SCAL' 'CENTRE' 'COMP' NOMCOM RYCEL ;
'SINON' ;
RYN = RYN 'ET'
('KCHT' $DM 'SCAL' 'CENTRE' 'COMP' NOMCOM RYCEL) ;
'FINSI' ;
'FIN' BL1 ;
YCEL = ('KCHT' $DM 'SCAL' 'CENTRE' 1.0) '-' YTOT ;
NOMCOM = PGAZ . 'ESPNEULE' ;
A0 = A0 '+' (('EXTRAIRE' 1 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A1 = A1 '+' (('EXTRAIRE' 2 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A2 = A2 '+' (('EXTRAIRE' 3 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A3 = A3 '+' (('EXTRAIRE' 4 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A4 = A4 '+' (('EXTRAIRE' 5 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
RTOT = RTOT '+' ((PGAZ . NOMCOM . 'R') '*' YCEL);
TN = PN '/' (RN '*' RTOT) ;
T2 = TN '*' TN ;
T3 = T2 '*' TN ;
T4 = T3 '*' TN ;
T5 = T4 '*' TN ;
ETHER = (A0 * TN) '+' ((A1 '/' 2.0) * T2) '+'
((A2 '/' 3.0) * T3) '+' ((A3 '/' 4.0) * T4) '+'
((A4 '/' 5.0) * T5) ;
GN = RN '*' VN ('MOTS' 'SCAL' 'SCAL') ('MOTS' 'UX ' 'UY ')
('MOTS' 'UX ' 'UY ') ;
GN = 'KCHT' $DM 'VECT' 'CENTRE' GN ;
ECIN = 0.5D0 '*' ('PSCAL' GN VN ('MOTS' 'UX ' 'UY ')
('MOTS' 'UX ' 'UY ')) ;
REN = (RN '*' ETHER) '+' ECIN ;
REN = 'KCHT' $DM 'SCAL' 'CENTRE' REN ;
'FINPROC' GN REN RYN ;
************************************************************
************************************************************
* *
***** Gas properties *****
* *
************************************************************
************************************************************
PGAZ = 'TABLE' ;
PGAZB = 'TABLE' ;
* Polynomial degree of specific heats
PGAZ . 'NORD' = 4 ;
PGAZB . 'NORD' = 0 ;
* Species explicitly treated in the Euler Equations
PGAZ . 'ESPEULE' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' ;
PGAZB . 'ESPEULE' = PGAZ . 'ESPEULE' ;
* Species non explicitly treated
PGAZ . 'ESPNEULE' = 'N2 ';
PGAZB . 'ESPNEULE' = PGAZ . 'ESPNEULE' ;
* Gas species properties
PGAZ . 'H2 ' = 'TABLE' ;
PGAZ . 'H2O ' = 'TABLE' ;
PGAZ . 'N2 ' = 'TABLE' ;
PGAZ . 'O2 ' = 'TABLE' ;
*----------------
PGAZB . 'H2 ' = 'TABLE' ;
PGAZB . 'H2O ' = 'TABLE' ;
PGAZB . 'N2 ' = 'TABLE' ;
PGAZB . 'O2 ' = 'TABLE' ;
* R (J/Kg/K) and molar masses
mH2 = 2. '*' 1.00797E-3 ;
mO2 = 2. '*' 15.9994E-3 ;
mH2O = mh2 '+' (0.5 '*' mo2) ;
mN2 = 2 '*' 14.0067E-3 ;
RGAS = 8.31441 ;
*-------------------------------------
PGAZ . 'H2 ' . 'R' = RGAS '/' mh2 ;
PGAZ . 'H2O ' . 'R' = RGAS '/' mh2o ;
PGAZ . 'N2 ' . 'R' = RGAS '/' mn2 ;
PGAZ . 'O2 ' . 'R' = RGAS '/' mo2 ;
*-------------------------------------
PGAZB . 'H2 ' . 'R' = RGAS '/' mh2 ;
PGAZB . 'H2O ' . 'R' = RGAS '/' mh2o ;
PGAZB . 'N2 ' . 'R' = RGAS '/' mn2 ;
PGAZB . 'O2 ' . 'R' = RGAS '/' mo2 ;
*-------------------------------------
* Polynomials regressions coefficients
PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 0.000862203836
-2.37281455E-07 1.84701105E-11 ;
PGAZ . 'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
-1.82753232E-08 2.44485692E-12 ;
PGAZ . 'N2 ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
8.78233606E-09 -3.05514485E-13 ;
PGAZ . 'O2 ' . 'A' = 'PROG' 575.012333 0.350522002 -0.000128294865
2.33636971E-08 -1.53304905E-12;
*-------------------------------------------------------
* Formation enthalpies at 0K (J/Kg))
*
* h_i(0K) = h_i(T0) '-' \int_0^{T0} cp_i(x) dx
* = h_i(T0) '-' (\int_0^{T0} cv_i(x) dx '+' R_i * T0)
*-------------------------------------------------------
PGAZ . 'H2 ' . 'H0K' = -4.195D6 ;
PGAZ . 'H2O ' . 'H0K' = -1.395D7 ;
PGAZ . 'N2 ' . 'H0K' = -2.953D5 ;
PGAZ . 'O2 ' . 'H0K' = -2.634D5 ;
*-------------------------------------
PGAZB . 'H2 ' . 'H0K' = -4.195D6 ;
PGAZB . 'H2O ' . 'H0K' = -1.395D7 ;
PGAZB . 'N2 ' . 'H0K' = -2.953D5 ;
PGAZB . 'O2 ' . 'H0K' = -2.634D5 ;
*-----------------------------------
*** Names of the passive scalars
*-----------------------------------
PGAZ . 'SCALPASS' = 'MOTS' 'H2IN' 'H2FI' 'K0 ' ;
*-----------------------------------------------
PGAZB . 'SCALPASS' = 'MOTS' 'H2IN' 'H2FI' 'K0 ' ;
****************************************************
****************************************************
* Initial conditions *
****************************************************
****************************************************
*** inside the domain
tg = 337.0D0 ;
pg = 100000.0D0 ;
xh2o = 0.26D0 ;
xh2 = 0.1D0 ;
xo2 = 0.1344D0 ;
xn2 = 1.0D0 - (xh2 + xo2 + xh2O) ;
uxg = 0.0 ;
uyg = 0.0 ;
*-------------------------
td = 337.0D0 ;
uxd = 0.0D0 ;
uyd = 0.0D0 ;
pd = 100000.0D0 ;
*---- values of K_0
K0IN = 0.41D0 ;
K0IG = 0.41D0 ;
K03 = 3.7D0 ;
*--------------------------------------------
* After burning: AIBC (same pressure)
*--------------------------------------------
tc = 1086.38864D0 ;
xh2c = zero ;
xo2c = 0.0888421D0 ;
xh2oc = 0.3789469D0 ;
xn2c = 1.0D0 '-' (xh2c '+' xo2c '+' xh2oc) ;
pc = pg ;
*---------------------------------------------
CP1 CV1 = CPCV tg PGAZ ;
*'-'*'-'*'-'*'-'*'-'*-
PGAZB . 'CP' = 'TABLE' ;
PGAZB . 'CP' . 'H2 ' = CP1 . 'H2' ;
PGAZB . 'CP' . 'H2O ' = CP1 . 'H2O' ;
PGAZB . 'CP' . 'O2 ' = CP1 . 'O2' ;
PGAZB . 'CP' . 'N2 ' = CP1 . 'N2' ;
*----------------------------------------
PGAZB . 'CV' = 'TABLE' ;
PGAZB . 'CV' . 'H2 ' = CV1 . 'H2' ;
PGAZB . 'CV' . 'H2O ' = CV1 . 'H2O' ;
PGAZB . 'CV' . 'O2 ' = CV1 . 'O2' ;
PGAZB . 'CV' . 'N2 ' = CV1 . 'N2' ;
*---------------------------------------------
* Creating CHAMPOINTS
*---------------------------------------------
*** Pressure
MAILC = 'DOMA' MDOMTOT 'CENTRE' ;
MAILCI = 'DOMA' MDOMI 'CENTRE' ;
MAILCR = 'DOMA' MDREST 'CENTRE' ;
MAILC2 = 'DOMA' MDSEC 'CENTRE' ;
MAILC3 = 'DOMA' MDTHIRD 'CENTRE' ;
*********************************
PINI = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' pc 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCR 1 'SCAL' pg 'NATU' 'DISCRET') ;
*** Temperature
TINI = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' tc
'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCR 1 'SCAL' tg
'NATU' 'DISCRET') ;
**** Velocity
WINI = 'MANUEL' 'CHPO' MAILC 2 'UX' uxg 'UY' uyg ;
*-------------------------------------------------------
* Combustion is complete, i.e. \ksi = 1 after combustion
*-------------------------------------------------------
CSIMAX = 'MANUEL' 'CHPO' MAILC 1 'SCAL' 1.0D0 ;
**** Molar fractions before combustion
XH21 = ('MANUEL' 'CHPO' MAILCI 1 'H2' xh2c
'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCR 1 'H2' xh2
'NATU' 'DISCRET') ;
XH2O1 = ('MANUEL' 'CHPO' MAILCI 1 'H2O' xh2oc
'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCR 1 'H2O' xh2o
'NATU' 'DISCRET') ;
XO21 = ('MANUEL' 'CHPO' MAILCI 1 'O2' xo2c
'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCR 1 'O2' xo2
'NATU' 'DISCRET') ;
XN21 = ('MANUEL' 'CHPO' MAILCI 1 'N2' xn2c
'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCR 1 'N2' xn2
'NATU' 'DISCRET') ;
XN = XH21 'ET' XH2O1 'ET' XO21 'ET' XN21 ;
**** Mass Fractions
MMOL = ('MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 1 'H2' mH2
'NATURE' 'DISCRET') 'ET'
('MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 1 'O2' mO2
'NATURE' 'DISCRET') 'ET'
('MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 1 'H2O' mH2O
'NATURE' 'DISCRET') 'ET'
('MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 1 'N2' mN2
'NATURE' 'DISCRET') ;
MTOT = 'PSCAL' MMOL XN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
YH21 = (XH21 '*' mH2) '/' MTOT ;
YO21 = (XO21 '*' mO2) '/' MTOT ;
YN21 = (XN21 '*' mN2) '/' MTOT ;
YH2O1 = (XH2O1 '*' mH2O) '/' MTOT ;
aaa = 'MAXIMUM' (YH21) ;
*----------------------------------------------
* Burned gas
* First we find how much of the burned gas
* produces using 1 mole of unburned gas
*----------------------------------------------
CHUN = 'MANUEL' 'CHPO' MAILC 1 'SCAL' 1.0D0 'NATU' 'DISCRET' ;
H2MAG = (0.5 '*' ('EXCO' 'H2' XN)) 'MASQUE' 'SUPERIEUR'
('EXCO' 'O2' XN) ;
O2MAG = CHUN '-' H2MAG ;
XO22a = 'MANUEL' 'CHPO' MAILC 1 'O2' zero ;
XH22a = XH21 '-' ('NOMC' 'H2' (2.0 '*' (XO21 '-' XO22a))) ;
XO22a = XO22a ;
XH22a = XH22a ;
XH22b = 'MANUEL' 'CHPO' MAILC 1 'H2' zero ;
XO22b = XO21 '-' ('NOMC' 'O2' (0.5 '*' (XH21 '-' XH22b))) ;
XH22 = (XH22a '*' H2MAG) '+' (XH22b '*' O2MAG) ;
XO22 = (XO22a '*' H2MAG) '+' (XO22b '*' O2MAG) ;
XH2O2 = XH2O1 '+' ('NOMC' 'H2O' (2 '*' (XO21 '-' XO22))) ;
XN22 = XN21 ;
*---------------------------------------------------
* Verifying that mtot2 = mtot
*---------------------------------------------------
MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
('MOTS' 'H2' 'O2' 'H2O' 'N2')
('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
'SI' (('MAXIMUM' (mtot2 '-' mtot) 'ABS') '>' zero) ;
'MESSAGE' 'problem1' ;
'FINSI' ;
***** xtot2 = mole of burned gas per mole of unburned gas
xtot2 = 'PSCAL' (XH22 '+' XO22 '+' XH2O2 '+' XN22)
CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL') ;
***** X now is the real molar fractions
XH22 = XH22 '/' xtot2 ;
XO22 = XO22 '/' xtot2 ;
XH2O2 = XH2O2 '/' xtot2 ;
XN22 = XN22 '/' xtot2 ;
xtot3 = 'PSCAL' (XH22 '+' XO22 '+' XH2O2 '+' XN22)
CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL') ;
'SI' (('MAXIMUM' (xtot3 '-' 1) 'ABS') '>' zero) ;
'MESSAGE' 'problem2' ;
'FINSI' ;
**** mtot2 = weight of 1 mole of burned gas
MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
('MOTS' 'H2' 'O2' 'H2O' 'N2')
('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
**** Mass fractions of burned gas after complete combustion
YH22 = (XH22 '*' MH2) '/' MTOT2 ;
YO22 = (XO22 '*' MO2) '/' MTOT2 ;
YN22 = (XN22 '*' MN2) '/' MTOT2 ;
YH2O2 = (XH2O2 '*' MH2O) '/' MTOT2 ;
*---------------------------------------------
* Conservative variables
*---------------------------------------------
RINI = (('NOMC' 'SCAL' yh21) '*' (PGAZ . 'H2' . 'R')) '+'
(('NOMC' 'SCAL' yo21) '*' (PGAZ . 'O2' . 'R')) '+'
(('NOMC' 'SCAL' yh2o1) '*' (PGAZ . 'H2O' . 'R')) '+'
(('NOMC' 'SCAL' yn21) '*' (PGAZ . 'N2' . 'R')) ;
RN = PINI '/' (RINI '*' TINI) ;
YN = YH21 '+' YO21 '+' YH2O1 ;
GN RETN RYN = CONS RN WINI PINI YN MDOMTOT PGAZ ;
*--------------------------------------
* Passive scalars
*--------------------------------------
YININ = 'MANUEL' 'CHPO' MAILC 1 'H2' aaa 'NATU' 'DISCRET' ;
YFINN = 'MANUEL' 'CHPO' MAILC 1 'H2' zero 'NATU' 'DISCRET' ;
K0N = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' K0IG
'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILC2 1 'SCAL' K0IN
'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILC3 1 'SCAL' K03
'NATU' 'DISCRET') ;
*---------------------------------------
RSN = RN '*' (('NOMC' 'H2IN' YININ 'NATU' 'DISCRET') 'ET'
('NOMC' 'H2FI' YFINN 'NATU' 'DISCRET') 'ET'
('NOMC' 'K0' K0N 'NATU' 'DISCRET')) ;
*-------------------
***** Verification
*-------------------
VN P T Y S GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
*-----------------------
PCEL = 'MAXIMUM' PINI ;
TCEL = 'MAXIMUM' TINI ;
'SI'(('MAXIMUM' ((P '-' PINI) '/' PCEL) 'ABS') '>' (zero '*' 10.0)) ;
'MESSAGE' 'problem3' ;
'LISTE' ('MAXIMUM' ((P '-' PINI) '/' PCEL) 'ABS') ;
'FINSI' ;
'SI'(('MAXIMUM' ((T '-' TINI) '/' TCEL) 'ABS') '>' (zero '*' 10.0)) ;
'MESSAGE' 'problem4' ;
'LISTE' ('MAXIMUM' ((T '-' TINI) '/' TCEL) 'ABS') ;
'FINSI' ;
'SI' (('MAXIMUM' VN 'ABS') '>' zero) ;
'MESSAGE' 'problem5' ;
'LISTE' ('MAXIMUM' VN 'ABS') ;
'FINSI' ;
'SI'(('MAXIMUM' (Y '-' YN) 'ABS') '>' zero) ;
'MESSAGE' 'problem6' ;
'LISTE' ('MAXIMUM' (Y '-' YN) 'ABS') ;
'FINSI' ;
*--------------------------------------
* Graphics
*--------------------------------------
'SI' GRAPH ;
CHM_TN = 'KCHA' MDOMTOT 'CHAM' TINI ;
CHM_RN = 'KCHA' MDOMTOT 'CHAM' RN ;
CHM_PN = 'KCHA' MDOMTOT 'CHAM' PINI ;
CHM_VN = 'KCHA' MDOMTOT 'CHAM' WINI ;
CHM_H2 = 'KCHA' MDOMTOT 'CHAM' YN ;
CHM_YIN = 'KCHA' MDOMTOT 'CHAM' YININ ;
CHM_YFI = 'KCHA' MDOMTOT 'CHAM' YFINN ;
CHM_K0 = 'KCHA' MDOMTOT 'CHAM' K0N ;
'TRACER' CHM_TN MDOMTOT
'TITR' ('CHAINE' 'Temperature at t=' 0.0) ;
'TRACER' CHM_RN MDOMTOT
'TITR' ('CHAINE' 'RN at t=' 0.0) ;
'TRACER' CHM_PN MDOMTOT
'TITR' ('CHAINE' 'PN at t=' 0.0) ;
'TRACER' CHM_VN MDOMTOT
'TITR' ('CHAINE' 'VN at t=' 0.0) ;
'TRACER' CHM_H2 MDOMTOT
'TITR' ('CHAINE' 'H2 at t=' 0.0) ;
'TRACER' CHM_YIN MDOMTOT
'TITR' ('CHAINE' 'YIN at t=' 0.0) ;
'TRACER' CHM_YFI MDOMTOT
'TITR' ('CHAINE' 'YFI at t=' 0.0) ;
'TRACER' CHM_K0 MDOMTOT
'TITR' ('CHAINE' 'K0 at t=' 0.0) ;
'FINSI' ;
*-----------------------------------------------
* Parameters for the computation
*-----------------------------------------------
*** Upwind scheme
METO = 'SS' ;
**** Iterations
NITER = 100 ;
**** Final time (no restriction)
TFINAL = 10000.1 ;
**** Safety factor for the time step
SAFFAC = 0.7 ;
*** Second order ?
LOGSO = VRAI ;
LOGSO = FAUX ;
*** EPS_CC '-' epsilon in CREBCOM criterion
EPS_CC = 0.5 ;
**** "explosive" initial conditions?
LOGEXP = VRAI ;
LOGEXP = FAUX ;
*** Are there any Boundary conditions
LOGBC = VRAI ;
*-------------------------------
* Time intervals for the output
*-------------------------------
LISTT = 'PROG' 4.5 ;
*** Mesh size
DELTAXN = ('DOMA' MDOMTOT 'VOLUME') '**' 0.5 ;
DELX = 'MAXIMUM' DELTAXN ;
*-------------------------------------
TABLIM = 'TABLE' ;
pwin = (('MAXIMUM' RN) '*' 9.8) '*' 4.668 ;
TABLIM . 'CHPOUT' =
'MANUEL' 'CHPO' ('DOMA' MFRONTH 'CENTRE') 1 'PN' (pg '-' pwin) ;
TABLIM . 'H' = 50.0D0 ;
TABLIM . 'TNUL' =
'MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE') 1 'SCAL' tg ;
TABLIM . 'FREQ' = 1 ;
*----- gravity ----------
TABLIM . 'GRAV' = 'MANUEL' 'CHPO' ('DOMA' MDOMTOT 'CENTRE')
'UX' 0.0D0 'UY' (-9.8D0) ;
*----------------------------------------------
* Table for evolution of the variables with time
* at the fixed points
*----------------------------------------------
TABC = 'TABLE' 'CAPTEUR' ;
TABC . 'LPOINTS' = PCEL1 'ET' PCEL2 'ET' PCEL3 'ET' PCEL4
'ET' PCEL5 'ET' PCEL6 'ET' PCEL7 'ET' PCEL8 'ET' PCEL9 ;
TABC . 'LTPS' = 'PROG';
TABC . 'NOMC' = 'MOTS' 'C1' 'C2' 'C3' 'C4' 'C5' 'C6'
'C7' 'C8' 'C9' ;
TABC . 'COMP' = 'MOTS' 'PN' 'TN' ;
TABC . 'FREQ' = 10 ;
MSLIN2 = 'DOMA' MDLIN2 'SOMMET' ;
bnb = 'NBEL' MSLIN2 ;
VOLGOR = DX1 '*' DX1 '*' bnb ;
*-----------------------------------------------
* In MDOM1 we 'burn' the gas;
* In order to do so we shall apply the
* operator 'FLAM' in MDOM1 with eps=0
* and DELTAT >> DELTATC
*----------------------------------------------
YN = RYN '/' RN ;
SN = RSN '/' RN ;
K0N = 'EXCO' 'K0' SN ;
YININ = 'EXCO' 'H2IN' SN 'H2' ;
YFINN = 'EXCO' 'H2FI' SN 'H2' ;
DELTATC = 0.25 '*' ('MINIMUM' (DELTAXN '/' K0N)) ;
LMOT1 = 'MOTS' 'H2' 'O2' 'H2O' ;
**** Stocheom. coefficients
LCOEF = 'PROG' 1.0 0.5 -1.0 ;
*-------------------------------------------
* This to activate in case of constant volume
* combustion
*-------------------------------------------
'SI' LOGEXP ;
DELTARE DELTARY = 'FLAM' 'CREBCOM2' MDOMI PGAZ LMOT1 LCOEF
('REDU' RN MAILCI) ('REDU' YN MAILCI)
('REDU' YININ MAILCI) ('REDU' YFINN MAILCI)
('REDU' K0N MAILCI) ('REDU' DELTAXN MAILCI)
0.0 (1D10 '*' DELTATC) zero ;
RYN = RYN '+' DELTARY ;
YN = RYN '/' RN ;
RETN = RETN '+' DELTARE ;
'FINSI' ;
*-----------------------------------------------
* Names of the components
*-----------------------------------------------
LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2'
'H2O' 'H2IN' 'H2FI' 'K0' ;
****** names of the primitive variables
LISTINCP = 'MOTS' 'RN' 'UX' 'UY' 'PN' 'H2' 'O2'
'H2O' 'H2IN' 'H2FI' 'K0' ;
*--------------------------------------------
LMOT111 = 'MOTS' 'SCAL' ;
GRADR CCC COEFSCAL = 'PENT' MDOMTOT 'CENTRE'
'EULESCAL' 'NOLIMITE' LMOT111 RN ;
LMOT2 = 'MOTS' 'UX' 'UY' ;
GRADV CCC COEFVECT = 'PENT' MDOMTOT 'CENTRE'
'EULEVECT' 'NOLIMITE' LMOT2 GN ;
*--------------------------------------------
* Volumes of the rooms ;
*--------------------------------------------
VV2 = 'DOMA' MDSEC 'VOLUME' ;
CHUN = ('MANUEL' 'CHPO' MAILC2 1 'SCAL' 1.0 'NATU' 'DISCRET') ;
VV2S = ('XTY' CHUN VV2 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
*-------------------------------
VV3 = 'DOMA' MDTHIRD 'VOLUME' ;
CHUN = ('MANUEL' 'CHPO' MAILC3 1 'SCAL' 1.0 'NATU' 'DISCRET') ;
VV3S = ('XTY' CHUN VV3 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
*-------------------------------------------
LISTPS2 = 'PROG' ;
LISTPS3 = 'PROG' ;
LISTTS2 = 'PROG' ;
LISTTS2 = 'PROG' ;
*---- 'LISTE' for the discharge at the coridor
LISCOR1 = 'PROG' ;
*---- 'LISTE' for the averaged velocity at the coridor
LISCOR2 = 'PROG' ;
*---- 'LISTE' for the discharge at the vent
LISVEN1 = 'PROG' ;
*---- 'LISTE' for the averaged velocity at the vent
LISVEN2 = 'PROG' ;
*--------------------------------------------
TPS = 0.0 ;
*------------------------------------------
* Chemical time step
*------------------------------------------
DT_CHEM = SAFFAC '*' DELTATC ;
*-------------------------------------------
**** Temporal loop
*-------------------------------------------
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
'MESSAGE' ('CHAINE' 'SAFFAC = ' SAFFAC) ;
'MESSAGE' ;
*-----------------------------------
'REPETER' BLG ('DIME' LISTT) ;
*----------------------------------
TFINAL = 'EXTRAIRE' &BLG LISTT ;
*----------------------------------
VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
TNM1 = 'COPIER' TN ;
'TEMPS' 'ZERO' ;
'REPETER' BL1 NITER ;
*
**** Primitive variables
*
FREQC = TABC.'FREQ' ;
'SI' (((&BL1 '/' FREQC) '*' FREQC) 'EGA' &BL1) ;
******* The "coridor" **********************
RL2 = 'REDU' RN MSLIN2 ;
GL2 = 'REDU' GN MSLIN2 ;
GLX2 = 'EXCO' 'UX' GL2 ;
CHUN = ('MANUEL' 'CHPO' MSLIN2 1
'SCAL' (DX1) 'NATU' 'DISCRET') ;
*------- rashod discharge ----------------
GLXAV = ('XTY' CHUN GLX2 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
GLXAV = GLXAV '/' H1 ;
LISCOR1 = LISCOR1 'ET' ('PROG' GLXAV) ;
*---- redifining CHUN -------------------
CHUN = ('MANUEL' 'CHPO' MSLIN2 1
'SCAL' (DX1 '*' DX1) 'NATU' 'DISCRET') ;
*------- averaged density ----------------
RLAV = ('XTY' CHUN RL2 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
RLAV = RLAV '/' VOLGOR ;
*------- averaged momentum ---------------
GLXAV = ('XTY' CHUN GLX2 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
GLXAV = GLXAV '/' VOLGOR ;
*-----------------------------------------
* UXL = GLXAV '/' RLAV ;
CHVIT = GLX2 '*' ('INVERSE' RL2) ;
UXL = 'MAXIMUM' CHVIT ;
LISCOR2 = LISCOR2 'ET' ('PROG' UXL) ;
******* The second room **********************
*----------------------------------------------------
* Computation of averaged pressure and temperature
*----------------------------------------------------
RN2 = 'REDU' RN MAILC2 ;
GN2 = 'REDU' GN MAILC2 ;
RETN2 = 'REDU' RETN MAILC2 ;
RYN2 = 'REDU' RYN MAILC2 ;
RSN2 = 'REDU' RSN MAILC2 ;
*--- we average the conservative data
RN2S = ('XTY' RN2 VV2 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
RN2S = RN2S '/' VV2S ;
*----------------------
RETN2S = ('XTY' RETN2 VV2 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
RETN2S = RETN2S '/' VV2S ;
*-----------------------
GNX2 = 'EXCO' 'UX' GN2 'UX' ;
GNY2 = 'EXCO' 'UY' GN2 'UY' ;
*-------------------
GNX2S = ('XTY' GNX2 VV2 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
GNX2S = GNX2S '/' VV2S ;
*-------------------
GNY2S = ('XTY' GNY2 VV2 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
GNY2S = GNY2S '/' VV2S ;
*--------------------
RYNH2 = 'EXCO' 'H2' RYN2 'H2' ;
RYNO2 = 'EXCO' 'O2' RYN2 'O2' ;
RYNH2O = 'EXCO' 'H2O' RYN2 'H2O' ;
*-------------------
RYH2S = ('XTY' RYNH2 VV2 ('MOTS' 'H2')
('MOTS' 'SCAL')) ;
RYH2S = RYH2S '/' VV2S ;
*-------------------
RYO2S = ('XTY' RYNO2 VV2 ('MOTS' 'O2')
('MOTS' 'SCAL')) ;
RYO2S = RYO2S '/' VV2S ;
*-------------------
RYHO2S = ('XTY' RYNH2O VV2 ('MOTS' 'H2O')
('MOTS' 'SCAL')) ;
RYHO2S = RYHO2S '/' VV2S ;
*--------------------------------------
*--------------------------------------
R2S = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' RN2S 'NATU' 'DISCRET') ;
G2S = ('MANUEL' 'CHPO' MAILCI 1 'UX' GNX2S 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCI 1 'UY' GNY2S 'NATU' 'DISCRET') ;
RET2S = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' RETN2S 'NATU' 'DISCRET') ;
RY2S = ('MANUEL' 'CHPO' MAILCI 1 'H2' RYH2S 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCI 1 'O2' RYO2S 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCI 1 'H2O' RYHO2S 'NATU' 'DISCRET') ;
RS2S = ('MANUEL' 'CHPO' MAILCI 1 'H2IN' (aaa) 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCI 1 'H2FI' zero 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCI 1 'K0' K0IN 'NATU' 'DISCRET') ;
*--------------------------------------------
VS PS TS YS SS GAMS = 'PRIM' 'PERFTEMP' PGAZ R2S G2S RET2S RY2S RS2S ;
PN2S = 'MAXIMUM' PS ;
TN2S = 'MAXIMUM' TS ;
*-------------------------------
LISTPS2 = LISTPS2 'ET' ('PROG' PN2S) ;
LISTTS2 = LISTTS2 'ET' ('PROG' TN2S) ;
*-------------------------------------------------
******** IN THE BIGGEST ROOM ********************
RN2 = 'REDU' RN MAILC3 ;
GN2 = 'REDU' GN MAILC3 ;
RETN2 = 'REDU' RETN MAILC3 ;
RYN2 = 'REDU' RYN MAILC3 ;
RSN2 = 'REDU' RSN MAILC3 ;
*--- we average the conservative data
RN2S = ('XTY' RN2 VV3 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
RN2S = RN2S '/' VV3S ;
*----------------------
RETN2S = ('XTY' RETN2 VV3 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
RETN2S = RETN2S '/' VV3S ;
*-----------------------
GNX2 = 'EXCO' 'UX' GN2 'UX' ;
GNY2 = 'EXCO' 'UY' GN2 'UY' ;
*-------------------
GNX2S = ('XTY' GNX2 VV3 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
GNX2S = GNX2S '/' VV3S ;
*-------------------
GNY2S = ('XTY' GNY2 VV3 ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
GNY2S = GNY2S '/' VV3S ;
*--------------------
RYNH2 = 'EXCO' 'H2' RYN2 'H2' ;
RYNO2 = 'EXCO' 'O2' RYN2 'O2' ;
RYNH2O = 'EXCO' 'H2O' RYN2 'H2O' ;
*-------------------
RYH2S = ('XTY' RYNH2 VV3 ('MOTS' 'H2')
('MOTS' 'SCAL')) ;
RYH2S = RYH2S '/' VV3S ;
*-------------------
RYO2S = ('XTY' RYNO2 VV3 ('MOTS' 'O2')
('MOTS' 'SCAL')) ;
RYO2S = RYO2S '/' VV3S ;
*-------------------
RYHO2S = ('XTY' RYNH2O VV3 ('MOTS' 'H2O')
('MOTS' 'SCAL')) ;
RYHO2S = RYHO2S '/' VV3S ;
*--------------------------------------
*--------------------------------------
R2S = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' RN2S 'NATU' 'DISCRET') ;
G2S = ('MANUEL' 'CHPO' MAILCI 1 'UX' GNX2S 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCI 1 'UY' GNY2S 'NATU' 'DISCRET') ;
RET2S = ('MANUEL' 'CHPO' MAILCI 1 'SCAL' RETN2S 'NATU' 'DISCRET') ;
RY2S = ('MANUEL' 'CHPO' MAILCI 1 'H2' RYH2S 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCI 1 'O2' RYO2S 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCI 1 'H2O' RYHO2S 'NATU' 'DISCRET') ;
RS2S = ('MANUEL' 'CHPO' MAILCI 1 'H2IN' (aaa) 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCI 1 'H2FI' zero 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' MAILCI 1 'K0' K0IN 'NATU' 'DISCRET') ;
*--------------------------------------------
VS PS TS YS SS GAMS = 'PRIM' 'PERFTEMP' PGAZ R2S G2S RET2S RY2S RS2S ;
PN2S = 'MAXIMUM' PS ;
TN2S = 'MAXIMUM' TS ;
*-------------------------------------------------
LISTPS3 = LISTPS3 'ET' ('PROG' PN2S) ;
LISTTS3 = LISTTS3 'ET' ('PROG' TN2S) ;
'FINSI' ;
*--------------------------------------------------
*-------------------------------------------
VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN
TNM1 ;
*--------------------------------------------
* Boundary conditions
*--------------------------------------------
*** find new CP and CV *******
* CP1 CV1 = CPCV tg PGAZ ;
*------------------------------------------
'SI' LOGBC ;
'SI' (&BLG 'EGA' 1) ;
SHIN = 'EXCO' SN 'H2IN' ;
SHFI = 'EXCO' SN 'H2FI' ;
SHK0 = 'EXCO' SN 'K0' ;
RCHLIM RCHRES = 'KONV' 'VF' 'PERFMULT' 'CLIM' 'RESI'
MDOMTOT MFRONTH PGAZB LISTINCO LISTINCP
RN VN PN YN SHIN SHFI SHK0
(TABLIM . 'CHPOUT') 'OUTP' ;
'FINSI' ;
*------------------------------------
MAILLIM = 'EXTRAIRE' RCHLIM 'MAILLAGE' ;
MAIL111 = 'EXTRAIRE' RCHRES 'MAILLAGE' ;
*------------------------------------
NBB = 'NBEL' MAIL111 ;
TTT = 'REDU' TN MAIL111 ;
UNIF = 'MANUEL' CHPO ('DOMA' MDOMTOT 'CENTRE')
1 'SCAL' 1.0 ;
UNI1 = 'REDU' UNIF MAIL111 ;
NUNI = UNI1 '*' (1.0 '/' NBB) ;
**** computing average temperature at the boundary
TAV = 'XTY' TTT NUNI ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
*** find new CP and CV *******
CPN CVN = CPCV TAV PGAZ ;
*** Updating CP and CV for the boundary conditions
PGAZB . 'CP' . 'H2 ' = CPN . 'H2' ;
PGAZB . 'CP' . 'H2O ' = CPN . 'H2O' ;
PGAZB . 'CP' . 'O2 ' = CPN . 'O2' ;
PGAZB . 'CP' . 'N2 ' = CPN . 'N2' ;
*----------------------------------------
PGAZB . 'CV' . 'H2 ' = CVN . 'H2' ;
PGAZB . 'CV' . 'H2O ' = CVN . 'H2O' ;
PGAZB . 'CV' . 'O2 ' = CVN . 'O2' ;
PGAZB . 'CV' . 'N2 ' = CVN . 'N2' ;
* 'FINSI' ;
*--------------------------------
SHIN = 'EXCO' SN 'H2IN' ;
SHFI = 'EXCO' SN 'H2FI' ;
SHK0 = 'EXCO' SN 'K0' ;
RCHLIM RCHRES = 'KONV' 'VF' 'PERFMULT' 'CLIM' 'RESI'
MDOMTOT MFRONTH PGAZB LISTINCO LISTINCP
RN VN PN YN SHIN SHFI SHK0
(TABLIM . 'CHPOUT') 'OUTP' ;
MAILLIM = 'EXTRAIRE' RCHLIM 'MAILLAGE' ;
MAIL111 = 'EXTRAIRE' RCHRES 'MAILLAGE' ;
*--------------------------------------------
* We compute the debit et vitesse
*--------------------------------------------
'SI' (((&BL1 '/' FREQC) '*' FREQC) 'EGA' &BL1) ;
VOLUP = 'REDU' ('DOMA' MDOMTOT 'VOLUME')
MAIL111 ;
RLIM1 = 'EXCO' 'RN' RCHRES ;
DEBUP = 'XTY' VOLUP RLIM1 ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
DEBUP = (DEBUP '/' LW) '*' (-1.0) ;
LISVEN1 = LISVEN1 'ET' ('PROG' DEBUP) ;
*--------------------------------------------
RNUP = 'REDU' RN MAIL111 ;
GNUP = 'REDU' GN MAIL111 ;
GXUP = 'EXCO' 'UY' GNUP ;
CUN = ('MANUEL' 'CHPO' MAIL111 1
'SCAL' 1.0D0 'NATU' 'DISCRET') ;
*------------------------
GXUPS = 'XTY' GXUP VOLUP ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
RNUPS = 'XTY' RNUP VOLUP ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
VITUP = GXUPS '/' RNUPS ;
LISVEN2 = LISVEN2 'ET' ('PROG' VITUP) ;
'FINSI' ;
*--------------------------------------------
'FINSI' ;
*------------------------------------
*--------------------------------
TNM1 = 'COPIER' TN ;
'SI' LOGSO ;
GRADR ALR = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') RN 'GRADGEO' COEFSCAL ;
GRADP ALP = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PN 'GRADGEO' COEFSCAL ;
GRADV ALV = 'PENT' MDOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY') VN 'GRADGEO' COEFVECT ;
GRADY ALY = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'H2' 'O2' 'H2O') YN 'GRADGEO' COEFSCAL ;
GRADS ALS = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'H2IN' 'H2FI' 'K0') SN 'GRADGEO' COEFSCAL ;
ROF VITF PF YF SF = 'PRET' 'PERFTEMP' 2 1
MDOMTOT PGAZ
RN GRADR ALR
VN GRADV ALV
PN GRADP ALP
YN GRADY ALY
SN GRADS ALS ;
'SINON' ;
ROF VITF PF YF SF = 'PRET' 'PERFTEMP' 1 1
MDOMTOT PGAZ
RN
VN
PN
YN
SN ;
'FINSI' ;
*-----------------------------------------------
'SI' LOGBC ;
RESIDU DELTAT = 'KONV' 'VF' 'PERFTEMP' 'RESI' METO
MDOMTOT PGAZ LISTINCO ROF VITF PF YF SF MAILLIM ;
RESIDU = RESIDU 'ET' RCHRES ;
'SINON' ;
RESIDU DELTAT = 'KONV' 'VF' 'PERFTEMP' 'RESI' METO
MDOMTOT PGAZ LISTINCO ROF VITF PF YF SF ;
'FINSI' ;
*-----------------------------------------------
DT_CON = SAFFAC '*' DELTAT ;
*
**** The time step linked to tps
*
DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
*
**** Total time step
*
DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS DT_CHEM ) ;
*
**** Increment of the variables (convection)
*
RESIDU = DTMIN '*' RESIDU ;
DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU ('MOTS' 'UX' 'UY') ;
DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
DRYN = 'EXCO' ('MOTS' 'H2' 'O2' 'H2O') RESIDU
('MOTS' 'H2' 'O2' 'H2O') ;
DRSN = 'EXCO' ('MOTS' 'H2IN' 'H2FI' 'K0') RESIDU
('MOTS' 'H2IN' 'H2FI' 'K0') ;
TPS = TPS '+' DTMIN ;
RN = RN '+' DRN ;
GN = GN '+' DGN ;
RETN = RETN '+' DRETN ;
RYN = RYN '+' DRYN ;
RSN = RSN '+' DRSN ;
YN = RYN '/' RN ;
SN = RSN '/' RN ;
*
**** Increment of the variables (chemical source)
*
K0N = 'EXCO' 'K0' SN 'SCAL' ;
YININ = 'EXCO' 'H2IN' SN 'H2' ;
YFINN = 'EXCO' 'H2FI' SN 'H2' ;
DRETN DRYN = 'FLAM' 'CREBCOM2' MDOMTOT PGAZ LMOT1 LCOEF RN YN
YININ YFINN K0N DELTAXN EPS_CC DTMIN (zero '*' 1.0D3) ;
RYN = RYN '+' DRYN ;
RETN = RETN '+' DRETN ;
*---------------------------------------------------
* Rate of loosing energy due to boundary conditions
*---------------------------------------------------
* 'SI' (TPS '>' 3.0) ;
VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
HH = TABLIM . 'H' ;
PERT = HH '*' (TN '-' (TABLIM . 'TNUL')) ;
RETN = RETN '-' (DTMIN '*' PERT) ;
* 'FINSI' ;
*--------------------------------------------------
* Gravity forces
*--------------------------------------------------
LIMOT1 = 'MOTS' 'SCAL' 'SCAL' ;
LIMOT2 = 'MOTS' 'UX' 'UY' ;
LIMOT3 = 'MOTS' 'UX' 'UY' ;
ROG = RN '*' (TABLIM . 'GRAV') LIMOT1 LIMOT2 LIMOT3 ;
GN = GN '+' (DTMIN '*' ROG) ;
ROGU = 'PSCAL' ROG VN LIMOT3 LIMOT2 ;
RETN = RETN '+' (DTMIN '*' ROGU) ;
*---------------------------------------------------
*---------------------------------------------
* Capturing variables at the capteurs
*---------------------------------------------
'SI' (((&BL1 '/' FREQC) '*' FREQC) 'EGA' &BL1) ;
CENTOT = 'DOMA' MDOMTOT 'CENTRE';
PTS = TABC . 'LPOINTS' ;
***** M1 is the number of capteurs
M1 = 'NBEL' PTS ;
LTPS = TABC . 'LTPS' ;
NTPS = 'DIME' LTPS ;
**** N1 is the number of variables captured
L1 = TABC . 'COMP' ;
N1 = 'DIME' L1 ;
* On etudie uniquement la pression
* peut evoluer par la suite !!!
LOG1 = 'EXISTE' TABC 'NOMC' ;
* Initialisation des tables
'SI' (NTPS EGA 0 );
'REPETER' BOU1 M1;
**** extract the names of the capteurs
PT1 = PTS 'POIN' &BOU1;
'SI' (LOG1 'ET' (('DIME' (TABC.'NOMC')) 'EGA' M1 )) ;
CAPI = 'EXTR' (TABC.'NOMC') &BOU1;
'SINON' ;
CAPI = 'CHAINE' 'CAPT' &BOU1;
'FINSI';
PT2 = CENTOT 'POIN' 'PROCHE' PT1;
'SI' (&BOU1 'EGA' 1);
TABC.'RPOINTS' = PT2 ;
'SINON';
TABC.'RPOINTS' = TABC .'RPOINTS' 'ET' PT2;
'FINSI';
TABC . CAPI = 'TABLE';
'REPETER' BOU2 N1;
MOI = EXTR L1 &BOU2 ;
TABC . CAPI . MOI = 'PROG';
'FIN' BOU2;
'FIN' BOU1;
'FINSI';
*--------------------------------------------
'REPETER' BOU1 M1 ;
PT1 = TABC.'RPOINTS' 'POIN' &BOU1;
'SI' (LOG1 'ET' (('DIME' (TABC.'NOMC')) 'EGA' M1 )) ;
CAPI = 'EXTR' (TABC.'NOMC') &BOU1;
'SINON' ;
CAPI = 'CHAINE' 'CAPT' &BOU1;
'FINSI';
'REPETER' BOU2 N1;
MOI = EXTR L1 &BOU2 ;
*** captage de la pression
'SI' ('EGA' MOI 'PN');
X1 = 'EXTR' PN 'SCAL' PT1 ;
* 'LISTE' X1 ;
TABC . CAPI . MOI = (TABC . CAPI . MOI) 'ET' ('PROG' X1);
'FINSI';
*** captage de temperature
'SI' ('EGA' MOI 'TN');
X1 = 'EXTR' TN 'SCAL' PT1 ;
TABC . CAPI . MOI = TABC . CAPI . MOI 'ET' ('PROG' X1);
'FINSI';
'FIN' BOU2;
'FIN' BOU1;
'FINSI' ;
*****---- creating listreel of time steps
'SI' (((&BL1 '/' FREQC) '*' FREQC) 'EGA' &BL1) ;
TABC . 'LTPS' = (TABC . 'LTPS') 'ET' ('PROG' TPS);
'FINSI' ;
*---------------------------------------------
'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
'MESSAGE' ('CHAINE' 'ITER =' &BL1 ' TPS =' TPS) ;
'FINSI' ;
*-------------------------------------------------------
'SI' (TPS '>EG' TFINAL) ;
'QUITTER' BL1 ;
'FINSI' ;
'FIN' BL1 ;
*------------------------------------
*-----------------------------------
'FIN' BLG ;
'TEMPS' ;
*************************************************
* Posttreatment *
*************************************************
LISTP0 = 'PROG' 10 '*' 100000.0D0 ;
OVPRES = (TABC . 'C1' . 'PN') '-' LISTP0 ;
EVERP1 = 'EVOL' 'MANU' 'Niter' (TABC . 'LTPS' ) 'PC1'
OVPRES ;
*------------------------------------------------
LPGOOD = 'PROG' 1.3086 21.753 43.070 41.450 30.604 26.187
32.615 43.088 48.922 47.665 ;
DIFFER = OVPRES '-' LPGOOD ;
'SI' (('MAXIMUM' DIFFER) <EG 1.D-3) ;
'MESSAGE' ;
'MESSAGE' 'OVERPRESSURE IS OK' ;
'MESSAGE' ;
'SINON' ;
'MESSAGE' ;
'MESSAGE' 'THERE IS A PROBLEM' ;
'MESSAGE' ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer