Télécharger la plaquette Cast3M


fr en

Télécharger transport6EFMH.dgibi

Retour à la liste

Numérotation des lignes : 

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

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