Determining the 3D permeability of fibrous media using the Newton method

Composites: Part B 36 (2005) 609–618 www.elsevier.com/locate/compositesb

Determining the 3D permeability of fibrous media using the Newton method Daniel Z. Turner*, Keith D. Hjelmstad Department of Civil and Environmental Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Received 1 September 2004; accepted 29 January 2005 Available online 5 July 2005

Abstract This work explores the reliability and sensitivity of using the Newton method to determine the elements of a fiber mat’s permeability tensor. The goal of this work is to determine the extent to which experimental errors affect the numerical results and to provide general guidelines on which types of fiber materials are appropriate for this method. The sensitivity analysis reveals that employing the Newton method is a rather robust approach, applicable to a wide range of materials. The sensitivity analysis also reveals robustness with regard to experimental error as well. q 2005 Elsevier Ltd. All rights reserved. Keywords: Fibres; C. Computational modelling; E. Resin transfer molding (RTM); Physical properties

1. Introduction To model the flow of resin within a mold, using a FE approach, the permeability of the fibers must be known. Permeability is a measure of how easily resin flows within the fibers. A symmetric tensor, K, which describes the flow in all directions completely defines the three-dimensional permeability of the fibers. Much like material properties such as Young’s Modulus or Poisson’s Ratio, which govern a materials response to loading, K characterizes the flow of resin within the fibers. To accurately model flow, one must not only obtain numerical values for K, but also a sense of the value’s quality. Several experimental methods have been proposed in the literature to determine the permeability of a fiber preform [1–5,7,12,16,17]. Important differences to denote among them regard the dimensional capability of the technique and the numerical solution procedure. Many of the twodimensional approaches for determining the in-plane permeability [3–5,7] are based on a combination of Darcy’s law with the continuity equation. In most cases, the permeability is determined numerically through a least * Corresponding author. Tel: C1 217 244 4673; fax: C1 217 244 4974. E-mail address: [email protected] (D.Z. Turner).

1359-8368/$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compositesb.2005.01.009

squares fit of data or extruded from a linear relationship. As flow through the thickness of the preform is considered (three-dimensionally), the analysis becomes more complicated because the solution scheme may require an iterative numerical technique. With regard to these methods, the accuracy of the solution is no longer based solely on the quality of the experiment, but also the numerical method used to solve the equations. This work focuses on the solution of the governing equations, which results from conducting a three-dimensional determination of the permeability. Nedanov and Advani [12] have developed an experimental means of determining K which is quite robust with respect to experimental techniques. Their method avoids several problems associated with race tracking and non-uniform fiber volume fractions. The experimentally refined aspects of their method allow us to concentrate on a numerical evaluation. In their method, K is determined by injecting resin into a plate of fibers. A series of pictures are taken at certain times to capture the flow front of resin. By equating the flow of resin in the injection tube to the expanding ellipsoid of resin in the fibers, one obtains a system of non-linear equations to solve for the elements of K. The Newton method is used to solve the system of non-linear equations. The purpose of this work is to establish the usability and reliability of values for the permeability tensor determined

610

D.Z. Turner, K.D. Hjelmstad / Composites: Part B 36 (2005) 609–618

Nomenclature a b b C g g J k kx ky kz K m n P

f p Q r r V v x X xv y Y z Z

scale factor for x component radius of the injection port scale factor for y component constant used for simplification residual of governing equation scale factor for z component least squares fit function general permeability x component of permeability y component of permeability z component of permeability permeability tensor viscosity of the resin normal vector pressure

by using the Newton method. Our intention is to shed light on the sensitivity of the method to experimental errors and to provide a general framework for which the method is effective. As an example, we examine the data used in the paper by Nedanov and Advani and introduce experimental error by using a random number generator with a specified coefficient of variation. Although the numerical results of this work are specific to the Nedanov and Advani’s data, the method is applicable to a wide variety of fiber mats.

Darcy’s law states that the velocity of fluid flow, v is proportional to the viscosity, m and the permeability, K of the material. In three dimensions, Darcy’s law can be expressed as (1)

