኱つᶍ⛉Ꮫィ⟬ࢩࢫࢸ࣒㛵㐃᱌ෆ㸺኱つᶍ⛉Ꮫィ⟬ࢩࢫࢸ࣒㛵㐃ᴗົࡣࠊࢧ࢖ࣂ࣮ࢧ࢖࢚ࣥࢫࢭࣥࢱ࣮ᮏ㤋ෆࡢ᝟ሗ㒊᝟ሗᇶ┙ㄢࡀᢸᙜࡋ࡚࠸ࡲࡍࠋ㸼KWWSZZZVVFFWRKRNXDFMS㝵㟁ヰ␒ྕ䠄ෆ⥺䠅 
HPDLOಀ䞉ᐊྡ฼⏝┦ㄯᐊ

୺䛺䝃䞊䝡䝇ෆᐜ䝃䞊䝡䝇᫬ 㛫ᖹ᪥ィ⟬ᶵ฼⏝඲⯡䛻㛵䛩䜛┦ㄯ୍㝵  䡚
኱ุ䝥䝸䞁䝍䚸฼⏝⪅➃ᮎ➼䛾 䡚฼⏝ 

ྛ䝉䞁䝍䞊ᗈሗ䛾㜀ぴ䡚

Ṕ௦䛾኱ᆺィ⟬ᶵ➼䛾ᒎ♧䠖䡚䠖  
VRGDQ#FFWRKRNXDFMS┦ㄯဨ୙ᅾ᫬

฼⏝⪅ㄯヰᐊᒎ♧ᐊ ศᩓ䡶䢙䢇䢛䡩䡬䡼༤≀㤋 ྍど໬ᶵჾᐊ

 ୕ḟඖྍど໬䝅䝇䝔䝮䛾฼⏝䡚

⥲ົಀ  
VRP#FFWRKRNXDFMS⥲ົ䛻㛵䛩䜛䛣䛸䡚

఍ィಀ  
NDLNHL#FFWRKRNXDFMS఍ィ䛻㛵䛩䜛䛣䛸䚸㈇ᢸ㔠䛾ㄳồ䛻㛵䛩䜛䛣䛸䡚୕ ඹྠ◊✲ᨭ᥼ಀ   
UVVHF#FFWRKRNXDFMS㝵䠄䠅ඹྠ฼⏝ᨭ᥼ಀ
  ཷ㻌 ௜ XNHWXNH#FFWRKRNXDFMSඹྠ◊✲䚸ィ⟬ᶵ䝅䝇䝔䝮䛻㛵䛩䜛䛣䛸

䡚

฼⏝ᡭ⥆䛝䚸฼⏝┦ㄯ䚸ㅮ⩦఍䚸䝷䜲䝤䝷䝸䚸䜰䝥䝸䜿䞊䝅䝵䞁䛻㛵䛩䜛䛣䛸䡚

  
QHWVHF#FFWRKRNXDFMS䝛䝑䝖䝽䞊䜽䛻㛵䛩䜛䛣䛸䡚

ᅄ◊✲㛤Ⓨ㒊
㝵  ஬㝵 ➃ᮎᶵᐊ

 

3& ➃ᮎᶵ䠄; ➃ᮎ䠅䝛䝑䝖䝽䞊䜽ಀ㸦 㸧ෆࡣᮾ໭኱Ꮫෆࡢࡳࡢෆ⥺㟁ヰ␒ྕ࡛ࡍࠋ㟷ⴥᒣ࣭ᕝෆᆅ༊௨እ࠿ࡽࡣ㢌࡟  ࢆຍ࠼ࡲࡍࠋ
ᮏㄅ䛾ྡ๓䛂㻿㻱㻺㻭㻯䛃䛾⏤᮶㻌᫛࿴  ᖺ䛻ᮾ໭ᆅ༊䛾᭱ึ䛾㟁Ꮚィ⟬ᶵ䛸䛧䛶䚸ᮾ໭኱Ꮫ㟁Ẽ㏻ಙ◊✲ᡤ䛻䛚䛔䛶᏶ᡂ䛥䜜䛯䝟䝷䝯䝖䝻䞁ᘧィ⟬ᶵ䛾ྡ๓䛷 6(1$& 6(1GDL$XWRPDWLF&RPSXWHU 䛛䜙䛸䛳䛶࿨ྡ䛥䜜䛯䚹㻌

— 1 —

SENAC Vol. 49, No. 4（2016. 10）［共同研究成果］

多速さ格子ガス法実用化展開への手がかり―“テシャラの方法”から“マクスウェルの達人操作”へ―

松岡 浩
東北大学電気通信研究所(客員)
一般財団法人高度情報科学技術研究機構
菊池 範子
株式会社サイエンス・サービス

1.

はじめに

もう 20 年以上昔(1993 年)になるが、テシャラ(Chistopher M. Teixeira)は、彼の MIT における博士論文“Continuum Limit of Lattice Gas Fluid Dynamics[1]”において、CFD より格段にすぐれた多速さ格子ガス法を提案した。同じ格子点規模の流体シミュレーションを行うのに、CFDに比べて、1/10 のメモリサイズで 100 倍近い並列高速計算ができるという。筆者らがこれを知ったのは、10 年前のことであるが、以後、東北大学サイバーサイエンスセンターの SX シリーズを中心にしたベクトル型スーパーコンピュータを利用させて頂き、その実用化研究に取り組んできた。本稿では、テシャラの方法の概念と、筆者らが考える今後の方向性について紹介する。従って、学術論文の概要紹介とは異なり、一種の読み物として掲載させていただくものである。なお、
本稿の本文は松岡が執筆し、掲載したシミュレーション事例に関する計算及び可視化はすべて菊池が分担した。

２.

本稿の概要

格子ガス法は、流体が存在する空間中に格子を張り、多数の仮想粒子が、その格子点上のみで他の仮想粒子と衝突して進行の向きを変えながら格子点間を移動していく様子を平均化(疎視化)
して流体の挙動を模擬する方法である。このとき、仮想粒子がもつ質量・運動量・エネルギーが衝突の前後で保存されるような粒子衝突を想定する限り、その挙動は、自然界におけるある条件
下の流体挙動とかなり似かよったものになる。しかしながら、多少の違いが生じる。この主な原因は、仮想粒子がもちうる速度が、流体分子がもちうる速度とは異なり、制約があるからである。
すなわち、仮想粒子は、格子点が存在する位置の間しか移動できないので、速度の大きさも向きも離散的な値をとり、自然界の流体分子のように連続的な値をとることができない。この制約によって、格子ガス法が導くマクロな挙動は、連続流体を仮定している CFD が導く流体挙動とは多少異なったものになる。
“テシャラの方法”は、多速さ格子ガス法において、異なる速さをもつ仮
想粒子どうしの衝突頻度をうまく調節することによって、この“多少の差異”を解消する方法である。これにより、多速さ格子ガス法は、CFD を代替できる精度をもったシミュレーション手法となりうる。本稿では、多少の差異を解消するテシャラの方法の概念を身近な大きさのシミュレーション空間に対して考察してみる。このとき付記したシミュレーション事例は、概算によって数値の大きさを感覚的につかむことが目的である。
テシャラの方法は、個別の粒子を識別して、それに応じた操作を行い、自然界では生じないマ
クロな挙動を自己組織化により生成する。これは、負のエントロピーを生み出す“マクスウェルの悪魔”の仕事に似ている。ここでは、上記の“多少の差異”をピッタリとキャンセルする“達人技”のイメージを強調して、このような個別粒子に対する意図的な操作を“マクスウェルの達人操作”と名付ける。今後、この種のシミュレーション手法の実用化展開に当たって、自由な発想による“マクスウェルの達人操作”の導入がひとつの重要な手がかりになるであろう。

— 2 —

SENAC Vol. 49, No. 4（2016. 10）３.

格子ガス法による仮想世界と現実空間の対応

3.1

現実空間における４次元直交格子の配置

格子ガス法では、格子が張られた仮想世界の中を、仮想粒子が動きまわる。この仮想世界と現実世界の対応を実感できるように、あえて、具体的な数値と単位を使う。
まず、身近な大きさをもつ空間として、人ひとりがゆったり入れるぐらいの空間：“一辺(L)が2.56m の立方体（L×L×L）”（L＝2.56m）を考える。次に、この空間を“格子間隔(⊿ℓ)が 1cm の直交格子”（⊿ℓ＝1cm）で埋め尽くす。直交格子の x,y,z 軸は、立方体空間の各辺に平行になるよ
うに配置する。従って、３次元格子点は、x,y,z 方向にそれぞれ約 256 個並ぶので、３次元格子点の総数は全部で 256×256×256≒1678 万個となる。我々がシミュレーションで知りたい流体挙動の解像度(⊿L)は、もちろん目的によるが、ここでは、“手の指の長さ程度の解像度 8cm”（⊿L＝8cm）としておく。ただし、ここで注意が必要である。テシャラは、３次元シミュレーション用の直交格子として、対称性に優れた“４次元面心超立方体(FCHC: Face-Centered Hyper Cubic)格子”を用いている。従って、上記のすべての３次元格子点には、第４番目の自由度をもった座標を割りあてる必要がある。イメージ的には、各３次元格子点には、その格子点を中心にして半径 r が異なる多層の球殻が配置されており、仮想粒子は、半径 r が異なる球殻に別々に存在しうるものとする。こうして、４次元格子点の位置座標は（x,y,z,r）で与えられることになる。ここで、４次元目の座標位置 r も長さの単位で測ることができ、４次元格子点は、x,y,z 座標と同様に 1cm間隔で配置される。座標軸 r のとりうる範囲は、現実のシミュレーション対象である３次元物体の制約を受けないので任意性があるが、ここでは、r＝0,1,2,3cm の４点のみから成る“循環座標”としておく。すなわち、仮想粒子が、r＝3cm の位置にいて、さらに r のプラス向きに移動したい
ときは、r 座標は r＝0cm につながる。逆に、r＝0cm の位置にいて、さらに r のマイナス向きに移動したいときは、r 座標は r＝3cm につながると考える。

3.2

仮想粒子をめぐる物理量の設定

仮想粒子の速さについては、「静止粒子」、「遅い粒子」、「速い粒子」の３種類を仮定する。また、
��“仮想粒子の質量(m)はどれも同じで 1μg”(m ＝ 10 kg)とする。格子ガス法では、多数の仮想粒子が、同じ時間間隔(⊿τ秒)で時刻ステップがひとつ刻まれるたびに、格子点から格子点へと飛び移る。ここでは、この移動過程に要する時間⊿τを千分の 1 秒(⊿τ＝0.001sec)に設定する。
仮想粒子がもちうる速度については、注意が必要である。４次元 FCHC 格子を採用しているので、各格子点は４次元空間中に存在するため、格子点の間を飛び回る仮想粒子は、速度成分を４つも��� ＝(���� ,����� ,����� ,����� ) m/sec の形になる。ここで、
つことになる。すなわち、“４次元速度”�����
下付き添え字の「j」は、仮想粒子がもつ速度の大きさ(速さ)の違い、すなわち、エネルギーの違いを表す。j＝0 なら｢静止粒子｣、j=1 なら｢遅い粒子｣、j=2 なら｢速い粒子｣とする。また、「i」
は速度の向きの違いを表す。ただし、「静止粒子」の「i」は、単に複数の静止粒子の区別を表す。
仮想粒子の速さは、「格子点間距離 ℓ[1cm]を時刻ステップの刻み間隔⊿τ[0.001sec]で移動する速さｃ」のオーダーである。ここでは、この c＝ℓ/⊿τ＝1cm/0.001sec＝10 m/sec を基準にし��� を、次のように具体的に与えてしまう。エネルギーも求まる。
て、仮想粒子��がもちうる速度 �����

・「静止粒子（j=0,i=1～8）」がもちうる速度： 粒子がもつエネルギー(運動エネルギー)は、もちろんε� ＝0[J]。

������
��� ＝（0,0,0,0）m/sec、静止粒子は複数想定してもよい。通常、１個または６個の静止粒子を仮定することが多いが、ここでは、８個としておく。

・「遅い粒子（j=1,i=1～24）」がもちうる速度：

�
�
��
�
��
������
������
エネルギーは、ε� ＝ m(��� ･��� )＝mc ＝10 ･10 ＝10 [J]。
�— 3 —

多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —������
��� ＝（±c,±c,0,0）m/sec，（±c,0,±c,0）m/sec，（±c,0,0,±c）m/sec，（0,±c,±c,0）m/sec，（0,±c,0,±c）m/sec，（0,0,±c,±c）m/sec

・「速い粒子（j=2,i=1～24）」がもちうる速度：の 24 種類。
����������
������
エネルギーは、ε� ＝ �m(��� ･��� ) ＝ 2mc ＝2×10 [J]。

������
�
（0,±2c,0,0）m/sec，（0,0,±2c,0）m/sec，（0,0,0,±2c）m/sec，
�� ＝（±2c,0,0,0）m/sec，（±c,±c,±c,±c）m/secの 24 種類。

仮想粒子が、４次元格子点上のどこかに存在すれば、上記の速度をもって移動する限り、仮想粒子は設定された４次元 FCHC 格子から逸脱することはない。ただし、体系の周辺境界や固体境界等では、特別な移動規則を与える必要がある。また、仮想粒子は、格子点上でのみ、他の仮想粒子と衝突して移動の向きや速さを変えることができる。このとき、衝突の前後で、衝突に関与する仮想粒子合計の粒子数(質量)、運動量及びエネルギー(運動エネルギー)が保存される衝突規則を適用する。衝突には、｢遅い粒子｣どうし、｢速い粒子｣どうし、｢遅い粒子｣と｢速い粒子｣の衝突などのように、“仮想粒子数のエネルギー分布を変えない衝突”と、｢静止粒子｣と｢速い粒子｣が２個の｢遅い粒子｣と入れ替わる“仮想粒子数のエネルギー分布を変える衝突”の２種類がある。
最後に、重要な仮定として、「各４次元格子点においては、静止粒子を除き、複数の仮想粒子が同じ４次元速度をもつことはできない」という制約を課す。

3.3

マクロな挙動を観察する時空間スケール

仮想粒子の動きを平均化(以下、“疎視化”という。)してマクロな挙動を導出したい。まず、疎視化する範囲を決める。3.1 に前述したとおり、現実の３次元空間の中で求める解像度は、⊿L＝8cm とした。従って、間隔⊿ℓ＝1cm で４次元空間中に配置された格子点については、(x 方向 8cm)
×(y 方向 8cm)×(z 方向 8cm)×(r 方向 4cm)の４次元空間中に存在する 8×8×8×4＝2048 個の４次元格子点ごとに、そこに存在する仮想粒子の動きをまとめればよい。最後が「×4」になっているのは、第 4 番目の次元である r 方向には、４次元格子点が r=0,1,2,3 の位置に４つ存在するだ
けで、しかも、これらはすべて同じ３次元位置(x,y,z)に属しているからである。
また、時間についても、⊿τ＝0.001 秒という短い時間では、各仮想粒子は、⊿ℓ＝1cm のオー
ダーの距離しか移動できず、⊿L＝8cm の解像度が求められるマクロなオーダーでは大きな変化は生じない。特に、L＝2.56m の立方体内部における大きなスケールの挙動を追跡する場合は、⊿τ＝0.001 秒よりもかなり長い時間間隔⊿T 秒でスナップショットをとっていけばよいであろう。このスナップショット間隔⊿T の選択には、かなりの任意性があるが、ここでは、約 0.5 秒（⊿T＝0.512sec）とした。これは、仮想粒子の移動を 512 時刻ステップだけ刻める時間であり、仮想粒子が速さ c＝10m/sec で、立方体空間の一辺方向に他の仮想粒子と全く衝突せずに直進できたとしたら、二辺分の長さだけ進行できる時間である。例えば、初期状態から、⊿T＝0.512 秒ごとに 13
回のスナップショットを観察すれば、一辺が 2.56m の立方体空間における、約 6 秒間(6.144sec)のマクロな挙動を観察できることになる。

3.4

サブマクロな世界の記述

現実世界の流体を構成する流体粒子は、実際は、分子レベルの非常にミクロな世界を運動している。他方、仮想粒子が動き回る直交格子は⊿ℓ＝1cm の世界であり、現実の流体粒子は、そのよ
うな格子に動きを制約されることはない。この意味で、仮想粒子が動き回る世界は、現実にはありえない“仮想世界”である。しかし、この仮想世界における仮想粒子の動きを疎視化すれば、解像度⊿L＝8cm 程度のレベルで、現実のマクロな世界における流体挙動が得られることを期待しているわけである。ここで、仮想粒子が動き回る仮想世界は、スケール的に見ると、マクロな世界のすぐ下の“サブマクロな世界”であるということができる。

— 4 —

SENAC Vol. 49, No. 4（2016. 10）

サブマクロな世界の状態を完全に記述するには、各４次元格子点に、仮想粒子��が実際に存在するか否か？を示してやればよい。前述したように、各４次元格子点(x,y,x,r)には、①複数個の「静止粒子」（ここでは、これを仮に 8 個とした。）、②24 個のいろいろな向きの速度を持ちうる「遅い粒子」、③24 個のいろいろな向きの速度を持ちうる「速い粒子」が存在しうる。従って、ひとつの４次元格子点には、最大 8+24+24＝56 個の仮想粒子が存在しうる。ここでは、仮想粒子��が存在するか否か？を「1」または「0」の２値数(��� )で表現する。ただし、�� は、最大 56 個存在しうる仮想粒子の識別番号である。
��� (x, y, z, t): ４次元 FCHC 格子上の位置(x, y, z, t)にある４次元格子点において、時刻ステップ
ｔに仮想粒子��が存在する場合「1」、存在しない場合「0」の値をとる２値数。これにより、
格子ガス法において仮想粒子が動き回る４次元 FCHC 格子という“サブマクロな世界”の状態を完全に記述することができる。

3.5

マクロな状態量の導出

サブマクロな世界を記述する状態量 ��� から、マクロな世界を記述する状態量を導出する。はじめに、“質量密度ρ”と“運動量密度 P”を次式で ��� から定義する。
[kg/�� ]・“質量密度”： ρ(x, y, z, t)≡m･f ∑������ m･��� (x, y, z, t)
ここで、m：仮想粒子１個の質量[kg/個]

1μg＝10�� kg

f：４次元単位体積(１�� )あたりに存在する４次元格子点の数[個/�� ]。

４次元格子点は、４次元空間中に⊿ℓ＝1cm の間隔で並んでいるので、４次元の単位体積中には、f＝(1/⊿ℓ)･(1/⊿ℓ)･(1/⊿ℓ)･(1/⊿ℓ)＝10� 個/�������・“運動量密度”： ����������������������
P(x, y, z, t)≡m･f ∑���������� ･��� (x, y, z, t)}
�����
ここで、�
�� ：

[kg/㎥ sec]

仮想粒子�� がもつ４次元速度ベクトル[m/sec]

成分で書けば、 �� (x, y, z, t)＝m･f ∑����������� ･��� (x, y, z, t)}[kg/㎥ sec]�� (x, y, z, t)＝m･f ∑����������� ･��� (x, y, z, t)}[kg/㎥ sec]�� (x, y, z, t)＝m･f ∑����������� ･��� (x, y, z, t)} [kg/㎥ sec]�� (x, y, z, t)＝m･f ∑����������� ･��� (x, y, z, t)}[kg/㎥ sec]�����
ここで、���� , ���� , ���� , ���� は、各々��� のｘ成分,ｙ成分,ｚ成分,ｒ成分［m/sec］である。

