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