A Direct-Forcing Immersed Boundary Method for Incompressible Flows Based on Physics-Informed Neural Network

Abstract

The application of physics-informed neural networks (PINNs) to computational fluid dynamics simulations has recently attracted tremendous attention. In the simulations of PINNs, the collocation points are required to conform to the fluid–solid interface on which no-slip boundary condition is enforced. Here, a novel PINN that incorporates the direct-forcing immersed boundary (IB) method is developed. In the proposed IB-PINN, the boundary conforming requirement in arranging the collocation points is eliminated. Instead, velocity penalties at some marker points are added to the loss function to enforce no-slip condition at the fluid–solid interface. In addition, force penalties at some collocation points are also added to the loss function to ensure compact distribution of the volume force. The effectiveness of IB-PINN in solving incompressible Navier–Stokes equations is demonstrated through the simulation of laminar flow past a circular cylinder that is placed in a channel. The solution obtained using the IB-PINN is compared with two reference solutions obtained using a conventional mesh-based IB method and an ordinary body-fitted grid method. The comparison indicates that the three solutions are in excellent agreement with each other. The influences of some parameters, such as weights for different loss components, numbers of collocation and marker points, hyperparameters in the neural network, etc., on the performance of IB-PINN are also studied. In addition, a transfer learning experiment is conducted on solving Navier–Stokes equations with different Reynolds numbers.

1. Introduction

Over the past few years, machine learning (ML) has permeated into various research areas of fluid mechanics [1], e.g., reduced-order modeling [2,3], wake-type clustering and classification [4,5,6,7], development of turbulence closure model [8,9,10], flow optimization and active control [11,12,13,14,15,16], to name a few. The successful application of ML to these areas relies on the availability of large-scale data, which are obtained from CFD simulations or experimental observations.

More recently, the physics-informed neural network (PINN) was proposed to solve partial differential equations (PDEs) [17,18,19]. PINNs use automatic differentiation to compute the derivatives, and the residuals of equations are incorporated into the loss function. Unlike the conventional practice in ML, PINN can be trained without any simulation (or observation) data or with only a small amount of data. Nowadays, PINN has been successfully applied to the simulations of incompressible [20,21,22,23] and compressible [24] flows.

PINN is not meant to be a replacement of traditional CFD solvers but instead a complementary approach. For solving a standard forward problem (one PDE system with fixed parameters), the state-of-art PINN cannot compete with the traditional CFD methods in terms of accuracy, efficiency and robustness. However, a surrogate model based on PINN is a good choice in some non-ideal situations where traditional CFD methods cannot handle or are very expensive to use, such as data assimilation with numerical solution, solving ill-posed problems due to a lack of boundary conditions and solving inverse problems. Other advantages of PINN include: (a) closed-form approximated solutions with continuous derivatives, (b) implementation simplicity on high-level marching learning platforms.

PINN can be regarded as a mesh-free method for solving PDEs since an explicit mesh is not needed. In PINN, the no-slip condition at the fluid–solid interface is imposed by incorporating velocity penalties at the boundary points into the loss function. As that in traditional mesh-based numerical methods, the boundary points are required to conform to the surfaces of immersed objects.

The immersed boundary (IB) method is an alternative technique for handling complex and moving boundaries in mesh-based methods. It utilizes non-boundary-conforming meshes in numerical discretization, and the no-slip condition on the surface of the immersed object is enforced by adding a volume force to the momentum equation. This method has recently gained its popularity due to the simplicity of implementation.

The combination of PINN and the IB method can eliminate the boundary-conforming requirement in point generation. This eases the point generation process in problems with complex boundaries and avoids redistributing the points in moving-boundary problems. These benefits motivate the development of a novel PINN in conjunction with the IB method in this work. Such benefits were also the main motivation for incorporating IB method into a smoothed particle hydrodynamics (SPH) flow solver [25]. The proposed IB-PINN can be considered a complementary approach to the classical (mesh-based) IB schemes (rather than a replacement). Here, we emphasize that the mesh-free nature of PINN can barely be considered a potential advantage in such a combination. This is because mesh generation in the mesh-based IB method is significantly simplified and not a challenging task anymore.

