論文

気液二相流れ場におけるパラメータ同定解析小林 雅也1・倉橋 貴彦2・剱地 利昭3・衛藤 俊彦41

長岡技術科学大学大学院 工学専攻 機械工学分野(Mechanical Engineering, Master’s Program in Engineering, Graduate School of Nagaoka University of Technology)
2

長岡技術科学大学 技学研究院 機械系（Department of Mechanical Engineering, Institute of GIGAKU, Nagaoka University of Technology）
3

函館工業高等専門学校 生産システム工学科（Department of Production System Engineering, National Institute of Technology, Hakodate College）
4

環境都市工学科（Department of Civil Engineering, National Institute of Technology, Nagaoka College）

Parameter identification analysis in gas-liquid two-phase flow field
Masaya KOBAYASHI1, Takahiko KURAHASHI 2, Toshiaki KENCHI3
and Toshihiko ETO4
Abstract
In this paper, we present the results of optimization in a gas-liquid two-phase flow field based on theadjoint variable method. When performing optimization in a gas-liquid two-phase flow field, the advection
equation to determine the interface must be added to Lagrangian as an equality constraint. The adjoint
advection equation is then derived, but it is clear that the adjoint variables in the advection equation are notrelated to the other stationary conditions and need not be solved. The analytical conditions used in thisanalysis have the characteristics of advection dominance, and a center-difference approximation to theadvection term results in numerical instability. To solve this problem, the flow field equations and theadjoint flow field equation were discretized based on the SLG method. Numerical analysis was performed
using the programing language FreeFEM++.
Key Words : Two-phase flow, Parameter identification, FreeFEM++

１．はじめに場での最適化を行う必要があるが，これに関する研究は多くない．

製品設計の際，性能向上や効率の良い製品設計をそこで，本論文では，気液二相流れ場における最目的として，パラメトリック最適化および形状最適適化を行った結果を示す．本解析では，先に示した化やトポロジー最適化に代表されるノンパラメトリ

形状最適化問題から，問題としての難易度を落とし，ック最適化といった最適化手法が用いられる．これ

非定常非圧縮性粘性流体を対象とした気液二相流れ

ら手法の適用先として，液体輸送トラックや液体輸場の部分領域において，目的とする流速を得るため，送船舶の走行安定性向上を目的として，タンク形状随伴変数法に基づき液体の粘度を決定するパラメーの最適化や，タンクに取り付けられたバッフル板形タ同定解析を行った．気液二相流れ場における最適状の最適化が挙げられる．この場合，気液二相流れ化の際には，界面を決定するための移流方程式を等28式制約としてラグランジュ関数に加える必要がある．

Φ=

それに伴い，随伴移流方程式が導出されるが，移流方程式にかかる随伴変数がその他停留条件に関係し1 tan−1 (𝐻𝐻)+𝜋𝜋2(6)ないことが明らかとなり，解く必要が無いことが確

流れ場の支配方程式として式(7)のナビエ・ストーク認できる．また本解析で用いた解析条件は移流卓越

ス方程式と式(8)の連続の式を与える．連続の式に関の特性を持ち，移流項に対して中心差分的近似を行して，非圧縮条件による数値的不安定性を解消する
ため，非常に小さなペナルティ係数𝜀𝜀𝑝𝑝 = 1.0 × 10−8

った場合，数値的不安定性が生じる．これを解消するため，流れ場の方程式と随伴流れ場の方程式はをかけたペナルティ項を付与している．ここで，SLG 法 1) 2)に基づき離散化した．数値解析はプログ
ラミング言語 FreeFEM++3) 4)を用いて行った．

𝑢𝑢𝑖𝑖 [m/s]は流速，𝑝𝑝[Pa]は圧力，𝑔𝑔 = −9.8[m/s 2 ]は重力加速度，𝛿𝛿𝑖𝑖𝑖𝑖 はクロネッカーのデルタを表す．

２．パラメータ同定解析の定式化

𝜌𝜌２．１

気液二相流れ場における界面の捕捉

𝐷𝐷𝑢𝑢𝑖𝑖+ 𝑝𝑝,𝑖𝑖 − 𝜇𝜇𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗 − 𝜌𝜌𝜌𝜌𝛿𝛿𝑖𝑖2 = 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � in Ω(7)𝐷𝐷𝐷𝐷(8)𝑢𝑢𝑖𝑖,𝑖𝑖 + 𝜀𝜀𝑝𝑝 𝑝𝑝 = 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � in Ω

本解析ではVOF法5)を用いて気液二相流れ場の界気液二相流れ場の支配方程式に与える初期条件，境面を捕捉する．VOF法は界面位置をVOF関数Φと呼界条件はそれぞれ式(9)~(12)のように定義する．ここ� は既知の関数，𝑛𝑛𝑖𝑖 は境界における𝑖𝑖方向の単位で，(∙)

ばれるスカラー関数について，式(1)に示す移流方程式を解くことにより捕捉する．ここで，𝐷𝐷 ∙/𝐷𝐷𝐷𝐷は物質微分，𝑡𝑡0 は初期時刻，𝑡𝑡𝑓𝑓 は終端時刻，Ωは計算領域を法線ベクトル，Γ1 はディリクレ境界，Γ2 はノイマン境界を表す．

表す．

𝐷𝐷Φ
= 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � in Ω𝐷𝐷𝐷𝐷

� at 𝑡𝑡 = 𝑡𝑡0 in Ω
Φ=Φ𝑢𝑢𝑖𝑖 = 𝑢𝑢�𝚤𝚤 at 𝑡𝑡 = 𝑡𝑡0 in Ω𝑢𝑢𝑖𝑖 = 𝑢𝑢�𝚤𝚤 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ1(1)𝑇𝑇𝑖𝑖 = 𝑇𝑇�𝚤𝚤 = �−𝑝𝑝𝛿𝛿𝑖𝑖𝑖𝑖 + 𝜇𝜇𝑢𝑢𝑖𝑖,𝑗𝑗 �𝑛𝑛𝑗𝑗 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ2 (12)VOF関数Φは流体であれば1，気体であれば0となるように定義する．このとき，密度𝜌𝜌[kg/m3 ]と粘度２．２𝜇𝜇[Pa ∙ s]はVOF関数Φを用いて式(2)および式(3)のよ(2)(3)

随伴変数法に基づき液体の粘度𝜇𝜇𝑙𝑙 を更新する．ここで，𝑄𝑄は位置を決定するための重み係数であり式(14)𝐽𝐽 =るためデジタイザ法6)に基づき界面鋭敏化を行う．デジタイザ法では式(1)に示すVOF関数Φに関する移流方程式を解かずに，式(4)によってタンジェント変換された式(5)に示すVOF関数𝐻𝐻に関する移流方程式を解く．その後，式(6)によって逆タンジェント変換し1 𝑡𝑡𝑓𝑓2� � 𝑄𝑄�𝑢𝑢𝑖𝑖 − 𝑢𝑢𝑖𝑖(𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡) � 𝑑𝑑Ω 𝑑𝑑𝑑𝑑
2 𝑡𝑡0 Ω𝑄𝑄 = �1 on Ω𝑐𝑐0 on Ω ∖ Ω𝑐𝑐

図-1 に解析モデルを示す．

界面を決定する．

𝐷𝐷𝐷𝐷= 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � in Ω𝐷𝐷𝐷𝐷

評 価関 数 計 算領 域 Ωc に お い て ，目 的 と する 粘 度
𝜇𝜇𝑙𝑙(𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡) での流速𝑢𝑢𝑖𝑖(𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡) との差を最小とするようで表す．

本解析では，数値拡散による界面のぼやけを解消す

1𝐻𝐻 = tan �𝜋𝜋 �Φ − ��2

評価関数

評価関数𝐽𝐽を式(13)のように定義する．本解析ではうに決定される．ここで，𝜌𝜌𝑙𝑙 と𝜇𝜇𝑙𝑙 は液体の密度と粘度，
𝜌𝜌𝑔𝑔 と𝜇𝜇𝑔𝑔 は気体の密度と粘度を表す．
𝜌𝜌 = 𝜌𝜌𝑙𝑙 Φ + 𝜌𝜌𝑔𝑔 (1 − Φ)𝜇𝜇 = μl Φ + 𝜇𝜇𝑔𝑔 (1 − Φ)

(9)
(10)
(11)

(4)

(5)

29(13)(14)

+ � Φ∗ �𝑡𝑡𝑓𝑓 �𝛿𝛿Φ�𝑡𝑡𝑓𝑓 � 𝑑𝑑Ω − � Φ∗ (𝑡𝑡0 )𝛿𝛿Φ(𝑡𝑡0 ) 𝑑𝑑Ω
Ω𝑡𝑡𝑓𝑓

Ω∗
+ � �𝛿𝛿𝑢𝑢𝑖𝑖 �𝑝𝑝∗ 𝛿𝛿𝑖𝑖𝑖𝑖 + 𝜇𝜇𝑢𝑢𝑖𝑖,𝑗𝑗
� 𝑛𝑛𝑗𝑗 𝑑𝑑Γ 𝑑𝑑𝑑𝑑
𝑡𝑡0 Γ
𝑡𝑡𝑓𝑓

𝑡𝑡𝑓𝑓

+ � �𝑢𝑢𝑖𝑖∗ 𝑛𝑛𝑖𝑖 𝛿𝛿𝛿𝛿𝛿𝛿Γ 𝑑𝑑𝑑𝑑 + � �−𝜇𝜇𝑢𝑢𝑖𝑖∗ 𝛿𝛿𝑢𝑢𝑖𝑖,𝑗𝑗 𝑛𝑛𝑗𝑗 𝑑𝑑Γ 𝑑𝑑𝑑𝑑
𝑡𝑡0 Γ
𝑡𝑡𝑓𝑓

𝑡𝑡0

∗
+ � � 𝛷𝛷𝑢𝑢𝑖𝑖,𝑗𝑗
𝑢𝑢𝑖𝑖,𝑗𝑗 𝛿𝛿𝜇𝜇𝑙𝑙 𝑑𝑑Ω 𝑑𝑑𝑑𝑑
Ω𝑡𝑡0

図-1Γ

(16)

ラグランジュ関数第一変分の停留条件より，式(17)気液二相流れ場の解析モデルおよび式(18)の随伴流れ場方程式，式(19)の随伴移流２．３方程式と式(20)~(22)で表されるその初期条件および随伴方程式および感度の導出評価関数と気液二相流れ場決定のための支配方程境界条件が得られる．このとき，随伴VOF関数を決∗式よりラグランジュ関数𝐽𝐽 を式(15)に定義する．ここ定するための随伴移流方程式(式(19))が導出されるで，𝑢𝑢𝑖𝑖∗ ，𝑝𝑝∗ ，Φ∗ はそれぞれ，随伴流速，随伴圧力，がラグランジュ関数第一変分のその他停留条件に関随伴VOF関数である．

係しないため，方程式は解かずに方程式がゼロとなるように随伴VOF関数が求められたとする．

∗𝐽𝐽 = 𝐽𝐽
𝑡𝑡𝑓𝑓

+ � � 𝑢𝑢𝑖𝑖∗ �𝜌𝜌
𝑡𝑡0 Ω𝑡𝑡𝑓𝑓

𝐷𝐷𝑢𝑢𝑖𝑖+ 𝑝𝑝,𝑖𝑖 − 𝜇𝜇𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗 − 𝜌𝜌𝜌𝜌𝛿𝛿𝑖𝑖2 � 𝑑𝑑Ω 𝑑𝑑𝑑𝑑𝐷𝐷𝐷𝐷+ � � 𝑝𝑝∗ �𝑢𝑢𝑖𝑖,𝑖𝑖 + 𝜀𝜀𝑝𝑝 𝑝𝑝� 𝑑𝑑Ω 𝑑𝑑𝑑𝑑
𝑡𝑡0 Ω𝑡𝑡𝑓𝑓+ � � Φ∗ �𝑡𝑡0Ω𝐷𝐷Φ� 𝑑𝑑Ω 𝑑𝑑𝑑𝑑
𝐷𝐷𝐷𝐷−𝜌𝜌
(15)

ラグランジュ関数の第一変分𝛿𝛿𝐽𝐽∗ を式(16)に示す．
𝑡𝑡𝑓𝑓

𝑡𝑡0

Ω+ � � 𝛿𝛿𝑢𝑢𝑖𝑖∗ �𝜌𝜌
𝑡𝑡0 Ω
𝑡𝑡𝑓𝑓

𝐷𝐷𝑢𝑢𝑖𝑖+ 𝑝𝑝,𝑖𝑖 − 𝜇𝜇𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗 − 𝜌𝜌𝜌𝜌𝛿𝛿𝑖𝑖2 � 𝑑𝑑Ω 𝑑𝑑𝑑𝑑
𝐷𝐷𝐷𝐷+ � � 𝛿𝛿𝑢𝑢𝑖𝑖 �−𝜌𝜌
𝑡𝑡0 Ω
𝑡𝑡𝑓𝑓

𝐺𝐺𝑖𝑖 を示す．

𝐷𝐷𝑢𝑢𝑖𝑖∗∗+ 𝑝𝑝,𝑖𝑖∗ − 𝜇𝜇𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗
� 𝑑𝑑Ω 𝑑𝑑𝑑𝑑
𝐷𝐷𝐷𝐷+ � � 𝛿𝛿Φ∗ �𝑡𝑡0 Ω
𝑡𝑡𝑓𝑓
𝑡𝑡0 Ω
𝑡𝑡𝑓𝑓
𝑡𝑡0ΩΩ𝛺𝛺(𝑙𝑙+1)

２．４

on Γ2 (22)(23)

𝜇𝜇𝑙𝑙(𝑙𝑙)= 𝜇𝜇𝑙𝑙 − 𝛼𝛼𝐺𝐺𝑖𝑖

方程式の離散化(24)

本解析で解く必要のある方程式の多くは移流項を持つ．方程式を離散化する際，移流項に対して中心− (𝜌𝜌1 − 𝜌𝜌2 )𝑔𝑔𝛿𝛿𝑖𝑖2 � 𝑑𝑑Ω 𝑑𝑑𝑑𝑑Ω(21)

上付き(𝑙𝑙)は反復回数である．

𝐷𝐷𝑢𝑢𝑖𝑖− (𝜇𝜇1 − 𝜇𝜇2 )𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗
𝐷𝐷𝐷𝐷+ � 𝑢𝑢𝑖𝑖∗ �𝑡𝑡𝑓𝑓 �𝛿𝛿𝛿𝛿𝑖𝑖 �𝑡𝑡𝑓𝑓 � 𝑑𝑑Ω − � 𝑢𝑢𝑖𝑖∗ (𝑡𝑡0 )𝛿𝛿𝑢𝑢𝑖𝑖 (𝑡𝑡0 ) 𝑑𝑑Ω

in Ω (19)(20)

式(24)に粘度の更新式を示す．ここで，αは更新幅，𝐷𝐷Φ∗� 𝑑𝑑Ω 𝑑𝑑𝑑𝑑
𝐷𝐷𝐷𝐷+ � � 𝛿𝛿Φ �𝑢𝑢𝑖𝑖∗ (𝜌𝜌1 − 𝜌𝜌2 )𝑡𝑡𝑓𝑓

𝑡𝑡0𝐷𝐷Φ
� 𝑑𝑑Ω 𝑑𝑑𝑑𝑑
𝐷𝐷𝐷𝐷+ � � 𝛿𝛿Φ �−𝐷𝐷𝐷𝐷− (𝜌𝜌1 − 𝜌𝜌2 )𝑔𝑔𝛿𝛿𝑖𝑖2 = 0∗𝑢𝑢𝑖𝑖,𝑗𝑗 𝛿𝛿𝜇𝜇𝑙𝑙 𝑑𝑑Ω 𝑑𝑑𝑑𝑑
𝐺𝐺𝑖𝑖 = � � Φ𝑢𝑢𝑖𝑖,𝑗𝑗∗+ � � 𝛿𝛿𝛿𝛿�−𝑢𝑢𝑖𝑖,𝑖𝑖+ 𝜀𝜀𝑝𝑝 𝑝𝑝∗ � 𝑑𝑑Ω 𝑑𝑑𝑑𝑑
𝑡𝑡0 Ω𝑡𝑡𝑓𝑓

𝐷𝐷𝐷𝐷

(17)
(18)

式(23)に設計変数に対するラグランジュ関数の勾配+ � � 𝛿𝛿𝑝𝑝∗ �𝑢𝑢𝑖𝑖,𝑖𝑖 + 𝜀𝜀𝑝𝑝 𝑝𝑝� 𝑑𝑑Ω 𝑑𝑑𝑑𝑑
𝑡𝑡0 Ω
𝑡𝑡𝑓𝑓

in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � in Ω∗∗−𝑢𝑢𝑖𝑖,𝑖𝑖 + 𝜀𝜀𝑝𝑝 𝑝𝑝 = 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � in Ω𝐷𝐷Φ∗𝐷𝐷𝑢𝑢𝑖𝑖−+ 𝑢𝑢𝑖𝑖∗ (𝜌𝜌1 − 𝜌𝜌2 )− (𝜇𝜇1 − 𝜇𝜇2 )𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 �∗∗�𝑢𝑢𝑖𝑖 = 𝑢𝑢𝚤𝚤 at 𝑡𝑡 = 𝑡𝑡𝑓𝑓 in Ω�𝚤𝚤∗ in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ1𝑢𝑢𝑖𝑖∗ = 𝑢𝑢∗∗∗∗�𝑇𝑇𝑖𝑖 = 𝑇𝑇
𝚤𝚤 = �𝑝𝑝 𝛿𝛿𝑖𝑖𝑖𝑖 + 𝜇𝜇𝑢𝑢𝑖𝑖,𝑗𝑗 �𝑛𝑛𝑗𝑗 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 �𝛿𝛿𝐽𝐽∗ = � � 𝑄𝑄�𝑢𝑢𝑖𝑖 − 𝑢𝑢𝑖𝑖(𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡) �𝛿𝛿𝑢𝑢𝑖𝑖 𝑑𝑑Ω 𝑑𝑑𝑑𝑑
𝑡𝑡𝑓𝑓

𝐷𝐷𝑢𝑢𝑖𝑖∗∗+ 𝑝𝑝,𝑖𝑖∗ − 𝜇𝜇𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗+ 𝑄𝑄�𝑢𝑢𝑖𝑖 − 𝑢𝑢𝑖𝑖(𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡) � = 0𝐷𝐷𝐷𝐷

差分的近似で離散化を行った場合，数値的不安定になることが知られている 5)．従って，本解析ではこの数値的不安定性を解消するため SLG 法に基づき30し step4 に戻る．

方程式を離散化する．SLG 法では時間微分項と移流項を物質微分の形で表し，物質微分項を時間方向に離散化することで移流計算と非移流計算を分離する．

３．数値解析例

移流計算に関しては，あらかじめ特性曲線法に基づいて解を求め，求めた解と非移流計算は有限要素法に基づき離散化を行うことで数値的不安定性を解消

数値解析は解析条件を変えて case1，2，3 の 3 通することができる．式(25)~(29)にある時間ステップり行った．数値解析におけるプログラムはにおける離散化式を示す．上からナビエ・ストーク

FreeFEM++を用いて作成した．

ス方程式と連続の式，移流方程式，随伴方程式であ� は特性曲線上の上流点における物理る．ここで，(∙)

３．１

図-2 に有限要素メッシュ，図-3 に VOF 関数の初量，𝑤𝑤𝑖𝑖 ，𝑞𝑞，ℎは重み関数，上付きの𝑛𝑛は時間ステップ

期 分 布 と 評 価 関 数 計 算 領 域 Ωc を 示 す ． 一 辺 𝐿𝐿 をを表す．

� �𝜌𝜌𝑤𝑤𝑖𝑖 �Ω𝑢𝑢𝑖𝑖𝑛𝑛+1 − 𝑢𝑢�𝚤𝚤
∆𝑡𝑡

𝑛𝑛+1
� − 𝑤𝑤𝑖𝑖,𝑖𝑖 𝑝𝑝𝑛𝑛+1 + 𝜇𝜇𝑤𝑤𝑖𝑖,𝑗𝑗 𝑢𝑢𝑖𝑖,𝑗𝑗− 𝜌𝜌𝜌𝜌𝑤𝑤𝑖𝑖 𝛿𝛿𝑖𝑖2 � 𝑑𝑑Ω = 0𝑛𝑛+1+ 𝜀𝜀𝑝𝑝 𝑞𝑞𝑝𝑝𝑛𝑛+1 � 𝑑𝑑Ω = 0
� �−𝑞𝑞𝑢𝑢𝑖𝑖,𝑖𝑖Ω� �𝜌𝜌𝑤𝑤𝑖𝑖 �Ω

� ℎ�

Ω
�𝚤𝚤∗
𝑢𝑢𝑖𝑖∗𝑛𝑛 − 𝑢𝑢

∆𝑡𝑡

�𝚤𝚤
𝐻𝐻𝑖𝑖𝑛𝑛+1 − 𝐻𝐻∆𝑡𝑡

� 𝑑𝑑Ω = 0∗𝑛𝑛
� − 𝑤𝑤𝑖𝑖,𝑖𝑖 𝑝𝑝 ∗𝑛𝑛 + 𝜇𝜇𝑤𝑤𝑖𝑖,𝑗𝑗 𝑢𝑢𝑖𝑖,𝑗𝑗Ω

２．５𝐿𝐿=0.6[m]とする二次元ダム崩壊問題を設定した．解析条件は表-1 および表-2 に示す通りである．メッシ

ュは一辺を 160 分割するよう buildmesh 関数により生成した．流速，圧力，VOF 関数，随伴流速，随伴(25)圧力，随伴 VOF 関数はすべて三角形一次要素により補間した．境界条件として，全ての境界に対して，(26)流速をゼロとするノンスリップ境界条件を与えた．
特性曲線上の上流点における物理量は convect 関数(27)により求めた．なお，ナビエ・ストークス方程式と随伴方程式における移流流速𝑎𝑎𝑖𝑖 ，𝑎𝑎𝑖𝑖∗ は非線形性を回避するため二次精度の Adams-Bashforth 法により式(31)および式(32)のように近似した．

𝑛𝑛+ 𝑄𝑄𝑤𝑤𝑖𝑖 �𝑢𝑢𝑖𝑖𝑛𝑛 − 𝑢𝑢𝑖𝑖(𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡)
�� 𝑑𝑑Ω = 0 (28)∗𝑛𝑛− 𝜀𝜀𝑝𝑝 𝑞𝑞𝑝𝑝𝑛𝑛 � 𝑑𝑑Ω = 0� �−𝑞𝑞𝑢𝑢𝑖𝑖,𝑖𝑖解析条件(29)

𝑎𝑎𝑖𝑖 =3 𝑛𝑛 1 𝑛𝑛−1
𝑢𝑢 − 𝑢𝑢
2 𝑖𝑖 2 𝑖𝑖

31𝑎𝑎𝑖𝑖∗ = 𝑢𝑢𝑖𝑖∗𝑛𝑛+1 − 𝑢𝑢𝑖𝑖∗𝑛𝑛+222

解析の流れ

以下に計算フローを示す．
Step1.計算条件を入力する．
Step2.目的とする粘度で気液二相流れ場の解析を行い，目的とする粘度での流速を得る．
Step3.初期粘度を入力する．
Step4.与えられた粘度で気液二相流れ場の解析を行い，流速を得る．
Step5.評価関数を計算し，式(30)に示す収束判定条件
式を満たした場合，解析を終了する．満たさない場合は step6 に進む．ここで，𝜀𝜀は収束判定定数である．
𝐽𝐽(𝑙𝑙) − 𝐽𝐽(𝑙𝑙−1)
�� < 𝜀𝜀𝐽𝐽(0)

(30)

図-2

Step6.随伴流れ場の解析を行い，感度を計算する．
Step7.上記 step6 で求められた感度により粘度を更新

31有限要素メッシュ(31)(32)

たし 0.004[%]となった．設計変数である液体の粘度
反𝜇𝜇1 も初期値 1[Pa ∙ s]から目的とする粘度に近づき，復回数𝑙𝑙=62 で 0.002[Pa ∙ s]となった．

図-3

解析モデルの寸法とVOF関数の初期分布
case1での数値解析条件

表-1

𝑡𝑡𝑓𝑓

終端時間

𝛥𝛥𝛥𝛥

時間刻み幅
気体の密度

液体の粘度
(目的値)
液体の粘度

1.5 × 10−5𝜇𝜇𝑙𝑙(𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡)[Pa ∙ s]1.0 × 10−3𝜇𝜇𝑙𝑙

(初期値)

𝛼𝛼

更新幅
収束判定定数表-2[Pa ∙ s][kg/m3 ][Pa ∙ s]1.0 × 105

1.0 × 10−5Case2
時間刻み幅
気体の密度
気体の粘度
液体の密度
液体の粘度
(目的値)
液体の粘度
(初期値)
更新幅
収束 0 判定定数

３．２𝛼𝛼𝜀𝜀

表-31.0𝑙𝑙=030

1030

case1 での数値解析結果

評価関数 [%]

1.0𝜌𝜌𝑙𝑙𝜇𝜇𝑙𝑙

表-3 に case1 における解析結果を示す．

0.2𝜇𝜇𝑔𝑔𝜇𝜇𝑙𝑙(𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡)

Case3
10𝛥𝛥𝛥𝛥𝜌𝜌𝑔𝑔

case1 における液体の粘度の推移

1.0 × 103
1.0

𝜀𝜀𝑡𝑡𝑓𝑓図-5

1.25

case2,3での数値解析条件

終端時間

case1 における評価関数の推移

0.001

𝜌𝜌𝑔𝑔
𝜌𝜌𝑙𝑙

液体の密度[s]図-40.5[kg/m3 ]𝜇𝜇𝑔𝑔

気体の粘度[s]𝑙𝑙=621000.004

粘度𝜇𝜇1 [Pa ∙ s]
10.002

図-6~8 に𝜇𝜇𝑙𝑙 =0.001[Pa ∙ s] (目的値)，=1[Pa ∙ s] (初期50値)，=0.002[Pa ∙ s] (同定された値)での VOF 関数の0.1[s]ごとの時間経過を示す．

2.0 × 108

1.0 × 10−5

解析結果

・case1
図-4 に評価関数の推移，図-5 に液体の粘度の推移を示す．評価関数は，初期値 100[%]から粘度の更新とともに減少し，反復回数𝑙𝑙=62 で収束判定条件を満32

図-6 𝜇𝜇𝑙𝑙 =0.001[Pa ∙ s] (目的値)での VOF 関数の遷移

図-8 𝜇𝜇𝑙𝑙 =0.002[Pa ∙ s] (同定値)での VOF 関数の遷移

・case2，3

図-9 に評価関数の推移，図-10 に液体の粘度の推移を示す．case2 では初期値に目的値より小さな値を設定し，case3 では大きな値を設定した．解析結果はいずれの場合も目的値へと漸近した．これにより初期値によらず，パラメータの同定を行うことができていることがわかる．

図-9

図-7 𝜇𝜇𝑙𝑙 =1.0[Pa ∙ s] (初期値)での VOF 関数の遷移

33

case2，3 における評価関数の推移

図-10

case2，3 における液体の粘度の推移

表-4 および表-5 に case2，3 での解析結果を示す．
表-4

case2 での数値解析結果

評価関数 [%]
100𝑙𝑙=0𝑙𝑙=540.045
表-5𝑙𝑙=147

1029.665

case3 での数値解析結果

評価関数 [%]
𝑙𝑙=0

粘度𝜇𝜇𝑙𝑙 [Pa ∙ s]

1000.739

粘度𝜇𝜇𝑙𝑙 [Pa ∙ s]

図-11 𝜇𝜇𝑙𝑙 =30[Pa ∙ s] (目的値)での VOF 関数の遷移

5030.839

図-11~13 に𝜇𝜇𝑙𝑙 =30[Pa ∙ s] (目的値)，=10[Pa ∙ s] (case2初期値)，=50[Pa ∙ s] (case3 初期値)での VOF 関数の2[s]ごとの時間経過を示す．

４．結論
本論文では，気液二相流れ場でのパラメータ同定における定式化と解析結果について示した．定式化
では，随伴移流方程式が導出されたが，その他停留条件に関係しないことが明らかとなり，解く必要が無いことが確認できた．また，case1~3 の解析結果を通して理論とプログラムの妥当性についても評価することができた．

謝辞：本研究は，公益財団法人TAKEUCHI育英奨学会の助成金制度による援助を受けた．記して，謝意を表す．
図-12 𝜇𝜇𝑙𝑙 =10[Pa ∙ s] (case2 初期値)での VOF 関数の遷移34

図-13 𝜇𝜇𝑙𝑙 =50[Pa ∙ s] (case3 初期値)での VOF 関数の遷移

参考文献

1) 丸岡晃, 小保内啓太, 奥村弘, “移流拡散問題におけるSemi-Lagrange Galerkin 法”, ながれ 27 巻, pp.143-152,
2008.
2) 奥村弘, 丸岡晃, “自由界面問題に対する Semi-LagrangeGalerkin (SLG) 法 の 評 価 ”, 応 用 力 学 論 文 集 Vol.12,pp.155-162, 2009/8.
3) FreeFEM++, https://freefem.org/ (2022/9/5 参照).
4) 大塚厚二, 高石武史, “有限要素法で学ぶ現象と数理FreeFEM++数理思考プログラミング-”, 2015/3/5.
5) 日本計算工学会編, “有限要素法による流れのシミュレーション 第 3 版”, 丸善出版, 2017/7/30.
6) 矢部孝, 内海隆行, 尾形陽一, “CIP 法-原子から宇宙までを解くマルチスケール解法-”,森北出版, 2008/8/5.
(2022. 9. 8 受付)

35