A liquid injected from a point source into a semipermeable material forms an ellipsoidal flowfront [2]. If a coordinate system with base {n1,n2,n3} is introduced that coincides with the major axes of the ellipse then the permeability matrix K can be expressed in terms of the principal permeabilities as K Z kx n1 5n1 C ky n2 5n2 C kz n3 5n3

(2)

A closed form solution for flow in an anisotropic medium with a point source can be determined. Let k be an equivalent isotropic permeability determined as k Z ðdet KÞ1=3 Z ðkx ky kz Þ1=3

(3)

Through a change of variable, the ellipsoid of flow can be converted to an equivalent sphere of radius r(t) 2

2

2

X CY CZ Z r

2

where r is the radius and X, Y, and Z are the principal dimensions of the sphere. The principal dimensions of the ellipsoid are related to the principal dimensions of the sphere by the following ratios rffiffiffiffiffi rffiffiffiffiffi rffiffiffiffiffi ky kx kz ; zZZ ; y ZY (5) x ZX k k k Thus, by substituting back into the equation of a sphere, the equivalent sphere has the equation x 2 y 2 z2 r2 C C Z kx ky kz k

2. Theory and governing equations

1 v Z K KVP m

porosity of the preform constant pi flow rate radius of the flow front density of the resin volume of voids in resin bulb velocity of the resin resin bulb width along principle axis coordinate along principle axis 1 vector of iterative estimates for resin bulb width along principle axis coordinate along principle axis 2 resin bulb width along principle axis coordinate along principle axis 3

(4)

(6)

The permeabilities can be determined by measuring the dimensions of the ellipse at various times. Each measurement of the flow front position contains at least three points, one along the principal x-axis, one along the principal y-axis, and one along the principal z-axis. Let the radius of the resin bulb be r(t). The volume of the bulb is 2pr3/3. Let f be the volume fraction of voids. Then the volume of voids in the resin bulb VZ2p3f/3. Thus, the rate of change of the volume of the voids is dV dr Z 2pr 2 f dt dt

(7)

Since the velocity is radial, if the resin is assumed to be incompressible, we have by conservation of mass vðrÞ Z

Q 2pr2

(8)

where Q is the rate of flow. By employing Darcy’s law, Eq. (8) implies that K

k dp Q Z m dr 2pr2

(9)

D.Z. Turner, K.D. Hjelmstad / Composites: Part B 36 (2005) 609–618

where P(r) is the pressure. Integrating the equation gives ðr ð0 k Q K 2p dP Z dr (10) 2 m P1 br where we assume that the pressure has the (measured) value P, at radius b, which we will take to be the radius of the inlet tube, and the pressure is zero at the outside radius r of the bulb. Carrying out the indicated operations and simplifying gives QZ

k 2pbrP1 m r Kb

(11)

Since by definition, QZdV/dt, we have from Eq. (7) 2pr 2 f

dr k 2pbr Z P dt m ðr K bÞ 1

(12)

By rearranging Eq. (12), dividing both sides by 2p, and setting CZbP/fm ðr K bÞr dr Z Ck dt

(13)

Integration of Eq. (13) provides the time it takes to fill a hemisphere of radius r 3

611

and Eq. (15), in discrete form becomes gi Z 6Ckti K 2ri3 C 3bri2 K b3

(17)

where g is the residual that we are trying to minimize. Let kxZak, kyZbk, and kzZgk. Eq. (16) now becomes ri2 Z x2i bg C y2i ag C z2i ab

(18)

To reduce the number of variables in Eq. (18), observe from Eq. (3) that abg K 1 Z 0

(19)

Eq. (19) can be rearranged as gZ

1 ab

(20)

which is used to eliminate g from Eq. (18). This choice is made in anticipation that the method will be applied to systems where kx and ky have similar values, but kz is quite different (generally smaller). Eq. (18) becomes ri2 Z x2i

1 1 C y2i C z2i ab a b

(21)

2

r br K Z Ckt C e 3 2

(14)

where e is a constant of integration. At tZ0, rZb so eZ Kb3/6. Thus, the equation governing the evolution of the radius r is 6Ckt K 2r 3 C 3br 2 K b3 Z 0

