Télécharger la plaquette Cast3M


fr en

Télécharger dvispqt.dgibi

Retour à la liste

Numérotation des lignes : 

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

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