Bayesian Optimal Experimental Design for Race Tracking in Resin Transfer Moulding

Featured Application

Predict race-tracking strength in an RTM process using minimal pressure sensor measurements and position the sensors optimally throughout the preform in order to do so.

Abstract

A Bayesian inference formulation is applied to the Resin Transfer Moulding process to estimate bulk permeability and race-tracking effects using measured values of pressure at discrete sensor locations throughout a preform. The algorithm quantifies uncertainty in both the permeability and race-tracking effects, which decreases when more sensors are used or the preform geometry is less complex. We show that this approach becomes less reliable with a smaller resin exit vent. Numerical experiments show that the formulation can accurately predict race-tracking effects with few measurements. A Bayesian A-optimality formulation is used to develop a method for producing optimal sensor locations that reduce the uncertainty in the permeability and race-tracking estimates the most. This method is applied to two numerical examples which show that optimal designs reduce uncertainty by up to an order of magnitude compared to a random design.

1. Introduction

Resin Transfer Moulding (RTM) is a common advanced processing method for producing near-net-shaped composite parts. It consists of five key stages:

• Preform Manufacturing: The reinforcing fabric is manufactured, typically by stacking layers of fabric upon one another.
• Lay-up and Draping: The fabric preform is layed up or draped in the mould.
• Mould Closure: The fabric layers are compressed to the final thickness of the part.
• Resin Injection: Resin is injected into the mould through the injection gate(s) until the preform is fully saturated and the resin exits through a vent.
• Cure and Demoulding: A part is removed from the mould after the resin cures.

RTM is often used within the industry because of the relatively short production times and simple experimental setup to construct near-net-shaped parts. However, these advantages come with caveats.

An important stage in RTM is resin injection, where the flow of the resin must be carefully controlled to avoid incomplete fillings, in which air can be trapped inside the mould, forming macroscopic dry spots. Even very small regions without resin can cause catastrophic failure of the composite part [1,2,3,4]. Accurate predictions of resin flow allow for strategic design of gates, vents, and injection schemes to optimally fill the composite part, with minimal chance of dry spot formation [1].

The flow of the resin through the fabric preform depends on the permeability. Deviations in the permeability from what is expected often come from the lay-up and draping stage of manufacturing, when the preform is placed into the mould. Occasionally, there can be areas of very high permeability at the boundaries of the preform where the fabric density is lower. These are ‘race tracks’ for the resin as it fills the mould. The strength and existence of race tracking are influenced by a variety of factors, including fabric type, the preform manufacturing method, and the placement of the preform into the mould. It can vary from one part to the next in the same production run and is usually not repeatable [5].

The unpredictable nature of race tracks, in terms of their strength, can lead to significant variations in resin propagation, significantly changing the flow front movement pattern [6,7,8]. This can lead to dry spots in the composite parts, as the flow moves in an unpredictable way, trapping air. It is of great interest to be able to predict the effects of race tracking and mitigate the chance of dry spot formation, even in cases of high race tracking. When race tracking is predicted, one can use it as a method for advantageous resin distribution, rather than race tracking being an unpredictable problem.

Research around the effects of race tracking is not mainly focused on the prevention of race tracking, which is inevitable to some extent, but on understanding and predicting the race-tracking effect. Techniques from areas of statistics and machine learning have been incorporated into prediction, as well as methods that incorporate existing physical models of RTM [8,9,10,11,12].

One significant example, [5], addresses the problem of predicting race tracking online, as the resin is injected. This is an important problem as it allows manufacturers to accurately predict the flow of the resin, taking into account the specific race-tracking scenario for each part. Online automation tools could then be used to change the resin flow to ensure that there are no dry spots by introducing control actions during filling, such as introducing new injection locations, opening and closing exit vents, and changing injection pressures [13]. The method is based on discretising potential race-tracking effects for each susceptible region and creating a database of resin flow for each permutation of potential effects for each region. As a part is filled, sensors collect data that are then compared to the closest scenario in the database. This method is efficient because all the extensive computation can be completed offline. However, by considering only discrete scenarios and approximating flow based on the closest scenario, the reality could differ (possibly significantly) from what is predicted.

To allow for the systematic incorporation of uncertainties with race-tracking prediction, we consider this inverse problem in the Bayesian framework [14,15,16,17]. In this framework, the solution to the underlying statistical inverse problem is given by the posterior probability density. For nonlinear inverse problems with expensive forward models and high-dimensional parameters (as is the case for RTM inverse problems), fully characterising the posterior is typically not tractable. Consequently, we compute the Laplace approximation of the posterior, which requires only the maximum a posteriori (MAP) estimate, i.e., permeability and race-tracking strength, which maximise the posterior density and the approximate posterior covariance.

The idea of incorporating Bayesian inference into this problem has been explored previously [18,19]. These investigations focus on the more general problem of predicting a complete permeability field throughout the mould. The purpose of this would be to use pressure and flow front data created during filling to predict defects in the final part after it is produced but before any structural testing. They showed that a permeability and porosity field can be predicted using Bayesian inference, including predicting potential defects. However, extracting enough data with multiple pressure sensors and linear flow front sensors could be considered impractical and tedious in real manufacturing cases.

The current work uses a novel, lower-dimensional parameterisation compared to that formulation [18,19], where permeability throughout the domain is characterised by a homogeneous region, along with homogeneous race-tracking regions with a greater permeability. The purpose of this work is to estimate the race tracking that occurs during mould filling. The parameterisation presented here allows for the efficient estimation of race-tracking strength while neglecting other smaller spatially varying permeability effects. The simplification to a lower-dimensional parameterisation also allows for much faster computation, which opens the potential for online calculation in the future. In this paper, we use offline calculations to compute the optimal locations for pressure sensors, a problem that has not yet been considered to the best of the knowledge of the authors. These locations could be used to predict race tracking online in the future. It is the hope that an optimally designed experimental setup will allow practitioners to extract maximum insight into the permeability, and therefore flow pattern, as the mould is being filled. From there, control actions can be taken, as has been investigated before [5,13,20,21], to ensure that the resin flow is controlled and there are minimal dry spots.

2. Computational Model

2.1. Constitutive Equations

This paper is limited to considering the 2D case, which is sufficient to model the vast majority of net-shaped parts for which RTM is used. However, the computational model and all Bayesian formulations in this model can be applied to higher-dimensional formulations.

For the slow-moving thermoset resins that are typically used in RTM, Darcy’s law adequately describes the 2D flow through a porous medium [1,5,13,20,22,23,24]:

$$\begin{array}{cc}\hfill \mathit{q}& =-\frac{\mathit{K}}{\mu }\nabla p,\hfill \end{array}$$

(1)

where $\mathit{q}\in {\mathbb{R}}^2$ is the flux vector of resin, p is the pressure, $\mu$ is the fluid viscosity, and $\mathit{K}\in {\mathbb{R}}^{2\times 2}$ is the permeability tensor, which is symmetric and positive definite. Combining Darcy’s law with conservation of mass gives:

$$\begin{array}{cc}\hfill -\nabla ·\left(\frac{\mathit{K}}{\mu }\nabla p\right)& =0.\hfill \end{array}$$

(2)

We base the following notation on the notation used in [19] to describe this problem. The domain of the mould D can be partitioned at any point in time into $F\left(t\right)$ and $U\left(t\right)$, which represent the filled and unfilled regions of the domain, respectively. The boundary of D is $\delta D=\delta D_I\cup \delta D_N\cup D_O$, where $\delta D_I$ represents the inlet, $\delta D_N$ represents the impermeable walls of the mould, and $\delta D_O$ is the outlet vent. We consider the case where resin is injected from a gate at a fixed constant pressure ($p_{\mathrm{in}}$) and exits from a vent. At this vent, resin pressure is fixed at atmospheric pressure, taken to be 0 without loss of generality. The walls of the mould have zero normal flow. All unfilled regions are also fixed at zero pressure, neglecting the scenario in which air is sealed and compressed in the mould, which we are trying to prevent. In summary, the problem has the following boundary conditions, where $\mathit{n}\left(\mathit{x}\right)$ is the normal vector to the domain boundary at $\mathit{x}$:

$$\begin{array}{ccc}\hfill p(\mathit{x},t)& =p_{\mathrm{in}},\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{x}\in \delta D_I,\end{array}$$

(3)

$$\begin{array}{ccc}\hfill p(\mathit{x},t)& =0,\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{x}\in \delta D_O,\end{array}$$

(4)

$$\begin{array}{ccc}\hfill \nabla p(\mathit{x},t)·\mathit{n}\left(\mathit{x}\right)& =0,\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{x}\in \delta D_N,\end{array}$$

(5)

$$\begin{array}{ccc}\hfill p(\mathit{x},t)& =0,\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{x}\in U,\end{array}$$

(6)

and the following initial conditions:

$$\begin{array}{c}\hfill U\left(0\right)=D-\delta D_I,(F\left(0\right)=\delta D_I).\end{array}$$

(7)

2.2. Numerical RTM Model

The RTM simulation developed here uses a standard multigrid method [25,26,27]. This uses control volumes (CVs) that keep track of resin saturation, offset from the triangular finite elements that are used to calculate the pressure through Equation (2). These are illustrated in Figure 1. This paper will refer to these as ‘control volumes’ (CVs) and ‘elements’, respectively, to avoid confusion.

The algorithm can be divided into a series of steps:

• Solve finite element equations for pressures at the nodes.
• Calculate the flow between each control volume.
• Calculate the time to fill the next control volume.
• Step forward in time and propagate the resin front forward.

These steps are repeated until all control volumes are full. This algorithm is available in LIMS [28,29], which can address complex geometries and 3D flows. However, for this work, the goal was to explore the A-optimality formulation in combination with the Bayesian framework. Once demonstrated, this can be implemented in other RTM simulations. As such, our numerical simulations in this paper use an implementation of this algorithm written in Julia code. Each of the steps of the algorithm is discussed in more detail in the following.

2.2.1. Solve for Pressure

Equation (2) can be written in a weak form using the Galerkin finite element method, for a given mesh. The derivation of this is given in Appendix A. This gives the stiffness equation to solve:

$$\begin{array}{c}\hfill \mathit{E}\left(\mathit{K}\right)\mathit{p}=\mathit{f},\end{array}$$

(8)

where $\mathit{E}\left(\mathit{K}\right)$ and $\mathit{f}$ are the finite element stiffness matrix and the forcing vector, respectively, as defined in Appendix A. All nodes with respective CVs in $U\left(t\right)$ are fixed at $p=0$, as discussed in Section 2.1. Applying boundary conditions based on these time-varying sets alters $\mathit{E}\left(\mathit{K}\right)$ to $\mathit{A}(\mathit{K},t)$, which is defined as follows, using MATLAB notation:

$$\begin{array}{cc}\hfill \mathit{A}(i,:)=\left\{\begin{array}{c}\mathit{E}(i,:)\hfill \\ {\mathit{e}}_i^T\hfill \end{array}\right.& \begin{array}{c}C{V}_i\in \mathcal{F},\hfill \\ C{V}_i\in U\left(t\right)\mathrm{or}{n}_i\cup \delta D_I\cup \delta D_O),\hfill \end{array}\hfill \end{array}$$

(9)

where $C{V}_i$ represents the ith CV, $n_i$ represents the ith node, and ${\mathit{e}}_i^T$ is the standard unit basis vector. The vector $\mathit{f}$ does not change as $\mathcal{U}$ and $\mathcal{F}$ change and is defined as follows:

