Time-dependent properties of high-density polyethylene with wood/graphene nanoplatelets reinforcement

The effect of graphene nanoplatelets (GNPs) on the long-term performance of wood fiber/high-density polyethylene (HDPE) composite is investigated by using short-term creep tests with an efficient, faster data analysis approach. Previously, it was shown that the addition of GNPs at 15 wt% into HDPE reduces the visco-plastic (VP) strain developed during 2 h creep by (cid:1) 50%. The current study shows that 25 and 40 wt% wood content in HDPE reduce the VP strains developed during 2 h creep time by >75% with no noticeable effect of the increased wood content. However, further addition of GNPs results in more than 90% total reduction in the VP strains. The current study shows that the development of the VP strains in the hybrid composites follows Zapas model. Viscoelastic (VE) response of these composites is nonlinear and thus is described by Schapery's model. Parameters for VP and VE models are obtained from the creep experiments and were validated in a separate loading-unloading test sequence. Results show a very good agreement between experiments and predictions for the studied materials as long as the micro-damage is not present.


| INTRODUCTION
Wood polymer composites (WPCs) have been the material of choice in many applications ranging from construction to common consumer goods.They possess desirable properties such as being lightweight, offering ease of formability, and for the relative low cost of the constituents.Conversely, the hydrophilic nature of wood and the low durability of WPCs disqualify them for longterm use, especially in load-bearing applications.Both wood fibers (WF) and polymers are viscoelastic (VE) materials, [1,2] thus WPCs also possess VE characteristics in their properties with a large dependence on time, temperature, humidity, and load. [2]Under constant load, the material tends to creep and deform, and depending on the conditions, this deformation can become permanent leading to progressive damage and eventually to a catastrophic failure when reaching the critical strain after sufficiently long time. [2]igh-density polyethylene (HDPE) is one of the most widely used polymers in general and is one of the polyolefins used for the manufacturing of WPCs.It is a recyclable polymer known for its relatively good mechanical performance and chemical inertness together with the high hydrophobicity that makes it attractive for wide range of applications.However, similar to other polymers, HDPE undergoes creep deformation when subjected to constant load for elongated times.Several studies have been performed to characterize and model this behavior at different conditions due to its importance in the respective applications as neat polymer [3][4][5] or as a matrix in WPC. [6,7]10][11][12] Research has been conducted to study and characterize this behavior together with investigating ways to reduce or delay damage. [13]Increasing the filler/fiber content is found to have an inverse relation to the creep strains. [14]The larger the concentration of the filler in the matrix, the lower the development of creep strains and the longer the service life.The addition of modifiers to improve the adhesion between fibers and matrix is another way through which the creep resistance has been improved. [15]Better compatibility between the matrix and the reinforcement leads to a more efficient interface and improved load transfer reducing the time-dependent deformation. [16]Many factors play an important role in the creep response of WPC such as the type of polymer, the degree of crystallinity, the type of reinforcement and the compatibility of the constituents, and so forth. [13][19][20] Modification with carbonaceous additives at the nanoscale (such as graphene or carbon nanotubes) to achieve multifunctional materials with added thermal or electrical conductivities is relatively new field of research that is growing constantly.These functionalities require rather large amounts of the nanoparticles to form the conductive network in the thermoplastic matrix [17] which results in a great change of the polymer microstructure.They have also been reported to have a significant effect on the time-dependent performance of a variety of polymers [21][22][23][24] and advanced polymer composites of carbon and glass fibers. [25,26]ue to the impracticality to obtain experimental data for long-term performance, a common practice is to perform a short-term test and describe the captured material behavior with models to predict the behavior at different test conditions.Most of the available models for VE are simplistic in nature and are typically used to describe the performance of the material in the linear VE region only (i.e., under stresses below 50% of the ultimate tensile stress for the majority of materials that exhibit linear and nonlinear response [13] ).In these models, the constitutive response depends on time and is expressed based on the Boltzmann superposition principle in its integral or differential form.In order to extend the understanding of performance beyond this region and cover the materials that do not show linear VE response (such as HDPE [27] ), more generalized models have been developed.These models are adopted based on the principles of thermodynamics and the expansion of Gibb's free-energy theory for viscoelasticity relating the internal state variables and using the linear evolution laws with respect to thermodynamic forces.In these, functions are added describing the nonlinear material behavior with respect to applied stress, temperature, and so forth, and are experimentally identified.Models proposed by Findley [28] and Schapery [29] and the following modifications [30,31] are examples of this type of analysis approach where single creep compliance cannot describe the material behavior under different stresses.
In the current research, results of experiments performed to analyze the VE, and VP behavior of WPC modified with graphene nanoplatelets (GNPs) at different concentrations are presented.The efficiency and effects of the single type of nano-or micro-scale reinforcement as well as the synergistic effects of the combination of these two are investigated facilitated by the study and comparison of different concentrations of these reinforcements.An analysis approach presented in [29,32] and used in [12,24,33] has been used to separate the VE and VP responses and to allow studying the effect of the reinforcement on these behaviors individually.To the best of the authors' knowledge, no studies have been performed to characterize the extent to which these additives can affect the time-dependent performance of WPCs, or to separate their effect on the VP and VE behavior.The effect of interaction between the WF and the nanoparticles to alter the time-dependent response is still unclear, especially at high reinforcing content.Validation of the employed approach is done by comparison of loadingunloading test results with predications using parameters extracted from short-term creep tests.The results generated in this type of analysis are of extreme importance as they offer a deeper understanding of the behavior of an emerging class of material with potential for use in thermal management of buildings.Such information can also enrich the database proposed in [34] to correct the criteria of material selection based on long-term performance rather than static properties.

