## Abstract

Self-healing smart materials have emerged into the research arena and have been deployed in industrial and biomedical applications, in which the modelling techniques and predicting schemes are crucial for designers to optimize these smart materials. In practice, plastic deformation is coupled with damage and healing in these systems, which necessitates a coupled formulation for characterization. The thermodynamics of inelastic deformation, damage and healing processes are incorporated here to establish the coupled constitutive equations for healing materials. This thermodynamic consistent formulation provides the designers with the ability to predict the irregular inelastic deformation of glassy polymers and damage and healing patterns for a highly anisotropic self-healing system. Moreover, the lack of a physically consistent method to measure and calibrate the healing process in the literature is addressed here. Within the continuum damage mechanics (CDM) framework, the physics of damage and healing processes is used to introduce the healing effect into the CDM concept and a set of two new anisotropic damage–healing variables are derived. These novel damage–healing variables together with the proposed thermodynamic consistent coupled theory constitute a well-structured method for accurately predicting the degradation and healing mechanisms in material systems. The inelastic and damage response for a shape memory polymer-based self-healing system is captured herein. While the healing experimental results are limited in the literature, the proposed theory provides the mathematical competency to capture the most nonlinear responses.

## 1. Introduction

A novel category of smart materials with the ability to heal micro- and macro-scale damages have been developed recently. These smart materials may use a microencapsulated healing agent inside their matrix, in which micro- or macro-damages fracture the wall of the microcapsules and the released healing agent solidifies in the presence of catalysts and heals the microcracks and/or voids (White *et al.* 2001). Toohey *et al.* (2007) extended this idea by designing a microvascular network, within the body of the system, to feed the healing agent continuously into the material. However, this system suffers from clogging of the vascular network after the first round of healing, which limits the repeatability of the healing process. In particular, the performance of these systems in healing macro-scale damages is still an obstacle that must be overcome. This kind of healing methodology is herein referred to as a coupled damage–healing system to indicate that both the damage and the healing processes are activated concurrently in the system. The healing systems using external triggering to activate the healing mechanism are herein referred to as a decoupled damage–healing system. In these types of systems, the damage and the healing processes occur separately, where just one of these processes is active. These systems may use thermoplastic particles (TPs; Zako & Takano 1999), thermally reversible covalent bonds (Liu & Chen 2007; Plaisted & Nemat-Nasser 2007), confined shape recovery of shape memory polymers (SMPs; Li & John 2008) or biomimic two-step close-then-heal (CTH) mechanisms to molecularly heal structural length-scale damages (Li & Nettles 2010; Li & Uppu 2010; Nji & Li 2010; Xu & Li 2010).

Extensive demand for using this new branch of smart materials in biomedical and industrial applications necessitates the development of physically consistent constitutive equations to accurately predict their inelastic, damage and healing responses. For inelastic simulation, one may use the basic thermodynamic principles to develop the inelastic constitutive equations. Viscoplasticity theories developed in the thermodynamic framework can effectively describe the metal-like inelastic responses (Iwakuma & Nemat-Nasser 1984; Nemat-Nasser & Obata 1986; Chaboche 1991, 2008). Inelastic responses of glassy polymers are quite irregular as the elastic region is followed by yielding and softening regions, which are representations of intermolecular resistances against rotation of molecular chain segments in a polymeric network. That is followed by the strain hardening region, which results from network resistance against alignment along the principal direction of plastic deformations (Argon 1973; Boyce *et al.* 1988*a*). This irregular response has been modelled by many researchers using phenomenological approaches (Boyce *et al.* 1988*a*; Khan & Zhang 2001; Dupaix & Boyce 2007).

The inelastic deformation and damaging phenomena have been treated as a dissipative thermodynamic process in the literature, while here these are considered to be partially reversible processes where healing governing equations establish that reversibility. The dissipative nature of plasticity and damage has been investigated by a number of authors (Lemaitre 1985; Yazdani & Schreyer 1990; Hansen & Schreyer 1994; Bryant *et al.* 2008; Naderi *et al.* 2010). The modelling of the healing process has been investigated by Adam (1999) and Barbero *et al.* (2005), and Voyiadjis *et al.* (2011), who used a thermodynamic framework to investigate the healing process. In this respect, a uniaxial damage–healing variable is proposed (Voyiadjis *et al.* in press *a*).

