An Improved Immersed Boundary Method for Simulating Flow Hydrodynamics in Streams with Complex Terrains

: Three-dimensional (3D) computational ﬂuid dynamic (CFD) simulations have gained substantial popularity in recent years for stream ﬂow modelling. The complex terrain in streams is usually represented by a 3D mesh conforming to the terrain geometry. Such terrain-conforming meshes are time-consuming to generate. In this work, an immersed boundary method is developed in an existing terrain-conforming CFD model named U 2 RANS as an alternative, in which terrains are represented implicitly in the Cartesian background mesh. An improved two-layer wall function is proposed in the framework of the k - ε turbulence model, with the aim of producing accurate and smooth wall shear stress distribution and paving the way for future model development on sediment transport and scour modeling. The improvement overcomes the inherent discontinuity and nonlinearity of the two-layer velocity proﬁle, which causes error in the estimation of shear velocity. The new algorithm utilizes a distance control on the image point in immersed boundary method and a modiﬁcation of velocity prediction in the laminar layer. The improved immersed boundary method is tested with 1D, 2D, and 3D cases, and comparisons with ﬂume experiments show promising results. velocity proﬁles agree well with the experimental data for all three proﬁles and at the eight streamwise stations. The results show that the recovery from the cylinder is slow, even at the last measured location about eight cylinder diameters downstream ( x = 24 cm) where the wake e ﬀ ect of the cylinder is still signiﬁcant.


Introduction
Traditional three-dimensional (3D) modeling requires a 3D mesh to conform to domain geometry, termed the "terrain-conforming method" in this paper. Despite the widespread use of the terrain-conforming method in Computational Fluid Dynamics (CFD) models, generation of a high-quality mesh is still a challenging task in 3D modeling of flows over complex terrain. A very refined near-wall mesh must be used to resolve the boundary layer to produce an accurate solution. The stability and accuracy of the terrain-conforming method highly depend on the mesh quality. In addition, the mesh size increases rapidly with the increase of Reynolds numbers [1].
In this work, we aim to present a user-friendly 3D CFD model for practical applications in hydraulic engineering. An alternative way to the terrain-conforming method is the immersed boundary (IB) method, in which terrains are embedded in a background mesh. The boundary conditions on terrain surface are implicitly represented by modifying the governing equations of the flow. Such a method is also referred to as the terrain-embedding method in this work. Special numerical treatments are developed to implement the solid boundary conditions for turbulent flows near complex terrain and solid objects. For example, a discrete-forcing term can be added in the discretized Navier-Stokes equations to represent the embedded terrain [2,3]. The forcing term is calculated from the no-slip

Numerical Methods
U 2 RANS is a 3D CFD model using the unstructured mesh with arbitrarily shaped cells [14]. The flow is assumed incompressible, viscous, and Newtonian. The governing RANS equations are as follows: where t is time; ρ is the flow density; U is the mean velocity; τ is the turbulence stress; P is the mean pressure; υ is the kinematic viscosity; and g is the gravity acceleration. The Reynolds stress tensor τ is calculated with the standard k-ε model of [17]: where δ is the Kronecker delta (a unit tensor) and the eddy viscosity υ t is calculated as: where C µ = 0.09; k is the turbulence kinetic energy; and ε is the turbulence dissipation rate. The transport equations for k and ε are: where G = τ : ∇U is the turbulence generation rate; C ε1 = 1.44; C ε2 = 1.92; σ k = 1.0; and σ ε = 1.3. A numerical solution to the above flow equations involves the use of a mesh to cover the model domain, followed by discretization of the governing equations. With the IB method, only a background mesh is needed. The background mesh itself, however, can be unstructured and assume polyhedral shapes. Such a mesh is the most flexible and has the advantage of uniting various mesh topologies into a single formulation. The cell-centered scheme is adopted with all dependent variables located at the centroid of a mesh cell. The other alternatives include the cell-vertex scheme, with which all variables are at the cell's vertices, or a staggered scheme, with which the velocity and pressure are stored at difference locations. The governing equations are discretized using the finite-volume approach using the Gauss theorem. An advantage of the finite-volume method is that the conservation of any flow property can be achieved locally and globally. The detailed numerical method was referred to in [14], is and not repeated herein.

Immersed Boundary Implementation
The proposed IB method uses a discrete forcing approach, where the forcing term is cell-based for an unstructured mesh. The IB cells are cells cut by the immersed surface, and their cell centers are located on the fluid side (yellow cells shown in Figure 1). Three different points are identified: (1) the IB cell center (IB); (2) the hit point (HP), which is the intersection of the immersed surface with its normal line through the cell center; and (3) the image point (IP), which is the point on the extended line of the normal vector through the cell center in the fluid field. Three different characteristic lengths are identified: (1) wall distance y IB : the distance from the IB cell center to the corresponding hit point; (2) image distance y IP : the distance from the image point to the corresponding hit point; (3) and IB cell length l * : the minimum dimension of all IB cells. To enforce the turbulence model conditions at the IB cell centers, the following steps are carried out:

1.
The tangential flow velocity u IP at an image point is reconstructed by using a distance-based weighting procedure with the neighboring cells around the image point.

2.
Based on the log-law velocity profile, the dimensionless distance y + IB and y + IP are computed (iteratively) as: where E = 9.8 and κ = 0.41. 3.
The shear velocity u τ on the immersed surface is calculated as:

4.
The tangential flow velocity u IB , k IB , and ε IB at the IB cell center are calculated based on y + IB : where C k = −0.416 and B k = 8.366. The implementation of k IB and ε IB are similar to the wall functions for k and ε used in OpenFOAM [18].

5.
Fix the values of flow variables on the IB cell centers when solving the momentum equation and the transport equations for k and ε.

6.
Adjust the flux balance on the faces of IB cell centers on the solid side for mass conservation.
3. The shear velocity on the immersed surface is calculated as: 4. The tangential flow velocity , , and at the IB cell center are calculated based on : = + ,  The key to the IB treatment is the estimation of the shear velocity , that is, the calculation of . Considering the nonlinear and discontinuous nature of velocity profile between the laminar sublayer and log-law sublayer, the result of converges to different values if different velocity profiles are used. To keep the consistency when estimating the shear velocity of all IB cells, Equation (7) assumes that the image points are located within the log-law sublayer such that only the log-law velocity profile is used in the iteration. Generally, the image point distance is proportional to the wall distance (for example, = 3 ), as shown in Figure 1a [19,20]. However, the wall distance is arbitrarily distributed around the immersed surface, making it impossible to guarantee that the image point is located in the logarithmic sublayer. In addition, an extremely small value of may result in numerical instability and even divergence of the model. In this work, is set to be proportional to the minimum dimension * of each IB cell ( = 3 * ). Consequently, the image points are uniformly distributed along the immersed surface, as shown in The key to the IB treatment is the estimation of the shear velocity u τ , that is, the calculation of y + IP . Considering the nonlinear and discontinuous nature of velocity profile between the laminar sublayer and log-law sublayer, the result of y + IP converges to different values if different velocity profiles are used. To keep the consistency when estimating the shear velocity u τ of all IB cells, Equation (7) assumes that the image points are located within the log-law sublayer such that only the log-law velocity profile is used in the iteration. Generally, the image point distance y IP is proportional to the wall distance y IB (for example, y IP = 3y IB ), as shown in Figure 1a [19,20]. However, the wall distance y IB is arbitrarily distributed around the immersed surface, making it impossible to guarantee that the image point is located in the logarithmic sublayer. In addition, an extremely small value of y IP may result in numerical instability and even divergence of the model. In this work, y IP is set to be proportional to the minimum dimension l * of each IB cell (y IP = 3l * ). Consequently, the image points are uniformly distributed along the immersed surface, as shown in Figure 1b. The numerical instability due to the small value of y IP is avoided. A similar implementation has been used in [12] and [11] for the k-ω model and SA model.
Another problem in the near-wall treatment of the immersed boundary method comes from the inconsistency when the IB cell center and its corresponding image point are located in different sublayers of the velocity profile. Although all image points are designed to be in the log-law layer, the wall distance y IB is still an arbitrary value and the IB cell center may be in the laminar sublayer. Therefore, Equation (9), which assumes the IB cell center and its image point follow the same logarithmic velocity profile, is not applicable in this situation. To address this problem, a modified velocity profile is used when the IB cell center is in the laminar sublayer. We assume a linear velocity profile between the wall and the image point, such that the first derivative of the velocity with respect to the wall distance is still a constant (∂u/∂y = const.) (Figure 2a). This method was proposed in [11] to correct the mass flux on the cell boundary by using a slip velocity boundary condition. Here, it is used to extend the logarithmic velocity profile to the wall when the IB cell center is in the laminar sublayer. The eddy viscosity is also modified to be a constant between the wall and the image point ( Figure 2b). Thus, k IB , and ε IB are constant in this region. This modification is based on the balance of shear stress on the boundary: Water 2020, 12, x FOR PEER REVIEW 5 of 14 Another problem in the near-wall treatment of the immersed boundary method comes from the inconsistency when the IB cell center and its corresponding image point are located in different sublayers of the velocity profile. Although all image points are designed to be in the log-law layer, the wall distance is still an arbitrary value and the IB cell center may be in the laminar sublayer. Therefore, Equation (9), which assumes the IB cell center and its image point follow the same logarithmic velocity profile, is not applicable in this situation. To address this problem, a modified velocity profile is used when the IB cell center is in the laminar sublayer. We assume a linear velocity profile between the wall and the image point, such that the first derivative of the velocity with respect to the wall distance is still a constant ( ⁄ = const. ) (Figure 2a). This method was proposed in [11] to correct the mass flux on the cell boundary by using a slip velocity boundary condition. Here, it is used to extend the logarithmic velocity profile to the wall when the IB cell center is in the laminar sublayer. The eddy viscosity is also modified to be a constant between the wall and the image point ( Figure 2b). Thus, , and are constant in this region. This modification is based on the balance of shear stress on the boundary: The modified flow velocity , , and at the IB cell center are calculated as the following: In the following, this is called the modified wall model, and the one presented before is called the original model. The modified flow velocity u IB , k IB , and ε IB at the IB cell center are calculated as the following: where ∂u + /∂y + = 1/(κy + ). In the following, this is called the modified wall model, and the one presented before is called the original model.

