Numerical Modeling of Jet at the Bottom of Tank at Moderate Reynolds Number Using Compact Hermitian Finite Differences Method

Abstract

In this manuscript, the injection of a homogeneous jet in a numerical tank is considered to revolve around discussing the limitation of the direct numerical simulation (DNS), to resolve the equations governing the problem of a jet emitted from the bottom of a numerical tank. The investigation has been made in the context of an unsteady, viscous, and incompressible fluid. The numerical resolution of the equations governing the problem is made by the compact Hermitian finite differences method (HFDM) high accuracy $O\left(h^2,h^4\right)$ First, the numerical code used in this work is validated by comparing the profiles of the velocity components at the median of the lid-driven cavity with the results of the literature. Furthermore, to confirm the validity of the present numerical code, an evaluation of mesh domain sensitivity is assessed by comparing the numerical vertical velocity profiles for different steps of y-direction (flow direction) with the analytical solution. Afterward, the aim is to perform the nonlinear simulations of the Navier–Stokes equations in a large computational domain. Next, the goal is to characterize the instabilities associated with high Reynolds numbers when a jet is emitted from the bottom of the numerical tank.

1. Introduction

The importance of the numerical studies of the dynamics of jets for many industrial applications, has been attracting many engineers and researchers to investigate the instabilities of jet flows. Many issues are of interest to the original appearance of the excitability of turbulent structures. This is the reason why, researchers and engineers are classified into three main groups: first, are in interest to study the effect of nozzle geometry that has been paid attention, and proven significant influence on the jets profiles [1,2,3,4]. The second group focuses on discussing the effect of the initial conditions on the behavior of the flow of jets [5,6,7,8,9,10]. The third categories give importance to the effect of Reynolds number investigations. Based on the state-of-the-art, the study of the effect of the geometry of nozzles, and initial conditions on the stability of jets, are still lacking conversely to the investigations on the Reynolds number effects.

The Reynolds number which is a number that measures the ratio between convection and diffusion can characterize the nonlinearity of the equations governing the jet flows. Numerous studies that have substantially explored the influence of Reynolds number on jet instabilities. From the work of Reynolds [11], Crow and Champagne [12], O’Neill [13] have been in interest to investigate the onset of instability and the growth of structures in the shear layer. Further, many authors have shown that numerical instabilities are associated with increasing the Reynolds number. In the case of round jets, Bogey and Bailly [14] showed the appearance of vertical structure closer to the nozzle as Reynolds number increases, Dimotakis [9] studied the effect of Reynolds number on the flow structures, and investigate the critical value of Reynolds number which makes the flow properties a weak function of Reynolds number. Experimentally, Ricou [15] has shown that the weak dependence of Reynolds number is related to the nozzle exit diameter. Furthermore, the square jet flow characteristics have been investigated by Grinstein et al. [16], Grinstein and DeVore [17], Quinn and Militzer [18], Sankar et al. [19], Tsuchiya et al. [20] and more recently by Ghasemi et al. [21,22], proved that the initial vortex structures are the main difference between round and square jets.

Additionally, it is worth pointing out the methods used to resolve equations governing jet flows. First, the Reynolds Averaged Navier Stokes (RANS) approach may be the most employed approach for the reason of considerable turbulent models within the formulation of turbulence equations. Even though, in all engineering applications, it is necessary to validate at the first the turbulence models used to guarantee accurate solutions [23]. Even though, RANS mode is insufficient to predict transition flow behavior. In the second place, an alternative approach contains most of the turbulent energy and responsible for most of the momentum transfer was adopted, this approach called; large-eddy simulation (LES) is more accurate than the RANS approach. Nevertheless, the most expensive one for simulating the jet flows even it is applied only to low Reynolds number flows, is the direct numerical simulation (DNS) in which the full Navier–Stokes equations are solved numerically [24].

In light of the state-of-the-art above, the evolution and the distribution of the concentration of jets flows have been attracting attention in many engineering applications [25,26,27,28,29]. For this reason, the objective at first and foremost is to study accurately the dynamic phenomenon of jet with a concentration C emitted from the bottom of the numerical tank. Second, the goal is to perform to the nonlinear resolutions of the Navier-Stokes equations in the case of high recirculation of the velocity field inside a large geometric area using a direct numerical simulation (DNS). The numerical resolution of the equations describing the phenomenon is made by the Hermitian finite differences method (HFDM) with higher accuracy [30,31,32,33]. Our choice is based on the high resolution of this method for the discretization of evolution equations and has better computational performance than classical schemes. Finally, we assess the capability and limitations of the numerical method used in this work to treat the jet flows.

This manuscript is organized as follows. After discussing the goal behind our research in the first section, Section 2 is dedicated to the position of the problem. Afterward, the numerical method used to resolve the Navier–Stokes equations is clearly described in the third section. Next, we validate the numerical method used to resolve the equations governing the problem of this work. Shortly after, the ambition is to characterize the instabilities when a jet of concentration C is emitted from the bottom of the numerical tank and to act on the nonlinear resolutions of the Navier–Stokes equation in a large computational domain. Before finishing, an evaluation of mesh domain sensitivity is presented to confirm the validity of the present numerical scheme, by comparing the numerical vertical velocity profiles for different steps of y-direction (flow direction) with the analytical solution. Finally, the paper is finished by conclusions and perspectives.

2. Problem Formulation

In the present study, the analysis is focused on studying the injection of homogeneous jet emitted from a bottom of numerical tank by considering an unsteady, incompressible, and viscous fluid in a numerical tank. The schematic diagram of the numerical tank setup is shown in Figure 1. The numerical tank consists of a rectangular domain with length L = 7 m and height H = 5 m. The jet is emitted via an orifice posed at the central zone of a flat bottom having the width e = 1.62 m, the injection is assumed uniform in the y-direction from an unsteady source of concentration C. The Reynolds number $Re=\frac{e{U}_m}{\upsilon }$ which is the key parameter of this study is calculated using the exit velocity ($U_m$) at the nozzle, width of nozzle (e), and kinetic $\upsilon$ viscosity of the fluid motion.

The formulation of this problem is described by the Navier–Stokes equations along with the mass transport equation, in terms of vorticity function $\omega$, stream function $\psi$, and concentration C, these equations are written in the dimensionless form as:

$$\frac{\partial \omega }{\partial t}+\left(\overrightarrow{U}.\overrightarrow{\nabla }\right)\omega =\frac{1}{Re}∆\omega +\frac{1}{F{r}^2}\left[\overrightarrow{\nabla }C\wedge \frac{\overrightarrow{g}}{\parallel \overrightarrow{g}\parallel }\right].\overrightarrow{k}$$

(1)

$$∆\psi =-\omega$$

(2)

$$\overrightarrow{U}=\overrightarrow{\nabla }\wedge \left(-\psi \overrightarrow{k}\right)$$

(3)

$$\frac{\partial C}{\partial t}+\left(\overrightarrow{U}.\overrightarrow{\nabla }\right)C=\frac{1}{ReSc}∆C$$

(4)

To achieve a complete formulation of the problem, it is necessary to set the initial and boundary conditions that are defined as:

-

At the initial instant, all variables are zero except on the diffuser where C = 1;

-

The conditions on lateral boundaries ${\Gamma }_U$ and ${\Gamma }_D$ are expressed as:

$$\frac{\partial C}{\partial n}=0,\frac{\partial \left(\overrightarrow{U}\cdot \overrightarrow{n}\right)}{\partial n}=0,\overrightarrow{U}\cdot \overrightarrow{\tau }=0,\frac{\partial \omega }{\partial n}=0$$

-

The free surface ${\Gamma }_F$ is assumed to be flat, and therefore the condition is written as:

$$\frac{\partial C}{\partial n}=0;\frac{\partial \left(\overrightarrow{U}\cdot \overrightarrow{\tau }\right)}{\partial n}=0,\overrightarrow{U}\cdot \overrightarrow{n}=0,\omega =0$$

-

On the bottom ${\Gamma }_B$ of the numerical tank, two conditions are defined:

• At the nozzle ${\Gamma }_E$, the conditions are expressed as:

  $$\overrightarrow{U}\cdot \overrightarrow{n}=1;C=1.\omega =0$$
• Aside from the nozzle ${\Gamma }_B$, the conditions are expressed as:

  $$\frac{\partial C}{\partial n}=0,\overrightarrow{U}=\overrightarrow{0},\frac{{\partial }^2\psi }{\partial n^2}=-\omega$$

3. Solution Procedures

In this section, we present the numerical method used to resolve the Navier–Stokes equations and the transport equation. Our choice fell on the compact Hermitian method $O\left(h^2,h^4\right)$ to resolve the equations describing the flow considered in this work [30,31,32,33]. Consequently, the $O\left(h^4\right)$ scheme combined with a scheme of implicit alternate direction methods (ADI) is used to resolve the Poisson equation of the stream function, and $O\left(h^2\right)$ scheme combined with a scheme of implicit alternate direction methods (ADI) is used to resolve the transport equations.

First, to use the $O\left(h^4\right)$ scheme combined with the scheme of implicit alternate direction methods (ADI) to resolve the Poisson equation, it is necessary to transform the elliptic equation of the stream function into a parabolic equation by introducing a fictitious time $\tau$. For that, a term $\frac{\partial \psi }{\partial \tau }$ is added to the second member of the Poisson equation. Hence, the equation rewrites as:

$$\frac{{\partial }^2\psi }{\partial x^2}+\frac{{\partial }^2\psi }{\partial y^2}=-\omega +\frac{\partial \psi }{\partial \tau }$$

