Stochastic resin transfer molding process
Abstract
We consider onedimensional and twodimensional models of the stochastic resin transfer molding process, which are formulated as random moving boundary problems. We study their properties, analytically in the onedimensional case and numerically in the twodimensional case. We show how variability of time to fill depends on correlation lengths and smoothness of a random permeability field.
1 Introduction
Over the last few decades the use of fiberreinforced composite materials in aerospace, automotive and marine industries, in sports and other areas has seen a significant growth. The main benefits of composites include their lightweight and highperformance nature, as well as their flexibility to accommodate complex designs (see e.g. [4, 2, 10] and references therein). One of the main manufacturing processes for producing advanced composites is resin transfer molding (RTM), which belongs to the Liquid Composite Molding class of composite manufacturing processes. RTM has five main stages [4, 2, 23, 7]: (i) manufacturing of a reinforcing preform (e.g., from carbon fiber, glass fiber or other fabric); (ii) packing the preform in a closed mold, which has a cavity with the shape of the designed part; (iii) injecting resin into the mold cavity to fill empty spaces between fibers; (iv) resin curing (it can start during or after the injection stage); and (v) demolding, i.e., taking the solidified part, after completion of curing, from the mold. In this work we will consider the third stage (filling the preform by resin) only and we will, for simplicity, assume that resin curing starts after completion of the filling, or in other words, that the filling process is much faster than the curing, which is a common assumption for RTM. The injecting resin stage is crucial for getting the expected properties of a material.
It is widely accepted (see e.g. [13, 14, 18, 24, 26, 29, 40, 49, 54] and references therein) that composite manufacturing processes are accompanied by uncertainties. The origins of these uncertainties include (a) variability of fiber placements due to imperfections of stages (i)(ii) of RTM; (b) variability of fiber properties; (c) variability of resin properties; and (d) variability of environment during stages (i)(iv) of RTM. Deviations from the design caused by uncertainties can result in defects, which have two manufacturing consequences [2, 24]. First, to avoid compromising performance of the material due to possible defects, one usually uses more conservative designs, making composites thicker and hence more expensive and less lightweight. Second, defects (e.g. dry spots) can lead to a relatively large amount of scrap which increases the cost of the material. Being able to quantify these uncertainties is highly important for further advances of composite manufacturing.
The abovementioned variabilities usually cause uncertainty of permeability/hydraulic conductivity, which in its turn leads to variability in mold filling patterns and filling times in RTM. In this paper our main focus is on dependence of variability of filling times on properties of random fields which model uncertainty of hydraulic conductivity.
Estimation of a filling time for each particular design of a composite part is important for a successful manufacturing process. The filling time should be neither too short nor too long. On the one hand, it should be sufficiently long to allow adequate impregnation of the fibers. On the other hand, a too long filling time can lead to premature gelation (or even solidification) of the resin, which should be avoided as it is a source of defects. Further, a large variability of the filling time is, practically, highly undesirable as it affects robustness of the technological process, making it difficult for automation and standardization (e.g., of an operator’s guidance). Thus, understanding of factors influencing filling time variability is of importance for fiberreinforced composite manufacturing.
We start (Section 2) with considering a stochastic onedimensional movingboundary problem, where analytical analysis is possible. The hydraulic conductivity is modelled as a stationary lognormal random field. In particular, we observe that though the mean filling time does not depend on correlation length of hydraulic conductivity, the filling time variance (as a measure of variability) does depend on the correlation length. We also consider dependence of variation of the filling time on smoothness of the random hydraulic conductivity. In Section 3 we formulate two variants of a twodimensional movingboundary problem with a random hydraulic conductivity tensor. In some previous related works (see e.g. [42, 54]) on stochastic movingboundary problems, it was assumed that the random hydraulic conductivity is isotropic. However, permeability (and hence hydraulic conductivity) is anisotropic by design in most cases of practical interest (see e.g.[8, 37, 30]). Moreover, geometric variability (i.e., deviations from the design such as variation of gaps between fiber yarns, variation of width and angles of fiber yarns in comparison with design, etc.) was observed in experiments (see, e.g. [39, 18, 24, 26] and references therein), which gives further evidence for the need to model hydraulic conductivity as an anisotropic random field. The main difficulty in modeling 2 and 3dimensional resin transfer processes is the possibility of dry spot appearances, i.e., forming enclosed areas with moving fronts behind the main front. The enclosed areas contain air, which is a compressible medium, while resin can be considered incompressible. A full description of this phenomenon requires a twophase model involving both incompressible and compressible phases, which is too computationally demanding. Here we work with a reduced model in which the air entrapment is taken into account by modifying the boundary conditions on the internal fronts. A discretetype version of such a model was proposed in [23] (see also [2]). Here we formulate its continuous analogue. We note that though our study is motivated by technological processes in production of composite materials, the considered models are also useful for other porous media problems (see [8, 37, 42] and references therein). The twodimensional model is solved numerically using the control volumefinite element method (CV/FEM) which we present for completeness in Section 4. The results of our numerical study of the twodimensional stochastic movingboundary problem are given in Section 5, where we also experimentally examine convergence of the CV/FEM algorithm with the quantities of interest being the filling time and void content. A discussion and concluding remarks are in Section 6.
2 Onedimensional model
We start with studying a onedimensional movingboundary model for a stochastic RTM process. In what follows we assume that we have a sufficiently rich probability space Let be a random hydraulic conductivity defined on and taking values on the positive real semiline. Assuming the resin’s incompressibility, the moving boundary problem for the pressure of resin takes the form [2, 42]:
(2.1)  
where is the moving boundary, is the medium porosity, is a pressure on the inlet and is the pressure at the outlet .
We recall that the hydraulic conductivity can be expressed as
(2.2) 
where is the permeability of the medium and is the viscosity of the resin.
Remark
For simplicity we assume that the porosity in (2.1) is constant. Porosity can be assumed constant when its variability is significantly smaller than variability of hydraulic conductivity, which is often the case (see e.g. [8, 37]). It is possible to modify the arguments of this section for the case of random porosity but we do not consider it here. Also, for definiteness, we impose the constant pressure condition at the inlet but other boundary conditions (e.g. constant rate) can be also considered.
Remark
We note that without losing generality we can put in (2.1) (i.e., either it is vacuum on the outlet or is considered as the relative pressure) but for the sake of the twodimensional model considered in the next section, it is convenient to keep the parameter
Let
Assumption 2.1. We assume that the random field is such that
(i) for and all
(ii) the integral
exists for and a.e.
Proposition 1
Note that Assumption 2.1 implies that is a.s. continuous and strictly increasing, which guarantee existence of the inverse , needed for (2.3). We also remark in passing that if dynamics of the front are given then the nonlinear problem (2.1) becomes linear.
As we mentioned in the Introduction, from the application’s point of view, an important characteristic is the time to fill a piece of material of length i.e., the random time such that
(2.4) 
It follows from (2.3) that
(2.5) 
It is natural to assume that the mean of corresponds to the hydraulic conductivity intended by the design of a composite part and the perturbation of this mean is a stationary random field which models uncertainty due to the manufacturing process. Let us now consider the case of the hydraulic conductivity being a stationary lognormal random field, which is a commonly used assumption (see, e.g. [29, 18, 24, 44, 54]), i.e.,
(2.6) 
where and is a stationary Gaussian random field with zero mean and covariance function Note that
(2.7)  
The filling time according to the design is equal to
The first condition in Assumption 2.1 is obviously satisfied by from (2.6). To satisfy the second condition, it is sufficient to require that realizations of are continuous with probability one and we make the following assumption.
Assumption 2.2. We assume that the stationary Gaussian random field has zero mean and continuous covariance function such that for some and
(2.8) 
for all with
Under Assumption 2.2 the random field has continuous sample paths with probability one (see e.g. [1]).
For from (2.6), we obtain the following statistical characteristics of the filling time:
(2.9)  
To understand the behavior of and in terms of it is convenient to rewrite (2.9) via the mean and variance of :
(2.10)  
(2.11) 
It is interesting that the mean filling time depends only on Var and and does not depend on the covariance, and hence it does not depend on the correlation length and smoothness of the random field We pay attention to the interesting fact that grows linearly with Var
The variance of the filling time depends on covariance of . Then to understand its behavior with respect to correlation length and smoothness of let us consider some particular cases of
First let us look at the simple case of being independent of i.e., being just a Gaussian random variable with zero mean and variance which can be viewed as a perfectly correlated medium. Then
(2.12) 
i.e., grows cubically with increase of Var in this case.
Remark
In this case of a perfectly correlated medium we also have from (2.5):
(2.13)  
We note in passing that (2.13) is often used in experiments for estimating an effective macroscopic hydraulic conductivity (and hence effective permeability) via observing time to fill of samples of a material (see e.g. [26]). But it is not difficult to see that when the hydraulic conductivity is not perfectly correlated, (2.13) might not be a good way to estimate the mean and variance of hydraulic conductivity (cf. (2.10)(2.11)).
Now consider the following Matérn covariance function for the stationary Gaussian random field (see [25, 51, 41, 35]):
(2.14) 
where is the gamma function, is the modified Bessel function of the second kind, is variance, is a smoothness parameter, is a characteristic correlation length, and is a scaled distance function. In particular, we have
(2.15)  
(2.16)  
(2.17)  
(2.18) 
The parameter in the above covariance functions controls the degree of smoothness of sample paths of the random field. The random field with Matérn covariance function (2.14) has sample path continuous derivatives with probability one. Hence, with the exponential covariance function from (2.15), the random field has continuous (but not differentiable) sample paths with probability one; with from (2.16) – once differentiable sample paths with probability one; with from (2.17) – twice differentiable sample paths with probability one; and with from (2.18) – infinitely many times differentiable sample paths.
Note that we will use the four Matérn covariance functions in this section for being onedimensional and in the next sections for being twodimensional. The scaled distance of the form , with being the usual Euclidean distance, corresponds to an isotropic random field. In the twodimensional case considered later in the paper we also model the random conductivity as an anisotropic random field, appropriately choosing the scaled distance (see Section 5).
We recall (cf. (2.9)) that the mean filling time does not depend on the choice of covariance and hence, in particular, it is the same for all four Matérn covariance functions (2.15)(2.18). It is not difficult to show that, in the case of these Matérn covariance functions, variance of the filling time has the following properties: (i) it is increasing with growth of the correlation length ; (ii) for small correlation lengths relative to the size of the material grows linearly with (iii) for large correlation lengths is almost independent of (note that in the case of a perfectly correlated medium we had (2.12)); and (iv) for a fixed variance of the filling time grows with smoothness of the random field. We illustrate these properties in Fig. 2.1.
As an example, we also give the expansion of in the case of and small (expansions in small of statistical moments of the interface dynamics and of the pressure were derived in [22]):
where
To summarize, in the onedimensional case,

