Scheme for Evolutionary Navier-Stokes-Fourier System with Temperature Dependent Material Properties Based on Spectral/hp Elements

The computational algorithm for evolutionary Navier-Stokes-Fourier system with temperature dependent material properties (density, viscosity and thermal conductivity) is presented and tested. As a consequence of the thermal expansion, the velocity field is not divergence free. The system fully couples the momentum balance with non-linear advection-diffusion equation for temperature. The model is suitable for low speed flow simulations of the Newtonian, calorically perfect fluid, which obeys Fourier law. The algorithm is an extension of a well known semi-implicit scheme for incompressible Navier-Stokes equations in primitive variable formulation. Performance is tested on manufactured solution, what gives an overview to the temporal convergence. Comparison with experimental data is also presented.


Introduction
This work presents a numerical algorithm for the system of the Navier-Stokes equations coupled with the balance of internal energy where v = [u, v, w] T is the velocity vector (by setting w = const. = 0 we restrict to 2D problem), p is a variable related to the thermodynamic pressure, 1 T denotes the temperature, D = 1 2 ∇v + (∇v) T is the symmetric part of the rate of strain tensor, constant Re is the Reynolds number and constant Pr is the Prandtl number (for sake of simplicity we set Re = Pr = 1 for the testing on exact solution). The fluid is expected to be (calorically) perfect, 2 Newtonian, 3 whose heat flux obeys the Fourier law. 4 In system (1), we consider those fluids, which become nonhomogeneous in variable temperature fields due to temperature dependence of its material parameters, namely the density ρ = ρ(T ), dynamic viscosity μ = μ(T ) and thermal conductivity κ = κ(T ).
Without loss of generality, solving (2) instead of (1a), we avoid specification of the second viscosity coefficient λ.
The forcing terms f v , f T , may represent action of volumetric forces, e.g. gravity or viscous heating, but m is set zero in most of realistic situations. In case of testing of our algorithm on a given solution [v e , p e , T e ] T , we construct the forcing terms such, that Eqs. (2), (1b) and (1c) are satisfied.
Our computational scheme is developed for simulations based on the spectral/hp element approximation in spatial coordinates. We use the polynomial approximations of degree 15 in our tests, what eliminates the numerical error in spatial coordinates and we are getting an overview of error production, which belongs directly to the algorithm/discretisation in time. The high order spatial approximations also naturally include approximations of higher-order derivatives, what is utilized in the scheme.
Stationary and non-stationary models are denoted stat. and nonst., unspecified variability of a property is denoted var.

Algorithm
Our approach is inspired by the velocity-correction scheme with the high order pressure boundary condition (HOPBC) proposed for the incompressible Navier-Stokes equations in [7]. The constant property case, [7], is widely used for its efficiency and was already extended to problems with variable viscosity in [6].
Its modification was used also to the incompressible Navier-Stokes-Fourier system with temperature dependent viscosity and thermal conductivity in [10]. Efficiency of the approach comes from the implicit-explicit (IMEX) formulation, which allows decoupling of the system. The main contribution of the present work, which is a continuation of [10], is in extension to the problems with temperature dependent density. However, the velocity divergence cannot be further neglected in the momentum balance, what is the substantial difference from the previously discussed models and algorithms.

Decoupled System
We use the IMEX scheme in which the Backward difference formula (BDF) of order Q approximates the temporal derivative and a consistent extrapolation is applied to chosen terms (N) In (3), u is the searched solution, L denotes the terms solved implicitly, which we expect to be constant in time. Subscript n + 1 (or operator in square brackets with subscript) denotes evaluation at time t n+1 = t 0 +(n+1) t, where t is the discrete time step. Coefficients {α q } Q−1 q=0 , {β q } Q−1 q=0 and γ for particular Q can be found, e.g., in [10]. Henceforward, we use ' * ' in the superscript to denote extrapolation, N * ≡ N * := Q−1 q=0 β q N n−q . The extrapolated terms are evaluated using data from previous time steps, , what allows separate/decoupled solution of the (generalized) Navier-Stokes equations (2)-(1b) and the non-linear energy equation (1c).
Solution during one time step may be summarized to the scheme 1. Update μ, κ, ρ, ∇ · v, and HOPBC using already known values {v n−q }

Balance of Momentum and Mass
The scheme decouples solution of the Navier-Stokes system (2)-(1b) to the pressure-Poisson equation and an elliptic equation for velocity. The equation for pressure is derived as a projection to the irrotational space by application of the divergence operator to (2) where we applied (1b), identities ∇ × ∇ × v = ∇∇ · v − ∇ 2 v and ∇ · ∇× ≡ 0, ∂v/∂t was substituted by BDF andv = The temporal derivative of the density, which is extrapolated in (4), is approximated by Q-th order BDF We denote the extrapolation of the derivative approximation by superscript ' * * '. Note, that we have to specify the initial value ∂ρ ∂t 0 or both ρ 0 , ρ −1 to initialise the scheme of the lowest order Q = 1.
Our model assumes, that the density is entirely determined by the temperature distribution. Then, the divergence of velocity, whose forward estimate, [∇ · v] n+1 , is required in (4), follows from (1b) The forward estimate of velocity divergence is the crucial step in the proposed scheme. HOPBC is the natural boundary condition for (4). It is derived as projection of the momentum equation (2) to the direction of normal n to the domain boundary ∂ The forward estimate of velocity divergence follows from (6) again. Similarly to (4), we approximate the acceleration term ∂v ∂t by the BDF of Q-th order, whose initialisation requires value ∂v is circumvented in many realistic simulations, which begin from a constant fields.
The solution of (4) gives estimate/prediction ofp n+1 and we can solve (2) as an elliptic problem for v n+1 . However, in the case of temperature dependent viscosity and density, the algebraic system derived for operators with variable coefficients has time dependent matrices, whose direct solution is inefficient. To preserve efficiency of the scheme, we split such operators to the time independent part, which is solved implicitly using a direct method and a variable part, which is extrapolated together with the non-linear terms. We introduce material properties in form whereμ andκ are time-independent, while μ i = μ i (x, t) and κ i = κ i (x, t). The variable density ρ = ρ(T ) acts in our scheme as an inverse value, c.f. (10), so the splitting is done accordingly.
To demonstrate the splitting, we consider the second order operator with variable viscosity and density Only the term with time independent operator 1 ρ ∇· μ (∇v) T is solved implicitly, while we apply extrapolation to terms containing variable parameters 1 ρ i and μ i . This approach is valid forμ =μ(x),κ =κ(x), resp. 1 is constant in space, the constant operator simplifies to μ ρ ∇ 2 v (resp. ∇ 2 v if the properties are normalized toμ =κ =ρ = 1, what is the case of (1), the balance equations in form independent of physical units).
The final form of the equation for velocity becomes

