Télécharger la plaquette Cast3M


fr en

Télécharger benchmark_imst.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  benchmark_imst.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ***************************************************************
  5. *                                                             *
  6. * NOM         : benchmark_imst.dgibi                          *
  7. * DESCRIPTION : fiche de validation CASTEM2000                *
  8. *               Mecanique des Fluides                         *
  9. *               Convection naturelle laminaire                *
  10. *               en cavite rectangulaire (2D)                  *
  11. * FONCTIONS                                                   *
  12. *   TESTEES   : 'NS' algorithme semi-implicite                *
  13. *                    option  CENTREE                          *
  14. *                    approximation de Boussinesq              *
  15. *                    elements QUA8 'MACRO'                    *
  16. *                                                             *
  17. * AUTEUR      : Stephane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)      *
  18. *               e-mail : gounand@semt2.smts.cea.fr            *
  19. *                                                             *
  20. ***************************************************************
  21. *                                                             *
  22. * VERSION     : v1, 19/9/97, version initiale                 *
  23. * HISTORIQUE  : v1, 19/9/97, creation                         *
  24. * HISTORIQUE  : 15/06/2014 passage EQPR -> EQEX               *
  25. * HISTORIQUE  : 11/07/2023 fixer la pression en 1 point       *
  26. * HISTORIQUE  :                                               *
  27. * HISTORIQUE  :                                               *
  28. *                                                             *
  29. ***************************************************************
  30. *                                                             *
  31. * Priere de completer la partie HISTORIQUE apres modification *
  32. * du jeu de donnees afin de faciliter la maintenance          *
  33. *                                                             *
  34. ***************************************************************
  35. 'OPTI' echo 0 ;
  36. *
  37. * DESCRIPTION DETAILLEE :
  38. *
  39. *  --  IMST --  Convection naturelle laminaire à Prandtl (Pr) = 0
  40. *                Algorithme semi implicite
  41. *                operateur NS  (option CENTREE)
  42. *
  43. *              phi = 0
  44. *        ------------------------
  45. *        |                      |
  46. * T1 = 0 | H                    | T2 = A
  47. *        |          L           |
  48. *        ------------------------
  49. *              phi = 0
  50. *
  51. *   A = L/H = 4  l'allongement
  52. *   U = 0 sur les parois
  53. *   Le profil de temperature est impose T = T1 + (T2-T1)*(x/L)
  54. *
  55. *   La structure de l'ecoulement ne depend que
  56. *   du nombre de Grashof Gr :
  57. *   Gr = g beta (T2 - T1) H**3 A / nu**2
  58. *
  59. *      ^
  60. *      | Psi max (valeur max de la fonction de courant)
  61. *      |
  62. *      |
  63. *     -|                     o o
  64. *      |  unsteady ->      o   v
  65. *      |   (3 cells)     .     v
  66. *      | -----------    o      v         . o
  67. *      |  steady ->    .^      v      o
  68. *      |  (3 cells)   . ^      v  .
  69. *      | ----------- o  ^      o
  70. *     -|                ^   .
  71. *      |  steady -> .   ^ .       <-- steady (2 cells)
  72. *      |                o
  73. *      | (1 cell)  .
  74. *      |
  75. *      |          .
  76. *      |         o
  77. *     -|-------------------------------------------> Gr
  78. *               I     I     I
  79. *              1.E4  2.E4  3.E4
  80. **************************************************************
  81. * References : - Gounand, S. 1997
  82. *                Simulation numerique de cellules
  83. *                thermoconvectives.
  84. *                Rapport CEA/DMT 97/357.
  85. *              - Le Garrec, S. 1988
  86. *                Convection naturelle en cavite pour des
  87. *                fluides à faible nombre de Prandtl.
  88. *                Rapport CEA/DMT 88/216.
  89. *              - Le Garrec, S. & Magnaud, J.-P.
  90. *                Numerical Simulation of oscillatory
  91. *                convection in low Prandtl Fluids with
  92. *                TRIO-EF.
  93. *                A GAMM Workshop Edited by Bernard Roux
  94. *                Notes on Numerical Fluid Mechanics,Volume 27
  95. *                pp. 189-199. (Vieweg, Braunschweig 1990)
  96. *
  97. ***************************************************************
  98. *
  99. * JEU DE DONNEE :
  100. *
  101. prg = ('CHAINE' 'Benchmark Imst') ;
  102. **** Options
  103. * court = VRAI  : le cas-test s'effectue en une dizaine
  104. *                 de secondes.
  105. *                 Gr = 5000 => une cellule stationnaire
  106. *                 la comparaison s'effectue sur le profil
  107. *                 Uy(x) le long de la mediane horizontale.
  108. *                 On verifie egalement que l'ecoulement
  109. *                 est centrosymetrique.
  110. * court = FAUX  : maillage raffine ;
  111. *                 Gr = 40000 => bifurcation 3->2 cellules
  112. *                 environ 1h30 CPU avec les graphiques
  113. *
  114. * graph = VRAI  : on trace les graphiques (historiques, champs)
  115. *                 On a plus de calculs avec cette option
  116. *                 (fonction de courant, convergence...)
  117. *                 On mettra donc graph = faux pour la fiche
  118. *                 de validation
  119. *
  120. * inta  = VRAI si on est en interactif ;  0 sinon
  121. *
  122. * discconv : option de discretisation des termes de
  123. *            convection pour les equations de Navier-Stokes
  124. *            choix possibles : 'CENTREE',  'SUPG'...
  125. *            voir notice de 'NS'
  126. court = VRAI ;
  127. graph = FAUX ;
  128. inta =  FAUX ;
  129. typelem = 'QUA8' ;
  130. discconv = 'CENTREE' ;
  131.  
  132. DISCR  = 'MACRO'  ;
  133. KPRESS = 'CENTRE' ;
  134. * betap : parametre de stabilisation de la pression
  135. BETAP  = 1.       ;
  136.  
  137. 'OPTION' 'DIME' 2 'ELEM' typelem ;
  138.  
  139. 'SI' ('NON' inta) ;
  140.   'OPTION' 'TRAC' 'PSC';
  141. 'SINON' ;
  142. **od = 'TEXTE' 'OPTION DONN 3' ;
  143. 'FINSI' ;
  144.  
  145. *************************************************************
  146. **** Definition de quelques procedures :
  147. **    -> SUMMARY : affiche un resume des parametres
  148. **    -> CALCPSI : calcul de la fonction de courant
  149. **                 (appelee par 'EXEC' via la table 'EQEX')
  150. **    -> MAJHIST : mise à jour d'historiques
  151. **                 (appelee par 'EXEC' via la table 'EQEX')
  152. **    -> TRHIST  : trace des historiques
  153. **    -> TRCHAMP : trace des differents champs
  154. **                 (appelee par 'EXEC' via la table 'EQEX')
  155. ** NB : pour simplifier, on ne transmet pas les parametres à
  156. **      SUMMARY. C'est de la mauvaise programmation...
  157. **      Les autres procedures sont correctes
  158. **      sous toutes reserves
  159. *
  160. ** SUMMARY
  161. 'DEBPROC' SUMMARY ;
  162. 'MESSAGE' ('CHAINE' 'Programme :' ' ' prg) ;
  163. 'SAUTER' 'LIGNE' ;
  164. 'MESSAGE' ('CHAINE' '***  GEOMETRIE :') ;
  165. 'MESSAGE' ('CHAINE' 'Boite carree :' ' ' l ' x' h) ;
  166. 'MESSAGE' ('CHAINE' '      soient :' '  ' nl 'x' nh '='  (nl '*' nh) ' elements' ' ' typelem) ;
  167. 'MESSAGE' ('CHAINE' '      soient :' '  ' (nl '*' nh '*' 4) ' elements lineaires') ;
  168. 'SAUTER' 'LIGNE' ;
  169. 'MESSAGE' ('CHAINE' '***  PHYSIQUE :') ;
  170. 'MESSAGE' ('CHAINE' 'Rapport d_aspect : ' ' ' ('ENTIER' aspect)) ;
  171. 'MESSAGE' ('CHAINE' 'Nombre de Grashof :' ' '('ENTIER'  Gr) ) ;
  172. 'MESSAGE' ('CHAINE' 'Temperature de reference :' ' ' tref ' comprise entre 0 et '  ('ENTIER' aspect)) ;
  173.   'SAUTER' 'LIGNE' ;
  174. 'MESSAGE' ('CHAINE' '***  NUMERIQUE :') ;
  175. 'MESSAGE' ('CHAINE' 'Discretisation' ' ' discconv ' des termes de convection.');
  176. 'MESSAGE' ('CHAINE' ' Alfa (mult. pasdt.) =' ' ' alpha) ;
  177. 'MESSAGE' ('CHAINE' ' Beta (stab. pression) =' ' ' betap) ;
  178. 'SAUTER' 'LIGNE' ;
  179. rvpdt = rv . 'PASDETPS' ;
  180. 'MESSAGE' ('CHAINE' '***  TEMPS :') ;
  181. 'MESSAGE' ('CHAINE' 'Nb. iterations :' ' ' ((rvpdt . 'NUPASDT') '-' 1)) ;
  182. 'MESSAGE' ('CHAINE' 'Temps physique :' ' ' (rvpdt . 'TPS')) ;
  183. 'MESSAGE' ('CHAINE' 'Temps CPU      :' ' ' ('ENTIER' (thist . 'tpscpu')) ' ms') ;
  184. 'FINPROC' ;
  185.  
  186.  
  187. ** CALCPSI
  188. 'DEBPROC' CALCPSI ;
  189. 'ARGUMENT' trx*'TABLE   ' ;
  190.  as2 ama1  = 'KOPS' 'MATRIK' ;
  191. trv = (trx . 'EQEX') ;
  192. tinco = (trv . 'INCO') ;
  193. $dom = (trx . 'DOMZ') ;
  194. cdom = 'CONTOUR' ('DOMA' $dom 'MAILLAGE') ;
  195. vitesse  = 'KCHT' $dom 'VECT' 'SOMMET' (tinco . 'UN') ;
  196. sw = 'KOPS' vitesse 'ROT' $dom ;
  197. rk = 'EQEX' $dom 'OPTI' 'EF' 'IMPL' 'ZONE' $dom 'OPER' 'LAPN' 1.0D0
  198.         'INCO' 'PSI' 'ZONE' $dom 'OPER' 'FIMP' sw
  199.         'INCO' 'PSI' 'CLIM' 'PSI' 'TIMP' cdom 0.0D0 ;
  200. rk . 'INCO' = 'TABLE' 'INCO' ;
  201. rk . 'INCO' . 'PSI' = 'KCHT' $dom 'SCAL' 'SOMMET' 0.0D0 ;
  202. *-*-*-*- A remplacer par EXEC qd ca marchera...
  203. EXEC rk ;
  204. *-*-*-*-
  205. tinco . 'PSI' = (rk . 'INCO' . 'PSI') ;
  206.  
  207. 'FINPROC' as2 ama1 ;
  208.  
  209. ** MAJHIST
  210. 'DEBPROC' MAJHIST ;
  211. 'ARGUMENT' trx*'TABLE   ' ;
  212.  as2 ama1  = 'KOPS' 'MATRIK' ;
  213. trv = (trx . 'EQEX') ;
  214. tinco = (trv . 'INCO') ;
  215. th    = (trv . 'HISTOIRE') ;
  216. tp    = (trv . 'PASDETPS') ;
  217. vit    = (tinco . 'UN') ;  vitm1 = (tinco . 'UN-1') ;
  218. vitref = (tinco . 'UREF') ;
  219. psi    = (tinco . 'PSI') ;
  220. npdt = (tp . 'NUPASDT') ;
  221. * Calcul de l'erreur relative (norme Linfini) entre deux
  222. * champs de vitesse successifs
  223. errrel = ('MAXIMUM' (vit '-' vitm1) 'ABS') '/' vitref ;
  224. 'SI' ('MULT' npdt (trv . 'FIDT')) ;
  225.   'MESSAGE' ('CHAINE' 'Pdt : ' npdt '  |Un - Un_1| Linf = ' errrel) ;
  226. 'FINSI' ;
  227. tinco . 'UN-1' = 'COPIER' vit ;
  228. * Mise à jour des differentes listes
  229. mpsi = 'MAXIMUM' psi 'ABS' ;
  230. th . 'ltps'  = ((th . 'ltps') 'ET' ('PROG' (tp . 'TPS'))) ;
  231. th . 'lpsi'  = ((th . 'lpsi') 'ET' ('PROG' mpsi)) ;
  232. th . 'lerr'  = ((th . 'lerr') 'ET' ('PROG' errrel)) ;
  233. thl = (th . 'lpdt') ;
  234. thl . 'conv' = ((thl . 'conv') 'ET' ('PROG' (tp . 'DTCONV'))) ;
  235. thl . 'diff' = ((thl . 'diff') 'ET' ('PROG' (tp . 'DTDIFU'))) ;
  236. thl . 'min'  = ((thl . 'min')  'ET' ('PROG' (tp . 'DELTAT-1'))) ;
  237. 'FINPROC' as2 ama1 ;
  238.  
  239.  
  240. ** TRHIST
  241. 'DEBPROC' TRHIST ;
  242. 'ARGUMENT' th*'TABLE   ' ;
  243. * On enleve le premier indice des listes
  244. hltps = 'ENLEVER' (th . 'ltps') 1 ;
  245. hlpsi = 'ENLEVER' (th . 'lpsi') 1 ;
  246. hlerr = 'ENLEVER' (th . 'lerr') 1 ;
  247. thl = (th . 'lpdt') ;
  248. hlpdt  = 'ENLEVER' (thl . 'min')  1 ;
  249. hlconv = 'ENLEVER' (thl . 'conv') 1 ;
  250. hldiff = 'ENLEVER' (thl . 'diff') 1 ;
  251. * Pas de temps = fn(t)
  252. evpdt = 'EVOL' 'MANU' 'TEMPS' hltps 'PDT' hlpdt ;
  253. 'DESSIN' evpdt 'MIMA' 'TITR' 'Pas de temps' ;
  254. * LOG10(pas de temps) = fn(t)
  255. ltco = ('LOG' (hlconv '*' alpha)) '/' ('LOG' 10.0D0) ;
  256. ltdi = ('LOG' (hldiff '*' alpha)) '/' ('LOG' 10.0D0) ;
  257. ltmi = ('LOG' hlpdt) '/' ('LOG' 10.0D0) ;
  258. evco = 'COULEUR' ('EVOL' 'MANU' 'TEMPS' hltps 'dtconv' ltco) 'BLEU' ;
  259. evdi = 'COULEUR' ('EVOL' 'MANU' 'TEMPS' hltps 'dtdiff' ltdi) 'ROUG' ;
  260. evmi = 'COULEUR' ('EVOL' 'MANU' 'TEMPS' hltps 'dt' ltmi) 'BLAN' ;
  261. tabt = 'TABLE' ;  tabt . 'TITRE' = 'TABLE' ;
  262. tabt . 1 = 'TIRC' ;
  263. tabt . 2 = 'TIRL' ;
  264. tabt . 3 = 'NOLI' ;
  265. tabt . 'TITRE' . 1 = 'CHAINE' 'Pdt convection' ;
  266. tabt . 'TITRE' . 2 = 'CHAINE' 'Pdt diffusion' ;
  267. tabt . 'TITRE' . 3 = 'CHAINE' 'Pdt min(conv, diff)' ;
  268. 'DESSIN' (evco 'ET' evdi 'ET' evmi) 'MIMA' 'TITR' 'Pas de temps' 'LEGE' tabt ;
  269. * Max |Fn de courant| = fn(t)
  270. evpsi = 'EVOL' 'MANU' 'TEMPS' hltps 'PSI' hlpsi ;
  271. 'DESSIN' evpsi 'MIMA' 'TITR' 'Max |Fn de courant|' ;
  272. * Max |Erreur Linf Vit.| = fn(t) (convergence)
  273. everr = 'EVOL' 'MANU' 'TEMPS' hltps 'ERRLINF' (('LOG' hlerr) '/' ('LOG' 10.0D0)) ;
  274. 'DESSIN' everr 'MIMA' 'TITR' 'Log Erreur Linf Vit.' ;
  275. 'FINPROC' ;
  276.  
  277. ** TRCHAMP
  278. 'DEBPROC' TRCHAMP ;
  279. 'ARGUMENT' trx*'TABLE   ' ;
  280.  as2 ama1  = 'KOPS' 'MATRIK' ;
  281.   trv = (trx . 'EQEX') ;
  282.   'SI' ('NON' ('MULT' (trv . 'PASDETPS' . 'NUPASDT') (trv . 'NPTRC'))) ;
  283.     'QUITTER' TRCHAMP ;
  284.   'FINSI' ;
  285. tinco = (trv . 'INCO') ;
  286. $dom = (trx . 'DOMZ') ;
  287. dom  = ('DOMA' $dom 'MAILLAGE') ;
  288. cdom = 'CONTOUR' dom ;
  289. tphy = (trv . 'PASDETPS' . 'TPS') ;
  290. * Vitesses
  291. vitesse = (tinco . 'UN') ;  vref = (tinco . 'UREF') ;
  292. titv = ('CHAINE' 'Champ de vitesse (VECT SOMMET) à t = ' tphy) ;
  293. vv = 'VECTEUR' vitesse (0.2D0 '/' vref) 'UX' 'UY' 'JAUNE' ;
  294. 'TRACER' vv cdom 'TITR' titv ;
  295. * Fonction de courant
  296. 'OPTION' 'ISOV' 'SULI' ;
  297. psi = (tinco . 'PSI') ;
  298. titpsi = ('CHAINE' 'Fonction de courant (SCAL SOMMET) à t = ' tphy) ;
  299. 'TRACER' psi dom cdom 14 'TITR' titpsi ;
  300. * Pression
  301. * NB : On n'utilise pas 'ELNO' : on dessine une pression
  302. * constante sur un element.
  303. pression = (trv . 'INCO' . 'PRESSION') ;
  304. modelp  = 'MODELISER' dom 'THERMIQUE' ;
  305. elemp   = 'KCHA' $dom 'CHAM' pression ;
  306. titpres = ('CHAINE' 'Pression (SCAL CENTRE) à t = ' tphy) ;
  307. 'TRACER' elemp modelp cdom 'TITR' titpres ;
  308. 'FINPROC' as2 ama1 ;
  309. *
  310. **  Fin de la definition des procedures
  311. **************************************************************
  312.  
  313. **************************************************************
  314. **  Debut du 'main'
  315. *
  316. **** Definition des constantes
  317. * l, h, nl, nh : donnees pour la boite.
  318. * egeom : erreur geometrique pour 'ELIMINATION'
  319. * aspect : rapport d'aspect de la boite
  320. * nu, g, rho : parametres du fluide
  321. * tg, td : temperatures sur les murs gauches et droits
  322. h = 1.0D0 ;
  323. 'SI' court ;  nh = 3 ; 'SINON' ;  nh = 5 ;  'FINSI' ;
  324. aspect = 4.0D0 ;
  325. l = h '*' aspect ;
  326. 'SI' court ;
  327.   nl = 'ENTIER' (nh '*' 3) ;
  328. 'SINON' ;
  329.   nl = 'ENTIER' (nh '*' aspect) ;
  330. 'FINSI' ;
  331. egeom = 1.0D-4 ;
  332. nu = 1.0D0 ; g = (0.D0 -1.D0) ;
  333. tg = 0.D0  ; td = aspect ;
  334.  
  335. **** Creation des maillages
  336. lp  = l '/' 2.0D0     ;  hp  = h '/' 2.0D0 ;
  337. mlp = (-1.0D0) '*' lp ;  mhp = (-1.0D0 '*' hp) ;
  338. pmil = lp hp ;
  339. pA = 0.0D0 0.0D0 ;  pB = l 0.0D0 ;
  340. pC = l h         ;  pD = 0.0D0 h ;
  341. pDA = 0.0D0 hp   ;  pBC = l hp ;
  342. bas    = pA 'DROIT' nl pB ;  haut  = pD 'DROIT' nl pC ;
  343. rwall  = pB 'DROIT' nh pC ;  lwall = pA 'DROIT' nh pD ;
  344. median = pDA 'DROIT' nl pBC ;
  345. mt = 'DALLER' haut rwall bas lwall ;
  346. mt = chan mt QUAF ;
  347. $mt= MODE mt 'NAVIER_STOKES' DISCR ; DOMA $MT 'IMPR' ;
  348. mt = (DOMA $mt   'MAILLAGE') ;  'TASSER' mt ;
  349. haut  = 'CHANGER' haut SEG2   ;  bas   = 'CHANGER' bas SEG2 ;
  350. lwall = 'CHANGER' lwall SEG2 ;   rwall = 'CHANGER' rwall SEG2 ;
  351. median =  'CHANGER' median SEG2 ;
  352. cmt = (lwall 'ET' bas 'ET' rwall 'ET' haut) ;
  353. 'ELIMINATION'  (mt 'ET' bas 'ET' rwall 'ET' haut  'ET' lwall 'ET' cmt 'ET' median) egeom ;
  354.  
  355. **** Initialisation de la table des historiques
  356. *
  357. thist = 'TABLE' 'HIST' ;
  358. TABTPS= 'TEMPS' 'NOEC';
  359. thist . 'tpscpu' = TABTPS.'TEMPS_CPU'.'INITIAL' ;
  360. * Listes
  361. thist . 'ltps'  = 'PROG' ;
  362. thist . 'lpsi'  = 'PROG' ;
  363. thist . 'lerr'  = 'PROG' ;
  364. thist . 'lpdt'  = 'TABLE' 'LPDT' ; thl = (thist . 'lpdt') ;
  365. thl . 'conv'  = 'PROG' ;  thl . 'diff' = 'PROG' ;
  366. thl . 'min'   = 'PROG' ;
  367.  
  368. **** Initialisation des tables rv et rvp
  369. rv = 'EQEX' $mt 'OPTI' discconv 'ZONE' $mt 'OPER' 'NS'   'NU' 'GB' 'TN' 'TREF' 'INCO' 'UN' 'OPTI' 'CENTREE' 'ZONE' $mt 'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN' ;
  370. * Conditions aux limites
  371. rv = 'EQEX' rv  'CLIM' 'UN' 'UIMP' cmt 0.0D0 'UN' 'VIMP' cmt 0.0D0 ;
  372. 'SI' graph ;
  373. rv = 'EQEX' rv  'ZONE' $mt 'OPER' 'CALCPSI' 'ZONE' $mt 'OPER' 'MAJHIST' 'ZONE' $mt 'OPER' 'TRCHAMP' ;
  374. 'FINSI' ;
  375.  
  376. RVP = EQEX 'OPTI' 'EF' KPRESS
  377.            'ZONE' $MT  'OPER'  'KBBT'  -1. betap  'INCO' 'UN' 'PRES' ;
  378. pfpres = 'POIN' ('DOMA' $MT KPRESS) 'PROC' pmil ;
  379. ppfpres = 'MANU' 'POI1' pfpres ;
  380. rvp = 'EQEX' rvp  'CLIM' 'PRES' 'TIMP' ppfpres 0.0D0 ;
  381.  
  382.  
  383.     rvp.'METHINV'.TYPINV=1 ;
  384.     rvp.'METHINV'.IMPINV=0 ;
  385.     rvp.'METHINV'.NITMAX=300;
  386.     rvp.'METHINV'.PRECOND=3 ;
  387.     rvp.'METHINV'.RESID  =1.e-8 ;
  388.     rvp.'METHINV' . 'FCPRECT'=100 ;
  389.     rvp.'METHINV' . 'FCPRECI'=100 ;
  390.  
  391.   RV.'PROJ'     =RVP ;
  392.  
  393. rv . 'INCO' = 'TABLE' 'INCO' ;
  394. rv . 'HISTOIRE' = thist ;
  395.  
  396. **** Initialisation des inconnues
  397. * V -> 0 ; T gradient constant avec tg sur lwall, td sur rwall
  398. rv . 'INCO' . 'UN-1' = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0  0.D0) ;
  399. rv . 'INCO' . 'UN'   = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0  0.D0) ;
  400. rv . 'INCO' . 'PRES' = 'KCHT' $mt 'SCAL'  KPRESS   0.D0        ;
  401. rv . 'INCO' . 'TN' = 'KCHT' $mt 'SCAL' 'SOMMET' (tg '+' ((('COORDONNEE' 1 mt) '/' l) '*' (td '-' tg)) ) ;
  402.  
  403. **** Valeurs des parametres :
  404. *
  405. ** Physiques
  406. *
  407. * Gr : Nombre de Grashof
  408. * tref : temperature de reference (approx. de Boussinesq)
  409. * Il faut la dissymetriser pour observer la transition
  410. * 3 -> 2 cellules
  411. * uref : echelle de vitesse (cf. adimensionnement)
  412. * ampvit : facteur d'amplification pour le trace des vitesses
  413. *
  414. 'SI' court ;
  415. Gr = 5.0D3 ;  tref = (tg '+' td) '/' 2.0D0 ;
  416. 'SINON' ;
  417. Gr = 40.0D3 ;  tref = tg ;
  418. 'FINSI' ;
  419. uref = (Gr '**' 0.5D0) ;
  420. *
  421. ** Numeriques
  422. *
  423. * nbiter : nombre d'iterations
  424. * nptrc : frequence pour le trace des champs (.ps ou à l'ecran)
  425. * fidt  : frequence d'affichage des messages à l'ecran
  426. * alpha : multiplicateur du pas de temps
  427. *         pour les elements MACRO
  428. 'SI' court ;
  429. nbiter = 250 ;  nptrc = 125 ; fidt = 50 ;
  430. 'SINON' ;
  431. nbiter = 15000 ;  nptrc = 2500 ; fidt  = 100 ;
  432. 'FINSI' ;
  433. alpha = 0.9D0 ;
  434. *
  435. ** Remplissage des tables
  436. *
  437. rv . 'INCO' . 'NU' = nu ;
  438. rv . 'INCO' . 'GB' = g '*' Gr ;
  439. rv . 'INCO' . 'TREF' = tref ;
  440. rv . 'INCO' . 'UREF' = uref ;
  441. rv . 'ALFA' = alpha ;
  442. rv . 'FIDT'   = fidt ;
  443. rv . 'NPTRC'  = nptrc ;
  444. rv . 'ITMA'   = nbiter ;
  445. *
  446. **
  447. ***
  448. **** Execution
  449. SUMMARY ;
  450. 'SI' inta ;  'OPTION' 'DONN' 5 ;  'FINSI' ;
  451. EXEC rv  ;
  452. TABTP1= 'TEMPS' 'NOEC';
  453. thist . 'tpscpu' = (thist . 'tpscpu') '+' TABTP1.'TEMPS_CPU'.'INITIAL' ;
  454. 'SI' graph ;  TRHIST thist ; 'FINSI' ;
  455. ****
  456. ***
  457. **
  458. *
  459. ** Le cas-test donne-t-il un resultat correct ??
  460. 'SI' court ;
  461. *
  462. ** Verification de la centrosymetrie sur les vitesse
  463. *
  464.   vit  = (rv . 'INCO' . 'UN') ;
  465.   tvit = vit 'TOURNER' 180.0D0 pmil ;
  466.   vit =  'KCHT' $mt 'VECT' 'SOMMET' vit ;
  467.   'ELIMINATION' (('EXTRAIRE' tvit 'MAIL') 'ET' mt) egeom ;
  468.   tvit = 'KCHT' $mt 'VECT' 'SOMMET' tvit ;
  469.   dv = 'KOPS' vit '-' tvit ;
  470.   critere1 = ('MAXIMUM' dv 'ABS') '/' uref ;
  471. * Valeur de critere1 au 2/10/97 : 1.718764099177966D-15
  472.   limite1 = 1.0D-12 ;
  473. *
  474. ** Evolution vitesse verticale le long de la mediane
  475. *
  476.   evvm = 'EVOL' 'CHPO' vit 'UY' median ;
  477.   evym = 'EXTRAIRE' evvm 'ORDO' ;
  478.   list (evym *1.e-6);
  479. *
  480. * Resultat obtenu le 2/10/97
  481. * correspondant à critere2 = 1.507288760336422D-16
  482.   evvref = 'PROG' 0.0D0 -.2205372973655706D+02 -.1537657415971809D+02  -.9886959157032784D+01 -.4994589224716669D+01 -.1325630316341117D+01   .5065207872228873D+00  .9297558065189220D+00 .5908058903982639D+00   .2945842774751989D-13 -.5908058903982122D+00 -.9297558065189337D+00  -.5065207872229142D+00  .1325630316341034D+01 .4994589224716784D+01   .9886959157032821D+01  .1537657415971808D+02 .2205372973655705D+02  0.0D0 ;
  483. * Resultat obtenu le 15/06/2014
  484.   evvref = 'PROG'  -7.24189E-39  -22.049  -15.375  -9.8821  -4.9918  -1.3199  0.50926  0.93264  0.59143     -1.04084E-12  -0.59143     -0.93264 -0.50926 1.3199   4.9918  9.8821       15.375       22.049      7.24188E-39 ;
  485.   critere2 = ('MAXIMUM' (evym '-' evvref) 'ABS') '/' uref ;
  486.   limite2 = 1.0D-4  ;
  487. * pour tenir compte de l'affichage de 5 chiffres significatifs
  488. * Tps CPU (le 2/10/97) ; 1009 centiemes de secondes
  489. * (court = VRAI ;  graph = FAUX ;  inta = FAUX ;)
  490.   TABTP1= 'TEMPS' 'NOEC';
  491.   thist . 'tpscpu'= (thist . 'tpscpu') '+' TABTP1.'TEMPS_CPU'.'INITIAL';
  492. *
  493. ** Affichage des messages et des erreurs eventuelles
  494. *
  495.   'SAUTER' 2 'LIGNE' ;
  496.   'MESSAGE' 'Debriefing :' ;
  497.   'MESSAGE' '------------' ;
  498.   'SAUTER' 'LIGNE' ;
  499.   'MESSAGE' ('CHAINE' 'Erreur/Limite de centrosymetrie du champ de vitesses :') ;
  500.   'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere1 ' / ' limite1) ;
  501.   'MESSAGE' ('CHAINE' 'Erreur/Limite Vy sur la mediane horizontale') ;
  502.   'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere2 ' / ' limite2) ;
  503.   'MESSAGE' ('CHAINE' 'Temps CPU / Temps CPU (30/9/97) : ') ;
  504.   'MESSAGE' ('CHAINE' ('ENTIER' (thist . 'tpscpu')) ' / ~1000' ' centieme(s) de secondes') ;
  505.    mess '      ->' critere1 '  'limite1' 'critere2'   'limite2   ;
  506.   'SI'  ('OU' (critere1 '>' limite1) (critere2 '>' limite2)) ;
  507.     'ERREUR' 5 ;
  508.   'FINSI' ;
  509.   'MESSAGE' 'Tout s_est bien passe' ;
  510.   'OPTION' 'ECHO' 1 ;
  511. 'FINSI' ;
  512. SUMMARY ;
  513. 'SI' inta ;  'OPTION' 'DONN' 5 ;  'FINSI' ;
  514.  'FIN' ;
  515. *
  516. **  Fin du 'main'
  517. **************************************************************
  518.  
  519.  
  520.  

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