Numérotation des lignes :

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 variablesc     idimpara3 is the number of parameters for parahot3c     H66t is the tangent matrixc     rt1 and rc1 are the vectors built in DPRAN3D and used in DAMTG3D for calculation of tangent matrixc     esigma6v is the effective stress vectorc     lg and lerror are logicalc     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 stepc     parahot3(idimpara3-35): equivalent plastic strain in tension     at the end of the current stepc     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 Uc      -----------------------       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 id4C     do iloc=1,3C       do jloc=1,3C         do kloc=1,3C           do mloc=1,3C             if ((iloc.eq.jloc).and.(kloc.eq.mloc)) thenC               id4(iloc,jloc,kloc,mloc) = 1.0d0C             elseC               id4(iloc,jloc,kloc,mloc) = 0.0d0C             endifC           end doC         end doC       end doC     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,ic     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. tangentc     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 inc     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 tangentc     modulus Dt. Now, we have to include the damage contribution. c     APEX of Drucker-Prager or equivalent stress of Rankine = 0c     First, verify that we are not at the apex of the surfacesc     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      endifc     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 W4c     Calculation of p12, p23, p13c     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 pijc     [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] Qtensc     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,- * rc1c     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