Accueil | Cast3M

* fichier : unsat_lindiriEFMH.dgibi
************************************************************************
************************************************************************
********************************************************************
GRAPH = faux ;
******************** CAS TEST : unsat_lindiri.dgibi **********************
*
* Test de fonctionnement de DARCYSAT en 1D avec effet de gravité
* en regime transitoire.
*
*-------------------------------------------------------------------
*
*     A. GENTY - DM2S/SFME/LSET - 03/2009
*
* ==================================================================
* Infiltration d'eau a pression d'eau imposee (h2<=0) depuis la 
* surface dans un milieu 1D non sature limite par une surface 
* inferieure a pression d'eau (h0<=0) imposee.
*
*                    h = h2 (charge H = h2 + L) 
*            z=L    ---------------------------
*                   I                         I
*                   I                         I
*                   I                         I    
*                   I                         I   
*                   I hi = solution du regime I 
*                   I      permanent avec h1  I
*                   I      impose en z=L      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 obtenue
* par linearisation (Cf. Srivastava et Yeh, 1991, WRR, Vol 27, 
* N 5, pp 753-762) puis par separation des variables.
* ==================================================================
*
* - conditions initiales : 
*   > la pression initiale hi correspond a la solution du probleme
*     en regime permanent avec h=h1 impose en z=L.
*
* - conditions aux limites :
*   > la pression sur la frontiere inferieure du domaine est imposee
*     a h = h0 (charge H = h + z = h0)
*   > la pression sur la frontiere superieure du domaine est imposee
*     a h = h2 (charge H = h2 + L)
*   > 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.1 h sur une duree totale 8 h.
*
* -  Les options numeriques testees sont les suivantes:
*       EFMH   DECENTRE    TYPINV = 3  PRECOND = 5
*
* -  La précision de convergence demandée est de 1e-05 m
*
* ===================================================================
*                          RESULTATS ATTENDUS
*
* Le champs de pression sera compare a la solution analytique a 
* 8 heures par l intermediare des erreurs L_infini et L_2.
* Une erreur L2 inferieure a 1.2 10-2 est attendue.
 
