A simple and efficient outflow boundary condition for the incompressible Navier–Stokes equations

ABSTRACT Many researchers have proposed special treatments for outlet boundary conditions owing to lack of information at the outlet. Among them, the simplest method requires a large enough computational domain to prevent or reduce numerical errors at the boundaries. However, an efficient method generally requires special treatment to overcome the problems raised by the outlet boundary condition used. For example, mass flux is not conserved and the fluid field is not divergence-free at the outlet boundary. Overcoming these problems requires additional computational cost. In this paper, we present a simple and efficient outflow boundary condition for the incompressible Navier–Stokes equations, aiming to reduce the computational domain for simulating flow inside a long channel in the streamwise direction. The proposed outflow boundary condition is based on the transparent equation, where a weak formulation is used. The pressure boundary condition is derived by using the Navier–Stokes equations and the outlet flow boundary condition. In the numerical algorithm, a staggered marker-and-cell grid is used and temporal discretization is based on a projection method. The intermediate velocity boundary condition is consistently adopted to handle the velocity–pressure coupling. Characteristic numerical experiments are presented to demonstrate the robustness and accuracy of the proposed numerical scheme. Furthermore, the agreement of computational results from small and large domains suggests that our proposed outflow boundary condition can significantly reduce computational domain sizes.

Dirichlet and Neumann boundary conditions have been widely used for outflow conditions. A Dirichlet boundary condition with a velocity profile (e.g. the Poiseuille profile in a channel) requires pre-knowledge of the flow conditions at the outlet for specifying the velocity profile explicitly, while this naturally provides a global mass flux balance. A Neumann boundary condition does not need pre-knowledge of the flow conditions; however, it requires a large enough computational domain to prevent or reduce numerical errors at the boundaries. Because of this requirement, additional computational cost is expected for applying a Neumann boundary condition. To overcome the limitations of Dirichlet and Neumann boundary conditions at the outlet, a convective boundary condition is often used for the outflow boundary condition. Ol'shanskii and Staroverov (2000) proposed a convective boundary condition that is based on a first-order advection equation (Halpern & Schatzman, 1989) and a scaling method (Hagstrom, 1991). Typically, the advection velocity is assumed to be a constant or have a Poiseuille profile, but is somewhat arbitrary. Therefore, the advection velocity should be supplemented by a correction velocity at each time step to conserve mass flux. Nevertheless, once mass flux has been kept as a constant, the convective boundary condition becomes a fixed Dirichlet boundary condition. Another effective method derived from a wave equation is the nonreflecting boundary method (Hedstrom, 1979;Jin & Braza, 1993;Johansson, 1993;Rudy & Strikwerda, 1980;Sani & Gresho, 1994;Thompson, 1987), which performs well to minimize the spurious artifacts at the outlet. This kind of model is suitable for wake and jet flow with moderate and high Reynolds numbers. Recently, Dong et al. (2014) proposed an efficient outflow boundary condition that can allow for the influx of kinetic energy into the domain but prevent uncontrolled growth in the energy of the domain. To decouple the pressure and velocity computation, they developed a rotational velocity-correction type strategy. The proposed method is efficient and can maximize the domain truncation without adversely affecting the flow physics and enable simulations even with much higher Reynolds number.
In this paper, we propose an efficient and simple transparent boundary condition to simulate channel-type flows efficiently without trivial boundary treatment, especially for low Reynolds number. A pressure boundary condition is also proposed to match with the Navier-Stokes equation. In the numerical algorithm, a staggered marker-and-cell grid is used and temporal discretization is based on a projection method. The intermediate velocity boundary condition is consistently adopted for handling the velocity-pressure coupling. Characteristic numerical experiments are presented to demonstrate the robustness and accuracy of the numerical scheme. Furthermore, the agreement of computational results from small and large domains suggests that our proposed outflow boundary condition can significantly reduce computational domain sizes.
The paper is organized as follows: in Section 2, velocity and pressure boundary conditions for an incompressible Navier-Stokes equation are presented. Section 3 outlines the numerical method. Section 4 provides a variety of results. Finally, our conclusions are given in Section 5.

Governing equations and boundary conditions
The non-dimensional Navier-Stokes equation and the continuity equation for time-dependent incompressible viscous flows can be written as where Re = ρUL/μ is the Reynolds number, ρ and μ are fluid density and viscosity, respectively, U and L are characteristic velocity and length scales, respectively, p(x, t) is the pressure, u(x, t) = (u(x, t), v(x, t)) is the velocity, where u and v are its components in the x-and y-directions, respectively, x = (x, y) are Cartesian coordinates on the domain ⊂ R 2 , and t is the temporal variable. The domain boundary consists of ∂ = ∂ wall ∂ in ∂ out , where ∂ wall , ∂ in , and ∂ out are the wall, inlet, and outlet boundaries, respectively. Here, the velocities on the boundary are defined as is a given function, u| ∂ out is computed by an outflow boundary condition. Figure 1 shows a schematic representation of the inflow-outflow conditions. For simplicity of exposition, we shall describe the inlet and outlet in the x-direction. The other directions can be similarly defined.

Velocity boundary conditions
At the inlet, the boundary condition can be chosen based on the physical problem. In our paper, we consider inflow parallel to the wall. Thus for ∂ in , we use the following boundary condition: whereû(y, t) is a given x-directional inflow velocity function. Jin and Braza (1993) developed a nonreflecting outlet boundary condition for incompressible unsteady Navier-Stokes calculations: (4) Halpern and Schatzman (1989) proposed the convective boundary conditions where Ol'shanskii and Staroverov (2000) proposed a simple convective boundary condition: The computed velocity U is supplemented by a correction velocity to balance the mass flux, i.e.
The function U is formally chosen as constant or as a Poiseuille profile. If the outlet flow is unknown, the velocity U can be calculated by taking the integral of both sides of Equation (6): However, when the inflow is a steady flow, the constant flow rate at the inlet can lead to a zero velocity U at the outlet. Thus, the outflow boundary condition, Equation (6), will reduce to a Dirichlet boundary condition, i.e. u t = 0. It should be noted that, for unsteady flow, the outflow is clearly not constant. Sohankar et al. (1998) also proposed to set U as a convective velocity but to add a correction velocity to guarantee the balance between the inlet and outlet at each time step. The correction velocity is a constant over the entire outlet and is zero when the outflow reaches a steady state. However, the choice of the convective velocity is somewhat arbitrary. Note that, in this approach, to compare with other boundary conditions, we will perform several tests by using a convective boundary condition. The chosen convective velocity U is a constant Poiseuille profile when the outflow is in the steady state; otherwise, it will be computed by using Equation (7). The fluid flow at low Reynolds number is generally characterized by smooth and constant fluid motion because viscous forces are dominant. Therefore, we can assume that the flow is only moved by the convection and that its main direction is normal to the outflow boundary. The first assumption implies that the transparent boundary condition u t + u · ∇u = 0 is satisfied and the second assumption implies v = 0 because the outflow boundary is perpendicular to the x-direction in our paper. Therefore, to simulate the fluid flow at low Reynolds number without trivial boundary treatment, we propose to use a simple and efficient transparent boundary condition: These conditions are admissible provided that the outlet boundary is not located quite nearby the downstream. Note that in solving the incompressible Navier-Stokes equations at the inflow and outflow boundaries, both the velocity-divergence and the velocity-pressure forms should be considered to balance the mass flux. This is also a new aspect of this paper. The form (8) may suffer from some numerical difficulties, which will be described in Section 3.1. Thus, by using the divergence-free condition ∇ · u = 0, we can rewrite Equation (8) as

Pressure boundary conditions
Taking the inner product of both sides of Equation (1) with n, which is the unit normal vector to the domain boundary, we obtain n · ∇p = n · −u t − u · ∇u + 1 Re u .
On ∂ wall , we get Here we have used the no-slip boundary condition, i.e. u = 0. On ∂ in , by using Equation (10), we get Here, we have used the inflow boundary condition (3). For the outlet boundary, by substituting Equation (8) into Equation (10), we get Thus, the pressure boundary condition is directly derived by using the Navier-Stokes equation and outlet flow boundary conditions. The pressure condition is mainly balanced with the viscous stress terms at the boundary. All the pressure boundary conditions are of Neumann type. Thus, the method used with the projection method would be convenient and simple. It is noteworthy that the present form of the pressure boundary condition is more plausible because the flows in microfluidics are mostly viscous dominated. For high Reynolds number flows, the pressure boundary condition becomes a homogeneous Neumann boundary condition. This is because the viscous term (u xx + u yy )/Re is negligible.