In this work, within the thermodynamic framework, a novel potential function is introduced to simulate the irregular inelastic deformation of glassy polymers. This thermodynamic consistent formulation is extended to capture the damage and healing responses in polymer-based self-healing systems. These new potential functions provide mathematical competency to simulate irregular inelastic response of polymers within the thermodynamic framework. The lack of physically consistent measurement variables for self-healing systems is addressed in this paper, and two new anisotropic damage–healing tensors are established with the ability to capture the multi-axial case of anisotropic damage and healing. These variables are developed based on the physical description of the damage and healing processes within the continuum damage mechanics (CDM) framework. While recently developed self-healing systems are capable of healing microscale as well as macroscale damages (Li & John 2008; Nji & Li 2010), owing to the fact that the established damage–healing variables have been developed within the CDM framework, the length scale of the damages is restricted within the microscale in this work in order to preserve the continuity assumption of the CDM framework. Consequently, one may find that these variables are applicable to the microscale damaging and healing processes, where the state of material is considered to be a continuum. The performance of the proposed damage–healing variables and coupled constitutive equations are examined for an SMP-based self-healing system. The proposed theory constructs the bases of a physically consistent approach towards modelling the coupled inelastic, damage and healing processes in self-healing systems. The introduced measurement technique for healing is generalized to measure the healing based on changes in the elasticity modulus. This physically practical method is used in characterizing the healing process, such as identifying the required amount of healing agent, defining the required programming procedure or even identifying chemical compositions to obtain the desirable recovery stress or strain and healing efficiency. The thermodynamic framework is used in §2 to derive the coupled constitutive equations. In §3, two new damage–healing variables are derived using the CDM framework. Section 4 contains experimental observations regarding a macroscale healing process, and §5 represents the computational aspect and simulation results.

## 2. Thermodynamics of coupled viscoplasticity, damage and healing

The thermodynamic framework provides a physically consistent description for different mechanisms involved in materials deformation, such as elasto-viscoplastic deformations and damage–healing mechanisms. It is common to use phenomenological approaches to formulate the irregular inelastic deformation of glassy polymers (Khan & Zhang 2001; Krempl & Khan 2003); however, in this work, the thermodynamic framework is incorporated to describe different mechanisms associated with the irregular inelastic deformation of glassy polymers. A detailed discussion on these internal mechanism changes in polymeric networks during their inelastic deformations is given by Boyce and co-workers (Boyce *et al.* 1988*a*,*b*, 1989). Accordingly, observable and internal variables are used to describe each process in the thermodynamic framework. The following internal variables are assumed to define the thermodynamic equilibrium state: strain tensor *ϵ*_{ij}, plastic strain tensor , plastic kinematic hardening variable tensor *α*_{ij}, plastic isotropic hardening variable *p*, damage tensor *d*_{ij}, damage kinematic hardening variable tensor , damage isotropic hardening variable *d*^{I}, healing tensor *h*_{ij}, healing kinematic hardening/softening variable tensor and healing isotropic hardening/softening variable *h*^{I}. This formulation results in the most generalized coupled constitutive equations, which incorporate all possible phenomena for each process such as kinematic and isotropic hardening/softening. To derive constitutive equations for inelastic, damage and healing processes, one may use the Legendre transformation of the internal energy, *u*, in order to obtain the Helmholtz free energy (HFE) function as: *Ψ*=*u*−*Ts*, where *T* is the absolute temperature and *s* is the entropy. In general, the HFE is a function of mentioned thermodynamic fluxes and is expressed as
2.1

