# Deformation behavior of an epoxy resin subject to multiaxial loadings. part II: constitutive modeling and predictions.

1. INTRODUCTION

The results presented in the companion paper, Part I (1), indicated that the tested epoxy had a similar behavior as other thermoset polymers, such as highly nonlinear viscoelasticity and loading path dependency. These deformation characteristics have design implications. For example, because of time-dependent behavior, dimensional changes in polymeric structures occur during their service life, and this has to be considered during the design phase of such structures. Although experimental investigations provide very useful information for the design of polymeric structures, they are generally expensive and time consuming to perform, especially for long-time behavior testing.

With the advent of high-speed computation hardware and software, numerical methods are finding more and more applications in evaluating the behavior of polymeric structures. For example, finite element methods have been widely used for the stress analysis and strength evaluation of polymeric structures (2, 3, among others). However, the success of such applications is largely dependent on the accurate description of the deformation behavior of constituents involved, especially the properties of matrices since fibers can generally be considered to remain linear elastic within their application range.

Various types of constitutive formulations have been developed to describe the mechanical behavior of polymeric materials. They can generally be grouped into two categories, viz, differential and integral forms. In the category of integral formulations for viscoelastic materials, the free volume theory of Knauss and Emri (4) and the thermodynamic based model of Lou and Schapery (5) are among those well-known and often quoted. This type of integral representation Is convenient for the stress and strain analysis, and the experimental data required for the determination of material constants is not onerous. However, their applicability for complex stress states has not been fully explored (6). Differential formulation is the other type of constitutive modeling that is often used because of its simple form and easy adaptation to finite element analysis (7, 8, among others). Because of the complex features exhibited by polymers, no single constitutive model is currently available which can describe all aspects of polymer behavior, let alone different polymers (9).

The objective of this paper is to present an improved constitutive model which describes the mechanical behavior of epoxy resins. It has been shown in Part I of this paper that there are certain features of epoxy resins which are different from those of the elastic-plastic materials, e.g. time effect, hydrostatic stress dependency. Since these features are typical of polymers, an appropriate constitutive model should be capable of predicting these characteristics. Because of multiaxial nature experienced by polymeric structures in engineering applications, a multiaxial nonlinear viscoelastic constitutive model is required and this is the focus of the present investigation.

In this paper, an improved nonlinear viscoelastic constitutive model in differential form is described. The predictions of this model are then compared with the experimental results. The differential form model is not a superposition type model, in contrast to Schapery type integral models. For a comparative purpose, the predictions of an integral viscoelastic constitutive model by Lai and Bakker (10), modified from Schapery's model (5), is also included In this paper.

2. DEVELOPMENT OF NONLINEAR VISCOELASTIC CONSTITUTIVE MODEL

The differential model presented here is physically based and is rheological. It can be viewed as a combination of Kelvin (Voigt) elements in a uniaxial representation. This type of rheological model ensures that the description of the deformation process is reliable and that the conservation laws of energy are preserved. One disadvantage of this kind of constitutive relation has been the manner in which the material constants were determined. In certain cases, a trial-and-error method had to be employed. One of the examples is the differential model by Xia and Ellyin (7). On the other hand the advantage of the differential model is that it is relatively simple to insert into general-purpose finite element codes; thus facilitating its application. In this paper, the proposed model alleviates the aforementioned drawback by using a similar procedure advocated by Zhang and Moore (11) to calibrate the material constants.

2.1 General Description of the Model

In the current model, the total strain rate {[[epsilon].sub.t]} is assumed to be the sum of the elastic and creep strain rates, {[[epsilon].sub.e]} and {[[epsilon].sub.c]} respectively, i.e.,

{[[epsilon].sub.t]} = {[[epsilon].sub.e]} + {[[epsilon].sub.c]} (1)

The elastic strain rate is calculated through Hooke's law,

{[sigma]} = E[[A].sup.-1] {[[epsilon].sub.e]} (2)

where E is an elastic modulus which is assumed to be a constant and [A] is a matrix related to the value of Poisson's ratio, defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For a number of Kelvin (Voigt) elements connected in series, creep strain rate {[[epsilon].sub.c]} is the sum of the strain rate of each element {[[epsilon].sub.ci]}, i.e.

{[[epsilon].sub.c]} = [summation over (n/i=1)] {[[epsilon].sub.ci]} = [summation over (n/i=1)]([A]/[E.sub.i][[tau].sub.i]{[sigma]} - 1/[[tau].sub.i]{[[epsilon].sub.ci]}) (4)

In the above [[tau].sub.i] = [[eta].sub.i]/[E.sub.i] (i = 1, 2, ..., n) denotes the retardation time, [E.sub.i] is the spring stiffness and [[eta].sub.i] is the dashpot viscosity for the i-th Kelvin (Voigt) element, respectively. Thus, the combination of Eqs 1 and 4 represents a model of one linear spring and several nonlinear Kelvin elements connected in series in a unixial representation. It is to be noted that the retardation time [[tau].sub.i] in Eq 4 has a damped exponential character as in an exponential-type function. Its value determines a time duration after which contribution from the individual Kelvin element becomes negligible. Therefore, the number of the Kelvin elements adopted in the constitutive equation depends on the required time range. For simplicity, we introduce a time scale factor [alpha], and assume that

[[tau].sub.i] = [([alpha]).sup.i-1] [[tau].sub.1] (5)

In this way all [[tau].sub.i] are related through the scale factor [alpha]. A time span of order of n would be covered if n Kelvin elements were adopted and the value of a was taken to be 10.

The experimental data in Part I also show clearly the nonlinear viscoelastic behavior of the studied epoxy polymer. There are several ways to incorporate the non-linearity into the constitutive equation (7, 11 among others). The description of the nonlinear behavior in the current model was achieved by letting [E.sub.i] be functions of the current equivalent stress, [[sigma].sub.eq]. Furthermore, a single function form for all [E.sub.i] is assumed. i.e.

[E.sub.i] = [E.sub.1]([[sigma].sub.eq]) (6)

The values of constants (E, v, [alpha], [[tau].sub.1]) and the function form of [E.sub.1]([[sigma].sub.eq]) can be determined from the uniaxial creep curves at different stress levels following a routine procedure. This will be described in section 3.

2.2 Description of Unloading Behavior

For most polymers, as pointed out in many references (see e.g., Ward (12) and Pae and Bhateja (13)), the unloading behavior is usually significantly different from the loading response. This is also true for many other path- and history-dependent engineering materials. For example, for elastoplastic materials, unloading usually indicates an abrupt change from plastic loading to elastic unloading. In an elastoplastic constitutive model (14) three different loading cases (elastic loading, monotonic plastic loading and plastic reloading) have been distinguished in order to correctly describe the complex plastic behavior of materials. For the viscoelastic materials, it Is also imperative that the constitutive model should have the ability to distinguish the loading and unloading cases. As such, one first has to have a criterion to distinguish unloading from loading. In plasticity, the mechanism associated with unloading is well understood and the criterion is well developed. The concept of unloading in plasticit y is closely related to the yield surface (including the subsequent yield surfaces). In viscoelasticity, there is no such clear concept akin to the yield surface. However, it is possible to introduce a current loading surface, which is defined as a stress envelope on which the current stress resides. If [[sigma].sup.t.sub.ij] is the current stress point at time t, [[sigma].sup.t.sub.ij] is the stress point at which the last switch of loading/unloading takes place and d[[sigma].sup.t.sub.ij] is the stress increment at time t, then

