Experimental validation of a new adaptable LCM mold filling software

Abstract

Resin Transfer Molding (RTM) is a manufacturing process for fiber reinforced polymer composites where dry fibers are placed inside a mold and resin is injected under pressure. During mold design, filling simulations can study different manufacturing concepts (i.e. placement of injection gates and vents) to guarantee complete filling of the part and avoid air entrapment where flow fronts converge. In this work, a novel software tool LCMsim, which was implemented by the authors, is benchmarked against other tools and real-world flow experiments. Its development was driven by two ideas: Easy-of-use for the mold engineer and maximum flexibility for the researcher. Two experiments were used for validation. In the first, zones with different preform properties were present and in the second, race-tracking was enforced. Flow fronts from LCMsim and experiment agree with 7% error and simulated flow fronts from LCMsim and the commercially available software PAM-RTM agree with 3% error.

GRAPHICAL ABSTRACT

Keywords:

• Resin transfer molding (RTM)
• liquid composite molding (LCM)
• filling simulation
• computational fluid dynamics (CFD)
• shell mesh
• finite area method
• flow front tracking
• resin flow

1. Introduction

1.1. LCM process

In all liquid composite molding (LCM) technologies dry fiber preforms are placed into a mold which is subsequently closed or sealed with a flexible bag by drawing a vacuum. Resin then enters into the mold until the preform is fully wetted with resin. As soon as the resin is cured the mold is opened and the part is extracted. In Resin Transfer Molding (RTM) resin is injected into a matched mold under pressure which gives excellent surface finish on both sides and high fiber volume fractions (55-65%) can be obtained. The mold walls are considered rigid and the preform is stationary during injection. In contrast to RTM, in VARI vacuum pulls resin through the preform on a single-sided mold with vacuum bagging. For a complete process description see e.g. Strong [1] or Niu [2].

1.2. LCM filling simulations

During mold design, filling simulations can study different manufacturing concepts (i.e. placement of injection gates and vents) to guarantee complete filling of the part and avoid air entrapment where flow fronts converge.

Required input for filling simulations are the mold cavity geometry, preform properties (e.g. permeability, porosity and thickness values which describe the fibrous preform regions), process parameters (e.g. initial cavity pressure and injection pressure) and fluid properties (e.g. resin viscosity).

During pre-processing the physical domain (i.e. the mold cavity) is transferred to a numerical domain (i.e. the mesh) and preform properties and initial conditions for the flow variables are assigned to the cells of the mesh. The solver calculates how the values of the flow variables in the cells change in the course of the filling. Different physical models of the filling can be implemented in the solver. The physical models have to describe all relevant issues of the filling process. During post-processing the resin flow front at different time instances is visualized with contour plots.

Before a simulation is started, parameters which describe the flow through porous media (permeability and porosity) must be determined. In-plane permeability is characterized experimentally from linear or radial flow experiments where the flow front location is tracked over time and the entries of the in-plane permeability tensor are computed according to specifically developed mathematical algorithms, see May et al. [3]. Experiments for different porosity levels are performed. Porosity is the fraction of the volume of voids over the total volume. For fibrous preforms inside a mold, the volume averaged porosity can be calculated from cavity thickness t, number of plies n, mass density ρ_f and areal weight A_f of the preform material as $\epsilon =1-(n{A}_f)/({\rho }_{f}t).$ The second term is the ratio of the preform volume over the total volume, see for example Neitzel et al. [4].

1.3. Requirements for LCM software

General purpose CFD software packages such as ANSYS® Fluent® or OpenFOAM as well as specialized software packages such as PAM-RTM®, RTM-Worx®, LIMS® or myRTM can be used for filling simulations.

The models implemented in the different software packages mentioned above perform well for most scientific and engineering applications, but there are applications where no satisfactory solution is found, e.g. for Compression-RTM where an equation of state to model fluid compressibility must be considered. See e.g. Bickerton and Abdullah [5] for a process description.

Since the required physical model to describe all relevant issues of the filling process often is not known in advance, a functional requirement for a new filling software is maximum flexibility for the user. Maximum flexibility is only possible if the code is open-source and code modification is possible for people with basic CFD and programming knowledge, similar to creating new solvers in the open-source CFD package OpenFOAM. Typical code modifications are the implementation of different boundary conditions (control pressure level or control opening and closing of ports) or adding additional equations (e.g. temperature or degree-of-cure equations since the resin viscosity is known to depend on temperature and degree of cure) or extracting the solver section (i.e. the loop for the time evolution) and integrate in separate code. The latter was already done for a control study with three inlet ports to keep the flow front straight in a rectilinear flow with reinforcement patches, see Stieber et al. [6]. If different models can be implemented easily, as a side result, comparisons between different models are possible also for complex geometries and not only for linear and radial test cases where analytical formulas are available.

In addition, mold engineers need an easy-to-use software to determine the optimal location of gates and vents. This is possible, if a shell mesh of the cavity with specified sets for different preform regions is available and the software can be started from a GUI with only preform-, process- and ports-relevant input data.

Commercial software packages (with the advantage of a user-friendly GUI, an extensive documentation, a large number of industrial case studies,…and with the obvious limitation in customizability) cannot fulfill these functional requirements. Non-commercial software packages from the list above are the open source CFD software package OpenFOAM and myRTM, a free industry standard for filling simulations. OpenFOAM was already used for implementing and testing new filling models, see for example Seuffert et al. [7] or Sebastian et al. [8], but OpenFOAM comes without GUI and requires both CFD and programming knowledge to perform appropriate LCM filling simulations. myRTM can only simulate RTM filling. The predicted filling pattern is correct but the predicted filling time differs considerably from that of the actual component, see Barandun et al. [9] and myRTM [10]. In contrast to RTMsim and LCMsim, myRTM includes a mesher for surface geometries and preform properties can be assigned graphically. The main differences between LCMsim and the pre-existing RTMsim is the capability to perform well for all relevant viscosity and permeability ranges, the capability to simulate VARI filling in addition to RTM and the customizability for other process variants. Table 1 compares the new LCMsim module with the pre-existing RTMsim module.