次に、“マクロな速度”は、“質量密度”と“運動量密度”から次式で定義しておく。
����������������������
y, z, t)／ρ(x, y, z, t)
・“マクロな速度”： ����������������������
U(x, y, z, t)≡P(x,

��
�����
＝〔m･f ∑���������� ･��� (x, y, z, t)}〕／{m･f ∑���� ��� (x, y, z, t)}�������＝ 〔∑��
�������� ･��� (x, y, z, t)}〕／{∑���� ��� (x, y, z, t)}[m/sec]・・・mや f に依存しない。

��
成分で書けば、�� (x, y, z, t)＝{∑����������� ･��� (x, y, z, t)}}／{∑���� ��� (x, y, z, t)}[m/sec]���� (x, y, z, t)＝{∑����������� ･��� (x, y, z, t)}}／{∑���� ��� (x, y, z, t)}[m/sec]���� (x, y, z, t)＝{∑����������� ･��� (x, y, z, t)}}／{∑���� ��� (x, y, z, t)}[m/sec]���� (x, y, z, t)＝{∑����������� ･��� (x, y, z, t)}}／{∑���� ��� (x, y, z, t)}[m/sec]

多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —

3.6

— 5 —

マクロな速度の大きさの事例評価“仮想格子流体”が一様に流れている状況を考察する。このため、ｘ＝0 の面に存在する全ての４次元格子点から＋ｘの向きに、４次元速度が(ｃ,0,0,-ｃ)（ｃ,0,0,ｃ）（2ｃ,0,0,0）（ただし、ｃ＝10m/sec）となる３種類の仮想粒子を適当な頻度(確率 1/4)で送りこんでみよう。このときの初期条件としては、立方体内部のすべての４次元格子点において、1/4 の存在確率でランダ
ムに仮想粒子が存在しているものとする。このとき、質量密度は、次式の値となり、数値的には、常温常圧の空気の密度にほぼ等しい。もちろん、ｍの値の設定の仕方で、どんな物質のどんな状態の質量にも、値を対応させることができる。
��
kg×10� 個/�� ×(56×1/4)＝1.4[kg/�� ]ρ(x, y, z, t) ＝ m･f ∑������ �� (x, y, z, t)＝10

また、ｙ方向とｚ方向には周期的境界条件を課し、ｘ方向出口では、出口直前に位置する格子
点状態 ��� を出口直後に位置する格子点状態 ��� に上書きする方法で境界条件を与えた。

時刻ステップ＝0 の初期状態では、いろいろな向きの速度をもつ仮想粒子が、等方的に同じ確率(0.25)で各格子点に存在する。従って、最初は、マクロな速度はゼロである。しばらくすると、＋ｘ向きの運動量を得て＋ｘ向きの“流れ”が発生する。⊿τ＝0.001 秒の時刻刻みを約 1000 回(1024 回)繰り返した１秒後の時点では、多少のゆらぎはあるが、マクロな速度としては、全体が一様に＋ｘ向きに流れている。シミュレーション結果は省略するが、このときの仮想粒子の存否
データ：��� (x,y,z,r)からｘ向きのマクロな速度を算出してみると、約 0.5m/sec（数値データの意味は省略するが、110×10／2300≒0.48m/sec）であることがわかる。
一般に、４次元格子点において、仮想粒子がもつ速度の向きをそろえてある向きのマクロな速度をできるだけ大きくしようとしても、その最適な向きに存在できる仮想粒子は高々１個である。
また、最適な向きをあきらめて似たような向きの速度をもつ仮想粒子を増やしていくと、最適な向きとはますます異なった周辺向きの速度をもった仮想粒子の存在場所しか空いていなくなる。
従って、個々の仮想粒子の速度をある特定の向きに束ねるには限界がある。この結果、マクロな速度の大きさＵm/sec は、ｃ＝10m/sec よりもかなり小さい値になる。すなわち、Ｕ<<ｃという関係が多くの場合に成立している。

４.波型せん断流の減衰解析
4.1

はじめに

第３章で、格子ガス法による仮想世界を現実世界に対応させて、４次元速度が(ｃ,0,0,-ｃ)（ｃ,0,0,ｃ）（2ｃ,0,0,0）（ただし、ｃ＝10m/sec）となる３種類の仮想粒子を確率 1/4 の頻度で＋ｘ向きに送りこむと、１秒後には、一辺が 2.56m の立方体中における静止した“仮想流体”が、
約 0.5m/sec のマクロな速度を得ることを予想した。格子ガス法によるサブマクロな世界では、質量、運動量、エネルギーが保存されている。従って、それを疎視化して得られた“仮想流体”も、これらが保存されることで、実際の流体に似た挙動を示すであろう。しかし、このときの質量密度が、たまたま常温常圧の空気の密度に等しいからと言って、この“仮想流体”が常にマクロな世界における実際の空気の流れを模擬できていると即決はできない。ここでは、もう少し複雑な状況、すなわち、粘性が関与する状況を考察してみる。

4.2

波型せん断流の事例計算

ここでも、まずは、疎視化したときに「＋ｘ向きにマクロな速度 �� (x, y, z, t)＝約 0.5m/sec で一様に流れているマクロな状態」を再現できる“仮想格子流体”を考える。次に、ある時点で、瞬間的に適当な力を加えて、マクロな速度のｙ方向成分 �� (x, y, z, t) を波型にしてみる。これによってｙ方向に粘性によるせん断力が働くことになる。この瞬間の３次元画像を図１に示す。こ— 6 —

SENAC Vol. 49, No. 4（2016. 10）の状態から時刻ステップ⊿τで時刻刻みを進めていくと、マクロな速度のｙ成分 ܷ௬ の速度勾配が減少しながら、マクロな速度のｘ成分 ܷ௫ に乗って、＋ｘ向きに移動する“波型せん断流の減衰過程”がスナップショットで得られる。この減衰から、粘性の定量値等を求めることができる。

図１

瞬間的に生じさせた波型せん断流（ParaView による可視化）

シミュレーションによって得られた“波型せん断流の減衰過程”を図２に示す。図２には、比較のために、格子点間隔が⊿ℓ＝1cm の場合以外に、⊿ℓ＝2cm と粗くした場合、及び、⊿ℓ＝0.5cm、0.25cm と細かくした場合の結果を一緒に示してある。また、図２において、各ケース名の下に付記したνの数値は、動粘性係数の大ざっぱな概算値である。この算出方法は、4.5 に後述する。
一般に、格子点間隔を短くしていけば、仮想粒子の平均自由行程が短くなり、粘性は小さくなる。実際、図２から、格子点間隔が短い「ケース４」に近づくほど、波型せん断流の減衰の仕方
が緩慢になっており、粘性が小さくなっていることを見てとれる。
ケース１：
格子点間隔⊿ℓ＝2cm の場合
ν≒0.04 [㎡/秒]

格子点数：210 万
計算時間：
SX-ACE 32 ノード
約 14 秒
ケース２：
格子点間隔⊿ℓ＝1cm の場合

200

T=0

100
0

T=1
T=2
1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31

-100

200

T=0
T=1

100
0
-100
-200

T=4
T=5

-200

ν≒0.0264 [㎡/秒]

格子点数：1678 万
計算時間：
SX-ACE 32 ノード
約 130 秒

T=3

T=2
1 3 5 7 9 1113151719212325272931

T=3
T=4
T=5

多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —

ケース３：

200

格子点間隔⊿ℓ＝0.5cm の場合

T=0
T=1

100

ν≒0.0216 [㎡/秒]0

格子点数：1 億 3422 万

計算時間：
SX-ACE 32 ノード
約 27 分 30 秒

-100

T=2
1 3 5 7 9 1113151719212325272931

200

格子点間隔⊿ℓ＝0.25cm の場合
ν≒0.0056 [㎡/秒]
格子点数：10 億 7374 万

計算時間：
SX-ACE 32 ノード
約 8 時間 43 分
図２

-100

T=4

T=0
T=1

100
0

T=3
T=5

-200

ケース４：

— 7 —

T=2
1 3 5 7 9 1113151719212325272931

-200

T=3
T=4
T=5

波型せん断流の減衰過程（マクロな速度の y 方向成分 �� の変化）

（横軸は x 方向の位置[m]、縦軸はマクロな速度のｙ方向成分[m/sec]で、値は相対値である。）（T=0：静止状態、T=1：一様＋X 向き流れ、T=2：�� に対して波型速度分布を瞬間的に付与、T=2 以降：減衰。）

4.3“仮想格子流体”が導くマクロな挙動が満足する方程式

サブマクロな世界における４次元 FCHC 格子上の仮想粒子の衝突・移動という“仮想格子流体”のモデルは、それを疎視化してマクロなレベルで眺めると、現実の流体挙動に似た状態を再現できる。しかし、実際には、“多少の差異”が生じる。すなわち、実際の流体現象を正確に表現できると言われているナビエ・ストークス方程式の解との間には、ある差異が生じる。
テシャラの論文にも示されているが、ここで想定している“仮想格子流体”を疎視化した場合に、マクロなレベルで得られる方程式を導出してみると、次式の形になる。これは、ナビエ・ス
トークスとは少し異なる。
��＋�Ｕ･∇Ｕ ≒ νΔＵ��
これを、4.2 の事例計算で仮定した初期条件：
�� (�� �� �� �) = �� 、�� (�� �� �� �) = �� ��� ��、�� (�� �� �� �) = 0の下で解くと、次の過渡変化の理論解が得られる。
これが、今回想定した“仮想格子流体”の動きを疎視化した後に得られるマクロな挙動である。
�� (�� �� �� �) = ����� (�� �� �� �) = �� ���� � ･���(�� � ��� ��)

4.4

�� (�� �� �� �) = 0“仮想格子流体”に関するパラメータ�の評価“仮想格子流体”が満足する方程式に含まれるパラメータ�の値については、シミュレーショ
ン結果として得られた“波型せん断流の減衰過程”における位相の解析から、定量的な評価を行うことができる。

— 8 —

SENAC Vol. 49, No. 4（2016. 10）

上述の �� (�� �� �� �) の理論解によれば、スナップショット取得時刻Ｔにおける振幅ピーク位置
は、������ � ��� �Ｔ�＝1 を満たすｘで与えられる。実際、�� � ��� �� ＝��� から一番最初のピーク
位置を知ることができ、�������� ＝�� ��＋���� であることがわかる。同様にして、スナップショッ
ト取得時刻Ｔ＋⊿Ｔにおける最初のピーク位置は、�������(��⊿�) ＝�� �(� � ⊿�)＋���� である。そのピーク位置の差は、�������(��⊿�) － �������� ＝�� �⊿� となり、この値は k やＴには依存しない。ここで、�� はシミュレーション結果から測定(�� ＝0.48m/sec)でき、⊿Ｔは既知(⊿Ｔ＝0.512 秒)であるから、ピーク位置の差の測定値から � の値を求めることができる。
シミュレーション結果より �������(��⊿�) － �������� ≒0.16ｍであるから、次の値を得る。

�＝�������(��⊿�) － ��������
�� ⊿�

＝

�����

(���������)�(��������)≒ 0.65

さて、前述した“仮想格子流体”が満足する方程式：
��＋�Ｕ･∇Ｕ = νΔＵ��は、�=1 の場合、“一定密度・一定温度の場合の非圧縮性流体に関するナビエ・ストークス方程式”に等しくなる。すなわち、今回計算した“仮想格子流体”の挙動は、�≠1 であるので、実在の流体の挙動とは � の差に応じた差異がある。後述するように、テシャラは、静止粒子と速い粒子が衝突して２個の遅い粒子になる衝突の頻度を制御することによって、この差異を解消する方法を考案した。

4.5“仮想格子流体”に関するパラメータ�の評価“仮想格子流体”が満足する方程式に含まれるパラメータλについても、シミュレーション結果として得られた“波型せん断流の減衰過程”における振幅の解析から、定量的な評価を行うこ�とができる。前述の �� (�� �� �� �) の理論解によれば、波型せん断流の振幅ピーク値の減衰は、� ��� �に従う。ここで、νは“仮想格子流体の動粘性係数[㎡/秒]”、k は波数[1/m]、t は経過時間[秒]である。
“波数”とは、波型せん断流の初期条件である正弦波の単位長さあたりの位相進行量である。今回の波型せん断流の初期条件では、立方体を構成するＸ方向の一辺 L=2.56m の長さの中に２つの波を配置したので、k＝4π/L＝π/0.64[1/m]であることがわかる。
�
さて、あるスナップショット取得時刻Ｔにおける振幅ピーク値の理論解は �� � ��� � で、その次�のスナップショット取得時刻Ｔ＋⊿Ｔにおける振幅ピーク値は �� � ��� (��⊿�) である。従って、そ�の減衰比は、�� やＴに依存せず、� ��� ⊿� で与えられることがわかる。ここで、k と⊿Ｔは既知であるから、シミュレーション結果から測定できる振幅ピーク値の減衰比から仮想格子流体の動粘性係数νの値を求めることができる。
いま、格子点間隔⊿ℓ＝1cm の場合を例にとる。波型せん断流の減衰過程における最初のスナッ
プショット取得時刻Ｔでの最初の振幅の相対ピーク値は約 150 で、Ｔ＋⊿Ｔでの相対ピーク値は
約 100 になっている。従って、波型せん断流の振幅ピーク値の減衰比は、100/150≒0.67 である。
�
これより、0.67≒� ��� ⊿� からνの値を概算することができ、ν≒0.0264[㎡/秒]という値を得る。これは、空気の動粘性係数のレベルよりもかなり大きい。従って、テシャラの方法によって、
�=1 の補正がなされ、これが何かの条件下における実在流体の挙動を模擬できていたとしても、粘性値を見ると、それは空気ではなく、うすい“はちみつ”ぐらいに対応するものだとわかる。

4.6“仮想格子流体”を考えるアプローチからの視点

さて、岩波講座「現代工学の基礎」の“計算熱流体力学”に“物理現象のメゾスコピックモデル”に関して以下の記述[2]がある。

多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —— 9 —「流体のように、多数の分子とその相互作用でできている系は、分子間の相互作用と相互作用の間の相関がすぐに薄れてしまうという特性(分子カオス)をもつ。このため全体としてのマクロな挙動は、相互作用のミクロな詳細には依存しなくなる。相互作用の詳細は、マクロな係数には影響するが、方程式の代数構造には現れなくなる。このことを逆に考えれば、多数の分子相互作用の集積として現れる物理現象を、連続体近似でもなければ分子レベルでもない異なるレベルでモデル化することが可能であることがわかる。マクロなレベルに疎視化したときに同じマクロな挙動を示すように仮想的なレベルで相互作用をモデル化できれば、それをマクロな現象のモデルとして使えることになる。
」（cf.計算熱流体力学、笠木伸英、松本洋一郎、大橋弘忠、岩波講座、現
代工学の基礎、p108）
上記の記述の下線部を前述の“波型せん断流の解析”の場合にあてはめると、以下のことに対応していることがわかる。すなわち、“仮想格子流体”を疎視化して得られるマクロな挙動も、“実在流体”も、その支配方程式の代数構造は、
��＋�Ｕ･∇Ｕ = νΔＵ��という共通の形をしている。違いは、マクロな係数 � のみに現れたわけである。

