白黒羊飼いの庭

unityで作るゲームの製作メモでしたが、最近は機械学習とかデザインとか何でもやってます

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

移転しました。
今後の更新は下記で行います。

www.shirokurohitsuji.studio






adjoint Shape Optimization Foam とは

久しぶりにOpenFOAMの話題です。

流体力学や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は全く関係ない基礎的な部分になります。

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


Adjoint Shape Optimization Foam の計算の流れ

目的関数と制約条件

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


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


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


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

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


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

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

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

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


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

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

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


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


と書けます。

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

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

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

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

解析領域

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

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

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

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

コラム:乱流|投稿一覧



(\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]  \displaystyle{ \dfrac{\partial L}{\partial \alpha_{i}} = \textbf{u}_{i} \cdot \textbf{v}_{i} V_{i}}


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

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


[7]  \displaystyle{ -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]  \displaystyle{ \nabla \cdot \textbf{u} = \frac{\partial J_{\Omega}}{\partial p}}


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

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


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