Télécharger damtg3d.eso

Retour à la liste

damtg3d
C DAMTG3D   SOURCE    CB215821  16/04/21    21:16:14     8920
      SUBROUTINE DAMTG3D(parahot3,idimpara3,H66t,rt1,rc1,
     .     esigma6v,lg,ntot,lerror,i1,i2,i3,i4,i5,i6,H66,
     .     lfulldamage,lcp)
 
c     SUBROUTINE added by THG
 
c     parahot3 is the array that contains the material parameters and state variables
c     idimpara3 is the number of parameters for parahot3
c     H66t is the tangent matrix
c     rt1 and rc1 are the vectors built in DPRAN3D and used in DAMTG3D for calculation of tangent matrix
c     esigma6v is the effective stress vector
c     lg and lerror are logical
c     ntot is the integration point (global number for the whole structure): counter
 
c     ==========================================================================================
c     !   DAMAGE BOX of the 3D plastic-damage concrete model by T. Gernay                      !
c     !   This subroutine calculates                                                           !
c     !                             the nominal stress at the end of the time step             !
c     !                             the nominal consistent algorithmic tangent modulus Ct      !
c     !                                                                                        !
c     !   The effective stress esigma6v and the eff. tangent modulus H66t given by plastic box !
C     !                                                                                        !
c     ==========================================================================================
 
      IMPLICIT REAL*8 (A-B,D-H,O-Z)
      implicit integer (I-K,M,N)
      implicit logical (L)
      implicit character*10 (C)
 
      dimension H66t(6,6)
c     Tangent matrix (output)
      dimension parahot3(idimpara3)
c     parahot3              : Material prop. and various info. at elevated temp. (input and output)
c     parahot3(idimpara3-34): equivalent plastic strain in compression at the end of the current step
c     parahot3(idimpara3-35): equivalent plastic strain in tension     at the end of the current step
c     parahot3(idimpara3-6:idimpara3-1) = effective stress obtained by dpran3D
      dimension esigma6v(6)
c     Effective stress (input)
 
      dimension rc1(1,6),rt1(1,6)
      dimension v1(3),v2(3),v3(3),rcossigmapr(3,3)
      dimension p1(3,3),p2(3,3),p3(3,3)
      dimension pt1(3,3,3,3),pt2(3,3,3,3),pt3(3,3,3,3)
      dimension projt(3,3,3,3),w4(3,3,3,3)
      dimension Qtens(3,3,3,3),Qcomp(3,3,3,3)
      dimension esigmat(6),esigmac(6)
      dimension sigma(6),projtm(6,6),projcm(6,6)
      dimension p12(3,3),p23(3,3),p13(3,3)
      dimension pp12(3,3,3,3),pp23(3,3,3,3),pp13(3,3,3,3)
      dimension W4m(6,6),H66temp(6,6)
      dimension sigrt(6,6),sigrc(6,6)
      dimension p121(3,3),p231(3,3),p131(3,3)
      dimension p122(3,3),p232(3,3),p132(3,3)
      dimension damtensor(6,6), damtensor2(6,6),rt1mat(3,3)
      dimension rc1mat(3,3),esigmatmat(3,3),esigmacmat(3,3)
      dimension sigrttens(3,3,3,3),sigrctens(3,3,3,3),H66(6,6)
      dimension Qtensm(6,6),Qcompm(6,6),rc1t(6)
      dimension U(6,6)
*     dimension id4(3,3,3,3)
 
c     Matrix of unity I, or U
c      -----------------------
 
      do jloc=1,6
        do iloc=1,6
          if (iloc.eq.jloc) then
            U(iloc,jloc) = 1.0d0
          else
            U(iloc,jloc) = 0.0d0
          endif
        enddo
      enddo
 
