Introduction. Submerged combustion apparatus (SCA) are thermal engineering devices for the direct heating and evaporation of solutions. The mixture is heated via direct contact with a jet of hot flue gases generated in the combustion chamber. The key advantage of this heating method is the absence of heat-transfer surfaces, where salt precipitation and scale deposition typically occur. During operation of the SCA, the flue gas jet contacts the liquid and breaks into bubbles, which bubble through the liquid layer and transfer heat to it. Turbulent flow occurs in the liquid due to the high-velocity inlet of the flue gas jet. Therefore, the flows in the apparatus are calculated using one of the turbulence models. Several classes of such models are known [1].

Figure 1 shows a general classification scheme for turbulence models. Most of them belong to the RANS class (Reynolds-averaged Navier-Stokes equations). The essence of averaging is to replace the true velocity with the ensemble-averaged velocity [2]:

 (1)

Here, v — true velocity; 〈〉 — average velocity; v′ — velocity fluctuation. The Navier–Stokes equation takes the form:

 (2)

  (3)

Equation (3) represents the Reynolds stress tensor. Equation (2) remains open until the tensor components are determined. RANS class models define the Reynolds stress tensor.

A subclass of models using the Boussinesq hypothesis [3], is based on the proposition of proportionality of the stress tensor and the strain rate tensor. In this case, the proportionality coefficient is turbulent viscosity, and its determination becomes the key to describing the liquid motion. A number of models have been proposed for calculating turbulent viscosity, differing in the number of differential equations (DEs) used. The following models are known:

– without DE (zero-order models described in [4] and [5]);

– with one DE [6];

– with two and more DE (k – ε [7], k – ω [8], k – ε – n² [9]).

An alternative turbulence modeling within the RANS framework involves a departure from the Boussinesq hypothesis [10]. From a mathematical standpoint, the Reynolds stresses represent the second-order single-point moments of the velocity fluctuations, and the transport equations for these quantities can be rigorously derived from the Navier–Stokes expressions:

 (4)

Here, Dijk — diffusion term, Pij — generation term, Fij — diffusion term, eij — dissipative term. Determination of the relationship between the terms of equation (4) is the essence of the problem of constructing Reynolds stress transfer models. A set of algebraic and differential equations can be used.

Eddy-resolving models, instead of the stress averaging procedure, implement a variable filtering procedure [11]:

 (5)

Large-scale turbulence eddies are resolved directly using filtered variables. Small-scale eddies, smaller than the grid cell size, are resolved using subgrid-scale models. The most well-known of these is large-eddy simulation (LES). In the limit of a vanishing grid size, LES reduces to direct numerical simulation (DNS) of turbulence, which is currently used only for fundamental research.

The combination of RANS and eddy resolution is the basis for a class of hybrid models. Among them, we note the detached eddy simulation (DES) model [12], as well as the group of scale-adaptive simulation (SAS) models [13]. This diversity, on the one hand, allows for the flexible selection of a turbulent model suitable for a specific computational case, and on the other hand, requires careful analysis and verification of the models.

For numerical experiments, an abstract representation of the SCA is used — a cell with a submerged jet and a free surface [14]. Its diagram is shown in Figure 2.

From a computational fluid dynamics (CFD) perspective, this cell is a challenging computational case. The complexity stems from five factors:

Recommendations for describing these factors contradict each other1. A multi-velocity continuum is proposed for modeling bubble flows, and a single-velocity continuum for modeling the free surface. Locally unstable flows require the use of eddy-resolving turbulence models, while laminar and turbulent regions require resolution of the laminar-turbulent transition. Researchers of gas-lift [15] and bubbling [16] systems use the k – ε model, but in such systems, there is no high-speed injection of the gas-liquid jet. Thus, it is highly challenging to make a definitive choice of turbulence model for the cell under consideration based solely on previous publications. At the same time, the national economy requires reliable technology for evaporating concentrated brines. These brines are either waste from potash fertilizer production or an intermediate product from underground salt leaching.

Materials and Methods. The following turbulence models were considered for the analysis:

