# Gridap.jl: Help with a Coupled PDE

Looking over the various FEM PDE solvers (those written in julia and other languages), few of them have tutorials that deal with fracture mechanics. While there are different formulations that could be used, I want to implement the cohesive-zone model (CZM) because it can readily handle both energy-based and strength-based crack propagation.

Assuming linear elasticity (since this is the simplest continuum model), the cohesive-zone model can be written in the following (coupled) weak form:
\int_{\Omega^+} \sigma : \varepsilon(v^+) d\Omega = \int_{\Gamma_F^+} F^+ \cdot v^+ d\Gamma + \int_{\Gamma_C^+} t^+ \cdot \delta(v^+) d\Gamma \,\,\, \forall v^+ \in U
\int_{\Omega^-} \sigma : \varepsilon(v^-) d\Omega = \int_{\Gamma_F^-} F^- \cdot v^- d\Gamma + \int_{\Gamma_C^-} t^- \cdot \delta(v^-) d\Gamma \,\,\, \forall v^- \in U ,
where \Omega^+ and \Omega^- make up the two elastic portions of \Omega. \Gamma_C splits \Omega in two and the surfaces can separate from each other. \Gamma_F^{+/-} is a boundary where forces are applied to the body (if they exist).

The two equations above are coupled by:
t^+ = -t^- = t(\delta) \, (\text{on } \Gamma_C) ,
where t(\delta) is the traction-separation law that determines the interface stresses and separation between the \Gamma_C^+ and \Gamma_C^- surfaces (\Gamma_C = \Gamma_C^+ \cup \Gamma_C^-). Importantly, t(\delta) will go to zero at some separation distance and the surfaces can separate (crack propagation).

With that in mind, I was wondering how this might be implemented in Gridap.jl. Looking at fluid-structure interaction tutorial (11 Fluid-Structure Interaction · Gridap tutorials), it appears as though all the necessary building blocks are required to solve this problem (e.g. CellField).

Does @fverdugo or anyone else in the Gridap.jl community have thoughts on how this might be implemented? It might require a new element type being written, and I have a CZM UEL (user-defined element) from Abaqus (written in Fortran) that I could convert to julia.

1 Like

Most of the pieces should be there. I think the main missing part is the computation of the distance between integration points between both sides of the crack. In any case, you can do this yourself.

Thank you for the suggestion. I will look look through the documentation to try and figure that out.

Do you think it is necessary to write a new finite element? Or might it be acceptable to just ignore the 1-element-thick crack propagation path? I wasn’t sure whether or not that might be a problem for a non-lnear solver.

Is there a way to simply provide the local stiffness matrix for each element? I am trying to implement the cohesive-zone model, which means traction-free boundaries will be generated once a crack propagates. Furthermore, a moment correction (appendix) is needed to be made for physically-accurate calculations.

If not, do you have a separate suggestion on how I might implement this traction-free boundary (and moment correction) into the standard element? I’m wondering if a new set of basis functions might be required.

Accessing the local matrices is very easy. Just evaluate your bilinear form with the test/trial finite element bases, and then index the result with a triangulation object. It will return an array with the local matrices on the cells of the given triangulation.

This is an example based on the 2nd Gridap tutorial.

Hope it helps!
Francesc

using Gridap
using Gridap.FESpaces
u(x) = x[1] + x[2]
f(x) = 0
domain = (0,1,0,1)
partition = (4,4)
model = CartesianDiscreteModel(domain,partition)
order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
V0 = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags="boundary")
U = TrialFESpace(V0,u)
degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)
a(u,v) = ∫( ∇(v)⊙∇(u) )*dΩ
b(v) = ∫( v*f )*dΩ
du = get_trial_fe_basis(U)
dv = get_fe_basis(V)
cell_mat = a(du,dv)[Ω]
cache = array_cache(cell_mat)
ncells = length(cell_mat)
for cell in 1:ncells
mat = getindex!(cache,cell_mat,cell) # This is the local cell matrix at cell 'cell'
end