概要

　構造最適化は、あらかじめ決められた設計・境界条件から、性能を最大限に発揮できる最適な構造や形状を求める方法です。その中で、トポロジー最適化があります。

　土木分野（建築分野？）でも使われており、土木設計などに適応できると考えております。（どう適応すればいいかわかりませんが。。。）

今回の記事では、簡単に、トポロジー最適化について説明したいと思います。

参考文献

以下の書籍を参考にしています。

この書籍にはMatlab コードが付随しており、トポロジー最適化の計算ができます。

平均コンプライアンスと全ポテンシャルエネルギー

　
　平衡状態において、剛性を最大化するような適切な構造を求めることを考える。そのためには、平均コンプライアンスと呼ばれる量を定義する必要がある。
　
　あらかじめ指定された設計領域内 $D$ に存在する弾性体を考える。弾性体には、物体力 $\bf b$ が作用しており、境界（Dirichlet境界条件） $\Gamma_u$ は完全に固定され、境界（Neumann境界条件） $\Gamma_t$ には、表面力 $\bf t$ が作用しているとする。

仮想原理から、

\begin{align}
\int_D \mathcal{d}\Omega \ \ {\bf \varepsilon}({\bf v}): \chi_{\Omega} {\bf E} : {\bf \varepsilon}({\bf u})
= \int_D \mathcal{d}\Omega \ \ \chi_{\Omega} {\bf b} \cdot {\bf v} 
+ \int_{\Gamma_t} \mathcal{d}\Gamma \ \ \chi_{\Omega} {\bf t} \cdot {\bf v} 
 \\
 {\bf u} \in V^D \ , \  \forall{\bf v} \in V^D 
\end{align}

とかけ、双１次形式 $a({\bf u},{\bf v})$ と１次形式 $l({\bf v})$ を導入し

\begin{align}
a({\bf u},{\bf v}) &= l({\bf v}) \ , \ {\bf u} \in V^D \ , \  \forall{\bf v} \in V^D  \\
a({\bf u},{\bf v}) &= \int_D \mathcal{d}\Omega \ \ {\bf \varepsilon}({\bf v}): \chi_{\Omega} {\bf E} : {\bf \varepsilon}({\bf u}) \\
l({\bf v}) &= \int_D \mathcal{d}\Omega \ \ \chi_{\Omega} {\bf b} \cdot {\bf v} + \int_{\Gamma_t} \mathcal{d}\Gamma \ \ \chi_{\Omega} {\bf t} \cdot {\bf v} 
\end{align}

と書く。 $\bf E$ は弾性テンソルであり ${\bf \varepsilon}$ はひずみテンソルである。$\chi_{\Omega}$ は、特性関数を表し

\begin{align}
\chi_{\Omega}({\bf x}) = \left\{
\begin{array}{ll}
1 & ({\bf x} \in \Omega_d)\\
0 & ({\bf x} \in D \setminus \Omega_d )
\end{array}
\right.
\end{align}

と書ける。$\Omega_d$ は最適構造となる領域を表し、特性関数は$\Omega_d$ 以外の領域ではゼロになる。実際の計算においては、上記の特性関数を使うことはなく、密度を用いた緩和法(SIMP)を用いる。具体的な方法は、後々紹介する。

強制変位がゼロの場合

ここで、境界において強制変位はゼロであるとする。平均コンプライアンス$l_{mc}$を

\begin{align}
l_{mc}= l({\bf v}) &= \int_D \mathcal{d}\Omega \ \ \chi_{\Omega} {\bf b} \cdot {\bf v} + \int_{\Gamma_t} \mathcal{d}\Gamma \ \ \chi_{\Omega} {\bf t} \cdot {\bf v} 
\end{align}

と定義する。平均コンプライアンスを最小化することは、弾性体に物体力と表面力を作用させた状態において、最小の変形をする構造を得ることになり、これは剛性の最大化になる。

　また、${\bf v}={\bf u}$ ならば

\begin{align}
l({\bf u}) = a({\bf u},{\bf u}) = 2 \times \mbox{ひずみエネルギー}
\end{align}

となり平均コンプライアンスの最小化は、ひずみエネルギーの最小化と等価になる。さらに、全ポテンシャルエネルギーは、

\begin{align}
 F({\bf v})= \frac{1}{2}a({\bf v},{\bf v})-l({\bf v}) \ , \  \forall{\bf v} \in V^D 
\end{align}

となり、平衡状態において

\begin{align}
 F({\bf u}) = \underset{{\bf v} \in V^D }{\text{minimize}} F({\bf v})= \frac{1}{2}a({\bf u},{\bf u})-l({\bf u}) = -\frac{1}{2}l({\bf u}) 
\end{align}

となるので、平均コンプライアンスの最小化は、全ポテンシャルエネルギーの最大化である。

　従って、目的関数は、

\begin{align}
 \underset{\mbox{設計変数} }{\text{maximize}}　F({\bf u}) =  \underset{\mbox{設計変数} }{\text{maximize}}　\left\{\underset{{\bf v} \in V^D }{\text{minimize}} F({\bf v}) \right\}
\end{align}

で与えられる。

強制変位が作用する場合

　弾性体には、物体力 $\bf b$ と表面力 $\bf t$ が作用せず、強制変位${\bf g}$が作用するとする。平均コンプライアンス$l_{mc}$を

\begin{align}
l_{mc}= \int_{\Gamma_u} \mathcal{d}\Gamma \ \ {\bf t}({\bf u}) \cdot {\bf g} 
\end{align}

と定義する。平均コンプライアンスを最大化すると、強制変位による表面力${\bf t}({\bf u})$ が最大化され剛性も最大化される。よって、強制変位がゼロの場合と設定が逆になり、

\begin{align}
 \underset{\mbox{設計変数} }{\text{maximize}}　l_{mc}=  \underset{\mbox{設計変数} }{\text{maximize}}　\int_{\Gamma_u} \mathcal{d}\Gamma \ \ {\bf t}({\bf u}) \cdot {\bf g} =  \underset{\mbox{設計変数} }{\text{maximize}}　a({\bf u},{\bf u})
\end{align}

と書ける。しかし、全ポテンシャルエネルギーは、

\begin{align}
 F({\bf u}) = \underset{{\bf v} \in V^D }{\text{minimize}} F({\bf v})= a({\bf u},{\bf u})
\end{align}

であるので、平均コンプライアンスの最大化は、全ポテンシャルエネルギーの最大化である。従って、目的関数は、

\begin{align}
 \underset{\mbox{設計変数} }{\text{maximize}}　F({\bf u}) =  \underset{\mbox{設計変数} }{\text{maximize}}　\left\{\underset{{\bf v} \in V^D }{\text{minimize}} F({\bf v}) \right\}
= \underset{\mbox{設計変数} }{\text{maximize}}　a({\bf u},{\bf u})
\end{align}

で与えられ、強制変位がゼロの場合と設定が同じになる。

　全ポテンシャルエネルギーを目的関数とすれば、強制変位の場合にも表面力が作用する場合にも計算することができる。

最適化問題の定式化

　最適設計を行うための大きな目的は、必要としない余肉を落とす為にあるので、それを反映できる制約条件は必要となる。総体積に関する制約として、

\begin{align}
g = \int_D \mathcal{d}\Omega \rho - \Omega_0
\end{align}

と設定する。$\rho$ は体積密度で設定変数であり、$\Omega_0$は体積の上限である。密度は、ゼロから$1$までの範囲をとる。全ポテンシャルエネルギーを目的関数とすると、最適化問題は以下のように定式化される。

\begin{align}
 &\underset{\rho }{\text{maximize}}　F({\bf u}) \\

\mbox{制約条件} &\\
g &= \int_D \mathcal{d}\Omega \rho - \Omega_0 \leq 0 \\
a({\bf u},{\bf v}) &= l({\bf v}) \ , \ {\bf u} \in V^D \ , \  \forall{\bf v} \in V^D  \\
0 &\leq \rho \leq 1
\end{align}

　弾性体には、境界（Dirichlet境界条件） $\Gamma_u$ は完全に固定され、境界（Neumann境界条件） $\Gamma_t$ には、表面力 $\bf t$ が作用している場合の、KKT条件を考える。ラグランジアンを

\begin{align}
L = -F({\bf u})+ a({\bf u},{\bf v}) -l({\bf v}) 
+ \Lambda\left( \int_D \mathcal{d}\Omega \rho - \Omega_0\right)
+ \int_D \mathcal{d}\Omega\left[\lambda_0(-\rho)+ \lambda_1(\rho-1)\right]
\end{align}

と定義する。${\bf v},\Lambda,\lambda_0,\lambda_1$ はラグランジュ未定乗数である。変分を行うと、

\begin{align}
\delta L &= -\frac{\partial F({\bf u}) }{\partial {\bf u}} \delta {\bf u} 
-\frac{\partial F({\bf u}) }{\partial \rho} \delta \rho 
+\frac{\partial a({\bf u},{\bf v}) }{\partial {\bf u}} \delta {\bf u} \\
&+ \Lambda \int_D \mathcal{d}\Omega \delta \rho 
+ \int_D \mathcal{d}\Omega\left[\lambda_0(-\delta\rho)+ \lambda_1\delta\rho\right] \\

&=  -\frac{\partial a({\bf u},{\bf u}) }{\partial {\bf u}} \delta {\bf u} 
-\frac{1}{2}\int_D \mathcal{d}\Omega \ \ {\bf \varepsilon}({\bf u}): \frac{\partial \chi_{\Omega} {\bf E}  }{\partial \rho}: {\bf \varepsilon}({\bf u})\delta \rho 

+\frac{\partial a({\bf u},{\bf v}) }{\partial {\bf u}} \delta {\bf u} \\
&+ \Lambda \int_D \mathcal{d}\Omega \delta \rho 
+ \int_D \mathcal{d}\Omega\left[\lambda_0(-\delta\rho)+ \lambda_1\delta\rho\right] \\
\end{align}

　
自己随伴 ${\bf v}={\bf u}$ とすると、１項目と3項目は消えて、

\begin{align}
\delta L 
&=  -\frac{1}{2}\int_D \mathcal{d}\Omega \ \ {\bf \varepsilon}({\bf u}): \frac{\partial \chi_{\Omega} {\bf E}  }{\partial \rho}: {\bf \varepsilon}({\bf u})\delta \rho 
+ \Lambda \int_D \mathcal{d}\Omega \delta \rho 
+ \int_D \mathcal{d}\Omega\left[\lambda_0(-\delta\rho)+ \lambda_1\delta\rho\right] \\
\end{align}

また、

\begin{align}
\frac{\delta L}{\delta \rho} = 0
\end{align}

よりKKT条件

\begin{align}
& -\frac{1}{2} {\bf \varepsilon}({\bf u}): \frac{\partial \chi_{\Omega} {\bf E}  }{\partial \rho}: {\bf \varepsilon}({\bf u})+\Lambda-\lambda_0+\lambda_1 = 0　 \\

&a({\bf u},{\bf v}) = l({\bf v}) \ , \ {\bf u} \in V^D \ , \  \forall{\bf v} \in V^D  \\

&\Lambda \geq 0 \ , \  g=\int_D \mathcal{d}\Omega \rho - \Omega_0 \leq 0 \ , \
\Lambda\left( \int_D \mathcal{d}\Omega \rho - \Omega_0\right) = 0 \\

&\lambda_0 \geq 0 \ , \ -\rho \leq 0  \ , \ \lambda_0(-\rho)=0 \\
&\lambda_1 \geq 0 \ , \ \rho-1 \leq 0  \ , \ \lambda_1(\rho-1)=0 \\

\end{align}

となる。

最適性基準法

　上記で説明したKKT条件から最適解を得る方法として、最適性基準法がある。以下、下付き添え字$i,j,k,\cdots$を設計変数 $i,j,k,\cdots$番目の番号とする。側面制約

\begin{align} 
\ -\rho_i \leq 0  \ , \ \rho_i-1 \leq 0  
\end{align}

が非活性である、つまり、等号が成り立っていないとすると

\begin{align}
&\rho_i \neq 0  \ , \ \lambda_{0i}(-\rho_i)=0 \Rightarrow \lambda_{0i} = 0 \\
&\rho_i \neq 1  \ , \ \lambda_{1i}(\rho_i-1)=0 \Rightarrow \lambda_{1i} =0 \\
\end{align}

となりラグランジュ未定乗数$\lambda_{0i},\lambda_{1i}$ はゼロになる。すると、

\begin{align}
 -\frac{1}{2} {\bf \varepsilon}({\bf u}): \frac{\partial \chi_{\Omega} {\bf E}  }{\partial \rho_i}: {\bf \varepsilon}({\bf u})+\Lambda_i &= 0 \\
\therefore \ \ \frac{{\bf \varepsilon}({\bf u}): \frac{\partial \chi_{\Omega} {\bf E}  }{\partial \rho_i}: {\bf \varepsilon}({\bf u})}{\Lambda_i}　 &= 1
\end{align}

が得られる。$2$は$\Lambda$ に取り込んだ。上式を満たすように設計変数とラグランジュ未定乗数が得られれば最適解の候補になる。

　最適解を探索する方法として、上付き添え字をステップ数とし、

\begin{align}
A_i^{(k)}&=\frac{{\bf \varepsilon}({\bf u}): \frac{\partial \chi_{\Omega} {\bf E}}{\partial \rho_i^{(k)}}: {\bf \varepsilon}({\bf u})}{\Lambda_i^{(k)}}　 \\

\rho_i^{(k+1)} &= \left(A_i^{(k)} \right)^{\eta} \rho_i^{(k)}
\end{align}

と設計変数を更新していく。$A_i^{(k)}$ は１になるように更新されるので、設計変数は最適解の候補となる。$\Lambda^{(k)}_i$ の更新の方法は、２分法を使い決定していく。

　側面制約が非活性であることを仮定しているので、設計変数の上限値と下限値を超えてしまうことがある。上限値と下限値の間に収まるように、ムーブリミット$\varsigma$を設けて

\begin{align}
\rho_i^{L(k)} &= \max \{\rho_i^{(k)}-\varsigma,0 \}  \\
\rho_i^{U(k)} &= \min \{\rho_i^{(k)}+\varsigma,1 \}  \\
\end{align}

上限値と下限値を更新し、最終的に設計変数を

\begin{align}
\rho_i^{(k+1)} = 
\begin{cases}
\rho_i^{L(k)} & \left(A_i^{(k)} \right)^{\eta}\rho_i^{(k)}  \leq \rho_i^{L(k)} \\
\left(A_i^{(k)} \right)^{\eta} \rho_i^{(k)} &   \rho_i^{L(k)} <  \left(A_i^{(k)} \right)^{\eta}\rho_i^{(k)} <  \rho_i^{U(k)}   \\
\rho_i^{U(k)}  &  \rho_i^{U(k)} \leq  \left(A_i^{(k)} \right)^{\eta} \rho_i^{(k)}
\end{cases}
\end{align}

と更新する。

設計変数の設定

　設計変数は、各要素に設定する場合と各節点に設定する方法がある。設定変数を各要素に設定すると、隣り合う要素同士では、互いに独立した値を取り得るため、チェッカーボードパターンを許容する。（チェッカーボードのような最適解になる。）
　設定変数を各節点に設定するとチェッカーボードパターンは回避できるが、グレースケールを生じやすい。今回の記事では、設計変数を節点に設定する方法を紹介する。

　FEM（有限要素法）で計算する際、各要素において密度を計算する必要があり、各節点における密度の値を使い要素における密度を補間する。その際は、FEMの計算で使われる形状関数 $N_j$ を用いて
　

\begin{align}
\sum_{j=1}^{N_e} N_j\rho_j
\end{align}

と補間する。これは、FEMの計算において要素剛性行列${\bf K_e}$を計算する際に、各要素の任意の点$(\xi,\eta)$のヤング率$E$を、設計変数である密度を使い

\begin{align}
E(\xi,\eta) = \left(\sum_{j=1}^{N_e} N_j\rho_j  \right)^pE^0 
\end{align}

と計算する。$p$ は、ペナライゼーションパラメーターである。平面ひずみを仮定する場合、弾性テンソルは、

\begin{align}
{\bf E} = \left(\sum_{j=1}^{N_e} N_j\rho_j\right)^p
\frac{E^0}{1-\nu^2}

\begin{pmatrix}
1   & \nu &  0 \\
\nu & 1   &  0 \\
0   & 0   & \frac{1-\nu}{2}  \\
\end{pmatrix}
\end{align}

として計算を行う。

計算手順

Step1 ： FEMの計算

　境界条件などを設定し、FEMの計算し変位$U$を求める。
FEMの計算は、密度を考慮した要素剛性行列${\bf K_e}$と要素荷重節点ベクトル${\bf F_e}$
　

\begin{align}
{\bf K_e} &= \int_{\Omega_e} \mathcal{d}\Omega \ \ {\bf B}^T{\bf E}{\bf B} \\
{\bf F_e} &= \int_{\Omega_e} \mathcal{d}\Omega \ \ {\bf N}^T {\bf b} + \int_{\Gamma_e} \mathcal{d}\Gamma \ \ {\bf N}^T {\bf t} 
\end{align}

から、全体の剛性行列${\bf K}$と荷重節点ベクトル${\bf F}$をつくり

\begin{align}
{\bf K}{\bf U} = {\bf F}
\end{align}

上記の連立方程式を解き、変位${\bf U}$ を求める。${\bf B}$ は形状関数の微分である。

Step２：全ポテンシャルエネルギーの計算

目的関数である全ポテンシャルエネルギーは、

\begin{align}
F({\bf U}) = -\frac{1}{2}l_{mc} = -\frac{1}{2}{\bf F}^T{\bf U}
= -\frac{1}{2} {\bf U}^T{\bf K}{\bf U} 
\end{align}

である。ステップ$k+1$番目と$k$番目の全ポテンシャルエネルギーの差が、ある閾値$\epsilon$より小さい、つまり、

\begin{align}
\|F({\bf U})^{(k+1)} -  F({\bf U})^{(k)} \| \leq \epsilon
\end{align}

ならば収束したと判断する。

Step３：設計感度の計算

　Step2で収束しない場合、最適性基準法を使い、設計変数(密度)を更新しなければならない。そのためには以下の量を計算する必要がある。

\begin{align}
{\bf \varepsilon}({\bf u}): \frac{\partial \chi_{\Omega} {\bf E}}{\partial \rho_i^{(k)}}: {\bf \varepsilon}({\bf u})
\end{align}

剛性行列を使い、上記の式を表すと

\begin{align}
{\bf U}^T \frac{\partial {\bf K}}{\partial \rho_i} {\bf U} 
\end{align}

と書け、節点$i$が含まれる要素$e$の要素剛性行列に対して

\begin{align}
\frac{\partial {\bf K_e}}{\partial \rho_i} &= \int_{\Omega_e} \mathcal{d}\Omega \ \ {\bf B}^T\frac{\partial{\bf E}　}{\partial \rho_i} {\bf B} 
\\

\frac{\partial{\bf E}　}{\partial \rho_i} &= pN_i\left(\sum_{j=1}^{N_e} N_j\rho_j\right)^{p-1}
\frac{E^0}{1-\nu^2}

\begin{pmatrix}
1   & \nu &  0 \\
\nu & 1   &  0 \\
0   & 0   & \frac{1-\nu}{2}  \\
\end{pmatrix}
\end{align}

を計算する。

Step４：最適性基準法による設計変数の更新

　最適性基準法を使い、設計変数(密度)を更新する。Step1 に戻り、収束するまで計算を行う。

計算結果

　シミュレーションの計算領域は、横幅 $1.2m$、高さ $0.8m$、格子間隔 $0.01m$ と設定した。物性値の値は、ヤング率 $1.0 GPa$、ポアソン比 $0.3$ とし、ペナライゼーションパラメーターは $p=3$ とした。
　
　境界条件は、左側面の変位はx,y方向どちらも固定する。そして、右側面の中心部から下方向に $1.0Gpa$ で引っ張る。体積の上限値は $\Omega_0=0.3$ とした。また、密度の下限値をゼロにしてしまうと計算ができなくなる可能性があるので、下限値は $0.01$ とした。以下の図は、約200ステップの計算結果である。

ヒートマップは、密度の値を表す。

まとめ

　今回の記事では、トポロジー最適化を用いて、全エネルギーポテンシャルが最大となるような最適設計をおこなった。設計変数である密度を使い、最適な構造や形状・形態を表現した。

密度ではなく、レベルセット関数を用いた形状表現も可能であり、機会があればレベルセット法によるトポロジー最適化の方法も紹介したい。