- Email: [email protected]
A numerical model to predict fiber tow saturation during liquid composite molding
Composites Science and Technology 63 (2003) 1725–1736 www.elsevier.com/locate/compscitech
A numerical model to predict fiber tow saturation during liquid composite molding Pavel Simacek, Suresh G. Advani* Department of Mechanical Engineering and Center for Composite Materials, University of Delaware, Newark, DE 19711, USA Received 22 March 2002; received in revised form 20 June 2002; accepted 25 June 2002
Abstract A dual scale porous medium contains two distinct scales of pores. We will consider the case in which small scale pores are ordered within well-defined sub-regions. Typical examples of such media are textile preforms used in various composite-manufacturing processes. Rigorously, the phenomenon can be modeled by using Darcy’s law to describe flow through porous medium and mass conservation for the flow within the larger pores. The smaller pores can be included within these equations as a sink term. This approach, though straightforward, poses implementation difficulties. In this paper, we suggest an alternative way to model this concept. We use the standard finite element/control volume approach and model the ‘‘internal’’ sink term by appending extra onedimensional elements to control volumes associated with the control volumes of discretized part geometry. This approach offers two advantages over previously attempted schemes: (1) the problem to be solved remains linear and flow can be calculated explicitly within the time domain and (2) existing simulation packages for RTM filling simulation will be able to incorporate saturation effects to simulate the flow in dual scale media. To illustrate this point we present the implementation of the developed algorithm within the frameworks of an existing simulation package LIMS. The implementation is capable of providing saturation data during filling of arbitrarily shaped part and can capture the influence of saturation on filling pressure or flow-rate. This should prove useful in determining the time needed to completely saturate the fiber tows, which is crucial for part performance and also allow us to explain the variation in injection pressure and difficulties in ‘‘assigning’’ single scale permeability for such porous medium. # 2003 Elsevier Ltd. All rights reserved. Keywords: Liquid composite molding
1. Introduction A dual scale porous medium can be defined as a porous medium that contains two distinct scales of pores. We will consider the case in which small scale pores are concentrated within dispersed, well-defined sub-regions encompassed by the larger pores. This is true when the porous medium has an ordered structure. Typical examples of such media are textile preforms used in various composite-manufacturing processes. The fluid flow in such media will saturate the larger pores faster than the smaller pores. To describe the physics of such flows one may still use the traditional governing equations such as the Darcy’s law and continuity equation. The Darcy’s law can be expressed as
* Corresponding author. Tel.: +1-302-831-8975; fax: +1302-8313619. E-mail address: [email protected] (S.G. Advani). 0266-3538/03/$ - see front matter # 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0266-3538(03)00155-6
hv i ¼
K rp
ð1Þ
Here hvi is the volume averaged flow velocity, rp is the pressure gradient, is the viscosity of the fluid and K describes the permeability of the porous medium. If the porous medium was homogenous with one intrinsic permeability, the continuity (mass conservation) equation results in r hv i ¼ 0
ð2Þ
The substitution of Eq. (1) in the continuity Eq. (2), results in the following governing equation [1]: K r rp ¼ 0 ð3Þ However, linear injection experiments with woven preforms have shown that the flow rate and the pressure drop relationship for such preforms is not linear as forecasted by Eq. (1) [2,3]. One must modify the Eq. (2)
1726
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
to account for impregnation of fiber tows when one has such a dual scale porous medium. Many structured (i.e. woven or braided) fibrous preforms have essentially the earlier-mentioned dual-scale geometrical structure [4,5]. This structure consists of the system of fiber tows on the macroscale level and the system of fibers on the microscale level as shown in Fig. 1. The empty pores between the fiber tows are of the order of millimeters and the space between the fibers is of the order of microns. The single scale approach essentially neglects any flow on the microscale, i.e. the details at the fiber tow level. This seems to be reasonable assumption for saturated flows, when the fluid has already filled all the pores [2,6,7]. The situation with unsaturated flows when some of the pores are empty, in the proximity of the moving flow-front during the filling process, is different. As the preform consists of two scales of pores (L1 and L2), the pressure driven flow within the preform is bound to fill the larger pores preferentially. Thus it must be considered as a dual scale process. This notion introduces the accepted terms of ‘‘saturated’’ and ‘‘unsaturated’’ permeability together with the concept of ‘‘partially saturated length’’ [8]. The macroscopic geometry of fibrous preform governs the flow at the ‘‘global’’ or macro- scale. However, on the microscopic level, the fiber tows contain both fibers and empty spaces that are being filled by the resin. The flow into these pores between the fibers continues even after the ‘‘macroscopic’’ flow front has moved beyond them. The partially saturated region is created due to this dual scale nature and is the region where the
macroscopic front has reached and occupied the spaces between the fiber tows but the pores between the fibers in a tow are not yet completely filled. The difference between the idealized flow pattern and the dual scale filling process is shown in Fig. 2. Note that, for the given fluid flow-rate, the macroscopic flow-front may move faster when dual-scale effects are incorporated. This phenomenon is, somewhat artificially, treated by using the ‘‘unsaturated’’ permeability [8]. This behavior is of practical significance. The complete saturation of fiber tows is imperative for part performance. To achieve complete saturation, the injection should continue to bleed out of the vents until the resin occupies the regions not only between the fiber tows but also within the tows. The duration to continue to bleed the resin through the vent is hard to estimate, as little systematic modeling or experimental data is available [9,10]. One can model the flow physics by introducing a subtle change in the governing equation. One can still retain the description of the macro-flow using Darcy’s equation. However, the mass conservation or the continuity equation should account for the loss of resin to the fiber tows at the microscale. This can be done by the incorporation of a sink term q within the continuity equation. Thus Eq. (3) will now assume the form [10]: K r rp ¼ qðp; sÞ ð4Þ
Fig. 1. Dual scale pores in textile preform.
Fig. 2. Idealized and dual scale flow in textile preforms.
The sink term magnitude will depend on the pressure surrounding the fiber tow and the fraction of the pores filled by the resin within a fiber tow. This is usually denoted by a saturation index s, which represents the fraction of empty spaces within a fiber tow occupied by the resin. When one uses Eq. (3) for the situation described by Eq. (4) to measure permeability, it leads to the claims of variable permeability, which are, some-
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
what artificially, solved by the introduction of ‘‘unsaturated’’ permeability within the range of the ‘‘unsaturated length’’. Note that, once the fiber tows within preform are fully saturated (s=1), the sink term q(p,s) will disappear. Then, Eq. (4) is reduced to the original form (3). When the fiber tows saturate at the same rate as the spaces between the fiber tows, then Eq. (2) describes the correct physics. However, if the tows saturate at a slower rate, Eq. (4) must be used. The form describing q(p,s) will vary depending on the perform type and possibly also on some other factors such as the filling history. As an approximation, the fiber tow itself can be considered a Darcian medium with permeability KT and the dependence of flow q on surrounding pressure p will be linear. Dependence of q on the fiber tow saturation s will vary depending on the tow cross-section and the filling history. Accurate description may be very difficult, especially as the properties of individual fiber tows may vary. We will simplify this relation by a step function: q ¼C:p; s < 1 0; s ¼ 1
ð5Þ
Thus, once the tow is saturated, the governing Eq. (4) returns to the traditional form (3). The value of constant C can be characterized by matching global flow rate vs. pressure drop behavior of the preform or from the geometry and dimensions of fiber tows. The dimensions of fiber tows are usually very difficult to measure accurately. It is difficult to include the Eqs. (4) and (5) in existing numerical simulations. Hence, we suggest an alternative approach, which is easier to model, is consistent with the physics and can be easily integrated with numerical models for fluid motion in single scale porous medium.
2. Modeling approach 2.1. FE/CV method As the modeling of impregnation of the fluid into porous medium involves a moving boundary, the finite element/control volume (FE/CV) solution scheme to simulate the filling process has served well to capture this physics efficiently [1]. The solution domain is meshed with a fixed finite element mesh. Control volumes are associated with each mesh node or alternatively every element [Fig. 3 (a)]. Each control volume has a fill factor associated with it. This factor ranges between zero (empty CV) and one (filled CV) and designates how much of the porous volume is already filled with the fluid [Fig. 3 (b)]. The pressure in the empty control volumes is known to be equal to that of the vent and the pressure in the filled control volumes is evaluated by the finite element method, using Eq. (3). Then, the flow Qij between indi-
1727
vidual control volumes i and j is determined using the computed pressure field [Fig. 3 (c)]. ð K ð6Þ Qij ¼ ni : rpdsij sij Once the flow rates are known, flow is advanced by explicit integration in time domain. The time step is selected so as to fill at least one additional control volume. This changes the fluid domain and hence the boundary conditions. The pressure solution is sought for the new domain and this process is repeated until the complete porous medium is saturated. This approach is relatively mature and widespread [1,11–16], fairly simple and quite robust, and its utilization stretches into the fields of optimization and control [17–19]. More importantly, if ‘‘incremental’’ system matrix decomposition is used as outlined in [20], the solution is very fast, as the set of system equations is solved only once for the complete filling cycle. The efficiency comes with some restrictions. The incremental algorithm [20] requires the system matrix not to change during the process, which necessarily means that the process must remain linear. If this is not the case, other efficient solution techniques are available [21–26]. For the global flow modeling, the porous fiber tows can be replaced by an impermeable medium [6,7]. This will provide the proper model for Darcy’s equation. Then, the sink term [Eq. (4)] can be added and removed only after it has consumed the proper amount of resin to saturate fiber tows within the appropriate volume. To accomplish this objective rigorously, one would need to modify the governing equation to Eq. (4). This involves change in whatever simulation package is used. This change may be unmaking of the FEM/CV algorithm described earlier, as the sink term depends on the unknown pressure and changes with fiber tow saturation (and time). Depending on the actual form of q(p,s) it may create non-linearity in the set of governing equations. Moreover, even when the form is carefully selected to circumvent this, the invariance of system matrix with time is violated and the very efficient algorithm [20] cannot be applied. We introduce an alternative approach. Instead of providing the sink term in the domain internally by modifying the equations, it may be added externally in the form of extra volume associated with each unit of the discretized model geometry, such an finite elements or control volumes. The additional volume represents the porous volume in fiber tows and must be compensated by reducing the original volume. To do so, the porosity might be decreased (fiber volume fraction is increased), so as to make the fiber tows impermeable. The extra volume is appended to nodes or control volumes in the form of 1D (runner) elements, as shown in Fig. 4.
1728
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
Fig. 3. FEM/CV approach to filling simulation: (a) domain discretization (b) flow-front and fill factors (c) flow between control volumes.
This approach is relatively easy to accomplish with CV/FEM model described earlier by following these steps: 1. One-dimensional (bar) elements representing volume within fiber tows are connected to each original ‘‘master’’ node/control volume. This operation creates ‘‘slave’’ nodes and control volumes. This process is shown in Fig. 4. 2. The porous volume within fiber tows has to be extracted from the original elements. Thus, either the volume or the fiber volume fraction within the original mesh is modified. Increasing the fiber volume fraction seems to be the easiest possibility, especially for three-dimensional meshes. 3. Simulation is run as usual. The fill factor in ‘‘master’’ control volumes represents fill factor in
macro pores between the fiber tows, while the fill factor in ‘‘slave’’ control volumes, connected by bar elements, represents the saturation of fiber tows. The mesh is changed as shown in Fig. 5. The only current limitation is that one must use bar (1D) elements to describe flow within fiber tows. However, it should be possible to add two-dimensional ‘‘slave’’ elements along the element edges, though this will make the evaluation of ‘‘equivalent’’ properties more difficult. The model described earlier represents a sink term with magnitude independent of fiber tow saturation (below s=1) and linearly dependent on pressure in that region. The sink vanishes once the fiber tow is saturated (s reaches 1). The stepwise approximation may cause some numerical problems, but in our case they turn out to be only in the form of minimal pressure oscillations. This model might be sufficient for approximate modeling of the saturation effects. If a more complex behavior is desired, it may be approximated using a more complex ‘‘slave’’ structure consisting of a chain of bar elements instead of a single ‘‘slave’’ node and a single element. However, it might prove difficult to build such a model and provide adequately accurate parameters for all of the newly created ‘‘slave’’ elements.
3. Implementation We will implement the simple model shown in Fig. 5 and obtain the proper ‘‘equivalent’’ parameters for the ‘‘slave’’ elements used by the model. A few details need to be addressed for this implementation. They are: Fig. 4. Addition of fiber tow internal volume externally by attachment of one-dimensional ‘‘slave’’ elements to all the nodes of an original ‘‘master’’ element (quad in this case).
The resin flow in regions where the fiber tows are saturated must be the same as when saturation effects are not considered in order to use the
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
1729
Fig. 5. Changes in mesh to model saturation of dual-scale porous medium.
commonly available ‘‘saturated’’ permeability value. This can be easily accomplished by maintaining dimensions and permeability of the original mesh intact. The dimensions and permeability of the bar elements should reflect the estimated fiber tow dimensions and permeability. The porous volume in ‘‘slave’’ control volume must correspond to the porous volume in fiber tows. The total porous volume must be conserved. The ‘‘slave’’ node/control volume adds some free volume and even the size of ‘‘master’’ control volume is influenced by the presence of a bar element. Thus, we can change the fiber volume fraction in ‘‘master’’ control volume, its dimensions or both. In order to maintain the permeability value on the level of ‘‘saturated’’ permeability we change the fiber volume fractions in connected elements as described in the next section.
sonable accuracy. Next, we will assume we have some estimate of transverse fiber tow permeability KTtt. This value is usually obtained by measuring the transverse permeability of an aligned fiber bed. However, as this value is difficult to measure accurately, we can lump it with other constitutive parameters we will use in our predictions. Later, it will be shown that it is possible to reduce the number of unknown parameters to one, and this parameter can be found by matching the global experimental results with the simulated ones. We assume the fiber cross-section is a long rectangle with rounded corners and dimensions approximately 2.b 2.h as shown in Fig. 6. This matches the shape in well-compacted performs as seen from Fig. 7. We also assume that the fiber tow saturates by transverse flow from the resin surrounding the tows. For such situations, the resin will largely infiltrate the tows from the upper and top sides and the flow rate within the depicted volume (Fig. 6) can be approximately described by Darcy’s law as: Qin ffi
2:KTtt :2b:L1 :p :h
4. Calculation of input parameters The geometry of the original part (and mesh) is usually known. This may be specified as a three-dimensional part or as a shell, i.e. spatial surface with specified thickness. Also, the fiber volume fraction vf within the mold is known and some estimate of ‘‘saturated’’ preform permeability is available. These data are the usual inputs for any simulation. In order to characterize the medium as a dual-scale medium we need to know additional data concerning the fiber tows. We may safely assume that the fiber volume fraction within tows vfT is known within rea-
Fig. 6. Assumed fiber tow cross-section.
ð7Þ
1730
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
cross-section A, permeability K and fiber volume fraction vf. They are as follows: K ¼ KTtt vf ¼ vfT l ¼ 2:h A ¼ 4:b:Li
Fig. 7. Fiber tow cross sections in compacted preforms.
The length L1 is related to the unit fiber tow volume fraction vT (space occupied by fiber tows in the unit volume, see Fig. 6) as follows: vT ¼
4:b:h:L1 Volume
Ai ¼
The sink term should disappear when fiber tow is saturated, i.e. the complete porous volume within a fiber tow is filled. The porous volume within all fiber tows in a unit volume of material, vp, is equal to vP ¼ vT : 1 vfT ð10Þ This value provides the volume for ‘‘slave’’ finite elements/ control volumes that the resin should fill for these tows to ‘‘saturate’’.
5. ‘‘Effective’’ properties for modified mesh In this section we will describe the transition from original mesh to augmented mesh for saturation evaluation as shown in Fig. 4. The finite elements mesh (elements and nodes) is used to evaluate the pressure field. The rest of the computations are performed on control volume mesh with control volumes centered around individual nodes, i.e. each control volume consists of parts of several elements. The calculation of equivalent properties is based on element data. For each element, its volume Vi and its fiber volume fraction vfi are known. The volume of fiber tows within the element may be computed [see Eq. (9)] as vfi VTi ¼ Vi ð11Þ vfT If we can determine the fiber tow dimensions b and h (Fig. 6), we can compute the length of fiber tows within element Li as Li ¼
ViT 4:b:h
However, we must connect ‘‘equivalent’’ elements as contribution to nodes. To do this, the area A must be divided between the nodes. Easiest way is to divide it evenly:
ð8Þ
The unit fiber tow volume fraction vT can be also evaluated from the preform fiber volume fraction vf and the fiber volume fraction within fiber tows vfT as: vf vT ¼ ð9Þ vfT
ð12Þ
Now we can easily determine all the properties of ‘‘equivalent’’ 1D element for this element: length l,
ð13Þ
4:b:Li NN
ð14Þ
where NN is number of nodes in the element. If mesh is irregular, more complex distribution schemes, consistent with volume integration used in the program, should be used. As mentioned earlier, all the individual element contributions in nodes may be simply merged by addition of the cross-sections as long as the perform characteristics vft, KttT and H remain constant. Otherwise we may either leave several 1D elements connected to each node or create another equivalent element by considering all the individual contributions in parallel to satisfy mass balance. Since there are four parameters to adjust, this will not be a unique process. The last element parameter to adjust is its own fiber volume fraction vfi to be changed into vfie. This is needed to conserve total (macro- and micro-) porous volume. Alternatively, one can modify the element dimensions, but that would be more demanding approach. The conservation of porous volume can be expressed by the equation: Vi : 1 vfi ¼ Vi : 1 vfie þ 8:h:b:Li : 1 vfT ð15Þ Substituting from Eqs. (12) and (13) we obtain vfi Vi : 1 vfi ¼ Vi : 1 vfie þ 2:Vi : 1 vfT vfT
ð16Þ
Thus, vfie
1 vfT ¼ vfi : 1 þ 2: vfT
ð17Þ
Assuming vfT is 0.8, we cannot apply Eq. (12) if vfi is over 66% since Eq. (14) would result in vfie over 1. However, one usually does not have vfi over 66% (or vfT over 0.8) due to packing limitations. Hence, this will not cause any numerical artifact. 5.1. Fiber tow parameters The physical parameters of fiber tows cannot be determined accurately. In reality, the only accurately known value is the volume fraction vfT, which is usually near the maximum packing fraction of around 0.8. The
1731
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
dimensions b, h and the transverse permeability KttT are only known within certain ranges or order of magnitude and may show significant deviations. The significance of this uncertainty can be reduced by lumping these effects into a single fitting parameter if we are interested only in flow rate into the sink and its volume. The volume fractions are well known, so the volumes are likely to be estimated accurately. This provides us with a good approximation of the factor b.h.L1 [see Eq. (6)]. Substituting this into Eq. (5) we can obtain the following equation for flow: Qin ffi
KTtt vT :Volume :p : h2
ð18Þ
If we use the equivalent parameters for the ‘‘slave’’ 1D elements as given in Eq. (13) that we use directly in the simulation, Eq. (18) reduces to Qin ffi
K 4:vT :Volume :p : ð2hÞ2
ð19Þ
This means that in order to provide ‘‘equivalent’’ simulation data we have a single ‘‘uncertain’’ parameter K/(2h)2, which one can be estimated by matching onedimensional pressure profile during an experiment with one-dimensional flow simulation. This parameter will ‘‘lump’’ all the non-uniformities in fiber tow cross-section and variations in flow rates into fiber tows.
6. Examples In this section we will present two examples of filling simulation that include saturation effects. The first example is a simple one: constant flow linear injection. It is, nonetheless, a very important problem because it verifies experimental findings and may be used to characterize the only parameter needed, K/(2h)2. We will also examine the injection pressure and the extent of unsaturated regions within the part. The second example simulates the filling and saturation of fiber tows in a box. This demonstrates not only the ability to predict fill time, but also the ability to predict the time to reach complete saturation as a function of the parameter K/ (2h)2. This will allow us to estimate the duration to bleed resin out of the vent.
If there is no saturation effect, the pressure at the injection gate should build up linearly. In experiments with fiber preforms, the inlet pressure is usually nonlinear to some extent. Our interest essentially is to show that we can duplicate the non-linear pressure development as recorded [8], at least in the qualitative sense when the saturation effects are included. We will use the K/(2h)2 parameter to adjust the nonlinear behavior. The other parameters for the simulation are listed in Table 1 and are close to parameters used to conduct one-dimensional permeability experiments. The value of K/(2h)2 parameter was varied from 1.10 6 to 1.10 12. This would correspond (as h is usually of the order of 1.10 3 m) to the fiber tow permeability range 1.10 12 to 1.10 18 m2. The realistic value should lie somewhere among the smaller values in this range. The inlet pressure dependence on the time is shown in Fig. 8 for a variety of K/(2h)2 values. This set of result shows how the pressure behavior depends on the permeability of the fiber tows. If this permeability is very large, the fiber tows are being filled almost immediately. This is the case for K/ (2h)2=1.10 6, which is unrealistically high value. The small pressure oscillations are an artificial result of abrupt closing of individual ‘‘slave’’ control volumes. The flow front progress is slower as the tows behind the flow-front are immediately saturated and the pressure grows linearly. The slope of the curve is determined by the ‘‘saturated’’ permeability. On the opposite end of spectrum is extremely low value of fiber tow permeability, for example K/ (2h)2=1.10 12. Pressure develops linearly again, since no fiber tows get impregnated in the timeframe of the filling process. However, pressure grows faster because the flow-front moves more rapidly. The change of slope of this line relative to the previous one will depend on preform fiber volume fraction vf (and the porosity of fiber tows vfT, which does not vary). If there are no fiber tows (vf=0), the lines should coincide. ‘‘Unsaturated’’ permeability—if one defines it as a value of permeability when fiber tows do not get saturated at all—could be determined from this curve. Note, however, that the usual concept of unsaturated permeability is rather abstract and it is aimed at reproducing the typical flow behavior rather than describing the flow under particular (extreme) conditions.
6.1. Linear injection We will examine a simple linear injection. Here, the resin is injected with constant flow rate into a narrow mold containing a stationary, dual-scale preform. This is a case of limited interest for manufacturing but is very significant for experimental work, as it will allow one to determine permeability and the parameter K/(2h)2. The injection is discontinued when the macro-flow reaches the end of the mold.
Table 1 Constant flow linear injection parameters Length
1.5 [m]
Width Preform permeability (Kpref) Thickness Fiber volume fraction Fiber tow volume fraction
0.2 [m] 1.10 10 [m2] 0.02[m] 0.5 0.80
1732
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
Fig. 8. Pressure vs. time for constant flow rate injection as predicted by simulation for various values of K/h2.
The more interesting result lies in-between these extremes and covers more realistic estimates of fiber tow permeability. In these cases, the curve follows the upper line, as fiber tows do not saturate immediately. After a while, the fiber tows start to saturate, flow front slows down and the pressure curve converges to the lower line (large fiber tow permeability). Depending on the actual value of fiber tow permeability this might create the impression of smooth, non-linear pressure development as is the case for K/(2h)2=3.10 11 in the figure. The point where the curve joins this line corresponds to the instant when the fiber tows at the inlet are saturated. As the flow into the fiber tows for a given pressure depends linearly on K/(2h)2 [Eq. (18)], one can collapse these curves on each other with proper scaling. The curves obtained correspond well to the experimental data we have [8] in qualitative sense. The attempt to adjust the parameters to match the experimental curves quantitatively is, however, beyond the scope of this paper. The question of fiber tow saturation at the end of the process is also important. This may be addressed by plotting the ‘‘unsaturated length’’ when the resin reaches the end of the mold. We define unsaturated length as the length of the domain behind the flow-front where fiber tow saturation is less than 100%. Fig. 9 shows this concept on the contour plot of saturation value in the sample for Kpref=1.10 10 [m2] and K/ (2h)2=4.10 8. The figure demonstrates that even for realistic input values, there is a small region that is not saturated yet when the macro-flow arrives at the end of the mold. The size of unsaturated length depends on both K/(2h)2 and
Kpref, and it may increase dramatically as the fiber tow permeability decreases. Fig. 10 shows this dependence. The length of unsaturated region, non-dimensionalized by square root of global permeability, is plotted as a function of K/(2h)2. The curves corresponding to different values of global permeability, Kpref should fold on top of each other [27]. They indeed coincide within the limits of accuracy imposed by our computations. The discrepancy in values for the largest fiber tow permeability can be easily explained by interpolation, as the unsaturated region is smaller than an element size in this case. It may be seen that the length of unsaturated region might be significant within the estimated range of values. In practice, the injection is stopped when the resin reaches the end of the mold. Significant part of the preform may be still left unsaturated at this point. After curing, the unsaturated region will have higher void content. This reduces the part quality and may even render the part useless. This result suggests the need for additional studies with saturation modeling. 6.2. Filling of a box mold To demonstrate the applicability of described algorithm for non-trivial simulations, which are generally of significant interest for practical manufacturing, filling of a five-sided box was selected. The geometry we used was from unrelated experimental studies conducted previously. It is a box with missing bottom side shown in Fig. 11. Fig. 11(a) is in the view consistent with the view in which we show the
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
1733
Fig. 9. Contour plot of saturation values for one-dimensional mold filling simulation with K/h2=4.10 8 when resin reaches the end of the mold.
Fig. 10. Fraction of unsaturated length (at the end of the injection) as a function of K/h2.
results and the injection location. Fig. 11(b) is a view from the bottom. Its dimensions are small: 0.1333 0.2032 0.1016 m. The wall thickness is 0.003175 m. The fiber volume fraction is 0.4. Global permeability Kbulk is isotropic. Two values were used, the original value of 1.1833.10 8 m2 and 1.1833.10 9 m2. There are racetracking channels along edges (narrow elements in Fig. 11) simulated by reduction of fiber volume fraction to 0.2 and increase of permeability to 1.6566.10 7 m2 and 1.6566.10 8 m2 respectively for the two cases. The injection is under constant pressure of 100 kPa from the location shown in Fig. 11. Overall, this is a reasonably complex part. Fig. 12 shows a snapshot of saturation during the filling of the mold containing the dual scale preform. The parameter values used are given in the figure. The global permeability is 1.18.10 8 m2, the fiber tow para-
Fig. 11. Geometry of the box mold: (a) view from top (b) view from bottom.
1734
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736 Table 2 Time to fill and time to saturate depending on Ktow/h2 and Kbulk for the box mold Ktow/h2
3
1.10 1.10 4 1.10 5 1.10 6 1.10 7 1.10 8 1.10 9
Fig. 12. Fiber tow saturation during filling of a box structure.
Kbulk=1.18.10 8 m2
Kbulk=1.18.10 9 m2
tfill [s]
tsaturate [s]
tfill [s]
tsaturate [s]
4.665 4.664 4.65 4.54 4.114 3.954 3.938
4.667 4.67 4.67 4.952 7.516 33.149 289.662
46.65 46.65 46.64 46.48 45.4 41.14 39.54
46.65 46.68 46.7 46.96 49.52 75.16 331.5
meter Ktow/h2 is 1.10 6. Saturation decreases from 1 (saturated, dark region) to 0 (unsaturated, light color). Macroscopic flow-front is shown in Fig. 12 as a dark line. The line labeled ‘‘Full Saturation’’ is the boundary of the fully saturated region. All tows from the injection
Fig. 13. Unsaturated zone (darker) as a function of fiber tow K/h2—when all the macro pores are filled.
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
gate to that boundary are completely saturated. Obviously, there is a significant lag between the macroand micro-saturated boundary. The most important implication of this example is the necessity to keep resin flowing long enough to saturate micro-pores within fiber tows in the complete part, not just to fill the macro-pores as is usually done. Table 2 shows time to fill the macro pores (tfill) and time to fully saturate the preform (tsaturate) as a function of Ktow/h2 and Kbulk. As can be seen, decreasing the fiber tow permeability can increase the time to saturate the part beyond any boundaries. The associated reduction in time to fill the macro voids is finite and corresponds directly to the percentage of micro-pores in tows relative to their total volume (which is about 16% in this case). All this can be easily explained using the governing equations. In complex parts such as the one presented there is no ‘‘unsaturated length’’ to measure. One may, instead, introduce concept of ‘‘unsaturated area’’ along similar lines. In the case of three-dimensional part this will become an ‘‘unsaturated volume.’’ Fig. 13 shows how the ‘‘unsaturated area’’ in blue (dark) increases as Ktow/ h2 decreases. Fig. 13 shows the mold at the instant when it is ‘‘conventionally filled,’’ i.e. when the resin filled all macro-space between fiber tows. Obviously, one might still obtain a significant area that is not fully saturated. As expected, lower permeability of fiber tows can dramatically influence the ‘‘unsaturated’’ region.
7. Conclusions We have introduced a method to account for the fiber tow saturation for dual scale porous medium when two distinct scales of pores are present. The implementation that accounts for proper conservation of mass balance is demonstrated in a moving boundary environment in which resin impregnates a dual scale fiber preform. It has been shown that the results do agree with experimental findings for 1D linear injection case. This case may be used to characterize the lower scale of pores with respect to the higher scale by matching the experimental results of the injection pressure with the simulation for a constant flow-rate injection. Within the simulation one can vary the properties of the smaller scale (which are difficult to measure) until a suitable agreement is found. This approach can be used to model filling of arbitrary complex domain and one can predict the size and location of the partially saturated domain for such dual scale porous medium. Alternatively, it can predict times that one should allow the resin to bleed through the vent to obtain complete saturation.
1735
Acknowledgements This work was sponsored by ONR under contract number N00014-00-C-0333.
References [1] Bruschke MV, Advani SG. A finite element/control volume approach to mold filling in anisotropic porous media. Polymer Composites 1990;69. [2] Parnas RS, Salem AJ, Sadiq TAK, Wang H-P, Advani SG. The interaction between micro-and macroscopic flow in RTM performs. Composite Structures 1994;27(1-2):93–107. [3] Sadiq T, Parnas R, Advani SG. Experimental investigation of flow in resin transfer molding. In: International SAMPE Electronics Conference, v 24, Simultaneus Convening of the 3rd International SAMPE Metals Conference and the 24th International SAMPE Technical Conference. Part 1 (of 2), 1992, p. 660–74. [4] Slade J, Pillai K, Advani SG. Investigation of unsaturated flow in woven, braided and stitched fiber mats during mold-filling in resin transfer molding. Polymer Composites 2001;22(4):491–505. [5] Ken Han K, Lee CW, Rice BP. Measurements of the permeability of fiber preforms and applications. Composites Science and Technology 2000;60:2435–41. [6] Markicevic B, Papathanasiou TD. The hydraulic permeability of dual porosity fibrous media. Journal of Reinforced Plastics and Composites 2001;20:871–80. [7] Papathanasiou TD. Flow across structured fiber bundles: a dimensionless correlation. International Journal of Multiphase Flow 2001;27:1451–61. [8] De Parseval Y, Pillai KM, Advani SG. A simple model for the variation of permeability due to partial saturation in dual scale porous media. Transport in Porous Media 1997;27:243–64. [9] Phelan Jr FR, Leung Y, Parnas RS. Modeling of microscale flow in unidirectional fibrous porous media. Journal of Thermoplastic Composite Materials 1994;7(3):208–18. [10] Pillai KM, Advani SG. Numerical simulation of unsaturated flow in woven fiber preforms during the resin transfer molding process. Polymer Composites 1998;19(1):71–80. [11] Hieber CA, Shen SF. A finite element/finite difference simulation of the injection mold filling process. J Non-Newtonian Fluid Mech 1980;7:1–31. [12] Osswald TA, Tucker CL. Polymer Eng and Sci 1988;28:413–20. [13] Kang-Moon-Koo X, Lee-Woo-II X, Yoo-Jung-Yu X, Cho-SungMin XSimulation of mold filling process during resin transfer molding. Journal of Materials Processings and Manufacturing Science 1995;3(3):297–313. [14] Lee LJ, Young WB, Lin RJ. Mold filling and cure modeling of RTM and SRIM processes. Composite Structure 1994;1–2(27):109–20. [15] Phelan FR. Flow simulation of the resin transfer molding process. In: Proceedings of the 7th Technical Conference of the American Society for Composites. Technomic Publ. Co. Inc, 1992, p. 90–102. [16] Chang W, Kikuchi N. Analysis of non-isothermal mold filling process in resin transfer molding (RTM) and structural reaction injection molding (SRIM). Computational Mechanics 1995;1(16): 22–35. [17] Nielsen DR, Pitchumani R. Closed-loop flow control in resin transfer molding using real-time numerical process simulations. Composites Science and Technology 2002;62:283–98. [18] Bickerton S, Stadtfeld HC, Steiner KV, Advani SG. Design and application of actively controlled injection schemes for resintransfer molding. Composites Science and Technology 2001;61: 1625–37.
1736
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
[19] Liu B, Bickerton S, Advani SG. Modeling and simulation of rtm:gate control, venting and dry spot prediction. Composites Part A 1996;27A:135–41. [20] Maier RS, Rohaly TF, Advani SG, Fickie KD. A fast numerical method for isothermal resin transfer mold filling. International Journal of Numerical Methods in Engineering 39 1996. [21] Trochu F, Boudreault JF, Gao DM, Gauvin R. Three-dimensional flow simulations for the resin transfer molding process. Materials and Manufacturing Processes 1995;1(10):21–6. [22] Lewis RW, Usmani AS, Cross JT. Finite element modeling of mold filling. In: Onate E, Periaux J, Samuelsson A, editors. Finite Elements in the 900 s. Barcelona: Springer-Verlag/CIMNE; 1991. [23] Voller VR, Chen YF. Prediction of filling times of porous cavities. International Journal Numerical Methods 10 1995.
[24] Gallez X, Advani SG. Numerical simulations for impregnation of fiber preforms in composite manufacturing: a fast method for three-dimensional geometries. In SPE ANTEC Conference 1997. [25] Ngo ND, Mohan RV, Chung PW, Tamma KK. Recent developments encompassing non-isothermal/isothermal liquid composite molding process modeling/analysis: physically accurate, computationally effective, and affordable simulations and validations. Journal of Thermoplastic Composite Materials 1998;6:493– 532. [26] Gauvin R, Trochu F. Key issues in numerical simulation for liquid composite molding processes. Polymer Composites 1998;3: 233–40. [27] Acheson JA, Simacek P, Advani SG. The implications of fiber compaction and saturation on fully coupled VARTM simulation. Composites Part A (in press).
A numerical model to predict fiber tow saturation during liquid composite molding Pavel Simacek, Suresh G. Advani* Department of Mechanical Engineering and Center for Composite Materials, University of Delaware, Newark, DE 19711, USA Received 22 March 2002; received in revised form 20 June 2002; accepted 25 June 2002
Abstract A dual scale porous medium contains two distinct scales of pores. We will consider the case in which small scale pores are ordered within well-defined sub-regions. Typical examples of such media are textile preforms used in various composite-manufacturing processes. Rigorously, the phenomenon can be modeled by using Darcy’s law to describe flow through porous medium and mass conservation for the flow within the larger pores. The smaller pores can be included within these equations as a sink term. This approach, though straightforward, poses implementation difficulties. In this paper, we suggest an alternative way to model this concept. We use the standard finite element/control volume approach and model the ‘‘internal’’ sink term by appending extra onedimensional elements to control volumes associated with the control volumes of discretized part geometry. This approach offers two advantages over previously attempted schemes: (1) the problem to be solved remains linear and flow can be calculated explicitly within the time domain and (2) existing simulation packages for RTM filling simulation will be able to incorporate saturation effects to simulate the flow in dual scale media. To illustrate this point we present the implementation of the developed algorithm within the frameworks of an existing simulation package LIMS. The implementation is capable of providing saturation data during filling of arbitrarily shaped part and can capture the influence of saturation on filling pressure or flow-rate. This should prove useful in determining the time needed to completely saturate the fiber tows, which is crucial for part performance and also allow us to explain the variation in injection pressure and difficulties in ‘‘assigning’’ single scale permeability for such porous medium. # 2003 Elsevier Ltd. All rights reserved. Keywords: Liquid composite molding
1. Introduction A dual scale porous medium can be defined as a porous medium that contains two distinct scales of pores. We will consider the case in which small scale pores are concentrated within dispersed, well-defined sub-regions encompassed by the larger pores. This is true when the porous medium has an ordered structure. Typical examples of such media are textile preforms used in various composite-manufacturing processes. The fluid flow in such media will saturate the larger pores faster than the smaller pores. To describe the physics of such flows one may still use the traditional governing equations such as the Darcy’s law and continuity equation. The Darcy’s law can be expressed as
* Corresponding author. Tel.: +1-302-831-8975; fax: +1302-8313619. E-mail address: [email protected] (S.G. Advani). 0266-3538/03/$ - see front matter # 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0266-3538(03)00155-6
hv i ¼
K rp
ð1Þ
Here hvi is the volume averaged flow velocity, rp is the pressure gradient, is the viscosity of the fluid and K describes the permeability of the porous medium. If the porous medium was homogenous with one intrinsic permeability, the continuity (mass conservation) equation results in r hv i ¼ 0
ð2Þ
The substitution of Eq. (1) in the continuity Eq. (2), results in the following governing equation [1]: K r rp ¼ 0 ð3Þ However, linear injection experiments with woven preforms have shown that the flow rate and the pressure drop relationship for such preforms is not linear as forecasted by Eq. (1) [2,3]. One must modify the Eq. (2)
1726
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
to account for impregnation of fiber tows when one has such a dual scale porous medium. Many structured (i.e. woven or braided) fibrous preforms have essentially the earlier-mentioned dual-scale geometrical structure [4,5]. This structure consists of the system of fiber tows on the macroscale level and the system of fibers on the microscale level as shown in Fig. 1. The empty pores between the fiber tows are of the order of millimeters and the space between the fibers is of the order of microns. The single scale approach essentially neglects any flow on the microscale, i.e. the details at the fiber tow level. This seems to be reasonable assumption for saturated flows, when the fluid has already filled all the pores [2,6,7]. The situation with unsaturated flows when some of the pores are empty, in the proximity of the moving flow-front during the filling process, is different. As the preform consists of two scales of pores (L1 and L2), the pressure driven flow within the preform is bound to fill the larger pores preferentially. Thus it must be considered as a dual scale process. This notion introduces the accepted terms of ‘‘saturated’’ and ‘‘unsaturated’’ permeability together with the concept of ‘‘partially saturated length’’ [8]. The macroscopic geometry of fibrous preform governs the flow at the ‘‘global’’ or macro- scale. However, on the microscopic level, the fiber tows contain both fibers and empty spaces that are being filled by the resin. The flow into these pores between the fibers continues even after the ‘‘macroscopic’’ flow front has moved beyond them. The partially saturated region is created due to this dual scale nature and is the region where the
macroscopic front has reached and occupied the spaces between the fiber tows but the pores between the fibers in a tow are not yet completely filled. The difference between the idealized flow pattern and the dual scale filling process is shown in Fig. 2. Note that, for the given fluid flow-rate, the macroscopic flow-front may move faster when dual-scale effects are incorporated. This phenomenon is, somewhat artificially, treated by using the ‘‘unsaturated’’ permeability [8]. This behavior is of practical significance. The complete saturation of fiber tows is imperative for part performance. To achieve complete saturation, the injection should continue to bleed out of the vents until the resin occupies the regions not only between the fiber tows but also within the tows. The duration to continue to bleed the resin through the vent is hard to estimate, as little systematic modeling or experimental data is available [9,10]. One can model the flow physics by introducing a subtle change in the governing equation. One can still retain the description of the macro-flow using Darcy’s equation. However, the mass conservation or the continuity equation should account for the loss of resin to the fiber tows at the microscale. This can be done by the incorporation of a sink term q within the continuity equation. Thus Eq. (3) will now assume the form [10]: K r rp ¼ qðp; sÞ ð4Þ
Fig. 1. Dual scale pores in textile preform.
Fig. 2. Idealized and dual scale flow in textile preforms.
The sink term magnitude will depend on the pressure surrounding the fiber tow and the fraction of the pores filled by the resin within a fiber tow. This is usually denoted by a saturation index s, which represents the fraction of empty spaces within a fiber tow occupied by the resin. When one uses Eq. (3) for the situation described by Eq. (4) to measure permeability, it leads to the claims of variable permeability, which are, some-
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
what artificially, solved by the introduction of ‘‘unsaturated’’ permeability within the range of the ‘‘unsaturated length’’. Note that, once the fiber tows within preform are fully saturated (s=1), the sink term q(p,s) will disappear. Then, Eq. (4) is reduced to the original form (3). When the fiber tows saturate at the same rate as the spaces between the fiber tows, then Eq. (2) describes the correct physics. However, if the tows saturate at a slower rate, Eq. (4) must be used. The form describing q(p,s) will vary depending on the perform type and possibly also on some other factors such as the filling history. As an approximation, the fiber tow itself can be considered a Darcian medium with permeability KT and the dependence of flow q on surrounding pressure p will be linear. Dependence of q on the fiber tow saturation s will vary depending on the tow cross-section and the filling history. Accurate description may be very difficult, especially as the properties of individual fiber tows may vary. We will simplify this relation by a step function: q ¼C:p; s < 1 0; s ¼ 1
ð5Þ
Thus, once the tow is saturated, the governing Eq. (4) returns to the traditional form (3). The value of constant C can be characterized by matching global flow rate vs. pressure drop behavior of the preform or from the geometry and dimensions of fiber tows. The dimensions of fiber tows are usually very difficult to measure accurately. It is difficult to include the Eqs. (4) and (5) in existing numerical simulations. Hence, we suggest an alternative approach, which is easier to model, is consistent with the physics and can be easily integrated with numerical models for fluid motion in single scale porous medium.
2. Modeling approach 2.1. FE/CV method As the modeling of impregnation of the fluid into porous medium involves a moving boundary, the finite element/control volume (FE/CV) solution scheme to simulate the filling process has served well to capture this physics efficiently [1]. The solution domain is meshed with a fixed finite element mesh. Control volumes are associated with each mesh node or alternatively every element [Fig. 3 (a)]. Each control volume has a fill factor associated with it. This factor ranges between zero (empty CV) and one (filled CV) and designates how much of the porous volume is already filled with the fluid [Fig. 3 (b)]. The pressure in the empty control volumes is known to be equal to that of the vent and the pressure in the filled control volumes is evaluated by the finite element method, using Eq. (3). Then, the flow Qij between indi-
1727
vidual control volumes i and j is determined using the computed pressure field [Fig. 3 (c)]. ð K ð6Þ Qij ¼ ni : rpdsij sij Once the flow rates are known, flow is advanced by explicit integration in time domain. The time step is selected so as to fill at least one additional control volume. This changes the fluid domain and hence the boundary conditions. The pressure solution is sought for the new domain and this process is repeated until the complete porous medium is saturated. This approach is relatively mature and widespread [1,11–16], fairly simple and quite robust, and its utilization stretches into the fields of optimization and control [17–19]. More importantly, if ‘‘incremental’’ system matrix decomposition is used as outlined in [20], the solution is very fast, as the set of system equations is solved only once for the complete filling cycle. The efficiency comes with some restrictions. The incremental algorithm [20] requires the system matrix not to change during the process, which necessarily means that the process must remain linear. If this is not the case, other efficient solution techniques are available [21–26]. For the global flow modeling, the porous fiber tows can be replaced by an impermeable medium [6,7]. This will provide the proper model for Darcy’s equation. Then, the sink term [Eq. (4)] can be added and removed only after it has consumed the proper amount of resin to saturate fiber tows within the appropriate volume. To accomplish this objective rigorously, one would need to modify the governing equation to Eq. (4). This involves change in whatever simulation package is used. This change may be unmaking of the FEM/CV algorithm described earlier, as the sink term depends on the unknown pressure and changes with fiber tow saturation (and time). Depending on the actual form of q(p,s) it may create non-linearity in the set of governing equations. Moreover, even when the form is carefully selected to circumvent this, the invariance of system matrix with time is violated and the very efficient algorithm [20] cannot be applied. We introduce an alternative approach. Instead of providing the sink term in the domain internally by modifying the equations, it may be added externally in the form of extra volume associated with each unit of the discretized model geometry, such an finite elements or control volumes. The additional volume represents the porous volume in fiber tows and must be compensated by reducing the original volume. To do so, the porosity might be decreased (fiber volume fraction is increased), so as to make the fiber tows impermeable. The extra volume is appended to nodes or control volumes in the form of 1D (runner) elements, as shown in Fig. 4.
1728
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
Fig. 3. FEM/CV approach to filling simulation: (a) domain discretization (b) flow-front and fill factors (c) flow between control volumes.
This approach is relatively easy to accomplish with CV/FEM model described earlier by following these steps: 1. One-dimensional (bar) elements representing volume within fiber tows are connected to each original ‘‘master’’ node/control volume. This operation creates ‘‘slave’’ nodes and control volumes. This process is shown in Fig. 4. 2. The porous volume within fiber tows has to be extracted from the original elements. Thus, either the volume or the fiber volume fraction within the original mesh is modified. Increasing the fiber volume fraction seems to be the easiest possibility, especially for three-dimensional meshes. 3. Simulation is run as usual. The fill factor in ‘‘master’’ control volumes represents fill factor in
macro pores between the fiber tows, while the fill factor in ‘‘slave’’ control volumes, connected by bar elements, represents the saturation of fiber tows. The mesh is changed as shown in Fig. 5. The only current limitation is that one must use bar (1D) elements to describe flow within fiber tows. However, it should be possible to add two-dimensional ‘‘slave’’ elements along the element edges, though this will make the evaluation of ‘‘equivalent’’ properties more difficult. The model described earlier represents a sink term with magnitude independent of fiber tow saturation (below s=1) and linearly dependent on pressure in that region. The sink vanishes once the fiber tow is saturated (s reaches 1). The stepwise approximation may cause some numerical problems, but in our case they turn out to be only in the form of minimal pressure oscillations. This model might be sufficient for approximate modeling of the saturation effects. If a more complex behavior is desired, it may be approximated using a more complex ‘‘slave’’ structure consisting of a chain of bar elements instead of a single ‘‘slave’’ node and a single element. However, it might prove difficult to build such a model and provide adequately accurate parameters for all of the newly created ‘‘slave’’ elements.
3. Implementation We will implement the simple model shown in Fig. 5 and obtain the proper ‘‘equivalent’’ parameters for the ‘‘slave’’ elements used by the model. A few details need to be addressed for this implementation. They are: Fig. 4. Addition of fiber tow internal volume externally by attachment of one-dimensional ‘‘slave’’ elements to all the nodes of an original ‘‘master’’ element (quad in this case).
The resin flow in regions where the fiber tows are saturated must be the same as when saturation effects are not considered in order to use the
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
1729
Fig. 5. Changes in mesh to model saturation of dual-scale porous medium.
commonly available ‘‘saturated’’ permeability value. This can be easily accomplished by maintaining dimensions and permeability of the original mesh intact. The dimensions and permeability of the bar elements should reflect the estimated fiber tow dimensions and permeability. The porous volume in ‘‘slave’’ control volume must correspond to the porous volume in fiber tows. The total porous volume must be conserved. The ‘‘slave’’ node/control volume adds some free volume and even the size of ‘‘master’’ control volume is influenced by the presence of a bar element. Thus, we can change the fiber volume fraction in ‘‘master’’ control volume, its dimensions or both. In order to maintain the permeability value on the level of ‘‘saturated’’ permeability we change the fiber volume fractions in connected elements as described in the next section.
sonable accuracy. Next, we will assume we have some estimate of transverse fiber tow permeability KTtt. This value is usually obtained by measuring the transverse permeability of an aligned fiber bed. However, as this value is difficult to measure accurately, we can lump it with other constitutive parameters we will use in our predictions. Later, it will be shown that it is possible to reduce the number of unknown parameters to one, and this parameter can be found by matching the global experimental results with the simulated ones. We assume the fiber cross-section is a long rectangle with rounded corners and dimensions approximately 2.b 2.h as shown in Fig. 6. This matches the shape in well-compacted performs as seen from Fig. 7. We also assume that the fiber tow saturates by transverse flow from the resin surrounding the tows. For such situations, the resin will largely infiltrate the tows from the upper and top sides and the flow rate within the depicted volume (Fig. 6) can be approximately described by Darcy’s law as: Qin ffi
2:KTtt :2b:L1 :p :h
4. Calculation of input parameters The geometry of the original part (and mesh) is usually known. This may be specified as a three-dimensional part or as a shell, i.e. spatial surface with specified thickness. Also, the fiber volume fraction vf within the mold is known and some estimate of ‘‘saturated’’ preform permeability is available. These data are the usual inputs for any simulation. In order to characterize the medium as a dual-scale medium we need to know additional data concerning the fiber tows. We may safely assume that the fiber volume fraction within tows vfT is known within rea-
Fig. 6. Assumed fiber tow cross-section.
ð7Þ
1730
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
cross-section A, permeability K and fiber volume fraction vf. They are as follows: K ¼ KTtt vf ¼ vfT l ¼ 2:h A ¼ 4:b:Li
Fig. 7. Fiber tow cross sections in compacted preforms.
The length L1 is related to the unit fiber tow volume fraction vT (space occupied by fiber tows in the unit volume, see Fig. 6) as follows: vT ¼
4:b:h:L1 Volume
Ai ¼
The sink term should disappear when fiber tow is saturated, i.e. the complete porous volume within a fiber tow is filled. The porous volume within all fiber tows in a unit volume of material, vp, is equal to vP ¼ vT : 1 vfT ð10Þ This value provides the volume for ‘‘slave’’ finite elements/ control volumes that the resin should fill for these tows to ‘‘saturate’’.
5. ‘‘Effective’’ properties for modified mesh In this section we will describe the transition from original mesh to augmented mesh for saturation evaluation as shown in Fig. 4. The finite elements mesh (elements and nodes) is used to evaluate the pressure field. The rest of the computations are performed on control volume mesh with control volumes centered around individual nodes, i.e. each control volume consists of parts of several elements. The calculation of equivalent properties is based on element data. For each element, its volume Vi and its fiber volume fraction vfi are known. The volume of fiber tows within the element may be computed [see Eq. (9)] as vfi VTi ¼ Vi ð11Þ vfT If we can determine the fiber tow dimensions b and h (Fig. 6), we can compute the length of fiber tows within element Li as Li ¼
ViT 4:b:h
However, we must connect ‘‘equivalent’’ elements as contribution to nodes. To do this, the area A must be divided between the nodes. Easiest way is to divide it evenly:
ð8Þ
The unit fiber tow volume fraction vT can be also evaluated from the preform fiber volume fraction vf and the fiber volume fraction within fiber tows vfT as: vf vT ¼ ð9Þ vfT
ð12Þ
Now we can easily determine all the properties of ‘‘equivalent’’ 1D element for this element: length l,
ð13Þ
4:b:Li NN
ð14Þ
where NN is number of nodes in the element. If mesh is irregular, more complex distribution schemes, consistent with volume integration used in the program, should be used. As mentioned earlier, all the individual element contributions in nodes may be simply merged by addition of the cross-sections as long as the perform characteristics vft, KttT and H remain constant. Otherwise we may either leave several 1D elements connected to each node or create another equivalent element by considering all the individual contributions in parallel to satisfy mass balance. Since there are four parameters to adjust, this will not be a unique process. The last element parameter to adjust is its own fiber volume fraction vfi to be changed into vfie. This is needed to conserve total (macro- and micro-) porous volume. Alternatively, one can modify the element dimensions, but that would be more demanding approach. The conservation of porous volume can be expressed by the equation: Vi : 1 vfi ¼ Vi : 1 vfie þ 8:h:b:Li : 1 vfT ð15Þ Substituting from Eqs. (12) and (13) we obtain vfi Vi : 1 vfi ¼ Vi : 1 vfie þ 2:Vi : 1 vfT vfT
ð16Þ
Thus, vfie
1 vfT ¼ vfi : 1 þ 2: vfT
ð17Þ
Assuming vfT is 0.8, we cannot apply Eq. (12) if vfi is over 66% since Eq. (14) would result in vfie over 1. However, one usually does not have vfi over 66% (or vfT over 0.8) due to packing limitations. Hence, this will not cause any numerical artifact. 5.1. Fiber tow parameters The physical parameters of fiber tows cannot be determined accurately. In reality, the only accurately known value is the volume fraction vfT, which is usually near the maximum packing fraction of around 0.8. The
1731
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
dimensions b, h and the transverse permeability KttT are only known within certain ranges or order of magnitude and may show significant deviations. The significance of this uncertainty can be reduced by lumping these effects into a single fitting parameter if we are interested only in flow rate into the sink and its volume. The volume fractions are well known, so the volumes are likely to be estimated accurately. This provides us with a good approximation of the factor b.h.L1 [see Eq. (6)]. Substituting this into Eq. (5) we can obtain the following equation for flow: Qin ffi
KTtt vT :Volume :p : h2
ð18Þ
If we use the equivalent parameters for the ‘‘slave’’ 1D elements as given in Eq. (13) that we use directly in the simulation, Eq. (18) reduces to Qin ffi
K 4:vT :Volume :p : ð2hÞ2
ð19Þ
This means that in order to provide ‘‘equivalent’’ simulation data we have a single ‘‘uncertain’’ parameter K/(2h)2, which one can be estimated by matching onedimensional pressure profile during an experiment with one-dimensional flow simulation. This parameter will ‘‘lump’’ all the non-uniformities in fiber tow cross-section and variations in flow rates into fiber tows.
6. Examples In this section we will present two examples of filling simulation that include saturation effects. The first example is a simple one: constant flow linear injection. It is, nonetheless, a very important problem because it verifies experimental findings and may be used to characterize the only parameter needed, K/(2h)2. We will also examine the injection pressure and the extent of unsaturated regions within the part. The second example simulates the filling and saturation of fiber tows in a box. This demonstrates not only the ability to predict fill time, but also the ability to predict the time to reach complete saturation as a function of the parameter K/ (2h)2. This will allow us to estimate the duration to bleed resin out of the vent.
If there is no saturation effect, the pressure at the injection gate should build up linearly. In experiments with fiber preforms, the inlet pressure is usually nonlinear to some extent. Our interest essentially is to show that we can duplicate the non-linear pressure development as recorded [8], at least in the qualitative sense when the saturation effects are included. We will use the K/(2h)2 parameter to adjust the nonlinear behavior. The other parameters for the simulation are listed in Table 1 and are close to parameters used to conduct one-dimensional permeability experiments. The value of K/(2h)2 parameter was varied from 1.10 6 to 1.10 12. This would correspond (as h is usually of the order of 1.10 3 m) to the fiber tow permeability range 1.10 12 to 1.10 18 m2. The realistic value should lie somewhere among the smaller values in this range. The inlet pressure dependence on the time is shown in Fig. 8 for a variety of K/(2h)2 values. This set of result shows how the pressure behavior depends on the permeability of the fiber tows. If this permeability is very large, the fiber tows are being filled almost immediately. This is the case for K/ (2h)2=1.10 6, which is unrealistically high value. The small pressure oscillations are an artificial result of abrupt closing of individual ‘‘slave’’ control volumes. The flow front progress is slower as the tows behind the flow-front are immediately saturated and the pressure grows linearly. The slope of the curve is determined by the ‘‘saturated’’ permeability. On the opposite end of spectrum is extremely low value of fiber tow permeability, for example K/ (2h)2=1.10 12. Pressure develops linearly again, since no fiber tows get impregnated in the timeframe of the filling process. However, pressure grows faster because the flow-front moves more rapidly. The change of slope of this line relative to the previous one will depend on preform fiber volume fraction vf (and the porosity of fiber tows vfT, which does not vary). If there are no fiber tows (vf=0), the lines should coincide. ‘‘Unsaturated’’ permeability—if one defines it as a value of permeability when fiber tows do not get saturated at all—could be determined from this curve. Note, however, that the usual concept of unsaturated permeability is rather abstract and it is aimed at reproducing the typical flow behavior rather than describing the flow under particular (extreme) conditions.
6.1. Linear injection We will examine a simple linear injection. Here, the resin is injected with constant flow rate into a narrow mold containing a stationary, dual-scale preform. This is a case of limited interest for manufacturing but is very significant for experimental work, as it will allow one to determine permeability and the parameter K/(2h)2. The injection is discontinued when the macro-flow reaches the end of the mold.
Table 1 Constant flow linear injection parameters Length
1.5 [m]
Width Preform permeability (Kpref) Thickness Fiber volume fraction Fiber tow volume fraction
0.2 [m] 1.10 10 [m2] 0.02[m] 0.5 0.80
1732
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
Fig. 8. Pressure vs. time for constant flow rate injection as predicted by simulation for various values of K/h2.
The more interesting result lies in-between these extremes and covers more realistic estimates of fiber tow permeability. In these cases, the curve follows the upper line, as fiber tows do not saturate immediately. After a while, the fiber tows start to saturate, flow front slows down and the pressure curve converges to the lower line (large fiber tow permeability). Depending on the actual value of fiber tow permeability this might create the impression of smooth, non-linear pressure development as is the case for K/(2h)2=3.10 11 in the figure. The point where the curve joins this line corresponds to the instant when the fiber tows at the inlet are saturated. As the flow into the fiber tows for a given pressure depends linearly on K/(2h)2 [Eq. (18)], one can collapse these curves on each other with proper scaling. The curves obtained correspond well to the experimental data we have [8] in qualitative sense. The attempt to adjust the parameters to match the experimental curves quantitatively is, however, beyond the scope of this paper. The question of fiber tow saturation at the end of the process is also important. This may be addressed by plotting the ‘‘unsaturated length’’ when the resin reaches the end of the mold. We define unsaturated length as the length of the domain behind the flow-front where fiber tow saturation is less than 100%. Fig. 9 shows this concept on the contour plot of saturation value in the sample for Kpref=1.10 10 [m2] and K/ (2h)2=4.10 8. The figure demonstrates that even for realistic input values, there is a small region that is not saturated yet when the macro-flow arrives at the end of the mold. The size of unsaturated length depends on both K/(2h)2 and
Kpref, and it may increase dramatically as the fiber tow permeability decreases. Fig. 10 shows this dependence. The length of unsaturated region, non-dimensionalized by square root of global permeability, is plotted as a function of K/(2h)2. The curves corresponding to different values of global permeability, Kpref should fold on top of each other [27]. They indeed coincide within the limits of accuracy imposed by our computations. The discrepancy in values for the largest fiber tow permeability can be easily explained by interpolation, as the unsaturated region is smaller than an element size in this case. It may be seen that the length of unsaturated region might be significant within the estimated range of values. In practice, the injection is stopped when the resin reaches the end of the mold. Significant part of the preform may be still left unsaturated at this point. After curing, the unsaturated region will have higher void content. This reduces the part quality and may even render the part useless. This result suggests the need for additional studies with saturation modeling. 6.2. Filling of a box mold To demonstrate the applicability of described algorithm for non-trivial simulations, which are generally of significant interest for practical manufacturing, filling of a five-sided box was selected. The geometry we used was from unrelated experimental studies conducted previously. It is a box with missing bottom side shown in Fig. 11. Fig. 11(a) is in the view consistent with the view in which we show the
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
1733
Fig. 9. Contour plot of saturation values for one-dimensional mold filling simulation with K/h2=4.10 8 when resin reaches the end of the mold.
Fig. 10. Fraction of unsaturated length (at the end of the injection) as a function of K/h2.
results and the injection location. Fig. 11(b) is a view from the bottom. Its dimensions are small: 0.1333 0.2032 0.1016 m. The wall thickness is 0.003175 m. The fiber volume fraction is 0.4. Global permeability Kbulk is isotropic. Two values were used, the original value of 1.1833.10 8 m2 and 1.1833.10 9 m2. There are racetracking channels along edges (narrow elements in Fig. 11) simulated by reduction of fiber volume fraction to 0.2 and increase of permeability to 1.6566.10 7 m2 and 1.6566.10 8 m2 respectively for the two cases. The injection is under constant pressure of 100 kPa from the location shown in Fig. 11. Overall, this is a reasonably complex part. Fig. 12 shows a snapshot of saturation during the filling of the mold containing the dual scale preform. The parameter values used are given in the figure. The global permeability is 1.18.10 8 m2, the fiber tow para-
Fig. 11. Geometry of the box mold: (a) view from top (b) view from bottom.
1734
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736 Table 2 Time to fill and time to saturate depending on Ktow/h2 and Kbulk for the box mold Ktow/h2
3
1.10 1.10 4 1.10 5 1.10 6 1.10 7 1.10 8 1.10 9
Fig. 12. Fiber tow saturation during filling of a box structure.
Kbulk=1.18.10 8 m2
Kbulk=1.18.10 9 m2
tfill [s]
tsaturate [s]
tfill [s]
tsaturate [s]
4.665 4.664 4.65 4.54 4.114 3.954 3.938
4.667 4.67 4.67 4.952 7.516 33.149 289.662
46.65 46.65 46.64 46.48 45.4 41.14 39.54
46.65 46.68 46.7 46.96 49.52 75.16 331.5
meter Ktow/h2 is 1.10 6. Saturation decreases from 1 (saturated, dark region) to 0 (unsaturated, light color). Macroscopic flow-front is shown in Fig. 12 as a dark line. The line labeled ‘‘Full Saturation’’ is the boundary of the fully saturated region. All tows from the injection
Fig. 13. Unsaturated zone (darker) as a function of fiber tow K/h2—when all the macro pores are filled.
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
gate to that boundary are completely saturated. Obviously, there is a significant lag between the macroand micro-saturated boundary. The most important implication of this example is the necessity to keep resin flowing long enough to saturate micro-pores within fiber tows in the complete part, not just to fill the macro-pores as is usually done. Table 2 shows time to fill the macro pores (tfill) and time to fully saturate the preform (tsaturate) as a function of Ktow/h2 and Kbulk. As can be seen, decreasing the fiber tow permeability can increase the time to saturate the part beyond any boundaries. The associated reduction in time to fill the macro voids is finite and corresponds directly to the percentage of micro-pores in tows relative to their total volume (which is about 16% in this case). All this can be easily explained using the governing equations. In complex parts such as the one presented there is no ‘‘unsaturated length’’ to measure. One may, instead, introduce concept of ‘‘unsaturated area’’ along similar lines. In the case of three-dimensional part this will become an ‘‘unsaturated volume.’’ Fig. 13 shows how the ‘‘unsaturated area’’ in blue (dark) increases as Ktow/ h2 decreases. Fig. 13 shows the mold at the instant when it is ‘‘conventionally filled,’’ i.e. when the resin filled all macro-space between fiber tows. Obviously, one might still obtain a significant area that is not fully saturated. As expected, lower permeability of fiber tows can dramatically influence the ‘‘unsaturated’’ region.
7. Conclusions We have introduced a method to account for the fiber tow saturation for dual scale porous medium when two distinct scales of pores are present. The implementation that accounts for proper conservation of mass balance is demonstrated in a moving boundary environment in which resin impregnates a dual scale fiber preform. It has been shown that the results do agree with experimental findings for 1D linear injection case. This case may be used to characterize the lower scale of pores with respect to the higher scale by matching the experimental results of the injection pressure with the simulation for a constant flow-rate injection. Within the simulation one can vary the properties of the smaller scale (which are difficult to measure) until a suitable agreement is found. This approach can be used to model filling of arbitrary complex domain and one can predict the size and location of the partially saturated domain for such dual scale porous medium. Alternatively, it can predict times that one should allow the resin to bleed through the vent to obtain complete saturation.
1735
Acknowledgements This work was sponsored by ONR under contract number N00014-00-C-0333.
References [1] Bruschke MV, Advani SG. A finite element/control volume approach to mold filling in anisotropic porous media. Polymer Composites 1990;69. [2] Parnas RS, Salem AJ, Sadiq TAK, Wang H-P, Advani SG. The interaction between micro-and macroscopic flow in RTM performs. Composite Structures 1994;27(1-2):93–107. [3] Sadiq T, Parnas R, Advani SG. Experimental investigation of flow in resin transfer molding. In: International SAMPE Electronics Conference, v 24, Simultaneus Convening of the 3rd International SAMPE Metals Conference and the 24th International SAMPE Technical Conference. Part 1 (of 2), 1992, p. 660–74. [4] Slade J, Pillai K, Advani SG. Investigation of unsaturated flow in woven, braided and stitched fiber mats during mold-filling in resin transfer molding. Polymer Composites 2001;22(4):491–505. [5] Ken Han K, Lee CW, Rice BP. Measurements of the permeability of fiber preforms and applications. Composites Science and Technology 2000;60:2435–41. [6] Markicevic B, Papathanasiou TD. The hydraulic permeability of dual porosity fibrous media. Journal of Reinforced Plastics and Composites 2001;20:871–80. [7] Papathanasiou TD. Flow across structured fiber bundles: a dimensionless correlation. International Journal of Multiphase Flow 2001;27:1451–61. [8] De Parseval Y, Pillai KM, Advani SG. A simple model for the variation of permeability due to partial saturation in dual scale porous media. Transport in Porous Media 1997;27:243–64. [9] Phelan Jr FR, Leung Y, Parnas RS. Modeling of microscale flow in unidirectional fibrous porous media. Journal of Thermoplastic Composite Materials 1994;7(3):208–18. [10] Pillai KM, Advani SG. Numerical simulation of unsaturated flow in woven fiber preforms during the resin transfer molding process. Polymer Composites 1998;19(1):71–80. [11] Hieber CA, Shen SF. A finite element/finite difference simulation of the injection mold filling process. J Non-Newtonian Fluid Mech 1980;7:1–31. [12] Osswald TA, Tucker CL. Polymer Eng and Sci 1988;28:413–20. [13] Kang-Moon-Koo X, Lee-Woo-II X, Yoo-Jung-Yu X, Cho-SungMin XSimulation of mold filling process during resin transfer molding. Journal of Materials Processings and Manufacturing Science 1995;3(3):297–313. [14] Lee LJ, Young WB, Lin RJ. Mold filling and cure modeling of RTM and SRIM processes. Composite Structure 1994;1–2(27):109–20. [15] Phelan FR. Flow simulation of the resin transfer molding process. In: Proceedings of the 7th Technical Conference of the American Society for Composites. Technomic Publ. Co. Inc, 1992, p. 90–102. [16] Chang W, Kikuchi N. Analysis of non-isothermal mold filling process in resin transfer molding (RTM) and structural reaction injection molding (SRIM). Computational Mechanics 1995;1(16): 22–35. [17] Nielsen DR, Pitchumani R. Closed-loop flow control in resin transfer molding using real-time numerical process simulations. Composites Science and Technology 2002;62:283–98. [18] Bickerton S, Stadtfeld HC, Steiner KV, Advani SG. Design and application of actively controlled injection schemes for resintransfer molding. Composites Science and Technology 2001;61: 1625–37.
1736
P. Simacek, S.G. Advani / Composites Science and Technology 63 (2003) 1725–1736
[19] Liu B, Bickerton S, Advani SG. Modeling and simulation of rtm:gate control, venting and dry spot prediction. Composites Part A 1996;27A:135–41. [20] Maier RS, Rohaly TF, Advani SG, Fickie KD. A fast numerical method for isothermal resin transfer mold filling. International Journal of Numerical Methods in Engineering 39 1996. [21] Trochu F, Boudreault JF, Gao DM, Gauvin R. Three-dimensional flow simulations for the resin transfer molding process. Materials and Manufacturing Processes 1995;1(10):21–6. [22] Lewis RW, Usmani AS, Cross JT. Finite element modeling of mold filling. In: Onate E, Periaux J, Samuelsson A, editors. Finite Elements in the 900 s. Barcelona: Springer-Verlag/CIMNE; 1991. [23] Voller VR, Chen YF. Prediction of filling times of porous cavities. International Journal Numerical Methods 10 1995.
[24] Gallez X, Advani SG. Numerical simulations for impregnation of fiber preforms in composite manufacturing: a fast method for three-dimensional geometries. In SPE ANTEC Conference 1997. [25] Ngo ND, Mohan RV, Chung PW, Tamma KK. Recent developments encompassing non-isothermal/isothermal liquid composite molding process modeling/analysis: physically accurate, computationally effective, and affordable simulations and validations. Journal of Thermoplastic Composite Materials 1998;6:493– 532. [26] Gauvin R, Trochu F. Key issues in numerical simulation for liquid composite molding processes. Polymer Composites 1998;3: 233–40. [27] Acheson JA, Simacek P, Advani SG. The implications of fiber compaction and saturation on fully coupled VARTM simulation. Composites Part A (in press).