# On a hybrid continuum-kinetic model for complex fluids

A. Chertock Department of Mathematics, North Carolina State University, USA;
<EMAIL_ADDRESS>P. Degond Institut de Mathématiques de Toulouse;
UMR5219; Université de Toulouse; CNRS; UPS; F-31062 Toulouse Cedex 9, France;
<EMAIL_ADDRESS>G. Dimarco Department of Mathematics and
Computer Science $\&$ Center for Modeling, Computing and Statistics of
University of Ferrara, Italy<EMAIL_ADDRESS>M. Lukáčová-Medvid’ová
Institute of Mathematics, University of Mainz, Germany<EMAIL_ADDRESS>A. Ruhi Institute of Mathematics, University of Mainz, Germany;
<EMAIL_ADDRESS>
###### Abstract

In the present work, we first introduce a general framework for modelling
complex multiscale fluids and then focus on the derivation and analysis of a
new hybrid continuum-kinetic model. In particular, we combine conservation of
mass and momentum for an isentropic macroscopic model with a kinetic
representation of the microscopic behavior. After introducing a small scale of
interest, we compute the complex stress tensor by means of the Irving–Kirkwood
formula. The latter requires an expansion of the kinetic distribution around
an equilibrium state and a successive homogenization over the fast in time and
small in space scale dynamics. For a new hybrid continuum-kinetic model the
results of linear stability analysis indicate a conditional stability in the
relevant low speed regimes and linear instability for high speed regimes for
higher modes. Extensive numerical experiments confirm that the proposed
multiscale model can reflect new phenomena of complex fluids not being present
in standard Newtonian fluids. Consequently, the proposed general technique can
be successfully used to derive new interesting systems combining the macro and
micro structure of a given physical problem.

Keywords: multiscale simulations; hybrid method; kinetic equations;
homogenization; scale separation; Newtonian and non-Newtonian flows; fluid
dynamics; complex fluids

MSC: 76Nxx, 82C40, 76A05, 76M25

##  1 Introduction

Many important fluid flow problems are entirely multiscale: microscopic
processes strongly influence macroscopic behavior of the fluid and need to be
taken into account in order to accurately describe fluid dynamics. Typical
examples are granular [33] and high-speed rarefied flows [6], the plastic
deformation in materials [31], the viscoelastic [29] and biological type of
fluids [37]. For this reason, in the last few decades, there has been a huge
interest in both modeling and numerical simulations of problems associated
with multilevel physical models, which are able to incorporate multiscale
effects in different ways.

For the Newtonian fluids there have been many rigorous theoretical studies of
hydrodynamic limits and the relationship between microscopic molecular
dynamics and/or mesoscopic kinetic models of the Boltzmann type with
macroscopic models such as the compressible Euler or Navier-Stokes equations,
see, e.g., [2, 10, 22, 21, 28] and the references therein. On the other hand,
for complex fluids theoretical understanding is certainly less developed, and
more research is needed. More precisely, compared to standard fluids, the
challenge associated with complex fluids lies in an accurate determination of
rheological relations that are typically obtained from physical or
computational experiments. Consequently, in many situations, e.g., soft
matters or colloid-polymer mixtures, their full analytical description is not
available.

In order to take small scale effects into account we can apply either
mesoscopic kinetic models or directly microscopic models, such as molecular
dynamics or dissipative particle dynamics, to reconstruct time evolution of
macroscopic quantities. However, as is well known, an obvious drawback of
meso- and microscopic descriptions, despite their higher accuracy, is
sometimes prohibitively high computational costs, which limit their direct
application in many practical situations. To overcome this disadvantage and
make large-scale simulations possible the so-called hybrid multiscale methods
have been developed in the literature. The latter combine the advantages of
both descriptions: the accuracy of microscopic models to be able to consider
complex rheological relationships with the efficiency of macroscopic models
based on using classical conservation laws.

A prototype of hybrid multiscale methods is the heterogeneous multiscale
method proposed by E, Enquist et al. [16, 14, 15, 17, 19, 34], see also [8,
38, 36, 39, 40] for its application for complex polymeric fluids. In this
context let us also refer to triple-decker atomistic-mesoscopic-continuum
method [20], the seamless multiscale methods [35, 18], the equation-free
multiscale methods [25, 26] or the internal-flow multiscale method [4, 5]. In
[27] a overview of multiscale flow simulations using particles is presented.
For classical gas dynamics, similar ideas were employed, for instance, in [11,
9, 7, 32].

However, such general methods can have several disadvantages: the precise
definition of the area in which the microscopic model should be applied is
problem-dependent and needs to be specified. If the microscopic description is
applied (almost) everywhere, the computational cost may become even larger
than the cost required to numerically solve a micro- or mesoscale model.
Clearly, bridging the large range of dynamically coupled scales is a
fundamental challenge that was and still is a driving force in the development
of new mathematical algorithms.

The goal of the present paper is to derive and analyse a new hybrid continuum-
kinetic model for complex fluids. We start first by introducing a general
methodological approach, which allows us to derive a class of new multiscale
models. The basic idea is based on two simple ingredients. First, we assume
that the problems we are studying have fast and slow scale dynamics. In
addition, we assume a scale separation between the microscopic (fast) and
macroscopic (slow) dynamics and thus, the phenomena may be imagined to act at
different domain scales. Second, we assume that the effects of the fast scale
dynamics on microscopic scale can be captured at the macroscopic level, at
least approximately, by homogenization of the microscopic properties of the
fluid over a finite size domain. In this work, we consider a prototype
situation of isentropic flows governed by the conservation of mass and
momentum at the continuum level. In order to model rheology of a complex
fluid, the non-Newtonian stress tensor is obtained by an upscaling
homogenizing procedure using the kinetic relaxation type equations. We will
study properties of the derived hybrid multiscale model using a linear
stability analysis and successively conduct several numerical experiments to
illustrate the accuracy and efficiency of the proposed hybrid model.

The rest of the paper is organized in the following way. In Section 2, we
derive the hybrid continuum-kinetic model including the non-Newtonian stress
tensor. The latter is obtained by the Irving-Kirkwood formula which represents
upscale microscopic effects. In Section 3, we perform a linear stability
analysis which shows that for the regimes of interest, the low speed flows,
the model is linearly stable. In Section 4, several numerical examples are
conducted that demonstrate the validity of the proposed model for different
prototype situations arising in complex flows. Finally, Section 5 is dedicated
to the discussion of the obtained results and future developments.

##  2 The hybrid model approach

In this section, we first detail the general framework which will be used to
derive a class of hybrid multiscale models. We then focus on the derivation
and analysis of a prototype case, a hybrid continuum-kinetic model for complex
fluids. This is realized by homogenization over a fixed size cell of the
microscopic domains composed by the fluid molecules. The latter process leads
to the Irving-Kirkwood formula for complex macroscopic stress tensor.

We start by considering the following general setting. We assume that a
microscopic process used to describe time evolution of the state of the system
of interest is known. For instance, molecular dynamics or kinetic mesoscopic
equations are able to provide such information accurately enough. We also
assume that we have at our disposal a macroscopic model in which a missing
information will be provided by means of a microscopic model. These two
models, micro/meso and macro, can be related through a reconstruction
(upscaling) and a projection (downscaling) operator which permits to commute
from one system to the other. The upscaling operator averages the micro
properties of the fluid up to a coarser description, while the downscaling
operator uses the coarser description to obtain the unknown variables at the
microscopic level. A classical example of the above operators are
downscaling/upscaling operators between kinetic and macroscopic description.
In this case, the upscaling procedure is obtained by the integration of the
distribution function multiplied by the so-called collision invariants over
the velocity space while the downscaling operator is recovered from the
knowledge of the macroscopic variable defining a so-called equilibrium
distribution. In such a setting, our aim is to be able to give a description
of the state of the system by working on a given macroscopic grid, defined
a-priori, and through the use of a macroscopic model bringing some information
from the microscopic/mesoscopic dynamics. Schematically, the proposed method
works as depicted in Figure 2.1. We upscale the microsolver information, which
is determined by time evolution of the unknown $\bm{u}$, through a
homogenization on a box of fixed size $[-\alpha,\alpha]^{d}$ with $d$ the
spatial dimension and $\alpha$ the characteristic length of the microscopic
variation. The size of the box depends on the problem under consideration.
This micro information is successively used into the macroscopic solver, with
unknown $\bm{U}$, over the macroscopic grid to update the solution of the
complex flow.

In what follows, we describe the details of such procedure in the case, where
the macroscopic model is represented by the isentropic fluid equations for
complex fluids while the microscopic model is the BGK kinetic equation [3].
From the depicted scenario, it is possible to imagine alternative types of
combination of micro and macro dynamics which will be discussed in future
investigations.

Figure 2.1: Sketch of the hybrid multiscale approach.

###  2.1 A prototype case and derivation of the model

In this section, we consider a fluid which, by hypothesis, can be described at
the macroscopic level by the compressible isentropic fluid equations. In this
case, the system is governed by two equations describing the conservation of
mass and momentum. This system reads

$\begin{split}&\partial_{t}\rho+\nabla_{\bm{x}}\cdot(\rho\bm{u})=0,\\\
&\partial_{t}(\rho\bm{u})+\nabla_{\bm{x}}\cdot(\rho\bm{u}\otimes\bm{u})=\nabla_{\bm{x}}\cdot\mathbb{T},\end{split}$
(2.1)

where $\rho$ and $\bm{u}$ stand for the fluid density and velocity,
respectively. The so-called Cauchy stress tensor $\mathbb{T}$ describes
specific rheological fluid properties. All unknowns are functions of space
$\bm{x}\in\mathbb{R}^{d}$ and time $t>0$. For inviscid fluids,
$\mathbb{T}=-p\mathbb{I}$, where $\mathbb{I}$ is the identity tensor. This
expression leads to the Euler equations. Considering isentropic fluids,
pressure $p=p(\rho)$ is a given function of density $\rho.$

If viscous effects are taken into account, the Cauchy stress tensor reads

$\mathbb{T}=-p\mathbb{I}+\mathbb{S},$

where $\mathbb{S}$ stands for the viscous stress. For Newtonian fluids, the
latter is given by the Newton rheological rule (NRR)

$\mathbb{S}=\mu\left(\nabla_{\bm{x}}\bm{u}+\nabla_{\bm{x}}\bm{u}^{T}-\frac{2}{d}\nabla_{\bm{x}}\cdot\bm{u}\mathbb{I}\right)+\lambda\nabla_{\bm{x}}\cdot\bm{u}\mathbb{I},$

with $\mu>0$ and $\lambda\geq 0$ being the constant shear and bulk viscosity
coefficients, respectively. This relation leads to the compressible Navier-
Stokes equations. In complex fluids, however, the rheological relation for the
Cauchy stress tensor is more general and typically obtained by computational
or physical experiments. Here, we propose instead that the Cauchy stress
tensor takes into account microscopic effects in a homogenized way. This new
model is derived from considering the kinetic equations underpinning fluid
models of the type (2.1) as a microscopic model.

More specifically, we make the hypothesis that the so-called BGK (Bhatnagar-
Gross-Krook) equation is a suitable microscopic model which can be upscaled to
yield a missing information on a complex Cauchy stress tensor at the
macroscopic level. The considered BGK equation reads [3]

$\partial_{t}f+\bm{v}\cdot\nabla_{\bm{x}}f=Q(f),$ (2.2)

where $f$ is the probability density function of fluid molecules at position
$\bm{x}$ having velocity $\bm{v}\in\mathbb{R}^{d}$ at time $t$ and $Q(f)$ is a
collision operator (modelling the molecular interactions) to be described
later on. This hypothesis can be relaxed giving rise to similar models with,
however, different coefficients.

We assume now as depicted in Figure 2.1 that $f$ has slow variations at the
domain scale $D$ and fast variations at scale $\alpha\ll D$. Given that
hypothesis, we consequently set
$f=\tilde{f}(\bm{x},\frac{\bm{x}}{\alpha},\bm{v},t)$ where
$\tilde{f}(\bm{x},\bm{y},\bm{v},t)$ is 2-periodic with respect to $\bm{y}$,
with unit cell $[-1,1]^{d}$. We also suppose that the collision operator has
magnitude $1/\alpha$, i.e. that the microscopic spatial effect are balanced by
the collision dynamics among molecules. Substituting this representation into
(2.2), yields the following microscopic model

$\partial_{t}\tilde{f}+\bm{v}\cdot\nabla_{\bm{x}}\tilde{f}+\frac{1}{\alpha}\bm{v}\cdot\nabla_{\bm{y}}\tilde{f}=\frac{1}{\alpha}Q(\tilde{f}).$
(2.3)

Now, we introduce two new functions, one for the density
$\tilde{\rho}=\tilde{\rho}(\bm{x},\bm{y},t)$ and one for the vector velocity
$\tilde{\bm{u}}=\tilde{\bm{u}}(\bm{x},\bm{y},t)$ to be precisely defined later
on. $Q(\tilde{f})$ is the BGK operator, i.e. a relaxation operator towards a
given Maxwellian distribution $M_{(\tilde{\rho},\tilde{\bm{u}})}$ of parameter
$\tilde{\rho}$ and $\tilde{\bm{u}}$. We observe at this stage that the way in
which $\tilde{\rho}$ and $\tilde{\bm{u}}$ will be defined will permit to close
the resulting system of macroscopic equations. Denoting by $T$ the (constant)
temperature and $R$ the universal gas constant,
$M_{(\tilde{\rho},\tilde{\bm{u}})}$ is given by

$M_{(\tilde{\rho},\tilde{\bm{u}})}(\bm{v})=\frac{\tilde{\rho}}{(2\pi
RT)^{d/2}}\exp{\left(-\frac{|\bm{v}-\tilde{\bm{u}}|^{2}}{2RT}\right)},$ (2.4)

and the BGK operator is given by

$Q(\tilde{f})=-\frac{1}{\tilde{\tau}}(\tilde{f}-M_{(\tilde{\rho},\tilde{\bm{u}})}),$
(2.5)

where $\tilde{\tau}$ is a relaxation time. Finally, (2.3) takes the form

$\alpha\big{(}\partial_{t}\tilde{f}+\bm{v}\cdot\nabla_{\bm{x}}\tilde{f}\big{)}+\bm{v}\cdot\nabla_{\bm{y}}\tilde{f}=-\frac{1}{\tilde{\tau}}(\tilde{f}-M_{(\tilde{\rho},\tilde{\bm{u}})}).$
(2.6)

The microscopic model (2.6) can be used to provide a missing detailed
information for the macroscopic model (2.1). To this end, the Cauchy stress
tensor $\mathbb{T}$ can be reconstructed from $\tilde{f}$ by an upscaling
homogenization procedure through the Irving-Kirkwood (or Kramer) formula [23]:

$\mathbb{T}(\bm{x},t)=-\frac{1}{2^{d}}\int\displaylimits_{\left[-1,1\right]^{d}}\int_{\mathbb{R}^{d}}(\bm{v}-\tilde{\bm{u}}(\bm{x},\bm{y},t))\otimes(\bm{v}-\tilde{\bm{u}}(\bm{x},\bm{y},t))\tilde{f}(\bm{x},\bm{y},\bm{v},t)d\bm{v}d\bm{y}$
(2.7)

with $\tilde{f}$ being the solution of (2.6). In other words, we measure the
microscopic effects by homogenization of the kinetic model over the
microscopic box of size $[-\alpha,\alpha]^{d}$. Problem (2.6) can be viewed as
a cell problem for $\bm{y}$ in the unit cell $\left[-1,1\right]^{d}$ and
$\bm{v}\in{\mathbb{R}}^{d}$. The size of this box is left as a free parameter,
which may depend on the system under consideration. An interesting setting,
which will be numerically explored in Section 4, consists of considering the
size of the box into relation with the macroscopic description of the flow at
the numerical level. This can be done, for instance, by fixing a ratio between
the box where the microscopic effects are measured and the size of the mesh
employed at the macroscopic level.

In order to have a well posed problem, we finally supplement (2.6) with the
periodic boundary conditions for the probability density function:

$\tilde{f}(\bm{x},\bm{y}+\bm{n},\bm{v},t)=\tilde{f}(\bm{x},\bm{y},\bm{v},t),\quad\forall\bm{n}\in(2\mathbb{Z})^{d}.$

Our aim now is to find an approximate solution of the cell problem by means of
the Chapman-Enskog expansion [6] for small relaxation times $\tilde{\tau}$.
This will allow us to upscale the microscopic effects at the macroscopic level
without resorting to the numerical resolution of the microscopic model which
would lead to an expensive computation. This step is performed in the next
section.

###  2.2 Chapman-Enskog expansion for the prototype hybrid model

In this part, we present a perturbation analysis of the distribution function
$\tilde{f}$ over the box $[-1,1]^{d}.$ We assume in the rest of the paper that
$\tilde{\tau}=O(\varepsilon),$ where $\varepsilon\ll 1$ is the Knudsen number,
i.e. the ratio of the relaxation parameter $\tilde{\tau}$ to a characteristic
hydrodynamic temporal scale. The latter represents the macroscopic scale of
interest in our analysis. More precisely, we write
$\tilde{\tau}=\varepsilon\hat{\tau}$, where $\hat{\tau}=O(1)$ when
$\varepsilon\rightarrow 0$ in the following. We proceed by introducing the
Chapman-Enskog expansion [6] that will be truncated at the first-order terms

$\tilde{f}^{\varepsilon}=\tilde{f}^{(0)}+\varepsilon\tilde{f}^{(1)}+{\mathcal{O}}(\varepsilon^{2}),$
(2.8)

where $\tilde{f}^{\varepsilon}$ is the solution of (2.6) when $\tilde{\tau}$
is replaced by $\varepsilon\hat{\tau}$. The so-called cell problem (2.3) is
then rewritten, using the Einstein summation convention under the above
scaling as

$\Big{(}\alpha\big{(}\partial_{t}\tilde{f}^{\varepsilon}+v_{i}\partial_{x_{i}}\tilde{f}^{\varepsilon}\big{)}+v_{i}\partial_{y_{i}}\tilde{f}^{\varepsilon}\Big{)}(\bm{x},\bm{y},\bm{v},t)=-\frac{1}{\varepsilon\hat{\tau}}\Big{[}\tilde{f}^{\varepsilon}(\bm{x},\bm{y},\bm{v},t)-M_{(\tilde{\rho}^{\varepsilon}(\bm{x},\bm{y},t),\tilde{\bm{u}}^{\varepsilon}(\bm{x},\bm{y},t))}(\bm{v})\Big{]}.$
(2.9)

Substituting (2.8) into (2.9) and equating the same powers of $\varepsilon$
yields

$\tilde{f}^{(0)}=M_{(\tilde{\rho},\tilde{\bm{u}})},\quad\tilde{f}^{(1)}=-\hat{\tau}\Big{(}\alpha\big{(}\partial_{t}\tilde{f}^{(0)}+v_{i}\partial_{x_{i}}\tilde{f}^{(0)}\big{)}+v_{i}\partial_{y_{i}}\tilde{f}^{(0)}\Big{)}.$

We proceed now, using the first order expansion (2.9), to the computation of
the Cauchy stress tensor appearing in equation (2.1) in the above described
setting. To this end, we denote

$\displaystyle\mathcal{T}_{ij}$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{d}}(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})\tilde{f}^{(0)}d\bm{v}+\varepsilon\int_{\mathbb{R}^{d}}(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})\tilde{f}^{(1)}d\bm{v}+{\mathcal{O}}(\varepsilon^{2})$
$\displaystyle=:$
$\displaystyle\mathcal{T}_{ij}^{(0)}+\mathcal{T}_{ij}^{(1)}+{\mathcal{O}}(\varepsilon^{2}),$

the contribution to the stress tensor coming from the Chapmann-Enskog
expansion and we consider in the following computations each term separately.

First, using the equation of state for perfect gases
$\tilde{p}=RT\tilde{\rho}$, we obtain that the leading term in (2.2) can be
written as

$\displaystyle\mathcal{T}^{(0)}_{ij}$
$\displaystyle=\int_{\mathbb{R}^{d}}(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})\tilde{f}^{(0)}d\bm{v}$
(2.11a)
$\displaystyle=\int_{\mathbb{R}^{d}}(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})\frac{\tilde{\rho}}{(2\pi
RT)^{d/2}}\exp{\left(-\frac{|\bm{v}-\tilde{\bm{u}}|^{2}}{2RT}\right)}d\bm{v}=\tilde{\rho}RT\delta_{ij}=\tilde{p}\delta_{ij}$
or in the matrix form as $\mathcal{T}^{(0)}=\tilde{p}\mathbb{I}.$ (2.11b)

The second term in (2.2) reads:

$\displaystyle\mathcal{T}^{(1)}_{ij}$ $\displaystyle=$
$\displaystyle\varepsilon\int_{\mathbb{R}^{d}}(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})\tilde{f}^{(1)}d\bm{v}$
(2.12) $\displaystyle=$
$\displaystyle-\hat{\tau}\varepsilon\int_{\mathbb{R}^{d}}(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})\Big{(}\alpha\big{(}\partial_{t}\tilde{f}^{(0)}+v_{\ell}\partial_{x_{\ell}}\tilde{f}^{(0)}\big{)}+v_{\ell}\partial_{y_{\ell}}\tilde{f}^{(0)}\Big{)}d\bm{v}.$

We analyze now each term of equation (2.12) separately. We start by applying
the chain rule first to

$\partial_{{y_{\ell}}}\tilde{f}^{(0)}=\partial_{\tilde{\rho}}\tilde{f}^{(0)}\partial_{{y_{\ell}}}\tilde{\rho}+\partial_{\tilde{u}_{m}}\tilde{f}^{(0)}\partial_{y_{\ell}}\tilde{u}_{m}.$
(2.13)

We then observe that we can write the following relation

$\displaystyle v_{\ell}\frac{\partial\tilde{f}^{(0)}}{\partial y_{\ell}}$
$\displaystyle=v_{\ell}\frac{\partial\tilde{f}^{(0)}}{\partial\tilde{\rho}}\frac{\partial\tilde{\rho}}{\partial
y_{\ell}}+v_{\ell}\frac{\partial\tilde{f}^{(0)}}{\partial\tilde{u}_{m}}\frac{\partial\tilde{u}_{m}}{\partial
y_{\ell}}$ (2.14)
$\displaystyle=v_{\ell}\frac{\tilde{f}^{(0)}}{\tilde{\rho}}\frac{\partial\tilde{\rho}}{\partial
y_{\ell}}+\frac{v_{\ell}(v_{m}-\tilde{u}_{m})}{RT}\tilde{f}^{(0)}\frac{\partial\tilde{u}_{m}}{\partial
y_{\ell}}=:A_{1}+A_{2},$

to separate different contributions to the stress tensor, i.e. the one,
denoted by $A_{1}$, is due to the mass variation and the other one, denoted by
$A_{2}$, is due to the variation of the mean velocity in the cell
$[-1,1]^{d}$. Now, by introducing the notation
$\langle\cdot\rangle\equiv\int_{\mathbb{R}^{d}}\cdot d\bm{v}$ to indicate
integration over the velocity space, one can compute the moments for the terms
$A_{1}$ and $A_{2}$ defined in (2.14) separately. To that aim, with the
notation
$\tilde{A}_{k}=\langle(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})A_{k}\rangle$,
with $k=1,\,2$, for the second order moments of $A_{1}$ and $A_{2}$ we obtain:

$\displaystyle\tilde{A}_{1}$ $\displaystyle:=\
\frac{1}{\tilde{\rho}}\frac{\partial\tilde{\rho}}{\partial
y_{\ell}}\langle(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})v_{\ell}\tilde{f}^{(0)}\rangle$
$\displaystyle\qquad=\frac{1}{\tilde{\rho}}\frac{\partial\tilde{\rho}}{\partial
y_{\ell}}\langle(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})(v_{\ell}-\tilde{u}_{\ell})\tilde{f}^{(0)}\rangle+\frac{1}{\tilde{\rho}}\frac{\partial\tilde{\rho}}{\partial
y_{\ell}}\tilde{u}_{\ell}\langle(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})\tilde{f}^{(0)}\rangle$
$\displaystyle\qquad=\frac{\partial\tilde{\rho}}{\partial
y_{\ell}}\tilde{u}_{\ell}RT\delta_{ij},$ $\displaystyle\tilde{A}_{2}$
$\displaystyle:=\ \frac{1}{RT}\frac{\partial\tilde{u}_{m}}{\partial
y_{\ell}}\langle(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})(v_{m}-\tilde{u}_{m})v_{\ell}\tilde{f}^{(0)}\rangle$
$\displaystyle\qquad=\frac{1}{RT}\frac{\partial\tilde{u}_{m}}{\partial
y_{\ell}}\langle(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})(v_{m}-\tilde{u}_{m})(v_{\ell}-\tilde{u}_{\ell})\tilde{f}^{(0)}\rangle+\frac{\tilde{u}_{m}}{RT}\frac{\partial\tilde{u}_{m}}{\partial
y_{\ell}}\langle(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})(v_{m}-\tilde{u}_{m})\tilde{f}^{(0)}\rangle$
$\displaystyle\qquad=\frac{\partial\tilde{u}_{m}}{\partial
y_{\ell}}\tilde{\rho}RT\big{(}\delta_{ij}\delta_{m\ell}+\delta_{im}\delta_{j\ell}+\delta_{i\ell}\delta_{jm}\big{)}.$

We now analyze the term involving the time derivative
$\frac{\partial\tilde{f}^{(0)}}{\partial t}$ in (2.12). Similarly to (2.13) we
have

$\frac{\partial\tilde{f}^{(0)}}{\partial
t}=\frac{\tilde{f}^{(0)}}{\tilde{\rho}}\frac{\partial\tilde{\rho}}{\partial
t}+\frac{v_{m}-\tilde{u}_{m}}{RT}\tilde{f}^{(0)}\frac{\partial\tilde{u}_{m}}{\partial
t},$

and analogously to (2.14) with $\bm{y}$ replaced by $\bm{x}$ we have

$v_{\ell}\frac{\partial\tilde{f}^{(0)}}{\partial
x_{\ell}}=v_{\ell}\frac{\tilde{f}^{(0)}}{\tilde{\rho}}\frac{\partial\tilde{\rho}}{\partial
x_{\ell}}+\frac{v_{\ell}(v_{m}-\tilde{u}_{m})}{RT}\tilde{f}^{(0)}\frac{\partial\tilde{u}_{m}}{\partial
x_{\ell}}.$

This leads to the following contribution for what concerns the time derivative
of the distribution $\tilde{f}^{0}$

$\begin{split}&\tilde{A}_{3}:=\Big{\langle}(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})\frac{\partial\tilde{f}^{(0)}}{\partial
t}\Big{\rangle}=\Big{\langle}(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})\left(\frac{\tilde{f}^{(0)}}{\tilde{\rho}}\frac{\partial\tilde{\rho}}{\partial
t}+\frac{v_{m}-\tilde{u}_{m}}{RT}\tilde{f}^{(0)}\frac{\partial\tilde{u}_{m}}{\partial
t}\right)\Big{\rangle}=\\\
&=\Big{\langle}(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})\frac{\tilde{f}^{(0)}}{\tilde{\rho}}\frac{\partial\tilde{\rho}}{\partial
t}\Big{\rangle}+\Big{\langle}(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})\frac{v_{m}-\tilde{u}_{m}}{RT}\tilde{f}^{(0)}\frac{\partial\tilde{u}_{m}}{\partial
t}\Big{\rangle}=RT\frac{\partial\tilde{\rho}}{\partial
t}\,\delta_{ij}=\frac{\partial\tilde{p}}{\partial t}\,\delta_{ij},\end{split}$

while the contribution to the stress tensor coming from the space variation of
the quantities $\tilde{\rho}$, $\tilde{u}_{m}$ at the macroscopic scale
$\bm{x}$ is

$\tilde{A}_{4}:=\Big{\langle}(v_{j}-\tilde{u}_{j})(v_{i}-\tilde{u}_{i})v_{\ell}\frac{\partial\tilde{f}^{(0)}}{\partial
x_{\ell}}\Big{\rangle}=\frac{\partial\tilde{\rho}}{\partial
x_{\ell}}\tilde{u}_{\ell}RT\delta_{ij}+\frac{\partial\tilde{u}_{m}}{\partial
x_{\ell}}\tilde{\rho}RT\big{(}\delta_{ij}\delta_{m\ell}+\delta_{im}\delta_{j\ell}+\delta_{i\ell}\delta_{jm}\big{)}.$

Collecting all these terms, we obtain the following form of
$\mathcal{T}^{(1)}_{ij}=-\hat{\tau}\varepsilon\left(\tilde{A}_{1}+\tilde{A}_{2}+\alpha(\tilde{A}_{3}+\tilde{A}_{4})\right)$:

$\displaystyle\mathcal{T}^{(1)}_{ij}$ $\displaystyle=$
$\displaystyle-\varepsilon\hat{\tau}RT\tilde{\rho}\Big{[}\frac{\partial\tilde{u}_{m}}{\partial
y_{\ell}}\left(\delta_{ij}\delta_{m\ell}+\delta_{im}\delta_{j\ell}+\delta_{i\ell}\delta_{jm}\right)+\delta_{ij}\frac{\tilde{u}_{\ell}}{\tilde{\rho}}\frac{\partial\tilde{\rho}}{\partial
y_{\ell}}\Big{]}$ (2.15)
$\displaystyle-\varepsilon\alpha\hat{\tau}\Big{[}\frac{\partial\tilde{p}}{\partial
t}\delta_{ij}+RT\tilde{\rho}\Big{(}\frac{\partial\tilde{u}_{m}}{\partial
x_{\ell}}\left(\delta_{ij}\delta_{m\ell}+\delta_{im}\delta_{j\ell}+\delta_{i\ell}\delta_{jm}\right)+\delta_{ij}\frac{\tilde{u}_{\ell}}{\tilde{\rho}}\frac{\partial\tilde{\rho}}{\partial
x_{\ell}}\Big{)}\Big{]}.$

In the matrix form, the contribution to the stress tensor given by the
truncated Chapmann-Enskog expansion can be written as

