Adaptation of the hybrid fictitious domain-immersed boundary method for Reynolds-averaged turbulence modeling
June 04, 2026
Abstract
Engineering practice often calls for shape or topology optimization (TO) of fluid defining components, while the ever-increasing computing power allows the optimized cost functions to be based on computational fluid dynamics (CFD). However, a common bottleneck in CFD-based TO frameworks is the requirement for frequent remeshing. In order to alleviate this bottleneck, we propose an adaptation of an immersed boundary (IB) method variant, the hybrid fictitious domain-immersed boundary method, to leverage Reynolds-averaged Navier-Stokes (RANS) equations and wall function. The main contribution of the present work lies in the design and open-source implementation of the IB-aware steady-state solution of the RANS equations via the SIMPLE algorithm in the OpenFOAM library. For the most common two-equation RANS models, Reynolds numbers from \(10^1\) to \(10^6\), and several benchmarks, such as flow over a backwards facing step or an Ahmed body, the framework gives results consistent with the standard body-fitted CFD. Furthermore, given the intended application in TO, special emphasis is placed on the robustness and applicability of the approach to general geometries, which is tested on a NACA profile under various angles of attack.
computational fluid dynamics, immersed boundary method, Reynolds averaged simulation, wall functions
1 Introduction↩︎
An immersed boundary method (IBM) [1] is a numerical tool used in computational fluid dynamics (CFD) that enables flows to be simulated around complex geometries without the need to create complex geometry-conforming meshes. Instead of a complex mesh, a simple one is used, with the geometry represented by a scalar indicator field and, in a direct forcing IBM, by additional force terms [2]. If the geometry changes, regenerating the scalar field is much cheaper than re-meshing. Thus, utilizing IBMs is very advantageous in applications that require geometry optimization. However, geometry optimization of a real-life component needs the IBM to be connected to feasible turbulence models, e.g., models based on Reynolds-averaged Navier-Stokes (RANS) equations.
The connection of IBMs to Reynolds-averaged turbulence models is a challenging task where the most limiting factor is the modeling of the boundary layer [3]. While having a non-conforming mesh, IBMs cannot leverage anisotropic mesh refinement near geometry walls, and therefore often rely on wall modeling approaches. In [4], the velocity at the immersed boundary is reconstructed by linear interpolation and its tangential component is corrected by wall functions. A different two-layer approach was proposed by [5], where a thin-boundary layer equation is solved near the IB and RANS equations in the rest of the domain.
The inclusion of wall modeling approaches showed promising results, but in flows with high Reynolds numbers, spurious oscillations of pressure and wall shear stress occurred at the immersed boundary [6]. [7] suggested avoiding these by interpolating the friction velocity rather than the raw velocity at the IB. Alternatively, [8] proposed a pre-processing step to optimize the positions of interpolation points based on an estimate of the boundary layer thickness. Moreover, [9] discussed another reason for instability arising from a not tight enough connection between the momentum and turbulence variables. When the local velocity profile is used to calculate the friction velocity, which in turn is used to correct the velocity at the immersed boundary, limited feedback is gained from the turbulence transport equations.
Another issue arises when solving the incompressible form of the RANS equations in direct forcing IBMs. As described in a review by [3], the forcing terms first impose correct velocity at the immersed boundary, but second, a solenoidal field is enforced by the incompressibility condition, and the velocity is changed indifferently of the IB. Several remedies for the problem were suggested, including additional iteration loops [10], Lagrange multipliers [11], or mass sources and sinks [12].
Overall, there is no standard approach to the connection of IBM and RANS, and of the many there are, a limited number of them is open-source. In this work, we aim to provide our own version of an IBM-RANS solver, but the source code is fully open and built upon the OpenFOAM library, incorporating as much of the standard CFD practices as possible. In particular, the IBM variant used is the hybrid fictitious domain-immersed boundary method (HFDIB), originally developed by [13] and further extended to include the discrete element method (DEM) for arbitrarily shaped particles (openHFDIB-DEM) by [14] and [15]. Here presented is the connection of HFDIB to RANS models (openHFDIBRANS), which implements the steady-state incompressible RANS equations with four two-equation turbulence models: the \(k\)-\(\omega\), \(k\)-\(\varepsilon\), \(k\)-\(\omega\) SST and realizable \(k\)-\(\varepsilon\) model.
To address the problems described above, our openHFDIBRANS solver utilizes a mixture of the approaches mentioned, IB-adjusted OpenFOAM practices, and custom additions. For wall modeling at the immersed boundary, wall functions are used, and to avoid instability of the solution, the friction velocity is evaluated using turbulence variables rather than local velocity profiles. Moreover, in high Reynolds number flows, the velocity is not reconstructed by interpolation but is corrected solely by the turbulence viscosity, gaining a tight connection between the momentum and turbulence variables. Last but not least, the incompressibility condition is only enforced in the fluid region of the computational domain to avoid conflict with the IBM forcing. This may result in a loss or gain of mass in the first iterations, but in steady-state, it leads to convergence. The solver performance is showcased and tested in a series of benchmarks, including the backward-facing step benchmark, flow over a smooth cylinder, or the Ahmed body.
2 Solver description↩︎
Here, the proposed variant of the immersed boundary (IB) method and its coupling with the Reynolds-averaged turbulence models are presented. The method stems from the hybrid fictitious domain-immersed boundary (HFDIB) method, an IB variant initially introduced by [13] and further refined by [14] for applications in particle-resolved direct numerical simulation of particle-laden flows; see [15], [16]. Following the naming logic from the aforementioned publications, we denote the here-presented method as openHFDIBRANS. The openHFDIBRANS solver is implemented using the open-source C++ library OpenFOAM [17]. Its complete source codes and examples are available from https://github.com/techMathGroup/openHFDIBRANS.
In the following, we first introduce the general governing equations solved. Afterwards, details on the IB treatment are provided. Particular emphasis is placed on distinguishing between resolved and unresolved boundary layer scenarios and evaluating wall functions at the immersed boundary.
2.1 Governing equations↩︎
The governing equations utilized in openHFDIBRANS are the steady-state Reynolds-averaged Navier-Stokes (RANS) equations complemented by two-equation turbulence closure models. The implemented form of the RANS equations is \[\label{eq:hfdibrans} \begin{array}{*1{>{\displaystyle}c}} \nabla \cdot (\boldsymbol{u} \otimes \boldsymbol{u}) = \nabla\cdot\left\{ \left(\nu + \nu_\mathrm{t}\right) \left[ \nabla\boldsymbol{u} + \nabla \boldsymbol{u}^{\text{T}} \right] \right\} - \nabla \tilde{p} + \boldsymbol{f}_\mathrm{ib}\,, \\ [0.3cm] (1 - \lambda)\left(\nabla \cdot \boldsymbol{u}\right) = 0\,, \end{array}\tag{1}\] where \(\boldsymbol{u}\) is the average velocity, \(\nu\) the kinematic viscosity, \(\nu_\mathrm{t}\) is the turbulent viscosity, and \(\tilde{p}\) the average kinematic pressure. The effect of the immersed boundary is included in the term \(\boldsymbol{f}_{\mathrm{ib}}\) and multiplicative factor \((1 - \lambda)\). In particular, the immersed boundary-induced source term \(\boldsymbol{f}_\mathrm{ib}\) simulates the effect of the solid body onto the flow momentum. Second, in the multiplicative factor, \(\lambda\) is a phase indicator field and by multiplying the conservation equation by \((1 - \lambda)\), its solution is limited to the fluid part of the domain.
To complement the RANS equations, four turbulence closure models are available in openHFDIBRANS. Namely, the implemented models are the \(k\)-\(\varepsilon\) model by [18] and [19], the \(k\)-\(\omega\) model by [20], the Menter’s \(k\)-\(\omega\) shear stress transport (SST) model [21] and the realizable \(k\)-\(\varepsilon\) model by [22].
In openHFDIBRANS, all turbulence models are implemented in a similar way. Therefore, the \(k\)-\(\omega\) model was chosen to illustrate the implementation details. The equations of the \(k\)-\(\omega\) model [20] are included in the following form \[\label{eq:hfdibkomega} \begin{array}{*1{>{\displaystyle}c}} \nabla \cdot (\boldsymbol{u}\, k) = \nabla \cdot \left[ \left(\nu + \frac{\nu_\mathrm{t}}{\sigma_{k_1}} \right) \nabla k \right] + P_k - \beta^*\, k\, \omega + S_\mathrm{ib}\,,\\[0.5cm] \nabla \cdot (\boldsymbol{u}\, \omega) = \nabla \cdot \left[ \left( \nu + \frac{\nu_\mathrm{t}}{\sigma_{\omega_1}} \right) \nabla \omega \right] + C_{\alpha 1}\, \frac{\omega}{k}\, P_k - C_{\beta 1}\, \omega^2 \,,\\[0.5cm] \end{array}\tag{2}\] where \(k\) is the turbulence kinetic energy, \(\omega\) the specific dissipation of \(k\), \(P_k\) the production of \(k\), and \(S_\mathrm{ib}\) is another immersed boundary-induced source term. The rest of the symbols represent model constants [20].
The steady-state governing equations 1 and 2 are solved on \(\Omega \subset \mathbb{R}^3\), an open and connected set representing the spatial computational domain. Note that since immersed boundary approach is of interest here, \(\Omega\) is assumed to comprise both the solid phase \(\Omega_\mathrm{s}\) and the fluid phase \(\Omega_\mathrm{f}\).
2.2 Immersed boundary-induced source terms↩︎
In the governing equations, two immersed boundary-induced source terms are present: \(\boldsymbol{f}_\mathrm{ib}\) in the momentum equation 1 \(_1\) and \(S_\mathrm{ib}\) in the conservation equation of \(k\) 2 \(_1\). These are computed as \[\label{eq:qibalpha} \boldsymbol{f}_{\mathrm{ib}} = \alpha_{\boldsymbol{u}}\,\left[\,\mathcal{M}(\boldsymbol{u}_{\mathrm{ib}}) + \nabla \tilde{p}\,\right],\quad S_\mathrm{ib} = \alpha_k\;\mathcal{N}(k_\mathrm{ib})\,,\tag{3}\] where \(\mathcal{M}\) and \(\mathcal{N}\) are operators that stem from the governing equations as \[\label{eq:operators} \begin{array}{*1{>{\displaystyle}c}} \mathcal{M} (\boldsymbol{u}) = \nabla \cdot (\boldsymbol{u} \otimes \boldsymbol{u}) - \nabla \cdot\left\{ (\nu + \nu_\mathrm{t}) \left[ \nabla\boldsymbol{u} + \nabla \boldsymbol{u}^{\text{T}} \right] \right\}\,, \\ [0.3cm] \mathcal{N}(k) = \nabla \cdot (\boldsymbol{u}\, k) - \nabla \cdot \left[ \left(\nu + \frac{\nu_\mathrm{t}}{\sigma_{k_1}} \right) \nabla k \right] - P_k+ \beta^*\, k\, \omega\,. \end{array}\tag{4}\] Moreover, the IB-induced source terms are evaluated based on masking fields \(\alpha_{\boldsymbol{u}}\) and \(\alpha_k\) that limit the activity of the source terms to relevant parts of \(\Omega\), and imposed fields \(\boldsymbol{u}_\mathrm{ib}\) and \(k_\mathrm{ib}\). In cells where the IB-induced source terms are active, the values of the imposed fields are enforced during the iterative solution process [13].
2.3 Masking fields construction↩︎
To construct the masking fields \(\alpha_{\boldsymbol{u}}\) and \(\alpha_k\), the solid bodies need to be projected on the spatial domain \(\Omega\). For the description, we will focus solely on discrete setting. Thus, let \(\Omega^h\approx\Omega\) be a discrete spatial domain, see Figure 1a. Based on the geometry of the solid body, \(\Omega^h\) can be divided into three subdomains \[\label{eq:lamomega} \Omega^h = \Omega^h_\mathrm{f} \cup \Omega^h_\mathrm{sf} \cup \Omega^h_\mathrm{s}\,,\tag{5}\] where \(\Omega^h_\mathrm{f}\) contains cells that are fully immersed in fluid, \(\Omega^h_\mathrm{s}\) contains those fully immersed in solid, and \(\Omega^h_\mathrm{sf}\) those intersected by the fluid-solid interface, see Figure 1b. The domain division is represented by an indicator scalar field \(\lambda\), see Figure 1c. With \(\Omega^h_P\in\Omega^h\) being a computational cell, the \(\lambda\) field is defined as \[\label{eq:lambda} \lambda = \left\{ \begin{array}{rc} 0 & \text{if}\; \Omega^h_P \in \Omega^h_\mathrm{f} \\ [0.2cm] \tilde{\lambda} \in (0,1) & \text{if}\; \Omega^h_P \in \Omega^h_\mathrm{sf} \\ [0.2cm] 1 & \text{if}\; \Omega^h_P \in \Omega^h_\mathrm{s} \end{array} \right.,\quad \tilde{\lambda} = \frac{1}{2}\left[1 - \text{tanh}\left( \frac{\sigma}{V(\Omega^h_P) ^{\frac{1}{3}}} \right) \right]\,,\tag{6}\] where \(V(\Omega^h_P)\) is the volume of the cell and \(\sigma\) is the signed perpendicular distance from the cell center to the solid surface. The distance \(\sigma\) is positive when the cell center is in the fluid and negative when in the solid. When the cell center is exactly at the fluid-solid interface, \(\sigma\) is zero.
In openHFDIBRANS, the \(\lambda\) field is used to sort cells into three groups. The groups are visualized in Figure 1d. Each cell can be
an in-solid cell (\(\Omega^h_\mathrm{ic}\)) if \(\lambda \geq 0.5\),
a boundary cell (\(\Omega^h_\mathrm{bc}\)) if \(\lambda \in (0,0.5)\), or \(\lambda = 0\) and there is an in-solid cell as the cell neighbor, or
a free-stream cell (\(\Omega^h_\mathrm{fc}\)) when \(\lambda = 0.0\) and none of the cell neighbors are in-solid cells.
Based on the cell groups, the masking fields that limit the area of effect of the sources terms 3 are defined as \[\label{eq:alphauk} \alpha_{\boldsymbol{u}} = \alpha_k = \left\{ \begin{array}{rl} 1 & \text{if}\;\Omega^h_P \in \Omega^h_\mathrm{ic} \\ [0.2cm] 1 & \text{if}\;\Omega^h_P \in \Omega^h_\mathrm{bc}\;\text{and}\;y^+_\mathrm{ib} \leq y^+_\mathrm{lam} \\ [0.2cm] 0 & \text{elsewhere} \end{array} \right.\,,\tag{7}\] where \(y^+_\mathrm{ib}\) is the \(y^+\) value evaluated with respect to the immersed boundary, and \(y^+_\mathrm{lam} \approx 11.0\) is a constant commonly used as an approximate border between the viscous sublayer and the logarithmic region of the fluid boundary layer [23].
2.3.0.1 Treatment of production and dissipation of \(k\)
Note that specific IB-induced source terms were added only to equations for transport of linear momentum 1 \(_1\) and of turbulence kinetic energy 2 \(_1\). On the other hand, the immersed boundary also affects the values of production and dissipation of \(k\), i.e., of \(P_k\) and \(\omega\), respectively. But, for \(\omega\), this cannot be accounted for by introducing another IB-induced source term, since \(\omega\) goes to infinity when approaching smooth walls [20]. Instead, imposed fields \(\omega_\mathrm{ib}\) and \(P_\mathrm{ib}\) are calculated and assigned directly to \(\omega\) and \(P_k\) in cells with active masking field \(\alpha_\omega\). This masking field is defined as \[\label{eq:alphaomega} \alpha_\omega = \left\{ \begin{array}{rl} 1 & \text{if}\;\Omega^h_P \in (\Omega^h_\mathrm{ic} \cup \Omega^h_\mathrm{bc}) \\ [0.2cm] 0 & \text{elsewhere} \end{array} \right.\,.\tag{8}\] For the production of \(k\), the direct assignment of \(P_\mathrm{ib}\) is simple, since no transport equation for \(P_k\) is considered and its value is calculated explicitly before the construction of both turbulence model equations 2 . On the other hand, the assignment of \(\omega_\mathrm{ib}\) to the \(\omega\) field is carried out by direct manipulation with the discretization matrix of the \(\omega\) equation. In the affected cells, the \(\omega\) equation is then not solved. Note that the same approach is commonly adopted by OpenFOAM to enforce Dirichlet boundary conditions near walls for \(\omega\) [17].
2.4 Calculation of \(y^+\) at the immersed boundary↩︎
To switch between the variants of masking fields 7 , it is necessary to compute \(y^+_\mathrm{ib}\) at the immersed boundary, i.e., to decide if the boundary layer around IB is resolved. For a boundary cell \(\Omega^h_P \in\Omega_{\mathrm{bc}}^h\), \(y^+_\mathrm{ib}\) is computed as \[\label{eq:yplus} y^+_\mathrm{ib} = \frac{y_\mathrm{eff}\;u_\tau}{\nu},\quad y_\mathrm{eff} = 0.5\left[y_\perp + 0.5\,V(\Omega_P^h)^{1/3}\right],\tag{9}\] where \(y_\perp\) is the perpendicular distance from the cell center to the solid surface and \(V(\Omega_P^h)\) is the cell volume.
For a boundary cell \(\Omega^h_P\in\Omega_{\mathrm{bc}}^h\) and its in-solid neighbor \(\Omega^h_N\in\Omega_{\mathrm{ic}}^h\) that shares at least one vertex with \(\Omega^h_P\), the perpendicular distance \(y_\perp\) is computed as \[\label{eq:yortho} y_\perp = \left\{ \begin{array}{cl} |\sigma_P| & \text{if}\;\lambda_P \in (0,0.5) \\ [0.3cm] \boldsymbol{n}_N \cdot \left[P - \left(N + |\sigma_N|\boldsymbol{n}_N\right)\right] & \text{if}\;\lambda_P = 0\;\text{and}\;\lambda_N \in \left[0.5, 1.0\right) \\ [0.3cm] \|P - \langle v_{PN} \rangle\| & \text{if}\;\lambda_P = 0\;\text{and}\;\lambda_N = 1.0 \end{array} \right.,\tag{10}\] where \(\langle v_{PN} \rangle\) is the average of the shared vertices and \(\sigma\) and \(\boldsymbol{n}\) are computed as \[\label{eq:yorthoaux} \sigma = V(\Omega^h)^{\frac{1}{3}}\,\text{tanh}^{-1}\left(1 - 2\,\lambda \right),\quad \boldsymbol{n} = \frac{-(\nabla \lambda)}{\|(\nabla \lambda)\|}\,.\tag{11}\] To find the in-solid neighbor \(\Omega_N^h\), all neighbors of the boundary cell \(\Omega_P^h\) are checked and the one with the lowest value of \(-\boldsymbol{n}_P \cdot (P-N)\) is chosen, where \(P\) and \(N\) are the cell centers and \(\boldsymbol{n}_P\) is the outer unit normal to the solid surface constructed in \(P\).
Nevertheless, the perpendicular distance \(y_\perp\) is not used directly for the computation of \(y^+_\mathrm{ib}\) 9 . Instead, an effective perpendicular distance \(y_\mathrm{eff}\) is utilized, because it better represents that the boundary cell (or neighboring in-solid cell) should act as if cut by the solid surface. This is illustrated in Figure 3 for three different cases of solid surface position. Moreover, \(y_\mathrm{eff}\) helps smooth out \(y^+_\mathrm{ib}\). For instance, if cell pairs (a) and (c) from Figure 3 exist next to each other, the difference in \(y_\mathrm{eff}\) is much smaller than in \(y_\perp\) and the inclusion of \(y_\mathrm{eff}\) leads to a smoother transition in \(y^+_\mathrm{ib}\).
2.5 Calculation of imposed fields↩︎
The imposed fields \(\boldsymbol{u}_\mathrm{ib}\), \(k_\mathrm{ib}\), \(\omega_\mathrm{ib}\), and \(P_\mathrm{ib}\) are calculated differently for different cell groups, visualized in Figure 1d. In the in-solid cells \(\Omega_{\mathrm{ic}}^h\), the imposed fields are prescribed as \[\label{eq:imval} \boldsymbol{u}_\mathrm{ib} = \boldsymbol{u}_\mathrm{s},\quad k_\mathrm{ib} = 0,\quad \omega_\mathrm{ib} = \displaystyle\max_{\Omega^h_\mathrm{bc}\, \cup\;\Omega^h_\mathrm{fc}}(\omega^\mathrm{old}),\quad P_\mathrm{ib} = 0\,,\tag{12}\] where \(\boldsymbol{u}_\mathrm{s}\) is the velocity of the solid and \(\omega^\mathrm{old}\) are the values of \(\omega\) from the previous iteration or the initial guess. In the boundary cells \(\Omega_{\mathrm{bc}}^h\), the calculation of the imposed fields is based on the value of \(y^+_\mathrm{ib}\). The values of \(\boldsymbol{u}_\mathrm{ib}\) and \(k_\mathrm{ib}\) are required only when \(y^+_\mathrm{ib} \leq y^+_\mathrm{lam}\), since only then the corresponding masking fields are active in the boundary cells, see 7 . The calculation itself is then done using polynomial interpolation.
![Figure 4: Construction of the interpolation polynomial for a boundary cell \Omega^h_P. (a) Interpolation points located along the direction of the solid surface normal \boldsymbol{n}_P. The points P_1 and P_2 are in the free stream, S is on the solid surface, and P in the boundary cell center. (b) Example of a quadratic interpolation profile. Adapted from [24].](/article/2606.06135/main-figure3.png "fig:")
The interpolation polynomial is constructed in the same way as in openHFDIB [14], [15]. The construction process is illustrated in Figure 4. In particular, the interpolation is carried out in a normal direction to the solid surface, see Figure 4a, and the polynomial, see the example in Figure 4b, is constructed based on the values from the interpolation points located in the free stream and the value at the solid surface. More details on the interpolation process can be found in [14].
For calculations in openHFDIBRANS, the values in the free stream are taken from the previous iteration or the initial guess. Regarding the surface values, the velocity is set equal to the velocity of the solid \(\boldsymbol{u}_\mathrm{s}\) and \(k\) is calculated using wall functions [25]. The utilized form of the wall functions is the one used in OpenFOAM [17], i.e., \[\label{eq:kws} k = k_\mathrm{vis}\,u_\tau^2,\quad k_\mathrm{vis} = \frac{2400\,C_f}{C^2_\mathrm{eps2}},\quad C_f = \frac{1}{(y^+_\mathrm{ib} + C)^2} + \frac{2\,y^+_\mathrm{ib}}{C^3} - \frac{1}{C^2}\,,\tag{13}\] where \(C_\mathrm{eps2}\) and \(C\) are model constants [20] and \(u_\tau\) is the friction velocity calculated as \[\label{eq:frictionvel} u_\tau = C_{\nu}^{\,0.25}\, \sqrt{k_P}\tag{14}\] where \(C_\nu\) is another model constant, and \(k_P\) is the value of \(k\) in cell \(\Omega^h_P\) from the previous iteration or the initial guess.
The values of the imposed fields \(\omega_\mathrm{ib}\) and \(P_\mathrm{ib}\) are needed regardless of the \(y^+_\mathrm{ib}\) value. The corresponding masking field is always active, see 8 . The calculation of the imposed fields is done using wall functions [25] in the form used in OpenFOAM [17] which is \[\label{eq:omegaws} \begin{array}{rcl} \omega_\mathrm{ib} & = & \displaystyle\frac{u_\tau^2}{\nu}\cdot\left\{ \begin{array}{rc} \displaystyle\frac{6}{\beta_1\,(y^+_\mathrm{ib})^2} & \text{if}\;\;y^+_\mathrm{ib} \leq y^+_\mathrm{lam}\\ [0.5cm] \displaystyle\frac{1}{\kappa\,\sqrt{C_\nu}\,y^+_\mathrm{ib}} & \text{if}\;\;y^+_\mathrm{ib} > y^+_\mathrm{lam} \end{array} \right.\!\!\!\!,\\[1.0cm] P_\mathrm{ib} & = & \left\{ \begin{array}{rc} P_k & \text{if}\;\;y^+_\mathrm{ib} \leq y^+_\mathrm{lam}\\ [0.5cm] \displaystyle\frac{\kappa\,\nu\,y^+_\mathrm{ib}\,\|\boldsymbol{n}_\mathrm{ib}\cdot\nabla\boldsymbol{u}\|^2}{\log^2(E\,y^+_\mathrm{ib})} & \text{if}\;\;y^+_\mathrm{ib} > y^+_\mathrm{lam} \end{array} \right.\!\!\!\!, \end{array}\tag{15}\] where \(\beta_1\), \(\kappa\) and \(E\) are model constants [20].
2.5.1 SIMPLE-based solution strategy with immersed boundary enforcement↩︎
In openHFDIBRANS, the RANS equations 1 and the turbulence model equations 2 are solved using the SIMPLE algorithm [26] that is extended by additional steps. The extended solver loop is sketched in Figure 5.
First, before the momentum predictor [27], the values of \(\boldsymbol{u}_\mathrm{ib}\) are calculated. Then the momentum equation is solved and a correction on the resulting velocity field is applied as \[\label{eq:corru} \boldsymbol{u} = \boldsymbol{u} + \alpha_{\boldsymbol{u}}(\boldsymbol{u}_\mathrm{ib} - \boldsymbol{u})\,.\tag{16}\] By this correction, the values of \(\boldsymbol{u}_\mathrm{ib}\) are directly assigned to cells where the masking field \(\alpha_{\boldsymbol{u}}\) is active. This correction is important to ensure convergence in steady-state. The reason for this is that, for steady-state solution, the equation matrix has to be underrelaxed to ensure diagonal dominance [17]. The underrelaxation adds auxiliary sources evaluated based on the values of \(\boldsymbol{u}\) from the previous iteration. The correction 16 ensures that theses sources are as correct as possible. With the corrected velocity field, the pressure corrector and velocity update steps follow as usual [27].
After the velocity is updated, the turbulence model part of the solver is entered. There, first the values of the imposed fields \(\omega_\mathrm{ib}\) and \(P_\mathrm{ib}\) are calculated. The imposed field \(P_\mathrm{ib}\) can already be assigned to \(P_k\). Then, the \(\omega\) equation is discretized and the values of \(\omega_\mathrm{ib}\) are enforced by direct manipulation with the discretization matrix. Afterward, the \(\omega\) equation is solved and the values of \(k_\mathrm{ib}\) are computed. Here, similarly as for velocity, a correction on the \(k\) field is applied and the solver loop ends.
3 Results and Discussion↩︎
Several verification and validation tests were prepared to compare openHFDIBRANS performance with standard geometry-conforming simulations run with the simpleFoam (sF) solver available from OpenFOAM [17]. Six test case scenarios are presented in this section, namely
2D pipe flow,
3D pipe flow,
2D flow over a backward facing step,
2D flow over a cylindrical obstacle,
2D flow over a NACA-0009 profile, and
3D flow over an Ahmed body.
In each scenario, different aspects of openHFDIBRANS are showcased and tested. These range from sensitivity to Reynolds number, over interface position and angle, to parallelization. Moreover, a study of the openHFDIBRANS convergence, or, e.g., evolution of the system matrix stiffness is included. Lastly, all 2D test scenarios are available as tutorials in https://github.com/techMathGroup/openHFDIBRANS.
3.1 2D flow in a pipe↩︎
First, a 2D simulation of pipe flow was used to investigate the openHFDIBRANS sensitivity to the flow Reynolds number, chosen turbulence model, and small variations of the solid-fluid interface position. The computational domain and mesh used for simpleFoam are given in Figure 6a. Note that the computational domain is in fact 3D, but has only one cell in the third direction. The domain counterpart used for openHFDIBRANS is illustrated in Figure 6b together with the \(\lambda\) field.
3.1.0.1 Boundary conditions
At the pipe inlet, the flow was defined by a zero-gradient boundary condition for pressure and Dirichlet boundary condition for velocity \(\boldsymbol{u} = (u_\mathrm{i},0,0)^\mathrm{T}\) and turbulence variables. The inlet values of the turbulence variables were calculated according to [28] as \[\label{eq:bcspkoein} k_\mathrm{i} = \frac{3}{2}\,u_\mathrm{i}^2\,I_\mathrm{i}^2,\quad \varepsilon_\mathrm{i} = C_\nu\,\frac{k_\mathrm{i}^\frac{3}{2}}{\ell_\mathrm{t}},\quad \omega_\mathrm{i} = \frac{\varepsilon_\mathrm{i}}{C_\nu\,k_\mathrm{i}}\,,\tag{17}\] where \(I_\mathrm{i} = 2.0\,\%\) was the turbulence intensity and \(\ell_\mathrm{t}\) the turbulent length-scale determined as a tenth of the pipe inlet width.
At the pipe outlet, a Dirichlet boundary condition for pressure \(p_\mathrm{o} = 0\) and zero-gradient condition for velocity and turbulence variables were used. For walls, the zero-gradient condition for pressure, no-slip boundary condition for velocity, and wall functions for turbulence variables were considered. In openHFDIBRANS, at walls that happen to be inside the solid phase, i.e., when the wall neighboring cells have \(\lambda > 0.0\), all boundary conditions were set to zero-gradient.
3.1.1 Flow Reynolds number↩︎
First, the 2D pipe flow simulation was run with different flow Reynolds numbers which ranged from \(10^1\) to \(10^6\). The fluid kinematic viscosity was fixed at \(\nu = 10^{-6}\,\mathrm{m}^{2}\,\mathrm{s}^{-1}\) and the Reynolds number was varied by changing the inlet velocity \(u_\mathrm{i}\). For all tested \(\mathrm{Re}\), the same computational meshes were used. Consequently, the \(y^+\) value ranged from \(y^+ \ll 1\) for low \(\mathrm{Re}\) to \(y^+ \lesssim 10^3\) for high \(\mathrm{Re}\). The simulation results computed with the \(k\)-\(\omega\) turbulence model are depicted in Figure 7.
The fields compared in Figure 7 are the \(x\) component of velocity \(u_x\), the turbulent viscosity \(\nu_\mathrm{t}\), and the \(x\)-\(y\) component of the Reynolds stress tensor divided by the fluid density \(\Sigma_{xy}\). All fields were sampled along the \(y\) direction through the middle of the geometry, as indicated by the dashed black line in Figure 6. Moreover, the spatial coordinate and the field values were normalized as \[\label{eq:normal} \begin{array}{*1{>{\displaystyle}c}} \tilde{y} = \left[y - \min{(y_\mathrm{sF})}\right] / \left[\max{(y_\mathrm{sF})} - \min{(y_\mathrm{sF})}\right]\\[0.2cm] \varphi^\mathrm{rel} = \varphi/\max(\|\varphi_\mathrm{sF}\|),\quad \varphi = \{u_x,\,\nu_\mathrm{t},\,\Sigma_{xy}\}\,, \end{array}\tag{18}\] where sF stands for simpleFoam.
| \(\mathrm{Re}\) | \(10^1\) | \(10^2\) | \(10^3\) | \(10^4\) | \(10^5\) | \(10^6\) |
|---|---|---|---|---|---|---|
| \(\text{avg}\!\left[\ell_2(u^\mathrm{rel}_x)\right]\) | 2.053e-04 | 1.997e-04 | 2.021e-04 | 4.257e-03 | 1.881e-03 | 3.040e-03 |
| \(\sum \ell_2(u^\mathrm{rel}_x)\) | 4.517e-03 | 4.394e-03 | 4.447e-03 | 9.366e-02 | 4.138e-02 | 6.688e-02 |
| \(\text{avg}\!\left[\ell_2(\nu^\mathrm{rel}_\mathrm{t})\right]\) | 6.344e-04 | 8.690e-04 | 8.643e-04 | 7.516e-03 | 2.526e-03 | 2.907e-03 |
| \(\sum \ell_2(\nu^\mathrm{rel}_\mathrm{t})\) | 1.396e-02 | 1.912e-02 | 1.902e-02 | 1.653e-01 | 5.558e-02 | 6.396e-02 |
| \(\text{avg}\!\left[\ell_2(\Sigma^\mathrm{rel}_\mathrm{xy})\right]\) | 1.935e-03 | 3.285e-03 | 3.220e-03 | 2.932e-02 | 1.806e-02 | 2.395e-02 |
| \(\sum \ell_2(\Sigma^\mathrm{rel}_\mathrm{xy})\) | 4.256e-02 | 7.226e-02 | 7.085e-02 | 6.451e-01 | 3.973e-01 | 5.269e-01 |
The comparison of the field profiles showed that for all values of \(\mathrm{Re}\), openHFDIBRANS is capable of providing results that are very close to simpleFoam. To support this, an \(l_2\) norm was calculated from the difference of the profiles. Its average and total value computed for each field is given in Table 1.
| turbulence model | \(k\)-\(\omega\) SST | \(k\)-\(\varepsilon\) | realizable \(k\)-\(\varepsilon\) |
|---|---|---|---|
| \(\text{avg}\!\left[\ell_2(u^\mathrm{rel}_x)\right]\) | 3.469e-03 | 4.175e-03 | 5.593e-03 |
| \(\sum \ell_2(u^\mathrm{rel}_x)\) | 7.633e-02 | 9.185e-02 | 1.230e-01 |
| \(\text{avg}\!\left[\ell_2(\nu^\mathrm{rel}_\mathrm{t})\right]\) | 1.739e-02 | 8.542e-03 | 5.322e-03 |
| \(\sum \ell_2(\nu^\mathrm{rel}_\mathrm{t})\) | 3.826e-01 | 1.879e-01 | 1.171e-01 |
| \(\text{avg}\!\left[\ell_2(\Sigma^\mathrm{rel}_\mathrm{xy})\right]\) | 5.996e-02 | 4.056e-02 | 3.637e-02 |
| \(\sum \ell_2(\Sigma^\mathrm{rel}_\mathrm{xy})\) | 1.319e+00 | 8.924e-01 | 8.000e-01 |
3.1.2 Different turbulence models↩︎
Next, the simulation of pipe flow with \(\mathrm{Re} = 10^6\) was recomputed with different turbulence models. In addition to the \(k\)-\(\omega\) model, openHFDIBRANS can operate with the \(k\)-\(\varepsilon\), \(k\)-\(\omega\) SST and realizable \(k\)-\(\varepsilon\) models. For each model, a graphical comparison of results from openHFDIBRANS and simpleFoam is given in Figure 8 and complemented by the values of the \(l_2\) norm in Table 2. The comparison showed greater differences in the \(\nu_\mathrm{t}\) profiles compared to \(k\)-\(\omega\), see Figure 7. To investigate further, a qualitative comparison of the whole \(\nu_\mathrm{t}\) fields is given in Figure 9. This revealed that the \(\nu_\mathrm{t}\) fields given by some models, namely \(k\)-\(\omega\) SST and realizable \(k\)-\(\varepsilon\), contain abrupt changes that were not caught precisely by openHFDIBRANS. Nevertheless, for all turbulence models, the qualitative behavior of \(\nu_\mathrm{t}\) is comparable between openHFDIBRANS and simpleFoam.
| interface | (a) | (b) | (c) | (d) | (e) | (f) |
|---|---|---|---|---|---|---|
| \(\text{avg}\!\left[\ell_2(u^\mathrm{rel}_x)\right]\) | 3.968e-03 | 3.040e-03 | 1.827e-03 | 2.354e-03 | 5.099e-03 | 7.907e-03 |
| \(\sum \ell_2(u^\mathrm{rel}_x)\) | 8.730e-02 | 6.688e-02 | 4.019e-02 | 5.179e-02 | 1.122e-01 | 1.740e-01 |
| \(\text{avg}\!\left[\ell_2(\nu^\mathrm{rel}_\mathrm{t})\right]\) | 3.637e-03 | 2.907e-03 | 3.390e-03 | 4.156e-03 | 6.963e-03 | 9.614e-03 |
| \(\sum \ell_2(\nu^\mathrm{rel}_\mathrm{t})\) | 8.001e-02 | 6.396e-02 | 7.459e-02 | 9.143e-02 | 1.532e-01 | 2.115e-01 |
| \(\text{avg}\!\left[\ell_2(\Sigma^\mathrm{rel}_\mathrm{xy})\right]\) | 3.775e-03 | 2.395e-02 | 2.691e-02 | 3.119e-02 | 4.046e-02 | 4.675e-02 |
| \(\sum \ell_2(\Sigma^\mathrm{rel}_\mathrm{xy})\) | 8.305e-02 | 5.269e-01 | 5.919e-01 | 6.862e-01 | 8.902e-01 | 1.029e+00 |
3.1.3 Position of solid-fluid interface↩︎
Lastly, the 2D pipe flow test was used to investigate the effect of small variations of the solid-fluid interface position on the openHFDIBRANS results. The different interface variations that were considered are illustrated in Figure 10. For openHFDIBRANS, the interface variations result in different values of the \(\lambda\) field, but the original domain from Figure 6b is still used. For simpleFoam, the width of the middle part of the original domain from Figure 6a was adjusted to match the variations in interface position.
For each interface variation, the pipe flow simulation was run with \(\mathrm{Re} = 10^6\) and the \(k\)-\(\omega\) turbulence model. The resulting field profiles are compared in Figure 11. As can be expected, the comparison showed that the better the IB is aligned with actual mesh faces, i.e., interface positions (b), (c) and (d), the better openHFDIBRANS works. For worse aligned IB interfaces (a), (e), and (f), the openHFDIBRANS estimates are burdened with a higher error. Still, the overall robustness of openHFDIBRANS with respect to the position of IB; and, consequently, to variations in \(y_\mathrm{eff}\) seems satisfactory.
3.2 3D pipe flow↩︎
As a second test scenario, a 3D pipe flow was chosen. Following the design of the 2D domains presented in Figure 6, three 3D variants were created, see Figure 12. With each 3D variant, a simulation of flow with \(\mathrm{Re} = 10^6\) was run while using the \(k\)-\(\omega\) turbulence model. A comparison of results from simpleFoam and openHFDIBRANS is given in Figure 13 and Table 4.
| 3D pipe variation | (a) | (b) | (c) |
|---|---|---|---|
| \(\text{avg}\!\left[\ell_2(u^\mathrm{rel}_x)\right]\) | 4.858e-03 | 3.927e-03 | 8.885e-03 |
| \(\sum \ell_2(u^\mathrm{rel}_x)\) | 1.069e-01 | 1.139e-01 | 2.665e-01 |
| \(\text{avg}\!\left[\ell_2(\nu^\mathrm{rel}_\mathrm{t})\right]\) | 3.624e-03 | 5.972e-03 | 1.367e-02 |
| \(\sum \ell_2(\nu^\mathrm{rel}_\mathrm{t})\) | 7.972e-02 | 1.732e-01 | 4.102e-01 |
| \(\text{avg}\!\left[\ell_2(\Sigma^\mathrm{rel}_\mathrm{xy})\right]\) | 1.399e-02 | 1.237e-02 | 1.083e-02 |
| \(\sum \ell_2(\Sigma^\mathrm{rel}_\mathrm{xy})\) | 3.077e-01 | 3.587e-01 | 3.250e-01 |
For all the tested 3D variants, the agreement between simpleFoam and openHFDIBRANS is satisfactory. The biggest deviations appear for the last variant which has the worse alignment between IB and mesh faces, see Figure 12c (right). Yet, even with the worst aligned IB, the \(\ell_2\) norms of the profile errors, given in Table 4, are not large in magnitude.
3.3 Flow over a backward facing step↩︎
As the next test scenario, the backward facing step benchmark was chosen. A detailed description of the benchmark can be found in [29]. In addition, experimental data by [30] are available.
The flow was simulated as 2D and the utilized computational domains, the computational mesh and \(\lambda\) field are given in Figure 14. The mesh resolution of the domain was kept such that \(y^+ \in (30,200)\). Furthermore, the boundary conditions were prescribed in the same way as in the pipe flow test scenarios, see Section 3.1. Lastly, following the description in [29], the fluid kinematic viscosity was set to \(\nu = 1.56\,10^{-5}\cdot\mathrm{m}^2\,\mathrm{s}^{-1}\), and the inlet velocity \(u_\mathrm{in}\) was calculated such that \(\mathrm{Re} = 36\,000\) based on the step height \(h_\mathrm{step} = 0.0127\,\mathrm{m}\).
| \(x/h_\mathrm{step}\) | \(-4.0\) | \(1.0\) | \(2.0\) | \(4.0\) | \(6.0\) | \(10.0\) |
|---|---|---|---|---|---|---|
| \(\text{avg}\!\left[\ell_2(u^\mathrm{rel}_x)\right]\) | 4.540e-03 | 8.304e-04 | 8.665e-04 | 5.580e-04 | 5.686e-04 | 3.266e-04 |
| \(\sum \ell_2(u^\mathrm{rel}_x)\) | 3.542e-01 | 7.224e-02 | 7.538e-02 | 4.855e-02 | 4.947e-02 | 2.841e-02 |
| \(\text{avg}\!\left[\ell_2(\nu^\mathrm{rel}_\mathrm{t})\right]\) | 7.583e-03 | 1.990e-03 | 2.543e-03 | 2.596e-03 | 2.342e-03 | 2.072e-03 |
| \(\sum \ell_2(\nu^\mathrm{rel}_\mathrm{t})\) | 5.914e-01 | 1.731e-01 | 2.212e-01 | 2.258e-01 | 2.037e-01 | 1.803e-01 |
| \(\text{avg}\!\left[\ell_2(\Sigma^\mathrm{rel}_\mathrm{xy})\right]\) | 1.816e-02 | 4.038e-03 | 4.292e-03 | 3.613e-03 | 3.584e-03 | 3.022e-03 |
| \(\sum \ell_2(\Sigma^\mathrm{rel}_\mathrm{xy})\) | 1.417e+00 | 3.513e-01 | 3.734e-01 | 3.143e-01 | 3.118e-01 | 2.629e-01 |
In Figure 15 and Table 5, the results from simpleFoam and openHFDIBRANS are compared. For the comparison, data were sampled at several \(x/h_\mathrm{step}\) positions indicated by dashed lines in Figure 14d. Note that the position of the step is at \(x/h_\mathrm{step} = 0.0\).
The compared profiles are in good agreement. Only openHFDIBRANS underestimates the increase in \(\nu_\mathrm{t}\) behind the step. To investigate this, a qualitative comparison of the \(u_x\) and \(\nu_\mathrm{t}\) fields in the vicinity of the step is given in Figure 16. There, the differences in the \(\nu_\mathrm{t}\) field are apparent, but the effect of these on the \(u_x\) field is small and the qualitative characteristics of the flow are captured correctly.
![Figure 17: Evolution of the (a) wall shear stress \tau_\mathrm{w} and (b) skin friction coefficient C_\mathrm{f} sampled along the bottom wall in the backward facing step benchmark. In (b), the x range is adjusted to match the available experimental data from [30].](/article/2606.06135/main-figure16.png "fig:")
Furthermore, the evolution of the wall shear stress \(\tau_\mathrm{w} = \nu \frac{\partial u_x}{\partial y}\big|_{y = y_\mathrm{wall}}\) and skin friction coefficient \(C_\mathrm{f} = \frac{\tau_\mathrm{w}}{0.5\,u_\mathrm{in}^2}\) along the bottom wall was examined. Comparison of simpleFoam, openHFDIBRANS and experiment [30] is given in Figure 17. At the domain beginning, openHFDIBRANS underestimates the magnitude of the peak in \(\tau_\mathrm{w}\). Nevertheless, this effect does not translate to the along the step and the two profiles coincide. Hence, we suspect that the lower \(t_\mathrm{w}\) magnitude at the beginning is caused by a non-identical domain transition. Creating an identical transition would require adding an additional domain segment with a jump down as is in domains in Figures 6 and 12. However, this would be against the benchmark description given in [29] and [30], so we accepted the presented version as the best trade-off.
Near the position of the step, the profiles given by simpleFoam and openHFDIBRANS agree nicely. Also, the agreement with experiment is acceptable given that the models do not operate with resolved boundary layer since \(y^+\) is kept in the logarithmic region. Moreover, the \(C_f\) profile can be used to calculate the reattachment length \(\ell_\mathrm{reat}\). In [30], the length was located as the point with zero skin friction which was \(\ell_\mathrm{reat} = (x_\mathrm{reat} - \ell_\mathrm{step})/h_\mathrm{step} = 6.26 \pm 0.10\). When calculated based on simpleFoam data, the reattachment length is \(4.46\), and based on openHFDIBRANS, it is \(4.43\).
3.4 2D flow over a cylindrical obstacle↩︎
To evaluate the openHFDIBRANS ability to estimate drag coefficients, a test scenario with flow over a cylindrical obstacle was prepared. The utilized simulation domains together with the \(\lambda\) field used by openHFDIBRANS are given in Figure 18. The inlet and outlet boundary conditions were set the same as in the pipe flow tests, see Section 3.1, and the kinematic viscosity was \(\nu = 10^{-6}\,\mathrm{m}^2\,\mathrm{s}^{-1}\).
Meshes with three different resolutions were created. The meshes had around \(80\), \(180\) and \(320\) thousands of cells. With each mesh resolution, a series of simulations was run with Reynolds numbers ranging from \(10^{0}\) to \(10^7\) where \(\mathrm{Re} = \|\boldsymbol{u}_\mathrm{in}\|\,2\,r/\nu\) with \(\boldsymbol{u}_\mathrm{in}\) being the inlet velocity and \(r\) the radius of the cylinder. From the simulation results, the drag coefficient was calculated as \[\label{eq:dragcoeff} \begin{array}{*1{>{\displaystyle}c}} C_\mathrm{d} = \frac{2}{\rho\,A_\mathrm{ref}\,\|\boldsymbol{u}_\mathrm{in}\|^2} \int_{\Gamma_\mathrm{sf}} \boldsymbol{f}_\mathrm{tot}\cdot\boldsymbol{v}_\mathrm{d}\,\mathrm{d}S,\quad A_\mathrm{ref} = 2\,r\,\Delta z,\quad \boldsymbol{v}_\mathrm{d} = (1, 0, 0)^\mathrm{T},\\[0.3cm] \boldsymbol{f}_\mathrm{tot} = \boldsymbol{f}_\mathrm{n} + \boldsymbol{f}_\mathrm{t},\quad \boldsymbol{f}_\mathrm{n} = \rho\,\boldsymbol{n}\,S\,(\tilde{p} - \tilde{p}_\mathrm{ref}),\quad \boldsymbol{f}_\mathrm{t} = -\rho\,S\,\tau_\mathrm{w}\,, \end{array}\tag{19}\] where \(\rho\) is the fluid density, \(\Delta z\) the thickness of the computational domain, \(\boldsymbol{n}\) the surface normal, and \(\tilde{p}_\mathrm{ref}\) a reference kinematic pressure. Next, \(\Gamma_\mathrm{sf}\) represents the fluid-solid interface which is formed, in simpleFoam, by cell faces, and in openHFDIBRANS, by the immersed boundary.
The evolution of the drag coefficient with respect to \(\mathrm{Re}\) is given in Figure 19. As a reference, a compilation of experiments from [31] is added. Compared to the reference, both solvers are unable to resolve the drag crisis. Yet, this is probably caused by the use of RAS turbulence models. As was shown in [32], more advanced two layer models are required to capture the boundary layer transition. In the rest of the \(\mathrm{Re}\) range, solvers follow the experimental trends. When compared to each other, openHFDIBRANS estimates slightly lower \(C_\mathrm{d}\) values. The biggest difference is in the slope transition region around \(\mathrm{Re} = 10^2\), but there the flow around the cylinder tends to be unsteady and the RAS results are merely approximate. Outside of the unsteady region, the agreement of openHFDIBRANS and simpleFoam is satisfactory and the solver behavior does not variate much for different mesh resolutions.
3.5 Flow over a NACA-0009 airfoil↩︎
As the next test scenario, the simulation of flow over a NACA-0009 airfoil was chosen. A view of the computational mesh used by simpleFoam is given in Figure 20a, and of the \(\lambda\) field used by openHFDIBRANS in Figure 20b. The simulation was 2D and the number of cells was \(\sim\!10^{5}\).
In Figure 20a, the domain boundaries are also distinguished. At the inlet boundaries, the velocity and turbulence variables were prescribed via Dirichlet boundary conditions. For pressure, the zero-gradient boundary condition was used. The inlet values of the turbulence variables were calculated according to [33]. At the outlet, the zero-gradient boundary condition for velocity and turbulence variables, and the Dirichlet boundary condition for pressure were used. The flow was simulated as incompressible with fluid kinematic viscosity \(\nu \approx 10^{-6}\;\mathrm{m^2\,s}^{-1}\) and Reynolds number \(\mathrm{Re} \approx 2\cdot 10^{6}\) with respect to the airfoil cord length.
| \(\gamma\) | \(0^{\circ}\) | \(2^{\circ}\) | \(4^{\circ}\) |
|---|---|---|---|
| \(\text{avg}\!\left[\ell_2(u^\mathrm{rel}_x)\right]\) | 2.035e-03 | 3.020e-03 | 2.998e-03 |
| \(\sum \ell_2(u^\mathrm{rel}_x)\) | 5.250e-01 | 7.792e-01 | 7.735e-01 |
| \(\text{avg}\!\left[\ell_2(\nu^\mathrm{rel}_\mathrm{t})\right]\) | 2.634e-02 | 2.848e-02 | 2.491e-02 |
| \(\sum \ell_2(\nu^\mathrm{rel}_\mathrm{t})\) | 6.796e+00 | 7.348e+00 | 6.428e+00 |
| \(\text{avg}\!\left[\ell_2(\Sigma^\mathrm{rel}_\mathrm{xy})\right]\) | 7.902e-03 | 9.526e-03 | 9.643e-03 |
| \(\sum \ell_2(\Sigma^\mathrm{rel}_\mathrm{xy})\) | 2.039e+00 | 2.458e+00 | 2.488e+00 |
A qualitative comparison of the simulation results is given in Figure 21. Moreover, Figure 22 and Table 6 features a comparison of field profiles sampled along the dashed line given in Figure 20b. In both figures, it is apparent that openHFDIBRANS is capable of capturing the important features of the flow. However, it also tends to overestimate the size of the wake behind the airfoil.
To quantify the effect of larger wake, the resulting flow fields were used to calculate the drag and lift coefficients. The coefficients were calculated as \[\label{eq:dlcoeffs} \begin{array}{*1{>{\displaystyle}c}} C_\mathrm{d} = \frac{2}{\rho\,A_\mathrm{ref}\,\|\boldsymbol{u}_\mathrm{in}\|^2} \int_{\Gamma_\mathrm{sf}} \boldsymbol{f}_\mathrm{tot}\cdot\boldsymbol{v}_\mathrm{d}\,\mathrm{d}S,\quad C_\mathrm{\ell} = \frac{2}{\rho\,A_\mathrm{ref}\,\|\boldsymbol{u}_\mathrm{in}\|^2} \int_{\Gamma_\mathrm{sf}} \boldsymbol{f}_\mathrm{tot}\cdot\boldsymbol{v}_\mathrm{\ell}\,\mathrm{d}S\\[0.3cm] \boldsymbol{v}_\mathrm{d} = (\cos(\gamma),\sin(\gamma),0)^\mathrm{T},\quad \boldsymbol{v}_\mathrm{\ell} = (-\sin(\gamma), \cos(\gamma),0)^\mathrm{T}\,, \end{array}\tag{20}\] where the calculation of the total force \(\boldsymbol{f}_\mathrm{tot}\) is given in 19 . The reference area was calculated as \(A_\mathrm{ref} = \ell_\mathrm{c}\cdot \Delta z\) where \(\ell_\mathrm{c}\) is the chord length and \(\Delta z\) the thickness of the computational domain.
The evolution of the drag and lift coefficients with respect to the angle of attack \(\gamma\) is given in Figure 23. To have a reference, the evolution of the lift coefficient according to the inviscid thin airfoil theory [34] was added. For low magnitudes of \(\gamma\), both solvers nicely agree with \(C_\mathrm{\ell}\) given by the theory. For increasing \(\gamma\) magnitudes, the solvers slowly deviate from the theory yet this is expectable given that the solvers work with a viscous fluid. Nonetheless, openHFDIBRANS deviates faster which reflects its tendency to overestimate the size of the wake behind the airfoil. The same effect is visible on the evolution of the drag coefficient where the openHFDIBRANS profile has mostly larger slope magnitudes. But apart from the effect of the larger wake, the profile shapes are captured correctly.
3.6 Flow over an Ahmed body↩︎
As the last test scenario, the simulation of flow over an Ahmed body [35] was chosen. In Figure 24, the Ahmed body geometry is given together with its \(\lambda\) field representation, positioning within the computational domain and domain description. The simulation was constructed as 3D and the computational mesh consisted of \(\sim 2\cdot10^{7}\) cells.
At the inlet to the domain, the values of velocity and turbulence variables were prescribed, and the zero-gradient condition was used for pressure. At the outlet, the pressure value was prescribed and the zero-gradient condition was used for the rest of the variables. Next, the boundary below the Ahmed body was treated as a wall with the no-slip condition for velocity and zero-gradient condition for the other variables. At the last remaining boundaries (top, left, and right), the slip boundary condition was used for velocity and zero-gradient for the rest. The fluid kinematic viscosity was \(\nu = 1.5\cdot10^{-5}\;\mathrm{m}^2\,\mathrm{s}^{-1}\) and the flow Reynolds number \(\mathrm{Re} \approx 3\cdot10^4\) with respect to the Ahmed body length.
A comparison of the resulting \(u_x\) and \(\nu_\mathrm{t}\) fields is given in Figure 25. Similarly to the test scenario with NACA-0009, the openHFDIBRANS solver was able to correctly predict the flow patterns, but the size of the recirculation region was overestimated.
4 Conclusion↩︎
In this work, we propose our custom solver that combines the hybrid fictitious domain-immersed boundary method with Reynolds-averaged turbulence models. In particular, the HFDIB is a direct forcing IBM, and it was adjusted for the solution of the steady-state RANS equations with two-equation turbulence models and wall functions. The main strengths of the custom solver, called openHFDIBRANS, include its open-source implementation using the OpenFOAM library, along with a stability oriented design of wall functions and the incompressibility condition at the immersed boundary. The solver performance was tested in a series of benchmark cases, focusing on, e.g., its sensitivity to the flow Reynolds number, interface position and curvature, or parallelization. In all test cases, a standard body-fitted simulation run by simpleFoam served as a reference. Our openHFDIBRANS solver achieved very good agreement with the reference in many of the tested scenarios, including the prediction of lift and drag coefficients as well as wall shear stress. However, in some test cases, it tended to slightly overestimate the size of the wake behind the simulated solid objects. Still, this effect behaves predictably, and we believe that our approach is ready for usage in geometry optimizations, where the consistency of results computed for different geometries is of the most importance. In our future studies, we plan to utilize this approach in the topology optimization of engineering devices such as biosensors or molds for die casting.
The solver source code and illustrative examples are available at https://github.com/techMathGroup/openHFDIBRANS.
Acknowledgments↩︎
The work was financially supported by the Czech Science Foundation project 25-17815S. Also, the authors acknowledge the financial support provided by the Ministry of Education, Youth, and Sports of the Czech Republic via the project No. CZ.02.01.01/00/23_020/0008501 (METEX), co-funded by the European Union, and the institutional support RVO:61388998.
Nomenclature↩︎
| \(d\) | diameter, \([\mathrm{L}]\) |
| \(\boldsymbol{f}\) | momentum source term, \([\mathrm{L\,T^{-2}}]\) |
| \(I\) | turbulence intensity, \([\mathrm{-}]\) |
| \(k\) | turbulence kinetic energy, \([\mathrm{L}^2\,\mathrm{T}^{-2}]\) |
| \(\ell\) | length, \([\mathrm{L}]\) |
| \(\mathcal{M}\) | operator of the momentum transport, \([\mathrm{-}]\) |
| \(\boldsymbol{n}\) | outer unit normal, \([\mathrm{-}]\) |
| \(\mathcal{N}\) | operator of the turbulence kinetic energy transport, \([\mathrm{-}]\) |
| \(\tilde{p}\) | average kinematic pressure, \([\mathrm{L}^2\,\mathrm{T}^{-2}]\) |
| \(P_k\) | production of turbulence kinetic energy, \([\mathrm{T}^{-1}]\) |
| \(S\) | turbulence kinetic energy source term, \([\mathrm{T}^{-1}]\) |
| \(\boldsymbol{u}\) | average velocity, \([\mathrm{L\,T^{-1}}]\) |
| \(u_\tau\) | friction velocity, \([\mathrm{L\,T^{-1}}]\) |
| \(v\) | vertex, \([\mathrm{-}]\) |
| \(V\) | volume, \([\mathrm{L}^3]\) |
| \(y^+\) | normalized wall distance, \([\mathrm{-}]\) |
| \(y_\perp\) | perpendicular distance, \([\mathrm{L}]\) |
Greek letters↩︎
| \(\alpha\) | masking field, \([\mathrm{-}]\) |
| \(\gamma\) | angle of attack, \([]\) |
| \(\epsilon\) | tolerance, \([\mathrm{-}]\) |
| \(\varepsilon\) | dissipation of turbulence kinetic energy, \([\mathrm{L}^2\,\mathrm{T}^{-3}]\) |
| \(\lambda\) | solid phase indicator field, \([\mathrm{-}]\) |
| \(\nu\) | kinematic viscosity, \([\mathrm{L}\,\mathrm{T}^{-2}]\) |
| \(\sigma\) | signed perpendicular distance, \([\mathrm{L}]\) |
| \(\Sigma\) | kinematic Reynolds stress tensor, \([]\) |
| \(\omega\) | specific dissipation of turbulence kinetic energy, \([\mathrm{T}^{-1}]\) |
| \(\Omega^h\) | discrete spatial domain, \([\mathrm{-}]\) |
| \(\Omega^h_P\) | computational cell, \([\mathrm{-}]\) |
Subscripts and superscripts↩︎
| \(\mathrm{atm}\) | atmospheric |
| \(\mathrm{bc}\) | boundary cells |
| \(\mathrm{eff}\) | effective |
| \(\mathrm{f}\) | fluid |
| \(\mathrm{fc}\) | free-stream cells |
| \(\mathrm{i}\) | inlet |
| \(\mathrm{ib}\) | immersed boundary |
| \(\mathrm{ic}\) | in-solid cells |
| \(\mathrm{lam}\) | laminar |
| \(\mathrm{log}\) | logarithmic region |
| \(\mathrm{o}\) | outlet |
| \(\mathrm{old}\) | previous iteration or initial guess |
| \(\mathrm{rel}\) | relative |
| \(\mathrm{s}\) | solid |
| \(\mathrm{t}\) | turbulent |
| \(\mathrm{sf}\) | solid-fluid interface |
| \(\mathrm{vis}\) | viscous sublayer |
Abbreviations↩︎
| CFD | computational fluid dynamics |
| DEM | discrete element method |
| HFDIB | hybrid fictitious domain-immersed boundary method |
| RANS | Reynolds-averaged Navier-Stokes equations |
| sF | simpleFoam |