As shown in equation (2.1), it is assumed that isotropic hardening fluxes are scalars (no distortion in yield surfaces). Kinematic and isotropic hardening variables show the shift in the centre and the change in the size of their yield surfaces, respectively. The damage tensor, *d*_{ij}, and healing tensor, *h*_{ij}, are captured by the damage–healing variables that are defined in §3. Other fluxes such as *d*^{K}_{ij}, *d*^{I}, and *h*^{I} provide enough flexibility to map the simulation of the damage and healing on their respective experimental results. The conjugate thermodynamic forces for each flux are then obtained as
2.2where *ρ* is the material density and *σ*_{ij} is the Cauchy stress tensor, *X*_{ij} is the deviatoric backstress tensor, which is defined as , and *δ*_{ij} is the Kronecker delta. Scalar *y*^{pI} is the isotropic hardening conjugate force, and are, respectively, second-order tensors of conjugate forces for the damage and healing. In equation (2.2), and are second-order tensors of the damage and healing kinematic hardening conjugate forces, respectively. However, *y*^{dI} and *y*^{hI} are, respectively, the damage and healing isotropic hardening conjugate forces (scalars). The scalar isotropic hardening conjugate forces, *y*^{dI} and *y*^{hI}, represent uniform changes in size of the damage and healing yield surfaces. The HFE is then decomposed as follows:
2.3where the first term is named the elastic HFE, while hardening function ℵ captures the hardening/softening effects. The relationship between the damaged–healed fourth rank elastic modulus tensor, *E*_{ijkl}, based on the virgin elastic modulus, , and the damage–healing variables is derived in §3. A generalized quadratic form of the hardening function, ℵ, is used here,
2.4where *c*_{i} (*i*=1–6) are material constants. In order to derive the flow rules associated with plasticity, damage and healing processes, the power of dissipation along with the HFE definition is used to obtain the following expression:
2.5Upon substituting the thermodynamic conjugate forces from equation (2.2) into equation (2.5), one may obtain the most general form of the dissipative power function, including the kinematic and isotropic hardening terms for each process, as follows:
2.6In order to accomplish the derivation of the constitutive equations, the associative yield surfaces for plasticity, damage and healing should be incorporated. The classical von Mises yield surface, *f*^{p}, is used for the plastic yield surface while damage, *f*^{d}, and healing, *f*^{h}, yield surfaces, which capture the kinematic and isotropic hardening effects in both damage and healing processes, are defined as
2.7where *s*_{ij} is the deviatoric stress tensor, which is defined as *s*_{ij}=*σ*_{ij}−(1/3)*σ*_{kk}*δ*_{ij}. Tensors and contain material constants in order to capture the growth of the damage and healing surfaces as well as the initiation of damage and healing and they are defined as follows (Voyiadjis *et al.* 2011):
2.8where *λ* is the Lamé constant, *ν* and *ν*′ are the damage and healing thresholds, respectively, and *ζ* and *ζ*′ are, respectively, the damage- and healing-related material constants. The irregular inelastic response of polymeric materials and the nonlinear nature of the damage and healing processes force the constitutive equations to be non-associative. In the case of inelastic response of polymers, the associative flow rules, which are widely used in simulation of viscoplastic responses of metals (Chaboche 2008), do not have mathematical competency to simulate the softening and subsequent hardening responses. This deficit was compensated by phenomenological approaches in the literature (Khan & Zhang 2001). The thermodynamic framework is used here to simulate the irregular inelastic response of polymers through a non-associative formulation. This idea is extended for the damage and healing processes to obtain the most generalized form of constitutive equations for these processes. The basic idea of this generalized thermodynamic-based theory was extracted from the well-known Chaboche viscoplasticity theory, in which the back stress evolution laws are incorporated to obtain the metallic inelastic responses and viscoelasticity responses of polymer-based material systems (Chaboche 1997). Here, in order to obtain the mathematical competency to capture different hardening and softening zones associated with the viscoplastic deformation of glassy polymers, the following potential functions are proposed, which provide an extensive flexibility to obtain different nonlinear responses for each phenomenon:
2.9where *a*_{(r)}, *b*_{(r)}, *c*_{(r)}, *Λ*_{(r)}, *z*_{(r)} and *β*_{(r)} are viscoplasticity-related material con- stants. Adding prime and double prime to all of these material constants identifies the damage- and healing-related material constants, respectively. Furthermore, *k*_{1(r)},*k*_{2(r)} and *k*_{3(r)} are viscoplasticity, damage and healing material constants, respectively. Tensors and *y*^{hK,s}_{ij(r)} show the saturation values of their respective conjugate tensors, and . Subscript *r* has the values of 1,2,…,*n*,*m*, and *p* in which *n*,*m* and *p* indicate, respectively, the total number of viscoplasticity, damage and healing material constants. The ‘|.|’ shows effective values and it is defined as for the second rank tensor (other effective values are obtained in the same way). As mentioned before, equation (2.9)_{1} provides the mathematical competency to simulate the irregular viscoplastic response of glassy polymers within the thermodynamic framework. This task is accomplished by spectral decomposition of the total strain and simulates each spectrum with a set of material constants. The physical basis of this spectral decomposition relies on distinctive resistance mechanisms which are active at each strain level. For example, intermolecular resistance at initial yielding produces a softening response, while network resistance at higher strains results in a strain hardening response (Argon 1973; Boyce *et al.* 1988*a*). Consequently, one may assume that a different thermodynamic process governs each spectrum and model it with a specific potential function. Taking into account the history of loading and variable changes ensure the continuity of the evaluated internal variables. In such a way, a thermodynamic consistent model that captures the irregular inelastic response of glassy polymers is derived. This decomposition approach appears to be quite precise in simulating irregular inelastic deformation of glassy polymers, as shown in §4. The material constants for each process need to be identified through matching the simulations with their respective experimental results. A detailed discussion on the curve fitting process to obtain the viscoplasticity-related material constants can be found in Voyiadjis & Abu Al-Rub (2003). Owing to the fact that the inelastic responses in glassy polymers are quite nonlinear, one may find that conventional curve fitting techniques, such as least-square methods (Voyiadjis & Abu Al-Rub 2003), are not applicable here. Then other strategies may be adopted for finding the material parameters; for example, in each strain spectrum, the material constants can be changed in a systematic way until the gap between the simulation and the experimental data passes a predefined limit (Voyiadjis *et al.* in press *b*). In the case of damage and healing processes, the same technique may be adopted to find the material constants, while in some cases with regular rate dependencies the least-square methods, described by Voyiadjis & Abu Al-Rub (2003) in detail, are applicable. In many engineering applications, the graphical curve fitting techniques result in enough accuracy and they are used in this work to find the material constants. It is shown in §5 that only seven material constants are required to capture the inelastic deformation in each strain spectrum zone, while the damage potential function requires only five material constants to precisely simulate the damage response. Finally, the coupled constitutive equations are obtained by minimizing the following functional:
2.10Applying the necessary derivatives, ∂*γ**/∂*σ*_{ij}=0, and to minimize *γ**, results in the following generalized coupled constitutive equations for a multi-axial case of loading:
2.11To simplify the derivation of kinematic and isotropic evolution laws, one may assume that hardening/softening mechanisms for inelasticity, and damage and healing processes are uncoupled from each other; the evolution laws for these internal variables are obtained as follows:
2.12The Lagrangian multipliers are obtained by applying consistency conditions: for plasticity, ; for damage, ; and for healing, . The relation between the damage and healing conjugate forces and the Cauchy stress tensor *σ*_{ij} at the damaged state is given by Voyiadjis *et al.* (2011),
2.13where *Q*_{ijkl} is a damage–healing transformation tensor introduced in §3 (equation (3.8)). Equation (2.13) renders the state of damage and healing explicitly as a function of Cauchy stress and damage and healing variables. These equations provide coupling between damage, healing and inelastic processes.

