Télécharger dvispassiLM.dgibi

Retour à la liste

Numérotation des lignes :

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

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