Sources Esope | Cast3M

huitch
C HUITCH    SOURCE    CHAT      10/12/16    21:15:59     6819
C=======================================================================
C=        MODELE HYPERELASTIQUE 8Chaines                               =
C=                 EN GRANDES TRANSFORMATIONS                          =
C=            CONTRAINTES PLANES   Modele INCOMPRESSIBLE (direction 3) =
C=            DEFORMATIONS PLANES  Modele Quasi INCOMPRESSIBLE         =
C=            TRIDIMENSIONNEL      Modele Quasi INCOMPRESSIBLE         =
C=  Contribution de Laurent Gornet - Ecole Centrale de Nantes (2009)   =
C=======================================================================
C=                                                                     =
C=  Exemple d'utilisation d'un modele UMAT en grandes transformations  =
C=                                                                     =
C=  Pour plus d'informations, voir la presentation de L. Gornet lors   =
C=  du Club Cast3m 2009, disponible sur le site Web de Cast3m.         =
C=                                                                     =
C= Note : Actuellement en grandes deformations dans PASAPAS, le modele =
C=        ne peut contenir que des modèles type UMAT. On ne peut pas   =
C=        "melanger" les derivees objectives dans Cast3m.              =
C=======================================================================
      SUBROUTINE HUITCH (STRESS, STATEV, DDSDDE, STRAN, DSTRAN,
     &                   TIME, DTIME, TEMP, DTEMP, PREDEF, DPRED,
     &                   NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS,
     &                   DFGRD0, DFGRD1, KSTEP, KINC,SSE,spd,scd,rpl,
     &                   DDSDDT,DRPLDE,DRPLDT,cmname,coords,drot,PNEWDT,
     &                   celent,NOEL, NPT, LAYER, KSPT )
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER*16  CMNAME
 
      INTEGER  NDI, NSHR, NTENS, NSTATV, NPROPS, KSTEP, KINC
     &         NOEL, NPT, LAYER, KSPT
      REAL*8   STRESS(NTENS), STATEV(NSTATV),SSE, SPD, SCD, RPL,
     &         DDSDDE(NTENS,NTENS), STRAN(NTENS), DSTRAN(NTENS),
     &         DDSDDT(NTENS), DRPLDE(NTENS), DRPLDT,
     &         TIME(2), DTIME, TEMP, DTEMP, PREDEF(*), DPRED(*),
     &         PROPS(NPROPS), DFGRD0(3,3), DFGRD1(3,3),
     &         COORDS(3),DROT(3,3),PNEWDT,CELENT
      REAL*8   CG11,CG22,CG12,CG33, CInv11,CInv22,CInv12, DLT,
     &         I1,I2, dWdI1,dWdI2, phyd, Coe1,Coe2,Coe3
     &         F11,F21,F12,F22, S11,S12,S21,S22
     &         F13,F31,F23,F32,F33
     &         CGg1,CGg2,CGg3,CGg4,CGg5,CGg6
     &         C2Gg1,C2Gg2,C2Gg3,C2Gg4,C2Gg5,C2Gg6
     &         dWisodI1bar,dWisodI2bar,eg1,eg2,Gam1,Gam2
      REAL*8 Cc11,Cc22,Cc12,Cc33
      REAL*8  I1B,I2B,cinv33,detF,mdetF,Jm2s3,Jm1s3
      REAL*8 S33,dwvdj,zz,C1,C2
C===============================================
C 1 DECLARER VOTRE MATERIAU 8Chaines :
      REAL*8 nkT,VN,klm,lmp,mpo
C
      REAL*8 zero,one,two,three,four,six
      Parameter(zero=0.d0,one=1.d0,two=2.d0,three=3.d0,
     & four=4.d0,six=6.d0)
      PARAMETER (cm5s3=-1.66666666666666666666666666666666666)
      PARAMETER (cm7s3=-2.33333333333333333333333333333333333)
      PARAMETER (cm2s3=-0.66666666666666666666666666666666666)
      PARAMETER (cm1s3=-0.33333333333333333333333333333333333)
      PARAMETER (cm4s3=-1.33333333333333333333333333333333333)
c      print*,'GORNET', ndi
      if(ndi.ne.-2) go to 100
* formulation 2D contraintes planes
C  PARAMETRES MATERIAU DU MODELE 8Chaines :
C ===============================================
C= Dans le cas du modele
C= la densite d'energie de deformation est definie par :
C=    W = f(I1,I2) I3= 1
C ===============================================
C  3 PARAMETRES  MATERIAU DU MODELE HARTSMITH :
C ===============================================
      nkt= PROPS(3)
      vn = PROPS(4)
