sandbox/cselcuk/README
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
Navier-Stokes
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).
Test cases
Convergences
Rotation
- Freely rotating cylinder in shear flow
- Stokes flow: freely rotating sphere in simple shear flow, shear along x, rotation along z
- Stokes flow: freely rotating sphere in simple shear flow, shear along y, rotation along x
- Stokes flow: freely rotating sphere in simple shear flow, shear along z, rotation along y
- Sphere with inital rotational velocity with gravity
Translation
Stencils and fictititous domains
- Colision test 2D
- Colision test 3D
- 2-2D overlapping fictitious-domains
- 3-2D overlapping fictitious-domains
- 3-2D Another overlapping fictitious-domains
Periodicity
- 2D-advection of a solid cylinder in the x-direction in a bi-periodic domain
- 2D-advection of a solid cylinder in the y-direction in a bi-periodic domain
- 2D-advection of a solid cylinder in the (x+y)-directions in a bi-periodic domain
- 3D-advection of a solid sphere in the (x+z)-directions in a tri-periodic domain
- 3D-advection of a solid sphere in the z-direction in a tri-periodic domain