Nonlinear Schapery viscoelastic material model for thermoplastic polymers

The nonlinear mechanical behavior of viscoelastic materials is modeled based on the Schapery integral model containing internal variables. In this context, a new approach for the strain-dependent material properties is introduced considering a one-dimensional formulation with strain-dependent nonlinear functions for an oscillatory load case. In addition to the viscoelastic storage and loss modulus, the higher order harmonic oscillations in the stress response are computed and compared to experimental data from Fourier transform rheology of Polyamide 6 (PA6). The com-parison reveals a good agreement between the predictions of the nonlinear model and the experimental data for the higher harmonic intensity I 2 = 1 .


| INTRODUCTION
The application of polymer-based material systems in the automotive and aircraft sectors has increased significantly in recent years due to their potential as low-weight materials, combining low weight with high specific stiffness and strength, their processability, and their relatively low cost. Discontinuous and continuous fiber reinforced composites, for instance, allow to construct and design structural parts that are tailored to specific industrial applications. 1 In such composites, thermosets and thermoplastics are often used as polymeric matrix. In terms of thermoplastic polymers, the solid material behavior is characterized by pronounced thermo-viscoelastic material properties. [2][3][4] Weidenmann and co-workers investigated the influence of the strain-rate on the elastic modulus for a long fiber reinforced composite. 5 Extensive work was published to experimentally characterize pure and reinforced polymers in terms of their strength, 6 damage, 7 or crack behavior. 8 Moreover, the influence of temperature, 9 crystallization, 10 and aging 11 on the material properties has been investigated. There are many approaches to describe the thermomechanical behavior of thermoplastic-based composites. Attempts have been made to predict the effective viscoelastic properties of fiber reinforced composites, 12 the effective thermoelastic properties, 13 or the damage behavior. 14 The material behavior of viscoelastic materials is in general rate-and temperature-dependent. In linear viscoelasticity, the Boltzmann-superposition that describes the stressstrain relationship by a single integral is usually applied. 15 Numerous approaches exist for modeling nonlinear viscoelastic material behavior in the solid state. According to Brinson and Brinson,15 common approaches are given by nonlinear mechanical models that introduce nonlinear spring and damper elements, or nonlinear creep power law approaches. To fit experimental data obtained by one-step uniaxial creep tests or tests with a constant strain-rate, multiple-or singleintegral approaches were introduced. Since the Boltzmannsuperposition is no longer applicable in the nonlinear regime in case of multiple step stress, varying stress, or multiaxial stress loadings, the multiple-as well as the single-integral approaches integrate modified superposition concepts. Green and Rivlin 16 developed an integral approach, in which the stress-strain relation is expressed as a sum of multiple integrals. Leaderman 17 developed a nonlinear viscoelastic single-integral model as an extension of the linear Boltzmann-superposition. Applications of these models can be found for different polymeric materials. 18,19 However, there is no distinct standard approach for modeling nonlinear viscoelasticity.
The time-dependent stress response of polymers under a sinusoidal strain becomes nonlinear under sufficiently large strain amplitudes. 20 Nonlinearity here means that the stress response contains beside sine and cosine terms at the fundamental angular frequency, ω, also higher order terms at integer multiples of ω, the so-called higher harmonics. Consequently, if the strain is a sine, the stress is not only a sine shifted by the phase angle. The storage E 0 and loss E 00 moduli lose in parts their physical meaning as characteristic linear material properties. 20,21 The nonlinear contributions can be detected and quantified with high precision via Fourier transform analysis of the stress.
In this work, we first motivate the necessity of nonlinear viscoelastic modeling of thermoplastics based on isochronous stress-strain curves for neat Polyamide 6 (PA6, Firestone, CLM200-001). Following Section 2, the nonlinear viscoelastic material model based on the Schapery model is introduced and the nonlinear approach for the stress-strain relation and the evolution equation are presented in Section 3. Onedimensional elongational strains are considered. A comparison of the proposed nonlinear viscoelastic approach with experimental results for the storage and loss modulus as well as for higher harmonics is presented and discussed in Section 4, followed by concluding remarks in Section 5.

