Turbulent Flow Simulations

Abstract

This chapter consists of four sections; Introduction, Boundary and initial conditions, Setting of nonuniformly spaced grid, and Simulation approach. The fields of the fluctuating velocity describe the turbulent flows. Such fluctuations blend the transported quantities, like species concentration, energy, and momentum and make the transported quantity fluctuations in addition. Because such fluctuations have a small scale and high frequency, so they are very computationally costly for simulation straightforward in the virtual engineering computations. Alternatively, the instantaneous accurate governing equations can be time-averaged, ensemble-averaged, or in different way handled for removing the small scales, causing a modified group of equations, which are computationally less costly for solving.

1. Introduction

The updated evolution in technology of computer has offered a fresh insight into the turbulent flow phenomenon. Particularly, a numerical analysis approach named direct numerical simulation (DNS), which is clearly depended upon the physical laws with no turbulence model, is anticipated to become too beneficial due to the Navier–Stokes equation reliability and versatility. Nevertheless, the DNS used to the turbulent flows with intricate shapes of boundary or large Reynolds (Re) nos. Puts a great challenge to the engineers, who are doing different efforts for overcoming the relevant problems. The grids that are staggered have traditionally been utilized for the computation, and the grids that are non-staggered are more often utilized in the recent time. For example, investigations on the organized grid systems have progressed quickly owing to their expansion to the generalized coordinate regimes. The use of DNS to the fields of the turbulent flow with the free surfaces, like seas and rivers, was initially tried by Lam and Banerjee in 1992. From that time, such method has supplied significant information on the turbulent structures and mechanism of the mass transport close to the free surfaces that are hard for investigating experimentally. The whole of such investigations have utilized either finite difference method or spectrum method utilizing a staggered grid, and there have been no try for using a normal grid for a method.

2. Boundary and initial conditions

2.1 Periodic boundary conditions

In the hydrodynamically and thermally fully (time) developed conditions as well as for the periodically fully developed conditions, where the flow uniformity in the stream-wise direction is possible to be assumed, the periodic boundary conditions (BCs) can be employed on the whole variables along the stream-wise direction, involving temperature and pressure. The channel flow is the highly often utilized geometries in the DNS, where the periodic BCs are entailed. The turbulent flow simulations in a channel having an infinite width (z) and length (x) are conducted via approximating the domain infinite extent in the span-wise and stream-wise directions with a finite length computational domain (named unit cell) and via using the periodic BCs in such directions. These circumstances theoretically impose the depending variable fields to become the periodic along the chosen direction. Practically, their use relies upon the specific approach itself: for example, the periodicity is inherent to the spectral approaches’ principle functions, whereas in the finite volume/finite difference methods, it is performed via entailing the computational stencils’ circular topology at the domain double terminals. A particular care must be considered when selecting the unit cell extent that must be highly sufficient for encompassing the flow biggest scales. When such states are met, the flow inside the computational domain could be regarded as accurate flow depicting in the infinite channel.

2.2 Inflow open boundary conditions

The periodic BCs could be used in the stream-wise direction to calculate variables if fully developed circumstances are satisfied. In case of the attention is a spatially evolving flow, the approach of DNS needs re-sorting to a method of the inflow turbulence creation for the fields of temperature and velocity. Recently, the methods of the inflow turbulence creation for the field of velocity have been categorized and surveyed by Wu in 2017. A short subsidiary periodic domain, where the turbulence is calculated in parallel to the principal computation, is giving outcome results utilized as inflow state for target simulation. Such kind of inflow (BCs) is recognized as a re-cycling approach. Fabricated approaches represent a substitute trial for fully modeling the inflow.

2.3 Outflow open boundary conditions

Contriving proper BCs at the open boundaries if the average flow exists, the domain can be a hard work. The highly familiar difficulties of the outflow open BCs are that they do not make genuine reflection of the velocity and of the transported scalars and they try to make wiggles at the boundary domain and might influence the solution algorithm stability. Then, the outflow BCs have to depict the physical domain continuation after calculating the domain edge, and the transport equations extrapolated from the domain interior are frequently enforced in this form:

$$\frac{\partial \varphi }{\partial t}+U_c\frac{\partial \varphi }{\partial x}=0$$ E1