５.マクロな係数の差異を解消するテシャラの方法
5.1

はじめに

マクロな係数の差異を解消する“テシャラの方法”の概念を述べる。本来は、サブマクロな世界における仮想粒子の運動から、マクロな挙動が導出される全過程を詳細に追跡する必要がある。
しかし、その導出過程はかなり長く複雑で、この部分を厳密に記載した参考論文([1]、[3]、[4]等)も当然存在するので、テシャラの方法が導入されるところまでだけの“大筋”を示す。具体的には、「局所平衡状態における仮想粒子分布を導出する式(後述の⑧式)」までを記述し、以後の過程は省略した。なお、“大筋”も、筆者なりの独自解釈で簡素化してあり、文献どおりではない。

5.2

サブマクロな世界における仮想粒子分布の記述

時刻 t に４次元空間内のある格子点位置(x,y,z,r)に到着した仮想粒子の“到着粒子分布”を����������������������������
�(���� �� �� �) 、衝突後の“出発粒子分布”を ������������������������������(�� �� �� �� �) と置く。これらの状態ベクトルの ji 成分

をそれぞれ ��� (�� �� �� �� �)、�� �� (�� �� �� �� �)と書き、仮想粒子 ji が存在するときは「1」、存在しないときは「0」と定義する。ここで、仮想粒子 ji とは、エネルギーが j で速度の向きが i の仮想粒子を意味する。ただし、静止粒子の場合(j=0)は速度の向きの区別はないので、i は、単に粒子の区別を表す。今、ji は 56 通りを考えるので、状態ベクトルは 56 ビットの“ビット列”で表現でき����������� (�＝1る。一般に、任意の 56 ビットのビット列について、全ての可能なビットパターンを、�(�)��～ 2 のいずれか)と書き、その ji 成分を ��� (�)と表現すれば、次式が成り立つ。
��� (�� �� �� �� �) ＝����〔��� (�) × ( �������������������������
�(�� �� �� �� �) �� �����������
�(�) )〕

���

・・・①式

����������������������������
������������������������������ �� �� �) �� ������������
�(�) ) は、２つのベクトル �(���� �� �� �) と ������������
�(�) が一致したと
ただし、上式で、( �(��きのみ｢1｣で、一致しないときは｢0｣になる２値数を表すものとする。

5.3“仮想粒子衝突規則”の表現

各格子点位置に到着した個々の仮想粒子は、そこに到着した他の仮想粒子と衝突を起こし、運動量やエネルギーを交換する。この結果、“到着粒子分布”は、別の“ビット列”で表現できる“出発粒子分布”に変換される。このときの具体的な変換内容を、“仮想粒子衝突規則”と呼び、���������� � �(�)����������) によって表す。ξ(�(�)���������� � �(�)����������) は、生じうるすべての入力ビット列のひとつひとつ
ξ(�(�)
������������
�(�) （�＝1～ 2 のいずれか）に対して、衝突後、それをただひとつの出力ビット列 �����������
�(�) （�＝1

— 10 —

SENAC Vol. 49, No. 4（2016. 10）���������� � ����������
～ 2�� のいずれか）に変換する規則を与える。ξ(�(�)�(�)) の値は、変換してよい入力ビット列���������� と出力ビット列 �(�)���������� が与えられたときのみ｢1｣になり、それ以外は｢0｣になる２値数である。
�(�)�(�)（�＝1～ 2�� ）に対して次式が成り立つ。
従って、任意の出力ビット列 �������������́ (�� �� �� �� �) �� ���������
�(�) � ＝ �� �����������������������

������������� � �����������ξ ��(�)�(�)� × � ����������������������������
�(�� �� �� �� �) �� �����������
�(�) ��

・・・②式“仮想粒子衝突規則”は、概念的には、2�� × 2�� の巨大な数値表になる。なお、すべての仮想粒子 ji が衝突を起こさず、素通りする場合も衝突の特別な場合である。この場合、ビット列に変化はないので、数値表の対角位置に｢1｣が並ぶ。また、数値表における｢1｣｢0｣の位置は、確率的に変動してもよい。すなわち、同じ入力ビット列に対して、確率的に異なるビット列を出力する“仮想粒子衝突規則”も可能である。いずれにしても、a,b＝1～ 2�� としたとき、任意の a に対して、
�

���

���

���������� � ����������
ξ(�(�)
�(�)) ＝ 1 ・・・③式 が成立する。

5.4

衝突過程における仮想粒子分布の変化

時刻 � での格子点位置(x,y,z,r)における“出発粒子分布”�����������������������������́ (�� �� �� �� �) の ji 成分を求める。
�́ �� (�� �� �� �� �)＝����〔��� (�) × ( ��������������������������́ (�� �� �� �� �) �� ����������
�(�) )〕

���

＝�

���

〔��� (�) × �

���

���

���

（∵ 前記①式の関係を適用）���������� � ����������
�(�)� × ( ����������������������������
�(�� �� �� �� �) �� �����������
�(�) )�〕�ξ ��(�)・・・④式 （∵ 前記②式を適用）

���������� � ����������
�(�)� × ( ����������������������������
�(�� �� �� �� �) �� �����������
�(�) )�〕�ξ ��(�)－④式から①式を引くことで、衝突過程における仮想粒子 ji の粒子分布の変化を求める。
�́ �� (�� �� �� �� �) － ��� (�� �� �� �� �)

＝�
＝�

���

〔��� (�) × �

���

���
���
���

�

���

���

－ �

＝�

���

���

�

���
���

���

���

�

���

〔��� (�) × ��

���

���

〔��� (�) × ( �������������������������
�(�� �� �� �� �) �� �����������
�(�) )〕

���

���

���

����������������������������
���������� � ������������ �� �� �) �� ������������
�(�) )〕 （∵ 値が 1 の③式を代入しても同じ。
）
ξ(�(�)
�(�))� × ( �(������������������������������
���������� � ����������
���� (�) × ξ ��(�)�(�)� × ( �(���� �� �� �) �� �����������
�(�) )� （∵ 上式の第２項で、�＝�，�＝�, �＝� とおく。
）
〔��� (�) × ��

���

���

�

����������������������������
���������� � ����������
���� (�) × ξ ��(�)�(�)� × ( �(���� �� �� �) �� �����������
�(�) )�

���

���

－ �

＝�

���

���

���

���������� � ����������
�(�) )〕
�(�))� × ( ����������������������������
�(�� �� �� �� �) �� �����������
ξ(�(�)

����������������������������
���������� � ����������
�(��� (�)－��� (�)) × ξ ��(�)�(�)� × ( �(���� �� �� �) �� �����������
�(�) )�

・・・⑤式

������������������������������ �� �� �) �� �����������
�(�) ) の具体的な評価式を考える。通常は、次式を用いる。
ここで、⑤式に含まれる ( �(�����(�� �� �� �� �) �� ���������
�(�) ) ＝ ∏��･�1 � ��� (�� �� �� �� �)�
( ��������������������������〔���� (�� �� �� �� �)�

� (�)����� (�)〕 ・・・⑥式

このとき、右辺に現れる ��� と ��� はいずれも｢1｣または｢0｣の値しかとらない２値数であり、「1� ＝1，1� ≡1，0� ＝0、0� ≡1」であることに注意すれば、右辺は、��� と ��� の全ビットが一致したときのみ｢1｣で、それ以外のときは｢0｣になることがわかる。このような“一致判定”は、２値数を扱う限りは、論理演算式で行う方が簡— 11 —

多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —単である。しかし、ここでは、後で集合平均をとるので、��� が確率として｢0 と 1 の間の実数値｣をとった場
合にもそのまま拡張できる評価式を用いたい。このため、⑥式の右辺の式表示を用いる。

5.5“仮想粒子衝突規則”の適用頻度を変えるテシャラの方法“テシャラの方法”を用いる場合、“仮想粒子衝突規則”の入力ビット列を作る前に、“到着粒子分布”に対してある特別な操作を行う。具体的には、“エネルギー ｊ の粒子が存在することを隠すマスク” �� と“エネルギー ｊ の粒子が存在しないことを隠すマスク” �� を導入する。�� と �� は２値数であり、通常は｢1｣で何の効果も発揮しない。マスク効果を発揮させるときのみ、値を「0」にする。シミュレーション計算では、ある一定の確率で「0」の値をとらせる。この場合、⑥式の表示は、次のとおりになる。
���(�� �� �� �� �) �� ���������
�(�) ) ＝∏��･��� (� � ��� (�� �� �� �� �))�( ����������������������
�������� ��� (�� �� �� �� �)�

� (�)����� (�)]・・・⑦式もし �� ＝0 であれば、静止粒子が存在(��� =1)していても、上式右辺で �� ��� =0 なので、入力ビット列ではその存在を検知できない。従って、「静止粒子」＋「速い粒子」⇒「遅い粒子」＋「遅い粒子」 という衝突変化は生じなくなる。また、�� ＝0 であれば、静止粒子が存在していない(��� =0)場合でも、上式右辺で�� (� � ��� )=0 なので、入力ビット列では存在していないことを検知できない。従って、上記と逆の「遅い粒子」＋「遅い粒子」⇒「静止粒子」＋「速い粒子」 という衝突変化は生じなくなる。最終的には、“「�� が 1 に���� 」の比”： �� ≡ ��� /����� が上記の衝突変化のバランスを決定することになる確率��� 」と「�� が 1 になる確率�
なる。すなわち、�� の値を調節して、局所平衡状態における仮想粒子数のエネルギー分布を制御できる。

5.6

局所平衡状態における仮想粒子分布を導出する式

��� (�� �� �� �� �) や �� �� (�� �� �� �� �) は、サブマクロな世界にある４次元 FCHC 格子上の時空間で、目���������� →まぐるしく｢1｣または｢0｣の値で変動している。また、5.3 に述べたとおり、“仮想粒子衝突規則” ξ(�(�)����������
�(�)) も、その数値表における｢1｣や｢0｣の位置が確率的に変化してもよい。ここで、このような確率的に変動する２値数を含む⑤式について、統計力学にいう“集合平均”をとる。詳しくは省略するが、局所平衡状態では「⑤式の集合平均＝０」の式が成り立つとして、これから、マクロな速度がＵであるときの仮想粒子分布を導出する。具体的に⑤式にテシャラの⑦式を代入して、集合平均をとると次式を得る。

0 ＝ �
＝�

���
���

���
���

�

���

���

���

�

�����������������������
����������������������������������
���������� → ����������
�(�)� × ( ����������������������������
�(�) )�
�(��� (�)－��� (�)) × ξ ��(�)�(�� �� �� �� �) �� ������������ (�)����� (�)����������������������
����������) ∏��������������������
���������� → �(�)� ������������������ �� ･���� ･(� � ����� (�) � ��� (�)}ξ(�(�)〕 ・・・⑧式
�� (�� �� �� �� �))��������� ･��� (�� �� �� �� �)�

���

ここで、集合平均をとるとき、これまで２値数として扱っていたいくつかの数は、0 以上 1 以下の実数値をとりうる数としての扱いに変更する必要が生じる。上式では、この数に上付きバーを付して区別した。ま���������� → ����������
����������������������������
�(�) )、�� と ��� 、 �� と (� � ��� ) の確率的変動は互いに独立なので、た、ξ(�(�)�(�)) と ( �(���� �� �� �) �� �����������
この積の集合平均をとると、その結果は、それぞれの集合平均の積になる。
⑧式から長い式変形過程を経て最終的に「マクロな速度がＵm/sec の局所平衡状態に対応する“仮想粒子数のエネルギー分布”」が以下のとおり求まる。ただし、u ≡(Ｕ/ｃ)<< 1 とする。
��������
��� (u)＝�＋����������
�����(�)��� �� ＋��(�)＋�� �(�)＋���� ･��� と定義するが、非熱流体解析用３速さモデルの場合は、�� = �� = 1 として、
ここで、�� ≡ ��� /��� のみを調節パラメータとして用いる。上式から一般に、�� の値を操作して、仮想粒子数のエネ
ルギー分布を変化させられることがわかる。マクロな質量、運動量、運動量フラックス等はすべてこの粒子分布から導出されるので、ナビエ・ストークス方程式等のマクロな支配方程式を正しく満足するように �� の値を調節して、マクロな係数の差異を解消することができる。

— 12 —

６.

SENAC Vol. 49, No. 4（2016. 10）“マクスウェルの達人操作”の展開

6.1 “マクスウェルの達人操作”
分子や原子の動きを個別に自由に操って、負のエントロピーが生じるような過程を生じさせる
想像上の小さな生物？で“マクスウェルの魔物”と言われるものがある。テシャラが ݃＝1 を達成するために用いた方法は、数多く存在する仮想粒子の衝突過程の中から、静止粒子が介在する衝突に注目し、静止粒子の“存在”や“非存在”をマスクによって隠すことによって、意図的にその発生頻度(確率)をある適切な値に制御するものであった。これによって、自然界では生じないようなある特別な“仮想粒子数の局所平衡エネルギー分布”をサブマクロな世界で実現し、それを疎視化したマクロな世界では、݃＝1 となる絶妙なバランスを得ることができる。ここでは、“多少の差異”をピッタリとキャンセルする“達人技”のイメージを強調して、このような個別粒子に対する意図的な操作を“マクスウェルの達人操作”と名付けることにする。
テシャラの方法は、“マクスウェルの達人操作”のひとつであると見なすことができる。本章では、この“マクスウェルの達人操作”のさらなる展開について考える。
6.2 負の粘性の発現
仮想粒子の衝突規則は、通常、すべての効果を合計したとき、“正の粘性”を発現させる[4]。
もし、“負の粘性”を発現させる特別な方法を見つければ、これを“正の粘性”を発現させる通常の衝突規則と確率的に組み合わせて、それを疎視化したときにマクロな世界で発現する“粘性”の値を自由に制御できるであろう。このことは、格子ガス法において、格子点数を増やさなければいけない２つの理由（①解像度を上げたい、②粘性を小さくしたい）のうち後者を解消し、格子点数は、もはや、求める解像度の要求からだけ決めることができるようになる。このことは、結果的に、必要とされる計算機パワーの大幅な節約を可能にする。
Daniel H. Rothman は、1989 年の論文[5]“Negative-Viscosity Lattice Gases”において、既に、“負の粘性”を発現させるひとつの方法を見出し、２次元の流体シミュレーションで粘性が実際に小さくなる効果を確認している。このとき用いられた方法は、仮想粒子に「赤」と「青」の色をつけてその違いに応じた処理を行う“着色法”の応用である。仮想粒子に対して通常の衝突をさせた後に、その格子点近傍の赤色粒子と青色粒子がもつ運動量の差の勾配が強調される向きに出発する仮想粒子に赤色を着色する。ここでは、“マクスウェルの達人”は、個々の粒子の動きをみて色を塗る着色屋さんの仕事をしている。
さて、粘性は、速度勾配があるところで、その効果を現す。例えば、ある短い時間、流線がある場所で局所的にある弧を描いたとする。いま、この流線上に位置する流体塊の重心速度がＵm/sec であり、弧の少し外側の流体塊の重心速度はＵより少し大きい速度Ｕ+εで、弧の少し内側の流体塊の重心速度はＵより少し小さい速度
Ｕ－εであったと仮定する。このとき、速度Ｕの向き(流線の接線
方向プラス向き)に横軸をとり、流線に垂直な方向に縦軸をとって、その近傍にある格子点間で、仮想粒子がどのようにやりとりされているか？を想像すると
図３のようになる。
図３

局所平衡的に速度勾配が維持された場所の仮想粒子分布(相対確率)

— 13 —

多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —図３では、議論を簡単にするため、２次元平面上の流れを想定し、仮想粒子の移動可能な向き
は、右上、右下、左上、左下の４つのみを仮定した。注目している流線Ｕ上の真ん中の格子点における、局所平衡状態時の仮想粒子の到着分布及び衝突後の出発分布が、図３中央の黄色点線枠
内に記載されている。すなわち、到着分布では、右上に向かう仮想粒子が左下の格子点から相対確率６の頻度で到着、右下に向かう仮想粒子が左上の格子点から相対確率７の頻度で到着、左上に向かう仮想粒子が右下の格子点から相対確率３の頻度で到着、左下に向かう仮想粒子が右上の格子点から相対確率３の頻度で到着していることを示す。また、これらの仮想粒子が衝突して移動の向きを変えた後が出力分布であり、仮想粒子は右上の格子点に向かって相対確率７の頻度で出発、右下の格子点に向かって相対確率６の頻度で出発、左上の格子点に向かって相対確率２の頻度で出発、左下の格子点に向かって相対確率３の頻度で出発していることを示している。
ここで、a. 速い側から到着する仮想粒子の寄与は、加速の寄与：減速の寄与＝７：２であるから、正身＋５の加速効果。
b. 遅い側から到着する仮想粒子の寄与は、加速の寄与：減速の寄与＝６：３であるから、正身＋３の加速効果。
c. 速い側へ出発する仮想粒子の寄与は、加速の寄与：減速の寄与＝２：７であるから、正身－５の減速効果。
（↑これは、速い側の流体に対しては、正身＋５の加速効果を与えている。
）
d. 遅い側へ出発する仮想粒子の寄与は、加速の寄与：減速の寄与＝３：６であるから、正身－３の減速効果。
（↑これは、遅い側の流体に対しては、正身＋３の加速効果を与えている。
）となり、流線上の運動量は収支がつり合って局所的な平衡状態を保持している。
以上の考察から、速度勾配が維持されている局所平衡状態では、仮想粒子の衝突前後における速度分布は、流線軸に対して反転されていることがわかる。
“負の粘性”が発現している状態とは、速度勾配がどんどん急激化していく変化途上の状態であるから、上記で考察したつり合い状態を破って、例えば、以下のような状態を作り出せば、“負の粘性”が発現するものと考えられる。
とりあえず、a.と b.は同じとして、c.と d.を次の e.と f.に変える。
e. 速い側へ出発する仮想粒子の寄与を、加速の寄与：減速の寄与＝１：８とし、正身－７の減速効果にする。
（↑これは、速い側の流体に対しては、正身＋７の加速効果を与える。）
f. 遅い側へ出発する仮想粒子の寄与を、加速の寄与：減速の寄与＝４：５とし、正身－１の減速効果にする。
（↑これは、遅い側の流体に対しては、正身＋１の加速効果を与える。）
流線を挟んで、速い側の流体に対する正身の加速効果は以前より増大し、遅い側の流体に対する正身の加速効果は以前より減少しているので、より大きな運動量勾配を作りだす方向に平衡状態が変化していくであろう。また、このときの新しい平衡状態に対応して、はじめに仮定した a.
と b.の状態も更新されていくことになる。以上の考察を図示すると、図４のようになる。
負粘性発現のために修正した

出発確率(相対値)

到着確率(相対値)

出発確率(相対値)

７

６
３

２

８
１
流線の向き

２

３
７
図４

６

４

５

負粘性を発現させるための仮想粒子の出発分布（相対確率）

— 14 —

SENAC Vol. 49, No. 4（2016. 10）

すなわち、“負粘性”を発現するには、通常の衝突規則によって得られた各３次元格子点における出発仮想粒子の速度分布(確率)を図４のように凹凸強調の修正をしてやればよいことがわかる。
図５は、筆者らが行った負粘性発現シミュレーションである。前述の“一辺が 2.56m の立方体”の３次元空間(ただし、全方向に周期的境界条件を適用)を想定した。
“負の粘性”の発現によって、静止状態から自発的に渦が発生して拡大していく様子がわかる。筆者らの方法では、「同じ３次元
格子点の異なる４次元格子点間で、それぞれの４次元格子点に存在する仮想粒子群の重心速度を測定し、その差(揺らぎ幅)を拡大する反発力(揺動力)を与える仮想粒子の速度ペアを見つけてそれらを反転させる」という操作を行ったが、これが“マクスウェルの達人”がした仕事である。

(a) 負粘性発現による渦の拡大過程（Z 軸に垂直な断面上）図５(b)ある時点の３次元可視化図
4 次元 FCHC 格子上で負粘性発現衝突規則を適用し疎視化して得たマクロな速度

6.3“マクスウェルの達人”のさらなる進化
Robert P. Bosch, Jr.は、彼の MIT における 1995 年の論文[6]“A Multi-grid Algorithm for
Lattice Gases”において、格子ガス法におけるマルチグリッド法を提案している。これは、CFD
における BCM(Building Cube Method)に似たものであり、物体に近くて解像度を上げたい領域では細かな立方格子を用い、物体から遠く離れた解像度をあまり必要としない領域では、粗い立方多速さ格子ガス法実用化展開への手がかり—“テシャラの方法”から“マクスウェルの達人操作”へ —— 15 —

格子を用いる方法である。これも、必要とされる計算機パワーを大幅に節約する方法である。両者の格子が接続される領域では、細かい格子の世界のみを動き回る軽い質量の仮想粒子８個と、粗い格子の世界のみを動き回る重い質量の仮想粒子 1 個が相互に変換される。しかし、同じ向きの速度をもつ軽い仮想粒子が８個集まるのを待っている間、質量や運動量が保存しなくなる。この過不足分の保存量を補うため、速度���で正の質量 m をもつ通常の仮想粒子に対して、速度は同�
じ���を有するが、負の質量:－m と負のエネルギー:－ � mv � をもつ“質量反転粒子”を導入している。
ここでの“マクスウェルの達人”の仕事は、２つの異なる格子点間隔の接続領域にやってくる仮
想粒子を観察して、必要に応じた数の“質量反転粒子”を生成または消滅させることである。
このように、格子ガス法における仮想粒子衝突規則には、それを疎視化したときに発現するマクロな挙動を支配する DNA が集約されており、衝突規則に対して“マクスウェルの達人”が行う操作を自由に発想することによって、まだまだ未知の可能性があることを感じさせる。さらに、
今や、“達人”を超える“人工知能（AI）”の時代が到来している。実のところ、仮想粒子の衝突規則に対するいろいろな操作が、それを疎視化したとき、どのようなマクロな挙動になって現れるのか？その因果関係を知ることは非常に難しい。しかし、他方、仮想粒子の衝突規則は、到着粒子の存在分布を表すビット列を、出発粒子の存在分布を表すビット列に対応させるものであるから、簡単にニューラルネットで入出力を扱うことができる。また、このとき、周辺格子点における仮想粒子の存在分布の影響も考慮する必要があれば、そこのビット列を入力に追加するだけでよい。いずれにしても、このとき用いるニューラルネットの基本構造は、流体中にあるすべての格子点において共通のはずである。また、格子ガス法の最大の特徴は、誤差が蓄積しないビッ
ト演算を使うことができるので、どんな激しい変化を示す流体挙動であってもその計算を安定に行い答えを出すことができるという点にある。もし、乱流について、深層学習を行ったらどうなるであろうか？人工知能を手に入れた“マクスウェルの達人”は、乱流状態における仮想粒子衝
突規則の操作方法を事前に知らなくても、乱流の測定データという自然界の教師データから何か
本質的なものを直接学びとるかもしれない。
これからの“AI コンピューティング時代”が楽しみである。

謝辞
本稿で述べた今後の研究展望は、これまで長期にわたり東北大学サイバーサイエンスセン
ターのスーパーコンピュータ(特に、ＳＸシリーズのベクトルコンピュータ)を利用すること
によってはじめて得られたものである。利用にあたっては同センター関係各位のご親切なご指導とご協力をいただき、心から深く感謝する次第である。
参考文献[1] Chistopher M.Teixeira, “Continuum Limit of Lattice Gas Fluid Dynamics”, MIT, 1993
[2] 笠木伸英, 松本洋一郎, 大橋弘忠, 岩波講座,現代工学の基礎, 計算熱流体力学, p108
[3] Uriel Frisch, Dominique d'Humiè res, Brosl Hasslacher, Pierre Lallemand, Yves Pomeau,
Jean-Pierre Rivet, “Lattice Gas Hydrodynamics in Two and Three Dimensions”, Complex
Systems,１(1987), pp.649-707, 1987
[4] Michel Hé non, “Viscosity of a Lattice Gas”, Complex Systems 1(1987), pp.763-789, 1987
[5] Daniel H. Rothman, “Negative-Viscosity Lattice Gases”, Statistical Physics, Vol.56,
Nos.3/4(1989), pp.517-524, 1989
[6] Robert P.Bosch,Jr.,“A Multi-grid Algorithm for Lattice Gases”,MIT,1995

— 16 —

SENAC Vol. 49, No. 4（2016. 10）［共同研究成果］

電流の偏りに注目したアンテナ放射効率の解析

東京農工大学

有馬 卓司，宇野 亨
大学院工学研究院 先端電気電子部門

身の回りは電波を利用した機器であふれている．無線電力伝送など新しい技術が現れており，
今後ますます電波の利用が多くなると考えられる．これら機器においては，電気信号もしくは電力の送受信はアンテナを用いて行われる．アンテナ自身に効率があり放射効率と呼ばれる．放射効率の良いアンテナでは電波をほぼ損失無く送れるのに対し，放射効率の悪いアンテナでは電波の多くを損失してしまう．このようにアンテナの放射効率は非常に重要なパラメータである．アンテナの放射効率は材質にも左右されるが，良導体を用いていても条件によっては効率が低下することが知られている．放射効率が低下するメカニズムを解明できれば，その改善に寄与することができ，省エネルギー化や携帯情報端末の電池のもちを改善することができる．本研究では，アンテナに導体が近づくことで，電流分布に偏りが生じることに注目し，偏りと放射効率の関係について解析的に評価している．解析には数値電磁解析手法の一つである FDTD 法を用いた．

1.

はじめに

現在身の回りの多くの電子機器にワイヤレス機能が実装されている．これらは無線 LAN を利用してインターネットに接続され，非常に多くの情報を利用することができ，生活が一段と豊かになった感がある．今後も IoT（Internet of Things）の普及によりますますワイヤレス機器が増えてくると考えられる．これら機器においては，情報の送受信はアンテナを用いて行われる．このためアンテナは電波と情報のインターフェイスとして用いられるため，アンテナの性能により通信の質が左右されると言っても過言ではない．アンテナ自身に効率があり放射効率と呼ばれる．
放射効率は，入力電力と放射した電力の比で定義される．そのため放射効率の良いアンテナでは電波をほぼ損失無く送れるのに対し，放射効率の悪いアンテナでは電波の多くを損失してしまい放射させることができない．以上よりアンテナの放射効率は非常に重要なパラメータである．
筆者らのグループではこれまでアンテナの放射効率低下メカニズム解明に関する研究を継続して行ってきた．たとえば，文献[1,2]では良導体のみで構成したアンテナにおいてもその条件により，アンテナ放射効率の低下が起きる事を報告してきた．また文献[3]では，アンテナを構成する
導体板の表面の粗さがその導電率に及ぼす影響について報告してきた．この研究では，導体自体の導電率は高いのであるが，その表面の粗さによって導電率が低下することを指摘している．この原因は，アンテナの表面上に電流が集中するので，その電流が表面粗さの影響を大きく受けることと考えられる．本研究では，アンテナの近くに導体が接近すると，アンテナ上の電流分布に偏りが生じる事による放射効率低下のメカニズムの解明を行った．これは，偏りにより電流分布がより表面に分布する事が原因と考えられる．解析には，電磁界解析の代表的な手法である FDTD
法（Finite Difference Time Domain method）を用いた[4,5]．この手法は，時間領域のマクスウ
エルの方程式を直接差分し時間領域で解析する手法である．FDTD 法は比較的容易に複雑なモデルを解析できることが強みであるが，解析モデルによっては比較的多くの計算時間が必要となる．
本研究の解析モデルにおいてもアンテナ内の電流分布を精密に計算するため多くの計算時間が必要になる．そこで， FDTD 法の解析時間の高速化について東北大学サイバーサイエンスセンターの多大な協力を得た．

— 17 —

電流の偏りに注目したアンテナ放射効率の解析

２.

解析モデルと解析手法

仮想閉曲面 S

解析領域

J
Ms，Js

図１．解析モデル

図２．遠方界解散手法

本研究では，図 1 に示す解析モデルを用いた．このモデルでは，一辺の幅が 10mm で全長 150mmのアンテナの近傍に同じ長さの導体棒が配置されたモデルである．この導体棒とアンテナの間隔
d を変化させることでアンテナ上の電流分布に偏りを生じさせその際の放射効率の変化を確認した．
FDTD 法は解析モデルを微小領域に分けて解析するが，本解析ではその大きさを 1mm とした．このため，アンテナは各幅方向に 10 分割されているのでアンテナ上の電流分布を詳細に解析出来る．
アンテナの放射効率は次式で定義される．

η�

����

(1）

���

η :放射効率��� :入力電力���� :放射電力

FDTD 法では入力電力は給電部の電圧と電流を１周期積分する事により求まる．

Pin 

1 T
V ( t )I( t )dt (2)T 0 0放射電力はアンテナから放射される電波のうち遠方に届く遠方界と呼ばれる成分が寄与するので，遠方界を求める必要がある．アンテナの遠方界は図 2 に示すようにアンテナを仮想的な閉曲面で囲い，閉曲面上の等価電磁流を解析する．解析した電磁流を数値積分することにより遠方界での電界成分 E,Eを求めることが出来る．求まった，E,E より，指向性関数 D を用いて最終的に

D( , )  E 2( , )  E 2( , )
2D( , )1
Prad    
sin d d dt
T 0 0 0Z0
T 2 

(3)(4)(4)式より，放射電力を求めることが出来る．以上より求めた入力電力と放射電力より(1)式の計算を行い効率を求める．

— 18 —

３.

SENAC Vol. 49, No. 4（2016. 10）

電流の偏りと放射効率の関係

current density[A/m

2]

0.6

0.4

0.2

0
−1

図 3．座標系1

図 4．アンテナ単体の電流分布3current density[A/m

2]

2]

2

current density[A/m0y/(ny/2)

1

0

2

1

0
−1

0
y/(ny/2)1

図 5．電流分布（d=15mm）−10y/(ny/2)1

図 6．電流分布（d=10mm）

図 3 に座標系を示す．図 4 にアンテナ単体の電流分布を示す．図 4 は横軸がアンテナの中心軸上の電流分布であり，半分の幅が 1 となるように正規化している．解析した周波数は 1GHz の連続波を給電部により給電している．当然ながらアンテナ単体では電流が対照的に流れており偏りは生じていないことが分かる．次にアンテナ近傍に導体棒を配置した際の電流分布を解析した結果を示す．図 5 は，アンテナと導体棒の間隔 d を 15mm とした際の電流分布，図 6 は，アンテナと導体棒の間隔 d を 10mm とした際の電流分布である．横軸は，図 4 と同じく中心軸上の電流分布であるが，アンテナ単体の際と比べて電流分布の大きさが，左に多く偏っていることが分かる．また
その偏りは，アンテナと導体棒の間隔 d が小さくなるとより顕著に偏りが生じている事が分かる．
次に，この際の入力電力を精度よく計算する手法について述べる．図 7 から図 9 に入力電力 Pinの時間変化を示す．このように，導体棒が近傍にあると入力電力がなかなか収束しないことが分かる．これではいつまでも計算を終わらせることができない．そこで本研究では収束値を予測する手法として ARMA と呼ばれるアルゴリズムを導入した[6]．ARMA アルゴリズムでは計算値を次式で示す有理関数で近似する手法である．
H�z� �

���� �� � �� � �� � �� � �� � � �� � ��
�(5)
� � �� � �� � �� � �� � � �� � ������

この際 H は Y と X の関数になっているが，本研究では H を放射効率としているので，Y は放射
電力，X は入射電力である．z は z 変換のパラメータである．

— 19 —

電流の偏りに注目したアンテナ放射効率の解析

図 7．電流分布（単体）

図 8．電流分布（d=15mm）

図 9．電流分布（d=15mm）

図 10．放射効率

上述した放射電力と入射電力の計算手法を駆使して解析したアンテナの放射効率を図 10 に示す．この解析では，ARMA アルゴリズムを組み込んでいる．図 10 の横軸は，アンテナと導体の距離であり，縦軸はアンテナの放射効率である．放射効率はアンテナと導体の距離に依存することが分かる．特に，本研究の例では，15mm 付近で効率が大きく落ち込んでいることが分かる．本研
究のアンテナの長さ等は一例であるが，この例ではアンテナの 15mm 近くに導体棒があると顕著な放射効率の低下が起きることが示唆された．この条件を回避することでアンテナの効率低下を防ぐことができると考えられる．

５.

まとめ

本研究では，アンテナ近傍に導体棒が有る際に起きる電流の偏りとアンテナの放射効率の関係について解析を行いその特性を明らかにした．解析には FDTD 法を用いた．解析の結果アンテナの電流が偏ると放射効率が低下することが示唆された．また，ある条件において顕著に放射効率の低下が起きることを示した．本研究で用いたアンテナは一辺の幅が 10mm で全長 150mm のアンテナであったが，この例ではアンテナの 15mm 近く導体棒があると顕著な放射効率の低下が起きることが示唆された．この条件を回避することでアンテナの効率低下を防ぐことができると考えられる．
今後さらなる検討により，アンテナの放射効率低下が起きる条件やその詳細なメカニズムを明らかにすることによりアンテナの放射効率が低下しない条件を明らかにする事が出来ると考える．
そのためにはさらなる解析が必要である．

謝辞
本研究を遂行するにあたり，東北大学サイバーサイエンスセンターのスーパーコンピュータを利用した．この恩恵にあずかり本研究が遂行できた．また，研究にあたっては同センター関係各位に有益なご指導とご協力をいただき遂行できた．

— 20 —

SENAC Vol. 49, No. 4（2016. 10）

参考文献[1] 山本雄太，宇野亨有，馬卓司，何一偉 "FDTD 法におけるアンテナ放射効率計算法に関する検討", 電子情報通信学会技術研究報告，AP2011-71，pp.37-42[2] 山本雄太，宇野亨，有馬卓司 "FDTD 法における遠方界を用いたアンテナ放射効率計算法に関する研究"，2012 電子情報通信学会全国大会論文集，B-1-175, p.175[3] 北住亮祐，宇野亨，有馬卓司，"表面粗さによるアンテナ放射効率低下に関する検討, "電子情報通信学会技術報告, vol.114, no.94, AP2014-130, pp. 13-16, 2014 年 11 月
[4] 宇野亨，FDTD 法による電磁界およびアンテナ解析，コロナ社，1998.
[5] 宇野亨，何一偉，有馬卓司，数値電磁界解析のための FDTD 法，コロナ社，2016.
[6] Takuji Arima, Akie Kuriyama and Toru Uno, "ARMA/FDTD Technique for Reflection
Properties of Periodic Structures at Oblique Incidence and Developing Wide Angle EBG
Application, " Proc. 2015 IEEE 4th Asia-Pacific Conference on Antennas and Propagation,
Bali Island, Indonesia, pp.8-10, 2015.7.1

— 21 —

SENAC Vol. 49, No. 4（2016. 10）［研究成果］— 情報科学研究科グループ利用制度を利用して —Shape Optimization による非定常流れ場の安定性制御
Shape Optimization Problem Improving Hydrodynamics Stabilityon the Non Stationary Navier-Stokes Problem
中澤嵩
東北大学大学院情報科学研究科数学連携推進室
本研究では，運動エネルギー最小化問題を行い，初期形状と最適形状における線形安定性を比較する．その際，初期形状を 2 次元 Cavity 流れとし，レイノルズ数Re = 11500において非定常
Navier-Stokes 問題を主問題とする．数値計算の結果，初期形状と比較して最適形状における臨界レイノルズ数Re� が減少したことが数値的に示された．運動エネルギーを最小化しているにも拘らずRe� が減少するというこの結果は，流体物理の観点では非常に興味深い．今後，更なる詳細な検証が求められる．

1. はじめに形状最適化問題は，偏微分方程式が定義された空間における境界の形状を最適化することである．初期には，1903 年に J. Hadamard によって薄膜の振動数最大化問題として行われた．その後，特にフランスの応用数学者によって線形弾性体が定義された空間における，散逸エネルギー最小化問題が数学的に解析され，形状最適化に必要な framework が数学的に構築された．
有限要素法により支配方程式の弱形式を導き，Lagrange 未定乗数法により Lagrange 関数を定義する．その際，Lagrange 乗数は有限要素法における trial function と結果的に同じものとなる．そして，この Lagrange 関数に対する物質微分をとることで，密度関数の形状微分が得られる．ここでの形状微分は，固定された座標において領域変動に対する密度関数の変動を意味しているため，Fréchet 微分そのものになる．そして，随伴変数法により，Lagrange 関数の物質微分を評価することで目的関数を最小化するための感度（または領域摂動）が得られる．このような解析を一般に感度解析と呼ぶ．
1973 年には，流体力学における要請から，O. Pironneau によって Stokes 問題[2]が，1974 年に定常 Navier-Stokes 問題[3]が定義された空間における散逸エネルギー最小化問題に対する感度解析が行われた．
近年になり，計算機の発達に伴い数値的に流体機器の最適設計を行うことが求められ，形状最適化問題に対する数値計算手法の構築が行われた[4]．ここでは，数値的に評価された感度（または領域摂動）を逐次的に初期の領域に足していくことで，目的関数を最小化するというものであった．しかし，形状更新の度に境界の滑らかさが落ちていく“波打ち現象”が観察されたため，境界の平滑化が必要となった．その際，大きく分けて Parametric 的または Non-Parametric 的な手法が提案された．特に，Non-Parametric 的手法として代表的な手法は，H. Azegami によって提案された H1 勾配法である[5-7]．この H1 勾配法では，領域を線形弾性体と仮定して，感度を境
界で Dirichlet 条件として定義し，領域摂動を平滑化させることを可能にし，誤差解析も行われている．そして，定常 Navier-Stokes 問題に対して，この H1 勾配法を用いた研究が E. Katamine et— 22 —

SENAC Vol. 49, No. 4（2016. 10）

al. [8, 9]及び T. Nakazawa et al. [10-12]によって行われた．更に，非定常 Navier-Stokes 問題に対しては Y. Iwata et al. [13]によって行われた．
流体工学において流体機器の最適設計は実用上，非常に重要な課題であるため，形状最適化問
題が広範囲に活用されてきた分野でもある．特に，航空工学では，2 次元翼断面[14]や 3 次元翼[15]の形状最適化が行われている．その際，揚力と抗力の比である揚抗比が目的関数として用いられている．しかしながら，航空機・高速列車は高速で飛行・走行するため機体の表面において壁面乱流が発生し，この乱流によって飛行・走行安定性を著しく悪化させたり，乱流そのものによる抵抗が予想される．そこで，乱流発生を抑制するために，流れ場の安定性を改善するような形状最適化問題を構築する必要がある．しかし，国内外問わず，そのような研究は著者が知る限り，
ほとんど行われていない．
本研究では，流れ場の安定性を改善するような形状最適化問題を構築することを目的として，散逸エネルギーを目的関数と定義し，非定常 Navier-Stokes 問題を主問題として用いた形状最適化問題を解き，初期形状と最適形状における線形安定性を比較する．その際，初期形状には 2 次
元 Cavity 流れを採用する．本論文は以下のような構成となっている．セクション 2 では，初期形状を定義し，領域の変動を数学的に定義する．そして，主問題である非定常 Navier-Stokes 問題を導入し，形状最適化問題を構築する．そして，線形安定性解析の概要を説明する．セクション
3 では，数値計算結果を示し，セクション 4 では本研究の考察を行う．

2. 問題の設定
本セクションでは，運動エネルギーを目的関数と定義し，非定常 Navier-Stokes 問題を主問題とした際の形状最適化問題と，流れ場の線形安定性解析について述べる．
セクション 2.1 では，初期形状Ωと領域変動�を定義し，汎関数�とその密度関数�に対する物質微分をとることで，汎関数�の形状微分公式を導出する．
セクション 2.2 では，主問題である非定常 Navier-Srokes 問題の弱形式を導き，その際，特性曲線法を導入する．
セクション 2.3 では，Lagarange 未定乗数法を用いるかたちで，運動エネルギーを目的関数，
非定常 Navier-Srokes 問題を主問題とする，Lagarange 関数を定義する．その際，セクション 2.2で弱形式を導出する際に用いた trial function は Lagarange 乗数法に相当する．そして，�の形状

微分公式により，Lagarange 関数の目的関数を最小化させるような領域摂動に対する変分をとり，感動を評価する．
セクション 2.4 では，本研究に対する線形安定性解析について定式化する．

2.1 初期形状と領域変動
ユークリッド空間における，2 次元領域Ω� ∈ ℝ� を非有界な凸領域とする．そして，Ω ⊂ Ω� 初期

形状とし，位置ベクトルを� ∈ ℝ� と記述する．その際，

— 23 —

Shape Optimization による非定常流れ場の安定性制御

Ω = {� = (�, �) ∈ ℝ� | � � � � 1, � � � � 1},Γ��� = {� = (�, �) ∈ ℝ� | � � � � 1, � = 1},�����Γ���� = �Ω � Γ��� .

次に，ℝ� 値の関数である Lipschitz transformation �� Ω � �(Ω)による領域変形を考える．この写像は，微小パラメーター|�| ≪ 1 を用いて，以下のように記述できると仮定する．
� = �� � �� � �(� � )ここで，�� は恒等写像であり�� (�) = �となる．また，�は領域変動である．そして，�の関数空間は M. Kimura[16]の Proposition 1.41 と eq. (1.18)より，

�

�(Ω) = �� ∈ � �,� (Ω; ℝ� ); �(��� )� � ����� � � �

� (�)������� ⊂ Ω� �.
� 1, �(Ω)

これから，汎関数の形状微分公式を導く．まず， � = {�� , �� , � } ∈ �� (Ω)をΩにおける物理量を表す実数値関数とする．そして，�(�, Ω)を密度関数として以下のような汎関数を定義する．
�(�� , �, Ω) = � �(�, Ω) ���また，変形領域における汎関数は以下のように記述することが可能である．

�(�, �(�), �(Ω)) = ��(�(�), �(Ω)) �(��)�(�)この変動領域�(Ω)における�(�, �(�), �(Ω))を，初期領域Ωにおいて記述する．
�の Jacobi 行列は，
(��� )�

��= ����� � � ����� � � �(� � )=����(�� , �� )
� �(� � )�(�� , �� )と記述することが可能であり，Jacobi 行列はdet((��� )� ) = 1 � �� � � � �(� � )となる．よって，以下の結果が得られる．

�(��)= det((��� )� )�� = �1 � �� � � � �(� � )���— 24 —

SENAC Vol. 49, No. 4（2016. 10）

次に，密度関数���(�)�に対して物質微分をとると，
�(�(�), �(Ω)) = �|���� + ���� + (� � �)��|���� + �(� � ).

ここで，( )�は物質微分を(�)は形状微分を表す．更に，形状微分は固定された座標系における領域変動であるので，Fréchet 微分と同意である．

���(�), �(Ω)�と�(��)を�(�, �(�), �(Ω))に代入することで，

���, �(�), �(Ω)� = � �� + ���� + (� � �)�� + �(� � )�|���� �� + �� � � + �(� � )����

= � �� + ���� + (� � �)� + �� � �� + �(� � )�|���� �� .
�となり，���, �(�), �(Ω)� � �(�� , �, Ω)
���
�
lim

�)
= ��(�, �, Ω= � ��� + (� � �)� + �� � ��|���� ���= � ��� + � � (��)�|���� ���= � �� (�, Ω) �� + � �(�, �Ω)� � ��� .
���

ここで， �は境界上の外向き単位法線ベクトルである．また，本研究ではΓ���� が設計境界となっているので，�の関数空間を以下のように再定義する．

�(Ω) = {� � �(Ω)� � = � �� Γ���� }.
そのため，
��(�, �, Ω � Γ���� ) = �� (�� , �, Ω) + �� (�, �, Γ���� ),�� (�� , �, Ω) = � �� (�, Ω) ��,��� (�, �, Γ���� ) = �������(�, Γ���� )� � ��� .

�� (�� , �, Ω)と�� (�, �, Γ���� )は，それぞれ�, �に対して形状微分（Fréchet 微分）をとっており，�� (�, Ω)は微分の連鎖則を用いることで評価することが可能である．

— 25 —

Shape Optimization による非定常流れ場の安定性制御

2.2 形状更新
本研究において，目的汎関数�(�� , �, Ω)の最小化問題を考え，その時の領域変動�を決定する．

その際，�(�� , �, Ω)の変分である��(�, �, Ω ∪ Γ���� )において，�が未知変数であり�(�, Γ���� )�が trial

function という関係になっている．そこで，�を形状最適化問題の目的関数，�, �を任意の微小パ
ラメーター，として，以下のように反復させて目的関数を最小化させるアルゴリズムを採用する．
a.
b.
c.
d.
e.
f.
g.

���
.
Set ���� , ���� � 0, � ��� , � ��� , � ��� , Ω��� , Γ����

Define the functional ���� , � � , Ω� �.

�
Derive the material derivative ����� , � � , Ω� ∪ Γ�����.

�
Evaluate the sensitivity � � �� � , Γ������.

Obtain the new domain Ω��� to reshape Ω� by Ω��� � �� (Ω� ) � ��� .
Operate �� gradient method [5-7].

Judge convergence.

-If the terminal condition �� ��� � � � � � � � ℝ is satisfied, then stop.
-Otherwise, replace � � 1 with � and return to b.

以下，便宜上，形状更新の反復�を�, �, �, �, �, Ω, Γ���� には明記しない．

2.3 主問題
本研究では，非定常 Navier-Stokes 方程式の初期値境界値問題を解くが，計算時間を短縮する
ために定常 Navier-Stokes 方程式の境界値問題の解を初期条件とする．そして，レイノルズ数，
速度ベクトルと圧力を，それぞれRe, � � �(�), � � �(�)とし，

�(�) � �� � �� (Ω; ℝ� )| � � � on Γ���� , � � (1�� � (� � 1)� , 0) on Γ��� �,
�(�) � {� � �� (Ω; ℝ)| ��(�)� �� � 0 �n Ω�.

Problem 1 (Stationary Navier-Stokes Problem)� � �(Ω), � � �(�), � � �(�)のとき，(�, �): Ω → ℝ� × ℝを求めよ．
(� � �)� � ��� �1
��, � � � � 0 �n Ω,Re� � � on Γ���� ,� � (1�� � (� � 1)� , 0)on Γ��� .

� � �(�) � �� � �� (Ω; ℝ� )� � � � on �Ω�, � � �(�)を�, �に対する trial function とすると，Problem 1 の弱形式は，

— 26 —

SENAC Vol. 49, No. 4（2016. 10）

�

�(�)

�{(� � �)�} � � �1�� � �� − (� � �)� − (� � �)�� �� = 0,Re� �(�, �) � �(�, �, �) � �(�, �) � �(�, �) = 0

1
�� � �� ��, �(�, �, �) = � {(� � �)�} � ��� , �(�, �) = − � (� � �)��� .
� Re
�
�

�(�, �) = �

ここで Problem 1 の解を�� ∈ �(�)とする．

Problem 2 (Non Stationary Navier-Stokes Problem)� ∈ �(Ω), � ∈ �(�), � ∈ �(�)のとき，(�, �): Ω × (0, �) → ℝ� × ℝを求めよ．
�� 1− �� � �� = 0, � � � = 0 �n Ω × (0, �)�� Re
� = � on Γ���� × (0, �)� = (1�� � (1 − �)� , 0) on Γ��� × (0, �)� = �� �n Ω �� � = 0.

そして，Problem 2 の弱形式は，� = {0, � , �}, ∀(�, �) ∈ �(�) × �(�)において，��, �� � �(�, �) � �(�, �) � �(�, �) = 0.
���

時間方向の離散化には有限差分法を用い，移流項の取り扱いには特性曲線法を用いる[18]．その際，��を時間刻み，� = ����を全時間ステップとして，時間を� � = ���(� ∈ {0, � , �})と記述す��

る．その結果，(�, �)における物質微分 の 1 時近似は，�(�, �) = �� を既知関数とすると，
���� − ���� �������(�, �) − �(� − ��, � − ��)(�, � � ) =� �(��) =� �(��),��������� = � − ������ (�), and ���� ���� = ���� (��� ).

最終的に，�� = �� を初期条件，{(�� , �� )}���� ∈ �(�) × �(�)を求めるための弱形式は，
��� − ���� ����, �� � �(�� , �� ) � �(�� , �� ) � �(�� , �� ) = 0.
��� � � � �
その際，{(�� , �� )}���� ∈ �(�) × �(�). ここで，(�, �, �, �) = {(� , � , � , � )}��� , �� , �� ∈ ℝとして，
� ,�以下のように汎関数��� � を定義する．
� ,�
��� � (�� , �, �, �, �, Ω) =��1�� − ���� ����� ��, �� � �(�� , �� ) � �(�� , �� ) � �(�� , �� )�.
���� − ��
����

Shape Optimization による非定常流れ場の安定性制御— 27 —

2.4 形状最適化問題
非定常 Navier-Stokes 問題を主問題とした際の，時間平均運動エネルギー最小化問題��11
� � � �� � �� ����=2 �
�� � ��
����

を定義する．

Problem 3 (Shape Optimization Problem Minimizing the Momentum Energy)� � �(Ω), �, � � �(�) � �(�)のとき，
��

11� � � �� � �� ��� | (�� , �� ) in ������� �, {(�� , �� )}�
min �� =��� in ������� ��.
���(�)
�� � ��
2 �
����

Lagrange 乗数法を用いて Problem 3 に対する Lagrange 関数（汎関数）を以下のように定義す
る．

��

1
1
� ,�
�(�� , �, �, �, �, Ω) =� � � �� � �� ��� � ��� � (�� , �, �, �, �, Ω)�� � ��
2 �
��

����

1
� � � � (�, �, �, �, Ω)��=�� � �������そして，形状微分公式を用いると汎関数の物質微分は，
��(�, �, �, �, �, Ω) = �� (�� , �, �, �, �, Ω) � �� (�, �, �, �, �, ����� ),と書ける．更に，微分の連鎖則を用いると�� (�, �, �, �, �, Ω)は，

�� (�� , �, �, �, �, Ω) = ��,� (�� , �, �, �, �, Ω) � ��,� (�� , �, �, �, �, Ω).

�, �に対して形状微分（Fréchet 微分）をとった際の汎関数を��,� (�� , �, �, �, �, Ω)，�, �に対しては��,� (�� , �, �, �, �, Ω)，�に対しては�� (�, �, �, �, �, Ω)，と記述し，具体的には，
��1� � � � (��, ��, �, �, Ω)��,��,� (�� , �, �, �, �, Ω) =�� � ��
�����

— 28 —

SENAC Vol. 49, No. 4（2016. 10）��1��,� (�� , �, �, �, �, Ω) =� � � � (�, �, ��, ��, Ω)��,�� − ��
�
����

��

1
� �� � � (�, �, �, �, ����� )� � ���� .
�� (�, �, �, �, �, ����� ) =�� − �����������

そして，Problem 3 対に対して随伴変数法を用いると，
��,� (�� , �, �, �, �, Ω) = 0 for ∀(�� , �� ) � �(�) × �(�),��,� (�� , �, �, �, �, Ω) = 0 for ∀(�� , �) � �(�) × �(�),となる．��,� (�� , �, �, �, �, Ω) = 0から Problem 2 を解くことになる．��,� (�� , �, �, �, �, Ω) = 0から，(�, �)に対する随伴変数(�, �)を求める必要がある．

Problem 4(Adjoint Problem for Problem 2)�� � �(Ω)とし，{(�� , �� )}���� � �(�) × �(�)が与えられたとき，(�, �): Ω × (�, 0) → ℝ × ℝを求め

よ．

−

�� 1
− �� � �� � � � � = 0, � � � = 0 �� Ω × (�, 0)
�� Re
� = �� �� Ω �� � = ���.

� = � o� �Ω × (�, 0),� � �
Problem 2，Problem 4 から求められた{(�� , �� )}���� � �(�) × �(�)と{(� , � )}��� � �(�) ×

�(�)を��(�, �, �, �, �, ����� )に代入すると，Problem 3 に対する感度は，��1� �� � � (�, �, �, �, ����� )� � ������(�, �, �, �, �, ����� ) =�� − ������������ � (�, �, �, �, ����� ) =1(��� : ��� )|��������
Re

��(�, �, �, �, �, ����� ) = 0となる場合の local minimizer である�を見つけることが出来，�を任意の微小パラメーターとすれば，�(Ω) = �� (Ω) � ��によって最適な形状が得られる．しかし，本研究において，主問題は非線形偏微分発展方程式であるので数学的に� � (�, �, �, �, ����� )を評価することは困難である．そこで，セクション 2.3 で記述した形状更新の反復で� � (�, �, �, �, ����� )を数値的に評価する．

— 29 —

Shape Optimization による非定常流れ場の安定性制御

2.5 線形安定性解析
線形安定性理論に基づき，固有値方程式を導出する．以下，固有値と固有関数を，それぞれ
�(�) × ��(�)とし，
� , �̂ ) � �λ, (��(�) = {�� � �� (�(Ω); ℂ� )| � = � �n �Ω},���(�) = {�̂ � �� (�(Ω); ℂ)| ��(�)�̂ �� = 0 in Ω}.

Problem 5(Linear Stability Analysis)� , �̂ ): Ω → ℂ × ℂ� を求めよ．
� � �(Ω)とし，�� � �(�)が与えられたとき，λ, (��−� + (�� ⋅ ∇)�� + (�� ⋅ ∇)���1
� + ��
� = 0, � ⋅ � = 0 in Ω��Re

� = � �n �Ω.

�(�) × ��(�)を(��, ��) � �� , �̂ )に対する trial function とすると Problem 5 の弱形式は，
∀(�
��

�(�)

� , �) + �(�� , �, �) + �(�� , �� , �) + �(��, �̂ ) + �(�� , ��) = 0.
� ⋅ ��� + �(�3. 数値計算スキーム
本研究における数値計算については，無償の有限要素法ソフトウェアである Freefen++[18]を用いる．数値計算アルゴリズムは以下の通りである．
I.
II.

Prepare a stationary velocity �� by solving Problem 1 with the Newton-Raphson

method at Re = 11500.

Set the stationary velocity �� for the initial domain and conditions, and put � = 0.

III.

Solve Problem 2 and Problem 4, Evaluate the sensitivity,

IV.

Reshape Ω by �(Ω) = �� + �.

V.

Judge convergence.

VI.

-If the terminal condition �� ��� − � � � � 0.01 is satisfied, then proceed to Ⅶ and stop.

-Otherwise, replace � + 1 with � and return to Ⅲ.

Solve Problem 5, and compare linear stability in the initial and optimal domain.

上記の数値計算アルゴリズムを進めるにあたって，空間及び時間方向の離散化を行う．空間方
向には有限要素法により三角形分割を行う．その際，速度，圧力には P2/P1 要素を用いる．また，
時間方向には時間刻み�� により有限差分法を用いている．空間及び時間方向の離散化精度につい— 30 —

SENAC Vol. 49, No. 4（2016. 10）ては以下に記述する．空間方向の離散化精度を担保するために，Γ��� において�� = (1,0)としてNode と Element を増加させ，Problem 1，Problem 5 を解き最大固有値を比較し，Table1 に示した．次に，Node と Element が(8281, 16200)の結果と，[19-23]との結果を比較した．特に，[23]との比較では相対誤差がそれぞれ 1.73%，2.86%となった．そのため，本研究では，Node とElement を(8281, 16200)とする．次に，時間刻み���については，0.001 として，Problem 1，
Problem 2 を解き運動エネルギーを Fig. 1 に示した．Fig. 1 の結果から，��は十分に小さいと仮定する．
Table1. ��� and Imag[�] depending on (♯ of Nodes, ♯ of Elements)(♯ of Nodes, ♯ of Elements)(2601, 5000)

Imag[λ]

Re�

7310

±3.02

7515

±2.95

(5041, 9800)

7594±2.916(6501, 12800)

7741

±2.89

(8281, 16200)

7890

±2.86

(3721, 7200)

Table2. ��� and Imag[�] in [19-23].

References

Imag[λ]Re�

[19]7763.4±2.8634i

[20]

8000

±2.8356i[21]8375.0±2.7640i

[22]

7890

±2.7646i[23]8026.6±2.82559i
0.04734

Momentum Energy

Momentum Energy

0.05

0.049

0.048

0.047

0

100

200

300

400

500

600

time

(a)�� = (0,1000).

700

800

900 1000

0.04733

0.04732

600

601

602

603

604

605

606

607

time

(b) � = (600,61).

Fig. 1. The momentum energy, (a) � = (�, ����) and (b) � = (���, ��).

608

609

610

— 31 —

Shape Optimization による非定常流れ場の安定性制御

4. 数値計算結果
セクション 3 で示した，数値計算アルゴリズムに従がって形状最適化問題を解いた．Fig. 2 と� の流れ関数の分布である．Fig.
Fig. 3 は，それぞれ初期形状と最適形状における定常解�� の擾乱�

3 からΓ���� が領域の内部に向かって移動していることがわかる．Fig. 4 は，初期形状と最適形状の最大固有値の実部を図示しており，最適形状の方がRe� が小さくなった．Fig. 5 は，Re = 11500における固有値の分布を示しており，Fig. 4 の結果を補強している．

(a)

(b)

Fig. 2. Stream functions of (a) � and (b) � in the initial domain.
�(a)Fig. 3. Stream functions of (a) �� and (b) � in the optimal domain.
7

Ω0
Ω1

0.06
0.05
0.04
0.03
0.02
0.01
0

5
4
3
2
1

-0.01
-0.02
6500

Ω0
Ω1

6

Imaginary parts

Real parts of the eigenvalues

0.07

(b)

0

7000

7500

Re

8000

8500

9000

Fig. 4. Real parts of the leading eigenvalues

-0.2

-0.15

-0.1

-0.05

Real parts

0

0.05

0.1

Fig. 5. Spectrum of the eigenvalues
at Re = 11500.

— 32 —

SENAC Vol. 49, No. 4（2016. 10）

5. 考察
本研究では，運動エネルギー最小化問題を行い，初期形状と最適形状における線形安定性を比較する．その際，初期形状を 2 次元 Cavity 流れとし，レイノルズ数Re = 11500において非定常
Navier-Stokes 問題を主問題とする．数値計算には無償の有限要素法ソフトウェアである
Freefen++[18]を用いた．

数値計算の結果，Γ���� が領域の内部に向かって移動していた．更に，初期形状と比較して最適形状における臨界レイノルズ数Re� が減少したことが数値的に示された．一般に，流れ場は剪断応
力が大きななった際に不安定になるとされている．目的関数を最小化する過程でΓ���� が領域の内部に移動することによってΓ���� 上で剪断応力が大きくなり，流れ場が不安定になったと考えられる．その結果，運動エネルギーを最小化しているにも拘らずRe� が減少した．

本研究は，2次元Cavity流れを対象としているが，将来的には，3次元の大規模高精度計算により複雑形状内流体の形状最適化を行う予定である．しかし，そのためには多くの計算コストが必要となるため，今後は東北大学サイバーサイエンスセンターのスーパーコンピュータを利用することで実現することが可能であると考えている．

謝辞
今後の研究にのために，東北大学サイバーサイエンスセンターのスーパーコンピュータ関係各位に，Freefem++をインストールして頂いた．これにより， Problem 2 を実行するに至った．現
在は，2 次元空間における数値解析しか行ってないが将来的には 3 次元空間の解析に役立てていきたい．更に，情報科学研究科とサイバーサイエンスセンターとの申し合わせにより研究室に負担金が生じることなく大規模な計算機が利用できる仕組みがあるので，大規模高精度な計算に積極的に利用したい．これは，著者の今後の研究にとって非常に有意義なことであり，感謝申し上
げたい．

参考文献
[1] J. Hadamard, Memoire sur un probleme d’analyse relatif a l’equilibre des plaques elastiques
encastrees, Memoire des savants etragers, Oeuvres de J. Hadamard, CNRS, Paris, 1968
[2] O. Pironneau, On optimum profiles in Stokes flow, JFM, 59, 117-128, 1973
[3] O. Pironneau, On optimum design in fluid mechanics, JFM, 64, 97-110, 1974
[4] O. Pironneau, Optimal shape design for elliptic systems, Springer-Verlag, New York, 1984.
[5] H. Azegami, S. Kaizu, M. Shimoda, E. Katamine, Irregularity of shape optimization prob-lemsand an improvement technique, Computer Aided Optimization Design of Structures V, Chap. 4,
Sec. 4, edited by Hernandez, S. and Brebbia, C. A., Computational Mechanics Publications,Southampton, 309-326, 1997
[6] H. Azegami, Z. Wu, Domain Optimization analysis in linear elastic problems: Approach using
traction Method, JSME Int. J. Series A, 39, 272-278, 1996

Shape Optimization による非定常流れ場の安定性制御— 33 —[7] H. Azegami, Regularized Solution to Shape Optimization Problem, Transactions of the Japan
Society for Industrial and Applied Mathematics, 23, 83-138, 2014[8] E. Katamine, Y. Nagatomo, H. Azegami, Shape optimization of 3D viscous flow fields, Inverse
Problems in Sci. and Eng., 17, 105-114, 2009.
[9] E. Katamine, H. Azegami, T. Tsubata, S. Itoh, Solution to shape optimization problems ofviscous flow fields, Int. J. of Com. Fluid Dyn., 19, 45-51, 2005.
[10] T. Nakazawa, Eiji Katamine, Akiyasu Tomoeda, Hideyuki Azegami，An application of shape
optimization method to a room design based on anincompressible fluid model, RIMS講究録別冊「トポロジー最適化理論とその現実問題への応用」, B(54), 1-15, 2015.
[11] T. Nakazawa, Increasing the Critical Reynolds number by maximizing dissipation energy
problem, Proc. of the 5th ICJWSF2015, Editor Antonio Segalini, Chapter 75, 2016[12] T. Nakazawa, Shape Optimization of the maximizing problem of the dissipation energy and its
effect on Hydrodynamic stability, Proc. of the 15th ETC2015, 477, Netherland, 2015[13] Y. Iwata, H. Azegami, T. Aoyama, E. Katamine, Numerical solution to shape optimization
problems for non-stationary Navier–Stokes problems, JSIAM Letters, 2, 37-40, 2010[14] D. Li, R. Hartman, Adjoint-based airfoil optimization with discretization error control, Int. J.
Numer. Meth. Fluids, 77, 1-17, 2015[15] S. Schmidt, C. Ilic, V. Schulz, N. Gauger, Three-dimensional large-scale aerodynamic shape
optimization based on shape calculus. AIAA Journal, 51, 2615-2627, 2013[16] M. Kimura, Lecture notes Volume IV, Topics in Mathematical Modeling, 1-38, 2008[17] H. Notsu, M. Tabata, Error estimates of a pressure-stabilized characteristics finite element
scheme for Oseen equations, J. Sci. Comput, 65, 940-955, 2015[18] Hecht, F. New development in FreeFem++. J. Numer. Math. 20, no. 3-4, 251–265. 65Y15, 2012[19] Poliashenko M, Aidun CK. A direct method for computation of simple bifurcations, J.l ofCom. Physics, 121, 246–260, 1995

[20] Fortin A, Jardak M, Gervais JJ, Pierre R. Localization of Hopf bifurcations in fluid flow
problems. Int. J. for Num. Methods in Fluids, 24, 1185–1210, 1997[21] Tiesinga G, Wubs FW, Veldman AEP. Bifurcation analysis of incompressible flow in adriven cavity by the Newton–Picard method, J. of Comp and Appl. Math., 140:751–772,
2002[22] Cadou JM, Potier-Ferry M, Cochelin B. A numerical method for the computation ofbifurcation points in fluid mechanics. Euro. J. of Mech. Fluids B, 25, 234–254, 2006[23] V. B. L. Boppana and J. S. B. Gajjar, Global flow instability in a lid-driven cavity, Int. J.
Numer. Meth. Fluids, 2009

— 34 —

SENAC Vol. 49, No. 4（2016. 10）

[報 告]
［報 告］

学際大規模情報基盤共同利用・共同研究拠点（JHPCN）
第8回シンポジウム報告
江川隆輔
スーパーコンピューティング研究部
平成28年7月14日，15日の両日，学際大規模情報基盤共同利用・共同研究拠点第8回シンポジウムが，例年通り東京品川The Grand
Hallにおいて開催されました．
学際大規模情報基盤共同利用・共同研究は，東北大学，北海道大学，東京大学，東京工業大学，名古屋大学，京都大学，大阪大学，
九州大学にそれぞれ附置するスーパーコンピュータを持つ8つの共同利用の施設を構成拠点とする「学際大規模情報基盤共同利用・共同研究拠点」が, 超大規模数値計算系応用分野，超大規模データ処
理系応用分野，超大容量ネットワーク技術分野，およびこれらの技術分野を統合した大規模情報システム関連研究分野，更には分野間に亘る複合分野の研究者らと取り組む学際的な共同利用・共同研究です．今年度から，これらの一
般共同研究枠に加えて，国際共同研究，企業共同研究，萌芽的研究枠が設置され，今後将来に向けて，我が国の学術・研究基盤の更なる高度化と恒常的な発展が大いに期待されております．これまで平成22年度から27年度の6年間において233件の課題が学際大規模情報基盤共同利用・共同研究として採択されており，そのうち当センターとの共同研究課題は，38件ありました．平成28年度は国際共同研究3件を含む39採択課題のうち，6件が当センターとの共同研究（うち1件が国際共同研
究）として採択されております．
本シンポジウムでは，昨年度採択された36件に及ぶ研究課題の成果発表，および今年度採択された39件の研究課題のポスター発表を通して積極的な議論が交わされたばかりでなく，拠点および拠点を構成する各センターにおけるそれぞれの共同研究への取り組みを確認することができました．
本報告では，今年度の東北大学との共同研究として採択された課題の成果は，JHPCNウェブページ
（http://jhpcn-kyoten.itc.u-tokyo.ac.jp/ja/adoption.php）にて公開されております．ご興味のある方には是非，先進的，且つ多様な共同研究成果をご覧いただければと存じます．
来年度も引き続き，学際大規模情報基盤共同利用・共同研究拠点の研究公募行う予定でおります．
ご興味のあるかたは奮ってご応募ください．当センターの計算機科学を専門とする教員との共同研
究の可能性を検討したい，手続き方法が分からない等，本応募に関して不明な点などあります場合
は，お気軽に当センターまでお問い合わせください．
【JHPCNに関する問い合わせ窓口】
joint̲research@cc.tohoku.ac.jp

— 35 —

SENAC Vol. 49, No. 4（2016. 10）

[報 告]
［報 告］

サイバーサイエンスセンターオープンキャンパス報告
小松 一彦
スーパーコンピューティング研究部

今年も東北大学 オープンキャンパスが開催されました．進学を考えている学生や，東北大学の研究や活動に興味・関心を持っている一般の方々に向けて，年に1回開催されており，今年は7月27日と28日の2日間に渡り開催されました．サイバーサイエンスセンターでは，「スパコンを見に行こう！」「ネットワークを知ろう！」の2つのキャッチフレーズを掲げ，大規模科学技術計算を支える最新型スーパーコンピュータ，計算結果を分かりやすく可視化する3次元可視化システム，
東北の大学ネットワークを支えるネットワーク機器，大規模計算機の歴史的資料を展示している
コンピュータ博物館を一般に公開しました．また，ネットワーク研究部，スーパーコンピューテ
ィング研究部，先端情報技術研究部の研究成果を展示しました．
津波被害予測，熱中症予測などの大規模シミュレーションを支えるスーパーコンピュータや，
大学のインフラを支えるネットワークへの関心も高く，1700人以上の方々にご来場いただきました．計算機専用棟(サイバーサイエンスセンター2号館)に導入された最新型スーパーコンピュータ
SX-ACEや関連設備を見学し，スーパーコンピュータの大きさやその音，冷却するための冷風などを肌で感じていただきました．また，50インチモニタ12面からなる可視化システムを用いて，スーパーコンピュータでの計算結果の3次元可視化結果を体験していただきました．コンピュータ博
物館では，歴史的価値の高いスーパーコンピュータやネットワーク機器，資料などを見学して，スーパーコンピュータの歴史を巡りました．これらを通じて，スーパーコンピュータやネットワークが研究だけでなく，実生活にどのように役立っているのかを感じてもらえたと思います．研
究開発部の展示では，情報セキュリティ・スーパーコンピュータ・サイバー医療に関する研究活動だけでなく，大学生活や大学での講義などについても，職員や学生と交流しておりました．来
年度のオープンキャンパスも一般に公開する予定です．みなさまのお越しを心よりお待ちしております．

— 36 —

SENAC Vol. 49, No. 4（2016. 10）[報 告] 計算科学・計算機科学人材育成のためのスーパーコンピュータ無償提供利用報告［報 告］計算科学・計算機科学人材育成のためのスーパーコンピュータ無償提供利用報告

東北大学大学院データサイエンスプログラムにおける
疑似プロジェクト演習
中尾光之・瀬川悦生・山田和範
東北大学大学院情報科学研究科
1. 概要
2016 年 7 月 27 日から 29 日の 3 日間において，東北大学大学院情報科学研究科が主催する
データサイエンスプログラムにて，講義の一環として DSP Training camp II および Big data
challenge を行いました．この実習では，修士課程の学生がトレーニングを受ける対象，博士課程の学生がトレイナーとして，現実社会に存在する，ビッグデータを用いた疑似的な研究プ
ロジェクトに取り組みました．
授業のテーマのひとつがビッグデータであり，大量のデータを扱う必要がありましたので，
サイバーサイエンスセンターによって提供されるスーパーコンピュータ無償提供制度を利用させていただきました．
実際の実習では，受講者 18 名が 3 つのチームに分かれ，計算言語学，生命科学，経済学に関連する問題を，提供していただいたスパコンをはじめとする計算機を駆使することで解決しました．今回，実習で扱った研究テーマは答えが用意されていない発見的なものでしたが，チームによっては，これまでに明らかにされていなかったような結果を出すことができ，実習が終わった現在においても，さらに発展的に研究を進めています．
以下では，その内容について簡単にご報告いたします．
2. 取り組んだ問題取り組んだ問題は以下の 3 つです．
A) 言語学の問題：1 年間にわたり取得された日本語のツイッターの時系列データから，1 日
毎における社会的な雰囲気を抽出し，ベクトル化，その社会的な雰囲気が，日本の主要会
社の株価，日経平均 225 のリーディングインジケーターとなり得るかどうかを検証しました．また，リーディングインジケーターとなり得る場合，その時系列データを利用して，株価の変動の予測器を構築することとしました．
B) 生命科学の問題：タンパク質立体構造中にディスオーダー領域と呼ばれる，水中で他の確固たる構造部位より変動が大きい部位があります．この部位は何らかの生物学的な機能を

東北大学大学院データサイエンスプログラムにおける疑似プロジェクト演習— 37 —C) 有する場合が多く，タンパク質立体構造科学的に重要な領域です．この領域を，これまでに得られた既存のデータから予測する予測器を，深層学習法を用いて構築しました．
D) 経済学の問題：株価の変動を意味するボラティリティという指標があります．ボラティリ
ティが大きいということは，株式取引を行う際に，その分，損をする可能性は大きくなる
ものの大きなキャピタルゲインを得ることもできることを意味します．その点において，
ボラティリティを精度よく予測するモデルは重要です．この問題では，これまでに提案されているいくつかのボラティリティのモデルを日経平均に適用し，どれが最も良くボラティリティを説明するのかを検証しました．
3. 実習の様子および成果
各チームとも，チーム内で役割分担を決め，効率的に仕事を進めました．特に扱うデータ量が大きかったのは，言語学の問題を扱ったチーム A でした．圧縮されたツイッターのデータだ
けで約 700GB です．その処理にはスパコンの 1 ノードをフル活用したにもかかわらず，約 21
時間もの計算時間が必要でした．このような処理はサイバーサイエンスセンターからスパコンを提供していただかなければ為し得なかったことでありましたので，今回使わせていただいて良かったと思います．その他のチームも慣れないスパコンのジョブ管理システムの使い方に戸惑いながらも，何とか使いこなし，与えられた問題のみならず各チームが発展的な研究成果を出しました．
また，得られた全ての研究成果は2016 年 8 月 9 日に東北大学とオハイオ州ケースウエスタンリザーブ大学との共同ワークショップにて発表しました．
4. 最後に本実習は，サイバーサイエンスセンタ

