Télécharger la plaquette Cast3M


fr en

Télécharger transport2.dgibi

Retour à la liste

Numérotation des lignes : 

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

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