! User element subroutine for the phase field model for fracture

! Quasi-Newton, suitable for fatigue and static fracture

! The code is distributed under a BSD license? ? ?

? ? ??

! If using this code for research or industrial purposes, please cite:

! P.K. Kristensen and E. Martínez-Pa?eda. Phase field fracture modelling

! using quasi-Newton and a new adaptive step scheme.

! Theoretical and Applied Fracture Mechanics 107, 102446 (2020)

! doi: 10.1016/j.tafmec.2019.102446

? ? ??

! Emilio Martínez-Pa?eda (e.martinez-paneda@imperial.ac.uk)

! Imperial College London

? ? ? module kvisual

? ? ? implicit none

? ? ? real*8 UserVar(70000,11,4),Ac

? ? ? integer nelem

? ? ? save

? ? ? end module

? ? ? subroutine uel(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars,

? ? ?1 props,nprops,coords,mcrd,nnode,u,du,v,a,jtype,time,dtime,

? ? ?2 kstep,kinc,jelem,params,ndload,jdltyp,adlmag,predef,npredf,

? ? ?3 lflags,mlvarx,ddlmag,mdload,pnewdt,jprops,njpro,period)

? ? ? use kvisual

? ? ? include 'aba_param.inc' !implicit real(a-h o-z)

? ? ??

? ? ? dimension rhs(mlvarx,*),amatrx(ndofel,ndofel),props(*),svars(*),

? ? ?1 energy(*),coords(mcrd,nnode),u(ndofel),du(mlvarx,*),v(ndofel),

? ? ?2 a(ndofel),time(2),params(*),jdltyp(mdload,*),adlmag(mdload,*),

? ? ?3 ddlmag(mdload,*),predef(2,npredf,nnode),lflags(*),jprops(*)

? ? ? parameter(ndim=2,ntens=4,ninpt=4,nsvint=12)

? ? ??

? ? ? dimension wght(ninpt),dN(nnode,1),dNdz(ndim,nnode),stran(ntens),

? ? ?2 dNdx(ndim,nnode),b(ntens,nnode*ndim),ddsdde(ntens,ntens),

? ? ?3 stress(ntens),dstran(ntens),statevLocal(nsvint)

? ? ??

? ? ? data wght /1.d0, 1.d0, 1.d0, 1.d0/

? ? ? ??

!? ? ?initialising

? ? ? do k1=1,ndofel

? ? ? ?rhs(k1,1)=0.d0

? ? ? end do

? ? ? amatrx=0.d0

? ? ??

!? ? ?find number of elements? ? ? ? ??

? ? ? if (dtime.eq.0.d0) then

? ? ? ?if (jelem.eq.1) then

? ? ? ? nelem=jelem

? ? ? ?else

? ? ? ? if (jelem.gt.nelem) nelem=jelem?

? ? ? ?endif?

? ? ? endif? ? ??

? ? ? if (jelem.eq.1) then

? ? ? Ac=0.d0

? ? ? endif

!? ? ?reading parameters

? ? ? xlc=props(3)

? ? ? Gc=props(4)

? ? ? xk=props(5)

? ? ? kFlagF=props(6)

? ? ? if (kFlagF.eq.1) then

? ? ? ? alphT=Gc/(2.d0*6.d0*xlc)

? ? ? else

? ? ? ? alphT=1.d10

? ? ? endif

? ? ? do kintk=1,ninpt

!? ? ?evaluate shape functions and derivatives

? ? ? ?call kshapefcn(kintk,ninpt,nnode,ndim,dN,dNdz)? ? ??

? ? ? ?call kjacobian(jelem,ndim,nnode,coords,dNdz,djac,dNdx,mcrd)

? ? ? ?dvol=wght(kintk)*djac

? ? ? ?

!? ? ?form B-matrix

? ? ? ?b=0.d0

? ? ? ?do inod=1,nnode

? ? ? ? b(1,2*inod-1)=dNdx(1,inod)

? ? ? ? b(2,2*inod)=dNdx(2,inod)

? ? ? ? b(4,2*inod-1)=dNdx(2,inod)

? ? ? ? b(4,2*inod)=dNdx(1,inod)

