Télécharger la plaquette Cast3M


fr en

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. ************ 
  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.  
  348.  
  349.  
  350.  
  351.  
  352.  
  353.  

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