C Deformation totale stockage perso dans var 4 5 6
C ******************** GREEN LAGRANGE **************
C
       F11 = dfgrd1(1,1)
       F21 = dfgrd1(2,1)
       F12 = dfgrd1(1,2)
       F22 = dfgrd1(2,2)
C*******************
C       print*,F11
C       print*,F21
C       print*,F12
C       print*,F22
C
C ******************** CC = FT * F
C
        cc11 =  F11**2 + F21**2
        cc22 =  F12**2 + F22**2
        cc12 =  (F11*F12)+(F21*F22)
C
C        print*,cc11
C        print*,cc22
C        print*,cc22
C ******************** CC-1
           DLT = CC11*CC22 - CC12*CC12
               cinv11 =  1.0D0/DLT*CC22
           cinv22 =  1.0D0/DLT*CC11
           cinv12 = -1.0D0/DLT*CC12
C
C***********************************************************
C det J =1  ==>c
            cc33=1.0D0/(cc11*cc22-cc12**2)
C trace C
            i1= cc11+cc22+cc33
            i2= .50D0*(i1**2-cc11**2-cc22**2-cc33**2-2*cc12**2)
C===============================================================
C 2 INTEGRER VOTRE MATERIAU 8Chaines : DWDI1 = ?, DWDI2 = ?
C Constantes materielles
C              W = f(i1,i2)
C234567
      klm = 33.D0/(1050.D0*(vN**2))*(I1**2)
      lmp = 19.D0*4.0D0/(7000.D0*(vN**3))*(I1**3)
      mpo = 519.D0*5.0D0/(673750.D0*(VN**4))*(I1**4)
C
      dWdI1 = nkT*(0.5D0+(0.1D0/vN)*I1+klm+lmp+mpo)
      dWdI2 = 0.d0
C===============================================================
C p : pression hydrostatique obtenu ici a partir de S33 = 0.
            p=2.0D0*((dwdi1+i1*dwdi2)-dwdi2*cc33)*cc33
C******************** CONTRAINTES PK2    ******************
C
        S11=  (2.0D0*((dwdi1+i1*dwdi2)-dwdi2*cc11)-p*cinv11)
        S22=  (2.0D0*((dwdi1+i1*dwdi2)-dwdi2*cc22)-p*cinv22)
        S12=  (2.0D0*(-dwdi2*cc12)-p*cinv12)
C
C
C******************** CONTRAINTES Cauchy     ******************
C A FAIRE PK2 a transformer par  cauchy  = j-1  F S FT
C ici incompressible j =1 a modifier par F et FT,
C***************************************************************
C234567
      stress(1)=(S11*F11**2)+(2.0D0*F11*F12*S12)+(F12**2*S22)
      stress(2)=(S11*F21**2)+(2.0D0*F21*F22*S12)+(F22**2*S22)
      stress(4)=(F21*F11*S11)+(F21*F12*S12)+(F22*F11*S12)
     . +(F22*F12*S22)
      RETURN
 100  CONTINUE
*  formulation 3D massive
      if(ndi.ne.2)  go to 200
C  PARAMETRES MATERIAU :
C =======================
C= Dans le cas du modele Mooney Rivlin compressible,
C= l'energie de deformation est decomposee en deux termes decouples :
C= - la partie isochorique ou incompressible Wiso, fonction des inva-
C=   riants du tenseur de Cauchy-Green modifie ;
C= - la partie purement volumique Wvol, dependant de J=det(F).
C= Pour le present modele, nous avons :
C= - Wiso = Wiso(I1bar,I2bar) = Coe2 * (I1bar-3.) + Coe3 * (I2bar-3.)
C= - Wvol = Wvol(J) = 1/Coe1 * (J-1)*(J-1)
C= Coe2 et Coe3 : coefficients du materiau
C= le module de cisaillement est egal a mu = 2*Coe2
C= bbar : tenseur de Cauchy-Green gauche modifie
C= Par definition, bbar = J**(-2/3)*(F.Ft) = J**(-2/3)*CGg
C= I1bar : 1er invariant de bbar (= trace(bbar))
C= Dans le cas de la quasi-incompressibilite, c.a.d. J proche de 1,
C= Wvol peut etre interpretee comme une fonction de penalisation.
C*    Youn = PROPS(1)
C*    Pois = PROPS(2)
c      Coe1 = PROPS(5)
c      Coe2 = PROPS(3)
c      Coe3 = PROPS(4)
C ===============================================
C  2 PARAMETRES  MATERIAU DU MODELE 8Chaines :
C ===============================================
      nkt= PROPS(3)
      vn = PROPS(4)
      Coe1 = PROPS(5)
