Next Article in Journal
A New Approach Utilizing Aza-Michael Addition for Hydrolysis-Resistance Non-Ionic Waterborne Polyester
Next Article in Special Issue
Design and Performance Evaluation of Polymer Matrix Composite Helical Springs
Previous Article in Journal
Sensor Fusion for Simultaneous Estimation of In-Plane Permeability and Porosity of Fiber Reinforcement in Resin Transfer Molding
Previous Article in Special Issue
Eco-Friendly Hybrid PLLA/Chitosan/Trichoderma asperellum Nanomaterials as Biocontrol Dressings against Esca Disease in Grapevines
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Institute of Science and Innovation in Mechanical and Industrial Engineering (INEGI), Rua Dr. Roberto Frias, 4200-465 Porto, Portugal
2
Departamento de Engenharia Mecânica, Faculdade de Engenharia, Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal
*
Authors to whom correspondence should be addressed.
Polymers 2022, 14(13), 2653; https://doi.org/10.3390/polym14132653
Submission received: 30 May 2022 / Revised: 14 June 2022 / Accepted: 27 June 2022 / Published: 29 June 2022

Abstract

:
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.

Graphical Abstract

1. 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.
Two families of polymeric matrices are commonly used in FRPs, namely thermoplastics and thermosets. While thermoplastics are processed by melting the polymer, and cooling it down to its final configuration, thermoset-based composites are manufactured with a non-reversible thermo-curing process where the polymer is transformed from a gel to a monolithic material. The latter have been the preferred option in aerospace applications due to their easier processing characteristics.
Inadequate manufacturing processes in composite materials generate low-quality products and performance reduction in as-manufactured conditions due to residual stresses, voids formation, or incomplete curing. This performance reduction can lead to design failure because material properties were overestimated, or because the geometry does not fit the requirements [5,6]. Linking of composites manufacturing simulation with in-service performance is a recent research topic that has become affordable due to increased computing power, but we still lack a complete understanding of how it influences the final performance.
Focusing on carbon fiber-based thermoset composite systems (CFRPs), micro-residual stresses appear at the microstructural level due to the mismatch between the thermal expansion coefficients of the fiber and the polymeric matrix, and due to the chemical shrinkage, that takes place during the curing process [7,8,9,10,11]. Additional residual stresses can be obtained in composite laminates due to the anisotropy of the plies, but they are out of the scope of this work, which will focus on the micro- (constituents) scale of polymer matrix composites.
The curing process is inherently a thermo-chemical phenomenon for thermoset polymers where an exothermic chemical reaction starts at a given temperature and finishes when the polymerization process ends. One of the most commonly used approximation models for polymer curing was proposed by Kamal [12]. This model relates the curing state in the material with the time and temperature evolution. Currently, several works can be found in the literature regarding the experimental characterization of the curing behavior of polymeric resins, including epoxies [13,14,15,16,17,18,19,20]. Many of them use differential scanning calorimetry (DSC) and dynamical mechanical analysis (DMA) to identify the curing behavior and the evolution of the elastic modulus.
Experimental works to measure actual residual stresses at micro-scale level are not easy to perform due to the reduced length scale at which they occur. For instance, Minakushi [21] performs fiber-optic based tests to measure cure shrinkage in fiber-reinforced laminates, obtaining local strain measurements in the fiber direction and through the thickness, while Seers et al. [22] presents a recent review of the measurement techniques of residual stresses in composites, most of them relying on macro-measurements to obtain indirect measures of the micro-residual stresses. Similar difficulties arise regarding the characterization of the mechanical properties of the material at micro-scale level due to size effects [23], especially for the fracture toughness.
Numerical simulation of the properties evolution with curing is a current research field. To the best of the authors’ knowledge, Ersoy et al. [8] is among the first studies to deal with these simulations, considering only the elastic part of the material. Yuan et al. [24,25] uses viscoelastic constitutive modeling for the elastic modulus evolution during curing and representative volume elements (RVEs) consisting of one fiber or multiple fibers coupled with a multi-scale approach to capture temperature profiles from the macro-scale.
A similar approach was followed by Hui et al. [26], who also use a multi-scale approach to analyze the heat transfer effects of the curing process at the micro-scale level. They include an extended Drucker–Prager model to account for the polymer failure response. Danzi et al. [10] instead used a more advanced material constitutive model that accounts for the plastic response and failure of the polymer matrix, a model that was previously proposed by Melro et al. [27] and has been used by several authors to simulate the mechanical response of epoxies. All these works converge in the usage of RVE models [28,29] to analyze the mechanical response of the composite material at the micro-scale level with the appropriate boundary conditions (BCs) to guarantee solution consistency [30].
Another relevant topic relates to the methodologies used to estimate the effective macro-scale mechanical properties and the imposed BCs. Volumetric average measures of stress and strain fields are the most commonly used technique to extract the effective strain and stress response [28]. Second order homogenization theories are also available in the literature [31] but are not considered in this work because they require the usage of strain gradient measures, which falls outside the scope of this work.
Considering first order approximations, several authors have studied the influence of the BCs in the effective strain–stress response, concluding that the periodic boundary conditions (PBCs) give the best compromise between accuracy and computational simplicity [32,33]. However, the usage of PBCs has limitations in the simulation of strain localization phenomena such as damage propagation, because they over-constrain the strain field at the boundaries [28,32,34]. Nevertheless, PBCs provide accurate estimates up to failure initiation, which is the focus of this work.
Here, a micromechanical model is developed using the commercial software Abaqus® and user-material (UMAT) subroutines to simulate the curing process and estimate the micro-residual stresses and its influence on the mechanical performance of an aerospace grade epoxy resin. A distinction between residual (macro-scale) and micro-residual (micro-scale) stresses is made because the micro-residual stresses are developed inside the RVE with traction-free boundary conditions.
Using the methodology proposed by Danzi et al. [10] that couples the polymer constitutive formulation proposed by Melro et al. [27] with the curing kinetics model, an enhanced temperature dependence of the material properties, especially for the elastic properties, is included in this work to appropriately model the shrinkage and the cooling processes.
The effective strain and stress measures are extracted directly from the PBCs output, i.e., the strain is obtained from the strain of the master node, and the stress from the force-conjugate measure of the applied strain, instead of using volumetric averages that are not appropriate when strain localization is achieved. The predictions obtained are then compared with experimental data obtained from standardized composite macro-scale tests.
A detailed analysis of the RVE response during the curing process is also performed using three different RVE geometries. Instantaneous measures of the effective macro-strain are extracted to see effective material expansion and contraction, besides the typical stress analysis. This analysis shows how the constituents are subjected to local stresses while keeping the macro effective stress equal to zero.
Additionally, due to the uncertainties of the material parameters at micro-scale level, a sensitivity analysis is performed in this work, to assess the influence of the thermal expansion coefficient, chemical shrinkage, resin elastic modulus and cure temperature on the residual stresses and the posterior effective mechanical properties, to complement the overall analysis of the micro-residual stresses.
The following section describes the materials, the constitutive modeling approach, the finite element procedures and the micro-residual stress analysis methodology. Then, the results are analyzed to highlight the micro-residual stresses distribution inside the RVE, and to link their influence on the effective mechanical properties. Finally, conclusions are drawn based on the discussion of the principal findings of this study.

2. Materials and Methods

2.1. Materials Selection

The fiber-reinforcement system selected for the experiments and simulations consists on the Toray T700 standard modulus carbon fiber. The parameters used in the current analysis are taken from the datasheet [35] and other previous works [10,28,36]. These parameters are summarized in Table 1, where C T E L and C T E T are, respectively, the longitudinal and transverse coefficients of thermal expansion, E L and E T are, respectively, the longitudinal and transverse elastic moduli, v L is the longitudinal Poisson’s ratio, G L and G T are the longitudinal and transverse shear moduli.
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 C T E 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 l m 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 ).

2.2. Constitutive Models

2.2.1. Carbon Fibers

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 ) )
ε 11 ε 22 ε 33 ε 12 ε 13 ε 23 = 1 / E L v L / E L v L / E L 0 0 0 v L / E L 1 / E T v T / E T 0 0 0 v L / E L v T / E T 1 / E T 0 0 0 0 0 0 1 / 2 / G L 0 0 0 0 0 0 1 / 2 / G L 0 0 0 0 0 0 1 / 2 / G T σ 11 σ 22 σ 33 σ 12 σ 13 σ 23

2.2.2. 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 t o t , 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.
φ = H ( t , T ) H t o t
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):
d φ d t = ( k 1 ( T ) + k 2 ( T ) φ m ) ( 1 φ ) n
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 ) .
k 1 ( T ) = A 1 exp ( Δ E 1 R · T ( t ) )
k 2 ( T ) = A 2 exp ( Δ E 2 R · T ( t ) )
where:
  • A 1 ; A 2 : Reaction velocities [1/s].
  • Δ E 1 ; Δ E 2 : Activation energies [J/mol]
  • m ; n : Fitting exponents
  • R : Universal gas constant [8.31432 J/mol.K]
  • T ( t ) : Temperature
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).
E ( φ ) = E m · S ( φ )
v ( φ ) = S ( φ )   v c + ( 1 S ( φ ) ) v 0
α ( φ ) = α m · S ( φ )
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 E ( 0 ) = 10 3 · E m .
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):
σ = [ D e ] · ε e
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:
Φ ( σ , σ C , σ T ) = 6 · J 2 + 2 · I 1 · ( S y S y + ) 2 · S y · S y + = 0
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
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:
α = 9 2 ( 1 2 v p 1 + v p )
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:
Δ ε ¯ p = k Δ ε p : Δ ε p
where k ensures that the equivalent plastic strain is the same to that obtained in the uniaxial test:
k = 1 1 + 2 v p 2
The yield stress hardening curve is approximated with a potential law as:
S y / + = S y 0 + K s · ( ε ¯ p ) n c
where S y 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):
D m = σ 11 2 + σ 22 2 + σ 33 2 2 E m ( 1 d m ) v m E m ( σ 11 σ 22 + σ 22 σ 33 + σ 33 σ 11 ) + 1 + v m E m ( 1 d m ) ( σ 12 2 + σ 13 2 + σ 23 2 ) + D m P
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 m P is not considered in the subsequent damage formulation.
The damage initiation criterion Φ m d 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.
Φ m d = 3 J 2 ˜ S u · S u + + I 1 ˜ · ( S u S u + ) S u · S u +
The damage evolution condition is given by Equation (19). The damage parameter r m ranges from 1, to satisfy the condition F m d < 0 when no damage is present, to + when the material is completely damaged.
F m d = Φ m d r m 0
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.
d m = 1 exp ( A m ( 3 7 + 2 r m 2 ) ) 7 + 2 r m 2 2
The damage dissipated energy Y m can be obtained differentiating Equation (17), as shown in Equation (21).
Y m = D m d m = σ 11 2 + σ 22 2 + σ 33 2 2 E m ( 1 d m ) 2 + 1 + v m E m ( 1 d m ) 2 ( σ 12 2 + σ 13 2 + σ 23 2 )
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 l m divided by the characteristic element size l e to ensure a mesh-independent damage law.
W f m = 1 + Y m · d m r m d r m = G f l m l e

2.2.3. Thermo-Curing Coupling Strategy

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, ε t h is thermal expansion induced strains and ε s h are the chemical contraction-induced strains.
ε t = ε e + ε p + ε t h + ε s h
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.
D e ( φ , T ) = D e 0 · S ( φ ) · β ( T )
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 .
β ( T ) = ( 1 + β E · log ( T T * ) )
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.
σ ˙ = σ ε · ε ˙ + σ T · T ˙
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 :
σ T · T ˙ = ( D e T · ε e + D e · ε e T ) · T ˙ = ( D e T · ε e ) · T ˙
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].
ε e = ( D e ) 1 · σ
The stress increment is calculated using a trapezoidal integration scheme as shown in Equation (29). Δ T n + 1 means only positive changes of temperature.
σ n + 1 = σ n + 1 2 · ( D e n + D e n + 1 ) · Δ ε n + 1 + ( D e T · ε e ) · Δ T n + 1

2.3. 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 ¯ f 2 (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.
L d = n f · π · d f 2 4 · V f
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 strain-displacement definition, applied to opposite faces as given in Equation (31).
ε i j = d u i d x j = u i A u i B L j
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.
{ u i } f { u i } + f L f · ε i f = 0
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.
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 displacement-controlled 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.

2.4. 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 micro-mechanics 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.1 d 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 over-estimations 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.
A sensitivity analysis was performed to analyze the influence of the chemical shrinkage, thermal expansion coefficient, elastic modulus and cure temperature in the micro-residual 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.

2.5. 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.
  • Longitudinal shear tests following ASTM D3518.
  • Transverse shear tests following ASTM D5379.
The coupons were extracted from composite plates manufactured following the curing cycle presented in Table 4 and trimmed by CNC machining.

3. Results and Discussion

3.1. 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.
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.
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.

3.2. 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, Figure 7 and Figure 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 σ P r e s s u r e = 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 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 above 0.9.

3.3. 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.
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.
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 + U T 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.

3.4. 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.
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.
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.

4. 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.

Author Contributions

Conceptualization: P.T.G. and A.A., methodology: P.T.G. and A.A., software: P.T.G. and A.A., validation: P.T.G. and A.A., formal analysis: P.T.G., A.A. and N.R., writing—review and editing: P.T.G., A.A., N.R. and L.P., funding acquisition: N.R. and L.P., project administration: N.R. All authors have read and agreed to the published version of the manuscript.

Funding

The current research work has been funded by the project UIDP/50022/2020 New Generation of Cryogenic Propulsion Systems for Future Launchers from LAETA—Laboratório Associado de Energia, Transportes e Aeronáutica, with the support of FCT/MCTES and PIDDAC (Programa de Investimentos e Despesas de Desenvolvimento da Administração Central).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Barbero, J.E. Introduction to Composite Materials Design, 3rd ed.; CRC Press: Boca Raton, FL, USA, 2017; ISBN 9781315296494. [Google Scholar]
  2. Mouritz, A.P. Fibre–Polymer Composites for Aerospace Structures and Engines, 1st ed.; Woodhead Publishing Ltd.: Cambridge, UK, 2012; ISBN 978-1-85573-946-8. [Google Scholar]
  3. Li, C.; Yin, X.; Liu, Y.; Guo, R.; Xian, G. Long-term service evaluation of a pultruded carbon/glass hybrid rod exposed to elevated temperature, hydraulic pressure and fatigue load coupling. Int. J. Fatigue 2020, 134, 105480. [Google Scholar] [CrossRef]
  4. Bazli, M.; Abolfazli, M. Mechanical Properties of Fibre Reinforced Polymers under Elevated Temperatures: An Overview. Polymers 2020, 12, 2600. [Google Scholar] [CrossRef] [PubMed]
  5. Struzziero, G.; Teuwen, J.J.E.; Skordos, A.A. Numerical optimisation of thermoset composites manufacturing processes: A review. Compos. Part A Appl. Sci. Manuf. 2019, 124, 105499. [Google Scholar] [CrossRef]
  6. D’Mello, R.J.; Maiarù, M.; Waas, A.M. Virtual manufacturing of composite aerostructures. Aeronaut. J. 2016, 120, 61–81. [Google Scholar] [CrossRef] [Green Version]
  7. Baran, I.; Cinar, K.; Ersoy, N.; Akkerman, R.; Hattel, J.H. A Review on the Mechanical Modeling of Composite Manufacturing Processes. Arch. Comput. Methods Eng. 2017, 24, 365–395. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  8. Ersoy, N.; Garstka, T.; Potter, K.; Wisnom, M.R.; Porter, D.; Clegg, M.; Stringer, G. Development of the properties of a carbon fibre reinforced thermosetting composite through cure. Compos. Part A Appl. Sci. Manuf. 2010, 41, 401–409. [Google Scholar] [CrossRef]
  9. Lavoie, J.A.; Adolfsson, E. Stitch Cracks in Constraint Plies Adjacent to a Cracked Ply. J. Compos. Mater. 2001, 35, 2077–2097. [Google Scholar] [CrossRef]
  10. Danzi, F.; Fanteria, D.; Panettieri, E.; Mancino, M.C.C. A numerical micro-mechanical study on damage induced by the curing process in carbon/epoxy unidirectional material. Compos. Struct. 2019, 210, 755–766. [Google Scholar] [CrossRef]
  11. Danzi, F.; Fanteria, D.; Panettieri, E.; Palermo, M. A numerical micro-mechanical study of the influence of fiber–matrix interphase failure on carbon/epoxy material properties. Compos. Struct. 2017, 159, 625–635. [Google Scholar] [CrossRef] [Green Version]
  12. Kamal, M.R. Thermoset characterization for moldability analysis. Polym. Eng. Sci. 1974, 14, 231–239. [Google Scholar] [CrossRef]
  13. O’Brien, D.J.; Sottos, N.R.; White, S.R. Cure-dependent Viscoelastic Poisson’s Ratio of Epoxy. Exp. Mech. 2007, 47, 237–249. [Google Scholar] [CrossRef]
  14. Billotte, C.; Bernard, F.M.; Ruiz, E. Chemical shrinkage and thermomechanical characterization of an epoxy resin during cure by a novel in situ measurement method. Eur. Polym. J. 2013, 49, 3548–3560. [Google Scholar] [CrossRef]
  15. Nawab, Y.; Casari, P.; Boyard, N.; Jacquemin, F. Characterization of the cure shrinkage, reaction kinetics, bulk modulus and thermal conductivity of thermoset resin from a single experiment. J. Mater. Sci. 2013, 48, 2394–2403. [Google Scholar] [CrossRef] [Green Version]
  16. Garschke, C.; Parlevliet, P.P.; Weimer, C.; Fox, B.L. Cure kinetics and viscosity modelling of a high-performance epoxy resin film. Polym. Test. 2013, 32, 150–157. [Google Scholar] [CrossRef]
  17. Chen, D.-L.; Chiu, T.-C.; Chen, T.-C.; Chung, M.-H.; Yang, P.-F.; Lai, Y.-S. Using DMA to Simultaneously Acquire Young’s Relaxation Modulus and Time-dependent Poisson’s Ratio of a Viscoelastic Material. Procedia Eng. 2014, 79, 153–159. [Google Scholar] [CrossRef] [Green Version]
  18. Péron, M.; Sobotka, V.; Boyard, N.; Le Corre, S. Bulk modulus evolution of thermoset resins during crosslinking: Is a direct and accurate measurement possible? J. Compos. Mater. 2017, 51, 463–477. [Google Scholar] [CrossRef] [Green Version]
  19. Feng, J.; Guo, Z. Temperature-frequency-dependent mechanical properties model of epoxy resin and its composites. Compos. Part B Eng. 2016, 85, 161–169. [Google Scholar] [CrossRef]
  20. Shah, D.U.; Schubel, P.J. Evaluation of cure shrinkage measurement techniques for thermosetting resins. Polym. Test. 2010, 29, 629–639. [Google Scholar] [CrossRef] [Green Version]
  21. Minakuchi, S. In situ characterization of direction-dependent cure-induced shrinkage in thermoset composite laminates with fiber-optic sensors embedded in through-thickness and in-plane directions. J. Compos. Mater. 2015, 49, 1021–1034. [Google Scholar] [CrossRef]
  22. Seers, B.; Tomlinson, R.; Fairclough, P. Residual stress in fiber reinforced thermosetting composites: A review of measurement techniques. Polym. Compos. 2021, 42, 1631–1647. [Google Scholar] [CrossRef]
  23. Chevalier, J.; Brassart, L.; Lani, F.; Bailly, C.; Pardoen, T.; Morelle, X.P. Unveiling the nanoscale heterogeneity controlled deformation of thermosets. J. Mech. Phys. Solids 2018, 121, 432–446. [Google Scholar] [CrossRef]
  24. Yuan, Z.; Wang, Y.; Yang, G.; Tang, A.; Yang, Z.; Li, S.; Li, Y.; Song, D. Evolution of curing residual stresses in composite using multi-scale method. Compos. Part B Eng. 2018, 155, 49–61. [Google Scholar] [CrossRef]
  25. Yuan, Z.; Zhang, B.; Yang, G.; Yang, Z.; Tang, A.; Li, S.; Li, Y.; Zhao, P.; Wang, Y. Multi-Scale Modeling of Curing Residual Stresses in Composite with Random Fiber Distribution into Consideration. Appl. Compos. Mater. 2019, 26, 983–999. [Google Scholar] [CrossRef]
  26. Hui, X.; Xu, Y.; Zhang, W. An integrated modeling of the curing process and transverse tensile damage of unidirectional CFRP composites. Compos. Struct. 2021, 263, 113681. [Google Scholar] [CrossRef]
  27. Melro, A.R.; Camanho, P.P.; Andrade Pires, F.M.; Pinho, S.T. Micromechanical analysis of polymer composites reinforced by unidirectional fibres: Part I—Constitutive modelling. Int. J. Solids Struct. 2013, 50, 1897–1905. [Google Scholar] [CrossRef] [Green Version]
  28. Melro, A.R.; Camanho, P.P.; Andrade Pires, F.M.; Pinho, S.T. Micromechanical analysis of polymer composites reinforced by unidirectional fibres: Part II—Micromechanical analyses. Int. J. Solids Struct. 2013, 50, 1906–1915. [Google Scholar] [CrossRef] [Green Version]
  29. Catalanotti, G. On the generation of RVE-based models of composites reinforced with long fibres or spherical particles. Compos. Struct. 2016, 138, 84–95. [Google Scholar] [CrossRef] [Green Version]
  30. Barbero, E.J. Finite Element Analysis of Composite Materials Using AbaqusTM.; CRC Press: Boca Raton, FL, USA, 2013; ISBN 9781466516632. [Google Scholar]
  31. Yvonnet, J.; Auffray, N.; Monchiet, V. Computational second-order homogenization of materials with effective anisotropic strain-gradient behavior. Int. J. Solids Struct. 2020, 191–192, 434–448. [Google Scholar] [CrossRef]
  32. Coenen, E.W.C.; Kouznetsova, V.G.; Geers, M.G.D. Novel boundary conditions for strain localization analyses in microstructural volume elements. Int. J. Numer. Methods Eng. 2012, 90, 1–21. [Google Scholar] [CrossRef]
  33. Otero, F.; Oller, S.; Martinez, X. Multiscale Computational Homogenization: Review and Proposal of a New Enhanced-First-Order Method. Arch. Comput. Methods Eng. 2018, 25, 479–505. [Google Scholar] [CrossRef] [Green Version]
  34. Otero, F.; Oller, S.; Martinez, X.; Salomón, O. Numerical homogenization for composite materials analysis. Comparison with other micro mechanical formulations. Compos. Struct. 2015, 122, 405–416. [Google Scholar] [CrossRef] [Green Version]
  35. Toray Composite Materials America, Inc. T700 Technical Data Sheet. 2018. Available online: https://www.toraycma.com/resources/data-sheets/ (accessed on 2 June 2021).
  36. Soden, P. Lamina properties, lay-up configurations and loading conditions for a range of fibre-reinforced composite laminates. Compos. Sci. Technol. 1998, 58, 1011–1022. [Google Scholar] [CrossRef]
  37. Bai, X.; Bessa, M.A.; Melro, A.R.; Camanho, P.P.; Guo, L.; Liu, W.K. High-fidelity micro-scale modeling of the thermo-visco-plastic behavior of carbon fiber polymer matrix composites. Compos. Struct. 2015, 134, 132–141. [Google Scholar] [CrossRef] [Green Version]
  38. Murakami, S. Continuum Damage Mechanics; Solid Mechanics and Its Applications; Springer: Dordrecht, The Netherlands, 2012; Volume 185, ISBN 978-94-007-2665-9. [Google Scholar]
  39. de Souza Neto, E.A.; Peri, D.; Owen, D.R.J. Computational Methods for Plasticity; John Wiley & Sons, Ltd.: Chichester, UK, 2008; ISBN 9780470694626. [Google Scholar]
  40. Dai, J.; Xi, S.; Li, D. Numerical Analysis of Curing Residual Stress and Deformation in Thermosetting Composite Laminates with Comparison between Different Constitutive Models. Materials 2019, 12, 572. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  41. Melro, A.R.; Camanho, P.P.; Pinho, S.T. Influence of geometrical parameters on the elastic response of unidirectional composite materials. Compos. Struct. 2012, 94, 3223–3231. [Google Scholar] [CrossRef]
  42. Tan, W.; Martínez-Pañeda, E. Phase field predictions of microscopic fracture and R-curve behaviour of fibre-reinforced composites. Compos. Sci. Technol. 2021, 202, 108539. [Google Scholar] [CrossRef]
Figure 1. Example of a RVE.
Figure 1. Example of a RVE.
Polymers 14 02653 g001
Figure 2. Model solution procedure.
Figure 2. Model solution procedure.
Polymers 14 02653 g002
Figure 3. Randomly generated RVE geometries used for the micro-residual stress analysis. (a) RVE1, (b) RVE2 and (c) RVE3.
Figure 3. Randomly generated RVE geometries used for the micro-residual stress analysis. (a) RVE1, (b) RVE2 and (c) RVE3.
Polymers 14 02653 g003
Figure 4. Results from DMA analysis. (a) Storage and loss modulus (b) Curing state and shift function.
Figure 4. Results from DMA analysis. (a) Storage and loss modulus (b) Curing state and shift function.
Polymers 14 02653 g004
Figure 5. (a) Shift function for storage modulus and (b) Elastic properties evolution with curing level.
Figure 5. (a) Shift function for storage modulus and (b) Elastic properties evolution with curing level.
Polymers 14 02653 g005
Figure 6. Coupons after testing. (a) Transverse tensile tests, (b) Longitudinal shear tests and (c) Transverse shear tests.
Figure 6. Coupons after testing. (a) Transverse tensile tests, (b) Longitudinal shear tests and (c) Transverse shear tests.
Polymers 14 02653 g006
Figure 7. RVE micro-residual von Mises stress distribution for (a) RVE1, (b) RVE2 and (c) RVE3.
Figure 7. RVE micro-residual von Mises stress distribution for (a) RVE1, (b) RVE2 and (c) RVE3.
Polymers 14 02653 g007
Figure 8. RVE micro-residual hydrostatic stress distribution for (a) RVE1, (b) RVE2 and (c) RVE3.
Figure 8. RVE micro-residual hydrostatic stress distribution for (a) RVE1, (b) RVE2 and (c) RVE3.
Polymers 14 02653 g008
Figure 9. RVE average hydrostatic and von Mises stress evolution during the curing process.
Figure 9. RVE average hydrostatic and von Mises stress evolution during the curing process.
Polymers 14 02653 g009
Figure 10. Effective RVE transversal and longitudinal strain evolution during curing process.
Figure 10. Effective RVE transversal and longitudinal strain evolution during curing process.
Polymers 14 02653 g010
Figure 11. RVE volume change during curing process.
Figure 11. RVE volume change during curing process.
Polymers 14 02653 g011
Figure 12. Failure indexes for yielding and damage initiation.
Figure 12. Failure indexes for yielding and damage initiation.
Polymers 14 02653 g012
Figure 13. Transverse tensile stress-strain response comparison.
Figure 13. Transverse tensile stress-strain response comparison.
Polymers 14 02653 g013
Figure 14. Shear stress-strain response comparison. (a) Transverse direction and (b) Longitudinal direction.
Figure 14. Shear stress-strain response comparison. (a) Transverse direction and (b) Longitudinal direction.
Polymers 14 02653 g014
Figure 15. Damage patterns at failure for the different loading cases (SDV14 refers to the damage state variable).
Figure 15. Damage patterns at failure for the different loading cases (SDV14 refers to the damage state variable).
Polymers 14 02653 g015
Figure 16. Micro-residual stresses and failure index/cure level sensitivity to material and process parameters. (a) Thermal expansion sensitivity results, (b) Chemical shrinkage sensitivity results, (c) Cure temperature sensitivity results, (d) Matrix elastic modulus sensitivity results.
Figure 16. Micro-residual stresses and failure index/cure level sensitivity to material and process parameters. (a) Thermal expansion sensitivity results, (b) Chemical shrinkage sensitivity results, (c) Cure temperature sensitivity results, (d) Matrix elastic modulus sensitivity results.
Polymers 14 02653 g016
Figure 17. Tensile stress–strain response sensitivity to material and process parameters. (a) Thermal expansion sensitivity results, (b) Chemical shrinkage sensitivity results, (c) Cure temperature sensitivity results and (d) Elastic modulus sensitivity results.
Figure 17. Tensile stress–strain response sensitivity to material and process parameters. (a) Thermal expansion sensitivity results, (b) Chemical shrinkage sensitivity results, (c) Cure temperature sensitivity results and (d) Elastic modulus sensitivity results.
Polymers 14 02653 g017
Figure 18. Impact of the average residual stresses on the strength reduction.
Figure 18. Impact of the average residual stresses on the strength reduction.
Polymers 14 02653 g018
Table 1. Carbon fiber material properties.
Table 1. Carbon fiber material properties.
PropertiesValue
Fiber diameter (mm)7 × 10−3
Density (kg/m3)1800
Specific Heat (kJ/kg·K)0.752
Thermal Conductivity (W/m·K)9.38
C T E L (°C−1)−0.38 × 10−6
C T E T (°C−1)6.94 × 10−6
E L (GPa)230
E T (GPa)15
v L (n.d)0.2
G L (GPa)15
G T (GPa)7
Table 2. Resin system material properties.
Table 2. Resin system material properties.
PropertiesValue
Density (kg/m3)1310
Specific Heat (kJ/kg·K)0.679
Thermal Cond. (W/m·K)0.15 1
C T E (10−6 °C−1)61.0 1
Shrinkage (%)2.0 1
T g (°C)135 2
E m (MPa)2850
v m (n.d)0.33
v p (n.d)0.30 1
S y + (MPa)55
S y (MPa)81 3
S u + (MPa)65
S u (MPa)93 3
ε f (%)2.6
G f l m (N/mm)0.12 3
1 Estimated from literature for similar epoxies. 2 From the storage modulus onset. 3 Estimated from experimental tests with glass fiber prepregs. The compression values were extrapolated from transverse compression tests and the fracture energy corresponds to the interlaminar fracture toughness measured with DCB tests according to ASTM 5528.
Table 3. Loading cases to identify the effective composite properties.
Table 3. Loading cases to identify the effective composite properties.
TestDescriptionElastic Tensor ParametersYield Surface ParametersFailure Surface Parameters
Longitudinal uniaxial tension Polymers 14 02653 i001 E L ; v L --
Transverse uniaxial tension 1 Polymers 14 02653 i002 E T ;   v T Y + U T S + U T
Transverse shear Polymers 14 02653 i003 G T Y S T S S T
Longitudinal shear Polymers 14 02653 i004 G L Y S L S S L
1 Although v T is not required to identify the elastic tensor, it can be easily determined.
Table 4. Resin curing cycle.
Table 4. Resin curing cycle.
Time (min)023.2538.2550110130
Temperature (°C)25858512012020
Table 5. Set of parameters used in the sensitivity analysis.
Table 5. Set of parameters used in the sensitivity analysis.
ParameterValue
NominalCase 1Case 2Case 3Case 4
Shrinkage (mm/mm)−0.02−0.01−0.015−0.025−0.03
CTE (10−6 °C−1)60204080100
Elastic mod. (MPa)22802565285031353420
T cure (°C)12090105130-
Table 6. Fitting parameters for the curing kinetics model.
Table 6. Fitting parameters for the curing kinetics model.
A 1   ( 1 / s ) Δ E 1   ( kJ / mol ) A 2   ( 1 / s ) Δ E 2   ( kJ / mol ) n m
3.24 × 1015134,627001.000
Table 7. Ply material properties measured from experimental tests.
Table 7. Ply material properties measured from experimental tests.
Test StandardPropertyAverageStandard Deviation
ASTM D3039 E T (MPa)7313192
ASTM D3039 v L (n/d)0.01090.0044
ASTM D5379 G T (MPa)2812.873.9
ASTM D3518 G L (MPa)3215187
ASTM D3039 S + U T (MPa)42.914.15
ASTM D5379 S S T (MPa)22.342.17
ASTM D3518 S S L (MPa)21.411.577
Table 8. Effective mechanical properties comparison.
Table 8. Effective mechanical properties comparison.
Experimental
Results
Numerical PredictionsExperiments vs.
Cure Predictions (%)
CureNo-Cure
E 22 (MPa)7313 ± 1927151 ± 1147235 ± 552.3
v 21 (n/d)0.0109 ± 0.00440.0115 ± 0.00020.0115 ± 0.0002−5.1
v 23 (n/d)0.34 10.34 ± 0.010.34 ± 0.01−0.6
G 23 (MPa)2729 ± 1322674 ± 712733 ± 742.0
G 12 (MPa)3215 ± 1873353 ± 633460 ± 80−4.3
S + U T (MPa)42.9 ± 4.246.0 ± 2.051.9 ± 1.8−7.3
S S T (MPa)22.3 ± 2.225.0 ± 2.632.3 ± 0.7−12.2
S S L (MPa)21.4 ± 1.623.5 ± 1.226.8 ± 0.8−9.6
1 Estimated from shear modulus.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Gonçalves, P.T.; Arteiro, A.; Rocha, N.; Pina, L. Numerical Analysis of Micro-Residual Stresses in a Carbon/Epoxy Polymer Matrix Composite during Curing Process. Polymers 2022, 14, 2653. https://doi.org/10.3390/polym14132653

AMA Style

Gonçalves PT, Arteiro A, Rocha N, Pina L. Numerical Analysis of Micro-Residual Stresses in a Carbon/Epoxy Polymer Matrix Composite during Curing Process. Polymers. 2022; 14(13):2653. https://doi.org/10.3390/polym14132653

Chicago/Turabian Style

Gonçalves, Paulo Teixeira, Albertino Arteiro, Nuno Rocha, and Luis Pina. 2022. "Numerical Analysis of Micro-Residual Stresses in a Carbon/Epoxy Polymer Matrix Composite during Curing Process" Polymers 14, no. 13: 2653. https://doi.org/10.3390/polym14132653

APA Style

Gonçalves, P. T., Arteiro, A., Rocha, N., & Pina, L. (2022). Numerical Analysis of Micro-Residual Stresses in a Carbon/Epoxy Polymer Matrix Composite during Curing Process. Polymers, 14(13), 2653. https://doi.org/10.3390/polym14132653

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop