C DPRAN3D SOURCE CB215821 25/04/08 21:15:14 12227 SUBROUTINE dpran3d(parahot3,idimpara3,H66t,rt1,rc1,esigma, . deps,lg,ntot,lerror,i1,i2,i3,i4,i5,i6,n,lnoconv, . istep,lcp,lbroyden,H66,lfulldamage) 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 used in DAMTG3D for the calculation of the tangent matrix c esigma is the effective stress vector c deps is the incremental strains vector c lg and lerror are logical c ntot is the integration point (global number for the whole structure): counter c n, istep, lnoconv, lbroyden and lfulldamage are used for the algorithmic strategy (step subdivider) c H66 is the elastic matrix c =========================================================================================== c ! Plastic Box of the 3D plastic-damage concrete model by T. Gernay ! c ! Yield functions of Drucker Prager in compression and Rankine in tension ! c ! Calculate the effective stresses and the effective elastoplastic tangent ! c ! modulus Dt (implicit process) ! c ! ! c ! input ! c ! DEPS : Incremental strain vector (6 components) ! c ! ! c ! output ! c ! H66t : tangent stiffness matrix (6x6) ! c ! ! c ! input and output ! c ! PARAHOT3 : the first columns contain the material properties at elevated temperature ! c ! the last columns contain strains, stresses, .... ! c ! ! c ! ESIGMA : Effective stress vector (6 components) ! 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 deps(i6) c Incremental strain since the end of the previous time step (input) c i.e. since the last converged point dimension H66t(i6,i6) c Tangent matrix (output) dimension parahot3(idimpara3) c parahot3: Material prop. and various info. at elevated temp. (input and output) c parahot3(idimpara3-37,ntot): equivalent plastic strain in tension at the last converged step c parahot3(idimpara3-35,ntot): equivalent plastic strain in tension at the end of the current step c parahot3(idimpara3-36,ntot): equivalent plastic strain in compression at the last converged step c parahot3(idimpara3-34,ntot): equivalent plastic strain in compression at the end of the current step c parahot3(idimpara3-30:idimpara3-25,ntot) = strain at the end of the previous time step c parahot3(idimpara3-24:idimpara3-19,ntot) = plastic strains at the beginning of dpran3D c parahot3(idimpara3-18:idimpara3-13,ntot) = strain at the end of dpran3D c parahot3(idimpara3- 6:idimpara3- 1,ntot) = effective stress at the end of dpran3D c parahot3(idimpara3,ntot) = thermal strain dimension esigma(i6) c Effective stress (output) **** dimension EPSMEC(i6),EPSSIG(i6),esigmae(i6),esigmaf(i6) dimension EPSMEC(6),EPSSIG(6),esigmae(6),esigmaf(6) **** dimension H66(i6,i6),x(i2),d2gt(i6,i6),trav66b(i6,i6),trav6(i6) dimension H66(6,6),x(2),d2gt(6,6),trav66b(6,6),trav6(6) **** dimension trav66(i6,i6),Da(i6,i6),dft(i6) dimension trav66(6,6),Da(6,6),dft(6) **** dimension dfc(i6),dgc(i6),dgt(i6) dimension dfc(6),dgc(6),dgt(6) **** dimension trav1(1),DftD(i1,i6),DfcD(i1,i6),Ddgt(i6),Ddgc(i6) dimension trav1(1),DftD(1,6),DfcD(1,6),Ddgt(6),Ddgc(6) **** dimension trav16(i1,i6),rn1(i1,i6),rn2(i1,i6),rn3(i1,i6) dimension trav16(1,6),rn1(1,6),rn2(1,6),rn3(1,6) **** dimension rt1(i1,i6),rc1(i1,i6),vloc1(i1),vloc(i6) dimension rt1(1,6),rc1(1,6),vloc1(1),vloc(6) **** dimension Ddgtn1(i6,i6),Ddgtn2(i6,i6),Ddgcn3(i6,i6) dimension Ddgtn1(6,6),Ddgtn2(6,6),Ddgcn3(6,6) **** dimension Ddgcn4(i6,i6),rn4(i1,i6),dplas(6),d2gc(6,6),vloc2(6) dimension Ddgcn4(6,6),rn4(1,6),dplas(6),d2gc(6,6),vloc2(6) c Deviatoric projection operator P2 (stressdeviator = P2 * stress) c -------------------------------- dimension P2(6,6) data P2 / 2.0d0, -1.0d0, -1.0d0, 0.0d0, 0.0d0, 0.0d0, . -1.0d0, 2.0d0, -1.0d0, 0.0d0, 0.0d0, 0.0d0, . -1.0d0, -1.0d0, 2.0d0, 0.0d0, 0.0d0, 0.0d0, . 0.0d0, 0.0d0, 0.0d0, 6.0d0, 0.0d0, 0.0d0, . 0.0d0, 0.0d0, 0.0d0, 0.0d0, 6.0d0, 0.0d0, . 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 6.0d0 / c Matrix of unity I, or U c ----------------------- dimension U(6,6) data U / 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, . 0.0d0, 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, . 0.0d0, 0.0d0, 1.0d0, 0.0d0, 0.0d0, 0.0d0, . 0.0d0, 0.0d0, 0.0d0, 1.0d0, 0.0d0, 0.0d0, . 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0, 0.0d0, . 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0 / c Vector pi2 c ----------------------- dimension pi2(6) data pi2 / 1.0d0, . 1.0d0, . 1.0d0, . 0.0d0, . 0.0d0, . 0.0d0 / ***** * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS) IF(I1.GT.1) PRINT *, ' DPRAND3D - ERREUR I1 = ', I1, ' > 1 ' IF(I6.GT.6) PRINT *, ' DPRAND3D - ERREUR I6 = ', I6, ' > 6 ' ***** lerror = .false. lnoconv = .false. P2(1,1) = 2.0d0 P2(2,2) = 2.0d0 P2(3,3) = 2.0d0 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. r2 = 2. r10 = 10. Eo = parahot3(i1) poison = parahot3(i2) fc = parahot3(i3) ft = parahot3(i4) c Strain at peak stress epsc1 = parahot3(i5) c Compressive limit of elasticity fc0 = parahot3(i6) c Biaxial compressive strength fb = parahot3(i7) c Dilatancy parameter ag = parahot3(i8) c Damage to peak stress dc = parahot3(i9) c Compressive dissipated energy parameter rxc = parahot3(i10) c Crack energy in tension gt = parahot3(i11) c Model parameters rkappa1 = parahot3(i13) ac = parahot3(i14) bc = parahot3(i15) c small strain for starting the return mapping algorithm tol = 0.0001*epsc1 c Hardening parameter for tension at the beginning of the step rkappait= parahot3(idimpara3-37) c Hardening parameter for compression at the beginning of the step rkappaic= parahot3(idimpara3-36) c =================================================================== c ! ! c ! CALCULATION OF THE MATRIX OF ELASTICITY ! c ! ======================================= ! c =================================================================== c Matrix of elasticity, (6x6, full 3D) c -------------------- c [Eo*(1-nu)/((1+nu)*(1-2*nu)) Eo*nu/((1+nu)*(1-2*nu)) Eo*nu/((1+nu)*(1-2*nu)) 0 0 0 ] c [Eo*nu/((1+nu)*(1-2*nu)) Eo*(1-nu)/((1+nu)*(1-2*nu)) Eo*nu/((1+nu)*(1-2*nu)) 0 0 0 ] c H66 = [Eo*nu/((1+nu)*(1-2*nu)) Eo*nu/((1+nu)*(1-2*nu)) Eo*(1-nu)/((1+nu)*(1-2*nu)) 0 0 0 ] c [0 0 0 Eo/(1+nu) 0 0 ] c [0 0 0 0 Eo/(1+nu) 0 ] c [0 0 0 0 0 Eo/(1+nu) ] c Matrix of elasticity, (3x3, plane stress) c -------------------- c [ 1 nu 0 0 0 0] c H66 = Eo/(1-nu*nu) [nu 1 0 0 0 0] c [ 0 0 0 0 0 0] c [ 0 0 0 (1-nu)/2 0 0] c [ 0 0 0 0 0 0] c [ 0 0 0 0 0 0] c++ do jloc=i1,i6 do iloc=i1,i6 H66(iloc,jloc) = r0 end do end do if (.not.lcp) then c 3D model rloc = Eo/((r1+poison)*(r1-2.0d0*poison)) H66(i1,i1) = rloc*(r1-poison) H66(i2,i2) = rloc*(r1-poison) H66(i3,i3) = rloc*(r1-poison) H66(i1,i2) = rloc*poison H66(i1,i3) = rloc*poison H66(i2,i3) = rloc*poison H66(i2,i1) = rloc*poison H66(i3,i1) = rloc*poison H66(i3,i2) = rloc*poison rloc2 = Eo/(2.0d0*(r1+poison)) H66(i4,i4) =rloc2 H66(i5,i5) =rloc2 H66(i6,i6) =rloc2 else c PLANE STRESS rloc = Eo/(r1-poison*poison) H66(i1,i1) = rloc H66(i2,i1) = poison * (rloc) H66(i1,i2) = poison * (rloc) H66(i2,i2) = rloc H66(4,4) = rloc * ((1.0d0-poison)/(2.0d0)) deps(3) = 0.d0 deps(5) = 0.d0 deps(6) = 0.d0 endif c ============================================================= c ! ! c ! ELASTIC STRESS PREDICTOR ! c ! ======================== ! c ============================================================= c Total mechanical strain: EPSMEC[i,j] = epsmec[i-1,conv] + deps[i,j] - epspl[i-1,j] c epsplastic has to be subtracted because Sigma initial = 0. do iloc=i1,i6 EPSMEC(iloc) = parahot3(idimpara3-31+iloc) . + ( (DEPS(iloc) * istep)/n ) . - parahot3(idimpara3-25+iloc) c Calculation of the INSTANTANEOUS STRESS-DEPENDENT STRAIN EPSSIG(iloc) = EPSMEC(iloc) - parahot3(17+iloc) end do c The effective stresses are calculated as functions of the instantaneous c stress-dependent strain eps,sig because the transient creep strain is explicitly c calculated apart from these constitutive relationships. Still in parahot, the mechanical c strains are saved, not the instantaneous stress-dependent strains. if ( (.not.lnoconv).and.(n.gt.1) ) then c We are in a subdivision of the step (for algorithmic reasons) c The plastic strain considered in epsmec has to be adjusted rkappait = parahot3(idimpara3-35) rkappaic = parahot3(idimpara3-34) x1 = rkappait - parahot3(idimpara3-37) x2 = rkappaic - parahot3(idimpara3-36) call dplas3D(x1,x2,P2,pi2,ag,fc,fb,parahot3,idimpara3, . dplas,lcp) do iloc=i1,i6 EPSSIG(iloc) = EPSMEC(iloc)-dplas(iloc)-parahot3(17+iloc) end do endif c Elastic stress predictor: eff sigma[i,j] = Hel x EPSSIG[i,j] call mulab(H66,EPSSIG,ESIGMAE,i6,i6,i6,i1) * * print *,' prediction contraintes = ',esigmae * c ============================================================== c ! ! c ! CONCRETE COMPLETELY DAMAGED ! c ! =========================== ! c ============================================================== dc = parahot3(idimpara3-40) dt = parahot3(idimpara3-41) if ((dc.gt.0.98).or.(lfulldamage).or.(dt.gt.0.98)) then c Concrete completely damaged do iloc=1,6 ESIGMA(iloc) = 0.d0 end do do iloc=1,6 parahot3(idimpara3-13+iloc) = esigma(iloc) end do do jloc=1,6 do iloc=1,6 H66t(iloc,jloc) = 0.d0 end do end do c store the mechanical strain at the end of this iteration c mechanical strain = epsmec + epspl do iloc=i1,i6 parahot3(idimpara3-19+iloc) = EPSMEC(iloc) + . parahot3(idimpara3-25+iloc) end do return endif c ========================================================================== c ! ! c ! CALCULATE THE PLASTICITY FUNCTIONS WITH THE ELASTIC STRESS PREDICTOR ! c ! ==================================================================== ! c ========================================================================== c Calculation of the hardening functions tau_compr = fhardcc(rkappaic,fc,fc0,rkappa1, . ac,bc,lerror) tau_tension = fhardct(rkappait,ft) if (lerror) then write(2,*) '>Error returned from compressive hardening . to DPRAN3d' return endif c Calculation of the yield functions for the elastic stress predictor alpha = (fb-fc)/(2.*fb-fc) f_compr = DPfunc(esigmae,alpha) f_tension = Ranfunc(esigmae,i3,i6,parahot3, . idimpara3,lcp) c++ f_compr = f_compr - ( (1.0d0-alpha) * tau_compr ) f_tension = f_tension - tau_tension c ======================================================= c ! ! c ! TEST ELASTIC OR PLASTIC STEP ! c ! ============================ ! c ======================================================= c The plasticity functions must be satisfied with a precision of 10 Pa, i.e. 0.000010 N/mm² if ((f_compr.lt.-r10).and.(f_tension.lt.-r10)) then C *************** c ! ELASTICITY ! C *************** do iloc=i1,i6 ESIGMA(iloc) = ESIGMAE(iloc) end do do jloc=i1,i6 do iloc=i1,i6 H66T(iloc,jloc) = H66(iloc,jloc) end do end do c store the stress at the end of this iteration do iloc=i1,i6 parahot3(idimpara3-7+iloc) = ESIGMAE(iloc) end do c The equivalent plastic strain in tension has not changed parahot3(idimpara3-35) = parahot3(idimpara3-37) c The equivalent plastic strain in compression has not changed parahot3(idimpara3-34) = parahot3(idimpara3-36) c vectors for the calculation of H66t in subroutine damtg3d do iloc=i1,i6 rt1(i1,iloc) = r0 rc1(i1,iloc) = r0 end do else C *************** c ! PLASTICITY ! C *************** * PRINT *, ' >>>>>>>> ON ENTRE EN PLASTICITE' c RETURN MAPPING c -------------- c Calculation of Dlambdac and Dlambdat, the increments of the eq. plastic strain which brings c back on the plasticity surfaces c We have to solve the following non linear system of equations c rkt * function_t(Dlambdac,Dlambdat) + (1-rkt) * Dlambdat = 0 c rkc * function_c(Dlambdac,Dlambdat) + (1-rkc) * Dlambdac = 0 c with rkc and rkt = 1 or 0, depending whether the related surface of plasticity is violated if (f_compr.lt.-r10) then c only the tension surface is violated by the trial stress Dlambdac = r0 rkc = r0 rkt = r1 else if (f_tension.lt.-r10) then c only the compression surface is violated by the trial stress Dlambdat = r0 rkc = r1 rkt = r0 else c Both surfaces are violated by the trial stress rkc = r1 rkt = r1 endif c X=0 at the beginning of every SAFIR iteration, X being the vector (Dlambdat,Dlambdac) X(i1) = r0 X(i2) = r0 c =================================================================================== c ! ! c ! CALL TO SUBROUTINE ALGOPLAS TO SOLVE THE NONLINEAR SYSTEM IN DLAMBDAT, DLAMBDAC ! c ! =============================================================================== ! c =================================================================================== if (.not.lcp) then c 3D MODEL if (.not.lbroyden) then call Algoplas1(rkt,rkc,esigmae,H66,P2,pi2,alpha, . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot, . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6, . parahot3,idimpara3,deps,lnoconv,lcp,U) else call Algoplas2(rkt,rkc,esigmae,H66,P2,pi2,alpha, . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot, . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6, . parahot3,idimpara3,deps,lnoconv,lcp,U) endif if (lnoconv) then return endif if (lerror) then write(2,*) 'lerror is returned from ALGOPLAS to . DPRAN3D. Ntot =',ntot return endif c The increment of plastic multipliers Dlambdat, Dlambdac have been found c ========================================================== c ! ! c ! UPDATE OF THE HARDENING PARAMETERS AND THE STRESSES ! c ! =================================================== ! c ========================================================== Dlambdat = x(i1) Dlambdac = x(i2) if (n.gt.1) then c we are currently in a step subdivision (for algorithmic reason) c we take the value saved at the end of the previous step-subdivision rkappait = parahot3(idimpara3-35) rkappaic = parahot3(idimpara3-34) else c we take the value at the previous (converged) time step rkappait= parahot3(idimpara3-37) rkappaic= parahot3(idimpara3-36) endif c update of the hardening parameters rkappat = rkappait+Dlambdat rkappac = rkappaic+Dlambdac c Update of the stresses: esigmaf = F(Dlambdat,Dlambdac (and esigmae)) ag = parahot3(8) call sigf3D_implicit(esigmae,Dlambdat,Dlambdac,esigmaf,H66, c *************** . P2,pi2,ag,i6,fc,fb,parahot3,idimpara3,lerror,lcp, . esigmae,U) if (lerror) then return endif else c PLANE STRESS MODEL if (.not.lbroyden) then call Algoplas1cp(rkt,rkc,esigmae,H66,P2,pi2,alpha, . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot, . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6, . parahot3,idimpara3,deps,esigmaf,rkappat,rkappac, . lnoconv,lcp,U) else call Algoplas2cp(rkt,rkc,esigmae,H66,P2,pi2,alpha, . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot, . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6, . parahot3,idimpara3,deps,esigmaf,rkappat,rkappac, . lnoconv,lcp,U) endif if (lnoconv) then return endif if (lerror) then write(2,*) 'lerror is returned from ALGOPLAS to . DPRAN3D. Ntot =',ntot return endif Dlambdat = x(i1) Dlambdac = x(i2) if (n.gt.1) then c we are currently in a step subdivision (for algorithmic reason) c we take the value saved at the end of the previous step-subdivision rkappait = parahot3(idimpara3-35) rkappaic = parahot3(idimpara3-34) else c we take the value at the previous (converged) time step rkappait= parahot3(idimpara3-37) rkappaic= parahot3(idimpara3-36) endif c update of the hardening parameters rkappat = rkappait+Dlambdat rkappac = rkappaic+Dlambdac endif c store the stress at the end of this iteration do iloc=i1,i6 ESIGMA(iloc) = ESIGMAF(iloc) end do do iloc=i1,i6 parahot3(idimpara3-i7+iloc) = ESIGMAF(iloc) end do c store the accumulated plastic strains at the end of this iteration parahot3(idimpara3-35) = rkappat parahot3(idimpara3-34) = rkappac c THE EFFECTIVE STRESSES HAVE BEEN FOUND c ************************************** c ================================================================ c ! TANGENT MATRIX. (effective!) ! c ! ---------------------------- ! c ! Here, we calculate the effective elastoplastic tangent ! c ! modulus Dt. As this is a damage-plastic model, we have ! c ! then to correct this modulus Dt with the damage terms to ! c ! obtain the (nominal) consistent algorithmic tangent modulus ! c ! Ct. This will be done in DAMTG3D. Here, we calculate Dt. ! c ================================================================ c APEX of Drucker-Prager surface alpha = (fb-fc)/(2.*fb-fc) f_compr = DPfunc(esigmae,alpha) if ((ABS(f_compr).le.100.).and.(rkc.eq.1)) then do jloc=i1,i6 do iloc=i1,i6 H66t(iloc,jloc) = 0.d0 end do enddo do iloc=i1,i6 parahot3(idimpara3-19+iloc) = EPSMEC(iloc) + . parahot3(idimpara3-25+iloc) end do RETURN endif if (lcp) then c recalculate the full 3D matrix do jloc=i1,i6 do iloc=i1,i6 H66(iloc,jloc) = r0 end do end do rloc = Eo/((r1+poison)*(r1-2.0d0*poison)) H66(i1,i1) = rloc*(r1-poison) H66(i2,i2) = rloc*(r1-poison) H66(i3,i3) = rloc*(r1-poison) H66(i1,i2) = rloc*poison H66(i1,i3) = rloc*poison H66(i2,i3) = rloc*poison H66(i2,i1) = rloc*poison H66(i3,i1) = rloc*poison H66(i3,i2) = rloc*poison rloc2 = Eo/(2.0d0*(r1+poison)) H66(i4,i4) =rloc2 H66(i5,i5) =rloc2 H66(i6,i6) =rloc2 endif if (lcp) then tau_tension = fhardct(rkappat,ft) psi1 = tau_tension - (esigmaf(1)+esigmaf(2))/2 if (abs(psi1).le.fc*1.D-5) then do jloc=1,6 do iloc=1,6 H66t(iloc,jloc) = 0 enddo enddo do iloc=i1,i6 parahot3(idimpara3-19+iloc) = EPSMEC(iloc) + . parahot3(idimpara3-25+iloc) end do RETURN endif endif c First, calculation of Da, the effective algorithmic modulus c ************************************************** call ddsigrank(esigmaf,d2gt,i1,i2,i3,i4,i5,i6,lcp, . parahot3,idimpara3) do jloc=i1,i6 do iloc=i1,i6 TRAV66b(iloc,jloc) = Dlambdat*d2gt(iloc,jloc) c [6;6] end do end do call mulAB(H66,TRAV66b,TRAV66b,i6,i6,i6,i6) c [6;6] = [6;6] x [6;6] tau_c = fhardcc(rkappac,fc,fc0,rkappa1, . ac,bc,lerror) if (lerror) then write(2,*) '>Error returned from compressive hardening . to DPRAN3d' return endif call mulATB(pi2,esigmaf,vloc1,i6,i1,i6,i1) c [1] = [1;6] x [6;1] psi2 = (1d0-alpha)*tau_c - alpha*vloc1(i1) c First come the second order terms. call mulab(P2,ESIGMAF,TRAV6,i6,i6,i6,i1) c [6;1] = [6;6] x [6;1] c t t t c T66 = T6 T6 = P2 s s P2 call mulABT(trav6,trav6,trav66,i6,i1,i6,i1) c [6;6] = [6;1] x [1;6] test = 1.D+6 psi2lim=max(1.D-7*test,1.d-3) if (abs(psi2).le.psi2lim) then do jloc=i1,i6 do iloc=i1,i6 TRAV66(iloc,jloc) = r0 end do end do else do jloc=i1,i6 do iloc=i1,i6 TRAV66(iloc,jloc) = (P2(iloc,jloc)-TRAV66(iloc,jloc) . /(r2*psi2*psi2))/(r2*psi2) c [6;6] end do end do endif do jloc=i1,i6 do iloc=i1,i6 TRAV66(iloc,jloc) = Dlambdac*TRAV66(iloc,jloc) end do end do call mulAB(H66,TRAV66,TRAV66,i6,i6,i6,i6) c **** ***** [6;6]=[6;6]x[6;6] do jloc=i1,i6 do iloc=i1,i6 TRAV66(iloc,jloc) = U(iloc,jloc)+TRAV66(iloc,jloc) c TRAV66(iloc,jloc)=U(iloc,jloc)+TRAV66(iloc,jloc) . +TRAV66b(iloc,jloc) c TRAV66b requires the calculation of d2gt - we prefer not to take it end do end do c t t -1 c Da = ( I + Dlambdat H66 d2gt + Dlambdac H66 (P2/2psi2)-(P2 s s P2)/(4psi2³) ) H66 call invert9(TRAV66,i6,i6) c **** ******* [6;6] call mulab(TRAV66,H66,Da,i6,i6,i6,i6) c [6;6]=[6;6]x[6;6] c The effective algorithmic modulus Da has been obtained c Now, calculation of the effective elastoplastic tangent modulus Dt c ****************************************************************** c Now come the first order terms dtau1 = dfhardct(rkappat) dtau2 = dfhardcc(rkappac,fc,fc0,rkappa1,ac,bc) c Calculation of n1, n2, n3, n4 call dsigRank(esigmaf,dgt,i3,i6,parahot3,idimpara3, c **** ******** [6;1] . rkappait,lcp) rloc = 0.d0 do iloc=1,6 rloc = rloc + (dgt(iloc))**2 enddo if (rloc.eq.0.d0) rkt = 0.d0 do iloc=i1,i6 dft(iloc) = rkt*dgt(iloc) c [6;1] = [1] x [6;1] end do e1 = r1-rkt*(r1+dtau1) c [1] e2 = r1-rkc*(r1+dtau2) c [1] fc = parahot3(i3) c Compressive strength fb = parahot3(i7) c Biaxial compressive strength call dsigDrucker(fc,fb,esigmaf,dfc,alpha,P2, c **** **************** [6;1] . pi2,i1,i6,lcp) do iloc=i1,i6 dfc(iloc) = rkc*dfc(iloc) c [6;1] = [1]x[6;1] end do ag = parahot3(i8) call dsiggc(esigmaf,dgc,ag,P2,pi2,i1,i6,fc,fb,lcp) c **** **************** [6;1] call mulATB(dft,Da,trav16,i6,i1,i6,i6) c [1;6] = [1;6]x[6;6] do iloc=i1,i6 dftD(i1,iloc) = trav16(i1,iloc) end do call mulab(Da,dgt,Ddgt,i6,i6,i6,i1) c [6;1] = [6;6] x [6;1] call mulATB(dfc,Da,trav16,i6,i1,i6,i6) c [1;6] = [1;6]x[6;6] do iloc=i1,i6 dfcD(i1,iloc) = trav16(i1,iloc) end do call mulab(Da,dgc,Ddgc,i6,i6,i6,i1) c [6;1] = [6;6]x[6;1] c n1 call mulATB(dft,Ddgc,trav1,i6,i1,i6,i1) c **** ****** [1]=[1;6]x[6;1] dftDdgc = trav1(i1) c [1] do iloc=i1,i6 rn1(i1,iloc) = -dftDdgc*dfcD(i1,iloc) end do c [1;6]=[1]x[1;6] c n2 [1]=[1;6]x[6;1] call mulATB(dfc,Ddgc,trav1,i6,i1,i6,i1) c **** ****** dfcDdgc = trav1(i1) do iloc=i1,i6 rn2(i1,iloc) = (dfcDdgc-e2)*dftD(i1,iloc) end do c n3 [1]=[1;6]x[6;1] call mulATB(dft,Ddgt,trav1,i6,i1,i6,i1) c **** ****** dftDdgt = trav1(i1) do iloc=i1,i6 rn3(i1,iloc) = (dftDdgt-e1)*dfcD(i1,iloc) end do c n4 [1]=[1;6]x[6;1] call mulATB(dfc,Ddgt,trav1,i6,i1,i6,i1) c **** ****** dfcDdgt = trav1(i1) do iloc=i1,i6 rn4(i1,iloc) = -(dfcDdgt)*dftD(i1,iloc) end do c denominator [1] Denom = (dftDdgt-e1)*(dfcDdgc-e2) - (dfcDdgt*dftDdgc) c Calculation of rt1 and rc1 that are used in Mater 2 for calculation of Ct do iloc=i1,i6 rt1(i1,iloc) = (rn1(i1,iloc)+rn2(i1,iloc))/denom rc1(i1,iloc) = (rn3(i1,iloc)+rn4(i1,iloc))/denom end do do iloc=i1,i6 vloc(iloc) = rn1(i1,iloc) end do call mulABT(Ddgt,vloc,Ddgtn1,i6,i1,i6,i1) c **** ****** [6;6]=[6;1]x[1;6] do iloc=i1,i6 vloc(iloc) = rn2(i1,iloc) end do call mulABT(Ddgt,vloc,Ddgtn2,i6,i1,i6,i1) c **** ****** [6;6]=[6;1]x[1;6] do iloc=i1,i6 vloc(iloc) = rn3(i1,iloc) end do call mulABT(Ddgc,vloc,Ddgcn3,i6,i1,i6,i1) c **** ****** [6;6]=[6;1]x[1;6] do iloc=i1,i6 vloc(iloc) = rn4(i1,iloc) end do call mulABT(Ddgc,vloc,Ddgcn4,i6,i1,i6,i1) c **** ****** [6;6]=[6;1]x[1;6] do jloc=i1,i6 do iloc=i1,i6 TRAV66(iloc,jloc) = Ddgtn1(iloc,jloc) + Ddgtn2(iloc,jloc) . + Ddgcn3(iloc,jloc) + Ddgcn4(iloc,jloc) c [6;6] TRAV66(iloc,jloc) = TRAV66(iloc,jloc)/Denom c [6;6] end do end do do jloc=i1,i6 do iloc=i1,i6 H66t(iloc,jloc) = Da(iloc,jloc) - TRAV66(iloc,jloc) end do end do c END OF CALCULATION OF EFFECTIVE TANGENT MATRIX c ************************************************** endif c END OF PLASTICITY c store the mechanical strain at the end of this iteration c mechanical strain = epsmec + epspl do iloc=i1,i6 parahot3(idimpara3-19+iloc) = EPSMEC(iloc) + . parahot3(idimpara3-25+iloc) end do RETURN END