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

TL;DR

This paper explores how cohesive approaches can model damage and crack growth in fiber-reinforced materials to improve structural predictions.

Contribution

The study introduces cohesive approaches combined with modified finite element methods for accurate modeling of damage in fiber composites.

Findings

Cohesive elements effectively model energy release at crack fronts in fiber-reinforced materials.

Combining cohesive approaches with modified FEM improves accuracy in predicting structural damage.

The methods show promise for practical applications in structural material analysis.

Abstract

During the operation of structures, stress and deformation fields occur inside the materials used, which often ends in fatal damage of the entire structure. Therefore, the modelling of this damage, including the possible formation and growth of cracks, is at the forefront of numerical and applied mathematics. The finite element method (FEM) and its modification will allow us to predict the behaviour of these structural materials. Furthermore, some practical applications based on cohesive approach are tested. The main effort is devoted to composites with fibres and searching for procedures for their accurate modelling, mainly in the area where damage can be expected to occur. The use of the cohesive approach of elements that represent the physical nature of energy release in front of the crack front has proven to be promising not only in the direct use of cohesive elements, but also in combination with modified methods of standard finite elements.

Figures

Full text

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 under external conditions. The field of science, summarised 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 modernise procedures or relevant standards, but also to use new ones with more accurate numerical methods, which with increasing probability can predict the behaviour of materials even under complex operating conditions. The behaviour of bodies containing defects has been described and predicted using fracture mechanics.

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], and the second one considers the accumulation of damage during the 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 in 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 uses the long-time experience of authors in this research field, 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 of existing microcracks, for example, those created during manufacturing processes, before these propagate further. The concept of the stability criterion can be extended to the case of general stress concentrators, in which case, the conditions under which the crack starts to propagate from the considered concentrator has to be selected. 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 $[eqn]$ criterion (the $[eqn]$ criterion can be used equivalently). For nonlinear fracture mechanics, the $[eqn]$ 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 $[eqn]$ integral with a special cohesive element for dynamic crack propagation can be seen in [19].

From a 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 discretisations both in the Euclidean space, 3-dimensional in general, and on certain time variables, even for a seemingly static, simplified evaluation of fracture development. The conservation principles contain both the total strain tensor $[eqn]$ and the stress tensor $[eqn]$ , both represented by symmetric square matrices of order 3. Since the values of components $[eqn]$ of $[eqn]$ with $[eqn]$ depend on the choice of a Cartesian coordinate system, for the following considerations, it is useful to introduce also three invariants of $[eqn]$ (independent of the choice of Cartesian coordinates) $[eqn]$ , $[eqn]$ and $[eqn]$ ; the first (linear) one is $[eqn]$ ; for the derivation of their complete set, see [21]. Most computational tools $[eqn]$ are decomposed into four components, denoted as $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ , referring to elastic, plastic, creep and thermal ones; moreover, an appropriate incorporation of the damage process is expected.

In particular, in the simplified small deformation theory for isotropic materials, using the standard Kronecker symbol $[eqn]$ , for purely elastic deformation, the following equation

$$[eqn]$$

can be written for all $[eqn]$ ; this relation is the famous empiric Hooke’s law, containing a couple of material parameters $[eqn]$ : E is the Young’s modulus and $[eqn]$ is the Poisson’s ratio. To derive $[eqn]$ for $[eqn]$ , unlike (1), some appropriate decomposition of $[eqn]$ must be suggested: in most computational tools, $[eqn]$ is considered as a simple sum of $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ . Nevertheless, proper nonlinear formulations, namely, for the elastoplasticity theory, see [21], and in the more general context for the theory of structured deformation by [22], revisiting the original idea in [23], require more complicated multiplicative decompositions, with non-classical material characteristics and internal variables.

The irreversible component $[eqn]$ is activated (takes non-zero values) after a certain function $[eqn]$ for principal stresses reaches some prescribed value $[eqn]$ : e.g., for metals, the von Mises criterion with

$$[eqn]$$

