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

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