A crack-opening-dependent numerical model for self-healing cementitious materials

A new damage-healing model for self-healing cementitious materials is described. The model is formulated using results from a discrete ligament model and guided by the findings of a linked experimental study. Healing is simulated using the interaction of curing fronts propagating from opposing crack faces within a body of healing-agent. This approach accounts for the dependency of the healing response on the crack opening displacement (COD) and its rate. The new damage-healing cohesive-zone model is applied to an element with an embedded strong-discontinuity within a coupled finite-element code, which simulates healing-agent transport and mechanical behaviour. The model is validated using data from tests with different CODs and COD rates. The validations show that the coupled model represents the mechanical and flow behaviour of an autonomic self-healing system with good accuracy for a range of cracking configurations and load paths.


Introduction
Biomimetic construction materials with the ability to self-heal are being developed so that future structures do not suffer from the cracking-related durability problems that have afflicted a significant proportion of our existing infrastructure (Gardner et al., 2018).The potential of these materials 'to deliver major change, most noticeably in the built environment' has been highlighted in a recent Royal Society report (The Royal Society, 2021).Numerical models can play an important role in the development of biomimetic materials because, once calibrated, they allow a larger parameter space to be investigated than is practicable with experimentation and they permit materials to be tailored to specific applications.In addition, the ability to simulate the cracking and healing behaviour of self-healing cementitious materials is important when designing structures formed from these materials (Jefferson et al., 2018).Such simulations should address the size and disposition of cracks, healing time and healing efficiency, changes in permeability, and load-displacement responses, all of which are important considerations in the design of this type of cementitious structure.
This work focusses on the simulation of autonomic self-healing cementitious materials (SHCMs) with embedded healing-agents (Van Tittelboom and De Belie, 2013;De Belie et al., 2018;Xue et al., 2019), but the processes studied and the new relationships developed in this work are relevant to a wide range of biomimetic materials and healing systems.
Early self-healing (SH) material models considered only mechanical behaviour (Barbero et al., 2005;Schimmel and Remmers, 2006;Granger et al., 2007) and nearly all of these, as well as many more recent contributions, were based on continuum damage mechanics (CDM) (Lemaitre and Desmorat, 2005).CDM is a natural choice for SH material models because the damage variables used in these formulations provide convenient measures of the material available for healing.The extended theory has been named 'continuum damage healing mechanics' (CDHM) and models based on this framework normally use a healing variable (or tensor) as a multiplier on the damaged portion of material (Darabi et al., 2012;Pan et al., 2018).CDHM theories have also been applied in micromechanical formulations (Davies and Jefferson, 2017), which are particularly well-suited to the simulation of SHCMs containing microcapsules (Han et al., 2021).
Thermodynamic principles underpin most damage formulations, with constitutive relationships normally derived from thermodynamic free-energy potentials such that models always predict non-negative energy dissipation (Lemaitre and Desmorat, 2005).The same basic principles are used in CDHM except that thermodynamic potentials are expanded to incorporate one or more healing variable(s) (Barbero et al., 2005;Voyiadjis et al., 2012;Alsheghri and Abu Al-Rub, 2015).An important thermodynamic aspect of the healing process is that, in contrast to damage, it does not dissipate mechanical energy since healing-agents normally cure within cracks in a nominally stress-free state.This implies that under isothermal constant-strain conditions, the recovery of stiffness due to healing causes no change in either the stress or mechanical energy state within the surrounding continuum.Furthermore, when healing occurs under non-zero strain conditions, a permanent 'healing' strain develops that is associated with the hardening of healing-agents in micro-cracks (Mergheim and Steinmann, 2013;Davies and Jefferson, 2017;Jefferson et al., 2018).
A range of approaches has been used for simulating time-dependent healing.Some models treat healing as a similar process to damage that develops according to a thermodynamic driving force (Barbero et al., 2005;Voyiadjis et al., 2011).A number of investigators have employed rate dependent healing approaches (Voyiadjis et al., 2011;Darabi et al., 2012), while others have developed time dependent exponential healing functions (Mergheim and Steinmann, 2013;Oucif et al., 2018).Wool and O'Connor (1981) introduced a recovery variable based on a convolution integral in their polymer healing model.Mergheim and Steinmann (2013) also derived a convolution integral from their exponential healing function and developed a convenient two-level recursive scheme for evaluating healing variables in a transient analysis.These different approaches reflect variations in the host material, healing systems, healing-agents and the loading scenarios considered.
CDHM is applicable to the formation and healing of micro-cracks and is well-suited to SH materials with embedded micro-capsules that contain healing-agents (White et al., 2001;De Belie et al., 2018).Alternative self-healing systems, with embedded hollow fibres (Dry, 1994;Van Tittelboom and De Belie, 2013), brittle tubes (Joseph et al., 2010;Van Belleghem et al., 2018) or interconnected vascular networks (Minnebo et al., 2017;De Nardi et al., 2020), are used for healing macrocracks, since larger volumes of healing-agent are required to fill such cracks.The simulation of these SH material systems requires an approach that represents macro-cracks and their healing.Most SH models for discrete cracks use cohesive crack idealisations based on damage-healing formulations (Caggiano et al., 2017), which have a similar form to their CDHM counterparts.These models have been incorporated in finite element codes using discrete elements (Zhou et al., 2017), elements with embedded strong discontinuities (Zhang and Zhuang, 2018;Freeman et al., 2020) and the extended finite element method (Gilabert et al., 2017a;Gilabert et al., 2017b).
In all of the self-healing materials discussed above, the healing process involves autonomic healing compounds being transported to crack sites and subsequently curing within micro-or macro-cracks.Relatively few attempts have been made to simulate the associated flow processes (Jefferson et al., 2018;Mauludin and Oucif, 2019), although Gilabert et al. (2017a) and Gilabert et al. (2017b) simulated the flow of a polyurethane-based healing-agent into a concrete crack using computational fluid dynamics and Gardner et al. (2014Gardner et al. ( , 2017) )  In this paper, the only damage process considered is cracking due to mechanical loading, thus the terms damage and cracking are used synonymously.
The primary purpose of this article is to present a new damagehealing model for simulating healing in both static cracks and in cracks with transient crack opening displacements (CODs).The paper describes how results from a one-dimensional ligament model were used to determine the governing relationships of a homogenised model that simulates the behaviour of a representative zone of material.The novelty of the model lies in the way that, (i) homogenised damage and healing variables evolve to represent interspersed damaged and healed material, (ii) healing is simulated using the interaction between diffuse curing fronts emanating from opposing crack faces with a COD rate dependent factor, and (iii) generalised curing front variables are accumulated and updated when re-damage and re-healing occur.The model simulates the strong dependency of healing on the COD and its rate, and places no restrictions on the timing or rate of cracking or healing.This is important because cracks are rarely inanimate during healing, and any movement of the crack boundaries affects the healing process.Furthermore, cracks in almost all applications vary greatly in their dimensions as well as the rate at which they open and close during healing.The degree of healing and rate of healing can change by orders of magnitude with the COD or COD rate.Thus, models that assume instantaneous healing or use a single healing rate function with no allowance for the COD or its rate can give very inaccurate results.In addition, models that assume that healing and cracking processes always occur over separate discrete time intervals are very restricted in the range of problems that they can simulate accurately (Jefferson et al., 2018).
The development of the new damage-healing cohesive zone model described herein formed part of a wider programme of work on the characterisation and modelling of SHCMs.This work included an experimental study on the transport and mechanical processes that govern the behaviour of SHCMs (Selvarajoo et al., 2020a;2020b); the development of a computational transport model (Freeman and Jefferson, 2020) (Appendix C) to describe the release, flow and curing of healing-agents in the cracked cementitious host material; and the derivation and implementation a specialised finite element with an embedded strong-discontinuity (Freeman et al., 2020).It is noted that the constitutive model presented in outline in the latter article was a considerably simplified version of the model presented in this paper that did not fully account for COD or COD rate effects.The coupled finite element model formed from the components described in (Freeman and Jefferson, 2020) and (Freeman et al., 2020), along with the new damage-healing model described herein, was used for the simulations reported in this article.Although our models were developed using data from a specific SHCM, the modelling approach should be applicable to a wide range of self-healing systems and materials.In the remainder of this paper, the autonomic system used in this study and the associated governing processes are described, along with other preliminary and background information (Section 2).The homogenised damage-healing model and its derivation from a set of onedimensional discrete ligament model simulations is explained (Section 3).An approach for computing the healing parameter from interacting curing fronts is presented (Section 4).An overall solution framework is described and details of the algorithm used to implement the damagehealing formulation in a cohesive zone model are given (Section 5).A set of examples is presented that evaluates the ability of the coupled model to represent the transport and mechanical behaviour of SHCM specimens for a range of transient loading conditions.The sensitivity of the predicted response to variations in the damage-healing model's material properties is also shown (Section 6).The main conclusions from the work are then given (Section 7).A proof of the thermodynamic validity of the formulation and details of the transport model are provided in the Appendices.