is used frequently, with $[eqn]$ 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 $[eqn]$ can be evaluated from the associated plastic flow rule, following [21,24]; all additional material parameters are expected to be identified from well-prepared experiments. The superposition of $[eqn]$ and $[eqn]$ relies on the Kelvin–Voigt model typically; $[eqn]$ is evaluated as a function (in general nonlinear, especially in the well-known Norton power-law model of creep) of stress components or invariants, frequently from the deviatoric part $[eqn]$ of $[eqn]$ , i.e.,

$$[eqn]$$

with $[eqn]$ again, and its Euclidean norm S satisfying

$$[eqn]$$

Most simplified evaluation of $[eqn]$ can rely on the linear dependence of $[eqn]$ on the prescribed temperature change, using certain thermal expansion coefficient(s); 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, utilised by [26,27,28,29], derives the crack propagation from the C-integral, working with energy rate. For a 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, followed 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 the accumulation of damage is processed by [36]. As an example of using damage mechanics for ductile failure, the Gurson–Tvergaard–Needleman (GTN) model can be pointed out: see [37,38]. Based on the evaluation of the performed experiments, it is recommended to introduce two (or three) optional parameters $[eqn]$ , usually $[eqn]$ , $[eqn]$ , $[eqn]$ . A careful analysis of the influence of these parameters on the development of damage was carried out by [39]. Here, $[eqn]$ is the yield stress of the matrix material, $[eqn]$ is the principal (hydrostatic) stress and $[eqn]$ 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

$$[eqn]$$

cf. (3) and (4). For practical calculations, the volume fraction of voids $[eqn]$ in (5) is usually replaced by the effective volume fraction $[eqn]$ , introduced as

$$[eqn]$$

here, $[eqn]$ is the critical volume for coalesce of voids, $[eqn]$ is the volume of voids at ultimate damage and $[eqn]$ , as needed in (5). The effective value $[eqn]$ by (6) applies to (5) when $[eqn]$ is less than $[eqn]$ .

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 it 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 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 $[eqn]$ on particular points x of the e-th element (three-dimensional vectors of functions in general) and displacements at selected element nodes $[eqn]$ , utilising some standard shape functions $[eqn]$ , i.e.,

$$[eqn]$$

here, $[eqn]$ is the set of nodes corresponding to the standard e-th element.

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

$$[eqn]$$

including certain specialised shape functions $[eqn]$ and $[eqn]$ for the intrinsic version of XFEM, unlike the extrinsic version that adopts the original functions $[eqn]$ from (7). $[eqn]$ here is a certain discontinuous 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 (8) is identical with that of (7), the second one supports the evaluation of crack formation, working with some set of nodes $[eqn]$ , and the third one is the inclusion of the local phenomena in front of the crack, working with a still other set of nodes $[eqn]$ . The well-known setting with $[eqn]$ , frequently used in the two-dimensional crack propagation modelling, relies on the choice of special functions $[eqn]$ by [48].

Figure 1a illustrates the location of areas A, B and C, a related model area around the moving crack; the sets of nodes $[eqn]$ , $[eqn]$ and $[eqn]$ of these areas occur in (7) and (8). Figure 1b shows the real crack shape in forged steel. The progress of this approach, adopted for three-dimensional modelling at the expense of much more complicated functions $[eqn]$ , in the last years can be seen from [49].

In general, the method of discretisation 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 its own names and designations) already has quite a rich history with 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, utilising an upgrade of (7), presented as IXFEM, is introduced by [57] to reach inexpensive computations with good geometrical accuracy and surface mesh resolution independence. Another approach relies on the simplification of “smeared” damage by [58,59,60,61,62]: instead of the detailed analysis of the initiation and propagation of particular cracks, certain damage factor, or more such factors, are analysed to evaluate the partial loss of material stiffness. Both approaches can be coupled, even for finite deformations and dynamic fracture, relying on certain nonlocal regularisation of material properties and on advanced crack tracking procedures, cf. [63,64]. Examples and classic procedures focused on modelling the growth, namely, the formation of crack propagation, can be found in [65]. The purpose of this paper is to show the results of a study focused on the occurrence of damage heterogeneous materials. Attention is paid to the direct application of the finite element method to different types of materials in order to find critical parameters determining behaviour of materials at the damage process.

