The incompressible limit of compressible finitely extensible nonlinear bead-spring chain models for dilute polymeric fluids

TL;DR

This paper investigates the transition of dilute polymeric fluid models from compressible to incompressible regimes, demonstrating the limiting behavior of solutions as the Mach number approaches zero.

Contribution

It establishes the rigorous incompressible limit for a class of bead-spring chain models with nonlinear elastic springs in dilute polymeric fluids.

Findings

Weak solutions converge to incompressible models as Mach number tends to zero.

The study confirms the validity of incompressible approximation in dilute polymeric fluids.

Results are obtained for models with FENE spring potentials and slip boundary conditions.

Abstract

We explore the behaviour of global-in-time weak solutions to a class of bead-spring chain models, with finitely extensible nonlinear elastic (FENE) spring potentials, for dilute polymeric fluids. In the models under consideration the solvent is assumed to be a compressible, isentropic, viscous, isothermal Newtonian fluid, confined to a bounded open domain in $\mathbb{R}^3$, and the velocity field is assumed to satisfy a complete slip boundary condition. We show that as the Mach number tends to zero the system is driven to its incompressible counterpart.

Equations

$$\begin{array}{rl}{\partial }_t\rho +{\mathrm{div}}_x(\rho \mathbf{u})& =0\text{mbox}in\Omega \times (0,T],\\ {\partial }_t(\rho \mathbf{u})+{\mathrm{div}}_x(\rho \mathbf{u}\otimes \mathbf{u})+\frac{1}{{\mathrm{Ma}}^2}{\nabla }_{x}p(\rho )& =\frac{1}{\mathrm{Re}}{\mathrm{div}}_x\mathsf{S}(\mathbf{u})+\frac{1}{{\mathrm{Fr}}^2}\rho \mathbf{f}+\frac{1}{\mathrm{Re}}{\mathrm{div}}_x{\tau }_1-\frac{\tilde{\xi }}{{\mathrm{Ma}}^2}{\nabla }_x{\varrho }^2\text{mbox}in\Omega \times (0,T],\\ {\partial }_t\psi +{\mathrm{div}}_x(\mathbf{u}\psi )+\sum _{i=1}^K{\mathrm{div}}_{{\mathbf{q}}_i}(({\nabla }_{x}u){\mathbf{q}}_i\psi )& =\delta {\Delta }_x\psi +\frac{1}{4\mathrm{De}}\sum _{i=1}^K\sum _{j=1}^K{A}_{ij}{\mathrm{div}}_{{\mathbf{q}}_i}\left(M{\nabla }_{{\mathbf{q}}_j}\left(\frac{\psi }{M}\right)\right)\text{mbox}in\Omega \times D\times (0,T],\end{array}$$

$$\begin{array}{rl}\rho (\cdot ,0)& ={\rho }_0(\cdot )\ge 0\text{mbox}in\Omega ,\\ (\rho \mathbf{u})(\cdot ,0)& =({\rho }_0{\mathbf{u}}_0)(\cdot )\text{mbox}in\Omega ,\\ \psi (\cdot ,\cdot ,0)& ={\psi }_0(\cdot ,\cdot )\ge 0\text{mbox}on\Omega \times D.\end{array}$$

$$\mathsf{A}=[A_{ij}{]}_{i,j=1}^K,\text{mbox}asymmetricpositivedefinitematrix,withsmallesteigenvalue{a}_0\in {\mathbb{R}}_{>0},$$

$$\boldsymbol{\mathsf{S}}(\boldsymbol{u}):=\mu^{S}\bigg{(}\boldsymbol{\mathsf{D}}\boldsymbol{u}-\frac{1}{3}({{\mathrm{div}}_{x}}\,\boldsymbol{u})\,\boldsymbol{\mathsf{I}}\bigg{)}+\mu^{B}({{\mathrm{div}}_{x}}\,\boldsymbol{u})\,\boldsymbol{\mathsf{I}},$$

$$p(\rho ):=c_p{\rho }^{\gamma },$$

$$D:=B(0,b_1^{\frac{1}{2}})\times \cdots \times B(0,b_K^{\frac{1}{2}}),\text{mbox}where{b}_i>0,i=1,\dots ,K,K\ge 1,$$

$$\boldsymbol{\mathsf{\tau}}(\psi)(\boldsymbol{x},t):=\frac{1}{\rm Re}\boldsymbol{\mathsf{\tau}}_{1}(\psi)(\boldsymbol{x},t)-\frac{1}{{\rm Ma}^{2}}\bigg{(}\int_{D\times D}\gamma(\boldsymbol{q},\boldsymbol{q}^{\prime})\psi(\boldsymbol{x},\boldsymbol{q},t)\psi(\boldsymbol{x},\boldsymbol{q}^{\prime},t){\rm\,d}\boldsymbol{q}{\rm\,d}\boldsymbol{q}^{\prime}\bigg{)}\,\boldsymbol{\mathsf{I}}=\frac{1}{\rm Re}\boldsymbol{\mathsf{\tau}}_{1}(\boldsymbol{x},t)-\frac{\tilde{\xi}}{{\rm Ma}^{2}}\varrho^{2}(\boldsymbol{x},t)\,\boldsymbol{\mathsf{I}},$$

$$\gamma (\mathbf{q},{\mathbf{q}}^{\prime })=\tilde{\xi },\text{mbox}where\tilde{\xi }\in {\mathbb{R}}_{\ge 0},$$

$$\varrho (\mathbf{x},t):={\int }_D\psi (\mathbf{x},\mathbf{q},t)\mathrm{d}\mathbf{q},(\mathbf{x},t)\in \Omega \times [0,T].$$

$${\tau }_1(\psi ):=\frac{1-\beta }{\mathrm{De}}\left(\sum _{i=1}^K{\mathsf{C}}_i(\psi )-(K+1){\int }_D\psi \mathrm{d}\mathbf{q}\mathsf{I}\right),$$

$$\boldsymbol{\mathsf{C}}_{i}(\psi)(\boldsymbol{x},t):=\int_{D}\psi(\boldsymbol{x},\boldsymbol{q},t)\,U^{\prime}_{i}\bigg{(}\frac{1}{2}|\boldsymbol{q}_{i}|^{2}\bigg{)}\boldsymbol{q}_{i}\boldsymbol{q}^{\rm T}_{i}{\rm\,d}\boldsymbol{q},\quad i=1,\dots,K,$$

$$\boldsymbol{F}_{i}(\boldsymbol{q}_{i}):=U^{\prime}_{i}\bigg{(}\frac{1}{2}|\boldsymbol{q}_{i}|^{2}\bigg{)}\boldsymbol{q}_{i},\quad i=1,\dots,K,$$

$$M_i({\mathbf{q}}_i):=\frac{1}{Z_i}{\mathrm{e}}^{-U_i(\frac{1}{2}|{\mathbf{q}}_i{|}^2)},Z_i:={\int }_{D_i}{\mathrm{e}}^{-U_i(\frac{1}{2}|{\mathbf{q}}_i{|}^2)}\mathrm{d}{\mathbf{q}}_{\mathbf{i}},i=1,\dots ,K.$$

$$M(\mathbf{q}):=\prod _{i=1}^K{M}_i({\mathbf{q}}_{\mathbf{i}})\text{mbox}forall\mathbf{q}=({\mathbf{q}}_1^{\mathrm{T}},\dots ,{\mathbf{q}}_K^{\mathrm{T}}{)}^{\mathrm{T}}\in D=D_1\times \cdots \times D_K.$$

$$M(\boldsymbol{q})\nabla_{q_{i}}[M(\boldsymbol{q})]^{-1}=-[M(\boldsymbol{q})]^{-1}\nabla_{q_{i}}M(\boldsymbol{q})=U_{i}^{\prime}\bigg{(}\frac{1}{2}|\boldsymbol{q}_{i}|^{2}\bigg{)}\boldsymbol{q}_{i},$$

$${\int }_{D}M(\mathbf{q})\mathrm{d}\mathbf{q}=1.$$

$$c_{i1}[\mathrm{dist}({\mathbf{q}}_i,\partial D_i){]}^{{\theta }_i}\le M_i({\mathbf{q}}_i)\le c_{i2}[\mathrm{dist}({\mathbf{q}}_i,\partial D_i){]}^{{\theta }_i}\text{mbox}forall{\mathbf{q}}_i\in D_i,$$

$$c_{i3}\leq{\rm dist}(\boldsymbol{q}_{i},\partial D_{i})\,U^{\prime}_{i}\bigg{(}\frac{1}{2}|\boldsymbol{q}_{i}|^{2}\bigg{)}\leq c_{i4}\quad\mbox{ for }\boldsymbol{q}_{i}\in D_{i}\mbox{ with }i=1,\dots,K.$$

$$\int_{D_{i}}\left[1+\bigg{[}U_{i}\bigg{(}\frac{1}{2}|\boldsymbol{q}_{i}|^{2}\bigg{)}\bigg{]}^{2}+\bigg{[}U^{\prime}_{i}\bigg{(}\frac{1}{2}|\boldsymbol{q}_{i}|^{2}\bigg{)}\bigg{]}^{2}\right]M_{i}(\boldsymbol{q}_{i}){\rm\,d}\boldsymbol{q}_{i}<\infty,\quad i=1,\dots,K.$$

$$\mathrm{Re}:=\frac{{\rho }_0{L}_0{U}_0}{{\eta }_s+{\eta }_p},$$

$$\mathrm{De}:=\frac{\lambda U_0}{L_0}=\frac{{\zeta }_0{U}_0}{4H{L}_0},$$

$$\mathrm{Ma}:=\frac{U_0}{\sqrt{p_0/{\rho }_0}},$$

$$\mathrm{Fr}:=\frac{U_0}{\sqrt{L_0/f_0}},$$

$$\delta :=\frac{(l_0/L_0{)}^2}{(4(K+1){\zeta }_0{U}_0)/(4H{L}_0)}=\frac{(l_0/L_0{)}^2}{4(K+1)\mathrm{De}}$$

$$\mathbf{u}\cdot \mathbf{n}{|}_{\partial \Omega }=0,[\mathsf{S}(\mathbf{u})+{\tau }_1]\mathbf{n}\times \mathbf{n}{|}_{\partial \Omega }=0.$$

$$\left[\frac{1}{4\mathrm{De}}\sum _{j=1}^K{A}_{ij}M{\nabla }_{q_j}\left(\frac{\psi }{M}\right)-({\nabla }_x\mathbf{u}){\mathbf{q}}_i\psi \right]\cdot \frac{{\mathbf{q}}_i}{|{\mathbf{q}}_i|}=0\text{mbox}on\Omega \times \partial {\overline{D}}_i\times (0,T],\text{mbox}fori=1,\dots ,K,$$

$$\delta {\nabla }_x\psi \cdot \mathbf{n}=0\text{mbox}on\partial \Omega \times D\times (0,T],$$

$$\psi (\cdot ,\cdot ,0)={\psi }_0(\cdot ,\cdot )\ge 0\text{mbox}on\Omega \times D,$$

$$\begin{split}\partial_{t}\rho+{{\mathrm{div}}_{x}}\,(\rho\boldsymbol{u})&=0\quad\mbox{ in }\Omega\times(0,T],\\ \partial_{t}(\rho\boldsymbol{u})+{{\mathrm{div}}_{x}}\,(\rho\boldsymbol{u}\otimes\boldsymbol{u})+\frac{1}{\varepsilon^{2}}\nabla_{x}p(\rho)&={{\mathrm{div}}_{x}}\,\boldsymbol{\mathsf{S}}(\boldsymbol{u})+{{\mathrm{div}}_{x}}\,\boldsymbol{\mathsf{\tau}}_{1}-{\bar{\xi}}\nabla_{x}\varrho^{2}\quad\mbox{ in }\Omega\times(0,T],\\ \partial_{t}\psi+{{\mathrm{div}}_{x}}\,(\boldsymbol{u}\psi)+\sum_{i=1}^{K}{\rm div}_{\boldsymbol{q}_{i}}((\nabla_{x}\boldsymbol{u})\boldsymbol{q}_{i}\psi)&=\delta\Delta_{x}\psi+\frac{1}{4}\sum_{i=1}^{K}\sum_{j=1}^{K}A_{ij}{\rm div}_{\boldsymbol{q}_{i}}\left(M\nabla_{\boldsymbol{q}_{j}}\bigg{(}\frac{\psi}{M}\bigg{)}\right)\quad\mbox{ in }\Omega\times D\times(0,T],\\ \mbox{ with }\varrho(\boldsymbol{x},t)&:=\int_{D}\psi(\boldsymbol{x},\boldsymbol{q},t){\rm\,d}\boldsymbol{q},\quad(\boldsymbol{x},t)\in\Omega\times[0,T]\end{split}$$

$$\Vert {\mathbf{u}}_{0,\epsilon }{\Vert }_{L^2(\Omega )}\le c,$$

Full text

Keywords: Navier–Stokes–Fokker–Planck system, Mach number, compressible fluid, incompressible fluid, weak solutions, singular limit

1 Introduction

The purpose of this paper is to explore the singular limit of finite-energy global-in-time weak solutions, as the Mach number tends to zero, to a class of kinetic models of dilute polymeric fluids with noninteracting polymer chains, where the solvent is a compressible, isentropic, viscous isothermal Newtonian fluid, and the polymer molecules suspended in the solvent are idealised as linear bead-spring chains with finitely extensible nonlinear elastic (FENE) spring potentials. The existence of finite-energy global weak solutions to this class of models, which are coupled compressible Navier–Stokes–Fokker–Planck systems, has been shown in [3] in the case of a nonslip boundary condition for the velocity of the solvent. With minor modifications, the existence proof presented in [3] extends to the case of a complete slip (also referred to as Navier slip) boundary condition and, for the sake of simplicity of the exposition, it is the latter complete slip boundary condition that we shall assume hereafter. By performing a rigorous passage to the limit we show that, as the Mach number tends to zero, the compressible Navier–Stokes–Fokker–Planck system is driven to its incompressible limit, which is the coupled FENE-type bead-spring chain model arising from the kinetic theory of dilute solutions of polymeric liquids with noninteracting polymer chains. The limiting model involves the unsteady incompressible Navier–Stokes equations, with an elastic extra-stress tensor appearing on the right-hand side of the momentum equation. The elastic extra-stress tensor stems from the random movement of the polymer chains and is defined by the Kramers expression through the associated probability density function that satisfies a Fokker–Planck-type parabolic equation. The existence of global-in-time weak solutions to the limiting incompressible model was shown in [1] in the case of a nonslip boundary condition for the velocity of the solvent. Again, the existence proof presented in [1] extends to the case of a complete slip boundary condition. Our main result is therefore that the incompressible Navier–Stokes–Fokker–Planck system considered in [1], but with a complete slip boundary condition, is the singular limit, as the Mach number tends to zero, of the compressible model considered in [3], with a complete slip boundary condition.

Our proofs rely on combining the theoretical machinery developed in the monograph of Feireisl and Novotný [7], devoted to the rigorous analysis of various asymptotic limits of the compressible Navier–Stokes–Fourier system, with the analytical techniques involved in the existence proofs in [1], [3], and [9] for coupled incompressible and compressible Navier–Stokes–Fokker–Planck systems.

For fully macroscopic models of compressible non-Newtonian fluids there are relatively few results in the literature concerning the incompressible limit. In [12], Lei and Zhou studied the existence of global classical solutions of the incompressible Oldroyd-B system on the two-dimensional torus. The Oldroyd-B model arises in the description of viscoelastic fluid flow and consists of a coupling of the Navier–Stokes equations with a system of first-order partial differential equations for the elastic extra stress tensor. Using well-posedness and stability results for the compressible case they pass to the incompressible limit. This strategy is based on dispersive energy estimates for the compressible Oldroyd-B model. The outcome of their analysis is a global existence/uniqueness result for solutions to the two-dimensional incompressible Oldroyd-B model for small data and a justification of the incompressible Oldroyd-B model as a limit of the slightly compressible Oldroyd-B model. In [11], Lei studied the incompressible limit problem for the compressible Oldroyd-B model on a torus, and showed that compressible flows with well-prepared initial data converge to incompressible flows as the Mach number tends to zero. The case of bounded domains $\Omega\subset\mathbb{R}^{3}$ was explored by Guillopé, Salloum and Talhouk [10], where they proved the local and global existence of strong solutions to the weakly compressible Oldroyd-B model, and, with the aid of conformal coordinates, they also studied the incompressible limit problem with well-prepared initial data. Subsequently Fang and Zi [6] proved, using the techniques of Danchin [4] in scale-invariant Besov spaces on $\mathbb{R}^{d}$, that solutions to the compressible Oldroyd-B model with ill-prepared initial data (in which case strong time-oscillations of solutions need to be considered) converge to those of the corresponding incompressible Oldroyd-B model as the Mach number tends to zero. More recently, Ren and Ou [13] showed the local existence of strong solutions to an Oldroyd-B model for incompressible viscoelastic fluids in a bounded domain $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, via the incompressible limit: the main idea of the paper was to derive uniform estimates with respect to the Mach number for a linearised compressible Oldroyd-B system.

To the best of our knowledge no rigorous asymptotic results are available in the literature for the incompressible limit of global weak solutions (with large initial data) to kinetic models of dilute compressible polymeric fluids involving the (compressible) Navier–Stokes–Fokker–Planck system, and it is the study of this question that is our objective here. The paper is structured as follows: in the next two subsections of this introductory section we perform a nondimensionalisation of the compressible Navier–Stokes–Fokker–Planck system and we state our main result; Section 2 is then devoted to the proof of our main result. We end the paper, in Section 3, with concluding remarks concerning possible extensions.

1.1 Dimensionless form of the Navier–Stokes–Fokker–Planck system

We consider the following system of nonlinear partial differential equations consisting of the equations of continuity and balance of linear momentum, having the form of the compressible Navier–Stokes equations in which the elastic extra-stress tensor (the polymeric part of the Cauchy stress tensor) appears as a source term in the balance of linear momentum equation. In particular, for a given $T\in{\mathbb{R}}_{>0}$ and a bounded open domain $\Omega\subset{\mathbb{R}}^{3}$ with $\partial\Omega\in C^{2,\alpha}$, $\alpha\in(0,1)$, we denote by $\rho:\Omega\times[0,T]\to{\mathbb{R}}$ the density of the solvent, by $\boldsymbol{u}:\overline{\Omega}\times[0,T]\to{\mathbb{R}}^{3}$ the velocity of the solvent, which satisfies the conservation of linear momentum equation for a viscous compressible fluid, and $\psi:\overline{\Omega}\times\overline{D}\times[0,T]\to{\mathbb{R}}_{\geq 0}$ is the probability density function associated with the random movement of the polymer chains and satisfying a Fokker–Planck equation to be stated below. The configuration space domain $D:=D_{1}\times\cdots\times D_{K}\subset\mathbb{R}^{3K}$ is the Cartesian product of $K$ balanced convex open sets $D_{i}\subset\mathbb{R}^{d}$, $i=1,\dots,K$, with $q_{i}\in D_{i}$ denoting the orientation vector of the $i$-th spring in the bead-spring chain representing an idealisation of a polymer chain suspended in the solvent. The coupled Navier–Stokes–Fokker–Planck system under consideration has the following form:

[TABLE]

with the initial conditions

[TABLE]

In the above system $\mathrm{Ma}$ is the Mach number, $\mathrm{Re}$ is the Reynolds number, $\mathrm{Fr}$ is the Froude number, $\mathrm{De}$ is the Deborah number, $\tilde{\xi}$ and $\delta$ are positive real numbers, $A_{ij}$, $i,j=1,\dots,K$, in the Fokker–Planck equation (1.1)3, are the entries of

[TABLE]

called the Rouse matrix, or connectivity matrix (e.g. $\boldsymbol{\mathsf{A}}={\rm tridiag}[-1,2,-1]$ in the case of a topologically linear bead-spring chain). The function $p=p(\rho)$ is the pressure, $\boldsymbol{\mathsf{S}}$ is the Newtonian part of the viscous stress tensor, defined by

[TABLE]

where $\boldsymbol{\mathsf{I}}$ is the $3\times 3$ identity tensor, $\boldsymbol{\mathsf{D}}\boldsymbol{v}:=\frac{1}{2}(\nabla_{x}\boldsymbol{v}+(\nabla_{x}\boldsymbol{v})^{\rm T})$ is the rate of strain tensor, $\mu^{S}$ and $\mu^{B}$ are positive constants referred to as the shear viscosity and the bulk viscosity, respectively, and the function $M=M(\boldsymbol{q})$ is the Maxwellian, to be definedhttps://www.overleaf.com/project/5bd0b4be1ae53d317f17aefd below, associated with the model. In addition, the pressure $p$ is assumed to be related to the density $\rho$ of the solvent by the isentropic equation of state:

[TABLE]

where $c_{p}\in{\mathbb{R}}_{>0}$ and the constant $\gamma>3/2$. The above system is supplemented with boundary conditions, which will be stated below.

In a bead-spring chain model consisting of $K+1$ beads, linearly coupled with $K$ elastic springs to represent a polymer chain, the extra-stress tensor $\boldsymbol{\mathsf{\tau}}$ is given by a version of the Kramers expression depending on the probability density function $\psi$ of the random conformation vector $\boldsymbol{q}=(\boldsymbol{q}^{\rm T}_{1},\dots,\boldsymbol{q}^{\rm T}_{K})^{\rm T}\in D:=D_{1}\times\cdots\times D_{K}\subset{\mathbb{R}}^{3K}$ of the chain, where $\boldsymbol{q}_{i}$ represents the 3-component conformation/orientation vector of the $i$-th spring in the chain. Typically $D_{i}$ is the whole of ${\mathbb{R}}^{3}$ or a bounded open $3$-dimensional ball centred at the origin $\boldsymbol{0}\in{\mathbb{R}}^{3}$, for each $i=1,\dots,K.$ If $K=1$, then the model is referred to as the dumbbell model. Here we shall concentrate on finitely extensible nonlinear elastic (FENE) bead-spring chain models, with

[TABLE]

and $B(\mathbf{0},b^{\frac{1}{2}}_{i})$ is a bounded open ball in ${\mathbb{R}}^{3}$ of radius $b^{\frac{1}{2}}_{i}$, centred at $\mathbf{0}\in{\mathbb{R}}^{3}$.

On the right-hand side of (1.1)2, the 3-component vector function $\boldsymbol{f}$ is the nondimensional density of body forces and the elastic extra-stress tensor is of the form:

[TABLE]

where $\gamma:D\times D\to{\mathbb{R}}_{\geq 0}$ is a smooth, time-independent, $x$-independent and $\psi$-independent interaction kernel, which we take here for the sake of simplicity to be

[TABLE]

and the polymer number density is defined by

[TABLE]

Moreover, $\boldsymbol{\mathsf{\tau}}_{1}(\psi)$ appearing in (1.7) is the Kramers expression:

[TABLE]

where $(1-\beta):=\eta_{p}/(\eta_{s}+\eta_{p})$, with $\eta_{p}$ signifying the polymeric viscosity and $\eta_{s}$ the viscosity of the solvent. In the above

[TABLE]

where $U_{i}$ is the $i$-th spring potential, which we shall now fix and we shall also define the associated Maxwellian.

Let ${\mathcal{O}}_{i}\subset[0,\infty)$ denote the image of $D_{i}$ ($0\in{\mathcal{O}}_{i}$) under the mapping $\boldsymbol{q}_{i}\in D_{i}\mapsto\frac{1}{2}|\boldsymbol{q}_{i}|^{2}$ and let us consider the spring potential $U_{i}\in C^{1}({\mathcal{O}}_{i};{\mathbb{R}}_{\geq 0}),$ $i=1,\dots,K$. In the case of the FENE spring potential, $\mathcal{O}_{i}=[0,b_{i})$, $i=1,\dots,K$. We shall suppose that $U_{i}(0)=0$ and that $U_{i}$ is unbounded on ${\mathcal{O}}_{i}$ for $i=1,\dots,K$ (i.e., $\lim_{s\rightarrow b_{i}{-}}U_{i}(s)=+\infty$). The elastic spring-force $\boldsymbol{F}_{i}:D_{i}\subseteq{\mathbb{R}}^{3}\to{\mathbb{R}}^{3}$ of the $i$-th spring in the chain is then defined by

[TABLE]

and the partial Maxwellian $M_{i}$, associated with the spring potential $U_{i}$, is defined by

[TABLE]

The total Maxwellian $M_{i}$ in the model is then

[TABLE]

A straightforward calculation yields that, for $i=1,\dots,K$,

[TABLE]

and

[TABLE]

We shall suppose that for $i=1,\dots,K$ there exist constants $c_{ij}>0,$ $j=1,\dots,4$, and $\theta_{i}>1$ such that the spring potential $U_{i}$ and the associated Maxwellian $M_{i}$ satisfy

[TABLE]

[TABLE]

It then follows from the above that

[TABLE]

In the system (1.1) the Reynolds number Re and the Deborah number De are defined, respectively, by

[TABLE]

[TABLE]

with $\mathrm{De}$ characterising the elastic relaxation properties, $\lambda>0$ being the characteristic relaxation time, $U_{0}$ the characteristic speed of the fluid, $L_{0}$ the characteristic macroscopic lengthscale (e.g. the diameter of the flow domain $\Omega$), $H>0$ the spring constant, $\zeta_{0}>0$ the characteristic drag coefficient, and $\rho_{0}$ is the characteristic density; the Mach number and the Froude number are defined, respectively, by

[TABLE]

[TABLE]

where $p_{0}$ is the characteristic pressure, $f_{0}$ signifies the characteristic density of body forces; and

[TABLE]

is the polymeric diffusion coefficient, where $l_{0}$ is the characteristic microscopic lengthscale (e.g. the length of a typical polymer chain suspended in the solvent).

We supplement the system of partial differential equations under consideration with a complete slip boundary condition for the velocity field of the solvent, that is,

[TABLE]

Moreover, we impose the following boundary and initial conditions on solutions of the Fokker–Planck equation:

[TABLE]

[TABLE]

[TABLE]

where $\partial\overline{D}_{i}:=D_{1}\times\cdots\times D_{i-1}\times\partial D_{i}\times D_{i+1}\times\cdots\times D_{K},$ $\boldsymbol{q}_{i}$ is normal to $\partial D_{i}$, as $D_{i}$ is a bounded ball centred at the origin, and $\boldsymbol{n}$ is the unit outward normal vector to $\partial\Omega.$

Remark 1.1

The Navier–Stokes–Fokker–Planck system stated above is in dimensionless form. It is arrived at by introducing new variables of the form $X^{*}=\frac{X}{X_{0}}$ into the dimensionful version of the system, where $X_{0}>0$ is the characteristic value of $X$; for example, $L_{0}>0$ is a reference length, $T_{0}>0$ is a reference time, $U_{0}>0$ is a reference velocity, $\rho_{0}>0$ is a reference density of the solvent, proceeding similarly with the characteristic values of the other physical entities entering into the equations: $p_{0}>0$, $\mu^{S}_{0}>0$, $\mu^{B}_{0}>0$, $f_{0}>0$, $\eta^{p}>0$. In the above scaling the polymeric pressure term $p_{p}=\xi\varrho^{2}$ is treated as a part of a fluid pressure, and it is therefore rescaled by the reference value $p_{0}$.

In this paper we concentrate on the following choices of the various dimensionless numbers:

• •

${\rm Ma}=\varepsilon\ll 1$, i.e., the characteristic velocity is dominated by the speed of sound;

• •

${\mathrm{De}}=1$ and ${\rm Re}=1$;

• •

One can choose ${\rm Fr}=\sqrt{\varepsilon}$ (the case of low stratification), but for the sake of simplicity we shall not consider the influence of external forces in the momentum equation in this work, and will therefore set $\boldsymbol{f}=\boldsymbol{0}$;

• •

For the dimensionless coefficient $\tilde{\xi}$ we shall assume here that $\tilde{\xi}=\bar{\xi}\,{\rm Ma}^{2}=\bar{\xi}\,\varepsilon^{2}$, with some $\bar{\xi}>0$.

The system (1.1) therefore takes the following form: for a given $T\in{\mathbb{R}}_{>0}$ and a bounded open domain $\Omega\subset{\mathbb{R}}^{3}$, with $\partial\Omega\in C^{2,\alpha}$, $\alpha\in(0,1)$, we consider the following system of partial differential equations:

[TABLE]

for each $\varepsilon\in(0,1),$ supplemented with the boundary conditions (1.22)–(1.24), the initial conditions (1.2), and the initial data appearing in (1.2) assumed to satisfy the following properties:

[TABLE]

[TABLE]

[TABLE]

where $\overline{\varrho}$ is a static constant density satisfying the equilibrium state (static state) equation

[TABLE]

In the above we have used the following notation:

[TABLE]

Clearly, $\mathcal{F}^{\prime}(s)=\log(s)$ and $\mathcal{F}^{\prime\prime}(s)=\frac{1}{s}$ for $s>0$.

Remark 1.2** (Remark 1.3 in [3])**

Defining the polymer number density by (1.8), formally integrating the Fokker–Planck equation in (1.1) over $D$, and using the boundary condition (1.23), we infer that

[TABLE]

with $\varrho(x,0)=\varrho_{0}$ in $\Omega$ and $\frac{\partial\varrho}{\partial n}|_{\partial\Omega}=0$. Hence by the boundary and initial conditions (1.24), (1.25) we get that

[TABLE]

If ${{\mathrm{div}}_{x}}\,\boldsymbol{u}=0$ and $\varrho(\cdot,0)$ is constant, then $\varrho(\boldsymbol{x},t)$ is constant (and equal to $\varrho(\boldsymbol{x},0)\equiv\mathrm{Const.}>0$) for all $(\boldsymbol{x},t)\in\Omega\times(0,T]$, and thus

[TABLE]

In other words, the polymer number density is constant.

Our aim is to show by rigorous analysis that, as the Mach number converges to zero, the primitive system (NSFPε) converges to the following incompressible Navier–Stokes–Fokker–Planck system:

[TABLE]

where $\Pi$ is a new pressure function, $\boldsymbol{U}$ is the limiting divergence-free velocity field, the density of the solvent $\bar{\rho}$ is constant,

[TABLE]

is the Kramers expression for the polymeric extra stress tensor in the incompressible limit, with, again, $(1-\beta):=\eta_{p}/(\eta_{s}+\eta_{p})$, and

[TABLE]

Remark 1.3

One can formally arrive at the system (1.32) by inserting in the primitive system the expansions:

[TABLE]

and neglecting all terms of order $\varepsilon$ and higher. Our goal here is to deduce (1.32) from the primitive system (NSFPε) through a rigorous passage to the limit, as $\varepsilon\rightarrow 0$ (i.e, as $\mathrm{Ma}\rightarrow 0$).

1.2 Main result

For $r\in[1,\infty)$, let $L^{r}_{M}(\Omega\times D)$ denote the Maxwellian-weighted Lebesgue space equipped with the norm

[TABLE]

and let $H^{1}_{M}(\Omega\times D)$ denote the Maxwellian-weighted Sobolev space with the norm

[TABLE]

We also define

[TABLE]

Definition 1.1

We call a triple $(\rho_{\varepsilon},\boldsymbol{u}_{\varepsilon},\widehat{\psi}_{\varepsilon})$ a weak solution to the system (NSFPε) provided that:

The functions $\rho_{\varepsilon}$, $\boldsymbol{u}_{\varepsilon}$, $\widehat{\psi}_{\varepsilon}$ satisfy the following regularity properties:

[TABLE]

[TABLE]

where $v\in[1,\infty)$ and $s>1+\frac{3}{2}(K+1)$, with finite relative entropy and Fisher information,

[TABLE]

[TABLE]

• 2)

In addition,

[TABLE]

and

[TABLE]

• 3)

Moreover the following relations are satisfied:

[TABLE]

for some $s\in\left(\left(\frac{6\gamma}{6+\gamma}\right)^{\prime},\infty\right]$, where $\left(\frac{6\gamma}{6+\gamma}\right)^{\prime}$ is the Hölder conjugate of $\frac{6\gamma}{6+\gamma}$, and with $\rho_{\varepsilon}(\cdot,0)=\rho_{0,\varepsilon}$,

[TABLE]

with $(\rho\boldsymbol{u})(\cdot,0)=(\rho_{0}\boldsymbol{u}_{0})(\cdot),$ $\vartheta(\gamma):=\frac{2\gamma-3}{3}$ for $\frac{3}{2}<\gamma\leq 4$ and $\vartheta(\gamma):=\frac{5}{12}\gamma$ for $4\leq\gamma$; $r:=\max\left\{4,\frac{6\gamma}{2\gamma-3}\right\}$, and

[TABLE]

with $\widehat{\psi}_{\varepsilon}(\cdot,0)=\widehat{\psi}_{0,\varepsilon}(\cdot)$ and $s>1+\frac{3}{2}(K+1).$

Let $\boldsymbol{H}$ denote the standard Helmholtz projection onto the space of solenoidal (divergence-free) functions; that is,

[TABLE]

where $\Phi\in H^{1}(\Omega)/{\mathbb{R}}$ is the unique solution of the problem

[TABLE]

Definition 1.2** (Weak solution to the limiting system)**

We shall say that the couple $(\boldsymbol{U},\Psi)$ is a weak solution to the problem (1.32), provided that the following properties hold:

$\boldsymbol{U}\in L^{\infty}(0,T;L^{2}(\Omega)^{3})\cap L^{2}(0,T;H^{1}(\Omega)^{3})$, $\boldsymbol{U}\cdot\boldsymbol{n}|_{\partial\Omega}=\boldsymbol{0}$ for a.e. $t\in(0,T)$, ${{\mathrm{div}}_{x}}\,\boldsymbol{U}=0$ for a.e. $(\boldsymbol{x},t)\in\Omega\times(0,T)$ and $\widehat{\Psi}\in L^{v}(0,T;Z_{1})\cap H^{1}(0,T;M^{-1}(H^{s}(\Omega\times D))^{\prime}),$ where $v\in[1,\infty)$ and $s>1+\frac{3}{2}(K+1)$, with finite relative entropy and Fisher information, i.e.,

[TABLE]

[TABLE]

• 2)

Moreover the following relations are satisfied:

[TABLE]

where $\boldsymbol{U}_{0}\in L^{2}(\Omega\times(0,T))^{3}$, $\vartheta(\gamma)=\frac{2\gamma-3}{3}$ for $\frac{3}{2}<\gamma\leq 4$ and $\vartheta(\gamma)=\frac{5}{12}\gamma$ for $4\leq\gamma$; $r=\max\left\{4,\frac{6\gamma}{2\gamma-3}\right\}$, and

[TABLE]

with $\widehat{\Psi}(\cdot,0)=\widehat{\Psi}_{0}(\cdot)$ and $s>1+\frac{3}{2}(K+1).$ Here $\widehat{\Psi}_{0}\in L^{1}_{M}(\Omega\times D;\mathbb{R}_{\geq 0})$ and $\mathcal{F}(\widehat{\Psi}_{0})\in L^{1}_{M}(\Omega\times D;\mathbb{R}_{\geq 0})$, and $\int_{D}M\widehat{\Psi}_{0}\in L^{\infty}(\Omega;\mathbb{R}_{\geq 0})$.

Our main result is then the following theorem.

Theorem 1.1** (Main theorem)**

Suppose that $\Omega\subset{\mathbb{R}}^{3}$ is a bounded open domain with $\partial\Omega\in C^{2,\alpha}$, with some $\alpha\in(0,1)$. Let $p$ satisfy (1.5) with $\gamma>\frac{3}{2}$, let the stress tensor $\boldsymbol{\mathsf{S}}$ satisfy (1.4), and let the Rouse matrix satisfy (1.3). Let $(\rho_{\varepsilon},\,\boldsymbol{u}_{\varepsilon},\,\psi_{\varepsilon})$ be a weak solution triple to the system (NSFPε) emanating from the initial data satisfying (1.26)–(1.28) and subject to the boundary conditions (1.22)–(1.25). Then, by passing to suitable subsequences if necessary, we have that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Moreover, the couple $(\boldsymbol{U},\Psi)$ is a weak solution to (1.32) according to Definition 1.2 with initial datum $\boldsymbol{U}(0,\cdot)=\boldsymbol{H}[\boldsymbol{U}_{0}]$, where $\boldsymbol{U}_{0}$ and $\Psi_{0}$ are weak limits of (1.26) and (1.28), and $\widehat{\Psi}_{0}:=\Psi_{0}/M$.

2 Proof of the main result

Before embarking on the proof, a remark is in order concerning our choice of boundary condition for the velocity field.

Remark 2.1

As was noted in the Introduction, we have confined ourselves in this paper to the case of a complete slip boundary condition, which models an acoustically hard boundary. Then, if $\Omega\subset{\mathbb{R}}^{3}$ is a bounded domain, as is being assumed here, the gradient part of the velocity field associated to acoustic waves may exhibit fast oscillations in time, and the convergence of the velocity field $\boldsymbol{u}_{\varepsilon}$ in the limit of $\varepsilon\rightarrow 0$ is genuinely weak with respect to the temporal variable (see [7]).

In the case of a Dirichlet (no-slip) boundary condition for the velocity field, a boundary layer may appear because of the presence of viscosity and the specific geometrical properties of the boundary. The decay of the acoustic waves and consequently of the velocity field can then be deduced as in [5].

Our starting point is the following result concerning the existence of weak solutions to the primitive system. Its proof for the case of a homogeneous Dirichlet boundary condition for the velocity field can be found in [3, Lemma 6.5, 6.2, Theorem 6.1], and the arguments contained therein can be easily adapted to the case of complete slip boundary condition by replacing the function space of divergence-free three-component vector functions contained in $H^{1}_{0}(\Omega)$ throughout the proof in [3] by the function space of divergence-free three-component vector functions contained in $H^{1}(\Omega)$ with vanishing normal trace on $\partial\Omega$.

Proposition 2.1** (Existence of solutions to the primitive system)**

For each fixed $\varepsilon>0$ there exists a triple $(\rho_{\varepsilon},\boldsymbol{u}_{\varepsilon},\widehat{\psi}_{\varepsilon})$, which is a global weak solution to problem (NSFPε) in the sense of Definition 1.1. Furthermore, the solution triple $(\rho_{\varepsilon},\boldsymbol{u}_{\varepsilon},\widehat{\psi}_{\varepsilon})$ satisfies, for a.e. $t^{\prime}\in(0,T)$, the energy inequality

[TABLE]

2.1 Energy equality

In order to obtain uniform bounds on $(\rho_{\varepsilon},\boldsymbol{u}_{\varepsilon},\widehat{\psi}_{\varepsilon})$ with respect to $\varepsilon$, we shall prove a formal energy equality. Its derivation is in fact the same as in the proof of the existence of weak solutions in Proposition 2.1. For the moment we shall keep all characteristic numbers in their general form and will permit a nonzero density of body forces $\boldsymbol{f}$ so as to state the formal energy equality in its most general form. By taking the $L^{2}(\Omega)$ inner product of the continuity equation first with $\frac{1}{2}|\boldsymbol{u}|^{2}$ and then with $P^{\prime}(\rho)-P^{\prime}(\bar{\rho})\rho$, where $P(\rho)=\frac{p(\rho)}{\gamma-1}$ and where $\bar{\rho}$ is determined by (1.29), and taking the $L^{2}(\Omega)^{3}$ inner product of the momentum equation with $\boldsymbol{u}$ (by noting the boundary conditions and performing partial integration) we deduce that

[TABLE]

Now let us concentrate on the Fokker–Planck equation. Multiplying the Fokker–Planck equation in (1.1) by $\mathcal{F}^{\prime}(\widehat{\psi})$, integrating over $\Omega\times D$ and noticing that $\nabla_{q_{i}}M=-M\boldsymbol{q}_{i}U^{\prime}_{i}\left(\frac{|\boldsymbol{q}_{i}|^{2}}{2}\right)$, we deduce that

[TABLE]

Multiplying (2.3) by $\frac{1-\beta}{{\rm Wi\,Re}}$ and adding to (2.2) we deduce that

[TABLE]

It remains to deal with the last term on the left-hand side of (2.4). To this end we integrate the Fokker–Planck equation over $D$ and we get (1.30). After multiplying (1.30) by $\varrho$ and integrating over $\Omega$ we have

[TABLE]

Multiplying (2.5) by $\frac{\tilde{\xi}}{{\rm Ma}^{2}}$ and substituting the resulting expression into (2.4) gives

[TABLE]

In accordance with our assumptions on the characteristic numbers we shall take $\boldsymbol{f}=\boldsymbol{0}$, ${\rm Ma}=\varepsilon$, ${\mathrm{De}}=1$, ${\rm Re}=1$, and $\tilde{\xi}=\bar{\xi}\varepsilon^{2}$ in the above equality; thereby,

[TABLE]

Consequently, using also (1.3) we obtain

[TABLE]

Remark 2.2

We note that ${P}^{\prime\prime}(\rho)=p^{\prime}(\rho)/\rho>0$ for all $\rho\in(0,\infty)$. Thus ${P}$ is a strictly convex function of $\rho$ on $[0,\infty)$, and ${P}(\rho)-{P}^{\prime}(\bar{\rho})(\rho-\bar{\rho})-{P}(\bar{\rho})\approx c(\rho-\bar{\rho})^{2}$ provided that $\rho$ is close to $\bar{\rho}$.

2.2 Uniform estimates

In order to deduce uniform bounds from (2.7) we introduce, similarly as in [7], the decomposition into essential and residual parts of a measurable function $h$ as

[TABLE]

with $\chi\in C^{\infty}_{c}((0,\infty))$, $0\leq\chi\leq 1,\ \chi=1$ on the set ${\cal{O}}_{{\rm{ess}}}$, where

[TABLE]

The following estimates result from the energy estimate, with a constant $c$ independent of $\varepsilon$. According to (1.26)–(1.28) the expression on the right-hand side of (2.8) is bounded for $\varepsilon\to 0$. Thus we infer that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

By (2.15), together with (2.13), we have that

[TABLE]

The bound (2.13) and Sobolev embedding yield

[TABLE]

Hence, by an interpolation argument, we deduce that

[TABLE]

As $\tilde{P}$ is strictly convex, the relation (2.11) yields

[TABLE]

and, by virtue of (1.5) and (2.12),

[TABLE]

Next, as a direct consequence of (2.9), we have that

[TABLE]

and by Korn’s inequality ([7, Theorem 10.17]) and (2.21), (2.10) gives

[TABLE]

Moreover by (2.9), (2.11), (2.12), and (2.22) we have that

[TABLE]

Furthermore, by using the continuity equation we additionally infer that

[TABLE]

Next, we shall establish the necessary bounds on the extra stress tensor $\boldsymbol{\mathsf{\tau}}_{1}$. We deduce from (1.10), (1.13), and by noting that $M|_{\partial D}=0$, that

[TABLE]

As $\nabla_{{q_{i}}}\varphi=\nabla_{q_{i}}(\sqrt{\varphi}\,)^{2}=2\sqrt{\varphi}\,\nabla_{q_{i}}\sqrt{\varphi}$ for any sufficiently smooth nonnegative function $\varphi$, we have that

[TABLE]

for $r\in[1,2)$. Then, for $s\in[1,2]$ we get that, for any such function $\varphi$,

[TABLE]

where $v=\frac{s}{2-s}$ if $s\in[1,2)$ and $v=\infty$ if $s=2$. We deduce from (2.13), (2.25) and (2.15) that, for $i=1,\dots,K,$

[TABLE]

where $r\in[1,2),$ $s\in[1,2]$, and $v=\frac{s}{2-s}$ if $s\in[1,2)$, and $v=\infty$ if $s=2$. Next by (2.13) and interpolation we have

[TABLE]

where $\nu=\frac{3(v-2)}{2v}$, and $v\in(2,6]$. By (1.9), (2.13), (2.15), (2.17), (2.18) (2.26), (2.29) we deduce that

[TABLE]

[TABLE]

where $c$ is independent of $\varepsilon$.

Now let us show that

[TABLE]

with $s>1+\frac{3}{2}(K+1)$. Testing the Fokker–Planck equation (NSFPε)3 with $\varphi\in L^{2}(0,T;W^{1,\infty}(\Omega\times D))$ and since $\nabla_{x}\psi=2\sqrt{\psi}\,\nabla_{x}\sqrt{\psi}$ for sufficiently smooth $\psi$, we infer that

[TABLE]

The last inequality holds by (2.13), (2.15), (2.22) and as $H^{s}(\Omega\times D)\subset W^{1,\infty}(\Omega\times D)$ thanks to our assumption that $s>1+\frac{3}{2}(K+1)$. Therefore (2.31) holds.

Finally, we note that by choosing $\varphi=1$ in (1.44) and by the choice of the initial datum $\widehat{\psi}_{0,\varepsilon}$ (cf. (1.28)), we have that

[TABLE]

where again $c$ is independent of $\varepsilon$.

2.3 Convergence as $\varepsilon\rightarrow 0$

Obviously $\rho_{\varepsilon}-\bar{\rho}=[\rho_{\varepsilon}-\bar{\rho}]_{ess}+[\rho_{\varepsilon}-\bar{\rho}]_{res}$; hence, thanks to (2.19), (2.20), we have that

[TABLE]

Next, by (2.22) we have

[TABLE]

(by extracting a subsequence if necessarily). Using (2.9), (2.21) and Sobolev embedding we get

[TABLE]

[TABLE]

Moreover, from (2.9) and (2.20) it follows that

[TABLE]

Then, employing (2.34) and (2.35), we have from (1.42) that

[TABLE]

Since the limiting velocity $\boldsymbol{U}\in L^{2}(0,T;H^{1}(\Omega)^{3})$, we obtain that

[TABLE]

Thus in the limit of $\varepsilon\to 0$ the continuity equation is reduced to the divergence-free condition for the limiting velocity field $\boldsymbol{U}$. This justifies the use of solenoidal test functions for the momentum equation.

As a direct consequence of (2.35) and (2.39) we have

[TABLE]

Thanks to (2.35) and (2.36) we deduce that

[TABLE]

In the above ${\overline{\overline{\rho}\boldsymbol{U}\otimes\boldsymbol{U}}}$ denotes a weak limit of $\{\rho_{\varepsilon}\boldsymbol{u}_{\varepsilon}\otimes\boldsymbol{u}_{\varepsilon}\}_{\varepsilon>0}$ in $L^{2}(0,T;L^{q}(\Omega)^{3\times 3})$ for a certain $q>1$. Although we do not really expect that

[TABLE]

we will show in the next section that

[TABLE]

for any $\boldsymbol{w}\in C^{\infty}_{c}((0,T)\times\overline{\Omega})^{3}$, ${{\mathrm{div}}_{x}}\,\boldsymbol{w}=0$, $\boldsymbol{w}\cdot\boldsymbol{n}|_{\partial\Omega}=0$. The relation (2.42) may be interpreted as an expression of the fact that the difference ${{\mathrm{div}}_{x}}\,(\overline{\bar{\rho}\boldsymbol{U}\otimes\boldsymbol{U}}-{\bar{\rho}\boldsymbol{U}\otimes\boldsymbol{U}})$ is proportional to a gradient and that it can be therefore incorporated into the limiting pressure.

Next we turn our attention to passage to the limit $\varepsilon\rightarrow 0$ in the terms related to the probability density function $\widehat{\psi}_{\varepsilon}$. In particular, we will show the strong convergence of the sequence $\{\widehat{\psi}_{\varepsilon}\}_{\varepsilon>0}$ using Dubinskiǐ’s compactness theorem (cf. [2], for example,) stated in the next lemma.

Lemma 2.2

Let $X_{1}$ be a seminormed set, and let $X_{2}$ and $X_{3}$ be Banach spaces such that $X_{1}\subset\subset X_{2}\subset X_{3}$. Then, for $p_{1},\,p_{2}\in[1,\infty)$ the following embedding is compact

[TABLE]

We shall employ Lemma 2.2 with $X_{2}=L^{1}_{M}(\Omega\times D)$, $X_{3}=M^{-1}H^{s}(\Omega\times D)^{\prime}$,

[TABLE]

$p_{1}=v$ and $p_{2}=2$. We note that $X_{1}\subset\subset X_{2}$ (for the details, see [1, Section 5]). Moreover, $X_{2}\subset X_{3}$. Hence, thanks to (2.31) and (2.16), we have

[TABLE]

We infer also that

[TABLE]

[TABLE]

[TABLE]

By (2.43), (2.13), and (2.30) we obtain that

[TABLE]

We observe that, by (2.13), as $\varepsilon\to 0$ we have that

[TABLE]

Next, we note that

[TABLE]

Indeed, (1.40) and (2.43), together with (2.48) provide (2.49).

By (2.35), (2.37), (2.47) and (1.40), we may pass to the limit $\varepsilon\to 0$ in (1.43) with test functions $\boldsymbol{w}\in C^{1}_{c}([0,T);C^{\infty}(\overline{\Omega})^{3})$, $\boldsymbol{w}\cdot\boldsymbol{n}|_{\partial\Omega}=0$ and ${{\mathrm{div}}_{x}}\,\boldsymbol{w}=0$, and obtain that

[TABLE]

where we have assumed that $\boldsymbol{u}_{\varepsilon,0}\rightharpoonup\boldsymbol{U}_{0}$ in $L^{2}(\Omega)^{3}$. We note that, since in the above formulation (2.50) we have restricted ourselves to divergence-free test functions, the limiting term $\nabla_{x}\varrho^{2}$ does not appear because the gradient of a scalar function can be absorbed into the limiting pressure $\Pi$.

Next let us pass to the limit in (1.44) and finally obtain (1.48). Initially we fix the test function to $\phi\in C([0,T];C^{\infty}(\overline{\Omega\times D}))$. The first term of (1.44) converges to the first term of (1.48) thanks to (2.46). For the second term of (1.44) we observe that

[TABLE]

Hence, by noting that $|\sqrt{c_{1}}-\sqrt{c_{2}}|\leq\sqrt{|c_{1}-c_{2}|}$ for all $c_{1},\,c_{2}\in{\mathbb{R}}_{\geq 0}$ and by (2.16), we have that

[TABLE]

Thus, (2.43) implies that $I_{1}$ converges to zero as $\varepsilon\to 0$. Similarly, because $M^{\frac{1}{2}}\sqrt{\widehat{\Psi}}\,\nabla_{q_{j}}\phi\in L^{2}(0,T;L^{2}(\Omega\times D)^{3})$ for $j=1,\dots,K$, it follows form (2.45) that, as $\varepsilon\to 0$,

[TABLE]

Therefore the second term in (1.44) converges to the second term in (1.48). For the remaining terms in (1.44) we note that

[TABLE]

Next, by (2.22), using Sobolev embedding and (2.33), and recalling (1.40), we obtain that

[TABLE]

and therefore, by invoking (2.43) and (2.48)2, we deduce that $I_{3}$ converges to zero as $\varepsilon\to 0$. Similarly, by noting that (1.40) implies that $\int_{D}M\widehat{\Psi}\boldsymbol{q}_{i}\otimes\nabla_{q_{i}}\phi{\rm\,d}\boldsymbol{q}\in L^{2}(\Omega\times(0,T))^{3\times 3}$, $i=1,\dots,K$, we get by (2.35) that

[TABLE]

Therefore the last term of (1.44) converges to the last term in (1.48). Analogously, the third term of (1.44) converges to the third term in (1.48). In this way we obtain (1.48) for smooth test functions $\phi\in C([0,T];C^{\infty}(\overline{\Omega\times D}))$. In order to extend the class of test functions for the limiting equation we use the density of the function space $C([0,T];C^{\infty}(\overline{\Omega\times D}))$ in $L^{2}(0,T;H^{s}({\Omega\times D}))$, the embedding $H^{s}(\Omega\times D)\subset W^{1,\infty}(\Omega\times D)$, (2.31), (2.43), (2.14), (2.16), (1.40), and (2.35).

2.4 Convergence of the convective part of the momentum equation

In this section we concentrate on proving (2.42). The proof is based on performing a Helmholtz decomposition of the momentum $\rho_{\varepsilon}\boldsymbol{u}_{\varepsilon}=\boldsymbol{H}[\rho_{\varepsilon}\boldsymbol{u}_{\varepsilon}]+\boldsymbol{H}^{\perp}[\rho_{\varepsilon}\boldsymbol{u}_{\varepsilon}]$, the proof of the compactness of the solenoidal part in the decomposition, and the analysis of the acoustic equation governing the time-evolution of the gradient component in the decomposition.

Motivated by the discussion in Section 5.4.2 of [7], we begin by decomposing the convective term as follows:

[TABLE]

Let us emphasise that the solenoidal part of the momentum $\boldsymbol{H}[\rho_{\varepsilon}\boldsymbol{u}_{\varepsilon}]$ does not exhibit oscillations in time and in particular it converges a.e. on the set $\Omega\times(0,T)$. In order to see this, we choose a test function in the momentum equation (1.43) of the form $\boldsymbol{H}[\boldsymbol{w}]$, where $\boldsymbol{w}\in C^{\infty}_{c}(\overline{\Omega}\times[0,T))^{3}$, $\boldsymbol{w}\cdot\boldsymbol{n}|_{\partial\Omega}=0$. Thanks to the uniform estimates obtained in Section 2.2, by noting that the singular term is irrelevant as ${{\mathrm{div}}_{x}}\,\boldsymbol{H}[\boldsymbol{w}]=0$, we deduce that the family $\{\int_{\Omega}\boldsymbol{H}[\rho_{\varepsilon}\boldsymbol{u}_{\varepsilon}]\cdot\boldsymbol{w}\,{\rm d}\boldsymbol{x}\}_{\varepsilon>0}$ forms a bounded and equicontinuous sequence in $C([0,T])$. Therefore by the Arzelà–Ascoli theorem we have that

[TABLE]

for any $\boldsymbol{w}$ as above. By a density argument and noting (2.36), we infer that

[TABLE]

Next, we have that

[TABLE]

The first term on the right-hand side of this equality converges to zero weakly in $L^{1}(\Omega\times(0,T))$ by (2.34), (2.35) and compact embedding. The second term, by (2.35), compact embedding, (2.52), and since $\gamma>\frac{3}{2}$, converges weakly in $L^{1}(\Omega\times(0,T))$ to $\bar{\rho}|\boldsymbol{U}|^{2}$. Hence, as $\boldsymbol{U}=\boldsymbol{H}[\boldsymbol{U}]$, we deduce that

[TABLE]

This, together with (2.52), implies that

[TABLE]

Combining (2.53) and (2.37) we get

[TABLE]

Summarising (2.51), (2.54), and (2.55), we see that in order to prove (2.42) all that is left to be shown is that

[TABLE]

As our a-priori estimates do not provide any bound on the time-derivative of the gradient part of the velocity/momentum, in order to prove (2.56) we follow the strategy from [7], which is based on the observation that possible oscillations in time mutually cancel in the acoustic waves described by means of $\boldsymbol{H}^{\perp}[\rho_{\varepsilon}\boldsymbol{u}_{\varepsilon}]$ governed by the acoustic equation associated with our (NSFPε) system.

2.4.1 Acoustic equation

The analysis of the following acoustic equation will allow us to control the temporal oscillations of the gradient part of the momentum. Let us set, to this end,

[TABLE]

and

[TABLE]

We can then rewrite the continuity equation and the momentum equation with the extra stress tensor on its right-hand side in the form of Lighthill’s acoustic analogy (cf. eqs. (5.137) and (5.140) in [7]):

[TABLE]

supplemented with the boundary condition

[TABLE]

The equations (2.59), together with boundary condition (2.60), are understood in a weak sense; specifically, the integral identity

[TABLE]

holds for any $\varphi\in C_{c}^{\infty}(\overline{\Omega}\times(0,T))$ and

[TABLE]

holds for any $\boldsymbol{\varphi}\in C_{c}^{\infty}(\overline{\Omega}\times(0,T))^{3}$, $\boldsymbol{\varphi}\cdot\boldsymbol{n}|_{\partial\Omega}=0$.

By the uniform estimates (2.9), (2.20) we obtain that

[TABLE]

and by (2.21), (2.38) we deduce that

[TABLE]

Moreover, by combining (2.9), (2.19), (2.20), (2.22), (1.40), (2.30) we deduce that

[TABLE]

where

[TABLE]

Now we can follow the strategy developed in [7]. We reduce the problem to a finite number of modes, which are represented by eigenfunctions of the wave operator in (2.61), (2.62). It then transpires that the nonvanishing oscillatory part of the convective term can be represented as the gradient of a scalar function, and it can be therefore incorporated into the limiting pressure $\Pi$, whereby it is irrelevant in the incompressible limit of $\varepsilon\rightarrow 0$. The details of the spectral analysis of the wave operator appearing in (2.61), (2.62) are contained in the next section.

2.4.2 Spectral analysis of the wave operator

Next, as in Section 5.4.5 of [7], we now focus on the following eigenvalue problem associated with the operator appearing on the left-hand side of (2.61), (2.62):

[TABLE]

and reformulate it as the following homogeneous Neumann problem:

[TABLE]

for which we have a countable system of eigenvalues $0=\Lambda_{0}<\Lambda_{1}\leq\Lambda_{2}\leq\cdots$, with the associated system of real eigenfunctions $\{\zeta_{n}\}_{n=0}^{\infty}$ being a basis of the Hilbert space $L^{2}(\Omega)$. The complex eigenfunctions $\{\boldsymbol{\phi}_{\pm n}\}_{n=0}^{\infty}$ are determined through (2.66) as

[TABLE]

Moreover, we can decompose the Hilbert space $L^{2}(\Omega)^{3}=L^{2}_{\mathrm{div}}(\Omega)^{3}\oplus H_{\perp}(\Omega)$, where $H_{\perp}(\Omega)$ is the closure in the $L^{2}(\Omega)^{3}$ norm of the space spanned by $\{({-i}/{\sqrt{p^{\prime}(\bar{\rho})}})\boldsymbol{\phi}_{n}\}_{n=1}^{\infty}$, and where $L^{2}_{\mathrm{div}}(\Omega)^{3}$ denotes the subspace of divergence-free 3-component vector functions contained in $L^{2}(\Omega)^{3}$, with vanishing normal trace on $\partial\Omega$; that is, the closure in the $L^{2}(\Omega)^{3}$-norm of the set of all $\boldsymbol{\varphi}\in C_{c}^{\infty}(\Omega)^{3}$ s.t. ${{\mathrm{div}}_{x}}\,\boldsymbol{\varphi}=0$ in $\Omega$. In order to reduce the problem to a finite number of modes, we introduce the corresponding orthogonal projection

[TABLE]

and we set

[TABLE]

(since $\boldsymbol{P}_{N}$ commutes with $\boldsymbol{H}^{\perp}$). For any $\boldsymbol{\varphi}\in L^{2}(\Omega)^{3}$ we consider the Fourier coefficients

[TABLE]

and a scale of Hilbert spaces $H_{\alpha,\perp}(\Omega)\subset H_{\perp}(\Omega)$, $\alpha\in[0,1]$, with the norm $\|\cdot\|_{L^{2}_{\alpha,\perp}}$, defined by

[TABLE]

where $\{\Lambda_{n}\}_{n=1}^{\infty}\subset\mathbb{R}_{>0}$ is the family of nonzero eigenvalues associated with (2.67). We note that

[TABLE]

As $H_{0,\perp}=H_{\perp}$ and $H_{1,\perp}\subset L^{6}(\Omega)^{3}$, we deduce by an interpolation argument that there exists an $\bar{\alpha}\in(0,1)$ such that

[TABLE]

Let us return to (2.56) and rewrite it as follows:

[TABLE]

By (2.36), (2.70) and (2.71) we deduce that

[TABLE]

uniformly as $\varepsilon\to 0$, for any fixed $\boldsymbol{w}\in C^{\infty}_{c}(\overline{\Omega}\times[0,T))^{3}$, such that ${{\mathrm{div}}_{x}}\,\boldsymbol{w}=0$ and $\boldsymbol{w}\cdot\boldsymbol{n}|_{\partial\Omega}=0$.

We note that $\Lambda_{N}\to\infty$ as $N\to\infty$. By a duality argument and (2.71) we have also that

[TABLE]

Indeed, this bound is an immediate consequence of the following:

[TABLE]

where the penultimate inequality in this chain of inequalities follows by noting that $a_{n}[\boldsymbol{H}^{\perp}[\boldsymbol{\psi}]]=a_{n}[\boldsymbol{\psi}]$ for all $n\geq 1$, that $a_{n}[\boldsymbol{H}^{\perp}_{N}[\boldsymbol{\psi}]]=a_{n}[\boldsymbol{\psi}]$ for $1\leq n\leq N$ and $a_{n}[\boldsymbol{H}^{\perp}_{N}[\boldsymbol{\psi}]]=0$ for $n>N$, whereby $\|\boldsymbol{H}^{\perp}[\boldsymbol{\psi}]\|_{L^{2}_{1,\perp}}=\|\boldsymbol{\psi}\|_{L^{2}_{1,\perp}}$ and $\|\boldsymbol{H}^{\perp}_{N}[\boldsymbol{\psi}]\|_{L^{2}_{1,\perp}}\leq\|\boldsymbol{\psi}\|_{L^{2}_{1,\perp}}$.

With these arguments in hand, the proof of (2.56) reduces to showing that

[TABLE]

as $\varepsilon\to 0$, or, rather, thanks to (2.34), it suffices to show that

[TABLE]

as $\varepsilon\to 0$; this is done in the next subsection.

2.5 Weak limit in the convective term

In order to prove (2.72) let us choose a test function for (2.61) of the following ‘separated’ form

[TABLE]

where $\zeta_{n}$ and $\Lambda_{n}$ solve the eigenvalue problem (2.67), and for (2.62) as a test function we take

[TABLE]

Consequently, we obtain the following system of ordinary differential equations:

[TABLE]

for $n=1,2,\dots$, where $a_{n}[\boldsymbol{V}_{\varepsilon}]$ are the Fourier coefficients of $\boldsymbol{V}_{\varepsilon}$, defined by (2.69), $b_{n}[r_{\varepsilon}]:=\int_{\Omega}r_{\varepsilon}\zeta_{n}\,{\rm d}\boldsymbol{x},$ and

[TABLE]

Then, in terms of the Helmholtz projection, (2.73) reads as follows:

[TABLE]

where we set $[r_{\varepsilon}]_{N}:=\sum_{n=1}^{N}b_{n}[r_{\varepsilon}]\zeta_{n}$, and thanks to (2.74) we get that $\|L_{\varepsilon,N}\|_{L^{1}(\Omega\times(0,T))}\leq c$ for any fixed $N\geq 1$. Notice that (2.75) is satisfied in a strong sense, since $[r_{\varepsilon}]_{N}$ and $\boldsymbol{H}^{\perp}[\rho_{\varepsilon}\boldsymbol{u}_{\varepsilon}]$ are regular enough. By introducing the potential $\Phi_{\varepsilon,N}$ via

[TABLE]

we can reformulate (2.72) as

[TABLE]

where $\boldsymbol{w}\in C^{\infty}_{c}(\overline{\Omega}\times(0,T))^{3}$, ${{\mathrm{div}}_{x}}\,\boldsymbol{w}=0$, $\boldsymbol{w}\cdot\boldsymbol{n}|_{\partial\Omega}=0.$ Then by (2.75) we obtain that

[TABLE]

Therefore, thanks to (2.74), the right-hand side of (2.77) converges to zero as $\varepsilon\to 0$ for any fixed $\boldsymbol{w}$ as above. This step completes the proof of (2.72) and consequently we may replace $\overline{\bar{\rho}\boldsymbol{U}\otimes\boldsymbol{U}}$ by ${\bar{\rho}\boldsymbol{U}\otimes\boldsymbol{U}}$ in (2.50), and deduce that (1.47) is satisfied.

The formulation (1.47) with the incompressibility constraint (2.39) is supplemented by the boundary conditions

[TABLE]

Moreover, using (2.52) we can show that $\boldsymbol{U}\in C([0,T];L^{2}_{weak}(\Omega))$, and we thus deduce also that

[TABLE]

in other words, $\boldsymbol{U}(0,\cdot)=\boldsymbol{H}[\boldsymbol{U}_{0}]$.

The formulation (1.47) is satisfied for all functions $\boldsymbol{w}$ that are smooth, divergence-free, and $\boldsymbol{w}\cdot{\boldsymbol{n}}|_{\partial\Omega}=0$. However this family of test functions can be extended to the one mentioned in (1.47) by a density argument and the estimates obtained in Section 2.2.

3 Concluding remarks

We studied the behaviour of global-in-time weak solutions to a class of bead-spring chain models with finitely extensible nonlinear elastic (FENE) spring potentials for dilute polymeric fluids, and we proved that as the Mach number tends to zero the system is driven to its incompressible counterpart. Our analysis was performed under the assumption that the velocity of the solvent in which the polymer molecules are immersed satisfies a complete slip boundary condition. The corresponding passage to the limit in the case of a nonslip boundary condition for the velocity field is, as in the case in the compressible Navier–Stokes system (i.e. in the absence of coupling to the Fokker–Planck equation), more technical (cf. Ch.7 of [7], particularly the first paragraph of subsection 7.1.2) and will be considered elsewhere. The existence of global weak solutions to the coupled incompressible Navier–Stokes–Fokker–Planck system was proved in [1] on arbitrary bounded open Lipschitz domains $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$, while the existence proof in [3] for the corresponding compressible Navier–Stokes–Fokker–Planck system was restricted to bounded domains $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$, with $\partial\Omega\subset C^{2,\alpha}$, $\alpha\in(0,1)$. Using the ideas in [8] the $C^{2,\alpha}$ regularity of $\partial\Omega$ assumed in [3] can be relaxed to $\Omega$ being a bounded open Lipschitz domain. For simplicity and for the sake of consistency with the assumptions on $\Omega$ in [7] we have however restricted ourselves here to domains of class $C^{2,\alpha}$.

Acknowledgments

AWK was partially supported by a Newton Fellowship of the Royal Society and by the grant Iuventus Plus 0871/IP3/2016/74 of Ministry of Sciences and Higher Education RP.