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

Abstract

Three-dimensional (3D) computational fluid dynamic (CFD) simulations have gained substantial popularity in recent years for stream flow 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²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 profile, 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 modification of velocity prediction in the laminar layer. The improved immersed boundary method is tested with 1D, 2D, and 3D cases, and comparisons with flume experiments show promising results.

1. 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 condition and the moving velocity of the immersed surface if the boundary motion is considered. This method is referred to as the Kinematic IB scheme [4].

With the IB method, mesh generation becomes relatively simple, as only a background mesh is needed, and mesh quality is easy to control. A key issue of the IB method is the model accuracy in simulating turbulent flows, although it has been demonstrated for laminar flows. Another advantage of the IB method is that it can easily track terrain deformation, which is important in the simulation of sediment transport and scour. With the terrain-conforming methods, the mesh has to be regenerated to conform to the terrain or solid body movement in each morphological time step. The mesh regeneration is difficult to perform and may lead to numerical instability and loss of accuracy due to the deterioration of the mesh quality caused by large movement. With the IB method, terrain deformation is captured by recalculating the discrete-forcing term and no mesh changes are necessary.

A key variable in sediment transport and scour simulation is the wall shear stress on erodible beds. Existing algorithms of the IB method have demonstrated good performance in the prediction of flow velocity [5,6]. However, the accuracy and smoothness of the wall shear stress in the IB method need improvements. This issue is mainly caused by the poor prediction of the velocity profile between the near-wall cell and the terrain boundary [7]. The other reason is the mass conservation problem on the terrain boundary [7]. In addition, the implementation of a wall function in the IB method is a significant challenge, especially for high Reynolds number flows. To avoid high mesh resolution in the near-wall region, attempts have been made in the implementation of wall models in the large eddy simulation (LES) framework [8,9]. For example, Roman et al. proposed a Reynolds-Averaged Navier-Stokes (RANS)-like eddy viscosity to reconstruct the wall shear stress in LES [8].

Considering that there are drastically different time scales in hydrodynamics and morphodynamics (the latter is much slower than the former) [10], a wall model in the RANS framework is more suitable for real-world applications. Tamaki et al. proposed a Spalart-Allmaras (SA) wall model using the IB method, and modified the eddy viscosity to balance the shear stress [11]. Capizzano adopted a two-layer wall model in the near-wall region for k-ω and k-g turbulence models with a blending estimation of $\omega$ and $g$ [12]. Zhou applied a blending of Spalding’s law, Reichardt’s law, and log-law in the SST k-ω turbulence model [13]. The key of the most existing wall functions is to avoid the transition of the wall model between the viscous sublayer and the log-law sublayer by either using the blending turbulence variables or using continuous velocity profiles.

Despite the progresses made in the IB methodology and its applications in river hydraulics, mature models are yet to be developed to have the accuracy, stability, and efficiency needed for real-world applications. In particular, a reliable sediment transport and scour model would hinge on the accurate prediction of the bed shear stress. In this work, we present an improved IB algorithm which is then implemented into the terrain-conforming model of [14] (named U²RANS). The aim is to develop an IB-method based 3D model so that accurate flow and bed shear stress may be predicted. This paper extends the work presented in [15,16]. Specifically, a two-layer wall function is implemented in association with the k-ε turbulence model [17] and the IB method. The model is then verified to produce an accurate velocity field, as well as a smooth and accurate wall shear stress distribution. A key contribution of the present study is the improvement of the model accuracy through the control of the image point position, which is a point used in the IB model to reconstruct the near-wall flow field. In addition, a modification of the wall model in the laminar sublayer allows better handling of the discontinuity between the viscous sublayer and the logarithmic layer. As a result, a better estimation of wall shear stress can be achieved.

The new IB-method-based U²RANS model is validated and demonstrated using a number of flow cases in the following. As discussed before, due to the limitations of the terrain-conforming method, the terrain-embedding method is more suitable for modeling of real-world cases, particularly for sediment transport and scour modeling. The scour modeling is under development and will be reported in the future.

2. Numerical Methods

U²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:

$$\nabla \cdot \mathit{U}=0,$$

(1)

$$\frac{\partial \mathit{U}}{\partial t}+\left(\mathit{U}\cdot \nabla \right)\mathit{U}=\upsilon {\nabla }^2\mathit{U}+\nabla \cdot \mathit{\tau }-\frac{\nabla P}{\rho }+\mathit{g},$$

(2)

where $t$ is time; $\rho$ is the flow density; $\mathit{U}$ is the mean velocity; $\mathit{\tau }$ is the turbulence stress; $P$ is the mean pressure; $\upsilon$ is the kinematic viscosity; and $\mathit{g}$ is the gravity acceleration.

The Reynolds stress tensor $\mathit{\tau }$ is calculated with the standard k-ε model of [17]:

$$\mathit{\tau }={\upsilon }_t\left(\nabla \mathit{U}+{\left(\nabla \mathit{U}\right)}^{\mathrm{T}}\right)-\frac{2}{3}k\mathit{\delta },$$

(3)

where $\mathit{\delta }$ is the Kronecker delta (a unit tensor) and the eddy viscosity ${\upsilon }_t$ is calculated as:

$${\upsilon }_t=C_{\mu }\frac{k^2}{\epsilon },$$

(4)

where $C_{\mu }=0.09$; $k$ is the turbulence kinetic energy; and $\epsilon$ is the turbulence dissipation rate.

The transport equations for $k$ and $\epsilon$ are:

$$\frac{\partial k}{\partial t}+\nabla \cdot \left(k\mathit{U}\right)=\nabla \cdot \left(\left(\upsilon +\frac{{\upsilon }_t}{{\sigma }_k}\right)\nabla k\right)+G-\epsilon ,$$

(5)

$$\frac{\partial \epsilon }{\partial t}+\nabla \cdot \left(\epsilon \mathit{U}\right)=\nabla \cdot \left(\left(\upsilon +\frac{{\upsilon }_t}{{\sigma }_{\epsilon }}\right)\nabla \epsilon \right)+C_{\epsilon 1}\frac{\epsilon }{k}G-C_{\epsilon 2}\frac{{\epsilon }^2}{k},$$

(6)

where $G=\mathit{\tau }:\nabla \mathit{U}$ is the turbulence generation rate; $C_{\epsilon 1}=1.44$; $C_{\epsilon 2}=1.92$; ${\sigma }_k=1.0$; and ${\sigma }_{\epsilon }=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.

3. 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 ${\mathit{y}}_{\mathit{I}\mathit{B}}$: the distance from the IB cell center to the corresponding hit point; (2) image distance ${\mathit{y}}_{\mathit{I}\mathit{P}}$: the distance from the image point to the corresponding hit point; (3) and IB cell length ${\mathit{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:

• 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.
• Based on the log-law velocity profile, the dimensionless distance $y_{IB}^{+}$ and $y_{IP}^{+}$ are computed (iteratively) as:

  $$y_{IP,n+1}^{+}=\frac{y_{IP,n}^{+}+\frac{\kappa u_{IP}{y}_{IP}}{\upsilon }}{1+\ln \left(E{y}_{IP,n}^{+}\right)}; y_{IB}^{+}=y_{IP}^{+}\frac{y_{IB}}{y_{IP}},$$

  (7)

  where $E=9.8$ and $\kappa =0.41$.
• The shear velocity $u_{\tau }$ on the immersed surface is calculated as:

  $$u_{\tau }=\frac{y_{IP}^{+}\upsilon }{y_{IP}},$$

  (8)
• The tangential flow velocity $u_{IB}$, $k_{IB}$, and ${\epsilon }_{IB}$ at the IB cell center are calculated based on $y_{IB}^{+}$:

  $$u_{IB}^{new}=u_{IP}-\frac{1}{\kappa }{u}_{\tau }\log \left(\frac{y_{IP}}{y_{IB}}\right),$$

  (9)

  $$k_{IB}=\frac{C_k\log \left(y_{IB}^{+}\right)}{\kappa }+B_k,$$

  (10)

  $${\epsilon }_{IB}=\frac{C_{\mu }^{3/4}{k}_{IB}{}^{3/2}}{\kappa y_{IB}},$$

  (11)

  where $C_k=-0.416$ and $B_k=8.366$. The implementation of $k_{IB}$ and ${\epsilon }_{IB}$ are similar to the wall functions for $k$ and $\epsilon$ used in OpenFOAM [18].
• Fix the values of flow variables on the IB cell centers when solving the momentum equation and the transport equations for $k$ and $\epsilon$.
• Adjust the flux balance on the faces of IB cell centers on the solid side for mass conservation.

The key to the IB treatment is the estimation of the shear velocity $u_{\tau }$, 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_{\tau }$ 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}$ = 3$y_{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}$ = 3$l_{*}$). 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 $\left(\partial u/\partial y=\mathrm{const}.\right)$ (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 ${\epsilon }_{IB}$ are constant in this region. This modification is based on the balance of shear stress on the boundary:

$$\left(\upsilon +{\upsilon }_t\right)\left(\frac{\partial u}{\partial y}\right)=u_{\tau }^2=\mathrm{const}.$$

(12)

The modified flow velocity $u_{IB}$, $k_{IB}$, and ${\epsilon }_{IB}$ at the IB cell center are calculated as the following:

$$u_{IB}^{new}=\left\{\begin{array}{ll}{u}_{IP}-\frac{1}{\kappa }{u}_{\tau }\log \left(\frac{y_{IP}}{y_{IB}}\right)& \mathrm{if} y_{IB}^{+}>y_{Laminar}^{+}\\ u_{IP}-{\left(\frac{\partial u^{+}}{\partial y^{+}}\right)}_{IP}\left(y_{IP}-y_{IB}\right) u_{\tau }& \mathrm{if} y_{IB}^{+}\le y_{Laminar}^{+}\end{array}\right.,$$

(13)

$$k_{IB}=\left\{\begin{array}{c}\frac{C_k\log \left(y_{IB}^{+}\right)}{\kappa }+B_k \mathrm{if} y_{IB}^{+}>y_{Laminar}^{+} \\ \frac{C_k\log \left(y_{IP}^{+}\right)}{\kappa }+B_k \mathrm{if} y_{IB}^{+}\le y_{Laminar}^{+}\end{array}\right.,$$

(14)

$${\epsilon }_{IB}=\left\{\begin{array}{c}\frac{C_{\mu }^{3/4}{k}_{IB}{}^{3/2}}{\kappa y_{IB}} \mathrm{if} y_{IB}^{+}>y_{Laminar}^{+} \\ \frac{C_{\mu }^{3/4}{k}_{IP}{}^{3/2}}{\kappa y_{IP}} \mathrm{if} y_{IB}^{+}\le y_{Laminar}^{+}\end{array}\right.,$$

(15)

where $\partial u^{+}/\partial y^{+}=1/\left(\kappa y^{+}\right).$ In the following, this is called the modified wall model, and the one presented before is called the original model.

4. Result

The above IB method is implemented into the U²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.

4.1. 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 $R{e}_L=5\times {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 ${\mathrm{∆}}_{x_N}$ to the size of the start cell ${\mathrm{∆}}_{x_1}$ along the edge direction ($R={\mathrm{∆}}_{x_N}/{\mathrm{∆}}_{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.

Table 1. Modified wall model performances on a single-core computer.

Grid	CPU Time in Modified Wall Model (s)			Total CPU Time (s)	Percent (%)
	Definition of IB Cells	Image Stencil	Boundary Condition
$N_{y2}=100$	20.68	1.281	0.0004	31.963	68.7
$N_{y2}=200$	23.61	1.678	0.0003	41.943	60.3
$N_{y2}=300$	16.43	2.155	0.0004	52.361	54.6

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.

4.2. 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 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).

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.

4.3. 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.

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]:

$$\left|{\tau }_b\right|=u_{\tau }^2={\left|\frac{30\kappa U_1}{\ln \left(z_1-\eta \right)/k_s}\right|}^2,$$

(16)

where $U_1$ is the velocity at 5 mm above the bed, the distance $z_1-\eta$ is 5 mm, and $k_s=1$ mm is the bed roughness.

The wall shear stress $u_{\tau }^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.

5. 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²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.