Télécharger transport3.dgibi

Retour à la liste

Numérotation des lignes :

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

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