Characterising the linear viscoelastic behaviour of an injection moulding grade polypropylene polymer

Abstract

The linear viscoelastic behaviour of an injection moulding grade polypropylene is studied using theoretical and computational methods. Polypropylene has a variety of engineering applications as a component. However, it commonly exhibits viscoelastic deformations. This paper analyses the creep and recovery responses of the BJ368MO polypropylene copolymer using the Burgers and generalised Maxwell models. Within the linear viscoelastic regime, an experimental creep strain at $20\ \text{MPa}$ is used to determine the rheological constants of the models. These constants (springs and dashpots) are determined using a nonlinear least-squares curve fitting of the experimental creep. Then they are used to predict the creep and recovery responses of the polymer at three different stresses, $10\ \text{MPa}$, $12.5\ \text{MPa}$ and $15\ \text{MPa}$. The experiments are made using tensile specimens designed according to the ASTM D638-14standard. The theoretical evaluations are made using the creep and recovery equations derived from their constitutive. Whereas COMSOL Multiphysics software is used during the finite element (FE) analyses. The results of the theoretical and FE calculations are verified using creep and recovery experiments. Based on the validation analyses, both viscoelastic models showed lower deviations from the experimental results when a computational approach is used. In addition, the viscoelastic models are compared by evaluating the residuals of the creep and recovery strain predictions. The theoretical analyses showed better predictions at $12.5\ \text{MPa}$ and $15\ \text{MPa}$ stresses when the generalised Maxwell model is used. However, the improvements are attributed to the recovery predictions. When FE is used, the Burgers model showed lower mean absolute percentage errors (MAPEs) in all creep and recovery predictions. The model has a minimum of 6.37% error at the $10\ \text{MPa}$ stress and a maximum of 8.23% error at the $15\ \text{MPa}$. By comparison, the generalised Maxwell model showed a minimum of 9.24% error at $12.5\ \text{MPa}$ and a maximum of 12.8% error at $15\ \text{MPa}$ stresses. The novelty of this paper is on predicting the creep and recovery behaviour of the polymer using the FE and theoretical approaches in the linear viscoelastic regime. The findings suggest that the FE analyses using the Burgers viscoelastic material model provide better predictions, with all calculated errors falling below 10%.

1 Introduction

Polypropylene is a widely used polymer in engineering applications for its attractive physical properties including, low density, good toughness, and better fatigue resistance under elevated temperatures (Harutun 2003). These excellent properties make polypropylene widely versatile for automotive interior, medical equipment, fluid transport systems, and commodity applications (Moor 1996). Along with other smart and state-of-the-art materials, polypropylene is currently a biomedical engineering solution as drug delivery systems and implantable devices (Maddah 2016). The state-of-the-art in sustainable materials development focusses on the synthesis of polypropylene based bio-composites and bio-nanocomposites. This has enhanced the mechanical properties of the composites, increasing their applicability in automotive industries, structural reinforcements, and packaging (Visakh 2017). However, due to the broad range of applications, the polymer experiences stress due to static, dynamic, or thermal loadings leading to time-dependent deformations. Therefore, understanding the creep and recovery behaviour of the polymer is essential in order to avoid failures caused by different loadings.

There are a few studies that discuss the linear viscoelastic behaviour of engineering polymers and their composites. Daver et al. (2016) used creep-recovery experiments to characterise the viscoelastic responses of polyolefin-rubber nanocomposites developed for additive manufacturing. Based on their experimental results, the authors used analytical and numerical techniques to model the viscoelastic behaviours using extrapolations and intercepts on creep curves. Martins et al. (2015) used curve-fitting on experimental creep data to determine the rheological constants for a polylactic acid (PLA) polymer blended with polycaprolactone. The authors determined the rheological constants using short-term experimental data. Yang Yang et al. (2006) investigated the creep resistance of polyamide 66 with titanium dioxide nanofillers. The long-term creep behaviours of Burgers and Findley power law material models are determined using the nonlinear curve fitting function on OriginPro. They set the initial values of the rheological parameters manually to study the simulated curve asymptote to their experimental results. Dean and Broughton (Dean and Broughton 2007) studied the influence of retardation time parameters casting nonlinear viscoelasticity on rectangular polypropylene samples under uniaxial tension and compression. Similarly, the creep and relaxation behaviours of polypropylene are studied using finite element analyses (Resan et al. 2015). In many cases, the viscoelastic parameters that are determined geometrically from the creep curve by taking strain rates at $t\ =\ 0$ and $t\ =\infty$ lack accuracy. The linear viscoelastic parameter analyses are discussed for different materials (Keenan et al. 2013; Ngudiyono et al. 2019; Bharadwaj et al. 2017). The FE studies by Ngudiyono et al. (2019) and Bharadwaj et al. (2017) used the Prony series to model the time-dependent shear moduli of the materials.

The viscoelastic behaviour of polypropylene has been studied by taking into account the different morphological compositions or additives for reinforcing specific mechanical properties. Wang et al. (2018), have analysed the creep and recovery properties of injection moulded isotactic polypropylene (iPP). They used the Burgers model for fitting the creep behaviour and a Weibull distribution function for the recovery. Kurt and Kasgoz (2021) have studied the effect of different molecular weight and their distributions on the creep behaviour of the same composite. They used the four parameters Burgers model to fit the creep. The effects of adding different crystals and porosities have been studied by Wu et al. (2019). Furthermore, several reinforced polypropylene composites have been studied lately using both experimental and complementary theoretical modelling, including materials such as wood (Lee et al. 2004; Homkhiew et al. 2013), silica (Rosa et al. 2018), glass fibre (Berecz et al. 2012) and e-glass (Zhai et al. 2018). The creep behaviour of polypropylene blended with other polymers Dian et al. (2020), Houshyar et al. (2005), Banik et al. (2007), Martins et al. (2015) and organic materials such as fibres (Militký and Jabbar 2015; Tiwari et al. 2021; Hao et al. 2014) and clay (Drozdova et al. 2009) are studied using experimental, theoretical and computational methods.

Traditionally, rheological constants that characterise viscoelastic deformations are determined by finding the slopes at instantaneous elasticity, long-term creep, and recovery. However, it is impossible to find the tangents to creep curves (Findley et al. 1976; Riande et al. 2000). This paper aims to predict the creep and recovery behaviour of the BJ368MO polypropylene copolymer at different stresses within the linear viscoelastic regime. Two comparable viscoelastic models, the Burgers and the generalised Maxwell models, will be used to estimate the creep and recovery behaviour of the polymer theoretically and computationally. In order to predict the behaviour, the rheological constants of the viscoelastic models will be determined first by fitting creep equations to experimental data. The regression analyses will be made using the nonlinear least-squares curve fitting on MATLAB. Finally, the creep and recovery estimations will be compared with the experimental findings to identify the approach and viscoelastic model that predict the phenomena better. The flowchart in Fig. 1 summarises the methodology we followed in this paper.

Fig. 1

The flowchart of methodology used to determine the suitable approach and model for predicting the linear viscoelastic behaviour of BJ368MO polypropylene copolymer

2 Linear viscoelasticity

2.1 Linear viscoelastic models

