* 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 = 'ORDO' (orig 'ET' hycen) ; 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 ; *-------------------------------------------------------------------- * 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.e11 ; 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 'H' Profsite 'NATURE' 'DISCRET' ; Hi_bo = 'MANU' 'CHPO' ('DOMA' Mbo 'CENTRE') 1 'H' hc_bo 'NATURE' 'DISCRET' ; H_ini = Hi_site 'ET' Hi_bo ; Flux_ini = 'MANU' 'CHPO' ('DOMA' modhyb 'FACE') 1 'FLUX' 0. 'NATURE' 'DISCRET' ; * --------------------------- * = 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 ; *- conditions initiales SATUR. 'TEMPS' . 0 = Tini ; SATUR. 'TRACE_CHARGE' . 0 = Th_ini ; 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 . 'TYPDISCRETISATION'= 'EFMH' ; 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 = 'MAXIMUM' ('ABS' (lp - lp2)); 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' 'S (s.d.)' ; TDES = TRACHIS SATUR 'CHARGE' LCENC LT ('MOTS' 'H') ('MOTS' ' ') 'PREF' ' ' ; DESTRA TDES 'MIMA' 'AXES' 'TITX' 'z (m)' 'TITY' 'Pw (m)' ; 'FINSI' ; 'FIN';