Télécharger la plaquette Cast3M


fr en

Télécharger transport6.dgibi

Retour à la liste

Numérotation des lignes : 

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

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