? ? ? ?end do? ? ? ? ? ? ? ? ? ? ?

? ? ? ?

!? ? ?compute from nodal values

? ? ? ?phi=0.d0

? ? ? ?do inod=1,nnode

? ? ? ? phi=phi+dN(inod,1)*u(ndim*nnode+inod)

? ? ? ?end do? ?

? ? ? ?if (phi.gt.1.d0) phi=1.d0

? ? ? ? ? ?

!? ? ?compute the increment of strain and recover history variables

? ? ? ?dstran=matmul(b,du(1:ndim*nnode,1))

? ? ? ?

? ? ? ?call kstatevar(kintk,nsvint,svars,statevLocal,1)?

? ? ? ?

? ? ? ?stress=statevLocal(1:ntens)

? ? ? ?stran(1:ntens)=statevLocal((ntens+1):(2*ntens))

? ? ? ?phin=statevLocal(2*ntens+1)

? ? ? ?Hn=statevLocal(2*ntens+2)

? ? ? ?alphn=statevLocal(2*ntens+3)

? ? ? ?alphBn=statevLocal(2*ntens+4)

? ? ? ?if (dtime.eq.0.d0) phin=phi

? ? ??

!? ? ?update crack length

? ? ? ?coordx=0.d0

? ? ? ?if (phi.ge.0.95d0) then

? ? ? ? do i=1,nnode

? ? ? ? ?coordx=coordx+dN(i,1)*coords(1,i)

? ? ? ? enddo

? ? ? ? if (coordx.gt.Ac) then

? ? ? ? ?Ac=coordx

? ? ? ? endif

? ? ? ?endif

!? ? ?call umat to obtain stresses and constitutive matrix?

? ? ? ?call kumat(props,ddsdde,stress,dstran,ntens,statevLocal)

? ? ? ?stran=stran+dstran

? ? ? ?

?!? ? ?compute strain energy density? ? ??

? ? ? ?Psi=0.d0

? ? ? ?do k1=1,ntens

? ? ? ? Psi=Psi+stress(k1)*stran(k1)*0.5d0

? ? ? ?end do

? ? ? ?alph=Psi*((1-phi)**2+xk)

!? ? ?Update fatigue history variable

? ? ? ?if ((alph.ge.alphn).and.(dtime.gt.0.d0)) then

? ? ? ? alphB = alphBn+abs(alph-alphn)

? ? ? ?else

? ? ? ? alphB=alphBn

? ? ? ?endif? ? ??

? ? ? ? ? ?

!? ? ?enforcing Karush-Kuhn-Tucker conditions

? ? ? ?if (Psi.gt.Hn) then

? ? ? ? H=Psi

? ? ? ?else

? ? ? ? H=Hn

? ? ? ?endif

? ? ? ?if (alphB.lt.alphT) then

? ? ? ? Fdeg= 1.d0

? ? ? ?else

? ? ? ? Fdeg=(2.d0*alphT/(alphB+alphT))**2.d0

? ? ? ?endif

? ? ? ?statevLocal(1:ntens)=stress(1:ntens)

? ? ? ?statevLocal((ntens+1):(2*ntens))=stran(1:ntens)

? ? ? ?statevLocal(2*ntens+1)=phi

? ? ? ?statevLocal(2*ntens+2)=H

? ? ? ?statevLocal(2*ntens+3)=alph

? ? ? ?statevLocal(2*ntens+4)=alphB? ? ?

? ? ? ?call kstatevar(kintk,nsvint,svars,statevLocal,0)

? ? ? ?

? ? ? ?amatrx(1:8,1:8)=amatrx(1:8,1:8)+

? ? ?1 dvol*(((1.d0-phi)**2+xk)*matmul(matmul(transpose(b),ddsdde),b))

? ? ? ??

? ? ? ?rhs(1:8,1)=rhs(1:8,1)-

? ? ?1 dvol*(matmul(transpose(b),stress)*((1.d0-phi)**2+xk))? ? ? ?

? ? ? ? ? ?

? ? ? ? amatrx(9:12,9:12)=amatrx(9:12,9:12)

? ? ?1 +dvol*(matmul(transpose(dNdx),dNdx)*Gc*xlc*Fdeg

? ? ?2 +matmul(dN,transpose(dN))*(Gc/xlc*Fdeg+2.d0*H))? ?

? ? ? ? rhs(9:12,1)=rhs(9:12,1)

? ? ?1 -dvol*(matmul(transpose(dNdx),matmul(dNdx,u(9:12)))*Fdeg

? ? ?2 *Gc*xlc+dN(:,1)*((Gc/xlc*Fdeg+2.d0*H)*phi-2.d0*H))? ? ? ?

? ? ? ? ? ? ? ? ? ?

! output

? ? ? ?UserVar(jelem,1:4,kintk)=statevLocal(1:ntens)*((1.d0-phi)**2+xk) ! Stress

? ? ? ?UserVar(jelem,5:9,kintk)=statevLocal((ntens+1):(2*ntens+1)) ! Strains and phi

? ? ? ?UserVar(jelem,10:11,kintk)=statevLocal((2*ntens+3):(2*ntens+4)) ! alpha and \bar{alpha}

? ? ??

? ? ? end do? ? ? ?! end loop on material integration points

? ? ??

? ? ? RETURN

? ? ? END

? ? ??

? ? ? subroutine kshapefcn(kintk,ninpt,nnode,ndim,dN,dNdz)

c

? ? ? include 'aba_param.inc'

c

? ? ? parameter (gaussCoord=0.577350269d0)

? ? ? dimension dN(nnode,1),dNdz(ndim,*),coord24(2,4)

? ? ??

? ? ? data? coord24 /-1.d0, -1.d0,

? ? ?2? ? ? ? ? ? ? ? 1.d0, -1.d0,

? ? ?3? ? ? ? ? ? ? ?-1.d0,? 1.d0,

? ? ?4? ? ? ? ? ? ? ? 1.d0,? 1.d0/? ? ??

? ? ?

!? ? ?2D 4-nodes

!? ? ?determine (g,h,r)

? ? ? g=coord24(1,kintk)*gaussCoord

? ? ? h=coord24(2,kintk)*gaussCoord

!? ? ?shape functions?

? ? ? dN(1,1)=(1.d0-g)*(1.d0-h)/4.d0

? ? ? dN(2,1)=(1.d0+g)*(1.d0-h)/4.d0

? ? ? dN(3,1)=(1.d0+g)*(1.d0+h)/4.d0

? ? ? dN(4,1)=(1.d0-g)*(1.d0+h)/4.d0

!? ? ?derivative d(Ni)/d(g)

? ? ? dNdz(1,1)=-(1.d0-h)/4.d0

? ? ? dNdz(1,2)=(1.d0-h)/4.d0

? ? ? dNdz(1,3)=(1.d0+h)/4.d0

? ? ? dNdz(1,4)=-(1.d0+h)/4.d0

!? ? ?derivative d(Ni)/d(h)

? ? ? dNdz(2,1)=-(1.d0-g)/4.d0

? ? ? dNdz(2,2)=-(1.d0+g)/4.d0

? ? ? dNdz(2,3)=(1.d0+g)/4.d0

? ? ? dNdz(2,4)=(1.d0-g)/4.d0

? ? ??

? ? ? return

? ? ? end?

? ? ? subroutine kjacobian(jelem,ndim,nnode,coords,dNdz,djac,dNdx,mcrd)

!? ? ?Notation: djac - Jac determinant; xjaci - inverse of Jac matrix?

!? ? ?dNdx - shape functions derivatives w.r.t. global coordinates

? ? ? include 'aba_param.inc'

? ? ? dimension xjac(ndim,ndim),xjaci(ndim,ndim),coords(mcrd,nnode),

? ? ?1 dNdz(ndim,nnode),dNdx(ndim,nnode)

? ? ? xjac=0.d0

? ? ? do inod=1,nnode

? ? ? ?do idim=1,ndim

? ? ? ? do jdim=1,ndim

? ? ? ? ?xjac(jdim,idim)=xjac(jdim,idim)+

? ? ?1? ? ? ? dNdz(jdim,inod)*coords(idim,inod)? ? ??

? ? ? ? end do

? ? ? ?end do?

? ? ? end do