実習の様子．各々の端末からスパコンに SSH アクセスをしています．

ーのスーパーコンピュータを利用する
ことではじめて実現することができました．また，実習中には技術的なサポートを頂きました
ことを感謝申し上げます．

— 38 —

SENAC Vol. 49, No. 4（2016. 10）

[報 告]
［報 告］
プレスリリース報告
平成 28 年 7 月 25 日、以下の内容でプレスリリースしましたのでご報告いたします。

気象データと連動した熱中症リスク評価システムの開発― 幼児 60 分で熱中症リスク、真夏にアスファルトの歩行 ―名古屋工業大学
東北大学サイバーサイエンスセンター
一般財団法人 日本気象協会【概要】
名古屋工業大学、東北大学サイバーサイエンスセンター、一般財団法人日本気象協会の共同研究
グループは、熱中症リスク評価シミュレーション技術に、気象予報データと経験から得られた数式を融合させたデータを組み込み、現実的な条件（例えば、アスファルト、運動場など）での熱中症リス
ク評価システムを開発しました。本システムでは、日本の夏場を想定し、屋外を歩行した場合、幼児の体温上昇は成人よりも早く、それに伴い発汗量も多く、60 分間で熱中症（軽度の脱水症状）になるリスクを有する結果を得ることができました。今後、個人属性を考慮した適切な熱中症リスク評価
技術を活用することで、今まで以上に場面に応じた発症数の低減に貢献することが期待できます。
【これまでの研究・開発の経緯】
名古屋工業大学、東北大学サイバーサイエンスセンター、一般財団法人日本気象協会の共同研究グ
ループは、1）乳幼児や高齢者などの個人特性を考慮した熱中症リスク評価のための複合物理・シス
テムバイオロジー統合シミュレーション技術を東北大学サイバーサイエンスセンターが有するスー
パーコンピュータ SX-ACE に効率的に実装、高速化、2）気象予報データと融合させることによる、個人特性を考慮した 3 時間後の熱中症のリスクを 10 分で評価する技術の開発に成功してきました（2015年 7 月 21 日プレスリリース）。
【新規性】
本リリースでは、これまで開発してきたシステムを改良し、アスファルト、運動場など気温・湿度の測定データより経験的に導出した数式と一般財団法人日本気象協会が提供する気象予測データを併用することにより、ある特定の環境下での気象データ（外気温、湿度）を活用し、熱中症リスク評価に成功しました。また、システムバイオロジーあるいは気象データの効率的な読み込みなどの高速化、最適化を実施しました。
開発したシステムを、歩道(アスファルト)上で 60 分間の歩行をした場合を想定したリスク評価に応用しました。具体的には、周辺の実際の気温・湿度の測定データ（特に、アスファルトからの輻射熱を考慮するために、高さ方向で温度が変化）に基づく（図１参照）推定値を用いて、成人男性および幼児の計算モデル実験を実施しました。その結果、外気温が約 34℃の環境下で、幼児の体温上昇
値は 1.12℃、成人の値 0.61℃に比べ、約 2 倍となりました（体表面温度分布．図２参照）。また 60分間の総発汗量は成人体重の 0.34%だったのに対し、幼児は体重の 2.3%に達し、初期の脱水症状になるリスクを有することがわかりました。これは成人と幼児の体形および生理的相違のみならず、歩道からの照り返しによる成人と幼児での周辺温度の相違によるものです。この手法を用いれば、ある特定の環境における将来のリスク評価が可能となります。熱中症高リスク群（高齢者、幼児など）と若者の感覚の違いなどを把握でき、周囲の気配りなどを促すのに利用することが期待されます。

— 39 —

気象データと連動した熱中症リスク評価システムの開発〜幼児 60 分で熱中症リスク、真夏にアスファルトの歩行〜

13:00
13:10
13:20
13:30
13:40

地面からの高さ[cm]

150

100

50

0

35

36

37

38

39

。

測定気温[ C]

40

41

42

図 1：アスファルト上での気温測定例（2016 年
7 月 14 日．名古屋）。高さ方向の気温は、1.5mにおける測定値から推定可能。

図 2：1 時間路上を歩行後の体表面温度【今後の展開】
本技術により、より現実に即した環境下でのリスク評価が可能となります。一例として、一般財団
法人日本気象協会の気象予測を用いて、2016 年 7 月 26 日から 28 日までの 12：00～13:00 において成人および幼児が屋外を散歩した場合の、体温上昇および発汗量予測値を表１に示します。
表１. 2016 年 7 月 26 日～28 日までの成人および幼児の体温上昇、発汗量予測値(12:00～13:00)
7 月 26 日
成人

7 月 27 日

体温上昇(温上0.31
総発汗量(g)129
幼児
体温上昇(℃)
0.47
総発汗量(g)193（参考） WBGT
26*日照時間の割合を 100%と仮定して算出した値。本例は仙台

7 月 28 日
0.31
130
0.48
195
26

0.33
140
0.55
209
27

【この件に関するお問合せ】
熱中症のリスク評価について
名古屋工業大学
電気・機械工学専攻
教授 平田晃正
Tel: 052-735-7916
e-mail: ahirata@nitech.ac.jp

スーパーコンピュータ・プログラムの最適化について
東北大学サイバーサイエンスセンター准教授 江川隆輔
Tel: 022-795-3416
e-mail:egawa@cc.tohoku.ac.jp
東北大学情報部情報基盤課
総務係 佐藤恵美子
Tel: 022-795-3407
e-mail:som@cc.tohoku.ac.jp【参考情報】
日本気象協会ホームページ：http://www.jwa.or.jp/news/2015/07/post-000536.html(平成 27 年 7 月 23 日公開)
【謝辞】本研究の一部は、学際大規模情報基盤共同利用・共同研究拠点の支援によるものです。

— 40 —

SENAC Vol. 49, No. 4（2016. 10）

[報 告]
［報 告］
プレスリリース報告
平成 28 年 9 月 30 日、以下の内容でプレスリリースしましたのでご報告いたします。

血行状態モニタリング装置「魔法の鏡」の開発に成功
東北大学 サイバーサイエンスセンター 先端情報技術研究部
東北大学革新的イノベーション推進機構
科学技術振興機構（ＪＳＴ）【ポイント】●遠隔・非接触的に人体の皮膚表面の血行状態をリアルタイムにわかりやすく動画像で表示します。
●心拍数ばかりでなく、身体の 2 か所の脈波伝搬時間を連続算出します。
●外付け Web カメラや内蔵カメラが付いたパソコンで動作します。
●家庭の鏡や自動車内のドライブレコーダなどへ応用することで、日常的な体調管理ができると期待されます。
国立大学法人東北大学（以下、東北大学）は、平成 25 年度 国立研究開発法人 科学技術振興機構のセンターオブイノベーション（COI）創出プログラム注)「さりげないセンシングと日常人間ドックで実現する理想自己と家族の絆が導くモチベーション向上社会創生拠点」（代表:
高山卓三・末永智一）の支援を受け、血行状態モニタリング装置「魔法の鏡」の開発に成功しました。
本装置は、2016 年 10 月 1 日（土）に仙台国際センターで開催される第 39 回日本高血圧学会「市民公開講座」会場入り口で展示されます。
注) センター・オブ・イノベーション（COI）プログラム
国立研究開発法人 科学技術振興機構による公募型研究開発プログラム。ハイリスクではあるが実用化の期待が大きい異分野融合・連携型の基盤的テーマに対し、集中的な支援を行っている。産学が連携する研究開発チームを形成し、イノベーションを推進する積極的な社会実装を目指しています。本プログラムは、文部科学省の施策革新的イノベーション創出プ
ログラム（COI STREAM）に包摂され、一体的に推進するものです。COI STREAM では、現在潜在している将来社会のニーズから導き出されるあるべき社会の姿、暮らしの在り方を見据えたビジョンを設定しています。COI プログラムでは、このビジョンを基に１０年後を見通した革新的な研究開発課題を特定した上で、既存分野・組織の壁を取り払い、企業だけでは実現できない革新的なイノベーションを産学連携で実現します。
連絡先：東北大学 サイバーサイエンスセンター 先端情報技術研究部
教授