Balance of Energy
The energy equation with temperature dependent thermal conductivity is strongly non-linear. We split the diffusion operator to the time independent and the variable part, following the technique shown for the velocity-correction (10). We setκ = 1 for simplicity and the discretized energy equation (1c) gets form

Temporal Convergence on Manufactured Solution and Application
Our convergence tests are based on the method of manufactured solutions, alternative to estimates of numerical analysis on simplified system. This approach lacks generality, because we always restrict to particular data and some representative of the solution space, but we get rough convergence estimate for unrestricted equation system (1), while proving also the correctness of method implementation.
As an exact solution, we take a smooth functions v e :
The incompressible Navier-Stokes equations define the pressure up to a constant value and only the boundary condition for velocity is needed. In this sense, we set the Dirichlet boundary condition for velocity on whole ∂ . However, the pressure-Poisson equation (4) requires setting a boundary condition as a consequence of the decoupling. We set HOPBC (7) at ∂ and solve the fully Neumann problem, which defines the solution up to a constant value, which we set by fixing the solution to zero in one of the grid points. The boundary condition for pressure is an artificial element of the computational scheme and its existence is related to the splitting error. The boundary condition for energy equation (11) is of Dirichlet type for whole ∂ .
We present the first and second order schemes in time in the convergence tests. The technique is applicable to higher-order schemes as well. A multi step schemes use data from multiple time steps, what complicates its initialisation. We apply the first order BDF method for initialisation of the second order scheme. The first order scheme needs data of only one backward time step, but the time step must be appropriately shortened.
As mentioned already, the acceleration in HOPBC (7) and the term ∂ρ ∂t in (4) require an initial value or one other backward value for proper initialisation also in case of the first order scheme, what is in contradiction to standard initial conditions for system (1), which require only the initial values. However, setting the correct values for calculation of the first time step is crucial for the final accuracy of the solution.
Finally, we trace appropriate norms of difference between the exact and computed solutions on a set of computations with time steps t = t/2 n , t = 0.2, n = 0, . . . , 9 for t ∈ [0 : 1].
We use the power laws for approximation of dependence of material parameters on temperature The temporal convergence of the above scheme for α m = α k = α r = 0.1 , β m = β k = β r = 2 is shown in Fig. 1.
A detail view of error production, Fig. 2, shows, that the dominant error production arises at the grid point, which was used to set the unknown constant for the Neumann problem.
The scheme was successfully applied in a 2D simulation of flow around the heated cylinder and the results were compared with experimental data [15], where the dependence of the vortex shedding frequency (Strouhal number St) on the wall temperature of the cylinder, T W , was observed. Figure 3 shows the substantial difference in results between the model neglecting the thermal expansion, [10], and the present one. Fig. 4 shows value range and structure of velocity divergence in a chosen realistic simulation.

Conclusion
The numerical scheme proposed for the Navier-Stokes-Fourier system with variable parameters allows to solve the highly complex mathematical model, which has an impact to understanding the processes connected with the heat exchange, transport and energy storage in fluids.
The computational scheme for a fluid flows influenced by temperature as modelled by system (2), (1b), (1c) was developed and tested. The scheme was primarily constructed for spatial discretisations based on spectral/hp finite elements and presented results were obtained after implementation to the Nektar++ framework [2], modified version 3.3.
We did not impose restrictions to the type of functional dependency of the material parameters on temperature. Graph of error convergence in L ∞ norm, Fig. 1, results from testing on a manufactured solution and shows a good convergence properties of the scheme, what is promising for applications.
Considered model neglects compressibility in the sense of direct dependence of density on pressure, but the velocity field is not divergence free as a consequence of the thermal expansion. A forward estimate of velocity divergence is needed in the proposed scheme and its successful approximation is one of the main contributions presented in this work. For these reasons, the scheme is unique among numerical schemes based on the finite element approximations in space.
Proposed scheme is an extension of the efficient semi-implicit solver for Incompressible Navier-Stokes system [7] and it is suitable for a fast and highly accurate simulations of problems on long time intervals. The present results inspire implementation of high order BDF schemes and extension of the solver to 3 spatial coordinates.
Derivation of the scheme includes a number of sub-steps, whose detail description is beyond the scope of this article and will be published separately, together with extension of the scheme for energy equation with variable density and further testing of performance as dependent on various physical parameters in the equations.
Also the results from application of the scheme to computations of a physically realistic problem and comparison of its results with experimental data exhibit a good coincidence and will be presented with detail description in a separate article.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.