Physically plausible and conservative solutions to Navier–Stokes equations using physics-informed CNNs

1.   Introduction

Partial differential equations (PDEs) are ubiquitous in scientific fields, such as mechanics, engineering, and economics, and are often solved using numerical methods. Although numerical methods, including the finite difference method (FDM), finite volume method (FVM), and finite element method (FEM), have made successful progress in solving PDEs over the past few decades, their high computational overhead still makes the simulation difficult in many scenarios. For instance, conventional numerical solvers are computationally expensive for complex systems with multiscale/multiphysics features, e.g., optimization^[1], inverse problems^[2], and uncertainty quantification^[3].

In recent years, with the development of deep learning, solving PDEs using neural networks has attracted researchers’ interest. Neural networks show strong competitiveness compared with numerical methods in solving PDEs because of their generalization ability, lower implementation cost, and lower computational cost^[4]. Since the actual values of PDE solutions are expensive to acquire in real-world scenarios, it is not practical to train neural network models with a large amount of labeled data. In addition, in the absence of physical constraints, neural networks may produce predictions that violate the laws of physics and do not generalize well. Raissi et al.^[2] proposed a physics-informed neural network (PINN) for solving PDEs. PINN uses the automatic differentiation technique^[5] to calculate partial derivatives and uses the equation residuals and boundary/initial conditions as penalty terms to construct the neural network’s loss function. By combining physical information with neural networks, PINN effectively reduces the amount of labeled data required for training. However, the slow convergence of PINN limits its performance in solving PDEs. The physics-informed convolutional neural network (PICNN), a variant of PINN, achieves the better results on a series of PDEs due to better scalability and faster convergence of convolutional neural networks (CNNs). Besides, the parameter-sharing property of CNNs makes PICNN more effective in learning spatial dependencies. PICNN computes partial derivatives using central difference convolutional filters, whose parameters are derived by the finite difference method and encode PDEs into the loss function^{[6, 7]}. Compared to PINN, PICNN has the following three advantages. First, PICNN is more efficient in computing higher-order partial derivatives. PINN uses automatic differentiation to calculate partial derivatives, which makes the computational overhead increase exponentially as the order of the PDE increases. PICNN uses numerical methods to calculate partial derivatives, which makes the computational overhead increase linearly as the order of the PDE increases. Second, PICNN imposes boundary conditions instead of using them as penalty terms of the loss function, so that the boundary conditions can be enforced. Third, PICNN leverages CNNs for modeling nonlinear dynamics, which is capable of well portraying the local solution details^[8].

Although PICNN has achieved promising results in solving many PDEs, its performance is still unsatisfactory for equations containing convection terms such as convection-diffusion equations and Navier–Stokes equations. Navier–Stokes equations are widely used in the field of computational fluid dynamics and numerical heat transfer. The Navier–Stokes equations are characterized by strong directionality in the convection terms they contain, which describes the directional motion of the fluid. Under the action of convection, the disturbance at a certain point can only be transmitted downstream but not in reverse. Existing PICNN-based methods^{[9, 10]} discretize the equations using a central difference scheme, which in some cases leads to implausible solutions, i.e., oscillations of the solutions. Besides, the discrete scheme derived from the FDM does not have conservation properties, which can further lead to inaccurate predictions.

In this study, we try to properly combine numerical methods with CNNs to ensure that more accurate physical information is encoded into the neural networks so that the predicted solutions of the Navier–Stokes equations are more in line with physical laws. First, we use the FVM to derive the discrete schemes of the conserved PDEs to ensure that the resulting discrete equations are conserved. The use of discrete equations with conservative forms can describe physical systems more accurately than equations with nonconservative forms, resulting in solutions that are more in line with physical laws^[11]. Second, we use the second-order upwind difference scheme to discretize the Navier–Stokes equations for constructing the loss function. The second-order upwind difference scheme can obtain information from the upstream properly, which avoids the introduction of downstream information into the discrete equations. Third, we hard impose boundary conditions to allow the PICNN to train and predict without any labeled data.

2.   Related work

In many scientific computing scenarios, it is difficult to obtain labeled data (e.g., coordinates and corresponding physical quantities). Using purely data-driven deep learning methods to solve scientific computing problems is neither cost-effective nor able to obtain solutions that obey the laws of physics. By combining prior knowledge with neural networks, solutions can be learned with less data^{[2, 12–14]} or even without any data^{[6, 15, 16]}. Based on PINN^[2], a series of PINN variants have been proposed to address more challenging problems, such as UQPINN^[17], PPINN^[18], fPINN^[19], vPINN^[20], and B-PINNs^{[13, 21]}.