Table 1. Comparison of the new LCMsim module with the pre-existing RTMsim module.

	RTMsim module	LCMsim module
Advantages	• GUI or text input file for assignment of preform and process parameters • Open source • Validation and verification cases for RTM available, see radial flow test case below and Obertscheider and Fauster [11]	GUI or text input file for assignment of preform and process parameters Open source Prepared for code extensions/modifications for people with basic CFD or programming knowledge Flow physics model allows for RTM and VARI simulations and also other process variants
Limitations	• Shell mesh required, no mesher for surfaces included • Flow front propagation too slow for high viscosity and small permeability, see linear flow test case below • Only RTM simulations	Shell mesh required, no mesher for surfaces included Not yet run-time optimized
Radial flow test case from Isoldi et al. [12] to show similar filling pattern for typical RTM filling test case	Preform size $0.6\times 0.6\times 0.003$ m^3, central injection gate (grey) with 0.02 m diameter, 35,000 Pa difference between injection and initial cavity pressure, dynamic viscosity 0.06 Pas, porosity 0.7, isotropic permeability $3·{10}^{-10}$ m^2. Same mesh for all simulations. Contour plot of the filling fraction after 200 s (yellow/blue is filled/unfilled):
Linear flow test case to show differences for increased viscosity and decreased permeability	Preform size $0.2\times 0.06\times 0.003$ m^3, line injection at the left (grey), 35,000 Pa difference between injection and initial cavity pressure, dynamic viscosity 0.1 Pas, porosity 0.7, isotropic permeability $3·{10}^{-11}$ m^2. Same mesh for all simulations. Flow front position at 0.17 m calculated with analytical formula, see for example Isoldi et al. [12]. Contour plot of the filling fraction after 1000s:

Notes: The first two lines list advantages and limitations. The third line shows the results of a typical radial RTM flow test case. The fourth line shows the results of a linear RTM flow test case with increased viscosity and decreased permeability. RTMsim gives reliable results for most RTM test cases, particularly if only the filling pattern is considered. LCMsim in contrast to RTMsim performs well for all relevant viscosity and permeability ranges (see the difference in the linear flow) and can simulate RTM, VARI and other process variants.

1.4. Description of new LCM software

The existing Julia module RTMsim described in Obertscheider and Fauster [11] for iso-thermal RTM filling simulations was extended with a new fluid model and with new functionalities. The extended Julia module LCMsim can be downloaded from the GitHub repository https://github.com/obertscheiderfhwn/LCMsim. Installation instructions and instructions on how to run a simulation are given in a README file on the repository. Figure 1 shows the GUI with explanations and typical values.

Figure 1. LCMsim GUI filled with parameters for the radial flow test case from Table 1. All parameters in SI units. The GUI is divided into four major sections as highlighted, with the most important fields of input data mentioned in the describing text.

The work at hand describes the flow physics model which is implemented in LCMsim in detail and presents a validation of LCMsim for RTM. Real-world RTM experiments under real process conditions, i.e. permeameter experiments with zones with different in-plane permeability and porosity levels, are the basis for the RTM assessment.

2. Materials and methods

2.1. Physical model

RTM is a manufacturing process suitable for producing thin-walled fiber reinforced polymer composites (FRPC) parts. The thickness of the part is much smaller than the overall dimensions as depicted in Figure 2. Dry reinforcement fibers are placed inside a mold and resin is injected under pressure into the fibrous preform. On a macroscopic level, the physical quantities porosity and permeability are used to describe resin flow through the fibrous media.

Figure 2. A schematic of a thin-shell RTM geometry (in-plane dimension much greater than the out-of-plane dimensions, $d_{\text{ip}}\gg d_{\text{oop}}$) for manufacturing a flat plate with linear flow and line injection gate (a) and with radial flow and single point injection (b).

For usual textiles the out-of-plane permeability $k_{\text{oop}}$ is one or two orders of magnitude smaller than the in-plane permeability $k_{\text{ip}},$ see for example Simacek and Advani [13] or Sirtautas et al. [14]. For thin-walled parts considered in this study the thickness $d_{\text{oop}}$ is two orders of magnitude smaller than the characteristic in-plane dimension $d_{\text{ip}}.$ This gives an aspect ratio (1) $$AR=\frac{d_{\text{oop}}}{d_{\text{ip}}}\sqrt{\frac{k_{\text{ip}}}{k_{\text{oop}}}}$$(1) between 0.01 and 0.1. Only for higher aspect ratio (for example thick laminates or a flow distribution medium on top of the preform) three-dimensional modeling of the flow is desirable [13]. Consequently, transverse flow is neglected and the flow is modeled as two-dimensional, considering in-plane flow contributions only.

In addition one assumes slip boundary conditions at the walls of the cavity since the flow is described on a macroscopic level. In the absence of a porous preform a flow profile with zero velocity at the top and bottom cavity walls and maximum velocity in the middle would develop. Since the cavity is filled with the preform, flow takes place in small flow channels with potential cross-flow but is described by a superficial velocity and the assumption of a plain flow front is valid. Figure 3 illustrates this difference. The physical velocity through the small channels is related to the superficial velocity via the porosity (fraction of the void volume over the total volume). The superficial velocity is smaller than the physical such that the same volumetric flow rate moves through a cross sectional area, see e.g. Tucker [15] or Nield and Bejan [16].

Figure 3. Flow velocity profile through an empty cavity (a) and assumed volume-averaged (superficial) flow velocity profile through a porous cavity (b).

In contrast to the original model from Obertscheider and Fauster [11], a more realistic model was implemented in the novel LCMsim software module. The implemented model results in an algorithm which (a) is explicit in time and (b) avoids pressure-velocity coupling iteration during one time step.

Initially the porous cavity is either evacuated or filled with air (air density according to cavity pressure) and pressurized resin is injected into the fibrous preform.

Fluid flow is described by (i) a compressible continuity equation for the mixture density, (ii) an equation-of-state which describes the pressure build-up with increasing mixture density and (iii) a momentum equation with dynamic viscosity of the resin in the Darcy term. The in-plane velocity vector $\mathit{u}=(u,v)$ is described in the system of the orthotropic in-plane permeability. For tracking the flow front a volume-of-fluid equation is solved.

The physical quantities in the governing equations are functions of space x and time t where $\mathit{x}\in \Omega \subset {\mathbb{R}}^3$ and Ω is the fluid body also called the fluid cavity. The governing equations¹ are (2) $$\frac{\partial \epsilon \rho }{\partial t}+\nabla ·\left(\rho \mathit{u}\right)=0$$(2) (3) $$\frac{\partial \rho \mathit{u}}{\partial t}+\nabla ·\left(\rho \mathit{uu}\right)=-\nabla p+S_u$$(3) where

• $\rho =c{\rho }_{\mathrm{r}\mathit{esin}}+(1-c){\rho }_{\mathrm{a}ir}$ is the mass density of the mixture,

• ε is the porosity of the preform,

• ${\rho }_{\mathrm{r}\mathit{esin}}$ is resin mass density at injection pressure,

• ${\rho }_{\mathrm{a}ir}$ is the air mass density at initial cavity pressure,

• the binary fraction function c is equal to c = 0 (i.e. empty) if $\rho <\overline{\rho }$ and c = 1 (i.e. completely filled) if $\rho \ge \overline{\rho }$ and the threshold value is $\overline{\rho }=({\rho }_{\mathrm{r}\mathit{esin}}+{\rho }_{\mathrm{a}ir})/2,$

• u is the superficial velocity in a local coordinate system which is related to the physical velocity via the porosity ε given by the stationary fibrous media,

• $\mathit{uu}$ is the dyadic product,

• p is the fluid pressure,

• S_u is a source term which contains the pressure loss from flow through a porous medium according to Darcy’s law, i.e. $S_u=-\mu {\mathsf{K}}^{-1}\mathit{u}$ with dynamic viscosity μ and permeability tensor $\mathsf{K},$

• the contribution $\rho \mathbf{g}$ of the gravitational acceleration in the momentum equation is neglected, since in the currently considered applications it is assumed that gravity pressure effects are negligible compared to the pressure differences resulting from injection and cavity pressure,

• the contributions of the viscous stress tensor in the momentum equation are neglected, since in many porous media flow problems the viscous stress is much smaller than the solid-fluid drag, see Advani [17],

• the equation of state (see Figure 4) models the pressure rise as the mixture mass density increases. It is assumed that the pressure build up is slow at small mixture density (i.e. zero slope $\partial p/\partial \rho$). Alternatively, the real compressibility of the resin can be used for the calculation of the slope at injection pressure. A more realistic, reduced resin compressibility will result in much smaller time steps.

Figure 4. The equation of state models the pressure build-up for the air-resin mixture. A quadratic polynom fit through the two mass density/pressure pairs for air at initial cavity pressure and resin at injection pressure is used. The shown curve is for air density 1.225 kg/m³ at 100,000 Pa and resin density 960 kg/m³ at 191,000 Pa which is used in the validation cases.

2.2. Implementation

In order to numerically solve the flow model described by Equation (2)(2) $$\frac{\partial \epsilon \rho }{\partial t}+\nabla ·\left(\rho \mathit{u}\right)=0$$(2) and (3)(3) $$\frac{\partial \rho \mathit{u}}{\partial t}+\nabla ·\left(\rho \mathit{uu}\right)=-\nabla p+S_u$$(3) , respectively, with the equation of state from Figure 4, the computational domain (time and space) is discretized. Discretization of the governing equations on an unstructured mesh follows the ideas presented in Versteeg and Malalasekera [18].

The temporal domain (i.e. the simulation time) is split into a finite number of time steps where values for the physical quantities on the spatial domain are calculated in a time-marching manner. An explicit method with adaptive time stepping is used.

The spatial domain is defined on the cavity’s mid-surface and described by a shell mesh in a format similar to the NASTRAN bulk data format, see NASTRAN [19]. The mid-surface model can be curved and cells can have edges where more than two cells are connected to each other such as for handling T-junctions. A general mesh is illustrated in Figure 5(a). The injection gate is part of the computational domain and injection pressure and resin fraction is assigned to the cells belonging to this region. Vents can be specified optionally. If no vents are specified (similar to software myRTM) the simulated filling is not influenced as long as the flow front does not reach the fictitious vent position. Initial and boundary conditions are assigned. In all cells except in those belonging to the injection gates and vents the initial velocity is zero, the resin fraction is zero and the pressure is equal to the initial cavity pressure which is an input parameter. In all cells belonging to the injection gates the fluid pressure is equal to injection pressure and resin fraction is equal to one. In all cells belonging to the optional vents the fluid pressure is equal to initial cavity pressure and resin fraction is zero. The values in the cells for the injection gates and vents are not changed during a simulation. If injection pressure changes during the filling process or injection gates are added or removed, the simulation can be continued from the previous filling state with changed boundary conditions.

