第 27	  回 理論懇シンポジウム	
「理論天文学・宇宙物理学と境界領域」	

カノニカル乱流︎の	  大規模直接数値計算 	
石原卓	  (名古屋大学大学院工学研究科)	

2014年 12月25日(木)	  国立天文台 三鷹キャンパス すば︎る棟大セミナー室	自己紹介	  – 名古屋大学大学院工学研究科計算理工学専攻(1997-­‐)	  – 日本学振振興会 未来開拓学術研究推進事業「計算科学」	  • 地球規模流動現象解明のための計算科学(1997-­‐2002)	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  （代表：金田行雄 名古屋大学教授）	  – 乱流の直接数値計算(DNS)(格子点数5123,	  10243)	  VPP500,	  VPP5000＠名古屋大学	  – 地球シミュレータ (with	  横川、板倉、宇野)	  • ゴードン・ベル賞(SC2002)	  • 乱流の世界最大規模DNS(格子点数40963)	  Annu.	  Rev.	  Fluid	  Mech.	  (2009)	  – 「京」コンピュータ	  • カノニカル乱流の大規模直接数値シミュレーション	  	  	  	  	  	  	  	  	  	  	  	  (with	  横川、森下、宇野、金田)	  – 乱流の世界最大規模DNS	  (格子点数122883)	

共立出版(2012)	  「乱流」について	社会的関心の高い多くの問題が「乱流」現象	  ² 理論的／数学的研究が困難	  ² 実験による詳細な把握が困難	  ² 性質の理解が不十分	  「乱流」の特徴	  ² 決定論的な方程式、強い非線形性	  ² 巨大自由度(多スケール)	  計算機の有効な活用	  ² 巨大自由度が扱える	  ² 詳細な可視化と解析が可能	  ² 様々な「実験」が可能	  フォン・ノイマンの夢(1949)	  「乱流」現象は計算科学	  により解明・予測できる！	
3	

乱流の普遍的な統計法則	
• コルモゴロフの考え（K41）	  – 非線形性が十分に強い（レイノルズ数が大きい）
発達した乱流の十分に小さいスケールでは	  普遍的な統計法則が成立する	
非線形性の強さ	

E(k) = C εlog	  E(k)	  ε :平均エネルギー散逸率	  ν :	  動粘性率	  η = (ν 3 / ε )1/4 :	  コルモゴロフ長	

k −5/3

慣性小領域	  Re=UL/ν	:	  レイノルズ数	U3/L	  	  :平均エネルギー注入率	

2/3

エネルギー保有領域	  O(L)	

(5E(k) = ν ε

散逸領域	  1/4) f (kη)O(η)	Log	  k	  

乱流の自由度と計算コスト	
L	

非線形性の強さ	

Re:	  レイノルズ数	
Re=UL/ν	U	

時速100kmの車の後流	

Re=106-­‐107	

L

3/ 4~
Re
スケール比：	η3

η=(ν3/ε)1/4	
ν:動粘性係数，ε:平均エネルギー散逸率	⎛ L ⎞
9/ 4~~
Re自由度 	 ⎜ ⎟⎝ η ⎠35	計算コスト	 ~ Re

Recent	  development	  of	  supercomputers	

K	  Computer	
10PFlops	
Earth	  Simulator	
40TFlops		  		  	Re	

10years	

hTp://www.top500.org/staVsVcs/perfdevel/	

x1000	

10years	x10	

History of representative DNS
Homogeneous Isotropic Turbulence under periodic BC

Challenge for higher Re	
Yokokawa,	  et	  al(2002)	  Earth	  Simulator	  Jimenez et al.(1993)

10000	  Number of grid points

40963

Caltech Delta machine1283
Kerr(1985)

1000	  643

5123	  VPP	  Cray-1S NCAR

1Siggia(1981)
Cray-1 NCAR

100	  

243204835122883

