* 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 = RESU ('ABS' (lp - lp2));
err1 = err1 / (RESU (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';
© Cast3M 2003 - All rights reserved.
Disclaimer