## 3. Anisotropic damage and healing variables

Kachanov (1958) made a pioneering contribution in describing the material degradation based on damage density and calibrating the area reduction owing to microscale damages. This concept was later generalized to a multi-axial anisotropic case by Murakami (1988). This approach is depicted in figure 1, where figure 1*a* depicts the real state of damaged material that is assumed to be decomposed to a fictitious effective configuration in figure 1*b*, which carries the load, and a fictitious fully damaged configuration, figure 1*c*, which cannot sustain load. Murakami (1988) defined the damage variable tensor, *ϕ*_{ij}, to represent the transformation between the real damaged area vector, *dAn*_{i} (figure 1*a*), and the effective fictitious area vector, (figure 1*b*), as follows:
3.1

The pure damaged area vector, (figure 1*c*), is then obtained as a function of the damage variable tensor, *ϕ*_{ij}, and the area vector, *dAn*_{j}, in the real damaged configuration. Equation (3.1) is rearranged to obtain
3.2

The physics of the healing process is used here to find a variable to calibrate the healing mechanism. During the healing process, some of the microscale damages are healed and then the effective area, which carries the load, is increased. To represent this phenomenon, it is assumed that the fictitious fully damaged configuration (figure 1*c*) undergoes the healing process. This process for an anisotropic healing case is shown in figure 2, where figure 2*a* shows the pure damaged state without load capacity, while figure 2*b* shows the healed configuration and figure 2*c*,*d* shows the fictitious effective healed configuration and remaining damages after accomplishing the healing process, respectively.

