スーパーコンピュータを用いた平行平板間流れ場の圧力に関する解析
下向 建秀 ££Ý柳谷 貴之 £大泉 健治 Ý神田 英貞 £

会津大学・コンピュータ理工

東北大学情報シナジーセンター

概要
平行平板間流れ場の圧力に関する子細な解析を試みた。 本論文では、 レ イノルズ  数による変化及び格子変化による影響について述べる。 精密な数値計算により、 境界層方程式、ベルヌーイの法則とは若干異なる結果を得た。 その差異とは、 流れの発展初期の段階においては、 壁側の圧力は中心線上のそれよりも低い値を示すというものである。 ただし 、 この結果は  数が  以下で顕著に現れている。 格子変化による影響は 、 全体ではさほど の差異は見られなかった。 なお、 この計算は流れ関数と渦度表示の方程式に対し 、 差分法で行った。 この方程式系の特徴は 、 変数の少なさ及び連続の式を厳密に満足する点にある。 しかし 、 格子幅を細密にとり、 また、 反復法を用いて解くため、 計算時間が長いという問題点がある。 高いベクトル演算機能をもつスーパーコンピュータを用い、 また、 計算アルゴ リズ
ムを見直すことにより、 計算時間の削減を試みた。



レ イノルズ数  ¾ ¼ 
定常状態までの助走区間の長さ無次元での定常状態までの助走区間の長さ
平行平板間の幅
方向の速度方向の速度

時間刻み
圧力入口から  方向への圧力の落ち込み 流れ関数
渦度   はじめにはじめに、 解析条件について述べる。 平行平板間の非圧縮性流体解析の体系化はすでに理論上において、 また、 数値計算においてもなされている。 本研究の主目的は、 二次元平行平板の入り口で、 一様流度分布が与えられた場合の、 流れ場における圧力分布に関して高精度の解析を行うことである。 特に、 流れ方向に垂直な断面上での圧力分布について研究する。
元来、 平行平板の壁側と中心線上での圧力分布に関しての議論では、 ベルヌーイの法則によると壁側の圧力の方が中心線上のそれよりも高いとされ 、 また、 境界層方程式を用いた解析によると、 壁側と中心線上での圧力に差は無いとされている。 本研究では、 近年発達の著し
い高速計算機を用い、 また、 解析方法としては、 完全な形の 

方程式により、解析を行う。

実際の計算としては、 まず、 平行平板内の流れについて、 入口で時刻    以降一様な流速で入る流れの定常流への移行を計算した。 以下に図  を示す。 入口で  方向への速度を ¼と与え、方向の速度は  とする。

Initial uniform
velocity profile

Channel wall

Fully developed
parabolic profile

hyx0Centerline

Entrance length(xe)Channel wall

Dimensionless entrance length(Le=xe/(2h Re))図      平行平板間流れにおいて、 入口での一様な速度は、 放物型に増加する。 これは、 壁面付近においての粘性力の働きによる。 圧力に関しては、  方向、 方向の両方の解析を行った。
格子間隔を小さくとることにより、 壁側付近と中心線上での圧力分布の微細な変化の解析を行った。

 基礎方程式
非圧縮性流体の運動を支配する方程式は 、 質量保存を表す連続の方程式  及び運動量の保存を表すナビエーストークス方程式  である。

     ¾ 
この時、  は圧力であり、  はレ イノルズ数である。
式の圧力項の消去及び 、 ベク







トル解析の公式よりを得る。 ただし 、   



は渦度で、
との関係式及び渦度輸送方程式

が得られる。 次に、 % 式を
·½#$    ¾ 

%

に関して、 陰的に解くオイラー陰解法を考えると、

     ·½     ·½&!

" 

¾¾である。 ! 式及び 、 速度  、  と流れ関数

より、 流れ関数ー渦度の方程式







 

 ¾ ·½  ¾ ·½  ¾
¾





'となる。 圧力分布は、 以下で与えられる ()。

