Validation of a code_aster / MFront cohesive model for the delamination of laminated composites (DCB)

1. Objective

Delamination is the dominant failure mode of laminated composites, and its predictive simulation is a major industrial concern: it makes it possible to size, for instance, bonded joints, aircraft wing stiffeners or structures under impact without multiplying the number of tests. This article shows that such simulations are now within reach using open-source tools: a cohesive law written in MFront [1] is integrated into code_aster through interface finite elements.

The case study chosen here is the benchmark of Alfano & Crisfield [2]: a Double Cantilever Beam (DCB) test [3] on a carbon-fiber laminated composite specimen. Its value lies in possessing an analytical solution of the force-displacement response during crack propagation, predicted by Linear Elastic Fracture Mechanics (LEFM) [4]. The objective is therefore to compare this reference with the code_aster + MFront simulation, and to show that the latter reproduces it faithfully.

2. The Double Cantilever Beam (DCB) problem [3]

2.1 Geometry and data

The specimen is a unidirectional laminate (along $x$) pre-cracked over an initial length $a_0$. Each arm is loaded by an imposed opening displacement $\pm u$, to which a reaction $P$ corresponds.

Quantity	Symbol	Value	Unit
Geometry
Total length	$L$	100	mm
Width	$b$	20	mm
Arm thickness	$t$	1.5	mm
Total thickness	$h=2t$	3	mm
Pre-crack length	$a_0$	30	mm
Material (ply)
Longitudinal modulus	$E_{11}$	135 300	MPa
Shear modulus	$G_{12}$	5 200	MPa
Poisson's ratio	$\nu_{12}$	0.24	-
Material (interface)
Critical energy release rate (mode I)	$G_c$	0.28	N/mm
Cohesive strength	$\sigma_c$	1.7	MPa

3. The reference analytical solution (LEFM)

Since Griffith [4], the fracture of brittle and quasi-brittle materials has been described by an energy balance: a crack propagates only if the elastic energy released by its advance is sufficient to create the new surfaces. Irwin [5] then recast this criterion in a form directly usable in engineering, by introducing the energy release rate $G$, the energy available per unit of cracked surface, and its critical threshold $G_c$ beyond which the crack propagates.

The DCB test admits a closed-form analytical solution based on this energy approach, so that $G$ takes the form:

$$G = \frac{P^2}{2b}\,\frac{\mathrm{d}C}{\mathrm{d}a}\tag{1}$$

where $C = 2u/P=\delta/P$ is the compliance (the inverse of the stiffness), directly measurable experimentally.

A crack always makes the structure more compliant ($\mathrm{d}C/\mathrm{d}a > 0$), so $G$ is always positive: energy is always available to drive cracking, the only question being whether there is enough of it to reach $G_c$.

Crack propagation. Combining the Griffith criterion $G = G_c$ (an empirically determined energy) with beam theory applied to the two arms, the reaction $P$ takes the form:

$$\boxed{\; P(u) = \frac{\lambda}{\sqrt{u}}, \qquad \lambda(G_c) = 62.1\ \mathrm{N\,mm^{1/2}} \;}\tag{2}$$

Physically, as the crack advances, the structure becomes more compliant: the required load decreases while the imposed opening increases, hence a response that drops continuously. The theory places initiation at $R_i = 68.8$ N for $u_i = 0.81$ mm, and crack breakthrough ($a=L$) at a theoretical arm displacement $u_{max} = 9.04$ mm.

The loading is applied under displacement control: $G \propto \delta^2/a^4$ then decreases sharply with $a$, so that propagation arrests on its own as long as the arms are not opened further. This control therefore makes propagation intrinsically stable and removes the need for more advanced numerical methods (arc-length or dynamic time integration).

4. The cohesive zone model

Whereas LEFM assumes a pre-existing crack and an infinite stress at the tip, the cohesive zone model, introduced by Barenblatt [6] and Dugdale [7], places along the crack path a law relating the traction to the opening of the crack faces. This law is expressed through two parameters, each tied to a different criterion: a strength $\sigma_c$ that governs initiation, and an energy $G_c$ that governs propagation and recovers the Griffith criterion. In the bilinear law used here, the traction increases linearly with the opening up to the peak $\sigma_c$, then decreases to zero; the area under this traction-separation curve is exactly $G_c$, and the model recovers LEFM as $\sigma_c \to \infty$.

$$G_c = \int_{0}^{\delta_c} \sigma(\delta)\,\mathrm{d}\delta = \tfrac{1}{2}\,\sigma_c\,\delta_c \tag{3}$$

The interface is discretized by zero-thickness elements connecting two layers of coincident but distinct nodes. The kinematic quantity is not a strain but a displacement jump across the interface:

$$\delta = u^+-u^-\tag{4}$$

and the dual variable is not a stress tensor but a traction vector:

$$\mathbf{\sigma} = t_n \tag{5}$$

With $t_n$ the normal traction component (mode I), $\sigma_c$ denotes the critical stress, $\delta_0$ its associated opening and $\delta_c$ the critical opening. The normal law comprises four branches:

$$\sigma(\delta) = \begin{cases} \alpha K\,\delta & \delta < 0 & \text{(compression, penalty)}\\[4pt] K\,\delta & 0 \le \delta \le \delta_0 & \text{(elastic)}\\[4pt] (1-d)\,K\,\delta & \delta_0 < \delta^{\max} < \delta_c & \text{(damage)}\\[4pt] 0 & \delta^{\max} \ge \delta_c & \text{(failure)} \end{cases} \tag{6}$$

We introduce $d$, the damage variable, which here represents the progressive degradation of the interface before its complete failure. This variable varies continuously in $[0,1]$: $0$ for an undamaged interface, $1$ for a fully failed interface, and it degrades the interface cohesion according to:

$$(1-d)\,K\,\delta \tag{7}$$

with

$${\ d(\delta) = \frac{1}{\eta}\left(1 - \frac{\delta_0}{\delta}\right)\ } \quad \text{with} \quad \eta \equiv 1 - \frac{\delta_0}{\delta_c} \tag{8}$$

Moreover, the irreversibility of damage requires that $d$ be a monotonically increasing function of the maximum opening reached $\delta^{\max}$, which is enforced numerically as:

$$d \leftarrow \min\!\Big(\max\big(d^{\,\text{prev}},\, d(\delta)\big),\ 1\Big) \tag{9}$$

One parameter governs the numerical balance, the cohesive zone length [6,7] $\ell_{cz} \sim E_{11} \,G_c/\sigma_c^2$, which the mesh must resolve. Since $\ell_{cz} \propto 1/\sigma_c^2$, the choice of $\sigma_c$ trades off physical fidelity (high strength, short fracture zone) against numerical robustness (low strength, more extended zone). The value used here, $\sigma_c = 1.7$ MPa, favors robustness.

The numerical aspects, namely the integration of the MFront material law and the derivation and computation of the tangent operator, are not addressed here, as they are already the subject of dedicated articles elsewhere.

5. Implementation in code_aster with MFront

The full workflow relies on free tools: CAD and meshing in Salome, material law written in MFront, computation in code_aster, post-processing in ParaVis (the module adapting ParaView within Salome).

Cohesive law: The bilinear law of the cohesive zone model is written in MFront, compiled on the fly and integrated into code_aster as the joint behavior. Writing the law separately makes it possible to modify it without recompiling code_aster.

Mesh: The interface at the fracture zone is discretized by zero-thickness elements connecting two layers of duplicated, coincident nodes (3D_JOINT modeling). The arms, for their part, are modelled with a 3D model.

Boundary conditions: This is the point that determines the success of the computation. A DCB is a compliant structure, held together only by the intact joint ahead of the crack. We therefore impose only an opening displacement along $z$ on the loaded edges, with respective values $+u$ and $-u$, while leaving the sections free to rotate. Only a few individual nodes are constrained to remove the remaining rigid-body modes.

Solution. The computation is carried out as a non-linear analysis with adaptive increments and iterations, under displacement control.

6. Results and validation

With $\sigma_c = 1.7$ MPa, the curve from the simulation is compared with the analytical reference $P=62.1/\sqrt{u}$.

The agreement on the propagation branch is excellent (0.29 %): this is the signature of the model's energy consistency, which dissipates exactly $G_c$. The more compliant initial phase and the lowered peak ($\sigma_c = 1.7$ MPa) are not implementation flaws, but the direct consequence of $\ell_{cz} \gg L$: the whole interface undergoes damage from the very start of loading, without a localised crack tip. This is the behavior documented in [2] for low strengths. A higher $\sigma_c$ would restore the peak, at the cost of more delicate convergence.

7. Conclusion

A 3D finite element model of mode-I delamination, built entirely with free tools (Salome-Meca, MFront, code_aster, ParaVis), quantitatively reproduces the benchmark of Alfano & Crisfield [2], with agreement of the order of one percent on crack propagation relative to the analytical solution. This agreement relies on three well-identified ingredients: boundary conditions that preserve the compliance of the structure, a cohesive law that fully releases the traction at failure, and a solution strategy suited to the softening damage branch.

This result shows that code_aster covers a need usually reserved for commercial codes: the predictive simulation of delamination. The same approach extends naturally to mixed mode, to multi-delamination and to real bonded structures.

References

[1] T. Helfer et al., Introducing the open-source MFront code generator, Computers & Mathematics with Applications, 70 (2015) 994–1023.

[2] G. Alfano, M. A. Crisfield, Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues, International Journal for Numerical Methods in Engineering, 50 (2001) 1701–1736.

[3] ASTM D5528, Standard Test Method for Mode I Interlaminar Fracture Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites.

[4] A. A. Griffith, The phenomena of rupture and flow in solids, Philosophical Transactions of the Royal Society A, 221 (1921) 163–198.

[5] G. R. Irwin, Analysis of stresses and strains near the end of a crack traversing a plate, Journal of Applied Mechanics, 24 (1957) 361–364.

[6] G. I. Barenblatt, The mathematical theory of equilibrium cracks in brittle fracture, Advances in Applied Mechanics, 7 (1962) 55–129.

[7] D. S. Dugdale, Yielding of steel sheets containing slits, Journal of the Mechanics and Physics of Solids, 8 (1960) 100–104.