According to the presented physical description of the healing process, a new second rank anisotropic healing variable tensor, *h*_{ij}, is proposed as follows:
3.3where *h*_{ij} captures the transformation between the real damaged area vector, *dAn*_{i} (figure 3*a*), and the fictitious healed area vector, (figure 3*b*). Figure 3 shows the overall mapping procedure between the real damaged and fictitious healed effective configurations.

Owing to the complex mapping techniques required in direct measurement of damage and healing variables based on defect density, one may use an indirect measurement method such as calibrating damage and healing based on elastic modulus changes during these processes (Lemaitre & Dufailly 1987). A fourth-order anisotropic damage variable tensor, *κ*_{ijkl}, to capture the elastic modulus change is defined as follows:
3.4where superscripts in equations (2.1) and (2.2) indicate the two different mathematical tensorial expressions for the damage tensor when normalizing it with reference to the inverse of the undamaged elasticity tensor, , and is the damaged elastic modulus. In order to measure the healing indirectly, a new fourth rank healing variable tensor, *h*′_{ijkl}, is introduced to measure the elastic modulus changes after accomplishing the healing process as follows:
3.5where is the elastic modulus of the healed material in which represents no healing case (where 0_{ijkl} is the fourth rank zero tensor). One may assume that the maximum healing is obtained when all produced damages are cured, and shows that the elastic modulus is fully recovered. Substituting from equation (3.4) into equation (3.5) results in the following expression for the healed elastic modulus :
3.6

As shown in figure 3, two kinds of stresses are dealt with in CDM, the Cauchy stress tensor, *σ*_{ij}, which is the stress state in the real damaged configuration (figure 3*a*), and the effective stress tensor, , which is a representation of the stress state in the effective healed configuration (figure 3*b*). In general, two approaches are available to derive the relation between these two stress tensors. One may equate the strains in both configurations, , and the second method is the so-called equivalent elastic energy (EEE) method, which is more reliable. In general, the EEE method relies on energy equivalence between the real and effective configurations that incorporate more physical phenomena than a simple equivalence of strains. The EEE approach has been extensively used in the literature to derive the relations between effective and real configurations (Lemaitre 1985; Yazdani & Schreyer 1990; Hansen & Schreyer 1994). Making use of the second approach, which is stated between real damaged and fictitious healed configurations, results in
3.7

A fourth rank damage–healing transformation tensor, *Q*_{ijkl}, is defined to represent the transformation between the Cauchy stress tensor, *σ*_{ij}, in the real state and the Cauchy stress tensor , in the effective configuration. It is expressed by the following relationship:
3.8

