Identification of the heat transfer frequency response in pulsating laminar and subcritical flow across a cylinder

The steady-state heat transfer from a cylinder in cross-flow is a prototype problem in thermo-fluiddynamics. However, in many applications such as the Rijke tube, the flow may fluctuate. This work analyses the phenomenon combining numerical simulation with system identification. Direct numerical simulation of laminar flow and Large Eddy Simulation at subcritical flow at Reynolds number equal to 3900 are used, respectively. Fluctuations of the inlet velocity in the simulation are excited over a wide range of frequencies. Time series of unsteady heat release and velocity are post-processed to identify dynamic models, which may be represented as transfer functions. They accurately describe the dynamic behavior and can be used for further modeling.


Introduction
The heat transfer between a cylinder and a steady fluid stream perpendicular to the axis of the cylinder has been studied in great detail. However, in some applications the flow is not steady, but fluctuating. In most cases, like it is in this work, the motivation was the Rijke tube [1,2]. A tube in which a constantly heated wire is exciting acoustic oscillations. This arrangement is the simplest setup of a thermo-acoustic system often used to study the phenomenon.
Literature on the dynamic behavior of the heat transfer to perturbations of the flow is scarce. The first and probably most notable publication was written by Lighthill [3]. The fundamental result was that the heat transfer behaves like a first order time lag system with time constant 0.2 d/u 0,∞ where d is the diameter of the wire and u 0,∞ the unperturbed mean flow velocity. Gersten [4], following Lighthill's approach, also calculated the response of the heat release to velocity perturbations. Again using the boundary layer equations Telionis and Romaniuk [5] solved for perturbations on a linearly decelerating flow numerically. For creeping flow at very low Reynolds numbers Bayly [6] was able to derive an analytical expression for the fluctuating heat release based on the Oseen equations. Kwon and Lee [7] also approached the problem numerically and with regards to the Rijke tube. Instead of solving the boundary layer equations they computed the linearized incompressible Navier-Stokes equations in frequency domain.
Although other researchers have approached this problem, Lighthill's expression for the time lag is still widely used to quantify the dynamic behavior. Today's computational methods allow us to use a different approach to this problem. Direct numerical simulation for laminar flow (Re = 0.4, 4 and 40) and Large Eddy simulation for subcritical flow (Re = 3900) is used to gather data and identify a data based black-box model for the heat transfer response to velocity perturbations. Three different Reynolds numbers for laminar flow are chosen to represent the regimes of creeping flow, flow without recirculation and flow with steady recirculation zone at the lee side of the cylinder respectively. The Reynolds number of 3900 represents subcritical flow and is used in many studies on cylinders in cross-flow so validation data is available in great detail. This method was already successfully applied to the laminar problem by Föller et al. [8]. Figure 1 depicts a sketch of the problem of fluctuating heat release rate from cylinder to fluid in pulsating cross-flow. The flow around a cylinder can be described by the Navier-Stokes equations. Instead of using the full set, some assumptions are made a priori. Fluctuations are assumed to travel instantaneously i.e. the cylinder is acoustically compact. For frequencies and mean flow velocities considered here this is a valid approximation. The temperatures of the cylinder and the far-field are assumed constant yielding a constant temperature difference ∆T . Also, fluid properties, namely density ρ, viscosity ν, heat conductivity λ and heat capacity c p are assumed to be constant. Invoking these assumptions, the incompressible Navier-Stokes equations can be used.

Problem description
Characteristic scales are the cylinder diameter d and the unperturbed free stream velocity u 0,∞ . The dimensionless temperature is θ ≡ T /∆T . Dimensionless numbers, namely Reynolds, Prandtl and Nusselt number, can thus be given as where n denotes the wall-normal direction of the cylinder surface scaled with d, and α is the heat transfer coefficient. Reynolds and Nusselt number can be seen as the linear combination of a steady value (index 0 ) and a fluctuation ( ). A transfer function linking these fluctuations can be defined as The Strouhal number Sr is a non-dimensional expression for the frequency of the fluctuation.

The CFD/SI approach
The method to combine CFD and system identification was first proposed by Polifke et al. [9] and has since been applied to several problem sets like laminar and turbulent flames, aeroacoustics and thermo-acoustics [10,11]. Figure 2 presents the procedure schematically. In the CFD/SI approach CFD simulations are conducted to acquire data which is subsequently postprocessed using methods from the field of system identification (SI) to achieve a model for the   On the mean flow velocity at the inlet a fluctuation in time is imposed. This fluctuation is created as random perturbations but fulfills some requirements in order to yield good results. The variation is less than 30% in order to yield a linear response [12]. Power spectral density (PSD) is high for the range of frequencies that is of interest 0 < Sr < 40. As the heat release behaves like a low-pass filter, higher frequencies only have small response amplitudes and will not be regarded further.
During the numerical computation the temporal development of the velocity fluctuation is acquired at a plane upstream of the cylinder to provide an input for the identification. The heat release at the cylinder is extracted as the wall normal temperature gradient which serves as output. The identification eventually yields a linear model that can be represented as frequency response (transfer function, equation 2).

CFD simulation
For the simulations, the PIMPLE algorithm, implemented in the OpenFOAM framework [13] (version 2.3.0), is used. It solves momentum, temperature and pressure correction equations for every time step. This is achieved by computing several outer corrections (SIMPLE loop) and two inner corrections (PISO loop), until a specified initial residual of 10 −6 is reached for pressure, velocity and temperature.

Direct numerical simulation of laminar flow
For the laminar flow cases a 2D representation of the domain was chosen and only half of the cylinder had to be considered due to symmetry [14]. A grid independence study was performed to determine domain and cell size and assess the mesh. All domain borders except the symmetry line parallel to the flow passing through the cylinder were placed at a distance of 50d to the cylinder to ensure that the solution close to the cylinder is not influenced. Grid independence for steady-state flow was achieved using a block structured grid with 72 square cells along the perimeter of the half cylinder ( figure 4 shows the unrefined mesh with 36 cells). In addition to grid independence at steady-state conditions the response to high frequency oscillations was assessed. The half-cylinder is discretized by 144 cells at Re = 0.4 and 4 and by 288 cells  The flow enters the domain at the inlet with specified mean velocity and superimposed fluctuation at constant inlet temperature. Pressure gradient is set to zero. Boundaries parallel to the flow are symmetry conditions, where all normal gradients vanish. The cylinder is represented as a non-slip wall with a fixed temperature ∆T higher than the inlet temperature. Pressure gradients vanish at the cylinder and the relative pressure is given as zero at the outlet. Momentum and energy are allowed to leave the domain under zero gradient conditions. An accurate solution for every time step is achieved for small time steps determined by a maximum allowed CFL number of CFL = ∆t u/∆x ≤ 0.2. The duration of the simulation corresponds to the flow passing the cylinder at least 200 times.

Large Eddy Simulation of subcritical flow
In contrast to direct numerical simulations where no further modeling is introduced Large Eddy Simulations can be performed for turbulent flows. The advantage of LES is the vast reduction of computational cost considering that cell size and time step in a DNS have to in the order of the smallest turbulent length and time scales. In this work a dynamic one equation eddy sub-grid scale model [15] was chosen to model the unresolved portion of turbulence.
The computational domain was chosen smaller than in the laminar case and spans now 15d in the directions perpendicular to the cylinder axis, except for the outlet border which is 30d downstream. This is justified by the fact that the boundary layer, i.e. the part of the flow-field strongly influenced by the flow is smaller than in the laminar case. The direction along the cylinder axis is resolved in the domain spanning πd [16,17].
The 3D mesh for LES is very similar to the 2D mesh for the laminar simulations shown in figure 4. It was adapted to the smaller domain size and the span-wise direction was discretized. Three refinement steps in concentric regions with decreasing radius around the cylinder ensure that cells are refined in the direction of the cylinder axis, yielding a total of 216 cells at the cylinder in span-wise direction. Additionally, the grid close to the cylinder experiences a grading in radial direction so the cells form layers around the cylinder with a non-dimensional wall distance of y + ≈ 1 for the innermost cell at an aspect ratio of 6. The circumference of the full cylinder is discretized using 216 cells leading to a total of 3.7 · 10 6 cells.
The computational grid was assessed for steady-state flow and simulation results were compared to a mesh a total 16.6 · 10 6 cells (refined by a factor of 5/3) and to experimental data. The drag coefficient of the cylinder for the mesh used here (c d = 1.029) was about 5% higher than simulated with the finer mesh (c d = 0.98) and 4% higher than reported in experiments (c d = 0.99) which corresponds to one standard deviation of the fluctuating quantity. Gradual Although the comparatively high Reynolds number leads to a very small Stokes layer at high frequencies, at least one cell was within a distance δ to the cylinder. Boundary conditions were chosen the same as in the 2D laminar cases. Additionally, the boundaries in span-wise direction are specified as periodic boundaries. To reduce computational time a maximum CFL number of CFL = 0.5 is applied to determine the time step size.

Linear system identification
System identification is the exercise of finding a model that describes a certain process on the basis of data. It has been widely used in contexts of electrical and chemical engineering and economics. An overview over the methods can be found in various textbooks [18,19,20].
In this case a linear dynamic model is sought from the time series acquired from CFD simulations. An input-output structure was chosen to represent the model because it is convenient for single input single output (SISO) systems. A valid model can only be found under the condition that the process is linear and time invariant (LTI) which is the case as long as only small perturbation amplitudes (up to 30% of the mean flow) are considered.
The model structure chosen in this work is often called a Box-Jenkins (BJ) model. In discrete time t = kT s (T s is the sampling time) it takes the form with the polynomials B(q) = b 0 + b 1 q −1 + b 2 q −2 + · · · + b (n b −1) q −(n b −1) and F (q) = 1 + f 1 q −1 + f 2 q −2 + · · · + f n f q −n f for the deterministic part and C(q) = 1 + c 1 q −1 + c 2 q −2 + · · · + c n b q −nc and D(q) = 1 + d 1 q −1 + d 2 q −2 + · · · + d n d q −n d to determine a model for the noise e[k]. q is the backward-shift operator (q −n x[k] = x[k − n]) and n b , n c , n d and n f are the orders of the polynomials i.e. the model orders. Figure 3 shows this type of model as a diagram. In the case C = D = 1 it is known as output error (OE) model where the noise is just passed on to the output and no disturbance model is produced. In both cases, a prediction error method is used to determine the parameters b 0 , b 1 , etc. An iterative algorithm applying a least squares approach minimizes a cost function which is essentially the squared sum of the differences between the known output from numerical simulations and predictions of the model. The identified models are validated subsequently to determine whether they correctly describe the dynamics of the process.
The system identification toolbox of the commercial software MATLAB [21] was used to perform post-processing and identification. For laminar flow represented by the cases up to a Reynolds number of Re = 40, OE models are identified. This model yields the best results for cases where the signal to noise ratio is very high.
Several test were performed on the prediction errors and the model itself after identification. In a final cross-validation test the output of the model is directly compared to the output of CFD simulations with the same input. All models for laminar flow perform better than 95% normalized root mean square fit using model orders of n b = 5 and n f = 4 for Re = 0.4 and 4 and n b = 4 and n f = 3 for Re = 40.
In the regime where turbulent fluctuations are present the full BJ model structure was chosen also fitting a model for the turbulent noise and the influence of wake structures. The model orders that proved best in this case were n b = 4, n f = 3, n c = 7, and n d = 7. A cross-validation is not possible here, because the stochastic fluctuations can not be predicted and always impede a comparison with time series from turbulent simulation.

Results
The results can be depicted as transfer function representations of the identified linear models. Figure 5 shows amplitude and phase of the complex frequency response for Re = 0.4, 4, 40 and 3900. In the regime of creeping flow (Re = 0.4) the steady-state gain is about 0.26 and the amplitude drops to small values rapidly with increasing frequency. Phase also shows a steep gradient at low frequencies and then approaches a value above −π/2. Gain values increase with Reynolds number as shown for Re = 4 and Re = 40. Additionally, a peak becomes visible in the region 0 < Sr < 1 appearing at higher frequencies and more pronounced for higher Reynolds numbers.
Similarly, a small phase lead occurs for very low frequencies at Re = 4 and stronger for Re = 40. At those Reynolds numbers, after the initial phase lead, a growing phase lag develops but with a gradient less steep than in the case of creeping flow. It also approaches a value above −π/2, but faster than in the case of creeping flow at high frequencies.
In subcritical flow at Re = 3900 the overall behavior is very similar to Re = 40. The steadystate gain is higher but no peak occurs at low frequencies. Also, no phase lead is visible and the high frequency limit is higher than in the laminar cases. Initial drops in amplitude and phase are steeper than in the case of Re = 40.

