Incremental Numerical Approach for Modeling the Macroscopic Viscoelastic Behavior of Fiber-Reinforced Composites Using a Representative Volume Element

Abstract

The objective of this study is to describe the stress relaxation behavior of an epoxy-based fiber-reinforced material. An existing incremental formulation of an orthotropic linear viscoelastic material behavior was adapted to Voigt notation and to the special case of an isotropic material. Virtual relaxation tests on a representative volume element were performed, and the behavior of individual components of the relaxation tensor of the transversely isotropic composite material was determined. The study demonstrated that in the case of only one viscoelastic material, each component of the relaxation tensor can be described in terms of a scalar form factor and the behavior of the neat resin. The developed method was implemented in an incremental finite element model (FEM) analysis to calculate the stress relaxation on the macroscopic ply level. A validation of the approach has shown a promising agreement up to a limit below the glass transition temperature of 15 °C in longitudinal and 35 °C in transverse direction. This study therefore demonstrates a novel way to incrementally describe the macroscopic viscoelastic behavior of materials with a single viscoelastic component with good controllability for engineering purposes.

1. Introduction

Fiber composite materials are known to have the potential to replace other materials in areas subjected to high mechanical loads. These materials can be particularly attractive for saving weight in moving masses. However, their potential can only be exploited if knowledge gaps in production and use are closed. This was the aim of the project within which this study was conducted. The research activities were focused on the influence of residual stresses on the fatigue strength of a composite leaf spring for the automotive industry.

During the manufacturing process of a thermoset material (e.g., compression resin transfer molding CRTM), its structure is exposed to chemical shrinkage and thermal contraction. Depending on the geometry of the part, especially if free deformation is restricted, this can lead to the formation of process-induced stresses. Analogous to a mechanical relaxation test, these stresses will relax with time, especially at higher temperatures and a low curing degree. In a typical manufacturing process for fiber-reinforced polymers, the material is in exactly this state at the beginning of the process, which means at comparably high curing temperatures compared to the actual glass transition temperature. When one estimates the extent of the residual stresses, the stress relaxation has to be considered in order not to make an overly high and thus conservative prediction.

For this reason, the introduction of a homogenization model to describe the relaxation effect of the composite based on the behavior of its individual components is advantageous. Isolated measurements of the components of the relaxation tensor are difficult to perform since a unidirectional strain is challenging to produce in a mechanical test. In contrast, an isotropic material can be described in terms of only two parameters, so the elimination of the lateral contraction is not necessary. Therefore, the relaxation behavior can also be measured in a multiaxial strain state, e.g., uniaxial tension. Furthermore, the introduction of a homogenization model from an isotropic to an orthotropic material allows the consideration of different fiber volume contents.

Within the current project, a virtual process simulation was set up based on earlier models [1,2], in which the modeling of the relaxation effect was the last in a series of implemented physical effects. It was assumed that for some loads, a lack of knowledge could lead to significant inaccuracy in the predicted final performance. To remedy this knowledge gap, the present study was performed, which investigated the magnitude of relaxation of isolated loads on microscopic unit cells. Its aim was to derive a method that can adapt to varying properties and contents of the individual components of a composite and is able to provide input values for the manufacturing process simulation implemented in an incremental FEM analysis.

1.1. State of the Art

The available approaches used to model process-dependent residual stresses in fiber-reinforced materials show a trend in which increasing accuracy comes with increased complexity and computational effort.

The most drastic simplification of the manufacturing process is an elastic model that considers the total thermal and/or chemical shrinkage together with the stiffness at a reference temperature, usually room temperature. The amount of thermal shrinkage is thus calculated as the difference between the reference temperature and a defined stress-free temperature (usually the curing temperature) multiplied by a coefficient of thermal expansion [3]. Since the material’s stiffness at the curing temperature is lower than at room temperature, this approach obviously leads to an overestimation of the stresses that occur in reality. This has the beneficial effect that an assessment of a manufactured part is conservative, which is generally required for engineering purposes.

To increase the accuracy, incrementally elastic models, for example, the CHILE (cure hardening instantaneously linear elastic [4]) model, have been developed. They are especially popular due to their ability to be easily implemented in an incremental FEM analysis [5,6]. This approach derives the occurring stress from the integration of the occurring strain with the actual stiffness over time:

(1)$\sigma \left(t\right)={{\displaystyle \int }}_0^{t}E\left(T,\mathsf{\alpha }\right)\frac{\mathrm{d}\epsilon }{\mathrm{d}\mathsf{\tau }}\mathrm{d}\mathsf{\tau }$

A further increase in the accuracy of the predicted stresses can be achieved if the viscoelastic effects are considered. It was found in [7] that the elastic approaches overestimate the occurring stresses up to 20% in some specific cases.

In linear viscoelastic approaches, the elastic modulus is replaced with a time-dependent modulus so that the formulation becomes:

(2)$\sigma \left(t\right)={{\displaystyle \int }}_0^{t}E\left(t-\tau \right)\frac{\partial \epsilon }{\partial \mathsf{\tau }}\mathrm{d}\mathsf{\tau }$

This integral, known as the convolution integral, takes an easier form if it is transformed into the Laplace domain. An example of viscoelastic approaches that take advantage of the Laplace transform can be found in [8], where a material model for woven fabrics was developed. The validation was carried out for the distortion of a plate with an asymmetrical layout. Other approaches, especially for the case of epoxy resin systems, were developed in [9,10].

A drawback of the methods involving Laplace transformation is that they are computationally expensive and because, in some cases, no analytical solution for the back-transformation can be found. To avoid these problems, approaches in the time domain have also been developed. Zocher et al. [11] developed an incremental formulation for a linear viscoelastic material whose stiffness tensor entries can be described in terms of the rheological model known as the generalized Maxwell model. Since this model is known to provide a good approximation for the relaxation behavior of a wide range of materials and at the same time it is possible to be implemented in an incremental FEM analysis, this approach is very promising. An example of the application of this formulation can be found in [12], where the material was modeled to relax instantly in the rubbery phase and show almost elastic behavior in the glassy phase during the curing of an epoxy-based thermoset composite material.

The presented viscoelastic approach was integrated in a holistic approach of a virtual manufacturing process [1]. It aims to include all relevant physical aspects leading to process-induced stresses and consists of a thermochemical and a subsequent mechanical part. As mentioned earlier, one area of uncertainty is the so-called virtual material characterization of the relaxation properties of a composite material from the properties of its individual components. While there are well-established concepts for elastic homogenization, there are few methods available for viscoelastic homogenization that are user-friendly and allow a fast assessment of the accuracy. In the following, we first present the available methods for homogenization of the elastic properties and then those for the viscoelastic properties.

1.1.1. Homogenization of Elastic Properties

In comparison with viscoelastic behavior, there are several widely used homogenization concepts for elastic behavior. They can be differentiated into empirical, analytical and numerical approaches.

A well-known empirical relationship was developed, namely the Halpin–Tsai equations [13,14]. All the mechanical properties can be described in terms of two equations as follows:

(3)$\frac{p}{p_m}=\frac{1+\xi \eta v_f}{\left(1-\eta v_f\right)}{E}_m where \eta =\frac{p_f/p_m-1}{p_f/p_m+{\xi }^{\prime }}$

where p is an arbitrary property of the composite, p_m and p_f denote a property of the matrix and fiber, respectively, and the parameter ξ′ is an adjustment parameter that needs to be determined to fit to the experimental data. In order to correct its accuracy for shear stiffness, ref. [15] formulated a correction factor for this formulation.

Analytical approaches usually use drastic simplifications of the composite material. A well-known approach is the rule of mixtures, which calculates the homogenized stiffness from differently arranged springs representing the stiffness of the matrix and the fiber, respectively. By arranging the springs in parallel, it is possible to describe the fiber-dominated properties such as the stiffness in the fibers’ direction, and arranging the springs in series can be used to describe matrix-driven properties such as the transverse stiffness or shear stiffness. A more accurate micromechanical model based on this concept is known as the rule of mixture developed by Chamis [16], which rose to eminence in the aerospace industry.

In terms of numerical homogenization, three approaches are commonly used for simulating composites: a micromechanics-based approach, an equivalent homogeneous material (EHM) based approach and a combination of the first two methods. Each approach has advantages and disadvantages described in the literature [17,18,19,20]. The micromechanical approach describes the material structure in detail, and it is possible to investigate local defects [18]. However, the computational effort is very high, because a much better mesh quality is required in comparison to the EHM model. The EHM method reduces computational effort by homogenizing the microscopic properties on the macro level, but it neglects local effects [21,22]. It is possible to combine the two models [18,20,23] to receive the advantages of both.

Microscale modeling involves consideration of uncertainties in the composite structure, which requires stochastic models [24]. In microscale models the composite properties are defined from all constituent materials. Basic approaches include elasticity theories based on the repetition of a representative unit cell, assuming a perfect bond between fibers and matrix. However, an advanced viscoelastic models could be implemented to describe fiber-media interaction and the fiber packing process [25]. The focus of the micro-scale study is on the fiber composition, geometry, and orientation. Such micromechanical theories are usually validated with experiments [24,26,27].

Multiscale techniques solve local problems on the microscale and take this information to the macro level through homogenization techniques [28]. The material properties of the macroscale model could be obtained from micromechanical analyses. Mechanical behaviors of heterogeneous materials are often described with representative volume elements (RVE) [29]. Hill’s theory assumes that the RVE should contain enough fibers to provide a statistical representation of the heterogeneous materials. One advantage of the RVE method is that the composite material replaced by a homogeneous media retains the anisotropic material properties. The effective properties derived from the RVE represent the material on the macroscopic scale, which is commonly known as the micro–meso–macro principle [30].

1.1.2. Homogenization of Viscoelastic Properties

Consideration of the time-dependent behavior of the individual components and their interaction brings additional complexity into the homogenization of a composite‘s properties. A common approach uses the principle that if the convolution integral in Equation (2) is transformed into the Laplace space, it takes the form of basic multiplication. Afterwards, a micromechanical homogenization approach (as well as classical methods from elastic homogenization) can be applied for a specified fiber volume content:

(4)${\widehat{\Gamma }}_c\left(s\right)=MM\left({\widehat{\Gamma }}_f\left(s\right),{\widehat{\Gamma }}_m\left(s\right),v_f\right)$

Here, MM denotes the micromechanical homogenization model, $\widehat{\Gamma }$ the viscoelastic response in Laplace space of the fiber, matrix, and composite respectively and $v_f$ the fiber volume fraction.

Another common approach for the homogenization in Laplace space is known as the asymptotic homogenization method (AHM). It considers the composite as a periodic assembly of RVE’s. One of the first applications of this method was carried out by [31]. With a focus on the individual indices of the creep or relaxation tensor, an early evaluation was presented in [32], where a micromechanical model was also developed directly in the Laplace domain. A more recent example of this method was given by [33], where a Prony’s series was used in order to model the individual phases’ relaxation in the time domain. The results for low fiber volume contents can be found in a later study [34].

As mentioned earlier, the main challenge of the application of Laplace transform is to obtain analytical formulas for the inverse Laplace transform to generate the viscoelastic solution back in the time domain. Therefore, a drawback of this method is the mathematical complexity which limits its wide application by different users.

For this reason, other researchers have applied homogenization in the time domain. In [35], quasi-elastic solutions were developed for cases of elementary loads summarized in a k-factor for each of them. While these factors can be used to describe the time-dependent behavior of the creep tensor indices, they cannot be used for the relaxation tensor indices. This limitation arises because there are no widely used elementary load cases for producing a unidirectional strain, since contraction in the lateral direction is hard to restrict in a mechanical test. However, it is remarkable that some indices increase over time, which is also shown in the current investigation. A similar approach was recently followed by [36], who decomposed the strain tensor to represent characteristic load cases. A solution for determining the stress relaxation of the composite material for these load cases was developed and validated. Although it is a very promising approach, it is hard to implement due to the high number of parameters and solution uncertainty for specific load cases.

1.1.3. Limits of the Linear Viscoelastic Assumption

The assumption that linear viscoelastic material behavior describes the homogenized composite’s behavior might be inaccurate because local stresses are much higher on the microlevel compared with the meso-level or ply level. On the microscale, the material shows greater stress by a factor of 3–12 depending on the stiffness ratio between the polymer and the fiber. It might exceed linear behavior (an assumption for the application of linear viscoelasticity) even at low load levels.

This is especially visible in the transverse direction. Brauner et al. [37] simulated the stress distribution in a unit cell by considering thermal and chemical shrinkage, temperature-dependent effects, micro-yielding, and degradation and relaxation. An example considering just thermal strain is shown in Figure 1. An applied nominal strain of 0.48% can lead to a maximum local strain of 2.93%, which is an amplification of about 6.

1.2. Approach

In this study, an incremental formulation of the linear viscoelastic response [11], was used to analyze the response of isolated components of the orthotropic relaxation tensor with the use of an RVE.

It has been suggested that in the case of only one time-dependent material, the response can be described as proportional to the response of that material because the decay time may not change. For this reason, the possibility of expressing the resulting time-dependent behavior of the individual indices in terms of a scalar form factor, together with the behavior of the neat resin system, was investigated.

These factors were integrated in [11] in order to provide an incremental description of the composite material. The presented approach was implemented by using user subroutines in the commercial finite element program Abaqus (6.13, Maastricht, The Netherlands). A validation was performed by measuring an epoxy-based glass fiber reinforced composite material, demonstrating its simplicity and computational advantage. Even if the linear viscoelastic region might be exceeded at the microscale, the suitability of the homogenization approach assuming this behavior can be assessed by the validation on the macroscale.

The present study therefore provides a novel approach for the description of a composite with only one linear viscoelastic phase, which is applicable to a wide range of materials. The form of expressing the elements of the relaxation tensor in terms of form factors offers a good controllability because each factor indicates the deviation of the response of the element to which it belongs to the response of the neat resin system. The incremental form of this description in the time domain also represents an extraordinarily efficient variant compared to established methods using Laplace transformation, since only six status variables have to be transmitted for each time step.

2. Materials and Methods

The research flow of the present work was structured according to the proposed approach: First, the formulation of [11] was adapted for implementation in a commercial finite element program (Section 2.1). A description for the orthotropic case has been simplified to the isotropic case, which can represent a viscoelastic resin system (Section 2.2). The elastic and viscoelastic properties of the resin system were determined in a three-point bending test so that this data could be used to calibrate the viscoelastic phase of an RVE (Section 2.4). Next, the response of the elements of the relaxation tensor to isolated load cases was determined and it was evaluated whether these can be described using a factor and a master curve (Section 2.3). Since this turned out to be promising, a formulation for describing the macroscopic behavior from these factors was developed and a validation was carried out to estimate the accuracy of the approach. For this purpose, a virtual copy of the 3PB experiment was created (Section 2.5) and compared with experimental data from the composite (Section 2.4).

2.1. Incremental Orthotropic Linear Viscoelastic Formulation

The general stress–strain relationship for an orthotropic linear viscoelastic material can be expressed by the constitutive equation:

(5)${\sigma }_{ij}\left(x_k,\xi \right)={{\displaystyle \int }}_0^{\xi }{C}_{ijkl}\left(x_k,\xi -{\xi }^{\prime }\right)\frac{\partial {\epsilon }_{kl}\left(x_k,{\xi }^{\prime }\right)}{\partial {\xi }^{\prime }}\mathrm{d}{\xi }^{\prime }$

This respects the history of the material by assuming that each response to a loading event can be superposed on all other responses. A common method for the description of the relaxation behavior of a linear viscoelastic material is to use a rheological model. One such model known to be suitable for describing the relaxation behavior of a wide range of materials is named the generalized Maxwell model, also known as Wiechert model [38]. It consists of a number of so-called branches, each consisting of a spring and a dashpot, which are connected in parallel (Figure 2).

The generalized Maxwell model has the beneficial property that the stress it describes during a relaxation test takes the form of a Prony series. Now, by assuming that the behavior of the elements of the stiffness tensor in Equation (5) can be represented in terms of a generalized Maxwell model with M branches, its value at time ${\xi }_{n+1}$ (following the notation of [11].) can be expressed as a Prony series as follows:

(6)$C_{ijkl}\left(x_k,{\xi }_{n+1}-{\xi }^{\prime }\right)=C_{ijk{l}_{\infty }}+{\displaystyle \sum }_{m=1}^{M_{ijkl}}{C}_{ijk{l}_m}{\mathrm{e}}^{-\left({\xi }_{n+1}-{\xi }^{\prime }\right)/{\rho }_{ijk{l}_m}}$

Comparing this description with the rheological model, one can see that spring elements that are connected to a dashpot contribute the less to the actual stress the longer it has been that strain was applied. This is commonly known as fading memory behavior.

Equation (6) was inserted in Equation (5) and converted to an incremental form by [11]. That means that the actual stress is the sum of the stress increment and the stress of the previous time step:

(7)${\sigma }_{ij\left(t+\Delta t\right)}={\sigma }_{ij\left(t\right)}+\Delta {\sigma }_{ij}$

It consists of an immediate increment, a stress increment due to thermal strain and a stress increment ($\mathsf{\Delta }{\sigma }_{ij}^{\mathrm{R}}$) that considers relaxation of the stress. For small time steps, the immediate increment can be understood as the elastic increment and $\mathsf{\Delta }{\sigma }_{ij}^{\mathrm{R}}$ as the viscoelastic part.

(8)$\mathsf{\Delta }{\sigma }_{ij}=C_{ijkl}^{\prime }\mathsf{\Delta }{\epsilon }_{kl}-{\beta }_{ij}^{\prime }\mathsf{\Delta }\mathsf{\Theta }+\mathsf{\Delta }{\sigma }_{ij}^{\mathrm{R}}$

In contrast to [11], here the temperature load considered ($\mathsf{\Delta }\mathsf{\Theta }$) was integrated with the mechanical load, $\mathsf{\Delta }{\epsilon }_{kl}$. The simplified incremental description of linear viscoelastic behavior then takes the form

(9)$C_{ijkl}^{\prime }\equiv C_{ijk{l}_{\infty }}+\frac{1}{\Delta \xi }{\displaystyle \sum }_{m=1}^{M_{ijkl}}{\eta }_{ijk{l}_m}\left(1-{\mathrm{e}}^{-\Delta \xi /{\rho }_{ijk{l}_m}}\right) \left(\mathrm{no} \mathrm{sum} \mathrm{on} i,j,k,l\right)$

(10)$\Delta {\sigma }_{ij}^{\mathrm{R}}=-{\displaystyle \sum }_{k=1}^3{\displaystyle \sum }_{l=1}^3{\displaystyle \sum }_{m=1}^{M_{ijkl}}\left(1-{\mathrm{e}}^{-\Delta \xi /{\rho }_{ijk{l}_m}}\right)S_{ijk{l}_m}\left({\xi }_n\right) \left(\mathrm{no} \mathrm{sum} \mathrm{on} i,j,k,l\right)$

(11)$\begin{array}{cc}\hfill S_{ijk{l}_m}\left(x_k,{\xi }_n\right)=& {\mathrm{e}}^{-\mathsf{\Delta }\xi /{\rho }_{ijkl}m}{S}_{ijk{l}_m}\left(x_k,{\xi }_n-\mathsf{\Delta }\xi \right)\hfill \\ & +{\eta }_{ijk{l}_m}{R}_{\epsilon ,kl}\left(1-{\mathrm{e}}^{-\mathsf{\Delta }\xi /{\rho }_{ij{k}_m}}\right) \left(\mathrm{no} \mathrm{sum} \mathrm{on} i,j,k,l\right) \hfill \end{array}$

where:

(12)$R_{\epsilon ,kl}\equiv \mathsf{\Delta }{\epsilon }_{kl}/\mathsf{\Delta }\xi$

The variable S represents a state variable that must be stored in the previous time step and passed to the current step. For the orthotropic case, each branch m of the Maxwell model requires a state variable for each index of the stress tensor ij. This leads to M × 36 state variables, which is cumbersome for practical use.

For presenting the result of fitting the rheological model to the measurement data, it is beneficial to introduce a parameter that describes the relationship of one spring element to the sum of all spring elements:

(13)${\mathsf{\phi }}_{\mathrm{m}}=\frac{k_{\mathrm{m}}}{k_{\mathrm{tot}}}$

2.2. Adaption to Isotropic Material Behavior

The homogenization for modeling a unit cell was carried out by considering the individual phases as isotropic and only the matrix’s material properties to be time-dependent. This was achieved by expressing the stress–strain relation in terms of the bulk and shear modulus, where only the shear modulus was time-dependent, similar to the approach of [39]. This means, that we can assume that the viscous part of the deformation is incompressible and thus the volumetric change is purely elastic. Hence, the constitutive stress–strain relationship is written as:

(14)${\sigma }_{ij}=K{\epsilon }_{el,vol}I+2G{\epsilon }_{d,ij}$

Here, the strain in the first term of the right-hand side is the volumetric change, derived from the trace of the strain tensor as:

(15)${\epsilon }_{el,vol}=trace\left({\epsilon }_{el}\right)$

The second is the deviatoric strain, which is defined as:

(16)${\epsilon }_d={\epsilon }_{el}-\frac{1}{3}{\epsilon }_{el,vol}I$

If we consider the shear modulus to be time-dependent; it can be written in terms of the linear viscoelastic constitutive equation. The stress–strain relationship for an isotropic material can be described by:

(17)${\sigma }_{ij}=K{\epsilon }_{el,vol}I+2{{\displaystyle \int }}_0^t\mathrm{G}\left(t-t^{\prime }\right)\frac{\partial {\epsilon }_{d,ij}}{\partial t^{\prime }}d{t}^{\prime }$

Therefore, the shear modulus is derived from the measurement results by assuming the bulk modulus to be constant:

(18)$G\left(t\right)=\frac{3KE\left(t\right)}{9K-E\left(t\right)}$

Similar to the elements of the stiffness tensor in Equation (6), it is possible to describe the relaxation behavior of most polymer materials with the generalized Maxwell model and thus in terms of a Prony series as follows:

(19)$G\left(x_k,{\xi }_{n+1}-{\xi }^{\prime }\right)=G_{\infty }+{\displaystyle \sum }_{m=1}^{\mathrm{M}}{G}_m{\mathrm{e}}^{-\left({\xi }_{n+1}-{\xi }^{\prime }\right)/{\rho }_m}$

Inserting the constitutive equation of the isotropic material in Equation (17) and the Prony series in Equation (19) in the procedure of Zocher [11], the general incremental form (Equations (8)–(11)) can then written in simplified form as:

(20)$\mathsf{\Delta }{\sigma }_{\alpha }=\mathrm{K}\mathsf{\Delta }{\epsilon }_{el,vol}{I}_{\alpha }^V+2G^{\prime }{R}_{\alpha }^{-1}\mathsf{\Delta }{\epsilon }_{d,\alpha }+\mathsf{\Delta }{\sigma }_{\alpha }^{\mathrm{R}}$

where $I_{\alpha }^V$ denotes the identity tensor in Voigt notation (Appendix A). Similarly, R is a transformation vector from index to Voigt notation. In Equation (2D), it is also known as the Reuter matrix (Appendix A). The variable G′ is defined as:

(21)$G^{\prime }\equiv G_{\infty }+\frac{1}{\Delta \xi }{\displaystyle \sum }_{m=1}^{\mathrm{M}}{\eta }_m\left(1-{\mathrm{e}}^{-\Delta \xi /{\rho }_m}\right)$

This term can be linearized while maintaining the conservative character of this approach. For small time steps applying the rule of de L’Hospital, it takes the form of the shear modulus:

(22)$\underset{\Delta \xi \to 0}{\lim }{G}^{\prime }=G_{\infty }+{\displaystyle \sum }_{m=1}^{\mathrm{M}}\frac{{\eta }_m}{{\rho }_m}=G_0$

By replacing the stiffness tensor C_ijkl with the time-dependent scalar 2G, Equations (10) and (11) can be written as:

(23)$\Delta {\sigma }_{\alpha }^{\mathrm{R}}=-{\displaystyle \sum }_{m=1}^{\mathrm{M}}\left(1-{\mathrm{e}}^{-\Delta \xi /{\rho }_m}\right)S_{{\alpha }_m}\left({\xi }_n\right)$

(24)$S_{{\alpha }_m}\left(x_k,{\xi }_n\right)=e^{-\Delta \xi /{\rho }_m}{S}_m\left(x_k,{\xi }_n-\Delta \xi \right)+{\eta }_{m}2R_{\alpha }^{-1}{R}_{\epsilon ,\alpha }\left(1-e^{-\Delta \xi /{\rho }_m}\right)$

2.3. Adaptation to Homogenized Composite Behavior

A representative volume element (RVE) was set up, which consisted of two isotropic materials: a linear elastic fiber and a linear viscoelastic matrix. The fiber’s properties were chosen the same as for physical tests in Section 2, while the matrix was modeled by the method described in the previous section and the measured parameter values shown in Section 3. The compression modulus of the neat resin was chosen to be 2.78 GPa, derived from the elastic modulus and Poisson’s ratio in the elastic state.

The geometries of different textile structures and fiber volume fractions can be created efficiently by using the available material modelers of commercial FEM software. The incremental formulation of linear viscoelastic materials was implemented in a user material subroutine of Abaqus 6.13. It aimed to find a description for each element of the orthotropic relaxation tensor:

(25)$\left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\sigma }_4\\ {\sigma }_5\\ {\sigma }_6\end{array}\right]=\left[\begin{array}{cc}\begin{array}{ccc}{\mathsf{\Psi }}_{11}& {\mathsf{\Psi }}_{21}& {\mathsf{\Psi }}_{21}\\ {\mathsf{\Psi }}_{12}& {\mathsf{\Psi }}_{22}& {\mathsf{\Psi }}_{23}\\ {\mathsf{\Psi }}_{12}& {\mathsf{\Psi }}_{23}& {\mathsf{\Psi }}_{22}\end{array}& \\ 0& \begin{array}{ccc}{\mathsf{\Psi }}_{44}& & \\ & {\mathsf{\Psi }}_{44}& \\ & & {\mathsf{\Psi }}_{66}\end{array}\end{array}\right]{\left[\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ {\epsilon }_3\\ {\epsilon }_4\\ {\epsilon }_5\\ {\epsilon }_6\end{array}\right]}_0$

To do this, virtual load cases were set up, that created a unidirectional strain state. It consisted of 3024 elements of type C3D15, which is a 15-node quadratic triangular prism element. The fiber-matrix interface was modeled as perfectly bonded. In Figure 3, one can see the direction of the applied loads and the qualitative behavior of the set-up micromechanical models. To measure the response one could either evaluate the volume averaged stress or the reaction force. Here, as an indicator the reaction force was chosen (R_f). The surface on which it was measured is also indicated in Figure 3. The boundary conditions are listed in Table A1 in Appendix B.

The resulting normalized values of the relaxation tensor at a temperature of 110 °C are shown in Figure 4. They can be compared with the normalized relaxation modulus of the neat resin, which is indicated by a dotted line. As shown by Meder et al. [35], some components appear to have an increasing value over time.

All these functions depend on one single time-dependent variable G(t). The decrease of each value has the same decay time, so all the curves can be expressed in terms of a scalar form factor and a master curve. Here, the form factor is defined as the ratio of the amount of normalized relaxation of the actual index and the neat resin reference at a reference time (Equation (26)). The respective elements of the relaxation tensor can then be described in terms of the elastic index, the form factor and the normalized relaxation modulus of the neat resin (Equation (27)):

(26)$f_{\alpha \beta }\left(t_{ref}\right)=\frac{1-{\overline{\mathsf{\Psi }}}_{\alpha \beta }\left(t_{ref}\right)}{1-\overline{G}\left(t_{ref}\right)}$

(27)${\mathsf{\Psi }}_{\alpha \beta }\left(t\right)=C{\left(0\right)}_{\alpha \beta }{\overline{\mathsf{\Psi }}}_{\alpha \beta }\left(t\right)=C_{\alpha \beta }\left(0\right)\left(\left(1-f_{\alpha \beta }\right)+f_{\alpha \beta }\overline{G}\left(t\right)\right)$

Inserting this approach in Equation (6) (for the case of a relaxation test $C{\left(t\right)}_{\alpha \beta }={\mathsf{\Psi }}_{\alpha \beta }\left(t\right)$), the Prony series description of the stiffness tensor indices becomes:

(28)$C_{\alpha \beta }\left({\xi }_{n+1}-{\xi }^{\prime }\right)=C_{\alpha \beta }\left(0\right)\left[\left(1-f_{\alpha \beta }\right)+f_{\alpha \beta }{G}_{\infty }+{\displaystyle \sum }_{m=1}^{\mathrm{M}}{f}_{\alpha \beta }{\overline{G}}_m{\mathrm{e}}^{-\left({\xi }_{n+1}-{\xi }^{\prime }\right)/{\rho }_m}\right]$

As in the isotropic case, this description of the stiffness tensor is inserted in the approach of Zocher [11] in Voigt notation. The respective form of Equation (A1) becomes:

(29)$\mathsf{\Delta }{\sigma }_{\alpha }=C_{\alpha \beta }^{\prime }\Delta {\epsilon }_{\beta }+\Delta {\sigma }_{\alpha }^R$

Similar to the isotropic case, the first parameter on the right-hand side was approximated for small time steps without losing conservativity:

(30)$\underset{\Delta \xi \to 0}{\lim }{C}_{\alpha \beta }^{\prime }=C_{\alpha \beta }\left(0\right)$

Since the decay time $\rho$ does not depend on β, the sums in the adapted form of (1C) can be interchanged:

(31)$\Delta {\sigma }_{\alpha }^{\mathrm{R}}=-{\displaystyle \sum }_{\beta =1}^6{\displaystyle \sum }_{m=1}^M\left(1-{\mathrm{e}}^{-\frac{\Delta \xi }{{\rho }_m}}\right)S_{\alpha {\beta }_m}\left({\xi }_n\right) =-{\displaystyle \sum }_{m=1}^{\mathrm{M}}\left(1-{\mathrm{e}}^{-\Delta \xi /{\rho }_m}\right){\displaystyle \sum }_{\beta =1}^6{S}_{\alpha {\beta }_m}\left({\xi }_n\right)$

By comparing Equations (6) and (28), one can notice that only the part in the sum operator is time-dependent. It can be understood that Equation (11) changes to:

(32)$S_{\alpha {\beta }_m}\left({\xi }_n\right)={\mathrm{e}}^{-\Delta \xi /{\rho }_m}{S}_{\alpha {\beta }_m}\left({\xi }_n-\Delta \xi \right)+C_{\alpha \beta }\left(0\right)f_{\alpha \beta }{\overline{\eta }}_{\mathrm{m}}{R}_{\epsilon }\left(1-{\mathrm{e}}^{-\Delta \xi /{\rho }_m}\right)$

With the aim of reducing the number of state variables, an auxiliary state variable $S^{\prime }$ is defined. It consists of the sum over $\beta$ of the initial state variables:

(33)$S_{{\alpha }_m}^{\text{'}}≔{\displaystyle \sum }_{\beta =1}^6{S}_{\alpha {\beta }_m}\left({\xi }_n\right)$

Equation (31) then takes the following form:

(34)$\Delta {\sigma }_{\alpha }^{\mathrm{R}}=-{\displaystyle \sum }_{m=1}^{\mathrm{M}}\left(1-{\mathrm{e}}^{-\Delta \xi /{\rho }_m}\right)S_{{\alpha }_m}^{\text{'}}$

By writing Equation (32) in terms of $S^{\prime }$, we obtain:

(35)$S_{{\alpha }_m}^{\text{'}}\left({\xi }_n\right)={\mathrm{e}}^{-\Delta \xi /{\rho }_m}{S}_{{\alpha }_m}^{\text{'}}\left({\xi }_n-\Delta \xi \right)+{\overline{\eta }}_m\left(1-{\mathrm{e}}^{-\Delta \xi /{\rho }_m}\right){\displaystyle \sum }_{\beta =1}^6{C}_{\alpha \beta }\left(0\right)f_{\alpha \beta }{R}_{\epsilon ,\beta }$

For small time steps, this can be approximated by:

(36)$S_{{\alpha }_m}^{\text{'}}\left({\xi }_n\right)={\mathrm{e}}^{-\Delta \xi /{\rho }_m}{S}_{{\alpha }_m}^{\text{'}}\left({\xi }_n-\Delta \xi \right)+\frac{{\overline{\eta }}_m}{{\overline{\rho }}_m}{\displaystyle \sum }_{\beta =1}^6{C}_{\alpha \beta }\left(0\right)f_{\alpha \beta }{R}_{\epsilon ,\beta }\Delta \xi$

We can then insert the relationship shown in Equation (12):

(37)$S_{{\alpha }_m}^{\text{'}}\left({\xi }_n\right)={\mathrm{e}}^{-\Delta \xi /{\rho }_m}{S}_{{\alpha }_m}^{\text{'}}\left({\xi }_n-\Delta \xi \right)+\frac{{\overline{\eta }}_m}{{\overline{\rho }}_m}{\displaystyle \sum }_{\beta =1}^6{C}_{\alpha \beta }\left(0\right)f_{\alpha \beta }\Delta {\epsilon }_{\beta }$

For the expression of the spring element of the normalized shear modulus, it is helpful to insert the parameter introduced in Equation (13). It represents the normalized amount of the decay described by one of the branches of the generalized Maxwell element, which was measured for the neat resin system.

(38)$S_{{\alpha }_m}^{\text{'}}\left({\xi }_n\right)={\mathrm{e}}^{-\Delta \xi /{\rho }_m}{S}_{{\alpha }_m}^{\text{'}}\left({\xi }_n-\Delta \xi \right)+{\phi }_m{\displaystyle \sum }_{\beta =1}^6{C}_{\alpha \beta }\left(0\right)f_{\alpha \beta }\Delta {\epsilon }_{\beta }$

This method of describing the relaxation tensor indices with form factors allows reduction of the number of state variables of M × 6.

The concept presented here only works if there is no more than one individual material of the composite modeled to show the viscoelastic behavior. If there is more than one viscoelastic material included, the components of the relaxation tensor will have different decay times and cannot be transformed to a master curve with a form factor.

2.4. Generation of Experimental Data

The resin system investigated in this study was provided by Huntsman (Basel, Switzerland) and has the commercial name XB3585. It is a fast-curing epoxy resin for automotive applications that can achieve short curing-cycle times. The system has a final glass transition temperature of 120–130 °C and a flexural modulus of 2.8–3.1 GPa.

The glass fiber reinforcement textile was provided by Kümpers Composites (Salzbergen, Germany) as a unidirectional non-crimping fabric. The number of layers was chosen for a targeted fiber volume content of 60%. A summary of the number of layers, their orientation and the dimensions of the samples are presented in Table 1. Because of the limited maximum force of the test device (15 N), the specimens for the measurement in longitudinal dimensions were manufactured with a lower thickness. The textile consisted of S-glass fibers provided by 3B (Hoeilaart, Belgium). They are listed as SE 2020 2400 dTex and are expected to provide a tensile modulus of 81 GPa.

The relaxation of the specimens was measured with a 3-point bending test in a dynamic mechanical analyzer (DMA, Q800) from TA Instruments (New Castle, DE, USA) at different temperatures. The fully cured specimens were tested at isothermal temperatures ranging from 125 °C to 50 °C in steps of −5 °C. For each temperature step, an initial displacement of 0.1% strain was applied and the stress relaxation was recorded for 60 min. The relaxation modulus of the measured data according to the software’s internal manual was calculated as:

(39)$E\left(t\right)=K_S\left(t\right)·\frac{L^3}{6I}\left[1+\frac{6}{10}\ast \left(1+v_P\right){\left(\frac{z}{L}\right)}^2\right]$

where K_s is the measured relationship between strain and force, L is the support width, z is the sample thickness, I is the moment of inertia and v_P is Poisson’s ratio. The support width was chosen to be 50 mm.

2.5. Validation and Application on the Meso-Level/Ply Level

In order to validate the developed model, a virtual copy of the three-point bending test was built, whose simulation results were compared with the experimental data from the DMA. It consists of 896 elements of type C3D20R which is a 20-node quadratic brick. As for the micromechanical models, the incremental formulations of the transversely isotropic linear viscoelastic material behavior were implemented in a user subroutine. In Figure 5, the meshing and typical stress distribution during load application are shown. The derived relaxation stiffness of the virtual test was calculated as for the physical test using Equation (39) and the displacement and total reaction force. The supports were modeled by constraining the degree of freedom in a vertical direction. This is sufficient for the consideration here since the modeling of nonlinear effects of the supports is not intended.

3. Results

This section is structured to provide the best possible overview of the workflow carried out: First, the experimentally generated data of the three-point bending test are presented, divided into the measurement of the elastic modulus and the relaxation modulus. The data from the pure resin sample are used to calibrate the RVE model. Because this paper focuses on viscoelastic homogenization, the measured elastic modulus of the fiber-reinforced samples is used for the calibration of the validation model, while the measured relaxation modulus is used as the validation. In Section 3.3 the results of the RVE modeling are shown and the form factors are derived. This is followed by the validation in Section 3.4.

3.1. Experimentally-Derived Elastic Modulus

In Figure 6 the storage and loss modulus of the dynamic DMA measurements of the neat resin sample are presented. One can see that after an approximately linear decrease at lower temperatures, the storage modulus drops dramatically, starting at about 110 °C.

The corresponding result for the composite samples can be seen in Figure 7. While the storage modulus in the longitudinal direction remains stable until a temperature up to 115 °C, the modulus in transverse direction shows a similar behavior as the neat resin system. It is remarkable that the loss modulus of the material in the transverse direction shows an increase even at lower temperatures.

In the rheological model, which was used to describe the relaxation behavior of the neat resin, the elastic modulus is the sum of all springs of the generalized Maxwell model. In order to derive an analytical description, a customized expression was fitted to the measurement data. It consists of a linear part and a step function:

(40)$E^{\prime }=a+b·x+\frac{c}{1+exp\left(\frac{T-d}{e}\right)} \left[MPa\right]$

The parameters obtained after fitting this expression by linear regression to the measurement data are listed in Table 2.

The macroscopic model also requires a continuous description of the modulus of elasticity. For this purpose, the data of the composite samples was fitted in the same way. The resulting values are given in Table 2.

3.2. Experimentally-Derived Relaxation Modulus

The relaxation modulus of the neat resin system was measured at different constant temperatures from 50 to 125 °C. The isotropic viscoelastic description presented in Section 2.2 with two branches (M = 2) was fitted to the data using a linear regression algorithm. The measurement data and the result of the fitting are shown in Figure 8.

The fitting of the isotropic viscoelastic description resulted in a set of values for the parameters of the rheological model. They are shown in Figure 9. With increasing temperatures, the material approaches a rubbery and liquid state, and the decay time, indicated by the parameter ${\mathsf{\rho }}_i$, decreases. Similarly, the amount of total relaxation increases, visible in the increase of ratio of spring elements connected to a branch with a damper element ($k_1/k_{tot}$ resp. $k_2/k_{tot}$). k_tot denotes the sum of all springs in the rheological model, which was set to the same value as the modeled elastic modulus from the previous measurement.

To validate the macroscopic model later, the relaxation behavior of the composite material was measured. The result is shown in Figure 10. It can be seen that the relaxation modulus in the fiber direction shows hardly any relaxation up to about 100 °C. At higher temperatures, the relaxation modulus drops and reaches a level about 40% lower than the elastic value at room temperature. The relaxation modulus in the transverse direction shows a decrease at room temperature and continuously decreases at higher temperatures up to a reduction greater than 95% at 125 °C.

3.3. Results of RVE Homogenization

The normalized elements of the relaxation tensor ${\overline{\mathsf{\Psi }}}_{\alpha \beta }$ obtained from the micromechanical model were compared with the behavior of the neat resin with Equation (26). The resulting form factors at a reference temperature of 110 °C are shown in Table 3.

In order to demonstrate the validity of the assumption that all elements ${\overline{\mathsf{\Psi }}}_{\alpha \beta }$ can be expressed in terms of $\overline{G}$ and a scalar value $f_{\alpha \beta }$, the normalized relaxation tensor indices were transformed with Equation (41) and plotted together with the neat resin reference in Figure 11.

(41)${\overrightarrow{\mathsf{\Psi }}}_{\alpha \beta }\left(t\right):=\frac{{\overline{\mathsf{\Psi }}}_{\alpha \beta }\left(t\right)-1}{f_{\alpha \beta }}+1$

The largest relative deviation can be found for index ${\overrightarrow{\mathsf{\Psi }}}_{44}$ with 20% followed by index ${\overrightarrow{\mathsf{\Psi }}}_{66}$ with 10%. For all other indices, the maximum relative deviations are less than 2%.

3.4. Validation of Macroscopic Model

For an assessment of the quality of the relaxation model, the measured modulus and the modeled relaxation modulus derived from a three-point bending test were compared. The result is shown in Figure 12 and Figure 13 in normalized form. The comparison in normalized form was chosen in order to assess the quality of the relaxation model independent of the quality of the elastic modeling.

Data for temperatures above 115 °C were excluded, since the normalized measurement data do not show a trend towards increasing relaxation anymore. This can be explained by the measurement setup, for which the fixed sampling rate was not fast enough to capture the decay at the beginning of the measurement. This is indicated by the decrease in k₂/k_tot in Figure 9.

The modeling in the longitudinal direction is in a good agreement with the measurement data, while the modeling in the transverse direction can be seen to underestimate the amount of relaxation.

4. Discussion

Although the approach of describing the response of the individual elements of the stiffness tensor in terms of the behavior of the pure resin system and a scalar factor leads to an error of less than 2% for most elements, a larger deviation of up to 20% can be seen for the index ${\overrightarrow{\mathsf{\Psi }}}_{44}$ and ${\overrightarrow{\mathsf{\Psi }}}_{66}$. As a consequence, the model currently creates a maximum error for shear deformation of about the same order of magnitude. The deviation could be caused by the setup of the boundary conditions for this load case which is visible in the uneven surfaces of the 2–3 plane. This could be increased by the introduction of periodic boundary conditions as described by [25] in detail, where also a user-defined material subroutine in Abaqus was used.

The presented incremental method for the macroscopic time-dependent description of fiber composite materials represents an efficient alternative to established methods for FEM simulations. By introducing scalar form factors for the indices of the relaxation tensor, it allows the use of Mx6 state variables instead of M × 36 to describe an orthotropic material. This reduces the computational effort.

An increase of the accuracy of the model could be achieved by the introduction of temperature-dependent form factors. Furthermore, in order to derive a continuous description of the temperature dependency of the relaxation, the presented homogenization is able to introduce the time-temperature superposition principle with the use of the reduced time variable ξ. The time shift factors can be found by creating a master curve using the experimental data in Figure 8.

Linear viscoelastic material behavior is only an appropriate assumption for some polymers in the solid state at relatively low load levels. The higher the temperature, the more the state of the polymer changes from the glassy to the rubbery state. In measurements of the elastic modulus, an indicator of this change is the increase of the loss modulus. In the measurements taken here, this was the case for the neat resin at about 115 °C (Figure 6), for the composite in the longitudinal direction about the same (Figure 7), but in transverse direction it can be seen that there is an increase of the loss modulus already at lower temperatures which is an indicator for increased nonlinear behavior. This can be explained by the strain amplification effect on the microlevel, which reaches of the elastic limit of the material in some places at lower nominal load levels. One can assume that the strain amplification effect is more pronounced in the transverse direction, since there are stiffness jumps in the loading direction.

The strain amplification also changes the time-dependent behavior of the composite from linear viscoelastic to non-linear viscoplastic behavior. This was especially pronounced in the matrix-dominated directions (transverse), where the amount of relaxation in the measurement—especially at higher temperatures—exceeds the model’s prediction which assumes linear viscoelastic behavior. This means that the model provides a conservative assessment of the residual stresses, which is generally required for engineering purposes.

5. Conclusions

An incremental form to describe the linear viscoelastic response of isotropic and orthotropic materials was derived. Both forms were successfully implemented in an incremental FEM analysis.

With the use of a representative volume element it was shown that in the case of only one viscoelastic phase, the response of an element of the relaxation tensor can be expressed in terms of a scalar factor and the response of the viscoelastic phase. In this case, six form factors for a transversely isotropic material were derived. By inserting these factors into the incremental form for an orthotropic material, the macroscopic description of transversely isotropic material was converted to an extraordinarily easy form shown in Equations (29), (30), (34) and (38). Since the form factors describe the deviation of an index of the relaxation tensor of the composite in comparison with the neat resin, the impact of the virtual parameter derivation can be assessed easily, which is beneficial for engineering purposes. The result is a linear viscoelastic model with good handling and efficiency.

The model was validated through simulation of a three-point bending test and comparing the simulated results with the measurements in the longitudinal and transverse directions of a unidirectional fiber-reinforced material. A promising agreement was shown for up to 15 °C below the glass transition temperature. Since at these temperatures it was seen that most of the stress relaxes and the decay times are small, the present model was assessed to sufficiently describe the induced stresses of a typical manufacturing process of a thermoset material. Towards higher temperatures and, thus, increasing deviation of the model from the measured relaxation modulus, the model assessment still stays conservative due to the linear viscoelastic assumption, which is beneficial for component design.

Author Contributions

N.G.: Writing—original draft, conceptualization, methodology; I.Z.: Writing—review and editing; C.B.: Formal analysis, investigation, supervision, project administration. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

Special thanks goes to KATZ, Kunststoff Ausbildungs- und Technologie-Zentrum (Aarau, Switzerland) for the use of their facilities.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Nomenclature

$\mathsf{\alpha }$	Degree of cure
${\beta }_{ij}^{\prime }$	Second-order tensor of relaxation moduli relating stress to thermal strain
$C_{ijkl}$	Fourth-order tensor of orthotropic relaxation moduli relating stress to mechanical strain
$C_{ijk{l}_{\infty }}$	Response of stiffness tensor element after infinite time
$C_{ijk{l}_m}$	Value of mth Prony series element of stiffness tensor
$E^{\prime }$	Storage modulus
$E\left(t\right)$	Time dependent tensile modulus
$E_{\perp }$	Tensile modulus in transverse direction
$E_{\Vert }$	Tensile modulus in longitudinal direction
${\epsilon }_d$	Deviatoric strain
${\epsilon }_{kl}$	Cauchy strain
${\epsilon }_{vol}$	Volumetric strain
$f_{\alpha \beta }$	Form factor, the normalized amount of relaxation compared to the neat resin
$G\left(t\right)$	Time dependent shear modulus isotropic material
$I$	Identity tensor
$K$	Bulk modulus
$k$	Maxwell model spring Parameter constant
$K_S$	Measured force divided by measured displacement
M	Number of Branches of Maxwell Model
$n$	Number of time step
$\phi$	Relation of spring element of rheological model to the sum of all spring elements
${\mathsf{\Psi }}_{kl}$	Relaxation tensor element, special case of $C_{ijkl}$ for the case of a constant strain load
${\overline{\mathsf{\Psi }}}_{\alpha \beta },\overline{G}$	Normalized relaxation function (value at t = 0 is 1)
$R$	Reuter matrix, transformation from index to Voigt notation
$R_{\mathsf{\epsilon }}$	Strain increment per reduced time
$R_f$	Reaction force
${\rho }_{ijk{l}_m}$	Decay time of branch m, ratio of viscosity η and spring stiffness k
$S_{ijk{l}_m}$	State variable of branch m
${\sigma }_{ij}$	Cauchy stress tensor
$\Delta {\sigma }_{ij}^{\mathrm{R}}$	Stress increment considering the change of stress due to linear viscoelastic behavior
$T$	Temperature
$\mathsf{\Delta }\mathsf{\Theta }$	Temperature strain
$x_k$	Spatial position
${\xi }_{n+1}$	Reduced time

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Von Mises strain distribution for the application of only thermal strain in the micromodel [37] (the blue part is the fiber).

Figure 2. Generalized Maxwell model or Wiechert [38] model with M = 2 branches.

Figure 3. Definition of the boundary conditions for the derivation of the form factors defined in Equation (25).

Figure 4. Normalized values of the entries of the relaxation tensor at 110 °C.

Figure 5. Geometry and qualitative stress distribution of the implemented validation model. Virtual copy of mechanical test setup. u denotes the applied deflection and Rf the reaction force in the middle of the specimen.

Figure 6. Measurement and fit of the storage and loss modulus of neat resin.

Figure 7. Measurement and fit of the storage and loss modulus of the UD-composite: (a) longitudinal; (b) transverse.

Figure 8. Comparison of the relaxation data and the relaxation model described by the continuous temperature-dependent parameters.

Figure 9. Resulting values of the element parameters of the rheological model after fitting to measurement data. It shows the decay time [Forumla omitted. See PDF.] (a): branch 1, (c): branch 2 and the spring constant [Forumla omitted. See PDF.] compared to the sum of all spring constants [Forumla omitted. See PDF.] (b): branch 1, (d): branch 2.

Figure 10. Measured relaxation modulus: (a) longitudinal direction; (b) transverse direction.

Figure 11. Relaxation tensor components transferred with the derived form factor. The transferred functions are in a good agreement with the reference function of the neat resin reference at 110 °C.

Figure 12. Comparison of the modeled and measured normalized relaxation modulus of the composite in the longitudinal direction.

Figure 13. Comparison of the modeled and measured normalized relaxation modulus of the composite in the transverse direction.

Table 1

Information of the samples.

Type	Thickness (mm)	Number of Layers	Width (mm)	Length (mm)
Neat resin	4 ± 0.06	-	15	65
Longitudinal	1 ± 0.04	2
Transverse	4 ± 0.06	8

Table 2

Fitting parameters of the composite’s elastic modulus.

Direction	a	b	c	d	e
Neat	2.846 × 102	−7.105	2.936 × 103	1.301 × 102	6.385
Longitudinal	3.176 × 104	−6.417	9.770 × 103	1.252 × 102	4.027
Transverse	2.595 × 103	−45.14	1.548 × 104	1.217 × 102	8.583

Table 3

Results for the form factor f (110 °C).

Index	${\mathit{f}}_{\mathit{\alpha }\mathit{\beta }}\left({\mathit{t}}_{\mathit{r}\mathit{e}\mathit{f}}\right)$
${\mathsf{\Psi }}_{11}$	2.307 × 10−2
${\mathsf{\Psi }}_{22}$	3.592 × 10−1
${\mathsf{\Psi }}_{44}$	9.161 × 10−1
${\mathsf{\Psi }}_{66}$	6.250 × 10−1
${\mathsf{\Psi }}_{12}$	6.453 × 10−3
${\mathsf{\Psi }}_{23}$	1.342 × 10−1

Appendix A

R is called the Reuter matrix:(A1)$$R=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ & & & 2& 0& 0\\ & sym& & & 2& 0\\ & & & & & 2\end{array}\right]$$

Furthermore, the identity matrix in Voigt notation becomes:(A2)$$I_{\alpha }^V=\left[\begin{array}{c}1\\ 1\\ 1\\ 0\\ 0\\ 0\end{array}\right]$$

Appendix B

Table A1

Boundary conditions.

	Ψ11	Ψ22	Ψ12	Ψ23	Ψ44	Ψ55
Xfront		Ux = 0		Ux = 0		Ux = 0
Xback	Ux = 0	Ux = 0	Ux = 0	Ux = 0		Ux = 0
Yfront	Uy = 0		Uy = 0		Uy, Uz = 0	Ux,uy = 0
Yback	Uy = 0	Uy = 0	Uy = 0	Uy = 0	Encastre	Encastre
Zfront	Uz = 0	Uz = 0	Uz = 0	Uz = 0	Uz = 0
Zback	Uz = 0	Uz = 0	Uz = 0	Uz = 0	Uz = 0
Displacement Plane	Xfront	Yfront	Xfront	Yfront	Yfront	Yfront
Loading Direction	X	Y	X	Y	X	Y
Reaction Plane	Xfront	Yfront	Yfront	Zfront	Yfront	Yfront

References

1. Brauner, C. Analysis of Provess Induced Distortions and Residual Stresses in Composite Structures; Logos Verlag: Berlin, Germany, 2013.

2. Hein, R.; Wille, T.; Gabtni, K.; Dias, J.P. Prediction of Process-Induce Distortions and Residual Stresses of a Composite Suspension Blade. Defect and Diffusion Forum; Trans Tech. Publications Ltd.: Bäch, Switzerland, 2015; Volume 362, pp. 224-243. [DOI: https://dx.doi.org/10.4028/www.scient]

3. Hahn, H.T.; Pagano, N.J. Curing Stresses in Composite Laminates. J. Compos. Mater.; 1975; 9, pp. 91-106. [DOI: https://dx.doi.org/10.1177/002199837500900110]

4. Johnston, A. An Integrated Model of the Development of Process-Induced Deformation in Autoclave Processing of Composite Structures. Ph.D. Thesis; University of British Columbia: Vancouver, BC, Canada, May 1997.

5. Bogetti, T.A.; Gillespie, J.W., Jr. Process-Induced Stress and Deformation in Thick-Section Thermoset Composite Laminates. J. Compos. Mater.; 1992; 26, pp. 626-660. [DOI: https://dx.doi.org/10.1177/002199839202600502]

6. White, S.R.; Hahn, H.T. Process Modeling of Composite Materials: Residual Stress Development during Cure. Part I. Model Formulation. J. Compos. Mater.; 1992; 26, pp. 2402-2422. [DOI: https://dx.doi.org/10.1177/002199839202601604]

7. Weitsman, Y. Residual Thermal Stresses Due to Cool-Down of Epoxy-Resin Composites. J. Appl. Mech.; 1979; 46, pp. 563-567. [DOI: https://dx.doi.org/10.1115/1.3424606]

8. Wijskamp, S. Shape Distorsions in Composite Forming. Ph.D. Thesis; University Twente: Twente, The Netherlands, June 2004.

9. White, S.R.; Kim, Y.K. Process-Induced Residual Stress Analysis of As4/3501-6 Composite Material. Mech. Compos. Mater. Struct.; 1998; 5, pp. 153-186. [DOI: https://dx.doi.org/10.1080/10759419808945897]

10. Prasatya, P.; McKenna, G.B.; Simon, S.L. A Viscoelastic Model for Predicting Isotropic Residual Stresses in Thermosetting Materials: Effects of Processing Parameters. J. Compos. Mater.; 2001; 35, pp. 826-848. [DOI: https://dx.doi.org/10.1177/a037322]

11. Zocher, M.; Groves, S.; Allen, D. A three-dimensional finite element formulation for thermoviscoelastic orthotropic media. Int. J. Numer. Methods Eng.; 1997; 40, pp. 2267-2288. [DOI: https://dx.doi.org/10.1002/(SICI)1097-0207(19970630)40:12<2267::AID-NME156>3.0.CO;2-P]

12. Svanberg, M.; Holmberg, A. Prediction of shape distortions Part I. FE-implementation of a path dependent constitutive model. Compos. Part A Appl. Sci. Manuf.; 2004; 35, pp. 711-721. [DOI: https://dx.doi.org/10.1016/j.compositesa.2004.02.005]

13. Halpin, J.C. Effects of Environmental Factors on Composite Materials; Air Force Materials Lab Wright-Patterson: Montgomery County, OH, USA, 1969.

14. Halpin, J.; Kardos, J. The Halpin-Tsai Equations: A Review. Polym. Eng. Sci.; 1976; 16, pp. 344-352.

15. Whitney, J.M. Shear Correction Factors for Orthotropic Laminates under Static Load. J. Appl. Mech.; 1973; 40, pp. 302-304. [DOI: https://dx.doi.org/10.1115/1.3422950]

16. Chamis, C.C. Simplified composite micromechanics equations for hygral, thermal and mechanical properties. Proceedings of the Annual Conference of the Society of the Plastics Industry (SPI) Reinforced Plastics/Composites Institute; Houston, TX, USA, 1 January 1983.

17. David Müzel, S.; Bonhin, E.P.; Guimarães, N.M.; Guidi, E.S. Application of the Finite Element Method in the Analysis of Composite Materials: A Review. Polymers; 2020; 12, 818. [DOI: https://dx.doi.org/10.3390/polym12040818] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32260389]

18. Dandekar, C.R.; Shin, Y.C. Modeling of machining of composite materials: A review. Int. J. Mach. Tools Manuf.; 2012; 57, pp. 102-121. [DOI: https://dx.doi.org/10.1016/j.ijmachtools.2012.01.006]

19. Arola, D.; Ramulu, M. Orthogonal cutting of fiber-reinforced composites: A finite element analysis. Int. J. Mech. Sci.; 1997; 39, pp. 597-613. [DOI: https://dx.doi.org/10.1016/S0020-7403(96)00061-6]

20. Barbero, E.J. Finite Element Analysis of Composite Materials Using Ansys; 2nd ed. CRC Press: Boca Raton, FL, USA, 2014.

21. Mkaddem, A.; Demirci, I.; Mansori, M. El A micro–macro combined approach using FEM for modelling of machining of FRP composites: Cutting forces analysis. Compos. Sci. Technol.; 2008; 68, pp. 3123-3127. [DOI: https://dx.doi.org/10.1016/j.compscitech.2008.07.009]

22. Yang, L.; Yan, Y.; Liu, Y.; Ran, Z. Microscopic failure mechanisms of fiber-reinforced polymer composites under transverse tension and compression. Compos. Sci. Technol.; 2012; 72, pp. 1818-1825. [DOI: https://dx.doi.org/10.1016/j.compscitech.2012.08.001]

23. Venu Gopala Rao, G.; Mahajan, P.; Bhatnagar, N. Machining of UD-GFRP composites chip formation mechanism. Compos. Sci. Technol.; 2007; 67, pp. 2271-2281. [DOI: https://dx.doi.org/10.1016/j.compscitech.2007.01.025]

24. Tenek, L.T.; Argyris, J. Finite Element Analysis for Composite Structures. Solid Mechanics and Its Applications; Gladwell, G.M. Springer: Dordrecht, The Netherlands, 1998; Volume 59.

25. Wang, Z.; Smith, D.E. Numerical analysis on viscoelastic creep responses of aligned short fiber reinforced composites. Compos. Struct.; 2019; 229, 111394. [DOI: https://dx.doi.org/10.1016/j.compstruct.2019.111394]

26. Jones, R.M. Mechanics of Composite Materials; 2nd ed. CRC Press: Boca Raton, FL, USA, 1998.

27. Hinton, M.J.; Kaddour, A.S.; Soden, P.D. Failure Criteria in Fibre-Reinforced-Polymer Composites; Elsevier: Amsterdam, The Netherlands, 2004.

28. Stier, B.; Simon, J.-W.; Reese, S. Numerical and experimental investigation of the structural behavior of a carbon fiber reinforced ankle-foot orthosis. Med. Eng. Phys.; 2015; 37, pp. 505-511. [DOI: https://dx.doi.org/10.1016/j.medengphy.2015.02.002] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25765189]

29. Tian, W.; Qi, L.; Zhou, J.; Liang, J.; Ma, Y. Representative volume element for composites rein-forced by spatially randomly distributed discontinuous fibers and its applications. Compos. Struct.; 2015; 131, pp. 366-373. [DOI: https://dx.doi.org/10.1016/j.compstruct.2015.05.014]

30. Lu, Z.; Yuan, Z.; Liu, Q. 3D numerical simulation for the elastic properties of random fiber composites with a wide range of fiber aspect ratios. Comput. Mater. Sci.; 2014; 90, pp. 123-129. [DOI: https://dx.doi.org/10.1016/j.commatsci.2014.04.007]

31. Schapery, R.A. Stress Analysis of Viscoelastic Composite Materials. J. Compos. Mater.; 1967; 1, pp. 228-267. [DOI: https://dx.doi.org/10.1177/002199836700100302]

32. Barbero, E.J.; Luciano, R. Micromechanical formulas for the relaxation tensor of linear viscoelastic composites with transversely isotropic fibers. Int. J. Solids Struct.; 1995; 32, pp. 1859-1872. [DOI: https://dx.doi.org/10.1016/0020-7683(94)00233-M]

33. Otero, J.; Díaz, R.G.; González, O.L.C.; Sabina, F.; Berger, H.; Böhlke, T. Asymptotic and numerical homogenization methods applied to fibrous viscoelastic composites using Prony’s series. Acta Mech.; 2020; 231, pp. 2761-2771. [DOI: https://dx.doi.org/10.1007/s00707-020-02671-1]

34. Rodríguez-Ramos, R.; Otero, J.A.; Cruz-González, O.L.; Guinovart-Díaz, R.; Bravo-Castillero, J.; Sabina, F.J.; Padilla, P.; Lebon, F.; Sevostianov, I. Computation of the relaxation effective moduli for fibrous viscoelastic composites using the asymptotic homogenization method. Int. J. Solids Struct.; 2020; 190, pp. 281-290. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2019.11.014]

35. Meder, G. Zur exakten und näherungsweisen Berechnung unidirektionalverstärkter Kunststoffe. Mater. Werkst.; 1981; 12, pp. 366-374. [DOI: https://dx.doi.org/10.1002/mawe.19810121008]

36. Chen, Y.; Zhao, Z.; Li, D.; Guo, Z.; Dong, L. Constitutive modeling for linear viscoelastic fiber-reinforced composites. Compos. Struct.; 2021; 263, 113679. [DOI: https://dx.doi.org/10.1016/j.compstruct.2021.113679]

37. Brauner, C.; Block, T.B.; Purol, H.; Herrmann, A.S. Microlevel manufacturing process simulation of carbon fiber/epoxy composites to analyze the effect of chemical and thermal induced residual stresses. J. Compos. Mater.; 2012; 46, pp. 2123-2143. [DOI: https://dx.doi.org/10.1177/0021998311430157]

38. Wiechert, E. Ueber Elastische Nachwirkung. Ph.D. Thesis; Königsberg University: Königsberg, Germany, 1889.

39. COMSOL. Documentation; Linear Viscoelastic Materials. 2022; Available online: https://doc.comsol.com/5.5/doc/com.comsol.help.sme/sme_ug_theory.06.26.html (accessed on 1 August 2022).