- Email: [email protected]
PEEK biocomposite under quasi-static tensile load
Materials Science and Engineering A 382 (2004) 341–350
Modeling of the mechanical behavior of HA/PEEK biocomposite under quasi-static tensile load J.P. Fan, C.P. Tsui∗ , C.Y. Tang Department of Industrial and Systems Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong, PR China Received 9 January 2004; received in revised form 29 April 2004
Abstract A cylindrical three-phase unit-cell model was established for predicting mechanical properties of hydroxyapatite (HA)-reinforced polyetheretherketone (PEEK) biocomposite. The model consists of an elastic–brittle HA spherical particle, an elasto-plastic matrix, and a very thin interphase region between the particle and the matrix. Perfect bonding was initially assumed among three phases before occurrence of particle–matrix debonding. For simulation of the particle–matrix debonding process of HA/PEEK, a damage evolution equation was incorporated into the material properties of the interphase region. A ductile damage evolution in the matrix was also taken into account, and was treated as a criterion for the failure of the biocomposite. A user subroutine for the damage-coupled material properties for the matrix and the interphase region was built and incorporated into a finite element code named ABAQUS. The influence of interfacial adhesion on the non-linear stress–strain relation of HA/PEEK has also been simulated by varying the tensile strength of the interphase region. The strength of the interphase region in the composite could be estimated by a comparison between the predicted and experimental results. By using the unit-cell modeling technique incorporated with the stress-based debonding and the composite failure criterion and selection of appropriate interphase strength, the numerical simulation was found to be close to the experimental results from available literatures. In vivo performance of HA/PEEK as implant materials was also discussed. © 2004 Elsevier B.V. All rights reserved. Keywords: Unit-cell model; HA/PEEK biocomposite; Finite element method (FEM); Interface debonding; Matrix damage; Particle volume fraction (PVF)
1. Introduction Metallic materials, such as stainless steel and cobalt– chromium alloys, have been successfully used as bone implants in orthopedic areas for many years. Recently, the newly developed polymers have become feasible for replacing these conventional materials due to their advantageous properties and high performance [1–4]. Among these polymers, polyetheretherketone (PEEK) is a good choice because of its excellent mechanical properties. Owing to its favorable structural and load-bearing functions, PEEK has gained much attention in the area of aerospace and marine industries [5]. In addition, it possesses exceptionally good chemical behavior, fatigue resistance, high-temperature durability, and excellent wear properties. From the processing viewpoint, PEEK can be easily formed using conventional plas-
∗
Corresponding author. Tel.: +852 2766 6635; fax: +852 2362 5267. E-mail address: [email protected] (C.P. Tsui).
0921-5093/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.msea.2004.04.078
tic processing methods, such as compression and injection molding, facilitating a rapid and economical production of implant devices [6,7]. Additional benefits of PEEK including its capacity to be repeatedly sterilized and shaped readily by machining and heat contouring render its suitability to become a biomaterial [8]. However, stiffness of PEEK cannot satisfy the requirement as a load-bearing implant material for replacing human bone. If the stiffness of the implant material such as the metallic one is too high, the newly generated living tissue will experience a stress-shielding problem [9] so that fracture healing can hardly function normally. On the other hand, if the stiffness of the implant material like some polymers is too low, the undesirable callus will form at the fracture site [10]. Hence, there is a need to develop an implant material that matches the elastic modulus of bone in order to minimize the stress-shielding effect as well as formation of the callus during the process of fracture healing. Therefore, incorporation of a second-phase hard material into the polymer matrix becomes inevitable for enhancing mechanical properties of the resulting polymer compos-
342
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
ites. Since the bone tissue of mammals has been found to contain 58% of calcium hydroxyapatite (HA), great efforts have been made to develop phosphate ceramics as a potential implant material [11]. Artificially synthesized HA has proven to be a successful substitute for the natural one. This material is not only biocompatible and non-toxic, but also exhibits unique osteoinductive properties [12]. HA/PEEK has been conceived as a biomaterial for bone replacement on the basis of producing appropriate mechanical compatibility as well as necessary biocompatibility. Besides considering whether the material stiffness can satisfy the requirement for load-bearing applications, the ultimate tensile strength of HA/PEEK must be adequate to prevent the material from catastrophic fracture. The implant made of biocomposite should undergo a sufficient geometrical deformation to show its ductile behavior. Some experimental observations from other investigators [12] have revealed that ductility of the biocomposites is greatly reduced by increasing HA particles. In the meantime, the ultimate tensile strength was found to decrease due to the presence of a poor interfacial adhesive bond. It is so because the increase in particle volume fraction (PVF) results in an increase in the area of debonded interface, which in turn causes a weakening of the ultimate tensile strength. Therefore, for HA/PEEK and all other potential composites utilizing HA particles, prediction of their mechanical properties using certain technique shows to be crucial for designing materials with certain desirable properties. Elastic properties and failure processes of HA particles-reinforced biocomposites have been investigated [1,2]. However, due to constraint of the linear elastic analysis employed in those investigations, their predictions for ductile failure processes may not attain a high accuracy. Hence, the elasto-plastic large deformation analysis is required for predicting the material properties of biocomposites. Recently, experiments have been performed by Bakar and his co-workers [3,4,13,14] to investigate the mechanical properties of HA/PEEK biocomposite for possible usage as bone replacement. Micro-structural analysis was carried out in their investigations to correlate the structure–property relationship of the composite. The dependency of the tensile properties, such as elastic modulus, strength, strain to fracture, and micro-hardness, on PVF of the HA particles were also presented. However, the failure mechanism of this biocomposite could not be completely revealed in their investigations. Thus, theoretical and numerical techniques are required for prediction of material properties. In this paper, a three-phase unit-cell model incorporating a large deformation and elasto-plastic constitutive equation has been built for predicting the mechanical properties of HA particles-reinforced biocomposite using finite element method (FEM). Interfacial debonding around the interphase region was considered in the formulation of the model by using the vanishing finite element technique [15]. Damage evolution in the matrix was also taken into account for serving as a composite failure criterion. The failure mode
transformation would be revealed from the damage distribution. The model was also applied to predict the non-linear stress–strain relation of the HA/PEEK biocomposites under the influence of different strength values of the interphase region. The optimum interphase strength would be estimated by a comparison between the predicted and experimental results. The model with the selected interphase strength, the debonding criterion, and the composite failure criterion would be further applied to predict the tensile strength and strain to fracture of HA/PEEK.
2. Elasto-plastic constitutive equations The HA particles in this study, which are brittle and with high modulus, could be accurately characterized by a pure elastic constitutive law. The polymer matrix was assumed to undergo elasto-plastic deformation with isotropic hardening, while the interphase region was assumed to have the same material properties as the matrix except for the damage evolution equation and tensile strength. Moreover, elastic strains of the HA particles were assumed to be very small with respect to unity, whereas their rotations could be arbitrarily large. The behavior of the HA particles can then be reasonably described by the following hypoelastic law e Jij = Cijkl Dkl
(1)
where σ J is the Jaumann derivative of the Cauchy stress tensor, D denotes the strain rate tensor. Ce expresses the fourth-order elastic constants tensor, which is represented by the Lame’s constants λ and µ in an isotropic material as e Cijkl = λδij δkl + µ(δik δjl + δil δjk )
(2)
where δ is the Kronecker delta. The elastic strains were assumed to be small in the matrix and the interphase region. Thus, the constitutive law used for the matrix and the interphase region is similar to Eq. (1), but the elastic strain rate is replaced by subtraction of plastic strain rate from the total strain rate, namely, p
e Jij = Cijkl (Dkl − Dkl )
(3)
where Dp is the plastic strain rate tensor. Assuming that the flow rule is derived from a plastic potential, it is easily found that the plastic strain rate Dp is coaxial with the deviator σ of the Cauchy stress tensor. P Dkl =
3p˙ σ 2σeq kl
(4)
where p˙ is the equivalent plastic strain rate and σeq is the von Mises equivalent stress. They can be expressed as p p p p p p p˙ = 2/3[(˙ε11 − ε˙ 22 )2 + (˙ε22 − ε˙ 33 )2 + (˙ε33 − ε˙ 11 )2 p2
p2
p2
+ 6(˙ε12 + ε˙ 23 + ε˙ 31 )]1/2
(5)
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
and
σeq
343
Substituting Eq. (10) into Eq. (9) gives [(σ11 − σ22 )2 + (σ22 − σ33 )2 + (σ33 2 + τ 2 + τ 2 )]1/2 − σ11 )2 + 6(τ12 23 31 = √ 2
eq
(6)
This constitutive law for the matrix and the interphase region can be represented by an one-dimensional relation between ε˙ p and σeq . In this study, the power hardening law is used as n ε σ = σ0 (7) for ε > ε0 ε0 where n denotes the hardening index, σ0 is the initial yield stress and ε0 is the corresponding strain such that ε0 = σ0 /E.
3. General description of dynamic explicit formulation Finite element programs based on an explicit dynamic formulation have been proven to be a very powerful tool in industrial application because of its numerical efficiency and other advantages, such as low memory requirement, no iteration convergent problem, easy treatment of contact condition and availability of the vanishing element treatment, as compared with implicit methods. The general form of the principle of virtual work is described as follows [16] σij δεij d + ρu¨ i δui d − ρfi δui d ρgi δui dA = 0 (8) − Af
where u and u¨ represent the displacement and accelerate vectors, r denotes mass density, f and g are the body and surface force vectors. δu and δε denote the virtual displacement and strain, respectively. Ω represents the volume of the body under study, and Af denotes the surface which is subject to an external applied force. The left-hand side of the above equation involves the terms for internal work, inertia work, work done by the body force, and work exerted by the surface force. Meshing the body with finite elements and introducing material behavior model, element shape function, and dynamics of rigid body, a finite element equation of motion is obtained and expressed in the following form at the time step n Mij u¨ nj = Rni − Fin
Mij un+1 = eq Fin j
(11)
where eq M and eq F denote the equivalent mass and equivalent applied force matrix, respectively. If eq M is lumped as a diagonal matrix, matrix inversion is not needed and a solution can be obtained directly by solving a linear equation. The lumping scheme is computationally economical because the matrix inversion consumes large computing time. Often in dynamic analysis, use of the lumping mass achieves more accurate results than that of consistent mass [17]. The central difference method has a selective convergence according to the magnitude of t, and the accuracy and convergence are linearly proportional to the square of t [17]. After the time increment of every element, the magnitude of global time increment can be determined by using the following equation: tn+1 = α t,
where
t = min{ t1 , t2 , . . . , tN } (12)
where N is the total element number and ti is the time increment of the ith element. The safety constant, a, is often selected to be less than 0.9. The critical time increment is determined such that tc = Ls /C, where Ls is the characteristic length, which is the given element area divided by the largest edge. Moreover, √the propagation speed C is determined such that C = Et /ρ, where Et is the tangent modulus and r represents the material density.
4. Ductile damage model for the matrix Under excessive external applied load, micro-voids or micro-cracks may appear in the polymer matrix. These micro-defects, which lead to material degradation, are non-reversible processes. In the present study, a damage variable, D, was used to characterize the degradation of the material’s load-bearing capability. Fig. 1 shows the damaged body of a representative volume element (RVE). A0 is the nominal intersection area and Aeff is the effective sectional
(9)
where M is the mass matrix, Rn is the overall external applied force vector and F n is the nodal force vector. In order to obtain the solution at the time step n + 1 when all of the variables at the time step n are known, the central difference method for the time discretization of the acceleration is introduced: u¨ n =
u˙ n+1/2 − u˙ n−1/2 un+1 − 2un + un−1 = t 2 t 2
(10)
Fig. 1. Representative volume element for a damaged material.
344
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
area without the presence of micro-defects and their mutual interaction. Thus, a damage variable, D, is expressed as D=1−
Aeff A0
(13)
According to this definition, the damage variable is between 0 and 1. Generally, D = 0 represents the case for an undamaged material and D = 1 denotes the case for material failure, but for most instances, the critical damage value cannot exceed a certain critical value Dcr = 1 [18]. In this study, a potential function ϕ for a ductile damage process in the polymer matrix is taken as follows [19] S −Y 2 p˙ ϕ= (14) α/n 2 S D p2/n
Fig. 2. Morphology of flame spherical HA particles.
where α and S are material-dependent constants. The damage ˙ can be obtained as evolution D ∂ϕ Y p˙ ˙ = D = (15) α/n ∂(−Y) S D p2/n The strain energy density release rate is in the form as Y=
2 R σeq v
2E(1 − D)2
(16)
where the triaxiality function Rv is described by Rv = 2 2 3 (1 + ν) + 3(1 − 2ν)(σm /σeq ) . Thus, the damage evolution ˙ D may be calculated for any loading process defined by the histories of the accumulated plastic strain rate p, ˙ the von Mises equivalent stress σeq and the hydrostatic stress σm . The damage-coupled material behaviors of the elasto-plastic polymer matrix were incorporated into ABAQUS/Explicit finite element program through a user-defined material subroutine: VUMAT sub-module.
5. Unit-cell modeling 5.1. Formation of a unit-cell model Bakar et al. [4] showed in their thermal spraying experiment that the HA particle is spherical in shape with a smooth surface as shown in Fig. 2. Although these spherical particles could hardly be evenly dispersed in a polymer composite, the micro-structure of the composite was assumed to be periodical in their investigation so that it could be represented by repeated three-dimensional hexagonal cells. As the use of three-dimensional hexagonal cells involves a larger amount of calculating time as compared with that of two-dimensional cells, a two-dimensional axisymmetric cylindrical cell [20], which has been extensively used to analyze the mechanical behaviors of particulate-reinforced composites [21,22], was used in the present study. The simplification process of the unit-cell model is illustrated as shown in Fig. 3. In spite of the gap between the exact cell and its approximated one, the simulation using the cylindri-
Fig. 3. Simplification process for the unit-cell model: (a) hexagonal array; (b) hexagonal cell model; and (c) cylindrical cell model.
cal cell only results in a relatively low error with a limited influence on the overall response of the composite [23]. Each representative cell contains a hard particle situated at its centre, which is surrounded by the polymer matrix cylinder. Due to the assumed axisymmetry for the packing of the spherical particles, only a half of the representative cell is analyzed, as shown in Fig. 4. The half height of the cylindrical cell is set equal to its radius. The volume fraction, Vp , of the particle is related to its radius. In this analysis, the constituent material properties, the particle shape, and the interphase thickness were kept fixed while only PVF was varied up to 40%. The material properties of each constituent phase are summarized in Table 1 [1,4].
Fig. 4. Geometry of the axisymmetric cylindrical cell model.
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350 Table 1 Material properties of HA and PEEK Materials
E(Gpa)
ν
ρ (kg/m3 )
PEEK HA
3.2 85
0.42 0.3
1291 3160
A third-phase material, the interphase region with a thickness of 1% of the cell radius, was introduced between the particle and the matrix [24]. Perfect bonding was initially assumed between the particle and the interphase region, and between the interphase region and the polymer matrix before occurrence of particle–matrix debonding. When damage and failure took place within the interphase region, interfacial debonding would occur and result in a degradation of the strength and stiffness of the composite. The cell model was then converted into a finite element cell model constructed by axisymmetric four-node linear elements. The finite element mesh of the cell model containing 10 vol.% particles is shown in Fig. 5, which clearly illustrates that the model is composed of the particle, the interphase region and the polymer matrix. 5.2. Boundary conditions The influence of adjacent material deformation was taken into account by using boundary conditions. By symmetry, there exist following displacement boundary conditions: uz = 0
on
z=0
(17)
ur = 0
on
r=0
(18)
where uz and ur represent the axial and radial displacement, respectively. To ensure compatibility among all periodic representative cells, the line at r = R0 of the cylindrical cell is required to remain straight and parallel to their initial state after deformation. Thus, the normal force acting on the edge at r = R0 must be set to zero such that σr dA = 0 on r = R0 (19) A1
where A1 is the hoop area of the cylinder, and all shearing stresses on all boundaries were set to zero.
Fig. 5. Finite element cell model with 10% particle volume fraction.
345
In order to satisfy this requirement, an arbitrary point P, which does not belong to any part of the model and is allowed to freely move in any directions, was created and related to all nodes on the line r = R0 . Thus, displacement of all element nodes on the line r = R0 in the r-direction are the same and equal to that of the point P. (ur )r=R0 = (ur )P
(20)
With this method, no external force is applied onto this edge. Hence, the boundary conditions set in Eq. (19) is satisfied. The velocity boundary condition is forced on the upper edge at z = H0 . In order to minimize the inertia effect, the loading velocity, v, was taken as 1 m/s.The overall uniaxial true stress, sth , and logarithm strain response, elog , of the composite are σth =
1 Tz (1 + ε) and A
εlog = ln (1 + ε)
(21)
where ε denotes the nominal strain in the axial direction and Tz is the traction applied in the same direction. A stands for current area of the cell normal to the axial direction During small-range deformation, no interfacial debonding and matrix plastic deformation take place. The overall deformation remains within the linear elastic range. Thus, the Young’s modulus of the composite, Ec , can be estimated by Ec =
(σz )z=H0 (εz )ave
(22)
where (σz )z=H0 = Tz /A is the nominal stress, and (εz )ave = (uz )z=H0 /H0 is the nominal strain. In small deformation, the nominal stress and strain are more or less the same as that of the true ones. 5.3. Particle–matrix debonding criterion It was revealed, from the tensile fracture surface of HA/PEEK samples [4], that extensive interfacial debonding between the HA particle and the polymer had occurred (Fig. 6). Thus, the composite was supposed to undergo interfacial debonding during deformation at certain critical stress. For the sake of simulating the particle–matrix debonding process, a very thin interphase region with a thickness of 1% of the cell radius was incorporated between the particle and the matrix. The interphase was initially assumed to have the same material properties as the polymer matrix except for the damage evolution and tensile strength. Perfect bonding was initially assumed between the particle and interphase and between the interphase and the polymer matrix before occurrence of the interfacial debonding. When radial strain in the interphase region exceeds a critical value, damage would be developed. The damage evolution, Di , in the interphase region is defined as Di =
σn σc
(23)
346
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
Fig. 6. Tensile fracture surface of HA/PEEK biocomposite.
where σ n and σ c are radial and critical stresses, respectively. The radial stress is defined as 1 1 σn = (σr + σz ) + (σr − σz ) cos(2ϕ) + τrz sin(2ϕ) (24) 2 2 where angle ϕ is shown in Fig. 4. During the damage process of the interphase region, the stiffness will degrade. When the radial stress σ c reaches its critical value, the failure was assumed to occur in the interphase region. When all of the material points in a certain element fail, all stress components in this element would be set to be zero. As this element looses its load-bearing capability, it would be removed from the mesh through the vanishing finite element technique [15] for simulation of the interfacial debonding process. All of these procedures have been accomplished in ABAQUS/Explicit analysis through a user-defined material subroutine: VUMAT sub-module. 6. Results and discussions 6.1. Young’s modulus In order to verify the proposed numerical model, Chow’s analytical solution [25] for the elastic modulus of a composite system with ellipsoidal particles based on Eshelby approach is used for comparison. For spherical particle, the Young’s modulus is reduced to (Kp /Km − 1) 2(Gp /Gm − 1) Ec = Em 1 + Vp + Vp 3B 3A (25) where Vp is PVF. K and G are bulk and shear moduli. The subscripts p and m represent the particle and the matrix, respectively. The material parameters are given by Kp A=1+ − 1 (1 − Vp )α, Km Gp B =1+ − 1 (1 − Vp )β (26) Gm
Fig. 7. Young’s modulus of HA/PEEK vs. volume percentage of the HA particles.
in which 1 1 + νm α= 3 1 − νm
and
β=
2 15
4 − 5νm 1 − νm
(27)
where νm is Poisson’s ratio of the matrix. The dependency of the predicted Young’s modulus of HA/PEEK on the HA content is plotted in Fig. 7. The predicted Young’s modulus increases with an increase in HA content. This means that the elastic modulus of the composite has been enhanced by the addition of HA particles. The results generated from both the numerical prediction and the analytical Chow’s formula are close to the experimental data [13,14]. In addition to the Chow’s formula, our proposed model could provide a simulation of stress–strain distribution and damage process in the composite. It can be clearly observed that HA/PEEK with 40% PVF of HA has Young’s modulus 3.5 times greater than the unfilled PEEK. Young’s modulus of cortical bone tissue is about 7–30 GPa [4]. Therefore, the stiffness of HA/PEEK with 40% PVF is comparable to that of the cortical bone. 6.2. Stress–strain relation The interaction between the HA particle and the PEEK matrix is not mechanical interlocking only. Experimental observation of the fracture surface of HA/PEEK [26] showed that most HA particles have fibrous PEEK residue adhering to the surface of the particle. Hence, there are some forms of interfacial adhesion between the particle and the matrix. In order to study the influence of the interfacial adhesion on the mechanical behaviors of the composite like HA/PEEK, tensile strength of the interphase region was varied between 40% and 80% of that of the polymer matrix. These composites were supposed to undergo interfacial debonding as usual. In addition, two extreme cases were also taken into account in this study for comparison. One of them is that the interphase region was assumed to have the same tensile strength as that of the matrix and to have no interfacial
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
347
Fig. 8. Stress–strain relation of HA/PEEK with 10% PVF and different interphase strength.
Fig. 9. Tensile strength of HA/PEEK vs. volume percentage of the HA particles.
debonding during deformation. The other extreme case is that the composite was assumed to have no interfacial bonding or strength so that all filled particles have been fully debonded before loading and thus contain only voids equal to the volume fraction of the fully debonded particles. Fig. 8 shows the effects of different interphase strength on the predicted stress–strain relation of HA/PEEK with 10% PVF. Within elastic deformation range, there is no distinguished effect of different interphase strengths on the stiffness of the composites. During plastic deformation, the stress of the composite with stronger interphase strength becomes significantly higher than that with weaker interphase strength. Hence, a proper choice of coupling agent is necessary for enhancing the interfacial strength and adhesion as well as the stiffness and strength of the composite. In addition, the non-linear stress–strain relation of the composite with the interphase region equal to 40% of the matrix strength is in good agreement with the experimental data [4]. Thus, the interphase strength was set to be about 40% of the matrix strength in this study.
shown in Fig. 9 because of the interfacial debonding. The numerically predicted results are in good agreement with the experimental data [4] and also close to the results of our previous model with an interphase having a different damage evolution equation [29]. During the interfacial debonding process, particles usually do not fully debond from the polymer matrix. Thus, there is a certain stress transfer from the matrix to the particle; the particle can still help to carry part of overall loading. As most of the loading becomes principally carried by the remaining matrix, the tensile strength of the composite with partially debonded particles is reduced. This explains why the numerically predicted results are much higher than that of the empirical model, which assumed no interfacial adhesion.
6.3. The ultimate tensile strength In most cases, the ultimate tensile strength could be enhanced by addition of particles, especially in metallic-matrix composite [27]. Unfortunately, this advantage fails to happen in biocomposite due to the poor interfacial bond. A simple empirical model that assumed no interfacial adhesion taking place between the matrix and the particle is expressed as [28] 2/3
σc = σm (1 − 1.21Vp )
(28)
where sc and sm are the tensile strength of the composite and the matrix, respectively. The predicted ultimate tensile strength decreases with increasing PVF and falls below that of the unfilled PEEK as
6.4. The strain to fracture In addition to consideration of the interfacial debonding in this study, the matrix damage was also taken into account. The critical damage value Dcr was taken to be 0.6. This value was obtained by fitting the experimental results. When the damage in any elements of the matrix reaches this critical value, catastrophic failure was assumed to happen in the composite. The logarithm strain determined using Eq. (21) is defined as the strain to fracture. Fig. 10 shows that there is an inverse relation between the strain to fracture and PVF. The strain to fracture was found to drop drastically before 10% PVF and then goes down gradually. Thus, the addition of the HA particles results in reduction of the ductility of the composite. Moreover, brittle failure dominates after 20% PVF. The predicted results are close to the experimental data [4,14]. The deformed finite element meshes at the moment of composite failure are shown in Fig. 11. Among the composites under study, only the composite with 10% PVF experiences a large deformation, while the composites with a
348
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
6.5. In vivo performance of HA/PEEK as implant materials
Fig. 10. Strain to fracture of HA/PEEK vs. volume percentage of the HA particles.
higher PVF experiences brittle fracture. The damage fields for 10% and 40% PVF are shown in Fig. 12. The damage path, along which the failure takes place, in the composite with 40% PVF is much shorter than that with 10% PVF. This explains why the composite with 40% PVF has undergone a brittle fracture.
Particulate-reinforced bicomposites like HA/PEEK, HAPEX, and HA/PLLA have been extensively developed for orthopedic applications, particularly as load-bearing implant materials [13]. PEEK is one of the high-performance polymer matrices suitable for biomedical applications because of its good chemical and fatigue resistance, high-temperature durability, good wear properties, and ability to be repeatedly sterilized before implantation [4]. An in vivo investigation has been investigated by Abu Bakar et al. [13] to study the biological response of the HA/PEEK implant materials inserted into pigs up to 24 months. Effects of varying porosity, pore size distribution, and volume content of HA particles in the implant materials on the biological response were also studied in their investigations. During early stages of implantation, fibroblast cells that promote vascularization were found. At longer implantation periods, osteroid and osteocytes were formed within lamellar bone to enable osteoblastic activites. Their in vivo investigations suggested that HA/PEEK had favorable bioactivity and biocompatibility suitable for load-bearing implant applications with long-term performance. Thus, our proposed model can be applied by material designers to design an implant
Fig. 11. Interfacial debonding at the moment of composite failure: (a) 10% PVF; (b) 20% PVF; (c) 30% PVF; and (d) 40% PVF.
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
349
Fig. 12. Damage distribution in the polymer matrix: (a) 10% PVF; and (b) 40% PVF.
material with desired mechanical properties and enough load-bearing capability.
7. Conclusions High-performance composites, such as particulatereinforced biocomposite like HA/PEEK, have got great development for load-bearing orthopedic applications. Prediction of their mechanical properties becomes important for material design and application. The cylindrical unit-cell model has been successfully applied to predict the mechanical properties of HA/PEEK biocomposites, including elastic modulus, ultimate tensile strength, stress–strain relation, and strain to fracture. The initiation of the particle–matrix debonding process has been taken into account on the basis of the stress-based criterion, which was only applied to the interphase region of the unit-cell model. Matrix damage evolution was also considered in this model to serve as a criterion for composite failure. A user-defined VUMAT material subroutine describing the damage-coupled elasto-plastic constitutive behavior of the intephase region and the matrix has been successfully established and incorporated into the ABAQUS finite element code for prediction of the mechanical properties. The mechanical behaviors of the particulate biocomposites have been found to be much influenced by the interfacial adhesion, which was successfully characterized by using different strengths of the interphase region in this study. Thus, the strength of the interphase region in the composite with 10% PVF could be estimated by comparison between the predicted and experimental results. The cell model incorporated with the proposed debonding, the failure criterion, and the selected interphase strength is found to be successful in predicting not only the trend of the mechanical properties of the HA/PEEK composites in relation with the HA content but also the non-linear constitutive relation and the ductile-to-brittle failure transition of HA/PEEK. Generally, the Young’s modulus increases with the increase in the HA content while the tensile strength and the strain
to fracture show a decrease with the increasing particle concentration. The predicted results are close to the experimental data from available literatures. This implies that the use of the stress-based debonding criterion as well as the composite failure criterion has been justified. Thus, the cylindrical unit-cell model is well suitable for prediction of the non-linear constitutive relation of the HA/PEEK composites. With this model, the interfacial debonding process and the damage path in the composite could be visualized for the material designer.
Acknowledgements The authors would like to thank the Research Grant Council of Hong Kong for funding support of this project (PolyU5176/00E) and Prof. Khor for providing us the SEM pictures of the HA/PEEK biocomposite.
References [1] F.J. Guild, W. Bonfield, Biomaterials 14 (1993) 985. [2] F.J. Guild, W. Bonfield, J. Mater. Sci.: Mater. Med. 9 (1998) 497. [3] M.S. Abu Bakar, P. Cheang, K.A. Khor, J. Mater. Process. Technol. 89–90 (1999) 462. [4] M.S. Abu Bakar, P. Cheang, K.A. Khor, Mater. Sci. Eng. A345 (2003) 55. [5] F.N. Cogswell, D.C. Leach, Plast. Rubber Compos. Process. Appl. 18 (1992) 249. [6] A. Lustiger, F.S. Uralil, G.M. Newaz, Polym. Compos. 11 (1990) 65. [7] J.R. Sarasua, P.M. Remiro, J. Mater. Sci. 30 (1995) 3501. [8] A.J. Barton, R.D. Sagers, W.G. Pitt, J. Biomed. Mater. Res. 30 (1996) 403. [9] J.V. Sloten, L. Labey, R.V. Audekercke, G.V. Perre, Biomaterials 19 (1998) 1455. [10] K.A. Jockisch, S.A. Brown, T.W. Bauer, K. Merritt, J. Biomed. Mater. Res. 26 (1992) 133. [11] G. Hienderauer, T. McGee, R. Kudej, Am. Ceram. Soc. Bull. 70 (1991) 1010. [12] G. Zambonin, M. Grano, Biomaterials 16 (1995) 397. [13] M.S. Abu Bakar, M.H.W. Cheng, S.M. Tang, S.C. Yu, K. Liao, C.T. Tan, K.A. Khor, P. Cheang, Biomaterials 24 (2003) 2245.
350
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
[14] M.S. Abu Bakar, P. Cheang, K.A. Khor, Compos. Sci. Tech. 63 (2003) 421. [15] G.C. Li, X.W. Ling, H. Shen, Int. J. Plast. 16 (2000) 39. [16] K.L. Bathe, Finite Element Procedure in Engineering Analysis, Prentice-Hall, New Jersey, 1982. [17] ABAQUS/Explicit: User’s Examples and Theory Manuals, Hibbitt, Karlsson & Sorensen Inc., 2001. [18] J. Lemaitre, Continuum Damage Mechanics Theory and Applications, Spinger-Verlag, 1987. [19] S. Chandrakanth, P.C. Pandey, Eng. Fract. Mech. 50 (1995) 457. [20] V. Tvergaard, On localization in ductile materials containing spherical voids, Int. J. Fract. 18 (1982) 237. [21] J. Llorca, A. Martin, J. Ruiz, M. Elices, Metall. Transact. A, Phys. Metall. Mater. Sci. 24 (1993) 1575.
[22] D.B. Zahl, R.M. Mc Meeking, Acta Metall. 39 (1991) 1117. [23] G. Bao, R.M. Meeking, J.W. Hutchinson, Acta Metall. Mater. 39 (1991) 1871. [24] W. Wu, K. Sadeghipour, K. Boberick, G. Baran, Mater. Sci. Eng. A332 (2002) 362. [25] T.S. Chow, J. Polym. Sci. 16 (1978) 959. [26] S.M. Tang, P. Cheang, M.S. Abu Bakar, K.A. Khor, K. Liao, Int. J. Fatigue 26 (2004) 49. [27] L. Farrissey, S. Schmauder, M. Dong, E. Soppa, M.H. Poech, P. McHugh, Comput. Mater. Sci. 15 (1999) 1. [28] G. Levita, A. Marchetti, A. Lazzeri, Polym. Compos. 10 (1989) 39. [29] J.F. Fan, C.P. Tsui, C.Y. Tang, C.L. Chow, Biomaterials, 2004 (in print).
Modeling of the mechanical behavior of HA/PEEK biocomposite under quasi-static tensile load J.P. Fan, C.P. Tsui∗ , C.Y. Tang Department of Industrial and Systems Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong, PR China Received 9 January 2004; received in revised form 29 April 2004
Abstract A cylindrical three-phase unit-cell model was established for predicting mechanical properties of hydroxyapatite (HA)-reinforced polyetheretherketone (PEEK) biocomposite. The model consists of an elastic–brittle HA spherical particle, an elasto-plastic matrix, and a very thin interphase region between the particle and the matrix. Perfect bonding was initially assumed among three phases before occurrence of particle–matrix debonding. For simulation of the particle–matrix debonding process of HA/PEEK, a damage evolution equation was incorporated into the material properties of the interphase region. A ductile damage evolution in the matrix was also taken into account, and was treated as a criterion for the failure of the biocomposite. A user subroutine for the damage-coupled material properties for the matrix and the interphase region was built and incorporated into a finite element code named ABAQUS. The influence of interfacial adhesion on the non-linear stress–strain relation of HA/PEEK has also been simulated by varying the tensile strength of the interphase region. The strength of the interphase region in the composite could be estimated by a comparison between the predicted and experimental results. By using the unit-cell modeling technique incorporated with the stress-based debonding and the composite failure criterion and selection of appropriate interphase strength, the numerical simulation was found to be close to the experimental results from available literatures. In vivo performance of HA/PEEK as implant materials was also discussed. © 2004 Elsevier B.V. All rights reserved. Keywords: Unit-cell model; HA/PEEK biocomposite; Finite element method (FEM); Interface debonding; Matrix damage; Particle volume fraction (PVF)
1. Introduction Metallic materials, such as stainless steel and cobalt– chromium alloys, have been successfully used as bone implants in orthopedic areas for many years. Recently, the newly developed polymers have become feasible for replacing these conventional materials due to their advantageous properties and high performance [1–4]. Among these polymers, polyetheretherketone (PEEK) is a good choice because of its excellent mechanical properties. Owing to its favorable structural and load-bearing functions, PEEK has gained much attention in the area of aerospace and marine industries [5]. In addition, it possesses exceptionally good chemical behavior, fatigue resistance, high-temperature durability, and excellent wear properties. From the processing viewpoint, PEEK can be easily formed using conventional plas-
∗
Corresponding author. Tel.: +852 2766 6635; fax: +852 2362 5267. E-mail address: [email protected] (C.P. Tsui).
0921-5093/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.msea.2004.04.078
tic processing methods, such as compression and injection molding, facilitating a rapid and economical production of implant devices [6,7]. Additional benefits of PEEK including its capacity to be repeatedly sterilized and shaped readily by machining and heat contouring render its suitability to become a biomaterial [8]. However, stiffness of PEEK cannot satisfy the requirement as a load-bearing implant material for replacing human bone. If the stiffness of the implant material such as the metallic one is too high, the newly generated living tissue will experience a stress-shielding problem [9] so that fracture healing can hardly function normally. On the other hand, if the stiffness of the implant material like some polymers is too low, the undesirable callus will form at the fracture site [10]. Hence, there is a need to develop an implant material that matches the elastic modulus of bone in order to minimize the stress-shielding effect as well as formation of the callus during the process of fracture healing. Therefore, incorporation of a second-phase hard material into the polymer matrix becomes inevitable for enhancing mechanical properties of the resulting polymer compos-
342
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
ites. Since the bone tissue of mammals has been found to contain 58% of calcium hydroxyapatite (HA), great efforts have been made to develop phosphate ceramics as a potential implant material [11]. Artificially synthesized HA has proven to be a successful substitute for the natural one. This material is not only biocompatible and non-toxic, but also exhibits unique osteoinductive properties [12]. HA/PEEK has been conceived as a biomaterial for bone replacement on the basis of producing appropriate mechanical compatibility as well as necessary biocompatibility. Besides considering whether the material stiffness can satisfy the requirement for load-bearing applications, the ultimate tensile strength of HA/PEEK must be adequate to prevent the material from catastrophic fracture. The implant made of biocomposite should undergo a sufficient geometrical deformation to show its ductile behavior. Some experimental observations from other investigators [12] have revealed that ductility of the biocomposites is greatly reduced by increasing HA particles. In the meantime, the ultimate tensile strength was found to decrease due to the presence of a poor interfacial adhesive bond. It is so because the increase in particle volume fraction (PVF) results in an increase in the area of debonded interface, which in turn causes a weakening of the ultimate tensile strength. Therefore, for HA/PEEK and all other potential composites utilizing HA particles, prediction of their mechanical properties using certain technique shows to be crucial for designing materials with certain desirable properties. Elastic properties and failure processes of HA particles-reinforced biocomposites have been investigated [1,2]. However, due to constraint of the linear elastic analysis employed in those investigations, their predictions for ductile failure processes may not attain a high accuracy. Hence, the elasto-plastic large deformation analysis is required for predicting the material properties of biocomposites. Recently, experiments have been performed by Bakar and his co-workers [3,4,13,14] to investigate the mechanical properties of HA/PEEK biocomposite for possible usage as bone replacement. Micro-structural analysis was carried out in their investigations to correlate the structure–property relationship of the composite. The dependency of the tensile properties, such as elastic modulus, strength, strain to fracture, and micro-hardness, on PVF of the HA particles were also presented. However, the failure mechanism of this biocomposite could not be completely revealed in their investigations. Thus, theoretical and numerical techniques are required for prediction of material properties. In this paper, a three-phase unit-cell model incorporating a large deformation and elasto-plastic constitutive equation has been built for predicting the mechanical properties of HA particles-reinforced biocomposite using finite element method (FEM). Interfacial debonding around the interphase region was considered in the formulation of the model by using the vanishing finite element technique [15]. Damage evolution in the matrix was also taken into account for serving as a composite failure criterion. The failure mode
transformation would be revealed from the damage distribution. The model was also applied to predict the non-linear stress–strain relation of the HA/PEEK biocomposites under the influence of different strength values of the interphase region. The optimum interphase strength would be estimated by a comparison between the predicted and experimental results. The model with the selected interphase strength, the debonding criterion, and the composite failure criterion would be further applied to predict the tensile strength and strain to fracture of HA/PEEK.
2. Elasto-plastic constitutive equations The HA particles in this study, which are brittle and with high modulus, could be accurately characterized by a pure elastic constitutive law. The polymer matrix was assumed to undergo elasto-plastic deformation with isotropic hardening, while the interphase region was assumed to have the same material properties as the matrix except for the damage evolution equation and tensile strength. Moreover, elastic strains of the HA particles were assumed to be very small with respect to unity, whereas their rotations could be arbitrarily large. The behavior of the HA particles can then be reasonably described by the following hypoelastic law e Jij = Cijkl Dkl
(1)
where σ J is the Jaumann derivative of the Cauchy stress tensor, D denotes the strain rate tensor. Ce expresses the fourth-order elastic constants tensor, which is represented by the Lame’s constants λ and µ in an isotropic material as e Cijkl = λδij δkl + µ(δik δjl + δil δjk )
(2)
where δ is the Kronecker delta. The elastic strains were assumed to be small in the matrix and the interphase region. Thus, the constitutive law used for the matrix and the interphase region is similar to Eq. (1), but the elastic strain rate is replaced by subtraction of plastic strain rate from the total strain rate, namely, p
e Jij = Cijkl (Dkl − Dkl )
(3)
where Dp is the plastic strain rate tensor. Assuming that the flow rule is derived from a plastic potential, it is easily found that the plastic strain rate Dp is coaxial with the deviator σ of the Cauchy stress tensor. P Dkl =
3p˙ σ 2σeq kl
(4)
where p˙ is the equivalent plastic strain rate and σeq is the von Mises equivalent stress. They can be expressed as p p p p p p p˙ = 2/3[(˙ε11 − ε˙ 22 )2 + (˙ε22 − ε˙ 33 )2 + (˙ε33 − ε˙ 11 )2 p2
p2
p2
+ 6(˙ε12 + ε˙ 23 + ε˙ 31 )]1/2
(5)
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
and
σeq
343
Substituting Eq. (10) into Eq. (9) gives [(σ11 − σ22 )2 + (σ22 − σ33 )2 + (σ33 2 + τ 2 + τ 2 )]1/2 − σ11 )2 + 6(τ12 23 31 = √ 2
eq
(6)
This constitutive law for the matrix and the interphase region can be represented by an one-dimensional relation between ε˙ p and σeq . In this study, the power hardening law is used as n ε σ = σ0 (7) for ε > ε0 ε0 where n denotes the hardening index, σ0 is the initial yield stress and ε0 is the corresponding strain such that ε0 = σ0 /E.
3. General description of dynamic explicit formulation Finite element programs based on an explicit dynamic formulation have been proven to be a very powerful tool in industrial application because of its numerical efficiency and other advantages, such as low memory requirement, no iteration convergent problem, easy treatment of contact condition and availability of the vanishing element treatment, as compared with implicit methods. The general form of the principle of virtual work is described as follows [16] σij δεij d + ρu¨ i δui d − ρfi δui d ρgi δui dA = 0 (8) − Af
where u and u¨ represent the displacement and accelerate vectors, r denotes mass density, f and g are the body and surface force vectors. δu and δε denote the virtual displacement and strain, respectively. Ω represents the volume of the body under study, and Af denotes the surface which is subject to an external applied force. The left-hand side of the above equation involves the terms for internal work, inertia work, work done by the body force, and work exerted by the surface force. Meshing the body with finite elements and introducing material behavior model, element shape function, and dynamics of rigid body, a finite element equation of motion is obtained and expressed in the following form at the time step n Mij u¨ nj = Rni − Fin
Mij un+1 = eq Fin j
(11)
where eq M and eq F denote the equivalent mass and equivalent applied force matrix, respectively. If eq M is lumped as a diagonal matrix, matrix inversion is not needed and a solution can be obtained directly by solving a linear equation. The lumping scheme is computationally economical because the matrix inversion consumes large computing time. Often in dynamic analysis, use of the lumping mass achieves more accurate results than that of consistent mass [17]. The central difference method has a selective convergence according to the magnitude of t, and the accuracy and convergence are linearly proportional to the square of t [17]. After the time increment of every element, the magnitude of global time increment can be determined by using the following equation: tn+1 = α t,
where
t = min{ t1 , t2 , . . . , tN } (12)
where N is the total element number and ti is the time increment of the ith element. The safety constant, a, is often selected to be less than 0.9. The critical time increment is determined such that tc = Ls /C, where Ls is the characteristic length, which is the given element area divided by the largest edge. Moreover, √the propagation speed C is determined such that C = Et /ρ, where Et is the tangent modulus and r represents the material density.
4. Ductile damage model for the matrix Under excessive external applied load, micro-voids or micro-cracks may appear in the polymer matrix. These micro-defects, which lead to material degradation, are non-reversible processes. In the present study, a damage variable, D, was used to characterize the degradation of the material’s load-bearing capability. Fig. 1 shows the damaged body of a representative volume element (RVE). A0 is the nominal intersection area and Aeff is the effective sectional
(9)
where M is the mass matrix, Rn is the overall external applied force vector and F n is the nodal force vector. In order to obtain the solution at the time step n + 1 when all of the variables at the time step n are known, the central difference method for the time discretization of the acceleration is introduced: u¨ n =
u˙ n+1/2 − u˙ n−1/2 un+1 − 2un + un−1 = t 2 t 2
(10)
Fig. 1. Representative volume element for a damaged material.
344
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
area without the presence of micro-defects and their mutual interaction. Thus, a damage variable, D, is expressed as D=1−
Aeff A0
(13)
According to this definition, the damage variable is between 0 and 1. Generally, D = 0 represents the case for an undamaged material and D = 1 denotes the case for material failure, but for most instances, the critical damage value cannot exceed a certain critical value Dcr = 1 [18]. In this study, a potential function ϕ for a ductile damage process in the polymer matrix is taken as follows [19] S −Y 2 p˙ ϕ= (14) α/n 2 S D p2/n
Fig. 2. Morphology of flame spherical HA particles.
where α and S are material-dependent constants. The damage ˙ can be obtained as evolution D ∂ϕ Y p˙ ˙ = D = (15) α/n ∂(−Y) S D p2/n The strain energy density release rate is in the form as Y=
2 R σeq v
2E(1 − D)2
(16)
where the triaxiality function Rv is described by Rv = 2 2 3 (1 + ν) + 3(1 − 2ν)(σm /σeq ) . Thus, the damage evolution ˙ D may be calculated for any loading process defined by the histories of the accumulated plastic strain rate p, ˙ the von Mises equivalent stress σeq and the hydrostatic stress σm . The damage-coupled material behaviors of the elasto-plastic polymer matrix were incorporated into ABAQUS/Explicit finite element program through a user-defined material subroutine: VUMAT sub-module.
5. Unit-cell modeling 5.1. Formation of a unit-cell model Bakar et al. [4] showed in their thermal spraying experiment that the HA particle is spherical in shape with a smooth surface as shown in Fig. 2. Although these spherical particles could hardly be evenly dispersed in a polymer composite, the micro-structure of the composite was assumed to be periodical in their investigation so that it could be represented by repeated three-dimensional hexagonal cells. As the use of three-dimensional hexagonal cells involves a larger amount of calculating time as compared with that of two-dimensional cells, a two-dimensional axisymmetric cylindrical cell [20], which has been extensively used to analyze the mechanical behaviors of particulate-reinforced composites [21,22], was used in the present study. The simplification process of the unit-cell model is illustrated as shown in Fig. 3. In spite of the gap between the exact cell and its approximated one, the simulation using the cylindri-
Fig. 3. Simplification process for the unit-cell model: (a) hexagonal array; (b) hexagonal cell model; and (c) cylindrical cell model.
cal cell only results in a relatively low error with a limited influence on the overall response of the composite [23]. Each representative cell contains a hard particle situated at its centre, which is surrounded by the polymer matrix cylinder. Due to the assumed axisymmetry for the packing of the spherical particles, only a half of the representative cell is analyzed, as shown in Fig. 4. The half height of the cylindrical cell is set equal to its radius. The volume fraction, Vp , of the particle is related to its radius. In this analysis, the constituent material properties, the particle shape, and the interphase thickness were kept fixed while only PVF was varied up to 40%. The material properties of each constituent phase are summarized in Table 1 [1,4].
Fig. 4. Geometry of the axisymmetric cylindrical cell model.
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350 Table 1 Material properties of HA and PEEK Materials
E(Gpa)
ν
ρ (kg/m3 )
PEEK HA
3.2 85
0.42 0.3
1291 3160
A third-phase material, the interphase region with a thickness of 1% of the cell radius, was introduced between the particle and the matrix [24]. Perfect bonding was initially assumed between the particle and the interphase region, and between the interphase region and the polymer matrix before occurrence of particle–matrix debonding. When damage and failure took place within the interphase region, interfacial debonding would occur and result in a degradation of the strength and stiffness of the composite. The cell model was then converted into a finite element cell model constructed by axisymmetric four-node linear elements. The finite element mesh of the cell model containing 10 vol.% particles is shown in Fig. 5, which clearly illustrates that the model is composed of the particle, the interphase region and the polymer matrix. 5.2. Boundary conditions The influence of adjacent material deformation was taken into account by using boundary conditions. By symmetry, there exist following displacement boundary conditions: uz = 0
on
z=0
(17)
ur = 0
on
r=0
(18)
where uz and ur represent the axial and radial displacement, respectively. To ensure compatibility among all periodic representative cells, the line at r = R0 of the cylindrical cell is required to remain straight and parallel to their initial state after deformation. Thus, the normal force acting on the edge at r = R0 must be set to zero such that σr dA = 0 on r = R0 (19) A1
where A1 is the hoop area of the cylinder, and all shearing stresses on all boundaries were set to zero.
Fig. 5. Finite element cell model with 10% particle volume fraction.
345
In order to satisfy this requirement, an arbitrary point P, which does not belong to any part of the model and is allowed to freely move in any directions, was created and related to all nodes on the line r = R0 . Thus, displacement of all element nodes on the line r = R0 in the r-direction are the same and equal to that of the point P. (ur )r=R0 = (ur )P
(20)
With this method, no external force is applied onto this edge. Hence, the boundary conditions set in Eq. (19) is satisfied. The velocity boundary condition is forced on the upper edge at z = H0 . In order to minimize the inertia effect, the loading velocity, v, was taken as 1 m/s.The overall uniaxial true stress, sth , and logarithm strain response, elog , of the composite are σth =
1 Tz (1 + ε) and A
εlog = ln (1 + ε)
(21)
where ε denotes the nominal strain in the axial direction and Tz is the traction applied in the same direction. A stands for current area of the cell normal to the axial direction During small-range deformation, no interfacial debonding and matrix plastic deformation take place. The overall deformation remains within the linear elastic range. Thus, the Young’s modulus of the composite, Ec , can be estimated by Ec =
(σz )z=H0 (εz )ave
(22)
where (σz )z=H0 = Tz /A is the nominal stress, and (εz )ave = (uz )z=H0 /H0 is the nominal strain. In small deformation, the nominal stress and strain are more or less the same as that of the true ones. 5.3. Particle–matrix debonding criterion It was revealed, from the tensile fracture surface of HA/PEEK samples [4], that extensive interfacial debonding between the HA particle and the polymer had occurred (Fig. 6). Thus, the composite was supposed to undergo interfacial debonding during deformation at certain critical stress. For the sake of simulating the particle–matrix debonding process, a very thin interphase region with a thickness of 1% of the cell radius was incorporated between the particle and the matrix. The interphase was initially assumed to have the same material properties as the polymer matrix except for the damage evolution and tensile strength. Perfect bonding was initially assumed between the particle and interphase and between the interphase and the polymer matrix before occurrence of the interfacial debonding. When radial strain in the interphase region exceeds a critical value, damage would be developed. The damage evolution, Di , in the interphase region is defined as Di =
σn σc
(23)
346
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
Fig. 6. Tensile fracture surface of HA/PEEK biocomposite.
where σ n and σ c are radial and critical stresses, respectively. The radial stress is defined as 1 1 σn = (σr + σz ) + (σr − σz ) cos(2ϕ) + τrz sin(2ϕ) (24) 2 2 where angle ϕ is shown in Fig. 4. During the damage process of the interphase region, the stiffness will degrade. When the radial stress σ c reaches its critical value, the failure was assumed to occur in the interphase region. When all of the material points in a certain element fail, all stress components in this element would be set to be zero. As this element looses its load-bearing capability, it would be removed from the mesh through the vanishing finite element technique [15] for simulation of the interfacial debonding process. All of these procedures have been accomplished in ABAQUS/Explicit analysis through a user-defined material subroutine: VUMAT sub-module. 6. Results and discussions 6.1. Young’s modulus In order to verify the proposed numerical model, Chow’s analytical solution [25] for the elastic modulus of a composite system with ellipsoidal particles based on Eshelby approach is used for comparison. For spherical particle, the Young’s modulus is reduced to (Kp /Km − 1) 2(Gp /Gm − 1) Ec = Em 1 + Vp + Vp 3B 3A (25) where Vp is PVF. K and G are bulk and shear moduli. The subscripts p and m represent the particle and the matrix, respectively. The material parameters are given by Kp A=1+ − 1 (1 − Vp )α, Km Gp B =1+ − 1 (1 − Vp )β (26) Gm
Fig. 7. Young’s modulus of HA/PEEK vs. volume percentage of the HA particles.
in which 1 1 + νm α= 3 1 − νm
and
β=
2 15
4 − 5νm 1 − νm
(27)
where νm is Poisson’s ratio of the matrix. The dependency of the predicted Young’s modulus of HA/PEEK on the HA content is plotted in Fig. 7. The predicted Young’s modulus increases with an increase in HA content. This means that the elastic modulus of the composite has been enhanced by the addition of HA particles. The results generated from both the numerical prediction and the analytical Chow’s formula are close to the experimental data [13,14]. In addition to the Chow’s formula, our proposed model could provide a simulation of stress–strain distribution and damage process in the composite. It can be clearly observed that HA/PEEK with 40% PVF of HA has Young’s modulus 3.5 times greater than the unfilled PEEK. Young’s modulus of cortical bone tissue is about 7–30 GPa [4]. Therefore, the stiffness of HA/PEEK with 40% PVF is comparable to that of the cortical bone. 6.2. Stress–strain relation The interaction between the HA particle and the PEEK matrix is not mechanical interlocking only. Experimental observation of the fracture surface of HA/PEEK [26] showed that most HA particles have fibrous PEEK residue adhering to the surface of the particle. Hence, there are some forms of interfacial adhesion between the particle and the matrix. In order to study the influence of the interfacial adhesion on the mechanical behaviors of the composite like HA/PEEK, tensile strength of the interphase region was varied between 40% and 80% of that of the polymer matrix. These composites were supposed to undergo interfacial debonding as usual. In addition, two extreme cases were also taken into account in this study for comparison. One of them is that the interphase region was assumed to have the same tensile strength as that of the matrix and to have no interfacial
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
347
Fig. 8. Stress–strain relation of HA/PEEK with 10% PVF and different interphase strength.
Fig. 9. Tensile strength of HA/PEEK vs. volume percentage of the HA particles.
debonding during deformation. The other extreme case is that the composite was assumed to have no interfacial bonding or strength so that all filled particles have been fully debonded before loading and thus contain only voids equal to the volume fraction of the fully debonded particles. Fig. 8 shows the effects of different interphase strength on the predicted stress–strain relation of HA/PEEK with 10% PVF. Within elastic deformation range, there is no distinguished effect of different interphase strengths on the stiffness of the composites. During plastic deformation, the stress of the composite with stronger interphase strength becomes significantly higher than that with weaker interphase strength. Hence, a proper choice of coupling agent is necessary for enhancing the interfacial strength and adhesion as well as the stiffness and strength of the composite. In addition, the non-linear stress–strain relation of the composite with the interphase region equal to 40% of the matrix strength is in good agreement with the experimental data [4]. Thus, the interphase strength was set to be about 40% of the matrix strength in this study.
shown in Fig. 9 because of the interfacial debonding. The numerically predicted results are in good agreement with the experimental data [4] and also close to the results of our previous model with an interphase having a different damage evolution equation [29]. During the interfacial debonding process, particles usually do not fully debond from the polymer matrix. Thus, there is a certain stress transfer from the matrix to the particle; the particle can still help to carry part of overall loading. As most of the loading becomes principally carried by the remaining matrix, the tensile strength of the composite with partially debonded particles is reduced. This explains why the numerically predicted results are much higher than that of the empirical model, which assumed no interfacial adhesion.
6.3. The ultimate tensile strength In most cases, the ultimate tensile strength could be enhanced by addition of particles, especially in metallic-matrix composite [27]. Unfortunately, this advantage fails to happen in biocomposite due to the poor interfacial bond. A simple empirical model that assumed no interfacial adhesion taking place between the matrix and the particle is expressed as [28] 2/3
σc = σm (1 − 1.21Vp )
(28)
where sc and sm are the tensile strength of the composite and the matrix, respectively. The predicted ultimate tensile strength decreases with increasing PVF and falls below that of the unfilled PEEK as
6.4. The strain to fracture In addition to consideration of the interfacial debonding in this study, the matrix damage was also taken into account. The critical damage value Dcr was taken to be 0.6. This value was obtained by fitting the experimental results. When the damage in any elements of the matrix reaches this critical value, catastrophic failure was assumed to happen in the composite. The logarithm strain determined using Eq. (21) is defined as the strain to fracture. Fig. 10 shows that there is an inverse relation between the strain to fracture and PVF. The strain to fracture was found to drop drastically before 10% PVF and then goes down gradually. Thus, the addition of the HA particles results in reduction of the ductility of the composite. Moreover, brittle failure dominates after 20% PVF. The predicted results are close to the experimental data [4,14]. The deformed finite element meshes at the moment of composite failure are shown in Fig. 11. Among the composites under study, only the composite with 10% PVF experiences a large deformation, while the composites with a
348
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
6.5. In vivo performance of HA/PEEK as implant materials
Fig. 10. Strain to fracture of HA/PEEK vs. volume percentage of the HA particles.
higher PVF experiences brittle fracture. The damage fields for 10% and 40% PVF are shown in Fig. 12. The damage path, along which the failure takes place, in the composite with 40% PVF is much shorter than that with 10% PVF. This explains why the composite with 40% PVF has undergone a brittle fracture.
Particulate-reinforced bicomposites like HA/PEEK, HAPEX, and HA/PLLA have been extensively developed for orthopedic applications, particularly as load-bearing implant materials [13]. PEEK is one of the high-performance polymer matrices suitable for biomedical applications because of its good chemical and fatigue resistance, high-temperature durability, good wear properties, and ability to be repeatedly sterilized before implantation [4]. An in vivo investigation has been investigated by Abu Bakar et al. [13] to study the biological response of the HA/PEEK implant materials inserted into pigs up to 24 months. Effects of varying porosity, pore size distribution, and volume content of HA particles in the implant materials on the biological response were also studied in their investigations. During early stages of implantation, fibroblast cells that promote vascularization were found. At longer implantation periods, osteroid and osteocytes were formed within lamellar bone to enable osteoblastic activites. Their in vivo investigations suggested that HA/PEEK had favorable bioactivity and biocompatibility suitable for load-bearing implant applications with long-term performance. Thus, our proposed model can be applied by material designers to design an implant
Fig. 11. Interfacial debonding at the moment of composite failure: (a) 10% PVF; (b) 20% PVF; (c) 30% PVF; and (d) 40% PVF.
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
349
Fig. 12. Damage distribution in the polymer matrix: (a) 10% PVF; and (b) 40% PVF.
material with desired mechanical properties and enough load-bearing capability.
7. Conclusions High-performance composites, such as particulatereinforced biocomposite like HA/PEEK, have got great development for load-bearing orthopedic applications. Prediction of their mechanical properties becomes important for material design and application. The cylindrical unit-cell model has been successfully applied to predict the mechanical properties of HA/PEEK biocomposites, including elastic modulus, ultimate tensile strength, stress–strain relation, and strain to fracture. The initiation of the particle–matrix debonding process has been taken into account on the basis of the stress-based criterion, which was only applied to the interphase region of the unit-cell model. Matrix damage evolution was also considered in this model to serve as a criterion for composite failure. A user-defined VUMAT material subroutine describing the damage-coupled elasto-plastic constitutive behavior of the intephase region and the matrix has been successfully established and incorporated into the ABAQUS finite element code for prediction of the mechanical properties. The mechanical behaviors of the particulate biocomposites have been found to be much influenced by the interfacial adhesion, which was successfully characterized by using different strengths of the interphase region in this study. Thus, the strength of the interphase region in the composite with 10% PVF could be estimated by comparison between the predicted and experimental results. The cell model incorporated with the proposed debonding, the failure criterion, and the selected interphase strength is found to be successful in predicting not only the trend of the mechanical properties of the HA/PEEK composites in relation with the HA content but also the non-linear constitutive relation and the ductile-to-brittle failure transition of HA/PEEK. Generally, the Young’s modulus increases with the increase in the HA content while the tensile strength and the strain
to fracture show a decrease with the increasing particle concentration. The predicted results are close to the experimental data from available literatures. This implies that the use of the stress-based debonding criterion as well as the composite failure criterion has been justified. Thus, the cylindrical unit-cell model is well suitable for prediction of the non-linear constitutive relation of the HA/PEEK composites. With this model, the interfacial debonding process and the damage path in the composite could be visualized for the material designer.
Acknowledgements The authors would like to thank the Research Grant Council of Hong Kong for funding support of this project (PolyU5176/00E) and Prof. Khor for providing us the SEM pictures of the HA/PEEK biocomposite.
References [1] F.J. Guild, W. Bonfield, Biomaterials 14 (1993) 985. [2] F.J. Guild, W. Bonfield, J. Mater. Sci.: Mater. Med. 9 (1998) 497. [3] M.S. Abu Bakar, P. Cheang, K.A. Khor, J. Mater. Process. Technol. 89–90 (1999) 462. [4] M.S. Abu Bakar, P. Cheang, K.A. Khor, Mater. Sci. Eng. A345 (2003) 55. [5] F.N. Cogswell, D.C. Leach, Plast. Rubber Compos. Process. Appl. 18 (1992) 249. [6] A. Lustiger, F.S. Uralil, G.M. Newaz, Polym. Compos. 11 (1990) 65. [7] J.R. Sarasua, P.M. Remiro, J. Mater. Sci. 30 (1995) 3501. [8] A.J. Barton, R.D. Sagers, W.G. Pitt, J. Biomed. Mater. Res. 30 (1996) 403. [9] J.V. Sloten, L. Labey, R.V. Audekercke, G.V. Perre, Biomaterials 19 (1998) 1455. [10] K.A. Jockisch, S.A. Brown, T.W. Bauer, K. Merritt, J. Biomed. Mater. Res. 26 (1992) 133. [11] G. Hienderauer, T. McGee, R. Kudej, Am. Ceram. Soc. Bull. 70 (1991) 1010. [12] G. Zambonin, M. Grano, Biomaterials 16 (1995) 397. [13] M.S. Abu Bakar, M.H.W. Cheng, S.M. Tang, S.C. Yu, K. Liao, C.T. Tan, K.A. Khor, P. Cheang, Biomaterials 24 (2003) 2245.
350
J.P. Fan et al. / Materials Science and Engineering A 382 (2004) 341–350
[14] M.S. Abu Bakar, P. Cheang, K.A. Khor, Compos. Sci. Tech. 63 (2003) 421. [15] G.C. Li, X.W. Ling, H. Shen, Int. J. Plast. 16 (2000) 39. [16] K.L. Bathe, Finite Element Procedure in Engineering Analysis, Prentice-Hall, New Jersey, 1982. [17] ABAQUS/Explicit: User’s Examples and Theory Manuals, Hibbitt, Karlsson & Sorensen Inc., 2001. [18] J. Lemaitre, Continuum Damage Mechanics Theory and Applications, Spinger-Verlag, 1987. [19] S. Chandrakanth, P.C. Pandey, Eng. Fract. Mech. 50 (1995) 457. [20] V. Tvergaard, On localization in ductile materials containing spherical voids, Int. J. Fract. 18 (1982) 237. [21] J. Llorca, A. Martin, J. Ruiz, M. Elices, Metall. Transact. A, Phys. Metall. Mater. Sci. 24 (1993) 1575.
[22] D.B. Zahl, R.M. Mc Meeking, Acta Metall. 39 (1991) 1117. [23] G. Bao, R.M. Meeking, J.W. Hutchinson, Acta Metall. Mater. 39 (1991) 1871. [24] W. Wu, K. Sadeghipour, K. Boberick, G. Baran, Mater. Sci. Eng. A332 (2002) 362. [25] T.S. Chow, J. Polym. Sci. 16 (1978) 959. [26] S.M. Tang, P. Cheang, M.S. Abu Bakar, K.A. Khor, K. Liao, Int. J. Fatigue 26 (2004) 49. [27] L. Farrissey, S. Schmauder, M. Dong, E. Soppa, M.H. Poech, P. McHugh, Comput. Mater. Sci. 15 (1999) 1. [28] G. Levita, A. Marchetti, A. Lazzeri, Polym. Compos. 10 (1989) 39. [29] J.F. Fan, C.P. Tsui, C.Y. Tang, C.L. Chow, Biomaterials, 2004 (in print).