C  GRADIENT DE LA TRANSFORMATION (FIN DU PAS) :
C ==============================================
      F11 = DFGRD1(1,1)
      F12 = DFGRD1(1,2)
      F13 = DFGRD1(1,3)
      F21 = DFGRD1(2,1)
      F22 = DFGRD1(2,2)
      F23 = DFGRD1(2,3)
      F31 = DFGRD1(3,1)
      F32 = DFGRD1(3,2)
      F33 = DFGRD1(3,3)
 
C  JACOBIEN DE LA TRANSFORMATION (FIN DU PAS) :
C ==============================================
      detF =  F11*F22*F33 - F12*F21*F33 + F12*F23*F31
     &      + F13*F32*F21 - F13*F31*F22 - F23*F32*F11
 
C  TENSEUR DES DEFORMATIONS DE CAUCHY-GREEN GAUCHE
C =================================================
C= Tenseur de Cauchy-Green gauche CGg = F.Ft
C G11 G22 G33  G12 G13 G23
      CGg1 = F11*F11 + F12*F12 + F13*F13
      CGg2 = F21*F21 + F22*F22 + F23*F23
      CGg3 = F33*F33 + F31*F31 + F32*F32
      CGg4 = F11*F21 + F12*F22 + F13*F23
      CGg5 = F11*F31 + F12*F32 + F13*F33
      CGg6 = F21*F31 + F22*F32 + F23*F33