Substituting from equation (3.8) and from equation (3.6) into equation (3.7) provides the relationships between *h*′_{ijkl}, and *Q*_{ijkl} as follows:
3.9Stating equilibrium between the real and effective configurations results in the relation between stresses in these two states, as shown in the following equation:
3.10where *I*_{ijkl} is the fourth-order identity tensor. The fourth rank damage effect tensor, *M*_{ijkl}, and the healing effect tensor, *H*_{ijkl}, are defined in the following. One may find that *H*_{ijkl}=0_{ijkl} shows the unhealed configuration in which . And *H*_{ijkl}=*I*_{ijkl} represents the fully healed state where . Consequently, the following tensorial representations for *M*_{ijkl} and *H*_{ijkl} are proposed:
3.11

Comparing equation (3.10) and equation (3.8), one may find
3.12The relationships between the second rank damage variable tensor, *ϕ*_{ij}, and the healing variable tensor, *h*_{ij}, with the fourth rank damage tensor, *d*_{ijkl}, and the healing variable tensor, *h*′_{ijkl}, are obtained by substituting *Q*_{ijkl} from equation (3.12) into equation (3.9),
3.13Upon substituting *M*_{ijkl} and *H*_{ijkl} from equation (3.11) into equation (3.13), the relationships between all introduced damage–healing variables are established. Consequently, by knowing one set of these damage–healing variables, the other set can be obtained using equation (3.13).

## 4. Experimental observations

The performance of the proposed coupled constitutive equations is examined for uniaxial compression on a thermoset polystyrene SMP (Veriflex, CRG Industries)-based self-healing system containing TPs (Abifor 501; Abifor Company; Nji & Li 2010). This study followed a bioinspired two-step CTH self-healing system proposed by Li & Uppu (2010). It uses the confined shape recovery of the SMP matrix for the purpose of closing macroscopic cracks and the dispersed TPs for the purpose of healing molecularly. In order to make the SMP matrix smart, a three-step thermomechanical programming is conducted, including compression at a temperature above the glass transition temperature (step 1), cooling while maintaining the compressive strain constant (step 2) and removal of the applied stress at a temperature well below the glass transition temperature (step 3). The compression test and sample preparation follow the American Society for Testing and Materials (ASTM) standard (ASTM-C365 2005). Figure 4 shows the three-step strain-controlled compression programming for the pure SMP. Figure 4 also has shows step 4, the free shape recovery (Li & Nettles 2010). It is shown that this system heals structural length-scale damage (crack) repeatedly, efficiently and molecularly (). The SMP sample is programmed under the aforementioned programming steps before the structural-level damage is induced into the sample. Figure 5*a* shows the schematic of a damaged configuration at room temperature in which a three-point single edge notch bend test is used to introduce a macroscale crack into the system. The spheres indicate the dispersed TPs in the SMP matrix. Figure 5*b* shows the heated configuration, from room temperature to a temperature above the SMP glass transition temperature, *T*_{g}, where the shape memory property of the SMP is activated and, under confined boundary conditions, the crack is sealed. Figure 5*c* illustrates the state of the sample when the temperature reaches above the melting point of TPs, in which the molten TP molecules diffuse into the crack surfaces and fill the gap. Figure 5*d* shows the state of the SMP after the cooling down process where the crack is fully healed in this configuration. A schematic of the confined shape recovery of the SMP is depicted in figure 5*e*, where the configuration of the rigid rod and sample is shown inside the rigid cylindrical fixture (Li & Uppu 2010). The confined recovery starts by heating the set-up to a temperature above *T*_{g}, while the rigid fixture and rod retain the overall shape of the sample. Consequently, the activated SM property, with external confinement by the rod and the rigid cylinder, forces the sample to expand and fill the open internal space and the cracks within the sample are closed. As shown in figure 6, it is experimentally confirmed that the TPs in this step diffuse into the SMP matrix and provide the desired molecular-level healing (). Experimental observations confirm the efficiency of this healing mechanism in which molecular entanglements are gained after healing and the mechanical properties of the self-healing system are fully recovered (Nji & Li 2010; Voyiadjis *et al.* 2011). Figure 6 shows scanning electron microscopy (SEM) and transmission electron microscopy (TE) results for a macroscale crack, produced by a three-point bending test (ASTM-D5045 2007). Figure 6*a* shows an SEM image of a macroscale crack in the self-healing system, where the opening of the crack, crack path and crack tip are highlighted on the SEM image. Figure 6*b* shows the healed configuration where the crack path has almost disappeared, and the TEM image shows the interface between the TPs and the SMP matrix after diffusion.

