Phase field modelling of microscopic failure in composite laminates

Abstract

The previous study of the authors on a phase field method for modelling progressive failure in multi-phase materials is further extended in this study to capture microscopic failure processes in composite laminates. The auxiliary phase field model for the regularization of material interface is re-formulated. In the computational model, a small part of the model is built from microscopic length scale to explicitly consider the microstructure of the material while the other area is built from macroscopic length scale to improve the efficiency. Cohesive elements are employed to link the two mesh parts with different length scales and also to capture the intralaminar failure. The existing ABAQUS subroutine “UEL” has been rewritten for the present study due to the introduction of cohesive elements and a BFGS solver is used to improve the iteration stability. The thickness effect of a thin ply made laminate is considered and the crack suppression effect is investigated.

Keywords

Phase field method delamination matrix crack thickness effect cohesive model

Introduction

Reinforced composite laminates are now widely used in industrial applications. A deep insight of the failure mechanism at early stages in composite laminates and the left carrying capability is desired. Essentially, composite laminates are multiphase materials. Even for a single lamina, it is composed of fiber, matrix, material interface, voids, defects, etc. From the point of view of a larger length scale, a laminate is made of several laminas by stacking them together with a designed sequence. Although various well designed experiments are available, most of them are conducted and observed from macroscopic length scale, a precise experimental observation of a complete failure process from microscopic length scale is however rather demanding.^1–3

In addition to experimental tests, virtual tests that rely on faithful numerical methods may serve as a good alternative. Virtual tests do not require the preparation of specimens and test set up, the job is ran with a computer and thus less costly. The kernel of a virtual test is the numerical method it uses, because it directly decides the modelling accuracy and efficiency. On one hand, we need a full understanding of the global response of laminates from macroscopic length scale. And the relative modelling scheme is more or less well developed. With the increasing demand from engineering, the local failure behavior of the material from microscopic length scale must be taken into consideration. Homogenization schemes used in most of the modelling approaches have prevented it from digging into the mechanical response of the microstructure and thus may bring significant error. So far, a reliable numerical method for the modelling of tensile failure in composite laminate from different length scales is still desired. We thus develop a phase field modelling scheme for this purpose in the present study.

Phase field model for fracture was developed on the basis of the variational principle proposed by Francfort and Marigo.⁴ In the variational principle, the total system potential energy functional has further included the energy dissipated for generating new crack surfaces. The dissipated energy is calculated by the critical energy release rate multiplied by the new crack surface area and thus makes the phase field model in consistence with the Griffith’s fracture theory.⁵ The phase field value can be solved with given boundary and loading conditions, and crack evolution including, crack nucleation, propagation, branching, etc, can be described by the evolution of the phase field value.⁶ In this way, complex failure processes can modelled without any crack tracking strategy. Due to the advantages of the phase field modelling, it attracts great attention. So far, phase field model for fracture has been widely used for ductile fracture problems,^7–11 brittle fracture problems,^12–21 anisotropic fracture problems,^22–27 cohesive fracture problems,^28–34 and dynamic fracture problems.^35–40 Phase field models are also developed to model material failure process in composite materials.^41–45 Phase field is normally computationally expensive for solving the nonlinear equations, a multi-scale modelling scheme in Patil et al.¹⁸ for higher efficiency was thus proposed. Only the area around the crack tip needs to use refined mesh in which the phase field equations are solved, whereas other areas are meshed using coarse meshes. The multi-scale phase field modelling scheme was then further enhanced and extended for more complicated problems i.e., composite materials, heterogeneous materials.^46–48 Recently, a few new adaptive mesh refinement strategies were developed for the localizing gradient damage model, and the refinement criterion and data transfer schemes are also studied for higher efficiency.⁴⁹

In the previous study, the authors developed a few phase field models for modelling failure process in composites.^50–53 Among the contributions, we introduced an auxiliary interface phase field model in Zhang et al.⁵³ to regularize the interface between inclusion and matrix. Such that the failure process in composite at microscopic length scale involving inclusion debonding and matrix cracking can be modelled effectively. However, the phase field model was developed mostly for modelling damage/failure in a small portion of material in lamina from microscopic length scale, it is inapplicable when failure processes in laminate are investigated. Specially, in a computational model involving multiple length scales, the existing microscopic phase field model needs further enhancement to account for other failure mechanisms such as delamination and the interaction between delamination and matrix cracking.

In this study, the existing modelling scheme is further enhanced to predict microscopic failures in composite laminates. A computational model with multi-grid meshes is introduced in detail, namely, a refined mesh region with randomly distributed fibers from microscopic length scale, and a coarse mesh region with homogenized material properties from macroscopic length scale are used. Cohesive elements are used to link the two parts of the meshes, and also to capture interlaminar failure between layers. The phase field model for the regularization of the material interface between inclusion and matrix is simplified for the fiber reinforced composite and an analytical solution is obtained.

Phase field model of cohesive model based cracks

The basic concept of the phase field model for fracture is depicted in Figure 1, and the fundamental equations are summarized as follows

\{\begin{cases} ω (d) \nabla σ + f = 0 & i n Ω \\ ω (d) σ \cdot n = t & i n \partial Ω_{t} \\ ω' (d) H + G_{s} \frac{1}{c_{0}} [\frac{1}{l_{0}} α' (d) - 2 l_{0} Δ d] = 0 & i n Ω \end{cases}

(1)

Figure 1.

Explicit and smeared crack models. (a) Explicit model, (b) Smeared model.

In the above equations, $σ$ and $ε$ are the stress and strain vectors, and $f$ and $t$ are vectors of body force and surface traction densities, respectively. $n$ is the outer normal vector. The scalar variable $d$ is the value of the fracture phase field which represents material damage state and $Δ d$ is the Laplacian of the phase field. The degradation function $ω (d)$ reflects material degradation due to damage. $α (d)$ is the geometric function and $c_{0} = 4 \int_{0}^{1} \sqrt{α (d)} d d$ . $G_{s}$ is the material’s critical energy release rate (CERR). The energy component

H = \max_{τ \in [0, t]} ψ^{+} (ε (x, τ))

(2)

is a historical value adopted to prevent self-healing.

ψ^{+}

is the split potential energy density whereas

H

is the largest value of

ψ^{+}

over the loading history. Here, we use a spectrum energy decomposition scheme⁶ to split potential energy into tensile and compressive parts, only the tensile part of the energy

ψ^{+}

is used to drive the evolution of the phase field, so as to prevent material from fracturing under compressive state. In the phase field model for fracture, a crack is smeared over its nearby region. And the internal length scale

l_{0}

in equation (1) reflects the width of the smeared zone. According to Wu’s study,³¹ it is known that a special combination of the geometric function

α (d)

and the degradation function

ω (d)

can be used to embed cohesive crack models in the phase field model. The general forms of the functions are given by

\begin{array}{l} ω (d) = \frac{{(1 - d)}^{p}}{{(1 - d)}^{p} + Q (d)}, \\ Q (d) = a_{1} d + a_{1} a_{2} d^{2} + a_{1} a_{2} a_{3} d^{3} + \dots \end{array}

(3)

α (d) = 2 d - d^{2}

(4)

where

p > 0

and the coefficients

a_{i}

are calculated from material properties and the cohesive constitutive relations. For the widely used cohesive zone model with linearly softening relationship, the coefficients are specified by

p = 2, a_{2} = - \frac{1}{2}, a_{3} = 0

(5)

It is noteworthy that for some other common cohesive zone models, the coefficients can also be derived and are available in Wu.³¹

An auxiliary phase field for fiber-matrix interface

In the modeling of transverse cracking in composite laminates from microscopic length scale, matrix cracking and fiber-matrix interface debonding are normally considered key failure mechanisms. When modelling such failure processes using the phase field model, CERR is used as an important material property. Although in reality material is continuous over matrix, fiber and the interface, and the material property CERR is thus also continuous, the CERRs of matrix and the interface are normally considered two different material properties. And experimental test procedures were designed to measure the two CERRs. Thus, there is a sudden change of CERR (see Figure 2(b)) without using any regularization scheme. Obviously, the material sudden change (or material heterogeneity) has made the modelling with phase field method challenging. Especially if we want to model failure processes in both the matrix and the interface.

Figure 2.

An auxiliary phase field and the distribution of material property (CERR), (a) is the one-dimensional distribution of the auxiliary phase field, (b) is the material property distribution without using the auxiliary phase field regularization, and (c) is the material property distribution using the auxiliary phase field regularization.

In the previous study,⁵³ the authors have proposed an interface regularization scheme to eliminate the material heterogeneity by introducing an auxiliary phase field termed as the interface phase field (IPF). The previous regularization theory was proposed mostly for concrete or materials in which the shape of the inclusion could be arbitrary. And an image processing scheme by identifying the value of the picture pixels was introduced to distinguish different material phases.

In the present investigation of fiber reinforced composite, the cross section of a fiber is simply a circular. Hence, the auxiliary IPF model can be reformulated, and the identification of the material phases can be conducted through an analytical but simple way. The IPF function $η (x)$ in a 1 D domain can be re-written as