These models can be considered generally applicable, as they are all included in the engineering handbook for chemical industry workers [17]. The zero-order model, obviously, is not intended for solving this class of problems. It was included in the experiment for indirect verification of the numerical model. Verification will be considered successful if the calculation using the zero-order model yields unphysical results or end with an error.

The models under consideration implement different approaches to turbulence modeling:

Furthermore, these models are integrated into the following computational packages: CFX, Fluent, FlowVision, and Logos.

The equations of the turbulence models under consideration were included in a general physicomathematical model of the hydrodynamic system under study. The general model describes a multiphase flow in which the carrier phase is a liquid, and the dispersed phase consists of gas bubbles and solid particles. The Euler–Euler approach [18], was used to calculate the gas-phase flow, whereas the Euler–Lagrange approach [19] was applied for the solid phase. Interphase interaction was calculated using empirical relationships:

Gas bubbles with a constant diameter of 3 mm are introduced into the combustion chamber nozzle exit. Gas density is variable, which is a consequence of strong temperature variations in the gas phase along the flow line.

Bubble coalescence was not calculated. Air was chosen as the gas, as the thermal properties of flue gases were close to those of the flue gas at an excess air ratio of 1.7. The solid phase inlet is located at the hypocenter of the combustion chamber, some distance from the nozzle exit. This distance is selected so that the inlet region is located within the gas-liquid jet, in the zone of maximum thermal stress. In the actual process, the liquid evaporates precisely at this point, which consequently results in salt precipitation. An open boundary allows liquid to escape from the computational domain, and gas to penetrate into it and partially escape. In this section of the boundary, the parameters of the gas phase beyond the computational domain are efficiently set: all walls are impermeable and adiabatic, plus no-slip conditions. Additionally, for the bottom wall (Fig. 1), a condition for complete absorption of the solid phase is introduced. This simulates the settling of solid particles on the bottom of the apparatus.

The equations of the physicomathematical model were solved by the finite volume method with the computational fluid dynamics (CFD) software package Ansys CFX 2019 R3.

Description of the mathematical model. The mathematical model of the cell is based on the Eulerian-Eulerian and Eulerian-Lagrangian equations. The equations of continuity, conservation of momentum, and energy are presented below:

 (6)

 (7)

 (8)

The viscous stress tensor from (7) and (8) is calculated using the corresponding equations of the turbulence models mentioned above.

Law of motion for Lagrangian particles:

  (9)

Forces in (7) and (9):

In the presented equations, α — volume fraction; ρ — density, kg/m3; v — velocity, m/s; t — time, s; p — pressure, Pa; H — specific enthalpy, J/kg; λ — thermal conductivity, W/m∙K; T — temperature, K; Q — amount of heat transferred from phase j, J/m3; cp — heat capacity at constant pressure, J/kg∙K; d — diameter, m; m — mass, kg. Vectors , , , represent the acceleration due to gravity, m/s², specific force of interphase interaction, N/m3, specific force of added mass, N/m3, and the force of gravity, N; τ — viscous stress tensor. Subscripts i, j indicate arbitrary phases. Subscripts g, l, s correspond to the gas, liquid, and solid phases, respectively. CVM — coefficient of added mass, set to 0.5. CD — drag coefficient calculated using the Schiller–Naumann correlation for solid particles or the Ishii–Zuber correlation for gas bubbles.

The boundary and initial conditions described below are added to the basic equations.

The first is the combustion chamber nozzle exit (the conditions correspond to the experimental setup [23]):

Second — open border:

Third — liquid inlet:

Fifth — initial conditions:

The free surface of the liquid is bounded by plane z = 0.8 m.

The following dimensions of the computational domain were adopted: length and width – 1.3 m, height – 1 m. The number of computational grid nodes in three dimensions was 130:130:100. A uniform tetrahedral grid was used. To verify the independence of the calculation results from the grid cell size, numerical experiments were conducted using grids of varying sizes. The results of these experiments are presented in Table 1.

As can be seen from Table 1, the average energy of turbulent pulsations (a key parameter for this study) for grids with 15- and 10-mm cells differs by 0.001 J/kg, which is 2% of the obtained values. Therefore, both the 15- and 10-mm grids can be used for further calculations. However, the 10-mm grid demonstrated greater calculation stability, so it was chosen for further work.

