Télécharger transport4.dgibi

Retour à la liste

Numérotation des lignes :

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

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