Abstract
Being motivated by the technological applications of bioabsorbable polymeric materials in the fields of biomechanics and medicine, this paper presents a simple but efficient extension of Lemaitre's elastoplastic damage model by incorporating a chemicalbased (hydrolysis) degradation term. The aim is to allow the simulation of devices subjected to both mechanical and chemical environments. The model applicability is tested by a set of numerical finiteelement examples. The encouraging results show expected adequate coupling between the elastoplastic and chemical damages. Although the model is presently restricted to linear kinematics, the basic idea can be extended to finite strains.
Bioabsorbable polymers; elastoplastic damage; hydrolytic degradation
Simple extension of Lemaitre's elastoplastic damage model to account for hydrolytic degradation
Eduardo A. Fancello^{I,II,}^{*} * Author email: fancello@grante.ufsc.br ; Leandro P. Lindenmeyer^{I}; Carlos R. M. Roesler^{I}; Gean V. Salmoria^{III,II}
^{I}GRANTEGroup of Analysis and Mech. Design; Mech. Engineering Department
^{II}LEBmLab. of Biomechanical. Eng.; University Hospital
^{III}CIMJECTPlastic Molding Lab.; Mechanical Engineering Department Federal Univ. of Santa Catarina, Brazil
ABSTRACT
Being motivated by the technological applications of bioabsorbable polymeric materials in the fields of biomechanics and medicine, this paper presents a simple but efficient extension of Lemaitre's elastoplastic damage model by incorporating a chemicalbased (hydrolysis) degradation term. The aim is to allow the simulation of devices subjected to both mechanical and chemical environments. The model applicability is tested by a set of numerical finiteelement examples. The encouraging results show expected adequate coupling between the elastoplastic and chemical damages. Although the model is presently restricted to linear kinematics, the basic idea can be extended to finite strains.
Keywords: Bioabsorbable polymers, elastoplastic damage, hydrolytic degradation
1 INTRODUCTION
Medical implants made of bioabsorbable polymeric materials have been comprehensively investigated as efficient substitutes for their metallic counterparts. Costs, mechanical behavior, and their ability of being absorbed by the body once their mechanical function has ceased are the main features that motivate the use of these polymers. Among the most used bioabsorbable polymers developed for surgical procedures, the most common are polyglycolic acid (PGA) and polylactic acid (PLA) copolymers.
The process of degradation of bioabsorbable implants is represented by the irreversible phenomenon called hydrolysis that comprises the progressive breakage of the chemical bonds of the polymeric net due to its contact with bodily fluids. During this process, an initial decrease in the molecular weight is observed, followed by a loss in mechanical performance; finally, a loss of mass is observed. Due to this sequential degradation, as the mechanical properties of the bioabsorbable devices worsen, stresses are progressively transferred to the healed new tissue formed around them. In this way, they are designed to abandon their mechanical functions before being fully absorbed. This final absorption is the last step in the degradation process, which is found in biochemical reactions.
The hydrolytic degradation of the polymer depends on the intrinsic properties of the material, and it is related to the diffusion coefficient of water, degradation rate, body size, degree of crys tallinity, and mechanical stress state.
Degradation can occur in the following two ways (Chen et al., 2011). In the polymer medium, if the diffusion of water is slower than the degradation, rate, then surface degradation is observed, which progresses from the surface to the core. Conversely, when the water diffusivity allows for a rapid soaking of the entire device, then bulk degradation is observed, which acts on all the material points almost simultaneously. In the context of the present work, we assume that the latter functions as the primary agent of degradation for the PGA and PLA polymers.
Several studies have been performed to investigate the behavior of bioabsorbable polymers considering the effects of hydrolysis and mechanical stresses. Not aiming to present a systematic review, we cite some examples below.
An early study by Miller and Williams (1984) revealed that the virgin material of PLA polymers used in surgical sutures was subjected to strain values between 25% and 50% of their rupture strain values. The stress relaxation after 1421 days was evaluated, revealing that the rate of hydrolysis clearly increased with the magnitude of the applied strain. Chu (1985) reported a relationship between the applied strain and the degradation rate of PLA fibers immersed in a solution. When prestress was applied to the material in an aqueous medium, microvoids appeared, which increased the water penetration and accelerated the degradation process. Zhong et al. (1993) analyzed the influence of a strain value of 4% on the degradation parameters when using PGA and PLA suture fiberssoaked in a solution of water and hydrogen peroxide. After two weeks of immersion, they concluded that the degradation rate was considerably higher in the presence of strain than its absence. Smit et al. (2008) evaluated the effect of axial compressive loading on the degradation of a spine cage made of 70/30 poly (l, dllactic acid) (PLDLLA). The results showed that the mechanical resistance of the material was dependent on the load, temperature, and humidity. Soares et al. (2008, 2009, 2010) presented constitutive models for bioabsorbable polymers that underwent degradation under applied strain. They proposed a scalar field that described the local state of degradation (isotropic damage). This field is related to an evolution law that depends on the state of the strain/stress. Muliana et al. (2009) and Muliana and Rajagopal (2011,2012) analyzed the different relationships between hydrolytic degradation, viscoelastic behavior, and water diffusion. Basic hypotheses suggest that the Young's modulus of the material decreases with degradation. In a very recent paper, Khan and ElSayed (2013) includes a damage type internal variable within a variational viscoelastic model to account for degradation phenomenon in stents. It is interesting to note that none of the above mentioned works consider the combined influence of the damages induced by mechanical and chemical sources. For example, it is possible to find this coupling in vascular stents, because these devices are subjected to severe mechanical action (plastic diametral expansion) during the implantation procedure.
The present study aims to test a simple extension of the classical Lemaitre's isotropic elasto plastic damage model to include a timedependent term that accounts for the dynamics of a chemical degradation processes. On the basis of the experimental observations, it can be shown that this term depends on the stress state. The model combines the effects of elastic and plastic strains, mechanical damage, and chemical (hydrolytic) damage. Although this model is restricted to infinitesimal kinematics, it can be extended to finite deformations.
Section 2 presents a short but selfcontained mathematical formulation of the model. The notation used in this section is obtained from the book by de Souza Neto et al. (2008) where certain details about the technical operations can be found.
Section 3 illustrates the behavior of the proposed model. To this end, a set of numerical tests are conducted to evaluate the single constitutive equation as well as its inclusion in finiteelement examples. It is worth mentioning, however, that the results do not represent any particular identified material, and the intention of these results is to demonstrate the capabilities of the proposed model. A final section is added to summarize the obtained results.
2 HYDROLYTICPLASTIC DAMAGE MODEL
Let us consider an isotropic damage variable that is split into two terms: the first is the classical ductile damage of the Lemaitre's model and the other is called hydrolytic damage. While the evolution of the former is directly related to the existence of plastic strains, the latter is dependent on time (chemical action) and stress state. Now, let us define the conventional external and internal state variables, namely, the total strain ε, plastic strain ε^{p}, equivalent plastic strain ^{p}, plastic damage d^{p}, and hydrolytic damage d^{h}. The elastic strain and total damage can be defined as
For simplicity, we assume that both plasticity and damage phenomena are isotropic in nature. In this case, the free energy can be additively decomposed into two contributions: the elastic potential y^{e}, depending on the elastic strain and scalar damage, and the plastic one ψ^{p}, depending on the accumulated plastic strain ^{p}only:
The elastic potential is related to the linear isotropic material and is dependent on the bulk modulus K and shear modulus G:
In (3), D is the isotropic elasticity tensor, ε_{v} = tr [ε^{e}], the volumetric strain, and = ε^{e}, the deviatoric elastic strain. Considering the linear constitutive relationship it is also possible to express this energy in terms of the stress: where p =tr [σ] is the hydrostatic pressure, s = σ  pI is the deviatoric stress, and q = = is the equivalent von Mises stress. Their relationships with the corresponding effective values are
On the basis of thermodynamical principles, the conjugated forces can be obtained from the derivative of ψ with respect to the state variables:
Y^{p} is the energy release rate due to an increment in the plastic damage. The new conjugate force Y^{h} results in exactly the same conjugate value of Y^{p}; however, this value is now related to the hydrolytic damage. Despite the fact that Y^{h} and Y^{p} share the same value by definition, both these variables are expressed separately.
The von Mises yield function is given by
where σ_{y}(k) is the hardening law.
To define the evolution rules of the internal variables, a common formalism can be used to assume the existence of a dissipative potential ψ which ensures the automatic satisfaction of the ClausiusPlanck inequality (nonnegative dissipation process) on the basis of its convexity properties. In the present case, the potential ψ is exactly the same as that of the Lemaitre's damage model along with a new term ψ^{h}:
In (9), we distinguish the yield function ϕ plastic potential ψ^{p}, and the (new) hydrolytic potential ψ^{h}. Since ψ ≠ ϕ, the flow model is nonassociative.
In our case, we use the simplified Lemaitre and Chaboche plastic potential (de Souza Neto et al., 2008):
where r and s denote the material constants. The Heaviside step function H is used to incorporate a plastic threshold for the accumulated plastic strain ^{p} below which no plastic damage rate is allowed.
The evolution of the internal variables ε^{p}, ^{p}, and d^{p} are intrinsically related to the plastic flow phenomenon and can be defined based on the dissipation potential ψ as follows:
The N, H, and V values defined in (11), (12), and(13), respectively, are associated with the direction of evolution of the internal variables, while the Lagrange multiplier rate accounts for their amplitudes.
In the present proposition, the hydrolytic damage is not dependent on the plastic flow, unless it is indirectly related through the stress state. Therefore, ^{h} is independent of . Moreover, the hydrolytic damage is a timedependent process driven by the kinetics of the chemical reactions. In chemical kinetics, a classical expression that accounts for the degradation rate is given below (de Paula & Atkins, 2011, pp.1328):
where k is a rate constant depending on the temperature T, activation energy E_{a}, gas constant R, and a preexponential factor A . Moreover, the exponent n is known as the order of the reaction. In the present case, since the temperature and activation energy are assumed to be constant, we assume that the value of k depends on the strain energy Y^{h}. Then, following the structure of (14), we define the evolution of d^{h} as follows:
The hydrolytic dissipative potential consistent with (15) is then expressed as
where l, m, g are material constants. The first simple expression arises when taking n = 1, i.e., a firstorder chemical reaction, which provides an evolution that is conceptually similar to that of Soares et al. (2009). An important value in (15) is g that guarantees the existence of the degradation rate even in the case of absence of stresses.
Although similar, the expressions in (13) and (15) behave significantly differently; in the case of the former, the damage rate increases when the damage itself increases, which is the opposite behavior of the latter. In (13), the damage evolution is a timeindependent process since it follows the plastic multiplier ; in (15), the hydrolytic damage is essentially dependent on the time scale of the chemical process.
2.1 Incremental updating
The time integration of the previous equations over the time step [t_{n}, t_{n}_{ + 1}] defines the corresponding incremental constitutive algorithm. Here, we use the standard full implicit elastic predictor plastic corrector technique, and the presentation is similar to the procedures shown in de Souza Neto et al. (2008) and Simo and Hughes (1998). We consider the set of state variables at time t_{n} (which is completely known) as well as the total strain ε_{n}_{ + 1} = ε_{n} + Δε at time t_{n}_{ + 1}. The increment of the plastic deformation and damage comes from (1115) and Δγ = Δ t:
In (20), the evolution law of (15) was substituted after considering n = 1. The computation of the elastic strain and corresponding effective stresses is given by
The consistence condition in (8) at t_{n }_{+ 1} takes the form
2.2 Elastic predictor (ϕ_{n + 1} < 0)
In this step, we assume that the strain increment is purely elastic, i.e., Δγ = 0. In this case, the socalled trial quantities can be defined as
Since Δγ = 0, the damage increment, if it exists, can be attributed only to the hydrolysis phenomenon. Therefore,
Isolating from (29), we get
which allows the computation of , , and as follows:
Finally, the yield function can be evaluated for the trial values as follows:
If the condition < 0 is satisfied, then the step is purely (damage) elastic (Δγ = 0) and the variables are updated using the corresponding trial values:
Otherwise, > 0; then, a plastic increment Δγ > 0 exists. In this case, the new updated variables are determined to satisfy the condition ϕ_{n + 1} = 0 , which corresponds to the plastic corrector step.
Plastic corrector (Δγ > 0)
In this step, the following set of equations is solved:
From (33), it is posible to write
which may be substituted in (35). Similarly, replacing (34) in (36), the system is reduced to the following pair of equations:
where
Isolating w_{n}_{ + 1} = 1  d_{n}_{ + 1} from (38) and introducing it in (39), a single nonlinear equation in Δγ can be obtained:
Note that this expression differs from the classical Lemaitre's model in the last term of (41). The solution of this single nonlinear equation is obtained by using the NewtonRaphson method. In de Souza Neto et al. (2008), a convenient initial value Δγ_{(0)} is proposed as follows:
Once Δγ is obtained, all the remaining variables are updated:
Finally, the above algorithm σ_{n}_{ + 1} = (ε_{n}_{ + 1}, ε_{n}, , d_{n}) of the incremental constitutive problem yields the derivative.
providing the consistent tangent elastoplastic tensor. Clearly, a closed expression of this derivative is obtainable, but not included in this paper.
3 NUMERICAL EXAMPLES
In this section, we present a set of examples aiming to demonstrate the behavior and applicability of the proposed model. The material parameters used, however, do not correspond to any specific material and are taken just as examples. It is worth mentioning that the parameters related to the hydrolytic damage are clearly dependent on the unit of time, since it is a phenomenon that is significantly influenced by the time scale in which the chemical phenomenon occurs. Here, since no identification process was considered, the time unit used will be generically called a "day." The values shown in Table 1 were used for all the further studied cases and are intended to resemble those of a polymeric material.
3.1 Uniaxial tensile test
3.1.1 Straincontrolled tensile test: Case 1
In this case, we consider a single material point subjected to an axial tensile test controlled by axial strain. Three loads are applied that are linearly varied from zero to the maximum strain values of 0.05, 0.03, and 0.01 during a time interval of [0,120] days. Figure 1 shows the curves of the von Mises stress versus time and damage versus time for all the load conditions, which can be used to compare the behavior of the proposed plastichydrolytic damage model (DEPH) with a purely plastic damage (DEP) model. As expected, the hydrolytic damage curves exhibit higher values in the cases with higher stress values. In the case of the maximum strain of 0.01, no plastic strain occurs and only hydrolytic damage is visible.
3.1.2 Straincontrolled tensile test: Case 2
In this example, a similar tensile test is conducted; however, we consider that the maximum strain is applied for a short interval of 1/24 day and then maintained constant for up to 120 days. Figure 2 shows the curves of the von Mises stress versus time and damage versus time for all the load conditions. Again, the higher the load, the higher is the initial damage rate. During the initial loading (1/24 day), the hydrolytic damage is imperceptible and only plastic damage is visible.
3.1.3 Straincontrolled tensile test: Case 3
The same maximum strains are linearly applied during the first 60 days and then linearly unloaded during the next 60 days. Figure 3 shows the stress component (σ_{x}) versus strain curve and damage versus time curve.
3.1.4 Stresscontrolled tensile test: Case 4
Similarly to the previous cases, in this case, the tensile test is controlled by the applied axial stresses with maximum values of 20 MPa and 26 MPa that are linearly applied during a short time step of 1/24 day and then maintained constant. Figure 4 shows the curves of strain versus time and damage versus time. Note that in both these instances, the applied loads are less than the plastic yield value of 50 MPa. However, due to hydrolytic damage, plasticity is finally achieved, thereafter driving the material point to collapse.
3.1.5 Stresscontrolled tensile test: Case 5
In this case, the stress values of 18 MPa and 20 MPa are linearly applied for 60 days and unloaded in the next 60 days. Figure 5 shows the curves of the von Mises stress versus time and damage versus time. Again, in both these instances, the applied loads are less than the plastic yield value of 50 MPa. For the load of 20 MPa, hydrolytic damage occurs before the unloading and collapse. For the load of 18 MPa, despite the increasing damage, unloading is faster and the material point is able to withstand the load.
3.2 Axisymmetric sample
3.2.1 Displacementcontrolled test: Case 6
Consider the axisymmetric sample shown in Figure 6. An axial load is applied by controlling the displacement of the top surface of the sample and imposing symmetry boundary conditions on the lower surface. Fournode Lagrangian finite elements (with 2D solid revolution) were used to simulate this case. The mesh is also shown in Figure 6.
The prescribed displacement increases linearly from zero to 0.15 mm within the time interval of [0,120] days. Figures 7 and 8 show the von Mises stress distribution and damage distribution for the two instances: 40 and 94 days, respectively. From Figure 7, it is evident that strong relaxation occurs in the instance when hydrolysis is taken into account. Further, the stress distribution changes significantly; the maximum value emanates from the notch and appears in the interior of the domain.
3.2.2 Forcecontrolled test: Case 7
The same sample as that used in Case 6 is subjected to a controlled axial force that grows linearly from zero to 1021 N within the time interval of [0,120] days. After that, progressive unloading is undertaken at the same rate, reaching (theoretically) zero value at 120 days. This maximum force value was set up so that the virgin material could reach the maximum stress value of 50 MPa. Figures 9 and 10 show the von Mises stress distribution and damage distribution for two different times: 60 and 79 days, respectively. A behavior similar to that in the previous case is evident; due to the presence of hydrolysis, the stresses at the notch relax and the maximum von Mises stress value is found to exist in the interior of the domain.
It is noteworthy that since this is a forcecontrolled test, the sample cannot withstand the load if hydrolysis is considered. Figure 11 shows the axial strain at the notch in the sample. The unloading at 60 days for both the instances is clearly visible; however, since hydrolysis keeps damaging the sample, plastic collapse is achieved at 79 days (Figure 15, loss of convergence).
3.2.3 Plate Sample: Case 8
In this case, we simulate a plate sample using the 8node Lagrangian brick elements (Figure 12), which can be used to represent a 3D sample discretized with solid elements. In this case, a distributed force is progressively applied on the top surface of the plate; the maximum value of 250.75 N is achieved at 60 days. After that, the unloading phase is initiated, reaching the zero load value at 120 days. Figures 13 and 14 show the distribution of the axial stress component u_{y} and the damage at two times: 60 and 120 days. As expected, the behavior is similar to that in
the cylindrical case; the damage due to hydrolysis significantly modifies the stress distribution at the notch. During the unloading phase, it is possible to visualize the competition between the remaining resistance of the specimen and the unloading tensile force. At the end of the process, the local compressive forces remain near the notch region for the hydrolytic model.
4 FINAL REMARKS
Being motivated by the technological applications of bioabsorbable polymeric materials, the objective of this paper is the proposal of a simple extension of the classical Lemaitre's elastoplastic damage model that includes the representation of damage induced by chemical phenomena (such as hydrolytic degradation). In a simple and efficient manner, this model demonstrates the relationship between the damage evolution due to mechanical and chemical actions. The applicability of the model is analyzed by undertaking a set of numerical tests. Although the model is presently restricted to linear kinematics, the idea can be extended to finite strains.
A set of specific technical remarks are highlighted below:
Hydrolytic damage interferes with all the state variables at all the material points subjected to chemical activity. Moreover, this damage appears even in the absence of stresses due to the material parameter g.
The evolution law for the hydrolytic damage in which the term (1  d)^{n}; n > 0 appears is essentially different than that of the mechanical damage. Because of the former, the damage rate decreases with increasing damage, behavior that is opposite to that given by the model of mechanical damage.
Mechanical damage is independent of time since it undergoes plastic evolution. Hydrolytic damage is essentially dependent on the time scale of the chemical process.
The results of the investigated examples are consistent with the expectations for different situations. At the beginning of the load, the evolution of the damage is driven only by the hydrolysis process. When the material reaches the plastic flow, the degradation rate increases significantly due to the coupled action between both the hydrolytic and plasticitydriven damages.
In the cases when the material point is subjected to a straincontrolled process, the expected stress relaxation depends on the degradation rate. On the other hand, when the material is subjected to a stresscontrolled process, the degradation forces the material to undergo plastic flow followed by collapse.
In all the investigated cases, hydrolytic damage increased faster than elastoplastic damage. This, however, should be attributed to the arbitrary choice of the material parameters in the examples.
The proposition in (2) regarding free energy decoupling, although conventional, cannot be considered a unique alternative. For example, in the studies of Brunig (2002, 2004) and the references therein, the damage follows a kinematic description similar to that of plastic strains and the free energy is additively decoupled in elastic (reversible) and inelastic (damage and plastic) contributions. This approach would certainly modify the definitions of the conjugate forces Y^{h} and Y^{p} in (7). In spite of these alternatives, the present work has aimed to maintain the same structure of the original Lemaitre's elastoplastic damage proposition (Lemaitre, 1985).
Summarily, it is worth mentioning that the actual representation of a bioabsorbable material deserves attention in several aspects that have not been considered here. Among others, one such aspect is the modification of the material properties due to the chemical actions: in the present simple damage model, the yield limit and elastic modulus of the virgin material decrease as long as the damage increases. The process of aging due to chemical actions may produce embrittlement and/or stiffening of the elastic behavior instead of the softening action observed in the present case.
Acknowledgments
The authors would like to thank the Brazilian agencies CNPq, CAPES and FINEP for the financial support offered for this research.
Received in 12 May 2013
In revised form 30 Jul 2013
 Y. Chen, S. Zhou, and Q. Li(2011). Mathematical modeling of degradation for bulkerosive polymers: Applications in tissue engineering scaffolds and drug delivery systems. ActaBiomaterialia, 7:11401149.
 C.C. Chu (1985).Strainaccelerated hydrolytic degradation of synthetic absorbable sutures. Surgical research, recent developments. Proceedings of the First Annual Scientific Session of the Academy of Surgical Research, 111115.
 M. Brunig (2002). Numerical analysis and elasticplastic deformation behavior of anisotropically damaged solids. International Journal of Plasticity, 18:12371270.
 M. Brunig (2004) An isotropic continuum damage model: Theory, and numerical analyses. Latin American Journal of Solid and Structures, 1: 185218.
 E. A. de Souza Neto, D. Peric, and D. R. J. Owen (2008). Computational methods for plasticity: Theory and applications. WILEY.
 K.A. Khan, T. ElSayed (2013). A Phenomenological constitutive model for the nonlinear viscoelastic responses of biodegradable polymers. Acta Mechanica 224(2):287305.
 J. Lemaitre (1985). Coupled Elastoplasticity and damage constitutive equations. Com. Meth. Appl. Mech. Engng, 51:3149.
 N.D. Miller and D.F. Williams (1984). The in vivo and in vitro degradation of poly(glycolicacid) suture material as a function of applied strain. Biomaterials, 5:365368.
 A. Muliana, K.R. Rajagopal, and S.C. Subramanian (2009). Degradation of an elastic composite cylinder due to the diffusion of a fluid. Journal of Composite Materials, 43(11 ):12251249.
 A.H. Muliana and K.R. Rajagopal (2011).Changes in the response of viscoelastic solids to changes in their internal structure. Acta Mechanica, 217(34):297316.
 A. Muliana and K.R. Rajagopal (2012). Modeling the response of nonlinear viscoelastic biodegradable polymeric stents. International Journal of Solids and Structures, 49(78):9891000.
 J. V. Simo and T. J. R. Hughes (1998). Computational Inelasticity. Springer.
 T. H.Smit, T A P Engels, P I J M Wuisman, and L E Govaert (2008). Timedependent mechanical strength of 70/30 poly(L,DLLactide) shedding light on the premature failure of degradable spinal cages. Spine, 33:1418.
 J. S. Soares, J. E. Moore, JR., and K. R. Rajagopal (2008).Constitutive framework for biodegradable polymers with applications to biodegradable stents. ASAIO Journal, 295301.
 J. S.Soares, K. R. Rajagopal, and J. E. Moore Jr. (2010).Deformationinduced hydrolysis of a degradable polymeric cylindrical annulus. Biomechanics and Modeling in Mechanobiology, 9:177186.
 S. P. Zhong, P. J. Doherty, and D. F. Williams (1993). The effect of applied strain on the degradation of absorbable suture in vitro. Clinical Materials, 14:183189.

*
Author email:
Publication Dates

Publication in this collection
03 Feb 2014 
Date of issue
Oct 2014
History

Accepted
30 July 2013 
Received
12 May 2013