Fast, Accurate, Reliable 3D
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.
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.
Simulation of flow in RTM and Vacuum Infusion on a macroscopic level is well described by Darcys Law. It is our starting point too and pretty straightforward to implement [KOOR02]. However, for a simulation code Darcys 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.
The flow of resin through reinforcement is well described by Darcys Law [TUCK94]. With gravity terms included, Darcys Law generalized to three dimensions and the continuity equations are:
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 Darcys 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:
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:
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.
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:
In addition density is assumed to be constant (the resin behaves as an incompressible fluid), which leads to the well-known Navier-Stokes equations:
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):
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 areas where 3D volume elements are typically needed and available - contain reinforcement.
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):
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:
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:
Integration of equation (3) with gravity forces included over the thickness leads to a similar result for the RTM shell element:
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:
Inspection shows that in this very particular case, the GHS formulation is equivalent to the RTM formulation with an apparent permeability of k = H2/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:
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.
Application of the Galerkin Finite-Element method and partial differentiation results in the following set of non-linear equations:
Here, fi 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 fi and bi represent body forces (gravity) and boundary conditions; bi is only nonzero in nodes where pressure or flowrate is prescribed. Note that the stiffness matrix coefficients Kij 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):
Because only the gradient of the pressure enters the equations, absolute pressure has to be fixed by prescribing it somewhere on the boundary.
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 Tadmors 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 Kij contain the viscosity which depends on temperature and degree of reaction (and in the most general case also on shear-rate).
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:
In this equation, x3 is the local coordinate in thickness direction, perpendicular to the in-plane pressure gradient. The degree of conversion is calculated from:
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).
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:
The time-derivative is discretised using an implicit (unconditionally stable) first order backward difference scheme:
Integration of the evolution equation for the degree of conversion, multiplied with weight functions w(x,t), results in:
In contrast with the equations for the pressure problem (and temperature conduction), the stiffness matrix is asymmetric and the Galerkin method (choice wi = Di 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 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 wi are chosen such that Mij and Kij become zero if pi > pj. The weighing functions are defined as follows:
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.
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:
With Tic the convected temperature, which is the solution of:
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( Tic, ai ). 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.
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.
For the viscosity, the model proposed by Castro & Macosko [CAST82] has been implemented in RTM-Worx:
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.
The energy released by the chemical reaction is proportional to the rate of conversion (by definition):
The reaction rate Ra 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]:
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.
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):
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.
flow model |
|beams and shells
|Products and Services
7443 EH Nijverdal
+31 548 612217