# FE simulation of traction separation relations

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

• You might want to show the actual equations you are trying to solve. A lot of people know PDEs, but few of them will know what "cohesive zone modelling using separation laws" actually means. Aug 30, 2023 at 16:47
• @WolfgangBangerth Added Aug 31, 2023 at 2:35
• @aaquib I am in the process of writing a DG code for CZM using dealii so it's possible. Aug 31, 2023 at 7:10
• If @CuteCompute already has a piece of code in deal.II, then you should contact them of course! Separately, step-8 and step-17 in the deal.II tutorial solve the linear elasticity equation, and that isn't very difficult at all. In your case, the equation is nonlinear and so you should look at step-15 for nonlinear solvers (or perhaps step-77). Sep 1, 2023 at 3:07
• I recommend searching the repositories of open source software to see what already exists. For example, here is one open issue for deal.II: github.com/dealii/dealii/issues/7153 Sep 1, 2023 at 15:29

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")

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()