Télécharger contactd_fmm.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : contactd_fmm.dgibi
  2. *********************************************************************
  3. * *
  4. * Propagatrion d'une discontinuité de contact. *
  5. * *
  6. * Methode implicite sans matrice *
  7. * *
  8. * BECCANTINI A., SFME/LTMF, Fevrier 2003 *
  9. * *
  10. * Boundary conditions imposed on the border without using ghost *
  11. * cells. *
  12. * Boundary conditions correctly taken into account into the *
  13. * explicit part. *
  14. * They are not recomputed at each Jacobi iteration. *
  15. * Dual time step strategy. *
  16. * Automatic time step strategy *
  17. * BDF2 *
  18. * *
  19. * Juin 2006 : - changement des conditions aux limites pour NS *
  20. * (elle peuvent venir de la partie convective) *
  21. * - possibilité de utiliser SGS *
  22. * *
  23. *********************************************************************
  24.  
  25.  
  26. *************************************************************************
  27. *************************************************************************
  28. * I) PROCEDURES *********************************************************
  29. *************************************************************************
  30. *************************************************************************
  31.  
  32. *************************************************************************
  33. ******PROCEDURE PROCPT **************************************************
  34. *************************************************************************
  35. *
  36. *
  37. 'DEBPROC' PROCPT ;
  38. 'FINPROC' ;
  39.  
  40. *************************************************************************
  41. *********** PROCEDURE PNSSM : ********************
  42. *********** SOLUTION OF THE NAVIER-STOKES EQUATIONS ********************
  43. *********** FMM ********************
  44. *************************************************************************
  45. *
  46. * RVX . 'RESULTS' : table containing the results
  47. *
  48. * RVX . 'FREQI' : frequency of presenting the results
  49. *
  50. * RVX . 'MODEL' : model object
  51. * RVX . 'RN0' : density
  52. * RVX . 'GN0' : qdm
  53. * RVX . 'RET0' : total energy
  54. *
  55. * RVX . 'PGAS' : table containing the gas model
  56. * RVX . 'PGAS' . 'GAMN' : gamma
  57. * RVX . 'PGAS' . 'MU' : dynamic viscosity (kg/m^3 x m^2/s)
  58. * RVX . 'PGAS' . 'LAMBDA' : heat diffusion (W/K/m)
  59. *
  60. *
  61. * RVX . 'GRAVITY' : gravity
  62. *
  63. * RVX . 'LISTCONS' : name of the conservative variables
  64. * RVX . 'LISTPRIM' : name of the primitive variables
  65. *
  66. * RVX . 'LISTERR' : name of the error variables
  67. *
  68. * RVX . 'METHOD' : numerical scheme
  69. * RVX . 'CUTOFF' : cut off speed
  70. *
  71. * RVX . 'SPACEA' : space accuracy
  72. * RVX . 'LIMITER' : limiter type
  73. * RVX . 'TIMEA' : time accuracy
  74. *
  75. * RVX . 'T0' : initial time
  76. * RVX . 'TFINAL' : final time
  77. * RVX . 'DTPS' ('CFL') : time step or CFL number
  78. *
  79. * RVX . 'DCFL' : CFL number for dual time
  80. *
  81. * Error criteria for dual time loop:
  82. * RVX . 'NDTITER' : number of iterations
  83. * RVX . 'RELERR' : Logical (Relative error or absolute error)
  84. * RVX . 'EPSDT' : error at which dual time iterations are stopped
  85. *
  86. * RVX . 'NJAC' : Jacobi iterations
  87. *
  88. * RVX . 'PROLIM' : table called by the procedure PROLIM,
  89. * procedure to compute boundary conditions
  90. *
  91. * RVX . 'DIFTIMP' : boundary condition on temperature
  92. * RVX . 'DIFGTIMP' : boundary condition on gradient of temperature
  93. * RVX . 'MDIFCT' : mesh in which we impose the temperature
  94. * arising from the convective BC
  95. * RVX . 'DIFVIMP' : boundary condition on speed
  96. * RVX . 'DIFGVIMP' : boundary condition on gradient of speed
  97. * RVX . 'MDIFCV' : mesh in which we impose the velocity
  98. * arising from the convective BC
  99. * RVX . 'DIFTAUI' : boundary condition on constraint tensor
  100. * RVX . 'DIFQIMP' : boundary condition on heat flux
  101. *
  102. *
  103.  
  104. 'DEBPROC' PNSSM ;
  105. 'ARGUMENT' RVX*TABLE ;
  106. *
  107. *
  108. 'SI' ('EXISTE' RVX 'RESULTS') ;
  109. 'MESS' 'Table RESULTS already exists' ;
  110. 'SINON' ;
  111. RVX . 'RESULTS' = 'TABLE' ;
  112. RVX . 'RESULTS' . 'TPS' = RVX . 'T0' ;
  113. RVX . 'RESULTS' . 'RN' = 'COPIER' (RVX . 'RN0') ;
  114. RVX . 'RESULTS' . 'GN' = 'COPIER' (RVX . 'GN0') ;
  115. RVX . 'RESULTS' . 'RET' = 'COPIER' (RVX . 'RET0') ;
  116. RVX . 'RESULTS' . 'LISTLINF' = 'PROG' ;
  117. RVX . 'RESULTS' . 'LISTITDT' = 'LECT' ;
  118. RVX . 'RESULTS' . 'LISTITER' = 'LECT' ;
  119. RVX . 'RESULTS' . 'NITER' = 0 ;
  120. 'FINSI' ;
  121.  
  122. MDOMINT = RVX . 'MODEL' ;
  123. NELT = 'NBEL' ('DOMA' MDOMINT 'CENTRE') ;
  124. *
  125. ***** Physical properties
  126. *
  127. GAMN = RVX . 'PGAS' . 'GAMN' ;
  128. GAMSCA1 = 'MAXIMUM' GAMN ;
  129. GAMSCA2 = 'MINIMUM' GAMN ;
  130. 'SI' ('EGA' GAMSCA2 GAMSCA1 0.0001) ;
  131. GAMSCAL = GAMSCA1 ;
  132. 'SINON' ;
  133. 'MESSAGE' ;
  134. 'MESSAGE' 'Gamma is not constant' ;
  135. 'ERREUR' 21 ;
  136. 'FINSI' ;
  137. *
  138. MU = (RVX . 'PGAS' . 'MU') ;
  139. LAMBDA = (RVX . 'PGAS' . 'LAMBDA') ;
  140. R = (RVX . 'PGAS' . 'R') ;
  141. CV = R '/' (GAMSCAL '-' 1.) ;
  142.  
  143. *
  144. LISTINCO = RVX . 'LISTCONS' ;
  145. LISTPRIM = RVX . 'LISTPRIM' ;
  146. LISTERR = RVX . 'LISTERR' ;
  147. * Names of the gradient of temperature
  148. LMGTEMP = 'MOTS' 'P1DX' 'P1DY' ;
  149. * Names of the gradient of speed
  150. LMGVIT = 'MOTS' 'P1DX' 'P1DY' 'P2DX' 'P2DY' ;
  151.  
  152. *
  153. * Upwind scheme
  154. *
  155. METO = RVX . 'METHOD' ;
  156.  
  157. * Space accuracy (1 or 2) and limiter
  158. * Time accuracy (1 or 2)
  159.  
  160. ORDESP = RVX . 'SPACEA' ;
  161. TYPELIM = RVX . 'LIMITER' ;
  162. ORDTPS = RVX . 'TIMEA' ;
  163.  
  164. * Initial/final time
  165. * Deltat or CFL
  166. TPS = RVX . 'RESULTS' . 'TPS' ;
  167. TFINAL = RVX . 'TFINAL' ;
  168. 'SI' ('EXISTE' RVX 'DTPS') ;
  169. 'SI' ('EXISTE' RVX 'CFL') ;
  170. 'MESSAGE' 'DTPS or CFL ???' ;
  171. 'ERREUR' 21 ;
  172. 'FINSI' ;
  173. DTPS = RVX . 'DTPS' ;
  174. 'SINON' ;
  175. CFL = RVX . 'CFL' ;
  176. 'FINSI' ;
  177.  
  178. * Dual time iterations
  179. NDT = RVX . 'NDTITER' ;
  180. * Relative error
  181. EPSDT = RVX . 'EPSDT' ;
  182. * Jacobi iterations
  183. * NJAC = RVX . 'NJAC' ;
  184. * Cut off speed
  185. ICO = RVX . 'CUTOFF' ;
  186. * To compute the diffusive cut-off
  187. DELTAX = 'DOMA' MDOMINT 'XXDIEMIN' ;
  188. USDELTAX = 'INVERSE' DELTAX ;
  189. *
  190. **** Conservative variables
  191. *
  192. *
  193. MOT1 = 'EXTRAIRE' LISTINCO 1 ;
  194. 'SI' ('EGA' ('VALE' 'DIME') 2) ;
  195. NOMMOM = 'EXTRAIRE' LISTINCO
  196. ('LECT' 2 3 ) ;
  197. NOMVEL = 'MOTS' 'UX' 'UY' ;
  198. MOT2 = 'EXTRAIRE' LISTINCO 4 ;
  199. 'SINON' ;
  200. NOMMOM = 'EXTRAIRE' LISTINCO
  201. ('LECT' 2 3 4) ;
  202. NOMVEL = 'MOTS' 'UX' 'UY' 'UZ' ;
  203. MOT2 = 'EXTRAIRE' LISTINCO 5 ;
  204. 'FINSI' ;
  205.  
  206. RN0 = 'REDU' (RVX . 'RESULTS' . 'RN') ('DOMA' MDOMINT 'CENTRE') ;
  207. GN0 = 'REDU' (RVX . 'RESULTS' . 'GN') ('DOMA' MDOMINT 'CENTRE') ;
  208. RET0 = 'REDU' (RVX . 'RESULTS' . 'RET') ('DOMA' MDOMINT 'CENTRE') ;
  209. *
  210. **** Primitive variables
  211. *
  212. VN0 PN0 = 'PRIM' 'PERFMONO' RN0 GN0 RET0 GAMN ;
  213. TN0 = PN0 '/' (R '*' RN0) ;
  214.  
  215. MOTRN = 'EXTRAIRE' LISTPRIM 1 ;
  216. 'SI' ('EGA' ('VALE' 'DIME') 2) ;
  217. MOTVN = 'EXTRAIRE' LISTPRIM
  218. ('LECT' 2 3 ) ;
  219. MOTPN = 'EXTRAIRE' LISTPRIM 4 ;
  220. 'SINON' ;
  221. MOTVN = 'EXTRAIRE' LISTPRIM
  222. ('LECT' 2 3 4) ;
  223. MOTPN = 'EXTRAIRE' LISTPRIM 5 ;
  224. 'FINSI' ;
  225. *
  226. ********************************************************************
  227. **** Coeff to compute gradients for convective term (MCHSCA, MCHVEC)
  228. **** and for diffusive terms (MCHDIT, MCHVEC)
  229. ********************************************************************
  230. *
  231. * Boundary conditions have to be taken into account
  232. *
  233.  
  234. RCHLIM RESLIM = PROLIM (RVX . 'PROLIM') MDOMINT LISTINCO LISTPRIM
  235. RN0 VN0 PN0 GAMN ;
  236.  
  237. MAILLIM = 'EXTRAIRE' RCHLIM 'MAILLAGE' ;
  238.  
  239. *
  240. * If only wall conditions (Maillim is then empty)
  241. *
  242.  
  243. 'SI' (MAILLIM 'EGA' 0) ;
  244.  
  245. CHPVID PIPI = 'KOPS' 'MATRIK' ;
  246. MAILLIM = 'DIFF' ('DOMA' MDOMINT 'CENTRE')
  247. ('DOMA' MDOMINT 'CENTRE') ;
  248.  
  249. * Convective geometric coefficients
  250.  
  251. GRADRN ALRN MCHSCA = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' 'NOLIMITE'
  252. ('MOTS' 'SCAL') RN0 ;
  253.  
  254. GRADVN ALVN MCHVEC = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' 'NOLIMITE'
  255. NOMVEL GN0 'CLIM' (RVX . 'DIFVIMP') ;
  256.  
  257. * Diffusive geometric coefficients
  258.  
  259. GRADTN MCHDIT = 'PENT' MDOMINT 'FACE' 'DIAMAN2' ('MOTS' 'SCAL')
  260. LMGTEMP TN0 (RVX . 'DIFTIMP') (RVX . 'DIFGTIMP') ;
  261.  
  262. GRADVN MCHDIV = 'PENT' MDOMINT 'FACE' 'DIAMAN2' NOMVEL LMGVIT
  263. VN0 (RVX . 'DIFVIMP') (RVX . 'DIFGVIMP') ;
  264.  
  265. *
  266. * If other imposed Boudary conditions we create false Chpoints
  267. *
  268.  
  269. 'SINON' ;
  270.  
  271. * Convective false champoints
  272.  
  273. SCALBC = 'MANUEL' 'CHPO' MAILLIM 1 'SCAL' 1.0
  274. 'NATU' 'DISCRET' ;
  275.  
  276. 'SI' (('VALEUR' 'DIME') 'EGA' 2) ;
  277. VECTBC = 'MANUEL' 'CHPO' MAILLIM
  278. 2 'UX' 0.0 'UY' 0.0 'NATU' 'DISCRET' ;
  279. * Same name as NOMVEL
  280. 'SINON' ;
  281. VECTBC = 'MANUEL' 'CHPO' MAILLIM
  282. 3 'UX' 0.0 'UY' 0.0 'UZ' 0.0 'NATU' 'DISCRET' ;
  283. 'FINSI' ;
  284.  
  285. VECTBC = VECTBC '+' (RVX . 'DIFVIMP') ;
  286.  
  287. * Convective geometric coefficients
  288.  
  289. GRADRN ALRN MCHSCA = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' 'NOLIMITE'
  290. ('MOTS' 'SCAL') RN0 'CLIM' SCALBC ;
  291.  
  292. GRADVN ALVN MCHVEC = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' 'NOLIMITE'
  293. NOMVEL VN0 'CLIM' VECTBC ;
  294. 'FINSI' ;
  295.  
  296. *
  297. * Diffusive false champoints
  298. *
  299. * We have several possibility
  300. * We decide to impose the temperature given by the convective boundary
  301. * condition on (RVX . 'MDIFCT') and the velocity given by the
  302. * convective boundary condition on (RVX . 'MDIFCV')
  303. *
  304. * We check that (RVX . 'MDIFCT') belongs to MAILLIM
  305. *
  306. NN1 = 'NBEL' MAILLIM ;
  307. NN2 = 'NBEL' (RVX . 'MDIFCT') ;
  308. 'SI' (('NEG' NN1 0) 'ET' ('NEG' NN2 0)) ;
  309. *
  310. *** If both meshes contain at least one element NN1 is equal
  311. * to the number of elements belonging to the intersection.
  312. * Otherwise NN1 = 0
  313. *
  314. NN1 = 'NBEL' ('INTERSECTION' (RVX . 'MDIFCT') MAILLIM) ;
  315. 'FINSI' ;
  316. 'SI' ('NEG' NN1 NN2) ;
  317. 'MESSAGE' 'Problem in MDIFCT' ;
  318. 'ERREUR' 21 ;
  319. 'FINSI' ;
  320. SCALBC = 'MANUEL' 'CHPO' (RVX . 'MDIFCT') 1 'SCAL' 1.0
  321. 'NATU' 'DISCRET' ;
  322. SCALBC = SCALBC '+' (RVX . 'DIFTIMP') ;
  323.  
  324. NN2 = 'NBEL' (RVX . 'MDIFCV') ;
  325. 'SI' (('NEG' NN1 0) 'ET' ('NEG' NN2 0)) ;
  326. NN1 = 'NBEL' ('INTERSECTION' (RVX . 'MDIFCV') MAILLIM) ;
  327. 'FINSI' ;
  328. 'SI' ('NEG' NN1 NN2) ;
  329. 'MESSAGE' 'Problem in MDIFCV' ;
  330. 'ERREUR' 21 ;
  331. 'FINSI' ;
  332. 'SI' (('VALEUR' 'DIME') 'EGA' 2) ;
  333. VECTBC = 'MANUEL' 'CHPO' (RVX . 'MDIFCV')
  334. 2 'UX' 0.0 'UY' 0.0 'NATU' 'DISCRET' ;
  335. 'SINON' ;
  336. VECTBC = 'MANUEL' 'CHPO' (RVX . 'MDIFCV')
  337. 3 'UX' 0.0 'UY' 0.0 'UZ' 0.0 'NATU' 'DISCRET' ;
  338. 'FINSI' ;
  339. *
  340. VECTBC = VECTBC '+' (RVX . 'DIFVIMP') ;
  341. *
  342. * Diffusive geometric coefficients
  343.  
  344. GRADTN MCHDIT = 'PENT' MDOMINT 'FACE' 'DIAMAN2' ('MOTS' 'SCAL')
  345. LMGTEMP TN0 SCALBC (RVX . 'DIFGTIMP') ;
  346. *
  347. GRADVN MCHDIV = 'PENT' MDOMINT 'FACE' 'DIAMAN2' NOMVEL LMGVIT
  348. VN0 VECTBC (RVX . 'DIFGVIMP') ;
  349. *
  350. ***************************************************************
  351. *** After each dual time loop, we could display
  352. * the evolution of the error in the dual time loop
  353. * the evolution of the dual time step (the safety factor)
  354. ***************************************************************
  355. *
  356.  
  357. LISTLINF = RVX . 'RESULTS' . 'LISTLINF' ;
  358. LISTITDT = RVX . 'RESULTS' . 'LISTITDT' ;
  359. LISTITER = RVX . 'RESULTS' . 'LISTITER' ;
  360. *
  361. 'MESSAGE' ;
  362. 'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
  363. 'MESSAGE' ;
  364. *
  365. 'TEMPS' 'ZERO' ;
  366. *
  367. ************************************************************************
  368. ************************************************************************
  369. **** Temporal loop *****************************************************
  370. ************************************************************************
  371. ************************************************************************
  372.  
  373. RN_N1M1 = 'COPIER' RN0 ;
  374. GN_N1M1 = 'COPIER' GN0 ;
  375. RET_N1M1 = 'COPIER' RET0 ;
  376.  
  377. AA = 'DIME' LISTITER ;
  378. 'SI' (AA > 0) ;
  379. PTITER = 'EXTRAIRE' LISTITER AA ;
  380. 'SINON' ;
  381. PTITER = 0 ;
  382. 'FINSI' ;
  383.  
  384. DUSDT = 0.0D0 ;
  385.  
  386. 'REPETER' BLITER (RVX . 'NITER') ;
  387.  
  388. PTITER = PTITER '+' 1 ;
  389. *
  390. **** Personal procedure
  391. *
  392. PROCPT RVX ;
  393. *
  394. *
  395. **** _N1M = (t^n,\tau^m)
  396. * _N1M1 = (t^n,\tau^{m+1})
  397. *
  398. *
  399. ************************************************************************
  400. ****** Loop on dual time***********************************************
  401. ************************************************************************
  402. *
  403. *
  404. *** DUSDT0 is the increment of DUSDT in the previous (physical) time
  405. * iteration.
  406. *
  407. DUSDT0 = DUSDT ;
  408. DUSDT = 0.0D0 ;
  409.  
  410. 'REPETER' BLDT NDT ;
  411.  
  412. RN_N1M = RN_N1M1 ;
  413. GN_N1M = GN_N1M1 ;
  414. RET_N1M = RET_N1M1 ;
  415. *
  416. **** Primitive variables
  417. *
  418.  
  419. VN_N1M PN_N1M = 'PRIM' 'PERFMONO' RN_N1M GN_N1M RET_N1M GAMN ;
  420. *
  421. *** Boundary conditions
  422. *
  423.  
  424. RCHLIM RESLIM = PROLIM (RVX . 'PROLIM') MDOMINT LISTINCO LISTPRIM
  425. RN_N1M VN_N1M PN_N1M GAMN ;
  426.  
  427. *
  428. ****** First/second order reconstruction
  429. *
  430. 'SI' (ORDESP 'EGA' 2) ;
  431. *
  432. NNLIM = 'NBNO' MAILLIM ;
  433. 'SI' (NNLIM 'EGA' 0) ;
  434. RNLIM = CHPVID ;
  435. PNLIM = CHPVID ;
  436. VNLIM = RVX . 'DIFVIMP' ;
  437. 'SINON' ;
  438. RNLIM = ('EXCO' MOTRN RCHLIM 'SCAL') ;
  439. PNLIM = ('EXCO' MOTPN RCHLIM 'SCAL') ;
  440. VNLIM = ('EXCO' MOTVN RCHLIM MOTVN)
  441. '+' (RVX . 'DIFVIMP') ;
  442. 'FINSI' ;
  443. GRADRN ALRN0 = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' TYPELIM
  444. ('MOTS' 'SCAL') RN_N1M
  445. 'CLIM' RNLIM 'GRADGEO' MCHSCA ;
  446. GRADPN ALPN0 = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' TYPELIM
  447. ('MOTS' 'SCAL') PN_N1M
  448. 'CLIM' PNLIM 'GRADGEO' MCHSCA ;
  449. GRADVN ALVN0 = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' TYPELIM
  450. ('MOTS' 'UX' 'UY') VN_N1M
  451. 'CLIM' VNLIM 'GRADGEO' MCHVEC ;
  452.  
  453. * 'SI' (&BLDT < NLCB) ;
  454. * ALRN0 = 'COPIER' ALRN ;
  455. * ALPN0 = 'COPIER' ALPN ;
  456. * ALVN0 = 'COPIER' ALVN ;
  457. * 'SINON' ;
  458. * 'SI' (&BLDT 'EGA' NLCB) ;
  459. * 'MESSAGE' ;
  460. * 'MESSAGE' 'On gele les limiteurs!!!' ;
  461. * 'MESSAGE' ;
  462. * 'FINSI' ;
  463. * 'FINSI' ;
  464.  
  465. ROF VITF PF GAMF = 'PRET' 'PERFMONO' 2 1
  466. MDOMINT
  467. RN_N1M GRADRN ALRN0
  468. VN_N1M GRADVN ALVN0
  469. PN_N1M GRADPN ALPN0
  470. GAMN ;
  471. 'SINON' ;
  472. ROF VITF PF GAMF = 'PRET' 'PERFMONO' 1 1
  473. MDOMINT
  474. RN_N1M VN_N1M PN_N1M GAMN ;
  475. 'FINSI' ;
  476.  
  477. RESIDU DELTAT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
  478. MDOMINT LISTINCO ROF VITF PF GAMF MAILLIM
  479. * ICO (MU '*' ('INVERSE' RN_N1M) '*' USDELTAX) ;
  480. ICO ICO ;
  481. RESIDU = RESIDU '+' RESLIM ;
  482.  
  483. *
  484. **** La gravite
  485. *
  486.  
  487. RESGRA = 'FIMP' 'VF' 'GRAVMONO' 'RESI' LISTINCO
  488. RN_N1M GN_N1M1 (RVX . 'GRAVITY') ;
  489.  
  490. RESIDU = RESIDU '+' RESGRA ;
  491. *
  492. **************************************
  493. **** Diffusive terms *****************
  494. **************************************
  495. *
  496. TN_N1M = PN_N1M '/' (R '*' RN_N1M) ;
  497. *
  498. * Computation of the CHPOINTS at the boundary of the domain
  499. *
  500. NN1 = 'NBEL' MAILLIM ;
  501. 'SI' (NN1 > 0) ;
  502. RNLIM = 'EXCO' MOTRN RCHLIM 'SCAL' 'NATU' 'DISCRET' ;
  503. PNLIM = 'EXCO' MOTPN RCHLIM 'SCAL' 'NATU' 'DISCRET' ;
  504. VNLIM = ('EXCO' MOTVN RCHLIM MOTVN 'NATU' 'DISCRET') ;
  505. TNLIM = (PNLIM '/' (RNLIM '*' R) ) ;
  506. VNLIM = ('REDU' VNLIM (RVX . 'MDIFCV')) '+'
  507. (RVX . 'DIFVIMP') ;
  508. TNLIM = ('REDU' TNLIM (RVX . 'MDIFCV')) '+'
  509. (RVX . 'DIFTIMP') ;
  510. 'SINON' ;
  511. VNLIM = (RVX . 'DIFVIMP') ;
  512. TNLIM = (RVX . 'DIFVIMP') ;
  513. 'FINSI' ;
  514. *
  515. GRADTN = 'PENT' MDOMINT 'FACE' 'DIAMAN2' ('MOTS' 'SCAL')
  516. LMGTEMP TN_N1M TNLIM (RVX . 'DIFGTIMP')
  517. 'GRADGEO' MCHDIT ;
  518. *
  519. * 'LIST' (RVX . 'DIFGVIMP') ;
  520. GRADVN = 'PENT' MDOMINT 'FACE' 'DIAMAN2' NOMVEL
  521. LMGVIT VN_N1M VNLIM (RVX . 'DIFGVIMP')
  522. 'GRADGEO' MCHDIV ;
  523. *
  524. * NOMVEL = 'MOTS' 'UX' 'UY' ;
  525. *
  526. *
  527. ICACCA RESIDI DTCACCA = 'LAPN' 'VF' 'PROPCOST' 'RESI' 'EXPL'
  528. MDOMINT MU LAMBDA CV RN_N1M VN_N1M TN_N1M GRADVN GRADTN
  529. LISTINCO 'VIMP' VNLIM 'TAUI' (RVX . 'DIFTAUI')
  530. 'QIMP' (RVX . 'DIFQIMP') ;
  531.  
  532. RESIDU = RESIDU '+' RESIDI ;
  533. *
  534. ***** Spectral radious of the viscosity matrix
  535. *
  536. COEFV1 = (0.0 '*' RN_N1M) '+' 1.0 ;
  537. COEFV1 = COEFV1 '*' (lambda '*' (gamscal '-' 1.)) ;
  538. COEFV1 = COEFV1 '/' (R '*' RN_N1M) ;
  539. COEFV2 = (0.0 '*' RN_N1M) '+' ((4.0 '/' 3.0) '*' MU) ;
  540. COEFV2 = COEFV2 '/' RN_N1M ;
  541. COEFV = 0.5 '*' (COEFV1 '+' COEFV2) ;
  542. COEFV = COEFV '+' (0.5 '*' ((COEFV1 '-' COEFV2) 'ABS')) ;
  543. * 'MESSAGE' ;
  544. * 'MESSAGE' 'Je fais n importe' ;
  545. * COEFV = COEFV '*' 0.25 ;
  546. *
  547. ****** Residuum for dual tims stepping also involved the
  548. * variation of the conserved variables with respect
  549. * to time
  550.  
  551. 'SI' ((&BLITER 'EGA' 1) 'OU' (ORDTPS 'EGA' 1)) ;
  552. RESIDU = RESIDU '-' DUSDT ;
  553. 'SINON' ;
  554. RESIDU = RESIDU '-' ((1.5 '*' DUSDT) '-' (0.5 '*' DUSDT0)) ;
  555. 'FINSI' ;
  556.  
  557. *
  558. *** Time step at the first iteration/jacobi iteration
  559. *
  560.  
  561. 'SI' (&BLDT 'EGA' 1) ;
  562. 'SI' ('EXISTE' RVX 'CFL') ;
  563. DTPS = (RVX . 'CFL') '*' 2.0D0 '*' DELTAT ;
  564. 'SINON' ;
  565. DTPS = RVX . 'DTPS' ;
  566. 'FINSI' ;
  567. DTPS = 'MINIMUM' ('PROG' DTPS ((TFINAL '-' TPS) '*' 1.001)) ;
  568. TPS = TPS '+' DTPS ;
  569. *
  570. NJAC = 'ENTIER' ('MINIMUM' (RVX . 'NJACITER')) ;
  571. NJAC0 = NJAC ;
  572. 'SINON' ;
  573. NJAC = 'IPOL' (('LOG' ERRINF) '/' ('LOG' 10))
  574. (RVX . 'NJACLERR') (RVX . 'NJACITER') ;
  575. NJAC = 'MAXIMUM' ('PROG' NJAC0 NJAC) ;
  576. NJAC = 'ENTIER' NJAC ;
  577. NJAC0 = NJAC ;
  578. 'FINSI' ;
  579.  
  580. *
  581. *** JACOBI
  582. *
  583.  
  584. *
  585. **** CFL dual
  586. *
  587. 'SI' (&BLDT 'EGA' 1) ;
  588. SAFFACD = ('MINIMUM' (RVX . 'DCFL')) '*' 2 ;
  589. 'SINON' ;
  590. SAFFACD = ('IPOL' (('LOG' ERRINF) '/' ('LOG' 10))
  591. (RVX . 'DCFLERR') (RVX . 'DCFL')) '*' 2 ;
  592. 'FINSI' ;
  593.  
  594. 'SI' ((&BLITER 'EGA' 1) 'OU' (ORDTPS 'EGA' 1)) ;
  595. * DUN IPRO = 'DETO' (RVX . 'TYPEJAC')
  596. DUN IPRO = 'KONV' 'VF' 'PMON1FMM' (RVX . 'TYPEJAC')
  597. LISTINCO MDOMINT RESIDU RN_N1M GN_N1M RET_N1M GAMN ICO
  598. DTPS SAFFACD
  599. NJAC 'CLIM' LISTPRIM RCHLIM COEFV ;
  600. 'SINON' ;
  601. * DUN IPRO = 'DETO' (RVX . 'TYPEJAC')
  602. DUN IPRO = 'KONV' 'VF' 'PMON1FMM' (RVX . 'TYPEJAC')
  603. LISTINCO MDOMINT RESIDU RN_N1M GN_N1M RET_N1M GAMN ICO
  604. (DTPS '/' 1.5)
  605. SAFFACD NJAC 'CLIM' LISTPRIM RCHLIM COEFV ;
  606. 'FINSI' ;
  607. * 'FINSI' ;
  608.  
  609. 'SI' (IPRO 'NEG' 0) ;
  610. 'MESSAGE' ;
  611. 'MESSAGE' 'Probleme dans FMM' ;
  612. 'MESSAGE' ;
  613. 'ERREUR' 21 ;
  614. 'FINSI' ;
  615.  
  616. *
  617. **** We compute DUSDT for the future loop
  618. *
  619.  
  620. DUSDT = DUSDT '+' (DUN '/' DTPS) ;
  621.  
  622. *
  623. **** We evaluate the conservative variables at t^{n+1}, \tau^{m+1}
  624. *
  625.  
  626. DRN = 'EXCO' MOT1 DUN 'SCAL' ;
  627. DGN = 'EXCO' NOMMOM DUN NOMVEL ;
  628. DRET = 'EXCO' MOT2 DUN 'SCAL' ;
  629.  
  630. RN_N1M1 = RN_N1M '+' DRN ;
  631. GN_N1M1 = GN_N1M '+' DGN ;
  632. RET_N1M1 = RET_N1M '+' DRET ;
  633.  
  634. ERRINF = 'MAXIMUM' DUN 'ABS' LISTERR ;
  635. LISTLINF = LISTLINF 'ET' ('PROG' ERRINF) ;
  636. LISTITDT = LISTITDT 'ET' ('LECT' &BLDT) ;
  637. LISTITER = LISTITER 'ET' ('LECT' PTITER) ;
  638.  
  639. 'SI' ((&BLDT 'EGA' 1) 'OU'
  640. (((&BLDT '/' (RVX . 'FREQI')) '*' (RVX . 'FREQI'))
  641. 'EGA' &BLDT)) ;
  642. 'MESSAGE' ;
  643. 'MESSAGE'
  644. ('CHAINE' 'ITER =' PTITER ' TPS =' TPS
  645. ' DTITER =' &BLDT ' LINF =' ERRINF
  646. ' DCFL =' SAFFACD ' NJAC =' NJAC) ;
  647. 'MESSAGE' ;
  648. 'FINSI' ;
  649. *
  650. *
  651. *** Update of RVX . 'RESULTS'
  652. *
  653. RVX . 'RESULTS' . 'RN' = RN_N1M1 ;
  654. RVX . 'RESULTS' . 'GN' = GN_N1M1 ;
  655. RVX . 'RESULTS' . 'RET' = RET_N1M1 ;
  656. RVX . 'RESULTS' . 'LISTITDT' = LISTITDT ;
  657. RVX . 'RESULTS' . 'LISTITER' = LISTITER ;
  658. RVX . 'RESULTS' . 'LISTLINF' = LISTLINF ;
  659. *
  660. 'SI' (RVX . 'RELERR') ;
  661. * Relative error
  662. 'SI' (&BLDT 'EGA' 1) ;
  663. ERRINF0 = ERRINF ;
  664. 'SINON' ;
  665. 'SI' (ERRINF < (EPSDT '*' ERRINF0)) ;
  666. 'QUITTER' BLDT ;
  667. 'FINSI' ;
  668. 'FINSI' ;
  669. 'SINON' ;
  670. 'SI' (ERRINF < EPSDT) ;
  671. 'QUITTER' BLDT ;
  672. 'FINSI' ;
  673. 'FINSI' ;
  674.  
  675.  
  676. 'FIN' BLDT ;
  677.  
  678. ************************************************************************
  679. ****** End of the loop on dual time*************************************
  680. ************************************************************************
  681. *
  682. *** Update of RVX . 'RESULTS'
  683. *
  684. RVX . 'RESULTS' . 'TPS' = TPS ;
  685. RVX . 'RESULTS' . 'NITER' = (RVX . 'RESULTS' . 'NITER') '+' 1 ;
  686. *
  687. 'SI' (TPS '>EG' TFINAL) ;
  688. 'QUITTER' BLITER ;
  689. 'FINSI' ;
  690.  
  691. 'FIN' BLITER ;
  692.  
  693. TCPU = 'TEMPS' 'NOEC' ;
  694. RVX . 'RESULTS' . 'TCPU' = TCPU ;
  695.  
  696. 'FINPROC' ;
  697.  
  698. *************************************************************************
  699. ****** FIN PROCEDURE PNSSM **********************************************
  700. *************************************************************************
  701.  
  702. *************************************************************************
  703. ******PROCEDURE PROLIM **************************************************
  704. *************************************************************************
  705. *
  706. 'DEBPROC' PROLIM ;
  707.  
  708. 'ARGUMENT' RVX*'TABLE' MDOMINT*'MMODEL' LISTINCO*'LISTMOTS'
  709. LISTPRIM*'LISTMOTS'
  710. RN*'CHPOINT' VN*'CHPOINT' PN*'CHPOINT' GAMN*'CHPOINT' ;
  711.  
  712. *
  713. * Récupération du nombres de CL imposées
  714. *
  715.  
  716. NCL = RVX . 'N' ;
  717.  
  718. *
  719. * Initialisation de RESLIM et de RCHLIM
  720. *
  721.  
  722. RESLIM VIDE = 'KOPS' MATRIK ;
  723. RCHLIM VIDE = 'KOPS' MATRIK ;
  724.  
  725. *
  726. * Boucle d'évaluation des résidus et des CHPOINTS aux limites
  727. *
  728.  
  729. 'REPETER' BCLIM NCL ;
  730.  
  731. MOTCLI = 'EXTRAIRE' (RVX . 'CLN' . &BCLIM) 1 ;
  732. MODELI = RVX . 'MODELN' . &BCLIM ;
  733. CHPOLI = RVX . 'CHPOLN' . &BCLIM ;
  734.  
  735. RCHLI RESLI = 'KONV' 'VF' 'PERFMONO' 'CLIM' 'RESI'
  736. MDOMINT MODELI LISTINCO LISTPRIM
  737. RN VN PN GAMN
  738. CHPOLI MOTCLI ;
  739.  
  740. RCHLIM = RCHLIM '+' RCHLI ;
  741. RESLIM = RESLIM '+' RESLI ;
  742.  
  743. 'FIN' BCLIM ;
  744.  
  745. 'RESPRO' RCHLIM RESLIM ;
  746. 'FINPROC' ;
  747.  
  748. *************************************************************************
  749. ****** FIN PROCEDURE PROLIM *********************************************
  750. *************************************************************************
  751.  
  752. *************************************************************************
  753. ***** PROCEDURE POUR CALCULER LES VARIABLES CONSERVATIVES ***************
  754. *************************************************************************
  755.  
  756. 'DEBPROC' CONS ;
  757. 'ARGUMENT' RN*'CHPOINT' VN*'CHPOINT' PN*'CHPOINT' GAMN*'CHPOINT' ;
  758.  
  759. RVN = RN '*' VN ('MOTS' 'SCAL' 'SCAL') ('MOTS' 'UX' 'UY')
  760. ('MOTS' 'UX' 'UY') ;
  761. CELL = 'PSCAL' RVN VN ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
  762. RECIN = 0.5 '*' CELL ;
  763. REIN = PN '/' (GAMN '-' 1.0) ;
  764. RET = RECIN '+' REIN ;
  765.  
  766. DETR CELL ;
  767. DETR RECIN ;
  768. DETR REIN ;
  769. 'RESPRO' RVN RET ;
  770.  
  771. 'FINPROC' ;
  772.  
  773. *************************************************************************
  774. *************************************************************************
  775. ****************** FIN PROCEDURES **************************************
  776. *************************************************************************
  777. *************************************************************************
  778.  
  779.  
  780. *************************************************************************
  781. *************************************************************************
  782. ****************** MAILLAGE ********************************************
  783. *************************************************************************
  784. *************************************************************************
  785.  
  786. TYEL = 'QUA4' ;
  787.  
  788. 'OPTION' 'DIME' 2 'ELEM' TYEL 'ISOV' 'SULI'
  789. 'ECHO' 1 'TRAC' 'X' ;
  790.  
  791. * GRAPH = VRAI ;
  792. GRAPH = FAUX ;
  793. *
  794.  
  795. RAF = 1 ;
  796. L = 1. ;
  797.  
  798. NX = 50 '*' RAF ;
  799. DX = L '/' NX ;
  800. * NY = 4 ;
  801. 'MESSAGE' ;
  802. 'MESSAGE' 'MODIF' ;
  803. 'MESSAGE' ;
  804. NY = 1 ;
  805.  
  806. LIGM = (0.0 0.0) 'DROIT' NY (0.0 (NY '*' DX)) ;
  807. LIGG = LIGM 'PLUS' ((-1 '*' NX '*' DX) 0.0) ;
  808. LIGD = LIGM 'PLUS' ((NX '*' DX) 0.0) ;
  809. DOM1 = LIGG 'REGLER' NX LIGM ;
  810. DOM2 = LIGM 'REGLER' NX LIGD ;
  811. DOMINT = DOM1 'ET' DOM2 ;
  812. CDOM = 'DIFF' ('CONTOUR' DOMINT) (LIGG 'ET' LIGD) ;
  813. *
  814. **** Creation of MODELS
  815. *
  816.  
  817. MDOMINT = 'MODELISER' DOMINT 'EULER' ;
  818. MDOM1 = 'MODELISER' DOM1 'EULER' ;
  819. MDOM2 = 'MODELISER' DOM2 'EULER' ;
  820. MLIGG = 'MODELISER' LIGG 'EULER' ;
  821. MLIGD = 'MODELISER' LIGD 'EULER' ;
  822. MCDOM = 'MODELISER' CDOM 'EULER' ;
  823.  
  824. TDOMINT = 'DOMA' MDOMINT 'VF' ;
  825. TDOM1 = 'DOMA' MDOM1 'VF' ;
  826. TDOM2 = 'DOMA' MDOM2 'VF' ;
  827. TLIGG = 'DOMA' MLIGG 'VF' ;
  828. TLIGD = 'DOMA' MLIGD 'VF' ;
  829. TCDOM = 'DOMA' MCDOM 'VF' ;
  830.  
  831. QDOMINT = TDOMINT . 'QUAF' ;
  832. QDOM1 = TDOM1 . 'QUAF' ;
  833. QDOM2 = TDOM2 . 'QUAF' ;
  834. QLIGD = TLIGD . 'QUAF' ;
  835. QLIGG = TLIGG . 'QUAF' ;
  836. QCDOM = TCDOM . 'QUAF' ;
  837.  
  838. 'ELIMINATION' QDOMINT (1.0D-3 '/' RAF) QDOM1 ;
  839. 'ELIMINATION' QDOMINT (1.0D-3 '/' RAF) QDOM2 ;
  840. 'ELIMINATION' QDOMINT (1.0D-3 '/' RAF) QLIGG ;
  841. 'ELIMINATION' QDOMINT (1.0D-3 '/' RAF) QLIGD ;
  842. 'ELIMINATION' QDOMINT (1.0D-3 '/' RAF) QCDOM ;
  843.  
  844. 'SI' GRAPH ;
  845. 'TRACER' DOMINT 'TITRE'
  846. ('CHAINE' 'Domaine, nbel=' ('NBEL' DOMINT)) ;
  847. 'FINSI' ;
  848.  
  849. *************************************************************************
  850. *************************************************************************
  851. **************** FIN MAILLAGE ******************************************
  852. *************************************************************************
  853. *************************************************************************
  854.  
  855. *************************************************************************
  856. *************************************************************************
  857. **************** INITIAL CONDITIONS *************************************
  858. *************************************************************************
  859. *************************************************************************
  860.  
  861. * Noms des variables
  862.  
  863. MOTRN = 'RN' ;
  864. MOTGNX = 'RUX' ;
  865. MOTGNY = 'RUY' ;
  866. MOTGN = 'MOTS' MOTGNX MOTGNY ;
  867. MOTVNX = 'UX' ;
  868. MOTVNY = 'UY' ;
  869. MOTVN = 'MOTS' MOTVNX MOTVNY ;
  870. MOTRET = 'RETN' ;
  871. MOTPN = 'PN' ;
  872.  
  873. LISTINCO = 'MOTS' MOTRN MOTGNX MOTGNY MOTRET ;
  874. LISTPRIM = 'MOTS' MOTRN MOTVNX MOTVNY MOTPN ;
  875.  
  876. GAMAIR = 1.4 ;
  877.  
  878. *
  879. ***** Left state
  880. *
  881.  
  882. ROL = 1.4 ;
  883. PL = 1.0 ;
  884. UL = 0.01 ;
  885.  
  886. HTL = ((GAMAIR '/' (GAMAIR '-' 1.0)) '*' (PL '/' ROL)) '+'
  887. (0.5 '*' UL '*' UL) ;
  888. SL = PL '/' (ROL '**' GAMAIR) ;
  889.  
  890. *
  891. ***** Right state
  892. *
  893.  
  894. ROR = 1.4 '/' 10. ;
  895. PR = PL ;
  896. UR = UL ;
  897.  
  898. *
  899. *** Boundary Conditions
  900. *
  901. *
  902. *
  903. * Subsonic inlet.
  904. * We impose the total enthalpy and the entropy.
  905. * The velocity direction is supposed to be 90 degrees
  906. *
  907. *
  908. CHPINL = 'MANUEL' 'CHPO' ('DOMA' MLIGG 'CENTRE') 2 'HT' HTL
  909. 'S' SL ;
  910. *
  911. * Subsonic outlet
  912. * We impose the pressure
  913. *
  914. CHPOUT = 'MANUEL' 'CHPO' ('DOMA' MLIGD 'CENTRE') 1 'PN' PR ;
  915. *
  916. * Point where we directly compute the boundary conditions
  917. *
  918. MAILLIM = ('DOMA' MLIGG 'CENTRE') 'ET' ('DOMA' MLIGD 'CENTRE') ;
  919. *
  920. RN0 = ('MANU' 'CHPO' ('DOMA' MDOM1 'CENTRE') 1 'SCAL' ROL
  921. 'NATURE' 'DISCRET') '+'
  922. ('MANU' 'CHPO' ('DOMA' MDOM2 'CENTRE') 1 'SCAL' ROR
  923. 'NATURE' 'DISCRET') ;
  924.  
  925. VN0 = ('MANU' 'CHPO' ('DOMA' MDOM1 'CENTRE') 2 'UX' UL
  926. 'UY' 0.0 'NATURE' 'DISCRET') 'ET'
  927. ('MANU' 'CHPO' ('DOMA' MDOM2 'CENTRE') 2 'UX' UR
  928. 'UY' 0.0 'NATURE' 'DISCRET') ;
  929.  
  930. PN0 = ('MANU' 'CHPO' ('DOMA' MDOM1 'CENTRE') 1 'SCAL' PL
  931. 'NATURE' 'DISCRET') 'ET'
  932. ('MANU' 'CHPO' ('DOMA' MDOM2 'CENTRE') 1 'SCAL' PR
  933. 'NATURE' 'DISCRET') ;
  934.  
  935. GAMN = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1 'SCAL' GAMAIR
  936. 'NATURE' 'DISCRET' ;
  937.  
  938. CSON = GAMN '*' (PN0 '/' RN0) ;
  939. CSON = CSON '**' 0.5 ;
  940.  
  941. GN0 RET0 = CONS RN0 VN0 PN0 GAMN ;
  942.  
  943. VN1 PN1 = 'PRIM' 'PERFMONO' RN0 GN0 RET0 GAMN ;
  944.  
  945. ERRO = 'MAXIMUM' (PN1 '-' PN0) 'ABS' ;
  946.  
  947. 'SI' (ERRO > 1.0D-6) ;
  948. 'MESSAGE' 'Problem in the ic file!!!'
  949. 'ERREUR' 5 ;
  950. 'FINSI' ;
  951. *
  952. **** Plot of IC
  953. *
  954. GRAVITE = 'MANUEL' 'CHPO' (TDOMINT . 'CENTRE') 2 'UX' 0.0 'UY' 0.0 ;
  955. * 'SI' (FAUX 'ET' GRAPH) ;
  956. 'SI' (GRAPH) ;
  957. CHM_RN = 'KCHA' MDOMINT 'CHAM' RN0 ;
  958. CHM_PN = 'KCHA' MDOMINT 'CHAM' PN0 ;
  959. CHM_VN = 'KCHA' MDOMINT 'CHAM' VN0 ;
  960. CHM_MN = 'KCHA' MDOMINT 'CHAM'
  961. ((('PSCAL' VN0 VN0 MOTVN MOTVN) '**' 0.5) '/' CSON) ;
  962.  
  963. 'TRACER' CHM_RN MDOMINT
  964. 'TITR' ('CHAINE' 'RN at t=' 0.0);
  965. 'TRACER' CHM_PN MDOMINT
  966. 'TITR' ('CHAINE' 'PN at t=' 0.0);
  967. 'TRACER' CHM_VN MDOMINT
  968. 'TITR' ('CHAINE' 'VN at t=' 0.0);
  969. 'TRACER' CHM_MN MDOMINT
  970. 'TITR' ('CHAINE' 'Mach at t=' 0.0);
  971.  
  972. 'FINSI' ;
  973.  
  974. *************************************************************************
  975. *************************************************************************
  976. **************** FIN INITIAL CONDITIONS *********************************
  977. *************************************************************************
  978. *************************************************************************
  979.  
  980. *************************************************************************
  981. *************************************************************************
  982. ************** COMPUTATION OF THE SOLUTION ******************************
  983. *************************************************************************
  984. *************************************************************************
  985.  
  986. RV = 'TABLE' ;
  987. RV . 'FREQI' = 1 ;
  988. RV . 'MODEL' = MDOMINT ;
  989. *
  990. **** Conservative variables / primitive variables
  991. *
  992. *
  993. RV . 'LISTCONS' = 'MOTS' MOTRN MOTGNX MOTGNY MOTRET ;
  994. RV . 'LISTPRIM' = 'MOTS' MOTRN MOTVNX MOTVNY MOTPN ;
  995. *
  996. RV . 'RN0' = RN0 ;
  997. RV . 'GN0' = GN0 ;
  998. RV . 'RET0' = RET0 ;
  999. *
  1000. **** Gas property/gravity
  1001. *
  1002. RV . 'PGAS' = 'TABLE' ;
  1003. RV . 'PGAS' . 'GAMN' = GAMN ;
  1004. RV . 'PGAS' . 'MU' = 0.0D0 ;
  1005. RV . 'PGAS' . 'R' = 1.0D0 ;
  1006. RV . 'PGAS' . 'LAMBDA' = 0.0D0 ;
  1007. RV . 'GRAVITY' = (GRAVITE '*' 0.0) ;
  1008. *
  1009. * Table for BC (injection part)
  1010. *
  1011. RV . 'PROLIM' = 'TABLE' ;
  1012. RV . 'PROLIM' . 'N' = 2 ;
  1013. RV . 'PROLIM' . 'CLN' = 'TABLE' ;
  1014. RV . 'PROLIM' . 'CLN' . 1 = 'MOTS' INSU ;
  1015. RV . 'PROLIM' . 'CLN' . 2 = 'MOTS' OUTP ;
  1016. RV . 'PROLIM' . 'MODELN' = 'TABLE' ;
  1017. RV . 'PROLIM' . 'MODELN' . 1 = MLIGG ;
  1018. RV . 'PROLIM' . 'MODELN' . 2 = MLIGD ;
  1019. RV . 'PROLIM' . 'CHPOLN' = 'TABLE' ;
  1020. RV . 'PROLIM' . 'CHPOLN' . 1 = CHPINL ;
  1021. RV . 'PROLIM' . 'CHPOLN' . 2 = CHPOUT ;
  1022. *
  1023. * Table for BC
  1024. *
  1025. * RVX . 'DIFTIMP' : boundary condition on temperature
  1026. * RVX . 'MDIFCT' : mesh in which we impose the temperature
  1027. * arising from the convective BC
  1028. * RVX . 'DIFGTIMP' : boundary condition on gradient of temperature
  1029. * RVX . 'DIFVIMP' : boundary condition on speed
  1030. * RVX . 'MDIFCV' : mesh in which we impose the velocity
  1031. * arising from the convective BC
  1032. * RVX . 'DIFGVIMP' : boundary condition on gradient of speed
  1033. * RVX . 'DIFTAUI' : boundary condition on constraint tensor
  1034. * RVX . 'DIFQIMP' : boundary condition on heat flux
  1035. *
  1036. CHPVID MATVID = 'KOPS' 'MATRIK' ;
  1037. RV . 'DIFTIMP' = CHPVID ;
  1038. RV . 'DIFGTIMP' = 'MANUEL' 'CHPO' (TCDOM . 'CENTRE') 2 'P1DX' 0.0
  1039. 'P1DY' 0.0 ;
  1040. RV . 'DIFQIMP' = 'MANUEL' 'CHPO' (TCDOM . 'CENTRE') 2 'UX' 0.0
  1041. 'UY' 0.0 ;
  1042. RV . 'MDIFCT' = ('DOMA' MLIGG 'CENTRE') 'ET'
  1043. ('DOMA' MLIGD 'CENTRE') ;
  1044. RV . 'DIFVIMP' = CHPVID ;
  1045. RV . 'DIFGVIMP' = 'MANUEL' 'CHPO' (TCDOM . 'CENTRE') 4 'P1DX' 0.0
  1046. 'P1DY' 0.0 'P2DX' 0.0 'P2DY' 0.0 ;
  1047. RV . 'DIFTAUI' = 'MANUEL' 'CHPO' (TCDOM . 'CENTRE') 3 'TXX' 0.0
  1048. 'TXY' 0.0 'TYY' 0.0 ;
  1049. RV . 'MDIFCV' = ('DOMA' MLIGG 'CENTRE') 'ET'
  1050. ('DOMA' MLIGD 'CENTRE') ;
  1051. *
  1052. **** Numerical parameters
  1053. *
  1054. *
  1055. * Variable to compute L2 error
  1056. *
  1057. RV . 'LISTERR' = 'MOTS' MOTGNX MOTGNY ;
  1058. RV . 'LISTRHO' = 'MOTS' MOTRN ;
  1059. * RV . 'LISTERR' = 'MOTS' MOTRET ;
  1060. *
  1061. * Upwind scheme
  1062. *
  1063. * RV . 'METHOD' = 'RUSANOLM' ;
  1064. RV . 'METHOD' = 'AUSMPLM' ;
  1065. * RV . 'METHOD' = 'HLLCLM' ;
  1066. *
  1067. * Low-Mach Cut off
  1068. *
  1069. CO = UL '*' 10.D-2 ;
  1070. RV . 'CUTOFF' = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1 'SCAL' CO
  1071. 'NATURE' 'DISCRET' ;
  1072. *
  1073. * Reconstruction/limiter
  1074. * Time accuracy (1 or 2)
  1075. * Iterations
  1076. * Final time
  1077. RV . 'SPACEA' = 2 ;
  1078. * RV . 'LIMITER' = 'NOLIMITE' ;
  1079. RV . 'LIMITER' = 'LIMITEUR' ;
  1080. RV . 'TIMEA' = 2 ;
  1081. *
  1082. **** Physical time
  1083. *
  1084. RV . 'T0' = 0 ;
  1085. * RV . 'CFL' = 20. ;
  1086. RV . 'TFINAL' = L '/' (2 '*' UL) ;
  1087. * RV . 'CFL' = 20. ;
  1088. RV . 'DTPS' = 0.5 '*' (DX '/' UL) ;
  1089. RV . 'NITER' = -1 ;
  1090. *
  1091. **** Dual time
  1092. *
  1093. * Safety factor for the dual time step
  1094. * Max. Dual time iterations
  1095. * Relative error
  1096. *
  1097. RV . 'DCFLERR' = 'PROG' 16. -16. ;
  1098. RV . 'DCFL' = 'PROG' 1.d6 1.d6 ;
  1099. RV . 'NDTITER' = 100 ;
  1100. RV . 'RELERR' = FAUX ;
  1101. RV . 'EPSDT' = 1.0D-10 ;
  1102. *
  1103. **** Jacobi iterations
  1104. *
  1105. * RV . 'TYPEJAC' = 'PJACO' ;
  1106. RV . 'TYPEJAC' = 'LJACOFB' ;
  1107. RV . 'NJACLERR' = 'PROG' 16 -16 ;
  1108. RV . 'NJACITER' = 'PROG' 15 15 ;
  1109. *
  1110. **** Parameters for PROCPT
  1111. *
  1112. RV . 'PROCPT' = 'TABLE' ;
  1113. RV . 'PROCPT' . 'RMAX' = 'PROG' ;
  1114. RV . 'PROCPT' . 'RMIN' = 'PROG' ;
  1115. RV . 'PROCPT' . 'RMED' = 'PROG' ;
  1116. RV . 'PROCPT' . 'PMAX' = 'PROG' ;
  1117. RV . 'PROCPT' . 'PMIN' = 'PROG' ;
  1118. RV . 'PROCPT' . 'PMED' = 'PROG' ;
  1119. RV . 'PROCPT' . 'TPS' = 'PROG' ;
  1120. *
  1121. PNSSM RV ;
  1122.  
  1123. *************************************************************************
  1124. *************************************************************************
  1125. ************** FIN COMPUTATION OF THE SOLUTION **************************
  1126. *************************************************************************
  1127. *************************************************************************
  1128.  
  1129. *************************************************************************
  1130. *************************************************************************
  1131. ************** POST TREATMENT *******************************************
  1132. *************************************************************************
  1133. *************************************************************************
  1134.  
  1135. *
  1136. RN = RV . 'RESULTS' . 'RN' ;
  1137. GN = RV . 'RESULTS' . 'GN' ;
  1138. RET = RV . 'RESULTS' . 'RET' ;
  1139. GAMN = RV . 'PGAS' . 'GAMN' ;
  1140.  
  1141. *
  1142. *** Convergence evolution inside of each iteration
  1143. *
  1144.  
  1145. LISTITER = RV . 'RESULTS' . 'LISTITER' ;
  1146. LISTITDT = RV . 'RESULTS' . 'LISTITDT' ;
  1147. LISTLINF = RV . 'RESULTS' . 'LISTLINF' ;
  1148.  
  1149. 'SI' GRAPH ;
  1150. NITERE = 'DIME' LISTITER ;
  1151. I1 = 1 ;
  1152. AA = 'PROG' ;
  1153. BB = 'PROG' ;
  1154. CC = 'PROG' ;
  1155. DD = 'PROG' ;
  1156. 'REPETER' BLITER NITERE ;
  1157. I0 = I1 ;
  1158. I1 = 'EXTRAIRE' LISTITER &BLITER ;
  1159. 'SI' (I1 'EGA' I0) ;
  1160. AA = AA 'ET' ('PROG' ('EXTRAIRE' LISTITDT &BLITER)) ;
  1161. BB = BB 'ET' ('PROG' ('EXTRAIRE' LISTLINF &BLITER)) ;
  1162. 'SINON' ;
  1163. everr = 'EVOL' 'MANU' 'niter' AA 'Log(Linf)'
  1164. (('LOG' (BB '+' ('PROG' ('DIME' BB) '*' 1.0D-20))) '/'
  1165. ('LOG' 10.)) ;
  1166. 'DESSIN' everr 'TITRE' ('CHAINE' 'Convergence at iter '
  1167. (I1 '-' 1)) 'NCLK' ;
  1168. AA = 'PROG' ;
  1169. BB = 'PROG' ;
  1170. CC = 'PROG' ;
  1171. DD = 'PROG' ;
  1172. 'FINSI' ;
  1173. 'FIN' BLITER ;
  1174. everr = 'EVOL' 'MANU' 'niter' AA 'Log(Linf)'
  1175. (('LOG' (BB '+' ('PROG' ('DIME' BB) '*' 1.0D-20))) '/'
  1176. ('LOG' 10.)) ;
  1177. 'DESSIN' everr 'TITRE' ('CHAINE' 'Convergence at iter '
  1178. (I1)) 'NCLK' ;
  1179. 'FINSI' ;
  1180. *
  1181. *
  1182. **** The mesh
  1183. *
  1184.  
  1185. 'SI' GRAPH ;
  1186. 'TRACER' DOMINT 'TITR' 'Maillage' ;
  1187. 'FINSI' ;
  1188.  
  1189. *
  1190. **** Initial conditions
  1191. *
  1192.  
  1193. 'SI' GRAPH ;
  1194.  
  1195. RN0 = RV . 'RN0' ;
  1196. GN0 = RV . 'GN0' ;
  1197. RET0 = RV . 'RET0' ;
  1198. VN PN0 = 'PRIM' 'PERFMONO'
  1199. RN0 GN0 RET0 GAMN 'TRICHE' ;
  1200.  
  1201. CHM_RN = 'KCHA' MDOMINT 'CHAM' RN0 ;
  1202. CHM_VN = 'KCHA' MDOMINT 'CHAM' VN ;
  1203. CHM_PN = 'KCHA' MDOMINT 'CHAM' PN0 ;
  1204.  
  1205. 'TRAC' CHM_RN MDOMINT ('CONTOUR' DOMINT)
  1206. 'TITR' ('CHAINE' 'rho at t=' 0.0) ;
  1207. 'TRAC' CHM_VN MDOMINT ('CONTOUR' DOMINT)
  1208. 'TITR' ('CHAINE' 'v at t= ' 0.0) ;
  1209. 'TRAC' CHM_PN MDOMINT ('CONTOUR' DOMINT)
  1210. 'TITR' ('CHAINE' 'p at t= ' 0.0) ;
  1211.  
  1212. 'FINSI' ;
  1213.  
  1214. *
  1215. **** The 2D graphics
  1216. *
  1217.  
  1218. NOMVEL = 'MOTS' 'UX' 'UY' ;
  1219. VN PN = 'PRIM' 'PERFMONO' RN GN RET GAMN 'TRICHE' ;
  1220. CN2 = GAMN '*' (PN '/' RN) ;
  1221. VN2 = 'PSCAL' VN VN NOMVEL NOMVEL ;
  1222. MACHN2 = VN2 '/' CN2 ;
  1223. MACHN = MACHN2 '**' 0.5 ;
  1224.  
  1225. HTN = (GAMN '/' (GAMN '-' 1.0)) '*' (PN '/' RN) ;
  1226. ECIN = 0.5 '*' ('PSCAL' VN VN NOMVEL NOMVEL) ;
  1227. HTN = HTN '+' ECIN ;
  1228.  
  1229. SN = PN '/' (RN '**' 1.4) ;
  1230.  
  1231. TPS = RV . 'RESULTS' . 'TPS' ;
  1232. 'SI' ('EGA' ('TYPE' TPS) 'CHPOINT') ;
  1233. TPS = 'MINIMUM' TPS ;
  1234. 'FINSI' ;
  1235.  
  1236. 'SI' GRAPH ;
  1237.  
  1238. CHM_RN = 'KCHA' MDOMINT 'CHAM' RN ;
  1239. CHM_VN = 'KCHA' MDOMINT 'CHAM' VN ;
  1240. CHM_PN = 'KCHA' MDOMINT 'CHAM' PN ;
  1241. CHM_MN = 'KCHA' MDOMINT 'CHAM' MACHN ;
  1242. CHM_HTN = 'KCHA' MDOMINT 'CHAM' HTN ;
  1243. CHM_SN = 'KCHA' MDOMINT 'CHAM' SN ;
  1244.  
  1245. 'TRAC' CHM_RN MDOMINT ('CONTOUR' DOMINT)
  1246. 'TITR' ('CHAINE' 'rho at t=' TPS) ;
  1247. 'TRAC' CHM_VN MDOMINT ('CONTOUR' DOMINT)
  1248. 'TITR' ('CHAINE' 'v at t= ' TPS) ;
  1249. 'TRAC' CHM_PN MDOMINT ('CONTOUR' DOMINT)
  1250. 'TITR' ('CHAINE' 'p at t= ' TPS) ;
  1251. 'TRAC' CHM_MN MDOMINT ('CONTOUR' DOMINT)
  1252. 'TITR' ('CHAINE' 'Mach at t= ' TPS) ;
  1253. 'TRAC' CHM_HTN MDOMINT ('CONTOUR' DOMINT)
  1254. 'TITR' ('CHAINE' 'ht at t= ' TPS ' hl =' HTL) ;
  1255. 'TRAC' CHM_SN MDOMINT ('CONTOUR' DOMINT)
  1256. 'TITR' ('CHAINE' 'sn at t= ' TPS ' hl =' SL) ;
  1257.  
  1258. VECN = 'VECTEUR' VN (1. '/' RAF) 'UX' 'UY' 'JAUNE' ;
  1259. 'TRACER' DOMINT VECN ('CONTOUR' DOMINT) 'TITR'
  1260. ('CHAINE' 'v at t= ' TPS) ;
  1261.  
  1262. 'OPTION' 'ISOV' 'LIGN' ;
  1263. RNV = 'ELNO' TDOMINT ('KCHT' TDOMINT 'SCAL' 'CENTRE' RN) ;
  1264. 'TRACER' DOMINT RNV ('CONTOUR' DOMINT) 15 'TITRE' 'ro';
  1265. PNV = 'ELNO' TDOMINT ('KCHT' TDOMINT 'SCAL' 'CENTRE' PN) ;
  1266. 'TRACER' DOMINT PNV ('CONTOUR' DOMINT) 15 'TITRE' 'p' ;
  1267.  
  1268. 'FINSI' ;
  1269.  
  1270. *
  1271. **** Evolution objects
  1272. *
  1273. *
  1274. * Position of the contact
  1275. *
  1276.  
  1277. TYPELIM = RV . 'LIMITER' ;
  1278.  
  1279. POS = TPS '*' UL ;
  1280. LISTPOS = 'PROG' (POS '-' (DX '/' 2)) (POS '+' (DX '/' 2)) ;
  1281.  
  1282. TITOLO = 'CHAINE' 'Nbel = ' ('NBEL' DOMINT) ', OT = '
  1283. (RV . 'TIMEA') ', ' RV . 'METHOD' ;
  1284. 'SI' VRAI ;
  1285. TITOLO = 'CHAINE' TITOLO ', OE = 2, ' TYPELIM ;
  1286. 'SINON' ;
  1287. TITOLO = 'CHAINE' TITOLO ', OE = 1 ' ;
  1288. 'FINSI' ;
  1289.  
  1290. 'OPTION' 'ELEM' QUA4 ;
  1291.  
  1292. AA = 'POIN' ('COORDONNEE' 1 (TDOMINT . 'CENTRE')) 'MINI' ;
  1293. DCEN = AA 'ET' ('DIFF' (TDOMINT . 'CENTRE') AA) ;
  1294. 'ORDONNER' DCEN ;
  1295. LIGEVOL = 'QUELCONQUE' DCEN 'SEG2' ;
  1296.  
  1297. LISTX = 'EXTRAIRE' ('EVOL' 'CHPO' ('COORDONNEE' 1 LIGEVOL) 'SCAL'
  1298. LIGEVOL) 'ORDO' ;
  1299.  
  1300. RNEX = (ROL * (('COORDONNEE' 1 (TDOMINT . 'CENTRE')) 'MASQUE'
  1301. 'INFERIEUR' POS)) '+'
  1302. (ROR * (('COORDONNEE' 1 (TDOMINT . 'CENTRE')) 'MASQUE'
  1303. 'SUPERIEUR' POS)) ;
  1304. 'SI' GRAPH ;
  1305. EVR = 'EVOL' 'CHPO' RN 'SCAL' LIGEVOL ;
  1306. LISTR = 'EXTRAIRE' EVR 'ORDO' ;
  1307. EVR = 'EVOL' 'MANU' 'x' LISTX 'Density' LISTR ;
  1308. EVPOS = 'EVOL' 'ROUG' 'MANU' 'x' LISTPOS ' ' ('PROG'
  1309. ('MINI' LISTR) ('MAXIMUM' LISTR)) ;
  1310. EVREX = 'EVOL' 'CHPO' RNEX 'SCAL' LIGEVOL ;
  1311. LISTR = 'EXTRAIRE' EVREX 'ORDO' ;
  1312. EVREX = 'EVOL' 'VERT' 'MANU' 'x' LISTX 'Density' LISTR ;
  1313. 'DESSIN' (EVR 'ET' EVPOS 'ET' EVREX) 'MIMA' 'TITRE' TITOLO 'GRILL'
  1314. 'YBOR' 0.0 1.6 ;
  1315.  
  1316. EVP = 'EVOL' 'CHPO' PN 'SCAL' LIGEVOL ;
  1317. LISTP = 'EXTRAIRE' EVP 'ORDO' ;
  1318. EVP = 'EVOL' 'MANU' 'x' LISTX 'Pressure' LISTP ;
  1319. EVPOS = 'EVOL' 'ROUG' 'MANU' 'x' LISTPOS ' ' ('PROG'
  1320. ('MINI' LISTP) ('MAXIMUM' LISTP)) ;
  1321. 'DESSIN' (EVP 'ET' EVPOS) 'MIMA' 'TITRE' TITOLO 'GRILL' ;
  1322.  
  1323. EVH = 'EVOL' 'CHPO' HTN 'SCAL' LIGEVOL ;
  1324. LISTH = 'EXTRAIRE' EVH 'ORDO' ;
  1325. EVH = 'EVOL' 'MANU' 'x' LISTX 'Total enthalpie' LISTH ;
  1326. EVPOS = 'EVOL' 'ROUG' 'MANU' 'x' LISTPOS ' ' ('PROG'
  1327. ('MINI' LISTH) ('MAXIMUM' LISTH)) ;
  1328. 'DESSIN' (EVH 'ET' EVPOS) 'MIMA' 'TITRE'
  1329. ('CHAINE' TITOLO ', hl =' HTL) 'GRILL' ;
  1330.  
  1331. EVS = 'EVOL' 'CHPO' SN 'SCAL' LIGEVOL ;
  1332. LISTS = 'EXTRAIRE' EVS 'ORDO' ;
  1333. EVS = 'EVOL' 'MANU' 'x' LISTX 'Entropie' LISTS ;
  1334. EVPOS = 'EVOL' 'ROUG' 'MANU' 'x' LISTPOS ' ' ('PROG'
  1335. ('MINI' LISTS) ('MAXIMUM' LISTS)) ;
  1336. 'DESSIN' (EVS 'ET' EVPOS) 'MIMA' 'TITRE'
  1337. ('CHAINE' TITOLO ', sl =' SL) 'GRILL' ;
  1338.  
  1339. 'FINSI' ;
  1340.  
  1341. *
  1342. **** Test de convergence
  1343. *
  1344. AA = 'EXTRAIRE' LISTLINF ('DIME' LISTLINF) ;
  1345. 'SI' (AA > 1.0D-8) ;
  1346. 'ERRE' 'Probleme en KONV' ;
  1347. 'FINSI' ;
  1348.  
  1349. ERRO = 'ABS' (RN '-' RNEX) ;
  1350. ERRO = 'EVOL' 'CHPO' ERRO LIGEVOL ;
  1351. ERRO = 'PRIM' ERRO ;
  1352. ERRO = 'EXTRAIRE' ERRO 'ORDO' ;
  1353. ERRO = 'EXTRAIRE' ERRO ('DIME' ERRO) ;
  1354. 'SI' (ERRO > 1.0D-1) ;
  1355. 'ERRE' 'Probleme en KONV' ;
  1356. 'FINSI' ;
  1357.  
  1358. 'FIN' ;
  1359.  
  1360.  
  1361.  
  1362.  
  1363.  

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