<figure><img src="image_1.png" /></figure>

HAL
open science

# A multi-scale analysis of local stresses during the cure of a composite tooling material

Emmanuel Lacoste, K. Szymanska, S. Terekhina, Sylvain Fréour, Frédéric Jacquemin, M. Salvia

► To cite this version:

Emmanuel Lacoste, K. Szymanska, S. Terekhina, Sylvain Fréour, Frédéric Jacquemin, et al.. A multi-scale analysis of local stresses during the cure of a composite tooling material. International Journal of Material Forming, 2013, 6 (4), pp.467-482. 10.1007/s12289-012-1100-0 . hal-01004853

HAL Id: hal-01004853

https://hal.science/hal-01004853

Submitted on 12 Feb 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A multi-scale analysis of local stresses development during the cure of a composite tooling material

E. Lacoste · K. Szymanska · S. Terekhina · S. Fréour · F. Jacquemin · M. Salvia

E. Lacoste · S. Fréour · F. Jacquemin
Institut de Recherche en Génie Civil et Mécanique (UMR CNRS
6183), LUNAM Université - Université de Nantes - Ecole Centrale
de Nantes,
37 Bd de l'Université, BP 406, 44602 Saint-Nazaire, France

K. Szymanska · S. Terekhina · M. Salvia
Laboratoire de Tribologie et Dynamique des Systèmes (UMR
CNRS 5513), Ecole Centrale de Lyon,
36 Avenue Guy de Collongue,
69134 Ecully Cedex, France

**Abstract** This paper is dedicated to the cure of an in-plane isotropic carbon-epoxy tooling material presenting a specific mesostructure. Eshelby-Kröner self-consistent model (EKSC) is used to achieve a two-steps scale transition procedure, allowing relating microscopic to macroscopic properties of the material, and estimating its multi-scale mechanical states. This procedure is used to predict the local residual stresses due to thermal and chemical shrinkage of the resin, depending on the manufacturing process conditions. An experimental investigation provides the BMI resin cure kinetics and mechanical properties as a function of the temperature and conversion degree. The consequences of these evolutions on the local mechanical states are investigated and discussed.

**Keywords** Composite material · Residual stress · Multiscale · Cure process · Thermoset polymer

## Introduction

General purpose

In the last decade, composite materials based on carbon fiber and thermoset matrices have been more and more

involved in the design of mechanical parts, particularly for aeronautical applications. Their high strength-to-weight ratio as well as corrosion and fatigue resistance make them an interesting alternative to metals, which are still commonly used. However, as a counterpart of this technological breakthrough, the engineering of composite parts yields several scientific challenges. Among them is the computation of internal stress at the microscopic level, due to the properties mismatch between the fibers and the matrix. As curing is considered, the strong chemical and thermal shrinkage of the matrix is counteracted by the fibers, which are far stiffer and less dilative, yielding self-compensated stresses ("residual stresses") in the constituents. These stresses only exist at a microscopic level; nevertheless, several authors showed that they significantly affect the macroscopic stress-strain and failure behaviour. They can result in fiber-matrix debonding, or in increased fibers waviness (microbuckling), both of these yielding a change of properties for the ply. A good review of this subject can be found in a series of bibliographical papers from Parlevliet et al. [35-37]. The present study is focused on the computation of these local stresses created during cure, when the material undergoes homogeneous loads.

## Review of predictive models

Various experimental approaches have been developed in order to measure the local stresses and strains in composite materials. An extensive description of these methods is given by Parlevliet et al. [36]. On the other hand, scale transition models proved to be an interesting approach for predicting both the local stresses, and the effective mechanical properties of composite plies.

A first class of models uses a statistical description of the microstructure, accounting for the mechanical properties, morphologies and volumic ratio of the phases constituting
the material, while neglecting the geometrical arrangement of the individual inclusions. The first historical examples of this approach are the well-known Tsaï-Hahn's rules of mixtures [46], which gave an easy-understandable estimate of the properties of unidirectional composite plies. Hashin and Shtrikman [17] also suggested rigorous bounds for composite material properties, which are still often used for achieving such calculations. An extension of this approach, the class of so-called "mean-field" models, became trendier in the last decade. Among them, Mori-Tanaka (MT) and Eshelby-Kröner Self-Consistent (EKSC) estimates [3, 23] are the most used. In exchange for their more pronounced mathematical complexity, these models can take into account the morphologies of the heterogeneities (inclusions) constituting the material (the interested reader can refer, for instance, to Jacquemin et al. [20]). Moreover, they can also be used to predict in a systematic and robust way, the mechanical states (stresses and strains) experienced by the very constituents of the materials. A comparison between these two models, focused on the application to the modeling of composite materials, can be found in [12]. In particular, predicting the residual micro-stresses during the cure process, which is the aim of the present study, can be achieved through both these approaches. However, since unidirectional structures are the most frequently considered, and relatively well represented by rules of mixtures, self-consistent methods are rarely used.

Another class of models takes advantage of the finite elements method (FEM) for a direct representation of the microstructure's geometrical complexity. This "full-field" approach provides a more detailed estimate of the local mechanical states. The implementation and computation costs required by FEM is nevertheless much more important than their mean-field counterparts. However, such an approach is not the topic of the present study, but the interested reader can, for instance, refer to [15, 16, 21].

## Purpose of the study

A semi-analytic self-consistent model will be used to carry out a multi-scale analysis of the residual cure stress in a specific composite material, which mechanical response cannot be appropriately described by classical rules of mixtures, because of a very specific microstructure. An experimental investigation enabled to characterise the mechanical behaviour and the cure kinetics of the resin constituting this material. These properties will be used as input data in EKSC model, in order to compute the composite's effective properties.

Residual stress at the microscopic level will also be computed in the case of a simplified cure process, illustrating the most interesting feature of the model. Three typical boundary conditions, approximately representing mould-

part interactions, will be considered here: (i) a stress-free cure of the material, (ii) a cure with fixed in-plane displacements, and (iii) a cure with prescribed thermal in-plane displacements. These conditions would approximately match those applied to a thin part cured into different kinds of moulds, respectively made of (i) polymer, (ii) Invar, or (iii) metal (steel). The stresses, experienced by the material, during the cure process, and after removal from the mould, will be computed at macroscopic and microscopic scales, through the multi-scale analysis.

## The scale transition procedure

