Stochastic resin transfer molding process

We consider one-dimensional and two-dimensional models of the stochastic resin transfer molding process, which are formulated as random moving boundary problems. We study their properties, analytically in the one-dimensional case and numerically in the two-dimensional case. We show how variability of time to fill depends on correlation lengths and smoothness of a random permeability field.

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 fiber-reinforced composite manufacturing.
We start (Section 2) with considering a stochastic one-dimensional moving-boundary problem, where analytical analysis is possible. The hydraulic conductivity is modelled as a stationary log-normal 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 two-dimensional moving-boundary problem with a random hydraulic conductivity tensor. In some previous related works (see e.g. [42,54]) on stochastic moving-boundary 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 3-dimensional 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 two-phase 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 discrete-type 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 two-dimensional model is solved numerically using the control volume-finite element method (CV/FEM) which we present for completeness in Section 4. The results of our numerical study of the two-dimensional stochastic moving-boundary 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. One-dimensional model. We start with studying a one-dimensional moving-boundary model for a stochastic RTM process. In what follows we assume that we have a sufficiently rich probability space (Ω, F, P ). Let K(x) = K(x; ω) be a random hydraulic conductivity defined on [0, x * ] × Ω and taking values on the positive real semi-line. Assuming the resin's incompressibility, the moving boundary problem for the pressure of resin p(t, x) takes the form [2,42]: x) = 0, 0 < x < L(t), t > 0, (2.1) p(0, x) = p 0 , x ∈ (0, x * ], where L(t) is the moving boundary, κ > 0 is the medium porosity, p I is a pressure on the inlet x = 0, and p 0 ≤ p I is the pressure at the outlet x = x * . We recall that the hydraulic conductivity K can be expressed as where k is the permeability of the medium and µ is the viscosity of the resin.
Remark 2.1. 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 2.2. We note that without losing generality we can put p 0 = 0 in (2.1) (i.e., either it is vacuum on the outlet or p is considered as the relative pressure) but for the sake of the two-dimensional model considered in the next section, it is convenient to keep the parameter p 0 .
It is not difficult to prove the following proposition (see similar statements e.g. in [2,42] Note that Assumption 2.1 implies that G(x), x ≥ 0, is a.s. continuous and strictly increasing, which guarantee existence of the inverse G −1 (·), needed for (2.3). We also remark in passing that if dynamics of the front L(t) 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 x * , i.e., the random time τ such that It follows from (2.3) that It is natural to assume that the mean of K(x) 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 K(x) being a stationary log-normal random field, which is a commonly used assumption (see, e.g. [29,18,24,44,54]), i.e., where K 0 > 0 and ϕ(x) = ϕ(x; ω), (x, ω) ∈ [0, x * ] × Ω, is a stationary Gaussian random field with zero mean and covariance function r(x). Note that The filling time according to the design is equal to The first condition in Assumption 2.1 is obviously satisfied by K(x) from (2.6). To satisfy the second condition, it is sufficient to require that realizations of ϕ(x) are continuous with probability one and we make the following assumption.
Assumption 2.2. We assume that the stationary Gaussian random field ϕ(x) has zero mean and continuous covariance function r(x) such that for some C > 0 and α, δ > 0 : for all x with |x| < δ.
Under Assumption 2.2 the random field ϕ(x) has continuous sample paths with probability one (see e.g. [1]).
For K(x) from (2.6), we obtain the following statistical characteristics of the filling time: To understand the behavior of Eτ and Varτ in terms of K(x), it is convenient to re-write (2.9) via the mean and variance of K(x): It is interesting that the mean filling time depends only on VarK(x) and EK(x) and does not depend on the covariance, and hence it does not depend on the correlation length and smoothness of the random field K(x). We pay attention to the interesting fact that Eτ grows linearly with VarK(x).
The variance of the filling time Varτ depends on covariance of K(x). Then to understand its behavior with respect to correlation length and smoothness of K(x), let us consider some particular cases of ϕ(x).
First let us look at the simple case of ϕ(x) being independent of x, i.e., being just a Gaussian random variable with zero mean and variance σ 2 which can be viewed as a perfectly correlated medium. Then i.e., Varτ grows cubically with increase of VarK in this case.

Remark 2.3.
In this case of a perfectly correlated medium we also have from (2.5): 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 ϕ(x) (see [25,51,41,35]): where Γ(·) is the gamma function, K ν (·) is the modified Bessel function of the second kind, σ 2 is variance, ν > 0 is a smoothness parameter, λ > 0 is a characteristic correlation length, and d(x, λ) is a scaled distance function. In particular, we have The parameter ν in the above covariance functions controls the degree of smoothness of sample paths of the random field. The random field ϕ(x) with Matérn covariance function (2.14) has ν − 1 sample path continuous derivatives with probability one. Hence, with the exponential covariance function r(x) from (2.15), the random field ϕ(x) has continuous (but not differentiable) sample paths with probability one; with r(x) from (2.16) -once differentiable sample paths with probability one; with r(x) from (2.17) -twice differentiable sample paths with probability one; and with r(x) from (2.18) -infinitely many times differentiable sample paths.
Note that we will use the four Matérn covariance functions in this section for x being one-dimensional and in the next sections for x being two-dimensional. The scaled distance of the form d(x, λ) = |x|/λ, with |x| being the usual Euclidean distance, corresponds to an isotropic random field. In the two-dimensional case considered later in the paper we also model the random conductivity as an anisotropic random field, appropriately choosing the scaled distance d(x, λ) (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 Varτ 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 x * /λ >> 1, Varτ grows linearly with λ; (iii) for large correlation lengths x * /λ << 1, Varτ is almost independent of λ (note that in the case of a perfectly correlated medium we had (2.12)); and (iv) for a fixed λ, variance Varτ of the filling time grows with smoothness of the random field. We illustrate these properties in Fig. 2  As an example, we also give the expansion of Varτ in the case of ν = 1/2 and small σ > 0 (expansions in small σ of statistical moments of the interface dynamics and of the pressure were derived in [22]):  To summarize, in the one-dimensional case, • The mean filling time Eτ does not depend on the correlation length λ or on smoothness of the random hydraulic conductivity K(x). It decreases with increase of the mean EK of the hydraulic conductivity and with increase of (p I − p 0 ); and it linearly increases with increase of variance VarK of the hydraulic conductivity and quadratically with increase of the length x * . • The standard deviation of filling time √ Varτ grows as √ λ for small correlation lengths x * /λ >> 1 and is almost independent of λ for large correlation lengths x * /λ << 1. 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 K(x) 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 Eτ 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 fiber-reinforced composite manufacturing as explained in the Introduction.
3. Two-dimensional model. In this section we formulate the two-dimensional analog of the one-dimensional model (2.1). To represent a mold, consider an open two-dimensional domain D with the boundary ∂D = ∂D I ∪ ∂D N ∪ ∂D O , where ∂D I is the inlet, ∂D N is the perfectly sealed boundary, and ∂D O is the outlet. Let K(x, y) = K(x, y; ω) be a random second-order hydraulic conductivity tensor defined onD × Ω. The moving-boundary problem for the pressure of resin p(t, x, y) takes the form (cf. [2,42]): where V (t, x, y) is the velocity of the moving boundary Γ(t) in the normal direction andn(x, y) andn(t, x, y) are the unit outward normals to the corresponding boundaries, D t ∈ D is the time-dependent domain bounded by the moving boundary Γ(t) and the appropriate parts of ∂D, κ > 0 is the medium porosity, p I is a pressure at the inlet ∂D I , and p 0 ≤ p I is the pressure at the outlet ∂D O . Remark 2.1 is applicable here.
Behavior of the two-dimensional model (3.1) is considerably more complex than of the one-dimensional 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 K(x, y) = K 1 I constant everywhere in the domainD (here I is the 2 × 2 unit matrix) except a relatively small region D low which has (again constant) hydraulic conductivity K(x, y) = K 2 I, where K 2 K 1 (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 one-dimensional model (2.1) and Proposition 2.1, it is not difficult to estimate that it is sufficient to have K 2 /K 1 < 1/9 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., p 0 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 two-phase 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 discrete-type formulation of this model, here we give its PDE formulation.
To describe this modification of (3.1), assume that at time t ≥ 0 there are (t) entrapments with closed boundaries Γ i (t) and volumes v i (t) formed at times τ i , i = 1, . . . , (t), behind the main front Γ 0 (t). Then we can write the model as x, y) are the velocities of the moving boundaries Γ i (t) in the normal direction and n(x, y) andn i (t, x, y) are the unit normals to the corresponding boundaries, D t ∈ D is the time-dependent domain bounded by Γ i (t), i = 0, . . . , (t), and the appropriate parts of ∂D, the rest of the notation is as in (3.1). The volume v i (t) of a void is computed as where H 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 two-dimensional 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 τ i , its volume v i (t) remains larger than or equal to p 0 v i (τ i )/p I and, consequently, the number of voids (t) is a non-decreasing function. Hence we have the following inequalities for a void's volume: and for the void content at a fixed time T : In the case of constant hydraulic conductivity, K(x) = K, (such a problem often called the Hele-Shaw problem or the quasi-stationary 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 non-constant K(x) can have singular-type 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 k(x, y) (and hence the hydraulic conductivity K(x, y)) is anisotropic [8,30,39,18,24,26]. Consequently, it is important to model the hydraulic conductivity K(x, y) as a random tensor. The principal-axis transformation of the hydraulic conductivity tensor gives We assume that K xx (x, y) = K xx (x, y; ω) and K yy (x, y) = K yy (x, y; ω) can be modeled as independent log-normal random fields and the angle θ(x, y) = θ(x, y; ω) is a Gaussian field independent of K xx (x, y) and K yy (x, y). We propose that the means of K xx (x, y), K yy (x, y) and θ(x, y) 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 interface-tracking 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 volume-of-fluids 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 volume-of-fluids methods [43], marker particle methods (see [38] and references therein), and boundary element methods (see e.g. [53] and references therein). Fixed-grid 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 sub-volumes 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.  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 CV i , 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, f i , called the fill factor, is assigned to each control volume CV i . The fill factor represents the ratio of the volume of fluid to the total volume of the control volume. The fill factor f i takes values from 0 to 1: f i = 1 for saturated region, f i = 0 for an empty region and 0 < f i < 1 for a partially filled region.
If 0 ≤ f i < 1, 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  By a void control volume, we understand an unsaturated control volume which is not 13 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, from vent and two pointers, iter and last iter, pointing at the first element and the last nonzero element of the array, respectively. This process is done in four steps (see  Step 1 Add all unsaturated control volumes where the vent is located into the array from vent, and set two pointers iter and last iter (Figure 4.2 (b)).
Step 2 The loop is carried out over the neighbors of a control volume CV i pointed at by the pointer iter to detect first visited neighbors whose fill factors are less than 1. Each time the neighbor is chosen to be added to from vent, move last iter to the next element ( Figure 4.2 (c)).
Step 3 After all neighbors of CV i are checked, advance iter to the next element ( Now the void control volumes are all unsaturated control volumes that do not belong to the from vent formed by Algorithm 4.1. If there are void control volumes at the current time step t k , then each individual void is represented by a connected set of void control volumes. Starting from any of the void control volumes, CV i , a search process (analogous to Algorithm 4.1) is used to find all void control volumes connected to CV i through void control volumes. The total volume of each void j is computed during the search. If the void is formed for the first time at the time t k then we store its volume inv * j ; otherwise, assign it tov j (t k ). Then the pressure value,p void j (t k ), in each void (including its boundary) at the time t k is determined using the ideal gas law as in (3.2): At each time step t k of the CV/FEM, the fully-saturated control volumes form the solution domain and the finite element method is used to approximate the pressure fieldp in the solution domain. The corresponding boundary conditions on the inlet ∂D I and the perfectly sealed boundary ∂D N are imposed. Further, the pressure on the outlet ∂D O and in all-partially filled or empty control volumes not belonging to voids is set to p 0 while the pressure inside voids (i.e., for all unsaturated control volumes belonging to the voids) is set to p 0 if the void is formed at this time step and top void j (t k ) otherwise.
After computing the approximate pressure fieldp, 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: where u 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 e i . The midpoints of the three sides of e i are denoted as a, b, and c and point o is the center of the element. The vectors n 1 , n 2 , and n 3 are unit normal vectors perpendicular to the surfaces oa, ob, and oc, respectively. The local flux into each subelement associated with a node x i in element e j , Q e j ,x i , i = 1, 2, 3, is calculated as where H is the thickness of the mold. The total fluxes entering into the control volume CV i is then calculated by assembling the local fluxes: where E i is the set of elements containing the node x i . Having computed all the fluxes, we calculate the new fill factor f i (t k+1 ) = f i (t k + δt) of the unsaturated control volume i as where V i is the volume of control volume i. The time increment δt k is calculated so that at least one control volume is filled during the current time step and where the minimum is taken over all partially-filled 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 t k . Note that the element on which the minimum is reached becomes filled at t k+1 = t k + δt k . 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 τ h is considered as the approximate filling time obtained with the mesh size h (h is the length of hypotenuse of the triangular element). To summarize, the CV/FEM is presented in Algorithm 4.2.
Step 2. Identify the saturated domain by the fill factors and detect void control volumes using Algorithm 4.1.
Step 3. Compute the pressure fieldp in the saturated domain with the appropriate boundary conditions by the FEM.
Step 4. Calculate the pressure gradients, ∇p, in the neighborhood of the flow front by differentiating the element shape functions and then compute the velocity field, u, using Darcy's law (4.2).
Step 6. Find the size of the time step δt k by (4.6).
Step 7. Update the fill factors with (4.5) and advance the time: t k+1 = t k + δt k and set k = k + 1.
Step 8. Repeat from Step 2 for the newly-filled domain until the mold is completely saturated or there is an equilibrium of pressure between the interior and the exterior of all voids and the rest of the mold is filled.
To simulate the stochastic model (3.2), (3.4), we need to generate the random hydraulic conductivity tensor K(x, y) at the centers of triangular elements, which requires to sample from stationary Gaussian distributions according to the proposed stochastic model for K(x, y) 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 block-regular grids.  K(x, y). The principal conductivity values, K xx and K yy , and the angle θ(x, y) are assumed to be independent stationary random fields, with K xx and K yy being log-normal and θ(x, y) being Gaussian. In our experiments the three Matérn covariance functions (2.15)-(2.17) for the Gaussian random fields log K α , α = xx, yy, and θ, are used with the following scaled distance In Sections 5.2-5.4 we will denote the means of K xx , K yy and θ by µ Kxx , µ Kyy and µ θ , respectively, and the standard deviations of K xx , K yy and θ by σ Kxx , σ Kyy and σ θ , respectively. In each particular experiment the covariance function and the correlation lengths λ x and λ y 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 3000 independent runs of (3.2), (3.4) in all the Monte Carlo experiments presented in Sections 5. 2-5.4. In all our experiments the resin is injected into a rectangular mold of size 20cm×10cm ×0.1cm as shown in Figure 5.1. The horizontal direction is viewed as the x-direction and vertical as the y-direction. Since in this work we assume that there is no flow in the thickness direction, we use two-dimensional elements for modeling the flow. In what follows we denote the mesh size by h, 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.

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 K(x, y) 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 v(t).
To estimate the error of the simulated filling time τ h obtained with the mesh size h, we use a reference solution τ h * obtained on a grid with the small mesh size h * = 10/128 cm. Three different shapes of the flow front determined by different choices of the hydraulic conductivity K(x, y) tensor are considered: (a) straight line parallel to y-axis (K xx = K yy = 10 −8 m 2 /sec·Pa); (b) cup-shaped front (K xx (x, y) = 10 −9 + 10 −8 (1 − sin(πy/0.1)) m 2 /sec·Pa, K yy = 10 −10 m 2 /sec·Pa); and (c) cap-shaped front (K xx (x, y) = 10 −9 + 10 −8 sin(πy/0.1) m 2 /sec·Pa, K yy = 10 −10 m 2 /sec·Pa). The convergence of the resin filling time τ h with respect to the mesh resolution h is shown in Fig. 5.2. We observe that |τ h * − τ h |/τ h * 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 cap-shaped 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 1% can substantially affect mechanical properties of the material, e.g., decrease of strength up to 30% in bending, 9% in torsional shear, 8% 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 (1cm ≤ x ≤ 3cm and 4cm ≤ y ≤ 6cm), where the permeability value is 100 times lower than that of the rest of the mold, 10 −7 m 2 /sec·Pa (see also the corresponding discussion in Section 3). In this setting only a single void can form. We look at the void volumev(t) at the time t * when it is formed and at the time t end when the mold is fully saturated except the void.

5.2.
Pseudo-1D flow. Now we consider the model (3.2), (3.4) with relatively high mean horizontal conductivity µ Kxx and low mean vertical conductivity µ Kyy (i.e., µ Kxx >> µ Kyy ), which results in limited movement of flow in the vertical direction. We also choose the horizontal correlation length λ x to be much greater than the vertical correlation length λ y (i.e., 10  λ x >> λ y ). The parameters of the random conductivity tensor K(x, y) used in this subsection are listed in Table 5.3. We conduct 9 separate numerical experiments using the 3 different values of λ x for each of the three different Matérn covariance functions with different smoothness ν.
The results of the experiments are presented in Tables 5.4 and 5.5. We observe that the mean filling time Eτ shows almost no dependence on either the horizontal correlation length λ x or on smoothness ν of the random hydraulic conductivity field K(x, y), while the variance of the filling time increases with increase of λ x and ν. This is consistent with the 1D flow properties studied in Section 2 and one can conclude that when µ Kxx >> µ Kyy , a good qualitative prediction about the filling time τ can be made using the analytically solvable one-dimensional problem (2.1). Note that the filling time according to the design here and in the corresponding one-dimensional case (see Section 2) coincide as expected but it is not so for Eτ (see Remark 5.1 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  Table 5.4: Mean filling time Eτ (in sec) and its 95% confidence interval for the pseudo-1D flow. The filling time according to the design is equal to 2.8 sec.
Remark 5.1. We note that though the filling times according to the design in the considered two-dimensional case here and in the one-dimensional case are both equal to 2.8 sec, there is a considerable difference between Eτ in the two cases: in the pseudo-1D flow Eτ ≈ 34 sec, while in the true 1-D case (see Section 2) Eτ ≈ 4.6 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 pseudo-1D flow so that K yy =0 and θ = 0. Suppose in our domain discretization we have l horizonal strips via which resin is propagating. Each strip has its own time-to-fill τ i on each realization of the random field K(x). Then for τ = max 1≤i≤l τ i we have Eτ ≥ Eτ i . It is an interesting probabilistic mini-problem to study the relationship between τ and τ i , in particular between Eτ and Eτ i , as well as the limiting case l → ∞, but we do not consider it here. We emphasize that in the deterministic case with K xx >> K yy the one-dimensional model gives good predictions for the two-dimensional 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 pseudo-1D 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.  Table 5.5: Variance of the filling time Varτ (in sec 2 ) and its 95% confidence interval for the pseudo-1D flow.
Unlike the pseudo-1D results presented in Section 5.2, where the mean filling time Eτ is independent of the correlation length and smoothness, here Eτ increases with growth of the correlation length as shown in Table 5.6 and slightly grows with increase of smoothness. In other words, 2D flow moves slower on a more homogeneous porous medium with less spatial variability than in the pseudo-1D case. We also observe a fast increase of variance Varτ with an increasing spatial correlation length in Table 5.7.
5.4. Two-dimensional flow with anisotropic mean. The third stochastic numerical experiment uses the following parameters of the random hydraulic conductivity: µ Kxx = 1e-8m 2 /sec · Pa, µ Kyy = 2e-8m 2 /sec · Pa, σ Kxx = σ Kyy = 8e-9m 2 /sec · Pa, µ θ = 0, σ θ = 0.0356 radian, and λ x = λ y = {0.01m, 0.03m, 0.05m}. In other words, here (in comparison with Section 5.3) we consider a model in which random hydraulic conductivity has an anisotropic mean. Recall (see Section 3) that we proposed that the means of K xx (x, y), K yy (x, y) and θ(x, y) correspond to the hydraulic conductivity intended by the design of a composite part and that in most cases of practical interest permeability k(x, y) (and hence the hydraulic conductivity K(x, y)) is anisotropic by design. Note that the mean of permeability in the x-direction (K xx ) is smaller than that in the y-direction (K yy ) and that the 1D model from Section 2 does not apply here (see the discussion in Section 5.2).
The results of the numerical experiments are presented in Tables 5.8 and 5.9. We observe that the mean filling time Eτ of the 2D mean anisotropic flow increases as the correlation length increases and also that an increasing spatial correlation length leads to an increase of Varτ .
Eτ λ x = λ y ν = 1/2 ν = 3/2 ν = 5/2 0.01 35 We note that in all three stochastic cases considered in Sections 5.2-5.4 the mean filling time is considerably higher that the filling time according to the design, which further emphasises the importance of stochastic modeling of RTM.
Remark 5.2. The MATLAB codes for the CV/FEM used in our experiments are available at https://github.com/parkmh/MATCVFEM.
6. Discussion and summary. In this work we considered stochastic one-dimensional and two-dimensional moving-boundary problems suitable for modeling RTM processes. The PDE formulation of the main two-dimensional model which takes into account compressible air entrapment in voids behind the main front is somewhat novel. The one-dimensional problem has an analytical solution while the two-dimensional model requires a numerical method for computing quantities of interest. Following the common practice in the RTM community [2], we use the control volume-finite element method (CV/FEM) to simulate the two-dimensional problem. We test accuracy of the CV/FEM algorithm for three particular cases of deterministic hydraulic conductivity, and we experimentally observed approximately its 2nd order convergence for the filling time and a lower order, approximately 1st order, convergence for the volume of void. We note that these tests were done in the case of infinitely smooth (in the case of the filling time) or piece-wise constant (in the case of the volume of void) hydraulic conductivity while less regular random field are used to model hydraulic conductivity. For future study, it is of interest to look at dependence of behavior of the CV/FEM algorithm on the smoothness of hydraulic conductivity and on various observables.
We studied properties of a stochastic one-dimensional moving-boundary problem with the hydraulic conductivity being modelled as a stationary log-normal random field. In particular, we observe that the mean filling time does not depend on correlation length or on smoothness of hydraulic conductivity, while the filling time variance (as a measure of its variability) does depend on both the correlation length and smoothness. A similar conclusion is made about the two-dimensional model's behavior when random permeability in the direction from inlet to outlet is much bigger than in the perpendicular direction (the case of pseudo-1D flow). However, we discovered that in other cases of random permeability the mean filling time does depend on correlation length and on smoothness. The important consequence of these conclusions is the observed (often high) sensitivity of the mean and variance of the filling time to changes in correlation length of the permeability as well as in its smoothness. This highlights the importance of conducting laboratory experiments from which covariance of the permeability field can be reconstructed. Further, sensitivity of filling time to smoothness of permeability serves as a warning that stochastic modeling of permeability via homogenization procedures needs to be done with a very careful choice of scales.
Among the main objectives of this paper was to attract attention of the UQ community to challenges posed by stochastic moving-boundary problems, which are highly relevant to modern technological processes in production of composite materials as well as to other porous media problems. The challenges include (i) establishing existence and uniqueness results for two and three dimensional moving-boundary problems of the type considered in this work; (ii) numerical analysis for the corresponding CV/FEM algorithm; (iii) development of faster sampling techniques, e.g. using a multi-level Monte Carlo approach [17,27,28] and/or polynomial chaos expansions [16,48,47]; (iv) computational experiments with complex geometry using outcomes of (iii); (v) comparing the CV/FEM with other numerical approaches to moving boundary problems, e.g. level sets methods [38]; (vi) design of laboratory experiments and collection and analysis of the corresponding data to recover characteristics of random permeability and to find a stochastic model of the conductivity consistent with experimental data (for recent research in this direction see e.g. [2,24,26]); and (vii) considering a two-phase model involving both incompressible (resin) and compressible (air) phases and comparing its properties with the ones for the reduced model studied here.