Viscoelastic model with complex rheological behavior (VisCoR): incremental formulation

Abstract A thermo-rheologically complex linear viscoelastic material model, accounting for temperature and degree of cure (DoC), is developed starting with series expansion of the Helmholtz free energy and systematically implementing simplifying assumptions regarding the material behavior. In addition to the temperature and DoC dependent shift factor present in rheologically simple models, the derived novel model contains three cure and temperature dependent functions. The first function is identified as the rubbery modulus. The second is a weight factor to the transient integral term in the model and reflects the current temperature and cure state, whereas the third function is under the sign of the convolution integral, thus affecting the “memory” of the material. An incremental form of this model is presented which, due to improved approximation inside the time increment, has better numerical convergence than most of the similar forms. Parametric analysis is performed simulating stress development in a polymer, geometrically constrained in the mold during curing and cool-down. The importance of using proper viscoelastic model is shown, and the role of parameters in the model is revealed and discussed. Graphical Abstract


Introduction
The modern composites industry is a rapidly evolving and manufacturing techniques and processes to produce optimized reinforced polymer composite parts are also being refined at a similar pace. Modelling of these processes using numerical methods is yet to mature. Prediction of process-induced residual stresses and deformations that occur during composite manufacturing requires special attention. It is vital to obtain a final composite component that matches up to specifications dimensionally and meet design criteria for strength.
On the micro-scale process-induced residual stresses occur primarily due to different behavior with temperature (T) and degree of cure (DoC) of the resin and the fiber [1]. In most major manufacturing processes, initially the resin is an uncured liquid. During the manufacturing, the resin undergoes phase changes from a liquid to a solid. These phase changes induce volumetric shrinkage of the resin as it solidifies; change its mechanical properties and thermal expansion coefficient. Together with external constraints and influences, like interaction between the tool and the part, thermal gradients, it leads to stress build-up in the composite part. This in turn causes undesired shape distortions and the final part deviating from the design specifications. Hence, simulation of these phenomena is vital in order to employ methods for compensation and corrective methods for minimization of shape distortions. For example, the post-cure time and temperature has a significant influence on the spring-in angles [2].
Over the years, several models have been developed [3][4][5][6][7][8] to understand (describe) and predict these phenomena for manufacturing techniques like resin transfer molding (RTM), autoclave processes, filament winding, and so on. The numerical methods developed so far broadly fall under two categories; one assuming that the material behavior is elastic and the other assuming viscoelastic (VE) behavior. Both approaches have their advantages and disadvantages depending on the application and available computational resources. While the number of input parameters for elastic models is small and the model can be easily implemented in a numerical code, the prediction accuracy is sometimes questionable. VE models provide a more accurate description of the material behavior and hence render better predictions. However, the characterization costs are greater, and their usage is computationally expensive. Hence, the quest challenge in the field of prediction of processinduced deformations and residual stresses has always been to develop models, ideally viscoelastic, with minimized characterization effort and low computational costs without sacrificing accuracy in the process. The VE behavior depends on temperature, but the material behavior also changes with the DoC. The combination of T-and DoC-dependence makes computation rather expensive and an analysis of shape distortions in composites using VE models is currently not widespread.
One of the still most frequently used elastic approaches is the Cure Hardening Instantaneously Linear Elastic (CHILE) model [7]; the properties are linear elastic within each material phase. However, this model is incapable to account for the so-called "frozen-in" deformations that occur during cure and procedures of their release.
Svanberg and Holmberg in [8] proposed so-called "linear elastic path-dependent model" was derived as a limiting case from the incremental form of the linear VE model formulated by Zocher et al. [9]. In fact, this model is "pseudo-viscoelastic"; the behavior in the rubbery state is assumed to have instant stress relaxation whereas in the glassy state there is no relaxation at all. Due to its simple implementation and relatively inexpensive computation costs, it has been widely used in various commercial tools for process simulations of industrial relevance [10][11][12][13].
One of the shortcomings of the path-dependent model is that it cannot account for viscoelastic relaxation in glassy state at high temperatures: there is no stress relaxation or strain recovery below glass transition temperature T g : Therefore, the strains are "frozen-in". They may be suddenly completely released when the material "jumps" to the rubbery state at temperature larger than T g : In reality, even in the glassy state at high temperature the VE relaxation is very fast and in the rubbery state the relaxation is gradual, not instant. Therefore, this model sometimes overestimates the built-in residual stresses, particularly when modelling high-temperature composites. The above discussion shows that VE models are required. Among several attempts to simulate the VE behavior, we can mention Machado et al. [14] who tested their incremental thermo-rheologically simple VE constitutive model against various thermoplastic composite materials. Benavente et al. [2] implemented a VE model in Abaqus to study relaxation during cure of an epoxy reinforced with 3 D interlock woven fabric.
In the aforementioned VE models, the DoC and T affect only the shift factor. The time shift is based on the time-temperature-cure-superposition (TTCS) principle. According to TTCS, high temperature and/or low DoC shorten the VE relaxation times. Kim and White [15] characterized the temperature shift factor at various DoCs. Simon et al. [16] used the DiBenedetto function to describe the glass transition temperature, T g as a function of the DoC. The T g was then used as the reference temperature in the temperature shift factor. O'Brien et al. [17] analysed the effect of cure on VE parameters using in experiments combination of parallel plate rheometry and three-point bending tests in a DMA.
Recent new studies [18,19] proposed to improve VE models by including DoC and T-dependence of the glassy and rubbery stiffness. In [20,21], a very general and thermodynamically consistent VE model with DoC and T dependences was developed following the derivation scheme suggested by Schapery [22]. This constitutive model allows accounting for the glassy modulus change with DoC. Ding et al. [23] have implemented complex VE behavior using a temperature dependent factor as a multiplier for Prony series. The multiplier is identified as a function of the ratio of stiffness at a given temperature to the reference temperature. Zhang et al. [24] used a very different formulation: in addition to reduced time they introduced also DoC dependent relaxation times, replacing the combination of them with one parameterthe Deborah number. They argue that the inverse of the Deborah number comprehensively reflects the effects of the DoC, the shift factors, and the relaxation times on the viscoelasticity. Using extreme values of the Deborah number leads to the model presented in [8].
Most of the above discussed models, except for [20,21], represent thermo-rheologically complex behavior as empirical, aimed to introduce flexibility in the model without checking for thermodynamic consistency.
In this paper, we present a thermodynamically consistent, rheologically complex VE model, which is derived from Helmholtz free energy, systematically using simplifying assumptions. The obtained material model accounts for the rubbery modulus and glassy modulus dependence on DoC and T, represented in the model by three DoC and T dependent functions. One the functions is the current value of the rubbery modulus, another one is a weight factor in front of the convolution integral characterizing T-and DoC state in the current instant and the last (the one under the sign of the heredity integral) defines (weights) the significance of different material development stages and temperature changes during the whole history. The material model is written in an incremental form with fast convergence (second order convergence rate), achieved by using improved approximation of functions in the time increment (the main effect comes from the improved approximation of the highly changing shift factor). Using the incremental form, the stress development in a geometrically constrained thermo-rheologically complex polymer during curing with following cool-down is simulated, and the significance of different parameters is revealed.

Thermo-rheologically complex viscoelastic material model (VisCoR)
The VE model complex thermo-rheological behavior in this paper is based on following assumptions: The cross-linking and the T changes may have the same time scale as the process time. Therefore, time derivatives of material functions in the model cannot be neglected. The material is linear viscoelastic and viscoplasticity is not included. For the sake of clarity, the presentation in this paper is one-dimensional (1D). Geometrically linear material behavior (small deformations) is assumed. Free thermal expansion and curing shrinkage strains are accounted for but are not explicitly included in derived expressions. This means that strain e in following expressions is, in fact, e ¼ e applied À e free thermal À e free curing (1) The derivation of this model, using the thermodynamics-based routine suggested by Schapery [22] for nonlinear VE materials, is outlined in Appendix 1; the Helmholtz free energy function is expanded with respect to applied strain e and internal state variables related to viscoelasticity.
With a being degree of cure and T being temperature, the derived stress-strain relationship is written as The reduced time w is defined by expression In the Prony series (3), C i are constants and s i are relaxation times. In (2), E r a, T ð Þ can be identified as the rubbery modulus; h 1 a, T ð Þ is a weight factor in front of the convolution integral and the stress value in a particular time instant depends on the value of this function in that instant; h 2 a, T ð Þ "shapes the fading memory", determining (weighing) for the current instant the importance of events that happened during the whole history. These dependencies do not exist in thermo-rheologically simple material models. Hence, assuming E r as independent on DoC and T and taking h 1 ¼ h 2 ¼ 1, equations (2)-(4) are reduced to stress-strain relationship for a thermo-rheologically simple linear VE material [9].
The presented model and its derivation are very similar to the model presented in [20]. The most general eq. (12) in [20] is identical to (2)-(4). Follow-up assumptions used in [20] leads to eq. (16), (17) that are different: a. the first term (rubbery modulus) is assumed independent on DoC and T; b. the function in front of the hereditary integral is taken under the sign of the integral, without stating that it depends on the current time instant t (not on integration variable s); c. without any justification the third function h 2 1: We will show in Section 7 that assuming for this function even a weak linear dependence on DoC, the derivative (the slope of the dependence) has significant effect on stress predictions.
Temperature dependence of these functions is not considered in [20,21] but its incorporation does not require any changes in the model.
Among other modeling efforts, we can mention [15], ref. eq. (9), where a rheologically complex model was suggested, with the rubbery modulus E r depends on DoC. In [15], all coefficients, C i in Prony series depend on DoC, which leads to glassy modulus dependence on DoC. They assumed that even the relaxation times s i depend on DoC. Finally, this model was simplified, leaving only the relaxation time dependence on DoC. In [17], O'Brien also assumed DoC dependent coefficients in Prony series, which results in DoC dependent rubbery and glassy shear moduli. They have also assumed that the relaxation times, s i depend on DoC and that this dependence is the same for all of them. Obviously, the latter dependence can be included in the shift factor. Zhang et al. [24] also assumed relaxation times dependent on DoC, but in a rheologically simple model. They combined the DoC dependent shift factor, a and the relaxation times, s i in one parameter called Deborah number. Ding in [23] assumed DoC and T dependent rubbery and glassy modulus, but with the limitation that they both have the same dependence.
The generalization in the latter papers is not systematic; i.e. replacing constants with functions dependent on DoC and T: The approach presented in Appendix 1 is consistent with thermodynamics; making a sequence of assumptions regarding the type of material and analyzing consequences. For example, introducing "reduced time" the change of all VE internal state variables is described by differential equations with constant coefficients. Solving these equations, we obtain solution in Prony series with coefficients and relaxation times that are constants. The potential physical relaxation time dependence on DoC is already included in the shift factor which defines the reduced time.
Arbitrarily replacing some constants with functions does not ensure consistency with basic principles and models may be unreliable in more complex cases.

Derivation of the incremental form
One path to overcome the storage problems in numerical calculations has been by replacing the integral representation with incremental form where recursive expressions have been derived for calculation of parts of the integral. This approach, followed also in the present paper, leads to significant memory savings because calculated values only in two last consecutive time instants have to be saved.
The second path [20,21] is solving directly the very basic differential equations obtained from principles of thermodynamics. These equations, that are used deriving the presentation with the heredity integral, are solved using different finite difference methods [25].

Stress increment
The list of notations used in the incremental formulation is given in Appendix 2. We substitute in (2) the transient compliance expression (3) and rewrite (2) for time instants t k and t kþ1 The integral in (6) can be written as a sum of an integrals from 0 to t k and from t k to t kþ1 Subtracting (5) from expression (7), we after some rearrangements obtain In (8), the following notation is used S k i and DH k i defined by (9, 10) are analysed in Sections 2.5 and 2.6. Closer investigation of expression (11) shows that Q kþ1 i is not an independent function Using (12) and (13) (8) can be written in a different form This form will be used in all following derivations.

Assumption about parameter change in the interval
We assume that in the small interval Dt kþ1 the VE parameters are linear functions of time s: The shift factor a is given by We assume linear dependence of We also assume a linear time dependence of the strain change.
In (16) and (17), H is the Heaviside step function. The linear approximation used in each increment leads to calculation scheme that converges much faster, see Section 6.2, than common schemes where constant values of parameters in the increment are used.
In contrast to most of other representations, the time increment does not have to be constant. For example, large regions where the parameter change and the strain change is linear can be treated in one-time step. Zocher et al. [9] used slightly different incremental formulation presenting formulas for constant reduced time increments instead of time increments, assuming that the strain change is linear with respect to the reduced time. This form of presentation has its advantages because it automatically adjusts the time step to the shift factor value: in the rubbery state with low values of the shift factor the time step becomes very low. Better accuracy is automatically achieved, but one could question whether it is necessary in the region before gelation and in the rubbery region where stresses are usually very low. Apparently, the selection depends on the aim of the simulation and the formulated targets.

Expression for Dw kþ1
Assumption that Dw kþ1 is a linear function of Dt kþ1 is not correct for finite increments and assumptions made in Section 2.2. According to (4) Substituting (15) in (18) and integrating, we obtain Direct substitution of d kþ1 a defined by (15) in (19) yields Expressions (19) and (20) are not suitable for numerical calculations in regions where a is not changing or changing slowly (d kþ1 a 0). Expansion of ln 1 þ x ð Þ% x can be used leading to Keeping also the second term in expansion of (20), a more accurate expression is obtained 2.4. Linear approximation of h 2 , e and the derivative of their product The integral in (2) contains the term d ds h 2 e ð Þ: From (16) and (17) (keeping in mind that HÃH ¼ H) The derivative of the Heaviside step function is the Dirac's delta-function. Hence, The numerical value of an integral involving delta-function is equal to the function under the sign of the integral at s ¼ t k : This means that after integration all terms containing s À t k will give null. The remaining (relevant) terms are Substituting expression (25) into (26) we obtain Expressions for integrals DI kþ1 i and DN kþ1 i are derived in Appendices 3 and 4.

Recursive expression for
can be written as a sum of two integrals (29) The second integral in (29) is DH k i defined by (27) and (A3.8), (A4.3). According to definition (10), the first integral in (29) is S kÀ1 i : Therefore (29) can be rewritten as

Incremental form for implementation in software
The equations described in Section 3 can be rewritten in a different form: In (38), In (39), Dr k R depends on strain in time instant t k and before.
The reduced time is calculated using (32), DÎ

Numerical instability
For small values of Dw k s i numerical accuracy problems can appear when calculating 1 À exp À Dw k s i : Therefore, if Dw k s i < d where d is small (for example 0.000001), we have to use first expansion terms in 1 À exp Àx ð Þ % x À x 2 2 : This approximation has to be used in (39). In addition, substituting this expression in the above expressions (36) and (37) we obtain If in addition to the above, Da k a kÀ1 is also very small, approximate expressions (21) and (22) can be used to calculate Dw k leading to

Results and discussion
In this section, we will focus on verification of the VE constitutive model and analyze results for one test case to provide insight into the behavior of materials that follow the proposed constitutive relation.

Material properties
We have used mechanical properties of the epoxy system LY5052, see Table 1, used for simulations in [8]. The input data and the temperature profiles used are not critical for the discussion in this paper. Parameters for DC in (3) are given in Table 2. The values are computed using weights taken from [21].
The temperature during curing is constant, T ¼ 120 C, and the DoC is a linear function of time, with a ¼ 1 at t ¼ t 1 (t 1 ¼ 100 sec).
The glass transition temperature, T g does not explicitly enter the used model: parameters in the model depend on DoC and temperature. To gain an additional insight, the development of the T g with time is calculated using which is obtained from the DiBenedetto expression. Figure 1 presents the temperature ramp and the calculated development of DoC and T g : The shift factor is presented as a product of temperature and DoC dependent functions [18,19] a T, The temperature shift factor a T ¼ 1 at T ref ¼ 30 C; the DoC shift factor a C ¼ 1 when the material is fully cured. Hence, the shift factor (47) is equal to one for a fully cured material at 30 C: The following dependences of parameters are assumed in the rheologically complex model (for simplicity referred to as "VisCoR" in contrast to the "simple" model which means "rheologically simple") In (51)-(52), k and m are constants. These forms are arbitrarily assumed, the aim is to illustrate the sensitivity of calculated stress distributions with respect to these parameters.

Loading case
We analyze the curing of a rod of a neat resin, which during curing and the following cool-down is constrained at the ends, i.e.
If not constrained, free chemical shrinkage and thermal expansion following the temperature change would take place during curing. These changes are assumed linear with respect to DoC and T, with values of coefficients in different regions given in Table 1. From (1) and (53) follows These conditions are close to conditions when a polymer is cured in a rigid mold with the difference that our simulation will be one-dimensional. At t gel ¼ 50 sec, when the degree of cure is 0.5 (gelpoint according to Table 1), the polymer transforms from a liquid to a rubber like solid. From this point, tensile stresses start to form due to constraint on chemical shrinkage, but since relaxation in the rubbery state is very fast, the stress is quite small. This is not necessarily the same when the constraint is 3-D, when the bulk modulus, which may be large even in the rubbery state, is the stress-governing property [9].
When the degree of cure increases the glass transition temperature, T g also increases according to Eq. (46). At t vit ¼ 96 sec, T g reaches the cure temperature and the polymer transforms from a rubber state to glassy state. In our VE model, this transition does not require sudden jump-wise changes of properties: they are changing with respect to T and DoC according to defined laws. The chemical shrinkage coefficient does not change in the glassy state, but the stress relaxation is much slower and, therefore, stresses from 96 sec increase faster. At t 1 ¼ 100 sec the rod is fully cured, and the glass transition temperature is 136 C: From t 1 until 150 sec, the temperature is reduced with constant rate from 120 to 20 C: Since thermal shrinkage is not allowed, tensile thermal stresses are rising (they are partially relaxed when the temperature is still high and the relaxation is fast). The stress Figure 1. Temperature ramp during the test, increase of glass transition temperature (T g ) and degree of cure (DoC).

Convergence studies
Incremental formulations have partially resolved the storage problem, but the computation time (the convergence rate of different schemes with respect to the time increment) is an additional problem. The convergence rate depends on the used approximation regarding the behavior of the DoC and T dependent nonlinearity functions in the time increment. In our incremental scheme, the function h 2 a, T ð Þ under the sign of hereditary integral is rather slowly changing and a linear dependence on time t in the increment is assumed.
However, this is not the major part of improving the convergence rate: the change of the reduced time w in the increment with monotonously changing T and DoC is much faster. The linear approximation of the shift factor and analytical integration, used in the paper, leads to logarithmic expression (20) for the reduced time increment. Using (20), which leads to quadratic convergence, instead of having the reduced time constant in the time increment, see example in [27], the time increment could be increased 100 times keeping the same relative error.
It is noteworthy, that a comprehensive analysis of convergence of incremental schemes with several different approximations of the shift factor was presented in [25]. Instead of rewriting the Schapery's viscoelastic model in incremental form as we did,  they took one step back and solved numerically the basic differential equations from thermodynamics analysis leading to Schapery's model. They have used a 4 th order Runge-Kutta finite difference scheme. It was shown in [25] that assuming the shift factor as piecewise constant over the time increment t 2 t k , t kþ1 ½ , the convergence rate is improved comparing with a rougher assumption that the shift factor is constant [26] (in [26] the reduced time w becomes a linear function of time in the increment). The refined approximation leads to quadratic convergence rate. The described studies were for stress dependent shift factors, not for DoC and T dependency, where dependence is significantly stronger.
The following example, using the VisCoR model for a "simple" case, h 1 ¼ h 2 ¼ 1, E r ¼ 28 MPa, explains some of the features of the approximations in the incremental formulation. The "simple" case is sufficient due to features of the VisCoR model and specifics of the loading case considered. In the example analyzed, because the aforementioned parameters are kept constant, they do not affect convergence rate: a. the stress in a certain time instant depends on the values of E r and h 1 in this instant, not on how these functions were changing with time making the time increment irrelevant. b. the model assumes linearity of h 2 , inside the increment. For a slowly changing function, this approximation does not require small time increments. More than that, for the assumed material properties we used (see (52)), a linear temperature dependence of this function and linear T change, makes the linear approximation exact making the value of the time increment irrelevant.
The loading case described in Section 6.1.2 is simulated, and the stress increase curve is shown in Figure  2. The present model assumes linear change of the shift factor in the increment, see Section 2.2, with a slope which is based on shift factor values at both ends of the increment. This approximation gives much faster convergence in simulations than assuming it constant. According to Table 3, where the peak values of the calculated stress are shown, taking Dt k ¼ 0:1 sec in the VisCoR model leads to higher accuracy than Dt k ¼ 0:01 sec used in the formulation in [8,10,21] where the shift factor value in the beginning of the increment is used: the number of time steps can be significantly reduced using the present formulation. Most of the published incremental procedures are formulated for constant time step Dt k ¼ Dt (except from Zocher's model [9], which is formulated for a constant reduced time step Dw k ¼ Dw). In the present model, the time step can be arbitrary changed during the simulation. For example, if we use Dt k ¼ 10:0 sec in the rubbery region and Dt k ¼ 0:1 sec in the glassy region, the predicted maximum stress is 11.75468 MPa, which is the same as using constant time step Dt k ¼ 0:1 sec over the whole region.
The path-dependent model [8] is in fact just a particular case of the VisCoR model with assumed extreme values of the total shift factor given by This implies that in the rubbery region (t < t vit ) when a ! 0, exp À Dw kþ1 s i ! 0 and DÎ k i ! 0: In the glassy region (t ! t vit ) when a ! 1, Dw kþ1 ! 0 and DÎ k i ! 1: Thus, the equations (35) and (39) can be simplified and rewritten as Equations (56) and (57) are exactly as described by the path-dependent model in one dimension and hence proving that the path-dependent model is in fact a limiting case of the VisCoR model.
The VisCoR model was used with the shift factor in (55) to reproduce results in [8] for the same loading case. A large difference was obtained in stress between the VisCoR model in "simple" case and the path-dependent model [8], see Figure 2a. The stress in the latter is significantly overestimated, the reason being an absence of any stress relaxation in this model in the glassy state. In fact, even in the glassy state, if the temperature is high, the stress relaxation is very fast and during cool-down in the range from 120 C to 80 C the stress is not at all proportional to the temperature change. The evolution of the parameters DH k i , S k i and DC k in the VisCoR model formulation in the simple case are also shown in Figure 2  VisCoR model with assumptions formulated in Section 2.2. In following figures, one curve, called "simple", is obtained assuming rheologically simple behavior of the polymer in the VisCoR model (h 1 ¼ h 2 ¼ 1, E r ¼ 28 MPa). Other curves in the same figure will correspond to rheologically complex behavior with only one of the three parameters changing during the simulation.
6.3.1. Effect of rubbery modulus dependence on DoC and temperature First, see Figure 3, we analyze the influence of the rubbery modulus E r which increases with DoC and reduces with increasing temperature according to the assumed law in (50). It is assumed that due to DoC change the rubbery modulus can increase by 50% and due to 100 C temperature increase E r would decrease by almost 30%.
In the VisCoR case, stress is systematically lower than in the "simple" case. This is understandable, since E r in the VisCoR case is always smaller than in the simple case. In the final stage of the cooldown, E r is almost the same as in the simple case (fully cured material at low temperature). Recalling that in eq. (2) E r is not linked with the convolution integral, we can expect that E r affects stress at a certain instant only through its value at this instant.
Indeed, the relative difference between predictions in the simple and complex case is the largest in the beginning of the stress build-up in the rubbery region. The stress difference is the largest (0.114 MPa) when the T g equals the curing temperature (t vit ¼ 96 s). The difference holds until t ¼ t 1 ¼ 100 s. After that, stresses are rising due to the cooldown. Then the difference between the two modelling cases reduces and at t ¼ 150s both values coincide (11.75456 MPa), confirming the above conclusion that only the instantaneous value of E r affects stresses at that particular instant. Since this parameter is relatively small, its effect on stresses in the glassy state at room temperature is almost negligible.

Effect of the scaling factor h 1 on viscolelastic behavior
The next parameter in the VisCoR model is h 1 T, a ð Þ: According to assumption (51), this parameter at high temperature has lower values than one: assuming 100 C temperature increase and k ¼ 0.002, h 1 will decrease by 20%. At very low DoC the h 1 is lower by 20%. We assume now that it is the only parameter changing (h 2 ¼ 1, E r ¼ 28 MPa). Since h 1 T, a ð Þ is a multiplier to the integral in (2), we can consider it as an instantaneous "weight factor" to all transient viscoelastic effects developed in the "simple" modelling case during the loading history. This interpretation explains the equality of the stress according to the "simple" and the "VisCoR" modelling in the final point (the value of the integral is the same in both models and h 1 ð150sÞ ¼ 1), see Figure 4. This type of behavior (without explanation) was observed in the path-dependent code in [21], where the same DoC and temperature dependent factor was implemented for the rubbery modulus and for all terms in Prony series. There is almost no effect of h 1 on stresses at t < t 1 : this parameter "weights" transient effects and, since in the rubbery region stress relaxation is very fast, slight variation of the h 1 function does not change the stress.
In Figure 4, in the time interval t > t 1 , the polymer is fully cured and a ¼ 1: Hence, h 1 changes only because of the temperature change. In a different thermo-mechanical loading ramp, for example, with cool-down before the full cure state is reached, the stress peak at T ¼ 20 C would be lower than in Figure 4. Reasons: a) since even in the final instant h 1 < 1, we may expect, see (51), lower stress than in the "simple" modeling in Figure 4; b) since a C < 1, the stress relaxation during the cool-down would be faster than in the "simple" case, additionally reducing the final stress.
6.3.3. Effect of h 2 T, a ð Þ on stress build-up Finally, we consider h 2 T, a ð Þ, which is the main generalization of the model in [20], changing it according to (52). At high temperature, this parameter is smaller: assuming 100 C temperature increase, the decrease of h 2 would be 10-50% for values m ¼ 0.001; m ¼ 0.003 and m ¼ 0.005. According to (52) the DoC can change h 2 by 20%. The meaning of this parameter and the temperature and DoC dependence of it is currently rather unclear and the results below may shed some light on it. In a test where DoC and T would be kept constant, h 2 ¼ const < 1: Then this parameter can be taken outside the integral (2) and its role is exactly the same as the role of h 1 described above: reduced peak stress in Figure 4. The situation dramatically changes, see Figure 5, if this parameter is changing with time according to (52). Results in Figure 5 show that in a fully cured state the derivative of h 2 , represented by coefficient m in (52), has an extremely strong effect on stress, even if the values of h 2 T, a ð Þ < 1 are not changing much. The peak value of the stress increases with m, almost two times exceeding the value predicted by the "simple" model. It has to be stressed that this T and DoC dependent function in thermorheologically complex material models for manufacturing process simulation has never been introduced and discussed before. There are no experimental data available proving or disproving its necessity in the model for a particular material and the functional dependence on T and DoC is entirely unknown. For, example, if with decreasing T it would decrease, the effect on stress may reverse: the stress would be lower than in the "simple" behavior model. Tests can be designed to study this parameter in the material model for a particular polymer. In designing a proper test, simulations of expected response are important. Example could be the effect of this parameter in relaxation test with changing temperature or with changing DoC.

Conclusions
A viscoelastic (VE) material model describing complex thermo-rheological behavior dependent on temperature and degree of cure (DoC) has been developed from expansion of the Helmholtz free energy in terms of strain and viscoelastic internal state variables. Systematic sequence of assumptions regarding material behavior has been applied, resulting in a model (for clarity presented in this paper in 1-D formulation) which, in addition to shift factor introduced as in a thermo-rheologically simple model, contains three temperature and DoC dependent functions. One of them can be interpreted as a modulus in rubbery state, E r : One of the other two, h 1 , is a weight factor in front of the convolution integral and it represents the current instant, whereas the last one, h 2 , which has never been introduced in other models is under the sign of the integral and its derivative with time affects the fading memory of the material.
To allow its implementation in commercial FEM codes, this model is rewritten in an incremental form. A more accurate than usually interpolation in the time increment is used, reaching high accuracy with relatively large time steps. The time step can be different in different stages of cure, temperature and loading.
Stress build-up in a polymer which is cured and then cooled down while being geometrically constrained in the mold is simulated showing that a proper VE model even in a thermo-rheologically simple formulation predicts much lower stresses that the very popular quasi-VE path-dependent model [7]. In the analyzed example, the rubbery modulus E r affects stresses only in the rubbery region but not the final stress value in the glassy state at room temperature. Similarly, the function h 1 has a limited effect on stress development as long as the temperature is high and/or the cure is not complete. However, the predicted stress is very sensitive with respect to the last function h 2 and its derivative.