Discussion
The results for all Reynolds numbers basically feature the expected shape of a low pass filter indicating the typical behavior of a system with time lag. Flow regions influenced by the presence of a cylinder are much larger at low Reynolds numbers as are diffusive effects. Thus, time scales are larger than in the case of higher Reynolds numbers and gain drops at even smaller frequencies. For cases Re 1 overall time scales are similar and within the order of 0.2 d/u 0,∞ predicted by Lighthill whose theory relies on the existence of a laminar boundary layer.
The peak in amplitude occurs at frequencies where the size of the unsteady boundary layer (Stokes layer) is similar to the size of the steady-state boundary layer (Prandtl layer). This suggests that an interaction between both layers amplifies the gain compared to the quasisteady state. The larger the Reynolds number, the more pronounced this effect becomes. There might be a possibility that a peak also exists at lower Reynolds numbers but less distinct and shifted to even lower frequencies. To capture this in the system identification, very long time series would be necessary in combination with very high model orders. The latter would have other negative effects as over-fitting might occur and is therefore not desirable. For subcritical flow no peak is visible. Most likely this is due to the fact that large amounts of turbulent noise and the vortex shedding, which effects the heat release at about Sr ≈ 2.4, complicate the fit of the deterministic model even if also the noise is modeled. The interval indicating 95% confidence for the frequency response is as large as [0.48 0.59] at the lowest frequency but becoming smaller quickly with rising Sr.
The low frequency limit i.e. the quasi-steady regime is a good indicator for the model performance and can be assessed approximately analytically [22]. For the heat exchange between fluid and cylinder in steady cross-flow, Sparrow et al. [23] proposed a correlation, valid for 1 ≤ Re ≤ 10 5 of the form for constant fluid parameters. Assuming, that this expression also holds in the quasi-steady limit (Sr → 0) an approximation for the limit of the transfer function can be found from a linearization. In the case of small amplitudes the transfer function can be approximated by  [24] which was developed for 0.02 < Re < 44 and linearized as above. The quantities can be compared to the results from the identification procedure.
The low frequency limit is challenging using the CFD/SI approach to obtain a linear model because simulations have to be run for a long time to satisfyingly represent low frequencies. However, this comparison shows that the agreement with known correlations in the quasi-steady limit is quite good.
The phase of the frequency response experiences a small phase lead at low frequencies. Again, this is not captured by the identification for Re = 0.4 because the effect will be even less intense as at Re = 4. At subcritical flow, phase lead is also not visible for the same reason as the peak in gain does not occur. This phase lead, however, is due to a phase lead in the velocity fluctuations in the boundary layer which is overcome by inertia at higher frequencies. The phase is directly linked to the time scales that govern the process indicating again larger times at smaller Re. In the limit of high frequencies, where the gain is already low, identification of the phase is less precise. For all models phase values approach a value above −π/2. Assuming a fist order time lag, a limit of −π/2 would be reached for Sr → ∞. In this case the lee side of the cylinder has a small contribution which acts with opposite phase with respect to the rest of the cylinder. Thus this theoretical limit is never reached for the overall heat release rate.

Conclusions
Models for the response of the heat release driven by a temperature difference between fluid and a cylinder in cross-flow to velocity perturbations were identified in unprecedented detail. The overall behavior found by some researchers was confirmed but several new insights could be gained from this study. The response is severely depending on the Reynolds number of the mean flow at lower Reynolds numbers represented by Re = 0.4, 4 and 40. It was also possible to identify a model at subcritical mean flow conditions at Re = 3900. Although much more demanding in terms of computational cost and to the tools used for identification also the dynamic heat release behavior in this flow regime could also be captured.
Models as presented can easily be transformed between time domain and frequency domain representations as well as between discrete time and continuous time. Therefore, those models can be used not only to study the dynamic heat transfer behavior of bluff bodies in general, but also as an input in other applications. Among them e.g. network models to determine the thermo-acoustic stability of a Rijke tube or domestic boilers.