| ISOCHRONOUS STRESS-STRAIN CURVES
The question arises whether the material response of a polymer under certain loads is linear or nonlinear. In terms of a linear constitutive response, the material response is independent of the applied strain load. If the viscoelastic material response depends not only on time but also on the load applied, the material behaves nonlinearly. The nonlinear viscoelastic behavior of polymers can be determined by creep (relaxation) tests at different stress (strain) levels. 15 The stress is then depicted over the strain at different times. When the so-called isochronous stress-strain curves at different times are linear, the material behaves viscoelastic linearly. In contrast, when the isochrones begin to deviate from linearity at a certain stress level, the material shows a nonlinear viscoelastic behavior.
In this work, creep tests at different stress levels are performed for Polyamide 6 (PA6) using the dynamicmechanical analysis (DMA) testing device GABO Eplexor ® 500 N. In Figure 1, the isochrones for PA6 at f gMPa are depicted. From this diagram, the isochronous stress-strain curves for a specific time leave the linear slopes for increasing stresses indicating a nonlinear viscoelastic response.

| Schapery model
For linear viscoelastic material behavior, the time-dependent stress response σ t ð Þ to a uniaxial strain excitation ε t ð Þ is given by a convolution integral over the viscoelastic relaxation modulus and the derivative of the strain 22 time-dependent viscoelastic response. This relation is known as the Boltzmann-superposition. An extension of the Boltzmann-superposition for nonlinear material behavior can be made using the Schapery model that has been developed for solid viscoelastic polymers. The Schapery model introduces straindependent material properties, 22 ξ t ð Þ and ξ 0 τ ð Þ are reduced times given by The functions h 0 ε t ð Þ ð Þ and h 1 ε t ð Þ ð Þ affect the elastic and viscoelastic parts of the stress response. The function h 2 ε t ð Þ ð Þ modifies the strain excitation in the convolution integral. The reduced time ξ t ð Þ allows time-stresssuperposition with strain-dependent shift factors a ε ε t ð Þ ð Þ. If all of the strain-dependent properties are equal to unity, the linear Boltzmann-superposition in the form of Equation (1) is regained. All of the above equations can also be stated in a strain-explicit form with stressdependent properties.
An equivalent formulation of the Schapery model is given by Banks et al. 23 Instead of an integral equation, time-dependent internal variables as a measure of stress are used. The function ΔE t ð Þ is assumed to be of the form which is a frequent approach for the relaxation modulus in linear viscoelasticity. Then, the stress-strain relation can also be stated as The viscoelastic stress behavior is calculated from a number of internal variables ε eα . Every internal variable follows an evolution equation, which in general is an ordinary differential equation with a nonlinear dependence on the strain excitation. Like in Equation (2), nonlinearity is represented by the strain-dependent functions h 0 ε t ð Þ ð Þ, h 1 ε t ð Þ ð Þ, and h 2 ε t ð Þ ð Þ. Time-stress-superposition can be considered by changing the relaxation times in the evolution equations along with the strain. For further generalization, the strain-dependent properties h 1 ε t ð Þ ð Þ, h 2 ε t ð Þ ð Þ, and the shift factors a ε ε t ð Þ ð Þ can be set separately for each internal variable as h 1α ε t ð Þ ð Þ, h 2α ε t ð Þ ð Þ, and a εα ε t ð Þ ð Þ. However, only in case of an equal choice of the straindependent functions for all elements, the equivalence between the integral and the differential formulation is ensured. Setting all strain-dependent functions to unity will lead to a linear model.