c     Fourth identity tensor id4
C     do iloc=1,3
C       do jloc=1,3
C         do kloc=1,3
C           do mloc=1,3
C             if ((iloc.eq.jloc).and.(kloc.eq.mloc)) then
C               id4(iloc,jloc,kloc,mloc) = 1.0d0
C             else
C               id4(iloc,jloc,kloc,mloc) = 0.0d0
C             endif
C           end do
C         end do
C       end do
C     end do
 
      i1 = 1
      i2 = 2
      i3 = 3
      i4 = 4
      i5 = 5
      i6 = 6
      i7 = 7
      i8 = 8
      i9 = 9
      i10=10
      i11=11
      i12=12
      i13=13
      i14=14
      i15=15
      i16=16
      i17=17
      i18=18
      i19=19
      i20=20
      i21=21
      i22=22
      i23=23
      i24=24
      i25=25
      i26=26
      i27=27
      i28=28
      i29=29
      i30=30
      i31=31
      i32=32
      i33=33
      i34=34
      i35=35
      i36=36
      i37=37
      i38=38
      i39=39
 
      r0 = 0.
      r1 = 1.
 
c     ================================================
c     !                                              !
c     !    EIGEN VALUES OF THE EFFECTIVE STRESSES    !
c     !    ======================================    !
c     ================================================
 
c     Effective stresses
      do iloc=i1,i6
        ESIGMA6V(iloc) = parahot3(idimpara3-i7+iloc)
      end do
 
 
c     Eigen values and eigen vector
      call sigmapr3ETC(esigma6v,S1,S2,S3,rcossigmapr,lcp)
      do iloc=i1,i3
        v1(iloc) = rcossigmapr(iloc,i1)
        v2(iloc) = rcossigmapr(iloc,i2)
        v3(iloc) = rcossigmapr(iloc,i3)
      end do
 
c     ===============================================
c     !                                             !
c     !           UPDATE DAMAGE VARIABLES           !
c     !           =======================           !
c     ===============================================
 
      rkappac = parahot3(idimpara3-34)
      rkappat = parahot3(idimpara3-35)
 
c     Tensile damage
      at = parahot3(16)
      dt = 1.0d0 - ((0.5d0)*exp(-at*rkappat)+(0.5d0)*
     .     exp(-6.0d0*at*rkappat))
      rdti = parahot3(idimpara3-41)
      if (dt.lt.rdti) then
        dt = rdti
      endif
      parahot3(idimpara3-33) = dt
 
 
c     Compressive damage
      rkappaic = parahot3(idimpara3-36)
      rdeltakc = rkappac - rkappaic
      rdci = parahot3(idimpara3-40)
      rdamprev = LOG(1.0d0-rdci)
      dc = 1.0d0 - exp((rdamprev)-parahot3(i14)*rdeltakc)
      if (dc.lt.parahot3(idimpara3-40)) dc = parahot3(idimpara3-40)
      parahot3(idimpara3-32) = dc
 
c     ===================================================================
c     !                                                                 !
c     !    PROJECTION OPERATORS OF THE EFFECTIVE STRESSES P+ AND P-     !
c     !    ========================================================     !
c     ===================================================================
 
c     p1(i3,i3) = v1 X v1 (tensorial product)
      call mulABT(v1,v1,p1,i3,i1,i3,i1)
c         [3;3]=[3;1]x[1;3]
      call mulABT(v2,v2,p2,i3,i1,i3,i1)
c         [3;3]=[3;1]x[1;3]
      call mulABT(v3,v3,p3,i3,i1,i3,i1)
c         [3;3]=[3;1]x[1;3]
 
c     P+ calculated as the sum = H(sig,i) p,i X p,i
c     First we calculate the PT,i = p,i X p,i
      do mloc=i1,i3
        do kloc=i1,i3
          do jloc=i1,i3
            do iloc=i1,i3
              pt1(iloc,jloc,kloc,mloc) = p1(iloc,jloc)*p1(kloc,mloc)
c             [3;3;3;3]        =   [3;3] x [3;3]    tensorial product
              pt2(iloc,jloc,kloc,mloc) = p2(iloc,jloc)*p2(kloc,mloc)
              pt3(iloc,jloc,kloc,mloc) = p3(iloc,jloc)*p3(kloc,mloc)
            end do
          end do
        end do
      end do
 
      if (S1.gt.r0) then
        h1 = r1
      else
        h1 = r0
      endif
      if (S2.gt.r0) then
        h2 = r1
      else
        h2 = r0
      endif
      if (S3.gt.r0) then
        h3 = r1
      else
        h3 = r0
      endif
 