For the mentioned self-healing SMP embedded with TPs, the elastic modulus of the SMP sample during a compression loading at room temperature increases from an initial 680 to 845 MPa upon loading and unloading in different strain ranges. One may consider these changes in the elasticity tensor as an indication of the state of the damage inside the material system (Lemaitre & Dufailly 1987). The same strategy may be adopted to calibrate the state of the healing in the self-healing material systems in which the elasticity tensor after accomplishing the healing process can be compared with the damaged elasticity tensor. The differences between the elasticity tensors in the healed and damaged configurations illustrate the magnitude of the healing process. This experimental procedure for calibration of the healing process in the self-healing materials is discussed in detail by Voyiadjis *et al.* (2011).

## 5. Computational aspect and simulation results

Generally speaking, the computational aspects of the coupled problem can be divided into two distinctive modules. One of them ensures that the basic solid mechanics governing equations are satisfied, while boundary conditions and geometry of structure are enforced. Basically, these are equilibrium, strain–stress and compatibility relations that result in an initial boundary value problem (IBVP). The computed state of inelastic strain, damage and healing variables is iteratively corrected (Newton–Raphson technique) to satisfy this IBVP along with the initial and boundary conditions (Mendelson & Manson 1957; Voyiadjis *et al.* in press *b*). The second computational module deals with the flow rules and nonlinear governing equations of inelastic strain, damage and healing variables. An incremental solution, which provides history dependency in a coupled inelastic–damage–healing problem, updates the internal variables. These updated values are used during the iteration process of the IBVP solution and this task continues until the IBVP solution is converged. Different solution algorithms may be incorporated to solve the nonlinear governing equations of inelasticity, damage and healing, including iterative return-mapping (Simo & Taylor 1986) or non-iterative methods (Sivakumar & Voyiadjis 1997). A comprehensive description of solution algorithms for coupled as well as uncoupled inelastic–damage–healing problems can be found in Voyiadjis *et al.* (in press *b*). The return-mapping technique proposed by Mendelson & Manson (1957) is generalized here for all the three processes.

To compute the damage inside the materials, the proposed viscoplasticity theory is used to simulate the half-cycle strain-controlled compression test of the SMP where damage and plastic deformation are coupled. Consequently, the Cauchy stress, *σ*_{ij}, and the stress rate, , in the real damaged state can be introduced into the damage computational module from the viscoplasticity solution. Once the damage criterion is satisfied, following each inelastic strain increment, a damage increment is computed while the updated values of stress and stress rate are incorporated. The computed damage values are then used to update the elastic modulus for the next load increment. This computational procedure allows full coupling between inelastic and damage computation (Voyiadjis *et al.* in press *b*).

In a coupled self-healing system, the healing computational module receives the updated inelastic deformation and damages from the inelastic and damage computational modules. In the case of a coupled healing system, after each increment of the inelastic deformation and damage, the healing computation module computes the amount of the healing. However, in an uncoupled self-healing system, after a certain amount of damage, all final inelastic deformations and damage variables are introduced into the healing computational module and the healing is computed based on these values. Finally, the updated elastic modulus is introduced to the viscoplasticity module for the next load step of computations. This computational procedure allows each mechanism to follow its respective governing equations while full coupling is considered.

