Télécharger fluxtotalVF.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : fluxtotalVF.dgibi
  2. *
  3. ********************* CAS TEST : transport1.dgibi ********************
  4. *
  5. GRAPH = 'N' ;
  6. 'OPTI' 'ECHO' 1 ;
  7. CRIT1 = 1.D-6 ;
  8. 'SAUT' 'PAGE' ;
  9. *
  10. *---------------------------------------------------------------------
  11. * TEST TRANSPORT1
  12. * CALCUL DARCY ISOTROPE TRANSPORT.
  13. * Transport d'un front.
  14. *
  15. * dT
  16. * -- + div (uT - Kgrad(T)) = 0
  17. * dt
  18. *
  19. * Ce test permet de vérifier le bon fonctionnement des opérateurs
  20. * utilisés afin de résoudre l'équation de transport par diffusion
  21. * convection d'un champ scalaire passif par la méthode d'éléments
  22. * finis mixtes hybrides implanté dans CASTEM2000.
  23. *
  24. * Le maillage ne comporte qu'un élément dans la direction transverse
  25. * à l'écoulement puisque le phénomène étudié est 1D. Dans le sens de
  26. * l'écoulement, le pas de discrétisation est constant.
  27. *
  28. * Les conditions initiales correspondent à l'entrée du front dans le
  29. * domaine. Comme le champ de vitesse choisi est positif, la trace de
  30. * concentration est égale à 1 à gauche du domaine.
  31. *
  32. * Les conditions aux limites imposent la trace de concentration à
  33. * gauche du domaine égale à un via une condition de flux total
  34. * (soit U C = U pour avoit C = 1) ; sur les autres frontières ce sont
  35. * les conditions aux limites naturelles.
  36. *
  37. * Le schéma en temps utilisé est le schéma d'Euler implicite
  38. * (teta-méthode avec teta=1.0).
  39. *
  40. * Données physiques :
  41. * VK : Coefficient de diffusion ................. 1.D-2
  42. * VX1 : Vitesse suivant x......................... 1.
  43. * VY1 : Vitesse suivant y......................... 0.
  44. * T0 : Concentration initiale.................... 0.
  45. * T1 : Concentration en entrée................... 1.
  46. *
  47. * Données temporelles :
  48. * TMIN : Temps initial............................. 0.
  49. * TMAX : Temps final............................... 30.
  50. * DELTAT: Pas de temps.............................. 1.
  51. * T0 : Concentration initiale.................... 0.
  52. * T1 : Concentration en entrée................... 1.
  53. *
  54. * Parametres de maillage :
  55. * L : longueur du domaine ....................... 100.
  56. * H : hauteur du domaine ........................ 1.
  57. * INUMX : Nombre d'éléments en x .................... 100
  58. * INUMY : Nombre d'éléments en y .................... 1 car 1D
  59. *
  60. *---------------------------------------------------------------------
  61. *
  62. *------------------
  63. * Options generales
  64. *------------------
  65. *
  66. 'OPTI' 'DIME' 2 'ELEM' 'QUA4' ;
  67. 'OPTI' 'ISOV' 'SURF' ;
  68. *
  69. *
  70. *=========
  71. * MAILLAGE
  72. *=========
  73. *
  74. *
  75. *- Création des points supports du contour du domaine, et des droites
  76. *- passant par les centres et les faces pour le post-traitement.
  77. *
  78. L = 100.D0 ;
  79. LS2 = L / 2.D0 ;
  80. H = 10.D0 ;
  81. HS2 = H / 2.D0 ;
  82. X0 = 0.D0 ;
  83. X1 = X0 + L ;
  84. Y0 = 0.D0 ;
  85. Y1 = Y0 + H ;
  86. INUMX = 100 ;
  87. INUMY = 1 ;
  88. INUM1 = INUMX - 1 ;
  89. Y01 = Y0 + Y1 * 0.5D0 ;
  90. DX = X1 - X0 / INUMX ;
  91. DX1 = DX / 2.D0 ;
  92. XG = X0 + DX1 ;
  93. XD = X1 - DX1 ;
  94. *
  95. A1 = X0 Y0 ;
  96. A3 = X1 Y0 ;
  97. D1 = X0 Y1 ;
  98. D3 = X1 Y1 ;
  99. B1 = X0 Y01 ;
  100. B3 = X1 Y01 ;
  101. C1 = XG Y01 ;
  102. C3 = XD Y01 ;
  103. P6 = LS2 Y01 ;
  104. *
  105. *- Création des DROITES frontieres
  106. *
  107. DRBAS = A3 'DROI' INUMX A1 ;
  108. DRGAU = A1 'DROI' INUMY D1 ;
  109. DRHAU = D1 'DROI' INUMX D3 ;
  110. DRDRO = D3 'DROI' INUMY A3 ;
  111. PELIM = DX1 / (5. * INUMX) ;
  112. *
  113. *- Creation maillage GEOMETRIQUE
  114. *
  115. PTOT1 = 'DALL' DRBAS DRGAU DRHAU DRDRO ;
  116. PTOT2 = 'ORIE' PTOT1 ;
  117. *
  118. *- Creation maillage HYBRIDE y compris sous-objets (cond. limites)
  119. *
  120. DRMID = B1 'DROI' INUMX B3 ;
  121. DRMIC = C1 'DROI' INUM1 C3 ;
  122. EXT1  = 'MANU' 'POI1' B1 ;
  123. *
  124. QFTOT= 'CHAN' PTOT2 QUAF ;
  125. QFGAU= 'CHAN' DRGAU QUAF ;
  126. QFDRO= /CHAN' DRDRO QUAF ;
  127. ELIM PELIM (QFTOT ET QFGAU ET QFDRO ET DRMID ET DRMIC ET EXT1) ;
  128. *
  129. *================
  130. * INITIALISATIONS
  131. *================
  132. *
  133. * ----------------
  134. * = MODELISATION =
  135. * ----------------
  136. MODHYB = MODE QFTOT 'DARCY' 'ANISOTROPE' ;
  137. MODDRO = MODE QFDRO 'DARCY' 'ANISOTROPE' ;
  138. MODGAU = MODE QFGAU 'DARCY' 'ANISOTROPE' ;
  139. CHYB1 = 'DOMA' MODHYB 'SURFACE' ;
  140. CHYB2 = 'DOMA' MODHYB 'NORMALE' ;
  141.  
  142.  
  143.  
  144. *
  145. * ---------------------
  146. * = Donnees physiques =
  147. * ---------------------
  148. *
  149. T0 = 0.D0 ; CTT1 = 1.D0 ;
  150.  
  151. T1 = (H ) / INUMY;
  152.  
  153. VX1 = 1.D0 ; VY1 = 0.D0 ;
  154. VK = 1.D-2 ;
  155. MATI2 = 'MATE' MODHYB DIRECTION (1. 0.)
  156. 'K11' VK 'K21' 0.D0 'K22' VK;
  157. MATI2 = KCHA MODHYB MATI2 'CHPO';
  158. *
  159. SPEED = 'MANU' 'CHPO' ('DOMA' MODHYB 'FACE') 2
  160. 'VX' VX1 'VY' VY1 'NATURE' 'DISCRET' ;
  161. *
  162. * -----------------------
  163. * = Donnees transitoire =
  164. * -----------------------
  165. * TETA : Parametre de le theta-méthode
  166. * TMAX : Temps final
  167. * TSUP : Temps pour conditions aux limites
  168. * DELTAT : Pas de temps
  169. *
  170. TETA = 0.00D0 ;
  171. TMIN = 0.D0 ;
  172. TMAX2 = 15.00D0 ;
  173. TMAX = 30.00D0 ;
  174. TSUP = 1.2D0 * TMAX ;
  175. DELTAT = 0.5D0 ;
  176. *
  177. LICALC = 'PROG' TMIN 'PAS' (DELTAT/2.D0) TMAX2 ;
  178. LICALC = LICALC ET ('PROG' (TMAX2 + DELTAT) 'PAS' DELTAT TMAX) ;
  179. LISAUV = 'PROG' TMAX ;
  180. *
  181. * ------------------------
  182. * = Conditions initiales =
  183. * ------------------------
  184. TFHY0 = 'MANU' 'CHPO' ('DOMA' MODHYB 'FACE') 1 'TH' T0
  185. 'NATURE' 'DISCRET' ;
  186. TFHY1 = 'MANU' 'CHPO' EXT1 1 'TH' T1 'NATURE' 'DISCRET' ;
  187. TFHYB = TFHY0 'ET' TFHY1 ;
  188. TCHYB = 'MANU' 'CHPO' ('DOMA' MODHYB 'CENTRE') 1 'H' T0
  189. 'NATURE' 'DISCRET' ;
  190. *
  191. * --------------
  192. * = T imposée =
  193. * --------------
  194. BBGAU = 'BLOQ' ('DOMA' MODGAU 'CENTRE') 'TH' ;
  195. EEGAU = 'DEPI' BBGAU T1 ;
  196. TT1 = 'KCHT' MODGAU SCAL CENTRE (-1.* T1);
  197. TTT1 = 'CHAR' TT1 ('EVOL' 'MANU' ('PROG' 0. TSUP) ('PROG' 1. 1.)) ;
  198. CHAGAU = 'CHAR' EEGAU ('EVOL' 'MANU' ('PROG' 0. TSUP) ('PROG' 1. 1.)) ;
  199. *
  200. * ---------------------------
  201. * = Table DARCY_TRANSITOIRE =
  202. * ---------------------------
  203. *
  204. *
  205. *
  206.  
  207. *-- Table de transport :
  208. Transp = 'TABLE';
  209. Transp . 'MODELE' = MODHYB ;
  210. Transp.'TEMPS' = 'TABLE';
  211. Transp.'CONCENTRATION' = 'TABLE';
  212. Transp.'FLUXDIFF' = 'TABLE';
  213. Transp.'FLUXCONV' = 'TABLE';
  214. Transp.'CARACTERISTIQUES' = MATI2;
  215. *Transp.'POROSITE' = MaPor;
  216. Transp.'CONVECTION' = SPEED ;
  217. Transp.'VITELEM' = kcht MODHYB VECT CENTRE (1.D0 0.);
  218. Transp.'VITELEM' = nomc (mots UX UY) (mots VX VY)
  219. Transp . 'VITELEM';
  220.  
  221. * Conditions initiales :
  222. Transp.'TEMPS'. 0 = TMIN ;
  223. Transp.'CONCENTRATION'. 0 = TCHYB;
  224. Transp.'FLUXDIFF'. 0 = 0.D0 * TFHY0;
  225. Transp.'FLUXCONV'. 0 = 0.D0 * TFHY0;
  226.  
  227. * Conditions aux limites :
  228. *Transp . 'TRACE_IMPOSE' = TTT1 ;
  229. Transp . 'FLUXTOT_IMP' = TTT1 ;
  230.  
  231.  
  232. * Paramètres numériques :
  233. Transp.'THETA_DIFF' = 1.D0;
  234. Transp.'THETA_CONVECTION' = 1.0D0;
  235. Transp.'LUMP' = FAUX;
  236. Transp.'TYPDISCRETISATION' = 'VF';
  237. Transp.'DECENTR' = FAUX;
  238.  
  239. TABRES = table METHINV;
  240. TABRES . 'TYPINV' = 1;
  241. TABRES . 'PRECOND' = 3;
  242.  
  243. Transp . 'METHINV' = TABRES;
  244.  
  245.  
  246. Transp.'TEMPS_CALCULES' = LiCalc;
  247. Transp.'TEMPS_SAUVES' = LiSauv;
  248.  
  249.  
  250. Transp . INTCONC = TABLE;
  251. Transp . INTCONC . 0 = 0.D0 * TCHYB;
  252. * ==========
  253. *Transp . 'DECROISSANCE' = 0.01D0;
  254. * | CALCUL |
  255. * ==========
  256.  
  257.  
  258. *=======================
  259. * Resolution transitoire
  260. *=======================
  261. *
  262. TRANSGEN TRANSP ;
  263. *
  264. *
  265. *=================
  266. * POST-TRAITEMENT
  267. *=================
  268.  
  269. TRANS2 = TABLE TRANSP;
  270.  
  271. *
  272. *--------------------------------------------------------------------
  273. * Dans chaque cas on trace
  274. * La trace de concentration le long de DRMID
  275. * La concentration le long de DRMIC
  276. *--------------------------------------------------------------------
  277. * Tests de NON-REGRESSION :
  278. * Principe du maximum
  279. * Position du centre de gravité du champ de concentration
  280. *--------------------------------------------------------------------
  281. *
  282. * Critères numériques
  283. CFL = VX1 * DELTAT / DX ;
  284. PEK = VX1 * DX / (2. * VK) ;
  285. FOU = 2. * VK * DELTAT / (DX * DX) ;
  286. 'SAUT' 1 'LIGNE' ;
  287. 'MESS' ' Critères numériques ' ;
  288. 'MESS' ' CFL ' CFL ;
  289. 'MESS' ' PECLET ' PEK ;
  290. 'MESS' ' FOURIER ' FOU ;
  291. 'MESS' ' ' ;
  292. *
  293. XC YC = 'COOR' (DOMA MODHYB 'CENTRE') ;
  294. ISOR1 = INDEX ( TRANS2 . 'TEMPS') ;
  295. NTSOR = DIME ISOR1 ;
  296. NTSO1 = NTSOR - 1 ;
  297. IOK = FAUX ;
  298. IRESU = 1 ;
  299. *
  300. *-----------------------
  301. REPETER VISURESU NTSO1 ;
  302. *-----------------------
  303. *
  304. IRESU = IRESU + 1 ;
  305. INDI1 = ISOR1.IRESU ;
  306. TTRA = TRANS2 . 'TEMPS' . INDI1 ;
  307. * Principe du maximum
  308. EPR1 = TRANS2 . 'CONCENTRATION' . INDI1 ;
  309. PMA = 'MAXI' EPR1 ; PMI = 'MINI' EPR1 ;
  310. VTEST = PMA + PMI / 1.D0 ;
  311. 'SAUT' 1 'LIGNE' ;
  312. 'MESS' ' Principe du maximum ' ;
  313. 'MESS' ' Max et min théorique ' T1 T0 ;
  314. 'MESS' ' Max et min par maille ' PMA PMI ;
  315. 'MESS' ' Précision demandée ' CRIT1 ;
  316. 'MESS' ' ' ;
  317. * Centre de gravité
  318. M0 = 'KCHA' MODHYB 'CHAM' EPR1 ;
  319. M1 = 'INTG' MODHYB M0 ;
  320. *
  321. XV1 = XC * EPR1 ;
  322. MXV1 = 'KCHA' MODHYB 'CHAM' XV1 ;
  323. XV2 = 'INTG' MODHYB MXV1 ;
  324. XV3 = XV2 / M1 ;
  325. YV1 = YC * EPR1 ;
  326. MYV1 = 'KCHA' MODHYB 'CHAM' YV1 ;
  327. YV2 = 'INTG' MODHYB MYV1 ;
  328. YV3 = YV2 / M1 ;
  329. XTEST = TTRA / 2.D0 ;
  330. *
  331. 'SAUT' 1 'LIGNE' ;
  332. 'MESS' ' Centre de gravité du nuage ' ;
  333. 'MESS' ' THEORIQUE ' XTEST HS2 ;
  334. 'MESS' ' OBTENU ' XV3 YV3 ;
  335. 'MESS' ' TAILLE D UNE MAILLE ' DX ;
  336. 'MESS' ' ' ;
  337. 'SAUT' 1 'LIGNE' ;
  338. * Tests de non regression
  339. 'SI' ( VTEST 'NEG' CTT1 CRIT1 ) ;
  340. IOK = VRAI ;
  341. 'FINS' ;
  342. 'SI' ( XV3 'NEG' XTEST DX ) ;
  343. IOK = VRAI ;
  344. 'FINS' ;
  345. 'SI' ( YV3 'NEG' HS2 CRIT1 ) ;
  346. IOK = VRAI ;
  347. 'FINS' ;
  348. *
  349. 'SI' ('NEG' GRAPH 'N' ) ;
  350. LTI2 = 'CHAINE' 'Front 1D-h temps ' TTRA ;
  351. 'TITR' LTI2 ;
  352. AV2 = 'EVOL' 'ROUG' 'CHPO' EPR1 'H' DRMIC ;
  353. 'DESS' AV2 'MIMA' ;
  354. 'FINS' ;
  355. *
  356. *-------------
  357. FIN VISURESU ;
  358. *-------------
  359. *
  360. *
  361. 'SI' ( IOK ) ;
  362. mess 'erreur' ((VTEST - T1) / T1);
  363. mess 'erreur' ((XV3 - XTEST) / XTEST);
  364. mess 'erreur' ((YV3 - HS2) / HS2);
  365. 'ERRE' 5 ;
  366. 'SINO' ;
  367. 'ERRE' 0 ;
  368. 'FINSI' ;
  369. *
  370. 'FIN' ;
  371.  
  372.  
  373.  
  374.  
  375.  

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