Introduction

Computational Fluid Dynamics (CFD) plays a critical role in the aerospace industry as it allows us to optimize the aerodynamic characteristics of aircraft, space, and other flying machines. For example, it helps develop efficient airfoils, wings, and control surfaces to reduce drag and improve lift. CFD is used to study flow patterns and combustion processes in gas turbine engines and rocket propulsion systems. It helps to optimize engine design, improve fuel efficiency and reduce emissions. CFD is used to analyze and predict heat transfer phenomena such as conduction, convection, and radiation in aerospace systems. It is important for thermal management and to ensure the structural integrity of components exposed to high temperatures.

Overall, computational fluid dynamics has revolutionized the aerospace industry, allowing engineers to gain valuable insight into complex fluid flow phenomena and optimize designs before creating costly physical prototypes. CFD has significantly reduced development time and costs while improving the safety, efficiency, and productivity of aerospace systems.

One of the main driving forces behind the growth of computational fluid dynamics was the aerospace industry. Over the past 40 years, it has evolved from a useful method of analysis to a mainstream design tool. In companies like Boeing, much of the early wing design work is done almost exclusively using CFD1,2.

It is known that turbulence is a problem of classical physics to be solved The importance of this problem lies in the fact that the vast majority of flows occurring in nature and in various technological processes are precisely of a turbulent nature. To date, there are several approaches to the mathematical modeling of turbulence. The most common is the Reynolds approach. Based on this approach, a Reynolds-averaged Navier–Stokes system of equations (RANS) is obtained. However, as is known, this system of equations is not closed. To close the resulting system of equations, a large number of different mathematical models were proposed. These models are based on the hypotheses of Boussinesq3, Kolmogorov4, Prandtl5, Karman6, etc. The NASA turbulence database7 provides a comparative analysis of various semi-empirical models. From this analysis, we can conclude that the models of Spalart and Allmaras8, and Menter k − ω SST9,10,11 have the highest ratings. To date, these models have been used to obtain numerical solutions to many important practical problems12,13,14.

An important step in the development of computational methods is the verification of the created mathematical models in wind tunnels by correcting the data obtained, excluding boundary induction15,16,17,18

However, at present, despite the fact that RANS methods are widely used, there are hydrodynamic problems the solution to which cannot give satisfactory results. These include the problem of transition from laminar to turbulent regime and separated flows.

Recently, due to the rapid development of computer technology, direct methods of turbulence simulation (DNS, LES) have become increasingly popular. These methods have high accuracy, but require large computational resources. Therefore, it will take some time to use them in solving engineering problems. The so-called hybrid RANS/LES methods, called the methods of detached-eddy simulation (DES) of vortices14, received good development. The essence of this method is that near solid surfaces, where high resolution of computational cells is required, the RANS model is used, and far from the walls, the LES model is used. The approach significantly saves computational resources and gives high-accuracy results19.

Recently, the two-fluid model of turbulence has become increasingly popular20,21. This turbulence model is based on the dynamics of two fluids, which, unlike the Reynolds approach, leads to a closed system of equations. These articles show that the two-fluid model is a low-Reynolds one and capable of describing complex anisotropic turbulent flows. In22, a two-fluid turbulence model was used to solve the problem of the transverse flow around a square cylinder. Comparison with experimental data showed high accuracy of the model.

Up to now, the numerical implementation of the two-fluid turbulence model has been conducted using proprietary computational codes. However, the model becomes more important if it is implemented using well-known software packages. To date, special programs such as ANSYS Fluent, Solidworks, Comsol Multiphysics, etc. can be used to simulate an airfoil23,24,25,26.

ANSYS Fluent is a widely used computational fluid dynamics (CFD) software package developed by ANSYS Inc. It is a powerful tool for modeling and analyzing fluid flow, heat transfer and related phenomena. ANSYS Fluent uses the control volume method to solve fluid dynamics equations.

