Télécharger pakzad.eso

Retour à la liste

pakzad
C PAKZAD    SOURCE    OF166741  25/11/04    21:15:58     12349          
       SUBROUTINE PAKZAD(DEPST,NSTRS,NCOMAT,NVARI,
     .            MFR1,IB,IGAU,
     .            XMAT,SIG0,VAR0,SIGF,VARF,DEFP,KERRE,DSIGT,
     .            SUCC1,SUCC2)
 
***************************************************************
*
*  Modèle d'argile de PAKZAD
*
***************************************************************
*
*_________________________________________________________________
*
*    ENTREES :
*    ---------
*
* DEPST = INCREMENT DE DEFORMATIONS TOTALES
* NSTRS = NBRE DE COMPOSANTES DES DEFORMATIONS
* NCOMAT= NBRE DE CARACTERISTIQUES MECANIQUES DU MATERIAU
* NVARI = NBRE DE VARIABLES INTERNES
* MFR1  = NUMERO DE LA FORMULATION
* IB    = NUMERO DE L ELEMENT COURANT
* IGAU  = NUMERO DU POINT COURANT
* SIG0(NSTRS) = CONTRAINTES AU DEBUT DU PAS D'INTEGRATION
* VAR0(NVARI) = VARIABLES INTERNES AU DEBUT DU PAS DE TEMPS
* XMAT(NCOMAT) = CARACTERISTIQUES MECANIQUES DU MATERIAU
* SUCC1  SUCCION AU DEBUT DU PAS
* SUCC2  SUCCION A LA FIN DU PAS
*
*    SORTIE :
*    --------
*
* SIGF(NSTRS) = CONTRAINTES FINALES
* VARF(NVARI) = VARIALES INTERNES A LA FIN DU PAS D'INTEGRATION
* DEFP(NSTRS) = INCREMENT DE DEFORMATION PLASTIQUE A LA FIN
*               DU PAS D'INTEGRATION
* ============================================================
*       ICI IL FAUT PROGRAMMER EN FORTRAN PUR
*=============================================================
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
 
      DIMENSION SIG0(*),DEPST(*),VAR0(*),XMAT(*),SIGF(*),
     &          VARF(*),DEFP(*),DSIGT(*)
      DIMENSION RSIG0(6),RDEPS0(6),RSIGF(6),RDEFP(6),RSIG1(6)
      DIMENSION DEVT(6),SIGT(6),DEFT(6)
      REAL*8 LS0,LS1
*
*  Adaptation de l'option de calcul vers le 3D massif de SIG0 a RSIG0
*====================================================================
*
      IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31)  THEN
*
*---> 1 formulation massive
*---> 2 formulation quasi incompressible
*---> MASSIF 3D
*
        IF (NSTRS .EQ. 6) THEN
          DO 110 I=1,NSTRS
            RSIG0(I)=SIG0(I)
            RDEPS0(I)=DEPST(I)
 110      CONTINUE
        ELSE IF ( NSTRS .EQ. 4 .AND. ((IFOUR .EQ. 0)
     &    .OR.(IFOUR .EQ. -1))) THEN
*
*---> Calcul en mode deformations planes ou axisymetrique
*
          DO 115 I=1,NSTRS
            RSIG0(I)=SIG0(I)
            RDEPS0(I)=DEPST(I)
 115      CONTINUE
          RSIG0(5)=0.D0
          RSIG0(6)=0.D0
          RDEPS0(5)=0.D0
          RDEPS0(6)=0.D0
        ENDIF
      ELSE
        KERRE = 99
        RETURN
      ENDIF
 
*   Passage des deformations de cisaillement exprimées
*   en GAMA aux déformations de cisaillement exprimées
*   en déformations
      DO 116 I=4,6
         RDEPS0(I)=RDEPS0(I)*0.5D0
 116  CONTINUE
 
*   Données du materiau
*===========================================================
*
      YOUN0=XMAT(1)
      XNU0=XMAT(2)
      XN0=XMAT(5)
      XKA0=XMAT(6)
      XGA0=XMAT(7)
      XPA0=XMAT(8)
      XPC0=XMAT(9)
      XM0=XMAT(10)
      XBET1=XMAT(11)
      XA0=XMAT(12)
      XPCR0=XMAT(13)
      XS0=XMAT(14)
      XM1=XMAT(15)
      XM2=XMAT(16)
      XM3=XMAT(17)
      XBET2=XMAT(18)
      XTAU0=XMAT(19)
 
