OpenFOAM の adjoint Shape Optimization Foam を考える #02境界条件

今回は、adjoint法を用いてトポロジー最適化を行う際の境界条件の設定について考えてみます。

前回の「OpenFOAM の adjoint Shape Optimization Foam を考える #01手法」に引き続き、「PhD course in CFD with OpenSource software」を大幅に参考にしています。

シリーズ

1. OpenFOAM の adjoint Shape Optimization Foam を考える #01手法
2. OpenFOAM の adjoint Shape Optimization Foam を考える #02境界条件 （今回）

速度と圧力の境界条件

速度と圧力の境界条件を以下のように設定した想定で考えます。

このとき、adjoint 速度 $u$、adjoint 圧力 $q$ についても考えてみます。
まず、壁面と流入部での境界条件は以下のようになります。

[1]　$\textbf{u}_t = 0$

[2]　$u_n = -\frac{\partial J_\Omega}{\partial p}$

[3]　$\textbf{n}\cdot\nabla q = 0$

添字「$t$」は “tangential” すなわち接線方向成分であること、添字「$n$」は “normal” すなわち法線方向成分であることを示します。

同様に、流出部分での境界条件はこうなります。

[4]　$q = \textbf{u}\cdot\textbf{v} + u_n v_n + \nu (\textbf{n}\cdot\nabla ) u_n + \frac{\partial J_\Gamma}{\partial v_n}$

[5]　$0 = v_n \textbf{u}_t + \nu (\textbf{n}\cdot\nabla ) \textbf{u}_t + \frac{\partial J\Gamma}{\partial \textbf{v}_t}$

$\frac{\partial J_\Gamma}{\partial v_n}$ と $\frac{\partial J_\Gamma}{\partial \textbf{v}_t}$ はそれぞれ、$\frac{\partial J_\Gamma}{\partial \textbf{v}}$ の法線成分と接線成分です。

流出部分での一様流

今回は、流出部分での流れが均一な一様流になるような条件を考えてみます。
コスト関数はこのようになります。

[6]　$J := \int_{\Gamma_O} { \frac{c}{2} (\textbf{v} – \textbf{v}^d)^2} d\Gamma_O$

ここで、 $c$ は、単位を統一するための定数です。
また、 $\textbf{v}^d$ は、流出部分に設定したい速度です。

これを先ほどの境界条件[1]～[5]に反映させると、壁面と流入部では、

[7]　$\textbf{u} = 0$

[8]　$\textbf{n}\cdot\nabla q = 0$

となり、流出部では、

[9]　$q = \textbf{u}\cdot\textbf{v} + u_n v_n + \nu (\textbf{n}\cdot\nabla ) u_n + c(v_n – v_n^d)$

[10]　$0 = v_n \textbf{u}_t + \nu (\textbf{n}\cdot\nabla ) \textbf{u}_t + c(\textbf{v}_t – \textbf{v}_t^d)$

となります。

最急降下法

この目的関数を満たすために、最急降下法を使います。
最急降下法とは、今自分がいる位置から見て、最も良さそうな方向に進む手法です。
今回は、目的関数の最小化が目標なので、今いる位置から最も急な下り坂（＝勾配）に向かって進みます。
式で表すと、

[11]　$\textbf{p}_k = -\nabla f (\textbf{x}_k)$

となります。
ここで、 $\textbf{x}_k$ は、 $k$ ステップ目の現在地、 $f$ は最小化したい関数です。

$\textbf{x}_k$ の微小変化によって、コスト関数 $f$ の値は小さくなります。

[12]　$\dfrac{\partial L}{\partial \alpha_{i}} = \textbf{u}_{i} \cdot \textbf{v}{i} V_{i}$

[12] で定義された感度について考えると、最急降下法を適用することで次のように porosity (空隙率) $\alpha$ を更新できます。

[13]　$\alpha_{n+1} = \alpha_n – \textbf{u}_{i} \cdot \textbf{v}_{i} V_{i} \delta$

上式において、 $\alpha_{n+1}$ は次のステップの porosity、 $\alpha_n$ は現時点での porosity であり、 $\delta$ はセルのステップ長さ、 $V_i$ はセルのステップ体積です。

OpenFOAM における、 adjointShapeOptimizationFoam での最急降下の実装でも、安定性を向上させるために緩和係数が使用されます。

実装の最終的な形式は、緩和係数 $\gamma$ を用いて、

[14]　$\alpha_{n+1} = \alpha_n (1-\gamma) + \gamma min( max(( \alpha_n – \textbf{u}_{i} \cdot \textbf{v}_{i} V_{i} \delta), 0), \alpha_{max})$

と書くことができます。
$\alpha_{max}$ は、自分で決めることができる porosity の上限値です。
porosity は負にならず、 $\alpha_{max}$ より大きくならないように制限するために最小関数minと最大関数maxを使っています。

また、 $\alpha_{max}$ は各ステップ長を制限し、積 $\textbf{u}_{i} \cdot \textbf{v}_{i} V_{i}$ の値が大きなマイナスの数になってしまったときの影響を低減する作用もあります。

まとめ

今回は比較的理解しやすい部分でした。
次は、ここまで見てきたことがOpenFOAMの中でどのようにコードとして表現されているのかを見てみます。