A shell model for resin flow and preform deformation in thin-walled composite manufacturing processes

The paper proposes a novel approach to model the in-plane resin flow in deformable thin-walled fiber preforms for liquid composite molding processes. By ignoring the through-thickness flow in large scale thin-walled components, the 3-D resin flow is simplified to an in-plane flow inside the preform by a specialized divergence theorem. Shell kinematics are used to describe the fiber preform deformation, and the compressible flow is modeled in the context of the free surface flow in porous media. For simplicity and efficiency, the normal stretch, which is driven by the internal fluid and applied external pressure, represents the fiber preform expansion and compression. As compared with full 3-D models, the proposed shell model significantly reduces the problem size, while it still represents the primary physical phenomena during the process. The effects of neglecting the through-thickness flow are illustrated in a numerical example that compares the flow for a set of preforms with different thickness. The model is demonstrated from the numerical example of the mold filling in a doubly curved thin-walled fiber preform. Due to the applied vacuum and the consequent resin flow motion, the relevant deformation of the preform is observed.


Introduction
The class of liquid composite molding (LCM) processes has been widely employed for manufacturing fiber reinforced polymer composite materials (FRPCMs) and helps manufacturers to carve out a niche amid the keen market competition. Since the mid-1980s, the automotive industries started to utilize the resin transfer molding (RTM) method to produce high volume production net shape structural components. Then the vacuum assisted resin transfer molding (VARTM) process sprung up in marine, energy and aerospace industries. The VARTM process can reduce the emission of volatile organic compounds, and produce high-quality FRPCM parts with flexible, handy and lowcost tooling. However, the challenges of defects also appear Ragnar Larsson ragnar.larsson@chalmers.se 1 Division of Material and Computational Mechanics, Department of Industrial and Materials Science, Chalmers University of Technology, SE-412 96, Göteborg, Sweden during LCM processes, e.g., dry spots, spring-in, microvoids and thickness variations [1]. The resin flow distribution arrangement dramatically influences the whole filling process, so it is keen to model and simulate the process numerically instead of relying on trial and error physically.
Between the late 80s and early 90s, Chan and Hwang [2] started to solve the pressure distribution in the RTM process based on the Darcy's law, which describes the fluid transfer in porous media. Fracchia et al. [3] employed the control volume finite element method (CVFEM) together with the concept of volume of fraction (VOF) to simulate the mold filling in the RTM process using conforming finite elements. Since then, many results following this approach have been published, e.g., [4,5]. The modified CVFEM, e.g., [6][7][8], have been proposed. The boundary element method (BEM), e.g., Um and Lee [9], the level set method, e.g., Soukane and Trochu [10] or Gantois et al. [11] and the smoothed particle hydrodynamics (SPH) method, e.g., [12,13] are also among the methods of LCM process modeling. What's more, Remacle et al. [14] presented a high-order model using the discontinuous Galerkin method (DGM). Wu and Larsson [15] proposed a homogenized flow model based on the theory of porous media (TPM) to simulate the wet-out RTM process. Another interesting study has been done recently by Dammann and Mahnken [16] who used the phase-field models to model the RTM process.
The research is not just limited to simulate the mold filling flow, but it is also extended to model the fiber preform response coupled with the resin flow. Niaki et al. [17] developed a three-phase integrated flow-stress model. Li and Tucker [18] reported a method to model the consolidation problem. Besides, Wysocki et al. [19][20][21] have modeled the coupled resin flow/perform deformation problems. As an extension, Rouhi et al. [22][23][24][25] presented a series of works based on the continuum mechanics framework.
In contrast to the rigid upper mold used in RTM processes, a plastic membrane is placed on top of the preform in the VARTM process. Thus, the preform deforms due to the difference between the internal and external pressure, which leads to thickness variations in final products. Due to this complexity, the VARTM process is more complicated to model than the RTM process. To simulate the VARTM process, Blais et al. [26] modeled the process as Stokes-Darcy coupled problem. Andriamananjara et al. [27] considered the effect of capillary pressure in their model. Besides, various studies on both resin flow and preform response were reported by researches e.g., [28][29][30][31][32][33]. However, to solve the full 3-D fluid and solid problems, it demands enormous computational efforts. What's more, it is also a challenge to mesh thin parts and keep the elemental skewness, aspect ratio and warpage in proper quality. On the other hand, the process optimization is an interesting topic for discussion. Hsiao et al. [34] built an optimization framework of flow distribution arrangement by using genetic algorithms. From the industry point of view, the speed of the optimization is vital, which relies on the computational efficiency of the process model.
To simplify the LCM process model for thin-walled FRPCMs, it is assumed that the resin flow is confined to the in-plane of the preform, and that the preform deforms solely along the normal of the mold. Consequently, the preform deformation is represented by a normal stretch variable. Based on the packing law [35,36], an explicit solution of the stretch is derived in terms of the fluid and atmosphere pressure. As a result, the LCM process is simplified to the 2-D in-plane resin flow in the 3D deformable preform. Also, in-plane elements are used based on shell kinematics. The resulting shell model for resin flow and preform deformation in thin-walled process applications significantly reduces the number of degrees of freedom, while the primary physical phenomena of the process can still be represented. The model can be applied at the preliminary process design stage to help industries improving both efficiency and quality of the production. Lastly, the validity of assumptions and the capabilities of the model are illustrated through two numerical examples.

Mass and momentum conservations for non-saturated porous media
Following the developments in [15], we model the wet-out process of the RTM based on the theory of porous media. A two phase porous media is considered, which consists of the solid fiber preform phase and the homogenized resin/gas fluid phase. Let n s indicate the volume fraction of solid phase, whereas n f represents the volume fraction of the (homogenized) fluid phase. So that n s and n f are interrelated through n s + n f = 1 , (1a) where the liquid volume fraction ϕ l and the gas volume fraction ϕ g further subdivide the fluid volume fraction n f . We introduce the saturation degree, 0 ≤ ξ ≤ 1, to indicate the case of partial saturation. We thus express the liquid and gas volume fraction in Eq. 1b in terms of n f and ξ as The dry (n f = ϕ g ) and full-saturated (n f = ϕ l ) cases are represented as ξ → 0 and ξ → 1, respectively. As shown in Fig. 1, the dry and full-saturated parts are separated by a thin boundary, which indicates the flow front position where the gradient of saturation degree exists. Using where D s • /Dt =• and D f • /Dt denote material time derivatives related to the solid reference configuration B 0 and fluid reference configuration B f 0 , respectively. Regarding solid particles, we introduce the deformation gradient F and its Jacobian J as where ∇ X is the material gradient operator with respect to the solid phase. From the principle of mass conservation, we find that the solid contents M s = J n s ρ s and fluid contents M f = J n f ρ f conserve aṡ where ρ s and ρ f are the intrinsic densities of the solid and fluid phases, respectively. In particular, the stationarity of M s in Eq. 5a yields n s ρ s J = n s 0 ρ s 0 , whereby the solid volume fraction is governed by n s = J −1 n s 0 , where • 0 represents the variable in the reference configuration.
Moreover, based on the arguments of homogenization, see [15], we obtain the homogenized fluid density as where ρ l denotes the intrinsic density of the liquid (resin) that is assumed as in-compressible; ρ g is the intrinsic density of the compressible gas. Thus the mixture fluid density is considered as compressible.
Considering the fluid mass flux ρ f v f , it can be expressed as a linear combination of the liquid and gas mass fluxes, scaled by the saturation degree, Multiplying v on both sides of Eq. 6, we obtain ρ f v = ξρ l v+(1−ξ)ρ g v; subtracting this by Eq. 7, and introducing the homogenized relative fluid velocity v rf := v f − v, we obtain In view of (2), Eq. 8 can be further elaborated as where the Darcian velocities v df , v dl and v dg are defined as Upon combining the relations (5) and Eq. 6 together with the saturation constraint Eq. 1a, we obtain the total mass conservation of the solid-fluid mixture as (11) which is named as pressure equation. Similarly, as to the liquid mass, e.g., M l = J ϕ l ρ l , we obtain the balance equation aṡ Given the fact that the liquid resin phase is incompressible, Eq. 12 can be formulated as an evolution equation of the saturation degree, see also [22] and [15], (13) Finally, from quasi-static momentum conservation of the mixture porous media, the linear momentum balance yields (14) whereρ = n s ρ s + n f ρ f , and the total stressσ relates to the effective (constitutive) stress σ of the fiber network and the fluid pressure p, viz., Terzaghi effective stress principle, From ref. 15, it implies that the deviatoric parts of total stress and effective stress are same, i.e.,σ dev = σ dev , whereby the total pressurep is obtained as the summation of the fluid pressure p and the effective pressure p e aŝ As to the fluid stress response, it is assumed that the fluid is ideal with negligible shear stress. Thereby the intrinsic fluid stress is represented by the fluid pressure p, which is homogenized between the pressures of the liquid resin and the gas. It follows from the homogenization in [15], that the (mixture) fluid pressure is an interaction in the degree of saturation ξ between the intrinsic liquid and gas pressures (p l and p g ) written as

Preform deformation and resin flow analysis
In this section, we consider the LCM process of a deformable thin-walled preform resting on a fixed lower surface. The preform undergoes the atmospheric pressure p a through a flexible membrane, as shown in Fig. 2a. Given the nature of a thin-walled preform, it is considered as a single director shell surface 0 with the unit normal N as indicated in both Figs. 1 and 2b, see also [37]. Under the pressure loading in Fig. 2a, we shall assume that the preform can compress or expand only along the normal N via the thickness stretch defined as where h is the current thickness of the preform whereas h 0 is the undeformed thickness. It is also assumed that the fluid pressure p and the saturation degree ξ keep constant along the thickness direction. We thereby simplify the 3-D Darcy flow to the in-plane flow in an arbitrary curved surface .

Geometry of the thin-walled pressure loaded preform
In order to define the geometry of the thin-walled pressure loaded preform, its current geometry B is assumed to deform solely with the normal stretch λ in Eq. 18 so that where the initial position vector θ 1 , θ 2 defines material points on the lower surface 0 of the preform. Next, based on the assumption that only normal direction deformation is allowed, the current geometry B[λ] is defined in terms of the normal stretch on 0 as In particular, the initial undeformed geometry B 0 = B [1] is obtained as From Eq. 21, we find that where the co-variant basis vectors (in-plane G α and out-ofplain N ) are defined as where α = 1, 2 and G 3 = G 3 = N, and the curvature of the preform K αβ is defined as Moreover, in Eq. 23 we introduced the contra-variant basis vectors G i (i = 1, 2, 3), which are defined from the identity relationship where G ij is the metric tensor. Finally, the infinitesimal area element dA on the 0 and volume element dV in the B 0 are formulated in the curve-linear coordinates as Likewise, from the current geometry in Eq. 20 and the curve-linear coordinates of the preform, we identify the deformation gradient from the linearization (27) where the co-variant basis vectors are identified as In the following, let us assume that bending effects of the preform can be neglected corresponding to (λN) ,α ≈ 0 and K αβ ≈ 0, see [37]. Consequently, we can obtain g α = where 1 and 1 = 1 + N ⊗ N are the in-plane and 3-D identity tensor respectively. From Eq. 29, it follows that the normal stretch λ is the Jacobian of the deformation gradient λ = J = det F .

In-plane resin flow in thin-walled preform
To describe 2-D LCM wet-out process of the thin-walled preform, we consider the fluid resin flow in B 0 as an inplane flow that is projected to the surface 0 . To this end, the following form of the divergence theorem is considered for the fluid flow where ∇ is the in-plane gradient operator, N is the outward unit normal of the boundary line 0 on the surface 0 (in Fig. 1), and v df is the in-plane fluid Darcy velocity. Please note that the second order tensors and vectors in 0 satisfy the orthogonality properties, e.g., v df ·N = 0, and N ·N = 0, see also [38]. In view of Eq. 30, we therefore obtain the in-plane pressure gradient Grad p as

Fiber packing law
Considering the mechanical behavior of the deformable fiber preform during the LCM process, a major mechanism is the fiber packing induced by applied pressure. Thus, we consider the Toll's model [35], which is an exponential law in terms of the fiber volume fraction n s . In the present model, we limit the analysis to hyper-elastic with the stored free energy defined as where a slight modification from [35] As the fiber volume fraction n s changes, (increasing when J < 1 and decreasing when J > 1), the packing effect is characterized by the internal contact variations, which is also illustrated in Fig. 3. Hence, given the effective stress Fig. 3 The normalized fiber preform free energy curves in terms of Jacobian. When the parameter m is set to 15 and various initial fiber volume fractions n s 0 are chosen, the fiber packing law exhibits that the free energy exponentially increases during compressing, and linearly increases as expanding principle in Eq. 16, we obtain the total pressure asp = p e [J ] + p. Moreover, in view of the kinematic assumptions of the preform in Eq. 29, the normal stretch can be expressed as λ = J = det[F ].

In-plane Darcy flow
As to the Darcy flow, we postulate that the 3-D flow is simplified to an effective 2-D flow by ignoring the throughthickness flow, whose transformation is defined by Eq. 30 as where K α is the permeability of the preform relating to the liquid/gas transportation. As suggested by Gebart [39], in turn, K α is defined in terms of the intrinsic permeabilityk , the viscosity of the liquid resin/gas μ α and the relative permeabilitties k rα as according to Burdine [40] are defined as Regarding the intrinsic permeability tensork, for simplicity, we consider the isotropic permeability represented bŷ k =k1. Here,k is the deformation dependent isotropic permeability, whose deformation dependence is assumed to follow the Kozeny-Carman equation, see also [39]. Thereby, the permeability is obtained in terms of the fiber volume fraction n s aŝ where r s is the particle (or fiber bundle) radius, and C is the Kozeny constant. From that, we can derive the relation between the deformed preform permeabilityk and the undeformed preform permeabilityk 0 aŝ We thus conclude that the in-plane Darcy law for the liquid transportation is expressed as v dl = −K l p ∇p − K l ξ ∇ξ with, K l p = K l and, and for the gas transportation we obtain v dg = −K g p ∇p − K g ξ ∇ξ with, K g p = K g and, K g ξ = K g 1 + ξ log ,ξ p c p c .
Combining Eqs. 10, 39 and 40, we now obtain the homogenized fluid Darcy velocity as where the mixture permeabilities K f p and K f ξ are obtained as

Capillary pressure and universal gas law
The capillary pressure in Eqs. 38 to 40 represents the difference between the gas and liquid pressure, we apply the phenomenological model by Brooks [41] where p ent > 0 denotes the entry pressure and n b > 0 controls the shape of the capillary curve. By combining Eqs. 43 and 17, the intrinsic pressures of liquid and gas are obtained completely in terms of fluid pressure p and saturation degree ξ as As to the compressible gas density ρ g , the universal gas law shows where R is the universal gas constant and T is the absolute temperature.

Boundary value problems
To solve the present process model, we introduce the mass balance relations Eqs. 11 and 13 and the quasistatic momentum balance Eq. 14 in weak form. The weak forms are formulated in terms of the shell kinematics and the coordinate systems introduced in "Geometry of the thin-walled pressure loaded preform". From "Fiber packing law", we also notice that it suffices to represent deformation gradient F and its Jacobian J = det[F ] by the stretch field λ ∈ L 2 [ 0 ]. So the present model aims to solve the fields of fluid pressure p, saturation degree ξ and stretch λ.

Weak form of momentum balance
The weak form of the momentum balance Eq. 14 is equivalent to the principle of virtual work, which is formulated in the reference configuration B 0 and in the absence of the gravitational force as B 0Ŝ : wheret 1 is the nominal traction vector acting on a surface element dA ∈ 0 . In particular, N ·t 1 := −p a is the prescribed atmospheric pressure at the top surface of the preform, as shown in Fig. 2a. We also introduced the symmetric total second Piola-Kirchhoff stressŜ from the pull-back transformation of the total Cauchy stressŜ = J F −1 ·σ · F −t . By implementing the Terzaghi stress principle Eq. 15, we obtain where C := F t · F is the right Cauchy-Green deformation tensor and S is the effective (constitutive) second Piola-Kirchhoff stress. As mentioned in Eq. 29, the deformation gradient simplifies to consequently, the packing law Eq. 32 yields the effective stress response as Moreover, we also obtain the explicit expressions for the right Cauchy-Green deformation C = F t · F and its inverse C −1 from Eq. 48 where the last equality was obtained using the Sherman-Morrison formula.
In view of Eqs. 49 and 50, the virtual work in Eq. 46 is now worked out as Hence, we find that the applied pressure p a balances with the effective and the fluid pressures, p e and p. From this balance relationship together with the constitutive relation Eq. 33, relating the effective pressure p e and the stretch λ, we obtain whereby the normal stretch field can be resolved explicitly from the hyper-elastic fiber packing law as (53)

Weak formulation of mass balance
Let's consider the in-plane flow in "In-plane Darcy flow" on the surface 0 that is bounded by 0 . We introduce the test functions for the fields of fluid pressures and saturation degrees as where we prescribe the fluid pressure and the saturation degree along the boundaries ∂ p 0 and ∂ ξ 0 , which satisfy that ∂ p 0 ∪ ∂ ξ 0 ⊂ 0 . By using the relation J = λ in Eqs. 11 and 13, we thus obtain the weak form of the pressure and saturation equation projected onto the surface 0 as where Q = N · v df , H = N · v dl , which are the in-plane fluxes across the boundary 0 .

FE interpolations
The weak forms Eqs. 55 and 56 are solved by the finite element method, which is stabilized by the streamline upwind/Petrov-Galerkin (SUPG) method, as discussed in [15]. The domain 0 is discretized by bilinear four-node element. The primary variables p and ξ are interpolated as nodal summation forms, where N I I =1,...,NNO are shape functions, and ξ I and p I are the time-dependent nodal saturation degree and fluid pressure respectively. Furthermore, we obtain the discretized in-plane gradient from Eq. 57 as According to the Eq. 53, we can write out the elemental stretch λ el as Hence, the nodal stretch can be calculated from the average of neighboring elements, where NE denotes the number of elements that is adjacent to λ I . Thus, the field λ is interpolated as

Numerical examples
Two representative examples are presented to validate the 2-D assumptions and demonstrate the model. Table 1 lists the parameters used in both examples. At the infusion inlet, the pressure is set to p 0 = 1 atm and the corresponding saturation degree is given as ξ 0 = 1.0. At the outlet, the boundary conditions are determined by p 1 = 3.2 mbar and ∂ x ξ = 0; the initial values are set to 0 ξ = ξ(x, 0) = 0.001 and 0 p = p(x, 0) = 3.2 mbar. Ideally, the initial field 0 ξ is very close to zero, but numerical singularity problems appear when 0 ξ is chosen too small. Thus the current value of 0 ξ has been chosen small enough, while it avoids numerical singularity.

Comparison of the resin flow in thickand thin-walled fiber preforms
To verify the in-plane flow assumption in thin-walled preforms, we studied two different models when the preform is considered fixed with λ = 1. The considered preform is a curved tunnel that starts and ends with a 0.2 meters long straight path to steady the flow; there is an 0.8 meters radius curve in between that introduces an unevenly distributed through-thickness flow as shown in Fig. 4. Following [15], a 2-D model is first made for capturing both the in-plane and through-thickness resin flow of the curved preform in Fig. 4a. Second, a 3-D shell model is built for the same geometry; however; because there is no flow along the transverse and through-thickness direction of the preform, the flow can be considered as 1-D in the longitudinal direction, as shown in Fig. 4b. The infusion starts from the inlet and ends at the outlet as depicted in Fig. 4. In order to study the influence of the through the thickness flow as induced by the curved preform, a global saturation degree measureξ is introduced to represent overall significance of the through the thickness flow effect. The measureξ is defined as where the NNO is the number of nodes, and the ξ i denotes the saturation degree at the ith node. By comparing theξ in the 2-D plane and shell models, we thus measure the importance of the through-thickness flow for various thickness to length ratios, t/L, of the preform. The thickness t varies among 0.05, 0.1, 0.15 and 0.2 meters. For each t/L case, the global saturation degreẽ ξ is computed and compared between two models. For the very thick preform t/L = 12%, Fig. 4a shows the unevenly distributed through-thickness flow for the 2-D plane model, i.e., due to the curved geometry the inner radius flow is faster than the outer radius flow. Figure 5 shows that the percent error, δ = ξ plane − ξ shell /ξ shell × 100%, decreases linearly as the preform thickness decreases. When the ratio t/L is reduced below 6%, the through-thickness flow effect is significantly reduced. For example, a thin preform of t/L = 3% yields a difference of 1.5%, which justifies the in-plane flow assumption made for the present thin-walled case.

LCM process of doubly curved thin-walled preform
To investigate the capabilities of the proposed model for simulating the infusion of a deformable thin-walled preform, we consider a LCM process example as shown in Fig. 6. The black edges including the hole are impervious, in contrast, the blue edges are the resin inlet and outlet. The curvature R and side length l of the surface 0 are 1 meter. The inlet is 0.25 meters wide and locates in the middle of the northwest edge, and the opening has the diameter d = 0.2 meters located at the center.
We choose four paths to represent the results. The first one is the line from the middle of the inlet to the outlet middle, viz., Center; the second path is the "northeast" edge in the Fig. 6, viz., Side; the next one is the line that equally splits the partition between the Center and Side, viz., Middle; and last, the "northwest" edge is named as Front.
The mesh size convergence study is based on the difference of the global saturation degreeξ of the coarse, regular and fine meshes relating to the finest mesh. The time step t = 1 × 10 −3 seconds is selected and the infusion ends at 150 seconds. Because the number of time steps is fixed, the error with respect to the spatial discretization e m is defined as the root-mean-square error (RMSE) as, where n is the total number of time steps;ξ mesh i, t denotes the global saturation degree of the i-th mesh scheme at t-th second. Figure 7a shows that the convergence rate decreases as the element number increases. The error dropping from the coarse to the regular mesh is much higher than the dropping from the regular to the fine mesh. Especially, the error reduction from the regular to the fine mesh is limited small. Thus, the fine mesh is "fine enough" for the given case. Figure 7b shows the convergence study of the time step size when the fine mesh is chosen. Since the mesh size is fixed in this study, the error is compared with the nodal saturation degree ξ of the time steps 1, 2 and 3 relating to the time step 4. With respect to the time discretization, the error e t is defined as, where NNO equals to the total node numbers of the fine mesh scheme;ξ time step n, i represents the i-th node of n-th time step size case. Similar to the convergence study of the  mesh size, the error drops significantly from 0.05 seconds to 0.01 seconds; however, the error decreases much slower when the time step size decreases from 0.01 seconds to 0.001 seconds. Because of this, the 0.001 seconds is chosen as the time step size to save computational effort while still keeping a good accuracy.
It is also interesting to track the flow in the deformable preform. Figure 8 plots the resin flow patterns during the infusion process. The red regions represent the fullsaturated parts of the preform; and the blue regions indicate the dry parts; the gradients between red and blue illustrate the process zone, where the flow front locates. The tiny white lines show the flow directions in each element.
At the beginning of the process (0 -50 seconds), the resin flow moves towards all directions simultaneously and forms a sector pattern. Once the flow reaches the hole (50 -140 seconds), the hole drags down the neighboring flow speed. As a result, the side flow gradually moves faster than the center flow. At 137 seconds, the outside flow has already caught up the inner flow around the hole. When the entire flow front has passed the hole, the flow front forms a nearly straight line as shown at 217 seconds. Figure 9 shows the pressure distribution. We can notice that there are local bands of low-pressure zones (blue bands) around flow fronts. The pressure drops to the bottom where the global minimum pressure locates, then just after the flow front, it rises again and forms a local pressure peak.
We also plot the preform profiles along different paths in Fig. 10, by using the finest mesh and 1 × 10 −3 seconds as time step size. Due to the vacuum, the preform is  compressed when the process just begins; after that, the resin moves into and inflates the preform up. From the front section, we can see that the inflation occurs in the vicinity of the inlet. This significant expansion happens immediately after the resin moves in. As the process going, the resin flow moves from the center to sides, which makes the swollen part expanding outwards gradually; and the depth of the swollen part reducing in the meanwhile. What is more, around the inlet, the inflation speed consistently slows down as the pressure approaches the atmospheric pressure.
The center, middle and side section mainly exhibit similar deformation profiles. The longitudinal profiles compress first then expand. The left side parts (x < 0) have larger expansion than the right side parts (x > 0). At 50, 100, and 200 seconds, there are local jumps around the flow front regions in the middle and side section. This behavior reflects low pressure bands in Fig. 9. After 200 seconds, the local jumps disappear and profiles turn smoother.

Conclusion
In this paper, we have presented a novel finite element based model for the LCM process of the thin-walled FRPCM components. The model contributes a shell element in the context of the theory of porous media to model both the preform deformation and the in-plane resin flow. In this fashion, the full 3-D fluid-structure interaction problem is reduced to a 2-D porous media problem, which significantly reduces the number of the degrees of freedom. The model solves for three primary variables: 1) the saturation degree ξ and 2) the fluid pressure p are solved from the mass balance relations that describe the homogenized Darcy flow transporting in the porous media; 3) the normal stretch λ is then obtained from the explicit formulation Eq. 53, which is derived from the linear momentum balance. The coupled Eqs. 55 and 56 are solved using the Streamline-Upwind/Petrov-Galerkin stabilized finite element method following the staggered approach, as proposed in [15]. The in-plane flow assumption has been justified through the comparison between the 2-D plane and shell models. The convergence studies for the doubly curved specimen shows that the method is converging with respect to the saturation degree for reasonable spatial and temporal discretizations. The model shows the ability to predict different mold filling settings, e.g., multiple resin inlets and different outlet positions. Besides, if the stretch field is prescribed, the same model can also be applied to the RTM process. Thus, this study proposes a simulation tool to optimize the wetout process in the deformable thin-walled fiber preform for large-scale FRPCM components, which can substitute for the costly physical trial and error approaches.