Mesoscale Process Modeling of a Thick Pultruded Composite with Variability in Fiber Volume Fraction

Pultruded fiber-reinforced polymer composites are susceptible to microstructural nonuniformity such as variability in fiber volume fraction (Vf), which can have a profound effect on process-induced residual stress. Until now, this effect of non-uniform Vf distribution has been hardly addressed in the process models. In the present study, we characterized the Vf distribution and accompanying nonuniformity in a unidirectional fiber-reinforced pultruded profile using optical light microscopy. The identified nonuniformity in Vf was subsequently implemented in a mesoscale thermal–chemical–mechanical process model, developed explicitly for the pultrusion process. In our process model, the constitutive material behavior was defined locally with respect to the corresponding fiber volume fraction value in different-sized representative volume elements. The effect of nonuniformity on the temperature and cure degree evolution, and residual stress was analyzed in depth. The results show that the nonuniformity in fiber volume fraction across the cross-section increased the absolute magnitude of the predicted residual stress, leading to a more scattered residual stress distribution. The observed Vf gradient promotes tensile residual stress at the core and compressive residual stress at the outer regions. Consequently, it is concluded that it is essential to take the effects of nonuniformity in fiber distribution into account for residual stress estimations, and the proposed numerical framework was found to be an efficient tool to study this aspect.