Result
The above IB method is implemented into the U 2 RANS model. In this section, a number of turbulent flow cases are selected to verify the improved IB method. In particular, model accuracy is examined and discussed.

Turbulent Flow over Flat Plate
The flow over a flat plate was used to investigate a 2D boundary layer without pressure gradient. The case is used to verify the IB method implementation with the proposed wall function. The length L of the plate is 2 m. The Reynolds number based on the plate length is Re L = 5 × 10 6 and incoming flow velocity is U = 2.5 m/s. The upper boundary is 0.1L away from the flat plate. The Cartesian mesh is refined near the leading edge and the plate. The mesh arrangement is shown in Figure 3. A short slip surface (0.1L) is added in front of the leading edge using the symmetric boundary condition. The flat plate is modeled as an immersed boundary, and the wall boundary condition is applied on it. We tested four different mesh sizes by using four different cell expansion ratios, R y = 1, 5, 10, 20, in the y-direction to change the mesh size of the first cell touching the wall. The cell expansion ratio, R, is that of the size of the end cell ∆ x N to the size of the start cell ∆ x 1 along the edge direction (R = ∆ x N /∆ x 1 ). Different mesh sizes change the value of y IB and y IP . For both original and modified wall models, the computed local friction coefficient (C f ) along the plate compares well with the experimental data from [21] (Figure 4a,b); the velocity profiles at 0.9L are in good agreement with the log-law (Figure 4c,d).
In addition, as y + IB decreases, the simulation results converge to the experimental data. The results show that the proposed IB algorithm with an original or modified wall model is insensitive to the image point distance y IP in the log-law layer.
To further investigate the stability of the algorithm and its dependence on wall distance, the mesh above the flat plate is rotated, as shown in Figure 5. The top line of the grid has the same height to maintain the same water depth throughout the length. The bottom of the mesh is rotated, and the downstream end of the bottom is moved down by 0.0025L. The red line represents the position of the flat plate. With this configuration, the mesh lines are not aligned with the flat plate. As the grid rotates (stretches), the wall distance is not uniformly distributed along the plate. Other details of the mesh arrangement are shown in Figure 5. N y2 is the number of mesh refinement in the y-direction used in the region between the bottom line and the line 0.1L away from the bottom. We changed the number of N y2 (N y2 = 100, 200, 300) to show the mesh independence of this method. The cell expansion ratio in the y-direction in the refinement zone is 1 to maintain the minimum dimension l * constant for each IB cell, such that the image point distance is the same; but the wall distance distribution is nonuniform over the plate. Figure 6a shows the numerical results of the local friction coefficient using the original wall model. The predicted results are comparable with the experiment. However, there are some small, semi-periodic oscillations due to the change of wall distance along the plate, especially when the value of y IB decreases and the IB cell center is in the laminar sublayer. To suppress the oscillation from the nonuniform wall distance, the modified wall model is applied and the simulated results are plotted in Figure 6b. The modified wall model can predict the local friction coefficient more accurately and the oscillation is greatly reduced. As the mesh is refined, the local friction coefficient converges to the experimental data. In addition, velocity profiles in Figure 6b and d at 0.9L agree well with the log-law. Both the original and modified wall models give a good prediction of velocity.
the computed local friction coefficient ( ) along the plate compares well with the experimental data from [21] (Figure 4a,b); the velocity profiles at 0.9L are in good agreement with the log-law ( Figure  4c,d). In addition, as decreases, the simulation results converge to the experimental data. The results show that the proposed IB algorithm with an original or modified wall model is insensitive to the image point distance in the log-law layer.   To further investigate the stability of the algorithm and its dependence on wall distance, the mesh above the flat plate is rotated, as shown in Figure 5. The top line of the grid has the same height to maintain the same water depth throughout the length. The bottom of the mesh is rotated, and the downstream end of the bottom is moved down by 0.0025L. The red line represents the position of the flat plate. With this configuration, the mesh lines are not aligned with the flat plate. As the grid rotates (stretches), the wall distance is not uniformly distributed along the plate. Other details of the mesh arrangement are shown in Figure 5.
is the number of mesh refinement in the y-direction used in the region between the bottom line and the line 0.1L away from the bottom. We changed the number of ( = 100, 200, 300) to show the mesh independence of this method. The cell expansion ratio in the y-direction in the refinement zone is 1 to maintain the minimum dimension * constant for each IB cell, such that the image point distance is the same; but the wall distance distribution is nonuniform over the plate. Figure 6a shows the numerical results of the local friction coefficient using the original wall model. The predicted results are comparable with the experiment. However, there are some small, semi-periodic oscillations due to the change of wall distance along the plate, especially when the value of decreases and the IB cell center is in the laminar sublayer. To suppress the oscillation from the nonuniform wall distance, the modified wall model is applied and the simulated results are plotted in Figure 6b. The modified wall model can predict the local friction coefficient more accurately and the oscillation is greatly reduced. As the mesh is refined, the local friction coefficient converges to the experimental data. In addition, velocity profiles in Figure 6b and d at 0.9L agree well with the log-law. Both the original and modified wall models give a good       Table 1 provides the computational cost in one iteration of the IB method with the modified wall model. The calculation of the modified wall model includes the finding of IB cells, identification of image stencil for reconstruction, and implementation of wall function. All simulations were performed on a Dell Precision Tower 5810 with an Intel Xeon CPU ES-1620. The computational cost of the modified wall model is about 60% of the total CPU time in one iteration for each case. However, considering the boundary is immobile in this work, the IB cells and image stencil are only calculated in the first iteration, and remain unchanged in the following iterations. The percentage of the computational cost of the modified wall model greatly decreases over the whole simulation time. In addition, the modified wall model avoids the computational cost required to fully resolve the near-wall flow in the laminar layer.

Turbulent Flow around a Cylinder over Scoured Beds
The proposed IB algorithm is verified next for its capability to simulate a case with an instream structure-a turbulent flow around an instream cylinder over scoured beds. The modified wall model is used. The simulated results are compared with the flume experiment by Jensen et al. [22]. Figure 7 shows the three bed profiles representing three scour phases observed in the experiment by Mao [23]. The cylinder has a diameter of 3 cm and is placed above three scoured bed profiles. The mean inlet The proposed IB algorithm is verified next for its capability to simulate a case with an instream structure-a turbulent flow around an instream cylinder over scoured beds. The modified wall model is used. The simulated results are compared with the flume experiment by Jensen et al. [22]. Figure 7 shows the three bed profiles representing three scour phases observed in the experiment by Mao [23]. The cylinder has a diameter of 3 cm and is placed above three scoured bed profiles. The mean inlet flow velocity is 0.2 m/s. The computational domain is 1.1 m in the streamwise direction (x-direction). The flow depth at the uneroded bed at the exit is maintained at 0.245 m (y-direction). Despite 2D flow in nature, the modeling is carried out in a 3D model domain with the dimension of 0.03 m along the cylinder (z-direction).

Figure 7.
Three cases are simulated corresponding to three bed profiles; the top, middle, and bottom profiles correspond to profiles 1, 3, and 5, respectively, of the experiments in [22]. The simulated pressure contours are also displayed.
Using profile 3 as an example, the Cartesian background mesh has 0.1 cm resolution near the cylinder and scour bed, which is the finest resolution similarly used by Smith and Foster [24]. The background mesh has a total of 355,515 3D cells (71,103 2D cells in the xy plane and 5 cells along z). The cylinder boundary and bed profiles are treated as the immersed boundaries. The fluid and IB cells near the immersed boundaries are shown in Figure 8. The meshes for the two other profiles are similar.  Figure 9 shows the comparisons between simulation results and measured data from experiments. It is seen that the computed streamwise (x) velocity profile approaching the cylinder is near-logarithmic, although a constant velocity boundary condition is applied at the inlet. The IB method-simulated velocity profiles agree well with the experimental data for all three profiles and at the eight streamwise stations. The results show that the recovery from the cylinder is slow, even at the last measured location about eight cylinder diameters downstream (x = 24 cm) where the wake effect of the cylinder is still significant.  [22]. The simulated pressure contours are also displayed.
Using profile 3 as an example, the Cartesian background mesh has 0.1 cm resolution near the cylinder and scour bed, which is the finest resolution similarly used by Smith and Foster [24]. The background mesh has a total of 355,515 3D cells (71,103 2D cells in the xy plane and 5 cells along z). The cylinder boundary and bed profiles are treated as the immersed boundaries. The fluid and IB cells near the immersed boundaries are shown in Figure 8. The meshes for the two other profiles are similar. structure-a turbulent flow around an instream cylinder over scoured beds. The modified wall model is used. The simulated results are compared with the flume experiment by Jensen et al. [22]. Figure 7 shows the three bed profiles representing three scour phases observed in the experiment by Mao [23]. The cylinder has a diameter of 3 cm and is placed above three scoured bed profiles. The mean inlet flow velocity is 0.2 m/s. The computational domain is 1.1 m in the streamwise direction (x-direction). The flow depth at the uneroded bed at the exit is maintained at 0.245 m (y-direction). Despite 2D flow in nature, the modeling is carried out in a 3D model domain with the dimension of 0.03 m along the cylinder (z-direction).

Figure 7.
Three cases are simulated corresponding to three bed profiles; the top, middle, and bottom profiles correspond to profiles 1, 3, and 5, respectively, of the experiments in [22]. The simulated pressure contours are also displayed.
Using profile 3 as an example, the Cartesian background mesh has 0.1 cm resolution near the cylinder and scour bed, which is the finest resolution similarly used by Smith and Foster [24]. The background mesh has a total of 355,515 3D cells (71,103 2D cells in the xy plane and 5 cells along z). The cylinder boundary and bed profiles are treated as the immersed boundaries. The fluid and IB cells near the immersed boundaries are shown in Figure 8. The meshes for the two other profiles are similar.  Figure 9 shows the comparisons between simulation results and measured data from experiments. It is seen that the computed streamwise (x) velocity profile approaching the cylinder is near-logarithmic, although a constant velocity boundary condition is applied at the inlet. The IB method-simulated velocity profiles agree well with the experimental data for all three profiles and at the eight streamwise stations. The results show that the recovery from the cylinder is slow, even at the last measured location about eight cylinder diameters downstream (x = 24 cm) where the wake effect of the cylinder is still significant.  Figure 9 shows the comparisons between simulation results and measured data from experiments. It is seen that the computed streamwise (x) velocity profile approaching the cylinder is near-logarithmic, although a constant velocity boundary condition is applied at the inlet. The IB method-simulated velocity profiles agree well with the experimental data for all three profiles and at the eight streamwise stations. The results show that the recovery from the cylinder is slow, even at the last measured location about eight cylinder diameters downstream (x = 24 cm) where the wake effect of the cylinder is still significant.

Turbulent Flow over 3D Dunes
3D dunes were used to verify the model performance with the turbulence model in a 3D form, especially the prediction of wall shear stress. The modified wall model is used in this case for a more accurate and smooth wall shear prediction. The simulation results are compared with the experiments of [25]. In the experiment, 14 fixed 3D dunes were placed on the bottom of the flume sequentially and experimental data were measured on the 11th and 12th numbered dunes. Only six dunes were simulated to reduce the computational cost, and the data were collected from the last two dunes, as shown in Figure 10a. The length of the simulation domain is 5.0 m in the x-direction and the width of the domain is 0.9 m in the y-direction. The bed elevation contour of the 3D dunes is shown in Figure 10b. The incoming velocity is 0.261 m/s in the x-direction and the water depth is 0.561 m in the z-direction.

