A preform deformation and resin flow coupled model including the cure kinetics and chemo-rheology for the VARTM process

The present paper deals with preform deformation and resin flow coupled to cure kinetics and chemo-rheology for the VARTM process. By monitoring the coupled resin infusion and curing steps through temperature control, our primary aim is to reduce the cycle time of the process. The analysis is based on the two-phase porous media flow and the preform deformation extended with cure kinetics and heat transfer. A novel feature is the consideration of temperature and preform deformation coupled to resin viscosity and permeability in the VARTM process. To tackle this problem, we extend the porous media framework with the heat transfer and chemical reaction, involving additional convection terms to describe the proper interactions with the resin flow. Shell kinematics is applied to thin-walled preforms, which significantly reduces the problem size. The proposed finite element discretized system of coupled models is solved in a staggered way to handle the partially saturated flow front under non-isothermal conditions efficiently. From the numerical example, we conclude that the cycle time of the VARTM infusion process can be shortened over 68% with the proper temperature control. Moreover, the proposed framework can be applied to optimize the processing parameters and check the compatibility of a resin system for a given infusion task.


Introduction
Nowadays, the vacuum-assisted resin transfer molding (VARTM) process is extensively used to produce fiberreinforced polymer composites. It has the benefits of the enclosed volatile organic compounds emission, high flexibility, scalability, and affordability [1]. The VARTM process mainly consists of two stages, 1) the liquid thermoset resin infusion, and 2) the resin curing stage. The viscous resin is infused into the pre-compacted porous fiber preform by the pressure difference between atmospheric pressure and the vacuum inside the flexible bag. In this stage, void formation and thickness variation are the main challenges. Those are primarily related to various process parameters, such as the arrangement of inlets and vents, What is more, Simacek & Advani [4] and Zhou, Alms & Advani [5] contributed methods of modeling dual-scale flow. To consider the deformation of the flexible vacuum bag, Loudad, Saouab, Beauchene, Agogue & Desjoyeaux [6] reported a multi-layer model. Chevalier, Bruchon, Moulin, Liotier & Drapier [7] integrated capillary actions into a general bi-fluid model. Moreover, Wu & Larsson [8] suggested a shell model that is capable of simulating both the resin flow and the lubrication and relaxation mechanisms of the preform during the VARTM process. These models have reported good agreements with the initial experiments, e.g., [9]. However, the models, as mentioned above, are limited to the isothermal condition. The onset of curing during the infusion stage has therefore been ignored.
To model and simulate the resin curing stage, many researchers have proposed a variety of cure kinetic and chemo-rheological models. Kamal & Sourour [10] and Bogetti & Gillespie [11] were the pioneers who proposed the most frequently used cure kinetics models. In addition, Stanko & Stommel [12] and Tao, Pinter, Antretter, Krivec & Fuchs [13] applied the Model Free Kinetics approach to model the curing under both non-isothermal and isothermal heating conditions. The curing processes of different processes and composite materials have been studied. Ruiz & Trochu [14] presented a reaction kinetics model for the Resin Transfer Molding (RTM) process. Muc, Romanowicz & Chwałl [15] reported optimization of the RTM curing process. Benavente, Marcin, Courtois, Lévesque & Ruiz [16] reported a model to predict the geometrical distortion during cure cycles. Moreover, the cure kinetics model for the autoclave molding process has been reported by Maji & Neogi [17]. As to the out-of-autoclave process, Dong, Zhao, Zhao & Yu [18] developed a fast curing model. Garschke, Parlevliet, Weimer & Fox [19] studied the curing stage of the out-of-autoclave (OOA) toughened resin film infusion (RFI) process. Guo On the other hand, concerning the viscosity, the Castro-Macosko model [23] is one of the most extensively wellknown chemo-rheological models. Researchers propose many chemo-rheological models as well. Henne [29] performed a nonlinear-regression analysis of viscosity and degree of cure data to quantify the dependence of viscosity on temperature and extent of cure reaction. Geissberger, Maldonado, Bahamonde, Keller, Dransfeld & Masania [30] reported a time-dependent chemo-rheological model without the knowledge of the cure kinetics.
The studies mentioned above either aimed at the infusion stage or the curing process independently; however, in fact, the resin flow development, heat transfer, and resin cure are strongly interrelated at both stages. The exothermic reactions increase resin temperature -additionally change the viscosity, which in turn affect the flow filling pattern. Moreover, the degree of cure upon the completion of the infusion stage contributes to the initial conditions for the curing stage. The reduction of the process cycle time may benefit from the adequate control of the temperature and initialization of polymerization. In those regards, Celle, Drapier & Bergheau [31] modeled the non-isothermal fluid flow during the RFI process omitting the changes of resin viscosity. Keller, Dransfeld & Masania [32] proposed a framework to simulate the flow and heat transfer during the Compression Resin Transfer Molding (CRTM) process. Shi & Dong [33] and Yang, Lu, Xian & Liu simulated the filling and curing processes in the non-isothermal RTM process. So far, the holistic simulation of the VARTM process has rarely been reported. Grujicic, Chittajallu & Walsh [34] and Ma, Athreya, Mehta, Barpanda & Shafi [35] simulated the non-isothermal VARTM process of large and thick composite parts. However, the consideration of the multi-phase flow, the capillary effect, and mathematical models of preform deformations are absent in those works. Based on the previous shell model [8], we include the cure kinetics, chemo-rheological model, and heat transfer in this study to simulate both the filling and curing stages in the non-isothermal VARTM process. By virtue of the simplified shell model, we propose a system of coupled models that can efficiently simulate and predict the flow pattern, pressure, temperature, degree of cure, viscosity, and preform deformation in the course of nonisothermal VARTM process for the thin-walled preform. Through the numerical experiment, we also show the potential of initiating polymerization during the infusion stage to optimize the total process cycle time. In the numerical experiment, the thermoset resin and carbon NCF preform are studied initially. Besides, the emerging reactive thermoplastic resin is considered. Because the viscosity of the reactive thermoplastic liquid methyl-methacrylate resin [36] is as low as the thermoset polymer at room temperature, it can be processed directly by the VARTM process. On the other hand, the polymerized thermoplastics can be melted. Based on this characteristic, the numerical experiment discusses the possibility of accelerating the VARTM process using reactive thermoplastic resin under high temperature.
The paper is organized as follows: in "Governing equations", we present the mathematical models, followed by numerical simulation in "Numerical simulation". Furthermore, the results from the simulation are compared and discussed in "Results and discussion". Finally, the concluding remarks are made in "Conclusion".