(5)

considering $\tau =k∆\tau$, and designating ${\psi }_{i,j}^{n+\frac{1}{2}}$ and ${\psi }_{i,j}^{n+1}$ the values of the stream function $\psi$ at each fictitious half step, where i and j denote respectively the x and y-direction indices, and n is the time step. The discrete scheme of Equation (5) is treated in two half steps: the first step in x-direction and the second step in y-direction

$${\left(\frac{{\partial }^2\psi }{\partial x^2}\right)}_{i,j}^{n+\frac{1}{2}}+{\left(\frac{{\partial }^2\psi }{\partial y^2}\right)}_{i,j}^n=-{\omega }_{i,j}+\frac{1}{∆{\tau }_x}\left({\psi }_{i,j}^{n+\frac{1}{2}}-{\psi }_{i,j}^n\right)$$

(6)

$${\left(\frac{{\partial }^2\psi }{\partial x^2}\right)}_{i,j}{}^{n+\frac{1}{2}}+{\left(\frac{{\partial }^2\psi }{\partial y^2}\right)}_{i,j}{}^{n+1}=-{\omega }_{i,j}+\frac{1}{∆{\tau }_y}\left({\psi }_{i,j}^{n+1}-{\psi }_{i,j}^{n+\frac{1}{2}}\right)$$

(7)

at the first step following the x-direction, the value of ${\left(\frac{{\partial }^2\psi }{\partial y^2}\right)}_{i,j}^n$ in Equation (8) is assumed to be known, and hence we can write:

$${\left(\frac{{\partial }^2\psi }{\partial x^2}\right)}_{i,j}^{n+\frac{1}{2}}=S_i+\frac{1}{∆{\tau }_x}{\psi }_{i,j}^{n+\frac{1}{2}}$$

(8)

where

$$S_i=-{\omega }_{i,j}-{\left(\frac{{\partial }^2\psi }{\partial y^2}\right)}_{i,j}^n-\frac{1}{∆{\tau }_x}{\psi }_{i,j}^n$$

(9)

Taking account, the Hermitian formula [30,31,32,33] following the x-direction, that is written as:

$${\left(\frac{{\partial }^2\psi }{\partial x^2}\right)}_{i+1,j}+10{\left(\frac{{\partial }^2\psi }{\partial x^2}\right)}_{i,j}+{\left(\frac{{\partial }^2\psi }{\partial x^2}\right)}_{i-1,j}=\frac{12}{{\left(∆x\right)}^2}\left({\psi }_{i+1,j}-2{\psi }_{i,j}+{\psi }_{i-1,j}\right)$$

(10)

and the Equations (11)–(13) that are deduced by induction from the Equation (8)

$${\left(\frac{{\partial }^2\psi }{\partial x^2}\right)}_{i-1,j}^{n+\frac{1}{2}}=S_{i-1}+\frac{1}{∆{\tau }_x}{\psi }_{i-1,j}^{n+\frac{1}{2}}$$

(11)

$$10{\left(\frac{{\partial }^2\psi }{\partial x^2}\right)}_{i,j}^{n+\frac{1}{2}}=10S_i+\frac{10}{∆{\tau }_x}{\psi }_{i,j}^{n+\frac{1}{2}}$$

(12)

$${\left(\frac{{\partial }^2\psi }{\partial x^2}\right)}_{i+1,j}^{n+\frac{1}{2}}=S_{i+1}+\frac{1}{∆{\tau }_x}{\psi }_{i+1,j}^{n+\frac{1}{2}}$$

(13)

We can simply write the discrete set of Poisson equation at the first half-step, following the x-direction is written as:

$$\left(\frac{12}{{\left(∆x\right)}^2}-\frac{1}{∆{\tau }_x}\right){\psi }_{i-1,j}^{n+\frac{1}{2}}-\left(\frac{24}{{\left(∆x\right)}^2}-\frac{10}{∆{\tau }_x}\right){\psi }_{i,j}^{n+\frac{1}{2}}+\left(\frac{12}{{\left(∆x\right)}^2}-\frac{1}{∆{\tau }_x}\right){\psi }_{i+1,j}^{n+\frac{1}{2}}=S_{i-1}+10S_i+S_{i+1}$$

(14)

Thereafter, the x-direction allowed us to know ${\psi }_{i,j}^{n+\frac{1}{2}}$ and ${\left(\frac{{\partial }^2\psi }{\partial x^2}\right)}_{i,j}^{n+\frac{1}{2}}$ for the use case in the second half-step following the y-direction. Hence, from the Equation (7) we can deduce the formula:

$${\left(\frac{{\partial }^2\psi }{\partial y^2}\right)}_{i,j}^{n+1}=S_j+\frac{1}{∆{\tau }_y}{\psi }_{i,j}^{n+1}$$

(15)

