Numerical Simulation of Two-Phase Flow in Liquid Composite Moulding Using VOF-Based Implicit Time-Stepping Scheme

Abstract

The filling stage in injection/infusion moulding processes plays a key role in composite manufacturing that can be influenced by the inlet and vent ports. This will affect the production of void-free parts and the desirable process time. Flow control is usually required in experiments to optimise such a stage; however, numerical simulations can be alternatively used to predict manufacturing-induced deficiencies and potentially remove them in the actual experiments. This study uses ANSYS Fluent software to model flow-front advancement during the impregnation of woven fabrics. A developed technique is applied by creating tracking points (e.g., on-line monitor) in the direction of the flow to report/collect data for flow-front positions as a function of time. The study adopts the FVM-VOF-based two-phase flow model together with an implicit time-stepping scheme, i.e., a dual-time formulation solution method with a preconditioned pseudo-time derivative. Initially, three time-step sizes, 5 s (small), 25 s, and 50 s (large), are evaluated to examine their impact on numerical saturation lines at various fabric porosities, 40%, 50%, and 60%, for a two-dimensional (2D) rectangular mould, and predictions are then compared with the well-known analytical Darcy. This is followed by a three-dimensional (3D) curved mould for a fillet L-shaped structure, wherein the degree-of-curvature of fibre preforms is incorporated using a User-Defined Function (UDF) to tailor the impregnation process. The developed approach shows its validation (1–5.7%) with theoretical calculations and experimental data for 2D and 3D cases, respectively. The results also stress that a shorter computational time can be achieved with a large time-step size while maintaining the same level of accuracy.

1. Introduction

During liquid composite moulding (LCM) processes, a liquid resin is injected/infused at a constant pressure/flow rate to impregnate a fibrous reinforcement [1,2]. This serves to fully saturate the placed/laid fibre preforms on the mould. This means that the fibre preforms are dry (unsaturated) during the pre-filling stage (c.f., Figure 1). The dry fibrous reinforcement needs to be processed with a thermosetting resin in order to produce a composite part. This process stage is referred to as the mould-filling stage, in which the injected/infused resin flow starts to impregnate the fibre preforms. The quality of the processed fibre–reinforcement composites would depend on this crucial process stage. This is due to the fact that gelation or early curing, as well as the emergence of voids, could occur within the impregnation process. To avoid these issues, the knowledge and control of filling time and flow-front advancement is required. For this reason, it is important to optimise the mould-filling process in an effective way to obtain good-quality composite parts/components. This optimisation can be done with the traditional trial and error method; however, this method shows drawbacks from the perspective of time and cost [3]. Recently, there have been alternative proposals to optimise such a process, and one of the most efficient methods is to use the available sophisticated numerical tools. This can include isothermal/non-isothermal conditions, but this study will focus on isothermal mould-filling processing.