2. Materials and Methods

Three groups of materials to which the cohesion approach can be applied were selected for the following modelling:

• Forged steel;
• Ceramics with fibres;
• Cementitious composites.

2.1. Forged Steel

If the possibilities of cohesive approaches for modelling the behaviour of materials had been used, it would have been necessary to have the following input data available based on experiments: (i) The standard true stress—true strain curve obtained by tensile tests; (ii) allowing crack growth through the selection of initialisation criteria; and (iii) traction–separation law representing micromechanic processes. However, it can also be used to analyse certain wider aspects of material or structural behaviour, which can be verified using experiments. The individual parts of the curve show the switching on of individual mechanisms that are dominant at a given moment in 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 extends over the boundaries of the classical elements. Thus, a potential crack can only propagate where cohesive elements have been used and nowhere else. The 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, general concept of cohesion approach is schematically presented in Figure 3. Figure 2b shows that the zero stress value is not used in real modelling, which improves numerical stability.

The experimental material was supplied as block-shaped castings with dimensions of 100 × 55 × 260 mm. These blocks were annealed and then cut into blanks with dimensions of 30 × 55 × 260 mm, from which test bodies for the three-point bending test and cylindrical test bodies of the circular cross-section for the static uniaxial tensile test were made. The initial measured length $[eqn]$ = 26 mm and the diameter of the tested length d = 6 mm. During the static tensile test, an extensometer MTS 634-12F-51 was used; during the test, both elongation and contraction were recorded at a crossbar speed of 2 mm per minute. The chemical composition of the material in weight percent was as follows: 0.09 C, 1.18 Mn, 0.37 Si, 0.01 P, 0.12 Cr, 0.29 Cu, 0.29 Ni, 0.03 Mo and 0.028 Al. Yield strength at room temperature was 250 MPa. Prominent localisation of plastic damage for forged steel is presented in Figure 4, examples of GTN modelling can be seen in Figure 5, and Figure 6 was used for the 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 by the displacement of standard elements, as can be seen in Figure 7. 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, and (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 work 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 [66,67,68,69,70].

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 $[eqn]$ $[eqn]$ content of the spinel is expressed as MgOn $[eqn]$ $[eqn]$ with $[eqn]$ .

For the following experiments, a plate with a thickness of 5 mm and a square dimension of 50 × 50 mm was made. The reinforcing fibres were placed into the middle plane of the plate periodically with nominal fibre to fibre at a 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; the results are presented in Figure 8 and images of the microstructure are shown in Figure 9.

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, such as high temperatures and corrosive environments, where elements are subjected to high cyclic contact stresses, the modelling and life prediction are 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 [71]. Similarly, the influence of boundary phase manipulation and the effect of grain bridging on the strength and toughness were illustrated in [72]. To understand crack growth mechanisms in ceramic materials, resistance curve (R-curve) is often used, and it helps to apply the given type of ceramics in a targeted manner in technical practice by [73]. However, more recently, it has been determined that some mechanisms cannot adequately account for the steep-rising crack growth resistance curves (R-curves). The development of advanced methods for determining the bridging stresses is documented in [74,75,76].

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 the 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. [77,78,79,80,81]. The properties of the resulting structure using both destructive and non-destructive methods could be determined 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 obtaining information related to particle orientation is complicated. 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 of 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. However, this means that subsequently, the deformation, 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 initiates a chemical reaction with the cement.

And where can most problems be found in the 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 as well its different mechanisms of deformation in stress and tension, more in [82,83], 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 10 and Figure 11. For more information on cementitious composites for engineering applications see the article [84]; for mechanical properties of fibre concrete, see [85]. The article focused on the preparation of composites with magnetically aligned metal fibres [86] and the authors’ practical application shows the way to the preparation of construction elements with directional properties.

2.4. Numerical Adjustment

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

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

$$[eqn]$$

Microstructural properties can be characterised 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 [87,88,89,90,91].

Several steps are necessary due to the numerical implementation: (i) remove of the stress singularity at $[eqn]$ and (ii) incorporate the initial fracture strength $[eqn]$ to improve the traction–separation law, as can be seen in Figure 7. $[eqn]$ represents the total fracture energy without $[eqn]$ . Following parameter $[eqn]$ in (10), see Figure 7, indirectly describes the peak used in traction–separation law

$$[eqn]$$

where parameter $[eqn]$ is in the range $[eqn]$ . 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, only 2D is presented. 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 (11), which generates the final stiffness matrix 12 × 12 for 2D and 48 × 48 in the generalised case for 3D. A convenient design of the cohesive element is initiated by the following set of equations: (11) for nodal displacements, (12) for displacement jumps, (13) for nodal displacement jumps on adjacent layers, (14) for displacement jumps, via shape functions, (15) for corresponding shape functions, up to the desirable results (16) and (17).

In more detail, let us consider the displacements (doubled) between adjacent nodes in cohesive elements

$$[eqn]$$

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:

$$[eqn]$$

The opening between adjacent layers belonging to paired nodes is defined:

$$[eqn]$$

where $[eqn]$ is the unity matrix having 6 rows and 6 columns, $[eqn]$ is a line vector with 6 columns. By default, the opening is interpolated to the integration points using shape functions. Let $[eqn]$ be the shape function for node pair $[eqn]$ where $[eqn]$ determines position in the local coordinate system $[eqn]$ . The relative displacement between the nodes within the elements is then given as

$$[eqn]$$

where $[eqn]$ is special matrix with 2 rows and 6 columns containing the quadratic shape function. For the 2D element, the shape of this matrix is as follows:

$$[eqn]$$

We obtain

$$[eqn]$$

where $[eqn]$ has 2 rows and 12 columns and $[eqn]$ 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 $[eqn]$ and the deformation state is predicted by the vector $[eqn]$ ; the coordinates of the reference surface $[eqn]$ are then calculated by linear interpolation:

$$[eqn]$$

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 fibres 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-dependent problems is detailed in [92,93,94]. This method has a number of modifications, even existing under different names; progress in recent years can be seen by comparing the work of [95] and modification via the boundary element method by [96] or more in [97]. A procedure for randomly oriented fibres is presented in [98]; an example of a real distribution can be found in the case of reinforced concrete. However, the models describing the formation of microcracks [99,100] when the stress field is modified by the nonlocal Eringen’s model [101] come out differently. The relationship between the matrix and the fibres realising the microscopic behaviour of the material can be described by a cohesive model, more in [102,103,104,105].

2.5. Computational Algorithms and Their Convergence Properties

In all sections above (as well as this one), our attention is paid to various kinds of geometric and material characteristics connected with the initiation and development of micro- and macrofractures. For reasonable computations, such characteristics must be inserted with some variational principles, such as those of classical thermomechanics, compatible with [20]: namely, 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 $[eqn]$ in the 3-dimensional Euclidean space $[eqn]$ , supplied with a Cartesian coordinate system $[eqn]$ , with the Lipschitz boundary $[eqn]$ , consisting of 2 disjoint parts $[eqn]$ (supports) and $[eqn]$ (with active surface loads). The displacements $[eqn]$ of $[eqn]$ for $[eqn]$ , i.e., of all points of $[eqn]$ , are redistributed in time t from a finite time interval $[eqn]$ ; it is useful to relate them to the reference configuration of $[eqn]$ for the starting time $[eqn]$ ; the upper dots will then denote partial derivatives with respect to t. Displacement rates $[eqn]$ are needed, assuming the Cauchy initial conditions $[eqn]$ for some $[eqn]$ from the standard Sobolev space of test functions $[eqn]$ . The fully dynamic (not only quasi-static) formulation works also with accelerations $[eqn]$ with the additional Cauchy initial conditions $[eqn]$ , with $[eqn]$ again. In the small deformation theory, we are allowed to introduce the symmetric strain tensor $[eqn]$ for each $[eqn]$ . Taking some volume loads $[eqn]$ and surface loads $[eqn]$ into account, variable in $[eqn]$ , the conservation of linear momentum, compatible with [106], reads

$$[eqn]$$

for all $[eqn]$ and $[eqn]$ where the symbols $[eqn]$ , $[eqn]$ and $[eqn]$ can be understood as scalar products in $[eqn]$ , $[eqn]$ and $[eqn]$ , $[eqn]$ being certain stress tensor. The aim is to find $[eqn]$ for any $[eqn]$ and corresponding v, or also a, respectively.

A constitutive equation is presented as the dependence of $[eqn]$ on u and/or v here. In linear elasticity, we can take $[eqn]$ 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’s modulus and the Poisson’s ratio). Some energy dissipation which can be implemented using structural damping in the Kelvin parallel viscous form $[eqn]$ , containing a damping factor $[eqn]$ , is here considered as a scalar for simplicity. Alternatively, mass damping takes $[eqn]$ in (18) instead of $[eqn]$ , working with another damping factor $[eqn]$ . The risk of cracking in tension (unlike compression) can be handled applying certain damage factor $[eqn]$ , compatible with [58], scalar in the simplest case, with values between 0 and 1, depending on certain invariants of $[eqn]$ and/or $[eqn]$ (to guarantee the objectivity). This motivates the modification of (18) to the form

