数理解析研究所講究録
第 80 巻 1970 年 45-53

Burgers Model の数値解

細川巖

航技研

山本稀義

Burgers model の一様な乱れの相関量を求めるために行なわれた数値計算について報告する。

乱流の相関量の計算は既に報告された様に 1) 特性汎数の初期値問題の計算に帰着するが, この中に含まれる汎関数積分のモンテカルロ法による数値計算を行い, 乱れのエネルギー, 相関々数澄よびエネ

ルギースペクトルの時間的経過を求めた。

91.

$\cdot$Burge rs mode 1の相関量

Burgers 方程式は無次元形て(1.1)$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\frac{1}{R}\frac{\partial^{2}u}{\partial x^{2}}$と与えられる。 この場合 $u,$$x_{0},$$t_{0}=x_{0}/u_{0}$と Cole3)$x,$$t$の単位はそれぞれ初期の速度の自乗平均である。従って $R=u_{0}x_{0}/\nu$ ($\nu$$u_{0}$, 乱流の初期スケール: 動粘性係数) である。上式については Ho pf$\emptyset$

によって初期値問題の解が次の様に与えられている。

$u(x, t)= \frac{\int_{arrow\infty}^{\infty}(x-x’)exp\{-\frac{R}{2}\int_{0}^{x’}u(x’’,0)dx^{\nu}-\frac{R(x-x^{J})^{2}}{4t}\}dx’}{t\int_{arrow}\infty exp\{-\frac{R}{2}\int_{0}u(\mathscr{K},0)dx^{v}-\frac{R(x-x)}{4t}\}dx,\infty x’2}$(1.2)

われわれはここで, 初期の乱れは一様で, 相関々数 $Q$ (拗が Gauss 型 $Q(x)=e^{arrow x^{2}}$ で与えられる場合を考える。 この時 $u(x, 0)$ を次の様なフーリェ級数形$u(x, 0)= \sum^{\infty}\sqrt{\sigma(k_{n})}\underline{1}$($\Lambda_{n}$co $sk_{n}x+M_{n}sink_{n}x$ )$\sqrt{\pi}n=1$$k_{n}=(n- \frac{1}{2})\triangle k$

で,$\sigma(k)$は–45-$\triangle k$(1.3)

(1.4)$\sigma(k)=\int_{arrow\infty}Q(x)e^{arrow ikx}dx\infty$で表わすと, 時刻$t$の相関々数は次の様になる事が示される 1)。

$<u(x_{1}, t)u(x_{2}, t)>= \int_{-\infty}\cdot\cdot\int_{-\infty}u(x_{1}\infty\cdot\cdot\infty, t)u(x_{2}, t)\prod_{n=1}^{\infty}exp\{-\frac{1}{2}(\Lambda_{n}^{2}+M_{n}^{2})\triangle k\}$(1.5)$\cross d(\Lambda_{n}\tau/^{1\frac{\triangle k}{2\pi}}\ulcorner)d(M_{n}\sqrt{\frac{\triangle k}{2\pi}})$

integrand の $u(x, t)$ は (13) によ. り$x$と$t$をパラメータとする{\Lambda n$\sqrt{\triangle k},$$M_{n}\sqrt{\triangle k}$} の関数である。上式の無限回繰返し積分はもちろん実行困難な \sim ので, 実際はフーリェ級数形を有限の $N$ 項($\triangle k$も有限) で打切って, 本質的には 2 $N$ 次元の積分で近似する。われわれはこの式を用いて, 乱れのエネルギー, 相関々数の数値計算を行なうが, その計算方法を次に述べる。

\S 2.

数値計算の方法

第一に, $N$ 有限の近似が $Narrow\infty$ で極限値を能率よく表わすかどうかが問題である。 われわれは $N=5$から始めて, $N=15$ の程度で収束の徴候を見た。 この問題は現在の数値計算法の根本問題であるので,もう少し深く検討を要するが, 現報告では $N=15$ および $N=30$ のデータを中心にまとめた。
$N$

有限であっても (1.5) は高度な多重積分であるから, モンテカルロ法を使うのが便利である。 この$t$場合, (1.5) の形から見て標準正知分布に従う確率変数として $\{\Lambda_{n}\sqrt{\triangle k}, M_{n}\sqrt{\triangle k}.\}$ を採用するど,estimator4) は $u^{N}(x_{1}, t)u^{N}(x_{2}, t)$ となる。 ($\text{事^{}\backslash }$を示す。 ) 換言すれば,

2 $N$ 個の標準正規乱数によって$u^{N}$は (13) を有限フーリェ級数で近似した$\{\Lambda_{n}\sqrt{\triangle k}, M_{n}\sqrt{\triangle k}\}$の組を作り, これから $u^{N}(x_{1}, t)u^{N}(x_{2}, t)$ を評価し, これを何度も繰返して estimator の多くのサン 7. ルについて平均をとれば, (15) の近似が得られるというわけである。 モンテカルロ法の精度は積分の次元に無関係に, サンプルの個数によって決まるので甚だ都合が良い。 どれ位のサンプルをとれば良い結果を得られるかについては,\S 5. で実例をもって示される。

標準正規乱数の発生についてはいくつかの方法があるが, ここでは比較的簡単でかつ良い精度を与える次のものを採用した$5)_{\circ}$

ま す, 二つの異る系列の一様乱数の対が次の方法で作られる。

$U_{\iota i+1},=U_{1,i}(2^{12}+1)+1,$$(mod2^{35})$$-46-$(2.1)

$(mod2^{35})$$U_{2,i+1}=(U_{2}, i(2^{7}+1)+0.788\cross 2^{35})(2^{7}+1)+0.788\cross 2^{35},$

ここで最初の出発値としては $U_{1,1}=233362477003,$ $U_{2,1}=212312312323i\partial^{i}$ とられる。 こ

一様乱数の対は最大値$\circ$$2^{35}$$[0,1]$ に分布する一様乱数で規格されて,$U_{1},$$U_{2}$に直されたのち標準正規乱数に変換されるが, その変換には次の 2 通りの方法が有る。
$X_{1}=(-210gU_{1})^{\frac{1}{2}}cos2\pi U_{2}$(22)$X_{2}=(-2logU_{1}^{\frac{1}{)^{2}}}sin2^{\pi}U_{2}$

これらはさて,$U_{1}$$U_{2}$から同時に作られるので, 計算の際には $X_{1},$$X_{2}$が交互に採用された。

estimator の計算問題に移ろう。 まず(2.3)$u^{N}(x, 0)= \frac{1}{\pi^{1}/4n}\sum_{=1}^{N}earrow k_{n}^{2}/8(\Lambda_{n}cosk_{n}x+M_{n}sink_{n^{X)\triangle}}k$に寿いて\triangle $k$ を決める必要がある。 $N\triangle k=k_{c}$ は $u$ (紛を構成する波数成分の考慮され得る最大波数を示すが,$k_{c}$が余り小さいと, エネル羊- の相当部分を含む波を無視してしまうことになり, 物理的に考えて, 好結果は期待できない。無視する部分を全体の 1%程度に止めるには(2.4)$4_{c}^{\infty_{e}k^{2}/8}arrow dk/f_{0}^{o}earrow k^{2}/8dk=1/100$という方程式を解けばよいだろう。結果は$k_{c}\fallingdotseq 5.17$である。 これから$\triangle k$が $N$ に依存して確定することになる。
$\int_{0}^{x’}u^{N}(x’’, 0)dx’’$は解析的に言「算できて$\int_{0}u^{N}(d’, 0)dx’’=x’\frac{1}{\pi^{1}/4}\sum_{n=1}^{N}\frac{e^{arrow\frac{k_{n}^{2}}{8}}}{k_{n}}$となるが, 必ずしも無限小でない$\triangle k$( $\Lambda_{n}sink_{n}x’-M_{n}c$ os $k_{n}x’+M_{n}$ )$\triangle k$(25)を扱う今の場合には, これをむしろ$k^{2}$$\int_{0}^{x’}u^{N}(x’’, 0)dx’’=\frac{1}{\pi^{1}/4}\sum_{n=1}^{N}\int_{n^{-\frac{1}{2}\triangle^{\triangle}}}^{k_{n}+\frac{1}{2}}kkk_{\frac{e^{-}}{k}}^{-}\tau${$\Lambda(k)sinkx’-M(k)$co $skx’+M(k)$ } $dk$$\fallingdotseq\frac{1}{\pi^{1}/4}\sum_{r\approx 1}^{N}\frac{earrow k_{n}^{2}/}{k_{n}}8\{(\Lambda_{n}sink_{n}x’-M_{n}cosk_{n}d)\frac{sin\triangle kx\prime/2}{x\prime/2}+M_{n}\triangle k\}$という表現に代えた方が精度が良い。実際 (25) の表現を採ると,$x$が大きくなると estimatorのサンプルが甚だしく振動して, モンテカノレロ法の能率を著しく悪くさせることが分った。

-47-

(2.6)

(2.6) の$\sin\triangle kx/2/x/2$ の因子はこの振動を抑える働きをする。

(2.6) を (1.2) に入れ, 分母分子をシンプソン則で積分すれば $u(x, t)$ が計算できる。 シンプソン則による積分には次の二つの問題がある。一つは, 上限下限を精度を損わないように有限に止めることであり, 他の一つは,integrand 自体が指数関数的増加 ($e^{\pm R}$の程度であるから $R>10$ だと大変なことになる。/) を行うので, 適当な規格化による抑制が必要になることである。 われわれは次のような方法を取った。要点だけ述べると, かなり広い $X’$ の範闘で,いての何点か (2 $0\sim 500$ 点) につ

$x’$.

.

integrand の $exp$ の肩の関数を試みに計算し, その中の最大値をさがす。次に $exp$ の肩の関数からこの最大値を差し引くことによって inte.grand の規格化を行い, この最大値を与える点から $x’$の両側に向って, 予め決められたステップ間隔でシンプソン則の積分を計算し始める。そしてgrand が$10^{arrow 7}$って, 両遠方で,以下になるとそこで積分を止める。 Gauss の因子$e$inte-xp $\{-R(x-x’)^{2}/4’t\}$ によintegrand が単調に減衰することが分っているから, 上のようにして積分範囲を決めることは妥当と思われる。
次節で多くの計算結果をグラフによって示し , 付論する。

\S 3.

結果と討論

乱れのエネルギー $<u(0, t)^{2}>$ の計算結果を第 1 図に示す。横軸は無次元時間$t$で,縦軸はエネル

ギーの 2 倍である。計算は $N=15$ で, $R$ は 0.1,1 , 10, 100,1000 について行なわれた。 破線は(1. 1) の非線型項を零とした線型理論の結果で右側は $R=0.1$, 右側は $R=1$ の場合である。実線はMeecham & Siegel の Wiener–Hermite 展開法からの結果 $(R=100)$ であるから, $R$ が小さい場合 $(R=0.1,1)$ は,$t$$6)_{\circ}$

図が小さい所で線型理論と一致するが, t が大きくなると線型理論よりも早くエネルギ 7 は減衰する事がわかる。 これは非線型項による低波数領域より高波数領域へのエネルギー輸送の結果と考えられる。
$R$$R$が大きい場合 $(R=10,100,1000)$ は線型理論と著るしく違った特徴を示している。すなわち,を大きくして行った時に, エネルギー減衰曲線は, $R=\infty$ と置いて期待されるエネルギー減衰のない状態に近ずかずに, ある $=—-$ クな一つの減衰曲線に漸近する傾向を示している。 これは$Rarrow\infty$で有限なエネルギー散逸がなければ可能なことではない。 この説明としては次の事が考えられる。 Bu rge rs.

–48–

方程式 (1.1) の $R$ の大きい場合の漸近解として Burgers によって導かれた 7) 主角型の shock f rontの相一るが, 一つの shock

front の全散逸エネ

$Js$

ギー

$\frac{1}{R}\int(\frac{\partial u}{\partial x})^{2}d\alpha$は $R$ に無関係な有限値となる事がまずわかる。 さらに, このような相似解の多数から作られるアンサンブルを考えても,全散逸エネルギーは $R$ に関係する事がないから, この場合 $R$ の如何にかかわらず常に一定のエネルギー散逸が起るはずである。つまり$Rarrow\infty$で漸近的なエネルギー散逸の量が存在するはずである。 この事からエネルギー減衰のある漸近曲線が得られる事は不自然ではない。 われわれの場合, 解を Bu $rgers$の相似解に限定したわけではないので, 得られた結果はこの主張が広く一般的に成立する事を示している。
また, Meecham&Siegelの結果との比較については $t=3$ 以上になると大きく違ってくる事がわかる。

第 2 図にエネルギーの計算 $(N=15,R=100)$ のモンテカルロ法の安定性の例を示す。横軸はモンテカルロ法の試行回数で, ほぼ 300 回の試行回数でほとんど一定値に近ずいている事がわかる。

第 3 図に相関々数の計算結果を示す。横軸は無次元の距離ある。 実線は初期の相関$-X^{2}$$e$$x$で, 縦軸は$<u(0, t)u(x, t)>$ での曲線を示し, 計算は $N=30$ で行なわれた。 $R=100$ の場合$t$は1と 3 の計算が行なわれた。 この図は時間の経過とともにエネルギー ( $x=0$ の値) が減衰し, 相関が広がっていく事を示している。また $R$ の違いによる相関々数の比較のため $R=1,$ $t=1$ の計算が行なわれた。 この場合のエネルギーの絶対量はほぼ $R=100,$ $t=3$ の場合と同じであるが, 明らかに $x=$$0$の近傍の勾配は違っており, この結果はエネルギースペクトルの高波数成分の相違となって現われている。 (第 5 図) 破線は Meecham& Siegel の $R=100,$ $l=1$ の結果で, われわれの結果と構造上の違いがはっきりみられる。

相関々数の計算におけるモンテカルロ法の安定性の一例を第 4 図に示す。 これは第 3 図の $R=100$ ,$t=1$ の場合で, 横軸は試行回数である。

第 3 図で計算された各点を曲線で結び, フーリェ変換して得られるエネルギースペクトルを第 5 図に
/7-\す。横軸は波数$k$

は Meecham&で,縦軸はエネルギースペクトル $E$ (k) である。実線は $t=0$ のスペクトルで, 破線Siegel の結果である。第 3 図から計算される値はいくぶんばらつきが有るので,図では代表的な点を示してある。 この図からわかる事は次の通りである。 $R=100$ の場合, 初期の相-49-関$e^{arrow x^{2}}$に対応するいわゆる釣鐘型のスペクトルからも出発しても, $t=1$ ではエネルギースペクトルに既に変化が起り, スペクトルの高波数領域にかなりの範囲にわたって$k^{arrow 2}$領域が見られる。 さらに$t$$=3$ では, 全エネルギーの絶対量が小さくなるが, スペクトル曲線にはまだかなりこの領域が保たれている。従って$k^{arrow 2}-$スペクトルを予言する),
inertial subrange の理論$6$$8$),$9$), $10$ )はここで実験的裏付けを得たと言えるだろう。 $R=-100$ 0-で $t=1$ の場合の計算も行なわれたが, 図 1 の結果から期待されるように, $R=100,$ $t=1$ と重なる結果が得られた。 従って, 示されているでは $R=100$$k^{arrow 2}$

領域ぱ$J$,$k$の範囲$t=1$ のスペクトルは $Rarrow\infty$ の漸近スペクトルと考えてもよい。

en 等 11) の数値実験の結果でもみられるが, われわれの場合とは定量的には一致して$g$いない。彼等の数値実験はやはり Hopf と $Co1e$ の厳密解を使い, 乱数によって初期条件のアンサンフルを作り出している点で, われわれの計算と似た所が有る。 しかし彼等のアンサンブルの要素の数(モンテカルロ法の試行回数に相当) は”severa1” であるからわれわれのものと比べて大きな違いが有り, 統計的性質が正しく導かれているかどうか大いに疑問が有る。 (この欠点を禰うために彼等はspace average という別の操作を導入している。 ) また, 初期条件のアンサンブルは我々のもの(正規型 1)) とかなり違っているので, 結果に違いが有っても不思議はない。
$R=1.0,$ $t=1$ のスペクトルは $R=100,$$t=3$ の場合に比較して初期のスペクトルに近い型をしている。両方とも初期のエネルギーからほぼ同程度のエネルギーを失った状態であるが, 明らかに後者の方が非線型項による低波数領域から高波数領域へのエネルギー輸送が活発に行なわれている事を示している。
またエネルギースペクトルの $k=0$ の点のエネルギーは保存すると考えられるが (Loitsianski定数に対応する$6)_{o}$), 計算された $k=0$ の値は 5%以内で一致し, モンテカルロ法の精度と同程度である。

\S 4.

結び

以上は, 特性汎関数の初期値問題の解としての汎関数積分を直接に扱うことにょって, 乱れを魏直的に追求した始めての試みとして興味が有るかも知れない。汎関数積分という概念は Wiener 積分以外
はいまだ数学者に広く認められておらず, 著者によって形式的に誘導された諸公式が, 本当に意味のあ–50-あるものかどうか疑問視する向きもなきにしも有らずであるが, 現在の計算によって, 少くとも Bur-gers model に対しては, 以上のような物理的に意味のある具体的な結果が出たからには, われわれの形式解が基本的に意味のあるものである事を疑う理由はほとんどなくなったといえるであろう。 もちろん, モンテカルロ法のために, 精度は場合によって十分とは言いがたいが, この方の欠陥はモンテカルロ法の技術そのものの refinement によって克服されることが期待される$4)_{\circ}$

Burgers model の乱れの物理的性質については, 前節で詳しく討論したが, この研究によって明らかになった新事実は,$Rarrow\infty$に於ける乱れの挙動についてであった。 $R=1$ を境として, 乱れのエネルギーの減衰のパターンは明らかに変り, $R=10$ で既にあるユニークな減衰曲線に収束する傾向を見せてくる。 Burgers 方程式で $R=\infty$ として得られる無減衰の解のアンサンブルと上の漸近的乱れのアンサンブルとの間には, 大きい不連続の溝があることが判明した。 Saw wave 展開にょって Bur-gers model の乱れを取扱う巽 10) の議論も定量的な面はとも角,$Rarrow\infty$におけるエネルギーの減衰の漸近的行動を与えている事は興味深い。 これは Burgers mod el だけの特徴であるか, もっと一般的な性質であるか, 予断は出来ないようである。漸近曲線について $N=30$ の場合の計算を目下検討しつつあるので, これが終ればなお確信のある結論が出されると思われる。

文献
1)

細川

巖: 数理解析研究所講究録47 (1968) 46, I. Hosokawa.:

Phys. Fluids 11 (1968)2052.
2)

E. Hopf: Commun. Pure Appl. Math. 3 (1950) 201.

3)

J. D. Cole: Quart. Appl. Math. 9 (1951) 225.

4)J. M. Hammersley and D. C. Handscomb: Monte Carlo Methods, (John Wiley&Sons,Inc. 1966).
5)

R. Kronmal: J. Assoc. Comput. Machinery 11 (1964) 357.

6)

W. C. Meecham and A. Siegel: Phys. Fluids 7 (1964) 1178.

7)

J. M. Burgers: Advances in Applied Mechanics 1 (1948) 171 (Academic Press).

8)

J. M. Burgers: Proc. Acad. Sci. Amsterdam 53247.

9)

R. H. Kraichnan: Phys. Fluids 11 (1968) 265.

10)

巽11)

D. T. Jeng, R. Foerster, S. Haaland and W. C. Meecham: Phys. Fluids 9 (1966) 2114.

友正: 数理解析研究所講究録47 (1968) 57.

-51-

臓
$\circ\iota$

図$H$$\triangleleft_{/}i$$\sim$$*$$1$$\aleph$$\rangle_{0}$$Q$ア$-S2-$

$f$.第1図

乱れのエネルギー減衰$\triangle R=l()()()$.

$R=J()(),$$(1v=15)$$\cross R=l()$$M_{(((}\cdot/l\iota m((’\cdot\backslash \{)^{\backslash }j_{((/(}’.’/,$. ▲ $R=l$ ,$—-$$\bullet$$R=0.1$ ,線 fFf!-J」里論

第 3 図相関々数

第2図-

モンテカルロ法の安定性エネルギーの計算 $(R=lm$$-S3-$$N=30$ )