Engineering components made from polymers are widely used in constant stress environments. Even if under constant loading, the components experience time-dependent deformation due to their viscoelastic properties. This time-dependent deformation is referred to as creep (Findley et al. 1976). The linear viscoelasticity approach can be used to analyse the strain response of viscoelastic materials under constant stresses. This strain response can be modelled using constitutive equations which govern the exhibited viscous and elastic natures (Riande et al. 2000). For a linear viscoelastic material, the stress-strain relations in the elastic (spring) and dashpot components are presented as (Brazel and Rosen 2012)

$$ \begin{aligned}[b] &\sigma =E\varepsilon ,\\ &\sigma = \eta \dot{\varepsilon }. \end{aligned} $$

(1)

$E$ and $\eta$ are the elastic and dashpot constants that are generally referred to as rheological constants. Analytically, the strain responses of polymers are evaluated using a variety of spring-dashpot configurations. Some of the basic viscoelastic models are Maxwell, Kelvin/Voigt, or Zener (Standard Linear Solid) (Gutierrez-Lemini 2014; Kelly 2013). However, these models are inadequate to predict the materials’ creep, recovery, or relaxation behaviours under constant or step loading profiles (Crawford and Martin 2020). A better viscoelastic representation can be made using complex viscoelastic models. In this regard, the Burgers and generalised Maxwell models are mostly used to analyse the creep and recovery behaviour of the materials (Maddah 2016).

The stress-strain constitutive equation of a linearly viscoelastic material under constant loading in shear or uniaxial tension is presented as (Findley et al. 1976)

$$ \sum _{i=0}^{n} p_{r} \frac{\partial ^{r}}{\partial t^{r}} \sigma = \sum _{i=0}^{m} q_{r} \frac{\partial ^{r}}{\partial t^{r}} \varepsilon . $$

(2)

The coefficients, $p_{r}$ and $q_{r}$, are varying combinations of material’s rheological constants. The subscripts $m$ and $n$ are interrelated and represent the number of combinations of the viscoelastic models, whereas $r$ denotes the order of differential. Equation (2) is the constitutive equation that can be modified to fit viscoelastic models, such as Burgers and generalised Kelvin-Voigt.

The Burgers viscoelastic model has instantaneous and delayed elasticities in which a single Maxwell model is connected in series with one Kelvin-Voigt; see Fig. 1. As it is a four-parameter model, the constitutive equation of the Burgers model can be given as

$$ p_{o} \sigma + p_{1} \dot{\sigma } + p_{2} \ddot{\sigma } = q_{1} \dot{\varepsilon } + q_{2} \ddot{\varepsilon }. $$

(3)

The creep of the Burgers viscoelastic model under a step loading can be derived from Eq. (3). The step loading is given using the Heaviside function as

$$ \sigma = \sigma _{o} H \left ( t \right ). $$

(4)

The total strain of the Burgers viscoelastic model presented in Fig. 2 is given as the algebraic sum of strains from the individual element.

Fig. 2

The Burgers viscoelastic model represented with a Kelvin and a Maxwell model connected in series

The total strain is

$$ \varepsilon = \varepsilon _{1} + \varepsilon _{2} + \varepsilon _{3}. $$

(5)

The strains can be presented as functions of the material’s rheological constants and applied stress. That is,

$$ \varepsilon _{1} = \frac{\sigma }{E_{1}},\qquad \dot{\varepsilon }_{2} = \frac{\sigma }{\eta _{1}},\qquad \dot{\varepsilon }_{3} + \frac{E_{2}}{\eta _{2}} \varepsilon _{3} = \frac{\sigma }{\eta _{2}}. $$

(6)

The Laplace transform of stress-strain relations in Eq. (6) can be replaced in Eq. (5) to determine the coefficients of the constitutive form presented in Eq. (3). With the substitutions, Eq. (3) becomes

$$ \sigma + \frac{1}{E_{1} E_{2}} \left ( \eta _{1} E_{2} + \eta _{1} E_{1} + \eta _{2} E_{1} \right ) \dot{\sigma } + \frac{\eta _{1} \eta _{2}}{E_{1} E_{2}} \ddot{\sigma } = \eta _{1} \dot{\varepsilon } + \frac{\eta _{1} \eta _{2}}{E_{2}} \ddot{\varepsilon }. $$

(7)

The strain response of the Burgers model is determined by solving the above second-order differential equation using the appropriate stress and strain initial conditions (Appendix A). The creep of the Burgers viscoelastic model is given as

$$ \varepsilon \left ( t \right ) = \sigma _{o} \left [ \frac{1}{E_{1}} + \frac{1}{\eta _{1}} t+ \frac{1}{E_{2}} \left ( 1- e^{- \frac{E_{2}}{\eta _{2}} t} \right ) \right ]. $$

(8)

In Eq. (8), the first term on the right-hand side is the instantaneous elasticity, whereas the last one shows the delayed elasticity.

The generalised Maxwell model is another viscoelastic representation of materials with a number of Maxwell models connected in parallel; see Fig. 3.

Fig. 3

The generalised Maxwell viscoelastic model with Maxwell branches connected in parallel

The viscous (fluid) or elastic (solid) behaviours of a viscoelastic material are alternatively represented by connecting an isolated dashpot or spring to the generalised Maxwell model. The generalised Maxwell model presented in Fig. 3 shows the instantaneous and delayed elasticities with multiple retardation times. The stress-strain relation for the generalised Maxwell model is given as (Gutierrez-Lemini 2014)

$$ \sigma = \left [ \sum _{i=1}^{N} \frac{E_{i} \eta _{i} \frac{\partial }{\partial t}}{E_{i} + \eta _{i} \frac{\partial }{\partial t}} \right ] \varepsilon . $$

(9)

Equation (9) is further manipulated to get the constitutive in a standard form. This can be approached by multiplying both sides of Eq. (9) by the minimum common denominator of the right-hand side expression. This leads to

$$ \left [ \prod _{i=1}^{N} \left ( E_{i} + \eta _{i} \frac{\partial }{\partial t} \right ) \right ] \sigma = \left [ \sum _{i=1}^{N} \left ( E_{i} \eta _{i} \frac{\partial }{\partial t} \right ) \prod _{j=1}^{N} \left ( E_{j} + \eta _{j} \frac{\partial }{\partial t} \right ) \right ] \varepsilon ,\quad i\neq j. $$

(10)

For a four-parameter model, the above constitutive equation reduces to

$$ \sigma + \left ( \frac{\eta _{1}}{E_{1}} + \frac{\eta _{2}}{E_{2}} \right ) \dot{\sigma } + \frac{\eta _{1} \eta _{2}}{E_{1} E_{2}} \ddot{\sigma } = \left ( \eta _{1} + \eta _{2} \right ) \dot{\varepsilon } + \eta _{1} \eta _{2} \left ( \frac{1}{E_{1}} + \frac{1}{E_{2}} \right ) \ddot{\varepsilon }. $$

(11)

The second-order differential form presented in Eq. (11) is solved using the Laplace transform with appropriate initial conditions (Appendix A). The creep of the four-element generalised Maxwell model is

$$ \varepsilon \left ( t \right ) = \sigma _{0} \left [ \frac{t}{\eta _{1} + \eta _{2}} + \frac{e^{-at}}{E_{1} + E_{2}} + \frac{1}{\eta _{1} + \eta _{2}} \left ( \frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}} - \frac{1}{a} \right ) \left ( 1- e^{-at} \right ) \right ]. $$