| Materials
High-density polyethylene (HDPE, MG9647S) is supplied by Borealis AG in a pellet form having a density of 0.964 g/cm 3 .A masterbatch of graphene in HDPE (heXo HDPE1-V20/35) is purchased from NanoXPLOre having 35 wt% GNPs.According to the manufacturer's technical data sheet, [35] the graphene platelets in the masterbatch have an average of 40 layers of graphene sheets with a thickness of 20 nm and a lateral dimension of 50 μm.The graphene platelets are functionalized at the edges (by the supplier) to enhance compatibility with the polymer.WF from Scandinavian Wood Fiber AB is a sawdust of spruce and pinewood with 75% of its particles in the size range 200-400 μm (based on the manufacturer's description and the confirmed measurement in supplementary materials in Ref. [17]).Maleic anhydride-grafted HDPE (MAPE), E265, provided by DuPont, is used as a compatibilizer.The masterbatch was used as received without any modifications or further characterizations to the nanoplatelets.

| Methods
The composites are manufactured by melt blending process in a co-rotating twin-screw extruder followed by compression molding.The extruder was run with zone temperatures of 180-200 C and a rotation speed of 300 rpm.To achieve a better dispersion, WPC was extruded for another round with 120 rpm.Details of the manufacturing parameters can be found elsewhere [17] and are presented in a sketch in the supplementary information document (Figure S1).Examination of the quality and distribution of the reinforcement is presented in the supplementary information through scanning electron microscopy (SEM) images of test-fracture surfaces of the samples (Figure S2 in the supplementary information document).The formulations, constituents, and sample notations of the composites used for this study are given in Table 1.The composite samples are named adapting the following notations "XWPCY" where X refers to the content of WF in the composites (in wt%) and Y represents the fraction of the nanoparticles (when present) in the composite (also in wt%), while WPC refers to wood polymer composite.It is worth noting that these composites have been developed for advanced applications [17,36] where thermal and electrical conductivities were target properties besides the mechanical performance.This justifies the selection of the high content of the reinforcement and the relevance of the studied properties to the targeted applications.For comparison, properties of the neat polymer and other nanocomposites analyzed in previous work are also presented when needed (full set of results are available in Ref. [24]).Specimens with dimensions 200 mm in length, 15 mm in width, and 4 mm in nominal thickness were cut by waterjet followed by thorough drying for 8 h at 80 C before testing.
Instron 3366 universal tester equipped with 10 kN load cell and pneumatic grips was utilized for the test in a displacement control mode on 100 mm gauge length of the specimen.Strains were measured using a clip-on extensometer (Instron 2620-601) having a base of 50 mm placed at the center of the specimen.A detailed description of the different tests employed in this study is provided in the following sections.
1. Short-term creep tests were performed on a single specimen of each material formulation.One creep test consists of multiple loading ramps at a rate of 3.5 mm/min (based on extensometer) to specific stress for increasing time intervals interrupted by strain recovery steps following the procedure described in [21,24] as shown in Figure 1A.Specimens were loaded to the predefined stress, held at that stress for a time t 1 where creep strain is observed, then unloaded to almost zero load level (4 N) to recover the reversible strains and held at that load for a time 8t 1 .The stresses were selected to cover around 25%-75% of the maximum stress in the tensile test as can be seen in values of Table 2.The load was not set to 0 N during the unloading due to the technical limitation of the machine.However, 4 N equals to a very small stress value that is considered ineffective in preventing the full recovery of reversible strains.Creep loading intervals (Δt 1 , Δt 2 , Δt 3 , Δt 4 ) of 10, 20, 30, and 60 min were used in one test amounting to a total of 2 h.These intervals were used to obtain VP strain dependence on time.The dependence on stress is obtained by loading the specimen to higher stresses for the same time intervals.A minimum of three stress levels are needed on each material formulation to obtain the VP parameters.Any unrecovered strains at the end of each recovery step are attributed to permanent deformation (VP strains) and were subtracted from the subsequent strains for the analysis of viscoelasticity.In fact, the presence of damage in form of fiber debond or matrix cracking could also contribute to the permanent deformations and increase the permanent irreversible strain (IR-strain).Whether the damage is interfering with VP strains or not is not fully established at this point and will be discussed further in the article.In the following, when only VP strains are discussed independently from damage, they are referred to as such (VP strains), otherwise, they are referred to as IR-strains. 24 Cyclic loading-unloading-recovery ramps (further referred to as cyclic loading test) were applied to investigate the development of damage and its effect on the material response.The test was designed to load the specimen in multiple successive cycles with incremental strains (see Figure 1B).Each cycle consists of loading the specimen to a specific strain with a constant crosshead speed of 2 mm/min, unloading it to 5 N, and waiting for the specimen to recover the reversible strains.The time of recovery was Â4 the time required for the loading and unloading ramps.
After each straining cycle, the stiffness of the material was measured by introducing an additional step where the specimen is loaded to 0.2% strain, unloaded, and left to recover.The modulus was then calculated from the slope of the unloading part of the stress-strain curve in the low stress region between 0.05% and 0.15% strain.Any unrecovered strains at the end of the recovery time were considered irreversible.Damage development was assessed by normalizing the stiffness after each maximum loading to the initial stiffness of the material, that is, E(σ max )/E 0 .