Autonomic healing system and experimental observations.
We studied the processes that govern the behaviour of an autonomic cementitious self-healing system in a series of experiments (Selvarajoo et al., 2020a(Selvarajoo et al., , 2020b)).The system studied comprised a vascular network embedded in concrete specimens with PC20 cyanoacrylate (CA) as the healing-agent.CA was chosen because it is a relatively fast acting agent that allowed us to study simultaneous cracking-healing processes in tests of modest duration (i.e.1-30 min).
The experiments undertaken to study the mechanical behaviour (Selvarajoo et al., 2020a) of the system included notched beam tests (SF and SO series, see Fig. 1a) and direct tension tests on doubly notched cube specimens (DT1 and DT2 series, Fig. 1b).The experimental procedures for these tests are outlined in Section 6.The SF and DT series were designed so that healing took place under fixed crack conditions.By contrast, the COD increased at a constant rate in the SO series tests so that cracking and healing occurred simultaneously.A summary of the tests considered in Section 6 is given in Table 1.
The experimental programme of tests used to study the transport and curing characteristics of the system (Selvarajoo et al., 2020b) considered the capillary flow behaviour of CA in static macro-cracks, the sorption of CA into concrete specimens through natural crack surfaces, the development and progress of CA curing fronts adjacent to a concrete substrate and the dynamic flow characteristics of CA in capillary channels.These transport and curing experiments, as well as those undertaken by others (Ferrara et al., 2018), showed that healing-agent is released from supply channels, flows into the crack and surrounding micro-cracked zone, and then cures gradually.The experiments showed that a curing front develops in the body of healing-agent adjacent to each crack wall, and then progresses towards the centre of the crack, with the curing front becoming more diffuse over time, as explained in Section 4.

