* fichier : alonso.dgibi
******************************************************************
* *
* Modèle d'ALONSO: argile *
* *
* Un cube d'argile de 50 cm de cote est soumis à des *
* deplacements imposes de 10.cm sur 3 de ses faces. *
* L'état de contraintes est volumique. *
* *
* Ce calcul est réalisé dans le domaine saturé de l'argile *
* (succion nulle) *
* *
******************************************************************
'OPTION' 'ECHO' 0 'DIME' 3 'ELEM' cub8 'TRACER' psc ;
* ajout de option epsilon lineaire pour la precision des test!
OPTION epsilon lineaire;
graph0 = vrai ;
ecri0 = vrai ;
elast0 = faux ;
pb0 = faux ;
xlim0 = 1.E-7 ;
*
* Maillage : un cube soumis à un déplacement
* isotrope
*
n1 = 1 ;
n2 = 1 ;
n3 = 1 ;
l1 = .5 ;
l2 = .5 ;
l3 = .5 ;
p1 = 0. 0. 0. ;
p2 = l1 0. 0. ;
p3 = l1 l2 0. ;
p4 = 0. l2 0. ;
d1 = 'DROITE' p1 n1 p2 ;
d2 = 'DROITE' p2 n2 p3 ;
d3 = 'DROITE' p3 n1 p4 ;
d4 = 'DROITE' p4 n2 p1 ;
sur1 = 'DALLER' d1 d2 d3 d4 'PLAN' ;
mail0 = sur1 'VOLU' n3 'TRAN' ( 0. 0. l3 ) ;
*
coo1 coo2 coo3 = 'COOR' mail0 ;
*
sur1 = coo3 'POIN' 'EGSUPE' l3 ;
sur2 = coo2 'POIN' 'EGSUPE' l2 ;
sur3 = coo1 'POIN' 'EGSUPE' l1 ;
sur4 = coo1 'POIN' 'EGINFE' 0. ;
sur5 = coo2 'POIN' 'EGINFE' 0. ;
sur6 = coo3 'POIN' 'EGINFE' 0. ;
*
*mail0 = sur1 'ET' sur2 'ET' sur3 'ET' sur4 'ET'
* sur5 'ET' sur6 'ET' mail0 ;
'ELIM' mail0 1.E-5 ;
'SI' graph0 ;
'TITR' ' Maillage ' ;
'TRAC' mail0 'CACH' ;
'FINSI' ;
*
* Initialisations diverses
*
* Longueur du déplacement imposé sur les côtés 1, 2 et 3 (m) :
depm0 = -.1 ;
tfin0 = 10. ;
EV0 = 'EVOL' 'MANU' 'temps' ( 'PROG' 0. tfin0 ) 'y' ( 'PROG' 0. 1. ) ;
* Liste des temps calculés
dt0 = tfin0 / 10. ;
tem0 = 'PROG' 0. 'PAS' dt0 tfin0 ;
*
* modèle et matériau
* ==================
* {Notation d'Alonso : "A constitutive model for partially saturated soils"
* E. Alonso, A. Gens, A. Josa (1990), Géotechnique 40, n°3}
* [notations de Gens : "Expertise de la problématique THM-Gaz dans le
* cadre de l'étude de la faisabilité d'un stockage en
* profondeur. UPC" Rapport ANDRA : C RP 0UPC 00-001.
* p' = contrainte nette isostatique]
*
*-- Lois de comportement :
* module d'Young tangent équivalent initial [E] (Pa).
* Il sert uniquement pour initialiser le premier pas de temps.
* N'importe quelle valeur supérieure à la pente à l'origine convient.
youn0 = 2.1E11 ;
* module de cisaillement équivalent initial [G] {G} (Pa)
* Alonso donne 10 MPa !
gtan0 = 150.E6 ;
* coefficient de Poisson [nu] (s.d.) :
nu0 = .3 ;
* Coef. d'écrouissage de traction en succion [k de P_s = - k s] {idem} (s.d.)
* Alonso donne 0.6 !
kk0 = 1.E-2 ;
* Pente q/p à l'état critique [M] {M} (s.d.) :
* Alonso donne 1 !
m0 = 1.2 ;
* paramètres de la loi lambda(s) = lambda(0) x [(1-r) exp (-beta s) + r]
* Compressibilité en pression isostatique sur un chemin plastique :
* charge d'un ech. normalement consolidé
* [lambda(0)] {idem} (s.d.) :
xlam1 = 5. / 111. ;
* [r] {r} :
* Alonso donne 0.75 !
xt0 = .5 ;
* [beta] {beta} (Pa-1) :
tau0 = 1.E-6 ;
* Compressibilité élastique initiale en p' à saturation [k_i^0] {kapa} :
* Il manque alpha_is pour définir k_i [Gens] {Alonso donne kapa constant}.
ktan0 = 1. / 111. ;
* Pente d'écrouissage en succion [lambda_s] {lambda_s} (s.d.)
* charge d'un ech. normalement consolidé
* [variable chez Gens] {constant chez Alonso}
lams0 = 5. / 111. ;
* Compressibilité élastique en succion [k_s] (s.d.)
* charge d'un ech. surconsolidé ou décharge
* Alonso donne 0.008 au lieu de 0.009. Gens le fait varier
ks0 = 1. / 111. ;
* pression de référence pour l'évolution de P0(s) [P^c] {P^c} (Pa)
xp0 = .1E6 ;
*-- Etat initial de l'échantillon
* Contrainte de compression isostatique initiale [p] {sigma_m} (Pa)
xp1 = 1.E6 ;
* succion initiale [s(t=0)] {s(t=0)} (Pa)
succ0 = 0. ;
* indice des vides initial [e(t=0)] {e(t=0)} (s.d.)
e00 = 0. ;
* pression de consolidation fictive (à saturation) initiale [P0 ou Pc(t=0)] {idem} (Pa)
xpc0 = 20.E6 ;
* limite initiale de plasticité en succion (Succion Increase) [s_0] {s_0} (Pa) ?
s00 = 1.E7 ;
*
'SI' elast0 ;
* Calcul élastique
'MESS' ' Calcul élastique' ;
mod0 = 'MODE' mail0 'MECANIQUE' 'ELASTIQUE' ;
mat0 = 'MATE' mod0 'YOUN' youn0 'NU' nu0 ;
'SINON' ;
* Calcul sur le modèle argileux ALONSO
'MESS' ' Calcul sur le modèle argileux ALONSO ' ;
mod0 = 'MODE' mail0 'MECANIQUE' 'ELASTIQUE' 'PLASTIQUE' 'ALONSO' ;
mat0 = 'MATE' mod0 'YOUN' youn0 'NU' nu0 'KS' ks0
'LAMS' lams0 'M' m0 'KK' kk0 'PC0' xpc0
'P0' xp0 'LAM1' xlam1 'T0' xt0 'TAU' tau0
'GTAN' gtan0 'KTAN' ktan0 'E0' e00 ;
*
* Volumes élémentaires de la structure
*
chel0 = 'MANU' 'CHML' mail0 'SCAL' 1. ;
vol0 = 'INTG' 'ELEM' mod0 chel0 ;
*
* Variables internes initiales
*
q01 = 'MANU' 'CHML' mod0
* 'VOLU' vol0 'PC0' xpc0 'SLIM' s00 'SUCC' succ0
'VOLU' vol0 'SLIM' s00 'SUCC' succ0
'TYPE' 'VARINTER' 'STRESSES' ;
q00 = 'ZERO' mod0 'VARINTER' ;
qq00 = q00 + q01 ;
*
* Contraintes initiales
*
p01 = 'MANU' 'CHML' mod0
'SMXX' (-1. * xp1) 'SMYY' (-1. * xp1) 'SMZZ' (-1. * xp1)
'SMXY' 0. 'SMXZ' 0. 'SMYZ' 0.
'TYPE' 'CONTRAINTES' 'STRESSES' ;
*
'FINSI' ;
*
* Conditions aux limites
*
cl1 = 'BLOQ' sur1 'UZ' ;
cl2 = 'BLOQ' sur2 'UY' ;
cl3 = 'BLOQ' sur3 'UX' ;
cl4 = 'BLOQ' sur4 'UX' ;
cl5 = 'BLOQ' sur5 'UY' ;
cl6 = 'BLOQ' sur6 'UZ' ;
cl0 = cl1 'ET' cl2 'ET' cl3 'ET' cl4 'ET' cl5 'ET' cl6 ;
*
* Chargement
*
dep1 = 'DEPI' cl1 depm0 ;
dep2 = 'DEPI' cl2 depm0 ;
dep3 = 'DEPI' cl3 depm0 ;
cha1 = 'CHAR' 'DIMP' dep1 EV0 ;
cha2 = 'CHAR' 'DIMP' dep2 EV0 ;
cha3 = 'CHAR' 'DIMP' dep3 EV0 ;
cha0 = cha1 'ET' cha2 'ET' cha3 ;
*
* Appel à PASAPAS
*
ta1 = 'TABLE' ;
ta1 .'MODELE' = mod0 ;
ta1 .'CARACTERISTIQUES' = mat0 ;
ta1 .'CHARGEMENT' = cha0 ;
ta1 .'BLOCAGES_MECANIQUES' = cl0 ;
'SI' ( 'NON' elast0 ) ;
ta1 .'VARIABLES_INTERNES' = 'TABLE' ;
ta1 .'VARIABLES_INTERNES'. 0 = qq00 ;
ta1 .'CONTRAINTES' = 'TABLE' ;
ta1 .'CONTRAINTES'. 0 = p01 ;
'FINSI' ;
ta1 .'TEMPS_CALCULES' = tem0 ;
*
TMASAU=table;
ta1 . 'MES_SAUVEGARDES'=TMASAU;
TMASAU .'DEFTO'=VRAI;
TMASAU .'DEFIN'=VRAI;
TMASAU .'DEFLO'=VRAI;
PASAPAS ta1 ;
*
'SI' elast0 ;
*
* Post traitement
*
dim0 = 'DIME' ta1 .'TEMPS' - 1 ;
sigm0 = ta1 .'CONTRAINTES'. dim0 ;
depl0 = ta1 .'DEPLACEMENTS'. dim0 ;
*
sigxx = 'EXCO' sigm0 'SMXX' 'SCAL' ;
sigyy = 'EXCO' sigm0 'SMYY' 'SCAL' ;
sigzz = 'EXCO' sigm0 'SMZZ' 'SCAL' ;
pres0 = -1. * ( sigxx + sigyy + sigzz ) / 3. ;
*
def1 = 'DEFO' mail0 depl0 1. 'ROUG' ;
def0 = 'DEFO' mail0 depl0 0. ;
*
'SI' graph0 ;
'TITR' ' Deformee ' ;
'TRAC' ( def1 'ET' def0 ) ;
'FINSI' ;
*
'SI' ecri0 ;
'MESS' ' sigxx = ' ( 'MAXI' sigxx ) ;
'MESS' ' sigyy = ' ( 'MAXI' sigyy ) ;
'MESS' ' sigzz = ' ( 'MAXI' sigzz ) ;
'MESS' ' pres0 = ' ( 'MAXI' pres0 ) ;
'FINSI' ;
*
'SINON' ;
*
dim0 = ( 'DIME' ta1 .'TEMPS' ) - 1 ;
*
ktan0 = 'MAXI' ( 'EXCO' mat0 'KTAN' 'SCAL' ) ;
xlam1 = 'MAXI' ( 'EXCO' mat0 'LAM1' 'SCAL' ) ;
xks0 = 'MAXI' ( 'EXCO' mat0 'KS' 'SCAL' ) ;
lams0 = 'MAXI' ( 'EXCO' mat0 'LAMS' 'SCAL' ) ;
xt0 = 'MAXI' ( 'EXCO' mat0 'T0' 'SCAL' ) ;
tau0 = 'MAXI' ( 'EXCO' mat0 'TAU' 'SCAL' ) ;
p00 = 'MAXI' ( 'EXCO' mat0 'P0' 'SCAL' ) ;
*
depl0 = ta1 .'DEPLACEMENTS'. dim0 ;
sigm0 = ta1 .'CONTRAINTES'. dim0 ;
var0 = ta1 .'VARIABLES_INTERNES'. dim0 ;
*
sigxx = 'EXCO' sigm0 'SMXX' 'SCAL' ;
sigyy = 'EXCO' sigm0 'SMYY' 'SCAL' ;
sigzz = 'EXCO' sigm0 'SMZZ' 'SCAL' ;
pres0 = -1. * ( sigxx + sigyy + sigzz ) / 3. ;
*
def1 = 'DEFO' depl0 mail0 1. 'ROUG' ;
def0 = 'DEFO' depl0 mail0 0. ;
'SI' graph0 ;
'TITR' ' Deformee ' ;
'TRAC' ( def1 'ET' def0 ) ;
'FINSI' ;
*
err_t1 = 'PROG' 0. ;
err_p1 = 'PROG' 0. ;
pres1 = 'PROG' xp1 ;
treps1 = 'PROG' 0. ;
pc1 = 'PROG' xpc0 ;
trepp1 = 'PROG' 0. ;
pres2 = 'PROG' xp1 ;
trepp2 = 'PROG' 0. ;
succ1 = 'PROG' succ0 ;
slim1 = 'PROG' s00 ;
succ00 = succ0 ;
xlam2 = xlam1 * ((1. - xt0) *
( 'EXP' ( -1.* tau0 * succ0)) + xt0) ;
exp0 = ( xlam1 - ktan0 ) / ( xlam2 - ktan0 ) ;
pcs0 = p00 * ( ( xpc0 / p00 ) ** exp0 ) ;
pcs1 = 'PROG' pcs0 ;
*
'REPETER' bloc0 dim0 ;
i0 = &bloc0 ;
'SI' ecri0 ;
'MESS' ' Calcul de restitution au pas ' i0 ;
'FINSI' ;
*
depl0 = ta1 .'DEPLACEMENTS'. i0 ;
sigm0 = ta1 .'CONTRAINTES'. i0 ;
var0 = ta1 .'VARIABLES_INTERNES'. i0 ;
def0 = ta1 .'DEFORMATIONS_INELASTIQUES'. i0 ;
*
sigxx = 'EXCO' sigm0 'SMXX' 'SCAL' ;
sigyy = 'EXCO' sigm0 'SMYY' 'SCAL' ;
sigzz = 'EXCO' sigm0 'SMZZ' 'SCAL' ;
pres0 = 'MAXI' ( -1. * ( sigxx + sigyy + sigzz ) / 3. ) ;
*
eps0 = 'EPSI' mod0 depl0 ;
epxx = 'EXCO' eps0 'EPXX' 'SCAL' ;
epyy = 'EXCO' eps0 'EPYY' 'SCAL' ;
epzz = 'EXCO' eps0 'EPZZ' 'SCAL' ;
treps0 = 'MAXI' ( -1. * ( epxx + epyy + epzz ) ) ;
*
pc0 = 'MAXI' ( 'EXCO' var0 'VPC0' 'SCAL' ) ;
*
defxx = 'EXCO' def0 'EIXX' 'SCAL' ;
defyy = 'EXCO' def0 'EIYY' 'SCAL' ;
defzz = 'EXCO' def0 'EIZZ' 'SCAL' ;
trepp0 = 'MAXI' ( -1. * ( defxx + defyy + defzz ) ) ;
*
pres1 = pres1 'ET' ( 'PROG' pres0 ) ;
pres01 = pres0 ;
treps1 = treps1 'ET' ( 'PROG' treps0 ) ;
pc1 = pc1 'ET' ( 'PROG' pc0 ) ;
trepp1 = trepp1 'ET' ( 'PROG' trepp0 ) ;
trepp01 = trepp0 ;
*
succ0 = 'MAXI' ( 'EXCO' var0 'SUCC' 'SCAL' ) ;
slim0 = 'MAXI' ( 'EXCO' var0 'SLIM' 'SCAL' ) ;
pc0 = 'MAXI' ( 'EXCO' var0 'VPC0' 'SCAL' ) ;
xlam2 = xlam1 * ((1. - xt0) *
( 'EXP' ( -1.* tau0 * succ0)) + xt0) ;
exp0 = ( xlam1 - ktan0 ) / ( xlam2 - ktan0 ) ;
pcs0 = p00 * ( ( pc0 / p00 ) ** exp0 ) ;
succ1 = succ1 'ET' ( 'PROG' succ0 ) ;
slim1 = slim1 'ET' ( 'PROG' slim0 ) ;
pcs1 = pcs1 'ET' ( 'PROG' pcs0 ) ;
*
* Solution analytique
*
treps0 = -1. * treps0 ;
p_test0 = xp1 * ( 'EXP' ( -1. / ktan0 * treps0 ) ) ;
xlam2 = xlam1 * ((1. - xt0) *
( 'EXP' ( -1.* tau0 * succ0)) + xt0) ;
rapp0 = ( xpc0 / p00 )
** ( ( xlam1 - ktan0 ) / ( xlam2 - ktan0 ) );
rapp0 = rapp0 * p00 / xp1 ;
'SI' ( p_test0 '>' xpc0 ) ;
trepp0 = ( ( treps0 / ktan0 ) + ( 'LOG' rapp0 ) ) ;
trepp0 = trepp0 * ( xlam2 - ktan0 ) * ktan0 / xlam2 ;
p_th0 = xp1 * ( 'EXP' ( -1. / ktan0 * ( treps0 - trepp0 ) ) ) ;
'SINON' ;
trepp0 = 0. ;
p_th0 = p_test0 ;
'FINSI' ;
trepp0 = -1. * trepp0 ;
trepp02 = trepp0 ;
*
trepp2 = trepp2 'ET' ( 'PROG' trepp0 ) ;
pres2 = pres2 'ET' ( 'PROG' p_th0 ) ;
*
* Calcul des erreurs
*
'SI' ( trepp01 '>' 1.E-6 ) ;
err_t0 = ( 'ABS' ( trepp01 - trepp02 ) ) / ( 'ABS' trepp01 ) ;
err_t0 = err_t0 * 100. ;
'SINON' ;
err_t0 = 0. ;
'FINSI' ;
err_p0 = ( 'ABS' ( pres01 - p_th0 ) ) / ( 'ABS' pres01 ) ;
err_p0 = err_p0 * 100. ;
'SI' ecri0 ;
'MESS' ' Erreur sur la deformation en % ' err_t0 ;
'MESS' ' Erreur sur la pression en % ' err_p0 ;
'FINSI' ;
err_t1 = err_t1 'ET' ( 'PROG' err_t0 ) ;
err_p1 = err_p1 'ET' ( 'PROG' err_p0 ) ;
*
'SI' ((err_t0 '>' xlim0) 'OU' (err_p0 '>' xlim0)) ;
pb0 = vrai ;
'FINSI' ;
*
'FIN' bloc0 ;
*
* Traces
*
'SI' graph0 ;
evol1 = 'EVOL' 'MANU' 'trace' treps1 'pression' pres1 ;
evol2 = 'EVOL' 'MANU' 'trace' treps1 'PC0' pc1 ;
evol21 = 'EVOL' 'MANU' 'trace' treps1 'PCS' pcs1 ;
evol3 = 'EVOL' 'MANU' 'trace' trepp1 'pression' pres1 ;
evol4 = 'EVOL' 'MANU' 'trace' trepp2 'pression' pres2 ;
evol5 = 'EVOL' 'MANU' 'trace' treps1 'succion' succ1 ;
evol6 = 'EVOL' 'MANU' 'trace' treps1 'succion' slim1 ;
evol7 = 'EVOL' 'MANU' 'trace' treps1 'erreur' err_t1 ;
evol8 = 'EVOL' 'MANU' 'trace' treps1 'erreur' err_p1 ;
TAB = 'TABLE' ;
TAB. 1 = 'MARQ PLUS ' ;
TAB. 2 = 'MARQ CROI ' ;
TAB. 3 = 'MARQ TRIA ' ;
TAB.'TITRE' = 'TABLE' ;
TAB.'TITRE'. 1 = MOT 'PRESSION' ;
TAB.'TITRE'. 2 = MOT 'PC' ;
TAB.'TITRE'. 3 = MOT 'PCS' ;
'DESS' ( evol1 'ET' evol2 'ET' evol21 )
'TITR' ' Trace de la deformation, pression ' 'LEGE' TAB ;
TAB = 'TABLE' ;
TAB. 1 = 'MARQ PLUS ' ;
TAB. 2 = 'MARQ CROI ' ;
TAB.'TITRE'= 'TABLE' ;
TAB.'TITRE'. 1 = MOT 'CALCUL' ;
TAB.'TITRE'. 2 =MOT 'THEORIE' ;
'DESS' ( evol3 'ET' evol4 )
'TITR' ' Comparaison calcul-theorie ' 'LEGE' TAB ;
TAB = 'TABLE' ;
TAB. 1 = 'MARQ PLUS' ;
TAB. 2 = 'MARQ CROI' ;
TAB.'TITRE'='TABLE' ;
TAB.'TITRE'. 1 =MOT 'SUCCION' ;
TAB.'TITRE'. 2 = MOT'LIMITE' ;
'DESS' ( evol5 'ET' evol6 )
'TITR' ' Succion ' 'LEGE' TAB ;
TAB = 'TABLE' ;
TAB. 1 = 'MARQ PLUS' ;
TAB. 2 = 'MARQ CROI' ;
TAB.'TITRE' = 'TABLE' ;
TAB.'TITRE'. 1 = MOT 'DEFORMATION' ;
TAB.'TITRE'. 2 = MOT 'PRESSION' ;
'DESS' ( evol7 'ET' evol8 )
'TITR' ' Erreur en % ' 'LEGE' TAB ;
'FINSI' ;
*
'FINSI' ;
*
'SI' pb0 ;
'ERRE' 5 ;
'SINON' ;
'ERRE' 0 ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer