九 州 大 学 大 学 院 総 合 理 工 学 報 告第27巻 第2号 229-240頁 平 成 17 年 9 月

Engineering Sciences Report, Kyushu University
Vol.27, No.2 pp.229-240, SEP. 2005

カットセル法を用いた海洋モデルについて
及川 正道*1†・経塚 雄策*2
(平成17年7月29日 受理)

Ocean Model using Cut-Cell Method
Masamichi OIKAWA, Yusaku KYOZUKA†

E-mail of corresponding author: oikawa@esst.kyushu-u.ac.jp

The ocean model using Cut-Cell method has been developed to simulate flow field around topographic
features with steep slope precisely. The model can describe flow field along topography including vertical flow
more naturally than an existing z-coordinate model.
Key words: Cut-Cell Method, Finite-Volume Method, Numerical Simulation, Ocean Model

1. 緒

言

より実際のものに近づき、計算精度が向上することが数値シミュレーションにおいて、急勾配で高低差の期待される。

大きい海山のような海底地形を取り扱う際には、種々モデルは、z座標を用いているMECモデル9)10)をベーの問題が生じる。σ座標系を用いたモデルでは、水平スとしており、MECモデルにCut-Cell法を用いるため

圧力勾配の計算誤差から偽の流れが現れる 。z座標系の、変数及びサブルーチンを付加することにより拡張を用いたモデルでは、海底地形を階段上に近似するこしている（以降は、MECのオリジナルモデルを単に1)とから、一般的に海底近くの流れの精度は低く、勾配

MECモデルと呼び、MECにCut-Cell法を適用したモデが大きい海底地形では局所的に強い鉛直流が現れる 。

ルをCut-Cellモデルと呼ぶことにする）。作成した海洋

航空機・自動車等の流体機械の設計等で使われる数モデルを用いてガウシアン型海山の周囲の流れを計算

値流体力学(CFD)では、任意形状の物体を取り扱えるした結果、従来のモデルを用いて計算したものに比べように様々な手法が検討されている。物体形状に適合て精度が向上していることを確認した。

2)するよう生成される格子（境界適合格子）を有限要素

本件では、第2節と第3節で本モデルのベースとなっ法や有限体積法を適用して用いる方法や、直交格子をたMECモデルと作成したCut-Cellモデルの概要につい

用 い た VCE (Virtual Cell Embedding) 法 3) 、 Immersedてそれぞれ説明し、第4節においてCut-Cellモデルを用

Boundary

法4)、CIP (Cubic-Interpolated Propagation)

法5)、

いて行った海山周りの流動場の計算結果を、MECモデ

仮想境界法 (Virtual Boundary)6)などが用いられていルを用いた際の計算結果と共に示す。最後に、第5節にる。しかしながら物体境界の取り扱いが厳密化される

おいて、本研究のまとめを行う。

ものの、境界適合格子では任意形状を扱う際の格子生成の計算コストが大きく、直交格子を用いた手法では、

2.

MEC モデルの概要

計算処理が複雑になり３次元問題を取り扱う際の計算

MECモデルは、沿岸海洋の物理環境を解くための数コストが大きくなる。そこで、格子生成の計算コスト値モデルであり、静水圧モデル、Full 3Dモデル、そしが小さく、計算処理も比較的単純な直交格子を用いたてその結合ルーチンからなる複合モデルである。本モ

Cut-Cell法 7),8) を海洋モデルに適用することにした。

デルでは、この内静水圧近モデルのみを使用した。以Cut-Cell法は、海底地形の形状を平面または曲面で近下にMECモデルの静水圧近似部分について簡単に説

似し、計算セル内において流体が占める領域にのみ有明する。

限体積法を適用し支配方程式を離散化している。これ

静水圧部分の支配方程式は、Boussinesq近似を適用しにより、海底近傍の流れが従来のz座標モデルに比してたNavie-Stokes型の運動方程式、連続の式、水温・塩分の移流・拡散式からなる。鉛直方向の運動方程式には*1 大気海洋環境システム学専攻博士課程*2 総合理工学研究院 流体環境理工学部門230

カットセル法を用いた海洋モデルについて

静水圧近似が適用されている。運動方程式及び連続の3.

式について、以下に示す。

従来のz座標系を用いた海洋モデルでは、Fig.1(a)に

1 ∂P
∂u
∂u
∂u
∂u
+u
+v
+w
= fv −
ρ 0 ∂x∂t∂x∂y∂z⎛ ∂ 2u ∂ 2u ⎞ ∂⎛ ∂u ⎞+ AM ⎜⎜ 2 + 2 ⎟⎟ + K M ⎜ ⎟∂z∂y ⎠⎝ ∂x ⎠⎝ ∂x

示すように海底地形を階段状に近似している。この場
合、海底地形近傍の流れの精度の低下は避けがたく、(1)に海底地形を平面または曲面で近似することにより、海底地形近傍の流れの精度を改善するものである。この場合、1つの計算セルの内部に流体領域と海底地形な(2)

どの固体領域が共存することになる。支配方程式の離散化は、計算セルの流体領域に対して有限体積法を適用することによって行っている。

1 ∂P
−g
ρ ∂z

(3)

∂u ∂v ∂w++=0
∂x ∂y ∂z(4)0=−

局所的に強い流れが現れることも指摘されている2)。
Cut-Cell法は、z座標系においてFig.1(b)に示すよう

1 ∂P
∂v
∂v
∂v
∂v
+u
+v
+w
= − fu −ρ 0 ∂y∂t∂x∂y∂z⎛ ∂ 2v ∂ 2v ⎞ ∂⎛ ∂v ⎞+ AM ⎜⎜ 2 + 2 ⎟⎟ +KM ⎜ ⎟∂z∂x∂y⎝ ∂x ⎠⎝⎠

Cut-Cell 法

上式で使われている変数は以下の通りである。
u: 流速（x方向）(m/s)

v: 流速（y方向）(m/s)

w: 流速（z方向）(m/s)

P: 圧力 (N/m2)

ρ: 密度 (kg/m3)(a) Step topography（ρ0は定数で代表密度を表す。）
f: コリオリパラメーター(s-1)

g: 重力加速度 (m/s2)

AM: 水平渦動粘性係数 (m2/s)

KM: 鉛直渦動粘性係数 (m2/s)

水位の時間発展の式は、連続の式を海底から海面まで積分することにより得られる。これによって解くべき未知量は水位・水平方向流速・水温・塩分となり、鉛直方向流速及び圧力については他の諸量から自動的(b) Cut Cellに決定される。これらの式をstaggered mesh (Arakawa Cgrid)上で離散化している。移流項は一次上流差分、拡Fig.1 The representation of a smoothly varying bottom

散項は二次中央差分となっている。水位と水平方向流(broken line) in a z-coordinate model.

速の時間発展の式の組み合わせには、Forward-Backwardスキームが用いられている。水温・塩分の時間発展の式はオイラーの陽解法によって解かれる。粘性及び拡散は、それぞれ渦動粘性係数及び渦動拡散係数によって表現される。水平方向の渦動粘性・拡散係数についてはリチャードソンの4/3乗則が適用され、格子間隔の4/3乗に比例するように与えられている。鉛直方向の渦動粘性・拡散係数については勾配型リチャードソン数を用いた成層化関数により求めら
れている。

3.1

有限体積法

ここで、モデルで使われている有限体積法について
以下に示すスカラー量φの保存式を例にして説明する。

r∂φ +∇⋅F = 0∂t
ここで、Fはφに関するフラックスを表す。

(5)

平成17年

九 州 大 学 大 学 院 総 合 理 工 学 報 告

第27巻

第2号

231

(5)式の両辺について体積がVである計算セル内で体積(11)式を、計算セル内で空間積分することにより積分積分を行う。第2項には、ガウスの発散定理を適用し、型の支配方程式が得られる。

体積積分を面積分に変換する。

r r∂
φdV + ∫A F ⋅ n dA = 0
∫
V
∂t

r
r
r r
∂u
∫ ∂t dV = − ∫ div(u ⊗ u )dV − ∫ div( pI )dV
V
V
V
I .J

(6)

r
1
∫ div( grad (u ))dV
Re V
r r
r r r
= − ∫ (u ⊗ u ) ⋅ n dA − ∫ ( pI ) ⋅ ndA
I .J

きの単位ベクトルを表している。(6)式は、体積内のφSの変化が計算セルの表面を通してのフラックスの出入

S

(12)

r r
1
+
grad (u ) ⋅ ndA∫
Re S

によって決定されることを意味している。
(6)式を計算セルへ適用することを考える。I番目のセル

上式(12)に、前節で示した有限体積法による離散化をについて(6)式は、以下のように表現される。

∂
φ I + ∑ FI , J AI , J = 0
∂t
J

I,J

+

ここで、nは計算セル表面の面積要素dAに垂直で外向

VI

I ,J

適用する。
(7)

r
r r
r r
r
∂u
= −∑ (u ⊗ u ) m ⋅ ndAm − ∑ Pm I ⋅ ndAm
∂t
m
m
r
r
1
+
∑ grad (u ) m ⋅ ndAm
Re mV

ここで、上式第2項はJにより識別される計算セルの表面要素に関する和である。また、φIとFI,Jは、それぞれ

φ の計算セル体積内の平均値とFの計算セル境界面内(13)の平均値であり、以下の式で表される。

φI =

1
∫ φdV
VI V

FI , J =

I

r r
1
F
⋅ n dA∫
AI , J A
I ,J

3.3
(8)

変数配置

本モデルで新規に導入する変数をTable1にまとめ、運動量方程式に関わる主要な変数の配置をFig.2に示す。

(9)
Table1以上のように、計算セルにおいて固体領域を除き変List of variables.

uc

velocity in x direction (m/s)

vcvelocity in y direction (m/s)

Asurface area of fluid part in x direction (m2)Bsurface area of fluid part in y direction (m2)3.2 支配方程式の離散化

C

surface area of fluid part in z direction (m2)有限体積法を適用することにより支配方程式の離散

V

volume of fluid part (m3)形した流体領域に対して、任意形状に適用可能な有限体積法を用いることによって、支配方程式の離散化が可能となる。

化を行う。ここでは、例として以下に示す非圧縮性流z体のNavie Stokes方程式の離散化を行う。

∂u i ∂ (u i u j )∂P1 ∂ 2 ui
+
=−
+
∂t
∂x j
∂x j Re ∂x j ∂x jyCi,j,k+1(10)

fluid partvi,j+1,k

ここで、uは速度、Pは圧力、Reはレイノルズ数を表している。i,jに関してはアインシュタインの縮約を取る
Ai,j,kものとする。

r
r
r
r r
∂u
1
= − div(u ⊗ u ) − div ( pI ) +div( grad (u ))∂tRe⎪⎧⊗ : Tensor product⎨r⎪⎩ I : Unit Tensor

uci,j,kvci,j,kui,j,k微分型の(10)式をテンソル形の発散型に書き換えると次式のようになる。

x

wi,j,k+1
Bi,j+1,k
ui+1,j,k

Ai+1,j,k

Pi,j,kvi,j,kBi,j,kwi,j,k(11)

solid partCi,j,kFig.2

A cell being cut by solid.

232

カットセル法を用いた海洋モデルについて

計算セル内の灰色の領域は固体領域を表しており、そ求められる境界面上の速度変数である。例としてuc(x)れ以外が流体領域を表している。海底地形を考慮したについては、補間係数qを用いて、以下のように表され計算セルにおける境界面積(A,B,C)と体積(V)の計算方る。

法を付録A.1に記す。

uci(,xj),k = qi − 2, j ,k uci − 2, j ,k + qi −1, j ,k uci −1, j , k

変数配置はMECモデルでは、staggered配置であるが、本モデルでは、計算セル中心に水平速度変数(uc,vc)を(15)+ q i , j ,k uci , j ,k + qi +1, j ,k uci +1, j , k

加えることによりcollocate配置としている。collocate配

本モデルは、格子間隔を変化させて用いるvariable

置では、格子毎に変数値が激しく変動する計算モード

meshの使用を想定しており、格子間隔に応じて補間係を起こしやすいとされるため、海洋モデルでは、計算数を定める必要がある。不等格子間隔の場合のQUICKの安定性が高いstaggered配置を使うことが多い。しか補間係数の求め方について付録A2を参照されたい。

し、staggered配置の場合には、有限体積法を適用する

格子間隔が一定の場合には以下のようになる。

変数が分散配置されるため、変数毎に有限体積法を適用する計算セルを設定し、個別に体積、表面積を評価しなければならず、計算が煩雑になる。本モデルでは、変数がセル中心に集中配置され計算セルを変数同士で共用でき、有限体積法の適用が容易となるcollocate配置を選択した。計算モードを抑制するためMECモデルで使われているstaggered配置された速度変数も併用する こ と と し た。 本 モ デ ル で新 し く 加 え た速 度変数(uc,vc)は運動量方程式を、MECモデルで使用される速度変数(u,v,w)は連続の式を満たすものとしている。
3.4

差分式

以下に、移流、圧力及び粘性の各項の具体的な差分式をFig.2に示した計算セルに基づき求める。

Case 1: u i , j ,k > 0
1
6
uci(,xj),k = − × uci − 2, j ,k + × uci −1, j ,k8
83
+ × uci , j ,k + 0 × uci +1, j ,k
8

(16a)

Case 2 : u i , j , k < 0

uci(,xj),k = 0 × uci − 2, j ,k +3
× uci −1, j ,k
8

(16b)

6
1
+ × uci , j ,k − × uci +1, j , k8
8
3.4.2 圧力項
圧力項は以下のように、差分表示される。

prs _ termr r
= ∑ Pm I ⋅ ndS m

3.4.1 移流項
移流項は以下のように、差分表示される。

m= [ Ai , j , k P 1

adv _ termr rr= ∑ (u ⊗ u ) m ⋅ n dS m

i + , j ,k2

− Ai −1, j , k P 1r+ ( Ai , j ,k − Ai −1, j , k ) Pp ] ⋅ im= [ Ai , j , k uc i(,xj), k u i , j , k − Ai −1, j , k uc i(−x1), j , k u i −1, j , k+ [ Bi , j , k Pr+ ( Bi , j , k − Bi , j −1, k ) Pp ] ⋅ j+ B i , j , k uc

( y)i, j ,kv i , j , k − Bi , j −1, k uc

+ C i , j , k uc

(z)i , j ,kwi , j , k − C i , j , k −1uc

( y)i , j −1, k

1i, j + ,k2

v i , j −1, k(z)i , j , k −1rwi , j , k −1 ] ⋅ i+ [ Ai , j , k vc i(,xj), k u i , j , k − Ai −1, j , k vc i(−x1), j , k u i −1, j , k+ Bi , j , k vc

( y)i, j ,kv i , j , k − Bi , j −1, k vc

( y)
i , j −1, ki − , j ,k2

− Bi , j −1, k P

1
i , j − ,k2(17)Pi-1/2,j,k, Pi+1/2,j,k, Pi,j-1/2,k, Pi,j-1/2,kは、計算セルの境界面上の圧力値を表し、これらの値は、隣り合う計算セル内v i , j −1, kr+ C i , j , k vc i(,zj), k wi , j , k − C i , j , k −1 vc i(,zj), k −1 wi , j , k −1 ] ⋅ jの値との1次の線形内挿により求められる。また、Ppwherer⎧⎪i : unit vector pointing to x direction⎨r⎪⎩ j : unit vector pointing to y direction

おり、当該計算セルと隣接する計算セルの圧力値を用は、流体と固体領域の境界面における圧力値を表して
いて1次の線形外挿により求められる。その際に当該セ(14)ルと隣接するセルの中心点を通る鉛直平面内における固体領域境界面の代表点の位置(dx,dz)のみを考慮し、移流項の差分式(14)の上添え字(x),(y)が付いた速度

この鉛直平面に垂直な方向の固体領域の位置は考慮し変数は、計算セル中心の速度変数を用いて、それぞれxないものとする。以下に、例としてx方向の圧力勾配の方向とy方向に3次精度のQUICK (quadratic upstream

求め方を、Fig.3に示す計算セルに基づいて説明する。

interpolation for convective kinematics) 補間を適用して

平成17年

九 州 大 学 大 学 院 総 合 理 工 学 報 告

dz = delz k −Ci,j,k+1
Ai,j,kPi-1,j,kPw

Ai+1,j,kdz

Pi,j,k

Pp

z

Pz dx
Ci,j,k

solid
x
Fig.3

Peb

第2号

233

(C i , j , k +1 ⋅ delz k − Vi , j , k )C i , j , k +1 − C i , j , k(21)

流体領域境界面上の圧力値Pebは、1次の線形内挿によdelzk

り次式で求められる。

Pi+1,j,k

Pe

第27巻

Peb =

Pee

(delxi ⋅ Pee + delxi +1 ⋅ Pz )delxi + delxi +1(22)

固体領域境界面上の圧力値Ppは、1次の線形外挿によdelxi

Arrangement of variables for calculation of apressuregradien term.

り求められる。

Pp = Pz − ( Pe − Pz ) ⋅x方向の圧力勾配は次式で表される。

(2dx − delxi )delxi + delxi +1(23)

ここで、dxは計算セルにおいて固体流域に覆われている境界面に相対する境界面から固体領域代表点まで

Pressure Gradient in X direction=

Ai , j , k ( Pe − Pw ) + ( Ai +1, j , k − Ai , j , k )( Peb − Pp )(18)のx方向の水平距離を表している。固体領域のx方向の幅は、計算セル内における固体領域の体積を計算セル

Vi , j , k

ここで、Vは計算セルの流体領域の体積を表している。

側面で固体領域が占める面積で割ることにより求めら

また、上式中で用いられている境界面における圧力値れる。以下の式で示すように、dxは計算セルのx方向のは以下のように求められる。

幅(delxi)から固体領域の幅を引くことにより求められx方向の2つの境界面上で、相対する領域が2つとも流る。

体となる領域（面積は、Ai,j,kで表される）における圧力値PwとPeは以下のように当該セルと隣接セルの圧力

dx = delxi −( Ai +1, j , k ⋅ delxi − Vi , j , k )値を用いて1次の線形内挿として求められる。

Pw =
Pe =

(delx i −1 ⋅ Pi , j , k + delx i ⋅ Pi −1, j , k )delx i −1 + delx i(delx i ⋅ Pi +1, j , k + delx i +1 ⋅ Pi , j , k )delx i + delx i +1

3.4.3

(24)

Ai +1, j , k − Ai , j , k

粘性項

(19a)

粘性項は以下のように、差分表示される。

(19b)

vis _ termrr
= ∑ grad (u ) m ⋅ ndS mm

x方向の2つの境界面上で、相対する領域のうちどちらかが固体となる領域における圧力値は、当該セルと隣= [ Ai , j , k grad (u ) x i + 1 , j , k − Ai −1, j , k grad (u ) x i − 1 , j , k2

2

接セルのそれぞれにおいて圧力定義点から固体領域鉛+ ( Ai , j , k − Ai −1, j , k ) grad (u ) ]直中心位置まで密度を用いて外挿して求められた2つ+ Bi , j , k grad (u ) y i , j + 1 , k − Bi , j −1, k grad (u ) y i , j − 1 , kの圧力値PzとPeeを用いて求められる。ただし、計算セ+ ( Bi , j , k − Bi , j −1, k ) grad (u )ル内では鉛直方向の圧力変化のみを考慮している。

Pz = Pi , j , k + grv ⋅ ri , j , k ⋅ (dz − 0.5 ⋅ delz k )(20a)

Pee = Pi +1, j , k + grv ⋅ ri +1, j , k ⋅ (dz − 0.5 ⋅ delz k )

(20b)

x
p

2

2

y
p

+ Ci , j , k grad (u ) z i , j , k + 1 − Ci , j , k −1 grad (u ) z i , j , k − 1
2
2
r
2
2
+ (Ci , j , k − Ci , j , k −1 ) ⋅ αuci , j , k uci , j , k + vci , j , k ] ⋅ i+ [ Ai , j , k grad (v) x i + 1 , j , k − Ai −1, j , k grad (v ) x i − 1 , j , k2

2

ここで、grvは重力加速度定数、ri,j,kは海水密度(kg/m )、+ ( Ai , j , k − Ai −1, j , k ) grad (v) xp ]dzは計算セル上部境界面から固体境界の代表点（Fig.3+ Bi , j , k grad (v) y i , j + 1 , k − Bi , j −1, k grad (v) y i , j − 1 , kの緑丸）までの距離を表している。固体領域の計算セ+ ( Bi , j , k − Bi , j −1, k ) grad (v)3

ル下面からの高さは、計算セル内における固体領域の体積を計算セル下面で固体領域が占める面積で割ることにより求められる。以下の式で示すように、dzは計
算セルの高さ(delzk)から固体領域の高さを引くことにより求められる。

2

2

y
p

+ Ci , j , k grad (v) z i , j , k + 1 − Ci , j , k −1 grad (v ) z i , j , k − 1
2
2
r
2
2
+ (Ci , j , k − Ci , j , k −1 ) ⋅ αvci , j , k uci , j , k + vci , j , k ] ⋅ j

(25)

234

カットセル法を用いた海洋モデルについて

粘性項は、上式のとおり、隣接する計算セルとの流体領域内の境界面と計算セル内の固体領域の境界面に

zt i(,nj+1) = zt i(,nj) + Δt[C i , j ,nk −1 wi(,nj),nk −1

おける水平速度の勾配とその境界面積との積の和にな+ Ai −1, j ,nk u i(−n1), j ,nk − Ai , j ,k u i(,nj),nkっている。ただし、鉛直方向の固体領域境界面（面積+ Bi , j −1,nk vi(,nj)−1, nk − Bi , j , k vi(,nj), nk ] / Ci , j , nk

(29)

はCi,j,k+1-Ci,j,k）では、海底からの摩擦力を受けるものとした。(25)式では、海底摩擦係数をαで表している。
例として、Fig.4の計算セルに基づいて水平速度成分
uがx方向の境界面を通して受ける応力を示す。

を参照）に従って、オイラーの陽解法により以下のよ
うに行われる。

Ai,j,kAi+1,j,kui,j,kuxi,j,k
uxp

z

solid
x

Fig.4

ui+1,j,kuci,j,kuxi+1,j,kdx
Ci,j,kGu i(,nj), k = −[adv _ term] x + [vis _ term] x

uci+1,j,k
delzk

)(

Vi , j , k

)

(26)

固体領域境界面の速度勾配のx方向成分uxpを、no slip
条件を仮定し以下のように表している。

ui +1, j , k(27)

上式において、dxは計算セルにおいて固体流域に覆われている境界面に相対する境界面から固体領域代表点までの距離を表しているが、圧力項の差分式で用いたものと同じで(24)式で表される。

数u,vを用いて連続の式により、以下のように計算される。

中心に配置された速度変数(uci-1, uci)を用いた境界面上における線形の内挿値に、2つの計算セルを一体とみなして計算される移流項と粘性項に、2つの計算セルの圧力から求められる局所的な境界面上の圧力勾配を足し合わせることにより、以下のように計算される。

ui(,nj+, k1) =(uci(,nj), k ⋅ delxi −1 + uci(−n1), j , k ⋅ delxi )delxi −1 + delxi⎧⎪ (Gu+ Gui(−n1), j , k ) grv ( Pi ,( nj ,)k − Pi (−n1), j , k ) ⎪⎫ (31)
+ Δt ⎨−⎬
0.5(delxi −1 + delxi ) ⎪⎭⎪⎩ Vi , j , k + Vi −1, j , k(n)i, j ,k

計算セルの中心に配置された速度変数では、圧力勾ているので、圧力値が1格子毎に激しく変化（計算モード）しても、それに対する抑制機構が働かない。
そこで、上式のように境界面に配置された速度変数の時間積分の際に、局所的な圧力勾配を加えてることに

wi , j ,k = C i , j , k −1 wi , j , k −1

よって、計算モードに対する抑制機構が強まるものと(28)

海面変位ztは、最上位層(k=nk)の計算セルにおいて、鉛直流速と同様にstaggered配置された速度変数u,vを用いて、以下のように計算される。

ップの値は、その境界面で接している2つの計算セルの圧力を隣接する2つの計算セルの圧力値の内挿を用い鉛直流速wはstaggered配置された境界面上の速度変+ Bi , j −1, k vi , j −1,k − Bi , j , k vi , j , k ] / C i , j , k

ることを示している。

配を求める際に、(19)式に示したように、境界面上の3.4.4 連続の式

+ [ Ai −1, j , k u i −1, j , k − Ai , j ,k u i , j , k

は時間ステップ数、そして下付き添え字xはx成分であ境界面上に配置された速度変数u(v)の次の時間ステ

ここで、KMは水平渦動粘性係数を表す。

dx

(30b)

Vi , j , k

表している。Δtは時間ステップ間隔、上付き添え字(n)∂ 2u= KM ×KM∂x 2
Ai , j , k u ix+1, j , k − uc ix, j , k + Ai +1, j , k − Ai , j , k u ix+1, j , k − u px

u px =(Gu i(,nj), k − [ prs _ term] x )（(14)式）、粘性項（(25)式）及び圧力項（(17)式）を

delxi

) (

uc i(,nj+, k1) = uc i(,nj), k + Δt

(30a)

adv_term,vis_term及びprs_termは、それぞれ移流項

Arrangement of variables for calculation of aviscosity term.

(

計算セルの中心に配置された速度変数uc(vc)の時間方向の積分は、有限体積法を適用した離散式（(13)式

Ci,j,k+1

uci-1,j,k

3.4.5 時間積分期待される。
3.5 Cut-Cell 法の問題点
Cut-Cell法において、計算セルにおける流体領域の割合が小さい場合、計算が不安定になる要因となる。
Adcroft2)らは、Cut-Cell法を用いる際に、計算の安定性を確保するためには、下式に示すように流体領域体平成17年

九 州 大 学 大 学 院 総 合 理 工 学 報 告

積を境界面の平均で割ることにより得られる実効的な

第27巻

第2号

(1) 計算領域

格子幅が、CFL条件を満たすように設定されることが・x軸方向: 80km

望ましいとしている。

・y軸方向: 60km

Δx eff ≈Δy eff ≈235

・z軸方向: 1000m （鉛直方向、深さ）Vi , j , k
0.5( Ai +1, j , k + Ai , j , k )(32a)(2) 格子Vi , j , k(32b)0.5( Bi , j +1, k + Bi , j , k )・格子数：(x:80)×(y:60)×(z:12)・格子間隔：

本モデルでは、海底地形を考慮して計算セルを構成

中心部：等間隔1km
格子数：(x:20)×(y:20)

する際（付録A.1参照）に、極端に小さい流体領域を持
つ計算セルを作らないために、使用する水深と計算セ

周辺部：不等間隔中心部から周辺部へ向けて1.1倍ルの上面（または下面）の距離が計算セルの高さの2の比で等比級数的に拡張

割未満の場合は、水深を上面（または下面）の位置に鉛直方向：12層変更するものとしている。
また、計算セルにおいて流体領域が占める割合が半

20m,30m,50m,100m×9層分以下の場合は、その計算セルの上部のものと合体し、一体のものとして通常の計算を行うこととした。計算(3) 海山形状により得られる変数値は合体した計算セル同士で共用

水深(m)は以下の関数によって表される。

するものとした。

depth = 1000
⎛ ( x − x0 ) 2 − ( y − y 0 ) 2 ⎞⎟− 900 exp⎜⎜ −⎟5000 2⎝⎠

4. 数 値 計 算 例作成したモデルの精度を検証するため、本モデルとMECモデルを用いてFig.5に示すような平坦な海底面にガウシアン型の海山を置いた簡易地形を使った数値(x0,y0): 海山山頂の座標位置 (Fig.5参照)
x0=y0=30km

計算を行い、結果を比較・検討した。
(4) 境界条件
西境界：流入North

流速値0.3m/s

水平・鉛直方向共に一様
東境界：流出

流速値0.3m/s

水平・鉛直方向共に一様
南北境界：Free Slip
East(5) 物理パラメータ
コリオリパラメータ： f=7.11×10-5 (s-1)
水平渦動粘性係数： 50 m2/sec （固定）
鉛直渦動粘性係数： 0.001 m2/sec （固定）

West
South

水温： 21℃（固定）

塩分： 34.8‰
Fig.5 Perspective plot of the seamount geometry, and thehorizontal numerical grid.

4.1 計算条件
計算条件を以下のとおり設定した。Cut-Cell法とMEC モ デ ル の 計 算 の 際 に は 同 じ 条 件 を 用 い た 。
Cut-Cell法とMECモデルでの海底地形の近似状態の違いを水平断面と鉛直断面についてFig.6に示す。

（固定）(6) 計算パラメータ
時間ステップ：

5 sec.

積分時間：

10 days

(33)

236

カットセル法を用いた海洋モデルについて(a) Vertical section at x=30km through the seamount top
Fig.6(b) Horizontal section at 550m deep

The representation of the seamount bottom in Cut-Cell model and MEC model.

l(a) Cut-Cell(c) MEC

depth=150m

depth=150m
Fig.7(b) Cut-Cell

depth=550m(d) MEC

depth=550m

Plain view of horizontal velocity vectors.

平成17年

九 州 大 学 大 学 院 総 合 理 工 学 報 告

第27巻

第2号

237

(a) Cut-Cell

depth=150m(b) Cut-Cell

depth=550m(c) MEC

depth=150m(d) MEC

depth=550m

Fig.8

Plain view of vertical velocity distributions.

4.2 計算結果

Fig.8に示した海山近傍の鉛直流でも、MECモデルで以下にCut-Cell法とMECモデルの計算結果を示す。

は、特に上流側（西側）に表れる上昇流の分布が分散

Fig.7には海山近傍における水平流速分布を、Fig.8

的であるのに対して、Cut-Cell法では、上流側でも上昇には鉛直流速分布を、それぞれ深度150mと550mの場合流の分布が滑らかになっている。これは、MECモデルについ示す。

では、海底地形を階段状に近似しているために、局所

Fig.7に示した海山近傍の水平方向の流れは、Cut-Cell

的に鉛直方向の強い流れが表れやすい傾向にあるため法の場合には、海山の地形に沿った流れが滑らかに表であると考えられる。Cut-Cell法では、海底地形の近似

現されているのに対して、MECモデルでは、海山の南の精度が高いために、鉛直方向の流れを滑らかに表現

北側面において、東向きの流れが強く現れ、海山の地することに成功している。2つのモデルで、下流側（東形に沿った流れの成分が弱まり、海山の後方の広い範側）の鉛直流分布が異なっているが、これは、主にFig.7

囲の流れが弱くなっている。これは、Fig.6(b)に示したに示したように2つのモデルで海山後方の流れの分布

ように、MECにおける海底地形の近似が粗く、格子のが違っているためであると考えられる。

並びに対して斜め方向の流れが表れにくい状況にあるためであると考えられる。一方、Cut-Cell法では、MECモデルと同じ格子数にも関わらず海底地形に沿った流れを滑らかに表すことができている。

238

5. 結

カットセル法を用いた海洋モデルについて

言

付

録

A1. 計算セルの体積及び表面積の求め方

急勾配を持った海底地形の周りでも、精度が高い流れの数値計算が可能となるようにMECモデルに有限計算セルの流体領域の体積及び表面積の計算方法を

体積法を用いたCut-Cell法を適用した。このモデルは、以下の鉛直方向に並ぶ2つの計算セルについて示す。上

格子生成の計算コストも低く、複雑な計算方式も必要

側を計算セル1、下側を計算セル2とする。

としない。ただし、実効的に計算セルを縮小するため、一般的に従来のz座標系モデルに比べ、計算する際の時

S4S3間ステップを小さくする必要がある。また、安定したB1

h4

計算を効率的に行うために計算セルの実効的な体積がh3S1小さくなり過ぎないように水深を微調整したり、小さS2A3くなった計算セルを他の計算セルと一体的に取り扱うA4といった工夫が必要である。

hm

z1

流体領域

Cut-Cell法とMECモデルを用いてガウシアン型の海h2

山を配した海底地形を用いた数値計算を行い、計算結果により、海山近傍の流れの精度が水平方向及び鉛直

M4 T
h1

D1

C1 B2

方向ともにCut-Cell法を適用することにより改善されることが確認された。また、海山近傍の流れの精度に

M3
C2

M1

M2

よる海山後方の流動場への影響が大きいことが確認さ
地形近傍の流動場の計算精度の向上に寄与するものと
参 考 文 献1) Mellor, G.L., Ezer, T., and Oey, L.-Y., “The Pressure GradientConundrum of Sigma Coordinate Ocean Models,” J. Atmos.
Oceanic Technol., 11, 1126-1134 (1994).
2) Adcroft, A., Hill, Chris, and Marchall, “Representation of
Topography by Shaved Cells in a Height Coordinate Ocean
Model,” Monthly Weather Review, 125, 2293-2315 (1997).
3) Landsberg, A. M., et al.: “The Virtual Cell Embedding
Method : A Simple Approach for Gridding Complex
Geometries,” AIAA-97-1982 (1997).
4) Ye, T., Mittal, R., Udaykumar, H.S., and Shvy, W., “Cartesian
Grid Method for Viscous Incompressible Flows with
Complex Immersed Boundaries,” AIAA-99-3312, 545-557(1999).
5) 矢部孝, 肖鋒, “固体・液体・気体の統一解放と CIP 法(2)”,
数値流体力学, 7, 103-114.
6) Saki, E.M., and Bringen, S., “Numerical Simulation of a
Cylinder in Uniform flow: Application of a Viertual Boundary
Method,” J. Comput. Phys., 123, 450-465 (1996).
7) 寺本進, 藤井孝蔵, “直交格子法による三次元物体周りの流れ解析”, 第 12 回数値流体力学シンポジウム講演論文集, 299-300.
8) 雷康斌, 岩田正子, 野田茂穂, 姫野龍太郎, “V-CAD デー
タを直接利用する二次元の任意形状流れ場の数値解析”,
ものつくり情報技術統合化研究（第２回）: Integrated
Volume-CAD System Research, 理化学研究所ものつくり
情報技術統合化研究プログラム
9) 日 本 造 船 学 会 , 海 洋 環 境 研 究 委 員 会 , “MEC Ocean
Model オペレーションマニュアル Version 1.1”, (2003)
10) 濱 田 孝 治 , “MEC Ocean Model に つ い て ”, 海 と 空 ,80, 27-36 (2004).

z2

A2

れ、本モデルで用いたCut-Cell法が海山等の複雑な海底期待される。

D2

固体領域

A1

(a) Side view of the two cells

M4

M3

x
s
B2

y

D1
dx1

C1

dy3 ds1

dy1

M1

D2
C2
dy2

dx2
ds2

dy4

M2

(b) Horizontal interface between the two cells
Fig.A1

Definition of the points in the program.

A1.1 定数・変数
A1.1.1 計算セル基本定数x,y: 水平方向格子幅

z1-2 : 鉛直方向格子幅
h1-4 : 水深（計算セル1の上面を基準とする）

平成17年

九 州 大 学 大 学 院 総 合 理 工 学 報 告

第27巻

第2号

239

+(z1-hm){x·y/4-[(x-dx1-dx2)(dy1+dy2)
A1.1.2 主要変数2

+dx1·dy1+dx2·dy2]/2}/2-x·y·z1/4

2 1/2

s=0.5(x +y )hm=(h1+h2+h3+h4)/4

A1.2.3 計算セル 1

ds1=s[1-(z1-hm)/(h1-hm)]

計算セル1における流体領域の体積

ds2=s[1-(z1-hm)/(h2-hm)]=[流体領域全体の体積]-[計算セル2の流体領域の体積]

dx1=0.5x·ds1/s
dx2=0.5x·ds2/s

A1.3 表面積dy1=0.5y·ds1/s
dy2=0.5y·ds2/s

A1.3.1 計算セル 1 側面

dy3= y(h1-z1)/(h1-h4)

境界面[S1S2-M1M2]=x·z1

dy4= y(h2-z1)/(h2-h3)境界面[S2S3-M2D2]=y·z1-(y-dy4)(z1-h3)/2
境界面[S3S4-A3A4]=(h3+h4)x/2

A1.2 流体領域の体積
A1.2.1 全体（計算セル 1+計算セル 2）

境界面[S4S1-D1M1]=y·z1-(y-dy3)(z1-h4)/2
A1.3.2 計算セル 2 側面

流体領域の体積（計算セル1+計算セル2）

境界面[M1M2-A1A2]=(h1+h2-2dz1)x/2=立体[S1S2B1-A1A2T]+立体[S2S3B1-A2A3T]

境界面[M2D2-A2]=dy4(h2-z1)/2+立体[S3S4B1-A3A4T]+立体[S4S1B1-A4A1T]

境界面[M1D1-A1]=dy3(h1-z1)/2立体[S1S2B1-A1A2T]の体積= (h1+h2+hm)x·y/12

A1.3.3 鉛直方向境界面
境界面[S1S2S3S4]=x·y
境界面[D1C1C2D2-M2M1]=[(x-dx1-dx2)(dy1+dy2)立体[S2S3B1-A2A3T]の体積+dx1·dy1+dx2·dy2]/2= (h2+h3+hm)x·y/12立体[S3S4B1-A2A3T]の体積= (h3+h4+hm)x·y/12立体[S4S1B1-A3A4T]の体積= (h4+h1+hm)x·y/12
A1.2.2 計算セル 2

A2. QUICK 補間係数の決定
QUICK (quadratic upstream interpolation for convective
kinematics)は、セル境界の物理量φi+1/2の内挿法として、
その境界の前後を挟んでいる2つのセルの値の他に、風上側にもう1つ離れたセルの値を用いて2次の多項式による内挿を行うものである。
下図において、格子間隔が等しい場合

計算セル2における流体領域の体積(dxi-1=dxi=dxi+1)、補間式は以下のようになる（流れは、=立体[A1-C1D1M1] + 立体[A2-C2D2M2]

iが増加する方向へ向いているものとする）。物理量φi-1,
φi, φi+1は、セル中心で定義されているものとする。

+ 立体[A1A2-C1C2M2M1]
立体[A1-C1D1M1]の体積=dx1·dy3(h1-z1)/6

φ 1 =
i+

2

−φi −1 + 6φi + 3φi +1
8(A1)

立体[A2-C2D2M2]の体積=dx2·dy4(h2-z1)/6

φi-1

φi

dxi-1

dxi

(φi+1/2)

φi+1

立体[A1A2-C1C2M2M1]=立体[S1S2B1-A1A2T]+立体[T-B2C1C2]-立体[S1S2B1-M1M2B2]=(h1+h2+hm)x·y/12

Fig.A2

Positions of variables.

dxi+1

240

カットセル法を用いた海洋モデルについて作成したモデルでは、variable meshを用いているので、格子間隔が等しくない場合のQUICK補間係数が必要となる。そこで下式において、精度が3次以上となるように係数(a,b,c)を定めることを考える。

φ 1 = aφi −1 + bφi + cφ i +1
i+

(A2)

Δx −1 = 0.5dx i −1 + dx i

(A4)

Δx 0 = 0.5dx i
Δx1 = 0.5dx i +1

(A2.3)式を(A2.2)式に代入し補間係数に関してまと
めると、1次と2次の誤差項が”0”になるという条件のも2

φi-1、φi及びφi+1をφi+1/2の回りでテーラー展開すること
により次式が得られる。

∂φ(Δx −1 ) 2 ∂ 2φ i +1 / 2
φ i −1 = φ 1 − Δx −1 ⋅ i +1 / 2 +
+K
⋅
i+
2
∂x
∂x 2
2

φ i = φ 1 − Δx 0 ⋅i+

2∂φ i +1 / 2 (Δx 0 ) 2 ∂ 2 φ i +1 / 2
+K
+
⋅
2
∂x
∂x 2

∂φ
(Δx1 ) 2 ∂ 2φ i +1 / 2
φ i +1 = φ 1 + Δx1 ⋅ i +1 / 2 +
+K
⋅
i+
2
∂x
∂x 2
2
(A3)
ここで、Δx-1、Δx0及びΔx1は、それぞれφi-1、φi及びφi+1の定義点とφi+1/2の定義点の距離を表しており、次式で与えられる。

と以下の行列方程式が得られる。

⎡ 1⎢ − Δx−1⎢⎢⎣(Δx −1 ) 2

1− Δx 0(Δx 0 ) 2

1 ⎤ ⎡ a ⎤ ⎡1 ⎤Δx1 ⎥⎥ ⎢⎢b ⎥⎥ = ⎢⎢0⎥⎥(Δx1 ) 2 ⎥⎦ ⎢⎣ c ⎥⎦ ⎢⎣0⎥⎦(A5)

これを、以下のように解くことによって等しくない
格子間隔の組み合わせに対する3次精度の補間係数を求めることができる。

⎡a ⎤ ⎡ 1⎢b ⎥ = ⎢ − Δx−1⎢ ⎥ ⎢⎢⎣ c ⎥⎦ ⎢⎣(Δx −1 ) 2

1− Δx 0(Δx 0 ) 2

1 ⎤Δx1 ⎥⎥(Δx1 ) 2 ⎥⎦

−1

⎡1 ⎤
⎢0 ⎥
⎢ ⎥
⎢⎣0⎥⎦

(A6)