In several situations, a consistent value for U_c that ensures the mass conservation is utilized. The interpretation for the convective velocity U_c can be entailed in Eq. (1) and depended upon the configuration specific flow. For instance, some researchers, who dealt with the turbulent jets and plumes and proposed using the exponential function for U_c (y). In the survey on the open outflow circumstances by the approach was detailed as a promising method for the buoyant, turbulent plumes. The approach suggested by some authors was depended upon a 1-D advection diffusion expression at the outflow, where the term of diffusion was added to Eq. (1) and a phase velocity was presented as well as U_c .

2.4 Boundary conditions for the thermal field

The geometry of turbulent channel is also highly often utilized for the investigations of the heat transfer of near-wall. The highly exact method is a conjugate heat transfer, where the heat conduction within the actual heated walls is considered. The related investigation for the uses of liquid metal is the DNS that conducted a conjugate heat transfer (DNS) at low Prandtl no. (Pr = 0.01). Majority of studies of (DNS) heat transfer were conducted with the simple thermal BCs with no solid walls, and dimensionless temperature at the contact plane of fluid–solid was set to zero:

$$\theta \left(y=\pm h\right)=0\left(\text{non}-\text{fluctuating thermal}\mathrm{BC}\right)$$ E2

Such kind of thermal BC is a too good approximation of actuality if the ratio of thermal activity $K=\sqrt{{\left(\mathit{\lambda \rho }{c}_p\right)}_f/{\left(\mathit{\lambda \rho }{c}_p\right)}_w}$ moves to zero, whereas the thickness of the heated wall and parameter ( $G={\alpha }_f/{\alpha }_w$ ) stay finite $\left(\alpha =\lambda /\rho c_p\right)$ . In such situation, the temperature fluctuations created in fluid are not penetrating into solid wall. Such kind of perfect thermal BC was represented as a non-fluctuating BC. A case study of this system is the solid and fluid air/metal combination. For the combination of water/steel, thermal BC is yet near to the non-fluctuating BC when the heated wall is sufficiently thick; nevertheless, this approximation does not succeed if the thickness of metal wall is small. The thermal BC that allows the peak penetration of the turbulent temperature fluctuations from fluid to solid is expressed as.

$${⟨\theta \left(y=\pm h\right)⟩}_{x,z,t}=0{\left(\frac{\partial {\theta }^{\prime }}{\partial y}\right)}_{y=\pm h}=0\left(\text{fluctuating thermal}\mathrm{BC}\right)$$ E3

The average dimensionless temperature for fluctuating temperature (BC) at wall is yet set to zero, whereas the zero derivative is suggested for the temperature ( ${\theta }^{\prime }$ ) fluctuating component. Such kind of thermal BC is performed if the ratio (K) of thermal activity moves to infinity; experiments can be repeated with a water flume, which is heated with a too thin metal foil. In 2001, when initial analysis of BCs (Eqs. 2 and 3) in a turbulent flow of channel was conducted by Tiselj et al. (2001b), such BCs were represented as “isothermal” BC (Eq. 2) and “isoflux BC for a dimensionless temperature” (Eq. 3). Particularly, the “isoflux” term results much confusion, since it was often blended with the standard “isoflux” BC for a physical temperature. Such pair of BCs must not be wrong: the aspect with a fixed volumetric heat origin within the solid walls matches to the “isoflux” BC for the physical temperature. The explained BC in such situation can be either non-fluctuating (Eq. 2) or fluctuating (Eq. 3). Hence, Tiselj and Cizelj [1] recommended new names for BCs (Eqs. 2 and 3). A significant characteristic of various thermal BCs (Eqs. 2 and 3) is that in the passive scalar approximation, they possess approximately negligible effect on the coefficient of heat transfer. In other words, the average temperature profiles close to the wall are specially the explained thermal (BC) independent. Certainly, an important discrepancy can be noted in the temperature fluctuations and the other turbulent statistics. For actual fluid–solid regimes, the thermal BC is common within the two bounded BCs of Eqs. (2) and (3). For the composition of steel and liquid sodium with (K) equals (1), the thermal (BC) is approximately in the middle within the non-fluctuating and fluctuating temperature (BC). Thus, so as to investigate the temperature fluctuation details with their penetration inside the solid wall, the method of a conjugate heat transfer is advised for the heat transfer (DNS) in the liquid metals.