Governing equations
In this section, we outline the flow model, the energy equation, the cure kinetics, and chemo-rheological models. Concerning the non-isothermal VARTM process, we assume 1. the preform is considered as thin-walled and the through-thickness flow is neglected [8]; 2. the heat transfer between the preform and resin/air mixture fluid are instantaneous [37]; 3. the viscosity dissipation energy in the energy equation is neglected [33]; 4. the gravity is neglected; 5. the infusion is stopped when the mold is filled.

Model of resin infusion in the deformable preform
The flow front position is identified using the concept of saturation degree. Let n s and n f represent the volume fraction of fibers and the mixture fluid of liquid and air, respectively, which satisfies n s + n f = 1 . (1) The mixture fluid is further divided by the saturation degree, 0 ≤ ξ ≤ 1, as where ϕ l := ξn f denotes the liquid volume fraction, and ϕ g := (1 − ξ)n f is the air volume fraction. The homogenized state variables can be expressed following the homogenization approach from [3] where χ stands for the density -ρ, viscosity -μ, thermal conductivity -k, and specific heat capacity -c p . Two equations from the mass balance relation have been derived from [8], i.e., the saturation equation and the pressure equation, In Eqs. 4 and 5, J is the Jacobian of the deformation gradient tensor F , which is defined as where I is the second order identity tensor, N is the normal of the preform top surface, and λ := h h 0 = J denotes the ratio between the deformed preform thicknesses h and the un-deformed preform thicknesses h 0 [8]. The Darcy velocity of the homogenized fluid, v df , is expressed as wherek is the intrinsic permeability tensor of the preform. The relative permeabilities of resin, k l r , and gas, k g r , are obtained from [38], respectively Equations 8 and 9 are derived from the Brooks and Corey's [39] capillary pressure model, where p c is the capillary pressure, the pressure of air is represented by p g and the pressure of resin is denoted as p l . The entry pressure p ent and the parameter n b control the capillary pressure curve.

Fiber preform compaction model
From [8] and [9], we obtain a constitutive relation between the fluid pressure and the fiber preform stretch, where k s E and m are parameters in the Toll's packing law [40], the atmospheric pressure is represented by p a , the volume fraction of the un-deformed fiber preform is denoted as n s 0 . We emphasize that Eq. 11 is based on an exponential hyper-elastic law that describes the fiber packing in terms of the fiber volume fraction n s , cf. [8].
As long as the preform compacts or expands during the VARTM process, n s and n f are varying accordingly, Thus, the permeability variations can be expressed using the Kozeny-Carman equation aŝ wherek ij is the in-plane principle permeabilities, where i = 1, 2, and j = 1, 2.

Cure kinetics and chemo-rheological models
The Kamal and Sourour [10] autocatalytic model is modified with the diffusion factor to improve the cure kinetics at the final cross-linking stage as where the function G denotes the rate of change of degree of cure, c is the degree of cure, T is the resin temperature. c d denotes the diffusion control factor, and c c is the critical degree of conversion, n 1 and n 2 are the overall reaction orders. k 1 and k 2 are the specific rate constants given by Arrhenius expressions, where A i is the pre-exponential factor, E i is the activation energy, and R is the universal gas constant. Given the fact that the cure reaction is moving as the fluid is transporting in the porous fiber preform, the degree of cure for a given temperature, T , depends on both time and positions of the flow, i.e., c := c [x, t]. We derive (16) from the mass balance of the cured resin (see Appendix A) following the theory of porous media as When the preform is fixed, we have v s = 0. Thus (16) reduces to the form reported in [37,41,42] as ϕ l ∂c ∂t the change of polymer (17) where D e is the effective diffusivity tensor, and D D is the hydrodynamic dispersion tensor, which are defined as where D = 8.2 × 10 −12 m 2 /s is the molecular diffusion coefficient of the resin [43], I is the second order identity tensor. where c D are longitudinal and transverse dispersion coefficients, respectively [44]. The Peclet number is defined as P e c = |v dl |d f /D, where d f = 7 μm is the fiber diameter.
From the Castro-Macosko [23] model, we obtain the viscosity as a function of the temperature and degree of cure, where μ ∞ is the viscosity at zero cure and T = ∞, E μ is the activation energy, c g is the degree of cure at gelation, and a and b are the parameters fitted from the dynamic and isothermal rheology tests.

Energy equation
We formulate the heat transfer equation considering the multi-phase flow and the deformable preform as (see Appendix A for detailed derivations) (21) where the homogenized specific heatĉ p = n s ρ s c s p + n f ρ f c f p . If the preform is fixed, the preform velocity turns to zero, v s = 0. Therefore, Eq. 21 reduces tô where T is the temperature, H is the reaction heat of resin, and c denotes the resin degree of cure. k f and k s are the conductivity of the fluid and preform, respectively.  Figure 1 shows the setup of the numerical simulation. The carbon NCF preform is 350 mm in length, 70 mm in width, and 6.5 mm thick. The thermoset resin inlet and outlet are marked by the red and blue lines, respectively. The black edges are impermeable and perfectly insulated, and the green centerline is selected to collect and represent results. Because of the symmetry, we only model half of the preform.

Numerical simulation
Equations 4, 5, 21, 16 and 20 form a coupled initialboundary value problem system. The system is decoupled by a staggered approach (shown in Fig. 2). In each time step, Eq. 4 is first solved for the degree of saturation. The updated degree of saturation is then used to solve the pressure from Fig. 1 The setup of the numerical simulation of the VARTM process. The inlet width is half of the length of the shorter edge Eq. 5. Next, Eq. 21 is capable of updating the temperature, and the degree of cure is solved from Eq. 16. Finally, the viscosity is computed from Eq. 20. Upon the end of the infusion stage, the curing process is solved without (4) and (5).
Equations 4, 5, 21 and 16 are solved by the Galerkin finite element method and stabilized by the streamline upwind Petrov-Galerkin method [3,45], when necessary. The backward Euler method is used for the time integration. The initial and boundary values are set as follows We discretize the preform to a number of four-node bilinear elements. Different numbers of elements and time step sizes are chosen to study the convergence in "Results and discussion". Tables 1, 2 and 3 list the parameters used in the simulation. Specifically, the parameters of the NCF fiber preform are obtained from the experiment in [9]. The parameters of cure kinetics and rheology of the thermoset resin are obtained from [33] and [46]. The temperature dependent heat conductivity and specific heat capacity of fiber and resin are obtained from [32]. The in-house code is developed in C with the Intel MKL PARDISO library and OpenMP library. The meshes are generated in Abaqus CAE, and the post processing is done by ParaView. We also remark that in the present development of the code distribution media has not been implemented.
The mesh size convergence study is based on the difference of state variables between mesh schemes 1 to 4 and the scheme 5, when the time step size is fixed at size-5. The error concerning the spatial discretization em RMS is defined as the root-mean-square (RMS) error as, with i = 1, 2, 3, 4 , (23) where NT is the total number of time steps; χ mesh−i, t denotes the variable of the i-th mesh scheme at the t-th second. Once em RMS s are computed, the normalized error em is computed as RMS,4 with i = 1, 2, 3, 4 .
In Eq. 23, χ is computed as the integration of the variables over the whole preform as where stands for the variable of temperature T , degree of cure c, viscosity μ, saturation degree ξ , and fluid pressure p in turns, B is the domain of the preform. with n = 1, 2, 3, 4 , (26) where NNO equals to the total node numbers of the mesh scheme 5. ς size−n, i represents the value of temperature T , degree of cure c, viscosity μ, saturation degree ξ , and fluid pressure p at the i-th node of n-th time step size. The normalized error et is computed as  Figure 3 shows that the fastest decrease in error happens when the element number changes from 1000 to 2250. As the element numbers increase, the error continues to reduce; however, the rate of change of error reduction becomes smaller and smaller. Thus, we choose the 6250 elements as the mesh scheme in this study since the error is lower than 10%.
The result of the convergence study of the time step size is depicted in Fig. 4. The error reduces faster when denser time steps are utilized, nevertheless the total error reduction is about 1%. That means that the smaller time step size will not contribute very much to the error control. Fig. 3 Convergence of staggered approach due to mesh refinements after 70 seconds of infusion when the initial temperature is 373 K. The time step size is fixed to 1 × 10 −3 seconds Because of this, the 0.01 seconds is chosen as the time step size to save computational effort while still keeping good accuracy.

Process simulation
The exothermic infusion process highly depends on the ambient temperature. To study the effects of the initial temperature, we design five different cases. The initial temperatures, 0 T , start from 298 K, end at 398 K, and we select the temperature every 25 K in between. Figure 5 shows the filling time and the maximum degree of cure for each given initial temperature. When the initial temperature varies from 298 K to 373 K, the filling out time drops 68.5%. On the contrary, the degree of cure increase three Fig. 4 Convergence of staggered approach due to time step sizes after 70 seconds of infusion when the initial temperature is 373 K. The mesh-5 scheme is selected times. In this regard, the infusion process is accelerated by increasing the initial temperature. However, when the initial temperature keeps rising to 398 K, the resin gels at around 39 seconds. At this moment, the thermoset resin can not move forward anymore, and the process stops finally. However, if we consider the reactive thermoplastic resin instead of conventional thermoset resin, e.g., reactive methyl-methacrylate [36], the in-situ polymerized resin can be melted and flow again due to the heat from exothermic chemical reactions. In this regard, we find out that the resin flow through the whole preform much longer than the other temperatures. It takes about 278 seconds to fill out the preform with a very high degree of cure (0.839) at the end, as the semitransparent parts shown in Fig. 5. However, this theoretical hypothesis of the active thermoplastic resin need further experimental validation as a future work. In conclusion, there is a window of the processing temperature that can shorten the time of the infusion stage. If the initial temperature is too low, the effect of the accelerating  process is small; however, if the temperature is too high, the resin gels before it flows through the whole preform. Nevertheless, regarding the advanced reactive thermoplastic resin, the total infusion time is unacceptable. Figure 6 shows the flow front positions during the process at different initial temperatures. When the temperature varies from 298 K to 373 K, curves exhibit similar shape. The higher temperature gives a higher tangent of the curve, i.e., the flow front moves faster. Apart from that, at 398 K, the flow front moves fastest at the beginning, and it significantly slows down after 35 seconds until the resin flow stops at around 39 seconds. This is because the high initial temperature results in a faster curing process and the gelation of resin lead the viscosity to increase significantly. As mentioned before, the stopped resin can flow again under the heating from the chemical reactions of the reactive thermoplastic resin. The dashed line in Fig. 6 indicates the flow front positions of the reactive thermoplastic resin. It is obvious to notice the "elbow" shape at around 39 seconds, which implies the dramatic decreasing of the infusion speed after the polymerization starts.
In the VARTM process, a fundamental problem is the preform deformation. As reported in [9], the relaxation and lubrication mechanisms dominate the fiber preform deformation modes. Figure 7 shows the preform deformation along the centerline when the initial temperature equals to 373 K. The regions nearby the inlet have the largest thickness. Due to the lubrication mechanism, the minimum thickness parts locate around the flow front, and the widths of them overlap the partially-saturated regions. Upon the end of the infusion stage, the initial flatted preform turns into the shape of a wedge.
To study the development of the flow front in the course of the VARTM process, we plot the changes in temperature throughout the filling stage in Fig. 8a. The changes of flow front temperatures are limited small except the case of 398 K. The dashed line shows the rising of the flow front temperature due to the reactions of the reactive thermoplastic resin. The dashed line tells that the rate of change of the temperature changes from high to low at around 130 seconds since the amount of polymer that can be reacted reduces. When the whole preform is filled out, the jump of temperature is about 43 K. This large jump of the temperature accelerates the curing process and slows down the infusion process. Furthermore, the temperature profile of the case 0 T =373 K along the centerline is shown in Fig. 8b. During the process, the temperature profiles show the shape of an upside-down parabola at all four times. What is more, the maximum temperature at a specific time locates in the middle of the length of the filled-out path. This  Similarly, the flow front degree of cure and the profile of degree of cure along the centerline are plotted in Fig. 9. From Fig. 9a, we can notice that the flow front degree of cure develops linearly when the initial temperature is low (298 K, 323 K, and 348 K). When the initial temperature locates at a middle range (373 K), the degree of cure development represents non-linear behavior. At a high initial temperature, e.g., 398 K, the fully non-linear curing process is exhibited as the dashed line shows. By the control of the initial temperature, we can find a suitable process window of the temperate to reduce the cycle time and limit the resin polymerization. On the other hand, for a given initial temperature, 0 T =373 K, the profile of the centerline degree of cure is shown in Fig. 9b. The highest degree of cure occurs at the fully saturated flow front. The curves of the degree of cure represent the "cliff" shape fronts. Ideally, if it is the fully saturated case, the "cliff" will be straight and upright. However, the existence of the partially-saturated flow front in this study makes the steep slop "cliff" observed instead.
Last but not least, the curves of viscosity at the flow front and along the centerline are shown in Fig. 10a and b, respectively. The flow front resin viscosity shows similar linear and non-linear behaviors like Fig. 9a The viscosity increases to the peak very fast when the initial temperature equals to 398 K. The time matches with the "elbon point on the purple line shown in Fig. 6.
Just because of this high viscosity, the flow moves hardly forward in case of the thermoset resin. After this time, the Fig. 10 a The changes in viscosity at flow fronts for different 0 T s during the infusion stage. b The profiles of the viscosity along the centerline when the initial temperature is 373 K at 10, 25, 50, and 70 seconds after the beginning of the infusion temperature keeps rising due to the exothermic reactions, and the viscosity starts reducing for a long time, if the reactive thermoplastic resin is chosen. As the resin fills out the mold, the cure accelerates and leads the viscosity to increase again as the dashed line shows. Figure 10b is similar to Fig. 9b, but the curve at 70 seconds changes from the convex to concave. Then position of the change point corresponds with the peak temperature location in Fig. 8b. The reduction of temperatures gives rise the slower resin polymerization and viscosity rise.
Next, we plot the contour plots from the 2-D simulation in Fig. 11. From Fig. 11a, we can see the temperature at the center of the width is lower than the edges. Furthermore, the differences in the temperature between the center and edges in the width direction are more extensive in the top half part, but very small in the lower half part, where is far away from the inlet. This pattern also shows in Fig. 11b regarding the degree of cure. Moreover, the lower part in Fig. 11bis almost in the same color, i.e., the resin is stacked at the same degree of cure in those regions. The color gradients only exist in the top half part, which means that the rate of change of the degree of cure is higher close to the inlet and lower when it is far from the inlet.

Conclusion
In this study, the coupled resin infusion and curing stages are modeled for the non-isothermal reactive VARTM process. Restriction is made to thin-walled fiber-reinforced polymer composite components, which significantly reduces the problem size. The implementation is thus made based on the shell model [8], extended by heat transfer, cure kinetics, and rheological models to investigate the effects of the exothermic polymerization during the VARTM process. In this development, two-phase porous media is used to model the inplane resin infusion flow. For the flow front, the capillary effect combined with the relative permeability concept were considered to describe the competition between the resin and remaining air inside the fiber preform. In this context, heat transfer and cure kinetics equations are extended to consider the fiber preform deformation. When the highly convective non-isothermal flow of the reactive resin migrates into the preform, the resin curing generates heat that accelerates the process of polymerization. Accordingly, the resin viscosity and the flow velocity increase or decrease depending on the temperature and the degree of cure. On the other hand, the compaction and extension of the preform change the structure and permeability of the preform, which affects the infusion speed. In order to handle the coupled problem, a FE-based staggered solution procedure is considered. The corresponding convergence studies are also carried out.
By applying the VARTM process in different environmental temperatures, we discovered the existence of a temperature window, allowing for the infusion process speed up via the proper temperature control. When the temperature is too low, the infusion stage cannot speed up too much; however, when the temperature is too high, the process is facing stagnation due to gelation. As to the type of resin, emerging reactive thermoplastic resin with low viscosity at room temperature can replace the traditional thermoset resins in the RTM and VARTM processes. If we do so, based on the characteristics of thermoplastic resins, the partially polymerized resin can melt due to heating from chemical reactions. However, even though the process can continue, the infusion process slows down due to the high viscosity. With the aid of the presented framework, we can optimize the process parameters and the design of the mold, such as the initial temperature, the vacuum level, the position of the inlet and outlet. What is more, the proposed method can check the compatibility of a resin system for a given infusion task without relying on the experience and the "trial and erron approach.