f([[sigma].sup.t.sub.ij] - [[sigma].sup.0.sub.ij] = [R.sup.2] (7)

defines a current loading surface, and [([partial]f/[partial][[sigma].sub.ij]).sub.[[sigma].sub.ij] = [[sigma].sup.t.sub.ij]] represents the direction of the normal of the loading surface at the current stress point. A criterion thus can be introduced to define the change from loading to unloading or vice versa as follow (see Fig. 1):

* if [partial]f/[partial][[sigma].sub.ij] d[[sigma].sub.t.sub.ij] [greater than or equtl to] 0, no chang in loading or unloading;

* if [partial]f/[partial][[sigma].sub.ij] d[[sigma].sub.t.sub.ij] < 0, either change from loading to unloading or vice versa.

The functional form off could be taken to be the same form as the von Mises surface, i.e.

f = 3/2 [S.sub.ij] [S.sub.ij] = 3/2 ([S.sub.ij] - [S.sup.0.sub.ij])([S.sub.ij] - [S.sup.0.sub.ij]) = [R.sup.2]

and [R.sup.2] = 3/2 ([S.sup.t.sub.ij] - [S.sup.0.sub.ij]([S.sup.t.sub.ij] - [S.sup.0.sub.ij]). (8)

where [S.sub.ij] = [[sigma].sub.ij] - [[sigma].sub.kk]/3 is the deviatoric part of the stress tensor.

If loading starts from a virgin state, then [[sigma].sup.0.sub.ij] = 0. At the instant at which one switches from loading to unloading or vice versa, [[sigma].sup.0.sub.ij] is up-dated to the stress point at that instant.

Having defined the criterion for the unloading, we now turn to the specification of the material constants and the modulus function for the unloading case. Phenomenological observation on the loading-unloading stress-strain curves has indicated that the slope of the loading branch decreases with the increasing stress while the trend is opposite in the unloading branch. The slope is closely related to the value of [E.sub.1]([[sigma].sub.eq]). In the loading case, the [E.sub.1]([[sigma].sub.eq]) will be a decreasing function of the applied stress which is consistent with the experimental observations (see Figs. 2 and 13). If the same function of [E.sub.1]([[sigma].sub.eq]) was used for the unloading case, then a trend for the slope change would be observed which is contrary to the experimental data. As a first approximation, it will be assumed that [E.sub.1] remains the same during the entire unloading process. In addition, to satisfy the continuity with the loading modulus function, it will be assumed that the co nstant value is equal to the value of [E.sub.1] at the switch stress point, [[sigma].sup.0.sub.ij]), where unloading takes place. All other material constants remain the same for both loading and unloading regimes. It will be seen later on that by using such a simple unloading rule, the model does predict very well the various test results including unloading cases.

2.3 Different Behaviors In Tension and Compression

For the Epon 826 epoxy, it was shown in the Part I that hydrostatic pressure has no significant effect on the deformation behavior of this epoxy. However, the material does show different behavior in uniaxial tension and compression. It was also shown in Part I that the difference in tensile and compressive deformation behavior can be adequately described by a Stassi equation. For a multiaxial stress state, the corresponding expression is as follow:

3[J.sub.2] + (R - 1)[[sigma].sup.*] [I.sub.1] = R[[sigma].sup.*2] (9)

where [I.sub.1] = [[sigma].sub.1] + [[sigma].sub.2] + [[sigma].sub.3] is the first invariant of the stress tensor, and [J.sub.2] = [S.sub.ij] [S.sub.ij]/2 is the second invariant of the deviatoric stress. In the above [[sigma].sup.*] is a uniaxial tensile yield stress in the original Stassi equation (15) and R is the ratio of a uniaxial compressive yield stress to a uniaxial tensile yield stress.

If [[sigma].sup.*] in Eq 9 is set to be the equivalent stress, [[sigma].sub.eq], then one gets,

[[sigma].sub.eq] = (R - 1)[I.sub.1] + [square root of ([(R - 1).sup.2] [I.sup.2.sub.1] + 12 [RJ.sub.2])]/2R (10)

which is different from the von Mises equivalent stress, [[sigma].sub.eq] = [square root of (3[J.sub.2])]. The equivalent stress expressed in Eq 10 is used in the current model for calculating the modulus function, [E.sub.1]([[sigma].sub.eq]). Note that when R = 1 in Eq 10, [[sigma].sub.eq] then reduced to the von Mises equivalent stress.

3. DETERMINATION OF MATERIAL CONSTANTS AND MODULUS FUNCTION

The proposed constitutive model contains a total of five material constants (E, v, [alpha], [[tau].sub.1], R) and a modulus function, [E.sub.1] ([[sigma].sub.eq]). Except for R, which has to be determined by a comparison between a uniaxial tensile curve with that of a compressive one, all other constant and function can be obtained from the experimental data of the stepped creep tests, (see Fig. 3). Assuming that the stress [[sigma].sub.0] is instantaneously reached from the zero stress state, the uniaxial strain can be reformulated from Eqs 1-6 as

[epsilon] = [[sigma].sub.0]/E + [[sigma].sub.0]/[E.sub.1]([[rho].sub.0]) [summation over (n/i=1)] (1-[e.sup.-t/[[alpha].sup.i-1][[tau].sub.1]]) (11)

The elastic modulus, E, can be determined from the stress-strain response in the initial loading stage of the stepped creep tests. If a fast loading rate is adopted in the tests, then a stable linear relationship between the stress and strain is observed. The value of E then is the slope of these straight lines. The Poisson's ratio, v, is determined from the ratio of the transverse strain to the axial strain in the uniaxial creep tests. As shown in Fig. 8 of Part I, a relatively stable value can be obtained. The time scale factor, [alpha], is usually taken to be ten. As for the number of the nonlinear Kelvin elements, a value of n = 5 to 6 will suffice. Finally, the [E.sub.1]([[rho].sub.eq]) and [[tau].sub.1] are determined from the experimental creep curves at different stress levels. For each creep stress [[rho].sub.0], Eq 11 is used to fit the corresponding creep curve by a nonlinear least squares fit procedure. Thus stress dependent functions [E.sub.1]([[rho].sub.0]) and [[tau].sub.1]([[rho].sub.0]) can be obtained. The fitting results show that [[tau].sub.1] is generally weakly dependent on the stress, therefore, it could be set to be the average value of all [[tau].sub.1]([[rho].sub.0]) values.

Following the above described procedure, all the material constants and the modulus function in the constitutive model can easily be calibrated. The material constants and the modulus function for the Epon 826 epoxy polymer were thus determined to be:

E = 3400 MPa,

v = 0.42, [alpha] = 10, [[tau].sub.1] = 6.116, R = 1.15 and

[E.sub.1][sigma] = 1.055 x [10.sup.5]e - [sigma] - 22.764/18.000 MPa (12)

Note that six Kelvin elements were used in the calculations. The individual values of [E.sub.1] ([[sigma].sub.0]) from fitting the creep curves at different stress levels are shown in Fig. 2, as well as the final fitting curve for the function [E.sub.1]([sigma]).

4. MODEL PREDICTIONS

The current differential form of nonlinear viscoelastic constitutive model has been used to predict the deformation behavior of Epon 826 epoxy polymer for the tests reported in Part I. For the sake of comparison, the predictive results of an integral form of constitutive model (10) are also included. A short description of the integral model and the calibrated material constants and functions is provided herein.

The constitutive law can be expressed as:

[[epsilon].sub.ij] = 1/2 [J.sub.0][g.sub.0][s.sub.ij](t) + 1/2 [g.sub.1] [[integral].sup.t.sub.0] [DELTA]J [[phi](t) - [phi]([tau])] d([g.sub.2][s.sub.ij])/d[tau] d[tau] + 1/9 [B.sub.0][g.sub.0] [[delta].sub.ij] [[sigma].sub.kk] (t)

+ 1/9 [[delta].sub.ij] [g.sub.1] [[integra].sup.t.sub.0] [DELTA]B [[psi](t) - [psi]([tau])] d([g.sub.2][[sigma].sub.kk])/d[tau] d[tau] (13)

where [J.sub.0] and [B.sub.0] are instantaneous shear and bulk compliances; [DELTA]J and [DELTA]B are transient shear and bulk compliances, respectively; [psi] and [psi]' are reduced times, [psi] = [[integral].sup.t.sub.0] d[tau]/[a.sub.[sigma]] and [psi]' = [[integral].sup.[tau].sub.0] d[tau]/[a.sub.[sigma]]; [g.sub.0], [g.sub.1], [g.sub.2] and [a.sub.[sigma]] are functions of multiaxial stress state. According to Lai and Bakker (10), [g.sub.0] and [a.sub.[sigma]] are functions of hydrostatic stress, [[sigma].sub.kk], and [g.sub.1] and [g.sub.2] are functions of the von Mises equivalent stress, [[sigma].sub.eq] = [square root of (3/2 [s.sub.ij][s.sub.if])].

[J.sub.0], [B.sub.0], [DELTA]J and [DELTA]B are related to the instantaneous tensile compliance [D.sub.0] and the transient tensile compliance [DELTA]D by the following relations,

[J.sub.0] = 2(1 + v)[D.sub.0] [DELTA]J([psi]) = 2(1 + v)[DELTA]D([psi])

[B.sub.0] = 3(1 - 2v)[D.sub.0] [DELTA]B([psi]) = 3(1 - 2v)[DELTA]D([psi]) (14)

where v is the same Poisson's ratio as defined in the differential model. [DELTA]D is described by a simple power law as [DELTA]D([psi]) = [D.sub.1][[psi].sup.n], [D.sub.1] and n are material constants to be determined from experiments.

As pointed out by Schapery (16), all the above material constants and functions can be determined by using a set of uniaxial creep-recovery test data at different stress levels (Fig. 3). The calibrated material constants and function relations of the integral constitutive model for the Epon 826 are listed in Table 1.

Creep Teats

The predicted creep and recovery curves by the differential constitutive model is compared with the experimental data in Fig. 3. It is to be noted that in the current model a switch from loading to unloading occurs at the instant when the creep test is terminated. It can be seen that the agreement is fairly good. At the highest stress level, the predictive accuracy deteriorates as the nonlinear viscoelastic relations do not account for the damage which may have occurred at high stress levels. The prediction by the integral model are similar to those shown in Fig. 3.

Since the model parameters for the two constitutive representations were determined from the uniaxial tensile creep tests, it is imperative to explore whether these models can accurately predict the responses for other loading conditions. The prediction for a shear creep test at a stress level of 28.95 MPa is shown in Fig. 4. The predicted shear creep curve by the differential model is very close to the experimental data while that of the integral model is not as good.

Multiaxial Behavior

The predicted stress envelopes for a fixed octahedral shear strain of [[epsilon].sub.oct] = 2.5% at three different octahedral strain rates of [10.sup.-3] [s.sup.-1], [10.sup.-4] [s.sup.-1] and [10.sup.-5] [s.sup.-1] are depicted in Fig. 5. It is observed that both the differential and integral models can predict the difference in behavior in tension and compression loading of this epoxy. However, the predicted stress envelopes by the present differential model are in much better agreement with the experimental data than that of the integral model.

The predicted stress envelopes for three strain values with the same octahedral shear strain rate of [10.sup.-4] [s.sup.-1] are shown in Fig. 6. One notes differences between the predictions of the two constitutive models. The predicted stress envelopes by the differential model are closer to the experimental data than the integral model predictions. The latter overestimates the stress envelopes for all three strain levels.

Proportional Loading Path

Figure 7 compares the experimental data of combined axial and shear proportional test (path OC) at an octahedral shear strain rate of [10.sup.-5] [s.sup.-1] with the predicted stress responses by the two models. It is seen that during the early loading, the predictions of both models fit fairly well the experimental data. However, with the increase in the stress level, the prediction by the differential model are in much better agreement with the experimental data than that of the integral model for both axial and shear stresses.

Proportional Loading Path.

The predicted axial and shear stresses by the two models for the non-proportional paths OBC are compared with the experimental data as shown in Fig. 8. Again, it is seen that the differential constitutive model predicts stress responses that are very close to the experimental results, while the prediction by the integral representation deviates significantly. At the proportional loading stage (Path OB) the integral model predicts a higher shear stress response and the same trend extends to the non-proportional stage (Path BC). In contrast, the predicted axial stress response is lower for the Path BC.

The predicted stress responses by the two models for another non-proportional paths OAC are shown in Fig. 9, together with the corresponding experimental data. The predicted trends for this loading path are similar to those shown in Fig. 8, i.e., the differential model predicts the stress variation very close to the experimental data while that predicted by the integral constitutive representation deviates significantly. In the proporitonal loading case (Path OA) the integral model predicted a higher axial stress response and the same trend is also observed in the non-proportional Path AC. A reversal trend is seen with respect to the predicted shear stress for the path AC where the predicted values are lower than the experimental ones.

From the above comparison, it is clear that the predictions by the differential constitutive model are in fairly good agreement with the experimental data, both qualitatively and quantitatively. However, in the case of the integral model only a qualitative agreement is observed.

Transient Loading Paths

Uniaxial Tensile Loading With Different Strain Rates

Figure 10 depicts the predictions for uniaxial tensile stress-strain curves at three different octahedral shear strain rates of [10.sup.-3] [s.sup.-1], [10.sup.-4] [s.sup.-1] and [10.sup.-5] [s.sup.-1]. The comparison of the predictions with the experimental data shows that at low stress levels, the prediction by both the differential and the integral models is close to each other and to the experiments. However, at high stress level, the integral model predicts much higher stress values. The prediction by the differential model is fairly good, especially for the two lower strain rates.

Uniaxial Compressive Loading With Different Strain Rates

The predictions for uniaxial compressive stress-strain curves at the same three octahedral shear strain rates as those for uniaxial tensile tests, are shown in Fig. 11. The comparison of the predicted values with the experimental data indicates that at low stress levels, the prediction by both the differential and the integral models is close to each other and to the experiments. However, at high stress level, the integral model predicts much higher stress values. The prediction by the differential model is closer to the experimental data. However, both models cannot capture the yield and subsequent strain softening behavior indicated by the experiments, since there are no damage mechanisms included in both model formulations.

Pure Shear Loading With Different Strain Rates

Figure 12 shows the predictions for pure shear stress-strain curves at three octahedral shear strain rates of [10.sup.-3] [s.sup.-1], [10.sup.-4] [s.sup.-1] and [10.sup.-5] [s.sup.-1]. The differential model predicts a very good agreement with the experimental data, especially at the two lower strain rates. The predictions for the three strain rates by the integral model are all higher than those of the experiments, especially at higher stress levels.

Equibiaxial Tensile Loading With Different Strain Rates

For the equibiaxial loading case, the predicted axial and hoop stress responses for each strain rate are the same and collapse to one curve. Again the differential model predicts the experimental data better than the integral model. The integral model predicts much higher stress responses for all three strain rates ([10.sup.-3], [10.sup.-4], [10.sup.-5] [s.sup.-1]) at high stress levels.

From the above comparison, it is clear that both models can accurately predict the deformation behavior at low stress levels. With the increase in stress level, the predicted values by the two models are different from each other. In general, the integral model tends to overestimate the stress responses for the strain-controlled tests and underestimate the strain responses for the stress-controlled tests (data for the latter are not presented here). The predictions by the differential model are in very good agreement with the experimental data for most loading cases, even at very high stress levels.

Unloading Behaviors

In the following, the ability of the current differential model to predict the unloading behavior by using the newly proposed rule for the modulus [E.sub.1], is examined. For the sake of comparison, the results by using the same loading modulus function [E.sub.1]([[sigma].sub.eq]) for both loading and unloading cases are also presented for the uniaxial tensile loading/unloading test.

Uniaxial Tensile Loading-Unloading

The predicted unloading behavior for the uniaxial tensile test at an octahedral shear strain rate of [10.sup.-5] [s.sup.-1] is shown in Fig. 13. This figure reveals that the strong nonlinear and viscoelastic response of the material is well predicted over the complete range of the test program. It is to be noted that the prediction of the unloading branch deviates significantly, if the same modulus function is applied without distinguishing between loading/unloading cases, see dotted line in Fig. 13.

Proportional Loading/Unloading Along Path O-C-O

The predicted behavior for the proportional path O-C-O with an octahedral shear strain rate of [10.sup.-5] [s.sup.-1] is shown in Fig. 14. It is noted that the proposed model can predict fairly well the biaxial (shear/axial) loading/unloading behavior of the proportional path O-C-O.

Non-Proportional Loading/Unloading Path

Figure 15 shows that the current differential model prediction for both the axial and shear loading/unloading behavior following the non-proportional path O-A-C-O. Similar trend is observed for the non-proportional path O-B-C-O. It is seen that the loading/unloading response is well predicted both qualitatively and quantitatively, by the current model.

5. CONCLUSIONS

A nonlinear differential viscoelastic constitutive model has been developed. The model can distinguish between loading and unloading cases by introducing a current loading surface and a corresponding switch rule. For the loading case, a stress-dependent modulus function is defined; for the unloading case the modulus is assumed to be constant and is equal to the value of the modulus function at the loading/unloading switch instant. Another feature of the model is the inclusion of a new definition of equivalent stress based on the Stassi equation to characterize the multiaxial stress state especially the difference in behavior in tension and compression loading of the epoxy resin. Furthermore, the model is relatively simple, it contains only five material constants, and a modulus function. The material constants and the modulus function can be easily determined from the uniaxial tensile creep curves at different stress levels through a well-defined routine procedure.

The model predictions have been compared with the comprehensive experimental data of Epon 826 epoxy polymer including creep, recovery, various transient loading/unloading paths and multiaxial proportional or nonproportional paths. The model predictions are in very good agreement with the experimental results not only qualitatively but also quantitatively.

This model is a phenomenological constitutive model focused on the mulitiaxial nonlinear viscoelastic behavior of epoxy polymers. The predictive ability of the model has also been checked for another type of epoxy polymer, ColdCure epoxy (Epon 828). see (17). However, the deformation behavior of polymers is influenced by many other factors such as cross-link density, additives, cure cycles, temperature, aging, moisture, etc. The prediction capability of the model on the influence of the aforementioned factors needs to be further investigated.

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

[FIGURE 13 OMITTED]

[FIGURE 13 OMITTED]

[FIGURE 14 OMITTED]

[FIGURE 15 OMITTED]

ACKNOWLEDGMENT

The work presented here is part of a general investigation of the behavior of polymeric composites under complex loading and adverse environment. This research is supported by the TCPL/NSERC Senior Industrial Research Chair Program (F. E.) and NSERC research grants to F. E. and Z. X.

REFERENCE

(1.) Y. Hu, Z. Xia, and F. Ellyin, Polym. Eng. Sci., this issue.

(2.) Y. Chen, Z. Xia, and F. Ellyin. J. Compos. Mater., 35, 522 (2001).

(3.) S. Yi and H. H. Hilton. Trans. ASME J. Eng. Mater. Technol., 119, 266 (1997).

(4.) W. G. Knauss and I. Emri, Polym. Eng. Sci., 27, 86 (1987).

(5.) Y. C. Lou and R. A. Schapery, J. Compos. Mater., 5, 208 (1971).

(6.) Y. Hu, Z. Xia, and F. Ellyin, Polym. Polym. Compos., 8. 157 (2000).

(7.) Z. Xia and F. Ellyin, Polym. Polym Compos., 6, 75(1998).

(8.) O. C. Zienkiewicz, M. Watson, and I. P. King. Int. J. Mech. Sci., 10, 807 (1968).

(9.) R. A. Schapery, Int. J. Solids Struct., 37, 359 (2000).

(10.) J. Lai and A. Bakker, Comput. Mech. 18, 182 (1996).

(11.) C. Zhang and D. Moore, polym. Eng. Sci., 37, 414 (1997).

(12.) I. M. Ward. Mechanical Properties of Solid Polymers, Wiley, New York (1983).

(13.) K. D. Pae and S. K. Bhateja, J. Macromol. Sci., Rev. Macromol. Chem., C13, 1 (1975).

(14.) F. Ellyin, Z. Xia, and J. Wu, Comput. Struct. 56, 283 (1995).

(15.) F. Stassi-D'Alia, Meccanica, 4, 349 (1969).

(16.) R. A. Schapery, Polym. Eng. Sci., 9, 295 (1969).

(17.) Y. Hu, "Multiaxial Behavior and Viscoelastic Constitutive Modeling of Epoxy Polymers," PhD thesis, Department of Mechanical Engineering, University of Alberta (2002).

ZIHUI KIA *

* Corresponding author. E-mail: Zthul.Xia@ualberta.ca

The results presented in the companion paper, Part I (1), indicated that the tested epoxy had a similar behavior as other thermoset polymers, such as highly nonlinear viscoelasticity and loading path dependency. These deformation characteristics have design implications. For example, because of time-dependent behavior, dimensional changes in polymeric structures occur during their service life, and this has to be considered during the design phase of such structures. Although experimental investigations provide very useful information for the design of polymeric structures, they are generally expensive and time consuming to perform, especially for long-time behavior testing.

With the advent of high-speed computation hardware and software, numerical methods are finding more and more applications in evaluating the behavior of polymeric structures. For example, finite element methods have been widely used for the stress analysis and strength evaluation of polymeric structures (2, 3, among others). However, the success of such applications is largely dependent on the accurate description of the deformation behavior of constituents involved, especially the properties of matrices since fibers can generally be considered to remain linear elastic within their application range.

Various types of constitutive formulations have been developed to describe the mechanical behavior of polymeric materials. They can generally be grouped into two categories, viz, differential and integral forms. In the category of integral formulations for viscoelastic materials, the free volume theory of Knauss and Emri (4) and the thermodynamic based model of Lou and Schapery (5) are among those well-known and often quoted. This type of integral representation Is convenient for the stress and strain analysis, and the experimental data required for the determination of material constants is not onerous. However, their applicability for complex stress states has not been fully explored (6). Differential formulation is the other type of constitutive modeling that is often used because of its simple form and easy adaptation to finite element analysis (7, 8, among others). Because of the complex features exhibited by polymers, no single constitutive model is currently available which can describe all aspects of polymer behavior, let alone different polymers (9).

The objective of this paper is to present an improved constitutive model which describes the mechanical behavior of epoxy resins. It has been shown in Part I of this paper that there are certain features of epoxy resins which are different from those of the elastic-plastic materials, e.g. time effect, hydrostatic stress dependency. Since these features are typical of polymers, an appropriate constitutive model should be capable of predicting these characteristics. Because of multiaxial nature experienced by polymeric structures in engineering applications, a multiaxial nonlinear viscoelastic constitutive model is required and this is the focus of the present investigation.

In this paper, an improved nonlinear viscoelastic constitutive model in differential form is described. The predictions of this model are then compared with the experimental results. The differential form model is not a superposition type model, in contrast to Schapery type integral models. For a comparative purpose, the predictions of an integral viscoelastic constitutive model by Lai and Bakker (10), modified from Schapery's model (5), is also included In this paper.

2. DEVELOPMENT OF NONLINEAR VISCOELASTIC CONSTITUTIVE MODEL

The differential model presented here is physically based and is rheological. It can be viewed as a combination of Kelvin (Voigt) elements in a uniaxial representation. This type of rheological model ensures that the description of the deformation process is reliable and that the conservation laws of energy are preserved. One disadvantage of this kind of constitutive relation has been the manner in which the material constants were determined. In certain cases, a trial-and-error method had to be employed. One of the examples is the differential model by Xia and Ellyin (7). On the other hand the advantage of the differential model is that it is relatively simple to insert into general-purpose finite element codes; thus facilitating its application. In this paper, the proposed model alleviates the aforementioned drawback by using a similar procedure advocated by Zhang and Moore (11) to calibrate the material constants.

2.1 General Description of the Model

In the current model, the total strain rate {[[epsilon].sub.t]} is assumed to be the sum of the elastic and creep strain rates, {[[epsilon].sub.e]} and {[[epsilon].sub.c]} respectively, i.e.,

{[[epsilon].sub.t]} = {[[epsilon].sub.e]} + {[[epsilon].sub.c]} (1)

The elastic strain rate is calculated through Hooke's law,

{[sigma]} = E[[A].sup.-1] {[[epsilon].sub.e]} (2)

where E is an elastic modulus which is assumed to be a constant and [A] is a matrix related to the value of Poisson's ratio, defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For a number of Kelvin (Voigt) elements connected in series, creep strain rate {[[epsilon].sub.c]} is the sum of the strain rate of each element {[[epsilon].sub.ci]}, i.e.

{[[epsilon].sub.c]} = [summation over (n/i=1)] {[[epsilon].sub.ci]} = [summation over (n/i=1)]([A]/[E.sub.i][[tau].sub.i]{[sigma]} - 1/[[tau].sub.i]{[[epsilon].sub.ci]}) (4)

In the above [[tau].sub.i] = [[eta].sub.i]/[E.sub.i] (i = 1, 2, ..., n) denotes the retardation time, [E.sub.i] is the spring stiffness and [[eta].sub.i] is the dashpot viscosity for the i-th Kelvin (Voigt) element, respectively. Thus, the combination of Eqs 1 and 4 represents a model of one linear spring and several nonlinear Kelvin elements connected in series in a unixial representation. It is to be noted that the retardation time [[tau].sub.i] in Eq 4 has a damped exponential character as in an exponential-type function. Its value determines a time duration after which contribution from the individual Kelvin element becomes negligible. Therefore, the number of the Kelvin elements adopted in the constitutive equation depends on the required time range. For simplicity, we introduce a time scale factor [alpha], and assume that

[[tau].sub.i] = [([alpha]).sup.i-1] [[tau].sub.1] (5)

In this way all [[tau].sub.i] are related through the scale factor [alpha]. A time span of order of n would be covered if n Kelvin elements were adopted and the value of a was taken to be 10.

The experimental data in Part I also show clearly the nonlinear viscoelastic behavior of the studied epoxy polymer. There are several ways to incorporate the non-linearity into the constitutive equation (7, 11 among others). The description of the nonlinear behavior in the current model was achieved by letting [E.sub.i] be functions of the current equivalent stress, [[sigma].sub.eq]. Furthermore, a single function form for all [E.sub.i] is assumed. i.e.

[E.sub.i] = [E.sub.1]([[sigma].sub.eq]) (6)

The values of constants (E, v, [alpha], [[tau].sub.1]) and the function form of [E.sub.1]([[sigma].sub.eq]) can be determined from the uniaxial creep curves at different stress levels following a routine procedure. This will be described in section 3.

2.2 Description of Unloading Behavior

For most polymers, as pointed out in many references (see e.g., Ward (12) and Pae and Bhateja (13)), the unloading behavior is usually significantly different from the loading response. This is also true for many other path- and history-dependent engineering materials. For example, for elastoplastic materials, unloading usually indicates an abrupt change from plastic loading to elastic unloading. In an elastoplastic constitutive model (14) three different loading cases (elastic loading, monotonic plastic loading and plastic reloading) have been distinguished in order to correctly describe the complex plastic behavior of materials. For the viscoelastic materials, it Is also imperative that the constitutive model should have the ability to distinguish the loading and unloading cases. As such, one first has to have a criterion to distinguish unloading from loading. In plasticity, the mechanism associated with unloading is well understood and the criterion is well developed. The concept of unloading in plasticit y is closely related to the yield surface (including the subsequent yield surfaces). In viscoelasticity, there is no such clear concept akin to the yield surface. However, it is possible to introduce a current loading surface, which is defined as a stress envelope on which the current stress resides. If [[sigma].sup.t.sub.ij] is the current stress point at time t, [[sigma].sup.t.sub.ij] is the stress point at which the last switch of loading/unloading takes place and d[[sigma].sup.t.sub.ij] is the stress increment at time t, then

f([[sigma].sup.t.sub.ij] - [[sigma].sup.0.sub.ij] = [R.sup.2] (7)

defines a current loading surface, and [([partial]f/[partial][[sigma].sub.ij]).sub.[[sigma].sub.ij] = [[sigma].sup.t.sub.ij]] represents the direction of the normal of the loading surface at the current stress point. A criterion thus can be introduced to define the change from loading to unloading or vice versa as follow (see Fig. 1):

* if [partial]f/[partial][[sigma].sub.ij] d[[sigma].sub.t.sub.ij] [greater than or equtl to] 0, no chang in loading or unloading;

* if [partial]f/[partial][[sigma].sub.ij] d[[sigma].sub.t.sub.ij] < 0, either change from loading to unloading or vice versa.

The functional form off could be taken to be the same form as the von Mises surface, i.e.

f = 3/2 [S.sub.ij] [S.sub.ij] = 3/2 ([S.sub.ij] - [S.sup.0.sub.ij])([S.sub.ij] - [S.sup.0.sub.ij]) = [R.sup.2]

and [R.sup.2] = 3/2 ([S.sup.t.sub.ij] - [S.sup.0.sub.ij]([S.sup.t.sub.ij] - [S.sup.0.sub.ij]). (8)

where [S.sub.ij] = [[sigma].sub.ij] - [[sigma].sub.kk]/3 is the deviatoric part of the stress tensor.

If loading starts from a virgin state, then [[sigma].sup.0.sub.ij] = 0. At the instant at which one switches from loading to unloading or vice versa, [[sigma].sup.0.sub.ij] is up-dated to the stress point at that instant.

Having defined the criterion for the unloading, we now turn to the specification of the material constants and the modulus function for the unloading case. Phenomenological observation on the loading-unloading stress-strain curves has indicated that the slope of the loading branch decreases with the increasing stress while the trend is opposite in the unloading branch. The slope is closely related to the value of [E.sub.1]([[sigma].sub.eq]). In the loading case, the [E.sub.1]([[sigma].sub.eq]) will be a decreasing function of the applied stress which is consistent with the experimental observations (see Figs. 2 and 13). If the same function of [E.sub.1]([[sigma].sub.eq]) was used for the unloading case, then a trend for the slope change would be observed which is contrary to the experimental data. As a first approximation, it will be assumed that [E.sub.1] remains the same during the entire unloading process. In addition, to satisfy the continuity with the loading modulus function, it will be assumed that the co nstant value is equal to the value of [E.sub.1] at the switch stress point, [[sigma].sup.0.sub.ij]), where unloading takes place. All other material constants remain the same for both loading and unloading regimes. It will be seen later on that by using such a simple unloading rule, the model does predict very well the various test results including unloading cases.

2.3 Different Behaviors In Tension and Compression

For the Epon 826 epoxy, it was shown in the Part I that hydrostatic pressure has no significant effect on the deformation behavior of this epoxy. However, the material does show different behavior in uniaxial tension and compression. It was also shown in Part I that the difference in tensile and compressive deformation behavior can be adequately described by a Stassi equation. For a multiaxial stress state, the corresponding expression is as follow:

3[J.sub.2] + (R - 1)[[sigma].sup.*] [I.sub.1] = R[[sigma].sup.*2] (9)

where [I.sub.1] = [[sigma].sub.1] + [[sigma].sub.2] + [[sigma].sub.3] is the first invariant of the stress tensor, and [J.sub.2] = [S.sub.ij] [S.sub.ij]/2 is the second invariant of the deviatoric stress. In the above [[sigma].sup.*] is a uniaxial tensile yield stress in the original Stassi equation (15) and R is the ratio of a uniaxial compressive yield stress to a uniaxial tensile yield stress.

If [[sigma].sup.*] in Eq 9 is set to be the equivalent stress, [[sigma].sub.eq], then one gets,

[[sigma].sub.eq] = (R - 1)[I.sub.1] + [square root of ([(R - 1).sup.2] [I.sup.2.sub.1] + 12 [RJ.sub.2])]/2R (10)

which is different from the von Mises equivalent stress, [[sigma].sub.eq] = [square root of (3[J.sub.2])]. The equivalent stress expressed in Eq 10 is used in the current model for calculating the modulus function, [E.sub.1]([[sigma].sub.eq]). Note that when R = 1 in Eq 10, [[sigma].sub.eq] then reduced to the von Mises equivalent stress.

3. DETERMINATION OF MATERIAL CONSTANTS AND MODULUS FUNCTION

The proposed constitutive model contains a total of five material constants (E, v, [alpha], [[tau].sub.1], R) and a modulus function, [E.sub.1] ([[sigma].sub.eq]). Except for R, which has to be determined by a comparison between a uniaxial tensile curve with that of a compressive one, all other constant and function can be obtained from the experimental data of the stepped creep tests, (see Fig. 3). Assuming that the stress [[sigma].sub.0] is instantaneously reached from the zero stress state, the uniaxial strain can be reformulated from Eqs 1-6 as

[epsilon] = [[sigma].sub.0]/E + [[sigma].sub.0]/[E.sub.1]([[rho].sub.0]) [summation over (n/i=1)] (1-[e.sup.-t/[[alpha].sup.i-1][[tau].sub.1]]) (11)

The elastic modulus, E, can be determined from the stress-strain response in the initial loading stage of the stepped creep tests. If a fast loading rate is adopted in the tests, then a stable linear relationship between the stress and strain is observed. The value of E then is the slope of these straight lines. The Poisson's ratio, v, is determined from the ratio of the transverse strain to the axial strain in the uniaxial creep tests. As shown in Fig. 8 of Part I, a relatively stable value can be obtained. The time scale factor, [alpha], is usually taken to be ten. As for the number of the nonlinear Kelvin elements, a value of n = 5 to 6 will suffice. Finally, the [E.sub.1]([[rho].sub.eq]) and [[tau].sub.1] are determined from the experimental creep curves at different stress levels. For each creep stress [[rho].sub.0], Eq 11 is used to fit the corresponding creep curve by a nonlinear least squares fit procedure. Thus stress dependent functions [E.sub.1]([[rho].sub.0]) and [[tau].sub.1]([[rho].sub.0]) can be obtained. The fitting results show that [[tau].sub.1] is generally weakly dependent on the stress, therefore, it could be set to be the average value of all [[tau].sub.1]([[rho].sub.0]) values.

Following the above described procedure, all the material constants and the modulus function in the constitutive model can easily be calibrated. The material constants and the modulus function for the Epon 826 epoxy polymer were thus determined to be:

E = 3400 MPa,

v = 0.42, [alpha] = 10, [[tau].sub.1] = 6.116, R = 1.15 and

[E.sub.1][sigma] = 1.055 x [10.sup.5]e - [sigma] - 22.764/18.000 MPa (12)

Note that six Kelvin elements were used in the calculations. The individual values of [E.sub.1] ([[sigma].sub.0]) from fitting the creep curves at different stress levels are shown in Fig. 2, as well as the final fitting curve for the function [E.sub.1]([sigma]).

4. MODEL PREDICTIONS

The current differential form of nonlinear viscoelastic constitutive model has been used to predict the deformation behavior of Epon 826 epoxy polymer for the tests reported in Part I. For the sake of comparison, the predictive results of an integral form of constitutive model (10) are also included. A short description of the integral model and the calibrated material constants and functions is provided herein.

The constitutive law can be expressed as:

[[epsilon].sub.ij] = 1/2 [J.sub.0][g.sub.0][s.sub.ij](t) + 1/2 [g.sub.1] [[integral].sup.t.sub.0] [DELTA]J [[phi](t) - [phi]([tau])] d([g.sub.2][s.sub.ij])/d[tau] d[tau] + 1/9 [B.sub.0][g.sub.0] [[delta].sub.ij] [[sigma].sub.kk] (t)

+ 1/9 [[delta].sub.ij] [g.sub.1] [[integra].sup.t.sub.0] [DELTA]B [[psi](t) - [psi]([tau])] d([g.sub.2][[sigma].sub.kk])/d[tau] d[tau] (13)

where [J.sub.0] and [B.sub.0] are instantaneous shear and bulk compliances; [DELTA]J and [DELTA]B are transient shear and bulk compliances, respectively; [psi] and [psi]' are reduced times, [psi] = [[integral].sup.t.sub.0] d[tau]/[a.sub.[sigma]] and [psi]' = [[integral].sup.[tau].sub.0] d[tau]/[a.sub.[sigma]]; [g.sub.0], [g.sub.1], [g.sub.2] and [a.sub.[sigma]] are functions of multiaxial stress state. According to Lai and Bakker (10), [g.sub.0] and [a.sub.[sigma]] are functions of hydrostatic stress, [[sigma].sub.kk], and [g.sub.1] and [g.sub.2] are functions of the von Mises equivalent stress, [[sigma].sub.eq] = [square root of (3/2 [s.sub.ij][s.sub.if])].

[J.sub.0], [B.sub.0], [DELTA]J and [DELTA]B are related to the instantaneous tensile compliance [D.sub.0] and the transient tensile compliance [DELTA]D by the following relations,

[J.sub.0] = 2(1 + v)[D.sub.0] [DELTA]J([psi]) = 2(1 + v)[DELTA]D([psi])

[B.sub.0] = 3(1 - 2v)[D.sub.0] [DELTA]B([psi]) = 3(1 - 2v)[DELTA]D([psi]) (14)

where v is the same Poisson's ratio as defined in the differential model. [DELTA]D is described by a simple power law as [DELTA]D([psi]) = [D.sub.1][[psi].sup.n], [D.sub.1] and n are material constants to be determined from experiments.

As pointed out by Schapery (16), all the above material constants and functions can be determined by using a set of uniaxial creep-recovery test data at different stress levels (Fig. 3). The calibrated material constants and function relations of the integral constitutive model for the Epon 826 are listed in Table 1.

Creep Teats

The predicted creep and recovery curves by the differential constitutive model is compared with the experimental data in Fig. 3. It is to be noted that in the current model a switch from loading to unloading occurs at the instant when the creep test is terminated. It can be seen that the agreement is fairly good. At the highest stress level, the predictive accuracy deteriorates as the nonlinear viscoelastic relations do not account for the damage which may have occurred at high stress levels. The prediction by the integral model are similar to those shown in Fig. 3.

Since the model parameters for the two constitutive representations were determined from the uniaxial tensile creep tests, it is imperative to explore whether these models can accurately predict the responses for other loading conditions. The prediction for a shear creep test at a stress level of 28.95 MPa is shown in Fig. 4. The predicted shear creep curve by the differential model is very close to the experimental data while that of the integral model is not as good.

Multiaxial Behavior

The predicted stress envelopes for a fixed octahedral shear strain of [[epsilon].sub.oct] = 2.5% at three different octahedral strain rates of [10.sup.-3] [s.sup.-1], [10.sup.-4] [s.sup.-1] and [10.sup.-5] [s.sup.-1] are depicted in Fig. 5. It is observed that both the differential and integral models can predict the difference in behavior in tension and compression loading of this epoxy. However, the predicted stress envelopes by the present differential model are in much better agreement with the experimental data than that of the integral model.

The predicted stress envelopes for three strain values with the same octahedral shear strain rate of [10.sup.-4] [s.sup.-1] are shown in Fig. 6. One notes differences between the predictions of the two constitutive models. The predicted stress envelopes by the differential model are closer to the experimental data than the integral model predictions. The latter overestimates the stress envelopes for all three strain levels.

Proportional Loading Path

Figure 7 compares the experimental data of combined axial and shear proportional test (path OC) at an octahedral shear strain rate of [10.sup.-5] [s.sup.-1] with the predicted stress responses by the two models. It is seen that during the early loading, the predictions of both models fit fairly well the experimental data. However, with the increase in the stress level, the prediction by the differential model are in much better agreement with the experimental data than that of the integral model for both axial and shear stresses.

Proportional Loading Path.

The predicted axial and shear stresses by the two models for the non-proportional paths OBC are compared with the experimental data as shown in Fig. 8. Again, it is seen that the differential constitutive model predicts stress responses that are very close to the experimental results, while the prediction by the integral representation deviates significantly. At the proportional loading stage (Path OB) the integral model predicts a higher shear stress response and the same trend extends to the non-proportional stage (Path BC). In contrast, the predicted axial stress response is lower for the Path BC.

The predicted stress responses by the two models for another non-proportional paths OAC are shown in Fig. 9, together with the corresponding experimental data. The predicted trends for this loading path are similar to those shown in Fig. 8, i.e., the differential model predicts the stress variation very close to the experimental data while that predicted by the integral constitutive representation deviates significantly. In the proporitonal loading case (Path OA) the integral model predicted a higher axial stress response and the same trend is also observed in the non-proportional Path AC. A reversal trend is seen with respect to the predicted shear stress for the path AC where the predicted values are lower than the experimental ones.

From the above comparison, it is clear that the predictions by the differential constitutive model are in fairly good agreement with the experimental data, both qualitatively and quantitatively. However, in the case of the integral model only a qualitative agreement is observed.

Transient Loading Paths

Uniaxial Tensile Loading With Different Strain Rates

Figure 10 depicts the predictions for uniaxial tensile stress-strain curves at three different octahedral shear strain rates of [10.sup.-3] [s.sup.-1], [10.sup.-4] [s.sup.-1] and [10.sup.-5] [s.sup.-1]. The comparison of the predictions with the experimental data shows that at low stress levels, the prediction by both the differential and the integral models is close to each other and to the experiments. However, at high stress level, the integral model predicts much higher stress values. The prediction by the differential model is fairly good, especially for the two lower strain rates.

Uniaxial Compressive Loading With Different Strain Rates

The predictions for uniaxial compressive stress-strain curves at the same three octahedral shear strain rates as those for uniaxial tensile tests, are shown in Fig. 11. The comparison of the predicted values with the experimental data indicates that at low stress levels, the prediction by both the differential and the integral models is close to each other and to the experiments. However, at high stress level, the integral model predicts much higher stress values. The prediction by the differential model is closer to the experimental data. However, both models cannot capture the yield and subsequent strain softening behavior indicated by the experiments, since there are no damage mechanisms included in both model formulations.

Pure Shear Loading With Different Strain Rates

Figure 12 shows the predictions for pure shear stress-strain curves at three octahedral shear strain rates of [10.sup.-3] [s.sup.-1], [10.sup.-4] [s.sup.-1] and [10.sup.-5] [s.sup.-1]. The differential model predicts a very good agreement with the experimental data, especially at the two lower strain rates. The predictions for the three strain rates by the integral model are all higher than those of the experiments, especially at higher stress levels.

Equibiaxial Tensile Loading With Different Strain Rates

For the equibiaxial loading case, the predicted axial and hoop stress responses for each strain rate are the same and collapse to one curve. Again the differential model predicts the experimental data better than the integral model. The integral model predicts much higher stress responses for all three strain rates ([10.sup.-3], [10.sup.-4], [10.sup.-5] [s.sup.-1]) at high stress levels.

From the above comparison, it is clear that both models can accurately predict the deformation behavior at low stress levels. With the increase in stress level, the predicted values by the two models are different from each other. In general, the integral model tends to overestimate the stress responses for the strain-controlled tests and underestimate the strain responses for the stress-controlled tests (data for the latter are not presented here). The predictions by the differential model are in very good agreement with the experimental data for most loading cases, even at very high stress levels.

Unloading Behaviors

In the following, the ability of the current differential model to predict the unloading behavior by using the newly proposed rule for the modulus [E.sub.1], is examined. For the sake of comparison, the results by using the same loading modulus function [E.sub.1]([[sigma].sub.eq]) for both loading and unloading cases are also presented for the uniaxial tensile loading/unloading test.

Uniaxial Tensile Loading-Unloading

The predicted unloading behavior for the uniaxial tensile test at an octahedral shear strain rate of [10.sup.-5] [s.sup.-1] is shown in Fig. 13. This figure reveals that the strong nonlinear and viscoelastic response of the material is well predicted over the complete range of the test program. It is to be noted that the prediction of the unloading branch deviates significantly, if the same modulus function is applied without distinguishing between loading/unloading cases, see dotted line in Fig. 13.

Proportional Loading/Unloading Along Path O-C-O

The predicted behavior for the proportional path O-C-O with an octahedral shear strain rate of [10.sup.-5] [s.sup.-1] is shown in Fig. 14. It is noted that the proposed model can predict fairly well the biaxial (shear/axial) loading/unloading behavior of the proportional path O-C-O.

Non-Proportional Loading/Unloading Path

Figure 15 shows that the current differential model prediction for both the axial and shear loading/unloading behavior following the non-proportional path O-A-C-O. Similar trend is observed for the non-proportional path O-B-C-O. It is seen that the loading/unloading response is well predicted both qualitatively and quantitatively, by the current model.

5. CONCLUSIONS

A nonlinear differential viscoelastic constitutive model has been developed. The model can distinguish between loading and unloading cases by introducing a current loading surface and a corresponding switch rule. For the loading case, a stress-dependent modulus function is defined; for the unloading case the modulus is assumed to be constant and is equal to the value of the modulus function at the loading/unloading switch instant. Another feature of the model is the inclusion of a new definition of equivalent stress based on the Stassi equation to characterize the multiaxial stress state especially the difference in behavior in tension and compression loading of the epoxy resin. Furthermore, the model is relatively simple, it contains only five material constants, and a modulus function. The material constants and the modulus function can be easily determined from the uniaxial tensile creep curves at different stress levels through a well-defined routine procedure.

The model predictions have been compared with the comprehensive experimental data of Epon 826 epoxy polymer including creep, recovery, various transient loading/unloading paths and multiaxial proportional or nonproportional paths. The model predictions are in very good agreement with the experimental results not only qualitatively but also quantitatively.

This model is a phenomenological constitutive model focused on the mulitiaxial nonlinear viscoelastic behavior of epoxy polymers. The predictive ability of the model has also been checked for another type of epoxy polymer, ColdCure epoxy (Epon 828). see (17). However, the deformation behavior of polymers is influenced by many other factors such as cross-link density, additives, cure cycles, temperature, aging, moisture, etc. The prediction capability of the model on the influence of the aforementioned factors needs to be further investigated.

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

[FIGURE 13 OMITTED]

[FIGURE 13 OMITTED]

[FIGURE 14 OMITTED]

[FIGURE 15 OMITTED]

Table 1 Integral Model Constants and Functions. [D.sub.0] = 2.86 X [10.sup.-4] [(MPa).sup.-1], [D.sub.-1] = 9.819 X [10.sup.-6] [(MPa).sup.-1], n = 0.145 [g.sub.i] = [C.sub.i] + [A.sup.i] (i = 0, 1, 2) [([sigma]/[[sigma].sub.o].sup. [B.sub.i]] [MATHEMATICAL EXPRESSIONS NOT [MATHEMATICAL EXPRESSIONS NOT REPRODUCIBLE IN ASCII] REPRODUCIBLE IN ASCII] [C.sub.0] = 1.0, [C.sub.1] = 1.0, [C.sub.2] = 1.0, [C.sub.[sigma]] = 0.2097, [[sigma].sub.0] = 80.0 [g.sub.i] = [C.sub.i] + [A.sup.i] [([sigma]/[[sigma].sub.o].sup. [B.sub.i]] [MATHEMATICAL EXPRESSIONS NOT REPRODUCIBLE IN ASCII] [C.sub.0] = 1.0, [A.sub.0] = 0.3250, [C.sub.1] = 1.0, [A.sub.1] = 2.6629, [C.sub.2] = 1.0, [A.sub.2] = 3.7015, [C.sub.[sigma]] = 0.2097, [A.sub.[sigma]] = 0.7901, [[sigma].sub.0] = 80.0 [[sigma].sub.1] = 18.97 [g.sub.i] = [C.sub.i] + [A.sup.i] [([sigma]/[[sigma].sub.o].sup. [B.sub.i]] [MATHEMATICAL EXPRESSIONS NOT REPRODUCIBLE IN ASCII] [C.sub.0] = 1.0, [B.sub.0] = 4.462 [C.sub.1] = 1.0, [B.sub.1] = 2.047 [C.sub.2] = 1.0, [B.sub.2] = 2.947 [C.sub.[sigma]] = 0.2097, [B.sub.[sigma]] = 13.845 [[sigma].sub.0] = 80.0

ACKNOWLEDGMENT

The work presented here is part of a general investigation of the behavior of polymeric composites under complex loading and adverse environment. This research is supported by the TCPL/NSERC Senior Industrial Research Chair Program (F. E.) and NSERC research grants to F. E. and Z. X.

REFERENCE

(1.) Y. Hu, Z. Xia, and F. Ellyin, Polym. Eng. Sci., this issue.

(2.) Y. Chen, Z. Xia, and F. Ellyin. J. Compos. Mater., 35, 522 (2001).

(3.) S. Yi and H. H. Hilton. Trans. ASME J. Eng. Mater. Technol., 119, 266 (1997).

(4.) W. G. Knauss and I. Emri, Polym. Eng. Sci., 27, 86 (1987).

(5.) Y. C. Lou and R. A. Schapery, J. Compos. Mater., 5, 208 (1971).

(6.) Y. Hu, Z. Xia, and F. Ellyin, Polym. Polym. Compos., 8. 157 (2000).

(7.) Z. Xia and F. Ellyin, Polym. Polym Compos., 6, 75(1998).

(8.) O. C. Zienkiewicz, M. Watson, and I. P. King. Int. J. Mech. Sci., 10, 807 (1968).

(9.) R. A. Schapery, Int. J. Solids Struct., 37, 359 (2000).

(10.) J. Lai and A. Bakker, Comput. Mech. 18, 182 (1996).

(11.) C. Zhang and D. Moore, polym. Eng. Sci., 37, 414 (1997).

(12.) I. M. Ward. Mechanical Properties of Solid Polymers, Wiley, New York (1983).

(13.) K. D. Pae and S. K. Bhateja, J. Macromol. Sci., Rev. Macromol. Chem., C13, 1 (1975).

(14.) F. Ellyin, Z. Xia, and J. Wu, Comput. Struct. 56, 283 (1995).

(15.) F. Stassi-D'Alia, Meccanica, 4, 349 (1969).

(16.) R. A. Schapery, Polym. Eng. Sci., 9, 295 (1969).

(17.) Y. Hu, "Multiaxial Behavior and Viscoelastic Constitutive Modeling of Epoxy Polymers," PhD thesis, Department of Mechanical Engineering, University of Alberta (2002).

ZIHUI KIA *

* Corresponding author. E-mail: Zthul.Xia@ualberta.ca

Printer friendly Cite/link Email Feedback | |

Author: | Kia, Zihui; Hu, Yafei; Ellyin, Fernand |
---|---|

Publication: | Polymer Engineering and Science |

Date: | Mar 1, 2003 |

Words: | 5392 |

Previous Article: | Deformation behavior of an epoxy resin subject to multiaxial loadings. Part I: experimental investigations. |

Next Article: | Scratch visibility of polymers measured using optical imaging. |

Topics: |