Use of Cohesive Approaches for Modelling Critical States in Fibre Reinforced Structural Materials

1. Introduction

The question of ensuring the safety of structural elements regardless of their size and predicting their service life is increasingly being asked in connection with the development of advanced materials, devices, structures and their components not only in the field of engineering or construction. This can depend on the occurrence of defects that can arise during the production phase or during their operation and loading by external conditions. The field of science, summarized as fracture mechanics, combining continuum mechanics with materials engineering, tries to predict the formation and further development of defects in structures. It is a complex defect-stress-material relationship, not included in both classical and advanced approaches to the modelling of composite structures by [1]. To understand relationships and extend lifespan, it is not only necessary to modernize procedures or relevant standards, but also to use new ones more accurate numerical methods, which with increasing probability can predict the behaviour of materials even under complex operating conditions. Using fracture mechanics, we describe and predict the behaviour of bodies containing defects.

There are two basic approaches to derive the conditions at the moment of initiation of unstable crack propagation. The first one uses the weakest link theory, following [2], the second one considers the accumulation of damage during life cycle, as discussed by [3]. Failure of structural materials is a continuous process in which the phases of plastic deformation, nucleation and initiation of cracks are intertwined. The final stage of failure development of bodies, which is the subject of fracture mechanics investigation, is crack propagation (unstable or stable). Attention is paid to understanding this behaviour of structural materials to the deterministic approaches and to the computational modelling of such processes for the case of several representative types of composites, unlike various types of probabilistic considerations by [4,5,6,7], or discrete element techniques by [8,9], or even those relying on techniques of artificial intelligence as [10]. This article used the long-time experience of both authors in this research field, which is confronted with new models, methods and results developed in recent years.

As a result of the fracture, the structure breaks into two or more parts. We define a crack as a violation of the cohesion of bodies along a surface bounded by a curve which is either closed or ends on the surface of the body. An important problem is to define criteria that will unambiguously determine the state when existing microcracks, for example those created during manufacturing processes, decide to propagate further. The concept of the stability criterion can be extended to the case of general stress concentrators, in that case, we look for the conditions under which the crack starts to propagate from the considered concentrator. In general, it is possible to formulate stability criteria and propagation criteria on the basis of energy by [11,12,13], or using the strain energy density, crack driving force, G, J integral, or based on the stress and strain at the crack tip (stress intensity factor concept, crack opening). For linear fracture mechanics, the more frequently used stability criterion is the $K_{IC}$ criterion (the $G_{IC}$ criterion can be used equivalently). For nonlinear fracture mechanics then the $J_{IC}$ criterion by [14] was formulated and with connection to cohesive energy by [15], or the crack tip opening displacement (CTOD) condition, with a series of criteria for combined stress, suggested by [16,17,18]. Connection of the $J_{IC}$ integral with special cohesive element for dynamic crack propagation can be seen in [19].

From the mathematical point of view, the complete system of partial differential equations of evolution, both in its classical differential form and in its variational or weak integral form, relies on the conservation principles of classical thermomechanics, supplied by appropriate constitutive equations, as presented by [20]. Since such formulations work with function spaces of infinite dimension, thus most computational approaches to real engineering problems need some discretizations both in the Euclidean space, 3-dimensional in general, and on certain time variable, even for a seemingly static simplified evaluation of fracture development. The conservation principles contain both the total strain tensor $\epsilon$ and the stress tensor $\sigma$, both represented by symmetric square matrices of order 3. Since the values of components ${\sigma }_{ij}$ of $\sigma$ with $i,j\in \{1,2,3\}$ depend on the choice of a Cartesian coordinate system, for the following considerations it is useful to introduce also 3 invariants of $\sigma$ (independent of the choice of Cartesian coordinates) ${\sigma }_1$, ${\sigma }_2$ and ${\sigma }_3$; the first (linear) one is ${\sigma }_1={\sigma }_{11}+{\sigma }_{22}+{\sigma }_{33}$, for the derivation of their complete set see [21], Part 1.5. In most computational tools $ϵ$ is decomposed into 4 components, denoted as ${ϵ}^e$, ${ϵ}^p$, ${ϵ}^c$ and ${ϵ}^{\theta }$, referring to elastic, plastic, creep and thermal ones; moreover an appropriate incorporation of damage process is expected.

In particular, in the simplified small deformation theory for isotropic materials, using the standard Kronecker symbol ${\delta }_{ij}=1-\mathrm{sgn}|i-j|$, for purely elastic deformation we are allowed write

$${\epsilon }_{ij}={\epsilon }_{ij}^e=\frac{1+\nu }{E}{\sigma }_{ij}-\frac{\nu }{E}{\sigma }_1{\delta }_{ij}$$

(1)

for all $i,j\in \{1,2,3\}$; this relation is the famous empiric Hook law, containing the couple of material parameters $(E,\nu )$: E is the Young modulus, $\nu$ the Poisson ratio. To derive $\epsilon$ for $\epsilon \ne {\epsilon }^e$, unlike (1), some appropriate decomposition of $\epsilon$ must be suggested: in most computational tools $\epsilon$ is considered as a simple sum of ${\epsilon }^e$, ${\epsilon }^p$, ${\epsilon }^c$ and ${\epsilon }^{\theta }$. Nevertheless, proper nonlinear formulations, namely for the elastoplasticity theory by [21], Part 6.1, and in the more general context for the theory of structured deformation by [22], revisiting the original idea of [23], require more complicated multiplicative decompositions, with non-classical material characteristics and internal variables.

