Accueil | Cast3M

* fichier :  transsat.dgibi
************************************************************************
************************************************************************
GRAPH = vrai ;
******************** CAS TEST : transsat.dgibi ************************
*
*
* Test de fonctionnement de DARCYSAT sur un probleme multizone.
* ===================================================================
* Infiltration d'eau dans une barrière ouvragée dans son site d'acceuil.
*
* -  condition initiale :
*    > BO : désaturation uniforme correspondant a une succion de 7827,42 m ;
*    > SITE saturé sous 500m d'eau
* - les paramètres des lois de perméabilité diffèrent entre le SITE et la BO
*
* ===================================================================
*    Les options de modélisation déclarées dans la table transmise à
*    la procedure DARCYSAT sont les suivantes :
*
* -  les effets gravitationnels ne sont pas pris en compte : pression =
*    charge (indice GRAVITE absent ; valeur par defaut : FAUX);
*
* -  Le pas de temps est automatique (indice TEMPS_CALCULES absent)
*    sont fournis :
*    > le pas de temps initial (indice 'DT_INITIAL'),
*    > le nombre d'itérations recherchées par pas de temps (indice 'NITER')
*    > le nombre de pas de temps maximum (indice 'NPAS')
*    > le nombre de CFL initial (indice 'CFL')
*
* -  La méthode de pénalisation est désactivée (indice XI absent) ;
*
* -  Homogénéisation : les propriétés physiques sont calculées à partir
*    des charges aux éléments (indice HOMOGENEISATION de valeur CENTRE)
*
*
************************************************************************
'OPTION' 'ECHO' 1 ;
'SAUTER' 'PAGE';
*
 
TITRE 'Saturation d une barriere ouvragee : transsat.dgibi'   ;
'OPTION' 'DIME' 2 'ELEM' 'QUA4';
'OPTION' 'ISOV' 'LIGN' 'TRACER' 'PSC' ;
 
*--------------------------------------------------------------------
* Création maillage
*
*- Discrétisation :
ENX = 1      ;
ENY = 100    ;
'DENS' (3./eny) ;
 
*- Création des points et des droites
*
A0 = -.5   6.D0 ;  B0 = .5   6.D0 ;
A1 = -.5   4.D0 ;  B1 = .5   4.D0 ;
A2 = -.5   1.D0 ;  B2 = .5   1.D0 ;
A3 = -.5   0.D0 ;  B3 = .5   0.D0 ;
 
*- Création des droites
*
AB0 = 'DROIT' ENX A0 B0 ;
AB1 = 'DROIT' ENX A1 B1 ;
AB2 = 'DROIT' ENX A2 B2 ;
AB3 = 'DROIT' ENX A3 B3 ;
 
*- Création des surfaces
*
BO      = 'REGLER' AB3 AB2 ;
BO      = 'COULEUR' BO 'VERT' ;
BO      = 'INVERSE' BO ;
SITE    = AB2 'REGLER' DINI (3./eny) DFIN 0.2 AB1 'REGLER' 10 AB0 ;
SITE    = 'COULEUR' SITE 'BLEU' ;
SITE    = 'INVERSE' SITE ;
MASSIF0 = BO 'ET' SITE;
 
* entrée
ENTR    = 'COULEUR' ('INVERSE' AB0) 'ROUGE' ;
 
'SI' GRAPH  ;
  'TRACER' (MASSIF0 'ET' ENTR) ;
'FINSI' ;
 
*- Création des maillages contenant tous les points
*************************************************************************
* IMPORTANT : ne pas concaténer de maillages QUAF                       *
*************************************************************************
*
QMASSIF = 'CHANGER' MASSIF0 'QUAF' ;
QSITE   = 'CHANGER' SITE    'QUAF' ;
QBO     = 'CHANGER' BO      'QUAF' ;
QFENTR  = 'CHANGER' ENTR    'QUAF' ;
 
'ELIMINATION' (QMASSIF 'ET' QFENTR 'ET' QBO 'ET' QSITE) 0.00001 ;
 
*- Modèles
*************************************************************************
* IMPORTANT : Il faut concaténer les modèles et ne surtout pas créer un *
* modèle total à une seule zone                                         *
*************************************************************************
MBO     = 'MODELE' QBO     'DARCY' 'ISOTROPE'        ;
MSITE   = 'MODELE' QSITE   'DARCY' 'ISOTROPE'        ;
MENTR   = 'MODELE' QFENTR  'DARCY' 'ISOTROPE'        ;
MODHYB  =  MSITE 'ET' MBO;
 
CEENTR  = 'DOMA' MENTR  'CENTRE' ;
HYCEN   = 'DOMA' MODHYB 'CENTRE' ;
HYFAC   = 'DOMA' MODHYB 'FACE';
 
*- Création ligne de suivi pour le post-traitement
* ligne des points centres (cas 1D)
*
* on ordonnance de bas en haut les points centres
orig  = 0. 0. ;
lcen1 = orig 'ET' hycen ;
'ORDO' lcen1 ;
lcen2 = 'DIFF' lcen1 ('MANU' 'POI1' orig) ;
* et on crée la ligne de seg2
lcenc = 'QUELCONQUE' lcen2 'SEG2' ;
 
*- ordonnées
*
NDIME = 'VALEUR' 'DIME' ;
ZCC   = 'COOR' HYCEN NDIME ;
ZFF   = 'COOR' HYFAC NDIME ;
 
*--------------------------------------------------------------------
*
* Paramètres
*
unjour   = 24 * 3600.      ;
unmois   = 30 * unjour     ;
unan     = unjour * 365.25 ;
Tini     = 0.              ;
Tfin     = 1000 * unan     ;
 
*--- profondeur du site
Profsite = 500.            ;
T        = 50 + 273.15     ;
hc_bo    = -7827.42        ;
RHOg     = 9810. ;
poro_bo  = 0.3 ;
 
*--- coefficient d'emmagasinement
* C = poro * rhoeau * g *
*    (compressibilite_eau - compressibilite_solide + alpha/poro)
* avec alpha = 3 (1 - 2 Nupoisson) / Eyoung
* pour FOCA :  nu = 0.25 ; E = 500 MPa
* poro = 0.3 ; RHOg = 9810 ; compres_eau = 5E-10 Pa-1 ; compres_solide = 0
emma_sit = 3.83E-6 ;
nu_bo    = 0.4   ;
Young_bo = 150.E6 * 2 * (1 + nu_bo) ;
alph_bo  = 3 * (1 - (2 * nu_bo)) / Young_bo  ;
emma_bo  = 1.D0 *poro_bo * RHOg * (5E-10 - 0 + (alph_bo/poro_bo))     ;
emma0    = ('MANU' 'CHPO' ('DOMA' msite 'CENTRE') 1
                   'SCAL' emma_sit 'NATURE' 'DISCRET')
      'ET' ('MANU' 'CHPO' ('DOMA' mbo 'CENTRE')  1
                   'SCAL' emma_bo  'NATURE' 'DISCRET') ;
emma0     = 'KCHT' modhyb 'SCAL' 'CENTRE' 'COMP' 'SCAL' emma0  ;
 
*emma0 = 1.D-2 *  emma0 ; 
 
 
*--------------------------------------------------------------------
* Conditions aux limites
* 1) flux nul sur la frontiere du colis => rien a faire
* 2) flux nuls en latéral => rien a faire
* 3.1) pression imposee en limite de site Pliq = h avec h =500m
BENTR = 'BLOQUE' CEENTR 'TH'  ;
EENTR = 'MANUEL' 'CHPO' CEENTR 1 'TH'  Profsite  ;
 
* chargement des CLs
LICALC = 'PROG' 0.D0  1.e10 ;
LIST1  = 'PROG' 2 * 1.D0 ;
VALI0  = 'CHAR' 'THIM' EENTR  ('EVOL' 'MANU' LICALC LIST1) ;
 
*--------------------------------------------------------------------
*- initialisation des inconnues
*  (doit être compatible avec les conditions aux limites)
*
Thi_site = 'MANU' 'CHPO' ('DOMA'  Msite 'FACE') 1
                  'TH'  Profsite 'NATURE' 'DISCRET' ;
Thi_bo   = 'MANU' 'CHPO' ('DOMA'  Mbo 'FACE') 1
                  'TH' hc_bo 'NATURE' 'DISCRET'   ;
Th_ini   = 'KCHT' modHYB 'SCAL' 'FACE' 'COMP' 'TH' Thi_site Thi_bo;
*
Hi_site  = 'MANU' 'CHPO' ('DOMA' Msite 'CENTRE') 1
                  'SCAL'  Profsite 'NATURE' 'DISCRET' ;
Hi_bo    = 'MANU' 'CHPO' ('DOMA' Mbo 'CENTRE') 1
                  'SCAL' hc_bo 'NATURE' 'DISCRET'    ;
