On stochastic inlet boundary condition for unsteady simulations

. The paper deals with the stochastic generation of synthesized turbulence, which may be used for a generating of an inlet boundary condition for unsteady simulations, e.g. Direct Numerical Simulation (DNS) or Large Eddy Simulation (LES). Assumptions for the generated turbulence are isotropy and homogeneity. The described method produces a stochastic turbulent velocity ﬁeld using the synthesis of a ﬁnite sum of random Fourier modes. The calculation of individual Fourier modes is based on known energy spectrum of turbulent ﬂow, and some turbulent quantities, e.g. turbulent kinetic energy and turbulent dissipation rate. A division of wave number range of the energy spectrum determines directly the number of Fourier modes, and has a direct impact on accuracy and speed of this calculation. Therefore, this work will examine the inﬂuence of the number of Fourier modes on a conservation of the ﬁrst and second statistical moments of turbulent velocity components, which are prespeciﬁed. It is important to ensure a su ﬃ cient size of a computational domain, and a su ﬃ cient number of cells for meaningful comparative results. Dimensionless parameters characterizing the resolution and size of the computational domain according to a turbulent length scale will be introduced for this purpose. Subsequently, the su ﬃ cient values of this parameters will be shown for individual numbers of Fourier modes.


Introduction
It is well known that the inlet boundary condition (BC) has a significant effect on solution, especially for unsteady simulations. Any input unphysicality may thus lead to larger or smaller error.
If we perform a steady-state simulation, such as Reynold Averaged Navier-Stokes (RANS) simulation, it may only be necessary to impose the mean velocity profile with mean turbulent quantities at the turbulent inflow. Despite the difficulties associated with setting of appropriate characteristics of turbulent flow, it is significantly easier than to generate turbulence itself, which is necessary for transient simulation such as LES and DNS.
As we mentioned before, the correct specification of the inlet BC is critical due to their strong influence on the solution. There exists many different approaches to create inflow turbulence. The first, very often used, approach is based on extracting of data directly from a precursor DNS of channel flow. This method is accurate if the Reynolds number (Re) of DNS is relevant. In another case, if the Re of the DNS differs from that which we need, it is very difficult to re-scale the DNS fluctuations. Therefore, this approach is applicable only for low Re due to the enormous computational demands of DNS simulation. The second approach to create inflow turbulence uses a mapping of resolved fluctuations at the plane downstream of an inlet plane, back on the inlet plane. Both of these two approaches are time-consuming, because we need firstly to perform the precursor DNS simulation, or we have to wait for sufficiently well developed turbulence due to the mapping. Moreover, the mapping area itself causes an enlargement of computational domain and thus it increases a coma e-mail: niedoba@fme.vutbr.cz putational effort. Hence, the idea of the method, which generates directly the turbulent inlet BC, evolves naturally. Otherwise, the turbulence is exactly governed by Navier-Stokes equations, and it seems to be a difficult task to generate physically realistic turbulent fluctuations independently of these equations. There is fortunately the third class of methods, the stochastic methods, which are able to preserve the spatial and temporal correlations together with the first and second statistical moments, and prescribed energy spectra.
The stochastic method described bellow is based on the approach devised by Kraichnan [1] and improved by Karweit et al. [2], where the spatially correlated turbulent velocity field is defined as a finite sum of discrete Fourier modes. This method was originally used for generating a turbulent velocity field from the turbulent quantities obtained by RANS simulation, and subsequent computing of acoustic source terms as a right hand side of an appropriate propagation equations of sound waves. The first formulation of this model, Stochastic Noise Generation and Radiation (SNGR) approach, was applied to subsonic jet noise calculation by Béchara et al. [3]. In a further development, Bailly et al. [4] introduced a time dependent term into the Fourier modes. A different way of introducing time dependency based on an asymmetric time filter was presented by Billson et al. [5].
A further development of this method for generating a turbulent inlet BC was proposed by Davidson [6].