$\displaystyle\mathcal{T}^{(1)}$ $\displaystyle=$
$\displaystyle-\varepsilon\hat{\tau}RT\tilde{\rho}\Big{[}\nabla_{\bm{y}}\tilde{\bm{u}}+\nabla_{\bm{y}}\tilde{\bm{u}}^{T}+\Big{(}\nabla_{\bm{y}}\cdot\tilde{\bm{u}}+\big{(}\frac{\tilde{\bm{u}}}{\tilde{\rho}}\cdot\nabla_{\bm{y}}\big{)}\tilde{\rho}\Big{)}\mathbb{I}\Big{]}$
(2.16)
$\displaystyle-\varepsilon\alpha\hat{\tau}\Big{[}\partial_{t}\tilde{p}\,\mathbb{I}+RT\tilde{\rho}\Big{[}\nabla_{\bm{x}}\tilde{\bm{u}}+\nabla_{\bm{x}}\tilde{\bm{u}}^{T}+\Big{(}\nabla_{\bm{x}}\cdot\tilde{\bm{u}}+\big{(}\frac{\tilde{\bm{u}}}{\tilde{\rho}}\cdot\nabla_{\bm{x}}\big{)}\tilde{\rho}\Big{)}\mathbb{I}\Big{]}=$
$\displaystyle=:$ $\displaystyle\mathcal{T}^{(1,1)}+\mathcal{T}^{(1,2)}.$

Combining now (2.7), (2.11b), (2.16), and discarding the terms
${\mathcal{O}}(\varepsilon^{2})$, yields

$\mathbb{T}(\bm{x},t)=-\frac{1}{2^{d}}\int_{[-1,1]^{d}}\left(\mathcal{T}^{(0)}+\mathcal{T}^{(1)}\right)d\bm{y}$

and after substituting plugging this into (2.1) we finally obtain

$\displaystyle\partial_{t}\rho+\nabla_{\bm{x}}\cdot(\rho\bm{u})=0,$ (2.17)
$\displaystyle\partial_{t}(\rho\bm{u})+\nabla_{\bm{x}}\cdot(\rho\bm{u}\otimes\bm{u})=$
$\displaystyle\hskip
28.45274pt\nabla_{\bm{x}}\cdot\Big{[}\frac{1}{2^{d}}\int\limits_{[-1,1]^{d}}\Big{\\{}-\tilde{p}\mathbb{I}+\varepsilon\hat{\tau}RT\tilde{\rho}\Big{(}\nabla_{\bm{y}}\tilde{\bm{u}}+\nabla_{\bm{y}}\tilde{\bm{u}}^{T}+\big{[}\nabla_{\bm{y}}\cdot\tilde{\bm{u}}+\big{(}\frac{\tilde{\bm{u}}}{\tilde{\rho}}\cdot\nabla_{\bm{y}}\big{)}\tilde{\rho}\big{]}\mathbb{I}\Big{)}$
$\displaystyle\hskip
56.9055pt\varepsilon\hat{\tau}\alpha\big{\\{}\partial_{t}\tilde{p}\,\mathbb{I}+RT\tilde{\rho}\Big{(}\nabla_{\bm{x}}\tilde{\bm{u}}+\nabla_{\bm{x}}\tilde{\bm{u}}^{T}+\big{[}\nabla_{\bm{x}}\cdot\tilde{\bm{u}}+\big{(}\frac{\tilde{\bm{u}}}{\tilde{\rho}}\cdot\nabla_{\bm{x}}\big{)}\tilde{\rho}\big{]}\mathbb{I}\Big{)}\big{\\}}\Big{\\}}d\bm{y}\Big{]}.$

We proceed in the next section by proposing a closing strategy for the system
(2.17). This will be established by relating the perturbations
$\tilde{\bm{u}},\tilde{\rho},\tilde{p}$ to the mean velocity $\bm{u}$ and
density $\rho$ of the fluid at the macroscopic scale. With this aim, we will
assume that $\tilde{\bm{u}}$ and $\tilde{\rho}$ have approximate polynomial
variations w.r.t. $\bm{y}$ in the unit cell. More precisely, we assume that
$\tilde{\rho}(\bm{x},\bm{y},t)\approx\rho(\bm{x}+\alpha\bm{y},t)$ and
similarly for $\tilde{\bm{u}}\approx\bm{u}(\bm{x}+\alpha\bm{y},t).$ We will
use the Taylor expansion and successive truncation with respect to the
parameter $\alpha$.

###  2.3 System closure

In order to get a closed system of equations we proceed with inserting our
assumption on a linear variation of the perturbed quantities $\tilde{\bm{u}}$
and $\tilde{\rho}$ into (2.17). We assume that both $\rho$ and $\bm{u}$ are
sufficiently regular to be Taylor-expanded from the center of the box. In our
expansion, we keep the terms up to order ${\mathcal{O}}(\alpha^{2})$. Thus,
$\rho(\bm{y},t),\bm{u}(\bm{y},t)$ can be written as:

$\displaystyle\tilde{\rho}(\bm{x},\bm{y},t)=\rho(\bm{x},t)+\alpha(\bm{y}\cdot\nabla_{\bm{x}})\rho(\bm{x},t)+\frac{\alpha^{2}}{2}(\bm{y}\cdot\nabla_{\bm{x}})^{2}\rho(\bm{x},t)+{\mathcal{O}}(\alpha^{3}),$
(2.18)
$\displaystyle\tilde{\bm{u}}(\bm{x},\bm{y},t)=\bm{u}(\bm{x},t)+\alpha(\bm{y}\cdot\nabla_{\bm{x}})\bm{u}(\bm{x},t)+\frac{\alpha^{2}}{2}(\bm{y}\cdot\nabla_{\bm{x}})^{2}\bm{u}(\bm{x},t)+{\mathcal{O}}(\alpha^{3}).$

Substituting (2.18) into (2.11a) and (2.15), integrating over the microscopic
box with respect to $\bm{y}$ and noting that odd terms with respect to
$\bm{y}$ cancel by antisymmetry, we obtain:

$\displaystyle-\frac{1}{2^{d}}\int\limits_{[-1,1]^{d}}\mathcal{T}^{(0)}\,d\bm{y}$
$\displaystyle=-\frac{1}{2^{d}}\int\limits_{[-1,1]^{d}}\tilde{\rho}(\bm{x},\bm{y},t)RT\mathbb{I}\,d\bm{y}$
(2.19)
$\displaystyle=-\frac{1}{2^{d}}\int\limits_{[-1,1]^{d}}\left[\rho+\alpha(\bm{y}\cdot\nabla_{\bm{x}})\rho+\frac{\alpha^{2}}{2}(\bm{y}\cdot\nabla_{\bm{x}})^{2}\rho\right]RT\mathbb{I}\,d\bm{y}+{\mathcal{O}}(\alpha^{3})$
$\displaystyle=-\left(\rho+\frac{\alpha^{2}}{6}\Delta_{\bm{x}}\rho\right)RT\mathbb{I}+{\mathcal{O}}(\alpha^{3}).$

Now taking the $\bm{y}$-derivative of the second equation in (2.18) relative
to the expansion of the velocity $\tilde{\bm{u}}$, we have

$\nabla_{\bm{y}}\tilde{\bm{u}}=\alpha\left(\nabla_{\bm{x}}\bm{u}+\alpha(\bm{y}\cdot\nabla_{\bm{x}})\nabla_{\bm{x}}\bm{u}\right)+{\mathcal{O}}(\alpha^{3}),$

with similar relations holding true for $(\nabla_{\bm{y}}\tilde{\bm{u}})^{T}$,
$(\nabla_{\bm{y}}\cdot\tilde{\bm{u}})$, and $(\nabla_{\bm{y}}\tilde{\rho})$.
Using again that odd terms in $\bm{y}$ cancel by antisymmetry, we get

$\begin{split}&-\frac{1}{2^{d}}\int\limits_{[-1,1]^{d}}\mathcal{T}^{(1,1)}\,d\bm{y}=\frac{\varepsilon\hat{\tau}RT}{2^{d}}\int\limits_{[-1,1]^{d}}\tilde{\rho}\big{(}\nabla_{\bm{y}}\tilde{\bm{u}}+(\nabla_{\bm{y}}\tilde{\bm{u}})^{T}+(\nabla_{\bm{y}}\cdot\tilde{\bm{u}})\,\mathbb{I}\big{)}\,d\bm{y}+\frac{\varepsilon\hat{\tau}RT}{2^{d}}\int_{[-1,1]^{d}}(\tilde{\bm{u}}\cdot\nabla_{\bm{y}})\tilde{\rho}\,\mathbb{I}\,d\bm{y}\\\
&=\frac{\varepsilon\hat{\tau}RT\alpha}{2^{d}}\int\limits_{[-1,1]^{d}}\Big{[}\rho+\alpha(\bm{y}\cdot\nabla_{\bm{x}})\rho\Big{]}\times\\\
&\phantom{mmmmmmmmmmmmm}\Big{[}\nabla_{\bm{x}}\bm{u}+(\nabla_{\bm{x}}\bm{u})^{T}+(\nabla_{\bm{x}}\cdot\bm{u})\,\mathbb{I}+\alpha(\bm{y}\cdot\nabla_{\bm{x}})\Big{(}\nabla_{\bm{x}}\bm{u}+(\nabla_{\bm{x}}\bm{u})^{T}+(\nabla_{\bm{x}}\cdot\bm{u})\mathbb{I}\Big{)}\Big{]}\,d\bm{y}\\\
&+\frac{\varepsilon\hat{\tau}RT\alpha}{2^{d}}\int\limits_{[-1,1]^{d}}\Big{[}\bm{u}+\alpha(\bm{y}\cdot\nabla_{\bm{x}})\bm{u}\Big{]}\cdot\Big{[}\nabla_{\bm{x}}\rho+\alpha(\bm{y}\cdot\nabla_{\bm{x}})(\nabla_{\bm{x}}\rho)\Big{]}\,\mathbb{I}d\bm{y}+{\mathcal{O}}(\alpha^{3})\\\
&=\varepsilon\hat{\tau}RT\alpha\Big{\\{}\rho\Big{[}\nabla_{\bm{x}}\bm{u}+(\nabla_{\bm{x}}\bm{u})^{T}+(\nabla_{\bm{x}}\cdot\bm{u})\,\mathbb{I}\Big{]}+\bm{u}\cdot\nabla_{\bm{x}}\rho\,\mathbb{I}\Big{\\}}+{\mathcal{O}}(\alpha^{3}).\end{split}$

Using similar computations for the term $\mathcal{T}^{(1,2)}$ in (2.16), we
have

$\displaystyle-\frac{1}{2^{d}}\int\limits_{[-1,1]^{d}}\mathcal{T}^{(1,2)}\,d\bm{y}$
$\displaystyle=\frac{1}{2^{d}}\varepsilon\hat{\tau}RT\alpha\int\limits_{[-1,1]^{d}}\Big{\\{}\partial_{t}\tilde{\rho}\,\mathbb{I}+\tilde{\rho}\big{(}\nabla_{\bm{x}}\tilde{\bm{u}}+(\nabla_{\bm{x}}\tilde{\bm{u}})^{T}+(\nabla_{\bm{x}}\cdot\tilde{\bm{u}})\,\mathbb{I}\big{)}+(\tilde{\bm{u}}\cdot\nabla_{\bm{x}})\tilde{\rho}\,\mathbb{I}\Big{\\}}\,d\bm{y}$
$\displaystyle=\frac{1}{2^{d}}\varepsilon\hat{\tau}RT\alpha\int\limits_{[-1,1]^{d}}\Big{\\{}\Big{(}\partial_{t}\rho+\alpha(\bm{y}\cdot\nabla_{\bm{x}})\frac{\partial\rho}{\partial
t}\Big{)}\,\mathbb{I}+\Big{[}\rho+\alpha(\bm{y}\cdot\nabla_{\bm{x}})\rho\Big{]}\,\Big{[}\nabla_{\bm{x}}\bm{u}+(\nabla_{\bm{x}}\bm{u})^{T}$
$\displaystyle\hskip
56.9055pt+(\nabla_{\bm{x}}\cdot\bm{u})\,\mathbb{I}+\alpha(\bm{y}\cdot\nabla_{\bm{x}})\Big{(}\nabla_{\bm{x}}\bm{u}+(\nabla_{\bm{x}}\bm{u})^{T}+(\nabla_{\bm{x}}\cdot\bm{u})\mathbb{I}\Big{)}\Big{]}$
$\displaystyle\hskip
56.9055pt+\Big{[}\bm{u}+\alpha(\bm{y}\cdot\nabla_{\bm{x}})\bm{u}\Big{]}\cdot\Big{[}\nabla_{\bm{x}}\rho+\alpha(\bm{y}\cdot\nabla_{\bm{x}})(\nabla_{\bm{x}}\rho)\Big{]}\,\mathbb{I}\Big{\\}}d\bm{y}+{\mathcal{O}}(\alpha^{3})$
$\displaystyle=\varepsilon\hat{\tau}RT\alpha\Big{\\{}\partial_{t}\rho\,\mathbb{I}+\rho\Big{[}\nabla_{\bm{x}}\bm{u}+(\nabla_{\bm{x}}\bm{u})^{T}+(\nabla_{\bm{x}}\cdot\bm{u})\,\mathbb{I}\Big{]}+\bm{u}\cdot\nabla_{\bm{x}}\rho\,\mathbb{I}\Big{\\}}+{\mathcal{O}}(\alpha^{3}).$

We finally assume that $\varepsilon=\sqrt{\alpha}$. This means that the
spatial inhomogeneity due to the microscopic interactions of order $\alpha$
are much smaller and specifically equal to the square of the relaxation
length, i.e. the mean free-path which is of order $\varepsilon$ as stated at
the beginning of Section 2.2. Thus, we drop all terms of order $\alpha^{2}$ in
front of terms of order $\varepsilon\alpha$. This is the case of the term
involving $\Delta_{x}\rho$ in (2.19). With this closure and using that
$\partial_{t}\rho=-\nabla_{\bm{x}}\cdot(\rho\bm{u})$ the system (2.1) is
written:

$\displaystyle\partial_{t}\rho+\nabla_{\bm{x}}\cdot(\rho\bm{u})=0,$ (2.20)
$\displaystyle\partial_{t}(\rho\bm{u})+\nabla_{\bm{x}}\cdot(\rho\bm{u}\otimes\bm{u})+\nabla_{\bm{x}}p$
$\displaystyle\hskip 62.59596pt=\alpha^{3/2}\hat{\tau}RT\
\nabla_{\bm{x}}\cdot\Big{\\{}2\rho\Big{(}\nabla_{\bm{x}}\bm{u}+(\nabla_{\bm{x}}\bm{u})^{T}\Big{)}+\nabla_{\bm{x}}\cdot(\rho\bm{u})\,\mathbb{I}\Big{\\}},$

with $p=\rho RT.$ We note that when terms of order $\alpha^{3/2}$ are
neglected, we recover the standard Euler equations.

######  Remark 1.

  * •

Using alternative closure relations with respect to the one defined in (2.18)
leads to alternative systems of balance laws which may be used to describe
different phenomena related to the microscopic dynamics.

  * •

Keeping additional terms in the Taylor expansion (2.18)/considering a
different scaling relation with respect to

$\varepsilon=\sqrt{\alpha}$

leads also to a different balance laws of dispersive type: third order terms
appears.

  * •

Different kinetic models can be applied instead of the BGK relaxation equation
(2.5). For example, molecular dynamics interacting through short range
potential can be used giving rise, for instance, to alternative Piola-
Kirchhoff stress tensors.

##  3 Linear stability analysis

In this section, we study the linear stability of system (2.20). Let us assume
$\rho=\rho_{0}+\beta\rho_{1}+\cdots$, and
$\bm{u}=\bm{u}_{0}+\beta\bm{u}_{1}+\cdots$ with $\beta\ll 1$, where
$(\rho_{0},\bm{u}_{0})$ represents a uniform steady state and
$(\rho-\rho_{0},\bm{u}-\bm{u}_{0})$ is a small perturbation from this
equilibrium state. We then substitute the above expressions into (2.20) and,
dropping terms of order $\beta^{2}$ or more, we obtain a system of linearized
equations for $\rho_{1}(\bm{x},t)$ and $\bm{u}_{1}(\bm{x},t)$. The scope of
this section is to investigate the stability properties of the resulting
linear system. The corresponding equations for $(\rho_{1},\bm{u}_{1})$ are

$\displaystyle\partial_{t}\rho_{1}+\rho_{0}\nabla_{\bm{x}}\cdot\bm{u}_{1}+\bm{u}_{0}\cdot\nabla_{\bm{x}}\rho_{1}=0,$
(3.1)
$\displaystyle\rho_{0}\left(\partial_{t}\bm{u}_{1}+(\bm{u}_{0}\cdot\nabla_{\bm{x}})\bm{u}_{1}\right)+RT\nabla_{\bm{x}}\rho_{1}$
$\displaystyle\hskip 62.59596pt=\alpha^{3/2}\hat{\tau}RT\
\nabla_{\bm{x}}\cdot\Big{\\{}2\rho_{0}\Big{(}\nabla_{\bm{x}}\bm{u}_{1}+(\nabla_{\bm{x}}\bm{u}_{1})^{T}\Big{)}+\Big{(}\rho_{0}\nabla_{\bm{x}}\cdot\bm{u}_{1}+\bm{u}_{0}\cdot\nabla_{\bm{x}}\rho_{1}\Big{)}\,\mathbb{I}\Big{\\}}.$

Using the Fourier transform in both space and time, we can write

$\begin{bmatrix}\rho_{1}\\\
\bm{u}_{1}\end{bmatrix}(\bm{x},t)=\begin{bmatrix}\bar{\rho}_{1}\\\
\bar{\bm{u}_{1}}\end{bmatrix}e^{i(\bm{k}\cdot\bm{x}-\omega t)}$ (3.2)

with $\bar{\rho}_{1}$ and $\bar{\bm{u}}_{1}$ being the Fourier coefficients
and $\bm{k}\in\mathbb{R}^{3}$, $\omega\in\mathbb{C}$. Using such
transformation, one can infer that the system is stable about the stationary
solution $(\rho_{0},\bm{u}_{0})$ if and only if the imaginary part of $\omega$
is nonpositive for all non-trivial solution (3.2). Moreover, one can state
that the system is stable if and only if it is stable for all
$(\rho_{0},\bm{u}_{0})\in[0,\infty)\times{\mathbb{R}}^{d}$.

For simplicity, before proceeding with the computations, we redefine the
following quantities $RT\equiv T$, $\alpha^{3/2}\hat{\tau}\equiv\tau$ and
$(\bar{\rho}_{1},\bar{\bm{u}}_{1})\equiv(\rho,\bm{u})$. Substituting now (3.2)
in (3.1) and using these redefined quantities, we get

$\displaystyle-i\omega\rho+\rho_{0}i(\bm{k}\cdot\bm{u})+i(\bm{k}\cdot\bm{u}_{0})\rho=0,$
$\displaystyle\rho_{0}\Big{(}-i\omega\bm{u}+i(\bm{k}\cdot\bm{u}_{0})\bm{u}\Big{)}+Ti\bm{k}\rho$
$\displaystyle\hskip 62.59596pt=\tau
Ti\bm{k}^{T}\Big{\\{}2\rho_{0}\Big{(}i\bm{k}\otimes\bm{u}+i\bm{u}\otimes\bm{k}\Big{)}+\Big{(}\rho_{0}i(\bm{k}\cdot\bm{u})+i(\bm{k}\cdot\bm{u}_{0})\rho\Big{)}\mathbb{I}\Big{\\}}.$

Simplifying the above expressions, we obtain

$\displaystyle\Big{(}-\omega+\bm{k}\cdot\bm{u}_{0}\Big{)}\rho+\rho_{0}\bm{k}\cdot\bm{u}=0,$
(3.3)
$\displaystyle\Big{(}-\omega+\bm{k}\cdot\bm{u}_{0}\Big{)}\bm{u}+\frac{T}{\rho_{0}}\Big{(}1-i\tau(\bm{k}\cdot\bm{u}_{0})\Big{)}\bm{k}\rho-2i\tau
T\Big{(}|\bm{k}|^{2}\bm{u}+(\bm{k}\cdot\bm{u})\bm{k}\Big{)}=0.$

Now, we want to analyze equation (3.3). To this aim, let us first assume that
$\bm{k}=0$. This choice directly implies $\omega\not=0$ otherwise, the
perturbed solution is just like the unperturbed equation: constant in time and
spatially uniform. From the above hypothesis, we immediately obtain
$-\omega\rho=0$ and $-\omega\bm{u}=0$. This means $\rho=0$ and $\bm{u}=0$ and
thus, the perturbed solution is again constant in time and space, i.e. there
does not exist any non-trivial solution to the system (3.3) in this setting.
Thus, in the sequel, we assume that $\bm{k}\not=0$.

Let thus, under the hypothesis $\bm{k}\not=0$, the orthogonal projection onto
$\\{\bm{k}\\}^{\perp}$ be denoted as $\mathtt{P}_{\bm{k}^{\perp}}$. Applying
$\mathtt{P}_{\bm{k}^{\perp}}$ to the second equation in (3.3), yields

$\Big{(}-\omega+\bm{k}\cdot\bm{u}_{0}-2i\tau
T|\bm{k}|^{2}\Big{)}\mathtt{P}_{\bm{k}^{\perp}}\bm{u}=0.$ (3.4)

It should be observed that the expression inside the parentheses in (3.4) can
not be zero, otherwise, by taking the imaginary part, one concludes that
$k=0$, which is a contradiction. Hence, the perturbed velocity is in the
direction of $\bm{k}$ and thus $\mathtt{P}_{\bm{k}^{\perp}}\bm{u}=0$. As a
consequence, we can write

$\bm{u}=v\frac{\bm{k}}{|\bm{k}|},$ (3.5)

where $v\in\mathbb{R}$. After substituting (3.5) into equation (3.3), the
system for $(\rho,v)$ reads

$\displaystyle\Big{(}-\omega+\bm{k}\cdot\bm{u}_{0}\Big{)}\rho+\rho_{0}|\bm{k}|v=0,$
$\displaystyle\Big{(}-\omega+\bm{k}\cdot\bm{u}_{0}\Big{)}v+\frac{T}{\rho_{0}}\Big{(}1-i\tau(\bm{k}\cdot\bm{u}_{0})\Big{)}|\bm{k}|\rho-4i\tau
T|\bm{k}|^{2}v=0.$

We introduce now the angle $\theta$ between $\bm{u}_{0}$ and $\bm{k}$, hence
$\bm{k}\cdot\bm{u}_{0}=ku_{0}\cos\theta$, where we denote $k=|\bm{k}|$ and
$u_{0}=|\bm{u}_{0}|$. This leads to

$\displaystyle\Big{(}-\omega+ku_{0}\cos\theta\Big{)}\rho+\rho_{0}kv=0,$
$\displaystyle\frac{T}{\rho_{0}}\Big{(}1-i\tau
ku_{0}\cos\theta\Big{)}k\rho+\Big{(}-\omega+ku_{0}\cos\theta-4i\tau
Tk^{2}\Big{)}v=0,$

or, in matrix form:

$\begin{bmatrix}-\omega+ku_{0}\cos\theta&\rho_{0}k\\\
\frac{T}{\rho_{0}}\Big{(}1-i\tau
ku_{0}\cos\theta\Big{)}k&-\omega+ku_{0}\cos\theta-4i\tau
Tk^{2}\end{bmatrix}\begin{bmatrix}\rho\\\ v\end{bmatrix}=0$

Thus, there exists a non trivial solution if and only if the determinant of
the above system is equal to zero. This gives

$\left|\begin{array}[]{llcc}-\omega+ku_{0}\mbox{cos}\theta&\rho_{0}k\\\
\frac{T}{\rho_{0}}\Big{(}1-i\tau
ku_{0}\mbox{cos}\theta\Big{)}k&-\omega+ku_{0}\mbox{cos}\theta-4i\tau
Tk^{2}\end{array}\right|=0.$

Letting $X=-\omega+ku_{0}\cos\theta$, the previous equation leads to the
following characteristic equation

$X^{2}-4i\tau Tk^{2}X-Tk^{2}(1-i\tau ku_{0}\cos\theta)=0.$ (3.6)

Let us consider the set of solutions of (3.6). We want to guarantee that the
imaginary part of $\omega$ remains non-positive such that the linear stability
holds true. To analyse (3.6), we refer to the Routh-Hurwitz criterion [24].
This criterion is summarized as follows: let a polynomial equation be of the
form

$(a_{0}+ib_{0})\mu^{n}+(a_{1}+ib_{1})\mu^{n-1}+\cdots+(a_{n}+ib_{n})=0.$ (3.7)

Then, all solutions of the above equation satisfy $\mbox{Re}(i\mu)<0$ if and
only if

$\displaystyle(-1)\Delta_{2}$ $\displaystyle=$
$\displaystyle-\left|\begin{array}[]{cc}a_{0}&a_{1}\\\
b_{0}&b_{1}\end{array}\right|>0,$ $\displaystyle(-1)^{2}\Delta_{4}$
$\displaystyle=$
$\displaystyle\left|\begin{array}[]{cccc}a_{0}&a_{1}&a_{2}&a_{3}\\\
b_{0}&b_{1}&b_{2}&b_{3}\\\ 0&a_{0}&a_{1}&a_{2}\\\
0&b_{0}&b_{1}&b_{2}\end{array}\right|>0,$ $\displaystyle\vdots$
$\displaystyle(-1)^{n}\Delta_{2n}$ $\displaystyle=$
$\displaystyle(-1)^{n}\left|\begin{array}[]{ccccccccc}a_{0}&a_{1}&\cdots&a_{n-1}&a_{n}&0&\cdots&\cdots&0\\\
b_{0}&b_{1}&\cdots&b_{n-1}&b_{n}&0&\cdots&\cdots&0\\\
0&a_{0}&\cdots&a_{n-2}&a_{n-1}&a_{n}&0&\cdots&0\\\
0&b_{0}&\cdots&b_{n-2}&b_{n-1}&b_{n}&0&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&a_{0}&a_{1}&a_{2}&a_{3}&\cdots&a_{n}\\\
0&0&\cdots&b_{0}&b_{1}&b_{2}&b_{2}&\cdots&b_{n}\end{array}\right|>0.$

Now, observing that $\textrm{Re}(iX)=\textrm{Im}(\omega)$, the stability
conditions follow immediately. We then apply the Routh-Hurwitz criterion
directly to (3.6). To this aim, we compare the coefficients of (3.6) with
those of (3.7). This gives

$\begin{array}[]{ll}a_{0}=1,&b_{0}=0\\\ a_{1}=0,&b_{1}=-4\tau Tk^{2}\\\
a_{2}=-Tk^{2},&b_{2}=\tau Tk^{3}u_{0}\cos\theta.\end{array}$ (3.11)

Using the above quantities in the Routh-Hurwitz criterion we get

$-\Delta_{2}=4\tau Tk^{2}>0$

and

$\Delta_{4}=\tau^{2}T^{2}k^{6}(16T-u_{0}^{2}\cos^{2}\theta).$

Thus, we can conclude that if $u_{0}^{2}\geq 16T$, then, $\Delta_{4}\geq 0$,
i.e. the model is stable for all values of $\theta$. On the other hand, if
$u_{0}^{2}>16T$, then, there is a threshold $\cos\theta_{c}$ defined by

$\cos^{2}\theta_{c}=\frac{16T}{u_{0}^{2}},$

such that the model is stable if $|\cos\theta|\leq\cos\theta_{c}$ and unstable
if $|\cos\theta|>\cos\theta_{c}$. In this case, the model is unstable for
waves propagating in directions close to the direction of the unperturbed
velocity.

To summarize the above stability analysis we have showed that when the flow is
subsonic the system is certainly linearly stable. Moreover, even for
supersonic flows the system remains stable, in fact for the model considered
the sound speed corresponds to $c=T$. Only for hypersonic flows (Mach number
$M>4$) we have the appearance of linearly unstable modes. However, for the
applications we have in mind, we always consider regimes where the fluid
velocity is small compared with the thermal velocity, and so, the stability
criterion is always verified.

##  4 Numerical methods and experiments

In this section, we study the behavior of the hybrid continuum-kinetic model
(2.20) derived in Section 2 by comparing it with the isentropic Euler
equations, the isentropic Navier-Stokes equations, and the BGK model. In what
follows, we shall refer to our hybrid multiscale model (2.20) as HMM. We
expect the HMM to improve the results of the standard macroscopic model and to
approach the behavior of the BGK equation, at least in specific regimes that
have been used in its derivation.

In the rest of this section we first describe the numerical methods used for
computations and then illustrate the performance of the derived hybrid
continuum-kinetic model in a number of experiments that clearly illustrate
that the HMM is capable to better describe the physics of complex fluids in
several regimes.

###  4.1 The BGK model and its numerical discretization

The HMM has been derived using the Chapmann-Enskog expansion starting from a
steady state BGK model, cf. Section 2.2. It is therefore natural to expect the
HMM to be close to the underlying kinetic equation in several regimes of the
relaxation parameter $\tau$. Consequently, we consider the following time
dependent BGK equation

$\partial_{t}f+\bm{v}\cdot\nabla_{\bm{x}}f=\frac{1}{\tau}(M_{(\rho,\bm{u})}-f),$
(4.1)

where as assumed in the derivation of the HMM, the Maxwellian distribution
(2.4) has constant temperature. In particular, we consider $T=1$ and also fix
the gas constant $R=1$. The shape of this equilibrium distribution is given in
equation (2.4). Now, the mean velocity $\bm{u}$ and the density $\rho$ are
those obtained from integration of the distribution $f$ in velocity space:

$\rho=\int_{\mathbb{R}^{d}}fd\bm{v},\qquad\bm{u}=\int_{\mathbb{R}^{d}}\bm{v}fd\bm{v}.$

For the sake of comparison with the other macroscopic models, we choose and
fix the dimension of the velocity space as $d=2$ both for the BGK as well as
for the HMM. Then, in order to numerically approximate (4.1) (see [13] for
details), we first replace the unbounded velocity space with a suitable
sufficiently large bounded set. This implies the truncation of the tails of
the distribution function, which normally lives in a non-compactly supported
set. Successively, we replace our continuous model by a so-called Discrete
Velocity Model (DVM) by discretizing this new bounded space by means of a
finite number of discrete points representing the discrete velocities that the
particles can assume. The result of this procedure is that the continuous BGK
model is replaced by $N$ linear transport equations coupled through a suitable
discretization of the relaxation operator $(M_{(\rho,\bm{u})}-f)$. We now
introduce the method and the notations, taking inspiration from [30]. We work
on a Cartesian grid $\mathcal{V}$ with

$\mathcal{V}=\left\\{\bm{v}_{\bm{k}}=\bm{k}\Delta\bm{v}+a,\
\bm{k}=(k^{(1)},k^{(2)}),\ \bm{a}=(a_{1},a_{2})\right\\},$

where $\bm{a}$ is an arbitrary vector, $\Delta\bm{v}$ is a constant mesh size
in velocity and where the components of the index $\bm{k}$ have some given
bounds $K^{(1)},K^{(2)}$. In this setting, the continuous distribution
function $f$ is replaced by the vector $f_{\mathcal{K}}(\bm{x},t)$ of size
$N$, where $N$ is chosen as a compromise between accuracy and computational
cost. Each component of this vector is assumed to be an approximation of the
distribution function $f$ at location $\bm{v}_{\bm{k}}$:

$f_{\mathcal{K}}(\bm{x},t)=(f_{\bm{k}}(\bm{x},t))_{\bm{x}},\qquad
f_{\bm{k}}(\bm{x},t)\approx f(\bm{x},\bm{v}_{\bm{k}},t).$

Thus, the discrete ordinate kinetic model consists of the following system of
ODEs to be solved

$\partial_{t}f_{\bm{k}}+\bm{v}_{\bm{k}}\cdot\nabla_{\bm{x}}f_{\bm{k}}=\frac{1}{\tau}(M_{\bm{k}}-f_{\bm{k}}),$
(4.2)

with $M_{\bm{k}}\approx M(\bm{x},\bm{v}_{\bm{k}},t)$ being a suitable
approximation of $M_{\rho,\bm{u}}(\bm{x},\bm{v}_{\bm{k}},t)$.

The system (4.2) is discretized in space using standard WENO approaches of
order three and we do not detail them here. With respect to the time
discretization, it should be observed, that the Maxwellian distribution in
(4.2) depends on the distribution function $f$ through its moments, cf. (2.4),
and hence the time integration of this ODE system is implemented using an
implicit-explicit (IMEX) Runge-Kutta method; see, e.g., [12] and references
therein. A general formulation of the IMEX Runge-Kutta method for (4.2) can be
written as

$\displaystyle F_{\bm{k}}^{(i)}$ $\displaystyle=f_{\bm{k}}^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\bm{v}_{\bm{k}}\cdot\nabla_{\bm{x}}F_{k}^{(j)}+\Delta
t\sum_{j=1}^{\nu}a_{ij}\frac{1}{\tau}\left(M^{(j)}_{\bm{k}}-F_{\bm{k}}^{(j)}\right),$
(4.3) $\displaystyle f_{\bm{k}}^{n+1}$ $\displaystyle=f_{\bm{k}}^{n}-\Delta
t\sum_{i=1}^{\nu}\tilde{w}_{i}\bm{v}_{\bm{k}}\cdot\nabla_{\bm{x}}F_{\bm{k}}^{(i)}+\Delta
t\sum_{i=1}^{\nu}w_{i}\frac{1}{\tau}\left(M_{\bm{k}}^{(i)}-F_{\bm{k}}^{(i)}\right),\
$

where the matrices $\tilde{A}=(\tilde{a}_{ij})$, $\tilde{a}_{ij}=0$ for $j\geq
i$ and $A=(a_{ij})$ are $\nu\times\nu$ matrices such that the resulting scheme
is explicit in $\bm{v}_{\bm{k}}\cdot\nabla_{\bm{x}}f$, and implicit in
$(M_{\rho,\bm{u}}-f)$. Here, we use the so-called second order in time
ARS(2,2,2) scheme [1], for which the coefficient vectors
$\tilde{w}=(\widetilde{w}_{1},..,\tilde{w}_{\nu})^{T}$,
$w=(w_{1},..,w_{\nu})^{T}$ are determined by the following double Butcher
tableau:

$\begin{array}[]{c|ccc}0&0&0&0\\\ \gamma_{1}&\gamma_{1}&0&0\\\
1&\gamma_{2}&1-\gamma_{2}&0\\\
\hline\cr&\gamma_{2}&1-\gamma_{2}&0\end{array}\qquad\begin{array}[]{c|ccc}0&0&0&0\\\
\gamma_{1}&0&\gamma_{1}&0\\\ 1&0&1-\gamma_{1}&\gamma_{1}\\\
\hline\cr&0&1-\gamma_{1}&\gamma_{1}\end{array}$

with $\gamma_{1}=1-1/\sqrt{2}$ and $\gamma_{2}=1-1/(2\gamma_{1})$. The above
scheme belongs to a particular class of IMEX methods for which the implicit
tableau is simply diagonally implicit, i.e. it is such that $a_{ij}=0$ if
$j>i$. This permits a direct evaluation of the implicit terms without
resorting to the inversion of non linear systems despite the nonlinearity of
the function which defines the equilibrium state $M_{\rho,\bm{u}}$. Indeed,
let us remark that the stage evaluation (4.3) can be rewritten as

$F_{\bm{k}}^{(i)}=\displaystyle f_{\bm{k}}^{n}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\bm{v}_{\bm{k}}\cdot\nabla_{\bm{x}}F_{\bm{k}}^{(j)}+\Delta
t\sum_{j=1}^{i-1}a_{ij}\frac{1}{\tau}\left(M^{(j)}_{\bm{k}}-F_{\bm{k}}^{(j)}\right)+\Delta
t\frac{a_{ii}}{\tau}\left(M^{(i)}_{\bm{k}}-F_{\bm{k}}^{(i)}\right),$ (4.4)

where the only implicit term is the diagonal factor
$\frac{a_{ii}}{\tau}\left(M^{(i)}_{\bm{k}}-F_{\bm{k}}^{(i)}\right)$, in which
$M_{\bm{k}}^{(i)}$ depends on the density and momentum of the distribution
function $(\rho,\bm{u})$. These macroscopic quantities can be obtained from
equation (4.4) by discrete integration in velocity space against
$\bm{\phi}_{\bm{k}}=(1,\bm{v}_{\bm{k}})$:

$\displaystyle(\rho^{(i)}(\bm{x}),\bm{u}^{(i)}(\bm{x}))$
$\displaystyle=\sum_{\bm{k}}(1,\bm{v}_{\bm{k}})F^{(i)}_{\bm{k}}(\bm{x})\,\Delta\bm{v}:=\langle\bm{\phi}_{\bm{k}}F_{\bm{k}}^{(i)}\rangle_{\mathcal{K}}$
(4.5)
$\displaystyle=\langle\bm{\phi}_{\bm{k}}f_{\bm{k}}^{n}\rangle_{\mathcal{K}}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\langle\bm{\phi}_{\bm{k}}\bm{v}_{\bm{k}}\cdot\nabla_{x}(F_{\bm{k}}^{(j)})\rangle_{\mathcal{K}}.$

As a consequence from the calculation performed in (4.5),
$(\rho^{(i)}(x),\bm{u}^{(i)}(x))$, and thus $M_{\bm{k}}^{(i)}$, can be
explicitly evaluated and then the scheme (4.3)-(4.4) is, in fact, explicitly
solvable.

The choice of this specific IMEX Runge-Kutta ODE solver reflects the facts
that we need to compare different model acting at different scales and the
scheme should allow to handle consistently the passage from the kinetic to the
fluid equations. In particular, when the scaling parameter $\tau\to 0$, the
scheme is formally equivalent to a discretization of the inviscid equation
(2.1). For small but non zero values of $\tau$, one can also expect the scheme
to be close to the new HMM. An analysis of the above scheme can be performed
showing that indeed it possesses the property of being consistent with the
limit macroscopic inviscid model (2.1) when $\tau\to 0$; see, e.g., [12].

###  4.2 A numerical method for the hybrid continuum-kinetic model (HMM)

We continue by introducing the numerical method for the HMM. We give the
details of the two-dimensional discretization since the following numerical
experiments will be restricted to two space dimensions. Let us fix, as for the
BGK case, the gas constant $R=1$ as well as the temperature $T=1$, set also
$\alpha^{3/2}\hat{\tau}\equiv\tau_{HMM}.$ HMM equations (2.20) rewritten in
the vector form read

${\bm{Q}}_{t}+{\bm{F}}_{x}+{\bm{G}}_{y}=\mathbf{S},$ (4.6)

where the vectors of conservative variables $\mathbf{Q}$, the fluxes
$({\bm{F}}(\mathbf{Q}),{\bm{G}}(\mathbf{Q})$ and the source terms
${\bm{S}}(\mathbf{Q})$ are

$\displaystyle{\bm{Q}}=\begin{bmatrix}\rho\\\ \rho u_{x}\\\ \rho
u_{y}\end{bmatrix},$ $\displaystyle{\bm{F}}(\mathbf{Q})=\begin{bmatrix}\rho
u_{x}\\\ \rho u_{x}^{2}+p\\\ \rho
u_{x}u_{y}\end{bmatrix},\quad{\bm{G}}(\mathbf{Q})=\begin{bmatrix}\rho u_{y}\\\
\rho u_{y}u_{x}\\\ \rho u_{y}^{2}+p\end{bmatrix},$ (4.7)
$\displaystyle{\bm{S}}(\mathbf{Q})=\tau_{HMM}\begin{bmatrix}0\\\
(\xi_{xx})_{x}+(\xi_{yx})_{y}\\\ (\xi_{xy})_{x}+(\xi_{yy})_{y}\\\
\end{bmatrix}.$

Here, $\bm{u}=(u_{x},u_{y})$ is the macroscopic velocity and

$\displaystyle\xi_{xx}=4\rho\frac{\partial u_{x}}{\partial
x}+\frac{\partial(\rho u_{x})}{\partial x}+\frac{\partial(\rho
u_{y})}{\partial y},\ \xi_{yy}=4\rho\frac{\partial u_{y}}{\partial
y}+\frac{\partial(\rho u_{y})}{\partial y}+\frac{\partial(\rho
u_{x})}{\partial x},$ (4.8)
$\displaystyle\xi_{xy}=\xi_{yx}=2\rho\left(\frac{\partial u_{x}}{\partial
y}+\frac{\partial u_{y}}{\partial x}\right),$

We also use the notation $\bm{x}=(x,y)\in\Omega$, where $\Omega$ is the
computational domain. We assume that
$\Omega=[x_{\min},x_{\max}]\times[y_{\min},y_{\max}]$ is paved with
$N_{x}\times N_{y}$ uniform cells of size $\Delta x\times\Delta y$. A cell is
labelled by two indices ${i,j}$, one for each direction, while when referring
to an interface between the cells we use respectively $({i+\frac{1}{2},j})$ or
$({i,j+\frac{1}{2}})$. Using this notation, the cell center is located at
point ${\bm{x}}_{i,j}=(x_{i},y_{j})$ while a face center lies at point

${\bm{x}}_{i+\frac{1}{2},j}=\left(\frac{1}{2}(x_{i}+x_{i+1}),y_{j}\right),\
{\bm{x}}_{i,j+\frac{1}{2}}=\left(x_{i},\frac{1}{2}(y_{j}+y_{j+1})\right).$

We denote any generic cell-centered quantity $m_{i,j}$, i.e. in the following
density, momentum, pressure and mean velocity. A conservative finite volume
scheme is adopted along with the first order explicit Euler scheme for time
integration, and we define the generic explicit operator $F[m_{i,j}^{n}]$
applied to $m_{i,j}^{n}$ as

$F[m_{i,j}^{n}]=m_{i,j}^{n}-\frac{\Delta t}{\Delta
x}\left(f_{i+\frac{1}{2},j}^{m}-f_{i-\frac{1}{2},j}^{m}\right)-\frac{\Delta
t}{\Delta y}\left(g_{i,j+\frac{1}{2}}^{m}-g_{i,j-\frac{1}{2}}^{m}\right),$

where the numerical fluxes $f_{i+\frac{1}{2},j}^{m}$ and
$g_{i,j+\frac{1}{2}}^{m}$ are of Rusanov type:

$\displaystyle f_{i+\frac{1}{2},j}^{m}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(f(m_{i+1,j}^{n})+f(m_{i,j}^{n})\right)-\frac{1}{2}|\lambda_{i+1/2,j}^{x,\text{max}}|\left(m_{i+1,j}^{n}-m_{i,j}^{n}\right),$
(4.9) $\displaystyle g_{i,j+\frac{1}{2}}^{m}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(g(m_{i,j+1}^{n})+g(m_{i,j}^{n})\right)-\frac{1}{2}|\lambda_{i,j+1/2}^{y,\text{max}}|\left(m_{i,j+1}^{n}-m_{i,j}^{n}\right)$
(4.10)

with the eigenvalues at the interfaces given by

$\displaystyle|\lambda_{i+\frac{1}{2},j}^{x,\text{max}}|=\max\left(|\lambda_{i+1,j}^{x,n}|,|\lambda_{i,j}^{x,n}|\right),\
|\lambda_{i,j+\frac{1}{2}}^{y,\text{max}}|=\max\left(|\lambda_{i,j+1}^{y,n}|,|\lambda_{i,j}^{y,n}|\right).$

In this setting, the density $\rho^{n}_{i,j}$ can be directly computed as

$\rho_{i,j}^{n+1}=F[\rho^{n}_{i,j}],$

and the components of the momentum equations satisfy

$\displaystyle(\rho{u})_{i,j}^{n+1}=F[(\rho{u})_{i,j}^{n}]+\frac{\Delta
t}{\Delta
x}\left((\xi_{xx})_{i+\frac{1}{2},j}^{n}-(\xi_{xx})_{i-\frac{1}{2},j}^{n}\right)+\frac{\Delta
t}{\Delta
y}\left((\xi_{yx})_{i,j+\frac{1}{2}}^{n}-(\xi_{yx})_{i,j-\frac{1}{2}}^{n}\right),$
$\displaystyle(\rho{v)}_{i,j}^{n+1}=F[(\rho{v})_{i,j}^{n}]+\frac{\Delta
t}{\Delta
x}\left((\xi_{xy})_{i+\frac{1}{2},j}^{n}-(\xi_{xy})_{i-\frac{1}{2},j}^{n}\right)+\frac{\Delta
t}{\Delta
y}\left((\xi_{yy})_{i,j+\frac{1}{2}}^{n}-(\xi_{yy})_{i,j-\frac{1}{2}}^{n}\right).$

According to the definitions (4.8) the computation of the stress components
$(\xi_{xx},\xi_{xy}$, $\xi_{yx},\xi_{yy})$ requires the knowledge of the
discrete velocity gradients. Here, they are computed on each boundary face of
the control volume using a trapezoidal quadrature rule, that is,

$\nabla_{\bm{x}}\bm{u}_{i+\frac{1}{2},j}=\frac{1}{4}\left(\nabla_{\bm{x}}\bm{u}_{i+\frac{1}{2},j-\frac{1}{2}}+\nabla_{\bm{x}}\bm{u}_{i+\frac{1}{2},j+\frac{1}{2}}+\nabla_{\bm{x}}\bm{u}_{i+\frac{1}{2},j-\frac{1}{2}}+\nabla_{\bm{x}}\bm{u}_{i+\frac{1}{2},j+\frac{1}{2}}\right),$
(4.11)

where the two-dimensional gradient of the velocity field $\bm{u}(t,\bm{x})$ is
obtained at each corner defined by subscript index
$({i+\frac{1}{2},j+\frac{1}{2}})$ as follows:

$\nabla_{\bm{x}}\bm{u}_{i+\frac{1}{2},j+\frac{1}{2}}=\frac{1}{2}\left[\begin{array}[]{l}\frac{\bm{u}_{i+1,j}-\bm{u}_{i,j}}{{\Delta
x}}+\frac{\bm{u}_{i+1,j+1}-\bm{u}_{i,j+1}}{{\Delta x}}\\\\[10.0pt]
\frac{\bm{u}_{i,j+1}-\bm{u}_{i,j}}{{\Delta
y}}+\frac{\bm{u}_{i+1,j+1}-\bm{u}_{i+1,j}}{{\Delta y}}\end{array}\right].$

Finally, to increase the accuracy in space, the cell-centred quantities
$m_{i,j}$ used in the definition of the numerical fluxes (4.9)-(4.10) are
replaced by a high order polynomial interpolation through the so-called WENO
reconstruction of order three.

####  4.2.1 A numerical method for the Navier-Stokes model

We briefly recall here the Navier-Stokes (NS) model which can obtained from
the standard Chapmann-Enskog expansion [6] in the case of isentropic flows and
is used for comparison purposes in the rest of the section. The NS model reads

$\displaystyle\partial_{t}\rho+\nabla_{\bm{x}}\cdot(\rho\bm{u})=0,$ (4.12)
$\displaystyle\partial_{t}(\rho\bm{u})+\nabla_{\bm{x}}\cdot(\rho\bm{u}\otimes\bm{u})+\nabla_{\bm{x}}p=\varepsilon
RT\
\nabla_{\bm{x}}\cdot\left(\rho\Big{(}\nabla_{\bm{x}}\bm{u}+(\nabla_{\bm{x}}\bm{u})^{T}\Big{)}\right).$

The numerical method follows the path of the discretization of the HMM model.
Thus, rewriting the system (4.12) as (4.6), where ${\bm{Q}},{\bm{F}},{\bm{G}}$
have the same meaning of (4.7) while the source term becomes

${\bm{S}}(\mathbf{Q})=\varepsilon RT\begin{bmatrix}0\\\
(\xi^{NS}_{xx})_{x}+(\xi^{NS}_{yx})_{y}\\\
(\xi^{NS}_{xy})_{x}+(\xi^{NS}_{yy})_{y}\\\ \end{bmatrix}.$

with

$\displaystyle\xi^{NS}_{xx}=2\rho\frac{\partial u_{x}}{\partial x},\
\xi^{NS}_{yy}=2\rho\frac{\partial u_{y}}{\partial y},\
\xi^{NS}_{xy}=\xi^{NS}_{yx}=\rho\left(\frac{\partial u_{x}}{\partial
y}+\frac{\partial u_{y}}{\partial x}\right).$

We then proceed, as done previously, using a first order explicit Euler scheme
in time and the same Rusanov flux with WENO reconstruction for the hyperbolic
fluxes and the same second order discrete velocity gradients (4.11) are used
for the viscous terms. This ends the description of the model and of the
numerical scheme.

###  4.3 Highly oscillating fluid

We start by considering a highly oscillating initial condition for a flow in
the computational domain $\Omega=[0,1]\times[0,2]$ and assume the periodic
boundary conditions to be imposed in both directions. The temperature $T$ as
well as the constant $R$ are set to one. The domain is paved with $N_{x}\times
N_{y}=64\times 128$ cells in the physical space while for what concerns the
BGK model we discretize the velocity space with $20\times 20$ cells with
$v_{x,max}=v_{y,max}=5$ and $v_{x,min}=v_{y,min}=-5$, where
$\bm{v}=(v_{x},v_{y})$. The initial distribution function is assumed to be at
local equilibrium, i.e. $f(\bm{x},\bm{v},0)=M_{\rho,\bm{u}}(\bm{x},\bm{v},0)$.
The initial density and velocities are

$\rho(x,y,0)=1+0.2\cos(10\pi x)\sin(12\pi
y),\quad\bm{u}(x,y,0)=(u_{x},u_{y})=(1,0).$

In Figure 4.1, we plot the highly oscillating density and the velocity
profiles after $100$ time iterations when the isentropic Euler model is used.
In Figure 4.2, we present the same result at a fixed value of $x=0.5$. In
particular, we compare the results obtained for the isentropic Euler, the BGK
and the HMM models using different values of the scaling parameter $\tau$. In
the BGK model we choose $\tau=0.001$, $\tau=0.005$ and $\tau=0.01$, while in
the HMM we set $\tau_{HMM}=\tau/3$. In Figure 4.2 we can clearly see that
oscillations are damped as fast as the scaling parameter becomes larger. In
these three tested situations, the HMM seems to be able to describe a kinetic
regime both for the case of the density as well as for the case of the mean
velocity $u_{x}(x,y,t)$. For sake of comparison, in Figure 4.3, we plot the
same density and velocity profiles at fixed $x=0.5$ where we added the results
of the Navier-Stokes model for $\varepsilon=\tau$. As one can see, for small
values of the damping parameters both NS and HMM give results very close to
those of the BGK model, while for $\tau=0.01$ the NS model tends to
overestimate the damping for the density and fails in describing the velocity
profile. Instead, the HMM model provides a very good description of the
density and while overdamping the velocity still is able to follow the
behavior of the BGK model.

Figure 4.1: Highly oscillating fluid. Profiles of density (left) and first
component of the velocity $u_{x}(x,y,t)$ (right) at final time for the
isentropic Euler equations.

Figure 4.2: Highly oscillating fluid. Comparison of the density (left) and
velocity $u_{x}(x,y,t)$ (right) profiles for $x=0.5$. The results for the
Euler, the BGK and the HMM equations with respectively $\tau=0.001$,
$\tau=0.005$ and $\tau=0.01$ are shown.

Figure 4.3: Highly oscillating fluid. Comparison of the density (left) and
velocity $u_{x}(x,y,t)$ (right) profiles for $x=0.5$. The results for the the
BGK, the Navier-Stokes and the HMM equations with respectively $\tau=0.001$,
$\tau=0.005$ and $\tau=0.01$ are shown.

###  4.4 Perturbed Couette flow

In this section, we consider a two-dimensional Couette-type flow in the domain
$\Omega=[0,1]\times[0,2]$ perturbed by high frequency waves. The setup is the
following. The boundary conditions in the $x$-direction are periodic while in
the $y$-direction we impose the so-called no-slip boundary conditions for
which the normal to the walls velocity is fixed to zero while the tangential
velocity is set equal to the velocity of the walls. The spatial domain is
discretized as before with $N_{x}\times N_{y}=64\times 128$ cells, while the
velocity space is approximated by using $20\times 20$ cells in the domain
$[-5,5]^{2}$. Local equilibrium for the distribution function is supposed at
time $t=0$. The initial density and velocities are

$\rho(x,y,0)=1+0.2\cos(10\pi x)\sin(12\pi y),\
\bm{u}(x,y,0)=(u_{x},u_{y})=(y(1+0.125\sin(8\pi x)),0).$

In other words, we simulate a steady state Couette flow with an added high
frequency perturbations both in the density and the velocity. In particular,
to the top wall and in the $y$-direction it is imposed an oscillatory velocity
in the $x$-direction which reads

$u_{y}(x,2,t)=1+0.25\sin(8\pi x).$

The temperature and the universal gas constant are set $T=1$ and $R=1$ as
before. In Figure 4.4, we show the density and the velocities profiles after
$100$ time iteration with a time step of $\Delta t=6.25\ 10^{-4}$. At this
time, the perturbed waves have not been completely damped out: it is still
possible to observe how the oscillations in velocity and density modify the
solution as the intensity of the relaxation parameter grows. In particular,
Figure 4.4 shows the solutions obtained with $\tau=0$, $\tau=0.005$ and
$\tau=0.01$ for the BGK model, while for the HMM we take $\tau_{HMM}=\tau/3$
as for the first test. In Figure 4.5, we present the same results at $x=0.5$
comparing the isentropic Euler, the BGK and the HMM solutions. In this
picture, we can clearly see that oscillations are damped when moving far from
the regime of validity of the isentropic Euler equations and that the HMM
captures the BGK solution very well for all regimes, considered at least for
the density and the first component of the velocity field. Finally, in Figure
4.6, we show a comparison between the HMM and the NS model using again the BGK
model as reference, all presented at $x=0.5$. For the specific case of the NS
model, we take $\varepsilon=\tau$ while all the others numerical parameters
are equal to the ones of the HMM model. From this figure we see that for small
values of $\tau$ both models capture very well the reference solution given by
the BGK equation, while for larger $\tau$, NS tends to overestimate the
damping for what concerns the density profile. Both models is capable of
capturing the first component of the velocity field $u_{x}$, while both fail
in the description of the second component $u_{y}$ when $\tau=0.01$. Let
observe anyway that $u_{y}$ is affected by waves, which are at least one order
of magnitude smaller than the perturbation waves acting on the other
macroscopic quantities and thus the error introduced in the solution remains
very small.

Figure 4.4: Perturbed Couette flow. Profiles of density and velocity for the
Euler and the BGK equations. From top to bottom: Density, x-velocity and
y-velocity profiles. Left $\tau=0$, middle $\tau=0.005$, right $\tau=0.01$.

Figure 4.5: Perturbed Couette flow. Comparison of the density (left),
$x$-velocity (middle) and $y$-velocity (right) profiles for $x=0.5$. The
results for the Euler, the BGK and the HMM equations with respectively
$\tau=0.001$, $\tau=0.005$ and $\tau=0.01$ are shown.

Figure 4.6: Perturbed Couette flow. Comparison of the density (left),
$x$-velocity (middle) and $y$-velocity (right) profiles for $x=0.5$. The
results for the the BGK, the Navier-Stokes and the HMM equations with
respectively $\tau=0.001$, $\tau=0.005$ and $\tau=0.01$ are shown.

###  4.5 A vortex test problem

In this last section, we consider a Taylor-Green-type problem. The initial
data are given by

$\rho(x,y,0)=1+0.1\cos(8\pi x)\sin(8\pi y),\
\bm{u}(x,y,0)=(u_{x},u_{y})=(\cos(x)\sin(y),-\cos(y)\sin(x)).$

The temperature is $T=1$, the universal gas constant $R=1$, and the
computational domain is $\Omega=[0,2\pi]^{2}$ discretized with
$N_{x}=N_{y}=128$ points in each direction. The velocity space is discretized
using the same parameters as in the previous tests. We take a time step
$\Delta t=2\ 10^{-4}$ and set a final time of $T_{f}=0.6$. The initial data
are represented in Figure 4.7 on the top left, while the final solution for
the isentropic Euler is reported in the same figure on the top right and
illustrates the vortices formation. The bottom image reports a three-
dimensional view of the same solution for the density. In Figure 4.8, we
present comparisons between the isentropic Euler, the BGK and our HMM at the
final time for the density and the velocities and for three different values
of the relaxation parameter $\tau=0.005$, $\tau=0.01$ and $\tau=0.025$ with
$\tau_{HMM}=\tau/3$. Once again, the obtained results indicate that the HMM is
able to capture microscopic structures with a sufficiently high accuracy for
the regimes under considaration.

Figure 4.7: Vortex flow. Profiles of density and velocity for the Euler and
the BGK equations. From top to bottom: Density initial data, density Euler
solution, density Euler solution three dimensional view and density BGK
solution with $\tau=0.01$. Red arrows represent the velocity vector field.

Figure 4.8: Vortex flow. Comparison of the density (left), $x$-velocity
(middle) and $y$-velocity (right) profiles for $x=0.5$. The results for the
Euler, the BGK and the HMM equations with respectively $\tau=0.005$,
$\tau=0.01$ and $\tau=0.025$ are shown.

##  5 Conclusions

In this work, we have introduced a new hybrid multiscale model coupling the
continuum and kinetic descriptions. We have applied the Chapman-Enskog
expansion on a suitably scaled stationary BGK equation and then upscaled the
kinetic contribution over a microscopic box. These microscopic regions are
located at each points where the fluid equation are solved. The probability
density function used to study the kinetic evolution is related to the unknown
density and velocity of the system through truncated Taylor expansions. Based
on a scaling assumption, related to the problem under consideration, a new
multiscale continuum-kinetic model has been derived. A linear stability
analysis shows that new hybrid multiscale model is conditionally linearly
stable in the strongly supersonic regime and unconditionally stable in the
subsonic and mildly supersonic one. Several numerical examples demonstrate
that the hybrid multiscale model is more accurate than standard fluid models
and represents complex flows at different regimes more precisely.

In future work we would like to derive different hybrid multiscale models
following the same general strategy outlined in this work: couple a
microscopic description with a macroscopic one through an upscaling of the
microscopic information obtained by homogenization of the micro quantities
over micro sized boxes.

Acknowledgment:

The work of A. Chertock was supported in part by NSF grants DMS-1818684. P.
Degond holds a visiting professor association with the Department of
Mathematics, Imperial College London, UK. The work of G. Dimarco was supported
by the Italian Ministry of Instruction, University and Research (MIUR) under
the PRIN Project 2017 (No. 2017KKJP4X). M. Lukáčová-Medvid’ová and A. Ruhi
were supported by the German Science Foundation (DFG) under the Collaborative
Research Center TRR 146 Multiscale Simulation Methods for Soft Matter Systems
(Project C5). M. Lukáčová-Medvid’ová is grateful to the Gutenberg Research
College and Mainz Institute of Multiscale Modelling for supporting her
research.

## References

  * [1] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri, Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations, Appl. Numer. Math., 25 (1997), pp. 151–167.  Special issue on time integration (Amsterdam, 1996).
  * [2] C. Bardos, F. Golse, and C. D. Levermore, Fluid dynamic limits of kinetic equations. II. Convergence proofs for the Boltzmann equation, Comm. Pure Appl. Math., 46 (1993), pp. 667–753.
  * [3] P. L. Bhatnagar, E. P. Gross, and M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev., 94 (1954), pp. 511–525.
  * [4] M. K. Borg, D. A. Lockerby, and J. M. Reese, Fluid simulations with atomistic resolution: a hybrid multiscale method with field-wise coupling, J. Comput. Phys., 255 (2013), pp. 149–165.
  * [5] , A hybrid molecular–continuum method for unsteady compressible multiscale flows, J. Fluid Mech., 768 (2015), p. 388–414.
  * [6] C. Cercignani, Mathematical methods in kinetic theory, Plenum Press, New York, second ed., 1990.
  * [7] A. Crestetto, N. Crouseilles, G. Dimarco, and M. Lemou, Asymptotically complexity diminishing schemes (ACDS) for kinetic equations in the diffusive scaling, J. Comput. Phys., 394 (2019), pp. 243–262.
  * [8] R. Datta, L. Yelash, F. Schmid, F. Kummer, M. Oberlack, M. Lukáčová-Medviďová, and P. Virnau, Shear-thinning in oligomer melts-molecular origin and applications, Polymers, 13 (2021), p. 2806.
  * [9] P. Degond, G. Dimarco, and L. Pareschi, The moment-guided Monte Carlo method, Internat. J. Numer. Methods Fluids, 67 (2011), pp. 189–213.
  * [10] P. J. Dellar, Two routes from the Boltzmann equation to compressible flow of polyatomic gases, Prog. Comput. Fluid Dyn., 8 (2008), pp. 84–96.
  * [11] G. Dimarco and L. Pareschi, Fluid solver independent hybrid methods for multiscale kinetic equations, SIAM J. Sci. Comput., 32 (2010), pp. 603–634.
  * [12] G. Dimarco and L. Pareschi, Asymptotic preserving implicit-explicit Runge-Kutta methods for nonlinear kinetic equations, SIAM J. Numer. Anal., 51 (2013), pp. 1064–1087.
  * [13] , Numerical methods for kinetic equations, Acta Numer., 23 (2014), pp. 369–520.
  * [14] W. E, Principles of multiscale modeling, Cambridge University Press, Cambridge, 2011.
  * [15] W. E and B. Engquist, The heterognous multiscale methods, Commun. Math. Sci., 1 (2003), pp. 87–132.
  * [16] W. E, B. Engquist, X. Li, W. Ren, and E. Vanden-Eijnden, Heterogeneous multiscale methods: A review, Commun. Comput. Phys., 2 (2007), pp. 367–450.
  * [17] W. E and X. Li, Analysis of the heterogeneous multiscale method for gas dynamics, Methods Appl. Anal., 11 (2004), pp. 557–572.
  * [18] W. E and J. Lu, Seamless multiscale modeling via dynamics on fiber bundles, Commun. Math. Sci., 5 (2007), pp. 649–663.
  * [19] W. E, P. Ming, and P. Zhang, Analysis of the heterogeneous multiscale method for elliptic homogenization problems, J. Amer. Math. Soc., 18 (2005), pp. 121–156.
  * [20] D. A. Fedosov and G. E. Karniadakis, Triple-decker: interfacing atomistic-mesoscopic-continuum flow regimes, J. Comput. Phys., 228 (2009), pp. 1157–1171.
  * [21] F. Golse and L. Saint-Raymond, The incompressible Navier-Stokes limit of the Boltzmann equation for hard cutoff potentials, J. Math. Pures Appl. (9), 91 (2009), pp. 508–552.
  * [22] H. Grad, Asymptotic equivalence of the Navier-Stokes and nonlinear Boltzmann equations, in Proc. Sympos. Appl. Math., Vol. XVII, Amer. Math. Soc., Providence, R.I., 1965, pp. 154–183.
  * [23] J. H. Irving and J. G. Kirkwood, he statistical mechanical theory of transport processes. IV. The equations of hydrodynamics, J. Chem. Phys., 18 (1950), pp. 508–552.
  * [24] Y. Ishida and T. Yamamoto, Linear And Nonlinear Rotordynamics: a modern treatment with applications, John Wiley & Sons, 2013.
  * [25] I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevrekidis, O. Runborg, and C. Theodoropoulos, Equation-free, coarse-grained multiscale computation: enabling microscopic simulators to perform system-level analysis, Commun. Math. Sci., 1 (2003), pp. 715–762.
  * [26] I. G. Kevrekidis and G. Samaey, Equation-free multiscale computation: algorithms and applications, Annu. Rev. Phy.s Chem., 60 (2009), pp. 321–344.
  * [27] P. Koumoutsakos, Multiscale flow simulations using particles, vol. 37 of Annu. Rev. Fluid Mech., Palo Alto, 2005, pp. 457–487.
  * [28] C. D. Levermore and N. Masmoudi, From the Boltzmann equation to an incompressible Navier-Stokes-Fourier system, Arch. Ration. Mech. Anal., 196 (2010), pp. 753–809.
  * [29] F.-H. Lin, C. Liu, and P. Zhang, On hydrodynamics of viscoelastic fluids, Comm. Pure Appl. Math., 58 (2005), pp. 1437–1471.
  * [30] L. Mieussens, Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamics, Mathematical Models and Methods in Applied Sciences, 10 (2000), pp. 1121–1149.
  * [31] J. Nečas and I. Hlaváček, Mathematical theory of elastic and elasto-plastic bodies: an introduction, vol. 3 of Studies in Applied Mechanics, Elsevier, 1980.
  * [32] G. A. Radtke, J.-P. M. Péraud, and N. G. Hadjiconstantinou, On efficient simulations of multiscale kinetic transport, Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 371 (2013), p. 20120182.
  * [33] K. K. Rao and P. R. Nott, An introduction to granular flow, Cambridge Series in Chemical Engineering, Cambridge University Press, 2008.
  * [34] W. Ren, Analytical and numerical study of coupled atomistic-continuum methods for fluids, J. Comput. Phys., 227 (2007), pp. 1353–1371.
  * [35] , Seamless multiscale modeling of complex fluids using fiber bundle dynamics, Commun. Math. Sci., 5 (2007), pp. 1027–1037.
  * [36] W. Ren and W. E, Heterogeneous multiscale method for the modeling of complex fluids and micro-fluidics, J. Comput. Phys., 204 (2005), pp. 1–26.
  * [37] S. E. Spagnolie, ed., Complex fluids in biological systems, Biological and Medical Physics, Biomedical Engineering, Springer, New York, 2015\.  Experiment, theory, and computation.
  * [38] S. Stalter, L. Yelash, N. Emamy, A. Statt, M. Hanke, M. Lukáčová-Medviďová, and P. Virnau, Molecular dynamics simulations in hybrid particle-continuum schemes: Pitfalls and caveats, Comput. Phys. Commun., 224 (2018), pp. 198–208.
  * [39] F. Tedeschi, G. G. Giusteri, L. Yelash, and M. Lukáčová-Medviďová, A multi-scale method for complex flows of non-Newtonian fluids, Math. Eng., 4(6) (2022), pp. 1–22.
  * [40] S. Yasuda and R. Yamamoto, Multiscale modeling and simulation for polymer melt flows between parallel plates, Phys. Rev. E, 81 (2010), p. 036308.