# Optimal Shape Design for Poisson Equation in OpenFOAM

Author: - 23 July 2018

We study the shape design problem through the minimization of the cost functional

$\mathcal{J} \left( \theta, y \right) = \frac{1}{2} \int_{\Omega \left( \theta \right )} \left( y - y_d \right) ^2 \, \mathrm{d} \Omega,$

where $y \in L^2 \left( \Omega \right)$ is the state variable, $y_d \in L^2 \left( \Omega \right)$ is a target function, and $\theta \left( \mathbf{x} \right)$ the normal displacement to a reference boundary,

$\Gamma \left( \theta \right) = \left\{ \mathbf{x} + \theta \left( \mathbf{x} \right) \mathbf{n} \left( \mathbf{x} \right) : \mathbf{x} \in \Gamma_{0} \right\}.$

The problem is subject to the following elliptic PDE in the domain $\Omega \left( \theta \right ) \subset \mathbb{R}^d$ with Dirichlet boundary conditions on $\Gamma \left( \theta \right ) = \Gamma_w \cup \Gamma_s \left( \theta \right )$,

$\begin{cases} -\Delta y = f & \text{in } \Omega \left( \theta \right ),\\ y = y_{w} & \text{in } \Gamma_w,\\ y = y_{s} & \text{in } \Gamma_s \left( \theta \right ), \end{cases}$

and $f \in L^2 \left( \Omega \right)$.

We consider the set of admissible domains whose measure is fixed,

$\mathcal{U}_{ad} = \left\{ \Omega \left( \theta \right ) : \lvert \Omega \left( \theta \right ) \rvert = \Omega_0 \right\},$

and aim at finding the optimal shape that minimizes the cost function. In order to do so, we use the steepest descent method with descent direction given by

$\delta \theta ^{\left( n \right)} = - \left[ \frac{1}{2} \left( y^{\left( n \right)} - y_d \right) ^2 + \frac{\partial v}{\partial n}^{\left( n \right)} \frac{\partial y}{\partial n}^{\left( n \right)} + q^{\left( n \right)} \right],$

where $v$ is solution of the adjoint problem

$\begin{cases} - \Delta v = y - y_d & \text{in } \Omega \left( \theta \right ), \\ v = 0 & \text{on } \Gamma \left( \theta \right ) = \Gamma_w \cup \Gamma_s \left( \theta \right ). \end{cases}$

The normal displacement field is updated at every iteration,

$\theta^{\left( n + 1 \right)} = \theta^{\left( n \right)} + \epsilon \delta{\theta^{\left( n \right)}},$

with $\epsilon$ sufficiently small. The Lagrange multiplier is computed in order to ensure that the volume contraint is fulfilled, thus

$q^{\left( n \right)} = - \frac{1}{\left|\Gamma_s \left( \theta \right )\right|} \int_{\Gamma_s \left( \theta \right )} \left[ \frac{1}{2} \left( y^{\left( n \right)} - y_d \right) ^2 + \frac{\partial v}{\partial n}^{\left( n \right)} \frac{\partial y}{\partial n}^{\left( n \right)} \right] \mathrm{d} \Gamma.$

## Mesh Motion Solver

The solutions to the primal and adjoint problems are commonly approximated by means of numerical methods, such as the finite element method (FEM) or the finite volume method (FVM). In order to apply these techniques, the domain under study must be tessellated with a mesh. Nevertheless, applying the displacement directly to the controlled boundary will deteriorate the surrounding elements after a few iterations and the computation will crash if the interior nodes of the domain are not reallocated. In order to avoid this, the domain can be re-meshed after a number of iterations. However, this can be very expensive as a completely new mesh must be generated. A commonly used alternative is to move the interior nodes of the mesh according to the displacements prescribed on the boundary. By doing this the number of elements and the nodes connectivities remain the same, only the mesh nodes positions are updated. The Solid Body Rotation Stress method and the Laplacian smoothing included in the OpenFOAM library have been used.

### Linear Elasticity

Mesh motion can be achieved by treating the mesh as an elastic body and solving the equations of Linear Elasticity for solids with prescribed displacements on the domain boundary,

$\begin{cases} \partial_j \left( G \partial_j u_i \right) + \partial_j \left( G \partial_i u_j \right) + \partial_i \left( \lambda \partial_j u_j \right) = 0 & \text{in } \Omega \left( \theta \right ), \\ u_i = 0 & \text{on } \Gamma_w, \\ u_i = \delta \theta n_i & \text{on } \Gamma_s \left( \theta \right ). \end{cases}$

### Solid Body Rotation (SBR) Stress

The Linear Elasticity model fails when dealing with rotating meshes. This can be mitigated by selecting the material properties in a proper manner,

$\begin{cases} 2 \partial_j \left( \gamma \left( \mathbf{x} \right) \partial_j u_i \right) + \partial_j \left[ \gamma \left( \mathbf{x} \right) \left( \partial_i u_j - \partial_j u_i - \delta_{ij} \partial_k u_k \right) \right] = 0, & \text{in } \Omega \left( \theta \right ), \\ u_i = 0 & \text{on } \Gamma_w, \\ u_i = \delta \theta n_i & \text{on } \Gamma_s \left( \theta \right ). \end{cases}$

### Laplacian Smoothing

A simple and widely used practice is to solve a Laplace equation with the prescribed displacements as boundary conditions,

$\begin{cases} \Delta u_i = 0 & \text{in } \Omega \left( \theta \right ), \\ u_i = 0 & \text{on } \Gamma_w, \\ u_i = \delta \theta n_i & \text{on } \Gamma_s \left( \theta \right ). \end{cases}$