H_ini    = kcht modhyb SCAL CENTRE (Hi_site 'ET' Hi_bo) ;
H_ini    = nomc 'H' H_INI ;
Flux_ini = 'MANU' 'CHPO' ('DOMA' modhyb 'FACE') 1
                  'FLUX' 0. 'NATURE' 'DISCRET' ;
 
*                                         ---------------------------
*                                         = Table DARCY_TRANSITOIRE =
*                                         ---------------------------
*- initialisation table
SATUR                     = 'TABLE' ;
SATUR. 'TEMPS'            = 'TABLE' ;
SATUR. 'CHARGE'           = 'TABLE' ;
SATUR. 'FLUX'             = 'TABLE' ;
 
*- données géométriques
SATUR. 'SOUSTYPE'         = 'DARCY_TRANSATUR' ;
SATUR. 'MODELE'           = MODHYB ;
 
*- conditions initiales
SATUR. 'TEMPS' . 0        = Tini  ;
SATUR. 'CHARGE' . 0       = H_ini  ;
SATUR. 'FLUX' . 0         = Flux_ini             ;
 
*- conditions aux limites et chargements
SATUR. 'TRACE_IMPOSE'       = VALI0  ;
 
 
*- définition de la loi de perméabilité
SATUR.'LOI_PERMEABILITE'  = 'TABLE' 'MULTIZONE' ;
 
* pour le site
LoiP                      = 'TABLE' 'PUISSANCE';
LoiP . 'ALPHA'            = RHOg / 1.D7 ;
LoiP . 'PERMSAT'          = 1.D-12  ;
LoiP . 'MODELE'           = MSITE ;
SATUR. 'LOI_PERMEABILITE'. 'SITE' = LoiP ;
 
* pour la BO
LoiP                      = 'TABLE' 'PUISSANCE';
LoiP .'ALPHA'             = RHOg / (7.D6);
LoiP .'PERMSAT'           = 6.D-14  ;
LoiP .'MODELE'            = MBO ;
SATUR.'LOI_PERMEABILITE'. 'BO' = LoiP ;
 
*- définition de la teneur en eau
SATUR.'LOI_SATURATION'    = 'TABLE' 'MULTIZONE' ;
 
* pour le site
LoiS                      = 'TABLE' 'VAN_GENUCHTEN';
LoiS . 'PORO'             = 0.05 ;
LoiS . 'TERESIDU'         = 0. ;
LoiS . 'NEXP'             = 1.7  ;
LoiS . 'MEXP'             = 1 - (1 / 1.7) ;
LoiS . 'BHETA'            = RHOg / 10.D6 ;
LoiS . 'MODELE'           = MSITE ;
SATUR.'LOI_SATURATION' . 'SITE' = LoiS ;
 
'SI' GRAPH  ;
*-  Vérification du choix du dp minimum pour le calcul de la capacité.
* droite support des variables
  dx  = 'DROIT' (0. -7827.4 ) 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 dum = HT_PRO (SATUR.'LOI_SATURATION'.'SITE') 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. ;
s1 t1 dum = HT_PRO (SATUR.'LOI_SATURATION'.'SITE') ('KOPS' zc '-'  dp);
* représentation de la capacité
  c1  = (t0 - t1) / dp;
  ev1 = 'EVOL' 'ROUGE' 'CHPO' c1 'SCAL' dx ;
  evc = 'EVOL' 'TURQ'  'MANU' px ('EXTR' ev1 'ORDO') ;
  'DESSIN' evc 'TITX' 'Pc(m)' 'TITY' 'Capa(1/m)'
               'TITRE' 'Capacite capillaire SITE' ;
FINSI ;
 
* pour la BO
LoiS                      = TABLE 'VAN_GENUCHTEN';
LoiS . 'PORO'             = poro_bo ;
LoiS . 'TERESIDU'         = 0. ;
LoiS . 'NEXP'             = 1. / (1. - 0.27)  ;
LoiS . 'MEXP'             = 0.27 ;
LoiS . 'BHETA'            = RHOg / 17.E6 ;
LoiS . 'MODELE'           = MBO ;
SATUR.'LOI_SATURATION' . 'BO' = LoiS ;
 
'SI' GRAPH  ;
*-  Vérification du choix du dp minimum pour le calcul de la capacité.
*--- calcul de la teneur en eau pour p = zc
* calcul de la teneur en eau pour la pression zc
  s0 t0 dum = HT_PRO (SATUR.'LOI_SATURATION'.'BO') 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. ;