Recent studies have shown that CNNs are suitable for solving large-scale, high-dimensional spatiotemporal problems due to their parameter-sharing features^[22]. For steady-state PDEs, Zhu et al.^{[15, 23]} designed a physics-informed convolutional encoder-decoder structure to learn solutions without using any labeled data. Gao et al.^[9] designed a geometric adaptation strategy, which transformed the irregular geometric domain into a regular rectangular domain through a coordinate transformation so that the irregular domain can be processed by CNN. Geneva et al.^[6] proposed the AR-DenseED method, which used a PICNN-based deep autoregressive model to predict dynamical PDEs. Ren et al.^[10] designed a surrogate model based on the AR-DenseED method to learn the time evolution of dynamical partial differential equations.

Existing studies have achieved promising results in solving fluid and heat transfer problems using PICNN^[24–28]. Most of these methods use the FDM-based central difference scheme to discretize the convection term. Bar-Sinai et al.^[29] introduced the FVM to CNN, enabling the model to generate coefficients for approximating spatial derivatives. However, the method proposed in Ref. [29] is data-driven and does not work in scenarios where labeled data are expensive and scarce.

3.   Methodology

3.1   Navier–Stokes equations

Eq. (1) represents the steady-state Navier–Stokes equations in a two-dimensional Cartesian coordinate system, where $\rho$ is the fluid density, $u$ and $v$ are the flow velocities in the x-axis and y-axis directions, respectively, $\mu$ is the viscosity, and $p$ is the pressure. The terms $\dfrac{\partial{(\rho u u)}}{\partial{x}}$, $\dfrac{\partial{(\rho v u)}}{\partial{y}}$, $\dfrac{\partial{(\rho u v)}}{\partial{x}}$, and $\dfrac{\partial{(\rho v v)}}{\partial{y}}$ in the first two rows of Eq. (1) are convection terms, which have strong directionality and describe the directional motion of the fluid. The terms $\dfrac{\partial}{\partial{x}}\left( {\mu\dfrac{\partial{u}}{\partial{x}}} \right)+\dfrac{\partial}{\partial{y}}\left( {\mu \dfrac{\partial{u}}{\partial{y}}} \right)$ and $\dfrac{\partial}{\partial{x}}\left( {\mu\dfrac{\partial{v}}{\partial{x}}} \right)+\dfrac{\partial}{\partial{y}}\left( {\mu \dfrac{\partial{v}}{\partial{y}}} \right)$ in the first two rows of Eq. (1) are the diffusion terms, which describe the irregular thermal motion of molecules.

In Eq. (1), the first two rows represent the momentum conservation equations and the third row represents the mass conservation equation. To use CNN to solve the Navier–Stokes equations, we first need to divide the domain into grids and use the coordinates of the grid points as the input to the CNN. The output of the CNN is set to the values of the unknown function at the corresponding coordinates. Then, we need to embed the equations into a loss function, which is why we call it physics-informed CNN (PICNN). More details will be introduced in the following subsections.

3.2   Second-order upwind difference scheme

To embed the Navier–Stokes equations into the loss function of PICNN, we first need to discretize the equations. This subsection presents the procedure for discretizing the Navier–Stokes equations using the second-order upwind difference scheme derived by FVM. Eq. (1) is the conserved form of the Navier–Stokes equations, and FVM is more suitable than the FDM for discretizing equations in conserved form. A detailed analysis of the limitations of existing FDM-based methods and the central difference scheme are presented in Appendix A.1.1 and A.1.2.

First, we discretize the momentum equation in the horizontal direction (i.e., the first row in Eq. (1)). As shown in Fig. 1, the upper, lower, left, and right boundaries of the control volume of point $P$ are $n$, $s$, $w$, and $e$, and the upper, lower, left, and right adjacent points of point P are N, S, W, and E, respectively. By integrating the equation over the control volume of $P$, we can obtain Eq. (2). Further expanding Eq. (2), we can get Eq. (3).

Figure  1.  Domain meshing of two-dimensional equations.

DownLoad: Full-Size Img PowerPoint

We define the flow rate on interface $e$ as $F_e=(\rho u)_e \Delta y$, and similarly the flow rates on interfaces $w$, $n$, and $s$ are defined as $F_w=(\rho u)_w\Delta y$, $F_n=(\rho v)_n\Delta x$, and $F_s=(\rho v)_s\Delta x$, respectively. We define the diffusion conductance on interfaces $e$, $w$, $n$, and $s$ as $D_e=\dfrac{\mu \Delta y}{(\delta x)_e}$, $D_w=\dfrac{\mu \Delta y}{(\delta x)_w}$, $D_n=\dfrac{\mu \Delta x}{(\delta y)_n}$, and $D_s=\dfrac{\mu \Delta x}{(\delta y)_s}$, respectively^[30]. In this work, we divide the domain equally so that $(\delta x)_e=(\delta x)_w=(\delta y)_n=(\delta y)_s=\Delta x=\Delta y$. Thus we can get $D_e=D_w=D_n=D_s$, which we write as $D$ for brevity. We discretize the $\dfrac{\partial u}{\partial x}$ and $\dfrac{\partial u}{\partial y}$ at the interface using linear profile, i.e., discretize them using the central difference. We use the mean value of the pressure of the adjacent grid points to the left and right to represent the pressure values $p_e$ and $p_w$ on the interface. By processing Eq. (3) as described above, we can obtain Eq. (4).

The value of velocity $u_e$ at the interface $e$ in Eq. (4) is shown in Eq. (5). The values of $u_w$, $u_n$, and $u_s$ are similar to $u_e$, that is, we take two points in the upwind direction and combine them according to the coefficients in Eq. (5). This is where the second-order upwind difference scheme comes in.

By bringing Eq. (5) into Eq. (4), we can replace all the velocities at the boundaries in Eq. (4) with the velocities at the grid points, which are the predicted values of PICNN. Note that the velocity at the boundary in Eq. (5) varies with its sign. To simplify the equation, we use a compact form as shown in Eq. (6).

By substituting the compact form given by Eq. (6) into Eq. (4), we can obtain the second-order upwind scheme of the momentum equation in the horizontal direction, as shown in Eq. (7), where $P_\Delta=F/D$ and $\delta x$ represents the grid size. $P_\Delta$ is the Péclet number^[30], which represents the relative magnitude of convection and diffusion. Péclet number can help us analyze the performance of different discrete schemes.

The second-order upwind scheme of the momentum equation in the vertical direction (i.e., the second row in Eq. (1)) is similar to the equation in the horizontal direction and is detailed in Appendix A.2.2. We discretize the mass equation (i.e., the third row in Eq. (1)) using the central difference scheme because it does not contain convection terms. The derivation of the discrete scheme of the mass equation is presented in Appendix A.2.1. And we give the derivation of the central difference scheme in Appendix A.2.3 for comparison.

3.3   Construction of loss function

Eq. (7) gives the discrete scheme of the Navier–Stokes equations at grid point $P$, which is formed by combining the predicted values at point $P$ and its neighbors. When the predicted value of PICNN at point $P$ is close to its real value, the value calculated by the left side of Eq. (7) should be close to $0$. We set the left side of Eq. (7) as the function ${\cal{F}}_u$, and ${\cal{F}}_{u}(P)$ represents the value of the function ${\cal{F}}_u$ at point $P$, that is, the degree to which the predicted value of PICNN satisfies the first equation in the Navier–Stokes equations. Similarly, we set the left side of the discrete scheme of the vertical direction equation (see Eq. (A9)) to be the function ${\cal{F}}_{v}$, and the left side of the discrete scheme of the mass equation (see Eq. (A8)) to be the function ${\cal{F}}_\text{mass}$.

The loss function consists of three parts, which respectively represent the degree of satisfaction of the predicted value with the three equations in the Navier–Stokes equations. For each part, we compute the ${\cal{F}}$ values at all grid points and compute the MSE value, as shown in the first row of Eq. (8), where $P_i$ represents the grid point, $\varOmega$ represents the definition domain, and $n$ represents the number of grid points. Because the ${\cal{F}}(P_i)$ is compared with $0$, we abbreviate $({\cal{F}}(P_i)-0)^2$ as ${\cal{F}}(P_i)^2$. Finally, we set the total loss function ${\cal{L}}$ to be the sum of the three partial loss functions, as shown in the second row of Eq. (8).

