オープン CAE 学会 技術報告
Technical Report of the Open CAE Society of Japan, OpenCAEtech2024-001

Copyright© The Open CAE Society of Japan

OpenFOAM による気液混相流解析の理解を深めるための例題集の作成
毛管現象 capillaryRise 例題の拡張
好田菜々1 中川慎二 1†1
富山県立大学

Tutorials for a better understanding of Gas-Liquid multiphase simulation with
OpenFOAM
Extension of capillary phenomena tutorial capillaryRise
Nana YOSHIDA* Shinji NAKAGAWA*†*
Toyama Prefectural University

Abstract
OpenFOAM is powerful toolkit for simulating multiphase flow. It is quite useful but difficult tounderstand and use correctly for beginners. We developed a series of tutorials which will teach beginners
how to setup cases with gas-liquid multiphase flow using VOF method. A standard tutorial capillaryRise
is modified and extended. The relation between pressure boundary conditions without hydrostaticpressure and liquid height on the outside of simulation region is discussed. Capillary phenomena with
non-wettable condition is simulated and a downward movement is reproduced.
Keywords: OpenFOAM, VOF, interFoam, tutorials, capillary phenomenon

1.

はじめに

OpenFOAM では気液混相流の解析が可能であり，現在広く使用されている．最も基本的な気液混相流解析
用ソルバ interFoam では，VOF 法により非混和性 2 流体を扱い，表面張力と壁面での接触角を考慮することができる．しかし，境界条件として圧力を直接与えずに静圧から静水圧成分を除去した値を与える必要があるなど，初心者が正しく理解して使うことが難しいことがある．
本報では OpenFOAM の interFoam ソルバの標準例題である capillaryRise 例題（毛管現象）を取り上げ，それを徐々に発展させた複数の例題を報告する．なお，計算には，OpenFOAM 7 を使用する．capillaryRise 例
題を徐々に発展させた複数の例題を実行・比較することにより，OpenFOAM の設定の物理的意味を初心者が理解しやすくなることを目指す．さらに，OpenFOAM でのポスト処理自動化方法も例題に組み込み，その使い方が学べるものとする．

2.

capillaryRise 例題について

2.1. capillaryRise 例題の概要
ここでは，OpenFOAM の標準例題である capillaryRise（毛管上昇）例題について説明する．密度  の液体の自由表面に細管を立てたとき，細管内の自由表面は表面張力  によって，周囲の自由表面よりも上昇（または下降）する．このような現象を毛管現象という．細管材料と自由表面との接触角  が 90°より小さい（ぬれ易い）場合は上昇，90°より大きい（ぬれ難い）場合には下降する．
図 1(a)に計算対象の模式図，図 1(b)に境界面の名称を示す．なお，図 1(a)では，計算領域（赤い破線で囲んだ領域）以外に周囲の液位（液体と気体との界面の位置）なども表示した．本報では，3.3 節で説明する
濡れにくい場合の毛管現象にも対応できるよう，流体を水銀に変更する例題も作成した．液体が必ずしも水ではないことから，水位ではなく「液位」という言葉を使用する．capillaryRise 例題では，平行な 2 枚の板を間隔 W = 1 mm で立てた状態を想定している．解析領域（幅 W×高さ H ×奥行きZ）は 1×20×1 mm3 の直方体領域であり，奥行方向には同じ状態が続く 2 次元問題である．このとき，重力による力と表面張力による力との釣り合いから，平板間液面の周囲液面からの理論的な上昇量 hth は次式から求められる．
† E-mail address of corresponding author: nakagawa-lab-pub@ml.pu-toyama.ac.jp(Received December 26, 2023; accepted April 5, 2024; published April 15, 2024)
1

オープン CAE 学会 技術報告
Technical Report of the Open CAE Society of Japan, OpenCAEtech2024-001

Copyright© The Open CAE Society of Japanℎth = 2𝜎cos𝜃⁄𝜌𝑊𝑔(1)

標準例題では，液体を水として  = 1000 kg/m3，重力加速度 g = 10 m/s2,  = 0.070711 N/m,  = 45 °としていることから，上昇する液面の高さは hth = 0.010 m となる．なお，気体は空気であり， = 1.0 kg/m3，動粘性係数ν = 1.48×10-5 m2/s としている．

(a) schematic overview(b) boundaries
Fig. 1 Schematic diagrams of OpenFOAM standard tutorial capillaryRise.
表 1 に境界条件を示す．流入口（inlet）では，圧力 p_rgh は 0 kg/m/s2 で固定されている．ここで指定す
る圧力によって周囲の液面高さ（あるいは，inlet 面が周囲液面から沈み込んだ深さ）を指定することになる．
流入口（inlet）では，流体は自由に出入りできるが，alpha.water の VOF 値が 1 に指定されているため，流入する流体の種類は水となる．
計算領域を，幅・高さ・奥行き方向にそれぞれ ( 20, 400, 1 ) に等分割してメッシュを作成する．メッシュ
は blockMesh により生成する．
Table 1 Boundary conditions of standard tutorial capillaryRise.
boundary
inletatmosphere

U

p_rgh

type pressureInletOutletVelocity;value uniform (0 0 0);

type fixedValue;
value uniform 0;

type inletOutlet;
value uniform 1;
inletValue uniform 1;

type pressureInletOutletVelocity;value uniform (0 0 0);

type fixedValue;
value uniform 0;

type zeroGradient;type noSlip;type
fixedFluxPressure;

type
constantAlphaContactAngle;
theta0 45;
limit gradient;value uniform 0;

type empty;

type empty;

type empty;

walls

frontAndBackalpha.water

図 2 に，OpenFOAM の標準例題 capillaryRise の結果の一例（時刻 0.5 秒）を示す．図 2(a)は計算領域全体
での体積分率 alpha.water の 2 次元分布である．水の体積分率（VOF 値；alpha.water）が 1（赤）は水，0（青）は空気であることを示す．界面の形状をわかりやすくするために，一部を拡大した図を併記する．接触角を 45°としているため，気液界面は下に凸な曲線となる．図 2(b)のグラフは，計算領域中心での高さ方向の水の体積分率分布であり，Y = 9.2 mm 付近に気液界面が存在することから，表面張力による周囲液面からの上昇量 h = 9.2 mm となる．理論上の上昇量は hth = 10 mm であることから，理論上の液面の Y 座標 Yth= 10 mm となる．理論値に比べると，シミュレーション結果では表面張力による界面の変位が小さい．
図 2(c)は，計算領域中心における圧力 p および静水圧成分を除去した圧力 p_rgh の高さ方向分布である．

2

オープン CAE 学会 技術報告
Technical Report of the Open CAE Society of Japan, OpenCAEtech2024-001

Copyright© The Open CAE Society of Japan

下に凸な気液界面が存在するため，上部の空気に対して，界面を通した液内の圧力は低下する．その低下した圧力は，高さが小さくなるほど，静水圧の効果によって直線的に増加していき，計算領域底面では 0
kg/m/s2 となる．静水圧を除去した圧力も，境界条件で指定しているため計算領域底面では 0 kg/m/s2 である．
計算領域内の液位は一定とならず，上下振動が継続する．その様子を図 3 に示す．これは，両壁面および壁面間を 4 等分する位置での液位の経時変化を表す．

(a) alpha.water distribution(b) alpha.water profile(c) pressure profiles
Fig. 2 Simulation result of standard tutorial capillaryRise; profiles at x = 0.50mm (center).

Fig. 3 Simulation result of standard tutorial capillaryRise; history of liquid height at various X locations.
2.2. capillaryRise 例題の課題
ある程度の時間が経過後には定常状態に達するはずであるが，図 3 では液面が振動を続ける様子が確認された．先行研究[1]では，OpenFOAM において液面が振動を続けることが報告されている．この振動は時間微分クラス ddtScheme の fvcDdtPhiCoeff メソッドを改良する[2]ことで低減されることを確認した．しかし，対象が変わると不安定になることがあるため，ソースコードへの改良が取り消された[3]．この挙動の問題については，今後の改良によって修正が望まれ，本稿では取り扱わない．
図 1(a)では，計算領域（赤い破線で囲んだ領域）以外に周囲の液位などを表示した．しかし，計算領域外部に存在する液位については例題には明示されていないため，初心者はこれらの状況を明確に意識していな
いことがある．あるいは，なぜ周囲の液位が計算領域底部にあるのか説明できないことがある．本稿では周囲液位と境界条件との関係を明確にすることを主な目的とする．
周囲液位と境界条件との関係を明確にすることは，ぬれ性の低い（ぬれ難い）条件でのシミュレーションを実現するために欠くことができない．しかし，capillaryRise 例題の物性値と接触角だけを変更しても，ぬれ難い場合の毛細管現象を再現することはできない．そこで，周囲液位と境界条件との関係を明確にした後に，物性値を水から水銀に変更して例題を作成することとした．
p と p_rgh の関係性や水位の変化を確認したい場合に，図 2，3 のようにグラフ化することでより数値が明確になり視覚的に結果を確認しやすい．しかし，capillaryRise 例題にはグラフを自動作成できるスクリプトは用意されていない．効率的に結果を確認できるようにするため，このようなポスト処理を自動化すること3

オープン CAE 学会 技術報告
Technical Report of the Open CAE Society of Japan, OpenCAEtech2024-001

Copyright© The Open CAE Society of Japanとした．このような自動化は OpenFOAM を活用するためには重要であるが，初心者には理解が難しい．複数のポスト処理方法を使用する方法を例題に含めることで，初心者の学習に役立つ例題にした．

3.

capillaryRise 例題を発展させた新たな例題の作成

新たに作成した例題の概略を表 2 に示す．ポスト処理を自動化するという機能の追加と，計算領域外の液位を境界条件でどのように設定するかを理解する例題などを作成した．
3.1. ポスト処理の自動化
複数の例題を実行しながら結果を比較し，OpenFOAM の操作を理解するためには，効率的に結果をチェックする仕組みをつくることが望まれる．そのため，結果の確認に役立つ情報の抽出とグラフ化を実行する仕組みをつくった．2.1 節で示したグラフは，この仕組みによって生成したものである．
利用した OpenFOAM の機能（function objects）は，指定位置での分布情報を取出す sampling，計算領域内の変数の体積積分を実行する volIntegrate（volFieldValue の中の操作の 1 つ）
，計算領域内の気液界面の位置を取出す interfaceHeight である．
これらの機能により，計算結果からグラフ化に必要なデータファイルが生成される．そのファイルを
gnuplot によりグラフ化するスクリプトを作成した．付録 A にポスト処理自動化のための設定ファイルを示す．
Table 2 List of developed tutorials.
classification

summary(A)automatic post-process

Average liquid height in domain is automatically obtained andits history is shown as a graph. User can learn how to use
functionObjects.

(B)symmetry boundary condition

Symmetric boundary condition is used to reduce computational
time.

(C)

mesh

Effect of mesh number and coarseness is examined.

(D)liquid height on the outside of domain

Boundary conditions are changed to manipulate a liquid heighton the outside of a simulation region.
This is needed when non-wetting condition is considered.

(E)

coordinate and p_rgh condition

The value of p_rgh at the boundary depends on the definition of
coordinate.

(F)

contact angle

Contact angle is changed from wetting condition to non-wetting
condition. Liquid height will be compared with the theoreticalvalue.

3.2. 座標と p_rgh の関係
2.1 節において，流入口（inlet）で圧力 p_rgh を 0 kg/m/s2 に指定したことで，周囲の液面高さ（あるいは，inlet 面が周囲液面から沈み込んだ深さ）を指定することになると説明した．ここで注意が必要なのは，
0 kg/m/s2 という値だけで周囲の液位が決まるのではなく，座標の取り方に依存することである．座標の取り方と p_rgh の関係を確認するための例題を作成した．本節で説明する例題は表 2(D)である．この例題では，周囲の高さの位置を知るために，接触角を 90°に設定した．接触角を 90°にすることで，表面張力が働いても，平行平板間の液の高さは周囲と同じになる．そのため，計算領域内の液位から，周囲流体の液位を知ることができる．
標準例題と同じ計算領域の大きさに対して，Y 座標の取り方が異なる場合を考える．図 4 に 3 種の座標の定義を示す．いずれの場合も周囲の液面の高さが計算領域の中央に位置すると考えるが，その Y 座標 Yref は異なる．標準例題では，計算領域最下部（流入境界）の Y 座標が 0 であった．図 4 の(1)から(3)では，物理的な条件は同じであるが，Y 座標の原点がことなる．このように座標を変更した場合，標準例題で与えた入口での境界条件（p_rgh = 0 kg/m/s2）で計算を実行すると，計算が破綻する．それは，周囲の液位が正しく4オープン CAE 学会 技術報告
Technical Report of the Open CAE Society of Japan, OpenCAEtech2024-001

Copyright© The Open CAE Society of Japan

指定できないためである．
周囲液位の指定を正確にするためには，p_rgh の定義について明確にする必要がある．p_rgh は式(2)で求めることができる．
(2)

𝑝_𝑟𝑔ℎ = 𝑝 − 𝜌𝑔ℎh は液面からの重力方向の距離に相当する．今回の例題では，重力の向きは Y 軸の負の方向としている．
図 4 の(a)から(c)の 3 つの場合について，h として与える値と，そこから計算した p_rgh を表 3 に示す．式(2)より，計算領域最下部での圧力は計算領域外での液面位置 Yref とそこでの圧力 pref とから，
(3)𝑝_𝑟𝑔ℎ = 𝑝ref − 𝜌𝑔𝑌refで求めることができる．圧力は大気圧を基準として考えてよく，pref = 0 kg/m/s2 としている．この例題では，
3.3 節で説明する濡れにくい場合の毛管現象にも対応できるよう，流体を水銀に変更し，密度を水銀の値  =13546 kg/m3 とした．このように境界条件を設定することで，正しく計算が実行できることを確認するための例題を作成した．
図 5 に計算結果の一例として，高さ方向の圧力分布を示す．いずれの場合も，計算領域上部の空気部分では圧力が 0 kg/m/s2 となる．接触角を 0°としているため，気液界面は湾曲せずに水平となる．そのために，
液面が曲率を持つ標準例題とは異なり，界面での急激な圧力変化は生じない．気液界面の位置（計算領域の中央の高さ）より下では，圧力 p は深さと共に線形に増加し，計算領域最下面では p = ρgh = 1354.6 kg/m/s2となる．この値は，計算領域外部に仮定した液面の位置から想定される値と等しい．
図 5 に示した静水圧を除去した圧力 p_rgh は，座標の取り方によって値が異なる．基準としたい液面の位
置に Y 座標の原点を合わせた(b)の場合には，標準例題と同様に，液体内部で p_rgh が 0 kg/m/s2 となる．しかし，そうでない場合には，液体内部で一様な値をとるものの，その値は 0 にはならないことに注意が必要
である．

Fig. 4 Various coordinate system with same simulation object. Yin and Yatm are the Y-coordinates of inlet andatmosphere boundary.
Table 3 Height used in boundary value calculation and resulting p_rgh at inlet boundary.
(a)𝑌ref [m]p_rgh = 𝑝ref − 𝜌𝑔𝑌ref[kg/m/s2]

(b)

(c)

-5.0×10-30

15.0×10-3-677.302031.95

オープン CAE 学会 技術報告
Technical Report of the Open CAE Society of Japan, OpenCAEtech2024-001

Copyright© The Open CAE Society of Japan

Fig. 5 p and p_rgh profiles with various coordinate system.
3.3. ぬれにくい場合の毛管現象（capillarySink）
大気中に置かれた水銀表面に平行平板を立てる場合を考える．板と水銀表面との接触角が 135°であり90°より大きいとすると，平板間の液面は周囲よりも低くなる．ここでは，このような状態を再現する例題を作成した．本節で説明する例題は表 2(F)である． = 13546 kg/m3， = 0.4821 N/m とした．このような問題を解くためには，前節で説明したように，計算領域外の液位を任意に指定することが必要となる．
座標の取り方は標準例題と同じ（inlet 境界で Y = 0 mm）とし，周囲の液位を Yref = 10 mm とした場合の結果を図 6 と図 7 に示す．初期の液位は 10 mm としたが，大きな接触角によって液位が下降することが確認できる．
理論的な平板間液面の周囲液面からの上昇量 hth は式(1)で計算でき，この例題では負の値となる．これは，周囲の液面高さの位置 Yref から hth 下がった位置に液面があることを示す．よって，Yth = Yref + hth = 0.010 –0.00503 = 0.00497 m が理論上の液面の Y 座標 Yth となる．図 7 の液位は変動を続け，その値 Y が理論値 Yth =5 mm よりも大きい時間が多い．Y が大きいことは，表面張力によって界面が移動する距離 h の大きさが小さいことを意味する．（Y = Yref + h であり，h は負の値となる．）2.1 節の標準例題でも同様に表面張力による移動量が小さい結果であった．この問題については，OpenFOAM での界面の表現方法などが要因であると推測している．

(a) alpha distribution and profile(b) pressure profiles
Fig. 6 Simulation result of tutorial capillarySink.

6

オープン CAE 学会 技術報告
Technical Report of the Open CAE Society of Japan, OpenCAEtech2024-001

Copyright© The Open CAE Society of Japan

Fig. 7 Simulation result of capillarySink. Average liquid height is obtained from liquid volume in the region by dividing
by bottom area.

4.

まとめ

OpenFOAM の代表的なソルバ interFoam に付属する例題 capillaryRise を対象として，それを徐々に発展させた複数の例題を作成した．境界条件として圧力を直接与えずに静圧から静水圧成分を除去した値を与える必要があるなど，初心者には理解が難しい境界値設定に着目し，その詳細を理解するための例題や，ガラス細管内の水銀の動きを模擬する例題を作成した．他にも，OpenFOAM でのポスト処理自動化方法も例題に組
み込み，その使い方が学べる例題を作成した．
本報で紹介した以外にも，対称性を考慮した計算モデルの作成方法や，メッシュの生成方法を変えた例題も作成した．これらの例題をまとめて公開する[4]．

参考文献[1] 稲田翔也, 中川慎二, 西田樹生, 清家美帆, 毛細管現象による OpenFOAM（interFoam 系ソルバ）の検証. オープン CAE シンポジウム 2018, A35, 2018.
[2] OpenFOAM/OpenFOAM-6, ddtScheme::fvcDstPhiCoeff: Improved formulation providing betterstability/accuracy balance, GitHub, https://github.com/OpenFOAM/OpenFOAM6/commit/da787200a6b208cf3fc4dfaa48b4438a3a038c1b, (accessed 2024-02-27).
[3] OpenFOAM/OpenFOAM-dev, ddtScheme::fvcDdtPhiCoeff: Reverted to previous flux-normalised scheme,
https://github.com/OpenFOAM/OpenFOAM-dev/commit/b162c4818b9bbdb824db9e75baffef6ac7c5a7a0,(accessed 2024-02-27).
[4] snaka-dev/capillaryRiseExtended, GitHub, https://github.com/snaka-dev/capillaryRiseExtended, (accessed2024-02-27).

7

オープン CAE 学会 技術報告
Technical Report of the Open CAE Society of Japan, OpenCAEtech2024-001

Copyright© The Open CAE Society of Japan

付録 A ポスト処理自動化のための設定ファイル
ポスト処理を自動化するために追加した設定ファイルを Code 1 から Code 3 に示す．
Code 1 Setting file for interfaceHeight.
1
2
3
4
5
6
7
8
9
10
11
12
13
14

type
libs
alpha
locationswriteControl
writeIntervalinterfaceHeight;("libfieldFunctionObjects.so");alpha.water;// width x 0-0.001; height 0-0.02; depth z 0-0.001((0
0 0.0005) //xMin wall
(0.00025 0 0.0005) //(0.0005 0 0.0005) //center(0.00075 0 0.0005) //(0.0010 0.0005) //xMax wall);
adjustableRunTime; //writeTime;
0.01; // this entry is needed for runTime/adjustableRunTime

Code 2 Setting file for volFieldValue.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

volIntegralAlpha{//typevolRegion; // for OpenFOAM 4
type
volFieldValue; // for OpenFOAM 6, v1806
libs ("libfieldFunctionObjects.so");enabledtrue;
writeControl
writeInterval

adjustableRunTime; //writeTime;
0.01; // this entry is needed for runTime/adjustableRunTime

log
writeFields

true; //false;
false; //true;

regionType//nameall; //cellZone;
outlet; // this entry is needed for cellZone regionType.

operation

volIntegrate;//sum, weightedSum, sumMag, average, weightedAverage;//volAverage, weightedVolAverage, volIntegrate, min, max CoV//weightField

rhoPhi;fields(alpha.water);}

Code 3 Setting file for sample.
1
2
3
4
5
6
7
8
9
10
11
12
13

type sets;
libs("libsampling.so");interpolationScheme cellPointFace;
setFormat raw;
sets
(
x05mm
{
type
axis

lineFace;xyz;

8オープン CAE 学会 技術報告
Technical Report of the Open CAE Society of Japan, OpenCAEtech2024-001
14
15
16
17
18
19
20
21
22
23
24
25

start
end

(0.0005 00.0005);(0.0005 0.02 0.0005);});
fields(Ualpha.water
p
p_rgh
);

9

Copyright© The Open CAE Society of Japan