(12)

In Eq. (12), the coefficient $a$ is the inverse of the retardation time, $a= \frac{\left ( \eta _{1} + \eta _{2} \right ) E_{1} E_{2}}{\eta _{1} \eta _{2} \left ( E_{1} + E_{2} \right )}$.

2.2 Viscoelastic material model on COMSOL

Nowadays, FE software provides possibilities to define mechanical behaviours of materials that contribute to accurate material models in computational analyses. Our previous work on the influence of stiffener geometry on flexural properties of additive manufactured beams used the COMSOL structural mechanics for modelling and analyses (Silas et al. 2020). COMSOL has suitable features to define viscoelastic material models, such as generalised Maxwell, Burgers and standard linear solids (SLS). In the FE calculations, the rheological constants are used in the governing equations. For the generalised Maxwell model, the viscoelastic deformation is treated using constitutive equations of the individual Maxwell branches (COMSOL 2019). The relationship in terms of the total strain, $\epsilon$, the $m$th branch strain, $\epsilon _{m}$ and its strain rate, $\dot{\epsilon }_{m}$ is

$$ \tau _{m} \dot{\epsilon }_{m} + \epsilon _{m} =\epsilon . $$

(13)

In Eq. (13), $\tau _{m}$ is a relaxation time measured in the frequency domain and is further given as the ratio of

$$ \tau _{m} = \frac{\eta _{m}}{G_{m}}. $$

(14)

In Eq. (14), $\eta_{m}$ is the dashpot, and $G_{m}$ is the stiffness of the branch. The viscoelastic stress of the generalised Maxwell model, $\sigma _{q}$, is the sum of the stresses in each Maxwell branch, $\sigma _{m}$. That is,

$$ \sigma _{q} = \sum _{m=1}^{N} \sigma _{m}. $$

(15)

Similarly, a second-order ODE is used as a governing constitutive equation to relate the stress and strain tensors in the Burgers model (COMSOL 2019).

For a viscoelastic material model under normal loading, the stiffnesses are defined using elastic moduli instead of the shear ones (Lopes et al. 2014). In addition, in linear viscoelasticity, Hooke’s law can be modified to include the viscoelastic stresses to govern the behaviour of materials (COMSOL 2019; Ward and Sweeney 2004).

2.3 Curve fitting

The curve fitting function in MATLAB is used to perform the regression analysis. The fit function can be used to perform linear or nonlinear regressions to fit curves to experimental data. A parametric curve fitting to specific data is also possible by introducing custom equations that define the fit. The curve fitting accuracy can be enhanced by assessing and improving the goodness-of-fit and confidence intervals. The code for the curve fitting function can be automatically generated and exported to a workstation for the manipulations (MATLAB 2004). In this paper, MATLAB codes will be used for the nonlinear least-square fits. The regression analyses will be made on the experimental creep results to determine the rheological constants of the two viscoelastic models.

3 Experimental methods

The creep experiments are made using injection moulded tensile test specimens prepared from the BJ368MO polypropylene. The specimens were produced on IntElect 100 injection moulding machine that has the state-of-the-art control features for high precision manufacturing. During the sample manufacturing, the nozzle temperature was kept at 220 °C. The temperatures at the different zones and the mould were set as; zone $3 = 220$ °C, zone $2 = 215$ °C, zone $1= 210$ °C, $\text{feed} = 30$ °C and $\text{mould} = 30$ °C. These temperature profiles are digitally controlled with tight tolerances of ±1 °C. The injection speed is set at $13\ \text{mm}/\text{s}$. For this injection speed, the injection pressure readings range between the maximum and minimum values of $500\ \text{bar}$ and $183\ \text{bar}$. The holding switchover is controlled via screw position. When the screw stroke reaches $8\ \text{ccm}$, the injection profile changes to packing by applying a constant pressure of $300\ \text{bar}$ for 3 seconds. After packing, the screw plasticises at $50\ \text{rpm}$ until it reaches the screw-back position of $32\ \text{ccm}$. A two-cavity mould is used to manufacture the samples with a cycle time of 18 seconds.

The test specimens are designed and manufactured according to the ASTM D638-14 standard for type I samples (ASTM D638-14 2014). The creep and recovery experiments are made on recently developed X-series high precision and digitally controlled tester by Testometric. The X350-20 tester has a controllable test speed range between 0.0001–1000 mm/s. The strain is carefully measured using the TE3442 digital extensometer acquired from Epsilon technology. Figure 4 shows the sample gripping and extensometer attachment during testing.

Fig. 4

Experimental setup of the creep and recovery experiments. The tests are made using the X350-20 machine from Testometric. The strains are measured using higher precision clip-on extensometer, TE3542

In order to apply the theoretical and computational methods, the creep and recovery responses should be within the linear viscoelastic regime. The linear viscoelastic creep and recovery responses of polypropylene composites have been experimentally studied in the past (Wang et al. 2018; Kurt and Kasgoz 2021; Wu et al. 2019; Lee et al. 2004; Homkhiew et al. 2013; Rosa et al. 2018). However, the experimental studies vary on the testing methods, the equipment used, and the testing conditions. Using the universal testing machines, uniaxial tension (Wu et al. 2019; Lee et al. 2004; Zhai et al. 2018; Houshyar et al. 2005; Drozdova et al. 2009) and flexural (Militký and Jabbar 2015; Tiwari et al. 2021) tests are made to study creep of polypropylene composites. Nano-indentation tests (Rosa et al. 2018; Dian et al. 2020), a dynamic oscillatory rheometer (Kurt and Kasgoz 2021) and a dynamic mechanical analyser (Wang et al. 2018; Houshyar et al. 2005) are used to study the effects of chemical modification on the mechanical performance of polypropylene. Similarly, the creep and recovery behaviours of the polypropylene composites are studied at different temperatures (Wang et al. 2018; Banik et al. 2007; Militký and Jabbar 2015). The different experimental approaches followed by the different authors have produced a wide range of values for the rheological constants.

We conducted several creep and recovery experiments before the curve fitting and the viscoelastic evaluations. The material’s viscoelastic responses are experimentally investigated at different stresses ranging between $10\ \text{MPa}$ and $35\ \text{MPa}$. Figure 5 presents the linear and nonlinear viscoelastic responses observed during the creep and recovery experiments. The nonlinear viscoelastic responses are exhibited at higher stresses, leading to failures at early stage of the creep. Thus, the study focusses on the responses at $20\ \text{MPa}$ and below. The creep response that is used for the curve fitting is made at the $20\ \text{MPa}$ loading; see Fig. 5. Each creep experiment is made for a length of time, $t=3600$ seconds. During each experiment, the load function was defined using the step tensile testing method. The machine is programmed with the load sequence defined as $\sigma =20\ \text{MPa}$ at $t=\ 0$ and $\sigma =\ -20\ \text{MPa}$ at $t=\ 3600$ seconds. Practically, it is not possible to apply the instantaneous stress changes during the experiments. Therefore, we used a higher travel speed of $350\ \text{mm}/s$ during initial loading in the step load sequence. As a result, it took below 0.5 seconds for the tester to reach the constant stress. This stress was maintained with a $25\ \text{mm}/s$ holding speed until the end of the creep experiment.