The mean filling time does not depend on the correlation length or on smoothness of the random hydraulic conductivity . It decreases with increase of the mean of the hydraulic conductivity and with increase of ; and it linearly increases with increase of variance of the hydraulic conductivity and quadratically with increase of the length .

The standard deviation of filling time grows as for small correlation lengths and is almost independent of for large correlation lengths . For a fixed it also grows with increase of smoothness of the random field.
The important consequence of these observations is significance of the correlation length for variability of the filling time. Dependence of variability of on smoothness of serves as a warning that stochastic modeling of permeability via homogenization procedures needs to be done carefully. We also note that the mean filling time is larger than the filling time expected from the design. Further, standard deviation of the filling time is comparable with the mean filling time, i.e. variability of the filling time is high, which can cause problems in fiberreinforced composite manufacturing as explained in the Introduction.
3 Twodimensional model
In this section we formulate the twodimensional analog of the onedimensional model (2.1). To represent a mold, consider an open twodimensional domain with the boundary , where is the inlet, is the perfectly sealed boundary, and is the outlet. Let be a random secondorder hydraulic conductivity tensor defined on . The movingboundary problem for the pressure of resin takes the form (cf. [2, 42]):
(3.1)  
where is the velocity of the moving boundary in the normal direction and and are the unit outward normals to the corresponding boundaries, is the timedependent domain bounded by the moving boundary and the appropriate parts of , is the medium porosity, is a pressure at the inlet and is the pressure at the outlet . Remark 2 is applicable here.
Behavior of the twodimensional model (3.1) is considerably more complex than of the onedimensional model (2.1). Let us start with an illustrative example (see also e.g. [2, Ch. 8]). Consider a rectangular piece of material which has hydraulic conductivity constant everywhere in the domain (here is the unit matrix) except a relatively small region which has (again constant) hydraulic conductivity , where (see Fig. 3.1). In this case the resin can race around the low permeability region and the front becomes discontinuous, creating a macroscopic void behind the main front as demonstrated in Fig. 3.2. Based on the onedimensional model (2.1) and Proposition 1, it is not difficult to estimate that it is sufficient to have for appearance of a void. Macroscopic voids are one of the main defects in composites leading to scrap and failures (see e.g. [2] and references therein). Possible discontinuities of the front cause difficulties in both analytical analysis of (3.1) and its numerical approximation as we discuss further in Sections 4 and 5.1.
In practice one does not have a deep vacuum (i.e., cannot be assumed negligible) and air is entrapped in macrovoids [2]. To take into account air entrapment, one can replace the model (3.1) by a twophase model with one phase (resin) being incompressible (as it is in (3.1)) and the other (air) being compressible. But such a model is computationally expensive, while to find an optimal design for a composite part’s production (e.g. optimal locations of vents and inlets), especially taking into account uncertainties, one needs to run model simulations very many times. To this end a simplified model is considered [23, 2], in which the air entrapment is taken into account by modifying the boundary conditions on the internal fronts. Such a model is widely used in the engineering community working on advanced composite manufacturing, including related commercial software (see, e.g. [2] and references therein). In the simplified model it is assumed [23, 2] that the pressure in a void increases according to the ideal gas law, i.e., that pressure inside the void multiplied by its volume remains constant when the void shrinks during the filling process (we assume that temperature remains constant during the process). In [23, 2] there is a discretetype formulation of this model, here we give its PDE formulation.
To describe this modification of (3.1), assume that at time there are entrapments with closed boundaries and volumes formed at times behind the main front Then we can write the model as
(3.2)  
where are the velocities of the moving boundaries in the normal direction and and are the unit normals to the corresponding boundaries, is the timedependent domain bounded by , , and the appropriate parts of , the rest of the notation is as in (3.1). The volume of a void is computed as
where is a fixed thickness of the material, which is assumed to be small so that flow through thickness can be neglected, i.e., that the twodimensional model is a good approximation for the three dimensional flow.
It is clear that in the model (3.2) once a void is formed at , its volume remains larger than or equal to and, consequently, the number of voids is a nondecreasing function. Hence we have the following inequalities for a void’s volume:
(3.3) 
and for the void content at a fixed time :
In the case of constant hydraulic conductivity, , (such a problem often called the HeleShaw problem or the quasistationary Stefan problem) existence and uniqueness of (3.1) have been established both locally and globally and in the classical and weak senses, see [12, 36, 15, 6, 9, 52] and also references therein. Note that in this case no voids can be formed and (3.1) and (3.2) coincide. The models (3.1) and (3.2) with nonconstant can have singulartype behavior when voids form and, in the case of (3.1), also when they collapse. Since there is no collapse of voids in (3.2), behavior of its solutions is less singular than for (3.1) and in this sense (3.2) can be viewed as a regularization of (3.1). We are not aware that questions concerning existence and uniqueness of solutions to (3.1) and (3.2) have been considered in the literature, and they are an interesting and important topic for further study, especially taking into account wide use of such models in applicable sciences. Local existence and uniqueness of solutions to (3.1) and (3.2) can potentially be addressed under some regularity of the data similarly to [33, 34].
In most cases of practical interest permeability (and hence the hydraulic conductivity ) is anisotropic [8, 30, 39, 18, 24, 26]. Consequently, it is important to model the hydraulic conductivity as a random tensor. The principalaxis transformation of the hydraulic conductivity tensor gives
(3.4) 
where is the rotation matrix
We assume that and can be modeled as independent lognormal random fields and the angle is a Gaussian field independent of and We propose that the means of and correspond to the hydraulic conductivity intended by the design of a composite part and the perturbation of these means are stationary random fields modeling uncertainty arising during the manufacturing process. Properties of the model (3.2), (3.4) are discussed in Section 5 below based on its simulation by the CV/FEM algorithm described in the next section.
4 Numerical algorithm
In this section, for completeness of exposition, we give an implementation of the interfacetracking control volume finite element method (CV/FEM) with a fixed grid [2, 20, 45], in a form suitable for the considered stochastic model (3.2), (3.4). CV/FEM is a volumeoffluids technique [20, 43]. It is widely used in the simulation of the RTM filling process [2, 23, 5, 32]. There are a number of alternatives to CV/FEM, including level set methods [38], other volumeoffluids methods [43], marker particle methods (see [38] and references therein), and boundary element methods (see e.g. [53] and references therein). Fixedgrid CV/FEM is currently the method of choice in the RTM community due to its computational efficiency [2] and we follow this common RTM practice here. At the same time, it is of interest in future work to compare CV/FEM with modern level set methods.
Let us turn to the CV/FEM description. The whole computational domain (an empty mold) is first discretized using triangular elements, and then each element is further divided into three subvolumes by connecting the center point and the midpoints of the edges of triangle. Each node is surrounded by a control volume that is composed of all of the subvolumes associated with that node. Note that the number of control volumes is equal to the total number of nodes.
Figure 4.1 shows two different types of tessellation of the spatial domain and the corresponding control volume subdivisions. Suppose that the two horizontal layers divided by a thick line in the middle domain as in Fig. 4.1 have two different permeability values. It is not difficult to see (and was checked experimentally) that in the case of the discretization as in Fig. 4.1(a) the exchange of permeability values between the two layers can result in different filling times due to the discretization asymmetry: the control volume , which has a boundary edge, has two subelements from the bottom layer and one from the top layer. In order to avoid this bias, we choose unbiased (i.e. symmetric) control volumes as shown in Fig. 4.1(b) for our numerical experiments.
In the CV/FEM, to track the interface, a scalar parameter, , called the fill factor, is assigned to each control volume . The fill factor represents the ratio of the volume of fluid to the total volume of the control volume. The fill factor takes values from 0 to 1: for saturated region, for an empty region and for a partially filled region. If , we will say that the control volume is unsaturated. The flow front can be reconstructed based on the nodes that have partially filled control volumes, i.e., those with (see e.g. [20, 3]). In this paper we are not interested in reconstruction of the front and hence we do not consider it here. Finding all parts of the unsaturated regions requires a void detection algorithm. Such an algorithm was introduced in [23] and here we present its implementation. For simplicity, we assume that there is a single vent in the mold but it is not difficult to generalize the algorithm to the case of many vents.
By a void control volume, we understand an unsaturated control volume which is not connected to a vent. To determine which control volumes belong to voids, we first find all unsaturated control volumes which are connected to the vent (i.e., which are not voids) through unsaturated control volumes. To do this, we introduce an array of size equal to the number of nodes, fromvent and two pointers, iter and lastiter, pointing at the first element and the last nonzero element of the array, respectively. This process is done in four steps (see Algorithm 4.1), and Figure 4.2 illustrates it with an example.
Now the void control volumes are all unsaturated control volumes that do not belong to the fromvent formed by Algorithm 4.1. If there are void control volumes at the current time step , then each individual void is represented by a connected set of void control volumes. Starting from any of the void control volumes, , a search process (analogous to Algorithm 4.1) is used to find all void control volumes connected to through void control volumes. The total volume of each void is computed during the search. If the void is formed for the first time at the time then we store its volume in ; otherwise, assign it to . Then the pressure value, , in each void (including its boundary) at the time is determined using the ideal gas law as in (3.2):
(4.1) 
At each time step of the CV/FEM, the fullysaturated control volumes form the solution domain and the finite element method is used to approximate the pressure field in the solution domain. The corresponding boundary conditions on the inlet and the perfectly sealed boundary are imposed. Further, the pressure on the outlet and in allpartially filled or empty control volumes not belonging to voids is set to while the pressure inside voids (i.e., for all unsaturated control volumes belonging to the voids) is set to if the void is formed at this time step and to otherwise.
After computing the approximate pressure field , we calculate the velocity field for completely filled control volumes adjacent to unsaturated control volumes at the centroid of each element using the Darcy’s law:
(4.2) 
where is the superficial fluid velocity.
Assuming that the fluid velocity is constant throughout each element, we compute the local flux into each subelement. Figure 4.3 illustrates the local flux calculation in the triangular element . The midpoints of the three sides of are denoted as , , and and point is the center of the element. The vectors , , and are unit normal vectors perpendicular to the surfaces , , and , respectively. The local flux into each subelement associated with a node in element , , , is calculated as
(4.3)  
where H is the thickness of the mold. The total fluxes entering into the control volume is then calculated by assembling the local fluxes:
(4.4) 
where is the set of elements containing the node .
Having computed all the fluxes, we calculate the new fill factor of the unsaturated control volume as
(4.5) 
where is the volume of control volume . The time increment is calculated so that at least one control volume is filled during the current time step and
(4.6) 
where the minimum is taken over all partiallyfilled control volumes and their neighboring empty control volumes (i.e., all the control volumes which are in a neighborhood of the moving front) at time . Note that the element on which the minimum is reached becomes filled at .
This simulation process continues for every time step until one of two conditions is met: either (1) the entire mold is filled or (2) there is an equilibrium of pressure between the interior and the exterior of all voids and the rest of the mold is filled. The corresponding time is considered as the approximate filling time obtained with the mesh size ( is the length of hypotenuse of the triangular element). To summarize, the CV/FEM is presented in Algorithm 4.2.
To simulate the stochastic model (3.2), (3.4), we need to generate the random hydraulic conductivity tensor at the centers of triangular elements, which requires to sample from stationary Gaussian distributions according to the proposed stochastic model for at the end of Section 3. To generate the required Gaussian random fields, we exploit the block circulant embedding method from [31], which is an extension of the classical circulant embedding method [11, 46] from regular grids to blockregular grids.
5 Numerical experiments
In this section, we present results of numerical experiments which aim at (i) examining convergence of Algorithm 4.2 for the model (3.2) with the quantities of interest being the filling time and void content (Section 5.1); and (ii) studying how variability of permeability affects the filling time in the RTM process (Sections 5.25.4). To experimentally study convergence in Section 5.1, we use the deterministic version of the model (3.2), i.e., when the hydraulic conductivity is a given deterministic tensor. In the stochastic experiments in Sections 5.25.4 we consider the model (3.2), (3.4) with various parameters of the random hydraulic conductivity . The principal conductivity values, and , and the angle are assumed to be independent stationary random fields, with and being lognormal and being Gaussian. In our experiments the three Matérn covariance functions (2.15)(2.17) for the Gaussian random fields , , and , are used with the following scaled distance
(5.1) 
In Sections 5.25.4 we will denote the means of , and by , and , respectively, and the standard deviations of , and by , and , respectively. In each particular experiment the covariance function and the correlation lengths and will be chosen to be the same for all three random fields. To compute expectation and variance of the filling time , we exploit the Monte Carlo technique using independent runs of (3.2), (3.4) in all the Monte Carlo experiments presented in Sections 5.25.4.
In all our experiments the resin is injected into a rectangular mold of size 20cm10cm0.1cm as shown in Figure 5.1. The horizontal direction is viewed as the direction and vertical as the direction. Since in this work we assume that there is no flow in the thickness direction, we use twodimensional elements for modeling the flow. In what follows we denote the mesh size by , which is the length of hypotenuse of the triangular element. The physical properties of the preform are chosen in all the experiments, except the ones with a void in Section 5.1, as listed in Table 5.1.
Porosity  

