Télécharger la plaquette Cast3M


fr en

Télécharger rdem_surf1Daxi.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier : rdem_surf1Daxi.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ****************************************************************
  5. ****    APPROCHE VF "Cell-Centred Formulation" pour la      ****
  6. ****    combustion.                                         ****
  7. ****    MODELE RDEM                                         ****
  8. ****                                                        ****
  9. ****    OPERATEURS 'PRIM', PRET, KONV                       ****
  10. ****    PROPAGATION D'UNE FLAMME DANS UN TUBE               ****
  11. ****    We verify that S_flame = S_tube                     ****
  12. ****                                                        ****
  13. ****    S. KUDRIAKOV  SFME/LTMF                             ****
  14. ****************************************************************
  15. *
  16.  'OPTION' 'ECHO' 1 'DIME' 2
  17.           'ELEM' 'QUA4' 'TRAC' 'X' ;
  18.  'OPTION' 'MODE' 'AXIS' ;
  19.  
  20.  RAF = 1 ;
  21. ************************************************
  22. *     Constructing mesh                        *
  23. ************************************************
  24. *---------------------------------
  25. ******** refinement level ********
  26. *---------------------------------
  27.  EPSI = 0.154 '/' 10000.0D0 ;
  28. *----------------------------------
  29. **  radius ************************
  30. *----------------------------------            
  31.  R = 0.15 '/' 2.0D0 ;
  32. *-----length of the obstacle-------- 
  33.  LD = 0.03015D0 ;
  34. *--------------------------------------------------------
  35. **** number of points along the 'REST' of the radius 
  36. *--------------------------------------------------------
  37.  NRAD = 4 '*' RAF ;          
  38. ******* base of the tube
  39. *!!!!!!!!!!!!!!!!!!!!
  40.  A0 = 0.0 0.0 ;
  41.  A1 = R 0.0  ;
  42. 'SI' ('EGA' ('VALEUR' 'MODE') 'AXIS') ;
  43.   A0 = EPSI 0.0 ;
  44. 'FINSI' ;  
  45. *********************************
  46.   DX = R '/' NRAD ;
  47. **************************************************** 
  48. **************************************************** 
  49.  A0A1 = 'DROIT' NRAD A0 A1 ;
  50. *----------------------------------------
  51.  IBASE = A0A1 ;
  52.  
  53. * 'TRACER' IBASE ;
  54. *-------------------------------------------------------
  55. ****** Length of the first section (without obstacles)
  56. *-------------------------------------------------------
  57.  LFIRST = 3.0 ;
  58.  NFIRST = 'ENTIER' ((LFIRST '/' R) '*' NRAD) ;  
  59. *---------------------------- 
  60.  B0 = A0 'PLUS' (0.0 LFIRST) ;
  61.  B1 = A1 'PLUS' (0.0 LFIRST) ;
  62. *--------------------------
  63.  A0B0 = 'DROIT' NFIRST A0 B0 ;
  64.  B0B1 = 'DROIT' NRAD B0 B1 ;
  65.  B1A1 = 'DROIT' NFIRST B1 A1 ;
  66. *------------------------------
  67.  DOMTOT = 'DALLER' A0B0 B0B1 B1A1 ('INVERSE' A0A1) 'PLANE' ;
  68. *********************************************************
  69. ***** Extracting ignition region  ***********************
  70. *********************************************************
  71.  LRIG = R ;
  72. *----------------
  73.  YYAXI = 'COORDONNEE' 2 DOMTOT ;
  74.  YY_I = YYAXI 'POIN'  'EGINFE' LRIG ;
  75.  TIGN = DOMTOT 'ELEM'  'APPUYE' 'STRICTEMENT' YY_I ;
  76.  
  77. * 'TRACER' TIGN ;
  78.  
  79.  'OUBLIER'  YYAXI ;
  80.  'OUBLIER'  YY_I ;
  81. *************************************************
  82. ***** adapting the names ************************
  83. *************************************************
  84.     DOM1 = TIGN ;
  85.     DOM2 =  'DIFF'  DOMTOT DOM1 ;
  86.     DOMLIM = 'CONTOUR' DOMTOT ;
  87.     DOMTOT = DOMTOT 'COULEUR' 'VERT' ;
  88. *********************************************
  89. **** Extracting a line for posttreatment ****
  90. *********************************************
  91.  'OPTION' 'ELEM' 'SEG2' ;
  92.   vev1 = (0.3 '*' DX) (0.25 '*' DX) ;
  93.   vev2 = (0.3 '*' DX) (LFIRST '-'(0.25 '*' DX)) ;  
  94.   vev3 = (1.0 '*' DX) (LFIRST '-'(0.25 '*' DX)) ;
  95.   vev4 = (1.0 '*' DX) (0.25 '*' DX) ;  
  96.  LIMCON = 'MANUEL' 'QUA4'
  97.     (A0 'PLUS' vev1)
  98.     (A0 'PLUS' vev2)
  99.     (A0 'PLUS' vev3)
  100.     (A0 'PLUS' vev4) ;
  101.  LIMCON = LIMCON 'COULEUR' 'ROUG' ;
  102.  
  103. ***************************************************
  104. ****** Making models !!!!!!!!!!!!!!!!! ************
  105. ***************************************************
  106.      MDNS = 'EULER' ;
  107. *---------------------------------
  108.  $DOMTOT = 'MODELISER' DOMTOT 'EULER' ;
  109.  $DOMLIM = 'MODELISER' DOMLIM 'EULER' ;
  110.  $DOM1   = 'MODELISER' DOM1 'EULER' ;
  111.  $DOM2   = 'MODELISER' DOM2 'EULER' ;
  112. * 
  113.  TDOMTOT = 'DOMA' $DOMTOT 'VF' ;
  114.  TDOMLIM = 'DOMA' $DOMLIM 'VF' ;
  115.  TDOM1   = 'DOMA' $DOM1   'VF' ;
  116.  TDOM2   = 'DOMA' $DOM2   'VF' ;
  117. * 
  118.  QDOMTOT = TDOMTOT  . 'QUAF' ;
  119.  QDOMLIM = TDOMLIM  . 'QUAF' ;
  120.  QDOM1   = TDOM1    . 'QUAF' ;
  121.  QDOM2   = TDOM2    . 'QUAF' ;
  122. *
  123.  'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOMLIM ;
  124.  'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOM1 ;
  125.  'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOM2 ;
  126. *--------------------------------------------
  127. *--------------------------------------------
  128.  
  129.  LIGEVO = (TDOMTOT . 'CENTRE') 'INCL' LIMCON ;
  130. *  P1 = 'POIN' (COOR 1 LIGEVO) 'MINIMUM' ;
  131. *  'ORDONNER' LIGEVO P1 ;
  132.  'ORDONNER' LIGEVO ;
  133.  LIGEVO = 'QUELCONQUE' 'SEG2' LIGEVO 'COULEUR' 'VERT';
  134. *************************************************
  135. *******   Volumes !!!!!!!!!!!   *****************
  136. *************************************************    
  137. ***---  Volume of the total domain
  138.  VTOT = 'DOMA' $DOMTOT 'VOLUME' ;
  139.  CENDOM = 'DOMA'  $DOMTOT 'CENTRE' ;
  140.  CHUN = ('MANUEL' 'CHPO' CENDOM 1 'SCAL' 1.0 'NATU' 'DISCRET') ;
  141.  VTOT = ('XTY' CHUN VTOT ('MOTS' 'SCAL')
  142.        ('MOTS' 'SCAL')) ;
  143. **************************************************
  144. ******* Pressure transducers *********************
  145. **************************************************
  146.  PBAS = A0 ;
  147.  PTOP = B0 ;
  148. *--------------------
  149.  PC1 =  ('DOMA' $DOMTOT 'CENTRE') 'POIN' 'PROC' PBAS ;
  150.  PC2 =  ('DOMA' $DOMTOT 'CENTRE') 'POIN' 'PROC' PTOP ;
  151.  
  152.  PCTOT = PC1 'ET' PC2 ;
  153. *-------------------------------------
  154.  GRAPH = VRAI ;
  155. * GRAPH = FAUX ;
  156. * Upwind scheme
  157.  METO = 'SS' ;
  158. ***************************************** 
  159. **** vitese fondamentale   **************
  160. *****************************************
  161.  K0 = 10.0D0 ;
  162. ***********************************************
  163. **** Computational parameters *****************
  164. ***********************************************
  165.  NITER = 10000000 ;
  166.  TFINAL = 1.0D-3 ;
  167.  SAFFAC = 1.0 ;
  168.  LOGSO = VRAI ;
  169. *************************************
  170. **** table for pressure transducers *
  171. *************************************
  172.  PCEN = 'TABLE' ;
  173.  PCEN.'TIME' = 'PROG' ;
  174.  PCEN.'PRESSION' = 'TABLE' ;
  175.  NPP = 'NBEL' PCTOT ;
  176. ***** two pressure transducers 
  177.  PCEN.'PRESSION'. 1 = 'PROG' ;
  178.  PCEN.'PRESSION'. 2 = 'PROG' ;
  179. *******************************************
  180. *** table for velocity along centerline  **
  181. *******************************************
  182.  VCEN = 'TABLE' ;
  183.  VCEN . 'TIME' = 'PROG' ;
  184.  IPOST = 0 ;
  185. **** PRIMCONS ****
  186. ******************
  187.  'DEBPROC' PRIMCONS ;
  188.  'ARGUMENT' PGAS*TABLE TN1*'CHPOINT ' TN2*'CHPOINT '
  189.    PN1*'CHPOINT ' PN2*'CHPOINT ' VN1*'CHPOINT ' VN2*'CHPOINT ' ;
  190. *
  191. * ETHER = int_0^T cv(T') dT'           T <  TMAX
  192. *       = int_0^TMAX cv(T') dT' '+'
  193. *         cv(TMAX)                     T >= TMAX
  194.  ESP1 = 'EXTRAIRE' (PGAS . 'SPECIES') 1 ;
  195. * DY1 = y_i - y_f for the species 1 
  196.  DY1 = (('EXTRAIRE' (PGAS . 'MASSFRA') 1) '-'
  197.    ('EXTRAIRE' (PGAS . 'MASSFRA') 2)) ;
  198.  COEF1 = ('EXTRAIRE' (PGAS . 'CHEMCOEF') 1) '*'
  199.     (PGAS . ESP1 . 'W') ;
  200.  YFINPH1 = 1.0 ;
  201.  YFINPH2 = 1.0 ;
  202.  'SI' (COEF1 > 0) ;
  203.     YPH2 = 'EXTRAIRE' (PGAS . 'MASSFRA') 2 ;
  204.     YPH1 = YPH2 '+' DY1 ;
  205.  'SINON' ;
  206.     YPH1 = 'EXTRAIRE' (PGAS . 'MASSFRA') 2 ;
  207.     YPH2 = YPH1 '-' DY1 ;
  208.  'FINSI' ;
  209.  YFINPH1 = YFINPH1 '-' YPH1 ;
  210.  YFINPH2 = YFINPH2 '-' YPH2 ;
  211.  PRYPH1 = 'PROG' YPH1 ;
  212.  PRYPH2 = 'PROG' YPH2 ;
  213.  'REPETER' BLESP (('DIME' (PGAS . 'SPECIES')) '-' 2) ;
  214.     ESP = 'EXTRAIRE' (PGAS . 'SPECIES') (&BLESP '+' 1)  ;
  215.     COEF = ('EXTRAIRE' (PGAS . 'CHEMCOEF') (&BLESP '+' 1)) 
  216.        '*' (PGAS . ESP . 'W') ;
  217.     DY = (DY1 * (COEF '/' COEF1)) ;
  218.     'SI' (COEF > 0) ;
  219.        YPH2 = 'EXTRAIRE' (PGAS . 'MASSFRA') (&BLESP '+' 2) ;
  220.        YPH1 = YPH2 '+' DY ;
  221.     'SINON' ;
  222.        YPH1 = 'EXTRAIRE' (PGAS . 'MASSFRA') (&BLESP '+' 2) ;
  223.        YPH2 = YPH1 '-' DY ;
  224.     'FINSI' ;
  225.     PRYPH1 = PRYPH1 'ET' ('PROG' YPH1) ;
  226.     PRYPH2 = PRYPH2 'ET' ('PROG' YPH2) ;
  227.     YFINPH1 = YFINPH1 '-' YPH1 ;
  228.     YFINPH2 = YFINPH2 '-' YPH2 ;
  229.  'FIN' BLESP ;
  230.  PRYPH1 = PRYPH1 'ET' ('PROG' YFINPH1) ; 
  231.  PRYPH2 = PRYPH2 'ET' ('PROG' YFINPH2) ;
  232. * 'LISTE' PRYPH1 ;
  233. * 'LISTE' PRYPH2 ; 
  234. *
  235.  TMAX = (PGAS . 'TMAX') ;
  236. * TCAL1 = MIN TN1, TMAX 
  237.  TCAL1 = 0.5D0 '*' ((TMAX '+' TN1) '-' ('ABS' (TN1 '-' TMAX))) ;
  238.  DTN1 = TN1 '-' TCAL1 ;
  239. * TCAL1 = MIN TN1, TMAX 
  240.  TCAL2 = 0.5D0 '*' ((TMAX '+' TN2) '-' ('ABS' (TN2 '-' TMAX))) ;
  241.  DTN2 = TN2 '-' TCAL2 ;
  242. *
  243. * Internal energy (J/kg in SI)
  244. *
  245.  ETHER1 = 0.0 ;
  246.  CV1 = 0.0 ;
  247.  ETHER2 = 0.0 ;
  248.  CV2 = 0.0 ;
  249.  FUNTN1 = 1.0 ;
  250.  FUNTN2 = 1.0 ;
  251.  'REPETER' BLPO ((PGAS . 'NORD') '+' 1) ;
  252.     'REPETER' BLESP ('DIME' (PGAS . 'SPECIES')) ;
  253.        ESP = 'EXTRAIRE' (PGAS . 'SPECIES') &BLESP ;
  254.        YCEL1 = 'EXTRAIRE' PRYPH1 &BLESP ;
  255.        YCEL2 = 'EXTRAIRE' PRYPH2 &BLESP ;
  256.        AA = 'EXTRAIRE' (PGAS . ESP . 'A') &BLPO ;
  257.        DCV1 = (AA * YCEL1 * FUNTN1) ;
  258.        DCV2 = (AA * YCEL2 * FUNTN2) ;
  259.        CV1 = CV1 '+' DCV1 ;
  260.        CV2 = CV2 '+' DCV2 ;
  261.        ETHER1 = ETHER1 '+' (DCV1 * TCAL1 '/' (&BLPO)) ;
  262.        ETHER2 = ETHER2 '+' (DCV2 * TCAL2 '/' (&BLPO)) ;
  263.     'FIN' BLESP ;
  264.     FUNTN1 = FUNTN1 '*' TCAL1 ;
  265.     FUNTN2 = FUNTN2 '*' TCAL2 ;
  266.  'FIN' BLPO ;
  267.  ETHER1 = ETHER1 '+' (CV1 '*' DTN1) ;
  268.  ETHER2 = ETHER2 '+' (CV2 '*' DTN2) ;
  269. *
  270. * Formation energy/enthalpy (J/kg in SI) and gas constant (J/kg/K)
  271. *
  272.  EFORM1 = 0.0 ;
  273.  EFORM2 = 0.0 ;
  274.  RGAS1 = 0.0 ;
  275.  RGAS2 = 0.0 ;
  276.  'REPETER' BLESP ('DIME' (PGAS . 'SPECIES')) ;
  277.     ESP = 'EXTRAIRE' (PGAS . 'SPECIES') &BLESP ;
  278.     YCEL1 = 'EXTRAIRE' PRYPH1 &BLESP ;
  279.     YCEL2 = 'EXTRAIRE' PRYPH2 &BLESP ;
  280.     EFORM1 = EFORM1 '+' (YCEL1 * (PGAS . ESP . 'H0K')) ;
  281.     EFORM2 = EFORM2 '+' (YCEL2 * (PGAS . ESP . 'H0K')) ;
  282.     RGAS1 = RGAS1 '+' (YCEL1 *  (PGAS . 'RUNIV') '/' 
  283.       (PGAS . ESP . 'W')) ;
  284.     RGAS2 = RGAS2 '+' (YCEL2 *  (PGAS . 'RUNIV') '/' 
  285.       (PGAS . ESP . 'W')) ;
  286.  'FIN' BLESP ;
  287. *
  288. * Computation of the conservative variables
  289. *
  290.  RN1 = PN1 '/' (RGAS1 '*' TN1) ;
  291.  RN2 = PN2 '/' (RGAS2 '*' TN2) ;
  292.  GN1 = RN1 * VN1 ;
  293.  GN2 = RN2 * VN2 ;
  294.  LVEL = 'MOTS' 'UX' 'UY' 'UZ' ;
  295.  ECIN1 = 0.5D0 '*' ('PSCAL' VN1 VN1 LVEL LVEL) ;
  296.  ECIN2 = 0.5D0 '*' ('PSCAL' VN2 VN2 LVEL LVEL) ;
  297.  RETN1 = RN1 '*' (ETHER1 '+' ECIN1 '+' EFORM1) ; 
  298.  RETN2 = RN2 '*' (ETHER2 '+' ECIN2 '+' EFORM2) ;
  299. *
  300.  'RESPRO' RN1 RN2 GN1 GN2 RETN1 RETN2 ;
  301.  'FINPROC' ;
  302. ***************************************************
  303. ********* end of personal subroutines *************
  304. ***************************************************
  305.  
  306.  
  307. ****************************************************************
  308. ****************************************************************
  309. *****      Initial conditions and gas properties.          *****
  310. ****************************************************************
  311. ****************************************************************
  312. *
  313. *************************************************
  314. **** The table for the properties of the gas ****
  315. *************************************************
  316. *
  317.  PGAS = 'TABLE' ;
  318. *
  319. **** Order of the polynomial order for cv = cv(T) 
  320. *    For T > TMAX, cv(T) = cv(Tmax)
  321. * 
  322.  PGAS . 'TMAX' = 6000.0 ;
  323.  PGAS . 'NORD' = 4 ;
  324. *
  325. **** Species involved in the mixture (before or after
  326. *    the chemical reaction)
  327. *
  328.  PGAS . 'SPECIES' = 'MOTS' 'H2  ' 'O2  ' 'H2O ' 'N2  ' ;
  329. *
  330. *
  331. **** Coefficient of the chemical reaction.
  332. *    Note that for the first species this coefficient should be positive
  333. *    Normal, we take it equal to 1.
  334. *
  335. *    H2 '+' 0.5 O2 ---> H2O
  336. *
  337.  PGAS . 'CHEMCOEF' = 'PROG' 1.0 0.5 -1.0 0.0 ;
  338. * 
  339. **** Mass fraction of the first species before and after the combustion
  340. *    Final mass fractions of the species with positive coefficients.
  341. *    Final mass fractions of the species with non-positive coefficient.
  342. *    The mass fraction of the last species is not given.
  343. *  meme que 'species' YH2_ini, YH2_fin, YO2_fin, YH2O_ini 
  344.  PGAS . 'MASSFRA' = 'PROG' 1.03333E-02  0.964039E-11  1.48501E-01
  345.                            0.127442E-10 ;
  346. *
  347. **** Coef with the gas properties
  348. *
  349.  PGAS .  'H2  ' = 'TABLE'  ;
  350.  PGAS .  'H2O ' = 'TABLE'  ;
  351.  PGAS .  'N2  ' = 'TABLE'  ;
  352.  PGAS .  'O2  ' = 'TABLE'  ;
  353. *
  354. **** Runiv (J/mole/K)
  355. *
  356.  PGAS . 'RUNIV' = 8.31441 ;
  357. *
  358. **** W (kg/mole). Gas constant (J/kg/K = Runiv/W)
  359. *
  360.  PGAS .  'H2  ' . 'W' =  2. * 1.00797E-3 ;
  361.  PGAS .  'O2  ' . 'W' =  2. * 15.9994E-3 ;
  362.  PGAS .  'H2O ' . 'W' = (PGAS .  'H2  ' . 'W' ) '+'
  363.     (0.5 * (PGAS .  'O2  ' . 'W' )) ;
  364.  PGAS .  'N2  ' . 'W' = 2 * 14.0067E-3 ;
  365. *
  366. **** Polynomial coefficients 
  367. *
  368.  PGAS .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836
  369.                                -2.37281455E-07 1.84701105E-11 ;
  370.  PGAS .  'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
  371.                               -1.82753232E-08 2.44485692E-12 ;
  372.  PGAS .  'N2  ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
  373.                               8.78233606E-09 -3.05514485E-13 ;
  374.  PGAS .  'O2  ' . 'A' = 'PROG' 575.012333  0.350522002 -0.000128294865
  375.                               2.33636971E-08 -1.53304905E-12;
  376. *
  377. **** Formation enthalpies (energies) at 0K (J/Kg)
  378. *
  379.  PGAS .  'H2  ' . 'H0K' = -4.195D6 ;
  380.  PGAS .  'H2O ' . 'H0K' = -1.395D7 ;
  381.  PGAS .  'N2  ' . 'H0K' = -2.953D5 ;
  382.  PGAS .  'O2  ' . 'H0K' = -2.634D5 ;
  383. *
  384. *************************************************
  385. **** The initial conditions *********************
  386. *************************************************
  387. *
  388.  eps = 1.0D-64 ;
  389. * 
  390.  T1   = 300.0D0 ;
  391. * alph11 = 1.0 '-' 1.0D-5 ;
  392.  alph11 = 1.0D-5 ;
  393.  alph12 =  1.0 '-' 1.0D-5 ;
  394.  alph21 = 1.0 '-' alph11 ;
  395.  alph22 = 1.0 - alph12 ;
  396.  T2   = 1325.1D0 ;
  397.  un1  = 0.0 ;
  398.  un2  = 0.0 ;
  399.  ut1  = 0.0 ;
  400.  ut2  = 0.0 ;
  401.  pre1 = 1.D5 ;
  402.  pre2 = 1.D5 ;
  403. *!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  404. *** these lists are needed for flame evolution with time
  405. *** along the centerline
  406.  DFLAM = 'PROG' LRIG  ;
  407.  DFLAM1 = 'PROG'  LRIG ;
  408.  TFLAM = 'PROG' 0.0 ;
  409.  TFLAM1 = 'PROG'  0.0 ;
  410. *!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  411. *
  412.  CHVN1  = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 2 'UX' un1
  413.    'UY' ut1) ;
  414.  CHVN2  = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 2 'UX' un2
  415.    'UY' ut2) ;
  416. *
  417.  CHTN1  = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' T1) ;
  418.  CHTN2  = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' T2) ;
  419. *
  420.  CHPN1  = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' pre1) ;
  421.  CHPN2  = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' pre2) ;
  422. ********************************************************
  423. ********************************************************
  424.  CHAL1  = ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL'
  425.    alph11) '+'
  426.    ('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 1 'SCAL'
  427.    alph12) ;
  428.  CHAL2  = ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL'
  429.    alph21) '+'
  430.    ('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 1 'SCAL'
  431.    alph22) ;
  432. ********************************************************   
  433. ********************************************************
  434.  CHRN1 CHRN2 CHGN1 CHGN2 CHRET1 CHRET2 = PRIMCONS
  435.    PGAS CHTN1 CHTN2 CHPN1 CHPN2 CHVN1 CHVN2 ;
  436. * 
  437.  RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS CHAL1 CHAL2
  438.    (CHAL1 * CHRN1) (CHAL2 * CHRN2) (CHAL1 * CHGN1)
  439.    (CHAL2 * CHGN2) (CHAL1 * CHRET1) (CHAL2 * CHRET2)
  440.     CHTN1 CHTN2 EPS ;
  441.  TN1M = COPIER TN1 ;
  442.  TN2M = COPIER TN2 ;
  443. *
  444.  ERRO = 'MAXIMUM' (PRE1 '-' ('REDU' PN1 (TDOM1 . 'CENTRE'))) 'ABS' ;
  445.  ERRO = ERRO '/' (pre1) ;
  446.  'SI' (ERRO > 1.0D-6) ;
  447.      'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
  448.      'ERREUR' 21 ;
  449.  'FINSI' ; 
  450.  ERRO = 'MAXIMUM' (PRE2 '-' ('REDU' PN2 (TDOM2 . 'CENTRE'))) 'ABS' ;
  451.  ERRO = ERRO '/' (pre2) ;
  452.  'SI' (ERRO > 1.0D-6) ;
  453.      'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
  454.      'ERREUR' 21 ;
  455.  'FINSI' ;
  456.  ERRO = 'MAXIMUM' (T1 '-' ('REDU' TN1 (TDOM1 . 'CENTRE'))) 'ABS' ;
  457.  ERRO = ERRO '/' (T1) ;
  458.  'SI' (ERRO > 1.0D-6) ;
  459.      'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
  460.      'ERREUR' 21 ;
  461.  'FINSI' ;
  462.  ERRO = 'MAXIMUM' (T2 '-' ('REDU' TN2 (TDOM2 . 'CENTRE'))) 'ABS' ;
  463.  ERRO = ERRO '/' (T2) ;
  464.  'SI' (ERRO > 1.0D-6) ;
  465.      'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
  466.      'ERREUR' 21 ;
  467.  'FINSI' ;
  468. **************************************************
  469. **** initial conditions **************************
  470. **************************************************
  471. *** density
  472. * DENSM = (CHAL1 * CHRN1) '+'  (CHAL2 * CHRN2) ; 
  473. * CHM_RN   =  'KCHA' $DOMTOT 'CHAM' DENSM ;
  474. * 'TRACER'  CHM_RN  $DOMTOT  
  475. *        'TITR'  ('CHAINE'  'RN at t='  0.0) ;
  476. **************************************************
  477. **** temperat
  478. * TEMPM = (CHAL1 * TN1) '+'  (CHAL2 * TN2) ;   
  479. * CHM_TN   =  'KCHA' $DOMTOT 'CHAM' TEMPM ;
  480. * 'TRACER'  CHM_TN  $DOMTOT  
  481. *        'TITR'  ('CHAINE'  'TN at t='  0.0) ;
  482.  
  483. * 'OPTION' DONN 5 ;   
  484. *
  485. **** Parameter for the time loop 
  486. *
  487. *    Names of the residuum components
  488. *
  489.  LISTINCO = ('MOTS' 'ALF1' 'RN1' 'RNX1' 'RNY1' 'RET1'
  490.           'ALF2' 'RN2' 'RNX2' 'RNY2' 'RET2') ;
  491. *
  492. **** BC
  493. *
  494.  KONLIM = 'DIFF' (TDOMLIM . 'CENTRE') (TDOMLIM . 'CENTRE') ;
  495.  CHPVID CHMVID = 'KOPS' 'MATRIK' ;
  496.  RESLIM = 'COPIER' CHPVID ;
  497. ************************************************************
  498. *   mass of burnt gase
  499. ************************************************************
  500.    MBURB = (CHAL2 * CHRN2) '*' ('DOMA' $DOMTOT 'VOLUME') ;
  501.    FMBURB = 'MAXIMUM'  ('RESULT' MBURB) ;
  502.    TPSM1 = 0.0 ;
  503.    DENUNB = (CHAL1 * CHRN1) ;
  504.    LSFLAM = 'PROG' R ;
  505.   'SI' ('EGA' ('VALEUR' 'MODE') 'AXIS') ;
  506.      LSFLAM = 'PROG'  ((R '**' 2.0) '/' 2.0D0) ;
  507.    'FINSI' ;
  508.  
  509. *   'LISTE'  FMBURB ;
  510. *   'LISTE'  ((CHAL2 * CHRN2)) ; 
  511. * 
  512. ****************************************************************
  513. ****************************************************************
  514. ****************************************************************
  515. *****           The computation                           ******
  516. ****************************************************************
  517. ****************************************************************
  518.  
  519.  AL1   = CHAL1          ;
  520.  AL2   = CHAL2          ;
  521.  ARN1  = CHAL1 * CHRN1  ;
  522.  ARN2  = CHAL2 * CHRN2  ;
  523.  AGN1  = CHAL1 * CHGN1  ;
  524.  AGN2  = CHAL2 * CHGN2  ;
  525.  AREN1 = CHAL1 * CHRET1 ;
  526.  AREN2 = CHAL2 * CHRET2 ;
  527. *
  528.  AL10   =  'COPIER' AL1   ;
  529.  AL20   =  'COPIER' AL2   ;
  530.  ARN10  =  'COPIER' ARN1  ;
  531.  ARN20  =  'COPIER' ARN2  ;
  532.  AGN10  =  'COPIER' AGN1  ;
  533.  AGN20  =  'COPIER' AGN2  ;
  534.  AREN10 =  'COPIER' AREN1 ;
  535.  AREN20 =  'COPIER' AREN2 ;
  536.  TN10   =  'COPIER' TN1   ;
  537.  TN20   =  'COPIER' TN2   ;
  538.  
  539. *
  540. **** Geometrical coefficient to compute grad(alpha)/|grad(alpha)|
  541. *
  542.  
  543.  EPSGR = 1.0D-6 ;
  544.  CHPL1 = CHPVID ;
  545.  CHPL2 = 'MANUEL' 'CHPO' (TDOMLIM . 'CENTRE') 2
  546.          'P1DX' 0.0 'P1DY'  0.0 ;
  547.  GRALP1 GRAL = 'PENT' $DOMTOT 'FACE'
  548.      'DIAMAN2'  ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY')
  549.      AL1  CHPL1 CHPL2 ;
  550. * 'SI' VRAI ;
  551. *    V1 = 'VECTEUR'
  552. *       ('NOMC' ('MOTS' 'P1DX' 'P1DY') GRALP1 ('MOTS' 'UX' 'UY')) ;
  553. *    'TRACER' DOMTOT V1 DOMLIM ;
  554. * 'FINSI' ;
  555.  GRALP1 = GRALP1 '+' EPSGR ;
  556.  GRALP1 = GRALP1 '/' (('PSCAL' GRALP1 GRALP1 ('MOTS' 'P1DX' 'P1DY')
  557.     ('MOTS' 'P1DX' 'P1DY')) '**' 0.5) ;
  558. * 'SI' VRAI ;
  559. *    V1 = 'VECTEUR'
  560. *       ('NOMC' ('MOTS' 'P1DX' 'P1DY') GRALP1 ('MOTS' 'UX' 'UY')) ;
  561. *    'TRACER' DOMTOT V1 DOMLIM ;
  562. * 'FINSI' ;
  563. *
  564. **** Geometrical coefficient to compute gradients
  565. *
  566.  GRADR0 ALR0 COEFSCAL =  'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
  567.           ('MOTS' 'SCAL') AL10 ;
  568.  GRADR0 = GRADR0 * 0.0 ;
  569.  ALR0   = ALR0   * 0.0 ;
  570.  
  571.  GRADV0 ALV0 COEFVECT =  'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
  572.           ('MOTS' 'UX' 'UY') AGN10 ;
  573.  GRADV0 = GRADV0 * 0.0 ;
  574.  ALV0   = ALV0   * 0.0 ;
  575.  
  576.  
  577. *
  578. **** Temporal loop
  579. *
  580.  
  581.  TPS = 0.0 ;
  582.  
  583.  'MESSAGE' ;
  584.  'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
  585.  'MESSAGE' ('CHAINE' 'SAFFAC      = ' SAFFAC) ;
  586.  'MESSAGE' ;
  587.  
  588.  'TEMPS' 'ZERO' ;
  589.  'REPETER' BLITER NITER ;
  590.  
  591. *
  592. **** Primitive variables
  593. *
  594.  
  595.     RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1 AL2
  596.        ARN1 ARN2 AGN1 AGN2 AREN1 AREN2 TN1M TN2M EPS ;
  597.  
  598.     TN1M = 'COPIER' TN1 ;
  599.     TN2M = 'COPIER' TN2 ;
  600.  
  601. **** Gradient of alpha
  602.  
  603.     GRALP1 = 'PENT' $DOMTOT 'FACE'
  604.     'DIAMAN2'  ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY')
  605.         AL1  CHPL1 CHPL2 'GRADGEO' GRAL  ; 
  606. *    'SI' VRAI ;
  607. *       V1 = 'VECTEUR'
  608. *          ('NOMC' ('MOTS' 'P1DX' 'P1DY') GRALP1 ('MOTS' 'UX' 'UY')) ;
  609. *       'TRACER' DOMTOT V1 DOMLIM ;
  610. *    'FINSI' ;
  611.     GRALP1 = GRALP1 '+' EPSGR ;
  612.     GRALP1 = GRALP1 '/' (('PSCAL' GRALP1 GRALP1 ('MOTS' 'P1DX' 'P1DY')
  613.        ('MOTS' 'P1DX' 'P1DY')) '**' 0.5) ;
  614.  
  615.     'SI' LOGSO ;
  616.  
  617.        GRADA1 ALA1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  618.           ('MOTS' 'SCAL') AL1  'GRADGEO' COEFSCAL ;
  619.  
  620.        GRADR1 ALR1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  621.           ('MOTS' 'SCAL') RN1  'GRADGEO' COEFSCAL ;
  622.  
  623.        GRADP1 ALP1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  624.           ('MOTS' 'SCAL') PN1  'GRADGEO' COEFSCAL ;
  625.  
  626.        GRADV1 ALV1 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
  627.           ('MOTS' 'UX' 'UY') VN1   'GRADGEO'  COEFVECT ;
  628.  
  629.        GRADA2 ALA2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  630.           ('MOTS' 'SCAL') AL2  'GRADGEO' COEFSCAL ;
  631.  
  632.        GRADR2 ALR2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  633.           ('MOTS' 'SCAL') RN2  'GRADGEO' COEFSCAL ;
  634.  
  635.        GRADP2 ALP2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  636.           ('MOTS' 'SCAL') PN2  'GRADGEO' COEFSCAL ;
  637.  
  638.        GRADV2 ALV2 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
  639.           ('MOTS' 'UX' 'UY') VN2   'GRADGEO'  COEFVECT ;
  640.  
  641.        CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
  642.          'PRET' 'DEM'  $DOMTOT
  643.           AL1  GRADA1  ALA1
  644.           AL2  GRADA2  ALA2
  645. *          AL1  GRADA1  ALR0
  646. *          AL2  GRADA2  ALR0
  647.           RN1  GRADR1  ALR1
  648.           RN2  GRADR2  ALR2 
  649.           VN1  GRADV1  ALV1
  650.           VN2  GRADV2  ALV2
  651.           PN1  GRADP1  ALP1
  652.           PN2  GRADP2  ALP2 ;
  653.  
  654.     'SINON' ;
  655.  
  656.        CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
  657.          'PRET' 'DEM'  $DOMTOT
  658.           AL1  GRADR0  ALR0
  659.           AL2  GRADR0  ALR0
  660.           RN1  GRADR0  ALR0
  661.           RN2  GRADR0  ALR0
  662.           VN1  GRADV0  ALV0
  663.           VN2  GRADV0  ALV0
  664.           PN1  GRADR0  ALR0
  665.           PN2  GRADR0  ALR0 ;
  666.  
  667.     'FINSI' ;
  668.  
  669.     RESIDU DELTAT SURFF = 'KONV' 'VF' 'DEM' 'RESI' METO 'CONS'
  670.           $DOMTOT PGAS LISTINCO AL1 AL2 CHFAL1 CHFAL2 CHFRN1 CHFRN2
  671.           CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS KONLIM ;
  672.  
  673.  
  674.     RESIDU = RESIDU '+' RESLIM ;
  675. *********************************************************
  676. *** axisymmetric case                                ****
  677. *   Contribution of pressure                         ****
  678. *   Warning: at first order, ALP initialised to 0    ****
  679. *********************************************************
  680.   'SI' ('EGA' ('VALEUR' 'MODE') 'AXIS') ;
  681. ******** first order for while  
  682.        ALP1 = 0.0 '*' ALP1 ;
  683.        ALP2 = 0.0 '*' ALP2 ;
  684.     RESIP1 = 'FIMP' 'VF' 'AXIS' 'RESI' $DOMTOT
  685.        ('MOTS' 'RN1' 'RNX1' 'RNY1' 'RET1')
  686.        (AL1 '*' PN1) GRADP1 ALP1 ;
  687.     RESIP2 = 'FIMP' 'VF' 'AXIS' 'RESI' $DOMTOT
  688.        ('MOTS' 'RN2' 'RNX2' 'RNY2' 'RET2')
  689.        (AL2 '*' PN2) GRADP2 ALP2 ;       
  690.     RESIDU  = RESIDU 'ET' RESIP1 'ET' RESIP2 ;
  691.   'FINSI' ;    
  692. ********************************************************
  693. ********************************************************    
  694. *    'REPETER' BL1 ('DIME' LISTINCO) ;  
  695. *       mot1 = 'EXTRAIRE' LISTINCO &BL1 ;
  696. *       valre1 = 'MAXIMUM' ('EXCO' RESLIM mot1) 'ABS' ;
  697. *       tit1 = 'CHAINE' 'Component ' mot1 ', reference value ' valre1 ;
  698. *       evaa = 'EVOL' 'CHPO' RESIDU mot1  LIGEVO ;
  699. *       'DESSIN' evaa 'TITRE' tit1 ;
  700. *    'FIN' BL1 ;
  701.  
  702. *    'OPTION' DONN 5 ;
  703. *   Problem with RNY1, RNY2
  704.  
  705.     DT_CON = SAFFAC '*' DELTAT ;     
  706. *
  707. **** The time step linked to tps
  708. *
  709.  
  710.     DTTPS = (TFINAL '-' TPS) * (1. '+' 1.0D-6) ;
  711.  
  712. *
  713. **** Total time step
  714. *
  715.  
  716.     DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ;
  717.  
  718. *
  719. **** Increment of the variables (predictor)
  720. *
  721. *    RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1 AL2
  722. *       ARN1 ARN2 AGN1 AGN2 AREN1 AREN2 CHTN1 CHTN2 EPS ;
  723.  
  724.     RESIDU = (0.5 * DTMIN) '*' RESIDU ;
  725.  
  726.     DALP1  = 'EXCO' 'ALF1' RESIDU 'SCAL' ;
  727.     DRN1   = 'EXCO' 'RN1'  RESIDU 'SCAL' ;
  728.     DGN1   = 'EXCO' ('MOTS' 'RNX1' 'RNY1') RESIDU ('MOTS' 'UX' 'UY') ;
  729.     DRET1  = 'EXCO' 'RET1' RESIDU 'SCAL' ;
  730.     DALP2  = 'EXCO' 'ALF2' RESIDU 'SCAL' ;
  731.     DRN2   = 'EXCO' 'RN2'  RESIDU 'SCAL' ;
  732.     DGN2   = 'EXCO' ('MOTS' 'RNX2' 'RNY2') RESIDU ('MOTS' 'UX' 'UY') ;
  733.     DRET2  = 'EXCO' 'RET2' RESIDU 'SCAL' ;
  734.  
  735.     AL1B   = AL1   '+'  DALP1 ;
  736.     ARN1B  = ARN1  '+'  DRN1  ;
  737.     AGN1B  = AGN1  '+'  DGN1  ;
  738.     AREN1B = AREN1 '+'  DRET1 ;
  739.     AL2B   = AL2   '+'  DALP2 ;
  740.     ARN2B  = ARN2  '+'  DRN2  ;
  741.     AGN2B  = AGN2  '+'  DGN2  ;
  742.     AREN2B = AREN2 '+'  DRET2 ;
  743.  
  744. *    
  745. **** Corrector
  746. *
  747.  
  748.     RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1B AL2B
  749.        ARN1B ARN2B AGN1B AGN2B AREN1B AREN2B TN1M TN2M EPS ;
  750.  
  751. **** Gradient of alpha
  752.  
  753.     GRALP1 = 'PENT' $DOMTOT 'FACE'
  754.     'DIAMAN2'  ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY')
  755.         AL1B  CHPL1 CHPL2 'GRADGEO' GRAL  ; 
  756. *    'SI' VRAI ;
  757. *       V1 = 'VECTEUR'
  758. *          ('NOMC' ('MOTS' 'P1DX' 'P1DY') GRALP1 ('MOTS' 'UX' 'UY')) ;
  759. *       'TRACER' DOMTOT V1 DOMLIM ;
  760. *    'FINSI' ;
  761.     GRALP1 = GRALP1 '+' EPSGR ;
  762.     GRALP1 = GRALP1 '/' (('PSCAL' GRALP1 GRALP1 ('MOTS' 'P1DX' 'P1DY')
  763.        ('MOTS' 'P1DX' 'P1DY')) '**' 0.5) ;
  764.  
  765.     'SI' LOGSO ;
  766.  
  767.        GRADA1 ALA1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  768.           ('MOTS' 'SCAL') AL1B  'GRADGEO' COEFSCAL ;
  769.  
  770.        GRADR1 ALR1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  771.           ('MOTS' 'SCAL') RN1  'GRADGEO' COEFSCAL ;
  772.  
  773.        GRADP1 ALP1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  774.           ('MOTS' 'SCAL') PN1  'GRADGEO' COEFSCAL ;
  775.  
  776.        GRADV1 ALV1 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
  777.           ('MOTS' 'UX' 'UY') VN1   'GRADGEO'  COEFVECT ;
  778.  
  779.        GRADA2 ALA2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  780.           ('MOTS' 'SCAL') AL2B  'GRADGEO' COEFSCAL ;
  781.  
  782.        GRADR2 ALR2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  783.           ('MOTS' 'SCAL') RN2  'GRADGEO' COEFSCAL ;
  784.  
  785.        GRADP2 ALP2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  786.           ('MOTS' 'SCAL') PN2  'GRADGEO' COEFSCAL ;
  787.  
  788.        GRADV2 ALV2 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
  789.           ('MOTS' 'UX' 'UY') VN2   'GRADGEO'  COEFVECT ;
  790.  
  791.  
  792.  
  793.        CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
  794.          'PRET' 'DEM'  $DOMTOT
  795.           AL1B  GRADA1  ALA1
  796.           AL2B  GRADA2  ALA2
  797.           RN1  GRADR1  ALR1
  798.           RN2  GRADR2  ALR2 
  799.           VN1  GRADV1  ALV1
  800.           VN2  GRADV2  ALV2
  801.           PN1  GRADP1  ALP1
  802.           PN2  GRADP2  ALP2 ;
  803.  
  804.     'SINON' ;
  805.  
  806.        CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
  807.          'PRET' 'DEM'  $DOMTOT
  808.           AL1B  GRADR0  ALR0
  809.           AL2B  GRADR0  ALR0
  810.           RN1  GRADR0  ALR0
  811.           RN2  GRADR0  ALR0
  812.           VN1  GRADV0  ALV0
  813.           VN2  GRADV0  ALV0
  814.           PN1  GRADR0  ALR0
  815.           PN2  GRADR0  ALR0 ;
  816.  
  817.     'FINSI' ;
  818.  
  819.     RESIDU DELTAT SURFP = 'KONV' 'VF' 'DEM' 'RESI' METO 'CONS'
  820.           $DOMTOT PGAS LISTINCO AL1B AL2B CHFAL1 CHFAL2 CHFRN1 CHFRN2
  821.           CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS KONLIM ;
  822.  
  823.     RESIDU = RESIDU '+' RESLIM ;
  824. *********************************************************
  825. *** axisymmetric case                                ****
  826. *   Contribution of pressure                         ****
  827. *   Warning: at first order, ALP initialised to 0    ****
  828. *********************************************************
  829.   'SI' ('EGA' ('VALEUR' 'MODE') 'AXIS') ;
  830. ******** first order for while  
  831.        ALP1 = 0.0 '*' ALP1 ;
  832.        ALP2 = 0.0 '*' ALP2 ;
  833.     RESIP1 = 'FIMP' 'VF' 'AXIS' 'RESI' $DOMTOT
  834.        ('MOTS' 'RN1' 'RNX1' 'RNY1' 'RET1')
  835.        (AL1B '*' PN1) GRADP1 ALP1 ;
  836.     RESIP2 = 'FIMP' 'VF' 'AXIS' 'RESI' $DOMTOT
  837.        ('MOTS' 'RN2' 'RNX2' 'RNY2' 'RET2')
  838.        (AL2B '*' PN2) GRADP2 ALP2 ;       
  839.     RESIDU  = RESIDU 'ET' RESIP1 'ET' RESIP2 ;
  840.   'FINSI' ;      
  841. *    
  842.  
  843.  
  844.     RESIDU = DTMIN '*' RESIDU ;
  845.  
  846.     DALP1  = 'EXCO' 'ALF1' RESIDU 'SCAL' ;
  847.     DRN1   = 'EXCO' 'RN1'  RESIDU 'SCAL' ;
  848.     DGN1   = 'EXCO' ('MOTS' 'RNX1' 'RNY1') RESIDU ('MOTS' 'UX' 'UY') ;
  849.     DRET1  = 'EXCO' 'RET1' RESIDU 'SCAL' ;
  850.     DALP2  = 'EXCO' 'ALF2' RESIDU 'SCAL' ;
  851.     DRN2   = 'EXCO' 'RN2'  RESIDU 'SCAL' ;
  852.     DGN2   = 'EXCO' ('MOTS' 'RNX2' 'RNY2') RESIDU ('MOTS' 'UX' 'UY') ;
  853.     DRET2  = 'EXCO' 'RET2' RESIDU 'SCAL' ;
  854.  
  855.     TPS   = TPS   '+'  DTMIN ;
  856.     AL1   = AL1   '+'  DALP1 ;
  857.     ARN1  = ARN1  '+'  DRN1  ;
  858.     AGN1  = AGN1  '+'  DGN1  ;
  859.     AREN1 = AREN1 '+'  DRET1 ;
  860.     AL2   = AL2   '+'  DALP2 ;
  861.     ARN2  = ARN2  '+'  DRN2  ;
  862.     AGN2  = AGN2  '+'  DGN2  ;
  863.     AREN2 = AREN2 '+'  DRET2 ;
  864.  
  865.     'SI' (((&BLITER '/' 10) '*' 10) 'EGA' &BLITER) ;
  866.         'MESSAGE' ('CHAINE' 'ITER =' &BLITER '  TPS =' TPS) ;
  867. *!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  868. *!!!!!    HERE                   !!!!!!!
  869. *!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  870.          IPOST = IPOST '+' 1 ;
  871. *         'LISTE'  DRN2 ;
  872. ************************************************
  873. ******* pressure transducers *******************
  874. ************************************************
  875.           PN = (AL1 '*' PN1) '+' (AL2 '*' PN2) ;
  876.           PPC = 'REDU' PN  PCTOT ;
  877.           PCEN.'TIME' = PCEN.'TIME' 'ET' ('PROG' TPS) ;
  878.         'REPETER' BOUC NPP ;
  879.            PP1 = 'EXTRAIRE' PPC 'SCAL' ('POIN' PCTOT &BOUC);
  880.           PCEN.'PRESSION'. &BOUC = (PCEN.'PRESSION'. &BOUC) 'ET'
  881.                                            ('PROG' PP1)  ;
  882.         'FIN' BOUC ;
  883. *************************************************
  884. ***** velocity along centerline *****************
  885. *************************************************
  886.          VN = (AL1 '*' VN1) '+' (AL2 '*' VN2) ;
  887.          VYV = 'EXCO' 'UY' VN ;
  888.          EVVYV  = ('EVOL' 'CHPO'  VYV LIGEVO) 'COULEUR' 'ROUG' ;
  889. *         'DESSIN' (EVVYV) 'GRILL' 'NCLK' ;
  890.          VCEN . IPOST = ('EVOL' 'CHPO'  VYV LIGEVO) ;
  891. *         'LISTE'  (VCEN .  IPOST) ;        
  892. *******************************************************                  
  893. ****** making variable ksi = \alpha_1 * (1-\alpha_1)
  894. *******************************************************
  895.           EVALP = 'EVOL' 'CHPO'  AL1 LIGEVO ;
  896. *          'DESSIN'  EVALP 'NCLK' ;
  897.           LCOOR = 'EXTRAIRE' EVALP 'ABSC' ; 
  898.           LALP1 = 'EXTRAIRE' EVALP 'ORDO' ;           
  899.           LUNI = 'PROG' ('DIME' LALP1) * 1.0 ; 
  900.           VKSI = (-1.0D0 '*' LALP1) '+' LUNI ;
  901.           VKSI = LALP1 '*' VKSI ;
  902.           EVKSI = 'EVOL' 'MANUEL' LCOOR VKSI ;
  903. ************************************************
  904. ****** new way of computing flame location *****
  905. ************************************************
  906.           EVKIS = 'EVOL' 'MANUEL' VKSI LCOOR ;
  907.           EVKI = 'EXTRAIRE' EVKIS 'APRE' 0.23D0 ;
  908.           EVK = 'EXTRAIRE'  EVKI 'ORDO' ;
  909.         'SI' (('DIME' EVK) '>' 1) ;  
  910.           X1S = 'EXTRAIRE'  EVK 1 ;
  911.           X2S = 'EXTRAIRE'  EVK ('DIME' EVK) ;
  912.           DISTFL1 = (X1S '+' X2S) '/' 2.0D0 ;
  913.           DFLAM1 = DFLAM1 'ET' ('PROG' DISTFL1) ;
  914.           TFLAM1 = TFLAM1 'ET' ('PROG' TPS) ;
  915.         'FINSI' ;
  916. **************************************************         
  917. ********* the flame is where the maximum of ksi
  918. **************************************************
  919.           OBJET2 DISTFL OBJET4  =  'MAXIMUM' EVKSI ;
  920.           DFLAM = DFLAM 'ET' ('PROG' DISTFL) ;
  921.           TFLAM = TFLAM 'ET' ('PROG' TPS) ;
  922.           EVDFLAM = 'EVOL' 'MANUEL' TFLAM DFLAM ;
  923. *********************************************************
  924. ********* Flame surface storage     *********************
  925. *********************************************************
  926.           LSFLAM = LSFLAM 'ET' ('PROG' SURFF) ;
  927.           EVSFLAM = 'EVOL' 'MANUEL' TFLAM LSFLAM ;
  928. *          'DESSIN'  EVSFLAM 'GRILL' 'NCLK' ;
  929. *-------------------------------------
  930. *      CHM_TN   =  'KCHA' $DOMTOT 'CHAM' AL1 ;
  931. *      CNT1 = 'CONTOUR' ('DOMA' $DOMTOT 'MAILLAGE') ;
  932. *     'TRACER'  CHM_TN  $DOMTOT CNT1
  933. *        'TITR'  ('CHAINE'  'TN at t='  TPS) 'NCLK' ;
  934. *    'FINSI' ;    
  935. *!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!         
  936. *!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!        
  937.     'FINSI' ;
  938.  
  939.     'SI' (TPS > TFINAL) ;
  940.        'QUITTER' BLITER ;
  941.     'FINSI' ;
  942.  
  943.  'FIN' BLITER ;
  944.  
  945. ****************************************************
  946. *****   Test de non-regression    ******************
  947. ****************************************************
  948. ***** exact flame surface
  949.   SEXACT = (R '**' 2.0) '/' 2.0D0 ;
  950. ***** maximum of computed flame surface  
  951.   SMCOM = 'MAXIMUM' LSFLAM ;
  952. *  'LISTE' SMCOM ;
  953. ****** relative error
  954.   ERREL = (SEXACT '-' SMCOM) '/' SEXACT ;
  955.   ERREL = ('ABS' ERREL) '*' 100.0D0 ;
  956.  
  957. *  'LISTE' ERREL ;
  958.  
  959.   'SI' (ERREL '>' 0.02 ) ;
  960.         'MESSAGE' 'Problem 1';
  961.         'ERREUR' 5 ;
  962.    'FINSI' ;
  963.  
  964.  
  965.  'FIN' ;
  966.  
  967.  
  968.  
  969.  
  970.  
  971.  
  972.  
  973.  
  974.  
  975.  
  976.  
  977.  
  978.  
  979.  
  980.  
  981.  
  982.  
  983.  
  984.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales