Discontinuous Galerkin Finite Element Investigation on the Fully-Compressible Navier – Stokes Equations for Microscale Shock-Channels

Microfluidics is a multidisciplinary area founding applications in several fields such as the aerospace industry. Microelectromechanical systems (MEMS) are mainly adopted for flow control, micropower generation and for life support and environmental control for space applications. Microflows are modeled relying on both a continuum and molecular approach. In this paper, the compressible Navier–Stokes (CNS) equations have been adopted to solve a two-dimensional unsteady flow for a viscous micro shock-channel problem. In microflows context, as for the most gas dynamics applications, the CNS equations are usually discretized in space using finite volume method (FVM). In the present paper, the PDEs are discretized with the nodal discontinuous Galerkin finite element method (DG–FEM) in order to understand how the method performs at microscale level for compressible flows. Validation is performed through a benchmark test problem for microscale applications. The error norms, order of accuracy and computational cost are investigated in a grid refinement study, showing a good agreement and increasing accuracy with reference data as the mesh is refined. The effects of different explicit Runge–Kutta schemes and of different time step sizes have also been studied. We found that the choice of the temporal scheme does not really affect the accuracy of the numerical results.


Introduction
It was 1959 when Richard Feynman gave his famous lecture at the meeting of the American Physical Society at Caltech called "There's Plenty of Room at the Bottom", where he proposed two challenges with a prize of $10,000 each: the first one was to design and build a tiny motor, while the second one was to write the entire Encyclopaedia Britannica on the head of a pin.
Nowadays, his speech is considered as the foundation of modern nanotechnology, since he highlighted the possibility to encode a number of pieces of information in very small spaces, hence producing small and compact devices [1].All those extremely small devices having characteristic length of less than 1 mm but more than 1 micron are called microelectromechanical systems (MEMS) and, as the name suggests, they combine both electrical and mechanical components [2].MEMS are small devices made of miniaturized structures, sensors, actuators and microeletronics and their components are between 1 and 100 micrometers in size.In recent years, several MEMS have been designed and developed, from small sensors to measure pressure, velocity and temperature, to micro-heat engines and micro-heat pumps and their numerical investigations are indispensable.
From a historical point-of-view, a pioneer experimental work on shock wave propagation in a low-pressure small-scale shock-tube was carried out by Duff (1959) [3], where a non-linear attenuation of the shock wave propagation for a certain diaphragm pressure ratio was observed.Other experimental works were performed by Roshko (1960) [4] and Mirels (1963Mirels ( , 1966) ) [5,6] confirming the strong attenuation of the shock wave and the acceleration of the contact surface, which propagates behind the shock wave in the classic shock-tube test case.The time interval between the shock wave and the contact surface measured at a certain point-which is also known as flow duration-rarefaction effects and thermal creeping were explained in depth.It is important to note that experimental and numerical studies on shock waves in different fields of engineering sciences attracted researchers over the past seventy years [3][4][5][6][7][8][9][10][11][12][13][14].However, researchers paid particular attention to microscale shock waves and rarefied gas dynamics recently, especially for aerospace applications.This is due to the fact that microengines are used in the development of aerospace propulsion systems, because of their reduced size and achievable high power density.One of the greatest difficulty in the design process of microengines is that the fast heat loss results in low efficiency of these microdevices.Therefore, researchers devoted attention to carrying out experimental and numerical works on shock wave propagation and formation in micro shock-tubes and channels for MEMS applications [15][16][17][18].
In the aerospace industry, microfluidics is becoming more and more popular having applications mainly in aerodynamics, micropropulsion, micropower generation and in life support and environmental control for space applications.For instance, MEMS can be adopted for flow control problems for both free and wall bounded shear layers flows.In 1998, Smith et al. [15] studied experimentally the control of separated flow on unconventional airfoils using synthetic jet actuators to create a "virtual aerodynamic shaping" of the airfoil in order to modify the airfoil characteristics.Microfluidics is also used through fluidic oscillators in order to produce high-frequency perturbations for example to decrease jet-cavity interaction tones [16].MEMS-based devices are adopted in the aerospace industry for the sake of turbulent boundary layer control.In fact, the small sizes of those systems (high density devices) allow to study near-wall flow structures [17].In space applications, micropropulsive devices are designed and developed for miniaturized satellites, mainly used for global positioning systems or to serve generic platforms [18].A detailed review on the application of microfluidics related devices in the aerospace industrial sector can be found in [18].
The flow behavior at those microscales is in general characterized by a granular nature for liquids and a rarefied behavior for gases; the walls "move", hence the classical no slip boundary conditions adopted in the macro regime fails.In agreement with [19], it is possible to classify the main differences among macrofluidics and microfluidics in the following list: noncontinuum effects, surface-dominated effects, low Reynolds number effects and multiscale and multiphysics effects.Furthermore, it is also observed that the diffusivity effects play an important role at this scales (see, e.g., [20]), especially when compared with the transport effects of the flow.Dealing with gases in micro devices, it is common practice to classify different flow regimes through the dimensionless Knudsen number Kn.Let λ be the mean-free path , which is the average distance traveled by a molecule between two consecutive collisions; denoting with the characteristic length of the generic problem considered (e.g., the hydraulic diameter for a channel flow problem), the Knudsen number is defined as Kn = λ .
Microfluidics can be modeled with two different approaches.The first one is the continuum model, and the flow is considered as a continuous and indivisible matter, while in the molecular model, the fluid is seen as a set of discrete particles.These models are valid in specific flow regimes determined by the Knudsen number and, when Kn increases, the validity of the continuum approach becomes questionable and the molecular approach should be adopted, as briefly sketched in Figure 1.
When the continuum approach is adopted, the fully-compressible Navier-Stokes (CNS) equations must be numerically solved.In the literature, this is usually performed adopting finite volume solvers.In the present work, the authors investigate how the discontinuous Galerkin finite element method (DG-FEM) performs applied to compressible flows at microscale levels in the slip flow regime (low Knudsen number).This method is selected because it takes advantages from the classical finite element method (FEM) and the finite volume method (FVM) since discontinuous polynomial functions are used and a numerical flux is defined among cells to reconstruct the solution.To verify the DG-FEM code, due to a lack of experimental data in microfluidics, the Zeitoun's test case [21] is adopted, which consists of a mini viscous shock channel problem numerically solved.In particular, they adopted the following models: the CNS equations in a FVM context, DSMC (Direct Simulation Monte Carlo) for the Boltzmann equation and the kinetic model BGKS (Bhatnagar-Gross-Krook with Shakhov equilibrium distribution function) model.In our work, the open source MATLAB code-developed by Hesthaven and Warburton [22]-has been adopted, modified and further improved.

Compressible Navier-Stokes Equations
Consider a generic domain Ω ⊂ IR d being d = 1, 2, 3 the dimension, provided with a sufficiently regular boundary ∂Ω ⊂ IR d−1 oriented by outward pointing normal unit vector n.On a two-dimensional (d = 2) cartesian reference system characterized by unit vectors i and j, the position vector is x = xi + yj.Consider a gas with vector velocity field u = ui + vj, density ρ, pressure p and total energy E. All the properties considered are both space and time dependent, e.g., u = u(x, y, t).The fully-compressible set of governing equations made of the continuity equation, Navier-Stokes momentum equations and energy conservation form a set of m partial differential equations which can be written in a vectorial form as ∂w ∂t In the equation above, w(x, t) = (ρ, ρu, ρv, E) T is the vector of conserved variables; f c (w) = (ρu, ρu 2 + p, ρuv, (E + p)u) T and g c (w) = (ρv, ρuv, ρv 2 + p, (E + p)v) T are the convective fluxes in the x and y directions, respectively; and f v (w, ∇ ⊗ w) = (0, τ xx , τ xy , τ xx u + τ xy v) T and g v (w, ∇ ⊗ w) = (0, τ xy , τ yy , τ xy u + τ yy v) T are the viscous fluxes in the x and y directions, respectively.The terms τ ij are the entries of the second-order viscous stress tensor τ τ τ.The latter is related to the velocity field according to the Navier-Stokes hypothesis for Newtonian, isotropic, viscous fluid through the following formulation: ] the strain rate tensor and I the identity matrix.The viscous stress tensor entries are The total energy is linked to the other fluid properties through the following equation of state (EOS) for a calorically ideal gas: E = p γ−1 + 1 2 ρ|u| 2 , where γ is the specific heat capacity ratio and Consider the compressible Navier-Stokes equations written in compact form (1) where the viscous fluxes are taken on the on the left hand side by ∂w ∂t if one defines F(w) as a m × d matrix having as columns the differences among the convective and viscous flux vectors, respectively, in the x and y direction Equation ( 5) can be expressed as In particular, w(x, t) :

Discontinuous Galerkin Finite Element Method (DG-FEM) Formulation
The physical domain is approximated by the computational domain Ω h which consists of an unstructured grid made of K geometry conforming non-overlapping elements D k , with k = 1, . . ., K. A non-negative integer N is introduced for each element k and let IP N be the space of polynomials of global degree less than or equal to N. The following discontinuous finite element approximation space is introduced [23]: being L 2 (Ω h ) the Hilbert space of square integrable functions on Ω h .Using DG-FEM, the vector of conserved variables w(x, t) is approximated by a function w h (x, t), which is the direct sum of K local polynomial solution w k h (x, t) by Analogously, one has which means that F F h = F(w h , ∇ ⊗ w h ).The local solution is expressed as a polynomial of order N through a nodal representation as being k i the multidimensional interpolating Lagrange polynomial defined by grid points x i on the element D k and N p the number of terms within the expansion which is related to the order of polynomial N through the relation N p = (N+1)(N+2)

2
. From this perspective, recalling Equation ( 7), the residual is formed as The residual can vanish requiring that it is orthogonal to all test functions φ h (x) ∈ V h on all the K grid elements Using the Gauss' theorem, it can be easily shown that the latter reduces to From the RHS of the last equation, one can observe that the solution at the element interfaces is multiply defined, thus, it is possible to refer to a solution F * h to be determined.Reconsidering the flux vectors of the matrix F, and considering that the normal vector is defined as n = nx i + ny j, one has which gives the following weak form: The numerical flux indicated with the superscript ' * ' is computed through the local Lax-Friedrich flux as where λ in general represents the local maximum of the directional flux Jacobian and an approximate local maximum linearized acoustic wave speed can be given [22] by Note that, even if this flux has a dissipative nature, hence strong shock wave in supersonic regime can have a smeared trend, it gives accurate results for subsonic and weakly supersonic flows.For a generic quantity, the superscripts "−" and "+" here indicate, respectively, an interior and exterior information, i.e., if the quantity is taken at the internal or external side of the face of the element considered.The symbols {{•}} and [[•]] are the average and the jump along a normal n, which are defined, for a generic vector v, as {{v}}

Temporal Integration Schemes
Considering the semi-discrete problem, written in the form of a system of ordinary differential equations (ODEs), the corresponding initial-value problem when initial conditions are given at time L(•) is the elliptic operator.Since the flow is strongly characterized by flow discontinuities, the strong stability-preserving Runge-Kutta (SSP-RK) schemes are adopted because they do not introduce spurious oscillations.Referring to Gottlieb et al. [24], the optimal second-order, two-stage and third-order three-stage SSP-RK schemes are expressed as 2 nd -order, two-stage SSP-RK : 3 rd -order, three-stage SSP-RK : Gottlieb et al. [24] showed that it is not possible to design a fourth-order, four-stage SSP-RK where all the coefficients are positive.The classical fourth-order four-stage explicit RK method (ERK4) might be adopted, however the main disadvantage of this approach is its high computational effort since for each time step, four arrays must be stored in the memory.A valid alternative to this method, is given by the low storage explicit Runge-Kutta (LSERK) scheme, firstly introduced in 1994 in [25].
The fourth-order LSERK is defined by (5) . ( The coefficients a i , b i and c i are listed in Table 1.As the formula above shows, different from the classical ERK4, in this case, only one additional storage level is required.However, the LSERK requires five stages instead of four.The time step size ∆t that ensures a stable solution is computed [22] as where a is the local speed of sound, which, using the ideal gas law, reads a = √ γRT = γ p ρ , while the geometrical factor ϑ is computed as ϑ = , where F scale is a matrix having dimension N f aces × K and its entries are the ratio of surface to volume Jacobian of face i on element k.From Equation (24), one can observe that, with very high order polynomials (N >> 1), this time step restriction becomes impracticable; furthermore, the time step decreases as the dynamic viscosity µ increases, hence for highly viscous fluids, this expression for the time step might be unfeasible.

Slope Limiting Procedure
Due to strong flow discontinuities, the solution might be affected by spurious unphysical oscillations, hence, slope limiters are added to the existing code in order to properly model and catch the large gradients in the flow field.In particular, van Albada type slope limiter suitable for DG-FEM, throughly described by Tu and Aliabadi in [26], are adopted.

Benchmark Test Problem with Its Initial and Boundary Conditions
The viscous shock wave propagation is studied in a microchannel characterized by characteristic length (hydraulic diameter) equal to H = 2.5 mm.The viscous shock channel problem of Zeitoun et al. [21] is characterized by a driver and a driven chamber, quantities referred to these states are denoted with subscripts 4 and 1 and summarized in Table 2.The driver chamber is characterized by a higher pressure and density.The gas used in both chambers is Argon (Ar) and the main fluid properties of this gas, considered in standard condition, are listed in Table 2c.A sketch of the microchannel in the cartesian reference system (x, y) is given in Figure 2. Geometric information, initial conditions and the main flow properties of the argon are given in Table 2.  Since the continuum approach is adopted, the rarefaction effects are usually taken into account imposing at wall the following conditions: • velocity slip boundary condition; • temperature jump boundary condition.
The small area where thermodynamic disequilibria occur is called Knudsen layer, having thickness of order of the mean free path λ.A generic form of the slip boundary condition is proposed by Maxwell.Let u slip be the fictitious velocity required to predict the velocity profile out of the layer, the slip velocity can be expressed as being u f the fluid velocity, n and s the normal and parallel directions to the wall, and σ u the tangential momentum accommodation coefficient which denotes the fractions of molecules absorbed by the walls due to the wall roughness, condensation and evaporation processes [27].For microchannels, accurate values of σ u are in the range 0.8-1.0[28].For the temperature jump [21], a condition is imposed by where T s is the temperature that must be computed at wall that takes into account the gas rarefied conditions, T wall is the reference wall temperature, Pr is the Prandtl number and σ T is the thermal accommodation coefficient.In the literature, different empirical and semi-analytical expressions are available for λ and they are based on the way the force exerted among molecules is defined.In this work, the inverse power law (IPL) model is used, firstly introduced in 1978 by Bird in [29].The model is based on a description of the mean free path based on the repulsive part of the force.It defines λ as where k 2 is a coefficient which varies according to the model taken into account.According to the Maxwell Molecules (MM) model, this constant is equal to Hence, if the temperature variation at walls is neglected, the slip boundary conditions becomes and the temperature jump The slip boundary condition and temperature jump at wall are conditions desired for high Knudsen number regimes, whereas rarefied condition of the gas becomes predominant.However, since the present work focuses on the investigation of low Knudsen number regimes, the no-slip boundary condition case is used as initial approach.
The shock channel problem consists of two chambers at high (on the left, denoted with number 4) and low (on the right, denoted with number 1) pressure separated by a diaphragm in a known position x d .When the diaphragm is instantaneously removed, i.e., when t > 0, due to the initial pressure difference, a combination of different wave patterns arises.The flow considered is originally at rest (u = 0) and the initial conditions of the problem are given by u(x, y, 0) = 0, (31) Numerical values of the quantities above are summarized in Table 2b.

Results and Discussion
To verify the MATLAB code and to validate the numerical results achieved, the benchmark test problem of Zeitoun et al. [21] on the investigation of viscous shock waves are considered, because this is one of the most frequently used benchmark problem for microscale applications.In their work, the viscous shock channel problem is solved at micro scales adopting three different approaches: compressible Navier-Stokes (CNS) equations with slip and temperature jump BCs using the CARBUR solver, the statistical Direct Simulation Monte Carlo (DSMC) method for the Boltzmann equation and the kinetic model Bhatnagar-Gross-Krook with the Shakhov equilibrium distribution function (BGKS).

Grid Convergence Study
The validation of the numerical results achieved with DG-FEM is performed through a grid convergence study using four grid levels.The grid levels are indicated with the index i, respectively, equal to 4, 3, 2 and 1.Let ∆x and ∆y be the mesh widths, respectively, in the x and y directions and N x and N y the number of grid points.The mesh widths in both directions have same length.A constant refinement ratio R = ∆x i+1 ∆x i = ∆y i+1 ∆y i among all grid levels equal to 2 is considered.The required data for the mesh refinement study are listed in Table 3, whereas the quantity h is the dimensionless grid spacing which is the ratio among the grid spacing of the i-th grid level considered and the grid spacing of the finer mesh, defined as h i = ∆x i ∆x 1 = ∆y i ∆y 1 , i = 1, . . ., 4.
Table 3. Grid levels adopted in the mesh refinement study.The simulations are performed setting the Knudsen number equal to 0.05 and the order of polynomials is kept equal to 1 for all the grid levels.The results are considered at the output time T f = 80 µs: before this time, the acoustic waves are propagating in the channel without considering the reflection at lateral walls.Figure 3a,b shows, respectively, the dimensionless density ρ/ρ 1 and temperature T/T 1 extracted at the centerline of the channel (y = H) plotted against the dimensionless x/H coordinate.Qualitatively speaking, referring to the density profile in Figure 3a, one can see that the accuracy of the numerical results achieved increases when the mesh used is finer; in particular, with the coarse and medium mesh, the profile is very diffusive yielding to an incorrect and inaccurate representation of the discontinuities in the flow field.In fact, the density in the rarefaction wave region is over predicted and, as a result, the positions of the contact wave and of the shock wave are imprecise.The real flow physics is matched when the fine and finer mesh are considered, since less numerical diffusion can be observed and, as a result, the density jumps are properly caught.Analogous considerations can be done for the dimensionless temperature profile shown in Figure 3b.Firstly, one can observe that the accuracy quickly increases as the mesh is refined, in fact, for instance, the results achieved adopting the coarse mesh do not match at all the jumps in the flow field observed by Zeitoun et al. [21].Furthermore, considering the finer grid level, it is observed that the position of the contact wave is properly achieved using DG-FEM, however the whole jump in temperature is slightly bigger than the one observed in the reference data (the relative error observed is approximately 10.6%).This produces a small under prediction of the shock wave position (relative error equal to 3.8%).

Mesh
The numerical results are validate also in terms of streamwise velocity u (Figure 3c) in x = 25H, which represents the position immediately before the shock wave (pre-shock state).The velocity is made dimensionless using the speed of sound in the driven chamber defined as a 1 = √ γRT 1 and plotted against the dimensionless y/H coordinate (for the first half of the channel's height).The velocity profile obtained with the coarse and medium mesh under predicts the outcomes given by Zeitoun et al. [21], while bigger values are achieved adopting the fine and finer mesh; however, when y/H ≈ 0.3, the profiles obtained overcome the reference data.The numerical results obtained for the stream wise velocity are quite accurate even if they are studied in the no slip fluid flow regime.However, as the Knudsen number increases, the slip condition becomes mandatory [19].To confirm that the numerical results achieved grid convergence, a simulation on an additional grid level 5-which is the finest mesh-is performed to compare the results between the finer and the finest mesh (see Figure 4).The refinement ratio is kept equal to 2, so the finest grid level is characterized by N x = 1536 and N y = 96 grid points in the x and y directions, respectively, and mesh widths ∆x = ∆y = 0.52 × 10 −1 mm. Figure 4 shows that the mesh further refinement does not improve significantly the already obtained accuracy of the numerical simulation results, which means that the further refinement of the mesh compared to the grid level 4 could increase the computational cost without significant further improvement of the achieved accuracy.For the sake of a quantitative analysis, the L 0 , L 1 and L 2 norms of the absolute error between the numerical results and the reference data given by Zeitoun et al. [21] are computed.Figure 5 shows the logarithmic plots of the L 0 (Figure 5a), L 1 (Figure 5b) and L 2 (Figure 5c) norms of the absolute error between the results achieved with DG-FEM and reference data given by [21] against the inverse of the dimensionless grid spacing h.Regarding the density and the temperature, the three plots clearly show that as the mesh is refined, the error drops in terms of all the norms considered.The same trend can be observed for the velocity profile, however, going from the fine to the finer grid level, the error increases after log(1/h) = 0.5.The reason is that the present investigation focuses on the low Knudsen number flow regime, where rarefaction effects start to become important but still not dominant.Rarefaction effects are taken into account in a continuum approach-i.e., using compressible Navier-Stokes equations-imposing wall slip boundary conditions producing hence a different velocity profile [19].The slip boundary condition becomes a mandatory requirement as the Knudsen number increases, which would yield to more physical and accurate results as confirmed by other authors in [19,21] when other continuum based numerical approaches were employed as well.This behavior is also met in the qualitative discussion above, since it is seen that the velocity profile achieved through the finer mesh overcomes the reference data when y/H ≈ 0.3.Furthermore, the lowest error norms are observed for the stream wise velocity profile extracted in x = 25H.For the sake of a complete analysis, the error norms are also summarized in Table 4.
Table 4. L 0 , L 1 and L 2 norms of the absolute error between the results achieved with DG-FEM and reference data in [21].Results presented for density (a), temperature (b) in the centerline of the microchannel, and stream wise velocity (c) in the cross section x = 25H.Results obtained at the final time T f = 80 µs with Kn = 0.05.

Coarse Medium
Fine Finer A convergence test is performed in order to understand if the formal order of accuracy matches (or not) the observed order of accuracy.Hence, within this approach, one can understand if the discretization error is reduced at the expected rate [30].The formal order of accuracy can be achieved from a truncation error analysis, and, in the FEM approach, it is equal to N. In particular, in this case, N = 1 since first-order polynomials are considered.The observed order of accuracy P is computed from the numerical outputs on systematically refined grids [30].The observed order of accuracy is based on the trend of the error.Consider two generic grid levels i and i + 1, being the i-th level the finer among them, and let e abs and e (i+1) abs be the absolute errors for these grid levels.The observed order of accuracy, based on the L p norms of the errors is defined by , with p = 0, 1, In this work, different grid levels are taken into account along with different observed orders of accuracy for each flow property and for each norms.Figure 6 shows three logarithmic plots of the observed order of accuracy adopting the L 0 , L 1 and L 2 norms of the absolute error between the results achieved with DG-FEM and reference data given by [21] against the dimensionless grid spacing h.The simulations are performed using the first-order polynomial representation N = 1 and, for a sake of clarity, the observed order of accuracy for each quantity, for each norm, is also listed in Table 5. Figure 6a shows that, regarding the density at the centerline of the microchannel, for the L 0 and L 2 norms, the observed order of accuracy increases as the mesh is refined, reaching the values 1.076 and 1.104, respectively, which are higher than the theoretical order N = 1.The same trend is not shared by the L 1 norm, which exhibits a maximum between the medium and fine grid level.Concerning the temperature profile (Figure 6b), all the values are below the theoretical order N = 1 and the observed order increases as the mesh is refined for the L 1 and L 2 norms, while the L 0 norm shows a maximum between the medium and fine mesh.The stream wise dimensionless velocity profile in Figure 6c shows a different trend, which means that the observed order of accuracy P decreases as the mesh is refined, as also confirmed by the previous discussions about Figures 3c and 5  Table 5. Observed order of accuracy P for density, temperature and u velocity based on L 0 , L 1 and L 2 norms of the absolute error between numerical solution obtained with DG-FEM and numerical data in [21].

Effect of Different Time Integration Schemes
The effect of different time integration schemes is investigated.The following explicit Runge-Kutta schemes are taken into account: second-order, two-stage SSP-RK, third-order, three-stage SSP-RK and fourth-order LSERK.The limiting procedure is applied for each stage of the methods.The investigation is performed adopting both the medium and fine mesh and Figure 7 shows the results obtained for density and streamwise velocity.Referring to Figure 7a,b, where the density profile is plotted, respectively, adopting the medium and fine mesh with the three RK schemes presented, no big differences are observed.On the right of the figures, some zooms are shown: a unique trend is not observed, with the medium mesh the accuracy of the second-order, two-stage SSP-RK is always between the other two schemes.When the fine mesh is used, the differences among the schemes become even smaller and the third-order, three-stage SSP-RK seems to hold an average trend among the other methods.Broadly speaking, the same trend can be observed for the stream wise velocity profiles in Figure 7c,d , where the medium and fine grids are used, respectively.When the medium mesh is used, the third-order, three-stage SSP-RK scheme gives the most accurate profile, while using the fine mesh the results achieved are very similar.In particular, the second-order, two-stage SSP-RK is more accurate in the first part of the profile and the fourth-order LSERK in the second.
Of course, one can see that, due to the small differences among the outputs obtained, a qualitative analysis cannot determine correctly which scheme yields to the most accurate results.Hence, as previously done, the L 0 , L 1 and L 2 norms of the absolute error are computed and listed in Table 6 for the medium and fine mesh.On the one hand, regarding the medium mesh, it is possible to observe that the third-order, three-stage SSP-RK scheme is the least accurate, since the lowest error norms are gotten.This behavior is met for all norms, for all the physical quantities considered.On the other hand, when the fine mesh is adopted, a unique trend is not observed: the fourth-order LSERK scheme is more accurate for temperature and density (except for the L 0 norm), while the second-order, two-stage SSP-RK is the most accurate for the stream wise velocity.To understand exactly how those time integration schemes perform, the relative difference among the previous error norms is compared.In particular, the relative difference among different temporal schemes using the L p norms are defined by where the superscripts I I , I I I and IV , respectively, indicate the error norms computed using the second-order, two-stage SSP-RK, third-order, three-stage SSP-RK and fourth-order LSERK.These quantities are collected, in percentage, in Table 7. From the table, one can see that the relative difference is bigger when the medium mesh is adopted and this trend is emphasized when the stream wise velocity is taken into account.However, generally speaking, when the fine mesh is considered, the relative differences has as order of magnitude 1%, which can be considered a negligible outcome.For this reason, it is possible to conclude that when a fine (or even a finer) mesh is considered, the choice of the temporal scheme do not really affect the accuracy of the numerical results.7. Relative difference among L 0 , L 1 and L 2 norms of the absolute error achieved using second-order, two-stage SSP-RK, third-order, three-stage SSP-RK and fourth-order LSERK scheme.Results obtained using the medium and fine mesh density and temperature in the microchannel's centerline and stream wise velocity in the cross section x = 25H.Results obtained at the final time T f = 80 µs with Kn = 0.05 and N = 1.

Effect of Different Time Step Sizes
The effect of different time step sizes ∆t on the accuracy of the numerical results is investigated.The time step size is computed using Equation ( 24) to get a stable solution according to [22].The following time step sizes are considered: 4∆t, 2∆t, ∆t, ∆t/2 and ∆t/4.The numerical simulations are performed using the fine mesh with second-order, two-stage SSP-RK scheme.Changing the time step size still gives stable and bounded solutions and relevant changes in terms of accuracy are not observed.The L 0 , L 1 and L 2 norms of the absolute error among reference data and numerical solutions are computed and it is observed that they share same precision up to the fourth digit.As for the temporal integration schemes study, the relative difference among the error norms is computed as: , for j = 1, . . ., 4 and p = 0, 1, 2. (37) The index j indicates a time step size level, in particular j = 1 corresponds to 4∆t and j = 5 to ∆t/4.Table 8 shows these relative differences and it can be observed that the order of magnitude, in percentage, is between 10 −7 and 10 −3 .For this reason, it is possible to conclude that time step sizes do not affect the accuracy of the numerical solution.
Table 8.Relative difference among L 0 , L 1 and L 2 norms of the absolute error achieved using different time step sizes.Results obtained with fine mesh and second-order, two-stage SSP-RK scheme and shown for density and temperature in the centerline of the microchannel and stream wise velocity in the cross section x = 25H.Results obtained at the final time T f = 80 µs with Kn = 0.05 and N = 1.

Conclusions
The two-dimensional unsteady fully-compressible Navier-Stokes equations for microfluidic problems are numerically solved adopting the nodal discontinuous Galerkin finite element method (DG-FEM) in space and different explicit Runge-Kutta (RK) schemes for the temporal integrations.The equations are solved at microscale level considering a miniaturized version of the shock-channel problem as test case.Unstructured meshes are generated and a MATLAB code developed by Hesthaven and Warburton [22] is adopted, modified and further improved.The numerical results are validated with reference data provided by Zeitoun et al. [21] where CNS equations in a FVM context, DSMC for the Boltzmann equation and the kinetic model BGKS (Bhatnagar-Gross-Krook with Shakhov equilibrium distribution function) models are adopted.A mesh refinement study is performed using four grid levels.The study showed that, as the mesh is refined, more accurate results are achieved.This trend is always confirmed for the density and temperature profiles, while for the streamwise velocity profile, the error slightly increases from the fine to the finer mesh.In fact, when the finer grid is adopted, it is observed that the method overestimates the velocity profile.A slip boundary condition implementation could improve this result.The observed order of accuracy based on different error norms for different flow properties is also computed and it is shown that, when the first-order of polynomial N = 1 is adopted (theoretical order of one), the following results are achieved: for the density and temperature, as the mesh is refined, the observed order of accuracy increases reaching (and overcoming for the density) the theoretical order; a different trend is achieved for the velocity profile, i.e., the order is higher as the mesh is coarser.The reason behind this behavior is the same as the previous one for the error norms.An additional grid level is then considered, and it is seen that the mesh further refinement does not improve significantly the already obtained accuracy.In addition to this, the effect of different temporal schemes (explicit Runge-Kutta) is investigated considering second-order, two-stage SSP-RK; third-order, three-stage SSP-RK; and fourth-order LSERK.A qualitative and quantitative analysis is done computing error norms and their relative differences, however no big differences among the schemes are observed as the mesh is refined and it is concluded that the choice of the temporal scheme does not really affect the accuracy of the numerical results, especially when fine meshes are used.Finally, different time step sizes have also been considered, and the solution remains stable and the accuracy unaffected.

Figure 1 .
Figure 1.Different flow regimes in function of the Knudsen number Kn.

Figure 2 .
Figure 2. The sketch of the viscous shock-channel of Zeitoun et al. [21] for microfluidic applications.

Figure 3 .
Figure 3. Dimensionless density (a); temperature (b) profiles in the centerline of the microchannel; and stream wise velocity u (c) in the cross section x = 25H using four grid levels at the final time T f = 80 µs with Kn = 0.05 in comparison with reference data of Zeitoun et al. [21].

Figure 4 .
Figure 4. Dimensionless density (a); temperature profiles (b) in the centerline of the microchannel; and stream wise velocity u (c) in the cross section x = 25H on the finer and the finest grid levels at the final time T f = 80 µs with Kn = 0.05.

Figure 5 .
Figure 5. L 0 , L 1 and L 2 norms of the absolute error between the results achieved with DG-FEM and reference data in [21] against the inverse of the dimensionless grid spacing h.Results presented for density and temperature in the centerline of the microchannel and stream wise velocity in the cross section x = 25H.Results obtained at the final time T f = 80 µs with Kn = 0.05. .

Figure 6 .
Figure 6.Logarithmic plots of the observed order of accuracy P using the L 0 , L 1 and L 2 norms of the absolute error between the results achieved with DG-FEM and reference data in [21] against the dimensionless grid spacing h.Results presented for: density (a); temperature (b) in the centerline of the microchannel; and streamwise velocity (c) in the cross section x = 25H.Results obtained at the final time T f = 80 µs with Kn = 0.05.

Figure 7 .
Figure 7. Dimensionless density profile in the centerline of the microchannel (a,b); and streamwise velocity u in the cross section x = 25H (c,d) using the: medium (a,c); and fine (b,d) mesh at the final time T f = 80 µs with Kn = 0.05 and N = 1.Comparison between different explicit Runge-Kutta schemes: second-order, two-stage SSP-RK, third-order, three-stage SSP-RK, and fourth-order LSERK.Validation with reference data of Zeitoun et al. [21].

Table 1 .
Coefficients a i , b i and c i used for the low storage five-stage fourth-order explicit Runge-Kutta method.

Table 2 .
Zeitoun's test case: geometric information and output time for the simulation (a); flow properties in the driver and driven chambers (b); and Argon properties in standard condition (c).