Injection pressure  MPa 
Vent pressure  MPa 
5.1 Convergence of the CV/FEM
In this section we deal with the deterministic version of the model (3.2), i.e., when the hydraulic conductivity is a given deterministic tensor. We apply Algorithm 4.2 to this model and examine its convergence looking at two quantities of interest: the filling time (it is not random in this experiment, of course) and volume of a void .
To estimate the error of the simulated filling time obtained with the mesh size , we use a reference solution obtained on a grid with the small mesh size cm. Three different shapes of the flow front determined by different choices of the hydraulic conductivity tensor are considered: (a) straight line parallel to yaxis ( m/secPa); (b) cupshaped front ( m/secPa, m/secPa); and (c) capshaped front ( m/secPa, m/secPa). The convergence of the resin filling time with respect to the mesh resolution is shown in Fig. 5.2. We observe that converges approximately quadratically in all three cases. At the same time, we see that the geometry of the flow front has an influence on the accuracy. The most accurate results are seen in the case of the flat flow front (Figure 5.2 (a)) and the least accurate results are seen in the case of the capshaped front (Figure 5.2 (c)).
Void formation in RTM processes has received considerable interest due to its effect on the degradation of physical and mechanical properties of the composite. According to [21, 19], even a void content of just can substantially affect mechanical properties of the material, e.g., decrease of strength up to in bending, in torsional shear, in impact, etc. Consequently, we view that it is important to look at convergence of the CV/FEM for the void content, despite not considering void formation in further experiments here.
Let us look at accuracy of the CV/FEM with presence of a void. In order to have a void in the solution of (3.2), we incorporate a low permeability patch ( and ), where the permeability value is 100 times lower than that of the rest of the mold, m/secPa (see also the corresponding discussion in Section 3). In this setting only a single void can form. We look at the void volume at the time when it is formed and at the time when the mold is fully saturated except the void.
8  2.04e08  2.10e08  4.07e08 

