Télécharger la plaquette Cast3M


fr en

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

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