数理解析研究所講究録
第 1800 巻 2012 年 216-225

216

窓関数を用いたスペクトル法による非周期的な場の時間発展の数値シミュレーション法の検討

大阪府立大学工学研究科

榑林 聖子 (Seiko Kurebayashi), 村上 洋一 (Youichi Murakami)

Graduate School of Engineering,Osaka Prefecture University
1

はじめに偏微分方程式の数値解法のひとっであるスペクトル法は，周期境界条件のもとでは非常に精度よく空間微分を計算できる．しかし境界が不連続な場合，フーリエ級数で関数を表すと不連続点の近傍で激しく振動し級数の収束性が悪くなり (ギブスの現象), 実用的でなくなる．本研究は，境界が不連続な初期条件においてもギブスの現象を回避してスペクトル法を用いて精度よく計算する手法の提案を目的としている．

現在この手法には，周辺領域法，ウィンドウ法が提案されている．Shlatter et al[1] は 2 次元
非圧縮の $Navier^{-}Stokes$ 方程式に周辺領域法，ウィンドウ法を適用して両者の比較を行っている．

本研究では，ウィンドウ法を用いて Burgers 方程式を計算した．まず，窓関数を掛け合わせた新しい関数を導入し，Burgers 方程式を修正した方程式を用いてシミュレーションを行った．
Burgers 方程式は，非線形項と散逸項を持つにもかかわらず，非定常な厳密解が存在する．したがって，本研究ではウィンドウ法の有効性を定量的に調べることができる．

2 数値計算法
2. lBurgers 方程式の修正$u(x, t)$ について
領域 [0,L] において非周期的 (左右の境界で値が不連続) な関数$u(x, t)$ を考える．の Burgers 方程式は，(1)$\frac{\partial u}{\partial t}=-\frac{1}{2}\frac{\partial u^{2}}{\partial x}+v\frac{\partial^{2}u}{\partial x^{2}}$で与えられる．式 (1) をスペクトル法を用いて精度良く計算するために，窓関数灰 x) によって $u(x, t)$を人工的に周期的な関数$\tilde{u}(x, t)$にする．Burgers 方程式を修正し，新しい関数$\tilde{u}(x, t)$についてのBurgers 方程式を導出する．
まず，$\tilde{u}(x, t)=W(x)u(x, t)$ ,(2)$\overline{u^{2}}(x, t)=W(x)\{u(x, t)\}^{2}$(3)を定義する．ここで，W(x) は両端近傍で急速に に減少する関数である．
$0$$\frac{\partial\tilde{u}}{\partial t}=\frac{\partial}{\partial t}(Wu)=W\frac{\partial u}{\partial t}$

(4)

より，
$\frac{\partial u}{\partial t}=\frac{1}{W}\frac{\partial\tilde{u}}{\partial t}$.

(5)

217次に，(6)$\frac{\partial}{\partial x}(\frac{1}{2}\overline{u^{2}})=\frac{\partial}{\partial x}(\frac{1}{2}Wu^{2})=W\frac{\partial}{\partial x}(\frac{1}{2}u^{2})+\frac{1}{2}u^{2}\frac{dW}{dx}$より，$\frac{\partial}{\partial x}(\frac{1}{2}u^{2})=\frac{1}{W}\{\frac{\partial}{\partial x}(\frac{1}{2}\overline{u^{2}})-\frac{1}{2}u^{2}\frac{dW}{dx}\}$.

(7)

最後に，
(8)

$\frac{\partial^{2}\tilde{u}}{\partial x^{2}}=\frac{\partial^{2}}{\partial x^{2}}(Wu)=\frac{1}{w}\frac{dW}{dx}(\frac{\partial\tilde{u}}{\partial x}-u\frac{dW}{dx})+W\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial}{\partial x}(u\frac{dW}{dx})$より，$\frac{\partial^{2}u}{\partial x^{2}}=\frac{1}{w}\{\frac{\partial^{2}\tilde{u}}{\partial x^{2}}-\frac{1}{w}\frac{dW}{dx}(\frac{\partial\tilde{u}}{\partial x}-u\frac{dW}{dx})-\frac{\partial}{\partial x}(u\frac{dW}{dx})\}$.

(9)以上得られた式 (5), (7), (9) を式 (33) に代入し，方程式全体に窓関数 $W(x)$ をかけると，

.

(10)

$\frac{\partial\tilde{u}}{\partial t}=-\frac{\partial}{\partial x}(\frac{1}{2}\overline{u^{2}})+\frac{1}{2}u^{2}\frac{dW}{dx}+v\{\frac{\partial^{2}\tilde{u}}{\partial x^{2}}-\frac{1}{W}\frac{dW}{dx}(\frac{\partial\tilde{u}}{\partial x}-u\frac{dW}{dx})-\frac{\partial}{\partial x}(u\frac{dW}{dx})\}$

式 (10) のすべての項は周期化されている．以上より，$\tilde{u}(x)$について修正した Burgers 方程式が得られた．

22

数値スキームについてまず，(1)(2) 式より，実空間で，$\tilde{u}(x, t),\overline{u^{2}}(x, t)$を作る (windowing).

式 (10) の各項を，それぞれフーリエ級数に展開すると，$\tilde{u}(x, t)=\sum_{k=-\infty}^{\infty}a_{k}(t)exp[-ikx]$

,

(11)
,

$\frac{1}{2}\{\tilde{u}(x, t)\}^{2}=\sum_{k=-\infty}^{\infty}b_{k}(t)exp[-ikx]$$u \frac{dW}{dx}=\sum_{k=-\infty}^{\infty}c_{k}(t)exp[-ikx]$

(12)

,

(13)

$\frac{1}{2}u^{2}\frac{dW}{dx}+\frac{v}{W}(\frac{dW}{dx})^{2}u=\sum_{k=-\infty}^{\infty}d(t)exp[-ikx]$

(14)

,

.

(15)

$-v \frac{1}{W}\frac{dW\partial\tilde{u}}{dx\partial x}=\sum_{k=-\infty}^{\infty}e_{k}(\tau)exp[-ikx]$

以上を式 (10) に代入すると，各フーリエ成分について，(16)$\frac{da_{k}}{dt}=ikb_{k}(t)-vk^{2}a_{k}(t)+ivkc_{k}(t)+d_{k}(t)+e_{k}(t)$が得られた．式 (16) を時間について前進オイラー法を用いて差分化すると，$a_{k}(t+\Delta t)=a_{k}(t)+\Delta t\{ikb_{k}(t)-vk^{2}a_{k}(t)+ivkc_{k}(t)+d_{k}(t)+e_{k}(t)\}$

これを実空間に戻して，$\tilde{u}(t+\Delta t)$

.

(17)

が得られる．次に，$u(x, t+\Delta t)=(1-W)u(0)+\tilde{u}(x, t+\Delta t)$(18)$u(x, t+\Delta t)$ を得る (dewindowing).
から，この操作をすることで，もとの境界条件を取り入れるこ218とになる．以上を繰り返す．

23

窓関数と領域について

今回は、 以下に示す窓関数 [1] を採用した．
(19).
図 1, 2 に $W(x)$ と $dW(x)/dx$ の概形を示す．この関数はなめらかで，何回でも微分可能であり，微$W(x)=10^{-a^{n}|2(x-x_{L})/(x_{R}-x_{L})-1|^{\mathfrak{n}}}$分の回数が増えると，極値の数と極値の絶対値がともに増大する．窓関数の中心の平らな領域が正しくシミュレーションしたい領域 (物理領域 l) である．2 つのパラメータ $a,$ はそれぞれ物理領域$\Gamma$$n$の幅と，窓関数の 1 から への落ち方を決める．は大きいほど平らな領域はせまくなり，は大き
$0$

$n$

$a$

いほど落ち方が急になる．

3 数値計算の結果
以下の左右の境界で値が不連続な初期条件，$u(x)=-2v t\frac{k}{2}+\frac{k}{2}\tanh(\frac{kx+x_{0}}{2})\}$

(20)

,

(21)

$u(x)=-2v \frac{k_{1}\exp(k_{1}x+x_{1})+k_{2}\exp(k_{2}x+x_{2})}{1+\exp(k_{1}x+x_{1})+\exp(k_{2}x+x_{2})}$のもとでシミュレーションを行$A\searrow$得られた数値解と厳密解を比較した．(20) 式は 1 つの衝撃波が振幅の幅によってきまる一定の速度で進行する場合，(21) 式は 2 つの衝撃波がそれぞれの振幅によってきまる一定速度で進行していき，やがて融合し，最終的に 1 つの衝撃波として進行していく$v=0.01,$ $k=-10,$ $x=50,$場合の初期条件である．初期条件のパラメータは，$60,$$k_{2}=-20,$$x_{2}=100$ ,$k_{1}=-10,$$x_{1}=$窓関数のパラメータは $a=1.18,$ $n=40$ とした．

図 3 に衝撃波が 1 つの場合，図 4 に衝撃波が 2 つの場合の時間発展の様子を示す．図 3, 図 4より，どちらの数値解も厳密解と比較して，正しくシミュレーションすることができた．
しかし，図 3 の $u=0.2$ 付近を拡大して見ると (図 5), 振幅が徐々に下がっていくのがわかる．
この誤差をできるだけ小さくするために，空間分割数 , 時間刻み , 窓関数のパラメータ との与える影響を調べた．図 6 より， を変えても大きな変化は見られなかった． については，$dt=0.1$ では差がかなり大きくなり正しくないが， を小さくすれば誤差も小さくなるというわけ
$rrx$

$dt$

$dt$

$nx$$n$$a$$dt$ではなかった (図 7). 窓関数のパラメータは8 ,$)$

$n$

$a$

を変えても誤差の入り方は大きくは変わらないが (図を変えると誤差の入り方に変化が見られた (図 9).

以上より，パラメータを変えることによって誤差を大きく改善することはできなかった．領域$t=100$ までの各時刻のピークの位置 (図 10) での誤差を定量的に調べた結果を図 11 にを拡大して，示す．
$t=100$ で振幅に対する誤差は 0.5%程度で十分小さいと言えるが，この誤差は時刻とともにほぼ線形に成長するのでシミュレーションには注意が必要である．

4 おわりにウィンドウ法による誤差は，十分小さいが，時刻とともに成長するので注意してシミュレーシ

219ョンしなければならない．今後この誤差を小さくしていく必要がある．(10) 式では周期境界条件のもとで髭だけで解いているが，各項を だけではなく も含んだ形に変更したり，dewindowing の方$u$$\tilde{u}$法を変更して誤差の入り方がどうなる力

5$\searrow$

現在検討している．

参考文献“A windowing method for periodic $inflow/outflow$ boundary
treatment of non-periodic flows”, Joumal of Computational Physics 206 (2005) pp.505-535[1] P.Schlatter, N.A.Adams, L.Kleiser,$\overline{\frac{x}{\geqq}}$

2$0$64810X

図1W(x) の概形$(x_{R}=0, x_{L}=10, a=1.18, n=40)$$\tau^{x}\backslash$$\hat{\frac{x}{\tau\geq}}$$0$264810X

図2dW(x)/ぬの概形$(x_{R}=0, x_{L}=10, a=1.18, n=40)$

220$t=0-$$t=10\cdots\cdots\cdots$$t=20\cdot’\cdot’\cdot’\cdot’\cdot t$$\overline{x\check{=}}$

$0$

2

6

48

10$x$

図 3(a)

厳密解による時間発展の様子 (衝撃波が 1 つの場合)$t=0-$$t=10\cdots\cdots\cdots$$t=20\cdot’\cdot|\sim\cdot\cdot$$\overline{\check{\supset}x}$$0$246810X

図 3(b)

数値解による時間発展の様子 (衝撃波が 1 つの場合),

221$t=0-$$t=5\ldots\ldots\ldots$$t=10\cdots\cdots\cdots\cdot\cdot$.

$t=15\prime\prime\prime’\cdot\prime\prime\prime\prime\prime\prime\prime\prime\prime$$\overline{\vee\supset x}$

$0$

2

64

810X

図 $4(a)$

厳密解による時間発展の様子 (衝撃波が 2 つの場合)$t=0-$.

$t=5\ldots\ldots\ldots$$t=10$$t=15$$\hat{\cross\check{3}}$

$0$

2

64

810$x$

図 4(b)

数値解による時間発展の様子 (衝撃波が 2 つの場合)$’\cdots\cdots$$\fallingdotseq\cdot\cdot 1$,

222
0.20005$t=0-$$t=10\cdots\cdots\cdots$0.2$t=20\cdot’\cdot’\cdots’\cdot$

01999501999$\overline{\vee\supset x}$

01998501998

01997501997
0199652

1$0$415678

3

X

図5

図 3(a) の $u=0.2$ 付近の拡大図$\hat{x\check{=}}$

1$0$

2

3

5

46X

図6$t=20$ で $nx$ を変化させた場合の $u(x)$

789,

223

$\hat{\check{3}x}$

2$0$

4

68

10$x$

図7$t=20$ で漉を変化させた場合の $u(x)$$\hat{\equiv x}$$0$123546フX

図8$t=20$ で $a$ を変化させた場合の $u(x)(n=40)$

89

224

$\overline{\vee x\supset}$$0$215437689X

図9$t=20$ で $n$ を変化させた場合の $u(x)(a=1.08)$$t=0-$$t=20\ldots\ldots\ldots$$t=40$$1\cdot\cdot$$t=60\prime\prime\prime\cdot u\prime\prime\prime\prime\prime\prime\cdot\cdot w\prime\prime$$t=80$$t=100$$\hat{x\check{3}}$

$0$

5

10

$|$15X

図 10

各時刻でのピーク位置

20$’.’..$ ,$’\cdot’.’\cdot’$.

225$\hat{*}$$\overline{\dot{\Phi}\rho}$

$0$

20

40

60$t$

図 11

振幅に対する誤差の時刻歴

80100