Télécharger la plaquette Cast3M


fr en

Télécharger villers_platten.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  villers_platten.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ***************************************************************
  5. *                                                             *
  6. * NOM         : villers_platten.dgibi                         *
  7. * DESCRIPTION : fiche de validation CASTEM2000                *
  8. *               Mécanique des Fluides                         *
  9. *               Convection naturelle laminaire +              *
  10. *               convection thermocapillaire (Marangoni)       *
  11. *               en cavité rectangulaire (2D)                  *
  12. * FONCTIONS                                                   *
  13. *   TESTEES   : 'NS' algorithme semi-implicite                *
  14. *                    option SUPG                              *
  15. *                    approximation de Boussinesq              *
  16. *                    éléments QUA8 'MACRO'                    *
  17. *               'TSCA' option SUPG                            *
  18. *               'TOIMP' pour imposer la contrainte            *
  19. *                       tangentielle de Marangoni en surface  *
  20. *               'KOPS' 'GRAD' pour le calcul de cette         *
  21. *                      contrainte (éléments SEG2)             *
  22. *                                                             *
  23. * AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)      *
  24. *               e-mail : gounand@semt2.smts.cea.fr            *
  25. *                                                             *
  26. ***************************************************************
  27. *                                                             *
  28. * VERSION     : v1, 09/10/97, version initiale                *
  29. * HISTORIQUE  : v1, 09/10/97, création                        *
  30. * HISTORIQUE  :                                               *
  31. * HISTORIQUE  :                                               *
  32. * HISTORIQUE  :                                               *
  33. *                                                             *
  34. ***************************************************************
  35. *                                                             *
  36. * Priere de completer la partie HISTORIQUE apres modification *
  37. * du jeu de donnees afin de faciliter la maintenance          *
  38. *                                                             *
  39. ***************************************************************
  40. *
  41. * DESCRIPTION DETAILLEE :
  42. *
  43. *  ----  Convection naturelle laminaire + convection
  44. *        thermocapillaire (Marangoni) sur la surface libre
  45. *        (supposée plane).
  46. *                Algorithme semi implicite
  47. *                opérateur NS  (option SUPG)
  48. *
  49. *              phi = 0
  50. *        - - - - - - - - - - - - 
  51. *        |                      |
  52. * T1 = 0 | H                    | T2 = A = L/H
  53. *        |          L           |
  54. *        ------------------------
  55. *              phi = 0
  56. *
  57. * U = 0 sur les parois sauf sur la surface ou on applique
  58. * la condition de Levich :
  59. *      rho nu (dUx/dz) = (dsigma/dT) (dT/dx).
  60. *
  61. * La structure de l'écoulement ne dépend que des 4 nombres
  62. * adimensionnés suivants :
  63. *   - le rapport d'allongement A = L / H ;
  64. *   - le nombre de Prandtl  : Pr = nu / kappa ;
  65. *   - le nombre de Rayleigh :
  66. *          Ra = (g beta (T2 - T1) h**3) / (kappa nu)
  67. *             = Gr * (A Pr)        Gr : nombre de Grashof
  68. *   - le nombre de Marangoni :
  69. *          Ma = -( dsigma/dT (T2 - T1) h) / (rho alpha nu)
  70. *             = Ma2 * (A Pr)      Ma2 : autre définition du
  71. *                                       nb de Marangoni
  72. *************************************************************
  73. * Références : - Gounand, S. 1997
  74. *                Simulation Numérique de cellules
  75. *                thermoconvectives. 
  76. *                Rapport CEA/DMT 97/357.
  77. *              - Villers, D. & Platten, J.K. 1991
  78. *                Coupled buoyancy and Marangoni convection
  79. *                in acetone : experiments and comparison
  80. *                with numerical simulations
  81. *                J. Fluid Mech. Vol. 234 pp. 487-510. 
  82. *                                                             
  83. ***************************************************************
  84. *
  85. * JEU DE DONNEE :
  86. *
  87. prg = ('CHAINE' 'Villers & Platten') ;
  88. **** Options
  89. * court = VRAI  : le cas-test s'effectue en une quinzaine
  90. *                 de secondes.
  91. *                 Les paramètres sont les suivants :
  92. *                        A = 4  ;  Pr = 4 ;
  93. *                       Ra = 0  ;  Ma = 1000
  94. *                  => une cellule stationnaire
  95. *                 La comparaison s'effectue sur le profil
  96. *                 Ux(y) le long de la médiane verticale.
  97. *                 Il doit se rapprocher de la solution
  98. *                 obtenue pour une cavité infinie, lorsque
  99. *                 l'effet Maragoni est la seule cause du mvt :
  100. *              
  101. *                   Ux(z) proportionnel à z(z-2/3)   
  102. *
  103. * court = FAUX  : maillage raffiné ; on utilise les
  104. *                 paramètres suivants :
  105. *                        A = 12 ;  Pr = 4.24 ;
  106. *                    Ra = 4300  ;  Ma = 7600
  107. *                 configuration à une cellule + (1-2)
  108. *                 rouleaux périphériques.
  109. *                 La comparaison s'effectue sur le profil
  110. *                 Ux(y) le long de la médiane verticale.
  111. *                 Il doit se rapprocher des solutions
  112. *                 numériques et expérimentales obtenue par
  113. *                 Villers et Platten (voir références) 
  114. *                 environ 2h CPU avec les graphiques
  115. *
  116. * graph = VRAI  : on trace les graphiques (historiques, champs)
  117. *                 On a plus de calculs avec cette option
  118. *                 (fonction de courant, convergence...)
  119. *                 On mettra donc graph = faux pour la fiche
  120. *                 de validation
  121. *
  122. * inta  = VRAI si on est en interactif ;  0 sinon
  123. *
  124. * discconv : option de discrétisation des termes de
  125. *            convection pour les équations de Navier-Stokes
  126. *            choix possibles : 'CENTREE',  'SUPG' ou 'SUPGCC' 
  127. *
  128. court = VRAI ;
  129. graph = FAUX ;
  130. inta  = FAUX ;
  131. typelem = 'QUA8' ;
  132. discconv = 'SUPG' ;
  133. 'OPTION' 'DIME' 2 'ELEM' typelem ;
  134.  
  135. 'SI' ('NON' inta) ;
  136.   'OPTION' 'TRAC' 'PS' ;
  137.   'OPTION' 'ECHO' 0 ;
  138. 'SINON' ;
  139. **od = 'TEXTE' 'OPTION DONN 3' ;
  140. 'FINSI' ;
  141.  
  142. *************************************************************
  143. **** Définition de quelques procédures :
  144. **    -> SUMMARY : affiche un résumé des paramètres
  145. **    -> CALCMAR : calcul de la contrainte de Marangoni
  146. **                 (cisaillement) sur la surface.
  147. **                 (appelée par 'EXEC' via la table 'EQEX') 
  148. **    -> MAJHIST : mise à jour d'historiques
  149. **                 (appelée par 'EXEC' via la table 'EQEX') 
  150. **    -> TRHIST  : tracé des historiques
  151. **    -> TRCHAMP : tracé des différents champs incluant 
  152. **                 la fonction de courant (qui est calculee)
  153. **                 (appelée par 'EXEC' via la table 'EQEX') 
  154. ** NB : pour simplifier, on ne transmet pas les paramètres à
  155. **      SUMMARY. C'est de la mauvaise programmation...
  156. **      Les autres procédures sont correctes
  157. **      sous toutes réserves
  158. *
  159. ** SUMMARY
  160. 'DEBPROC' SUMMARY ;
  161. 'MESSAGE' ('CHAINE' 'Programme : ' prg) ;
  162. 'SAUTER' 'LIGNE' ;
  163. 'MESSAGE' ('CHAINE' '***  GEOMETRIE :') ;
  164. 'MESSAGE' ('CHAINE' 'Boite carrée : ' l ' x' h) ;
  165. 'MESSAGE' ('CHAINE' '      soient :  ' nl 'x' nh '='  (nl '*' nh) ' éléments ' typelem) ;
  166. 'MESSAGE' ('CHAINE' '      soient :  ' (nl '*' nh '*' 4) ' éléments linéaires') ;
  167. 'SAUTER' 'LIGNE' ;
  168. 'MESSAGE' ('CHAINE' '***  PHYSIQUE :') ;
  169. 'MESSAGE' ('CHAINE' 'Rapport d_aspect : ' ('ENTIER' aspect)) ;
  170. 'MESSAGE' ('CHAINE' 'Nombre de Prandtl   : ' Pr ) ;
  171. 'MESSAGE' ('CHAINE' 'Nombre de Rayleigh  : ' ('ENTIER' Ra) ) ;
  172. 'MESSAGE' ('CHAINE' 'Nombre de Marangoni : ' ('ENTIER' Ma) ) ;
  173. 'MESSAGE' ('CHAINE' 'Nombre de Grashof   : '  Gr ) ;
  174. 'MESSAGE' ('CHAINE' 'Nombre de Marangon2 : '  Ma2) ;
  175. 'MESSAGE' ('CHAINE' 'Température de référence : ' tref
  176.                     ' comprise entre 0 et '  ('ENTIER' aspect)) ;
  177.   'SAUTER' 'LIGNE' ;
  178. 'MESSAGE' ('CHAINE' '***  NUMERIQUE :') ;
  179. 'MESSAGE' ('CHAINE' 'Discrétisation ' discconv ' des termes de convection.') ;
  180. 'MESSAGE' ('CHAINE' ' Alfa  = ' alpha) ;
  181. 'MESSAGE' ('CHAINE' ' Béta  = ' beta) ;
  182. 'SAUTER' 'LIGNE' ;
  183. rvpdt = rv . 'PASDETPS' ;
  184. 'MESSAGE' ('CHAINE' '***  TEMPS :') ;
  185. 'MESSAGE' ('CHAINE' 'Nb. itérations : ' ((rvpdt . 'NUPASDT')
  186.                                                        '-' 1)) ;
  187. 'MESSAGE' ('CHAINE' 'Temps physique : ' (rvpdt . 'TPS')) ;
  188. 'MESSAGE' ('CHAINE' 'Temps CPU      : ' ('ENTIER' (thist . 'tpscpu')) ' centième(s) de seconde') ;
  189. 'FINPROC' ;
  190.  
  191. ** CALCMAR
  192. 'DEBPROC' CALCMAR ;
  193. 'ARGUMENT' trx*'TABLE   ' ;
  194. trv = (trx . 'EQEX') ;
  195. tinco = (trv . 'INCO') ;
  196. $surf = (trx . 'DOMZ') ;
  197. clevich = (-1.0D0 '*' (tinco . 'CLEV')) ;
  198. tn  = (rv . 'INCO' . 'TN') ;
  199. * On restreint le champ de température à la surface
  200. tn2 = 'KCHT' $surf 'SCAL' 'SOMMET' tn ;
  201. tinco . 'TSURF' =  'KCHT' $surf 'VECT' 'CENTRE' (('KOPS' tn2 'GRAD' $surf) '*' clevich) ;
  202. 'FINPROC' trx ;
  203.  
  204. ** MAJHIST
  205. 'DEBPROC' MAJHIST ;
  206. 'ARGUMENT' trx*'TABLE   ' ;
  207. trv = (trx . 'EQEX') ;
  208. tinco = (trv . 'INCO') ;
  209. th    = (trv . 'HISTOIRE') ;
  210. tp    = (trv . 'PASDETPS') ;
  211. vit    = (tinco . 'UN') ;  vitm1 = (tinco . 'UN-1') ;
  212. vitref = (tinco . 'UREF') ;
  213. npdt = (tp . 'NUPASDT') ;
  214. * Calcul de l'erreur relative (norme Linfini) entre deux
  215. * champs de vitesse successifs
  216. errrel = ('MAXIMUM' (vit '-' vitm1) 'ABS') '/' vitref ;
  217. 'SI' ('MULT' npdt (trv . 'FIDT')) ;
  218.   'MESSAGE' ('CHAINE' 'Pdt : ' npdt '  |Un - Un_1| Linf = ' errrel) ;
  219. 'FINSI' ;
  220. tinco . 'UN-1' = 'COPIER' vit ;
  221. * Mise à jour des différentes listes
  222. th . 'ltps'  = ((th . 'ltps') 'ET' ('PROG' (tp . 'TPS'))) ;
  223. th . 'lerr'  = ((th . 'lerr') 'ET' ('PROG' errrel)) ;
  224. thl = (th . 'lpdt') ;
  225. thl . 'conv' = ((thl . 'conv') 'ET' ('PROG' (tp . 'DTCONV'))) ;
  226. thl . 'diff' = ((thl . 'diff') 'ET' ('PROG' (tp . 'DTDIFU'))) ;
  227. thl . 'min'  = ((thl . 'min')  'ET' ('PROG' (tp . 'DELTAT-1'))) ;
  228. 'FINPROC' trx ;
  229.  
  230.  
  231. ** TRHIST
  232. 'DEBPROC' TRHIST ;
  233. 'ARGUMENT' th*'TABLE   ' ;
  234. * On enlève le premier indice des listes
  235. hltps = 'ENLEVER' (th . 'ltps') 1 ;
  236. hlerr = 'ENLEVER' (th . 'lerr') 1 ;
  237. thl = (th . 'lpdt') ;
  238. hlpdt  = 'ENLEVER' (thl . 'min')  1 ;
  239. hlconv = 'ENLEVER' (thl . 'conv') 1 ;
  240. hldiff = 'ENLEVER' (thl . 'diff') 1 ;
  241. * Et meme les deux premiers pour l'erreur sur les vitesses
  242. hltp2 = 'ENLEVER' hltps 1 ;
  243. hler2 = 'ENLEVER' hlerr 1 ;
  244. * Pas de temps = fn(t)
  245. evpdt = 'EVOL' 'MANU' 'TEMPS' hltps 'PDT' hlpdt ;
  246. 'DESSIN' evpdt 'MIMA' 'TITR' 'Pas de temps' ;
  247. * LOG10(pas de temps) = fn(t)
  248. ltco = ('LOG' (hlconv '*' alpha)) '/' ('LOG' 10.0D0) ;
  249. ltdi = ('LOG' (hldiff '*' alpha)) '/' ('LOG' 10.0D0) ;
  250. ltmi = ('LOG' hlpdt) '/' ('LOG' 10.0D0) ;
  251. evco = 'COULEUR' ('EVOL' 'MANU' 'TEMPS' hltps 'dtconv' ltco) 'BLEU' ;
  252. evdi = 'COULEUR' ('EVOL' 'MANU' 'TEMPS' hltps 'dtdiff' ltdi) 'ROUG' ;
  253. evmi = 'COULEUR' ('EVOL' 'MANU' 'TEMPS' hltps 'dt' ltmi) 'BLAN' ;
  254. tabt = 'TABLE' ;  tabt . 'TITRE' = 'TABLE' ;
  255. tabt . 1 = 'TIRC' ;
  256. tabt . 2 = 'TIRL' ;
  257. tabt . 3 = 'NOLI' ;
  258. tabt . 'TITRE' . 1 = 'CHAINE' 'Pdt convection' ;
  259. tabt . 'TITRE' . 2 = 'CHAINE' 'Pdt diffusion' ;
  260. tabt . 'TITRE' . 3 = 'CHAINE' 'Pdt min(conv, diff)' ;
  261. 'DESSIN' (evco 'ET' evdi 'ET' evmi) 'MIMA' 'TITR' 'Pas de temps' 'LEGE' tabt ;
  262. * Max |Erreur Linf Vit.| = fn(t) (convergence)
  263. everr = 'EVOL' 'MANU' 'TEMPS' hltp2 'ERRLINF' (('LOG' hler2) '/' ('LOG' 10.0D0)) ;
  264. 'DESSIN' everr 'MIMA' 'TITR' 'Log Erreur Linf Vit.' ;
  265. 'FINPROC' ;
  266.  
  267. ** TRCHAMP
  268. 'DEBPROC' TRCHAMP ;
  269. 'ARGUMENT' trx*'TABLE   ' ;
  270.   trv = (trx . 'EQEX') ;
  271.   'SI' ('NON' ('MULT' (trv . 'PASDETPS' . 'NUPASDT') (trv . 'NPTRC'))) ;
  272.     'QUITTER' TRCHAMP ;
  273.   'FINSI' ;
  274. tinco = (trv . 'INCO') ;
  275. $dom = (trx . 'DOMZ') ;
  276. dom  = ($dom . 'MAILLAGE') ;
  277. cdom = 'CONTOUR' dom ;
  278. ** On calcule tout d'abord la fonction de courant
  279.   vitesse = (tinco . 'UN') ;
  280.   sw = (-1.) * ('KOPS' vitesse 'ROT' $dom) ;
  281.   rk = 'EQEX' $dom 'OPTI' 'EF' 'IMPL' 'ZONE' $dom 'OPER' 'LAPN' 1.0D0   'INCO' 'PSI' 'ZONE' $dom 'OPER' 'FIMP' sw      'INCO' 'PSI' 'CLIM' 'PSI' 'TIMP' cdom 0.0D0 ;
  282.   rk . 'INCO' = 'TABLE' 'INCO' ;
  283.   rk . 'INCO' . 'PSI' = 'KCHT' $dom 'SCAL' 'SOMMET' 0.0D0 ;
  284. *-*-*-*- A remplacer par EXEC qd ca marchera...
  285.   EXIC rk ;
  286. *-*-*-*-
  287.   psi = (rk . 'INCO' . 'PSI') ;
  288. ** On effectue une affinité si la cavité est trop allongée
  289. augmen = ((tinco . 'ASPECT') '>' 6.0D0) ;
  290. 'SI' augmen ;
  291.   dy = 'NOMC' 'UY'  (('COORDONNEE' 2 ($dom . 'SOMMET')) 'ET'  ('COORDONNEE' 2 ($dom . 'CENTRE'))) 'NATU' 'DISCRET' ;
  292.   dx = 'NOMC' 'UX' (dy '*' 0.0D0)  'NATU' 'DISCRET' ;
  293.   dv = 2.0D0 '*' (dx 'ET' dy) ;
  294.   svform = 'FORME' ;
  295.   'FORME' dv ;
  296. 'FINSI' ;
  297. ** Puis on trace les différents champs :
  298. tphy = (trv . 'PASDETPS' . 'TPS') ;
  299. * Vitesses
  300. vref = (tinco . 'UREF') ;
  301. titv = ('CHAINE' 'Champ de vitesse (VECT SOMMET) à t = '
  302.                  tphy) ;
  303. vv = 'VECTEUR' vitesse (2.0D0 '/' vref) 'UX' 'UY' 'JAUNE' ;
  304. 'TRACER' vv cdom 'TITR' titv ;
  305. * Fonction de courant
  306. 'OPTION' 'ISOV' 'SULI' ;
  307. titpsi = ('CHAINE' 'Fonction de courant (SCAL SOMMET) à t = '
  308.                  tphy) ;
  309. 'TRACER' psi dom cdom 14 'TITR' titpsi ;
  310. * température
  311. tn = (tinco . 'TN') ;
  312. tittn  = ('CHAINE' 'Champ de température (SCAL SOMMET) à t = '
  313.                  tphy) ;
  314. 'TRACER' tn  dom cdom 14 'TITR' tittn ;
  315. * Pression
  316. * NB : On n'utilise pas 'ELNO' : on dessine une pression
  317. * constante sur un élément.
  318. pression = (trv . 'PRESSION' . 'PRESSION') ;
  319. modelp  = 'MODELISER' dom 'THERMIQUE' ;
  320. elemp   = 'KCHA' $dom 'CHAM' pression ;
  321. titpres = ('CHAINE' 'Pression (SCAL CENTRE) à t = '
  322.                  tphy) ;
  323. 'TRACER' elemp modelp cdom 'TITR' titpres ;
  324. ** On rétablit le rapport d'aspect initial
  325. 'SI' augmen ;  'FORME' svform ; 'FINSI' ;
  326. 'FINPROC' ;
  327. *                 
  328. **  Fin de la définition des procédures
  329. **************************************************************
  330.  
  331. **************************************************************
  332. **  Début du 'main'
  333. *
  334. **** Définition des constantes
  335. * l, h, nl, nh : données pour la boite.
  336. * egeom : erreur géométrique pour 'ELIMINATION'
  337. * aspect : rapport d'aspect de la boite
  338. * nu, g  : paramètres du fluide
  339. * tg, td : températures sur les murs gauches et droits
  340. h = 1.0D0 ;
  341. 'SI' court ;
  342.   nh = 3 ;   aspect = 4.0D0 ;
  343. 'SINON' ;
  344.   nh = 15 ;  aspect = 12.0D0 ;
  345. 'FINSI' ;
  346. nl = 'ENTIER' (nh '*' 4) ;
  347. l = h '*' aspect ;
  348. egeom = 1.0D-4 ;
  349. g = (0.D0 -1.D0) ;
  350. tg = 0.D0  ; td = aspect ;
  351.  
  352. **** Création des maillages 
  353. lp  = l '/' 2.0D0     ;  hp  = h '/' 2.0D0 ;
  354. mlp = (-1.0D0) '*' lp ;  mhp = (-1.0D0 '*' hp) ;
  355. pmil = lp hp ;
  356. pA = 0.0D0 0.0D0 ;  pB  = l 0.0D0 ;
  357. pC = l h         ;  pD  = 0.0D0 h ;
  358. pAB = lp 0.0D0   ;  pCD = lp h ;
  359. bas     = pA  'DROIT' nl pB ;  haut  = pD 'DROIT' nl pC ;
  360. rwall   = pB  'DROIT' nh pC ;  lwall = pA 'DROIT' nh pD ;
  361. mediany = pAB 'DROIT' nh pCD ;
  362. mt = 'DALLER' haut rwall bas lwall ;
  363. $mt = 'DOMA' mt 'MACRO' 'IMPR' ;
  364. mt = ($mt . 'MAILLAGE') ;  'TASSER' mt ;
  365. haut  = 'CHANGER' haut SEG2   ;  bas   = 'CHANGER' bas SEG2 ;
  366. lwall = 'CHANGER' lwall SEG2 ;   rwall = 'CHANGER' rwall SEG2 ;
  367. mediany =  'CHANGER' mediany SEG2 ;
  368. enleelem = haut 'ELEM' 1 ;
  369. haut2 = 'DIFF' haut enleelem ;
  370. cmt = (lwall 'ET' bas 'ET' rwall 'ET' haut) ;
  371. $haut  = 'DOMA' haut  'INCL' $mt egeom ;
  372. $haut2 = 'DOMA' haut2 'INCL' $mt egeom ;
  373. $bas   = 'DOMA' bas   'INCL' $mt egeom ;
  374. 'ELIMINATION'  (mt 'ET' bas 'ET' rwall 'ET' haut  'ET' haut2 'ET' lwall 'ET' cmt 'ET' mediany) egeom ;
  375.  
  376. **** Initialisation de la table des historiques
  377. * 
  378. thist = 'TABLE' 'HIST' ;
  379. TABTPS = TEMP 'NOEC';
  380. thist . 'tpscpu' = TABTPS.'TEMPS_CPU'.'INITIAL';
  381. * Listes
  382. thist . 'ltps'  = 'PROG' ;
  383. thist . 'lerr'  = 'PROG' ;
  384. thist . 'lpdt'  = 'TABLE' 'LPDT' ; thl = (thist . 'lpdt') ;
  385. thl . 'conv'  = 'PROG' ;  thl . 'diff' = 'PROG' ;
  386. thl . 'min'   = 'PROG' ;
  387.  
  388. **** Initialisation des tables rv et rvp 
  389. rv = 'EQEX' $mt 'OPTI' discconv 'ZONE' $mt 'OPER' 'NS'   'CLPU' 'CBOU' 'TN' 'TREF' 'INCO' 'UN' 'OPTI' discconv 'ZONE' $mt 'OPER' 'TSCA' 'CLPT'  'UN' 0.0D0 'INCO' 'TN' ;
  390. rv = 'EQEX' rv  'ZONE' $haut 'OPER' 'CALCMAR'   'ZONE' $haut 'OPER' 'TOIMP' 'TSURF' 'INCO' 'UN' ;
  391. tbas = 'KCHT' $bas 'SCAL' 'SOMMET' (tg '+' ((('COORDONNEE' 1 bas) '/' l) '*' (td '-' tg)) ) ;
  392. * Conditions aux limites :
  393. * traitement special pour pA et pB communs respectivement
  394. * à lwall-bas et a bas-rwall sinon les temperatures valent
  395. * deux fois la valeur qu'on veut leur donner
  396. * On utilise 'MANUEL' car 'CLIM' ne prend que les objets de
  397. * type 'MAILLAGE'
  398. *
  399. rv = 'EQEX' rv 'CLIM' 'UN' 'UIMP' (lwall 'ET' rwall 'ET' bas) 0.0D0 'CLIM' 'UN' 'VIMP' (lwall 'ET' rwall 'ET' haut 'ET' bas) 0.0D0 'CLIM' 'TN' 'TIMP' lwall tg 'CLIM' 'TN' 'TIMP' ('MANUEL' 'POI1' pA)  (-1.0D0 '*' tg) 'CLIM' 'TN' 'TIMP' rwall td 'CLIM' 'TN' 'TIMP' ('MANUEL' 'POI1' pB)  (-1.0D0 '*' td) 'CLIM' 'TN' 'TIMP' bas   tbas ;
  400. 'SI' graph ;
  401. rv = 'EQEX' rv  'ZONE' $mt 'OPER' 'MAJHIST' 'ZONE' $mt 'OPER' 'TRCHAMP' ;
  402. 'FINSI' ;
  403. rvp = 'EQPR' $mt 'ZONE' $mt 'OPER' 'PRESSION' 0.D0 'PIMP' 0.D0 ;
  404. rv . 'INCO' = 'TABLE' 'INCO' ;
  405. rv . 'PRESSION' = rvp ;
  406. rv . 'HISTOIRE' = thist ;
  407.  
  408. **** Initialisation des inconnues
  409. * V -> 0 ; T gradient constant avec tg sur lwall, td sur rwall
  410. rv . 'INCO' . 'UN-1' = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0  0.D0) ;
  411. rv . 'INCO' . 'UN'   = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0  0.D0) ;
  412. rv . 'INCO' . 'TN' = 'KCHT' $mt 'SCAL' 'SOMMET' (tg '+' ((('COORDONNEE' 1 mt) '/' l) '*' (td '-' tg)) ) ;
  413.  
  414. **** Valeurs des paramètres :
  415. *
  416. ** Physiques
  417. *
  418. * Pr : Nombre de Prandtl
  419. * Ra : Nombre de Rayleign (Gr : Grashof)
  420. * Ma : Nombre de Marangoni (Ma2 : autre définition)
  421. * tref : température de référence (approx. de Boussinesq)
  422. * uref : échelle de vitesse (cf. adimensionnement)
  423. * ampvit : facteur d'amplification pour le tracé des vitesses
  424. 'SI' court ;
  425. Pr = 4.0D0 ; Ra = 0.0D0 ;  Ma = 1000.0D0 ;
  426. 'SINON' ;
  427. Pr = 4.24D0 ; Ra = 4300.0D0 ;  Ma = 7600.0D0 ;
  428. 'FINSI' ;
  429. Gr = Ra '/' (Pr '*' aspect) ;  Ma2 = Ma '/' (Pr '*' aspect) ;
  430. sMa2 = (Ma2 '**' 0.5D0) ;
  431. clpu = 1.0D0 '/' sMa2  ;
  432. cbou = Gr '/' Ma2 ;
  433. clpt = 1.0D0 '/' (sMa2 '*' Pr) ;
  434. clev = 1.0D0 ;
  435. uref = sMa2 ;
  436. tref = (tg + td) / 2.0D0 ;
  437. ampvit = 2.0D0 '/' uref ;
  438. *
  439. ** Numériques
  440. *
  441. * nbiter : nombre d'itérations
  442. * nptrc : fréquence pour le tracé des champs (.ps ou à l'écran)
  443. * fidt  : fréquence d'affichage des messages à l'écran
  444. * alpha : multiplicateur du pas de temps
  445. * beta  : paramètre de stabilisation de la pression
  446. *         pour les éléments MACRO
  447. 'SI' court ;
  448. nbiter = 250 ;  nptrc = 125 ; fidt = 10 ;
  449. alpha = 0.9D0 ;  beta = 1.0D0 ;
  450. 'SINON' ;
  451. *nbiter = 10 ;  nptrc = 500 ; fidt  = 100 ;
  452. nbiter = 10000 ;  nptrc = 2000 ; fidt  = 100 ;
  453. alpha = 0.5D0 ;  beta = 1.0D0 ;
  454. 'FINSI' ;
  455. *
  456. ** Remplissage des tables
  457. *
  458. rv . 'INCO' . 'ASPECT' = aspect ;
  459. rv . 'INCO' . 'CLPU' = clpu ;
  460. rv . 'INCO' . 'CBOU' = g '*' cbou ;
  461. rv . 'INCO' . 'TREF' = tref ;
  462. rv . 'INCO' . 'UREF' = uref ;
  463. rv . 'INCO' . 'CLPT' = clpt ;
  464. rv . 'INCO' . 'CLEV' = clev ;
  465. rv . 'ALFA' = alpha ;
  466. rv . 'PRESSION' . 'BETA' = beta ;
  467. rv . 'FIDT'   = fidt ;
  468. rv . 'NPTRC'  = nptrc ;
  469. rv . 'ITMA'   = nbiter ;
  470. *
  471. **
  472. ***
  473. **** Exécution
  474. SUMMARY ;
  475. 'SI' inta ;  'OPTION' 'DONN' 5 ;  'FINSI' ;
  476. EXEC rv  ;
  477. TABTPS = TEMP 'NOEC';
  478. thist . 'tpscpu' = (thist . 'tpscpu') '+' TABTPS.'TEMPS_CPU'.'INITIAL';
  479. 'SI' graph ;  TRHIST thist ; 'FINSI' ;
  480. ****
  481. ***
  482. **
  483. *
  484. vit  = (rv . 'INCO' . 'UN') ;
  485. evvm = 'EVOL' 'BLEU' 'CHPO' vit 'UX' mediany ;
  486. ** Le cas-test donne-t-il un résultat correct ??
  487. 'OPTION' 'ECHO' 0 ;
  488. 'SI' court ;
  489.   'OPTION' 'ECHO' 1 ;
  490. *
  491. ** Vérification :
  492. ** La composante horizontale de la vitesse le long de la
  493. ** médiane verticale doit etre prop. à z(z-2/3)
  494. *
  495.   cfpoly evpoly = @POMI evvm 2 ;
  496.   cfprop = (cfpoly . 2) ;
  497.   cflin  = (cfpoly . 1) '/' cfprop ;
  498.   cfcte  = (cfpoly . 0) '/' cfprop ;
  499.   cflana = -2.0D0 '/' 3.0D0 ;
  500. * Résultat obtenu le 14/10/97 :
  501.   cfpref = -.3128910770680292D+01 ;
  502.   critere1 = 'ABS' ((cfprop '-' cfpref) '/' cfpref) ;
  503.   critere2 = 'ABS' ((cflin  '-' cflana) '/' cflana) ;
  504.   critere3 = 'ABS' cfcte ;
  505. * Résultats obtenus le 14/10/97 :
  506. * cflref = -.6784204323293079D+00 ;
  507. * cfcref =  .5876882831320513D-02 ;
  508. * correspondant à :
  509. * critere1 = 2.83861D-16 
  510. * critere2 = 1.76306D-02
  511. * critere3 = 5.87688D-03
  512.   limite1 = 1.D-07 ;
  513.   limite2 = 2.D-02 ;
  514.   limite3 = 1.D-02 ;
  515.   'SI' graph ;
  516.     tabg = 'TABLE' ;  tabg . 'TITRE' = 'TABLE' ;
  517.     tabg . 1 = 'NOLI MARQ CROI' ;
  518.     tabg . 'TITRE' . 1 = 'CHAINE' 'Résultat Castem2000' ;
  519.     tabg . 'TITRE' . 2 = 'CHAINE' 'Parabole théorique' ;
  520.     evpoly = ('COULEUR' evpoly 'JAUN') ;
  521.     titg = 'CHAINE' 'Ux(z) sur la médiane verticale' ;
  522.     'DESSIN' (evvm 'ET' evpoly) 'MIMA' 'TITR' titg 'LEGE' tabg ;
  523.   'FINSI' ;
  524. * Tps CPU (le 10/10/97) ; 1797 centièmes de secondes
  525. * (court = VRAI ;  graph = FAUX ;  inta = FAUX ;)
  526.   TABTPS = TEMP 'NOEC';
  527.   thist .'tpscpu' = (thist . 'tpscpu') '+' TABTPS.'TEMPS_CPU'.'INITIAL';
  528. *
  529. ** Affichage des messages et des erreurs éventuelles
  530. *
  531.   'OPTION' 'ECHO' 0 ;
  532.   'SAUTER' 2 'LIGNE' ;
  533.   'MESSAGE' 'Debriefing :' ;
  534.   'MESSAGE' '------------' ;
  535.   'MESSAGE' 'Comparaison de ux(z) sur la médiane verticale à : ' ;
  536.   'MESSAGE' 'alpha (z**2 - 2/3 z + 0)' ;
  537.   'SAUTER' 'LIGNE' ;
  538.   'MESSAGE' ('CHAINE' 'Erreur/Limite sur le coefficient alpha :') ;
  539.   'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere1 ' / ' limite1) ;
  540.   'MESSAGE' ('CHAINE' 'Erreur/Limite sur le coefficient de z (2/3)') ;
  541.   'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere2 ' / ' limite2) ;
  542.   'MESSAGE' ('CHAINE' 'Erreur/Limite sur le coefficient constant (0)') ;
  543.   'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere3 ' / ' limite3) ;
  544.   'MESSAGE' ('CHAINE' 'Temps CPU / Temps CPU (10/10/97) : ') ;
  545.   'MESSAGE' ('CHAINE' ('ENTIER' (thist . 'tpscpu')) ' / ~1800' ' centième(s) de secondes') ;
  546.   'SI' ('OU' ('OU' (critere1 '>' limite1) (critere2 '>' limite2))  (critere3 '>' limite3)) ;
  547.     'ERREUR' 5 ;
  548.   'FINSI' ;
  549.   'MESSAGE' 'Tout s_est bien passé' ;
  550. 'SINON' ;
  551.   'OPTION' 'ECHO' 1 ;
  552. **
  553. ** Vérification pour le cas long :
  554. ** on compare ux(z) sur la médiane avec des résultats
  555. ** experimentaux (attention, l'expérience n'est pas 2D =>
  556. ** vitesse plus faible a cause des parois latérales)
  557. ** et une autre simulation numérique (effectuée avec les
  558. ** memes parametres que la notre)..
  559. ** Voir l'article de Villers et Platten pour plus de précisions
  560. **
  561.   zexp  = 'PROG' 0.0D0 0.1D0  0.2D0  0.35D0 0.57D0 0.75D0 0.97D0 1.2D0  1.5D0  1.8D0  2.0D0  2.2D0  2.5D0 ;
  562.   uxexp = 'PROG' 0.0D0 0.6D0  0.75D0 0.95D0 1.05D0 0.95D0 0.8D0 0.4D0  -0.4D0 -1.0D0 -2.0D0 -2.5D0 -3.6D0 ;
  563.   zpas = (2.5D0 '/' 32) ;
  564.   zsvp  = 'PROG' 0.0D0 'PAS' zpas 'NPAS' 32 ;
  565.   uxsvp = 'PROG' 0.0D0  0.22D0 0.4D0   0.6D0  0.75D0 0.86D0 0.98D0 1.05D0 1.13D0  1.15D0 1.16D0 1.15D0 1.13D0 1.08D0  1.02D0 0.92D0  0.81D0 0.7D0 0.55D0  0.4D0  0.2D0 0.01D0  -0.2D0 -0.43D0 -0.7D0 -0.95D0 -1.2D0 -1.5D0 -1.8D0 -2.15D0 -2.5D0 -2.85D0 -3.25D0 ;
  566. *
  567. * On remet l'évolution obtenu avec C2000 en dimensionné
  568. * 
  569. *          mu à 20°C      densité
  570.   nuacet = 0.3265D-3  '/' 870.0D0 ;
  571.   lref = 2.5D0 ;
  572.   uref = ((nuacet '*' 1.0D6) '/' lref) '*' sMa2 ;
  573.   zdim  = ('EXTRAIRE' evvm 'ABSC') '*' lref ;
  574.   uxdim = ('EXTRAIRE' evvm 'ORDO') '*' uref ;
  575. * Comparaison des profils :
  576. *  on les remet sur zdim par interpolation linéaire
  577.   uxexpi = 'IPOL' zdim zexp uxexp ;
  578.   uxvpi  = 'IPOL' zdim zsvp uxsvp ;
  579.   lmux = 'PROG' ('MAXIMUM' uxdim 'ABS') ('MAXIMUM' uxexpi 'ABS') ('MAXIMUM' uxvpi 'ABS') ;
  580.   echux = 'MINIMUM' lmux ;
  581.   critere1 = ('MAXIMUM' (uxdim '-' uxvpi) 'ABS')  '/' echux ;
  582.   critere2 = ('MAXIMUM' (uxdim '-' uxexpi) 'ABS') '/' echux ;
  583. * Résultats obtenus le 15/10/97 :
  584. * critere1 = .7530548488495176D-01
  585. * critere2 = .3228687516012702D+00
  586.   limite1 = 0.08D0 ;
  587.   limite2 = 0.35D0 ;
  588. * trace eventuel
  589.   'SI' graph ;
  590.     evexp  = 'EVOL' 'ROSE' 'MANU' 'Z' zexp 'UX' uxexp ;
  591.     evsvp  = 'EVOL' 'JAUN' 'MANU' 'Z' zsvp 'UX' uxsvp ;
  592.     evexpi = 'EVOL' 'ROUG' 'MANU' 'Z' zdim 'UX' uxexpi ;
  593.     evvpi  = 'EVOL' 'VERT' 'MANU' 'Z' zdim 'UX' uxvpi ;
  594.     evdim  = 'EVOL' 'TURQ' 'MANU' 'Z' zdim 'UX' uxdim ;
  595.     titd = 'CHAINE' 'Profil Ux(z) sur la médiane verticale' ;
  596.     tabd = 'TABLE' ;  tabd . 'TITRE' = 'TABLE' ;
  597.     tabd . 1 = 'MARQ PLUS TIRR TITR Experience' ;
  598.     tabd . 2 = 'MARQ CROI TIRC TITR Vp_num' ;
  599.     tabd . 3 = 'MARQ PLUS NOLI TITR Exp_ipol' ;
  600.     tabd . 4 = 'MARQ CROI NOLI TITR Vp_num_ipol' ;
  601.     tabd . 5 = 'MARQ ETOI NOLI TITR K2000' ;
  602.     tabd . 'TITRE' . 1 = 'CHAINE' 'Expérience' ;
  603.     tabd . 'TITRE' . 2 = 'CHAINE' 'Numérique (V&P)' ;
  604.     tabd . 'TITRE' . 3 = 'CHAINE' 'Exp. (interpolé)' ;
  605.     tabd . 'TITRE' . 4 = 'CHAINE' 'V&P  (interpolé)' ;
  606.     tabd . 'TITRE' . 5 = 'CHAINE' 'Castem 2000' ;
  607.     'DESSIN' (evexp 'ET' evsvp 'ET' evexpi 'ET' evvpi 'ET' evdim) 'LEGE' tabd 'TITR' titd ;
  608.   'FINSI' ;
  609. * Tps CPU (le 10/10/97) ; 7043 secondes
  610. * (court = FAUX ;  graph = VRAI ;  inta = FAUX ;)
  611.   TABTPS = TEMP 'NOEC';
  612.   thist .'tpscpu' = (thist . 'tpscpu') '+' TABTPS.'TEMPS_CPU'.'INITIAL';
  613. *
  614. ** Affichage des messages et des erreurs éventuelles
  615. *
  616.   'OPTION' 'ECHO' 0 ;
  617.   'SAUTER' 2 'LIGNE' ;
  618.   'MESSAGE' 'Debriefing :' ;
  619.   'MESSAGE' '------------' ;
  620.   'MESSAGE' 'Comparaison de ux(z) sur la médiane verticale' ;
  621.   'SAUTER' 'LIGNE' ;
  622.   'MESSAGE' 'par rapport à la simulation numérique de V&P' ;
  623.   'MESSAGE' 'Erreur / Limite :' ;
  624.   'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere1 ' / ' limite1) ;
  625.   'SAUTER' 'LIGNE' ;
  626.   'MESSAGE' 'par rapport à l_expérience de V&P' ;
  627.   'MESSAGE' 'Erreur / Limite :' ;
  628.   'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere2 ' / ' limite2) ;
  629.   'SAUTER' 'LIGNE' ;
  630.   'MESSAGE' ('CHAINE' 'Temps CPU / Temps CPU (15/10/97) (moy2 queue) : ') ;
  631.   'MESSAGE' ('CHAINE' ('ENTIER' ((thist . 'tpscpu') '/' 100.0D0)) ' / ~7040' ' secondes') ;
  632.   'SI' ('OU' (critere1 '>' limite1) (critere2 '>' limite2)) ;
  633.     'ERREUR' 5 ;
  634.   'FINSI' ;
  635.   'MESSAGE' 'Tout s_est bien passé' ;
  636. 'FINSI' ;
  637. SUMMARY ;
  638. 'OPTION' 'ECHO' 1 ;
  639. 'SI' inta ;  'OPTION' 'DONN' 5 ;  'FINSI' ;
  640. 'FIN' ;
  641. *
  642. **  Fin du 'main'
  643. **************************************************************
  644.  
  645.  
  646.  
  647.  
  648.  
  649.  
  650.  
  651.  
  652.  
  653.  
  654.  
  655.  
  656.  
  657.  
  658.  

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