2.5 The initial circumstances

Beyond the whole equations, computational domain and BCs have been identified; the initial circumstances are needed for beginning the simulations. Provided that the statistics for a long term are not dependent of the initial circumstances and given (CPU) time is not a bounding factor, it is frequently enough to begin with any type of primary velocity containing in minimum some perturbations of long enough wavelengths. Since such simulations’ initial stages are run on the highly rough grids, the consumption of CPU time should not become very high. In major situations, if the geometry stays alike, it is enough to utilize instantaneous fields determined from the simulations at the other Re or Pr nos. Beyond proper scaling. In state of computational expense is a bounding factor; thus the primary transient is shortened with a few more advanced kinds of the initial circumstances.

3. Setting of nonuniformly spaced grid

3.1 The grid system outline

The generation of computational grid can be divided into two types: a regular grid regime and a staggered grid regime in accordance with the evaluation points’ arrangement of pressure and velocity. The staggered types have been widely utilized in the spatial separation approaches of the analysis of flow either in the MAC or Simple approach since they govern the solution numerical vibration. Nevertheless, they need many interpolating processes, and calculating the loads is timorous, particularly if a generalized coordinate regime is employed. So, the grid regime choice is tightly related to the turbulent flow computation compactness [2]. On the other hand, the regular grids have been indicated the suffering from the taking place of solution spatial vibration. Nevertheless, they are suitable for the programming since the physical quantities are specified at similar point. The other advantages of these grids are the properness to improve the solution accuracy and the ease of conversion into a generalized coordinate regime that is influential in the deal with the intricate boundary difficulties. The regarded numerical vibration can be eliminated via the nonuniformly spaced grid adoption, and the active simulation can be recognized via utilizing a rough grid over the areas where the velocity spatial variations are small, with exception of the near-wall zone.

3.2 Influence of nonuniformly spaced grid

Because the nonuniformly spaced grid adoption, in general, results in lowered accuracy, here the influences of such kind on the numerical solution are debated from a 1-D linear convection problem referred by Eq. (4):

$$\begin{array}{c}\frac{\partial u}{\partial t}+U\frac{\partial u}{\partial x}=0,\left(U=1,-\infty <x<\infty \right)u\\ \left(x,0\right)=\left\{\begin{array}{l}\frac{1-cos\left(2x-\pi \right)}{2},\left(\frac{\pi }{2}\le x\le \frac{3\pi }{2}\right)\\ 0,\mathit{\text{otherwise}}\end{array}\right\}\end{array}$$ E4

As depicted in Figure 1, three kinds of grids are employed: a uniform-spaced grid, a smooth nonuniformly spaced grid (A) having a continuous metric ( $x_{\zeta }$ ), and a smooth nonuniformly spaced grid (B) having discontinuities in metric ( $x_{\zeta }$ ). Because the accurate solution of Eq. (4) is really similar to the primary waveform and does not vary with the time, the numerical solution beyond the (m-th) cycle is compared with the accurate solution via the periodic boundary condition adoption for the two terminals of x-axis. The scheme of the third-order Adams-Bashforth is employed for the time progress.

Figure 1.

Coordinate transformation.

The comparison of waveforms beyond 50 cycles of computation for a Courant no. of 0.1 in state of the uniform-spaced grid is presented in Figure 2. Graphs from (a) to (c) at the top depended upon a central difference scheme, while the graphs from (d) to (f) at the bottom side are based on an upwind difference scheme. This figure manifests that the nonuniformly spaced grid (A) causes somewhat a higher attenuation error than that for the uniform-spaced grid. That is likely resulted via the truth that the accuracy is reduced by one order if the grid spacing does not smoothly vary in the nonuniformly spaced grid, as it can be concluded from the second-order differential approximation truncation error which is inferred via the Taylor expansion. Graphs from (c) to (f) propose that such reduction of accuracy can be substituted via the higher-order finite difference approximation adoption. As for the nonuniformly spaced grid (B), not merely the error of attenuation but also the phase one is created in the whole situations. Since such a grid possesses similar smoothness as the nonuniformly spaced grid (A), the error of phase is regarded to be associated with the metric ( ${x_{\zeta }}^{\prime }$ S) discontinuities. Concerning the finite difference convection term scheme, the upwind difference is, in general, stated as a sum of the numerical viscosity and the central difference. So, the central differences in the graphs from (a) to (c) be the upwind differences in the graphs from (d) to (f) via the addition of truncation errors (numerical viscosity) of the second-, fourth-, and sixth-order differential, correspondingly. Therefore, the numerical viscosity influence can be studied via comparing the graphs at the top of Figure 2 with those at the bottom. While the numerical viscosity possesses a suppressing influence of the vibrating solution obtained in the central difference, the influence is highly deep in the lower-order upwind differences so as the solution is skewed. Even the third-order upwind difference that is frequently utilized in the practical calculation is not fitted for the analysis of turbulence in which too high accuracy is needed. The fifth-order upwind difference that entails a small reverse influence on the solution works as a filter with a high cut, which eliminates merely the components having a high wavenumber. In the central differences, even in state of high-order accuracy, the nonlinear instability is probably triggered via the dispersion which resulted via the truncation errors of the odd-numbered differential. So, it is regarded that the fifth-order upwind difference adoption is active in eliminating the aliasing errors that are resulted via the components having a high wavenumber created in a continuous way from the nonlinear terms and in the suppressing solutions, which include the numerical vibrations having no physical meaning.

Figure 2.

Waveform comparison for different types of grid spacing ( $c=0.1,x/2\pi =50$ ).

Therefore, in the current analysis of the study utilizing a regular grid in a generalized coordinate regime, a nonuniformly spaced grid with enough smoothness and a continuous metric is created. Also, the scheme of the fifth-order upwind difference is employed for the term of convection.

4. The approach of simulation

4.1 Numerical aspects

