Compressive response of self-healing polymer foams containing bilayered capsules: Coupled healing agents diffusion and stress simulations Journal of the Mechanics and Physics of Solids

Numerical modelling on the compressive response of self-healing polymer foams embedded with novel, bilayered alginate capsules was developed based on the coupled pore fluid diffusion and stress simulations. Micromechanical models were developed to link the damage variable to permeability as well as the saturation to the capillary pressure within damaged polymer foams. These micromechanical models were calibrated against experimental measurements and were implemented into the coupled simulations. To give physical insight into how the damage evolution coupled with the mass conservation, an illustrative example was presented for one- dimensional (1D) compression problem. Two-dimensional (2D) detailed finite element simulations were conducted to interpret the experimental findings. It was demonstrated that the nu- merical study could capture the main features of the self-healing process. The predicted healing efficiency has good agreement with that measured by the experiments. Based on the numerical models, parameter study was conducted to understand the effects of the key design parameters of the healing system.


Introduction
Self-healing materials, which can repair themselves automatically, have been emerged to provide resilient solutions for critical engineering structures. For self-healing polymer composites, the self-healing solutions can be achieved through intrinsic and extrinsic self-healing mechanisms (Blaiszik et al., 2010;Wang et al., 2015). Intrinsic self-healing is activated by the parent polymers under external stimuli, including various thermo-mechanical/chemical stimuli (Yu et al. 2019). Extrinsic self-healing can be achieved through releasing the prefilled healing agents within either vascular (Trask et al., 2007;Hansen et al., 2011) or capsular (Hia et al., 2016;Al-Mansoori et al., 2017;Sun et al., 2019) containers when damage occurs. Capsular self-healing systems are suitable for industrial-scale production with the advantages of easy fabrication, low cost, and versatility (Wang et al., 2015). The healing agents can be encapsulated through three approaches, i.e. dual-capsule system (Wang et al., 2017), capsules-catalysts system (Brown et al., 2005), and mono-capsules system (Caruso et al., 2007;Al-Mansoori et al., 2017). Owing to mutually reactive nature, the two-part healing agents must be stored separately in the healing system. This may make the healing system less efficient as the released epoxy resin and hardener may not be able to be well-mixed to achieve good healing effects. To overcome the difficulty, Cao et al. (2020) have recently reported an environment-friendly, multi-stage encapsulating process that can encapsulate the mutually reactive healing agents within single calcium-alginate capsules. The resulted capsules have a bilayered microstructure with the outer layer containing hardener and the inner layer containing epoxy resin, see Fig 1 (a). The capsules have been used to create self-healing polymer foams, see Fig. 1. Experimental study has been conducted to evaluate the self-healing performance of the self-healing system, including (1) cyclic quasi-static compression tests for foam samples; (2) quasi-static three-point bending for foam core sandwich beams; and (3) high-speed soft impact for foam core sandwich beams. These experimental evaluations have demonstrated that the bilayered capsule systems can achieve better self-healing effect compared to the dual-capsule system. The multicore-like internal microstructures of the bilayered capsules enable multiple healing events owing to release of the residual healing agents stored in the undamaged pores inside damaged capsules. In comparison, the capsules with core-shell structure (Yuan et al., 2008) can only show single healing event due to the complete release of healing agents after breaking of capsules.
For the purposes of interpretation of experimental findings and system optimisation, this paper aims to develop numerical modelling on the response of self-healing polymer foams containing bilayered capsules under cyclic compression. Although capsules based self-healing materials have been under intensive investigation for the last two decades, only limited work has been presented in the theoretical analyses and numerical modelling of self-healing systems (Balazs, 2007). Verberg et al. (2006) and Alexeev et al. (2007) modelled the motion of a fluid-driven microcapsules along the damaged adhesive substrate. The lattice Boltzmann model was used to simulate the releasing of nanoparticles from capsules, which had the capability of healing, and the lattice spring model was used to simulate the motion of elastic solids of capsules. The results suggested that the nanoparticles could heal the damage and enabled the motion of capsules along the substrate. Huh et al (2000) and Lee et al (2004) developed the three-dimensional computational model via Monte Carlo (MC) simulations to predict the healing effects of the polymer composite embedded with micro-particles. The model integrated the multilayer solids that combined crack, polymer matrix and micro-particles. The lattice spring model was employed to determine the mechanical properties of the polymeric system in the undamaged, damaged, and healed situation. The results showed that the particles in the polymer matrix potentially could enter into the localized damage regions and healed the damage. Mathematical models have been developed to simulate transportation events of healing agents. Huang and Ye (2014) proposed a model to predict the healing effect of self-healing cement via embedding capsules that were filled with water. When microcracks passed through capsules, the water was released from the capsules and triggered the hydration to seal the cracks. The constitutive behaviour of water transportation was expressed by the mass balance, and the hydration process was simulated via the Fick's second law of ion diffusion. Zemskov et al. (2012) developed a mathematical model for self-healing cracks in concrete. The healing agents were encapsulated within the spherical clay capsules and the healing effect was triggered when released healing agents encountered the bacteria that were pre-embedded in the concrete matrix. The diffusion of the released healing agent was solved by the Galerkin finite element method.
In this paper, the self-healing events are simulated via coupling the system equilibrium equations and mass conservation for the healing agents flowing inside cracked media. This approach has been used to simulate the mechanical responses of saturated or semisaturated porous media (Selvaduria and Suvorov, 2016). However, it has not been used to simulate the capsule based self-healing events. We consider that the closed-cell polymer foams are impermeable before damage. After onset of damage, the healing agents are released from the capsules and flow inside the cracked polymer foams to heal damage. To enable the simulations, micromechanical models were developed to link damage to permeability of the cracked polymer foams as well as the saturation to the capillary pressure within the damaged polymer foams. It was demonstrated, through the comparisons between the detailed two dimensional (2D) finite element simulations and experimental measurement, that the computational model could capture the main features of the self-healing process. The present model has the following advantages: 1. It incorporates the key concepts in Continuum Damage Mechanics (CDM) into the detailed simulations of self-healing events, which enable detailed comparison with experimental measurements. 2. The simulations can be easily implemented using the finite element solvers for the coupled pore fluid diffusion and stress analysis within the commercially available finite element package such as ABAQUS through user defined subroutines.
The paper is organised as follow: In Section 2, the main outcome of the experimental study is summarised; In Section 3, the framework of coupled pore fluid diffusion and stress simulations is described; In Section 4, the micromechanical models are developed to facilitate the coupled analysis; In Section 5, an illustrative example is given for one dimensional compression problem; and in Section 6, the results obtained by 2D finite element simulations are presented.

Bilayered alginate capsules
The capsules were formed through the chemical reaction when sodium alginate 1 solution encountered calcium chloride 1 solution. The two-part healing agents, i.e. rapid repair epoxy resin 2 and formulated amine hardener 2 , were encapsulated within single calciumalginate capsules using the two-stage encapsulation method developed by Cao et al. (2020). Briefly, in the first stage, the well-mixed epoxy resin-alginate emulsion was dropped into the calcium chloride solution to form single layer alginate capsules containing the epoxy resin; in the second stage, the dried single layer epoxy resin capsules were mixed with alginate-hardener emulsion. The emulsion covered epoxy resin capsules were dropped into the calcium solution to form the outer alginate-calcium layer containing the hardener agent. After completion, the diameters of the bilayered alginate capsules ranged from 3.5 mm to 4 mm. The scanning electron microscope (SEM) image showing the Internal structure of the bilayered alginate capsule is presented in Fig. 1 (a), which suggests that the bilayered microcapsule has multicore-like internal structure with pore diameters in the outer (hardener) around 150μm, and pore diameters in the inner (resin) around 70μm.

Polymer foams
The polymer foams were created by mixing PB250 epoxy resin 3 and DM03 hardener 3 , as detailed by Cao et al. (2020). The cured polymer foams had the closed-cell microstructure with the density ρ e = 0.33g /cm 3 , as the SEM image shown in Fig. 1 (c). The porosity of the foams, defined as the ratio of the volume of voids to the volume occupied by a bulk foam, was 0.67. To fabricate foams embedded with capsules ( Fig. 1 (b)), the bilayered calcium-alginate capsules were added to the mixture at the beginning of the foaming process, which can be uniformly distributed within the foams.

Experimental protocol
Cyclic quasi-static compressive tests of polymer foams were conducted using a screw-driven Instron Universal Testing Machine at the room temperature following the procedure defined by ASTM 1621 -04a, as described by Cao et al. (2020). The tests were conducted on cubic specimens of original edge length l f = 50 mm at a crosshead speed of 5 mm/min and with three repetitions for each type of foams, i.e. neat polymer foams and polymer foams embedded with bilayered capsules. The measured compressive force P f and the vertical displacement Δl of the crosshead were recorded by a 50KN load cell and Linear Variable Differential Transformer (LVDT), respectively. The engineering compressive stress σ c and strain ε c of foams can be calculated as σ c = P f /A f and ε c = |l a − l f | /|l f |, respectively, where l a is the height after compression, and A f = (l f ) 2 the original area of the cross section of the samples. Fig. 2 (a) shows the definition of a compressive loading cycle which includes a loading phase, an unloading phase and a healing phase. It is worth noting that the strain recovery during a healing phase is not attributed to the presence of the capsule based healing system as the strain recovery has been observed in the neat foam samples. After damage occurred, the healing agents were released from the capsules and  flowed into damaged areas to repair microcracks ( Fig. 1 (b)). To allow stable healing recovery, the interval time between two adjacent compression cycles was 24 hours. To achieve adequate incremental damage, the applied strain was controlled based on the relation ε n = ε * n− 1 + 0.075 for nth compressive loading cycle, where ε * n− 1 is the residual strain measured after the unloading phase in the (n-1)th compression cycle, and ε * 0 = 0. The damage evolution of a specimen is quantified by a scalar parameter of damage variable D, where E da denotes the elastic modulus of a damaged specimen measured at the unloading phases, and E vi the initial elastic modulus. To evaluate the self-healing effect, we define a healing cycle as the period from the unloading phase in the (n-1)th compression cycle to the loading phase in the nth compression cycle, highlighted by the solid curve in Fig. 2 (b). Hence, the healing efficiency ψ is defined as where E he is the elastic modulus after healing measured at the loading phase of nth compression cycle.

The key experimental results
The measured cyclic quasi-static compressive responses of the foam embedded with 10% volume fraction (VF) bilayered capsules and the neat foam are shown in Fig. 3 (a) and (b), respectively. The compressive strength and elastic modulus of the foam embedded with capsules measured at the first compression cycle are higher than those of the neat foam, which indicates that the capsules can provide reinforcing effect to the foam matrix. Fig. 3 (c) and (d) show the SEM micrographs of a microcrack completely healed by the released healing agents during the first healing cycle, which suggests that the current healing approach is effective in healing microcracks. To evaluate the healing effect, the self-healing efficiency ψ is calculated based on the experimental measurements. The values of healing efficiency ψ measured at selected healing cycles are shown in Fig. 3 (e) and (f) for the two types of foam samples, respectively. For the neat foam sample, the healing efficiency ψ is non-zero. The absorption of air during the unloading phase could be the reason (Michal Petrů and Ondřej Novák, 2017). Therefore, the healing effects can be evaluated based on the threshold value ψ s = 0.234 that is equivalent to the maximum value of the healing efficiency of the neat foam sample, i.e. there is healing effect if ψ > ψ s , otherwise, no healing effect. For the foam sample embedded with bilayered capsules, the healing efficiency ψ at both the 1 st and the 2 nd healing cycles is higher than ψ s , which indicates the presence of the multiple healing effect. The multiple healing effect diminishes in For the neat foam sample, the response measured by the cyclic compression test ( Fig. 3(b)) suggests that the cyclic compression causes degradation of elastic modulus, albeit there is negligible degradation to the plateau stress σ c . To understand the damage evolution during the cyclic compression, the relationships between the damage variable D and the plastic strain ε p of the tested specimens are shown in Fig. 3(g). Owing to the presence of the self-healing mechenism, the damage variables for the foam embedded with bilayered capsules are smaller than those of the neat foam sample. The numerial modelling of the self-healing polymer foams under compression will be discribed next.

Coupled pore fluid diffusion and stress simulations
As schematically shown in Fig. 4, a foam was initially modelled as an isotropic material without damage. After onset of damage, the healing agent flowed from capsules to the damaged areas, which was dictated by mass conservation. The foam was treated as the porous medium containing incompressible fluid phase (the healing agents), and the interaction between the solid phase and fluid phase of the porous medium was dictated by the Darcy's law. The system equilibrium equations and mass conservation for the healing agents flowing inside cracked media were coupled to solve the structural responses.

