流れ場のトポロジー最適化に関する研究
Study on the topology optimization of a flow field
精密工学専攻 12号 加藤 英祐 Eisuke Kato
1. はじめに
自動車が水の溜まった路面を走行するとき，ハイドロプ レーニング現象が生じることがある．これは，Fig.1のよ うに路面とタイヤとの隙間に水が入り，車が薄い水の膜の 上を滑るようになり，ハンドル操作やブレーキが利かなく なる現象である．ハイドロプレーニング現象はスピードの 出しすぎや，タイヤの磨耗など，路面に溜まった水の量が タイヤの排水能力を超えた場合に発生する．各タイヤメー カではハイドロプレーニング現象に関して，Fig.2のよう な実験や数値シミュレーションが行なわれている．Fig.2 を見ると，路面とタイヤの間に入り込んだ水がタイヤの 溝パターン（トレッドパターン）に沿ってを排出されてお り，トレッドパターンは排水能力に大きく影響していこと がわかる．排水能力の高いトレッドパターンを最適設計に よって
最適設計問題は，抗力を最小にする翼型設計などの形状 最適化と，本研究のようなトポロジー最適化とに大きく 二分される．形状最適化では，最初に設定した翼型物体の 境界が変形することのみを考え，翼型が二つに分かれると いった境界の増減は行えない．一方，トポロジー最適化は，
設定した領域全体を最適設計領域とし，物体・流体の有無 を一点で定義するため，境界の増減が可能となる．
本研究では，タイヤの接地面に見立てた計算領域を流れ る流体の最適化計算を行う．実際には，タイヤが回転する ことにより，接地面のトレッドパターンも時々刻々変化し ていくが，簡易化のために計算領域を固定し，流体を非圧 縮粘性流体の遅い定常流れとする．流れの解析には有限要 素法を用い，最適化手法には逐次凸関数近似法のMMA(3) を用いて近似化した最適流路設計問題を形成し，最急降下 法によって最適解を求める．
2. 流れ解析のための定式化 2.1 流れの支配方程式
本研究では非圧縮粘性流体を扱い，Fig.3のような厚さ 2ρの領域R を考える．流れの支配方程式は次の定常ス トークス方程式と，連続の方程式となる．
µ∇2u¯ =∇¯p−f¯ (1)
∇ ·u¯ = 0 (2)
ここで，µは流体の粘性係数，p¯は圧力，u¯ は流速ベクト ル，f¯ は物体力ベクトルである．
境界∂R全体で流速が指定されており，
¯
u= ¯g (3)
で表される．また，連続性より，
Z
∂R
¯
g·n= 0 (4)
が満たされる．ここで，nは外向き法線単位ベクトルで ある．
Fig. 1 Hydroplaning phenomenon
Fig. 2 Experiment of Hydroplaning phenomenon(8)
ここで，領域Rのx1，x2方向のスケールに比べてx3
方向の厚さ2ρは微小であり，x3方向の流体運動や圧力の 変化は無視できると仮定し，次式で定義する平均流速を用 いて式(1)，(2)をx3方向に平均化する．
u= 1 2ρ
Z ρ
−ρ
¯
udx3 (5)
また，ρを正の定数とする関数 ξ(x) =ξ(x1, x2, x3) =3
2 µ
1−x23 ρ2
¶
(6)
を用い，x1，x2方向の流速の分布を
u¯ =ξu (7)
と仮定することにより，平均化されたストークス方程式と 連続の方程式
µ∇2u+α(ρ)u=∇p−f (8)
∇ ·u= 0 (9)
となる．ここで，式(8)の左辺第2項は式(1)の左辺を平 均化したときに生じる項であり，
α(ρ) = 3µ
2ρ2 (10)
である．これにより，3次元流れの問題はFig.4に示され るような2次元流れ問題に帰着される．
Fig. 3 Illustration of a plane flow model
Fig. 4 Analytic area
また，このときの境界∂Ωでの条件として，
u=g (11)
が与えられる．
式(8)と多孔質体における一般化されたDarcyの法則 u=−κ
µ(∇p−f) (12)
とを比較すると，式(8)における粘性µ= 0としたとき，
α−1はダルシーの法則(12)で見られる透水係数κに相当 するパラメータとなる．透水係数は物質の流れ易さを表す 係数であり，値が大きいほど流れ易い事を意味する．
2.2 支配方程式の離散化
平均化されたストークス方程式(8)と連続の方程式(9) の弱形式はガラーキン法を用いて
Z
Ω
α(ρ)u·u∗+µ Z
Ω
∇u:∇u∗
− Z
Ω
p∇ ·u∗ = Z
Ω
f·u∗ (13)
− Z
Ω
p∗ ∇ ·u= 0 (14)
となる．ここで，u∗とp∗は重み関数である．弱形式(13)， (14)の空間方向の離散化には有限要素法を用いる．Fig.5 のように三角形を四つの小三角形に分割し，流速uを各小 三角形で1次近似，圧力pは大三角形で1次近似とする．
また，ρは小要素内一定と近似する．
次に，弱形式(13)，(14)は離散化すると，
· D −H HT 0
¸ ½ U P
¾
=
½ F 0
¾
(15)
(a) Velocity (b) Pressure Fig. 5 Finite elements
が得られる．ここで，Uは流速ベクトル，Pは圧力ベクト ル，Fは物体力ベクトルを表し，D，Hはそれぞれ粘性項，
圧力項に対応する係数行列である．
本研究ではLU分解法を用いて連立一次方程式(15)か ら流速U，圧力Pを求める．
3. 最適流路形状決定問題の定式化
トポロジー最適化問題は，目的関数に計算領域全体のポ テンシャルエネルギを用い，最小ポテンシャルエネルギ 問題と考えることができる．また，設計変数にρを用い，
ρ= 0でα=∞となり，固体領域を示し，ρ= 1でα= 0 となり，流体領域を示す．ρの分布により計算領域の流路 を決定する．最適化の手法はK Svanberg(3) が提案した MMAを用いる．
3.1 最適化問題の定式化
本研究における流体のポテンシャルエネルギを
φ(ρ) =1 2
Z
Ω
α(ρ)u·u +µ
Z
Ω
∇u:∇u− Z
Ω
f ·u (16)
と定義する．これを離散化すると，
Φ(ρ) =1
2UTDU−F˜TU˜ (17) が得られる．ここで，F˜ とU˜ はそれぞれ領域Ω内のみの 節点における外力ベクトルと流速ベクトルを表し，∂Ωの 節点値は含まれない．この目的関数は，第1項で消散エネ ルギを最小にし，第2項で外力の生じている場所で流速を 最大にすることを意味している．消散エネルギを最小にす ることは，抗力や領域内の圧力低下を最小にすることに等 しい．また，式(17)を式(15)で得られる流速，圧力を用 いて書き直すと，
Φ(ρ) =1 2
©UT(D˜g−F) +PTB˜gª
(18)
となる．ここで，g˜は境界節点の流速値を成分に持ち，Ω内 の節点に対応する成分は0で構成されるベクトルである．
次に，設計変数に関する制約条件は，
0≤ρ≤1 (19)
Z
Ω
ρdΩ≤γ|Ω| (20)
となる．式(19)はρ= 0で固体領域を表し，ρ= 1で流 体領域を表すことから，設計変数の下限上限条件を意味す
る．また，式(20)は面積制約条件であり，固定設計領域Ω 中を流体領域が占める割合の上限を意味し，γは0から1 の間の一定値である．本研究ではγ= 0.5を用いる．
これらを離散化し，まとめると，求める最適流路決定問 題Pは
P :
minρ Φ (ρ)
s.t. eTρ≤γ|Ω|
0≤ρ≤1
(21)
と表すことができる．ここで，eは小三角形要素の面積ベ クトルである．設計変数ρは小三角形要素内で一定と近似 したベクトルρで表す．
最終的な最適解はρが0（固体領域）か1（流体領域）
をとることを目的とするので，ρが中途半端な値をとるこ とのないようにαを選ぶ必要がある．また，ρ = 0のと きα =∞となってしまうエラーをなくすためのペナル ティーパラメータqを導入し，α(ρ)に次の補間関数を用 いることにする．
αq(ρ) =α+ (α−α)ρ1 +q
ρ+q (22)
この補間関数を用いることで，ρ= 0のとき，α(ρ)は∞で はなく，有限に大きな値をとることができる．また，ρ= 1 のとき，α(ρ)→0 の値をとることができるため，一般的 な2次元のストークス方程式
µ∇2u=∇p−f (23)
を保持することができる．本研究では，q= 0.1，µ= 1.0， α= 1.5µ/1002，α= 1.5µ/0.012を用いる.
3.2 MMA展開
逐次凸関数近似法の一つであるMMA（Method of Mov- ing Asymptotes）は構造最適化問題の手法としてSvan-
berg K(3)によって考案された．設計変数に媒介変数を導
入し，目的関数と制約条件を媒介変数に基づき1次のテイ ラー展開により線形近似する．この媒介変数を感度（目的 関数の勾配）の正負によって決定することにより，目的関 数と制約関数は凸関数に近似される．近似された目的関数 と制約条件から成る新たな最適化問題をMMAの反復計 算の中で毎ステップ形成することで，簡易化された問題を 解くことができる．
MMA展開の反復回数kにおける近似目的関数は Φ(k)(ρ) =r(k)+
Xm
i=1
qi(k)
ρi−l(k)i (24)
となる．ここで，mは小三角形要素数，ρiはi番目のρの 要素値を表す．パラメータli(k)，ρ(k)そしてρ(k)はMMA での繰り返しの間で更新される．また，パラメータr(k)と q(k)i は，
Φ(k)
³ ρ(k)
´
= Φ
³ ρ(k)
´
(25)
∂Φ(k)¡ ρ(k)¢
∂ρ(k)i = ∂Φ¡ ρ(k)¢
∂ρ(k)i (26)
を満たすように選ばれる．
また，MMAのパラメータ決定に必要な目的関数の勾 配は，
∂Φ (ρ)
∂ρi = 1
2α0(ρi)U(ρ)TDi(α)U(ρ) (27) となる．さらに，α0(ρi)≤0，U(ρ)TDi(α)U(ρ)≥0が 成り立ち，勾配は0または負となることがわかる．つまり，
目的関数は最適値に向かって改良される方向に向かう結果 になる．ここで，Di(α)はαに依存する要素マトリック スである．
制約条件も同様に近似することで，離散化された最適化 問題(21)はMMA展開を用いて，反復回数kとそのとき の設計変数ρ(k)における流速ベクトルU¡
ρ(k)¢
と圧力ベ クトルP¡
ρ(k)¢
から，次の近似化した最適流路設計問題 P(k)を形成する．
P(k):
minρ Φ(k)(ρ)
s.t. eTρ≤γ|Ω|
ρ(k)≤ρ≤ρ(k)
(28)
ここで，MMAの反復計算ついて簡単にその計算手順を 説明する．
Step1. 始点ρ(0)を与え，k= 0とする．
Step2. 反復点ρ(k)を与え，Φ¡ ρ(k)¢
と勾配∇Φ¡ ρ(k)¢ を計算する．
Step3. Step1.からの計算によって目的関数，制約関数を
近似することで，最適化問題Pを近似化した最適化 問題P(k)に置き換える．
Step4. P(k)を計算し，この最適解を次の反復点ρ(k+1) にする．k:=k+ 1としてStep1.に行く．
このように最適化問題P から反復kにおける近似した 最適化問題P(k)を得ることで，問題を容易にしてから解 くことができる．本研究では式(27)より勾配を用いるこ とができるため，最適化手法には勾配法を用いることにす る．勾配法の中でもアルゴリズムが簡単で大域収束性のあ る最急降下法を用いることとする．
以上，全体の計算の流れをまとめると，前章で述べた流 れ計算により求まる流速Uと圧力Pから全ポテンシャル
Fig. 6 The flow chart of the total calculation
エネルギΦと勾配∂Φ/∂ρを得る．MMAの反復計算に よって，近似された最適化問題P(k)の最適解ρ を更新し ていき，最適流路形状決定問題P の最適解を求める．
一連の計算の流れをFig.6 に示す．ここで，Fig.6にお ける収束条件は2回続けて目的関数の相互誤差が10−5 よ り小さくなることとする．
4. 数値計算例
計算モデルをFig.7に示す．計算領域は一辺の長さが1 の正方形領域とし，左を流入境界，右の中央1/3を流出境 界，上下を含んだその他の壁は全てno-slip境界条件とす る．流入境界と流出境界は共にディリクレ型境界条件とし て，流入中央を最大流入速度1，流出中央を最大流出速度 3としてそれぞれ分布をFig.7のような放物型とする．ま た，圧力の基準点として，流出境界に圧力p= 0を与え，
計算領域に働く外力f = 0とする．トポロジー最適化プロ グラムを用いてこの計算モデルで流入から流出までの流体 領域と固体領域の分布を明確にすることで最適流路形状を 求める．
計算領域を大三角要素数400，節点数221のメッシュ
Fig. 7 Design domain for the diffuser model
(a) 400 element mesh (b) 1600 element mesh Fig. 8 Optimal design of the diffuser model
Fig. 9 Gray zone of Fig.8(a)
に切ったときと，大三角要素数1600，節点数841のメッ シュに切ったときの最適流路形状はそれぞれFig.8の(a) と(b)に示す結果が得られた．白がρ= 1で流体領域を示 し，黒がρ = 0で固体領域を表している．どちらも流入 部から流路が狭まり，ほぼ一定の流路幅をもって流出部分 に向かっている傾向が得られた．また，要素数の違いによ り，流体領域と固体領域の境界の滑らかさが異なり，要素 数を多くすることで得られる形状はより滑らかになること がわかる．
また，Fig.8(a)の左下流入部分を拡大したFig.9を見る とわかるように，流入部分に若干グレーの領域が見られ る．この領域は流体領域と固体領域の中間の値を示してお り，トポロジー最適化においては好ましくない．このよう な領域が存在することで，境界をどこにするのが最適であ るかという設計上の問題点になる．Fig.9に示した流入部 分でのグレー領域は，境界流速が遅いために現れたと考え ることができる．このような境界部分において，より細か いメッシュを切る必要があるが，全体として白か黒のどち らかを示しており，境界を正確に捉えた最適形状を決定す る上で良好な結果が得られたといえる．
5. おわりに
本研究では，厚さに相当するパラメータρを設計変数と して，ポテンシャルエネルギΦを目的関数とした2次元定 常ストークス流れにおける最適流路形状決定のための，流 体トポロジー最適化プログラムを構築した．グレーの領域 がほとんどなく，流体領域と固体領域の境界が明確に判断 できる流路形状を得ることができた．
流体領域と固体領域とで，流れの支配方程式を同じもの として扱ったが，流体領域(ρ= 1)ではストークス方程式，
固体領域(ρ= 0)ではDarcyの法則といった，支配方程式 を可変にすることで，より正確なトポロジー最適化が行え ると考える．さらにはトレッドパターンのような複雑な流 入・流出境界条件での計算を行うなどによって発展が見込 まれる．
参考文献
(1) T. Borrvall and J. Petersson, Topology optimiza- tion of fluids in stokes flow. International Journal for Numerical Methods in Fluids, 2003;41:77-107.
(2) M. P. Bendsøe and N. Kikuchi, Generating opti- mal topologies in structural design using a homog- enization method. Computer Methods In Applied Mechanics And Engineering, 1998;71:197-224.
(3) K Svanberg, The method of moving asymptotes – a new method for structural optimization.Interna- tional Journal for Numerical Methods in Engineer- ing, 1987; 24:359-373.
(4) J Cahouet, JP Chabard, Some fast 3D solvers for the generalized Stokes problem.International Jour- nal for Numerical Methods in Fluids, 1988; 8:869- 895.
(5) James K. Guest and Jean H. Pr´evost, Topology optimization of creeping fluid flows using a Darcy- Stokes finite element.International Journal for Nu- merical Methods in Engineering, 2006;66:461-484.
(6) 今野 浩,山下 浩,非線形計画法 日科技連出版社, 1978.
(7) 久志 本茂,最適化問題の基礎 森北出版, 1979.
(8) http://www.yokohamatire.jp/check-de- smile/sp characteristic/index.html
