Sources Esope | Cast3M

elnlin
C ELNLIN    SOURCE    CB215821  26/06/25    21:15:08     12581          
       subroutine ELNLIN (da, uel, force, XMKel,
     &             delta_prec, Ftr, delta,deltamax)
 
        IMPLICIT INTEGER(I-N)
        IMPLICIT REAL*8(A-H,O-Z)
*
        DIMENSION da(20), XMKel(5,5), uel(5), force(5), Ftr(5)
        DIMENSION Res(3), AA(3,3),AA1(3,3),DVM(3)
*
        a=da(8)
        b=da(9)
        c=da(10)
        d=da(11)
        e=da(12)
        f=da(13)
        a9=da(15)
*
        A1=2.5D0
        epsilon2=1.D-12
C       Calcul de la limite élastique en décollement.
        uelu   =sqrt(uel(3)**2+uel(5)**2)
        forceu =sqrt(force(3)**2+force(5)**2)
        XMKelu =1.D0/sqrt(2.D0)*sqrt(XMKel(3,3)**2+XMKel(5,5)**2)
*
        if (a9.eq.1) then
          q_1=4.D0
          q_2=1.D0
        else
          if (a9.eq.2) then
            q_1=4.D0
            q_2=1.D0
          else
            q_1=6.D0
            q_2=2.D0
          endif
        endif
*
        fe=forceu**2.D0-(force(1)/q_1*EXP(-A1*force(1)))**2
*
        if (fe.le.0.D0) then
        call MATVE1 (XMKel, uel, 5,5, Ftr,2)
        delta=delta_prec
        else
*
        deltamax=1.D0
*
                ff1=(force(1)*(EXP(-A1*force(1))
     & +q_2*(delta_prec*deltamax)))/q_1-forceu
                gg1=-(((delta_prec/1.D0)/(1-(delta_prec/1.D0))
     & +LOG(1.D0-(delta_prec/1.D0)))
     & *force(1)*EXP(-A1*force(1)))
     & /(2.D0*q_1*XMKelu)+force(1)/XMKel(1,1)-uel(1)
        Xhh1=((delta_prec/1.D0)**2.D0*force(1)*EXP(-A1*force(1)))/q_1
     & +(1.D0-(delta_prec/1.D0))*(forceu-XMKelu*uelu)
        Res(1)=-ff1
        Res(2)=-gg1
        Res(3)=-Xhh1
        XNorm_Res=SQRT(Res(1)**2.D0+Res(2)**2.D0+Res(3)**2.D0)
                do while (XNorm_Res.gt.epsilon2)
                dff1dV=(EXP(-A1*force(1))+q_2*(delta_prec*deltamax))/q_1
     & -(A1*force(1)*EXP(-A1*force(1)))/q_1
                dff1dM=-1.D0
                dff1ddelta=1.D0*deltamax*q_2*force(1)/q_1
                dgg1dV=(((delta_prec/1.D0)/(1.D0-(delta_prec/1.D0))
     & +LOG(1.D0-(delta_prec/1.D0)))*A1*force(1)*EXP(-A1*force(1)))
     & /(2.D0*q_1*XMKelu)
     & -(((delta_prec/1.D0)/(1.D0-(delta_prec/1.D0))
     & +LOG(1.D0-(delta_prec/1.D0)))
     & *EXP(-A1*force(1)))/(2.D0*q_1*XMKelu)+1.D0/XMKel(1,1)
                dgg1dM=0.D0
                dgg1ddelta=-1.D0/1.D0*(((delta_prec/1.D0)*force(1)
     & *EXP(-A1*force(1)))/
     &          (2.D0*q_1*(1.D0-(delta_prec/1.D0))**2.D0*XMKelu))
                dhh1dV=((delta_prec/1.D0)**2.D0*EXP(-A1*force(1)))/q_1
     & -((delta_prec/1.D0)**2.D0*A1*force(1)*EXP(-A1*force(1)))/q_1
                dhh1dM=1.D0-(delta_prec/1.D0)
                dhh1ddelta=1.D0/1.D0*(2.D0*((delta_prec/1)*force(1)
     & *EXP(-A1*force(1)))/q_1-forceu+XMKelu*uelu)
                AA(1,1)=dff1dV
                AA(1,2)=dff1dM
                AA(1,3)=dff1ddelta
                AA(2,1)=dgg1dV
                AA(2,2)=dgg1dM
                AA(2,3)=dgg1ddelta
                AA(3,1)=dhh1dV
                AA(3,2)=dhh1dM
                AA(3,3)=dhh1ddelta
                call HINV33(AA,AA1)
                call MATVE1(AA1,Res,3,3,DVM,2)
                force(1)=force(1)+DVM(1)
                forceu=forceu+DVM(2)
                delta_prec=delta_prec+DVM(3)
*
                ff1=(force(1)*(EXP(-A1*force(1))
     & +q_2*(delta_prec*deltamax)))/q_1-forceu
                gg1=-(((delta_prec/1.D0)/(1.D0-(delta_prec/1.D0))
     & +LOG(1.D0-(delta_prec/1.D0)))
     & *force(1)*EXP(-A1*force(1)))
     & /(2.D0*q_1*XMKelu)+force(1)/XMKel(1,1)-uel(1)
        Xhh1=((delta_prec/1.D0)**2.D0*force(1)*EXP(-A1*force(1)))/q_1
     & +(1.D0-(delta_prec/1.D0))*(forceu-XMKelu*uelu)
        Res(1)=-ff1
        Res(2)=-gg1
        Res(3)=-Xhh1
        XNorm_Res=SQRT(Res(1)**2+Res(2)**2+Res(3)**2)
                enddo
        force(2)=XMKel(2,2)*uel(2)
        force(4)=XMKel(2,2)*uel(4)
        force(3)=uel(3)/uelu*forceu
        force(5)=uel(5)/uelu*forceu
        uelu1=sqrt(uel(3)**2+uel(5)**2)
       DO I=1,5
        Ftr(I)=force(I)
       ENDDO
        delta=delta_prec
        endif
        return
        end