1 Introduction

Large-eddy simulation (LES) has been increasingly applied to theoretical studies and engineering applications, particularly on flows with complex geometric boundaries. However, the LES suffers the prohibitive computational cost in resolving the turbulent boundary layer^[1]. The number of grid points required for the LES to resolve the viscous sublayer scales with Re_L^1.8[2], where Re_L is the Reynolds number based on the large-scale velocity and its integral length scale. More than 99% of the grid points have to be concentrated in the inner region of the turbulent boundary layer when Re_L>10^6[3]. The wall-modeled LES is able to significantly reduce the grid resolution, which bypasses the flow details near the wall to avoid fully resolving the inner region of the turbulent boundary layer. Different wall models have been proposed during the last decades^{[1, 3-4]}. The wall models were initially designed for the boundary-conformal mesh method, and recently extended to the non-boundary-conformal mesh method, e.g., the immersed boundary (IB) method^[5-6].

Tessicini et al.^[7] combined the IB method with wall models. They identified all grid points outside the solid bodies, and solved the flow equations down to the second grid point from the wall (see Fig. 1(a)). The flows between the second grid point and the wall were reconstructed based on an equilibrium stress balance model. The reconstruction can ensure the appropriate velocity profiles when the grid points are within the log layer, and thus extend the applicability of the IB method to turbulent flows. Cristallo and Verzicco^[8] improved the above procedure by solving the flow equations to the first grid point from the wall and computing the wall shear stress based on the equilibrium wall model. Choi et al.^[9] introduced the tangency correction to reconstruct the flows near the wall by using a power-law function without explicitly referring to the wall shear stress. Besides near-wall velocity reconstruction, Roman et al.^[10] introduced a near-wall eddy viscosity model with the consideration of near-wall damping effects. Yang et al.^[11] proposed a generalized off-wall boundary condition to combine wall models with the IB method, in which the conservation of the near-wall momentum was ensured by imposing velocity boundary conditions to convection terms and shear stress boundary conditions to diffusion terms. Yang et al.^[12] improved the wall-modeled LES with the IB method by referring to the integral wall model to account for the non-equilibrium effects. All these wall models are combined with the sharp-interface IB method, where the solid and fluid phases are treated distinctly^[5]. To implement the sharp-interface IB method, we need to identify the types of the points on the Eulerian grids, e.g., the fluid grid points, solid grid points, and interface grid points (see Fig. 1(a)). The type identification for the Eulerian grid points involves the famous problem of "point-in-polygon", which challenges the robustness and efficiency of the algorism in handling complex geometries and moving boundaries. In addition, the sharp-interface IB method suffers from spurious oscillations because of the spatial and temporal discontinuities of the flow near the IB, especially those associated with moving bodies^[5]. The issues of spatial and temporal discontinuities are not encountered in the diffuse-interface IB method because the diffuse-interface provides a smooth transition between the fluid phase and the solid phase^[5]. The diffuse-interface IB method does not need to identify the type of Eulerian grid points, which is efficient and robust, especially for the massively-parallel computing with domain decomposition strategy. However, the diffuse-interface IB method has rarely been applied to the LES of turbulent flows because of the difficulties in fully resolving the structures of near-wall flows and/or combining the wall models. It is still unclear whether the diffuse-interface IB method is applicable to the simulation of turbulent flows or not.

The advantages of the diffuse-interface method have been shown in the LES of wind turbines with rotor models. The rotor models (such as the actuator line model)^[13] for wind turbines share some similar ideas with the diffuse-interface IB method, where the effects of blades on the incoming wind are modeled by an effective body force, and are computed with the blade element method, parameterized airfoil geometry (e.g., chord length and twist angle), and aerodynamic (e.g., drag and lift coefficients) information instead of resolving the boundary layer flow around the blade. To model the nacelle of a wind turbine, Yang and Sotiropoulos^[14] recently developed an actuator surface model, and calculated the normal component of the force with the direct forcing IB method under the non-penetration condition. The tangential component of the force was computed by using the incoming wind speed and an empirical friction coefficient, because the grid was extremely coarse and could not even resolve the logarithmic layer. The diffuse-interface enables the implementation of the model without identifying the grid points outside the nacelles/at different sides of the blades (see Fig. 1(b)). The circumvention of identifying different grid points ensures the efficient coupling between the turbulent flows and the motion of wind turbines. Though the rotor models instead of the wall models are mainly used, the successful applications to different wind turbines^[15] indicate the potential to combine the wall models with the diffuse-interface IB method.

The aim of this work is to propose a novel method for combining the diffuse-interface IB method with wall models for the LES. We use wall models to reduce the computational cost while retain the advantages of the diffuse-interface IB method in handling complex geometries on Cartesian grids. The proposed method is validated by the benchmark simulation of the DARPA SUBOFF model, which is one of the well-documented turbulent benchmarks for complex geometries. The flows over the DARPA SUBOFF model involve both adverse and favorable pressure gradients, which are of great challenges for the wall-modeled LES. We apply the diffuse-interface IB method to the simulation of turbulent flows, and the results obtained are well comparable to the ones reported in the literature. The remainder of the paper is organized as follows. In Section 2, the flow equations and the numerical methods are described. In Section 3, the method for combining wall models with the diffuse-interface IB method and the implementations of different kinds of wall models are presented. In Section 4, the application of the proposed method to the simulation of flow around the DARPA SUBOFF model is presented. Finally, in Section 5, the summary and conclusions are given.

2 Numerical method

We solve the filtered Navier-Stokes (NS) equations for incompressible flows as follows:

(1)

(2)

where and are the filtered velocity components and pressure, respectively, f_i (i=1, 2, 3) are the body forces representing the boundary effects on the flows for the IB method, ν is the kinematic viscosity of the fluid, and is the sub-grid stress (SGS) given by

(3)

in which ν_t is the SGS viscosity, and is the strain rate tensor based on the resolved velocity field. The wall-adapting local eddy-viscosity (WALE) model^[16] is used for calculating ν_t, i.e.,

(4)

where C_w=0.6, and is the traceless symmetric part of the tensor defined by

(5)

In the above equation,

The flow equations are solved on an Eulerian grid without conforming to the immersed walls, for which geometry and kinematics are described by a set of Lagrangian points (see Fig. 1(b)). Equations (1) and (2) are solved by using a projection method as follows:

(6)

(7)

(8)

(9)

where the superscript n denotes the previous time step, are the intermediate velocity components, h_r contains the discretized advective and diffusive terms, and δp is the correction to the pressure. All the spatial derivatives are discretized by using a second-order center difference scheme on a staggered grid. The evaluation of the momentum forcing term f_i will be discussed in the next section. The details of the flow solver and the related validations can be found in Refs. [17] and [18].

3 Wall models for the diffuse-interface IB method 3.1 Difficulties in implementing wall models for the diffuse-interface IB method

The effective body forces, f_i^n+1/2 in Eq. (6), represent the effects of immersed walls on flows. To implement wall models in the framework of diffuse-interface IB methods, we first focus on the direct-forcing diffuse-interface IB method proposed by Vanella and Balaras^[19], which can be easily extended to other direct-forcing IB methods with diffuse-interfaces. According to Vanella and Balaras^[19], the effective body forces can be calculated by

(10)

where W_ij are interpolation operators based on the moving-least-squares with a linear basis^[19], and F_j^n+1/2 are the effective forces obtained by Lagrangian points, i.e.,

(11)

In the above equation, U_j^d are the desired velocity boundary conditions specified on the immersed walls, are the predicted velocities at the Lagrangian points, and

(12)

where is the predicted velocity on the Eulerian grid computed by using Eq. (6) without the effective body forces f_i^n+1/2, i.e.,

(13)

The effective body forces obtained from Eq. (10) have been shown to be able to give reasonable results for laminar flows as long as the grid is fine enough to resolve the near-wall flows. In the wall-modeled LES, the near-wall grids are not fine enough to resolve the viscous sublayer. The interpolation based on a linear basis (see Eq. (12)) is not valid anymore, which results in inaccurate predicted velocities and inaccurate effective body forces at the Lagrangian points. In the sharp-interface IB method, an explicit local velocity reconstruction based on wall models is used to account for the non-linear velocity distribution near the wall, for which the effective body force on an Eulerian grid can be directly obtained by

(14)

where u_i^r is the velocity near the wall reconstructed based on the wall model, and is the predicted velocity obtained from Eq. (13) at the same Eulerian point.

As discussed above, the sharp-interface IB method uses the velocity near the wall (see u_i^r and in Eq. (14)), which can be easily reconstructed using a wall model, to compute the effective body force on the Eulerian grid, while the effective body force in the diffuse-interface IB method is usually computed by using the velocity on the wall (see U_j^d and U_j^** in Eq. (11)). This makes the explicit local velocity reconstruction used in the wall model for the sharp-interface IB method cannot be directly used for the diffuse-interface IB method. This is one of the major challenges. The wall shear stress boundary condition is usually used in the wall-modeled LES as it can model the wall effects on the outflow with a coarse mesh. To apply wall models to the diffuse-interface IB method, a feasible approach is used to compute the effective body force based on the wall shear stress and/or introduce an explicit reconstruction near the diffuse-interface.

3.2 Combination of wall models with the diffuse-interface IB method

We propose a method for combining wall models with the diffuse-interface IB method, with which the effective body force is calculated based on the wall shear stress and a local flow reconstruction on a set of Lagrangian points near the wall. The combination is processed as follows.

(ⅰ) Decompose the velocities in a local orthogonal coordinate system

Wall models usually need the velocity in the tangential direction of the wall to compute the wall shear stress. A convenient way to implement wall models is to use the local orthogonal coordinate system with its coordinate origin fixed on the corresponding Lagrangian point on the wall (see Fig. 2), where the axis η points outside the solid body along the normal direction, the axis ξ is parallel to the tangential velocity near the wall, and the axis ζ goes along the direction perpendicular to the axes η and ξ according to the right-hand rule. The flow velocity at the near-wall Lagrangian point M (see Fig. 2) can be decomposed as follows:

(15)

where u_w is the velocity at Point S on the wall, e_η and e_ξ are the unit vectors of the local orthogonal coordinate system along the normal and tangential directions, respectively, and u_η and u_ξ are the normal and tangential components of the velocity in the local coordinate system, respectively. The velocity on the wall u_w is zero in the current work due to a stationary solid body. The velocity component u_ξ is zero (not shown in Eq. (15)) in the local orthogonal coordinate system because of the definition of the ξ-direction. Similarly, the effective body forces are decomposed in the local orthogonal coordinate system as follows:

(16)

where f_η and f_ξ are the normal and tangential components of the effective body forces on the Eulerian grid, respectively, and F_η and F_ξ are the normal and tangential components of the effective body forces on the Lagrangian grid, respectively. We calculate the normal and tangential components of the effective body forces separately.

(ⅱ) Compute the tangential components of the effective body forces

We integrate the momentum equation along the normal direction to link the effective body force to the wall shear stress, i.e.,

(17)

where Δ is the distance between Points M and S on the wall. We set Δ to be the same as the Eulerian grid length Δh, which is about 200 wall units in the current simulation. Equation (17) can be expressed in terms of the variables at the Lagragian points between Points S and M by using the mean value theorem for integrals, i.e.,

(18)

where Δ₁ and Δ₂ are points between Points S and M. The capital letters indicate the variables at the Lagrangian points. τ_w and τ_Δ are the shear stresses at the wall and Point M, respectively. For the wall-modeled LES of the turbulent flows over the stationary wall, the right-hand-side of Eq. (18) can be approximated by the dominant wall shear stress term as follows:

(19)

It is worth mentioning that Eq. (18) is consistent with the direct-forcing method proposed by Vanella and Balaras^[19] if the grid is fine enough to compute the spatial derivatives of the flows. The discretized form of Eq.(18) is

(20)

