Télécharger linekmanimp.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : linekmanimp.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ****************************************************************************
  5. *
  6. *
  7. * Date : 06/01/99
  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. * La solution est permanente.
  41. *
  42. * Cet écoulement admet une solution analytique car
  43. * les termes de transport sont négligés.
  44. *
  45. * Type de calcul : PERMANENT, PAS DE TERME DE TRANSPORT
  46. *
  47. * Type de méthode : IMPLICITE
  48. *
  49. * Type d'éléments : ELEMENTS QUA9 CENTREP1
  50. *
  51. * Opérateurs utilisés : KBBT LAPN MDIA
  52. *
  53. *
  54. *
  55. * Auteurs : Gilles Bernard-Michel
  56. * : Daniel Guilbaud pour la présente version implicite
  57. *
  58. *****************************************************************************
  59. *
  60.  
  61. *
  62. ********************* Les options utilisateur ******************************
  63. *
  64.  
  65. *
  66. ** Déformation du maillage.
  67. *
  68.  
  69. DEFORM = FAUX ;
  70.  
  71. *
  72. ** Sortie graphique.
  73. *
  74.  
  75. GRAPH = FAUX ;
  76.  
  77. *
  78. ** Test long ou court.
  79. *
  80.  
  81. COMPLET = FAUX ;
  82.  
  83. **************************************************************************
  84.  
  85. option dime 2 elem qua8 mode axis;
  86.  
  87. *
  88. ************************ On définit les propriétés physiques ***************
  89. *
  90. ** viscosité dynamique mu/rho.
  91. *
  92.  
  93. nu =1.E-1 ;
  94.  
  95. *
  96. ** La vitesse de rotation initiale.
  97. *
  98.  
  99. omeg = 1 ;
  100.  
  101. *
  102. ** Le nombre de Rossby.
  103. *
  104.  
  105. Rosb = 0.0 ;
  106.  
  107. *
  108. ** L'impulsion.
  109. *
  110.  
  111. ddom = 1 ;
  112.  
  113. *
  114. ********************* On traite les aspects géométriques ********************
  115. *
  116.  
  117. *
  118. ** Les paramètres de la géométrie.
  119. *
  120.  
  121. rint = 0.5 ;
  122. rext = 1.;
  123. hh = 1. ;
  124.  
  125. *
  126. ** Densité du maillage.
  127. *
  128.  
  129. SI (NON COMPLET);
  130. nptx = 10;
  131. nptz = 20;
  132. SINON;
  133. nptx = 30;
  134. nptz = 60;
  135. FINSI;
  136.  
  137. *
  138. ** Constantes caractérisant les différentes couches limites.
  139. *
  140.  
  141. ekma = nu / (abs omeg) ;
  142. ddek = ekma**0.5;
  143.  
  144. *
  145. ** Les sommets du bol
  146. *
  147.  
  148. p1 = rint 0. ;
  149. p2 = rext 0. ;
  150. p3 = rext hh ;
  151. p4 = rint hh ;
  152.  
  153. *
  154. ** Assemblage des bords du domaine.
  155. *
  156.  
  157. bas = d nptx p1 p2;
  158. cdro = d nptz p2 p3;
  159. haut = d nptx p3 p4;
  160. cgau = d nptz p4 p1;
  161.  
  162. *
  163. ** On maille.
  164. *
  165.  
  166. cnt = bas et cdro et haut et cgau;
  167. mt= bas cdro haut cgau daller ;
  168. tass mt ;
  169.  
  170. *
  171. ** On réoriente les éléments.
  172. *
  173.  
  174. mt = orienter mt ;
  175.  
  176. *
  177. ** On déforme le maillage avec un bruit gaussien.
  178. ** Sous réserve que l'option soit selectionnée.
  179. *
  180.  
  181. SI (DEFORM);
  182.  
  183. ** Amplitude des déformation (nombre arbitraire, 1 donne une déformation
  184. ** acceptable pour le maillage).
  185.  
  186. ampdef = 1. ;
  187.  
  188. mmt = CHANGE mt POI1;
  189. ccnt= CHANGE cnt POI1;
  190. inmt= diff mmt ccnt;
  191. denn = (0.03 * ampdef) / nptx;
  192. brbl1 = BRUIT BLAN GAUS 0 denn inmt;
  193. brbl2 = BRUIT BLAN GAUS 0 denn inmt;
  194. brbl = (NOMC 'UR' brbl1) ET (NOMC 'UZ' brbl2);
  195. DEPLAC mt PLUS brbl;
  196. FINSI;
  197.  
  198. *
  199. ** On trace le maillage obtenu.
  200. *
  201.  
  202. SI GRAPH ;
  203. titre 'tracé du maillage' ;
  204. trace mt ;
  205. FINSI ;
  206.  
  207. *
  208. ** On crée le domaine associé au maillage.
  209. *
  210.  
  211. * transformation des éléments en qua9
  212.  
  213. £mt = CHANGER mt QUAF ;
  214. £haut = CHANGER haut QUAF ;
  215. £cnt = CHANGER cnt QUAF ;
  216. £cgau = CHANGER cgau QUAF ;
  217. £bas = CHANGER bas QUAF ;
  218. £cdro = CHANGER cdro QUAF ;
  219.  
  220.  
  221.  
  222.  
  223. ELIM 1.e-6 (£mt et £cnt et £haut et £cgau et £bas et £cdro) ;
  224.  
  225. * formulation du modèle NAVIER_STOKES
  226.  
  227. $mt = MODE £mt 'NAVIER_STOKES' QUAF ;
  228.  
  229. mmt = doma $mt maillage ;
  230.  
  231. *
  232. ************************ On défini les conditions sur les vitesses **************
  233. *
  234.  
  235. *
  236. ** On se donne l'échelle de grandeur pour le tracage.
  237. *
  238.  
  239. uref= 0.2 ;
  240.  
  241. *
  242. ** Calcul de la solution exacte.
  243. *
  244. corx = coor 1 mmt ;
  245. cory = coor 2 mmt ;
  246. zba = cory / ddek;
  247. dzba = zba * 180. / PI ;
  248. mzba = -1.*zba;
  249. urdo = ddom * corx * (EXP mzba) * (SIN dzba) ;
  250. udom = NOMC 'UX' (kcht $mt SCAL SOMMET urdo) 'NATU' 'DISCRET';
  251. wzdo = -1. * ddom * ddek * (1. -
  252. ((EXP mzba) * ((SIN dzba) + (COS dzba))));
  253. wdom = NOMC 'UY' (kcht $mt SCAL SOMMET wzdo) 'NATU' 'DISCRET';
  254. upla = kcht $mt VECT SOMMET (udom ET wdom);
  255. vdom = ddom * corx * (EXP mzba) * (COS dzba) ;
  256. pdom = -2. * ddom * ekma * (EXP mzba) * (SIN dzba) ;
  257. zb = 1. / ddek;
  258. dzb = zb * 180. / PI ;
  259. mzb = -1 * zb;
  260. pdelt = 2. * ddom * ekma * (EXP mzb) * (SIN dzb) ;
  261. pdom = pdom + pdelt;
  262.  
  263. *
  264. ****************** Mise en place du calcul numérique ********************
  265. *
  266.  
  267.  
  268. * equations de qdm sur ur et uz + div V = 0
  269.  
  270. RV = EQEX $mt
  271. OPTI 'EF' 'IMPL' 'CENTREP1' 'SUPG'
  272. * il faut opti CENTREP1 pour KBBT
  273. ZONE $mt OPER KBBT (+1.) INCO 'UN' 'PN'
  274. ZONE $mt OPER LAPN nu INCO 'UN'
  275. ZONE $mt OPER MDIA 'cour' INCO 'VT' 'UN'
  276. ;
  277.  
  278. * equation scalaire sur ut
  279.  
  280. RV = EQEX RV
  281. OPTI 'EF' 'IMPL' 'SUPG'
  282. ZONE $mt OPER LAPN nu INCO 'VT'
  283. ZONE $mt OPER MDIA 'nusr' INCO 'VT'
  284. ZONE $mt OPER MDIA 'cout' INCO 'UN' 'VT'
  285. ;
  286.  
  287. * conditions aux limites
  288.  
  289. RV = EQEX RV
  290. 'CLIM' 'UN' 'UIMP' £cnt urdo
  291. 'UN' 'VIMP' (£cgau et £bas et £cdro) wzdo
  292. 'VT' 'TIMP' £cnt vdom
  293. ;
  294.  
  295.  
  296. * ________________
  297. *
  298. * INITIALISATION
  299. * ________________
  300.  
  301. RV.INCO = TABLE INCO ;
  302.  
  303. RV.INCO.'UN' = KCHT $mt VECT SOMMET (1.e-8 1.e-8) ;
  304. RV.INCO.'VT' = KCHT $mt SCAL SOMMET 1.e-8 ;
  305. RV.INCO.'PN' = KCHT $mt SCAL CENTREP1 0. ;
  306. RV.INCO.'cour' = KCHT $mt VECT SOMMET ( (-2.*omeg) 0. ) ;
  307. RV.INCO.'cout' = KCHT $mt VECT SOMMET ( (+2.*omeg) 0. ) ;
  308.  
  309. * terme -nu/r2
  310.  
  311. coorx = 'KCHT' $mt scal SOMMET corx ;
  312. nusr2 = KOPS ( +1. * nu ) '/' ( KOPS coorx '*' coorx ) ;
  313.  
  314. RV.INCO.'nusr' = KCHT $mt SCAL SOMMET nusr2 ;
  315.  
  316. * __________
  317. *
  318. * CALCUL
  319. * __________
  320.  
  321. EXEC RV ;
  322.  
  323.  
  324. *
  325. ** Calcul de l'erreur par rapport à la solution analytique.
  326. *
  327.  
  328. urdo = KCHT $mt SCAL SOMMET urdo ;
  329. wzdo = KCHT $mt SCAL SOMMET wzdo ;
  330. vdom = KCHT $mt SCAL SOMMET vdom ;
  331. pdom = KCHT $mt SCAL SOMMET pdom ;
  332.  
  333.  
  334. UNR = EXCO UX RV.INCO.'UN' ;
  335. UNR = KCHT $mt SCAL SOMMET UNR ;
  336. ERUNR = KOPS ( KOPS UNR '-' urdo) '/' ( (MAXI (ABS UNR) ) + 1.e-13) ;
  337. ERRUNR = MAXI ( ABS ERUNR ) ;
  338. MESS 'ERRUNR = ' ERRUNR ;
  339.  
  340. UNZ = EXCO UY RV.INCO.'UN' ;
  341. UNZ = KCHT $mt SCAL SOMMET UNZ ;
  342. ERUNZ = KOPS ( KOPS UNZ '-' wzdo) '/' ( (MAXI (ABS UNZ) ) + 1.e-13) ;
  343. ERRUNZ = MAXI ( ABS ERUNZ ) ;
  344. MESS 'ERRUNZ = ' ERRUNZ ;
  345.  
  346. VT = RV.INCO.'VT' ;
  347. ERVT = KOPS ( KOPS VT '-' vdom) '/' ( (MAXI (ABS VT) ) + 1.e-13) ;
  348. ERRVT = MAXI ( ABS ERVT ) ;
  349. MESS 'ERRVT = ' ERRVT ;
  350.  
  351. PNN = ELNO $mt (RV.'INCO'.'PN') CENTREP1 ;
  352. pnn = kcht $mt SCAL SOMMET PNN;
  353. ERPN = KOPS ( KOPS pnn '-' pdom) '/' ( (MAXI (ABS pnn) ) + 1.e-13) ;
  354. ERRPN = MAXI ( ABS ERPN ) ;
  355. MESS 'ERRPN = ' ERRPN ;
  356.  
  357. SI ( DEFORM ) ;
  358. erv = 2.e-3 ;
  359. erp = 1.e-1 ;
  360. SINON ;
  361. erv = 1.e-5 ;
  362. erp = 5.e-2 ;
  363. FINSI ;
  364.  
  365. BOOL1 = ERRUNR <EG erv ;
  366. BOOL2 = ERRUNZ <EG erv ;
  367. BOOL3 = ERRVT <EG erv ;
  368. BOOL4 = ERRPN <EG erp ;
  369.  
  370. BOOLT = BOOL1 et BOOL2 et BOOL3 et BOOL4 ;
  371.  
  372. SI BOOLT ;
  373. MESS 'TEST CORRECT' ;
  374. SINON ;
  375. MESS 'TEST FAUX' ;
  376. ERREUR 5 ;
  377. FINSI ;
  378.  
  379. *
  380. **************** Traitement graphique ***************************************
  381. *
  382.  
  383. SI GRAPH ;
  384.  
  385. *
  386. ** Tracé pour les vitesses dans le plan (rz).
  387. *
  388. UN = RV.INCO.'UN' ;
  389. ung1=vect un uref ux uy jaune ;
  390. titre 'vitesse VR VZ du calcul numérique' ;
  391. trac ung1 £mt titre 'vitesse VR VZ du calcul numérique' ;
  392. titre 'vitesse VT du calcul numérique' ;
  393. trac vt £mt ;
  394. titre 'pression numérique' ;
  395. trac pnn £mt ;
  396.  
  397. upla1=vect upla uref ux uy roug ;
  398. titre 'vitesse VR VZ du calcul analytique' ;
  399. trac upla1 £mt titre 'vitesse VR VZ du calcul analytique' ;
  400. titre 'vitesse VT du calcul analytique' ;
  401. trac vdom £mt ;
  402. titre 'pression analytique' ;
  403. trac pdom £mt ;
  404.  
  405. titre 'difference relative de vitesse UR' ;
  406. trac ERUNR £mt ;
  407.  
  408. titre 'difference relative de vitesse UZ' ;
  409. trac ERUNZ £mt ;
  410.  
  411. titre 'difference relative de vitesse VT' ;
  412. trac ERVT £mt ;
  413.  
  414. titre 'difference relative de pression' ;
  415. trac ERPN £mt ;
  416.  
  417. *
  418. ** Tracé pour la vitesse ortho-radiale en 3D.
  419. *
  420. VT = RV.INCO.'VT' ;
  421. vtx = kcht $mt SCAL SOMMET 0;
  422. vtx = NOMC 'UX' vtx 'NATU' DISCRET;
  423. vty = kcht $mt SCAL SOMMET 0.;
  424. vty = NOMC 'UY' vty 'NATU' DISCRET;
  425. * La vitesse est inversée car le maillage
  426. * est en réalité dans le plan rz et non xy.
  427. * Pour recreer la véritable orientation, on change
  428. * le signe de la vitesse ortho-radiale.
  429. vtz = NOMC 'UZ' (-1.*vt) 'NATU' DISCRET;
  430. uuu = kcht $mt VECT SOMMET (vtx ET vty ET vtz);
  431. ung2 = vect uuu uref ux uy uz jaune;
  432. titre 'vitesse azimutale 3D' ;
  433. trac ung2 mt titre 'vitesse azimutale 3D' ;
  434.  
  435. FINSI ;
  436.  
  437. fin ;
  438.  
  439.  
  440.  
  441.  
  442.  
  443.  
  444.  

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