η (x) = {\begin{cases} 1 at | x | \leq D_{i} / 2 \\ e^{(D_{i} / 2 - | x - O_{i} |) / l_{i}} otherwise \end{cases}

(6)

where

D_{i}

and

O_{i}

are the diameter and center coordinates of

i th

fiber, respectively. The distribution of the IPF is shown in Figure 2(a). Now, the zero thickness interface is transformed into a narrow band with a specific width which is affected by the internal length scale

l_{i}

. In the narrow band, the equivalent material property should be a combination of those of matrix and the interface. In the previous study, it has been proven that the following distribution gives reasonable numerical predictions, i.e.,

Θ_{e} (η) = \{\begin{matrix} Θ_{m} {(1 - η)}^{2} + Θ_{i} [2 η - η^{2}] & for & 0 \leq η < 1 \\ Θ_{f} & for & η = 1 \end{matrix}

(7)

where

Θ_{e}

Θ_{m}

Θ_{i}

Θ_{f}

are material properties in different material phases.

Θ

can be strength or CERR, and the subscripts “

e

m

i

and

f

“denote “equivalent”, “matrix”, “the interface”, and “fiber”, respectively. For the case where an overlapping exist between two and more adjacent IPFs, the effective IPF takes the one with the largest value, i.e.,

η (x) = \max (η_{1} (x), η_{2} (x), η_{3} (x), \dots)

(8)

The above one-dimensional (1 D) theory can easily be extended to a two-dimensional (2 D) case, and it is not discussed here for the sake of simplicity.

Generation of the local microscopic computational model

In the previous study,⁵³ the distribution of $η$ is obtained through solving the finite element equation. For this purpose, one has to set the value of $η$ to be 1 to satisfy the Dirichlet boundary condition. The image processing scheme was used to find out all the target nodes of fiber boundaries. The image processing scheme, however, could be complex and time consuming for the case with multiple fibers. As mentioned above, the complex image processing is unnecessary in this study. We first generate a set of randomly distributed points as fiber centers,⁵⁴ and set the value of $η$ for the surrounding points of each fiber center to one. The value of $η$ could be updated according to equation (8) if the value of a new IPF is larger. Although it is not limited to, it is assumed that all the fibers have the same size. It is noteworthy that for each IPF we do not set its value over the whole domain, this is because the value of IPF decays rapidly and is negligible in the area far away. In practice, we just set a circular area with a diameter of $2 D$ which is estimated to be adequate. Depending on the value of the internal length scale $l_{i}$ , the size of the area should be adjusted if necessary. The value of $η$ is continuous over the domain and hence the equivalent material properties are also continuous, the heterogeneous material is thus homogenized with graded properties. On the basis of various experimental observations, we assume that fiber does not sustain any damage or failure throughout the modelling. In practice, we use an artificial fracture toughness with large value for fibers. Although the material property sudden change still exists between the interface and fiber, it does not affect the phase field modelling.

Generation of the global computational model with multi-grid meshes

With the scheme for the generation of microscopic models introduced above, we further generate a multi-grid computational model. Taking the test specimen considered in ‘Validation’ section 8.2 as an example, the pre-processing scheme for other stacking sequences is similar, and details of the pre-processing are depicted in Figure 3. The $90^{o}$ plies are modelled from microscopic length scale while the $0^{o}$ plies are modelled from a macroscopic length scale. Hence, a free mesh strategy with isoparametric quadrilateral elements (marked as “bulk elements, macro-scale” in Figure 3) is used for the $0^{o}$ plies. In order to capture the whole process of transverse cracking in the $90^{o}$ plies, including crack nucleation (fiber-matrix and matrix), extension and merging, the material microstructure in these plies are considered in the computational model. Despite the shapes of the cross section of fiber, the microscopic model is built based on a structured mesh with rectangular (square) four-node regular elements (marked as “bulk-elements, micro-scale” in Figure 3). We thus firstly generate a few randomly distributed center points of the fibers and the fiber distribution is controlled such that the periodic condition can be satisfied. Then we set the values of the IPF at the nodes according to equation (8) so as to form the equivalent region of the interface. See the details of the regular mesh used to build the microstructure at lower right part of Figure 3. In order to improve the resolution of the fiber, we use dense meshes for the $90^{o}$ ply. Contours of the IPF in the $90^{o}$ plies with three different thicknesses are shown in Figure 4. In the figure, the interior areas of the fiber cross section are not shown, and the patterns in Figure 4 actually reflect position and distribution of the fibers.

Figure 3.

The detailed finite element model of the tested specimen.

Figure 4.

The IPF distributions in the $90^{o}$ plies (the red circles are the equivalent fibers interfaces).

In addition to matrix cracking and fiber debonding, delamination between adjacent plies should also be considered when investigating failure processes in composite laminates. According to the existing studies, cohesive element is considered one of the most proper and popular approaches for the modelling of delamination in composite laminates. Cohesive elements are placed along the interface between adjacent plies which essentially is the potential crack path of delamination. Moreover, the numerical predictions by using cohesive elements are very accurate in comparison with experimental observations as reported repeatedly in open literatures. So, we choose to use cohesive elements and insert them into the ply interfaces between the $0^{o}$ and $90^{o}$ plies in the pre-processing. Due to the limit of computational costs, only a part of the specimen can be modelled from microscopic length scale. Since we use different modelling approaches for matrix cracking and delamination, the interaction between the two methods should be addressed. In fact, the interaction between cohesive element and another modelling approach for bulk cracking⁵² has been discussed in other studies. In Lu et al.,⁵⁵ the authors reveal that a cohesive element must be partitioned artificially when a bulk crack reaches it, otherwise the expected initiation of delamination from that cohesive element may be delayed or even missed. Similar conclusions can also be found in literature.^56,57 In light of the conclusions of the existing studies, we suggest removing the cohesive element when the phase field crack reaches it. In this way, it provides an initiation (an entrance) for the delamination. The whole failure process is illustrated in Figure 5, one finds that the removal of the cohesive element is the key step to treat the interaction between a smeared crack model and the inserted cohesive elements. In order to realize the removal of a cohesive element, the material damage state of the neighboring bulk elements should be taken into account. For a cohesive element, there are two adjacent bulk elements, i.e., the upper one and the lower one, once the material in any of the two bulk elements are considered nearly failed, i.e.,

\max (d_{upper}, d_{lower}) > 0.99

(9)

where the subscripts “upper” and “lower” represent the upper element and the lower element, respectively, then we remove the cohesive element. From a mathematical perspective of view, the present treatment of element removal is essentially identical with the one in Lu et al.⁵⁵ if the present smeared crack tends to be an explicit one by using an infinitesimal internal length scale in an extremely dense mesh.

Figure 5.

Interaction between phase field model and cohesive element (“CE” represents cohesive element).

The fiber-matrix interface and the ply interface

In this study, two kinds of material interfaces are involved, i.e., the interface between adjacent plies and the fiber-matrix interface. We use different approaches to capture the interfacial failure processes, e.g., phase field model for the fiber-matrix interface debonding and cohesive element for delamination along the ply interface. Despite the fact that cohesive element is also proven proper and has been widely used for the modelling of fiber debonding, it is still recommended, like in the present modelling strategy, to use phase field model for the modelling of debonding. It is mostly because we consider a more complicated failure mechanism involving fiber debonding, matrix cracking and the interaction between them. Cohesive element might be satisfactory when modelling a single failure mechanism in material interfaces, fiber debonding or delamination, it brings difficulties when combined with other numerical methods for matrix cracking (phase field model in the present study). But in the proposed strategy, the fiber interface is smeared over the surrounding area, and both matrix cracking and fiber debonding are treated with the fracture phase field model. The interaction between fiber debonding and matrix cracking can also be captured automatically with the phase field model, and no more special treatment is required. It is noteworthy that the cohesive elements in between the plies are necessary to link the two mesh parts, the micro-scale dense mesh and the macro-scale coarse mesh. For the case we consider purely a micro-scale model, it is then unnecessary to use cohesive elements.

The present modelling strategy also brings advantages for the generation of computational models. It is known that one of the challenges in modelling heterogeneous material is the generation of computational models because, for example, the arbitrarily distributed fibers and the material property variation among different material phases. Normally, it requires extensive efforts being devoted to analytically build the model topology and identify the material phases. We find that the proposed strategy uses structured mesh with rectangular elements even for the fiber debonding, and the pre-processing is simplified and thus more convenient.

Finite element formulation and the BFGS solver

The formed computational model is solved numerically with finite element method (FEM). The fully coupled nonlinear equation set of the fracture phase field model can be specified by

\begin{array}{l} R^{u} = F^{ext} - \int_{Ω} ω (d) {B^{u}}^{T} σ d Ω = 0 \\ R^{d} = \int_{Ω} - ω^{'} (d) H {N^{d}}^{T} \\ - \frac{G_{s}}{c_{0}} [\frac{1}{l_{0}} α^{'} (d) {N^{d}}^{T} + 2 l_{0} {[B^{d}]}^{T} \nabla d] d Ω = 0 \end{array}

(10)

where

R = {[R^{u}, R^{d}]}^{T}

is the vector of the residual vectors of displacement and phase field fields.

B^{u}

and

B^{d}

are formed from the derivatives of the shape functions. Although it doesn’t have to, we assume the shape functions

N^{u}

and

N^{d}

of displacement and phase field elements are identical. In the Newton-Raphson iteration algorithm, the residual equation is linearized which leads to

R + \frac{\partial R}{\partial g} δ g = R - K δ g = 0, K = - \frac{\partial R}{\partial g} = [\begin{matrix} K_{u u} & K_{u d} \\ K_{d u} & K_{d d} \end{matrix}]

(11)

where

g = {[u, d]}^{T}

is the vector of nodal degrees of freedom and

K

is the tangential stiffness matrix. It is known that the Newton-Raphson monolithic algorithm is efficient but unstable for the fully coupled phase field problem of which the functional is non-convex with respect to the variation of

g

. Unless one of the two fields is fixed to make the functional convex, the Newton-Raphson algorithm is not recommended in practice due to the convergence issue. In this study, we use a BFGS (Broyden–Fletcher–Goldfarb–Shanno) quasi Newton algorithm in which the tangential stiffness matrix is modified such that

\tilde{K} (δ g) = δ R, δ R = R^{k} - R^{k - 1}

(12)

For this purpose, the stiffness matrix is specified by

\tilde{K} = {\tilde{K}}^{k - 1} - \frac{(K^{k - 1} δ g) {(K^{k - 1} δ g)}^{T}}{δ g^{T} K^{k - 1} δ g} + \frac{δ R δ R^{T}}{δ g^{T} δ R}

(13)

For the systems of symmetric stiffness matrix, the inverse of the modified stiffness matrix is explicitly available

{\tilde{K}}^{- 1} = (I - \frac{δ g δ R^{T}}{δ g^{T} δ R}) {[K^{k - 1}]}^{- 1} {(I - \frac{δ g δ R^{T}}{δ g^{T} δ R})}^{T} + \frac{δ g δ g^{T}}{δ g^{T} δ R}

(14)

In fact, the tangential stiffness matrix $K$ is not necessarily symmetric and positive-definite, it is recommended to use an uncoupled stiffness matrix as the initial tangential stiffness matrix

{\tilde{K}}^{0} = [\begin{matrix} K_{u u} & 0 \\ 0 & K_{d d} \end{matrix}]

(15)

The above stiffness matrix can be guaranteed symmetric and positive-definite in practice. In fact, the modified stiffness matrix ${\tilde{K}}^{- 1}$ should not be updated for each iteration step, it is recommended to update it after 8 iteration steps to improve the solving efficiency.

The present theory is implemented into the commercial finite element package ABAQUS through its users’ subroutine (UEL). Both the bulk elements (micro-scale, for the $90^{o}$ plies) and the cohesive elements are implemented using UEL, and for the other bulk elements (macro-scale, for the $0^{o}$ plies) we use the embedded plane strain elements “CPE4” and “CPE3” directly for convenience. The BFGS quasi Newton algorithm is also available in the package, and it can be activated by including the command line in the input file “*SOLUTION TECHNIQUE, TYPE=QUASI-NEWTON, REFORM KERNEL = 8”. More details about the implementation of the phase field method in ABAQUS and the programming examples using FORTRAN can be found in our previous study.⁵⁰

Verification and validation

Verification

It is considered a square shaped FRP composite sample with 24 randomly placed fibers (volume fraction 42%) as depicted in Figure 6. The width of the sample is $33.5 μm$ and the diameter of the fibers is $5 μ m$ . Material properties of the matrix are Young’s modulus $E_{m} = 3.9 GPa$ , Poisson’s ratio $υ_{m} = 0.37$ , critical energy release rate $G_{m} = 0.5 N/mm$ and material strength $σ_{\max}^{m} = 50 MPa$ . Fibers are modeled as linear elastic solid with Young’s modulus $E_{f} = 86.9 GPa$ and Poisson’s ratio $υ_{f} = 0.23$ . And the fracture properties of the fiber-matrix interface are critical energy release rate $G_{i} = 0.25 N/mm$ and strength $σ_{\max}^{i} = 25 MPa$ . All the material properties are obtained from Alfaro et al.⁵⁸ The boundary conditions as depicted in Figure 6 are considered in the simulation. A regular grid mesh of $670 \times 670$ elements with element size $h_{e} = 0.05 μm$ is used.

Figure 6.

A composite lamina subjected to tensile loading and the distribution of the fibers (unit: micron).

The microstructure and the contours of $η (x)$ of the investigated material is shown in Figure 6. The narrow bands in Figure 6 around the fibers are the smooth transition areas with equivalent interface properties. The predicted stress-strain curve converted from the load-displacement curve is depicted in Figure 7, in which four points A, B, C and D corresponding to the strain $ε = 0.0179$ , $0.0537$ , $0.537$ and $0.896$ are marked. The crack patterns corresponding to the four points are also depicted in Figure 7. For point “A”, it is the stage of crack initiation with slight damage, and the material of the structure is nearly intact. With the increasing of the external loading, the plate starts to sustain damage in the form of micro cracks initiated from the transition area close to the fiber borders with high stress concentration. The micro cracks in the transition area continue to propagate and eventually migrate into the matrix to merge into each other as shown in the crack pattern of point “B”. Due to material damage, the external loading starts to drop in this stage. With the further increment of external loading, the micro cracks keep evolving and a longer crack is eventually formed at point “C”. It can be seen that the reaction force has decreased to about 20% from the peak force at point “C” and it continues to decrease until the sample is nearly broken into two parts at point “D”. The strength of the material is about 45 MPa which is higher than that of the interface 25 MPa and lower than that of the matrix 50 MPa. With the consideration of the fiber-matrix interface which has a lower strength in the modelling, the present results suggest that the introduction of fiber has reduced the overall strength but has made the material more ductile. In Paris et al.,⁵⁹ the process of transverse failure in composite material was studied experimentally and the results are also depicted in the sub-figures of Figure 7. The failure patterns obtained from Paris et al.⁵⁹ are marked by “(a)”, “(b)”, “(c)” and “(d)”, respectively. In the present simulation results, it can be seen that all the experimentally observed failure stages have been captured with a correct sequence of occurrence.

Figure 7.

Modelling results of the specimen shown in Figure 6, stress-strain curve and failure mode.

Validation

An experiment conducted by Saito et al.,¹ is modelled to further validate the proposed method. In the experiment, a cross-ply composite laminate subjected to a tensile load was tested, and fiber-matrix debonding in the $90^{o}$ plies was observed and recorded. The cross-ply specimen has a stacking sequence of $[0^{o}_{2} / 90^{o}_{n} / 0^{o}_{2}]$ ( $n = 1, 2, 4$ ). The $0^{o}$ plies are made of ordinary prepreg with normal thickness and the $90^{o}$ plies are made of special thin prepreg. The thickness of a $0^{o}$ ply is $240 μ m$ while that of a single $90^{o}$ ply is $40 μ m$ . The experiment set up is illustrated in Figure 8, the specimen is subjected to a tensile loading acting on the two ends of it. The total thickness $t$ of the specimen depends on how many $90^{o}$ plies are used. For $n = 1$ , $2$ and $4$ , the corresponding total thickness of the specimen is $t = 1.00 mm$ , $1.04 mm$ and $1.12 mm$ , respectively. The finite element model of the test with detailed mesh strategy has been introduced in ‘Generation of the global computational model with multi-grid meshes’ section. For $n = 1$ , $2$ and $4$ , the total number of elements is 81543, 134677 and 241235, and the total number of nodes is 81515, 134799 and 241595, respectively. It is noteworthy that the $0^{o}$ plies were considered to have a homogeneous isotropic material property in Saito et al.,² and this hypothesis is also adopted in the study. In fact, $0^{o}$ plies in the conducted experiment are mainly used to prevent rapid transverse cracking in the $90^{o}$ plies, and the material details of the $0^{o}$ plies may be less important to the modelling results. Material properties of the laminate used in the modelling are listed in Table 1. Most of the material properties are available in literature^1,2 and the others are provided by the author of literature.^1,2

Figure 8.

The experimental set up of Saito’s test. (“N” is the number of the $90^{o}$ plies).

Table 1.

Material properties used in the modelling.^1,2

Property	Value
$90^{o}$ Ply fiber diameter	0.0025 mm
$90^{o}$ Ply fiber volume fraction	55%
$0^{o}$ Ply equivalent Young’s modulus	130 GPa
$0^{o}$ Ply equivalent Poisson’s ratio	0.33
$90^{o}$ Ply resin Young’s modulus	8.31 GPa
$90^{o}$ Ply resin Poisson’s ratio	0.45
$90^{o}$ Ply resin critical energy release rate (CERR)	0.04 N/mm
$90^{o}$ Ply resin cohesive strength	75 MPa
$90^{o}$ Ply fiber-matrix interface CERR	0.006 N/mm
$90^{o}$ Ply fiber-matrix interface cohesive strength	50 MPa
Ply interface CERR (normal)	0.04 N/mm
Ply interface cohesive strength (normal)	75 MPa
Ply interface CERR (shear)	0.12 N/mm
Ply interface cohesive strength (shear)	150 MPa

Due to the limit of computational costs, only a part of the specimen with a length of 0.12 mm is modelled in this study, see the right sub-figure of Figure 8. A parametric study with 0.15 mm, 0.18 mm, and 0.20 mm has been conducted, and the results suggest that 0.12 mm is adequate because modellings with other lengths give converged results.

The crack patterns at different loading stages are shown in Figure 9, in which the dark curves represent intralaminar failures (matrix-fiber debonding are matrix cracking) and the pink line segments are interlaminar failures (delamination). The intralaminar failure patterns are generated automatically using a “Tecplot” software with the calculation results. The intralaminar failure patterns are drawn artificially using a “Photoshop” software onto the former failure pattern based on the calculation results. The whole failure processing is divided into five stages as shown in Figure 9, i.e., the tensile strain reaches 0.5%, 0.8%, 1.0%, 1.5%, 2.0%. We term the five stages “stage I, II, III, VI and V” for simplicity. The three specimens mainly undergo an elastic stage from the start to “stage I” since not a clear damage pattern is observed. Intralaminar failure is found in “stage II” for all the specimens. Fiber-matrix debonding takes place first and kinks into the matrix, the cracks then merge with each other. Once reaching the ply interface, the previously formed intralaminar crack could then trigger a delamination. The first delamination triggered by the matrix crack is observed in the specimen with 4 layers of $90^{o}$ plies at “stage III”. The delamination is further extended with the increasing of the loading. Delamination triggered by matrix crack is also observed in the other two specimens at “stage IV”. However, the delamination crack in the specimen with 1 layer of $90^{o}$ ply is arrested by matrix cracks, and it is unlikely to propagate like in the other two thicker specimens.

Figure 9.

Numerical prediction of the failure process (black curves represent intralaminar failure while pink lines represent interlaminar failure).

The variation of crack density (which is defined by the total crack length in the thickness direction divided by the thickness of the $90^{o}$ plies) with respect to tensile strain is depicted in Figure 10. The present results agree with the experimental results.^1,2 The $90^{o}$ plies are fully penetrated by matrix crack if crack density reaches 1.0. The critical tensile strain of penetration for the thick-ply composite ( $n = 4$ ) is 1.0% and the thin-ply composite ( $n = 1$ ) is 1.4%. It is found that the thick-ply composite is easier to be penetrated then the thin-ply one. The crack propagation speed is also slowed down in the thin-ply composite ( $n = 1$ ), unlike in the other two cases, we don’t find a sudden increment of crack density. As reported in Saito et al.,² the $0^{o}$ plies suppress the crack propagation, and this crack suppression effect is more significant in thinner plies.

Figure 10.

Crack density vs. tensile strain, present prediction and the experimental results.²

Numerical predictions of crack pattern when crack density reaches 1.0 are depicted in Figure 11, in which the experimental observations are provided. In comparison with thick-ply composites, matrix crack is also shorter in thin-ply composites. It is known that stress intensity factors at the crack tips of longer cracks are normally larger and thus easier to trigger another crack (delamination in this case). And this explains why the delamination in Figure 9 for $n = 4$ at 2.0% strain looks more serious than the other cases. From the results, it is found that delamination is also suppressed in the thin-ply composite, and this could be another crack suppression mechanism of thin-ply composites.

Figure 11.

Crack paths obtained from the experimental test² and the present modelling.

Conclusion

As an extension of the previous study, the auxiliary interface phase field model is in fact a special case of the general model for the inclusion with arbitrary shape. The auxiliary interface phase field model is thus analytically derived and explicitly specified. Composite lamina with complex microstructure shapes can be modelled on a regular mesh with rectangular elements with the proposed model. The introduction of cohesive elements in the computational model has made it possible to capture the interlaminar failure. The proposed phase field modelling scheme is effective in capturing a complicated failure process with multiple failure mechanisms in composite laminates. In addition to the findings reported in Saito et al.,² the present numerical results suggest that the suppression of delamination could be another crack suppression mechanism of thin-ply composites. The crack suppression capability of a laminate is found to have a direct correlation with its capability in limiting the length of matrix cracks.

Footnotes

Acknowledgement

Professor Hiroshi Saito is sincerely appreciated for providing material properties and experimental data.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (No. 11872143) and the National Key Research and Development Program of China [No.2016YFB0200702].

ORCID iDs

Xiaofei Hu

Weian Yao

Appendix

References

Saito

Takeuchi

Kimpara

Experimental evaluation of the damage growth restraining in 90 layer of thin-ply CFRP cross-ply laminates. Adv Compo Mater 2012; 21: 57–66.

Saito

Takeuchi

Kimpara

A study of crack suppression mechanism of thin-ply carbon-fiber-reinforced polymer laminate with mesoscopic numerical simulation. J Compos Mater 2014; 48: 2085–2096.

Mortell

Tanner

McCarthy

CT.

In-situ SEM study of transverse cracking and delamination in laminated composite materials. Compos Sci Technol 2014; 105: 118–126.

Francfort

Marigo

JJ.

Revisiting brittle fracture as an energy minimization problem. J Mech Phys Solids 1998; 46: 1319–1342.

Griffith

AA.

The phenomena of rupture and flow in solids. Philos Trans R Soc London 1921; 221: 163–198.

Miehe

Welschinger

Hofacker

Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations. Int J Numer Meth Engng 2010; 83: 1273–1311.

Ambati

Gerasimov

De Lorenzis

Phase-field modeling of ductile fracture. Comput Mech 2015; 55: 1017–1040.

Alessi

Marigo

J-J

Vidoli

Gradient damage models coupled with plasticity: variational formulation and main properties. Mech Mater 2015; 80: 351–367.

Duda

Ciarbonetti

Sánchez

, et al. A phase-field/gradient damage model for brittle fracture in elastic–plastic solids. Int J Plast 2015; 65: 269–296.

10.

Borden

Hughes

Landis

, et al. A phase-field formulation for fracture in ductile materials: finite deformation balance law derivation, plastic degradation, and stress triaxiality effects. Comput Meth Appl Mech Eng 2016; 312: 130–166.

11.

Miehe

Hofacker

Schänzel

L-M

, et al. Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic–plastic solids. Comput Meth Appl Mech Eng 2015; 294: 486–522.

12.

Miehe

Hofacker

Welschinger

A phase field model for rate-independent crack propagation: robust algorithmic implementation based on operator splits. Comput Meth Appl Mech Eng 2010; 199: 2765–2778.

13.

Bourdin

Francfort

Marigo

JJ.

The variational approach to fracture. Berlin: Springer Science+Business Media B.V, 2008.

14.

Tanne

Bourdin

, et al. Crack nucleation in variational phase-field models of brittle fracture. J Mech Phys Solids 2018; 110: 80–99.

15.

Ambati

Gerasimov

Lorenzis

LD.

A review on phase-field models of brittle fracture and a new fast hybrid formulation. Comput Mech 2015; 55: 383–405.

16.

Msekh

Sargado

Jamshidian

Areias

, et al. Abaqus implementation of phase-field model for brittle fracture. Comput Mater Sci 2015; 96: 472–484.

17.

Nguyen

Yvonnet

Bornert

, et al. On the choice of parameters in the phase field method for simulating crack initiation with experimental validation. Int J Fract 2016; 197: 213–226.

18.

Patil

Mishra

Singh

IV.

An adaptive multiscale phase field method for brittle fracture. Comput Meth Appl Mech Eng 2018; 329: 254–288.

19.

Zhang

Vignes

Sloan

, et al. Numerical evaluation of the phase-field model for brittle fracture with emphasis on the length scale. Comput Mech 2017; 59: 737–752.

20.

Mesgarnejad

Bourdin

Khonsari

Validation simulations for the variational approach to fracture. Comput Meth Appl Mech Eng 2015; 290: 420–437.

21.

Huang

Nguyen

VP.

On the BFGS monolithic algorithm for the unified phase field damage theory. Comput Meth Appl Mech Eng 2020; 360: 112704.

22.

Bryant

Sun

A mixed-mode phase field fracture model in anisotropic rocks with consistent kinematics. Comput Meth Appl Mech Eng 2018; 342: 561–584.

23.

Peco

Millán

, et al. Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy. Int J Numer Meth Eng 2015; 102: 711–727.

24.

Teichtmeister

Kienle

Aldakheel

, et al. Phase field modeling of fracture in anisotropic brittle solids. Int J Non-Linear Mech 2017; 97: 1–21.

25.

Nguyen

Réthoré

Baietto

M-C.

Phase field modelling of anisotropic crack propagation. Eur J Mech – A/Solids 2017; 65: 279–288.

26.

Bleyer

Alessi

Phase-field modeling of anisotropic brittle fracture including several damage mechanisms. Comput Meth Appl Mech Eng 2018; 336: 213–236.

27.

Clayton

Knap

Phase field modeling of directional fracture in anisotropic polycrystals. Comput Mater Sci 2015; 98: 158–169.

28.

Paggi

Reinoso

Revisiting the problem of a crack impinging on an interface: a modeling framework for the interaction between the phase field approach for brittle fracture and the interface cohesive zone model. Comput Meth Appl Mech Eng 2017; 321: 145–172.

29.

Carollo

Reinoso

Paggi

A 3D finite strain model for intralayer and interlayer crack simulation coupling the phase field approach and cohesive zone model. Compos Struct 2017; 182: 636–651.

30.

Conti

Focardi

Iurlano

Phase field approximation of cohesive fracture models. Ann Linstitut Henri Poincaré 2016; 33: 1033–1067.

31.

JY.

A unified phase-field theory for the mechanics of damage and quasi-brittle failure. J Mech Phys Solids 2017; 103: 72–99.

32.

Nguyen

T-T

Yvonnet

Zhu

Q-Z

, et al. A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography. Comput Meth Appl Mech Eng 2016; 312: 567–595.

33.

Freddi

Iurlano

Numerical insight of a variational smeared approach to cohesive fracture. J Mech Phys Solids 2017; 98: 156–171.

34.

Verhoosel

De Borst

Eacute. A phase-field model for cohesive fracture. Int J Numer Meth Eng 2013; 96: 43–62.

35.

Ren

Zhuang

Anitescu

, et al. An explicit phase field method for brittle dynamic fracture. Comput Struct 2019; 217: 45–56.

36.

Borden

Verhoosel

Scott

, et al. Phase-field description of dynamic brittle fracture. Comput Meth Appl Mech Eng 2012; 217-220: 77–95.

37.

Bourdin

Larsen

Richardson

CL.

A time-discrete model for dynamic fracture based on crack regularization. Int J Fract 2011; 168: 133–143.

38.

Schlüter

Willenbücher

Kuhn

, et al. Phase field approximation of dynamic brittle fracture. Comput Mech 2014; 54: 1141–1161.

39.

Doan

Bui

Van Do

, et al. A rate-dependent hybrid phase field model for dynamic crack propagation. J Appl Phys 2017; 122: 115102.

40.

Doan

Bui

Duc

, et al. Hybrid phase field simulation of dynamic crack propagation in functionally graded glass-filled epoxy. Compos Part B: Eng 2016; 99: 266–276.

41.

Natarajan

Annabattula

RK.

Modeling crack propagation in variable stiffness composite laminates using the phase field method. Compos Struct 2018; 209: 424–433.

42.

Msekh

Silani

Jamshidian

Areias

, et al. Predictions of J integral and tensile strength of clay/epoxy nanocomposites material using phase field model. Compos Part B: Eng 2016; 93: 97–114.

43.

Espadas-Escalante

van Dijk

Isaksson

A phase-field model for strength and fracture analyses of fiber-reinforced composites. Compos Sci Technol 2019; 174: 58–67.

44.

Alessi

Freddi

Phase-field modelling of failure in hybrid laminates. Compos Struct 2017; 181: 9–25.

45.

Msekh

Cuong

, et al. Fracture properties prediction of clay/epoxy nanocomposites with interphase zones using a phase field model. Eng Fract Mech 2018; 188: 287–299.

46.

Patil

Mishra

Singh

, et al. A new multiscale phase field method to simulate failure in composites. Adv Eng Softw 2018; 126: 9–33.

47.

Patil

Mishra

Singh

IV.

A local moving extended phase field method (LMXPFM) for failure analysis of brittle materials. Comput Meth Appl Mech Eng 2018; 342: 674–709.

48.

Patil

Mishra

Singh

IV.

A multiscale framework based on phase field method and XFEM to simulate fracture in highly heterogeneous materials. Theoret Appl Fract Mech 2019; 100: 390–415.

49.

Sarkar

Singh

Mishra

BK.

Adaptive mesh refinement schemes for the localizing gradient damage method based on biquadratic-bilinear coupled-field elements. Eng Fract Mech 2020; 223: 106790.

50.

Zhang

Wang

, et al. An iteration scheme for phase field model for cohesive fracture and its implementation in abaqus. Eng Fract Mech 2018; 204: 268–287.

51.

Zhang

Bui

, et al. Phase field modeling of fracture in fiber reinforced composite laminate. Int J Mech Sci 2019; 161-162: 105008.

52.

Zhang

Feng

Bui

, et al. Modelling distinct failure mechanisms in composite materials by a combined phase field method. Compos Struct 2020; 232: 111551.

53.

Zhang

Yang

, et al. Modelling progressive failure in multi-phase materials using a phase field method. Eng Fract Mech 2019; 209: 105–124.

54.

Catalanotti

On the generation of RVE-based models of composites reinforced with long fibres or spherical particles. Compos Struct 2016; 138: 84–95.

55.

Chen

Tan

VBC

, et al. Adaptive floating node method for modelling cohesive fracture of composite materials. Eng Fract Mech 2018; 194: 240–261.

56.

Tay

TE.

Modelling delamination migration using virtual embedded cohesive elements formed through floating nodes. Compos Struct 2018; 204: 500–512.

57.

Zhi

Tay

TE.

Interrogating failure mechanisms of notched composites through a discrete crack modeling approach. Compos Sci Technol 2020; 196: 108203.

58.

Alfaro

MVC

Suiker

ASJ

de Borst

Transverse failure behavior of fiber-epoxy systems. J Compos Mater 2010; 44: 1493–1516.

59.

Paris

Correa

Mantic

Kinking of transversal interface cracks between fiber and matrix. J Appl Mech – Trans ASME 2007; 74: 703–716.