Flow-front advancement has been discussed by a considerable amount of research [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]. This includes experimental [4,5,6,7,8,9,10,11,12] and numerical studies [5,8,13,22], in which the latter are subject to validation of the adopted numerical tool with the available experimental data or theoretical models. Experimental works for monitoring/tracking the flow progression used various types of techniques; these would include cameras, ultrasonic sensors, optic sensors, dielectric sensors, and pressure transducers. On the other hand, the developed numerical models involve a wide range of approaches that are based on the control volume finite element method (CVFEM), finite volume method (FVM), finite element method (FEM), and finite difference method (FDM). These discretisation methods allow the solving of multi-phase or multi-physics flow problems, and this could apply to resin flow in LCM processes. Hoes et al. [4] conducted an experimental work that involved filling front tracking at a constant pressure injection. The authors placed electric sensors together with pressure transducers in the bottom flat mould plate, wherein these sensors were located on straight lines and predefined in the connected data collection computer. This method enabled the automatic observation of flow-front positions as a function of time. Binetruy et al. [5] investigated, experimentally and theoretically, the interaction between micro- and macro-flow through unidirectional woven fabrics, and their impacts on the flow-front progression during the impregnation phase. The authors used a flat Plexiglass (transparent) piece on the top of the mould to allow visual observations of the filling process. The study highlighted that, for unsaturated (dry) tows, a lag of the intra-tow flow progression would lead the macro-flows to cause transverse fluid penetration into the yarns. This study concluded that the micropores could affect flow-front behaviour and slow the filling process. A contribution by Luce et al. [6] was presented by performing mould-filling experiments for multi-ply fibre preforms that consisted of different fabric architectures, such as 3-D woven and random mat (RM). Two cameras were placed on the top and side-view of the mould to monitor and characterise the in-plane and through-thickness advancing flow behaviour. The study showed that the flow-front advancement through RM behaves uniformly with rapid progression. However, this was different for the 3-D woven ply, in which the flow encountered progress resistance caused by the transverse tows, besides the emergence (presence) of out-plane flows resulting from advancing macro-infiltration. The application of ultrasonics was employed by Schmachtenberg et al. [7] to obtain the on-line monitoring of the flow-front propagation in the RTM manufacturing process. The generated experimental data stressed that signals, a transmission time, and an amplitude could explain the arriving flow-front, and hence achieve the sensitive optimisation of the RTM filling process. The work by Pierpaolo et al. [8] embedded dielectric sensors to capture pressure data of the infused unidirectional flow in LCM, as well as monitoring the progressing flow front. The results were validated against numerical modelling, with satisfaction in terms of pressure profile and flow progression. Their study incorporated multiscale modelling, in which a microstructure of the used fabric (2-D woven), was scanned to obtain a meso-scale image, whereby the meso-permeability was identified and introduced in the macro-scale simulation to validate the approach with the experimental work. At present, a broad range of sophisticated software packages are available to model mould-filling processes in LCM, including RTM and VARTM (see

Table 1. Selection of numerical contributions that evaluated simulation tools for flow-front modelling.

References	Fabric Architecture	Injection Method	Flow Modelling	Computational Approach
Tan et al. [3]	Bi-axial	Unidirectional	Dual-scale	FE/CV-PoreFlow
Simacek et al. [13]	UD	Unidirectional	Dual-scale	FE/CV-LIMS
Oliveira et al. [14]	Fibre mats	Unidirectional	−	Darcy-based PAM-RTM
Grossing et al. [15]	UD/Triaxial NCF	Radial	Dual-scale	FVM-VOF OpenFoam
Sas et al. [16]	UD	Unidirectional	−	FEM-LSM COMSOL
Wei et al. [22]	Fibre mats	Unidirectional	−	FVM-VOF Moldex3D

). These would include, but are not limited to, LIMS [13], PAM-RTM [14,15], OpenFoam [15], and COMSOL [16]. Simacek et al. [13] implemented an algorithmic code based on the finite element/control volume (FE/CV) method using the simulation tool LIMS to model tow saturation during the mould-filling process. The numerical study considered the dual-scale flow, and emphasised that the developed model is able to capture the filling progression as well as the required time to fully impregnate and saturated the fibre preforms. They also added that this approach can be applied to arbitrary complex shapes. Voller et al. [17] followed CVFEM, but using Fortran to track flow-front positions as a function of time. The work analysed the filling front with structured/unstructured mesh grids for fibre–reinforcement mats. The results showed good agreement with the well-known analytical Darcy for various time-step sizes. Multiscale modelling of bi-axial fabrics was developed by Tan et al. [3] to predict the filling front during RTM processes. The numerical model adopted the FE/CV approach to simulate a coupled macro- and micro-flow problem. The PoreFlow program was used in the study, which is based on the Fortran modular package, and hence showed its capability of simulating resin impregnation for dual-scale porous media. This was validated with experimental data obtained by the same authors. Grössing et al. [15] performed a numerical mould-filling analysis via two different simulation tools, i.e., OpenFoam and PAM-RTM. The study evaluated both tools and highlighted that non-porous and porous zones can be both modelled by OpenFoam, whereas this was different with Darcy-based PAM-RTM, whereupon porous zones can only be modelled. The authors also argued that the issue with PAM-RTM could be resolved by assigning non-porous regions as race-tracking zones with %100 porosity value. The numerical study was based on an RTM experimental work that involved two types of fibre preforms, such as NCF and UD. Both numerical results agreed well with the flow-front experimental data.