Flow and curing processes
In Freeman and Jefferson (2020) we presented a transport and curing model for this healing system that was used for the simulations reported in Section 6.The transport component of the model is summarised in Appendix C. The curing model component simulates a polymerisation reaction which is initiated by hydroxide ions that are transported from the cementitious substrate into and through the body of CA via a diffusion process.The rate of progress of the polymerisation, or curing, front gradually slows as the cured material forms an increasingly impermeable barrier to the further migration of OH -ions (Tomlinson et al., 2006).We showed that the mean position of the curing front (z f (t c )) at curing time t c is well represented by the following equation: where z c0 is the critical curing depth and τ is the curing time parameter.
By solving an advection-diffusion equation for the degree of cure  (ϕ), we showed that the distribution of ϕ across a crack may be simulated with the following equation: ) ) (2 where x is defined in Fig. 2, z c1 is the curing front constant and z c2 is the wall factor.Fig. 2a shows a comparison between a predicted ϕ distribution obtained using (2) and measured values obtained from the curing front experiments described in Selvarajoo et al. (2020b).The parameters used in the comparison shown in Fig. 2a and further background to the curing front function is provided in Appendix C. A schematic diagram illustrating diffuse curing fronts progressing from opposing crack faces is presented in Fig. 2b.Throughout this work we have assumed that, within a filled crack, opposing curing fronts are symmetric.

Observed mechanical behaviour and applicability of damage mechanics
Fig. 3 shows characteristic responses from the beam and direct tension tests in which cracks had fixed CODs during the healing period (DT and SF series), and from beam tests that maintained a constant CMOD rate during healing (SO series).The former (Fig. 3a and b) show that the re-cracking response of healed specimens has the same characteristic softening behaviour as that of the virgin material, although the plots also show that the peak strength and fracture energy of the virgin and healed material can be markedly different from each other.It is well-known that this type of fracture response can be replicated with a cohesivezone damage model (Elices et al., 2002;De Borst et al., 2012).
The responses of the SO specimens were very different from those of the fixed crack tests, as illustrated in Fig. 3c.In these specimens, cracking and healing progressed simultaneously and the material recracked and re-healed continuously throughout the experiments.This type of behaviour is far more challenging to represent, particularly when the healing dependency on the COD and COD rate is taken into account.A new cohesive zone damage-healing model to represent this behaviour is described in Sections 3 and 4.
In the description of the model, parameters related to healed material are designated with the subscript h.The model distinguishes between the first and subsequent occurrences of damage and healing, the former being designated 'virgin' events (subscript v) and the latter with the prefix 're-' (subscript r).

Cohesive-zone crack-plane model
The cohesive-zone model uses the concept of a 'crack-plane', which is defined as the mid surface of a narrow band that contains a macrocrack and/or a number of micro-cracks (Jefferson and Mihai, 2015).The crack-plane model relates the crack-plane tractions (σ) to the relative-displacements across the crack plane (ũ), with the three components of both the σ and ũ vectors coinciding with the directions r to r 3 (r 1 being normal to the crack-plane mid-surface, and r 2 & r orthogonal in-plane directions).Prior to cracking, the relationship between σ and ũ is defined by the following elastic constitutive relationship: where ke is an elastic crack-plane constitutive matrix that depends on the elastic moduli of the cementitious material (E c ) and the assumed width of the cracking zone (h fpz ) (Appendix A). ũe is the elastic relativedisplacement vector.The inelastic component of ũ is obtained from The components of σ are also equal to the normal and shear components of the transformed total stress tensor (i.e. This type of cohesive zone model may be applied to interface finite elements as well as to embedded discontinuities within finite elements (De Maio et al., 2021).Such models may also be applied to crack planes within finite elements using a smeared model, with the equivalent relative-displacements being computed from the transformed strains multiplied by an element characteristic length (Jefferson and Mihai, 2015).

Damage-healing model
Our approach to developing a damage-healing model for use in our coupled finite element program was first to select a viable theoretical framework based on our experimental observations (see Section 2.3) and then to use a discrete ligament model to explore damage-healing characteristics over a representative area of crack.The results from a series of ligament model simulations were used to guide the development of a homogenised crack-plane damage-healing model.
The ligament and homogenised damage-healing models discussed in this section are independent of the COD and assume that cracks fill instantaneously with healing agent.Here, we concentrate on understanding the damage, re-damage, healing and re-healing processes over a representative area.COD and COD rate effects are addressed using the overlapping curing front model described in Section 4, along with the method used for relating the curing front variables to the homogenised healing variables.

Damage and healing discrete ligament model
Using the ideas from statistical damage theories (Krajcinovic, 1996), we assume that virgin and healed strengths of nominal ligaments (f v (ξ) and f h (ξ) respectively) of a representative cracking-healing area (A r ) may be defined in terms of a relative area variable ξ ∈ [0, 1].These where ke 1 and ke h1 are elastic normal stiffnesses (Appendix A) for virgin and healed material respectively, and A r is nominally assumed to be ~(5d agg ) 2 , where d agg is the coarse aggregate particle size (see Gitman et al. (2007) for a discussion on representative volumes).
The proportion of damaged material is given by the damage variable experienced by an increment of material over the time interval t 0 to t.Furthermore, the value of the relative area variable ξ for ζ = ζ ξ (ξ) is equal to the damage variable ω(ζ).This damage variable can be described as a cumulative probability function P b (ζ) which gives the proportion of material that is likely to be damaged when the relativedisplacement parameter is equal to ζ (Krajcinovic, 1996).The corresponding damage variable for healed material is denoted ω h (ζ h ) and functions for ζ ef , ω and ω h are given in Appendix A.
The link between ξ and the effective relative-displacement parameters is illustrated in Fig. 4 for a particular damage-healing scenario that involves damage, healing, re-damage and re-healing.
When an open fractured ligament heals, healing-agent fills and cures in the open space between the crack faces.This introduces a healing relative-displacement ũhξ (ξ, t) that is equal to ũ at the time of healing.This means that that there is no change in stress in a ligament due to healing alone, which is important from a thermodynamic point of view (Appendix B).Taking this healing relative-displacement into account, the traction vector (σ hξ (ξ, t)) on a healed material increment may be written as follows: where kh is the elastic crack-plane constitutive matrix for the healed material.
Considering infinitesimal ligaments (dξ) of undamaged and damagehealed material, the total traction vector by may be obtained by integrating across the representative area, as follows: ) ) dξ  where H ω and H h are Heaviside functions indicating damage and healing respectively, with H ω = 0 and H h = 0 indicating undamaged and unhealed ligaments respectively, and H ω = 1 and H h = 1 denoting damaged and healed ligaments respectively; noting that an infinitesimal ligament (dξ) is damaged or re-damaged when ) ) exceeds ζ ξ (ξ) or ζ ξh (ξ) respectively.Furthermore, H h changing from 1 to 0 signifies a re-damage event.
If ũ changes with time during the healing process, the stress and crack-healing displacement will vary across A r .To explore this situation, the traction history across A r has been evaluated over a time interval t 0 to t n at n t discrete times (t i ) using the following discrete approximation of (6): where n l is the number of discrete ligaments, Δξ = 1/n l , subscript j denotes the value of a variable at ξ = ξ j and Δt = (t n − t 0 )/(n t − 1); the time subscript i has been dropped from the right-hand-side terms for clarity.
The scenario chosen to illustrate the features of the damage-healing process is a band of cementitious material with the geometric and material properties given in Table A1 (Appendix A).A schematic of the system with 10 ligaments is shown in Fig. 5a.A normal displacement (ũ 1 ) is applied to all ligaments at a constant rate, from time 0 to t max , until it reaches the maximum value u lim .This creates a state of uniaxial tension in the band and so only scalar stress and relative-displacement measures are considered.It is assumed that there is a fixed time interval (Δt h ) between a ligament fracturing and re-healing (Table A1).Two cases are considered, which have the following maximum displacement and duration: case 1 u lim = 0.25 mm and t max = 240 s: and case (ii) u lim = 0.25 mm and t max = 60 s.
The overall traction-relative-displacement response of the system is given for both cases in Fig. 5b.The proportions of material that are damaged and healed at three selected times for case 1 are shown in Fig. 5c-e.This figure also shows the value of a parameter r ef (ξ, t) , which provides a measure of how close a material increment is to the re-damage surface.ref = 1 indicates material that is on the re-damage surface andr ef = 0 is associated with unstressed material.A convergence study showed that n l = n t = 4000 gave tractions results that were converged to within 1%, which was considered adequate for the present purposes.
The behavioural characteristics exhibited by the ligament model in this solution, in ξ space, include; I. the healing front (h v ) progresses from ξ = 0 towards ξ = 1, with damaged material becoming interspersed with healed material in the range 0 to h v as the applied displacement increases; II. it appears that discrete blocks of healed material (resembling waves) form within the range 0 to h v , and the number and size of these blocks tend to increase and decrease respectively as the damage-healing process proceeds; although the blocks would not necessarily manifest in this way in reality, since they artefacts of the increasing-strength representation in ξ space; III.ũhξ may vary in the range 0 to h v in an irregular manner; IV. r ef may vary from 0 to 1 within a block of healed material, implying a considerable variation in the propensity of different ligaments to re-damage with further displacement.