Equilibrium equations
Consider a closed cell foam containing bilayered micro-capsules as the healing agent carriers, as schematically shown in Fig. 4 at the reference and deformed configurations without and with damage. As the closed cell foam without damage is nearly impermeable, the healing agent is not able to flow inside the foam sample. Hence, the effect of healing agents to structural response of the foam is negligible, see Fig. 4 (b). Let Ω(t) represents the volume occupied by a bulk foam material at the current configuration with the boundary S Ω (t) at time t under Cartesian coordinates x 1 − x 2 − x 3 ; σ the effective stress tensor (Cauchy stress), i.e. σ = 1 Ω ∫ Ωs σ s dΩ, where Ω s denotes the volume occupied by the skeletons of the foams; σ s the micromechanical variation of the stress field within Ω s ; and f the body force. The equilibrium equation of the system without damage reads where ∇ = ∂/∂x with x representing the spatial coordinates at the current configuration. Throughout the paper, the inertial effect to the equilibrium equation is ignored as the whole process is quasi-static. For the foam with damage ( Fig. 4 (c)), microcracks at damaged locations provide the micro-channels that allow healing agents to flow through. Let Ω c represents the total volume of microcracks at the current configuration, Ω v the total volume of voids within the foam before onset of damage, and Ω w the volume of the healing agents that flows through the microcracks (Ω w ≤ Ω v ). The concentration of microcracks ξ is defined as ξ = Ω c /Ω, with ξ ≤ 0.03 as measured by the current research. The porosity n p of the foam can be calculated as n p = Ω v /Ω + ξ with ξ = 0 for an undamaged foam. For the foam with original porosity of 0.67 as considered in the current research, the effect of ξ to n p is negligible. Hence, throughout the numerical study, it is reasonable to assume n p ≈ Ω v /Ω. Let saturation s denotes the ratio of the volume saturated with healing agents to the total volume of voids, i.e.s = Ω w /Ω v ,s ∈ [0,1]. With presence of the healing agents, the equilibrium equation of the system reads where ρ w is the density of the healing agents; g the gravity acceleration; and σ the total stress (Cauchy stress) at a material point (Skempton, 1960;Lade and De Boer, 1997), which can be defined as where p w is the pressure of healing agents; p a the air pressure trapped within the closed cells; I the unit tensor; and χ the Bishop parameter that is a function of saturation s (Bishop, 1959). To simplify the problem, we ignore the effect of the air trapped in the closed cells within the foam (Hou et al., 2015). Hence, Eq. (3.3) can be rewritten as Gray and Schrefler (2001) suggested that the Bishop parameter χ could be replaced by the saturation s. Therefore, Eq. (3.4) can be further rewritten as

Mass conservation for the healing agents flowing in cracked media
The mass conservation dictates that the rate of mass of the healing agents flowing into a volume equates the rate of mass increase of the healing agents stored at the volume. Consider a volume Ω f (t) at the current configuration at time t with surface S f (t).Mass conservation can be described mathematically as ∫ where v w and ρ w are the velocity and density of the healing agents, respectively; J the determinant of deformation gradient F,J = detF; and n the outward normal to surface S f . Using the divergence theorem, Eq. (3.6) can be rewritten as The volume rate of healing agents flowing through a surface with area A s ,Q(mm 3 /s), is assumed to be governed by the Darcy's law (Bear, 1972) where k(mm 2 ) is the permeability of the damaged foam which will be determined in Section 4.1;η(Pa • s) the dynamic viscosity; and ∂p w /∂x the pressure gradient (Pa /mm).The negative sign in Eq. (3.8) indicates that the fluid flows from the region of high pressure to the region of low pressure. The averaged velocity of the healing agents flowing through the effective unit area can be defined as (Bear, 1972) Hence, Eq. (3.7) can be rewritten as The coupled pore fluid diffusion and stress analysis can be conducted via simultaneously solving Eqs. (3.2) and (3.10). To complete the formulation, the functional relationships between the damage variable of a damaged foam D and the concentration of microcracks ξ; the damage variable D and the permeability k; and the saturation s and the capillary pressure of the healing agents p w are described in Section 4.

Permeability of a cracked foam
The bulk permeability k of a cracked foam needs to be determined to enable the simulations of healing agents flowing inside a cracked foam (Eq. (3.10)). Fig. 5 shows a SEM image of a microcrack within a damaged foam sample, which contains fractured hollow cells. The microcrack provided a micro-channel that allowed the flowing of the released healing agents within the damaged area. This observation motivated the following micromechanical model that could estimate the value of bulk permeability k. To simplify the problem, we assume that all microcracks are of cylindrical geometry. It is a reasonable assumption as observed through SEM images, such as Fig. 5. The effect of the shape of the microcracks will be explored in the future study. Consider an isolated crack of cylindrical geometry and constant cross section embedded in a cubic volume as shown in Fig. 5. Let r and l represent the radius and length of a microcrack, respectively; r ∈ [0, R] with R representing the possible maximum crack radius; and l ∈ [0, L], with L representing the possible maximum crack length. In the Cartesian coordinate system x 1 − x 2 − x 3 at the current configuration, the projection of the crack in x 1 − x 2 plane has an angle to axis-x 2 , φ(0 ≤ φ ≤ 2π); and the crack has an angle to axis-x 3 , θ( − π /2 ≤ θ ≤ π /2). Experimental measurement suggests that it is reasonable to assume that both angles follow uninform distributions, i.e. θ ∼ U[− π /2, π /2] and φ ∼ U [0,2π].This assumption dictates that the bulk permeability k of a cracked foam containing a random ensemble of cracks is isotropic. In the following analysis, to simplify the problem, we consider the scenario that the bulk flow through the foam is only in the x 3 -direction.
For an incompressible and Newtonian fluid in the laminar flow flowing through a cylindrical crack of a constant cross section, the volume flow rate Q s can be related to the pressure drop Δp w through the Hagen-Poiseuille law, i.e.
where n 0 is a unit vector along the axis of the crack; and Δp w the pressure difference between the two ends of the crack. As the fluid flows from the high pressure p h to low pressure p l , the negative sign comes from the definition of Δp w = p l − p h < 0. As the length of the majority of microcracks is smaller than 10 mm, it is reasonable to assume Δp l ≈ dpw dl . In the Cartesian coordinate system, Eq. (4.1) can be rewritten as As the bulk flow through the foam is assumed to be only in the x 3 -direction, we have ∂p w /∂x 1 = 0 and ∂p w /∂x 2 = 0. Eq. with where e i (i = 1, 2,3) is the unit vector in x i -direction, and |v w | the length of the vector v w . For a random ensemble of cracks saturated with healing agents, the average flow velocities can be calculated as where f 1 (r) denotes the probability density function of the radii of cracks, and 〈•〉 the mean value of a random variable. Eq. (4.6) suggests that the average flow velocity 〈v w3 〉 for the ensemble of saturated microcracks is not related to the length of the cracks l and the angles φ.The mean value μ r and variance σ 2 r of the radii of the cracks can be calculated as Substituting Eq. (4.9) through the third equation of Eq. (4.6), we have Recall the Darcy's law given in Eq. (3.9) and the bulk flow through the volume is only in the x 3 -direction, i.e.〈v w3 〉e 3 ≡ v w . For the ensemble of fully saturated microcracks,sn p ≡ ξ. Hence, we have where k s denotes a coefficient that is introduced to take into account the effect of saturation with k s = 1 for fully saturated media. Nguyen and Durso (1983) suggested k s = s 3 . To calculate the permeability of a cracked foam via the Eq. (4.11) at the nth compressive loading cycle,μ rn , σ rn and ξ n can be estimated based on histograms of radii and lengths of the microcracks measured from the tested samples using the following formula where N n is the total number of microcracks in the foam sample at the nth loading cycle.

Capillary pressure and saturation of cracked foams
For porous media such as the polymer foams considered in the paper, capillary pressure plays an important role owing to the high porosity of the porous media. In the following analysis, the functional relation between saturation s and capillary pressure p w in the porous media will be established. Consider the isolated microcrack of cylindrical geometry and the constant cross section embedded in a cubic volume, as described in Section 4.1 and shown in Fig. 5. The capillary pressure in the microcrack can be calculated as (Fanchi, 2002) where τ s is the surface tension, and ω the liquid-solid contact angle on the walls of a micro-channel. For the epoxy system, τ s = 0.00003475N /mm and ω = 43 0 55 ′′ (Cheng and Lin, 2008). Let r * n , (r * n = r n /R n ; 0 ≤ r * n ≤ 1), denote the normalised radius, and l * n , (l * n = l n /L n ; 0 ≤ l * n ≤ 1), the normalised length of a microcrack at the nth loading cycle, where R n and L n are the maximum radius and length among all cracks within a foam sample at the nth loading cycle. In the following analysis, we assume r * n and l * n follow the Beta Distribution, i.e.
where (α r * n , β r * n ) and (α l * n , β l * n ) are the shape factors for f 1 (r * n ; α r * n , β r * n ) and f 2 (l * n ; α l * n , β l * n ) respectively. The change of saturation ds n at the nth compressive loading cycle induced by the healing agents flowing into all microcracks of radius r * n can be calculated as where μ l * n is the mean normalised length of microcracks. Eq. (4.17) can be further written as where μ r * n is the mean normalised radius of microcracks. The total volume of microcracks Ω c can be calculated as where μ r * 2 n is the second moment normalised radius Substituting Eq. (4.21) through Eq. (4.20) and from (4.16), we have Let p * n represent the normalised capillary pressure in the nth loading cycle,p * n = p w (r n ) /p w (R n ). From Eq. (4.15), we have Integrating Eq. (4.24), we can establish the functional relation between capillary pressure and saturation for the nth loading cycle.

The model for damage recovery after healing
The recovery of damage after healing at a material point can be related to the saturation of the healing agents at the point during healing. Here, we assume the recovery of damage as a linear function of the saturation, i.e.
where D re denotes the residual damage variable post healing; and s n re the saturation which enables the complete recovery of damage at the nth compression cycle. s n re was obtained through calibration against experimental measurement, See Appendix C.

Data fitting and model calibration
The foam samples were examined by an optical microscope after each compressive loading cycle to measure the distribution and geometry of the microcracks. The measurements were used for data fitting and model calibration for establishment of the functional relationships between (1) the damage variable D and the concentration of microcracks ξ; (2) the damage variable D and the permeability k; and (3) the saturation s and the capillary pressure of the healing agents p w . In addition, the evolution of damage variable D against the plastic strain ε p was also obtained through the curve fitting against experimental measurement.

The damage variable D and plastic strain ε p
As shown in Fig. 3 (g), the damage vairable D envolves as a function of plastic strain ε p . The functional relation can be discribed as a power law fuction through data fitting to the experimental data measured for the neat foam where C 1 andC 2 are constants obtained from curve fitting as listed in Table 1.

The damage variable D and the crack concentration ξ
The experimental measurements suggest that the functional relation of cracks concentration ξ and damage variable D follows an exponential equation as given below (4.28) where C 3 andC 4 are constants obtained from curve fitting as listed in Table 1. The curve fitting of Eq. (4.28) is plotted in Fig. 6 (a) in comparison with experimental measurements. Feng and Yu (2010) employed an analytical model and FE analysis to investigate the relation between the damage evolution and the cracks growth of quasi-brittle isotropic materials. The results indicated that the damage showed an exponential growth with the increase of microcrack density. Farrokhabadi et al. (2013) investigated the damage evolution of laminated composites under fracture tests experimentally and numerically. The damage variable and crack density again exhibited an exponentially functional relationship. The trend of the functional relation in Eq. (4.28) is consistent with the findings from these existing research.

The damage variable D and the permeability k
For the ensemble of fully saturated microcracks (k s = 1), based on the Eq. (4.11), the permeability k can be estimated via measuring the lengths and radii of the microcracks from the 2D micrographs and cracks concentration ξ at each loading cycle. Eq. (4.11) also suggests that permeability k is nearly proportional to microcracks concentration ξ as the mean value μ r and variance σ 2 r of the radii of the microcracks have small variations at different loading cycles. As the damage variable D follows the exponential relation with respect to the concentration of microcracks ξ(Eq. (4.28)), the functional relation between the damage variable D and the permeability k can also be described as an exponential equation, i.e. k = C 5 e (C4D) (4.29) where C 5 is a constant obtained from curve fitting as listed in Table 1. The fitted relationship of Eq. (4.29) is plotted in Fig. 6 (b).

The saturation s and the capillary pressure of the healing agents p w
The constitutive relation between saturation s and capillary pressure p w can be described via the Beta distribution function of the normalised radii of the microcracks within the foam samples (Eq. (4.16)). The histograms and the Beta distribution fitting of the normalised radii of the microcracks after each loading cycle,Beta(α r * n , β r * n ), are presented in Fig. 7.
Based on Eq. (4.25) and measured values of α r * n , β r * n , ξ n and n p , the functional relations between s and the capillary pressure p w for each loading cycle are plotted in Fig. 8. As the difference of these functional relations among each loading cycles is small, we use a single exponential relationship to represent the functional relation between saturation s and capillary pressure p w , i.e. s = C 6 e C7pw (4.30) where C 6 and C 7 are the constants from curve fitting as listed in Table 1.

An illustrative example -one-dimensional (1D) compression problem
To give the physical insight into how the damage evolution coupled with the mass conservation, an illustrative example is given for one dimensional compression problem. Consider a foam containing two bilayered capsules under uniaxial cyclic compression, as schematically in Fig. 9 (a), the healing agents is flowing inside the damaged areas from the embedded bilayered capsules. Under the Cartesian coordinates x 1 − x 3 at the current configuration, we take a small region of thickness h s within the foam for analysis, as schematically shown in Fig. 9 (b). For ease of calculations, the following assumptions are made, i.e.
(1) During the calculations, the following pore pressure boundary conditions are assumed p w (t, x 3 = h s ) = p 0 for ∀t ≥ 0 (5.1)  where p 0 denotes the pore pressure imposed by the healing agents released from the capsules, and p i the initial capillary pore pressure. It is reasonable to assume the p i = 0.00065 MPa for the dried porous media (Jung et al., 2016). The assumptions implicate that we only consider the situation before the healing agents reach the bottom surface.
(2) The foam follows Hooke's law in the elastic region with initial elastic modulus E vi , and linear hardening in the plastic region with hardening modulus H when the compressive stress exceeds the compressive strength σ c . The evolution of the damage variable D of the foam only affects the value of elastic modulus, i.e. E da = (1 − D)E vi (Eq. (2.1)), and has negligible effect on the values of the compressive strength σ c and hardening modulus H. This assumption is consistent with experimental measurement of foam samples under uniaxial compression (Fig. 3), i.e. the damage of foam has negligible effect on the plateau stress. (3) The thickness of selected region h s is much smaller than the thickness of foam domain h f , i.e. h s ≪ h f . Therefore, the strain rate of the selected region during the deformation could be assumed to be a constant value ±C. Hence, a compressive loading cycle starting from time instant t o includes the following temporal phases (Fig. 2 (a)), i.e.
• The loading phase, including • The healing phase where ε e is the elastic strain.
(4) Owing to small thickness h s , the effect of saturation s to permeability k is ignored, i.e. k s = 1.

1D Constitutive behaviours
The response in elastic loading phase t ∈ [t o , t o +T e ] is governed by the Hook's Law, i.e.
where s o and ε o p are the saturation of healing agents and plastic strain at t = t o . In the plastic loading phase t ∈ (t o + T e ,t o + T e + T p ], the response follows linear hardening relation, i.e. where σ is the stress rate, ε the strain rate, and ε p the plastic strain rate. Combining the two equations in Eqs. (5.5) and (4.27), we have Eq. (5.6) can be solved using the fifth-order Runge-Kutta method based ordinary differential equation solver ode45 in MATLAB with the initial condition ε p (t = t o +T e ) = ε o p and ε = − C. The calculated plastic strain ε p can be fitted by a polynomial function of time t, i.e.
where C 8 and C 9 are two constants obtained from curve fitting. In the unloading phase t ∈ (t o + T e + T p , t o + 2T e + T p ), there is no further development of plastic strain ε p . Hence, we have where ε L p is the plastic strain at t = t o + T e + T p . The damage recovery model (Eq. (4.26)) is applied for healing phase, i.e. t = t o + 2T e + T p , to enable recovery of elastic modulus, i.e.
where s L is the saturation of the healing agents at t = t o + 2T e + T p .

Mass conservation for the healing agents
When damage occurs, the foam matrix becomes permeable. Under 1D cyclic compression as shown in Fig. 9 with X 3 -axis as the loading direction, a material point located at the coordinates (X 1 , X 2 , X 3 ) under the initial configuration will deformed to the position (x 1 , x 2 , x 3 ) under the current configuration. The expansion of the mass conservation (Eq. (3.10)) for the healing agents flowing in a cracked media can be written as The deformation gradient is given by , and the determinant of F is In Eq. (5.11), the material time derivative operator can be written in the form where v f is the velocity field in the either solid phase (∂x 3 /∂t) or fluid phase (v w ). Here, we use the s = C 6 e C7pw (Eq. (4.30)) to describe the relation between saturation and capillary pressure at the selected region. Combining Eqs. (4.30), (5.11), (5.13) and (5.14), we have As the strain rate of foam is a constant C throughout the selected region, we have ∂ε /∂x 3 = 0, ∂ε /∂t = − C in loading phase, and ∂ε /∂t = C in the unloading phase. Combining with the Darcy's law in Eqs. (3.9), (5.15) can reduce to

±C
(1 + ε) n p C 6 e C7pw + n p C 6 C 7 e C7pw As damage variable D is related to plastic strain ε p through Eq. (4.27), which is independent to the spatial coordinates (Eq. (5.7)).
In the unloading phase t ∈ (t o + T e + T p ,t o + 2T e + T p ), the value of permeability remain unchanged. Eq. (5.18) can be written as At the healing phase t = t o + 2T e + T p , combining with Eqs. (4.26) and (5.18) can be written as where D L re denotes the damage variable at t = t o + 2T e + T p . The Eqs. (5.19), (5.21)-(5.23) were solved using forward Euler finite difference (FD) method as described below ∂p w ∂t where Δt was a very small time incremental that enabled convergence of solutions. From Eq. (5.25), the Eqs. (5.19), (5.21)-(5.23) can be rewritten as the where A 1 = C 5 e C4D /ηn p C 6 C 7 e C7pw for the loading phase, A 1 = C 5 e C4D(ε L p ) /ηn p C 6 C 7 e C7pw for the unloading phase, and A 1 = C 5 e C4D L re /ηn p C 6 C 7 e C7pw for the healing phase;A 2 = C 5 e C4D /C 6 e C7pw n p η for the loading phase, A 2 = C 5 e C4D(ε L p ) /C 6 e C7pw n p η for the unloading phase, and A 2 = C 5 e C4D L re /C 6 e C7pw n p η for the healing phase; and A 3 = C C7(1− ε) in the loading, and A 3 = − C C7(1− ε) in the unloading phase and healing phases. Substituting Eq. (5.25) into Eq. (5.26), we have In the forward Euler FD calculations, we assumed the small region had height h s = 10 mm; the foam had porosity n p = 0.67, the initial elastic modulus E vi = 300 MPa, hardening modulus H = 10 MPa, and compressive strength σ c = 6.12 MPa; the released healing agents had the dynamic viscosity η = 0.000001 MPa•s; and strain rate C = 0.001/s . From curve fitting, we have C 8 = − 4.0 ×10 − 7 /s 2 and C 9 = 1.1 × 10 − 3 /s. For comparison purpose, the 1D problem was also solved by finite element (FE) simulations using the commercially available FE package ABAQUS Standard. The FE discretization was conducted using 252 four-node plane strain elements, i.e. CPE4P in ABAQUS notation, which were based on the FE formulation of coupled stress analysis and fluid diffusion within a porous media as described in Section 6. The FE model employed in the simulation is shown in Fig. 9 (c).

Results and discussion
The numerical simulations were conducted to simulate the response of the foam sample under the cyclic compression with the temporal phrases defined in Table 2. Fig. 10 shows the time histories of total strain, plastic strain, capillary pressure, damage variable, and permeability at the selected locations, i.e. X 3 = hs 2 and 3hs 4 , obtained via the FE and FD calculations, respectively. For the FD calculations, numerical experiment suggested that the converged results could be achieve when Δt ≤ 0.000833s. As the strain rate is a constant, the time histories of overall strain, plastic strain, damage variable, and permeability are identical at the two selected locations. However, temporal responses of the capillary pressure are distinct. The numerical results obtained by the FD calculations have good agreements with those obtained by the FE calculations. The discussions on the numerical results are presented as follow.
The closed-cell polymer foam is impermeable during the elastic loading phase in the 1st compression cycle (t = 0 ∼ 15s). The healing agents are not able to flow inside the foam. The damage variable, capillary pressure, and permeability remain unchanged. During the plastic loading phase (t = 15 ∼ 40s), the damage occurs and the damage variable D increases based on Eq. (4.27). The  permeability k increases accordingly with damage variable based on Eq. (4.29), which increases the ability of healing agents flowing inside the foam. Once the damage occurs, the values of the capillary pressure p w at the two locations increase dramatically. According the Eqs. (5.19), (5.21)-(5.23), the temporal variation of the capillary pressure p w relates to compressive deformation and spatial variation of the capillary pressure. During the unloading phase (t = 40 ∼ 55s), the elastic strain is recovered. The plastic strain, damage variable and permeability remain constant values. The capillary pressure p w keeps increasing, which may indicate the healing agents can still flow into the damaged areas owing to spatial variation of the capillary pressure. At the healing phase (t = 55s), both the damage variable and permeability decrease owing to the healing effect. However, as the damage is only partially healed, the healing agents can still flow into the damaged areas with increased capillary pressure. During the elastic loading phases in the following compressive loading cycles (t = 55 ∼ 70s,t = 110 ∼ 125s), the plastic strain, damage variable, and permeability remain constant as no further damage develops. The capillary pressure increases owing to the flowing of healing agents into the unhealed areas. At the beginning of the plastic loading phases (75 s or 110 s), the damage variable and permeability recover to the values before healing in the previous compressive loading cycle. This is consistent to experimental observation: the healing system can only have healing effect on elastic response, and the healed microcracks reopen under the plastic deformation.

Finite element simulations of two-dimensional (2D) problems
6.1. The finite element formulation for the coupled pore fluid diffusion and stress analysis

Weak form of the equilibrium equations
For foams without damage, the structural response (Eq. (3.1)) can be analysed using finite element formulation through the principle of virtual work for the volume Ω at the current configuration at time t (updated Lagrangian formulation) ∫ where t is the traction force on the boundary S Ω (Fig. 4); f the body force;v the spatial velocity; D the effective rate of deformation, i.e. D = 1 /Ω ∫ Ωs D s dΩ, where D s the micromechanical variation of the strain field within the volume occupied by the skeletons of the foams Ω s .
For the foams with damage, microcracks at damaged locations may provide the microchannels that allow healing agents flow. Eq. (6.1) can be rewritten as ∫ where D denotes the total effective rate of deformation at a material point, which can be calculated as where K b denotes bulk modulus of the foams.

Weak form of mass conservation for the healing agents flowing in a cracked media
The equivalent weak form of Eq. (3.7) can be written as ∫ where δp w is an arbitrary, continuous and variational variable, which can be interpreted as virtual pore pressure in the liquid phase. In the reference configuration, we have ∫ where Ω f0 is the volume in the reference configuration. Using backward Euler integration, Eq. (6.5) can be rewritten as ∫ which, over the current volume Ω f , is ∫ Using the divergence theorem Eq. (6.7) can be rewritten ∫ The coupled pore fluid diffusion and stress finite element analysis is achieved by solving Eqs. (6.2) and (6.8) simultaneously. Under the finite element formulation, the nodal variables are {v, p w }, which can be achieved by adding an additional degree of freedom representing the pressure in the liquid phase p w to displacement based finite element formulation. In this paper, the FE solver for the coupled pore fluid diffusion and stress analysis within the commercially available finite element package ABAQUS Standard was employed to conduct the FE analysis, as detailed below.