Synthesized turbulence
The method for generating a homogeneous isotropic turbulence having the proper spatial characteristics was pro-posed in [1]. The turbulent velocity field is defined using the Fourier mode approach.
Consider a Fourier decomposition of a turbulent homogeneous isotropic field u = (u 1 , u 2 , u 3 ) at given point where κ is a wave vector and j is the unit imaginary number ( j 2 = −1). Assuming that u is a real function, and using the even property of cosine function, we may approximate the equation (1) as a finite sum of N discrete Fourier modes whereû n , ψ n , and σ n are the amplitude, phase, and direction, respectively, of the nth Fourier mode associated with the wave vector κ n . Each wave vector κ n is picked randomly on a sphere of radius specified by a wave number κ n = ||κ n || to ensure isotropy. Hence, the wave vector κ n may be characterized by spherical coordinates (κ n , ϕ n , θ n ). Now we define the choice of wave numbers κ n . The highest wave number is dependent on a mesh resolution, i.e. κ max = π/Δ, where Δ is equal to grid spacing, when we have an equidistant grid. In the case of non-equidistant grid, we may decide to pick a mean value of grid spacing as a Δ. The smallest wave number is defined by relation κ min = κ e /p, where κ e = A9π/(55L t ) is the wave number corresponding to the most energy containing eddies [6], p is the parameter which should be larger then one to make the smallest wave number smaller then κ e . The numerical constant A will be introduced later, and L t is a turbulent length scale. In [7] it is found that L t 0.1δ in is a suitable choice, where δ in indicates an inlet boundary layer thickness. Now, we may divide the wave number range, from κ min to κ max , into N equally large segments of size Δκ. The center of nth segment then corresponds to the nth wave number κ n .
We also assume an incompressibility of the turbulent field and hence, the relation is satisfied [3]. Therefore in a spectral space, the unit vector σ n is always perpendicular to the wave vector κ n , and its direction on the plane perpendicular to κ n is determined by the polar angle α n . As a consequence of isotropy, homogeneity, and incompressibility, the angles ϕ n , θ n , α n and ψ n are chosen randomly with probability distributions given in table 1 [8]. Table 1. Probability distributions of the random angels For a complete description of the turbulent velocity field u (x), see eq. (2), we have to determine the amplitudeû n of nth Fourier mode. Firstly, we introduce a turbulent kinetic energy where u rms is a root mean square value of the turbulent velocity field u . A relation between k andû n should be established using the equations (2) and (4) as in [3].
For the homogeneous isotropic turbulence, the energy spectrum E(κ) satisfies the following equality [9] ∞ 0 E(κ)dκ = k Now its clear, from the relations (5) and (6), that our approximation given by (2) causes that we can approximate the integral of energy spectrum by the finite sum of squares of Fourier mode amplitudes. Hence, the amplitude of nth Fourier mode is given bŷ The modified Von Kármán spectrum is employed to simulate the isotropic energy spectrum [3] (8) where κ η = 1/4 ν −3/4 is the Kolmogorov wave number, A is the earlier mentioned numerical constant set to fulfill the relation (6). Thus [8], For completeness, we should note that the symbols ν, and denotes a kinematic viscosity, and a turbulent dissipation rate, respectively. The latter quantity should be calculated as follows [4] = Now, we define a time sequence {u (x, t)} T t=0 of the random turbulent velocity field u (x) (see eq. (2)), where T is a total time. This time sequence corresponds to the generation of fluctuations at each time step. Until, the individual members u (x, t) of the sequence are independent of each other, the time correlation will be zero, which is unphysical.
For introducing of time dependency, we adopt the approach proposed in [5], where the field at current time step is computed as a weighted sum of the field at previous step, and the random field generated by (2). Thus, the time correlated random turbulent velocity field where Δt denotes a time step, a = exp(−Δt/T t ) is the weighting coefficient by which, the turbulent time scale T t may be prescribed, and the choice of the coefficient b = √ 1 − a 2 ensures a preservation of the root mean square values of the random turbulent velocity fields [6], i.e. v rms = u rms .

Comparative characteristics
As a consequence of isotropy and homogeneity, the Reynolds stress tensor u i u j is diagonal [3] and hence the following equalities should be verified.
where denotes a statistical mean, and δ i j is the Kronecker delta symbol. It is clear that with the increasing number of Fourier modes N, and time steps N t = T/Δt, the difference between the left-hand sides (simulation) and the right-hand sides (theory) of (12) is reduced. Note that in the result section, sec. 4, we will compare the normalized values of (12) averaged over all the axis directions. Now, we introduce two dimensionless parameters on which these equalities (12) also depends. The first one is the resolution parameter where Δx is a grid step. This parameter determines the number of cells used at the distance of the one turbulent length scale. The second characteristic is the size parameter where n p is the number of cells at inlet plane in one axis direction (will be illustrated further). The size parameter indicates the number of resolved large scale eddies (turbulent length scales) along one axis of a computational domain.

