Abstract
The primary objective of this research is to develop a multiscale simulation framework to arrive at more realistic estimates of cure-induced residual stresses in the vicinity of the fiber-matrix interface in thick thermoset composites. The methodology involves simulations at the part level—employing homogenized rendering of the composite using micromechanics approach—within a finite element framework to obtain part-level temperature and degree-of-cure gradients and strains, and imposition of this information as boundary conditions at the mesoscale simulations, employing microstructural representative volume elements (RVE). A simple implementation of the multiscale framework, involving simulations at the part as well as the RVE levels, is demonstrated in the context of a thick, unidirectional continuous-glass-fiber-reinforced thermoset composite. The trends in the mesoscale residual stresses estimated by employing different RVE-level thermal and thermomechanical boundary conditions—displaying different degrees of coupling between the global and part-level simulations—are then examined. Significant differences are observed in the estimates of mesolevel cure-induced residual stress evolution obtained from simulations with a conventional symmetric RVE and those obtained by employing the multiscale approach involving detailed boundary conditions that realistically account for global thermal and mechanical strain histories.
1. Introduction and Background
During fabrication and cure of fiber-reinforced thermoset-matrix composites, chemical shrinkage of the thermoset matrix, the differential thermal shrinkage of the fiber and the matrix, and the constraints arising due to mold-composite interaction result in residual stresses in the resulting parts. In thick composites, additional large cross-thickness gradients in temperature and extent of cure (and the resulting resin property evolution) develop due to the exothermic nature of the cure reaction and the low thermal conductivity of the resin. Combined with microstructural defects and inhomogeneities, such large gradients in temperature and cure can result in residual stresses that can be a significant fraction of the matrix strength and/or the interfacial shear strength [1, 2] and may be comparable to or higher than the stresses encountered during service conditions [3–6]. The phenomenological aspects of residual stress development during thermoset composite cure have been enunciated in detail in several earlier reviews (e.g., [7]). At the length scales of the part, residual stresses result in shape distortion defects, such as warpage and spring-in. At the length scales of the reinforcement, residual stresses may become comparable to the strength of the matrix/interfacial adhesive strength and impact the strength of the composite as a whole by developing microcracks, delamination, and fiber-matrix debonding [4, 5]. Cure process models provide a valuable tool to understand the effect of composite microstructure, part geometry, and processing parameters on the development of residual stresses; such models also provide a relatively inexpensive virtual means of optimizing cure processing to minimize the residual stresses. Process models for composite cure are available that address residual stress development at the part level [6, 8–18] or at the “meso” levels involving the microstructure or fiber-matrix interface [19–23]. However, a systematic framework for bridging the information at the two scales has not been demonstrated; such bridging of information is critical for obtaining realistic estimates of residual stresses at the fiber matrix interface and for arriving at an understanding of the efficacy of reinforcement. This study is focused on demonstrating a multiphysics multiscale modeling approach—in the context of a thick unidirectional glass-fiber-reinforced thermoset composite—that effectively channels information from part-level simulations in arriving at estimates of residual stresses at the fiber-matrix interface. In the remaining portion of this section, the current state of the art in terms of process models for cure induced residual stresses is summarized (Sections 1.1 and 1.2) in order to highlight the need for a detailed multiphysics multiscale modeling approach; the objectives addressed in this study and the multiscale modeling framework are then presented (Section 1.3).
1.1. Background: Simulation of Residual Stresses at Part Level
Mathematical modeling of residual stresses at the length scales of the part or laminate provides insights into the interplay between the various thermochemical and mechanical phenomena driving the shape distortions. Typically, these models involve three aspects: first, the coupled solution of cure kinetics and heat transfer equations to obtain the temperature and degree of cure profiles in the part; this is especially important for thick composites [8–10]; second, the development of models for resin, laminate, and composite property evolution with cure and temperature; and third, the solution of the structural mechanics equations to estimate the induced strains and residual stresses in the part. Models for composite property evolution during cure can be based on direct experimental measurements on the composite itself during cure [11–14]. Alternatively, the properties of the matrix resin are measured as a function of the degree of cure and temperature. Assuming the fiber properties to be constant, the properties of the composite are then estimated either using self-consistent field micromechanics theories or 3D laminate theory [6, 12, 15–17]. Estimation of residual stresses and shape distortion effects in the composite part is obtained either by employing the composite laminate theory [5, 6, 15], finite element methods [5, 11, 14, 16, 17], or closed-form analytical solutions [18].
1.2. Background: Simulation of Residual Stresses at Reinforcement Level
At the smaller length scales corresponding to those of the reinforcement, the composite local microstructure becomes of great significance. The matrix and reinforcement phases are treated as separate entities in residual stress modeling at the mesoscales. This allows for detailed description of the matrix property evolution during cure and better estimation of residual stresses at the matrix-fiber interface. Consequently, this framework may also be employed to study localized damage initiation and propagation. Typically, mesoscale models for continuous fiber-reinforced composites are constructed by employing a representative volume element (RVE) which best represents the repeating microstructural unit of the composite. The RVE may be modified in order to study the effects of random spatial heterogeneities in orientation, aggregation, and distribution of the reinforcement or defects, such as a void. Alternatively, an actual composite micrograph may be digitized and employed as the geometry to carry out the simulations. The accuracy and relevance of previous attempts of modeling residual stresses at the mesoscale have been limited by the extent of the details of cure phenomenology and the composite microstructure captured by the models.
Zhang et al. [19] modeled the evolution of cure-induced residual stresses in cross-ply laminates by employing a thermoviscoelastic model of a 3D RVE accounting for both 0° and 90° orientations. This model however only accounted for residual stress build-up during cooling from the cure temperature and therefore did not account for cure-induced shrinkage. The model considered temperature-dependent modulus and cure-independent glass transition temperature and captured the build-up and relaxation of residual stresses. Andersson et al. [20] modeled the stress development during cool-down from cure temperature by employing the digitized micrograph of an actual composite as the simulation geometry. While the model employed a real microstructure, thermoelastic properties of the resin were assumed to be temperature independent. Also, cure-induced shrinkage and the stress development during cure were not accounted for. Rosso et al. [21] estimated mesolevel cure-induced residual stresses by employing a plane-strain thermoelastic model of a unidirectional carbon fiber composite with hexagonal packing; however, they incorporated cure-shrinkage in an ad hoc manner by modifying the coefficient of thermal expansion of the resin only for a fixed temperature range (thereby not truly accounting for cure-dependent gelation effects). More recently, Zhao and coworkers have reported attempts to capture the cure phenomenology at the RVE level in greater detail with separate accounting for cure-induced shrinkage and thermal expansion effects in unidirectional composites by employing thermoelastic models [22] and thermoviscoelastic models [23]. These RVE simulations involved modeling of isothermal cure and cool-down stages of the thermal cycle. Chemical shrinkage was assumed to follow a linear dependence on conversion. Zhao et al. [22, 23] ignored the interdependence of cure and temperature, which might be a valid assumption for isothermal curing systems with low heats of reaction. However, such an assumption would not be valid for an RVE that is located in the core of a thick composite with a resin that has a very high heat of reaction. The RVE/unit cell studies summarized above involve simulations on idealized microstructures in isolation. Thermomechanical boundary conditions do not take into account part level gradients and the location of the microstructural feature within the part.
1.3. Motivation and Objective
The above discussion makes it clear that, on the one hand, part-level models do not provide information about stress distribution within the different phases of the composite and therefore may not be reliably interpreted to explain the efficacy of reinforcement. On the other hand, while the existing fiber-level models can provide information relevant to matrix microcrack development and fiber-matrix debonding, often they are based on simplistic or inaccurate assumptions regarding cure-dependent property evolution; more significantly, they do not account for the part-level variations in processing conditions and property evolution, which can be significant for thick composites. Therefore, the primary objective of this research is to develop and demonstrate a comprehensive multiphysics multiscale modeling framework that systematically combines the information from global processing conditions and composite microstructural details, along with detailed models for cure-dependent material properties, to arrive at more realistic estimates of residual stresses in the vicinity of the fiber-matrix interface.
The methodology involves simulations at the part level—employing homogenized rendering of the composite using micromechanics approach—to obtain part level gradients and strains and imposition of this information as boundary conditions at the mesoscale simulations employing microstructural RVEs. The implementation of the methodology is facilitated by multiphysics modeling strategies which have enabled detailed and comprehensive accounting of cure phenomenology (see, e.g., [24]). The multiscale framework demonstrated here may be broadly divided into three steps: the first step involves thermomechanical characterization of reinforcing fibers, measurement of cure-kinetics and cure-dependent properties of the matrix, and development of correlations for cure-dependent property-evolution of the composite as a whole. The second step involves feeding the property-evolution models into part-level cure simulations to arrive at estimates of cure and temperature gradients as well as spatiotemporal variation of stresses and strains within the part. The final step involves channeling the part-level information of cure, temperature, and strain transients as boundary conditions for mesolevel cure simulations to estimate residual stresses at the fiber-matrix interface.
In the subsequent, a simple implementation of the multiscale modeling framework is demonstrated in the context of a thick unidirectional continuous-glass-fiber-reinforced thermoset polyester matrix composite; this polyester system has been characterized by Bogetti and Gillespie [9, 15]. The implementation of the first step of the multiscale modeling framework is described in Section 2; in this section the correlations derived by Bogetti and Gillespie [15] for calculation of the average properties of the unidirectional composite using self-consistent field micromechanics theory are summarized. The implementation of the second step of the framework involving the part-level simulation of cure is described in Section 3. Section 4 describes the implementation of the final step of the multiscale framework in which information regarding estimated global variations in temperatures, conversion, and the associated part-level deformation is incorporated in the mesolevel simulations carried out using a 2D RVE. Different degrees of coupling between the part-level and RVE-level simulations were employed to clearly establish the expediency of employing the multiscale modeling framework. The estimates of mesolevel stresses obtained by employing RVEs with thermal or thermomechanical boundary conditions that are more truly representative of the prevalent magnitudes and gradients of global processing variables are then compared with those obtained using a conventional symmetric RVE.
2. Definition of Fiber, Matrix and Composite Properties
In this section, the material properties of the fiber and matrix components of the composite, as well as the composite properties, are described. The composite system, composed of a polyester matrix unidirectionally reinforced with continuous glass fibers, was studied by Bogetti and Gillespie [9, 15]. The properties of the resin, reinforcing fiber, and the composite are described subsequently in terms of three principal directions. “1” is the direction of alignment of the fiber reinforcement, “2” is the in-plane transverse direction, and “3” is the out of plane direction with respect to the composite. The composite itself is assumed to be aligned with the global coordinate system such that 1, 2, and 3 directions with respect to the composite correspond to the global -, the global -, and the global -axes, respectively, as shown schematically in Figure 1.
The glass-fiber reinforcement was modeled as an isotropic elastic material [9, 15] with a constant (cure and temperature independent) elastic modulus, , shear modulus, , and Poisson’s ratio, . The coefficient of thermal expansion of the fibers was also assumed to be constant (temperature independent). The expressions for these variables, as derived in [15], are listed in Table 1, and the fiber thermomechanical properties are listed in Table 2. The matrix was modeled as an isotropic elastic material with cure-dependent elastic modulus using a rule-of-mixtures type relationship [15] (see Table 1). This relationship interpolates the cure-dependent modulus of the resin between the limiting moduli of the fully uncured resin, , and the fully cured resin, , using a function of the instantaneous conversion , which involves the degree of cure corresponding to the gel point of the resin, , and the limiting conversion, , corresponding to the onset of diffusion controlled phenomena, beyond which the modulus does not develop further with conversion. The function is defined below [15]:
In this study, following [15], these model parameters were assigned the following values: = 0 (modulus evolution right from the beginning) and = 1.0 (no diffusion controlled regime). Also, the temperature dependence of the modulus of the fully cured resin, as well as any vitrification effects during the evolution of the modulus, were ignored in this analysis. The expressions of the constant isotropic Poisson’s ratio and the cure-dependent shear modulus of the resin are also shown in Table 1. The thermal conductivity of the resin was assumed to be constant (temperature and cure independent) and isotropic. The evolution of cure-induced shrinkage was assumed to follow a linear dependence on the degree [15] of cure as shown below:
In the above equation, is the overall composite volume change at the completion of cure per unit initial volume. The limiting resin moduli and the resin Poisson’s ratio, , the coefficient of thermal expansion of the resin, , and are listed in Table 3.
Based on the fiber thermomechanical properties and the evolution of the matrix resin properties during cure (refer to Table 1), the evolution of the transversely isotropic thermomechanical properties of the composite with the progress of cure was modeled using the self-consistent field micromechanics model for unidirectional continuous fiber reinforced composites as described by Bogetti and Gillespie [15]. The expressions for the composite variables as derived by Bogetti and Gillespie [15] are listed in Table 4.
3. Multiscale Framework: Part Level Simulation
The part-level finite element simulations, with the homogeneous rendering of the composite properties described in the earlier section, were carried out by employing a 3D planar geometry of a cuboidal shape, with 12.7 cm length, 12.7 cm width, and 1.27 cm thickness. As schematically shown in Figure 2, this geometry represents one half of the thickness of the symmetric composite with the bottom face (z = 0) representing the core and the symmetry plane (the symmetric representation ignores any differential mechanical constraints arising due to the contact between the top and bottom composite surfaces and the autoclave/supporting frame and assumes that both the composite surfaces are exposed to the autoclave thermal cycle). The length and width of the composite geometry are aligned with the global - and global -axes, respectively, which are also the composite longitudinal and in-plane directions, respectively (cf. Figure 2).
(a)
(b)
3.1. Kinetics and Heat Transfer
The cure kinetics model for the composite, described by Bogetti and Gillespie [9, 15] based on a calorimetric study of the cure, is summarized below:
In the above equation, the subscript denotes that the kinetic parameters are estimated for the composite as a whole. The overall reaction order for the cure of the glass/polyester composite is 2. The preexponential factor, , the activation energy, , and the exponents, and , are listed in Table 5. R is the universal gas constant. Transient heat conduction in the curing system is described as shown below:
In the above equation,ρ, , and are, respectively, the density, specific heat, and thermal conductivity of the composite, and these are also listed in Table 5. The rate of heat evolution during cure is in turn governed by the overall enthalpy of the cure reaction and the instantaneous rate of reaction:
The kinetic model (3) was solved with an initial concentration of 1.0 mol/L. The heat transfer problem (4) and (5) was solved with a surface convective boundary condition at the top (z = 1.27 cm) surface of the composite which is exposed to the autoclave temperature cycle (cf. Figure 2(a)):
The time-dependent autoclave thermal cycle is presented in Figure 3. On all other edges, including the bottom surface, which is representative of the composite core, a thermal insulation boundary condition was employed:
Following Bogetti and Gillespie [9, 15], was set to 0, and was set to 1.
3.2. Laminate Level Structural Mechanics Problem Definition
The total strain developed within the composite during cure can be divided into two parts: mechanical strains and expansional strains that are caused by thermochemical effects during cure:
The expansional strains in the composite are the resultant of the expansional strains in the fiber and those in the matrix. The expansional strains in the fiber are entirely composed of thermal strains induced by thermal expansion/shrinkage associated with temperature change:
In the above equation, is the coefficient of linear thermal expansion (CLTE), in a given principal direction i (with respect to the composite reinforcement), and is the reference temperature. On the other hand, in the resin matrix, the expansional strains are composed of thermal strains, as well as shrinkage strains brought about by cure induced shrinkage during chemical crosslinking [15]:
Based on the expressions for thermochemical strains in the fiber and the matrix resin, the transversely isotropic expansional strains in the composite as derived using the self-consistent field model approach [15] are shown below:
Since the composite is assumed to remain linearly elastic throughout its cure, the Cauchy mechanical/residual stress tensor is proportional to the mechanical strain tensor [25]:
In the above equation, C is a function of the effective composite modulus and Poisson’s ratio of the material. The momentum balance equation takes the form shown in (13) in the absence of body forces [25]:
The momentum balance was implemented in the geometry shown in Figure 2(b). The composite was assumed to be infinite in the longitudinal x- and the in-plane transverse y- directions. At the four faces that are normal to these directions, a straight-face or equal-normal-displacement boundary condition was imposed as an ad hoc measure to maintain periodicity (cf. Figure 2(b)):
The thermochemomechanical material constants required to evaluate the composite engineering properties and the expansional strains are listed in Table 2 for the fiber and Table 3 for the matrix resin. For the calculations, the reference temperature, , was taken to be 293 K. The coupled cure-kinetics, heat-transfer, and structural mechanics problem defined above was implemented in COMSOL Multuphysics Version 3.5a [26]. As shown in Figure 2, the evolution of the key process variables, such as temperature, conversion, and strains, was tracked at two locations within the composite, one on the top surface [point A] and the other within the core of the composite [point B].
3.3. Part-Level Estimations
Figure 4 shows the estimated temperature profiles at the core (point B) and the surface (point A) of the composite during the course of cure. The corresponding conversions of the cross-linking reaction at the surface and the core are plotted in Figure 5(a). It should be noted that the long cure cycle, lasting for nearly six hours, is typical of cure of thick composites. As seen in Figure 4, the temperature transient at the part surface precisely follows the autoclave temperature cycle. For approximately the first 9000s, the temperature rise in the core clearly trails that at the surface owing to the low thermal conductivity of the composite. For the same initial duration, the degrees of cure are lower in the core compared to those at the surface, as can be seen from Figure 5(a). As the cure progresses, due to the low thermal conductivity of the composite, the exothermic heat in the core cannot be dissipated fast enough, and the temperature at the core shows a rapid increase, crossing the temperature curve at the surface at around 9000s. This in turn leads to faster reaction at the core and, therefore, more heat due to the reaction exotherm. The maximum temperature reached at the core, at about 9840s, is almost 40°C higher than that imposed on the surface. The estimates in Figure 4 are similar to those obtained by Bogetti and Gillespie (cf. [9]), verifying the COMSOL implementation of the thermochemical coupling model. At the timescales corresponding to the temperature maximum at the core (cf. Figure 4), the conversion at the core also shows a rapid increase, significantly exceeding that at the surface and rapidly reaching complete conversion at the time corresponding to the temperature maximum (cf. Figure 5(a)). After the temperature maximum corresponding to 9840s, as the conversions reach limiting values at the core, the reaction rates in the core decline. This results in the exhaustion of the exothermic heat and a gradual drop of the temperatures. The core once again trails the temperatures at the surface subsequently. As seen in Figure 5(a), the conversions at the surface require substantially longer times to reach the limiting values compared to the core.
(a)
(b)
The evolution of composite moduli at the corresponding locations within the composite is plotted in Figure 5(b). The composite moduli clearly follow the cure transients at the surface and the core. At times corresponding to the temperature maximum associated with the cure exotherm, the moduli in the core rapidly exceed those at the surface, reaching the limiting values. It is also clear that the longitudinal modulus, , dominated by the glass fiber reinforcements, is substantially higher than the transverse in-plane and out-of-plane moduli ( and , resp.) and is also close to its limiting value even at the beginning of the cure cycle. The cure of matrix clearly adds only a small percentage of the total magnitude of , whereas it amounts to a significant portion of the limiting values of and .
The evolution of the principal residual stresses in the composite at the surface (point A, cf. Figure 2) and the core (point B in Figure 2) is plotted in Figures 6(a) and 6(b), respectively. A quick inspection of Figures 6(a) and 6(b) reveals that, in the initial as well as the final stages of the cure, the longitudinal and in-plane transverse principal stresses typically assume similar magnitudes that are significantly higher than the out-of-plane stresses. For the first 1800s, during the first temperature ramp, the in-plane transverse and the out-of-plane stresses are negligible compared to the longitudinal stresses. During this first positive temperature ramp, the chemical reaction, and consequently, the crosslink shrinkage has not set in. Therefore, the composite has a tendency to expand. While this expansion can occur freely in the in-plane transverse and out-of-plane directions, the composite is inherently constrained in the longitudinal direction due to the alignment of the continuous fiber reinforcements in the x-direction, as also evident from the high magnitudes of even at the beginning of the cure cycle. Due to the fact that the temperature at the core significantly trails that at the surface during this initial period, the tendency to expand is also relatively higher at the surface compared to the core. Consequently, driven by the inherent fiber constraints, the surface develops a compressive longitudinal stress, and this causes the development of tensile longitudinal stresses in the core. These thermal gradient driven initial longitudinal stresses at the surface and the core die out, as the temperature is equilibrated between the surface and the core during the first temperature hold that lasts up to 5000s. The second temperature ramp imposed between about 5000 and 9000s is not as steep as the first temperature ramp. Moreover, the chemical reaction has still not reached high enough levels of conversion for crosslink shrinkage to set in a significant manner. Therefore, no significant additional stresses may be observed in any direction in the composite between 3000s and 9000s.
(a)
(b)
As the cure exotherm sets in at the core around 9000s, resulting in rapid increase in the degree of cure, the matrix in the core shows a tendency to shrink. This shrinkage occurs relatively unconstrained in the out-of-plane direction in the absence of any mechanical constraints but is inherently locally constrained by the presence of fibers in the longitudinal and in-plane-transverse directions. As a result, the composite rapidly develops significant magnitude of tensile stresses in these directions at the core, further aided by the increase in the composite moduli. The magnitudes of the in-plane-transverse tensile stresses are slightly higher than the longitudinal tensile stresses. These higher magnitudes of the in-plane transverse stresses arise from the satisfaction of the uniform-edge-displacement periodicity boundary condition (14), which clearly incurs a higher penalty in the in-plane transverse direction compared to the longitudinal direction (which is inherently constrained due to the high magnitudes of ). Of course, it should be noted that, in the absence of any external constraints, the tensile stress peak in the core is primarily driven by the large temperature and conversion gradients that set in along the thickness of the composite during cure. As the temperature and degree of cure gradients decrease, there is a drop in the tensile longitudinal and in-plane-transverse tensile stresses from their peak values, but the tensile stresses do not fully relax, because as the cure exotherm dies out, the temperatures decrease, resulting in a further tendency to shrink and the development of residual tensile stresses. These residual tensile stresses finally disappear as the composite relaxes during the third temperature ramp up to 14400s.
As mentioned earlier, the stresses in the current simulation are primarily driven by the gradients in the temperature and degree of cure along the thickness of the composite, and not by any externally imposed constraints. Consequently, the longitudinal and in-plane-transverse stresses developed at the surface are typically of the opposite nature as compared to those at the core. At times corresponding to the tensile stress peak in the core, the surface develops longitudinal and in-plane-transverse compressive stresses. However, the longitudinal and in-plane-transverse stresses at the surface do not develop the intermediate compressive stress maximum analogous to the tensile stress maximum observed in the core. This is because the evolution of conversion and the consequent cross-link shrinkage is more gradual at the surface (refer to Figure 5). During the final cool-down, the surface undergoes more rapid temperature loss compared to the core. This causes a greater tendency for the material at the surface to shrink compared to that at the core. However, as the moduli are fully developed by this time, the resulting stiffness prevents free shrinkage, resulting in the development of residual tensile stresses on the surface, with the core developing compressive stresses. Once again, as the temperatures get finally equilibrated at the end of the cure cycle, the stresses are also minimal at the end of cure cycle.
The simulated evolution of principal strains in the longitudinal (x-) direction, the in-plane transverse (y-) direction, and the out-of-plane (z-) direction in the composite are plotted in Figures 7(a), 7(b), and 7(c), respectively. As seen in Figure 7(b), the in-plane transverse strains, , evolve in exactly the same manner in the surface (point [A] in Figure 2) and the core (point [B] in Figure 2) of the composite, as governed by the uniform-edge-displacement boundary condition (refer to (14)) imposed on all the side edges. Up to about 9000s, the strains are tensile in nature, indicative of free expansion in the in-plane transverse direction due to the ramping of temperature. The strains rapidly become compressive at times corresponding to the cure exotherm in the core and remain compressive thereafter for the remaining duration of cure. The magnitudes of compressive strains increase further during the final cool-down of the fully cured composite. By contrast the out-of-plane strains, (cf. Figure 7(c)), follow different transients in the surface (point [A] in Figure 2) and the core (point [B] in Figure 2) of the composite. As in the case of in-plane strains, the out-of-plane strains are also tensile for the first 9000s. In the core, the strains rapidly turn compressive at times corresponding to the cure exotherm. The compressive strains in the core are further enhanced up to about 12700s, due to the cool-down of the exotherm, and get partially relieved due to the further heating of the composite. The final cool-down further enhances the compressive strains, and the final state of strain in the composite core is compressive out-of-plane strain. In the surface, the out-of-plane strains reach the same final magnitude; however, the build-up of strains is more gradual, and follows the development of cure transient up to about 15000s. Once again, the final cool-down results in substantial addition to the compressive out-of-plane strains in the composite surface. Compared to and , the principal strains in the longitudinal direction are significantly lower (cf. Figure 7(a)). This may be attributed to the inherent longitudinal constraint imposed by the continuous fibers aligned in the x-direction. Also, like , evolve in exactly the same manner at the surface and the core of the composite due to the uniform-edge-displacement boundary conditions on the side edges. Whether this remains true in a true periodic boundary scenario (nonstraight edge periodicity using multipoint constraints on the edges) remains to be verified; this is the subject of ongoing investigations and will be reported subsequently.
(a)
(b)
(c)
4. Multiscale Framework: RVE Level Simulation
The 2D unit cell representation of the unidirectional composite microstructure is shown in Figure 8(a). The structural mechanics problem within the RVE was solved using the plane strain approximation, with the longitudinal (x-) direction considered perpendicular to the unit cell plane (this choice is further appropriate because, in the part-level analyses, the in-plane-transverse and out-of-plane residual strains, and , were found to be significantly higher than the longitudinal strains, ). In the RVE, the composite is represented using a 100μm by 100μm square. The fiber is represented as a circle of area the same as the overall volume fraction of the fibers in the composite. The remaining portion of the domain represents the matrix. Perfect adhesion was assumed between the matrix and the fiber. The evolution of cure-induced residual stresses at the reinforcement scale was tracked at three points within the matrix portion of the RVE, as indicated in Figure 8(a).
(a)
(b)
(c)
(d)
4.1. Kinetics and Heat Transfer
In the RVE level simulations, the kinetics was modeled in the same way as it was carried out at the part level (cf. (3)). The heat transfer equation (cf. (4)) however did not include the exothermic heat generation term. It was found that, given the small volume of the matrix resin accounted for in the RVE, relative to the length of the edge on which the thermal boundary condition was imposed, any exothermic heat generated was instantaneously equilibrated. Instead, the thermal boundary conditions for the RVE (discussed in Section 4.2) contained information regarding the temperature transients (and the cure-exotherm), obtained from the part-level simulations, at any given location along the thickness of the composite.
4.2. Impact of Thermomechanical Boundary Conditions
The impact of part-level variations in cure-induced deformation and temperature and cure transients on the development of residual stresses at the fiber-matrix interface was investigated by employing three types of thermomechanical boundary conditions at the RVE-level. In the first scenario, for which the thermomechanical boundary conditions are shown in Figure 8(b), mechanical symmetry (equal-normal-displacement at all points on an edge) was employed on all the external edges of the RVE, while subjecting it to the autoclave thermal cycle, a scenario that is equivalent to the boundary conditions employed by Zhao et al. [22, 23]. In the second scenario, for which the thermomechanical boundary conditions are shown in Figure 8(c), the evolution of residual stresses was studied in a mechanically symmetric RVE (equal-normal-displacement at all points on an edge) subjected to the temperature transient estimated for the core in the part-level composite cure simulations (cf. Figure 4). In the third scenario, for which the thermomechanical boundary conditions are shown in Figure 8(d), the RVE boundary conditions contain information about the temperature transient as well as the mechanical strain history at a chosen location in the composite part. The third scenario presents the most comprehensive accounting of the spatial variation of part level cure- and temperature-transients as well as deformation history. The trends in residual stress evolution within an RVE located at the surface and the core of the composite are then compared and contrasted for these three boundary condition scenarios.
Case 1: Simulations with Symmetric RVE and the Autoclave Thermal Cycle. The simplest thermomechanical boundary condition investigated was the global autoclave thermal cycle (corresponding to the temperature profile of the composite surface, point [A] in Figure 2) imposed on a symmetric RVE. This is the typical manner in which RVE simulations are carried out (see, e.g., [22, 23]), assuming the applicability of the imposed thermal cycle in the RVE irrespective of its location within the composite and ignoring any global deformations (hence the symmetric RVE). In this case, the thermomechanical boundary conditions are as shown schematically in Figure 8(b). On one of the outer edges of the RVE, the autoclave temperature cycle was imposed. The remaining outer edges were taken to be insulated. It should be noted here that this is only one of the several ways in which the RVE level heat transfer problem may be solved. The temperature boundary condition may be imposed on any one of the four outer edges of the RVE. The temperature boundary condition may also be simultaneously imposed on all the four outer edges of the RVE. Due to the very low thermal mass associated with the 100μm by 100μm RVE, the temperatures get equilibrated within the RVE almost instantaneously. In order to maintain a mechanically symmetric RVE while satisfying periodicity, an equal-normal-displacement boundary condition was imposed on all the outer edges of the RVE as shown in Figure 8(b).
The contour plot of the simulated maximum principal stress within the RVE at the end of the cure cycle (at 24000s) for this scenario is shown in Figure 9(a). It is clear that, at the end of cure, a significant portion of the matrix phase of the RVE develops extremely high magnitudes of tensile stresses. The tensile stresses estimated within the matrix are significantly higher than those within the fiber and tend to be highest at the fiber-matrix interface (the stresses at fiber-matrix interface are not very well resolved with the ten contour levels plotted in Figure 9(a)). This is due to dual constraints on the matrix offered by the significant mismatch of the CLTE between the fiber and the matrix and the controlled displacement boundary conditions imposed on the outer edges of the RVE. The maximum tensile stresses estimated within the matrix are of the same order as the tensile strength of the fully cured matrix. The simulated evolution of the maximum principal stresses within the matrix of the RVE, at three locations (indicated in Figure 8(a)), are plotted in Figure 9(b). The symmetric mechanical boundary conditions result in identical evolution of the predominantly compressive maximum principal stresses at the (50, 100) and (100, 50) locations within the matrix. The maximum principal stresses are tensile at the matrix-rich location (100, 100), driven by constrained shrinkage of the matrix phase. At all the three locations within the matrix, the increase in the magnitude of the maximum principal stresses between 10000s and 15000s is consistent with the increase in degree of cure up to its completion during the same time period. The second ramp in the stresses happens during the final cool down of the composite. As seen from the contour plot in Figure 9(a), this simulation predicts almost complete rupture of the matrix at all locations within the composite (since spatial variation of cure, temperature, and strain are not accounted for) due to the extremely high tensile stresses estimated with the symmetric RVE; clearly, this is not realistic.
(a)
(b)
Case 2: Simulations with Symmetric RVE and Part-Level Core Thermal Transient. The impact of part-level variation in temperature transients on the RVE level stresses is explored in Figure 10. Figure 8(c) depicts the modified thermal boundary condition for one of the edges of the RVE, in which the temperature transients at the core are imposed as the thermal boundary condition. Thermal insulation is applied on the remaining outer edges. In this scenario too, in order to maintain a mechanically symmetric RVE while satisfying periodicity, an equal-normal-displacement boundary condition was imposed on all the outer edges of the RVE. The simulated evolution of the maximum principal stresses in this scenario at the various locations within the matrix (cf. Figure 8(a)) is shown in Figures 10(a) and 10(b).
(a)
(b)
In Figure 10(a), the estimated maximum principal stress evolution at the (100, 100) location (matrix rich) within the RVE for this scenario is compared with the stress transient obtained from the simulation with the symmetric RVE subjected to the surface thermal transient (shown in Figure 9(b)). From an inspection of Figure 10(a), the impact of the cure exotherm within the core at around 9800s is obvious; the rapid completion of the cure reaction results in the generation of tensile stresses much more rapidly within the core RVE compared to the surface RVE. At times corresponding to the cure exotherm, the rapid evolution of constrained shrinkage in the matrix-rich corners of the RVE drives a steep increase in compressive residual stresses at the locations (100, 50) and (50, 100), as seen in Figure 10(b). The choice of symmetric mechanical boundary conditions results in identical evolution of stress transients at (100, 50) and (50, 100). It is worth noting that even though the evolution of the maximum principal stresses is different in the RVE simulations with autoclave and core thermal transients, the final magnitudes of the stresses at the end of cure cycle are identical at the respective locations within the matrix for the two scenarios. In fact, the contour plot of maximum principal stresses at the end of cure cycle (24000s) is identical for the two RVE simulation scenarios. This may be attributed to the fact that, in the part-level simulations (which provide the information for the RVE-level temperature boundary conditions), the final temperature (refer to Figure 4) and the degree of cure (refer to Figure 5) within the surface as well as the core are the same at the end of the thermal cycle.
Case 3: RVE Simulations with Part-Level Thermal and Strain Histories. In a complete implementation of the multiscale modeling involving the translation of global variations in thermal history and mechanical strains, the RVE level thermomechanical boundary conditions were implemented as shown in Figure 8(d). The surface and core thermal histories were implemented in the same manner as discussed in the descriptions of Cases 1 and 2, respectively. Additionally, the information about the in-plane transverse () and out-of-plane () principal strains in the composite at the surface and the core (cf. Figures 7(b) and 7(c)) were also imposed in the form of displacement boundary conditions on the RVE, as described in Figure 8(d). In Figure 11(a), the estimated maximum principal stress evolution at the (100, 100) location, simulated with the detailed thermomechanical boundary conditions for the surface and core, is compared with the stress transient obtained from the simulation with the symmetric RVE subjected to the autoclave thermal cycle (shown in Figure 9(b)). Figures 11(b) and 11(c) are analogous to Figure 11(a) and compare the respective transients at the (50, 100) and (100, 50) locations within the RVE.
(a)
(b)
(c)
(d)
From Figure 11(a), it is clear that, at (100, 100) (matrix-rich location), the simulated maximum principal stresses upon the completion (and for the most part) of the cure cycle remain tensile in nature, but the magnitudes of these tensile stresses are almost 50% lower as compared to those calculated in a symmetric RVE subjected to the autoclave temperature cycle. Again, even though the part-level out-of-plane residual strain () transients are very different at the surface and the core of the part, the RVE level stresses predicted by accounting for these strains, while evolving differently, end up at the same magnitude at the end of the cure cycle. This may be explained from an inspection of Figures 7(b) and 7(c), in which the estimates of the part-level in-plane-transverse and out-of-plane strains at the surface and the core match at the end of the cure cycle (since the evolution of residual stresses and strains at the part level were driven by the cure and temperature gradients along the thickness of the curing composite, and not by any externally imposed constraints). In the presence of external constraints, the final values of residual stresses and strains in the surface and the core of the part would have been different, and this would also have reflected in the stress evolution at the RVE level.
In contrast to the identical evolution of stresses at the matrix-poor locations, (100, 50) and (50, 100) in Figures 9 and 10 (obtained by accounting for only thermal history at the RVE level), the stresses at these two locations evolve differently when part-level mechanical strain history is also taken into account, as evident from inspection of Figures 11(b) and 11(c), respectively. The magnitudes of compressive stresses developed at these locations at the end of cure cycle are almost 100–150% higher compared to those predicted in a symmetric RVE subjected to the autoclave temperature cycle (due to the additional compressive strains imparted by accounting for part-level residual strains).
These aspects discussed in the foregoing are also clear from an inspection of Figure 11(d), which shows the contour plot of the simulated maximum principal stress within the surface RVE, simulated with the thermomechanical boundary conditions equivalent to the part-level temperature and residual strain transients at the part surface, at the end of the cure cycle (at 24000s). A quick comparison of this contour plot with that shown in Figure 9(a) reveals that accounting for global strain history results in significantly lower estimates for the tensile stresses within the matrix-rich regions of the RVE and a more widespread compressive state of stress in the matrix phase compared to the scenario in which only thermal transients are taken into account. Once again, it is worth noting that the contour plot obtained with the core RVE at the end of cure cycle is identical to that obtained for surface RVE.
It is now clear that not accounting for global strain history in the RVE level simulations results in a gross overestimation of the tensile stresses in the matrix-rich areas. At other locations it also results in the prediction of significant tensile stresses in the matrix, as compared to compressive stresses when the strain history is accounted for. Thus, simulations with a symmetric RVE with thermal boundary conditions, without accounting for the global strain history, may result in significant overestimation of the matrix damage during cure.
5. Summary and Future Work
The multiscale modeling framework for estimation of mesolevel cure-induced residual stresses has been demonstrated in the context of a thick, unidirectional, continuous-glass-fiber-reinforced thermoset polyester composite. The part-level simulation involves a homogenized rendering of the composite, with the average properties of the composite calculated using self-consistent field micromechanics theory. The 2D RVE simulations are conducted by employing three different types of thermomechanical boundary conditions that involve different degrees of coupling with the part-level simulations: a mechanically symmetric RVE, in which cure is simulated using a thermal cycle that corresponds to the globally imposed autoclave-thermal cycle, a conventionally employed choice of boundary conditions that is similar to Zhao et al. [22, 23]; a mechanically symmetric RVE, in which the thermal boundary conditions are more representative of its location within the composite; and an RVE, in which both the local temperature evolution and the global mechanical strain evolution are accounted for in the thermomechanical boundary conditions. Based on comparison of these scenarios, the following may be concluded.(1)Simulations with a symmetric RVE result in a gross overestimation of the tensile stresses in the matrix rich areas compared to the simulations involving the imposition of global strain history at the RVE-level. At other locations within the matrix phase of the RVE, simulations with symmetric RVEs result in the prediction of significant tensile stresses, while simulations accounting for global strain history result in compressive stresses. Thus, in this scenario, simulations with a symmetric RVE with only thermal boundary conditions, without accounting for the global strain history, may result in significant overestimation of the matrix damage during the cure.(2)From the current implementation, it may be concluded that, even in the case of thick composites with a significant exotherm, the final state of stress in the surface and core RVEs at the end of cure cycle is not significantly different, especially if the global residual strains are primarily driven by the cure and temperature gradients. On the other hand, the presence of externally imposed global mechanical constraints and modified periodic boundary conditions may result in significant differences in the final state of stress in the surface and core RVEs.
Having demonstrated the potential value of the multiscale modeling framework in arriving at estimates of residual stresses at the fiber-matrix interface, the following research problems are currently being investigated or proposed to further bring forth the applicability of such multiscale framework to scenarios of high practical relevance. The results from these studies will be presented in future.(1)As mentioned in the foregoing, the presence of externally imposed global mechanical constraints (such as mold part interaction) and modified periodic boundary conditions at the part level can result in significant differences in the final state of stress in the surface and core RVEs (the latter was not observed in the current study since mold constraints were not modeled). It is therefore important to implement the multiscale modeling approach to scenarios involving different degrees of part-level constraints during global cure processing, to understand how the part-level cure, temperature, and deformation transients would govern mesolevel residual stress development. Modification of periodic boundary conditions at the edge to allow for nonstraight edges would further improve the accuracy of residual stress estimation both at the part-level and mesolevels. These periodic boundary conditions need to be systematically arrived at.(2)In the current implementation, the matrix resin was simulated as a cure-dependent elastic material. In such rendition, the instantaneous properties of the matrix are governed only by the instantaneous degree of cure, and are not affected by the processing history. On the other hand, it was observed by Patham [24] that, by employing a more comprehensive cure- and temperature-dependent viscoelastic model for the thermoset resin, subtle details of the interplay between the processing history and spring-back behavior may be brought forth, which are not available from elastic models. Therefore it would be interesting to investigate the effect of viscoelasticity of the matrix resin in governing the residual stress development both at the part-level and at the mesolevels. Of course this would require the solution of another significant research problem: the modification of self-similar-field micromechanics theory for estimation of composite properties to account for matrix viscoelasticity.(3)Another aspect which was not accounted for in detail in the current implementation was the properties of the fiber-matrix “interphase,” since perfect adhesion was assumed at the interface. Of course that would not be the case in real scenarios since a finite interfacial region of steep property gradients would exist between the fiber and matrix which would significantly govern the residual stresses at the mesoscales. This interphase would in turn be significantly impacted not only by the processing conditions, but also by the nature of reinforcement, for example, glass fiber, carbon fiber, or Kevlar fiber, and the interfacial modification (sizing) employed during composite formulation. The multiscale modeling framework offers rich possibilities for defining fundamental investigations centered on this aspect.(4)Going further, the model will be employed to investigate the interplay between the microstructural heterogeneities, such as, resin rich areas, fiber rich areas, and voids at different locations along the thickness of the composite, and the globally imposed processing conditions, in the development of residual stresses. The generated information from this modeling framework would help in identifying potential locations within a thick composite where presence of certain defects might be most detrimental to the strength of the composite.
Conflict of Interests
The authors hereby declare that they have no conflict of interests.
Acknowledgments
The authors acknowledge motivating discussions with Hamid G. Kia, N. Ramakrishnan, Arun M. Kumar, and Rajesh S. Raghavan during the course this work. Constructive feedback from the reviewers of this paper is also appreciated.