Numerical Analysis of Micro-Residual Stresses in a Carbon/Epoxy Polymer Matrix Composite during Curing Process

The manufacturing process in thermoset-based carbon fiber-reinforced polymers (CFRPs) usually requires a curing stage where the material is transformed from a gel state to a monolithic solid state. During the curing process, micro-residual stresses are developed in the material due to the different chemical–thermal–mechanical properties of the fiber and of the polymer, reducing the mechanical performance of the composite material compared to the nominal performance. In this study, computational micromechanics is used to analyze the micro-residual stresses development and to predict its influence on the mechanical performance of a pre-impregnated unidirectional CFRP made of T700-fibers and an aeronautical grade epoxy. The numerical model of a representative volume element (RVE) was developed in the commercial software Abaqus® and user-subroutines are used to simulate the thermo-curing process coupled with the mechanical constitutive model. Experimental characterization of the bulk resin properties and curing behavior was made to setup the models. The higher micro-residual stresses occur at the thinner fiber gaps, acting as triggers to failure propagation during mechanical loading. These micro-residual stresses achieve peak values above the yield stress of the resin 55 MPa, but without achieving damage. These micro-residual stresses reduce the transverse strength by at least 10%, while the elastic properties remain almost unaffected. The numerical results of the effective properties show a good agreement with the macro-scale experimentally measured properties at coupon level, including transverse tensile, longitudinal shear and transverse shear moduli and strengths, and minor in-plane and transverse Poisson’s ratios. A sensitivity analysis was performed on the thermal expansion coefficient, chemical shrinkage, resin elastic modulus and cure temperature. All these parameters change the micro-residual stress levels and reduce the strength properties.


Introduction
Fiber-reinforced polymers (FRPs) are commonly used in aerospace and high-performance applications due to their higher strength-to-weight ratio. They are made with a fiber reinforcement embedded in a polymeric matrix. This reinforcement ranges from short fibers, continuous long fibers and fabrics. The most common fiber reinforcement materials are glass fibers (GFRPs), carbon fibers (CFRPs) and other organic fibers such as Kevlar, Dyneema, etc. [1]. Composite materials possess some advantages related to traditional (metallic) materials in addition to the weight reduction, including toughening for impact, fatigue resistance, corrosion resistance, electromagnetic transparency, erosion and wear resistance, acoustic and vibration damping, low thermal expansion, among others [1][2][3][4]. However, the most important advantage of FRPs is their tailoring ability to fit design requirements.

Properties Value
Fiber diameter (mm) 7 × 10 −3 Density (kg/m 3 ) 1800 Specific Heat (kJ/kg·K) 0.752 Thermal Conductivity (W/m·K) 9.38 The polymer matrix used in this study corresponds to an aeronautic grade epoxy resin. Table 2 shows the material parameters used for model setup, where CTE is the coefficient of thermal expansion, T g is the glass transition temperature, E m is the elastic modulus, v m is the elastic Poisson's ratio, v p is the plastic Poisson's ratio, S y+ and S y− are the tensile and compressive yield strengths and S u+ and S u− are the tensile and compressive ultimate strengths, ε f is the failure strain at S u+ and G f lm is the fracture toughness. The mechanical properties for the pure resin were extracted from the datasheet, from dog-bone tensile tests and from the literature [10,27]. The properties reported in Table 2 correspond to the nominal properties of the matrix in as-manufactured conditions. (Degree of cure ϕ = 1). The fibers are modeled using a transversely isotropic elastic constitutive relation. This is a common approach for RVE modeling of composites whose transverse failure is controlled by the polymer matrix and its interface with the fibers. Equation (1) shows the elastic compliance matrix considering that direction 1 is the fiber direction. The Poisson's ratio v T can be estimated from G T = E T /(2(1 + v T )) Polymers 2022, 14, 2653

Polymer Matrix
The curing process in thermoset resins is usually an exothermic reaction process. It is normally characterized by the total curing reaction enthalpy H tot , which is the total heat released by the curing reaction. If the instantaneous reaction enthalpy is defined as H(t, T), hence the curing process variable ϕ is defined as given in Equation (2). This parameter ranges from ϕ = 0 when curing is not initiated, and ϕ = 1 when the cured state is fully achieved.
The model proposed by Kamal [12] is used to describe the evolution of the curing state ϕ(t) with a similar approach to the one used by Danzi et al. [10]. The curing kinetics is given in Equation (3): The temperature dependent functions k 1 (T) and k 2 (T) are given in Equations (4) and (5). They can represent different reaction processes such as a cure for k 1 (T), and a post-cure for k 2 (T).
where: The parameters of the curing kinetics model are calibrated from differential scanning calorimetry tests.
On the other hand, it is also necessary to relate the mechanical properties evolution with the curing state. A phenomenological approach [37] is proposed to link the evolution of the mechanical properties with the curing state as given in Equation (6). The shift function S(ϕ) ranges from S(0) = 0 when no curing is achieved and, S(1) = 1 when the cure is completely achieved. This equation uses an arctan function because it approximates the actual experimental data shape and fulfills the boundary limit values (0 to 1).