Fig. 5

Experimentally studied creep and recovery responses of the BJ368MO polypropylene at different stresses. The creep load function is defined for $\sigma =20\ \mbox{MPa}$ whereas, the creep & recovery load function is defined for the other stresses

The initial loading conditions of the creep and recovery experiments for the other stresses are defined similarly to the $20\ \text{MPa}$ creep test. During the recovery, the loads are removed, and the strain recoveries are studied for 1200 seconds. Each experiment is repeated five times, and the arithmetic mean of the responses is used for the viscoelastic analyses.

3.1 Curve fitting for rheological parameters

The rheological constants of the Burgers and generalised Maxwell models are determined using a custom equation fit on MATLAB. The creep of the two viscoelastic models presented in Eqs. (8) and (12) are used for fitting the experimental data. Figure 6 (a) and (b) present the fittings made that correlate the models with the experiment. The nonlinear least-squares fit showed good agreement with acceptable goodness-of-fit characteristics. The details of these fits are given. See Appendix B.

Fig. 6

Curve fitting of the linear viscoelastic models on the experimental creep obtained at stress, $\sigma =20\ \mbox{Mpa}$ a) Burgers model fit. b) Generalised Maxwell model fit

The values of the rheological constants determined for the viscoelastic models are presented in Table 1.

Table 1 Rheological constants determined by curve fitting

3.2 Viscoelastic modelling analyses on COMSOL

Based on the rheological constants presented in Table 1, the material’s creep and recovery behaviour is studied computationally using COMSOL. The responses of the Burgers and generalised Maxwell viscoelastic material models are studied for the $10\ \text{MPa}$, $12.5\ \text{MPa}$ and $15\ \text{MPa}$ stress loadings. From the solid mechanics physics, the Burgers and generalised Maxwell viscoelastic material models are selected during the analyses. The rheological constants and other physical properties from the manufacturer’s data are used to define the material model (Borealis 2019). The loads are defined using a piecewise function, which follows the experimental step tensile load sequence. The step loading introduces nonlinearity leading to nonconverging solutions due to the residuals between successive iterations. In order to avoid this, the piecewise load function is smoothed using a continuous second-order derivative. The modelling and analyses started with building the ASTM D 368-14 type I sample using the geometry tool on the 3D space dimension. Tetrahedral mesh with a user-controlled meshing sequence is assigned to the model. The minimum element size is set at $0.1\ \text{mm}$, which generated closer to half a million degrees of freedom to be solved. In each study, boundary loads are applied for the first 3600 seconds and removed during the next 1200 seconds. Time-dependent studies are conducted, and the mesh refinements are approached by sequentially decreasing the element size until the solutions showed convergences. Figure 7 shows the mesh and boundary settings of the FE model.

Fig. 7

Finite element modelling and boundary conditions. a) Tetrahedral mesh. b) Boundary load in a uniaxial direction

4 Results

The linear viscoelastic deformation of the BJ368MO polypropylene copolymer is studied using experimental, computational and theoretical methods. The rheological constants determined via regression analyses are used to evaluate the creep and recovery of the material at different stresses. The responses are studied using theoretical, and FE approaches. The Burgers and generalised Maxwell viscoelastic models are alternatively assigned to the material during the analyses. The results of the theoretical and FE evaluations are compared to the experimental findings.

4.1 Theoretical results

The creep and recovery responses of the two viscoelastic models are theoretically calculated. Employing the superposition principle, Eqs. (8) and (12) are modified to evaluate the recovery strain after the load removal (Riande et al. 2000). The evaluations are made for the $10\ \text{MPa}$, $12.5\ \text{MPa}$ and $15\ \text{MPa}$ stresses. Figure 8 (a) – (c) present the results of the creep and recovery evaluations compared to the experimental ones.

Fig. 8

Theoretically calculated creep and recovery strains of the Burgers and generalised Maxwell models presented in comparisons with the experimental results at a) $\sigma _{o} =10\ \mbox{MPa}$, b) $\sigma _{o} =12.5\ \mbox{MPa}$, c) $\sigma _{o} =15\ \mbox{MPa}$

The instantaneous deformation of the material is closely similar in all cases. However, with the increase in time, the creep predictions by the models deviated from the experimental results. In contrast, the recovery predictions of the models showed variations for each loading. Generally, the rheological constants can be attributed to the deviations observed. The creep and recovery of the two viscoelastic models have constant, linear and exponential terms with different combinations of rheological constants. The retardation times of the two models have similar orders of magnitude. Therefore, the exponent terms of the models in Eqs. (8) and (12) behave similarly during the creep. However, at each viscoelastic branch, the dashpots of the Burgers model are higher than the generalised Maxwell ones. These influence the linear terms of the creep strains predicted by the models. Similarly, the recovery strains predicted by the models have terms related to previous load histories and exponentially decaying expressions which depend on the rheological constants. In addition, unlike the practical experimentations, the theoretical models exhibit instantaneous recoveries which are derived from the principles of their constitutive.

Comparing the theoretical estimations, the Burgers model predicts the creep and recovery of the material better at $\sigma _{o} =10\ \text{MPa}$. However, when the stress increases to $12.5\ \text{MPa}$ and $15\ \text{MPa}$, the generalised Maxwell becomes a better model. The improvements of the model’s accuracy arise from the long-term strain recovery. Considering only the creep, the theoretical prediction by the Burgers model is better in all stress loadings.

4.2 Finite element analyses results

The creep and recovery analyses of the polymer are made on COMSOL using the Burgers and the generalised Maxwell viscoelastic material models. The viscoelastic analyses predicted the material’s response for the $10\ \text{MPa}$, $12.5\ \text{MPa}$, and $15\ \text{MPa}$ stress loadings. The time-dependent strains in the loading direction are studied for the length of time, $t=\ 4800$ seconds, with a 1 second increment. From the graphic window, the average surface strains of both viscoelastic models for the $15\ \text{MPa}$ stress are presented as examples; see Figs. 9 (a) and (b).

Fig. 9

Surface strain in the local x-axis for the $15\ \mbox{MPa}$ stress. a) Burgers viscoelastic material model. b) Generalised Maxwell viscoelastic material model

The FE creep and recovery predictions by the two models and the experimental investigations are compared; see Figs. 10 (a) – (c).

Fig. 10

FE creep and recovery predictions of the Burgers and generalised Maxwell models presented in comparison with the experimental findings at, a) $\sigma _{o} =10\ \mbox{MPa}$, b) $\sigma _{o} =12.5\ \text{MPa}$, c) $\sigma _{o} =15\ \text{MPa}$