The present work is motivated to contribute a simple, accurate technique using ANSYS Fluent to model resin flow advancement in the RTM/VARTM processes. The FVM-VOF-based multiphase flow approach is adopted to monitor the mould-filling process and to report flow-front positions as a function of time. This approach allows the prediction of the required time to fully saturate/impregnate the fibre preforms in simple and complex shapes. During the mould-filling simulation, the saturated, partially saturated, and unsaturated regions can be observed and located. The numerical simulations fit well with the analytical predictions and the experimental data for a rectilinear/channel flow injection.

2. Numerical Simulation Approach

2.1. Volume of Fluid (VOF)

The free-surface flows can be modelled by the prominent technique VOF, whereby the two immiscible fluids’ interface position is tracked and located simultaneously. This VOF model can be solved either using implicit or explicit time formulation (see Section 2.2) throughout structure/unstructured fixed meshes—a Eulerian-based approach. Such formulations (interpolation schemes) will allow discretisation of the volume fraction (tracking of the interface) by applying widely adopted options (in ANSYS Fluent)—geometric reconstruction (a piecewise-linear approach) and Modified HRIC (High-Resolution Interface Capturing)—for explicit and implicit, respectively [23,24,25]. This is to obtain face fluxes (fluid advection) for all filled and empty cells, and also near-interface (partially filled) cells [23,24,25]. Momentum and continuity equations are applied to each cell throughout the domain, in which the flow motion is solved by Navier–Stokes equations (N–S) based on the FVM discretisation (control-volume-based) approach. For RTM/VARTM filling simulation, the volume fraction in each control volume (computational grid cell) denotes either one of the phases (resin or air), or a mixture of both phases (partial saturation). This approach follows a transient, laminar, incompressible, multiphase flow problem with Newtonian behaviour and an isothermal filling condition. Thus, the continuity/transport equation can be written as follows:

$$\frac{\partial s_i}{\partial t}+\nabla ·\left(s_i\mathbf{u}\right)=0$$

(1)

where $\mathbf{u}$ is the volume-averaged velocity, t is the time, and $s_i$ denotes a phase volume fraction such that $s_f$ is the fluid (e.g., resin), while $s_a$ defines air (or empty) regions. This form assumes a constant fluid density and applies to fixed/stationary control volumes in this study. $s_f$ varies from 0 to 1 in each grid cell, and this would indicate full saturation for cells with a value of 1, partial saturation for cells with a value ranging from $0<s_f<1$, and unsaturation for cells with a value of 0. Since the resin flow impregnates a porous material, the momentum equation used by the simulation tool incorporates a source term to allow porous media modelling. Equation (2) demonstrates the momentum equation after applying the above-mentioned assumptions, while Equation (3) expresses the source/sink term:

$$\frac{\partial }{\partial t}\left(\rho \mathbf{u}\right)+\nabla \left(\rho \mathbf{uu}\right)=-\nabla p+\mu {\nabla }^2\mathbf{u}+S$$

(2)

$$S=-\frac{\mu }{\mathbf{K}}\mathbf{u}$$

(3)

Here, $\nabla p$ is the pressure gradient, $\rho$ is the fluid density, $\mu$ is the fluid viscosity, S is the sink term, and $\mathbf{K}$ is the permeability tensor. This filling modelling is performed on a rectangular mould (2-D RTM) in a macro-scale flow problem. The porous medium permeability is isotropic; therefore, the in-plane permeability can be denoted as $K_{xx}=K_{yy}=K$. The permeability value was computed in a previous work by the current authors at a dual-scale flow (inter- and intra-tow flow) through a plain-weave fabric unit cell at various aggregate porosity values: 40%, 50%, and 60%. Due to the fact that the unit cell can be taken as a representative volume element (RVE) of the woven fabric ply, the obtained in-plane dual-scale permeability is inputted in the source term to enable true viscous resistance of the woven model at a macroscopic level. The developed approach used the implicit scheme (see Section 2.2) provided by ANSYS Fluent, owing to the fact this method is preferred for slow flow movement problems, i.e., creeping flows, quasi-static motion, and static/dynamic motion [26]. In this method, the time-step size is independent of the results, and can be large to obtain a shorter computational time [17,26]. This work developed a flow-front tracking technique by creating points on a straight line (flow direction) at certain locations. Each point will report resin flow volume fraction data as a function of time; this would be similar to real-time monitoring RTM/VARTM experiments wherein sensors are embedded/placed (in fabrics or on mould).