The irreversible component ${\epsilon }^p$ is activated (takes non-zero values) after certain function $\mathcal{F}({\sigma }_1,{\sigma }_2,{\sigma }_3)$ reaches some prescribed value ${\sigma }_{*}$: e. g. for metals the von Mises criterion with

$$\mathcal{F}({\sigma }_1,{\sigma }_2,{\sigma }_3)=\sqrt{\frac{1}{2}\left({({\sigma }_1-{\sigma }_2)}^2+{({\sigma }_2-{\sigma }_3)}^2+{({\sigma }_3-{\sigma }_1)}^2\right)}$$

(2)

is used frequently, with ${\sigma }_{*}$ considered as the yield strength of material in simple tension; the left-hand side of (2) can be then understood as an equivalent stress. Then ${\epsilon }^p$ can be then evaluated from the associated plastic flow rule, following [21], Part 6.4, and [24]; all additional material parameters are expected to be identified from well-prepared experiments. The superposition of ${\epsilon }^c$ and ${\epsilon }^e$ relies on the Kelvin parallel viscoelastic model typically; ${\epsilon }^c$ is evaluated as a function (in general nonlinear, especially in the well-known Norton power-law model of ccreep) of stress components or invariants, frequently from the deviatoric part ${\sigma }^{\prime }$ of $\sigma$, i. e.

$${\sigma }_{ij}^{\prime }={\sigma }_{ij}-\frac{{\sigma }_1}{3}{\delta }_{ij}$$

with $i,j\in \{1,2,3\}$ again. The most simplified evaluation of ${ϵ}^{\theta }$ can rely on the linear dependence of ${ϵ}^{\theta }$ on the prescribed temperature change, using certain thermal expansion coefficient(s); the deeper analysis results in multi-physical models including both mechanical deformation and thermal transfer, which exceeds the scope of this article.

In the field of both linear and nonlinear fracture mechanics, special curvilinear integrals are introduced as its substantial ingredients, namely for the analysis of crack tip phenomena. The above mentioned energetic contour path integral, called J-integral, comes from [14,25]. An alternative way, utilized by [26,27,28,29], derives the crack propagation from the C-integral, working with energy rate. For better understanding of such individual types of curvilinear integrals the theoretical Warp3D software system manual [30] can be recommended. However, in the case of multiaxial stresses, the one-parameter fracture mechanics, referring to the analysis of the stress around the crack front, may be insufficient; the multi-parameter fracture mechanics, following by [31,32], work with such concepts as T-stresses or Q-parameters.

In the history of FEM-based numerical simulations, some simple approaches to modelling of crack propagation should be mentioned, namely the element deletion approach by [33] and the node release method by [34], upgraded as the continuous remeshing technique of [35]. Another model of accumulation of damage is processed by [36]. As an example of using damage mechanics for ductile failure, Gurson - Tvergaard - Needleman (GTN) model can be pointed out: see [37,38]. Based on the evaluation of the performed experiments, it is recommended introducing two (or three) optional parameters $q_1,q_2,q_3$, usually $q_1=1.5$, $q_2=1$, $q_3=q_1^2$. A careful analysis of the influence of these parameters on the development of damage was carried out by [39]. Here, ${\sigma }_{YS}$ is the yield stress of the matrix material, ${\sigma }_m$ is the principal (hydrostatic) stress and f is the volume fraction of voids. The energy balance of such so-called complete model by [40], implemented as a user procedure in Abaqus software, as documented by [41], gives an equation for plastic potential describing plastic flow in porous materials

$$\frac{2S^2}{3{\sigma }_{YS}^2}+2q_1\mathfrak{F}cosh\frac{3q_2^2{\sigma }_m}{2{\sigma }_{YS}}-q_3{\mathfrak{F}}^2=1$$

(3)

where, using the deviatoric stress components,

$$S=\sqrt{{\sum }_{i,j=1}^3{\sigma }_{ij}^{\prime }{\sigma }_{ij}^{\prime }}.$$

For practical calculations, the volume fraction of cavities $\mathfrak{F}$ in (3) is usually replaced by the effective volume fraction ${\mathfrak{F}}^{*}$, introduced as

$${\mathfrak{F}}^{*}={\mathfrak{F}}_c+\frac{{\mathfrak{F}}_u^{*}-{\mathfrak{F}}_c}{{\mathfrak{F}}_F-{\mathfrak{F}}_c}(\mathfrak{F}-{\mathfrak{F}}_c);$$

here ${\mathfrak{F}}_c$ is the critical volume for coalesce of voids, ${\mathfrak{F}}_F$ is the volume of voids at ultimate damage and ${\mathfrak{F}}_u^{*}=1/q_1$, as needed in (3). This effective value ${\mathfrak{F}}^{*}$ applies to (3) when ${\mathfrak{F}}_c$ is less than $\mathfrak{F}$.

FEM is widely used in the numerical analysis of initial and boundary value problems for differential equations, however, as already indicated in the article, the FEM mesh may not always be ideal for modelling crack propagation. And this is one of the most important interests in the problems of mechanics of solid bodies. The first models were established on a weak (deformation) discontinuity that could pass through the finite element mesh, using the variational principle by [42]. Other authors and researchers considered a strong (displacement) discontinuity by modification of the principle of virtual work (which also applies to models with the traction separation law), as [43,44], including the stability and convergence of such problems, as [45], and improving the accuracy of such modelling, namely [46,47].

In the case of a strong discontinuity, the displacement consists of regular and amplified components, where the improved components are given over a jump over a discontinuous surface by [44]. There is a modification of the basic equation, well-known from the standard finite element method (FEM), for the relationship between displacements $u^{ϵ}\left(x\right)$ on particular points x of the $ϵ$-th element (3-dimensional vectors of functions in general) and displacements at selected element nodes $u^i$, utilizing some standard shape functions ${\mathcal{N}}_i\left(x\right)$, i. e.

$$u^{ϵ}\left(x\right)=\sum _{i\in {\mathcal{E}}_A}{\mathcal{N}}_i\left(x\right)u^i;$$

(4)

here ${\mathcal{E}}_A$ is the set of nodes corresponding to the standard e-th element.

If we consider a modification of FEM such as XFEM (extended finite element method), some more terms must be added to slightly modified displacements inside element end displacements at element nodes. It is the second term that points to the technique of penetration into the element and the crack movement, the third term the failure (decohesion) and the possible direction of movement according to the preferred criteria. So (4) for the case of crack growth for 2-dimensional modelling (not limited to one element) is modified into the form

$$u^{ϵ}\left(x\right)=\sum _{i\in {\mathcal{E}}_A}{\mathcal{N}}_i\left(x\right)u^i+\sum _{j\in {\mathcal{E}}_B}{\mathcal{N}}_j\left(x\right){\mathcal{H}}_j\left(x\right)+\sum _{k\in {\mathcal{E}}_C}{\mathcal{N}}_k\left(x\right)\sum _{m=1}^{\mathcal{M}}{\varphi }_k^m\left(x\right)c_k^m,$$

(5)

including certain specialized shape functions ${\mathcal{N}}_j\left(x\right)$ and ${\mathcal{N}}_k\left(x\right)$ for the intrinsic version of XFEM, unlike the extrinsic version adopting the original functions ${\mathcal{N}}_i\left(x\right)$ from (4); $\mathcal{H}\left(x\right)$ here is certain discontinous function with values 0 and 1, interpretable as the Heaviside function in the local crack coordinate system. Whereas the first additive right-hand-side term of (5) is identical with that of (4), the second one supports the evaluation of crack formation, working with some set of nodes ${\mathcal{E}}_B$, and the third one the inclusion of the local phenomena in front of the crack, working with a still other set of nodes ${\mathcal{E}}_C$. The well-known setting with $\mathcal{M}=4$, used in the two-dimensional crack propagation modelling frequently, relies on the choice of special functions ${\varphi }_k^m\left(x\right)$ by [48]. Figure 1 illustrates the location of areas A, B and C, related to a crack; the sets of nodes ${\mathcal{E}}_A$, ${\mathcal{E}}_B$ and ${\mathcal{E}}_C$ ot these areas occur in (5) and (4). The progress of this approach, adopted to the three-dimensional modelling at the expense of much more complicated functions ${\varphi }_k^m\left(x\right)$, in the last years can be seen from [49].

In general, the method of discretization in time and the extended finite element method (XFEM) can be used for the adaptive change and / or enrichment of the set of basis functions near singularities. This computational approach (involving its numerous variants with own names and designations) already has quite a rich history with the remarkable progress in recent years, as documented by [50,51,52,53,54,55], including the detailed analysis of convergence properties for various enrichment strategies by [56]. Namely an improvement of XFEM, utilizing an upgrade of (4), presented as IXFEM, is introduced by [57] to reach inexpensive computations with good geometrical accuracy and surface mesh resolution independence. Another approach works relies on the simplification of “smeared” damage by [58,59,60,61,62]: instead of the detailed analysis of initiation and propagation of particular cracks, certain damage factor, or more such factors, are evaluated to evaluate the partial lost of material stiffness. Both approaches can be coupled, even for finite deformations and dynamic fracture, relying on certain non-local regularization of material properties and on advanced crack tracking procedures, cf. [63,64].

2. Materials and Methods

Three groups of materials where this procedure can be applied and the behavior of these materials can be characterized by an unambiguous traction separation law were chosen as test materials for the use of the cohesive approach:

• forget steel;
• ceramics with fibres;
• cementitious composites.

2.1. Forget Steel

If we want to use the possibilities of cohesive approaches for modelling the behaviour of materials, it is necessary to have the following input data available based on experiments: (i) The standard true stress - true strain curve obtained by tensile tests. (ii) Allow crack growth including selection of initialization criteria. (iii) Traction separation law representing micromechanism processes. However, it can also be used to analyse certain wider aspects of material or structural behaviour, which can be verify using experiments. The individual parts of the curve show the switching on of individual mechanisms that are dominant at a given moment for the response of materials to external conditions. In a numerical model the cohesive elements are surrounded by the classical elements. When the cohesive elements are damaged the crack is extended over the boundaries of the classical elements. Thus, a potential crack can only propagate where we used cohesive elements and nowhere else. Shape of the traction-separation law is defined according to many authors for various types of materials and damage as can be seen in Figure 2, prominent localization of plastic damage for forged steel is in Figure 3, examples of GTN modelling can be seen in Figure 4, whereas and Figure 5 was used for determination of basic model parameters.

