Télécharger Marangoni2.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : Marangoni2.dgibi
  2. * _____________________
  3. *
  4. * CAS TEST EFFET MARANGONI
  5. * UTILISATION DE LA LOI DE CHAN, CHEN, MAZUMDER
  6. * Journal of Heat Transfert, vol 110, 1988
  7. * _____________________
  8. * OPTIONS TESTEES PLAN QUAF
  9.  
  10. ***************************************************************
  11. ** DESCRIPTION DU TEST :
  12. **
  13. ** Ce test a pour but d'étudier l'opérateur toimp au travers de l'effet
  14. ** Marangoni.
  15. **
  16. ** Une poche de metal liquide est soumise au bombardement d'un flux de
  17. ** chaleur surfacique ayant une loi parabolique. La loi de CHAN décrit
  18. ** ce type d'écoulement. Cette loi donne le profil de température du
  19. ** liquide ainsi que le profil des vitesses a la surface de la poche
  20. ** liquide.
  21. **
  22. ** Lors de ce test, nous imposons directement sur la surface du bain de
  23. ** liquide le profil de température obtenue grâce a la loi de CHAN. Ce
  24. ** profil de température permet de calculer les forces de Marangoni, seul
  25. ** terme moteur de l'écoulement. Le but est de retomber sur le profil
  26. ** théorique des vitesses donné par la loi de CHAN.
  27. **
  28. ** Ref : Journal of Heat Transfert, vol 110, 1988
  29. ** "Asymptotic solution for thermocapillary flow at high and low
  30. ** Prandtl numbers due to concentrated surface heating."
  31. ** CHAN, CHEN, MAZUMDER.
  32. **
  33. ***************************************************************
  34.  
  35. ***********************
  36. ** OPTIONS A TESTER :
  37. ***********************
  38. option DIME 2 ELEM qua8 MODE plan ; MOAXI = FAUX ;
  39. *option DIME 2 ELEM qua8 MODE axis; MOAXI = VRAI ;
  40. MACRO = QUAF ;
  41.  
  42. ************************************************************
  43. ************************************************************
  44.  
  45. *DIMPR = VRAI ;
  46. DIMPR = FAUX ;
  47. GRAPH = FAUX ;
  48. COMPLET=FAUX ;
  49. option trac psc ;
  50.  
  51. * _________________________________
  52. *
  53. * différents parametres du calcul
  54. * _________________________________
  55.  
  56. * Nombre de mailles
  57. NMAIL = 10 ;
  58.  
  59. * nombre maximal d'itérations
  60. NITER = 8 ; TERR=0.35;
  61. 'SI' COMPLET ; NITER = 50 ; TERR=0.05 ; 'FINSI' ;
  62.  
  63. * coefficient de relaxation pour le champ de vitesse
  64. ALFAU = 0.5 ;
  65.  
  66. * précision requise
  67. PREC = 1.e-3 ;
  68.  
  69. * _______________________________________________________
  70. *
  71. * géométrie du creuset : tronc de cone (dimension en m)
  72. * _______________________________________________________
  73. *
  74.  
  75. * rayon du bas du domaine
  76. rbas = 0.0199 ;
  77.  
  78. * rayon du haut du domaine
  79. rhaut = 0.020 ;
  80.  
  81. * hauteur de liquide
  82. hliq = 0.005 ;
  83.  
  84. * épaisseur de la couche limite
  85. eclm = 1.e-4 ;
  86.  
  87. * rayon moyen du faisceau électonique
  88. rbe = 0.010 ;
  89.  
  90. * __________________________________________________________________
  91. *
  92. * Caractéristiques du faisceau d'électrons à répartition d'énergie
  93. * parabolique : puissance PCANON et "rayon" RFOCAL
  94. * __________________________________________________________________
  95.  
  96. PCANON = 1.e3 ;
  97. RFOCAL = rbe ;
  98.  
  99.  
  100. * ______________________
  101. *
  102. * constantes physiques
  103. * ______________________
  104.  
  105. * constante des gaz parfaits
  106. RG = 8.32 ;
  107.  
  108. * constante de STEPHAN-BOLTZMAN
  109. SIGMA = 5.67E-8 ;
  110.  
  111.  
  112.  
  113. *
  114. * __________________________________________________
  115. *
  116. * propriétés physiques du fluide incompressible
  117. * __________________________________________________
  118.  
  119. * viscosité dynamique (kg/m/s)
  120. mu = 1.3e-3 ;
  121.  
  122. * masse volumique (kg/m3)
  123. rho0 = 6600. ;
  124.  
  125. * conductivité thermique (W/m/K)
  126. lambda = 5. ;
  127.  
  128. * chaleur spécifique (J/kg/K)
  129. Cp = 270 ;
  130.  
  131. * tension de surface (N/m)
  132. sig0 = 0.34 ;
  133.  
  134. * 1/sig0 dsig0/dT (1/K)
  135. gam0 = 5.88e-4 ;
  136.  
  137. * viscosité cinématique (m2/s)
  138. nu = mu / rho0 ;
  139.  
  140. * effet MARANGONI
  141. lev = 1. * gam0 * sig0 / rho0 ;
  142.  
  143. * ______________________
  144. *
  145. * CREATION DU MAILLAGE
  146. * ______________________
  147.  
  148. p0 = 0. 0. ;
  149. p1 = rhaut 0. ;
  150. p2 = rbas (-1. * hliq) ;
  151. p3 = 0. (-1. * hliq) ;
  152.  
  153. p01 = rbe 0. ;
  154. p32 = rbe (-1. * hliq) ;
  155.  
  156. p03 = 0. (-10. * eclm) ;
  157. p12 = rhaut (-10. * eclm) ;
  158.  
  159. d0 = rhaut / NMAIL ;
  160. d1 = rhaut / NMAIL ;
  161. d2 = rhaut / NMAIL ;
  162. surli0 = (p01 droi dini d1 dfin d1 p0 ) coul roug ;
  163. surli1 = (p1 droi dini d2 dfin d1 p01) coul roug ;
  164. surli = surli1 et surli0 ;
  165.  
  166. d0 = hliq / 9 / NMAIL ;
  167. d1 = hliq / 2 / NMAIL ;
  168. d2 = hliq * 3 / NMAIL ;
  169. smm0 = (p0 droi dini d0 dfin d1 p03) coul jaun ;
  170. smm1 = (p03 droi dini d1 dfin d2 p3 ) coul jaun ;
  171. smm = smm1 et smm0 ;
  172.  
  173. scon10 = (inve (smm0 proj 'CYLI' (1. 0.) 'DROI' p1 p2) ) coul vert ;
  174. scon11 = (inve (smm1 proj 'CYLI' (1. 0.) 'DROI' p1 p2) ) coul vert ;
  175. scon1 = scon11 et scon10 ;
  176.  
  177. dz0 = (smm0 proj 'CYLI' (1. 0.) 'DROI' p01 p32) coul bleu ;
  178. dz1 = (smm1 proj 'CYLI' (1. 0.) 'DROI' p01 p32) coul bleu ;
  179. dz = dz0 et dz1 ;
  180.  
  181. zc = (-1. * (rhaut - rbe ) * hliq ) / ( rhaut - rbas) ;
  182.  
  183. scon21= (inve (surli1 proj 'CONI' (rbe zc) 'DROI' p2 p3) ) coul vert ;
  184. scon20= (inve (surli0 proj 'CYLI' (0. -1.) 'DROI' p2 p3) ) coul vert ;
  185. scon2 = scon20 et scon21 ;
  186.  
  187. dr1 = (surli1 proj 'CONI' (rbe zc) 'DROI' p12 p03) coul rose ;
  188. dr0 = (surli0 proj 'CYLI' (0. -1.) 'DROI' p12 p03) coul rose ;
  189. dr = dr1 et dr0 ;
  190.  
  191. scon = scon1 et scon2 ;
  192.  
  193. elim 1.e-6 ( surli1 et surli0 et scon10 et scon11 et scon20
  194. et scon21 et dr0 et dr1 et dz0 et dz1 et smm0 et smm1 ) ;
  195.  
  196. meta1 = daller surli1 dz0 dr1 scon10 ;
  197. meta2 = daller (inve dr1) dz1 scon21 scon11 ;
  198. meta3 = daller (inve dz0) surli0 smm0 dr0 ;
  199. meta4 = daller (inve dz1) (inve dr0) smm1 scon20 ;
  200.  
  201. meta = meta1 et meta2 et meta3 et meta4 ;
  202.  
  203.  
  204. oublier meta1 ;
  205. oublier meta2 ;
  206. oublier meta3 ;
  207. oublier meta4 ;
  208. *trac meta ;
  209. *
  210. ** On réoriente les éléments.
  211. *
  212.  
  213. meta = 'ORIENTER' meta ;
  214. tass meta ;
  215. nbel1 = nbel meta ; list nbel1 ;
  216.  
  217. * transformation des éléments en qua9
  218.  
  219. £META = CHANGER META MACRO ;
  220. £SURLI = CHANGER SURLI MACRO ;
  221. £SCON = CHANGER SCON MACRO ;
  222. £SMM = CHANGER SMM MACRO ;
  223.  
  224. ELIM 1.e-6 (£META ET £SURLI ET £SCON ET £SMM) ;
  225.  
  226. * formulation du modèle NAVIER_STOKES
  227.  
  228. $META = MODE £META 'NAVIER_STOKES' MACRO ;
  229. $SURLI = MODE £SURLI 'NAVIER_STOKES' MACRO ;
  230. $SCON = MODE £SCON 'NAVIER_STOKES' MACRO ;
  231. $SMM = MODE £SMM 'NAVIER_STOKES' MACRO ;
  232.  
  233. mmeta = doma $META maillage ;
  234. ssurli = doma $SURLI maillage ;
  235.  
  236.  
  237.  
  238. * _______________________________________________
  239. *
  240. * SYSTEME D'EQUATIONS ET CONDITIONS AUX LIMITES
  241. * _______________________________________________
  242.  
  243.  
  244.  
  245. * equations de qdm sur ur et uz + div V = 0
  246.  
  247. RV = EQEX $META
  248. OPTI 'EF' 'IMPL' CENTREP1 'SUPG'
  249. ZONE $META OPER KBBT 1 (-1.) INCO 'UN' 'PN'
  250. ZONE $META OPER NS 1. 'UN' nu INCO 'UN'
  251. OPTI 'CENTREE'
  252. ZONE $SURLI OPER TOIMP 'tau' INCO 'UN'
  253. ;
  254.  
  255.  
  256. * conditions aux limites
  257. RV = EQEX RV
  258. CLIM UN UIMP £SMM 0.
  259. UN VIMP £SURLI 0.
  260. ;
  261.  
  262. * ________________
  263. *
  264. * INITIALISATION
  265. * ________________
  266.  
  267. RV.INCO = TABLE INCO ;
  268.  
  269. RV.INCO.'UN' = KCHT $META VECT SOMMET (1.e-8 1.e-8) ;
  270. RV.INCO.'PN' = KCHT $META SCAL CENTREP1 0. ;
  271. RV.INCO.'tau' = KCHT $SURLI VECT CENTRE (0. 0. ) ;
  272.  
  273.  
  274. ********************
  275. * terme tau
  276. * grad T < 0 ,lev > 0 => Marangoni < 0 donc opposé à la qdm
  277. ***********************
  278.  
  279. * Courbe de temperature selon CHAN, CHEN, MAZUMDER
  280. 'SI' ( MOAXI ) ;
  281. bb = 0.70185 ;
  282. cc = 0.35093 ;
  283. f1 = 0.89872 ;
  284. SINON ;
  285. bb = 1. ;
  286. cc = 0.33333 ;
  287. f1 = 1. ;
  288. 'FINSI' ;
  289.  
  290. I0 = PCANON /(PI * (RFOCAL**2.)) ;
  291. q0 = I0 ;
  292. q1 = q0 '/' rbe '/' rbe ;
  293. DTmax = bb * (cc ** -.25) * q0 * (q1 ** -.25)
  294. * ((2 * (rho0 ** 2) * (Cp ** 3) * lev * mu) ** -.25) ;
  295. 'LISTE' DTmax ;
  296. XK = ssurli coor 1 ;
  297. Ttheo = 'KOPS' DTmax * ('KOPS' 1 '-' ('KOPS' (cc / bb / rbe / rbe)
  298. * ('KOPS' XK ** 2 ) ) ) ;
  299. Ttheo = 'KOPS' Ttheo '+' 1200. ;
  300. TK = KCHT $SURLI SCAL SOMMET Ttheo ;
  301.  
  302.  
  303. **************procedure gradient******
  304. DEBPROC GRADS FK*'CHPOINT ' LIGNK*'MAILLAGE' ;
  305. * procédure de calcul de df/ds
  306. * gradient d'une fonction FK suivant l'abscisse curviligne de
  307. * la ligne LIGNK
  308. NELIGNK = NBEL LIGNK ;
  309. REPETER BLIGNK NELIGNK ;
  310. ELEMI = ELEM LSURLI &BLIGNK ;
  311. PO1I = POIN ELEMI INITIAL ;
  312. PO2I = POIN ELEMI FINAL ;
  313. XK1I = COOR 1 PO1I ;
  314. XK2I = COOR 1 PO2I ;
  315. ZK1I = COOR 2 PO1I ;
  316. ZK2I = COOR 2 PO2I ;
  317. FK1I = EXTR FK 'SCAL' PO1I ;
  318. FK2I = EXTR FK 'SCAL' PO2I ;
  319. DXKI = XK2I - XK1I ;
  320. DZKI = ZK2I - ZK1I ;
  321. DSKI = ( ( DXKI * DXKI ) + ( DZKI * DZKI ) ) ** 0.5 ;
  322. DFKI = FK2I - FK1I ;
  323. SI ( DXKI &lt;EG 0. ) ;
  324. SIGNI = -1. ;
  325. SINON ;
  326. SIGNI = 1. ;
  327. FINSI ;
  328. DFDI = 0.5 * SIGNI * DFKI / DSKI ;
  329. CDFDI = 'MANU' 'CHPO' (PO1I ET PO2I ) 1 SCAL ( PROG 2 * DFDI )
  330. 'NATURE' DISC ;
  331. SI ( &BLIGNK EGA 1 ) ;
  332. DFDS = CDFDI ;
  333. SINON ;
  334. DFDS = DFDS ET CDFDI ;
  335. FINSI ;
  336. FIN BLIGNK ;
  337. FINPROC DFDS ;
  338. *************** fin de la procedure *******
  339.  
  340.  
  341. ** calcul de l'effet Marangoni
  342. LSURLI = CHANGER SURLI SEG2 ;
  343. GTK = GRADS TK LSURLI ;
  344. MARA = KOPS lev '*' GTK ;
  345. MARA = 'NOMC' UX MARA ;
  346. RV.INCO.'tau' = KCHT $SURLI VECT CENTRE MARA ;
  347. MARA = KCHT $SURLI VECT CENTRE MARA ;
  348. MARAN = ELNO $SURLI MARA ;
  349.  
  350. * __________________
  351. *
  352. * BOUCLE ITERATIVE
  353. * __________________
  354.  
  355.  
  356. REPETER BLOCKI NITER ;
  357.  
  358. MESS 'ITERATION INTERNE N°' &BLOCKI ;
  359.  
  360. UN0 = COPI RV.INCO.'UN' ;
  361.  
  362.  
  363.  
  364. *
  365. * CALCUL D'UNE ITERATION
  366. *
  367. RV.METHINV.FCPRECT = 1 ;
  368. RV.METHINV.FCPRECI = 1 ;
  369. EXEC RV ;
  370.  
  371. *
  372. * TESTS DE CONVERGENCE
  373. *
  374. UNR = EXCO UX RV.INCO.'UN' ;
  375. UNZ = EXCO UY RV.INCO.'UN' ;
  376. URMAX = MAXI UNR ;
  377. UZMAX = MAXI UNZ ;
  378.  
  379. MESS 'URMAX = ' URMAX ;
  380. MESS 'UZMAX = ' UZMAX ;
  381.  
  382.  
  383. UR0 = UN0 EXCO 'UX' scal ;
  384. UR0 = KCHT $META SCAL SOMMET UR0 ;
  385. UNR = KCHT $META SCAL SOMMET UNR ;
  386. ERUNR = KOPS ( KOPS UNR '-' UR0 ) '/' ( (MAXI (ABS UNR) ) + 1.e-13) ;
  387. ERRUNR = MAXI ( ABS ERUNR ) ;
  388. MESS 'ERRUNR = ' ERRUNR ;
  389.  
  390. UZ0 = EXCO UY UN0 ;
  391. UZ0 = KCHT $META SCAL SOMMET UZ0 ;
  392. UNZ = KCHT $META SCAL SOMMET UNZ ;
  393. ERUNZ = KOPS ( KOPS UNZ '-' UZ0 ) '/' ( (MAXI (ABS UNZ) ) + 1.e-13) ;
  394. ERRUNZ = MAXI ( ABS ERUNZ ) ;
  395. MESS 'ERRUNZ = ' ERRUNZ ;
  396.  
  397.  
  398. ERR = ERRUNR ; VAR = MOTS 'UNR' ;
  399.  
  400. SI ( ERRUNZ > ERR ) ;
  401. ERR = ERRUNZ ; VAR = MOTS 'UNZ' ;
  402. FINSI ;
  403.  
  404.  
  405. VAR1 = EXTR 1 VAR ;
  406. MESS 'erreur relative maximale sur ' VAR1 ERR ;
  407.  
  408. * RELAXATION
  409.  
  410. RV.INCO.'UN' = KOPS ( KOPS ALFAU '*' RV.INCO.'UN')
  411. '+' ( KOPS (1. - ALFAU) '*' UN0 ) ;
  412.  
  413.  
  414. SI (ERR &lt;EG PREC) ;
  415. * convergence des itérations internes
  416. MESS 'CONVERGENCE EN ' &BLOCKI 'ITERATIONS' ;
  417. QUITTER BLOCKI ;
  418. FINSI ;
  419.  
  420. SI (&BLOCKI EGA NITER ) ;
  421. MESS 'PAS DE CONVERGENCE EN ' NITER 'ITERATIONS' ;
  422. FINSI ;
  423.  
  424. ******************************
  425. ** FIN DE LA BOUCLE ITERATIVE
  426. *******************************
  427. FIN BLOCKI ;
  428.  
  429. *************************************
  430. ** COURBE THEORIQUE DE Ur SELON CHAN
  431. *************************************
  432. XK = ssurli coor 1 ;
  433. URtheo = ((cc * q1 * 2 * lev / mu / Cp) ** .5) * f1 ;
  434. URtheo = 'KOPS' URtheo * XK ;
  435. URtheo = KCHT $SURLI SCAL SOMMET URtheo ;
  436. evURtheo = evol roug chpo URtheo (inve ssurli ) ;
  437.  
  438. ******************
  439. ** IMPRESSIONS
  440. ******************
  441. 'SI' (DIMPR) ;
  442.  
  443. titre 'temperature imposee le long de surli.' ;
  444. evtk = evol roug chpo TK (inve ssurli ) ;
  445. 'SI' GRAPH ; dess evtk ; 'FINSI' ;
  446.  
  447. titre 'vitesses calculees (ampl= 0.001) (m/s)' ;
  448. UN = RV.INCO.'UN' ;
  449. VIT1 = VECT UN 0.001 UX UY ROUG ;
  450. 'SI' GRAPH ; trac VIT1 mmeta ; 'FINSI' ;
  451.  
  452. titre 'vitesse radiale le long de surli (m/s) (en rouge, la theorie)' ;
  453. evux = evol chpo UN UX (inve ssurli ) ;
  454. 'SI' GRAPH ; dess (evux 'ET' evURtheo) 'LEGE' ; 'FINSI' ;
  455.  
  456. titre 'vitesse radiale theorique le long de surli (m/s)' ;
  457. 'SI' GRAPH ; dess evURtheo ; 'FINSI' ;
  458.  
  459. titre 'contrainte de Marangoni le long de surli (Pa)' ;
  460. evmaran = evol chpo MARAN UX (inve ssurli ) ;
  461. evmaran = -1. * evmaran ;
  462. evmaran = rho0 * evmaran ;
  463. 'SI' GRAPH ; dess evmaran ; 'FINSI' ;
  464.  
  465. * affichage de la composante imposee tau
  466. nomper = 'EXTRAIRE' 3 (rv . 'LISTOPER') ;
  467. notable= 'MOT' ('TEXTE' ('CHAINE' 3 nomper)) ;
  468. mess 'Operateur ' nomper ;
  469. msi mai=('TEXTE' nomper) (rv . notable) ;
  470. msir = redu msi surli ;
  471. ev1 = evol vert chpo msir 1UN (INVE surli) ;
  472. 'SI' ( MOAXI ) ;
  473. XK = ssurli coor 1 ;
  474. XK = KCHT $SURLI SCAL SOMMET XK ;
  475. MARAN = 'EXCO' MARAN UX ;
  476. MARAN = KCHT $SURLI SCAL SOMMET MARAN ;
  477. marx = KOPS (KOPS MARAN * (-1. * rhaut / 2 / NMAIL )) * (2 * pi * XK) ;
  478. SINON ;
  479. MARAN = 'EXCO' MARAN UX ;
  480. marx = KOPS MARAN * (-1. * rhaut / 2 / NMAIL ) ;
  481. 'FINSI' ;
  482. titre 'effet marangoni par maille (la theorie en rouge)' ;
  483. evmarx = evol roug chpo marx (inve ssurli ) ;
  484. 'SI' GRAPH ; dess (evmarx et ev1) 'LEGE' ; 'FINSI' ;
  485.  
  486. 'FINSI' ;
  487.  
  488. ************
  489. titre 'vitesse radiale le long de surli (m/s) (en rouge, la theorie)' ;
  490. UN = RV.INCO.'UN' ;
  491. evux = evol chpo UN UX (inve ssurli ) ;
  492. 'SI' GRAPH ; dess (evux 'ET' evURtheo) 'LEGE' ; 'FINSI' ;
  493.  
  494. **** TEST DE FONCTIONNEMENT ***********
  495. UNR = EXCO UX RV.INCO.'UN' ;
  496. UNR = 'REDU' UNR SURLI ;
  497. UNR = KCHT $SURLI SCAL SOMMET UNR ;
  498. ERtheo = KOPS ( KOPS UNR - URtheo ) '/' ( (MAXI (ABS UNR) ) + 1.e-13) ;
  499. ERRtheo = MAXI ( ABS ERtheo ) ;
  500. MESS 'erreur de la vitesse radiale par rapport à la théorie =' ERRtheo ;
  501.  
  502. 'SI' (ERRtheo < TERR) ;
  503. 'ERREUR' 0 ;
  504. 'SINON' ;
  505. 'ERREUR' 5 ;
  506. 'FINSI' ;
  507.  
  508. ***************************************
  509.  
  510. 'FIN' ;
  511.  
  512.  
  513.  
  514.  
  515.  
  516.  
  517.  
  518.  
  519.  
  520.  
  521.  
  522.  
  523.  

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