To the best of the authors’ knowledge, PINNs, in conjunction with IB method, are still barely explored. Balam et al. recently presented a PINN for solving two-dimensional elliptic equations with singular forces on an embedded interface [26]. The IB method involved in their work fell into the ’continuous-forcing’ sub-type. The forcing terms in the equations were explicitly computed via a jump condition at the interface, together with a smoothed delta function for spreading the singular force. Thus, no major modification to the original framework of PINN was needed in their work.

In this paper, we consider the ‘direct-forcing’ sub-type of IB methods. Unlike the one in [26], in this variant of IB methods, the forcing term is determined indirectly via a no-slip constraint imposed on the embedded interface. Due to the lack of an explicit expression for the forcing term, the loss components that account for the equation residual and non-slip condition on the embedded interface need to be reformulated.

The paper is organized as follows. In Section 2, we present an introduction to the direction-forcing IB method and the architecture of neural network used in this study. The results of the simulations using the proposed IB-PINN are presented in Section 3. Finally, the summary of the present work and suggestions on future directions are provided in Section 4.

2. Methods

2.1. Direct-Forcing Immersed Boundary Method

We consider the incompressible Navier–Stokes equations with a volume force added to the momentum equation. The governing equations can be written in a dimensionless form as:

$$\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}·\nabla )\mathbf{u}=-\nabla p+\frac{1}{Re}{\nabla }^2\mathbf{u}+\mathbf{f},$$

(1)

$$\nabla ·\mathbf{u}=0,$$

(2)

$$\mathbf{f}(\mathbf{x},t)={\int }_s\mathbf{F}\left(\xi \right(s,t\left)\right)\delta (\xi -\mathbf{x})ds,$$

(3)

$$\mathbf{u}\left(\xi (s,t)\right)={\int }_{\mathbf{x}}\mathbf{u}\left(\mathbf{x}\right)\delta (\mathbf{x}-\xi )d\mathbf{x}=0,$$

(4)

where $Re$ is the Reynolds number, $\mathbf{u}$, p and $\mathbf{f}$ denote the velocity, pressure, volume force (density), respectively. The spatial variable $\mathbf{x}$ denotes the position in the computational domain (which includes the region occupied by the immersed object). $\xi$ represents an embedded fluid–solid interface. s is the coordinate along the interface for parameterizing it. $\mathbf{F}$ represents the singular force that is added to the fluid–solid interface. Equation (3) is the formula for the transfer between the surface force and the volume force. Equation (4) is the no-slip constraint that must be satisfied at the surface of a stationary object.

For incompressible Navier–Stokes equations, the pressure p acts as a set of Lagrange multipliers, which enforce the divergence-free condition (Equation (2)). In direct-forcing IB methods, the surface (and volume) forces can also be regarded as a set of Lagrange multipliers that enforce the no-slip condition (Equation (4)) on the embedded interface. In numerical implementations, the forcing term and pressure can be obtained sequentially [27,28] or simultaneously in a coupled way [29,30] via splitting the system (Equations (1)–(4)).

In conventional mesh-based discretization, $\mathbf{u}$, p and $\mathbf{f}$ are defined at the grid points, while $\mathbf{F}$ is defined at a set of marker points that represent the interface. The integral operations in Equations (3) and (4) are replaced by summation, and the delta function is replaced by a smoothed delta function with compact support [31].

2.2. Physics-Informed Neural Network

Instead of using conventional mesh-based methods, in this work, a deep neural network (DNN) is used to approximate the solution of incompressible Navier–Stokes equations. The architecture of the DNN used here is similar to those proposed in some previous references [17,20,21,22] (see Figure 1). The DNN is composed of multiple hidden layers, and it can be expressed as:

$$\begin{array}{c}{\mathbf{z}}_0=(\mathbf{x},t),\hfill \\ {\mathbf{z}}_{\ell }={\sigma }_{\ell }(W_{\ell }{\mathbf{z}}_{\ell -1}+{\mathbf{b}}_{\ell }),1\le \ell \le L-1\hfill \\ {\mathbf{z}}_{\ell }=W_{\ell }{\mathbf{z}}_{\ell -1}+{\mathbf{b}}_{\ell },\ell =L,\hfill \end{array}$$