which provides an explicit equation for ri in terms of the parameters a and b. A least squares fit of several data points is performed on Eq. (17) as follows: Let

(15) Jðk; a; bÞ h

3. Analysis In the experiment, several measurements of the principal dimensions x(t) and y(t) of the elliptical flow front were determined from photos taken at different times tZti, iZ 1,.,N (see Fig. 1). A measurement of the z(t) dimension was also taken at the final time when the fluid reached the bottom of the preform. Let xiZx(ti), yiZy(ti) and ziZz(ti). Eq. (6) can be written at tZti as x2i y2 z2 r2 C i C i Z i kx ky kz k

(16)

N 1X g2 ðk; a; bÞ 2 iZ1 i

(22)

To find the minimum of J the gradient is taken and set to zero as N vJ X vg Z gi i Z 0 vk vk iZ1

(23)

N vJ X vg Z gi i Z 0 va va iZ1

(24)

N vJ X vg Z gi i Z 0 vb vb iZ1

(25)

These three equations constitute three equations in three unknowns, which can be solved using the Newton method [8]. The residual equation becomes f ðxÞ Z VJ Z

N X

gi Vgi Z 0

(26)

iZ1

where xZ{k,a,b}. Given an estimate xi, an improved estimate can be determined as Fig. 1. Experimental setup.

xvC1 Z xv K ½Vf ðxv ÞK1 f ðxv Þ

(27)

612

D.Z. Turner, K.D. Hjelmstad / Composites: Part B 36 (2005) 609–618

where Vf ðxÞ Z

4. Results and verification N X

gi V2 gi C Vgi 5Vgi

(28)

iZ1

The first derivative terms are computed as vgi Z Cti vk

(29)

  vgi vg vri x2i 2 Z K 3ðri K bÞ zi b K vri vai va a2

(30)

  vgi vg vr y2 Z i i K 3ðri K bÞ z2i a K i2 vb vri vb b

(31)

where from Eq. (21) we have   vri 1 x2i 2 Z z bK 2 2ri i va a

(32)

  vri 1 y2 Z z2i a K i2 2ri vb b

(33)

The second derivative terms are computed in similar fashion as  2  2 v2 gi 3 x2 x ZK z2i b K i2 K 6ðri K bÞ i3 (34) 2ri vava a a    v2 g i 3 y2i x2i 2 2 ZK z aK 2 zi b K 2 K 3ðri K bÞðz2i Þ 2ri i vavb a b (35)  2  2 v2 gi 3 x2 y ZK z2i a K i2 K 6ðri K bÞ i3 2ri vbvb b b

(36)

v2 gi Z0 vkva

(37)

v2 gi Z0 vkvb

(38)

v2 gi Z0 vkvk

(39)

k, a, and b are determined by successfully iterating until the residual is with the tolerance set in the program. The above method was implemented in a Matlab 6.1 code, Determine3DK.m, which is provided in [15].

To verify that Determine3DK.m works correctly, a sample set of data was created by back-calculating the x, y, and z points from a chosen solution. The data was then run in the program to compare the program’s converged solution to the chosen solution from which the data was back-calculated. The chosen solution and the backcalculated data are listed in Table 1. The first column contains the arbitrarily chosen parameters (from which the second column was calculated). The converged solution varies only slightly from the chosen solution suggesting that Determine3DK.m functions correctly. The variation between the converged solution and the chosen solution is due to the precision reported for the converged solution. The verification performed above shows that Determine3DK.m functions correctly for an exact data set, but the original intention of the program was to obtain the solution to an experimental data set, which includes experimental error. The method described in Section 3 is based on a least squares fit of several data points (see Eq. (22)). A least squares fit should account for reasonable experimental error and approximate the best solution to the given data. To further verify Determine3DK.m, another set of data was run for which a solution was provided by Nedanov and Advani [12]. Nedanov and Advani’s data and results are listed in Table 2. The results from Determine3DK.m are also listed in Table 2. Determine3DK.m’s converged solution is quite close to the results provided by Nedanov and Advani. The variation between Determine3DK.m’s solution and the results from Nedanov and Advani are listed in the last column of Table 2. Since the program calculates g from a and b directly, the 3.1% variation in g stems from the 1.5% variation in a and the 1.47% variation in b. Because the governing equations and the formulation described in this work are the same as that used by Nedanov and Advani, one would expect the results to be exactly the same, but the results are slightly different. The variation arises from the method used for solving the system of equations. Nedanov and Advani employ a built-in solver from Maple V to compute the elements of the symmetric K tensor, kx, ky, and kz. Also, they do not scale the elements, as in Eq. (18), but solve for the values directly. Because the goal of this work is to explore the use of the Newton method to determine the permeability tensor, a solver, which incorporates the Newton method, was written into Determine3DK.m. It is unknown whether Maple V uses the Newton method to solve non-linear systems of equations.

