Télécharger la plaquette Cast3M


fr en

Télécharger transport1EFMH.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  transport1EFMH.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; sur les autres frontières ce sont
  34. * les conditions aux limites naturelles.
  35. *
  36. * Le schéma en temps utilisé est le schéma d'Euler implicite
  37. * (teta-méthode avec teta=1.0).
  38. *
  39. * Données physiques :
  40. *    VK    : Coefficient de diffusion .................    1.D-4
  41. *    VX1   : Vitesse suivant x.........................    1.
  42. *    VY1   : Vitesse suivant y.........................    0.
  43. *    T0    : Concentration initiale....................    0.
  44. *    T1    : Concentration en entrée...................    1.
  45. *
  46. * Données temporelles :
  47. *    TMIN  : Temps initial.............................    0.
  48. *    TMAX  : Temps final...............................   30.
  49. *    DELTAT: Pas de temps..............................    1.
  50. *    T0    : Concentration initiale....................    0.
  51. *    T1    : Concentration en entrée...................    1.
  52. *
  53. * Parametres de maillage :
  54. *    L     : longueur du domaine .......................  100.
  55. *    H     : hauteur du domaine ........................    1.
  56. *    INUMX : Nombre d'éléments en x ....................  100
  57. *    INUMY : Nombre d'éléments en y ....................    1 car 1D
  58. *
  59. *---------------------------------------------------------------------
  60. *
  61. *------------------
  62. * Options generales
  63. *------------------
  64. *
  65. 'OPTI' 'DIME' 2 'ELEM' 'QUA4' ;
  66. 'OPTI' 'ISOV' 'SURF'          ;
  67. *
  68. *
  69. *=========
  70. * MAILLAGE
  71. *=========
  72. *
  73. *
  74. *- Création des points supports du contour du domaine, et des droites
  75. *- passant par les centres et les faces pour le post-traitement.
  76. *
  77. L     = 100.D0          ;
  78. LS2   = L / 2.D0        ;
  79. H     = 20.D0            ;
  80. HS2   = H / 2.D0        ;
  81. X0    = 0.D0            ;
  82. X1    = X0 + L          ;
  83. Y0    = 0.D0            ;
  84. Y1    = Y0 + H          ;
  85. INUMX = 100             ;
  86. INUMY = 1               ;
  87. *INUM1 = INUMX - 1       ;
  88. *Y01   = Y0 + Y1 * 0.5D0 ;
  89. DX    = X1 - X0 / INUMX ;
  90. DY    = H '/' INUMY     ;
  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 = DY / (100.D0) ;
  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) ;
  128. *
  129. *================
  130. * INITIALISATIONS
  131. *================
  132. *
  133. *                                                    ----------------
  134. *                                                    = MODELISATION =
  135. *                                                    ----------------
  136. MODHYB = MODE QFTOT 'DARCY' 'ANISOTROPE' ;
  137. MODDRO = MODE QFDRO 'DARCY' 'ISOTROPE' ;
  138. MODGAU = MODE QFGAU 'DARCY' 'ISOTROPE' ;
  139. CHYB1 = 'DOMA' MODHYB 'SURFACE'          ;
  140. CHYB2 = 'DOMA' MODHYB 'NORMALE'          ;
  141.  
  142.  
  143. xcoor ycoor = 'COORDONNEE' ('DOMA' modhyb centre);
  144.  
  145.  
  146. *
  147. *                                               ---------------------
  148. *                                               = Donnees physiques =
  149. *                                               ---------------------
  150. *
  151. T0    = 0.D0 ; T1  = 1.D0    ;
  152. VX1   = 1.D0 ; VY1 = 0.D0    ;
  153. VK    = 1.D-4               ;
  154. MATI2 = MANU 'CHPO' ('DOMA' MODHYB 'CENTRE')  'K' VK ;
  155. *MATI2 = KCHA MODHYB MATI2 'CHAM';
  156. *
  157. SPEED = 'MANU' 'CHPO' ('DOMA' MODHYB 'FACE')  2
  158.         'VX' VX1 'VY' VY1 'NATURE' 'DISCRET' ;
  159. SPEEDC = 'MANU' 'CHPO' ('DOMA' MODHYB 'CENTRE')  2
  160.         'VX' VX1 'VY' VY1 'NATURE' 'DISCRET' ;
  161. MOT1  = 'MOTS' 'VX' 'VY'                     ;
  162. MOT2  = 'MOTS' 'UX' 'UY'                     ;
  163. FLU1  = 'PSCA' SPEED CHYB2 MOT1 MOT2         ;
  164. FLU2  = CHYB1 * FLU1                         ;
  165. FLU3  = 'NOMC' 'FLUX' FLU2                   ;
  166. *
  167. *                                             -----------------------
  168. *                                             = Donnees transitoire =
  169. *                                             -----------------------
  170. *  TETA   : Parametre de le theta-méthode
  171. *  TMAX   : Temps final
  172. *  TSUP   : Temps pour conditions aux limites
  173. *  DELTAT : Pas de temps
  174. *
  175. TETA   = 1.00D0       ;
  176. TMIN   = 0.D0         ;
  177. TMAX   = 30.00D0      ;
  178. TSUP   = 1.2D0 * TMAX ;
  179. DELTAT = 0.5D0        ;
  180. *
  181. LICALC = 'PROG' TMIN 'PAS' DELTAT TMAX ;
  182. LISAUV = 'PROG' TMAX                   ;
  183. *
  184. *                                             ------------------------
  185. *                                             = Conditions initiales =
  186. *                                             ------------------------
  187. *TCHYB =  'MANU' 'CHPO' ('DOMA' MODHYB 'CENTRE') 1 'H' 1.D0
  188. *                                             'NATURE' 'DISCRET' ;
  189. TCHYB = 'MASQUE' xcoor SUPERIEUR ((L/2.D0) '-' (dx/2.D0));
  190. TCHYB = TCHYB * (
  191.     'MASQUE' xcoor INFERIEUR ((L/2.D0) '+' (dx/2.D0)));
  192. TCHYB = TCHYB * (
  193.     'MASQUE' ycoor INFERIEUR ((H/2.D0) '+' (dy/2.D0)));
  194. TCHYB = TCHYB * (
  195.     'MASQUE' xcoor SUPERIEUR ((L/2.D0) '-' (dx/2.D0)));
  196. TCHYB = 'NOMC' 'H' TCHYB;
  197. *
  198. *                                                       --------------
  199. *                                                       = T imposée  =
  200. *                                                       --------------
  201. *opti donn 5;
  202.  
  203. *
  204. *                                          ---------------------------
  205. *                                          = Table DARCY_TRANSITOIRE =
  206. *                                          ---------------------------
  207. *
  208. *
  209. *
  210.  
  211. *-- Table de transport :
  212. Transp                         = 'TABLE';
  213. Transp . 'MODELE'              = MODHYB ;
  214. Transp.'TEMPS'                 = 'TABLE';
  215. Transp.'CONCENTRATION'         = 'TABLE';
  216. Transp.'FLUXDIFF'              = 'TABLE';
  217. Transp.'FLUXCONV'              = 'TABLE';
  218. Transp.'CARACTERISTIQUES'      = MATI2;
  219. Transp.'POROSITE'              = 1.D0;
  220. Transp.'COEF_RETARD'           = 1.D0;
  221. Transp.'CONVECTION'            =  NOMC 'SCAL' FLU3 ;
  222. Transp.'CONVECTION'            =  Transp.'CONVECTION' /
  223.                            (doma MODHYB SURFACE) ;
  224. Transp . 'CONVECTION'        = Transp . CONVECTION *
  225.                            (doma Modhyb NORMALE);
  226. Transp.'VITELEM'            =   SPEEDC ;
  227. *Transp . 'ALPHAL'           = MANU 'CHPO' (doma MODHYB centre)
  228. *                              'SCAL' 1.D0;
  229. *Transp . 'ALPHAT'           = MANU 'CHPO' (doma MODHYB centre)
  230. *                              'SCAL' 0.D0;
  231.  
  232. * Conditions initiales :
  233. Transp.'TEMPS'. 0              = TMIN ;
  234. Transp.'CONCENTRATION'. 0      = TCHYB;
  235. Transp.'FLUXDIFF'. 0           = 0.D0 * FLU3 ;
  236. Transp.'FLUXCONV'. 0           = 0.D0 * FLU3;
  237.  
  238. * Conditions aux limites :
  239. Transp . 'TRACE_IMPOSE' = CHAR ('MANU' 'CHPO' ('DOMA' MODGAU 'CENTRE') 1
  240.               'H' 1.D0  'NATURE' 'DISCRET')
  241.               ('EVOL' 'MANU' ('PROG' 0. TSUP) ('PROG' 1. 1.))       ;
  242. *Transp . 'FLUXTOT_IMP' =  TTT1       ;
  243.  
  244.  
  245. * Paramètres numériques :
  246. Transp.'THETA_DIFF'            = 1.D0;
  247. Transp.'THETA_CONVECTION'      = 1.D0;
  248. Transp.'LUMP'                  = FAUX;
  249. Transp.'TYPDISCRETISATION'     = 'EFMH';
  250. Transp.'DECENTR'          = VRAI;
  251.  
  252. TABRES = table METHINV;
  253. TABRES . 'TYPINV' = 1;
  254. TABRES . 'PRECOND' = 3;
  255.  
  256. Transp . 'METHINV' = TABRES;
  257.  
  258.  
  259. Transp.'TEMPS_CALCULES'        = LiCalc;
  260. Transp.'TEMPS_SAUVES'          = LiSauv;
  261.  
  262.  
  263. Transp . INTCONC = TABLE;
  264. Transp . INTCONC . 0 = 0.D0 * TCHYB;
  265. *                                                            ==========
  266. *                                                            | CALCUL |
  267. *                                                            ==========
  268.  
  269.  
  270. *=======================
  271. * Resolution transitoire
  272. *=======================
  273. *
  274. TRANSGEN TRANSP ;
  275.  
  276.  
  277.  
  278.  
  279.  
  280. CFL = VX1 * DELTAT / DX                   ;
  281. PEK = VX1 * DX / (2. * VK)                ;
  282. FOU = 2. * VK * DELTAT / (DX * DX)        ;
  283. *
  284. XC YC = 'COOR' (DOMA MODHYB 'CENTRE')   ;
  285. ISOR1 = INDEX ( TRANSp . 'TEMPS') ;
  286. NTSOR = DIME ISOR1                ;
  287. NTSO1 = NTSOR - 1                 ;
  288. IOK   = FAUX                      ;
  289. IRESU = 1                         ;
  290. *
  291. *-----------------------
  292. REPETER VISURESU NTSO1 ;
  293. *-----------------------
  294. *
  295. IRESU = IRESU + 1                        ;
  296. INDI1 = ISOR1.IRESU                      ;
  297. TTRA  = TRANSP . 'TEMPS' . INDI1         ;
  298. *                                              Principe du maximum
  299. EPR1  = TRANSP . 'CONCENTRATION' . INDI1        ;
  300. PMA   = 'MAXI' EPR1 ; PMI = 'MINI' EPR1  ;
  301. VTEST = PMA + PMI   ;
  302. 'SAUT' 1 'LIGNE'                         ;
  303. 'MESS' '    Max et min théorique  ' T1 T0     ;
  304. 'MESS' '    Max et min par maille ' PMA PMI   ;
  305. 'MESS' '    Précision demandée    ' CRIT1     ;
  306. *                                              Centre de gravité
  307. M0 = 'KCHA' MODHYB 'CHAM' EPR1 ;
  308. M1 = 'INTG' MODHYB M0         ;
  309. *
  310. XV1   = XC * EPR1               ;
  311. MXV1  = 'KCHA' MODHYB 'CHAM' XV1 ;
  312. XV2   = 'INTG' MODHYB MXV1      ;
  313. XV3   = XV2 / M1                ;
  314. YV1   = YC * EPR1               ;
  315. MYV1  = 'KCHA' MODHYB 'CHAM' YV1 ;
  316. YV2   = 'INTG' MODHYB MYV1      ;
  317. YV3   = YV2 / M1                ;
  318. XTEST = TTRA / 2.D0             ;
  319. *
  320. 'SAUT' 1 'LIGNE'                            ;
  321. 'MESS' '    THEORIQUE           ' XTEST HS2     ;
  322. 'MESS' '    OBTENU              ' XV3 YV3       ;
  323. 'MESS' '    TAILLE D UNE MAILLE ' DX            ;
  324. 'SAUT' 1 'LIGNE'                            ;
  325. 'SI' ( VTEST 'NEG' T1 CRIT1 )           ;
  326.      IOK = VRAI                         ;
  327.      'MESSAGE' 'vtest';
  328. 'FINS'                                  ;
  329. 'SI' ( XV3 'NEG' XTEST DX    )  ;
  330.      IOK = VRAI                 ;
  331.      'MESSAGE' 'XV3';
  332. 'FINS'                          ;
  333. 'SI' ( YV3 'NEG' HS2  CRIT1 )   ;
  334.      IOK = VRAI                 ;
  335.      'MESSAGE' 'YV3';
  336. 'FINS'                          ;
  337. *
  338. 'SI' ('NEG' GRAPH 'N' )                        ;
  339.    LTI2 = 'CHAINE' 'Front 1D-h temps ' TTRA    ;
  340.    'TITR' LTI2                                 ;
  341.    iepr1 = 'ELNO' modhyb epr1;
  342.    AV2  = 'EVOL' 'ROUG' 'CHPO' iEPR1 'SCAL'  DRhau ;
  343.    'DESS' AV2 'MIMA'                           ;
  344. 'FINS'                                         ;
  345. *
  346. *-------------
  347. FIN VISURESU ;
  348. *-------------
  349. *
  350. *
  351. 'SI' ( IOK )  ;
  352.      'ERRE' 5 ;
  353. 'SINO'        ;
  354.      'ERRE' 0 ;
  355. 'FINSI'       ;
  356. *
  357. 'FIN' ;
  358.  
  359.  
  360. *
  361. *
  362. *=================
  363. * POST-TRAITEMENT
  364. *=================
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  

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