COMSOL Multiphysics is a powerful software package for modeling physical phenomena in various disciplines, including fluid dynamics, heat transfer, structural mechanics, electromagnetism, and chemical reactions. COMSOL Multiphysics uses the Finite Element Method to solve hydrodynamic equations. COMSOL Multiphysics offers several advantages over ANSYS Fluent and Solidworks, particularly in terms of versatility, multiphysics capabilities, and customization options. In COMSOL Multiphysics, the users have the option to define and solve their own partial differential equations (PDEs) using the Custom PDE functions. This feature enables the users to model and simulate specific physical phenomena that may not be addressed by the pre-built physics modules.

An attempt is done in this article to solve this problem and the following goals are posed:

  1. 1.

    Using the Custom PDE functions to implement a two-fluid turbulence model in the COMSOL Multiphysics software package.

  2. 2.

    Validation of the two-fluid turbulence model and verification of the computational algorithm on a number of simple test problems, such as flows around a flat plate with a zero pressure gradient, a DSMA661 airfoil with an angle of attack of 0 degrees and a NACA 4412 airfoil with an angle of attack of 13.87 degrees.

  3. 3.

    Compare the obtained results with the results of the well-known SST turbulence model (built into the COMSOL Multiphysics program) and experimental data from the NASA Turbulence Modeling Resource (TMR) website7.

Two-fluid turbulence model.

The description of this model is presented in several publications by one of the authors of this article20,21,22. The main equations for studying the tasks posed are the hydrodynamic equations of the two-fluid model24 for an incompressible medium

$$ \left\{ \begin{gathered} \frac{{\partial V_{j} }}{{\partial x_{j} }} = 0, \hfill \\ \frac{{\partial V_{i} }}{\partial t} + \frac{{\partial V_{j} V_{i} }}{{\partial x_{j} }} + \frac{\partial p}{{\rho \partial x_{i} }} = \frac{\partial }{{\partial x_{j} }}\left[ {\nu (\frac{{\partial V_{i} }}{{\partial x_{j} }} + \frac{{\partial V_{j} }}{{\partial x_{i} }}) - \vartheta_{j} \vartheta_{i} )} \right], \hfill \\ \frac{{\partial \vartheta_{i} }}{\partial t} + \frac{{\partial V_{j} \vartheta_{i} }}{{\partial x_{j} }} = - \vartheta_{j} \frac{{\partial V_{i} }}{{\partial x_{j} }} + \frac{\partial }{{\partial x_{j} }}\left[ {\tilde{\nu }_{ji} (\frac{{\partial V_{i} }}{{\partial x_{j} }} + \frac{{\partial V_{j} }}{{\partial x_{i} }})} \right] + F_{si} + F_{fi} , \hfill \\ \tilde{\nu }_{ji} = 3\nu + 2\left| {\frac{{\vartheta_{i} \vartheta_{j} }}{{def(\vec{V})}}} \right|\;\;\;i \ne j,\;\;\;\tilde{\nu }_{ii} = 3\nu + \frac{{\sum\limits_{k = 1}^{3} {\vartheta_{k} \vartheta_{k} \left| {\frac{{\partial \vartheta_{k} }}{{\partial x_{k} }}} \right|} }}{{def(\vec{V})\sum\limits_{k = 1}^{3} {\left| {\frac{{\partial \vartheta_{k} }}{{\partial x_{k} }}} \right|} }}, \hfill \\ \vec{F}_{f} = - K_{f} \vec{\vartheta },\;\;\;\vec{F}_{s} = C_{s} rot\vec{V} \times \vec{\vartheta }. \hfill \\ \end{gathered} \right. $$
(1)

In the given system of equations \(V_{i}\) is the component of the average flow velocity, \(\vartheta_{i}\) is the component of the relative velocity, \(\tilde{\nu }_{ji}\) is the molar viscosity tensor, \(p\) is the pressure, \(\rho\) is the density of the medium, \(\nu\) is the molecular viscosity, Kf is the friction coefficient, \(C_{s}\) is the coefficient at the Saffman force, \({\text{def}}(\vec{V})\) is the strain rate, determined as:

