C SJONL2 SOURCE PASCAL 21/09/24 21:15:01 11108 SUBROUTINE SJONL2(sigi, depst, vari, xmat, trac, jtrac , > tras, jtras, trat, jtrat, sigf, varf, > depsp, ireturn) c --- SJONL2 --- c computing plastic flow in joint element c under 2D plane (strain) conditions c c former version 13-06-00 ylp. c update 15-03-01 ylp. c a shear unloading parameter has been added --> beta * 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 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(2), ddepsp(2) dimension xmat(*) c xmat 1 --> Ks elastic shear modulus c 2 --> Kt elastic traction modulus c 3 --> rho? mass per unit volume c 4 --> alpha? thermal coefficient c 5 --> edge location of the edge of the friction cone c 6 --> cplg make damage dependencies c 7 --> beta control shear unloading dimension trac(2,jtrac), tras(2,jtras), trat(2,jtrat) c 1 --> strain c 2 --> stress c WARNING original version c dimension TRAS(2,NTRAT), TRAT(2,NTRAS) ????? c_______________________________________________________________________ c * save iloop * data iloop/0/ * iloop = iloop + 1 c do i=1,2 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 c look at xmat(6) and deduce elastic/plastic... c 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 Icompr = 4 Istrain = 1 Istress = 2 c unloading constitutive parameters c the unloading shear stiffness is computed by using a linear c interpolation between the elastic stiffness and the c damaged stiffness... Beta is equal to 0 to obtain a c simple elastic unloading and is equal to 1 to get the c damaged stiffness Beta = xmat(7) 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(Ipress) = depst(Ipress)/INCR ddepst(Ishear) = depst(Ishear)/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 c [1 c if ( JTP .ne. 0 ) then Aipt = vari(Ipress) c \--> initial plastic tensile strain c print *,'initial plastic tensile strain',Aipt 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 check if Aipt has a positive value or zero c [2 if ( ireturn .ne. 0 ) return c 2] c \--> ERROR c print *,'tensile stress RTS', RTS c endif c 1] c [1 c if ( JSP .ne. 0 ) then Aips = vari(Ishear) c \--> initial plastic shear strain c print *,'initial plastic shear strain', Aips c compute the initial cohesion residual c print *,tras(Istress,1) call yofxcu(Aips, tras, jtras, RSS, dRSS, ireturn) c [2 if (ireturn .ne. 0) return c \--> ERROR c 2] c print *,'initial cohesion residual',RSS c endif c 1] c if ( JCP .eq. 0 ) then Aipc = vari(Icompr) c \--> initial plastic compression strain c compute the initial compressive strength call yofxcu(Aipc, trac, jtrac, RCS, dRCS, ireturn) c print *,'Aipc RCS ',Aipc, RCS c [2 if (ireturn .ne. 0) return c \--> ERROR c 2] c endif c 1] c______________________________________________________________________ c 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 Aitt = vari(Itotal) c \--> initial total tensile strain c compute final total tensile strain... Aftt = Aitt + ddepst(Ipress) * write(6,*) 'Aitt, Aftt =',Aitt, Aftt 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 ** mise 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' * 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 [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 damage 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 ) c \--> damaged stiffness 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 damage 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... c print *,ubound(trac) 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 Ksd = RSS / (Aips + TS0 / Ks0) c [3 if ( RSS .eq. TS0 ) then c \--> still in the elastic regime Ks = Ksd else c \--> non linear regime... c the stiffness is computed by using an c interpolation between the elastic stiffness c and the damaged stiffness Ks = ( ( 1.0D0 - Beta ) * Ks0 ) + ( Beta * Ksd ) endif c 3] endif c update the shear damage variable if (JST .eq. 0) Dshear = RSS / TS0 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(Ipress) + ( ddepst(Ipress) * Kt * Dshear) endif c 2] Tfs = sigi(Ishear) + ( ddepst(Ishear) * Ks * Dtensi) 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(Ipress) + ( ddepst(Ipress) * Kc ) endif c 2] Tfs = sigi(Ishear) + ( ddepst(Ishear) * Ks * Dcompr) endif c 1] c print *,'stress state elastic prediction' c print *,'tension compression shear' c print *,Tft, Tfc, Tfs 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 c 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(Tfs) .gt. (RSS*Chom) ) then c \--> YIELDING IN SHEAR c print *,'YIELDING IN SHEAR' c compute the total plastic shear strain... Afps = Aips + ((abs(Tfs)/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) c print *,'updated shear stress ',Tfs endif c 5] else c \--> already got some damage due to tension Afpstmp = Aips+ (abs(Tfs)/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(Tfs) .gt. (Chom*RSSTMP) ) then c \--> YIELDING IN SHEAR Afps = Afpstmp c print *,'YIELDING IN SHEAR 2' RSS = RSSTMP Tfs = sign(Chom*RSS , Tfs) c print *,'updated shear stress ',Tfs 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(Tfs) .gt. (RSS*Chom) ) then c \--> YIELDING IN SHEAR c print *,'YIELDING IN SHEAR' c compute the total plastic shear strain... Afps = Aips + ((abs(Tfs)/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) c print *,'updated shear stress ',Tfs endif c 5] else c \--> already got some damage due to compression Afpstmp = Aips+ (abs(Tfs)/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(Tfs) .gt. (Chom*RSSTMP) ) then c \--> YIELDING IN SHEAR Afps = Afpstmp c print *,'YIELDING IN SHEAR 2' RSS = RSSTMP Tfs = sign(Chom*RSS , Tfs) c print *,'updated shear stress ',Tfs 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 ! c 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(Ipress) = Tft sigi(Ipress) = Tft else sigf(Ipress) = Tfc sigi(Ipress) = Tfc endif c 1] sigf(Ishear) = Tfs sigi(Ishear) = Tfs 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(Ipress) = varf(Ipress) - vari(Ipress) else ddepsp(Ipress) = varf(Icompr) - vari(Icompr) endif c 1] ddepsp(Ishear) = varf(Ishear) - vari(Ishear) vari(Ipress) = Afpt vari(Ishear) = Afps vari(Icompr) = Afpc vari(Itotal) = Aftt depsp(Ipress) = depsp(Ipress) + ddepsp(Ipress) depsp(Ishear) = depsp(Ishear) + ddepsp(Ishear) c END OF MAIN LOOP c 0] enddo c print *,'output variables :' c print *,' --> stress state - normal shear' c print *, sigf(Ipress), sigf(Ishear) c print *,'internal variables - normal - shear - total' c print *, varf(Ipress), varf(Ishear), varf(Itotal) c print *,'plastic strain - normal - shear' c print *, depsp(Ipress), depsp(Ishear) c And, that's all folks !!!! c print *,'sjonl2 ended normaly...' end