A Review on Cementitious Self-Healing and the Potential of Phase-Field Methods for Modeling Crack-Closing and Fracture Recovery

Improving the durability and sustainability of concrete structures has been driving the enormous number of research papers on self-healing mechanisms that have been published in the past decades. The vast developments of computer science significantly contributed to this and enhanced the various possibilities numerical simulations can offer to predict the entire service life, with emphasis on crack development and cementitious self-healing. The aim of this paper is to review the currently available literature on numerical methods for cementitious self-healing and fracture development using Phase-Field (PF) methods. The PF method is a computational method that has been frequently used for modeling and predicting the evolution of meso- and microstructural morphology of cementitious materials. It uses a set of conservative and non-conservative field variables to describe the phase evolutions. Unlike traditional sharp interface models, these field variables are continuous in the interfacial region, which is typical for PF methods. The present study first summarizes the various principles of self-healing mechanisms for cementitious materials, followed by the application of PF methods for simulating microscopic phase transformations. Then, a review on the various PF approaches for precipitation reaction and fracture mechanisms is reported, where the final section addresses potential key issues that may be considered in future developments of self-healing models. This also includes unified, combined and coupled multi-field models, which allow a comprehensive simulation of self-healing processes in cementitious materials.


Introduction
Concrete is characterized by its high compressive strength, a wide availability of its raw materials, and simple production methods, which is the main reason that it became the most commonly used construction material in the world [1,2]. However, its low tensile strength is the main reason that various types of cracks can occur in a concrete element that may adversely affect its service life [3]. While under internal, external, or environmental load, open or closed micro-and/or meso cracks may develop inside a concrete element that may successively result in a loss of structural integrity [4]. Open surface cracks may also allow water or hazardous substances to enter and thereby severely impairing its durability [5,6]. Therefore, improving the durability of concrete structures, asks for a limitation or reduction of the number of cracks where self-healing strategies could be solution.
When considering the number and type of experiments required to study the performance of self-healing concrete, it turns out that optimizing self-healing mechanisms through extensive experimental studies is a very demanding task. However, this task becomes more doable when employing numerical simulation models. However, most existing models do not incorporate physically/chemically driven boundary movements for an accurate simulation of solid-liquid interfaces. To overcome these difficulties, phase-field (PF) methods have been proposed as a powerful tool for handling moving interfaces caused by phase transitions [54][55][56]. In conventional numerical models for phase transformations and microstructural evolutions, interfaces are considered to be infinitely sharp and have to be schematized explicitly [57][58][59]. It leads to incompatibilities that makes calculations very complex and difficult to implement in a computer program. Contrarily, PF methods are based on thermodynamic principles and assume a diffuse interface, which makes them suitable for solving complex morphological evolutionary processes. The evolution of the "field", over time and space, is controlled by the nonlinear Cahn-Hilliard diffusion equation and its relaxation by the Allen-Cahn equation [60,61]. For concrete, a self-healing mechanism is physically almost similar to a dissolution and/or precipitation principle that evolves at the cracked surfaces. It makes a PF modeling approach very suitable for solving this type of moving interface problems at cracked surfaces, caused by phase transformations.
This article provides a review on existing models to simulate self-healing in cracked concrete, with emphasis on PF methods. After the introduction in Section 1, the currently available self-healing methods for concrete are reported in Section 2. In Section 3, the possibility of using PF methods for simulating self-healing in concrete is presented and discussed. Then, in Section 4, the basic equations of a PF method are presented. Next, in Sections 5 and 6 existing PF techniques for precipitation and fracture in concrete are reported, respectively. Finally, items that should be addressed in self-healing models along with future research priorities and a concluding discussion on the whole article is given in Section 7.

Self-Healing Mechanisms in Concrete
In general, self-healing processes in cement-based materials can be divided into two categories: (1) autogenous self-healing and (2) autonomous self-healing [9,62,63]. Autogenous self-healing involves only the original components of a concrete. These components may, due to their specific chemical compositions, promote crack healing under favorable environmental conditions, driven by chemical reactions or transitions [10,64,65]. However, autonomous self-healing processes can only take place with the help of healing additives, such as microcapsules that may contain healing agents like polymers or bacterial spores [14,66]. Autogenous healing mechanisms have a limited healing capacity, typically only being able to heal cracks of about 100-150 µm in width [67]. In contrast to this, autonomous mechanisms can easily heal cracks up to 300 µm and sometimes even more than 1 mm [67]. These self-healing mechanisms are described below in detail.

