概要
構造最適化は、あらかじめ決められた設計・境界条件から、性能を最大限に発揮できる最適な構造や形状を求める方法です。その中で、トポロジー最適化があります。
土木分野でも研究されており、土木設計などに適応できると考えております。例えば、コンクリートを生成する際に大量の二酸化炭素を放出します。最適なコンクリート構造物を設計することで（不要なコンクリート部分を削減することで）、コンクリートの生成量を抑えることができ、二酸化炭素を減らすことができます。
他には、コンクリートは圧縮力に強いが引張力に弱いことが知られており、引張応力が生じる箇所には steel を挿入して引張応力に耐えられるように設計します。トポロジー最適化を用いることで steel を挿入する最適な箇所（引張応力が生じる箇所）を知ることができます。
前回の記事では、設計変数を節点に設定したトポロジー最適化を紹介しました。
今回の記事では、設計変数を要素に設定したトポロジー最適化について紹介します。
参考文献
以下の書籍を参考にしています。
最適化問題の簡単な例
横幅 $L=6$m, 高さ $x$ m ,幅 $b=x/3$m の構造物があり、構造物上部から等分布荷重 $q = 5 kN/m^2$ が掛かっているとする。構造物の体積が最小となるように高さ $x$ を決めたいとする。ただし、許容できるたわみは最大で $u_{\mbox{max}} = 0.03$ m とする。
均質材料の場合、たわみ$u(x)$は
\begin{align}
u(x) = \frac{5}{384}\frac{qL^4}{EI} \ ,\ I = \frac{bx^3}{12} =\frac{x^4}{36}
\end{align}
と書くことができる。$E$ はヤング率、$I$ は慣性モーメントである。許容できるたわみを最大 $u_{\mbox{max}} = 0.03$ m とすると、このときの高さ$x_{\mbox{max}}$は
\begin{align}
x_{\mbox{max}} &= \left(\frac{5\times36}{384}\frac{qL^4}{Eu_{\mbox{max}}} \right)^{1/4} \\
&= 0.241 \mbox{m}
\end{align}
なので、高さは $0.241$ m より低くなってはならない。
この問題を最適化問題として定式化すると、目的関数$f(x)$を構造物の体積として
\begin{align}
f(x) = bxL=\frac{x^2L}{3}
\end{align}
不等式制約 $g(x)$
\begin{align}
g(x) = 0.241-x \leq 0
\end{align}
を満たすように最適化問題を解くことになる。目的関数から、$x$ のとりうる範囲のなかで最小となる値を選べばよく、最小となる体積は、$x =0.241 m$ のとき $0.116m^3$ である。
FEMを用いたトポロジー最適化の定式化
目的関数の最適化は、平均コンプライアンスの最小化として定義され、これは剛性最大化を意味する。FEMの基礎方程式は、
\begin{align}
{\bf K}(x) {\bf U} ={\bf F}
\end{align}
と書くことができる。${\bf K}$は剛性行列であり、 ${\bf U}$ は変位ベクトル、${\bf F}$ は荷重ベクトルである。$x$ は設計変数であり、剛性行列は設計変数に依存する。
平均コンプライアンスは、
\begin{align}
c = f(x)= {\bf U}^T {\bf K}(x) {\bf U}
= \sum_{e=1}^{N_e} {\bf u}_e^T {\bf k}_e(x_e) {\bf u}_e
\end{align}
と定義され、平均コンプライアンスは変位の２乗に比例する。$x_e$ は要素 $e$ における設計変数であり、$u_e,K_e$ は要素 $e$ における変位と剛性行列を表す。$N_e$ を要素数とする。平均コンプライアンスを小さくすることは、変位を小さくする構造を探すことを意味し、これは剛性が最大になるように最適化することと同じである。（変位が小さいならば硬い構造と言える。）平均コンプライアンスを目的関数として、最小化を行う。
最適化を行う際は、不等式制約として
\begin{align}
g(x) = V -\beta V^0 = \sum_{e=1}^{N_e} x_e v_e -\beta\sum_{e=1}^{N_e} v_e \leq 0
\end{align}
を課す。$v_e$ は要素における体積である。$\beta$ は体積の減率を表し、例えば、$\beta=0.5$ の場合、計算初期の体積$V^0$を$50$%以下に削減するよう最適化を行う。
各要素の剛性行列に含まれるヤング率$E_e$は、設計変数$x_e$と物理的なヤング率$E^0$を使い、
\begin{align}
E_e(x_e) = x_e^p E^0
\end{align}
と表現する。$p$ はペナルティ指数であり$p>1$ である。この方法におけるトポロジー最適化をSIMPと呼ぶ。設計変数の範囲は、
\begin{align}
0< x_e^L \leq x_e \leq x_e^U=1
\end{align}
である。（$x_e = 0$ だと剛性行列が特異になるため、下限値を$x_e^L=10^{-3}$に設定する）
ラグランジュ関数と停留条件
上記の最適化問題をラグランジュ関数で表すと
\begin{align}
L &= f(x) + \Lambda g(x) + \sum_{e=1}^{N_e}\lambda_e (x_e^L-x_e)
+ \sum_{e=1}^{N_e}\gamma_e (x_e-x_e^U) \\
&= {\bf U}^T {\bf K}{\bf U} + \Lambda\left( \sum_{e=1}^{N_e} (x_e-\beta) v_e \right) + \sum_{e=1}^{N_e}\lambda_e (x_e^L-x_e)
+ \sum_{e=1}^{N_e}\gamma_e (x_e-x_e^U) \\
& \Lambda \geq 0 \ , \ \lambda_e \geq 0 \ , \ \gamma_e \geq 0
\end{align}
となる。$\Lambda, \lambda_e ,\gamma_e$ はラグランジュ未定乗数である。ラグランジュ関数を設計変数で微分し停留条件を求めると、
\begin{align}
\frac{\partial L}{\partial x_e} &= \frac{\partial {\bf U}^T}{\partial x_e} {\bf K} {\bf U} + {\bf U}^T \frac{\partial {\bf K}}{\partial x_e} {\bf U} +
{\bf U}^T {\bf K} \frac{\partial {\bf U}}{\partial x_e} +
\Lambda v_e -\lambda_e + \gamma_e = 0 \\
\end{align}
となる。荷重は設計変数に依らないとすると
\begin{align}
&\frac{\partial }{\partial x_e}\left({\bf K}{\bf U}-{\bf F} \right)=0 \\
&\therefore {\bf K}\frac{\partial {\bf U}}{\partial x_e} = -\frac{\partial {\bf K} }{\partial x_e}{\bf U}
\end{align}
と書ける。剛性行列は対称行列とすれば停留条件は、
\begin{align}
\frac{\partial L}{\partial x_e} &= -{\bf U}^T \frac{\partial {\bf K}}{\partial x_e} {\bf U} +
\Lambda v_e -\lambda_e + \gamma_e = 0 \\
\end{align}
と書ける。また、剛性行列の設計変数の微分は、
\begin{align}
\frac{\partial f(x)}{\partial x_e} = -{\bf U}^T \frac{\partial {\bf K}}{\partial x_e} {\bf U} = -px_e^{p-1} {\bf u}^T_e {\bf k}_e {\bf u}_e
\end{align}
である。
$\lambda_e,\gamma_e$ に関しては、$x_e^L-x_e < 0 $ または $x_e-x_e^U < 0$ のとき非活性とよばれ（等式制約が成立していない）、$\lambda_e=0 ,\gamma_e=0 $ となる。
フィルタリング
トポロジー最適化は、チェッカーボードの問題が生じる。これは、設計変数を要素に設定すると、隣り合う要素同士では独立した値を取るため、計算結果がチェッカーボードのようになる。
チェッカーボード問題を克服する有名な方法としてフィルタリングがある。フィルタリングは、目的関数の要素 $e$ における設計変数の微分を
\begin{align}
\frac{\partial \tilde{f}(x)}{\partial x_e} = \frac{1}{x_e \sum_{i=1}^{N_e}H_{ei}} \sum_{i=1}^{N_e}H_{ei} x_i \frac{\partial f(x)}{\partial x_i} \\
H_{ei} = \begin{cases}
r_{\min} -\Delta_{ei} & \Delta_{ei} \leq r_{\min} \\
0 & r_{\min} < \Delta_{ei}
\end{cases}
\end{align}
と書き換える。$\Delta_{ei}$ は要素 $e$ と $i$ の中心からの距離である。
（重み付き平均の形にすることで隣り合う要素の値に相関性を持たせることができる？）
設計変数の更新
最適基準法を用いる場合、
\begin{align}
G_e \equiv -\frac{\frac{\partial \tilde{f}(x)}{\partial x_e}}{ \Lambda v_e }=1
\end{align}
となるように $\Lambda$ を決め $x_e$ を更新する。$x_e$ の更新は、$k$ を繰り返しステップ数とすると
\begin{align}
x_e^{(k+1)} =
\begin{cases}
M_x^{-} & \left[G_e^{(k)} \right]^{0.5} x_e^{(k)} \leq M_x^{-} \\
\left[G_e^{(k)} \right]^{0.5}x_e^{(k)} &
M_x^{-} < \left[G_e^{(k)} \right]^{0.5}x_e^{(k)} < M_x^{+} \\
M_x^{+} & M_x^{+}\leq \left[G_e^{(k)} \right]^{0.5} x_e^{(k)} \\
\end{cases}
\end{align}
\begin{align}
M_x^{-} &= \max \{(1-\mu_x)x_e^{(k)},x_e^L \} \\
M_x^{+} &= \max \{x_e^U,(1+\mu_x)x_e^{(k)} \} \\
\end{align}
となるように更新する。$\mu_x$ は move limit と呼ばれ $\mu_x =0.2$ に設定する。
収束条件としては、
\begin{align}
\max \left| x_e^{(k+1)} - x_e^{(k)} \right| \leq 10^{-2}
\end{align}
を課し、設計変数が更新しなくなったら計算を終了する。
計算結果
以下の図のように、幅 $6$m 、高さ $1$m のコンクリートに上端中央部から集中荷重 $200$ kN の荷重をかける。
ヤング率は $28$Gpa 、ポアソン比は $0.2$ とし平面応力とした。体積の減率は $\beta = 0.2$ とした。フィルタリングで用いるパラメータ $r_{\min}$ は、最大となる要素辺の長さの $1.5$ 倍とした。
以下がトポロジー最適化の結果である。参考書と類似した結果が得られた。
まとめ
今回の記事では、トポロジー最適化を用いてコンクリートの最適設計を行った。設計変数を要素ごとに定義するとチェッカーボード問題が生じるが、フィルタリングをすることで、この問題を避けた。
コンクリートは、引張に弱いので引張応力がかかるところには steel を挿入する。次回の記事では、最適な箇所に steel を挿入するトポロジー最適化について紹介したい。