where H_R includes all the discretized advection term, the diffusion term, and the pressure gradient in the NS equations. Equation (20) reduces to Eq. (11) as Point M moves towards the wall (Δ, Δ₁, Δ₂ → 0), which recovers the direct-forcing method for laminar flows. If the grid resolution is far from resolving the flow details around the wall (or too coarse to resolve the log-law region of the boundary layer), the effective body forces on the Lagrangian points (see Eq. (18)) can be computed by an integral model for specific flows. The proposed model then reduces to some kinds of actuator type models, e.g., the actuator line/surface model for wind turbines. In this sense, Eq. (18) reduces to the results of Refs. [13] and [14] for the LES with rotor models.

(ⅲ) Compute the normal components of the effective body forces

The normal component of the effective body forces is computed at a Lagrangian point Q which has a distance Δ₁ (0 < Δ₁ < Δ) away from the wall, i.e.,

(21)

where U_η^**(Δ₁) is the intermediate normal velocity component near the wall in the predictor step interpolated from the intermediate velocities based on the moving-least-square reconstruction. U_η(Δ₁) is the desired normal velocity component near the wall, and can be computed from an interpolation using the velocity at Point P (see Fig. 2). In order to satisfy the impenetrability condition and the zero normal derivative of the wall-normal velocity component at the wall, a parabolic distribution of the wall-normal velocity between Points S and P is assumed^[10]. Therefore,

(22)

where Δ_p is the distance from Point P to the wall.

Note that both the tangential and normal components of the effective body forces are computed at the auxiliary point Q near the wall instead of Point S on the wall (see Fig. 2). The Lagrangian points near the wall, which are denoted as Q in Fig. 2, are the new Lagangrian points introduced in this work as a surrogate wall to combine wall models with the diffuse-interface IB method. We find that the setup of Δ₁ = Δh/4.0 gives acceptable results for all the simulation in the current work. We will show the effects of Δ₁ (Δ₁∈ [0, Δh/2.0]) in the next section. The effective body forces computed at the auxiliary Lagrangian point Q are applied to the Eulerian grids by using Eq. (10) in the same way as in the diffuse-interface IB method.

3.3 Implementation of the wall models

The wall shear stress can be computed based on the unsteady thin-boundary-layer equation (TBLE)^[20-21] as follows:

(23)

where

(24)

and ν and ν_t are the fluid kinematic viscosity and eddy viscosity, respectively. is the filtered near-wall pressure from the outer flow. is detected from a point along the normal direction of the wall (hereinafter referred to as the probe point) as Point P shown in Fig. 2. Note that the probe point P is not necessary to be coincident with Point M, since the probe point P is used by the wall model while Point M is an auxiliary point to introduce the forcing point Q for the IB method. ũ is the filtered velocity vector. e_ξ is the unit vector of the local orthogonal coordinate system along the tangential direction, and η indicates the direction normal to the wall.

We investigate two simplified models of Eq. (23). The first one has the source term in Eq. (23) as S=0, which corresponds to the boundary layer equation with the local equilibrium hypothesis (hereinafter referred to as the EB model). The second one has the source term with , which accounts for the non-equilibrium effect caused by the pressure gradient (hereinafter referred to as the NEB). For the EB and NEB models, Eq. (23) reduces to an ordinary differential equation, which can be integrated from the probe point down to the wall to give a closed-form expression for the wall shear stress, i.e.,

(25)

where ũ(Δ_p) is the filtered velocity at the probe point P.

Following the work of Duprat et al.^[22], we use the following definition of eddy viscosity which takes the pressure gradients into account:

(26)

where the mixing scaling for the normalized wall distance is defined as η^*=ηu_τp/ν, where is a combination of the frictional velocity and the velocity based on the pressure gradient . The parameter in Eq. (26) quantifies the preponderant effect between the shear stress and the pressure gradient. κ =0.41, β =0.78, and A=18 are constant according to the work of Duprat et al.^[22]. If the effect of the pressure gradient is negligible, Eq. (26) will recover to the van Driest formula^[23], which gives the velocity profile for the flow with zero pressure gradient.

To determine η^* and thus the wall-layer eddy viscosity ν_t in Eq. (26), the friction velocity u_τ is required, which depends on the wall shear stress. In the present implementation, is evaluated using the instantaneous τ_w from the previous time step. In this sense, the simplified models given by Eq. (25) are algebraic, and the wall shear stress can be reconstructed by the numerical integration of Eq. (25) under the boundary condition provided by the probe point in each time step.

We also implement the Werner-Wengle wall model^[24] (hereinafter refer to as the WW model), which is based on the assumption of a power law velocity profile outside the viscous sublayer, i.e.,

(27)

where u⁺=|u_ξ|/u_τ with u_ξ being the tangential velocity, η⁺=ηu_τ/ν, A₁=8.3, and B=1/7. The wall shear stress can be computed by

(28)

where u_ξ(Δ_p) is the tangential velocity at the probe point P, and . Δ_p is the distance between the probe point and the wall.

4 Wall-modeling for the LES of flows around the DARPA SUBOFF model 4.1 Numerical setup

We validate the proposed method by using the simulation of flows over the DARPA SUBOFF model. The axisymmetric bare hull of the SUBOFF model is used in the simulation. The hull consists of a streamlined forebody, a parallel middle body, and a stern with contraction in the radial direction (see Fig. 3). The hull has a maximum diameter D and a total length L, satisfying L/D = 8.6. The model is fixed in uniform flows at 0 angle of attack and 0 yaw angle. The Reynolds number based on the velocity of the uniform upstream flow and the length of the model is Re_L = U_∞L/ν = 1.2 × 10⁶.

Simulation is conducted in a domain of [-4.3D, 4.3D]×[-4.3D, 4.3D]×[-2.6D, 23.2D], according to the setup of Posa and and Balaras^[25]. The uniform upstream flow is specified at the inlet, and the free-convection boundary condition is set at the outlet. The proposed method along with the diffuse-interface IB method is used to represent the effect of hull on the flows with different wall models. The free-slip boundary condition is used at other boundaries. The flow is initially developed by the direct-forcing IB method without wall models, and then is restarted with wall models.

The computational domain is discretized by a block-structured Eulerian mesh with about 53 million grid points. The Eulerian grid length is Δh=0.033 6D, which corresponds to Δh⁺≈ 225 near the middle body. The geometry of the SUBOFF is given by a set of Lagrangian points distributed on the model surface. The averaged distance between the Lagrangian points is approximately equal to the Eulerian grid length. The time step is dynamically adjusted to keep the maximum Courant-Friedrichs-Lewy number at 0.3. The distance from the immersed surface to the probe point P varies from Δ_p=2.0 Δh to Δ _p=1.0 Δh. The semi-width of the interpolation stencil is Δ_r=1.2 Δh, and the probe point is carefully selected to keep the undesired pollution away from the diffuse-interface.

4.2 Distributions of the velocity and pressure

The instantaneous distribution of the streamwise velocity on a symmetric plane obtained with the NEB model is shown in Fig. 4(a). The velocity distributions obatined with the EB and WW models are similar. The low speed regions are concentrated around the stagnation point, boundary layer, and stern. There is no apparent flow separation at the streamlined forebody and parallel middle body, and the flows are almost axisymmetric (see Figs. 4(b) and 4(c)). The disturbances to the streamwise flow increase near the stern (see the distributions of the steamwise velocity in Figs. 4(d) and 4(e) and the vorticity magnitude in Fig. 5(b)). The eddies move downstream and form vortex structures in the wake (see Fig. 5).

The interaction between the upstream flows and the hull results in a deficit for the velocity profile in the wake. The time-averaged streamwise velocity profiles by the NEB wall model in the intermediate wake are shown in Fig. 6. The comparison is provided in the self-similar coordinates, which are the maximum velocity deficit u₀ and half-wake width l₀.

The velocity defects located at three different points in the wake are self-similar, since they nearly collapse onto one single curve at the scaled vertical distances. The mean velocity profiles of the current simulation are in good agreement with the similarity law and the results in Refs. [25] and [27].

