* fichier : cacul.dgibi ************************************************************************ * Section : Fluides Transitoire ************************************************************************ 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) * ((INTG DfT_evol) ) * 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) * ((INTG KfH_evol) ) * 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 INTG * (-1.) / dt prog ; 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 INTG * (-1.D0) / dt prog ; *-------------------------- 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' 'DECENTRE' ; 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';