Télécharger dvispassi2.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : dvispassi.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. * CAVITE CARREE ­ VAHL­DAVIS Methode de projection implicite *
  5. * Malvina Renesson Aout 1999 *
  6. ***************************************************************
  7. GRAPH=VRAI ;
  8. GRAPH=FAUX ;
  9. COMPLET=VRAI;
  10. COMPLET=FAUX ;
  11. TCPT = FAUX;
  12. TKPR = VRAI;
  13. TRESOU = VRAI;
  14. IMPARA = FAUX;
  15.  
  16. Si (NON Complet);
  17. itmax= 80 ;
  18. DISCR='LINE' ;
  19. *DISCR='MACRO';
  20. DISCR='QUAF';
  21. d1 = 0.04; d2 = 0.12;
  22. d1 = 0.02; d2 = 0.06;
  23. DT=.25;
  24. Ra = 1.e6 ;
  25. Sinon ;
  26. *Ra = 1.e5 ;
  27. *itmax= 600 ;
  28. Ra = 1.e6 ;
  29. itmax= 600 ;
  30. DISCR='LINE' ;
  31. d1 = 0.01; d2 = 0.07;
  32. DT=.25;
  33.  
  34. Finsi ;
  35.  
  36. KSUPG ='CENTREE' ;
  37. KPRESS='MSOMMET' ;
  38. KPRESS='CENTRE';
  39. KPRESS='CENTREP1';
  40.  
  41. ******************************
  42. *ESTIMATION DE LA CONVERGENCE*
  43. ******************************
  44. DEBPROC CALCUL;
  45. ARGU RVX*'TABLE' MDZ/'MMODEL' MDZ/'TABLE';
  46. RV = RVX.'EQEX';
  47. DD = RV.PASDETPS.'NUPASDT';
  48. NN = DD/5;
  49. LO = (DD-(5*NN)) EGA 0;
  50. SI (LO);
  51. UN = RV.INCO.'UN';
  52. UNM1 = RV.INCO.'UNM1';
  53. unx = kcht $mt scal sommet (exco 'UX' un);
  54. unm1x = kcht $mt scal sommet (exco 'UX' unm1);
  55. uny = kcht $mt scal sommet (exco 'UY' un);
  56. unm1y = kcht $mt scal sommet (exco 'UY' unm1);
  57. ERRX = KOPS unx '-' unm1x;
  58. ERRY = KOPS uny '-' unm1y;
  59. ELIX = MAXI ERRX 'ABS';
  60. ELIY = MAXI ERRY 'ABS';
  61. ELIX = (LOG (ELIX + 1.0E-10))/(LOG 10.);
  62. ELIY = (LOG (ELIY + 1.0E-10))/(LOG 10.);
  63. MESSAGE 'ITER' RV.PASDETPS.'NUPASDT' 'ERREUR LINF' ELIX ELIY;
  64. RV.INCO.'UNM1' = KCHT $MT 'VECT' 'SOMMET' (RV.INCO.'UN');
  65. IT = PROG RV.PASDETPS.'NUPASDT';
  66. ER = PROG ELIY;
  67. RV.INCO.'IT' = (RV.INCO.'IT') ET IT;
  68. RV.INCO.'ER' = (RV.INCO.'ER') ET ER;
  69. FINSI;
  70. as2 ama1 = 'KOPS' 'MATRIK' ;
  71. 'FINPROC' as2 ama1 ;
  72.  
  73. ************
  74. * MAILLAGE *
  75. ************
  76. OPTI ELEM QUA8;
  77. crit = 0.0000001;
  78. xxll = 1.;
  79. p1 = 0. 0.;
  80. p15 = (xxll/2.) 0.;
  81. p2 = xxll 0.;
  82. p25 = xxll 0.5;
  83. p3 = xxll 1.;
  84. p35 = (xxll/2.) 1.;
  85. p4 = 0. 1.;
  86. p45 = 0. 0.5;
  87.  
  88. bas1 = p1 d dini d1 dfin d2 p15;
  89. bas2 = p15 d dini d2 dfin d1 p2;
  90. bas = (bas1 et bas2);elim crit bas;
  91. cdro = p2 d dini d1 dfin d2 p25 d dini d2 dfin d1 p3;
  92. haut = p3 d dini d1 dfin d2 p35 d dini d2 dfin d1 p4;
  93. cgau = p4 d dini d1 dfin d2 p45 d dini d2 dfin d1 p1;
  94. coupx = bas plus (0. 0.5) ;
  95. coupy=cgau plus (0.5 0.) ;
  96. cnt = bas et cdro et haut et cgau;
  97. mt = bas cdro haut cgau daller;
  98. bord = cont mt;
  99. *mt = surf bord;
  100. elim (mt et coupx) 1.e-5;
  101. Si GRAPH ; trace mt ; Finsi ;
  102. *option donn 5;
  103. orienter mt;
  104. ********
  105. * DOMA *
  106. ********
  107.  
  108. Mmt = chan QUAF mt ;
  109. $mt = MODE Mmt 'NAVIER_STOKES' DISCR ;
  110. mt = doma $mt maillage;
  111.  
  112. Mcgau=chan cgau QUAF ;
  113. $cgau=mode Mcgau 'NAVIER_STOKES' DISCR ;
  114.  
  115. Mcdro=chan cdro QUAF ;
  116. $cdro=mode Mcdro 'NAVIER_STOKES' DISCR ;
  117.  
  118. Mcnt=chan cnt QUAF ;
  119. $cnt=mode Mcnt 'NAVIER_STOKES' DISCR ;
  120.  
  121. Mcoupx = chan coupx 'QUAF' ;
  122. $coupx=mode Mcoupx 'NAVIER_STOKES' DISCR ;
  123.  
  124. Mcoupy = chan coupy 'QUAF' ;
  125. $coupy=mode Mcoupy 'NAVIER_STOKES' DISCR ;
  126.  
  127. Elim (Mmt et Mcgau et Mcdro et Mcnt et Mcoupx
  128. et Mcoupy) 1.e-5 ;
  129. cgau=doma $cgau maillage ;
  130. cdro=doma $cdro maillage ;
  131. cnt=doma $cnt maillage ;
  132. coupx=doma $coupx maillage ;
  133. coupy=doma $coupy maillage ;
  134.  
  135. *option donn 5;
  136. **************
  137. * PARAMETRES *
  138. **************
  139. Pr = 0.71;
  140. Gr = Ra/Pr;
  141. NU = 1/(Gr**0.5);
  142. ALF= NU/Pr;
  143. gb = kcht $mt vect centre (0. -1.);
  144. uref = 0.2 ;
  145.  
  146. ***********************
  147. * CREATION DES TABLES *
  148. ***********************
  149. RV = EQEX 'OMEGA' 1. 'NITER' 1 ITMA itmax 'FIDT' 1
  150. 'OPTI' 'EF' 'IMPL' KSUPG
  151. 'ZONE' $MT 'OPER' CALCUL
  152. 'ZONE' $MT 'OPER' 'NS' 1. 'UN' NU 'GB' 'TN' 0.5 'INCO' 'UN'
  153. 'ZONE' $MT 'OPER' 'TSCAL' 1. 'UN' ALF 0. 'INCO' 'TN'
  154. 'OPTI' 'CENTREE'
  155. ZONE $mt OPER DFDT 1. 'UN' DT INCO 'UN'
  156. ZONE $mt OPER DFDT 1. 'TN' DT INCO 'TN'
  157. ;
  158. ****
  159.  
  160. RV = EQEX RV
  161. 'CLIM' 'UN' 'UIMP' cnt 0. 'UN' 'VIMP' cnt 0.
  162. 'TN' 'TIMP' cgau 1.
  163. 'TN' 'TIMP' cdro 0.;
  164.  
  165. rv.'METHINV'.TYPINV=1 ;
  166. rv.'METHINV'.IMPINV=0 ;
  167. rv.'METHINV'.NITMAX=400;
  168. rv.'METHINV'.PRECOND=3 ;
  169. rv.'METHINV'.RESID =1.e-10;
  170. rv. 'METHINV' . 'FCPRECT'=1 ;
  171. rv. 'METHINV' . 'FCPRECI'=1 ;
  172.  
  173.  
  174. ****
  175. * La cavité est fermée il faut imposer la pression en 1 point !
  176. prep1=doma $mt kpress;
  177. bcp=elem prep1 POI1 (lect 1 ) ;
  178.  
  179. RVP = EQEX
  180. 'OPTI' 'EF' KPRESS
  181. ZONE $mt OPER KBBT -1. INCO 'UN' 'PRES'
  182. CLIM PRES TIMP bcp 0. ;
  183.  
  184. rvp.'METHINV'.TYPINV=1 ;
  185. rvp.'METHINV'.IMPINV=0 ;
  186. rvp.'METHINV'.NITMAX=300;
  187. rvp.'METHINV'.PRECOND=3 ;
  188. rvp.'METHINV'.RESID =1.e-8 ;
  189. rvp.'METHINV' . 'FCPRECT'=100 ;
  190. rvp.'METHINV' . 'FCPRECI'=100 ;
  191.  
  192.  
  193.  
  194. RV.INCO = TABLE INCO;
  195. RV.INCO.'UN' = kcht $MT VECT SOMMET (0. 0.);
  196. *RV.INCO.'TN' = kcht $MT SCAL SOMMET (1.-corsx);
  197. RV.INCO.'TN' = kcht $MT SCAL SOMMET 0.5 ;
  198. RV.INCO.'PRES' = kcht $MT SCAL KPRESS 0. ;
  199. RV.INCO.'UNM1' = kcht $MT VECT SOMMET (1.E-5 1.E-5);
  200. RV.INCO.'IT' = PROG 1;
  201. RV.INCO.'ER' = PROG 0.;
  202. RV.INCO.'GB' = gb ;
  203. RV.'PROJ' = RVP;
  204.  
  205. RV.'TCPT' = TCPT ;
  206. RV.'TKPR' = TKPR ;
  207. RV.'TRESOU'= TRESOU;
  208. RV.'IMPARA'= IMPARA;
  209. temps;
  210. EXEC RV;
  211. temps ;
  212.  
  213. *************
  214. * RESULTATS *
  215. *************
  216. *option sauv 'SAV';
  217. *sauv Pr Ra Gr NU ALF gb uref RV;
  218. *************
  219. OPTI ISOV SULI;
  220.  
  221.  
  222. un=(RV.INCO.'UN') * (Pr*(Gr**0.5)) ;
  223. unch = vect un (0.3/(Pr*(Gr**0.5))) UX UY JAUNE;
  224.  
  225. u1y = exco un 'UX' ;
  226. u2x = exco un 'UY' ;
  227. evu=evol 'CHPO' u2x coupx ;
  228. evv=evol 'CHPO' u1y coupy ;
  229.  
  230. EVOL4 = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG---E---inf'
  231. (RV.INCO.'ER');
  232.  
  233. GRADT = KOPS RV.INCO.'TN' 'GRAD' $MT;
  234. DTDX = KCHT $MT 'SCAL' 'CENTRE' (EXCO 'UX' GRADT);
  235. DTDX = (-1)*(ELNO $MT DTDX);
  236.  
  237. EVOLd = EVOL 'CHPO' DTDX SCAL cdro ;
  238. EVOLg = EVOL 'CHPO' DTDX SCAL cgau ;
  239. evol1=evold et evolg ;
  240. LISTE EVOL1;
  241. TAB1 = TABLE;
  242. TAB1.1 = 'MARQ LOSA';
  243. TAB1.2 = 'MARQ CROI';
  244.  
  245. Si GRAPH ;
  246. DESS EVOL1 'XBOR' 0. 1. 'GRIL' TAB1 'TITRE' ' Nusselt' ;
  247. trac RV.INCO.'TN' mt (cont mt) 14 'TITRE' ' Isothermes' ;
  248. trace unch mt (cont mt) 'TITRE' ' Vitesses' ;
  249. dess evu 'TITRE' 'Coupe ox à y=1/2 (uy)' ;
  250. dess evv 'TITRE' 'Coupe oy à x=1/2 (ux)' ;
  251. dess EVOL4 'XBOR' 0. itmax 'YBOR' -10.0 0.0;
  252. Finsi;
  253.  
  254. Nusd = (redu DTDX cdro);
  255. Nusg = (redu DTDX cgau);
  256.  
  257. Numg = somt (Nusg*(doma $cgau 'XXDIAGSI')) ;
  258. Numd = somt (Nusd*(doma $cdro 'XXDIAGSI')) ;
  259.  
  260. un2=KOPS un '*' (Pr*(Gr**0.5));
  261. sw = KOPS un 'ROT' $mt ;
  262. rk = EQEX $mt 'OPTI' 'EF' 'IMPL'
  263. ZONE $mt OPER LAPN -1. INCO 'PSI'
  264. ZONE $mt OPER FIMP sw INCO 'PSI'
  265. 'CLIM' 'PSI' 'TIMP' cnt 0.;
  266. rk.'INCO'=table 'INCO' ;
  267. rk.'INCO'.'PSI'=kcht $mt scal sommet 0. ;
  268. EXec rk ;
  269. psi=rk.'INCO'.'PSI';
  270. psi2=kops psi '*' (Pr*(Gr**0.5));
  271.  
  272. Si GRAPH ;
  273. trac psi mt CNT 14 'TITRE' 'Fonction de courant ';
  274. Finsi;
  275.  
  276. opti 'ECHO' 0 ;
  277. MESS '----------- Ra = ' Ra '------------------------------------';
  278. MESS ' NUSSELT MOYEN NUSSELT MAX NUSSELT MIN' ;
  279. MESS 'CAST3M : Gauche ' NUMG ' ' (MAXI Nusg) ' ' (MINI Nusg) ;
  280. MESS 'CAST3M : Droite ' NUMD ' ' (MAXI Nusd) ' ' (MINI Nusd) ;
  281. Si (EGA Ra 1.e5);
  282. MESS 'Vahl Davis : Ra 1.e6 4.519 7.717 0.729';
  283. Finsi ;
  284. Si (EGA Ra 1.e6);
  285. MESS 'Vahl Davis : Ra 1.e6 8.8 17.925 0.989';
  286. Finsi ;
  287. MESS '---------------------------------------------------------------';
  288. MESS ' Umax Vmax ' ;
  289. MESS 'CAST3M : '
  290. (MAXI (exco un 'UX')) (MAXI (exco un 'UY'));
  291. Si (EGA Ra 1.e5);
  292. MESS 'Vahl Davis : Ra 1.e6 34.73 68.59 ' ;
  293. Finsi ;
  294. Si (EGA Ra 1.e6);
  295. MESS 'Vahl Davis : Ra 1.e6 64.63 219.36 ' ;
  296. Finsi ;
  297. MESS '---------------------------------------------------------------';
  298. MESS ' Psi max ' ;
  299. MESS 'CAST3M : ' (MAXI PSI) ;
  300. Si (EGA Ra 1.e5);
  301. MESS 'Vahl Davis : Ra 1.e6 9.612 ' ;
  302. Finsi;
  303. Si (EGA Ra 1.e6);
  304. MESS 'Vahl Davis : Ra 1.e6 16.75 ' ;
  305. Finsi;
  306. MESS '---------------------------------------------------------------';
  307. opti 'ECHO' 1 ;
  308.  
  309. Si (NON COMPLET) ;
  310.  
  311. ern=abs( NUMD - 7.835)/7.835 ; ernm=abs((MAXI Nusg) - 15.6)/15.6 ;
  312. eru=abs((MAXI (exco un 'UY')) - 238.5)/238.5;
  313. erp=abs((MAXI PSI) - 31.7)/31.7;
  314. mess 'ern = ' ern ' ernm = ' ernm ' eru = ' eru ' erp = ' erp ;
  315. Si (ern > 0.005); erreur 5 ; Finsi ;
  316. Si (ernm > 0.005); erreur 5 ; Finsi ;
  317. Si (eru > 0.005); erreur 5 ; Finsi ;
  318. Si (erp > 0.005); erreur 5 ; Finsi ;
  319.  
  320. Finsi ;
  321.  
  322. FIN ;
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  
  340.  
  341.  
  342.  
  343.  
  344.  
  345.  
  346.  
  347.  
  348.  

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