(5)

where ${\mathbf{z}}_{\ell }$ denotes the hidden variables of the ℓth layer, and $W_{\ell }$ and ${\mathbf{b}}_{\ell }$ denote the weight matrix and bias vector of the ℓth layer, respectively. ${\sigma }_{\ell }$ is a nonlinear activation function. Throughout this work, the hyperbolic tangent function is used as the activation function. The Xavier initialization technique is used to initialize the network [32]. The DNN takes the spatial-temporal variables $\{\mathbf{x},t\}$ as the inputs. The outputs of the last layer are used to approximate the solution $\{\mathbf{u},p,\mathbf{f}\}$. Please note that the surface force $\mathbf{F}$, which is an auxiliary variable for obtaining the volume force in some mesh-based IB methods, is not used here. The DNN is implemented using the TensorFlow machine learning platform, and the source code can be found in https://github.com/huangyi89/IB-PINN (accessed 12 January 2022).

PINN solves a PDE system by converting it into an optimization problem in which the loss function is minimized through iteratively updating the weights and biases. Generally speaking, the total loss for a PDE system is composed of equation residual loss, initial condition loss and boundary condition loss. In solving a PDE system in the framework of direct-forcing IB methods, special care needs to be taken in constructing the loss component associated with the non-slip boundary condition at the embedded interface. In this study, we propose incorporating velocity penalties at some marker points that represent the embedded interface into the total loss. Obtaining the forcing term by adding a velocity penalty to the loss function is analogous to obtaining the pressure by adding the residual of the continuity equation into the loss function in the PINNs proposed by Sun et al. [20] and Jin et al. [21]. This analogy is rooted in the fact that both the forcing term and pressure are Lagrangian multipliers for imposing constraint on the system [29,30]. In addition to the velocity penalty term, the total loss in the proposed IB-PINN also includes a force penalty term that is used to ensure compact distribution of volume force. To this end, the force penalties are applied at the collocation points, which are exterior to and beyond certain distance from the immersed interface. The compact distribution of the forcing term is a necessary condition for accurately mimicking the dynamic effect of an immersed object on fluid flow. In conventional mesh-based IB methods, such a condition is automatically satisfied by using an interpolation kernel function (i.e., smoothed delta function) with compact support [31].

As a result, the loss function used for training IB-PINN can be expressed as:

$$J=\lambda J_{PDE}+\alpha J_b+\beta J_{IB}+\gamma J_f,$$

(6)

$$J_{PDE}=\frac{1}{N_r}\left(\sum _{i=1}^{N_r}{∥{\mathbf{r}}_m\left({\mathbf{x}}_r^i,t\right)∥}^2+\sum _{i=1}^{N_r}{∥r_c\left({\mathbf{x}}_r^i,t\right)∥}^2\right),$$

(7)

$$J_b=\frac{1}{N_{b_1}}\sum _{i=1}^{N_{b_1}}{∥\mathbf{u}\left({\mathbf{x}}_b^i,t\right)-\widehat{\mathbf{u}}\left({\mathbf{x}}_b^i,t\right)∥}^2+\frac{1}{N_{b_2}}\sum _{i=1}^{N_{b_2}}{∥p\left({\mathbf{x}}_b^i,t\right)-\widehat{p}\left({\mathbf{x}}_b^i,t\right)∥}^2,$$

(8)

$$J_{IB}=\frac{1}{N_{IB}}\sum _{i=1}^{N_{IB}}{∥\mathbf{u}\left({\mathbf{x}}_{IB}^i,t\right)∥}^2,$$

(9)

$$J_f=\frac{1}{N_{r_1}}\sum _{i=1}^{N_{r_1}}{∥\mathbf{f}\left({\mathbf{x}}_r^i,t\right)∥}^2.$$

(10)

Here J denotes the total loss function, which is the summation of four components corresponding to the residual loss ($J_{PDE}$), the outer boundary loss ($J_b$), the boundary loss associated with the immersed boundary ($J_{IB}$) and the force compactness loss ($J_f$). As a proof-of-concept, only steady problems are considered in this work. Thus, the initial condition loss is neglected in Equation (6). $\lambda$, $\alpha$, $\beta$ and $\gamma$ are the weights for the aforementioned four loss components, respectively.

The equation residual loss can be further decomposed into the loss for the momentum Equation (${\mathbf{r}}_m$) and the loss for the continuity Equation ($r_c$). In computing the equation residual loss, the partial differential operations are performed using automatic differentiation (AD). $\widehat{\mathbf{u}}$ and $\widehat{p}$ denote the velocity and pressure that are specified on the outer boundaries, respectively. $N_r$ and $N_{IB}$ denote the numbers of collocation points and marker points, respectively. $N_{b_1}$ and $N_{b_2}$ denote the numbers of boundary points where velocity and pressure are specified, respectively. $N_{r_1}$ denotes the number of collocation points on which force penalties are applied.

The locations of the three sets of points, namely collocation points, boundary points and marker points, are illustrated in Figure 2. The collocations points are located inside the computational domain (including the area occupied by the immersed object) and on the outer boundary. The boundary points are located on the outer boundary where data for velocity or pressure are available. In this work, the boundary points and the collocation points coincide with each other on the outer boundary. Such a choice of boundary points is just for the simplicity of implementation. A different set of boundary points that does not belong to the collocation points can also be chosen [26]. The marker points are located inside the domain for representing the fluid–solid interface (i.e., immersed boundary). Please note that the marker points do not belong to the set of collocation points and are only used for computing the loss. The force penalties are applied at the collocations points that are located in the gray area (the area outside the cylinder and beyond a distance of Δ from the immersed boundary). Please note that no force penalties are applied at the collocation points that lie inside the immersed object.

3. Results

3.1. Problem Description and Numerical Settings

We use the proposed IB-PINN to simulate fluid flow past a circular cylinder that is placed in a channel [22]. Since only a steady case is considered here, the Reynolds number based on the diameter of the cylinder is set to 40, which is below the critical value for the onset of vortex shedding. The computational domain and boundary conditions are illustrated in Figure 3. At the inlet, a parabolic profile of the streamwise velocity is specified, and the crosswise velocity component is set to zero. Mathematically, the parabolic velocity profile is defined as

$$u(0,y)=\frac{4(H-y)y}{H^2},$$

(11)

where H denotes the height of the channel. At the outlet, the pressure is set to zero. The no-slip condition is imposed on the top and bottom walls of the channel.

A total of 72,365 collocation points (with 2564 outer boundary points) and 251 marker points are used in the simulation. The collocation points are placed at the vertices of a uniformly distributed Cartesian mesh with a grid width of h = 0.025. Please note that this mesh is generated for arranging the collocation points only and is not used explicitly in the IB-PINN simulations. The marker points are uniformly distributed on the surface of the cylinder. As a usual practice in IB methods, the spacing between two neighboring marker points is comparable with that between two neighboring collocation points. The distance Δ in Figure 2, which represents the supporting width of volume force, is set to 2 h.

The DNN used in the simulations is composed of six hidden layers, of which the number of layer neurons is 100. The weights for the loss components are: λ = 1.0; α = 2.0; β = 2.0; γ = 1.0. In the training process, we follow the strategy adopted in some previous works, such as [21,22], where certain numbers of iterations are first performed using the Adam optimizer; after that, the L-BFGS-B optimizer is used to finetune the DNN to obtain a more accurate solution. In this work, 10⁴ iterations in the Adam optimizer are first performed with a fixed learning rate of 5 × 10⁻⁴ before the L-BFGS-B optimizer is used. The training with the L-BFGS-B optimizer is terminated based on the increment tolerance (or if the total number of iterations reaches 10⁵). The training is performed on Tesla-V100 16G/32G GPU.

A reference solution is obtained by using a conventional mesh-based IB method. This numerical method is based on the discrete stream-function formulation for solving incompressible Navier–Stokes equations [28]. The spatial discretization and temporal advancement in the flow solver are of second-order accuracy. For detailed descriptions and validations of the mesh-based IB method, please refer to [28,33,34]. The uniform rectangular mesh aforementioned (with h = 0.025) is used in the simulation. The smoothed delta function used in the mesh-based IB method has a supporting width of 4 h (the half supporting width is equal to Δ) [35].

An additional reference solution is also obtained by performing a simulation using Ansys Fluent. This software uses an ordinary CFD method (with body-fitted mesh and velocity-pressure formation). The resolution of the simulation is comparable with that of the mesh-based IB method. The numerical settings of this simulation are described in Appendix A.

3.2. Predicted Solution vs. Reference Solutions

The convergence history of the training process is shown in Figure 4. It is seen that after the training is terminated, the total loss is reduced to a value less than 4 × 10⁻⁵, while the maximum value among the four loss components is around 2 × 10⁻⁵. The low values of the final loss indicate the convergence of training.

The contours of the streamwise and crosswise velocity components obtained using Fluent, mesh-based IB method and IB-PINN are displayed in Figure 5. The contours for the difference between the two IB solutions are also shown in the figure. The comparisons indicate that the three solutions are in excellent agreement qualitatively. Actually, an evident difference does exist in the interior of the cylinder between the solutions of the two IB methods. However, the interior flows in IB methods are merely fictitious and irrelevant to the prediction accuracy. Thus, in plotting the difference between the two IB solutions, data for the interior flows are removed, and the interior regions are filled with white color.

The streamlines for the three solutions are shown in Figure 6. Again, it is seen that the flow fields exterior to the cylinder are in excellent agreement. Evident difference in the fictitious flows obtained by the two IB methods can be visualized in the interior of the cylinder. This difference is mainly due to the fact the force compactness constraint is imposed inside the cylinder in the mesh-based IB method but not in IB-PINN.

The profiles of streamwise velocity along the horizontal centerline and vertical line passing through the center of the cylinder are shown in Figure 7. The quantitative agreement between the three solutions can be seen clearly. Again, larger discrepancies are only found on the segment inside the cylinder (corroding to the region where fictitious flows are produced in the two IB methods).

The pressure contours obtained using IB-PINN and Fluent are compared in Figure 8. The distributions of the pressure coefficient along the cylinder surface are also compared in Figure 9. A good agreement between the two solutions is clearly demonstrated. The slight under-prediction of the pressure coefficient near the front stagnation point in IB-PINN is associated with the smearing effect of diffusive-interface IB methods. The pressure prediction in IB-PINN can be improved by adopting a local refinement strategy in the generation of collocation points. Here, the pressure field for the mesh-base IB method is not available. This is because the mesh-based IB solver is based on the discrete stream-function formulation in which the pressure is eliminated.

The length and position of the separation bubble behind the cylinder, the separation angle and the drag coefficient are listed in

Table 1. The length and position of the separation bubble, the separation angle and the drag coefficient obtained using Ansys Fluent, mesh-based IB method and IB-PINN.

	$\mathit{l}/\mathit{D}$	$\mathit{a}/\mathit{D}$	$\mathit{b}/\mathit{D}$	${\mathit{\theta }}_{\mathit{s}}$	${\mathit{C}}_{\mathit{d}}$
Fluent	1.19	0.45	0.48	47.7°	2.17
IB	1.19	0.46	0.48	46.7°	2.21
IB-PINN	1.17	0.45	0.48	46.7°	2.18

for comparison. For the IB methods, the drag coefficient is computed by

$$C_d=\left(-\sum _{j=1}^{N_r}{f}_x{h}^2\right)/\left(\frac{1}{2}\rho U_{max}^{2}D\right),$$

(12)

where $f_x$ denotes the x-component of volume force $\mathbf{f}$; $\rho$, $U_{max}$ and D are (non-dimensional) density, maximum velocity and diameter, respectively, (i.e., $\rho$$U_{max}^{2}D=1$). For the solution of Fluent, the drag coefficient is computed by the integration of pressure and friction forces on the cylinder surface. From this table, it can be seen that the results from the three simulations also agree very well.

The distributions of the volume force obtained in the two simulations of IB methods are compared in Figure 10. It is interesting to note that although the velocity fields obtained in the two simulations are very similar, the two volume forces differ significantly. The dissimilarity in the distributions of volume force is partially due to the lack of the force compactness constraint inside the cylinder in the IB-PINN simulation. Another reason for the discrepancy between the two volume forces is that no interpolation kernel function is used in IB-PINN, while in the mesh-based IB method, an interpolation kernel function is used to regulate the distribution of volume force.