Introduction
Composite materials' structure governs their effective properties, and the structure is governed by the processing steps during manufacture [1]. Here, the fiber arrangement impacts several properties of the composite part (e.g., electrical, thermal, and mechanical) as well as process-induced effects. The mechanical performance, which is considered a mechanical property, is directly influenced by the process-induced residual stress [2,3]. Therefore, the coupled nature of processing, structure, and residual stress gives critical information about composite materials' performance. In this paper, the relationship between the meso-structure (associated with the fiber arrangement) and residual stress in a pultruded profile is the main focus.
Pultrusion is an increasingly popular production technique to manufacture fiberreinforced composite profiles [4]. The process is a highly efficient way to produce composite profiles due to its continuous and automatized nature. In a pultrusion process, the resin impregnated fiber reinforcement is pulled through a heated die to produce constant crosssectional composite profiles [5].
Pultruded profiles are subject to residual stresses due to non-uniform heating/cooling, non-uniform fiber distribution, cure shrinkage, coefficient of thermal expansion mismatch, elastic moduli evolution gradient through the cross-section, resulting in geometrical distortions [6][7][8], and pre-mature cracks [9]. The locked in residual stresses should be taken into account to prevent unprecedented failure or undesirable product rejection due to geometrical distortions [2]. Since the cure reaction takes place while the profile is moving during the pultrusion process, a specialized numerical framework is needed to analyze pultrusion. This is a subject that has been studied extensively in the literature.
In one of the earliest works, a thermo-chemical model of pultrusion was carried out with a 1D Lagrangian model that followed the profile throughout the process [10]. Later, full 2D and 3D methods were developed for thermo-chemical modeling [11][12][13]. The impregnation flow of the resin system has also been coupled with thermo-chemical models [14]. This analysis is essential for resin injection pultrusion. In [15,16], 2D and 3D mechanical models in the Lagrangian frame were coupled with a 3D thermo-chemical model in the Eulerian frame to predict the residual stress formation in pultrusion. Accordingly, in a recent paper, a fully coupled 3D Eulerian model was developed [17]. Using these approaches, numerous studies have been carried out in the literature to predict the residual stress in pultruded components with different cross-sectional shapes, reinforcement configurations, and for different process parameters [6,7,9,15,18,19]. These numerical models showed that the residual stress locked in a pultruded profile could reach critical levels for the failure initiation and the profile's ultimate strength.
The process models of process-induced stress in pultrusion have been validated experimentally. For example, the spring-in of an L-shaped pultruded profile and warpage in a hollow pultruded profile were analyzed numerically and experimentally in [6,7], respectively. In later research, Ref. [20] reported the spring-in values of L-shaped profiles produced with different pulling speeds. In addition to shape distortions, experimental validation of process-induced stresses has been conducted as well. The tensile residual stress at the core of a thick square profile was validated via hole drilling experiments in [21,22]. Here it was found that the measured and the predicted residual stress values can reach up to 10-15% of the estimated transverse tensile strength. Therefore, it was argued that it is essential to consider the residual stress in the initial design steps.
The structural variability should also be considered in the design steps as V f is one of the most prominent parameters defining the performance of a fiber-reinforced polymer composite (FRPC). The structural variability in FRPC can be found in the form of fiber waviness, misalignment, wrinkling, intra and inter-tow resin-rich regions, etc. [23][24][25][26]. Especially in thick composites, consolidation results in ply thickness variability, which causes fiber volume fraction (V f ) variability in the through thickness direction [27]. The internal variability in fiber distribution has also consequences on residual stress and processinduced deformations. For example, in [28], a V f gradient was introduced by a resin bleeder, which resulted in spring-in of a composite laminate. The effect of fiber wrinkling on the spring-in value was shown in [29]. Pultruded composites have structural variability as well. The V f values from different locations within I-shaped pultruded profiles were observed via optical microscopy and image processing in [30], and the fiber distribution was reported as neither uniform nor totally random. Local and global variability of fiber volume fraction in pultruded profiles differs for the reinforcement type or geometry of the profile [31,32]. In the resin bath pultrusion, the resin impregnated rovings are consolidated by a solid die and the excess resin is bleed out via squeeze flow [33]. This peripheral compaction, as schematically shown in Figure 1, can be hypothesized as one of the sources of V f gradient through the cross-section. In addition to this global V f variability, a micro/meso-scale variability in V f is inevitable due to intra-roving resin-rich layers and randomly distributed filament within the inter-roving regions. Detailed analyses of fiber misalignment, fiber mat misalignment, resin-rich areas, and void topology in pultrusion products have been performed via optical micrography, Fourier transformation, and computational tomography [34][35][36]. In addition to the characterization of the variability in fiber distribution, its effects on a pultruded bridge deck profile's ultimate properties/performance were investigated in [37]. The effect of non-uniform fiber distribution on the cure degree and temperature evolution through the process was studied numerically in [38]. Local V f variability has been implemented in using a pixel-based finite element method [39]. Here, the elastic material properties within a selected area were defined with respect to each finite element's corresponding V f . In our previous work [40], we have developed a simplistic model and investigated the effect of V f nonuniformity through the cross-section of a pultruded profile. There has been limited studies to unravel the influence of local V f distribution on the residual stresses for pultruded thick composites. It is a scientific challenge that requires the actual V f distribution of the whole cross-section of a thick pultruded composite and implementation of location dependent material properties with a proper RVE selection. Indeed, the microstructural variability through the whole structure is mostly overlooked in the literature, which leads to ruling out the effects caused by the processing-structure relationship. Moreover, to the best of the author's knowledge, the correlation of V f distribution with the local residual stress evolution has not been studied yet.
In the present work, the V f nonuniformity in a unidirectional (UD) fiber-reinforced pultruded profile, which is a consequence of compaction at the entrance of the die as hypothesized, is analyzed in depth, and a mesoscale process modeling framework, which takes the nonuniformity into account, is demonstrated. we build upon our previous outcomes and improve our model by introducing; (i) the effect of patch size on the temperature, cure degree, and local variations of residual stress and its magnitude, (ii) location dependent residual stress concerning local V f . We address these points by doing the following. First, the fiber distribution across the whole cross-section of a relatively thick pultruded profile is presented. The V f dependent thermal, chemical, and mechanical material property evolutions are calculated and implemented into a 2D Lagrangian thermo-chemical-mechanical process model. The difference in predicted temperature and cure degree fields at the exit of the die for different RVE sizes representing the fiber arrangement is presented. The evolution of residual stresses on various locations through the cross-section is inspected to reveal the influence of V f on the predicted residual stress. Location dependent residual stress distribution for the corresponding V f values is analyzed in depth. The critical regions under different loading conditions are discussed concerning processing-structure-property relationships.

Materials and Methods
A unidirectional glass fiber-reinforced polyester-based pultruded profile (19.5 × 19.5 mm 2 ) was investigated in this study. Samples were cut using a diamond saw with water cooling. Before embedding the sample into cold mounting epoxy, they were kept in a vacuum oven overnight at 45 • C. The samples were ground and polished using 500, 1000, 2000, 4000 grit size grinding papers and 9, 3, 1 µm solutions in an automated polishing machine 'Struers Tegramin-30' (Struers, Cleveland, OH, USA). An optical microscope, 'Keyence VHX 1000' (Keyence, Osaka, Japan), was used to capture the micrographs. Ring illumination was used during micrography as it results in higher contrast due to the opaque resin system and transparent glass fibers. Images were taken with 400× magnification with automatic stitching. The stitched images taken from the microscope were stitched again manually to yield a single image of the whole cross-sectional area of 19.5 × 19.5 mm 2 . The investigated cross-sectional area can be seen in Figure 2a.
The V f distribution through the cross-section was calculated using an image processing tool in Matlab (R2019b). First, individual fibers were detected via the circle finding function in Matlab (R2019b) , 'imfindcircles'. Subsequently, overlapped circles were identified and eliminated using an in-house code. The V f was calculated from the area fraction of fibers to the total area. The average fiber volume fraction values were calculated for different-sized RVEs throughout the cross-section. A complete cross-sectional view of the sample together with RVEs of the selected dimensions can be seen in Figure 2. The labels used to represent each RVE size case and the corresponding number of patches through the cross-section together with the RVE edge lengths can be seen in Table 1.

Numerical Modeling Framework
A 2D Lagrangian modeling framework was used to solve the coupled thermo-chemical-mechanical problem in the present study. In this modeling approach, a 2D material frame, representing a cross-section of the profile, moves along the pulling direction starting from the die entrance until the profile is cooled down. Similar to [10], the heat transfer in the fiber direction was neglected. This simplified 2D model allows using a fine mesh discretization, as it yields a fast computation time. This is essential to implement the V f variability with a high resolution throughout the cross-section, which is why this approach was chosen. A schematic representation of the modeling framework can be seen in Figure 3. As seen in Figure 3, a temperature boundary condition was implemented on the outer surfaces of the profile. In other words, the equations were only solved for the composite material. The 2D transient heat conduction equation is given in Equation (1).
where ρ c is the density, Cp c is the specific heat capacity of the composite, T is the temperature, t is time, X and Y are the in-plane directions, and k X = k Y are the thermal conductivities. Finally, q is the heat source term coming from the exothermic heat reaction of thermoset resin. A 2D generalized plane strain formulation was used to estimate the displacement field of the cross-section. The incremental mechanical strain and stress components followed Equations (2) and (3).
where ∆ mech ij , ∆ tot ij , ∆ therm ij and ∆ chem ij are the mechanical, total, thermal and chemical strain increments, respectively. ∆σ ij is the incremental stress component and J is the Jacobian matrix. The entries in J can be seen in Equation (23) in [15].

Chemical and Thermo-Mechanical Material Models of Polyester Resin
The resin's chemical structure evolves via curing, which results in an exothermic reaction, where the physical state of resin changes from a liquid to a rubbery and then a solid state. Therefore, the temperature and cure degree dependent constitutive behavior of resin is essential in process modeling. The resin system used in this profile was characterized in [41]. The corresponding material models, namely the cure kinetics and CHILE (cure hardening instantaneous linear elastic) models, are briefly presented in this study. Readers are referred to [41,42] for detailed information about the characterization of the material systems used in pultrusion.
The cure kinetics was modeled using an autocatalytic model with an Arrhenius type temperature dependency. Therefore, the cure degree (α) evolution and heat generation terms were captured by Equation (4).
where A 0 is the pre-exponential constant, E a is the activation energy, m and n are the reaction orders. α is the instantaneous cure degree, R is the gas constant, and T is the absolute temperature. The total exothermic heat reaction was measured as 175 ± 15 kJ/kg. The cure kinetics model parameters can be seen in Table 2. A CHILE model was used to implement the cure degree and temperature dependent elastic modulus evolution in the process model. The model used in the present study is shown in Equation (5) as: where T C1 , T C2 , T C3 and E 0 r , E 1 r , E ∞ r are the critical temperatures and the associated elastic moduli, respectively. A e and K e are the constants describing an exponential rise in stiffness during cure. T * is equal to the difference between the instantaneous temperature and the glass transition temperature of the resin (T * = T g − T). The instantaneous glass transition temperature is modeled with a linear relationship to the degree of cure from 0 to 135 • C (T g = 135 × α). The constants used in the elastic modulus model can be found in Table 3. Table 3. Parameters used in the CHILE model to characterize the instantaneous elastic modulus of the resin system (reproduced from [41] with the permission). Other than cure kinetics and the elastic modulus evolution of resin through the process, the thermal expansion and cure shrinkage behaviors of the resin system are also essential to predicting process-induced stresses. As shown in [43], the coefficient of thermal expansion (CTE) value of a thermoset resin is generally higher above T g than below T g . In this study, the CTE of polyester resin was set to 72 ppm/ • C below T g . Accordingly, the CTE in the rubbery state (above T g ) was set to 180 ppm/ • C, which is 1.5 times higher than the CTE in its glassy state. In [7], experimental and numerical analysis of an L-shaped pultruded profile produced with a similar resin system showed a good agreement with the experimentally measured spring-in value when a 6% of volumetric shrinkage was chosen. Therefore, the total volumetric shrinkage value was defined as 6% in the present study.

Thermal and Thermo-Mechanical Constitutive Material Properties
The thermal and mechanical constitutive material properties of the composite material were calculated based on the corresponding fiber volume fraction and constituents' properties (the glass fiber and polyester resin). The thermal and mechanical properties of the glass fiber and polyester resin used in this study are presented in Table 4. Table 4. Thermal and mechanical material properties of the glass fiber and resin system used in this study (reproduced from [7] with the permission).
Polyester resin 1100 1830 0.17 Density, specific heat capacity, and thermal conductivity are the required properties to solve the thermal equations. The density, the specific heat capacity, and the conductivity of representative volume elements were calculated via rule of mixture as a commonly accepted and a simplistic approach for process modeling of composites. The thermal conductivity transverse to the fiber direction was calculated using the relationship in Equation (6) [44].
where k f and k r are the thermal conductivity of the fibers and resin, respectively. ω f and ω r are the weight fractions of the fiber and resin. The mechanical constitutive behavior, the equivalent thermal expansion/shrinkage, and the chemical shrinkage are required to estimate the evolution of process-induced stress.
Here, the V f dependent elastic behavior of the composite as well as the thermal expansion and cure shrinkage coefficients were captured via the self-consistent field micromechanics (SCFM) by Bogetti and Gillespie [45] which is a commonly used micromechanics approach for continuous UD fiber-reinforced composites.

Coupled Temperature Displacement Model
The coupled thermo-chemical-mechanical model was developed in Abaqus. As shown in Figure 4, a temperature boundary condition was applied to the outer nodes of the 2D Lagrangian frame that was moved along the pulling direction. The die was considered to be rigid. The interaction between the die and the composite was defined as frictionless in the tangential direction. It was defined as hard contact in the normal direction. The separation between the profile and the die was allowed, which prevented the suppression of cure-related shrinkage. After exiting the die, which corresponds to 400 s after the initial entrance for a 150 mm/min pulling speed and a 1 m die, cooling of the profile was simulated as convective heat transfer towards the outer surroundings. The ambient temperature was defined as 20 • C, where the convection coefficient was defined as 10 W/(m 2 K). The simulation was stopped when the profile was cooled down to 20 • C.
A set of subroutines were used to implement the material behavior into the process model. The elastic modulus evolution was defined in UMAT. USDFLD was used to import the V f distribution into the ABAQUS model. UEXPAN was used to define the V f dependent coefficient of thermal expansion and cure shrinkage. The heat generation through exothermic cure reaction depending on V f was included via HETVAL. V f dependent specific heat capacity, thermal conductivity, and density were defined in the user interface in the form of tabular data in which the V f was taken as a field variable. In total, 4080 elements were used in the model, which corresponds to approximately 64 × 64 elements through the cross-section. The mesh structure used in this study can be seen in Figure 4. The total number of elements was same for each RVE size case, i.e., the mesh density was same for each model. The element type was set to '8-node generalized plane strain thermally coupled quadrilateral, biquadratic displacement, bilinear temperature, hybrid, linear pressure, reduced integration' (CPEG8RHT).

Fiber Distribution through the Cross-Section
The estimated V f distribution maps for different RVE sizes are depicted in Figure 5. The average fiber radius is 14.5 µm with a standard deviation of 1.5 µm. The smallest RVE edge length used in this analysis is 312.5 µm, 22 times larger than the fiber radius. In accordance with the literature, this RVE size is large enough to represent the local mechanical properties [46]. The V f distributions show that there is a V f gradient through the cross-section, with increasing fiber content towards the center of the profile. A similar observation was reported before in the literature for a profile manufactured using the resin injection pultrusion process [38].
The average V f value through the cross-section of the profile is 58%. The overall span and the variation of the corresponding V f values for RVEs increase with the decreasing RVE size [47]. The minimum V f values for 'Case E' and 'Case A' are 21.0 and 54.0%, respectively. The maximum V f values are 89.4 and 60.5% respectively for the smallest and the largest windows. The coefficient of variation reduces from 11.48 to 3.59% for 'Case E' to 'Case A'. The distribution of V f values and the coefficient of variation for different RVE sizes can be seen in Figure 6.

Effect of Nonuniformity on the Temperature and Cure Degree Evolution
In a conventional pultrusion process, the profile is heated from the outside via conduction from the pultrusion die. Consequently, the temperature rises at the profile's outer regions first and then in the core region. In addition to the heating die, internal heat generation occurs during the exothermic cure reaction, which affects the thermal history, particularly in the core of the profile. In the thinner profiles, the temperature distribution can be almost uniform throughout the cross-section [15]. On the other hand, for thicker profiles, it is possible to have a more non-uniform temperature distribution due to the exothermic heat reaction. With the internal heat generation's help, the temperature at the center of the profile can overshoot and reach a higher temperature than the die temperature. This macro-scale temperature nonuniformity is related to the process parameters and exists profiles with a uniform fiber distribution [15].
The thermal conductivity of glass fiber is higher than the conductivity of the polyester resin system. It is vice versa for the specific heat capacity. With increasing fiber content, the equivalent thermal conductivity and specific heat capacity of an RVE increases and decreases, respectively [44]. Consequently, the higher fiber content at the profile's outer regions enhances the heat transfer from the die to the core of the profile. On the other hand, a higher fiber volume content means lower exothermic heat generation due to lower resin content. These phenomena influence the overall thermal history of the profile during the process. Temperature distributions through the cross-section at the exit of the die for different window sizes are shown in Figure 7. At the die exit, the maximum temperature of the uniform case and 'Case E' were 145.9 and 149 • C, respectively. Throughout the process, the maximum temperatures observed in these extreme cases were 146.1 and 149.6 • C, respectively. Therefore, it can be said that the experimentally observed nonuniformity has a negligible effect on the temperature distribution. Temperature evolution through the process is strongly affected by the mold temperature, geometry, resin kinetics, pulling speed. Still, it can be said that the maximum temperature values predicted in this work were in the range of the reported values (∼140-170 • C)in the literature for polyester resin [7,48,49].
Cure degree evolution throughout the process depends on the thermal history as given in Equation (4). Here, a higher heat generation due to lower fiber content amplifies the local cure degree evolution. The cure degree distributions at the exit of the die can be seen in Figure 7. The nonuniformity had a more prominent effect on the cure degree distribution than the temperature distribution. The RVEs with relatively higher fiber content concerning the adjacent RVEs at the core region had a lower cure degree. This phenomenon can be attributed to the difference in local internal heat generation. It is worth noting that the minimum cure degree value within the profile for each case was 0.99 at the end of the process (after cool-down). At this cure degree value, the profile can be considered to be fully cured.

Effect of Nonuniformity on Global Residual Stress Distribution
Residual stress formation in a composite material is a consequence of the thermochemical-mechanical history of the material itself and external sources. Therefore, the residual stress field reveals the combined effects of the temperature, cure degree, and elastic modulus evolution through the process. The predicted stress fields in the transverse plane for the process conditions presented in this study are shown in Figure 8. A similar trend was reported in the literature for square profiles as tensile residual stress in the core and the compressive at the outer sections [15,16,49]. In each RVE size case in this study, the core region was mainly in tension, which was equilibrated with compressive stresses at the outer region. However, the local stress values were scattered with decreasing window size. In addition to the observed scatter, local maximum and minimum stress values within the profile also increased with decreasing window size. The maximum tensile residual stress value observed in the 'Case E' was 20.58 MPa, whereas it was 3.95 MPa for the uniform case. Similarly, the maximum compressive stress predicted within the cross-section was −43.53 MPa, whereas it was −8.86 MPa in the uniform case. The increasing trend of the stress values with decreasing RVE size can be seen in Figure 8. The boundaries of the core and outer regions were defined to compare the average stresses between different cases. The area where the stress was higher than 1.5 MPa in the uniform case was labeled as the "core region". The area where the stress was lower than −2.0 MPa was labeled as the "outer region". The corresponding core and outer regions are depicted in Figure 8. These areas were defined separately for the X and Y directions. The average stress and the V f values in the corresponding areas for different patch sizes are summarized in Table 5. The corresponding average V f values can be seen in parentheses.
The average V f values relatively low for the core, and high for the outer regions. Together with the temperature and cure degree nonuniformity, which is mainly caused by the nature of the process as seen in Figure 7, relatively higher V f at the outer region has an amplifying effect on the residual stress. This is because a higher V f increases the elastic modulus, resulting in higher stiffness at the outer region during the process [45]. On the other hand, relatively lower V f at the core region means higher thermal contraction during cooling and higher cure shrinkage due to polymerization of unsaturated polyester resin. Due to the outside-in heating of the profile, the core region cures later than the outer regions. Therefore, the stiffer outer shell and higher thermo-chemical contraction in the non-uniform V f cases result in higher average tensile stress at the core compared to the uniform V f case. As tabulated in Table 5, the average V f decreased from 0.58 to 0.55 from the uniform case to non-uniform cases. It is worth pointing out that the predicted average residual stress at the core region showed a sudden increase from the uniform case to 'Case A'. The average tensile stress value at the core was almost equal for each non-uniform V f case. This sudden increase and stabilization after 'Case A' in the average predicted stress indicates the effect of global nonuniformity on the average tensile stress in the core. The compressive average stress values within the outer regions were converged after 'Case B', where the V f value is stabilized too. Moreover, slight differences between the average stress values in X and Y directions were caused by the asymmetry in the V f distribution with respect to the center axis.

Local Residual Stresses
The stress evolution within the selected RVEs was analyzed in detail to discuss the effect of local variability in V f on the residual stress formation. Two RVEs from the core region and two RVEs from the outer region were selected as the comparison cases. One out of these two RVEs at the core and outer regions is in a relatively highly packed position (high V f ), and the other one is in a position with relatively lower V f . They were picked up from the locations closer to each other to minimize the location dependent processinduced stress difference between them. These RVE locations on the cross-section and the corresponding stress evolution in the X and Y directions within these RVEs are shown in Figure 9. The stress evolution curves at the same integration points are plotted for each case to reveal the effect of RVE size on the local stresses.
The first two rows in Figure 9 show the stress evolution within the RVEs located in the core region. At the core, process-induced stresses are supposed to be in tension. The V f of the first RVE was 0.22 for 'Case E' due to a local resin pocket, as shown in the corresponding micro-graph. The V f of the second RVE was 0.79 for 'Case E', which is a highly packed RVE. For the uniform case, integration points in both RVEs had similar stress in both directions, which were around 3.3 MPa. Stress evolution in the RVE with the lower V f shows that the predicted stress value increases with decreasing patch size. Still, there were exceptions such as 'Case A' and 'Case B'. The predicted stress value for the former case was higher than the latter despite the larger RVE size. This can be explained with the estimated V f values of these RVEs. Another critical aspect in the first two rows of Figure 9 is the difference between the predicted stress values in the X and Y directions. This difference was related to the anisotropy in the neighbor RVEs' V f values. On the other hand, stress evolution in an RVE with relatively high V f shows that the predicted stress values may even revert to compression as the patch size decreases. As a general explanation of the stress evolution within these two RVEs in the core region, a relatively low V f value results in enhancing the final tensile residual stress locally. On the other hand, a higher local V f value introduces compression to the final stress state.
Similarly, the local stress evolution at the outer region is shown in Figure 9. For the uniform distribution case, both integration points present at the depicted locations had almost zero stress in the X-direction since they are close to a free edge. In the uniform case, both have compressive stress in the Y-direction, which can be considered to be the equilibrating component of the tensile stress at the core. Stress evolution trends with decreasing RVE size in Figure 9 showed similar behavior as in the core. In other words, in a highly packed RVE, the predicted stress value increased in compression with decreasing window size.
On the other hand, stress evolution of the RVE shown at the bottom row in Figure 9, which has relatively lower fiber content due to a resin-rich layer, turned out to be in tension.
The results indicate that relatively lower V f alters the stress value in tension, and higher V f alters the local stress in compression, in general. Stress values mainly differ during cooling, as shown in Figure 9. Therefore, the relation between V f and the process-induced stress can be attributed to the thermal expansion coefficient, which is directly related to V f and dominant on the thermal strain components during cooling.
The predicted in-plane principal stress values are shown in scatter plots for each case in Figure 10 to reveal the correlation between V f and the predicted values throughout the cross-section. The general trends in stress distribution with respect to V f are similar for each case. However, a broader range of V f values for smaller RVE sizes results in higher stress values. The limit values for 'Case E' were around 24 MPa and −38 MPa, which corresponded to 21% and 89% fiber contents, respectively. The measured transverse tensile strength for glass fiber-reinforced polyester pultruded profiles are at a level of 50 MPa, and the compressive strength values are at a level of 70 MPa in the literature [50,51]. The predicted values are near critical levels. Therefore, the residual stresses should be taken into account when evaluating the performance of pultruded profiles. Moreover, here we have a location dependent V f gradient instead of a random distribution, which would play an essential role in the critical locations depends on the loading scenario.

Processing-Structure-Property Relationship
Investigating the reasons behind the V f gradient throughout the cross-section of this particular pultruded profile is not within the scope of this study. Still, this formation can be claimed to be formed because of the processing steps. Compaction from the outer edges at the beginning of the die can be considered one of the possible reasons for this gradient. In the end, this experimentally observed V f gradient promotes the tensile residual stress at the core due to relatively lower V f values. Similarly, this gradient promotes the compressive residual stress at the outer regions due to relatively higher V f as explained for the local stress evolution in Figure 9. From a structural aspect, the V f value affects the failure behavior of composites. Therefore, superposed effects of the location dependent residual stress and loading-induced stresses, both being related to the location dependent V f , would be decisive on the mechanical performance.
The predicted maximum and minimum in-plane principal stress distributions on different locations through the cross-section, namely the core, neutral and outer regions, are depicted in Figure 11 for 'Case E' to emphasize the processing-structure-property relationships. The previously defined core and outer regions, shown in Figure 8, were defined separately for X and Y components to see the direction dependent average stress predictions. These regions were united as depicted in Figure 11 for the in-plane principal stress analysis. The core region is colored red, the neutral region is colored green, and the outer region is colored blue in Figure 11. The slopes of the mean lines of maximum and minimum in-plane principal stresses were similar for the core and the neutral region, as shown in Figure 11a,b,d,e. However, the slope differentiated in the outer region as depicted in Figure 11c,f. The mean line comparisons in Figure 11g,h clearly show that the average V f value and residual stress values deviated for different regions. The critical regions for tension are primarily located in the core region, and the critical regions for compression are mostly located in the outer regions. Figure 11. Maximum in-plane principal stresses vs. V f for the core (a), the neutral (b) and the outer regions (c). The minimum in-plane principal stresses vs. V f for the core (d), the neutral (e) and the outer regions (f). Best linear fits and confidence intervals of the maximum (g) and the minimum (h) in-plane principal stresses.

Conclusions
In this study, we investigated the V f distribution in a thick UD fiber-reinforced pultruded profile and its effects on the residual stress formation. These were done by characterizing the cross-sectional V f distribution with optical light microscopy and using the local V f dependent material behavior in a 2D numerical process model developed for mesoscale thermo-chemical-mechanical modeling of pultrusion.
The optical light microscopy analysis revealed that the profile had a lower fiber content in the core compared to the outer region. When we used this information in the mesoscale model, the magnitude of the overall residual stress in the core (tension) and the outer region (compression) increased more than 50%.
By implementation of the full field V f nonuniformity in the mesoscale, we also captured the scatteredness in the residual stress field. The stress field appeared more scattered when the RVE size became smaller. In some areas, the stress values even shifted from tension to compression, but overall, the local regions with high fiber content appeared to exhibit more compressive residual stress.
We exemplified the importance of analyzing the processing-structure-property relationship for a pultruded composite profile in this study. In conclusion, the results underline that the coupled nature of the local fiber volume fraction and the residual stress formation can be of critical importance when accessing the structural integrity of pultruded profiles. As, a huge effort is put on simulating the structural performance of pultruded composites [52]. We shed a light on the mesoscale nonuniformity in Vf and its effects on residual stress field, which can be further used for the structural analysis in the presence of this mesoscale nonuniformity.
The proposed image analysis and coupling with the thermo-chemical-mechanical model might not be directly used for pultruded composites consisting of continuous filament mat or woven fabrics. Hence, a 3D microstructural analysis and a modification in the proposed process simulation would be necessary for those types of reinforcements.
The mesoscale process modeling scheme developed in this study can also be used in a multiscale process modeling framework in future studies. The location dependent spatial fiber distribution can be implemented into advanced microscale methods to estimate the fiber distribution dependent thermal [53] and mechanical [46] properties which can be used in a multiscale process modeling scheme. Location dependent coupling between micromeso-macro levels in a multi-scale analysis would give valuable insight into pultruded profiles' structural limits.

Data Availability Statement:
The data presented in this study are available on request from the corresponding authors.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: