Télécharger dvisp.dgibi

Retour à la liste

Numérotation des lignes :

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

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