s1 t1 dum = HT_PRO (SATUR.'LOI_SATURATION'.'BO') ('KOPS' zc '-'  dp);
* représentation de la capacité
  c1  = (t0 - t1) / dp;
  ev1 = 'EVOL' 'ROUGE' 'CHPO' c1 'SCAL' dx ;
  evc = 'EVOL' 'TURQ'  'MANU' px ('EXTR' ev1 'ORDO') ;
  'DESSIN' evc 'TITX' 'Pc(m)' 'TITY' 'Capa(1/m)'
               'TITRE' 'Capacite capillaire B.O.' ;
FINSI ;
SATUR. 'COEF_EMMAGASINEMENT' = emma0    ;
 
* données numériques
SATUR . 'ITMAX' = 40;
SATUR. 'TEMPS_FINAL'      = 1.e10 ;
SATUR. 'HOMOGENEISATION'  = 'CENTRE' ;
SATUR. 'NPAS'             = 100 ;
SATUR . 'RESIDU_MAX'      = 1.D-4;
SATUR. 'DT_INITIAL'       = 12 * 3600. ;
SATUR. 'MESSAGES'         = 1 ;
SATUR . 'LUMP' = FAUX    ;
 
SATUR . 'LUMP'   = FAUX;
SATUR . 'TYPDISCRETISATION'   = 'VF' ;
SATUR . 'DECENTR'  = FAUX;
TABRES = table METHINV;
TABRES . 'TYPINV' = 3;
TABRES . 'PRECOND' = 3;
 
SATUR . 'METHINV' = TABRES;
SATUR . 'METHINV' . RESID = 1.D-20;
 
SATUR. 'TEMPS_SAUVES'    = ('PROG' 1 'PAS' 1 41.)*36000. ;
SATUR. 'TEMPS_CALCULES'    = ('PROG' 1 'PAS' 1 41.)*36000. ;
 
 
* ===============
* Calcul
* ===============
DARCYSAT SATUR ;
 
* ===============
* Post-traitement
* ===============
* Le test est effectué en vérifiant la solution
 
*- solution après 40 pas de temps
lp = 'PROG'  -7827.4 -7827.4 -7827.4 -7827.4 -7827.4 -7827.4 -7827.4
             -7827.4 -7827.4 -7827.4 -7827.3 -7827.2 -7826.9 -7826.3
             -7825.2 -7823.1 -7819.2 -7812.2 -7800.1 -7779.9 -7747.0
             -7695.1 -7615.8 -7498.0 -7328.3 -7090.4 -6765.7 -6334.1
             -5773.9 -5064.7 -4189.3 -3141.8 -1944.8 -1272.1 -1188.6
             -1101.0 -1009.7 -915.24 -818.60 -720.85 ;
lp = lp 'ET' ('PROG'
             -623.36 -527.60 -435.04 -346.88 -263.88 -186.13 -113.58
             -47.363 14.082  73.448  130.79  185.41  236.60  283.74
             326.27  363.79  396.07  423.06  444.92  462.03  474.90
             484.16  490.51  494.64  497.16  498.60  499.36  499.72
             499.88  499.95  499.98  499.99  500.00  500.00  500.00
             500.00);    
 
 
 
id   = 40 ;
ev   = 'EVOL' 'CHPO' SATUR.'CHARGE'.id 'H' LCENC ;
lp2  = 'EXTR' ev 'ORDO' ;
err1 = SOMM ('ABS' (lp - lp2));
err1 = err1 / (SOMM (ABS lp));
mess 'régression : ' err1 ;
 
'LISTE' lp2 ;
 
 
'SI' (err1 > 5.D-1);
   mess  'erreur anormale après ' (@ARR id 0) ' pas : ' err1 ;
*   'ERREUR' 5;
'SINON' ;
  'ERREUR' 0;
'FINSI' ;
 
* ===============
* Post-traitement
* ===============
 
'SI' GRAPH  ;
*--- liste des indices des temps postraites
   NN        = (('DIME' SATUR.TEMPS) - 1) / 3 ;
   LT        = 'LECT' 0 'PAS' 3 (3 * NN)  ;
   TDES      = TRACHIS SATUR 'SATURATION' LCENC LT 
                       ('MOTS' 'SCAL') ('MOTS' ' ') 'PREF' ' '   ;
   DESTRA TDES 'MIMA' 'AXES' 'TITX' 'z (m)' 'TITY' 'Pw (m)' ;
   TDES      = TRACHIS SATUR 'CHARGE' LCENC LT 
                       ('MOTS' 'H') ('MOTS' ' ') 'PREF' ' '   ;
   DESTRA TDES 'MIMA' 'AXES' 'TITX' 'z (m)' 'TITY' 'Pw (m)' ;
'FINSI' ;
 
'FIN';