- Email: [email protected]
plastic mechanical properties of coatings via combined nanoindentation and non-linear finite element modeling
Progress in Organic Coatings 47 (2003) 312–323
Multi-parameter models of the viscoelastic/plastic mechanical properties of coatings via combined nanoindentation and non-linear finite element modeling Timothy C. Ovaert a,∗,1 , Byung Ro Kim b , Jianjun Wang c a
Department of Aerospace and Mechanical Engineering, University of Notre Dame, 365 Fitzpatrick Hall, Notre Dame, IN 45556 5637, USA b L.G. Chemical Ltd., Daejeon, South Korea c ABB Robotics, Windsor, CT, USA
Abstract Coatings deteriorate from a variety of failure mechanisms. To improve coating durability and/or enable structure-performance correlations, it is necessary to develop more advanced methods of mechanical characterization. For complex multifunctional coatings, multi-parametric constitutive models that simultaneously account for elastic, viscoelastic, and plastic mechanical properties should be used, especially when mechanical properties in the (macro-scale) bulk state differ from properties that occur as a result of thin-film application, post-treatment processes, or aging effects. The nanoindentation creep experiment combined with non-linear finite element modeling of nanoindentation is an effective tool for characterizing the properties of such coatings. Three- and four-parameter viscoelastic/plastic finite element models, implemented using the ABAQUSTM commercial finite element software, have been developed to simulate the isotropic indentation response of coatings. Unified constitutive models where both plastic and viscoelastic deformation are considered simultaneously have not been published previously within the indentation modeling literature. The parameters are determined by an optimization program that automatically matches the load vs. indentation deformation plot from the nanoindentation experiment, with the load vs. indentation deformation plot obtained by the finite element simulation. The computed parameters become a unique “thumbprint” for a particular coating. These parameters may then be used as input data for more complex simulations, for example, capable of computing stress and strain fields, strain energy dissipation, residual stress, and residual strain during particulate scratching; or various other forms of mechanical loading. © 2003 Elsevier B.V. All rights reserved. Keywords: Nanoindentation; Multi-parameter models; Finite element modeling
1. Introduction Coatings are used extensively in a wide variety of system and component designs. It has become increasingly important for both manufacturers and end users to obtain more detailed information on coating mechanical properties, which influences their resistance to stress, deformation, fracture, wear, etc. While the uniaxial tensile test is a commonly-used method for measuring mechanical properties of fully-cured materials in bulk form, it is difficult to utilize this method for thin coatings since accurate intrinsic properties (such as Young’s modulus) are often difficult to obtain, due to test method anomalies. Completely cured coatings are often se∗ Corresponding author. Tel.: +1-574-631-9371. E-mail address: [email protected] (T.C. Ovaert). 1 Department of Aerospace and Mechanical Engineering, University of Notre Dame, Notre Dame, IN 46556, USA.
0300-9440/$ – see front matter © 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0300-9440(03)00145-0
lected based on their resistance to scratching, substrate adhesion capabilities, and resistance to deformation, all of which may be assessed to some degree by mechanical characterization. In addition, related properties may also be determined during the curing phase, where coatings transform from the liquid to the solid phase. For these reasons, the nanoindentation technique, which was originally developed to measure hardness, is widely used to characterize the mechanical properties of many coating materials. Nanoindentation is conceptually similar to most other indentation hardness measurements, however, since an extremely low load is applied to the coating (requiring specialized load/displacement transducers and vibration and acoustic isolation), the properties may be measured without influence from the substrate material, even at sub-micron coating thicknesses. A typical nanoindentation load vs. time curve, applied to a sample material, along with its corresponding load vs. indentation depth (the maximum indentation depth at the center of the
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
Fig. 1. Typical load vs. time (a) and load vs. indentation depth (b) curves during nanoindentation of inelastic material, showing unloading slope, S; maximum indentation depth at maximum load, hmax ; final indentation depth at zero load, hf ; and extrapolated indentation depth [2], hc .
indenter) curve may be seen in Fig. 1. It is evident that the material “creeps” at maximum load, indicative of an inelastic response while being held at the maximum load (Fig. 1(b)). The determination of mechanical properties via the nanoindentation test began with Doerner and Nix’s work in the 1980s [1]. Using a pyramid-shaped (Berkovich) diamond tip, they assumed that the unloading curve was controlled by a purely elastic process and that the contact area remained constant during the initial unloading phase. When these assumptions are satisfied, the contact area may be computed from the indentation depth, and Young’s modulus of the sample material, E, may be calculated using the following equations: √ π S Er = (1) √ 2 A0 1 − ν02 1 − ν2 1 = + Er E E0
(2)
A0 = 24.5h2
(3)
where S is the slope of the unloading curve (as shown in Fig. 1(b)), A0 the contact area, h indentation depth, Er the reduced elastic modulus, E and ν are the Young’s modulus and Poisson’s ratio of the sample material, respectively, and
313
E0 and ν0 the similar parameters for the indenter (usually diamond). Oliver and Pharr [2], however, realized that many materials did not follow the linear curve during unloading, as proven by experiments with different materials. They claimed that the slope should be measured continuously during loading and unloading to achieve more precise properties since the load vs. indentation depth curve was non-linear even in the initial stages. Unlike previous work, they proposed an extrapolated indentation depth, hc (Fig. 1(b)), which lies between hf and hmax to calculate the contact area. Compared to Doerner and Nix, their results produced a better estimate of Young’s modulus, which meant that their depth, hc , was more appropriate for contact area estimation. Though Young’s modulus is an important mechanical property, other properties such as the yield stress and hardening index are also necessary to fully describe the behavior of elasto-plastic materials. Cheng and Cheng [3] performed dimensional analyses with elasto-plastic solids exhibiting work hardening. Using scaling relationship and FE models, empirical formulas from the loading and unloading curves were derived for a conical indenter. Based on these empirical relations, they determined the effective yield stress, hardening index, and Young’s modulus. The methods that utilize the unloading curve to determine Young’s modulus are powerful and simple, but they are questionable for some coating materials because the unloading curve is influenced by viscoelastic effects (time-dependent relaxation). Even when an accurate contact area is measured experimentally, the unloading curve itself may not be purely elastic and the slope (S) varies with the unloading rate. To overcome this problem, Loubet et al. [4] developed a depth-sensing indentation (DSI) technique with the idea that a sinusoidal load applied to the sample yields a sinusoidal displacement for a viscoelastic material. As a result, storage and loss moduli may be estimated from the amplitude ratio and phase angle between the load and displacement, according to the following relation: √ π K Er = (4) √ 2 A0 where Er is the (viscoelastic) storage modulus and K the contact stiffness. Er will be equal to E for purely elastic materials, and lower when viscoelastic behavior is present. This method is applicable for pure viscoelastic materials, however, the issue of how to obtain an accurate contact area, A0 , still remains. As many researchers have focused on viscoelastic simulations, plastic deformation has often been overlooked in the modeling of coating materials. However, it is common for permanent deformation to remain after indentation; thus plasticity should be accounted for in describing the behavior of inelastic materials. Cheng et al. [5] studied the spherical indentation problem of a viscoelastic body, which was modeled as a standard solid model composed of two springs and a dashpot. In a subsequent paper [6], they added plastic behavior to one of
314
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
the springs on the same model. Young’s modulus was obtained from the unloading curve, yield strength from Tabor’s approximation [7] and Johnson’s cavity model [8], and viscoelastic properties were determined from the method by the previous paper [5]. Though all types of deformation were included in the calculation, they did not consider any coupling effect. That is, viscoelastic properties were obtained under the assumption of no plastic deformation, and yield strength was obtained by only considering one (elastic) spring element. Unified constitutive models where both plastic and viscoelastic deformations are considered simultaneously have not been published previously within the indentation modeling literature. There are, however, some related resources in the field of study in rheology [9–11]. Bardenhagen et al. [11] proposed a modified standard solid model consisted of a Maxwell model and an elasto-plastic spring. They used the material constants within the constitutive equations not only for pure mechanical property description but also for an explanation of the effects of temperature, humidity, etc. Their results are quite useful in describing material behavior and are readily implemented in finite element algorithms. Recently, Kim [12] and Wang [13] followed the method of Bardenhagen et al. [11] and proposed five-parameter [12], four-parameter [13], and three-parameter [13] models to describe viscoelastic/plastic mechanical behavior. The main objective of this paper is to present the results of their three- and four-parameter finite element models. The constants in the models are determined by an iterative approach using a combination of nanoindentation experiments with numerical simulations. Unlike Bardenhagen et al., several constants were omitted (for temperature, humidity, and a non-linear dashpot) to decrease the number of constants. For both parametric models, Poisson’s ratio is an input variable. In addition, the four-parameter model utilizes the
Ramberg–Osgood relation for modeling elasto-plastic behavior, while the three-parameter model utilizes an input yield stress and linear hardening parameter.
2. Experimental The nanoindentation method is well established and the details are presented elsewhere [2,14]. The maximum load capability of the nanoindenter system is 10 mN with less than 1 nN load resolution. The transducer and motion control systems are housed in an enclosure with acoustic and vibration suppression systems. The indenter tip used in these studies is a 1 m radius hemispherical diamond. There are various types of indentation tests that may be performed, each having unique loading profiles (e.g., ramp, creep, relaxation, sinusoidal). In this investigation, creep tests were performed, with constant loading and unloading rates (2 N/s) and a set holding time (25, 50, 75, or 100 s). This corresponded to a maximum load of 50, 100, 150, or 200 N, as shown in Fig. 2. Creep tests were used mainly because they facilitated the corresponding FE modeling, and strain rates are too low to induce viscoelastic heating (internal friction), which may occur with high frequency sinusoidal tests. The sample used in this study is a polyester-based material, fabricated by an air knife coating process on a vinyl substrate of 25 m thickness. The polyester was then end-capped with acrylic acid to yield a 100% solid ultraviolet curable resin. The UV cure was a total of 2.0 doses (J/cm2 ) divided approximately 33/67 between a nitrogen atmosphere and an air atmosphere. Poisson’s ratio was input as 0.49, based on bulk material estimates. Fig. 3 presents the indentation depth vs. time and load vs. indentation depth curves obtained after applying the loading
Fig. 2. Load vs. time curves for polyester-based material used in this study.
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
315
Fig. 3. Indentation depth vs. time (a) and load vs. indentation depth (b) for polyester-based material.
functions from Fig. 2 on the polyester sample. These curves were reasonably consistent even though data were obtained at different locations and with different maximum loads. A reasonable amount of creep implies that deformation of the sample is dependent on time as well as load, and creep cannot be overlooked in mechanical modeling. Permanent deformation remained after loads were completely removed, verifying that coupled plastic deformation should also be considered when describing the mechanical behavior of the sample. While load is proportional to the indentation depth to the power of 1.5 for the Hertzian elastic contact, the shape of the loading curve in Fig. 3(b) is linear. This fact verifies the influence of plastic and/or time-dependent deformation on the results. The fact that the slope of the loading curve was nearly identical for every test suggests that the transition
from the elastic region to the plastic region is smooth, and that it is difficult to pinpoint a distinct yield stress. This behavior is typical of tough, rubbery polymers and other types of coatings.
3. Four-parameter finite element indentation model It can be concluded from the above nanoindentation experiments that a specific constitutive material model is required that includes both viscoelastic and plastic behavior. To accomplish that, four-parameter constitutive equation was developed and implemented in a finite element model. The model essentially consists of a linear dashpot in parallel with an elasto-plastic spring, as shown in Fig. 4. The main task in the finite element formulation is to determine the tangent
316
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
the above equation (for each stress component) after rearranging may be rewritten as d 1 2 (t) σij
= (ηv + 2ηd ) ε11 + ηv ε22 + ηv ε33 d + σ11 |t t, . . .
d 1 2 (t) σ12
(10)
d = 2ηd ε12 + σ12 |t t, . . .
(11)
As a result, the Jacobian matrix becomes Fig. 4. Schematic representation of four-parameter model.
stiffness matrix which relates stress to strain. Since the two elements in Fig. 4 (essentially a Kelvin–Voigt model) are connected in parallel, the total stress, σ, is the sum of the stresses in the dashpot, σ d , and the elasto-plastic spring, σ s , while the strains in both cases, ε, are the same: σ =σ +σ d
s
(5)
ε = εd = εs
(6)
Here, the superscript ‘d’ refers to a dashpot and ‘s’ the spring. The equations illustrate that the total tangent stiffness matrix necessary for the FE simulation may be obtained by adding the dashpot and the elasto-plastic spring stiffness matrices, which may be determined independently. d ∂σ11 d ∂σ22 ∂σ d
ηv + 2ηd ηv ηv 2 33 = d t 0 ∂σ12 0 d ∂σ23 0 d ∂σ31 For the dashpot, whose constant is η, the one-dimensional governing equation is σ d = η˙ε
(7)
Here, ε˙ is the strain rate and η the linear damping constant. This equation may be extended to a three-dimensional form with the analogy of linear elasticity. Assuming the same Poisson’s ratio in both the volumetric and deviatoric terms, the governing differential equation for the dashpot in three dimensions has the following form: σijd = ηv ε˙ kk δij + ηd ε˙ ij
(8)
where the subscript ‘v’ refers to the volumetric term and ‘d’ the deviatoric term, and δij the Kronecker delta. With the introduction of a simple central difference operator, Eq. (8) may be written as d 1 2 (t) σij
= ηv εkk + 2ηd εij − σijd |t t
(9)
d ∂σ11 2 = (ηv + 2ηd ) ∂ε11 t
(12)
d ∂σ11 2ηv = ∂ε22 t
(13)
d ∂σ12 4ηd = ∂ε12 t
(14)
The final form of the constitutive equation for the dashpot is then expressed as d ∂εkl ∂σijd = Cijkl
(15)
If 3 × 3 stress and strain tensors are changed into 6 × 1 vectors, the constitutive equation for the dashpot has the following form: ηv
ηv
0
0
ηv + 2ηd
ηv
0
0
ηv
ηv + 2ηd 0
0 2ηd
0
0 0
0
0
0 2ηd
0
0
0
0
0 ∂ε11 ∂ε22 0 0 ∂ε 33 ∂ε 0 12 ∂ε23 0 ∂ε 31 2ηd
(16)
The finite element implementation of the plastic response utilizes a rate-independent model with the von Mises yield criterion and the associated flow rule. With its simplicity, the isotropic hardening model is used since it is not necessary to consider repeated loading and the Bauschinger effect since there was no repeated loading on the test region. For isotropic materials, the total stress and strain components may be divided into their respective volumetric and deviatoric parts. Since the volumetric component has no influence on plastic deformation, this separation in the above equations has advantages in treating a plasticity problem. With an associated flow rule, the von Mises yield function Φ may be expressed as follows: Φ = 21 τij τij − 13 σ¯ 2 (ep ) = 0
(17)
where σ¯ 2 (ep ) is yield stress and ep the equivalent plastic strain. After mathematical manipulation, the total
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
317
constitutive equation in tensor notation for the elasto-plastic spring may be expressed as
1 nˆ ⊗ nˆ s s s s ∂σ = K II + 2G J − II − G s ∂ε 3 G + (1/3)σ¯ (18) In the above equation, Ks and Gs are the bulk and shear moduli of the elasto-plastic spring, I is the fourth-order unit tensor expressed as (1/2)(δik δjl + δil δjk ), the other fourth-order tensor is defined as a cross product of the two second-order unit tensors, I ⊗ I, nˆ is the unit normal vector to the yield surface, and a prime denotes the derivative with respect to the equivalent plastic strain. Since the von Mises stress and equivalent plastic strain are not only influenced by yield stress and hardening index but also affect each other, this is a non-linear equation that must be solved via an iterative procedure. Finally, the total tangential stiffness matrix is the sum of the linear dashpot and elasto-plastic spring stiffness matrices: ∂σij = Cijkl ∂εkl p
v Cijkl = Cijkl + Cijkl
(19)
and m instead of yield strength and hardening index. As a result, the number of unknown parameters required to obtain the total stiffness matrix remains at four: η, E, A, and m, plus Poisson’s ratio which is assumed to be a known quantity.
(20)
3.1. Determination of the four parameters
Therefore, the analysis yields four unknowns (parameters): one material constant for the dashpot, η, and three constants for the elasto-plastic spring (stiffness, E; yield stress, Y; and hardening index, H). Many tough, rubbery organic coatings do not exhibit a distinct yield stress, and the transition from the elastic region to the plastic region is smooth, without any noticeable inflection points in the loading curve, as shown in Fig. 5. Thus, the Ramberg–Osgood relation is substituted for the yield stress and hardening index in modeling the elasto-plastic spring:
s m σs σ ε= (21) +A E E The first term on the right-hand side is the elastic deformation term and the second term models plastic deformation. This equation introduces the non-dimensional constants A
Fig. 5. Ramberg–Osgood elasto-plastic stress–strain behavior.
Fig. 6. Axisymmetric finite element model.
The normal indentation problem was simulated by a spherical indenter pressing against the polymer sample that was assumed to be an isotropic body. Since the tip velocity is very low, it is not necessary to consider impact in the analysis. This contact problem is highly non-linear, and thus the multi-purpose finite element code ABAQUSTM was used. Fig. 6 illustrates the finite element model mesh for the axisymmetric (two-dimensional) indentation problem. The parameters for the sample were obtained at four different loads by matching the experimental test results with the FE simulation results. Here, each experimental load vs. indentation depth curve was examined separately using the four-parameter model. The polymer sample was modeled mainly with isoparametric axisymmetric linear quadrilateral elements, and triangular elements were also used to satisfy compatibility. Since the accuracy of the finite element analysis depends on the mesh size, it is necessary to refine the portion of the mesh near the expected contact zone to obtain satisfactory results for stress and displacement. This is a normal procedure for most all contact problems of this sort. Mesh refinement was accomplished by subdividing a previously-used quadrilateral element into two or four elements. Then, the compatibility between both refined and unrefined regions was satisfied with triangular elements at the boundary. Approximately 5000 nodes and 5000 elements were employed for modeling the polymer sample. A 2 × 2 integration scheme was used for the four-node elements, and three integration points were implemented for the triangular elements. Friction was not considered in the analysis since it has only a very small effect on the normal displacement during indentation. After the constitutive equation was
318
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
implemented in an ABAQUSTM user subroutine (UMAT), the material parameters for the sample were determined by an iterative process. Table 1 summarizes the parameters obtained from each of the four loading curves. Fig. 7 shows the numerical simulation results acquired with the parameters provided in Table 1 and the experimental results shown in Fig. 3. The agreement between the numerical simulations and the experiments was satisfactory, which verifies the feasibility of using the four-parameter constitutive equation to model the time-dependent viscoelastic/plastic material behavior. The results in Table 1 illustrate that the spring constant (or stiffness) decreases as the sample is indented deeper.
This is due to a surface “skin effect” that occurs during curing with UV light. Since the intensity of the light was highest at the surface and decreases exponentially beneath the surface, additional cross-linking occurred at the surface. Although isotropic material properties have been assumed thus far for the FE simulations, one or more discrete layers or a continuous gradient in stiffness down to some predetermined depth could be modeled to obtain more precise results. The increase in the dashpot constant η revealed that the viscosity (or damping) of the sample was not constant but depended on the applied load. According to rheology theory, the relationship between viscosity and pressure may be
Fig. 7. Indentation depth vs. time (a) and load vs. indentation depth (b) experimental (solid) with corresponding finite element model fits (dashed) for polyester-based material based on parameters in Table 1 and Poisson’s ratio = 0.49.
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323 Table 1 Material parameters determined from the numerical simulations with the four-parameter model Maximum load (N)
Maximum indentation depth (nm) E (GPa) η (GPa s) A m
50
100
150
200
161.2
327.6
498.1
697.7
0.455 3.103 20 3.0
0.345 3.654 28 3.0
0.290 4.688 23 3.0
0.224 6.205 8 3.0
expressed as (after Barus [16]): η = η0 exp(αP)
(22)
where η0 is the viscosity at ambient pressure, α the pressure coefficient of viscosity, and P the hydrodynamic pressure. Though this equation is normally used for fluids, it may also be applied to model the solid-state sample. To determine η0 (the material parameter), one would input several values (two minimum) of η and P (maximum normal pressure at the surface, obtained from the finite element simulation) and thereby determine α and η0 . The introduction of the non-linear dashpot may yield more precise results. The Ramberg–Osgood (power law) relationship for elasto-plastic stress–strain behavior yields two nondimensional parameters, A and m. As can be seen in Table 1, m was unchanged in the four simulations. On the other hand, A increased slightly until the 100 N maximum load and then decreased for the 150 and 200 N simulations. In a power law relation such as Ramberg–Osgood, m represents the log ‘slope’ and A the ‘intercept’ on a log–log curve. One possible interpretation for the lack of change in m is that
319
the polymer matrix in all four simulations was the same. Thus, the physics of polymer chain movement (slip) on a molecular level would also be the same. Therefore, m represents essentially a micro-structural term that accounts for plastic flow on a micro-structural level. On the other hand, A has more to do directly with the contribution of plastic strain to the total strain and is itself a (non-dimensional) strain term. The reason for the change in A is more perplexing, however, it may relate to the way in which polymer chains are “packed” or compacted together under load. At low loads, the polymer chains will maintain their separation more readily, thus localized deformation may occur without the effect of neighboring chains. At higher loads, however, the relative amount of chain compaction is greater, and the total number of polymer chains under load is greater as well. Thus, the material behaves more like a ‘compacted’ elastic solid, with a lower relative contribution of plastic strain to the total stress–strain behavior, especially during the holding and unloading sequences.
4. Three-parameter finite element indentation model The second term in the Ramberg–Osgood equation (Eq. (21)) may be written as
m Y εp = A (23) E where Y is now a classical yield stress. Fig. 8 shows a typical hardening curve (dashed line) for a tough, rubbery polymer material. In this case, E = 100, A = 1000, and m = 3. Here, the vertical axis is yield stress, Y, and the horizontal axis is the plastic strain, !p .
Fig. 8. Yield stress (arbitrary units) vs. plastic strain for Ramberg–Osgood plasticity model and its approximation.
320
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
elastic, the estimated Young’s modulus would be larger than the actual Young’s modulus. Therefore, elastic estimation of the unloading curve yields an upper bound on E. 4.2. Viscoelastic indentation Fig. 9. Schematic representation of simplified three-parameter model.
In the actual finite element implementation of the fourparameter model discussed above, the smooth continuous Ramberg–Osgood curve is replaced by smaller piecewise line segments, as shown by the solid line. The initial yield stress, Y, is assumed to be close to zero (although it does not have to be in practice). To simplify the four-parameter model further into a three-parameter version, the Ramberg–Osgood equation is approximated by a linear hardening model. Letting the hardening parameter be H, which is the slope of the first line segment in Fig. 8 (for the given values of E, η, and m, H = 40), Eq. (21) may now be written as σ σ ε= + for loading (24) E H σ ε= for unloading (25) E The four-parameter model is now simplified to a threeparameter model, as shown in Fig. 9. When indenting with this model, during the loading period, the spring E is responsible for elastic deformation and the spring H for plastic deformation. The dashpot constant, η, and the effective spring constant of E and H (equal to EH/(E + H)) may be used to determine a viscoelastic time constant, τ. During the unloading period, only spring E is effective. 4.1. Elastic indentation For some simple model configurations, analytical solutions exist. Although one cannot use these solutions exclusively for a viscoelastic/plastic model, these solutions can yield valuable information on the indentation process. For example, when the indented material is purely elastic, the analytical relation between indentation load, P, and indentation depth, h, is given by the well-known Hertzian formula for point contact: h1.5 =
9P √ 16E R
(26)
where R is the indenter radius and E the sample Young’s modulus (assuming a diamond (rigid) indenter). For time-dependent materials, however, the unloading curve is not only influenced by elastic recovery but also by viscoelastic recovery. The viscoelasticity will impede the speed of elastic recovery. As a result, the actual unloading curve of our model will have a steeper slope than that of a purely elastic one. That is, if one treats the unloading segment as purely
If one neglects any plastic deformation, then the loading and holding segments are a viscoelastic problem, which also has an analytical solution: t Pmax 9 t − τ 1 − exp − √ τ 16E R t0 during loading
h1.5 (t) =
h1.5 (t) =
Pmax 9 √ 16E R t0
t t0 − t × t0 + τ exp − − τ exp τ τ during holding
(27)
(28)
where τ = η/E . Although both equations (27) and (28) may be used to find E and η, in Eq. (27), the linear term, which is elastic deformation, will be dominant. The estimation using Eq. (27) will be further influenced by “noise” during loading. Therefore, Eq. (28) is more applicable since the holding period is more greatly dominated by viscoelastic effects. When Eqs. (27) and (28) are used to fit the finite element simulation curve, it is found that E should be treated as the effective spring constant EH/(E + H), which is the combination of the elastic spring, E, and plastic hardening spring, H, in the simplified three-parameter model. An example can be seen in Fig. 10. A sample experimental data curve is shown in yellow. Applying the elastic estimation to the loading and unloading portions, one obtains an elastic modulus E = 123 MPa for loading, and E = 807 MPa for unloading. The viscoelastic estimation using the loading period yields E = 109 MPa, the time constant τ = 8.6 s, and η = 949 MPa s; while the estimation using the holding period gives E = 88 MPa, t = 19.3 s, and h = 1701 MPa s. Taking E = EH/(E + H) = 88 MPa, η = 1701 MPa s, different E [MPa] values (92, 163, 235, 306, 378, 449, 521, 592, 664, 735, and 807 MPa corresponding to curves 1–11, respectively) yield the results shown in Fig. 10. From Fig. 10, it can seen that the viscoelastic analytical solution is a reasonably good approximation for the effective spring constant, E , and the viscosity constant, η. Different elastic modulus (E) values yield different unloading curves, with the loading and holding portions remaining nearly the same. Note that the elastic modulus E does influence the loading and especially the holding time, although the effect is small. This influence may be explained in terms of plasticity. A “similar” plastic spring will give a larger displacement than an elastic spring under the same loading conditions.
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
321
Fig. 10. Experimental and finite element simulation curves for simplified three-parameter model.
This phenomenon may be verified by the pure tension test. Plasticity tends to distribute itself amongst all of the strain components even under the simple tension or compression test, while elastic deformation will tend to proceed more in the loading direction. Therefore, under the same loading conditions, pure elastic deformation will result in a smaller displacement than an elasto-plastic deformation. In the above example, E remains constant. Smaller E (smoother unloading) will give larger H values. A larger plastic component of E means larger loading displacement, since the “elastic equivalence” (e.g., E ) of E becomes smaller. This in turn results in a larger viscoelastic time constant (η/E ). The larger time constant will act to prevent viscoelastic recovery during holding, which implies a smaller holding displacement. This is believed to be the reason why the FE simulations yield a shorter holding period than the experimental curve when using the above estimation method. For both the four- and three-parameter models, any changes to these parameters in their respective simulations will have a large impact on the model fit. For example, changing the spring constant E has a large effect on the loading and unloading curve slopes. Changing the dashpot constant η has a large impact on the holding time. Changing the plastic parameters A, m, or H has a large effect on the final value of h.
5. Future work and application to multifunctional coatings Identification of mechanical properties is an important step in obtaining structure–property relationships of coating materials. The use of nanoindentation in conjunction with advanced finite element modeling provides a powerful tool for constructing multi-parameter representations of coating
properties, which is not limited to the typical elastic property approximations that are widely used. In order to facilitate this, one of the goals of this research is to develop stand-alone finite element simulations of the nanoindentation process, that are capable of running on PC platforms without the use of commercial finite element software. This will aid in software deployment across a wider range of potential end users. Parametric constitutive equations may also be modified as in [11] to incorporate the effects of other variables, such as temperature, humidity, or other micro-structural (or multifunctional) features. In many cases, these additional parameters may be de-coupled in the FE simulations (for example, humidity may not affect E, however, it could affect η) to facilitate greater computational efficiency. In addition, many coatings are subject to external degradation mechanisms due to abrasion, erosion, UV deterioration, etc. Utilizing the constitutive equations presented here, another future goal is to develop a stand-alone PC-based finite element model of the point scratching configuration, as shown in Fig. 11. The use of an elasto-plastic constitutive equation will more accurately determine the residual stresses, residual strains, total strain energy expended during per unit scratch length, permanent deformation, etc., during a single scratch or repeated scratches, than one would obtain using a basic elastic or elastic–perfectly plastic model which is more common. This information may then be used to assess the effects of coating functional groups and/or obtain correlations with coating in-service field performance data over time, such as gloss or reflectivity changes, chemically-induced micro-structural changes, corrosion resistance, etc. Coatings that must dry or cure over time exhibit drastic changes in properties during the curing or drying phase. For example, viscosity may change by 3–4 orders of magnitude
322
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
Fig. 11. Three-dimensional finite element model of hemispherical point scratching (after Subhash and Zhang [15]).
during cure, and there is a natural “gap” that exists when trying to bridge between liquid and solid “mechanical” behavior. Indentation methods and FE simulations may be used to examine coating property changes when going from the liquid to the solid phase, facilitated by the rheological methodologies utilized here. This information is important, for example, when assessing coating property changes as a function of time when curing or drying on inclined surfaces.
6. Conclusions 1. Though the nanoindentation technique is one of the most widely used methods for determining material properties of very thin layers, the results from elastic approximations are prone to produce errors in time-dependent materials due to their viscoelastic nature. The DSI technique, which was developed to explain the time-dependent behavior of materials, is limited because plastic deformation is not considered. Therefore, it is necessary to deal with viscoelastic and plastic properties simultaneously. 2. A commercial FEM package, ABAQUSTM , was employed for contact modeling with new constitutive equations. The constitutive equations were formulated with a conventional tangent stiffness method and implemented on an user-defined subroutine. The agreement between the numerical results and the experimental results was satisfactory.
3. Since elasto-plastic mechanical parameters of polymeric materials cannot be obtained analytically nor by experiments alone, they were acquired by matching nanoindentation test results with FE simulation results. 4. Increased computer processor speed, processor quantity, and memory storage capabilities at low cost facilitates the employment of additional parametric constants for describing coating behavior. For example, the introduction of non-linear elasticity in both the spring and dashpot may yield improved results. Other parameters that account for coating structure may also be used in the analysis, facilitating the development of multifunctional coatings.
Acknowledgements The authors would like to acknowledge the generous financial support of Armstrong World Industries Inc. Special thanks go to Drs. Jeffrey Ross, Sunil Ramachandra, and Thomas Garrett. References [1] M.F. Doerner, W.D. Nix, A method for interpreting the data from depth-sensing indentation instruments, J. Mater. Res. 1 (1986) 601– 609. [2] W.C. Oliver, G.M. Pharr, An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments, J. Mater. Res. 7 (1992) 1564–1583.
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323 [3] Y.T. Cheng, C.M. Cheng, Scaling approach to conical indentation in elastic–plastic solids with work hardening, J. Appl. Phys. 84 (1998) 1284–1291. [4] J.L. Loubet, B.N. Lucas, W.C. Oliver, Some measurements of viscoelastic properties with the help of nano-indentation, in: D.T. Smith (Ed.), Proceedings of the International Workshop on Instrumented Indentation, NIST SP 896, 1995, pp. 31–34. [5] L. Cheng, X. Xia, W. Yu, L.E. Scriven, W.W. Gerberich, Flat-punch indentation of viscoelastic material, J. Polym. Sci., Part B 38 (2000) 10–22. [6] L. Cheng, L.E. Scriven, W.W. Gerberich, Viscoelastic analysis of micro- and nano-indentation, Mater. Res. Soc. Symp. Proc. 522 (1998) 193–198. [7] D. Tabor, The Hardness of Metals, Clarendon Press, Oxford, 1951. [8] K.L. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, 1985. [9] H. Ghoneim, Y. Chen, A viscoelastic–viscoplastic constitutive equation and finite element implementation, Comput. Struct. 17 (1983) 499–509.
323
[10] C. Miehe, J. Keck, Superimposed finite elastic–viscoelastic–plastoelastic stress response with damage in filled rubbery polymers. Experiments, modelling and algorithmic implementation, J. Mech. Phys. Solids 48 (2) (2000) 323–365. [11] S.G. Bardenhagen, M.G. Stout, G.T. Gray, Three-dimensional, finite deformation, viscoplastic constitutive models for polymeric materials, Mech. Mater. 25 (1997) 235–253. [12] B.R. Kim, Modeling of particle interaction with polymer multi-layered structures, Ph.D. Thesis, The Pennsylvania State University, Pennsylvania, 1999. [13] J. Wang, Material property identification of polymer thin films under the indentation test, Ph.D. Thesis, The Pennsylvania State University, Pennsylvania, 2001. [14] W.C. Oliver, C.J. McHargue, S.J. Zinkle, Thin film characterization using a mechanical properties microprobe, Thin Solid Films 153 (1987) 185–196. [15] G. Subhash, W. Zhang, Investigation of the overall friction coefficient in single-pass scratch test, Wear 252 (2002) 123–134. [16] C. Barus, Isothermals, isopiestics, and isometrics relative to viscosity, Am. J. Sci. 45 (1893) 87–96.
Multi-parameter models of the viscoelastic/plastic mechanical properties of coatings via combined nanoindentation and non-linear finite element modeling Timothy C. Ovaert a,∗,1 , Byung Ro Kim b , Jianjun Wang c a
Department of Aerospace and Mechanical Engineering, University of Notre Dame, 365 Fitzpatrick Hall, Notre Dame, IN 45556 5637, USA b L.G. Chemical Ltd., Daejeon, South Korea c ABB Robotics, Windsor, CT, USA
Abstract Coatings deteriorate from a variety of failure mechanisms. To improve coating durability and/or enable structure-performance correlations, it is necessary to develop more advanced methods of mechanical characterization. For complex multifunctional coatings, multi-parametric constitutive models that simultaneously account for elastic, viscoelastic, and plastic mechanical properties should be used, especially when mechanical properties in the (macro-scale) bulk state differ from properties that occur as a result of thin-film application, post-treatment processes, or aging effects. The nanoindentation creep experiment combined with non-linear finite element modeling of nanoindentation is an effective tool for characterizing the properties of such coatings. Three- and four-parameter viscoelastic/plastic finite element models, implemented using the ABAQUSTM commercial finite element software, have been developed to simulate the isotropic indentation response of coatings. Unified constitutive models where both plastic and viscoelastic deformation are considered simultaneously have not been published previously within the indentation modeling literature. The parameters are determined by an optimization program that automatically matches the load vs. indentation deformation plot from the nanoindentation experiment, with the load vs. indentation deformation plot obtained by the finite element simulation. The computed parameters become a unique “thumbprint” for a particular coating. These parameters may then be used as input data for more complex simulations, for example, capable of computing stress and strain fields, strain energy dissipation, residual stress, and residual strain during particulate scratching; or various other forms of mechanical loading. © 2003 Elsevier B.V. All rights reserved. Keywords: Nanoindentation; Multi-parameter models; Finite element modeling
1. Introduction Coatings are used extensively in a wide variety of system and component designs. It has become increasingly important for both manufacturers and end users to obtain more detailed information on coating mechanical properties, which influences their resistance to stress, deformation, fracture, wear, etc. While the uniaxial tensile test is a commonly-used method for measuring mechanical properties of fully-cured materials in bulk form, it is difficult to utilize this method for thin coatings since accurate intrinsic properties (such as Young’s modulus) are often difficult to obtain, due to test method anomalies. Completely cured coatings are often se∗ Corresponding author. Tel.: +1-574-631-9371. E-mail address: [email protected] (T.C. Ovaert). 1 Department of Aerospace and Mechanical Engineering, University of Notre Dame, Notre Dame, IN 46556, USA.
0300-9440/$ – see front matter © 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0300-9440(03)00145-0
lected based on their resistance to scratching, substrate adhesion capabilities, and resistance to deformation, all of which may be assessed to some degree by mechanical characterization. In addition, related properties may also be determined during the curing phase, where coatings transform from the liquid to the solid phase. For these reasons, the nanoindentation technique, which was originally developed to measure hardness, is widely used to characterize the mechanical properties of many coating materials. Nanoindentation is conceptually similar to most other indentation hardness measurements, however, since an extremely low load is applied to the coating (requiring specialized load/displacement transducers and vibration and acoustic isolation), the properties may be measured without influence from the substrate material, even at sub-micron coating thicknesses. A typical nanoindentation load vs. time curve, applied to a sample material, along with its corresponding load vs. indentation depth (the maximum indentation depth at the center of the
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
Fig. 1. Typical load vs. time (a) and load vs. indentation depth (b) curves during nanoindentation of inelastic material, showing unloading slope, S; maximum indentation depth at maximum load, hmax ; final indentation depth at zero load, hf ; and extrapolated indentation depth [2], hc .
indenter) curve may be seen in Fig. 1. It is evident that the material “creeps” at maximum load, indicative of an inelastic response while being held at the maximum load (Fig. 1(b)). The determination of mechanical properties via the nanoindentation test began with Doerner and Nix’s work in the 1980s [1]. Using a pyramid-shaped (Berkovich) diamond tip, they assumed that the unloading curve was controlled by a purely elastic process and that the contact area remained constant during the initial unloading phase. When these assumptions are satisfied, the contact area may be computed from the indentation depth, and Young’s modulus of the sample material, E, may be calculated using the following equations: √ π S Er = (1) √ 2 A0 1 − ν02 1 − ν2 1 = + Er E E0
(2)
A0 = 24.5h2
(3)
where S is the slope of the unloading curve (as shown in Fig. 1(b)), A0 the contact area, h indentation depth, Er the reduced elastic modulus, E and ν are the Young’s modulus and Poisson’s ratio of the sample material, respectively, and
313
E0 and ν0 the similar parameters for the indenter (usually diamond). Oliver and Pharr [2], however, realized that many materials did not follow the linear curve during unloading, as proven by experiments with different materials. They claimed that the slope should be measured continuously during loading and unloading to achieve more precise properties since the load vs. indentation depth curve was non-linear even in the initial stages. Unlike previous work, they proposed an extrapolated indentation depth, hc (Fig. 1(b)), which lies between hf and hmax to calculate the contact area. Compared to Doerner and Nix, their results produced a better estimate of Young’s modulus, which meant that their depth, hc , was more appropriate for contact area estimation. Though Young’s modulus is an important mechanical property, other properties such as the yield stress and hardening index are also necessary to fully describe the behavior of elasto-plastic materials. Cheng and Cheng [3] performed dimensional analyses with elasto-plastic solids exhibiting work hardening. Using scaling relationship and FE models, empirical formulas from the loading and unloading curves were derived for a conical indenter. Based on these empirical relations, they determined the effective yield stress, hardening index, and Young’s modulus. The methods that utilize the unloading curve to determine Young’s modulus are powerful and simple, but they are questionable for some coating materials because the unloading curve is influenced by viscoelastic effects (time-dependent relaxation). Even when an accurate contact area is measured experimentally, the unloading curve itself may not be purely elastic and the slope (S) varies with the unloading rate. To overcome this problem, Loubet et al. [4] developed a depth-sensing indentation (DSI) technique with the idea that a sinusoidal load applied to the sample yields a sinusoidal displacement for a viscoelastic material. As a result, storage and loss moduli may be estimated from the amplitude ratio and phase angle between the load and displacement, according to the following relation: √ π K Er = (4) √ 2 A0 where Er is the (viscoelastic) storage modulus and K the contact stiffness. Er will be equal to E for purely elastic materials, and lower when viscoelastic behavior is present. This method is applicable for pure viscoelastic materials, however, the issue of how to obtain an accurate contact area, A0 , still remains. As many researchers have focused on viscoelastic simulations, plastic deformation has often been overlooked in the modeling of coating materials. However, it is common for permanent deformation to remain after indentation; thus plasticity should be accounted for in describing the behavior of inelastic materials. Cheng et al. [5] studied the spherical indentation problem of a viscoelastic body, which was modeled as a standard solid model composed of two springs and a dashpot. In a subsequent paper [6], they added plastic behavior to one of
314
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
the springs on the same model. Young’s modulus was obtained from the unloading curve, yield strength from Tabor’s approximation [7] and Johnson’s cavity model [8], and viscoelastic properties were determined from the method by the previous paper [5]. Though all types of deformation were included in the calculation, they did not consider any coupling effect. That is, viscoelastic properties were obtained under the assumption of no plastic deformation, and yield strength was obtained by only considering one (elastic) spring element. Unified constitutive models where both plastic and viscoelastic deformations are considered simultaneously have not been published previously within the indentation modeling literature. There are, however, some related resources in the field of study in rheology [9–11]. Bardenhagen et al. [11] proposed a modified standard solid model consisted of a Maxwell model and an elasto-plastic spring. They used the material constants within the constitutive equations not only for pure mechanical property description but also for an explanation of the effects of temperature, humidity, etc. Their results are quite useful in describing material behavior and are readily implemented in finite element algorithms. Recently, Kim [12] and Wang [13] followed the method of Bardenhagen et al. [11] and proposed five-parameter [12], four-parameter [13], and three-parameter [13] models to describe viscoelastic/plastic mechanical behavior. The main objective of this paper is to present the results of their three- and four-parameter finite element models. The constants in the models are determined by an iterative approach using a combination of nanoindentation experiments with numerical simulations. Unlike Bardenhagen et al., several constants were omitted (for temperature, humidity, and a non-linear dashpot) to decrease the number of constants. For both parametric models, Poisson’s ratio is an input variable. In addition, the four-parameter model utilizes the
Ramberg–Osgood relation for modeling elasto-plastic behavior, while the three-parameter model utilizes an input yield stress and linear hardening parameter.
2. Experimental The nanoindentation method is well established and the details are presented elsewhere [2,14]. The maximum load capability of the nanoindenter system is 10 mN with less than 1 nN load resolution. The transducer and motion control systems are housed in an enclosure with acoustic and vibration suppression systems. The indenter tip used in these studies is a 1 m radius hemispherical diamond. There are various types of indentation tests that may be performed, each having unique loading profiles (e.g., ramp, creep, relaxation, sinusoidal). In this investigation, creep tests were performed, with constant loading and unloading rates (2 N/s) and a set holding time (25, 50, 75, or 100 s). This corresponded to a maximum load of 50, 100, 150, or 200 N, as shown in Fig. 2. Creep tests were used mainly because they facilitated the corresponding FE modeling, and strain rates are too low to induce viscoelastic heating (internal friction), which may occur with high frequency sinusoidal tests. The sample used in this study is a polyester-based material, fabricated by an air knife coating process on a vinyl substrate of 25 m thickness. The polyester was then end-capped with acrylic acid to yield a 100% solid ultraviolet curable resin. The UV cure was a total of 2.0 doses (J/cm2 ) divided approximately 33/67 between a nitrogen atmosphere and an air atmosphere. Poisson’s ratio was input as 0.49, based on bulk material estimates. Fig. 3 presents the indentation depth vs. time and load vs. indentation depth curves obtained after applying the loading
Fig. 2. Load vs. time curves for polyester-based material used in this study.
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
315
Fig. 3. Indentation depth vs. time (a) and load vs. indentation depth (b) for polyester-based material.
functions from Fig. 2 on the polyester sample. These curves were reasonably consistent even though data were obtained at different locations and with different maximum loads. A reasonable amount of creep implies that deformation of the sample is dependent on time as well as load, and creep cannot be overlooked in mechanical modeling. Permanent deformation remained after loads were completely removed, verifying that coupled plastic deformation should also be considered when describing the mechanical behavior of the sample. While load is proportional to the indentation depth to the power of 1.5 for the Hertzian elastic contact, the shape of the loading curve in Fig. 3(b) is linear. This fact verifies the influence of plastic and/or time-dependent deformation on the results. The fact that the slope of the loading curve was nearly identical for every test suggests that the transition
from the elastic region to the plastic region is smooth, and that it is difficult to pinpoint a distinct yield stress. This behavior is typical of tough, rubbery polymers and other types of coatings.
3. Four-parameter finite element indentation model It can be concluded from the above nanoindentation experiments that a specific constitutive material model is required that includes both viscoelastic and plastic behavior. To accomplish that, four-parameter constitutive equation was developed and implemented in a finite element model. The model essentially consists of a linear dashpot in parallel with an elasto-plastic spring, as shown in Fig. 4. The main task in the finite element formulation is to determine the tangent
316
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
the above equation (for each stress component) after rearranging may be rewritten as d 1 2 (t) σij
= (ηv + 2ηd ) ε11 + ηv ε22 + ηv ε33 d + σ11 |t t, . . .
d 1 2 (t) σ12
(10)
d = 2ηd ε12 + σ12 |t t, . . .
(11)
As a result, the Jacobian matrix becomes Fig. 4. Schematic representation of four-parameter model.
stiffness matrix which relates stress to strain. Since the two elements in Fig. 4 (essentially a Kelvin–Voigt model) are connected in parallel, the total stress, σ, is the sum of the stresses in the dashpot, σ d , and the elasto-plastic spring, σ s , while the strains in both cases, ε, are the same: σ =σ +σ d
s
(5)
ε = εd = εs
(6)
Here, the superscript ‘d’ refers to a dashpot and ‘s’ the spring. The equations illustrate that the total tangent stiffness matrix necessary for the FE simulation may be obtained by adding the dashpot and the elasto-plastic spring stiffness matrices, which may be determined independently. d ∂σ11 d ∂σ22 ∂σ d
ηv + 2ηd ηv ηv 2 33 = d t 0 ∂σ12 0 d ∂σ23 0 d ∂σ31 For the dashpot, whose constant is η, the one-dimensional governing equation is σ d = η˙ε
(7)
Here, ε˙ is the strain rate and η the linear damping constant. This equation may be extended to a three-dimensional form with the analogy of linear elasticity. Assuming the same Poisson’s ratio in both the volumetric and deviatoric terms, the governing differential equation for the dashpot in three dimensions has the following form: σijd = ηv ε˙ kk δij + ηd ε˙ ij
(8)
where the subscript ‘v’ refers to the volumetric term and ‘d’ the deviatoric term, and δij the Kronecker delta. With the introduction of a simple central difference operator, Eq. (8) may be written as d 1 2 (t) σij
= ηv εkk + 2ηd εij − σijd |t t
(9)
d ∂σ11 2 = (ηv + 2ηd ) ∂ε11 t
(12)
d ∂σ11 2ηv = ∂ε22 t
(13)
d ∂σ12 4ηd = ∂ε12 t
(14)
The final form of the constitutive equation for the dashpot is then expressed as d ∂εkl ∂σijd = Cijkl
(15)
If 3 × 3 stress and strain tensors are changed into 6 × 1 vectors, the constitutive equation for the dashpot has the following form: ηv
ηv
0
0
ηv + 2ηd
ηv
0
0
ηv
ηv + 2ηd 0
0 2ηd
0
0 0
0
0
0 2ηd
0
0
0
0
0 ∂ε11 ∂ε22 0 0 ∂ε 33 ∂ε 0 12 ∂ε23 0 ∂ε 31 2ηd
(16)
The finite element implementation of the plastic response utilizes a rate-independent model with the von Mises yield criterion and the associated flow rule. With its simplicity, the isotropic hardening model is used since it is not necessary to consider repeated loading and the Bauschinger effect since there was no repeated loading on the test region. For isotropic materials, the total stress and strain components may be divided into their respective volumetric and deviatoric parts. Since the volumetric component has no influence on plastic deformation, this separation in the above equations has advantages in treating a plasticity problem. With an associated flow rule, the von Mises yield function Φ may be expressed as follows: Φ = 21 τij τij − 13 σ¯ 2 (ep ) = 0
(17)
where σ¯ 2 (ep ) is yield stress and ep the equivalent plastic strain. After mathematical manipulation, the total
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
317
constitutive equation in tensor notation for the elasto-plastic spring may be expressed as
1 nˆ ⊗ nˆ s s s s ∂σ = K II + 2G J − II − G s ∂ε 3 G + (1/3)σ¯ (18) In the above equation, Ks and Gs are the bulk and shear moduli of the elasto-plastic spring, I is the fourth-order unit tensor expressed as (1/2)(δik δjl + δil δjk ), the other fourth-order tensor is defined as a cross product of the two second-order unit tensors, I ⊗ I, nˆ is the unit normal vector to the yield surface, and a prime denotes the derivative with respect to the equivalent plastic strain. Since the von Mises stress and equivalent plastic strain are not only influenced by yield stress and hardening index but also affect each other, this is a non-linear equation that must be solved via an iterative procedure. Finally, the total tangential stiffness matrix is the sum of the linear dashpot and elasto-plastic spring stiffness matrices: ∂σij = Cijkl ∂εkl p
v Cijkl = Cijkl + Cijkl
(19)
and m instead of yield strength and hardening index. As a result, the number of unknown parameters required to obtain the total stiffness matrix remains at four: η, E, A, and m, plus Poisson’s ratio which is assumed to be a known quantity.
(20)
3.1. Determination of the four parameters
Therefore, the analysis yields four unknowns (parameters): one material constant for the dashpot, η, and three constants for the elasto-plastic spring (stiffness, E; yield stress, Y; and hardening index, H). Many tough, rubbery organic coatings do not exhibit a distinct yield stress, and the transition from the elastic region to the plastic region is smooth, without any noticeable inflection points in the loading curve, as shown in Fig. 5. Thus, the Ramberg–Osgood relation is substituted for the yield stress and hardening index in modeling the elasto-plastic spring:
s m σs σ ε= (21) +A E E The first term on the right-hand side is the elastic deformation term and the second term models plastic deformation. This equation introduces the non-dimensional constants A
Fig. 5. Ramberg–Osgood elasto-plastic stress–strain behavior.
Fig. 6. Axisymmetric finite element model.
The normal indentation problem was simulated by a spherical indenter pressing against the polymer sample that was assumed to be an isotropic body. Since the tip velocity is very low, it is not necessary to consider impact in the analysis. This contact problem is highly non-linear, and thus the multi-purpose finite element code ABAQUSTM was used. Fig. 6 illustrates the finite element model mesh for the axisymmetric (two-dimensional) indentation problem. The parameters for the sample were obtained at four different loads by matching the experimental test results with the FE simulation results. Here, each experimental load vs. indentation depth curve was examined separately using the four-parameter model. The polymer sample was modeled mainly with isoparametric axisymmetric linear quadrilateral elements, and triangular elements were also used to satisfy compatibility. Since the accuracy of the finite element analysis depends on the mesh size, it is necessary to refine the portion of the mesh near the expected contact zone to obtain satisfactory results for stress and displacement. This is a normal procedure for most all contact problems of this sort. Mesh refinement was accomplished by subdividing a previously-used quadrilateral element into two or four elements. Then, the compatibility between both refined and unrefined regions was satisfied with triangular elements at the boundary. Approximately 5000 nodes and 5000 elements were employed for modeling the polymer sample. A 2 × 2 integration scheme was used for the four-node elements, and three integration points were implemented for the triangular elements. Friction was not considered in the analysis since it has only a very small effect on the normal displacement during indentation. After the constitutive equation was
318
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
implemented in an ABAQUSTM user subroutine (UMAT), the material parameters for the sample were determined by an iterative process. Table 1 summarizes the parameters obtained from each of the four loading curves. Fig. 7 shows the numerical simulation results acquired with the parameters provided in Table 1 and the experimental results shown in Fig. 3. The agreement between the numerical simulations and the experiments was satisfactory, which verifies the feasibility of using the four-parameter constitutive equation to model the time-dependent viscoelastic/plastic material behavior. The results in Table 1 illustrate that the spring constant (or stiffness) decreases as the sample is indented deeper.
This is due to a surface “skin effect” that occurs during curing with UV light. Since the intensity of the light was highest at the surface and decreases exponentially beneath the surface, additional cross-linking occurred at the surface. Although isotropic material properties have been assumed thus far for the FE simulations, one or more discrete layers or a continuous gradient in stiffness down to some predetermined depth could be modeled to obtain more precise results. The increase in the dashpot constant η revealed that the viscosity (or damping) of the sample was not constant but depended on the applied load. According to rheology theory, the relationship between viscosity and pressure may be
Fig. 7. Indentation depth vs. time (a) and load vs. indentation depth (b) experimental (solid) with corresponding finite element model fits (dashed) for polyester-based material based on parameters in Table 1 and Poisson’s ratio = 0.49.
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323 Table 1 Material parameters determined from the numerical simulations with the four-parameter model Maximum load (N)
Maximum indentation depth (nm) E (GPa) η (GPa s) A m
50
100
150
200
161.2
327.6
498.1
697.7
0.455 3.103 20 3.0
0.345 3.654 28 3.0
0.290 4.688 23 3.0
0.224 6.205 8 3.0
expressed as (after Barus [16]): η = η0 exp(αP)
(22)
where η0 is the viscosity at ambient pressure, α the pressure coefficient of viscosity, and P the hydrodynamic pressure. Though this equation is normally used for fluids, it may also be applied to model the solid-state sample. To determine η0 (the material parameter), one would input several values (two minimum) of η and P (maximum normal pressure at the surface, obtained from the finite element simulation) and thereby determine α and η0 . The introduction of the non-linear dashpot may yield more precise results. The Ramberg–Osgood (power law) relationship for elasto-plastic stress–strain behavior yields two nondimensional parameters, A and m. As can be seen in Table 1, m was unchanged in the four simulations. On the other hand, A increased slightly until the 100 N maximum load and then decreased for the 150 and 200 N simulations. In a power law relation such as Ramberg–Osgood, m represents the log ‘slope’ and A the ‘intercept’ on a log–log curve. One possible interpretation for the lack of change in m is that
319
the polymer matrix in all four simulations was the same. Thus, the physics of polymer chain movement (slip) on a molecular level would also be the same. Therefore, m represents essentially a micro-structural term that accounts for plastic flow on a micro-structural level. On the other hand, A has more to do directly with the contribution of plastic strain to the total strain and is itself a (non-dimensional) strain term. The reason for the change in A is more perplexing, however, it may relate to the way in which polymer chains are “packed” or compacted together under load. At low loads, the polymer chains will maintain their separation more readily, thus localized deformation may occur without the effect of neighboring chains. At higher loads, however, the relative amount of chain compaction is greater, and the total number of polymer chains under load is greater as well. Thus, the material behaves more like a ‘compacted’ elastic solid, with a lower relative contribution of plastic strain to the total stress–strain behavior, especially during the holding and unloading sequences.
4. Three-parameter finite element indentation model The second term in the Ramberg–Osgood equation (Eq. (21)) may be written as
m Y εp = A (23) E where Y is now a classical yield stress. Fig. 8 shows a typical hardening curve (dashed line) for a tough, rubbery polymer material. In this case, E = 100, A = 1000, and m = 3. Here, the vertical axis is yield stress, Y, and the horizontal axis is the plastic strain, !p .
Fig. 8. Yield stress (arbitrary units) vs. plastic strain for Ramberg–Osgood plasticity model and its approximation.
320
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
elastic, the estimated Young’s modulus would be larger than the actual Young’s modulus. Therefore, elastic estimation of the unloading curve yields an upper bound on E. 4.2. Viscoelastic indentation Fig. 9. Schematic representation of simplified three-parameter model.
In the actual finite element implementation of the fourparameter model discussed above, the smooth continuous Ramberg–Osgood curve is replaced by smaller piecewise line segments, as shown by the solid line. The initial yield stress, Y, is assumed to be close to zero (although it does not have to be in practice). To simplify the four-parameter model further into a three-parameter version, the Ramberg–Osgood equation is approximated by a linear hardening model. Letting the hardening parameter be H, which is the slope of the first line segment in Fig. 8 (for the given values of E, η, and m, H = 40), Eq. (21) may now be written as σ σ ε= + for loading (24) E H σ ε= for unloading (25) E The four-parameter model is now simplified to a threeparameter model, as shown in Fig. 9. When indenting with this model, during the loading period, the spring E is responsible for elastic deformation and the spring H for plastic deformation. The dashpot constant, η, and the effective spring constant of E and H (equal to EH/(E + H)) may be used to determine a viscoelastic time constant, τ. During the unloading period, only spring E is effective. 4.1. Elastic indentation For some simple model configurations, analytical solutions exist. Although one cannot use these solutions exclusively for a viscoelastic/plastic model, these solutions can yield valuable information on the indentation process. For example, when the indented material is purely elastic, the analytical relation between indentation load, P, and indentation depth, h, is given by the well-known Hertzian formula for point contact: h1.5 =
9P √ 16E R
(26)
where R is the indenter radius and E the sample Young’s modulus (assuming a diamond (rigid) indenter). For time-dependent materials, however, the unloading curve is not only influenced by elastic recovery but also by viscoelastic recovery. The viscoelasticity will impede the speed of elastic recovery. As a result, the actual unloading curve of our model will have a steeper slope than that of a purely elastic one. That is, if one treats the unloading segment as purely
If one neglects any plastic deformation, then the loading and holding segments are a viscoelastic problem, which also has an analytical solution: t Pmax 9 t − τ 1 − exp − √ τ 16E R t0 during loading
h1.5 (t) =
h1.5 (t) =
Pmax 9 √ 16E R t0
t t0 − t × t0 + τ exp − − τ exp τ τ during holding
(27)
(28)
where τ = η/E . Although both equations (27) and (28) may be used to find E and η, in Eq. (27), the linear term, which is elastic deformation, will be dominant. The estimation using Eq. (27) will be further influenced by “noise” during loading. Therefore, Eq. (28) is more applicable since the holding period is more greatly dominated by viscoelastic effects. When Eqs. (27) and (28) are used to fit the finite element simulation curve, it is found that E should be treated as the effective spring constant EH/(E + H), which is the combination of the elastic spring, E, and plastic hardening spring, H, in the simplified three-parameter model. An example can be seen in Fig. 10. A sample experimental data curve is shown in yellow. Applying the elastic estimation to the loading and unloading portions, one obtains an elastic modulus E = 123 MPa for loading, and E = 807 MPa for unloading. The viscoelastic estimation using the loading period yields E = 109 MPa, the time constant τ = 8.6 s, and η = 949 MPa s; while the estimation using the holding period gives E = 88 MPa, t = 19.3 s, and h = 1701 MPa s. Taking E = EH/(E + H) = 88 MPa, η = 1701 MPa s, different E [MPa] values (92, 163, 235, 306, 378, 449, 521, 592, 664, 735, and 807 MPa corresponding to curves 1–11, respectively) yield the results shown in Fig. 10. From Fig. 10, it can seen that the viscoelastic analytical solution is a reasonably good approximation for the effective spring constant, E , and the viscosity constant, η. Different elastic modulus (E) values yield different unloading curves, with the loading and holding portions remaining nearly the same. Note that the elastic modulus E does influence the loading and especially the holding time, although the effect is small. This influence may be explained in terms of plasticity. A “similar” plastic spring will give a larger displacement than an elastic spring under the same loading conditions.
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
321
Fig. 10. Experimental and finite element simulation curves for simplified three-parameter model.
This phenomenon may be verified by the pure tension test. Plasticity tends to distribute itself amongst all of the strain components even under the simple tension or compression test, while elastic deformation will tend to proceed more in the loading direction. Therefore, under the same loading conditions, pure elastic deformation will result in a smaller displacement than an elasto-plastic deformation. In the above example, E remains constant. Smaller E (smoother unloading) will give larger H values. A larger plastic component of E means larger loading displacement, since the “elastic equivalence” (e.g., E ) of E becomes smaller. This in turn results in a larger viscoelastic time constant (η/E ). The larger time constant will act to prevent viscoelastic recovery during holding, which implies a smaller holding displacement. This is believed to be the reason why the FE simulations yield a shorter holding period than the experimental curve when using the above estimation method. For both the four- and three-parameter models, any changes to these parameters in their respective simulations will have a large impact on the model fit. For example, changing the spring constant E has a large effect on the loading and unloading curve slopes. Changing the dashpot constant η has a large impact on the holding time. Changing the plastic parameters A, m, or H has a large effect on the final value of h.
5. Future work and application to multifunctional coatings Identification of mechanical properties is an important step in obtaining structure–property relationships of coating materials. The use of nanoindentation in conjunction with advanced finite element modeling provides a powerful tool for constructing multi-parameter representations of coating
properties, which is not limited to the typical elastic property approximations that are widely used. In order to facilitate this, one of the goals of this research is to develop stand-alone finite element simulations of the nanoindentation process, that are capable of running on PC platforms without the use of commercial finite element software. This will aid in software deployment across a wider range of potential end users. Parametric constitutive equations may also be modified as in [11] to incorporate the effects of other variables, such as temperature, humidity, or other micro-structural (or multifunctional) features. In many cases, these additional parameters may be de-coupled in the FE simulations (for example, humidity may not affect E, however, it could affect η) to facilitate greater computational efficiency. In addition, many coatings are subject to external degradation mechanisms due to abrasion, erosion, UV deterioration, etc. Utilizing the constitutive equations presented here, another future goal is to develop a stand-alone PC-based finite element model of the point scratching configuration, as shown in Fig. 11. The use of an elasto-plastic constitutive equation will more accurately determine the residual stresses, residual strains, total strain energy expended during per unit scratch length, permanent deformation, etc., during a single scratch or repeated scratches, than one would obtain using a basic elastic or elastic–perfectly plastic model which is more common. This information may then be used to assess the effects of coating functional groups and/or obtain correlations with coating in-service field performance data over time, such as gloss or reflectivity changes, chemically-induced micro-structural changes, corrosion resistance, etc. Coatings that must dry or cure over time exhibit drastic changes in properties during the curing or drying phase. For example, viscosity may change by 3–4 orders of magnitude
322
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323
Fig. 11. Three-dimensional finite element model of hemispherical point scratching (after Subhash and Zhang [15]).
during cure, and there is a natural “gap” that exists when trying to bridge between liquid and solid “mechanical” behavior. Indentation methods and FE simulations may be used to examine coating property changes when going from the liquid to the solid phase, facilitated by the rheological methodologies utilized here. This information is important, for example, when assessing coating property changes as a function of time when curing or drying on inclined surfaces.
6. Conclusions 1. Though the nanoindentation technique is one of the most widely used methods for determining material properties of very thin layers, the results from elastic approximations are prone to produce errors in time-dependent materials due to their viscoelastic nature. The DSI technique, which was developed to explain the time-dependent behavior of materials, is limited because plastic deformation is not considered. Therefore, it is necessary to deal with viscoelastic and plastic properties simultaneously. 2. A commercial FEM package, ABAQUSTM , was employed for contact modeling with new constitutive equations. The constitutive equations were formulated with a conventional tangent stiffness method and implemented on an user-defined subroutine. The agreement between the numerical results and the experimental results was satisfactory.
3. Since elasto-plastic mechanical parameters of polymeric materials cannot be obtained analytically nor by experiments alone, they were acquired by matching nanoindentation test results with FE simulation results. 4. Increased computer processor speed, processor quantity, and memory storage capabilities at low cost facilitates the employment of additional parametric constants for describing coating behavior. For example, the introduction of non-linear elasticity in both the spring and dashpot may yield improved results. Other parameters that account for coating structure may also be used in the analysis, facilitating the development of multifunctional coatings.
Acknowledgements The authors would like to acknowledge the generous financial support of Armstrong World Industries Inc. Special thanks go to Drs. Jeffrey Ross, Sunil Ramachandra, and Thomas Garrett. References [1] M.F. Doerner, W.D. Nix, A method for interpreting the data from depth-sensing indentation instruments, J. Mater. Res. 1 (1986) 601– 609. [2] W.C. Oliver, G.M. Pharr, An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments, J. Mater. Res. 7 (1992) 1564–1583.
T.C. Ovaert et al. / Progress in Organic Coatings 47 (2003) 312–323 [3] Y.T. Cheng, C.M. Cheng, Scaling approach to conical indentation in elastic–plastic solids with work hardening, J. Appl. Phys. 84 (1998) 1284–1291. [4] J.L. Loubet, B.N. Lucas, W.C. Oliver, Some measurements of viscoelastic properties with the help of nano-indentation, in: D.T. Smith (Ed.), Proceedings of the International Workshop on Instrumented Indentation, NIST SP 896, 1995, pp. 31–34. [5] L. Cheng, X. Xia, W. Yu, L.E. Scriven, W.W. Gerberich, Flat-punch indentation of viscoelastic material, J. Polym. Sci., Part B 38 (2000) 10–22. [6] L. Cheng, L.E. Scriven, W.W. Gerberich, Viscoelastic analysis of micro- and nano-indentation, Mater. Res. Soc. Symp. Proc. 522 (1998) 193–198. [7] D. Tabor, The Hardness of Metals, Clarendon Press, Oxford, 1951. [8] K.L. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, 1985. [9] H. Ghoneim, Y. Chen, A viscoelastic–viscoplastic constitutive equation and finite element implementation, Comput. Struct. 17 (1983) 499–509.
323
[10] C. Miehe, J. Keck, Superimposed finite elastic–viscoelastic–plastoelastic stress response with damage in filled rubbery polymers. Experiments, modelling and algorithmic implementation, J. Mech. Phys. Solids 48 (2) (2000) 323–365. [11] S.G. Bardenhagen, M.G. Stout, G.T. Gray, Three-dimensional, finite deformation, viscoplastic constitutive models for polymeric materials, Mech. Mater. 25 (1997) 235–253. [12] B.R. Kim, Modeling of particle interaction with polymer multi-layered structures, Ph.D. Thesis, The Pennsylvania State University, Pennsylvania, 1999. [13] J. Wang, Material property identification of polymer thin films under the indentation test, Ph.D. Thesis, The Pennsylvania State University, Pennsylvania, 2001. [14] W.C. Oliver, C.J. McHargue, S.J. Zinkle, Thin film characterization using a mechanical properties microprobe, Thin Solid Films 153 (1987) 185–196. [15] G. Subhash, W. Zhang, Investigation of the overall friction coefficient in single-pass scratch test, Wear 252 (2002) 123–134. [16] C. Barus, Isothermals, isopiestics, and isometrics relative to viscosity, Am. J. Sci. 45 (1893) 87–96.