Télécharger la plaquette Cast3M


fr en

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

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