Clearance gap flow : simulations by discontinuous Galerkin method and experiments

Compressible viscous fluid flow in a narrow gap formed by two parallel plates in distance of 2 mm is investigated numerically and experimentally. Pneumatic and optical methods were used to obtain distribution of static to stagnation pressure ratio along the channel axis and interferograms including the free outflow behind the channel. Modern developing discontinuous Galerkin finite element method is implemented for numerical simulation of the fluid flow. The goal to make progress in knowledge of compressible viscous fluid flow characteristic phenomena in minichannels is satisfied by finding a suitable approach to this problem. Laminar, turbulent and transitional flow regime is examined and a good agreement of experimental and numerical results is achieved using γ − Reθt transition model.


Introduction
Fluid flow in narrow channels and gaps of various crosssections is a phenomenon occuring in many technical applications, such as clearance gap flow in control valves of steam turbines, tip leakage flow in gas and steam turbines, clearance gaps in screw compressors etc.There is a need to understand fluid flow characteristic phenomena in these types of channels and be able to simulate it properly.A simple straight narrow gap formed by two parallel plates is chosen as a basic simplification of these applications.The distance of the plates, and hence the characteristic dimension of the channel, is varied from 1 mm to 5 mm, length of the channel is 100 mm and width also 100 mm, so that side walls don't influnce flow in the centre of the channel.A few classifications of small channels according to their size can be considered.Following Mehendale et al. (2000), [1], we can classify the considered gap as a compact passage, which is a category between mesoand microchannels on one side and conventional channels on the other side.Probably more used classification has been set by Kandlikar and Grande (2003), [2].They define microchannels up to 200 μm, minichannels from 200 μm to 3 mm and conventional channels from 3 mm, while the characteristic dimension is diameter in case of circular tubes and the smaller dimension of cross-section in case of rectangular channels.With respect to this classification and the emphasis on the gap of 2 mm we will call the channels under investigation minichannels, with only partial overlap into conventional channels category.
The objective of this study is to examine properties of air flow in minichannels, especially turbulent/laminar regime of the flow, and simulate it properly.It is supposed to assist in future research of more complicated cases of curved or narrowing minichannels and microchannels, in which the authors are engaged for a time, [3].The flow a e-mail: prausova@ntis.zcu.czregime is influenced by many aspects as properties of the fluid, velocity of the flow, roughness of the walls, outer turbulence intensity and of course dimensions of the region.Computation of Reynolds number with respect to height of the channel gives values of the order of 10 4 , which suggests turbulent regime of the flow.However, while in macrochannels the flow passes from laminar to turbulent regime for Reynolds numbers of the order of 10 3 , in microchannels viscous effects are predominant and the flow remains laminar even for relatively high values of Reynolds number, [4].Situation in minichannels, which stand between these two categories, needs to be further studied.
Another subject of interest in the presented work is the free jet flow behind the minichannel, where the air leaves the channel and continues into a free space.First set of experiments and numerical simulations carried out on the minichannel showed that in spite of relatively high Reynolds number the flow in the minichannel is rather laminar, but it is obvious that the flow behind the outflow is turbulent and the unsteady phenomenons present here are influencing also the flow in the end part of the minichannel.Transition from laminar to turbulent flow regime occurs.
Experimental and numerical research of the outlined issue is done to make progress in knowledge of the compressible viscous fluid flow characteristic phenomena in minichannels.

Governing equations
To simulate the flow numerically, the problem is simplified into two dimensions.Such a simplification is enabled thanks to the width of the measured channel.The following mathematical model in dimensionless form is con- sidered.Favre-averaged system of compressible Navier-Stokes equations is written as where ¯ , ũi , p, ẽ are dimensionless mass-or time-averaged values of density, velocity, pressure and energy.Massaveraged stress tensor and Reynolds stress tensor are given by tij = 2μ S i j , τ i j = 2μ t S i j − 2 3 δ i j k ¯ , where S i j = 1 2 ∂x k .This system is completed by two-equation k-ω model of turbulence with cross-diffusion term according to Wilcox (2006), [5], composed by two transport equations for turbulence kinetic energy k and specific dissipation rate ω, in dimensionless form written as where i, j = 1, 2 and closure coefficients are set as α = 13 25 , β * = 0.09, β = 0.0708, σ = 0.5, and turbulent viscosity is given by Joining the transport equations of the turbulence model to the system of Favre averaged compressible Navier-Stokes equations, the following system of equations is obtained where w is a vector of conservative variables, f s (w) are the inviscid numerical fluxes, f v s (w, ∇w) are the viscid numerical fluxes and p(w, ∇w) is a production term.The vectors are written as In case of modeling transition from laminar to turbulent flow, the Langtry and Menter's γ −Re θt model, [6], [7], is used to complete the above mentioned equations.The model uses local variables and is based on two transport equations, one for intermittency γ and one for momentum thickness Reynolds number R e θt , which sets the transition onset.The transport equation for intermittency in dimensionless form is given as where the production and destruction terms are S is the strain rate magnitude and Ω is the vorticity magnitude.Empirical correlation for the transition onset F onset is based on vorticity Reynolds number Re ν and viscosity ratio R T and is given by following relations where y is the distance from the nearest wall and .

The empirical transport equation for momentum thickness
Reynolds number R e θt can be written as where the source term is defined is the time scale and function F θt is The boundary layer thickness δ is computed using where Ū is absolute value of velocity.The empirical dependence of parameter Re θc , the length of the transition F length and the transition onset Re θt on R e θt are used as given by Langtry and Menter (2009), [6].
The first two terms on the right-hand side in the transport equation for turbulence kinetic energy k (production and dissipation terms) are multiplied by intermittency γ.While multiplying the dissipation term, the value of intermittency should be within 0.1 and 1. Transport equations for γ and Re θt are added to the system of equations given above.

Discretization
The spatial discretization of equations ( 1) is performed using the discontinuous Galerkin finite element method (DGFEM), [8], [9], [10]. Let . ., K} be a triangulation of the computational domain Ω ∈ R 2 and Γ = Γ int ∪ Γ b be the set of inner Γ int and boundary Γ b edges.The solution of equations ( 1) is considered in the space of vector functions and P q (Ω k ) is the space of polynomials of q-th order at the most.
To discretize the second order derivatives using DGFEM, the auxiliary variable θ r = ∂W ∂x r for the gradient is used.The system of equations ( 1) can be rewritten into the system of first order equations Multiplying this system by a test function v ∈ S h , integrating it over the element Ω k and using the Green's theorem, the following integral identity is obtained The key part of the method is an approximation of curve integrals.They are evaluated in the same way as it is done in the case of finite volume method.The inviscid flux is evaluated as where H is Lax-Friedrichs flux, [11], n = [n 1 , n 2 ] is vector of outer normal, values w − resp.w + are the left resp.right limits of the vector w on the edge Γ and w b is a boundary value of w computed from boundary conditions.Similarly, the value of viscous flux and the value of w on the edge for gradient equation are evaluated using the appropriate numerical fluxes.In present paper, the approximation known as local discontinuous Galerkin (LDG) method was chosen, [12].
Let w k be a part of solution w corresponding to control element Ω k and ϕ k,i be i-th basis function of the space P q (Ω k ).The m components w m k , m = 1, 2, ...4, of vector w k can be expressed as a linear combination Inserting it into the integral identity and evaluating the volume and curve integrals using the Gauss integration rules of appropriate order, the following system of ordinary differential equations is obtained where vector W consist of coefficients of linear combination (2).Taylor polynomials are chosen as the basis functions of the space P q (Ω k ) where x s and y s are coordinates of element center.Time integration of equation ( 3) is performed using the third order Runge-Kutta method in the form Computational efficiency of explicit schemes is directly affected by the maximal size of the global time step, which is limited by the CFL condition of stability where λ xk , λ yk is maximum wave speed in x, y direction, Re is Reynolds number, Δl k is diameter of control element Ω k and p is the order of polynomial approximation.A local time-stepping improvement of explicit method was implemented to enhance computational efficiency, [13].
Numerical simulations are performed using computational environment Matlab and programming language Java.After ensuring independence of the results on the mesh resolution, the computational mesh consisting of approximately 65000 nodes was used.The use of Galerkin method increases the number of degrees of freedom up to almost 200000.Capturing of effects occuring in boundary layer is accomplished by sufficient refining of the mesh near the channel walls and in the free jet.

Experimental Facility
The experiments were carried out in the Aerodynamics Laboratory of the Institute of Thermomechanics AS CR, v.v.i. in Nový Knín.The laboratory is equipped with the wind tunnel of the indraft type and intermittent operation.
For the purposes of minichannel flow examination, a special measurement facility was designed and manufactured.The model of the facility can be seen in Figure 1.The facility allows to perform both optical and pneumatic measurements of the flow through channels of length l up to 200 mm and allows smooth adjustment of the slit height h in range from 0 to 5 mm.The proper height was set up using feeler gauge.Since sensitivity of optical methods depends on the width of the minichannel and also the assumption of the two dimensional flow needs to be satisfied reasonably, the width w was chosen 100 mm.The air enters the test section directly from the laboratory and flows through the minichannel formed by two parallel steel plates with polished surface (upper and lower walls) and by the side walls with glass windows (Figure 2).Then the air continues to the settling chamber, where the back pressure p B is maintained during the measurement using the regulating slide valve downstream the settling chamber.The settling chamber is connected to the vacuum storage by the pipe with the spherical fast-acting gate.Scheme of the set up of the experiment is depicted in Figure 3.

Methods of Measurements
The dimensions of minichannels are considerably limiting for any experimental investigation.Application of classical pressure probes as well as thermoanemometric probes is problematic since the dimension of the probes might be comparable to the height of the minichannel.Use of these techniques could cause significant disturbances in examined flow field, resulting in unreliability of the measured data.For the optical measurements, special care must be devoted to setting the parallelism between channel walls and the beam of the light.This is important particularly for identification of the interference fringes near the wall.Used measurement methods will be briefly described in further paragraphs.
Pneumatic measurements were carried out using static pressure taps at the centreline of the top wall of the minichannel.Exact position of the pressure taps can be seen in the Figure 4. Pressure taps near the exit are spaced more densely in order to capture expected higher pressure gradient in the streamwise direction in this region.It can be also noted that in the exit region the pressure taps are drilled with a certain angle of the minichannel axis.This arrangement ensures that the downstream pressure taps are not influenced by those which are positioned upstream and allows denser spacing.Also the pressure in the settling chamber and the barometric pressure was measured.
For the pressure data acquisition two 16 channel pressure scanners PSI 9116 were used.The scanner integrates 16 silicon piezoresistive pressure sensors and microprocessor that performs a correction for sensor zero, span, linearity and thermal errors.The manufacturer guarantees system accuracy of up to ±0.05% of the full scale.The connection between the pressure scanners and pressure taps was realised using the plastic tubes of diameter 0.5 mm and length up to 0.5 m.
Total pressure measurements were not carried out, due to the small dimensions of the channel that does not allow to use Pitot probe inside the channel for most of the investigated dimensions.
The optical methods are very useful for minichannel measurements since it allows to perform the measurements without the intrusion of any probe in to the flow.It is very important because the probe dimensions might be of the same order as a characteristic dimension (height) of the minichannel, and thus an intrusion of the probe can affect the flow and bias the measurements.
Aerodynamic Laboratory in Nový Knín is equipped with the Mach-Zehnder type interferometer.The scheme of the optical path is given in Figure 5.During the exposition, the light from the light source goes through the monochromatic filter and field lens to make a collimated beam of monochromatic light.Lasers, spark light generators or LED diodes can be used as a light source.In this case the spark light source was used providing sufficient light intensity and sufficiently short exposition time (2 • 10 −6 s) which allows to take a picture even when the flow is unsteady.The beam is then split by the beam splitter forming a reference and a sample beam.The sample beam further goes through test section, where the measured model is placed, meanwhile, the reference beam is reflected by the mirror to the upper branch of interferometer, where the compensating glasses are placed in order to retain the same optical path as for the sample beam.Then both sample and reference beams meet again and are reflected to the CCD camera through some other optical elements.
When there is no flow in the test section and the mirrors are adjusted exactly in line and parallel, then both beams have the same phase shift and no interference patterns can be observed on the screen (camera).This set-up is known as an Infinite-Fringe technique.During the measurement there are significant density gradients and thus gradients of index of refraction in the test section, causing the phase shift of the sample beam.If there is a phase shift between sample and reference beams, the interference patterns are formed and it is possible to take the interferogram using the digital camera.
The interference fringes visible in the interferograms represent the lines of the constant refraction index.Evaluation of the flow parameters from the interferograms is based upon Lorentz-Lorenz equation, which relates the density ρ and the absolute index of refraction n of the transparent medium as follows where r is a specific refraction.It is experimentally verified that for the changes of density, pressure, temperature and even for some changes of state of the transparent media, specific refraction can be considered as a constant, [14].
Since the index of refraction for the air n .= 1, it is possible to rewrite the relation (4) into equation which can be simplified to Gladston-Dale relation where K is the Gladston-Dale constant.
From geometrical optics and considering assumption of negligible influence of beam curvature on the optical path length, the following relation can be derived ρ S is a known reference density, i is a number of interference fringe from the reference one and C is a constant involving wavelength of the used light, Gladston-Dale constant and width of the test section Knowing the density at some point of the flow field, the distribution of the density can be evaluated from interferograms using equation (7).As a reference density ρ S the stagnation density of the air was assumed in the entrance plane of the minichannel, i.e. the first interference fringe has a number i = 1.The example of an interferogram can be seen in Figure 6.

Results of the experiment
At first, the experimental investigation was carried out for a minichannel of the height h = 2 mm.Both optical and pneumatic measurements were performed for various pressure ratios π (ratio of the back pressure to stagnation pressure).
In Figure 7 the comparison of pressure ratio p i p 01 (ratio of the static pressure to stagnation pressure p 01 ) for pressure ratio π = 0.188 obtained using pneumatic and optical methods can be seen.It is necessary to be aware of the fact that these two methods should not provide exactly same results, because for the computation of the pressure from optical measurements the isoentropic relation between pressure and density p = p 01 ρ ρ 01 κ was used, and thus the value of static pressure is lower than the one obtained using pneumatic method.The pressure ratio π = 0.188 was chosen as a representative for an illustration.For all the cases when the flow is choked the plots of pressure are very similar.Another set of experimental data was obtained using only pneumatic methods for minichannel of the height 0.5, 1, 2, 4, and 5 mm and for the same pressure ratios π.The resulting plots of the pressure ratio distribution along the minichannel axis for various heights and pressure ratio π = 0.188 can be observed in Figure 8.As can be seen, the narrower the gap, the steeper course of the pressure ratio was measured.Near the minichannel exit the measured pressure ratio is almost the same for all heights when the flow is choked.It is of particular importance that for minichannel heights of 4 and 5 mm the static pressure measured by the pressure tap located at x = 10 mm is lower than the pressure measured by the next pressure tap downstream.This might be caused by the radius of the convergent part of the minichannel and following discontinuous change in curvature at the very beginning of the channel with constant cross-section.

Numerical results in comparison with experiment
The gap of height h = 2 mm and pressure ratio π = 0.188 was chosen as a representative.First laminar and fully turbulent simulations were performed.Results are presented in Figure 9 in form of the pressure ratio p i /p 01 (static pressure to stagnation pressure) along the minichannel axis.As can be seen from the graph, neither laminar nor fully turbulent aproach to this case is suitable and as stated in the introduction, transition from laminar to turbulent flow regime has to be simulated.In this case Figure 9 shows very good agreement of the experimental and numerical results, which was achieved by use of the γ − Re θt model of transition.Distribution of Mach number along the minichannel axis is displayed in Figure 10.The air flow behind the minichannel is also a subject of interest in this study.Good agreement of the numerical results with the experimental ones confirmes suitability of the used mathematical model and the chosen numerical method for the simulation.Since interference fringes are directly related to the density isolines (through relations ( 6), ( 7) and ( 8)), comparison of the interferograms and computed density contour lines in Figures 11 and 12 is possible.They show corresponding sequence of shock waves, which is formed in the outflow behind the minichannel.

Summary
Flow conditions in the narrow gaps classified as minichannels have been studied both experimentally and numerically.Modern developing discontinuous Galerkin finite element method was implemented for spatial discretization of both the system of Favre averaged compressible Navier-Stokes equations and the transport equations of turbulence and transitional model.Presented results for minichannel of height 2 mm and pressure ratio π = 0.188 show good agreement between experimental measurement and numerical simulation using the γ − Re θt transitional model.Preliminary simulations have been run for another pressure ratios and minichannel widths and will be the aim of the subsequent study, together with giving more precision to the measurements.
During the measurements some drawbacks of the measurement facility emerged and must be solved for further investigation.The most serious one is very laborious setting of the minichannel height and the parallelism of the walls.Even though the smallest feeler gauge used for the adjustment of the height had a thickness of 0.05 mm, the actual inaccuracies might be higher due to compliance of the whole structure.It turned out that the flow is very sensitive on the accuracy of the minichannel walls adjustment and its parallelism.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences, 2015

Fig. 1 .
Fig. 1.The measurement facility depicted in isometric projection (without the right wall).

Fig. 2 .Fig. 3 .
Fig. 2. Schematic view on the minichannel formed by two steel plates with marked main dimensions l, w and h (windows serving as a side walls are not depicted).experimental facility settling chamber pipe sliding valve vacuum storage

Fig. 4 .
Fig. 4. Scheme of the upper wall of the minichannel of the length 100 mm with plotted positions of pressure taps and a detailed view of the pressure tap.

Fig. 5 .
Fig. 5. Schematic depiction of the light beam passing through the Mach-Zehnder interferometer.

Fig. 6 .
Fig. 6.Interferogram of the flow inside a minichannel of the height 2 mm and pressure ratio π = 0.188 (the image is outstretched in the vertical direction).

Fig. 7 .Fig. 8 .
Fig.7.Ratio of the static pressure to stagnation pressure p 01 along the minichannel axis for pressure ratio π = 0.188, obtained using pneumatic and optical methods.

Fig. 9 .
Fig. 9. Distribution of static to stagnation pressure ratio along the minichannel axis for pressure ratio π = 0.188.

Fig. 10 .
Fig. 10.Distribution of Mach number along the minichannel axis for pressure ratio π = 0.188; speed of sound is achieved at about 85% of the channel length.

Fig. 11 .
Fig. 11.Interferogram of the flow at the end of the minichannel and in the free outflow, for minichannel of height of 2 mm and pressure ratio π = 0.188.