$$S_j=-{\omega }_{i,j}-{\left(\frac{{\partial }^2\psi }{\partial x^2}\right)}_{i,j}^{n+\frac{1}{2}}-\frac{1}{∆{\tau }_y}{\psi }_{i,j}^{n+\frac{1}{2}}$$

(16)

Taking account, the Hermitian formula [30,31,32,33] following the y-direction, that is written as:

$${\left(\frac{{\partial }^2\psi }{\partial y^2}\right)}_{i,j+1}+10{\left(\frac{{\partial }^2\psi }{\partial y^2}\right)}_{i,j}+{\left(\frac{{\partial }^2\psi }{\partial y^2}\right)}_{i,j-1}=\frac{12}{{\left(∆y\right)}^2}\left({\psi }_{i,j+1}-1{\psi }_{i,j}+{\psi }_{i,j-1}\right)$$

(17)

and the Equations (18)–(20) that are deduced by induction from the Equation (15):

$${\left(\frac{{\partial }^2\psi }{\partial y^2}\right)}_{i,j+1}^{n+1}=S_{j+1}+\frac{1}{∆{\tau }_y}{\psi }_{i,j+1}^{n+1}$$

(18)

$$10{\left(\frac{{\partial }^2\psi }{\partial y^2}\right)}_{i,j}^{n+1}=10S_j+\frac{10}{∆{\tau }_y}{\psi }_{i,j}^{n+1}$$

(19)

$${\left(\frac{{\partial }^2\psi }{\partial y^2}\right)}_{i,j-1}^{n+1}=S_{j-1}+\frac{1}{∆{\tau }_y}{\psi }_{i,j-1}^{n+1}$$

(20)

Simply, the discrete set of Poisson equation at the second half-step, following the y-direction is written as:

$$\left(\frac{12}{{\left(∆y\right)}^2}-\frac{1}{∆{\tau }_y}\right){\psi }_{i,j-1}^{n+1}-\left(\frac{24}{{\left(∆y\right)}^2}-\frac{10}{∆{\tau }_y}\right){\psi }_{i,j}^{n+1}+\left(\frac{12}{{\left(∆y\right)}^2}-\frac{1}{∆{\tau }_y}\right){\psi }_{i,j+1}^{n+1}=S_{j-1}+10S_j+S_{j+1}$$

(21)

Moreover, based on the first-order derivatives of the Hermitian formula [30,31,32,33], the calculation of $\psi$ allowed us to calculate the components u and v of the velocity, for the sake of clarity:

-

The definition $\frac{\partial \psi }{\partial x}=v$, is used to calculate the component v using the formula

$${\left(\frac{\partial \psi }{\partial x}\right)}_{i+1,j}^n+4{\left(\frac{\partial \psi }{\partial x}\right)}_{i,j}^n+{\left(\frac{\partial \psi }{\partial x}\right)}_{i-1,j}^n=\frac{3}{∆x}\left({\psi }_{i+1,j}^n-{\psi }_{i-1,j}^n\right)$$

(22)

-

The definition $\frac{\partial \psi }{\partial y}=-u$, is used to calculate the component u using the formula

$${\left(\frac{\partial \psi }{\partial y}\right)}_{i,j+1}^n+4{\left(\frac{\partial \psi }{\partial y}\right)}_{i,j}^n+{\left(\frac{\partial \psi }{\partial y}\right)}_{i,j-1}^n=\frac{3}{∆y}\left({\psi }_{i,j+1}^n-{\psi }_{i,j-1}^n\right)$$

(23)

Furthermore, to resolve the transport equations in C and ω, we have chosen the $o\left(h^2\right)$ scheme combined with the implicit alternate direction method (ADI), the precision at the order 2 is applied to space variables and the order 1 to time variables. For more details, the transport Equations (1) and (4) are expressed in $\phi$ variable that is written in a conservative form as:

$$\frac{\partial \phi }{\partial t}+\frac{\partial \left(u\phi \right)}{\partial x}+\frac{\partial \left(v\phi \right)}{\partial y}=K\left(\frac{{\partial }^2\phi }{\partial x^2}+\frac{{\partial }^2\phi }{\partial y^2}\right)+S$$

(24)

where $\phi \equiv C$ or $\omega$:

-

in the case of $\phi =\omega$, this satisfies $S=\frac{1}{F{r}^2}\left[\overrightarrow{\nabla }C\wedge \frac{\overrightarrow{g}}{g}\right].\overrightarrow{k}$, and $K=\frac{1}{Re}$,

-

in the case of $\phi =C$, this satisfies S=0 and $K=\frac{1}{ReSc}$.

Let:

• ${\delta }_x$ denotes the discretized operator corresponding to $\frac{\partial }{\partial x}$
• ${\delta }_y$ denotes the discretized operator corresponding to $\frac{\partial }{\partial y}$
• ${\delta }_x^2$ denotes the discretized operator corresponding to $\frac{{\partial }^2}{\partial x^2}$
• ${\delta }_y^2$ denotes the discretized operator corresponding to $\frac{{\partial }^2}{\partial y^2}$