The creep and recovery results of the FE analyses generally showed similar responses during the instantaneous and long-term creep. However, higher deviations from the experimental results are locally observed at initial loadings. During the creep and recovery experiments, the initial loadings were made at a test speed of $350~\text{mm}/\text{s}$. With this travel speed, the tester reached the constant load sequence within 0.5 seconds in all cases. However, in FE analyses, abrupt changes in boundary load lead to nonlinearly, resulting in nonconvergence (Lopes et al. 2014). In order to resolve this, the piecewise load functions are smoothed at loading and unloading transitions. A second-order continuous in time smoothing function is used to modify the step loadings slightly. The initial creep responses of the FE analyses deviated from the experimental ones due to this load smoothing. During the recovery, the return speed of the tester was reasonably lowered in the experiments. As a result, the recovery of the instantaneous elasticity in both experimental and FE analyses are similar. However, the models exhibited quick recovery predictions for the $12.5\ \text{MPa}$ and $15\ \text{MPa}$ loadings. Comparing the models, the Burgers model showed good agreement with the experimental results in all loadings. The generalised Maxwell model comparatively has higher deviations in creep. However, the recovery predictions of the model are better.

The mean absolute percentage errors (MAPE) of the theoretical and FE analyses are further investigated. These deviation analyses determine the margins of acceptability of the results predicted by the two viscoelastic models. The MAPE is calculated as

$$ \mathit{MAPE}= \frac{1}{n} \sum _{t=0}^{n} \left \vert \frac{\varepsilon _{e x_{t}} - \varepsilon _{t}}{\varepsilon _{ex_{t}}} \right \vert *100\%. $$

(16)

In Eq. (16), $n$ is the number of data points used, $\varepsilon _{ex_{t}}$ is the experimental strain, and $\varepsilon _{t}$ is the calculated strain using either the theoretical or FE methods. Table 2 presents percentage errors of the two viscoelastic models as compared to the experimental approach.

Table 2 Mean absolute percentage errors of the viscoelastic models using analytical and FE

The MAPEs in Table 2 show the deviations of the theoretical and FE strain predictions by the models. The residual analyses indicate that the theoretical creep and recovery estimation by the Burgers model has its smallest deviation of 13.11% for the $10\ \text{MPa}$ stress loading. The prediction errors showed increments when the load changes to $12.5\ \text{MPa}$ and $15\ \text{MPa}$; see Table 2. However, the majority of the model’s prediction errors came from the recovery strain response. Considering only the creep, the model’s prediction error reduced to approximately 7.45% in all loading cases. Unlike the Burgers model, the generalised Maxwell theoretical creep and recovery prediction errors do not show proportional increments with the loading. The lowest prediction error of the model is 12.16% for the $12\ \text{MPa}$ loading. Relatively, the model showed lower prediction errors for the $12.5\ \text{MPa}$ and $15\ \text{MPa}$ loading. The improvements are due to lower deviations during strain recovery. The creep prediction of the model showed a range of MAPE between 8% and 9%.

The deviations of the FE creep and recovery analyses are also evaluated. The percentage error of creep and recovery prediction by the Burgers model is below 10% in all three cases. Like the theoretical analyses, the model exhibited increments in percentage deviation with the load increments. On the other hand, the generalised Maxwell showed a higher percentage deviation in all the cases. This model has a maximum deviation of 12.8% for the $15\ \text{MPa}$ loading. Unlike the Burgers, the percentage deviations of the generalised Maxwell model do not have a linear correlation to the stress increments. With minimum residuals, the FE approach represents the linear viscoelastic behaviour of the polymer better than the theoretical counterpart.

Two viscoelastic models are used to predict the creep and recovery behaviours of the polymer theoretically and computationally. The evaluations are compared with the results of the experimental investigation; see Figs. 8 and 10. The experimental approach is used for validating the analyses by the methods. The MAPEs of all analyses are computed to compare and identify the better approach and model that predicts the viscoelastic behaviour of the material. Based on the results of all analyses, the Burgers viscoelastic model showed good agreement with experimental results when the FE method is used.

5 Conclusion

The linear viscoelastic behaviour of the BJ368MO polypropylene copolymer is studied using theoretical and FE methods. In both methods, the material’s viscoelastic deformation is predicted using the Burgers and generalised Maxwell models. The rheological constants of the two viscoelastic models are determined using nonlinear curve fitting of the experimental creep made at $20\ \text{MPa}$ stress. The determined rheological constants are used in theoretical and FE approaches to predict the creep and recovery of the material at different stresses. Generally, the results of the FE analyses are better than the theoretical ones for both viscoelastic models. Furthermore, the FE creep and recovery prediction using the Burgers viscoelastic model showed minimum percentage errors for all stress loadings.

The comparison between the models’ theoretical creep and recovery results indicated that the Burgers model was better at the $10\ \text{MPa}$ loading. For the $12.5\ \text{MPa}$ and $15\ \text{MPa}$ loadings, the generalised Maxwell model became the better model. If the comparison is made solely on the creep, the Burgers model predicts better, with deviation errors falling close to 7.45% for all stress loadings. For both models, the errors in recovery prediction contributed highly to the increment of the overall prediction errors. One of the reasons for the higher deviations is the models’ instantaneous recovery during load removal. Experimentally, abrupt load removal is not possible.

The results of the FE creep and recovery analyses showed improved predictions by the viscoelastic models. The Burgers model showed the smallest deviation errors falling below 10% in all stress loadings. The maximum deviation error of the model is 8.23% at $15\ \text{MPa}$ stress. The FE calculations using the generalised Maxwell model showed smaller prediction errors than its theoretical counterpart. However, it exhibited higher prediction errors at all stress loadings when compared with the Burgers FE results. The higher prediction errors of the model are attributed to the values of its rheological constants as well as its spring-dashpot layout. Even if both models have four elements, their layouts are different. This affects the material’s creep and recovery predictions. The generalised Maxwell model is suitable for modelling the relaxation behaviour of materials. Using parallel Maxwell branches, materials’ viscoelastic relaxation at different times can be determined (Baumgaertel and Winter 1992; Jalocha et al. 2015).

Characterisation studies of the creep behaviour of polypropylene polymer are made in Dean and Broughton (2007) and Resan et al. (2015). However, the studies generally focus on the creep and relaxation responses using a smaller number of spring-dashpot elements. Kasgoz et al. (2018) studied temperature-dependent creep and relaxation of isotactic polypropylene. However, the study solely focussed on predicting the creep of the polymer using time-temperature superposition. Other works have focussed on the mix of polypropylene with additives such as polymers (Houshyar et al. 2005; Banik et al. 2007; Martins et al. 2015), organic materials (Homkhiew et al. 2013; Militký and Jabbar 2015; Tiwari et al. 2021; Hao et al. 2014) and many others (Drozdova et al. 2009; Loghman and Shayestemoghadam 2016). This work proposed a method for predicting the creep and recovery behaviour of polypropylene using two different viscoelastic models and approaches. With minimum residuals, the FE analysis using the Burgers model is better in representing the viscoelastic behaviour of the polymer. Using the proposed method, the creep and recovery behaviour of the polymer in the linear viscoelastic regime can be estimated with acceptable accuracies. This will avoid the long-term creep and recovery experimentations, which are often costly. Finally, this work will open the door to further research on new polymeric materials or blends in the linear and nonlinear viscoelastic regimes.

Appendices

Appendix A

Creep of a Burgers viscoelastic model for an applied constant stress:

Boundary load information,

$$ \sigma = \sigma _{o},\quad \dot{\sigma } = \ddot{\sigma } =0. $$

The constitutive equation of the Burgers model is first derived using the Laplace transform and its inverse on Eqs. (4) and (5) (Kelly 2013). The strain solution follows a similar approach; however, it utilises the initial conditions given above.