3.4   Boundary encoding

The solutions of the steady-state PDEs are determined by both the equations and the boundary conditions. A Dirichlet boundary condition gives the value of the solution at the boundary, and a Neumann boundary condition gives the derivative or partial derivative of the solution at the boundary. In this work, we adopt the encoding method of boundary conditions proposed by Gao et al.^[9] As shown in Fig. 2 (left), Dirichlet boundary conditions can be imposed by applying constant padding, which is independent of the internal field and does not vary during each iteration. For Neumann boundary conditions, we can use the finite difference method to calculate the values at the boundary points from the predicted value of the interior points adjacent to the boundary points. Because the Neumann boundary conditions depend on the internal field, the padding values are different at each iteration. By imposing the boundary conditions hard, we can guarantee that they are enforced during training. Therefore, the PICNN model can be trained without any labeled data when the boundary conditions are determined.

Figure  2.  Hard-imposition of boundary conditions (BC). The left side represents the Dirichlet boundary conditions and the right side represents the Neumann boundary conditions.

DownLoad: Full-Size Img PowerPoint

The first row of Eq. (9) is the commonly used fluid boundary condition at the outlet, which is a Neumann boundary condition. By using the FDM, we can obtain the numerical relationship of the boundary point and its adjacent interior point, as shown in the second row of Eq. (9). Furthermore, we can obtain $v_\text{boundary}=v_\text{internal}$, which we can use for boundary padding.

3.5   Model training and prediction

Fig. 3 illustrates our PICNN architecture for solving Navier–Stokes equations. The input of the model is image-like data of 2 channels, which are composed of x-coordinates and y-coordinates separately. Considering the order of magnitude difference in the values of different variables predicted by the model (e.g., the streamwise velocity component can be orders of magnitude greater than the spanwise velocity component), we use sub-CNNs to predict them separately. Refs. [16, 31] demonstrated the superiority of using subnets for multivariate regression in data-driven and data-free scenarios. We use three sub-CNNs to predict velocity components and pressure ($u$, $v$, and $p$), each with the same four-layer convolutional network structure. Each layer of the sub-CNN has trainable filters with a kernel size of $5 \times 5$, and 2D convolutional operations with padding of 2 and stride of 1. It is important to emphasize that the use of sub-CNNs is not to avoid the use of normalization. If necessary, sub-CNNs and normalization can be used at the same time. In our cases, which are represented in Section 4, the input data range is between 0–1, so we did not use normalization.

Figure  3.  The PICNN architecture for Navier–Stokes equations.

DownLoad: Full-Size Img PowerPoint

For the first three convolutional layers, each layer is followed by a ReLU activation function and a batch normalization layer. The last convolutional layer is followed by an upsampling layer to expand the final output to the desired size. The variables predicted by the sub-CNNs are used to calculate the loss function after boundary encoding. After minimizing the PDE residuals, the prediction from the PICNN can approximate the true values.

The predicting process is similar to the training process, taking the coordinate values and passing them through the already trained sub-CNNs. The final predictions can be obtained after boundary encoding the outputs of the sub-CNNs.

4.   Experiments

Datasets and evaluation metric. We evaluate the accuracy and physical plausibility of our proposed methods on steady-state Navier–Stokes equations under different scenarios, including flow through a cavity, convective heat transfer, and a lid-driven cavity flow. We use the numerical solutions calculated by COMSOL, a commercial numerical calculation software, as the ground truth for comparison. The relative error metric is defined as $\text{Error}=||\tilde{u}-u||_{L2}/||\tilde{u}||_{L2}$, where $\tilde{u}$ is the numerical solution and $u$ is the predicted solution.

Baselines. The proposed method (we abbreviate it as SUS-FVM) is compared to the CD-FVM method, which is a PICNN model with a central difference scheme derived from the FVM to build the loss function. We also compare our proposed method with existing convolutional filter-based PICNN methods^[9], including methods based on the FDM-derived central difference scheme (CD-FDM) and second-order upwind difference scheme (SUS-FDM), to demonstrate that our method achieves higher prediction accuracy by using FVM. In addition, we also compare the PINN method, which is an multilayer perceptron (MLP) model based on automatic differentiation.