c     Projection operators of the effective stresses P+ [3;3;3;3]
      do mloc=i1,i3
       do kloc=i1,i3
          do jloc=i1,i3
            do iloc=i1,i3
              projt(iloc,jloc,kloc,mloc) = h1*pt1(iloc,jloc,kloc,mloc)+
     .        h2*pt2(iloc,jloc,kloc,mloc) + h3*pt3(iloc,jloc,kloc,mloc)
            end do
          end do
        end do
      end do
 
c     Transform the [3;3;3;3] tensor into [6;6] matrix (Voigt notation)
c     projtm is the tensile projection matrix [6;6] of the effective stress
      do jloc=i1,i3
        do iloc=i1,i3
          projtm(iloc,jloc) = projt(iloc,iloc,jloc,jloc)
        end do
      end do
      do iloc=i1,i3
        projtm(iloc,i4) = 0.5d0*(projt(iloc,iloc,i1,i2) +
     .                    projt(iloc,iloc,i2,i1))
        projtm(iloc,i5) = 0.5d0*(projt(iloc,iloc,i1,i3) +
     .                    projt(iloc,iloc,i3,i1))
c       5th component is eps,xz
        projtm(iloc,i6) = 0.5d0*(projt(iloc,iloc,i2,i3) +
     .                    projt(iloc,iloc,i3,i2))
c       6th component is eps,yz
      end do
      do jloc=i1,i3
        projtm(i4,jloc) = 0.5d0*(projt(i1,i2,jloc,jloc) +
     .                    projt(i2,i1,jloc,jloc))
        projtm(i5,jloc) = 0.5d0*(projt(i1,i3,jloc,jloc) +
     .                    projt(i3,i1,jloc,jloc))
        projtm(i6,jloc) = 0.5d0*(projt(i2,i3,jloc,jloc) +
     .                    projt(i3,i2,jloc,jloc))
 
      end do
      projtm(i4,i4) = 0.5d0*(projt(i1,i2,i1,i2) +
     .                projt(i1,i2,i2,i1))
      projtm(i4,i5) = 0.5d0*(projt(i1,i2,i1,i3) +
     .                projt(i1,i2,i3,i1))
      projtm(i4,i6) = 0.5d0*(projt(i1,i2,i2,i3) +
     .                projt(i1,i2,i3,i2))
      projtm(i5,i4) = 0.5d0*(projt(i1,i3,i1,i2) +
     .                projt(i1,i3,i2,i1))
      projtm(i5,i5) = 0.5d0*(projt(i1,i3,i1,i3) +
     .                projt(i1,i3,i3,i1))
      projtm(i6,i5) = 0.5d0*(projt(i2,i3,i1,i3) +
     .                projt(i2,i3,i3,i1))
      projtm(i6,i4) = 0.5d0*(projt(i2,i3,i1,i2) +
     .                projt(i2,i3,i2,i1))
      projtm(i5,i6) = 0.5d0*(projt(i1,i3,i2,i3) +
     .                projt(i1,i3,i3,i2))
      projtm(i6,i6) = 0.5d0*(projt(i2,i3,i2,i3) +
     .                projt(i2,i3,i3,i2))
 
c     Identity matrix [6,6] U used to calculate P- and Q-
c       for P-: not really necessary as sig- = sig - sig+
c       for Q-: necessary for calculation of consistent algorit. tangent
c     Projection operators of the effective stresses P-
      do jloc=i1,i6
        do iloc=i1,i6
          projcm(iloc,jloc) = U(iloc,jloc) - projtm(iloc,jloc)
        end do
      end do
 
c     ===================================================================
c     !                                                                 !
c     !    POSITIVE AND NEGATIVE PARTS OF THE EFFECTIVE STRESS TENSOR   !
c     !    ==========================================================   !
c     ===================================================================
 
c     projtm has to be transposed
      call mulATB(projtm,esigma6v,esigmat,i6,i6,i6,i1)
c            [6;1] = [6;6] x [6x1]
      do iloc=i1,i6
        esigmac(iloc) = esigma6v(iloc) - esigmat(iloc)
      end do
 
c     Save the negative part of the effective stress tensor to be used in
c     the transient creep strain calculation
      do iloc=i1,i6
        parahot3(23+iloc) = esigmac(iloc)
      end do
 
c     ========================================
c     !                                      !
c     !           NOMINAL STRESSES           !
c     !           ================           !
c     ========================================
 
c     SIGMA = (1-dt) * SIG,eff,tension + (1-dc) * SIG,eff,comp
 
c     concrete completely damaged
      if ((lfulldamage).or.(dc.gt.0.98).or.(dt.gt.0.98)) then
        do iloc=1,6
          SIGMA(iloc) = 0.d0
        end do
        do iloc=i1,i6
          parahot3(idimpara3-i13+iloc) = sigma(iloc)
        end do
        do jloc=1,6
          do iloc=1,6
           H66t(iloc,jloc) = 0.d0
          end do
        end do
        return
      endif
 
      dt = parahot3(idimpara3-33)
      do iloc = i1,i6
        SIGMA(iloc) = (1.0d0-dt)*esigmat(iloc) +
     .                (1.0d0-dc)*esigmac(iloc)
      end do
 
      do iloc=i1,i6
        parahot3(idimpara3-i13+iloc) = sigma(iloc)
      end do
 
c     ***************************************************************
c     ========================================
c     !                                      !
c     !            TANGENT MODULUS           !
c     !            ===============           !
c     ========================================
 
c     Calculation of (nominal) consistent algorithmic tangent modulus Ct.
c     In DPRan3D we have calculated the effective elastoplastic tangent
c     modulus Dt. Now, we have to include the damage contribution.
 
c     APEX of Drucker-Prager or equivalent stress of Rankine = 0
c     First, verify that we are not at the apex of the surfaces
c     If we are, Dt = 0 and in this case we should have Ct = 0 too
      rloc = 0.d0
      do jloc=1,6
        do iloc=1,6
          rloc = rloc+ABS(H66t(iloc,jloc))
        enddo
      enddo
      if (rloc.lt.1.e+4) then
        do jloc=1,6
          do iloc=1,6
           H66t(iloc,jloc) = 0.d0
          end do
        end do
        return
      endif
c     End of test APEX
 
      do iloc=i1,6
        rt1(i1,iloc) = ( (at/2.)*EXP(-at*rkappat)+(6.*at/2.)*
     .                EXP(-6.*at*rkappat) ) * rt1(i1,iloc)
        rc1(i1,iloc) = ( parahot3(14)*EXP(-parahot3(14)*rkappac) )
     .                                       * rc1(i1,iloc)
      end do
 
c     Calculation of W4
c     Calculation of p12, p23, p13
c     p12(i3,i3) = (1/2) (v1 X v2 + v2 X v1) (tensorial product)
      call mulABT(v1,v2,p121,i3,i1,i3,i1)
c          [3;3]=[3;1]x[1;3]
      call mulABT(v2,v1,p122,i3,i1,i3,i1)
      call mulABT(v2,v3,p231,i3,i1,i3,i1)
      call mulABT(v3,v2,p232,i3,i1,i3,i1)
      call mulABT(v1,v3,p131,i3,i1,i3,i1)
      call mulABT(v3,v1,p132,i3,i1,i3,i1)
      do jloc=i1,i3
        do iloc=i1,i3
          p12(iloc,jloc) = (1.0d0/2.0d0) * (p121(iloc,jloc) +
     .                                      p122(iloc,jloc))
          p23(iloc,jloc) = (1.0d0/2.0d0) * (p231(iloc,jloc) +
     .                                      p232(iloc,jloc))
          p13(iloc,jloc) = (1.0d0/2.0d0) * (p131(iloc,jloc) +
     .                                      p132(iloc,jloc))
        end do
      end do
 
