Analysis of a model for the dynamics of microswimmer suspensions [original]

Received: 26 November 2020 Revised: 30 April 2021 Accepted: 24 June 2021

DOI: 10.1002/mma.7673

RESEARCH ARTICLE

Analysis of a model for the dynamics of microswimmer

suspensions

Etienne Emmrich Lukas Geuter

Institut für Mathematik, Technische

Universität Berlin, Berlin, Germany

Correspondence

Etienne Emmrich, Institut für

Mathematik, Technische Universität

Berlin, Straße des 17. Juni 136, 10623

Berlin, Germany.

Email: [email protected]

Communicated by: T. Roubíˇ

cek

Funding information

Deutsche Forschungsgemeinschaft,

Grant/Award Number: SFB 910

In this paper, a model that was recently derived in Reinken et al. to describe

the dynamics of microswimmer suspensions is studied. In particular, the global

existence of weak solutions, their weak–strong uniqueness, and a connection to

a different model that was proposed in Wensink et al. is shown.

KEYWORDS

active fluid, existence, relative energy, weak solution, weak–strong uniqueness

MSC CLASSIFICATION

35Q35; 76A05

1INTRODUCTION

Microswimmers are small, self-propelling particles, for example, bacteria like Bacillus subtilis, algae like Chlamydomonas

reinhardtii, or artificial nanorods. Suspended in a liquid, the interaction between these particles themselves and between

them and the surrounding fluid gives rise to interesting phenomena such as active turbulence. For a high density of

microswimmers, such a suspension can be regarded as an example of an active fluid.

There are several approaches to capture the dynamics of an active fluid starting with adaptions of a model proposed in

Vicsek et al.1Later approaches include the derivation of continuum limit hydrodynamic equations as in Toner and Tu2or

adaptions of liquid crystal models as in Thampi and Yeomans.3Another model we want to mention is the one suggested

in Wensink et al,4where a Toner–Tu-like equation is supplemented with a fourth-order Swift–Hohenberg term.

The model that we want to study here was recently derived in Reinken et al5with the goal of it being able to incorporate

short-range interactions favoring alignment as well as long-range hydrodynamic interactions. The authors start from

overdamped Langevin equations for a generic microscopic model and couple them with a Stokes equation supplemented

by an ansatz for the stress tensor.

Our paper is arranged as follows. First, we introduce the equations that are the object of our analysis. Then—under a

simplifying assumption—we prove the existence of weak solutions via a Galerkin approximation. The main part of the

paper is then devoted to prove that such weak solutions obey a relative energy inequality. Such a relative energy inequality

can be employed to a variety of uses, such as showing the stability of equilibria (e.g., in Feireisl6), deriving a posteriori

estimates for modeling errors (see Fischer7), or as the basis of a generalized concept of solution (e.g., in Lasarzik8).

In our case, we will apply the relative energy inequality to prove the weak–strong uniqueness of weak solutions as well

as the convergence of those to strong solutions of another problem as one parameter tends to zero.

Throughout this paper, we denote by c>0 a generic constant that does not depend on any changing quantity.

Furthermore, let Ω⊂R3be a bounded domain of class 𝒞2and T>0 some fixed time.

This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium,

provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

wileyonlinelibrary.com/journal/mmaMath Meth Appl Sci. 2021;44:14041–14058.

14041

EMMRICH AND GEUTER

Spaces of vector- or matrix-valued functions are denoted by bold letters, for example, Lp(Ω) ∶= Lp(Ω; R3)for the

spaces of integrable functions, and Wk,p(Ω) ∶= Wk,p(Ω; R3)as well as Hk(Ω) ∶= Wk,2(Ω; R3)for Sobolev spaces. The

Sobolev space associated with homogeneous Dirichlet boundary conditions is denoted by H1

0(Ω) ∶= closH1𝒞∞

c(Ω; R3),

where 𝒞∞

c(Ω; R3)denotes the set of infinitely many times differentiable functions that are compactly supported in Ω.We

write 𝒞∞

c,𝜎(Ω; R3)for all such functions that are in addition solenoidal and then denote by Lp

𝜎(Ω) and H1

0,𝜎(Ω) the closure

of 𝒞∞

c,𝜎(Ω; R3)with respect to the standard norms of Lp(Ω) and H1(Ω), respectively.

The spaces for time-dependent continuous and continuously differentiable functions mapping [0, T] into a Banach

space Xare denoted by 𝒞([0,T]; X)and 𝒞1([0,T]; X), respectively, while the corresponding spaces of Bochner-integrable

and weakly differentiable functions are denoted by Lp(0, T;X)andW1, p(0, T;X), respectively. We often omit the time

interval (0, T)andthedomainΩin this notation and write, for example, Lp(W1,q). Also, for the sake of brevity, we generally

do not write out the time dependence of functions under the integral.

For an arbitrary Banach space X, we denote its dual by X*and the dual pairing between X*and Xby ⟨·,·⟩. With a slight

abuse of notation, however, the dual pairing between the spaces Lp(Ω) and Lq(Ω),whereqis the conjugate exponent to

p, is denoted by (·,·). We use the same notation for the inner product on the Hilbert space L2(Ω) × L2(Ω).

The space H2(Ω) ∩ H1

0,𝜎(Ω) is equipped with the norm ||Δ·||L2(cf. Boyer & Fabrie9, Proposition IV.5.9).

2MODEL

We want to study the following equations derived in Reinken et al.5Consider the initial-boundary value problem

−Δu+𝜇1Δ2p−𝛾1Δp+𝜆1(p·∇)p+∇𝜋1=0,

𝜕tp+𝜇2Δ2p−𝛾2Δp+𝜆2(p·∇)p+𝛼|p|2p+𝛽p

+(u·∇)p+𝜅(∇u)symp−(∇u)skwp+∇𝜋2=0,

∇·u=∇·p=0,

(1a)

on Ω×(0, T)and

u=p=Δp=0on𝜕Ω×(0,T),

p(·,0)=p0in Ω,(1b)

where 𝛼,𝜆1,𝜆2,𝜇1,𝜇2>0, and 𝛽,𝛾1,𝛾

2,𝜅 ∈R. In contrast to other models, both the influence of the polar ordering of the

microswimmers themselves, described by the vector field p∶

Ω×[0,T]→R3, as well as the one from the velocity field

u∶

Ω×[0,T]→R3of the suspension fluid is taken into account. The latter is incorporated into the model by coupling

the Stokes equation with the equation for the polar ordering parameter p. The strength of the coupling is described by a

dimensionless parameter 𝜀, meaning there exist 𝜀, 

𝜆1,𝜇1>0, and 𝛾1∈Rsuch that

𝛾1=𝜀𝛾1,𝜆

1=𝜀

𝜆1,𝜇

1=𝜀𝜇1.(2)

Then for 𝜀=0, both equations decouple (see Section 5).

Since both the solvent as well as the whole fluid are assumed to be incompressible, it follows that uand pare supposed

to be solenoidal vector fields. The Lagrange multipliers for these divergence-free constraints are denoted by 𝜋1and 𝜋2,

respectively. The velocity field of the whole fluid is then given as the sum u+𝜈pof both vector fields, where 𝜈>0 denotes

a reference velocity of the swimmers. The unknown so-called pressure 𝜋∶

Ω×[0,T]→Renters the equations via the

relation 𝜋=𝜋1+𝜈𝜋2and is the Lagrange multiplier for the divergence-free constraint ∇·(u+𝜈p)=0 on the velocity field

of the whole fluid. For a certain constant c1>0, one can think of c1|p|2as the so-called active pressure and of 𝜋1−c1|p|2

as the effective pressure of the active fluid (cf. Dunkel et al.10). We consider the problem in a variational framework with

solenoidal test functions and will not be dealing with the terms 𝜋,𝜋1,and𝜋2since they vanish in the weak formulation.

Unfortunately, we are not able to show the appropriate a priori estimates to prove the existence of weak solutions for the

general case above. The problematic term here is 𝜅(∇u)symp. If we proceed in a standard fashion to derive a priori estimates

and multiply the second equation by p, integrate and perform an integration by parts, we end up with the expression

−𝜅∫T

0∫Ω

(p·∇)p·udxdt.

14042

EMMRICH AND GEUTER

Since by the first equation uisofthesameorderasΔp, we seem not to be able to estimate and absorb this term into the

left-hand side. Therefore, we restrict ourselves to the simpler case 𝜅=0.

3EXISTENCE OF WEAK SOLUTIONS

We start by giving the definition of a weak solution to (1). Note that the term Δ2poccurs in the Stokes equation, meaning

that Δuhas at most the regularity of Δ2p. Hence, the weak formulation of the second equation necessitates a very weak

formulation for the Stokes equation. For this, by A−1∶L2

𝜎(Ω) →H2(Ω) ∩ H1

0,𝜎(Ω), we denote the inverse of the Stokes

operator, meaning that for 𝜑∈L2

𝜎(Ω), the vector field A−1𝜑∶= vis the unique solution of

−Δv+∇𝜋=𝜑in Ω,

∇·v=0inΩ,

v=0on𝜕Ω

in the weak sense. For the existence and properties of this operator, see Temam.11, Chapter I, § 2.6

Definition 1 (Weak solution). Let p0∈H2(Ω)∩H1

0,𝜎(Ω) be given. A pair (u,p)∈L2(L2

𝜎)×(L∞(L2

𝜎)∩L4(L4

𝜎)∩L2(H2∩

0,𝜎)) such that 𝜕tp∈L4

3((H2∩H1

0,𝜎)∗)is called weak solution to (1) if the equations

∫T

(u,𝜑)dt−𝜇1∫T

(Δp,ΔA−1𝜑)dt+𝛾1∫T

(Δp,A−1𝜑)dt

−𝜆1∫T

((p·∇)p,A−1𝜑)dt=0

(3)

and

∫T

0⟨𝜕tp,𝜓⟩dt+𝜇2∫T

(Δp,Δ𝜓)dt+𝛾2∫T

(∇p,∇𝜓)dt+𝜆2∫T

((p·∇)p,𝜓)dt

+𝛼∫T

(|p|2p,𝜓)dt+𝛽∫T

(p,𝜓)dt+∫T

((u·∇)p,𝜓)dt

2∫T

((p·∇)𝜓−(𝜓·∇)p,u)dt=0

(4)

hold for all test functions 𝜑, 𝜓 ∈𝒞∞

c,𝜎(Ω × (0,T); R3)and the initial condition p(0)=p0is fulfilled in the weak sense,

that is,

p(t)⇀p0in L2

𝜎(Ω) as t→0.

Remark 1. In order to justify that this weak formulation is well defined, we only consider the nonlinear terms. Since

pis divergence-free, the identity

(p·∇)p=∇·(p⊗p),

where ⊗denotes the dyadic product, holds almost everywhere in Ω×(0, T). Performing an integration-by-parts and

using Hölder's inequality as well as the continuous embedding H2(Ω) → W1,6(Ω) then yields

|((p·∇)p,w)|≤c||p||2

L12

5||w||H2

for any w∈H2(Ω) ∩ H1

0,𝜎(Ω). Therefore, we can estimate the norm of the functional w→ ((p·∇)p,w)by

||(p·∇)p||(H2∩H1

0,𝜎 )∗≤c||p||2

L12

14043

EMMRICH AND GEUTER

almost everywhere in (0, T). Hence,

∫T

0||((p·∇)p,A−1𝜑)||dt≤∫T

0||(p·∇)p||(H2∩H1

0,𝜎 )∗||A−1𝜑||H2∩H1

0,𝜎 dt

≤c||p||2

L4(L12

5)||𝜑||L2(L2),

where we also used the boundedness of A−1∶L2

𝜎(Ω) →H2(Ω) ∩ H1

0,𝜎(Ω) (see Temam12, Section 2.5). For the remaining

nonlinear terms, applications of Hölder's inequality yield

∫T

0|((p·∇)p,𝜓)|dt≤||p||L2(L∞)||∇p||L2(L6)||𝜓||L∞(L6

5),

∫T

0|||(|p|2p,𝜓)|||dt≤||p||3

L4(L4)||𝜓||L4(L4),

∫T

0|((u·∇)p,𝜓)|dt≤||u||L2(L2)||∇p||L2(L6)||𝜓||L∞(L3),

∫T

0|((p·∇)𝜓−(𝜓·∇)p,u)|dt≤||p||L2(L∞)||∇𝜓||L∞(L2)||u||L2(L2)

+||𝜓||L∞(L3)||∇p||L2(L6)||u||L2(L2).

(5)

Together with the continuous embedding H2(Ω) → L∞(Ω), we see that this notion of a weak solution is well defined

for all test functions 𝜓∈L4(H2∩H1

0,𝜎)∩L∞(H1

0,𝜎)and 𝜑∈L2(L2

𝜎).

We can prove the existence of such a solution via a Galerkin approximation.

Theorem 1. Let p0∈H2(Ω) ∩ H1

0,𝜎(Ω). Then there exists a weak solution to (1) in the sense of Definition 1.

Proof. First, we choose a Galerkin basis {wm}m∈N⊂H2(Ω) ∩ H1

0,𝜎(Ω) consisting of eigenfunctions of the Stokes

operator A∶H2(Ω) ∩ H1

0,𝜎(Ω) ⊂L2

𝜎(Ω) →L2

𝜎(Ω).Foranyn∈N, let Wnbe the finite dimensional subspace of

H2(Ω) ∩ H1

0,𝜎(Ω) spanned by w1,…,wnand Πn∶L2

𝜎(Ω) →Wnthe L2(Ω)-orthogonal projection onto Wn.Forthe

initial value p0,n∶= Πnp0, we are then looking for a solution (un,pn)∈ 𝒞1([0,T]; Wn)×𝒞1([0,T]; Wn)to the discretised

problem

(un,v)−𝜇1(Δpn,ΔA−1v)+𝛾1(Δpn,A−1v)−𝜆1(pn·∇)pn,A−1v)=0,(6a)

(𝜕tpn,w)+𝜇2(Δpn,Δw)+𝛾2(∇pn,∇w)+𝜆2((pn·∇)pn,w)+𝛼(|pn|2pn,w)

+𝛽(pn,w)+((un·∇)pn,w)+1

2((pn·∇)w−(w·∇)pn,un)=0,

pn(0)=p0,n

(6b)

in (0, T)andforallv,w∈Wn. In order to show the existence of a solution to this system, we fix n∈Nand represent

each function pn∶[0,T]→Wnvia the corresponding basis of Wn,thatis,

pn(t)=

∑

i=1

Pni(t)wi.

We make use of the inverse J−1∶(L2

𝜎(Ω))∗→L2

𝜎(Ω) of the Riesz isomorphism on the Hilbert space L2

𝜎(Ω) in order

to define

un=Π

nJ−1gpn,(7)

where gpn∈(L2

𝜎(Ω))∗is given as

⟨gpn,v⟩=𝜇1(Δpn,ΔA−1v)−𝛾1(Δpn,A−1v)+𝜆1((pn·∇)pn,A−1v)

14044

EMMRICH AND GEUTER

for any pn∈Wnand v∈L2

𝜎(Ω). In order to transform the problem to an autonomous system of ordinary differential

equations, we let y∶[0,T]→Rnand y0∈Rnbe defined as

(y(t))i=Pni(t),(y0)i=(p0,n,wi)

for i=1,…,n. Then the discretized problem (6) can equivalently be written as

y′=f(y)on [0,T],

y(0)=y0,

where the right-hand side f∶Rn→Rnis given as

f(y)=−𝜇2(Δpn,Δw)−𝛾2(∇pn,∇w)−𝜆2((pn·∇)pn,w)−𝛼(|pn|2pn,w)

−𝛽(pn,w)−((un·∇)pn,w)−1

2((pn·∇)w−(w·∇)pn,un)

for i=1,…,n, as well as

pn=

∑

i=1

yiwi,

and unas in (7). Note that the mass matrix in this case is just the identity since the eigenfunctions wi

are orthonormal with respect to the L2-inner product. Then fdoes not depend on tand is continuous with

respect to y. Therefore, the Peano theorem (see Hale13, Chapter I, Theorem 5.1) provides a maximally continued solution

(un,pn)∈ 𝒞1([0,Tn); Wn)×𝒞1([0,Tn); Wn)to (6), where Tn≤Tmay depend on the dimension nof Wn.Wenow

derive a priori estimates to show that there is no blow-up of these solution; hence, they exist globally in time.

To this end, we first test Equation (6b) with pn.Sinceunand pnare both divergence-free, this results in

dt||pn||2

L2+𝜇2||Δpn||2

L2+𝛾2||∇pn||2

L2+𝛼||pn||4

L4+𝛽||pn||2

L2=0in(0,T).

Here, we also used the identity

(𝜕tpn,pn)=1

dt||pn||2

L2in (0,T)

for continuously differentiable functions. If 𝛾2,𝛽≥0, this already gives the desired a priori estimate for p.Weset

a−∶= {0,if a≥0

−a,if a>0.

for any a∈R. Performing an integration-by-parts and using Cauchy–Schwarz's as well as Young's inequality then

leads to

𝛾−

2||∇pn||2

L2≤(𝛾−

2)2

2𝜇2||pn||2

L2+𝜇2

2||Δpn||2

L2,

and by another application of Young's inequality as well as the continuous embedding L4(Ω) → L2(Ω), it follows that

((𝛾−

2)2

2𝜇2

+𝛽−)||pn||2

L2≤|Ω|

2𝛼((𝛾−

2)2

2𝜇2

+𝛽−)2

+𝛼

2||pn||4

L4.

By absorbing the norms of pninto the left-hand side, we arrive at

dt||pn||2

L2+𝜇2||Δpn||2

L2+𝛼||pn||4

L4≤|Ω|

𝛼((𝛾−

2)2

2𝜇2

+𝛽−)2

14045

Loading more pages...