Numerical experiments were conducted in steady-state and transient modes. The time step in the transient mode was 0.1 s. In both modes, the single-continuum (homogeneous) and multi-continuum (heterogeneous) formulations were alternately applied for turbulence modeling. The remaining physical quantities were computed using the heterogeneous formulation.

Research Results. The following data were obtained as a result of the computational experiments:

In the tables, the type of data that could not be obtained during the experiment is designated as n/a. Symbol х indicates cases where the calculation was not completed due to a solver error.

Figures 3 and 4 present the liquid velocity fields with the gas bubble streamlines overlaid on them.

Additionally, the accumulation of solids at the bottom of the cell in a transient mode was studied. The mass of solid particles touching the lower boundary of the computational domain over the entire duration of the numerical experiment (60 s) was calculated. The data are presented in Table 6.

As can be seen, the mass of the settled solid phase is in the range of 1.18–1.38 kg, indicating a low scatter of values. The difference between the two SST model formulations is insignificant (6%). The homogeneous SST and DES models yielded the smallest difference in results (3%), despite the differences in flow patterns.

Discussion. A comparison of the specific turbulent kinetic energy values (Table 2) shows that they are all of the same order of magnitude, ranging from tenths to hundredths of Joules per kilogram. The homogeneous formulation yields higher fluctuation energy values than the heterogeneous one.

The average liquid velocities (Table 3) are also of the same order of magnitude — tenths of a meter per second. The spread of pulsation energies and velocities is small.

The liquid volume values (Table 4) in the computational domain are close to each other for the steady-state mode. The exception is the heterogeneous SST model, with a volume value 0.04 m³ lower than the average for the other models. In the transient mode, the DES and homogeneous SST models yield volume values of approximately 1.19 m³, while the other two yield values of approximately 1.3 m³. This indicates a low degree of conservatism for the first two models, meaning they may yield unacceptable results over long-time scales.

The distribution of mean temperature values (Table 5) shows a different trend. Thus, in the transient analysis, the difference between the two SST model formulations is minimal (0.7°C), while the difference between the SST and eddy-resolving models is several degrees. The steady-state mode is characterized by a spread in mean temperature values. The minimum (28.2°C) is given by the k – ω model, and the maximum (58.7°C) — by the heterogeneous k – ε.

It should be noted that all calculations using the zero-order model resulted in errors, as expected. This is entirely consistent with expectations and indicates indirect verification of the numerical model.

As can be seen from Tables 2–5, twelve out of twenty-four calculations (excluding calculations with the zero-order model) resulted in errors. This indicates low computational stability for the numerical model containing both a free surface and a dense bubble flow. This instability is due to the need to calculate the position of the free surface in a heterogeneous formulation. As mentioned earlier, in a homogeneous formulation in a number of numerical experiments, only turbulence is modeled. To simulate flows with a free surface, it is recommended to use a complete homogeneous formulation2. However, the need to calculate bubble flow precludes the use of a homogeneous approach.

The steady-state calculation mode turned out to be more stable than the transient one. This can be explained by the accumulation of numerical errors in the transient solution [24]. This effect is especially noticeable for multiphase problems, where a larger number of equations (a multiple of the number of phases) are solved. Thus, the steady-state mode can be recommended for modeling stabilized modes, when the time dependence of variables is not significant.

Examining the velocity fields obtained in the transient mode (Fig. 3), it can be observed that the DES and LES models produce a flow structure with a greater number of eddies than the SST. The free surface of the liquid turned out to be strongly deformed in the calculations with the DES and SST models within the homogeneous formulation.

Examining the velocity fields obtained in the steady-state mode (Fig. 4), we see that none of the models produces significant free surface deformation, even the SST. The latter yields a flow pattern with a slightly more complex free surface relief. The result obtained using the k – ω model is worth noting. Here, the free surface is smoothest, and the velocity field is symmetrical, which is essentially reminiscent of the reference flow pattern in the SCA [25]. The remaining models gave similar results in both homogeneous and heterogeneous formulations (Fig. 4 а–f, h): a two-eddy flow pattern with a predominant left-to-right flow direction across the computational domain.