The distribution of the time-averaged pressure in a symmetric plane is shown in Fig. 7. The pressure reaches the maximum at the front stagnation point, and decreases along the forebody. A nearly constant pressure is achieved at the parallel middle body. The adverse pressure gradient is generated at the front part of the stern.

The distribution of the time-averaged pressure coefficient on the meridian of the hull is shown in Fig. 8, where the time-averaged pressure coefficient is computed in terms of

(29)

where 〈p〉 is the time-average pressure. p_∞ and U_∞ are the free-stream pressure and velocity, respectively. We have calculated the pressure coefficient with different wall models. Figure 8 shows the pressure coefficient predicted by using the WW, EB, and NEB models, respectively, where the numerical results obtained by Posa and and Balaras^[25] are wall-resolved LES with full appendages while the results obtained by Posa and and Balaras^[25] and Jimenez et al.^[27] are obtained from the side opposite to the sail. The distributions of the pressure coefficient are not sensitive to different wall models for most locations except in the trailing edge, where the pressure coefficient obtained by the NEB model is higher and agrees better with the data of Jemenez et al.^[27] and Huang et al.^[26]. Overall, all the wall models give an acceptable prediction for the distribution of the pressure coefficient.

4.3 Distribution of the skin-friction coefficient

The time-averaged skin-friction coefficient is defined as follows:

(30)

where 〈τ_w〉 is the wall shear stress.

Figure 9 shows the distributions of the time-averaged skin-friction coefficients predicted by different wall models, where the probe point has a distance of Δ_p=2.0 Δh from the wall.

From Fig. 9, we can see that both the WW model and the EB model give the nearly constant skin-friction coefficient on the parallel middle body. However, the peaks of the skin-friction on the stern predicted by the WW model and the EB model are shifted toward the downstream, and neither the WW model nor the EB model obtains the peak of the skin-friction coefficient on the forebody. The NEB model obtains both of the peaks near the forebody and the stern, though the skin-friction coefficient is over-predicted.

The over-prediction of the skin-friction coefficient can be circumvented by moving the probe point toward the wall. The distance between the probe point and the wall is set to be far away enough from the wall to reduce the effect of the diffuse-interface on the wall model and to be close enough to the wall to ensure that the wall models are able to capture the flow features. We use Δ _p=2.0 Δh in the above simulation, and have the results of Δ_p⁺ ranging from 300 to 400 near the forebody and middle body of the hull. Grid refinements on the Eulerian meshes can move the probe point towards the wall. However, the grid refinements near the wall usually cause a dramatic increase in the computational cost. Instead of refining the Eulerian grid, we move the probe point P towards the wall by reducing Δ_p from 2.0 Δh to 1.5 Δh, 1.2 Δh, and 1.0 Δh, respectively. The prediction of the skin-friction coefficient is improved when the probe point P moves toward the wall (see Fig. 10). The skin-friction coefficient predicted by the non-equilibrium model is comparable to the wall-resolved LES of Posa and Balaras^[25] and the rescaled experimental measurement of Huang et al.^[26] when the probe point P has a distance of Δ_p < 1.2 Δh (corresponding to a Δ_p⁺ ranging from 150 to 240 near the forebody and middle body of the hull).