| Cyclic loading test
Figure 2A shows the time-dependent response of the WPC at different loading steps.The response appears in a form of hysteresis loops where the unloading curves follow different paths than those of the loading but approaches the original point after sufficient recovery time has passed.Any unrecoverable strains after that time are removed during data redaction and the offset stress-strains curves are plotted for each subsequent step.Upon the next loading step, the path and slope of the curves are identical to the previous steps at the low-stress regions.At regions of higher stresses, the curves slightly shift downwards indicating degradation in stiffness due to the introduction of damage.An enlarged part of the curves showing the initial slopes corresponding to the stiffness is presented in the supplementary information document (Figure S3). Figure 2B shows the hysteresis loops of the studied material systems as they are being loaded to around The small peaks correspond to the 0.2% strain step for the modulus calculation, and the larger peaks are the loading ramps, which were stepwise increased from 0.4% to 6%.
20 MPa.With the addition and increase of the fraction of reinforcements, the curves are shifting upwards with a decrease in the size (area) of the loop and a reduction in the deformation strains.This is an indication of the increased stiffness and the increased ability of the material to store energy.In other words, the addition of the stiff particles increased the fraction of the material that contributes to the elastic response compared to that of viscous response (typically from the amorphous part of the polymer) that tends to dissipate the energy as it flows at the same applied stress.Utilizing the change in slope of hysteresis loops, finding the modulus after each applied stress/strain would give an idea about the state of damage in the material subjected to a certain load level.This is depicted in Figure 2C,D where the ratio between modulus at damaged state (increased stress or strain) to the initial modulus is plotted against applied stress and strain.These data will be used in the evaluation of VP strain limits in the upcoming discussion.Observation of F I G U R E 3 IR-strains as a result of the maximum applied stress in the cyclic loading test.
the normalized moduli with respect to the applied stress/ strain shows that all materials share a general trend.No change or slight steady decrease of modulus up to a certain level of the load was observed, followed by a significant stiffness drop until failure occurs.A decrease in stiffness of about 8% (for 25WPC) and 13% (for 40WPC) is observed at failure.With the increase in the amount of the reinforcement the material gets stiffer but also more brittle which leads the material to fail at lower strains.For example, the sample 40WPC7.5 is heavily reinforced and sustains very high loads but, at the same time, the probability for damage initiation increases significantly especially at the polymer-reinforcement interfaces.Moreover, the damage may develop in a more severe way than other materials until a critical defect size is reached that causes ultimate failure.For all the formulations studied, the noticeable stiffness degradation indicated by this test occurs at stresses higher than 15-17.5 MPa which suggests that below these stress levels, all IR-strains could be attributed to the viscoplastic deformation.Beyond these stress values, a possible interaction between damage and VP strains could occur.The development of IR-strains as a result of the max applied stress is shown in Figure 3, where the materials show a trend of reduced IR-strains with increased amounts of reinforcement and that the lowest IR-strains occur for sample 40WPC7.5.Once again, the materials show a rapid increase in the permanent deformation around the 20 MPa applied stress due to the expected interaction between the damage and plastic deformation of the polymer chains.As the applied stress increases, WF contributes more efficiently to reducing these IRstrains (in sample 40WPC) compared to the higher or equivalent efficiency of the nanoparticles at lower stress levels (sample 25WPC10).Even though these two materials have different loading histories in the experiments, the trend has been validated by simulations as well (see the supplementary information document -Figure S4).

