Télécharger la plaquette Cast3M


fr en

Télécharger gridturb.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  gridturb.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. **********************************************************
  5. * GRID TURBULENCE : convection of homogeneous turbulence *
  6. * Analysis of the K-Epsilon TURBULENCE MODEL             *
  7. * Mohammadi/Pironneau p. 74     (Wiley)                  *
  8. * H. PAILLERE/TTMF/AVRIL 1997 (à l'origine 1/2 explicite)*
  9. * Version implicite MARS 2000                            *
  10. * Teste K-Epsilon en implicite méthode de projection     *
  11. **********************************************************
  12. * WARNING : Il faut que Nt/Nu > 1
  13.  
  14. *opti trace 'PSC';
  15. GRAPH = VRAI ;
  16. GRAPH = FAUX ;
  17. KPRES='MSOMMET';
  18. DISCR='LINE';
  19. KSUPG='SUPGDC';
  20. *
  21.         D = 1 ;
  22.         NITMA =200 ;
  23.         LCFL = PROG 0.1 0.5 1. 10.;
  24.         LCFL = PROG 1. 10. 50. 100.;
  25.         lerkm= prog -1.5 -2. -4. -4. ;
  26.         lerem= prog -1.5 -2. -4. -4. ;
  27.  
  28. NX = D*40 ;
  29. NY = D ;
  30.  
  32. opti elem qua8 ;
  33. opti isov suli ;
  34.  
  35. ******************************************************************
  36. * PROCEDURE POUR ESTIMER LA CONVERGENCE VERS L'ETAT STATIONNAIRE *
  37. ******************************************************************
  38.  
  39. DEBPROC CALCUL ;
  40. ARGU RVX*'TABLE' ;
  41.  
  42. RV = RVX.'EQEX' ;
  43.  
  44. DD = RV.PASDETPS.'NUPASDT' ;
  45. NN = DD/10 ;
  46. LO = (DD-(10*NN)) EGA 0 ;
  47. SI ( LO ) ;
  48.  
  49. EN   = RV.INCO.'EN' ;
  50. ENM1 = RV.INCO.'ENM1' ;
  51. KN   = RV.INCO.'KN' ;
  52. KNM1 = RV.INCO.'KNM1' ;
  53.  
  54. ERRE = KOPS EN '-' ENM1 ;
  55. ELIE = MAXI ERRE 'ABS' ;
  56. ELIE = (LOG (ELIE + 1.0E-20))/(LOG 20.) ;
  57. ERRK = KOPS KN '-' KNM1 ;
  58. ELIK = MAXI ERRK 'ABS' ;
  59. ELIK = (LOG (ELIK + 1.0E-20))/(LOG 20.) ;
  60. MESSAGE 'ITER ' RV.PASDETPS.'NUPASDT' '   ERREUR LINF ' ELIK ELIE ;
  61. RV.INCO.'ENM1' = KCHT $DOMTOT 'SCAL' 'SOMMET' (RV.INCO.'EN') ;
  62. RV.INCO.'KNM1' = KCHT $DOMTOT 'SCAL' 'SOMMET' (RV.INCO.'KN') ;
  63. IT = PROG RV.PASDETPS.'NUPASDT' ;
  64. ERK = PROG ELIK ;
  65. ERE = PROG ELIE ;
  66. RV.INCO.'IT' = (RV.INCO.'IT') ET IT ;
  67. RV.INCO.'ERK' = (RV.INCO.'ERK') ET ERK ;
  68. RV.INCO.'ERE' = (RV.INCO.'ERE') ET ERE ;
  69. FINSI ;
  70.  as2 ama1  = 'KOPS' 'MATRIK' ;
  71.  FINPROC as2 ama1 ;
  72.  
  73. ***********************
  74. * GEOMETRY OF PROBLEM *
  75. ***********************
  76.  
  77. A1 = 0.0 0.0 ;
  78. A2 = 1.0 0.0 ;
  79. A3 = 1.0 0.02;
  80. A4 = 0.0 0.02;
  81.  
  82. A1A2 = A1 'DROI' NX A2 ;
  83. A2A3 = A2 'DROI' NY A3 ;
  84. A3A4 = A3 'DROI' NX A4 ;
  85. A4A1 = A4 'DROI' NY A1 ;
  86.  
  87. DOMTOT = 'DALLER' A1A2 A2A3 A3A4 A4A1 'PLAN' ;
  88. DOMTOT = CHAN DOMTOT QUAF ;
  89. A4A1   = CHAN A4A1   QUAF ;
  90. A2A3   = CHAN A2A3   QUAF ;
  91. A1A2   = CHAN A1A2   QUAF ;
  92. A3A4   = CHAN A3A4   QUAF ;
  93. ELIM (DOMTOT et A4A1 et A2A3 et A3A4 et A1A2) 1.E-9 ;
  94. MCOTE  = A1A2 et A3A4;
  95.  
  96. $DOMTOT = 'MODE' DOMTOT 'NAVIER_STOKES' DISCR ;
  97. DOMA $DOMTOT IMPR ;
  98. $INLET  = 'MODE' A4A1 'NAVIER_STOKES' DISCR ;
  99.  INLET  = 'DOMA' $INLET 'MAILLAGE' ;
  100. $OUTLET = 'MODE' A2A3 'NAVIER_STOKES' DISCR ;
  101.  OUTLET = 'DOMA' $OUTLET 'MAILLAGE' ;
  102.  poutlet= 'DOMA' $OUTLET 'MSOMMET'  ;
  103. $COTE   = 'MODE' MCOTE 'NAVIER_STOKES' DISCR ;
  104. COTE    = DOMA $COTE MAILLAGE ;
  105.  
  106. DOMTOT = DOMA $DOMTOT MAILLAGE ;
  107. CNT = CONT DOMTOT ;
  108.  
  109. *****************
  110. * PHYSICAL DATA *
  111. *****************
  112.  
  113. Uinlet = 1.71;
  114. Kinlet = 0.013;
  115. Einlet = 0.029;
  116. deltax = (coor 1 A2)/nx ;
  117. Nu = 1.E-3 ;
  118.  
  119. *****************************************
  120. * ANALYTICAL SOLUTION FOR K(x) AND E(x) *
  121. *****************************************
  122.  
  123. C2 = 1.92 ;
  124. X1 = COOR 1 (DOMA $DOMTOT SOMMET) ;
  125. TINLET = KINLET/EINLET;
  126. TAN=(((C2 - 1.)*X1) + (UINLET*TINLET))/UINLET;
  127. M=1. - C2 ;
  128. KAN=KINLET * (( 1. - ((M*X1)/UINLET/TINLET))**(1./M));
  129. EAN=KAN*(INVE TAN) ;
  130.  
  131.  
  132. ********************
  133. * SET OF EQUATIONS *
  134. ********************
  135.  
  136. RV = EQEX $DOMTOT 'ITMA' NITMA  'OMEGA' 1. 'NITER' 1 'FIDT' 20
  137.         'ZONE' $DOMTOT 'OPER' CALCUL
  138.         'OPTI' KSUPG 'IMPL' EF
  139.         'ZONE' $DOMTOT 'OPER' 'NS' 'Ro' 'UN' 'MUF' 'INCO' 'UN'
  140.  'ZONE' $DOMTOT 'OPER' 'KEPSILON' 'Ro' 'UN' 'Mu' 'DT' 'INCO' 'KN' 'EN'
  141.  OPTI EFM1 'CENTREE'
  142.  'ZONE' $DOMTOT 'OPER' 'DFDT' 'Ro' 'UN' 'DT'            'INCO' 'UN'
  143.          ;
  144.  
  145. RV = EQEX RV
  146.         'CLIM' 'UN' 'UIMP' INLET Uinlet
  147.                'UN' 'VIMP' INLET 0.
  148.                'UN' 'VIMP' COTE  0.
  149.         'CLIM' 'KN' 'TIMP' INLET Kinlet
  150.                'EN' 'TIMP' INLET Einlet
  151.  ;
  152.  
  153. RVP = EQEX 'OPTI' 'EF' KPRES
  154.         ZONE $DOMTOT 'OPER' KBBT  (-1.)     INCO 'UN' 'PRES'
  155.  CLIM 'PRES' TIMP poutlet 0.
  156. ;
  157.  
  158. rvp.'METHINV' . 'FCPRECT'=300 ;
  159.  
  160.  
  161.  
  162. RV.INCO = TABLE INCO ;
  163. RV.INCO.'UN' = KCHT $DOMTOT VECT SOMMET (uinlet 0.) ;
  164. RV.INCO.'KN' = KCHT $DOMTOT SCAL SOMMET Kinlet;
  165. RV.INCO.'EN' = KCHT $DOMTOT SCAL SOMMET 1.e-3   ;
  166. RV.INCO.'ENM1' = KCHT $DOMTOT SCAL SOMMET Einlet ;
  167. RV.INCO.'KNM1' = KCHT $DOMTOT SCAL SOMMET Kinlet ;
  168. RV.INCO.'MUF' = KCHT $DOMTOT SCAL SOMMET 1.E-10 ;
  169. RV.INCO.'Ro' = KCHT $DOMTOT SCAL SOMMET 1. ;
  170. RV.INCO.'Mu' = KCHT $DOMTOT SCAL SOMMET Nu ;
  171. RV.INCO.'PRES' = KCHT $DOMTOT SCAL KPRES 0. ;
  172. RV.INCO.'IT' = PROG 1 ;
  173. RV.INCO.'ERK' = PROG 0.;
  174. RV.INCO.'ERE' = PROG 0.;
  175. RV.INCO.'EAN'=EAN;
  176. RV.INCO.'KAN'=KAN;
  177.  
  178.  RV.'PROJ' = RVP ;
  179.  
  180. EVOLKA = EVOL (BLEU) 'CHPO' KAN 'SCAL' A1A2 ;
  181. EVOLEA = EVOL (BLEU) 'CHPO' EAN 'SCAL' A1A2 ;
  182.  
  183. NCFL = DIME LCFL ;
  184. REPETER BCFL NCFL;
  185. CFL = extr LCFL &BCFL ;
  186. dt = CFL*deltax/Uinlet ;
  187. mess ' CFL=' CFL ' deltax ' deltax ' dt = ' dt ;
  188. RV.INCO.'DT' = dt;
  189.  
  190. RV.INCO.'IT' = PROG 1 ;
  191. RV.INCO.'ERK' = PROG 0.;
  192. RV.INCO.'ERE' = PROG 0.;
  193. RV.'PASDETPS'.'NUPASDT'= 1;
  194.  
  195. EXEC RV ;
  196.  
  197. EVOLKA = EVOLKA et (EVOL (ROUG) 'CHPO' RV.INCO.'KN' 'SCAL' A1A2) ;
  198. EVOLEA = EVOLEA et (EVOL (ROUG) 'CHPO' RV.INCO.'EN' 'SCAL' A1A2) ;
  199.  
  200. EVOLK = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
  201.            (RV.INCO.'ERK') ;
  202. EVOLE = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
  203.            (RV.INCO.'ERE') ;
  204. Si (EGA &BCFL 1);
  205. EVOLER = EVOLK et EVOLE;
  206. Sinon ;
  207. EVOLER = EVOLK et EVOLE et EVOLER ;
  208. Finsi ;
  209. * TESTS D'ERREUR *
  210.  
  211. erkm= extr &bcfl lerkm;
  212. erem= extr &bcfl lerem;
  213. list rv.inco.'ERK' ;
  214. mess ' Sur K' (MINI RV.INCO.'ERK') erkm;
  215. list rv.inco.'ERE' ;
  216. mess ' Sur E' (MINI RV.INCO.'ERE') erem;
  217. SI ( (MINI RV.INCO.'ERE') > erem ) ;
  218. *        ERREUR 5 ;
  219. FINSI ;
  220. SI ( (MINI RV.INCO.'ERK') > erkm ) ;
  221. *        ERREUR 5 ;
  222. FINSI ;
  223. FIN BCFL ;
  224.  
  225. SI ( GRAPH ) ;
  226. TRACE DOMTOT ;
  227. TAB1 = TABLE ;
  228. TAB1.1 = 'TIRC MARQ TRIA REGU' ;
  229. TAB1.2 = 'TIRC MARQ TRIA REGU' ;
  230. TAB1.3 = 'TIRL MARQ LOSA REGU' ;
  231. TAB1.4 = 'TIRL MARQ LOSA REGU' ;
  232. TAB1.5 = 'TIRM MARQ TRIB REGU' ;
  233. TAB1.6 = 'TIRM MARQ TRIB REGU' ;
  234. TAB1.7 = 'TIRR MARQ ETOI REGU' ;
  235. TAB1.8 = 'TIRR MARQ ETOI REGU' ;
  236.  
  237. TAB1.'TITRE' = TABLE ;
  238. REPETER BCL NCFL ;
  239. TAB1.'TITRE'. (2*&bcl-1) = (CHAI 'Cv_ k : CFL=' (EXTR LCFL &bcl));
  240. TAB1.'TITRE'. (2*&bcl)   = (CHAI 'Cv_ e : CFL=' (EXTR LCFL &bcl));
  241. FIN BCL ;
  242.  
  243. titl = chai 'Convergence Linf : Cfl = ' cfl ;
  244. DESS EVOLER 'XBOR' 0. (2.*NITMA) 'YBOR' -15.0 0.0 'GRIL'
  245.  'TITR' titl  LEGE TAB1 ;
  246.  
  247. TRACE RV.INCO.'KN' DOMTOT CNT 'TITR' 'KN' 20 ;
  248. TRACE RV.INCO.'EN' DOMTOT CNT 'TITR' 'EN' 20 ;
  249.  
  250.  
  251. TAB1 = TABLE ;
  252. TAB1.1 = 'MARQ CROI ' ;
  253. TAB1.2 = 'TIRC MARQ TRIA REGU' ;
  254. TAB1.3 = 'TIRL MARQ CARR REGU' ;
  255. TAB1.4 = 'TIRM MARQ LOSA REGU' ;
  256. TAB1.5 = 'TIRM MARQ ETOI REGU' ;
  257. TAB1.'TITRE' = TABLE ;
  258. TAB1.'TITRE'. 1 = 'MOT' 'k : ANALYTIQUE' ;
  259. REPETER BCL NCFL ;
  260. TAB1.'TITRE'. (1+&bcl) = (CHAI 'k : CFL=' (EXTR LCFL &bcl));
  261. FIN BCL ;
  262.  
  263. mess ' CFL=' CFL ;
  264. titk = chai 'Grid Turbulence K : Cfl = ' cfl ;
  265. tite = chai 'Grid Turbulence E : Cfl = ' cfl ;
  266. DESS EVOLKA LEGE TAB1 'TITR' titk 'GRIL';
  267.  
  268. TAB1.'TITRE'. 1 = 'MOT' 'epsilon ANALYTIQUE' ;
  269. REPETER BCL NCFL ;
  270. TAB1.'TITRE'. (1+&bcl) = (CHAI 'e : CFL=' (EXTR LCFL &bcl));
  271. FIN BCL ;
  272. DESS EVOLEA LEGE TAB1 'TITR' tite 'GRIL';
  273. FINSI ;
  274.  
  275. ERRK = KOPS RV.INCO.'KN' '-' KAN ;
  276. ERRK = KOPS ERRK '*' ERRK ;
  277. ERR2 = 0. ;
  278. REPETER BLOC1 (NBNO (DOMA $DOMTOT MAILLAGE)) ;
  279. P1 = (doma $DOMTOT SOMMET) POIN &BLOC1 ;
  280. ERR1 = 'EXTR' ERRK 'SCAL' P1 ;
  281. ERR2 = ERR2 + ERR1 ;
  282. FIN BLOC1 ;
  283. ERR2 = ERR2/(nbno (doma $DOMTOT MAILLAGE)) ;
  284. ERR2 = ERR2 ** 0.5 ;
  285. ERR2K = ERR2/Kinlet ;
  286. MESSAGE 'ERREUR SUR K EN NORME L2 = ' ERR2K ;
  287.  
  288. ERRE = KOPS RV.INCO.'EN' '-' EAN ;
  289. ERRE = KOPS ERRE '*' ERRE ;
  290. ERR2 = 0. ;
  291. REPETER BLOC1 (NBNO (DOMA $DOMTOT MAILLAGE))  ;
  292. P1 = (DOMA $DOMTOT SOMMET) POIN &BLOC1 ;
  293. ERR1 = 'EXTR' ERRE 'SCAL' P1 ;
  294. ERR2 = ERR2 + ERR1 ;
  295. FIN BLOC1 ;
  296. ERR2 = ERR2/(NBNO (DOMA $DOMTOT MAILLAGE)) ;
  297. ERR2 = ERR2 ** 0.5 ;
  298. ERR2E = ERR2/Einlet ;
  299. MESSAGE 'ERREUR SUR E EN NORME L2 = ' ERR2E ;
  300.  
  301.  
  302. FIN ;
  303.  
  304.  
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  

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