吉澤

誠

〒980-8579 仙台市青葉区荒巻字青葉 6-6-05
電話：022-795-7128
E-mail：yoshizawa@cc.tohoku.ac.jp

血行状態モニタリング装置「魔法の鏡」の開発に成功— 41 —【開発の背景】
少子高齢化に伴う医療費・介護費の増大に対処するためには、病気に至る前の未病対策が重要で、家庭や職場において常時健康状態のチェックが可能であることが理想的です。
これを目的として、既にリストバンド型のウェアラブルセンサなどで常時身体の状況を記録する装置が市販されています。しかしそれらのほとんどは、接触式センサによって心拍数や活動度などの健康に関連すると思われる既存の指標を得るものです。しかし、特別なセンサを常時身に付けることは煩わしく、毎日意識して機器を操作する必要があるような健康管理法は習慣化しにくいものです。
これに対して、洗面所や脱衣所にある鏡は多くの人々がほぼ毎日利用します。特に、洗面所や鏡台の鏡を使って自分の顔を毎日入念にチェックする人が少なくありません。また、風呂場などの脱衣所の鏡では衣服を必ず脱ぐために顔に限らず全身をチェックすることが可能です。
一方、自分自身の体調を客観的かつ簡単に知ることは容易ではありません。特に、多くの人々が悩んでいる冷え性・肩こり・血行不良・むくみ・だるさ・食欲不振などの不定愁訴のような自律神経に関係する症状は、自覚することはできてもこれを客観的に把握することは困難です。
【目的】
そこで本研究では、ビデオカメラとコンピュータを内蔵した鏡型ディスプレイの前に立つだけで、自律神経指標に基づいたその日の健康予報を使用者に直感的で分かりやすく表示するツールとして、
いわば「魔法の鏡」を実現することを目指しています。

【方法】
われわれはこれまで、脈波信号の計測によく利用される光電脈波計の代わりに、ビデオカメラで撮影した身体映像から、皮下の血液中のヘモグロビンが吸収する緑色信号に基づいて、脈波伝搬時間の推定を行う技術を開発し、これが血圧変動と正の相関をすることを明らかにしています[1]。本研究では、この技術に基づいて上図のような鏡型ディスプレイを構築しました[2]。
まず、身体映像の領域をモザイク状の小領域に分割し、各領域の緑色信号のうち心拍周波数近傍の成分が強いものだけを対象として選択し、心拍変動に無相関な運動や周辺光変化による雑音成分をリ
アルタイムにキャンセルするアルゴリズムを開発しました。さらに、映像脈波情報から脈波伝搬時間

— 42 —

SENAC Vol. 49, No. 4（2016. 10）

あるいは血行状態を推定するために、身体の異なる 2 箇所の領域間の信号の位相差を抽出します。各
モザイク領域を変遷する 2 次元的な血行パターンを、顔などの実映像に重畳して表示し、これを利用者が見ることにより、その変化と自分自身の体調の変化とを比較することができるようになります。
これにより、洗面所や鏡台などに取り付けた「魔法の鏡」では、毎朝身支度を整えたり化粧をしたりする際に、自分の顔に重畳する血行状態を観察することが可能になります。また、風呂場などの脱衣所に取り付けた「魔法の鏡」では、ビデオカメラを利用者の後方に取り付けることにより、普段は見えない利用者の肩や背中の血行状態が観察できるようになります。
【仕様】
 動作コンピュータ:

Windows マシン

 ビデオカメラ：USB 接続 Web カメラ，コンピュータ内蔵カメラ，SONY 製 Playstation Eye【今後の展開】
今後は、心拍数変動や脈波振幅変動などから得られる自律神経系指標を表示する予定です。これと併せて身体の血行状態を毎日チェックすることで、自覚した体調の変化との関係を利用者自身で学習して行くことができ、その日の体調を予測するための手掛かりを把握できるようになると予想されます。また、次のような応用が期待されます。
 インターネットを介した遠隔からの見守り用体調管理
 自動車内体調監視
 アスリートの体調推定
 テレビ視聴者の感情分析 デジタルサイネージの評価
 ロボットによる応対者の体調・感情推定【文献】[1] Norihiro Sugita, Kazuma Obara, Makoto Yoshizawa, et al.: Techniques for estimating bloodpressure variation using video images, Proc. of 37th Annual International Conference of
the IEEE Engineering in Medicine and Biology Society (EMBC 2015), (2015).
[2] Makoto Yoshizawa, Norihiro Sugita, Makoto Abe, et al.: Blood perfusion display based onvideo pulse wave, Proc. of 38th Annual International Conference of the IEEE Engineeringin Medicine and Biology Society (EMBC 2016), (2016).

— 43 —

SENAC Vol. 49, No. 4（2016. 10）[Web 版大規模科学計算システムニュースより]
大規模科学計算システムニュースに掲載された記事の一部を転載しています。

http://www.ss.cc.tohoku.ac.jp/tayori/

コンパイラのバージョンアップについて（No.227）
2016 年 8 月 28 日に FORTRAN90/SX コンパイラ、Intel コンパイラのバージョンアップを行いますのでお知らせいたします。詳細につきましては、リリースメモをご覧ください。
対象システムコンパイラ・ライブラリ名