3.3. Influences of Some Parameters on the Performance of IB-PINN

In this section, we study how the performance of IB-PINN is affected by the variations of some parameters. The IB-PINN with the numerical settings aforementioned is regarded as the baseline model to compare with. The value of final loss after training and the relative $L_2$ error of velocity are used to quantify the convergence behavior and prediction accuracy, respectively. Here, the relative $L_2$ error is defined as

$${\epsilon }_{\phi }={{∥\widehat{\phi }-\phi ∥}_2/\phantom{{∥\widehat{\phi }-\phi ∥}_2∥\widehat{\phi }∥}∥\widehat{\phi }∥}_2,$$

(13)

where $\phi$ represents the two velocity components $(u,v)$ obtained by IB-PINN. $\widehat{\phi }$ represents an accurate solution, which is obtained using the conventional IB method, on a very fine mesh of h = 0.005. Please note that the fictitious velocity inside the circular cylinder is excluded in the computation of relative $L_2$ error.

The effects of the DNN architecture on the convergence of training and velocity error are tested by using DNNs with nine different configurations. To be more specific, three numbers of hidden layers (5, 6, 7) and three numbers of neurons per layer (90, 100, 110) are considered. The results obtained by the nine DNNs are summarized in

Table 2. Final loss and velocity error obtained in nine DNNs with different configurations.

		90	100	110
5	Final Loss	7.3764 × 10${}^{-5}$	8.1033 × 10${}^{-5}$	5.4935 × 10${}^{-5}$
	${\epsilon }_u$(%)	0.28	0.27	0.26
	${\epsilon }_v$(%)	0.64	0.84	2.65
6	Final Loss	5.8943 × 10${}^{-5}$	3.8017 × 10${}^{-5}$	4.4175 × 10${}^{-5}$
	${\epsilon }_u$(%)	0.28	0.16	0.25
	${\epsilon }_v$(%)	0.87	0.62	0.82
7	Final Loss	3.4572 × 10${}^{-5}$	4.1086 × 10${}^{-5}$	3.2451 × 10${}^{-5}$
	${\epsilon }_u$(%)	0.20	0.22	0.20
	${\epsilon }_v$(%)	1.18	0.96	0.60

. The low values of final loss indicate good convergence in training all DNNs. As to the prediction accuracy, the smallest velocity error is achieved in DNN with 6 hidden layers and 100 neurons per layer. This implies that for a fixed number of collocation points, a DNN with a more complex architecture (i.e., larger size) does not necessarily provide a more accurate prediction. Similar findings were also reported in some previous studies on PINNs [20,21,22,26].

Figure 11 shows the influence of the number of Adam iterations (before L-BFGS-B) on the convergence of training. From Figure 11a, it is seen that the final loss is not very sensitive to the number of Adam iterations. A low value of final loss (less than 5 × 10⁻⁵) can be achieved if the number of Adam iterations lies in the range of 1.0 × 10⁴ to 2.5 × 10⁴. From the convergence histories shown in Figure 11b, it can be seen that a much larger loss declining rate can be achieved in the L-BFGS-B optimizer in comparison with the Adam optimizer. However, since L-BFGS-B is a local optimizer, a certain number of iterations in Adam is necessary at the beginning to avoid becoming stuck on the local optimum. To keep the total number of iterations (in both optimizers) as small as possible while maintaining a low final loss, 1.0 × 10⁴ iterations in Adam before switching to L-BFGS-B is an appropriate choice. Due to the non-convexness of the loss function in this study, optimization can get stuck in local minima in a few cases with certain numbers of collocation points. It is found that the optimization can escape the local minima by adjusting the number of Adam iterations.

The effects of four weights in the loss function on the convergence of training and prediction accuracy are tested by changing one weight once a time while keeping the other three fixed. The influences of the weights on the value of the final loss and prediction error are shown in Figure 12. It is seen that choosing appropriate weights is crucial for ensuring the convergence of training and achieving high accuracy in prediction. A rule that provides guidance to the selection of weights can be easily inferred from the figure. If it is required that the low final loss is less than 1.0 × 10⁻⁴ and the relative velocity errors are less than 1%, $\lambda$ and $\gamma$ should lie in range of [0.1, 1.0], while $\alpha$ and $\beta$ should lie in the range of [2.0, 100.0].

