A study of the diffusive properties of a modified compressible Navier-Stokes model

The aim of this study is to provide further validation for the weakly well-posed modified compressible Navier-Stokes system proposed in Svärd (Phys A 506:350–375, 2018) when applied to ideal gases. We do so by considering sound attenuation, both theoretically and numerically for argon and oxygen, and make comparisons with experimental values in the literature. Furthermore, we compute shock profiles for argon and nitrogen, and compare them with experiments in the literature. Our numerical simulations have revealed problems when using experimental attenuation data, as presented in the available literature, for validation and determination of diffusion coefficients. However, comparisons with the shock data, suggest that the modified system may benefit from an additional heat diffusive term. In view of these and previously published validation tests, the model proposed in Svärd (Phys A 506:350–375, 2018) is equally accurate as the standard compressible system. However, with more complete experimental information for the attenuation case at hand, it might be possible to further improve the accuracy by more precise determination of the diffusion coefficients. We propose a tentative adjustment of the model that may be tested/validated, if more detailed experimental information becomes available.


Introduction
The compressible Navier-Stokes equations, for a heat conductive and viscous compressible ideal gas, take the form, where x T = (x, y, z) and Ω is the spatial domain. E = 1 2 |v| 2 + e is the total energy; e = c v T is the internal energy and T the temperature; e = p −1 , where and p are the density and pressure. v T = (u, v, w) are the Cartesian velocity components. Furthermore, the adiabatic exponent = c p ∕c v , R is the gas constant and the heat conductivity. The stress tensor is given by: ( ) ij = ij = − 2 3 v k,k ij + (v i,j + v j,i ) + v k,k ij . > 0 is the dynamic viscosity and since we are focusing on ideal gases, we assume that the bulk viscosity = 0.
The Navier-Stokes Eq. (1) have been the foundation of fluid dynamics for almost two centuries and since the 1950's there has been an ever increasing effort to design numerical methods that solve these equations. For smooth solutions, linear theory has (1) x ∈ Ω, t ∈ (0, T], t (E) + ∇ ⋅ (Ev + pv) = ∇ ⋅ v + ∇ ⋅ ( ∇T), p = RT, ideal gas law, been succesful and ever more accurate methods have been proposed. However, many phenomena are fundamentally non-linear and for those the endeavour to design numerical schemes has been severely hampered by the lack of well-posedness results for the Navier-Stokes equations; it is not known what properties a numerical scheme should enforce on its solutions. For example, it is shown in [11] that common state-of-the-art schemes for (1) may generate stable but erroneous solutions. In short, there is a dire need for a proof of well-posedness to advance numerical simulations further. Mathematical and numerical challenges aside, the standard Navier-Stokes equations have recently been questioned for physical reasons and various modifications of them have been proposed. Perhaps the most well-known is the modification proposed by Brenner, [3,4]. Quite remarkably this system was subsequently shown to admit weak solutions in [8]. More recently, Reddy et al. derived a class of modifications that was shown to produce more accurate results than the standard version (1) with respect to Rayleigh-Brillouin scattering experiments, see [15].
Herein, we will consider the modification proposed in [19], where the following line of arguments were put forth: The standard derivation of (1) is carried out for a mass element in a Lagrangian reference frame, which allows the stress tensor to act on the mass element. This is analogous to the equations for deformable solids. (See also the original derivation [18], where it is clear that this is the inspiration for the introduction of the viscous stress tensor.) However, in a solid, mass elements are natural as molecules are essentially fixed in relation to each other and the forces modelled by the stress tensor are intermolecular, caused by electromagnetism. In a gas, however, molecules mix or diffuse since they travel a finite distance before they hit another molecule. That is they are not locked in fixed positions in relation to other molecules but there is a random motion superposed on their mean velocities. It is this diffusive mechanism that gives rise to the macroscopic viscosity. Hence, and as pointed out in [19], diffusion is the primary property in a gas and viscosity is a secondary property. Since molecules that mix by diffusion have to transfer their mass, momentum and energy, and not just a subset of these three, the following alternative system, modelling diffusion, was proposed in [19]: where is a diffusion coefficient and a heat diffusion coefficient. (The notation is the same as in (1), apart from the diffusive coefficients, and both systems must also be augmented with boundary conditions. See [19] and [20]).
The diffusion coefficient is proposed to take the form where 0 , 1 are constants. Clearly, these must be chosen to make the model accurate. It was proposed in [19], that 0 is taken as the dynamic viscosity. This results in a −1 dependence on , which is consistent with a diffusive coefficient according to kinetic theory.
Taking 0 as the dynamic viscosity (and 1 = = 0 ) implies that (2) reduces to the incompressible Navier-Stokes equations for = constant and (2) would thus exactly reproduce the Blasius layer in the incompressible limit. Hence this choice suggests that, 1 is so small that solutions remain within the accuracy limits of the experiments that have validated the Blasius layer. Furthermore, the system (2) with = 1 = 0 and = 0 ∕ where 0 is the dynamic viscosity, has been validated in the compressible regime for a wide range of Mach and Reynolds numbers, in [6] for NACA0012 and in [16,17] for cylinders. A remarkable agreement with the standard system (1) was demonstrated in all cases. These results also suggest that 1 has to be orders of magnitude smaller than 0 .
However, 1 can not be zero. If 1 = 0 , the diffusion would effectively dissappear as → ∞ . Although, it is physically reasonable that the diffusion is smaller in a dense gas due to a shorter mean free path, it must not be allowed to decrease to zero. If the diffusion does vanish, the model cannot prevent infinite growth of the density, which is unphysical. The term, 1 ( 1 << 1 ), in (3) prevents this from happening (see [20]), while under normal conditions is still essentially ∼ 0 ∕ . Finally, in [20] the system (2) with 0 < 1 << 0 and = 1 4 r T 3 was considered. Although the latter term could be interpreted as a radiation diffusion, it is mainly a technical assumption that is necessary to control pressure and r > 0 should be very small. In stark contrast to the standard system (1), it was shown in [20] that, with these coefficients, (2) admits weak solutions.
The purpose of this paper is to further validate (2). To this end, we consider sound attenuation. We derive the attenuation coefficents for linear sound waves and carry out numerical simulations for the full non-linear systems for both argon and oxygen. We compare the results with experimental results available in literature. Next, we simulate argon and nitrogen shocks numerically using both systems and compare with experimental results available in the literature. Finally, we discuss our findings.

Remark
We emphasize that since (2) has as many physical parameters as (1), we do not expect (2) to be systematically more accurate than (1). The difference between the two systems lies in the weak well-posedness of (2) (see [20]), which is a necessary requirement for truly predictive simulations.

Sound attenuation
A sound wave is an irrotational and adiabatic (isentropic) wave of sufficiently small amplitude such that non-linear effects are negligible. Diffusive processes absorb energy resulting in a decaying amplitude of sound waves. In [13], the absorption coefficient for sound waves governed by the Navier-Stokes equations is derived and we refer to it for more details. Briefly, the attenuation is modelled as the reduction of energy, E, of the sound wave, i.e. the kinetic and internal energy of the deviation from the equilibrium. Since the wave is isentropic, one arrives at the following relation for small amplitude waves, where S = c v log(p∕ ) is the specific entropy. By assuming a velocity wave of the form, u 0 cos(kx − t) , and by using linearised thermodynamic relations (see [13]), one arrives at, which is the mean value of the energy dissipation. The absorption coefficient Γ is obtained by dividing |(Ē) t | by 2cE SW where E SW = 1 2 u 2 0 L is the integral of the total energy of the sound wave, L is the domain length and c the speed of sound. The amplitude of the sound wave decays as exp(−2Γt) . (In [13] the decay in x rather than t is computed. They differ by a scaling with the sound speed). For (1), one obtains f is the frequency, the wave number and = 0 the bulk viscosity.
Turning to the modified system, the same exercise results in, By taking = 0 ∕ in (3) and = 0 , (5) becomes Using the thermodynamic relations R = c p − c v and = c p ∕c v , it is straighforward to deduce that, For an ideal gas with = 5∕3 and Pr = 2∕3 (recall that noble gases are close to being ideal), one obtains, This difference between the two models was first noted in [14]. 1

Numerical experiments
Although the linear analysis used to obtain (4) and (6) is asymptotically valid for small amplitude waves, we have chosen to compute the attenuation numerically. This allows us to quantify the amplitude for which waves actually behave linearly. The sound waves have to be integrated for a long time, which calls for a method with very high accuracy in both space and time to keep numerical errors at a minimum without the need to use excessive computational resources. Since the problem is periodic and we intend to study linear waves, we opt for a Fourier spectral method [2] The initial conditions for both argon and oxygen are, where A is the amplitude. The domain length is L = 1 . Unless otherwise stated, the domain is discretized with 12 points, which leads to an accuracy of at least 4 significant digits in the attenuation coefficients. (See also remark below.) The simulations were run for 10 4 periods for p 0 = 10 5 and 10 3 periods for p 0 = 10 3 . The attenuation of the kinetic energy We ran the simulations with the amplitude, A = 10 −6 . In the Tables 1 and 2 the attenuation coefficients at the time intervals 10-40%, 40-70% and 70-100% are listed. (The first 10% of the simulations are omitted on account of initial transients.) We observe that the numerical coefficients are close to the theoretical values (7) for both oxygen and argon. The simulation with higher background pressure appears slightly more accurate than the lower.
We conclude that the attenuation coefficients obtained from the simulations of the full non-linear models match the coefficients obtained with linear theory, but only for extremely small amplitude waves. Hence, the theoretical discrepancy (8), which is only valid asymptotically for small amplitudes, has been formally confirmed. Furthermore, Greenspan's experiments for argon (See Fig. 2 in [10]) shows very good agreement with the standard model (1) in the continuum regime (corresponding to r ∈ (7, 70) in [10]).

Remark
The attenuation coefficients reported for the last third in Tables 1 and 2 are expected to be the most accurate since they have the smallest amplitude. They match the theoretical coefficient with 4 significant digits, which implies that the Fourier spectral scheme with 12 points in space is at least that accurate.
In view of Greenspan's measurements, it looks like the model (2) would need to be adjusted. The data in [10], match the Navier-Stokes well but there is no discussion on data analysis or experimental errors. Neither is the complete experimental setup given. The amplitudes of the waves and the distance between the measuring points are not mentioned. Similar experiments were run in [5] at 300K with ambient pressure. Unlike Greenspan's, their data points are somewhat scattered, but they do favour the standard model. Furthermore, the waves in Tables 1 and 2 behave according to linear theory but their amplitudes are minuscule. (Well below the auditory threshold.) If instead, we run the larger, but still very small, initial amplitude A = 5 ⋅ 10 −5 (still below the auditory threshold), the results are very different. (Here, we use 61 grid points in space to ensure that the distortion is not caused by numerical errors.) In Table 3, the attenuation coefficients for argon are listed. Their decays are depicted in Fig. 1 along with the shape of the wave at the end time ( T = 10.4).
A number of observations can be made: The computed attenuation coefficients are all too large, which is caused by weakly non-linear effects (the deformed shape in Fig. 1b) that increase the dissipation. Note that initially the wave is sinusoidal and the attenuation is closer to the theoretical and, as it deforms, the attenuation increases. This effect is not accounted for in (4) and (6) where linearity is assumed. In the last third, the actual (computed) attenuation coefficient of the modified system (2) is closest to the theoretical value for the standard system (1). (See also Fig. 1a). If this is the actual situation in the experiments we use as reference, then the experiments could equally well point to the modified system (2) as being the more accurate model! Our simulations clearly demonstrate that only extremely small amplitude waves behave linearly. The difficulty in measuring small amplitudes was pointed out in [7]. Thus, there is a risk that experiments overestimate attenuation due to non-linear effects. Indeed, [7] even recorded a growth of the amplitude in some cases showing that it is challenging to measure amplitudes accurately.
Furthermore, it is easier to measure attenuation for high frequencies and Greenspan used 11 MHz [10]. On the other hand, in [5] it is remarked that high frequencies are excessively attenuated due to boundary layers and turbulence, and they therefore carry out experiments at 1-3 MHz. These observations also suggest that Greenspan's experiments may model a situation similar to Fig. 1, where the attenuation is over-estimated due to non-linearities.
Our point here is not to discredit the experiments; we merely note that there are uncertainties when comparing models to experiments since there is a possibility that the attenuation measured in experiments is somewhat exaggerated. Had we known the  [1] amplitudes and travel times in the experiments, we could have assessed the non-linear effects and more decisively conclude whether or not (2) needs adjustments. Given that the current attenuation data does not establish that the waves are truly linear, one can not claim that (1) provides a better match than (2). Nevertheless, if we assume that attenuation data favour (1), then (2) lacks some diffusive mechanism. Otherwise (1) is too diffusive.

