Télécharger la plaquette Cast3M


fr en

Télécharger gtkl.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  gtkl.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. **********************************************************
  5. * GRID TURBULENCE : voir aussi gridturb.dgibi            *
  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 2007                            *
  10. * Teste K-L 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. rv.'ALGO_KEPSILON'=mots 'KL';
  145.  
  146. RV = EQEX RV
  147.         'CLIM' 'UN' 'UIMP' INLET Uinlet
  148.                'UN' 'VIMP' INLET 0.
  149.                'UN' 'VIMP' COTE  0.
  150.         'CLIM' 'KN' 'TIMP' INLET Kinlet
  151.                'EN' 'TIMP' INLET Einlet
  152.  ;
  153.  
  154. RVP = EQEX 'OPTI' 'EF' KPRES
  155.         ZONE $DOMTOT 'OPER' KBBT  (-1.)     INCO 'UN' 'PRES'
  156.  CLIM 'PRES' TIMP poutlet 0.
  157. ;
  158.  
  159. rvp.'METHINV' . 'FCPRECT'=300 ;
  160.  
  161.  
  162.  
  163. RV.INCO = TABLE INCO ;
  164. RV.INCO.'UN' = KCHT $DOMTOT VECT SOMMET (uinlet 0.) ;
  165. RV.INCO.'KN' = KCHT $DOMTOT SCAL SOMMET Kinlet;
  166. RV.INCO.'EN' = KCHT $DOMTOT SCAL SOMMET 1.e-3   ;
  167. RV.INCO.'ENM1' = KCHT $DOMTOT SCAL SOMMET Einlet ;
  168. RV.INCO.'KNM1' = KCHT $DOMTOT SCAL SOMMET Kinlet ;
  169. RV.INCO.'MUF' = KCHT $DOMTOT SCAL SOMMET 1.E-10 ;
  170. RV.INCO.'Ro' = KCHT $DOMTOT SCAL SOMMET 1. ;
  171. RV.INCO.'Mu' = KCHT $DOMTOT SCAL SOMMET Nu ;
  172. RV.INCO.'PRES' = KCHT $DOMTOT SCAL KPRES 0. ;
  173. RV.INCO.'IT' = PROG 1 ;
  174. RV.INCO.'ERK' = PROG 0.;
  175. RV.INCO.'ERE' = PROG 0.;
  176. RV.INCO.'EAN'=EAN;
  177. RV.INCO.'KAN'=KAN;
  178. rv.inco.'Echl'=KAN**1.5*(inve EAN);
  179.  
  180.  RV.'PROJ' = RVP ;
  181.  
  182. EVOLKA = EVOL (BLEU) 'CHPO' KAN 'SCAL' A1A2 ;
  183. EVOLEA = EVOL (BLEU) 'CHPO' EAN 'SCAL' A1A2 ;
  184.  
  185. NCFL = DIME LCFL ;
  186. REPETER BCFL NCFL;
  187. CFL = extr LCFL &BCFL ;
  188. dt = CFL*deltax/Uinlet ;
  189. mess ' CFL=' CFL ' deltax ' deltax ' dt = ' dt ;
  190. RV.INCO.'DT' = dt;
  191.  
  192. RV.INCO.'IT' = PROG 1 ;
  193. RV.INCO.'ERK' = PROG 0.;
  194. RV.INCO.'ERE' = PROG 0.;
  195. RV.'PASDETPS'.'NUPASDT'= 1;
  196.  
  197. EXEC RV ;
  198.  
  199. EVOLKA = EVOLKA et (EVOL (ROUG) 'CHPO' RV.INCO.'KN' 'SCAL' A1A2) ;
  200. EVOLEA = EVOLEA et (EVOL (ROUG) 'CHPO' RV.INCO.'EN' 'SCAL' A1A2) ;
  201.  
  202. EVOLK = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
  203.            (RV.INCO.'ERK') ;
  204. EVOLE = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
  205.            (RV.INCO.'ERE') ;
  206. Si (EGA &BCFL 1);
  207. EVOLER = EVOLK et EVOLE;
  208. Sinon ;
  209. EVOLER = EVOLK et EVOLE et EVOLER ;
  210. Finsi ;
  211. * TESTS D'ERREUR *
  212.  
  213. erkm= extr &bcfl lerkm;
  214. erem= extr &bcfl lerem;
  215. list rv.inco.'ERK' ;
  216. mess ' Sur K' (MINI RV.INCO.'ERK') erkm;
  217. list rv.inco.'ERE' ;
  218. mess ' Sur E' (MINI RV.INCO.'ERE') erem;
  219. SI ( (MINI RV.INCO.'ERE') > erem ) ;
  220. *        ERREUR 5 ;
  221. FINSI ;
  222. SI ( (MINI RV.INCO.'ERK') > erkm ) ;
  223. *        ERREUR 5 ;
  224. FINSI ;
  225. FIN BCFL ;
  226.  
  227. SI ( GRAPH ) ;
  228. TRACE DOMTOT ;
  229. TAB1 = TABLE ;
  230. TAB1.1 = 'TIRC MARQ TRIA REGU' ;
  231. TAB1.2 = 'TIRC MARQ TRIA REGU' ;
  232. TAB1.3 = 'TIRL MARQ LOSA REGU' ;
  233. TAB1.4 = 'TIRL MARQ LOSA REGU' ;
  234. TAB1.5 = 'TIRM MARQ TRIB REGU' ;
  235. TAB1.6 = 'TIRM MARQ TRIB REGU' ;
  236. TAB1.7 = 'TIRR MARQ ETOI REGU' ;
  237. TAB1.8 = 'TIRR MARQ ETOI REGU' ;
  238.  
  239. TAB1.'TITRE' = TABLE ;
  240. REPETER BCL NCFL ;
  241. TAB1.'TITRE'. (2*&bcl-1) = (CHAI 'Cv_ k : CFL=' (EXTR LCFL &bcl));
  242. TAB1.'TITRE'. (2*&bcl)   = (CHAI 'Cv_ e : CFL=' (EXTR LCFL &bcl));
  243. FIN BCL ;
  244.  
  245. titl = chai 'Convergence Linf : Cfl = ' cfl ;
  246. DESS EVOLER 'XBOR' 0. (2.*NITMA) 'YBOR' -15.0 0.0 'GRIL'
  247.  'TITR' titl  LEGE TAB1 ;
  248.  
  249. TRACE RV.INCO.'KN' DOMTOT CNT 'TITR' 'KN' 20 ;
  250. TRACE RV.INCO.'EN' DOMTOT CNT 'TITR' 'EN' 20 ;
  251.  
  252.  
  253. TAB1 = TABLE ;
  254. TAB1.1 = 'MARQ CROI ' ;
  255. TAB1.2 = 'TIRC MARQ TRIA REGU' ;
  256. TAB1.3 = 'TIRL MARQ CARR REGU' ;
  257. TAB1.4 = 'TIRM MARQ LOSA REGU' ;
  258. TAB1.5 = 'TIRM MARQ ETOI REGU' ;
  259. TAB1.'TITRE' = TABLE ;
  260. TAB1.'TITRE'. 1 = 'MOT' 'k : ANALYTIQUE' ;
  261. REPETER BCL NCFL ;
  262. TAB1.'TITRE'. (1+&bcl) = (CHAI 'k : CFL=' (EXTR LCFL &bcl));
  263. FIN BCL ;
  264.  
  265. mess ' CFL=' CFL ;
  266. titk = chai 'Grid Turbulence K : Cfl = ' cfl ;
  267. tite = chai 'Grid Turbulence E : Cfl = ' cfl ;
  268. DESS EVOLKA LEGE TAB1 'TITR' titk 'GRIL';
  269.  
  270. TAB1.'TITRE'. 1 = 'MOT' 'epsilon ANALYTIQUE' ;
  271. REPETER BCL NCFL ;
  272. TAB1.'TITRE'. (1+&bcl) = (CHAI 'e : CFL=' (EXTR LCFL &bcl));
  273. FIN BCL ;
  274. DESS EVOLEA LEGE TAB1 'TITR' tite 'GRIL';
  275. FINSI ;
  276.  
  277. ERRK = KOPS RV.INCO.'KN' '-' KAN ;
  278. ERRK = KOPS ERRK '*' ERRK ;
  279. ERR2 = 0. ;
  280. REPETER BLOC1 (NBNO (DOMA $DOMTOT MAILLAGE)) ;
  281. P1 = (doma $DOMTOT SOMMET) POIN &BLOC1 ;
  282. ERR1 = 'EXTR' ERRK 'SCAL' P1 ;
  283. ERR2 = ERR2 + ERR1 ;
  284. FIN BLOC1 ;
  285. ERR2 = ERR2/(nbno (doma $DOMTOT MAILLAGE)) ;
  286. ERR2 = ERR2 ** 0.5 ;
  287. ERR2K = ERR2/Kinlet ;
  288. MESSAGE 'ERREUR SUR K EN NORME L2 = ' ERR2K ;
  289.  
  290. ERRE = KOPS RV.INCO.'EN' '-' EAN ;
  291. ERRE = KOPS ERRE '*' ERRE ;
  292. ERR2 = 0. ;
  293. REPETER BLOC1 (NBNO (DOMA $DOMTOT MAILLAGE))  ;
  294. P1 = (DOMA $DOMTOT SOMMET) POIN &BLOC1 ;
  295. ERR1 = 'EXTR' ERRE 'SCAL' P1 ;
  296. ERR2 = ERR2 + ERR1 ;
  297. FIN BLOC1 ;
  298. ERR2 = ERR2/(NBNO (DOMA $DOMTOT MAILLAGE)) ;
  299. ERR2 = ERR2 ** 0.5 ;
  300. ERR2E = ERR2/Einlet ;
  301. MESSAGE 'ERREUR SUR E EN NORME L2 = ' ERR2E ;
  302.  
  303.  
  304.  FIN ;
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  

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