The interpretations of the optimized ranges of weights observed in this figure are as follows. First, weights much larger than one should be assigned to the boundary losses. This is because the boundary losses are usually one order of magnitude smaller than the other two loss components in the training process. Small weights for the boundary losses can result in the domination of equation and force components in the total loss. This makes the optimization hard to converge since the boundary conditions are never properly satisfied. Second, when the weights for the PDE and force losses are too small, a prediction error rises even if the training converges very well. Here, the large discrepancy between predicted and true solutions is linked with the violations of the governing equation and force compactness condition.

The convergence of the prediction error with respect to the number of collocation points is also examined. In the refinement and coarsening process, auxiliary rectangular meshes with different resolutions are generated for arranging the collocation points. The variation of the velocity error with the number of collocation points is shown in Figure 13a. It is seen that with the increasing number of collocation points, an (almost) converged solution can be reached when the number of collocation points is larger than 7.2 × 10⁴. The relative $L_2$ velocity errors in the converged solution are in the order of 10⁻³ to 10⁻². If the number of collocation points is increased further from 7.2 × 10⁴, the prediction error almost levels off (even increases slightly). This can be explained by the fact that the prediction accuracy is ultimately limited by the optimization error in training the DNN if the number of collocation points becomes sufficiently large. The supporting evidence for this explanation can be found in Figure 13b. It is seen that the value of final loss does not change significantly if the number of collocation points is larger than 1.8 × 10⁴.

Similarly, the effect of the number of marker points on the prediction accuracy is tested by changing the number of marker points from 8 to 314 while keeping the number of collocation points fixed to 7.2 × 10⁴. From Figure 14, it is seen that the prediction accuracy is less sensitive to the number of marker points. Thus, the requirement of marker-point resolution for representing the cylinder surface is not very stringent. In this test, 125 marker points on the cylinder surface correspond to a resolution that is compatible with that of the collocation points. However, 31 marker points are sufficient for achieving an accurate prediction.

The influence of imposing a force compactness constraint inside the cylinder is also examined. In other words, force penalties are also applied to the collocation points that are in the interior of the cylinder and beyond a distance of Δ from the immersed boundary. We notice that in the mesh-based IB method used for obtaining the reference solution, half of the support in the interpolation kernel function lies inside the cylinder and spans a width of Δ. Thus, the distribution of volume force in the IB-PINN with such a constraint is supposed to be more akin to that in the mesh-based IB method. The final loss, velocity errors and drag coefficient obtained by IB-PINNs with and without such a constraint are compared in

Table 3. Final loss, velocity error and drag coefficient obtained in IB-PINNs with and without the interior force compactness constraint.

	Final Loss	${\mathit{\epsilon }}_{\mathit{u}}$(%)	${\mathit{\epsilon }}_{\mathit{v}}$(%)	${\mathit{C}}_{\mathit{d}}$
With interior constraint	2.5749 × 10 −3	3.10	6.77	1.96
Without interior constraint	3.8017 × 10 −5	0.16	0.62	2.21

. It is seen that imposing such a constraint makes the training process hard to converge and thus leads to a much larger prediction error.

3.4. Transfer Learning

In the current setting of IB-PINN, we only consider solving Navier–Stokes equations with a fixed Reynolds number. For a prediction at a different Reynolds number, the DNN needs to be re-trained. The large re-training overhead is an inherent drawback that severely limits its application in optimization problems.

Transfer learning is a powerful technique in machine learning to accelerate convergence and save training time. In transfer learning, a base network is first trained on one task. It is then used as a starting point for training the target network on a different but related task. This technique has been successfully applied to the fields of computer vision and natural language processing. In one recent work by Chen et al. [36], it was shown that for solving N-S equations at different Reynolds numbers using PINN, the re-training benefited from the transfer learning technique.

In this study, a transfer learning experiment is conduced on solving N-S equations with different Reynolds numbers using IB-PINN. We consider four tasks of solving N-S equations with Reynolds numbers Re = 10, 20, 30, 40. We first train the base DNN for the task of Re = 40. In training the other three targeted DNNs for the tasks of Re = 10, 20, 30, two approaches are adopted. In the original approach, the three DNNs are trained from scratch (with Xavier initialization). In the transfer learning approach, the three DNNs are initialized by copying the layers from the trained base DNN.

The total number of iterations at convergence, the final loss and the prediction errors in three targeted DNNs trained by using the two aforementioned approaches are compared in

Table 4. A comparison of the performances in re-training the DNNs using original and transfer learning approaches.

	Original				Tranfer Learning
	# Iterations	Final Loss	${\mathit{\epsilon }}_{\mathit{u}}$(%)	${\mathit{\epsilon }}_{\mathit{v}}$(%)	# Iterations	Final Loss	${\mathit{\epsilon }}_{\mathit{u}}$(%)	${\mathit{\epsilon }}_{\mathit{v}}$(%)
Re = 10	45011	8.7173 × 10 −4	1.62	2.50	38258	2.7832 × 10 −4	0.79	1.36
Re = 20	81341	2.1812 × 10 −4	0.47	1.23	25420	7.3278 × 10 −5	0.27	0.82
Re = 30	53915	4.7878 × 10 −5	0.19	0.78	11500	5.6311 × 10 −5	0.24	0.78

. The reference solutions for Re = 10, 20, 30 are obtained by using the conventional IB method on a very fine mesh of h = 0.005. It is observed that for predictions at different Reynolds numbers, the re-training of DNNs benefits from the transfer learning approach. The reduction in the amount of training time is substantial, especially when the Reynolds numbers in the targeted DNN and base DNN are close. This observation is consistent with the finding reported by Chen et al. [36].

4. Summary

In this paper, we propose using a direct-forcing IB method in combination with a PINN to solve steady incompressible Navier–Stokes equations. In the proposed IB-PINN, the inputs are the spatial coordinates, and the outputs are the velocity, pressure and volume force. The velocity penalties at the IB marker points and force penalties at some collocation points are incorporated into the total loss function to enforce the no-slip condition and ensure a compact distribution of the volume force.

Some differences in the implementations of the proposed IB-PINN and the mesh-based IB method are highlighted here. First, the surface force and interpolation kernel function are not used in IB-PINN. Second, in IB-PINN, the pressure and volume force are two hidden states that can be obtained automatically via imposing constraints without the extra computational cost of splitting the Navier–Stokes equations using the projection approach [29].

The laminar flow past a circular cylinder that is placed in a channel is simulated using the proposed IB-PINN method, and the result is compared with the reference solution obtained using a conventional mesh-based IB method. The comparison indicates that the results from IB-PINN and from the mesh-based IB method are in excellent agreement.

A study is carried out to investigate the influences of some parameters on the performance of IB-PINN. If the number of collocation points is fixed, the neural network with an appropriate size can achieve the best accuracy. Both the convergence of training and prediction error are found to be sensitive to the weights of different components in the loss function. It is found that with increasing number of collocation points, a converged solution can be achieved. The prediction accuracy is found to be less sensitive to the resolution of marker points for a given number of collocation points.

A transfer learning experiment is conducted on solving Navier–Stokes equations with different Reynolds numbers. It is found that the re-training DNNs can benefit from the transfer learning approach in terms of convergence acceleration.

There are several future research avenues. First, the effect of the supporting width in the force compactness constraint on the performance of IB-PINN needs to be explored systematically. Second, an adaptive refinement strategy needs to be employed to enhance the solution accuracy of IB-PINN with respect to collocation points’ distribution. Third, aiming at fluid flow problems with unsteadiness, three-dimensionality and high Reynolds numbers, the computational efficiency of IB-PINN needs to be greatly improved. We notice that the parallel PINNs developed based on the strategies of parallel in-time [37] and domain-decomposition [38] are two promising candidates for efficiently solving time-dependent PDEs. In addition, the proposed IB-PINN can also be implemented with a parametric setting in which the Reynolds number is one input [20]. Under such circumstances, re-training the DNN is not needed for predictions at different Reynolds numbers.