$$[eqn]$$

As discussed by [106], to guarantee the existence of $[eqn]$ in appropriate Bochner–Sobolev spaces, based on the method of discretisation in time and on the properties of several types of Rothe sequences, some regularisation for the evaluation of $[eqn]$ is necessary, as that using the Eringen’s approach by [101]; for more details to regularisation, see [107]. More complicated approaches can be classified by their evaluation formulae of scalar or more complicated $[eqn]$ and by their use of additional internal variables. Similar results can be derived even with (i) $[eqn]$ added to $[eqn]$ in (19), working with a thermal expansion coefficient $[eqn]$ , $[eqn]$ being the absolute Kelvin temperature, a priori known here, supplied by some reference value $[eqn]$ , and (ii) $[eqn]$ added to $[eqn]$ . 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 (19), a contribution of a similar term introduced on internal interfaces (instead of $[eqn]$ ) can be considered, with the carefully defined analogue of g dependent on jumps of u and/or v on interfaces.

Another approach to such a coupled failure model is presented by [108]: 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 $[eqn]$ in (18), as well as still open problems in existence and convergence theory. Following [20], in accordance with [109], which can be seen as a substantial generalisation of [110], the thermodynamic formulation can be based on the evaluation of the power of dissipation, from the specific Helmholtz free energy, or from the specific Gibbs energy, respectively. This can cover, using appropriate internal variables, following [111], a rather large class on nonlinear thermodynamically admissible models. In the simplest case of [112], the dissipation potential $[eqn]$ can be reduced to a certain set of characteristic functions. The local well-posedness of such formulation is studied in [113], and the recent techniques for estimation of such potentials via crack tips are discussed in [114].

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 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 $[eqn]$ (J-integral resistance) curve, sometimes in literature known as the $[eqn]$ , where $[eqn]$ refers to crack length differences, cf. Figure 12; for the critical review of $[eqn]$ curve analysis, see [115].

