— 21 —SENAC Vol. 49, No. 4（2016. 10）［研究成果］— 情報科学研究科グループ利用制度を利用して —Shape Optimization による非定常流れ場の安定性制御
Shape Optimization Problem Improving Hydrodynamics Stabilityon the Non Stationary Navier-Stokes Problem
中澤嵩
東北大学大学院情報科学研究科数学連携推進室
本研究では，運動エネルギー最小化問題を行い，初期形状と最適形状における線形安定性を比較する．その際，初期形状を 2 次元 Cavity 流れとし，レイノルズ数Re = 11500において非定常
Navier-Stokes 問題を主問題とする．数値計算の結果，初期形状と比較して最適形状における臨界レイノルズ数Re� が減少したことが数値的に示された．運動エネルギーを最小化しているにも拘らずRe� が減少するというこの結果は，流体物理の観点では非常に興味深い．今後，更なる詳細な検証が求められる．

1. はじめに形状最適化問題は，偏微分方程式が定義された空間における境界の形状を最適化することである．初期には，1903 年に J. Hadamard によって薄膜の振動数最大化問題として行われた．その後，特にフランスの応用数学者によって線形弾性体が定義された空間における，散逸エネルギー最小化問題が数学的に解析され，形状最適化に必要な framework が数学的に構築された．
有限要素法により支配方程式の弱形式を導き，Lagrange 未定乗数法により Lagrange 関数を定義する．その際，Lagrange 乗数は有限要素法における trial function と結果的に同じものとなる．そして，この Lagrange 関数に対する物質微分をとることで，密度関数の形状微分が得られる．ここでの形状微分は，固定された座標において領域変動に対する密度関数の変動を意味しているため，Fréchet 微分そのものになる．そして，随伴変数法により，Lagrange 関数の物質微分を評価することで目的関数を最小化するための感度（または領域摂動）が得られる．このような解析を一般に感度解析と呼ぶ．
1973 年には，流体力学における要請から，O. Pironneau によって Stokes 問題[2]が，1974 年に定常 Navier-Stokes 問題[3]が定義された空間における散逸エネルギー最小化問題に対する感度解析が行われた．
近年になり，計算機の発達に伴い数値的に流体機器の最適設計を行うことが求められ，形状最適化問題に対する数値計算手法の構築が行われた[4]．ここでは，数値的に評価された感度（または領域摂動）を逐次的に初期の領域に足していくことで，目的関数を最小化するというものであった．しかし，形状更新の度に境界の滑らかさが落ちていく“波打ち現象”が観察されたため，境界の平滑化が必要となった．その際，大きく分けて Parametric 的または Non-Parametric 的な手法が提案された．特に，Non-Parametric 的手法として代表的な手法は，H. Azegami によって提案された H1 勾配法である[5-7]．この H1 勾配法では，領域を線形弾性体と仮定して，感度を境
界で Dirichlet 条件として定義し，領域摂動を平滑化させることを可能にし，誤差解析も行われている．そして，定常 Navier-Stokes 問題に対して，この H1 勾配法を用いた研究が E. Katamine et— 22 —

SENAC Vol. 49, No. 4（2016. 10）

al. [8, 9]及び T. Nakazawa et al. [10-12]によって行われた．更に，非定常 Navier-Stokes 問題に対しては Y. Iwata et al. [13]によって行われた．
流体工学において流体機器の最適設計は実用上，非常に重要な課題であるため，形状最適化問
題が広範囲に活用されてきた分野でもある．特に，航空工学では，2 次元翼断面[14]や 3 次元翼[15]の形状最適化が行われている．その際，揚力と抗力の比である揚抗比が目的関数として用いられている．しかしながら，航空機・高速列車は高速で飛行・走行するため機体の表面において壁面乱流が発生し，この乱流によって飛行・走行安定性を著しく悪化させたり，乱流そのものによる抵抗が予想される．そこで，乱流発生を抑制するために，流れ場の安定性を改善するような形状最適化問題を構築する必要がある．しかし，国内外問わず，そのような研究は著者が知る限り，
ほとんど行われていない．
本研究では，流れ場の安定性を改善するような形状最適化問題を構築することを目的として，散逸エネルギーを目的関数と定義し，非定常 Navier-Stokes 問題を主問題として用いた形状最適化問題を解き，初期形状と最適形状における線形安定性を比較する．その際，初期形状には 2 次
元 Cavity 流れを採用する．本論文は以下のような構成となっている．セクション 2 では，初期形状を定義し，領域の変動を数学的に定義する．そして，主問題である非定常 Navier-Stokes 問題を導入し，形状最適化問題を構築する．そして，線形安定性解析の概要を説明する．セクション
3 では，数値計算結果を示し，セクション 4 では本研究の考察を行う．

2. 問題の設定
本セクションでは，運動エネルギーを目的関数と定義し，非定常 Navier-Stokes 問題を主問題とした際の形状最適化問題と，流れ場の線形安定性解析について述べる．
セクション 2.1 では，初期形状Ωと領域変動�を定義し，汎関数�とその密度関数�に対する物質微分をとることで，汎関数�の形状微分公式を導出する．
セクション 2.2 では，主問題である非定常 Navier-Srokes 問題の弱形式を導き，その際，特性曲線法を導入する．
セクション 2.3 では，Lagarange 未定乗数法を用いるかたちで，運動エネルギーを目的関数，
非定常 Navier-Srokes 問題を主問題とする，Lagarange 関数を定義する．その際，セクション 2.2で弱形式を導出する際に用いた trial function は Lagarange 乗数法に相当する．そして，�の形状

微分公式により，Lagarange 関数の目的関数を最小化させるような領域摂動に対する変分をとり，感動を評価する．
セクション 2.4 では，本研究に対する線形安定性解析について定式化する．

2.1 初期形状と領域変動
ユークリッド空間における，2 次元領域Ω� ∈ ℝ� を非有界な凸領域とする．そして，Ω ⊂ Ω� 初期

形状とし，位置ベクトルを� ∈ ℝ� と記述する．その際，

— 23 —

Shape Optimization による非定常流れ場の安定性制御

Ω = {� = (�, �) ∈ ℝ� | � � � � 1, � � � � 1},Γ��� = {� = (�, �) ∈ ℝ� | � � � � 1, � = 1},�����Γ���� = �Ω � Γ��� .

次に，ℝ� 値の関数である Lipschitz transformation �� Ω � �(Ω)による領域変形を考える．この写像は，微小パラメーター|�| ≪ 1 を用いて，以下のように記述できると仮定する．
� = �� � �� � �(� � )ここで，�� は恒等写像であり�� (�) = �となる．また，�は領域変動である．そして，�の関数空間は M. Kimura[16]の Proposition 1.41 と eq. (1.18)より，

�

�(Ω) = �� ∈ � �,� (Ω; ℝ� ); �(��� )� � ����� � � �

� (�)������� ⊂ Ω� �.
� 1, �(Ω)

これから，汎関数の形状微分公式を導く．まず， � = {�� , �� , � } ∈ �� (Ω)をΩにおける物理量を表す実数値関数とする．そして，�(�, Ω)を密度関数として以下のような汎関数を定義する．
�(�� , �, Ω) = � �(�, Ω) ���また，変形領域における汎関数は以下のように記述することが可能である．

�(�, �(�), �(Ω)) = ��(�(�), �(Ω)) �(��)�(�)この変動領域�(Ω)における�(�, �(�), �(Ω))を，初期領域Ωにおいて記述する．
�の Jacobi 行列は，
(��� )�

��= ����� � � ����� � � �(� � )=����(�� , �� )
� �(� � )�(�� , �� )と記述することが可能であり，Jacobi 行列はdet((��� )� ) = 1 � �� � � � �(� � )となる．よって，以下の結果が得られる．

�(��)= det((��� )� )�� = �1 � �� � � � �(� � )���— 24 —

SENAC Vol. 49, No. 4（2016. 10）

次に，密度関数���(�)�に対して物質微分をとると，
�(�(�), �(Ω)) = �|���� + ���� + (� � �)��|���� + �(� � ).

ここで，( )�は物質微分を(�)は形状微分を表す．更に，形状微分は固定された座標系における領域変動であるので，Fréchet 微分と同意である．

���(�), �(Ω)�と�(��)を�(�, �(�), �(Ω))に代入することで，

���, �(�), �(Ω)� = � �� + ���� + (� � �)�� + �(� � )�|���� �� + �� � � + �(� � )����

= � �� + ���� + (� � �)� + �� � �� + �(� � )�|���� �� .
�となり，���, �(�), �(Ω)� � �(�� , �, Ω)
���
�
lim

�)
= ��(�, �, Ω= � ��� + (� � �)� + �� � ��|���� ���= � ��� + � � (��)�|���� ���= � �� (�, Ω) �� + � �(�, �Ω)� � ��� .
���

ここで， �は境界上の外向き単位法線ベクトルである．また，本研究ではΓ���� が設計境界となっているので，�の関数空間を以下のように再定義する．

�(Ω) = {� � �(Ω)� � = � �� Γ���� }.
そのため，
��(�, �, Ω � Γ���� ) = �� (�� , �, Ω) + �� (�, �, Γ���� ),�� (�� , �, Ω) = � �� (�, Ω) ��,��� (�, �, Γ���� ) = �������(�, Γ���� )� � ��� .

�� (�� , �, Ω)と�� (�, �, Γ���� )は，それぞれ�, �に対して形状微分（Fréchet 微分）をとっており，�� (�, Ω)は微分の連鎖則を用いることで評価することが可能である．

— 25 —

Shape Optimization による非定常流れ場の安定性制御

2.2 形状更新
本研究において，目的汎関数�(�� , �, Ω)の最小化問題を考え，その時の領域変動�を決定する．

その際，�(�� , �, Ω)の変分である��(�, �, Ω ∪ Γ���� )において，�が未知変数であり�(�, Γ���� )�が trial

function という関係になっている．そこで，�を形状最適化問題の目的関数，�, �を任意の微小パ
ラメーター，として，以下のように反復させて目的関数を最小化させるアルゴリズムを採用する．
a.
b.
c.
d.
e.
f.
g.

���
.
Set ���� , ���� � 0, � ��� , � ��� , � ��� , Ω��� , Γ����

Define the functional ���� , � � , Ω� �.

�
Derive the material derivative ����� , � � , Ω� ∪ Γ�����.

�
Evaluate the sensitivity � � �� � , Γ������.

Obtain the new domain Ω��� to reshape Ω� by Ω��� � �� (Ω� ) � ��� .
Operate �� gradient method [5-7].

Judge convergence.

-If the terminal condition �� ��� � � � � � � � ℝ is satisfied, then stop.
-Otherwise, replace � � 1 with � and return to b.

以下，便宜上，形状更新の反復�を�, �, �, �, �, Ω, Γ���� には明記しない．

2.3 主問題
本研究では，非定常 Navier-Stokes 方程式の初期値境界値問題を解くが，計算時間を短縮する
ために定常 Navier-Stokes 方程式の境界値問題の解を初期条件とする．そして，レイノルズ数，
速度ベクトルと圧力を，それぞれRe, � � �(�), � � �(�)とし，

�(�) � �� � �� (Ω; ℝ� )| � � � on Γ���� , � � (1�� � (� � 1)� , 0) on Γ��� �,
�(�) � {� � �� (Ω; ℝ)| ��(�)� �� � 0 �n Ω�.

Problem 1 (Stationary Navier-Stokes Problem)� � �(Ω), � � �(�), � � �(�)のとき，(�, �): Ω → ℝ� × ℝを求めよ．
(� � �)� � ��� �1
��, � � � � 0 �n Ω,Re� � � on Γ���� ,� � (1�� � (� � 1)� , 0)on Γ��� .

� � �(�) � �� � �� (Ω; ℝ� )� � � � on �Ω�, � � �(�)を�, �に対する trial function とすると，Problem 1 の弱形式は，

— 26 —

SENAC Vol. 49, No. 4（2016. 10）

�

�(�)

�{(� � �)�} � � �1�� � �� − (� � �)� − (� � �)�� �� = 0,Re� �(�, �) � �(�, �, �) � �(�, �) � �(�, �) = 0

1
�� � �� ��, �(�, �, �) = � {(� � �)�} � ��� , �(�, �) = − � (� � �)��� .
� Re
�
�

�(�, �) = �

ここで Problem 1 の解を�� ∈ �(�)とする．

Problem 2 (Non Stationary Navier-Stokes Problem)� ∈ �(Ω), � ∈ �(�), � ∈ �(�)のとき，(�, �): Ω × (0, �) → ℝ� × ℝを求めよ．
�� 1− �� � �� = 0, � � � = 0 �n Ω × (0, �)�� Re
� = � on Γ���� × (0, �)� = (1�� � (1 − �)� , 0) on Γ��� × (0, �)� = �� �n Ω �� � = 0.

そして，Problem 2 の弱形式は，� = {0, � , �}, ∀(�, �) ∈ �(�) × �(�)において，��, �� � �(�, �) � �(�, �) � �(�, �) = 0.
���

時間方向の離散化には有限差分法を用い，移流項の取り扱いには特性曲線法を用いる[18]．その際，��を時間刻み，� = ����を全時間ステップとして，時間を� � = ���(� ∈ {0, � , �})と記述す��

る．その結果，(�, �)における物質微分 の 1 時近似は，�(�, �) = �� を既知関数とすると，
���� − ���� �������(�, �) − �(� − ��, � − ��)(�, � � ) =� �(��) =� �(��),��������� = � − ������ (�), and ���� ���� = ���� (��� ).

最終的に，�� = �� を初期条件，{(�� , �� )}���� ∈ �(�) × �(�)を求めるための弱形式は，
��� − ���� ����, �� � �(�� , �� ) � �(�� , �� ) � �(�� , �� ) = 0.
��� � � � �
その際，{(�� , �� )}���� ∈ �(�) × �(�). ここで，(�, �, �, �) = {(� , � , � , � )}��� , �� , �� ∈ ℝとして，
� ,�以下のように汎関数��� � を定義する．
� ,�
��� � (�� , �, �, �, �, Ω) =��1�� − ���� ����� ��, �� � �(�� , �� ) � �(�� , �� ) � �(�� , �� )�.
���� − ��
����

Shape Optimization による非定常流れ場の安定性制御— 27 —

2.4 形状最適化問題
非定常 Navier-Stokes 問題を主問題とした際の，時間平均運動エネルギー最小化問題��11
� � � �� � �� ����=2 �
�� � ��
����

を定義する．

Problem 3 (Shape Optimization Problem Minimizing the Momentum Energy)� � �(Ω), �, � � �(�) � �(�)のとき，
��

11� � � �� � �� ��� | (�� , �� ) in ������� �, {(�� , �� )}�
min �� =��� in ������� ��.
���(�)
�� � ��
2 �
����

Lagrange 乗数法を用いて Problem 3 に対する Lagrange 関数（汎関数）を以下のように定義す
る．

��

1
1
� ,�
�(�� , �, �, �, �, Ω) =� � � �� � �� ��� � ��� � (�� , �, �, �, �, Ω)�� � ��
2 �
��

����

1
� � � � (�, �, �, �, Ω)��=�� � �������そして，形状微分公式を用いると汎関数の物質微分は，
��(�, �, �, �, �, Ω) = �� (�� , �, �, �, �, Ω) � �� (�, �, �, �, �, ����� ),と書ける．更に，微分の連鎖則を用いると�� (�, �, �, �, �, Ω)は，

�� (�� , �, �, �, �, Ω) = ��,� (�� , �, �, �, �, Ω) � ��,� (�� , �, �, �, �, Ω).

�, �に対して形状微分（Fréchet 微分）をとった際の汎関数を��,� (�� , �, �, �, �, Ω)，�, �に対しては��,� (�� , �, �, �, �, Ω)，�に対しては�� (�, �, �, �, �, Ω)，と記述し，具体的には，
��1� � � � (��, ��, �, �, Ω)��,��,� (�� , �, �, �, �, Ω) =�� � ��
�����

— 28 —

SENAC Vol. 49, No. 4（2016. 10）��1��,� (�� , �, �, �, �, Ω) =� � � � (�, �, ��, ��, Ω)��,�� − ��
�
����

��

1
� �� � � (�, �, �, �, ����� )� � ���� .
�� (�, �, �, �, �, ����� ) =�� − �����������

そして，Problem 3 対に対して随伴変数法を用いると，
��,� (�� , �, �, �, �, Ω) = 0 for ∀(�� , �� ) � �(�) × �(�),��,� (�� , �, �, �, �, Ω) = 0 for ∀(�� , �) � �(�) × �(�),となる．��,� (�� , �, �, �, �, Ω) = 0から Problem 2 を解くことになる．��,� (�� , �, �, �, �, Ω) = 0から，(�, �)に対する随伴変数(�, �)を求める必要がある．

Problem 4(Adjoint Problem for Problem 2)�� � �(Ω)とし，{(�� , �� )}���� � �(�) × �(�)が与えられたとき，(�, �): Ω × (�, 0) → ℝ × ℝを求め

よ．

−

�� 1
− �� � �� � � � � = 0, � � � = 0 �� Ω × (�, 0)
�� Re
� = �� �� Ω �� � = ���.

� = � o� �Ω × (�, 0),� � �
Problem 2，Problem 4 から求められた{(�� , �� )}���� � �(�) × �(�)と{(� , � )}��� � �(�) ×

�(�)を��(�, �, �, �, �, ����� )に代入すると，Problem 3 に対する感度は，��1� �� � � (�, �, �, �, ����� )� � ������(�, �, �, �, �, ����� ) =�� − ������������ � (�, �, �, �, ����� ) =1(��� : ��� )|��������
Re

��(�, �, �, �, �, ����� ) = 0となる場合の local minimizer である�を見つけることが出来，�を任意の微小パラメーターとすれば，�(Ω) = �� (Ω) � ��によって最適な形状が得られる．しかし，本研究において，主問題は非線形偏微分発展方程式であるので数学的に� � (�, �, �, �, ����� )を評価することは困難である．そこで，セクション 2.3 で記述した形状更新の反復で� � (�, �, �, �, ����� )を数値的に評価する．

— 29 —

Shape Optimization による非定常流れ場の安定性制御

2.5 線形安定性解析
線形安定性理論に基づき，固有値方程式を導出する．以下，固有値と固有関数を，それぞれ
�(�) × ��(�)とし，
� , �̂ ) � �λ, (��(�) = {�� � �� (�(Ω); ℂ� )| � = � �n �Ω},���(�) = {�̂ � �� (�(Ω); ℂ)| ��(�)�̂ �� = 0 in Ω}.

Problem 5(Linear Stability Analysis)� , �̂ ): Ω → ℂ × ℂ� を求めよ．
� � �(Ω)とし，�� � �(�)が与えられたとき，λ, (��−� + (�� ⋅ ∇)�� + (�� ⋅ ∇)���1
� + ��
� = 0, � ⋅ � = 0 in Ω��Re

� = � �n �Ω.

�(�) × ��(�)を(��, ��) � �� , �̂ )に対する trial function とすると Problem 5 の弱形式は，
∀(�
��

�(�)

� , �) + �(�� , �, �) + �(�� , �� , �) + �(��, �̂ ) + �(�� , ��) = 0.
� ⋅ ��� + �(�3. 数値計算スキーム
本研究における数値計算については，無償の有限要素法ソフトウェアである Freefen++[18]を用いる．数値計算アルゴリズムは以下の通りである．
I.
II.

Prepare a stationary velocity �� by solving Problem 1 with the Newton-Raphson

method at Re = 11500.

Set the stationary velocity �� for the initial domain and conditions, and put � = 0.

III.

Solve Problem 2 and Problem 4, Evaluate the sensitivity,

IV.

Reshape Ω by �(Ω) = �� + �.

V.

Judge convergence.

VI.

-If the terminal condition �� ��� − � � � � 0.01 is satisfied, then proceed to Ⅶ and stop.

-Otherwise, replace � + 1 with � and return to Ⅲ.

Solve Problem 5, and compare linear stability in the initial and optimal domain.

上記の数値計算アルゴリズムを進めるにあたって，空間及び時間方向の離散化を行う．空間方
向には有限要素法により三角形分割を行う．その際，速度，圧力には P2/P1 要素を用いる．また，
時間方向には時間刻み�� により有限差分法を用いている．空間及び時間方向の離散化精度につい— 30 —

SENAC Vol. 49, No. 4（2016. 10）ては以下に記述する．空間方向の離散化精度を担保するために，Γ��� において�� = (1,0)としてNode と Element を増加させ，Problem 1，Problem 5 を解き最大固有値を比較し，Table1 に示した．次に，Node と Element が(8281, 16200)の結果と，[19-23]との結果を比較した．特に，[23]との比較では相対誤差がそれぞれ 1.73%，2.86%となった．そのため，本研究では，Node とElement を(8281, 16200)とする．次に，時間刻み���については，0.001 として，Problem 1，
Problem 2 を解き運動エネルギーを Fig. 1 に示した．Fig. 1 の結果から，��は十分に小さいと仮定する．
Table1. ��� and Imag[�] depending on (♯ of Nodes, ♯ of Elements)(♯ of Nodes, ♯ of Elements)(2601, 5000)

Imag[λ]

Re�

7310

±3.02

7515

±2.95

(5041, 9800)

7594±2.916(6501, 12800)

7741

±2.89

(8281, 16200)

7890

±2.86

(3721, 7200)

Table2. ��� and Imag[�] in [19-23].

References

Imag[λ]Re�

[19]7763.4±2.8634i

[20]

8000

±2.8356i[21]8375.0±2.7640i

[22]

7890

±2.7646i[23]8026.6±2.82559i
0.04734

Momentum Energy

Momentum Energy

0.05

0.049

0.048

0.047

0

100

200

300

400

500

600

time

(a)�� = (0,1000).

700

800

900 1000

0.04733

0.04732

600

601

602

603

604

605

606

607

time

(b) � = (600,61).

Fig. 1. The momentum energy, (a) � = (�, ����) and (b) � = (���, ��).

608

609

610

— 31 —

Shape Optimization による非定常流れ場の安定性制御

4. 数値計算結果
セクション 3 で示した，数値計算アルゴリズムに従がって形状最適化問題を解いた．Fig. 2 と� の流れ関数の分布である．Fig.
Fig. 3 は，それぞれ初期形状と最適形状における定常解�� の擾乱�

3 からΓ���� が領域の内部に向かって移動していることがわかる．Fig. 4 は，初期形状と最適形状の最大固有値の実部を図示しており，最適形状の方がRe� が小さくなった．Fig. 5 は，Re = 11500における固有値の分布を示しており，Fig. 4 の結果を補強している．

(a)

(b)

Fig. 2. Stream functions of (a) � and (b) � in the initial domain.
�(a)Fig. 3. Stream functions of (a) �� and (b) � in the optimal domain.
7

Ω0
Ω1

0.06
0.05
0.04
0.03
0.02
0.01
0

5
4
3
2
1

-0.01
-0.02
6500

Ω0
Ω1

6

Imaginary parts

Real parts of the eigenvalues

0.07

(b)

0

7000

7500

Re

8000

8500

9000

Fig. 4. Real parts of the leading eigenvalues

-0.2

-0.15

-0.1

-0.05

Real parts

0

0.05

0.1

Fig. 5. Spectrum of the eigenvalues
at Re = 11500.

— 32 —

SENAC Vol. 49, No. 4（2016. 10）

5. 考察
本研究では，運動エネルギー最小化問題を行い，初期形状と最適形状における線形安定性を比較する．その際，初期形状を 2 次元 Cavity 流れとし，レイノルズ数Re = 11500において非定常
Navier-Stokes 問題を主問題とする．数値計算には無償の有限要素法ソフトウェアである
Freefen++[18]を用いた．

数値計算の結果，Γ���� が領域の内部に向かって移動していた．更に，初期形状と比較して最適形状における臨界レイノルズ数Re� が減少したことが数値的に示された．一般に，流れ場は剪断応
力が大きななった際に不安定になるとされている．目的関数を最小化する過程でΓ���� が領域の内部に移動することによってΓ���� 上で剪断応力が大きくなり，流れ場が不安定になったと考えられる．その結果，運動エネルギーを最小化しているにも拘らずRe� が減少した．

本研究は，2次元Cavity流れを対象としているが，将来的には，3次元の大規模高精度計算により複雑形状内流体の形状最適化を行う予定である．しかし，そのためには多くの計算コストが必要となるため，今後は東北大学サイバーサイエンスセンターのスーパーコンピュータを利用することで実現することが可能であると考えている．

謝辞
今後の研究にのために，東北大学サイバーサイエンスセンターのスーパーコンピュータ関係各位に，Freefem++をインストールして頂いた．これにより， Problem 2 を実行するに至った．現
在は，2 次元空間における数値解析しか行ってないが将来的には 3 次元空間の解析に役立てていきたい．更に，情報科学研究科とサイバーサイエンスセンターとの申し合わせにより研究室に負担金が生じることなく大規模な計算機が利用できる仕組みがあるので，大規模高精度な計算に積極的に利用したい．これは，著者の今後の研究にとって非常に有意義なことであり，感謝申し上
げたい．

参考文献
[1] J. Hadamard, Memoire sur un probleme d’analyse relatif a l’equilibre des plaques elastiques
encastrees, Memoire des savants etragers, Oeuvres de J. Hadamard, CNRS, Paris, 1968
[2] O. Pironneau, On optimum profiles in Stokes flow, JFM, 59, 117-128, 1973
[3] O. Pironneau, On optimum design in fluid mechanics, JFM, 64, 97-110, 1974
[4] O. Pironneau, Optimal shape design for elliptic systems, Springer-Verlag, New York, 1984.
[5] H. Azegami, S. Kaizu, M. Shimoda, E. Katamine, Irregularity of shape optimization prob-lemsand an improvement technique, Computer Aided Optimization Design of Structures V, Chap. 4,
Sec. 4, edited by Hernandez, S. and Brebbia, C. A., Computational Mechanics Publications,Southampton, 309-326, 1997
[6] H. Azegami, Z. Wu, Domain Optimization analysis in linear elastic problems: Approach using
traction Method, JSME Int. J. Series A, 39, 272-278, 1996

Shape Optimization による非定常流れ場の安定性制御— 33 —[7] H. Azegami, Regularized Solution to Shape Optimization Problem, Transactions of the Japan
Society for Industrial and Applied Mathematics, 23, 83-138, 2014[8] E. Katamine, Y. Nagatomo, H. Azegami, Shape optimization of 3D viscous flow fields, Inverse
Problems in Sci. and Eng., 17, 105-114, 2009.
[9] E. Katamine, H. Azegami, T. Tsubata, S. Itoh, Solution to shape optimization problems ofviscous flow fields, Int. J. of Com. Fluid Dyn., 19, 45-51, 2005.
[10] T. Nakazawa, Eiji Katamine, Akiyasu Tomoeda, Hideyuki Azegami，An application of shape
optimization method to a room design based on anincompressible fluid model, RIMS講究録別冊「トポロジー最適化理論とその現実問題への応用」, B(54), 1-15, 2015.
[11] T. Nakazawa, Increasing the Critical Reynolds number by maximizing dissipation energy
problem, Proc. of the 5th ICJWSF2015, Editor Antonio Segalini, Chapter 75, 2016[12] T. Nakazawa, Shape Optimization of the maximizing problem of the dissipation energy and its
effect on Hydrodynamic stability, Proc. of the 15th ETC2015, 477, Netherland, 2015[13] Y. Iwata, H. Azegami, T. Aoyama, E. Katamine, Numerical solution to shape optimization
problems for non-stationary Navier–Stokes problems, JSIAM Letters, 2, 37-40, 2010[14] D. Li, R. Hartman, Adjoint-based airfoil optimization with discretization error control, Int. J.
Numer. Meth. Fluids, 77, 1-17, 2015[15] S. Schmidt, C. Ilic, V. Schulz, N. Gauger, Three-dimensional large-scale aerodynamic shape
optimization based on shape calculus. AIAA Journal, 51, 2615-2627, 2013[16] M. Kimura, Lecture notes Volume IV, Topics in Mathematical Modeling, 1-38, 2008[17] H. Notsu, M. Tabata, Error estimates of a pressure-stabilized characteristics finite element
scheme for Oseen equations, J. Sci. Comput, 65, 940-955, 2015[18] Hecht, F. New development in FreeFem++. J. Numer. Math. 20, no. 3-4, 251–265. 65Y15, 2012[19] Poliashenko M, Aidun CK. A direct method for computation of simple bifurcations, J.l ofCom. Physics, 121, 246–260, 1995

[20] Fortin A, Jardak M, Gervais JJ, Pierre R. Localization of Hopf bifurcations in fluid flow
problems. Int. J. for Num. Methods in Fluids, 24, 1185–1210, 1997[21] Tiesinga G, Wubs FW, Veldman AEP. Bifurcation analysis of incompressible flow in adriven cavity by the Newton–Picard method, J. of Comp and Appl. Math., 140:751–772,
2002[22] Cadou JM, Potier-Ferry M, Cochelin B. A numerical method for the computation ofbifurcation points in fluid mechanics. Euro. J. of Mech. Fluids B, 25, 234–254, 2006[23] V. B. L. Boppana and J. S. B. Gajjar, Global flow instability in a lid-driven cavity, Int. J.
Numer. Meth. Fluids, 2009