Télécharger transport6VF.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : transport6VF.dgibi
  2. ******************** CAS TEST : transport6.dgibi ************************
  3. *
  4. GRAPH = 'O';
  5. 'OPTION' ECHO 1 ;
  6. CRIT1 = 5.D-2 ;
  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' = 'VF';
  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. 'TITRE' 'Limite de solubilite';
  638. MaSol2 = 'KCHA' MoTransp MaSol1 'CHAM';
  639. 'TRAC' MoTransp MaSol2 LTrame;
  640. FINSI;
  641.  
  642. * Tracé des profils de conc, prec, diss le long des lignes choisies :
  643. * -------------------------------------------------------------------
  644. SI TrLEvol;
  645.  
  646. IndLEvol = 'INDEX' LEvol;
  647. TabEvlC = table;
  648. SI Solub;
  649. TabEvlP = table;
  650. TabEvlL = table;
  651. Finsi;
  652.  
  653. 'REPETER' Bcl ('DIME' LEvol);
  654. * pour chaque ligne i où suivre les valeurs :
  655. i = IndLEvol.&Bcl;
  656.  
  657. * détermination du nom de la composante à suivre et évolution de
  658. * la limite de solubilité :
  659. si ('EGA' (LEvol.i.'P') 'CENTRE');
  660. mot1 = 'H';
  661. finsi;
  662. si ('EGA' (LEvol.i.'P') 'FACE');
  663. mot1 = 'TH';
  664. finsi;
  665. si ('EGA' (LEvol.i.'P') 'SOMMET');
  666. mot1 = 'SCAL';
  667. finsi;
  668.  
  669. * Evolution des concentrations :
  670. TitDesC = 'CHAINE' 'Concentration ' (LEvol.i.'T') ;
  671. TabEvlC.i = TRACHIS Transp 'CONCENTRATION' (LEvol.i.'M') LiIndevl
  672. ('MOTS' mot1) ('MOTS' ' ') 'PREF' ' ' ;
  673.  
  674. si (SoluL et ('NEG' (LEvol.i.'P') 'FACE'));
  675. * Adjonction de la limite de solubilité :
  676. si ('EGA' (LEvol.i.'P') 'CENTRE');
  677. EvSatL = 'EVOL' 'BLAN' 'CHPO' MaSol1 'H' (LEvol.I.'M');
  678. finsi;
  679. si ('EGA' (LEvol.i.'P') 'SOMMET');
  680. EvSatL = 'EVOL' 'BLAN' 'CHPO' MaSolS 'SCAL' (LEvol.I.'M');
  681. finsi;
  682. EvSatL = 'EXTR' EvSatL 'COUR' 1;
  683. 'REPETER' bcl2 ('DIME' TabEvlC.i);
  684. j = ('DIME' TabEvlC.i) - &bcl2 + 1;
  685. TabEvlC.i.(j + 1) = TabEvlC.i.j;
  686. 'FIN' bcl2;
  687. TabEvlC.i.1 = table;
  688. TabEvlC.i.1 .'VALEUR' = EvSatL ;
  689. TabEvlC.i.1 .'LEGEND1' = 'Csat' ;
  690. TabEvlC.i.1 .'LEGEND2' = ' ';
  691. finsi;
  692.  
  693. * DESTRA (TabEvlC.i) 'MIMA' 'AXES'
  694. * 'TITR' TitDesC 'TITX' (LEvol.i.'L') 'TITY' 'C (mol/l)'
  695. * 'YBOR' 0. CSat1 ;
  696.  
  697. SI (Solub et ('NEG' (LEvol.i.'P') 'FACE'));
  698. * Evolution des précipités :
  699. TitDesP = 'CHAINE' 'Precipite ' (LEvol.i.'T') ;
  700. TabEvlP.i = TRACHIS Transp 'PRECIPITE' (LEvol.i.'M') LiIndevl
  701. ('MOTS' mot1) ('MOTS' ' ') 'PREF' ' ' ;
  702. DESTRA (TabEvlP.i) 'MIMA' 'AXES' 'TITR' TitDesP
  703. 'TITX' (LEvol.i.'L') 'TITY' 'P (mol/dm3)' ;
  704. * Evolution des dissolutions :
  705. TitDesL = 'CHAINE' 'Dissolution ' (LEvol.i.'T') ;
  706. TabEvlL.i = TRACHIS Transp 'DISSOLUTION' (LEvol.i.'M') LiIndevl
  707. ('MOTS' mot1) ('MOTS' ' ') 'PREF' ' ' ;
  708. DESTRA (TabEvlL.i) 'MIMA' 'AXES' 'TITR' TitDesL
  709. 'TITX' (LEvol.i.'L') 'TITY' 'L (mol/dm3/an)' ;
  710. Finsi;
  711.  
  712. 'FIN' bcl;
  713.  
  714. FINSI;
  715.  
  716. * Tracé des évolutions de conc, prec, diss aux points choisis :
  717. * -------------------------------------------------------------
  718. * (L'indice 0 est réservé à un point de la source)
  719. SI (TrPEvol ET (('DIME' ('INDEX' PEVol)) > 0));
  720.  
  721. IndPEvol = 'INDEX' PEvol;
  722. TabEvpC = table;
  723. SI Solub;
  724. TabEvpP = table;
  725. TabEvpL = table;
  726. Finsi;
  727.  
  728. 'REPETER' Bcl ('DIME' PEvol);
  729. * pour chaque point i où suivre les valeurs :
  730. i = IndPEvol.&Bcl;
  731.  
  732. * détermination du nom de la composante à suivre et évolution de
  733. * la limite de solubilité :
  734. si ('EGA' (PEvol.i.'P') 'CENTRE');
  735. mot1 = 'H';
  736. finsi;
  737. si ('EGA' (PEvol.i.'P') 'FACE');
  738. mot1 = 'TH';
  739. finsi;
  740. si ('EGA' (PEvol.i.'P') 'SOMMET');
  741. mot1 = 'SCAL';
  742. finsi;
  743.  
  744. * Evolution des concentrations :
  745. TitDesC = 'CHAINE' 'Concentration ' (PEvol.i.'T') ;
  746. TabEvpC.i = TRACHIT Transp 'CONCENTRATION'
  747. ('MANU' 'POI1' (PEvol.i.'M')) LiIndEvl
  748. ('MOTS' mot1) 'PREF' ' ' ;
  749.  
  750. si (SoluL et ('NEG' (PEvol.i.'P') 'FACE'));
  751. * Adjonction de la limite de solubilité :
  752. si ('EGA' (PEvol.i.'P') 'CENTRE');
  753. ValCSat = MaSol1 'EXTR' 'H' (PEvol.i.'M');
  754. finsi;
  755. si ('EGA' (PEvol.i.'P') 'SOMMET');
  756. ValCSat = MaSolS 'EXTR' 'SCAL' (PEvol.i.'M');
  757. finsi;
  758. EvSatP = 'EVOL' 'MANU' ('PROG' 0. TFin) ('PROG' 2*ValCSat);
  759. 'REPETER' bcl2 ('DIME' TabEvpC.i);
  760. j = ('DIME' TabEvpC.i) - &bcl2 + 1;
  761. TabEvpC.i.(j + 1) = TabEvpC.i.j;
  762. 'FIN' bcl2;
  763. TabEvpC.i.1 = table;
  764. TabEvpC.i.1 .'VALEUR' = EvSatP ;
  765. TabEvpC.i.1 .'LEGEND1' = 'Csat' ;
  766. TabEvpC.i.1 .'LEGEND2' = ' ';
  767. FINSI;
  768.  
  769. * DESTRA (TabEvpC.i) 'MIMA' 'AXES'
  770. * 'TITR' TitDesC 'TITX' 't (an)' 'TITY' 'C (mol/l)'
  771. * 'YBOR' 0. CSat1 'XBOR' Tini TFin ;
  772.  
  773. SI (Solub et ('NEG' (PEvol.i.'P') 'FACE'));
  774. * Evolution des précipités :
  775. TitDesP = 'CHAINE' 'Precipite ' (PEvol.i.'T') ;
  776. TabEvpP.i = TRACHIT Transp 'PRECIPITE'
  777. ('MANU' 'POI1' (PEvol.i.'M')) LiIndevl
  778. ('MOTS' mot1) 'PREF' ' ' ;
  779. * DESTRA (TabEvpP.i) 'MIMA' 'AXES' 'TITR' TitDesP
  780. * 'TITX' 't (an)' 'TITY' 'P (mol/dm3)'
  781. * 'XBOR' Tini TFin;
  782. * Evolution des dissolutions :
  783. TitDesL = 'CHAINE' 'Dissolution ' (PEvol.i.'T') ;
  784. TabEvpL.i = TRACHIT Transp 'DISSOLUTION'
  785. ('MANU' 'POI1' (PEvol.i.'M')) LiIndevl
  786. ('MOTS' mot1) ('MOTS' ' ') 'PREF' ' ' ;
  787. * DESTRA (TabEvpL.i) 'MIMA' 'AXES' 'TITR' TitDesL
  788. * 'TITX' 't (an)' 'TITY' 'L (mol/dm3/an)'
  789. * 'XBOR' Tini TFin;
  790. Finsi;
  791.  
  792. 'FIN' bcl;
  793.  
  794. FINSI; COMM 'SI TrPEvol';
  795.  
  796. * Tracé des champs de conc, prec, aux différents instants :
  797. * ---------------------------------------------------------
  798. * On déplace les supports des chpos pour tracer les fractions de
  799. * précipité à côté des concentrations :
  800. TrPrec = TrPrec ET Solub;
  801.  
  802. *-- On trace pour chaque temps choisi :
  803. SI (TrConc OU TrPrec);
  804.  
  805. *-- Echelles :
  806. SI TrConc;
  807. SI SoluL;
  808. CMinD = 'MINIMUM' MaSol1;
  809. CMaxD = 'MAXIMUM' MaSol1;
  810. SINON;
  811. CMinD = 0.;
  812. CMaxD = 0.;
  813. FINSI;
  814. FINSI;
  815.  
  816. 'REPETER' Bcl17 NTSauv;
  817. i = &Bcl17;
  818. Id = IndTSauv.i;
  819. SI TrConc;
  820. CMinD = 'MINIMUM'
  821. ('PROG' CMinD ('MINIMUM' Transp.Concentration.Id));
  822. CMaxD = 'MAXIMUM'
  823. ('PROG' CMaxD ('MAXIMUM' Transp.Concentration.Id));
  824. FINSI;
  825. 'FIN' Bcl17;
  826.  
  827. SI TrConc;
  828. CMinD = 'MINIMUM' ('PROG' CMinD (CMaxD * 1.e-2));
  829. CMinD = 'MAXIMUM' ('PROG' CMinD 0.);
  830. EchC0 = 'PROG' 0. 'PAS' (CMaxD / 13.) CMaxD;
  831. FINSI;
  832.  
  833. 'REPETER' Bcl3 NTGrph;
  834. Id = 'EXTR' LiIndgph &Bcl3;
  835. Tps = Transp.Temps.Id;
  836. TitGraph = ' ' ;
  837. SI TrConc;
  838. TitGraph = 'CHAINE' TitGraph 'Solute ';
  839. FINSI;
  840. SI (TrConc ET TrPrec);
  841. TitGraph = 'CHAINE' TitGraph 'et ';
  842. FINSI;
  843. SI TrPrec;
  844. TitGraph = 'CHAINE' TitGraph 'Precipite ';
  845. FINSI;
  846. TitGraph = 'CHAINE' TitGraph 'a t=' Tps ' ans ';
  847. 'TITRE' TitGraph;
  848.  
  849. SI TrConc;
  850. Conc2 = 'KCHA' MoTransp Transp.Concentration.Id 'CHAM';
  851. MoTot = MoTransp;
  852. Tot2 = Conc2;
  853. FINSI;
  854. SI TrPrec;
  855. Prec1 = 'KCHA' MoTransp Transp.Precipite.Id 'CHAM';
  856. MoTra2 Prec2 = MoTransp Prec1 'PLUS' Decal1;
  857. LTra2 = LTrame 'PLUS' Decal1;
  858. 'ELIMINATION' (ltra2 'ET' ('EXTR' Prec2 'MAIL')) epselim;
  859. SI (NON TrConc);
  860. Tot2 = Prec2;
  861. MoTot = MoTra2;
  862. SINON;
  863. Tot2 = Tot2 ET Prec2;
  864. MoTot = MoTot ET MoTra2;
  865. FINSI;
  866. FINSI;
  867.  
  868. SI (('MINIMUM' Tot2) < (('MAXIMUM' Tot2) - EpsPrec));
  869. SI (TrConc ET (NON TrPrec));
  870. SI (('MINIMUM' Conc2) < (('MAXIMUM' Conc2) - EpsPrec));
  871. 'TRACER' Motot conc2 echc0 LTrame;
  872. FINSI;
  873. SINON;
  874. SI TrConc;
  875. 'TRACER' MoTot Tot2 echc0;
  876. SINON;
  877. 'TRACER' MoTot Tot2 echc0 ltra2;
  878. FINSI;
  879. FINSI;
  880. FINSI;
  881. 'FIN' Bcl3;
  882. FINSI;
  883.  
  884. * Tracé des champs de dissolution :
  885. * ---------------------------------
  886. TrLib = TrLib ET Solub;
  887. SI TrLib;
  888. 'REPETER' Bcl22 (NTgrph - 1);
  889. Id = ('EXTR' LiIndgph (&Bcl22 + 1)) - 1 ;
  890. Tps = Transp.Temps.Id;
  891. Id2 = Id + 1;
  892. Tp2 = Transp.Temps.Id2;
  893. dt = Tp2 - Tps;
  894. TitGraph = 'CHAINE' 'Quantite dissoute entre '
  895. Tps ' et ' Tp2 ' annees (mol/m3/an)';
  896. 'TITRE' TitGraph;
  897. Libe2 = Transp.Dissolution.Id2;
  898. Libe3 = 'KCHA' MoTransp Libe2 'CHAM';
  899. SI ((('MAXIMUM' Libe3) - ('MINIMUM' Libe3)) > EpsPrec);
  900. 'TRAC' MoTransp Libe3 Ltrame ;
  901. FINSI;
  902. 'FIN' Bcl22;
  903. FINSI;
  904.  
  905. * Tracé du flux traversant la ligne choisie :
  906. * ------------------------------------------
  907. SI TrFlux ;
  908.  
  909. TitDess = 'CHAINE' 'Flux traversant ' DFluxTIT ;
  910. TabDess = 'TABLE';
  911. TabDess.1 = 'MARQ CARR' ;
  912. TabDess.2 = 'MARQ TRIA' ;
  913. TabDess.3 = 'MARQ PLUS' ;
  914. TabDess.'TITRE' = 'TABLE';
  915. TabDess.'TITRE' . 1 = 'Debit_Convectif';
  916. TabDess.'TITRE' . 2 = 'Debit_Diffusif';
  917. TabDess.'TITRE' . 3 = 'Debit_Total';
  918. EvTot = EvFlC ET EvFlD ET EvFl;
  919. 'DESSIN' EvTot 'MIMA' TabDess 'TITR' TitDess 'TITX' 't (an)' 'TITY'
  920. 'mol/l/an' 'LEGE' 'AXES' 'XBOR' TIni TFin;
  921. FINSI;
  922.  
  923. * ===========================
  924. * | TEST DE BON DEROULEMENT |
  925. * ===========================
  926. * Le test se fait en comparant le flux total traversant la ligne avec
  927. * une liste obtenue par un calcul préalable.
  928. *liflsav = 'PROG' .68457 1.1205 1.3157 1.3279 1.2973 1.1393 .79427
  929. * .48185 .34038 .24040 ;
  930.  
  931. liflsav = 'PROG' 5.69132E-02 .10321 .11576 .10992 9.835E-02 8.6297E-02
  932. 7.53200E-02 6.57577E-02 5.75179E-02 5.04195E-02 ;
  933.  
  934. *'LISTE' lifl;
  935.  
  936. err1 = 'MAXIMUM' ('ABS' (lifl - liflsav)) / ('MAXIMUM' liflsav);
  937. mess 'Erreur = ' err1 ;
  938.  
  939. si (err1 > CRIT1);
  940. 'ERREUR' 5;
  941. 'OPTION' donn 5 ;
  942. sinon;
  943. 'ERREUR' 0;
  944. finsi;
  945.  
  946. Fin;
  947.  
  948.  
  949.  
  950.  
  951.  
  952.  
  953.  
  954.  
  955.  
  956.  
  957.  
  958.  
  959.  
  960.  
  961.  

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