* ===================================================================
*
'OPTION' 'ECHO' 1 ;
'SAUTER' 'PAGE';
*
'TITRE' 'Transitoire homogene 1D : Dirichlet'   ;
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 ;
* pression imposee en haut initiale (m)
PSI1 = -1.1D0 ;
* pression imposee en haut t>0 (m)
PSI2 = -1.7D0 ; ;
*================================================
*-------------------------------------------------------------------
*--------------------- Création du maillage 2D ---------------------
*
*- Niveau de rafinement
*- Nombre de mailles = 50 * NIV1
*- Pas de temps = (Temps final) / (20 * (NIV1**2))
NIV1 = 2 ;
*- Densités de mailles
NDX1 = 1 ;
NDY1 = NIV1 * 50 ;
*
*- 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
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
*-----------------
*- Pression imposee initial P = HN1
HN1 = PSI1 + LONG1 ;
*- Pression imposee P = HN2
HN2 = PSI2 + LONG1 ;
HN2P = HN2 + LONG1 ;
* On impose la charge H = P + z (z=L en haut)
ESORH = 'MANU' 'CHPO' CESORH 1 'TH' HN2P 'NATURE' 'DISCRET';
VALI2 = 'CHAR' 'THIM' (ESORH)  ('EVOL' 'MANU' LICALC LIST1) ;
*
*--------------------------------------------------------------------
*----------------- initialisation des inconnues --------------------- 
*
*- charge initiale (en metres d'eau) dans le domaine (H = P + z)
* Pression h
ZAD1 = ZCC * ALPHA1 ;
*
L1 = ALPHA1 * LONG1 ;
CST1 = EXP(ALPHA1 * HN0) ;
CST2A = EXP(ALPHA1 * HN1) ;
CST2B = EXP(ALPHA1 * HN2) ;
CST3 = 1.0D0 - (EXP(-1.0D0*L1)) ;
CST4 = (CST1 - CST2A) / CST3 ;
CST5 = (CST1 - CST2B) / CST3 ;
*
EAP1 = ((EXP(-1. * ZAD1)) - 1.0D0) * CST4 + CST1 ;
PM   = (LOG(EAP1)) / ALPHA1 ;
*
* Traces de controle condition initiale en pression
* valeur aux noeuds
ZAD2 = ZSS * ALPHA1 ;
EAP2 = ((EXP(-1. * ZAD2)) - 1.0D0) * CST4 + CST1 ;
EAP3 = ((EXP(-1. * ZAD2)) - 1.0D0) * CST5 + CST1 ;
PM2  = LOG(EAP2) * (1. / ALPHA1) ;
PM3  = LOG(EAP3) * (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 ;
PINI2 = REDU PM3 LIN2 ;
EVIN1 = EVOL VERT 'CHPO' PINI1 LIN2 ;
EVIN2 = EVOL ROUG 'CHPO' PINI2 LIN2 ;
'SI' GRAPH  ;
DESS (EVIN1 ET EVIN0 ET EVIN2) ;
'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 ET VALI2 ;
*
* GBM - NOUVEAU
SATUR . 'LUMP'                = FAUX ;
SATUR . 'TYPDISCRETISATION'   = 'EFMH' ;
SATUR . 'DECENTR'             = VRAI ;
*
* 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 8 heures
TF1 = 8.0D0 * 3600.D0 ;
SATUR. 'TEMPS_FINAL'     = TF1 ;
SATUR. 'HOMOGENEISATION' = 'CHAINE' 'CENTRE' ;
SATUR. 'RESIDU_MAX'      = 1.D-5 ;
SATUR. 'NITER'           = 50 ;
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;
*
NBRPAS1 = (20. * (NIV1**2)) ;
DELTAT1 = TF1 / NBRPAS1 ;
SATUR . 'TEMPS_CALCULES' = ('PROG' 0. 'PAS' 1. NBRPAS1) * (DELTAT1) ;
*
*********************************************************************
*------------------------- Exécution procédure ----------------------
*********************************************************************
*
DARCYSAT SATUR ;
*
*********************************************************************
*
*-Legende des traces
TABI1 = TABLE ;
TABI1.'TITRE'= TABLE ;
TABI1.1      = 'REGU MARQ CARR' ;
TABI1.'TITRE' . 1  = MOT 'Initial 0 h' ;
TABI1.2      = 'REGU MARQ TRIA' ;
TABI1.'TITRE' . 2  = MOT 'Calcul 8 h' ;
TABI1.3      = 'REGU MARQ LOSA' ;
TABI1.'TITRE' . 3  = MOT 'Analytique 8 h' ;
TABI1.4      = 'REGU MARQ TRIB' ;
TABI1.'TITRE' . 4  = MOT 'Permanent' ;
*
*-------------------------------------------------------------------
*-------------------------- 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 ;
DESS (EVIN1 ET EVOL4 ET EVIN2) ;
'FINSI' ;
*
*--------------------------------------------------------------------------
* =========================================================================
* Solution analytique Richards lineaire Dirichlet
* u = exp(alpha.h)
* u(Z,t) = us(Z) + Somme ( Bn . sin(n.Pi.Z/L) . exp(-Z/2) . exp(-p.p.t) )
*--------------------------------------------------------------------------
*
* 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 ;
*
* Z addimentionnel 
ZCC1 = ZCCZ1 * ALPHA1 ;
*
*----------------------------------------------------------------------------
* Calcul du terme stationaire us(Z)
* us(Z) = TT1(Z) = exp(alpha*h0) - 
*                  (1-exp(-z))*(exp(alpha*h0)-exp(alpha*h2))/(1-exp(-L))
*----------------------------------------------------------------------------
CC1 = EXP(ALPHA1*HN0) ;
CC2 = EXP(ALPHA1*HN1) ;
CC3 = EXP(ALPHA1*HN2) ;
CC4 = EXP(-1.D0 * LNEW1) ;
CUZ1 = (CC3 - CC1) / (1.D0 - CC4) ;
TT1 = CC1 + (CUZ1 * (1.D0 - (EXP(-1.D0 * ZCC1)))) ;
CSTB1 = (CC2 - CC3) / (1.D0 - CC4) ;
BN1 = 2.D0 / LNEW1 ;
*
*----------------------------------------------------------------------------
* calcul du terme Somme
* Somme (Bn sin(n.Pi.z/L) exp(-z/2) exp(-p*p*t)
* Remarque: l'indice n=0 n'est pas calcule dans la somme
*----------------------------------------------------------------------------
* constantes trigo pour calcul
PI1 = PI ;
CSD1 = 180.D0 / PI1 ;
* 
* calcul de exp(-Z/2)
TLZ1 = KOPS ZCC1 '*' -0.5D0 ;
CSTE2 = EXP(TLZ1) ;
*
CSL1 = EXP(0.5D0 * LNEW1) ;
CSL2 = (CC3 - CC2) ;
NBOUC1 = 1000 ;
NNN1 = 1.D0 ;
SOM1 = 0.D0 ;
*---------------------------------------------
* calcul de la somme
*---------------------------------------------
REPETER BOU1 NBOUC1 ;
*
* Calcul des L_n et somme en sin(L_n.Z),
* avec ici L_n = n.Pi/L
*
  MESS 'Indice de INTG' NNN1 ;
*******************************
  LM1 = NNN1 * PI1 / LNEW1 ;
* 
*******************************
* calcul de SI1 = sin(L_n Z)
  MU1Z = KOPS ZCC1 '*' LM1 ;
* en degres
  MU2Z = MU1Z * CSD1 ;
  SI1 = SIN(MU2Z) ;
***************************************
* calcul de exp(-t p**2)
* avec p**2 = 1/4 + L_n**2
  LM2 =  LM1 * LM1 ;
  PCAR1 = 0.25D0 + LM2 ;
  TEX1 = -1.D0 * TAD1 * PCAR1 ;
  NP1 = EXP(TEX1) ;
***************************************
* calcul de Bn
* avec Bn = -(2/L) . (exp(alpha.h1)-exp(alpha.h2)) .
*            ((-1)**n) . L_n . exp(L/2) / p**2
  PMU1 = (-1.D0) ** NNN1 ;
  BN2 = (PMU1 * LM1) / (0.25D0 + LM2) ;
***************************************
* calcul produit total et somme
  COEF1 = SI1 * NP1 * BN2 ;
  SOM1 = SOM1 + COEF1 ;
  NNN1 = NNN1 + 1. ;
FIN BOU1 ;
***************************************
* fin de calcul du terme somme
*-----------------------------------------------
*
* Calcul de u
KTOT2 = BN1 * SOM1 * CSTE2 * CSL2 * CSL1 ;
KTOT1 = KTOT2 + TT1 ;
*
**************************************************************
*
* calcul de h associe a u
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 ET EVIN2) 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 'ERREUR L INFINI' ERINF1;
*
**********************
* 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 'ERREUR L2' ERL2 ;
*
**************************************************************
* Test de non regression sur erreur L2
ERMAX1 = 0.012 ;
SI (ERL2 > ERMAX1) ;
   mess 'ERREUR ANORMALE' ERL2 '>' ERMAX1 ;
  'ERREUR' 5;
SINON ;
  'ERREUR' 0;
FINSI ;
*
*-------------------------------------------------------------------
*
FIN;