Figure 5. (a) Annulus filler mesh with T-sections and curved regions. (b) and (c) Locally flat mesh.

For discretizing the equations on the shell mesh of the part’s mid-surface it is assumed that the geometry is locally flat. The neighboring cells are rotated about the common edges to lie in the plane of the considered cell (see Figure 5(b,c)) and the velocity vectors of the neighboring cells are transformed into the same coordinate system.

The solver section (which is the relevant section if a different physical model is implemented) consists of a for-loop over all interior and wall cells inside a while-loop for the time evolution. For every interior and wall cell, the following steps are performed:

1. The pressure and filling fraction gradients are evaluated.

2. The time evolution for new mass density, x- and y-velocities with the discretized continuum, x- and y-momentum equations in the cell coordinate systems is performed.

3. The pressure is evaluated according to the equation of state.

4. The filling fraction is evaluated.

Boundary conditions need not be updated since in the inlet and outlet cells pressure and mass density remain unchanged and the velocity values at the cell boundaries are evaluated in the numerical flux evaluation.

The discretized versions of the governing Equations (2)(2) $$\frac{\partial \epsilon \rho }{\partial t}+\nabla ·\left(\rho \mathit{u}\right)=0$$(2) and (3)(3) $$\frac{\partial \rho \mathit{u}}{\partial t}+\nabla ·\left(\rho \mathit{uu}\right)=-\nabla p+S_u$$(3) , respectively, give the numerical scheme (4) $${\rho }_i^{n+1}=({\epsilon }_i{\rho }_i^n-\frac{\Delta t}{V_i}\sum _j{\mathit{n}}_{i,j}·\underset{{(\rho \mathit{u})}_{i,j}^n}{\underbrace{\frac{1}{2}\left({\rho }_i^n+{\rho }_j^n\right) \frac{1}{2}\left({\mathit{u}}_i^n+{\mathit{u}}_i^n\right)}}{A}_{i,j})/{\epsilon }_i$$(4) (5) $$u_i^{n+1}=\left({\rho }_i^n{u}_i^n-\frac{\Delta t}{V_i}\sum _j{\mathit{n}}_{i,j}·{(\rho \mathit{u})}_{i,j}^n{u}_{i,j}^n{A}_{i,j}-{\left(\frac{\partial p}{\partial x}\right)}_i\Delta t\right)/\left({\rho }_i^{n+1}+({\mu }_i/k_1)\Delta t\right)$$(5) (6) $$v_i^{n+1}=\left({\rho }_i^n{v}_i^n-\frac{\Delta t}{V_i}\sum _j{\mathit{n}}_{i,j}·{(\rho \mathit{u})}_{i,j}^n{v}_{i,j}^n{A}_{i,j}-{\left(\frac{\partial p}{\partial y}\right)}_i\Delta t\right)/\left({\rho }_i^{n+1}+({\mu }_i/k_2)\Delta t\right)$$(6) for cell i where $j\in S_i$ represents the index of a cell neighbor and the combined index i, j represents the index of the cell interface between cell i and cell j. The Darcy source term is treated implicitly. Otherwise a much smaller time step than given by the CFL condition introduced in Courant et al. (20] would be required. The physical quantities $u_{i,j}^n$ and $v_{i,j}^n$ at the cell interface are evaluated using first-order upwinding. For flow from cell i to cell j, first order upwinding is e.g. $u_{i,j}^n=u_i^n.$ The j-sum runs over all neighboring cells. When the neighboring cell is a pressure inlet cell, the velocity in the factor without upwinding in the convective term is evaluated according to Darcy’s law.

Code snippets shown in the appendix illustrate how the different flow physics models are implemented and how an additional equation can be added.

3. Results and discussion

3.1. Running a simulation

The simulation is started with a set of process parameters and a set of material properties, see the README file on the GitHub repository https://github.com/obertscheiderfhwn/LCMsim, either from the command line or with the GUI. Figure 1 shows the GUI with the relevant input parameters. Contour plots for the filling fraction and the pressure at different time instances are generated during post-processing.

3.2. Simulation overview

The RTMsim module was developed for fast filling simulations and implements a simple fluid flow model, whereas LCMsim implements a more complex flow physics model which can be used for the simulation of all LCM processes. The RTMsim module was verified with three RTM test cases as reported in Obertscheider and Fauster [11]. The new model implementation provided in LCMsim was initially tested against these test cases. Table 1 shows the results of a radial and a linear flow experiment run with RTMsim and LCMsim. The radial flow is a typical RTM filling test case from Isoldi et al. [12] and gives well comparable contour plots of the filling faction after 200 s. Significant differences between LCMsim and the RTMsim occur in rectilinear flow with high viscosity and small permeability as shown in the linear flow experiment. This is a result of modeling the resin as compressible in RTMsim, whereas the resin in LCMsim is only slightly compressible. If the resin is compressible, the pressure build-up in the cells requires more material and this slows down the flow. The slower the flow, the stronger is this effect. Slower flow is a result of higher viscosity and lower permeability.

Table 2 shows the validity of the used mesh (mesh 3) for the first validation case. The flow front position after 120 s filling time shows numerical convergence. The flow front position is comparable for meshes 2, 3 and 4. This is quantified by the calculation of the filled area. For the second validation case a finer mesh away from the pressure injection gate is used because the flow around a patch is investigated there.

Table 2. Results of a mesh refinement study for validation case 1 in terms of well comparable contour plots of the filling fraction after 120 s (yellow/blue is filled/unfilled).

	Mesh 1	Mesh 2	Mesh 3	Mesh 4
Mesh density	Coarse mesh	Finer at inlet, coarser outside	Medium mesh	Fine mesh
Filling pattern (whole domain)
Filling pattern (detail)
Number of cells	2516	2668	5150	12,296
CPU time	3.6%	8.2%	100%	323%
Filled area	0.282	0.401	0.405	0.411

Notes: The pictures in the first line show the whole domain, the pictures in the second line show the mesh around the pressure injection gate. The cell size is reduced by a factor 1.5 between meshes 2 and 3 and between meshes 3 and 4. Furthermore, the number of cells, the computational time (relative to the reference of mesh 3) and the filled area as fraction of filled by total domain area are given.

The following general conclusions for simulations with LCMsim can be drawn from this mesh refinement study. If the mesh at the injection gates is too coarse (see detail view of mesh 1) the flow front propagation speed is influenced. In general, the mesh size must be finer around the pressure injection gates and can increase with the distance from the gates (see mesh 2). Depending on the purpose of a filling simulation, a coarser or a finer mesh must be chosen. If the actual filling time is important a fine mesh must be chosen and a numerical convergence study similar to the one shown in Table 2 must be performed. For simulations during mold design (e.g. to find locations for injection gates and vents) a coarser mesh (for example mesh 2) can be chosen to keep the simulation time short. Furthermore, it can be seen that the required computational time does not only depend on the number of cells, also on the size of the smallest cells. This is a result of the adaptive time step calculated with the CFL condition [20]. According to the CFL condition, the maximum time step is proportional to the cell size and inversely proportional to the flow speed in the cell. For every cell the maximum time step must be calculated and the minimum over all cells is then taken as next time step. Summing up, the smaller the cell size for the same flow speed (i.e. for the same considered case), the smaller the time steps for the simulation.

The need for fine resolution around the injection gate comes from the implemented non-linear equation of state shown in Figure 4. This is illustrated by the following thought experiment: Consider a linear flow through a sequence of equidistant cells, the first belonging to the inlet gate, all others belonging to the preform. At the beginning of the filling the pressure gradient between the inlet and the first cell inside the preform remains unchanged until the cell is almost completely filled. Consequently there is no flow out of the preform cell into neighboring preform cells (which are still at initial cavity pressure) for some time. Then, the first interior cell of the preform is split into several small cells and the cell belonging to the inlet is also replaced by a small cell. At the beginning of the filling the pressure gradient between the small inlet cell and the first small cell inside the preform is higher compared to the first setup with large cells resulting from the smaller spatial distance. If the flow reaches the last small cell the pressure gradient is equal to the pressure gradient between the two large cells from the first setup. Since the fluid velocity is proportional to the pressure gradient according to Darcy’s law, this thought experiment shows that the filling for finer cells is faster. If the equidistant cell size is below a certain limit (which must be shown by a numerical convergence study) the filling is independent of the cell size.

3.3. Flow experiments

The verification and validation benchmark presented in this study consists of two real-world flow experiments. All test cases are derived from an optical permeameter available at Montanuniversität Leoben, see Grössing et al. [21] and Fauster et al. [22]. For one of the test cases also a verification with a different software tool was performed.

Figure 6 shows the test rig which is used to run radial flow experiments. Directly on top of the working table the metal mold half is mounted. The upper mold half is made of glass and is framed with metal profiles. The actual mold cavity is specified through the cavity frame showing an inner dimension of 300 × 400 mm. After placing the stack of reinforcing materials of dimension 290 × 390 mm inside the cavity, the mold is closed by flapping the glass plate into a horizontal position. The radial flow experiment is then executed by injecting the test fluid into the cavity through a central injection point in the metal bottom mold half. Due to the thickness of the glass plate and the clamping frame used to tighten the glass plate, a maximum injection pressure of 4 bar can be applied. A camera system, mounted 1 m above the mold, is used to acquire an image sequence during the radial flow experiment. If the experiment is used for permeability characterization, the sequence images are evaluated by means of a digital image processing algorithm specifically developed for this application. The flow front is described by an elliptical geometry model. Thus, the radially advancing flow front is finally obtained in terms of the major and minor axes length characteristics.

Figure 6. Optical permeameter for radial flow experiments.

3.4. Validation with real-world flow experiments

The preforms for the validation experiments consist of layers of woven fabric Hexcel 1202. The preform is positioned such that the first principal flow direction is aligned with the long edge of the preform. The fabric was characterized beforehand for different numbers of layers in the described optical permeameter. Also mass density and dynamic viscosity of the used test fluid at lab temperature were determined before the experiments were conducted.

Initially the cavity is filled with air and plant oil is injected as test fluid. The fluid properties are air density 1.205 kg/m³ at initial cavity pressure and oil density 960 kg/m³ at injection pressure. The measured dynamic viscosity of the oil is 0.071 Pas in the first test case and 0.062 Pas in the second test case. The differences result from temperature variations in the lab.

3.4.1. Variation of preform properties

The setup is a model of the cavity in the optical permeameter (see Figure 6), i.e. a flat plate with 390 × 290 × 3.14 mm³. The main preform covers the leftmost section (245 mm wide) with the central injection port with 13 mm diameter. The remaining section (145 mm wide) is manipulated to show lower permeability and porosity. Figure 7 shows the top view of the preform model (mesh) with the central injection port and the reinforced section.

Figure 7. Mesh with highlighted central injection gate and section with different properties than the main preform (left) and dimensions (right) for the first validation test case.

The preform consists of 11 layers of fabric Hexcel 1202 in the main section and 14 layers in the section in the right end. The fabric properties are given in Table 3.

Table 3. Preform properties (number of layers N, porosity ε, orthotropic permeability k₁ and k₂) for the first validation case (variation of preform properties).

	N [–]	ε [–]	k1 [${10}^{-12}$ m^2]	k2 [${10}^{-12}$ m^2]
Main section	11	0.604	163	50.3
Manipulated section	14	0.468	28.6	3.4

Initial cavity pressure is ambient pressure and injection is with 100,000 Pa over-pressure. Due to pressure losses in the feed line the measured pressure difference between injection gate and cavity is 91,000 Pa.

Different propagation speeds occur in the different preform zones. This can be inspected if the flow front is smooth. Consequently, for this case a fine mesh (5150 cells) was chosen.

For the validation of the simulation the flow front propagation is compared with that of the flow experiments. Three experiments were conducted and all three showed a similar filling behaviour. A visual comparison of the filling after 120 s filling time for two different experiments is used to show the reproducibility of the experiments. Figure 8 shows the difference in filling after 120 s for two experiments. The relative error in the filling areas is 2%. This difference must be considered in the interpretation where flow fronts from experiment and simulation are compared.

Figure 8. Filling after 120 s for two different experiments with the same process, permeability and fluid parameters is shown in subplots (a) and (b). The difference in the filling is shown in (c).

Figure 9(c–f) compares the filling of one representative flow experiment and simulation. Pictures from similar time instances are transparently overlaid for visual comparison. Due to the reinforcement with lower permeability and lower porosity in the rightmost section the flow front propagation slows down there. In general, if propagation speed and filling pattern are compared, experiment and simulation show good agreement. The simulated filling leads the experimental filling in horizontal direction and trails in vertical direction. This can be explained with small discrepancies between the experimentally determined volume-averaged in-plane permeability values from those of the actual preform.

Figure 9. Inlet and patch areas (a), picture from the start of the flow experiment (b) and overlaid filling from simulation and experiment for different time instances (c–f) for the first validation case (variation of preform properties). The flow front propagation slows down in the zone of reduced porosity and in-plane permeability. Simulation and experiment show very good agreement.

Race tracking was expected in a vertical channel at the boundary between the two preform zones but was not observed in the experiments.² Since race-tracking is omnipresent in RTM, an additional experiment with enforced race-tracking was used for validation.

3.4.2. Flow manipulation by local variation of preform properties

The setup is again a model of the cavity in the optical permeameter, but with cavity thickness 3.0 mm. A 2.0 mm wide gap was realized around the edges of the patch, where race-tracking is enforced. The setup is illustrated in Figure 10.

Figure 10. Mesh (left) and dimensions (right) for the second validation test case. A 2 mm wide gap was realized around the edges of the patch, where race-tracking is enforced. This zone is highlighted in the mesh.

The preform consists of 11 layers of fabric Hexcel 1202 in the main section and 16 layers in the patch. The fabric properties are given in Table 4.

Table 4. Preform properties (number of layers N, porosity ε, orthotropic permeability k₁ and k₂) for the second validation case (flow manipulation by local variation of preform properties).

	N [–]	ε [–]	k1 [${10}^{-12}$ m2]	k2 [${10}^{-12}$ m2]
Main section	11	0.583	96.6	38.6
Patch	16	0.393	1.2	15.3
Race-tracking zone	0	0.960	1500	1500

Note: The permeability values in the main preform were slightly increased compared to the values from the preform characterization to match the flow front propagation in the main preform.

Initial cavity pressure is ambient pressure and injection is with 100,000 Pa over-pressure. Due to pressure losses in the feed line the measured pressure difference between injection gate and cavity is 91,000 Pa.

The expected flow front propagation is an ellipse in the main preform and fast flow through the race-tracking channel to create an entrapment in the patch. A similar mesh as for the first test case was chosen (5308 cells).

One of the best-known software tools for predicting the flow front advancement in RTM is PAM-RTM, see PAM-RTM [23]. Therefore, the predicted flow front advancements of LCMsim are compared to the results of PAM-RTM using the same shell mesh as well as identical material properties and process parameters.

The simulated flow front from both simulation tools are compared with that of the flow experiment. Again, three experiments were conducted and all three showed a similar filling behavior. Simulation results are compared with pictures of one representative flow experiment. Pictures from similar time instances are transparently overlaid for visual comparison in Figure 11. Both simulation tools render very similar filling patterns for the shown time instances, and simulations and experiment are in very good agreement. The overall flow front propagation is elliptical. As soon as the race-tracking channel is completely filled, air entrapment is formed in the patch area. However, the filling of the patch is much faster in both simulations compared with the experiment. From a mold engineering point of view it is sufficient to see that an entrapment occurs for a certain configuration and then the injection strategy (i.e. the location of gates and vents) must be changed. From a research point of view this phenomenon must be studied in more detail.

Figure 11. Flow front from experiment (first row), overlaid filling from LCMsim simulation and experiment (second row) and overlaid filling from PAM-RTM simulation and experiment (third row) for different time instances for the second validation case (flow manipulation by local variation of preform properties).

The two simulation tools perform well for applications with flow through preforms with different permeability and porosity values. This is demonstrated for example in validation case 1 of this study and in Grössing et al. [21]. Only in case of race-tracking (which is modeled by patches with much higher permeability and porosity), the filling from the experiment is different than the filling from both simulations tools. If the preform values are adjusted to match the flow in race-tracking channel direction, there is too much cross flow in both simulation tools. Since flow through porous media is described by Darcy’s law in both simulations tools, the results of validation case 2 suggest that there is an additional flow resistance between regions with significantly different porosity and/or permeability values. Summing up, both filling simulations describe the flow through porous media based on Darcy’s law, which is used because of its simplicity. However, it is only applicable under certain conditions (isothermal, steady and laminar flow of a Newtonian fluid through a homogeneous-porous medium). These assumptions are apparently violated in this test case at the transition of the preform to the patch (with significant differences in fiber volume content, porosity and permeability) and consequently simulation and experiment are no longer in good agreement there.

3.5. Discussion

In order to compare the simulated and experimentally determined flow fronts quantitatively, the filling after 120 s for the first test case and the filling after 60 s for the second test case are analyzed in more detail.

The flow front position is characterized by n equally distributed rays starting at the inlet center and ending at the intersection with the flow front. The relative error along ray i in direction θ_i is calculated by: (7) $$e_i=\frac{l_{i,1}-l_{i,2}}{l_{i,2}},$$(7) where $l_{i,1}$ and $l_{i,2}$ denote the length of ray i with respect to the first and second type of flow front, respectively. Figure 12 shows how the lengths along ray i are measured.

Figure 12. Sketch explaining the used variables for calculating the relative error e_i along ray i.

In the first test case, the flow front predicted by LCMsim is compared with the flow experiment. In the second test case, the flow front predicted by LCMsim is compared with that predicted by PAM-RTM. In both cases, n = 37 was chosen, i.e. the flow fronts are sampled with a spacing of ${10}^{°}.$ The root-mean-square (RMS) error in the flow front determined by the two different methods is then: (8) $$e=\sqrt{\frac{1}{n}\sum _{i=1}^n{e}_i^2}$$(8)

Since in case 1 the simulated filling leads the experimental filling in horizontal direction and trails in vertical direction, the flow fronts from experiment and LCMsim have a RMS error of 7%. The flow fronts from case 2 simulated with PAM-RTM and LCMsim exhibit a deviation of 3% in terms of the RMS error. The entrapment in case 2 is characterized by the enclosed area. Black/white pictures of equal size where the black area covers the entrapment are analyzed. The numbers of black pixels are used to calculate a relative error in the entrapment area. The entrapment area is nearly identical (0.5% difference).

Figure 13 shows the analyzed pictures. LCMsim and experiment have a sharp boundary for the flow front. For PAM-RTM the position of the flow front was taken between light and dark green in the contour plot which corresponds to a filling of 50%.

Figure 13. The flow front position is characterized by 36 equally distributed rays starting at the inlet center and ending at the intersection with the flow front. (a) and (d) describe the flow fronts from case 1 after 120 s simulated with LCMsim and measured in the experiment. (b) and (e) describe the flow fronts from case 2 after 60 s simulated with LCMsim and PAM-RTM. The entrapment is characterized by the enclosed area. (c) and (f) describe the entrapment in the patch from case 2 after 60 s simulated with LCMsim and PAM-RTM.

It has been confirmed that the LCMsim software is appropriate for performing RTM filling simulations under real-world conditions. Case 1 shows that the predicted flow front advancement simulated with LCMsim agrees with those from the experiment. Case 2 shows that both simulation tools, LCMsim and the commercially available PAM-RTM, render similar results, although different modeling equations are solved.

4. Conclusion and outlook

The work at hand presents a benchmark of the software LCMsim with a new flow physics model for filling simulations in RTM and how this model was implemented in the programming language Julia. The Julia module can be downloaded from a GitHub repository https://github.com/obertscheiderfhwn/LCMsim. In general, the software is easy-to-use for the tool engineer and gives maximum flexibility for scientific investigations. The new software tool can be customized for non-standard filling applications. For example, by adding new or customizing existing equations or boundary conditions. Furthermore, since the code is open source, code snippets can be extracted and integrated in separate codes, all intermediate results are available by creating new interfaces, and the simulation parameters can be modified during a run (for example, changing inlet pressure levels as function of flow front shape).

The software LCMsim implements a new fluid mixture model (continuity equation for the mixture density, momentum equations in a local cell coordinate system, equation of state which models the pressure build-up as function of mixture density, filling fraction determined from mixture density). The implementation results in a robust algorithm without pressure-velocity coupling during one time step, only takes a shell mesh as input and can be extended easily. In order to solve the governing equations on a shell mesh, the finite area method is generalized for curved surfaces with junctions.

Verification of the presented filling model was previously done with radial flow cases and by comparison with filling results from other software tools [11]. The study at hand successfully validates the new software with real world experiments under real process conditions for RTM. The software can deal with patches with different in-plane permeability and porosity levels and with race-tracking. Furthermore, it has been shown that LCMsim and PAM-RTM yield very similar results.

Based on the already presented verification and validation cases, the next step will be an extension of the software to cover particularities of the vacuum-assisted resin infusion (VARI) process. There, transverse flow from the flow distribution media into the fibrous preform will be implemented on the shell mesh. If needed, cells with two layers representing the flow distribution medium and the preform will be provided. Moreover, models for resin temperature and viscosity as well as cure kinetics will be added to the software. Apart from improving the physical model, the comparison in Section 1.3 and Table 1 shows the need of including a mesher and improving the computational performance.

Acknowledgments

The authors thank Oliver Rausch-Schott for running the flow experiments on the optical permeameter test rig at Montanuniversität Leoben.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The novel software tool LCMsim is openly available on GitHub at https://github.com/obertscheiderfhwn/lcmsim. Installation instructions and instructions on how to run a simulation are given in a README file on the repository. The input files input_case4a.txt and input_case4b.txt can be used to derive the data presented in the two test cases.

Notes

1 In this paper, mathematical symbols and units are noted according to the ISO 80000 standard.