In practice, the cohesive model uses two types of finite elements. The first type describes a standard continuum; the second type specifies the connection between these areas. The decoupling of cohesion elements is determined from the displacement of standard elements, as can be seen in Figure 6. However, in the case of composite materials modelled using damage mechanics, knowledge of the decisive micro-mechanisms affecting damage is essential, which leads to realistic predictions corresponding to experimental results. The direction of crack propagation is also important, in the case of crack propagation in the direction perpendicular to the strengthening fibres, the damage is then determined by these following micromechanisms: i) damage of matrix, ii) delamination of an interface fibre versus matrix, iii) fibre cracking, iv) fibre pull out including bridging. However, it is not easy to describe the behaviour of the fibre-matrix interface.

The 42CrMo4 steel, which was produced in the form of a forging with a cross-section of 500 × 500 mm, was chosen as an experimental material for the area of upper threshold values. The experiments were mainly focused on obtaining the fracture mechanical characteristics, especially of the surface layers of the forging. All the results of experimental works were obtained on test specimens, the longest dimension of which was always oriented parallel to the directions of the fibres in the initial semi-finished product. Due to the assumed calculations using the FEM, hardness tests were carried out with the aim of determining the gradient of properties and microstructures in the surface layers of the forging. In the case of numerical modelling, which is focused on monitoring crack growth, it is essential to describe the material as precisely as possible, especially in front of the crack tip. The observed hardness gradient in this area is not significant. For the above reasons, therefore, for the following calculations the material was considered homogeneous and isotropic. The necessary dependences of true stress versus true strain were obtained from tensile tests. Determination of individual parameters for the GTN model, relevant verification procedures for cohesion modelling and various approaches are listed in [65,66,67,68,69].

2.2. Ceramics with Fibres

Modern ceramics such as ceramics based on silicon nitride or SiC are resistant to high temperatures and have high abrasion resistance. The influence of grain bridging on the strength and fracture toughness was found, so the grains behave as short reinforcing fibres. The prediction of crack propagation through interface elements based on the fracture mechanics approach and cohesive zone model is investigated from various damage models.

A transparent matrix with carbon fibres covered with a layer of SiC was chosen for testing having final diameter of 0.15 mm. The matrix - polycrystalline transparent magnesium aluminate “spinel” ceramics were obtained by hot-pressing and hot isostatic pressing. The Al${}_2$O${}_3$ content of the spinel expressed as MgO.nAl${}_2$O${}_3$ with $n=1.5$.

For the following experiments a plate with a thickness of 5 mm and a square dimension of 50 x 50 mm was made. The reinforcing fibres were placed into the middle plane of the plate periodically with nominal fibre to fibre in distance 4 mm. The position of the fibres was very important because the experimental material was cut into rectangular shapes with respect to the own position of these fibres. Three configurations were used: samples without fibres (A), made in the same way as samples containing fibres, samples with fibres symmetrically placed in front of the tip with a sharp chevron notch (B), and finally samples with a fibres in the middle of the chevron notch (C); the properties of these microstructures are easy to spot as fracture toughness result in Figure 7, images of microstructure in Figure 8.

For the model case of crack propagation, the ideal direction is perpendicular to the reinforcing fibres, when practically all basic damage mechanisms are activated as are matrix cracking, delamination of interface fibres (grains) and matrix, fibre (grain) cracking and pull out. Under extreme operating conditions as high temperatures and corrosive environments, rolling are, where elements are subjected to high cyclic contact stresses, the modelling and life prediction is more complicated. It is well recognised that the properties of ceramics can be enhanced by suitably tailoring the microstructure based on realistic application and working conditions. Namely, the effects of tailoring the grain structure on the fracture toughness of silicon nitride were demonstrated in [70]. Similarly, the influence of boundary phase manipulation and the effect of grain bridging on the strength and toughness were illustrated in [71]. To understand crack growth mechanisms in ceramic materials resistance curve (R-curve) if often used, it helps to apply the given type of ceramics in a targeted manner in technical practice by [72]. However, more recently it has been determined that some mechanisms can not be adequately accounted for the steep-rising crack growth resistance curves (R-curves). The development of advanced methods for determining the bridging stresses can be documented in [73,74,75].

2.3. Cementitious Composites

Fibre-reinforced cementitious composites have appeared in the construction industry in recent years as modern materials (e.g. concrete reinforced with steel fibres), preventing tensile stress and deformations as sources of unwanted micro- and macro-cracking. The mechanical properties of such composites are determined by choice of volume fraction of the matrix and fibres; all this is dependent to technological procedures, e.g. special compaction and eventual hardening of the matrix and some procedures (e.g. compaction) and to the early age treatment, as well as by the bond - slip interface relations – cf. [76,77,78,79,80]. We can determine the properties of the resulting structure using both destructive and non-destructive methods. taken from the early-age matrix or obtained from the crushed part of the existing structure.

The volume fraction of particles can be evaluated accurately while any information related to particle orientation is missing. Moreover, such experiments with many constructions are not allowed by technical standards. This is obviously the impetus for using some reliable non-destructive or (at least) semi-destructive methods measuring, for example, the orientation of fibres or determining the proportion of fibres, usable in situ, manipulation of homogeneity and isotropy and determining other parameters and the overall structure of the material.

A crucial problem of the last few decades is the development of a computational model for the complex hydro-thermo-chemo-mechanical behaviour of cement-based composites, where the cement particles are activated by water and subsequently form the structure of the cement matrix. The rate of cement hydration is determined by diffusion; all the components of this building material react with each other, creating significant heat that is redistributed. However, this means that subsequently the deformation and stress and volume of the material are constantly changing, without the influence of external mechanical load. Every cement particle can be assumed to initiate hydration from the moment that it comes into contact with water. Water diffuses through the hydrate layer and reaches the surface cement particles and and initiates a chemical reaction with the cement.

And where most problems can be found of computational modelling of early-age behaviour of cement-based materials? Their origin is certainly in the non-periodic structure of the material, created as a result of relatively complicated chemistry response, its different mechanisms of deformation in stress and tension, more in [81,82], including a high probability of fracture, imprecisely defined boundary conditions according to the external environment and the necessity of coexistence with normal or pre-stressed reinforcement or other reinforcement particles, see Figure 9 and Figure 10.

2.4. Numerical Adjustment

As has already been mentioned several times in this text we think over a body having crack perpendicular to the direction of oriented fibres. In the case of fibre composites, where bridging stress ${\sigma }_{br}$ is the dominant mechanism determining crack propagation and arrest, a relation between the bridging stress ${\sigma }_{br}$ and fracture energy, then by means of fracture mechanics can predict crack growth and propagation. The law is identical in every point of bridging zone. Under shock loading, fibres are damaged, so the existence of a characteristic opening ${\delta }_0$ is assumed, which determines the moment when the bridging effect disappears.

To put it simply, it is energy, we have to act on body to crack growth which can be occurred at a given moment. The bridging law (e.g. in the case of materials with fibres) is accepted such as another material characteristic

$${\sigma }_{br}=\frac{\mathsf{\Delta }{J}_{ss}}{2{\delta }_0}\sqrt{\frac{{\delta }_0}{\delta }}.$$

Microstructural properties can be characterized using a geometric factor representative volume element (RVE) determining the range of action of the given mechanisms and thus also its representation using a cohesive element. Crack propagation through the element is controlled by the fracture mechanic and by the cohesive model. This model is simpler than the classical models and its parameters are given experimentally [83,84,85,86,87].

Several steps are necessary to do due to the numerical implementation: (i) remove of the stress singularity at ${\delta }_u=0$ and (ii) incorporate the initial fracture strength $J_0$ to improve traction separation law as can be seen in Figure 6, $J_{SS}$ represents the total fracture energy without $J_0$. The following parameter ${\sigma }_0$, see Figure 6, indirectly describes peak used in traction separation law

$$K_{\mathrm{fac}}=\frac{J_0}{{\sigma }_0\mathsf{\Delta }{u}_2}\left(\alpha +\frac{1}{\alpha }\right),$$

(6)

where parameter $\alpha$ is in the range $[1,100]$. An element consists of two quadratic line elements for 2D planar elements or two quadratic planar elements for 3D. To indicate the determination of the new cohesive element and to explain its formation, we will stick to 2D. The two surfaces of the contact element initially lie on top of each other in the initial phase and are in an unloaded deformation state. The relative displacements of the element surfaces produce normal and shear displacements depending on the constitutive equation. Now assume a quadratic linear element for 2D simulations. This element has 12 degrees of freedom and the nodal displacement vector is given by (7), which generates the final stiffness matrix 12×12 for 2D, 48×48 in the generalized case for 3D. A convenient design of the cohesive element is initiated in the following set of equations (8), (9), (10), (11), (12), (13).

Let us consider the displacements (doubled) between adjacent nodes in cohesive elements

$$d_N=[D_x^1,D_y^1,D_x^2,D_y^2,....D_x^6,D_y^6].$$

(7)

Appropriate 3D cohesive element has 48 degrees of freedom. For element numbering the standard conventions is used, then the opening of the connection element is determined as a difference in displacements between upper and lower nodes:

$$\mathsf{\Delta }u=u^{UPP}-u^{LOW}.$$

(8)

As follows, the opening between adjacent layers belonging to paired nodes is necessary to define:

$$\mathsf{\Delta }{\mathit{u}}_N=\mathsf{\Phi }{\mathit{d}}_N=[-I_{6x6}|I_{6x6}]{\mathit{d}}_N$$

(9)

where $I_{6x6}$ is unity matrix having 6 rows and 6 columns, ${\mathit{u}}_N$ is a line vector with 6 columns. By default, the opening is interpolated to the integration points using shape functions. Let $N_i\left(\xi \right)$ be the shape function for node pair $i\in \{1,2,3\}$ where $\xi$ determines position in the local coordinate system $-1<\xi <1$. The relative displacement between the nodes within the elements is then given as

$$\mathsf{\Delta }u\left(\xi \right)=\left[\begin{array}{c}\mathsf{\Delta }{u}_x\left(\xi \right)\\ \mathsf{\Delta }{u}_y\left(\xi \right)\end{array}\right]=H\left(\xi \right)\mathsf{\Delta }{\mathit{u}}_N$$

(10)

where $H\left(\xi \right)$ is special matrix having 2 rows and 6 columns containing the quadratic shape function. For 2D element the shape of this matrix is following:

$$H\left(\xi \right)=\left[\begin{array}{cccccc}{N}_1\left(\xi \right)& 0& N_2\left(\xi \right)& 0& N_3\left(\xi \right)& 0\\ 0& N_1\left(\xi \right)& 0& N_2\left(\xi \right)& 0& N_3\left(\xi \right)\end{array}\right].$$

(11)

As a result, we get

$$\mathsf{\Delta }\mathit{u}\left(\xi \right)=H\left(\xi \right)\mathsf{\Phi }{\mathit{d}}_N=B\left(\xi \right){\mathit{d}}_N$$

(12)

where $B\left(\xi \right)$ has 2 rows and 12 columns and $\mathsf{\Delta }u\left(\xi \right)$ only 2 rows and 1 column. In the case of using large deformations, which is common when solving practical problems, local deformations in both standard directions are also calculated. Then it is necessary to use the centre nodes in both layers. The coordinates of the default configuration are given by the vector $x_N$ and the deformation state is predicted by the vector ${\mathit{d}}_N$; the coordinates of the reference surface $x_N^R$ are then calculated by linear interpolation:

$$x_N^R=\frac{1}{2}\left[I_{6x6}\right|I_{6x6}](x_N+{\mathit{d}}_N).$$

(13)

It is widely known that the mechanical strength and crack resistance of fibre composites is high. Constructions made of these materials are then finer and the use of high-strength steels is significantly lower. Even damage resistance is often lower. Final damage is often initiated from surface or near-surface cracks, which fibers often prevent. A crucial problem is the design of such experiments, which are macroscopic in nature, but intended to estimate the microscopic behaviour of the material, e.g. random or deliberately preferred fibre directions. Some non-invasive procedures may be the way to apply for it, as are especially radiographic, tomographic and electromagnetic methods, working with a stationary magnetic field or with a harmonic electromagnetic field.

The procedure for using XFEM for time problems-dependent is detailed in [88,89,90]. This method has a number of modifications, even existing under different names; progress in recent years can be seen by comparing works [91] and modification via the boundary element method by [92] or more in [93]. A procedure for randomly oriented fibres is presented in [94], an example of a real distributions can be found in case of reinforced concrete. However, the models describing the formation of microcracks [95,96] when the stress field is modified by the non-local Eringen model [97] come out differently. The relationship between the matrix and the fibres realizing the microscopic behaviour of the material can be described by a cohesive model, more in [98,99,100,101].

2.5. Computational Algorithms and their Convergence Properties

In all sections above (including this one), our attention was paid to various kinds of geometric and material characteristics connected with the initiation and development of micro- and macrofracture. For reasonable computations, such characteristics must be inserted some variational principles, as those of classical thermomechanics, compatible with [20]: namely of mass, momentum and energy. From this point of view, to introduce a complete 3D model problem, let us consider a deformable body occupying a domain $\mathsf{\Omega }$ in the 3-dimensional Euclidean space ${\mathcal{R}}^3$, supplied with a Cartesian coordinate system $x=(x_1,x_2,x_3)$, with the Lipschitz boundary $\partial \mathsf{\Omega }$, consisting of 2 disjoint parts $\mathsf{\Theta }$ (supports) and $\mathsf{\Gamma }$ (with active surface loads). The displacements $u_i$ of $x_i$ for $i\in \{1,2,3\}$, i. e. of all points of $\mathsf{\Omega }$, are redistributed in time t from a finite time interval $\mathcal{I}$; it is useful to relate them to the reference configuration of $\mathsf{\Omega }$ for the starting time $t=0$; upper dots will then denote partial derivatives with respect to t. Displacement rates $v=\dot{u}$ are needed, assuming the Cauchy initial conditions $u\left(0\right)=u_0$ for some $u_0$ from the standard Sobolev space of test functions $\mathcal{V}=\{w\in W^{1,2}{\left(\mathsf{\Omega }\right)}^3:w=(0,0,0)\mathrm{on}\mathsf{\Theta }\}$. The fully dynamic (not only quasi-static) formulation works also with accelerations $a=\dot{v}$ will the additional Cauchy initial conditions $u\left(0\right)=v_0$, with $v_0\in \mathcal{V}$ again. In the small deformation theory we are allowed to introduce the symmetric strain tensor ${\epsilon }_{ij}\left(w\right)=(\partial w_i/\partial x_j+\partial w_j/\partial x_i)/2$ for each $i,j\in \{1,2,3\}$ . Taking some volume loads $f\in L^2{\left(\mathsf{\Omega }\right)}^3$ and surface loads $g\in L^2{\left(\mathsf{\Gamma }\right)}^3$ into account, variable in $t\in \mathcal{I}$, the conservation of linear momentum, compatible with [102], reads

$$(w,\rho a)+\left[\epsilon \right(w),\sigma (u,v\left)\right]=(w,f)+\langle w,g\rangle$$

(14)

for all $w\in \mathcal{V}$ and $t\in \mathcal{I}$ where the symbols $(.,.)$, $[.,.]$ and $\langle .,.\rangle$ can be understood as scalar products in $L^2{\left(\mathsf{\Omega }\right)}^3$, $L^2{\left(\mathsf{\Omega }\right)}^{3\times 3}$ and $L^2{\left(\mathsf{\Gamma }\right)}^3$, $\sigma (u,v)$ being certain stress tensor. Our aim is to find $u\in \mathcal{V}$ for any $t\in \mathcal{I}$ and corresponding v, or also a, respectively.

A constitutive equation is presented as the dependence of $\sigma$ on u and / or v here. In linear elasticity we can take $\sigma =C\epsilon \left(u\right)$ where C is a 4th order tensor of stiffness, containing (thanks to all symmetries) 21 independent material characteristics in general; in particular, for an isotropic medium the set of these characteristics reduces to 2 Lamé factors (or to the Young modulus and the Poisson ratio). Some energy dissipation which can be implemented using structural damping in the Kelvin parallel viscous form $\sigma =C\epsilon \left(u\right)+\alpha C\epsilon \left(v\right)$, containing a damping factor $\alpha$, here considered as scalar for simplicity. Alternatively, mass damping takes $\rho a+\beta \rho v$ in (14) instead of $\rho a$, working with another damping factor $\beta$. The risk of cracking in tension (unlike compression) can be handled applying certain damage factor $\mathcal{D}$, compatible with [58], scalar in the simplest case, with values between 0 and 1, depending on certain invariants of $\sigma$ and / or $\epsilon \left(u\right)$ (to guarantee the objectivity). This motivates the modification of (14) to the form

$$(w,\rho (a+\beta v\left)\right)+\left[\epsilon \right(w),(1-\mathcal{D}\left)C\right(\epsilon \left(u\right)+\alpha \epsilon \left(v\right)\left)\right]=(w,f)+\langle w,g\rangle .$$

(15)

As discussed by [102], to guarantee the existence of $u,v,a$ in appropriate Bochner-Sobolev spaces, based on the method of discretization in time and on the properties of several types of Rothe sequences, some regularization for the evaluation of $\mathcal{D}$ is necessary, as that using the Eringen approach by [97]; for more details to regularization see [103]. More complicated approaches can be classified by their evaluation formulae of scalar or more complicated $\mathcal{D}$ and by their use of additional internal variables. Similar results can be derived even with i) $\gamma \nabla (\vartheta -{\vartheta }^{\mathrm{ref}})$ added to $\epsilon \left(u\right)$ in (15), working with a thermal expansion coefficient $\gamma$, $\vartheta$ being the absolute Kelvin temperature, a priori known here, supplied by some reference value ${\vartheta }^{\mathrm{ref}}$, and ii) $\gamma \nabla \dot{\vartheta }$ added to $\epsilon \left(v\right)$. Moreover, for modelling macroscopic cohesive interfaces, implementing XFEM (or similar) techniques in computational algorithms, in addition to the last right-hand-side additive term of (15), a contribution of a similar term introduced on internal interfaces (instead of $\mathsf{\Gamma }$) can be considered, with the carefully defined analogue of g dependent on jumps of u and / or v on interfaces.

Another approach to such coupled failure model is presented by [104]: it incorporates both smeared and discrete crack models through an energy-consistent thermodynamic framework. However, for the detailed development of damage a less simplified model can be needed, with the evaluation of power of dissipation without non-physical ad hoc assumptions. One can expect much more complicated formulae at least for the evaluation of $\sigma$ in (14), as well as still open problems in existence and convergence theory. Following [20], in accordance with [105], which can be seen as a substantial generalization of [106], the thermodynamic formulation can be based on the evaluation of power of dissipation from the specific Helmholtz free energy, or from the specific Gibbs energy, respectively. This can cover, using appropriate internal variables, following [107], a rather large class on nonlinear thermodynamically admissible models. In the simplest case of [108] the dissipation potential $\phi$ can be reduced to certain set of characteristic functions the local well-possedness of such formulation is studied in [109], the recent techniques for estimation of such potentials via crack tips are discussed by [110].

3. Results and Discussion

The solved problem is focused on a body with a crack under the assumption of a deformation-controlled failure mechanism. This task actually represents the use of a hybrid a methodology combining numerical modelling, experiment and microscopic observation into one relatively complex whole. Prediction results for both the GTN model and the use of cohesion elements were verified using the experimentally obtained so-called $J-R$ (J-integral resistance) curve, sometimes in literature known as the $J-\mathsf{\Delta }a$ where $\mathsf{\Delta }a$ refers to crack length differences, cf. Figure 11; for the critical review of $J-R$ curve analysis see [111].

Briefly, main effort is being focused on the crack propagation modelling by both methods, see Figure 11. They use the capabilities of the Abaqus and Warp3D software systems and focus more on the influence of each parameters in the models for the correct prediction of crack propagation. Like a model material chosen was forged steel, where both could be used due to ductile failure presented methods. The previously indicated question arises whether the procedure tested for another class of materials can be applied to another, where the dimensional factor is an order of magnitude lower. The are e.g. composites reinforced with SiC fibres with a glass matrix. From a size point of view, we are ready below and here the proposition of the traction separation law requires even more understanding of physics essence of ongoing events. Therefore, a user procedure was developed in the Fortran language, which has been implemented into the commercial Abaqus system.

The identification technique was chosen to study the fracture behaviour of short cracks. This procedure makes it possible to monitor almost all stages of the failure and as well in the stage of crack propagation. The detail of very short crack (having length of several grains) is shown in Figure 8. For real modelling, a certain calibration of the traction law is necessary using the following procedure: a) we will try first increase ${\sigma }_0$ and decrease $\mathsf{\Delta }{u}_0$, then minimize bridging shape close to maximum stress value; b) we will try decrease ${\sigma }_0$ and increase $\mathsf{\Delta }{u}_0$, so we maximize bridging shape close to exponential behaviour part of this curve.

The following results deserve to be emphasized:

• Fibre bridging in crack growth realizes an increase in fracture toughness associated with an increased level of strength in ceramics, similar to the interaction of grains in ductile materials.
• When modelling behaviour commercial ceramic as is Si${}_3$N${}_4$, see Figure 12, the onset and initiation of the crack length $\mathsf{\Delta }a$ is slower when we introduce the effect of bridging into the model.
• It is very likely that early real bridging starts due to numerical oscillation and the obtained values of displacements after the crack initiation are smaller, see the shape of the traction separation law in Figure 6.
• Real determination of the shape of the separation curve generates J – R prediction, at least maximum stress ${\sigma }_0$ must be determined on the base of careful experimental procedures.

As indicated in the introductory section, this type of modelling requires collaboration with careful experiments and techniques that see more into the microstructure. An example can be non-destructive testing of material structure via image processing (2D X-ray, 3D tomographic) and stationary magnetic and non-stationary electromagnetic approaches. For fibre cement composites, where fibres are almost randomly distributed, which is often the situation in the construction industry, there is control over the volume fraction and orientation of fibres so far only possible in the production of fresh fibre concrete mix. Metal fibres will of course improve the mechanical properties of concrete, especially fracture toughness, compressive strength, impact strength and durability of structures, in addition fatigue strength, too.

A cement matrix reinforced with metal fibres was chosen for the following numerical tests. In practice, the most important case is cement composites containing short and intentionally or quasi-randomly oriented steel, sometimes ceramic or polymer fibres with primary suppression of some stress components, is introduced in [100], while a more detailed mathematical formulation is in [82]. Its numerical approach relies on a modified XFEM where one can use as a criterion for the formation of a crack, the cohesive traction separation law. The results in the case of implementation of the nonlocal constitutive stress-strain relation of the integral type are very interesting, see [58].

Then attention is paid in particular to Eringen model [97] for generating the multiplicative damage factor, related quasi-static analysis. In the first step the XFEM is used, then the stress in front of the crack front is recalculated according to the non-local approach, in the entire body according to the exponential law of violation. The following Figure 13 and Figure 14 present some comparative results.

In some complicated situations, where ignoring crack tip effects, the stability of the X-FEM discretizations is trivial for open cracks but remains a challenge if we constrain the crack to be closed (i.e., the bi-material problem [91]). The demand for lightweight and high-strength materials in aerospace, automotive and marine industries has necessitated the use of fibre reinforced polymeric composites in place of metal alloys [112]. Related to the mentioned issue is the issue of complex loading, or in the case of the aviation industry, also modelling the response to cyclic loading, [113]. Recently, there have even been articles focusing on radiation damage and modelling the formation of cracks inside solar cells, e.g. in [114]. Indeed, almost unexpected is the use of cohesion elements to estimate the response of human tissue when inserting deep needles into soft tissues in [115], or for the development of bone-inspired composites in [116].

4. Conclusions

The conclusions can be characterized by following:

• The procedure for implementing the cohesion element into the FEM system was indicated.
• With the advent of fibre composites in technical practice, it is essential to be able to predict or model the behaviour of these materials. Numerical methods solve not only new or modified procedures, including the existence of solutions, but the modelling result must clearly approach reality. The problem is that many of the input data are estimated, which increases the risk of a possible wrong prediction.
• An example of introduced numerical problems is the form of the traction separation law in cohesion models.
• Talking about the fineness of the FEM network can be counterproductive, it is necessary to start from the size of the RVE (representative volume element).
• Small modifications of XFEM with a focus on the applicability of these procedures were also tested on practical examples.
• For the modeling of microstructural behaviour using XFEM, it is often necessary to use, or rather to introduce, a real traction separation law.
• Careful determination of the traction separation law representing all phases of fibre reinforced composite behaviour enables a more accurate prediction of crack propagation predominantly in the initial phase of failure.
• In the case of cement composites, it is reasonable to use models that in a certain way average the stress field in front of the crack front.
• The combination of traction separation law and XFEM is a promising approach for crack propagation modelling, as a strong motivation for further research.