旧バージョン

新バージョン

リリースメモ*SX-ACEFORTRAN90/SX

Rev.520

Rev.530FORTRAN90/SX リリースメモ.pdf

ASL

24.0 リビジョン 2

24.0 リビジョン 3

ASL リリースメモ.pdf

Intel Fortran Compiler

15.0

Update 2

16.0

Update 3Intel Fortran Release Notes

Intel C++ Compiler

15.0

Update 2

16.0

Update 3Intel C++ Release Notes

MKL

11.2

Update 2

11.3

Update 3

MKL Release Notes

MPI

4.1

Update 3

5.1

Update 3

MPI Release Notes.pdf

LX406Re-2＊リリースメモは以下からご参照ください。
http://www.ss.cc.tohoku.ac.jp/tayori/tayori/227.html
なお、コンパイルコマンドに変更はありません。オプションの詳細は、sxman コマンドや PDF 版マニュアルで参照できます。参照方法は、以下をご覧ください。
http://www.ss.cc.tohoku.ac.jp/super/online_manual.html（共同利用支援係，共同研究支援係）

数式処理プログラム Mathematica のバージョンアップについて（No.228）
数式処理プログラム「Mathematica」のバージョンアップを行いますのでお知らせいたします。
Mathematica は Stephen Wolfram によって作られた、プログラミング言語を備えた数式処理システ
ムです。Mathematica の機能は、数値計算、記号計算、グラフィックスという 3 つに大別でき、この3 つが一体となって使いやすいインターフェイスを提供しています。
バージョン名バージョンアップ日
サービスホスト起動コマンド: Mathematica 10.4
: 2016 年 8 月 29 日（月）: front.cc.tohoku.ac.jp（並列コンピュータ）: mathematica （GUI 版）: math （テキスト版）

新機能の概要、機能の詳細、実行方法については開発元 HP、およびセンターHP をご覧ください。
開発元 HP
http://www.wolfram.com/mathematica/new-in-10/
センターHP
http://www.ss.cc.tohoku.ac.jp/application/mathematica.html（共同利用支援係）

— 44 —

SENAC Vol. 49, No. 4（2016. 10）

科学技術計算言語 MATLAB のバージョンアップについて（No.228）
科学技術計算言語「MATLAB」のバージョンアップを行いますのでお知らせいたします。
MATLAB は、アルゴリズム開発、データの可視化、数値計算を行うための高レベルなテクニカルコ
ンピューティング言語と対話型環境です。 MATLAB を利用することにより、C、C++、Fortran といった伝統的なプログラミング言語よりも短時間で科学技術計算の問題を解決することが可能です。
バージョン名バージョンアップ日
サービスホスト起動コマンド: MATLAB R2016a (Ver.9): 2016 年 8 月 29 日（月）: front.cc.tohoku.ac.jp（並列コンピュータ）: matlab（GUI 版）: matlab -nosplash -nodesktop -nodisplay（テキスト版）

新機能の概要、機能の詳細、実行方法については開発元 HP、およびセンターHP をご覧ください。
開発元 HP
http://www.mathworks.co.jp/products/new_products/latest_features.html
センターHP
http://www.ss.cc.tohoku.ac.jp/application/matlab.html（共同利用支援係）

MSC Software 社アプリケーションの利用について（No.228）

現在センターの並列コンピュータでは Marc/Mentat を提供していますが、以下のアプリケーションを別途研究室のワークステーションや PC 等にインストールして利用可能です。
対象者は東北大学に所属の方（教職員・学生等）で、センターの利用者番号が必要です。動作環境等については MSC Software 社の HP をご参照ください。
利用可能アプリケーション• Marc/Mentat（センターの並列コンピュータでも実行できます）• MSC Nastran• Patran• Dytran• FlightLoads• Sindaアプリケーション情報（MSC Software 社）
http://www.mscsoftware.com/
動作環境（MSC Software 社）
http://www.mscsoftware.com/support/platform-support

— 45 —

Web 版大規模科学計算システムニュースよりご利用方法については、共同利用支援係までお問合せください。
東北大学サイバーサイエンスセンター
情報部情報基盤課 共同利用支援係
Tel: (022)795-6251
E-mail: uketuke@cc.tohoku.ac.jp（共同利用支援係）

バッチリクエストの状態確認コマンドについて（No.229）
バッチリクエストの状態確認コマンドとして、qstat コマンドが利用できるようになりましたのでご案内いたします。
qstat コマンドは、投入されたリクエストの状態を表示します。状態は STT 項目に表示されます。
リソースに空きがなければ実行待ち状態（QUE）になり、順番が回ってくると自動的に実行状態（RUN）に入ります。リクエストのアサイン後に sstat コマンドで実行開始予定時刻を確認することが出来ます。
コマンドの説明および使用例は、SX-ACE および LX406Re-2 のプログラミング利用ガイドをご確認ください。
・SX-ACE プログラミング利用ガイド
http://www.ss.cc.tohoku.ac.jp/super/guide.html#バッチリクエスト
・LX 406Re-2 プログラミング利用ガイド
http://www.ss.cc.tohoku.ac.jp/scalar/guide.html#バッチリクエスト（共同利用支援係，共同研究支援係）

— 46 —

SENAC Vol. 49, No. 4（2016. 10）

平成 29 年度の「京」を含む HPCI システム共用計算資源の利用研究課題の募集について（No.230）

平成 29 年度「京」を含む HPCI システム利用研究課題の募集(年二回募集の第一回目)と募集説明会
についてご案内いたします。
課題募集開始申請受付開始申請受付締切： 平成 28 年 9 月 7 日(水)： 平成 28 年 10 月 6 日(木)： 平成 28 年 11 月 7 日(月)17:00(JST)押印済申請書の郵送期限 平成 28 年 11 月 21 日(月)
課題選定結果通知 ： 平成 29 年 2 月初旬
募集対象資源及び課題実施期間 ：募集対象資源
募集対象課題
一般課題
若手人材育成課題「京」

課題実施期間平成 29 年 4 月 1 日～
平成 30 年 3 月 31 日(1 年)

平成 29 年 4 月 1 日～
平成 29 年 9 月 30 日(半年)※「京」以外の HPCI システム 一般課題
平成 29 年 4 月 1 日～
一般課題
産業利用課題(実証利用)
平成 30 年 3 月 31 日(1 年)※「京」以外の HPCI 共用計算資源との同時利用、第 2 希望の選択はできません。
産業利用課題(実証利用)

なお、「京」一般課題(競争的資金等獲得課題)、「京」産業利用課題(ASP 事業実証利用、個別利用(成果非公開，有償))、「京」以外の HPCI 産業利用課題(個別利用(成果非公開、有償))、ならびにトライアル・ユースについては、随時応募を受け付けております。
（スーパーコンピューティング研究部，共同利用支援係，共同研究支援係）

コンパイラのバージョンアップについて（No.231）
2016 年 9 月 26 日に C++/SX コンパイラをバージョンアップいたします。詳細につきましては、
リリースメモをご覧ください。
対象システムSX-ACE

コンパイラ名C++/SX ※1

旧バージョン
Rev.102

新バージョン
Rev.110

リリースメモ*C++/SX リリース
メモ.pdf＊リリースメモは以下からご参照ください。
http://www.ss.cc.tohoku.ac.jp/tayori/tayori/231.html※1 適用される言語仕様の既定値の変更により、従来コンパイルできたソースプログラムがコンパ
イルできなくなることがあります。詳しくはリリースメモをご確認ください。
なお、コンパイルコマンドに変更はありません。オプションの詳細は、sxman コマンドや PDF 版マニュアルで参照できます。参照方法は、以下をご覧ください。
http://www.ss.cc.tohoku.ac.jp/super/online_manual.html（共同利用支援係，共同研究支援係）

— 47 —

SENAC Vol. 49, No. 4（2016. 10）― SENAC 執筆要項 ―１．お寄せいただきたい投稿内容
サイバーサイエンスセンターでは、研究者・技術者・学生等の方々からの原稿を募集しております。
以下の内容で募集しておりますので、皆さまのご投稿をお待ちしております。なお、一般投稿いただ
いた方には、謝礼として負担金の一部を免除いたします。
・一般利用者の方々が関心をもたれる事項に関する論説・センターの計算機を利用して行った研究論文の概要
・プログラミングの実例と解説・センターに対する意見、要望
・利用者相互の情報交換

２．執筆にあたってご注意いただく事項(1)原稿は横書きです。
(2)術語以外は、「常用漢字」を用い、かなは「現代かなづかい」を用いるものとします。
(3)学術あるいは技術に関する原稿の場合、200 字～400 字程度のアブストラクトをつけてください。
(4)参考文献は通し番号を付し末尾に一括記載し、本文中の該当箇所に引用番号を記入ください。
・雑誌：著者,タイトル,雑誌名,巻,号,ページ,発行年
・書籍：著者,書名,ページ,発行所,発行年
３．原稿の提出方法
原稿のファイル形式は Word を標準としますが、PDF での提出も可能です。サイズ*は以下を参照してください。ファイルは電子メールで提出してください。
－Word の場合－
・用紙サイズ：A4
・余白：上＝30mm 下＝25mm 左右＝25mm 綴じ代＝0
・標準の文字数（45 文字 47 行）＜文字サイズ等の目安＞
・表題＝ゴシック体 14pt 中央 ・副題＝明朝体 12pt 中央
・氏名＝明朝体 10.5pt 中央
・所属＝明朝体 10.5pt 中央
・本文＝明朝体 10.5pt
・章・見出し番号＝ゴシック体 11pt～12pt*余白サイズ、文字数、文字サイズは目安とお考えください。
４．その他(1)執筆者には、希望により本誌*（10 部以内の希望部数）と本誌 PDF 版を進呈します。
*2014 年末で、別刷の進呈は終了しました。
(2)一般投稿を頂いた方には謝礼として、負担金の一部を免除いたします。免除額は概ね 1 ページ
1 万円を目安とします。詳細は共同利用支援係までお問い合わせください。
(3)投稿予定の原稿が 15 ページを超す場合は共同利用支援係まで前もってご連絡ください。
(4)初回の校正は、執筆者が行って、誤植の防止をはかるものとします。
(5)原稿の提出先は次のとおりです。
東北大学サイバーサイエンスセンター内 情報部情報基盤課共同利用支援係
e-mail uketuke@cc.tohoku.ac.jp
TEL022-795-3406

— 48 —

SENAC Vol. 49, No. 4（2016. 10）

スタッフ便り
国際学術系無線 LAN ローミング基盤「eduroam(エデュローム)」の国内基盤である「eduroam JP」
が、2016 年夏をもって 10 周年を迎えました。eduroam JP は、平成 18 年度からの文部科学省特別教育研究経費事業「大学間連携のための全国共同電子認証基盤構築事業」のプロジェクトの一つとして始まったもので、東北大学が国内運用の責任校となり、サイバーサイエンスセンター・CSI
研究室がその運用や海外連携の中心的役割を果たしてきました。これまで実証実験の位置付けだったものを、本年 4 月にようやく、国立情報学研究所(NII)の正式サービス化まで漕ぎ着けることができました。学術認証フェデレーション(学認)と並ぶ、国内の新しい学術基盤の一つにまで育ったと自負しておりますが、ここまで実用化が進んだのは、先進環境を積極的に取り入れた約 160の参加機関や、協力会社のお陰です。
世界の eduroam は、単に普及促進のみならず、新技術の取り入れや新システムの開発、適用先の拡充など、現在でも進化が続いています。eduroam の仕組みに倣った、政府関係者向けの govroamというシステムも立ち上がりました(国内では未導入)。欧米では初等・中等教育機関にも導入例が増えており、先進 ICT 活用教育環境の整備のために、日本も早期追従すべきでしょう。国内では最近、学校無線 LAN から不正アクセスを許し、個人情報の漏洩につながった重大インシデントが発生しており、無線 LAN システムのセキュリティ対策も急務です。これらの課題に対処するため、私は今夏、小中高校向けの先進学校無線 LAN システムの啓発・普及を目的としたプロジェク
トを企画しました。詳細はこれから詰めますが、多くの学校やベンダ、ISP に協力いただけるよ
う、また、日本の先進 ICT 活用教育に貢献できるよう、頑張りたいと思います。（後藤英昭）
サイバーサイエンスセンターの向かい側に広がる新キャンパスで建設工事が進んでいます。すでに引っ越しが終わった新棟の 1 つに環境科学研究科があります。昼休みに散歩がてら気になっていた 1 階のサロンを覘いてみました。大学の施設とばかりと思っていましたが、実は仙台市で運営している「せんだい環境学習館 たまきさんサロン」というものでした。一般市民に開放している体験型スペースとなっており、各種講座も定期的に開催されているとのことでした。真新
しい椅子やテーブルでゆったりと本や DVD などを閲覧することもできます。貸出しも可能とのことで、さっそくカードを作って本を借りてきました。昼休みの楽しみが 1 つ増えました。
農学研究科の大きな新棟も大分完成してきたように見えます。他のキャンパスから引っ越して来られた方々に、当センターを以前より身近に感じていただき、新たにスーパーコンピュータ等を利用される方が増えるとうれしく思います。（K.S）

サイバーサイエンスセンター・情報部情報基盤課スタッフ
退職者のお知らせ[退職]
2016. 8.31

本郷 由美（総務係・事務補佐員）

SENAC 編集部会
小林広明 水木敬明 後藤英昭
江川隆輔 佐藤恵美子 吉田貴子
大泉健治 小野 敏 斉藤くみ子

編集・発行

印

開発中の青葉山キャンパス

刷

平成 28 年 10 月発行
東北大学
サイバーサイエンスセンター
仙台市青葉区荒巻字青葉 6-3
郵便番号 980-8578
東北大学生活協同組合
プリントコープࢩࢫࢸ୍࣒ぴィ⟬ᶵࢩࢫࢸ࣒ᶵ ✀ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢱ

6;$&(୪ิࢥࣥࣆ࣮ࣗࢱ

/;5H

ࢧ࣮ࣂ࡜࣍ࢫࢺྡ
ࣇࣟࣥࢺ࢚ࣥࢻࢧ࣮ࣂIURQWFFWRKRNXDFMS

66+ ࢔ࢡࢭࢫㄆド㘽⏕ᡂࢧ࣮ࣂ

NH\FFWRKRNXDFMSࢧ࣮ࣅࢫ᫬㛫฼⏝ࢩࢫࢸ࣒ྡ➼฼⏝᫬㛫ᖏ

ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢱ

㐃

⥆

㐠

㌿

୪ิࢥࣥࣆ࣮ࣗࢱ

㐃

⥆

㐠

㌿

ࢧ ࣮ ࣂ

㐃

⥆

㐠

㌿

ྍど໬ᶵჾᐊ

ᖹ᪥ 㹼

㤋ෆ฼⏝ᖹ᪥ 㹼

ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢱ㸦6;$&(㸧ࡢ฼⏝ᙧែ࡜ไ㝈್฼⏝ᙧែ฼⏝ࣀ࣮ࢻᩘ ᐇ⾜᫬㛫ไ㝈ͤ㸦⤒㐣᫬㛫㸧
ࠥ

つᐃ್㸸 㐌㛫᭱኱್㸸 ࣨ᭶ࠥ

つᐃ್㸸 ࣨ᭶᭱኱್㸸 ࣨ᭶ ᫬㛫

ࠥ ᫬㛫

ࠥ ᫬㛫

㏻ᖖ

↓ᩱ
ࢹࣂࢵࢢ

࣓ࣔࣜࢧ࢖ࢬไ㝈 T ࢜ࣉࢩࣙࣥ E ࢜ࣉࢩࣙࣥ

*%ࣀ࣮ࢻᩘ

฼⏝ࣀ࣮ࢻᩘ

V[*%
*%ࣀ࣮ࢻᩘI
GHEXJ

฼⏝ࣀ࣮ࢻᩘ

ͤࣀ࣮ࢻ௨ୖࢆ฼⏝ࡋࡓ୪ิᐇ⾜࡟ࡣ03,ࡢ฼⏝ࡀᚲ⏝୪ิࢥࣥࣆ࣮ࣗࢱ㸦/;5H㸧ࡢ฼⏝ᙧែ࡜ไ㝈್฼⏝ᙧែ฼⏝ࣀ࣮ࢻᩘͤᐇ⾜᫬㛫ไ㝈㸦⤒㐣᫬㛫㸧࣓ࣔࣜࢧ࢖ࢬไ㝈

㏻ᖖ

ࠥ

つᐃ್㸸 ࣨ᭶᭱኱್㸸 ࣨ᭶*%ࣀ࣮ࢻᩘ

T ࢜ࣉࢩࣙࣥ E ࢜ࣉࢩࣙࣥ
฼⏝ࣀ࣮ࢻᩘ
O[࢔ࣉࣜࢣ࣮ࢩࣙࣥ

࡞ࡋ

*%఍ヰᆺ㸦 ࢥ࢔ࡲ࡛㸧 ᫬㛫㸦&38 ᫬㛫ྜィ㸧

*%

ͤࣀ࣮ࢻ௨ୖࢆ฼⏝ࡋࡓ୪ิᐇ⾜࡟ࡣ03,ࡢ฼⏝ࡀᚲ⏝D