Modeling of Creep Behavior of Particulate Composites with Focus on Interfacial Adhesion Effect

Evaluation of creep compliance of particulate composites using empirical models always provides parameters depending on initial stress and material composition. The effort spent to connect model parameters with physical properties has not resulted in success yet. Further, during the creep, delamination between matrix and filler may occur depending on time and initial stress, reducing an interface adhesion and load transfer to filler particles. In this paper, the creep compliance curves of glass beads reinforced poly(butylene terephthalate) composites were fitted with Burgers and Findley models providing different sets of time-dependent model parameters for each initial stress. Despite the finding that the Findley model performs well in a primary creep, the Burgers model is more suitable if secondary creep comes into play; they allow only for a qualitative prediction of creep behavior because the interface adhesion and its time dependency is an implicit, hidden parameter. As Young’s modulus is a parameter of these models (and the majority of other creep models), it was selected to be introduced as a filler content-dependent parameter with the help of the cube in cube elementary volume approach of Paul. The analysis led to the time-dependent creep compliance that depends only on the time-dependent creep of the matrix and the normalized particle distance (or the filler volume content), and it allowed accounting for the adhesion effect. Comparison with the experimental data confirmed that the elementary volume-based creep compliance function can be used to predict the realistic creep behavior of particulate composites.


Introduction
The majority of lightweight design parts are realized using reinforced composites [1]. Knowledge of their reliable and relevant physico-mechanical materials data, such as stiffness, strength, or creep compliance, is indispensable for engineering and design purposes. Due to the viscoelastic nature of polymer matrices, the dimensional stability in structural applications, and the corresponding long-term behavior, are often estimated by either extrapolation of short-term data or under accelerated testing conditions using higher initial stresses, elevated temperatures, and relative humidity [2]. However, these investigations are time and energy-consuming, especially in the case of creep experiments. To overcome this, modeling as an analytical approach is currently used to predict the performance of the composites.
The creep behavior of composites depends on the creep response of a matrix and a dispersed phase (fibers or particles), the content of the dispersed phase, interface adhesion between filler and matrix, and experimental parameters such as temperature and initial stress [3]. During the creep experiment, the interface adhesion in the composites can be reduced (by, e.g., delamination), leading to accelerated creep. The creep compliances increase faster with time than in the case of pure viscoelastic creep. If the load transfer from matrix to filler is decreased, the creep behavior of composites is determined dominantly by the matrix creep.
The filler/matrix adhesion depends on the compatibility between the filler and the matrix. Bek et al. [4] compared the long-term creep compliances of neat and recycled polypropylene (PP) reinforced with untreated wood fibers. They found a superior creep resistance for neat PP composites and attributed it to a higher crystallinity and better interfacial adhesion due to trans-crystalline structures around the wood fibers occurring only for neat PP. Chandekar et al. [5] studied the effect of various chemical pre-treatments of jute fibers on the creep recovery behavior of PP-jute composites. The creep evaluation using the Burgers model revealed that creep strains decreased due to the chosen fiber pre-treatment in the following manner: untreated > potassium permanganate > alkali (NaOH) > silane > NaOH with maleic anhydride grafted PP. Similar results were presented by Marcovich and Villar [6] for the creep behavior of low-density polyethylene compounded with untreated wood flour and modified with both organic peroxide and MAH. They concluded that MAH increased the interfacial adhesion between filler and matrix. Furthermore, they found that the functionalization using organic peroxide hindered crystallization and led to composites with lower crystallinity and, thus, more pronounced creep behavior.
In general, the creep at constant load is divided into three ranges of deformation [7]. After an instantaneous elastic deformation due to applied initial stress, primary creep, characterized by a decreasing creep rate, occurs. After a certain time, the creep rate becomes almost constant, indicating the range of secondary creep. Finally, tertiary creep occurs with a fast-increasing creep rate leading to a fracture. To quantitatively describe the measured creep behavior of polymers and composites, several models were developed, Table 1. They can be divided into two categories [8]: Table 1. Some models of time-and stress-dependent creep compliances for polymers.

Models Creep Strain/Creep Compliance
Burgers with time t, initial stress σ 0 , stiffnesses of Burgers model E 1 , E 2 and viscosities of Burgers model η 1 , η 2

