Multiphase flows (MPFs) include motions of various phases and deformations of interfaces formed by the phases. Interactions among the phases are highly non-linear and complex. Although it is possible to predict MPFs by collecting their historical information, such as the flow rate, temperature, pressure, and positions of the phases, the mechanism is so complex that it is still difficult to build a general mathematical model to describe different kinds of MPFs or to solve those models efficiently and accurately. Predicting MPFs becomes a popular topic and many researchers have focused on it because they benefit a wide range of fields, e.g., the chemical industry, environmental protection, life sciences, and industrial process [1][2][3][4][5][6] .
In recent years, the development of neural network (NN) technology 7,8 promotes time series modeling, which also contributes to predicting MPFs. With the use of selforganizing NNs, Mi et al. 9 simulated the vertical bubbly, slug, churn, and annular flows for both the impedance of the multiphase admixture and the flow patterns. Based on the given experimental and theoretical velocity data, Valero and Bung 10 identified the characteristic shapes of the fluid phases and explored the behaviors of the air-water flow.
Rashid et al. 11 developed a radial basis function NN to predict MPFs under critical conditions and discovered the relationships between the upstream pressure and gas-liquid ratio, and between the choke bin size and liquid flow rate. Serra et al. 12 proposed a randomized hough transform with an NN to predict the void fraction with given image samples, and their results are consistent with the actual void fraction values in natural circulation-based systems of nuclear power plants. Xuan et al. 13 proposed a data-driven NN model to predict the separation-induced transition. Lino et al. 14 proposed a graph neural network to extrapolate the time evolution of the unsteady Eulerian fluid dynamics. More works on MPFs simulation with the help of NNs were summarized in Yan s 15 and Karniadakis s 16 reviews.
To model MPFs, parameters inside NN will be updated until flow characteristics predicted by NN are similar to or even the same as the given samples. Some popular optimization algorithms [17][18][19] can be applied to update the parameters inside NN. However, how to incorporate physical laws when optimizing the parameter space of NNs so that the predictions strictly meet physical constraints is still a problem. Consequently, NN can easily produce nonphysical results, especially when the number of samples is insufficient to train a complex NN with a large parameter space. Thus, recent studies are considering including physical constraints to improve the performance of NNs in the sense that not only available data are more effectively used, but also the burden of data acquisition is alleviated. The physical constraints can be enforced either implicitly or explicitly. For implicit constraints, penalties are added to the loss function and one obtains the so-called physics-informed neural networks (PINNs) 20 . For explicit constraints, predictions from NN are modified to satisfy the physical constraints exactly [21][22][23][24] .
Physics-informed neural networks (PINNs) aim to include the physical constraints implied in the available data, by adding penalties in the loss function 20 . This helps to reduce the hypothesis space of parameters, and thus requires fewer observations. Moreover, an NN model with desired physical properties is more acceptable and explainable in scientific problems because it has a better generalizable property in out-of-sample scenarios 25 .
Pang et al. 26 developed fractional physics-informed neural networks (fPINNs) to encode fractional-order partial differential equations (PDEs). Yang and Perdikaris 27 proposed probabilistic PINNs to approximate arbitrary conditional probability densities and to generate samples similar to the one given by PDEs. Dwivedi and Srinivasan 28 considered the advantages of PINNs and extreme learning machine, and proposed Physics Informed Extreme Learning Machine (PIELM) to speed up the network training process.
Leake and Mortari 29 transformed PDEs into unconstrained problems, and then solved them by training NNs. Goswami et al. 30 proposed a PINN algorithm to solve brittle fracture problems and optimize its loss by decreasing the variational energy of the system. Schiassi et al. 31 proposed a novel PINNs framework to satisfy DE constraints analytically but with unconstrained parameters in NNs, which was also applied to thermal creep flow 32 and Poiseuille flow simulations 33 . Qiu et al. 34 proposed PINNs for incompressible phase field method, which can capture the dynamic behavior precisely. Many existing studies have indicated that an NN with physical constraints contributes to improving prediction accuracy [35][36][37][38][39] , discovering PDEs and governing equations 40,41 , modeling inverse problems 42,43 , and solving uncertainty quantification 27,44 . Some other studies also focus on considering a general deep learning framework to learn diverse continuous non-linear operators [45][46][47] . Therefore, adding penalties from the physical constraints in the loss function becomes more and more acceptable and interpretable when developing NNs for physical problems.
In the present study, we proposed a PCNN to predict the temporal dynamics of MPFs.
Our method partly remedies a defect of PINNs that the physical constraints are not enforced exactly. Specifically to MPFs, the defect of PINNs may produce local voids, overfilling, mass loss, negative volume fractions, or negative density of the fluid mixture. In the present method, three physical constraints 24 on the order parameters (OPs), which locate individual phases, are exactly satisfied by correcting the prediction from the NNs, instead of adding penalties to the loss function. The first constraint is called the summation constraint, which requires the summation of the volume fractions to be unity. The second one is called the conservation constraint, which requires the mass (or volume) of each phase to be conserved if there are no external inputs. The last one is called the boundedness constraint, which requires the volume fractions produced by OPs to be between zero and one. We demonstrate in the present study that the nonphysical behaviors due to failures of satisfying the three physical constraints will be propagated into the training process and gradually accumulated. As a result, large errors are introduced as time goes on. Once we obtain the prediction of OPs, we proceed to compute the density of the fluid mixture, which later becomes the input to the velocity prediction with the constraint of momentum conservation.
The structure of the proposed PCNN for MPFs is as follows. The OPs at the new time level are first predicted from the NN model with the OPs at previous time steps. The time-dependent dynamics of the OPs are encoded by a recurrent neural network and decoded predicted OPs for the next time step. The predicted OPs are then corrected by the multiphase consistent and conservative boundedness mapping (MCBOM) algorithm 24 to strictly satisfy the four physical constraints mentioned, and the corrected OPs are the final prediction of the OPs at the new time level. As a result, one can update the density of the fluid mixture using the OPs at the new time level. Finally, the velocity at the new time level is predicted from the velocity at previous time levels and the density of the fluid mixture, using the NN model where the momentum conservation is included as a penalty in the loss function. Considering the momentum conservation as a penalty can shrink the model parameter space towards initial momentum with respect to the maximum likelihood estimates.
The remaining sections are structured as follows. In Section II, the problem of interest is defined. Section III introduces the proposed PCNN for MPFs. Then, in Section IV, we verify the feasibility and practicability of the proposed model through two applications.
Finally, conclusions are drawn in Section V.
In this section, we define the MPFs to be modeled, as well as their physical constraints and analyze the MPFs with time series.
In the present study, the incompressible MPFs are considered, where there are N (N 1) immiscible fluid phases. Each phase has a fixed density and viscosity, denoted by ρ p and µ p , respectively, for phase p, and every two-phase has a surface tension, denoted as σ p,q for phases p and q. Locations of the phases are determined from the order parameters {φ p } N p=1 , such that φ p = 1 represents pure phase p while phase p is absent at φ p = -1. Three physical constraints need to be strictly satisfied by the order parameters. The first one is called the summation constraint, i.e.,
The summation constraint requires that the summation of the volume fractions of all the phases is one so that local voids or overfilling is not produced.
The second constraint is the mass (or volume) conservation, i.e.,
where
is the domain of interest, ∂Ω is the domain boundary, n is the unit outward normal at the domain boundary, and v is the velocity. If there is no normal velocity at the domain boundary, the integral of φ p (1 p N) over the domain of interest does not change. As a result, the volume and mass of each phase, i.e., Ω 1+φ p 2 dΩ and Ω ρ p 1+φ p 2 dΩ, respectively, 1 p N, do not change with time either. The third constraint is the boundedness constraint, which requires
One can physically explain {φ p } N p=1 as volume fraction contrasts only when they are in [-1, 1]. This corresponds to the volume fractions
Here, the volume fraction contrast of Phase p is defined as φ p = C p -∑ N q=1,q =p C q , where {C p } N p=1 is the volume fractions of the phases.
The density of the fluid mixture is
From Eq. ( 4), if φ p is outside [-1, 1], the density of the fluid mixture can become negative.
Then, (ρv) is the momentum and should be conserved if there are no external forces at the domain boundary.
Time series is a sequence of data set arranged in chronological order 48,49 . The time series tool is applied to predict the dynamics of MPFs, or specifically, the order parameters and flow velocities of MPFs. Given the time series of the observed order parameters {φ i } m i=1 and flow velocity {v i } m i=1 , the goal is to predict the unobserved data in the future, i.e., {φ i } n i=m+1 and {v i } n i=m+1 . Here, m and n are the numbers of the observed time steps and total (including both observed and unobserved) time steps, respectively. The subscripts i is timestamp indexes for the data. To achieve the goal, we specify two projections, where one f 1 (•) from φ k to φ k+1 , and another f 2 (•) from v k to v k+1 such that:
where
Here ε 1 and ε 2 are both bias vectors {ε i } t+i i
To approximation f (•), we consider neural networks f (•):
where
) are parameters to minimize the difference between ground truth and prediction. In the next section, we will elaborate the way to design our neural networks and to incorporate the physical laws inherent in MPFs.
This section will elaborate on the proposed PCNN for multiphase flow prediction, including the NN model to predict the order parameters in Section III A, the boundedness mapping (MCBOM) from Huang et al. 24 to correct the order parameters in Section III B, and the second NN model to predict the velocity in Section III C.
The following operations are applied to φ p (1 p N) independently, and we skip the superscript for convenience. Given a set of order parameters {φ i } m i=1 as the input, a NN is developed to predict the order parameters {φ i } n i=m+1 . The order parameters are firstly encoded by an LSTM, whose workflow is shown in Figure 1(a). LSTM is effective to solve long-term dependent problems because it introduces the forget gate f i in each LSTM cell, see Figure 1(a), to control the circulation and loss of the sequential data. For each LSTM cell, the input of an LSTM cell includes the output of its previous cell, and all LSTM cells share the same structure and parameters.
LSTM cells are then connected by the hidden states and cell states to form a complete LSTM network. In the following content, we will formulate the structure of an LSTM cell and then explain how LSTM cells connect with subsequent parts. Inside the i th LSTM cell (1 i m -1), a forget gate f i ∈ [0, 1] is defined as a Sigmoid unit 50 ,
where S(•) is the sigmoid function, h i-1 is the hidden state in the previous LSTM cell.
Here, b, U, and W are the bias, input weight, and cycle weight, respectively, and their subscript f represents the forget gate. The same format is used in Eqs.( 8)- (11).
Then, an input gate I i , i.e.,
decides which information needs to be updated. To obtain the candidate value C i , a hypertangent activation function is applied,
Then, we obtain the current cell state C i , incorporating the previous cell state C i-1 , the candidate value C i in Eq. ( 9), the forget gate f i in Eq. ( 7), and the input gate I i in Eq. ( 8),
i.e.,
The cell states can be regarded as the memory of LSTM, which connects the previous, current, and future LSTM cells. Next, we move to an output gate o i , which decides the general state of the next hidden state. Here we use the sigmoid function as the activation function and project φ i to the output gate o i by
where h i-1 is the previous-step hidden state. Finally, hidden state h i is computed from the current cell state C i in Eq. ( 10) and the output gate o i in Eq. ( 11) with
The output range for sigmoid is [0, 1], while tanh is [-1, 1]. The input gate, forget gate, and output gate use sigmoid because they control the quantities between 0 and 1, where 0 means no flow and 1 means complete flow. The estimated means of cell state C i and hidden state h i are expected to be 0, so we consider tanh function. We have completed the operations in an LSTM cell. Once the hidden state h i is computed, we send current hidden state h i to a quasi encoder-decoder model. The representation r i is defined as a projection of hidden state h i and order parameter φ i :
where H(•) is an affine transformation using a non-linear activation (usually this is a sigmoid function):
Here we first concatenate vectors φ i and h i as a new vector, transpose the new vector, multiply it by the parameter matrix A 1 , and then add the offset parameter b 1 . We will do the same affine transformation in Eq. ( 16) but with different parameters. Then a commutative operation maps all the representations {r i } m-1 i=1 into a one-dimensional element R i with the following rules:
where ⊕ here is the mean operation. Then we use another affine transformation with a non-linear activation to map the elements R i into φ i+1 :
To update the parameters in the NN model, we defined a loss function to calculate the difference between the predicted φ and their observations. Given the observed order parameters φ i+1 and the prediction φ i+1 , the loss function L φ (θ 1 ) is defined as:
where N C here is the number of predicted order parameters in the domain Ω, and n C denotes the index of the grid cells. In practice, we minimize L φ (θ 1 ) by Adam optimizer 19 (see Algorithm 1).
After predicting {φ p k } N p=1 (m + 1 k n) from all the previous data before timestamp k by a NN model in Section III A, there is no guarantee that {φ p k } N p=1 satisfy the summation constraint, the mass conservation, and the boundedness constraint, which are the physical constraints on the order parameters, as discussed in Section II. Even though errors due to violating those constraints are small in a time step, they will accumulate and become significant as the computation proceeds. The multiphase consistent and conservative boundedness mapping (MCBOM) algorithm proposed in 24 Given a set of order parameters {φ p } N p=1 and a set of scalars {Φ p } N p=1 that represents total amounts of {φ p } N p=1 inside the domain of interest, the boundedness mapping includes three steps. The first one is the clipping step:
The second one is the re-scaling step
The last step is the conservation step
where W p,q is the weight function
and {B p } N p=1 are obtained from solving the linear system
In practice, the integral is approximated by the mid-point rule.
MCBOM corrects the prediction given by NN so that the properties of the order parameters are strictly satisfied. The resulting {φ p b } N p=1 has the following properties:
The cost of the MCBOM algorithm is negligible compared to the neural networks because it is an pointwise algebraic operator. The major cost of MCBOM comes from inverting an N by N symmetric and diagonal-dominant matrix, where N is the number of the phases that is usually much smaller than the parameter space of the NNs. Interested readers can refer to Huang et al. 24 for more details.
Once the order parameters are predicted from the first NN model in Section III A and corrected by MCBOM in Section III B, the density of the mixture is computed from Eq. ( 4). Thanks to satisfying the boundedness constraint of the order parameters, the density of the fluid mixture will stay in their physical interval as well, no matter how large the density ratio of the problems could be. The remaining task is to predict the velocity {v j } n j=m+1 given the observations {v i } m i=1 . This is accomplished by following the second NN model, whose structure is the same as the one in Section III A, but the input here becomes the velocity rather than the order parameters, i.e.,
Here (•) denotes the variables and parameters in the second NN model for velocity prediction.
The loss function of NN additionally considers the momentum conservation of the MPFs, which is different from the model in Section III A. The momentum of the multi-phase flow is
where M i is the total momentum of the multiphase flow at timestamp i, ∆Ω is the discrete volume. Recall that the density of the fluid mixture at timestamp i, i.e., ρ i , is computed from Eq. ( 4) once { φ p i } N p=1 are available. Finally, the loss function of NN for velocity prediction is
where λ is the Lagrange multiplier to enforce the momentum conservation. Here, ρ i+1 and v i+1 form a physics-informed constraint and penalize a nonphysical solution. Right now, the loss function in Eq. ( 26) not only includes the part to minimize the difference between the prediction and the ground truth, but also has the penalty due to violating the momentum conservation. A flowchart of the proposed PCNN in Sections III A-III C is shown in Figure 2. The methods to update the NN parameters and predictions are elaborated in Algorithm 1. Including the physical constraints help to alleviate the nonphysical behaviors produced simply by using the traditional NNs, and those effects will then be shown in the next section.
1 FC Layer FC Layer FC Layer FC Layer FC Layer Operation FC Layer LSTM …… …… …… …… LSTM LSTM Operation Operation FC Layer FC Layer FC Layer Operation …… Operation Operation 2 v FC Layer FC Layer FC Layer Algorithm 1 Update strategy using Adam optimization 19 .
Require: Decay rate: ρ φ and ρ v , momentum coefficient α.;
Require: Learning rate for order parameters and velocities: η, η .
Require: Exponential decay rates:
Require: Small constant for stabilization: . 1: Initialize time stamp: i = 1.
2: Initialize the first and second moment: s = 0, r = 0, s = 0, and r = 0. 3: while stopping criterion not met do 4:
6:
Update prediction:
Update loss: L φ (θ 1 ) ← 1
Eq. ( 17)
Compute gradient: g ← g -∇L φ (θ 1 ).
Update the first moment estimate:
Update the second moment estimate:
Correct the second moment estimation: r ← r 1-β 2 .
Update parameters: θ ← θη ŝ √ r+ .
Update loss:
Eq. ( 26).
Compute gradient: g ← g -∇L v (θ 2 ).
Update the first moment estimate: s ← β 1 • s + (1β 1 )g .
Update the second moment estimate: r ← β 2 • r + (1β 2 )g g .
Correct the first moment estimation: ŝ ← s 1-β 1 .
Correct the second moment estimation: r ← r 1-β 2 .
Update parameters: θ ← θη ŝ √ r + .
Update time step: i ← i + 1.
We will elaborate on two experiments here to demonstrate the effectiveness of the proposed PCNN. The first experiment is the horizontal shear layer problem and the other one is the dam break problem. The data sets are generated from the consistent and conservative phase-field method for incompressible MPFs 23 , including the bound-edness mapping algorithm 24 that is also introduced in Section III B. The sequential data has been spatially discretized, and the discrete data set is summarized in Table I, including the number of phases, simulation duration, number of time steps, and number of grids of the simulation. The dataset used in the present study is generated by the physical multiphase flow model in Huang et al. 24 which is a set of partial differential equations (PDEs), and the numerical methods therein. On the other hand, the present method is data-driven without solving any PDEs. The given data is on a uniform grid. The issue of using a random arrangement is that the total phasic mass and momentum in the entire domain, computed from Eqs. ( 22) and ( 25), are probably not accurate. If the given data on a random arrangement does not have that issue, the present model will still work.
Table I: A general summary of the data set. Dataset # phases # time steps Total time Grid size Horizontal shear layer 3 160 2.0 [128 × 128] Dam break 3 200 10.0 [512 × 128]
Horizontal shear flows happen under a variety of locations 52 , e.g., along coastal fronts, at the edges of energetic ocean currents, in oceans, and atmospheric currents around the terrain. Studies of the horizontal shear flows and predictions of their propagation trends can be widely applied to the hemodynamic analysis 53 , ocean circulation 54 , etc.
In this case, the numerical domain of the problem is in a [1 × 1] box with the periodic boundary condition. The domain is partitioned into 128 × 128 non-overlapping, uniform cells. The time step is set to be ∆t = 0.0125. The density of the fluids here are: ρ 1 = 50.0, ρ 2 = 10.0 and ρ 3 = 1.0. At t=0, phase 1 occupies the region of y ∈ (y 0 , y 2 ] and keeps stationary. phase 2 is in between y 1 and y 0 , and moves to the right with a unity speed.
phase 3 is at the rest of the domain and moves to the left with a unity speed. Here we have y 0 =0.5, y 1 =0.25, and y 2 =0.75. A sinusoidal vertical velocity is applied here to model the perturbation with an amplitude of 0.05 and wavelength of 2π. More details of the problem setup are available 23 .
For the NNs, the encoder consists of a single feed-forward layer, recurrent layers with 64 LSTM cells, and two hidden layers with 32 neurons each. The activation function is
ReLU. The learning rate is set to be 10 Where the order parameter of a phase is 1, only the corresponding phase occupies that location. Where the order parameter of a phase is -1, the corresponding phase is absent at that location. Here, we omit the comparison between the prediction results and the ground truth, because we cannot distinguish their differences according to the plots (the plots are almost the same).
We first illustrate the necessity of implementing MCBOM in Section III B, after obtaining predictions from the first NN model in Section III A. For a clear presentation, we only show the results of the first phase (phase 1) in Figure 4. For comparison, we present the results from the first NN model in Section III A, and the training loss converges, which, in our experiment, needs about 400 epochs. Figure 4(a) shows the mass change of phase 1 in the entire domain versus time. The ground truths of the results are obtained from the numerical solver in Huang et al. 24 . Therefore, the results we show are comparing to those with the method in Huang et al. 24 . It is clear that the mass of the phase is timeindependent, after applying MCBOM, while the mass is fluctuating from the output of the NN model in Section III A. In other words, mass gain or loss appears in the neural network prediction, which violates the principle of mass conservation. This issue is successfully addressed with the help of MCBOM. ful for long time prediction. Therefore, implementing MCBOM not only exactly enforces the physical constraints for the order parameters, but also produces more accurate and reliable long-time predictions. Figs. 6(a) and 6(b) show the momentum in the X and Y directions predicted by the proposed NN in Section III C or the traditional NN model without including the momentum conservation in the cost function. The momentum predicted from the NN is almost independent of time, and its value is very close to the ground truth. On the other hand, × 1] with the no-slip boundary condition. The number of grid points is [512×128]. The time step is ∆t = 5× 10 -2 . We used the same environment and equipment as in the HSL simulation, which took 5 days to train the neural network. The flows is at t = 0. More details of the problem setup are available 23 .
The PCNN here is identical to the one in the Section IV A except that it has 128 LSTM cells. Part of the predictive results is shown in Figure 7. not the case without the help of MCBOM. The results given by the NN model oscillate between 0.02 and 0.04. Figs. 9 shows the prediction of the momentum in the x direction.
Results predicted by the NN model will drop from -0.01 to -0.07. On the other hand, the prediction of the proposed NN including a physics-informed penalization of momentum conservation, is much closer to the validation data.
In the present study, a PCNN model is developed to predict MPFs. It consists of the first NN model to predict the order parameters, the MCBOM to correct the order parameters predicted by the NN model in order to exactly satisfy the mass conservation and the summation and boundedness constraints for the order parameters, and the second NN model to predict the flow velocity with a penalization of the momentum conservation using the density of the fluid mixture computed from the corrected order parameters.
We demonstrate that the prediction of the order parameters using only the NN model violates all the physical constraints. Hence, nonphysical behaviors, such as mass loss, local voids, or overfilling, are observed in our numerical experiments. After combining the NN and MCBOM, all those issues are thoroughly removed. Moreover, the accuracy and reliability of long-time predictions are improved when including the MCBOM. The momentum conservation is acceptably reproduced, although not exactly, by the proposed NN with physical constraints to predict the velocity. However, the traditional NN model fails to predict such a physical behavior.
In summary, we consider the NN model as the basic structure of the PCNN for MPFs, and enforce the physical constraints either implicitly by penalization or explicitly (or exactly) by MCBOM. The proposed PCNN model is an effective tool to predict the dynamics of the MPFs.