| Nonlinear approach
To describe the nonlinear viscoelastic behavior for a specific material, the stress-or strain-dependent properties need to be chosen appropriately. These properties are often described by discrete values for a limited number of stress-or strain-levels, for instance in creep tests. 18 By choosing analytical functions for the nonlinear properties instead, a model for arbitrary load cases can be obtained. In the following, only strain-dependent properties will be considered and the model will later be compared to DMA data for a single testing frequency. Consequently, time-stress-superposition can be neglected and the number of internal variables is set to one. Furthermore, the function h 1 ε t ð Þ ð Þ is assumed to be unity. For the two remaining strain-dependent properties h 0 ε t ð Þ ð Þ and h 2 ε t ð Þ ð Þ, several requirements need to be fulfilled. First, the functions need to reproduce the relevant nonlinear phenomena, that is, higher harmonic oscillations in the stress response for a DMA test. Next, both functions need to reach the asymptotic limit of the linear case for small strain amplitudes. The number of parameters for describing the nonlinearity should be kept as low as possible. Especially for the function h 2 ε t ð Þ ð Þ, a simple expression is desirable. Then, for some load cases such as relaxation, constant strain rate tests, or harmonic strain excitation, the evolution equation in Equation (6) can be solved analytically, which simplifies the parameter identification. Considering these requirements, h 0 ε t ð Þ ð Þ and h 2 ε t ð Þ ð Þ are chosen as with α, β, γ∈. The elastic part of the stress response contains linear and cubic terms in the strain. The right hand side of the differential evolution equation is a polynomial of degree three. The stress-strain relation and the evolution equation follow as The stress can be calculated by inserting the respective strain excitation into Equations (8) and (9) and solving the differential equation with appropriate boundary conditions.

| Dynamic strain input
If the strain is a harmonic oscillation of the form the stress response of the proposed nonlinear model is composed of a mean stress, a phase-shifted oscillation of the base frequency, and higher order harmonic oscillations. The stress response in steady state reads The nonlinear stress response can be written as a sum of the mean stress and in-and out-ofphase components or using intensities I n and phase angles δ n . 20,24 Using the cubic polynomial approach for the nonlinear properties, N ¼ 3 follows. Higher harmonics of two or three times the fundamental angular frequency ω occur. The coefficients for the oscillatory parts are obtained from the solution of the evolution Equation (9) with the strain input following Equation (10) The coefficients contain strain-independent contributions from the linear part of the model as well as straindependent nonlinear contributions. Additionally, the coefficients can be subdivided into contributions with or without a dependence on the loading frequency. The mean stress can be calculated directly. A mean elastic modulus E m is defined as the ratio of mean stress to mean strain For small mean strains, E m transitions to E 0 . If both mean strain and strain amplitude are small, the whole model reduces to the linear formulation.
Considering cyclic strain loads with controlled tension/tension (T=T) mode, the constant static strain ε m is larger than the strain amplitude ε 0 of the dynamic load. The ratio between the minimum ε min and the maximum ε max strain is described by the deformation ratio R. The oscillatory time-dependent elongational strain ε follows Equation (10). The complete stress response is represented by a Fourier series with N ! ∞. The magnitude of the generated higher harmonics (I n , n > 1) is typically normalized by the fundamental (I 1 ) and reported as I n =I 1 = I n=1 . Experimental data reveal, that the contributions I 2=1 and I 3=1 strongly predominate. 25 The normalized intensities of the higher harmonics should scale as power-law functions of ε 0 with exponents of n À 1. For example, I 2=1 and I 3=1 are expected to scale linearly and quadratically with respect to ε 0 , as observed and mathematically derived for suspensions, polymer melts, and solids under shear/torsion loads. 26,27 Experimentally, the higher harmonics can be detected with high sensitivity using high data acquisition rates and oversampling 28 as well as Fourier transform for the data analysis, resulting in the decomposition of a nonlinear stress waveform in the time domain into a Fourier series in the frequency domain. 20,21,24

