Universidad de Cantabria, E.T.S. de Ingenieros Industriales y de Telecomunicacíon

Departmento de Matemática Aplicada y Ciencias de la Computacíon

Avd. Los Castros 44, 39005 Santander, Spain

Žitná 25, CZ-115 67 Praha 1, Czech Republic

\(^1\) BCAM, Basque Center for Applied Mathematics

Mazarredo 14, E48009 Bilbao, Bizkaia, Spain.

\(^2\)IKERBASQUE, Basque Foundation for Science,

Plaza Euskadi 5, 48009 Bilbao, Bizkaia, Spain.

\(^3\)“Simion Stoilow" Institute of the Romanian Academy,

21 Calea Griviţei, 010702 Bucharest, Romania.

2020 Mathematics Subject Classification: 35Q35; 35B40; 74F10; 76D03.

Long time behavior, Fluid-structure interaction, Navier-Stokes equations

1 Introduction↩︎

Recently, the authors in [1], [2] studied the long time behaviour of a single rigid body immersed in a viscous fluid filling the entire space \(R^3\). Specifically, if the initial velocity of the system is small enough, the body approaches some final configuration as time goes to infinity. As pointed out in [2], the methodology could be adapted to the system of several rigid bodies provided the initial velocities are so small that their mutual contacts are excluded. Their approach is based on the concept of strong solution to the fluid structure system, therefore small data are necessary in view of the well known regularity issues related to the 3D incompressible Navier–Stokes system. Let us mention some available works related to the long time behaviour of the fluid-rigid body system in the literature. The authors in [3] analyzed the asymptotic behavior of a rigid ball interacting with an incompressible viscous fluid, establishing precise decay profiles for the coupled system in \(R^3\). In [4], the long-time dynamics of a rigid disk moving in a two-dimensional viscous fluid, proving convergence toward steady states under small initial data has been studied. In [5], the authors investigated the interaction of an anisotropic rigid body with a Poiseuille flow in an unbounded two-dimensional channel, characterizing the asymptotic velocity and orientation of the body over time. In [6], the authors constructed solutions with unbounded energy for the fluid–disk system and analyzed the long-time behavior for large initial data, highlighting nonlinear phenomena absent in small-data regimes. The article [7] is devoted to the analysis of large time behavior of a rigid body of arbitrary shape under the presence of external forces and torques and it shows that if the external forces vanish for large time in \(L^2\) sense then the global strong solution converges to zero as time approaches infinity. The large time behavior of solutions around nontrivial basic states for fluid-rigid ball interaction problem has been analyzed in [8]. In [9], the authors establish that if the initial data for the fluid-rigid body system is small under appropriate norm, then the kinetic energy of the system vanishes as \(t\rightarrow \infty\). Related results on the long-time behaviour of fluid-rigid body interaction systems under a spring-damper feedback force were obtained in [10]–[12].

We propose an alternative approach based on an appropriate concept of weak solutions. For the system of equations that describes the evolution of a family of rigid bodies immersed in a Navier-Stokes fluid contained in a bounded domain \(\Omega \subset R^3\), we introduce a concept of dissipative weak solution based on the original idea of Judakov [13] combined with a suitable form of energy balance. These weak solutions naturally include the case of strong solutions while having some more structure than the distributional solutions, namely having what can be described informally as a differential form of the energy inequality. We show that global–in–time dissipative solutions always exist for any finite energy initial data. Moreover, in the absence of driving forces, the system stabilizes to a static equilibrium for large time. The result holds for

• arbitrarily large initial kinetic energy of the system (no smallness assumption on the velocity);

• independently of possible collisions of several rigid bodies and/or collisions of the bodies with the domain boundary.

The above scenario should be compared with that of the Navier-Stokes fluid without rigid bodies and a possibly large initial velocity. Although the solution may not be regular, the velocity stabilizes exponentially to zero as long as external forces are absent.

1.1 Problem formulation↩︎

We consider a domain \(\Omega \subset R^3\) containing a viscous incompressible fluid and with several rigid bodies moving inside the fluid. The rigid bodies are connected compact sets \((\mathcal{S}^{sec:i})_{i = 1}^N\). The part of the domain \(\Omega\) occupied by the fluid at a time \(t\) is given by \[\Omega_{F}(t) = \Omega \setminus \cup_{i = 1}^N \mathcal{S}^i(t),\;and {we} set\;Q_{F} = \left\{ (t,x) \;\Big|\;t \in (0,T),\;x \in \Omega_{F}(t) \right\}.\] For the sake of simplicity, we set the fluid density \(\varrho_F\) to be \(1\). The motion of the fluid in \(Q_{F}\) is described by the fluid velocity \({\boldsymbol{u}}_F = {\boldsymbol{u}}_{F}(t,x)\) satisfying the Incompressible Navier–Stokes system: \[\begin{align} {\rm div}_x\mathbf{u}_{F} &= 0, \tag{1} \\ \partial_t \mathbf{u}_{F} + {\rm div}_x(\mathbf{u}_{F} \otimes \mathbf{u}_{F}) + \nabla_x\Pi &= {{\rm div}_x\mathbb{S} (\mathbb{D}_x\mathbf{u}_{F})} + {\boldsymbol{g}}, \tag{2} \\ {\mathbb{S}(\mathbb{D}_x\mathbf{u}) = \mu \mathbb{D}_x\mathbf{u}},\; \mathbb{D}_x\mathbf{u}& = \frac{ \nabla_x\mathbf{u}+ \nabla_x^t \mathbf{u}}{2} ,\;\mu > 0 \tag{3}, \end{align}\] where \(\Pi\) is the pressure and the function \({\boldsymbol{g}}\) denotes an external volume force. Denoting \(\varrho_{\mathcal{S}^i}\) the mass density of each body, we define the total mass \(m^{sec:i}\) and the center of mass \({\boldsymbol{h}}_{sec:i}(t)\) as \[m^i = \int_{\mathcal{S}^i(t)} \varrho_{\mathcal{S}^i}(t,x)\, {\rm d} x,\quad {\boldsymbol{h}}_{sec:i}(t)= \frac{1}{m^i} \int_{\mathcal{S}^i(t)} \varrho_{\mathcal{S}^i}(t,x)\, x\, {\rm d} x.\] Moreover, we consider the inertial tensor \(J^{sec:i}\) defined through its action on vectors \(a\) and \(b\) as: \[J^i a \cdot b = \int_{\mathcal{S}^i(t)} \varrho_{\mathcal{S}^i}(t,x)\, [a \times (x - {\boldsymbol{h}}_{sec:i}(t))] \cdot [b \times (x - {\boldsymbol{h}}_{sec:i}(t))]\, {\rm d} x.\] The associated “rigid” velocities at each point \(x\in \mathcal{S}^{sec:i}(t)\) are given by \[\label{m5} \mathbf{u}_{\mathcal{S}^i} (t,x) = {\boldsymbol{Y}}_i(t) + \mathbb{Q}_{sec:i}(t) ( x - {\boldsymbol{h}}_i(t)),\tag{4}\] where \[{\boldsymbol{Y}}_i(t) = \frac{{\rm d}}{\,{\rm d} t } {\boldsymbol{h}}_i (t) ,\;\mathbb{Q}_i(t) = \frac{ {\rm d}}{\,{\rm d} t } \mathbb{O}_i (t) \circ {\mathbb{O}_i}^{-1}(t),\;\mathbb{O}_i \in SO(3).\] Here \({\boldsymbol{Y}}_i\) denotes the translational velocity and \(\mathbb{Q}_{sec:i}\) is the angular velocity of the body. We suppose \({\boldsymbol{Y}}_i \in L^\infty_{\rm loc}([0, \infty); R^3)\), \(\mathbb{Q}_i \in L^\infty_{\rm loc} ([0, \infty); R^{3 \times 3})\), \(i=1,\dots, N\). The matrix \(\mathbb{Q}_{sec:i}\) is skew-symmetric, therefore it can be represented by a vector \(\omega^i\), \[\mathbb{Q}_{sec:i}(t)\big(x - {\boldsymbol{h}}_i(t)\big) = \omega^i \times (x - {\boldsymbol{h}}_i(t)).\] The balance of linear and angular momentum for the body \(\mathcal{S}^i(t)\) reads \[\begin{align} m^i \frac{d}{dt} {\boldsymbol{Y}}_i(t) &= \int_{\partial \mathcal{S}^i(t)} (\mathbb{S} (\mathbb{D}_x\mathbf{u}_{F}) - pI)\, n\, d\sigma + \int_{\mathcal{S}^i(t)} \varrho_{\mathcal{S}^i} {\boldsymbol{g}}\, {\rm d} x, \tag{5}\\[1ex] J^i(t) \frac{d}{dt} \omega^i(t) &= J^i \omega^i(t) \times \omega^i(t) + \int_{\partial \mathcal{S}^i(t)} (x - {\boldsymbol{h}}_i(t)) \times (\mathbb{S} (\mathbb{D}_x\mathbf{u}_{F}) - pI)\, n\, d\sigma \notag \\ &\quad + \int_{\mathcal{S}^i(t)} \varrho_{\mathcal{S}^i} (x - {\boldsymbol{h}}_i(t)) \times {\boldsymbol{g}}\, {\rm d} x, \qquad i=1,\dots,N. \tag{6} \end{align}\] The configuration of the rigid bodies at a time \(t \geq 0\) is parametrized by a family of affine isometries \((\sigma_i(t))_{i=1}^N\), such that \[\label{m4} \mathcal{S}^i (t) = \sigma^i [\mathcal{S}^i],\;\sigma^i(t)x = \mathbb{O}_{sec:i} (t) x + {\boldsymbol{h}}_{sec:i} (t),\;\mathbb{O}_i \in SO(3).\tag{7}\]

The state of the system is described by its mass density \(\varrho= \varrho(t,x)\) and the velocity \({\boldsymbol{u}}= {\boldsymbol{u}}(t,x)\). These quantities are defined throughout \(\Omega\) by prescribing their values in the fluid region and inside each rigid body, namely: \[\begin{align} \varrho(t, x) = \left\{ \begin{array}{l} 1 (= \varrho_F) \;for\; x \in \Omega_F(t), \\ \varrho_{\mathcal{S}^i} \;for\;x \in \mathcal{S}^i(t), \end{array} \right. \tag{8} \\ {\boldsymbol{u}}(t,x) = \left\{ \begin{array}{l} {\boldsymbol{u}}_F(t,x) \;for\; x \in \Omega_F(t), \\ {\boldsymbol{u}}_{\mathcal{S}^i}(t,x) \;for\;x \in \mathcal{S}^i(t) \end{array} \right. \tag{9} \end{align}\] for a.a. \(t > 0\). Moreover we assume no-slip boundary condition:

\[{{\boldsymbol{u}}_F |_{\partial \mathcal{S}_i} = {\boldsymbol{u}}_{\mathcal{S}_i},\;i = 1,\dots, N} ,\;{\boldsymbol{u}}|_{\partial \Omega} = 0.\] For the sake of simplicity, we have supposed the mass densities \(\varrho_{\mathcal{S}^i}\) of the rigid bodies are constant.

1.1.1 Mass conservation↩︎

The density \(\varrho\) and the velocity \({\boldsymbol{u}}\) satisfy the equation of continuity \[\label{m8} \int_0^\infty \int_{\Omega} \Big[ \varrho\partial_t \varphi + \varrho{\boldsymbol{u}}\cdot \nabla_x\varphi \Big] \ \,{\rm d} {x} \,{\rm d} t = - \int_{\Omega} \varrho_0 \varphi (0, \cdot) \ \,{\rm d} {x}\tag{10}\] for any \(\varphi \in C^1_c([0, \infty) \times \overline{\Omega})\). As the velocity \({\boldsymbol{u}}\) is solenoidal, DiPerna-Lions theory [14] asserts that for a given velocity field \[\label{m10} {\boldsymbol{u}}\in L^\infty ([0, \infty); L^2(\Omega; R^3)) \cap L^2(0, \infty; W^{1,2}_0 (\Omega; R^3)),\;{\rm div}_x{\boldsymbol{u}}= 0,\tag{11}\] equation 10 admits a renormalized solution \(\varrho\) unique in the class \[\label{m9} \varrho\in C([0, \infty; L^1(\Omega)) \cap L^\infty ((0, \infty) \times \Omega).\tag{12}\]

1.1.2 Momentum balance↩︎

The weak formulation of the momentum balance goes back to Judakov [13] with a slight modification introduced in [15], [16]. We require the integral identity \[\begin{align} \int_0^\infty &\int_{\Omega} \Big[ \varrho{\boldsymbol{u}}\cdot \partial_t \boldsymbol{\varphi}+ \varrho{\boldsymbol{u}}\otimes {\boldsymbol{u}}: \nabla_x\boldsymbol{\varphi}- \mathbb{S}(\mathbb{D}_x{\boldsymbol{u}}) : \mathbb{D}_x\boldsymbol{\varphi}\Big] \ \,{\rm d} {x} \,{\rm d} t \nonumber \\ &= - \int_0^\infty \int_{\Omega} \varrho{\boldsymbol{g}} \cdot \boldsymbol{\varphi} \ \,{\rm d} {x} \,{\rm d} t - \int_{\Omega} (\varrho_0 {\boldsymbol{u}}_0) \cdot \boldsymbol{\varphi}(0, \cdot) \ \,{\rm d} {x} \label{m11} \end{align}\tag{13}\] to hold for any test function \(\boldsymbol{\varphi}\in C^1_c([0, \infty) \times \Omega; R^3)\), \({\rm div}_x\boldsymbol{\varphi}= 0\), \[\label{m12} \mathbb{D}_x\boldsymbol{\varphi}(t, \cdot) = 0 \;on an open neighbourhood of\; \mathcal{S}^i(t),\;i = 1,\dots, N,\tag{14}\] where the velocity field \({\boldsymbol{u}}\) belong to the class 11 and satisfies the compatibility conditions 9 . In particular, \[\mathbb{D}_x\boldsymbol{\varphi}(t, \cdot) = 0 \;on\;\mathcal{S}^i (t) \;for a.a.\;t \in (0, \infty).\] Note that the weak formulations 10 , 13 , together with the compatibility condition 9 includes both the Navier-Stokes system satisfied in the fluid region as well as the validity of the momentum equations 5 , 6 on the rigid bodies.

1.1.3 Energy balance↩︎

In contrast with [15], [16], we require the energy balance to hold in the differential form \[\begin{align} - \int_0^\infty \partial_t \psi \int_{\Omega} \frac{1}{2} \varrho|{\boldsymbol{u}}|^2 \ \,{\rm d} {x} \,{\rm d} t &+ \int_0^\infty \psi \int_{\Omega} \mathbb{S}(\mathbb{D}_x{\boldsymbol{u}}) : \mathbb{D}_x{\boldsymbol{u}} \ \,{\rm d} {x} \,{\rm d} t \nonumber \\ &\leq \frac{1}{2} \psi (0) \int_{\Omega} \varrho_0 |{\boldsymbol{u}}_0|^2 \ \,{\rm d} {x} + \int_0^\infty \psi \int_{\Omega} \varrho{\boldsymbol{g}} \cdot {\boldsymbol{u}} \ \,{\rm d} {x} \,{\rm d} t \label{m13} \end{align}\tag{15}\] for any \(\psi \in C^1_c[0, \infty)\), \(\psi \geq 0\).

1.2 Dissipative weak solutions↩︎

The original concept of weak solution à la Judakov is based on a more restrictive interpretation of the momentum equation, where the test functions satisfy the compatibility condition \(\mathbb{D}_x\boldsymbol{\varphi}= 0\) only on the rigid bodies. Enlarging the class of test functions to 14 allows to accommodate possible contacts of several rigid bodies and/or rigid bodies with the boundary \(\partial \Omega\). Such a concept of weak solution was used in [15], [16] to develop the theory of global in time weak solutions including possible collisions in both the compressible and the incompressible case. In both [15] and [16], the energy inequality was imposed in the following form \[\begin{align} \int_{\Omega} \frac{1}{2} \varrho|{\boldsymbol{u}}|^2 (\tau, \cdot) \ \,{\rm d} {x} + & \int_0^\tau \int_{\Omega} \mathbb{S}(\mathbb{D}_x{\boldsymbol{u}}) : \mathbb{D}_x{\boldsymbol{u}} \ \,{\rm d} {x} \,{\rm d} t \nonumber \\ &\leq \int_{\Omega} \frac{1}{2} \varrho_0 |{\boldsymbol{u}}_0|^2 \ \,{\rm d} {x} + \int_0^\tau \int_{\Omega} \varrho{\boldsymbol{g}} \cdot {\boldsymbol{u}} \ \,{\rm d} {x} \,{\rm d} t \label{m14} \end{align}\tag{16}\] for a.a. \(\tau \geq 0\).

Our goal is to show that the construction of weak solutions presented in [16] can be modified to accommodate the differential version of the energy inequality 15 . This has an immediate impact on the long–time behaviour of global in time solutions. Whether or not the same is possible for the rigid bodies immersed in a viscous compressible fluid remains an open problem.

The paper is organized as follows. In Section 2, we revisit the construction of the weak solution in [16] and show how to accommodate the differential form of the energy inequality 15 . Applications to the long–time behaviour of dissipative solutions are discussed in Section 3.

2 Weak solutions to the fluid structure interaction problem revisited↩︎

Our goal is to show the existence of global–in–time dissipative solutions for any finite energy initial data \(\varrho_0\), \({\boldsymbol{u}}_0\) and any distribution of the rigid objects \[\label{rr1} \mathcal{S}^i,\;\mathcal{S}^i = {\rm cl}[{\rm int}[\mathcal{S}^i]],\; \mathcal{S}^i \;compact, connected,\;i = 1,\dots, N\tag{17}\] such that their initial positions satisfy \[\label{rr2} {\rm int}[\mathcal{S}^i(0)] \cap {\rm int}[\mathcal{S}^j(0)] = \emptyset \;as soon as\; i \ne j.\tag{18}\] To this end, we follow step by step the construction presented in [16] and show how to accommodate the strong form of the energy inequality 15 .

The strategy of [16] consists of several steps:

1. For smooth and initially well separated bodies and a smooth spatial domain \(\Omega\) construct a weak solution up to the time of the first contact, meaning either the contact of two or more rigid bodies or a contact of a rigid body with \(\partial \Omega\), see Proposition 2.3 and Corollary 2.4 in [16].

2. Extend the result to arbitrary compact bodies as in 17 and \(\Omega\) a general bounded open set, see [16]. This is done by means of an approximation procedure using the result of Step 1. The solution exist up to the time of the first contact.

3. Continue the solution after the time of the first contact considering the colliding rigid objects as a single body and repeating the arguments of Step 2. As there is a finite number of possible collision times, the process gives rise to a global in time weak solution. This has been explored in [16].

It is worth noting that at each step of the above procedure, one has to perform the limit in the convective term \[\varrho_n{\boldsymbol{u}}_n\otimes {\boldsymbol{u}}_n\to \varrho{\boldsymbol{u}}\otimes {\boldsymbol{u}}\] in a weak sense. In particular, one obtains the convergence of the weak time derivative of the kinetic energy \[\int_0^\infty \partial_t \psi \int_{\Omega} \frac{1}{2} \varrho_n |{\boldsymbol{u}}_n|^2 \ \,{\rm d} {x}\,{\rm d} t \to \int_0^\infty \partial_t \psi \int_{\Omega} \frac{1}{2} \varrho|{\boldsymbol{u}}|^2 \ \,{\rm d} {x}\,{\rm d} t .\] Thus the differential form of the energy inequality 15 remains valid for the limit solution as long as it holds for the approximate solutions constructed in Step 1. Our goal therefore reduces to finding a suitable construction of solutions with smooth data (up to the first contact) that would satisfy the energy balance 15 .

2.1 Approximate solutions with smooth data↩︎

In accordance with the previous discussion, we consider the situation of Step 1 of initially well separated smooth bodies and a smooth spatial domain \(\Omega\). Our goal is construct a weak solution up to the first contact satisfying the differential form of the energy inequality 15 .

We use the idea of San Martin, Starovoitov, and Tucsnak [17] and approximate the rigid bodies by a fluid of high viscosity. Specifically, adapting the strategy proposed in [18] we proceed as follows:

• As \(\Omega \subset R^3\) is bounded, we may suppose \[\Omega \subset (-L,L)^3\] for \(L > 0\) large enough. Thus we may consider the problem with the spatially periodic conditions, meaning on the flat torus \[\mathbb{T}^3 = \left( [-L,L]|_{\{ -L; L \} } \right)^3.\]

• We consider a system of equations

  \[\begin{align} \partial_t \varrho+ {\rm div}_x(\varrho[{\boldsymbol{u}}]_\delta ) &= 0, \tag{19} \\ \partial_t (\varrho{\boldsymbol{u}}) + {\rm div}_x(\varrho{\boldsymbol{u}}\otimes [{\boldsymbol{u}}]_\delta) + \nabla_xP &= {\rm div}_x([\mu_{\varepsilon}]_\delta \mathbb{D}_x{\boldsymbol{u}}) - \chi_\varepsilon{\boldsymbol{u}}+ \varrho{\boldsymbol{g}}, \tag{20} \\ \partial_t \mu_{\varepsilon} + {\rm div}_x(\mu_{\varepsilon} [{\boldsymbol{u}}]_\delta) &= 0,\;\tag{21}\\ {\rm div}_x{\boldsymbol{u}}&= 0, \tag{22} \end{align}\] where \[[h]_\delta = \sigma_\delta * h\] denotes a convolution in space with a family of regularization kernels with \[\sigma_{\delta}(x)=\frac{1}{\delta^3}\sigma\left(\frac{|x|}{\delta}\right),\;\sigma\in C_c^{\infty}(-1,1),\;\sigma(z)>0for -1<z<1, \sigma(z)=\sigma(-z),\;\int\limits_{-1}^{1} \sigma(z)\;dz=1.\]

  The quantity \(\chi_\varepsilon\) represents a penalization term, \[\chi_\varepsilon= \frac{1}{\varepsilon} \chi,\;\chi \in C^\infty(\mathbb{T}^3),\; \chi(x) = \left\{ \begin{array}{l} 0 \;for\;x \in \overline{\Omega}, \\ > 0 \;for\;x \in \mathbb{T}^3 \setminus \overline{\Omega}. \end{array} \right.\]

  The quantity \(\mu_\varepsilon\) provides a regularized characterization of the rigid body regions, whose evolution is governed by the transport equation 21 ; similarly the friction coefficient \(\chi_\varepsilon\) is a penalization of the no-slip boundary condition on \(\partial \Omega\), see [18] for details.

• System 19 –22 is supplemented with smooth initial data \[\begin{align} \varrho(0, \cdot) = \varrho_{0, \delta} = \left\{ \begin{array}{l} 1 (= \varrho_F) \;for\;x \in \mathbb{T}^3 \setminus \cup_{i=1}^N \mathcal{S}^i (0) ,\\ \in [0, \varrho_{S^i}] \;for\;x \in \mathcal{S}^i(0),\; {\rm dist}[x, \partial \mathcal{S}^i(0)] < \delta, \\ \varrho_{S^i} \;for\;x \in \mathcal{S}^i(0),\; {\rm dist}[x, \partial \mathcal{S}^i] \geq \delta,\;i = 1,\dots, N; \end{array}\right. \label{r5} \end{align}\tag{23}\] \[\begin{align} {\boldsymbol{u}}(0, \cdot) = {\boldsymbol{u}}_{0, \delta} = \left\{ \begin{array}{l} 0 \;for\;x \in \mathbb{T}^3 \setminus \Omega ,\\ {\boldsymbol{u}}_{0,F} \;for\;x \in \Omega \setminus \cup_{i=1}^N \mathcal{S}^i(0),\\ {\boldsymbol{u}}_{0,S^i} \;for\;x \in \mathcal{S}^i(0); \end{array}\right. \label{r6} \end{align}\tag{24}\] \[\begin{align} \mu(0, \cdot) = \mu_{0, \varepsilon, \delta} = \left\{ \begin{array}{l} 1 \;for\;x \in \mathbb{T}^3 \setminus \cup_{i = 1}^N \mathcal{S}^i ,\\ \in \left[ 0,\;\frac{1}{\varepsilon} \right] \;for\; x \in \mathcal{S}^i(0),\;{\rm dist}[x, \partial \mathcal{S}^i] < \delta, \\ \frac{1}{\varepsilon} \;for\;x \in \mathcal{S}^i(0),\;{\rm dist}[x, \partial \mathcal{S}^i(0)] \geq \delta,\;i = 1,\dots, N. \end{array}\right. \label{r7} \end{align}\tag{25}\]

As shown in [18], system 19 –22 , supplemented with the initial conditions 23 – 25 , admits a solution defined in \((0,T) \times \mathbb{T}^3\), see [18]. Moreover, performing successively the limits \(\varepsilon\to 0\), \(\delta \to 0\), the solutions converge to the desired solution of the fluid structure interaction problem as long as there are no contacts of rigid objects, meaning up to the first contact.

Our ultimate goal is therefore to show that the limit solution satisfies the energy balance in the differential form. To see this, we first record the energy balance for the original system: \[\frac{{\rm d}}{\,{\rm d} t } \int_{\mathbb{T}^3} \frac{1}{2} \varrho|{\boldsymbol{u}}|^2 \,{\rm d} {x} + \int_{\mathbb{T}^3} [\mu]_\delta |\mathbb{D}_x{\boldsymbol{u}}|^2 \,{\rm d} {x}+ \int_{\mathbb{T}^3} \chi_\varepsilon |{\boldsymbol{u}}|^2 \,{\rm d} {x}\leq \int_{\mathbb{T}^3} \varrho{\boldsymbol{g}} \cdot {\boldsymbol{u}}\,{\rm d} {x}.\] As a matter of fact, due to the regularization of the convective terms, the approximate solutions are regular enough for the energy balance to be derived directly form the equations. Using “distributional” formulation, we may also include the initial energy, \[\begin{align} - \int_0^T \partial_t \psi \int_{\mathbb{T}^3} \frac{1}{2} \varrho|{\boldsymbol{u}}|^2 \,{\rm d} {x}\,{\rm d} t &+ \int_0^T \psi \int_{\mathbb{T}^3} [\mu]_\delta |\mathbb{D}_x{\boldsymbol{u}}|^2 \,{\rm d} {x}\,{\rm d} t + \int_0^T \psi \int_{\mathbb{T}^3} \chi_\varepsilon |{\boldsymbol{u}}|^2 \,{\rm d} {x}\,{\rm d} t \nonumber \\ &\leq \psi(0) \int_{\mathbb{T}^3} \frac{1}{2} \varrho_{0,\delta} |{\boldsymbol{u}}_{0, \delta}|^2 \,{\rm d} {x} + \int_0^T \psi \int_{\mathbb{T}^3} \varrho{\boldsymbol{g}} \cdot {\boldsymbol{u}}\,{\rm d} {x}\,{\rm d} t \label{r8} \end{align}\tag{26}\] for any \(\psi \in C^1_c[0,T)\), \(\psi \geq 0\).

Thanks to the specific choice of initial data, we have \[\label{r9} \int_{\mathbb{T}^3} \frac{1}{2} \varrho_{0,\delta} |{\boldsymbol{u}}_{0, \delta}|^2 \,{\rm d} {x}= \int_{\Omega} \frac{1}{2} \varrho_{0,\delta} |{\boldsymbol{u}}_{0, \delta}|^2 \ \,{\rm d} {x}.\tag{27}\]

Next, because of the domain penalization term, \[\label{r10} \int_0^T \int_{\mathbb{T}^3 \setminus \Omega} |{\boldsymbol{u}}|^2 \stackrel{<}{\sim}\varepsilon\; \Rightarrow \;\int_0^T \int_{\mathbb{T}^3 \setminus \Omega} \varrho|{\boldsymbol{u}}|^2 \stackrel{<}{\sim}\varepsilon\tag{28}\] since the densities remain bounded uniformly in time by their initial values. Moreover, by the same token, the artificial viscosities \([\mu]_\delta\) remain bounded below, \[\label{r11} [\mu]_{\delta} \geq 1.\tag{29}\]

In view of the above observations, we may rewrite 26 as \[\begin{align} - \int_0^T \partial_t \psi \int_{\Omega} \frac{1}{2} \varrho|{\boldsymbol{u}}|^2 \ \,{\rm d} {x} \,{\rm d} t &+ \int_0^T \psi \int_{\Omega} |\mathbb{D}_x{\boldsymbol{u}}|^2 \ \,{\rm d} {x} \,{\rm d} t \nonumber \\ &\leq \psi(0) \int_{\Omega} \frac{1}{2} \varrho_{0,\delta} |{\boldsymbol{u}}_{0, \delta}|^2 \ \,{\rm d} {x} + \int_0^T \psi \int_{\Omega} \varrho{\boldsymbol{g}} \cdot {\boldsymbol{u}} \ \,{\rm d} {x} \,{\rm d} t + \omega\left(\varepsilon, \| \psi \|_{C^1[0,T]} \right) \label{r12} \end{align}\tag{30}\] where the error term satisfies \[\omega(\varepsilon, \| \psi \|_{C^1[0,T]}) \to 0 \;as\;\varepsilon\to 0,\] for any (fixed) \(\psi \in C^1_c[0,T)\), \(\psi \geq 0\).

At this stage, we can perform the double limit \(\varepsilon\to 0\), \(\delta \to 0\) in 30 . As shown in [18], the convective term \[\label{r13} \varrho_{\varepsilon, \delta} {\boldsymbol{u}}_{\varepsilon, \delta} \otimes {\boldsymbol{u}}_{\varepsilon, \delta} \stackrel{\varepsilon\to 0}{\to} \varrho_\delta {\boldsymbol{u}}_{\delta} \otimes {\boldsymbol{u}}_{\delta} \stackrel{\delta \to 0}{\to} \varrho{\boldsymbol{u}}\otimes {\boldsymbol{u}}\;weakly in\;L^q((0,T) \times \Omega; R^{3 \times 3})\;for some\;q > 1.\tag{31}\] In particular, we conclude by the convergence of the weak time derivative of the kinetic energy, \[\int_0^T \partial_t \psi \int_{\Omega} \frac{1}{2} \varrho_{\varepsilon, \delta} |{\boldsymbol{u}}_{\varepsilon, \delta}|^2 \ \,{\rm d} {x} \,{\rm d} t \stackrel{\varepsilon\to 0}{\to} \int_0^T \partial_t \psi \int_{\Omega} \frac{1}{2} \varrho_\delta |{\boldsymbol{u}}_\delta|^2 \ \,{\rm d} {x} \stackrel{\delta \to 0}{\to} \int_0^T \partial_t \psi \int_{\Omega} \frac{1}{2} \varrho|{\boldsymbol{u}}|^2 \ \,{\rm d} {x}.\]

We have shown the following result.

3 Long–time behaviour of dissipative solutions↩︎

The fact the dissipative solutions satisfy the energy balance in the differential form 15 has an immediate impact on their long–time behaviour.

3.1 System without external forcing↩︎

Consider the case \({\boldsymbol{g}} = 0\). First, by Korn–Poincaré inequality, the energy balance 15 gives rise to \[\begin{align} - \int_0^\infty \partial_t \psi \int_{\Omega} \frac{1}{2} \varrho|{\boldsymbol{u}}|^2 \ \,{\rm d} {x} \,{\rm d} t + C_{K,P} \int_0^\infty \psi \int_{\Omega} |{\boldsymbol{u}}|^2 \ \,{\rm d} {x} \,{\rm d} t \leq \frac{1}{2} \psi (0) \int_{\Omega} \varrho_0 |{\boldsymbol{u}}_0|^2 \ \,{\rm d} {x} \nonumber \end{align}\] for any \(\psi \in C^1_c[0, \infty)\), \(\psi \geq 0\) and a certain positive constant \(C_{K,P}\). Moreover, as the density is uniformly bounded, we get \[\label{L1} - \int_0^\infty \partial_t \psi \int_{\Omega} \frac{1}{2} \varrho|{\boldsymbol{u}}|^2 \ \,{\rm d} {x} \,{\rm d} t + \frac{2 C_{K,P}}{\overline{\varrho}} \int_0^\infty \psi \int_{\Omega} \frac{1}{2} \varrho|{\boldsymbol{u}}|^2 \ \,{\rm d} {x} \,{\rm d} t \leq \frac{1}{2} \psi (0) \int_{\Omega} \varrho_0 |{\boldsymbol{u}}_0|^2 \ \,{\rm d} {x}.\tag{32}\] for any \(\psi \in C^1_c[0, \infty)\), \(\psi \geq 0\). Relation 32 yields the exponential decay of the kinetic energy, \[\label{L2} \int_{\Omega} \varrho|{\boldsymbol{u}}|^2 (\tau, \cdot) \ \,{\rm d} {x} \leq \exp\left( - \frac{2 C_{K,P}}{\overline{\varrho}} \tau \right) \int_{\Omega} \varrho_0 |{\boldsymbol{u}}_0|^2 \ \,{\rm d} {x},\;\tau \geq 0.\tag{33}\]

Finally, we use the fact the mass densities of the rigid objects are bounded below away from zero and their interiors are non–empty. This yields the exponential decay for the corresponding rigid velocities, \[\label{L3} | {\boldsymbol{Y}}_i(\tau) | = \left| \frac{{\rm d}}{\,{\rm d} t } {\boldsymbol{h}}_i (\tau) \right| \leq C \exp (- \Lambda \tau), \;|\mathbb{Q}_i(\tau)| \leq C \exp (- \Lambda \tau),\;i = 1, \dots, N,\tag{34}\] for some \(\Lambda > 0\).

We have shown the following result.

3.2 Gravitational forcing↩︎

We conclude the paper by a short discussion of the asymptotic behaviour of the system driven by gravitational (potential) external force. Accordingly, we suppose \[{\boldsymbol{g}} = \nabla_xG,\;G = G(x) \;a smooth potential.\] In this case, the forcing term on the right–hand side of 15 can be written with the help of the equation of continuity as \[\int_0^\infty \psi \int_{\Omega} \varrho\nabla_xG \cdot {\boldsymbol{u}} \ \,{\rm d} {x} = - \int_0^\infty \partial_t \psi \int_{\Omega} \varrho G \ \,{\rm d} {x} + \psi(0) \int_{\Omega} \varrho_0 G \ \,{\rm d} {x}.\] Consequently, the energy inequality 15 takes the form \[\begin{align} - \int_0^\infty \partial_t \psi \int_{\Omega} \varrho\left( \frac{1}{2} |{\boldsymbol{u}}|^2 - G \right) \ \,{\rm d} {x} \,{\rm d} t &+ \int_0^\infty \psi \int_{\Omega} \mathbb{S}(\mathbb{D}_x{\boldsymbol{u}}) : \mathbb{D}_x{\boldsymbol{u}} \ \,{\rm d} {x} \,{\rm d} t \nonumber \\ &\leq \psi (0) \int_{\Omega} \varrho_0 \left( \frac{1}{2} |{\boldsymbol{u}}_0|^2 - G \right) \ \,{\rm d} {x} \label{L4} \end{align}\tag{35}\] for any \(\psi \in C^1_c[0, \infty)\), \(\psi \geq 0\).

The function \[\tau \mapsto \int_{\Omega} \varrho\left( \frac{1}{2} |{\boldsymbol{u}}|^2 - G \right)(\tau, \cdot) \ \,{\rm d} {x}\] is non–increasing and bounded from below. We may infer \[\label{L5} \int_{\Omega} \varrho\left( \frac{1}{2} |{\boldsymbol{u}}|^2 - G \right)(\tau, \cdot) \ \,{\rm d} {x} \to - \mathcal{E}_\infty \;as\;\tau \to \infty.\tag{36}\]

In addition, we have \[\label{L6} \int_0^\infty \int_{\Omega} |\mathbb{D}_x{\boldsymbol{u}}|^2 \ \,{\rm d} {x} \,{\rm d} t < \infty.\tag{37}\] Moreover, by virtue of the equation of continuity , we deduce \[\label{L7} \int_{\Omega} \varrho{G} (\tau_1, \cdot) \ \,{\rm d} {x} - \int_{\Omega} \varrho{G} (\tau_2, \cdot) \ \,{\rm d} {x} = \int_{\tau_1}^{\tau_2} \int_{\Omega} \varrho{\boldsymbol{u}}\cdot \nabla_xG \ \,{\rm d} {x} \,{\rm d} t ,\tag{38}\] where, furthermore, \[\begin{align} \left| \int_{\tau_1}^{\tau_2} \int_{\Omega} \varrho{\boldsymbol{u}}\cdot \nabla_xG \ \,{\rm d} {x} \,{\rm d} t \right| &\leq \int_{\tau_1}^{\tau_2} \int_{\Omega} \varrho|\nabla_xG|^2 \ \,{\rm d} {x} \,{\rm d} t + \int_{\tau_1}^{\tau_2} \int_{\Omega} \varrho|{\boldsymbol{u}}|^2 \ \,{\rm d} {x} \,{\rm d} t \nonumber \\ &\leq \int_{\tau_1}^{\tau_2} \int_{\Omega} \varrho|\nabla_xG|^2 \ \,{\rm d} {x} \,{\rm d} t + C \int_{\tau_1}^{\tau_2} \int_{\Omega} |\mathbb{D}_x{\boldsymbol{u}}|^2 \ \,{\rm d} {x} \,{\rm d} t . \label{L8} \end{align}\tag{39}\]

Combining 36 –39 we conclude \[\label{L9} \int_{\Omega} \varrho|{\boldsymbol{u}}|^2 (\tau, \cdot) \ \,{\rm d} {x} \to 0,\; \int_{\Omega} \varrho(\tau, \cdot) G \ \,{\rm d} {x} \to \mathcal{E}_\infty \;as\;\tau \to \infty.\tag{40}\] We have shown the following result.

In the full generality allowed by Theorem 6, it seems difficult to identify the final position of rigid bodies. More specific results in the simplified 2D case have been obtained in [19].

References↩︎

---

1. The work of M.B. was supported by the project "Ecuaciones en derivadas parciales motivadas por procesos de difusión y mecánica de fluidos: propiedades, asintóticas y homogeneización" (Ayuda financiada contrato Programa Gob. Cantabria -UC). The work of E.F. was partially supported by the Czech Sciences Foundation (GAČR), Grant Agreement 24–11034S. The Institute of Mathematics of the Academy of Sciences of the Czech Republic is supported by RVO:67985840. E.F. is a member of the Nečas Center for Mathematical Modelling. A.R is supported by the Grant RYC2022-036183-I funded by MICIU/AEI/10.13039/501100011033 and by ESF+. A.R, A.Z have been partially supported by the Basque Government through the BERC 2022-2025 program and by the Spanish State Research Agency through BCAM Severo Ochoa CEX2021-001142-S and through project PID2023-146764NB-I00 funded by MICIU/AEI/10.13039/501100011033 and cofunded by the European Union. A.Z. was also partially supported by a grant of the Ministry of Research, Innovation and Digitization, CNCS - UEFISCDI, project number PN-III-P4-PCE-2021-0921, within PNCDI III.↩︎