The behavior of the Laplacian smoothing can be improved by adding a non-uniform diffusivity term,

$\partial_k \left( \gamma \left( \mathbf{x} \right) \partial_k u_i \right) = 0 \quad \text{in } \Omega \left( \theta \right ).$

The diffusion field $\gamma \left( \mathbf{x} \right): \mathbb{R}^d \rightarrow \mathbb{R}$ decreases with the distance to the controlled boundary. A common choice is to make it depend on the inverse of the distance as

$\gamma \left( \mathbf{x} \right) = \frac{1}{d^m},$

so that the nodes next to the deforming boundaries move with similar displacements as those on the boundary.

## Getting Started

The solver must be compiled in the terminal. It is advisable to first clean previous compilations with

wclean


and then use

wmake


### Dynamic Mesh

The mesh motion solver is specified in the dictionary constant/dynamicMeshDict.

• For the Laplacian solver:
dynamicFvMesh   dynamicMotionSolverFvMesh;

motionSolverLibs ( "libfvMotionSolvers.so" );

solver          displacementLaplacian;

displacementLaplacianCoeffs
{
//diffusivity  	uniform;
//diffusivity     	inversePointDistance (deformedWall);
}

• For SBR Stress method:
motionSolver 	displacementSBRStress;

displacementSBRStressCoeffs
{
// diffusivity  	uniform;
// diffusivity  	directional (1 200 0);
// diffusivity  	motionDirectional (1 1000 0);
// diffusivity  	file motionDiffusivity;
// diffusivity	inversePointDistance (deformedWall);
}


### Prerequisites

OpenFOAM C++ library must be installed in order to compile the code.

The OpenFOAM distribution provided by the OpenFOAM Foundation was used.

## Running a Case

In order to run the solver move to the case folder poissonOptShapeFoamCase and type in the command line

./Allprepare

poissonOptShapeFoam


The shape optimization method described in the previous section has been tested with a simple example. The Poisson equation is posed in a two-dimensional circular domain with boundary

$\Gamma_w = \left\{ \left( x, y \right) \in \mathbb{R}^2: \sqrt{x^2 + y^2} = R_w \right\}.$

The reference geometry to be optimized is an inner hole with boundary given by

$\Gamma_s \left( 0\right) = \left\{ \left( x, y \right) \in \mathbb{R}^2 : \sqrt{\left( x - c_x \right)^2 + \left( y - c_y \right)^2} = R_s \right\}.$

The target function $u_d$ will be the analytical solution of the Poisson equation with the inner hole centered in the origin,

$u_{analytic}\left( r \right) = -\frac{f}{4}r^2 + c_1 \log r + c_2,$

where the integration constants are obtained from the boundary conditions,

\begin{align*} c_1 = & \frac{u_w - u_s + \displaystyle \frac{f}{4} \left( R_w^2 - R_s^2 \right)}{\log\left( \displaystyle \frac{R_w}{R_s} \right)}, \\ c_2 = & \frac{u_w + u_s}{2} + \frac{f}{8}\left( R_w^2 + R_s^2 \right) - \frac{c_1}{2} \log\left( R_w R_s \right). \end{align*}

When the hole center is displaced from the origin, the target state must be extended in the region $r \leq R_1$, thus

$u_d\left( r \right)= \begin{cases} u_s, \quad & r \leq R_s, \\ u_{analytic}, \quad & R_s < r \leq R_w. \end{cases}$

It is clear that for

$\Gamma^{*}_s = \left\{ \left( x, y \right) \in \mathbb{R}^2 : \sqrt{x^2 + y^2} = R_s \right\}$

the cost function equals zero, thus it is an optimal solution. The steepest descent algorithm has been coded in the OpenFOAM solver poissonOptShapeFoam with $\epsilon = 10^{-3}$.

### Warning

It might be needed to use

sed -i -e 's/\r\$//' filename


and

chmod +x filename


in order to be able to execute

./filename


## References

• O. Pironneau. Optimal shape design for elliptic systems. Springer Science & Business Media, 2012.
• J Simon. Diferenciación de problemas de contorno respecto del dominio. Technical report, Universidad de Sevilla, Facultad de Matemáticas, Departamento de Análisis Matemático, 1989.
• The OpenFOAM Foundation.
• Moving boundary problem based on calculated data, CFDonline, 2013.

Linear Elasticity mesh motion method:

• Andrew A. Johnson and Tayfun E. Tezduyar. Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces. Computer methods in applied mechanics and engineering, 119(1-2):73–94, 1994.
• Thomas D. Economon, Francisco Palacios, and Juan J. Alonso. Unsteady continuous adjoint approach for aerodynamic design on dynamic meshes. AIAA Journal, 53(9):2437–2453, 2015.
• George S. Eleftheriou and Guillaume Pierrot. Rigid motion mesh morpher: A novel approach for mesh deformation. In OPT-i, An International Conference on Engineering and Applied Sciences Optimization, Kos, Greece, pages 4–6, 2014.

Solid Body Rotation (SBR) Stress method:

• Richard P. Dwight. Robust mesh deformation using the linear elasticity equations. In Computational fluid dynamics 2006, pages 401–406. Springer, 2009.

Laplacian Smoothing:

• Peter Hansbo. Generalized laplacian smoothing of unstructured grids. International Journal for Numerical Methods in Biomedical Engineering, 11(5):455–464, 1995.
• Rainald Löhner and Chi Yang. Improved ale mesh velocities for moving bodies. Communications in numerical methods in engineering, 12(10):599–608, 1996.
• Hrvoje Jasak and Zeljko Tukovic. Automatic mesh motion for the unstructured finite volume method. Transactions of FAMENA, 30(2):1–20, 2006.