Modified Burgers
with relaxation time τ and fitting parameters a and b Findley power law with initial strain ε 0 , transient strain ε + , exponent n, reference time t 1 Findley modified power law with reference stress σ 1

Bailey-Norton
with coefficient A and exponents m and n + Power series with relaxation times τ i with relaxation times τ i

Rheological Models
They combine springs and dashpots in a suitable manner to reproduce measured behavior (Maxwell model, Kelvin/Voigt model, Burgers model, etc.). The model parameters are the stiffnesses of springs and viscosities of dashpots, which can partly be determined by experiments or fitting procedures.

Empirical Equations
They approximate the measured behavior and the model parameter are determined by fitting procedures.
The Burgers model is an in-series arrangement of the Maxwell model (elastic spring E 1 in series with a viscous dashpot η 1 ) and the Kelvin model (elastic spring E 2 parallel to viscous dashpot η 2 ) [9]. For the constant initial stress σ 0 , it reproduces qualitatively the creep behavior with the initial elastic strain and primary and secondary creep strains; Equation (1), Table 1. For long times, only the dashpot having viscosity η 2 determines the creep behavior in a time-linear manner, which does not display the real creep behavior of most materials. To give the Burgers model more generality and to account for both constant and tertiary creep strains, Sarabi [10] introduced a time-dependent viscosity η 2 , Equation (2).
Both the Findley power law model [7], Equation (3), and the Findley modified power law model [11], Equation (4), are empirical power functions describing creep with an elastic strain and a transient, non-elastic strain. Equation (3) was employed to describe creep for a variety of materials, such as polymers [12][13][14], metals [7], and concrete [15]. However, some materials show a non-linear viscoelastic behavior in cases of higher initial stresses. Thus, Equation (4) provided successful estimates of long-term creep on the basis of shortterm creep data at higher stress levels [11].
The models of Burgers and Findley are widely used to depict the creep behavior of polymers and polymer composites. In Table 1, further models are listed (Equations (5)-(7)) using exponential and power functions or power series. These models may provide suitable and better adjustments to measured creep data, as can be found in [16][17][18].
All models of Table 1 allow for fitting adaption of measured creep data. As a result, one can tabulate model parameters with respect to initial stresses and temperatures. An interface adhesion of polymer composites is a hidden parameter within other parameters. It means that the effect of the reduced and time-dependent interface adhesion on the creep behavior of polymer composites has not been considered yet. The aim of this study is to derive a time-dependent compliance function containing the adhesion term on the base of a cube in cube model to predict the creep behavior of reinforced polymers.

Measured Creep Compliances and Their Fitting Using Burgers and Findley Models
The measured creep compliance curves exhibit an initial elastic compliance followed by primary and partly secondary creep behavior within the measuring time of 1000 h; Figure 1. As expected, the creep compliances increase with time and initial stress and decrease with filler volume content. In the secondary creep range, the creep rates decrease as glass beads are a factor 20 stiffer than the poly(butylene terephthalate) (PBT) matrix, thus, exhibiting negligible creep for small initial stresses [19]. Only PBT GB30 at σ 0 = 22 MPa clearly reached the range of the secondary creep with the constant creep rate.
Burgers model and Findley power law model were used to fit the creep curves; Figure 1. They provide good fits up to 250 h, especially as long as the initial stresses remain small. The determined fitting parameters are summarized in Tables 2 and 3. In general, the Findley power law model performs better in the investigated time ranging from 250 to 1000 h, respectively. From the creep curves, it is obvious that the secondary creep has not been reached yet, Figure 1a,b,d-f. Except for PBT GB30 at 40% σ max , the Burgers model performed better as it shows nice linear creep, which can be attributed to a constant slippage rate of polymer chains in load direction [20], Figure 1c.
Fitting the measured creep compliances using the Burgers model shows that the model parameters E 1 , E 2 , η 1, and η 2 decrease with initial stresses (=softening effect due to more free volume) and increase with glass beads volume contents (=stiffening effect due to hard fillers), Table 2. In general, the values of the model parameters decrease if the upper time limit is increased from 250 to 1000 h. This behavior is an effect of fitting procedure as viscous flow behavior mainly addressed by a dashpot η 1 comes more into play for longer time ranges with a consequence of decreasing viscosities η 1 .
The parameter E 1 can be associated with Young's modulus of the material, and the E 1 values in Tables 3 and 4 are similar to Young's moduli of PBT determined by tensile tests [21]. The first data point of the creep experiment is not determined at the beginning but with a delay of 10 to 20 s. This means that it contains already small creep portions with the consequence that the start value of E 1 can be chosen up to 5% smaller as E 2 may partly compensate for that in the fitting process.
Fitting the measured creep compliances using the Burgers model shows that the model parameters E1, E2, η1, and η2 decrease with initial stresses (=softening effect due to more free volume) and increase with glass beads volume contents (=stiffening effect due to hard fillers), Table 2. In general, the values of the model parameters decrease if the upper time limit is increased from 250 to 1000 h. This behavior is an effect of fitting procedure as viscous flow behavior mainly addressed by a dashpot η1 comes more into play for longer time ranges with a consequence of decreasing viscosities η1.
The parameter E1 can be associated with Young's modulus of the material, and the E1 values in Tables 3 and 4 are similar to Young's moduli of PBT determined by tensile tests [21]. The first data point of the creep experiment is not determined at the beginning but with a delay of 10 to 20 s. This means that it contains already small creep portions with the consequence that the start value of E1 can be chosen up to 5 % smaller as E2 may partly compensate for that in the fitting process. There have been some trials to connect the model parameters with physical properties, e.g., E 2 was associated with the stiffness of amorphous polymer chains [22], which plays a role in the short-time creep stage. However, the moduli E 2 of Table 2 are far too large for amorphous polymers. Similarly, viscosities η 1 and η 2 were connected to damage of crystalline or oriented non-crystalline regions and characteristic times of damaging processes [23]. Nevertheless, as damage mechanisms in crystalline regions differ from those of amorphous regions with respect to physical processes and characteristic times, it is clear that all model parameters except E 1 are only fit parameters.
Fitting of the measured creep compliance curves using the Findley power law model provides values for the model parameters, Table 3, for the upper fitting time limits of 250 h and 1000 h. The parameter ε 0 represents the initial elastic strain, and as expected, it increases with the initial stress and decreases with the content of glass beads. The conversion to the initial modulus E 1 shows coincidence to the moduli determined for the Burgers model and Young's moduli from the literature within the accuracy of measurements.There have been some trials to connect the model parameters with physical properties, e.g., E 2 was associated with the stiffness of amorphous polymer chains [22], which plays a role in the short-time creep stage. However, the moduli E 2 of Table 2 are far too large for amorphous polymers. Similarly, viscosities η 1 and η 2 were connected to damage of crystalline or oriented non-crystalline regions and characteristic times of damaging processes [23]. Nevertheless, as damage mechanisms in crystalline regions differ from those of amorphous regions with respect to physical processes and characteristic times, it is clear that all model parameters except E 1 are only fit parameters.
Fitting of the measured creep compliance curves using the Findley power law model provides values for the model parameters, Table 3, for the upper fitting time limits of 250 h and 1000 h. The parameter ε 0 represents the initial elastic strain, and as expected, it increases with the initial stress and decreases with the content of glass beads. The conversion to the initial modulus E 1 shows coincidence to the moduli determined for the Burgers model and Young's moduli from the literature within the accuracy of measurements. There have been some trials to connect the model parameters with physical properties, e.g., E2 was associated with the stiffness of amorphous polymer chains [22], which plays a role in the short-time creep stage. However, the moduli E2 of Table 2 are far too large for amorphous polymers. Similarly, viscosities η1 and η2 were connected to damage of crystalline or oriented non-crystalline regions and characteristic times of damaging processes [23]. Nevertheless, as damage mechanisms in crystalline regions differ from those of amorphous regions with respect to physical processes and characteristic times, it is clear that all model parameters except E1 are only fit parameters.
Fitting of the measured creep compliance curves using the Findley power law model provides values for the model parameters, Table 3, for the upper fitting time limits of 250 h and 1000 h. The parameter ε0 represents the initial elastic strain, and as expected, it increases with the initial stress and decreases with the content of glass beads. The conversion to the initial modulus E1 shows coincidence to the moduli determined for the Burgers model and Young's moduli from the literature within the accuracy of measurements. The power coefficient n ranges from 0.35 to 0.55, meaning that the Findley model can describe only primary creep. As the n increases with the initial stresses, it provides enhanced creep. The introduction of glass beads leads to the smaller n indicating the reduction of the creep. Interestingly, n remains unchanged if the upper fitting time limit is increased from 250 h to 1000 h. This makes the Findley model a good tool for practitioners to predict creep behavior from short-term creep experiments [12,24].
In this evaluation, the determined transient strain ε + corresponds to the creep strain after 1 h, which should increase with the initial stress. This is only true for the neat PBT, whereas it decreases for the glass-beads filled PBT. Furthermore, the transient strains ε + of the glass bead-filled PBT exceed those of the neat PBT by a factor 2 to 3. This may mean that the glass bead-filled PBT exhibits an enhanced strain rate during the first hour of the creep before the creep equilibrium is established. The power coefficient n ranges from 0.35 to 0.55, meaning that the Findley model can describe only primary creep. As the n increases with the initial stresses, it provides enhanced creep. The introduction of glass beads leads to the smaller n indicating the reduction of the creep. Interestingly, n remains unchanged if the upper fitting time limit is increased from 250 h to 1000 h. This makes the Findley model a good tool for practitioners to predict creep behavior from short-term creep experiments [12,24].
In this evaluation, the determined transient strain ε + corresponds to the creep strain after 1 h, which should increase with the initial stress. This is only true for the neat PBT, whereas it decreases for the glass-beads filled PBT. Furthermore, the transient strains ε + of the glass bead-filled PBT exceed those of the neat PBT by a factor 2 to 3. This may mean that the glass bead-filled PBT exhibits an enhanced strain rate during the first hour of the creep before the creep equilibrium is established.
The evaluation of the creep compliance curves using empirical models always provides parameter sets depending on the initial stresses and materials compositions (e.g., the content of glass beads). It means that for calculation, one has to measure the creep compliances first with subsequent parameter evaluation before one is able to predict the creep behavior beyond the measured time or for other filler contents.
In the models of Burgers and Findley, the parameter E 1 represents the elastic initial response on the external load. Therefore, Young's modulus E 1 can be used as an input parameter that takes into account the filler content already before introducing it to the creep models. This can be achieved by the use of cube in cube models with various complexity, which were used for more than 60 years to calculate Young's moduli of particulate composites as a function of matrix modulus E M , filler modulus E F and filler volume content v F [25].
Young's modulus E 1 is the input parameter to most creep models. Therefore, it seems reasonable to use models to calculate it before introducing it to the creep models. To do so, cube in cube models with various complexity were used for more than 60 years to calculate Young's moduli of the particulate composites as the function of matrix modulus E M , filler modulus E F and filler volume content v F [25].

Determination of the Time Dependent Compliance Function of Glass Beads Filled Composites
The simplest cube in cube models was proposed by Paul [26] and Ishai-Cohen [27] in order to calculate Young's moduli of particle-filled composites. They assume a cubic elementary volume (EV) with edge length D + a containing a single inclusion of diameter D shown for glass beads in Figure 2a. For modeling purposes, the inclusion has to be transformed into a cube. Then the EV is separated into a matrix part and a composite part Due to the transformation, LC becomes (k D) and LM becomes (a + 1 − k) with the efficiency factor k. It considers that less than the maximum cross-section contributes to the stress transfer and thus depends on the particle shapes. For spheres, it is determined by the condition It is obvious that the strains of EV, matrix part, and composite part are different. The strain of EV εEV represents the macroscopic strain measured in a creep experiment, whereas the strains of matrix part εM and composite part εC represent "unknown" microscopic strains yet, which are related to each other for perfect interface adhesion according to [28] , with the normalized inclusion distance d to which a filler volume content vF is linked by 1 Thus, the microscopic strains are linked to the macroscopic strain by the dimensional ratio of the EV and the inclusion or by the filler volume content. However, if the cross-sectional ratio of filler inclusion and matrix exceeds 1 (identical cross-sections of filler particle and matrix), the related filler volume content vF exceeds 26%. This yields increasing shear stresses [28]; hence, the cross-sectional ratio represents a restricting/limiting parameter in this approach. In the creep experiment, the test bars are loaded with a constant force F0. This allows the expression of the strains in terms of initial stress σ0 and moduli of EV EEV, matrix part EM, and composite part EC. strain of EV (11) strain of matrix part (12) strain of composite part (13) Introduction of Equations (11)- (13) to Equation (9) yields During creep experiments, there may occur delamination between the matrix and filler depending on time and initial stress. This reduces the interface adhesion and the  [26] with the external load F 0 .
Due to the transformation, L C becomes (k D) and L M becomes (a + 1 − k) with the efficiency factor k. It considers that less than the maximum cross-section contributes to the stress transfer and thus depends on the particle shapes. For spheres, it is determined by the condition It is obvious that the strains of EV, matrix part, and composite part are different. The strain of EV ε EV represents the macroscopic strain measured in a creep experiment, whereas the strains of matrix part ε M and composite part ε C represent "unknown" microscopic strains yet, which are related to each other for perfect interface adhesion according to [28] with the normalized inclusion distance d to which a filler volume content v F is linked by Thus, the microscopic strains are linked to the macroscopic strain by the dimensional ratio of the EV and the inclusion or by the filler volume content. However, if the crosssectional ratio of filler inclusion and matrix exceeds 1 (identical cross-sections of filler particle and matrix), the related filler volume content v F exceeds 26%. This yields increasing shear stresses [28]; hence, the cross-sectional ratio represents a restricting/limiting parameter in this approach. In the creep experiment, the test bars are loaded with a constant force F 0 . This allows the expression of the strains in terms of initial stress σ 0 and moduli of EV E EV , matrix part E M , and composite part E C .
Introduction of Equations (11)- (13) to Equation (9) yields During creep experiments, there may occur delamination between the matrix and filler depending on time and initial stress. This reduces the interface adhesion and the load transfer to the filler particles. The structure of Equation (14) allows for the introduction of an adhesion factor k adh , or an adhesion function k adh (t) in case of a time dependency. The square of the adhesion function can be geometrically interpreted as the effectively available cross-section for stress transfer [28]. Thus, it is possible to adjust it between k adh = 0 (no adhesion) and k adh = 1 (perfect adhesion). After dividing Equation (14) by the initial stress σ 0 , one ends up with the time-dependent creep compliance of the EV The creep of the composites part can be considered elastic because the external load acts mainly on the filler as it is distributed between filler and matrix by the ratio E F :E M . This means that the matrix of the composites part experiences much less stress than the matrix part. Therefore, its viscoelastic behavior can be neglected here. Equation (15) follows that the time dependency of the composite creep J EV (t,σ 0 ) is completely determined by the matrix creep J M (t,σ 0 ) for the same given initial stress σ 0 and that for not too high initial stresses, the composite part contributes to the creep of particlefilled composite only elastically representing the time-independent constant term. Furthermore, quantitative information about the time-dependent interface adhesion is accessible via creep experiments as all parameters of Equation (15) are known except for k adh .
If the measured creep compliance curves of neat PBT are introduced in Equation (15), creep compliance curves of particle-filled composites can be calculated for any given adhesion factor. The calculated creep compliance curves for PBT GB20 and PBT GB30 are shown in Figure 3 as the lines for the adhesion factors 0, 0.2, 0.4, 0.6, and 1.
The adhesion factor k adh = 1 leads to the lowest creep compliance curves (solid lines), with values being smaller than measured creep compliance curves, Figure 3. With the lower adhesion factors, the calculated creep compliances achieve coincidence with the measured ones. Furthermore, the slope of the measured creep compliance curves exceeds the slope of the calculated ones. This means that the adhesion factor decreases with time. A pointwise evaluation provides the time-dependent adhesion function, Figure 3. Over 250 h, the adhesion factors k adh decrease from 0.55 to 0.36 (PBT GB20 at σ 0 = 17 MPa), 0.50 to 0.32 (PBT GB20 at σ 0 = 22 MPa), 0.53 to 0.45 (PBT GB30 at σ 0 = 11 MPa), and 0.57 to 0.34 (PBT GB30 at σ 0 = 22 MPa) and it decreases faster at the higher initial stresses. Between 250 and 1000 h, the adhesion factors tend towards an asymptotic final adhesion factor except for PBT GB30 at σ 0 = 22 MPa, which decreases linearly due to the secondary creep. The adhesion factors were found to decrease logarithmically with time. Therefore, the timedependent adhesion factors were fitted using the function k adh (t) = −∆k adh * ln(t) + k t=1 h adh (16) with the incremental slope ∆k adh and initial adhesion factor k t=1 h adh . The fit parameters are depicted in Table 4. For the fitting time limit of 250 h, both ∆k adh and k t=1 h adh increase with the initial stress. For 1000 h the k t=1 h adh also increases with initial stress. However, for the ∆k adh the decrease is found for PBT GB20 and the increase for PBT GB30, which can be attributed to the effects of the asymptotic behavior for long times. For PBT GB30 at 22 MPa, the segmental fitting procedure due to the occurrence of the secondary creep provides better fitting of the measured creep compliance.  The introduction of Equation (16) to Equation (15)   The introduction of Equation (16) to Equation (15)

Creep Tests
Creep tests were performed according to DIN EN ISO 899-1 using a Zwick 1411 creep stand (ZwickRoell, Ulm, Germany -software testXpert) equipped with optical strain measurements at 23 • C and 50% r.h. [29]. The initial stresses σ 0 were chosen to be 11, 17, and 22 MPa, respectively, for neat PBT. This corresponds roughly to yield strength ratios of 20, 30, and 40%. Initial stresses σ 0 were chosen to be 17 and 22 MPa for PBT GB 20 and 11 and 22 MPa for PBT 30 GB. Creep strains were measured for maximum of 1000 h. Subsequently, the creep strain curves were corrected for errors in the measured initial length due to non-vertical mounting and sliding in clamps using a method proposed in [30].

Creep Modeling and Determination of Adhesion Factors
Measured creep strain curves were converted to creep compliance curves and, in the first step, fitted with Burgers model, Equation (1), and Findley power law model, Equation (3), using an Excel solver in the least square approximations for 250 and 1000 h. In the second step, the model parameters were adjusted manually using creep compliance curves in a time logarithmic scale to compensate for the domination of the short-term range. The extended cube in cube model, Equation (15), was used to calculate the creep compliances J EV (t, σ 0 , v F ) of the glass bead-filled PBT composites in a pointwise manner to determine time-dependent interfacial adhesion factors. If the matrix creep compliance has been determined experimentally, all quantities in Equation (15) are known except for k adh (t). Thus, Equation (16) is used to determine the adhesion functions k adh (t) by fitting the time-dependent adhesion factors using Excel solver and then introduced to Equation (15).

Conclusions
It is known that the creep behavior of composites depends on filler/fiber volume content, initial stress as well as delamination between filler and matrix occurring at longlasting or high loads. The creep behavior was investigated on poly(butylene terephthalate) (PBT), and its glass beads reinforced composites. First, the measured creep compliance curves were evaluated for 250 and 1000 h using the empirical models of Burgers and Findley. It turned out that the Findley model performs better as long as there is only the primary creep. If secondary creep occurs, the Burgers model is more suitable. The values of the model parameters show that increasing filler volume content leads to stiffening and increasing initial stress to softening. For predicting purposes, the problem is that one gets different parameter sets for each initial stress, and changes in interface adhesion are hidden within them. Both models contain Young's modulus E 1 as a parameter. Therefore, Paul's cube in cube model, being an elementary volume approach (EV), was used as it offers the chance to introduce E 1 as a filler volume-dependent quantity. Further, the analysis showed that the separation of the EV into the matrix and composite parts leads to time-dependent compliance, which only depends on the creep compliance of the matrix and the mean distance between filler particles or filler volume content. In the composite part, the majority of the stress is transferred to the stiff filler particles. Thus, it contributes only with the constant value to the creep compliance as long as the initial stresses remain small with respect to the tensile strengths. However, the creep compliance of the composites part has to include interface adhesion properties-namely, perfect adhesion. Due to its mathematical structure, it allows for the introduction of the adhesion factor having values between 0 and 1. This enables the elucidation of the time-dependent adhesion factors. This approach offers (design) engineers the chance to predict the realistic creep behavior of particulate composites for any filler volume content if the creep compliance of the matrix for any given initial stress and the adhesion factor is available.

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.