*   Quelques initialisations
*=============================================================
*
*---> Succion
*  . Succion initiale
*  . succion finale
*  . Succion de désaturation
*
******************************
*   MODIFICATION  POUR LA SUCCION
*
      IF(SUCC1.LT.-1.E34.AND.SUCC2.LT.-1.E34)  THEN
*
*       cas à succion constante : on la prend dans var0
*
         SUCI0=VAR0(3)
         SUCF0=SUCI0
      ELSE
         SUCI0=SUCC1
         SUCF0=SUCC2
      ENDIF
*
      IF (VAR0(2).LE.1.D-10) VAR0(2)=XS0+XM1*(XPC0-XS0)
      SDT=VAR0(2)
*
*---> Pression de consolidation initiale dans le domaine saturé
*
      IF (VAR0(4).LE.1.D-10) VAR0(4)=XPC0
*
*---> On est en train d'activer la surface de contact ou de gonflement
*
      IF ((VAR0(21).LE.1.D-10).AND.(VAR0(22).LE.1.D-10)) THEN
*
*  On active la surface de contact
*
        IGONF=1
      ELSE
*
*  On active la surface de gonflement
*
        IGONF=0
      ENDIF
*
*---> On initialise les variables internes à leur valeur initiale
*
      DO 49 I=1,NVARI
        VARF(I)=VAR0(I)
 49   CONTINUE
*
*---> On initialise les déformations plastiques au cours du pas à
*     zéro
*
      DO 48 I=1,6
        RDEFP(I)=0.D0
 48   CONTINUE
 
*    Sous incrémentations : définition du nombre de sous
*    incréments
*=============================================================
*
      nmax0=1
      iter2=0
      crit0=XKA0*(RDEPS0(1)+RDEPS0(2)+RDEPS0(3))
      crit0=ABS(crit0)
      IF (crit0.GT.1.D-2) nmax0=NINT(100.D0*crit0)+1
      IF (nmax0.GT.10000) nmax0=10000
      DO 50 I=1,6
        RDEPS0(I)=RDEPS0(I)/nmax0
 50   CONTINUE
*
*   Début de la boucle de sous incrémentation
*=============================================================
*
 51   iter2=iter2+1
*
*---> Succion au début (SUC1) et à la fin (SUC2) du sous incrémént
*
      SUC1=SUCI0+(SUCF0-SUCI0)*(iter2-1)/nmax0
      SUC2=SUCI0+(SUCF0-SUCI0)*iter2/nmax0
*
*---> Variable de test de la plasticité
*   . iplas vaut 1 si on plastifie sur la surface de contact
*
      IPLAS=0
*
*---> Pression d'origine capillaire au début du sous incrément
*
      LS0=SUC1-XM1*VARF(4)+(XM1-1.D0)*XS0
      IF (LS0.LE.0.D0) THEN
*
*  Etat quasi saturé
*
        PSAT0=SUC1
      ELSE
*
*  Etat partiellement saturé
*
        XK1=(1.D0+(SUC1-SDT)/((XM2-1.D0)*SDT))**2.D0
        XK1=1/XK1
        PSAT0=SDT+(XK1**(.5D0))*(SUC1-SDT)
      ENDIF
*
*---> Contraintes effectives au début du sous incrément
*
      DO 9 I=1,6
        A=0.D0
        IF (I.LE.3) A=1.D0
        RSIG1(I)=RSIG0(I)-PSAT0*A
 9    CONTINUE
      PRES00=-1.D0/3.D0*(RSIG1(1)+RSIG1(2)+RSIG1(3))
*
*---> Déformations élastiques test et déformation plastique
*
      DO 10 I=1,6
        DEFT(I)=0.D0
 10   CONTINUE
      TREPEL0=RDEPS0(1)+RDEPS0(2)+RDEPS0(3)