The mean value theorem for integrals does not tell us the exact distance between the surrogate wall and the real wall, i.e., the exact value of Δ₁ in Eqs. (17) and (18). The exact value of Δ₁ depends on the distribution of the effective body force (or the velocity) near the wall. When the effective body force is linearly distributed within [0, Δh], we have Δ₁ = Δh/ 2.0. When the effective body force has a nonlinear distribution and is concentrated near the wall, 0 < Δ₁ < Δh/2.0. We use the setup of Δ₁ = Δh/4.0 in all the above simulation. We have investigated the effects of Δ₁ on the distributions of the pressure coefficient and the skin-friction coefficient in the simulation with the NEB model and Δ_p = 1.2 Δh. Neither the pressure coefficient nor the skin-friction coefficient sensitively depends on the distance between the surrogate wall and the real wall when 0 < Δ₁ < Δh/2.0, except some deviation around the stern in the case with Δ₁ = Δh/2.0, as shown in Fig. 11. The deviations of the pressure coefficient and skin-friction coefficient associated with the case of Δ₁ = Δh/2.0 should be caused by the interaction between the diffuse-interface and the wall model. The distance between the surrogate wall and the real wall (Δ₁) affects the normal velocity at the surrogate wall (see Eq. (22)) and the region on which the effective body forces are spread (see Eq. (10)). For a given probe point at Δ_p = 1.2 Δh and a fixed semi-width of the interpolation stencil Δ_r=1.2 Δh, a large Δ₁ will result in a large overlap domain in spreading the effective body forces (see Eq. (10)) and in computing the velocity at the probe point (see Eq. (12)), and thus increase the interaction between the wall models and the diffuse-interface. We suggest to set up Δ₁ to reduce the interaction between the diffuse-interface and the probe points, since the velocities given by the wall models are determined by the flow at the probe points.

5 Summary and conclusions

The prohibitive computational cost for resolving the inner region of the turbulent boundary layer blocks the applications of the LES to high-Reynolds number turbulent flows with complex boundaries. The wall-modeled LES with the IB method is expected to reduce the computational cost and handle the complex boundaries. Several approaches have been proposed for implementing wall models in the sharp-interface IB method, in which the velocities near the wall are reconstructed based on wall models. However, the implementation of the wall models in the diffuse-interface IB method has been little explored, though the diffuse-interface IB method has advantages in the efficiency and robustness in handling moving boundaries. The major challenge is that the diffuse-interface IB method does not have an explicit procedure for reconstructing the velocity near the wall.

We propose a method for implementing wall models in the diffuse-interface direct-forcing IB method. In the method, the effective body forces representing the effect of walls on the flows are computed based on the wall shear stress and a local velocity reconstruction at a set of Lagrangian points near the wall. We integrate the momentum equation along the wall-normal direction to link the tangential component of the effective body force for the IB method to the wall shear stress predicted by the wall models. We introduce a set of Lagrangian points near the wall to compute the normal component of the effective body force for the IB method by reconstructing the normal component of the velocity near the wall. The proposed method reduces to a classical direct-forcing IB method if the grid is fine enough to resolve all the spatial gradients of the flows near the wall. If the grid is far from resolving the flows near the wall, the proposed method can be reduced to the actuator type models, e.g., actuator line/surface models for wind turbines, by introducing integral models for the IB.

The proposed method has been validated by the benchmark simulation of flows around the DARPA SUBOFF by using the wall-modeled LES with three different wall models, i.e., the equilibrium stress balance model, the non-equilibrium stress balance model with the pressure gradient, and the WW model. The distribution of the streamwise velocity in the wake, the pressure coefficient, and the skin-friction coefficients on the wall predicted by different models are reported and compared with the wall-resolved LES and experimental results in the literature. An acceptable agreement is obtained for all the three models, while the non-equilibrium stress balance model with the pressure gradient term gives better predictions on the peaks of pressure and skin-friction coefficients.

We show that the diffuse-interface IB method combined with wall models has the capability in handling wall-modeled LES of turbulent flows. Though the proposed method is validated by a benchmark flow with both adverse and favorable pressure gradients, the capability of the proposed method in predicting separations and reattachments is still to be investigated. We combine a direct-forcing diffuse-interface IB method with different algebraic wall models, where the non-equilibrium effect associated with the pressure gradient is taken into account. However, the combination of the direct-forcing diffuse-interface IB method with the ordinary differential equation/partial differential equation (ODE/PDE) based wall models is expected for the LES of turbulent flows with strong non-equilibrium effects.

Acknowledgements

The author Xiaolei YANG would like to acknowledge the hospitality received at LNM during his visit where he accomplished this work. The computations are conducted on Tianhe-1 at the National Supercomputer Center in Tianjin.

Open Access

This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article's Creative Com- mons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.