Télécharger la plaquette Cast3M


fr en

Télécharger back_pression_2.dgibi

Retour à la liste

Numérotation des lignes : 

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

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