S(ϕ) =
arctan(β(ϕ − γ)) + arctan(β·γ) arctan(β·γ) + arctan(β(1 − γ)) The parameters β and γ can be obtained by fitting experimental data. β represents the slope and γ a phase shift to account for the delay of the property evolution with the actual curing state. The experimental data is obtained by DMA, but it can also be obtained by tensile tests at different curing states. Therefore, the instantaneous elastic modulus E(ϕ) is obtained by Equation (7), the Poisson's ratio by Equation (8) and the chemical shrinkage by Equation (9). where: v(0) = 0.5 Poisson's ratio at un-cured state v(1) = v m Poisson's ratio at cured state E(0) = 0 Elastic modulus at un-cured state E(1) = E m Elastic modulus at cured state α(0) = 0 Chemical shrinkage at un-cured state α(1) = α m Chemical shrinkage at cured state. The properties at completely cured state v m , E m , α m are the nominal properties of the cured material at a given temperature. To avoid numerical issues, it is considered that The initial elastic behavior is defined by a linear isotropic relation between the stress tensor, σ, and the elastic strain, ε e as given in Equation (10): where D e is the isotropic elastic tensor typically defined by the Young's modulus E, and the Poisson's ratio v.
A paraboloidal yield criterion Φ, Equation (11), is used to account for hydrostatic pressure effects: where S y− and S y+ are the compressive and tensile yield strengths of the material, respectively, J 2 is the second invariant of the deviatoric stress tensor, and I 1 is the first invariant of the stress tensor. The flow rule ψ is given in Equation (12). It is non-associative to account for the volumetric deformation: ψ(σ) = 3·J 2 + α/9·(I 1 ) 2 (12) where α corresponds to the material parameter that adjusts the volumetric component of the plastic flow. It is dependent of the plastic Poisson's ratio, and it is given by: The yield strengths are used to define the yield surface. Hardening is considered to affect both yield strengths and the increment of equivalent plastic strain is: where k ensures that the equivalent plastic strain is the same to that obtained in the uniaxial test: The yield stress hardening curve is approximated with a potential law as: where Sy 0 is the initial yielding stress, K s is the hardening coefficient and n c is the hardening exponent. One curve is used for the tensile part, and another curve is used for the compressive case.
The damage evolution model proposed by Melro et al. [27] enforces the damage to develop only in the material elastic modulus instead of the complete elastic tensor, an approach commonly followed by other continuum damage mechanics models usually implemented in a finite element method context [38,39].
The thermodynamic free energy potential D m is given in Equation (17): where σ corresponds to the stress tensor, E m and v m correspond to the undamaged material elastic modulus and Poisson's ratio, respectively, and d m corresponds to the damage variable. It is assumed that the damage process is de-coupled from the plastic dissipation process. Therefore, the plastic dissipation D P m is not considered in the subsequent damage formulation.
The damage initiation criterion Φ d m is given by Equation (18). This equation has a similar shape compared to the yield criterion, with the difference that S u− and S u+ correspond to compressive and tensile ultimate strengths, while I 1 and J 2 correspond to the effective un-damaged effective stress invariants in the material.
The damage evolution condition is given by Equation (19). The damage parameter r m ranges from 1, to satisfy the condition F d m < 0 when no damage is present, to +∞ when the material is completely damaged.
An exponential damage law is proposed to relate the damage parameter r m in the range [1, +∞] with the damage variable d m in the range [0, 1]. Equation (20) presents the exponential law where A m is a calibration parameter to regularize the volumetric damage dissipated energy in the element with the fracture energy.
The damage dissipated energy Y m can be obtained differentiating Equation (17), as shown in Equation (21).
The total dissipated energy due to the fracture process W f m can be calculated integrating the dissipated energy over the internal damage parameter as given in Equation (22). The total dissipated energy per unit volume should be equal to the surface fracture energy G f lm divided by the characteristic element size l e to ensure a mesh-independent damage law. The total strain mechanical strain ε t cam be decomposed as shown in Equation (23) where ε e is the elastic strain, ε p is the plastic strain, ε th is thermal expansion induced strains and ε sh are the chemical contraction-induced strains.
To couple the curing phenomena in the material with the mechanical response, the elastic constitutive tensor D e (ϕ, T) is assumed to be a function of the curing level and the temperature as shown in Equation (24), where S(ϕ) is the shift function and β(T) is the temperature dependence function of the elastic modulus.
The temperature dependence function β(T) is given in Equation (25) and it is a function of the coefficient β E for a given reference temperature T * . It was, originally proposed by Bai et al. [37]. β E is extracted by fitting the variation of the storage modulus with temperature obtained from a DMA test on a sample of cured resin. The obtained value for β E = 1.98.
Due to the nature of the curing problem, the equilibrium stress rate should be written in rate form similar to an hypoelastic law T , because the elastic constitutive tensor changes with time and temperature [10,40]. With this approach, material hardening due to curing and softening due to a temperature increase is appropriately considered, even at constant mechanical strain.
The instantaneous stress is approximated by σ = D e ·ε e , to be differentiated with temperature. In this procedure, the elastic strain is assumed to be constant during temperature variation ∂ε e /∂T = 0. For . T < 0, ∂D e /∂T > 0 which implies a stress increase for fixed elastic strain. This phenomenon is unrealistic, hence, for . T < 0, ∂D e /∂T = 0: Finally, during the curing process the elastic strain measure cannot be directly related with the initial/undeformed configuration. Hence, an instantaneous elastic strain measure can be calculated using the compliance tensor (inverse of elastic tensor) for the current stress state. With this approach the undeformed configuration is kept as the stress-free configuration which is not the initial configuration. This formulation is compatible with the plasticity theory but without damage. At damage onset, plasticity and curing are not allowed to develop [6].
The stress increment is calculated using a trapezoidal integration scheme as shown in Equation (29). ∆T n+1 means only positive changes of temperature.

Finite Element Model
The RVE geometry is defined using a constant fiber diameter d f and a square transverse section. Setting the number of fibers n f = n 2 f (to achieve a square RVE) and fiber volume fraction V f , the side size L d of the RVE can be explicitly calculated with Equation (30). Next, to generate the random distribution of fibers inside the RVE, the methodology proposed by Calatanotti [29] is implemented. Figure 1 shows an example of a typical RVE with n f = 25 and V f = 60%, where L a corresponds to the length of the RVE in the fiber direction.
Additional constraints for the RVE generation include a maximum fiber separation of 4 to avoid unrealistic resin-rich regions and a minimum fiber separation of 1.1 to avoid meshing difficulties. Lower fiber separations can also be achieved if required, but considerations regarding a good discretization of the inter-fiber regions led to the adoption of this value.
PBCs are implemented to guarantee continuity in the displacement field locally without over-constraining the boundary faces. The PBCs are derived from the strain-displacement definition, applied to opposite faces as given in Equation (31). The numerical implementation follows the procedure proposed by Barbero [30], where the nodal constrain equations are generated for each set of opposite faces, considering that edges and corners must be joined to guarantee the equations consistency. The general form of the constrain equations are given in Equation (32), where corresponds to the displacement in the direction at face , while corresponds to the displacement in the same direction and opposite face.
is the length of the RVE in the Additional constraints for the RVE generation include a maximum fiber separation of 4 × d f to avoid unrealistic resin-rich regions and a minimum fiber separation of 1.1 × d f to avoid meshing difficulties. Lower fiber separations can also be achieved if required, but considerations regarding a good discretization of the inter-fiber regions led to the adoption of this value.
PBCs are implemented to guarantee continuity in the displacement field locally without over-constraining the boundary faces. The PBCs are derived from the straindisplacement definition, applied to opposite faces as given in Equation (31).
The numerical implementation follows the procedure proposed by Barbero [30], where the nodal constrain equations are generated for each set of opposite faces, considering that edges and corners must be joined to guarantee the equations consistency. The general form of the constrain equations are given in Equation (32), where {u i } + f corresponds to the displacement in the direction i at face f , while {u i } − f corresponds to the displacement in the same direction and opposite face. L f is the length of the RVE in the direction f and ε i f is the imposed strain component. The loading in the RVE is imposed via a master node.
The effective composite properties, incl. elastic tensor components, yield strength and ultimate strength, are identified from the force-displacement response. The material parameters are calculated by enforcing that the load-displacement behavior of an equivalent monolithic RVE matches the full model RVE. Instead of volume average strain or stress measures, the effective strain corresponds to ε i f from Equation (32), and the force conjugate is obtained from the force reaction at the master node. The effective stress is obtained by dividing this force reaction by the associated cross-section area of the RVE. Table 3 shows the test loading cases required to identify the selected material parameters during the current study. The effective elastic tensor is approximated with a transversely isotropic tensor, that is defined by five elastic constants, namely E L , E T , v L , G L and G T . Additionally, yielding stresses and failure stresses were also obtained for this loading conditions to show the influence of the manufacturing parameters on the effective inelastic deformation and strength properties of the composite. In the current study, attention will be given to the transverse properties, which are dominated by the matrix response. Hence, the results for the longitudinal modulus E L will not be included in the analysis. Table 3. Loading cases to identify the effective composite properties.

Transverse shear
Longitudinal shear

Transverse shear
Longitudinal shear

Transverse shear
Longitudinal shear

Transverse shear
Longitudinal shear Although v T is not required to identify the elastic tensor, it can be easily determined.
The solution procedure to analyze the cure residual stresses begins with a thermomechanical analysis where the temperature curing cycle is imposed in the RVE domain. A pure mechanical analysis is considered because the temperature gradients inside the RVE are negligible. During this analysis step, the curing is simulated to determine the elastic properties evolution. The micro-residual internal stresses develop as a consequence of the interaction between the resin and the fibers. The PBCs are imposed in this step with a traction-free condition to allow the RVE to expand and shrink freely but conserving the periodicity. The constitutive model is completely implemented in a UMAT user subroutine, the shrinkage and thermal expansion via a UEXPAN user subroutine and the PBCs are implemented using a Python script to generate the input file. Figure 2 shows a flowchart of the solution procedure applied at each element integration point. The characteristic element size is given directly by the Abaqus solver in the UMAT subroutine, and it corresponds to the cubic root of the element volume.
After the curing step is completed, another analysis step is made to perform the test loading cases given in Table 3. The PBCs are kept the same but in this step a displacementcontrolled load is applied in accordance, keeping the remaining degrees of freedom tractionfree. To improve the computational efficiency, a restart point technique is used at the end of the curing step, to be used as the starting point of each loading case.
After the curing step is completed, another analysis step is made to perform the test loading cases given in Table 3. The PBCs are kept the same but in this step a displacementcontrolled load is applied in accordance, keeping the remaining degrees of freedom traction-free. To improve the computational efficiency, a restart point technique is used at the end of the curing step, to be used as the starting point of each loading case.

Micro-Residual Stress Analysis
The nominal curing cycle recommended by the material supplier, given in Table 4, is used to manufacture the coupons and to perform the micromechanical simulations. A set of three different statistically representative RVE geometries are used to perform all simulations, reducing the geometry dependence effects. Figure 3 shows the selected RVE geometries, namely RVE1, RVE2 and RVE3. After a mesh convergence analysis, the element size was set to be approximately 1/20 the fiber diameter for the matrix, and 1/10 for the fibers (to reduce computational costs). The average element size in the matrix is 0.00035 mm, which is defined in a compromise with the minimum fiber separation 1.1 to ensure at least two elements in the smaller gaps and avoid spurious artificial stresses due to badly shaped elements. Linear solid elements C3D8 from the Abaqus library were used. The matrix and the fiber are connected by tie constraints.
The fact that interface failure is not considered in these analyses can lead to overestimations in the effective strength prediction of the material, because the interface strength is typically lower than the bulk resin strength [28,42]. Although the current mod-

Micro-Residual Stress Analysis
The nominal curing cycle recommended by the material supplier, given in Table 4, is used to manufacture the coupons and to perform the micromechanical simulations. After a preliminary RVE size analysis considering RVEs with 9, 16, 25 and 36 fibers, an RVE with n f = 16 gives the best compromise between computational cost and convergence of the results convergence. This definition is consistent with previous computational micromechanics works [10,28], leading to an RVE size L d ≈ 4.58·d f . The RVE thickness was defined to guarantee at least 3 finite elements along this direction, i.e., L a ≈ 0.15·d f , because the RVE thickness has a minor effect on the transverse properties [41], the ones of interest in this study.
A set of three different statistically representative RVE geometries are used to perform all simulations, reducing the geometry dependence effects. Figure 3 shows the selected RVE geometries, namely RVE1, RVE2 and RVE3. After a mesh convergence analysis, the element size was set to be approximately 1/20 the fiber diameter for the matrix, and 1/10 for the fibers (to reduce computational costs). The average element size in the matrix is 0.00035 mm, which is defined in a compromise with the minimum fiber separation 1.1d f to ensure at least two elements in the smaller gaps and avoid spurious artificial stresses due to badly shaped elements. Linear solid elements C3D8 from the Abaqus library were used. The matrix and the fiber are connected by tie constraints.
The fact that interface failure is not considered in these analyses can lead to overestimations in the effective strength prediction of the material, because the interface strength is typically lower than the bulk resin strength [28,42]. Although the current modelling approach does not impede the consideration of interface failure, e.g., through cohesive interfaces, since no information regarding the interface properties for the current material system is available, and the evolution of the interface properties with curing state is, to the best of the authors' knowledge, unknown, it was decided to focus the analysis on the evaluation of the micro-residual stresses developed during the curing process and their effects on matrix failure only. elling approach does not impede the consideration of interface failure, e.g., through cohesive interfaces, since no information regarding the interface properties for the current material system is available, and the evolution of the interface properties with curing state is, to the best of the authors' knowledge, unknown, it was decided to focus the analysis on the evaluation of the micro-residual stresses developed during the curing process and their effects on matrix failure only. A sensitivity analysis was performed to analyze the influence of the chemical shrinkage, thermal expansion coefficient, elastic modulus and cure temperature in the microresidual stress state and mechanical properties. Table 5 shows the parameters used in this sensitivity analysis. Each parameter is arbitrarily ranged around the nominal value of the resin system.

Experimental Test Procedures
In order to calibrate the curing kinetics model, differential scanning calorimetry (DSC) tests were performed to un-cured resin samples following ASTM D3418. Next, to measure the instantaneous evolution of the storage modulus with curing, three DMA tests to un-cured resin samples are made, keeping the curing temperature fixed at 120 °C.
Additionally, coupon level tests were performed to measure the macro-mechanical properties of the corresponding composite system. The following tests were performed: • Transverse tensile tests following ASTM D3039.
The coupons were extracted from composite plates manufactured following the curing cycle presented in Table 4 and trimmed by CNC machining. A sensitivity analysis was performed to analyze the influence of the chemical shrinkage, thermal expansion coefficient, elastic modulus and cure temperature in the microresidual stress state and mechanical properties. Table 5 shows the parameters used in this sensitivity analysis. Each parameter is arbitrarily ranged around the nominal value of the resin system.

Experimental Test Procedures
In order to calibrate the curing kinetics model, differential scanning calorimetry (DSC) tests were performed to un-cured resin samples following ASTM D3418. Next, to measure the instantaneous evolution of the storage modulus with curing, three DMA tests to uncured resin samples are made, keeping the curing temperature fixed at 120 • C.
Additionally, coupon level tests were performed to measure the macro-mechanical properties of the corresponding composite system. The following tests were performed: • Transverse tensile tests following ASTM D3039.
The coupons were extracted from composite plates manufactured following the curing cycle presented in Table 4 and trimmed by CNC machining.

Experimental Results
The curing model parameters for the epoxy resin used in this study were obtained by fitting the experimental data of the curing process measured with DSC tests. The parameters are given in Table 6.  Figure 4a shows the results of one representative DMA curing test. It shows how the elastic modulus begins to grow after the temperature of 100 • C is reached while the curing level develops.
(1/s) (kJ/mol) (1/s) (kJ/mol) 3.24 × 10 15 134,627 0 0 1.00 0 Figure 4a shows the results of one representative DMA curing test. It shows how the elastic modulus begins to grow after the temperature of 100 °C is reached while the curing level develops.
Using the maximum obtained storage modulus and Equation (7) to obtain the instantaneous values of the shift function , Figure 4b shows the evolution of the shift function and the curing state of the material with time. It shows the storage modulus follows the curing state response with some delay.  Accordingly, the shift function given in Equation (7) is fitted, obtaining the parameters β = 4 and γ = 0.8. A comparison between the experimental relation and the relation from Equation (7) is given in Figure 5a. Figure 5b shows the estimated mechanical properties of the resin system using Equations (7)-(9), which will be implemented in the numerical material models.  Using the maximum obtained storage modulus and Equation (7) to obtain the instantaneous values of the shift function S(ϕ), Figure 4b shows the evolution of the shift function S(ϕ) and the curing state of the material with time. It shows the storage modulus follows the curing state response with some delay.
Accordingly, the shift function given in Equation (7) is fitted, obtaining the parameters β = 4 and γ = 0.8. A comparison between the experimental relation and the relation from Equation (7) is given in Figure 5a. Figure 5b shows the estimated mechanical properties of the resin system using Equations (7)-(9), which will be implemented in the numerical material models. rameters are given in Table 6.  15 134,627 0 0 1.00 0 Figure 4a shows the results of one representative DMA curing test. It shows how the elastic modulus begins to grow after the temperature of 100 °C is reached while the curing level develops.
Using the maximum obtained storage modulus and Equation (7) to obtain the instantaneous values of the shift function , Figure 4b shows the evolution of the shift function and the curing state of the material with time. It shows the storage modulus follows the curing state response with some delay.  Accordingly, the shift function given in Equation (7) is fitted, obtaining the parameters β = 4 and γ = 0.8. A comparison between the experimental relation and the relation from Equation (7) is given in Figure 5a. Figure 5b shows the estimated mechanical properties of the resin system using Equations (7)-(9), which will be implemented in the numerical material models. The tests coupons after mechanical testing are shown in Figure 6. Two transverse tensile coupons failed near the tabs; however, the measured strength has a low deviation, and the values are similar to the ones reported by the manufacturer and other test campaigns, making them suitable for the scope of this work. The failure modes for the longitudinal and transverse shear tests are within the expected behavior. level.
The tests coupons after mechanical testing are shown in Figure 6. Two transverse tensile coupons failed near the tabs; however, the measured strength has a low deviation, and the values are similar to the ones reported by the manufacturer and other test campaigns, making them suitable for the scope of this work. The failure modes for the longitudinal and transverse shear tests are within the expected behavior. A summary of the experimental results is given in Table 7

Residual Stress Analysis
During the curing process, the elastic modulus is developed at the same time the material expands due to the effects of the thermal expansion, and contracts due to chemical shrinkage. If the matrix material is completely free to deform it is not possible to generate stresses at the micro-mechanical level. However, in the case of fiber-reinforced composites, the fibers introduce local restrictions to the free deformation of the resin leading to the generation of local micro-residual stresses, the first stress generating mechanism. The second stress generating mechanism is the CTE incompatibility between the two constituents that becomes more relevant during the cooling-down stage.
To analyze the generated micro-residual stresses, Figures 7 and 8 show the stress distribution for the three different RVEs by means of the von Mises stress and the pressure stress. These measures were selected because they correspond the stress invariants that define the yielding and the failure criteria given in Equations (11) and (18). A summary of the experimental results is given in Table 7.

Residual Stress Analysis
During the curing process, the elastic modulus is developed at the same time the material expands due to the effects of the thermal expansion, and contracts due to chemical shrinkage. If the matrix material is completely free to deform it is not possible to generate stresses at the micro-mechanical level. However, in the case of fiber-reinforced composites, the fibers introduce local restrictions to the free deformation of the resin leading to the generation of local micro-residual stresses, the first stress generating mechanism. The second stress generating mechanism is the CTE incompatibility between the two constituents that becomes more relevant during the cooling-down stage.
To analyze the generated micro-residual stresses, Figures 7 and 8 show the stress distribution for the three different RVEs by means of the von Mises stress and the pressure stress. These measures were selected because they correspond the stress invariants that define the yielding and the failure criteria given in Equations (11) and (18).
The regions between the fibers concentrate higher deviatoric (von Mises) stresses in the matrix, being the most favorable regions to damage initiation and propagation. On the other hand, the hydrostatic stress distribution in the matrix is predominantly negative, which means that the bulk is in tensile state (remembering that σ Pressure = −1/3 trace(σ)), while the fibers remain in compression, as expected since the whole RVE is in global equilibrium and under traction-free boundary conditions. The regions between the fibers concentrate higher deviatoric (von Mises) stresses in the matrix, being the most favorable regions to damage initiation and propagation. On the other hand, the hydrostatic stress distribution in the matrix is predominantly negative, which means that the bulk is in tensile state (remembering that = −1/ 3 trace ), while the fibers remain in compression, as expected since the whole RVE is in global equilibrium and under traction-free boundary conditions.
To investigate the stress components evolution during the curing process, Figure 9 shows the von Mises stress and the hydrostatic stress components evolution as a function of time, and Figure 10 shows the effective RVE strain measured at the master nodes. Temperature is also superposed to ease the analysis. In this analysis, the stress measures correspond to a volumetric average over the matrix volume region. During the heating stage prior to the initiation of the curing in the polymer, no stress is generated and the RVE expands due to the combined effect of the positive CTE of the fibers and resin in the transverse direction, while it contracts in the longitudinal direction due to the negative CTE of the fibers.
Next, the curing process begins, and the matrix shrinks at a constant temperature (inside the plateau), until the curing process stops, achieving a stress equilibrium state with a smaller volume. A cooling stage follows, where additional straining is retrieved, developing more micro-residual stresses. The regions between the fibers concentrate higher deviatoric (von Mises) stresses in the matrix, being the most favorable regions to damage initiation and propagation. On the other hand, the hydrostatic stress distribution in the matrix is predominantly negative, which means that the bulk is in tensile state (remembering that = −1/ 3 trace ), while the fibers remain in compression, as expected since the whole RVE is in global equilibrium and under traction-free boundary conditions.
To investigate the stress components evolution during the curing process, Figure 9 shows the von Mises stress and the hydrostatic stress components evolution as a function of time, and Figure 10 shows the effective RVE strain measured at the master nodes. Temperature is also superposed to ease the analysis. In this analysis, the stress measures correspond to a volumetric average over the matrix volume region. During the heating stage prior to the initiation of the curing in the polymer, no stress is generated and the RVE expands due to the combined effect of the positive CTE of the fibers and resin in the transverse direction, while it contracts in the longitudinal direction due to the negative CTE of the fibers.
Next, the curing process begins, and the matrix shrinks at a constant temperature (inside the plateau), until the curing process stops, achieving a stress equilibrium state with a smaller volume. A cooling stage follows, where additional straining is retrieved, developing more micro-residual stresses. To investigate the stress components evolution during the curing process, Figure 9 shows the von Mises stress and the hydrostatic stress components evolution as a function of time, and Figure 10 shows the effective RVE strain measured at the master nodes. Temperature is also superposed to ease the analysis. In this analysis, the stress measures correspond to a volumetric average over the matrix volume region. During the heating stage prior to the initiation of the curing in the polymer, no stress is generated and the RVE expands due to the combined effect of the positive CTE of the fibers and resin in the transverse direction, while it contracts in the longitudinal direction due to the negative CTE of the fibers. The strain measure evolution with time and the degree of cure present trends similar to the strain measured experimentally with optical fibers by Minakushi [21]. Therefore, the matrix volumetric shrinkage is the precursor of the stress development during curing because it occurs at constant temperature, while the thermal contraction (CTE) is the stress precursor during cooling. The longitudinal strain deformation of the RVE is very low because it is dominated by the high stiffness of the fibers and their low negative CTE.   The strain measure evolution with time and the degree of cure present trends similar to the strain measured experimentally with optical fibers by Minakushi [21]. Therefore, the matrix volumetric shrinkage is the precursor of the stress development during curing because it occurs at constant temperature, while the thermal contraction (CTE) is the stress precursor during cooling. The longitudinal strain deformation of the RVE is very low because it is dominated by the high stiffness of the fibers and their low negative CTE. The chemical shrinkage for this resin is 2%, while the thermal contraction during cooling is only 0.6% (CTE multiplied by the temperature amplitude-from curing temperature to room temperature). However, the final volume reduction is around 0.9% estimated from the final volume of the RVE as presented Figure 11.   Next, the curing process begins, and the matrix shrinks at a constant temperature (inside the plateau), until the curing process stops, achieving a stress equilibrium state with a smaller volume. A cooling stage follows, where additional straining is retrieved, developing more micro-residual stresses.
The strain measure evolution with time and the degree of cure present trends similar to the strain measured experimentally with optical fibers by Minakushi [21]. Therefore, the matrix volumetric shrinkage is the precursor of the stress development during curing because it occurs at constant temperature, while the thermal contraction (CTE) is the stress precursor during cooling. The longitudinal strain deformation of the RVE is very low because it is dominated by the high stiffness of the fibers and their low negative CTE.
The chemical shrinkage for this resin is 2%, while the thermal contraction during cooling is only 0.6% (CTE multiplied by the temperature amplitude-from curing temperature to room temperature). However, the final volume reduction is around 0.9% estimated from the final volume of the RVE as presented Figure 11. If a full restrained analysis is performed, the volume contraction is enough to cause material failure, because the failure strain is lower than 0.7%. For the current example, no damage is retrieved from the numerical results, and only small plastic straining is found in the stress concentration zones.
Other results that can be extracted from the curing analysis are the failure index measurements in the RVE given by Equations (11) and (18). For this purpose, Figure 12 shows the maximum failure index retrieved from the RVE as a function of time. It is clear from the previous figures that the micro-residual stresses achieve representative values especially between the narrowest fiber gaps and these stresses lead to the development of plastic strains. For the current example, no damage is achieved but the failure index is If a full restrained analysis is performed, the volume contraction is enough to cause material failure, because the failure strain is lower than 0.7%. For the current example, no damage is retrieved from the numerical results, and only small plastic straining is found in the stress concentration zones.
Other results that can be extracted from the curing analysis are the failure index measurements in the RVE given by Equations (11) and (18). For this purpose, Figure 12 shows the maximum failure index retrieved from the RVE as a function of time. It is clear from the previous figures that the micro-residual stresses achieve representative values especially between the narrowest fiber gaps and these stresses lead to the development of plastic strains. For the current example, no damage is achieved but the failure index is above 0.9.
If a full restrained analysis is performed, the volume contraction is enough to cause material failure, because the failure strain is lower than 0.7%. For the current example, no damage is retrieved from the numerical results, and only small plastic straining is found in the stress concentration zones.
Other results that can be extracted from the curing analysis are the failure index measurements in the RVE given by Equations (11) and (18). For this purpose, Figure 12 shows the maximum failure index retrieved from the RVE as a function of time. It is clear from the previous figures that the micro-residual stresses achieve representative values especially between the narrowest fiber gaps and these stresses lead to the development of plastic strains. For the current example, no damage is achieved but the failure index is above 0.9.

Effective Mechanical Properties
To study the influence of the micro-residual stresses in the mechanical performance of the composite material, several simulations were performed using the selected RVEs subjected to different loading conditions as planned in Table 3, giving special attention to the properties dominated by the matrix behavior. Two sets of simulations are analyzed: the first corresponds to an analysis with the nominal material properties without considering the curing analysis and without the micro-residual stress effects, which is named "no-cure", and another set that is analyzed using the micro-residual stresses and the actual material condition coming from the curing step analysis, named "cure".

Effective Mechanical Properties
To study the influence of the micro-residual stresses in the mechanical performance of the composite material, several simulations were performed using the selected RVEs subjected to different loading conditions as planned in Table 3, giving special attention to the properties dominated by the matrix behavior. Two sets of simulations are analyzed: the first corresponds to an analysis with the nominal material properties without considering the curing analysis and without the micro-residual stress effects, which is named "no-cure", and another set that is analyzed using the micro-residual stresses and the actual material condition coming from the curing step analysis, named "cure".
The stress-strain response for the transverse tensile loading case is presented in Figure 13. A direct comparison between the three different RVEs for the "no-cure" and "cure" condition is made. First, the overall response of the RVE is linear until failure, although the matrix resin is subjected to some plastic deformation. The elastic modulus remains unaffected because the material is completely cured, but the transverse tensile strength is reduced by the micro-residual stresses as expected. The stress-strain response for the transverse tensile loading case is presented in Figure 13. A direct comparison between the three different RVEs for the "no-cure" and "cure" condition is made. First, the overall response of the RVE is linear until failure, although the matrix resin is subjected to some plastic deformation. The elastic modulus remains unaffected because the material is completely cured, but the transverse tensile strength is reduced by the micro-residual stresses as expected. A similar response is retrieved for the shear loading simulations in both transverse and longitudinal directions, as shown in Figure 14. The shear modulus is practically unaffected by the micro-residual stresses and only the shear strength is reduced.  RVE1 Cure S12 RVE2 Cure S12 RVE3 Cure S12 RVE1 No-cure S12 Figure 13. Transverse tensile stress-strain response comparison.
A similar response is retrieved for the shear loading simulations in both transverse and longitudinal directions, as shown in Figure 14. The shear modulus is practically unaffected by the micro-residual stresses and only the shear strength is reduced. A similar response is retrieved for the shear loading simulations in both transverse and longitudinal directions, as shown in Figure 14. The shear modulus is practically unaffected by the micro-residual stresses and only the shear strength is reduced.  Figure 15 shows the fracture pattern for each loading case for the "no-cure" and "cure" tests. A similar failure pattern is retrieved for both cases because it is given by the loading condition, with the only difference in the initiation point that changes the position of the crack bands. In the case of the cure analysis, the initiation point is influenced by the micro-residual stress levels, unlike the no-cure analysis.

Stress (MPa)
Strain (mm/mm) RVE1 Cure S12 RVE2 Cure S12 RVE3 Cure S12 RVE1 No-cure S12 RVE2 No-cure S12 RVE3 No-Cure S12  Figure 15 shows the fracture pattern for each loading case for the "no-cure" and "cure" tests. A similar failure pattern is retrieved for both cases because it is given by the loading condition, with the only difference in the initiation point that changes the position of the crack bands. In the case of the cure analysis, the initiation point is influenced by the micro-residual stress levels, unlike the no-cure analysis. To facilitate the comparison of results between both analysis sets, Table 8 shows a summary of the predicted effective mechanical properties and the experimental results. From the tensile tests, the transverse elastic modulus and the Poisson's ratio from the experimental measurements are similar to the numerical predictions, the transverse tensile strength presents a reduction of 11% from the "no-cure" condition although the experimental measured value is slightly lower. This difference can be attributed to the effect of the fiber-matrix interface and other manufacturing defects which are not considered in the current analysis [28,42]. Another relevant result that should be highlighted is the fact that the tensile strength of the composite is almost 30% lower than the bulk matrix strength, a tendency that is commonly highlighted in the literature [9,10].

Numerical Predictions
Experiments vs. To facilitate the comparison of results between both analysis sets, Table 8 shows a summary of the predicted effective mechanical properties and the experimental results. From the tensile tests, the transverse elastic modulus E 2 and the Poisson's ratio v 21 from the experimental measurements are similar to the numerical predictions, the transverse tensile strength S +UT presents a reduction of 11% from the "no-cure" condition although the experimental measured value is slightly lower. This difference can be attributed to the effect of the fiber-matrix interface and other manufacturing defects which are not considered in the current analysis [28,42]. Another relevant result that should be highlighted is the fact that the tensile strength of the composite is almost 30% lower than the bulk matrix strength, a tendency that is commonly highlighted in the literature [9,10]. A similar trend is observed for the predicted transverse and longitudinal shear moduli, which are similar to the experimental results, while the strength measures present a reduction of 22% and 12%, respectively, when compared with the "no-cure" condition. Similarly, the experimental measured values are lower than the numerical predictions. The effective properties estimated considering the cure condition are closer to the experimental results.
These results show evidence that the micro-residual stresses in the RVE after the curing process can reduce the material strength and should be considered when performing homogenization procedures. Regarding the fracture response, further work is required to understand the damage mechanisms and the appropriate size effects that must be considered, because the application of the material models and mechanical properties from the macro-scale experimental measurements are straightforward for the elastic and strength parameters, but not for the fracture properties.

Sensitivity Analysis Results
A sensitivity analysis of the relevant material and process parameters that influence the curing micro-residual stresses is presented in this section. After a preliminary study following previous literature results [6], the material parameters that have a direct influence in the residual stresses are the CTE, the chemical shrinkage and the matrix elastic modulus. On the other hand, from the process point of view, since thermal gradients and transient effects are negligible at the RVE scale, the only relevant parameter is the cure temperature because it has a direct influence in the final curing level. Figure 16 shows the results for the sensitivity analyses for each of the selected variables regarding the micro-residual stress level and the failure index. The failure index corresponds to the maximum value of the failure initiation criteria, given by Equation (18), over the matrix region. The results presented for the stress measure corresponds to the volumetric average over the entire matrix region of the RVE.
A clear tendency to increase the micro-residual stresses arises from increasing the CTE and the chemical shrinkage. For higher levels of micro-residual stresses, damage is achieved, with the failure index becoming higher than unity.
The elastic modulus also has a relevant influence on the micro-residual stresses. The tendency shows that micro-residual stresses increase with increasing elastic modulus, including scenarios of damage onset. It is worth to note that the strength properties are kept constant during the sensitivity analyses. transient effects are negligible at the RVE scale, the only relevant parameter is the cure temperature because it has a direct influence in the final curing level. Figure 16 shows the results for the sensitivity analyses for each of the selected variables regarding the micro-residual stress level and the failure index. The failure index corresponds to the maximum value of the failure initiation criteria, given by Equation (18), over the matrix region. The results presented for the stress measure corresponds to the volumetric average over the entire matrix region of the RVE.  The cure temperature has a different influence on the micro-residual stresses because it changes the degree of cure evolution and the final curing level in the matrix and not the mechanical response directly. This difference in the curing level changes the final elastic modulus, Poisson's ratio and overall mechanical behavior of the material. For lower temperatures, keeping the cure cycle time unchanged, the cure level is reduced to 20%, a very low value, but (as known) the final properties would not be suitable for a composite material. Figure 17 shows the effective stress-strain relations, which present the expected tendency following the micro-residual stress analyses from Figure 16. The identification of the curves corresponds to Table 5; for example, CTE1 corresponds to the CTE value of case 1. The increase in micro-residual stresses generates a direct reduction of the material tensile strength.
The cure temperature reduction leads to a reduction in the elastic modulus followed by a loss of material strength (although resin strength parameters are kept fixed). The variation of the elastic modulus has a direct effect in the composite stiffness, as expected, and leads to a tensile strength reduction due to the higher micro-residual stress levels.
The sensitivity results show that the strength reduction can be mainly attributed to the micro-residual stress state that was characterized by the average von Mises stress and hydrostatic stress. Using all the data obtained above and plotting the tensile strength against the von Mises stress and the hydrostatic component, as shown in Figure 18, two results are retrieved. First, the von Mises stress follows a linear relation with the average hydrostatic component for the curing micro-residual stresses state. Second, independently of the parameter that is being modified, the micro-residual stress effects hold, allowing us to conclude, within the scope of the current work, that the loss of material performance (tensile strength) is a direct consequence of the micro-residual stress levels. case 1. The increase in micro-residual stresses generates a direct reduction of the material tensile strength.
The cure temperature reduction leads to a reduction in the elastic modulus followed by a loss of material strength (although resin strength parameters are kept fixed). The variation of the elastic modulus has a direct effect in the composite stiffness, as expected, and leads to a tensile strength reduction due to the higher micro-residual stress levels. The sensitivity results show that the strength reduction can be mainly attributed to the micro-residual stress state that was characterized by the average von Mises stress and hydrostatic stress. Using all the data obtained above and plotting the tensile strength against the von Mises stress and the hydrostatic component, as shown in Figure 18, two results are retrieved. First, the von Mises stress follows a linear relation with the average hydrostatic component for the curing micro-residual stresses state. Second, independently of the parameter that is being modified, the micro-residual stress effects hold, allowing us to conclude, within the scope of the current work, that the loss of material performance (tensile strength) is a direct consequence of the micro-residual stress levels. Although it is recognized that additional work is required to appropriately understand the influence of process parameters on the mechanical response of the matrix resin at the micromechanical level, the difficulties associated to experimental work at this scale make the numerical insight a valuable tool to understand better the effective macro-mechanical response of composite materials, as demonstrated by the methodology proposed in this work.

Conclusions
The micro-residual stress analysis showed that stress concentrations are retrieved between the fibers in the narrower gaps, while the higher hydrostatic components are present in the resin-rich regions. This effect is clearly attributed to the constraining effect that the fiber arrangement imposes on the matrix, because the whole RVE is in traction-free Although it is recognized that additional work is required to appropriately understand the influence of process parameters on the mechanical response of the matrix resin at the micromechanical level, the difficulties associated to experimental work at this scale make the numerical insight a valuable tool to understand better the effective macro-mechanical response of composite materials, as demonstrated by the methodology proposed in this work.

Conclusions
The micro-residual stress analysis showed that stress concentrations are retrieved between the fibers in the narrower gaps, while the higher hydrostatic components are present in the resin-rich regions. This effect is clearly attributed to the constraining effect that the fiber arrangement imposes on the matrix, because the whole RVE is in traction-free condition.
During the curing process, the whole RVE experiences volume changes that follow the superposition of the CTE effects during heating, chemical shrinkage during curing, and the final thermal contraction during cooling down. This superposition gives a final volume reduction for the current example of 0.9%, approximately half of the chemical shrinkage (2%).
The micro-residual stresses have almost no influence in the effective stiffness of the composite system, but they have a more considerable effect in the strength properties, for example, reducing the transverse tensile strength by 11%. However, the experimental measured strength properties with macro-scale coupons are still slightly lower than the predicted properties considering the curing conditions.
The sensitivity analysis shows that increasing the CTE, chemical shrinkage and elastic modulus contributes to an increase of the micro-residual stresses. The cure temperature influences the final degree of cure of the material, degrading the elastic modulus.
A direct analysis of the transverse tensile strength as a function of the micro-residual stresses shows a linear relation between the von Mises stress and the hydrostatic stress. Therefore, the loss of performance can be attributed to the micro-residual stresses.