C SJONL3 SOURCE CHAT 06/06/01 21:20:08 5450 SUBROUTINE SJONL3(sigi, depst, vari, xmat, tras, jtras , > trat, jtrat, trac, jtrac, sigf, varf, > depsp, ireturn) c --- SJONL3 --- c computing plastic flow in joint element c under 3D plane (strain) conditions c c beat release 24-7-00 ylp. c sigi IN initial stress c depst IN total strain increment c vari IN initial internal variable c xmat IN material properties c tras IN shear plastic softening behaviour c jtras IN number of points c trat IN tensile plastic softening behaviour c jtrat IN number of points c trac IN compression plastic behaviour c jtrac IN number of points c sigf OUT final stress c varf OUT final internal variable c depsp OUT plastic strain increment c ireturn OUT routine control integer IMPLICIT REAL*8(A-H,K-Z) IMPLICIT INTEGER(I-J) dimension sigi(*), depst(*), vari(*), > sigf(*), depsp(*), varf(*) dimension ddepst(3), ddepsp(3) dimension xmat(*) dimension tras(2,jtras), trat(2,jtrat), trac(2,jtrac) c 1 --> strain c 2 --> stress c_______________________________________________________________________ c * save iloop * data iloop/0/ c print *,'debug -- sjonl3d -- check everything' c c do i=1,3 c print *,'depst ',depst(i) c print *,'vari ', vari(i) c print *,'sigi ', sigi(i) c print *,'sigf ', sigf(i) c print *,'varf ',varf(i) c print *,'depst ',depst(i) c enddo c do i=1,6 c print *,'xmat ', xmat(i) c enddo c print *,'tras ', tras ,' jtras ', jtras c print *,'trat ', trat ,' jtrat ', jtrat c print *,'trac ', trac ,' jtrac ', jtrac c c print *,'ireturn ', ireturn * iloop = iloop + 1 c look at xmat(6) and deduce elastic/ plastic... c ...coupled/uncoupled variables. c all variables Jxxx are integer used as logical with the c following rule (0 is true, 1 is false) c JCP : use plastic law in compression c JSP : use plastic law in shear c JTP : use plastic law in tension c JST : couple shear and tension c JSC : couple shear and compression call sjoncpl(xmat(6),JCP,JSP,JTP,JST,JSC) * print *,'DEBUG --> entree dans sjonl2' * print *,'DepsN DepsS ' ,depst(2), depst(1) * set some usefull variables ZERO = 0.0d0 ONE = 1.0d0 c set the initial location of the edge of the cone ETi = xmat(5) c index definition Ipress = 2 Ishear = 1 Itotal = 3 Ishea1 = 1 Ishea2 = 2 Inorm = 3 Icompr = 4 Istrain = 1 Istress = 2 c number of increments INCR = 5 c______________________________________________________________________ c c elastic initialisation of the joint stiffness --> K c K[z] c with z = t tensile ; s shear ; c compression Kt0 = xmat(Ipress) Ks0 = xmat(Ishear) Kc0 = Kt0 c print *,'Deps - tension ',depst(Ipress) c print *,'Deps - shear ',depst(Ishear) c split the displacement increment in few smaller increments ddepst(Inorm) = depst(Inorm)/INCR ddepst(Ishea1) = depst(Ishea1)/INCR ddepst(Ishea2) = depst(Ishea2)/INCR c STARTING MAIN LOOP c [0 do iflag =1 , INCR c______________________________________________________________________ c c initial residual strength c RTS --> tensile c RSS --> shear c c internal variables initialisation --> Alpha c A[xyz] c with x = i initial ; f final ; k intermediate computational value c y = p plastic ; t total c z = t tensile ; s shear Aipt = vari(Ipress) c \--> initial plastic tensile strain Aips = vari(Ishear) c \--> initial plastic shear strain Aipc = vari(Icompr) c \--> initial plastic compression strain c______________________________________________________________________ c c initial residual strength c RTS --> tensile c RSS --> shear c RCS --> compression c compute the tensile stress RTS corresponding to the initial c plastic tensile strain with evolution call yofxcu(Aipt ,trat, jtrat, RTS, dRTS, ireturn) c \--> Aipt initial plastic tensile strain c RTS OUT c dRTS = dTkt/dAkps useless value c [1 if (ireturn .ne. 0) return c \--> ERROR c 1] c compute the initial cohesion residual call yofxcu(Aips, tras, jtras, RSS, dRSS, ireturn) c [1 if (ireturn .ne. 0) return c \--> ERROR c 1] c comput the initial compression residual call yofxcu(Aipc, trac, jtrac, RCS, dRCS, ireturn) c [1 if (ireturn .ne. 0) return c \--> ERROR c 1] c______________________________________________________________________ c Aitt = vari(Itotal) c \--> initial total tensile strain c compute final total tensile strain... Aftt = Aitt + ddepst(Inorm) c ... and the initial tensile strength... TT0 = trat(Istress,1) c ... and the elastic compressive strength TC0 = trac(Istress,1) c print *,'TC0 ',TC0 c check one or two things if entering sjonl2 for first time c [1 * mis en commentaire par TC * if ( iloop .eq. 1) then c [2 * if ( TT0 .eq. ZERO ) then * print *,'cut-off tensile strength is zero' * print *,'denied option' * return * endif c 2] c [2 * if ( TT0 .gt. ETi ) then * print *,'WARNING cut-off ouside the friction cone...' * print *,'... tensile strength at the edge of the cone' * print *,'set to cut-off tensile strength' * return * endif c 2] * endif * fin de mise en commentaire par TC c 1] c initialization of the damage variables Dtensi = 1.0d0 Dshear = 1.0d0 Dcompr = 1.0d0 c WARNING the damage variables are understand here as : c D = 1 means no damage c D = 0 means full damage... c which is exactly the opposite of the usual definition !!! c compute the predicted elastic stress state c using the secant stiffness c c first, look at the tension modulus... c maximal (initial) tensile stress initialisation c N.B.: it is thus assumed that the behaviour in tension is elastic c plastic, and that the peak value features the end of the elastic c regime. c check whether tension or compression region is reached c [1 if ( Aftt .gt. ZERO ) then c \--> tensile region reached... c [2 if ( JTP .ne. 0 ) then c \--> use elastic law Kt = Kt0 else c \--> use damaged elastic predictor c [3 if ( RTS .eq. ZERO ) then c \--> zero residual tensile strength Kt = ZERO c \--> tensile modulus set to zero else c \--> still got residual tensile strength Kt = RTS / ( Aipt + TT0 / Kt0 ) endif c 3] c update tensile damage variable to enable coupling effect c [3 if ( JST .eq. 0) Dtensi = RTS / TT0 c 3] endif c 2] c print *,'Dtensi ' ,Dtensi JCT = 0 c JCT is an integer value used as a logical to control the location c of the stress state with respect to tension or compression: c tension --> JCT = 0 c compression --> JCT = 1 else c \--> compression region reached... c [2 if ( JCP .ne. 0 ) then c \--> use elastic law Kc = Kc0 else c \--> use damaged elastic predictor c [3 if ( RCS .eq. ZERO ) then Kc = ZERO else Kc = RCS / (Aipc + TC0 / Kc0 ) endif c 3] c update compression damage variable c Dcompr = Kc / Kc0 --> too fast c ... try something else... endif c 2] JCT = 1 endif c ]1 c ... then, look at the shear modulus TS0 = tras(Istress,1) c [1 if (JSP .ne. 0 ) then Ks = Ks0 else c [2 if ( RSS .eq. ZERO ) then c \--> zero residual shear strength c or opened bonds Ks = ZERO c \--> shear modulus set to zero c Such case shall only appear if the bonds are opened. c otherwise, residual cohesion is generally observed else c \--> still got residual shear strength Ks = RSS / (Aips + TS0 / Ks0) endif c update the shear damage variable if (JST .eq. 0) Dshear = RSS / TS0 c(Ks / Ks0) c print *,'Dshear ' , Dshear c 2] endif c 1] c [1 if (JCT .eq. 0) then c check if previous stress state was in compression c [2 if ( (Aitt * Aftt) .lt. ZERO ) then Tft = Kt * Aftt * Dtensi else Tft = sigi(Inorm) + ( ddepst(Inorm) * Kt * Dshear) endif c 2] Tfs1 = sigi(Ishea1) + ( ddepst(Ishea1) * Ks * Dtensi) Tfs2 = sigi(Ishea2) + ( ddepst(Ishea2) * Ks * Dtensi) Tfs0 = sqrt((Tfs1**2)+(Tfs2**2)) else c check if previous stress state was in tension c [2 if ( (Aitt * Aftt) .lt. ZERO ) then Tfc = Kc * Aftt * Dcompr else c Tfc = sigi(Ipress) + ( ddepst(Ipress) * Kc * Dshear) Tfc = sigi(Inorm) + ( ddepst(Inorm) * Kc ) endif c 2] Tfs1 = sigi(Ishea1) + ( ddepst(Ishea1) * Ks * Dcompr) Tfs2 = sigi(Ishea2) + ( ddepst(Ishea2) * Ks * Dcompr) Tfs0 = sqrt((Tfs1**2)+(Tfs2**2)) endif c 1] c print *,'stress state elastic prediction' c print *,'tension compression shear' c print *,Tft, Tfc, Tfs1, Tfs2 c______________________________________________________________________ c c before checking if any kind of yielding occured, update the c final internal variables in the case of simple elastic regime. Afpt = Aipt Afps = Aips Afpc = Aipc c______________________________________________________________________ c c now, check yielding c N.B.: note that each mode might be activated independently c without any restriction!!!! c first, look at the tension/compression mode... c [1 if ( JCT .eq. 0 ) then c [2 if ( JTP .eq. 0 ) then c [3 if ( abs(Tft) .gt. RTS*Dshear ) then c \--> YIELDING IN TENSION c print *,'yielding in tension...' c Afpt = Aftt - TT0 / Kt0 Afpt = Aipt + (Tft-RTS*Dshear)/Kt c print *,'Afpt ',Afpt c \--> updated plastic tensile strain * print *,'updated plastic tensile strain ', Afpt c compute the updated residual tensile strength... call yofxcu(Afpt, trat, jtrat, RTS ,dRTS, ireturn) c \--> Afpt updated plastic tensile strain c RTS updated residual tensile strength c dTkt dummy Tft = RTS * Dshear c \--> compute the final tensile stress c print *,'residual tensile strength ',Tft endif c 3] endif c 2] else c [2 if ( JCP .eq. 0 ) then c if ( abs(Tfc) .gt. abs(RCS)*Dshear) then c [3 if ( abs(Tfc) .gt. abs(RCS) ) then c \--> YIELDING IN COMPRESSION c print *,'yielding in compression' c Afpc = abs(Aftt) - TC0 / Kc0 c Afpc = Aipc + ((Tfc*-1.0d0)-RCS*Dshear)/Kc Afpc = Aipc + ((Tfc*(-1.0d0))-RCS)/Kc c \--> updated plastic compression strain c print *,'updated plastic compression strain ', Afpc c compute the updated residual compression strength call yofxcu(Afpc, trac, jtrac, RCS , dRCS, ireturn) c then, go back to the compression domain (negative stress) Tfc = RCS * (-1.0d0) c print *,'residual compressive strength ',RCS endif c 3] endif c 2] endif c 1] c ... finally, look at the shear mode... c compute the tensile value at the edge using an homothetic transform. c the cut-off value plays the role of a drag-stress c ETf = RTS + ETi - TT0 c ... definitely not relevant!!!! c ... try something else... like a simple homothetic rule related to c the location of the residual cut-off with regard to the initial c location of the edge of the friction cone ETf = ETi * RTS / TT0 c the behaviour in shear is assumed to follow a c homothetic rule with respect to the behaviour in pure shear c the coefficient of this homothetic transform is given by... c [1 if ( JSP .eq. 0 ) then c [2 if ( ETf .ne. ZERO ) then c [3 if ( JCT .eq. 0 ) then c \--> shear yielding with tension Chom = (ETf - Tft)/(ETf)*Dtensi c print *,'Coef. hom. ', Chom c WARNING the original version used c Chom = (RTS - Tft)/(RTS +Pnor) c xmat(5) or PNOR is a simple numerical value used is order to c avoid division by zero if the residual cohesion vanishes. c print *,'???', Dtensi, ONE c [4 if ( Dtensi .eq. ONE) then c \--> no tensile damage c [5 if ( abs(Tfs0) .gt. (RSS*Chom) ) then c \--> YIELDING IN SHEAR c print *,'YIELDING IN SHEAR' c compute the total plastic shear strain... Afps = Aips + ((abs(Tfs0)/Chom)-RSS)/Ks c print *,'plasti shear strain ', Afps c ...and the final cohesion residual call yofxcu(Afps,tras,jtras,RSS,dRSS,ireturn) c [6 if (ireturn .ne. 0) return c \--> ERROR c 6] Tfs = sign(Chom*RSS , Tfs) if (Tfs0 .eq. ZERO) then Tfs1 = ZERO Tfs2 = ZERO else Tfs1= Tfs*Tfs1/Tfs0 Tfs2= Tfs*Tfs2/Tfs0 endif c print *,'updated shear stress ',Tfs c print *,'Tfs1 Tfs2 ',Tfs1,Tfs2 endif c 5] else c \--> already got some damage due to tension Afpstmp = Aips+ (abs(Tfs0)/Chom-RSS)/(Ks*Dtensi) c print *, Afpstmp c [5 if ( Afpstmp .gt. ZERO ) then call yofxcu(Afpstmp,tras,jtras,RSSTMP,dRSS,ireturn) c [6 if (ireturn .ne. 0) return c \--> ERROR c 6] c print *,'RSSTMP', RSSTMP c [6 if ( abs(Tfs0) .gt. (Chom*RSSTMP) ) then c \--> YIELDING IN SHEAR Afps = Afpstmp print *,'YIELDING IN SHEAR 2' RSS = RSSTMP Tfs = sign(Chom*RSS , Tfs) if (Tfs0 .eq. ZERO) then Tfs1 = ZERO Tfs2 = ZERO else Tfs1= Tfs*Tfs1/Tfs0 Tfs2= Tfs*Tfs2/Tfs0 endif print *,'updated shear stress ',Tfs print *,'Tfs1 Tfs2 ',Tfs1,Tfs2 endif c 6] endif c 5] endif c 4] else c \--> check shear yielding in compression region Chom = (ETf - Tfc)/(ETf)*Dcompr c print *,'Coef. hom. ', Chom c print *,'???', Dcompr, ONE c [4 if ( Dcompr .eq. ONE) then c \--> no crushing damage c [5 if ( abs(Tfs0) .gt. (RSS*Chom) ) then c \--> YIELDING IN SHEAR c print *,'YIELDING IN SHEAR' c compute the total plastic shear strain... Afps = Aips + ((abs(Tfs0)/Chom)-RSS)/Ks c print *,'plasti shear strain ', Afps c ...and the final cohesion residual call yofxcu(Afps,tras,jtras,RSS,dRSS,ireturn) c [6 if (ireturn .ne. 0) return c \--> ERROR c 6] Tfs = sign(Chom*RSS , Tfs) if (Tfs0 .eq. ZERO) then Tfs1 = ZERO Tfs2 = ZERO else Tfs1= Tfs*Tfs1/Tfs0 Tfs2= Tfs*Tfs2/Tfs0 endif c print *,'updated shear stress ',Tfs c print *,'Tfs1 Tfs2 ',Tfs1,Tfs2 endif c 5] else c \--> already got some damage due to compression Afpstmp = Aips+ (abs(Tfs0)/Chom-RSS)/(Ks*Dcompr) c print *, Afpstmp c [5 if ( Afpstmp .gt. ZERO ) then call yofxcu(Afpstmp,tras,jtras,RSSTMP,dRSS,ireturn) c [6 if (ireturn .ne. 0) return c \--> ERROR c 6] c print *,'RSSTMP', RSSTMP c [6 if ( abs(Tfs0) .gt. (Chom*RSSTMP) ) then c \--> YIELDING IN SHEAR Afps = Afpstmp c print *,'YIELDING IN SHEAR 2' RSS = RSSTMP Tfs = sign(Chom*RSS , Tfs) if (Tfs0 .eq. ZERO) then Tfs1 = ZERO Tfs2 = ZERO else Tfs1= Tfs*Tfs1/Tfs0 Tfs2= Tfs*Tfs2/Tfs0 endif c print *,'updated shear stress ',Tfs c print *,'Tfs1 Tfs2 ',Tfs1,Tfs2 endif c 6] endif c 5] endif c 4] endif c 3] else c \--> the edge has reach the origin... no homothetic rule c can be computed ! print *,'EDGE ZERO -- ERROR' endif c 2] endif c 1] c______________________________________________________________________ c c final stage --> updating the output variables... c ... first the final stress... c [1 if ( JCT .eq. 0 ) then sigf(Inorm) = Tft sigi(Inorm) = Tft else sigf(Inorm) = Tfc sigi(Inorm) = Tfc endif c 1] sigf(Ishea1) = Tfs1 sigi(Ishea1) = Tfs1 sigf(Ishea2) = Tfs2 sigi(Ishea2) = Tfs2 c ... then the final internal variables... varf(Ipress) = Afpt varf(Ishear) = Afps varf(Icompr) = Afpc varf(Itotal) = Aftt c ... and finally, the plastic strain increments. c [1 if ( JCT .eq. 0 ) then ddepsp(Inorm) = varf(Ipress) - vari(Ipress) else ddepsp(Inorm) = varf(Icompr) - vari(Icompr) endif c 1] ddepsp(Ishea1) = varf(Ishear) - vari(Ishear) ddepsp(Ishea2) = varf(Ishear) - vari(Ishear) vari(Ipress) = Afpt vari(Ishear) = Afps vari(Icompr) = Afpc vari(Itotal) = Aftt depsp(Inorm) = depsp(Inorm) + ddepsp(Inorm) depsp(Ishea1) = depsp(Ishea1) + ddepsp(Ishea1) depsp(Ishea2) = depsp(Ishea2) + ddepsp(Ishea2) c END OF MAIN LOOP c 0] enddo c print *,'output variables :' c print *,' --> stress state - normal shear' c print *, sigf(Inorm), sigf(Ishea1), sigf(Ishea2) c print *,'internal variables - normal - comp - shear - total' c print *, varf(Ipress), varf(Icompr),varf(Ishear), varf(Itotal) c print *,'plastic strain - normal - shear' c print *, depsp(Inorm), depsp(Ishea1),depsp(Ishea2) c And, that's all folks !!!! c print *,'sjonl2 ended normaly...' end