* fichier : transport2EFMH.dgibi ******************** CAS TEST : transport2.dgibi ************************ * GRAPH = 'O'; 'OPTION' ECHO 1 ; CRIT1 = 1.D-3 ; 'SAUTER' 'PAGE'; 'TITRE' 'Validation du transport'; * Petit programme testant le transport en milieu poreux avec * - limite de solubilité diminuant exponentiellement vers l'aval. * équilibre chimique soluté/précipité instantané * (méthode : point fixe sur la dissolution) * * pour un cas test 2D simple : * Une petite source baignant dans un champ de vitesse uniforme. * * P4 *---------------* P3 * | | * | | * | | * P5 |---------------| P2 * | | * | | * | | * | | * P6 *----* P9 - - - * P11 * | | | * P7 *----* P8 - - - * P10 * | | * | | * P0 *---------------* P1 * 'SAUTER' 'PAGE'; * ============== * | PARAMETRES | * ============== * * Paramètres internes : * ===================== 'OPTION' 'DIME' 2 ELEM 'QUA4' ; 'OPTION' ISOV 'SURFACE' ; 'OPTION' 'TRAC' PSC ; * Paramètres de convivialité : * ============================ * Trace les limites de solubilité ? TrLimSol= Vrai ; * Trace les profils de concentration suivant certaines lignes ? TrLEvol = Vrai ; * Trace les évolutions de concentration en certaines points ? TrPEvol = Vrai ; * Trace les champs de concentration en soluté, précipité ? TrConc = vrai ; TrPrec = vrai ; * Trace la quantité de précipité dissout ? TrLib = vrai ; * Trace le flux à travers une ligne ? TrFlux = Vrai ; SI ('EGA' GRAPH 'N'); TrLEvol = faux ; TrPEvol = faux ; TrConc = faux ; TrPrec = faux ; TrLib = faux ; TrFlux = faux ; TrLimSol= faux ; FINSI; * Grandeurs physiques : * ===================== * VITESSE DE DARCY UXImp = 0.; UYImp = 1.; * POROSITE ValPor = 0.99; *-- Facteur de forme : FacForme = 1.; * LIMITE DE SOLUBILITE *-- Limite de solubilité en mol/l (global et fond): * Cette valeur ne sert que si la dissolution est instantanée ou * d'ordre 1. CSat1 = 1. ; CSat2 = 0.01 ; Distrib = 'expo'; * DIFFUSIVITE *-- Diffusivité moléculaire du traceur en m²/s DifMolec = 2.d0; * QUANTITES INITIALES *-- Concentration en marqueur dans la source à t=0 (mol/l d'eau): CActif0 = CSat1; *-- Précipité de marqueur dans la source à t=0 (mol/dm² de milieu): PActif0 = 10000. * CSat1 ; * TEMPS DE CALCUL * initial, final, et pas de temps : TIni = 0.; TFin = 100.; DelTaT = 10. ; * METHODES NUMERIQUES *-- Coef. d'implicitation (1 = implicite, 0 = explicite) * Diffusion : TetaDif = 1. ; * Convection : TetaCnv = 1. ; *-- Coefficient d'élimination : EpsElim = .001 ; *-- Précision de la machine : EpsPrec = 1.D-14; * ============ * | MAILLAGE | * ============ * Dimensions (m) : * ---------------- LX1 = 1.; LXM = 5.; LY1 = 4.; LY2 = 5.; LY3 = 8.; LYM = 12.; * Découpages : * ------------ N = 1; NX = 'ENTIER' (N * LXM); NY = 'ENTIER' (N * LYM + 1); NX1 = 'MAXIMUM' ('LECT' ('ENTIER' (NX * LX1 / LXM)) 1); NX2 = NX - NX1; NY1 = 'MAXIMUM' ('LECT' ('ENTIER' (NY * LY1 / LYM)) 1); NY2 = 'MAXIMUM' ('LECT' ('ENTIER' (NY * (LY2 - LY1) / LYM)) 1); NY3 = 'MAXIMUM' ('LECT' ('ENTIER' (NY * (LY3 - LY2) / LYM)) 1); NY4 = 'MAXIMUM' ('LECT' ('ENTIER' (NY * (LYM - LY3) / LYM)) 1); NY = NY1 + NY2 + NY3 + NY4; * Points : * -------- p0 = 0. 0.; p1 = LXM 0.; p2 = LXM LY3; p3 = LXM LYM; p4 = 0. LYM; p5 = 0. LY3; p6 = 0. LY2; p7 = 0. LY1; p8 = LX1 LY1; p9 = LX1 LY2; p10 = LXM LY1; p11 = LXM LY2; * Droites : * --------- 'OPTION' 'COULEUR' 'BLAN'; D4D5 = 'DROITE' P4 NY4 P5; D5D6 = 'DROITE' P5 NY3 P6; D6D7 = 'DROITE' P6 NY2 P7; D7D0 = 'DROITE' P7 NY1 P0; D2D5 = 'DROITE' P2 NX P5; 'OPTION' 'COULEUR' 'TURQ'; D0D1 = 'DROITE' P0 NX P1; D1D10 = 'DROITE' P1 NY1 P10; D10D11= 'DROITE' P10 NY2 P11; D11D2 = 'DROITE' P11 NY3 P2; D2D3 = 'DROITE' P2 NY4 P3; D3D4 = 'DROITE' P3 NX P4; D7D8 = 'DROITE' P7 NX1 P8; D8D9 = 'DROITE' P8 NY2 P9; D9D6 = 'DROITE' P9 NX1 P6; D10D8 = 'DROITE' P10 NX2 P8; D11D9 = 'DROITE' P11 NX2 P9; D2D5 = 'DROITE' P2 NX P5; D8D7 = 'INVERSE' D7D8; D8D10 = 'INVERSE' D10D8; D9D8 = 'INVERSE' D8D9; D6D9 = 'INVERSE' D9D6; D9D11 = 'INVERSE' D11D9; D5D2 = 'INVERSE' D2D5; D4D3 = 'INVERSE' D3D4; *-- Pour les conditions aux limites : LInf = D0D1; Lsup = D3D4; LDroit = D1D10 ET D10D11 ET D11D2 ET D2D3; LGau = 'INVERSE' (D4D5 ET D5D6 ET D6D7 ET D7D0); *-- Pour faire la somme des flux : * pour le tranfert d'activité (loin des conditions aux limites) : LFlux = D5D2; * Surfaces maillées : * ------------------- SActif = 'ORIENTER' ('COULEUR' ('DALLER' D7D8 D8D9 D9D6 D6D7) 'ROUG') ; SFond = 'ORIENTER' ('COULEUR' ('DALLER' D2D3 D3D4 D4D5 D5D2) 'VERT') ; SCoeur = 'ORIENTER' ('COULEUR' ( ('DALLER' D0D1 D1D10 (D10D8 ET D8D7) D7D0) ET ('DALLER' D8D10 D10D11 D11D9 D9D8) ET ('DALLER' (D6D9 ET D9D11) D11D2 D2D5 D5D6) ) 'TURQ' ) ; SGran = 'COULEUR' (SCoeur et SFond) 'TURQ' ; Mtot = SActif ET SGran; LTrame = 'COULEUR' ((CONT SActif) ET (CONT SGran)) 'BLAN'; * * Elements contenant tous les points de calcul QMtot = 'CHAN' Mtot 'QUAF' ; QActif = 'CHAN' SActif 'QUAF' ; QLInf = 'CHAN' LInf 'QUAF' ; QLSup = 'CHAN' LSup 'QUAF' ; QLFlux = 'CHAN' LFlux 'QUAF' ; ELIM 0.001 ( QMtot ET QActif ET QLInf ET QLSup ET QLFlux) ; * * Domaines : * ---------- * Modèle (Mo) : MoHydro = 'MODELE' QMtot DARCY ISOTROPE ; MoActif = 'MODELE' QActif DARCY ISOTROPE ; MoInf = 'MODELE' QLInf DARCY ISOTROPE ; MoSup = 'MODELE' QLSup DARCY ISOTROPE ; MoFlux = 'MODELE' QLFlux DARCY ISOTROPE ; CeFlux = 'DOMA' MoFlux 'CENTRE' ; * Longueurs, Normales et Orientations des arêtes du domaine ; Chyb1 = 'DOMA' MoHydro 'SURFACE'; Chyb2 = 'DOMA' MoHydro 'NORMALE'; VolChPro = 'DOMA' MoHydro 'VOLUME'; * Création des normales à la ligne sur laquelle on mesure les flux : NleFlux = 'MANU' 'CHPO' CeFlux 2 'UX' 0. 'UY' +1. NATURE DIFFUS; DFluxTIT = 'CHAINE' 'la ligne y=' 'FORMAT' '(F5.0)' ('COOR' 2 P2) 'm'; * =============== * | HYDRAULIQUE | * =============== * La donnée est U. Pour obtenir VDarcy et QF, on simule une résolution * de Darcy. On prend pour cela K = 1 et grad h = U : * Matériau associé (Ma) : MaHydro = 'MATERIAU' MoHydro 'K' 1.; * Matrice hybride élémentaire : MeHydro = 'MHYB' MoHydro MaHydro; * Trace de Concentration : FaHydro = 'DOMA' MoHydro 'FACE' ; XX YY = 'COOR' FaHydro ; Un = 'MANU' 'CHPO' FaHydro 1 'SCAL' 1.; THHydro = (-1.) * (UXImp * XX) - (UYImp * YY) + Un; THHydro = 'NOMC' THHydro 'TH'; * Concentration : CeHydro = 'DOMA' MoHydro 'CENTRE' ; XX YY = 'COOR' CeHydro ; Un = 'MANU' 'CHPO' CeHydro 1 'SCAL' 1.; HHydro = (-1.) * (UXImp * XX) - (UYImp * YY) + Un; HHydro = 'NOMC' HHydro 'H'; * Flux et vitesses : QFHydro = 'HDEB' MoHydro MeHydro HHydro THHydro ; VCDarcy = 'HVIT' MoHydro QFHydro ; * Points où l'on suivra la concentration : * ---------------------------------------- PEvol = 'TABLE'; * Centre du combustible : PEvol.1 = 'TABLE'; PEvol.1 .'T' = 'CHAINE' 'dans la source' ; PEvol.1 .'P' = 'CENTRE'; PEvol.1 .'M' = ( (.5 * LX1) ((LY2 + LY1) / 2.) ); PEvol.1 .'M' = CeHydro 'POINT' 'PROC' (PEvol.1 .'M'); * Juste après la transition aval : PEvol.2 = 'TABLE'; PEvol.2 .'T' = 'CHAINE' 'post-transition aval' ; PEvol.2 .'P' = 'CENTRE'; PEvol.2 .'M' = CeHydro 'POINT' 'PROC' (P5 'PLUS' (0.05 0.05)); * Lignes où capturer les évolutions des concentrations : * ------------------------------------------------------ LEvol = 'TABLE'; 'REPETER' Bcl8 ('NBEL' LGau); pt = CeHydro 'POINT' 'PROC' ((LGau 'POINT' &Bcl8) 'PLUS' (0 0.01)); SI (&Bcl8 > 2); LGauC = LGauC ET ('MANU' 'SEG2' pt0 pt); FINSI; SI (&Bcl8 'EGA' 2); LGauC = ('MANU' 'SEG2' pt0 pt); FINSI; pt0 = pt; 'FIN' Bcl8; 'OUBLIER' pt0; 'OUBLIER' pt; LEvol.1 = 'TABLE'; LEvol.1 .M = 'COULEUR' LGauC 'BLAN'; LEvol.1 .P = 'CENTRE'; LEvol.1 .L = 'Y (m)'; LEvol.1 .T = 'CHAINE' 'longitudinale' ; * ============= * | TRANSPORT | * ============= * Modèle propre au transport : * ---------------------------- MoTransp = 'MODELE' QMtot DARCY ANISOTROPE ; * Matrice Porosité cinématique (Type ChamElem) : * ---------------------------------------------- CeTransp = 'DOMA' MoTransp 'CENTRE' ; MaPor1 = 'MANU' 'CHPO' CeTransp 1 'CK' ValPor NATURE DIFFUS; MaPor = MaPor1; * Porosité * volume de maille * ---------------------------- MaPorVol= VolChPro * ('EXCO' 'CK' MaPor1 'SCAL'); * Temps : * ------- * Temps à calculer LiTCalc = 'PROG' TIni 'PAS' DeltaT Tfin; NTCalc = 'DIME' LiTCalc; 'REMPLACER' LiTCalc NTCalc Tfin; * Ecriture de ces temps sous forme de table ; TabTCalc = 'TABLE'; 'REPETER' Bcl26 NTCalc; Id = &Bcl26 - 1; Tps = 'EXTR' LiTCalc &Bcl26; TabTCalc.Id = Tps; 'FIN' Bcl26; * Temps sauvegardés LiTSauv = LiTCalc ; NTSauv = NTCalc ; * Temps de tracage des champs LiIndGph = 'LECT' 0 5 10 ; NTGrph = 'DIME' LiIndGph; * Temps de tracage des évolutions LiTEvol = litsauv ; liindEvl = 'LECT' 0 pas 1 (ntsauv - 1) ; NTEvol = 'DIME' litevol; * Temps servant à la création des chargements TSup = 1.2 * TFin; DelTDef = 0.001 * TSup; * ==================================== * | CONSTRUCTION DES DONNEES DE BASE | * ==================================== * Dissolution / Libération ; * -------------------------- *-- Matrice Limite de solubilité ; * Variation exponentielle décroissante de la source vers l'aval : XX YY = 'COOR' CeTransp ; Sup2 = 'POINT' YY 'EGINF' LY2; MaSol2 = 'MANU' 'CHPO' Sup2 1 'H' CSat1 NATURE DIFFUS; Sup3 = 'POINT' YY 'SUPERIEUR' LY2; YY3 = 'REDU' (YY - LY2) Sup3; a = ('LOG' (CSat2 / CSat1)) / (LYM - LY2); MaSol3 = CSat1 * ('EXP' (a * YY3)); MaSol3 = 'NOMC' MaSol3 'H' NATURE DIFFUS; MaSol1 = MaSol2 ET MaSol3; * Matériau associé au transport : * ------------------------------- * Terme diffusif pur (diffusivité effective) = Porosité * Diffusivité : DifEff = DifMolec * FacForme * MaPor1; MaTransp = NOMC 'K' DifEff ; * Conditions Initiales : * ---------------------- *-- Concentration aux centres des éléments (sert à définir la quantité * initiale d'élément) : Conc01 = 'MANU' 'CHPO' CeHydro 1 'H' 0. 'NATURE' 'DISCRET'; Conc02 = 'MANU' 'CHPO'( 'DOMA' MoActif 'CENTRE') 1 'H' CActif0 'NATURE' 'DISCRET'; Conc0 = Conc01 ET Conc02; Conc0 = 'NOMC' ('KCHT' MoHydro 'SCAL' 'CENTRE' ('NOMC' Conc0 'SCAL')) 'H'; *-- Flux diffusif initial (aux centres des faces) - décoratif puisque la * diffusion est implicite : QFDif0 = 'MANU' 'CHPO' FaHydro 1 'FLUX' 0. 'NATURE' 'DISCRET'; *-- Précipité : Prec01 = 'MANU' 'CHPO' ( 'DOMA' MoActif 'CENTRE') 1 'H' PActif0 'NATURE' 'DISCRET'; Prec02 = 'MANU' 'CHPO' CeHydro 1 'H' 0. 'NATURE' 'DISCRET'; Prec0 = Prec01 ET Prec02; Prec0 = 'NOMC' ('KCHT' MoHydro 'SCAL' 'CENTRE' ('NOMC' Prec0 'SCAL')) 'H'; * Conditions aux limites (type CHARGEMENT) : * ------------------------------------------ * Concentration imposée en amont et en aval, flux nul ailleurs. Ceinf = 'DOMA' Moinf 'CENTRE' ; Cesup = 'DOMA' Mosup 'CENTRE' ; CCL = 'MANU' 'CHPO' (Ceinf et Cesup) 1 'H' 0. 'NATURE' 'DISCRET'; CharCL = CCL ; ChgCL = 'CHARG' CharCL ('EVOL' 'MANU' ('PROG' 0. TSup) ('PROG' 1. 1.)); * ============== * | Résolution | * ============== *-- Table de transport : Transp = 'TABLE'; Transp.'TEMPS' = 'TABLE'; Transp.'TRACE_CONC' = 'TABLE'; Transp.'CONCENTRATION' = 'TABLE'; Transp.'FLUX' = 'TABLE'; Transp.'PRECIPITE' = 'TABLE'; Transp.'DISSOLUTION' = 'TABLE'; Transp.'SOUSTYPE' = 'DARCY_TRANSITOIRE'; Transp.'MODELE' = MoTransp; Transp.'CARACTERISTIQUES' = MaTransp; Transp.'POROSITE' = MaPor; Transp.'COEF_RETARD' = 10.D0; *Transp.'CONVECTION' = QFHydro; Transp.'CONVECTION' = (nomc SCAL (1.D0 * QFHydro)) / (doma MOTRANSP SURFACE) ; Transp . 'CONVECTION' = Transp . CONVECTION * (doma MoTRANSP NORMALE); Transp . 'VITELEM' = 'HVIT' MOtransp QFHYDRO ; Transp.'LIMITE_SOLUBILITE' = MaSol1; Transp.'COEF_DISSOLUTION' = 1.D8; * Conditions initiales : Transp.'TEMPS'. 0 = TIni; Transp.'CONCENTRATION'. 0 = Conc0; Transp . 'FLUXDIFF' = table; Transp.'FLUXDIFF' . 0 = 0.D0 * QFDif0; Transp . 'FLUXCONV' = table; Transp.'FLUXCONV' . 0 = 0.D0 * QFDif0; Transp.'PRECIPITE'. 0 = Prec0; * Transp.'DISSOLUTION'. 0 = Libe0; * Conditions aux limites : Transp.'TRACE_IMPOSE' = ChgCL; * Paramètres numériques : Transp.'THETA_DIFF' = TetaDif; Transp.'THETA_CONVECTION' = TetaCnv; Transp.'LUMP' = VRAI; Transp.'TYPDISCRETISATION' = 'EFMH'; Transp.'DECENTR' = VRAI; TABRES = table METHINV; TABRES . 'TYPINV' = 1; TABRES . 'PRECOND' = 3; Transp . 'METHINV' = TABRES; *ransp.'THETA_DISS' = TetaDis; SI ('EGA' Methode 'Point fixe Dissolution'); Transp.'EPSI_LIM' = EpsSolub; Transp.'ITMAX_LIM' = NbItMax; FINSI; Transp.'TEMPS_CALCULES' = LiTCalc; Transp.'TEMPS_SAUVES' = LiTSauv; * ========== * | CALCUL | * ========== TRANSGEN Transp; * =================== * | POST-TRAITEMENT | * =================== MESS ' '; MESS 'Post-Traitement :'; MESS '-----------------'; MESS ' '; IndTSauv = 'INDEX' (Transp.temps); * Préparation tracé : * ------------------- * Calcul du flux à travers une certaine ligne : * --------------------------------------------- * Les propriétés de la ligne sont stockées dans DFlux et NleFlux. * Le nom de la ligne doit etre mis dans DFluxTIT . NORHydro = 'DOMA' MoTransp 'NORMALE' ; Ori2 = 'PSCAL' ('REDU' NORHydro CeFlux ) NleFlux ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY'); FlInt = 0.; liabs2 = 'ENLEVER' litcalc 1; 'REPETER' Bcl21 (NTSauv - 1); i = &Bcl21 + 1; Id = IndTSauv.i; Tps = Transp.Temps.Id; Id1 = IndTSauv.(i - 1); Tp1 = Transp.Temps.(id - 1); dt = Tps - Tp1; *-- Flux diffusif et convectif à l'instant d'indice Id * (diffusion implicite TetaDif et convection implicite TetaCnv) : * méthode par produit scalaire flux * normale sur 1 ligne * Débit diffusif sortant : DebD1 = ((1. - TetaDif) * Transp.FLUXDIFF.Id1) + (TetaDif * Transp.FLUXDIFF.Id); DebD2 = 'REDU' DebD1 CeFlux ; DebD3 = 'RESULT' (DebD2 * Ori2); FlD = ('MAXIMUM' DebD3); Fl = FlD; * Débit convectif sortant : THCnv = ((1. - TetaCnv) * Transp.FLUXCONV.Id1) + (TetaCnv * Transp.FLUXCONV.Id); DebC1 = THCnv ; DebC2 = 'REDU' DebC1 CeFlux ; DebC3 = 'RESULT' (DebC2 * Ori2); FlC = ('MAXIMUM' DebC3); Fl = FlC + FlD; FlInt = FlInt + (Fl * dt); SI (i 'EGA' 2); LiFlC = 'PROG' FlC; LiFlD = 'PROG' FlD; LiFl = 'PROG' Fl; LiFlI = 'PROG' FlInt; SINON; LiFlC = 'INSERER' LiFlC (i - 1) FlC; LiFlD = 'INSERER' LiFlD (i - 1) FlD; LiFl = 'INSERER' LiFl (i - 1) Fl; LiFlI = 'INSERER' LiFlI (i - 1) FlInt; FINSI; 'FIN' Bcl21; EvFlC = 'EVOL' 'JAUN' 'MANU' 't (ans)' LiAbs2 'Debit_Conv' LiFlC; EvFlD = 'EVOL' 'BLEU' 'MANU' 't (ans)' LiAbs2 'Debit_Diff' LiFlD; EvFl = 'EVOL' 'VERT' 'MANU' 't (ans)' LiAbs2 'Debit' LiFl; EvFlI = 'EVOL' 'VERT' 'MANU' 't (ans)' LiAbs2 'Debit_int.' LiFlI; * Tracé : * ------ * Décalage des champs les uns par rapport aux autres : Decal1 = (1.2 * LXM) 0.; Decal2 = (2.4 * LXM) 0.; * Tracé des profils de conc, prec, diss le long des lignes choisies : * ------------------------------------------------------------------- SI TrLEvol; IndLEvol = 'INDEX' LEvol; TabEvlC = table; TabEvlP = table; TabEvlL = table; 'REPETER' Bcl ('DIME' LEvol); * pour chaque ligne i où suivre les valeurs : i = IndLEvol.&Bcl; * détermination du nom de la composante à suivre et évolution de * la limite de solubilité : si ('EGA' (LEvol.i.'P') 'CENTRE'); mot1 = 'H'; finsi; si ('EGA' (LEvol.i.'P') 'FACE'); mot1 = 'TH'; finsi; si ('EGA' (LEvol.i.'P') 'SOMMET'); mot1 = 'SCAL'; finsi; * Evolution des concentrations : TitDesC = 'CHAINE' 'Concentration ' (LEvol.i.'T') ; TabEvlC.i = TRACHIS Transp 'CONCENTRATION' (LEvol.i.'M') LiIndevl ('MOTS' mot1) ('MOTS' ' ') 'PREF' ' ' ; si ( ('NEG' (LEvol.i.'P') 'FACE')); * Adjonction de la limite de solubilité : si ('EGA' (LEvol.i.'P') 'CENTRE'); EvSatL = 'EVOL' 'BLAN' 'CHPO' MaSol1 'H' (LEvol.I.'M'); finsi; si ('EGA' (LEvol.i.'P') 'SOMMET'); EvSatL = 'EVOL' 'BLAN' 'CHPO' MaSolS 'SCAL' (LEvol.I.'M'); finsi; EvSatL = 'EXTR' EvSatL 'COUR' 1; 'REPETER' bcl2 ('DIME' TabEvlC.i); j = ('DIME' TabEvlC.i) - &bcl2 + 1; TabEvlC.i.(j + 1) = TabEvlC.i.j; 'FIN' bcl2; TabEvlC.i.1 = table; TabEvlC.i.1 .'VALEUR' = EvSatL ; TabEvlC.i.1 .'LEGEND1' = 'Csat' ; TabEvlC.i.1 .'LEGEND2' = ' '; finsi; * DESTRA (TabEvlC.i) 'MIMA' 'AXES' * 'TITR' TitDesC 'TITX' (LEvol.i.'L') 'TITY' 'C (mol/l)' * 'YBOR' 0. CSat1 ; SI ( ('NEG' (LEvol.i.'P') 'FACE'));; * Evolution des précipités : TitDesP = 'CHAINE' 'Precipite ' (LEvol.i.'T') ; TabEvlP.i = TRACHIS Transp 'PRECIPITE' (LEvol.i.'M') LiIndevl ('MOTS' mot1) ('MOTS' ' ') 'PREF' ' ' ; DESTRA (TabEvlP.i) 'MIMA' 'AXES' 'TITR' TitDesP 'TITX' (LEvol.i.'L') 'TITY' 'P (mol/dm3)' ; Finsi; 'FIN' bcl; FINSI; * Tracé des évolutions de conc, prec, diss aux points choisis : * ------------------------------------------------------------- * (L'indice 0 est réservé à un point de la source) SI (TrPEvol ET (('DIME' ('INDEX' PEVol)) > 0)); IndPEvol = 'INDEX' PEvol; TabEvpC = table; TabEvpP = table; TabEvpL = table; 'REPETER' Bcl ('DIME' PEvol); * pour chaque point i où suivre les valeurs : i = IndPEvol.&Bcl; * détermination du nom de la composante à suivre et évolution de * la limite de solubilité : si ('EGA' (PEvol.i.'P') 'CENTRE'); mot1 = 'H'; finsi; si ('EGA' (PEvol.i.'P') 'FACE'); mot1 = 'TH'; finsi; si ('EGA' (PEvol.i.'P') 'SOMMET'); mot1 = 'SCAL'; finsi; * Evolution des concentrations : TitDesC = 'CHAINE' 'Concentration ' (PEvol.i.'T') ; TabEvpC.i = TRACHIT Transp 'CONCENTRATION' ('MANU' 'POI1' (PEvol.i.'M')) LiIndEvl ('MOTS' mot1) 'PREF' ' ' ; si (('NEG' (PEvol.i.'P') 'FACE')); * Adjonction de la limite de solubilité : si ('EGA' (PEvol.i.'P') 'CENTRE'); ValCSat = MaSol1 'EXTR' 'H' (PEvol.i.'M'); finsi; si ('EGA' (PEvol.i.'P') 'SOMMET'); ValCSat = MaSolS 'EXTR' 'SCAL' (PEvol.i.'M'); finsi; EvSatP = 'EVOL' 'MANU' ('PROG' 0. TFin) ('PROG' 2*ValCSat); 'REPETER' bcl2 ('DIME' TabEvpC.i); j = ('DIME' TabEvpC.i) - &bcl2 + 1; TabEvpC.i.(j + 1) = TabEvpC.i.j; 'FIN' bcl2; TabEvpC.i.1 = table; TabEvpC.i.1 .'VALEUR' = EvSatP ; TabEvpC.i.1 .'LEGEND1' = 'Csat' ; TabEvpC.i.1 .'LEGEND2' = ' '; FINSI; * DESTRA (TabEvpC.i) 'MIMA' 'AXES' * 'TITR' TitDesC 'TITX' 't (an)' 'TITY' 'C (mol/l)' * 'YBOR' 0. CSat1 'XBOR' Tini TFin ; SI ( ('NEG' (PEvol.i.'P') 'FACE'));; * Evolution des précipités : TitDesP = 'CHAINE' 'Precipite ' (PEvol.i.'T') ; TabEvpP.i = TRACHIT Transp 'PRECIPITE' ('MANU' 'POI1' (PEvol.i.'M')) LiIndevl ('MOTS' mot1) 'PREF' ' ' ; * DESTRA (TabEvpP.i) 'MIMA' 'AXES' 'TITR' TitDesP * 'TITX' 't (an)' 'TITY' 'P (mol/dm3)' * 'XBOR' Tini TFin ; Finsi; 'FIN' bcl; FINSI; COMM 'SI TrPEvol'; * Tracé des champs de conc, prec, aux différents instants : * --------------------------------------------------------- * On déplace les supports des chpos pour tracer les fractions de * précipité à côté des concentrations : TrPrec = TrPrec; *-- On trace pour chaque temps choisi : SI (TrConc OU TrPrec); *-- Echelles : SI TrConc; CMinD = 'MINIMUM' MaSol1; CMaxD = 'MAXIMUM' MaSol1; SINON; CMinD = 0.; CMaxD = 0.; FINSI; FINSI; 'REPETER' Bcl17 NTSauv; i = &Bcl17; Id = IndTSauv.i; SI TrConc; CMinD = 'MINIMUM' ('PROG' CMinD ('MINIMUM' Transp.Concentration.Id)); CMaxD = 'MAXIMUM' ('PROG' CMaxD ('MAXIMUM' Transp.Concentration.Id)); FINSI; 'FIN' Bcl17; SI TrConc; CMinD = 'MINIMUM' ('PROG' CMinD (CMaxD * 1.e-2)); CMinD = 'MAXIMUM' ('PROG' CMinD 0.); EchC0 = 'PROG' 0. 'PAS' (CMaxD / 13.) CMaxD; FINSI; 'REPETER' Bcl3 NTGrph; Id = 'EXTR' LiIndgph &Bcl3; Tps = Transp.Temps.Id; TitGraph = ' ' ; SI TrConc; TitGraph = 'CHAINE' TitGraph 'Solute '; FINSI; SI (TrConc ET TrPrec); TitGraph = 'CHAINE' TitGraph 'et '; FINSI; SI TrPrec; TitGraph = 'CHAINE' TitGraph 'Precipite '; FINSI; TitGraph = 'CHAINE' TitGraph 'a t=' Tps ' ans '; 'TITRE' TitGraph; Conc2 = 'KCHA' MoTransp Transp.Concentration.Id 'CHAM'; MoTot = MoTransp; Tot2 = Conc2; SI TrPrec; Prec1 = 'KCHA' MoTransp Transp.Precipite.Id 'CHAM'; MoTra2 Prec2 = MoTransp Prec1 'PLUS' Decal1; LTra2 = LTrame 'PLUS' Decal1; 'ELIMINATION' (ltra2 'ET' ('EXTR' Prec2 'MAIL')) epselim; SI (NON TrConc); Tot2 = Prec2; MoTot = MoTra2; SINON; Tot2 = Tot2 ET Prec2; MoTot = MoTot ET MoTra2; FINSI; FINSI; 'OPTION' ISOV SURF; SI (('MINIMUM' Tot2) < (('MAXIMUM' Tot2) - EpsPrec)); SI (TrConc ET (NON TrPrec)); SI (('MINIMUM' Conc2) < (('MAXIMUM' Conc2) - EpsPrec)); 'TRACER' Motot conc2 echc0 LTrame; FINSI; SINON; SI TrConc; 'TRACER' MoTot Tot2 echc0; SINON; 'TRACER' MoTot Tot2 echc0 ltra2 ; FINSI; FINSI; FINSI; 'FIN' Bcl3; FINSI; * Tracé des champs de dissolution : * --------------------------------- TrLib = TrLib; SI TrLib; 'REPETER' Bcl22 (NTgrph - 1); Id = ('EXTR' LiIndgph (&Bcl22 + 1)) - 1 ; Tps = Transp.Temps.Id; Id2 = Id + 1; Tp2 = Transp.Temps.Id2; dt = Tp2 - Tps; TitGraph = 'CHAINE' 'Quantite dissoute entre ' Tps ' et ' Tp2 ' annees (mol/m3/an)'; 'TITRE' TitGraph; 'FIN' Bcl22; FINSI; * Tracé du flux traversant la ligne choisie : * ------------------------------------------ SI TrFlux ; TitDess = 'CHAINE' 'Flux traversant ' DFluxTIT ; TabDess = 'TABLE'; TabDess.1 = 'MARQ CARR' ; TabDess.2 = 'MARQ TRIA' ; TabDess.3 = 'MARQ PLUS' ; TabDess.'TITRE' = 'TABLE'; TabDess.'TITRE' . 1 = 'Debit_Convectif'; TabDess.'TITRE' . 2 = 'Debit_Diffusif'; TabDess.'TITRE' . 3 = 'Debit_Total'; EvTot = EvFlC ET EvFlD ET EvFl; 'DESSIN' EvTot 'MIMA' TabDess 'TITR' TitDess 'TITX' 't (an)' 'TITY' 'mol/l/an' 'LEGE' 'AXES' 'XBOR' TIni TFin; FINSI; * =========================== * | TEST DE BON DEROULEMENT | * =========================== * Le test se fait en comparant le flux total traversant la ligne avec * une liste obtenue par un calcul préalable. liflsav = 'PROG' 0.56061 1.0216 1.2444 1.3013 1.2365 1.0595 0.93058 0.64840 0.41501 0.29818; err1 = 'MAXIMUM' ('ABS' (lifl - liflsav)) / ('MAXIMUM' liflsav); 'MESSAGE' 'difference avec sol de ref' err1; si (err1 > CRIT1); 'ERREUR' 5; sinon; 'ERREUR' 0; finsi; 'LISTE' crit1; Fin;