FE simulation of traction separation relations

Question

I am trying to move completely open source. One of the things my research group does is cohesive zone modelling using traction separation laws, which I currently implement in ABAQUS. How easy would it be to write my own code for this? Our geometry is very simple, just two beams which we assume are linear elastic. Are there any open source libraries which can do the same?

Edit: I am adding the actual PDEs that are involved. This is basically a linear elasticity problem $\nabla.\tau=0$ in $\Omega$.

$\tau =\mathbf{C}\epsilon$ in $\Omega$.

$\epsilon=\frac{1}{2}\left(\nabla \vec{u}+\nabla{\vec{u}^T}\right)$ in $\Omega$.

$t=\tau \vec{n}$ on $\Gamma$

$t=\bar{t}$ on $\Gamma_t$

$\vec{u}=\vec{\bar{u}} $ on $\Gamma_u$

In addition to this, we have an extra condition

$t= f(\vec{u})$ on $\Gamma_c$

The overall boundary is composed of $\Gamma_u,\Gamma_t,\Gamma_c$. The main difficulty associated with this is the nature of the function $f$. It first increases with $\vec{u}$ and then drops to zero (in fact the most common traction-separation law is a bilinear). Once the traction falls to zero, that part of $\Gamma_c$ is essentially a part of $\Gamma_t$ with homogeneous BC. As far as I understand, the difficulty is with implementing the softening part of the cohesive law (where $f$ drops). Just attaching the ABAQUS manual for this to show what functionality is desired. http://130.149.89.49:2080/v6.8/books/usb/default.htm?startat=pt09ch31s01alm62.html

Answer 1

Cohesive zone couplings are very straightforward to implement in many modern open source finite element libraries.

For example in GetFEM you can use this code as a starting point:

import getfem as gf
gf.util_trace_level(1)

# Input data
NX = 60         # number of elements in horizontal direction
NY = 60         # number of elements in vertical direction
LX = 1.         # [mm] Length
LY = 1.         # [mm] Height

E = 210e3       # [N/mm^2]
nu = 0.3

tN = 4000 # [N/mm^2]
tS = 4000 # [N/mm^2]
delta0 = 0.001  # [mm]
deltaf = 0.003  # [mm]

eps_max = 0.01
steps = 30

disp_fem_order = 2       # displacements finite element order

# mesh generation
mesh = gf.Mesh("import", "structured",
               f"GT='GT_QK(2,2)';ORG=[{-LX/2:e},0];SIZES=[{LX:e},{LY:e}];NSUBDIV=[{NX},{NY}]")
# mesh regions
T_RG = 6 # top
BR_RG = 8 # bottom right
mesh.set_region(T_RG, mesh.outer_faces_in_box([-LX/2-1e-5,LY-1e-5],[LX/2+1e-5,LY+1e-5]))
mesh.set_region(BR_RG, mesh.outer_faces_in_box([-1e-5,-1e-5],[LX/2+1e-5,1e-5]))

# FEM
mfu = gf.MeshFem(mesh, 2)
mfu.set_classical_fem(disp_fem_order)
mfdir = mfu

mfout = gf.MeshFem(mesh)
mfout.set_classical_discontinuous_fem(2)

# Integration methods
mim9 = gf.MeshIm(mesh, 5)

# Model
md = gf.Model("real")
md.add_fem_variable("u", mfu)      # displacements field

md.add_initialized_data("K", E/(3.*(1.-2.*nu)))
md.add_initialized_data("G", E/(2*(1+nu)))
md.add_initialized_data("tN0", tN)
md.add_initialized_data("tS0", tS)
md.add_initialized_data("r", delta0/deltaf)
md.add_initialized_data("deltaf", deltaf)
md.add_macro("deveps", "Sym(Grad(u))-Div(u)/3*Id(2)")
md.add_linear_term(mim9, "(K*Div(u)*Id(2)+2*G*deveps):Grad(Test_u)")
md.add_macro("uN", "pos_part(min(1,u(2)/deltaf))")
md.add_macro("uS", "max(-1,min(1,u(1)/deltaf))")
md.add_macro("tN", "(uN-pos_part(uN-r)/(1-r))/r*tN0")
md.add_macro("tS", "(uS+(pos_part(uS-r)-pos_part(r-uS))/(1-r))/r*tS0")
md.add_nonlinear_term(mim9, "[tS;tN-neg_part(u(2))*1e4].Test_u", BR_RG)

# loading and boundary conditions
md.add_fem_data("dirichlet_data", mfu);
ibdir = md.add_Dirichlet_condition_with_multipliers(mim9, "u", mfdir, T_RG, "dirichlet_data")
dirmultname = md.mult_varname_Dirichlet(ibdir)

print("Displacement dofs: %i\nTotal model dofs: %i" % (mfu.nbdof(), md.nbdof()))
with open("gf_cohesive_zone_forces.dat", "w") as f:
  for step in range(steps):
    eps = eps_max*step/(steps-1.)
    print(f"{step}: Solve for applied average strain {eps:.4e}")
    md.set_variable("dirichlet_data", md.interpolation(f"{eps:.15g}*[0;X(2)]", mfdir))
    nit, conv =  md.solve("noisy", "lsolver", "mumps", "max_iter", 20, "max_res", 1e-9,
                          "lsearch", "simplest", "alpha max ratio", 1e9, "alpha min", 1.,
                          "alpha mult", 0.1, "alpha threshold res", 1e9)
    out = (mfu, md.variable("u"), "Displacements")
    for i,j in [[1,1],[2,2],[1,2]]:
      out += (mfout, md.interpolation(f"(K*Div(u)*Id(2)+2*G*deveps)({i},{j})", mfout),
              f"Cauchy Stress {i}{j}")
    mfout.export_to_vtu(f"gf_cohesive_zone_{step}.vtu", *out)
    f.write("step=%i average strain=%e reaction force=%e\n" %
            (step, eps, gf.asm_generic(mim9, 0, dirmultname+"(2)", T_RG, md)))
    f.flush()