16  6.84e08  7.09e08  1.37e07 
32  8.57e08  8.85e08  1.71e07 
64  9.23e08  9.54e08  1.86e07 
Table 5.2 shows that the inequality (3.3) holds in the discrete case. The volumes of void at two times, and , are both increasing as the mesh size becomes smaller. Two plots in Fig. 5.3 show that the volume of void converges with approximately first order in .
To conclude, we experimentally observed approximately 2nd order convergence of the CV/FEM Algorithm 4.2 for the filling time and a lower order, approximately 1st order, convergence for the volume of void.
5.2 Pseudo1D flow
Now we consider the model (3.2), (3.4) with relatively high mean horizontal conductivity and low mean vertical conductivity (i.e., ), which results in limited movement of flow in the vertical direction. We also choose the horizontal correlation length to be much greater than the vertical correlation length (i.e., ). The parameters of the random conductivity tensor used in this subsection are listed in Table 5.3. We conduct 9 separate numerical experiments using the 3 different values of for each of the three different Matérn covariance functions with different smoothness .
1e8  0 radian  
8e9  0.0356 radian  
1e9  0.01m, 0.03m, 0.05 m  
8e10  0.001m 
The results of the experiments are presented in Tables 5.4 and 5.5. We observe that the mean filling time shows almost no dependence on either the horizontal correlation length or on smoothness of the random hydraulic conductivity field , while the variance of the filling time increases with increase of and . This is consistent with the 1D flow properties studied in Section 2 and one can conclude that when , a good qualitative prediction about the filling time can be made using the analytically solvable onedimensional problem (2.1). Note that the filling time according to the design here and in the corresponding onedimensional case (see Section 2) coincide as expected but it is not so for (see Remark 5.2 below). We also plot typical sample densities for (see Fig. 5.4). We observe that the probability of filling time being in a range close to the designed time (2.8 sec here) is very low.
0.01  33.90.1  33.30.1  33.10.1 
0.03  33.80.2  33.40.2  33.40.2 
0.05  33.90.2  33.60.2  33.70.2 
Remark
We note that though the filling times according to the design in the considered twodimensional case here and in the onedimensional case are both equal to sec, there is a considerable difference between in the two cases: in the pseudo1D flow sec, while in the true 1D case (see Section 2) sec. The reason for this disparity of the mean filling times in 1D and 2D cases can be illustrated in the following way. Let us imagine the extreme case of the pseudo1D flow so that =0 and . Suppose in our domain discretization we have horizonal strips via which resin is propagating. Each strip has its own timetofill on each realization of the random field . Then for we have . It is an interesting probabilistic miniproblem to study the relationship between and , in particular between and , as well as the limiting case , but we do not consider it here. We emphasize that in the deterministic case with the onedimensional model gives good predictions for the twodimensional model and, being considerably simpler, it is often used in practice for this purpose. However, here we highlight that in the stochastic case there are considerable quantitative differences between statistical characteristics of the pseudo1D flow and the 1D flow despite the fact that qualitatively they are similar. This observation is important from the practical point of view. Further, in [22, 50] expansions of moments of the interface dynamics were derived using a dynamical mapping of the Cartesian coordinate system onto a coordinate system associated with the moving front. Making use of the ideas of [22] to obtain expansions for variance of filling times is an interesting problem (even if we assume that discontinuity of the front due to possible void formation can be neglected) for future research.
0.01  12.91.5  13.51.8  15.31.9 
0.03  18.81.5  20.01.6  20.41.5 
0.05  22.71.7  26.22.0  28.42.2 
5.3 Twodimensional isotropic flow
The second stochastic numerical experiment uses the following parameters of the random hydraulic conductivity: , , , radians, and {0.01m, 0.03m, 0.05m}. The focus of this experiment is again to investigate the impact of the correlation length and smoothness of the random field on the filling time.
0.01  37.80.2  39.40.2  39.80.2 
0.03  42.10.5  43.60.5  43.30.5 
0.05  43.30.6  45.20.7  45.40.7 
0.01  32.42.6  35.12.5  35.42.4 
0.03  17012  21617  21920 
0.05  30023  43241  42441 