2 最適形状決定問題

修士論文発表要旨(2009年度)

随伴法を用いた非圧縮性流体中の物体の最適形状決定

Shape Optimization of a Body Located in Incompressible Flow Using Adjoint Method

土木工学専攻  ２２号 坂本 雅人 Masato SAKAMOTO

1 はじめに

流体中に置かれた物体が受ける抗力の大きさは物体形 状に依存する.以前はある条件を満たす形状の設計は実験 や技術者の経験に基づいて行なわれきた.しかし実験や経 験で形状を決定した場合,理論的裏づけがなく,費用や時 間がかかるといった問題がある.そこで本研究では,最適 な形状を最適制御理論に基づいて数値的に決定する.抗力 最小化などの最適形状決定問題の取り組みは, Pironneau 1)らによってStokes近似の成立する範囲でその解法のア ルゴリズムと解の性質などが示され,現在では数値流体力 学の発展によって流れ場が精度良く解析できる事によっ て本格化している.汎用性のある解法として抗力や揚力を 直接目的関数に取り,物体表面の座標値を設計変数に使用 して形状決定を行なう取り組みがされており, 2),3),4),5) の中で定常,非定常流れでの形状同定に成功している.し かし,ここで得られている解については初期形状に依存し ている可能性がある.そこで本研究では,全体制御アルゴ リズムを用いて異なる初期形状より導出された解につい て比較検証し,最終形状と初期形状の依存を検証する.ま た物体を更新する勾配が流体力の影響を受けるため,全体 制御アルゴリズムを用いた場合,物体後方部分の勾配が前 方に比較し小さいという結果を得た. この問題により物 体後方の形状勾配は形状更新に反映され難くなる. そこ で本研究で提案する部分制御アルゴリズムを用いて後方 部分を独立に更新にさせる事によって,初期形状に依存し ない最終形状を導出する.このように本研究の目的は,最 適形状問題における初期形状の依存の検証と部分制御ア ルゴリズムを用いた最適形状の導出である.

2 最適形状決定問題

2.1

境界Γを有する２空間領域Ω ∈R²を定義し,領域Ω は非圧縮性流体で満たされているものとする. 支配方程 式には,以下の無次元化された非圧縮性Navier-Stokes方 程式を用いる.

˙

u_i+u_ju_i,j+p_,i−ν(u_i,j+u_j,i)_,j = 0 in Ω(1) u_i,i = 0 in Ω(2)

ここでu_iとpは流速と圧力であり,νは動粘性係数である.

2.2

本研究では直接に流体力を評価関数を用いる. 本研究 では抗力と揚力を最小とするような物体形状を解とする ため,評価関数は支配方程式より計算された流体力の二乗 和をもって定義される.

J = _tf

t0

1

2(q₁F₁²+q₂F₂²)dt, (3) ここで,q1とq2はそれぞれ抗力と揚力に対する重み関数 を表す.

2.3

支配方程式の補間関数には安定化気泡関数有限要素法 を用いる.この手法では,流速に関しては気泡関数要素を, 圧力に関しては一次要素を補間関数に用いる.

u_i= Φ_αu_αi= Φ₁u_1i+ Φ₂u_2i+ Φ₃u_3i+ Φ₄u˜_4i, (4)

˜

u_4i=u_4i−1

3(u_1i+u_2i+u_3i), (5) Φ₁=l₁, Φ₂=l₂, Φ₃=l₃, Φ₄= 27l₁l₂l₃, (6) p= Ψ_λp_λ= Ψ₁p₁+ Ψ₂p₂+ Ψ₃p₃, (7) Ψ₁=l₁, Ψ₂=l₂, Ψ₃=l₃, (8)

2

3 1

4

図1:気泡関数要素

2

3 1

図2:線形要素

2.4

面積一定条件を導入する事により,物体形状は全ての 計算過程において面積が一定となるように最適化される.

物体面積が一定である事は,有限要素メッシュの全要素の 面積の合計が一定である事と等価である. よって面積一 定条件は次の方程式で示される.

m

e=1

(a_e(X_i))^(l)−A₀= 0, (9)

ここで,a_e(x_i)は各要素の面積を,A₀は初期領域の面積 を表す.

2.5

本研究ではLagrange乗数法を用いて,拘束条件付きの 最小化問題を拘束条件なしの最小化問題に置き換える.

Lagrange乗数 を導入すると,拡張評価関数は式(6)で表

される.

J^∗ = _tf

t0

1

2(q₁F₁²+q₂F₂²)dt

− _tf

t0

Ωu^∗_i{u˙_i+u_ju_i,j+p_,i−ν(u_i,j+u_j,i)_,j}dΩdt +

_tf

t0

Ωp^∗u_i,idΩdt.

(10) そして、この拡張評価関数の第一変分を取ると、式(11) のようになる.

δJ^∗ = − _tf

t0

Ωδu^∗_i{u˙_i+u_ju_i,j+p_,i

− ν(u_i,j+u_j,i)_j}dΩdt +

_tf

t0

Ωδp^∗u_i,idΩdt

− _tf

t0

Ωδu_i{−u˙^∗_i −u_ju^∗_i,j+u_j,iu^∗_j+p^∗_,i

− ν(u^∗_i,j+u^∗_j,i)_,j}dΩdt +

_tf

t0

Ωδpu^∗_i,idΩdt

− _tf

t0

ΓD

δu_is_idΓdt− _tf

t0

ΓS

δu₁s₁dΓdt

− _tf

t0

ΓB

δu_is_idΓdt +

_tf

t0

ΓU

δt_iu^∗_idΓdt+ _tf

t0

ΓS

δt₂u^∗₂dΓdt +

_tf

t0

ΓB

δt₁(u^∗₁−q₁F₁)dΓdt +

_tf

t0

ΓB

δt₂(u^∗₂−q₂F₂)dΓdt,

−

Ωu_i(t_f)^∗δu_i(t_f)dΩ +

Ωu_i(t₀)^∗δu_i(t₀)dΩ, (11) 評価関数を最小とするための最適条件は,拡張評価関 数の第一変分をゼロとする事である. その条件を考慮す る事により,支配方程式,随伴方程式(12),(13),境界条件 (14),(15),(16),(17),終端条件(18),勾配(19)を得る事がで

きる.

−u˙^∗_i −u_ju^∗_i,j+u_j,iu^∗_j + p^∗_,i

−ν(u^∗_i,j+u^∗_j,i)_,j = 0 in Ω, (12) u^∗_i,i = 0 in Ω (13) u^∗_i = 0 on Γ_U(14) u^∗₁−q₁F₁= 0, u^∗₂−q₂F₂ = 0 on Γ_B(15) s₁= 0, u^∗₂ = 0 on Γ_S(16) s_i = 0 on Γ_D(17) u^∗_i(t_f) = 0 in Ω (18) grad(J^∗) = −s_ju_j,i (19) u^∗, p^∗を各条件を満足するように解き,式(19)より形状更 新に用いる勾配を計算する.

2.6

本研究では最小化手法として重み付き勾配法を用いる.

形状決定のパラメータは以下の式によって更新される.

X_i^(l+1)W_ij^(l)=X_i^(l)W_ij^(l)−grad(J^∗)^(l)_i , (20)

また, 形状更新には物体周りの節点全てを更新する全 体制御アルゴリズムと部分的に更新する部分制御アルゴ リズムの２点を適用する.

2.7

本研究では全体制御アルゴリズム、部分制御アルゴリ ズムの２つのアルゴリズムを用いて形状決定を行なう.部 分制御アルゴリズムでは勾配の小さくなる物体後方の勾 配を形状更新に反映させるため物体後方半分を部分制御 境界とする. 部分制御アルゴリズムでは部分的に制御し ただけでは制御境界で形状のギャップが発生してしまう ため、全体制御と部分制御を交互に行なう.

全体制御アルゴリズム 1.初期の物体座標の決定.

2.基礎方程式の状態量u⁽⁰⁾_i ,p⁽⁰⁾を解く.

3.随伴方程式の随伴変数u^∗(0)_i ,p^∗(0)を解く.

4.物体全体の座標をX_i^(l)式(8)によって計算する.

5.形状を更新し,計算領域のリメッシュをする.

6.基礎方程式の状態量u⁽⁰⁾_i ,p⁽⁰⁾を解く.

7.もし|J^(l+1)−J^(l)|< ならば計算を終了し, そうでなければ3へ.

部分制御アルゴリズム 1.初期の物体座標の決定.

2.基礎方程式の状態量u⁽⁰⁾_i ,p⁽⁰⁾を解く.

3.随伴方程式の随伴変数u^∗(0)_i ,p^∗(0)を解く.

4.物体表面の座標をX_i^(l)式(8)によって計算する.

5.形状全体を更新し,計算領域のリメッシュをする.

6.基礎方程式の状態量u⁽⁰⁾_i ,p⁽⁰⁾を解く.

7.もし|J^(l+1)−J^(l)|< ₁ならば8へ, そうでなければ形状を更新前に戻して8へ.

8.随伴方程式の随伴変数u^∗(0)_i ,p^∗(0)を解く.

9.部分制御部分の座標をX_i^(l)式(8)によって計算する.

10.形状を更新し,計算領域のリメッシュをする.

11.基礎方程式の状態量u⁽⁰⁾_i ,p⁽⁰⁾を解く.

12.もし|J^(l+1)−J^(l)|< ₂ならば13へ, そうでなければ形状を更新前へ戻して13へ.

13.全体,部分制御共に座標更新無しならば計算終了.

そうでなければ3へ.

3 数値解析結果

 数値解析例として,図３,４のような領域における物 体の流体力最小化問題を行なった. 全体制御を施す初期 形状として３つの初期形状を設定し,初期形状への依存を 比較検討する. Case 1は円, Case 2はRe=100での最終形 状, Case 3はNACA 4 Digits SeriesよりCase 2と同様な 幅厚比の形状とし,それぞれの初期形状より最終形状を求 める. 部分制御を施す初期形状をCase 4として,幅厚比 が1 : 7となるような楕円を設定し,全体制御を施した形 状と比較検討する. またレイノルズ数は250.0と設定し た.数値解析結果として,各ケースの初期形状と最終形状 を図５〜１２,流線図を図１３〜１６,評価関数,抗力,揚 力の時系列を図１７,１９,２０に示す.

円形状時の評価関数値を1.0と基準化すると, Case 1,2,3,4 はそれぞれ71.5%,77.5%,78.8%,78.8%減少させる事がで きた.それぞれの形状を比較すると, Case 1においては後 方部分が丸みを帯びており,初期形状である円に依存して いる. 物体後方部ではカルマン渦が発生し、物体の揚力 変動が大きい事がわかる. Case 2においては物体後方部 分の勾配が影響が小さく, Case 3,4に比べ物体後方で最大 幅が減少せず抗力が大きくなっている. Case 3では物体 後方部での渦が発生しない初期形状を選択したため、最 終形状においても物体後方部では渦が発生しない形状と なった. Case 4では部分制御アルゴリズムを適用したた め、物体前方だけでなく後方部も鋭くなり、渦、揚力共 にの発生しない安定した形状となった. 物体また最大厚 さの部分がCase 1,2では後方, Case 3では前方に配置し

ており,この点からも初期形状への依存がわかる. Case 4 においては初期形状である楕円に依存する事なく,後方部 分が鋭くなっている事がわかる. 評価関数はCase 3,4の 最終形状での値が同値となり、形状自体もほぼ同様の形 状となった.

4 結論

 最適制御理論に基づく形状決定のアルゴリズムを適 用しての定式化,解析を行なう事に成功した.また異なる 初期形状より評価関数を減少させるような最終形状を導 出する事に成功したが,最終形状は異なる形状となった.

今回の研究結果より,最終形状は初期形状に依存する事が わかり、この原因としては物体後方部分の勾配が前方よ り小さいため、物体後方の勾配が形状最適化に反映され ない事が挙げられる. 部分制御アルゴリズムを適用した ケースでは初期形状の楕円は両端が曲線状になっている にも関わらず,最終形状では両端が鋭い形状となり,初期 形状に依存しない最終形状を導出する事ができた. この 結果より最適形状問題における部分制御アルゴリズムの 有効性が証明された.

参考文献

[1] Pironneau,O. Consist Approximations, Automatic Dif- ferentiation and Domain Decomposition For Optimal Shape Design, International Series, Mathematical Sci- ence and Applications Vol. 16, Computational Methods for Control Application, pp. 167-178, 2001.

[2] Jameson,A. Aerodynamic Shape Optimization Using the Adjoint Method, Lectures at the Von Karman In- stitute, Brusseles, 2003.

[3] Nakajima,S. and M.Kawahara, Shape Optimization of a Body in Compressible Inviscid Flows, Comp. Meth.

Appl. Mech. Engrg., Vol. 197, pp. 4521-4530, 2008.

[4] Yagi,H. and M.kawahara, Optimal Shape Determina- tion of Body Located in Incompressible Viscous Fluid Flow, Comp. Meth. Appl. Mech. Engng, Vol. 196, pp.

5084-5091, 2007.

[5] Yoshida,H. and M.Kawahara, Shape Optimization of an Oscillating Body in Fluid Flow by Adjoint Equa- tion And ALE Finite Element Methods, Int. J. Numer.

Meth. Fluids, Vol. 22, pp. 229-239, 2006.

D

U = 0 V = 0 V = 0

V = 0

12.5D 12.5D

(-14.0, -4.5)

4.0D Ω

U = 1,V = 0

(14.0, -4.5) (14.0, 4.5) (-14.0, 4.5 )

A

B

C

D

図3:解析領域と境界条件

図4:有限要素メッシュ

図5: Case 1:初期形状 図6: Case 1:最終形状

図7: Case 2:初期形状 図8: Case 2:最終形状

図9: Case 3:初期形状

図10: Case 3:最終形状

図11: Case 4:初期形状 図12: Case 4:最終形状

X

Y

2 2.2 2.4

-0.2 0 0.2

図13: Case 1:後方の流線

X

Y

2 2.2 2.4

-0.2 0 0.2

図14: Case 2:後方の流線

X

Y

2 2.2

-0.2 0 0.2

図15: Case 3:後方の流線

X

Y

2 2.2 2.4

-0.2 0 0.2

図16: Case 4:後方の流線

Iteration

Performancevalue

0 100 200 300 400 500

0.2 0.4 0.6 0.8 1

Case 1 Case 2 Case 3 Case 4

図17:評価関数の時刻暦

X

Y

-2 -1 0 1 2

-1 -0.5 0 0.5

1 Circle

Case 1 Case 2 Case 3 Case 4

図18:最終形状の比較

Step

DragForce

500 1000 1500 2000

0.4 0.6 0.8 1 1.2 1.4 1.6

Circle Case 1 Case 2 Case 3 Case 4

図19:抗力の時刻暦

Step

LiftForce

0 500 1000 1500 2000

-1 -0.5 0 0.5 1

1.5 Circle

Case 1 Case 2 Case 3 Case 4

図20:揚力の時刻暦