Télécharger la plaquette Cast3M


fr en

Télécharger coneq.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  coneq.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *=---------------------------------------------------------------------
  5. *= Transport d'un cone par un champ de vitesse à rotationnel constant.
  6. *= Comparaison de schéma en temps EF implicites
  7. *=---------------------------------------------------------------------
  8. *
  9. *------------------
  10. * Options générales
  11. *------------------
  12. *
  13. GRAPH = 'N'                   ;
  14. ERR1=5.e-5     ;
  15. COMPLET = FAUX ;
  16. SI ( COMPLET ) ;
  17.   DIAM=2. ;
  18. SINON        ;
  19.   DIAM=0.5 ;
  20. FINSI ;
  21.  
  22. DISCR='QUAF';
  23.  
  24. 'OPTI' 'DIME' 2 'ELEM' 'QUA4' ;
  25. 'OPTI' 'ISOV' 'SULI'          ;
  26. 'OPTI' 'ECHO' 1 'TRACER' 'PSC' ;
  27. *
  28. *=========
  29. * MAILLAGE  
  30. *=========  
  31. *
  32. *
  33. *- Création des points supports des DROITES
  34. *
  35. L     = 2.D0 ; LS2   = L / 2.D0        ;
  36. H     = L    ; HS2   = H / 2.D0        ;
  37. X0    = -1.D0 * LS2     ; X1    = X0 + L          ;
  38. Y0    = -1.D0 * HS2     ; Y1    = Y0 + H          ;
  39. INUMX = 20              ;
  40. INUMY = INUMX           ;
  41. XNUMY = 'FLOT' INUMX    ;
  42. INUX1 = INUMX - 1       ;
  43. INUY1 = INUMY - 1       ;
  44. X01   = X0 + X1 * 0.5D0 ;
  45. Y01   = Y0 + Y1 * 0.5D0 ;
  46. DX    = X1 - X0 / INUMX ;
  47. DY    = Y1 - Y0 / INUMY ;
  48. DX1   = DX / 2.D0       ;
  49. DY1   = DY / 2.D0       ;
  50. *
  51. B1 = X0 Y0  ;
  52. B3 = X1 Y0  ;
  53. H1 = X0 Y1  ;
  54. H3 = X1 Y1  ;
  55. *
  56. *
  57. *- Création des DROITES frontieres
  58. *
  59. DRBAS = B1 'DROI' INUMX B3 ;
  60. DRDRO = B3 'DROI' INUMY H3 ;
  61. DRHAU = H3 'DROI' INUMX H1 ;
  62. DRGAU = H1 'DROI' INUMY B1 ;
  63. DRCOT = DRBAS 'ET' DRDRO 'ET' DRGAU 'ET' DRHAU ;
  64. PELIM = DX1 / (5.D0 * INUMX) ;
  65. *
  66. *- Creation maillage GEOMETRIQUE
  67. *
  68. PTOT1 = 'DALL' DRBAS DRDRO DRHAU DRGAU ;
  69. *
  70. *- Creation maillage HYBRIDE y compris sous-objets (cond. limites)
  71. *
  72. PTOTC = CHAN PTOT1 QUAF ;
  73. MDRCOT = CHAN DRCOT QUAF ;
  74. ELIM (PTOTC et MDRCOT) 1.E-5 ;
  75.  
  76. HYTOT = 'MODE' PTOTC 'NAVIER_STOKES' DISCR;
  77. $DRCOT= 'MODE' MDRCOT 'NAVIER_STOKES' DISCR ;
  78.  
  79. DRCOT=doma $drcot 'MAILLAGE' ;
  80.  
  81. XVOL= DOMA HYTOT 'VOLUME'  ;
  82. DOMA HYTOT 'IMPR' ;
  83. ptot1=DOMA HYTOT  'MAILLAGE';
  84.  
  85. *
  86. *
  87. *- Coordonnées des points par classe de points
  88. *
  89. XS YS = 'COOR' ( DOMA HYTOT  'MAILLAGE') ;
  90. *
  91. *  XCEN : Abscisse du sommet du cone
  92. *  YCEN : Ordonnée du sommet du cone
  93. *  RCEN : Rayon du cone
  94. *
  95. XCEN = 0.D0       ;
  96. YCEN =-0.5D0      ;
  97. RCEN = 3.0D0 * DX ;
  98.  
  99. *================
  100. * INITIALISATIONS
  101. *================
  102. *
  103. *                                            ------------------------
  104. *                                            = Conditions initiales =
  105. *                                            ------------------------
  106. * h(x,y) = f(r)
  107. * avec f(r) = max( h(1.-r/rcen) , 0.)
  108. * On se sert des égalités suivantes :
  109. *      m = min(f,g) = (f+g - abs(f-g)) / 2
  110. *      M = max(f,g) = (f+g + abs(f-g)) / 2
  111. *
  112. HAUT = 1.D0 ;
  113. *
  114. *- Charge aux sommets
  115. *
  116. CRS1  = (XS-XCEN)*(XS-XCEN) + ((YS-YCEN)*(YS-YCEN)) ** 0.5D0 ;
  117. CRS2  = 1.D0 - (CRS1 / RCEN) * HAUT                          ;
  118. CRS3  = 'ABS' CRS2                                           ;
  119. CRS4  = CRS2 + CRS3 / 2.D0                                   ;
  120. SI ( COMPLET ) ;
  121. daxh=drbas plus (0. 0.5)  ;
  122. daxv=drgau plus (1. 0.)  ;
  123. elim (daxh et daxv et ptot1) 1.e-5 ;
  124. SINON ;
  125. daxh=drbas plus (0. 1.)  ;
  126. daxv=drgau plus (0.5 0.)  ;
  127. elim (daxh et daxv et ptot1) 1.e-5 ;
  128. FINSI ;
  129. *
  130. *- Vitesses aux sommets
  131. *
  132. SPEEX = 'NOMC' 'UX' YS 'NATU' 'DISCRET' ;
  133. SPEEY = 'NOMC' 'UY' XS 'NATU' 'DISCRET' ;
  134. SPEEY = (-1.D0) * SPEEY                 ;
  135.  
  136. TN    = KCHT HYTOT  SCAL SOMMET  CRS4        ;
  137. evolth=evol 'CHPO' TN  daxh ;
  138. evoltv=evol 'CHPO' TN  daxv ;
  139. s2c= SOMT ((NOEL HYTOT TN)*XVOL) ;
  140. mess ' Solution exacte projetée ' s2c ;
  141.  
  142.  
  143. DEBPROC TEST ;
  144. ARGU OPDFDT*MOT  OPKONV*MOT BETA*FLOTTANT CFL*FLOTTANT
  145.      DIAM*FLOTTANT GRAPH*MOT TIT*MOT;
  146. DT=CFL*0.1 ;
  147. NBIT=enti (3.146 * DIAM / DT) + 1 ;
  148.  
  149.  coef = 1.e-20 ;
  150.  
  151.  RV =  'EQEX' HYTOT ITMA NBIT 'FIDT' 1 'DTI' 0.
  152.  'OPTI' 'EF' 'SEMI' BETA  OPKONV 'MATCONS' 'MTCV'
  153.  'ZONE' HYTOT 'OPER' 'KONV' 1.  'UN' COEF DT 'INCO' 'TN'
  154.  'CLIM' 'TN' 'TIMP' DRCOT  0.  ;
  155.  
  156.  Si ( EGA OPDFDT 'BDF2');
  157.  RV =  'EQEX' RV
  158.  'OPTI' 'EF' 'CENTREE' BDF2
  159.  'ZONE' HYTOT 'OPER' 'DFDT' 1. 'TNM' 'TNMM' DT 'UN' coef 'INCO' 'TN'
  160.  ;
  161.  Sinon ;
  162.  RV =  'EQEX' RV
  163.  'OPTI' 'EF'     OPDFDT
  164.  'ZONE' HYTOT 'OPER' 'DFDT' 1. 'TNM' DT 'UN' coef 'INCO' 'TN'
  165.  ;
  166.  FINSI ;
  167.  
  168.  RV.'INCO'         = TABLE 'INCO'                        ;
  169.  RV.'INCO'.'TN'    = KCHT HYTOT  SCAL SOMMET  0.          ;
  170.  RV.'INCO'.'TNM'    = KCHT HYTOT  SCAL SOMMET  CRS4       ;
  171.  RV.'INCO'.'TNMM'    = KCHT HYTOT  SCAL SOMMET  CRS4      ;
  172.  RV.'INCO'.'T0'    = KCHT HYTOT  SCAL SOMMET  0.          ;
  173.  RV.'INCO'.'UN'    = KCHT HYTOT  VECT SOMMET  SPEEX  SPEEY ;
  174.  
  175. SI  ('EGA' graph 'O' );
  176. u = vect RV.'INCO'.'UN' 1. ux uy ; 
  177.  trac u ptot1 'NCLK';
  178.  trac RV.'INCO'.'TNM'  ptot1 'NCLK';
  179. Finsi ;
  180.  
  181. rv.'OMEGA'= 1. ;
  182.  rv.'NITER'=1 ;
  183.  
  184. exec rv ;
  185.  
  186. maxc=maxi (RV.'INCO'.'TN') ; minc=mini (RV.'INCO'.'TN');
  187. mess 'MAXC=' maxc 'MINC=' minc ;
  188. TN=(RV.'INCO'.'TN') ;
  189. s2c= SOMT ((NOEL HYTOT TN)*XVOL) ;
  190. mess ' S2C= ' s2c ;
  191.  
  192. SI  ('EGA' graph 'O' );
  193.    P1 = 'PROG' -5. -5. 2. ;
  194.    MONTAGNE (RV.'INCO'.'TN') ptot1 2. P1 'TITRE' tit 'CACHE' ;
  195.    'TRACER' RV.'INCO'.'TN'  ptot1 ('PROG' -0.1 pas 0.1 0.9) drcot
  196.            'TITRE' TIT 'NCLK' ;
  197. FINSI ;
  198. FINPROC RV ;
  199.  
  200. * SCHEMA CRANK-NICHOLSON GENERALISE GALERKIN (CNG)
  201.  
  202. CFL=0.3 ;
  203. BETA=1.  ;
  204. OPDFDT='CENTREE' ;
  205. OPKONV='CNG' ;
  206. TIT= CHAI 'CNG (1/4 de TOUR)' ;
  207. rv = TEST OPDFDT OPKONV BETA  CFL DIAM GRAPH TIT;
  208. evolh=evol 'CHPO' (RV.'INCO'.'TN')  daxh ;
  209. evolv=evol 'CHPO' (RV.'INCO'.'TN')  daxv ;
  210. evolth=evolth et evolh ;
  211. evoltv=evoltv et evolv ;
  212.  
  213. SI(NON COMPLET ) ;
  214. lrr=prog
  215.    0.0000     -2.55800E-04 -1.72779E-03 -2.61253E-03 -2.39863E-03
  216.  -7.29778E-04  1.29162E-03  2.72106E-02  0.33912      0.70209    
  217.   0.92622      0.65468      0.33257      3.10799E-02 -7.50800E-03
  218.  -2.72013E-03 -6.27409E-03 -4.72605E-03 -1.68431E-03  4.86326E-03
  219.    0.0000;
  220.   lr=extr evolv 'ORDO';
  221.   list lr ;
  222.   ER=SOMM( abs (lr - lrr) ) ;
  223.   mess ' Erreur sur CNG ' er ;
  224.   si ( er > 1.e-4) ; erreur 5 ; finsi ;
  225. FINSI ;
  226.  
  227. * SCHEMA Tenseur Visqueux (TVISQ)
  228.  
  229. CFL=0.2 ;
  230. BETA=0.  ;
  231. OPDFDT='CENTREE' ;
  232. OPKONV='TVISQUEU' ;
  233. TIT= CHAI 'TVISQUEU (1/4 de TOUR)' ;
  234. rv = TEST OPDFDT OPKONV BETA  CFL DIAM GRAPH TIT;
  235. evolh=evol 'CHPO' (RV.'INCO'.'TN')  daxh ;
  236. evolv=evol 'CHPO' (RV.'INCO'.'TN')  daxv ;
  237. evolth=evolth et evolh ;
  238. evoltv=evoltv et evolv ;
  239.  
  240. SI(NON COMPLET ) ;
  241. lrr=prog
  242.    0.0000     -6.50725E-04  1.84987E-04 -3.96571E-04  1.01095E-03
  243.  -8.97916E-04 -2.64323E-03  1.68876E-02  0.31908      0.66562    
  244.   0.92950      0.68354      0.35128      3.29776E-02  3.66575E-03
  245.  -2.57892E-04  5.48000E-04  1.06891E-03 -4.26246E-04 -8.88270E-04
  246.    0.0000;
  247.   lr=extr evolv 'ORDO';
  248.   list lr ;
  249.   ER=SOMM( abs (lr - lrr) ) ;
  250.   mess ' Erreur sur Tenseur Visqueux ' er ;
  251.   si ( er > 1.e-4 ) ; erreur 5 ; finsi ;
  252. FINSI ;
  253.  
  254. * SCHEMA CRANK-NICHOLSON GALERKIN (CN)
  255.  
  256. CFL=1. ;
  257. BETA=0.5 ;
  258. OPDFDT='CENTREE' ;
  259. OPKONV='CENTREE' ;
  260. TIT= CHAI 'CN  (1/4 de TOUR)' ;
  261. rv = TEST OPDFDT OPKONV BETA  CFL DIAM GRAPH TIT;
  262. evolh=evol 'CHPO' (RV.'INCO'.'TN')  daxh ;
  263. evolv=evol 'CHPO' (RV.'INCO'.'TN')  daxv ;
  264. evolth=evolth et evolh ;
  265. evoltv=evoltv et evolv ;
  266.  
  267. SI ( NON COMPLET ) ;
  268. lrr=prog
  269.    0.0000     -1.03945E-03 -2.15865E-05  1.46772E-03 -8.37973E-04
  270.  -5.20953E-04  9.73026E-03  5.68738E-02  0.24455      0.61662    
  271.   0.85196      0.80630      0.39733      7.72915E-02 -5.45626E-02
  272.   3.51786E-03  2.85817E-03  9.89801E-04 -5.25209E-03 -7.05692E-03
  273.    0.0000;
  274.   lr=extr evolv 'ORDO';
  275.   list lr ;
  276.   ER=SOMM( abs (lr - lrr) ) ;
  277.   mess ' Erreur sur CN ' er ;
  278.   si ( er > 1.e-4 ) ; erreur 5 ; finsi ;
  279. FINSI ;
  280.  
  281. * SCHEMA CRANK-NICHOLSON PETROV-GALERKIN (CN-PG)
  282.  
  283. CFL=1. ;
  284. BETA=0.5 ;
  285. OPDFDT='SUPG' ;
  286. OPKONV='SUPG' ;
  287. TIT= CHAI 'CN-PG  (1/4 de TOUR)' ;
  288. rv = TEST OPDFDT OPKONV BETA  CFL DIAM GRAPH TIT;
  289. evolh=evol 'CHPO' (RV.'INCO'.'TN')  daxh ;
  290. evolv=evol 'CHPO' (RV.'INCO'.'TN')  daxv ;
  291. evolth=evolth et evolh ;
  292. evoltv=evoltv et evolv ;
  293.  
  294. SI (NON COMPLET ) ;
  295. lrr=prog
  296.    0.0000      2.74117E-08  5.34498E-07  7.58818E-07  2.94414E-05
  297.   7.57305E-04  8.80454E-03  5.71758E-02  0.24644      0.61474    
  298.   0.84870      0.78601      0.39482      8.58106E-02 -4.81594E-02
  299.   2.41366E-03  9.25329E-04 -5.22345E-05 -2.55730E-04 -1.64667E-04
  300.    0.0000;
  301.   lr=extr evolv 'ORDO';
  302.   list lr ;
  303.   ER=SOMM( abs (lr - lrr) ) ;
  304.   mess ' Erreur sur CN-PG ' er ;
  305.   si ( er > 1.e-4 ) ; erreur 5 ; finsi ;
  306. FINSI ;
  307.  
  308. * SCHEMA EULER IMPLICITE PETOV-GALERKIN (IMPL-PG)
  309.  
  310. CFL=1. ;
  311. BETA=1. ;
  312. OPDFDT='SUPG' ;
  313. OPKONV='SUPG' ;
  314. TIT= CHAI 'IMPL-PG  (1/4 de TOUR)' ;
  315. rv = TEST OPDFDT OPKONV BETA  CFL DIAM GRAPH TIT;
  316. evolh=evol 'CHPO' (RV.'INCO'.'TN')  daxh ;
  317. evolv=evol 'CHPO' (RV.'INCO'.'TN')  daxv ;
  318. evolth=evolth et evolh ;
  319. evoltv=evoltv et evolv ;
  320.  
  321. SI ( NON COMPLET ) ;
  322. lrr=prog
  323.    0.0000     -9.04250E-06  6.03244E-05 -1.72068E-04  2.18717E-03
  324.   2.76778E-02  8.96770E-02  0.19189      0.32630      0.45198    
  325.   0.51491      0.48040      0.35676      0.20231      8.34051E-02
  326.   1.93513E-02  9.73832E-04  2.47997E-04  4.06084E-05  2.72083E-05
  327.    0.0000;
  328.   lr=extr evolv 'ORDO';
  329.   list lr ;
  330.   ER=SOMM( abs (lr - lrr) ) ;
  331.   mess ' Erreur sur IMPL-PG ' er ;
  332.   si ( er > 1.e-4) ; erreur 5 ; finsi ;
  333. FINSI ;
  334.  
  335. * SCHEMA BDF2
  336.  
  337. CFL=1.  ;
  338. BETA=0.5 ;
  339. OPDFDT='BDF2' ;
  340. OPKONV='SUPG' ;
  341. TIT= CHAI 'BDF2  (1/4 de TOUR)' ;
  342. rv = TEST OPDFDT OPKONV BETA  CFL DIAM GRAPH TIT;
  343. evolh=evol 'CHPO' (RV.'INCO'.'TN')  daxh ;
  344. evolv=evol 'CHPO' (RV.'INCO'.'TN')  daxv ;
  345. evolth=evolth et evolh ;
  346. evoltv=evoltv et evolv ;
  347.  
  348. SI(NON COMPLET ) ;
  349. lrr=prog
  350.    0.0000      5.49984E-08  2.57798E-07 -9.63347E-07  1.52810E-05
  351.   2.97655E-04  4.82244E-03  3.65515E-02  0.16889      0.45454    
  352.   0.75960      0.84691      0.62439      0.24678     -3.80415E-02
  353.  -0.10517     -1.80764E-03 -1.91432E-03  1.90002E-04  4.46618E-04
  354.    0.0000;
  355.   lr=extr evolv 'ORDO';
  356.   list lr ;
  357.   ER=SOMM( abs (lr - lrr) ) ;
  358.   mess ' Erreur sur BDF2 ' er ;
  359.   si ( er > 1.e-4) ; erreur 5 ; finsi ;
  360. FINSI ;
  361.  
  362. si ('EGA' graph 'O' );
  363. titre 'Coupe ox à y=1/2 ' ;
  364. TAB1=TABLE;
  365. TAB1.'TITRE'=TABLE ;
  366. TAB1.1='TIRR MARQ CROI REGU';
  367. TAB1.'TITRE' . 1 = mot 'T initial ' ;
  368. TAB1.2='MARQ ETOI REGU';
  369. TAB1.'TITRE' . 2 = mot 'CNG  CFL 0.9';
  370. TAB1.3='MARQ ETOI REGU';
  371. TAB1.'TITRE' . 3 = mot 'TVISQ CFL 0.5';
  372. TAB1.4='MARQ CROI REGU ';
  373. TAB1.'TITRE' . 4 = mot 'CN  CFL 1. ';
  374. TAB1.5='MARQ LOSA REGU ';
  375. TAB1.'TITRE' . 5 = mot 'CNPG CFL 1.' ;
  376. TAB1.6='MARQ CARR REGU ';
  377. TAB1.'TITRE' . 6 = mot 'IMPLPG CFL 1.' ;
  378. TAB1.7='MARQ PLUS REGU ';
  379. TAB1.'TITRE' . 7 = mot 'BDF2 CFL 1.' ;
  380. DESS EVOLTH 'TITX' 'Ox'    LEGE  TAB1 ;
  381. titre 'Coupe oy à x=0   ' ;
  382. DESS EVOLTV 'TITX' 'Oy'    LEGE  TAB1 ;
  383. FINSI ;
  384.  
  385.  FIN ;
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  

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