* fichier : transport4.dgibi
******************** CAS TEST : transport4.dgibi ************************
*
GRAPH = faux ;
CRIT1 = 1.D-3 ;
'SAUTER' 'PAGE';
************************************************************************
*
* Petit programme testant le transport en milieu poreux avec
* - limite de solubilité diminuant linéairement vers l'aval.
* l'équilibre chimique soluté/précipité à cinétique linéaire
* (théta-méthode)
* - coefficient de retard constant
* (équilibre soluté/adsorbat instantané linéaire)
*
* 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
*
************************************************************************
'OPTION' 'ECHO' 0 ;
'SAUTER' 'PAGE';
'TITRE' 'Validation du transport';
* ==============
* | PARAMETRES |
* ==============
*
* Paramètres internes :
* =====================
'OPTION' 'ECHO' 0 'DIME' 2 'ELEM' 'QUA4' ;
'OPTION' 'ISOV' 'SURFACE' ;
'OPTION' 'TRAC' 'PSC' ;
* Paramètres de convivialité :
* ============================
* Trace les limites de solubilité ?
TrLimSol= vrai ;
* Trace le coefficient de dissolution ?
TrCoDis = Vrai ;
* Trace le coefficient de retard ?
TrRet = 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 ;
TrCoDis = faux ;
TrRet = faux ;
FINSI;
* Grandeurs physiques :
* =====================
* VITESSE DE DARCY
UXImp = 0.;
UYImp = 1.;
* POROSITE
ValPor = 1.;
*-- Facteur de forme :
FacForme = 1.;
* RETARD
*-- Coefficient de retard et matrice de retard
* (éventuellement scalaire) :
CofRetar = 10. ;
* 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.1 ;
* DIFFUSIVITE
*-- Diffusivité moléculaire du traceur en m²/s
DifMolec = 2.d0;
* NATURE DISSOLUTION
*-- Nature de la dissolution (si existe) et logiques associés :
*NaturDis = 'Pas de dissolution';
*NaturDis = 'Imposée';
NaturDis = 'Ordre 1';
*NaturDis = 'Instant_Pénal';
*NaturDis = 'Instant_PFDis';
Solub = 'NEG' NaturDis 'Pas de dissolution';
SoluA = 'EGA' NaturDis 'Imposée';
SoluP = 'EGA' NaturDis 'Ordre 1';
SoluI = ('EGA' NaturDis 'Instant_PFDis')
OU ('EGA' NaturDis 'Instant_Pénal');
SoluL = SoluP OU SoluI;
* PARAM. NUM. DISSOLUTION
*-- Grandeur(s) associée(s) à la dissolution :
si SoluP;
* Coefficient de dissolution
ValCoDis = LOG 2.;
finsi;
* 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):
SI Solub;
PActif0 = 10. * CSat1 * CofRetar;
FINSI;
* 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. ;
* Dissolution progressive
TetaDis = 1. ;
*-- Coefficient d'élimination :
EpsElim = .1 ;
*-- 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' ;
QFond = 'CHAN' SFond 'QUAF' ;
QActif = 'CHAN' SActif 'QUAF' ;
QGran = 'CHAN' SGran 'QUAF' ;
QInf = 'CHAN' LInf 'QUAF' ;
QDro = 'CHAN' LDroit 'QUAF' ;
QSup = 'CHAN' Lsup 'QUAF' ;
QGau = 'CHAN' LGau 'QUAF' ;
QLFlux = 'CHAN' LFlux 'QUAF' ;
ELIM 0.001 ( QMtot ET QActif ET QFond ET QGran ET QLFlux ET
QGau ET QInf ET QDro ET QSup) ;
* Domaines :
* ----------
* Modèle (Mo) :
MoHydro = 'MODELE' QMtot DARCY ISOTROPE ;
MoFond = 'MODELE' QFond DARCY ISOTROPE ;
MoActif = 'MODELE' QActif DARCY ISOTROPE ;
MoGran = 'MODELE' QGran DARCY ISOTROPE ;
MoInf = 'MODELE' QInf DARCY ISOTROPE ;
MoDro = 'MODELE' QDro DARCY ISOTROPE ;
MoSup = 'MODELE' QSup DARCY ISOTROPE ;
MoGau = 'MODELE' QGau DARCY ISOTROPE ;
MoFlux = 'MODELE' QLFlux DARCY ISOTROPE ;
* Longueurs, Normales et Orientations des arêtes du domaine ;
Chyb1 = 'DOMA' MoHydro 'SURFACE';
Chyb2 = 'DOMA' MoHydro 'NORMALE';
VolChPro = 'DOMA' MoHydro 'VOLUME';
CeHydro = 'DOMA' MoHydro 'CENTRE';
CeFlux = 'DOMA' MoFlux 'CENTRE';
* 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';
* 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' ;
* ===============
* | 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 :
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 ;
* =============
* | TRANSPORT |
* =============
* Modèle propre au transport :
* ----------------------------
MoTransp = 'MODELE' QMtot DARCY ISOTROPE ;
* Matrice Porosité cinématique (Type ChamElem) :
* ----------------------------------------------
CeTransp= 'DOMA' MoTransp 'CENTRE' ;
MaPor1 = 'MANU' 'CHPO' CeTransp 1 'CK' ValPor NATURE DIFFUS;
MaPor = 'KCHA' MoTransp MaPor1 'CHAM';
MaPorVol= VolChPro * ('EXCO' 'CK' MaPor1 'SCAL');
* Matrice coef. de retard (Type Scalaire car constant) :
* ----------------------------------------------
MaRet1 = CofRetar;
* Porosité * retard * volume de maille (x Ccentre = soluté + adsorbat par maille)
* ------------------------------------
PorVoRet = MaPorVol * MaRet1;
* 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 linéaire décroissance de l'amont vers l'aval :
XX YY = 'COOR' CeHydro;
MaSol4 = 'MANU' 'CHPO' CeHydro 1 'H' (CSat2 - CSat1)
NATURE DIFFUS;
MaSol1 = (YY * MaSol4 / LYM) + CSat1;
* Matrice Coefficient de dissolution :
* ------------------------------------
SI SoluP;
MaCoDis1 = 'MANU' 'CHPO' CeHydro 1 'SCAL' ValCoDis;
FINSI;
* Matériau associé au transport :
* -------------------------------
* Terme diffusif pur (diffusivité effective) = Porosité * Diffusivité :
DifEff = DifMolec * FacForme * MaPor;
MaTransp = 'MATERIAU' MoTransp 'K' DifEff ;
* Matrice hybride élémentaire :
* -----------------------------
MeTransp = 'MHYB' MoTransp MaTransp;
* Conditions Initiales :
* ----------------------
*-- Traces de concentration (quelconque car le transport est implicité).
TC0 = 'MANU' 'CHPO' FaHydro 1 'TH' 0. 'NATURE' 'DISCRET';
*-- 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é :
SI Solub;
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';
FINSI;
*-- Libération de soluté à partir du précipité (décoratif aussi):
SI Solub;
Libe0 = 'MANU' 'CHPO' CeHydro 1 'H' 0. 'NATURE' 'DISCRET';
FINSI;
* 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 'TH' 0.
'NATURE' 'DISCRET';
RigCL = 'BLOQUE' (CeInf et CeSup) 'TH';
CharCL = 'DEPIMPOSE' RigCL CCL ;
ChgCL = 'CHARG' CharCL
('EVOL' 'MANU' ('PROG' 0. TSup) ('PROG' 1. 1.));
* ===============================
* | Vérification des paramètres |
* | de CONVERGENCE |
* ===============================
MESS 'Vérification des paramètres de CONVERGENCE :';
MESS '--------------------------------------------';
* Cohérence dissolution progressive :
* -----------------------------------
* Pour éviter que la dissolution progressive ne soit si rapide que
* l'on dépasse régulièrement la limite de solubilité, au risque
* d'osciller, voire de diverger, vérifier que Codis * DeltaT < 1 :
SI SoluP;
Stietel = (1. - TetaDis) * MaCoDis1 * DeltaT;
SI ('NEG' (TYPE Stietel) 'CHPOINT ');
StietelM = Stietel;
SINON;
StietelM = 'MAXIMUM' Stietel;
FINSI;
MESS 'Le nombre de Stietel (dissolution progressive) atteint '
StietelM ;
SI (StietelM > 1.);
MESS 'Attention, il est supérieur à 1';
SI ('NEG' (TYPE MaCoDis1) 'CHPOINT ');
MinCoDis = MaCoDis1;
SINON;
MinCoDis = 'MINIMUM' MaCoDis1;
FINSI;
MESS 'Il faut prendre un pas de temps d au plus'
(1. / MinCoDis) 'ans. ';
SINON;
MESS 'Il est bien inférieur à 1';
FINSI;
MESS ' ';
FINSI;
* ==============
* | Résolution |
* ==============
*-- Table de transport :
Transp = 'TABLE';
Transp.'TEMPS' = 'TABLE';
Transp.'TRACE_CONC' = 'TABLE';
Transp.'CONCENTRATION' = 'TABLE';
Transp.'FLUX' = 'TABLE';
SI Solub;
Transp.'PRECIPITE' = 'TABLE';
Transp.'DISSOLUTION' = 'TABLE';
FINSI;
Transp.'SOUSTYPE' = 'DARCY_TRANSITOIRE';
Transp.'MODELE' = MoTransp;
Transp.'CARACTERISTIQUES' = MaTransp;
Transp.'POROSITE' = MaPor;
Transp.'CONVECTION' = QFHydro;
Transp.'COEF_RETARD' = MaRet1;
SI SoluL;
Transp.'LIMITE_SOLUBILITE' = MaSol1;
FINSI;
SI SoluP;
Transp.'COEF_DISSOLUTION' = MaCoDis1;
FINSI;
* Conditions initiales :
Transp.'TEMPS'. 0 = TIni;
Transp.'TRACE_CONC'. 0 = TC0;
Transp.'CONCENTRATION'. 0 = Conc0;
Transp.'FLUX'. 0 = QFDif0;
SI Solub;
Transp.'PRECIPITE'. 0 = Prec0;
Transp.'DISSOLUTION'. 0 = Libe0;
FINSI;
* Conditions aux limites :
Transp.'BLOCAGE' = RigCL;
Transp.'TRACE_IMPOSE' = ChgCL;
* Paramètres numériques :
Transp.'THETA_DIFF' = TetaDif;
Transp.'THETA_CONVECTION' = TetaCnv;
SI SoluP;
Transp.'THETA_DISS' = TetaDis;
FINSI;
Transp.'TEMPS_CALCULES' = LiTCalc;
Transp.'TEMPS_SAUVES' = LiTSauv;
* ==========
* | CALCUL |
* ==========
DARCYTRA 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.Flux.Id1)
+ (TetaDif * Transp.Flux.Id);
DebD2 = 'REDU' DebD1 CeFlux ;
DebD3 = 'RESULT' (DebD2 * Ori2);
FlD = ('MAXIMUM' DebD3);
Fl = FlD;
* Débit convectif sortant :
THCnv = ((1. - TetaCnv) * Transp.Trace_Conc.Id1)
+ (TetaCnv * Transp.Trace_Conc.Id);
DebC1 = ('EXCO' 'TH' THCnv 'SCAL') * Transp.Convection;
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é du champ de limite de solubilité :
* ----------------------------------------
SI (TrLimSol ET SoluL);
'OPTION' ISOV SURF;
'TITRE' 'Limite de solubilite';
MaSol2 = 'KCHA' MoTransp MaSol1 'CHAM';
'TRAC' MoTransp MaSol2 LTrame;
FINSI;
* Tracé du champ de coefficient de dissolution :
* ----------------------------------------------
SI (TrCoDis ET SoluP);
'TITRE' 'Coefficient de dissolution (1/an)';
MaCoDis2 = 'KCHA' MoTransp MaCoDis1 'CHAM';
'OPTION' ISOV SURF;
'TRAC' MoTransp MaCoDis2 LTrame;
FINSI;
* Tracé du champ de coefficient de retard :
* -----------------------------------------
SI (TrRet ET ('EXISTE' Transp COEF_RETARD));
'TITRE' 'Coefficient de Retard';
MaRet0 = Transp.COEF_RETARD;
SI ('EGA' (TYPE MaRet0) 'FLOTTANT');
MaRet1 = MANU 'CHPO' CeTransp 1 'SCAL' MaRet0;
SINON;
MaRet1 = MaRet0;
FINSI;
MaRet2 = 'KCHA' MoTransp MaRet1 'CHAM';
MaRet3 = 'KCHT' MoTransp 'SCAL' 'CENTRE' ('NOMC' MaRet1 'SCAL');
MaRet4 = 'ELNO' MoTransp MaRet3;
'OPTION' ISOV SURF;
Rmin = 'MINIMUM' MaRet1;
Rmax = 'MAXIMUM' MaRet1;
SI (('MINIMUM' MaRet1) < (('MAXIMUM' MaRet1) - EpsPrec));
Ech = 'PROG' Rmin 'PAS' ((Rmax - Rmin) / 4.) Rmax;
'TRAC' MoTransp MaRet2 LTrame Ech;
SINON;
Ech = 'PROG' 0. ('MAXIMUM' MaRet1);
'TRAC' MoTransp MaRet2 LTrame Ech;
FINSI;
FINSI;
* Tracé des profils de conc, prec, diss le long des lignes choisies :
* -------------------------------------------------------------------
SI TrLEvol;
IndLEvol = 'INDEX' LEvol;
TabEvlC = table;
SI Solub;
TabEvlP = table;
TabEvlL = table;
Finsi;
'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 (SoluL et ('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;
nn = 'DIME' TabEvlC.i ;
'REPETER' bcl2 nn ;
j = nn - &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 (Solub et ('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)' ;
* Evolution des dissolutions :
TitDesL = 'CHAINE' 'Dissolution ' (LEvol.i.'T') ;
TabEvlL.i = TRACHIS Transp 'DISSOLUTION' (LEvol.i.'M') LiIndevl
('MOTS' mot1) ('MOTS' ' ') 'PREF' ' ' ;
DESTRA (TabEvlL.i) 'MIMA' 'AXES' 'TITR' TitDesL
'TITX' (LEvol.i.'L') 'TITY' 'L (mol/dm3/an)' ;
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;
SI Solub;
TabEvpP = table;
TabEvpL = table;
Finsi;
'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 (SoluL et ('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);
nn = 'DIME' TabEvpC.i ;
'REPETER' bcl2 nn ;
j = nn - &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 (Solub et ('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 ;
* Evolution des dissolutions :
TitDesL = 'CHAINE' 'Dissolution ' (PEvol.i.'T') ;
TabEvpL.i = TRACHIT Transp 'DISSOLUTION'
('MANU' 'POI1' (PEvol.i.'M')) LiIndevl
('MOTS' mot1) 'PREF' ' ' ;
DESTRA (TabEvpL.i) 'MIMA' 'AXES' 'TITR' TitDesL
'TITX' 't (an)' 'TITY' 'L (mol/dm3/an)'
'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 ET Solub;
*-- On trace pour chaque temps choisi :
SI (TrConc OU TrPrec);
*-- Echelles :
SI TrConc;
SI SoluL;
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;
SI TrConc;
Conc2 = 'KCHA' MoTransp Transp.Concentration.Id 'CHAM';
MoTot = MoTransp;
Tot2 = Conc2;
FINSI;
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 ET Solub;
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;
Libe2 = Transp.Dissolution.Id2;
Libe3 = 'KCHA' MoTransp Libe2 'CHAM';;
SI ((('MAXIMUM' Libe3) - ('MINIMUM' Libe3)) > EpsPrec);
'OPTION' ISOV SURF;
'TRAC' MoTransp Libe3 Ltrame ;
FINSI;
'FIN' Bcl22;
FINSI;
* Tracé du flux traversant la ligne choisie (en mol/l/an) :
* ---------------------------------------------------------
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' .30290 .58543 .78354 .92082 1.0185 1.0894 1.0285
.86475 .69182 .54609 ;
*'LISTE' lifl ;
err1 = 'MAXIMUM' ('ABS' (lifl - liflsav)) / ('MAXIMUM' liflsav);
mess 'Erreur = ' err1 ;
si (err1 > CRIT1);
'ERREUR' 5;
sinon;
'ERREUR' 0;
finsi;
Fin;
© Cast3M 2003 - All rights reserved.
Disclaimer