The finite element model
To create the 2D FE models of foams embedded with bilayered alginate capsules, the capsules were assumed to have circular geometry of an identical radius, which were inserted into the foam matrix using the approach developed by Cao et al. (2019). Briefly, the location of a circular capsule within a foam matrix satisfies two conditions.
• the coordinates (x, y) of the centre of a capsule obey the following relation where r c denotes the radius of the capsule, l f the edge length of the FE model l f = 50mm.
• Capsule i does not overlap with an adjacent capsule j, i.e. ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ ( where n c is the total number of capsules. The four-node plane strain elements (CPE4P in ABAQUS notation) were employed in the FE simulations. Each node of the element had two translational degrees of freedom and one degree of freedom representing pore pressure of the healing agents. Numerical study suggested that the convergence of the numerical simulation could be achieved when the ratio between the maximum elemental edge length and the edge length of the cubic foam sample is less than 1/50. A total of 1391 elements were employed to discretise the solution domain containing both circular alginate capsules and foam matrix.
The isotropic constitutive foam model (Deshpande and Fleck, 2000) in conjunction with a damage model was employed to simulate the mechanical behaviour of the polymer foam, as described in Appendix A. A user defined field routine (USDFLD in ABAQUS notation) was developed to implement the damage evolution model (Eq. (A.5)). The experimental data obtained from uniaxial compression test on neat foam samples (without capsules) were used to calibrate the constitutive model. Fig. 3 (b) shows the comparison between the FE prediction and experimental measurement of a neat foam under cyclic compression. The FE prediction can capture the plastic behaviour and damage evolution, however, it could not capture the strain recovery of the real foam sample after completely unloaded, which might be caused by the recovery of the trapped air within the closed-cells during the unloading phase. To predict the self-healing behaviour, the USDFLD was used the exchange the damage variable (the field variable) with the main program of the ABAQUS to implement the functional relation of the damage variable D and the permeability k (Eq. (4.29)); the damage recovery model (Eq. (4.26)) was implemented into the USDFLD to update the damage variable for healing phase modelling.
As the current research focused on modelling the self-healing behaviour of the polymer foams, to simplify the problem, the constitutive behaviour of the bilayered alginate capsules was modelling as elasto-J2-plastic solids with isotropic hardening. The mechanical parameters of the material were determined through the single capsule compression test, see Appendix B. In the numerical simulations, the bilayered capsule was treated as a fully saturated solid with porosity n p = 0.7 measured via Thermogravimetric Analysis (TGA) (Cao et al., 2020). The permeability of the bilayered alginate capsules is subject to further experimental investigation. In the simulations, the bilayered alginate capsules were assumed to have a large permeability, i.e. k= 0.20mm 2 .
To simulate the compressive test, the 2D FE model of the sample was sandwiched between two rigid plates: one of the rigid plates was stationary, and another rigid plate was movable along the x 3 -axis direction. The movable rigid plate imposed pressure on the sample at a constant velocity (5mm/min for crosshead movement). A penalty contact approach was employed to simulate the interaction between all surfaces with a friction coefficient 0.2, and numerical trial tests suggested that the simulation results were not sensitive to the value of the friction coefficient employed in the calculations. The other parameters employed in the numerical simulations are same to those employed for the 1D simulation presented in Section 5. The outcomes of the 2D numerical simulations are presented next. Fig. 11 shows the predicted response of the foam sample embedded with 10% VF bilayered capsules under cyclic compression. Under the 2D coordinates x 1 − x 3 , the cyclic compressive load was imposed in x 3 -axis direction. The numerical results for the first two healing cycles are presented as the self-healing effect is more pronounced than the later stage. The experimental measurement is included for comparison purpose. It should be noted that the current numerical simulations were based on the 2D plane strain finite element formulation, which may introduce certain discrepancy in comparison with the measured 3D behaviours of the test samples. As  highlighted in Fig. 11 (a) for the stress-strain curves, the segments A-B, B-C, C-D(D'), and D-E represent the first plastic loading phase, first unloading phase, first healing phase, and second loading phase, respectively. Healing (damage recovery) takes place numerically at point D. The FE predictions again could not capture the strain recovery measured during phase C-D'. For the comparison purpose, the healing efficiency obtained via numerical predictions was also evaluated by the threshold value ψ s = 0.234, as shown in Fig. 11 (b). The healing efficiencies obtained by the FE predictions agree well with those obtained by the experimental measurements. Fig. 11 (c) presents an X-ray microcomputed tomography (μCT) image of a foam specimen embedded with bilayered capsules at the end of first compression cycle (Point D). The capsules, healed microcracks and unhealed microcracks are visible in the image as the white particles, white lines, and black lines, respectively. The white lines linking the two adjacent capsules indicated the flowing of the released healing agents through the microcracks. In comparison, Fig. 11 (d) shows the contour of saturation s obtained from numerical predictions at Point D, which shows the similar level of healing agent movement. As saturation s can be related to both healing agent movement (Eq. (4.30)) and recovery of damage (Eq. (4.26)) and damage variable D represents the level of damage within the material, the contours of damage variable D and saturation s at selected strains are presented in Fig. 12 (a) and (b), respectively. These strains correspond to Point A to Point E in Fig. 11 (a). During the elastic region in the first loading phase (before Point A), healing agents could not flow inside the foam matrix without damage. Damage occurred at the first plastic loading phase (A-B): the areas around the capsules had higher damage level owing to stress concentration. The adjacent capsules were linked by the microcacks that provided the microchannel for flowing of healing agents. The majority of the released healing agents concentrated surrounding capsules at the Point B. During the first unloading phase (B-C), the released healing agents flowed inside the damaged areas owing to the capillary pressure, while the contour of damage variable remained unchanged. After the unloading phase (Point C), the healing phase took place (C-D). The damage variable was partially recovered in the entire damaged areas while the contour of saturation remained unchanged. The second loading phase (D-E) resulted in further damage and more released healing agents.

Parameter study on the key design parameters of the healing system
The numerical study was conducted to evaluate the effects of the key design parameters of the healing system, i.e. the volume fraction (VF) of the capsules embedded in a foam sample V f ; the normalised elastic modulus of the capsules E * , defined as E * = E c /E f where E c is the elastic modulus of capsules, and E f the elastic modulus of foam matrix; the normalised averaged diameter of the capsules d * = d c /l f where d c is the averaged diameter of capsules, and l f the edge length of foam samples.

Effect of the volume fraction of capsules V f
To evaluate the effect of V f to self-healing efficiency ψ, the FE simulations were conducted for the foam samples embedded with 5%, 10%, and 15% VF capsules. The self-healing efficiencies of the foams embedded with different VF capsules have been measured by Cao et al (2020). The numerical results and the experimental measurements are compared in Fig. 13 (a), Fig. 11 (b), and Fig. 13 (b) for 5%, 10%, and 15% VF capsules, respectively, which show good agreements between the two methods. In the first healing cycle, the healing efficiencies ψ of the foams embedded with 10% and 15% VF capsules were similar and higher than that of the foam embedded with 5% VF capsules. In the second healing cycle, the healing efficiency of the foam embedded 10% VF capsules was smaller than that of the foam embedded with 15% VF of capsules; the healing efficiency of the foam embedded with 5% VF capsules fell to the threshold demonstrating that there was small healing effect for the second healing cycle. These results suggest that higher volume fraction of capsules could bring better multiple self-healing performance. The numerically predicted damage and saturation of the foams embedded with 5% VF bilayered capsules, 10% VF bilayered capsules, and 15% VF bilayered capsules are shown at selected strains in Fig. 14(a) and (b), Fig. 12 (a) and (b), and Fig. 14 (c) and (d), respectively. The presence of the lower volume fraction capsules, such as 5% VF bilayered capsules, led to more localised damage pattern owing to stress concentration. Higher volume fraction capsules were more efficient in delivery of the healing agents as the damage pattern was more spreaded, which provided more microchannels for transport of healing agents.

The normalised elastic modulus of the capsules E *
To evaluate the effects of relative elastic modulus E * to self-healing efficiency, the FE calculations were conducted for the foam samples with E * = 0.5, E * = 0.119 and E * = 0.073. During the calculations, we kept the elastic modulus of capsules at a fixed value and the volume fraction of the capsules was fixed at 10%. Hence, the higher value of E * is related to the softer foam matrix. The selfhealing efficiencies obtained by FE simulations in the first and second healing cycles are shown in Fig. 15 (a), which suggests lower E * could provide better self-healing effect. To understand the mechanism, the numerically predicted damage and saturation of the foams with E * = 0.5, E * = 0.119, and E * = 0.073 are presented in Fig. 16 (a) and (b), Fig. 12 (a) and (b), and Fig. 16 (c) and (d), respectively. Higher value of E * , say E * = 0.5, led to more spreaded but lower level damage within the matrix foam owing to stress concentration around capsules. On the contrary, lower value of E * , say E * = 0.073, led to higher level damage that concentrated around the capsules. The healing system was more efficient to foams with lower value of E * owing to the facts that (1) capsules subjected to greater deformation which enabled more heal agents to be released from the capsules; and (2) the saturation level of healing agents in the concentrated damage areas was higher than that of the foam sample with higher value of E * in the spread areas of damage.

The normalised averaged diameter of the capsules d *
To evaluate the effects of the normalised averaged diameter of the capsules d * , the FE calculations were conducted for the foam samples with d * = 0.04, d * = 0.07, and d * = 0.16. In these calculations, the edge length of the foam samples was fixed at l f = 50 mm and the volume fraction of the capsules was fixed at 10%. The self-healing efficiencies obtained by FE simulations in the first and second healing cycles are shown in Fig. 15 (b), which suggests smaller capsules may have better self-healing efficiency at either the first healing cycle (d * = 0.07) or the second healing cycle (d * = 0.04). To understand the mechanism, the numerically predicted damage and saturation of the foams with d * = 0.16, d * = 0.07, and d * = 0.04 are presented in Fig. 17 (a) and (b), Fig. 12 (a) and (b), and Fig. 17 (c) and (d), respectively. The bigger capsules, say d * = 0.16, led to higher level localised damage within the samples. On the contrary, smaller capsules, say d * = 0.04, led to lower level, spreaded damage. The smaller capsules made the healing system more efficient as it was easier for the healing agents to flow to the smaller damaged areas around capsules to heal the damage.

Concluding remarks
This study provides the framework in modelling capsule based self-healing polymer foams under cyclic compression. In conjunction with micromechanical models that link the damage variable D to permeability k as well as the saturation s to the capillary pressure of the healing agents p w , the coupled pore fluid diffusion and stress simulations can capture the main features included in the self-healing system. Parametric studies based on 2D plane strain finite element simulation were conducted to understand the effects of the key design parameters of the healing system, which suggested the healing efficiency of the system could be improved if • the system had higher volume fraction of embedded capsules; • the ratio of the elastic modulus of capsules to the elastic modulus of foam matrix decreased; • smaller capsules were chosen. Further research will employ three-dimensional FE simulations to improve the accuracy of the numerical simulations.

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.

B. Finite element modelling of the bilayered capsules
The bilayered capsules were modelled as elasto-J2-plastic solids. To obtain the mechanical properties, quasi-static single capsule compressive tests were conducted. The effective elastic modulus E e of the capsules can be estimated based on the elastic Hertzian contact between the bilayered capsule and stainless plates during compression, i.e. where F c and δ c are compressive force and displacement measured from the Instron Universal Testing Machine and a linear variable differential transformer (LVDT), respectively; d c the diameter of a capsules. The elastic modulus E c of the capsules can be calculated as where μ c is the Poisson's ratio of the capsules, μ c = 0.1 (Nguyen et al, 2009) The finite element simulations on the single circular capsule compression test were conducted in order to calibrate the material parameters. Four-node plane strain elements (CPE4P in ABAQUS notation) were employed in the simulations. The capsule domain consisted of 20 elements, which could achieve converged results. Curve fitting via finite element simulation on the single capsule compression test indicated that the compressive strength and elastic modulus were 2.09 MPa and 45.04 MPa, respectively. The comparison between experimental measurements and numerical simulations is plotted in Fig. S1.
C. Calibration of the parameter of s re in the damage recovery model s n re represents the saturation which enables the complete recovery of damage at the nth compression cycle. To calibrate the values of s n re , 2D Finite element calculations on a foam sample containing 10% volume fraction capsules were conducted in comparison with experimental measurement. The details of the FE model are described in Section 6.2. Figure S2 shows the normalised elastic modulus E n ,E n = E he /E vi , at beginning of each compression cycle obtained by numerical simulations at selected values of s n re . The experimental measurement is included for comparison. The comparison suggests that s 1 re =0.175 and s 2 re =0.375 give the best fits for the first compression cycle and the second compression cycle, respectively.