As part of this study, three-dimensional flow simulations of complex shapes are included to examine the flow-front prediction model for a real filling process of a composite component. An experimental work by Geng et al. [27] was selected, wherein a curved mould was used with a VARTM process type. Hence, the present work added a set of equations (see Equation (Section 2.1)) that accounts for the effect of curvature on resin flow through the curved regions during the resin impregnation of fibre preforms. This is done by using the User-Defined Functions (UDF) to interpret/compile these equations with a C-language-based code.

$$\begin{array}{c}\hfill D_c=2{sin}^{-1}\left(\frac{C}{2R_c}\right),\theta =\frac{D_c}{2}\end{array}$$

(4a)

$$\begin{array}{c}\hfill K_H\left(\theta \right)=K_{H}cos\left(\theta \right)\end{array}$$

(4b)

$$\begin{array}{c}\hfill K_V\left(\theta \right)=K_{V}sin\left(\theta \right)\end{array}$$

(4c)

The degree of curvature and the deflection angle are $D_c$ and $\theta$, respectively. Here, $R_c$ is the radius of curvature, and C is the cord length. Since part of the filling process is vertical, the permeability is set as a function of curvature for both vertical and horizontal resin flow progressions.

2.2. Time-Stepping Scheme

A temporal discretisation in transient simulations is defined by space and time. This means that an integration is required for each term in the differential equations with every time-step $\Delta t$. To apply such a time discretisation, there are two prominent schemes, i.e., implicit and explicit. These methods are to be carefully selected, since a diverse range of restrictions arise in each scheme in certain circumstances. The implicit method is applicable to slow motions (e.g., creeping flows), and it can handle long-duration analyses that vary from seconds to days [26]. On the other hand, the explicit case is suitable for rapid motions (e.g., drop tests, shocks, etc.), and it is restricted by the Courant–Friedrich–Lewy (CFL) condition, in which each grid cell uses the same time-step that depends on the courant number, cell volume, and the sum of outcoming fluxes [26]. Therefore, this work finds the implicit time-stepping a suitable technique for such a flow problem, and an evaluation was conducted to examine the impact of time-stepping size on the flow-front modelling. It is noteworthy that implicit time-stepping is a dual-time formulation that adopts a preconditioned pseudo-time (inner iterations) derivative at each physical time-step (Euler backward) to provide an accurate transient solution [26,28]. Figure 2 describes the time discretisation methods and how the time-dependent equations are computed differently throughout the domain. This shows that a function of a quantity/variable $f\left(\phi \right)$ is evaluated at a future time level $(n+1)$ for the implicit case, while $f\left(\phi \right)$ is computed at the current time level $\left(n\right)$ for the explicit time-stepping method. The general expression used in the current numerical simulation is demonstrated below.

$$\begin{array}{c}\hfill \frac{{\phi }^{n+1}-{\phi }^n}{\Delta t}=f\left(\phi \right)\end{array}$$

(5a)

$$\begin{array}{c}\hfill \left\{\begin{array}{c}n+1=t+\Delta t\\ n=t\\ n-1=t-\Delta t\end{array}\right.\end{array}$$

(5b)

2.3. Darcy’s Law for a Transient Flow

The general form of Darcy’s law (Equation (6)) for an isothermal filling flow problem can be integrated and simplified to suit transient flow conditions [29]. At a constant pressure injection for one-directional flow movement (e.g., x-coordinate), Equation (6) can be rewritten as shown in Equation (7). By further simplification and integration of Equation (7), the position of the flow-front can be obtained by Equation (8).

$$\mathbf{u}=-\frac{\mathbf{K}}{\mu }\nabla p$$

(6)

$$\frac{u_x}{{\varphi }_{\mathrm{o}}}=\frac{dx}{dt}=-\frac{K}{\mu {\varphi }_{\mathrm{o}}}\frac{\Delta p}{\Delta x}$$

(7)

$$x_f=\sqrt{\frac{2K{p}_o}{\mu {\varphi }_{\mathrm{o}}}{t}_f}$$

(8)

where $x_f$ is the flow-front position, $t_f$ is the flow-front time, ${\varphi }_o$ is the porosity of the medium, and since the pressure at the flow-front can be assumed as zero, the pressure difference $\Delta p$ becomes equal to the injection pressure $\left(p_o\right)$ [29]. Thus, the proposed approach in the present study can be subjected to a comparison analysis with the above-mentioned theoretical solution for validation. On this premise, the developed numerical model will be also feasible for complex structural shapes, and this would contribute to the knowledge and decision of LCM optimisation. This is in terms of the optimal cycle time and location of injection ports/vents, thereby achieving void-free and good-quality composite parts. Figure 3 illustrates the geometry and boundary conditions used in the numerical analysis for the 2-D rectangular mould filling involving the flow-front tracking points.

3. Results and Discussion

3.1. Two-Dimensional Rectangular Mould for Regular Shapes

The numerical simulations were performed at various fibre preform aggregate porosities $\left({\varphi }_o\right)$ $40\%$, $50\%$, and $60\%$, with in-plane dual-scale permeabilities $\left(K\left[{\mathrm{m}}^2\right]\right)$ $4.89\times {10}^{-11}$, $9.45\times {10}^{-11}$, and $2.09\times {10}^{-10}$, respectively. Filling front analysis was computed at a constant pressure injection ($p_o=10\left[\mathrm{kPa}\right]$), and with fluid flow properties $\rho =1300[\mathrm{kg}⁄{\mathrm{m}}^3]$, and $\mu =0.15[\mathrm{Pa}·\mathrm{s}]$. Three time step-sizes, 5 s, 25 s, and 50 s, were selected based on the generated grids/meshes for the transient multiphase flow through a 2-D macro-scale geometry. These different step sizes yielded significant factors in terms of the actual flow-front position or the so-called saturation lines and the computational time. It appeared that the partially filled control volumes (i.e., partially saturated region) could indicate the actual flow-front for all three time-step sizes. This can also be seen in similar numerical studies, such as [3], in which the flow-front was observed within the partially filled computational grid cells. Accordingly, the results stressed the fact that the flow-front progression with 40–50% resin volume fraction agreed well with the analytical solutions for all cases. Less computational time was observed at large time-steps such as 50 s; however, smaller step sizes are preferred (e.g., determined by grid cell size divided by fluid velocity value) [26]. By applying this, sharp flow-front behaviour could be seen in which the partial saturation zone was reduced (c.f., Figure 4). Consequently, filling front positions would be more manifested or perceivable, and an accurate on-line flow control could be attained. It should be noted that different implicit-scheme-based time-steps would offer the same required time for the mould-filling process; nonetheless, computational time processing would vary. For instance, to accomplish a fully saturated fibre preform with the geometry dimensions and boundary conditions illustrated in Figure 5 for the 50% porosity value case, $1600 \mathrm{s}$ is needed. This was achieved by all time-stepping sizes, $5 \mathrm{s}$, $25 \mathrm{s}$, and $50 \mathrm{s}$, but with different numbers of time-steps: 320, 64, and 32 respectively. This means that a shorter computational time can be obtained with the maximum time-stepping size, while maintaining the same accuracy of the output results, and this can be explained by Figure 5 together with

Table 2. Computational processing time with different time-stepping sizes for each medium porosity.

Medium Porosity [%]	Time-Stepping Size [s]	CPU [s]	Real Time [s]
40	5	19.75	481
	25	5.44	89
	50	2.97	57
50	5	14.67	319
	25	3.76	67
	50	2.2	40
60	5	8.14	148
	25	2.26	42
	50	1.55	26

. Therefore, when supposing that the user is interested in a shorter computational time to optimise the mould-filling process, it would be recommended to follow the maximum time-stepping in accordance with the mesh/grid element size. It is worth mentioning that an Intel Core i7-1165G7 Processor (Central Processing Unit CPU) was used to run the present numerical simulations. Figure 5 presents the computed results for various woven fabric aggregate porosities, 40%, 50%, and 60%, with different time-step sizes, $5 \mathrm{s}$, $25 \mathrm{s}$, and $50 \mathrm{s}$, for each porous medium case. The results showed perfect alignment with ones calculated by transient Darcy for all flow-front simulations. Therefore, the numerical approach proposed in this study shows its reliability and can be confidently used to optimise RTM/VARTM filling processes.

3.2. Three-Dimensional Curved Mould for Complex Shapes

For a complex shape, a study by Geng et al. [27] was selected to examine the current approach against the populated experimental data. Figure 6 shows a schematic diagram of the experimental setup that was used by [27], in which VARTM was adopted. Geng et al. [27] managed to perform experimental works for curved composite shapes using VARTM, and this was done with different bend angles, i.e., 180, 120, 90, and 60 degrees. Therefore, this paper assesses the developed numerical flow-front prediction for a three-dimensional complex shape during the mould-filling process. This is applied to curved non-crimp fabric (NCF) plates with 90 degree bending, and simulations are conducted on single and multiple layers based on the VARTM experiment by Geng et al. [27]. The permeabilities of fibre preforms are set to be $2.98\times {10}^{-11}$ $\left[{\mathrm{m}}^2\right]$ and $4.6\times {10}^{-11}$ $\left[{\mathrm{m}}^2\right]$, with fibre volume fractions $\left(V_f\right)$ $22\%$ and $40\%$ for 1-layer and 6-layer, respectively [27]. This is along with the flow properties 968 $[\mathrm{kg}⁄{\mathrm{m}}^3]$ and $0.35$ $[\mathrm{Pa}·\mathrm{s}]$ for density and viscosity, and a volume-average velocity range $\left(\mathbf{u}\right)$ 0.00025–0.00035 $[\mathrm{m}⁄\mathrm{s}]$ [27]. With such a curved composite part, the curvature region is said to be impacting the resin impregnation, as was thoroughly discussed by Alotaibi et al. [30], and hence it is considered in the present work. In such a case, the resin flow will be affected by the degree of curvature as long as it progresses within the curved zone—see Figure 7. This would show 90 degrees of curvature $\left(D_c\right)$ with a deflection angle impact range $(0^{\circ }$≤ $\theta$ ≤${45}^{\circ })$. Thereby, a set of equations is required to be incorporated with the flow equations—see Equation (7)—and this is done using the User-Defined Function (UDF) in ANSYS Fluent. The results give a good impression of the capability of the current numerical model, in which most of the tracking flow-front points (c.f., Figure 8) fit well, showing $\le 5.7\%$ discrepancy with the experimental data. This is in addition to the mould-filling time, wherein times of 760 s and 475 s were achieved for 1-ply and 6-ply, respectively, in comparison to 750 s and 470 s experimental fill time. Figure 9 portrays the resin flow progression contours/outlines for 1-ply and 6-ply NCFs, and with a more permeable fibre perform (e.g., 6-ply case), this would hasten the flow advancement. It should be noted that a hexahedron mesh topology is preferred when applying the relevant degree-of-curvature equations through the UDF, since it provides a sequential quadrilateral cell type, and this makes it compatible with the degree-of-curvature concept. The present numerical methodology proves its capability to monitor and predict resin flow advancement in RTM/VARTM processes that include regular and complex shapes, and it also offers a simple and accurate technique that permits efficient computational modelling.

4. Conclusions

A numerical technique was proposed using the data-track-point feature of ANSYS Fluent to predict the flow-front position as a function of time throughout the computational domain, as a means to simultaneously monitor filling progression during LCM processes (e.g., RTM/VARTM). The results stressed that the grid-independent solution (implicit VOF) of structured or unstructured control volumes does not impact the output of flow parameters, while allowing less computational time at larger time-step sizes. The simulation approach followed the FVM-VOF-based multiphase flow problem together with implicit time discretisation. The numerical model has been validated with Darcy’s law and an experimental work [27] for a transient flow at various aggregate porosities of dual-scale fabrics. It is noteworthy that the present work assumes no chemical conversion of the fluid flow (i.e., a constant viscosity) during the impregnation process. A future work will include a cure-temperature-time-dependent viscosity to characterise their impacts on the impregnation/saturation of complex (woven) fabric structures with heterogeneous permeability.