| Short-term creep tests
Representative creep strain responses for 60 min loading at 10 and 15 MPa creep stresses are shown in Figure 4A,B, respectively for the different studied formulations.It can be seen that all components of the strain response are affected by the addition of the reinforcement.The instantaneous strain is much lower for the 40WPC7.5 which is the result of higher stiffness of the sample.It was reduced from around 0.5% for the 25WPC to about 0.2% for 40WPC7.5-morethan 50% reduction.The maximum creep strain has been significantly reduced at the increase of the reinforcements, going from 1% to 0.4% for the same two mentioned formulations.The rate of development of the strain is also different going into a slower rate at both primary as well as secondary creep stages where the material with higher content of reinforcements reaches the semi-constant rate earlier with shorter transition region compared to other materials (see the inlet to the graphs presenting the strain rate analysis for the loading part of the creep curves).A more pronounced reduction is noticed at the secondary stage of creep where 40WPC7.5 shows the slowest development of creep strains and the largest recovered strains in the recovery part.
It is worth noting that samples 25WPC10 and 40WPC show similar responses at 10 MPa creep test, but they slightly differ at the 15 MPa.On the other hand, the addition of the extra 7.5 wt% of GNP to the 40WPC results in a total strain reduction of $62% (compared to sample 25WPC).This observation is similar to what has been introduced in the previous section.This implies that for small loads, increasing the amount of microreinforcement might have the same effect as the addition of the nanoparticles.However, there is a practical and technical limit for the maximum content of fibers (typically 20%-30% for short natural fibers [37] ) that can be added to the polymer while still being wetted and processed properly.After this value, mechanical properties show significant degradation and further improvement of the properties might only be possible by the addition of the nano-reinforcement.
For the same loading history, the addition of nanoparticles to WPC results in a significant decrease in the IR-strains compared to the WPC.The accumulated IRstrains in 25WPC at the end of the test were 0.57% (after loading to 7.5 MPa ! 10 MPa) which was reduced to 0.28% (a reduction of 50.3%) with the addition of 10 wt% GNPs.Similarly, the IR-strain was reduced by 60% at the addition of the 7.5 wt% GNP to the 40WPC sample (reduced from 1.58% to 0.63%) after loading to (10 MPa ! 15 MPa !17.5 MPa !20 MPa).Compared to previously tested neat HDPE, [24] the reduction is even higher both in creep strains and IR-strains due to the synergistic effect of the two types of reinforcements surpassing the effect of either of them present individually as can be seen in Figure 4C,D.
The development of IR-strains at the end of 2 h creep test as a function of the applied creep stress is presented in Figure 5. Compared to the IR-strains from the cyclic loading test, the materials show the same trend with much higher values at comparable stresses due to the effect of different loading histories.
The creep rupture of materials occurs at the last stage of creep.It can be seen from Figure 6A that only samples with 25 wt% WF (i.e., 25WPC and 25WPC10) have demonstrated that type of failure while the other formulations exhibited sudden brittle failure.Possibly, the advanced state of damage reached for highly reinforced composite could be one reason why failure occurs catastrophically with a sudden release of energy in contrast to the typical rapid progressive failure in creep.Moreover, it can be inferred that the small amount of polymer between the reinforcements is highly stressed and constrained which prevents it from deforming.This will result in suppressing matrix plasticity and promoting other modes of failure such as cracking and debonding at multiple sites that could eventually connect forming one larger crack of a critical size causing the material fracture.This is also evident by the reduced intensity of the tertiary creep for sample 25WPC10 compared to 25WPC F I G U R E 5 IR-strains as a result of the applied creep stress.
(circled in the figure).A complementary analysis of the development of strain rate is presented in the supplementary information document (Figure S5).A previous study [38] has correlated the development of the tertiary stage of creep in WPC to the level of stress in the creep test where low stress levels (15% of the ultimate tensile stress) do not allow for development of this creep stage while medium and high stresses (45%-60% of the ultimate tensile stress) causes the development of tertiary creep and a brittle fracture.In the current study, all materials have failed at high stress levels (>60% of tensile stress) which might suggest that the advanced stage of damage mentioned above is more responsible for the observed trends.
The improvement in the time-dependent response of the studied composites is caused by the large constrains offered by the rigid particles to the movement of the polymer chains.Homogeneously distributed particles modify the properties of the matrix leading to an overall improved response.Clusters and agglomerates have a negative impact on the properties preventing them from reaching the theoretically calculated reinforcing potential.Those have not been spotted often when the fracture surface was examined but could lead to even higher improvements if they were completely avoided.For thermoplastic materials produced by conventional extrusion, dispersion of the nanoparticles is challenging due to the high viscosity of the melt and the need for strong shear forces to properly distribute them, which simultaneously can damage the reinforcement or degrade the length of the polymer chains.However, incorporating the nanoparticles from a masterbatch in which they are well dispersed initially, leads to improved quality of composites and easier hybridization with other types of reinforcements.Moreover, the addition of nanoparticles helps increase the chances to further modify the composite when the technical limits of micro-sized reinforcement have been reached.It has been reported [39,40] that the nanoparticles contribute to the improved interface between the fiber and the matrix which subsequently results in an improvement in the macro-properties of the composites.[43][44] Different mechanisms are responsible for this improvement some of which are mechanical, others are chemical, or driven by polarity. [45]he nanoparticles used in this study are functionalized at the edges which promote their improved bonding to the functional hydroxyl groups at the surface of fibers or those of the compatibilizer, especially since the reinforcement percentage is quite high and large amounts of particle end up in the vicinity of the fiber surface.
Figure 6B depicts IR-strain as a function of creep strain which reveals a power function dependency as was previously reported by Starkova et al. [46] on polyolefins nanocomposites.This work presented a general formula in the form of a power function relating the creep strains (ϵ cr ) and the IR-strains (ϵ IR ) by ϵ IR = C ϵ cr n .A wide range of materials tested in creep for different creep/recovery times followed the same dependency with n ranging between 0.7 and 2. In the current study, the materials show the same dependency where the value of n lies at the higher end for all formulations with no clear correlation to the amount or type of reinforcement.Obviously, with the increase in the amount of reinforcement, the curves shift to the left and downwards, while they progress to the right with the increase in the creep stress and time.This information, however, concerns the overall IR-strains independent of whether it comes from viscoplasticity or damage which is favorable for the quick assessment of the structural stability of the material.On the other hand, stress rather than strain has a direct influencing effect on the amount of developed IR-strain.Therefore, in the next sections, a more detailed analysis of the strains is presented and modeled giving a more comprehensive understanding of the material's time-dependent response.

