Fast, Accurate, Reliable 3D
Reactive RTM Simulation

The contents of this page were presented by Arjen Koorevaar as a paper at the ISCM 2002 conference, organised by the Dutch Aerospace Laboratories (NLR) in Vollenhove, May 2002.

Abstract

The isothermal flow 2.5D simulation technology has matured to a point where accuracy and speed of the calculation are balanced and behaviour of the model is very predictable and easy to understand without detailed knowledge of the numerical model. RTM-Worx has been extended to include (1) calculation of temperatures and resin conversion and (2) 3D volume elements compatible with existing beam (1.5D) and shell (2.5D) elements. The added complexity and specific issues related to transport of energy by flow in a fixed mesh (Eulerian discretisation), variable viscosity (which tends to destabilise the calculations) and huge number of unknowns (3D) necessitated the development, implementation and extensive testing of new numerical methods. This paper details those methods and how we achieved the main design goals: balanced (high) speed and (high) numerical accuracy while the much more complex code behaves with a predictability comparable to the much more simple isothermal flow model.

Introduction

Simulation of flow in RTM and Vacuum Infusion on a macroscopic level is well described by Darcy’s Law. It is our starting point too and pretty straightforward to implement [KOOR02]. However, for a simulation code Darcy’s Law is not sufficient because it is only valid for flow through porous media. In practice there are tubes (spirals) and channels in the mould that do not contain reinforcement. While an equivalent permeability (based on comparison with analytical formulas for Poiseuille or Hele-Shaw flow) can be used in the isothermal case, we need the 3D-velocity field to be able to predict temperature and reaction rate. Because the RTM-Worx software shares the elements with applications for (gas-assisted) thermoplastic injection moulding and transfer moulding (for chips encapsulation for example), the formulation is more general than required for RTM and Vacuum Infusion.

3D Flow Through Porous Media

The flow of resin through reinforcement is well described by Darcy’s Law [TUCK94]. With gravity terms included, Darcy’s Law generalized to three dimensions and the continuity equations are:

(1)

(2)

Here, u is the local flux density (or superficial velocity), K is the permeability tensor, · is the resin viscosity, p is the pressure in the resin, Á_r is the (local) resin density and g is the gravity vector. Substitution of Darcy’s Law (1) in equation (2) results in a formulation of the continuity equation where the only unknown is the pressure (a 3D scalar quantity) within the fluid in the mould:

(3)

This is a 3D formulation of the so-called ‘pressure problem’ for flow through porous media.

With the assumption that thermal contact between fabric and resin is ideal (e.g. it is allowed to use a continuum approach in which resin and fabric properties are averaged), conservation of energy can be formulated as:

(4)

Where › is the effective conductivity tensor, with contributions – averaged in some, yet unspecified way - from reinforcement and resin. Note that pressure and energy equation are tightly coupled because the velocities in the energy equation directly follow from calculated pressures, for which the equation contains viscosity that depends on temperature and degree of conversion. In other words, the problem is highly non-linear and solving both equations simultaneously is not a trivial task.

Generalized Newtonian Model

The model for flow through sections without reinforcement is based on the so-called Generalised Newtonian Fluid. This means that elasticity and bulk viscosity of the fluid are neglected and viscosity only depends on pressure, temperature, degree of conversion and (for shear thinning/thickening fluids) on the deformation rate tensor:

(5)

(6)

(7)

In addition density is assumed to be constant (the resin behaves as an incompressible fluid), which leads to the well-known Navier-Stokes equations:

(8)

Conservation of energy can be described by an equation very similar to the one used before (eq. 4), but with inclusion of the dissipation term (which is actually not significant for RTM and Vacuum Infusion):

(9)

(10)

Solving those equations requires a formulation with both pressures and velocities as unknown quantities (or definition of a stream function), which makes it very difficult to combine the elements with the pressure-only elements that are used to solve Darcy flow. In addition, solving those equations is very time consuming and difficult to implement in a software targeted to users that not necessarily have the detailed knowledge needed to deal with the intricacies of solving the Navier-Stokes problem in 3D. Therefore, we only consider shell and beam elements for which a much more efficient model can be formulated. The advantages - to combine the remaining element types in any configuration and calculate at very high speed - largely outweighs this limitation. In practice, this turns out to be no limitation at all because the area’s where 3D volume elements are typically needed – and available - contain reinforcement.

Unified GHS/RTM Model Formulation

The restriction to thin-walled sections makes it possible to introduce a number of simplifications in the Navier-Stokes equations which lead to the so-called Generalised Hele-Shaw (GHS) model [HIEB80], better known as the 2½D model (see [BOSH88] for a detailed description):

• The pressure gradient in thickness direction can be neglected because it is much smaller than the pressure gradient in the (local) plane.

• Due to the previous assumption, velocity is zero in thickness direction.

• Conduction of heat in plane can be neglected because temperature differences in thickness direction are much larger compared to the temperature gradient in plane.

• Velocity and temperature are assumed to be symmetrical. While it is possible to formulate the model for asymmetrical cases, it is not very practical because it introduces additional difficulties, like in the case of branching flow [ZOET95].

As a consequence of those simplifications, the pressure calculation reduces locally to a 2D problem and convection and conduction terms in the energy equation can be solved decoupled from each other. However, in comparison with the model as formulated in [HIEB80], there are two major differences in the model presented here:

• Gravity forces are not neglected. They are significant for Vacuum Infusion because the pressure gradient is relatively small and the products typically quite large (boat hulls for example).

• The velocity in thickness direction is recovered from the mass balance and the resulting true 3D velocity field is used to solve the convection problem. This makes it possible to model phenomena like the fountain flow at the front (in the absence of reinforcement) and two-phase flow (like the combined flow of nitrogen gas and polymer in the gas-assisted injection moulding process).

Analytical integration of the Navier-Stokes equations over the thickness, leads to the following, much simpler, formulation for the pressure problem on a GHS shell element:

,

(11)

(12)

Integration of equation (3) with gravity forces included over the thickness leads to a similar result for the RTM shell element:

(13)

(14)

This analysis naturally results in a rule of mixtures for the permeability, for the case where the viscosity is constant. Note that for the non-isothermal case, it is not correct to average permeability and viscosity separately and use the quotient of the averages instead of S* (as was done in [BRUS92]) because it will give wrong results!

The major difference between equation (11) and (13) is that the flow conductivity S is a tensor, but only if the permeability of the reinforcement is not isotropic. In the general case, the directions of the pressure gradient and the velocity are not necessarily the same. If the permeability is isotropic (and equals k in all directions), S will be a scalar value. If in addition the viscosity is constant, the flow conductivity S in equations (12a) and (14a) simplifies to:

and

(15)

Inspection shows that in this very particular case, the GHS formulation is equivalent to the RTM formulation with an apparent permeability of k = H²/12. This is the previously mentioned ‘equivalent permeability’ and it is clear now why the GHS flow element can not be replaced with RTM flow elements in general. For 3D RTM volume elements, straightforward reformulation of equation (3) results in:

(16)

(17)

The formulation of the beam elements is very similar (derivation is left as an exercise to the reader), with the difference that the local tensors are one-dimensional. We now have a generalized model for the calculation of pressure (and resulting velocities) in 3D volume (Darcy flow only), 2D shell (Darcy and Hele-Shaw flow) and 1D beam (Darcy and Poiseuille flow) elements.

Finite Element Discretisation of the Pressure Problem

Application of the Galerkin Finite-Element method and partial differentiation results in the following set of non-linear equations:

		(18)
		(19)
		(20)

Here, f_i are the interpolation functions on the elements. In RTM-Worx, linear (line, triangle and tetrahedron) elements are used. The formulation of those elements can be found in [HUGH87] for example. The right-hand side vectors f_i and b_i represent body forces (gravity) and boundary conditions; b_i is only nonzero in nodes where pressure or flowrate is prescribed. Note that the stiffness matrix coefficients K_ij contain the viscosity and therefore depend on temperature, degree of reaction and (in the most general case) on shear rate.

Once the pressures have been obtained by solving eqs. 18-20, the resin flow velocities can be calculated (resin density assumed constant):

GHS 2D shell:		(21)
RTM 2D shell/3D tet:		(22)

Because only the gradient of the pressure enters the equations, absolute pressure has to be fixed by prescribing it somewhere on the boundary.

Flow Front Tracking, Control Volume Method

The mould filling process is a transient problem, and the free surface that separates the part that has been filled with resin from the still empty (dry) part has to be defined. Several solutions have been published; a good overview is given in [TUCK89]. In RTM-Worx, the so-called ‘Volume of Fluid’ method is used, introduced in [HIRT81]. It is based on similar ideas as Tadmor’s FAN method, developed for regular Finite Difference grids [TADM74]. It has been shown [WANG86] that a Galerkin FEM formulation on linear elements is identical to a Control Volume method.

The flow-front tracking method used in RTM-Worx, and the equivalence of the FEM and CV formulations used have been published before [KOOR02] and will be summarised briefly.

In all nodes, the filling state of the associated Control Volume is defined by a fill factor f which defines the fraction of the volume filled with resin relative to the amount of resin the Control Volume can contain. This fill factor is zero for empty nodes and one for full (wetted) nodes. All nodes that are partly filled (0 < f < 1) are front nodes.

For a given state of fill, the pressures are calculated with prescribed pressure in the front nodes: ambient pressure or in case of vacuum (assistance) the vacuum level. From the pressure solution, the flow rates at the front can be easily calculated by a matrix multiplication (eqs. 18-20) which is a direct consequence of the equivalence of the FEM and CV formulations. Multiplication of flow rates with the time step increment gives the amount of resin that enters each front node, which is used to update the fill factors. This is repeated until the mould is full.

Because this is an explicit scheme, only one node can be filled in each time step increment; which defines the size of the time step increment. The total number of time steps to fill the mould is therefore equal to the total number of nodes in the mesh.

For the isothermal case, the stiffness matrix K (eq. 18) only needs to be assembled once. However, in the non-isothermal reactive calculation it needs to be reassembled each time step because the coefficients K_ij contain the viscosity which depends on temperature and degree of reaction (and in the most general case also on shear-rate).

Adding a Dimension to Beams and Shells

Calculation of temperature, degree of conversion and pressure is coupled through the viscosity. Those equations are solved decoupled in an iterative scheme. The energy equation (4) can be simplified for the 2.5D shell elements, using the previously formulated assumptions:

(23)

In this equation, x₃ is the local coordinate in thickness direction, perpendicular to the in-plane pressure gradient. The degree of conversion is calculated from:

(24)

For storage of temperature and degree of conversion, beam and shell elements need an additional grid in radial and thickness direction respectively. The grid can be defined at the center of the elements, which has the advantage that for the convection only the topology is important [LIAN93], the dissipation and source terms can be calculated locally on the elements. However, this has the disadvantage that integration of the convection terms has to be done (partly) explicitly, which limits the size of the time step increment by the Courant number. Effectively, this means that choice for the time step is dominated by the smallest elements with the highest velocities: typically the resin supply tubes. This time step is typically much smaller than needed for sufficient accuracy and thus only slows down the calculation. In addition, results from several elements need to be averaged to define nodal values, which affects accuracy in the form of numerical (artificial) diffusion.

Therefore, temperature and degree of conversion are stored in a through thickness grid at the element vertices. This makes it possible to use a very efficient implicit method for integration of convection and diffusion terms in the energy equation. The viscosity and through thickness velocity however, are defined on a grid at the center of the elements because it is approximated by (or results in) a constant value in thickness direction on each element in the pressure problem. For the calculation of viscosity, the averages of values from the grid at the element vertices for temperature and degree of conversion are used (figure 2).

Reaction Kinetics, Convection

Now that the grid has been defined, the degree of conversion can be represented with element shape functions that interpolate it linearly on the elements (separately in each layer) according to the Finite Element Method:

,

(25)

The time-derivative is discretised using an implicit (unconditionally stable) first order backward difference scheme:

(26)

Integration of the evolution equation for the degree of conversion, multiplied with weight functions w(x,t), results in:

	(27)

In contrast with the equations for the pressure problem (and temperature conduction), the stiffness matrix is asymmetric and the Galerkin method (choice w_i = D_i for the weight functions) leads to a numerically unstable algorithm. A lot has been published on this subject, see [TUCK89] for an overview and more specifically about the solution of the convection problem in the 2½D problem [ZOET95]. The discretisation error manifests itself as artificial diffusion. In the case of the Galerkin formulation those diffusion terms are negative which results in oscillations in the solution. What is needed is a method that adds just enough diffusion to stabilize the calculation, which also implies that the added diffusion has to depend on the flow direction. Our design goals are:

• The calculation has to be robust, reliable and unconditionally stable.

• Refinement of the mesh must result in a more accurate solution.

• The size of the time step increment should not be limited.

• The algorithm has to be fast and use a modest amount of memory.

• Accuracy should match solution of pressures and velocities.

• Implementation should not be very complex and generalisable to 3D.

The first five are very important for the user-friendliness of the software. The user should not be concerned with the theory and implementation details, especially because instabilities can be hard to detect – it is difficult enough to interpret results from a transient non-isothermal reactive flow simulation. The last goal ensures that the chances for programming errors (bugs) are minimized, especially when future extensions will be implemented (and a lot of details about the current implementation are easily forgotten or overlooked).

The major source for inaccuracy in the convection equation is that field (continuous) information is represented by linear interpolation between a limited number of locations where values are stored in combination with a certain inaccuracy in the velocities used to transport those quantities. Regardless of the cleverness of the algorithms, there is only one effective way to obtain higher accuracy: increase the number of nodes. Refining the mesh also increases accuracy of calculated velocities. The algorithm should make optimal use of all available information, but produce sensible results on a course mesh too.

In RTM-Worx, the choice has been made for the so-called ‘Strong Implicit Upwind’ method [BAKHxx], based on a Petrov-Galerkin method in which the weighing functions w_i are chosen such that M_ij and K_ij become zero if p_i > p_j. The weighing functions are defined as follows:

• equal to D_i for elements upstream to the node (p_i < p_j for all nodes j);

• equal to zero for downstream elements (p_i > p_j for all nodes j);

• a weight function for all the other elements defined according to the formula below (eq. 28) in which the coefficients a_ki and b_ki are chosen such that M_ij and K_ij become zero for all downstream nodes in the element:

(28)

If the nodes are properly sorted, the stiffness matrix will become triangular: it only contains non-zero terms above or below and including the diagonal, which makes it trivial to solve the linearised system of equations. In RTM-Worx, the matrix is actually never assembled; all the equations are solved locally node-by-node. The amount of work is comparable to a simple upwind method, but this method is second-order accurate and does not introduce diffusion (e.g. inaccuracy) perpendicular to the direction of flow. The vector u is calculated consistently, with the same weighing functions as the convection terms.

The rate of conversion can be calculated explicitly; a second order Runge-Kutta method (also known as midpoint method) is used to update the degree of reaction from the values obtained by convecting the results from the previous time step.

Temperature calculation

In comparison with the reaction kinetics, the energy equation contains additional conduction terms (source terms are calculated explicitly using the results from the previous time step). The temperature is calculated with [CASP95] using a two-step algorithm:

(29)

With T_i^c the ‘convected’ temperature, which is the solution of:

(30)

In words, this comes down to a scheme in which the first step consists of moving all particles to their new location (after the front has been updated), followed by the remaining part as in a Lagrangian discretisation (in which the mesh is tied to the material). Decoupling convection from the remaining part of the calculation is a form of particle tracking. If, for example each particle that enters the mould was labeled, convection of the labels allows us to identify each particle at each instant in time (such a calculation, in which the particles are labeled with time of entry is actually a standard feature of the RTM-Worx non-isothermal reactive extension module).

The source term is calculated as š = f( T_i^c, a_i ). This only states that the properties of the resin depend on temperature and degree of cure; any model (including a zero value) can be used without consequences for the implementation or speed of the calculation.

Material models

To complete the model, constitutive equations are needed to define the dependence of viscosity on temperature and degree of conversion and the reaction kinetics. In RTM-Worx, the models for viscosity and reaction kinetics are easily extended or replaced (as explained in the last part of the previous paragraph). However, the models already available will suffice in most cases.

Viscosity

For the viscosity, the model proposed by Castro & Macosko [CAST82] has been implemented in RTM-Worx:

(31)

Several variants, simplifications actually, can be represented by this model by using specific values for some of the coefficients as summarized in the table below.

Reaction kinetics

The energy released by the chemical reaction is proportional to the rate of conversion (by definition):

(32)

The reaction rate R_a is nothing else than da/dt, the time derivative of the degree of conversion of a particle. RTM-Worx uses the following autocatalytic model for the reaction rate [KENN90]:

	(33)

This model contains a large number (9!) of coefficients and will therefore have sufficient flexibility to fit data from DSC measurements. A number of simplified models can also be represented by this model, as summarized in table 2.

In general, the best choice is the simplest model (with the least number of coefficients) that is capable of modeling the reaction kinetics with sufficient accuracy.

Implementation Notes