Ishihara,	  et	  al(2014)	  K	  Computer	  10243
Ishihara＆Kaneda (2001),
Gotoh&Fukayama (2001)
VPP5000/56NUCC

Yamamoto et al(1994)NWT
10	  323

Vincent&Meneguzzi (1991)

Orszag(1969)

2403

IBM Model 360-95

1	  1960	  1970	  1980	  1990	  Year

2000	  2010	  7	

DNS	  of	  box	  turbulence	  using	  Earth	  Simulator	
40TFlops	  @	  640	  node	• Fourier	  Spectral	  method	  (Yokokawa	  et	  al	  2002,	  Kaneda	  et	  al	  2003)	

!∂$!
2# + ν k − c(k) & û(k) = P(k) ⋅ (u × ω )(k)" ∂t%
ω = ∇×uk ⋅ uˆ (k ) = 0
Forcing:	  negaVve	  viscosity	  (to	  keep	  the	  total	  energy	  constant	  )	⎧c (k < 2.5)c(k ) = ⎨⎩0 otherwise	  • 1D	  decomposiVon	  16.4TFlops	  @512	  node	  Eﬃciency	  is	  50%!!	

~40963	地球シミュレータを用いた乱流DNS	  (最大格子点数40963)による結果	

RV365-FL41-10

ARI

12 November 2008

14:55

The	  higher	  the	  Reynolds	  numbers,	  	  the	  wider	  the	  inerVal	  range	

The	  higher	  the	  Reynolds	  numbers,	  	  the	  stronger	  the	  intermiTency	

Small-scale statisti
3.5

a

3.0

0

k5/3E(k)/ ε 2/3∏(k)/ εRλ = 167
257
471
732
1131

2.5

10

Rλ = 167
257
471
732
1131

2 /( ε r)2/3
b10 (δu Lr)(a)10–2

σP (ζ)

2.0

1

1.5
1.0

Rλ = 167
257
471
732
1131

– (δuRLrλ)3=/(1131ε r)732
Rλ = 167
471
257
471
257
732
167
1131

ζ = ∂u/∂xHigher	  Rλ	10–410–60.50Figure 3

0.001

0.01

kη0.1

10.1

10–8
10–30
100–20

100–10r/η01000
10Rλ = 675

2030

2/3 (solidI,	  Gotoh,	  Kaneda	  ARF2009	
I,	  Kaneda,	  et	  aversusl.	  JFM2007	
(a) Compensated energy spectrum
k5/3 E(k)/⟨ε⟩colored lines) and energy-flux (c)spectrum "(k)/⟨ε⟩(dashedcolored lines)kη by429
268at
direct numerical simulation at various Rλ ’s. The solid gray line is the approximation (Equationwith Ko measured in the runζ =8)∂v/∂xEnergy	  Spectrum	  in	  the	  InerVal	  Range	  • K41	  	  E(k) = CK ε2/3

