Télécharger la plaquette Cast3M


fr en

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

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