論文

Semi-Lagrange法を用いた特性曲線有限要素法による界面を有する非圧縮性粘性流れの解析- FreeFEM++を用いた数値解析 小林 雅也1・広瀬 智史2・倉橋 貴彦2・衞藤 俊彦3
1

長岡技術科学大学 機械創造工学課程（Department of Mechanical Engineering, Nagaoka University of Technology）
2

長岡技術科学大学大学院 機械創造工学専攻（Department of Mechanical Engineering, Graduate School of Nagaoka University of Technology）
3

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

Analysis of incompressible viscous flow with interface
based on the characteristic finite element method using the semi-Lagrange method- Numerical analysis using FreeFEM++ Masaya KOBAYASHI1, Satoshi HIROSE2, Takahiko KURAHASHI2
and Toshihiko ETO3
Abstract
In this paper, we present analysis results for the two-dimensional dam break problem based on thecharacteristic finite element method using the Semi-Lagrange method. The characteristic finite element
method consists of the method of characteristic and the finite element method, and the calculation process
is divided into two phases which are advection and non–advection phases considering the Semi-Lagrange
method. It is known that a stable numerical solution can be obtained even when advection is predominant
in case that the characteristic finite element method using the Semi-Lagrange method is employed. Thefluid analysis is performed by using FreeFem++. FreeFem++ provides the “convect function” that finds theupstream point of the characteristic curve. The analysis results using the “convect function” are also shownin the appendix of this paper.
Key Words : Semi-Lagrange Method, Characteristic Finite Element Method, FreeFem++, Dam-breakproblem１．

はじめに貨物の多くは貨物自動車(トラック)を用いて目的地まで運送される．貨物自動車は走行速度，積荷の種

陸上での貨物の運送が主流となっている日本では，
1類，積載状態などが要因となり走行安定性が低下しと Semi-Lagrange 法の概要を示す．

事故が起こりやすい．特に積荷が液体の場合は走行に伴い重心位置が移動しやすいため，走行安定性がより低くなることが報告されている 1)．また，地方へ貨物を運送する場合は，勾配やカーブの急な道・特性曲線有限要素法
Navier-Stokes 方程式における時間微分項と移流項(山道など)を通ることは避けられない．これらからは Lagrange 微分の形で式(3)左辺のように物質微分

貨物自動車の横転事故は多く発生している．

項として表される．特性法から，粒子軌道の出発点を 𝑋𝑋𝑖𝑖 ，終点を 𝑋𝑋𝑖𝑖′ とし，時間刻みを∆𝑡𝑡とすると式(3)

上記を考慮し，流体を積んだ貨物自動車について，数値計算を援用した構造最適化手法により走行安定右辺のように離散化することができる．

性を向上させる製品設計をすることが望ましく，そ𝐷𝐷𝑢𝑢𝑖𝑖 𝑢𝑢𝑖𝑖 (𝑋𝑋𝑖𝑖′ , 𝑡𝑡) − 𝑢𝑢𝑖𝑖 (𝑋𝑋𝑖𝑖 , 𝑡𝑡 − ∆𝑡𝑡)
=𝐷𝐷𝐷𝐷
∆𝑡𝑡のためにはタンク内における流体の挙動を精度良く数値シミュレーションする必要がある．しかし，流れの現象において，拡散に対して移流による影響が(3)

特性法により時間方向に離散化された物質微分項をGalerkin 有限要素法により求める方法が特性曲線有卓越する場合，移流項に対して中央差分的近似を行うと解が不安定になりやすいことが知られている 2)． 限要素法である．特性曲線有限要素法により解の空
そこで本論文では，移流が卓越する場合の解の不安間的特徴に方向性はなくなり，移流項が卓越する場定性を回避するため，Semi-Lagrange 法を特性曲線合でも安定した解を求めることができる．

有 限 要 素 法 3) に 組 み 込 ん だ SLG (Semi-LagrangeGalerkin) 法 4)-5)を用いて流体の数値シミュレーショ

・Semi-Lagrange 法
Semi-Lagrange 法とは Euler 法と Lagrange 法の折ンを行い，移流による影響が卓越する場合でも安定して解が得られることを示す．数値シミュレーショ
ンはプログラミング言語 FreeFem++6)-7)

衷法であり，粒子軌跡の終点 𝑋𝑋𝑖𝑖′ を固定された計算を用いて行格子の格子点となるようにし，粒子軌道の出発点う．

𝑋𝑋𝑖𝑖 を式(4)によって求める．ここで𝑎𝑎𝑖𝑖 は移流流速である．

２．

界面を有する非定常粘性流れ場の解析

２．１

𝑋𝑋𝑖𝑖 = 𝑋𝑋𝑖𝑖′ − �𝑡𝑡−∆𝑡𝑡

非定常粘性流れ場の支配方程式を用いて反復的に解くことが必要となる．しかし反復計算を行うと計算時間の増加を招くことから，本解析では移流流速を2次の Adams-Bashforth 法を用いここで，𝑢𝑢𝑖𝑖 は流速，𝑝𝑝 は圧力，𝜌𝜌 は密度，𝜇𝜇 は粘度，𝑔𝑔 は重力加速度，𝛿𝛿𝑖𝑖2 はクロネッカーのデルタて式(5)のように近似し，方程式を線形化した．

である．

２．２

(4)

厳密に解を求めるためにはNewton-Raphson 法など

支配方程式は式(1)，式(2)のように表される．

𝑢𝑢𝑖𝑖,𝑖𝑖 = 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � in Ω

𝑎𝑎𝑖𝑖 (𝑥𝑥𝑖𝑖 , 𝑡𝑡)𝑑𝑑𝑑𝑑

Navier-Stokes 方程式は移流項が非線形であるため，非 定 常 非圧 縮 粘性 流 れ場 の 支 配方 程 式と してNavier-Stokes 方程式，連続の式を考える．これらの𝜌𝜌𝑢𝑢̇𝚤𝚤 + 𝜌𝜌𝑢𝑢𝑗𝑗 𝑢𝑢𝑖𝑖,𝑗𝑗 + 𝑝𝑝,𝑖𝑖 − 𝜇𝜇𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗 = −𝜌𝜌𝜌𝜌𝛿𝛿𝑖𝑖2
in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 �in Ω𝑡𝑡

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

(1)

(5)

本解析では式(4)右辺第2項の積分は積分点を4として粒子の出発点を求めた（図-1）．Semi-Lagrange(2)法を用いて解析を行う場合，粒子の出発点 𝑋𝑋𝑖𝑖 は必ずしも計算格子の格子点と重なるとは限らない．そのため空間補間によって粒子の出発点における物理SLG法による支配方程式の離散化量を求める必要がある．

流れ場の解析において移流項が卓越する場合，解の空間的特徴に方向性が現れ中央差分で近似を行う
Galerkin 有限要素法では解が不安定になりやすい．
そこで本研究では Semi-Lagrange 法を特性曲線有限要素法に組み込んだ Semi-Lagrange Galerkin (SLG) 法を用いて解析を行った．以下に特性曲線有限要素法
2𝑢𝑢𝑖𝑖𝑛𝑛+1 − 𝑢𝑢�𝚤𝚤𝑛𝑛
𝑛𝑛+1
� �𝜌𝜌 �� 𝑤𝑤𝑖𝑖 − 𝑝𝑝𝑛𝑛+1 𝑤𝑤𝑖𝑖,𝑖𝑖 + 𝜇𝜇�𝑢𝑢𝑖𝑖,𝑗𝑗
𝑤𝑤𝑖𝑖,𝑗𝑗 �∆𝑡𝑡
Ωe+ 𝜌𝜌𝜌𝜌𝛿𝛿𝑖𝑖2 𝑤𝑤𝑖𝑖 � 𝑑𝑑Ωe𝑛𝑛+1�𝑛𝑛𝑗𝑗 𝑑𝑑Γe = 0− � 𝑤𝑤𝑖𝑖 �−𝑝𝑝𝑛𝑛+1 + 𝜇𝜇𝑢𝑢𝑖𝑖,𝑗𝑗
Γe𝑛𝑛+1� 𝑞𝑞𝑢𝑢𝑖𝑖,𝑖𝑖𝑑𝑑Ω𝑒𝑒 = 0Ωe

界における外向きの単位法線ベクトルを表す．
Galerkin 法に基づくと各要素において，𝑤𝑤𝑖𝑖 , 𝑞𝑞 , 𝑢𝑢𝑖𝑖 ,本研究では式(6)で示す3角形1次要素による内挿補𝑢𝑢�𝚤𝚤 , 𝑝𝑝 は式(13)のように近似される．

間により近似して求めた．

𝑒𝑒𝑤𝑤𝑖𝑖ℎ = 𝑁𝑁𝛼𝛼𝑒𝑒 𝑤𝑤𝑖𝑖𝑖𝑖
𝑞𝑞 ℎ = 𝑁𝑁𝛼𝛼𝑒𝑒 𝑞𝑞𝛼𝛼𝑒𝑒𝑒𝑒𝑢𝑢𝑖𝑖ℎ = 𝑁𝑁𝛼𝛼𝑒𝑒 𝑢𝑢𝑖𝑖𝑖𝑖

𝑢𝑢(𝑋𝑋1 , 𝑋𝑋2 ) = 𝑁𝑁1 (𝑋𝑋1 , 𝑋𝑋2 )𝑢𝑢1 + 𝑁𝑁 2 (𝑋𝑋1 , 𝑋𝑋2 )𝑢𝑢2+ 𝑁𝑁 3 (𝑋𝑋1 , 𝑋𝑋2 )𝑢𝑢3

𝑁𝑁 𝛼𝛼 (𝑥𝑥1 , 𝑥𝑥2 ) = 𝑎𝑎𝛼𝛼 + 𝑏𝑏 𝛼𝛼 𝑥𝑥1 + 𝑐𝑐 𝛼𝛼 𝑥𝑥2𝑒𝑒𝑢𝑢�𝚤𝚤ℎ = 𝑁𝑁𝑒𝑒𝛼𝛼 𝑢𝑢�𝚤𝚤𝚤𝚤(𝛼𝛼 = 1 , 2 , 3)1∆𝑒𝑒 = �(𝑥𝑥11 − 𝑥𝑥13 )(𝑥𝑥22 − 𝑥𝑥23 )2− (𝑥𝑥23 − 𝑥𝑥21 )(𝑥𝑥13 − 𝑥𝑥12 )�

ここで 𝑁𝑁𝛼𝛼𝑒𝑒 は要素 𝑒𝑒 における形状関数である．

式(13)を用いることで空間方向に離散化された重み付き残差式を得ることができ，これらすべてを重ね合わせた連立方程式を解くことで解を得ることができる．

(6)

２．３

位置をVOF関数 Φ と呼ばれるスカラー関数により表現する．図-2 のように液体であれば 1 ，気体で出発点における物理量 𝑢𝑢�𝚤𝚤𝑛𝑛 を用いると，時間方向に離散化された支配方程式 (式(7),式(8)) を得る．
𝜌𝜌 �

∆𝑡𝑡

𝑛𝑛+1
� + 𝑝𝑝,𝑖𝑖𝑛𝑛+1 − 𝜇𝜇𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗+ 𝜌𝜌𝜌𝜌𝛿𝛿𝑖𝑖2 = 0

𝑛𝑛+1𝑢𝑢𝑖𝑖,𝑖𝑖=0

VOF法による界面の決定

界面はVOF法を用いて決定した．VOF法は界面Semi-Lagrange 法によって求められた，粒子軌道の𝑢𝑢𝑖𝑖𝑛𝑛+1 − 𝑢𝑢�𝚤𝚤𝑛𝑛

(13)

𝑝𝑝ℎ = 𝑁𝑁𝛼𝛼𝑒𝑒 𝑝𝑝𝛼𝛼𝑒𝑒𝑎𝑎𝛼𝛼 =(𝛼𝛼, 𝛽𝛽, 𝛾𝛾) = (1,2,3), (2,3,1), (3,1,2)(12)

ここで， Ω は境界 Γ で囲まれた解析領域， 𝑛𝑛𝑗𝑗 は境

図-1 数値積分イメージ

1�𝑥𝑥 𝛽𝛽 𝑦𝑦 𝛾𝛾 − 𝑦𝑦 𝛽𝛽 𝑥𝑥 𝛾𝛾 �2∆𝑒𝑒1
�𝑦𝑦 𝛽𝛽 − 𝑦𝑦 𝛾𝛾 �𝑏𝑏 𝛼𝛼 =2∆𝑒𝑒
1�𝑥𝑥 𝛾𝛾 − 𝑥𝑥 𝛽𝛽 �𝑐𝑐 𝛼𝛼 =2∆𝑒𝑒(11)

あれば 0 となるようにVOF関数を定義する．

(7)
(8)

有限要素法を適用し，Navier-Stokes 方程式にかかる重み関数を 𝑤𝑤𝑖𝑖 ，連続の式にかかる重み関数を 𝑞𝑞 とすると重み付き残差式は式(9)，式(10)で表される．
� 𝑤𝑤𝑖𝑖 �𝜌𝜌 �Ωe𝑢𝑢𝑖𝑖𝑛𝑛+1 − 𝑢𝑢�𝚤𝚤𝑛𝑛
𝑛𝑛+1+ 𝜌𝜌𝜌𝜌𝛿𝛿𝑖𝑖2 � 𝑑𝑑Ωe� + 𝑝𝑝,𝑖𝑖𝑛𝑛+1 − 𝜇𝜇𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗∆𝑡𝑡

図-2 VOF関数による流体の判別(9)=0𝑛𝑛+1� 𝑞𝑞𝑢𝑢𝑖𝑖,𝑖𝑖𝑑𝑑Ωe = 0Ωe

領域内の各節点における液体の密度 𝜌𝜌𝑙𝑙 と粘度 𝜇𝜇𝑙𝑙 ，気体の密度 𝜌𝜌𝑔𝑔 と粘度 𝜇𝜇𝑔𝑔 はVOF関数を用いて式(14)，(10)式(15)のように表される．

𝜌𝜌 = 𝜌𝜌𝑙𝑙 Φ + 𝜌𝜌𝑔𝑔 (1 − Φ)式(8)，式(9)の弱形式は式(11)，式(12)で表される．

3𝜇𝜇 = 𝜇𝜇𝑙𝑙 Φ + 𝜇𝜇𝑔𝑔 (1 − Φ)

(14)
(15)

VOF 関数の支配方程式は移流方程式であり，式(16)Φ̇ + 𝑎𝑎𝑖𝑖 Φ,𝑖𝑖 = 02 −1/22 22‖𝑎𝑎𝑖𝑖 ‖� +�� �ℎ𝑒𝑒∆𝑡𝑡で表される．

𝑒𝑒𝜏𝜏𝑚𝑚= ��(16)(19)ℎ𝑒𝑒 は要素サイズであり本研究では各要素の 2 倍の面本解析では式(16)を SUPG 法に基づく安定化有限要積と同じ面積の正方形の 1 辺の長さとした．

素法 2)と Crank-Nicolson 法を用いて解析を行った．
以下にそれぞれの概要を示す．

・Crank-Nicolson 法
時間ステップ 𝑛𝑛 + 𝜃𝜃 (0 ≤ 𝜃𝜃 ≤ 1) を満たす 𝜃𝜃 にお・SUPG 法に基づく安定化有限要素法

いて Φ𝑛𝑛 ，Φ𝑛𝑛+1 を 𝑛𝑛 + 𝜃𝜃 まわりでテイラー展開し，移流項が卓越する場合，解が不安定になりやすい差をとると式(20)が得られる．
𝑑𝑑𝑢𝑢𝑛𝑛+1 − 𝑢𝑢𝑛𝑛 = ∆t 𝑢𝑢𝑛𝑛+𝜃𝜃
𝑑𝑑𝑑𝑑(∆𝑡𝑡)2 𝑑𝑑 2 𝑛𝑛+𝜃𝜃
22𝑢𝑢+{𝜃𝜃 − (1 − 𝜃𝜃) }2 𝑑𝑑𝑡𝑡 233(∆𝑡𝑡) 𝑑𝑑 𝑛𝑛+𝜃𝜃
𝑢𝑢−{𝜃𝜃 3 − (1 − 𝜃𝜃)3 }6 𝑑𝑑𝑡𝑡 3+⋯

ことは先ほど述べた．この問題を解決する手法として主に特性法と風上法が提案されている．特性法に基づいた解の不安定性の回避については 2.2 節で述べた通りである．今回は風上法に基づいて解の不安定性を回避する．SUPG 法に基づく安定化有限要素法とは，図-3 に示すように中央差分で近似される
Galerkin 法の重み関数に加えて移流項に対して後退(20)式(20)を変形すると式(21)が得られる．

差分で近似したときに導出される式(18)で表される
人工拡散項 𝜉𝜉𝑖𝑖ℎ を付加した SUPG 法の重み関数を用

𝑑𝑑 𝑛𝑛+𝜃𝜃 𝑢𝑢𝑛𝑛+1 − 𝑢𝑢𝑛𝑛𝑢𝑢=∆𝑡𝑡
𝑑𝑑𝑑𝑑∆𝑡𝑡 𝑑𝑑 2 𝑛𝑛+𝜃𝜃
22𝑢𝑢−{𝜃𝜃 − (1 − 𝜃𝜃) }2 𝑑𝑑𝑡𝑡 223(∆𝑡𝑡) 𝑑𝑑 𝑛𝑛+𝜃𝜃
𝑢𝑢+{𝜃𝜃 3 − (1 − 𝜃𝜃)3 }6 𝑑𝑑𝑡𝑡 3+⋯いることで，解の不安定性を解消する方法である．

(21)

式(21)において 𝜃𝜃 = 1/2 とするとき Crank-Nicolson法と呼ばれ時間 2 次精度となる．このとき， 𝑢𝑢 𝑛𝑛+1/2は式(22)で表される．

𝑢𝑢𝑛𝑛+1/2 =𝑢𝑢𝑛𝑛+1 + 𝑢𝑢𝑛𝑛
2

(22)

VOF 関数の支配方程式である移流方程式にかかる重み関数を 𝑤𝑤� とすると，SUPG 法に基づく安定化有限要素法と Crank-Nicolson 法によって離散化された重み付き残差式は式(23)で表される．

𝜕𝜕ΦhΦℎ 𝑛𝑛+1 − Φℎ 𝑛𝑛
� 𝑤𝑤 ℎ �+ 𝑎𝑎𝑖𝑖ℎ� 𝑑𝑑Ω𝑒𝑒∆𝑡𝑡𝜕𝜕𝑥𝑥𝑖𝑖Ωe𝜕𝜕𝑤𝑤 ℎ Φℎ 𝑛𝑛+1 − Φℎ 𝑛𝑛
𝜕𝜕Φh+ � 𝜏𝜏𝑚𝑚 �𝑎𝑎𝑖𝑖ℎ+ 𝑎𝑎𝑖𝑖ℎ��� 𝑑𝑑Ω𝑒𝑒 = 0𝜕𝜕𝑥𝑥∆𝑡𝑡
𝜕𝜕𝑥𝑥𝑖𝑖𝑖𝑖
Ωe

図-3 重み関数(1 次元 2 節点 1 次要素)

�𝚤𝚤ℎ = 𝑤𝑤 ℎ + 𝜉𝜉 ℎ𝑤𝑤𝑖𝑖𝑖𝑖

𝜕𝜕𝑤𝑤𝑖𝑖ℎ𝜉𝜉𝑖𝑖ℎ = 𝜏𝜏𝑚𝑚 𝑎𝑎𝑗𝑗ℎ𝜕𝜕𝑥𝑥𝑗𝑗(17)

𝑤𝑤 ℎ = 𝑁𝑁𝛼𝛼𝑒𝑒 𝑤𝑤𝛼𝛼𝑒𝑒Φℎ = 𝑁𝑁𝛼𝛼𝑒𝑒 Φ𝛼𝛼𝑒𝑒

(18)

(23)

(24)

これらすべてを重ね合わせた連立方程式を解くことで解を得ることができる．

ここで 𝜏𝜏𝑚𝑚 は時間の次元をもつ安定化パラメータであり一般的に式(19)で表され，本研究でも式(19)をSUPG 法に基づく安定化有限要素法では移流流速が用いた．

急激に変化する界面近傍では，オーバーシュート及4びアンダーシュートの挙動が見られることが知られ３．２

解析結果ている．そのため本研究では VOF 関数が 1 を超え解析結果を図-5 から図-10 に示す．結果よりダムる場合は VOF 関数を 1 とし，同様に 0 を下回る場合ブレイク時における水塊が壊れていく過程を適切には 0 とする処理を施した．

再現できていることがわかる．

３．

数値解析例

2 次元ダムブレイク問題において計算精度の検討を行った．
３．１

解析条件

解析モデルは図-4 に示す通り 1 辺 𝐿𝐿 = 0.6 [m] の剛体タンクを設定した．気液二相流体には水と空気の物性値を与え，初期条件として静止状態の幅 𝐿𝐿/4，高さ 𝐿𝐿/2 の水柱を与えた．流速と圧力，VOF 関数

は共に FreeFem++ 上の関数「buildmesh」を用いて 1
辺を 160 分割するよう 3 角形 1 次要素を生成した．
その他解析条件は表-1 に示す通りである．

また，すべての境界に対して slip 境界条件を与えた．

図-5 界面関数の分布 (経過時間 𝑡𝑡 = 0 s)

図-6 界面関数の分布 (経過時間 𝑡𝑡 = 0.1 s)
図-4 解析モデル
表-1 数値解析条件

節点数
要素数
時間ステップ 𝑡𝑡𝑓𝑓 [s]
時間刻み幅 ∆𝑡𝑡 [s]液体の密度 𝜌𝜌𝑙𝑙 [kg/m3]
液体の粘度 𝜇𝜇𝑙𝑙 [Pa∙ s]
気体の密度 𝜌𝜌𝑔𝑔 [kg/m3]
気体の粘度 𝜇𝜇𝑔𝑔 [Pa∙ s]
重力加速度 𝑔𝑔 [m/s2]

29346
58051
0.5
0.001
1.0 × 103
1.0 × 10−3
1.2
1.5 × 10−5
9.85

図-7 界面関数の分布 (経過時間 𝑡𝑡 = 0.2 s)

場合においても Semi-Lagrange Galerkin 法を用いて解析を行うことで安定して解が得られることが示せた．

謝辞：本研究を行うにあたり，長岡技術科学大学機械創造工学専攻2020年度修了の尾関優汰氏から多くの助言を受けた．記して謝意を表する．

付録
FreeFem++には上流点の位置とその位置における物理量を求める convect 関数が備わっている．図-11から図-16に convect 関数を用いた解析結果を示す．

図-8 界面関数の分布 (経過時間 𝑡𝑡 = 0.3 s)

convect 関数を用いた場合も適切に解析が行えていることを確認できる．

図-9 界面関数の分布 (経過時間 𝑡𝑡 = 0.4 s)
図-11 界面関数の分布 (経過時間 𝑡𝑡 = 0 s)

図-10 界面関数の分布 (経過時間 𝑡𝑡 = 0.5 s)

４. 結論

図-12 界面関数の分布 (経過時間 𝑡𝑡 = 0.1 s)

2 次元ダムブレイク問題において移流が卓越する
6図-13 界面関数の分布 (経過時間 𝑡𝑡 = 0.2 s)

図-16 界面関数の分布 (経過時間 𝑡𝑡 = 0.5 s)

参考文献1) 交通事故総合分析センター 事業用自動車事故調査委員会, “事業用自動車事故調査報告書”, 2015/12/7.
2) 日本計算工学会編, “有限要素法による流れのシミュレーション 第 3 版”, 丸善出版, 2017/7/30.
3) 田端正久, “ガレルキン・特性曲線有限要素法-理論と応用-”, 数理解析研究所講究録, 第 1791 巻, pp.76-86, 2012.
4) 丸岡晃, 小保内啓太, 奥村弘, “移流拡散問題におけるSemi-Lagrange Galerkin 法”, ながれ 27 巻, pp.143-152, 2008.
5) 奥村弘, 丸岡晃, “自由界面問題に対する Semi-Lagrange
Galerkin (SLG) 法 の 評 価 ”, 応 用 力 学 論 文 集 Vol.12,pp.155-162, 2009/8.
図-14 界面関数の分布 (経過時間 𝑡𝑡 = 0.3 s)

6) FreeFem++, https://freefem.org/ (2021/7/23 参照).
7) 大塚厚二, 高石武史, “有限要素法で学ぶ現象と数理FreeFem++数理思考プログラミング-”, 2015/3/5.
(2021. 8. 19 受付)

図-15 界面関数の分布 (経過時間 𝑡𝑡 = 0.4 s)

7