C= Tenseur de Cauchy-Green gauche CGg2 = (F.Ft)*(F.Ft)
C indices (1 a 6) =(11 22 33  12 13 23
      C2Gg1 = CGg1*CGg1 + CGg4*CGg4 + CGg5*CGg5
      C2Gg2 = CGg4*CGg4 + CGg2*CGg2 + CGg6*CGg6
      C2Gg3 = CGg5*CGg5 + CGg6*CGg6 + CGg3*CGg3
      C2Gg4 = CGg1*CGg4  + CGg4*CGg2 + CGg5*CGg6
      C2Gg5 = CGg1*CGg5 + CGg4*CGg6 + CGg5*CGg3
      C2Gg6 = CGg4*CGg5 + CGg2*CGg6 + CGg6*CGg3
C INVARIANTS : I1B= J**-1/3 I1, I2B= J**-4/3 I2
      I1=(CGg1+CGg2+CGg3)
      i2= .5D0*(i1**2-C2Gg1-C2Gg2-C2Gg3-2.D0*(C2Gg4+C2Gg5+C2Gg6))
      I1B = (detF**cm2s3)*I1
      I2B = (detF**cm4s3)*I2
 
C
            I1B = (detF**cm2s3)*I1
            I2B = (detF**cm4s3)*I2
C  CONTRAINTES DE CAUCHY (FIN DU PAS) :
C ======================================
C= Les contraintes de Cauchy STRESS se decomposent en deux termes :
C= - STRESS = SCvol + SCiso
C= - SCvol = phyd.Iden = (dWvol(J)/dJ).Iden
C=   avec Iden = tenseur identite d'ordre 2,
C=        phyd = pression hydrostatique,
C=   avec Gam1 =   (dWiso/dI1bar+I1bar.dWiso/dI2bar)
C=        Gam2 =  (dWiso/dI2bar)
C= - SCiso =  2. J**(-5/3)Gam1 CGg - 2. J**(-7/3) Gam2 C2Gg
C
C ENTRER LE MATERIAU ICI :
C===============================================================
C 2 INTEGRER VOTRE MATERIAU HARTSMITH : DWDI1 = ?, DWDI2 = ?
C Constantes materielles
C              W = f(i1,i2)
C
      klm = 33.D0/(1050.D0*(vN**2))*(I1**2)
      lmp = 19.D0*4.0D0/(7000.D0*(vN**3))*(I1**3)
      mpo = 519.D0*5.0D0/(673750.D0*(VN**4))*(I1**4)
C
      dWisodI1bar = nkT*(0.5D0+(0.1D0/vN)*I1+klm+lmp+mpo)
      dWisodI2bar = 0.d0
 
      Gam1 =  (dWisodI1bar+I1B*dWisodI2bar)
      Gam2 =  (dWisodI2bar)
 
      phyd = 2.0D0 * (detF-1.0D0) / Coe1
 
      eg1   = 2.0D0 * (detF**cm5s3) * Gam1
      eg2   = 2.0D0 * (detF**cm7s3) * Gam2
 
 
      STRESS(1) = eg1 * CGg1 - eg2 * C2Gg1 + phyd
      STRESS(2) = eg1 * CGg2 - eg2 * C2Gg2 + phyd
      STRESS(3) = eg1 * CGg3 - eg2 * C2Gg3 + phyd
      STRESS(4) = eg1 * CGg4 - eg2 * C2Gg4
      STRESS(5) = eg1 * CGg5 - eg2 * C2Gg5
      STRESS(6) = eg1 * CGg6 - eg2 * C2Gg6
      return
 200   if(ndi.ne.-1) go to 300
C Deformations Planes
C
      nkt= PROPS(3)
      vn = PROPS(4)
      Coe1 = PROPS(5)
C Deformation totale stockage perso dans var 4 5 6
C ******************** GREEN LAGRANGE **************
C
       F11 = dfgrd1(1,1)
       F21 = dfgrd1(2,1)
       F12 = dfgrd1(1,2)
       F22 = dfgrd1(2,2)
       F33 = dfgrd1(3,3)
       detF =F11*F22-F12*F21
C ******************** CC = FT * F
C
        cc11 =  F11**2 + F21**2
        cc22 =  F12**2 + F22**2
        cc12 =  (F11*F12)+(F21*F22)
C ******************** CC-1
           DLT = CC11*CC22 - CC12*CC12
               cinv11 =  1.0D0/DLT*CC22
           cinv22 =  1.0D0/DLT*CC11
           cinv12 = -1.0D0/DLT*CC12
           cinv33 =  1.0D0/DLT
C
C***********************************************************
C det J =1  ==>c
            cc33=1.0D0
C trace C  et CB
            i1= cc11+cc22+cc33
            i2= .50D0*(i1**2-cc11**2-cc22**2-cc33**2-2*cc12**2)
C
            I1B = (detF**cm2s3)*I1
            I2B = (detF**cm4s3)*I2
C===============================================================
C 2 INTEGRER VOTRE MATERIAU   : DWDI1B = ?, DWDI2B = ?
C attention nom !!!  dwdi1 =  DWDI1B et dwdi2 = DWDI2B
C              Wiso = f(i1b,i2b)
C
C===============================================================
C 2 INTEGRER VOTRE MATERIAU HARTSMITH : DWDI1 = ?, DWDI2 = ?
C Constantes materielles
C              W = f(i1,i2)
C
      klm = 33.D0/(1050.D0*(vN**2))*(I1**2)
      lmp = 19.D0*4.0D0/(7000.D0*(vN**3))*(I1**3)
      mpo = 519.D0*5.0D0/(673750.D0*(VN**4))*(I1**4)
C
      dWdI1 = nkT*(0.5D0+(0.1D0/vN)*I1+klm+lmp+mpo)
      dWdI2 = 0.d0
C===============================================================
                   dwvdj=  2.0D0 * (detF-1.0D0) / Coe1
C===============================================================
                   Jm2s3=detF**cm2s3
                   Jm1s3=detF**cm1s3
                   phyd = Jm1s3*dwvdj
C
C******************** CONTRAINTES PK2    ******************
C
C Epaisseur a ajouter dans les formules!!!
C234567
       S11=2.0D0*Jm2s3*((dwdi1+i1b*dwdi2)-dwdi2*cc11)+phyd*cinv11
       S22=2.0D0*Jm2s3*((dwdi1+i1b*dwdi2)-dwdi2*cc22)+phyd*cinv22
       S12=  2.0D0*Jm2s3*(-dwdi2*cc12)+phyd*cinv12
       S33=2.0D0*Jm2s3*((dwdi1+i1b*dwdi2)-dwdi2*cc33)+phyd*cinv33
C
C
C******************** CONTRAINTES Cauchy     ******************
C A FAIRE PK2 a transformer par  cauchy  = j-1  F S FT
C ici incompressible j =1 a modifier par F et FT,
C***************************************************************
      mdetF = 1.0D0/detF
C
      stress(1)=mdetF*((S11*F11**2)+(2.0D0*F11*F12*S12)+(F12**2*S22))
      stress(2)=mdetF*((S11*F21**2)+(2.0D0*F21*F22*S12)+(F22**2*S22))
      stress(4)=mdetF*((F21*F11*S11)+(F21*F12*S12)+(F22*F11*S12)
     . +(F22*F12*S22))
      stress(3)= mdetF*S33
      return
 300   kinc=-2
      return
      END