k−5/3∞, Π(k) ≡ ∫ k T (k#)dk# = ε• DNS	  of	  Turbulence	  in	  a	  Periodic	  Box	  with	  40963	  Grid	  Points	  E(k) ~ ε	  2/3

k −5/3−µ , µ ~ 0.10

Π(k) ~ ε• Similar	  slope	  – Hyper	  Viscosity	  SimulaVon	  Haugen	  &	  Brandenburg	  (2004)	  – DNS	  with	  40963	  grid	  points	  Donzis	  &	  Sreenivasan	  JFM	  (2010)	  InterpretaVon	  of	  the	  exponent	  	  	  	  	  …	  not	  yet	  Kaneda	  et	  al.	  PF	  (2003)	  乱流中の強い渦構造のReによる質的変化	
Re=3.6x104	  Rλ	  =O(103)	

Ishihara,	  Gotoh,	  Kaneda:	  ARF2009	

Ishihara,	  Kaneda,	  Hunt:	  FTC	  (2013)	
複雑な薄い剪断層	高Re	

L/λ=32	λ/η=	  66	  Re	  =	  uo	  L	  /ν	  Rλ	  =	  uo	  λ	  /ν	  	  	  L:	  積分長	  ~	  6	  L	

(外力のスケール)	  λ:	  テイラー長	  	  (慣性領域最小スケール)	  η:コルモゴロフ長	  Ishihara	  et	  al.	  JFM2007	

低Re	

0.7L	
大きさ	  ~	  O(L),	  厚さ	  ~	  O(λ)	 >>	  η		  (粘性のスケール)	

孤立した管状渦	

Re=3.2x102	  Rλ	  =O(102)	

長さ	  ~	  O(L)	  太さ	  ~	  O(10η)	
~	  O(λ)	

L/λ=3	λ/η=	  20	  ~	  6	  L	L	高Re乱流中の強い渦構造の統計解析	
Ishihara,	  Kaneda,	  Hunt:	  FTC	  (2013)	 One	  of	  acVve	  (high	  enstrophy)	  sub-­‐domains	

y	Based	  on	  one	  snapshot	  data	  of	  DNS	  	  of	  forced	  NS	  eq.	  with	  40963	  	  grid	  points	  	  (Kaneda	  et	  al,	  2003)	  Outside	
Inside	Outside	kmaxη=1,	  Rλ=1131	  y	512(=8x8x8)	  subdomains	

Enstrophy	  distribuVon	  on	  a	  cross	  secVon	x	z	z	x	
Analysis	  of	  CondiVonal	  staVsVcs	  z		  Bisset	  et	  al	  2002	  	  Westerweel	  et	  al	  2009	  	  da	  Silva	  &	  dos	  Reis	  2011	  xL	xR	x	CondiVonal	  averages	

Ishihara,	  Kaneda,	  Hunt:	  FTC	  (2013)	

ω 2 / 2 (ξ ) Ωω 2 / 2 (ξ ) Ω
4λ	Outside	−300η −150η

4λ	300ηw L (ξ ) w'
Outside	

Enstrophy	  in	  the	  layer	  	  ~	  10Ω	  	Outside	

Inside	Inside	ξ=x-­‐xL	 150ηλ	  ∝	  L	  Re	  -­‐1/2	

RL−300η −150η

ξ=x-­‐xR	
150η

w R (ξ ) w'4λ	

4λ	

Inside	Inside	Outside	

Existence	  of	  	  a	  layer-­‐scale	  shear	  	Δw=O(uo)	
−300η −150η ξ=x-­‐x 150η
L	300η

300η

−300η −150η

ξ=x-­‐xR	 150η

300η

Sharp	  interface	Sharp	  interface	  of	  internal	  thin	  shear	  layer	  (IKH2013)	 T/NT	  interface	  of	  external	  TBL	  (I,Ogasawara,Hunt2014)	

ω = 2 ω2

1/2

ω = 6 ω2

1/2

De-­‐correlaVon	  of	  the	  velocity	  ﬂuctuaVons	  across	  the	  sharp	  interface	  	  	  -­‐	  Blocking	  mechanism	  (Hunt	  &	  Durbin	  1999)	
CondiVonal	  cross	  correlaVons	  of	  u	Inside	

Outside	−300η −150η
ξ=x-­‐xR	 150η

300η

CondiVonal	  cross	  correlaVons	  of	  u	T	NT	Energy	  transfer	  T(x, k) and	  energy	  dissipaVon	  	  ε	  	  	  near	  the	  layer	

Ishihara,	  Kaneda,	  Hunt:	  FTC	  (2013)	

εLarge	  amplitude	  posiVve/negaVve	  (i.e.	  downscale/upscale)	  	  	  ﬂuctuaVons	  of	  T	  near	  the	  thin	  layer		  l :thickness	  of	  the	  layer	
<T(x,k)>Inside	  =	  <ε>Inside	  ~10	  <ε>	  	  for	  k	  >	  π/	  l10	  ~	  L/	  l ~ Rλ/100	
A	  net	  energy	  ﬂux	  from	  the	  larger	  scale	  moVons	  from	  outside	Data	  handling	  DNS	  with	  40963	  grid	  points	  	  High	  Enstrophy	  ranking		  	  (Kaneda	  et	  al,	  Phys.	  Fluids,	  2003)	  kmaxη=1,	  Rλ=1131(at	  a	  staVsVcally	  steady	  state)	  y	Divide	  a	  40963	  box	  into	  	  512(=8x8x8)	  boxes	  Ω=2174	

RANKING(new)	 nx	 ny	 nz	1	  2	  3	
4	
5	6	7	

Ω	  box	4	 5	 5	 6127 	
5	 6	 5	 5850 	
4	 6	 5	 5828 	5	 4	 5	 5639 	
6	 6	 5	 5360 	
5	 5	 5	 5091 	
4	 7	 5	 5036 	

x	509	
510	z	

Each	  box	  	  has	  	  3u,	  v,	  w,	  p		  512 	  grid	  points	  511	
512	

2	 7	 1	3	 8	 8	2	 7	 7	2	 7	 8	796 	
746 	
690 	623 	

Top10	  Enstrophy	  subdomains	①	②	③	

④	

4,5,5	

5,6,5	

4,6,5	

5,4,5	

⑤	

⑥	

⑦	

⑧	

6,6,5	

5,5,5	

4,7,5	

5,7,5	

⑨	

⑩	Clear	  shear	  layers	
Shear	  layers	  at	  	  subdomain	  boundaries	

4,8,5	

4,1,5	

No	  clear	  layer-­‐like	  structures	The	  most	  signiﬁcant	  shear	  layer	  (at	  t	  =	  t0)	  (In	  the	  most	  acVve	  subdomain	  of	  the	  40963	  DNS)	  ω = 6 ω2
ω =8 ω

0.7L	

1/2

2 1/2The	  most	  signiﬁcant	  shear	  layer	  	  from	  t	  =	  t0	  to	  t0+10τ (=t0+2.55λ/uo)	  0.7L	

The	  extreme	  points	  (	  	  ω	  	  	  >	   20 ω )	  	  in	  the	  most	  acVve	  subdomain	  2 1/2

Inside	  structure	  of	  the	  shear	  layers	
DistribuVon	  of	  the	  strong	  vorVces	  inside	  the	  layer	
Outside	Inside	Outside	Ishihara,	  Kaneda,	  Hunt:	  FTC	  (2013)	
(Δ=2π/4096~3η:	  grid	  spacing)	

ω2 2

w / uox	  Δw=O(uo)	
x/Δ	

x/Δ	

Thickness	  of	  the	  micro-­‐scale	  vorVces:	(insensiVve	  to	  their	  strength)	

∂w∂x

Very	  strong	  vorVcity	  of	  O(uo/10η)	  	  	>>	  uKol	  /η=1/τKol	  	  	  (K41)	
Velocity	  jump	  of	  O(uo)	  over	  distances	  of	  O(10η)	  	  	  	>>	  uKol	  ~	  uo	  Re	  -­‐1/4	  	  (K41)	
The	  layers	  may	  dominate	  the	  extreme	  point	  values	  of	  	  the	  staVsVcal	  distribuVons	  of	  dissipaVon,	  velocity	  and	  	  21	vorVcity	  ﬂuctuaVons	

High	  enstrophy	  regions	  in	  coarse-­‐grained	  data	  	  	Large	  clusters	  of	  the	  connected	  	  high	  enstrophy	  regions	Large	  clusters	  of	  the	  connected	  	  high	  enstrophy	  regions	

4,7,5	4,6,5	

5,7,5	

5,6,5	

6,6,5	

4,6,5	

4,5,5	

5,5,5	

6,5,5	

4,5,5	

4,4,5	

4,4,5	

5,4,5	

6,4,5	

高Re乱流中の薄い剪断層の特徴	
IKH	  FTC2013	

薄い剪断層	  強い管状渦が密に集合した渦の組織構造	  • 厚さ=O(λ∝Re-­‐1/2L),	  大きさ	  =O(L)	  • 強い剪断：O(λ)の層を横切って、O(uo)の速度差	  • 層の内部では非常に強い渦度と速度差(>>	  K41の予測)が生成	  • 剪断層は外部の速度ゆらぎをブロック・フィルタリング	  • 剪断層内部で強いエネルギー散逸とエネルギー輸送	  • 剪断層の周りでエネルギー輸送の強いゆらぎ（順カスケードと逆カスケード）	  DNS	  of	  turbulence	  in	  a	  periodic	  box	  on	  K	  computer	
11.28PFlops	  @	  	  88128	  node	• The	  same	  spectral	  method	  as	  used	  for	  Earth	  Simulator	  Forcing:	  negaVve	  viscosity	  (to	  keep	  the	  total	  energy	  constant	  )	⎧c (k < 2.5)c(k ) = ⎨⎩0 otherwise	  	  • A	  2D	  decomposiVon	  code	  for	  the	  K	  computer	  N	#	  of	  nodes	 TFlops	 eﬃciency	

6144	

96x64	

30.20	

3.84%	

8192	128x64	

32.93	

3.14%	

12288	 192x128	

70.46	2.24%	(double	  precision)	
Eﬃciency	  of	  50%	  was	  obtained	  on	  ES	  using	  512	  nodes.	

DNS	  with	  61443/81923/122883	  grid	  points	
kmaxη=1	3000	

Rλ- Vme	  history	

2500	

Rλ=1131	

40963	

61443	

81923	 122883	
Rλ=2314	Rλ=1482	

Rλ=1800	z	 5	
Re=1.5x10z	 4	
Re=8.9x10
Re=3.6x104	
(Ishihara	  et	  al	  JFM	  2007)	

Re=6.1x10z	 4	0	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  1	  	  	  	  	  0	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  		  	  	  1	  	  	  	  0	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  		  	  1	  	  	  	  	  	  	  T=L/uo~2	40963	

ω=6 ω

Rλ=1131	
Re=3.6x104	

~	  6	  L	61443	ω=6 ω

Rλ=1482	
Re=6.1x104	~	  6	  L	81923	

ω=6 ω

Rλ=1800	
Re=8.9x104	

~	  6	  L	61443	Rλ=1482	
Re=6.1x104	

81923	 Rλ=1800	
Re=8.9x104	

ω=6 ω

まとめ	• 地球シミュレータおよび「京」コンピュータを用いたカノニカル乱流の直接数値計算(DNS)では、各々最大格子点数40963、122883の高効率な計算を実現し、統計的準定常状態における最大レイノルズ数	  Re=3.6×104(Rλ=1131)、Re=1.5×105(Rλ=2316)	  の高Re乱流場のデータベースを構築した。	  
• エネルギースペクトル	  – Re>104,	  Rλ>600の乱流場はどれも E(k)~k-­‐5/3-­‐μ,	  (μ~0.1)	 となる波数領域	  – Re>>104ではその低波数側に真の慣性領域	  • 強い渦構造	  – Re<<104では孤立した管状渦	  – Re>104では渦の組織構造（薄い剪断層）	  剪断層：強い管状渦が密に集合した組織構造	  – 内部に非常に強い渦度、速度差が生成	  – Reが大きいほど微細な内部構造、シャープな境界	  課題	微惑星の形成過程	

・従来の問題点：乱流がダスト粒子の凝縮を阻害	  →微惑星が形成されない	  ・新たな理論(?)：高レイノルズ数乱流がダスト粒子の凝縮を促進	  高Re乱流DNSを用いた慣性粒子の運動のシミュレーション	
ストークス数	  St=tp/tη=0.01~1000	
慣性ありの粒子	

渦	慣性なしの粒子