Télécharger la plaquette Cast3M


fr en

Télécharger transport1VF_vs_EFMH.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  transport1_new.dgibi
  2. *
  3. ********************* CAS TEST : transport1.dgibi ********************
  4. *
  5. GRAPH = FAUX    ;
  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 verifier le bon fonctionnement des operateurs
  20. * utilises afin de resoudre l'equation de transport par diffusion
  21. * convection d'un champ scalaire passif par la methode d'elements
  22. * finis mixtes hybrides et VF implantes dans CASTEM2000.
  23. * Le decentrement est utilise. 
  24. *
  25. *
  26. *---------------------------------------------------------------------
  27. *
  28. *------------------
  29. * Options generales
  30. *------------------
  31. *
  32. 'OPTI' 'DIME' 2 'ELEM' 'QUA4' ;
  33. 'OPTI' 'ISOV' 'SURF'          ;
  34. *
  35. *
  36. *=========
  37. * MAILLAGE
  38. *=========
  39. *
  40. *
  41. *- Creation des points supports du contour du domaine, et des droites
  42. *- passant par les centres et les faces pour le post-traitement.
  43. *
  44. L     = 100.D0          ;
  45. H     = 20.D0            ;
  46. X0    = 0.D0            ;
  47. X1    = X0 + L          ;
  48. Y0    = 0.D0            ;
  49. Y1    = Y0 + H          ;
  50. INUMX = 141              ;
  51. INUMY = 17               ;
  52.  
  53. DX    = X1 - X0 / INUMX ;
  54. DY    = H '/' INUMY     ;
  55.  
  56. *
  57. A1    = X0 Y0   ;
  58. A3    = X1 Y0   ;
  59. D1    = X0 Y1   ;
  60. D3    = X1 Y1   ;
  61.  
  62. *
  63. *- Creation des DROITES frontieres
  64. *
  65. DRBAS = A3 'DROI' INUMX A1 ;
  66. DRGAU = A1 'DROI' INUMY D1 ;
  67. DRHAU = D1 'DROI' INUMX D3 ;
  68. DRDRO = D3 'DROI' INUMY A3 ;
  69. PELIM = DY / (100.D0) ;
  70. *
  71. *- Creation maillage GEOMETRIQUE
  72. *
  73. PTOT1  = 'DALL' DRBAS DRGAU DRHAU DRDRO ;
  74. PTOT2  = 'ORIE' PTOT1                   ;
  75. *
  76. *- Creation maillage HYBRIDE y compris sous-objets (cond. limites)
  77. *
  78. *DRMID = B1 'DROI' INUMX B3              ;
  79. *DRMIC = C1 'DROI' INUM1 C3              ;
  80. *EXT1  = 'MANU' 'POI1' B1                ;
  81. *
  82. QFTOT= 'CHAN' PTOT2  QUAF ;
  83. QFGAU= 'CHAN' DRGAU  QUAF ;
  84. QFDRO= 'CHAN' DRDRO  QUAF ;
  85.  ELIM PELIM (QFTOT ET QFGAU ET QFDRO) ;
  86. *
  87. *================
  88. * INITIALISATIONS
  89. *================
  90. *
  91. *                                                    ----------------
  92. *                                                    = MODELISATION =
  93. *                                                    ----------------
  94. MODHYB = MODE QFTOT 'DARCY' 'ANISOTROPE' ;
  95. MODDRO = MODE QFDRO 'DARCY' 'ISOTROPE' ;
  96. MODGAU = MODE QFGAU 'DARCY' 'ISOTROPE' ;
  97. CHYB1 = 'DOMA' MODHYB 'SURFACE'          ;
  98. CHYB2 = 'DOMA' MODHYB 'NORMALE'          ;
  99.  
  100.  
  101. xcoor ycoor = 'COORDONNEE' ('DOMA' modhyb centre);
  102.  
  103.  
  104. *
  105. *                                               ---------------------
  106. *                                               = Donnees physiques =
  107. *                                               ---------------------
  108. *
  109. T0    = 0.D0 ; T1  = 100.D0    ;
  110. VX1   = 1.D0 ; VY1 = 0.D0    ;
  111. VK    = 1.D-1               ;
  112. MATI2 = MANU 'CHPO' ('DOMA' MODHYB 'CENTRE')  'K' VK ;
  113. MATI2 = ('NOMC' K11 MATI2) '+' ('NOMC' K22 mati2) '+' ('NOMC' K21 (0.D0 * MATI2));
  114. *MATI2 = KCHA MODHYB MATI2 'CHAM';
  115. *
  116. SPEED = 'MANU' 'CHPO' ('DOMA' MODHYB 'FACE')  2 'VX' VX1 'VY' VY1 'NATURE' 'DISCRET' ;
  117. SPEEDC = 'MANU' 'CHPO' ('DOMA' MODHYB 'CENTRE')  2 'VX' VX1 'VY' VY1 'NATURE' 'DISCRET' ;
  118. MOT1  = 'MOTS' 'VX' 'VY'                     ;
  119. MOT2  = 'MOTS' 'UX' 'UY'                     ;
  120. FLU1  = 'PSCA' SPEED CHYB2 MOT1 MOT2         ;
  121. FLU2  = CHYB1 * FLU1                         ;
  122. FLU3  = 'NOMC' 'FLUX' FLU2                   ;
  123. *
  124. *                                             -----------------------
  125. *                                             = Donnees transitoire =
  126. *                                             -----------------------
  127. *  TETA   : Parametre de le theta-methode
  128. *  TMAX   : Temps final
  129. *  TSUP   : Temps pour conditions aux limites
  130. *  DELTAT : Pas de temps
  131. *
  132. TETA   = 1.00D0       ;
  133. TMIN   = 0.D0         ;
  134. TMAX   = 60.00D0      ;
  135. TSUP   = 1.2D0 * TMAX ;
  136. DELTAT = 1.0D-0        ;
  137. *
  138. LICALC = 'PROG' TMIN 'PAS' DELTAT TMAX ;
  139. LISAUV = 'PROG' TMAX                   ;
  140. *
  141. *                                             ------------------------
  142. *                                             = Conditions initiales =
  143. *                                             ------------------------
  144. *TCHYB =  'MANU' 'CHPO' ('DOMA' MODHYB 'CENTRE') 1 'H' 1.D0
  145. *              'NATURE' 'DISCRET' ;
  146.  
  147. INUMX1 = INUMX '+' 1;
  148. INUMY1 = INUMY '+' 1;
  149. TCHYB = 'MASQUE' xcoor SUPERIEUR (0.2D0 * L);
  150. TCHYB = TCHYB * ( 'MASQUE' xcoor INFERIEUR (0.3D0 * L));
  151. TCHYB = TCHYB * ( 'MASQUE' ycoor INFERIEUR (0.6D0 * H));
  152. TCHYB = TCHYB * ( 'MASQUE' ycoor SUPERIEUR (0.4D0 * H));
  153. TCHYB = 'NOMC' 'H' TCHYB;
  154.  
  155. * On rend la masse initiale independante du maillage.
  156. TCHYB = TCHYB / (maxi (resu (TCHYB * (doma modhyb volume))));
  157.  
  158.  
  159.  
  160. *
  161. *                                                       --------------
  162. *                                                       = T imposee  =
  163. *                                                       --------------
  164. *opti donn 5;
  165. CHCLIM = TABLE;
  166. CHCLIM . 'DIRICHLET' = 'MANU' 'CHPO' ('DOMA' MODDRO 'CENTRE') 1 'H' T0 'NATURE' 'DISCRET' ;
  167. CHCLIM . 'DIRICHLET' = CHCLIM . 'DIRICHLET' + ('MANU' 'CHPO' ('DOMA' MODGAU 'CENTRE') 1 'H' T1 'NATURE' 'DISCRET') ;
  168.  
  169. GEOL1 = TABLE;
  170. GEOL1 . 'CONCENTRATION'     = 'MANUEL' CHPO ('DOMA' modhyb 'CENTRE') 'H' 0.D0 ;
  171. GEOL1 . 'LUMP'              = FAUX                  ;
  172. GEOL1 . 'TYPDISCRETISATION' = 'EFMH'                ;
  173. GEOL1 . 'THETA_DIFFUSION'   = 1.0D0                 ;
  174. GEOL1 . 'THETA_CONVECTION'  = 1.0D0                 ;
  175. GEOL1 . 'DECENTREMENT'      = FAUX                  ;
  176. GEOL1 . 'DELTAT'            = 10.D15                ;
  177. * conduct en m/an
  178. GEOL1 . 'DIFFUSIVITE'       =  mati2 '/' vk      ;
  179. GEOL1 . 'SOLVEUR'           = 2                     ;
  180. GEOL1 . 'PRECONDITIONNEUR'  = 3                     ;
  181. GEOL1 . 'POROSITE'          = 'MANU' 'CHPO' (doma MODHYB CENTRE) 'CK' 0.D0 ;
  182. GEOL1 . 'CLIMITES'          = CHCLIM                ;
  183. GEOL1 . 'RECALCUL'          = VRAI                  ;
  184. *
  185.  
  186. *'OPTION' donn 5;
  187.  
  188. GEOL1 GEOL2 = TRANGEOL Modhyb GEOL1;
  189.  
  190.  
  191.  
  192.  
  193.  
  194.  
  195.  
  196.  
  197.  
  198. QFACE1 =  GEOL1 . 'FLUXDIFF'               ;
  199. VCENT1 = 'HVIT' MODHYB (nomc FLUX QFACE1)              ;
  200.  
  201.  
  202. *
  203. *                                          ---------------------------
  204. *                                          = Table DARCY_TRANSITOIRE =
  205. *                                          ---------------------------
  206. *
  207. *
  208. *
  209.  
  210. *-- Table de transport :
  211. Transp                         = 'TABLE';
  212. Transp . 'MODELE'              = MODHYB ;
  213. Transp.'TEMPS'                 = 'TABLE';
  214. Transp.'CONCENTRATION'         = 'TABLE';
  215. Transp.'FLUXDIFF'              = 'TABLE';
  216. Transp.'FLUXCONV'              = 'TABLE';
  217. Transp.'CARACTERISTIQUES'      = 1.D-10 * MATI2;
  218. Transp.'POROSITE'              = 1.D0;
  219. Transp.'COEF_RETARD'           = 1.D0;
  220. *Transp.'CONVECTION'            =  NOMC 'SCAL' FLU3 ;
  221. Transp.'CONVECTION'            =  (nomc SCAL (1.D0 * QFACE1)) / (doma MODHYB SURFACE) ;
  222. Transp . 'CONVECTION'        = Transp . CONVECTION * (doma Modhyb NORMALE);
  223. Transp.'VITELEM'            =   VCENT1 ;
  224. *Transp . 'ALPHAL'           = MANU 'CHPO' (doma MODHYB centre)
  225. *                              'SCAL' 1.D0;
  226. *Transp . 'ALPHAT'           = MANU 'CHPO' (doma MODHYB centre)
  227. *                              'SCAL' 0.D0;
  228.  
  229. * Conditions initiales :
  230. Transp.'TEMPS'. 0              = TMIN ;
  231. Transp.'CONCENTRATION'. 0      = TCHYB;
  232. Transp.'FLUXDIFF'. 0           = 0.D0 * FLU3 ;
  233. Transp.'FLUXCONV'. 0           = 0.D0 * FLU3;
  234.  
  235. * Conditions aux limites :
  236. Transp . 'TRACE_IMPOSE' = CHAR ('MANU' 'CHPO' ('DOMA' MODGAU 'CENTRE') 1 'H' 0.D0  'NATURE' 'DISCRET') ('EVOL' 'MANU' ('PROG' 0. TSUP) ('PROG' 1. 1.))       ;
  237. *Transp . 'FLUXTOT_IMP' =  TTT1       ;
  238.  
  239.  
  240. * Paramètres numeriques :
  241. Transp.'THETA_DIFF'            = 1.D0;
  242. Transp.'THETA_CONVECTION'      = 1.D0;
  243. Transp.'LUMP'                  = FAUX;
  244. Transp.'TYPDISCRETISATION'     = 'EFMH';
  245. Transp.'DECENTR'          = VRAI;
  246.  
  247. TABRES = table METHINV;
  248. TABRES . 'TYPINV' = 1;
  249. TABRES . 'PRECOND' = 3;
  250.  
  251. Transp . 'METHINV' = TABRES;
  252.  
  253.  
  254. Transp.'TEMPS_CALCULES'        = LiCalc;
  255. *Transp.'TEMPS_SAUVES'          = LiSauv;
  256.  
  257.  
  258. Transp . INTCONC = TABLE;
  259. Transp . INTCONC . 0 = 0.D0 * TCHYB;
  260. *                                                            ==========
  261. *                                                            | CALCUL |
  262. *                                                            ==========
  263.  
  264.  
  265. *=======================
  266. * Resolution transitoire
  267. *=======================
  268. *
  269. TRANSGEN TRANSP ;
  270. *
  271. *
  272. *=================
  273. * POST-TRAITEMENT
  274. *=================
  275.  
  276. ****************************************************************************
  277. 'SI' (GRAPH) ;
  278. titre 'Premier champ de concentration' ;
  279. trac (kcha TRANSP . CONCENTRATION . 1 MODHYB CHAM) modhyb;
  280. *
  281. NbTps = 'DIME' (Transp.'TEMPS') ;
  282. titre 'Dernier champ de concentration calcule' ;
  283. trac (kcha (TRANSP . CONCENTRATION . (NbTps - 1)) MODHYB CHAM) modhyb;
  284. 'FINSI' ;
  285.  
  286.  
  287.  
  288. LiMASS = 'PROG' ;
  289. LiTps = 'PROG' ;
  290. zozo = 0.D0;
  291. roro = 'PROG';
  292. litpm1 = 'PROG' ;
  293. *
  294.  
  295.  
  296.  
  297. 'REPETER' BclFlT (('DIME' (Transp.'TEMPS')) - 1 );
  298.   ICour   = (&BclFlT) ;
  299.   toto = 'KOPS' transp . CONCENTRATION . Icour  * transp . POROSITE;
  300.   toto = toto * ('DOMA' modhyb VOLUME);
  301. *  'LISTE' ('MAXIMUM' ('RESULT' toto));
  302.    LiTps   = LiTps 'ET' ('PROG' (Transp. 'TEMPS' . ICour)) ;
  303.   FlTCour = 'MAXIMUM' ('RESULT' toto) ;
  304. *  'LISTE' fltcour;
  305.   LiMASS   = LiMASS 'ET' ('PROG' FlTCour) ;
  306. 'FIN' BclFlT ;
  307.  
  308.  
  309. mailimg = 'DOMA' modgau centre;
  310. mailimd = 'DOMA' moddro centre;
  311.  
  312.  
  313. EvMASS = 'EVOL' 'ROUGE' 'MANUEL' 'Time (y)' Litps 'Total mass' (LiMASS) ;
  314.  
  315. 'SI' (graph) ;
  316. 'DESSIN' EvMASS 'TITRE' 'evolution de la masse dans les fractures' ;
  317. 'FINSI' ;
  318.  
  319.  
  320. LiFLT = 'PROG' ;
  321. LiTps = 'PROG' ;
  322. zozo = 0.D0;
  323. roro = 'PROG';
  324. litpm1 = 'PROG' ;
  325. 'REPETER' BclFlT ('DIME' (Transp.'TEMPS')) ;
  326.   ICour   = (&BclFlT) '-' 1 ;
  327.   LiTps   = LiTps 'ET' ('PROG' (Transp. 'TEMPS' . ICour)) ;
  328.   LesFlCo = ('REDU' (TRANSP . FLUXCONV . ICour '+' TRANSP . FLUXDIFF . ICour ) mailimg) ;
  329.   FlTCour = 'MAXIMUM' ('RESULT' LesFlCo) ;
  330. *  'LISTE' fltcour;
  331.   LiFLT   = LiFLT 'ET' ('PROG' FlTCour) ;
  332.   SI ( Icour < (('DIME' (Transp.'TEMPS')) '-' 1));
  333.   zozo = zozo '+' (fltcour * (((Transp. 'TEMPS' . (ICour+1))) '-' ((Transp. 'TEMPS' . ICour))));
  334.   roro = roro 'ET' ('PROG' zozo);
  335.   LiTpm1 = LiTpm1 'ET' ('PROG' (Transp. 'TEMPS' . ICour)) ;
  336.   'FINSI' ;
  337. 'FIN' BclFlT ;
  338. *
  339. XConvSY = 1.D0 ;
  340. LiTpsy = LiTps '/' (XConvSY) ;
  341. EvFLT = 'EVOL' 'ROUGE' 'MANUEL' 'Time (y)' LiTpsy 'Total FLux' (LiFLT '*' XconvSY) ;
  342. *
  343. ChTitr = 'CHAINE' 'Evolution du flux (1/year) en sortie (positif sortant)' ;
  344.  
  345. 'SI' (GRAPH) ;
  346. 'DESSIN' EvFLT 'TITRE' ChTitr ;
  347. 'FINSI' ;
  348.  
  349. EvFLTS = 'EVOL' 'JAUNE' 'MANUEL' 'Time (y)' LiTpm1 'Total FLux (1/y) sortant'  roro ;
  350.  
  351. ChTitr = 'CHAINE' 'Evolution du flux total integre dans le temps en sortie (positif sortant)';
  352.  
  353.  
  354.  
  355. *
  356. LiFLT = 'PROG' ;
  357. LiTps = 'PROG' ;
  358. *
  359. *  Pas de CORRECTION normale sortante à droite
  360.  
  361. zozo = 0.D0;
  362. roro = 'PROG';
  363. litpm1 = 'PROG' ;
  364. 'REPETER' BclFlT ('DIME' (Transp.'TEMPS')) ;
  365.   ICour   = (&BclFlT) '-' 1 ;
  366.   LiTps   = LiTps 'ET' ('PROG' (Transp. 'TEMPS' . ICour)) ;
  367.   LesFlCo =  ('REDU' (TRANSP . FLUXCONV . ICour '+' TRANSP . FLUXDIFF . ICour) mailimd)  ;
  368.   FlTCour = 'MAXIMUM' ('RESULT' LesFlCo) ;
  369.   LiFLT   = LiFLT 'ET' ('PROG' FlTCour) ;
  370.   SI ( Icour < (('DIME' (Transp.'TEMPS')) '-' 1));
  371.   zozo = zozo '+' (fltcour * (((Transp. 'TEMPS' . (ICour+1))) '-' ((Transp. 'TEMPS' . ICour))));
  372.   roro = roro 'ET' ('PROG' zozo);
  373.   LiTpm1 = LiTpm1 'ET' ('PROG' (Transp. 'TEMPS' . ICour)) ;
  374.   'FINSI' ;
  375. 'FIN' BclFlT ;
  376.  
  377. LiTpsy = LiTps '/' (XConvSY) ;
  378. EvFLT = 'EVOL' 'BLEU' 'MANUEL' 'Time (y)' LiTpsy 'Total FLux (1/y)' (LiFLT*xconvsy) ;
  379.  
  380. ChTitr = 'CHAINE' 'Evolution du flux total en entree (positif sortant)';
  381.  
  382. 'SI' (graph) ;
  383. 'DESSIN' EvFLT 'TITRE' ChTitr ;
  384. 'FINSI';
  385.  
  386. EvFLTE = 'EVOL' 'JAUNE' 'MANUEL' 'Time (y)' LiTpm1 'Total FLux (1/y) entree ' roro ;
  387.  
  388. ChTitr = 'CHAINE' 'Evolution du flux total integre dans le temps à lentree 'ET' sortie (positif sortant)';
  389. 'SI' (graph) ;
  390. 'DESSIN' (EvFLTE 'ET' evflts) 'TITRE' ChTitr ;
  391.  
  392. 'DESSIN' (evflts '+' evflte) TITRE 'flux total en entree ET sortie';
  393.  
  394. 'DESSIN' (evflts '+' evflte '+' evmass) TITRE 'conservation masse totale';
  395. 'FINSI' ;
  396.  
  397. * Projection sur des SUPPORTS COMMUNS
  398. ABSCi  = EXTR evflts 'ABSC' ;
  399. LRflteA= IPOL ABSCi evflte  ;
  400. LRmassA= IPOL ABSCi evmass  ;
  401. evflteA='EVOL' 'JAUNE' 'MANUEL' 'Time (y)' ABSCi 'Total FLux (1/y) entree ' LRflteA ;
  402. evmassA='EVOL' 'ROUGE' 'MANUEL' 'Time (y)' ABSCi 'Total mass' LRmassA ;
  403.  
  404. Bilan = evflts '+' evflteA '+' evmassA '-' 1.D0 ;
  405. LRE1= 'EXTR' Bilan 'ORDO';
  406. err1=('MAXI' 'ABS' LRE1) ;
  407. 'MESS' 'conservation de la masse' err1;
  408.  
  409. 'SI' (err1 > 1.D-11 ) ;
  410.    'MESS' 'mauvaise conservation de la masse:' err1 ;
  411.    'ERREUR' 5;
  412. 'SINON' ;
  413. *  pas d'erreur de conservation
  414.    ERCONSEF = VRAI;   
  415. 'FINSI'   ;
  416.  
  417.  
  418. **************post VF ********************************************
  419. * on sauve la derniere concentration EFMH
  420.  
  421. concEFMH = Transp . CONCENTRATION . (('DIME' transp . TEMPS) '-' 1);
  422.  
  423.  
  424.  
  425. *-- Table de transport :
  426. Transp                         = 'TABLE';
  427. Transp . 'MODELE'              = MODHYB ;
  428. Transp.'TEMPS'                 = 'TABLE';
  429. Transp.'CONCENTRATION'         = 'TABLE';
  430. Transp.'FLUXDIFF'              = 'TABLE';
  431. Transp.'FLUXCONV'              = 'TABLE';
  432. Transp.'CARACTERISTIQUES'      = 1.D-10 * MATI2;
  433. Transp.'POROSITE'              = 1.D0;
  434. Transp.'COEF_RETARD'           = 1.D0;
  435. *Transp.'CONVECTION'            =  NOMC 'SCAL' FLU3 ;
  436. Transp.'CONVECTION'            =  (nomc SCAL (1.D0 * QFACE1)) / (doma MODHYB SURFACE) ;
  437. Transp . 'CONVECTION'        = Transp . CONVECTION * (doma Modhyb NORMALE);
  438. Transp.'VITELEM'            =   VCENT1 ;
  439. *Transp . 'ALPHAL'           = MANU 'CHPO' (doma MODHYB centre)
  440. *                              'SCAL' 1.D0;
  441. *Transp . 'ALPHAT'           = MANU 'CHPO' (doma MODHYB centre)
  442. *                              'SCAL' 0.D0;
  443.  
  444. * Conditions initiales :
  445. Transp.'TEMPS'. 0              = TMIN ;
  446. Transp.'CONCENTRATION'. 0      = TCHYB;
  447. Transp.'FLUXDIFF'. 0           = 0.D0 * FLU3 ;
  448. Transp.'FLUXCONV'. 0           = 0.D0 * FLU3;
  449.  
  450. * Conditions aux limites :
  451. Transp . 'TRACE_IMPOSE' = CHAR ('MANU' 'CHPO' ('DOMA' MODGAU 'CENTRE') 1 'H' 0.D0  'NATURE' 'DISCRET') ('EVOL' 'MANU' ('PROG' 0. TSUP) ('PROG' 1. 1.))       ;
  452. *Transp . 'FLUXTOT_IMP' =  TTT1       ;
  453.  
  454.  
  455. * Paramètres numeriques :
  456. Transp.'THETA_DIFF'            = 1.D0;
  457. Transp.'THETA_CONVECTION'      = 1.D0;
  458. Transp.'LUMP'                  = FAUX;
  459. Transp.'TYPDISCRETISATION'     = 'VF';
  460. Transp.'DECENTR'          = VRAI;
  461.  
  462. TABRES = table METHINV;
  463. TABRES . 'TYPINV' = 1;
  464. TABRES . 'PRECOND' = 3;
  465.  
  466. Transp . 'METHINV' = TABRES;
  467.  
  468.  
  469. Transp.'TEMPS_CALCULES'        = LiCalc;
  470. *Transp.'TEMPS_SAUVES'          = LiSauv;
  471.  
  472.  
  473. Transp . INTCONC = TABLE;
  474. Transp . INTCONC . 0 = 0.D0 * TCHYB;
  475. *                                                            ==========
  476. *                                                            | CALCUL |
  477. *                                                            ==========
  478.  
  479.  
  480. *=======================
  481. * Resolution transitoire
  482. *=======================
  483. *
  484. TRANSGEN TRANSP ;
  485. *
  486. *
  487. *=================
  488. * POST-TRAITEMENT
  489. *=================
  490.  
  491. ****************************************************************************
  492. 'SI' (GRAPH) ;
  493. titre 'Premier champ de concentration' ;
  494. trac (kcha TRANSP . CONCENTRATION . 1 MODHYB CHAM) modhyb;
  495. *
  496. NbTps = 'DIME' (Transp.'TEMPS') ;
  497. titre 'Dernier champ de concentration calcule' ;
  498. trac (kcha (TRANSP . CONCENTRATION . (NbTps - 1)) MODHYB CHAM) modhyb;
  499. 'FINSI' ;
  500.  
  501.  
  502.  
  503. LiMASS = 'PROG' ;
  504. LiTps = 'PROG' ;
  505. zozo = 0.D0;
  506. roro = 'PROG';
  507. litpm1 = 'PROG' ;
  508. *
  509.  
  510.  
  511.  
  512. 'REPETER' BclFlT (('DIME' (Transp.'TEMPS')) - 1 );
  513.   ICour   = (&BclFlT) ;
  514.   toto = 'KOPS' transp . CONCENTRATION . Icour  * transp . POROSITE;
  515.   toto = toto * ('DOMA' modhyb VOLUME);
  516. *  'LISTE' ('MAXIMUM' ('RESULT' toto));
  517.    LiTps   = LiTps 'ET' ('PROG' (Transp. 'TEMPS' . ICour)) ;
  518.   FlTCour = 'MAXIMUM' ('RESULT' toto) ;
  519. *  'LISTE' fltcour;
  520.   LiMASS   = LiMASS 'ET' ('PROG' FlTCour) ;
  521. 'FIN' BclFlT ;
  522.  
  523.  
  524. mailimg = 'DOMA' modgau centre;
  525. mailimd = 'DOMA' moddro centre;
  526.  
  527.  
  528. EvMASS = 'EVOL' 'ROUGE' 'MANUEL' 'Time (y)' Litps 'Total mass' (LiMASS) ;
  529.  
  530. 'SI' (graph) ;
  531. 'DESSIN' EvMASS 'TITRE' 'evolution de la masse dans les fractures' ;
  532. 'FINSI' ;
  533.  
  534.  
  535. LiFLT = 'PROG' ;
  536. LiTps = 'PROG' ;
  537. zozo = 0.D0;
  538. roro = 'PROG';
  539. litpm1 = 'PROG' ;
  540. 'REPETER' BclFlT ('DIME' (Transp.'TEMPS')) ;
  541.   ICour   = (&BclFlT) '-' 1 ;
  542.   LiTps   = LiTps 'ET' ('PROG' (Transp. 'TEMPS' . ICour)) ;
  543.   LesFlCo = ('REDU' (TRANSP . FLUXCONV . ICour '+' TRANSP . FLUXDIFF . ICour ) mailimg) ;
  544.   FlTCour = 'MAXIMUM' ('RESULT' LesFlCo) ;
  545. *  'LISTE' fltcour;
  546.   LiFLT   = LiFLT 'ET' ('PROG' FlTCour) ;
  547.   SI ( Icour < (('DIME' (Transp.'TEMPS')) '-' 1));
  548.   zozo = zozo '+' (fltcour * (((Transp. 'TEMPS' . (ICour+1))) '-' ((Transp. 'TEMPS' . ICour))));
  549.   roro = roro 'ET' ('PROG' zozo);
  550.   LiTpm1 = LiTpm1 'ET' ('PROG' (Transp. 'TEMPS' . ICour)) ;
  551.   'FINSI' ;
  552. 'FIN' BclFlT ;
  553. *
  554. XConvSY = 1.D0 ;
  555. LiTpsy = LiTps '/' (XConvSY) ;
  556. EvFLT = 'EVOL' 'ROUGE' 'MANUEL' 'Time (y)' LiTpsy 'Total FLux' (LiFLT '*' XconvSY) ;
  557. *
  558. ChTitr = 'CHAINE' 'Evolution du flux (1/year) en sortie (positif sortant)' ;
  559.  
  560. 'SI' (GRAPH) ;
  561. 'DESSIN' EvFLT 'TITRE' ChTitr ;
  562. 'FINSI' ;
  563.  
  564. EvFLTS = 'EVOL' 'JAUNE' 'MANUEL' 'Time (y)' LiTpm1 'Total FLux (1/y) sortant'  roro ;
  565.  
  566. ChTitr = 'CHAINE' 'Evolution du flux total integre dans le temps en sortie (positif sortant)';
  567.  
  568.  
  569.  
  570. *
  571. LiFLT = 'PROG' ;
  572. LiTps = 'PROG' ;
  573. *
  574. *  Pas de CORRECTION normale sortante à droite
  575.  
  576. zozo = 0.D0;
  577. roro = 'PROG';
  578. litpm1 = 'PROG' ;
  579. 'REPETER' BclFlT ('DIME' (Transp.'TEMPS')) ;
  580.   ICour   = (&BclFlT) '-' 1 ;
  581.   LiTps   = LiTps 'ET' ('PROG' (Transp. 'TEMPS' . ICour)) ;
  582.   LesFlCo =  ('REDU' (TRANSP . FLUXCONV . ICour '+' TRANSP . FLUXDIFF . ICour) mailimd)  ;
  583.   FlTCour = 'MAXIMUM' ('RESULT' LesFlCo) ;
  584.   LiFLT   = LiFLT 'ET' ('PROG' FlTCour) ;
  585.   SI ( Icour < (('DIME' (Transp.'TEMPS')) '-' 1));
  586.   zozo = zozo '+' (fltcour * (((Transp. 'TEMPS' . (ICour+1))) '-' ((Transp. 'TEMPS' . ICour))));
  587.   roro = roro 'ET' ('PROG' zozo);
  588.   LiTpm1 = LiTpm1 'ET' ('PROG' (Transp. 'TEMPS' . ICour)) ;
  589.   'FINSI' ;
  590. 'FIN' BclFlT ;
  591.  
  592. LiTpsy = LiTps '/' (XConvSY) ;
  593. EvFLT = 'EVOL' 'BLEU' 'MANUEL' 'Time (y)' LiTpsy 'Total FLux (1/y)' (LiFLT*xconvsy) ;
  594.  
  595. ChTitr = 'CHAINE' 'Evolution du flux total en entree (positif sortant)';
  596.  
  597. 'SI' (graph) ;
  598. 'DESSIN' EvFLT 'TITRE' ChTitr ;
  599. 'FINSI';
  600.  
  601. EvFLTE = 'EVOL' 'JAUNE' 'MANUEL' 'Time (y)' LiTpm1 'Total FLux (1/y) entree ' roro ;
  602.  
  603. ChTitr = 'CHAINE' 'Evolution du flux total integre dans le temps à lentree 'ET' sortie (positif sortant)';
  604.              'SI' (graph) ;
  605. 'DESSIN' (EvFLTE 'ET' evflts) 'TITRE' ChTitr ;
  606.  
  607. 'DESSIN' (evflts '+' evflte) TITRE 'flux total en entree ET sortie';
  608.  
  609. 'DESSIN' (evflts '+' evflte '+' evmass) TITRE 'conservation masse totale';
  610. 'FINSI' ;
  611.  
  612. * Projection sur des SUPPORTS COMMUNS
  613. ABSCi  = EXTR evflts 'ABSC' ;
  614. LRflteA= IPOL ABSCi evflte  ;
  615. LRmassA= IPOL ABSCi evmass  ;
  616. evflteA='EVOL' 'JAUNE' 'MANUEL' 'Time (y)' ABSCi 'Total FLux (1/y) entree ' LRflteA ;
  617. evmassA='EVOL' 'ROUGE' 'MANUEL' 'Time (y)' ABSCi 'Total mass' LRmassA ;
  618.  
  619. Bilan = evflts '+' evflteA '+' evmassA '-' 1.D0 ;
  620. LRE2  = 'EXTR' Bilan 'ORDO';
  621. err2  =('MAXI' 'ABS' LRE2) ;
  622. 'MESS' 'conservation de la masse' err2;
  623.  
  624.  
  625.  
  626.  
  627. 'SI' ( err2 > 1.D-11 ) ;
  628.    'MESS' 'mauvaise conservation de la masse:' err2 ;
  629.    'ERREUR' 5;
  630. 'SINON' ;
  631. *  pas d'erreur de conservation
  632.    ERCONSVF = VRAI;
  633. 'FINSI'    ;
  634.  
  635.  
  636. * On sauve la concentration VF finale
  637.  
  638. concVF =  Transp . CONCENTRATION . (('DIME' transp . TEMPS) '-' 1);
  639.  
  640.  
  641. **************** COMPAR VF - EFMH **********************************
  642.  
  643. toto = concEFMH '-' concVF ;
  644. toto = 'KOPS' toto * toto;
  645. toto = ('DOMA' modhyb volume) * toto;
  646. toto = 'MAXIMUM' (('RESULT' toto));
  647. titi = ('DOMA' modhyb volume) * ('KOPS' concVF * concVF);
  648. titi = 'MAXIMUM' ('RESULT' titi);
  649. err3 = (toto '/' titi)**0.5D0;
  650.  
  651. 'MESSAGE' 'erreur relative L2' err3;
  652.  
  653. 'SI' (err3 > 1.D-2) ;
  654.    'MESS' 'Comparaison VF - EFMH mauvaise:' err3 ;
  655.    'ERREUR' 5 ;
  656. 'SINON' ;
  657.     erEFVF = VRAI;
  658. 'FINSI' ;
  659.  
  660.  
  661. **************** POST HYDRAU ****************************************
  662.  
  663. FLTotE  = 'RESULT' (('REDU' QFACE1 mailimd)) ;
  664. FLTotE = -1.D0 * FltotE;
  665. fltote = maxi fltote;
  666.  
  667. * FLux sortants
  668. FLTotS  = 'RESULT' (('REDU' QFACE1 mailimg)) ;
  669. fltots = maxi fltots;
  670.  
  671. * Affichage des resultats
  672.  
  673. 'MESSAGE' 'Resultats hydraulique : ' ;
  674. 'MESSAGE' ' ' ;
  675. 'MESSAGE' 'FLux entrant (total, rapporte à la surface) ' FLTotE  ;
  676. 'MESSAGE' 'FLux sortant (total, rapporte à la surface) ' FLTotS  ;
  677. 'MESSAGE' ' ' ;
  678. 'MESSAGE' 'Keq associes ' Keq1 Keq2 ;
  679.  
  680.  
  681. 'SI' (erefvf 'ET' erconsvf 'ET' erconsef) ;
  682.   'ERREUR' 0;
  683. 'SINON' ;
  684.   'ERREUR' 5;
  685. 'FINSI' ;
  686.  
  687.  
  688. 'FIN';
  689.  
  690.  
  691.  
  692.  
  693.  
  694.  
  695.  
  696.  
  697.  

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