| MATERIAL MODELING
It is assumed that the strain ε caused by stress ramp σ t ð Þ can be decomposed into VE strain (ε VE ) and VP strain (ε VP ) as: The linear form of Equation ( 1) may be considered as the first terms in the series expansion of a more complex relationship where the interaction of all possible mechanisms is expected.Strains ε VE and ε VP in Equation ( 1) are nonlinear with respect to stress.The damage state which evolves with load and time is characterized by a damage function d.The simplest assumption is that the damage development is an elastic process and hence d does not depend explicitly on time.Thus, damage function d depends only on the highest (most damaging) stress state experienced during the stress ramp.In the uniaxial loading case, which has been considered in the current study, it is the highest experienced axial tensile stress, σ max , and hence d can be found in Ref. [47] expressed as: where E 0 is the initial modulus at the undamaged state and E σ max ð Þ is the modulus of the damaged material after application of certain stress σ max For the current work, this parameter is defined to be 1 (damage is not present) based on the results from the cyclic loading test presented in Figure 2 (showing that damage starts at stresses higher than those used for the creep test).Therefore, the IR-strains at the end of each recovery period are attributed to only ε VP .Zapas [48] showed that VP strain development could be described by the function presented below: where C VP , m and M are fitting parameters, t Ã and σ Ã are arbitrarily chosen normalization constants.For this study, t Ã is selected to be equal to the total creep time (120 min) at one stress level, and σ Ã is selected to be 1 MPa to simplify calculations.Using the single specimen approach described in Ref. [33], VP strain has been calculated.Due to the previous loading history, calculations become slightly more difficult starting from the second loading step.The detailed derivation of this methodology can be found in Ref. [33].This approach is valid only if the scatter between the specimens is too small (as is the case in this study shown by the small standard deviations of the tensile test results in Table 2) so that the single specimen could be assumed representative for the material in general.However, validating with more than one specimen is advisable whenever possible or representativeness of the selected specimen should otherwise be confirmed. [49]The equations used for parameter identification are presented belowat the first stress level (σ 0 ), VP strain has been fitted using Equation ( 4) and at k-th stress level (σ k ¼ n k Á σ 0 ) the VP strain is fitted using Equation ( 5) In Figure 7A, experimental VP strains (represented by symbols) and modeled VP strain using Zapas model (represented by solid lines) for specimen 25WPC at different creep stresses are presented.All studied materials showed the same trends.The curves for the VP strain development for all materials can be found in the supplementary information document (Figure S6).It can be seen that the IR-strains are increasing with creep time and with increasing applied creep stress.The predictions are very accurate for the first creep stresses but at higher stresses deviates the model significantly underestimates the development of VP strains.The cyclic loading test indicated that the damage develops at stress levels higher than 17.5 MPa, thus the IR-strains in creep tests at the end of each recovery period were assumed as pure VP strain.However, the high deviation of experimental IR-strains from the VP model at high stresses could possibly be caused by the development of damage and its interference with VP strains.It is inferred that the different loading scenarios in the two types of the tests might result in different damage states especially for materials with high content of reinforcement.Previous investigations on hybrid reinforced materials with lower reinforcement content (up to 11%) did not show such deviations. [21]This suggests the need to find more reliable means to identify and evaluate the damage, where damage can be separated from VP strain to study the interaction between damage and the development of VP strains in more detail.A new approach is currently examined to advance the modeling procedure followed here and investigate the time effect on the development of VP strains.An example of the improved model by the new approach can be found in the supplementary information document (Figure S7) and will be discussed in more detail in an upcoming publication.This preliminary approach, however, did not apply to all materials in this study and needs further investigations.Nevertheless, it was possible to obtain VP parameters for the stress levels up to 15 MPa.The   obtained VP parameters are presented in Table 3. Data for neat polymer from the previous study are also presented for comparison.After VP parameters have been obtained, they can be used to remove VP strains from the total creep strains and obtain pure VE strains to perform VE analysis.The methodology for VP strain removal presented in Ref. [50] has been used.For the first loading steps, the VP is removed using Equation ( 6) and for the steps with previous loading history Equation (7) has been used.
In Equation ( 6), the superscript c refers to the current (in this case the first step), and in Equation ( 7) k refers to any k th step preceded by step k À 1.The result is the pure VE strains which can be separately analyzed.Figure 8A shows an example of the resulting pure VE strains for sample 25WPC at different stress levels together with the obtained creep compliances in Figure 8B.Each set of lines of the same type represent the different loading times within the same stress level.
As can be seen from Figure 8B, the creep compliance curves show nonlinear VE response with respect to stress even at relatively low-stress levels.Linear VE region (where compliance curves from the different stresses coincide with each other) could not be observed for this material as well as all other materials in this study.At higher stress levels, the individual lines of the compliance curves within the same stress level deviate from each other, thus indicating improper removal of VP-strain.This could be a result of damage that is developing and not taken into account in the model as was mentioned above.The nonlinear VE compliance indicates that analysis and modeling of viscoelasticity should be done using nonlinear VE models, for example, Schapery's. [30]Schapery's model for nonlinear VE was modified to account for VP strains and damage.In the onedimensional case, the material model is presented as: In Equation ( 8) ε el represents the initial elastic strain in undamaged material which, generally speaking, might be a nonlinear function of stress.The "Reduced time" ψ is introduced as ψ ¼ R t 0 dt a σ , and consequently ψ ¼ R τ 0 dτ a σ .Parameters g 1 , g 2 and the shift factor a σ are stress invariant-dependent as well as temperature-, humidity-, aging-, and so forth, dependent functions.In the current study, the environmental factors were kept constant, thus these functions are only stress-dependent.The transient part of the VE response is characterized by ΔS(ψ), which according to the theory by Schapery [30] does not depend on stress and has a form of Prony series, where C i are constants and τ i are retardation times.It is possible to apply the VE material model in Equation (8) to the creep strain and recovery strains separately by specifying the time of each ramp where t [0, t 1 ] is the time for creep and t > t 1 is the recovery time.
The expressions used to model the creep and recovery strains are then as in Equations ( 10) and (11), respectively.The detailed methodology for VE parameter and function identification has been presented in Ref. [51], Application of the expressions presented above results in the fitting curves presented in Figure 9.The model curves fit almost perfectly to the experimental data.The small deviations were observed at the higher stress levels as can be seen for sample 25WPC during loading and unloading presented in Figure 9A,B.Curves for all materials showed similar trends and are presented in the supplementary information document (Figure S10).Comparison for the loading and unloading response of all materials at 10 MPa stress level is presented in Figure 9C,D, respectively.
Figure 10 shows the visual presentation of the obtained parameters g 1 , g 2 , and a σ as a function of stress for sample 25WPC.The detailed fitting functions for these parameters for the studied materials are presented in Table 4 and the Proney Series coefficients in Table 5.
In theory, the parameters obtained from the analysis of viscoplasticity (Table 3), viscoelasticity (Tables 4 and  5), and damage obtained from the cyclic loading test can be used in Equation (1) to simulate an entire response of the material in this study at any loading profile below 15 MPa (where damage is not observed).In order to validate the obtained parameters and models, the curves from the cyclic loading test were Since these tests were performed in a displacement control mode (strain is an input), the material model in Equation ( 1) needs to be re-formulated so that it is a function of strain rather than stress.This formulation is not trivial and it has been done in Ref. [51], where the model has been written in an inverted-incremental form.Derived equations used for the simulations are: where, where the parameter functions, g 1 , g 2 , C i , and ψ have been introduced earlier and the expression D was introduced in [52] relating the VE parameters.
It was observed that finding a single region for all materials to obtain the elastic modulus is impossible, especially for the 25WPC, which exhibited highly non-linear behavior.This could be due to nonlinear viscoelasticity, which was observed even at very low stresses (as discussed in earlier paragraphs).The region for elastic modulus determination had a high effect on the simulation, as it is demonstrated in Figure 11A.It is demonstrated on two simulation curves by using modulus obtained from a region close to the origin (near 0% strain) and modulus obtained from the strain interval of 0.05%-0.15%strain.It can be observed that the use of region close to 0% strain shows a perfect fit with experimental data, while the use of region at higher strains leads to deviation from experimental data.These nuances will be further investigated in future papers.
For consistency reasons, a common region for elasticity was chosen between 0.05% and 0.15% of strain.Figure 11B shows the simulated curves of the stress-strain loops for all materials where the stress is around 15 MPa and damage is not present.The other loops are shown in the supplementary information document (in Figure S12), there again deviation is seen at higher stresses where damage is expected to initiate.It was also observed that in simulations outside the experimentally performed and studied creep tests, the deviations become more notable.This is an indication that the extrapolation of parameters outside the experimentally obtained values could lead to significant simulation inaccuracies.It is worth noting that all the issues occurring with applying the current approach on  the highly reinforced materials are the center for a future publication dedicated to advancing the modeling procedure and are not within the scope of this particular study.