$$ {\text{def}}(\vec{V}) = \sqrt {2S_{ij} S_{ij} } ,\;\;\;S_{ij} = \frac{1}{2}\left( {\frac{{\partial V_{i} }}{{\partial x_{j} }} + \frac{{\partial V_{j} }}{{\partial x_{i} }}} \right). $$
(2)

The coefficient of friction is found from the following relation:

$$ \,\,K_{f} = C_{1} \lambda_{\max } + C_{2} \frac{{\left| {d \cdot \vartheta } \right|}}{{d^{2} }}. $$
(3)

In this expression, d is the nearest distance to the solid wall, λmax is the real part of the largest root of the characteristic equation \(\det (A - \lambda E) = 0,\,\) where A is the matrix

$$ \begin{aligned} & A = \left| {\begin{array}{*{20}c} { - \frac{{\partial V_{1} }}{{\partial x_{1} }}} & { - \frac{{\partial V_{1} }}{{\partial x_{2} }} - C_{s} \zeta_{3} } & { - \frac{{\partial V_{1} }}{{\partial x_{2} }} + C_{s} \zeta_{2} } \\ { - \frac{{\partial V_{2} }}{{\partial x_{1} }} + C_{s} \zeta_{3} } & { - \frac{{\partial V_{2} }}{{\partial x_{2} }}} & { - \frac{{\partial V_{2} }}{{\partial x_{3} }} - C_{s} \zeta_{1} } \\ { - \frac{{\partial V_{3} }}{{\partial x_{1} }} + C_{s} \zeta_{2} } & { - \frac{{\partial V_{3} }}{{\partial x_{2} }} + C_{s} \zeta_{1} } & { - \frac{{\partial V_{3} }}{{\partial x_{3} }}} \\ \end{array} } \right|, \\ & \zeta = rot\overrightarrow {V} . \\ \end{aligned} $$
(4)

Constant patterns are \(C_{s} = 0.2,\;C_{1} = 0.7825,\;C_{2} = 0.306.\)

Consider a two-dimensional stationary solution to system (1). For the finite element method, the application of the standard Galerkin method will lead to a weak form:

$$ \left\{ \begin{gathered} \nu \int\limits_{\Omega } {\nabla {\mathbf{u}}:\nabla {\mathbf{v}}} + \int\limits_{\Omega } {\left( {{\mathbf{u}} \cdot \nabla {\mathbf{u}}} \right) \cdot {\mathbf{v}}} - \int\limits_{\Omega } {p\left( {\nabla \cdot {\mathbf{v}}} \right)} = \int\limits_{\Omega } {{\mathbf{f}} \cdot {\mathbf{v}}} , \hfill \\ \tilde{\nu }\int\limits_{\Omega } {\nabla {\tilde{\mathbf{u}}}:\nabla {\tilde{\mathbf{v}}}} + \int\limits_{\Omega } {\left( {{\mathbf{u}} \cdot \nabla {\tilde{\mathbf{u}}}} \right) \cdot {\tilde{\mathbf{v}}}} = \int\limits_{\Omega } {{\mathbf{f}} \cdot {\tilde{\mathbf{v}}}} , \hfill \\ \int\limits_{\Omega } {q\left( {\nabla \cdot {\mathbf{u}}} \right)} = 0. \hfill \\ \end{gathered} \right. $$
(5)

Here, the equations of system (5) are a weak form of the equation of motion for the averaged and relative velocities, as well as the equation of continuity. Here, \({\mathbf{v}}\), \({\tilde{\mathbf{v}}}\), and \({\mathbf{q}}\) are weight functions for the average velocity \(V\), relative velocity \(\vartheta\), and pressure \(p\), respectively.

Using the Finite Element Method (FEM) as a discretization method for finding solutions to turbulence models can be quite a challenge. The standard Galerkin problem statement deals with potential sources of numerical instabilities. This is, for example, the case when convection or reaction conditions prevail in the flow27,28. It is difficult to find such a stabilization that makes the solution of equations as reliable as possible. A stabilized FEM formulation can be created by adding grid-dependent, consistent, and numerically stabilizing terms to the standard Galerkin method. Various stabilization methods were proposed, many of which are based on the Petrov–Galerkin (PG) approach, which uses a modification of the standard Galerkin weight term. Examples of popular SG methods are the Streamline-Upwind/Petrov–Galerkin scheme29,30,31,32,33 (SUPG) and the pressure stabilization scheme/Petrov–Galerkin34 (SUPG). Over the years, the Petrov–Galerkin method has been developed and gave rise to new, more advanced methods. One of them is the Galerkin Least Squares (GLS) stabilization method. GLS is a general stabilization method applicable to a wide range of problems35,36,37,38,39. Its theoretical basis is that the test function should be chosen so as to minimize the squared residual of the equations. By working with matrices and not just with scalar equations, GLS can be formulated for systems of transport equations.

By adding the least squares term, we get the Galerkin/least squares formula for the two-fluid equations as shown below:

$$ \left\{ \begin{gathered} \nu \int\limits_{\Omega } {\nabla {\mathbf{u}}:\nabla {\mathbf{v}}} + \int\limits_{\Omega } {\left( {{\mathbf{u}} \cdot \nabla {\mathbf{u}}} \right) \cdot {\mathbf{v}}} - \int\limits_{\Omega } {p\left( {\nabla \cdot {\mathbf{v}}} \right)} + R_{GLS} = \int\limits_{\Omega } {{\mathbf{f}} \cdot {\mathbf{v}}} , \hfill \\ \tilde{\nu }\int\limits_{\Omega } {\nabla {\tilde{\mathbf{u}}}:\nabla {\tilde{\mathbf{v}}}} + \int\limits_{\Omega } {\left( {{\mathbf{u}} \cdot \nabla {\tilde{\mathbf{u}}}} \right) \cdot {\tilde{\mathbf{v}}}} + \tilde{R}_{GLS} = \int\limits_{\Omega } {{\mathbf{f}} \cdot {\tilde{\mathbf{v}}}} , \hfill \\ \int\limits_{\Omega } {q\left( {\nabla \cdot {\mathbf{u}}} \right)} = 0. \hfill \\ \end{gathered} \right. $$
(6)

where \(R_{GLS}\) and \(\tilde{R}_{GLS}\) are given as:

$$ \left\{ \begin{gathered} R_{GLS} = \tau \left( {{\mathbf{u}} \cdot \nabla {\mathbf{v}} - \nu \nabla^{2} {\mathbf{v}}} \right)\left( {{\mathbf{u}} \cdot \nabla {\mathbf{u}} - \nu \nabla^{2} {\mathbf{u}} - {\mathbf{f}} \cdot {\mathbf{v}}} \right) + \tau_{p} \left( {\nabla q\nabla p} \right), \hfill \\ \tilde{R}_{GLS} = \tilde{\tau }\left( {{\mathbf{u}} \cdot \nabla {\tilde{\mathbf{v}}} - \tilde{\nu }\nabla^{2} {\tilde{\mathbf{v}}}} \right)\left( {{\mathbf{u}} \cdot \nabla {\tilde{\mathbf{u}}} - \tilde{\nu }\nabla^{2} {\tilde{\mathbf{u}}} - {\mathbf{f}} \cdot {\tilde{\mathbf{v}}}} \right). \hfill \\ \end{gathered} \right. $$
(7)

Parameters \(\tau\) and \(\tilde{\tau }\) for stationary problems are determined by the following expression:

$$ \left\{ \begin{gathered} \tau = \min \left( {\frac{h}{{2\rho \left| {\mathbf{u}} \right|}},\frac{{h^{2} }}{48\nu }} \right), \hfill \\ \tilde{\tau } = \min \left( {\frac{h}{{2\rho \left| {\mathbf{u}} \right|}},\frac{{h^{2} }}{{48\;\min \left( {\tilde{\nu }_{11} ,\tilde{\nu }_{12} } \right)}}} \right). \hfill \\ \end{gathered} \right. $$
(8)

And for pressure, parameter \(\tau_{p}\) is determined by the following expression:

$$ \tau_{p} = \frac{{\left| {\mathbf{u}} \right|h}}{2}\min \left( {1,{\text{Re}}^{h} } \right). $$
(9)

where \({\text{Re}}^{h}\) is the Reynolds number of the element:

$$ {\text{Re}}^{h} = \frac{{\left| {\mathbf{u}} \right|h}}{\nu }. $$
(10)

Solution method and boundary conditions

COMSOL Multiphysics offers a range of solvers to solve various types of problems in physics. The choice of solver depends on the type of physics being modeled, the complexity of the problem, the sought-for accuracy, and the available computational resources. To solve the equations of the two-fluid turbulence model, a fully coupled approach was used with the direct solver algorithm (PARDISO). Newton's iterative method with a damping factor of 0.1 was used. The iterative process for the problem of flow around a flat plate with zero pressure gradient lasted up to 250 iterations, and for the remaining problems, the iterative process continued up to 350 iterations. The tolerance factor is 1, the residual factor is 1000.

The standard SST turbulence model uses standard COMSOL Multiphysics solvers. The following boundary conditions were set for the SST model:

$$ \omega_{\infty } = 10\frac{{U_{\infty } }}{L},\;\;\;k_{\infty } = 0.1\frac{{\nu_{\infty } U_{\infty } }}{L}. $$
(11)

where L is the characteristic size of the streamlined body. The remaining boundary conditions were specified in the standard way.

For the problem of an aerodynamic profile, the initial distribution of velocity and pressure was given by the potential field of velocities. Assuming an irrotational inviscid flow, \(\varphi\) the velocity potential is defined as

$$ {\mathbf{u}} = - \nabla \varphi . $$
(12)

The velocity potential must satisfy the continuity equation for an incompressible flow \(\nabla {\mathbf{u}} = 0\). The continuity equation can be represented as the Laplace equation

$$ \nabla \left( { - \nabla \varphi } \right) = 0. $$
(13)

which is the potential flow equation.

After calculating the velocity potential, the pressure can be approximated using the Bernoulli equation for stationary flows:

$$ p = - \frac{\rho }{2}\left| {\nabla \varphi } \right|^{2} . $$
(14)

Flat plate with zero pressure gradient

The main purpose of this experiment is to test the implementation of the two-fluid turbulence model in Comsol Multiphysics and compare the obtained results with the experimental data presented on the NASA website7. In the calculation, the plate length was L = 2 m with Reynolds number Re = 5,000,000 per unit length. In this case, the maximum thickness of the boundary layer is approximately 0.03 L, so the computational grid height was removed by distance y = L, which is sufficient to have little effect on the results. The turbulence intensity of the oncoming flow was 0.039%. The boundary conditions are shown in Fig. 1a, and an illustration of the computational grid and domain are shown in Fig. 1b. A computational grid of 69 × 49 in size, presented on the NASA website was used7.

Figure 1
figure 1

Flat plate with zero pressure gradient (a) boundary conditions and (b) computational grid and domain.

Full size image

Below are comparisons of the obtained numerical results with known experimental data. Figure 2 shows: a) the dependence of the friction coefficient along the plate; b) the dimensionless longitudinal flow velocity as a function of the dimensionless distance to the plate, as well as the results of Cole’s theorie40,41.

Figure 2
figure 2

Dependence of the coefficient of friction along the plate (a), profile of the dimensionless longitudinal flow velocity on the dimensionless distance to the plate (b).

Full size image

Here \(C_{f}\) is the coefficient of friction of the plate:

$$ C_{f} = \frac{2}{{\text{Re}}}\left( {\frac{{\partial V_{x} }}{\partial y}} \right). $$
(15)

Figure 3 shows the profiles of the dimensionless longitudinal velocity at different Reynolds numbers in two sections: (a) x = 0.97 m and (b) x = 1.97 m.

Figure 3
figure 3

Dimensionless longitudinal flow velocity at various Reynolds numbers.

Full size image

The solid line in Fig. 3 shows the results of numerical calculation for the dimensionless longitudinal flow velocity depending on the dimensionless distance to the plate. Dimensionless velocities and distance were determined by the following formulas:

$$ u^{ + } = \frac{{V_{x} }}{{u^{*} }},\;\;\;y^{ + } = {\text{Re}} yu^{*} ,\;\;\;u^{*} = \sqrt {0.5C_{f} } . $$
(16)

Figures 2 and 3 show that, for all Reynolds numbers, the model well describes both laminar and turbulent zones22.

Airfoil DSMA661

The DSMA661 airfoil is designed for low Reynolds flow. It was developed by the Delft University of Technology in the Netherlands. The DSMA 661 airfoil is commonly used in unmanned aerial vehicles (UAVs), especially those operating at low speeds. This airfoil has a relatively large maximum thickness of 16% of the chord length. It has moderate camber and is designed to provide good lift and low drag at low Reynolds numbers. The DSMA 661 airfoil is known for its stable and predictable behavior, making it suitable for a variety of UAV missions.

The DSMA661 airfoil42 provides another opportunity to evaluate the implementation of the proposed two-fluid model in Comsol Multiphysics. Four different two-dimensional grids were used in the work. Each coarser grid represents exactly every second point of the next finer grid, ranging from the 1121 × 193 fine grid to the coarsest 141 × 25 grid, presented on the NASA website7. There are 305 points on the fine grid along the trace from the trailing edge of the profile to the outflow boundary (39 points on the coarsest grid). The main results were taken for a grid of 561 × 97, and the remaining computational grids were used to check the grid convergence. The boundary conditions are shown in Fig. 4a, and an illustration of the computational grid and computational domain is shown in Fig. 4b.

Figure 4
figure 4

Airfoil DSMA661 (Model A). (a) boundary conditions; and (b) computational grid and area.

Full size image

For the problem at hand, the results of the well-known SST turbulence model are also obtained. The results are obtained with a free stream flow velocity of 18 m/s at zero angle of attack, which corresponds to a Reynolds number based on a chord of 1.2 million. For both models, the free flow turbulence intensity was 0.088%.

The distribution of the surface pressure coefficient on an airfoil is characterized by a change in pressure on its surface depending on the distance from a certain point. Typically, the analysis uses the surface pressure coefficient Cp, defined as the ratio of the pressure difference between a point on the surface of the profile and the pressure of the free flow to the dynamic pressure of the free flow.

$$ C_{p} = \frac{{p - p_{\infty } }}{{0.5\rho U_{0}^{2} }} $$
(17)

where p is the pressure at a point on the profile surface, P∞ is the pressure of the free flow, ρ is the density of the free flow, U0 is the velocity of the free flow.

The distribution of the surface pressure coefficient on an airfoil can be used to analyze its aerodynamic characteristics such as lift, drag coefficient, etc. Figure 5a shows the distribution of the surface pressure coefficient Cp of the DSMA661 profile at an angle of attack \(\alpha = 0^{0}\).

Figure 5
figure 5

Research results (a) pressure coefficient (b) friction coefficient.

Full size image

The distribution of the coefficient of skin friction on an airfoil is characterized by a change in the friction force on its surface depending on the distance from a certain point. The coefficient of skin friction Cf is defined as the ratio of the friction force acting on the surface of the profile to the dynamic pressure of the free flow.

$$ C_{f} = \frac{F}{{0.5\rho U_{0}^{2} S}} $$
(18)

where F is the friction force acting on the profile surface, S is the profile surface area oriented along the flow. Figure 5b illustrates the numerical and experimental results for Cf for both the top and bottom profile surfaces.

Figures 5 show that the results of both turbulence models practically coincide and are in good agreement with the experimental data.

Table 1 presents the errors in the deviation of numerical results from experimental data for Cf and Cp.

Table 1 Advances in LES of Complex Flows.
Full size table

Figures 6 and 7 show the longitudinal velocity \(U/U_{0}\) profiles along the top (Fig. 6) and along the bottom (Fig. 7) surfaces of the profile at different sections along the flow.

Figure 6
figure 6figure 6

Longitudinal velocity profiles on the top surfaces at x/c = 0.493, 0.593, 0.693, 0.793, 0.893, 0.94, 0.97, 0.99, 1.

Full size image
Figure 7
figure 7figure 7

Longitudinal velocity profiles on the bottom surfaces at x/c = 0.493, 0.593, 0.693, 0.793, 0.893, 0.94, 0.97, 0.99, 1.

Full size image

As can be seen from Figs. 6 and 7, the results of both models are close to the experimental results.

Below are the numerical results of the models and experiment for the profiles of longitudinal velocity \(U/U_{0}\) (Fig. 8) and turbulent stresses \(\overline{{u^{\prime}\vartheta^{\prime}}} /U_{0}^{2}\) (Fig. 9) in the wake after the profile at sections x/c = 1.01, 1.05, 1.20, 1.40, 1.80 and 2.19.

Figure 8
figure 8

Longitudinal flow velocity in the turbulent wake in sections x/c = 1.01, 1.05, 1.20, 1.40, 1.80 and 2.19.

Full size image
Figure 9
figure 9

Wake turbulent stresses in sections x/c = 1.01, 1.05, 1.20, 1.40, 1.80 and 2.19.

Full size image

Table 2 presents the errors in the deviation of the calculated results from the experimental data in the cross section x/c = 2.19 for U/U0.

Table 2 Deviation of tu.
Full size table

Figure 8 shows that the results of both models coincide with each other and with the experimental results for the average computational grid 561 × 97 with high accuracy.

For turbulent stress \(\overline{u^{\prime}\vartheta ^{\prime}} /U_{0}^{2}\), there is also a good agreement between the results of models and experiment.

Figure 10 shows the contours for longitudinal velocity and turbulent stresses.

Figure 10
figure 10

Loops for longitudinal velocity (m/s) and turbulent stresses (m2/s2).

Full size image

To check the grid convergence, Fig. 11 shows the numerical results when the computational grid is changed. The results in section x/c = 0.99 for the bottom surface of the profile are shown.

Figure 11
figure 11

Checking the grid convergence of turbulence models.

Full size image

In Fig. 11 shows that the results of both models do not depend on the resolution of the computational grid.

Figure 12 shows the change in the lift coefficient CL from the resolution of the computational grid.

Figure 12
figure 12

Change in the lift force coefficient CL from the resolution of the computational grid.

Full size image

It is also clear from Fig. 12 that the results of both models are less sensitive to the resolution of the computational grid.

Airfoil NACA 4412

The NACA 4412 airfoil is an airfoil with a maximum thickness of 12% of the chord length, located at 40% of the chord length. It is widely used in various fields, including aircraft wings and propeller blades. This section presents the validation of a two-fluid turbulence model for the NACA 4412 airfoil43,44. The case is considered for the Mach number M = 0.09, the angle of attack α = 13.87 ◦ and the Reynolds number Re = 1,520,000. 449 × 129 grid is presented in7. The boundary conditions are shown in Fig. 13a, and an illustration of the computational grid and domain is shown in Fig. 13b.

Figure 13
figure 13

NACA 4412 airfoil. (a) Boundary conditions and (b) computational grid and domain.

Full size image

Figure 14 shows the isolines of the flow velocity.

Figure 14
figure 14

Isolines of flow velocity.

Full size image

It can be seen from this figure that both models show close results to the experimental data.

Figure 15 shows the distribution of the surface pressure coefficient Cp on the NACA 4412 airfoil. From this figure it can be seen that the results of the two-fluid model are in better agreement with the experimental data. This can be explained by the fact that near the edge of the profile as shown in Fig. 14, a recirculating flow movement occurs, which causes anisotropic turbulence. The SST model uses the Boussinesq hypothesis, which is valid for isotropic turbulence. Therefore, under anisotropic turbulence, the SST results are somewhat worse. As for the two-fluid turbulence model, as shown in previous works, it is capable of describing anisotropic turbulence with great accuracy.

Figure 15
figure 15

Results for surface pressure coefficient.

Full size image

Table 3 presents the deviation of model results from experimental data for Cp.

Table 3 The relative.
Full size table

The controls of the airliners are mainly located on the edge of the profile. Therefore, the accuracy of calculating the distribution of the pressure coefficient near the profile edge is of great practical importance. In this regard, Fig. 16 shows the distribution of the pressure coefficient at trailing of the edge. From this figure it can be seen that the results of the two-fluid model are significantly better than the results of the SST model.

Figure 16
figure 16

Pressure coefficient distribution at trailing of the edge.

Full size image

Figure 17 shows the same results using different models for the NACA4412 airfoil at an angle of attack of 120 and Reynolds number Re = 1.64·106. These results are obtained from the work45.

Figure 17
figure 17

Pressure coefficient distribution near the profile edge.

Full size image

The results presented in Figs. 16 and 17 show that the accuracy of the two-fluid model is higher than that of the RANS models and no worse than that of the LES and DES models.

Figure 18 shows the longitudinal velocity \(U/U_{0}\) profiles along the upper surface of the NACA 4412 airfoil at different sections downstream.

Figure 18
figure 18

Longitudinal velocity profiles on the upper surface of the profile at x/c = 0.6753, 0.7308, 0.7863, 0.8418, 0.8973, 0.9528.

Full size image

Here, too, there is a correspondence between the results of the models and experimental measurements. It can be seen from the figure that the velocity profiles according to the results of the SST model slightly deviate from the experiment in all sections, and the results of the two-fluid model are in good agreement with the experimental data.

Table 4 presents the relative deviations of the model results from the experimental data in the cross section x/c = 0.6753 for U/U0.

Table 4 The Cp model.
Full size table

Figure 19 shows the turbulent stress \(\overline{{u^{\prime}\vartheta^{\prime}}} /U_{0}^{2}\) along the top surface for the NACA 4412 airfoil in different sections.

Figure 19
figure 19

Profiles of turbulent stresses on the upper surface in sections x/c = 0.6753, 0.7308, 0.7863, 0.8418, 0.8973, 0.9528.

Full size image

From Fig. 19 it is clear that the turbulent stresses correspond to the experimental data worse than the averaged values. This is due to the fact that they have small values and to improve compliance it is necessary to use calculation schemes of a higher order of accuracy.

Calculation time for problem NACA 4412: 420 iterations and 1740s in the SST model, and 380 iterations and 1588 s in the two-fluid model. All calculations were performed on a computer with a 2.5 GHz quad-core Intel i5-7300HQ processor, 16 GB DDR3 RAM, a 1024 GB hard drive, and Windows 7 (64-bit).

Conclusion

This article demonstrates the capability of the two-fluid turbulence model in the Comsol Multiphysics software package, which uses the finite element method. To stabilize the model equations, a stabilizer based on the Galerkin least squares method was used. For validation of the model and verification of the numerical algorithm, the problems of flow past a flat plate with a zero pressure gradient, as well as DSMA661 and NACA 4412 airfoils with angles of attack of 0 and 13.87 degrees, respectively, are considered. For the first time in the Comsol Multiphysics program, a two-fluid model was introduced and the results were obtained. To verify the created program, the results obtained are compared with the results of other models, as well as with experimental data. It is shown that the results of the two-fluid model correspond better to experimental data than the results of other RANS models and are no worse than the results of the LES and DES models. In addition, it is shown that the two-fluid model requires less computational resources than the SST model. Therefore, the two-fluid model can be recommended for calculating engineering problems of turbulent hydrodynamics.