| APPLICATION OF NONLINEAR VISCOELASTIC APPROACH TO EXPERIMENTAL DATA
The mechanical nonlinear behavior of PA6 under T=T mode was experimentally investigated in a previous publication, 25 revealing the detection and quantification of nonlinear contributions in the stress via Fourier transform. These experimental data are compared with the theoretical predictions of the developed model. For the model, R, ε 0 , and ε m are chosen according to the experiments. To identify an optimal set of parameters, a gradient-based nonlinear optimization with the experimental data and the coefficients from Equations (13) and (14) is conducted. The mean elastic modulus E m , the storage modulus E 0 , the loss modulus E 00 , and the intensities of the second and third order higher harmonics, I 2=1 and I 3=1 , are compared. The results are shown in Figures 2 and 3.
For small strain amplitudes up to approximately 0:3%, the material behavior is predominantly characterized by the linear viscoelastic properties E 0 and E 00 . In the experimental data, the loss modulus E 00 is constant, but the storage modulus E 0 significantly decreases for small strain amplitudes. As in the model both E 0 and E 00 have slightly decreasing values for β < 0 and small strain amplitudes, ε 0 ! 0, the predictions differ from the experiment. The nonlinear higher harmonics at those strain amplitudes are below the noise level of the machine and cannot be quantified, resulting in noisy experimental data, whereas very small I n=1 values < 10 À3 ð Þ are predicted by the model. For medium strain amplitudes, the nonlinear higher harmonic contributions increase above the noise threshold of the machine, so they can be measured and quantified. The intensity I 2=1 increases linearly with the strain amplitude and the intensity I 3=1 scales quadratically with the strain amplitude, both in the model and the experimental data. The model predicts the increase of I 2=1 well, although there are several outliers in the experimental data causing a fluctuating error. However, I 3=1 is estimated to be too small by a factor of 10 approximately. Similar I 2=1 and I 3=1 trends were found for the analysis of hyperelastic material models, such as the Neo-Hooke, Mooney Rivlin, Ogden, or Arruda-Boyce model, where the predicted I 2=1 fits well to experimental data and to the prediction of the here presented model, but underestimates I 3=1 by over a factor of 10. 29 The storage modulus E 0 decreases and the loss modulus E 00 increases in the experiments. Both trends are predicted by the model with an relative error of less than 25% for E 0 and E 00 , and an average error of approximately 11%, respectively. The mean strain increases proportionally with the strain amplitude. The mean elastic modulus E m over the entire strain range behaves similar to E 0 in the experimental data. The model predicts an initially constant values for E m and a decrease for larger strains. The maximum relative error is also at approximately 20%.

| CONCLUSIONS
For comparing a viscoelastic material model with experimental data, typically only moduli and/or the complete signal are available and nonlinear effects can hardly be determined. Simply fitting the theoretical stress curves into experimental data is a low-sensitive and low-key approach to detect and compare nonlinear contributions. Advances in experimental techniques allow to fully quantify with high sensitivity (signal to noise ratio in the range of 10 4 : 1 to 10 5 : 1) 20,24 a nonlinear stress response via Fourier transform, such that even small nonlinear contributions can be quantified and compared to theoretical predictions. Thus, the approach to Fourier transform the stress response of experimental and theoretically calculated data closes this gap.
In the work at hand, the nonlinear mechanical behavior of viscoelastic materials is addressed. Experimental data from creep tests reveal nonlinear effects for stresses greater than approximately 8 MPa. A mechanical model based on Schapery's theory in a stress-explicit formulation by Banks et al. 23 with a new approach for the strain-dependent material properties is proposed and formulated for the one-dimensional case. The model is compared with experimental data of the Fourier-transformed stress response for sinusoidal strains. Due to the simplicity of the approach, the model equations can be solved analytically for this loadcase. The resulting nonlinear stress response contains all experimentally identified contributions and correctly predicts their qualitative strain dependencies for strains up to 3%. With the optimal parameters set, the model can well-describe I 2=1 , but highly underestimates I 3=1 . A different, more complex, approach for the straindependent material properties could enhance the quality of the model predictions. This would also make solving the equations and identifying parameters more difficult, though. The applicability of the model and the transferability of its parameters considering other load cases has not yet been investigated, here. Nevertheless, to consider other load cases, such as relaxation or creep, the strain-dependent functions in the nonlinear approach can be adapted if needed.