Laplace transforms of the individual strains in Eq. (5) are

$$\begin{gathered} \mathcal{L} \left ( \varepsilon _{1} \right ) = \overline{\varepsilon }_{1} \left ( s \right ) = \frac{1}{E_{1}} \overline{\sigma } \left ( s \right ) \\ \mathcal{L} \left ( \varepsilon _{2} \right ) = \overline{\varepsilon }_{2} \left ( s \right ) = \frac{1}{\eta _{1}} \frac{\overline{\sigma } \left ( s \right )}{s} \\ \mathcal{L} \left ( \varepsilon _{3} \right ) = \overline{\varepsilon }_{3} \left ( s \right ) = \frac{1}{\eta _{2}} \frac{\overline{\sigma } \left ( s \right )}{\left ( s+ {E_{2}} / {\eta _{2}} \right )} \end{gathered}$$

Replacing individual Laplace transformed strains in Eq. (4) in Laplace form,

$$ \overline{\varepsilon } = \overline{\varepsilon }_{1} \left ( s \right ) + \overline{\varepsilon }_{2} \left ( s \right ) + \overline{\varepsilon }_{3} \left ( s \right ) $$

$$ \overline{\varepsilon } \left ( s \right ) = \frac{\overline{\sigma } \left ( s \right )}{E_{1}} + \frac{\overline{\sigma } \left ( s \right )}{s \eta _{1}} + \frac{\overline{\sigma } \left ( s \right )}{\eta _{2} \left ( s+ {E_{2}} / {\eta _{2}} \right )} $$

The total strain and stress are related in Laplace form. Rearranging,

$$ \overline{\varepsilon } (s) \left [ E_{1} \eta _{1} \eta _{2}^{2} s^{2} + E_{1} \eta _{1} E_{2} \eta _{2} s \right ] = \overline{\sigma } \left ( s \right ) [ \eta _{1} \eta _{2}^{2} s^{2} + \left ( \eta _{1} E_{2} \eta _{2} + E_{1} \eta _{2}^{2} + E_{1} \eta _{1} \eta _{2} \right ) s+ E_{1} E_{2} \eta _{2} ] $$

Taking the Laplace inverse of the above expression,

$$ \mathcal{L}^{-1} \left ( \overline{\varepsilon } (s) \left [ E_{1} \eta _{1} \eta _{2}^{2} s^{2} + E_{1} \eta _{1} E_{2} \eta _{2} s \right ] \right ) = E_{1} \eta _{1} \eta _{2}^{2} \ddot{\varepsilon } + E_{1} \eta _{1} E_{2} \eta _{2} \dot{\varepsilon } $$

$$\begin{aligned} &\mathcal{L}^{-1} \left ( \overline{\sigma } \left ( s \right ) \left [ \eta _{1} \eta _{2}^{2} s^{2} + \left ( \eta _{1} E_{2} \eta _{2} + E_{1} \eta _{2}^{2} + E_{1} \eta _{1} \eta _{2}^{2} \right ) s+ E_{1} E_{2} \eta _{2} \right ] \right ) \\ &\quad = \eta _{1} \eta _{2}^{2} \ddot{\sigma } + \left ( \eta _{1} E_{2} \eta _{2} + E_{1} \eta _{2}^{2} + E_{1} \eta _{1} \eta _{2} \right ) \dot{\sigma } + E_{1} E_{2} \eta _{2} \sigma \end{aligned}$$

Putting the above inverse in standard form,

$$ \frac{\eta _{1} \eta _{2}}{E_{2}} \ddot{\varepsilon } + \eta _{1} \dot{\varepsilon } = \frac{\eta _{1} \eta _{2}}{E_{1} E_{2}} \ddot{\sigma } + \left ( \frac{\eta _{1}}{E_{1}} + \frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{2}} \right ) \dot{\sigma } +\sigma $$

In order to solve the creep, the initial conditions and the boundary stress are considered, see given data, we consider the instantaneous response,

$$ \varepsilon _{o} = \frac{\sigma _{0}}{E_{1}} $$

The initial strain rate can be found by integrating the constitutive equation at the instant of load application. That is,

$$ \sigma =0,\qquad t= 0^{-};\qquad \sigma = \sigma _{o},\qquad t= 0^{+} $$

Integrating between the time range,

$$ \int _{0^{-}}^{0^{+}} \frac{\eta _{1} \eta _{2}}{E_{2}} \ddot{\varepsilon } + \int _{0^{-}}^{0^{+}} \eta _{1} \dot{\varepsilon } = \int _{0^{-}}^{0^{+}} \frac{\eta _{1} \eta _{2}}{E_{1} E_{2}} \ddot{\sigma } + \int _{0^{-}}^{0^{+}} \left ( \frac{\eta _{1}}{E_{1}} + \frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{2}} \right ) \dot{\sigma } + \int _{0^{-}}^{0^{+}} \sigma dt $$

$$ \frac{\eta _{1} \eta _{2}}{E_{2}} \dot{\varepsilon }_{0} + \eta _{1} \varepsilon _{0} = \frac{\eta _{1} \eta _{2}}{E_{1} E_{2}} \dot{\sigma }_{0} + \left ( \frac{\eta _{1}}{E_{1}} + \frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{2}} \right ) \sigma _{0} + \left . \sigma _{0} t \right \vert _{0^{-}}^{0^{+}} $$

The initial strain rate becomes

$$ \dot{\varepsilon }_{0} = \frac{\sigma _{o}}{\eta _{1} \eta _{2}} \left ( \eta _{1} + \eta _{2} \right ) $$

The Laplace transform of the constitutive form given in Eq. (6) is

$$ \mathcal{L} \left ( \eta _{1} \dot{\varepsilon } + \frac{\eta _{1} \eta _{2}}{E_{2}} \ddot{\varepsilon } \right ) =\mathcal{L} \left ( \sigma _{0} \right ) $$

$$ \frac{\eta _{1} \eta _{2}}{E_{2}} \left ( s^{2} \overline{\varepsilon } \left ( s \right ) -s \varepsilon _{0} - \dot{\varepsilon }_{0} \right ) + \eta _{1} \left ( s \overline{\varepsilon } \left ( s \right ) - \varepsilon _{0} \right ) = \frac{\sigma _{0}}{s} $$

Replacing the initial conditions and the boundary stress,

$$ \frac{\eta _{1} \eta _{2}}{E_{2}} \left ( s^{2} \overline{\varepsilon } \left ( s \right ) -s \frac{\sigma _{0}}{E_{1}} - \frac{\sigma _{o}}{\eta _{1} \eta _{2}} \left ( \eta _{1} + \eta _{2} \right ) \right ) + \eta _{1} \left ( s \overline{\varepsilon } \left ( s \right ) - \frac{\sigma _{0}}{E_{1}} \right ) = \frac{\sigma _{0}}{s} $$

Collecting similar terms,

$$ \frac{\eta _{1} \eta _{2}}{E_{2}} s^{2} \overline{\varepsilon } \left ( s \right ) +s \eta _{1} \overline{\varepsilon } \left ( s \right ) = \sigma _{0} \left [ \frac{1}{s} +s \frac{\eta _{1} \eta _{2}}{E_{2} E_{1}} + \frac{\eta _{1} + \eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}} \right ] $$

