数理解析研究所講究録
1539 巻 2007 年 201-206201

重合格子法を用いた回転物体周りの二次元流れの数値解析
Numerical Study of two-dimensional flow around aRotating Body using Overset Grid
お茶の水女子大学総合情報処理センター

佐藤祐子 (Yuko Sato)

Information, Media and Education Square, Ochanomizu University
お茶の水女子大学大学院人間文化研究科

水上洋子 (Yoko Mizukami)

Graduate School of Humanities and Science8, Ochanomizu University
お茶の水女子大学大学院人間文化研究科

河村哲也 (Tbtuya Kawamura)

Graduate School of Humanities and Science8, Ochanomizu University

1.

はじめに
近年, 風力エネルギーは環境にやさしく再生可能なエネルギー源の–つとして注目を集めており, 国内外において風力発電量は著しく増加してきている. その–方, 国内では, 台風等の強風時における風車翼損傷問題, 回らない風車問題, 騒音問題等の解決すべき課題も多く残っている. 風車の開発は既に実用段階に入っているが, これらの課題を解決するためには,$\mathrm{C}\mathrm{F}\mathrm{D}$(数値流体力学) による数値実験解析は, 非常に重要であると思われる.

著者らはこれまでに回転軸が風向に対し垂直な垂直軸型風車の数値シミュレーションを行ってきた $[1],[2],[3]$ . この風車は風向の変化による追従制御の必要がなく, 回転速度もあまり速くならないため騒音が少ないという特徴が挙げられる. 計算方法としては, 風車翼の形状に沿った単–の境界適合格子を用いて, 回転座標系の非圧縮性 NavierStokes 方程式を基礎方程式とし, 差分法で解くというものであるが, 回転座標系を単–の計算格子で取り扱うのには, いくつかの問題点があった. まずは, 遠方境界条件の取り扱いである. 回転座標系

では回転角度に応じて、 流入・流出条件を課す位置が移動するため, 流入では–様流, 流出
では自由流出等の条件を与えると, その境目で数値的な振動が起こりやすくなる. また, 遠方境界全体を–様流とする場合は, 十分に広い計算領域を取る必要性があるため, 非常に収束性が悪くなる. その上, 遠方境界の相対速度が非常に大きくなるため数値安定性も悪くなる. 次に, 風車翼が複数になる場合や回転軸の考慮をすると,単–の構造格子の生成が非常に困難になる.
そこで, 風車翼の近傍は今まで通り境界適合の細かい計算格子を用いて回転座標系で解き,風車翼から離れた領域では直交等間隔の粗い計算格子を用いて静止座標系で解くという重合格子法を実装することを考える. 本論文ではまず基礎的な検証を行うため, 実際の風車翼ではなく, 単体の NACA0012 翼を用いて数値シミュレーションを行うこととする.

202

2. 数値計算法
2. 1

基礎方程式

想定する物体の回転速度はあまり速くないため, 流れ場は非圧縮性く avier-Stokes 方程式で表される. 物体周りの流れは二次元流れとし, 物体近傍は回転座標系を用いた以下の方程式で計算を行う.
連続の式:$\frac{\partial U}{\partial X}+\frac{\partial V}{\partial \mathrm{Y}}=0$

運動方程式

,

:
,

$\frac{\partial U}{\partial t}+U\frac{\partial U}{\partial X}+V\frac{\partial U}{\partial \mathrm{Y}}-al^{2}X+2\omega V=-\frac{\partial p}{\partial X}+\frac{1}{Re}(\frac{\partial^{2}U}{\partial X^{2}}+\frac{\partial^{2}U}{\partial \mathrm{Y}^{2}})$$\frac{\partial V}{\partial t}+U\frac{\partial V}{\partial X}+V\frac{\partial V}{\partial \mathrm{Y}}-\omega^{2}\mathrm{Y}-2a\Pi=-\frac{\partial p}{\partial \mathrm{Y}}+\frac{1}{Re}(\frac{\partial^{2}V}{\partial X^{2}}+\frac{\partial^{2}V}{\partial \mathrm{Y}^{2}})$

ここで,$X,\mathrm{Y}$,および U,V は回転座標系での位置と速度, $a$ は物体の回転角速度,

$P$

は圧力,

Re は翼弦長 C および–様流速 uo に基づいたレイノルズ数である. 本論文では物体の回転角速度は–定とみなして計算を行う. これは物体の自励回転ではなく強制回転について計算を行うことを意味する.
また, 物体から遠方の領域では静止座標系を用いた以下の方程式で計算を行う.
連続の式:$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$

運動方程式

,:.

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{\partial p}{\partial x}+\frac{1}{Re}(\frac{\partial^{2}u}{\partial x^{?}\sim}+\frac{\partial^{2}u}{\phi\underline’})$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\phi}=-\frac{\partial p}{\partial y}+\frac{1}{Re}(\frac{\partial^{2}v}{\partial x^{2}}+\frac{\partial^{2}v}{\partial y^{2}})$

ここで,,x, y および u, 嫁よ静止座標系での位置と速度である. 回転座標系と静止座標系における位置と速度の関係は, 回転角度を$x=X\cos\theta+y\sin\theta$,$\theta$とすると, 以下の式で表される.

$X=x\cos\theta-y\sin\theta$

,

,

,

$y=-X\sin\theta+\mathrm{Y}\cos\theta$$\mathrm{Y}=x\sin\theta+y\cos\theta$$u=U\cos\theta+V\sin\theta+a\nu$ ,$U=u\cos\theta-v\sin\theta-a)\mathrm{Y}$$\nu=-U\sin\theta+V\cos\theta-a\alpha$,,$V=u\sin\theta+v\cos\theta+\alpha X$.

回転座標系で用いる計算格子と静止座標系で用いる計算格子の重ね合った部分では, 上記の関係式を用いて物理量の交換を行う.

203
これらの基礎方程式を, Fractional Step 法を用いて差分法により解く. 以下に回転座標系での運動方程式を用いて Fractional Step 法の流れを簡単に示す.
まず, 運動方程式の圧力項以外を用いて, 仮の速度 VIE を求める.
$\mathrm{V}^{k}=\mathrm{V}^{n}+\Delta t\{-(\mathrm{V}^{\prime l}\cdot\nabla)\mathrm{V}^{;\mathrm{t}}+\frac{1}{Re}\nabla^{2}\mathrm{V}^{i1}-\omega\cross(\omega\cross \mathrm{r})-2\omega\cross \mathrm{V}^{n}\}$

次の時間ステップ $n+1$ において連続の式が成り立つとすると,を求めることができる.

以下の式により圧力〆+1$\nabla^{2}p^{n+1}=\frac{\nabla\cdot \mathrm{V}^{l}}{\Delta t}$

最後に, 仮の速度 V* および圧力 p“+l を用いて, 次の時間ステップの速度 V731 を求める.
$\mathrm{V}^{n+1}=\mathrm{V}^{\cdot}-\Delta t\nabla p^{n+1}$

ここで, \Delta t は時間刻み幅,$\mathrm{r}=(X,\mathrm{Y})$, V=(U,V) とする

圧力の PoisSOn 方程式は, SOR法を用いて解く. 非線形項は 3 次精度の風上差分を用いて近似し, その他の空間微分は 2 次精度中心差分, 時間微分は 1 次精度前進差分を用いる.

2. 2 計算格子
NACA0012 翼の周りに$\mathrm{O}$型格子を生成する. これを内部格子とする. また, 周辺領域には直交格子を生成し, これを外部格子とする. これらを重ね合わせたものが全体の計算格子となる. 格子点数は内部格子 $96\cross 61$ , 外部格子 $199\cross 81$ である. 図 1 に計算格子を示す.

(a) 物体 (NACAOOI 2 翼) 近傍の計算格子(b) 計算格子全体
図1

計算格子

204
2. 3:

境界条件

流入は–様流とし, 物体表面は, 回転座標系においてすべりなしとする. 流出は速度, 圧力ともに外挿とする.

$-\circ D$$\mathrm{I}^{26\mathrm{C}}$$\mathrm{y}\llcorner \text{一}$$-\mathrm{c}:\text{翼弦長}$

図2

計算モデル

回転座標系で用いる内部格子と静止座標系で用いる外部格子の重ね合った部分では, 物理量の補間を行い, それぞれ相手の物理量を自分の境界条件として用いることとする. 補間式は,以下を用いる.
$f_{(1}=.

\frac{1}{1_{f}’r_{1}+1’\prime r_{2}\prime+1^{l},’’\prime r_{\tau}+1r_{4}}.(\frac{1}{r_{1}}f_{1}+\frac{1}{r_{2}}f_{2}+\frac{1}{r_{:}}f_{!}+\frac{1}{r_{4}}f_{4})$図3

格子間の物理量補間

3. 計算結果
レイノルズ数は翼弦長に基づき 10000 に固定して計算を行った. 翼は回転することにより迎角が変化することとなる. 図 4 は, 迎角が$0,5,10,20,2\tilde{\mathrm{o}}^{\mathrm{o}}$における無次元時間$\mathrm{T}=30$のときの圧力のシェーディングおよび等高線と, 速度ベクトル図を示す. 格子が重ね合った部

分の流れ場も概ね滑らかに解くことができている.

$\cdot:.\backslash \backslash$$::::::$$\oplus$:.
:.
$\mathrm{w}’[]\cdot:\mathrm{s}.\infty 000$$.u[]:\uparrow \mathfrak{g}0000$(a)$\alpha=0$

Deg.

205

(b)
ffld [\S g$\alpha=5$

Deg.
$arrow=$$\kappa 8!r_{\wedge’}’$$-arrow-arrow\simarrow\sim \mathfrak{g}r\vec{\sigma}\mathrm{t}99\vec{a}1\alpha\vec{c}_{r}^{-_{--}arrow--\vee^{-==}}arrowarrow-\wedge\equivarrowarrow’\sim_{<}’arrow\congarrowarrowarrow-\wedge\sim\approx\overline{\overline{\sim}}\underline{\sim}$$arrow\simarrow-arrowarrow-=’\sim-\wedge\prime\prime-arrowarrow-arrowarrow-arrow-arrow-\wedge$$\sim-arrowarrowarrow-arrow----arrow\approx--\frac{arrow-\prime}{arrow}\vee^{-\wedge^{-arrow\wedge}}\simarrowarrow-arrowarrow\frac{\prime_{arrow}\overline--arrow---arrowarrow}{}-\vee^{-\wedge-}\approx\frac{\approxarrow}{\simeq-}-\sim’\sim-\sim-arrow-\simeq-\wedgearrow’\simeqarrowarrow\sim\simarrow=_{arrow\sim}\wedge\wedgearrow-’arrow\sim$

図4

(c)

$\alpha=10$

Deg.

(d)

$\alpha=20$

Deg.

(e)

$\alpha=25$

Deg.

翼周りの流れ場$(\mathrm{T}=30)$

206
4.

まとめ
本研究では, 重合格子を用いて, NACA0012 翼周りの流れ場の数値シミュレーションを行った.

格子の重なり合った部分における圧力および速度の物理量の交換は, 概ね良好であるが, 翼端からの放出渦等の物理量の変化が大きい状態の移流を捉えることや, 回転速度が大きくなったときの数値振動に関しては,まだ問題が残されている.

今後は, 物体に働く揚力, 抗力, 回転力について詳細に評価検討を行う予定である. また, 実際の風車翼に合わせ, 複数の翼を配置した格子で計算を行っていきたい.
なお, 本研究の–部は科学研究費補助金 (基盤研究 (B)(2)) (16360478) の補助金を受けて行った.

参考文献[11 T.Kawamura, A.Shinohara, Y.Sato, M.Kan, \dagger ’Numerical study

dynamics around a(2002),$\mathrm{c}\mathrm{r}\mathrm{o}\mathrm{s}\mathrm{s}\cdot \mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{u}\mathrm{l}\mathrm{b}\mathrm{i}\mathrm{n}\mathrm{e}^{\prime 1}$of the flow andTheoretical and Applied Mechanics Vol.51,$\mathrm{p}\mathrm{p}231\cdot 240$

[21 Y.Sato, T.Kawamura “Numerical study of the blade dynamics for a$\mathrm{c}\mathrm{r}\mathrm{o}\mathrm{s}\mathrm{s}\cdot \mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}$turbine”

ACFD5 Proceedings, (2003)[31 Y.Sato, T.Kawamura, “Numerical study of the flow around a high torque vertical axis
wind turbine”, WWEC2004 abstracts, (2004), pp161