Télécharger la plaquette Cast3M


fr en

Télécharger transport3.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  transport3.dgibi
  2. ******************** CAS TEST : transport3.dgibi ************************
  3. *
  4. GRAPH = VRAI ;
  5. CRIT1 = 1.D-3 ;
  6.  
  7. ************************************************************************
  8. *
  9. * Petit programme testant le transport en milieu poreux avec
  10. *    - décroissance radioactive implicite
  11. *    - limite de solubilité (équilibre instantané) variant brutalement en aval
  12. *      (méthode : pénalisation)
  13. *     
  14. * pour un cas test 2D simple :
  15. * Une petite source baignant dans un champ de vitesse uniforme.
  16. *
  17. *           P4  *---------------* P3
  18. *               |               |
  19. *               |               |
  20. *               |               |
  21. *           P5  |---------------| P2
  22. *               |               |
  23. *               |               |
  24. *               |               |
  25. *               |               |
  26. *           P6  *----* P9 - - - * P11
  27. *               |    |          |
  28. *           P7  *----* P8 - - - * P10
  29. *               |               |
  30. *               |               |
  31. *           P0  *---------------* P1
  32. *
  33. ************************************************************************
  34. 'OPTION' 'ECHO' 0 ;
  35. 'SAUTER' 'PAGE';
  36. 'TITRE' 'Validation du transport';
  37. *                                                         ==============
  38. *                                                         | PARAMETRES |
  39. *                                                         ==============
  40. * 
  41. * Paramètres internes :
  42. * =====================
  43. 'OPTION' 'ECHO' 0 'DIME' 2 'ELEM' 'QUA4' ;
  44. 'OPTION' 'ISOV' 'SURFACE' ;
  45. 'OPTION' 'TRAC' 'PSC' ;
  46.  
  47. * Paramètres de convivialité :
  48. * ============================
  49. * Trace les limites de solubilité ?
  50. TrLimSol= Vrai ;
  51. * Trace les profils de concentration suivant certaines lignes ?
  52. TrLEvol = Vrai ;
  53. * Trace les évolutions de concentration en certaines points ?
  54. TrPEvol = Vrai ;
  55. * Trace les champs de concentration en soluté, précipité ?
  56. TrConc  = vrai ;
  57. TrPrec  = vrai ;
  58. * Trace la quantité de précipité dissout ?
  59. TrLib   = vrai ;
  60. * Trace le flux à travers une ligne ?
  61. TrFlux  = Vrai ;
  62.  
  63. SI ('EGA' GRAPH 'N');
  64.   TrLEvol = faux ;
  65.   TrPEvol = faux ;
  66.   TrConc  = faux ;
  67.   TrPrec  = faux ;
  68.   TrLib   = faux ;
  69.   TrFlux  = faux ;
  70.   TrLimSol= faux ;
  71. FINSI;
  72.  
  73. * Grandeurs physiques :
  74. * =====================
  75. *                                                      VITESSE DE DARCY
  76. UXImp = 0.;
  77. UYImp = 1.;
  78.  
  79. *                                                              POROSITE
  80. ValPor = 1.;
  81.  
  82. *-- Facteur de forme :
  83. FacForme = 1.;
  84.  
  85. *                                                          DUREE DE VIE
  86. * Durée de demi-vie de l'élément (années), convertie ensuite en terme
  87. * de décroissance tel que N = N0 * exp (- Lambda * t) :
  88. TVie = 10.0;
  89. Lambda = ('LOG' 2) / TVie;
  90. *Lambda = 0.;
  91. Decro  = (Lambda 'NEG' 0.);
  92.  
  93. *                                                  LIMITE DE SOLUBILITE
  94. *-- Limite de solubilité en mol/l (global et fond):
  95. *   Cette valeur ne sert que si la dissolution est instantanée ou
  96. *   d'ordre 1.
  97. CSat1 = 1. ;
  98. CSat2 = 0.1 ; 
  99.  
  100. *                                                           DIFFUSIVITE
  101. *-- Diffusivité moléculaire du traceur en m²/s
  102. DifMolec = 2.d0;
  103.  
  104. *                                                    NATURE DISSOLUTION
  105. *-- Nature de la dissolution (si existe) et logiques associés :
  106. *NaturDis = 'Pas de dissolution';
  107. *NaturDis = 'Imposée';
  108. *NaturDis = 'Ordre 1';
  109. NaturDis = 'Instant_Pénal';
  110. *NaturDis = 'Instant_PFDis';
  111.  
  112. Solub = 'NEG' NaturDis 'Pas de dissolution';
  113. SoluA = 'EGA' NaturDis 'Imposée';
  114. SoluP = 'EGA' NaturDis 'Ordre 1';
  115. SoluI = ('EGA' NaturDis 'Instant_PFDis')
  116.         OU ('EGA' NaturDis 'Instant_Pénal');
  117. SoluL = SoluP OU SoluI;
  118.  
  119. *                                               PARAM. NUM. DISSOLUTION
  120. *-- Grandeur(s) associée(s) à la dissolution : 
  121. Methode = 'Pénalisation';
  122. * Coefficient de pénalisation
  123. Penalty = 1.D7;
  124. * Nombre maxi d'itérations
  125. NbItMax = 50;
  126.  
  127. *                                                   QUANTITES INITIALES
  128. *-- Concentration en marqueur dans la source à t=0 (mol/l d'eau):
  129. CActif0 = CSat1;
  130.  
  131. *-- Précipité de marqueur dans la source à t=0 (mol/dm² de milieu):
  132. SI Solub;
  133.   PActif0 = 10. * CSat1;
  134. FINSI;
  135.  
  136. *                                                       TEMPS DE CALCUL
  137. * initial, final, et pas de temps :
  138. TIni   = 0.;
  139. TFin   = 10.;
  140. DelTaT = 1. ;
  141.  
  142. *                                                   METHODES NUMERIQUES
  143. *-- Coef. d'implicitation (1 = implicite, 0 = explicite)
  144. * Diffusion :
  145. TetaDif = 1. ;
  146. * Convection :
  147. TetaCnv = 1. ;
  148. * Décroissance :
  149. TetaDec = 1. ;
  150.  
  151. *-- Coefficient d'élimination :
  152. EpsElim = .1 ;
  153.  
  154. *-- Précision de la machine :
  155. EpsPrec = 1.D-14;
  156.  
  157. *                                                          ============
  158. *                                                          | MAILLAGE |
  159. *                                                          ============
  160. * Dimensions (m) :
  161. * ----------------
  162. LX1 = 1.;
  163. LXM = 5.;
  164. LY1 = 4.;
  165. LY2 = 5.;
  166. LY3 = 8.;
  167. LYM = 12.;
  168.  
  169. * Découpages :
  170. * ------------
  171. N   = 1;
  172. NX  = 'ENTIER' (N * LXM);
  173. NY  = 'ENTIER' (N * LYM + 1);
  174. NX1 = 'MAXIMUM' ('LECT' ('ENTIER' (NX * LX1 / LXM)) 1);
  175. NX2 = NX - NX1;
  176. NY1 = 'MAXIMUM' ('LECT' ('ENTIER' (NY * LY1 / LYM)) 1);
  177. NY2 = 'MAXIMUM' ('LECT' ('ENTIER' (NY * (LY2 - LY1) / LYM)) 1);
  178. NY3 = 'MAXIMUM' ('LECT' ('ENTIER' (NY * (LY3 - LY2) / LYM)) 1);
  179. NY4 = 'MAXIMUM' ('LECT' ('ENTIER' (NY * (LYM - LY3) / LYM)) 1);
  180. NY  = NY1 + NY2 + NY3 + NY4;
  181.  
  182. * Points :
  183. * --------
  184. p0  = 0.        0.;
  185. p1  = LXM       0.;
  186. p2  = LXM       LY3;
  187. p3  = LXM       LYM;
  188. p4  = 0.        LYM;
  189. p5  = 0.        LY3;
  190. p6  = 0.        LY2;
  191. p7  = 0.        LY1;    
  192. p8  = LX1       LY1;
  193. p9  = LX1       LY2;
  194. p10 = LXM       LY1;
  195. p11 = LXM       LY2;
  196.  
  197. * Droites :
  198. * ---------
  199. 'OPTION' 'COULEUR' 'BLAN';
  200. D4D5  = 'DROITE' P4     NY4     P5;
  201. D5D6  = 'DROITE' P5     NY3     P6;
  202. D6D7  = 'DROITE' P6     NY2     P7;
  203. D7D0  = 'DROITE' P7     NY1     P0;
  204. D2D5  = 'DROITE' P2     NX      P5;
  205. 'OPTION' 'COULEUR' 'TURQ';
  206. D0D1  = 'DROITE' P0     NX      P1;
  207. D1D10 = 'DROITE' P1     NY1     P10;
  208. D10D11= 'DROITE' P10    NY2     P11;
  209. D11D2 = 'DROITE' P11    NY3     P2;
  210. D2D3  = 'DROITE' P2     NY4     P3;
  211. D3D4  = 'DROITE' P3     NX      P4;
  212. D7D8  = 'DROITE' P7     NX1     P8;
  213. D8D9  = 'DROITE' P8     NY2     P9;
  214. D9D6  = 'DROITE' P9     NX1     P6;
  215. D10D8 = 'DROITE' P10    NX2     P8;
  216. D11D9 = 'DROITE' P11    NX2     P9;
  217. D2D5  = 'DROITE' P2     NX      P5;
  218. D8D7  = 'INVERSE' D7D8;
  219. D8D10 = 'INVERSE' D10D8;
  220. D9D8  = 'INVERSE' D8D9;
  221. D6D9  = 'INVERSE' D9D6;
  222. D9D11 = 'INVERSE' D11D9;
  223. D5D2  = 'INVERSE' D2D5;
  224. D4D3  = 'INVERSE' D3D4;
  225.  
  226. *-- Pour les conditions aux limites :
  227. LInf   = D0D1;
  228. Lsup   = D3D4;
  229. LDroit = D1D10 ET D10D11 ET D11D2 ET D2D3;
  230. LGau   = 'INVERSE' (D4D5 ET D5D6 ET D6D7 ET D7D0);
  231.  
  232. *-- Pour faire la somme des flux :
  233. * pour le tranfert d'activité (loin des conditions aux limites) :
  234. LFlux  = D5D2;
  235.  
  236. * Surfaces maillées :
  237. * -------------------
  238. SActif = 'ORIENTER' ('COULEUR' ('DALLER' D7D8 D8D9 D9D6 D6D7) 'ROUG') ;
  239. SFond  = 'ORIENTER' ('COULEUR' ('DALLER' D2D3 D3D4 D4D5 D5D2) 'VERT') ;
  240. SCoeur = 'ORIENTER' ('COULEUR' (
  241.                 ('DALLER' D0D1 D1D10 (D10D8 ET D8D7) D7D0)
  242.              ET ('DALLER' D8D10 D10D11 D11D9 D9D8)
  243.              ET ('DALLER' (D6D9 ET D9D11) D11D2 D2D5 D5D6)
  244.                       ) 'TURQ' ) ;
  245. SGran  = 'COULEUR' (SCoeur et SFond) 'TURQ' ;
  246. Mtot   = SActif ET SGran;
  247.  
  248. LTrame = 'COULEUR' ((CONT SActif) ET (CONT SGran)) 'BLAN';
  249.  
  250. *
  251. * Elements contenant tous les points de calcul
  252.  
  253. QMtot  = 'CHAN' Mtot   'QUAF' ;
  254. QFond  = 'CHAN' SFond  'QUAF' ;
  255. QActif = 'CHAN' SActif 'QUAF' ;
  256. QLFlux = 'CHAN' LFlux  'QUAF' ;
  257. QInf   = 'CHAN' LInf   'QUAF' ;
  258. QSup   = 'CHAN' Lsup   'QUAF' ;
  259. ELIM 0.001 ( QMtot ET QActif  ET QFond ET QLFlux ET QInf ET QSup) ;
  260.  
  261. * Domaines :
  262. * ----------
  263. * Modèle (Mo) :
  264. MoHydro = 'MODELE' QMtot DARCY ISOTROPE ;
  265. MoFond  = 'MODELE' QFond  DARCY ISOTROPE ;
  266. MoActif = 'MODELE' QActif DARCY ISOTROPE ;
  267. MoFlux  = 'MODELE' QLFlux DARCY ISOTROPE ;
  268. MoInf   = 'MODELE' QInf  DARCY ISOTROPE ;
  269. MoSup   = 'MODELE' QSup  DARCY ISOTROPE ;
  270. CeFlux =  DOMA MoFlux 'CENTRE' ;
  271. CeHydro = 'DOMA' MoHydro 'CENTRE' ;
  272.  
  273. * Longueurs, Normales et Orientations des arêtes du domaine ;
  274. Chyb1 = 'DOMA' MoHydro 'SURFACE';
  275. Chyb2 = 'DOMA' MoHydro 'NORMALE';
  276. VolChPro = 'DOMA' MoHydro 'VOLUME';
  277.  
  278. * Création des normales à la ligne sur laquelle on mesure les flux :
  279. NleFlux  = 'MANU' 'CHPO' CeFlux 2 'UX'  0. 'UY' +1. NATURE DIFFUS;
  280. DFluxTIT = 'CHAINE' 'la ligne y=' 'FORMAT' '(F5.0)'
  281.                          ('COOR' 2 P2) 'm';
  282.  
  283. * Points où l'on suivra la concentration :
  284. * ----------------------------------------
  285. PEvol   = 'TABLE';
  286. * Centre du combustible :
  287. PEvol.1 = 'TABLE';
  288. PEvol.1 .'T' = 'CHAINE' 'dans la source' ;
  289. PEvol.1 .'P' = 'CENTRE';
  290. PEvol.1 .'M' = ( (.5 * LX1) ((LY2 + LY1) / 2.) );
  291. PEvol.1 .'M' = CeHydro 'POINT' 'PROC' (PEvol.1 .'M');
  292. * Juste après la transition aval :
  293. PEvol.2 = 'TABLE';
  294. PEvol.2 .'T' = 'CHAINE' 'post-transition aval' ;
  295. PEvol.2 .'P' = 'CENTRE';
  296. PEvol.2 .'M' = CeHydro 'POINT'
  297.                'PROC' (P5 'PLUS' (0.05 0.05));
  298.  
  299. * Lignes où capturer les évolutions des concentrations :
  300. * ------------------------------------------------------
  301. LEvol = 'TABLE';
  302.  
  303. 'REPETER' Bcl8 ('NBEL' LGau);
  304.   pt = CeHydro 'POINT' 'PROC'
  305.        ((LGau 'POINT' &Bcl8) 'PLUS' (0 0.01));
  306.   SI (&Bcl8 > 2);
  307.     LGauC = LGauC ET ('MANU' 'SEG2' pt0 pt);
  308.   FINSI;
  309.   SI (&Bcl8 'EGA' 2);
  310.     LGauC = ('MANU' 'SEG2' pt0 pt);
  311.   FINSI;
  312.   pt0 = pt;
  313. 'FIN' Bcl8;
  314. 'OUBLIER' pt0;
  315. 'OUBLIER' pt;
  316.  
  317. LEvol.1    = 'TABLE';
  318. LEvol.1 .M = 'COULEUR' LGauC 'BLAN';
  319. LEvol.1 .P = 'CENTRE';
  320. LEvol.1 .L = 'Y (m)';
  321. LEvol.1 .T = 'CHAINE' 'longitudinale' ;
  322.  
  323. *                                                        ===============
  324. *                                                        | HYDRAULIQUE |
  325. *                                                        ===============
  326. * La donnée est U. Pour obtenir VDarcy et QF, on simule une résolution
  327. * de Darcy. On prend pour cela K = 1 et grad h = U :
  328.  
  329. * Matériau associé (Ma) :
  330. MaHydro = 'MATERIAU' MoHydro 'K' 1.;
  331.  
  332. * Matrice hybride élémentaire :
  333. MeHydro = 'MHYB' MoHydro  MaHydro;
  334.  
  335. * Trace de Concentration :
  336. FaHydro = 'DOMA' MoHydro 'FACE' ;
  337. XX YY   = 'COOR' FaHydro  ;
  338. Un      = 'MANU' 'CHPO' FaHydro 1 'SCAL' 1.;
  339. THHydro = (-1.) * (UXImp * XX) - (UYImp * YY) + Un;
  340. THHydro = 'NOMC' THHydro 'TH';
  341.  
  342. * Concentration :
  343. XX YY  = 'COOR' CeHydro ;
  344. Un     = 'MANU' 'CHPO' CeHydro 1 'SCAL' 1.;
  345. HHydro = (-1.) * (UXImp * XX) - (UYImp * YY) + Un;
  346. HHydro = 'NOMC' HHydro 'H';
  347.  
  348. * Flux et vitesses :
  349. QFHydro = 'HDEB' MoHydro  MeHydro HHydro THHydro ;
  350. VCDarcy = 'HVIT' MoHydro  QFHydro ;
  351.  
  352. *                                                          =============
  353. *                                                          | TRANSPORT |
  354. *                                                          =============
  355. * Modèle propre au transport :
  356. * ----------------------------
  357. MoTransp = 'MODELE' QMtot DARCY ISOTROPE ;
  358.  
  359. * Matrice Porosité cinématique (Type ChamElem) :
  360. * ----------------------------------------------
  361. CeTransp = 'DOMA' MoTransp 'CENTRE' ;
  362. MaPor1  = 'MANU' 'CHPO' CeTransp 1 'CK' ValPor NATURE DIFFUS;
  363. MaPor   = 'KCHA' MoTransp MaPor1 'CHAM';
  364.  
  365. * Porosité * volume de maille 
  366. * ---------------------------
  367. MaPorVol= VolChPro * ('EXCO' 'CK' MaPor1 'SCAL');
  368.  
  369. * Temps :
  370. * -------
  371. * Temps à calculer
  372. LiTCalc  = 'PROG' TIni 'PAS' DeltaT Tfin;
  373. NTCalc   = 'DIME' LiTCalc;
  374. 'REMPLACER' LiTCalc NTCalc Tfin;
  375.  
  376. * Ecriture de ces temps sous forme de table ;
  377. TabTCalc = 'TABLE';
  378. 'REPETER'  Bcl26 NTCalc;
  379.   Id  = &Bcl26 - 1;
  380.   Tps = 'EXTR' LiTCalc &Bcl26;
  381.   TabTCalc.Id = Tps;
  382. 'FIN' Bcl26;
  383.  
  384. * Temps sauvegardés
  385. LiTSauv = LiTCalc ;
  386. NTSauv  = NTCalc ;
  387.  
  388. * Temps de tracage des champs
  389. LiIndGph = 'LECT' 0  5  10 ;
  390. NTGrph   = 'DIME' LiIndGph;
  391.  
  392. * Temps de tracage des évolutions
  393. LiTEvol  = litsauv ;
  394. liindEvl = 'LECT' 0 pas 1 (ntsauv - 1) ;
  395. NTEvol   = 'DIME' litevol;
  396.  
  397. * Temps servant à la création des chargements
  398. TSup     = 1.2 * TFin;
  399. DelTDef  = 0.001 * TSup;
  400.  
  401. *                                  ====================================
  402. *                                  | CONSTRUCTION DES DONNEES DE BASE |
  403. *                                  ====================================
  404. * Dissolution / Libération ;
  405. * --------------------------
  406. *-- Matrice Limite de solubilité ;
  407. * Cas de deux zones :
  408.   MaSol2 = 'MANU' 'CHPO' CeHydro 1 'H' CSat1 NATURE DIFFUS;
  409.   MaSol3 = 'MANU' 'CHPO' (DOMA MoFond 'CENTRE')  1 'H'
  410.            (CSat2 - CSat1) NATURE DIFFUS;
  411.   MaSol1 = MaSol2 + MaSol3;
  412.   MaSol1 = 'NOMC'
  413.           ('KCHT' MoHydro 'SCAL' 'CENTRE' ('NOMC' MaSol1 'SCAL')) 'H';
  414.   MaSolS = 'ELNO' MoHydro MaSol1;
  415.  
  416. * Matériau associé au transport :
  417. * -------------------------------
  418. * Terme diffusif pur (diffusivité effective) = Porosité * Diffusivité :
  419. DifEff   = DifMolec * FacForme * MaPor;
  420. MaTransp = 'MATERIAU' MoTransp 'K' DifEff ;
  421.  
  422. * Matrice hybride élémentaire :
  423. * -----------------------------
  424. MeTransp = 'MHYB' MoTransp  MaTransp;
  425.  
  426. * Conditions Initiales :
  427. * ----------------------
  428. *-- Traces de concentration (quelconque car le transport est implicité).
  429. TC0 = 'MANU' 'CHPO' FaHydro 1 'TH' 0. 'NATURE' 'DISCRET';
  430.  
  431. *-- Concentration aux centres des éléments (sert à définir la quantité
  432. *   initiale d'élément) :
  433.  Conc01 = 'MANU' 'CHPO' CeHydro 1 'H' 0.      'NATURE' 'DISCRET';
  434.  Conc02 = 'MANU' 'CHPO'( 'DOMA' MoActif 'CENTRE') 1 'H' CActif0
  435.                                              'NATURE' 'DISCRET';
  436.  Conc0  = Conc01 ET Conc02;
  437.  Conc0  = 'NOMC'
  438.           ('KCHT' MoHydro 'SCAL' 'CENTRE' ('NOMC' Conc0 'SCAL')) 'H';
  439.  
  440. *-- Flux diffusif initial (aux centres des faces) - décoratif puisque la
  441. *   diffusion est implicite :
  442. QFDif0 = 'MANU' 'CHPO' FaHydro 1 'FLUX' 0. 'NATURE' 'DISCRET';
  443.  
  444. *-- Précipité :
  445. SI Solub;
  446.  Prec01  = 'MANU' 'CHPO' ( 'DOMA' MoActif 'CENTRE') 1 'H'
  447.                            PActif0 'NATURE' 'DISCRET';
  448.  Prec02  = 'MANU' 'CHPO' CeHydro 1 'H' 0.      'NATURE' 'DISCRET';
  449.  Prec0   = Prec01 ET Prec02;
  450.  Prec0   = 'NOMC'
  451.            ('KCHT' MoHydro 'SCAL' 'CENTRE' ('NOMC' Prec0 'SCAL')) 'H';
  452. FINSI;
  453.  
  454. *-- Libération de soluté à partir du précipité (décoratif aussi):
  455. SI Solub;
  456.   Libe0  = 'MANU' 'CHPO' CeHydro 1 'H' 0. 'NATURE' 'DISCRET';
  457. FINSI;
  458.  
  459. * Conditions aux limites (type CHARGEMENT) :
  460. * ------------------------------------------
  461. * Concentration imposée en amont et en aval, flux nul ailleurs.
  462. Ceinf= DOMA MoInf 'CENTRE' ;
  463. Cesup= DOMA Mosup 'CENTRE' ;
  464.  CCL    = 'MANU' 'CHPO' (Ceinf et Cesup) 1 'TH' 0.
  465.                'NATURE' 'DISCRET';
  466.  RigCL  = 'BLOQUE' (Ceinf et Cesup) 'TH';
  467.  CharCL = 'DEPIMPOSE' RigCL CCL ;
  468.  ChgCL  = 'CHARG' CharCL
  469.                   ('EVOL' 'MANU' ('PROG' 0. TSup) ('PROG' 1. 1.));
  470.  
  471. *                                       ===============================
  472. *                                       | Vérification des paramètres |
  473. *                                       |       de CONVERGENCE        |
  474. *                                       ===============================
  475. MESS 'Vérification des paramètres de CONVERGENCE :';
  476. MESS '--------------------------------------------';
  477.  
  478. * Critères sur la décroissance :
  479. * ------------------------------
  480. * La partie explicite de la décroissance peut rendre C < 0
  481. * Prévenir si        (1 - béta) * Lambda * DeltaT > 1
  482. * bien que ce ne soit pas un critère suffisant pour avoir C < 0.
  483. * De plus le calcul devient imprécis si DeltaT < TVie
  484. SI Decro;
  485.   Broc = (1. - TetaDec) * Lambda * DeltaT;
  486.   MESS 'Le nombre de Broc (décroissance) atteint ' Broc;
  487.   SI (Broc > 1.);
  488.     MESS 'ATTENTION, il est supérieur à 1';
  489.     MESS 'Il faut prendre un pas de temps d au plus ' 
  490.          (1. / (1. - TetaDec) / Lambda) ' ans. ';
  491.   SINON;
  492.     MESS 'Il est bien inférieur à 1';
  493.   FINSI;
  494.  
  495.   MESS 'De plus, dt/Tvie vaut au plus ' (DeltaT / TVie) ;
  496.   SI (DeltaT > TVie);
  497.     MESS 'ATTENTION, on rappelle que pour un calcul précis, il faut'
  498.          ' dt < TVie';
  499.   SINON;
  500.     MESS 'On a bien dt < TVie';
  501.   FINSI;
  502.   KDecr1 = Lambda * DeltaT;
  503.   Err1 =
  504.      ( (1. - ((1. - TetaDec) * KDecr1)) / (1. + (TetaDec * KDecr1)) )
  505.      - (EXP (-1. * KDecr1) ) * 100.;
  506.   KDecr2 = Lambda * DeltaT;
  507.   Err2 = 
  508.      ( (1. - ((1. - TetaDec) * KDecr2)) / (1. + (TetaDec * KDecr2)) )
  509.      - (EXP (-1. * KDecr2)) * 100.;
  510.   MESS 'Compte tenu des pas de temps, une mesure de l erreur '
  511.        'typique serait';
  512.   MESS ' comprise entre' Err1 'et' Err2 '% de C0';
  513.   MESS ' ';
  514.  
  515. FINSI;
  516.  
  517. *                                                        ==============
  518. *                                                        | Résolution |
  519. *                                                        ==============
  520. *-- Table de transport :
  521. Transp                         = 'TABLE';
  522. Transp.'TEMPS'                 = 'TABLE';
  523. Transp.'TRACE_CONC'            = 'TABLE';
  524. Transp.'CONCENTRATION'         = 'TABLE';
  525. Transp.'FLUX'                  = 'TABLE';
  526. SI Solub;
  527.   Transp.'PRECIPITE'           = 'TABLE';
  528.   Transp.'DISSOLUTION'         = 'TABLE';
  529. FINSI;
  530.  
  531. Transp.'SOUSTYPE'              = 'DARCY_TRANSITOIRE';
  532. Transp.'MODELE'                = MoTransp;
  533. Transp.'CARACTERISTIQUES'      = MaTransp;
  534. Transp.'POROSITE'              = MaPor;
  535. Transp.'CONVECTION'            = QFHydro;
  536. Transp.'DECROISSANCE'          = Lambda;
  537. SI SoluL;
  538.   Transp.'LIMITE_SOLUBILITE'   = MaSol1;
  539. FINSI;
  540.  
  541. * Conditions initiales :
  542. Transp.'TEMPS'. 0              = TIni;
  543. Transp.'TRACE_CONC'. 0         = TC0;
  544. Transp.'CONCENTRATION'. 0      = Conc0;
  545. Transp.'FLUX'. 0               = QFDif0;
  546. SI Solub;
  547.   Transp.'PRECIPITE'. 0        = Prec0;
  548.   Transp.'DISSOLUTION'. 0      = Libe0;
  549. FINSI;
  550.  
  551. * Conditions aux limites :
  552. Transp.'BLOCAGE'               = RigCL;
  553. Transp.'TRACE_IMPOSE'          = ChgCL;
  554.  
  555. * Paramètres numériques :
  556. Transp.'THETA_DIFF'            = TetaDif;
  557. Transp.'THETA_CONVECTION'      = TetaCnv;
  558. SI Decro;
  559.   Transp.'THETA_DEC'           = TetaDec;
  560. FINSI;
  561. SI ('EGA' Methode 'Pénalisation');
  562.   Transp.'PENALISATION'        = Penalty;
  563.   Transp.'ITMAX_LIM'           = NbItMax;
  564. FINSI;
  565. Transp.'TEMPS_CALCULES'        = LiTCalc;
  566. Transp.'TEMPS_SAUVES'          = LiTSauv;
  567.  
  568. *                                                            ==========
  569. *                                                            | CALCUL |
  570. *                                                            ==========
  571.  
  572. DARCYTRA Transp;
  573.  
  574. *                                                   ===================
  575. *                                                   | POST-TRAITEMENT |
  576. *                                                   ===================
  577.  
  578. MESS ' ';
  579. MESS 'Post-Traitement :';
  580. MESS '-----------------';
  581. MESS ' ';
  582.  
  583. IndTSauv = 'INDEX' (Transp.temps);
  584.  
  585. *                                                   Préparation tracé :
  586. *                                                   -------------------
  587. * Calcul du flux à travers une certaine ligne :
  588. * ---------------------------------------------
  589. * Les propriétés de la ligne sont stockées dans DFlux et NleFlux.
  590. * Le nom de la ligne doit etre mis dans DFluxTIT .
  591. NORHydro = 'DOMA' MoTransp 'NORMALE' ;
  592. Ori2     = 'PSCAL' ('REDU' NORHydro  CeFlux ) NleFlux
  593.                ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY');
  594. FlInt    = 0.;
  595. liabs2   = 'ENLEVER' litcalc 1;
  596. 'REPETER' Bcl21 (NTSauv - 1);
  597.   i   = &Bcl21 + 1;
  598.   Id  = IndTSauv.i;
  599.   Tps = Transp.Temps.Id;
  600.   Id1 = IndTSauv.(i - 1);
  601.   Tp1 = Transp.Temps.(id - 1);
  602.   dt  = Tps - Tp1;
  603.  
  604. *-- Flux diffusif et convectif à l'instant d'indice Id 
  605. * (diffusion implicite TetaDif et convection implicite TetaCnv) :
  606. * méthode par produit scalaire flux * normale sur 1 ligne
  607.  
  608. * Débit diffusif sortant :
  609.   DebD1 = ((1. - TetaDif) * Transp.Flux.Id1)
  610.         + (TetaDif * Transp.Flux.Id);
  611.   DebD2 = 'REDU' DebD1 CeFlux ;
  612.   DebD3 = 'RESULT' (DebD2 * Ori2);
  613.   FlD   = ('MAXIMUM' DebD3);
  614.   Fl    = FlD;
  615.  
  616. * Débit convectif sortant :
  617.   THCnv = ((1. - TetaCnv) * Transp.Trace_Conc.Id1)  
  618.         + (TetaCnv * Transp.Trace_Conc.Id);
  619.   DebC1 = ('EXCO' 'TH' THCnv 'SCAL') * Transp.Convection;
  620.   DebC2 = 'REDU' DebC1 CeFlux ;
  621.   DebC3 = 'RESULT' (DebC2 * Ori2);
  622.   FlC   = ('MAXIMUM' DebC3);
  623.   Fl  = FlC + FlD;
  624.  
  625.   FlInt = FlInt + (Fl * dt); 
  626.  
  627.   SI (i 'EGA' 2);
  628.     LiFlC = 'PROG' FlC;
  629.     LiFlD = 'PROG' FlD;
  630.     LiFl  = 'PROG' Fl;
  631.     LiFlI = 'PROG' FlInt;
  632.   SINON;
  633.     LiFlC = 'INSERER' LiFlC (i - 1) FlC;
  634.     LiFlD = 'INSERER' LiFlD (i - 1) FlD;
  635.     LiFl  = 'INSERER' LiFl  (i - 1) Fl;
  636.     LiFlI = 'INSERER' LiFlI (i - 1) FlInt;
  637.   FINSI;
  638. 'FIN' Bcl21;
  639.  
  640. EvFlC  = 'EVOL' 'JAUN' 'MANU' 't (ans)' LiAbs2 'Debit_Conv' LiFlC;
  641. EvFlD  = 'EVOL' 'BLEU' 'MANU' 't (ans)' LiAbs2 'Debit_Diff' LiFlD;
  642. EvFl   = 'EVOL' 'VERT' 'MANU' 't (ans)' LiAbs2 'Debit'      LiFl;
  643. EvFlI  = 'EVOL' 'VERT' 'MANU' 't (ans)' LiAbs2 'Debit_int.' LiFlI;
  644.  
  645. *                                                             Tracé :
  646. *                                                             ------
  647. * Décalage des champs les uns par rapport aux autres :
  648. Decal1 = (1.2 * LXM) 0.;
  649. Decal2 = (2.4 * LXM) 0.;
  650.  
  651. * Tracé du champ de limite de solubilité :
  652. * ----------------------------------------
  653. SI (TrLimSol ET SoluL);
  654.   'OPTION' ISOV SURF;
  655.   'TITRE'  'Limite de solubilite';
  656.   MaSol2 = 'KCHA' MoTransp MaSol1 'CHAM';
  657.   'TRAC' MoTransp MaSol2 LTrame;
  658. FINSI;
  659.  
  660. * Tracé des profils de conc, prec, diss le long des lignes choisies :
  661. * ------------------------------------------------------------------- 
  662. SI TrLEvol;
  663.  
  664.   IndLEvol = 'INDEX' LEvol;
  665.   TabEvlC  = table;
  666.   SI Solub;
  667.     TabEvlP  = table;
  668.     TabEvlL  = table;
  669.   Finsi;
  670.  
  671.   'REPETER' Bcl ('DIME' LEvol);
  672. *   pour chaque ligne i où suivre les valeurs :  
  673.     i = IndLEvol.&Bcl;
  674.  
  675. *   détermination du nom de la composante à suivre et évolution de
  676. *   la limite de solubilité :
  677.     si ('EGA' (LEvol.i.'P') 'CENTRE');
  678.       mot1 = 'H';
  679.     finsi;
  680.     si ('EGA' (LEvol.i.'P') 'FACE');
  681.       mot1 = 'TH';
  682.     finsi;
  683.     si ('EGA' (LEvol.i.'P') 'SOMMET');
  684.       mot1 = 'SCAL';
  685.     finsi;
  686.  
  687. *   Evolution des concentrations :
  688.     TitDesC   = 'CHAINE' 'Concentration ' (LEvol.i.'T') ;
  689.     TabEvlC.i = TRACHIS Transp 'CONCENTRATION' (LEvol.i.'M') LiIndevl
  690.                         ('MOTS' mot1) ('MOTS' ' ') 'PREF' ' ' ;
  691.  
  692.     si (SoluL et ('NEG' (LEvol.i.'P') 'FACE'));
  693. *     Adjonction de la limite de solubilité :
  694.       si ('EGA' (LEvol.i.'P') 'CENTRE');
  695.         EvSatL = 'EVOL' 'BLAN' 'CHPO' MaSol1 'H' (LEvol.I.'M');
  696.       finsi;
  697.       si ('EGA' (LEvol.i.'P') 'SOMMET');
  698.         EvSatL = 'EVOL' 'BLAN' 'CHPO' MaSolS 'SCAL' (LEvol.I.'M');
  699.       finsi;
  700.       EvSatL = 'EXTR' EvSatL 'COUR' 1;
  701.       nn = 'DIME' TabEvlC.i ;
  702.       'REPETER' bcl2 nn ;      
  703.         j = nn - &bcl2 + 1;
  704.         TabEvlC.i.(j + 1) = TabEvlC.i.j;
  705.       'FIN' bcl2;
  706.       TabEvlC.i.1 = table;
  707.       TabEvlC.i.1 .'VALEUR'  = EvSatL ;
  708.       TabEvlC.i.1 .'LEGEND1' = 'Csat' ;
  709.       TabEvlC.i.1 .'LEGEND2' = ' ';
  710.     finsi;
  711.  
  712.     DESTRA (TabEvlC.i) 'MIMA' 'AXES' 
  713.            'TITR' TitDesC 'TITX' (LEvol.i.'L') 'TITY' 'C (mol/l)'
  714.            'YBOR' 0. CSat1 ;
  715.  
  716.     SI (Solub et ('NEG' (LEvol.i.'P') 'FACE'));;
  717. *     Evolution des précipités :
  718.       TitDesP   = 'CHAINE' 'Precipite ' (LEvol.i.'T') ;
  719.       TabEvlP.i = TRACHIS Transp 'PRECIPITE' (LEvol.i.'M') LiIndevl
  720.                   ('MOTS' mot1) ('MOTS' ' ') 'PREF' ' ' ;
  721.       DESTRA (TabEvlP.i) 'MIMA' 'AXES' 'TITR' TitDesP 
  722.              'TITX' (LEvol.i.'L') 'TITY' 'P (mol/dm3)' ;
  723. *   Evolution des dissolutions :
  724.       TitDesL   = 'CHAINE' 'Dissolution ' (LEvol.i.'T') ;
  725.       TabEvlL.i = TRACHIS Transp 'DISSOLUTION' (LEvol.i.'M') LiIndevl
  726.                   ('MOTS' mot1) ('MOTS' ' ') 'PREF' ' ' ;
  727.       DESTRA (TabEvlL.i) 'MIMA' 'AXES' 'TITR' TitDesL 
  728.              'TITX' (LEvol.i.'L') 'TITY' 'L (mol/dm3/an)' ;
  729.     Finsi;
  730.  
  731.   'FIN' bcl;
  732.  
  733. FINSI;
  734.  
  735. * Tracé des évolutions de conc, prec, diss aux points choisis :
  736. * ------------------------------------------------------------- 
  737. * (L'indice 0 est réservé à un point de la source)
  738. SI (TrPEvol ET (('DIME' ('INDEX' PEVol)) > 0));
  739.  
  740.   IndPEvol = 'INDEX' PEvol;
  741.   TabEvpC  = table;
  742.   SI Solub;
  743.     TabEvpP  = table;
  744.     TabEvpL  = table;
  745.   Finsi;
  746.  
  747.   'REPETER' Bcl ('DIME' PEvol);
  748. *   pour chaque point i où suivre les valeurs :  
  749.     i = IndPEvol.&Bcl;
  750.  
  751. *   détermination du nom de la composante à suivre et évolution de
  752. *   la limite de solubilité :
  753.     si ('EGA' (PEvol.i.'P') 'CENTRE');
  754.       mot1 = 'H';
  755.     finsi;
  756.     si ('EGA' (PEvol.i.'P') 'FACE');
  757.       mot1 = 'TH';
  758.     finsi;
  759.     si ('EGA' (PEvol.i.'P') 'SOMMET');
  760.       mot1 = 'SCAL';
  761.     finsi;
  762.  
  763. *   Evolution des concentrations :
  764.     TitDesC   = 'CHAINE' 'Concentration ' (PEvol.i.'T') ;
  765.     TabEvpC.i = TRACHIT Transp 'CONCENTRATION'
  766.                ('MANU' 'POI1' (PEvol.i.'M')) LiIndEvl
  767.                ('MOTS' mot1) 'PREF' ' ' ;
  768.  
  769.     si (SoluL et ('NEG' (PEvol.i.'P') 'FACE'));
  770. *     Adjonction de la limite de solubilité :
  771.       si ('EGA' (PEvol.i.'P') 'CENTRE');
  772.         ValCSat = MaSol1 'EXTR' 'H' (PEvol.i.'M');
  773.       finsi;
  774.       si ('EGA' (PEvol.i.'P') 'SOMMET');
  775.         ValCSat = MaSolS 'EXTR' 'SCAL' (PEvol.i.'M');
  776.       finsi;
  777.       EvSatP = 'EVOL' 'MANU' ('PROG' 0. TFin) ('PROG' 2*ValCSat);
  778.       nn = 'DIME' TabEvpC.i ;
  779.       'REPETER' bcl2 nn ;      
  780.         j = nn - &bcl2 + 1;
  781.         TabEvpC.i.(j + 1) = TabEvpC.i.j;
  782.       'FIN' bcl2;
  783.       TabEvpC.i.1 = table;
  784.       TabEvpC.i.1 .'VALEUR'  = EvSatP ;
  785.       TabEvpC.i.1 .'LEGEND1' = 'Csat' ;
  786.       TabEvpC.i.1 .'LEGEND2' = ' ';
  787.     FINSI;
  788.  
  789.     DESTRA (TabEvpC.i) 'MIMA' 'AXES' 
  790.            'TITR' TitDesC 'TITX' 't (an)' 'TITY' 'C (mol/l)'
  791.            'YBOR' 0. CSat1 'XBOR' Tini TFin ;
  792.  
  793.     SI (Solub et ('NEG' (PEvol.i.'P') 'FACE'));;
  794. *     Evolution des précipités :
  795.       TitDesP   = 'CHAINE' 'Precipite ' (PEvol.i.'T') ;
  796.       TabEvpP.i = TRACHIT Transp 'PRECIPITE'
  797.                   ('MANU' 'POI1' (PEvol.i.'M')) LiIndevl
  798.                   ('MOTS' mot1) 'PREF' ' ' ;
  799.       DESTRA (TabEvpP.i) 'MIMA' 'AXES' 'TITR' TitDesP 
  800.              'TITX' 't (an)' 'TITY' 'P (mol/dm3)'
  801.              'XBOR' Tini TFin ;
  802. *     Evolution des dissolutions :
  803.       TitDesL   = 'CHAINE' 'Dissolution ' (PEvol.i.'T') ;
  804.       TabEvpL.i = TRACHIT Transp 'DISSOLUTION'
  805.                   ('MANU' 'POI1' (PEvol.i.'M')) LiIndevl
  806.                   ('MOTS' mot1) 'PREF' ' ' ;
  807.       DESTRA (TabEvpL.i) 'MIMA' 'AXES' 'TITR' TitDesL 
  808.              'TITX' 't (an)' 'TITY' 'L (mol/dm3/an)'
  809.              'XBOR' Tini TFin ;
  810.     Finsi;
  811.  
  812.   'FIN' bcl;
  813.  
  814. FINSI;    COMM 'SI TrPEvol';
  815.  
  816. * Tracé des champs de conc, prec, aux différents instants :
  817. * ---------------------------------------------------------
  818. * On déplace les supports des chpos pour tracer les fractions de 
  819. * précipité à côté des concentrations :
  820. TrPrec = TrPrec ET Solub;
  821.  
  822. *-- On trace pour chaque temps choisi :
  823. SI (TrConc OU TrPrec);
  824.  
  825. *-- Echelles :
  826.   SI TrConc;
  827.     SI SoluL;
  828.       CMinD = 'MINIMUM' MaSol1;
  829.       CMaxD = 'MAXIMUM' MaSol1;
  830.     SINON;
  831.       CMinD = 0.;
  832.       CMaxD = 0.;
  833.     FINSI;
  834.   FINSI;
  835.  
  836.   'REPETER'  Bcl17 NTSauv;
  837.     i  = &Bcl17;
  838.     Id = IndTSauv.i;
  839.     SI TrConc;
  840.       CMinD = 'MINIMUM'
  841.               ('PROG' CMinD ('MINIMUM' Transp.Concentration.Id));
  842.       CMaxD = 'MAXIMUM'
  843.               ('PROG' CMaxD ('MAXIMUM' Transp.Concentration.Id));
  844.     FINSI;
  845.   'FIN' Bcl17;
  846.  
  847.   SI TrConc;
  848.     CMinD = 'MINIMUM' ('PROG' CMinD (CMaxD * 1.e-2));
  849.     CMinD = 'MAXIMUM' ('PROG' CMinD 0.);
  850.     EchC0 = 'PROG' 0. 'PAS' (CMaxD / 13.) CMaxD;
  851.   FINSI;
  852.  
  853.   'REPETER' Bcl3 NTGrph;
  854.     Id  = 'EXTR' LiIndgph &Bcl3;
  855.     Tps = Transp.Temps.Id;
  856.     TitGraph = ' ' ;
  857.     SI TrConc;
  858.       TitGraph = 'CHAINE' TitGraph 'Solute ';
  859.     FINSI;
  860.     SI (TrConc ET TrPrec);
  861.       TitGraph = 'CHAINE' TitGraph 'et ';
  862.     FINSI;
  863.     SI TrPrec;
  864.       TitGraph = 'CHAINE' TitGraph 'Precipite ';
  865.     FINSI;
  866.     TitGraph = 'CHAINE' TitGraph 'a t=' Tps ' ans ';
  867.     'TITRE' TitGraph;
  868.  
  869.     SI TrConc;
  870.       Conc2 = 'KCHA' MoTransp Transp.Concentration.Id 'CHAM';
  871.       MoTot = MoTransp;
  872.       Tot2  = Conc2;
  873.     FINSI;
  874.     SI TrPrec;
  875.       Prec1 = 'KCHA' MoTransp Transp.Precipite.Id 'CHAM';
  876.       MoTra2 Prec2 = MoTransp Prec1 'PLUS' Decal1;
  877.       LTra2 = LTrame 'PLUS' Decal1;
  878.       'ELIMINATION' (ltra2 'ET' ('EXTR' Prec2 'MAIL')) epselim;
  879.       SI (NON TrConc);
  880.         Tot2  = Prec2;
  881.         MoTot = MoTra2;
  882.       SINON;
  883.         Tot2  = Tot2  ET Prec2;
  884.         MoTot = MoTot ET MoTra2;
  885.       FINSI;
  886.     FINSI;
  887.     'OPTION' ISOV SURF;
  888.  
  889.     SI (('MINIMUM' Tot2) < (('MAXIMUM' Tot2) - EpsPrec));
  890.       SI (TrConc ET (NON TrPrec));
  891.         SI (('MINIMUM' Conc2) < (('MAXIMUM' Conc2) - EpsPrec));
  892.           'TRACER' Motot conc2 echc0 LTrame;
  893.         FINSI;
  894.       SINON;
  895.         SI TrConc;
  896.           'TRACER' MoTot Tot2 echc0 ;
  897.         SINON;
  898.           'TRACER' MoTot Tot2 echc0 ltra2 ;
  899.         FINSI;
  900.       FINSI;
  901.     FINSI;
  902.   'FIN' Bcl3;
  903. FINSI;
  904.  
  905. * Tracé des champs de dissolution :
  906. * ---------------------------------
  907. TrLib = TrLib ET Solub;
  908. SI TrLib;
  909.   'REPETER' Bcl22 (NTgrph - 1);
  910.     Id  = ('EXTR' LiIndgph (&Bcl22 + 1)) - 1 ;
  911.     Tps = Transp.Temps.Id;
  912.     Id2 = Id + 1;
  913.     Tp2 = Transp.Temps.Id2;
  914.     dt  = Tp2 -  Tps;
  915.     TitGraph = 'CHAINE' 'Quantite dissoute entre ' 
  916.                Tps ' et ' Tp2 ' annees (mol/m3/an)';
  917.     'TITRE' TitGraph;
  918.     Libe2 = Transp.Dissolution.Id2;
  919.     Libe3 = 'KCHA' MoTransp Libe2 'CHAM';;
  920.     SI ((('MAXIMUM' Libe3) - ('MINIMUM' Libe3)) > EpsPrec);
  921.       'OPTION' ISOV SURF;
  922.       'TRAC' MoTransp Libe3 Ltrame ;
  923.     FINSI;
  924.   'FIN' Bcl22;
  925. FINSI;
  926.  
  927. * Tracé du flux traversant la ligne choisie (en mol/l/an) :
  928. * ---------------------------------------------------------
  929. SI TrFlux ;
  930.  
  931.   TitDess = 'CHAINE' 'Flux traversant ' DFluxTIT ;
  932.   TabDess   = 'TABLE';
  933.   TabDess.1 = 'MARQ CARR' ;
  934.   TabDess.2 = 'MARQ TRIA' ;
  935.   TabDess.3 = 'MARQ PLUS' ;
  936.   TabDess.'TITRE'     = 'TABLE';
  937.   TabDess.'TITRE' . 1 = 'Debit_Convectif';
  938.   TabDess.'TITRE' . 2 = 'Debit_Diffusif';
  939.   TabDess.'TITRE' . 3 = 'Debit_Total';
  940.   EvTot = EvFlC ET EvFlD ET EvFl;
  941.   'DESSIN' EvTot 'MIMA' TabDess 'TITR' TitDess 'TITX' 't (an)' 'TITY' 
  942.                  'mol/l/an' 'LEGE' 'AXES' 'XBOR' TIni TFin;
  943. FINSI;
  944.  
  945. *                                           ===========================
  946. *                                           | TEST DE BON DEROULEMENT |
  947. *                                           ===========================
  948. * Le test se fait en comparant le flux total traversant la ligne avec
  949. * une liste obtenue par un calcul préalable.
  950. liflsav = 'PROG' .66446  1.4155 1.4363 1.1050 .77084 .51835 .35170
  951.                  .24170 .16442 .13200 ;   
  952.  
  953. *'LISTE' lifl ;
  954.  
  955. err1 = 'MAXIMUM' ('ABS' (lifl - liflsav)) / ('MAXIMUM' liflsav);
  956. mess 'Erreur = ' err1 ;
  957.  
  958. si (err1 > CRIT1);
  959.   'ERREUR' 5;
  960. sinon;
  961.   'ERREUR' 0;
  962. finsi;
  963.  
  964. Fin;
  965.  
  966.  
  967.  
  968.  
  969.  
  970.  
  971.  
  972.  
  973.  
  974.  
  975.  
  976.  
  977.  
  978.  
  979.  
  980.  
  981.  

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