2 If race-tracking is present it is included in the simulation with a small zone with increased permeability and porosity. Typically, the permeability in the race tracking zone is set one or two orders of magnitude higher than the permeability in the main preform and the porosity is close to 1.

Appendix A:

code modifications

The code snippet in A1 shows the implementation of the numerical scheme as described by Equations (4)(4) $${\rho }_i^{n+1}=({\epsilon }_i{\rho }_i^n-\frac{\Delta t}{V_i}\sum _j{\mathit{n}}_{i,j}·\underset{{(\rho \mathit{u})}_{i,j}^n}{\underbrace{\frac{1}{2}\left({\rho }_i^n+{\rho }_j^n\right) \frac{1}{2}\left({\mathit{u}}_i^n+{\mathit{u}}_i^n\right)}}{A}_{i,j})/{\epsilon }_i$$(4) , (5)(5) $$u_i^{n+1}=\left({\rho }_i^n{u}_i^n-\frac{\Delta t}{V_i}\sum _j{\mathit{n}}_{i,j}·{(\rho \mathit{u})}_{i,j}^n{u}_{i,j}^n{A}_{i,j}-{\left(\frac{\partial p}{\partial x}\right)}_i\Delta t\right)/\left({\rho }_i^{n+1}+({\mu }_i/k_1)\Delta t\right)$$(5) and (6)(6) $$v_i^{n+1}=\left({\rho }_i^n{v}_i^n-\frac{\Delta t}{V_i}\sum _j{\mathit{n}}_{i,j}·{(\rho \mathit{u})}_{i,j}^n{v}_{i,j}^n{A}_{i,j}-{\left(\frac{\partial p}{\partial y}\right)}_i\Delta t\right)/\left({\rho }_i^{n+1}+({\mu }_i/k_2)\Delta t\right)$$(6) , respectively. During the time evolution (line 8) for every internal cell (index ind in line 12) the main steps are:

• The pressure gradient function (line 15) is called to evaluate ${(\frac{\partial p}{\partial x})}_i$ and ${(\frac{\partial p}{\partial y})}_i.$

• The numerical flow (e.g. ${\mathit{n}}_{i,j}·{(\rho \mathit{u})}_{i,j}^n{u}_{i,j}^n{A}_{i,j}$ for the x-velocity) at the boundaries of the considered cell i to all neighboring cells j (lines 26–67) is evaluated.

• The mass density and velocity components at the new time step (lines 69–78) are calculated according to the numerical schemes. Depending on the flow physics model, a different continuity equation is used. The original flow physics model corresponds to the blocks with if i_model==1 and the current flow physics model to the blocks with elseif i_model==2. The factors for porosity, permeability,…are equal to 1 and prepared the code for VARI with elseif i_model==3

• Depending on the flow physics model, the filling fraction and the fluid pressure are calculated (lines 80–105).

• In the evaluation of the cell boundary flow, first-order upwinding is used (e.g. lines 132–137 for the calculation of $u_{i,j}^n$ in the x-velocity flow ${\mathit{n}}_{i,j}·{(\rho \mathit{u})}_{i,j}^n{u}_{i,j}^n{A}_{i,j}$) depending on flow out of (with $\ge$) or iinto the considered cell. The contribution n_dot_rhou at boundaries to pressure inlet and outlet cells is calculated differently (lines 129 and 171) with the Darcy velocity (line 56) or the outflow velocity (line 58). This difference is the reason for two different functions. Different (higher-order) methods for the evaulation of the cell boundary flow can be added with the function argument i_method.

If an additional equation, e.g. for the degree-of-cure α, is added, the following steps must be performed:

• Additional parameters for the initial and inlet degree of cure (alpha_init = 0 and alpha_a) and for the cure kinetics model (e.g. k and n for the n-th order reaction model $S_{\alpha }=k{(1-\alpha )}^n$ are defined. Additional arrays alpha_new and alpha_old must be initialized similar to rho_new and rho_old. Instead of hard-coded parameters, other sections in the code can be modified to take these as input parameter. Directly before the time evolution section begins, the additional arrays and parameters are prepared:

• If the degree-of-cure equation is added to the new flow physics model, replace all occurances of if i_model==2 by if i_model==2 ǁ i_model==4 since the new flow physics model with degree-of-cure is identified with i_model = 4.

• After the equation of state section (after line 105) add a section inside an if i_model==4 block, where the degree-of-cure equation

(A1) $$\frac{\partial \epsilon \alpha }{\partial t}+\mathit{u}·\nabla \alpha =S_{\alpha }$$(A1)

in a conservative form (A2) $$\frac{\partial \epsilon \alpha }{\partial t}+\nabla ·\left(\alpha \mathit{u}\right)-\alpha \left(\nabla ·\mathit{u}\right)=S_{\alpha }$$(A2)

is solved. Therefore, the code snippet of the filling fraction equation of i_model = 1 which is also a transport equation is copied (lines 81–83). Gamma is replaced by alpha and a source function is added:

• Also in the functions for the calculation of the flow at the cell boundaries, the sections for gamma (lines 145–152 and 186–193) are copied and all gamma occurrences are replaced by alpha. In addition, alpha_P and alpha_A must be added to the vector arguments vars_P and vars_A. alpha_P = gamma_old[i_P]; and alpha_A = gamma_old[i_A]; are assigned just before the vectors are formed (before line 41). In the flow functions they must be assigned by alpha_P = vars_P[5]; and alpha_A = vars_A[5];.

• At the end of the flow definition (after line 25) the variables for alpha must be defined:

  •  F_alpha_num = Float64(0.0);F_alpha_num_add = Float64(0.0); and

  •  F_alpha_num1 = Float64(0.0);F_alpha_num1_add = Float64(0.0);

A1 Code snippets