* fichier : cacul.dgibi
opti echo 2;
**********************************************************************************
* LA SOLUTION ANALYTIQUE
**********************************************************************************
ALF = 0.08D0 * 100.D0 ;
KS = 5.85D-2 / 3600.D0 ;
TS = 0.3D0 ;
*TS = 1.D0;
TR = 0.055D0 ;
*TR = 0.0D0;
TSR = TS - TR ;
N = 2.0304D0 ; N1 = 1.D0/N ; MOINSN = (-1.D0)*N ;
M = -0.5075D0 ; M1 = 1.D0/M ;
BET = -0.029227D0 * 100.D0 ;
nb = 300 ;
H00 = -1.D-12 ;
Hinf = -1.15D0 ;
RA2PI = PI ** (-0.5D0) * 2.D0 ;
*opti trac ps ;
********************************************
DEBPROC KKR_pro TT*FLOTTANT ;
*K0 = (TT '-' TR) '/' TSR * ALF EXP * KS ;
K0 = (TT '-' TR) '/' TSR ** ALF * KS ;
FINPROC K0 ;
********************************************
********************************************
DEBPROC HHT_pro H1*FLOTTANT ;
T1 = H1 * BET ** N + 1.D0 ** m * TSR + TR ;
FINPROC T1 ;
********************************************
********************************************
DEBPROC TH_pro T1*FLOTTANT ;
H1 = T1 - TR / TSR ** M1 - 1.D0 ** N1 / BET ;
FINPROC H1 ;
********************************************
********************************************
DEBPROC DHDT_pro H1*FLOTTANT T1*FLOTTANT ;
DHDT0 = BET*(H1 '-' 1.D-13) ** MOINSN + 1. * (H1 '-' 1.D-13)
/ (T1 - TR) * N1 * M1 ;
FINPROC DHDT0 ;
********************************************
********************************************
DEBPROC KINPHI int00*FLOTTANT ;
phi0 = 0. ;
int0 = int00 ;
Int_p0 = prog int0 ;
phi_p0 = prog phi0 ;
repeter b1 (nb - 1);
d0 = extr D2_prog &b1 ;
phi0 = 2. * d0 / int0 + phi0 ;
phi_p0 = phi_p0 et (prog phi0) ;
int0 = int0 - phi0 ;
Int_p0 = Int_p0 et (Prog int0) ;
fin b1 ;
d0 = extr D2_prog NB;
phi0 = 2. * d0 / iNT0 + phi0 ;
phi_p0 = phi_p0 et (prog phi0) ;
FINPROC phi_p0 Int_p0 ;
********************************************
********************************************
DEBPROC KPROG T00*FLOTTANT TN0*FLOTTANT NB0*ENTIER MO1*MOT ;
DT = T00 - TN0 / NB0 ;
SI (EGA MO1 DEMI ) ;
MESS 'DEMI = ' MO1 ;
T = T00 - (dt*0.5) ;
h_pr0 = prog ;
K_pr0 = prog ;
D_pr0 = prog ;
C_pr0 = prog ;
T_pr0 = prog ;
SINON ;
MESS 'TOTA = ' MO1 ;
T = T00 ;
H = TH_pro T ;
* K = KR_pro H ;
K = KKR_pro T ;
DHDT = DHDT_pro H T ;
D = K*DHDT ;
C = 1.D0/(DHDT '+' 1.D-10) ;
h_pr0 = prog H;
K_pr0 = prog K;
D_pr0 = prog D;
C_pr0 = prog C;
T_pr0 = prog T;
T = T - DT ;
FINSI;
repe b0 nb0 ;
H = TH_pro T ;
* K = KR_pro H ;
K = KKR_pro T ;
DHDT = DHDT_pro H T ;
D = K*DHDT ;
C = 1./DHDT ;
h_pr0 = h_pr0 ET (prog H) ;
K_pr0 = K_pr0 ET (prog K) ;
D_pr0 = D_pr0 ET (prog D) ;
C_pr0 = C_pr0 ET (prog C) ;
T_pr0 = T_pr0 ET (prog T);
T = T - DT ;
FIN B0 ;
FINPROC h_pr0 K_pr0 D_pr0 C_pr0 T_pr0 ;
********************************************
********************************************
DEBPROC A x*FLOTTANT ;
si (x < 2.9 ) ;
X2 = (-1.D0)*X*X ;
A0 = X2 EXP * X * (1. - (ERF X) ** (-1.D0)) * RA2PI + (2.D0*X2) ;
sinon ;
X2 = X ** (-2.D0) * -0.5D0 ;
A0 = 108830.*X2+8162.*X2+706.*X2+74.*X2+10.*X2+2.*X2+1.;
*A0 = 8162.*X2+706.*X2+74.*X2+10.*X2+2.*X2+1.;
*A0 = 706.*X2+74.*X2+10.*X2+2.*X2+1.;
*A0 = 74.*X2+10.*X2+2.*X2+1.;
FINSI ;
FINPROC A0 ;
********************************************
********************************************
DEBPROC INTN phi*LISTREEL ;
PNM1 = EXTR phi NB ;
DD0 = EXTR D2_PROG NB ;
XX = PNM1 * 0.5D0 * (DD0 ** -0.5D0) ;
IC0 = PNM1*0.5D0 + (2.D0*DD0/PNM1*(A XX)) ;
FINPROC IC0 ;
********************************************
HN = HINF ;
TN = HHT_pro HN ; LIST TN ;
HN = TH_pro TN ; LIST HN ;
*KN = KR_pro HN ; LIST KN ;
KN = KKR_pro TN ; LIST KN ;
DHDTN = DHDT_pro HN TN ; LIST DHDTN ;
DN = KN*DHDTN ;
CN = 1./DHDTN ;
T0 = HHT_pro H00 ;
LIST (HN - (TH_pro (TN - 1.E-6)) / 1.E-6) ;
DT = T0 - TN / NB ;
h_prog K_prog D_prog C_prog T_prog = KPROG T0 TN NB 'TOTA';
TfH_evol = evol manu h_prog T_prog ;
*DESS TfH_evol MIMA
TITR 'teneur en eau fct(P)'
TITX 'pression (m)'
TITY 'T eau' ;
KfH_evol = evol manu h_prog K_prog ;
*DESS KfH_evol MIMA LOGY
TITR 'perméabilité relative fct(P)'
TITX 'pression (m)'
TITY 'KR (m/s)' ;
CfH_evol = evol manu h_prog C_prog ;
*DESS CfH_evol MIMA
TITR 'capacité fct(P)'
TITX 'pression (m)'
TITY 'C m' ;
DfH_evol = evol manu h_prog D_prog ;
*DESS DfH_evol MIMA LOGY
TITR 'diffusivité fct(P)'
TITX 'pression (m)'
TITY 'D m2/s' ;
TfH_evol = evol manu T_prog H_prog ;
*DESS TfH_evol MIMA
TITR 'teneur en eau fct(T)'
TITX 'teneur en eau'
TITY 'T eau' ;
KfH_evol = evol manu T_prog K_prog ;
*DESS KfH_evol MIMA LOGY
TITR 'perméabilité relative fct(T)'
TITX 'teneur en eau'
TITY 'KR (m/s)' ;
CfH_evol = evol manu T_prog C_prog ;
*DESS CfH_evol MIMA
TITR 'capacité fct(T)'
TITX 'teneur en eau'
TITY 'C m' ;
*DfH_evol = evol manu T_prog D_prog ;
*DESS DfH_evol MIMA LOGY
*TITR 'diffusivité fct(T)'
*TITX 'teneur en eau'
*TITY 'D m2/s' ;
*OPTI TRAC X ;
TN = HHT_pro Hinf ;
TN1 = extr T_prog (DIME T_prog) ;
T01 = extr T_prog 1 ;
db0 = 2.D0 / ((T01 - TN1) ** 2.D0) ;
x_prog = T_prog - (prog (DIME T_prog)*TN1) * D_prog ;
DfT_evol = evol manu T_prog x_prog ;
*-------------------------- Dmoy par D
Db = (-1.D0) * ((somm DfT_evol) extr 1) * db0 ;
*-------------------------- calcul de I 1/2 (j = 1)
*????????????????????????????????????
Int = Db/pi**0.5D0*2.D0*nb ; list Int ;
Y = T_prog - (prog (DIME T_prog)*TN1) * K_prog ;
KfH_evol = evol manu H_prog Y ;
*-------------------------- Dmoy par K
Db = (-1.D0) * ((somm KfH_evol) extr 1) * db0 ;
*-------------------------- calcul de I 1/2 (j = 1)
*????????????????????????????????????
Int = Db/pi**0.5*2.*nb ; list Int ;
*Int = 1.e-3/(ts - tn)*nb ;list Int ;
*Int = int * 2. ;
*-------------------------- calcul de D (I + 1/2)
h2_prog K2_prog D2_prog C2_prog T2_prog = KPROG T0 TN nb DEMI ;
*-------------------------- calcul de PHI (I + 1/2)
phi_prog Int_prog = KINPHI Int ;
I2_evol = evol manu T_prog phi_prog ;
I2mi = I2_evol somm * (-1.) / dt ; LIST I2mi ;
*I2_evol = evol manu phi_prog T_prog ;
*DESS I2_evol ;
*opti donn 5 ;
*-------------------------- calcul de ICHA (N - 1/2)
ICHA = INTN phi_prog ;
*-------------------------- DELTA (J=1)
DELTA = (EXTR Int_prog NB) - ICHA ;
*-------------------------- calcul de I 1/2 (j = 2)
int1 = int ;
Int = int1 - (DELTA * 0.5d0) ;
deltI = Int - Int1 ;
repeter b2 20 ;
*-------------------------- calcul de PHI (I + 1/2)
phi_prog Int_prog = KINPHI Int ;
I2_evol = evol manu T_prog phi_prog ;
I2mi = I2_evol somm * (-1.D0) / dt ;
*-------------------------- calcul de ICHA (N - 1/2)
ICHA = INTN phi_prog ;
mess
'iter :' &b2 'icha' ICHA 'DELTA' DELTA 'Itot' (maxi I2mi) 'int' int;
*-------------------------- DELTA (J=1)
delta1 = DELTA ;
DELTA = (EXTR Int_prog NB) - ICHA ;
*-------------------------- calcul de I 1/2
int1 = int ;
*int = int1 - (DELTA * delta1 / (DELTA + delta1)) ;
*int = int1 - (DELTA * deltI / (DELTA + delta1)) ;
Int = int1 - (DELTA * 0.5d0 ) ;
deltI = Int - Int1 ;
si ( (DELTA abs) < (ICHA*1.e-6))
quitter b2 ;
finsi;
menage ;
fin b2 ;
*'LISTE' h_prog;
*dess ( evol manu phi_prog T_prog);
*dess ( evol manu phi_prog h_prog);
nlist = 'DIME' h_prog;
nnlist = 2 * nlist;
xabs = 'PROG' 0.D0 PAS (1.D0 '/' nnlist) 1.D0;
* on inverse l'axe des x
*'ORDONNER' CROISSANT h_prog ;
*'LISTE' h_prog;
* il faut tronquer la derniere valeur de phi_prog et mettre
* un grand nombre à la place.
roro = 'ENLEVER' phi_prog nlist ;
roro = roro 'ET' ('PROG' 10.D0);
tabcool = 'TABLE' ;
tabcool . 1 = 'TURQ';
tabcool . 2 = 'TURQ';
tabcool . 3 = 'TURQ';
tabcool . 4 = 'TURQ';
***********************************************************************************
* LE CAS TEST
**********************************************************************************
******************** CAS TEST : cacul.dgibi ************************
*
*
* Test de fonctionnement de DARCYSAT en 2D sans effet de gravité.
* ===================================================================
* Infiltration d'eau dans une colonne horizontale de sable uniformément
* désaturé.
*
* - condition initiale : désaturation uniforme correspondant à une
* succion de 1,15 m ;
* - a l'instant initial, une extrémité de la colonne est noyée.
* La frontière est supposée rester a pression nulle ;
*
* ===================================================================
* Les options de modélisation declarées dans la table transmise à
* la procédure DARCYSAT sont les suivantes :
*
*
* - les effets gravitationnels ne sont pas pris en compte : pression =
* charge (indice GRAVITE absent ; valeur par défaut : FAUX);
*
* - Une liste de temps de sauvegarde est fournie en valeur de l'indice
* TEMPS_SAUVES ;
*
* - Le pas de temps est d'abord automatique (indice TEMPS_CALCULES absent)
* L'utilisateur fournit
* > le pas de temps initial (indice 'DT_INITIAL'),
* > le nombre d'itérations recherché par pas de temps (indice 'NITER')
* > le nombre de pas de temps (indice 'NPAS')
*
* - Le calcul est ensuite refais avec une liste de temps de calcul
* donnée à l'indice TEMPS_CALCULES ;
*
* - L'homogenéisation spatiale des propriétes physiques s'effectue
* par moyenne arithmétique des valeurs obtenues aux faces (indice
* HOMOGENEISATION de valeur DECENTRE)
*
* cas test tiré du rapport DMT 97/25 :
* "Implémentation dans CASTEM 2000 d'un modèle de transfert hydrique
* en milieu poreux non saturé"
*
*
********************* CAS TEST : cacul.dgibi ************************
'OPTION' 'ECHO' 1 ;
'SAUTER' 'PAGE';
*
'TITRE' 'infiltration horizontale dans le sable : cacul.dgibi' ;
'OPTION' 'DIME' 2 'ELEM' 'QUA4' ;
*'OPTION' 'ISOV' 'LIGN' ;
*'TRACER' 'PSC' ;
************************************************************************
* Fonction qui calcule la perméabilité en fonction de la
* saturation réduite.
*
* Loi puissance :
* Perméabilité : K = Ks S^B
*
* Paramètres physiques à définir dans la table LOI :
* --------------------------------------------------
* ALPHA : coef. B (s.d.)
* PERMSAT : coef. Ks, perméabilité à saturation (m/s)
*
* Entrée :
* --------
* LOI : table de données décrite ci-dessus
* SAT : saturation réduite
*
* Sortie :
* --------
* K1 : perméabilité totale en eau (m/s)
*
************************************************************************
************************************************************************
COMPLET = FAUX ;
GRAPH = FAUX ;
*--------------------------------------------------------------------
* Création maillage
*
*- Discrétisation :
ENX = 1 ;
ENY = 200 ;
'DENS' (1./eny) ;
*- Création des points et des droites
*
A0 = -0.5 1.0D0; B0 = 0.5 1.0D0;
dens (4./eny);
A1 = -0.5 0.75D0; B1 = 0.5 0.75D0;
A2 = -0.5 0.5D0 ; B2 = 0.5 0.5D0 ;
A3 = -0.5 0.0D0; B3 = 0.5 0.0D0 ;
*- Création des droites
*
AB0 = 'DROIT' ENX A0 B0 ;
AB1 = 'DROIT' ENX A1 B1 ;
AB2 = 'DROIT' ENX A2 B2 ;
AB3 = 'DROIT' ENX A3 B3 ;
*- Creation des surfaces
*
MASSIF0 = AB3 'REGLER' AB2 'REGLER' AB1;
MASSIF0 = MASSIF0 'REGLER' AB0 ;
ENY = 'NBEL' MASSIF0 ;
ENTR = 'COULEUR' ('INVERSE' AB0) 'ROUGE' ;
SORT = 'COULEUR' (AB3) 'ROUGE' ;
'SI' GRAPH ;
'TRACER' (MASSIF0 'ET' ENTR 'ET' SORT) ;
'FINSI' ;
*- Creation des maillages contenant tous les points
*
QFTOT = 'CHANGER' MASSIF0 'QUAF' ;
QFSORT = 'CHANGER' SORT 'QUAF' ;
QFENTR = 'CHANGER' ENTR 'QUAF' ;
'ELIMINATION' 0.00001 (QFTOT 'ET' QFSORT 'ET' QFENTR) ;
*- Modèles
*
MODHYB = 'MODELE' QFTOT 'DARCY' 'ISOTROPE' ;
MODENTR = 'MODELE' QFENTR 'DARCY' 'ISOTROPE' ;
MODSORT = 'MODELE' QFSORT 'DARCY' 'ISOTROPE' ;
CEENTR = 'DOMA' MODENTR 'CENTRE' ;
CESORT = 'DOMA' MODSORT 'CENTRE' ;
HYCEN = 'DOMA' MODHYB 'CENTRE' ;
HYFAC = 'DOMA' MODHYB 'FACE';
*- Création ligne de suivi pour le post-traitement et le test
* ligne des points centres (cas 1D)
*
'REPETER' BCL (ENY - 1) ;
IP = &BCL ;
JP = IP + 1 ;
PI = 'POINT' HYCEN IP ;
PJ = 'POINT' HYCEN JP ;
'SI' (IP 'EGA' 1);
LCENC = ('MANU' 'SEG2' PI PJ) ;
'SINON' ;
LCENC = LCENC 'ET' ('MANU' 'SEG2' PI PJ) ;
'FINSI' ;
FIN BCL ;
*--------------------------------------------------------------------
*- pression initiale (metre d'eau) dans le sable
HN = -1.15 ;
*--------------------------------------------------------------------
*- Conditions aux limites
* GBM modifie supprime bloque
*- frontière en limite du domaine de calcul (milieu désaturé)
ESORT = 'MANU' 'CHPO' CESORT 1 'TH' HN 'NATURE' 'DISCRET';
*- frontière mouillée
EENTR = 'MANU' 'CHPO' CEENTR 1 'TH' 0. 'NATURE' 'DISCRET';
*- chargement des CLs
LICALC = 'PROG' 0.D0 1.e20 ;
LIST1 = 'PROG' 2 * 1.D0 ;
VALI0 = 'CHAR' 'THIM' (ESORT et EENTR) ('EVOL' 'MANU' LICALC LIST1) ;
*--------------------------------------------------------------------
*- initialisation des inconnues
*
*- trace de charge d'eau
*TH0 = 'MANU' 'CHPO' ('DIFF' CEENTR HYFAC) 1 'TH' HN 'NATURE' 'DISCRET';
*TH0 = TH0 + SIMPE ;
*- charge d'eau
H0 = 'MANU' 'CHPO' HYCEN 1 'H' HN 'NATURE' 'DIFFUS';
*- flux
QFACE0 = 'MANU' 'CHPO' HYFAC 1 'FLUX' 0.D0 'NATURE' 'DISCRET' ;
* ---------------------------
* = Table DARCY_TRANSITOIRE =
* ---------------------------
*- initialisation table
SATUR = 'TABLE' ;
SATUR. 'TEMPS' = 'TABLE' ;
SATUR. 'CHARGE' = 'TABLE' ;
SATUR. 'FLUX' = 'TABLE' ;
SATUR . 'ITMAX' = 45;
*- données géommétriques
SATUR. 'SOUSTYPE' = 'DARCY_TRANSATUR' ;
SATUR. 'MODELE' = MODHYB ;
SATUR . 'XI' = 1.D-10;
SATUR . 'COEF_EMMAGASINEMENT' = 0.D0 ;
*- instant initial
SATUR. 'TEMPS' . 0 = 0. ;
SATUR. 'CHARGE' . 0 = H0 ;
SATUR. 'FLUX' . 0 = QFACE0 ;
*********************************************
*- conditions aux limites et chargements
SATUR . 'TRACE_IMPOSE' = VALI0 ;
SATUR . 'LUMP' = FAUX;
SATUR . 'TYPDISCRETISATION' = 'EFMH' ;
TABRES = table METHINV;
TABRES . 'TYPINV' = 2;
TABRES . 'PRECOND' = 5;
SATUR . 'METHINV' = TABRES;
SATUR . 'METHINV' . RESID = 1.D-20;
*******************************************
*******************************************
*- données physiques
* loi de perméabilité
LoiP = 'TABLE' 'PUISSANCE';
LoiP. 'ALPHA' = 8.D0 ;
LoiP. 'PERMSAT' = 5.85E-2 / 3600. ;
SATUR.'LOI_PERMEABILITE' = 'TABLE' LoiP ;
* loi de succion
LoiS = 'TABLE' 'VAN_GENUCHTEN';
LoiS. 'PORO' = 0.3 ;
*LoiS. 'PORO' = 1. ;
LoiS. 'TERESIDU' = 0.055 ;
*LoiS. 'TERESIDU' = 0.0 ;
LoiS. 'NEXP' = 2.0304 ;
*LoiS. 'NEXP' = 1.D0 ;
LoiS. 'MEXP' = 0.5075 ;
*LoiS. 'MEXP' = -1.D0 ;
LoiS. 'BHETA' = 0.029227 * 100. ;
SATUR.'LOI_SATURATION' = 'TABLE' LoiS ;
*- données numériques
'SI' (COMPLET) ;
SATUR. 'TEMPS_FINAL' = 10000.D0 ;
'REPETER' bbb 4 ;
tfin = 2500.D0 * &bbb ;
* var de Boltzman
znew = xabs '/' (tfin '**' 0.5D0);
ress = 'IPOL' znew roro h_prog;
'ORDONNER' CROISSANT ress;
'SI' ('EGA' &bbb 1) ;
courb = evol tabcool . 1 manu xabs ress ;
'SINON' ;
courb = courb 'ET'
('EVOL' tabcool . &bbb manu xabs ress);
'FINSI' ;
'FIN' bbb;
TAB1 = 'TABLE' ;
TAB1 . 'TITRE' = 'TABLE' ;
TAB1 . 'TITRE' . 1 = 'exacte 2500' ;
TAB1 . 'TITRE' . 2 = 'exacte 5000' ;
TAB1 . 'TITRE' . 3 = 'exacte 7500' ;
TAB1 . 'TITRE' . 4 = 'exacte 10000' ;
TAB1 . 1 = 'TIRR ';
TAB1 . 2 = 'TIRR ';
TAB1 . 3 = 'TIRR ';
TAB1 . 4 = 'TIRR ';
'SINON' ;
'REPETER' bbb 4 ;
tfin = 500.D0 * &bbb ;
* var de Boltzman
znew = xabs '/' (tfin '**' 0.5D0);
ress = 'IPOL' znew roro h_prog;
'ORDONNER' CROISSANT ress;
'SI' ('EGA' &bbb 1) ;
courb = evol tabcool . 1 manu xabs ress ;
'SINON' ;
courb = courb 'ET'
('EVOL' tabcool . &bbb manu xabs ress);
'FINSI' ;
'FIN' bbb;
TAB1 = 'TABLE' ;
TAB1 . 'TITRE' = 'TABLE' ;
TAB1 . 'TITRE' . 1 = 'exacte 500' ;
TAB1 . 'TITRE' . 2 = 'exacte 1000' ;
TAB1 . 'TITRE' . 3 = 'exacte 1500' ;
TAB1 . 'TITRE' . 4 = 'exacte 2000' ;
TAB1 . 1 = 'TIRR ';
TAB1 . 2 = 'TIRR ';
TAB1 . 3 = 'TIRR ';
TAB1 . 4 = 'TIRR ';
SATUR. 'TEMPS_FINAL' = 2000.D0 ;
'FINSI' ;
SATUR. 'HOMOGENEISATION' = 'CHAINE' 'TRACE' ;
SATUR. 'NPAS' = 4000 ;
SATUR. 'RESIDU_MAX' = 1.D-4 ;
SATUR. 'NITER' = 10 ;
SATUR. 'DT_INITIAL' = 1.D0 ;
'SI' COMPLET ;
SATUR. 'TEMPS_SAUVES' = ('PROG' 250 'PAS' 250 1000)*10. ;
'SINON' ;
SATUR. 'TEMPS_SAUVES' = ('PROG' 500 'PAS' 500 2000) ;
'FINSI' ;
*SATUR. 'TEMPS_CALCULES' = ('PROG' 1 'PAS' 1 10.)*360. ;
SATUR. 'MESSAGES' = 2 ;
'SI' GRAPH ;
*--------------------------------------------------------------------
*- Vérification du choix du dp minimum pour le calcul de la capacité.
* droite support des variables
dx = 'DROIT' (0. -1.15 ) 1000 (0. 0.) ;
zc = 'COOR' dx 2 ;
ev2 = 'EVOL' 'BLEU' 'CHPO' zc 'SCAL' dx ;
px = 'EXTR' ev2 'ORDO' ;
* calcul de la teneur en eau pour la pression zc
s0 t0 cap = HT_PRO (SATUR.'LOI_SATURATION') ZC ;
ev0 = 'EVOL' 'TURQ' 'CHPO' s0 'SCAL' dx ;
evt = 'EVOL' 'VERT' 'MANU' px (100. * ('EXTR' ev0 'ORDO')) ;
'DESSIN' evt 'TITX' 'Pc(m)' 'TITY' 'S(%)'
'TITRE' 'Loi capillaire S(Pc)' ;
* calcul de la teneur en eau pour la pression zc - dp
dp = 1.e-4 ;
s1 t1 cap = HT_PRO (SATUR.'LOI_SATURATION')
('KOPS' zc '-' dp) ;
* représentation de la capacité
c1 = (t0 - t1) / dp;
ev1 = 'EVOL' 'ROUGE' 'CHPO' c1 'SCAL' dx ;
evc = 'EVOL' 'JAUNE' 'MANU' px ('EXTR' ev1 'ORDO') ;
'DESSIN' evc 'TITX' 'Pc(m)' 'TITY' 'Capa(1/m)'
'TITRE' 'Capacite capillaire' ;
'FINSI' ;
* ===========================
* | 1er CALCUL |
* ===========================
*- fonctionnement avec temps automatiques
DARCYSAT SATUR ;
* Post-traitement
* ===============
LT = 'LECT' 1 PAS 1 ((dime SATUR.TEMPS) - 1) ;
liopt = 'MOTS' 'MIMA' 'AXES';
TDES = TRACHIS SATUR 'CHARGE' LCENC LT 'PREF' ' ' 'UNIT' 's' ;
TDES . 5 = 'TABLE' ;
TDES . 5 . VALEUR = 'EXTRAIRE' courb 'COUR' 1;
TDES . 5 . LEGEND2 = TAB1 . TITRE . 1;
TDES . 5 . LEGEND1 = ' ' ;
TDES . 6 = 'TABLE' ;
TDES . 6 . VALEUR = 'EXTRAIRE' courb 'COUR' 2;
TDES . 6 . LEGEND2 = TAB1 . TITRE . 2;
TDES . 6 . LEGEND1 = ' ' ;
TDES . 7 = 'TABLE' ;
TDES . 7 . VALEUR = 'EXTRAIRE' courb 'COUR' 3;
TDES . 7 . LEGEND2 = TAB1 . TITRE . 3;
TDES . 7 . LEGEND1 = ' ' ;
TDES . 8 = 'TABLE' ;
TDES . 8 . VALEUR = 'EXTRAIRE' courb 'COUR' 4;
TDES . 8 . LEGEND2 = TAB1 . TITRE . 4;
TDES . 8 . LEGEND1 = ' ' ;
'SI' GRAPH ;
DESTRA TDES liopt 'TITX' 'z (m)' 'TITY' 'Pw (m)' ;
'FINSI' ;
* On compare l'integrale de la saturation au temps final entre les
* solutions analytiques et non analytiques.
toto = 'PRIM' TDES . 8 . VALEUR;
inta = 'EXTRAIRE' toto 'ORDO';
inta = 'EXTRAIRE' ('DIME' inta) inta;
toto = 'PRIM' TDES . 4 . VALEUR;
intb = 'EXTRAIRE' toto 'ORDO';
intb = 'EXTRAIRE' ('DIME' intb) intb;
err = 'ABS' (inta '-' intb '/' inta);
'LISTE' err;
'SI' (err > 7.D-2) ;
'ERREUR' 5;
'FINSI' ;
'FIN';
© Cast3M 2003 - All rights reserved.
Disclaimer