Télécharger conem.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : conem.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='MACRO' ;
  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
  153. 'MATCONS' 'MTCV'
  154. 'ZONE' HYTOT 'OPER' 'KONV' 1. 'UN' COEF DT 'INCO' 'TN'
  155. 'CLIM' 'TN' 'TIMP' DRCOT 0. ;
  156.  
  157. Si ( EGA OPDFDT 'BDF2');
  158. RV = 'EQEX' RV
  159. 'OPTI' 'EF' 'CENTREE' BDF2
  160. 'ZONE' HYTOT 'OPER' 'DFDT' 1. 'TNM' 'TNMM' DT 'UN' coef 'INCO' 'TN'
  161. ;
  162. Sinon ;
  163. RV = 'EQEX' RV
  164. 'OPTI' 'EF' OPDFDT
  165. 'ZONE' HYTOT 'OPER' 'DFDT' 1. 'TNM' DT 'UN' coef 'INCO' 'TN'
  166. ;
  167. FINSI ;
  168.  
  169. RV.'INCO' = TABLE 'INCO' ;
  170. RV.'INCO'.'TN' = KCHT HYTOT SCAL SOMMET 0. ;
  171. RV.'INCO'.'TNM' = KCHT HYTOT SCAL SOMMET CRS4 ;
  172. RV.'INCO'.'TNMM' = KCHT HYTOT SCAL SOMMET CRS4 ;
  173. RV.'INCO'.'T0' = KCHT HYTOT SCAL SOMMET 0. ;
  174. RV.'INCO'.'UN' = KCHT HYTOT VECT SOMMET SPEEX SPEEY ;
  175.  
  176. SI ('EGA' graph 'O' );
  177. u = vect RV.'INCO'.'UN' 1. ux uy ;
  178. trac u ptot1 'NCLK';
  179. trac RV.'INCO'.'TNM' ptot1 'NCLK';
  180. Finsi ;
  181.  
  182. rv.'OMEGA'= 1. ;
  183. rv.'NITER'=1 ;
  184.  
  185. exec rv ;
  186.  
  187. maxc=maxi (RV.'INCO'.'TN') ; minc=mini (RV.'INCO'.'TN');
  188. mess 'MAXC=' maxc 'MINC=' minc ;
  189. TN=(RV.'INCO'.'TN') ;
  190. s2c= SOMT ((NOEL HYTOT TN)*XVOL) ;
  191. mess ' S2C= ' s2c ;
  192.  
  193. SI ('EGA' graph 'O' );
  194. P1 = 'PROG' -5. -5. 2. ;
  195. MONTAGNE (RV.'INCO'.'TN') ptot1 2. P1 'TITRE' tit 'CACHE' ;
  196. 'TRACER' RV.'INCO'.'TN' ptot1 ('PROG' -0.1 pas 0.1 0.9) drcot
  197. 'TITRE' TIT 'NCLK';
  198. FINSI ;
  199. FINPROC RV ;
  200.  
  201. * SCHEMA CRANK-NICHOLSON GENERALISE GALERKIN (CNG)
  202.  
  203. CFL=0.3 ;
  204. BETA=1. ;
  205. OPDFDT='CENTREE' ;
  206. OPKONV='CNG' ;
  207. TIT= CHAI 'CNG (1/4 de TOUR)' ;
  208. rv = TEST OPDFDT OPKONV BETA CFL DIAM GRAPH TIT;
  209. evolh=evol 'CHPO' (RV.'INCO'.'TN') daxh ;
  210. evolv=evol 'CHPO' (RV.'INCO'.'TN') daxv ;
  211. evolth=evolth et evolh ;
  212. evoltv=evoltv et evolv ;
  213.  
  214. SI(NON COMPLET ) ;
  215. lrr=prog
  216. 0.0000 2.18498E-04 -1.30556E-03 -1.17804E-04 5.72424E-04
  217. 1.98609E-04 -3.67838E-03 2.49501E-02 0.35472 0.66388
  218. 0.95009 0.66041 0.32277 1.85347E-02 1.33351E-02
  219. -1.43291E-02 5.74531E-03 -4.85268E-03 -5.00940E-03 3.95789E-03
  220. 0.0000;
  221. lr=extr evolv 'ORDO';
  222. list lr ;
  223. ER=SOMM( abs (lr - lrr) ) ;
  224. mess ' Erreur sur CNG ' er ;
  225. si ( er > 1.e-4 ) ; erreur 5 ; finsi ;
  226. FINSI ;
  227.  
  228. * SCHEMA Tenseur Visqueux (TVISQ)
  229.  
  230. CFL=0.2 ;
  231. BETA=0. ;
  232. OPDFDT='CENTREE' ;
  233. OPKONV='TVISQUEU' ;
  234. TIT= CHAI 'TVISQUEU (1/4 de TOUR)' ;
  235. rv = TEST OPDFDT OPKONV BETA CFL DIAM GRAPH TIT;
  236. evolh=evol 'CHPO' (RV.'INCO'.'TN') daxh ;
  237. evolv=evol 'CHPO' (RV.'INCO'.'TN') daxv ;
  238. evolth=evolth et evolh ;
  239. evoltv=evoltv et evolv ;
  240.  
  241. SI(NON COMPLET ) ;
  242. lrr=prog
  243. 0.0000 8.25958E-08 -2.63735E-06 2.29954E-05 -1.63805E-04
  244. 1.10280E-03 -6.85438E-03 2.29617E-02 0.31028 0.66016
  245. 0.90422 0.70314 0.34157 4.19636E-02 -3.59161E-03
  246. -5.01307E-06 7.59100E-04 -4.31731E-04 -9.34623E-05 2.18262E-05
  247. 0.0000;
  248. lr=extr evolv 'ORDO';
  249. list lr ;
  250. ER=SOMM( abs (lr - lrr) ) ;
  251. mess ' Erreur sur Tenseur Visqueux ' er ;
  252. si ( er > 1.e-4 ) ; erreur 5 ; finsi ;
  253. FINSI ;
  254.  
  255. * SCHEMA CRANK-NICHOLSON GALERKIN (CN)
  256.  
  257. CFL=1. ;
  258. BETA=0.5 ;
  259. OPDFDT='CENTREE' ;
  260. OPKONV='CENTREE' ;
  261. TIT= CHAI 'CN (1/4 de TOUR)' ;
  262. rv = TEST OPDFDT OPKONV BETA CFL DIAM GRAPH TIT;
  263. evolh=evol 'CHPO' (RV.'INCO'.'TN') daxh ;
  264. evolv=evol 'CHPO' (RV.'INCO'.'TN') daxv ;
  265. evolth=evolth et evolh ;
  266. evoltv=evoltv et evolv ;
  267.  
  268. SI ( NON COMPLET ) ;
  269. lrr=prog
  270. 0.0000 2.84874E-04 3.76178E-04 -4.74619E-04 -1.23049E-03
  271. 8.65360E-05 8.77272E-03 5.63920E-02 0.24097 0.60724
  272. 0.84842 0.79075 0.42250 7.93819E-02 -4.90737E-02
  273. 2.36598E-04 3.05821E-03 -2.56679E-03 -1.60211E-03 2.56096E-03
  274. 0.0000;
  275. lr=extr evolv 'ORDO';
  276. list lr ;
  277. ER=SOMM( abs (lr - lrr) ) ;
  278. mess ' Erreur sur CN ' er ;
  279. si ( er > 1.e-4 ) ; erreur 5 ; finsi ;
  280. FINSI ;
  281.  
  282. * SCHEMA CRANK-NICHOLSON PETROV-GALERKIN (CN-PG)
  283.  
  284. CFL=1. ;
  285. BETA=0.5 ;
  286. OPDFDT='SUPG' ;
  287. OPKONV='SUPG' ;
  288. TIT= CHAI 'CN-PG (1/4 de TOUR)' ;
  289. rv = TEST OPDFDT OPKONV BETA CFL DIAM GRAPH TIT;
  290. evolh=evol 'CHPO' (RV.'INCO'.'TN') daxh ;
  291. evolv=evol 'CHPO' (RV.'INCO'.'TN') daxv ;
  292. evolth=evolth et evolh ;
  293. evoltv=evoltv et evolv ;
  294.  
  295. SI (NON COMPLET ) ;
  296. lrr=prog
  297. 0.0000 -4.15500E-08 3.71752E-07 -1.51473E-06 -2.69314E-05
  298. 5.77943E-04 7.95387E-03 5.65109E-02 0.24563 0.60380
  299. 0.84205 0.77323 0.40991 8.19064E-02 -4.45954E-02
  300. -4.11115E-03 3.71337E-03 -2.85247E-04 -5.58008E-04 1.36503E-04
  301. 0.0000;
  302. lr=extr evolv 'ORDO';
  303. list lr ;
  304. ER=SOMM( abs (lr - lrr) ) ;
  305. mess ' Erreur sur CN-PG ' er ;
  306. si ( er > 1.e-4 ) ; erreur 5 ; finsi ;
  307. FINSI ;
  308.  
  309. * SCHEMA EULER IMPLICITE PETROV-GALERKIN (IMPL-PG)
  310.  
  311. CFL=1. ;
  312. BETA=1. ;
  313. OPDFDT='SUPG' ;
  314. OPKONV='SUPG' ;
  315. TIT= CHAI 'IMPL-PG (1/4 de TOUR)' ;
  316. rv = TEST OPDFDT OPKONV BETA CFL DIAM GRAPH TIT;
  317. evolh=evol 'CHPO' (RV.'INCO'.'TN') daxh ;
  318. evolv=evol 'CHPO' (RV.'INCO'.'TN') daxv ;
  319. evolth=evolth et evolh ;
  320. evoltv=evoltv et evolv ;
  321.  
  322. SI ( NON COMPLET ) ;
  323. lrr=prog
  324. 0.0000 2.72088E-06 5.67138E-05 -5.55864E-04 2.90601E-03
  325. 2.81410E-02 8.90610E-02 0.19109 0.32444 0.44825
  326. 0.51077 0.47795 0.35589 0.20087 8.24119E-02
  327. 1.96227E-02 5.88033E-04 -4.15394E-05 -1.40971E-05 1.93301E-05
  328. 0.0000;
  329. lr=extr evolv 'ORDO';
  330. list lr ;
  331. ER=SOMM( abs (lr - lrr) ) ;
  332. mess ' Erreur sur IMPL-PG ' er ;
  333. si ( er > 1.e-4 ) ; erreur 5 ; finsi ;
  334. FINSI ;
  335.  
  336. * SCHEMA BDF2
  337.  
  338. CFL=1. ;
  339. BETA=0.5 ;
  340. OPDFDT='BDF2' ;
  341. OPKONV='SUPG' ;
  342. TIT= CHAI 'BDF2 (1/4 de TOUR)' ;
  343. rv = TEST OPDFDT OPKONV BETA CFL DIAM GRAPH TIT;
  344. evolh=evol 'CHPO' (RV.'INCO'.'TN') daxh ;
  345. evolv=evol 'CHPO' (RV.'INCO'.'TN') daxv ;
  346. evolth=evolth et evolh ;
  347. evoltv=evoltv et evolv ;
  348.  
  349. SI(NON COMPLET ) ;
  350. lrr=prog
  351. 0.0000 1.54805E-10 -3.66235E-09 5.89492E-08 -8.20960E-07
  352. 3.73495E-06 3.06438E-04 6.61244E-03 6.84728E-02 0.34155
  353. 0.84425 1.1003 0.84443 0.27789 -0.14289
  354. -0.24111 -1.27424E-02 -6.61840E-03 5.55885E-03 1.59619E-03
  355. 0.0000;
  356. lr=extr evolv 'ORDO';
  357. list lr ;
  358. ER=SOMM( abs (lr - lrr) ) ;
  359. mess ' Erreur sur BDF2 ' er ;
  360. si ( er > 1.e-4 ) ; erreur 5 ; finsi ;
  361. FINSI ;
  362.  
  363. si ('EGA' graph 'O' );
  364. titre 'Coupe ox à y=1/2 ' ;
  365. TAB1=TABLE;
  366. TAB1.'TITRE'=TABLE ;
  367. TAB1.1='TIRR MARQ CROI REGU';
  368. TAB1.'TITRE' . 1 = mot 'T initial ' ;
  369. TAB1.2='MARQ ETOI REGU';
  370. TAB1.'TITRE' . 2 = mot 'CNG CFL 0.9';
  371. TAB1.3='MARQ ETOI REGU';
  372. TAB1.'TITRE' . 3 = mot 'TVISQ CFL 0.5';
  373. TAB1.4='MARQ CROI REGU ';
  374. TAB1.'TITRE' . 4 = mot 'CN CFL 1. ';
  375. TAB1.5='MARQ LOSA REGU ';
  376. TAB1.'TITRE' . 5 = mot 'CNPG CFL 1.' ;
  377. TAB1.6='MARQ CARR REGU ';
  378. TAB1.'TITRE' . 6 = mot 'IMPLPG CFL 1.' ;
  379. TAB1.7='MARQ PLUS REGU ';
  380. TAB1.'TITRE' . 7 = mot 'BDF2 CFL 1.' ;
  381. DESS EVOLTH 'TITX' 'Ox' LEGE TAB1 ;
  382. titre 'Coupe oy à x=0 ' ;
  383. DESS EVOLTV 'TITX' 'Oy' LEGE TAB1 ;
  384. FINSI ;
  385.  
  386. FIN ;
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  

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