Factoring out and simplifying,

$$ \frac{\eta _{1} \eta _{2}}{E_{2}} \overline{\varepsilon } \left ( s \right ) s(s+ \frac{E_{2}}{\eta _{2}} )= \sigma _{0} \left [ \frac{1}{s} +s \frac{\eta _{1} \eta _{2}}{E_{2} E_{1}} + \frac{\eta _{1} + \eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}} \right ] $$

$$ \frac{\eta _{1} \eta _{2}}{E_{2}} \overline{\varepsilon } \left ( s \right ) = \sigma _{0} \left [ \frac{1}{s^{2} (s+ \frac{E_{2}}{\eta _{2}} )} + \frac{\eta _{1} \eta _{2}}{E_{2} E_{1}} \frac{1}{(s+ \frac{E_{2}}{\eta _{2}} )} + \left ( \frac{\eta _{1} E_{1} + \eta _{2} E_{1} + \eta _{1} E_{2}}{E_{2} E_{1}} \right ) \frac{1}{s(s+ \frac{E_{2}}{\eta _{2}} )} \right ] $$

$$ \frac{\eta _{1} \eta _{2}}{E_{2}} \overline{\varepsilon } \left ( s \right ) = \sigma _{0} \left [ \frac{1}{s^{2} (s+ \frac{E_{2}}{\eta _{2}} )} + \frac{\eta _{1} \eta _{2}}{E_{2} E_{1}} \frac{1}{(s+ \frac{E_{2}}{\eta _{2}} )} + \left ( \frac{\eta _{1} \eta _{2} E_{1} + \eta _{2}^{2} E_{1} + \eta _{1} \eta _{2} E_{2}}{E_{2}^{2} E_{1}} \right ) \frac{\frac{E_{2}}{\eta _{2}}}{s(s+ \frac{E_{2}}{\eta _{2}} )} \right ] $$

Taking the inverse Laplace,

$$\begin{aligned} \frac{\eta _{1} \eta _{2}}{E_{2}} \varepsilon &= \sigma _{0} \left [ \frac{\eta _{2}}{E_{2}} t- \left ( \frac{\eta _{2}}{E_{2}} \right )^{2} \left ( 1- e^{- {E_{2}} / {\eta _{2}} t} \right ) + \frac{\eta _{1} \eta _{2}}{E_{2} E_{1}} e^{- {E_{2}} / {\eta _{2}} t}\right. \\ &\left.\quad {}+ \frac{\eta _{1} \eta _{2} E_{1} + \eta _{2}^{2} E_{1} + \eta _{1} \eta _{2} E_{2}}{E_{2}^{2} E_{1}} \left ( 1- e^{- {E_{2}} / {\eta _{2}} t} \right ) \right ] \end{aligned}$$

$$ \varepsilon = \sigma _{0} \left [ \frac{1}{\eta _{1}} t- \frac{\eta _{2}}{\eta _{1} E_{2}} \left ( 1- e^{- {E_{2}} / {\eta _{2}} t} \right ) + \frac{1}{E_{1}} e^{- {E_{2}} / {\eta _{2}} t} + \left ( \frac{1}{E_{2}} + \frac{\eta _{2}}{\eta _{1} E_{2}} + \frac{1}{E_{1}} \right ) \left ( 1- e^{- {E_{2}} / {\eta _{2}} t} \right ) \right ] $$

Simplifying the creep becomes

$$ \varepsilon = \sigma _{0} \left [ \frac{1}{E_{1}} + \frac{1}{\eta _{1}} t+ \frac{1}{E_{2}} \left ( 1- e^{- {E_{2}} / {\eta _{2}} t} \right ) \right ] $$

Creep of four element generalised Maxwell model under constant stress:

Initial and boundary load information,

$$ \sigma = \sigma _{o},\qquad \dot{\sigma } = \ddot{\sigma } =0. $$

The constitutive of the generalised Maxwell model given in Eq. (9) is customised for a four-element model. That is,

$$ \left [ \left ( E_{1} + \eta _{1} \frac{\partial }{\partial t} \right ) \left ( E_{2} + \eta _{2} \frac{\partial }{\partial t} \right ) \right ] \sigma = \left [ \eta _{1} E_{1} \frac{\partial }{\partial t} \left ( E_{2} + \eta _{2} \frac{\partial }{\partial t} \right ) + \eta _{2} E_{2} \frac{\partial }{\partial t} \left ( E_{1} + \eta _{1} \frac{\partial }{\partial t} \right ) \right ] \varepsilon $$

$$\begin{aligned} &\left [ \left ( E_{1} E_{2} + \left ( E_{1} \eta _{2} + E_{2} \eta _{1} \right ) \frac{\partial }{\partial t} + \eta _{1} \eta _{2} \frac{\partial ^{2}}{\partial t^{2}} \right ) \right ] \sigma \\ &\quad = \left [ \eta _{1} E_{1} E_{2} \frac{\partial }{\partial t} + E_{1} \eta _{1} \eta _{2} \frac{\partial ^{2}}{\partial t^{2}} + \eta _{2} E_{2} E_{1} \frac{\partial }{\partial t} + \eta _{1} \eta _{2} E_{2} \frac{\partial ^{2}}{\partial t^{2}} \right ] \varepsilon \end{aligned}$$

$$ E_{1} E_{2} \sigma + \left ( E_{1} \eta _{2} + E_{2} \eta _{1} \right ) \dot{\sigma } + \eta _{1} \eta _{2} \ddot{\sigma } = \left ( \eta _{1} E_{1} E_{2} + \eta _{2} E_{2} E_{1} \right ) \dot{\varepsilon } + \left ( E_{1} \eta _{1} \eta _{2} + \eta _{1} \eta _{2} E_{2} \right ) \ddot{\varepsilon } $$

Putting in standard form,

$$ \sigma + \left ( \frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}} \right ) \dot{\sigma } + \frac{\eta _{1} \eta _{2}}{E_{1} E_{2}} \ddot{\sigma } = \left ( \eta _{1} + \eta _{2} \right ) \dot{\varepsilon } + \eta _{1} \eta _{2} \left ( \frac{1}{E_{1}} + \frac{1}{E_{2}} \right ) \ddot{\varepsilon } $$

In order to solve the creep, the initial conditions and the boundary stress are considered, see given data, we consider the instantaneous response,

$$ \varepsilon _{o} = \frac{\sigma _{o}}{\left ( E_{1} + E_{2} \right )} $$

The initial strain rate can be found by integrating the constitutive equation at the instant of load application. That is,

$$ \sigma =0, t= 0^{-}; \sigma = \sigma _{o}, t= 0^{+} $$

Integrating between the time range,

$$\begin{aligned} &\int _{0^{-}}^{0^{+}} \eta _{1} \eta _{2} \left ( \frac{1}{E_{1}} + \frac{1}{E_{2}} \right ) \ddot{\varepsilon } + \int _{0^{-}}^{0^{+}} \left ( \eta _{1} + \eta _{2} \right ) \dot{\varepsilon } \\ &\quad = \int _{0^{-}}^{0^{+}} \frac{\eta _{1} \eta _{2}}{E_{1} E_{2}} \ddot{\sigma } + \int _{0^{-}}^{0^{+}} \left ( \frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}} \right ) \dot{\sigma } + \int _{0^{-}}^{0^{+}} \sigma dt \end{aligned}$$