¾   ¾   ¾ 



 "

¾



¾

 ¾*

流れ関数+渦度の計算及び 、 圧力分布の計算の差分方程式は、 双方ともヤコビ法を用いて解かれた。 また、 壁側の渦度は３点から１階微分を決めるので、 最も精度が良くなる様に選択
すると 、 以下の式となる。

        !

"&   "        

& ¾

同様に、 圧力分布に関しても、 ２階の精度を得られるように近似式を選択した。





InletChannel wall

OutletFlowJ0∆yh∆xj1Centerline1i

図  , i+1I0 -.  格子生成及び計算手順
生成した格子の概略図は図  に示す。 プログラムの実装は、 図 ! のフローチャートに沿った。 ここで、   、   は、 それぞれ  方向、方向の格子点の最大数である。 計算には、

東北大学情報シナジーセンター所有の /0 1+% スーパーコンピュータを使用した。 なお、計算は４節、５節で述べるように十分な高速化チューニングの後、 行った。

StartStart

Grid generation

Grid generation

Get initial value
ψ and ωGet initial value
ψ and ωSolve velocity
distribution
u and vset loop
Lp1

Solve velocity
distribution
u and vEq.(5)

Eq.(5)

Solve for
ω(n+1) from ω(n) Eq.(8)Solve for
ω(n+1) from ω(n) Eq.(8)Iteration methodsfor ψ(n+1) from Eq.(6)ψ(n) and ω(n+1)

Iteration methodsfor ψ(n+1) from Eq.(6)ψ(n) and ω(n+1)

Check| ψ(n+1)-ψ(n) |<εCheck| ψ(n+1)-ψ(n) |<εUpdate and check| ω(n+1)-ω(n) |< ε2Update and check| ω(n+1)-ω(n) |< ε2Solve pressure
distribution pEq.(9)Solve pressure
distribution pOutput
Output

改良前
図 ! 2.   3

-





改良後

Eq.(9)表  改良以前の計算結果
4 5  !"%$!%9!#*",67

,2867

8/

9 6

9  :

'%'9#'#'

!"'9#!!

"*9$'*"*9%"%"*;  %!#"$9%"'$!

表  改良  計算アルゴ リズムを見直した結果
4 5  

'"#'9$'"#%,67

,2867

8/

9 6

9  :

!9'#%'!

$#9$#"$*!

"%9#$#"*9!$'#!

;  "%#9*'$

高速化の手順プログラミングは 2  * を用いて行われ 、 コンパイラは 3* を用いた。 また、 計算モデルとして、 格子数を  方向に  、  方向に # とした平行平板流れ場を考える。 高速化の手順としては、

9計算アルゴ リズムの見直し9

ベクトル演算率の向上及び平均ベクトル長の改善!9バンクコンフリクト時間の削減を行った。

 計算アルゴリズムの見直し
計算アルゴ リズムは、 図 ! のフローチャートの  から - へと変更した。  では時間刻みを考慮し 、 全体のループ回数を指定出来る。 しかし 、 充分に誤差判定を満足しても計算を続けるため、 効率が悪い。 - では、 誤差判定基準を満足するまで時間発展させれば良いため、 計算効率は改善される。 なお、 同時にコンパイルオプション +7  の指定による自動並列化を行い、 "074 で実行した。 結果は、 表  から表  となった。 平均ベクトル長の増加及びバンクコンフリクトの削減により、 計算時間は短縮されたものの、 ベクトル演算率が極端に低下した。

 ベクト ル演算率の向上及び平均ベクト ル長の改善
まず、 多重ループにおいて配列データへのアクセスのアドレスが連続になる様に、 ループの入れ換えを行った。 この改良は、 キャッシュラインの有効利用を目的として施される。 結果として、データ参照が効率化され 、 メモリアクセスのオーバーヘッドが減少し 、 計算の高速化が期待できる。 図  で示している様に、 本プログラムでは   の最大値     は既知であり、 常に      であるため、 以下の例の様に入れ換えを行った。 プログラム全体では、 ! 個所が改良された。

              

次に 、 漸化式のマクロ演算により、 ベクトル演算率の向上を目指した。 下に実際改良を加えた部分を示す。

                                     と書かれた部分を 、 以下の形にする。 また 、 あらかじめ計算できる部分はループ 外で前もって計算し 、演算数（ 量）を減らした（  、 、 ）。
1  1+ <式

式 の形式は漸化式のマクロ演算にあたる。

                                        

この改良により、 表 ! の結果が示された。 ベクトル演算率が向上し 、 計算時間に大幅な改善がみられたものの、 バンクコンフリクト時間が大きく、 十分な高速化が達成されたことにはならない。 しかし 、 さらにいくつかの配列宣言の見直しとサブオプションにより、 表 "の結果にまで削減できた。

表 ! 改良  ベクトル演算率と平均ベクトル長の向上
4 5  '"#!$9!#"'#,67

,2867

8/

9 6

9  :

"!#9'*!%$9!"!

"%9##"*%"''9##*!

;  #%!9### まとめ
施したチューニング全体では、 計算時間はおよそ =# 程度まで削減できた。 図 " にその推移を示す。

CPU Times
35000
30000
25000
20000
15000
10000
5000

12

3
Tuning Steps

図 " 高速化チューニングと計算時間4表 " 改良 ! バンクコンフリクト時間を削減した結果
4 5  !"'9$#%$

,67

,2867

8/

9 6

9  :

"""9'$'#

!9!""%#!

"*9""'*"***9'$#;  
!9%!'計算結果
計算精度の比較として、 無次元での定常状態までの長さである  より求められる、 平板間の中心線上の速度達成率 *: 、 *': 、 **: 、 ***: の地点までの距離を表 # に示す。 なお、 これらの数値に関して、 先行研究との比較も同時に示す。 また、 である。

次に、 図の説明をする。 以下、 図 # から図 ' は、 格子数を  方向に  、  方向に  、

#、  、 # と変化させた場合の壁側の圧力の落ちこみ 7   を示している。 また、  数は ' 、  、 $ 、  とした。 なお、 無次元での格子幅を & &し、 と  の範囲で解析を行った。 格子の変化による圧力分布の微細な変化を確認する
  の範囲で表示する。 計算時間 新プログラム、 にて実行 は表 に示ため、   

1+%

$す。 また、 旧プログラム 1+" にて実行 の計算時間を表 % に示す。 格子数の関係により、  、     とした計算（ 表中< ）にて比較できる。 結果は!!'%    '!

!!

1+"



'!*

1+%となった。 実行環境が 1+" から 1+% へと変更されているものの、 アルゴ リズムの改善により、 約 *# 倍の計算時間の短縮を得た。

図 * から図  は、 格子数を  方向に  と固定し 、 る壁側と中心線上の 7   の変化を示した。

  の範囲で  数の変化によ表 # > 

  

*:*':**:**9*:'9#!

9!!!

9""

9%!



9#!

9!!

9"

9%"'

$

9#$

9!!$

9"%

9%!



9#!

9!!

9"

9%"%

 ? @  ()

9"%+++; ? 6   ()

+

9!"

9""

9%$

0+

9!"++A ? 8. ()+9!!#++@ (")++

9!*$+++

9""#+

B ? .

2 



()

? C  (*)

表 $ ,   074  D 1+%

E=F

8 074  '

=

%'

!'"'+

=#

$'

#'

+

=

!"%"#!''+=#

"!

!!%



=

##*'!* <+=#$'

*#'+

=

!'

!#$#%

+

=#

"''!

!""#"

$

=

!$'"

"%*

+

=#

"!$"

*!!'

+

=

"'!

!$"

+

=#

#%'$!

#*!



=

!

#$

+

=#

!%%%*

!$"'!

+

=

""'$

"#

+

=#

"%!#

#$!*

  

表 % ,   074  D 1+"074 

Pressure Drop

1
0.8



E=F

5+ 

6

 

#-=# # "# !!

-=#! " "# "-

=#   "$"-=# * "' %%=

 " #' "!

=

! " ! !'% ' <#=  "*# "#=

 # "** #

j=100
j=150
j=200
j=250

0.6
0.4
0.2
00

0.002

0.004

X

0.006

図 #   '

0.008

0.01

Pressure Drop

1
0.8

j=100
j=150
j=200
j=250

0.6
0.4
0.2
000.002

0.004

X

0.006

0.008

0.01

0.008

0.01

図 $   Pressure Drop

10.8

j=100j=150j=200
j=250

0.6
0.4
0.2
000.0020.004

X0.006

図 %   $Pressure Drop

10.8

j=100j=150
j=200
j=250

0.6
0.4
0.2
000.002

0.004

X

0.006

0.008

0.01

0.0015

0.002

Pressure Drop

図 '   

0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
00WallCenter Line

0.00050.001

X

図 *   'Pressure Drop

0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
00WallCenter Line

0.00050.001X0.0015

0.002

0.0015

0.002

Pressure Drop

図    

0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
00WallCenter Line

0.00050.001

X

図    $Pressure Drop

0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
00WallCenter Line

0.0005

0.001

0.002

0.0015X

図    むすび
流れ関数と渦度表示の方程式系に対し 、 高いベクトル化率での高速な計算を行うアルゴ リズムを実現した。 その結果、 計算時間はアルゴ リズムの改良で =#1+% にて実行 、 ハー
ド ウェアの変更 1+" から 1+% も含めると =*# にまで削減できた。 今後、 スーパーコ
ンピュータの発達により、 さらなる高速化が見込めると考えられる。 また、 計算結果より、 方向の格子変化は、 圧力分布の解析全体においては大きな影響をもたらすものではなかった。   での解析では、 平行平板の入口領域ではかなり大きな圧力勾配が 、 壁側と 中心線間にみられた。 なお、 この勾配は、  数の増加に従い小さくなっている。 中心線上の方が壁側よりも圧力が高いという結果となった。

謝辞
本研究は 、 東北大学情報シナジーセンターとの共同研究として行われた。 1+% 上でのプログラムの最適化及びベクトル化チューニングに関し 、 情報シナジーセンターより適切な助言を頂いた。 この場を借りて、 感謝したい。

参考文献() ;

F9 92. >

 6  C2 >G  C

9 9  C9 

3 7 7 9 $#+%$*$9  0

() 0

,90 .

7

9 3 2A9 9

9 #8

9 +"9(") @89 9@    /  8>C,/ F 3 C

CECC !$ C ($) B@90 

9 '*+*$,9

#/9/  2.   0 CECC 7



7



*'+%*

7



**'9

 2.9 7 9 /

9 *%+"- +C,/ 2.9 7 9 0

7 9 3 C,/ 2/  > 

***97 9 3 C,/ 299 !"+!!$

 @ 

7  7
"*$9 I , 3 5    0  C  3  8

9 '() 

-



--. @+7 /  >   7C,/ 2/>+9 #$92  >

2.  +***9>G   0    2.

 3*$#9/  > 

C,/ 2/>+9 #

(*) B 

 >  -- 8

9 "!+"* , ? /-9*%

3  ,  0   

7 

>

7 9 3 C,/ 2

9 '*+*$(%) B @9 0(')  @9E I , 3 5    0  >J 3  7 -
2/>+9 #3

9 E ,(#) B @9  6  B9  7 7*$9       (!) H 2.   E  3   ()  79 F99****'9() A

L9 89

C9E909/9 F

/G 3 +43  E ;7 J  3 F,/9

*%9,9  0 

9 "!+#

,92. -. 7  7

8;; 3 ,  E 5 -

8   /  

3 5.

 3 59

*%9            8.9 

79 C9
9!

8@   C-K K

2.   E  3 7  7

9 !!+!**$"9