Accueil | Cast3M

* 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 ;
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 ;  
ta1 .'HYPOTHESE_DEFORMATION'   = 'LINEAIRE' ;
* 
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' ;