|
This paper presented the effect on the long-term behavior of reinforcing highly nonlinear viscoelastic (VE) polymer (HDPE) with two types of reinforcement at different scales.It was observed that the nanoparticles contributed to the improvement of creep behavior of composites in general, but especially at low-stress levels.At higher stress levels, both nanosized particles (GNPs) and micron-sized fibers (wood fibers) contribute to the improved material properties.Time-dependent properties in terms of creep strain and accumulated Viscoplastic (VP) strains have been significantly reduced at all amounts of reinforcements.It was shown that mechanical response under creep loading increasingly improved with the increased amount of both types of reinforcement studied.The incorporation of large amounts of the two reinforcements not only reduces the irreversible strains (IR-strains) but also changes the mode of failure in the composites.Samples with less wood and more nanoparticles exhibited higher resistance to creep with a less catastrophic mode of failure at higher stresses.This is due to the participation of nanoparticles in mitigating the damage initiation at the fiber/ matrix interface by sharing larger loads in the matrix.These particles significantly restrict the movement of the polymer chains and thus reduce the extensive deformations in the matrix leading to improved time-dependent response.Depending on the application, exposure to low frequent loads for an extended time might require different compositions of the reinforcement compared to exposure to higher loads in a less frequent manner.It was also demonstrated that it is possible to separate and analyze the different modes of material behavior.Viscoelasticity, viscoplasticity, and damage were identified and characterized using appropriate models.Both materials -WPC and WPC modified with GNPs are shown to follow Zapas model in the development of VP strains.All types of studied composites showed highly nonlinear VE response even at low stresses as was proven by the creep compliance curves.Using simple short-term creep tests performed on a single specimen, it was possible to extract parameters to describe VP and VE.The parameters obtained from these analyses have been successfully used to simulate the experimental data of the cyclic loading test.The materials with high amounts of reinforcement might impose the need for advancing the modeling approach to cover the specificity of the materials and the interaction between damage and VP strains.
To date, it is not clear how this interaction happens, and no models are available to accurately describe it.A more accurate method for the elastic regions for highly nonlinear materials is required to avoid challenges in modeling their behavior.

F I G U R E 2
The hysteresis loops showing the time-dependent response of: (A) sample 25WPC at different strain levels; (B) the different materials around the same stress level (20 MPa); (C, D) the evolution of stiffness of materials with the applied max stress and strain, respectively.

F
I G U R E 4The strain response as a function of time for 60 min creep cycle at constant creep stress of 10 MPa (A) and 15 MPa (B) with inlets showing the strain rate analysis for example materials.Creep maximum strains (C) and accumulated IR-strains (D) with % change at the addition of reinforcement compared to the neat polymer at 10 and 15 MPa stress levels.Note that in (A) and (B) only the first 10,000 s of the test are shown to demonstrate the trends and full recovery of the reversible strain was possible at the end of the set recovery time.Curves are also shifted to start from origin independent of the accumulated VP strain from the previous step.

F
I G U R E 6 (A) failure of the materials at the different ultimate creep stresses.The circles highlight the rapid increase in strain prior to failure.Note that the sample 25WPC10 failed outside the extensometer gauge length and the curve is obtained from the displacement data instead.Curves for the strain analysis at failure are presented in the supplementary information; (B) IR-strain dependence on the maximum strain at creep test of 30 and 60 min compared to the trend in the work by Starkova et al. on PP/MWCNT.

F I G U R E 7
Development of IR-strains with time at different stress levels for 25WPC (A) and for all materials at 10 MPa creep stress (B).The R 2 values next to the fitting lines refer to the correlations between the experiment and model (coefficient of determinations).

F
I G U R E 8 Pure VE strains (elastic and time-dependent parts) after subtraction of VP strains (A), and creep compliance curves (B) for sample 25WPC.Note that the curves in (A) are presented for the first 10,000 s of the test step to demonstrate the trends.Full graph showing the complete recovery of the strains is presented in the supplementary information document (FigureS8) demonstrating complete recovery of the strains at the end of the test and the creep compliance curves of the other materials (FigureS9).

T A B L E 3
Viscoplastic parameters from Equation(3) I G U R E 9 Experimental and modeled VE strains for creep (A, C) and recovery (B, D) strains for 25WPC (A, B) and all materials at 10 MPa creep stress (C, D).

F I G U R E 1 0
Parameters of the nonlinear VE functions.Functions of the other materials are presented in the supplementary information document (Figure S11) T A B L E 4 Functions of the VE parameters used for fitting VE strains.

F I G U R E 1 1
Simulations of the experimental points below 15 MPa using the parameters obtained from the individual of the VE and VP response for 25WPC (A), and for all other materials (B) Tensile properties and the ranges of stresses selected for the creep test.
T A B L E 2 T A B L E 5 Coefficients of Prony series