# 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 $\rho \left({\partial }_{t}\mathbf{\text{u}}+\mathbf{\text{u}}\cdot \mathbf{\text{\nabla}}\mathbf{\text{u}}\right)=-\mathbf{\text{\nabla}}p+\mathbf{\text{\nabla}}\cdot \left(2\mu \mathbf{\text{D}}\right)+\rho \mathbf{\text{a}}+\mathbf{\text{\lambda}}\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}\Omega$ $\left(1-\frac{\rho }{{\rho }_{s}}\right)\left(M\left({\partial }_{t}\mathbf{\text{U}}-\mathbf{\text{g}}\right)\right)=-{\int }_{P\left(t\right)}\mathbf{\text{\lambda}}dv\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}P\left(t\right)$ $\left(1-\frac{\rho }{{\rho }_{s}}\right)\left(\mathbf{\text{I}}{\partial }_{t}\mathbf{\text{\omega}}+\mathbf{\text{\omega}}×\mathbf{\text{I}}\mathbf{\text{\omega}}\right)=-{\int }_{P\left(t\right)}\mathbf{\text{r}}×\mathbf{\text{\lambda}}dv\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}P\left(t\right)$ $\mathbf{\text{u}}-\left(\mathbf{\text{U}}+\mathbf{\text{\omega}}×\mathbf{\text{r}}\right)=0\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}P\left(t\right)$ $\mathbf{\text{\nabla}}\cdot \mathbf{\text{u}}=\mathbf{\text{0}}\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}\Omega$

with unknowns $\mathbf{\text{u}},\mathbf{\text{U}},\mathbf{\text{\omega}}$ and $\mathbf{\text{\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\left(t\right)$ with density ${\rho }_{s}$, inertia tensor $\mathbf{\text{I}}$ and mass $M$. The vector $\mathbf{\text{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:

• 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.

$\rho \left({\partial }_{t}\mathbf{\text{u}}+\mathbf{\text{u}}\cdot \mathbf{\text{\nabla}}\mathbf{\text{u}}\right)=-\mathbf{\text{\nabla}}p+\mathbf{\text{\nabla}}\cdot \left(2\mu \mathbf{\text{D}}\right)+\rho \mathbf{\text{a}}\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}\Omega$

$\mathbf{\text{\nabla}}\cdot \mathbf{\text{u}}=\mathbf{\text{0}}\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}\Omega .$

## 2nd subproblem: fictitious-domain problem

We are solving a fictitious-domain problem given by:

$\rho {\partial }_{t}\mathbf{\text{u}}=\mathbf{\text{\lambda}}\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}\Omega$ $\left(1-\frac{\rho }{{\rho }_{s}}\right)\left(M\left({\partial }_{t}\mathbf{\text{U}}-\mathbf{\text{g}}\right)\right)=-{\int }_{P\left(t\right)}\mathbf{\text{\lambda}}dv\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}P\left(t\right)$ $\left(1-\frac{\rho }{{\rho }_{s}}\right)\left(\mathbf{\text{I}}{\partial }_{t}\mathbf{\text{\omega}}+\mathbf{\text{\omega}}×\mathbf{\text{I}}\mathbf{\text{\omega}}\right)=-{\int }_{P\left(t\right)}\mathbf{\text{r}}×\mathbf{\text{\lambda}}dv\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}P\left(t\right)$ $\mathbf{\text{u}}-\left(\mathbf{\text{U}}+\mathbf{\text{\omega}}×\mathbf{\text{r}}\right)=0\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}P\left(t\right),$

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:

${\mathbf{\text{u}}}_{ns}^{n+1}={\mathbf{\text{u}}}_{ns}^{n}+\Delta t\left(-{\mathbf{\text{u}}}_{ns}^{n+1/2}\cdot \nabla {\mathbf{\text{u}}}_{ns}^{n+1/2}-\nabla {p}^{n+1}/\rho +\nabla \cdot \left[2\mu {\mathcal{\text{𝒟}}}_{v}^{n+1}\right]+{\mathbf{\text{a}}}^{n+1}\right),$

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

### Fictitious domain

The nunerical scheme for the fictitious domain problem reads: $\rho \left(\frac{{\mathbf{\text{u}}}_{fd}^{n+1}-{\mathbf{\text{u}}}_{ns}^{n+1}}{\Delta t}\right)-{\mathbf{\text{\lambda}}}^{n+1}=\mathbf{\text{0}}\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}\Omega$ $\left(1-\frac{\rho }{{\rho }_{s}}\right)\left(M\left(\left[\frac{{\mathbf{\text{U}}}^{n+1}-{\mathbf{\text{U}}}^{n}}{\Delta t}\right]-\mathbf{\text{g}}\right)\right)=-{\int }_{P\left(t\right)}{\mathbf{\text{\lambda}}}^{n+1}\phantom{\rule{0.333em}{0ex}}dv\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}P\left(t\right)$

$\left(1-\frac{\rho }{{\rho }_{s}}\right)\left(\mathbf{\text{I}}\left(\frac{{\mathbf{\text{\omega}}}^{n+1}-{\mathbf{\text{\omega}}}^{n}}{\Delta t}\right)+{\mathbf{\text{\omega}}}^{n+1}×\mathbf{\text{I}}{\mathbf{\text{\omega}}}^{n+1}\right)=-{\int }_{P\left(t\right)}\mathbf{\text{r}}×{\mathbf{\text{\lambda}}}^{n+1}\phantom{\rule{0.167em}{0ex}}dv\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}P\left(t\right)$ ${\mathbf{\text{u}}}_{fd}^{n+1}-\left({\mathbf{\text{U}}}^{n+1}+{\mathbf{\text{\omega}}}^{n+1}×\mathbf{\text{r}}\right)=0\phantom{\rule{0.333em}{0ex}}\mathrm{\text{over}}\phantom{\rule{0.333em}{0ex}}P\left(t\right).$