The model presented in this paper has been implemented in RTM-Worx. Actually, RTM-Worx has been the test vehicle for a lot of variants of this model and this is what turned out closest to our original goals. There are several other issues we had to deal with, which are beyond the scope of this paper but mentioned here for completeness (some of this material will be published in future papers):

• At the flow front in the GHS formulation (beam and shell elements), a phenomenon called ‘fountain flow’ occurs. This is caused by the flattening of the velocity field and has to be taken into account because the velocities in x₃ direction can not be neglected in this region. In RTM-Worx, no special measures have been taken to deal with the front; we simply included convection in x₃ direction everywhere using the velocities recovered from the mass balance.

• If the non-isothermal reactive calculation is used to do an isothermal simulation (e.g. constant viscosity and zero reaction rate), the results should be equal to an isothermal simulation. This is not trivial because the non-isothermal reactive module uses numerical integration through thickness where the isothermal module uses the analytical solution of those integrals. Therefore, special numerical integration schemes had to be developed to make sure that both models produce exactly the same results for the same problems.

• Nodes need to have their temperatures and degree of cure initialized when they become part of the flowfront. This was solved in a very elegant way using the generalized 3D-convection algorithm.

• Boundary conditions for temperature and degree of cure must be specified at injection nodes, but only as long as flow occurs. Again, this has been included as part of the convection algorithm. When the flow stops, the boundary condition is no longer needed and vanishes naturally because the convection terms become zero. This makes it possible to use one single model for both the filling and post-filling (curing) stages in the process.

• The material models contain terms that easily become extremely small or very large. Special care had to be taken to make sure that the software does not crash, even if the results do not make sense, and fails gracefully in case there are really no possibilities left to continue.

Conclusion

It may be clear that the non-isothermal reactive model for RTM and Vacuum Infusion with elements for Generalised Hele-Shaw (and Poiseuille) flow is far more complex than the 3D model limited to Darcy flow presented in an earlier paper [KOOR02]. However, all the ingredients are necessary to make an implementation useful for practical purposes. In a sense, the model presented here is the simplest one that includes all the phenomena important in RTM and Vacuum Infusion and combines that with straightforward easy to understand behaviour without the need study the model at the level of the implementation details in order to use it. In addition, the calculation speed is very high; a non-isothermal reactive calculation only takes about twice the CPU time compared to the (very fast) isothermal calculation (in the order of a minute on a regular PC).

Another result of this development is a significant improvement (and better understanding) of the isothermal flow module in RTM-Worx. The resulting modules form the basis of a true user-friendly practical software code that belongs in the toolbox of every engineer involved with development of RTM and/or Vacuum Infusion processes.

References

Symbols

A	m2	Cross-area
A1, A2	1/s	Fitted coefficients in kinetic model
bi	m3/s	Prescribed nodal flowrate (injection point)
B0, B1, B2	–, K-1, K-2	Fitted coefficients for extended kinetic model
Cij	J/(m3·K·s)	Conductivity matrix for temperature calculation
Cp	J/(kg·K)	Specific heat capacity
D	1/s	Deformation rate tensor
E1, E2	J/mol	Fitted coefficients in kinetic & viscosity models
f	–	Filled fraction of Nodal Control Volume, 0 ≤ f ≤ 1
	N/kg	Gravity vector
Hm	W/(m2·K)	Heat transfer coefficient in Biot boundary condition
I	–	Unity tensor
Kij	m3/(Pa·s)	Fluidity matrix (pressure problem)
M	–	Mass matrix in convection equation
Mij	J/(m3·K)	Mass matrix in temperature calculation
m,n	–	Fitted coefficients in kinetic/viscosity models.
	M	Normal, = 1.
P	Pa	Pressure
Qi	m3/s	Calculated residual flowrate in nodes
Ra	1/s	Reaction rate
ri	J/m3	Source term in energy equation
S	m3/(Pa·s)	Fluidity (averaged over thickness in 2½D model)
	J/(m3·s)	Specific energy released in chemical reaction
t	S	Time
T	K	Temperature
	m/s	Velocity vector
V	m3	Volume (VI for the Control Volume of node i)
wi	–	Modified weight function in Petrov-Galerkin FEM
z	M	Locale coordinate in thickness direction

Greek symbols

Indices

home |  abstract |  introduction |  flow model |  generalized newtonian
unified GHS/RTM |  finite elements |  control volumes |  beams and shells
reaction kinetics |  temperature calculation |  material models
implementation notes |  conclusion |  references |  symbols

© 1997-2022 Polyworx