Table 1 Sample verification data and results Arbitrarily chosen solution

Back-calculated data

a b g k

x y z t

3.55 2.88 0.0978 1.43!10K7

Converged solution from Determine3DK.m [4.32,5.38,6.13,6.73] [3.89,4.85,5.52,6.058] [1.116] [50,100,150,200]

a b g k

3.5516 2.8799 0.0978 1.4283!10K7

D.Z. Turner, K.D. Hjelmstad / Composites: Part B 36 (2005) 609–618

613

Table 2 Experimental verification data and results Data

x y z t

Results

[3.73,4.81,5.6,7.5,8.97] [3.58,4.61,5.54,6.69, 8.06] [1.49] [55.9,116.1,197.9,317.9,453.8]

Variation (%)

Nedanov and Advani

Determine3DK.m

a b g k

a b g k

5. Sensitivity and reliability In the following section, we evaluate the sensitivity and reliability of employing the Newton method in the solution scheme. In a sense, we are also presenting the effects of measurement errors on prediction errors as the measurement errors manifest themselves in the numerical computation. Thus, one must understand the effects of measurement errors in light of their consequences when the Newton method is used. To aid in this evaluation, we examine a number of aspects of sensitivity including sensitivity to initializing the search routine and the amount of error in the experimental measurements. To successfully employ a non-linear solver, one must initialize the search routine with values that are as close as possible to the solution values. If the initial estimates are too inaccurate, the solver may converge to an extraneous solution or not converge at all. Accurate data are also necessary. Although a least squares fit is applied to the experimental data, the data must also be accurate enough to lead to a converged solution. To test the sensitivity of Determine3DK.m several cases were run, varying the input parameters. The first set of cases provides a general range of initial estimates for which the program will converge. The second set of cases reveals the effect of experimental error on the converged solution. The first parametric study involved varying the initial estimates. For each case, one of the initial estimates varied

3.5546 2.8814 0.0976 1.4386!10K7

3.5010 2.8389 0.1006 1.4005!10K7

1.5 1.4 3.1 2.6

from K100 to 100% of the solution to see if Determine3DK.m would converge (and if it did, how many iterations it took). Fig. 2 shows the number of iterations required to reach the converged solution (given in Table 1) to the backcalculated exact data for given variations in the initial guesses. The results of the first parametric study reveal that Determine3DK.m converges for initial estimates within 15–20% of the exact solution. If the initial estimates are worse than 20% greater or less than the solution, the program will not converge (as least within a limit of 15 iterations). The results illustrate an interesting challenge to successfully employ the Newton method, that the experimentalist must have a pretty good idea of what the solution will be in order to determine it. Although observing the results of varying one set of data does not fully describe the program’s sensitivity in general, it provides a picture of the window of allowable inaccuracy. Fig. 2 also reveals another interesting feature about the Newton method; the initial estimate of k is almost inconsequential. Above a 5% variation, no matter how poor the initial estimate of k is, the number of iterations required to solve does not increase. To understand why the initial estimate of k does not affect the number of iterations, one must reexamine Eqs. (17) and (23). Observe that Eq. (17) is linear in k. The linearity in k leads to constant first derivatives and zero second derivatives. The linearity is evident in the convergence of the algorithm. The linearity of k in Eq. (17) suggests that Eq. (23) might be explicitly

Fig. 2. Number of iterations required for variations in initial estimates.

614

D.Z. Turner, K.D. Hjelmstad / Composites: Part B 36 (2005) 609–618

Fig. 3. Sensitivity of the solution set to randomly generated input data with varying degrees of imposed error.

determined (in terms of a and b). Then k could be eliminated by substitution in Eqs. (24) and (25). The second parametric study involves producing randomly generated error in the input data. The data used for this study are that of the exact solution and backcalculated data listed in Table 1. To simulate the effects of varying degrees of experimental error, 1500 or more iterations were run for increasingly varied input data. The variation was accomplished by adding or subtracting random amounts from the data listed in Table 1. The randomly generated error is quantified by the range or amplitude over which the error varies. This amplitude is represented by the coefficient of variation (CV) of the input data. For example, for the first set of 1500 iterations, the amplitude of the random error was set to a CV of 0.058, meaning that the input data listed in Table 1 for x, y, and z were randomly augmented 1500 times to create 1500 data sets for which the overall cluster of points for x, y, or z were within a CV of 0.058 from the data listed in Table 1. The elements x, y, and z in the data set were randomly varied simultaneously. The experimental error for the x, y and z data listed in Table 1 is considered zero because the data was back-calculated from an arbitrarily chosen solution set. These values serve as the mean from which the standard deviation and thus the CV are calculated. The goal of the second parametric study is to determine the sensitivity of the method employed in Determine3DK.m to experimental error. The desired sensitivity is represented by the CV of the solution for a given CV of input data. A large CV for the solution value resulting from a small CV in the input data (which represents the experimental error) would suggest that the method is highly sensitive to experimental error. The converse is also true. Fig. 3 shows the CV of the solution values for a given CV of input data. One can see that the relationship between the CV of the solution value to the CV of the input data is almost linear.

Also the CV for the solution is greater than the CV of the input data for any given point, but not by much. The two aforementioned points suggest that the method is quite robust regarding its sensitivity to experimental error. If the method were highly sensitive, the CV of the solution would increase rapidly for an increased CV of the input data. The key question regarding sensitivity concerns the response of the Determine3DK.m to error in the z measurement. The experimental procedure provides only one measurement for z (rather than the N values for x and y). Do errors in the z measurement lead to more substantial deviations in the solution values? Similar iterations run varying only z(1) reveal that the method is no more sensitive to z than it is to errors in the x and y data sets, probably because the time is large. Fig. 4 shows Determine3DK.m’s response to random variations in z only. Fig. 4 reveals that approximately the same relationship holds between the CV of the input values and the CV of the solution. The correlation supports Nedanov and Advani’s postulation that an accurate solution can be achieved with only one point for z and that our similar approach does not decrease the accuracy of our determined parameters. The parametric study described above (with random errors in x, y, and z) reveals not only the sensitivity of the method to experimental error, but also which solution values are the most sensitive to error. Fig. 5 shows three permutations of the solution values plotted against each other. Each point in the graphs represents a given solution for a particular iteration. The cluster of points in each graph represents 5000 iterations run for a CV of 1. The graphs on the left-hand side of Fig. 5 show the points for each converged solution. The graphs on the right-hand side depict the number of points in each grid box of the graphs on the left. To interpret the graphs, one must evaluate the shape of the darker cluster of points. The axis for which the cluster

D.Z. Turner, K.D. Hjelmstad / Composites: Part B 36 (2005) 609–618

615

Fig. 4. Sensitivity of the solution set to randomly generated data for z(1) only (curves for alpha and beta are identical).

