Optimal Shape Design for Poisson Equation in OpenFOAM

Download Code

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);
    diffusivity     	quadratic 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  	quadratic inverseFaceDistance (deformedWall);
    // diffusivity  	quadratic inverseDistance (deformedWall);
    diffusivity  	quadratic inversePointDistance (deformedWall);
    // 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.

Work Packages ≻ Optimal Control in OpenFOAM ≻

Optimal Shape Design for Poisson Equation in OpenFOAM

Mesh Motion Solver

Linear Elasticity

Solid Body Rotation (SBR) Stress

Laplacian Smoothing

Getting Started

Dynamic Mesh

Prerequisites

Running a Case

Warning

References

Related Posts

Interior Control of the Poisson Equation with the Conjugate Gradient Method in OpenFOAM

Optimal Shape Design in External Flow

Interior Control of the Heat Equation with the Steepest Descent Method and OpenFOAM

Interior Control of the Poisson Equation with the Steepest Descent Method in OpenFOAM