* fichier : srivastava1VF.dgibi
********************************************************************
GRAPH = faux ;
******************** CAS TEST : srivastava.dgibi **********************
*
* Test de fonctionnement de DARCYSAT en 1D avec effet de gravité
* en regime transitoire.
*
*-------------------------------------------------------------------
*
* A. GENTY - DM2S/SFME/MTMS - 10/2007
*
* ==================================================================
* Infiltration d'eau a debit impose depuis la surface dans un
* milieu 1D non sature limite par une surface inferieure a pression
* d'eau (h0<=0) imposee.
*
* qB imposé
* z=L ---------------------------
* I I
* I I
* I I
* I I
* I hi = solution du regime I
* I permanent avec qA I
* I impose I
* I (charge Hi = hi + z) I
* flux nul I I flux nul
* I I
* I I
* I I
* I I
* I I
* I I
* z=0 ---------------------------
* x=0 h = h0 (charge H = h0) x=x0
*
*
* La solution analytique (en pression) du probleme est decrite
* dans Srivastava et Yeh, 1991, WRR, Vol 27, N 5, pp 753-762.
* ==================================================================
*
* - conditions initiales :
* > la pression initiale hi correspond a la solution du probleme
* en regime permanent avec un debit impose qA.
*
* - conditions aux limites :
* > la pression sur la frontiere inferieure du domaine est imposee
* a h = h0 (charge H = h + z = h0)
* > le flux sur la frontiere superieure du domaine est imposee
* a qB
* > les autres faces exterieures du domaine sont imposees a flux nul
*
* ===================================================================
* Les options de modélisation declarées dans la table transmise à
* la procédure DARCYSAT sont les suivantes :
*
* - les effets gravitationnels sont pris en compte :
* FORCE_GRAVITE = 1.
*
* - Le calcul est réalisé avec l'option de calcul de pas de temps
* fixe dt = 0.5 h sur une duree totale 10 h (20 pas de temps).
*
* - Option numerique testee :
* VF CENTRE TYPINV = 3 PRECOND = 5
*
* - La précision de convergence du Picard demandée est de 1e-03 m
* ===================================================================
* RESULTATS ATTENDUS
*
* Le champs de pression sera compare a la solution analytique a
* 10 heures par l intermediare des erreurs L_infini et L_2.
* Une erreur inferieure a 5 10-2 est attendue.
* ===================================================================
*
'OPTION' 'ECHO' 1 ;
'SAUTER' 'PAGE';
*
'TITRE' 'Transitoire homogene 1D : Srivastava VF' ;
OPTI DIME 2 ELEM QUA4 ;
OPTI ISOV SULI ;
*
*================================================
* Parametres du calcul
*---------------------
* longueur milieu (m)
LONG1 = 1.D0 ;
* largeur du milieu (m) (libre car pheno 1D)
LARG1 = 1.D0 ;
* alpha (m-1) pour la loi exponentielle
ALPHA1 = 10.D0 ;
* permeabilite a saturation (m/s)
KSAT1 = 0.01D0 / 3600.D0 ;
* porosite (-)
PORO1 = 0.4D0 ;
* teneur en eau residuelle (-)
TRESID1 = 0.06D0 ;
* pression imposee en bas (m)
PSI0 = 0.D0 ;
* flux initial qA (m/s)
QA1 = 0.001D0 / 3600.D0 ;
* nouveau flux qB (m/s)
QB1 = 0.009D0 / 3600.D0 ;
*================================================
*-------------------------------------------------------------------
*--------------------- Création du maillage 2D ---------------------
*
*- Densités de mailles
NIV1 = 1 ;
NDX1 = 1 ;
NDY1 = NIV1 * 100 ;
*
*- Cotes du domaine rectangulaire (x0 et z0)
XCA1 = LARG1 ;
YCA1 = LONG1 ;
*
*- Création des points du domaine
A1 = 0. 0. ;
A2 = XCA1 0. ;
A3 = XCA1 YCA1 ;
A4 = 0. YCA1 ;
*
*- Création des lignes
LIN1 = DROIT NDX1 A1 A2 ;
LIN2 = DROIT NDY1 A2 A3 ;
LIN3 = DROIT NDX1 A3 A4 ;
LIN4 = DROIT NDY1 A4 A1 ;
*
*- Faces droite et gauche
LDROIT = LIN2 ;
LGAUCH = LIN4 ;
*
*- Création du domaine
MASSIF0 = DALLER LIN1 LIN2 LIN3 LIN4 PLAN ;
LHAUT = LIN3 COUL BLEU ;
LBAS = LIN1 COUL ROUG ;
*
*-Tracé du maillage avec faces sup et inf pour conditions limites
* et surface d'entree du debit
ELIM 0.00001 (MASSIF0 ET LBAS ET LHAUT ET LGAUCH ET LDROIT) ;
'SI' GRAPH ;
TRAC CACH (MASSIF0 ET LBAS ET LHAUT);
'FINSI' ;
*
***************************************************************************
*
*- Création des maillages contenant tous les points
QFTOT = 'CHANGER' MASSIF0 'QUAF' ;
QFSORB = 'CHANGER' LBAS 'QUAF' ;
QFSORH = 'CHANGER' LHAUT 'QUAF' ;
*
'ELIMINATION' 0.00001 (QFTOT 'ET' QFSORB 'ET' QFSORH) ;
*
*-------------------------------------------------------------------
*----------------------- Création du Modèle ------------------------
*
MODHYB = 'MODELE' QFTOT 'DARCY' 'ISOTROPE' ;
MODSORB = 'MODELE' QFSORB 'DARCY' 'ISOTROPE' ;
MODSORH = 'MODELE' QFSORH 'DARCY' 'ISOTROPE' ;
CESORB = 'DOMA' MODSORB 'CENTRE' ;
CESORH = 'DOMA' MODSORH 'CENTRE' ;
HYCEN = 'DOMA' MODHYB 'CENTRE' ;
HYFAC = 'DOMA' MODHYB 'FACE' ;
HYSOM = 'DOMA' MODHYB 'SOMMET' ;
*
ZSS = COOR 2 HYSOM ;
ZCC = COOR 2 HYCEN ;
ZFF = COOR 2 HYFAC ;
*-------------------------------------------------------------------
*--------------- Conditions aux limites (blocages) ------------------
*
LICALC = 'PROG' 0.D0 1.D8 ;
LIST1 = 'PROG' 2 * 1.D0 ;
*
*- surface inferieure
*--------------------
* Pression imposee P = HN0
HN0 = PSI0 ;
* On impose la charge H = P + z (z=0 en bas)
ESORB = 'MANU' 'CHPO' CESORB 1 'TH' HN0 'NATURE' 'DISCRET';
VALI0 = 'CHAR' 'THIM' (ESORB) ('EVOL' 'MANU' LICALC LIST1) ;
*
*- face superieure
*-----------------
*- Flux impose entrant
Q1 = -1.0D0 * QB1 ;
QCB1 = (Q1 * LARG1) / NDX1 ;
ESORH = 'MANU' 'CHPO' CESORH 1 'FLUX' QCB1 'NATURE' 'DISCRET';
VALI1 = 'CHAR' 'FLUX' (ESORH) ('EVOL' 'MANU' LICALC LIST1) ;
*
*--------------------------------------------------------------------
*----------------- initialisation des inconnues ---------------------
*
*- charge initiale (en metres d'eau) dans le domaine (H = P + z)
* Pression h
QAD1 = QA1 / KSAT1 ;
ZAD1 = ZCC * ALPHA1 ;
*
CST1 = EXP(ALPHA1 * HN0) - QAD1 ;
*
EAP1 = (EXP(-1. * ZAD1)) * CST1 + QAD1 ;
PM = (LOG(EAP1)) / ALPHA1 ;
*
* Traces de controle condition initiale en pression
* valeur aux noeuds
ZAD2 = ZSS * ALPHA1 ;
EAP2 = (EXP (-1. * ZAD2)) * CST1 + QAD1 ;
PM2 = LOG(EAP2) * (1. / ALPHA1) ;
*
CHP2 = ELNO MODHYB PM ;
'SI' GRAPH ;
TRAC CHP2 MASSIF0 ;
'FINSI' ;
CHP1 = KCHA MODHYB PM 'CHAM' ;
'SI' GRAPH ;
TRAC MODHYB CHP1 ;
'FINSI' ;
*
PIN1 = REDU CHP2 LIN2 ;
EVIN0 = EVOL BLEU 'CHPO' PIN1 LIN2 ;
'SI' GRAPH ;
DESS EVIN0 ;
'FINSI' ;
*
PINI1 = REDU PM2 LIN2 ;
EVIN1 = EVOL VERT 'CHPO' PINI1 LIN2 ;
'SI' GRAPH ;
DESS (EVIN1 ET EVIN0) ;
'FINSI' ;
*
* Charge H
HM = PM '+' ZCC;
*
H0 = nomc (HM) 'H' 'NATURE' 'DISCRET';
*
*- Flux
QFACE0 = MANU CHPO HYFAC 1 'FLUX' 0.D0 'NATURE' 'DISCRET' ;
*
*--------------------------------------------------------------------
*----------------------------- Calcul -------------------------------
*
* ---------------------------
* = Table DARCY_TRANSITOIRE =
* ---------------------------
*- initialisation table
SATUR = 'TABLE' ;
SATUR. 'TEMPS' = 'TABLE' ;
SATUR. 'TRACE_CHARGE' = 'TABLE' ;
SATUR. 'CHARGE' = 'TABLE' ;
SATUR. 'FLUX' = 'TABLE' ;
*
*- données géométriques
SATUR. 'SOUSTYPE' = 'DARCY_TRANSATUR' ;
SATUR. 'MODELE' = MODHYB ;
*
*- instant initial
SATUR. 'TEMPS' . 0 = 0. ;
SATUR. 'CHARGE' . 0 = H0 ;
SATUR. 'FLUX' . 0 = QFACE0 ;
*
*- flux impose
SATUR. 'FLUX_IMPOSE' = VALI1 ;
***********************************************************
* SEULS GROS CHGTS GBM
***********************************************************
*
*- conditions aux limites et chargements
* GBM - MODIFIE
SATUR . 'TRACE_IMPOSE' = VALI0 ;
*
* GBM - NOUVEAU
SATUR . 'LUMP' = FAUX ;
SATUR . 'TYPDISCRETISATION' = 'VF' ;
SATUR . 'DECENTR' = FAUX ;
*
* gravite
SATUR . 'FORCE_GRAVITE' = 1.D0 ;
*
************************************************************
*- données physiques
************************************************************
*- Loi de permeabilite relative (exponentielle de h)
LoiP = 'TABLE' 'EXPONENTIELLE';
LoiP. 'ALPHA' = ALPHA1 ;
LoiP. 'PERMSAT' = KSAT1 ;
LoiP. 'COEF_N' = 1.0 ;
LoiP. 'COEF_C' = 1.0 ;
SATUR.'LOI_PERMEABILITE' = 'TABLE' LoiP ;
*
*****************************************************************
*- Loi de succion (exponentielle de h)
LOIS = 'TABLE' 'EXPONENTIELLE' ;
LOIS.'PORO' = PORO1 ;
LOIS.'TERESIDU' = TRESID1 ;
LOIS.'COEF_N' = 1.0 ;
LOIS.'COEF_C' = 1.0 ;
LOIS.'BHETA' = ALPHA1 ;
SATUR . 'LOI_SATURATION' = 'TABLE' LOIS ;
*
*****************************************************************
*****************************************************************
*- données numériques
*
* Temps de calcul final 10 heures
TF1 = 10. * 3600. ;
SATUR. 'TEMPS_FINAL' = TF1 ;
SATUR. 'HOMOGENEISATION' = 'CHAINE' 'CENTRE' ;
SATUR. 'RESIDU_MAX' = 1.D-3 ;
SATUR. 'NITER' = 25 ;
SATUR. 'MESSAGES' = 2 ;
*
* SOLVEUR Grad conj, ILU0
TABRES = table METHINV;
TABRES . 'TYPINV' = 3;
TABRES . 'PRECOND' = 5;
*
SATUR . 'METHINV' = TABRES;
SATUR . 'METHINV' . RESID = 1.D-25;
*
SATUR . 'TEMPS_CALCULES' = ('PROG' 0. 'PAS' 1. 20.) * (3600. / 2.) ;
*
*********************************************************************
*------------------------- Exécution procédure ----------------------
*********************************************************************
*
DARCYSAT SATUR ;
*
*********************************************************************
*
*-Legende des traces
TABI1 = TABLE ;
TABI1.'TITRE'= TABLE ;
TABI1.1 = 'REGU MARQ CARR' ;
TABI1.'TITRE' . 1 = MOT '0 h' ;
TABI1.2 = 'REGU MARQ TRIA' ;
TABI1.'TITRE' . 2 = MOT '10 h' ;
TABI1.3 = 'REGU MARQ LOSA' ;
TABI1.'TITRE' . 3 = MOT 'Analytique 10 h' ;
*
*-------------------------------------------------------------------
*-------------------------- Extration des charges ---------------------
NDI1 = (('DIME' SATUR . TEMPS) '-' 1) ;
* Dates de traces
T44 = SATUR . TEMPS . NDI1 ;
CALCTH4 = SATUR . CHARGE . NDI1 ;
*
*-------------------------------------------------------------------
*------------------- Calcul et trace de la pression ----------------
CPRE4 = CALCTH4 '-' ('NOMC' 'H' ZCC) ;
*
CALOO4 = KCHA MODHYB CPRE4 'CHAM';
*
SAVCAL = NOMC SCAL CPRE4 ;
*
* Trace sur elements
'SI' GRAPH ;
TRAC MODHYB CALOO4 ;
'FINSI' ;
*
* Trace sur noeuds
CALNO4 = ELNO MODHYB CPRE4 ;
CONT1 = CONTOUR MASSIF0 ;
*
'SI' GRAPH ;
TRAC CALNO4 (MASSIF0) CONT1 ;
'FINSI' ;
*
* Trace 1D sur les noeuds
VL4D = REDU CALNO4 LIN2 ;
EVOL4 = EVOL VERT 'CHPO' VL4D LIN2 ;
'SI' GRAPH ;
DESS (EVIN1 ET EVOL4) LEGE TABI1 ;
'FINSI' ;
*
*--------------------------------------------------------------------------
* =========================================================================
* Solution analytique Srivastava and Yeh homogene
*--------------------------------------------------------------------------
*
* recuperation des coordonnees Z aux centres
ZCCZ1 = COOR 2 HYCEN ;
*
* constantes utiles
DATE1 = T44 ;
TAD1 = (ALPHA1 * KSAT1 * DATE1) / (PORO1 - TRESID1) ;
KZER1 = LoiP. 'PERMSAT' ;
ALPHA1 = ABS((LoiP. 'ALPHA')) ;
HZER1 = HN0 ;
LNEW1 = ALPHA1 * LONG1 ;
QBD1 = QB1 / KSAT1 ;
*
* Z addimentionnel
ZCC1 = ZCCZ1 * ALPHA1 ;
*
*----------------------------------------------------------------------------
* calcul du premier terme en exp(-z)
* TT1(Z) = QB - (QB - exp(alpha*h0))*exp(-z)
*----------------------------------------------------------------------------
CC1 = (QBD1 - (EXP(ALPHA1*HZER1))) ;
TT1 = QBD1 - (CC1 * (EXP(-1.D0 * ZCC1))) ;
*
*----------------------------------------------------------------------------
* calcul du second terme
* Remarque: l'indice n=0 n'est pas calcule dans la somme
*----------------------------------------------------------------------------
* constantes trigo pour calcul
PI1 = PI ;
LIST PI1 ;
CSD1 = 180.D0/PI1 ;
*
CSTE1 = -4.D0 * (QBD1 - QAD1) * (EXP(-0.25D0 * TAD1));
TLZ1 = 0.5D0 * (LNEW1 - ZCC1) ;
CSTE2 = EXP(TLZ1) ;
*
NBOUC1 = 1000 ;
NNN1 = 1.D0 ;
SOM1 = 0.D0 ;
* calcul de la somme
*---------------------------------------------
REPETER BOU1 NBOUC1 ;
*
* boucle de calcul des L_n et somme en sin(L_n.Z),
* partir de n=1.
* Calcul de la somme.
*
MESS 'Indice de Somme =' NNN1 ;
*******************************
* extraction des L_n avec un Newton
* Calcul de 2L_n = -tan(L_n L) par la methode de Newton-Raphson
* On cherche en fait les zeros de f(al) = al + (L/2)*tan(al)
* avec al = L_n L
*
* Point de depart (AL0) et erreur recherchee (ERCONV1)
EPSI1 = 1.E-4 ;
AL0 = (NNN1 - 0.5D0) * PI1 + EPSI1 ;
ERR1 = 1.D0 ;
ERCONV1 = 1.E-06 ;
REPETER BLOC1 ;
AL0D1 = AL0 * CSD1 ;
FAL1 = AL0 + ((0.5D0 * LNEW1) * ((SIN(AL0D1))/(COS(AL0D1)))) ;
DFAL1 = 1.D0 + ((0.5D0 * LNEW1) / ((COS(AL0D1))*(COS(AL0D1)))) ;
AL1 = AL0 - (FAL1 / DFAL1) ;
AL1D1 = AL1 * CSD1 ;
ERR1 = abs(AL1 + ((0.5D0 * LNEW1)*((SIN(AL1D1))/(COS(AL1D1))))) ;
AL0 = AL1 ;
SI (ERR1 < ERCONV1) ;
QUITTER BLOC1 ;
FINSI ;
FIN BLOC1 ;
Mess 'Iterations Newton =' (&BLOC1 - 1) ;
LM1 = AL0 / LNEW1 ;
*
*******************************
* calcul de SI1 = sin(L_n Z)
MU1Z = KOPS ZCC1 '*' LM1 ;
* en degres
MU2Z = MU1Z * CSD1 ;
SI1 = SIN(MU2Z) ;
*******************************
* calcul de SI2 = sin(L_n L)
LAL1 = LM1 * LNEW1 ;
* en degres
LALD1 = LAL1 * CSD1 ;
SI2 = SIN(LALD1) ;
***************************************
* calcul de exp(-t L_n**2)
LM2 = LM1 * LM1 ;
TEX1 = -1.D0 * TAD1 * LM2 ;
NP1 = EXP(TEX1) ;
***************************************
* calcul du denominateur
DEN1 = 1.D0 + (0.5D0 * LNEW1) + (2.D0 * LM2 * LNEW1) ;
***************************************
* calcul produit total et somme
COEF1 = (SI1 * SI2 * NP1) / DEN1;
SOM1 = SOM1 + COEF1 ;
NNN1 = NNN1 + 1. ;
FIN BOU1 ;
***************************************
* fin de calcul du terme somme
*-----------------------------------------------
*
* Calcul de K
KTOT2 = SOM1 * CSTE1 * CSTE2 ;
KTOT1 = KTOT2 + TT1 ;
*
**************************************************************
*
* calcul de h associe a K
VPHK2 = LOG(KTOT1) ;
VCH1 = VPHK2 / ALPHA1 ;
* Trace h sur elements
HFI1E = KCHA MODHYB VCH1 'CHAM';
'SI' GRAPH ;
TRAC MODHYB HFI1E ;
'FINSI' ;
* Trace h sur noeuds
HFI1N = ELNO MODHYB VCH1 ;
'SI' GRAPH ;
TRAC HFI1N (MASSIF0) CONT1 ;
'FINSI' ;
*
**************************************************************
VL8D = REDU HFI1N LIN2 ;
EVOL8 = EVOL BLEU 'CHPO' VL8D LIN2 ;
'SI' GRAPH ;
DESS (EVIN1 ET EVOL4 ET EVOL8) LEGE TABI1 ;
'FINSI' ;
**************************************************************
*
* Calcul de l'erreur relative par element (en pression)
error = ((SAVCAL '-' VCH1) '/' VCH1) 'ABS' ;
error = 'KCHA' modhyb error CHAM;
'SI' GRAPH ;
'TRACER' modhyb error;
'FINSI' ;
*
* Erreur Linfini
NUME1 = ABS(SAVCAL '-' VCH1) ;
DENU1 = ABS(VCH1) ;
DENU2 = MAXI (DENU1) ;
ERR1 = NUME1 / DENU2 ;
ER1C = KCHA MODHYB ERR1 'CHAM';
'SI' GRAPH ;
TRAC MODHYB ER1C ;
'FINSI' ;
*
ERINF1 = MAXI (ERR1) ;
MESS '________________________' ;
MESS 'ERREUR L INFINI =' ERINF1;
MESS '________________________' ;
*
**********************
* Calcul de l'erreur L2 par element (en pression)
ERRL21 = NUME1 ** 2.D0 ;
ERRL22 = DENU1 ** 2.D0 ;
ERRL2 = (ERRL21 '/' ERRL22) ** 0.5D0 ;
ERRL3 = 'KCHA' MODHYB ERRL2 CHAM ;
'SI' GRAPH ;
TRAC MODHYB ERRL3 ;
'FINSI' ;
*
* Calcul erreur L2
SURF1 = 'DOMA' MODHYB 'VOLUME' ;
ERNUM1 = ERRL21 * SURF1 ;
ERDEN1 = ERRL22 * SURF1 ;
INT1 = RESU ERNUM1 ;
INT2 = RESU ERDEN1 ;
ERL2 = (MAXI (INT1 / INT2)) ** 0.5 ;
MESS '________________________' ;
MESS 'ERREUR L2 =' ERL2 ;
MESS '________________________' ;
*
**************************************************************
* Test de non regression
ERMAX1 = 0.05 ;
SI (ERL2 > ERMAX1) ;
mess 'ERREUR ANORMALE' ERL2 '>' ERMAX1 ;
'ERREUR' 5;
SINON ;
'ERREUR' 0;
FINSI ;
*
*-------------------------------------------------------------------
*
FIN;
© Cast3M 2003 - All rights reserved.
Disclaimer