Turbulent Flow over 3D Dunes
3D dunes were used to verify the model performance with the turbulence model in a 3D form, especially the prediction of wall shear stress. The modified wall model is used in this case for a more accurate and smooth wall shear prediction. The simulation results are compared with the experiments of [25]. In the experiment, 14 fixed 3D dunes were placed on the bottom of the flume sequentially and experimental data were measured on the 11th and 12th numbered dunes. Only six dunes were simulated to reduce the computational cost, and the data were collected from the last two dunes, as shown in Figure 10a. The length of the simulation domain is 5.0 m in the x-direction and the width of the domain is 0.9 m in the y-direction. The bed elevation contour of the 3D dunes is shown in Figure 10b. The incoming velocity is 0.261 m/s in the x-direction and the water depth is 0.561 m in the z-direction.
The Cartesian background mesh is shown in Figure 11. The numbers of cells are 450, 90, and 70 in the x, y, and z directions, respectively. The mesh is refined near the dunes, especially where the bed elevation changes rapidly.
dunes were simulated to reduce the computational cost, and the data were collected from the last two dunes, as shown in Figure 10a. The length of the simulation domain is 5.0 m in the x-direction and the width of the domain is 0.9 m in the y-direction. The bed elevation contour of the 3D dunes is shown in Figure 10b. The incoming velocity is 0.261 m/s in the x-direction and the water depth is 0.561 m in the z-direction. The Cartesian background mesh is shown in Figure 11. The numbers of cells are 450, 90, and 70 in the x, y, and z directions, respectively. The mesh is refined near the dunes, especially where the bed elevation changes rapidly.  Figure 12 is a comparison of results from the IB method and the experiment. The streamwise velocities at different locations are well-predicted. The only noticeable deviation is downstream of the measured two dunes at y = 0 m. In this slice, the bed elevation has the largest slope, such that a long distance from the inlet is required for the flow to be fully developed. In the experiment, the measured dunes are the eleventh and twelfth. However, the measured dunes in the simulation are at the fifth and sixth, due to the limitation of computing capacity. Thus, the flow condition is slightly different.
Comparison of bed shear stress is shown in Figure 13. The simulation results show that the new IB algorithm can provide a smooth wall shear stress distribution. The normalized wall shear stress of the measured data (Figure 13a) is estimated using the velocity 5 mm above the bed [25]: where is the velocity at 5 mm above the bed, the distance − is 5 mm, and = 1 mm is the bed roughness.
The wall shear stress from the IB simulation is based on Equations (7) and (8) and shown in Figure 13b. Even with the distribution of IB cells and image points changing arbitrarily with the bed elevation, the predicted shear stress is smooth. A comparison with the experimental data in Figure  13a shows the existence of mismatch between the two. This is mainly because the accuracy of the simulation is limited by the computing capacity, such that the velocity prediction at the places with large slope deviates from the experiment. However, the general characteristics of the wall shear stress are captured by the numerical model. The wall shear stress is the highest at the crests of dunes and relatively small elsewhere. The smooth distribution of wall shear stress is very important for modeling the sediment transport and scours, which is currently under development. The Cartesian background mesh is shown in Figure 11. The numbers of cells are 450, 90, and 70 in the x, y, and z directions, respectively. The mesh is refined near the dunes, especially where the bed elevation changes rapidly.  Figure 12 is a comparison of results from the IB method and the experiment. The streamwise velocities at different locations are well-predicted. The only noticeable deviation is downstream of the measured two dunes at y = 0 m. In this slice, the bed elevation has the largest slope, such that a long distance from the inlet is required for the flow to be fully developed. In the experiment, the measured dunes are the eleventh and twelfth. However, the measured dunes in the simulation are at the fifth and sixth, due to the limitation of computing capacity. Thus, the flow condition is slightly different.
Comparison of bed shear stress is shown in Figure 13. The simulation results show that the new IB algorithm can provide a smooth wall shear stress distribution. The normalized wall shear stress of the measured data (Figure 13a) is estimated using the velocity 5 mm above the bed [25]: where is the velocity at 5 mm above the bed, the distance − is 5 mm, and = 1 mm is the bed roughness.
The wall shear stress from the IB simulation is based on Equations (7) and (8) and shown in Figure 13b. Even with the distribution of IB cells and image points changing arbitrarily with the bed elevation, the predicted shear stress is smooth. A comparison with the experimental data in Figure  13a shows the existence of mismatch between the two. This is mainly because the accuracy of the simulation is limited by the computing capacity, such that the velocity prediction at the places with large slope deviates from the experiment. However, the general characteristics of the wall shear stress are captured by the numerical model. The wall shear stress is the highest at the crests of dunes and relatively small elsewhere. The smooth distribution of wall shear stress is very important for modeling the sediment transport and scours, which is currently under development.  Figure 12 is a comparison of results from the IB method and the experiment. The streamwise velocities at different locations are well-predicted. The only noticeable deviation is downstream of the measured two dunes at y = 0 m. In this slice, the bed elevation has the largest slope, such that a long distance from the inlet is required for the flow to be fully developed. In the experiment, the measured dunes are the eleventh and twelfth. However, the measured dunes in the simulation are at the fifth and sixth, due to the limitation of computing capacity. Thus, the flow condition is slightly different.  Comparison of bed shear stress is shown in Figure 13. The simulation results show that the new IB algorithm can provide a smooth wall shear stress distribution. The normalized wall shear stress of the measured data (Figure 13a) is estimated using the velocity 5 mm above the bed [25]: where U 1 is the velocity at 5 mm above the bed, the distance z 1 − η is 5 mm, and k s = 1 mm is the bed roughness.

