Fictitious domain method

In numerical analysis, the fictitious domain method (FD) is a numerical technique designed to solve partial differential equations on complex geometries by embedding the physical domain into a larger and simpler computational domain. The method consists of extending the equations beyond the physical boundaries and enforcing the interface conditions through a distributed Lagrange multiplier in order to recover the correct solution within the original domain of interest.^[1]

Example of a fitted mesh for a circular obstacle

Example of a background unfitted mesh with an overlapping circular obstacle

This method belongs to the more general family of unfitted methods (also known as embedded or immersed), which allow solving interface problems on complex or evolving domains without generating a mesh that conforms to the domain’s boundaries.^[2] For this reason, the construction of two independent meshes is considered in the fictitious domain method: one for the unfitted background domain and another overlapping mesh that conforms to the interface, but which can move freely during the simulation. As with the immersed boundary method^[3] and the level-set method,^[4] the fictitious domain method is particularly effective for problems with evolving geometries since there is no need to rebuild the mesh but just to move the overlapping one.^[5]

The Dirichlet boundary conditions on immersed closed interfaces are imposed weakly^[6] by the introduction of the distributed Lagrange multiplier which enforces the boundary condition over the whole overlapping region.^[1]

A particular application of the fictitious domain method is in fluid–structure interaction (FSI) problems, where a fluid and an immersed solid interact one each other, exchanging forces through the interface.^[7] In this context, the solid is considered to be embedded within the computational fluid dynamics domain and the coupling conditions at the interface are enforced weakly over the whole fictitious fluid region.^[8]

Weak imposition of Dirichlet boundary conditions

Let ${\displaystyle \Omega _{f}\subset \mathbb {R} ^{d}}$ be a fluid domain in which we want to solve the Stokes equation for incompressible Newtonian fluids. We want to impose the Dirichlet boundary condition ${\displaystyle {\textbf {u}}={\textbf {g}}}$ on the boundary of an immersed circle described by the interface ${\displaystyle \Gamma =\partial \Omega _{c}}$. In this context, the fictitious domain method allows us to solve the Stokes equation over the whole domain ${\displaystyle \Omega =\Omega _{f}\cup \Omega _{c}}$, and to impose weakly the constraint ${\displaystyle {\textbf {u}}={\textbf {g}}}$ by means of the Lagrange multiplier ${\displaystyle \mathbf {\lambda } }$.

Velocity field around a fixed circular obstacle

The fictitious domain weak formulation reads:^[9]$${\displaystyle {\begin{cases}\displaystyle \int _{\Omega }\mu \nabla {\textbf {u}}:\nabla {\textbf {v}}-\int _{\Omega }p\nabla \cdot {\textbf {v}}+\langle \mathbf {\lambda } ,{\textbf {v}}\rangle _{\Omega _{c}}=0\quad &\forall {\textbf {v}}\in H_{0}^{1}(\Omega )\\[1ex]\displaystyle \int _{\Omega }q\ \nabla \cdot {\textbf {u}}_{f}=0\quad &\forall q\in L_{0}^{2}(\Omega )\\[1ex]\displaystyle \langle {\textbf {u}}-{\textbf {g}},\mathbf {\eta } \rangle _{\Omega _{c}}=0\quad &\forall \mathbf {\eta } \in \Lambda \end{cases}}}$$

The last equation of the system guarantees that, in a weak sense, the interface condition ${\displaystyle {\textbf {u}}={\textbf {g}}}$ is satisfied over the whole domain ${\displaystyle \Omega _{c}}$, and in particular on the interface ${\displaystyle \Gamma }$, while in ${\displaystyle \Omega _{f}}$ the equation is left untouched.^[6]

The Lagrange multiplier ${\displaystyle \mathbf {\lambda } }$ acts on the fluid as an external force defined only in ${\displaystyle \Omega _{c}}$ that enforces the fictitious fluid velocity ${\displaystyle {\textbf {u}}}$ to fulfil the constraint inside the circle.^[8]

Fluid-structure interaction application

As it is typically done in fluid-structure interaction problems, we consider a fluid domain ${\displaystyle \Omega _{f}\subset \mathbb {R} ^{d}}$, where the Navier–Stokes equations are usually solved, and ${\displaystyle {\mathcal {B}}\subset \mathbb {R} ^{d}}$ representing the region occupied by the structure.^[7] The fluid is described in an Eulerian framework while the solid is described in a Lagrangian formulation. This means that the solid domain ${\displaystyle {\mathcal {B}}}$ can be expressed as the image of a fixed reference configuration ${\displaystyle {\mathcal {B}}^{*}}$ through the Lagrangian map ${\displaystyle x={\text{X}}(s,t)}$.^[7]

The solid structure interacts with the fluid through its interface ${\displaystyle \Gamma =\partial {\mathcal {B}}}$, and from the mathematical point of view, it is required that the following kinematic condition is satisfied at the interface ${\displaystyle \Gamma }$:$${\displaystyle {\textbf {u}}(x,t)={\dfrac {\partial {\text{X}}(s,t)}{\partial t}}{\Big |}_{x={\text{X}}(s,t)}}$$where ${\displaystyle {\textbf {u}}}$ is the fluid velocity and ${\displaystyle x={\text{X}}(s,t)}$ is the position of the Lagrangian coordinate ${\displaystyle s}$ at time ${\displaystyle t}$. The latter kinematic condition can be seen as a specific Dirichlet boundary condition at the interface ${\displaystyle \Gamma }$ which involves both the unknowns ${\displaystyle {\textbf {u}}}$ and ${\displaystyle {\text{X}}}$ of the problem.^[7]

Therefore, in the spirit of the fictitious domain method, we can extend the fluid equation to the whole fixed computational domain ${\displaystyle \Omega =\Omega _{f}\cup {\mathcal {B}}}$ and impose weakly the kinematic coupling condition through the Lagrange multiplier ${\displaystyle \mathbf {\lambda } }$, requiring that the fictitious velocity ${\displaystyle {\textbf {u}}}$ matches the time derivative of the Lagrangian map ${\displaystyle {\text{X}}}$ in a weak sense inside the solid domain ${\displaystyle {\mathcal {B}}}$.^[7]

Numerical simulation of a heart valve

This requirement leads to the following weak formulation:^[7]$${\displaystyle {\begin{cases}a_{f}({\textbf {u}},{\textbf {v}})-b(p,{\textbf {v}})+\displaystyle \langle \mathbf {\lambda } ,{\textbf {v}}({\text{X}})\rangle _{{\mathcal {B}}^{*}}=F_{f}({\textbf {v}})\quad &\forall {\textbf {v}}\in H_{0}^{1}(\Omega )\\[1ex]b(q,{\textbf {u}})=0\quad &\forall q\in L_{0}^{2}(\Omega )\\[1ex]a_{s}({\text{X}},{\text{Y}})-\displaystyle \langle \mathbf {\lambda } ,{\text{Y}}\rangle _{{\mathcal {B}}^{*}}=F_{s}({\text{Y}})\quad &\forall {\text{Y}}\in H^{1}({\mathcal {B^{*}}})\\[1ex]\displaystyle \langle {\textbf {u}}({\text{X}})-{\dfrac {\partial {\text{X}}}{\partial t}},\mathbf {\eta } \rangle _{{\mathcal {B}}^{*}}=0\quad &\forall \mathbf {\eta } \in \Lambda \end{cases}}}$$

where ${\displaystyle a_{f}(\cdot ,\cdot )}$ and ${\displaystyle b(\cdot ,\cdot )}$ are the classical bilinear forms coming from the Navier–Stokes equations, while ${\displaystyle a_{s}(\cdot ,\cdot )}$ is the bilinear form associated with the elastic model describing the solid.^[7]

The kinematic condition is represented in a weak sense by the last equation, and the Lagrange multiplier ${\displaystyle \mathbf {\lambda } }$ acts as a body force on both the fluid and solid equations in order to fulfil this condition.^[8]

Numerical aspects

The fictitious domain method proves to be effective for the numerical simulation of moving obstacles, as it employs two independent meshes: a fixed background grid and an independent, overlapping mesh that describes the obstacle geometry.^[5] As we have seen in the previous applications, in order to weakly enforce the interface condition, we have to compute the duality pairing ${\textstyle \displaystyle \langle \mathbf {\lambda } ,{\textbf {v}}\rangle _{\Omega _{ob}}}$ over the obstacle domain ${\textstyle \Omega _{ob}}$.^[7][8][9]

The Lagrange multiplier ${\displaystyle \mathbf {\lambda } }$ is defined over ${\textstyle \Omega _{ob}}$ while the test function ${\displaystyle \mathbf {v} }$ is defined over the whole computational domain ${\displaystyle \Omega }$. Therefore, upon finite element method discretization, the discrete Lagrange multiplier ${\displaystyle \mathbf {\lambda } _{h}}$ is defined on the overlapping mesh, whereas the discrete test function ${\displaystyle \mathbf {v} _{h}}$ lives on the background mesh.^[10] As a result, computing the above coupling term requires the evaluation of basis functions defined on the background mesh at locations that belong to the overlapping mesh, which may be change over time. This task is computationally non-trivial and may introduce numerical errors and possibly wrong solutions.^[10]

This issue can be handled with different strategies:^[11]

• Augmented quadrature rule: sub-triangulate the overlapping region in order to adopt an enhanced quadrature formula capable of accurately computing the coupling term. While exact, this approach can be computationally expensive.^[11]
• Inexact quadrature rule: maintain the original quadrature used for the overlapping mesh, accepting a quadrature error that is controlled by the mesh size ${\displaystyle h}$.^[10]
• Collocation method: project ${\displaystyle \mathbf {\lambda } _{h}}$ over the discrete space of ${\displaystyle \mathbf {v} _{h}}$ in order to have both the functions defined on the same mesh.^[8]

See also

• Computational fluid dynamics
• Finite element method
• Fluid structure interaction
• Immersed boundary method
• Level set method
• Meshfree method