The DNS of the turbulent boundary layer)TBL(with the embedded instruments, like vortex generators, represents a step in intricacy [3], and it is likely with the NTS code [Owing to its combination of numerical dissipation, multi-block capability, and implicit numeric]. This code supplies the ability for the accurate computations with time with the structured multi-block overset grids and uses the scheme of an implicit second order in time flux-difference splitting. The flow is dealt as an incompressible. The inviscid fluxes in governing equations are approximated with the centered fourth-order accurate differences everywhere, with the exception of the devices’ stagnation points near vicinity where the third-order upwind approximations are utilized. Which is required for suppressing the odd–even pressure oscillations resulted via the vigorous gradients of pressure in such areas. The viscous fluxes are approximated with the centered second-order differences. At each step of time, the resulted finite-difference equations are solved with the employment of Gauss-Seidel relaxations via planes and sub-iterations in “pseudo-time.”

4.2 Turbulent inflow conditions

It is the more original aspect of the approach. It is preferred to possess a standard TBL that encounters the device, which controls the flow, without a long method zone. Some researchers revealed that the inflow conditions which depended on the random nos. Can simply take 20 $\delta$ for recovery; this would be too wasteful if the beneficial zone is of the order of 25 $\delta$ . Their method of recycling method, in some sense, wasted around 8 $\delta$ that is too smaller and could likely have been decreased more with no damage. Here, the method is vigorously inspired via LWS, but it is simpler, and a number of separate enhancements were done. The LWS used a various processing to the outer and inner zones of TBL and mixed the two expressions. That may seem optimum, but it needs re-scaling the velocity of friction as well as the thickness of BL and presents many arbitrary parameters and almost intricate expressions. This method failed to have the benefits of a pair of desirable facts. The first is that the near-wall turbulence re-creates itself too quickly than the outer-zone turbulence, and so a small damage is made via using the outer-zone scaling in every part. The second is that if the recycling plant is becoming too near the inflow that is favorable in respect to the cost of computing, the clashing between the outer-zone and the inner-zone scaling basically disappears. Also, a short recycling distance decreases the shifting in the span-wise length scales and the time scales, neither of which is taken into consideration (so, the recycling is better modified with shorter distances, till the distance is too short that the recycling and inflow planes are too near that no eddies development can occur). For such causes, the current method utilizes a single re-scaling. In a similar way, the corrections to the component (v) of the wall-normal velocity possess a too little influence and are deleted. Eventually, there is no subsidiary simulation here. The initial simulation yields its own recycling information.

The method conclude of utilizing the velocity field at a “recycling” plant (x = x_r ), which is partway within the domain, properly changed, to give the velocity field at the inflow (x = 0). Because the BL is somewhat thicker at (x_r ), the scaling in the wall-normal direction (y) is the minimum change, which is required. Let the thicknesses at (x = 0) and (x = x_r ) be ( ${\delta }_o$ ) and ( ${\delta }_r$ ). Thus, the inflow condition for the velocity vector (U) is:

$$U\left(0,y,z,t\right)=U\left(x_r,y{\delta }_r/{\delta }_0,z+\mathrm{\Delta }z,t-\mathrm{\Delta }t\right)$$ E5

where the span-wise shift, Δz, is presented to preserve the turbulence at the inlet and the recycling sections out of the phase (in the simulations debated below, the Δz value was put equal to half period of the span-wise). That is a small step, which is taken for disorganizing any lasting span-wise changes of the average flow that would otherwise be recycled and likely possess more time for damping via the span-wise diffusion. Utilizing the condition from the final step of (t – Δt) is highly suitable.

The procedure simplicity is visible in comparison with the LWS, but the fundamental concept is theirs. The generalization to flows with the gradients of pressure and the curved boundaries stays to be done. The operation is governed via a pair of two parameters; the difference of the (Re) no. from the inflow to recycling, $\mathrm{\Delta }{R}_x=x_r{U}_e/v$ , and ( ${\delta }_0/{\delta }_r$ ) ratio. By assuming the approximate scaling for the thickness of BL, the flow should be stabilized with the following magnitude for the thickness of inflow:

$$\delta =0.37\left(\frac{v}{U_e}\right).{\left[\frac{\mathrm{\Delta }{R}_s}{{\left({\delta }_r/{\delta }_0\right)}^{5/4}-1}\right]}^{4/5}$$ E6

This provides a fine control over the thickness of inflow (despite the factor (0.37) is not accurate, the little discrepancy in the rates of growth between the thickness of momentum s and the thickness of BL is ignored, and the turbulence is not the exact evolved condition), in which the response to the little variations in ( $\delta$ Rx) and ( ${\delta }_0/{\delta }_r$ ) can be predicated (as the flow is properly set). The primary conditions for this simulation have some of significance. For example, the random perturbations, which are very feeble or very insufficient in length scales, can too well “die out.” When it is, the simulation is becoming laminar and producing a higher thickness of BL than the anticipated with the turbulence, due to the various scaling than in the above equation, which will make the obvious failure. That is more possibly with a short distance of re-cycling. The ( $\delta$ R_x ) and ( ${\delta }_0/{\delta }_r$ ) little values costly affect the steady state but resist the turbulence maturation, which is primarily required. So, it can be useful to begin with the recycling plant farther downstream and then move it nearer to the inflow via decreasing ( $\delta$ R_x ) and ( ${\delta }_0/{\delta }_r$ ) in agreement. Instead of, when a device, which controls the flow, finishes the normal BL and bounds the recycling plant location, the simulation can begin with a further length that is attached to the upstream grid, which would later be removed beyond the turbulence is well established, for reducing the long simulation’s cost required for an adequate sample time. Moreover, the early part of the simulation does not require the complete length of the zone containing the recovery and device, particularly when the parameters of the recycling are being modified. In brief, the flexibility in re-beginning the simulations from each other with various domains will be too useful, and many approaches exist for surmounting the problems of the early system, based upon the preference of the user. As the procedure of recycling is properly setup, the solution loses the memory of the primary state processing. The method used here takes the advantages from another code, so as it is not a universal. It is as the following: a RANS solution is computed with the wanted thickness, and the perturbations is added to it. Such perturbations are determined from a code, which is utilized for initializing the uniform turbulence simulations. It gives turbulence volume of cubes shape with 0.063 size, each one has random phases that are positioned side-by-side upon the wall beyond their velocity field is multiplied via a shape function, so as they become zero at the wall and in the free stream. That acts properly, despite nothing is made for injecting the exact Reynolds shear stress. An appropriate solution is determined in some flow-through periods of times.