*
*---> Contraintes test: 2 cas suivant la valeur de XN0
*
      dif0=ABS(XN0-1.D0)
      IF (dif0.LE.1.D-5) THEN
*
*  Cas où XN0=1
*
        PRES1=PRES00*EXP(-1.D0*XKA0*TREPEL0)
        DO 11 I=1,6
          A=0.D0
          IF (I.LE.3) A=1.D0
          DEVT(I)=RSIG1(I)+PRES00*A
          DEVT(I)=DEVT(I)+2.D0*XGA0*PRES00*(RDEPS0(I)-TREPEL0/3.D0*A)
          SIGT(I)=DEVT(I)-PRES1*A
 11     CONTINUE
      ELSE
*
*  Autres cas
*
        XKA1=XKA0*XPA0*((PRES00/XPA0)**XN0)
        XGA1=XGA0*XPA0*((PRES00/XPA0)**XN0)
        PRES1=PRES00-(XKA1*TREPEL0)
        DO 21 I=1,6
          A=0.D0
          IF (I.LE.3) A=1.D0
          DEVT(I)=RSIG1(I)+PRES00*A
          DEVT(I)=DEVT(I)+2.D0*XGA1*(RDEPS0(I)-TREPEL0/3.D0*A)
          SIGT(I)=DEVT(I)-PRES1*A
 21     CONTINUE
      ENDIF
*
*---> Pressions de consolidation test à la succion SUC2 et module plastique BET1
* à la fin du sous incrément
*    PCT1 : pression de consolidation
*    PCT0 : pression de consolidation "fictive" dans le domaine
*           saturé ( = VARF(4))
*
      LS1=SUC2-XM1*VARF(4)+(XM1-1.D0)*XS0
      PCT0=VARF(4)
      IF (LS1.LE.0.D0) THEN
*
*  Etat quasi saturé
*
        PCT1=PCT0
        BET1=XBET1
      ELSE
*
*  Etat partiellement saturé
*
        XK2=(1.D0+(SUC2-SDT)/(XM3*SDT))**2.D0
        XK2=1/XK2
        PCT1=PCT0+(XK2**.5D0)*(SUC2-SDT)
        BET1=XBET2-((XBET2-XBET1)*EXP(-1.D0*XTAU0*(SUC2-SDT)))
*
*  Dans ce qui suit, on va négliger les variations de BET1 avec SDT
*  (et donc avec l'écoulement plastique)
*
      ENDIF
*
*---> Si on est en recharge, on vérifie la condition de gonflement avant
*     la condition de contact
*     A condition de ne pas déja avoir active la surface de charge
*     (IGONF=0)
*
        IF (((RDEPS0(1)+RDEPS0(2)+RDEPS0(3)).LE.0.D0)
     .      .AND.(IGONF.EQ.0)) THEN
          CALL PAKGON(XMAT,SIGT,VARF,SIGT,DEFT,VARF,
     .                RDEPS0,KERRE,IB,IGAU,PCT0,
     .                SUC2,RSIG1)
          PRES1=-1.D0/3.D0*(SIGT(1)+SIGT(2)+SIGT(3))
          DO 52 I=1,6
             A=0.D0
             IF (I.LE.3) A=1.D0
             DEVT(I)=SIGT(I)+PRES1*A
 52       CONTINUE
          IGONF=0
        ENDIF
*
*---> Contrainte équivalente
*
      STEST2=1.5D0*PROCON(DEVT,DEVT,6)
      IF (STEST2.LE.(1.D-10*YOUN0)) STEST2=0.D0
*
*---> Fonction de charge
*
      PHIT=Stest2+XM0*XM0*PRES1*(PRES1-PCT1)
*
*  Vérification du critère de plasticité
*=========================================================
*
*---> Erreur admise
*
      ERR0=1.D-7*ABS(PHIT)
      ERR0=MAX(ERR0,1.D-10*YOUN0)
*
*---> Critère de plasticité
*
      PHI0=PHIT
      iter0=0
*
 99   IF (PHI0.LE.ERR0) THEN
*
*  On est élastique par rapport à la surface de contact
*=========================================================
*
*
*---> On vérifie la condition de gonflement
*
      IF (IPLAS.EQ.0) THEN
*
*  On a pas plastifié sur la surface de contact
*  Si on est en décharge, on vérifie la condition de gonflement
*   après la condition de contact
*  Cas particulier, le premier pas du calcul, IGONF peut valoir 1 car
*   toutes les valeurs de des variables internes ne sont peut etre pas
*   initialisées de facon adéquate.
*   Or on doit vérifier la surface de gonflement systématiquement
*   au premier pas du calcul et ce jusqu'a ce qu'on active
*   la surface de gonflement ou la surface de contact
*
*
        IF (((RDEPS0(1)+RDEPS0(2)+RDEPS0(3)).GT.0.D0).OR.
     .       (IGONF.EQ.1)) THEN
          CALL PAKGON(XMAT,SIGT,VARF,SIGT,DEFT,VARF,
     .                RDEPS0,KERRE,IB,IGAU,PCT0,
     .                SUC2,RSIG1)
          IGONF=0
        ENDIF
        DO 53 I=1,6
          RDEFP(I)=RDEFP(I)+DEFT(I)
          RSIGF(I)=SIGT(I)
 53     CONTINUE
      ELSE
*
*  On a plastifié sur la surface de contact, on ré-initialise
*  les variables internes du gonflement
*    . le rayon est remis à sa valeur initiale
*    . le "Centre" de l'ellipse (P1,Q1) est mis à
*      la valeur courante des contraintes
*
       VARF(5)=XA0
       VARF(6)=XA0
       VARF(7)=PRES1
       VARF(8)=PRES1
       VARF(21)=0.D0
       VARF(22)=0.D0
       DO 399 I=1,6
         RDEFP(I)=RDEFP(I)+DEFT(I)
         RSIGF(I)=SIGT(I)
         VARF(8+I)=DEVT(I)
         VARF(14+I)=DEVT(I)
 399  CONTINUE
      ENDIF
*
*---> A la fin du sous incrémént: Est on saturé ou non ?
*
      LS1=SUC2-XM1*PCT0+(XM1-1.D0)*XS0
*
*---> Mise à jour de la succion de désaturation à la fin di sous incrément
*
      SDT=XS0+XM1*(PCT0-XS0)
*
*---> Pression d'origine capillaire à la fin de sous incrément
*
      IF (LS1.LE.0.D0) THEN
*
*  Etat quasi saturé
*
        PSAT0=SUC2
      ELSE
*
*  Etat partiellement saturé
*
        XK1=(1.D0+(SUC2-SDT)/((XM2-1.D0)*SDT))**2.D0
        XK1=1/XK1
        PSAT0=SDT+(XK1**(.5D0))*(SUC2-SDT)
      ENDIF
*
*---> Grandeurs à la fin de sous incrément
*
      DO 400 I=1,6
         A=0.D0
         IF (I.LE.3) A=1.D0
         RSIGF(I)=RSIGF(I)+PSAT0*A
         RSIG0(I)=RSIGF(I)
         DEFT(I)=RDEFP(I)-(RDEFP(1)+RDEFP(2)+RDEFP(3))/3.D0*A
 400  CONTINUE
      EPSPT=PROCON(DEFT,DEFT,6)
      IF (EPSPT.LE.1.D-10) EPSPT=0.D0
      EPSPT=(2.D0/3.D0*EPSPT)**.5D0
      VARF(1)=VAR0(1)+EPSPT
      VARF(2)=SDT
      VARF(3)=SUCF0
      VARF(4)=PCT0
*
      ELSE
*
*  On plastifie
*=========================================================
*
*---> Mise du flag IPLAS à 1: on plastifie en contact
*     Mise du flag IGONF à 1: on plastifie en contact
*
      IPLAS=1
      IGONF=1
*
*---> Condition de consistance: calcul du paramètre plastique
*
      IF (dif0.LE.1.D-5) THEN
*
*  Cas où XN0=1
*
        DF1=XM0*XM0*(2.D0*PRES1-PCT1)
        DF2=XM0*XM0*PCT1
        DF0=12.D0*XGA0*PRES00*Stest2
        denom0=DF0+XKA0*DF1*DF1*PRES1+PRES1*BET1*DF1*DF2
*
      ELSE
*
*  Autres cas
*
        DF1=XM0*XM0*(2.D0*PRES1-PCT1)
        DF2=XM0*XM0*PCT1
        DF0=12.D0*XGA1*Stest2
        denom0=DF0+XKA1*DF1*DF1+PRES1*BET1*DF1*DF2
      ENDIF
*
      dlam0=PHIT/denom0
*
*---> Mise à jour des contraintes
*
      IF (dif0.LE.1.D-5) THEN
*
*  Cas où XN0=1
*
        PRES1=PRES1*EXP(-1.D0*XKA0*DF1*dlam0)
        DO 12 I=1,6
          A=0.D0
          IF (I.LE.3) A=1.D0
          DEVT(I)=DEVT(I)-6.D0*XGA0*PRES00*DEVT(I)*dlam0
          SIGT(I)=DEVT(I)-PRES1*A
 12     CONTINUE
      ELSE
*
*  Autres cas
*
        PRES1=PRES1-(XKA1*DF1*dlam0)
        DO 22 I=1,6
          A=0.D0
          IF (I.LE.3) A=1.D0
          DEVT(I)=DEVT(I)-6.D0*XGA1*DEVT(I)*dlam0
          SIGT(I)=DEVT(I)-PRES1*A
 22     CONTINUE
      ENDIF
*
*---> Mise à jour de la pression de consolidation à la succion finale
*
      PCT1=PCT1*EXP(BET1*DF1*dlam0)
      PCT0=PCT0*EXP(XBET1*DF1*dlam0)
*
*---> Déformation élastiques et plastiques
*
      IF (dif0.LE.1.D-5) THEN
*
*  Cas où XN0=1
*
      TREPEL0=-1.D0/XKA0*LOG(PRES1/PRES00)
      DO 13 I=1,6
        A=0.D0
        IF (I.LE.3) A=1.D0
        DEFT(I)=(DEVT(I)-RSIG1(I)-PRES00*A)/(2.D0*XGA0*PRES00)
        DEFT(I)=RDEPS0(I)-DEFT(I)-TREPEL0*A/3.D0
 13   CONTINUE
      ELSE
*
*  Autres cas
*
      TREPEL0=-1.D0/XKA1*(PRES1-PRES00)
      DO 23 I=1,6
        A=0.D0
        IF (I.LE.3) A=1.D0
        DEFT(I)=(DEVT(I)-RSIG1(I)-PRES00*A)/(2.D0*XGA1)
        DEFT(I)=RDEPS0(I)-DEFT(I)-TREPEL0*A/3.D0
 23   CONTINUE
      ENDIF
*
*---> Mise à jour de la contrainte équivalente
*
      STEST2=1.5D0*PROCON(DEVT,DEVT,6)
      IF (STEST2.LE.(1.D-10*YOUN0)) STEST2=0.D0
*
*---> Mise à jour de la fonction de charge
*
      PHIT=Stest2+XM0*XM0*PRES1*(PRES1-PCT1)
      PHI0=ABS(PHIT)
*
*---> Test de convergence ou itération suivante
*
      iter0=iter0+1
      IF (iter0.LT.200) THEN
         GOTO 99
      ELSE
         PHI0=0.D0
         KERRE=460
         GOTO 99
      ENDIF
*
      ENDIF
*
*   Vérification des sous incrémentations
*========================================================================
*
      IF ((iter2.LT.nmax0).AND.(KERRE.EQ.0)) THEN
         GOTO 51
      ENDIF
*
*===========================================================================
*
*   Passage des deformations de cisaillement exprimées
*   en déformations aux déformations de cisaillement exprimées
*   en GAMA
      DO 117 I=4,6
         RDEFP(I)=RDEFP(I)*2.D0
 117  CONTINUE
*
*  Passage a l'option de calcul pour les contraintes
*  Grandeurs finales
*=========================================================
*
      IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31)  THEN
*---> MASSIF 3D
*---> Calcul axisymétrique ou contraintes planes
          DO 170 I=1,NSTRS
            SIGF(I)=RSIGF(I)
            DEFP(I)=RDEFP(I)
            DSIGT(I)=RSIGF(I)-RSIG0(I)
 170      CONTINUE
      ENDIF
*
      RETURN
      END
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales