OpenFOAM の adjoint Shape Optimization Foam を考える #01手法

adjoint Shape Optimization Foam とは

流体力学やCFDに深い理解がない状態で他サイトの記事を読んだりコードを見て書いていますので間違っている可能性が多分にあります。
お気づきの点がありましたら教えてください。
また、よくわからないままメモとして書いている部分もあります。

構造最適化手法の一種に「トポロジー最適化」というものがあります。
ある設計空間にトポロジー最適化を適用することによって不要な部分が削られ、構造にとって必要な部分だけが残ります。
これを実装しているのがOpenFOAMのadjoint Shape Optimization Foamです。

元になっているのは、C.Othmerの「Implementation of a Continuous Adjoint for Topology Optimization of Ducted Flows」です。

本記事では「PhD course in CFD with OpenSource software」を参考に複数回に分けて、順番に手法と実装内容、またそのカスタマイズ方法まで考えたいと思います。

今回は元になっている方程式や手法についてです。
OpenFOAMは全く関係ない基礎的な部分になります。

そもそもトポロジー最適化とは何なのかについてはググれば日本語の資料もたくさん出てくるので今回は省略します。

目的関数と制約条件

ある関数 $J$ の最小化を目的にトポロジー最適化を行うと考えると、目的関数は以下のように書けます。

[1]　$minimize J = J(\textbf{v},p,\alpha)$

ここで、$\textbf{v}$ は速度、$p$ は圧力、 $\alpha$ は設計変数です。
また制約条件として今回は非圧縮性の定常のナビエストークス方程式を解くことを考えます。

[2-1]　$(\textbf{v}\cdot\nabla)\textbf{v}+\nabla p-\nabla(2\nu D(\textbf{v}))+\alpha\textbf{v}=0$

[2-2]　$\nabla\cdot\textbf{v}=0$

ここでは、ダルシー則を用いてシンク項として設計変数である porosity (空隙率) $\alpha$を導入しています。
コスト関数の値を上昇させるセルはこの porosity を増加させることによって影響を減らします。

ひずみ速度テンソル $D(\textbf{v})$は $D(\textbf{v})=\frac{1}{2}(\nabla\textbf{v}+(\nabla\textbf{v})^{T})$ を示しています。

【参考】ひずみ速度テンソル

対称テンソルと反対称テンソル [物理のかぎしっぽ]

動粘性係数 $\nu$ は分子粘性係数と渦粘性係数（乱流粘性係数）を足し合わせたものです。

ここから、この制約条件はレイノルズ応力にブシネスク近似を用いたレイノルズ平均ナビエストークス方程式だということが言えます。

ラグランジュ乗数を用いてラグランジュ緩和を行うことで目的関数と制約条件を一つの式にすると、

[3]　$L := J+ \int_\Omega {(\textbf{u}, q)} Rd\Omega$

と書けます。

ラグランジュの未定乗数法

ラグランジュの未定乗数法とは、

$L(x,y,\lambda)=f(x,y)-\lambda g(x,y)$
を作ったときに、
$\frac{\partial L}{\partial x}= \frac{\partial L}{\partial y}=\frac{\partial L}{\partial\lambda}=0$ の解、もしくは$\frac{\partial g}{\partial x}=\frac{\partial g}{\partial y}=0$ の解である $ (\alpha,\beta)$ が極値を与える

という、等式制約付きの関数最大化、最小化問題を解く手法のことです。

【参考】ラグランジュの未定乗数法

ラグランジュの未定乗数法と例題 _ 高校数学の美しい物語

ラグランジュの未定乗数法の解説と直感的な証明 _ yunabe.jp

解析領域

ここで、 $\Omega$ とは、 flow domain（流体が通る解析領域）を指します。
$(\textbf{u}, q)$ の組はラグランジュ乗数で、それぞれ adjoint velocity (随伴速度) と adjoint pressure (随伴圧力) を示します。

このようにして評価関数の勾配をLagrange乗数法を用いて求める方法をadjoint法（随伴変数法）と言います。
微分の連鎖律で勾配を求める直接法は多制約問題に有利であるのに対して、adjoint法は多設計変数問題に有利な手法です。

$R$ は非圧縮性のレイノルズ平均モデルと[2-2]のような連続の式を指します。

【参考】レイノルズ平均モデル

もっと知りたい！ 熱流体解析の基礎71　第7章 乱流計算：7.3.1 レイノルズ平均モデル (1)｜技術コラム

$(\textbf{u}, q)R$ は、それぞれの制約条件にラグランジュ乗数を掛けたものが足しあわされて、全領域にわたって合計されることを示します。

適切な adjoint velocity と adjoint pressure の組である $(\textbf{u}, q)$ を選ぶため、$(\textbf{v}, p)$ を式から消去します。

[4]　$ \delta L = \delta_{\alpha}L + \delta_{\textbf{v}}L + \delta_{p}L$

[5]　$\delta_{\textbf{v}}L + \delta_{p}L = 0$

設計変数についての求める sensitivity（解析感度）を得るために [4] の $L$ の全変動を調べます。

[4] [5] から、セル $i$ の porosity についてのコスト関数の感度は、

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

と表すことができます。
ここで、 $V_{i}$ とは、セル $i$ の体積を表します。

また、 [5] は次の [7] [8] のような adjointナビエストークス方程式に変換されます。

[7]　$-2D(\textbf{u})\text{v} = -\nabla q + \nabla \cdot (2\nu D(\textbf{u})) -\alpha\textbf{u} – \frac{\partial J_{\Omega}}{\partial\textbf{v}} $

[8]　$\nabla \cdot \textbf{u} = \frac{\partial J_{\Omega}}{\partial p}$

なお、 adjoint velocity や adjoint pressure は速度や圧力という名前を冠してはいますが、物理的な意味での速度や圧力として解釈すべき値ではありません。
主変数である速度や圧力との類似性の強調、さらに主変数を解くのと同様の手法が適用できることを示すためにこのような名前を用いています。

flow domain からの影響を受けず、境界からのみ影響を受けるコスト関数の場合、 $J_{\Omega}$ を含む項はゼロになります。
$ \Omega$ は境界を除く領域の体積部分を示しています。
したがってこのとき、adjoint 運動方程式 [7] とadjoint圧力方程式 [8] は変化しません。

[9]　$ \int_{\Gamma}( \textbf{n}( \textbf{u} \cdot \textbf{v} ) + \textbf{u}( \textbf{v} \cdot \textbf{n} ) + 2\nu\textbf{n} \cdot D( \textbf{u} ) – q \textbf{n} + \frac{\partial J_{\Gamma}}{ \partial \textbf{v}} ) \delta \textbf{v} d\Gamma – \int_{\Gamma} 2\nu \textbf{n} \cdot D(\delta \textbf{v}) \cdot \textbf{u} d\Gamma = 0 $

[10]　$ \int_{\Gamma} ( \textbf{u} \cdot \textbf{n} + \frac{\partial J_{\Gamma}}{\partial p} ) \delta p d \Gamma = 0 $

その代わりに境界領域 $ \Gamma$ を用いて、境界条件によってコスト関数に影響を与えます。
[9] [10] は、 [7] [8] と同じ条件下のものです。
これらの adjoint 方程式 [7]-[10] と元のナビエストークス方程式を合わせた形が最適化問題の一般形になります。

この adjoint 方程式の導出には「凍結乱流」として知られる近似も含んでいます。
凍結乱流を用いると、乱流領域では正しい値を出さなくなってしまいますが簡単に近似できるようになります。

まとめ

adjoint法というのは、ラグランジュの未定乗数法を用いてある制約条件の元の目的関数を解く手法でした。
現時点で集めた資料をまとめただけですが、ひとまずこういうものだと思ってこの程度にしておき、今後コードを読み進めてさらに理解を深めていきたいです。