Télécharger la plaquette Cast3M


fr en

Télécharger condmixtesVF.dgibi

Retour à la liste

Numérotation des lignes : 

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

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