A specialised finite element for simulating self-healing quasi-brittle materials

A new specialised finite element for simulating the cracking and healing behaviour of quasi-brittle materials is presented. The element employs a strong discontinuity approach to represent displacement jumps associated with cracks. A particular feature of the work is the introduction of healing into the element formulation. The healing variables are introduced at the element level, which ensures consistency with the internal degrees freedom that represent the crack; namely, the crack opening, crack sliding and rotation. In the present work, the element is combined with a new cohesive zone model to simulate damage-healing behaviour and implemented with a crack tracking algorithm. To demonstrate the performance of the new element and constitutive models, a convergence test and two validation examples are presented that consider the response of a vascular self-healing cementitious material system for three different specimens. The examples show that the model is able to accurately capture the cracking and healing behaviour of this type of self-healing material system with good accuracy.

An alternative approach to the nodal based enrichment of XFEM, GFEM and Cut-FEM is to use an element-based enrichment, EFEM. In EFEM, the additional degrees of freedom can be eliminated at the element level using static condensation [49], which leads to a global system of equations comprising only the original degrees of freedom. A comparative study on the two approaches was presented by Oliver et al. [50], who found that EFEM was computationally more efficient than XFEM. Ortiz et al. [51] presented the first work incorporating a weak discontinuity to capture shear bands, whilst Dvorkin et al. [52] presented the first model incorporating a strong discontinuity to capture strain localisation. There are now a range of EFEM models that can be found in the literature [17,49,[51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67], with significant developments having been made. A full review of the first decade of research on EFEM is provided in Jirasek [68].
A formulation that allowed for non-constant displacement jumps in each element and a global or local approach, depending on whether or not static condensation was employed, was presented in Alfaiate et al. [49]. The former ensures continuity of the displacement jumps, at the cost of increasing the size of the global system of equations, whilst the latter requires the jumps to be averaged at coincident crack nodes. In Dias-da-Costa et al. [58][59][60] the generalised strong discontinuity approach GSDA was presented. The model employed a global approach and allowed for both non-constant displacement jumps and stretching. Linder and Armero [61,62] developed a new embedded discontinuity element that allowed for branching cracks in a T-shape pattern. The new element was employed in the simulation of a dynamic fracture benchmark problem and showed good agreement with experimental and numerical results. Dynamic effects were also considered in the work of Saksala et al. [64,65], who developed a 3D model combining continuum visco-damage for the pre-peak response with the embedded discontinuity model for the post-peak response. The model is rate-dependent and employs an explicit time integration to boost computational efficiency. Lu et al. [67] presented a model combining a multi-scale finite element method (MsFEM) [69] with the EFEM for discontinuities in saturated porous media (see also Lu et al. [66] for modelling cracks in solids). The model combines the ability of EFEM to capture localisations with the computational efficiency of MsFEM. Numerical examples were presented which showed that the proposed method gave accurate predictions when compared with EFEM, with a greater computational efficiency. A model for self-healing quasi-brittle materials was presented by Zhang and Zhuang [17]. The authors utilised an EFEM approach that allowed for crack opening and crack sliding. The model used a reduced integration for cracked elements and smeared the displacement jump across the crack band. The healing was introduced at the constitutive level into the traction-separation law, with an Arrhenius-type law being employed for the solidification of the healing agent.
In this study, the EFEM approach is employed and a specialised element is presented that takes healing into account at the element level. The element-based strong-discontinuity approach for simulating both cracking and healing is preferred to the smeared approach of Zhang and Zhuang [17] because the former provides a more accurate description of the crack geometry and crack displacements, which is particularly important in the present work, in which the mechanical model is coupled to a transport model for the flow of healing agents through the discrete cracks [19]. The approach adopted for this work ensures that the healing is consistent with the internal degrees of freedom introduced to represent the crack and allows for a variational derivation. A cohesive zone approach with a damage-healing constitutive model is employed, whilst equilibrium over the crack plane is satisfied in a weak sense using a local variation in the internal degrees of freedom. The element allows for different healing mechanisms and for overlapping cycles of damage and healing when combined with a suitable crack plane model, such as that employed in the present work [70].
The principal contributions of this work are, (i) a new element with healing degrees of freedom linked to a unique damage-healing constitutive model, (ii) the implementation of the element within a new coupled transport-mechanical computational framework, and (iii) validations using new experimental data and idealisations.
The layout of the remainder of this paper is as follows: • "Theoretical basis for the element" section describes the theoretical formulation of the element including the introduction of healing: • "Healing agent transport and curing model" section describes the healing agent transport and curing model employed in this study: • "Crack-plane constitutive model" section describes the crack-plane constitutive model employed in this study: • "Example problems" section presents a convergence test and two example problems concerning a self-healing concrete with a vascular healing system filled with cyanoacrylate, the first being based on a direct tension test and the second an L-shaped specimen: • "Conclusions" section gives some conclusions.

A quadrilateral finite element with an embedded strong discontinuity
The aim of this section is to present the formulation of a quadrilateral finite element (Q4) capable of describing embedded crack phenomena governed by any thermodynamically valid inelastic constitutive relationship. The starting point will be the general Q4 bilinear quadrilateral element to which is added a straight discontinuity, which crosses the element at an angle ψ , as illustrated in Fig. 1.
The problem domain is defined as Ω ∈ R 2 with boundary Γ . The total virtual work ( δΠ ) for this domain may be written as follows, noting that δ denotes a virtual quantity: where x is the Cartesian basis, u denotes the displacement vector; t ∈ [t 0 , T ] is the time., Γ D and Γ N are the complementary subsets of the boundary where the prescribed displacements u D (x, t) and tractions g D (x, t) are applied, ε(u) is the strain tensor, σ(ε) is the stress tensor, F B (x, t) and F ∂Ω (x, t) describe the external body and surface forces respectively, and n is a unit vector on the boundary.
If the domain Ω is divided into standard Q4 elements, it is assumed that within each elemental subdomain ( Ω e ) any displacement can be linearly interpolated from the values at the nodes (vector ū ) through the shape functions (i.e. u(x) = N(x) u e ). In this paper, direct tensor notation is used for Eq. (1) and the continuum constitutive equations but matrix notation is employed when the discretised form of equations and crack-plane constitutive relationships are considered.
A method for introducing a discontinuity into the element is now described. Let ϑ(x; x C , ψ) be the signed-distance equation describing the linear discontinuity (see Fig. 1), where x C stands for the central point of the crack within the element. Two subdomains Ω + k and Ω − k can be defined as x ∈ Ω + e ∀ ϑ(x) ≥ 0 and x ∈ Ω − e ∀ ϑ(x) < 0 . Then, using the Heaviside function H ϑ placed at this discontinuity, the elemental displacements can be expressed as the sum of two functions u e and [[u]] , the later representing the jump. In other words: In the present formulation, the continuum material is assumed to be elastic. The inelastic displacement associated with the discontinuity is extracted from a narrow band of material of finite width (See Fig. 2). This narrow band material containing the crack will be termed 'the crack-plane' , which is defined in the orthonormal (local) basis {� r, � s} , with r being normal to the crack-plane. The local crack-plane relative-displacement vector ( ũ ) is the sum of an inelastic part ( ⌣ u ), which contains the crack opening and sliding displacements, and an elastic part ( ũ e ). The total crackplane relative displacements, including those of the elastic band of material either side of the discontinuity, are illustrated in Fig. 2.  The stress in the continuum part of the element ( Ω + k ∪ Ω − k ) is required, and is given by: where D is the elasticity tensor and ε ϑ ( ⌣ w, x j ) is the strain in the continuum part of the element caused by the displacement jump at the discontinuity, noting that ⌣ w is defined below.
The second term requires the tractions and displacements across the discontinuity. In the present approach, the latter are derived from the following crack-plane vector, which contains crack opening and sliding displacement components, as well as the rotation at the mid-point of the crack within an element (i.e. at x = x C ); The opening and sliding displacement components at a position Z along the crack (see Fig. 1) are obtained from ⌣ w as follows; in which Λ Z = 1 0Z 0 1 0 and, in the present 2D case, Where, the variable Z ∈ 0.5 · [−l ϑ ,l ϑ ] spans the crack length l ϑ in the same direction as s ; and t hk is the out-of-plane thickness of the crack band.
The work conjugate vector to ⌣ w , which contains force and moment terms, is obtained by integrating the tractions along the length of the discontinuity, as follows: The tractions applied to the crack plane of material are equal to those applied across the discontinuity, thus σ = ⌣ σ and F = ⌣ F . In order to derive ε ϑ , an expression is needed for the displacement in the positive part of the elements (i.e. x ∈ Ω + ) due solely to the relative-displacement across the discontinuity (i.e. u ϑ (x) ). This is given by the following expression, which is based on the assumption that the rotation ( α ) is small (See Fig. 1): in which The operator T w (x; ψ,α, x C , H ϑ ) maps the rigid body motion induced by the discontinuity to the parent element nodes and the introduction of H ϑ ensures that vector u ϑ (x) is only non-zero in the positive part of the element. By considering only the rigid body motions, it is implicitly assumed the enhanced strain in the continuum is null.
This means that integrations over multiple element zones are avoided [59], which makes implementation more straightforward. However, this comes at the expense of restricting the shear jumps to being constant along the discontinuity in each element [59]. In some shear dominated problems, this restriction may mean that a finer mesh would be required to achieve the same level of accuracy compared to a solution using elements with a tangential relative stretching mode [63], but this issue is not believed to be significant for the problems considered in this paper.
ε ϑ is then obtained by applying Eq. (7) to the element nodal displacements u ϑ (x j )(j denoting the element local node number) and by employing the standard strain-displacement matrix B(x) derived from the element shape functions N (x) . This gives the following equation; The discontinuity is assumed to follow the basic elastic-damage cohesive crack constitutive relationship given in Jefferson et al. [71]. This is based on damage mechanics and gives the following relationship between the crack-plane tractions and relative displacements; 0k s ; k r and k s the normal and shear elastic moduli respectively for the crack-plane; ω ∈ [0, 1] is the crackplane damage variable, which depends on the damage evolution parameter ( ζ ), according to the functions given in Sect. 4. It is also noted that subscript e denotes an element quantity and superscript e refers to an elastic entity. Using Eq. (9) in Eq. (6) leads to the following relationship between the crack-plane displacement and force vectors: is the parametric coordinate along the discontinuity and J = dZ/dξ the resultant Jacobian.
The relationship between the crack-plane and inelastic displacement vectors is given by; in which I 3 is the rank 3 identity matrix. Making use of Eqs. (1), (3) and (11), the total virtual work may be written as the sum of the virtual work of the continuum part of the element and that of the discontinuity, as follows: Using (9) in (13) results in the following revised virtual work equation; Since the constitutive relationships (10) and (11) are written in terms of the total crack-plane displacement vector, it proves convenient to express (14) in terms of w rather then ⌣ w . This is accomplished by using Eqs. (9) and (12) Equilibrium between the tractions across the discontinuity and the stresses in the adjacent continuum is considered in a weak sense by assuming that the element nodal displacements are fixed and that there is a virtual change in w , which leads to; from which the following relationship between the crack-plane displacement vector and element nodal displacements can be derived; Using Eq. (18) in (16) gives; The element stiffness relationship may then be obtained by cancelling the common virtual displacements and rearranging Eq. (19), as follows: In this approach, equilibrium across the discontinuity and the crack-plane displacements are defined at the element level (see Eq. (18)). Whilst this means that there is no continuity of the crack-plane displacements across element boundaries (though the mean values at coincident nodes could be adopted in the constitutive law [49]), it ensures that the bandwidth of the global stiffness matrix remains constant [58]. We also note that the crack path used in the discrete flow computations is smoothed by averaging the crack openings from adjacent elements at common boundaries.

Extension to damage-healing
To account for healing, the local constitutive relationship (10) is modified by adding a term to account for the damaged material that has been restored. The proportion of material that is healed at time t is given by h (t) ∈ [0,ω] and the healing displacement ( ũ h ) is introduced to ensure that an increment of healing takes place under constant stress, which ensures that a healing step doesn't violate the second law of thermodynamics (See Jefferson et al. [20]). Various forms of healing model are described in Jefferson et al. [20] and a specific model (see "Healing agent transport and curing model" and "Crack-plane constitutive model" sections for an overview) is presented in Freeman and Jefferson [19] and Jefferson et al. [70]. The aim here is to show how a general form of local healing model may be incorporated into the element with an embedded strong discontinuity.
The elastic-damage-healing counterpart to Eq. (11) may be obtained from Eq. (19), as follows: in which The inelastic component of displacement is now given by: and the strain change due to the relative displacements across the discontinuity is; h . Since an increment of healing takes place under a zero change in stress, it may be assumed that;δ ⌣ w = I 3 −K e −1 K +K h δw and δε ϑ = M w δw The virtual work equation for the damage-healing is obtained by using Eqs. (22)(23)(24) in Eq. (14), as follows; The equivalence of Eq. (17), which is used to satisfy equilibrium between the crack and the continuum for the damage-healing case, is From which the following relationship may be derived; Using Eq. (27) in Eq. (25) gives; from which the following element stiffness relationship may be derived; in which and where w h is defined at the element level. In the following sections the specific healing agent transport, curing and crackplane constitutive models that are employed in this study are described. The authors (25)  emphasise that the element is general and could be readily combined with other forms of healing model.

Healing agent transport and curing model
The transport of the healing agent through the crack was calculated using the modified Lucas-Washburn equation of Gardner et al. [72][73][74] (see also Selvarajoo et al. [75]), that reads: where z denotes the cumulative flow length in the crack, the superior dot denotes a time derivative, P c is the capillary pressure, P app is the applied pressure, ρ is the density, g is the acceleration due to gravity, φ is the inclination of the crack, β s and β m are factors to allow for stick-slip of the meniscus and frictional dissipation at the meniscus respectively and η is the viscous resistance given by: where β w is a factor to allow for wall slip, k =w 2 r,C 12 is the permeability and µ is the dynamic viscosity of the healing agent.
The capillary pressure is given by the Young-Laplace equation: where γ denotes the surface tension and θ d is the dynamic contact angle. The dynamic contact angle is a function of the meniscus velocity and is described using the relationship of Jiang et al. [76]: where c 1 and c 2 are constants, and θ s denotes the static contact angle.
Cyanoacrylate cures through a polymerisation reaction in the presence moisture [75,77]. The moisture required for the reaction is transported into the glue from the substrate or the surrounding air. This leads to a reaction front which propagates into the glue, which, experimental evidence shows is diffuse in nature [75,78]. This can be described by [19]: where x denotes the position measured from the crack face, t is time, z c is a wall factor, z c1 is a diffusion coefficient and z(t) denotes the position of the reaction front.
The propagation of the reaction front is given by [19,75]: where z c0 is the critical reaction front depth and τ is the characteristic time.
The degree of mechanical healing found in Eq. (21) in the absence of re-damage-is given as the degree of cure at the centre of the crack:

Crack-plane constitutive model
In the present work, different constitutive models -that are coupled through the equilibrium of stresses and the crack initiation criterion [58] are employed for the continuum and the discrete cracks. The continuum is assumed to be linear elastic, whilst the discontinuity is assumed to follow the elastic-damage-healing cohesive crack constitutive relationship given in Eq. (21).

Damage-only model
The damage variable ( ω ) is a function of the damage evolution parameter ( ζ ), as follows: where u t = f t K where f t is the tensile strength of the material, K = E/h cp is the crackplane stiffness, E is Young's modulus of the continuum material (also the crack-plane band of elastic material), h cp is the assumed thickness of the crack-plane (taken as 10 mm u.n.o.), u m is the relative-displacement at the effective end of the softening curve (taken as 0.2 mm u.n.o.) and c 1 = 5 is a softening constant.
The evolution of ζ is governed by the following damage function:

Damage-healing model
The damage component of the model is unchanged from that described by Eqs. (37) and (38). Healing depends on the flow and curing of the healing agent, as explained in the previous section.
To allow for re-damage and re-healing, the evolution of the healing variable, found in the elastic-damage-healing cohesive crack constitutive relationship given in Eq. (21), is to be based on the following; in which the relative area (a) of the crack exposed to healing agent is given by; (35) in which Δa c is the incremental area of virgin filled crack, ∆a redam is the incremental area of re-damaged material and ∆a rec is the incremental area of re-filled cracks.
An important aspect of the damage-healing model is that an increment of healing causes no change in mechanical energy when it takes place in a static crack. This is equivalent to saying that healing agent cures in a stress-free state. Applying this principle to Eq. (21) for a healing update (Δh), with ũ remaining constant, leads to the following expression: from which the update ( ũ h ) may be derived to be:

Example problems
In this section, example problems are presented to demonstrate the performance of the model. Crack continuity was ensured using the algorithm of Alfaiate et al. [49,79] whilst U-turns in the crack path were prevented using the approach of Cervera et al. [80]. The latter two examples concern self-healing concrete specimens with embedded vascular networks. Depending upon the configuration of the channels, such systems can require a three-dimensional approach, since the transport of the healing agent may not be uniform across the width of the specimen. However, for the present examples, the authors believe that two-dimensional idealisations are adequate. This is justified by the fact that in the direct tension example, experimental observations showed that the healing agent had spread uniformly across the width of the crack [81,82]. Whilst in the L-shaped specimen example, the embedded channels are hypothetical and it is assumed that there are sufficient channels for the healing agent flow front to be effectively uniform across the specimen. The model parameters employed for the flow model can be seen in the Appendix.

Convergence test
This example considers a hypothetical prismatic singly notched specimen that is loaded in tension, as illustrated in Fig. 3. In this test it was assumed that the rate of transport of the healing agent in the crack was much faster than the rate of loading such that the crack was filled instantaneously with healing agent. The finite element meshes used in the analysis can be seen in Fig. 4, whilst the model parameters are given in Table 1. The number of elements used ranged from 36 for Mesh1 to 1296 for Mesh4, the time step employed was 0.5 s. A convergence tolerance of 0.01% for the out of balance force and L 2 iterative displacement norms was employed in this example and an average of 5 Newton iterations were required per time step.
(40) a = a + a c − a redam + a rec The load response curves from the simulations are given in Fig. 5. It can be seen from the figure that there is a significant difference in the load-response when healing is considered, in addition to there being a greater difference in the responses predicted for each of the meshes employed. To analyse the mesh convergence, the L 2 -error norm of the response curves was employed: where P and P r are the vectors of loads found in the response curves for a given solution and the reference solution respectively and the integral is taken over the response curve.
In the absence of an analytical solution to the problem, the solutions obtained with Mesh4 for the damage only and damage-healing cases were taken as the reference solutions. The results of the convergence study can be seen in Fig. 6, where the error norm is plotted against the element size, h e . The figure shows that the characteristics of the convergence behaviour of the model are the same for the damage only and damage-healing case, but that the error is larger when healing is considered.

Direct tension test
This example considers a set of direct tension tests on doubly notched concrete specimens (illustrated in Fig. 7), which contained embedded healing channels that were filled with cyanoacrylate (CA), presented in Selvarajoo [81] and Selvarajoo et al. [82]. The tests involved loading the specimen until a macro crack had formed and then opened to a given crack mouth opening displacement (CMOD). The CA supply was then released whilst the crack was then held at the CMOD value for a selected period of time. At the end of this 'healing period' , the specimen was loaded to until the CMOD clip gauge reached 0.3 mm, after which the specimens were unloaded. The test series considered two different crack openings and fixed healing periods of 0 s, 60 s and 600 s (and 1200 s for the 0.2 mm opening). The test set up can be seen in Fig. 7. The finite element meshes used in the analyses can be seen in Fig. 8. The model parameters are given in Table 2, where the subscript h indicates that these parameters are for the healed material. The (43) (P r − P)(P r − P)ds (P r )(P r )ds It can be noted that in this example, following experimental observations, the degree of healing was limited to 0.85. This is due to the fact that in some areas of the crack, the healing agent was in constant flux and as such did not stabilise and cure [82].
The comparison of the numerical simulations with the experimental data is given in Fig. 9. It can be seen from the figure that the numerical predictions accurately capture the experimental behaviour for both the 0.1 mm and 0.2 mm cases. The comparison of   the predicted and experimentally observed crack pattern can be seen in Fig. 10. It can be seen from the figure that there is a generally good match, although the experimental crack path is more tortuous. In this example, cracks propagated from the notches to the centre of the specimen, where they coalesced into one crack that crossed the centre of the specimen. The crack was fully formed prior to the release of healing agent in both cases and, as a result of this, the healing had no effect on the crack pattern. The numerical response curves for the two meshes are coincident, indicating mesh convergence.

L-shaped specimen test
This next example concerns the L-shaped specimen presented in Winkler et al. [83]. The test considered both plain and reinforced concrete specimens and both the load displacement curves and crack patterns were recorded. In the present study, the plain specimens were considered. To investigate the effect of healing, two hypothetical embedded channels filled with CA and placed vertically at a distances of 175 mm and 115 mm from the right-hand-side of the specimen were considered. The location of these hypothetical  (2020) 7:32 channels can be seen in Fig. 11. The finite element meshes used in the analyses can be seen in Fig. 12, whilst the model parameters are given in Table 3. The meshes include two structured meshes and one unstructured mesh that contains distorted elements. The number of elements was 770, 1158 and 1976 for Mesh1, Mesh2 and Mesh3 respectively, whilst the time step size used was 10 s. A convergence tolerance of 0.01% for the out of balance force and L 2 iterative displacement norms was employed in this example and an average of 5 Newton iterations were required per time step. The comparison of the numerical simulations with the experimental data can be seen in Figs. 13 and 14, whilst contour plots of the degree of healing in the crack at four points in the test (indicated by markers a-d in Fig. 13) is given in Fig. 15. It may be seen from the figures that the numerical predictions for the damage only case are able to accurately capture the experimental behaviour in terms of both the load-displacement curve and the predicted crack patterns. It is interesting to note the healing has a significant effect on the load-displacement curve, but that the predicted crack pattern remains unchanged. A key feature of the post peak response for the healing case is the two distinct post healed peak loads, each of which follows the release of CA from one of the embedded channels. The contour plots show that the degree of healing is greater than 80% throughout the crack at the first healed peak load, and that by the end of the test almost all of the healed material has been completely redamaged. The numerical response curves for the three meshes are in very close agreement, indicating mesh convergence. In addition to this, the predicted crack patterns are also in good agreement.

Conclusions
A novel specialised element for simulating self-healing in quasi-brittle materials was presented. The following conclusions can be drawn from this work: • The introduction of the healing variables at the element level ensures consistency with the internal degrees of freedom associated with the crack. • The equilibrium condition over the crack plane can be satisfied in a weak sense using the virtual work associated with the internal degrees of freedom. • The proposed element allows for different healing mechanisms and is readily coupled to a transport model that could be used to calculate the amount of available healing agent at a damage site. • The proposed element allows for overlapping cycles of damage and healing when combined with a suitable crack-plane model.     • The new finite element model is able to accurately capture the damage-healing behaviour of a vascular self-healing cementitious material system in terms of both the mechanical response and the observed crack patterns, as demonstrated with the validation examples.