Télécharger back_proj_2.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : back_proj_2.dgibi
  2. **************** CAS TEST : back_proj_2.dgibi ********************
  3. *
  4. * Ce test permet de vérifier le bon fonctionnement des opérateurs
  5. * utilisés pour la résolution des équations de NAVIER_STOKES en EF
  6. * par un algorithme de projection.
  7. *
  8. * Le cas étudié est celui d'un écoulement laminaire dans un canal
  9. * en présence d'une marche descendante. On teste la position du
  10. * point de réattachement par rapport au pied de la marche.
  11. *
  12. * Référence : K. Morgan, J. Périaux anf F. Thomasset, editors,
  13. * Analysis of laminar flow over a backward facing step, Vol9 of
  14. * Notes on numerical Fluid Mechanics, Vieweg, 1984.
  15. *
  16. *-------------------------------------------------------------------
  17. * Auteur : F.Dabbene (LTMF) 03/03
  18. *-------------------------------------------------------------------
  19. *
  20. 'SAUT' 'PAGE' ;
  21. 'OPTION' 'DIME' 2 'ELEM' 'QUA8' 'ECHO' 0 ;
  22. *
  23. *- Pilotage du calcul
  24. *
  25. * DISCR : Element en vitesse (LINE/MACRO/QUAF)
  26. * KPRES : Element en pression (CENTRE/CENTREP2/MSOMMET)
  27. * KSUPG : Méthode de décentrement (CENTREE/SUPG/SUPGDC)
  28. * KMASS : Matrice masse lumpée ou non (EFM1/EF)
  29. *
  30. * DT : Valeur du pas de temps
  31. * (En cas de problème de convergence, augmenter DT et/ou
  32. * prendre KMASS en EFM1)
  33. * EPS0 : Test d'arret
  34. * RAF : Taille de maille définie par 0.1*raf
  35. *
  36. * FREQ0 : Fréquence d'évaluation du résidu
  37. * GRAPH : Booleen pour les tracés
  38. * COMPLET : Booleen volume de calcul
  39. *
  40. DISCR = 'LINE' ;
  41. KPRES = 'MSOMMET' ;
  42. KSUPG = 'SUPG' ;
  43. KMASS = 'EFM1' ;
  44. GRAPH = FAUX ;
  45. COMPLET = FAUX ;
  46. 'SI' COMPLET ;
  47. DT = 3. ;
  48. FREQ0 = 20 ;
  49. EPS0 = 1.D-5 ;
  50. RAF = 1. ;
  51. 'SINO' ;
  52. DT = 3. ;
  53. FREQ0 = 20 ;
  54. EPS0 = 1.D-2 ;
  55. RAF = 2. ;
  56. 'FINSI' ;
  57. *
  58. *---------------------------------------------------------------------------
  59. * Recherche du point de réattachement (point où dUx/dy=0)
  60. *---------------------------------------------------------------------------
  61. * 1/ Après avoir calculé le gradient de Ux, on ne conserve que les valeurs
  62. * sur la frontière qui nous intéresse sous la forme d'une évolution.
  63. * 2/ On borne l'évolution au voisinage du point de réattachement afin
  64. * d'avoir une variation monotone sur l'intervale de dUx/dy.
  65. * 3/ On recherche par interpolation le zero de la fonction (dUx/dy(s),s)
  66. * La valeur obtenu est l'abscisse curviligne cherchée
  67. *
  68. * Remarques :
  69. * (i) 1/ permettrait de calculer le coeff de frottement à la paroi :
  70. * il suffirait de diviser EV1 par le bon coefficient (Re/2 ici)
  71. * (ii) Il est impératif que dUx/dy soit monotone sur l'intervale 2/ afin
  72. * qu'il y ait unicité du zero (principe des valeurs intermédiaires)
  73. *---------------------------------------------------------------------------
  74. *
  75. 'DEBPROC' attac ;
  76. * 1/
  77. Ux = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' 'COMP' 'UX' 0. RV.'INCO'.'UN' ;
  78. DUxDY = 'EXCO' 'UY' ('KOPS' Ux 'GRADS' $DOMTOT) 'SCAL' ;
  79. EV1 = 'EVOL' 'CHPO' DUxDY BOTTOM ;
  80. * 2/
  81. EV2 = 'EXTR' EV1 'APRE' 5. ;
  82. EV3 = 'EXTR' EV2 'AVAN' 10. ;
  83. * 3/
  84. LX1 = 'EXTR' EV3 'ABSC' ;
  85. LY1 = 'EXTR' EV3 'ORDO' ;
  86. LY1 = 'ORDO' LY1 ;
  87. Ymin = 'MINI' LY1 ;
  88. Ymax = 'MAXI' LY1 ;
  89. Delta = Ymin * Ymax ;
  90. 'SI' (Delta < 0.)
    Y0 = 0. ;
  91. 'SINON' ;
  92. 'MESS' 'Fonction 'NON' monotone' ;
  93. Y0 = Ymin ;
  94. * 'ERRE''Fonction 'NON' monotone' ;
  95. 'FINSI' ;
  96. V1 = 'IPOL' Y0 LY1 LX1 ;
  97. 'FINP' V1 ;
  98. *
  99. *==================================================================
  100. * Calcul du résidu basé sur la composante horizontale de la vitesse
  101. * et arrêt suivant un critère transmis
  102. *==================================================================
  103. * E/ : RVX : TABLE : TABLE des données créées par EQEX
  104. * ARG1 : Fréquence d'impression
  105. * ARG2 : Critère d'arrêt
  106. * ARG3 : Booleen de tracé
  107. * /S : MAT1 : MATRIK : Objet vide
  108. * /S : CHP1 : CHPO : Objet vide
  109. *------------------------------------------------------------------
  110. * MAT1 et CHP1 permettent d'assurer la compatibilité des opérateurs
  111. * de discrétisation avec les procédures personnelles
  112. *------------------------------------------------------------------
  113. 'DEBPROC' residu rvx*table ;
  114. RV = RVX . 'EQEX' ;
  115. FREQ = RVX . 'ARG1' ;
  116. EPS0 = RVX . 'ARG2' ;
  117. GRAPH = RVX . 'ARG3' ;
  118. NITER = RV . 'NITER' ;
  119. DD = RV . 'PASDETPS' . 'NUPASDT' ;
  120. NN = DD '/' FREQ ;
  121. L0 = 'EGA' (DD '-' (FREQ*NN)) 0 ;
  122. 'SI' L0 ;
  123. RANG0 = RV . 'PASDETPS' . 'NUPASDT' ;
  124. TIME0 = RV . 'PASDETPS' . 'TPS' ;
  125. UN0 = 'EXCO' 'UX' RV . 'INCO' . 'UN' 'SCAL' ;
  126. UNM0 = 'EXCO' 'UX' RV . 'INCO' . 'UN2' 'SCAL' ;
  127. ERR0 = ('MAXIMUM' ('ABS' (UN0 '-' UNM0))) '+' 1.D-20 ;
  128. ERR10 = ('LOG' ERR0 ) '/' ('LOG' 10.) ;
  129. 'MESSAGE' 'Résidu en vitesse suivant X au pas'
  130. RANG0 '(t=' TIME0 ') :' ERR0 ':' ERR10 ;
  131. RV . 'INCO' . 'IT' = RV . 'INCO' . 'IT' 'ET' ('PROG' RANG0) ;
  132. RV . 'INCO' . 'TI' = RV . 'INCO' . 'TI' 'ET' ('PROG' TIME0) ;
  133. RV . 'INCO' . 'ER' = RV . 'INCO' . 'ER' 'ET' ('PROG' ERR10) ;
  134. V1 = attac ;
  135. RV . 'INCO' . 'POSI' = RV . 'INCO' . 'POSI' 'ET' ('PROG' V1) ;
  136. Y1 = ('LOG' EPS0) '/' ('LOG' 10) ; Y2 = 0. ;
  137. 'SI' GRAPH ;
  138. EV1 = 'EVOL' 'MANU' (RV . 'INCO' . 'IT')(RV . 'INCO' . 'ER') ;
  139. 'DESSIN' EV1 'YBOR' Y1 Y2 'NCLK' ;
  140. 'FINSI' ;
  141. 'SI' ((ERR10 < Y1) 'ET' (DD > ('MAXI' ('LECT' 10 FREQ)))) ;
  142. RV . 'TFINAL' = RV . 'PASDETPS' . 'TPS' ;
  143. 'FINSI' ;
  144. 'FINSI' ;
  145. RV . 'INCO' . 'UN2' = 'COPIER' RV . 'INCO' . 'UN' ;
  146. mat1 chp1 = 'KOPS' 'MATRIK' ;
  147. 'FINP' mat1 chp1 ;
  148.  
  149. *
  150. *==========
  151. * Maillage
  152. *==========
  153. *
  154. *
  155. *-----------------------------------------------------------------------
  156. *
  157. *
  158. * < L1 X L2 - L1 >
  159. *
  160. * 7___________6__________________________________________________5
  161. * ^ ^
  162. * INLET H1
  163. * v H2 OUTLET
  164. * 1___________2
  165. * | v
  166. * 3__________________________________________________4
  167. * BOTTOM
  168. *
  169. *
  170. *-----------------------------------------------------------------------
  171. *
  172. *
  173. * L1 : Longueur de la section d'entrée (avant la marche)
  174. * L2 : Longueur de la totalité du dispositif
  175. * H1 : Hauteur de la section d'entrée
  176. * H2 : Hauteur de la section de sortie
  177. * d1 : Dimension caractéristique d'une maille
  178. L1 = 3.0 ;
  179. L2 = 22.0 ;
  180. H1 = 1.0 ;
  181. H2 = 1.5 ;
  182. d1 = 0.1 * raf ;
  183. *
  184. * H2-H1 : Hauteur de la marche servant à l'adimensionnalisation
  185. HDIM = H2 - H1 ;
  186. L1 = L1 / HDIM ;
  187. L2 = L2 / HDIM ;
  188. H1 = H1 / HDIM ;
  189. H2 = H2 / HDIM ;
  190. d1 = d1 / HDIM ;
  191. *
  192. * Points du maillage
  193. p1 = 0. (H2-H1) ;
  194. p2 = L1 (H2-H1) ;
  195. p3 = L1 0. ;
  196. p4 = L2 0. ;
  197. p5 = L2 H2 ;
  198. p6 = L1 H2 ;
  199. p7 = 0. H2 ;
  200. p8 = L2 (H2-H1) ;
  201. *
  202. * Section d'éntrée
  203. p1p2 = p1 'DROIT' 'DINI' d1 'DFIN' d1 p2 ;
  204. p2p6 = p2 'DROIT' 'DINI' d1 'DFIN' d1 p6 ;
  205. p6p7 = p6 'DROIT' 'DINI' d1 'DFIN' d1 p7 ;
  206. p7p1 = p7 'DROIT' 'DINI' d1 'DFIN' d1 p1 ;
  207. mesh1 = 'DALL' p1p2 p2p6 p6p7 p7p1 'PLAN' ;
  208. *
  209. * Section de sortie
  210. p6p2 = 'INVE' p2p6 ;
  211. p2p3 = p2 'DROIT' 'DINI' d1 'DFIN' d1 p3 ;
  212. p3p4 = p3 'DROIT' 'DINI' d1 'DFIN' d1 p4 ;
  213. p4p8 = p4 'DROIT' 'DINI' d1 'DFIN' d1 p8 ;
  214. p8p5 = p8 'DROIT' 'DINI' d1 'DFIN' d1 p5 ;
  215. p5p6 = p5 'DROIT' 'DINI' d1 'DFIN' d1 p6 ;
  216. mesh2 = 'DALL' (p6p2 'ET' p2p3) p3p4 (p4p8 'ET' p8p5) p5p6 'PLAN' ;
  217. *
  218. *
  219. *=======================
  220. * Modèles NAVIER_STOKES
  221. *=======================
  222. *
  223. *
  224. * Définition des équations
  225. * $DOMTOT : Modèle volumique défini sur le maillage complet
  226. *
  227. * Conditions aux limites en vitesse :
  228. * $INLET : Modèle surfacique défini à l'entrée fluide (Poiseuille)
  229. * $WALL : Modèle surfacique défini sur les murs (adhérence en paroi)
  230. *
  231. * Conditions aux limites en pression :
  232. * $OUTLET : Modèle surfacique défini à la sortie fluide (sert à imposer
  233. * la pression de sortie en pression continue MSOMMET)
  234. *
  235. * Post-traitement
  236. * $BOTTOM : Modèle surfacique défini sur le plancher après la marche
  237. * (sert à évaluer la position du point de réattachement)
  238. *
  239. DOMTOT = 'CHAN' 'QUAF' (mesh1 ET mesh2) ;
  240. $DOMTOT = 'MODE' DOMTOT 'NAVIER_STOKES' DISCR ;
  241. $INLET = 'MODE' p7p1 'NAVIER_STOKES' DISCR ;
  242. $OUTLET = 'MODE' (p4p8 'ET' p8p5) 'NAVIER_STOKES' DISCR ;
  243. $BOTTOM = 'MODE' p3p4 'NAVIER_STOKES' DISCR ;
  244. $WALL = 'MODE' (p1p2 ET p2p3 ET p3p4 ET p5p6 ET p6p7)
  245. 'NAVIER_STOKES' DISCR ;
  246. *
  247. * Elimination ad hoc
  248. * (En 2D, il faut éliminer les points centres des modèles surfaciques
  249. * avec les points faces des modèles volumiques à cause des MACROs)
  250. FDOMTOT = 'DOMA' $DOMTOT 'FACE' ;
  251. CINLET = 'DOMA' $INLET 'CENTRE' ;
  252. COUTLET = 'DOMA' $OUTLET 'CENTRE' ;
  253. CBOTTOM = 'DOMA' $BOTTOM 'CENTRE' ;
  254. CWALL = 'DOMA' $WALL 'CENTRE' ;
  255. 'ELIM' (FDOMTOT 'ET' CINLET 'ET' COUTLET 'ET' CWALL 'ET' CBOTTOM) EPS0 ;
  256. *
  257. * On écrase les anciens maillages afin d'éviter toute ambiguitée
  258. DOMTOT = 'DOMA' $DOMTOT 'MAILLAGE' ;
  259. INLET = 'DOMA' $INLET 'MAILLAGE' ;
  260. OUTLET = 'DOMA' $OUTLET 'MAILLAGE' ;
  261. BOTTOM = 'DOMA' $BOTTOM 'MAILLAGE' ;
  262. WALL = 'DOMA' $WALL 'MAILLAGE' ;
  263. *
  264. * Maillage pour d'éventuelles conditions aux limites en pression
  265. 'SI' ('EGA' kpres 'MSOMMET') ;
  266. OUTLETP = 'DOMA' $outlet KPRES ;
  267. 'SINON' ;
  268. OUTLETP = 'DOMA' $domtot KPRES ;
  269. 'FINSI' ;
  270. *
  271. *===========================
  272. * Description des équations
  273. *===========================
  274. *
  275. * Grandeurs adimensionnées
  276. Umax = 1.0 ;
  277. Re = 150. ;
  278. *
  279. * Profil de vitesse parabolique à l'entrée
  280. YINLET = 'COOR' 2 INLET ;
  281. YMAX = 'MAXI' YINLET ;
  282. YMIN = 'MINI' YINLET ;
  283. UIN = (YINLET '-' YMAX) '*' (YINLET '-' YMIN) ;
  284. UIN = UIN '*' (-4.0*Umax/((YMAX-YMIN)*(YMAX-YMIN))) ;
  285. UIN = 'NOMC' 'UX' UIN 'NATU' 'DISCRET' ;
  286. VIN = 'KCHT' $INLET 'SCAL' 'SOMMET' 'COMP' 'UY' 0. ;
  287. *
  288. * Description du système en vitesse-pression
  289. RV = 'EQEX' $DOMTOT 'ITMA' 5000 'ALFA' 1. 'FIDT' 20000
  290. 'OPTI' 'EF' 'IMPL' KSUPG KPRES
  291. 'ZONE' $DOMTOT 'OPER' residu FREQ0 EPS0 GRAPH
  292. 'ZONE' $DOMTOT 'OPER' 'NS' 1. 'UN' (1./Re) 'INCO' 'UN'
  293. * Conditions aux limites par défaut
  294. * 'ZONE' $OUTLET 'OPER' 'TOIM' (0. 0.) 'INCO' 'UN'
  295. 'OPTI' KMASS KSUPG
  296. 'ZONE' $DOMTOT 'OPER' 'DFDT' 1. 'UN' DT 'UN' (1./Re) 'INCO' 'UN'
  297. ;
  298. 'SI' ('EGA' kpres 'MSOMMET') ;
  299. RVP = 'EQEX' 'OPTI' 'EF' KPRES 'IMPL'
  300. 'ZONE' $DOMTOT 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES'
  301. 'CLIM' 'PRES' 'TIMP' OUTLETP 0.
  302. ;
  303. 'SINON' ;
  304. RVP = 'EQEX' 'OPTI' 'EF' KPRES 'IMPL'
  305. 'ZONE' $DOMTOT 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES'
  306. * En milieu ouvert, pression discontinue il est inutile
  307. * d'imposer la pression en un point
  308. * 'CLIM' 'PRES' 'TIMP' (OUTLETP 'ELEM' 1) 0.
  309. ;
  310. 'FINSI' ;
  311. *
  312. * Description des conditions aux limites
  313. RV = 'EQEX' RV
  314. 'CLIM' 'UN' 'UIMP' WALL 0. 'UN' 'VIMP' WALL 0.
  315. 'UN' 'UIMP' INLET UIN 'UN' 'VIMP' INLET VIN
  316. ;
  317. *
  318. * Déclaration des inconnues et initialisations (table INCO)
  319. RV . 'INCO' = 'TABLE' 'INCO' ;
  320. RV . 'INCO' . 'UN' = 'KCHT' $DOMTOT 'VECT' 'SOMMET' (0. 0.) ;
  321. RV . 'INCO' . 'PRES' = 'KCHT' $DOMTOT 'SCAL' KPRES 0. ;
  322. *
  323. * Champs supplémentaires pour la procédure residu
  324. RV . 'INCO' . 'IT' = 'PROG' ;
  325. RV . 'INCO' . 'TI' = 'PROG' ;
  326. RV . 'INCO' . 'ER' = 'PROG' ;
  327. RV . 'INCO' . 'POSI' = 'PROG' ;
  328. RV . 'INCO' . 'UN2' = 'KCHT' $DOMTOT 'VECT' 'SOMMET' (1.D-3 1.D-3) ;
  329. *
  330. * Algorithme de résolution : méthode de projection
  331. RV . 'PROJ' = RVP ;
  332. *
  333. EXEC RV ;
  334. *
  335. *=================
  336. * Post-traitement
  337. *=================
  338. *
  339. * Localisation du point de réattachement
  340. * Tracés de la vitesse et de la pression
  341. *
  342. CNT1 = WALL ;
  343. NLI0 = 14 ;
  344. *
  345. * Point de réattachement (point où dUx/dy=0)
  346. V1 = attac ;
  347. 'SAUT' 10 'LIGNE' ;
  348. 'MESS' 'ABSCISSE DU POINT DE REATTACHEMENT :' V1 ;
  349. 'SI' GRAPH ;
  350. EV1 = 'EVOL' 'MANU' RV . 'INCO' . 'TI' RV . 'INCO' . 'POSI' ;
  351. 'DESS' EV1 'MIMA' 'GRIL' 'TITR' 'Localisation reattachement = f(t)' ;
  352. *
  353. * Vitesse
  354. un = RV . 'INCO' . 'UN' ;
  355. vun = 'VECT' UN 0.1 'UX' 'UY' 'JAUNE' ;
  356. trace vun DOMTOT CNT1 'TITR' 'Vitesse' ;
  357. *
  358. * Pression
  359. pe = 'EXCO' (RV . 'INCO' . 'PRESSION') 'PRES' 'SCAL' ;
  360. 'SI' ('EGA' kpres 'MSOMMET') ;
  361. mp = 'DOMA' $DOMTOT 'MMAIL' ;
  362. trace pe mp ('CONTOUR' mp) ;
  363. 'SINO' ;
  364. pn = 'ELNO' $DOMTOT ('KCHT' $DOMTOT 'SCAL' kpres 0. pe) kpres;<'s0an>
  365. trace pn domtot CNT1 NLI0 'TITR' 'Pression' ;
  366. 'FINSI' ;
  367. 'FINSI' ;
  368. *
  369. *========================
  370. * Test de non régression
  371. *========================
  372. *
  373. 'SI' COMPLET ;
  374. TEST = ('ABS' (V1 '-' 5.8297)) < 0.0005 ;
  375. 'SINON' ;
  376. TEST = ('ABS' (V1 '-' 5.4447)) < 0.0005 ;
  377. 'FINSI' ;
  378. 'MESSAGE' 'Element Vitesse : ' DISCR ;
  379. 'MESSAGE' 'Element Pression : ' KPRES ;
  380. 'MESSAGE' 'Décentrement : ' KSUPG ;
  381. 'MESSAGE' 'Matrice masse : ' KMASS ;
  382. 'SI' TEST ;
  383. 'ERREUR' 0 ;
  384. 'SINON' ;
  385. 'ERREUR' 5 ;
  386. 'FINSI' ;
  387. *
  388. 'FIN' ;
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  

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