On computations of temperature dependent incompressible flows by high order methods

Parallel vortex shedding from the heated cylinder as dependent on the cylinder temperature is investigated in numerical simulation. The computational scheme for the system of the incompressible Navier-Stokes equations coupled with the heat equation is presented. The scheme is based on the mixed implicit-explicit solver of the third order and the spectral element method for the spatial approximation. The results are compared with the experiments for flow of air and water. The Strouhal-Prandtl-Reynolds relationship and evolution of the separation angle as dependence on the temperature is the final result of the computations.


Introduction
The vortex shedding from the circular cylinder is widely studied phenomena both experimentally and in the numerical simulations. This, geometrically simple, problem demonstrates various effects in the transition of the flow from the laminar to turbulent behaviour. From the wide literature concerning this problem we have already knowledge about the dependence of the vortex shedding frequency on the velocity magnitude and vortical structures in the wake. In the present study we restrict our attention to the flow regime, where the vortex shedding parallel to the cylinder occurs, and where a 2D numerical simulation is justified. This regime is achievable in the range of Reynolds numbers Re c < Re < 160, if Re c denotes the critical Reynolds number, separating the steady and unsteady flow. In this flow regime we study the influence of the cylinder wall heating to the vortex shedding frequency (described in terms of the Strouhal number S t). The Reynolds number is here defined as Re = |u ∞ | D ρ ∞ μ ∞ , where D stands for the cylinder diameter, μ is the dynamic viscosity, ρ density, |u| the characteristic velocity magnitude and subscript ∞ denotes values in the far-field/inlet. The Strouhal number is shedding.
The fluids considered are water and air. For our study is most important that these two substances exhibit opposite temperature dependence of its viscosities. We expect the flow to be incompressible in case of both studied fluids, because the velocity magnitudes satisfy M 2 1 (M is the Mach number). Effects of the buoyancy and viscous heating are neglected, since it has a minor influence to the flow in our setting. Therefore the influence of the flow by the temperature distribution is given by the temperature variation of the dynamic viscosity μ and the thermal conductiva e-mail: jpech@it.cas.cz ity λ. Frequency of the vortex shedding depends mostly on the viscosity, but variability in the thermal conductivity is also important, especially for preserving the correct local value of the Prandtl number (Pr = c p μ ∞ /λ ∞ ), which moderates the heat equation. For the isothermal case, when the cylinder wall has the same temperature as the inlet fluid, the results of numerical simulation using the scheme of [1] give excellent agreement with S t-Re relationship (see e.g. [2]). For the study of the heating effects to the S t-Re relation, the scheme for isothermal case is enhanced using technique already proposed in [3] and applied in the same manner to the heat equation as in [4]. In comparison to the previous results of [4] using the first order scheme, we present the third order scheme and optimise the computation in the spatial solver.
The flow parameters were chosen to correspond with those of [5], where an influence of heating to the Strouhal number was studied experimentally. The Strouhal number obtained from our computation in comparison to the values from the experiment then provides a basic information about the relevance of the computational results.

Physical model
A system of equations describing the motion and temperature distribution of an incompressible fluid with the temperature dependent viscosity and thermal conductivity may be written in the dimensionless form as where u ∈ [0 : 1] denotes velocity, p is dimensionless kinematic pressure, μ(T ) is the temperature dependent dy- This is an Open Access article distributed under the terms of the Creative Commons Attribution License 4.0, which permits distribution, and reproduction in any medium, provided the original work is properly cited.
XQUHVWULFWHG XVH To complete above mathematical system the functions λ(T ) and μ(T ) are approximated as power laws based on the data [6] and procedure suggested in [2], respectively. Exponents M and L, used in the computations, are listed in the table 1.
The whole system is fully coupled by the variable properties. The technique used in the numerical scheme decouples the system as the operator splitting method together with linearization and results in acceptably accurate and fast solver.

Numerical methods
The computational scheme for the system (1)-(3) originate from the high order splitting scheme [1] enhanced by the approach for variable viscosity ( [3]) and was already presented in its first-order version in [4].
The temperature dependent properties of viscosity and thermal conductivity are divided into a constant and variable (subscript s) part The terms with variable parts (μ s , λ s ) are left to the explicit evaluation together with the convective non-linearities. The constants μ ∞ and λ ∞ stay on the place as in the constantproperty system. Introducing a notation of a function value at the n-th time step, f (n Δt) = [ f ] n , we can write the explicit part of the scheme for the incompressible Navier-Stokes equations Evaluation of the viscosity follows [μ s (T )] n = μ s ([T ] n ) and coefficients α q , β q are listed in the table 2. The pressure Poisson equation as known from the fractional step methods follows and is accompanied with the high order pressure boundary condition (HOPBC), which has crucial impact on the final accuracy of the scheme The condition in the form 9 expects a constant velocity on the boundary, what is satisfied in our model. The second step (8) allows use of the spectral element method as the spatial solver. It also applies the incompressibility to the second intermediate velocity fieldû, which satisfŷ and finally forms the right hand side of the Helmholtz equation for the velocity at the new time step The temperature dependent thermal conductivity results in a strongly non-linear heat equation, which is linearized and splitted in the same manner as the equation of motion. The operator splitting is applied to perform linear terms implicitly and the non-linearities, containing convective term and variable part of the diffusion term, are evaluated explicitly. The explicit step of the third-order scheme then results in Implicit part has again the form of the Helmholtz equation 1 Pr Re and is solved by the spectral element method.

Results
The computational domain with dimensions measured from the center of the cylinder had 51D upstream, 100D downstream and 51D on both sides from the cylinder. A refinement towards the cylinder and the area of the wake was 02089-p.2 Here we denoted C the cylinder surface and T = −pI + μ(∇u + (∇u) T ) the stress tensor. The shear stress was computed as a tangential component of the velocity gradient in the direction of the normal n to the cylinder wall Values of C D , C L and positions of all detected zeros of the shear stress were saved in every time step.

Temperature dependence of St-Re relationship
If Re > Re c , a disturbances from initialisation of the computational process are not suppressed and the flow develops to the vortex shedding regime. History of C D and C L reflects the evolution of the flow and must be processed, since only the fully developed vortex street contains the frequency, which refers to the Strouhal number. Experimental results of [5] are denoted as exp, other notation is adopted from figure 2.
increment in case of air and increases in case of water.
Presently computed data coincide with those computed in [4] by the polynomial approximation up to order 49. The new results therefore support applicability of very high order spatial approximations and shows that distinction of the computation from the experiment is not given by the spatial approximation with mentioned settings.

Influence of heating to the separation angle
The separation angle Θ is measured on the arc defined by frontal stagnation point, cylinder center and the sepa-  sulting angles as dependent on the Reynolds number and temperature are plotted in figure 5 and figure 6. The frontal stagnation point very accurately preserves mean value of Θ = 0 and oscillates with moderate amplitude < 2 degrees. Time-averaged positions of the separation points (E1 and E2 in figure 4) appear symmetrically (Θ S 1 360 − Θ S 2 ). A difference between angles evaluated on the present mesh and mesh for very high order basis of [4] was observed. This deviation in results is up to 2 angular degrees and suggests a need of improvement of the present evaluation method, whose accuracy depends on lo-cal distribution of quadrature points, which is non-uniform over the elements.

Conclusion
There are still various settings in the computation, which need to be revised in relation to the results of the experiment. Values of the final Strouhal number are strongly influenced by the choice of the reference temperature and dependencies of μ and λ. Even though the buoyancy effects are negligible, the model expecting a constant density is limiting for the case of the heated flows even if the velocity magnitudes in the flow are small. Successes of the present computations are in confirmation of relevance of the previous very high order computations [4] and inclusion of the variable properties into the solver capable of the fast and accurate computations of the fluid flow. Results of the computations clearly shows the impacts of heating in both settings with properties of water and air. Accuracy of the computation is comparable with tolerance of the experimental data. Observed deviation of the absolute values of computation and experiment are left for wider discussion. Computations of the separation angle shows its dependence on temperature and confirms its dependence on the flow velocity. Observation of this quantity has many complications in experimental studies and the computation clarifies the phenomenon especially in its time dependence.