c     Calculation of PPijkl = pij X pij
c     [3;3;3;3]=[3;3] X [3;3]
      do mloc=i1,i3
        do kloc=i1,i3
          do jloc=i1,i3
            do iloc=i1,i3
              pp12(iloc,jloc,kloc,mloc)=p12(iloc,jloc)*p12(kloc,mloc)
              pp23(iloc,jloc,kloc,mloc)=p23(iloc,jloc)*p23(kloc,mloc)
              pp13(iloc,jloc,kloc,mloc)=p13(iloc,jloc)*p13(kloc,mloc)
            end do
          end do
        end do
      end do
 
c     Calculation of Q+           ! [3;3;3;3]
      do mloc=i1,i3
        do kloc=i1,i3
          do jloc=i1,i3
            do iloc=i1,i3
              Qtens(iloc,jloc,kloc,mloc) = r0
            end do
          end do
        end do
      end do
      if ((S1-S2).ne.r0) then
        rloc = (h1*S1) - (h2*S2)
        do mloc=i1,i3
          do kloc=i1,i3
            do jloc=i1,i3
              do iloc=i1,i3
                Qtens(iloc,jloc,kloc,mloc)=Qtens(iloc,jloc,kloc,mloc)
     .                  + 2.0d0 * (rloc/(S1-S2))
     .                          *  pp12(iloc,jloc,kloc,mloc)
              end do
            end do
          end do
        end do
      else
        continue
      endif
      if ((S2-S3).ne.r0) then
        rloc = (h2*S2) - (h3*S3)
        do mloc=i1,i3
          do kloc=i1,i3
            do jloc=i1,i3
              do iloc=i1,i3
                Qtens(iloc,jloc,kloc,mloc)=Qtens(iloc,jloc,kloc,mloc)
     .                  + 2.0d0 * (rloc/(S2-S3))
     .                          * pp23(iloc,jloc,kloc,mloc)
              end do
            end do
          end do
        end do
      else
        continue
      endif
      if ((S1-S3).ne.r0) then
        rloc = (h1*S1) - (h3*S3)
        do mloc=i1,i3
          do kloc=i1,i3
            do jloc=i1,i3
              do iloc=i1,i3
                Qtens(iloc,jloc,kloc,mloc)=Qtens(iloc,jloc,kloc,mloc)
     .                  + 2.0d0 * (rloc/(S1-S3))
     .                          * pp13(iloc,jloc,kloc,mloc)
              end do
            end do
          end do
        end do
      else
        continue
      endif
      do mloc=i1,i3
        do kloc=i1,i3
          do jloc=i1,i3
            do iloc=i1,i3
              Qtens(iloc,jloc,kloc,mloc) = projt(iloc,jloc,kloc,mloc)
     .                                  + Qtens(iloc,jloc,kloc,mloc)
c              [3;3;3;3]
            end do
          end do
        end do
      end do
 