Damage and healing homogenised model
The above characteristics, along with some experimental observa-tions and theoretical considerations, led to a set of assumptions on which the homogenised model is based.These assumptions, the relevant characteristics determined from the discrete model (I-IV) and associated model components are now considered.
(i) Assumption: a homogenised form of Eq. ( 6) exists that relates the averaged traction vector to the relative-displacement vector, and which depends on the scalar damage variable (ω) and a healing variable (h ∈ [0, ω]), as well as a homogenised healing relativedisplacement vector ũh , as follows: (ii) Characteristics I and II from the discrete model solution (Section 3.1) suggest that the re-damaged material should be expressed as a proportion of h v and that this may be related to the homogenised healing variable as follows: where ζ h is the homogenised effective healed relative-displacement.
Remark: The total amount of material that is re-damaged at any one time can grow and diminish, which implies that ω h is not a conventional irreversible damage variable since it reduces when re-healing occurs.
(iii) Assumption: a time-dependent healing function exists that gives the proportion of material that is healed at time t in the absence of any re-damage.The following exponential healing function, used by other investigators (Mergheim and Steinmann, 2013), is assumed for the purpose of comparing the responses computed from the homogenised and ligament models: where 〈x〉 denotes Macauley brackets such that 〈x〉 = x if x ≥ 0 and = 0 if x < 0, t co is the time when curing begins and τ h is the healing time parameter.
Remarks: a linear homogenised healing model may be more consistent with the fixed healing period used in the ligament model; however, the linear function is far less convenient than (10) because of the need to track the end of the healing period after each occurrence of re-damage.Also, it was found that the response from the linear model in the range t = 0 to τ h (before the tracking is required) was very similar to that obtained with (10) for the parameter ranges considered.A more accurate approach for computing healing, which accounts for COD and COD rate dependency, is presented in Section 4.
(iv) The healing variable may be derived from the following convolution integral (Mergheim and Steinmann, 2013): where a(s) represents the relative area of curing healing-agent within the crack at time s.
Remarks: h and a change when re-damage occurs.(11) also applies to h v but this variable does not reduce with re-damage.
(v) Characteristic III (Section 3.1) relates to the variation of ũhξ (ξ, t) in the ξ range (0 to h v ).A homogenised healing displacement vector (ũ h ) may be derived from (6), for a particular time t, as follows: This expression is not convenient for the homogenised model since ũhξ (ξ, t) is not available.We therefore use the condition that a positive healing rate and zero damage rate should result in a zero stress rate.This condition ensures that no spurious strain energy is created when healing occurs (Appendix B) and may be expressed by the following rate equation: where the superior dot denotes a time derivative.
An equation for the ũh rate may be derived from ( 13) by neglecting second order terms, as follows: and in the absence of re-damage, ũh = ∫ t t0 uh dt.Remark: an equivalent expression may be derived for ũh by considering the contribution of each infinitesimal healing increment to its value (i.e.d ũh = ( ḣũ(t)/h ) dt).This results in a convolution integral for ũh , which is used for proving the thermodynamic validity of the model (Appendix B).
(vi) Characteristic IV (Section 3.1) relates to the variation of r ef across ξ and the implied variation of ζ hξ .As with ũh , it is inconvenient to derive ζ h from an integral expression involving an approximation of ζ hξ across ξ.A more convenient, and consistent, approach is to derive a homogenised value for ζ h using the following condition: (vii) Assumption: healing-agent is supplied instantaneously to damaged or re-damaged material increments, which implies that all of the damaged material is undergoing curing i.e. a = ω.
Remark: this assumption is not made in the coupled model described in Section 4, in which crack filling is simulated explicitly.
(viii) Assumption: h reduces when re-damage occurs at the rate ḣ = − ωh h v .
In order to compare the response of the homogenised model ( 8) with that of the discrete model ( 7), an algorithm is needed to integrate (8) over a specified relative-displacement path for the same time period considered for the ligament model (Fig. 5).The algorithm employs the following standard recursive relationship for updating the healing variables, which is based on the convolution integral (11) (Simo and Hughes, 1998;Mergheim and Steinmann, 2013): where the value of a variable at time t i is denoted with subscript i.
It is assumed that, within a time increment, damage and healing may be computed sequentially, provided that the time-step is converged.We assume that a response is time-step converged when the response at any time does not change by more than 0.1% when the time-step is halved.
In the damage sub-step, ω and ω h may increase, whilst h v remains constant and h reduces in accordance with (9), with the adjusted value being denoted h p .In the healing sub-step, h and h v are updated, and ζ h and the dependent ω h are computed using ( 16) and ( 17).
The algorithm is presented in algorithm Box 1.
Algorithm 1. Homogenised crack-plane damage-healing model algorithm for comparison with ligament model