Briefly, the main effort is being focused on the crack propagation modelling by both methods, see Figure 12. They use the capabilities of the Abaqus and Warp3D software systems and focus more on the influence of each parameter in the models for the correct prediction of crack propagation. Like the model, the 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 one class of materials can be applied to another where the dimensional factor is an order of magnitude lower, e.g., composites reinforced with SiC fibres with a glass matrix. From a size point of view, the proposition of the traction–separation law requires even more understanding of the physics 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 a very short crack (having a length of several grains) is shown in Figure 9. For real modelling, a certain calibration of the traction law is necessary using the following procedure: (a) Try first to increase $[eqn]$ and decrease $[eqn]$ , then minimise the bridging shape close to the maximum stress value. (b) Try to decrease $[eqn]$ and increase $[eqn]$ , so maximise the bridging shape close to the exponential behaviour part of this curve.

The authors present their long-term experience with the issue and compare it with the state of knowledge in recent years. The following results deserve to be emphasised:

• Fibre bridging in crack growth realises 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 the behaviour of commercial ceramic as $[eqn]$ $[eqn]$ , see Figure 13, the onset and the initiation of the crack length $[eqn]$ is slower when the effect of bridging into the model is introduced, which is closer to an experiment.
• 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 7.
• Real determination of the shape of the separation curve generates J – R prediction, at least the maximum stress $[eqn]$ must be determined on the basis 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 [104], while a more detailed mathematical formulation is in [83]. 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 the 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’s model [101] for generating the multiplicative damage factor, related to quasi-static analysis. In the first step, the XFEM is used, then the stress in front of the crack is recalculated according to the nonlocal approach in the entire body according to the exponential law of violation. The following Figure 14 and Figure 15 present some comparative results.

For most numerical tests, a slightly modified method of the extended finite element method was used; this method uses new degrees of freedom for crack propagation within the standard element (extrinsic XFEM). The test material, in our case, a cement paste reinforced with metal fibres, had the following material characteristics: Young’s modulus E = 3.2 GPa, Poisson’s ratio $[eqn]$ = 0.3 and tensile strength 12 MPa. Using a nonlocal approach to calculate the stress value in front of the crack will significantly change the crack propagation rate in the model case. Accelerating the damage with a damage parameter higher than 0.2 creates the preconditions for multiple small cracks around the stress concentrator.

In some complicated situations, ignoring crack tip effects, the stability of the XFEM discretisations is trivial for open cracks but remains a challenge if we constrain the crack to be closed (i.e., the bi-material problem [95]). The demand for lightweight and high-strength materials in the aerospace, automotive and marine industries has necessitated the use of fibre-reinforced polymeric composites in place of metal alloys [116]. 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, [117]. Recently, there have even been articles focusing on radiation damage and modelling the formation of cracks inside solar cells, e.g., in [118]. Indeed, almost unexpected is the use of cohesion elements to estimate the response of human tissue when inserting deep needles into soft tissues in [119], or for the development of bone-inspired composites in [120].

4. Conclusions

The conclusions can be characterised by the 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. Not only must numerical methods solve new or modified procedures, including the existence of solutions, but the modelling result must also 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 numerical problems introduced 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 with 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 modelling 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.
• The combination of the traction–separation law and XFEM is a promising approach for crack propagation modelling and is a strong motivation for further research.
• Three groups of materials were tested: ductile steel, ceramic composites and cementitious composites reinforced with metal fibres. From a microstructural point of view, the materials have different RVE sizes but the procedures and modelling of crack initiation and propagation using XFEM and using the cohesive approach are very similar.
• The outlined path is probably promising for the future higher industrial deployment of these fibre composites. For example, in the field of building composites, the production of building components with directed fibres in the given component is expected, which will enable production with clear, predictable properties.
• The theory of sharp stress concentrators in front of the crack tip for the named types of materials is abandoned, and so-called smeared models with averaged stress are promising and could be a subject of future research, as presented in this paper.

Some more comments: Remarkable progress in engineering in the last decades is connected with the use of advanced materials, whose design for implementation in engineering structures cannot rely on long-time experience, justifying the simplified computational formulae. A proper physical formulation covering elastic, viscous and plastic deformation including damage for a sufficiently wide class of materials— including those three discussed in Section 2, taking all available information about material structure into account, together with relevant additional physical, chemical, etc. processes— is still a research challenge for the near future. It should be accompanied by the mathematical analysis of existence, uniqueness, smoothness and other important properties of solutions of related evolutionary problems for partial differential equations, as well as by the development of robust and effective computational algorithms, reflecting the progress in computer hardware and software, as in parallel and distributed systems, artificial intelligence, etc.