c     Calculation of the matrix [6;6] Qtens
c     Transform the [3;3;3;3] Qtens tensor into [6;6] matrix (Voigt notation)
      do jloc=1,3
        do iloc=1,3
          Qtensm(iloc,jloc) = Qtens(iloc,iloc,jloc,jloc)
        end do
      end do
      do iloc=1,3
        Qtensm(iloc,4) = 0.5d0*(Qtens(iloc,iloc,1,2)
     .                        + Qtens(iloc,iloc,2,1))
        Qtensm(iloc,5) = 0.5d0*(Qtens(iloc,iloc,1,3)
     .                        + Qtens(iloc,iloc,3,1))
        Qtensm(iloc,6) = 0.5d0*(Qtens(iloc,iloc,2,3)
     .                        + Qtens(iloc,iloc,3,2))
      end do
      do jloc=1,3
        Qtensm(4,jloc) = 0.5d0*(Qtens(1,2,jloc,jloc)
     .                        + Qtens(2,1,jloc,jloc))
        Qtensm(5,jloc) = 0.5d0*(Qtens(1,3,jloc,jloc)
     .                        + Qtens(3,1,jloc,jloc))
        Qtensm(6,jloc) = 0.5d0*(Qtens(2,3,jloc,jloc)
     .                        + Qtens(3,2,jloc,jloc))
      end do
      Qtensm(4,4) = 0.5d0*(Qtens(1,2,1,2) + Qtens(1,2,2,1))
      Qtensm(4,5) = 0.5d0*(Qtens(1,2,1,3) + Qtens(1,2,3,1))
      Qtensm(4,6) = 0.5d0*(Qtens(1,2,2,3) + Qtens(1,2,3,2))
      Qtensm(5,4) = 0.5d0*(Qtens(1,3,1,2) + Qtens(1,3,2,1))
      Qtensm(5,5) = 0.5d0*(Qtens(1,3,1,3) + Qtens(1,3,3,1))
      Qtensm(5,6) = 0.5d0*(Qtens(1,3,2,3) + Qtens(1,3,3,2))
      Qtensm(6,4) = 0.5d0*(Qtens(2,3,1,2) + Qtens(2,3,2,1))
      Qtensm(6,5) = 0.5d0*(Qtens(2,3,1,3) + Qtens(2,3,3,1))
      Qtensm(6,6) = 0.5d0*(Qtens(2,3,2,3) + Qtens(2,3,3,2))
 
      do jloc=1,6
        do iloc=1,6
          Qcompm(iloc,jloc) = U(iloc,jloc) - Qtensm(iloc,jloc)
        enddo
      enddo
      do iloc=1,3
        Qcompm(iloc,4) = 0.5d0*U(iloc,4) - Qtensm(iloc,4)
        Qcompm(iloc,5) = 0.5d0*U(iloc,5) - Qtensm(iloc,5)
        Qcompm(iloc,6) = 0.5d0*U(iloc,6) - Qtensm(iloc,6)
      enddo
      do jloc=1,3
        Qcompm(4,jloc) = 0.5d0*U(4,jloc) - Qtensm(4,jloc)
        Qcompm(5,jloc) = 0.5d0*U(5,jloc) - Qtensm(5,jloc)
        Qcompm(6,jloc) = 0.5d0*U(6,jloc) - Qtensm(6,jloc)
      enddo
      do jloc=4,6
        do iloc=4,6
          Qcompm(iloc,jloc) = 0.5d0*U(iloc,jloc) - Qtensm(iloc,jloc)
        enddo
      enddo
 
      do jloc=i1,6
        do iloc=i1,6
          w4m(iloc,jloc) = dt*Qtensm(iloc,jloc) + dc*Qcompm(iloc,jloc)
        end do
      end do
 
c     Calculation of (nominal) consistent algorithmic tangent modulus Ct.
c          [6;6]
c     Ct = (U - w4m) * Dt - sig,eff,+ * rt1 - sig,eff,- * rc1
c     In this subroutine Dt is called H66t (comes from DPRan3d)
c     and at the end, Ct has to be stored into H66t
 
c     Calculation of the product (U - w4m) * Dt
      do jloc=i1,i6
        do iloc=i1,i6
          w4m(iloc,jloc) = U(iloc,jloc) - w4m(iloc,jloc)
        end do
      end do
      call mulAB(w4m,H66t,H66temp,i6,i6,i6,i6)
c         [6;6]
      do j=1,6
        do i=1,6
          H66t(i,j) = H66temp(i,j)
        enddo
      enddo
 
c     esigmat and esigmac have been calculated above
      call mulAB(esigmat,rt1,sigrt,i6,i1,i1,i6)
c           [6;6] = [6;1] * [1;6]
      call mulAB(esigmac,rc1,sigrc,i6,i1,i1,i6)
 
      do jloc=i1,i6
        do iloc=i1,i6
          H66t(iloc,jloc) = H66t(iloc,jloc)
     .                       - sigrt(iloc,jloc)
     .                       - sigrc(iloc,jloc)
        end do
      end do
 
c    END OF TANGENT MODULUS
 
      RETURN
      END
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

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