$$\begin{array}{cc}\hfill \mathit{f}\left(i\right)=\left\{\begin{array}{c}{p}_{\mathrm{in}}\hfill \\ 0\hfill \end{array}\right.& \begin{array}{c}{n}_i\in \delta D_I,\hfill \\ \mathrm{otherwise}.\hfill \end{array}\hfill \end{array}$$

(10)

Equation (8) is never actually solved, but, instead, the following is solved every timestep:

$$\begin{array}{c}\hfill \mathit{A}(\mathit{K},t)\mathit{p}\left(t\right)=\mathit{f}.\end{array}$$

(11)

It is worth noting that $\mathit{A}(\mathit{K},t)$ is generally very sparse.

2.2.2. Calculate Flow

Given the finite element pressure distribution, the flow between CVs is calculated. The flow between filled CVs is not important, as the net flux will add to zero. $C{V}_i$ has a saturation ${\Psi }_i$ where

• ${\Psi }_i=0$ represents an empty CV.
• ${\Psi }_i=1$ represents a filled CV.
• $0<{\Psi }_i<1$ is a partially filled CV.

CVs that increase in saturation during the current timestep must be considered. To this end, we define the set X that contains all the values for i such that $C{V}_i$ has an adjacent CV that is filled and ${\Psi }_i<1$.

For a given finite element, we can express Darcy’s law for the flow rate $\mathit{q}$—Equation (2)—in terms of the nodal pressures:

$$\begin{array}{cc}\hfill \mathit{q}& =-\frac{1}{\mu }\mathit{K}\left[\begin{array}{c}{\displaystyle \frac{\partial p}{\partial x}}\\ {\displaystyle \frac{\partial p}{\partial y}}\end{array}\right],\left[\begin{array}{c}{\displaystyle \frac{\partial p}{\partial x}}\\ {\displaystyle \frac{\partial p}{\partial y}}\end{array}\right]={\mathit{J}}^{-1}\left[\begin{array}{c}\sum _{i=1}^3{p}_i{\displaystyle \frac{\partial {\varphi }_i}{\partial {\xi }_1}}\\ \sum _{i=1}^3{p}_i{\displaystyle \frac{\partial {\varphi }_i}{\partial {\xi }_2}}\end{array}\right].\hfill \end{array}$$

(12)

Here, $\mathit{J}$ is the finite element Jacobian, defined as follows, based on the three coordinates of the element $(x_1,y_1)$, $(x_2,y_2)$, and $(x_3,y_3)$:

$$\begin{array}{c}\hfill \mathit{J}=\left[\begin{array}{cc}{x}_1-x_3& y_1-y_3\\ x_2-x_3& y_2-y_3\end{array}\right].\end{array}$$

(13)

Note the lack of dependence of $\mathit{q}$ on the element coordinates—the linear elements for pressure give a constant flow rate throughout an element. Through a boundary C, the volume flow rate $Q_C$ can be obtained by integrating $\mathit{q}·\mathit{n}$ across the boundary. Within one element, there is no spatial dependence, so we have the following for $Q_C$:

$$\begin{array}{cc}\hfill Q_C& =(\mathit{q}·\mathit{n})\left|C\right|,\hfill \end{array}$$

(14)

where $\left|C\right|$ is the length of the boundary C. By summing $Q_C$ for all boundaries for the ith CV, the total flow $Q_i$ is calculated.

2.2.3. Time to Fill Next CV

Once the total flow $Q_i$ is found, the time to fill the CV can be calculated based on the current saturation. For all $i\in X$:

$$\begin{array}{cc}\hfill t_{\mathrm{fill}}^{\left(i\right)}& =\frac{\varphi (1-{\Psi }_i)A_i}{Q_i},\hfill \end{array}$$

(15)

where $\varphi$ is the porosity of the domain and $A_i$ is the area of $C{V}_i$. We then know the timestep to propagate to the next iteration:

$$\begin{array}{cc}\hfill \Delta t& =\underset{i\in X}{min}{t}_{\mathrm{fill}}^{\left(i\right)},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill i^{\prime }& =\underset{i\in X}{\mathrm{argmin}}{t}_{\mathrm{fill}}^{\left(i\right)}.\hfill \end{array}$$

(17)

2.2.4. Timestep

The time dependence in this problem comes only from the moving boundary, causing a change in $U\left(t\right)$ and $F\left(t\right)$. At the end of the iteration of this method, the following change takes place:

$$\begin{array}{cc}\hfill U(t+\delta t)& =U\left(t\right)-C{V}_i,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill F(t+\delta t)& =F\left(t\right)+C{V}_i.\hfill \end{array}$$

(19)

If there are still unfilled CVs, we return to solving for pressure (Section 2.2.1). This begins with a change from $\mathit{A}(\mathit{K},t)$ to $\mathit{A}(\mathit{K},t+\delta t)$ (Equation (9)).

2.2.5. Interpolation

For this investigation, it is of interest to be able to output the measurements that pressure sensors at specific locations would give at specified measurement times. The finite element solution defines the pressure at every location in the domain. However, a disadvantage of the algorithm described here is that the timesteps are dependent on the filling times. Therefore, extracting pressure measurements at specific times requires linear interpolation.

3. Bayesian Formulation

3.1. Permeability Vector

The unknown model parameters that represent the permeability are described by $\mathit{k}=[c_1,c_2,c_3,{\alpha }_1,{\alpha }_2,\dots ,{\alpha }_n]$. We use $c_1$, $c_2$, and $c_3$ to represent the log-Cholesky factorisation of the bulk permeability tensor, ${\mathit{K}}_{\mathrm{bulk}}$—see (20). ${\mathit{K}}_{\mathrm{bulk}}$ represents the permeability tensor for the majority of the domain, aside from the region experiencing race tracking. Specifically, we take:

$$\begin{array}{c}\hfill {\mathit{C}}^T\mathit{C}={\mathit{K}}_{\mathrm{bulk}},\mathit{C}=\left[\begin{array}{cc}{e}^{c_1}& c_2\\ 0& e^{c_3}\end{array}\right].\end{array}$$

(20)

Parameterising the bulk permeability in this way ensures that the bulk permeability tensor is symmetric and positive definite, regardless of the values of $c_1$, $c_2$, and $c_3$ [30]. The remaining parameters ${\alpha }_1,{\alpha }_2,\dots ,{\alpha }_n$ represent race tracking in various areas. For given values of these parameters, the permeability tensor of each area can be constructed as follows.

$$\begin{array}{cc}\hfill {\mathit{K}}_j& =exp\left({\alpha }_j\right){\mathit{K}}_{\mathrm{bulk}},j=[1,2,\dots ,n],\hfill \end{array}$$

(21)

where ${\mathit{K}}_j$ is the permeability tensor in a specific race-tracking region. The use of the exponential function is convenient, as it allows the parameters ${\alpha }_j$ to be normally distributed, but gives an exponential distribution in the race-tracking multipliers. This exponential distribution reflects typical race-tracking multipliers, which are often 1–5 times larger than the bulk permeability but can be 100 times higher in some instances [1,5,13,20]. The benefit of using normally distributed random parameters will be discussed in the following sections.

It should be noted that we have limited the stochastic variables in our model to only those parameterised by k. This includes the bulk permeability and the race-tracking strength of each region. We fix the porosity and race-tracking width parameters. The reasons for excluding each are different. Stochastic porosity has been considered in a previous study [19] but will be neglected here due to its typical low variation and because the focus of this research is specifically on race tracking. We neglect the race-tracking width, although stochastic in actuality, because the dual effect of the race-tracking width and strength can be replicated with just strength. As the interest in race tracking is in its effect rather than describing the exact source, this is an appropriate assumption.

3.2. Bayesian Inference

The objective of this investigation is to estimate the permeability vector $\mathit{k}$, which exists and produces a vector of measurements $\mathit{d}$, as the composite part is produced. For the RTM problem, the length of this vector is the number of measurement times multiplied by the number of pressure sensors. We assume that these data contain some error or noise $\mathit{e}$, so we have:

$$\begin{array}{c}\hfill \mathit{d}=\mathcal{G}\left(\mathit{k}\right)+\mathit{e},\mathit{e}\sim \mathcal{N}(\mathbf{0},{\mathbf{\Gamma }}_e),\end{array}$$

(22)

where ${\mathbf{\Gamma }}_e$ is a covariance matrix that describes the distribution of this error, and $\mathcal{G}$ represents the ‘forward model’. That is, $\mathcal{G}\left(\mathit{k}\right)$ produces noise-free measurements that would be taken when completing an RTM resin injection with permeability vector $\mathit{k}$. The computational model described in Section 2 is used every time $\mathcal{G}\left(\mathit{k}\right)$ is calculated.

This estimation of $\mathit{k}$ can be achieved by Bayesian inference. We treat the unknown $\mathit{k}$ as a random variable and determine the posterior probability distribution $\pi \left(\mathit{k}\right|\mathit{d})$. This has the added benefit of automatically quantifying the uncertainty in the parameters. The posterior can be found using Bayes’ rule:

$$\begin{array}{cc}\hfill \pi \left(\mathit{k}\right|\mathit{d})& ={\displaystyle \frac{\pi \left(\mathit{d}\right|\mathit{k}\left)\pi \right(\mathit{k})}{\pi \left(\mathit{d}\right)}}\hfill \\ \hfill & \propto \pi \left(\mathit{d}\right|\mathit{k}\left)\pi \right(\mathit{k}).\hfill \end{array}$$

(23)

Here, $\pi \left(\mathit{d}\right|\mathit{k})$ is the likelihood function, which is the probability of receiving the data for a given permeability vector. If we assume that the noise $\mathit{e}$ and the parameters $\mathit{k}$ are independent, and given Equation (22), this likelihood is [14]:

$$\begin{array}{cc}\hfill \pi \left(\mathit{d}\right|\mathit{k})& \propto exp\left(-\frac{1}{2}{(\mathcal{G}\left(\mathit{k}\right)-\mathit{d})}^T{\mathbf{\Gamma }}_e^{-1}(\mathcal{G}\left(\mathit{k}\right)-\mathit{d})\right)\hfill \\ \hfill & =exp\left(-\frac{1}{2}\left|\right|{\mathit{L}}_e(\mathcal{G}\left(\mathit{k}\right)-\mathit{d}){\left|\right|}^2\right).\hfill \end{array}$$

(24)

where ${\mathbf{\Gamma }}_e^{-1}={\mathit{L}}_e^T{\mathit{L}}_e$.

On the other hand, $\pi \left(\mathit{k}\right)$ is the prior that encodes the prior beliefs of $\mathit{k}$. Given the parameterisation described in Section 3.1, it is now sensible to describe the prior using a multivariate Gaussian distribution with the mean ${\mathit{\mu }}_k$ and covariance ${\mathbf{\Gamma }}_k$. This gives the following equation for the prior probability density:

$$\begin{array}{cc}\hfill \pi \left(\mathit{k}\right)& \propto exp\left(-\frac{1}{2}{(\mathit{k}-{\mathit{\mu }}_k)}^T{\mathbf{\Gamma }}_k^{-1}(\mathit{k}-{\mathit{\mu }}_k)\right),\hfill \\ \hfill & =exp\left(-\frac{1}{2}\left|\right|{\mathit{L}}_k(\mathit{k}-{\mathit{\mu }}_k){\left|\right|}^2\right),\hfill \end{array}$$

(25)

where ${\mathbf{\Gamma }}_k^{-1}={\mathit{L}}_k^T{\mathit{L}}_k$. Combining Equations (23)–(25) gives the following equation for the posterior:

$$\begin{array}{cc}\hfill \pi \left(\mathit{k}\right|\mathit{d})& \propto exp\left(-\frac{1}{2}\left|\right|{\mathit{L}}_e{(\mathcal{G}\left(\mathit{k}\right)-\mathit{d})\left|\right|}^2-\frac{1}{2}\left|\right|{\mathit{L}}_k(\mathit{k}-{\mathit{\mu }}_k){\left|\right|}^2\right)\hfill \end{array}$$

(26)

3.3. Exploring the Posterior

Although Equation (26) is the posterior, we need some way of interpreting this distribution. The posterior could be analysed using Markov Chain Monte Carlo (MCMC) methods [8], but this would involve evaluating the posterior hundreds of thousands of times. Calculating the posterior is a computationally expensive operation because it requires calculating $\mathcal{G}\left(\mathit{k}\right)$ and running the RTM model. The RTM model involves solving the system of equations for pressure at each timestep, where the number of timesteps scales with the number of nodes. Other research in this area focuses on the development of less computationally expensive methods [18,19].

We use a Laplace approximation to the posterior, which requires several orders of magnitude less computational time. The Laplace approximation works by approximating the posterior about the maximum a posteriori (MAP) estimate. We then take a Gaussian approximation of this estimate. However, the calculation of the MAP estimate is not trivial and requires maximising the posterior. This can be simplified by noting the following.

$$\begin{array}{cc}\hfill {\mathit{k}}_{\mathrm{MAP}}& =\underset{\mathit{k}}{\mathrm{argmax}}\left(\pi \left(\mathit{k}\right|\mathit{d})\right)\hfill \\ \hfill & =\underset{\mathit{k}}{\mathrm{argmin}}\left(-2ln\pi \left(\mathit{k}\right|\mathit{d})\right)\hfill \\ \hfill & =\underset{\mathit{k}}{\mathrm{argmin}}\left(\left|\right|{\mathit{L}}_e{(\mathcal{G}\left(\mathit{k}\right)-\mathit{d})\left|\right|}^2+\left|\right|{\mathit{L}}_k(\mathit{k}-{\mathit{\mu }}_k){\left|\right|}^2\right).\hfill \end{array}$$

(27)

Once the MAP estimate has been found, the Laplace approximation to the posterior can be found by considering the Taylor expansion of $\mathcal{G}\left(\mathit{k}\right)$ about ${\mathit{k}}_{\mathrm{MAP}}$:

$$\begin{array}{cc}\hfill \mathcal{G}\left(\mathit{k}\right)& \approx \mathcal{G}\left({\mathit{k}}_{\mathrm{MAP}}\right)+\mathcal{J}(\mathit{k}-{\mathit{k}}_{\mathrm{MAP}}),\hfill \end{array}$$

(28)

where $\mathcal{J}$ is the Jacobian of the RTM model evaluated at the MAP estimate. This Jacobian considers the derivative of all outputs of $\mathcal{G}\left(\mathit{k}\right)$, with respect to each of the components of $\mathit{k}$:

$$\begin{array}{c}\hfill \mathcal{J}=\left[\begin{array}{cccc}{\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial {\mathit{k}}_1}}& {\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial {\mathit{k}}_2}}& \cdots & {\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial {\mathit{k}}_n}}\\ {\displaystyle \frac{\partial {\mathcal{G}}_2}{\partial {\mathit{k}}_1}}& {\displaystyle \frac{\partial {\mathcal{G}}_2}{\partial {\mathit{k}}_2}}& \cdots & {\displaystyle \frac{\partial {\mathcal{G}}_2}{\partial {\mathit{k}}_n}}\\ & ⋮& \ddots & \\ {\displaystyle \frac{\partial {\mathcal{G}}_m}{\partial {\mathit{k}}_1}}& {\displaystyle \frac{\partial {\mathcal{G}}_m}{\partial {\mathit{k}}_2}}& \cdots & {\displaystyle \frac{\partial {\mathcal{G}}_m}{\partial {\mathit{k}}_n}}\end{array}\right],\end{array}$$