Autogenous Self-Healing
Autogenous self-healing has been extensively investigated in the last decades [9,10,65,68,69], mainly by using experimental techniques. Figure 1a shows three main categories: physical, chemical, and mechanical healing. The physical healing mechanism is the process where the crack surface inside a cement matrix absorbs water and causes volume expansion [70,71]. The chemical healing mechanism consists of two main reactions, namely, a further hydration of the still unhydrated cement clinker inside a concrete, generating additional Calcium Silicate Hydrates (C-S-Hs), and carbonation of the additionally formed portlandite [65,[72][73][74]. Finally, mechanical healing mechanisms refers to the filling of a crack with fine cement particles, which appear in a crack by water transport or diffusion [72]. The chemical mechanism is the primary and most promising healing method for hardened concrete at a young age [16]. Due to the relatively high content of unhydrated cement particles in these concretes, continuing hydration will still be possible and may result in a healing of cracks [18,64]. At later ages after crack initiation, the formation and growth of calcium carbonate crystals (CaCO 3 ) becomes the main healing mechanism [75]. Figure 1b,c shows the main healing products and their chemical components.

Mechanical causes
Fine particles: Broken of from fracture surface Originally in the water a b c Figure 1. The autogenous self-healing mechanisms, products, and their corresponding chemical composition. (a) Schematic representation of the mechanisms of autogenic self-healing. Reproduced with permission from the authors of [70]. Copyright 2013, Springer. (b) Morphology of healing products (GHP refers to the gel-like healing product and CHP refers to the crystal-like healing product). Reproduced with permission from the authors [76]. Copyright 2013, Elsevier. (c) Ratios of Ca/Si and Al/Si of healing products with time. Reproduced with permission from the authors of [76]. Copyright 2013, Elsevier.
To improve the effectiveness of autogenous crack repair, an improved self-healing method called Dissoluble Encapsulated Particles (DEP) has been proposed [11,12,77]. In this self-healing method a certain amount of cement in a concrete mixture remains unhydrated for a predefined period of time because of the pre-encapsulation of certain cement fractions which are covered with a thin membrane that can dissolve whenever it is affected by a crack (Figure 2). A crack in a cementitious surface may open the DEP membrane due to either (1) a dissolution mechanism caused by low pH-conditions, i.e., due to increased CO 2 ingress, or (2) by mechanical fracture. After this happened, the original unhydrated cement will be exposed to the local environmental temperature and humidity conditions causing the cement to react and finally close the crack [77].  [12]. (b) Initial state of microstructure by vol.−10% cement replacement by DEP [12]. (c) A high pH value will cause the DEP capsule to rupture, the healing agent will be released and a special hydration reaction with accompanying volume expansion will begin [77].

Self-Healing Based on Mineral Admixtures
Mineral admixtures are materials that are mixed in a concrete and react with water to form reaction products with an expanded volume to heal cracks developed in an already hardened concrete. With this healing mechanism [13,91,94,95], crack widths up to 120 µm can be repaired [67]. Depending on the type of mineral additives, three subcategories can be identified: (1) expansive additives, (2) geo-material based additives, and (3) chemical agents (crystalline additives) ( Figure 3). Expansive additives develop reaction products with an increased volume that can fill the cracks [96]. Commonly used are sulfoaluminate based expansive additives (C-S-A) [78]. The geo-material-based additives consist of silicon dioxide, sodium aluminum silicate hydroxide, and bentonite clay, which have the capacity to swell [79,97,98]. When this type of geo-material is exposed to water, its volume may increase 15-18 times its initial dry volume [79]. The most basic crystalline additive is tricalcium silicate (C 3 S), which is the main clinker component in cement and reacts with water to form calcium silicate hydrate C-S-H phases [26].

Self-Healing Based on Bacteria
A certain category of bacteria can be applied for healing cracks in concrete [28]. It results in a closed crack which is watertight and has a limited capacity to restore the mechanical strength of a concrete [80][81][82]. The maximum crack width that can be healed with this system are~150 µm [83], which is rather limited whenever compared with other healing systems. Figure 4 shows a schematic impression of a fractured concrete with microencapsulated bacterial spores and the results of previous experiments [80,81,99].
Bacteria provide an important reaction component in a self-healing mechanism, where they are enhancing the calcium carbonate CaCO 3 production, needed for crack closing [100]. During healing, the mechanism passes the following two sequential steps; (1) conversion of calcium lactate and (2) hydrolysis of urea through (ureolytic) bacterial metabolism. In the first mechanism, oxygen and water penetrate into the concrete interior through cracks where the bacteria are activated to convert calcium lactate into CaCO 3 crystals and CO 2 . Portlandite particles near the cracks will further react with CO 2 to produce more CaCO 3 which precipitates at the crack surfaces [81]. In the second mechanism, many components capable of producing organic urea (e.g., Bacillus cohnii, Sphaericus, Subtilis, Pasteurii, Megaterium, and Sporosarcina ureae) can act as a catalyst during the self-healing process [101]. As it undergoes demineralization, negatively charged bacterial cells take up components from the cell wall and then react to CaCO 3 precipitates [102].
The efficiency of the precipitates generated by bacterial induction is determined first by the available water content and moisture movement in the concrete matrix [103,104], and second by the concentration of calcium ions, the pH of the pore solution, the concentration of inorganic carbon and by the presence of nucleation sites [105,106]. The first three are available in the concrete matrix, while the last one is related to the type of bacteria used [82]. In addition, factors that affect the effectiveness of healing include (1) the type of carrier (direct [107], encapsulated [108] containers like clay and aggregates [109,110]) and (2) the concrete compatible chemical reactions taking place in producing CaCO 3 [86,111].

Self-Healing Based on Adhesive Agents
This method is based on injecting adhesives into a crack to induce manual healing [112,113]. The crack widths which can be healed with these systems vary from 50 µm up to 250-300 µm [67,114]. Adhesive agents can be divided into one-component and multicomponent systems. Commonly used one-component adhesive agents are polyurethane [115] and epoxy [116]. Multicomponent adhesives are methylmethacrylate [117] and ureaformaldehyde/epoxy [113]. Adhesive agents are encapsulated in spherical capsules [112], tubular-shaped capsules [117,118], and hollow fibers [119][120][121] that are mixed with fresh concrete ( Figure 5). When cracks occur, rupture of the encapsulation takes place, where the adhesive will be released into the crack by capillary action, initiating crack healing with time.
Recently, many problems in solid mechanics deal with the use of PF for describing fracture phenomena and to capture complex crack patterns [136][137][138][139][140]. Based on the present literature review, the following can be summarized.
• PF is an extremely powerful mathematical modeling scheme for accurately describing physical movements of phase boundaries. • PF was mainly employed for solving solidification dynamics, material phase changes/separations, growing phases driven by chemo-kinetics and transport phenomena, nucleation and coalescence processes between particles in micro-to-mesostructures. • PF has been successfully employed in fracture mechanics to capture the cracking response of brittle/ductile materials without the need for employing Discrete Crack Approaches (DCAs) and/or Smeared Crack Approaches (SCA).
Because of this, and as also supported by various state-of-the-art reports [55,56,[141][142][143][144], PF models can be employed for self-healing of brittle or plastic (ductile) materials in a fundamental and consistent way. It will combine the impact of two main phase changes that occur simultaneously in a self-healing mechanism, i.e., chemical reactions and fracture. Gradual changes from the fully-cracked (failure) to the uncracked configuration can be driven through the so-called Phase-Field order parameter (φ). It will provide a smooth transition of all relevant phenomena between the fully cracked configuration and the intact material phases: this strength and crack recoveries actually represent the self-healing process. The governing equations of the proposed unified model will be derived in the framework of thermodynamics concepts, in terms of kinematics and balance equations, dissipation inequality and constitutive laws. Particularly, the free energy will be considered as the sum of the contributions due to elasticity, reaction PF and fracture PF. The free energy of the system is described in a unified form over the entire phase transition region. In this regard, the advantage of the PF method over other competitive numerical methods is its enormous capability of capturing movements of interfaces, without the need for introducing any additional ad hoc technique, criteria and/or remeshing strategies, and also without any explicit tracking of the actual interface positions of these coupled processes. The governing equations of PF models for chemical/moisture reactions and fracture processes, associated with self-healing, as well as the coupling among them, can be formulated in a unified PF framework. The next sections report a review on the available formulations for a unified and coupled set of PF approaches for modeling reactions and fracture of self-healing mechanisms in concrete.

Main Equations of a Phase-Field Approach
The phase-field (PF) approach is a very powerful technique to simulate complex physical phenomena in multi-field environments. The main attributions of this approach are simplicity and generality. A popular PF application is a diffusion interface model that is frequently used to simulate phase transformation problems in materials research [145][146][147]. The classical PF method is formulated based on the theory of Ginzburg and Landau, elaborated in the 1950s [148]. Compared with the sharp interface model, the PF diffusion interface model has the important advantage that no boundary conditions are specified on the interface between the different domains ( Figure 6). A diffusive order parameter φ is a continuous function coordinate of time and space, which indicates each phase to convert between 0∼1 or −1∼1 within a thin translation layer [54,144]. Moreover, φ is controlled by a set of coupled partial differential equations that can be discretized and solved numerically by evolving the equations. Any phase transformation is driven by a reduction of the free energy of the system F, which can be described by a set of conserved c i and non-conserved φ i field variables. The domain of the model is the entire phase transition system. The free energy of the system consists of the energy contributions from the homogenous bulk phases F bulk and the diffuse interface region F int , according to [146] where f loc defines the local free energy density (including chemical, interfacial and elastic strain free energy density), while f int defines the diffusive interface energy density. From the computational point of view, monolithic or staggered algorithms can be computed to solve the problem unknowns, in which mechanical, chemical, interface, and phase-field variables are computed simultaneously or sequentially, respectively. For more details the interested reader is referred to the works in [136,149,150]. In those works, robust and efficient monolithic schemes were employed for the numerical implementation.

Evolution Equation
The generalized PF method is represented by the Ginzburg-Landau or Onsager kinetic equation combined with the well fitted Landau-or Redlich-Kister-type free energy density functionals, which are dependent on both conserved and non-conserved field variables [146]. The time-dependent evolution of the conserved field variables (chemical concentration) is defined using a modified Cahn-Hilliard equation [151], while the Allen-Cahn equation describes the transformations with non-conserved variables (e.g., crystal orientation, long-range order, crystal structure, and elastic strain) [152].
The Cahn-Hilliard equation is where c i is the conserved concentration field variable, M c is the kinetic coefficient of diffusion (associated mobility), t is the time and r is the spatial coordinate, ∇ is a vector of partial derivative operator, and δ denotes the variational derivation of the functional F.

The Allen-Cahn equation is
where φ i (r, t) are the i different structure field variables with i = 1, 2 ..., n, while L φ is the kinetic structure operators (order parameter mobility). Depending on the problem, L φ has different expressions [124,153,154].

Local Free Energy Function
The local free energy function is a key component in the PF model [155]. This function describes the free energy density of each bulk phase, whose coefficients are obtained from thermodynamic data [153]. The expression of the local free energy depends on the problem of interest. For example, a double-well form is often used for solidification [147,156]. When dealing with an electromigration problem, a double-obstacle potential is usually applied [55,157]. A crystalline energy function is used to describe an overlapped dislocation of an elastically anisotropic crystal [158][159][160]. When the problem is temperature-controlled, as in the melting and solidification processes of crystals, the local free energy function contains a temperature field [161,162]. In such cases, the phase-field is needed to be coupled with a temperature field [161][162][163][164][165]. Furthermore, a Landau-type polynomial potential can be applied for the treatment of a solid-state phase transformation [166][167][168][169][170][171]. Table 1 summarizes examples of the universal expressions, the graphs of the local free energies and existing phase-field applications. Table 1. Expressions, graphs, and applications of the local free energy. Double-well where A is the height of the potential energy between the two states at the minimum free energy.
dislocation system [158,159,202,203] spiral growth [160,204]  Potential with temperature field is the difference between the current temperature and the melting temperature; α is a positive constant.

Phase-Field Modeling of Precipitation Reaction Mechanisms
Self-healing of concrete can be numerically treated as a precipitation process of solutes at the solid-liquid crack interface [215,216], which is time-dependent and controlled by chemical reactions and diffusion [36,217]. When the rate of the chemical reactions at the interface is sufficiently high and there is no fluid flow, diffusion will be the only mechanism left for solute transport. The whole process is then a diffusion-controlled precipitation one [216]. However, when the chemical kinetics is slow enough, the precipitation process becomes chemically determined [218]. A review of existing models for self-healing that are based on chemical reactions show that these models are employing a reaction-diffusion process to describe the self-healing evolution [12,[31][32][33][34]. These models focus on two processes: (1) the diffusion mechanism where dissolved ions (e.g., calcium ions) are transferred from the concrete interior toward the surface of the crack, and (2) the precipitation of mineral ions reacting with, for example, carbon dioxide or carbonate ions to form calcium carbonate. They mostly consider how the chemical environment affects the formation of self-healing products and how to achieve agreement with experimental results [12,[32][33][34][35][36]41,45,219].
However, these models have several limitations. First, they only simulate chemical reactions in solution and do not explicitly account for the change of the initial solid phase boundary due to the dissolution of soluble minerals at the fracture surface. Reaction diffusion models only include precipitation reactions in solution and do not simulate the dissolution reactions of the solid phase with a solution. Second, these models only uniformly simulate the healing process at the crack and do not accurately simulate the change in micro-morphology of the crack. The change in crack morphology is directly influenced by the concentration of aqueous substances and precipitations inside the solution [111]. In return, the change in crack morphology does affects the local concentrations of aqueous substances and precipitations in the solution. This interaction between the two factors is not reflected by existing models.
A PF method can fill these gaps. Figure 7 shows schematically a potential application of a PF model for an autogenous self-healing mechanism. The solid-liquid phase distribution is described by an eigenfunction in the value range [0, 1]. The solid phase can be subdivided into an initial solid phase (φ 1 ) and a healing solid phase (φ 2 ), while φ 3 represents the solution phase. The solid-solid (I SS ) and solid-liquid (I SL ) interfaces are simulated continuously. In addition to the solution (D L ), diffusion constants are distinguished between the concrete (D S1 ) and the healing region (D S2 ) due to differences in the meso-and microstructures. Neumann boundary conditions (Zero composition flux) were applied at the top, bottom, left and right (the light gray part) boundary for the solute concentration c i and the order parameter φ i . The Dirichlet boundary condition (c 3 = 0.1 and φ 3 = 1) was applied at the right boundary (the blue part). The initial conditions are set based on the initial concentration in each phase. In this model, we chose to use the diffusion equation instead of Cahn-Hilliard equation because there is no phase separation. The Allen-Cahn equation is applied for solving the order parameter φ i .
This approach can accurately capture information about the alteration of the crack morphology due to solidification by the hydration reactions or the accumulation of precipitates [65,68,69]. With this, an overview of the PF approaches to the solute precipitation [180,220] and precipitation in binary alloys [179,221,222] is provided that are instructive for simulating self-healing mechanisms of concrete. The following models are presented in chronological order ( Table 2). Kinetic data of existing databases CALPHAD applied into the PF model.

Solute Precipitation
Solute precipitation is the process at which a solute changes from a liquid phase to a solid phase and precipitates outside its solution [230,231]. In fact, precipitates are mostly insoluble [232].

Xu-Meakin Model, 2008
Xu and Meakin [180,220] developed a PF model for studying the dynamics of liquid-solid interfaces due to precipitation and/or dissolution of phases, based on the Karma-Rappel model [154] for pure melt solidifications. Discontinuities in the solute concentration at the interface are explicitly considered. An additional term has been added to the solute diffusion equation to describe the discontinuity of the solute concentration gradient at the interface. In addition, a detailed asymptotic analysis was used to establish a connection between the sharp interface and the PF model by correlating the reaction rate parameter k with the microscopic PF parameters. This ensures that the PF model will converge to the corresponding sharp-interface limit. A modified solute diffusion equation is built up as follows, where the second additional term of the equation is corresponding to the discontinuity the solute concentration gradient at the interface. While the third additional term represents the net source or sink of the solute coming from the discontinuity in the solute concentration across the interface; D is the diffusion coefficient; A 1 and A 2 are two constants, which can be determined by the sharp-interface boundary conditions.

Noorden-Eck Model, 2011
Van Noorden and Eck [179] proposed a PF model for a precipitation and/or dissolution process. The model describes a single-phase free boundary problem with dynamic conditions at the moving boundary. The concentration on the precipitate side of the interface is specified, and the velocity normal to the interface is nonlinear dependent to the concentration on the other side of the interface. The evolution equation of φ and c is described according to where p(φ) is a double-well potential; f (c) is a rate function; k(φ) is an interpolation function; α, β, D, and ρ are physical parameters; and is the thickness of an interfacial layer. Redeker and Rohde [221,222] extended the Noorden-Eck model by incorporating curvature effects between two fluid phases to simulate precipitation in a porous medium. The model contains two immiscible fluids and one solid phase. Dissolved ions in one of the fluids can precipitate at the pore boundaries. Bringedal et al. [226] considered not only the diffusion of ions in the fluid phase, but also the effect of fluid flow on precipitation.

Metal Precipitation
Unlike solute precipitation, metal precipitation occurs in a supersaturated solid solution. Metals and metal oxides exist in the form of crystals. A crystal is a structure in which its atoms or molecules are arranged in an orderly fashion according to certain rules. A crystal is pure when all the components are just a single substance or a compound. If there is another substance involved that occupies the original atomic location and does not destroy the original structure, then this is a solid solution [233]. The original component is equivalent to a solvent and the foreign component is equivalent to a solute. As with a solution, when the solute in a solid solution is supersaturated in the solvent, it can no longer remain stable in the crystal structure and eventually precipitates [234].
The precipitate particles are generally metallic compounds, but may also be formed by aggregation of solute atoms in supersaturated solid solutions in a number of small solute-rich regions [235]. The precipitated particles act as barriers to dislocation movement, allowing significant increase in strength and hardness of most structural alloys of aluminum, magnesium, nickel, and titanium, as well as some steels and stainless steels [236]. The precipitation mechanisms of different binary and ternary alloys have been intensively studied by using PF models [182,237,238].

Wang-Chen Model, 1993
In the earlier study by Wang et al. [178], a PF model based on a microscopic kinetic model and elastic strain theory was developed to study the morphological evolution of the solid-state precipitation, controlled by transformation-induced elastic strain. The free energy of an inhomogeneous solid solution is given by the following equation, where φ(r, t) is the non-equilibrium single crystal sites of solute atoms, r is the crystal lattice site, W(r − r ) is the pairwise interaction energy of two atoms at the lattice site r and r , and k B is the Boltzmann's constant. The drawback of this model is that the matrix phase and the precipitates are iso-structurally treated. However, this assumption does not apply to the simulation of Al-Li alloy precipitation.

Rubin-Khachaturyan Model, 1999
Rubin and Khachaturyan [168] developed a 3D stochastic PF model for simulating the microstructural evolution of Ni-Al superalloys. This model considers the coherency strain in an elastic anisotropic system. The coarse grained stress-free free energy was expressed as where α ij and β ij (p) are the gradient coefficients, ∇ ic and ∇ jc denote the gradient terms of multi-composition profile c(r, t), ∇ iφ p and ∇ jφ p are the gradient terms of multi-component long-range order parameter φ(r, t), the specific free energy f (c, φ 1 , φ 2 , φ 3 ) is approximated by a polynomial, and the second integral term is the total strain energy functional based on the Fourier transform microelasticity method.

Chen-Ma Model, 2004
Chen et al. [177] designed a quantitative PF modeling scheme for multicomponent diffusion-controlled precipitate growth and dissolution in Ti-Al-V system in which the thermodynamic and kinetic data of existing databases CALPHAD was directly inserted into the PF model. The total Gibbs free energy is described as follows, where G m is the local molar Gibbs free energy; k i and k j are the gradient-energy coefficients for concentration and order parameter inhomogeneities, respectively; V m is molar volume. The temporal evolution of the composition is governed by Cahn-Hilliard diffusion equation on the basis of the phenomenological Fick-Onsager equations where M kj are chemical mobilities related to atomic mobilities.
Starting with the works of Bourdin et al. [258] and Miehe et al. [141], fracture processes were modeled explicitly by a PF approach. Due to its simplicity this methodology gained a wide interest and started to be used in the engineering community since 2010. From there on many scientist have worked in this field and developed PF approaches for finite elements methods (FEM), isogeometirc analysis (IGA), and recently also for the virtual element methods (VEM). The main driving force for these developments is the possibility to handle complex fracture phenomena within numerical methods in various dimensions. Thus, research on PF approaches is still actual and points in many different directions.
In this review article, the simulation of fracture processes in concrete is achieved by utilizing the continuum PF method, which is based on the regularization of sharp crack discontinuities. This avoids the use of complex discretization methods for crack discontinuities and can account for multi-branched cracks within a solid skeleton (e.g., hydrated cement paste, unhydrated clinker particles, and stones).
In particular due to the over-complicated geometry and content of concrete at multi-scales, in Figure 8 an example for PF modeling of water-induced failure mechanics in concrete microstructure is presented. In recent years, several brittle  and ductile [149, PF fracture formulations have been proposed in literature. These studies range from modeling 2D/3D small and large strain deformations, variational formulations, multi-scale/physics problems, mathematical analysis, different decompositions and discretization techniques with many applications in science and engineering.
All these examples demonstrate the potential of PF method for crack propagation.
The aforementioned PF approaches consider the fracture behavior of concrete, i.e., as a crack initiation and propagation. However, an important aspect in concrete is the treatment of the crack-closure effects. This response was firstly investigated in the works [325,326] for fatigue crack closure under cyclic tension. Thereby, the results indicate a fatigue crack, propagating under zero-to-tension loading may be partially or completely closed at zero load. A review of this physical phenomena can be seen in [327][328][329]. To the author's best knowledge, a PF approach for modeling crack closure is still an open issue. To this end, cohesive elements along the crack path will be coupled with the PF formulations to prevent overlapping of the crack faces. Another future direction is to use a contact scheme at the crack faces similar to the work developed by [330]. A further important aspect is the PF modeling of crack-closure induced by a self-healing mechanism (introduced in Section 3) in cementitious systems. These topics await investigation.

Fundamental Variational Formulations
In Griffith-type fracture formulations, the mechanical deformation denoted generally by "state" and the sharp crack surface Γ in a brittle elastic solid (e.g., cement paste) are determined by the incremental minimization problem developed by Francfort and Marigo [331] as where G c is the Griffith critical surface energy release and H(Γ) is the Hausdorff surface measure of the crack set Γ. In Equation (11), the functional E has a structure identical to that for image segmentation developed by Mumford and Shah [332]. It consists of the strain energy stored in the solid as well as the energy release due to fracture.

Regularized Variational Theory
The numerical evaluation of the sharp crack interface in the functional E (Equation (11)) is not suitable within a standard finite element framework, as outlined in the work of Bourdin et al. [258]. Therefore, a regularized crack interface using a specific regularization profile γ is introduced in the studies of Miehe et al. [141,297]. It is based on a geometric regularization of sharp crack discontinuities that is governed by a crack PF It characterizes locally for the initial condition φ = 0 the unbroken and for φ = 1 the fully broken state of the material. Thus, the critical fracture energy is approximated by (13) in terms of the crack surface density function per unit volume of the solid. The regularization is governed by a fracture length scale l f . Note that the limit for vanishing the fracture length scale l f → 0 gives the sharp crack surface Γ. Therefore, the minimization problem represented by Equation (11) can be expressed in the following form, defined in terms of the total work density function W as contains a degraded elastic work density and the crack energy release per unit volume. g(φ) is a degradation function defined as g(φ) = (1 − φ) 2 . It describes the degradation of the solid with the evolving crack phase-field φ, as depicted in Figure 8b-d.

Discussion and Conclusions
Based on the above literature review, it can be observed that PF methods have a great potential for simulating self-healing mechanisms in concrete. Therefore, it can be applied to solve problems that cannot be addressed by commonly applied models. It has the potential of an unprecedented breakthrough. As self-healing of concrete is a rather complex process, it is an interaction between physical, chemical and mechanical mechanisms. Obtaining a novel, versatile model for self-healing concrete is a multidisciplinary study involving civil engineering, materials science, and chemistry.
Many studies have been conducted in these fields using the PF approach, while it will be a great reference for the development of a self-healing PF model.
In future research, it would be recommended to include in the polynomial system of the PF approach the pore structure, concrete matrix, water dissolution, and hydration product phases at the crack front. In this way, the free-energy equations will combine hydration kinetics, crystallization kinetics, polymerization reaction kinetics, mass transport and chemical energies to provide a detailed description of the phase nucleation and growth mechanisms at the crack front. Coupling a reactive PF model with a fracture PF model allows to simulate the crack development and its mechanical self-healing recovery effects at different stages and under different environmental conditions. In order to achieve this goal, there are several self-healing mechanisms that need to be studied in great detail. Validating these models should be continuously done by comparing them with experimental results. The following potential future steps are identified: 1. Evolution of the pore structure at the crack surface: During the process of autonomous self-healing, soluble substances at the crack surface enter the solution and undergo various dissolution reactions, followed by hydration and carbonation crystallization reactions. Part of the solution will diffuse into the capillary pores of the concrete matrix, where crystallization and precipitation also occur. The growth of the cracked surface also forms a new pore structure, which further affects the diffusion and chemical reaction processes. Thus, the pore structure of the crack boundary is constantly changing with ongoing reaction. Its interaction with the crack morphology, reactant concentration, and mass transport needs to be investigated in the future.

Influencing factors and simulations for mechanical repair of cracks:
The fracture PF part is a combination of elastic and fracture energies. Elastic free energy will follow the classical assumptions while the fracture part will account for the fracture toughness, order formulation, evolution equations, and healing regain laws. Moreover, both are closely related to the packing density field. This is because the mechanical properties at fracture mainly depend on the solid-phase continuity. The mechanical properties are enhanced in a homogeneously dense position of the filler and, conversely, worse in the disconnected parts of the solid phase. The packing density field, in turn, is related to the mass transport. Therefore, a numerical transport-mechanical coupling strategy shall be developed to simulate the overall performance of the self-healing mechanism.
3. Evolution of crack healing morphology: The morphology of the crack greatly influences its local healing effect. At the crack tip, healing products are produced faster and more frequently because of the higher concentration of reactants. The movement of the crack tip is faster than at other locations. Thus the crack morphology changes continuously with the healing process. As the PF model avoids tracking the boundary conditions at the interface and instead simulates the evolution of the auxiliary field. Therefore, the evolution of the interfacial morphology is easier to simulate. In addition, the simulation of interfacial morphology will take into account the distribution of bacteria, adhesive agents and mineral admixtures. Therefore, the macroscopic representation of a crack healing morphology shall be simulated from a micro-level point of view.
4. Free energy to distinguish between various product phases: Self-healing products contain multiple substances (CSH, CH, or additional byproducts) that, although they have the same healing mechanism (aggregation, crystallization and precipitation), their chemical reaction kinetics are different. This affects the rate of healing of the cracks as a whole. Therefore, the free energies of the various product phases and the corresponding thermodynamic parameters will be distinguished in the future and reflected in specific simulations.

Determination of PF parameters:
A formulation for the determination of the PF parameters needs to be provided. Information on the PF parameters and their interrelationships will be obtained from thermodynamic and diffusion databases in combination with experimental data. Combined with the second law of thermodynamics and non-equilibrium thermodynamics, the self-diffusion, mutual diffusion, and chemical diffusion coefficients will be related to the diffusion mobility (M). The order parameter mobility (L) will be derived and their relationship to other phase-field parameters will be investigated. 6. Development of a three-dimensional model: As a self-healing process includes complex physical-chemical-mechanical processes, these mechanisms can only be accurately simulated in a fully three-dimensional system. Therefore, a three-dimensional simulation of the self-healing process need to be performed with realistic boundary conditions. The simulation results need to be verified and compared with 3D computed tomography scan (CT scan) results of concrete specimens.
In conclusion, the use of a PF method is feasible and has a significant application advantages in the field of self-healing concrete applications. Although this method still has a long way to go before it becomes a fully fledged simulation tool, these early studies are considered to be an important step towards reaching this goal.