MATHEMATICAL MODEL

sites of plumes have not been predicted. The convergence toward numerically stable and stationary plumes for low and high initial density of cells is revealed in [142].

In this chapter, a computational model based on the finite element method is pro-posed aiming at investigating the behavior of the two-dimensional chemotaxis–diffusion–

convection system.

5.2 Mathematical model

The mathematical model for the chemotaxis–diffusion–convection is proposed in [136]

and reads as follows:

ρ ∂u

∂t + u · ∇u



+ ∇p − µ∇²u = −n V_bg(ρ_b− ρ)j,

∇ · u = 0,

∂n

∂t + ∇ · [un − D_b∇n + S_dimr(c)n∇c] = 0,

∂c

∂t + ∇ · (uc − D_O∇c) = −n κr(c),

(5.1)

where u = (u, v) denotes the velocity field of water (solvent), p the pressure, ρ and µ the water density and viscosity, n the areal number density of bacteria (i.e., number of bacteria per unit area in a 2D space), V_b and ρ_b the volume and volumetric mass density of a bacterium, c the concentration of oxygen, V_bg(ρ_b − ρ)j the buoyancy force exerted by a bacterium on the fluid in the vertical direction (unit vector j), S_dim the dimensional chemotaxis sensitivity, D_b and D_O the bacterium and oxygen diffusivities, κ the bac-terium oxygen consumption rate, and r(c) the dimensionless cut-off function for oxygen concentration. The cut-off function r(c) is defined by the step function

r(c) = 1 if c > c^∗,

0 if c ≤ c^∗, (5.2)

where c^∗ = 0.3.

Bacteria are slightly denser than water and are diluted in the solvent, so that (ρ_b − ρ)/ρ  1 and n V_b  1, respectively. The consumption of oxygen is proportional to the bacterium population density n. Both n and c are advected by the fluid. When the oxygen concentration is lower than a threshold, the bacteria become quiescent [136], that is, they neither consume oxygen nor swim toward sites of higher oxygen concentrations. The dynamics of space filling, intercellular signaling, and quorum sensing are also neglected.

The system of equations (5.1) with the boundary conditions introduced in the previous papers (e.g. [135, 139, 142]) is solved in a two-dimensional rectangular container (Ω).

The top boundary (Γ_top) represents the interface between liquid and air. On the free surface the concentration of oxygen is equal to the air concentration of oxygen (c_air) and

5.2. MATHEMATICAL MODEL

Figure 5.1: Boundary conditions for the the chemotaxis–diffusion–convection system (5.1). The air–water interface, where the oxygen concentration is equal to that of air, is not crossed by bacteria; the fluid vertical velocity component equals zero and the fluid is assumed to be free of tangential stress. The container walls are impermeable to bacteria and oxygen; a no-slip condition is imposed.

5.3. COMPUTATIONAL MODEL

the free tangential stress condition as well as the absence of bacterium flux are prescribed (figure 5.1). It is therefore rational to prescribe the following boundary conditions

∂u

∂y = 0, v = 0, c = c_air,

S_dimn r(c) ∇c · n − D_b∇n · n = 0 on Γ_top,

(5.3)

where n is the unit outward normal vector. On the container walls (Γ_w), a no-slip bound-ary condition is prescribed and the fluxes of bacteria and oxygen are equal to zero:

u = 0, v = 0, ∇n · n = 0, ∇c · n = 0 on Γ_w. (5.4) A no-slip condition at the air–water interface would enable the formation of hydrody-namic instabilities. The effect of a moving boundary due to the advection caused by an incompressible fluid flow is to be explored.

5.3 Computational model

5.3.1 Scaling and setting for numerical simulations

The characteristic length is defined by the container height h and the characteristic bacterium density by the average of the initial bacterium population density defined as

n₀ := 1

|Ω|

Z

Ω

n₀(x, t) dx. (5.5)

This particular choice of the characteristic bacterium density gives an easy measure of the total number of bacteria in each simulation for different initial distributions of bacteria.

Dimensionless variables are defined as x⁰ = x

h, t⁰ = t

h²/D_b, n⁰ = n n₀, c⁰ = c

c_air, p⁰ = p

µD_b/h², u⁰ = u D_b/h.

(5.6)

Five dimensionless parameters given below characterize the hydrodynamic and chemo-taxis transport equations:

Pr_τ= ν

D_b, Ra_τ=g V_bn₀(ρ_b−ρ)L³

D_bµ , S=S_dimc_air D_b , H=κ n₀L²

c_airD_b, Le_τ=D_O D_b

(5.7)

where Pr_τ is the taxis Prandtl number, Ra_τ the taxis Rayleigh number (buoyancy-driven flow), Le_τ the taxis Lewis number, S the dimensionless chemotaxis sensitivity, and H the

5.3. COMPUTATIONAL MODEL

chemotaxis head. The taxis Prandtl, Rayleigh, and Lewis numbers are analogous to the respective heat and mass Prandtl, Rayleigh, and Lewis numbers in heat and mass transfer (table 5.2). The chemotaxis sensitivity (S) and head (H) characterize the chemotaxis sys-tem. Only Ra_τ and H depend on the characteristic length L and characteristic bacterium density n₀.

After removing the prime in dimensionless quantities, the set of dimensionless equa-tions becomes

∂u

∂t + u · ∇u − Pr_τ∇²u + Pr_τ∇p = −Ra_τPr_τn j,

∇ · u = 0,

∂n

∂t + u · ∇n − ∇²n + S ∇ · (r(c)n∇c) = 0,

∂c

∂t + u · ∇c − Le_τ∇²c = H n r(c).

(5.8)

The system (5.15) is solved in a rectangular domain Ω = [−`, `] × [0, 1] with the initial conditions:

u(x, 0) = u₀(x), n(x, 0) = n₀(x), c(x, 0) = c₀(x). (5.9) On the top of the domain Ω, the dimensionless boundary conditions are prescribed as

∂u

∂y = 0, v = 0, S r(c) n ∇c · n − D_b∇n · n = 0, c = 1, (5.10) while on the other boundaries, the dimensionless boundary conditions are imposed as follows

u = 0, v = 0, ∇n · n = 0, ∇c · n = 0. (5.11) In the following sections, the hydrodynamic system will refer to the first two equations in (5.15) and the chemotaxis system to the last two equations in (5.15).

5.3.2 Numerical methods

The governing equations in (5.15) are solved using the finite element method. The bi-quadratic quadrilateral elements for the primitive variables u and the bilinear quadrilateral elements for the primitive variables p are adopted to satisfy the LBB stability condition [84–87]. There are nine nodes in one biquadratic quadrilateral element and four nodes in one bilinear quadrilateral element. In each element, the function φ can be written as φ = Σ_iNⁱφ_i, where φ_iare the nodal unknowns.

To avoid the convective instability while solving the convection dominated flow equa-tions, a streamline upwind/Petrov-Galerkin (SUPG) method [143] is employed. The in-consistent Petrov-Galerkin weighted residual model developed in [144] yields the

follow-5.3. COMPUTATIONAL MODEL

In the previous inconsistent formulation, the test function W is only applied to the convective term so as to add a numerical stabilizing term. The main purpose of employing a SUPG model is to add an amount of streamline artificial viscosity. The test function W is therefore rewritten in two parts W = N + B. B is called the biased part. On each node i, the biased part is Bⁱ_j = τ u_j^∂N_∂xⁱ

j where j ∈ {1, 2} corresponds to the Cartesian coordinates and (u₁, u₂) = (u, v).

From the exact solution of the convection–diffusion equation in one dimension, fol-lowing the derivation in [144], the constant τ is determined as

τ = δ(γ₁) u h₁ + δ(γ₂) v h₂

2(u²+ v²) , (5.13)

where γ_j = ^uj_{2 k}^hj with (h₁, h₂) = (∆x, ∆y) being denoted as the grid sizes. Finally, the derived expression for δ(γ) is given below

δ(γ) =

For comparison purpose in section 5.4.4, the double diffusion system (5.23) and the Rayleigh–Bénard system (5.24) are solved using the software Freefem++ [1]. The code uses a finite element method based on the weak formulation of the problem. Taylor-Hood P2–P1 elements are chosen in FreeFem++. These elements are used together with a characteristic/Galerkin formulation to stabilize the convection terms.

5.3.3 Code validation for the coupled Navier-Stokes and Keller-Segel equations

In this validation study, the following dimensionless differential equations account-ing for the coupled Keller-Segel and incompressible viscous hydrodynamic equations are solved