(29)

where ${\mathcal{G}}_i$ is the ith component of the output vector from $\mathcal{G}$. Substituting this into Equation (26) gives an equation for the Laplace approximation of $\pi \left(\mathit{k}\right|\mathit{d})$, ${\pi }_{\mathrm{Laplace}}\left(\mathit{k}\right|\mathit{d})$:

$$\begin{array}{cc}\hfill {\pi }_{\mathrm{Laplace}}\left(\mathit{k}\right|\mathit{d})& \propto exp\left(-\frac{1}{2}\left|\right|{\mathit{L}}_e{(\mathcal{J}\mathit{k}-\mathit{z})\left|\right|}^2-\frac{1}{2}\left|\right|{\mathit{L}}_k(\mathit{k}-{\mathit{\mu }}_k){\left|\right|}^2\right),\hfill \end{array}$$

(30)

where $\mathit{z}=\mathit{d}+\mathcal{J}{\mathit{k}}_{\mathrm{MAP}}-\mathcal{G}\left({\mathit{k}}_{\mathrm{MAP}}\right)$, which does not depend on $\mathit{k}$. This probability density function describes a Gaussian distribution (that is, ${\pi }_{\mathrm{Laplace}}\sim \mathcal{N}({\mathit{k}}_{\mathrm{MAP}},{\mathbf{\Gamma }}_{\mathrm{post}})$), the covariance of which can be calculated as follows:

$$\begin{array}{cc}\hfill {\mathbf{\Gamma }}_{\mathrm{post}}& ={({\mathcal{J}}^T{\mathbf{\Gamma }}_e^{-1}\mathcal{J}+{\mathbf{\Gamma }}_k^{-1})}^{-1}.\hfill \end{array}$$

(31)

4. Numerical Examples

We now apply the outlined methodology to two example problems. For these, the data are generated by running the RTM model with a true permeability, which we denote by ${\mathit{k}}_{\mathrm{true}}$, which we are trying to estimate. Then, noise is added to these simulated pressures. It is assumed that the permeability in the bulk of the domain (${\mathit{K}}_{bulk}$) is homogeneous and that the permeability in the race-tracking regions is a scalar multiple of the bulk region.

We assume the data come in the form of pressure measurements taken at discrete sensors at predefined times. The pressure sensors are placed at optimal locations using the scheme described in the later section on Bayesian Optimal Experimental Design (Section 6). All pressure measurements are then normalised by dividing by the inlet pressure.

The programming language Julia was used to both execute the RTM simulations and to complete the Bayesian calculations in general. This allows for faster code execution compared to other interpreted languages, and we use the ForwardDiff package for automatic derivative calculation [31].

A Gauss–Newton algorithm [32] with random restarts is used to carry out the optimisation required to find the MAP estimate. We can reliably produce solutions with smaller residuals (i.e., a lower cost function described by Equation (27) than ${\mathit{k}}_{\mathrm{true}}$, with less than 10 restarts, for the small-scale models considered in this paper.

4.1. Model 1

The first investigated example is illustrated in Figure 2. Resin is injected into the mould from a 5 cm injection gate in the lower left corner of the domain at a constant pressure of $p_{\mathrm{in}}=400$ kPa. The only exit vent is at the top right corner of the domain, which is again fixed at 5 cm. There are three areas of potential race tracking, with different permeabilities at the top and bottom of the domain, as well as around a central block. In reality, there will be potential for race tracking along the vertical walls, but this is ignored here for simplicity.

As there are three race-tracking regions, the permeability vector is

$$\begin{array}{c}\hfill \mathit{k}={[c_1,c_2,c_3,{\alpha }_1,{\alpha }_2,{\alpha }_3]}^T,\end{array}$$

(32)

where:

$$\begin{array}{c}\hfill {\mathit{K}}_i=exp\left({\alpha }_i\right){\mathit{K}}_{\mathrm{bulk}},\end{array}$$

(33)

and ${\mathit{K}}_i$ represents the permeability tensor for region $i=1,2,3$, as shown in Figure 2.

Here, ${\mathit{k}}_{\mathrm{true}}$ was arbitrarily defined as follows:

$$\begin{array}{c}\hfill {\mathit{k}}_{\mathrm{true}}={[0,0,0,ln\left(10.0\right),ln\left(1.0\right),ln\left(50.0\right)]}^T.\end{array}$$

(34)

This represents a scenario where there is no race tracking along the bottom, 10 times higher permeability on the top, and 50 times higher permeability around the centre of the domain. All permeability values here have units of ${10}^{-11}{\mathrm{m}}^2$.

Error was added to the output of this model with the distribution $\mathcal{N}(\mathbf{0},{\mathbf{\Gamma }}_e)$. For this experiment, we set

$$\begin{array}{c}\hfill {\mathbf{\Gamma }}_e={\sigma }_e^2\mathit{I},{\sigma }_e=0.01.\end{array}$$

(35)

As with the error distribution, the prior for this experiment is arbitrary and can be adjusted to the beliefs of the implementer but should, in general, be created broadly so as not to introduce bias into the prediction. This experiment used the following prior, which roughly approximated our prior beliefs in the race-tracking and permeability values.

$$\begin{array}{cc}\hfill {\mathit{\mu }}_k& ={[0,0,-0.144,2,2,2]}^T,\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill {\mathbf{\Gamma }}_k& =\mathrm{diag}\left(\right[0.25,1,0.25,2,2,2\left]\right).\hfill \end{array}$$

(37)

The prior mean represents a permeability field that it is slightly more permeable in the x direction than in the y direction. Specifically, it represents the following tensors:

$$\begin{array}{c}\hfill {\mathit{K}}_{\mathrm{bulk}}=\left[\begin{array}{cc}1& 0\\ 0& \frac{3}{4}\end{array}\right],{\mathit{K}}_1,{\mathit{K}}_2,{\mathit{K}}_3=e^2{\mathit{K}}_{\mathrm{bulk}}\end{array}$$

(38)

4.2. Model 2

This example considers a more complicated scenario in which there are ten race-tracking regions. This time, we let resin exit through the entire right-hand boundary. The reasons for this are discussed in Section 7.1. There are now two central blocks, each with a different race-tracking strength along each of its four edges. The domain is illustrated in Figure 3.

The associated permeability vector is

$$\begin{array}{c}\hfill \mathit{k}={[c_1,c_2,c_3,{\alpha }_1,{\alpha }_2,\dots ,{\alpha }_{10}]}^T,\end{array}$$

(39)

where ${\alpha }_1,{\alpha }_2,\dots ,{\alpha }_{10}$ are defined in Equation (21) and their respective regions are labelled in Figure 3; ${\mathit{k}}_{\mathrm{true}}$ is taken as a random sample from the prior, which is defined as similar to Model 1, specifically:

$$\begin{array}{cc}\hfill {\mathit{\mu }}_k& ={[0,0,-0.144,2,2,\dots ,2]}^T,\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill {\mathbf{\Gamma }}_k& =\mathrm{diag}\left(\right[0.25,1,0.25,2,2,\dots ,2\left]\right).\hfill \end{array}$$

(41)

We use the same error distribution as in Model 1, given in Equation (35).

5. Results

5.1. Model 1

The posterior and prior for Model 1 are compared in Figure 4 when only two pressure sensors are used. We take measurements every 1000 s, for up to 15,000 s. When ${\mathit{k}}_{\mathrm{true}}$ is used, this gives a fill time of 10,500 s, which means that we measure both during filling and after the resin has reached the vent—a steady state pressure distribution.

The posterior encompasses the truth for all parameters as expected. The posteriors are generally broad due to the low amount of data used to develop the posterior, particularly in parameters such as ${\alpha }_1$, ${\alpha }_2$, and ${\alpha }_3$, which are perhaps less influential on the pressure measurements themselves. It is important to note that although Figure 4, Figure 5 and Figure 6 show each parameter’s marginal posterior individually, there is a multivariate correlation between the parameters developed by the posterior that is not visualised.

Figure 5 shows how the posterior varies when 10 pressure sensors are used. With the increase in information provided by the 10 sensors, the marginal posteriors have less uncertainty. However, the caveats of using the Laplace approximation are evident here, particularly for the posterior for ${\alpha }_2$, where the truth is not encompassed by the posterior. Although we can be sure that the unapproximated posterior would cover these truths, the Laplace approximation can sometimes be inaccurate. In addition, the posterior with 10 sensors has greater uncertainty than with 2 sensors for this parameter. One key reason for these inaccuracies in the Laplace approximation is discussed in Section 7.1.

These results are quantitatively summarised in

Table 1. Figure 4 and Figure 5 summarised in a table. Confidence intervals are 99% based on 3 standard deviations.

	$c_1$	$c_2$	$c_3$
Truth	$-0.141$	$0.531$	$-1.758$
Prior	$0.000\pm 1.500$	$0.000\pm 2.121$	$-0.144\pm 1.500$
Posterior (2 Sensor)	$-0.128\pm 0.034$	$0.537\pm 0.036$	$-1.731\pm 0.076$
Posterior (10 Sensor)	$-0.127\pm 0.020$	$0.533\pm 0.015$	$-1.747\pm 0.034$
	${\alpha }_1$	${\alpha }_2$	${\alpha }_3$
Truth	$3.737$	$2.824$	$2.189$
Prior	$2.000\pm 4.243$	$2.000\pm 4.243$	$2.000\pm 4.243$
Posterior (2 Sensor)	$3.751\pm 0.421$	$2.612\pm 0.166$	$2.093\pm 0.563$
Posterior (10 Sensor)	$3.653\pm 0.254$	$2.605\pm 0.196$	$2.067\pm 0.379$

, which shows how the prior, posterior, and truth compare when a 99% confidence interval is developed from the prior and posterior.

5.2. Model 2

The marginal priors and posteriors for the more complex Model 2 are compared in Figure 6. Because there are more parameters to estimate, we now take more data. We consider only the 10-sensor setup and take measurements every 100 s for a total of 10,000 s. We take measurements in a shorter time period because the mould takes only 7800 s to fill with ${\mathit{k}}_{\mathrm{true}}$.

The uncertainty in influential parameters such as the bulk permeability and the top and bottom race-tracking regions is very low because these parameters have a greater impact on domain pressure. The uncertainty in the less significant regions around the central blocks is much larger. The truth is encompassed by the posterior for all parameters, as expected.

6. Bayesian Optimal Experimental Design

This paper has so far shown that a lot of information about race tracking and permeability in the preform could be determined with even two sensors, taking 15 discrete measurements each. This section investigates where to optimally place sensors in order to bring posterior uncertainty to a minimum and how this varies when different amounts of sensors are used. The goal here is to allow practitioners to extract the maximum amount of information using minimal measuring equipment.

6.1. Method

We have $N\in \mathbb{N}$ sensors; we wish to place them within the mould so that when considering all potential permeability scenarios, the variance in parameters is minimised. One way to achieve this is to determine the covariance matrix for a variety of permeability scenarios with a sensor setup S. Then, we aim to minimise the sum of these covariance traces, which is known as A-optimality [33,34,35,36,37,38,39,40]. In essence, this is minimising the following cost function:

$$\begin{array}{c}\hfill \underset{S}{min}\sum _{i=1}^M\mathrm{tr}\left({\mathbf{\Gamma }}_{\mathrm{post}}^S\left({\mathit{k}}_i\right)\right),\end{array}$$

(42)

where M is the number of permeability scenarios considered and ${\mathbf{\Gamma }}_{\mathrm{post}}^S\left({\mathit{k}}_i\right)$ is the posterior covariance for the permeability scenario ${\mathit{k}}_i$ and sensor design S. A-optimality is appealing for this problem, where we could be interested in using the MAP permeability estimate to predict future resin flow and want minimal uncertainty in each individual estimate. Common alternatives, such as D-optimality or E-optimality [33], would consider the off-diagonal elements of the posterior covariance, which are not as relevant for this situation.

To simplify this problem, potential sensor locations will be limited to discrete locations spaced 0.025 m apart throughout the domain, resulting in 697 potential locations. When considering that the possible number of locations for N sensors is ${697}^N$, we will limit the optimisation to follow a greedy decision process [36,39]. That is, the optimal location for the first sensor will be found, and then the second sensor will be optimally placed, given that the first is fixed. The process will be repeated for all N sensors. Sensors are also prohibited from occupying the same position. An advantage of the greedy algorithm is that the optimal design for N sensors gives us the optimal design for all $n<N$ sensors.

Despite reducing the computation time significantly with a greedy algorithm, the computation for solving this problem is still extremely intensive. Determining the MAP estimate by optimisation took 1–2 h on a standard desktop computer (12th Gen Intel Core i7-12700 CPU, no GPU acceleration). For more complex geometries and finer meshes, this time can significantly increase. Solving this inverse problem M times for each sensor location will take $697M$ times longer for every sensor. Although this could be parallelised, this would still take too long (potentially years) to receive a solution in a reasonable amount of time. As a computationally feasible approximation [39], we use:

$$\begin{array}{c}\hfill {\mathbf{\Gamma }}_{\mathrm{post}}^{\prime }={({\mathcal{J}}_{\mathrm{true}}^T{\mathbf{\Gamma }}_e^{-1}{\mathcal{J}}_{\mathrm{true}}+{\mathbf{\Gamma }}_k^{-1})}^{-1},\end{array}$$

(43)

where ${\mathcal{J}}_{\mathrm{true}}^T$ is the Jacobian of the forward model evaluated at ${\mathit{k}}_{\mathrm{true}}$. This is in contrast to Equation (31), where now we use the Jacobian evaluated at ${\mathit{k}}_{\mathrm{true}}$ for the current sample, instead of evaluating it at the MAP estimate. Although this is an approximation, the optimal sensor locations using this approximation should not vary significantly from the true optima [39] and bring the problem to a solvable time frame (hours). Note that although this approximation does not make sense to use in a standard Bayesian inference scenario, it is sufficient for the purposes of judging a given sensor design using Equation (42).

Computation can be further shortened by noting that we can precompute ${\mathcal{J}}_{\mathrm{true}}$ for all potential sensor locations. Then, for a given sensor design S, we can simply extract the rows of the Jacobian that correspond to the relevant outputs, where the sensors have been placed. This can be computed for all M permeability scenarios.

6.2. Results

The optimal sensor locations for Model 1 are shown in Figure 7, calibrated using $M=200$ different permeability scenarios. The choice of M was decided by investigating the robustness of the optimal locations trained with various amounts of data. This is omitted here for clarity but is contained in Appendix C, along with a brief discussion on the robustness of the design. The sensors tend to be distributed around the race-tracking areas. For the top and bottom race-tracking areas, it seems to be beneficial to have sensors at the start and end of these regions. There are some sensors (most notably sensor 2) that are positioned to gain pressure measurements at the inlet, where the large pressures make the influence of noise less confounding.

Figure 8 illustrates how the objective function in Equation (42) changes with an increasing number of sensors, for Model 1. The results for (the more complex) Model 2 are given in Appendix B. The figure shows a monotonic decrease in objective function as more sensors are added, regardless of location, as expected. For the random design, there are points where adding a sensor gives no change in the objective. This is because the sensor is placed in a location with no information—for example, inside the central block. The optimal sensor locations lead to a significant reduction in uncertainty compared to using random sensor locations, with roughly five times lower variance in the permeability estimates. However, this difference would be expected to converge to zero as the number of sensors increases to the maximum [33,34,35,36,39,40]. It is also clear in Figure 8 that the largest variances come from the race-tracking parameters ${\alpha }_1$, ${\alpha }_2$, and ${\alpha }_3$ rather than the bulk parameters $c_1$, $c_2$, and $c_3$. This is in agreement with the posterior variances shown in Section 5. The eighth sensor appears to have greatly reduced the uncertainty in $c_1,c_2$, and $c_3$. However, the logarithmic scale of the graph emphasises this drop, and the sum of variances for all parameters still decreases, as seen in the upper graph (red line).

The MAP estimate for each parameter is shown in Figure 9, along with 99% confidence intervals. These intervals are based on the Laplace approximation to the posterior about the MAP estimate (Equation (31)), rather than the approximation defined in Equation (43). The MAP estimate tends to converge to the truth as more sensors are added, as expected. The variance tends to decrease as more sensors are added, as in Figure 8, but not monotonically. When calculating the Laplace approximation about the MAP estimate, there is no guarantee of strictly decreasing variances because the MAP estimate changes. There are some instances where the truth is not contained within the 99% confidence interval. This is due to a discrepancy between our Laplace approximation and the unapproximated posterior. One key reason for this discrepancy in this scenario is discussed in Section 7.1.

7. Discussion and Conclusions

We end this paper by discussing the limitations of this methodology and potential future work.

7.1. Non-Differentiability

The Laplace approximation of the posterior, along with all gradient-based optimisation techniques that are used to find the MAP estimate, relies on the RTM process being differentiable. In this section, we discuss how this assumption only holds under specific experiments.

When pressure measurements are extracted at fixed, equal time intervals within a predefined period, the pressure data from a given sensor can be split into three distinct phases:

• Before the resin reaches the sensor, the sensor gives 0 pressure (with noise).
• Once the resin reaches the sensor, the pressure increases from 0 in a decaying, quasi-logarithmic fashion.
• When the resin has almost filled the mould, the limited size of the exit vent means that the boundary condition on the resin changes from a Dirichlet boundary to a zero-flow Neumann boundary for the large majority of the boundary. This rapid change in the boundary condition causes a large and sudden increase in pressure throughout the mould. This change is sometimes so quick that it is ‘almost discontinuous’.

This is illustrated in Figure 10. The three phases are clearly visible here. The ‘almost discontinuity’ that appears in the posterior comes from the following: if, in the example shown in Figure 10, there was a measurement taken at t = 23,000 s, then a slight increase in any of the components in the permeability vector $\mathit{k}$ from either the bulk permeability or from the race tracking will cause the pressure measured here to rapidly change from roughly 0.2 to roughly 0.35. This effect is more prominent in models with a small exit vent, where there is a large change in boundary conditions as the resin reaches the end. The scenario in Figure 10 is one with a very small but realistic (5 mm wide) exit vent—smaller than what is used in the numerical examples in this paper. However, the effect still exists, to a lesser extent, for larger vents.

Figure 11 shows the true posterior (without using the Laplace approximation) for a much simpler formulation of the problem, where the only parameter to predict is $K=K_{xx}=K_{yy}$ for a homogeneous domain, with $K_{xy}=0$. Again, a central block is included, as in Section 4.1. The true value for this parameter is one. The inverse problem is solved with four different outlet boundary conditions: a single node, the top 0.1 m of the right boundary, the top 0.25 m of the right boundary, and the entire right boundary. As the outlet vent becomes smaller, a nondifferentiable point develops around $K=1.02$. $K=1.02$ is a point that exhibits the behaviour discussed above, where increasing this value a small amount causes one measurement to be taken after the resin has reached the end of the mould, instead of before. This point still exists even when the outlet is as large as 0.25 m. The roughness of the posterior shape aside from this discontinuity should be ignored, it is the product of using a coarse mesh in order to develop the true posterior within an acceptable time frame.

The non-differentiability causes two key limitations to our results:

• It increases the discrepancy between the Laplace approximation and the true posterior, as seen in Figure 11 and in the results for Model 1, where the true values were outside of the Laplace approximation.
• Finding the MAP estimate through optimisation techniques becomes much more difficult, as the function becomes nonsmooth.

For Model 1, the nondifferentiable points are far away from the MAP estimate, so there is no noticeable effect on the posterior. However, the optimisation to find the MAP estimate is long and difficult due to these nondifferentiable regions. For Model 2, we use a large exit vent to reduce the effects of this problem, speeding up the optimisation, and ensuring that the posterior is not affected.

For the numerical examples in this paper, we can minimise this issue by measuring only while the mould is filling. In a potential future case where we wish to estimate race tracking during the Resin Injection stage, it would not be useful to estimate after the mould has already been filled. However, when we consider more complicated geometries, there could be scenarios in which the resin reaches the end of a channel but the mould is not yet filled. So, the nondifferentiable issues mainly arise at the end of filling for our examples; this is not the case in general for all geometries.

Even in cases where our exact methodology is not used, the non-differentiability should be noted in all research where a Jacobian or gradient of the RTM process is considered, as differentiability is often assumed.

7.2. Conclusions

By representing the permeability throughout a preform as a homogeneous but anisotropic property, except for regions of race tracking, we are able to predict the race-tracking effect as well as the bulk permeability to a high level of certainty. A Bayesian inversion algorithm was used to achieve this and was demonstrated with two example scenarios that showed that the precision of the estimates increased with fewer race-tracking regions and when more sensors were used. Although these models are simple and only two-dimensional, the approach could be extended to three-dimensional, with many race-tracking regions.

In this paper, we only estimated the permeability once the mould had been filled and measurements had been collected throughout the filling process. However, in the future, this methodology could be applied through different stages of mould filling. This would allow researchers to estimate race tracking and then predict the flow front propagation pattern for the remainder of the filling. It would then be possible to take some control actions if we predict that the resin will propagate in an unfavourable way.

The uncertainty in the computed posterior distributions developed here considers only measurement error from sensors, where the measurement from the sensor is made noisy by the addition of randomly and independently sampled error. However, there are approaches that deal with another significant type of error—the approximation errors that arise due to the use of an approximate model, which differs from the true RTM process. There are approaches to deal with this so-called Bayesian approximation error [41,42,43] that could be investigated and integrated into the current framework in the future. In particular, this paper made several simplifying assumptions to the model, such as a piecewise homogeneous permeability and fixed porosity, which add to this Bayesian approximation error.

We have also used the Bayesian formulation to develop a methodology for optimising sensor placement within a mould to extract maximum insight. This investigation also revealed the decaying value of adding more sensors. The optimal sensor configuration showed significant decreases in uncertainty compared to random locations. In the future, this framework could be combined with an online permeability estimation tool for increased predictive power. The results show that using a large number of sensors with a tedious, disruptive, and potentially expensive setup is not necessary—very similar accuracy can be achieved with between five and seven sensors, provided that they are in optimal locations. The investigation has also shown that even two sensors can provide reasonably accurate estimates for race tracking for simple scenarios.