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

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