移流拡散方程式への Tangent 変換 CIP 法の適用性

(1)

移流拡散方程式への Tangent 変換 CIP 法の適用性

日大生産工（学部） ○ 家塚史仁日大生産工登坂宣好

１はじめに

連続体モデルとしての流れの数値シミュレーションの課題は，移流方程式の安定かつ高精度近似スキームの構築である．そのスキームの一つと

してCIP法(Cubic Interpolated Propagation

Method)が存在している¹⁾．CIP法の特徴は，物

理量のみならず微分量をも独立変数として用いることにより，補間精度を向上させていることにある．固体・液体・気体の三態の統一解法として，

また圧縮性，非圧縮性流体を統一的に扱える手法として，これまで流体力学をはじめ様々な分野で応用され成果を上げている．しかし，このCIP法の適用性に関する考察は少ない．

本論では，まず CIP 法の基本的な考えを示した後，2次元の非定常移流拡散方程式に対して同手法を適用し，2つの例について数値計算を行う．

また，長い時間ステップでの計算結果について，

スキームの適用性を検証する．

２ CIP 法基礎概念 2.1 CIP 法

移流方程式

) 0 (

0 = >

∂ = + ∂

∂

∂ u const

x u f t

f (1)

の解はよく知られているように，次式で表される．

) 0 , ( ) ,

(x t f x ut

f = − (2) これは初期のプロファイルが速度で移動することを表している．初期に図1の実線のような物理量のプロファイルがあるとする．

で時間移動するとき，厳密解は破線のようになるはずである．しかし，データは離散的に，格子点上でしか得られないので，セル内のプロファイルを単純に補間してしまうと点線のようになり，厳密解からはズレたものになってしまう．そこで式(1)を微分すると

) 0

(t< u

>0 u

∆t

) / (

0 g f x

t u g t

g = =∂ ∂

∂ + ∂

∂

∂ (3)

となり，gもまたで伝播することを表している．

従って，物理量同様，微分量も移流させれば，格子点上では移動後のプロファイルに図の太矢印

u

図1 CIP法の考え方

のような制限が加わり，メッシュ間のプロファイルが移動前に近くなることになる．

CIP法では格子点i−1,iの間のプロファイルを次の3次関数で補間する．

i i i i i i i

i x a x x b x x c x x d

F( )= ( − )³+ ( − )²+ ( − )+ (4)

4つの未知係数は関数が格子点上で与えられた値及び微分値を取るという次の条件

i i i

i b c d

a, , , F_i(x) , _i₋1

i f

f g_i,g_i₋₁

, 2

) 3 (

, )

( ) , (

, )

(

1 1 2

1 2

3 1

− −

−

= +

∆

−

∆

=

= +

∆

−

∆ +

∆

−

=

i i i i

i i

i i i i i i

i

i i i i

g c x b x dx a

x dF

f d x c x b x a x

F

g dx c

x dF

f d x F

(5)

から未知数が次のように決定できる．

2 . ) (

3

), (

2

1 2

1

3 1 2

1

x g g x

f b f

x f f x

g a g

i i i i i

i i i

i i

∆ + +

∆

= −

∆

− −

∆

= +

−

(6)

ただしは式(5)で既に与えてある．この結果，

次の時刻での値は，このプロファイルを

i

i d

c,

u∆tだけ遡ったものであるから

. 2

3

,

2 1

2 3 1

n i i i

n i

n i n i i i n i

g b a

g

f g b a f

+ +

=

+ + +

=

+ +

ξ ξ

ξ (7)

ここで，ξ =−u∆t．なお，u<0の場合にはi−1を +1

i に，−∆xを∆xに置き換えればよい．

Applicability to Advection-Diffusion Equation of CIP Method with Tangent Conversion

Fumihito IETSUKA and Nobuyoshi TOSAKA

(2)

2.2 Tangent 変換 CIP 法

CIP法でも十分精度のよい計算が行えるが，物体や液体の表面等の記述にはまだ不十分である．

そこで関数変換の方法を用いる．

Tangent 関数変換では式(1)でのを解く代わ

りに

f

(

12

tan )

(f = f −

H π

))

(8) の変換を行う．Hについても式(1)と同じ式が得られることはすぐに分かる．そこで，f 代わりに

の Hを

2 次元非定常移流拡散方程式

適用する

未知数として式(8)を使い，必要に応じてこの結果を逆変換し f を求める．この変換は物質の境界面を記述する密度関数のように，0または 1の指標で表される場合には非常に有効である．

３

次の2次元移流拡散方程式へCIP法を

．

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

∂ +∂

∂

= ∂

∂ +∂

∂

∂f

2 2 2

) 2

( ) (

y f x

f y

vf x

uf

t κ (9)

ここで，κ は拡散係数である．左辺を少し変形し

y G f x

f y

f v x f u

y v f u f

f + + ∂ =

∂

x t

⎟⎟≡

⎠

⎜⎜ ⎞

⎝

⎛

∂ +∂

∂ + ∂

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

∂ + ∂

∂

− ∂

∂

2 2 2

κ 2

(10)

としておく．CIP法では空間微分も必要になる．

x, f v f u

f G f v

f_x u _x _x ∂

∂ −

−

∂ =

∂ +

(11) x

y x

t ∂ ∂ ^x ^x ∂ ^y ∂

∂

y. f v y f u y G

v f x u f t f

y x

y y y y

∂

− ∂

∂

− ∂

∂ = + ∂

∂ + ∂

∂

∂ (12)

ここで，x, yの下付き添え字は各量での偏微分を流項と非移流項を分離して解く．

す x y x^,fy

)

と中間の値を求める．次に非移流相で表す．

CIP法では移

なわち，まず移流相を計算し

(

^fⁿ^,^fⁿ^,^fⁿ

)

^→

(

^f^*^,^f^* ^*

(

^*^, ^*^, ^*

)

^→

(

⁺¹^, ⁺¹^, yⁿ⁺¹

)

n x

n f f

f f

f

f _x _y

を求め，次のステップの値とする．

移流相としては

0 ,

0 , 0

∂ = + ∂

∂ + ∂

∂

= ∂

∂ + ∂

∂

∂ = + ∂ + ∂

∂ f

f v f u

∂

y v f x u f t f y v f x u f t f

y x t

y y x y

x x

(13) であり，CIP補間を用いて計算する．非移流相では以下の式を差分近似より求める．

,

, x

f v x f u t G

G f t

f _x

∂ −

−

∂ =

y. f v y f u t G

f

y x y y

∂

− ∂

∂

− ∂

∂ =

∂

(14) まず，移流相であるが，2次元空間での空間プロファイルは1次元同様に3次関数で補間する．

この補間関数の選び方には幾つかあるが，ここでは以下のような関数を採用した(A型CIP)．

2 2 , 1 1

, 1 2 1 , 2

1 , 0 2 , 0 3 , 0

0 , 0 0 , 1 2 0 , 2 3 0 , 3

, ( , )

XY C XY C Y X C

Y C Y C Y C

C X C X C X C y x F_i_j

+ +

+

+ + +

+ +

+

=

(15)

ここで，X =x−x_i,Y =y−y_i．これは，1次元CIP をそれぞれの方向に拡張しただけであり，項数が少なくてすむので，多次元への拡張には適している．ここにはCα,βの10個の未知係数が現れるが，

格子点(i,j),(i+1,j),(i,j+1)でそれぞれ与えられたの値を持ち，点での値を持つという10個の関係を用いると，次のように係数が決定できる．

y

x f

f

f, , (i+1,j+1) f_i₊₁_,_j₊₁

). (

) ( )

(

), ( )

( ) ( ) (

), (

) ( )

(

), ( ) 2 (

) (

3

), (

) (

2 ) (

, ,

,

), ( ) 2 (

) (

3

), (

) (

2 ) (

2 2 , 1 1 , 1 1 2 , 2

, ,

2 , , 2 , 0

3 , , 2

, , 3 , 0

, 1 , 0 ,

0 , 0 ,

0 , 1

, , 2

, , 0

, 2

3 , , 2

, , 0

, 3

u sign y x

v sign N u sign y x C L

u sign x

N v

sign y

M v sign u sign y x C L

v sign y x

u sign M v sign y x C L

y

v sign f

f y

f C f

v sign y

f f y

f C f

C f

C

x

u sign f

f x

f C f

u sign x

f f x

f C f

j yi jup yi j yi jup i

jup i j i j

yi jup yi

j yi j

i j

xi

j xi j xiup j

i j iup

j iup j i j

xi j xiup

×

∆

− ×

×

∆

= −

×

−∆

×

−∆

×

∆

= −

×

∆

− ×

×

∆

= −

∆ + +

∆

= +

×

∆

− −

∆

= +

=

∆ + +

∆

= −

×

∆

− −

∆

= +

(16) ただし，

. ,

,

, ,

, , ,

,

j yi j yiup

j xi jup xi

jup iup j iup jup i j i

f f N

f f M

f f f f L

−

=

−

=

+

−

=

(17)

ここで

⎩⎨

⎧

<

+

≥

= −

⎩⎨

⎧

<

+

≥

= −

) 0 ( 1

) 0 ( , 1

) 0 ( 1

v j

v jup j

u i

u

iup i (18)

y x

x ∂

∂

(3)

⎩⎨

⎧

<

−

= ≥

) 0 ( 1

) 0 ( ) 1

( x

x x

sign (19) という記号を導入した．これはの正負によりプログラムを場合分け(2 次元ならば上流が4 通り存在する)しなくとも良いようにするためである．

v u,

非移流相は以下の差分で計算する．式(21), (22) では，式(20)での微分がすでに用いられているので，をで置き換えて計算している．これにより，式(21), (22)ではの微分は必要なくなる．

G

y

x G

G , (fⁿ⁺¹− f^*)/∆t

G

2 , 2

2 2

1 , 1

* , , , 1 ,

* 1 ,

2

* 1 ,

* ,

* 1 , 2

* , 1

* ,

* , 1

* , 1 ,

y t u f u

x u f u

y f f f

x f f f

f f

j i j i j i j i j i j i

j i n

j i

⎭∆

⎬⎫

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

∆

− −

∆

− −

⎪⎩

⎪⎨

⎧

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

∆ + + −

=

− +

+

κ

(20)

2 , 2

2

, 1 ,

* 1 , ,

1 ,

* 1 ,

* , 1

* , 1 1 , 1 1 ,

* 1 , 1 ,

x t v f v

x t u f u

x f f f f f

f

j i j i j yi j

i j i j xi

j i j i n

j i n

j i j xi n

j xi

∆ ∆

− −

∆ ∆

− −

∆ +

− + −

=

− +

− + +

− + + +

(21)

2 . 2

2

1 , 1

* , , 1

, 1

* , ,

* 1 ,

* 1 , 1

1 , 1

1

* , , 1 ,

y t v f v

y t u f u

y f f f f f

f

j i j i j yi j

i j i j xi

j i j i n

j i n

j i j yi n

j yi

∆ ∆

− −

∆ ∆

− −

∆ +

− + −

=

− +

− + +

− + + +

(22) 2次元の場合にもTangent変換CIP法では，

上記の f の代わりに式(8)で変換されるH を用いることは，全く同様である．

４数値計算例 4.1 Skew Flow 問題

図2に問題設定を示す．メッシュに対して斜め一定方向の流れであり，常に u =1とする．

とし，20×20の同サイズの固定正方形格

子を使用する．流入境界は不連続(0 または 1)であり，流出境界は自然境界条件( する．

2通りの

10⁻6

κ =

)

=0

fn と

θ(= 45^°, 67.5^°)について，計算結果の立

面図を図3〜6に示す．

図2 Skew Flow問題

図3 θ= 45°(CIP法)

図4 θ= 45°(Tan変換CIP法)

図5 θ= 67.5°(CIP法)

図6 θ= 67.5°(Tan変換CIP法)

(4)

4.2 Rotating Cone 問題

図7に問題設定を示す．初期形状が文字“C”

になるよう 0 または 1 を与え，一定方向に回転するベクトル場

( ) (

u,v = −y,x

)

を与える．

とし，100×100の固定正方形格子を用いる．境界上では常に値を0とする．1回転後の立面図及び圧力図を図 8, 9に，100回転後の立面図を図 10に示す．

10⁻9

κ =

図7 (a) Rotating Cone問題 (b) 初期状態

図8 1回転後(CIP法) (a) 立面図，(b) 圧力図

(a)

図9 1回転後(Tan変換CIP法) (a) 立面図，(b) 圧力図

図10 100回転後(Tan変換CIP法)

５おわりに

本論では，CIP法を2次元非定常移流拡散方程式に適用し，2つの数値計算例について同手法の適用性を検討した．

CIP 法単独では二例とも若干オーバーシュートが起こり，また数値拡散もみられる．一方，

Tangent 変換CIP法では，何れの問題でも境界

をシャープに捕らえオーバーシュートもない．式 (8)のように簡単な変換でこれだけの効果が得られることは，非常に意味が大きい．しかし，Rota-

ting Cone問題での長い時間ステップでの結果は，

初期形状とは全く違ったものとなった．これは

Skew Flow 問題では流れの方向が常に一定であ

るのに対し，Rotating Cone問題では流れの方向が常に回転していることが原因と思われる．

今後は，より一般的な方程式への CIP 方の適用性を検討していきたい．

参考文献

1) 矢部孝，内海隆行，尾形陽一，CIP法，森北出版，(2003)．

2) Brooks, A.N., and Hughes, T.J.R., Streamline upwind / Petrov-Galelkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Meths. Appl. Mech.

Engrg., 32, (1982), pp.199-259.

(a)

(b) (a)

(b) (b)

移流拡散方程式への Tangent 変換 CIP 法の適用性