Remark
We also note in passing that the historical validation of (1) involves a certain freedom to slightly change Pr, (T) , and . For instance, we use both Pr = 0.68 and Pr = 0.661 for argon in this paper. These adjustments may change the outcomes by a few per cent, such that it is not exactly the same version that passes all validation tests as accurately. Thus, the difference to (2) also varies somewhat. The goal for both models should be a unique set of physical parameters, , and , respectively, for a given gas, and not different sets for different problems.
Finally, for oxygen we note that the theoretical, as well as the numerical, difference between the two models is about 6% , which is less than for monatomic gases. Opposite to the monatomic case, the modified model is more diffusive than the standard model. In [7], attenuation experiments for oxygen were performed and it is quite clear that the experimental errors exceed the model differences. However, we note that the high-frequency data is consistently higher than the theoretical value for (1). Hence, we can not use this case to determine the diffusive coefficients with greater certainty. We also stress that the data does not rule out that an extra diffusion term is necessary in (2). If so, then (1) also lacks a diffusive mechanism since (2) is already more diffusive.

Argon shocks
There are several articles on shock-tube experiments in the literature. Here, we use [1], where we have digitized data from Fig. 1 for Ma U = 1.55 and Ma U = 2.05 shocks. Since data is given in nondimensional form, we follow [9] and use the same initial data and physical constants.
For Argon, we have = 5∕3 . The upstream conditions (subscript U) for temperature, viscosity, sound speed, Mach number and density are: such that u U = Ma U . Since the sound speed c = √ RT and p = RT , we obtain R = −1 , the upstream pressure p U = 1 and the upstream mean free path The downstream conditions (subscript D) are given by the Rankine-Hugoniot conditions: The simulations are initialized with the up-and down stream states and a jump in between them at x = 0 . In the figures below, the down stream state, U , is to the left and the up stream state D to the right. Further, in [1] the density is non-dimensionalized with U and the x-axis is adjusted such that = 0.5 (nondimensionalized) is located at x = 0 . Hence, to compare results we translate the computational X-axis to x = (X − a)∕ U , where a is the appropriate constant that makes (x = 0) = 0.5.
In the Navier-Stokes model, we use the stress tensor with the bulk viscosity = 0 for monatomic gases. Furthermore, we use the same powerlaw model for the dynamic viscosity as in [9] as it gives a good match. In this non-dimensionalization, it is, for argon. (We have also tested Sutherland's formula but we do not report the results here as they are similar and do not change any conclusions.) Moreover, the Prandtl number is, Pr = 0.68 ; the heat capacity at constant pressure is, c p = ( − 1) −1 ; the heat conductivity is = c p Pr . For (2), we use 0 = (T) given by (12) in (3).

Numerical results
Since the problem is not periodic, we can not use the Fourier spectral method. Furthermore, the shock solutions are time independent, removing the need for an accurate time integration scheme.
Hence, a standard node-centred central finitevolume scheme without any artificial diffusion is sufficiently accurate on reasonably sized grids. At the boundaries, all variables are set by injecting the Rankine-Hugoniot values. The simulations are run, using the Euler-forward method in time, to steady state on a domain length of 40 U discretized with N = 801 points in space, such that the grid spacing is h ≈ 40∕N (since U ≈ 1).
The temporal residual (for density) is ≲ 10 −9 for the modified model (2) and ≲ 10 −8 for the Navier-Stokes model (1) (despite running more than twice as many steps in time). In all cases, the numerical errors are negligible (less than 1%).
We remark that the modified equations converge much faster to steady state than the standard model. This is not due to the modified system inherently being more diffusive. In view of the sound attenuation results, it is in fact less diffusive for argon. The reason is the structure of the diffusive terms, that makes the system less susceptible to amplify numerical errors. (See also [19,20]).
We compute shock waves with Ma U = 1.55 and Ma U = 2.05 with experimental data from [1]. The results are displayed in Fig. 2. Throughout, the standard model (1) is denoted as NS; the modified model (2) as modNS and experimental values with diamonds.
In view of Fig. 2, we note that the two models are equal downstream of the shock but the standard model is more accurate upstream.

Tentative update of the modified model
Despite uncertainties pertaining to experimental data, both the attenuation and shock cases seem to indicate that (2) may lack some diffusive mechanism. Hence, we investigate the possibility of adding diffusive terms. Any new terms are subject to a number of constraints: Remark The third constraint is included under the assumption that the attenuation data corroborating (1) are indeed truly linear waves; experimental data would otherwise be preferable to optimize against.
A naive rescaling of by some constant factor , such that = 0 + 1 , is not a viable option since it would make the model inconsistent with the incompressible Navier-Stokes and the Blasius solution. Furthermore, it was shown in [6,16] (with is not as accurate as = 1 that we have considered here. To gain some insight into the possibly missing diffusive mechanism, we return to the derivation of (2). We assume that we have field variables at hand. (The step from a particle cloud to field variables is discussed at length in [19]). We imagine an Eulerian box in the flow, with side length L. (See Fig. 3.) L   Fig. 3 Eulerian volume has to be significantly larger than the mean free path ("lambda" in the figure). The model is derived by using conservation arguments (for mass, momentum and energy) for the means inside the box, much like the derivation of the Euler equations. Indeed, the inviscid part is identical to the Euler equations and the diffusive terms are simply modelling particles crossing the volume boundary. As they do they bring their mass, momentum and (total) energy. That is why there is only one diffusion coefficient . The spatial arbitrariness of the box, leads to the set of PDEs (2).

Remark
The single diffusion coefficient can easily be verified using kinetic theory with a Maxwellian velocity distribution. See e.g. [12] where the dynamic viscosity is derived. The same principles can be used to derive the kinetic energy diffusion and the heat diffusion. For the latter, the coefficient is not the same as in (1). Instead, it leads to the heat diffusive term in (2). Finally, the same arguments can be applied to the mass of the molecules leading to the mass diffusion appearing in (2). All variables diffuse at the same rate since the diffusion coefficient is directly linked to the transport of molecules across a surface.
Since (2) is derived from the fundamental conservation principles, only the flow across the Eulerian box boundaries affect the mean values of the conservative variables , v and E inside the box. However, since diffusion takes place on the length scale of the mean-free path, there will be relaxation in the interior of the box due to collisions. Although, such collisions do not affect the conserved total energy, or any other conserved property, they might affect the balance between kinetic and internal energy, which in turn affect the dynamics of the system. Thus, the total energy may be subject to additional diffusion.
The only entropy consistent option, that does not directly affect the momentum and density equations, is to add a Fourier heat diffusion to the energy equation in (2). If we choose 0 such that the sound attenuation of (2) exactly matches (1) when r = 0 , we arrive at the heat diffusion coefficient This choice satisfies all the listed constraints.

Revisiting the validation cases
Verifying that (14) in (2), with r = 0 , leads to exactly the same attenuation coefficient as (1) is straightforward. Furthermore, we have recomputed the shock solutions with the extra diffusion term. The results are displayed in Fig. 4. The differences between (2) and (1) are now very small. The standard system is still slightly more diffusive downstream and slightly less diffusive upstream of the shock. Finally, we provide an example of the effect of the second diffusion coefficient, 1 . We choose 1 = 0.01 0 . The results are displayed in Fig. 5. The match with experimental data is somewhat improved but the effect is modest. Determining 1 with greater accuracy would require very accurate experimental data for intermediate Knudsen numbers. The Knudsen number in the shock layer is large and neither (1) nor (2) can be expected to be very accurate.

Diatomic gases
As mentioned above, the attenuation coefficient for (2) is larger than that of (1) for diatomic gases. However, experimental data is not accurate enough to favour one model over the other. Furthermore, we have run a Ma U = 1.53 nitrogen shock and compared the results with digitized experimental data provided in [1]. We have used (12) as the same exponent was claimed to be the most accurate also for nitrogen in [1]. Here, we can not choose the extra diffusion to match (1), since (2) is already slightly more diffusive. However, it may not be unreasonable that the model is updated in the same way, i.e., as (14) (with Pr = 0.72 and = 7∕5 for nitrogen). Hence, we have run the nitrogen shock using (2), with and without the extra diffusion given in (14), and the standard model (1). The results are shown in Fig. 6.
Without the extra diffusion the two systems are equally diffusive downstream and (2)  more diffusive everywhere. However, all the computed shock profiles are still too sharp in view of the experimental data.

Discussion
To further validate the modified Navier-Stokes model (2), we have studied its diffusive properties for sound attenuation in argon and oxygen, as well as Mach 1.55 and 2.05 shock waves in argon and a Mach 1.53 shock wave in nitrogen.
For the sound attenuation case, we derived the theoretical coefficients for (2) and (1). (Originally done in [14]). For argon, the attenuation of (2) is smaller and for oxygen it is larger. We also verified the attenuation rates numerically for both models. The good match between (1) and experimental data for argon (2) may be interpreted as (2) lacking diffusion. (For oxygen, the available experimental data is inconclusive with respect to the model differences.) However, our simulations also demonstrated that the linear attenuation only occurs at extremely small amplitudes, which makes it unclear if the experiments that display the best match with (1) truly are in the linear regime. To be able to use attenuation data to decide in (2), we would need, apart from the background state and frequency, also the wave amplitudes and travel times between measurements.
Turning to the shock simulations, both models predict too sharp shock profiles. For argon, (2) produces slightly sharper profiles than (1). Hence, we tentatively proposed to modify in (2) to provide further diffusion. The proposed adjustments imply that (2) and (1) attenuate sound waves at the same rate and their argon shock profiles become very similar. Furthermore, we demonstrated the effect of a small secondary diffusion ( 1 = 0.01 0 ), which, as expected, improved the profiles slightly. Nevertheless, we stress that the modification of requires further experimental justification as the sound attenuation experiments may have overpredicted the attenuation and shocks are anyway outside the validity range of (2) (and (1)).
Finally, we ran a Mach 1.53 nitrogen shock. First, with the original system (2), which performed slightly better than (1), and then with the extra heat diffusion term, which resulted in a small improvement.
In summary, for noble gases the evidence provide some support, although far from conclusive, for the addition of an extra heat diffusion term. However, its impact for real flow applications is probably minuscule as previous validation of (2) suggest, [6,16]. For diatomic gases, the experimental evidence in the attenuation case is inconclusive and it is unclear if (2) should be augmented with an extra Fourier diffusion, and if so, what size it should have. However, we recorded a slight improvement for the nitrogen shock with the extra heat diffusion term.
We conclude that, even without the modification suggested herein, the two models are, for all practical purposes [6,16], about as accurate models of compressible flows. However, (2) is easier to code, computationally less expensive and relaxes significantly faster to steady state solutions. But the single most profound advantage of (2) over (1) is its weak wellposedness, which is the bedrock of reliable and predictive simulations.
Funding Open access funding provided by University of Bergen (incl Haukeland University Hospital).

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, 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 article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article'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. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.