Computational domain
The random turbulent velocity field u(x, t) is generated on the inlet plane (y − z plane). The square equidistant mesh was chosen with n p = N s N r , see eq. (14), cells in the wall-normal (y) and the spanwise (z) direction. Due to the equidistance, the grid step Δx = L t /N r , see eq. (13), is the same for the whole domain. Hence, the domain size

Influence of resolution and size parameters
When testing the effect of the parameters N r and N s , we use the following setup. The total time is T = 0.4 s, and the time step is set to Δt = 2·10 −3 s, which determines the number of time steps N t = 200. As a consequence of homogeneous turbulence, we set the constant turbulent kinetic energy k = 13.5 m 2 s −2 specifying the root mean square of turbulent velocity component, see eq. (4), u rms = 3 ms −1 . When the turbulent length scale is set to L t = 0.05 m for the whole inlet plane, then the turbulent dissipation rate is also constant, i.e. = 992.04 m 2 s −3 , see eq. (10). The factor determining the smallest wave number κ min is p = 5, and the kinematic viscosity ν = 15.29 · 10 −6 m 2 s −1 is set for air at the temperature τ = 20 • C [10].  Figure 1 illustrates the dependence of the statistical value u rms ≡ u i u i /3 (using Einstein summation convention) on the parameters N r and N s for different numbers of Fourier modes N. We may see an enlargement of a suitable area (the area with maximum of two percent error compared to u rms ) with increasing number of Fourier modes. Specifically, the suitable area is approximately two times wider for N = 100, than the suitable area corresponding to N = 50. What may seem unexpected is the displacement of the lower boundary of the suitable area upward, which should denote the increasing of the minimal resolution N r with the increasing number of Fourier modes N. Furthermore, we can notice that from a certain value N s , the contours are nearly horizontal, i.e. the statistical value u rms is almost independent of the parameter N s .  we choose the smallest value of N from which the statistical values, u and u rms , have almost a constant behavior. Note that the computational complexity increases linearly with the number of Fourier modes N. Thus, the number of Fourier modes N = 200 seems to be an optimal choice. Further, the variance of the statistical signals in figure 2 is smaller for N t = 200 than for N t = 50, i.e. increasing of the number of time steps N t causes the smoothing. In addition, the signal of u rms slightly underestimates the prescribed value of u rms for N ≥ 180 approximately. It may be caused by moving the suitable area upward with increasing number of Fourier modes N (see figure 1) and thus, it causes that the point, N r = 30 and N s = 2, gets outside the suitable area. The dependency of the statistical values u and u rms on the number of time steps N t is illustrated by figure 3. Due to the previous result (see figure 2), we selected the number of Fourier modes N = 200. Furthermore, we chose the resolution and size parameter from upper part of the suitable area (see figure 1), i.e. N r = 50 and N s = 4. It is seen that the increasing of the number of time steps leads to a smaller difference between the statistical values and the theoretical values. We can also identify the minimum number of time steps needed to stabilize the statistical values. Thus, N t = 500 is sufficient to stabilize the statistical value u rms . On the other hand, u requires much more time steps for this purpose, i.e. at least N t = 5000.

Time correlation
It is well known that the random turbulent velocity field u (x) (see eq. (2)) is correlated only in space. Thus, we tried to establish a time correlation by relations (11), which is the conventional method to prescribe the turbulent time scale T t of the turbulent velocity field. Hence, the theoretical time correlation of v i should be equal to exp(−Δt/T t ) [6].
So, we set the turbulent time scale T t = 0.05 s, and then at the figure 4, we compared the time correlation of v 1 (x, t) by using the normalized auto-correlation function B norm ii (t ) = (1/v 2 rms ) v i (t) v i (t +t ) wheret is a time separation, with the theoretical time correlation. We can see that both time correlations are in good agreement.

Conclusion
This paper describes the stochastic method for generating random turbulent velocity field with the space-time correlation. The purpose was to set the physically correct inlet boundary condition for unsteady simulations. Above all, we tested the effect of various parameters on the conservation of the first and second statistical moments of the random turbulent velocity field. Specifically, the influence of the resolution parameter N r , and the size parameter N s was shown. We demonstrated that the suitable area is changing due to the changes in the number of Fourier modes N. Moreover, we are able to select the optimal value of N. The described stochastic method preserves the prescribed space-time correlation and thus it seems to be a promising method for the purpose of the inlet BC.
However, for some real tasks, the main disadvantage of this method may be the need to know the energy spectrum in advance.