$$ \eta _{1} \eta _{2} \left ( \frac{1}{E_{1}} + \frac{1}{E_{2}} \right ) \dot{\varepsilon }_{o} + \left ( \eta _{1} + \eta _{2} \right ) \varepsilon _{o} = \frac{\eta _{1} \eta _{2}}{E_{1} E_{2}} \dot{\sigma }_{o} + \left ( \frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}} \right ) \sigma _{o} + \left . \sigma _{0} t \right \vert _{0^{-}}^{0^{+}} $$

The initial strain rate becomes

$$ \dot{\varepsilon }_{o} = \sigma _{o} \left [ \left ( \frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}} \right ) - \frac{\eta _{1} + \eta _{2}}{E_{1} + E_{2}} \right ] \left ( \frac{E_{1} E_{2}}{\eta _{1} \eta _{2} \left ( E_{1} + E_{2} \right )} \right ) $$

The Laplace transform of the constitutive form given in Eq. (11) is

$$ \mathcal{L} \left ( \left ( \eta _{1} + \eta _{2} \right ) \dot{\varepsilon } + \eta _{1} \eta _{2} \left ( \frac{1}{E_{1}} + \frac{1}{E_{2}} \right ) \ddot{\varepsilon } \right ) = \mathcal{L} \left ( \sigma _{0} \right ) $$

$$ \left ( \eta _{1} + \eta _{2} \right ) \left ( s \overline{\varepsilon } \left ( s \right ) - \varepsilon _{0} \right ) + \frac{\eta _{1} \eta _{2} \left ( E_{1} + E_{2} \right )}{E_{1} E_{2}} \left ( s^{2} \overline{\varepsilon } \left ( s \right ) -s \varepsilon _{0} - \dot{\varepsilon }_{0} \right ) = \frac{\sigma _{0}}{s} $$

Replacing the initial conditions,

$$\begin{aligned} &\left (\eta _{1} + \eta _{2} \right ) \left( s \overline{\varepsilon } \left ( s \right ) - \frac{\sigma _{0}}{E_{1} + E_{2}} \right ) + \frac{\eta _{1} \eta _{2} \left ( E_{1} + E_{2} \right )}{E_{1} E_{2}} \left ( s^{2} \overline{\varepsilon } \left ( s \right ) -s \frac{\sigma _{0}}{E_{1} + E_{2}} - \sigma _{o} \left[ \left ( \frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}} \right )\right.\right. \\ &\left.\left.\quad {}- \frac{\left ( \eta _{1} + \eta _{2} \right )}{\left ( E_{1} + E_{2} \right )} \right ] \left ( \frac{E_{1} E_{2}}{\eta _{1} \eta _{2} \left ( E_{1} + E_{2} \right )} \right ) \right ) = \frac{\sigma _{0}}{s} \end{aligned}$$

Collecting similar terms,

$$\begin{aligned} &s \overline{\varepsilon } \left ( s \right ) \left [ \left ( \eta _{1} + \eta _{2} \right ) +s \frac{\eta _{1} \eta _{2} \left ( E_{1} + E_{2} \right )}{E_{1} E_{2}} \right ] \\ &\quad = \sigma _{0} \left [ \frac{1}{s} +s \frac{\eta _{1} \eta _{2}}{E_{1} E_{2}} + \frac{\eta _{1} + \eta _{2}}{E_{1} + E_{2}} + \left ( \frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}} \right ) - \frac{\eta _{1} + \eta _{2}}{E_{1} + E_{2}} \right ] \end{aligned}$$

$$ s \overline{\varepsilon } \left ( s \right ) \left [ \left ( \eta _{1} + \eta _{2} \right ) +s \frac{\eta _{1} \eta _{2} \left ( E_{1} + E_{2} \right )}{E_{1} E_{2}} \right ] = \sigma _{0} \left [ \frac{1}{s} +s \frac{\eta _{1} \eta _{2}}{E_{1} E_{2}} + \left ( \frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}} \right ) \right ] $$

Solving for the strain in Laplace form,

$$ \overline{\varepsilon } \left ( s \right ) = \frac{\sigma _{0}}{\eta _{1} \eta _{2} \left ( E_{1} + E_{2} \right )} \left [ \frac{E_{1} E_{2} + s^{2} \eta _{1} \eta _{2} +s \left ( \eta _{1} E_{2} + \eta _{2} E_{1} \right )}{s^{2} \left ( s+ \frac{\left ( \eta _{1} + \eta _{2} \right ) E_{1} E_{2}}{\eta _{1} \eta _{2} \left ( E_{1} + E_{2} \right )} \right )} \right ] $$

The strain in Laplace form can be arranged for computing the inverses individually.

Let

$$ a= \frac{\left ( \eta _{1} + \eta _{2} \right ) E_{1} E_{2}}{\eta _{1} \eta _{2} \left ( E_{1} + E_{2} \right )} $$

Hence the strain in Laplace form can be expressed as

$$\begin{aligned} \overline{\varepsilon } \left ( s \right ) &= \sigma _{0} \left [ \frac{1}{\eta _{1} \eta _{2} \left ( E_{1} + E_{2} \right )} \left [ \left ( E_{1} E_{2} \frac{1}{s^{2} \left ( s+a \right )} \right ) + \left ( \eta _{1} \eta _{2} \frac{1}{s+a} \right ) \right ]\right. \\ &\left.\quad {}+ \frac{\frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}}}{\eta _{1} + \eta _{2}} \left ( \frac{a}{s \left ( s+a \right )} \right ) \right ] \end{aligned}$$

Taking the inverse of Laplace on the above expression,

$$\begin{aligned} \mathcal{L}^{-1} \left ( \overline{\varepsilon } \left ( s \right ) \right ) &= \mathcal{L}^{-1} \left( \sigma _{0} \left[ \frac{1}{\eta _{1} \eta _{2} \left( E_{1} + E_{2} \right)} \left[ \left( E_{1} E_{2} \frac{1}{s^{2} \left ( s+a \right)} \right) + \left ( \eta _{1} \eta _{2} \frac{1}{s+a} \right ) \right] \right. \right. \\ &\left. \left. \quad {} + \frac{\frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}}}{\eta _{1} + \eta _{2}} \left ( \frac{a}{s \left ( s+a \right )} \right ) \right ] \right ) \end{aligned}$$

The creep becomes

$$ \varepsilon \left ( t \right ) = \sigma _{0} \left [ \frac{t}{\eta _{1} + \eta _{2}} + \frac{e^{-at}}{E_{1} + E_{2}} + \frac{1}{\eta _{1} + \eta _{2}} \left ( \frac{\eta _{2}}{E_{2}} + \frac{\eta _{1}}{E_{1}} - \frac{1}{a} \right ) \left ( 1- e^{-at} \right ) \right ] $$

Appendix B

Matlab Script of curve fitting for a Burgers viscoelastic model (color online):

Matlab Script of curve fitting for a four-element generalised Maxell viscoelastic model (color online):