Télécharger la plaquette Cast3M


fr en

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