Conclusions
This work proposes a new wall model algorithm for use by the IB method in association with the − turbulence model. It was implemented into the terrain-conforming U 2 RANS model of [14]. The aim was to develop a flexible IB-method-based CFD model so that turbulent flow field can be simulated accurately and smooth wall shear stress can be generated. A key contribution is to use a consistent velocity profile-log-law velocity profile-to estimate the wall shear velocity. In addition, a modification was proposed in the wall model for the laminar sublayer to suppress the oscillation caused by the small value of wall distance . In such a way, the discontinuity and nonlinearity of the velocity profile between the laminar sublayer and log-law sublayer are avoided. The proposed two-layer wall function has the same drawbacks of standard wall functions for non-equilibrium flow Figure 13. (a) Normalized wall shear stress from [26], |τ T | = 0.00046 m 2 /s; (b) normalized wall shear stress from the IB method.
The wall shear stress u 2 τ from the IB simulation is based on Equations (7) and (8) and shown in Figure 13b. Even with the distribution of IB cells and image points changing arbitrarily with the bed elevation, the predicted shear stress is smooth. A comparison with the experimental data in Figure 13a shows the existence of mismatch between the two. This is mainly because the accuracy of the simulation is limited by the computing capacity, such that the velocity prediction at the places with large slope deviates from the experiment. However, the general characteristics of the wall shear stress are captured by the numerical model. The wall shear stress is the highest at the crests of dunes and relatively small elsewhere. The smooth distribution of wall shear stress is very important for modeling the sediment transport and scours, which is currently under development.

Conclusions
This work proposes a new wall model algorithm for use by the IB method in association with the k-ε turbulence model. It was implemented into the terrain-conforming U 2 RANS model of [14]. The aim was to develop a flexible IB-method-based CFD model so that turbulent flow field can be simulated accurately and smooth wall shear stress can be generated. A key contribution is to use a consistent velocity profile-log-law velocity profile-to estimate the wall shear velocity. In addition, a modification was proposed in the wall model for the laminar sublayer to suppress the oscillation caused by the small value of wall distance y IB . In such a way, the discontinuity and nonlinearity of the velocity profile between the laminar sublayer and log-law sublayer are avoided. The proposed two-layer wall function has the same drawbacks of standard wall functions for non-equilibrium flow with a large pressure gradient, such as stagnation point and separated flow. However, it provides a relatively accurate prediction of a high Reynolds number and wall-bounded flow with an affordable computational cost. The new method was tested and validated with selected 1D, 2D, and 3D flow cases. Good results were obtained in each case. The proposed method produces an accurate and smooth wall shear stress distribution, paving the way for the next stage of model development for sediment transport and scour modeling.