extends the longest is the most sensitive direction. One can see from Fig. 5(b) and (c) that the axis for which the cluster extends the longest is the k-axis. Notice that in Fig. 5(b) and (c), the cluster is almost confined to the grid axes. This is due to the nature of the relationship between k and a or b. k is a measure of the general permeability of the material. a and b scale the permeability for a specific direction within the material. k and a or b are directly related to each other through Eqs. (21) and (17), but a and b are independent of each other. a and b scale the radius which was determined by k. Therefore, for high values of a, k is low and vice versa. The same applies to b. The identical nature of Fig. 5(b) and (c) is also due to the independence of both a and b. a and b scale the overall k of the material in the same manner. Therefore, they respond to random variation similarly. Notice also that for Fig. 5(a), the data points are not constrained to the axes. Fig. 5(a) further proves the independent nature of a and b. Although the locus of points is still in close proximity to the exact solution (1,1), they are vastly dispersed. Any given a is not dependent on b, but rather the r data points (which are derived from k). Why must we employ the Newton method if the permeabilities are independent? Although the permeabilities are independent of each other, they are still coupled through the radius of the expanding sphere. In other words, the permeabilities do not compete with each other, but collectively contribute to the general permeability of the material. Fig. 5 shows the following hierarchy of sensitivity, k is the most sensitive followed by a and b, which are equally sensitive. The higher sensitivity of k is shown by the locus of points in Fig. 5(b), for which the longer axis is along the k direction. The similar nature of the sensitivity of a and b is highlighted by Fig. 5(b) and (c), which are almost identical. This hierarchy of sensitivity is also supported by Fig. 4 which shows that for any given CV of input data, the CV of the solution for k

is always greater than that for a or b which are approximately equally sensitive. The preceding discussion reveals that one can know the a and b parameters with greater certainty than k. Since k represents the overall permeability of the fibers and a and b simply scale k to apply to the x or y direction, the experimentalist can know more accurately how to scale k than the value of k itself.

6. Utility To explore the probability of success with this numerical method in extreme cases, we examine two types of fiber materials. Unlike the typical fiber material for which azbzgz1, in material A (see Fig. 6), the permeabilities for directions perpendicular to the injection direction are very different. For example, let us consider materials for which a is several times greater than b. Such is the case for dense materials with all of the fibers running in one direction. As the resin enters the fibers, it has a much greater propensity to flow in one direction over the other. Material A is shown in Fig. 6. When these materials are injected with resin, the expanding resin bulb grows much faster in the direction of the fiber alignment. Even if the fibers are not uniformly aligned, the bulb grows faster in the direction in which the majority of fibers are aligned. In the extreme case, the permeability parallel to the fibers overwhelmingly controls the flow. Determine3DK.m converges for any scale of magnitude between a and b, thus the question of utility lies in the quality of the solution for materials such as material A. To explore the quality of the solution, 5000 samples with random data were run for material A to observe if the scatter plots changed shape. Fig. 7 plots the results from each iteration.

616

D.Z. Turner, K.D. Hjelmstad / Composites: Part B 36 (2005) 609–618

Fig. 5. (a) Solution values for a divided by the exact solution a (x-axis) vs. solution values for b divided by the exact solution b (y-axis). (b) Solution values for b divided by the exact solution b (y-axis) vs. solution values for k divided by the exact solution k (x-axis). (c) Solution values for a divided by the exact solution a (y-axis) vs. solution values for k divided by the exact solution k (x-axis).

Fig. 6. Material A.

D.Z. Turner, K.D. Hjelmstad / Composites: Part B 36 (2005) 609–618

617

Fig. 7. (a) Solution values for a divided by the exact solution a (x-axis) vs. solution values for b divided by the exact solution b (y-axis). (b) Solution values for b divided by the exact solution b (y-axis) vs. solution values for k divided by the exact solution k (x-axis). (c) Solution values for a divided by the exact solution a (y-axis) vs. solution values for k divided by the exact solution k (x-axis).

As seen in Fig. 7, the scatter plots are almost identical to the plots from the typical material. The only difference is the values of magnitude along the axes. The similarity of the shapes proves that the standard deviation of the solution in the extreme case is merely a scaled value of the standard deviation of the typical case. As a or b increases, the standard deviation (relative to the typical case) increases by a factor of 1.3 times the sought parameter. For example, the standard deviation of the typical case for a or b is 1.12, but the standard deviation for the case aZ20, bZ1 is (1.3)(20)Z26. The linear relationship between the standard deviations simply means that as the disparity between the parameters grows, the quality of the solution decreases

rapidly. In terms of utility, Determine3DK.m is functional for all cases, but the reliability of the solution for extreme cases is much less than for typical materials. In material B (see Fig. 8), the permeability parallel to the injection direction is many times greater than any perpendicular direction. The resin flow is characterized by a tubular cylinder that extends from the source for which the width varies slowly, but the height grows rapidly. For this type of flow, the radius is nearly constant. Thus the flow can be represented by a linear formulation for a fluid flowing in a channel with a constant permeability in one direction. Material B is merely an alternate representation of material A for which the injection is made parallel to the fibers.

618

D.Z. Turner, K.D. Hjelmstad / Composites: Part B 36 (2005) 609–618

Fig. 8. Material B.

Interestingly, the program again converges without limit in regard to extreme permeabilities parallel to the injection direction, g, but the standard deviation does not increase as rapidly. Whereas in material A the standard deviation increases by a factor of 1.3, in material B the standard deviation increases by a factor of 0.55. Why does the program converge with more accuracy for material A than for material B? Conceptually, this makes sense because in the case of material A, the resin is injected in a direction in which it cannot flow. Resin flow is almost prohibited for material A. Determine3DK.m is only useful for materials, such as material A, if the permeability is extremely high in all directions perpendicular to the injection, thus not violating our fundamental assumption that the flow bulb expands as an elliptical sphere. In this case, the flow remains bulb-like, even though it is quite shallow. The above discussion provides a guideline for which materials are appropriate for Determine3DK.m. The method is successful for all materials for which the resin flow does not violate our fundamental assumption that the fluid expands as an elliptical sphere. When this principle is violated, the quality of the solution degrades rapidly. Inappropriate materials would be those with dense fibers of uniform orientation perpendicular to the injection direction.

7. Conclusion The preceding discussion shows that the Newton method works quite well for determining the permeability of a fiber

material. Not only is it easy to employ, but it also provides highly reliable results for appropriate materials. The permeability of the fibers is the foundation on which accurate numerical modeling is built. Thus, the high level of confidence in obtaining the problem parameters translates to more confidence in accurately modeling the resin flow. This work also shows that including linear equations in the system of equations solved by the Newton method is appropriate, but not necessary. Linear equations can be solved implicitly and removed from the system, thus reducing the number of unknowns. We have determined the range in which one must initialize the Newton method in order to be successful. The initial estimates of the unknowns must be with 15% of the solutions for the program to converge. Lastly, we have revealed the importance of determining the z data point. Great care should be taken in determining z because of its direct impact on the solution.

References [1] Adams K, Rebenfeld L. Permeability characteristics of multiplayer fiber reinforcements Part I: experimental observations. Polym Compos 1991;2:179. [2] Ahn S, Lee W, Springer G. Measurement of the three-dimensional permeability of fiber preforms using embedded fiber optic sensors. J Compos Mater 1995;29:714–33. [3] Ferland P, Guittard D, Trochu F. Concurrent methods for permeability measurement in resin transfer molding. Polym Compos 1996;17(1): 149–58. [4] Gauvin R, Trochu F, Lemann Y, Diallo L. Permeability measurement and flow simulation through fiber reinforcement. Polym Compos 1996;17(1):34–42. [5] Rikard B, Lidstrom P. Measurement of in-plane permeability of anisotropic fiber reinforcements. Polym Compos 1996;17(1):43–51. [7] Han K, Lee C, Rice B. Measurement of the permeability of fiber preform and applications. Compos Sci Technol 2000;60:2435–41. [8] Hjelmstad K. Fundamentals of structural mechanics. Englewood Cliffs, NJ: Prentice-Hall; 1997 p. 369–72. [12] Nedanov P, Advani S. A method to determine 3D permeability of fibrous reinforcements. J Compos 2002;36:241–54. [15] Turner D. A finite element approach to anlayzing resin flow during injection for the GFRP cuff joint mold. Master’s Thesis. University of Illinois at Urbana-Champaign, Urbana, USA; 2003. [16] Woerdman D, Phelan F, Parnas R. Interpretation of 3-D permeability measurements for RTM manufacturing. Polym Compos 1995;16: 470–80. [17] Wang TJ, Wu CH, Lee JL. In-plane permeability measurement and analysis in liquid composite modling. Polym Compos 1994;15:278–88.