? ? ? djac=xjac(1,1)*xjac(2,2)-xjac(1,2)*xjac(2,1)

? ? ? if (djac.gt.0.d0) then ! jacobian is positive - o.k.

? ? ? ?xjaci(1,1)=xjac(2,2)/djac

? ? ? ?xjaci(2,2)=xjac(1,1)/djac

? ? ? ?xjaci(1,2)=-xjac(1,2)/djac

? ? ? ?xjaci(2,1)=-xjac(2,1)/djac

? ? ? else ! negative or zero jacobian

? ? ? ?write(7,*)'WARNING: element',jelem,'has neg. Jacobian'

? ? ? endif

??

? ? ? dNdx=matmul(xjaci,dNdz)?

? ? ? return

? ? ? end

c*****************************************************************

? ? ? subroutine kstatevar(npt,nsvint,statev,statev_ip,icopy)

c

c? ? ?Transfer data to/from element-level state variable array from/to

c? ? ?material-point level state variable array.

c

? ? ? include 'aba_param.inc'

? ? ? dimension statev(*),statev_ip(*)

? ? ? isvinc=(npt-1)*nsvint? ? ?! integration point increment

? ? ? if (icopy.eq.1) then ! Prepare arrays for entry into umat

? ? ? ?do i=1,nsvint

? ? ? ? statev_ip(i)=statev(i+isvinc)

? ? ? ?enddo

? ? ? else ! Update element state variables upon return from umat

? ? ? ?do i=1,nsvint

? ? ? ? statev(i+isvinc)=statev_ip(i)

? ? ? ?enddo

? ? ? end if

? ? ? return

? ? ? end

c*****************************************************************

? ? ? subroutine kumat(props,ddsdde,stress,dstran,ntens,statev)

c

c? ? ?Subroutine with the material model

c

? ? ? include 'aba_param.inc' !implicit real(a-h o-z)

? ? ??

? ? ? dimension props(*),ddsdde(ntens,ntens),stress(ntens),statev(*),

? ? ?+ dstran(ntens)

!? ? ?Initialization

? ? ? ddsdde=0.d0

? ? ? E=props(1) ! Young's modulus

? ? ? xnu=props(2) ! Poisson's ratio

? ? ??

!? ? ?Build stiffness matrix

? ? ? eg2=E/(1.d0+xnu)

? ? ? elam=(E/(1.d0-2.d0*xnu)-eg2)/3.d0

? ? ??

!? ? ?Update stresses

? ? ? do k1=1,3

? ? ? ?do k2=1,3

? ? ? ? ddsdde(k2,k1)=elam

? ? ? ?end do

? ? ? ?ddsdde(k1,k1)=eg2+elam

? ? ? end do

? ? ? ddsdde(4,4)=eg2/2.d0

? ? ??

? ? ? stress=stress+matmul(ddsdde,dstran)? ?

? ? ? return

? ? ? end

c*****************************************************************

? ? ? subroutine umat(stress,statev,ddsdde,sse,spd,scd,rpl,ddsddt,

? ? ?1 drplde,drpldt,stran,dstran,time,dtime,temp2,dtemp,predef,dpred,

? ? ?2 cmname,ndi,nshr,ntens,nstatv,props,nprops,coords,drot,pnewdt,

? ? ?3 celent,dfgrd0,dfgrd1,noel,npt,layer,kspt,jstep,kinc)

? ? ? use kvisual

? ? ? include 'aba_param.inc' !implicit real(a-h o-z)

? ? ? character*8 cmname

? ? ? dimension stress(ntens),statev(nstatv),ddsdde(ntens,ntens),

? ? ?1 ddsddt(ntens),drplde(ntens),stran(ntens),dstran(ntens),

? ? ?2 time(2),predef(1),dpred(1),props(nprops),coords(3),drot(3,3),

? ? ?3 dfgrd0(3,3),dfgrd1(3,3),jstep(4)

? ? ? ddsdde=0.0d0

? ? ? noffset=noel-nelem

? ? ? statev(1:nstatv-1)=UserVar(noffset,1:nstatv-1,npt)

? ? ? statev(nstatv)=Ac

? ? ? return

? ? ? end