Comparison between discrete and homogenised models
A comparison between the predicted responses from the ligament and homogenised models for cases 1 and 2 is shown in Fig. 5b.In both cases, the response of the homogenised model lies within the response envelope of the discrete model, which provides an indication that the assumptions used in to derive the homogenised model are reasonable.Matching the response of the ligament model is not an end within itself.The final arbiter of the accuracy of the model lies in its ability to reproduce a good range of experimental responses.Before these are considered in Section 6, the issue of the healing dependency on the COD and its rate is addressed.

The computation of healing from interacting curing fronts
As illustrated in Fig. 2, the degree of cure in an open crack varies across the body of healing-agent.Measurable mechanical healing commences when a solid phase of material first percolates across the body of healing-agent.This occurs after the symmetrical curing fronts from opposing crack surfaces overlap, with the degree of cure at the centre of the crack being 2ϕ m , where ϕ m = ϕ(w c /2, z f ).If the degree of healing were determined using a classical Evans-Avrami phase-change model (Wool and O'Connor, 1981), h f = 1 − e − 2ϕ m ; however, this function did not accurately represent the behaviour of our direct tension tests, as illustrated in Section 6.Rather, we found that once a threshold had been reached, the degree of healing is well-approximated using the assumption that the degree of healing is equal to ϕ m .Setting x = 0.5w c gives the following expression for the degree of healing in a stationary crack (h f stat ): Eq. ( 18) represents the degree of healing when the COD rate is zero throughout the healing period; however, when the COD changes, two important changes are required to capture COD rate effects.The first is to introduce a term (Θ ( ẇc ) that accounts for the disruption in the propagation of the polymerisation reaction caused by significant movement in the body of healing-agent.The second is to replace z f by a generalised curing front variable z that is accumulated in the manner described below.
By considering the ẇc rates in all of the tests (Selvarajoo et al., 2020a), we found that Θ ( ẇc was unity for nominally static cracks (i.e. ẇc was less than the threshold rate (w rt )) but approximately half this value when the rate exceeded the threshold.In the absence of detailed data around the threshold, we propose the following function that provides a smooth transition between the nominally static and dynamic crack rate domains: where the static and dynamic factors are f stat = 1 and f dyn = 0.45 respectively, w rt = 10 − 3 ζ m /τ h mm/s andw rnom = w rt /10.Eq. ( 18) is then replaced by the function that gives an equivalent degree healing for a given value of z: The generalised front variable is described into terms of the following convolution integral: where a(t c ) = A c (t c )/A cg is the relative area of a crack at a crack integration point that is undergoing curing, with A c being the actual area undergoing curing and A cg being the area associated with a crack numerical integration point or, in a general context, a representative crack area.a(t c ) depends on the curing time, the level of damage (ω) and the proportion of area filled with healing-agent (r fill ).Defining the relative fill area (a fill = ωr fill ) and making the assumption that an incremental component of healing-agent starts to cure immediately it is exposed to a new surface of crack wall, a accumulates as follows: As with (11), the convolution integral (21) may be solved over a finite time interval using the following recursive expression: The virgin generalised curing front variable z v , given in Eq. ( 24), is also required and this is used to compute the virgin healing front variable (h v ).
An increment of healing (Δh) is then computed using (20), as follows: and, where h i now represents the accumulated degree of healing at time t i .( 20) to ( 26) allow healing to be simulated for a crack with a timevarying fill-area, but further consideration is required to allow for recracking and re-healing.When re-cracking occurs, the accumulated curing front and healing variables must be reduced due to the loss of area (a red ).This is achieved by first reducing h by a red to give h r , as shown in ( 27), and then equating (20) to h r to find the consistent value of z (denoted z r ), as shown in ( 28), in which c h = atanh(2 hr Θ( ẇc) − 1).
This results in ( 23) being replaced with Eq. ( 29), as follows: Furthermore, it is necessary to account for the healed material present in the crack when re-damage and re-healing occur.This is accomplished by introducing the re-crack opening displacement (w r = w c − w h ) into the healing update (25), as calculated from the current total crack and healed crack displacements.
The accumulated virgin healing variable (h v ) and the current healing variable are then computed using the following expressions: where Δh v and Δh rh are the virgin and re-healing increments respectively.These steps form the basis of the damage-healing algorithm presented in Section 5.2.
The consistent set of damage, curing front and healing variables described in Sections 3 and 4 are summarised in Table 2.

Overall solution framework
We use the finite element method to solve the nonlinear coupled transport solid-mechanics problem using a staggered solution method.The governing equations for discrete and continuum healing-agent flow are given in Appendix C, and we use a conventional virtual work approach to derive the governing solid-mechanics system of equations.The discretised form of these equations is derived in the standard manner (Zienkiewicz et al., 2014) and may be expressed by the following matrix equation for the nonlinear system: where u g is the global nodal displacement vector, K g ( u g ) is the global nonlinear stiffness matrix and F g is the global vector of nodal loads.
) is assembled from the individual element stiffness matrices, which take either the standard form (K e ), for linear uncracked elements, or the special form with an embedded strong-discontinuity (K sd ) for cracked elements (Freeman et al., 2020).The nonlinear system of equations is solved using a standard Newton incremental iterative scheme (De Borst et al., 2012).This involves the nodal forces being applied in a set of time steps and then, at each time step, the out of balance force vector (Ψ) (33) being reduced in a series of iterations to a specified tolerance as measured by an L2 norm of Ψ.
where F λ (t) is the global applied force vector at time t, ∑ ASS denotes the element assembly process and B is the strain-displacement matrix.
The solution scheme for the coupled problem is given in Algorithm box 2, in which a sequential coupling procedure is employed that utilises sub-stepping for the transport problem (Freeman and Jefferson, 2020).
Set out-of-balance force vector to current load increment.for it = 1 to nt Enter iteration loop. δug= and χ e & χ m denote the convected coordinate of the healing agent meniscus in a discrete crack within an element and with reference to the supply channel respectively.

Healing update algorithm
The algorithm for updating the healing variables for the damagehealing model described in Sections 3 and 4 is presented in Algorithm Box 3.

Examples and the influence of healing material parameters
We assess the predictive performance of the coupled model using test data from our own experiments (examples 1 and 2) (Selvarajoo et al., 2020a), as well as data from a well-known concrete fracture test (example 3).The latter (Nooru-Mohamed, 1992) is re-imagined as a healing test in order to explore the healing response in a specimen with a more complex crack pattern than those that occur in examples 1 and 2.
The experiments considered in examples 1 and 2 all used the same healing-agent (PC20) and design concrete mix, although it was found that there were natural variations in the concrete properties between batches.The mechanical material properties are given in Table 3 and the flow properties, which were kept the same for all of the analyses, are given in Appendix C. The results reported in example 1 were obtained using the crackplane model alone and this example concentrates on basic damagehealing behaviour for different CODs and healing periods.The example also considers the sensitivity of the predicted responses to variations in the healing material parameters.The analyses reported in examples 2 and 3 were undertaken using the coupled finite element program developed for this work.All of the finite element analyses were two-dimensional and employed 4-noded bilinear isoparametric elements.The staggered incremental iterative Newton solution (Algorithm 2) was used to solve the nonlinear coupled transport-mechanical equations and a convergence tolerance of 10 -5 was used for L2 norms of the mass balance and residual force vectors.The number of time steps used for each analysis is given with each example and in every case a check was undertaken to ensure that the solutions were converged with respect to the time step size, using a tolerance 1% (based on changes in CMOD, reactions and meniscus flow position).

Example 1. Direct tension tests with different crack openings and healing periods
The direct tension (DT) experiments considered in this example used the doubly notched concrete cuboid specimens shown in Fig. 1b.In each test, the specimens were loaded until a crack opened to a pre-selected crack mouth opening displacement (CMOD) value (i.e.0.1 mm for DT1 and 0.2 mm for DT2).The test was then paused for a specified 'healing period' (t f ) whilst healing-agent was fed into the crack via a set of supply tubes and healing took place.Once the healing period had elapsed, the healing-agent supply was stopped and loading recommenced until a designated CMOD value had been attained.Here, we consider DT1 specimens with healing periods of 60 s and 600 s, and DT2 specimens with healing periods of 60 s and 1200 s.
In the experiments it was observed that the healing-agent did not cure over the full crack area, and to account for this we have used an effective curing area of 70% and 60% of the unnotched ligament for the DT1 and DT2 analyses respectively.For all but the final set (DT2_1200 secs), three tests were conducted with the same test parameters, in addition to a no-healing control test.Results from an example control test are included in the DT2_1200s graph.We have removed the results from those tests with known experimental problems (Selvarajoo et al., 2020a).A relevant experimental observation was that uncured healingagent was present in the crack after the supply was stopped for both of the 60 s healing period tests.This was simulated by assuming that the supply continued for 20 s after the supply tubes were clamped.This time corresponds approximately to when the outflow of healing-agent from the exit tubes was observed to cease.
In some of the experiments, the response in the first phase of loading exhibited an unloading-reloading loop (see Fig. 6b).This, along with other characteristics of the pre-healed softening response, are discussed in Selvarajoo et al. (2020a).The form of the cohesive-zone model presented in this paper does not account for this type of non-linear unloading-reloading response, but it would be possible to enhance the model to include such behaviour using an approach developed by the first author (Jefferson, 2002: Jefferson andMihai 2015).
A comparison between the experimental and numerical responses is given in Fig. 6 along with predicted healing time-histories.The average stress in the graphs is defined as the applied load divided by the  A.D. Jefferson and B.L. Freeman unnotched area (i.e.F L /5000 N/mm 2 ).Fig. 6 shows that the model is able to reproduce the major experimental trends with good accuracy; namely, (i) less healing occurs in a 0.2 mm crack than the 0.1 mm crack, and (ii) the tests with 60 s healing periods exhibit less healing than those with a t f of 600 s and 1200 s for DT1 and DT2 respectively.The numerical responses beyond the healing peak load are apparently less ductile than the experimental responses, particularly in the 60 s healingperiod cases.These differences are attributed to a flux in the uncured healing-agent beyond the end of the effective supply time that was not simulated in these crack-plane analyses.In addition, for the 60 s and 600 s DT1 cases, results from analyses undertaken using the Evans-Avrami (E-A) healing expression (Section 4) in place of ( 18) are presented.The results from these simulations, denoted 'Numerical E-A' in Fig. 6a and 6b, show that (18) provides a more accurate representation of the observed healing behaviour than the E-A equation for the partially cured 60 s case.As may be expected, there is less difference between the responses when the healing-agent is fully cured, as in the 600 s DT1 specimens.
We now explore the effect of varying the healing parameters (z c0 , τ and z c1 ) on the post-healing response.The DT1_60s case was selected because the healing-agent is partially cured at this time (60 s), which implies that the response should change with the healing parameters.Fig. 7 shows the effect on the average stress v CMOD response of halving and doubling the values of each of the parameters in turn, relative to those given in Table 3.The graphs show that the healing response is strongly dependent on τ and z c0 but that z c1 is of second order importance for the parameter range considered here.
The study presented here is contrasted with that we presented in Freeman et al. (2020).The example in the latter gave results for a single healing period and concentrated on showing the ability of the strong discontinuity element to represent crack-healing behaviour using a simplified damage-healing constitutive model.Here, we have explored the effects of healing time, as well as the sensitivity of the computed response to changes in the material parameters, for the cohesive zone model alone.

Example 2. Notched beam tests for a range of loading rates.
The second example is based on two series of tests on notched concrete beam specimens with embedded channels that were loaded in three-point bending, as illustrated in Fig. 1a.The channels were connected to healing-agent supply tubes that had the facility to be pressurised.
The experiments considered in this example comprise a set of tests from the SF series and three sets of tests from the SO series.In the former (SF tests), each specimen was loaded until the CMOD reached 0.15 mm.The healing-agent supply was then activated and the specimen held at this CMOD for 120 s.Loading was then restarted at a CMOD rate of 0.001 mm/s until the CMOD had reached 0.4 mm, at which point the specimen was unloaded.In the latter (SO tests), the healing-agent supply was activated from the start of the experiment and the load applied so as to maintain a constant CMOD rate once the primary crack had formed.The SO tests considered here are those with CMOD rates of 0.0005, 0.001 and 0.002 mm/s.All of the SO tests were continued until the CMOD reached 0.3 mm after which the specimens were unloaded.A delivery pressure of 0.5 bar was used for the SO tests; no excess pressure was applied in the SF tests.
The response of the specimens was simulated with the coupled finite element program using two meshes (see Fig. 8).The number of timesteps used for the SO and SF analyses was between 40 and 80, and the point load was applied using prescribed displacements that achieved the required CMOD rate once the primary crack had formed.A lower prescribed displacement rate was used in the initial stages to allow the prepeak response to be captured.The computed meniscus rise response and selected degree of healing results are presented in Fig. 9, and the experimental and numerical load-CMOD responses are given in Fig. 10.
The computed mechanical response for the SF case matches the experimental response with good accuracy, as shown in Fig. 10a.The results from the computations with the two different meshes are sufficiently close for the coarser mesh results to be considered mesh converged, and therefore only mesh 1 was used for the SO analyses.
The SF computed meniscus rise response (Fig. 9a) is not compared directly with any experimental data but the numerical results do have the same characteristics as those reported in Gardner et al. (2014) for CA in tapering concrete openings of similar magnitude.For this test, the healing period covered the time interval 100 to 220 s and the meniscus reached a plateau at a height of 60 mm at approximately 112 s.This means that healing took place under static flow and displacement conditions from 112 s to 220 s.The results show that this period was sufficient for full healing to take place.This conclusion is supported by the results of the full SF experimental series which showed that the experimental response did not change appreciably when the healing period was lengthened beyond 2 min (Selvarajoo et al., 2020a).
The numerical predictions for the three SO specimens (Fig. 10b-d) also match the experimental results well and show that the model is able to reproduce the main rate effect in the healing response i.e. that the first healing peak becomes less pronounced and less abrupt as the CMOD rate increases.Fig. 9 shows that the degree of healing through the depth at different times indicates that there is greater healing in the narrower sections but that there are some fluctuations in h near the top of the areas filled with healing-agent resulting from the damage-healing history at these locations.The tendency of the SO responses to reach a load plateau, above that of the corresponding control test, suggests that damage and healing rates balance each other in the latter stages of these tests.
The meniscus rise response (Fig. 9a) is relatively rapid in all cases.The upward movement of the meniscus, after the first plateau, reflects the extension of the primary crack as loading continues and the redamage of previously healed material.The latter creates a flow path for healing-agent to reach newly cracked sections.

Example 3. Nooru-Mohamed test 4B reimagined as a healing test
We used experiment 4B (Ref 46-05) from the test series of Nooru-Mohamed (1992) as the basis for the final example.Nooru-Mohamed's tests are well-known and have has been used by a number of investigators to validate their numerical models, e.g.(Cervera and Chiumenti, 2006).We first considered the original experiment with our numerical model to demonstrate its ability to reproduce the experimentally observed behaviour.We then reimagined the experiment as a healing test by assuming that a healing-agent supply system had been introduced into the specimen.The testing arrangement and imagined healing system are given in Fig. 11a, and the finite element mesh, boundary conditions and load path used for the analysis are given in Fig. 11b.As shown, a shear load (F s ) of 10kN was first applied to the upper part of the specimen.Then, while this shear load was maintained at the same value, the upper plate was restrained and displaced in a vertical upward direction with the associated normal load (F n ) being recorded.In the healing tests, it was assumed that the healing-agent was supplied at no excess pressure and activated once the shear load reached 10kN.A displacement rate of 2.7 × 10 − 5 mm/s was assumed for the upper surface, which is based on that used in the original experiment and 350 load steps were used in the solution.
The experimental and no-healing numerical responses are compared in Fig. 12a along with a comparison of the experimental and numerical crack patterns in Fig. 12b.The model is able to reproduce the experimental behaviour with good accuracy, in terms of both the load-displacement response and the crack path.The contrasting response of the imagined self-healing specimen is shown in Fig. 12a.This shows similar behaviour to that seen in the SO-series notched beams tests, except that there are some minor undulations after the first post-healing peak.The computed degree of healing and the extent of crack filling are shown in Fig. 12c at the first healing peak.At this stage, healing-agent had reached approximately 40% of each crack and some healing has occurred throughout the filled regions.
Although the computed healed response cannot be compared to experimental data, the example illustrates that the model works for multiple curved cracks and gives a credible solution.

Discussion
The examples show that the coupled model is able to reproduce the major trends observed in the experiments with good accuracy.In particular, the computed responses demonstrate that the new damagehealing model is able to capture the healing dependency on; (i) the COD, as evidenced by the contrasting numerical responses of specimens with different (nominally uniform) CODs (DT1 and DT2) and tapering cracks (SF); (ii) the healing time, which may be seen in the responses of the DT specimens with different healing periods; and (iii) the COD rate, as shown in the SO specimens tested at different displacement rates.The latter also illustrates the ability of the model to simulate the response of structural elements when damage and healing occur simultaneously over a sustained period of time.
The model does not reproduce the multiple healing peaks evident in the experimental SO-series responses.We believe that statistical variations in material strengths across the specimen and a more detailed representation of tortuous crack paths would be required to capture this behaviour.
The damage-healing model, which is the primary subject of this article, has been implemented in 3D although the overall simulations presented here are based on a two-dimensional idealisation.For full 3D simulations, a model that simulates healing-agent flow emanating from channels in curvilinear cracks is required.Such a model has been developed very recently by the authors and will be the subject of a forthcoming article; however, using the 3D flow model in the present examples makes negligible difference to the predicted response and the authors believe that including an account of this 3D flow model would detract from the main message of the paper.

Conclusions
Existing models built on continuum-damage-healing-mechanics principles are not able to simulate accurately simultaneous crackinghealing responses in systems that exhibit significant COD and CODrate effects.
A multi-ligament model of a representative tensile cracking zone, in which the ligament strengths are set according to an appropriate statistical function, can provide valuable insight into the interaction between damage and healing behaviour when healing takes place under continuous crack opening conditions.Simulations undertaken with such a ligament model, for a load case with a constant COD rate, show that discrete blocks of healed material may form over the representative crack area, with the number and size of these blocks tending to increase and decrease respectively as the damage-healing process proceeds.
A consistent set of homogenised damage-healing variables derived from the ligament model results are compatible with a cohesive zone formulation and capture the behavioural crack healing characteristics of a representative zone of SHCM.
The curing of an autonomic healing agent within a discrete crack progresses in the form of two diffuse curing fronts emanating from opposing crack faces.These nonlinear curing fronts may be assumed to be symmetric, with each front progressing in time according to an exponential function.The degree of crack healing may be related to the degree of cure in the overlapping zone between the opposing curing fronts.
A generalised curing front variable, expressed as a convolution integral and updated according to a two-level recursive scheme, is able to capture the curing characteristics within a crack during transient filling, healing, re-cracking and re-healing events.
The application of the curing-front approach in a damage-healing cohesive zone formulation produces a new model that simulates cracking and healing, and naturally accounts for the healing dependency on the COD and its rate.
The implementation of the new cohesive zone model in finite elements with embedded strong discontinuities provides an effective means of simulating the mechanical behaviour of structural components formed from SHCMs.
Overall, the behaviour of autonomic self-healing cementitious materials is governed by a set of interacting mechanical and transport processes.These may be simulated by a coupled finite element model that uses a Navier-Stokes discrete crack-flow model coupled to a mass balance equation for simulating matrix flow, and a combination of the strong-discontinuity element and new damage-healing cohesive model for mechanical behaviour.
A series of validation examples, based on recent experimental work undertaken by the authors' group, as well as data from the literature, show that the coupled model is able to represent the characteristic flow and mechanical behaviour of a vascular self-healing system in a cementitious structural element with good accuracy over a wide range of crack-healing conditions.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
In the present study, stability of the crack flow model is ensured using a continuous interior penalty, whilst a quasi-mass lumped storage matrix approach is employed for the matrix flow model.The transport model parameters used in the simulations are given in Table C2.Further background to the curing model Eqs.(1) and ( 2) are applicable to the curing of a body of CA when the thickness is relatively thin (<1mm) and the spread dimensions are at least an order of magnitude greater than the thickness.For larger thicknesses of healing-agent, thermal convection can become significant (Li et al., 2017) and a modification may be required to both equations (1) and (2).Since, the thickness of the cracks considered in this work do not exceed 0.4 mm, equations (1) and (2) are considered appropriate for the examples presented in Section 6.However, the thickness of the layer of CA used in our curing characterisation experiments (Selvarajoo et al., 2020b) was 6 mm and thus thermal convection is likely to have had an influence on the curing response in these tests.We have allowed for this in an approximate way by using a larger value of z c0 for the comparison shown in Fig. 2a than that used for the Section 6 examples.The material parameters used for the theoretical curve in Fig. 2 were as follows: τ = 60 s, z c0 = 0.3 mm, z c1 = 25 mm, z c2 = 0.0001 mm.
used a modified Lucas-Washburn equation to model the flow of healing-agents in discrete cracks.

Fig. 1 .
Fig. 1.Experimental arrangement for (a) SF and SO test series, and (b) DT test series.

Fig. 4 .
Fig. 4. Material states in terms the relative area variable (ξ), illustrated for the case in which the damage and re-damage evolution functions have the same material parameters.

Fig. 5 .
Fig. 5. (a) Schematic of ligament model.(b) Comparisons of ligament and homogenised model responses for cases 1 and 2. (c)-(e) Distribution of damaged and healed material from the case 1 ligament model solution at three separate displacements.
Note: max(a,b) = larger of a or b.

Algorithm 2 .
Staggered solution algorithm for transport and solid-mechanics equation systems, showing damage and healing phases t = t0; Fλ(t0) = 0 Initialise all cumulative variables.for i = 1 to n step Loop time steps.ti = ti− 1 + Δt; Read Fλ(ti)Set time variable and current loads.

Fig. 9 .
Fig. 9. (a) Meniscus rise responses from SO and SF analyses.(b) Meniscus position and degree of saturation at t = 343 s for SO_0.001.(c) Degree of healing in crack at t = 244 s and 343 s for SO_0.001.

Fig. 12 .
Fig. 12.(a) Experimental and numerical responses.(b) Final experimental and numerical crack patterns.(c) Central part (in height) of the finite element mesh showing the degree of healing at u n = 0.039 mm, F n = 9.66kN.

Table 1
Summary of experiments considered with the model.

Table 2
Summary of damage, healing and curing front variables.At damage sub-step:ζ h = max(ζ hprev ,

Table 3
Mechanical material parameters used for the analyses.