Numerical methods
We use a staggered grid, where the pressure field is stored at cell centres and velocities at cell interfaces (see Figure 2). Let h be a uniform mesh spacing. The centre of each cell, ij , is located at ( , N x and j = 1, . . . , N y . N x and N y are the numbers of cells in the x-and y-directions, respectively. The cell edges are located at ( At the beginning of each time step, given u n , we want to find u n+1 and p n+1 by solving the following discretized equations in time with a projection method (Chorin, 1968): Here, ∇ d and d denote the centred difference approximations for the gradient and Laplacian operators, respectively: The outline of the main procedure in one time step is as follows.
Step 1. Initialize u 0 to be the divergence-free velocity field.
Step 2. Update u n+1 and n · ∇ d p n+1 on the entire boundaries. The computational details will be given in Section 3.1.
Step 3. Solve an intermediate velocity fieldũ without the pressure gradient term, Step 4. Solve the Poisson equation for the pressure, The resulting linear system of Equation (17) is solved using a multigrid method. Note that Equation (17) is derived by applying the divergence operator to Equation (18): It should be pointed out that, to find the solution for pressure in Equation (17), we should know the boundary condition of the intermediate velocity. Our proposed intermediate velocity boundary condition is described in Section 3.2.
Step 5. Update the divergence-free velocity, These steps complete one time step. For more detail, see Li, Jung, Lee, Lee, and Kim (2012). In this study, we focus on introducing efficient outflow boundary conditions and pressure boundary conditions.

Boundary conditions for velocity and pressure gradient
In this section, we will give detailed descriptions of boundary conditions for the velocity and pressure gradient on the entire boundary. First, we will consider the velocity boundary condition. The velocities u n+1 on ∂ in and ∂ wall are given as a known function and zero, respectively. The outflow boundary condition Equation (8) can be discretized with an upwind scheme as Hereū n is computed using the upwind procedure: Note that Equation (19)  = 0 for every j = 1, ..., N y . This means that, once u is zero at the boundary, u will always be zero.
To avoid this trivial case, we use the conservation form of the convective terms in Equation (9). The discretized form can be written as where we used v = 0 on the outflow boundary condition. For the pressure boundary condition, we discretize Equation (12) on ∂ in as p n+1 At the wall, without loss of generality, we only describe the pressure in the y-direction. Discretizing Equation (11), we can get p n+1 Since the no-slip boundary condition is used at the wall, we can get v i,1/2 = 0. Meanwhile, on the staggered grid, v n i,−1/2 can be simply computed as −v n i,3/2 to match the no-slip boundary condition. Thus in our approach, Equation (11) can be reduced to n · ∇p n+1 = 0.
For ∂ out , the pressure boundary condition can be discretized as Note that, since there is no information for u N x + 3 2 ,j , we use a second-order one-sided approximation to u xx .

Boundary condition for intermediate velocity
Since we have used the projection method to discretize the time-dependent incompressible Navier-Stokes equation in time, it is necessary to define the boundary condition for the intermediate velocity field. In our approach, we can directly compute u * . Since we have calculated u n+1 and n · ∇ d p in Step 2 of Section 3, we can use them straightforwardly. From Equation (19), we havẽ On ∂ wall , we get For the inflow boundary, by recalling Equation (19), we havẽ (27) Here, we have used the Dirichlet boundary condition for velocity and the homogeneous Neumann boundary condition for pressure. On the outflow boundary, we obtaiñ Thus, the intermediate velocity boundary condition is consistently adopted for the velocity-pressure coupling.
With different projection methods, the gradient pressure term (∇ d p) is added (Kim, 2005;Kim & Moin, 1985;Poux, Glockner, & Azaïez, 2011;Tseng & Huang, 2014) or subtracted (Chorin, 1968;Li, Jung, Lee, Lee, & Kim, 2012;Li, Yun, Lee, Shin, Jeong, & Kim, 2013 Next, we will briefly describe the extension of our proposed method to an implicit scheme for the momentum equation (1): Then we decompose Equation (29) into two steps: Since the boundary conditions for velocity and pressure gradient are known in Section 3.2, we can straightforwardly compute the intermediate velocity values in the domain boundary as However, in that case the time step may be restricted by numerical stability, since explicit time integration for the diffusion term is considered on the boundary. One possible way to improve numerical stability is to replace d u n by dũ in Equations (22)-(24) and substitute them into Equation (32). After that, we can combine Equations (30) and (32) and find the solutions of the intermediate velocity. Onceũ is known, we can straightforwardly compute p n+1 and u n+1 by using Equations (17) and (19), respectively. Note that, in this paper, we use Equations (14) and (15) so that we only need to solve Equation (17 ). However, if the Reynolds number is much smaller, implicit or Crank-Nicolson schemes should be considered. It should also be noted that the projection method on staggered grids has the property that an arbitrary boundary condition for the pressure gradient can be used first, and then a corresponding boundary condition can be provided for the intermediate velocity through Step 5. The corresponding boundary condition should certainly make the outflow satisfy its boundary condition.

Mass flux correction algorithm
Undesirable mass flux loss in a global or local sense is possible because of boundary singularities in a finite difference framework. The proposed outflow conditions will leverage for unexpected mass flux loss by addressing the coupling effect between pressure and velocity at the singular point (see the circle in Figure 3(a)). The centred difference approximation for the velocity in Equations (21) and (28) at the singular point is not a good choice and may lead to flux loss because of oscillation of the values. Furthermore, the present boundary conditions require a transient time to achieve mass balance, which will be discussed in the next section. Thus, we propose a simple correction algorithm to maintain the global mass flux. The idea is summarized as follows.
The following should be pointed out.
(1) Our scheme satisfies the divergence-free condition with the mass flux correction algorithm because Equations (15) and (18) are used in the projection method at every time step.
(2) Since the mass correction method also leads to spurious velocities, our proposed outflow boundary condition does not require it if the outflow is parallel to the wall (see Figure 3 (b)), because in this case the mass flux loss is much smaller, which will be shown in the next section. (3) In the case of a singular point at the outlet, the correction algorithm is applied after several flow-through time periods t 0 , which satisfies t 0 ∂ in u dy = αL x L y . Here, α is the flow-through time. In our approach, we set α = 3. This is because transient flow fields may promote numerical instability. (4) In the case without a singular point at the outlet boundary, we can use the mass flux correction algorithm only once to save transient time, i.e. α = 1.

Numerical results
In this section, to demonstrate the performance of our proposed scheme we perform several numerical experiments for parabolic flow, convergence testing, unsteady flow in a channel, backward-facing step flow, and outflow with a block on the outlet.

Parabolic flow
The first numerical experiment is performed for parabolic flow through a channel. The channel ranges from 0 to 2 in the streamwise (x) direction and from 0 to 1 in the normal (y) direction. The initial velocity field consists of random perturbations with a maximum amplitude of 2; that is, u(x, y, 0) = rand(x, y) and v(x, y, 0) = rand(x, y), where rand(x,y) is a random number between −1 and 1. The inflow velocity is set as u(0, y, t) = 8(y − y 2 ) and v(0, y, t) = 0. A no-slip boundary condition is applied at the top and bottom walls; i.e. (u, v) = (0, 0). It is well known that, at the beginning of a random velocity field, the wake gradually attains a steady state for a low Reynolds number. At steady state, the shape of the flow is the same as that of the inflow. The evolution is run up to T = 7.813 with a time step of t = 4h 2 . Here h = 1/64 and Re = 100 are used. Note that, although u 0 does not satisfy the divergence-free condition, it converges rather quickly to a divergence-free velocity field in five time steps. It should be noted that, although when we begin with u 0 the projection method used can obtain a divergence-free velocity field in one time step, it will take many more iterations to ensure that the l 2 normal error of divergence of the velocity field is lower than 1E-6. Therefore, to reduce the computational cost, we divided each time step process into five steps with several iterations. Figure 4 shows the evolution of parabolic flow. The results show that our method performs well for simulating parabolic flow. To compare the velocity at the outlet, in Figure 5(a), we show the results obtained by our proposed models with and without a flux correction algorithm, along with those of the convective model and the nonreflecting model. Here U is set as a Poiseuille profile, since the outlet flow is known in this test. Furthermore, we calculate the discrete l 2 -norm of error, which is N y j=1 (u N x + 1 2 ,j − u 1 2 j ) 2 /N y .
The errors are 8.49E-3, 8.48E-3, 1.07E-3, and 2.26E-2 for our proposed models with and without a flux correction algorithm, the convective model, and the nonreflecting  model, respectively. These results indicate that the first three models can efficiently emulate parabolic flow. The reason for the failure of the nonreflecting model is that it is derived from a wave equation and not for parabolic flow. In Figure 5(b), we illustrate the flux loss rate, which is computed by the ratio of inlet and outlet fluxes. As can be seen, the convective model and our model satisfy flux conservation, whereas the nonreflecting model does not.
Note that our correction algorithm is performed after before α = 3 is introduced. The jump in flux rate as shown in Figure 5(b) is due to the correction algorithm after three time periods. The numerical jump disappears later. It should be pointed out that, without the mass correction algorithm, our model also maintains flux conservation after more fluid exits the outlet boundary, because our outflow boundary condition and pressure boundary conditions are defined by matching with the Navier-Stokes equation. Unless otherwise mentioned, we will not use our correction algorithm for the channels, which have no singular points.

Convergence test
We next perform spatial and temporal convergence tests of the proposed method. To obtain the spatial convergence rate, we perform a number of simulations with increasingly fine grids h = 1/2 n for n = 5, 6, 7, and 8 on the domain = (0, 2) × (0, 1). The initial condition and parameters are the same as in Section 4.1. It is easy to find that, at equilibrium, the constant pressure gradient equals u yy /Re, i.e. 16/Re. Thus the equilibrium solution of parabolic flow is given by Observing the evolution, which is described in Section 4.1, we can see that the outflow reaches a steady state after t = 7.813. Thus, in this section, we also run the simulations up to T = 7.813. We define the error of the numerical solution on a grid as the discrete l 2 -norm of the difference between the numerical solution u h and the exact solution, u exact : e h ij = u h ij − u exact ij . The other terms are defined similarly. The rates of convergence are defined as the ratio of successive errors: log 2 ( e h 2 / e h/2 2 ) and log 2 ( e h ∞ / e h/2 ∞ ). Here e 2 2 is a discrete l 2 norm and is defined as Ny j=1 e 2 ij /(NxNy), and e ∞ is a discrete l ∞ norm, which is defined as e ∞ = max(|e ij |). Since we refined the spatial and temporal grids by factors of 4 and 2, respectively, the ratio of successive errors should increase by a factor of 2. The errors and rates of convergence obtained using these definitions are given in Table 1. The results suggest that the scheme is almost second-order accurate in space, as expected from the discretization.
To obtain the convergence rate for temporal discretization, we fix the space step size as h = 1/64 and choose a set of decreasing time steps t = 1. 953E-3, 9.766E-4, 4.883E-4, and 2.441E-4. The errors and rates of convergence obtained using these definitions are given in Table 2. First-order accuracy with respect to time is observed, as expected from the discretization.

Backward-facing step flow
To examine the robustness and accuracy of the proposed method, we consider a backward-facing step flow. We then compare our results with the previous experimental (Armaly, Durst, & Pereira, 1983) and numerical (Barton, 1997;Kim & Moin, 1985) results.
The agreement in our computational results obtained from a larger domain and a smaller domain suggests that our proposed outflow boundary condition is efficient. It should be pointed out that, if the computational domain is smaller, then the numerical error increases (see Figure 9(a)). This is because our proposed method is simple, we have assumed that the flow only moves by convection, and that its main direction is normal to the outflow boundary.

Comparison with Neumann boundary condition
In this section, we will compare the results obtained by our proposed and Neumann outflow boundary conditions to show the efficiency and accuracy of our proposed method. The backward-facing step flow drawn Figure 10. Snapshots of streamlines for the proposed outflow boundary (a) and a Neumann outflow boundary (b). From left to right, they are results for = (0, 4) × (0, 1) and = (0, 7) × (0, 1), respectively.  in Section 4.3 is considered. The parameters and initial condition are chosen to be the same as in Section 4.3.
The comparisons with the two boundary conditions are drawn in Figure 10. From left to right, they are results for = (0, 4) × (0, 1) and = (0, 7) × (0, 1), respectively. As can be observed from Figure 10(a), the agreement between the results computed from the larger and smaller domains is good. This suggests that our proposed outflow boundary condition is more efficient than the Neumann boundary condition.
Furthermore, we perform a number of simulations with increasing domains = (0, L x ) × (0, 1) for L x = 4, 5, 6, and 7. Since there is no closed-form analytical solution for this problem, we consider a reference numerical solution, which is obtained in a large domain = (0, 8) × (0, 1). Therefore, we can compute the l 2 errors of velocities directly for each case. These l 2 errors are given in Table 3. The results suggest that, as the domain size increases, the l 2 errors converge. These results also suggest that our proposed outflow boundary condition is more accurate than the Neumann outflow boundary because the l 2 errors obtained by our proposed outflow boundary condition are much smaller than those obtained by the Neumann outflow boundary.

Efficiency of the proposed method
To investigate the efficiency of our proposed method, we measure CPU times needed to solve the following problem: u(x, y, 0) = rand(x, y), v(x, y, 0) = rand(x, y) on the domain = (0, 4) × (0, 1) with 2 n+2 × 2 n meshes for n = 4, 5, 6, 7, and 8. The inflow velocities are u(0, y, t) = max(24(1 − y)(y − 0.5), 0) and v(0, y, t) = 0. Here, Re = 500 is fixed. All simulations are run up to 1000 t with the time step t = h 2 . Note that each simulation is run until the maximum error of the residual is less than 10 −7 . Since we refined the spatial grids by a factor of four, computation times should be increased ideally by a factor of four owing to the discretization (17), which is solved by using a multigrid method, with no special computing operators for the inflow and outflow boundary conditions added. Figure 11 shows the total CPU time versus number of mesh grid (N x N y ). The straight fitting plot implies that the multigrid solver achieves O(N x N y ) efficiency. Furthermore, the computing time for the boundary conditions can be negligible because our method updates the boundary velocity explicitly. If an implicit update for the boundary condition is chosen, then the computational time is not negligible.

Three-dimensional flow
Our method can also be straightforwardly extended to three-dimensional flows. Here, we consider a backwardfacing step flow with an initial condition of a random velocity field on the computational domain = (0, 3) × (0, 1) × (0, 1). The inflow velocities are u(0, y, z, t) =

Block on the outlet
We now consider the flow geometry shown in Figure 13(a). A fluid enters the upper half of a rectilinear channel at the left domain boundary and exits the lower half of the channel at the right boundary. The inflow velocity is set as u(0, y, 0) = max(24(1 − y)(y − 0.5), 0) and v(0, y, 0) = 0 on = (0, L x ) × (0, 1) with h = 1/64, t = 4h 2 , L x = 4, and Re = 500. A block is set at y > 0.5 at the outflow boundary. Note that, in this case, there is a singularity point on the outflow boundary. Thus we should use the proposed mass correction method. Since there is no closed-form analytical solution for this problem, we consider a reference numerical solution, which is first obtained in a large domain = (0, L x + 2) × (0, 1), as shown in Figure 13(b), and then in the domain = (0, L x ) × (0, 1). The other parameters are chosen to be the same as in the above initial condition. Figure 14 shows contour plots of the stream functions at t = 3.052 and 42.724 for convective, nonreflecting, our proposed model, and a reference numerical solution, respectively. Here, U is computed by using Equation (7), since the outlet flow is unknown in this test. As can be seen, compared with the convective and proposed models, fluids cannot completely go across the outlet boundary by using the nonreflecting model. Since its outflow boundary condition is derived from a wave equation without flux correction, unphysical flow appears for lower Reynolds numbers. Also, the convective model can simulate inflow and outflow when a block is set on the outlet boundary, whereas the outlet flow boundary is reduced to a Dirichlet boundary though the outflow does not reach a steady flow. This is not a naturally occurring physical phenomenon. Compared with them, in our proposed model the outflow seems to be a parabolic-type flow, which is closer to the reference numerical solution. The reason for this is that our model allows the outflow to be adjusted on the basis of the current inflow and inside fluid, whereas the convective model works with a given artificial velocity and mass flux correction method. Figures 15(a)-15(d) show the contour plot of the stream function with three domain sizes L x = 2, 4, and 6 for convective, nonreflecting, our proposed models, Figure 14. Evolution of outflow with a block obtained by using our proposed model. From top to bottom, they are for convective, nonreflecting, our proposed models, and the reference numerical solution, which is obtained in a larger domain as shown in Figure 13   and the reference numerical solutions, respectively. In Figure 16, we compare the outlet velocity for the four cases with the three domain sizes.
As can be seen, the outflows obtained by using the nonreflecting model lead to flux loss with the increase in the domain size. Meanwhile, with the convective model, the outflow reaches a similar constant outflow profile for different domain sizes. Compared with the nonreflecting and convective models, our proposed model reaches a natural and efficient outflow, which seems to be a parabolic-type flow, and is closer to the reference numerical solution.

Conclusions
In this paper, we presented a simple and efficient outflow boundary condition for an incompressible Navier-Stokes equation. The proposed outflow boundary condition is based on the transparent equation, where a weak formulation is used. By matching with the Navier-Stokes equation adopted inside the domain, we exactly presented the pressure boundary condition. Furthermore, we proposed the pressure boundary for convective and nonreflecting models. In the numerical algorithm, a staggered marker-and-cell grid was used and temporal discretization was based on a projection method. The intermediate velocity boundary condition is consistently adopted to handle the velocity-pressure coupling. Although our numerical scheme is an explicit scheme in time and has first-order accuracy in time, we focus on proposing efficient outflow boundary conditions and pressure boundary conditions in this study. The extensions to implicit and Crank-Nicolson schemes are straightforward. Characteristic numerical experiments were presented to demonstrate the accuracy of our proposed model. Although a first-order equation is used on the outlet boundary, almost second-order accuracy is observed for our model in the convergence test, because we use a full second-order scheme to discretize the Navier-Stokes equation inside the domain. Compared with the convective model and our proposed model, the nonreflecting model is not a good choice for low Reynolds number flow and the outflow when a block is on the outlet boundary, because the nonreflecting model cannot maintain mass flux conservation and parabolic flow is not its solution. The convective model is simple and effective satisfying mass flux conservation, although once the mass flux of the inflow becomes constant, the convective model will be reduced to a Dirichlet boundary condition. Compared with these two methods, our proposed method allows the outlet flow to be adjusted on the basis of the current inflow and fluids inside the domain, when the Reynolds number is low. Our proposed method is simple compared to that of Dong et al. (2014). However, this comparison is in some ways unfair, because the method presented in Dong et al. (2014) can maximize the domain truncation without adversely affecting the flow physics even at higher Reynolds numbers, whereas our approach aims to propose an efficient and simple transparent boundary condition to simulate channeltype flows efficiently without trivial boundary treatment at lower Reynolds numbers. When the Reynolds number is high, our method requires a much larger computational domain to reduce numerical errors compared to that of Dong et al. (2014), owing to the simplicity of our proposed transparent boundary condition. In future, we will extend our work to simulate multiphase fluids.
Mass conservation should be satisfied for each fluid separately if the mass flux is to be conserved. Therefore, the inflow and outflow boundaries should be coupled with the phases in that case. In future work, we will investigate efficient inflow and outflow boundaries to conserve the total mass of each fluid.