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
Hi @James_Gorman
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.
Hi @James_Gorman,
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
@fverdugo,
I have been distracted with other projects, but I still want to formulate this fracture mechanics problem. As I return to looking into using Gridap.jl, there is some nuance to the FSI tutorial for boundary coupling compared to what I am wanting to do. While numerous codes ignore interfacial thickness for the cohesive-zone model, I do not want to ignore this because it doesnāt allow for compression. As such, there are ``twoāā interfaces where tractions need to be applied (although they are just equal and opposite).
With that in mind, do you have a suggestion of how to couple the two interfaces and ignore the single-layer of elements separating the two elastic bodies? It appears as though the traditional InterfaceTriangulation(), etc., is only meant for the interaction between two volumes/surfaces in direct contact at the boundary and not separated by some finite distance.
James
Hi @James_Gorman,
This is correct. The InterfaceTriangulation is to couple two volumes via a surface.
The implementation of the āthink interfaceā approach will depend on the actual formulation/way you consider. If your formulation can be written in terms of integrals of functions over surface/volumes, it will provably possible to implement it via the high-level interface of Gridap. If you express it in terms of nodes/cells, local vectors, etc, then you will need to do low level work (but also possible since all the code is written in Julia).
czm_fracture_mwe.jl (4.2 KB)
lap-shear6ātriārename_msh.jl (1.4 MB)
I thought I implemented the boundary conditions correctly as integrals in the āminimum working exampleā (czm_fracture_mwe.jl, while the necessary mesh needs to be renamed to lap-shear6ātri.msh), but unfortunately I may have ran into an error. When I run the attached script with the linear form (the presumed cause of the issue):
l((v_t,v_b)) = -ā«( (trac(disp_correction(v_t .- v_b), k_slope, interface_thickness)ā
v_t) )dĪ_cz_top + ā«( (trac(disp_correction(v_t .- v_b), k_slope, interface_thickness)ā
v_b) )dĪ_cz_bot
I get an error
LoadError: MethodError: no method matching getindex(::Gridap.CellData.OperationCellField{ReferenceDomain}, ::Int64)
Any suggestions on what I need to correct to fix the error? The script does do a rotation on the boundary conditions (perhaps part of what is causing the issue), but this is necessary because the interfaces will likely move prior to crack propagation.
@fverdugo, for the moment I decided to go with the low-level implementation of the fracture problem using the local nodes/vectors. I would like to make the code general to both 2D and 3D problems, but I do not seem to understand get_cell_coordinates() syntax because I only access 2D elements in a 3D mesh (quad elements, not cubic elements). Is there an option to explicitly cycle through volume elements for 3D problems?
Separately, do you know of a quicker option to sort through all the elements? In the 2D problem I am currently running through all the elements and flagging all the elements along the x-axis, but this approach wonāt work for a general fracture problem. I tried selecting a portion of the domain (like the āsolidā section in FSI Tutorial 11), but this still produced all the same elements from the imported mesh.