The variable $\phi$ is calculated in two half steps, we denote by ${\phi }^{n+\frac{1}{2}}$ the value of the variable calculated at the half step, and ${\phi }^{n+1}$ the value of the variable in the second half step. Based on ADI scheme, the discrete set of Equations (24) is written as:

• at the 1st half step:

  $$\frac{{\phi }^{n+\frac{1}{2}}-{\phi }^n}{∆t/2}-{\delta }_x\left(u^n{\phi }^{n+\frac{1}{2}}\right)+{\delta }_y\left(v^n{\phi }^n\right)=K\left({\delta }_x^2\left({\phi }^{n+\frac{1}{2}}\right)+{\delta }_y^2\left({\phi }^n\right)\right)+S$$

  (25)
• at the 2nd half step:

  $$\frac{{\phi }^{n+1}-{\phi }^{n+\frac{1}{2}}}{∆t/2}-{\delta }_x\left(u^n{\phi }^{n+\frac{1}{2}}\right)+{\delta }_y\left(v^n{\phi }^{n+1}\right)=K\left({\delta }_x^2\left({\phi }^{n+\frac{1}{2}}\right)+{\delta }_y^2\left({\phi }^{n+1}\right)\right)+S$$

  (26)

For greater clarity, the diffusion terms are discretized using central differences scheme, and the convection terms are treated by forward-backward differences method according to u and v velocities sign, this procedure ensures the preponderance of the main diagonal of the resulting matrices, for that matter we define

$$v_{i,j=}{\left(\frac{\partial \psi }{\partial x}\right)}_{i,j}^n{\alpha }_{i,j}=1\mathrm{if}{v}_{i,j}>0\mathrm{and}0\mathrm{elsewhere}$$

(27)

$$u_{i,j}=-{\left(\frac{\partial \psi }{\partial y}\right)}_{i,j}^n{\beta }_{i,j}=1\mathrm{if}{u}_{i,j}<0\mathrm{and}0\mathrm{elsewhere}$$

(28)

The discretization of the convection and diffusion operators are written as:

• at the 1st half step:

  $${\delta }_x\left(u^n{\phi }^{n+\frac{1}{2}}\right)=\frac{{\beta }_{i,j}\left(u_{i,j}{\phi }_{i,j}^{n+\frac{1}{2}}-u_{i-1,j}{\phi }_{i-1,j}^{n+\frac{1}{2}}\right)}{∆x}+\frac{\left(1-{\beta }_{i,j}\right)\left(u_{i+1,j}{\phi }_{i+1,j}^{n+\frac{1}{2}}-u_{i,j}{\phi }_{i,j}^{n+\frac{1}{2}}\right)}{∆x}$$

  (29)

  $${\delta }_y\left(v^n{\phi }^n\right)=\frac{\left(1-{\alpha }_{i,j}\right)\left(v_{i,j}{\phi }_{i,j}^n-v_{i,j-1}{\phi }_{i,j-1}^n\right)}{∆y}+\frac{{\alpha }_{i,j}\left(v_{i,j+1}{\phi }_{i,j+1}^n-v_{i,j}{\phi }_{i,j}^n\right)}{∆y}$$

  (30)

  $${\delta }_x^2\left({\phi }^{n+\frac{1}{2}}\right)=\frac{{\phi }_{i+1,j}^{n+\frac{1}{2}}-2{\phi }_{i,j}^{n+\frac{1}{2}}+{\phi }_{i-1,j}^{n+\frac{1}{2}}}{∆x^2}$$

  (31)

  $${\delta }_y^2\left({\phi }^n\right)=\frac{{\phi }_{i,j+1}^n-2{\phi }_{i,j}^n+{\phi }_{i,j-1}^n}{∆y^2}$$

  (32)
• at the 2nd half step:

  $${\delta }_x\left(u^n{\phi }^{n+\frac{1}{2}}\right)=\frac{\left(1-{\beta }_{i,j}\right)\left(u_{i,j}{\phi }_{i,j}^{n+\frac{1}{2}}-u_{i-1,j}{\phi }_{i-1,j}^{n+\frac{1}{2}}\right)}{∆x}+\frac{{\beta }_{i,j}\left(u_{i+1,j}{\phi }_{i+1,j}^{n+\frac{1}{2}}-u_{i,j}{\phi }_{i,j}^{n+\frac{1}{2}}\right)}{∆x}$$

  (33)

  $${\delta }_y\left(v^n{\phi }^{n+1}\right)=\frac{{\alpha }_{i,j}\left(v_{i,j}{\phi }_{i,j}^{n+1}-v_{i,j-1}{\phi }_{i-1,j}^{n+1}\right)}{∆y}+\frac{\left(1-{\alpha }_{i,j}\right)\left(v_{i,j+1}{\phi }_{i,j+1}^{n+1}-v_{i,j}{\phi }_{i,j}^{n+1}\right)}{∆y}$$

  (34)

  $${\delta }_x^2\left({\phi }^{n+\frac{1}{2}}\right)=\frac{{\phi }_{i+1,j}^{n+\frac{1}{2}}-2{\phi }_{i,j}^{n+\frac{1}{2}}+{\phi }_{i-1,j}^{n+\frac{1}{2}}}{∆x^2}$$

  (35)

  $${\delta }_y^2\left({\phi }^{n+1}\right)=\frac{{\phi }_{i,j+1}^{n+1}-2{\phi }_{i,j}^{n+1}+{\phi }_{i,j-1}^{n+1}}{∆y^2}$$

  (36)

Furthermore, the source term poses difficulties associated with gradient term, for that reason, the source term is treated by the expression giving the concentration using Green’s function, that is written as:

$$C\left(\overrightarrow{x},t\right)={{\displaystyle \iint }}_{\mathit{sup}\left(C\right)}\frac{Re}{4\pi t}C\left[{\overrightarrow{x}}^{\prime },t\right]\mathit{exp}\left(-\frac{Re}{4t}{\left(\overrightarrow{x}-{\overrightarrow{x}}^{\prime }\right)}^2\right)d\overrightarrow{x}\prime$$

(37)

where:

$$\left[\overrightarrow{\nabla }C\wedge \frac{\overrightarrow{g}}{g}\right].\overrightarrow{k}=\frac{\partial C}{\partial x}\left(\overrightarrow{x},t\right))$$

(38)

$$\overrightarrow{x}=x\overrightarrow{i}+y\overrightarrow{j}$$

(39)

After derivation of the expression giving the concentration using Green’s function

$$\begin{array}{c}\frac{\partial C\left(\overrightarrow{x},t\right)}{\partial x}=\frac{\partial }{\partial x}{{\displaystyle \iint }}_{\mathit{sup}\left(C\right)}\frac{Re}{4\pi t}{x}^{\prime }C\left[{\overrightarrow{x}}^{\prime },t\right]\mathit{exp}\left(-\frac{Re}{4t}{\left(\overrightarrow{x}-{\overrightarrow{x}}^{\prime }\right)}^2\right)d{\overrightarrow{x}}^{\prime }\\ ={{\displaystyle \iint }}_{\mathit{sup}\left(C\right)}\frac{Re}{4\pi t}C\left[{\overrightarrow{x}}^{\prime },t\right]\left(-\frac{Re}{2t}\left(x-x^{\prime }\right)\right)\mathit{exp}\left(-\frac{Re}{4t}{\left(\overrightarrow{x}-{\overrightarrow{x}}^{\prime }\right)}^2\right)d{\overrightarrow{x}}^{\prime }\\ =-\frac{Rex}{2t}{{\displaystyle \iint }}_{\mathit{sup}\left(C\right)}\frac{Re}{4\pi t}C\left[{\overrightarrow{x}}^{\prime },t\right]\mathit{exp}\left(-\frac{Re}{4t}{\left(\overrightarrow{x}-{\overrightarrow{x}}^{\prime }\right)}^2\right)d{\overrightarrow{x}}^{\prime }\\ +\frac{Re}{2t}{{\displaystyle \iint }}_{\mathit{sup}\left(C\right)}C\left[{\overrightarrow{x}}^{\prime },t\right]\left(Re\frac{x^{\prime }}{4\pi t}\right)\mathit{exp}\left(-\frac{Re}{4t}{\left(\overrightarrow{x}-{\overrightarrow{x}}^{\prime }\right)}^2\right)d{\overrightarrow{x}}^{\prime }\end{array}$$

(40)

Namely, we write:

$$\frac{\partial C}{\partial x}=-\frac{Re}{2t}xC\left(\overrightarrow{x},t\right)+\frac{Re}{2t}{{\displaystyle \iint }}_{\sup \left(C\right)}\frac{Re}{4\pi t}{x}^{\prime }C\left[{\overrightarrow{x}}^{\prime },t\right]\mathit{exp}\left(-\frac{Re}{4t}{\left(\overrightarrow{x}-{\overrightarrow{x}}^{\prime }\right)}^2\right)d{\overrightarrow{x}}^{\prime }$$

(41)

The passage to the discretized form is done by supposing that the particles cover the entire field, reads than for a particle i:

$${\left(\frac{\partial C}{\partial x}\right)}_i={{\displaystyle \iint }}_{P_i}\frac{\partial C}{\partial x}ds\approx -\frac{x_{i}Re\left|P_i\right|}{2t}{C}_i+\frac{Re}{2t}\left|P_i\right|{{\displaystyle \sum }}_{j=1}^{N_p}\frac{Re{x}_j{C}_j}{4\pi t}\mathit{exp}\left(-\frac{Re}{4t}{\left(\overrightarrow{x_i}-\overrightarrow{x_j}\right)}^2\right)$$

(42)

where Supp(C) stands for the support of the concentration function, $\left|P_i\right|$ is the surface of the particle with i index, $\overrightarrow{x_i}$ is the position of the particle, $\overrightarrow{x_j}$ and support particles, and $N_p$ is the total points number.

For more clarity, the Poisson equation and transport equations are constructed as linear systems, the associated matrices are written as $\left[A\right]\left\{\psi \right\}=\left\{B\right\}$, and $\left[C\right]\left\{\phi \right\}=\left\{D\right\}$, these algebraic systems are resolved by the application of the successive over-relaxation (SOR) iterative method [34] taking into account the convergence criterion

$$\left|{\tilde{\psi }}^{k+1}-{\tilde{\psi }}^k\right|\le {10}^{-4}$$

(43)

$$\left|{\tilde{\phi }}^{k+1}-{\tilde{\phi }}^k\right|\le {10}^{-4}$$

(44)

where ${\tilde{\psi }}^k$ and ${\tilde{\phi }}^k$ represent unknown vectors at $k^{th}$ iteration associated with (SOR) method.

4. Results and Discussion

4.1. Validation Test

The objective of this subsection is to validate the numerical method used in this work. For that, we consider a viscous, unsteady, incompressible fluid in a lid-driven cavity. For the sake of clarity, all walls are considered with no-slip or impervious boundary conditions, except the velocity at the upper wall set to unity (see Figure 2). To verify the validity of our results, we compare the profiles of the velocity components with the reference [35] at the median of the cavity, for different Reynolds numbers: 100, 400, and 1000.

The Figure 3 and Figure 4, outline the results respectively of the velocity profiles: u-components as a function of y by fixing the x-axis in (x = 0.5), and v-components as a function of x by fixing the y-axis in (y = 0.5) for the Reynolds numbers Re = 100, 400, 100, and 2000.

Based on the visualization of the velocity profiles presented in Figure 3 and Figure 4, we confirm that the results have shown a good agreement with the results obtained in the literature [36], which confirm the validation of the accuracy of the numerical code used in this work.

4.2. Jet Emitted from the Bottom of the Numerical Tank

The intention in this subsection is to analyze the dynamic properties of jet type flow, and to assess on the limitation of the compact Hermitian finite differences method (FDM) higher accuracy, to resolve the problem of jet with a concentration C emitted from a bottom of numerical tank (see Figure 1), For that, the Equations (1)–(4) are resolved using the Hermitian finite differences method in a dimensionless rectangular domain sized $\mathsf{\Omega }=$ (0, 7) × (0, 5).

First, in Figure 5 and Figure 6 we present the velocity fields, and the iso-contours of concentration for the different instant of time, by fixing the dimensionless numbers: Re = 100, Fr = 10, and Sc = 0.7, these results are presented for two reasons; in the first place, to assess the capability of the direct numerical simulation (DNS) to describe the evolution of jet flow in terms of concentration and velocity fields for low Reynolds number. In the second place, the objective is to analyze the validity of the proposed numerical approach and also for refining the description of the flow, particularly in the vicinity of the diffuser as the significance of this problem for many industrial applications.

Second, the goal is to assess the limitation and the response of the numerical code used in this work for high Reynolds number. For that, we present in the Figure 7 and Figure 8 the results respectively obtained by the stream function, and iso-values of vortex function by fixing the time at t = 1s, and varying Reynolds numbers from Re = 100 to Re = 1000.

The results presented as a form of velocity fields in the Figure 5, show the appearance of symmetrical vortex cells located on either side of the axis of jet for different instant of time. Further, the Figure 6 shows also symmetrical profiles of iso-contours of concentration function according to the axis of the jet over the time. Next, the results of the Figure 7 and Figure 8 show a dissymmetry of profiles of streamlines, and iso-vortex with respect to the axis of the jet at Reynolds number Re = 1000. In addition, our numerical experience show that flow is weak function of Reynolds number, when Reynolds number is higher than 1000. Further, the present numerical results have efficiently provided information about the associated eddy structures in the near-wall or impact region. Hence, it is quite possible to point out the influence of Reynolds number on eddy structures. In addition, the kinematic description of the impact region of the jet has not been the subject of many studies although significant transfers (heat/mass) occur only in this region. It is characterized by a significant vertical diffusion of momentum for the mean flow.

Furthermore, to check the accuracy associated with Reynolds number, it is worthy to present the mean absolute error (MAE) for different Reynolds numbers, these errors are given by the following formula:

$${\mathrm{MAE}}_{\mathsf{\omega }}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}\left|{\mathsf{\omega }}_{\mathrm{i}+1}-{\mathsf{\omega }}_{\mathrm{i}}\right|}{\mathrm{N}}$$

(45)

$${\mathrm{MAE}}_{\mathrm{C}}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}\left|{\mathrm{C}}_{\mathrm{i}+1}-{\mathrm{C}}_{\mathrm{i}}\right|}{\mathrm{N}}$$

(46)

$${\mathrm{MAE}}_{\psi }=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}\left|{\psi }_{\mathrm{i}+1}-{\psi }_{\mathrm{i}}\right|}{\mathrm{N}}$$

(47)

In the

Table 1. Mean absolute error (MAE) for different Reynolds numbers.

Mean Absolute Error (MAE)	Reynolds Numbers
	100	400	1000
${\mathrm{MAE}}_{\omega }$	4 × 10−5	3 × 10−4	5 × 10−3
${\mathrm{MAE}}_C$	5 × 10−5	2 × 10−4	7 × 10−3
${\mathrm{MAE}}_{\psi }$	2 × 10−5	4 × 10−4	9 × 10−3

we present the mean absolute error (MAE), for vortex function $\omega$, concentration function C, and stream function $\psi$, for Reynolds numbers $Re=100,400\mathrm{and}Re=1000$.

The results illustrated in the Table 1 that represents the mean absolute error for stream function, vortex function, and concentration function, show the onset of the numerical instabilities at Re = 1000, these instabilities come from the ill-conditioning effects of the algebraic systems when Reynolds number increases. For further clarity, the presented results show that the Hermitian finite differences method is suited to the problem of jet up to Reynolds number Re = 1000. Otherwise, in the case of Reynolds number higher than Re = 1000, the method appears inappropriate. Furthermore, the DNS using the Hermitian finite differences method allows us to find satisfactory results and a good representation of the jet flow but requires a refined mesh.

4.3. Evaluation of Mesh Domain Sensitivity

To ensure the validity of the present numerical scheme, in Figure 9 we provide a valuable insight into the grid domain sensitivity of the jet flow, by presenting the vertical velocity profiles for different steps of y-direction (direction of the jet flow), and compare them with the analytical solution of Bradbury (1965), these comparisons will be presented for different Reynolds numbers. To assess the grid sensitivity and evaluating mesh dependence to the accuracy of the jet flow, we use different two-dimensional grids mesh $\left(N_x\times N_y\right)$, and adopting a systematical increase of mesh points in y-direction, the step in y-direction is defined as $∆y=\frac{H}{N_y-1}$, where H is the height of computational domain. The different finite grids used in this validation test are 26 × 61. 26 × 75, and 26 × 90, leading respectively to different steps in y-direction $∆y_1=0.0847$, $∆y_2=0.0676$, and $∆y_3=0.051$.

The numerical simulations of vertical velocity profiles, using the different steps of y-direction $∆y_1$, $∆y_2,$ and $∆y_3$, have shown a good agreement with the analytical solution for different Reynolds numbers, which confirm the accuracy of the numerical procedure used in this work. Further, the proposed study shows that the center line velocity is independent of Reynolds number as pointed in the references [36,37].

5. Conclusions

In this paper, a numerical approach was proposed to resolve equations that have strong nonlinearities because of the high recirculation of the velocity in the vicinity of the boundary layers. The fluid flow is described by the Navier–Stokes equations that are resolved by finite difference method (FDM) high accuracy. First, the validation tests are ensured by comparing the profiles of the velocity components at the median of the lid-driven cavity with previous works of the literature, and assessing the mesh domain sensitivity by comparing the numerical velocity profiles for different steps of y-direction (flow direction) with the analytical solution. The presented numerical profiles have proven a good agreement with reference results, which confirm the validity of the numerical method used in this work.

Second, the objective was to revolve around discussing the behavior and the dynamic properties of jet type flow. The phenomena studied locally in terms of velocity fields and concentration made it possible to realize the importance of the interaction and the evolution of the flow over the time. Further, the instabilities mechanism of the flows associated with high Reynolds number are detected by the numerical simulations. Furthermore, finite differences method showed that it is well suited to the problems of jet flows at moderate Reynolds number. Moreover, the results of this work prove that the velocity profiles does not show any variation with the Reynolds number which confirms that the development of the jet is not affected by the Reynolds number.

To sum up, this research describes the limitation of the DNS using the HFDM to treat the problem of jet emitted from a bottom of numerical tank, the results presented in this research have made it possible to find satisfactory results, but there are difficulties related to the numerical instabilities associated with high Reynolds numbers. As a perspective, we endeavor to investigate the limitation of the Hermitian finite differences method to treat the 3D jet flow simulation. Besides, to inspect the capability of spectral methods [38], and meshless methods [39,40,41,42,43] to deal with the jet problem at higher Reynolds number.