Télécharger linekman.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : linekman.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ****************************************************************************
  5. *
  6. *
  7. * Date : 19/3/97
  8. *
  9. * Description: Cas test simulant l'écoulement dans une région
  10. * limitée par une plaque horizontale infinie en
  11. * rotation autour d'un axe perpendiculaire.
  12. *
  13. * l'ensemble du fluide tourne à une vitesse Omeg
  14. * tandis que la plaque tourne à la vitesse
  15. * (Omeg + ddom).
  16. *
  17. * Les calculs étant effectués en axi-symétrique,
  18. * seule la moitié de la région est représentée.
  19. * L'axe de rotation constitue l'axe des z de notre
  20. * étude. La plaque inférieure a pour équation z = 0.
  21. *
  22. * La région étant par ailleurs infinie, on est bien
  23. * entendu obligé de la limiter dans notre étude.
  24. * On met les conditions adéquates au bord.
  25. *
  26. *
  27. * Objet : Le but de ce test est de valider la mise
  28. * en place sous CASTEM 2000 de calculs en
  29. * "vrai axi-symétrique", c'est à dire pour
  30. * lesquels la composante ortho-radiale de
  31. * la vitesse n'est pas nulle mais par contre
  32. * toutes les dérivées par rapport à la composante
  33. * ortho-radiale du repère cylindrique sont nulles.
  34. *
  35. * Les résultats issus du calcul numérique seront
  36. * comparés aux résultats prédits par la théorie
  37. * (voir M. Ungarish, Hydrodynalics of Suspensions,
  38. * p. 47, Springer-Verlag, 1993): linear ekman layer.
  39. *
  40. * Cet écoulement admet une solution analytique.
  41. *
  42. *
  43. * Auteur : Gilles Bernard-Michel
  44. *
  45. * Révision : Ce fichier diffère du précédent car il n'y a
  46. * pas d'implicitation.
  47. *
  48. *****************************************************************************
  49. *
  50.  
  51. *
  52. ** Choix des éléments
  53. *
  54.  
  55. option dime 2 elem qua8 mode axis;
  56.  
  57. *
  58. ************************ Quelques procédures ********************************
  59. *
  60.  
  61. *
  62. ** Erreur Linfini entre deux Champoints.
  63. *
  64.  
  65. DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ;
  66. vmax = MAXI vit 'ABS' + 1.e-13;
  67. dimax = MAXI (vitp1 - vit) 'ABS' ;
  68. err = (dimax /vmax) ;
  69. FINPROC err ;
  70.  
  71. *
  72. ** Calcul de la norme L2 d'un champ.
  73. *
  74.  
  75. DEBPROC CALCNOL2 vit*'CHPOINT' mod*'MMODEL'
  76. mt*'MAILLAGE';
  77. opti dime 2 mode PLAN;
  78. vgau = CHANGE 'CHAM' vit mt 'RIGIDITE';
  79. usol2 = vgau * vgau;
  80. us2 = INTG mod usol2;
  81. opti dime 2 mode axis;
  82. FINPROC us2 ;
  83.  
  84. *
  85. ** Calcul de l'erreur L2 pour une composante de la vitesse.
  86. *
  87.  
  88. DEBPROC CALCURR vit*'CHPOINT' vitsol*'CHPOINT' mod*'MMODEL'
  89. mt*'MAILLAGE';
  90. opti dime 2 mode PLAN;
  91. vgau = CHANGE 'CHAM' vit mt 'RIGIDITE';
  92. vsolgau = CHANGE 'CHAM' vitsol mt 'RIGIDITE';
  93. udiff = vgau - vsolgau;
  94. udif2 = udiff*udiff;
  95. ud2 = INTG mod udif2;
  96. usol2 = vsolgau * vsolgau;
  97. us2 = INTG mod usol2;
  98. us2 = us2 + 1.e-14;
  99. err = ud2 / us2;
  100. err = err**0.5;
  101. opti dime 2 mode axis;
  102. FINPROC err ;
  103.  
  104. *
  105. ** Erreur L2 pour l'ensemble des composantes de la vitesse.
  106. *
  107.  
  108. DEBPROC CALCERTO vit1*'CHPOINT' vit2*'CHPOINT' vit3*'CHPOINT'
  109. vi1*'CHPOINT' vi2*'CHPOINT' vi3*'CHPOINT' mod*'MMODEL'
  110. mt*'MAILLAGE';
  111. vit1 = vit1 - vi1;
  112. vit2 = vit2 - vi2;
  113. vit3 = vit3 - vi3;
  114. nom1 = CALCNOL2 vit1 mod mt;
  115. nom2 = CALCNOL2 vit2 mod mt;
  116. nom3 = CALCNOL2 vit3 mod mt;
  117. nomi = nom1 + nom2 + nom3;
  118. deno1 = CALCNOL2 vi1 mod mt;
  119. deno2 = CALCNOL2 vi2 mod mt;
  120. deno3 = CALCNOL2 vi3 mod mt;
  121. denom = deno1 + deno2 + deno3 + 1.e-14;
  122. err = nomi/denom;
  123. err = err**0.5;
  124. FINPROC err;
  125.  
  126. *
  127. ************************ On définit les propriétés physiques ***************
  128. *
  129. ** viscosité dynamique mu/rho.
  130. *
  131.  
  132. nu =1.E-1 ;
  133.  
  134. *
  135. ** La vitesse de rotation initiale.
  136. *
  137.  
  138. omeg = 1 ;
  139.  
  140. *
  141. ** Le nombre de Rossby.
  142. *
  143.  
  144. Rosb = 0.0 ;
  145.  
  146. *
  147. ** L'impulsion.
  148. *
  149.  
  150. ddom = 1 ;
  151.  
  152. *
  153. ********************* Les options utilisateur ******************************
  154. *
  155.  
  156. *
  157. ** Critère de convergence implicite.
  158. *
  159.  
  160. crit = 0.07 ;
  161.  
  162. *
  163. ** Nombre max de boucles d'implicitation.
  164. *
  165.  
  166. MAXIMP = 1 ;
  167.  
  168. *
  169. ** Déformation du maillage.
  170. *
  171.  
  172. DEFORM = FAUX;
  173.  
  174. *
  175. ** Amplitude des déformation (nombre arbitraire, 1 donne une déformation
  176. ** acceptable pour le maillage).
  177. *
  178.  
  179. ampdef = 1. ;
  180.  
  181. *
  182. ** Sortie graphique.
  183. *
  184.  
  185. GRAPH = 'N';
  186.  
  187. *
  188. ** Test long ou court.
  189. *
  190.  
  191. COMPLET = FAUX;
  192.  
  193. *
  194. ** On calcule les erreurs L2.
  195. *
  196.  
  197. ERRCALC = VRAI;
  198.  
  199. *
  200. ** Periode de calcul des erreurs L2.
  201. *
  202.  
  203. period = 5;
  204.  
  205. *
  206. ********************* On traite les aspects géométriques ********************
  207. *
  208.  
  209. *
  210. ** Les paramètres de la géométrie.
  211. *
  212.  
  213. rint = 0.5 ;
  214. rext = 1.;
  215. hh = 1. ;
  216.  
  217. *
  218. ** Densité du maillage.
  219. *
  220.  
  221. SI (NON COMPLET);
  222. nptx = 3;
  223. nptz = 4;
  224. NBITER = 60;
  225. SINON;
  226. nptx = 30;
  227. npty = 60;
  228. NBITER = 15000;
  229. FINSI;
  230.  
  231. *
  232. ** Constantes caractérisant les différentes couches limites.
  233. *
  234.  
  235. ekma = nu / omeg;
  236. ddek = ekma**0.5;
  237.  
  238. *
  239. ** Les sommets du bol
  240. *
  241.  
  242. p1 = rint 0. ;
  243. p2 = rext 0. ;
  244. p3 = rext hh ;
  245. p4 = rint hh ;
  246.  
  247. *
  248. ** Assemblage des bords du domaine.
  249. *
  250.  
  251. bas = d nptx p1 p2;
  252. cdro = d nptz p2 p3;
  253. haut = d nptx p3 p4;
  254. cgau = d nptz p4 p1;
  255.  
  256. *
  257. ** On maille.
  258. *
  259.  
  260. cnt = bas et cdro et haut et cgau;
  261. mt= bas cdro haut cgau daller ;
  262. tass mt ;
  263.  
  264. *
  265. ** On trace le maillage obtenu.
  266. *
  267.  
  268. *trace mt ;
  269. *opti donn 5;
  270.  
  271. *
  272. ** On réoriente les éléments.
  273. *
  274.  
  275. mt = orienter mt ;
  276.  
  277. *
  278. ** On déforme le maillage avec un bruit gaussien.
  279. ** Sous réserve que l'option soit selectionnée.
  280. *
  281.  
  282. SI (DEFORM);
  283. mmt = CHANGE mt POI1;
  284. ccnt= CHANGE cnt POI1;
  285. inmt= diff mmt ccnt;
  286. denn = (0.03 * ampdef) / nptx;
  287. brbl1 = BRUIT BLAN GAUS 0 denn inmt;
  288. brbl2 = BRUIT BLAN GAUS 0 denn inmt;
  289. brbl = (NOMC 'UR' brbl1) ET (NOMC 'UZ' brbl2);
  290. DEPLAC mt PLUS brbl;
  291. FINSI;
  292.  
  293. *
  294. ** On crée le domaine associé au maillage.
  295. *
  296.  
  297. $mt = DOMA mt 'MACRO' impr;
  298. mt = $mt.maillage;
  299.  
  300. *
  301. ** On crée un modèle pour pouvoir intégrer des champs.
  302. *
  303. MODL1 = MODE mt 'MECANIQUE';
  304.  
  305. *
  306. ************************ On défini les conditions sur les vitesses **************
  307. *
  308.  
  309. *
  310. ** On se donne l'échelle de grandeur pour le tracage.
  311. *
  312.  
  313. uref= 0.2 ;
  314.  
  315. *
  316. ** le terme gb défini la présence de terme source
  317. ** pour l'équation radiale et son absence dans
  318. ** l'équation axiale.
  319. *
  320.  
  321. gb = (-1.0 0.0) ;
  322.  
  323. *
  324. ** Calcul de la solution exacte.
  325. *
  326.  
  327. corx = coor 1 mt ;
  328. cory = coor 2 mt ;
  329. zba = cory / ddek;
  330. dzba = zba * 180 / 3.141592654;
  331. mzba = -1*zba;
  332. urdo = ddom * corx * (EXP mzba) * (SIN dzba) ;
  333. udom = NOMC 'UX' (kcht $mt SCAL SOMMET urdo) 'NATU' 'DISCRET';
  334. wzdo = -1 * ddom * ddek * (1 -
  335. ((EXP mzba) * ((SIN dzba) + (COS dzba))));
  336. wdom = NOMC 'UY' (kcht $mt SCAL SOMMET wzdo) 'NATU' 'DISCRET';
  337. upla = kcht $mt VECT SOMMET (udom ET wdom);
  338. vdom = ddom * corx * (EXP mzba) * (COS dzba) ;
  339. pdom = -2. * ddom * ekma * (EXP mzba) * (SIN dzba) ;
  340. zb = 1. / ddek;
  341. dzb = zb * 180 / 3.141592654;
  342. mzb = -1 * zb;
  343. pdelt = 2. * ddom * ekma * (EXP mzb) * (SIN dzb) ;
  344. pdom = pdom + pdelt;
  345.  
  346. *
  347. ****************** Mise en place du calcul numérique ********************
  348. *
  349.  
  350. *
  351. ** On crée une table RX pour stocker les vitesses.
  352. *
  353.  
  354. RX = TABLE ;
  355.  
  356. *
  357. ** On initialise les CHPOINTS utilisés à t=0.
  358. *
  359.  
  360. RX.'UN' = kcht $mt VECT SOMMET (0. 0.) ;
  361. RX.'VT' = kcht $mt SCAL SOMMET 0. ;
  362.  
  363. *
  364. ** Création de la table RV associé à la résolution des équations
  365. ** radiales et axiales. On utilise la version Boussinesq pour avoir
  366. ** les termes sources aux sommets et non au centre.
  367. *
  368.  
  369. RV = EQEX $mt OPTI 'ALE' 'ITMA' 1 'ALFA' 0.8
  370. 'ZONE' $mt 'OPER' 'NS' nu gb 'SR' 0. 'WN' 'INCO' 'UN'
  371. 'ZONE' $mt 'OPER' 'TSCA' nu 'WN' 'ST' 'INCO' 'VT'
  372. ;
  373.  
  374.  
  375. RV = EQEX RV
  376. 'CLIM' 'UN' 'UIMP' cnt urdo
  377. 'UN' 'VIMP' (cgau et bas et cdro) wzdo
  378. 'VT' 'TIMP' cnt vdom
  379. ;
  380.  
  381. *
  382. ** On charge la pression (méthode de Choleski).
  383. *
  384.  
  385. RVP = EQPR $mt KTYPI 1
  386. ZONE $mt 'OPER' 'PRESSION' 0.
  387. * PIMP 0.
  388. ;
  389.  
  390. RV.'INCO' = TABLE 'INCO';
  391. RV.PRESSION = RVP;
  392.  
  393. *
  394. ** Les champs transportants sont nuls à tout t>=0.
  395. ** On initialise aussi les champs de vitesse.
  396. *
  397.  
  398. RV.INCO.'WN' = kcht $mt VECT SOMMET (0. 0.) ;
  399. RV.INCO.'UN' = copi RX.'UN';
  400. RV.INCO.'VT' = copi RX.'VT';
  401.  
  402. *
  403. ****************** Boucle itérative sur les pas de temps ***************
  404. *
  405.  
  406. indice = 0 ;
  407. nper = 1;
  408.  
  409. REPETER ITER NBITER;
  410.  
  411. indice = indice + 1 ;
  412.  
  413. *
  414. ** stocke les vitesses dans des variables intermédiaires.
  415. *
  416.  
  417. un = RX.'UN' ;
  418. vt = RX.'VT' ;
  419.  
  420. *
  421. ** Bouclage du calcul implicite des termes sources et du
  422. ** terme transportant.
  423. *
  424.  
  425. nimpl = 0 ;
  426.  
  427. REPETER IMPLIC MAXIMP ;
  428.  
  429. nimpl = nimpl + 1;
  430.  
  431. *
  432. ** On charge les termes sources implicites.
  433. *
  434.  
  435. uimpl = kcht $mt VECT SOMMET RV.INCO.'UN';
  436. vimpl = kcht $mt SCAL SOMMET RV.INCO.'VT';
  437. sr = 2. * omeg * vimpl ;
  438. uscal = exco RV.INCO.'UN' ux scal;
  439. st = -2. * omeg * uscal - (nu*vt/(corx*corx));
  440. stt = kcht $mt SCAL SOMMET st;
  441. st = noel $mt stt;
  442. RV.INCO.'SR' = kcht $mt SCAL SOMMET sr;
  443. RV.INCO.'ST' = kcht $mt SCAL CENTRE st;
  444.  
  445. *
  446. ** On reinitialise les tables de vitesses non implicites
  447. *
  448.  
  449. RV.INCO.'UN' = kcht $mt VECT SOMMET un ;
  450. RV.INCO.'VT' = kcht $mt SCAL SOMMET vt;
  451. RV.PRESSION =RVP ;
  452.  
  453. *
  454. ** On execute le calcul des équations radiales et axiales.
  455. *
  456.  
  457. EXEC RV ;
  458.  
  459. *
  460. ** On test la convergence de notre implicitation.
  461. ** Ce qui a un sens pour MAXIMP >= 2.
  462. *
  463.  
  464. SI (MAXIMP > 1);
  465.  
  466. err1 = CALCERR uimpl RV.INCO.'UN';
  467. err2 = CALCERR vimpl RV.INCO.'VT';
  468.  
  469. *
  470. ** On charge la premiere evolution.
  471. *
  472.  
  473. SI (nimpl EGA 1);
  474. errdebu = err1 ;
  475. errdebv = err2 ;
  476. FINSI;
  477.  
  478. *
  479. ** On effectue le test.
  480. *
  481.  
  482. bool1 = nimpl > 1;
  483. bool2 = err1 <EG (crit*errdebu);
  484. bool3 = err2 <EG (crit*errdebv);
  485. bool = bool1 et bool2 et bool3;
  486. SI (bool);
  487. QUITTER IMPLIC;
  488. FINSI;
  489.  
  490. SI ((nimpl EGA MAXIMP) ET (NON bool) ET (nimpl > 1)) ;
  491. mess 'erreur d' 'implicitation' err1 err2;
  492. FINSI;
  493. FINSI;
  494.  
  495.  
  496. FIN IMPLIC ;
  497.  
  498. *
  499. ** Calcul de l'erreur par rapport à la solution analytique.
  500. *
  501.  
  502. SI (ERRCALC ET (indice EGA (nper*period)));
  503. un1 = exco RV.INCO.'UN' ux;
  504. un2 = RV.INCO.'VT';
  505. un3 = exco RV.INCO.'UN' uy;
  506. pn = RV.PRESSION.PRESSION;
  507. pnn = elno $mt pn;
  508. pnn = kcht $mt SCAL SOMMET pnn;
  509. errl2u = CALCURR un1 urdo MODL1 mt;
  510. errl2v = CALCURR un2 vdom MODL1 mt;
  511. errl2w = CALCURR un3 wzdo MODL1 mt;
  512. errl2p = CALCURR pnn pdom MODL1 mt;
  513. errtot = CALCERTO un1 un2 un3 urdo vdom wzdo MODL1 mt;
  514.  
  515. *mess 'pas de temps' indice ;
  516. *mess 'valeur des derivees' errdebu errdebv;
  517. *mess 'erreur l2' errtot errl2u errl2v errl2w errl2p;
  518.  
  519. *
  520. ** Stockage des erreurs en listes.
  521. *
  522.  
  523. SI (nper EGA 1);
  524. lerl2u = PROG errl2u;
  525. lerl2v = PROG errl2v;
  526. lerl2w = PROG errl2w;
  527. lerl2p = PROG errl2p;
  528. lerl2t = PROG errtot;
  529. liter = PROG indice;
  530. SINON; SI (nper > 1);
  531. lerl2u = lerl2u ET (PROG errl2u);
  532. lerl2v = lerl2v ET (PROG errl2v);
  533. lerl2w = lerl2w ET (PROG errl2w);
  534. lerl2p = lerl2p ET (PROG errl2p);
  535. lerl2t = lerl2t ET (PROG errtot);
  536. liter = liter ET (PROG indice);
  537. FINSI;
  538. FINSI;
  539. nper = nper + 1;
  540. FINSI;
  541.  
  542. *
  543. ** On sauve les expressions obtenues pour les vitesses.
  544. *
  545.  
  546. RX.'UN' = copi RV.INCO.'UN' ;
  547. vt = RV.INCO.'VT';
  548. vtn = kcht $mt SCAL SOMMET vt;
  549. RX.'VT' = copi vtn ;
  550.  
  551. FIN ITER ;
  552.  
  553. *
  554. ******************** Test de l'erreur ************************
  555. *
  556.  
  557. SI (NON COMPLET);
  558. erdif1 = ABS (errtot - 2.51e-2);
  559. mess ' erdif1 = ' erdif1 ' errtot = ' errtot;
  560. bool1 = erdif1 <EG 1.e-4;
  561. erdif2 = ABS (errl2p - .69732);
  562. mess ' erdif2 = ' erdif2 ' errl2p = ' errl2p;
  563. bool2 = erdif2 <EG 1.e-4;
  564. bool3 = bool1 ET bool2;
  565. SI (NON bool3);
  566. mess 'houra';
  567. ERREUR 5;
  568. FINSI;
  569. SINON;
  570. bool1 = errtot < 0.02e-2;
  571. bool2 = errl2p < 0.01;
  572. bool3 = bool1 ET bool2;
  573. SI (NON bool3);
  574. ERREUR 5;
  575. FINSI;
  576. FINSI;
  577.  
  578. *
  579. **************** Traitement graphique ***************************************
  580. *
  581.  
  582. SI ('EGA' GRAPH 'O');
  583.  
  584. *
  585. ** Tracé pour les vitesses dans le plan (rz).
  586. *
  587.  
  588. ung1=vect un uref ux uy jaune ;
  589. trac ung1 mt ;
  590.  
  591. *
  592. ** Tracé pour la vitesse ortho-radiale en 3D.
  593. *
  594.  
  595. vtx = kcht $mt SCAL SOMMET 0;
  596. vtx = NOMC 'UX' vtx 'NATU' DISCRET;
  597. vty = kcht $mt SCAL SOMMET 0.;
  598. vty = NOMC 'UY' vty 'NATU' DISCRET;
  599. * La vitesse est inversée car le maillage
  600. * est en réalité dans le plan rz et non xy.
  601. * Pour recreer la véritable orientation, on change
  602. * le signe de la vitesse ortho-radiale.
  603. vtz = NOMC 'UZ' (-1.*vt) 'NATU' DISCRET;
  604. uuu = kcht $mt VECT SOMMET (vtx ET vty ET vtz);
  605. ung2 = vect uuu uref ux uy uz jaune;
  606. trac ung2 mt;
  607.  
  608. *
  609. ** tracé de l'erreur d'implicitation en fonction du
  610. ** pas de temps.
  611. *
  612.  
  613. UL2 = EVOL 'MANU' 'ITERATION' liter
  614. 'ERREUR L2 EN UR' lerl2u;
  615. VL2 = EVOL 'MANU' 'ITERATION' liter
  616. 'ERREUR L2 EN UTETA' lerl2v;
  617. WL2 = EVOL 'MANU' 'ITERATION' liter
  618. 'ERREUR L2 EN UZ' lerl2w;
  619. PL2 = EVOL 'MANU' 'ITERATION' liter
  620. 'ERREUR L2 EN P' lerl2p;
  621. TL2 = EVOL 'MANU' 'ITERATION' liter
  622. 'ERREUR L2 VITESSE' lerl2t;
  623. dessin ul2 logy;
  624. dessin vl2 logy;
  625. dessin wl2 logy;
  626. dessin tl2 logy;
  627. dessin pl2 logy;
  628. FINSI;
  629.  
  630. *
  631. ** Sauvegarde des données.
  632. *
  633.  
  634. *REPIX RV;
  635. *REPIX RA;
  636. *fsauv = '/test4/?????.sauv' ;
  637. *OPTI SAUV fsauv; SAUV;
  638.  
  639. FIN;
  640.  
  641.  
  642.  
  643.  
  644.  
  645.  
  646.  

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