Non-uniform subcell scheme

次式で評価する．

(q_r)^m_i+1/2,j,k∆θ_j∆z_k∆t

=sign

(qr)^m_i+1/2,j,k 

∆θ_j∆z_k 2







max



ri+1/2− (q_r)^m_i+1/2,j,k∆t ri+1/2 ,ri+1/2





2

−min



ri+1/2 − (q_r)^m_i+1/2,j,k∆t ri+1/2 ,ri+1/2





2



(4.2.34)

また，(q_θ)^m_i,j+1/2,k∆r_i∆z_k∆tは，

(qθ)^m_i,j+1/2,k∆ri∆zk∆t= 1 2

r_i+1/2² −r_i−1/2²

(qθ)^m_i,j+1/2,k∆t∆zk

r_i (4.2.35)

により評価する．(q_z)^m_i,j,k+1/2∆t(r²_i+1/2−r²_i−1/2)∆θ_j/2は，この離散式のまま評価する．

r_i+1/2

|(q_r/r)_i+1/2,j,k∆t|

|(q_r/r)_i−1/2,j,k∆t|

(q_r)_i+1/2,j,k (q_r)_i−1/2,j,k

∆r^sub ∆θ^sub

∆z^sub r^sub

r_i−1/2

r_i+1/2 −(q_r/r)_i+1/2,j,k∆t

axis r = 0

i+1/2 face i−1/2 face

transferred volumes, N_rp, N_rm non-transferred

volume, N_C

subcell interface

図4.2.2 界面セルに準備したサブセル

∆Ωi−1/2,j,k =−1 2

"

2ri−1/2− qr

r

i−1/2,j,k∆t

# qr

r

i−1/2,j,k∆t∆θj∆zk

である．この被輸送体積の定義は，ポアッソン方程式導出時における連続の式の扱い方と 同様である．これら被輸送体積中に含まれる体積率αを，セル内の界面位置を考慮して算 出する．界面はセル(i, j,k)内で平面と仮定する．平面の方程式は次式で与えられる．

n·xint +φint = 0 (4.2.39)

ここで，φint は界面‐セル中心間の符号付距離，xint はセル中心を始点とする界面上の点 への距離ベクトルである．nは界面法線であり，次式で計算する．

n= ∇φ

|∇φ| (4.2.40)

φは局所レベルセット関数である．φについては4.2.5節で述べる．

サ ブ セ ル を 利 用 し て φ_int を 求 め る ．図 4.2.2 に 示 す よ う に ，被 輸 送 体 積

∆Ωi+1/2,j,k, ∆Ωi−1/2,j,k お よ び 残 留 体 積 に 各 々 Nr p, Nrm, NC 個 の サ ブ セ ル を 準 備 す

る ．図 中 ∆r^sub, ∆θ^sub, ∆z^sub は r, θ, z 方 向 の サ ブ セ ル 幅 で あ る ．サ ブ セ ル 体 積 Ω^sub = 0.5(2r^sub−∆r^sub)∆r^sub∆θ^sub∆z^sub をすべてのサブセルについて加えると，セル体積 Ωi,j,k =0.5(2ri−1/2−∆ri)∆ri∆θj∆zkに一致する．

N_C

X

l=1

Ω_l^sub+

Nr p

X

m=1

Ω_m^sub+

N_rm

X

n=1

Ω^sub_n = Ωi,j,k (4.2.41)

cell center

subcell

d₍₁₎ d₍₂₎

d₍₄₎ d₍₃₎

d₍₈₎

d₍₇₎ d₍₅₎ d₍₆₎ interface

n

condition for d_(v) d_(v)^. n < d_(v+1)^. n

図4.2.3 距離ベクトルd(v)

ただし，r^sub は，rが小さい方のサブセル表面のr座標である(図4.2.2)．サブセルに含ま れる流体体積のセル体積Ωi,j,k に対する割合を，サブセル体積率α^sub と定義する．定義よ り，α^sub の値がとりうる範囲は0 < α^sub < Ω^sub/Ω_i,j,kである．α^sub はφ_int の関数である．

α^sub(φint)を次の区分一次関数で与える．

α^sub_m (φ_int)=max









min



 Ω^sub_m

4φintn−P₈

v=5d(v)

·n

Ωi,j,k

P₄

v=1 d(v)−P₈

v=5d(v)

·n , Ω^sub_m

Ωi,j,k



,0









 (4.2.42)

d_(v) はセル中心からサブセルmの頂点vへの距離ベクトルである(図4.2.3)．次式を満足 するように頂点番号を付す．

d(v)·n≤ d(v+1) ·n (4.2.43)

すべてのサブセル体積率 α^sub の総和が αi,j,k に一致するという制約条件から次式が成り 立つ．

f (φint)=

NC

X

l=1

α_l^sub(φint)+

N_{r p}

X

m=1

α_m^sub(φint)+

Nrm

X

n=1

α_n^sub(φint)−αⁿ_i,j,k= 0 (4.2.44) 上式を反復法により解いて，φint を得る．本研究ではブレント法[18, 19]を用いて上式を 解く．

式(4.2.44)の解φint により，式(4.2.38)右辺大括弧内の２項は各々以下の諸式で与えら

れる．

(αqr)_i+1/2,j,k∆t∆θ∆z=

Nr p

X

m=1

α^sub_m (φint)Ωi,j,k (4.2.45)

(αqr)_i−1/2,j,k∆t∆θ∆z=−

N_rm

X

n=1

α_n^sub(φint)Ωi,j,k (4.2.46)

これらの値を式(4.2.38)に代入して，αⁿ_i,j,kの値を更新する．

4.2.5 レベルセット関数

2.5節で提案した局所レベルセット関数の概念を３次元円柱座標系用に拡張する．界面 セルの距離関数は式(4.2.44)の解として得られる．

φi,j,k =φint (if 0< αi,j,k <1) (4.2.47)

界面セルの近傍セルのφ は，φint に基づいて与える．φはスカラー量であるため，第 2 章で説明したデカルト座標における方法をそのまま拡張できる．界面の存在しないセル (i, j,k)のφを次式で与える．

φ_i,j,k = P

l,m,n∈NWl,m,nφl,m,n

P

l,m,n∈NWl,m,n

(4.2.48)

ここで，W_l,m,n は重み関数であり，

Wl,m,n =( r_i,j,k− rl,m,n^−p if 0< αl,m,n <1

0 otherwise (4.2.49)

と定義した．式(4.2.48)中の N は，セル(i, j,k)を中心とした，重み付き平均をとる範囲 に含まれるセル数を意味する．また，式(4.2.49)中の pは経験的に与えるパラメーターで ある．

一方，次の移流方程式を解いてレベルセット関数の時間変化を求めることもできる．

∂φ

∂t +u· ∇φ=0 (4.2.50)

この場合，重み付き平均は用いない．上式にVerziccoの変数変換を適用し，連続の式を用 いて保存形に変形すると，次式を得る．

r∂φ

∂t + ∂qrφ

∂r + ∂qθφ

∂θ +r∂qzφ

∂z =0 (4.2.51)

移流方程式を解くと，距離関数としての性質(|∇φ|=1)が失われる．距離関数としての性 質を回復させるために，移流方程式に続けて，次の再初期化方程式を解く．

∂φ

∂τ = sign_ε(φ₀) (1− |∇φ|) (4.2.52)

ここで，τは擬似的な時間，φ0 は再初期化開始前のφの値，sign_ε は次式で定義される平 滑なsign関数である．

sign_ε(φ0)= φ0

q

φ²₀+h²

(4.2.53)

hはセル幅を表わす．ただし，移流方程式および再初期化方程式を解く際，φとαの記述 する界面位置が整合するように，界面セルのφの値をφint に固定する．

4.3 ３次元一般曲線座標系用 NSS

ドキュメント内 2007 2 (ページ 140-144)