Implementation. We train all the models using the Adam optimizer with 200k epochs. The learning rate is set as 0.01 and is reduced by 50% every 50k epochs. We use the kaiming normalization to initialize the weights of convolutional layers. Other details for training can be found in the following subsections. To demonstrate the advantages of our proposed method over the baselines introduced above, we only optimize the loss function while keeping the network structure unchanged for a fair comparison. For more details, check the codes and datasets in the Data availability section.

4.1   Flow through a cavity

We choose the spatial domain size as $[0,1]\times[0,1]$ and divide it into $50\times50$ ($\delta x=0.02$) grids evenly. We set the inlet on the left side of the domain and the outlet on the upper side. The lengths of inlet and outlet are both set to $d=0.3$. The remaining boundaries are set as no-slip walls on which the velocity is 0. The fluid at the inlet is set to a constant velocity of $init\_u=1$ and the pressure at the outlet is 0. The velocity at the outlet meets the Neumann boundary condition $\dfrac{\partial v}{\partial y}=0$. We select Navier–Stokes equations with Reynolds numbers ($Re$) of 30150, and 300 for evaluation. The Reynolds number is used to describe the flow of fluids, defined by $Re=\rho(init\_u)d/\mu$, where $\rho=1$, $init\_u=1$, $d=3$, and $\mu=0.01,\; 0.002, \;0.001$ for $Re=30,\;150,\;300$, respectively. We change the Reynolds number by adjusting the viscosity $\mu$ of the fluid. As the Reynolds number increases, the convection becomes stronger.

Table 1 shows the relative error of predicted pressure and velocity ($\sqrt{u^2+v^2}$) for the Navier–Stokes equations. We can find that as the Reynolds number increases from $Re=150$ to $Re=300$, the relative error of prediction for the CD-FVM increases dramatically, while the relative error of prediction for the SUS-FVM increases slightly. When the $Re$ is low ($Re=30$), CD-FVM has a higher prediction accuracy. And when the $Re$ is high ($Re=300$), the SUS-FVM scheme we used shows higher accuracy in the prediction of velocity and pressure.

Table  1.  Relative errors of predictions for the flow through a cavity.

		CD-FVM	SUS-FVM
$Re=30$	Pressure	0.0478	0.0675
	Velocity	0.0318	0.0393
$Re=150$	Pressure	0.0456	0.0855
	Velocity	0.0404	0.0389
$Re=300$	Pressure	0.0967	0.0839
	Velocity	0.1280	0.0486

 | Show Table

DownLoad: CSV

Fig. 4 shows the comparison of solutions between CD-FVM and SUS-FVM for the flow through a cavity. The fluid velocity is close to 0 in most areas of the domain, and only close to 1 on the way from the inlet to the outlet. SUS-FVM predicts the velocities well and has a higher prediction accuracy and better physical plausibility in the regions we care about (inlet, outlet, and pathways between them). To better illustrate the good prediction performance of SUS-FVM in local areas, we select a section parallel to the outlet ($y=0.96$) and a section perpendicular to the outlet ($x=0.50$) to observe the velocity distribution on them, as shown in Fig. 5. In Fig. 5, we can find that the predicted solutions obtained using SUS-FVM can perfectly capture the changes at locations with large velocity gradients. Compared to SUS-FVM, the solutions obtained using CD-FVM exhibit strong oscillations near the outlet, which extend vertically up to the position of $y=0.5$ (Fig. 5b). The predicted solutions from $x=0.45$ to $x=0.65$ in Fig. 5a more clearly show the difference between the solutions obtained using CD-FVM and the real values. In practice, we cannot tolerate the oscillating solutions, and the SUS-FVM scheme we use effectively reduces the oscillation of the predicted solutions, making PICNN’s prediction of the Navier–Stokes equations reliable.

Figure  4.  Comparison of solutions between CD-FVM and SUS-FVM for flow through a cavity (Re = 300).

DownLoad: Full-Size Img PowerPoint

Figure  5.  Comparison of velocity accuracy on section y=0.96 and section $x=0.50$ ($Re=300$).

DownLoad: Full-Size Img PowerPoint

Advantage of conservative schemes. Both CD-FDM and SUS-FDM use discrete schemes derived from nonconservative equations. The method we proposed and the CD-FVM use the discrete schemes derived from conservative equations. All the schemes are tested on dataset $Re=30$ and $Re=300$. From Table 2 we can find that the scheme derived from the conservative equations gives better predictions on both pressure and velocity than those derived from or using the nonconservative equations. When $Re=30$, the relative errors for CD-FDM and SUS-FDM are an order of magnitude larger than those for CD-FVM and SUS-FVM, respectively. When $Re=300$, the CD-FDM and SUS-FDM methods even diverse. The experimental results illustrate the advantages of deriving discrete schemes with conservation properties in prediction accuracy.

Table  2.  Relative errors of predictions for conservative and nonconservative schemes of the Navier–Stokes equations.

		CD-FVM	CD-FDM	SUS-FVM	SUS-FDM
$Re=30$	Pressure	0.0478	0.1392	0.0675	0.7760
	Velocity	0.0318	0.1078	0.0393	0.3281
$Re=300$	Pressure	0.0967	–	0.0839	–
	Velocity	0.1280	–	0.0486	–
Dash “–” means divergence and error cannot be calculated.

 | Show Table

DownLoad: CSV

Comparison with PINN methods. We contrast the proposed method with DeepXDE^[32] as shown in the Table 3. DeepXDE is a library for scientific machine learning and physics-informed learning, which includes the PINN method and integrates various strategies to improve the accuracy of PINN, including residual-based adaptive sampling^[33], gradient-enhanced PINN^[34], and PINN with multi-scale Fourier features^[35].

Table  3.  Relative errors of predictions for DeepXDE and SUS-FVM for the flow through a cavity.

		DeepXDE-s	DeepXDE-l	SUS-FVM
$Re=30$	Pressure	0.4420	0.3588	0.0675
	Velocity	0.2233	0.1715	0.0393
$Re=150$	Pressure	0.9956	0.5955	0.0855
	Velocity	0.6086	0.4062	0.0389
$Re=300$	Pressure	1.5442	0.6548	0.0839
	Velocity	0.7859	0.4522	0.0486

 | Show Table

DownLoad: CSV

We train the DeepXDE model with two data volumes separately. One uses a small amount of data (DeepXDE-s), the amount of which is the same as the amount of training data used by the SUS-FVM model, i.e., $49\times49$ interior points and $200$ boundary points. The other uses a large amount of data (DeepXDE-l), that is, $50000$ interior points and $5000$ boundary points. From Table 3, we can find that the prediction accuracy of DeepXDE-l is higher than that of DeepXDE-s, thanks to its amount of training data. However, the prediction error of SUS-FVM is an order of magnitude smaller than that of DeepXDE-l and DeepXDE-s. The experimental result shows that although a larger amount of data increases the generalization of PINN models, limited by the fully-connected neural networks and the soft imposition of boundary conditions, the prediction accuracy of PINN is not as good as that of PICNN.

4.2   Convective heat transfer

Compared with the Navier–Stokes equations, the convective heat transfer problems add a temperature equation, as shown in Eq. (10), where $T$ represents the temperature, $C_p$ represents the heat capacity at constant pressure, $k$ represents the thermal conductivity and $\rho$ represents the density of the medium. We set the inlet of the fluid to the left side of the domain and the outlet to the right side. The velocities $u$ and $v$ are calculated by commercial software on the Navier–Stokes equations with $Re=500$. By varying the value of thermal conductivity $k$, we obtain different datasets. For temperature, we choose the Dirichlet boundary conditions and set the temperature of the left wall to $293.15$ and the temperature of the remaining walls to $333.15$.

Table 4 shows the relative errors of predicted temperature for the convective heat transfer problems. Since the maximum velocity, in this case, is around 1, we can infer that when $k=10$ (data1), the grid Peclet number $P_\Delta=2$. As $k$ gradually decreases to 0.1 (data2 to data4), $P_\Delta$ gradually increases to 200. As $P_\Delta$ increases, the relative error of the prediction of the CD-FVM increases dramatically, while the relative error of the SUS-FVM increases slowly, which is an order of magnitude smaller than that of the CD-FVM. Using data3 as an example, as shown in Fig. 6, we can see that using CD-FVM for prediction yields solutions with strong oscillations near the outlet. We select a section perpendicular to the outlet ($y=0.64$) in Fig. 6 to further observe the temperature distribution on it, as shown in Fig. 7. The strongly oscillating solutions obtained using CD-FVM completely lose the physical plausibility, while the solutions obtained using SUS-FVM fit the real solutions precisely.

Table  4.  Relative errors of predictions ($T$) for the convective heat transfer problems.

	CD-FVM	SUS-FVM
Data1: $k=10$	0.0023	0.0027
Data2: $k=1$	0.0196	0.0042
Data3: $k=0.5$	0.0338	0.0076
Data4: $k=0.1$	0.0922	0.0293
$Re=500$, $\rho=1$, $C_p=1000$

 | Show Table

DownLoad: CSV

Figure  6.  Comparison of solution accuracy between SUS-FVM and CD-FVM schemes for the convective heat transfer problems (data3).

DownLoad: Full-Size Img PowerPoint

Figure  7.  Comparison of solution accuracy on section $y=0.64$ between SUS-FVM and CD-FVM schemes (data3).

DownLoad: Full-Size Img PowerPoint

Figure  8.  Comparison of predicted and true velocities for the lid-driven cavity flow.

DownLoad: Full-Size Img PowerPoint

4.3   Lid-driven cavity flow

Lid-driven cavity flow is a common benchmarking case in computational fluid dynamics. We use it to show the performance of our proposed method at high Reynolds numbers. We set the horizontal movement speed of the lip (upper boundary of the domain) to $1$, and set the other walls to be no-slip walls. By adjusting the viscosity $\mu$, we obtained three datasets with Reynolds numbers of $100$, $500$, and $1000$, respectively.

As shown in Table 5, SUS-FVM outperforms CD-FVM in prediction on all three datasets. When the Reynolds number is low ($Re=100$), the vortex in the cavity is not obvious, which is shown in Fig. 8, and both CD-FVM and SUS-FVM perform well in velocity prediction. When the Reynolds number is high ($Re=500, 1000$), the vortex in the cavity is obvious, which makes the prediction performance of CD-FVM and SUS-FVM for velocity decrease. But it can be found that the prediction accuracy of SUS-FVM decreases more slowly than that of CD-FVM. The results show that our proposed method is still effective in predicting the Navier–Stokes equations in the case of high Reynolds numbers and outperforms the existing central difference based methods.

Table  5.  Relative errors of predictions for the lid-driven cavity flow.

		CD-FVM	SUS-FVM
$Re=100$	pressure	0.4277	0.3986
	velocity	0.1082	0.1074
$Re=500$	pressure	0.2873	0.2125
	velocity	0.1644	0.1317
$Re=1000$	pressure	0.3737	0.3395
	velocity	0.2936	0.2768

 | Show Table

DownLoad: CSV

4.4   Discussion

In this work, our main study is solving the incompressible Navier–Stokes equations under laminar flow, i.e., Navier–Stokes equations under low Reynolds number, using physics-informed CNN. For turbulent flow scenarios, i.e., high Reynolds number, the prediction difficulty of Navier–Stokes equations is much higher than that of laminar flow scenarios. We need to optimize the neural network structure so that the model can predict the local details of turbulent flow well, which is our future research direction.

Furthermore, for compressible Navier–Stokes equations such as shock wave, our method is currently not suitable to solve such problems. There is a huge difference in fluid flow between shock wave and laminar flow. To better simulate the shock wave, it is necessary to accurately predict the fluid properties near the obstacle. However, the PICNN method we use is currently unable to describe the scene with obstacles in the domain of definition. It is a challenging problem and deserves further investigation in our future work.

5.   Conclusions

In this study, we focus on solving Navier–Stokes equations using PICNN. Aiming at the problem that PICNN will generate oscillating predictions when solving the Navier–Stokes equations, we propose to use the FVM-based second-order upwind difference scheme to construct the loss function. We give formulation derivations and conduct experiments on the Navier–Stokes equations under different scenarios to demonstrate that our proposed method effectively improves the physical plausibility and accuracy of the predictions.

Our proposed method is not limited to convolutional neural networks but applies to other network architectures using physically constrained loss functions. Because of the introduction of grid-scale, our method can solve PDEs under variable-step rectangular grid division. Furthermore, our method can be extended to graph neural networks to solve irregular spatial domains and grids. Besides, the proposed method can be used for the loss function calculation of PDEs solving models for different purposes, such as data-driven PDEs discovery models and boundary surrogate models, to improve the solution accuracy and physical plausibility. Finally, this work focuses on solving steady-state PDEs. For solving dynamic PDEs, we can further derive their discrete schemes and leverage suitable network architectures (e.g., LSTM) in future work.