Figure 7*a* represents changes in the elastic modulus of the SMP sample under compression tests with a strain rate of 0.002 (s^{−1}) in which the proposed damage theory is used to capture this behaviour. Figure 7*b* shows the compression test results of the SMP at room temperature with a strain rate of 0.002 (s^{−1}) and the simulation is obtained using the proposed plasticity theory. Tables 1 and 2 give, respectively, the required material constants for damage and inelastic deformation of the SMP sample. In order to show the capability of the proposed theory to simulate a vast range of irregular inelastic deformations of glassy polymers, additional experimental data, regarding the inelastic mechanical response of polyethylene terephthalate (PET), are presented here. Table 3 shows the material constants to simulate PET tension at room temperature; the simulation results are shown in figure 8 and the experimental results are reported by G'Sell *et al.* (2002). The inflection in experimental data points shows the respective softening and subsequent strain hardening region for PET, which is captured by the proposed model. In these tables, indicates the yield strain and and are, respectively, the final values of *X*^{s} at previous strain spectrum (Voyiadjis *et al.* in press *b*).

It is worthwhile indicating that precise inelastic simulations of material systems are crucial in the damage and healing computations where the stress state in the real damaged and/or healed configurations is used to compute the state of the damage and/or healing.

## 6. Concluding remarks

The mechanics of damage and healing processes is investigated for self-healing systems, and the bases of a physically consistent modelling of these phenomena are constituted. The lack of healing variables for calibrating the healing process is compensated by introducing two new anisotropic damage–healing variables that are defined based on the physics of the damage and healing processes. Thermodynamics of the plasticity, damage and healing are investigated, and coupled constitutive equations are developed in the most generalized forms. Kinematic and isotropic hardening effects for all processes are incorporated into the formulation, together with a new set of potential functions in order to achieve the mathematical competency in precisely simulating each phenomenon. The irregular inelastic deformation of the SMP and PET is captured using the proposed viscoplasticity theory. The developed viscoplasticity theory is based on a strain spectrum strategy where each strain spectrum is captured using a set of material constants. These strain spectrums are functions of strain rate and temperature and they should be evaluated experimentally. The performance of the proposed damage–healing variables and coupled constitutive equations is examined for an SMP-based self-healing system. The proposed anisotropic damage–healing variables provide the designers with the ability to measure the state of the damage and healing in a highly anisotropic self-healing system, and the constitutive equations perform quite well in simulating the plastic and damage responses of glassy polymers. As shown in this work, the proposed potential functions are capable of capturing the irregular inelastic deformation and the damage responses of the polymer-based material systems. However, availability of the experiments for the self-healing systems is limited in the literature; it is expected that the healing process shows a similar trend but inverse to that of the damaging process in which the proposed theory provides the mathematical competency to capture it. The healing test set-up is under preparation and will be presented in a forthcoming paper together with the corresponding healing experiments. This will be reported together with the corresponding simulations, using the proposed theory. It is worth noting that the proposed theory considers a continuum scale without incorporating the microstructural changes during each process. It is well known that such a macroscopic model, although it may provide good correlation with a specific loading condition, may not be capable of capturing a general case of loading in which loading histories deviate from the experiments used to fit the material parameters. This opens a new field to incorporate such microstructural changes into the governing equations of these coupled processes. The microstructural physics of the problem can be incorporated in CDM framework by introducing *fabric tensors*. Fabric tensors capture the crack and void distributions in the damaged material and they have been linked to the CDM concept (Voyiadjis & Kattan 2006, 2007*a*,*b*). This will be the work in a forthcoming paper in associating fabric tensors and micromechanics with the processes of both damage and healing.

## Acknowledgements

This study was sponsored by NSF under grant nos CMMI0900064 and HRD0932300. The financial support by NSF was greatly appreciated. G.Z.V. acknowledges the collaboration with Prof. Taehyo Park of the Hanyang University, Seoul, Korea, under the World Class University project funded by the National Research Foundation of Korea.

- Received May 21, 2011.
- Accepted August 1, 2011.

- This journal is © 2011 The Royal Society