### Description of the EKSC model

Eshelby-Kröner scale transition model is based on a representation of the material at two distinct scales: on the one hand, the «local» scale, denoted by the superscript ¹, where one observes the behaviour of each constituent, considered as an ellipsoidal and homogeneous inclusion with principal axes lengths {2 a₁, 2 a₂, 2 a₃} (also called Base Volume or BV); on the other hand, the macroscopic scale denoted by the superscript ¹, where one can observe the behaviour of the Homogeneous Equivalent Medium (or HEM). In the case of a chemical-thermo-chemio-elastic linear behaviour, the following law is written:

$$ \sigma^k = L^k : (\varepsilon^k - \alpha^k \cdot \Delta T - \eta^k \cdot \Delta \chi), \quad \text{where } k = i, I \quad (1) $$

In this relation, the stiffnesses are represented by the 4th-order tensors L, the Coefficients of Thermal Expansion (CTE) and Chemical Expansion (CCE) being referred to, by the 2nd-order tensors α and η. The temperature and cure degree increments are denoted by ΔT and Δχ respectively, whereas σ and ε stand for the stress and strain. Relation (1) can be modified to handle any inelastic strain (or eigenstrain). The scale transition relations are basically written as volume averages operations on stresses and strains. Hill [19] demonstrated, in a very general way, the equivalence between set (i.e. volume) averages and volume integrals. Hill's volume average relations over the mechanical states (also called, in the following of the present work: "consistency principles on mechanical states"), write as follows:

$$ \varepsilon^I = \langle \varepsilon^I \rangle, \quad (2) $$

$$ \sigma^I = \langle \sigma^I \rangle, \quad (3) $$

where the brackets <...> represent the arithmetical volume average. An extensive discussion about the choice of a mathematical method for achieving set averages can be found in [13]. In a fundamental work, Eshelby [6] studied
the behaviour of a single inclusion embedded in a homogeneous medium loaded at the infinite. He demonstrated that, if and only if the inclusion had an ellipsoidal shape, the local stresses and strains were homogeneous inside the BV, and fulfilled the following relation:

$$ \sigma^i - \sigma^I = -L^*: (\varepsilon^i - \varepsilon^I), \quad (4) $$

where $L^*$ is Hill's constraint tensor, which represents the interaction between the inclusion and its surrounding medium. It can be obtained from Eshelby tensor $S^I_{esh}$ or Morris tensor $E^I$, thanks to (Eq. 5), where $I$ is the 4th-order Identity tensor. Morris tensor is written as an integral that must be numerically computed, except for a few simplified cases. The interested reader can refer to [23, 24, 33], where expressions of Morris tensor and a detailed presentation of Eshelby-Kröner model, are provided.

$$ L^* = L^I : (S^I_{esh}^{-1} - I) = E^I^{-1} - L^I \quad (5) $$

In the EKSC scheme, the embedding medium is assumed to exhibit the macroscopic properties of the material. This implies that the homogenization relations providing the effective properties are implicit and must be solved numerically. Some models, such as the dilute approximation or the Mori-Tanaka model [29], can provide explicit expressions by considering the major phase (often the matrix phase) as the embedding medium. However, as these models neglect inter-particles interactions, they are not reliable for materials where the constituents have comparable volume fractions. The deviation between the effective properties, obtained by these two models (in auxiliary computations that will not be presented here), was not found extremely significant for the material considered in the present work (except for the in-plane Coulomb modulus). However, it was observed that the model had some influence over the predicted local mechanical states: Mori-Tanaka approximation slightly underestimates some local stresses components, in comparison with EKSC model.

Starting from the local and macroscopic constitutive responses, and using the scale transition relation given above (Eq. 4), one can express the local strains as a function of their macroscopic counterparts, as follows [10]:

$$ \varepsilon^i = (L^i + L^*)^{-1} : [(L^i + L^*) : \varepsilon^I + (L^i : \alpha^i - L^i : \alpha^i) \Delta T + (L^i : \eta^i - L^i : \eta^i) \chi] \quad (6) $$

Introducing Hill's average relation over the strains, and the macroscopic constitutive response, one then obtains the effective properties given by Eqs. (7) to (9). The model is self-consistent, i.e. leads to identical results (effective

properties and mechanical states) if written upon stresses instead of strains.

$$ L^I = L^I : \langle (L^i + L^*)^{-1} : (L^I + L^*) \rangle^{-1} \quad (7) $$

$$ \begin{aligned} \alpha^I &= L^{I-1} : \langle (L^i + L^*)^{-1} \rangle^{-1} \\ & \quad : \langle (L^i + L^*)^{-1} : L^i : \alpha^i \rangle \end{aligned} \quad (8) $$

$$ \begin{aligned} \eta^I &= L^{I-1} : \langle (L^i + L^*)^{-1} \rangle^{-1} \\ & \quad : \langle (L^i + L^*)^{-1} : L^i : \eta^i \rangle \end{aligned} \quad (9) $$

In the following, these three Eqs. (7 to 9) will be used to compute the effective mechanical properties of the composite, all along the cure process; the Eqs. (6) and (1) will then be used, in a time-differentiated version, to compute the local mechanical states in the constituents of the material.

### Application to the Hextool®

The present application is focused on a high-performance composite material, developed by Hexcel Composites©, for the design of composite parts to be used as a tooling material [18]. The Hextool® is made of unidirectional, rectangular shaped (50×8×0.15 mm³) reinforcing strips, being randomly arranged in the layout (see Fig. 1). The strips themselves are composed of AS4 fibers and M61 bismaleimid (BMI) matrix, set as a unidirectional ply.

The material presents itself as a thick layer (1.3 mm) exhibiting in-plane isotropy of its mechanical properties. In order to lighten the expressions, a “local” coordinate system $R_{xyz}$, oriented along the axis of the fibers constituting the strips, is introduced. This coordinate system is obtained by a rotation $\Theta$ around the 3-axis of the reference frame $R_{123}$ bound to the effective material. The overall fiber volume ratio is about 53 %, and thus the volume ratio of matrix is 47 % (the porosity is negligible). The matrix is assumed to be mainly contained in the strips, but a small fraction of the matrix (5 % of the total matrix volume) was also considered being located between the strips, to ensure the material’s consistency.

As a consequence, the structure of the material should be described at three different scales: microscopic, mesoscopic (intermediate) and macroscopic. The EKSC model will be used to perform a two-steps scale transition procedure (see Fig. 2). First, the effective properties of the reinforcing strip (mesoscopic scale) are estimated from those of the intra-reinforcements matrix and the carbon fibers (microscopic
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Schematic description of the studied material microstructure</figcaption></figure>

scale). Then, another homogenization procedure is achieved in order to find the behaviour of the material (macroscopic scale), from the properties of the extra-reinforcements matrix and those of the reinforcing strips (mescoscopic scale). The strong reinforcement ratio of the mesostructure (95 % of reinforcing strips) explains why the EKSC scheme was chosen, preferably to the MT model.

For the micro-meso scale transition, the fibers of infinite length, are supposedly aligned and have a volumic ratio $V^f= 56 \%$. The strip thereby exhibits a transverse isotropy of its properties. An analytical expression of Morris tensor, given in [11] for this specific configuration, enables a very fast computation of the reinforcing strips effective properties during the first step of the homogenization procedure (see Fig. 2). The non-negligible components of this tensor are prompted in Eq. 10; the negligible components were given a low constant value (resulting from a numerical evaluation of the tensor).

$$ 
\begin{aligned}
E_{22}^I &= E_{33}^I = \frac{3}{8L_{22}^I} + \frac{1}{4(L_{22}^I - L_{23}^I)} & E_{23}^I &= E_{32}^I = \frac{L_{22}^I + L_{33}^I}{8(L_{22}^I - L_{23}^I - L_{22}^I)} \\
E_{44}^I &= \frac{1}{8L_{22}^I} + \frac{1}{4(L_{22}^I - L_{23}^I)} & E_{55}^I &= \frac{1}{8L_{55}^I}
\end{aligned}
\quad (10) $$

For the meso-macro scale transition, the reinforcing strips, the volume fraction of which is $V^f=95 \%$, were supposed to be randomly oriented in the plane. This led to an in-plane isotropic behaviour at the macroscopic scale. Moreover, their morphology was assumed to correspond to an infinitely thin disc, in order to respect some fundamental hypothesis of the EKSC model, described in a previous paper with a similar material [26]. Those limits were

recently dispelled with a rigorous method, which will be described in future works. Moreover, extensive tests numerically demonstrated that this “thin-discs” approximation had neither significant influence on the predicted effective properties nor on the local mechanical states. This enabled also using an analytical expression of Morris tensor (which non-negligible terms are given in Eq. 11), historically established in [47]. This ensures a very fast computation of the macroscopic properties, at the second step of the homogenization procedure (see Fig. 2).

$$ E_{33}^I = \frac{1}{L_{33}^I} \quad E_{44}^I = E_{55}^I = \frac{1}{4L_{55}^I} \quad (11) $$

## Experimental investigation of the mechanical properties and cure kinetics of the BMI resin

The development of residual stresses in organic matrix composites is strongly related to the evolutions of the mechanical properties of the resin. In order to precisely describe those evolutions, the mechanical properties and cure kinetics of the M61 resin were experimentally investigated by S. Terekhina and K. Szymanska, respectively.

### Cure kinetics

### Models

The M61 resin is a toughened thermoset resin cured at 190° C, following Hexcel Composites’ recommendations [18].

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Principle of the two-steps scale transition procedure</figcaption></figure>
The polymerisation reaction is exothermic. As a consequence, one can define the conversion degree $\chi$ as the following ratio:

$$ \chi = \frac{\Delta H}{\Delta H_{tot}}, \qquad (12) $$

where $\Delta H$ is the enthalpy released since the beginning of the reaction, whereas $\Delta H_{tot}$ is the enthalpy released during the whole process. This conversion degree is a state variable; it allows the parameterization of the cure kinetics, but also of the thermal and mechanical properties of the resin. Enthalpy variations are easily measured using Differential Scanning Calorimetry (DSC), that's why this technique is widely used to experimentally investigate the cure kinetics of polymers.

The conversion rate $\dot{\chi}$ is generally defined as a function of both the absolute temperature $T$ and conversion degree:

$$ \dot{\chi} = \frac{d\chi}{dt} = f(\chi, T) \qquad (13) $$

Rigorous mechanistic descriptions of the cure kinetics can be established, from the knowledge of the elementary chemical reactions beneath the polymerization process. This can be accomplished, for a given resin formulation, using techniques such as near infrared spectroscopy [38]. However, as the resin formulation and required instruments are not always available, several phenomenological models have been proposed in the literature to describe the cure kinetics of arbitrary polymers from DSC measurements. One can quote the widely-used model of Kamal and Sourour [22], which fundamental relation writes:

$$ \dot{\chi} = (K_1 + K_2\chi^m)(1 - \chi)^n, \qquad (14) $$

where m and n are adjustment parameters, whereas $K_1$ and $K_2$ are Arrhenius-type functions of the temperature. Another model, proposed by Bailleul [2], enables a more accurate description of the cure kinetics, through:

$$ \dot{\chi} = K(T) \cdot f(\chi) = \exp\left(-\frac{E_a}{R \cdot T}\right) \cdot G(\chi), \qquad (15) $$

with $G(\chi)$ a polynomial expression, and $K(T)$ an Arrhenius-type function, depending on two parameters: the temperature and an activation energy. For all these models, the Arrhenius-type functions (e.g. the activations energies) are supposed independent of the conversion degree, and merely correspond to optimisation parameters. This lack of strong physical meaning is also the source of some inconsistencies occurring between the interpretations of isothermal and anisothermal data. To avoid those inconsistencies, Sbirrazzuoli and Vyazovkin [42] proposed to write the activation energy (in Eq. 15) as a function of the conversion degree ("Model Free Kinetics", MFK). This method was retained as

the most relevant to describe the complexity of the cure kinetics of the BMI resin.

## Experimental investigation

The cure kinetics were investigated through Differential Scanning Calorimetry (DSC), using a Mettler Toledo DSC1 STARe System, under a neutral atmosphere (nitrogen) for avoiding resin oxidation phenomena. Dynamic and isothermal runs were operated on 14.40±0.15 mg-weighting samples of M61 resin systems. The dynamic runs were performed at the following heating rates: 5 °C/min, 10 °C/min, 15 °C/min and 20 °C/min. Two isothermal runs were also performed at 190 °C and 220°, during four hours. After each (isothermal or dynamic) run, a second dynamic run was performed in order to establish the conversion residue, which is equal to 8.7 % for the 190 °C isothermal run, and null for all the others. The resulting conversion-temperature curves for the dynamic runs, as well as the conversion-time curves for the static runs, are given in Fig. 3a and b, respectively.

## Parameterisation of the cure kinetics

The activation energy and the conversion function were computed upon the five considered dynamic runs, using the "isoconversional method", as defined by Sbirrazzuoli and Vyazovkin. According to this definition, the activation energy was computed by minimizing the following functional over several heating programs (temperature instructions):

$$ \begin{aligned} \Phi(E_a(\chi)) &= \sum_{i=1}^{n} \sum_{j \neq i}^{n} \frac{J[E_a, T_i(t_\chi)]}{J[E_a, T_j(t_\chi)]} \\ \text{with } J[E_a, T_i(t_\chi)] &= \int_{t_\chi - \Delta\chi}^{t_\chi} \exp\left(-\frac{E_a(\chi)}{R \cdot T_i(t_\chi)}\right) \end{aligned} \qquad (16) $$

This operation was performed with a dedicated module of the Mettler Toledo DSC analysis system, and corresponds to the curve shown in Fig. 4a. The conversion functions were then deduced from each of the anisothermal runs, and plotted on Fig. 4b. Unexpectedly, these four curves where found significantly different, which can be attributed to the definition of the activation energy by an optimisation process. The conversion function obtained from the 5 °C/min run was arbitrarily chosen as a reference to establish the cure kinetics.

In order to simplify the implementation, the curves corresponding to the activation energy $E_a(\chi)$ and to the conversion function $G(\chi)$ were fitted with polynomial expressions. A special attention was paid to this operation,
<figure><img src="image_5.png" /><figcaption>**Fig. 3** (a) Conversion-temperature curves for dynamic runs (b) Conversion-time curves for static runs</figcaption></figure>

as the functions are highly contrasted. Moreover, a slight change in the activation energy has a strong influence on the conversion rate, as a result of the exponential form of the Arrhenius-type function.

## Description of the considered manufacturing process

The manufacturing process recommended by Hexcel Composites, shown on Fig. 6, was considered in the present work. A cure plateau of four hours at 190 °C (preceded by a short plateau at 150 °C) is applied to the part, which is then removed from the mould before being post-cured during sixteen hours at 220 °C. The MFK method was applied to this cure process, and the corresponding values of the conversion rate was numerically integrated (with the method of rectangles), in order to obtain the conversion degree reported on Fig. 5.

One can notice that the predicted conversion rates during the 190 °C and 220 °C isothermal phases are slightly different from those obtained with the corresponding DSC runs (Fig. 3b). The reason for these discrepancies is that the MFK model was optimised on dynamic runs, for which the

temperature stays above the glass transition temperature; this is not the case for isothermal runs, where the glass transition occurs and sharply slows down the conversion rate. In order to correctly model this effect, diffusion-controlled reaction kinetics should be applied (see, for instance, Sanford and McCullough [41], or Matsuoka et al. [27]).

Nevertheless, being given the hypotheses of this study, and the cure process considered here, only three values are assumed to significantly influence the predicted stress states: the temperature at which the gelation of the resin occurs (190 °C), the conversion degree after the first cure cycle (91.3 %) and the final conversion degree (100 %). The slight change of the conversion degree during the 190 °C–20 °C–220 °C cool-down and heat-up phase is supposed to have a nearly negligible influence over the predicted stress states. Consequently, the results can be uncoupled from the time-dependent cure kinetics, and related only to those critical “state indicators”. As those values are in good agreement with the results of the isothermal DSC runs (Fig. 3b), this cure kinetics model will be considered as acceptable for the purpose of this study.

<figure><img src="image_6.png" /><figcaption>**Fig. 4** (a) Activation energy (isoconversional method) (b) Conversion functions for the MFK method</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 5 Manufacturing process considered (temperature, conversion degree)</figcaption></figure>

## Mechanical properties of the BMI resin system

### Elastic modulus

The polymerisation of the resin along the cure process induces deep changes of its physical state and, subsequently, mechanical properties; from an initial liquid state, the resin turns into an amorphous solid, then into a glassy solid. The liquid-solid transition, named “gelation point”, is considered by most authors as the initiation of residual stresses [31, 32]. All these transformations yield extreme evolutions of the viscous and elastic moduli of the resin, which strongly affect the level of residual stresses created. However, in this study, only the long-term elastic response of the resin, and its evolutions (the toughening effect), will be considered; the viscous and history effects will be not be taken into account.

In order to model the toughening of the resin, most authors suppose a linear dependence of the Young's modulus on the conversion degree, from the gelation point (considering that the viscous stresses created before this point are instantly relaxed), with a constant Poisson's ratio. Boetti and Gillespie [4] use an enhanced version of the expression, historically provided by Dillman and Seferis [5], where a second-order term enables one to take into account the competition between the toughening of the resin and the viscous relaxation. However, Adolf and Chambers [1] showed that the percolation theory gives a more suitable (and physically-consistent) description of the equilibrium stiffness, in the amorphous domain, with the following expression for the Young's modulus (and, still, a constant Poisson's ratio):

$$E(T, \chi) = E_{\text{cured}}(T) \cdot \left( \frac{\chi^2 - \chi_{\text{gel}}^2}{1 - \chi_{\text{gel}}^2} \right)^{8/3}, \quad (17)$$

with $χ_{gel}$ the conversion degree at the gelation point and $E_{cured}(T)$ the modulus of the fully cured resin, at a given

temperature. The influence of glass transition temperature should be accounted for via a corrective term, such as the one proposed by Msallem et al. [31, 32]; however, this level of complexity could not be reached in this paper, due to the lack of experimental data on the stiffness of the resin in the amorphous state. Consequently, the expression given by Eq. 17 will be used in the following of the present work, but still keeping in mind the potential sensitiveness of the results to the glass transition temperature.

The evolution of the stiffness of the fully cured resin as a function of the temperature (namely, $E_{cured}(T)$) was determined through Dynamic Mechanical Analysis (DMA). DMA testing is a standard thermo-mechanical analysis technique for characterizing viscoelastic properties of polymers and polymer-based materials as a function of frequency and temperature [28]. It consists in the observation of the time-dependent behaviour of a material under dynamic periodic (sinusoidal) strain or stress. According to Etienne et al. [7], one can extrapolate the complex shear modulus $G^*$ from the stress and strain measurements: $G^*=G'+jG''$, where $G'$ is the storage (elastic component) modulus, $G''$ standing for the loss (viscous component) modulus. The loss (or damping) factor is denoted by $\tan \delta=G''/G'$, where $\delta$ represents the phase shift between stress and strain. In this paper, the tests were carried out at controlled strains in the linear domain of viscoelasticity of the material, and the corresponding stress was measured.

The DMA50 0.1 dB Metravib dynamic mechanical analyzer was used for measuring the temperature-dependent elastic modulus (instantaneous stiffness) of the fully cured resin. As the fully cured resin is in a vitrous state for the temperature range considered in this DMA test, the instantaneous and equilibrium (long-term) moduli were supposed equal. Dynamic testing was performed on a 26×4.1×1.5 mm³ rectangular bar, in the tensile/compressive mode, at a frequency of 1 Hz. The temperature ranged from -20 °C to 230 °C at a heating rate of 1.5°C/min, under a flow of nitrogen. Figure 6 sums up the evolution of the stiffness versus the temperature and conversion degree. A decrease of Young modulus E with rising of temperature was observed [45].

DMA tests were also performed on the uncured resin for estimating the gelation point. As the resin was initially in a liquid state, another standard was used: the tests were performed on the shear mode (sinusoidal shear strains) at a frequency of 1 Hz, using two plane plates (one rotating plate, the other being fixed). A heat ramp at 5 °C/min, from 24 °C to 350 °C, was applied at the same time, in order to polymerise the resin. The resulting log(G') versus time curve remains horizontal during the first half-an-hour of the experiment, then shows a sudden increase of the storage modulus, rapidly followed by a linear evolution. The gelation time was then defined as the intersection between the
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Young's modulus of the M61 resin, as a function of the temperature and conversion degree</figcaption></figure>

corresponding linear interpolation line, and the axis of time. By comparison with the 5 °C/min DSC runs, it was then easy to determine the conversion degree at this point (χ<sub>gel</sub>), which was found equal to 45 %.

### Thermal expansion

The CTE of the resin also experiences severe evolutions during the curing process: it increases with temperature and decreases with conversion degree. The evolution of the fully cured resin CTE was given by Hexcel for three ranges of temperatures. A linear data interpolation was achieved using the middle point of each temperature range (see Table 1).

The dependence upon the conversion degree could also have a significant effect on the results, as the resin undergoes temperature changes while still not being fully cured. This dependency was represented by a linear evolution of the CTE as a function of conversion degree, the CTE of the totally uncured resin having been fixed at a mean value of 220.10<sup>-6</sup>/K.

### Cure shrinkage

The polymerisation reaction corresponds to the creation of covalent bonds between the macromolecules, and a constriction of the amorphous network. This induces important bulk shrinkage of the resin, in the order of 3 % up to 9 %

(depending on the chemical formulation of the resin). This chemical shrinkage is comparable to the thermal shrinkage and thus, shall not be neglected. An estimation of 5.7 % for the bulk shrinkage of our BMI resin was obtained (unpublished results), using a dedicated equipment (a PVTα mold), also used by Msallem et al. This shrinkage was also supposedly proportional to the conversion coefficient, leading to a linear Coefficient of Chemical Expansion (CCE) η being equal to -1.67 %. The hypothesis of a constant CCE along the cure process is, yet, quite discussable from a physical point of view (see Rabearison et al. [39]).

### Effective properties and internal stresses computation during the cure process

#### Computation of the effective properties

The above presented scale transition procedure is used for computing the effective properties of the material during the cure process. The properties at room temperature (20 °C) of the fully cured Hextool, and of its constituents, are summarized in Table 2. A few experimental results can confirm the estimated effective properties.

The in-plane elastic modulus was determined through tensile tests (EN 2561). It varies between 41 and 45.7 GPa. The in-plane shear modulus was found equal to 14.4 GPa (Iosipescu test ASTM D5379), and the CTEs were also measured at room temperature (ASTM E831 and E228). Large relative uncertainties (in the order of 10 %) are observed upon the in-plane properties, as the sample sizes are probably too small to ensure a good representativity. Taking into account those uncertainties, one can say that the experimental and calculated values are close enough. Using the Mori-Tanaka model instead of the self-consistent approach gives a slightly better estimate of the shear moduli, but lower CTEs. The discrepancies between the numerical and experimental results can also be attributed to the out-of-plane waviness of the reinforcing strips (which deteriorate their in-plane apparent properties).

To follow the evolution of these properties during the manufacturing process, the temperature and conversion degrees presented on Fig. 5 are considered. The properties of the matrix are then computed and introduced in the homogenization scheme, in order to obtain the macroscopic

<table><caption>**Table 1** Coefficients of thermal expansion of the fully cured BMI resin</caption><thead><tr><th>20&lt;T&lt;70 °C</th><th>70&lt;T&lt;120 °C</th><th>120&lt;T&lt;170 °C</th><th>Interpolation [10<sup>-6</sup>/ °C]</th></tr></thead><tbody><tr><td>41.4×10<sup>-6</sup>/ °C</td><td>49.3×10<sup>-6</sup>/ °C</td><td>58.8×10<sup>-6</sup>/ °C</td><td>α<sub>cured</sub>(T)=33.285+0.1743×T</td></tr></tbody></table>
<table><caption>**Table 2** Thermo-chemio-elastic effective properties of the Hextool and its constituents</caption><thead><tr><th></th><th colspan="6">Mechanical moduli</th></tr><tr><th></th><th>E<sub>x</sub> [GPa]</th><th>E<sub>y</sub>,E<sub>z</sub> [GPa]</th><th>v<sub>xy</sub>, v<sub>xz</sub></th><th>G<sub>xy</sub>,G<sub>xz</sub> [GPa]</th><th>v<sub>yz</sub></th><th>G<sub>yz</sub> [GPa]</th></tr></thead><tbody><tr><td>M61 resin</td><td>3.5</td><td>3.5</td><td>0.40</td><td>1.6</td><td>0.40</td><td>1.6</td></tr><tr><td>AS4 fiber <sup>a</sup></td><td>207</td><td>20.7</td><td>0.20</td><td>27.6</td><td>0.05</td><td>6.9</td></tr><tr><td>UD strips</td><td>117.0</td><td>9.5</td><td>0.28</td><td>7.6</td><td>0.56</td><td>3.0</td></tr><tr><td></td><td>E<sub>1</sub>,E<sub>2</sub> [GPa]</td><td>E<sub>3</sub> [GPa]</td><td>v<sub>12</sub></td><td>G<sub>12</sub> [GPa]</td><td>v<sub>13</sub>,v<sub>23</sub></td><td>G<sub>13</sub>,G<sub>23</sub> [GPa]</td></tr><tr><td>Hextool</td><td>46.3</td><td>12.3</td><td>0.27</td><td>18.2</td><td>0.37</td><td>3.9</td></tr><tr><td>Exp. <sup>b</sup></td><td>41~45.7</td><td>X</td><td>X</td><td>14.4</td><td>X</td><td>X</td></tr><tr><td></td><td>CTE</td><td></td><td>CCE</td><td></td><td></td><td></td></tr><tr><td></td><td>α<sub>x</sub> [10<sup>-6</sup>/K]</td><td>α<sub>y</sub>, α<sub>z</sub> [10<sup>-6</sup>/K]</td><td>η<sub>x</sub> [10<sup>-6</sup>]</td><td>η<sub>y</sub>,η<sub>z</sub> [10<sup>-6</sup>]</td><td></td><td></td></tr><tr><td>M61 resin</td><td>38.2</td><td>38.2</td><td>-19000</td><td>-19000</td><td></td><td></td></tr><tr><td>AS4 fiber</td><td>-0.35</td><td>11.8</td><td>0</td><td>0</td><td></td><td></td></tr><tr><td>UD strips</td><td>0.32</td><td>27.4</td><td>-353</td><td>-9940</td><td></td><td></td></tr><tr><td></td><td>α<sub>1</sub>,α<sub>2</sub> [10<sup>-6</sup>/K]</td><td>α<sub>3</sub> [10<sup>-6</sup>/K]</td><td>η<sub>1</sub>,η<sub>2</sub> [10<sup>-6</sup>]</td><td>η<sub>3</sub> [10<sup>-6</sup>]</td><td></td><td></td></tr><tr><td>Hextool</td><td>2.98</td><td>42.5</td><td>-1323</td><td>-15894</td><td></td><td></td></tr><tr><td>Exp. <sup>b</sup></td><td>3.5~5</td><td>45~52</td><td>X</td><td>X</td><td></td><td></td></tr></tbody></table>

effective properties given in Fig. 7, for each time points. The method is very cost-effective, as the computation of one set of effective properties lasts less than ten seconds on a recent computer.

All properties are set to zero until the gelation of the matrix occurs (χ<sub>gel</sub>=45 %, t<sub>gel</sub>=129 min): i.e. the resin is considered as a liquid before this point. The sudden change of some properties at the gelation point is a numerical issue, and a side effect of the analytical expressions used for the Morris tensor. One can then observe a significant evolution of the properties as the matrix stiffens, followed by a stabilisation; weaker evolutions can also be observed during the heat-up and the cool-down steps of the process. The curves show that out-of-plane properties are clearly matrix-controlled, whereas in-plane properties are mainly related to the fibers' high stiffness and weak coefficients of expansion (particularly the in-plane modulus E<sub>1</sub>). However, this last observation should be moderated, as the behaviour of the matrix is strongly viscoelastic just beyond the gelation point.

### Computation of the local residual stresses

As the properties of the material and the constituents vary along the cure process, the Eqs. (1) and (6), giving the constitutive responses, must be written in a differential form, then integrated over time. Considering a small time step δt, the increments of temperature, conversion degree, strain and stress are denoted as follows: δT, δχ, δε and δσ. Moreover, the effect of a chemical

shrinkage is taken into account, which yields the following relations:

$$ \delta\sigma^k = L^k : (\delta\varepsilon^k - \alpha^k \delta T - \eta^k \delta\chi), \text{ where } k = i, I \quad (18) $$

$$ \delta\varepsilon^i = (L^i + L^*)^{-1} : [(L^i + L^*) : \delta\varepsilon^i + (L^i : \alpha^i - L^i : \alpha^i)\delta T + (L^i : \eta^i - L^i : \eta^i)\delta\chi] \quad (19) $$

Those increments are computed, using the mechanical properties computed at the beginning of the time step, then cumulated to obtain the mechanical states (with the same method as the evaluation of the conversion degree). The time increments are small enough to ensure a precision higher than 1 % over the integrated mechanical states.

### Results and discussion

For the range of this study, a fictive thin part is considered, and is affected idealised boundary conditions, representing mould-part interactions. This enables to free oneself of any structural or transient effects, and to obtain uniform temperature and macroscopic states fields. Consequently, time appears as a virtual parameter: the results only depend of a few “state indicators” described in “Description of the considered manufacturing process”, and of the constitutive response of the resin (which is itself described by state variables). Three typical cases are considered here: (i) a cure without external mechanical loading of the material (case #1); (ii) a cure with prevented in-plane macroscopic strains

<footnote><sup>a</sup> Elastic moduli provided by White and Kim [48], CTEs provided by Hexcel Composites [18]</footnote>

<footnote><sup>b</sup> Values provided by Hexcel Composites [18]</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 7** Evolutions of the material effective properties during the cure process</figcaption></figure>

(case #2); (iii) and a cure with imposed thermal strains in the plane of the ply (case #3). For these two last cases, the removal from the mould will also be considered, by releasing the macroscopic stresses. Any of these loads respects the in-plane isotropy of the material, so the local states are independent of the orientation Θ of the considered strip (see Lacoste et al. [26])

Process case #1

For this case, the macroscopic stresses are set to zero during
the whole cure process. The gelation point (t_gel=129 min) is

considered as the initiation of residual stresses and strains.
Figure 8 shows the evolution of stresses in the constituents,
at the microscopic and mesoscopic scale. The stresses fol-
low the evolution of the thermo-chemical shrinkage, with a
first peak after the first cure cycle, and another peak at the
end of the post-cure stand. The stresses are very weak when
the temperature reaches 220 °C, which roughly corresponds
to the stress-free temperature. During the first cure cycle, the
curves also show a neat decomposition of chemical and
thermal shrinkage, the second yielding nearly 300 % more
residual stresses than the first. This result is consistent with
results from other authors (see, for instance, Zhao et al.
[49]). However, during the post-cure cycle, the thermal expansion from 190 °C to 220 °C nearly nullifies the effect of the chemical shrinkage, so that the post-cure temperature can roughly be considered as a stress-free temperature. This result justifies the hypothesis of neglecting the viscoelastic relaxation at high temperature, but this should not be transposed to the viscoelastic effects occurring during the cooling. Also, the chemical shrinkage shall not be neglected, as it plays an important role over the strain states, particularly in the z-direction (out-of-plane).

The microscopic constituents (fibers and intra-reinforcement matrix) undergo more important stresses than mesoscopic elements. In particular, the intra-reinforcement matrix undergoes about 75 MPa in the x and y directions, which might lead to its failure or to fiber-matrix debonding. In order to determine the risk of matrix failure, a Tsai-Wu failure criterion is used, and the inverse of the strength factor is plotted on Fig. 8. The parameters of the criterion are taken from the engineering strength of a high-strength epoxy

matrix, given by Fiedler et al. [9]. This shows that residual stresses in the matrix can reach 64 % of its strength at failure, which is critical. Moreover, the self-consistent approach provides *average* (mean-field) values of the local stresses, so these may reach even stronger values in the actual material. This subject will be discussed later. On the contrary, the self-consistent approach gives a more realistic description for rigid elements (fibers and reinforcing strips), as they can be viewed as inclusions with well-defined morphologies. In the fiber direction, they experience compression stresses of –60 MPa and –170 MPa, respectively. These loads may not cause local degradations because of their high strength, but there is a risk of micro-buckling of these elements. Despite weaker stresses, there is a greater danger for reinforcing strips, because of their pre-existent waviness. This could also be detrimental for the surface state of the material and would increase the need of a surface treatment after fabrication.

<figure><img src="image_10.png" /><figcaption>**Fig. 8** Time dependent evolutions of the local stress states during the curing process, computed for the case #1.</figcaption></figure>
Process case #2

For this case, the macroscopic in-plane strains were set to zero during the first cure cycle, and the cumulated macroscopic stresses were released at the end of the cure cycle. The post-cure step was treated in the same way as in case #1 (stand-alone post-cure). Given the in-plane symmetry of the problem, the macroscopic strain increments of the Hextool are then given by Eq. 20. This results in tensile stress states in the Hextool during the cure phase.

$$ \left| \begin{array}{l} \delta\epsilon_{11}^I = \delta\epsilon_{22}^I = 0 \\ \delta\epsilon_{33}^I = (\alpha_3^I \cdot \delta T + \eta_3^I \cdot \delta\chi) + (\alpha_1^I \cdot \delta T + \eta_1^I \cdot \delta\chi) \times 2L_{13}^I/L_{33}^I \end{array} \right. \quad (20) $$

The resulting stresses in the composite and its constituents are depicted in Fig. 9. The stress states discontinuity corresponds to the removal of the material from the mould. It can be seen that the ultimate stresses experienced at the end of the manufacturing process #2 are very close to those observed for the stress-free cure. This shows that the

macroscopic stresses applied before releasing, combined to the slight evolution of the in-plane modulus of the composite during the cure, have little influence on the final result. One can also observe that the stress states experienced by the matrix before releasing are very close to those observed for the stress-free cure, which can be interpreted by the very weak stiffness of the matrix. The in-plane mechanical stresses applied to the composite are mainly transmitted to the rigid elements in the x-direction. This results to weaker stress states for those elements, which is beneficial for the material.

Process case #3

For this case, the macroscopic in-plane strains were assumed proportional to the temperature increments, the proportionality factor being equal to the CTE of steel ($\alpha_{steel} = 12.3 \times 10^{-6}/°C$). The strain increments in the composite are then computed via Eq. 21, and the resulting stress states are

<figure><img src="image_11.png" /><figcaption>Fig. 9 Time dependent evolutions of the local stress states during the curing process, computed for the case #2.</figcaption></figure>
given in Fig. 10. These stresses coarsely represent the inter-action between a composite part and a massive steel mould.

$$ \begin{gathered} |\delta\epsilon'_{11}| = |\delta\epsilon'_{22}| = \alpha^{\text{steel}} \cdot \delta T \\ |\delta\epsilon'_{33}| = (\alpha_3^I \cdot \delta T + \eta_3^I \cdot \delta\chi) + ((\alpha_1^I - \alpha^{\text{steel}}) \cdot \delta T + \eta_1^I \cdot \delta\chi) \times 2L_{13}^I / L_{33}^I \end{gathered} \quad (21) $$

The mould-part interactions result in tensile stress states during the cure plateau at 190 °C (due to the chemical shrinkage of the material), whereas strong compressive stresses occur during the cooling of the part.

At first hand, one must notice that those boundary conditions are unrealistically severe, as they yield very important compressive stresses in the composite (about -93 MPa). In any real manufacturing process, the adhesion of the mould and the part is not sufficiently strong to sustain those stresses, and debonding would probably occur.

The stress states in the matrix constituents are not strongly affected by the mould-part interactions, and their influence is mostly beneficial to the matrix, as they reduce the stress states experienced during the cure process. On the other hand, for the rigid elements, the external compressive load yields important stresses in the x-direction, before removal from the mould: -230 MPa in the reinforcing strips and -445 MPa in the fibers. These stresses are known to depend on the tooling material's CTE, but also on the cooling rate, due to viscoelastic effects [25]. This shows that the choice of a proper tooling material and cure cycle is fundamental to avoid material degradations during the cure cycles.

Indeed, with such levels of compression, one can fear a buckling of the rigid elements, particularly the reinforcing strips ("meso-buckling") which already have an important initial waviness. Rosen [40] gave a very coarse estimate of the microbuckling load, by simply relating the critical

<figure><img src="image_12.png" /><figcaption>Fig. 10 Time dependent evolutions of the local stress states during the curing process, computed for the case #3.</figcaption></figure>
compressive stress in the fiber to the transverse shear modulus of the matrix:

$$ \sigma_{\text{critical}}^{f} = \frac{1}{V^{f}V^{m}} G_{13}^{m} \qquad (22) $$

This estimate, quite reassuring, is nevertheless considered as excessive by most authors [34, 43, 44], since it does not take into account the initial reinforcement waviness. Many studies on this subject can be found in the literature (see, for instance, works by Grandidier et al. [14]), but cannot be easily transposed to this particular material because of its very specific mesostructure. We may say, as a comparison, that the experimental compressive strength of the composite (330 MPa, EN-2850-B test) is much higher than the loads considered in the present section. Moreover, the risk of mesobuckling shall be moderated considering that the stresses applied for this virtual process case are unrealistically severe.

## Discussion of the results

The numerical results obtained with the multiscale simulation could not (yet) be compared with experimental results, except for some of the mechanical properties. The random mesostructure makes it difficult to obtain a reliable measure of the local strains or stresses with classical methods (Fiber Bragg Gratings, incremental hole drilling method...).

Despite this lack of experimental confrontation, the simulation results highlight several interesting facts. The first is the effect of the chemical shrinkage, which has a significant influence over the local stresses after the first cure cycle, but which seems to be erased by the post-cure cycle (for our set of material parameters). Secondly, the evolution of the local properties affect the final stress state, and this final state doesn't seem to depend much on the history of the mechanical loads (or tool-part interaction) applied during the curing process (after removal from the mold, of course).

The results also highlight two main concerns about the strength of the material under residual stresses and, thereafter, about its durability under service-life loads: the risk of micro- and meso-buckling of the rigid elements, under compressive residual stresses, and the failure of the resin under tensile stresses. Unfortunately, both mechanisms are difficult to represent with the mean-field approach, used in the present work, which neglects the inter-particles interactions. As a consequence, the average values obtained for the local stresses do not provide sufficient information for treating phenomenon such as the failure or the progressive damage of the material. The computation of resin-reinforcement interface stresses would be a possible and easy enhancement to the model, as this is a corollary of the historical Eshelby's inclusion problem. However, this simplified problem does not take into account the perturbations induced by the

proximity of the inclusions. For instance, at the mesoscopic scale, the random shear stress experienced by the resin due to the disorientation of the UD strips, should not be neglected, as they are a key damage mechanism of the Hextool. In order to properly take into account the geometrical information upon the microstructure and mesostructure, several approaches can be considered. If the local structure corresponds to a periodic pattern, two complementary methods are widely used: the multilevel FE² method, which combines macroscopic and periodic degrees of freedom [8, 15, 16, 21], and the FFT-based method, which uses the Fast Fourier Transform to compute the Green's influence operator [30]. These methods are relevant for the UD strips, which microstructure is much more structured. However, in order to describe the mesostructure of the Hextool, with its random nature and the reinforcements' waviness, the use of a full-field approach seems unavoidable.

## Conclusions and perspectives

A scale transition procedure based on Eshelby-Kröner self-consistent model was established, in order to describe the multi-scale thermoelastic behaviour of a composite material exhibiting a very specific microstructure. This procedure was applied to the prediction of the time dependent evolutions of the effective properties and local stress states in the composite, during the cure process. An experimental investigation was led to determine the cure kinetics and the properties of the BMI resin along the cure process. The obtained results were taken into account in the homogenization procedure devoted to compute the composite's effective properties. These properties were found in good agreement with experimental results on the cured material.

The local stresses created by the thermal and chemical shrinkage could be computed for three processing conditions. The thermal shrinkage was found much more significant over local stresses states than the chemical shrinkage; nevertheless, this last shall not be neglected as it has a significant influence on out-of-plane strains. Also, the influence of viscoelastic relaxation during the cool-down periods remains to be characterized.

The influence of macroscopic loads on the local mechanical states was also investigated. At first hand, it seems that they do not fairly influence the residual stresses after removal from the mould, which can be explained by the relative stability of the in-plane stiffness during the cure. However, the macroscopic loading has a tremendous influence on the stresses experienced by the rigid elements along the cure process. It was shown that using a rigid mould with a strong CTE might result in important material degradations before removal from the mould, due to microbuckling effects.
36. Parlevliet PP, Bersee HEN, Beukers A (2007) Residual stresses in thermoplastic composites – a study of the literature. Part II: Experimental techniques. Compos Appl Sci Manuf 38:651–665

37. Parlevliet PP, Bersee HEN, Beukers A (2007) Residual stresses in thermoplastic composites – a study of the literature. Part III: Effects of thermal residual stresses. Compos Appl Sci Manuf 38:1581–1596

38. Pascault JP, Sautereau H, Verdu J, Williams R (2002) Thermosetting polymers. Marcel Dekker Inc. p. 157–196

39. Rabearison N, Jochum C, Grandidier JC (2009) A FEM coupling model for properties prediction during the curing of an epoxy matrix. Comput Mater Sci 45(3):715–724

40. Rosen VW (1965) Mechanics of composite strengthening. In: Fiber Composite Materials. American Society of Metals, Metals Park, Ohio, pp 37–75

41. Sanford WM, McCullough RL (1990) A free-volume-based approach to modeling thermoset cure behavior. J Polymer Sci, Part B: Polymer Phys 28:973–1000

42. Sbirrazzuoli N, Vyazovkin S (2002) Learning about epoxy cure mechanisms from isoconversional analysis of DSC data. Thermochim Acta 388(1–2):289–298

43. Schultheisz CR, Waas AM (1996) Compressive failure of composites, part I: testing and micromechanical theories. Prog Aerosp Sci 32(1):1–42

44. Schultheisz CR, Waas AM (1996) Compressive failure of composites, part II: Experimental studies. Prog Aerosp Sci 32(1):43–78

45. Terekhina S, Salvia M, Fouvry S (2011) Contact fatigue and wear behaviour of bismaleimide polymer subjected to fretting loading under various temperature conditions. Tribol Int 44(4):396–408

46. Tsai SW, Hahn HT (1980) Introduction to composite materials. Technomic, Westport

47. Welzel U, Fréour S, Mittemeijer EJ (2005) Direction-Dependent Elastic Grain-Interaction Models – A Comparative Study. Philos Mag 85:2391–2414

48. White SR, Kim KS (1998) Process-induced residual stress analysis of AS4/3501-6 composite material. Mech Compos Mater Struct 5:153–186

49. Zhao LG, Warrior NA, Long AC (2007) A thermo-viscoelastic analysis of process-induced residual stress in fibre-reinforced polymer-matrix composites. Mater Sci Eng, A 452–453:483–498