The liquid velocity fields correspond to regions where the volume fraction of liquid exceeds 85%. This allows assessing the state of the free surface. In all the illustrations, we see a similarly shaped region of the gas-liquid jet. Its range is virtually identical for all the models and computational modes considered. The jet diameter on the free surface varies slightly, except in cases where the free surface is highly curved (Figs. 3 a and 3 c). The liquid velocity vectors are directed tangentially to the gas phase flow lines. This is consistent with previously obtained results from studies of the hydrodynamic situation in the SCA [23].

Let us analyze the flow structure in the computational domain. The problem formulation itself is asymmetric: the liquid inlet is located on the right side of the computational domain, and the liquid and flue gas outlet is on the left. This reflects the design of real SCA, where a symmetric inlet and outlet configuration of the working media is impractical and even technically impossible. Consequently, the velocity field in the cell must also be asymmetric.  

The mirror asymmetry of all velocity fields, except for that obtained using the k – ω, was noted above. This is explained by the asymmetric eddy structure of the flow (three eddies cannot be symmetrical relative to the vertical cross-sectional plane). This is due to the use of eddy-resolving models. Therefore, the feature in question is caused by the inlet and outlet flows, on the right and left, respectively.

The SST model is a combination of k – ω and k – ε, and the contribution of each part is determined by a weighting coefficient [26]. The qualitative similarity of the SST and k – ε modeling results is explained by the prevailing weight of the k – ε component.

The EARSM model uses a different Reynolds stress averaging procedure and can therefore serve as a benchmark for comparing RANS-class models. According to Figure 4h, the flow structure reproduced by the EARSM model is close to that of the k – ε and SST models, confirming the validity of the results obtained using them. It can be concluded that the k – ω model yields an unphysical result, since strict symmetry of the velocity field cannot exist in an asymmetric problem.

It was previously noted that the homogeneous DES model reproduces the maximum asymmetry of the free surface of the liquid. The obtained inclination angle (about 30º) does not correspond to the real hydrodynamic situation in the SCA; therefore, the authors consider this unphysical result to be an extreme case and deem it inappropriate to use a homogeneous DES model.

When discussing the reliability of various models, it is worth highlighting the similarity between most of them in both qualitative and quantitative results. Let us consider three extreme cases:

The remaining models are generally equivalent. However, for steady-state conditions, the SST model is recommended. It has demonstrated the greatest computational stability and allows for resolving the laminar-turbulent transition. For cases where process dynamics must be considered, the LES and heterogeneous SST models are recommended.

Conclusion. Numerical simulations of the SCA were performed using various methods for resolving turbulent fluctuations. The results of the numerical experiments suggest a similar flow pattern for all turbulence models. Similarities were found across four criteria:

All models, except the heterogeneous RNG k – ε, k – ω and DES, demonstrate insignificant scatter in integral parameters such as turbulent pulsation energy, average velocity, temperature, and liquid volume. These same models reproduce the asymmetry of the forced liquid flow in the apparatus. The k – ω model proved insensitive to asymmetry, and DES reproduced a free surface with an inclination angle of 30º, which does not correspond to the hydrodynamics of the SCA.

Numerical experiments demonstrated low computational stability in a problem with both free surface and dense bubble flow. Half of the calculations failed due to a solver error. However, it is found that the steady-state mode is more stable. This allows us to recommend the steady-state mode for cases where parameter changes over time are not of interest. This happens, for example, when studying a steady-state mode or when there is no need to take into account the accumulation of Lagrangian particles. Furthermore, SST can be recommended for SCA calculations. This model has demonstrated the greatest computational stability and also allows for resolving the laminar-turbulent transition.

Further research involves a detailed study of the reasons why turbulence models reproduce such extreme results as elevated mean temperature or significant curvature of the free surface.

1. Ansys CFX Reference Guide. Release 2021 R1. Canonsburg: ANSYS Inc; 2021. 450 p.2. Ansys CFX Tutorials. Release 2024 R2. Canonsburg: Ansys Inc., 2024. 1002 p.