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 '='
  166. (nl '*' nh) ' éléments ' typelem) ;
  167. 'MESSAGE' ('CHAINE' ' soient : ' (nl '*' nh '*' 4)
  168. ' éléments linéaires') ;
  169. 'SAUTER' 'LIGNE' ;
  170. 'MESSAGE' ('CHAINE' '*** PHYSIQUE :') ;
  171. 'MESSAGE' ('CHAINE' 'Rapport d_aspect : ' ('ENTIER' aspect)) ;
  172. 'MESSAGE' ('CHAINE' 'Nombre de Prandtl : ' Pr ) ;
  173. 'MESSAGE' ('CHAINE' 'Nombre de Rayleigh : ' ('ENTIER' Ra) ) ;
  174. 'MESSAGE' ('CHAINE' 'Nombre de Marangoni : ' ('ENTIER' Ma) ) ;
  175. 'MESSAGE' ('CHAINE' 'Nombre de Grashof : ' Gr ) ;
  176. 'MESSAGE' ('CHAINE' 'Nombre de Marangon2 : ' Ma2) ;
  177. 'MESSAGE' ('CHAINE' 'Température de référence : ' tref
  178. ' comprise entre 0 et ' ('ENTIER' aspect)) ;
  179. 'SAUTER' 'LIGNE' ;
  180. 'MESSAGE' ('CHAINE' '*** NUMERIQUE :') ;
  181. 'MESSAGE' ('CHAINE' 'Discrétisation ' discconv ' des termes de
  182. convection.') ;
  183. 'MESSAGE' ('CHAINE' ' Alfa = ' alpha) ;
  184. 'MESSAGE' ('CHAINE' ' Béta = ' beta) ;
  185. 'SAUTER' 'LIGNE' ;
  186. rvpdt = rv . 'PASDETPS' ;
  187. 'MESSAGE' ('CHAINE' '*** TEMPS :') ;
  188. 'MESSAGE' ('CHAINE' 'Nb. itérations : ' ((rvpdt . 'NUPASDT')
  189. '-' 1)) ;
  190. 'MESSAGE' ('CHAINE' 'Temps physique : ' (rvpdt . 'TPS')) ;
  191. 'MESSAGE' ('CHAINE' 'Temps CPU : '
  192. ('ENTIER' (thist . 'tpscpu')) ' centième(s) de seconde') ;
  193. 'FINPROC' ;
  194.  
  195. ** CALCMAR
  196. 'DEBPROC' CALCMAR ;
  197. 'ARGUMENT' trx*'TABLE ' ;
  198. trv = (trx . 'EQEX') ;
  199. tinco = (trv . 'INCO') ;
  200. $surf = (trx . 'DOMZ') ;
  201. clevich = (-1.0D0 '*' (tinco . 'CLEV')) ;
  202. tn = (rv . 'INCO' . 'TN') ;
  203. * On restreint le champ de température à la surface
  204. tn2 = 'KCHT' $surf 'SCAL' 'SOMMET' tn ;
  205. tinco . 'TSURF' = 'KCHT' $surf 'VECT' 'CENTRE'
  206. (('KOPS' tn2 'GRAD' $surf)
  207. '*' clevich) ;
  208. 'FINPROC' trx ;
  209.  
  210. ** MAJHIST
  211. 'DEBPROC' MAJHIST ;
  212. 'ARGUMENT' trx*'TABLE ' ;
  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. npdt = (tp . 'NUPASDT') ;
  220. * Calcul de l'erreur relative (norme Linfini) entre deux
  221. * champs de vitesse successifs
  222. errrel = ('MAXIMUM' (vit '-' vitm1) 'ABS') '/' vitref ;
  223. 'SI' ('MULT' npdt (trv . 'FIDT')) ;
  224. 'MESSAGE' ('CHAINE' 'Pdt : ' npdt
  225. ' |Un - Un_1| Linf = ' errrel) ;
  226. 'FINSI' ;
  227. tinco . 'UN-1' = 'COPIER' vit ;
  228. * Mise à jour des différentes listes
  229. th . 'ltps' = ((th . 'ltps') 'ET' ('PROG' (tp . 'TPS'))) ;
  230. th . 'lerr' = ((th . 'lerr') 'ET' ('PROG' errrel)) ;
  231. thl = (th . 'lpdt') ;
  232. thl . 'conv' = ((thl . 'conv') 'ET' ('PROG' (tp . 'DTCONV'))) ;
  233. thl . 'diff' = ((thl . 'diff') 'ET' ('PROG' (tp . 'DTDIFU'))) ;
  234. thl . 'min' = ((thl . 'min') 'ET' ('PROG' (tp . 'DELTAT-1'))) ;
  235. 'FINPROC' trx ;
  236.  
  237.  
  238. ** TRHIST
  239. 'DEBPROC' TRHIST ;
  240. 'ARGUMENT' th*'TABLE ' ;
  241. * On enlève le premier indice des listes
  242. hltps = 'ENLEVER' (th . 'ltps') 1 ;
  243. hlerr = 'ENLEVER' (th . 'lerr') 1 ;
  244. thl = (th . 'lpdt') ;
  245. hlpdt = 'ENLEVER' (thl . 'min') 1 ;
  246. hlconv = 'ENLEVER' (thl . 'conv') 1 ;
  247. hldiff = 'ENLEVER' (thl . 'diff') 1 ;
  248. * Et meme les deux premiers pour l'erreur sur les vitesses
  249. hltp2 = 'ENLEVER' hltps 1 ;
  250. hler2 = 'ENLEVER' hlerr 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'
  259. ('EVOL' 'MANU' 'TEMPS' hltps 'dtconv' ltco) 'BLEU' ;
  260. evdi = 'COULEUR'
  261. ('EVOL' 'MANU' 'TEMPS' hltps 'dtdiff' ltdi) 'ROUG' ;
  262. evmi = 'COULEUR'
  263. ('EVOL' 'MANU' 'TEMPS' hltps 'dt' ltmi) 'BLAN' ;
  264. tabt = 'TABLE' ; tabt . 'TITRE' = 'TABLE' ;
  265. tabt . 1 = 'TIRC' ;
  266. tabt . 2 = 'TIRL' ;
  267. tabt . 3 = 'NOLI' ;
  268. tabt . 'TITRE' . 1 = 'CHAINE' 'Pdt convection' ;
  269. tabt . 'TITRE' . 2 = 'CHAINE' 'Pdt diffusion' ;
  270. tabt . 'TITRE' . 3 = 'CHAINE' 'Pdt min(conv, diff)' ;
  271. 'DESSIN' (evco 'ET' evdi 'ET' evmi)
  272. 'MIMA' 'TITR' 'Pas de temps'
  273. 'LEGE' tabt ;
  274. * Max |Erreur Linf Vit.| = fn(t) (convergence)
  275. everr = 'EVOL' 'MANU' 'TEMPS' hltp2 'ERRLINF'
  276. (('LOG' hler2) '/' ('LOG' 10.0D0)) ;
  277. 'DESSIN' everr 'MIMA' 'TITR' 'Log Erreur Linf Vit.' ;
  278. 'FINPROC' ;
  279.  
  280. ** TRCHAMP
  281. 'DEBPROC' TRCHAMP ;
  282. 'ARGUMENT' trx*'TABLE ' ;
  283. trv = (trx . 'EQEX') ;
  284. 'SI' ('NON' ('MULT' (trv . 'PASDETPS' . 'NUPASDT')
  285. (trv . 'NPTRC'))) ;
  286. 'QUITTER' TRCHAMP ;
  287. 'FINSI' ;
  288. tinco = (trv . 'INCO') ;
  289. $dom = (trx . 'DOMZ') ;
  290. dom = ($dom . 'MAILLAGE') ;
  291. cdom = 'CONTOUR' dom ;
  292. ** On calcule tout d'abord la fonction de courant
  293. vitesse = (tinco . 'UN') ;
  294. sw = (-1.) * ('KOPS' vitesse 'ROT' $dom) ;
  295. rk = 'EQEX' $dom 'OPTI' 'EF' 'IMPL'
  296. 'ZONE' $dom 'OPER' 'LAPN' 1.0D0 'INCO' 'PSI'
  297. 'ZONE' $dom 'OPER' 'FIMP' sw 'INCO' 'PSI'
  298. 'CLIM' 'PSI' 'TIMP' cdom 0.0D0 ;
  299. rk . 'INCO' = 'TABLE' 'INCO' ;
  300. rk . 'INCO' . 'PSI' = 'KCHT' $dom 'SCAL' 'SOMMET' 0.0D0 ;
  301. *-*-*-*- A remplacer par EXEC qd ca marchera...
  302. EXIC rk ;
  303. *-*-*-*-
  304. psi = (rk . 'INCO' . 'PSI') ;
  305. ** On effectue une affinité si la cavité est trop allongée
  306. augmen = ((tinco . 'ASPECT') '>' 6.0D0) ;
  307. 'SI' augmen ;
  308. dy = 'NOMC' 'UY' (('COORDONNEE' 2 ($dom . 'SOMMET'))
  309. 'ET' ('COORDONNEE' 2 ($dom . 'CENTRE')))
  310. 'NATU' 'DISCRET' ;
  311. dx = 'NOMC' 'UX' (dy '*' 0.0D0) 'NATU' 'DISCRET' ;
  312. dv = 2.0D0 '*' (dx 'ET' dy) ;
  313. svform = 'FORME' ;
  314. 'FORME' dv ;
  315. 'FINSI' ;
  316. ** Puis on trace les différents champs :
  317. tphy = (trv . 'PASDETPS' . 'TPS') ;
  318. * Vitesses
  319. vref = (tinco . 'UREF') ;
  320. titv = ('CHAINE' 'Champ de vitesse (VECT SOMMET) à t = '
  321. tphy) ;
  322. vv = 'VECTEUR' vitesse (2.0D0 '/' vref) 'UX' 'UY' 'JAUNE' ;
  323. 'TRACER' vv cdom 'TITR' titv ;
  324. * Fonction de courant
  325. 'OPTION' 'ISOV' 'SULI' ;
  326. titpsi = ('CHAINE' 'Fonction de courant (SCAL SOMMET) à t = '
  327. tphy) ;
  328. 'TRACER' psi dom cdom 14 'TITR' titpsi ;
  329. * température
  330. tn = (tinco . 'TN') ;
  331. tittn = ('CHAINE' 'Champ de température (SCAL SOMMET) à t = '
  332. tphy) ;
  333. 'TRACER' tn dom cdom 14 'TITR' tittn ;
  334. * Pression
  335. * NB : On n'utilise pas 'ELNO' : on dessine une pression
  336. * constante sur un élément.
  337. pression = (trv . 'PRESSION' . 'PRESSION') ;
  338. modelp = 'MODELISER' dom 'THERMIQUE' ;
  339. elemp = 'KCHA' $dom 'CHAM' pression ;
  340. titpres = ('CHAINE' 'Pression (SCAL CENTRE) à t = '
  341. tphy) ;
  342. 'TRACER' elemp modelp cdom 'TITR' titpres ;
  343. ** On rétablit le rapport d'aspect initial
  344. 'SI' augmen ; 'FORME' svform ; 'FINSI' ;
  345. 'FINPROC' ;
  346. *
  347. ** Fin de la définition des procédures
  348. **************************************************************
  349.  
  350. **************************************************************
  351. ** Début du 'main'
  352. *
  353. **** Définition des constantes
  354. * l, h, nl, nh : données pour la boite.
  355. * egeom : erreur géométrique pour 'ELIMINATION'
  356. * aspect : rapport d'aspect de la boite
  357. * nu, g : paramètres du fluide
  358. * tg, td : températures sur les murs gauches et droits
  359. h = 1.0D0 ;
  360. 'SI' court ;
  361. nh = 3 ; aspect = 4.0D0 ;
  362. 'SINON' ;
  363. nh = 15 ; aspect = 12.0D0 ;
  364. 'FINSI' ;
  365. nl = 'ENTIER' (nh '*' 4) ;
  366. l = h '*' aspect ;
  367. egeom = 1.0D-4 ;
  368. g = (0.D0 -1.D0) ;
  369. tg = 0.D0 ; td = aspect ;
  370.  
  371. **** Création des maillages
  372. lp = l '/' 2.0D0 ; hp = h '/' 2.0D0 ;
  373. mlp = (-1.0D0) '*' lp ; mhp = (-1.0D0 '*' hp) ;
  374. pmil = lp hp ;
  375. pA = 0.0D0 0.0D0 ; pB = l 0.0D0 ;
  376. pC = l h ; pD = 0.0D0 h ;
  377. pAB = lp 0.0D0 ; pCD = lp h ;
  378. bas = pA 'DROIT' nl pB ; haut = pD 'DROIT' nl pC ;
  379. rwall = pB 'DROIT' nh pC ; lwall = pA 'DROIT' nh pD ;
  380. mediany = pAB 'DROIT' nh pCD ;
  381. mt = 'DALLER' haut rwall bas lwall ;
  382. $mt = 'DOMA' mt 'MACRO' 'IMPR' ;
  383. mt = ($mt . 'MAILLAGE') ; 'TASSER' mt ;
  384. haut = 'CHANGER' haut SEG2 ; bas = 'CHANGER' bas SEG2 ;
  385. lwall = 'CHANGER' lwall SEG2 ; rwall = 'CHANGER' rwall SEG2 ;
  386. mediany = 'CHANGER' mediany SEG2 ;
  387. enleelem = haut 'ELEM' 1 ;
  388. haut2 = 'DIFF' haut enleelem ;
  389. cmt = (lwall 'ET' bas 'ET' rwall 'ET' haut) ;
  390. $haut = 'DOMA' haut 'INCL' $mt egeom ;
  391. $haut2 = 'DOMA' haut2 'INCL' $mt egeom ;
  392. $bas = 'DOMA' bas 'INCL' $mt egeom ;
  393. 'ELIMINATION' (mt 'ET' bas 'ET' rwall 'ET' haut
  394. 'ET' haut2 'ET' lwall 'ET' cmt 'ET' mediany) egeom ;
  395.  
  396. **** Initialisation de la table des historiques
  397. *
  398. thist = 'TABLE' 'HIST' ;
  399. TABTPS = TEMP 'NOEC';
  400. thist . 'tpscpu' = TABTPS.'TEMPS_CPU'.'INITIAL';
  401. * Listes
  402. thist . 'ltps' = 'PROG' ;
  403. thist . 'lerr' = 'PROG' ;
  404. thist . 'lpdt' = 'TABLE' 'LPDT' ; thl = (thist . 'lpdt') ;
  405. thl . 'conv' = 'PROG' ; thl . 'diff' = 'PROG' ;
  406. thl . 'min' = 'PROG' ;
  407.  
  408. **** Initialisation des tables rv et rvp
  409. rv = 'EQEX' $mt
  410. 'OPTI' discconv
  411. 'ZONE' $mt 'OPER' 'NS' 'CLPU' 'CBOU' 'TN' 'TREF'
  412. 'INCO' 'UN'
  413. 'OPTI' discconv
  414. 'ZONE' $mt 'OPER' 'TSCA' 'CLPT' 'UN' 0.0D0
  415. 'INCO' 'TN' ;
  416. rv = 'EQEX' rv
  417. 'ZONE' $haut 'OPER' 'CALCMAR'
  418. 'ZONE' $haut 'OPER' 'TOIMP' 'TSURF'
  419. 'INCO' 'UN' ;
  420. tbas = 'KCHT' $bas 'SCAL' 'SOMMET'
  421. (tg '+' ((('COORDONNEE' 1 bas)
  422. '/' l)
  423. '*' (td '-' tg)) ) ;
  424. * Conditions aux limites :
  425. * traitement special pour pA et pB communs respectivement
  426. * à lwall-bas et a bas-rwall sinon les temperatures valent
  427. * deux fois la valeur qu'on veut leur donner
  428. * On utilise 'MANUEL' car 'CLIM' ne prend que les objets de
  429. * type 'MAILLAGE'
  430. *
  431. rv = 'EQEX' rv
  432. 'CLIM' 'UN' 'UIMP' (lwall 'ET' rwall 'ET' bas) 0.0D0
  433. 'CLIM' 'UN' 'VIMP' (lwall 'ET' rwall 'ET' haut 'ET' bas) 0.0D0
  434. 'CLIM' 'TN' 'TIMP' lwall tg
  435. 'CLIM' 'TN' 'TIMP' ('MANUEL' 'POI1' pA) (-1.0D0 '*' tg)
  436. 'CLIM' 'TN' 'TIMP' rwall td
  437. 'CLIM' 'TN' 'TIMP' ('MANUEL' 'POI1' pB) (-1.0D0 '*' td)
  438. 'CLIM' 'TN' 'TIMP' bas tbas ;
  439. 'SI' graph ;
  440. rv = 'EQEX' rv
  441. 'ZONE' $mt 'OPER' 'MAJHIST'
  442. 'ZONE' $mt 'OPER' 'TRCHAMP' ;
  443. 'FINSI' ;
  444. rvp = 'EQPR' $mt
  445. 'ZONE' $mt 'OPER' 'PRESSION' 0.D0
  446. 'PIMP' 0.D0 ;
  447. rv . 'INCO' = 'TABLE' 'INCO' ;
  448. rv . 'PRESSION' = rvp ;
  449. rv . 'HISTOIRE' = thist ;
  450.  
  451. **** Initialisation des inconnues
  452. * V -> 0 ; T gradient constant avec tg sur lwall, td sur rwall
  453. rv . 'INCO' . 'UN-1' = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0 0.D0) ;
  454. rv . 'INCO' . 'UN' = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0 0.D0) ;
  455. rv . 'INCO' . 'TN' = 'KCHT' $mt 'SCAL' 'SOMMET'
  456. (tg '+' ((('COORDONNEE' 1 mt)
  457. '/' l)
  458. '*' (td '-' tg)) ) ;
  459.  
  460. **** Valeurs des paramètres :
  461. *
  462. ** Physiques
  463. *
  464. * Pr : Nombre de Prandtl
  465. * Ra : Nombre de Rayleign (Gr : Grashof)
  466. * Ma : Nombre de Marangoni (Ma2 : autre définition)
  467. * tref : température de référence (approx. de Boussinesq)
  468. * uref : échelle de vitesse (cf. adimensionnement)
  469. * ampvit : facteur d'amplification pour le tracé des vitesses
  470. 'SI' court ;
  471. Pr = 4.0D0 ; Ra = 0.0D0 ; Ma = 1000.0D0 ;
  472. 'SINON' ;
  473. Pr = 4.24D0 ; Ra = 4300.0D0 ; Ma = 7600.0D0 ;
  474. 'FINSI' ;
  475. Gr = Ra '/' (Pr '*' aspect) ; Ma2 = Ma '/' (Pr '*' aspect) ;
  476. sMa2 = (Ma2 '**' 0.5D0) ;
  477. clpu = 1.0D0 '/' sMa2 ;
  478. cbou = Gr '/' Ma2 ;
  479. clpt = 1.0D0 '/' (sMa2 '*' Pr) ;
  480. clev = 1.0D0 ;
  481. uref = sMa2 ;
  482. tref = (tg + td) / 2.0D0 ;
  483. ampvit = 2.0D0 '/' uref ;
  484. *
  485. ** Numériques
  486. *
  487. * nbiter : nombre d'itérations
  488. * nptrc : fréquence pour le tracé des champs (.ps ou à l'écran)
  489. * fidt : fréquence d'affichage des messages à l'écran
  490. * alpha : multiplicateur du pas de temps
  491. * beta : paramètre de stabilisation de la pression
  492. * pour les éléments MACRO
  493. 'SI' court ;
  494. nbiter = 250 ; nptrc = 125 ; fidt = 10 ;
  495. alpha = 0.9D0 ; beta = 1.0D0 ;
  496. 'SINON' ;
  497. *nbiter = 10 ; nptrc = 500 ; fidt = 100 ;
  498. nbiter = 10000 ; nptrc = 2000 ; fidt = 100 ;
  499. alpha = 0.5D0 ; beta = 1.0D0 ;
  500. 'FINSI' ;
  501. *
  502. ** Remplissage des tables
  503. *
  504. rv . 'INCO' . 'ASPECT' = aspect ;
  505. rv . 'INCO' . 'CLPU' = clpu ;
  506. rv . 'INCO' . 'CBOU' = g '*' cbou ;
  507. rv . 'INCO' . 'TREF' = tref ;
  508. rv . 'INCO' . 'UREF' = uref ;
  509. rv . 'INCO' . 'CLPT' = clpt ;
  510. rv . 'INCO' . 'CLEV' = clev ;
  511. rv . 'ALFA' = alpha ;
  512. rv . 'PRESSION' . 'BETA' = beta ;
  513. rv . 'FIDT' = fidt ;
  514. rv . 'NPTRC' = nptrc ;
  515. rv . 'ITMA' = nbiter ;
  516. *
  517. **
  518. ***
  519. **** Exécution
  520. SUMMARY ;
  521. 'SI' inta ; 'OPTION' 'DONN' 5 ; 'FINSI' ;
  522. EXEC rv ;
  523. TABTPS = TEMP 'NOEC';
  524. thist . 'tpscpu' = (thist . 'tpscpu') '+' TABTPS.'TEMPS_CPU'.'INITIAL';
  525. 'SI' graph ; TRHIST thist ; 'FINSI' ;
  526. ****
  527. ***
  528. **
  529. *
  530. vit = (rv . 'INCO' . 'UN') ;
  531. evvm = 'EVOL' 'BLEU' 'CHPO' vit 'UX' mediany ;
  532. ** Le cas-test donne-t-il un résultat correct ??
  533. 'OPTION' 'ECHO' 0 ;
  534. 'SI' court ;
  535. 'OPTION' 'ECHO' 1 ;
  536. *
  537. ** Vérification :
  538. ** La composante horizontale de la vitesse le long de la
  539. ** médiane verticale doit etre prop. à z(z-2/3)
  540. *
  541. cfpoly evpoly = @POMI evvm 2 ;
  542. cfprop = (cfpoly . 2) ;
  543. cflin = (cfpoly . 1) '/' cfprop ;
  544. cfcte = (cfpoly . 0) '/' cfprop ;
  545. cflana = -2.0D0 '/' 3.0D0 ;
  546. * Résultat obtenu le 14/10/97 :
  547. cfpref = -.3128910770680292D+01 ;
  548. critere1 = 'ABS' ((cfprop '-' cfpref) '/' cfpref) ;
  549. critere2 = 'ABS' ((cflin '-' cflana) '/' cflana) ;
  550. critere3 = 'ABS' cfcte ;
  551. * Résultats obtenus le 14/10/97 :
  552. * cflref = -.6784204323293079D+00 ;
  553. * cfcref = .5876882831320513D-02 ;
  554. * correspondant à :
  555. * critere1 = 2.83861D-16
  556. * critere2 = 1.76306D-02
  557. * critere3 = 5.87688D-03
  558. limite1 = 1.D-07 ;
  559. limite2 = 2.D-02 ;
  560. limite3 = 1.D-02 ;
  561. 'SI' graph ;
  562. tabg = 'TABLE' ; tabg . 'TITRE' = 'TABLE' ;
  563. tabg . 1 = 'NOLI MARQ CROI' ;
  564. tabg . 'TITRE' . 1 = 'CHAINE' 'Résultat Castem2000' ;
  565. tabg . 'TITRE' . 2 = 'CHAINE' 'Parabole théorique' ;
  566. evpoly = ('COULEUR' evpoly 'JAUN') ;
  567. titg = 'CHAINE' 'Ux(z) sur la médiane verticale' ;
  568. 'DESSIN' (evvm 'ET' evpoly) 'MIMA' 'TITR' titg
  569. 'LEGE' tabg ;
  570. 'FINSI' ;
  571.  
  572. * Tps CPU (le 10/10/97) ; 1797 centièmes de secondes
  573. * (court = VRAI ; graph = FAUX ; inta = FAUX ;)
  574. TABTPS = TEMP 'NOEC';
  575. thist .'tpscpu' = (thist . 'tpscpu') '+' TABTPS.'TEMPS_CPU'.'INITIAL';
  576. *
  577. ** Affichage des messages et des erreurs éventuelles
  578. *
  579. 'OPTION' 'ECHO' 0 ;
  580. 'SAUTER' 2 'LIGNE' ;
  581. 'MESSAGE' 'Debriefing :' ;
  582. 'MESSAGE' '------------' ;
  583. 'MESSAGE' 'Comparaison de ux(z) sur la médiane verticale à : ' ;
  584. 'MESSAGE' 'alpha (z**2 - 2/3 z + 0)' ;
  585. 'SAUTER' 'LIGNE' ;
  586. 'MESSAGE' ('CHAINE'
  587. 'Erreur/Limite sur le coefficient alpha :') ;
  588. 'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere1 ' / ' limite1) ;
  589. 'MESSAGE' ('CHAINE'
  590. 'Erreur/Limite sur le coefficient de z (2/3)') ;
  591. 'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere2 ' / ' limite2) ;
  592. 'MESSAGE' ('CHAINE'
  593. 'Erreur/Limite sur le coefficient constant (0)') ;
  594. 'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere3 ' / ' limite3) ;
  595. 'MESSAGE' ('CHAINE' 'Temps CPU / Temps CPU (10/10/97) : ') ;
  596. 'MESSAGE' ('CHAINE' ('ENTIER' (thist . 'tpscpu')) ' / ~1800'
  597. ' centième(s) de secondes') ;
  598. 'SI' ('OU' ('OU' (critere1 '>' limite1)
  599. (critere2 '>' limite2))
  600. (critere3 '>' limite3)) ;
  601. 'ERREUR' 5 ;
  602. 'FINSI' ;
  603. 'MESSAGE' 'Tout s_est bien passé' ;
  604. 'SINON' ;
  605. 'OPTION' 'ECHO' 1 ;
  606. **
  607. ** Vérification pour le cas long :
  608. ** on compare ux(z) sur la médiane avec des résultats
  609. ** experimentaux (attention, l'expérience n'est pas 2D =>
  610. ** vitesse plus faible a cause des parois latérales)
  611. ** et une autre simulation numérique (effectuée avec les
  612. ** memes parametres que la notre)..
  613. ** Voir l'article de Villers et Platten pour plus de précisions
  614. **
  615. zexp = 'PROG' 0.0D0
  616. 0.1D0 0.2D0 0.35D0 0.57D0 0.75D0 0.97D0
  617. 1.2D0 1.5D0 1.8D0 2.0D0 2.2D0 2.5D0 ;
  618. uxexp = 'PROG' 0.0D0
  619. 0.6D0 0.75D0 0.95D0 1.05D0 0.95D0 0.8D0
  620. 0.4D0 -0.4D0 -1.0D0 -2.0D0 -2.5D0 -3.6D0 ;
  621. zpas = (2.5D0 '/' 32) ;
  622. zsvp = 'PROG' 0.0D0 'PAS' zpas 'NPAS' 32 ;
  623. uxsvp = 'PROG' 0.0D0 0.22D0 0.4D0 0.6D0 0.75D0 0.86D0
  624. 0.98D0 1.05D0 1.13D0 1.15D0 1.16D0 1.15D0
  625. 1.13D0 1.08D0 1.02D0 0.92D0 0.81D0 0.7D0
  626. 0.55D0 0.4D0 0.2D0 0.01D0 -0.2D0 -0.43D0
  627. -0.7D0 -0.95D0 -1.2D0 -1.5D0 -1.8D0 -2.15D0
  628. -2.5D0 -2.85D0 -3.25D0 ;
  629. *
  630. * On remet l'évolution obtenu avec C2000 en dimensionné
  631. *
  632. * mu à 20°C densité
  633. nuacet = 0.3265D-3 '/' 870.0D0 ;
  634. lref = 2.5D0 ;
  635. uref = ((nuacet '*' 1.0D6) '/' lref) '*' sMa2 ;
  636. zdim = ('EXTRAIRE' evvm 'ABSC') '*' lref ;
  637. uxdim = ('EXTRAIRE' evvm 'ORDO') '*' uref ;
  638. * Comparaison des profils :
  639. * on les remet sur zdim par interpolation linéaire
  640. uxexpi = 'IPOL' zdim zexp uxexp ;
  641. uxvpi = 'IPOL' zdim zsvp uxsvp ;
  642. lmux = 'PROG' ('MAXIMUM' uxdim 'ABS')
  643. ('MAXIMUM' uxexpi 'ABS')
  644. ('MAXIMUM' uxvpi 'ABS') ;
  645. echux = 'MINIMUM' lmux ;
  646. critere1 = ('MAXIMUM' (uxdim '-' uxvpi) 'ABS') '/' echux ;
  647. critere2 = ('MAXIMUM' (uxdim '-' uxexpi) 'ABS') '/' echux ;
  648. * Résultats obtenus le 15/10/97 :
  649. * critere1 = .7530548488495176D-01
  650. * critere2 = .3228687516012702D+00
  651. limite1 = 0.08D0 ;
  652. limite2 = 0.35D0 ;
  653. * trace eventuel
  654. 'SI' graph ;
  655. evexp = 'EVOL' 'ROSE' 'MANU' 'Z' zexp 'UX' uxexp ;
  656. evsvp = 'EVOL' 'JAUN' 'MANU' 'Z' zsvp 'UX' uxsvp ;
  657. evexpi = 'EVOL' 'ROUG' 'MANU' 'Z' zdim 'UX' uxexpi ;
  658. evvpi = 'EVOL' 'VERT' 'MANU' 'Z' zdim 'UX' uxvpi ;
  659. evdim = 'EVOL' 'TURQ' 'MANU' 'Z' zdim 'UX' uxdim ;
  660. titd = 'CHAINE' 'Profil Ux(z) sur la médiane verticale' ;
  661. tabd = 'TABLE' ; tabd . 'TITRE' = 'TABLE' ;
  662. tabd . 1 = 'MARQ PLUS TIRR TITR Experience' ;
  663. tabd . 2 = 'MARQ CROI TIRC TITR Vp_num' ;
  664. tabd . 3 = 'MARQ PLUS NOLI TITR Exp_ipol' ;
  665. tabd . 4 = 'MARQ CROI NOLI TITR Vp_num_ipol' ;
  666. tabd . 5 = 'MARQ ETOI NOLI TITR K2000' ;
  667. tabd . 'TITRE' . 1 = 'CHAINE' 'Expérience' ;
  668. tabd . 'TITRE' . 2 = 'CHAINE' 'Numérique (V&P)' ;
  669. tabd . 'TITRE' . 3 = 'CHAINE' 'Exp. (interpolé)' ;
  670. tabd . 'TITRE' . 4 = 'CHAINE' 'V&P (interpolé)' ;
  671. tabd . 'TITRE' . 5 = 'CHAINE' 'Castem 2000' ;
  672. 'DESSIN' (evexp 'ET' evsvp 'ET' evexpi 'ET' evvpi 'ET' evdim)
  673. 'LEGE' tabd 'TITR' titd ;
  674. 'FINSI' ;
  675. * Tps CPU (le 10/10/97) ; 7043 secondes
  676. * (court = FAUX ; graph = VRAI ; inta = FAUX ;)
  677. TABTPS = TEMP 'NOEC';
  678. thist .'tpscpu' = (thist . 'tpscpu') '+' TABTPS.'TEMPS_CPU'.'INITIAL';
  679. *
  680. ** Affichage des messages et des erreurs éventuelles
  681. *
  682. 'OPTION' 'ECHO' 0 ;
  683. 'SAUTER' 2 'LIGNE' ;
  684. 'MESSAGE' 'Debriefing :' ;
  685. 'MESSAGE' '------------' ;
  686. 'MESSAGE' 'Comparaison de ux(z) sur la médiane verticale' ;
  687. 'SAUTER' 'LIGNE' ;
  688. 'MESSAGE' 'par rapport à la simulation numérique de V&P' ;
  689. 'MESSAGE' 'Erreur / Limite :' ;
  690. 'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere1 ' / ' limite1) ;
  691. 'SAUTER' 'LIGNE' ;
  692. 'MESSAGE' 'par rapport à l_expérience de V&P' ;
  693. 'MESSAGE' 'Erreur / Limite :' ;
  694. 'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere2 ' / ' limite2) ;
  695. 'SAUTER' 'LIGNE' ;
  696. 'MESSAGE'
  697. ('CHAINE' 'Temps CPU / Temps CPU (15/10/97) (moy2 queue) : ') ;
  698. 'MESSAGE' ('CHAINE' ('ENTIER' ((thist . 'tpscpu')
  699. '/' 100.0D0)) ' / ~7040'
  700. ' secondes') ;
  701. 'SI' ('OU' (critere1 '>' limite1)
  702. (critere2 '>' limite2)) ;
  703. 'ERREUR' 5 ;
  704. 'FINSI' ;
  705. 'MESSAGE' 'Tout s_est bien passé' ;
  706. 'FINSI' ;
  707. SUMMARY ;
  708. 'OPTION' 'ECHO' 1 ;
  709. 'SI' inta ; 'OPTION' 'DONN' 5 ; 'FINSI' ;
  710. 'FIN' ;
  711. *
  712. ** Fin du 'main'
  713. **************************************************************
  714.  
  715.  
  716.  
  717.  
  718.  
  719.  
  720.  
  721.  
  722.  
  723.  
  724.  
  725.  
  726.  

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