Accueil | Cast3M

* fichier : gaculVF.dgibi
************************************************************************
************************************************************************
GRAPH = FAUX ;
******************** CAS TEST : gacul.dgibi ************************
*
*
* Test de fonctionnement de DARCYSAT en 2D avec effet de gravite
* ==============================================================
* Infiltration d'eau dans une colonne verticale 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 sont pris en compte (indice GRAVITE de
*    valeur VRAI);
*
* -  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')
*
* -  Une sous-relaxation de 0,8 est prise en compte (indice SOUSRELAX)
*
* -  L'homogenéisation spatiale des propriétes physiques s'effectue
*    par moyenne arithmétique des valeurs obtenues aux faces (indice
*    HOMOGENEISATION de valeur DECENTRE)
*
*!!!!!!!!!!!!!!!!!!!!!!!!
* -  La précision de convergence demandée est de 0,5 mm (indice RESIDU_MAX)
* 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 : gacul.dgibi ************************
 
'OPTION' 'ECHO' 1 ;
'SAUTER' 'PAGE';
*
'TITRE' 'infiltration verticale dans le sable : gacul.dgibi';
'OPTION' 'DIME' 2 'ELEM' 'QUA4' ;
'OPTION' 'ISOV' 'LIGN' 'TRACER' 'PSC' ;
 
*--------------------------------------------------------------------
* Création maillage
*
*- Discrétisation :
ENX = 1      ;
ENY = 200    ;
'DENS' (1./eny) ;
 
*- Création des points et des droites
*
A0 = -.5   1.0D0;  B0 = .5  1.0D0;
A1 = -.5   0.5D0;  B1 = .5  0.5D0;
A2 = -.5   0.25D0;  B2 = .5  0.25D0 ;
A3 = -.5   0.0D0;  B3 = .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 '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';
 
*- ordonnées
*
NDIME = 'VALEUR' 'DIME' ;
ZCC   = 'COOR' HYCEN NDIME ;
ZFF   = 'COOR' HYFAC NDIME ;
 
*- 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
 
*- frontière en limite du domaine de calcul (milieu désaturé)
ESORT = 'NOMC' (('REDU' ZFF CESORT) + HN) 'TH' 'NATURE' 'DISCRET';
 
 
*- frontière mouillée
EENTR = 'NOMC' ('REDU' ZFF CEENTR) 'TH' '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
*  (doit être compatible avec les conditions aux limites)
*
 
*- trace de charge d'eau
z      = 'REDU' (ZFF + HN) ('DIFF' CEENTR HYFAC) ;
 
 
*- charge d'eau
H0     = 'NOMC' (ZCC + HN) 'H' 'NATURE' 'DISCRET';
 
*- 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' = 20;
 
*- données géommétriques
SATUR. 'SOUSTYPE'         = 'DARCY_TRANSATUR' ;
SATUR. 'MODELE'           = MODHYB              ;
 
*- 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'   = 'VF' ;
 
* GBM MODIFIE
 
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'              = 7.3 ;
LoiP.'PERMSAT'            = 5.85E-2 / 3600. ;
SATUR.'LOI_PERMEABILITE'  = 'TABLE' LoiP ;
 
* loi de succion
LoiS                      = 'TABLE' 'VAN_GENUCHTEN';
LoiS. 'PORO'              = 0.3 ;
LoiS. 'TERESIDU'          = 0.055 ;
LoiS. 'NEXP'              = 2.0304  ;
LoiS. 'MEXP'              = 0.5075 ;
LoiS. 'BHETA'             = 0.029227 * 100. ;
SATUR.'LOI_SATURATION'    = 'TABLE' LoiS ;
 
* gravité (rho_w g,  en mètres d'eau, vaut 1)
SATUR. 'FORCE_GRAVITE'    = 1.;
 
*- données numériques
SATUR. 'SOUS_RELAXATION'  = 1.0D0;
SATUR. 'TEMPS_FINAL'      = 3000. ;
SATUR. 'HOMOGENEISATION'  = 'CHAINE' 'DECENTRE' ;
SATUR. 'NPAS'             = 1000 ;
SATUR. 'RESIDU_MAX'       = 1.D-2 ;
SATUR. 'NITER'            = 12 ;
SATUR. 'DT_INITIAL'       = 5. ;
SATUR. 'TEMPS_SAUVES'     = ('PROG' 1 'PAS' 1 6.)*500. ;
*SATUR. 'MESSAGES'         = 0 ;
 
 
'SI' GRAPH ;
*--------------------------------------------------------------------
*--- Visualisation de la loi de capacité capillaire.
* 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
  TAB = 'TABLE'  ;
  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' ;
 
*                                           ===========================
*                                           |         CALCUL          |
*                                           ===========================
DARCYSAT SATUR ;
 
* Post-traitement
* ===============
*-- Test de non régression
* Le test est effectué en vérifiant la solution à 1 heure de temps physique
 
*- solution a 1 heure de temps physique
lp = 'PROG' 35*-1.1500 -1.1501 -1.1496 -1.1516 -0.95203 -0.40387
            -0.22041 -0.15351 -0.11533 -8.84629E-02 -6.74391E-02
            -4.98338E-02 -3.43395E-02 -2.01365E-02 -6.63732E-03 ;
 
*- construction de la pression aux centres
SATUR.'PRESSION' = 'TABLE'  ;
'REPETER' bb (dime SATUR.TEMPS) ;
   id = &bb - 1 ;
   HN = SATUR.'CHARGE'.id ;
   SATUR.'PRESSION'.id = (nomc HN SCAL 'NATURE' 'DISCRET') - ZCC ;
'FIN' bb ;
 
*ev   = 'EVOL' 'CHPO' SATUR.'PRESSION'. 4 'SCAL' LCENC ;
*lp2  = 'EXTR' ev 'ORDO' ;
*err1 = 'MAXIMUM' ('ABS' (lp - lp2)) ;
*mess 'régression : ' err1 ;
 
*-- Tracé (tous les temps)
  LT        = 'LECT' 0  PAS 1 ((dime SATUR.TEMPS) - 1)  ;
  liopt     = 'MOTS' 'MIMA' 'AXES';
  TDES      = TRACHIS SATUR 'CHARGE' LCENC LT 'PREF' ' ' 'UNIT' 's' ;
'SI' GRAPH ;
  DESTRA TDES liopt 'TITX' 'z (m)' 'TITY' 'Pw (m)' ;
'FINSI' ;
 
toto = 'PRIM' TDES . 6 . VALEUR;
inta = 'EXTRAIRE' toto 'ORDO';
inta = 'EXTRAIRE'  ('DIME' inta) inta;
 
err = 'ABS' (inta '+' 0.235 '/' 0.235);
'LISTE' err;
 
'SI' (err > 2.D-2) ;
   'ERREUR'  5;
'FINSI'  ;
 
 
'FIN';