# Can’s sandbox

You can find here documentations, examples and test cases related to my implementation of a fictitious-domain method within basilisk.

# Distributed Lagrange multiplier/fictitious domain method for rigid particle laden flows

We are looking for a numerical approximation for the solution of the following equations \displaystyle \rho\left(\partial_t{\mathbf{u}} +\mathbf{u}\cdot \mathbf{\nabla}\mathbf{u}\right) = -\mathbf{\nabla} p + \mathbf{\nabla} \cdot\left(2\mu\mathbf{D}\right) + \rho\mathbf{a} + \mathbf{\lambda}~\text{over}~\Omega \displaystyle \left(1-\frac{\rho}{\rho_s}\right)\left(M\left(\partial_t{\mathbf{U}}-\mathbf{g}\right)\right)=-\int_{P(t)} {\mathbf{\lambda}} dv~\text{over}~P(t) \displaystyle \left(1-\frac{\rho}{\rho_s}\right)\left(\mathbf{I}\partial_t{\mathbf{\omega}} + \mathbf{\omega} \times\mathbf{I}\mathbf{\omega}\right) = -\int_{P(t)}\mathbf{r}\times{\mathbf{\lambda}} dv~\text{over}~P(t) \displaystyle \mathbf{u}-\left(\mathbf{U}+\mathbf{\omega}\times \mathbf{r}\right)=0~\text{over}~P(t) \displaystyle \mathbf{\nabla}\cdot\mathbf{u} = \mathbf{0} ~\text{over}~\Omega

with unknowns \mathbf{u}, \mathbf{U}, {\mathbf{\omega}} and {\mathbf{\lambda}} being respectively the fluid velocity field, the particle’s translational velocity, the particle’s rotational velocity and the Lagrange multipliers. The particle occupies the domain P(t) with density \rho_s, inertia tensor \mathbf{I} and mass M. The vector \mathbf{g} is the gravity acceleration here.

From a numerical point of view, this is a complicated problem to solve direclty. One has to deal with three main difficulties:

• an advection-diffusion equation
• the imcompressibility contraint with the pressure p as unknown
• the particle’s solid body’s constraint with the Lagrange multipliers as unknow.

The classic strategy is to split (and decouple) the problem into subproblems and solve them successively. We chose here a two-steps time-spliting.

## 1st sub-problem: Navier-Stokes

We are using the centred(.h) solver for this problem.

\displaystyle \rho\left(\partial_t{\mathbf{u}}+\mathbf{u}\cdot {\mathbf{\nabla}}\mathbf{u}\right) = -{\mathbf{\nabla}}p + {\mathbf{\nabla}}\cdot\left(2\mu\mathbf{D}\right) + \rho\mathbf{a}~\text{over}~\Omega

\displaystyle {\mathbf{\nabla}}\cdot\mathbf{u} = \mathbf{0} ~\text{over}~\Omega.

## 2nd subproblem: fictitious-domain problem

We are solving a fictitious-domain problem given by:

\displaystyle \rho \partial_t\mathbf{u} = {\mathbf{\lambda}}~\text{over}~\Omega \displaystyle \left(1-\frac{\rho}{\rho_s}\right)\left(M\left(\partial_t\mathbf{U}-\mathbf{g}\right)\right)=-\int_{P(t)}{\mathbf{\lambda}} dv~\text{over}~P(t) \displaystyle \left(1-\frac{\rho}{\rho_s}\right)\left(\mathbf{I}\partial_t {\mathbf{\omega}} + {\mathbf{\omega}} \times\mathbf{I}{\mathbf{\omega}}\right) = -\int_{P(t)}\mathbf{r}\times{\mathbf{\lambda}} dv~\text{over}~P(t) \displaystyle \mathbf{u}-\left(\mathbf{U}+{\mathbf{\omega}}\times \mathbf{r}\right)=0~\text{over}~P(t),

This leads to a saddle-point problem which is solved with an iterative solver (Uzawa, conjugate gradient algorithm).

## Numerical scheme

Basilisk handles the resolution with its own temporal scheme. It reads:

\displaystyle \mathbf{u}_{ns}^{n+1} = \mathbf{u}_{ns}^{n} + \Delta t (-\mathbf{u}_{ns}^{n+1/2} \cdot \nabla \mathbf{u}_{ns}^{n+1/2}-\nabla p^{n+1}/\rho + \nabla \cdot [2 \mu \mathcal{D}^{n+1}_v] + \mathbf{a}^{n+1}),

\displaystyle \mathbf{\nabla}\cdot\mathbf{u}_{ns}^{{n+1}} =\mathbf{0}

### Fictitious domain

The nunerical scheme for the fictitious domain problem reads: \displaystyle \rho\left(\frac{\mathbf{u}_{fd}^{n+1}-\mathbf{u}_{ns}^{n+1}}{\Delta t}\right) - \mathbf{\lambda}^{n+1} = \mathbf{0} ~\text{over}~\Omega \displaystyle \left(1-\frac{\rho}{\rho_s}\right)\left(M\left(\left[\frac{\mathbf{U}^{n+1}-\mathbf{U}^n}{\Delta t}\right] -\mathbf{g}\right)\right)=-\int_{P(t)}\mathbf{\lambda}^{n+1}~dv~\text{over}~P(t)

\displaystyle \left(1-\frac{\rho}{\rho_s}\right)\left(\mathbf{I}\left(\frac{\mathbf{\omega}^{n+1}-\mathbf{\omega}^n}{\Delta t}\right) + \mathbf{\omega}^{n+1} \times\mathbf{I}\mathbf{\omega}^{n+1}\right) = -\int_{P(t)}\mathbf{r}\times\mathbf{\lambda}^{n+1}\, dv~\text{over}~P(t) \displaystyle \mathbf{u}_{fd}^{n+1}-\left(\mathbf{U}^{n+1}+\mathbf{\omega}^{n+1}\times \mathbf{r}\right)=0~\text{over}~P(t).