኱つᶍ⛉Ꮫィ⟬ࢩࢫࢸ࣒㛵㐃᱌ෆ㸺኱つᶍ⛉Ꮫィ⟬ࢩࢫࢸ࣒㛵㐃ᴗົࡣࠊࢧ࢖ࣂ࣮ࢧ࢖࢚ࣥࢫࢭࣥࢱ࣮ᮏ㤋ෆࡢ᝟ሗ㒊᝟ሗᇶ┙ㄢࡀᢸᙜࡋ࡚࠸ࡲࡍࠋ㸼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 䛛䜙䛸䛳䛶࿨ྡ䛥䜜䛯䚹㻌SENAC Vol. 48, No. 4（2015. 10）— 1 —［共同研究成果］

結晶性分子ジャイロスコープの構造と回転動力学の密度汎関数強束縛法によるシミュレーション
菅野 学,1 小林倫仁,1 Wilfredo Credo Chung,1,2 山崎 馨,1,3 瀬高 渉,4 河野裕彦 1
1
東北大学大学院理学研究科化学専攻
2
Department of Chemistry, De La Salle University−Manila, Philippines
3
北海道大学大学院理学研究院化学部門
4
首都大学東京大学院都市環境科学研究科分子応用化学域

動的機能分子の 1 つとして有望視される分子ジャイロスコープは外部骨格（固定子）によって保護された環（回転子）を持ち, それが両者の間の化学結合を軸として回転する. 近年合成されたフッ素置換回転子を有する結晶性分子ジャイロスコープ（RotF）は, テラヘルツ波による内部回転の直接励起が可能と考えられている. その実現に向けて, 我々は密度汎関数法の 1/10 から 1/100の計算負荷で同程度の計算精度が期待できる密度汎関数強束縛（DFTB）法を用いて, RotF の結晶構造と回転動力学のシミュレーションを行った. まず, DFTB 法の計算プログラム DFTB+を並列コンピュータ NEC LX 406Re-2 に実装し, その並列化性能を検証した. 次に, 内部回転の有効ポテ
ンシャルエネルギー曲線の探索と, 熱運動の分子動力学シミュレーションを実行した. 結果の詳細な解析から, RotF の活性化エネルギーの値や室温における回転頻度を予測した.

1. 序論
新規な動的機能を発現する分子の創製およびその制御はナノテクノロジーの中心的課題である.
ナノサイズの分子は単分子でも多数の自由度を有するが, それらがファンデルワールス力で結び付いた固体（分子性結晶）においては隣接する分子が互いに相関し, 超多次元系を形成する. その超多次元空間の中から如何にして有効な反応座標を抽出し, ポテンシャルエネルギー曲面の鞍点に相当する遷移状態を制御するかという問題は, 分子の動的機能を高効率で操作するための方法論と直結する重要な学術的テーマである.
21 世紀に入り, 米国の Garcia-Garibay らは外部骨格（固定子）によって保護された環（回転子）が結合を軸として回転する分子ジャイロスコープの概念を提案し, 実際に種々の結晶性分子ジャイロスコー
固定子プを合成した. 図 1 にその代表例と核磁気共鳴（Nuclear magnetic resonance; NMR）スペクトル測定
回転子から見積もられた回転の活性化エネルギーを示した
回転軸[1-3]. 嵩高い固定子ほど結晶中で回転子の周囲に大きな空間を確保し, 活性化エネルギーを低下（遷移状態を安定化）させる.
瀬高らは固定子の設計によって遷移状態の安定性を調節するこのアイディアを発展させ, 大きく頑強
12.8 kcal/mol< 5 kcal/mol
なジシロキサンかご状骨格を持つ分子ジャイロスコ
ープ（図 2 で X = H としたもの, 以下 RotH と表記） 図 1 Garcia-Garibay らが合成した結晶性分子ジャイロスコープの例. 数値はを合成した[4]. 結合軸の両端を架橋したジシロキサ
実測された回転の活性化エネルギー
ン鎖が周囲の影響を軽減し, フェニレン環は室温で[1-3].
NMR の時間スケールよりはるかに速く回転する.

— 2 —

SENAC Vol. 48, No. 4（2015. 10）

河野らは RotH の結晶構造を理論計算で再現することに成功し, 活性化エネルギーが図 1 の従来型と比べて非常に低い（0.7 kcal/mol）ことを示した[5].
瀬高らが合成した分子ジャイロスコープの複屈折性が温度に依存して劇的に変化する様子が偏光顕微鏡観察により確認されている[6]. このことから, 環の内部回転を瞬時に誘起できれば, 液晶よりも高速に動作する光学材料としての応用が見込まれる. しかし, 温度だけでは回転子の高速駆動は不可能である. そこで我々は, 固定子の設計による遷移状態の安定化に加えて, 電磁波などの外部摂動によっ
回転子て回転子を直接駆動すれば, 有効反応座標に沿った内部回
固定子
転運動を誘起して遷移状態を高速に通過させることができるとの着想に至った. この方針の下, 瀬高らはフェニレン
X = H or F
si = SiMe2環の片側をフッ素で置換した分子ジャイロスコープ（図 2
図 2 ジシロキサンかご状骨格をで X = F としたもの, 以下 RotF と表記）の合成に成功した持つ結晶性分子ジャイロスコー[7]. フッ素置換によって回転子−固定子間の立体反発が大プ. 無置換体（X = H）[4]を RotH,
きくなり, RotH と比べて回転の活性化エネルギーは増加す
フッ素置換体（X = F）[7]を RotFると予測される. しかし, それと引き替えに RotF の回転子と呼ぶ.には大きな電気双極子モーメントが発生するため, 外部電場（特にテラヘルツ波）による回転駆動が原理的に可能である.
テラヘルツ波を利用した高速内部回転駆動の将来的な実現に向けて, まずは RotF の活性化エネ
ルギーの値や回転の動力学的機構を理論計算で解明する必要がある. ナノサイズの分子が集合した結晶性分子ジャイロスコープのシミュレーションを実行するためには, 以下に挙げる 3 つの条件を満たす計算法が必要である.
(ⅰ) 単位格子を無限に繰り返す周期境界条件の下での計算が行える.
(ⅱ) 計算負荷が低く, ピコ秒からナノ秒に及ぶ分子動力学シミュレーションが実行できる.
(ⅲ) 回転子−固定子間および分子間に働くファンデルワールス力までも高い精度で評価できる.
我々はこれら全ての条件を満たす密度汎関数強束縛（Density-functional-based tight-binding; DFTB）
法[8-11]を採用した. 固体物理学の分野で開発された DFTB 法は, 密度汎関数（Density functional
theory; DFT）法に基づく半経験的手法である. 上述のように, 我々は既に DFTB 法を RotH のシミュレーションに適用してその有用性を確認している[5]. 本研究では, DFTB 法の代表的な計算プ
ログラム DFTB+[12]を東北大学サイバーサイエンスセンター所有の並列コンピュータシステム
NEC LX 406Re-2 に実装した. その並列化性能をテストした後に, RotF の結晶構造と回転動力学のシミュレーションを行った.

2. 密度汎関数強束縛（DFTB）法
最初に, DFTB 法の理論的枠組みを概説する. DFTB 法では, 系の電子密度 ρ を中性原子で構成される電荷の偏りが無い参照系の電子密度 ρ0 と電子密度揺らぎ δρ の和で表す. そして, ρ の汎関数である DFT 法の Kohn-Sham エネルギーを ρ0 の近傍で δρ の多項式として展開する. 1 次の項は 0となり, 系のエネルギーは参照系の Kohn-Sham ハミルトニアンと原子間反発から得られる 0 次の主要項, および δρ に依存する 2 次以上の補正項に分けられる.
電荷の偏りを無視して 0 次解を評価する手法を Non-self-consistent-charge (Non-SCC) DFTB 法[8,9]と呼ぶ. 強束縛近似（化学の用語では LCAO 近似）に基づいて参照系の Kohn-Sham 軌道を原子軌道の線形結合で表すと, Kohn-Sham 方程式は次の強束縛（永年）方程式に帰着される.

det  H 0   S   0(1)

結晶性分子ジャイロスコープの構造と回転動力学の密度汎関数強束縛法によるシミュレーション

— 3 —

0
H 0   H 
 は原子軌道間の Kohn-Sham ハミルトニアン行列, S   S  は重なり積分（内積）行列である. DFTB 法は, これらの行列要素を事前に DFT 法で見積もられたパラメーターとすること
で, 高速かつ DFT 法に近い精度で結晶の計算を可能にする. (1)式を解いて Kohn-Sham 軌道の軌道
エネルギー  i  と展開係数 ci  が得られ, 系のエネルギーE はE  2 f i  i  Erep

(2)iと与えられる. 2fi は Kohn-Sham 軌道 i の平均占有数であり, 通常 fi   0,1 は Fermi-Dirac 分布関数から決定される. Erep は原子間反発エネルギーであり, 解析的 2 体ポテンシャルの和として記述され, DFT 計算に合うようにフィッティングされる. 電荷の偏りを無視するこの Non-SCC DFTB 法
はナノカーボンなどの等核分子や非極性分子の計算に適している.
一方, 系のエネルギーに対して 2 次の補正項も考慮する手法を Self-consistent-charge (SCC)
DFTB 法[10]と呼ぶ. SCC DFTB 法では, (1)式の代わりに次の行列式を解く.

det  H   S   0(3)

H   H   は参照系の Kohn-Sham ハミルトニアン行列 H0 に各原子の電子密度揺らぎの間に働く
Coulomb 相互作用を加えたものである.

1
0H   H  S    IK   JK  qK
2
K

(4)

添字 I と J はそれぞれ原子軌道 μ と ν が局在する原子を指す. ΔqK は原子 K の Mulliken 電荷[13]の揺らぎ（中性原子との差）, γIK や γJK は原子の DFT 計算から決まる Hubbard パラメーター[14]と原子間距離に依存する関数である. (3)式と(4)式を自己無撞着的に解くためには, 反復計算が必要である. まず{ΔqK}の初期値を設定して(4)式の H を評価し, それを(3)式に代入して解く. 得られた
展開係数 ci から新しい{ΔqK}を算出し, その値が収束するまで(3)式と(4)式を解く過程を繰り返す. 最終的な軌道エネルギー  i  を(2)式に代入して系のエネルギーE が求められる. SCC DFTB 法
は主に反復計算の数だけ Non-SCC DFTB 法よりも計算負荷が高く, 典型的には 10 倍以上に達する.
それでもパラメーターを使用しない DFT 法と比べればはるかに低負荷であり, 異核分子や極性分
子の電荷分布を適切に記述することが期待できる.

 

3. MPI 並列版 DFTB+の並列化性能
DFTB 法に基づく量子化学計算パッケージの 1 つである DFTB+[12]は, ナノ材料などの大規模系への応用に向けた並列化が進められている. 本章では, 最近テスト版が公開された MPI 並列版
DFTB+の並列化性能を検証するために SiC 単結晶と RotF を対象として行った計算の結果を報告する. その準備として, 次節では RotF の結晶構造について説明する.

3.1

RotF の結晶構造

瀬高らが X 線結晶構造解析を用いて温度 273 K で観測した RotF の単位格子を図 3(a)に描画した. これは 4 分子で構成される直方晶（斜方晶）であり, 互いに直交する 3 つの格子ベクトルの長さは図中に記した通りである. 空間群 Pbcn に属し, 格子内のいずれか 1 つの分子に反転, 21 らせん, 対角映進の 3 種類の対称操作を施すことで他の 3 つの分子が得られる. 図 3(b)は格子内の 1分子を抜き出したもの, 図 3(c)はそれを固定子と回転子をつなぐ Si−C 結合（回転軸）に沿って眺めたものである. 一見すると, フェニレン環に 4 個のフッ素が置換され, 分子全体として C2 対称性を持っているように見える. RotF には図 3(d)に示した互いに完全に等価な安定構造 A と B が存

— 4 —

SENAC Vol. 48, No. 4（2015. 10）

在し, 結晶中の各分子が 1:1 の存在比でランダムに A または B の構造にある. X 線結晶構造解析では, 平均的な構造として 2 つの安定構造が重なって観測されている. 1 分子に含まれる原子数は195 個, 単位格子の 4 分子では 780 個に達する.
固定子を成す 3 本の鎖に特徴的な Si−O−Si（ジシロキサン）結合はほぼ直線的である. 図 3(c)の左右の鎖は捻れて広がっているのに対し, 中央の鎖は C2 軸方向へ伸びている. 一方, 回転子に着目すると, 環とフッ素の C−F 結合は環平面に対してわずかに C2 軸方向へ傾く.

図 3 X 線結晶構造解析から得られた温度 273 K における RotF の結晶構造[7]. (a) 4 分子で構成される単位格子. (b) 格子内の 1 分子. (c) 回転軸に沿って眺めた 1 分子. (d) 等価な 2 つの安定構造 Aと B の回転子近傍の拡大図. 観測された結晶構造は安定構造 A と B が重なった平均的な構造である.

3.2

SiC 単結晶と RotF を対象としたテスト計算

並列コンピュータシステム NEC LX 406Re-2 に MPI 並列版 DFTB+を実装した. コンパイルにはIntel Fortran Compiler XE 14.0 を使った. また, MPI ライブラリとして Intel MPI Library 4.1 Update 3
Build 048 を用い, MPI 並列用の線形代数ライブラリには Intel MKL Version 11.1 Update 1 に収録されている ScaLAPACK を指定した.
まず, SiC 単結晶をテストの対象とした. 原子数 1,000（原子価軌道数 No = 4,000）, 2,000（No =8,000）, 4,000（No = 16,000）の 3 つの系に周期境界条件を課して, SCC DFTB 法で電子エネルギー計算を行った. 計算終了までの経過時間（Wall clock time）とコア数の両対数グラフを図 4 に示した. コア数の増加と共に経過時間が期待通りにほぼ線形減少することが分かる.
次に, RotF に関しても経過時間のコア数依存性を調べた. 図 3 の結晶構造を基にして, 4 分子が

— 5 —

結晶性分子ジャイロスコープの構造と回転動力学の密度汎関数強束縛法によるシミュレーション

同じ安定構造にある単位格子を用意した. 図 5(a)は 4分子が全て安定構造 A にある場合, 図 5(b)は全て安定構造 B にある場合の単位格子を明示したものである.
両者は等価な結晶を成すので, どちらを計算に用いても同じ結果が得られる. 2×2×2 = 8 個の単位格子の集まりに対して周期境界条件を適用した. 総原子数は6,240, 総軌道数は 13,632 である. SiC 単結晶の 16,000軌道の場合と同様に, コア数の増加につれて経過時間が減少した（図 4 の○印）.
以上のように, 16,000 軌道程度の計算では, 384 コア（16 ノード）まで線形スケールすることが確認された.
しかし, 384 コアを使用した場合の並列化効率はわず
か 19 %であった（48 コアの場合は 50 %）. この原因を究明するために, 大規模計算パフォーマンス解析ツール Score-P[15]および Scalasca[16]を用いてボトルネック解析を行った. その結果, ScaLAPACK 収録の分割
統治法（Divide-and-conquer algorithm）に基づくハミル
トニアン行列の対角化サブルーチンにおいてデータの通信・同期を行う際に遅延が発生し, 計算負荷が使用されているコアに対して均等に分散されていないことが判明した. この通信･同期の遅延は, 子プロセス間のデータ送受信の順番が効率的に行われていないことが原因と考えられる. 現在, 子プロセス間のデータ送受信等を中心としたプログラムコードの最適化を進めている.

図 4 NEC LX 406Re-2 において測定した, 周期境界条件を課した SiC 単結晶の電子エネルギー計算の経過時間（Wallclock time）とコア数の関係 . 原子数1,000（原子価軌道数 No = 4,000）, 2,000（No = 8,000）, 4,000（No = 16,000）の結果をそれぞれ△, □, ◇で示した. ○はRotF の 2×2×2 = 8 個の単位格子の集まりに周期境界条件を適用した例. 実線
はプロセス数を n 倍にしたときに経過
時間が 1/n になる理想値を表している.

図 5 4 分子が同じ安定構造にある単位格子. (a) 全て安定構造 A の場合. (b) 全て安定構造 B の場
合. 回転子の向きが分かり易いように, 各分子のフッ素を点線で囲んだ. (a)を単位格子とする結晶と(b)を単位格子とする結晶は等価である.

4. RotF の結晶構造と回転動力学のシミュレーション
フッ素置換の影響で活性化エネルギーがどの程度増加するかを評価するために, RotF の有効反応座標に沿ったポテンシャルエネルギー曲線を求めた. また, 熱運動で環が回転して安定構造間— 6 —

SENAC Vol. 48, No. 4（2015. 10）を移動する動的挙動を探る分子動力学シミュレーションを実行した. 以下では, 比較的少数のコ
アで済む（並列化効率が高い）単位格子 1 個に周期境界条件を課した計算の結果を報告する.

4.1

結晶構造最適化

図 5 の 4 分子が同じ安定構造（A または B）にある単位格子 1 個に周期境界条件を課し, 格子定数は実験値に固定して SCC DFTB 法で構造最適化を行った. 図 3 の安定構造 A と B に対応す
る等価な 2 つの SCC DFTB 最適構造を図 6 に示した. 図 3 と比較すると, 固定子は C2 対称性を失い, Si−O−Si 結合が折れ曲がっている. 回転子では C2 軸に垂直だった環平面が傾いて C−F 結合とほぼ平行になった.

図 6 4 分子が同じ安定構造にある単位格子を SCC DFTB 法で最適化した結果. 格子内の 1 分子を抜き出した. (a)は安定構造 A, (b)は安定構造 B に対応しており, 両者は等価である. 紫色の 4 原子で回転軸周りの二面角を定義した.
回転子の配向を定量化する尺度として, 図 6 の紫色で指定した 4 原子で回転軸周りの二面角を定義し, その値を表 1 にまとめた. 構造 A と B は等価であるが, 前者の二面角は実験値に近いのに対し, 後者はややその差が大
きい. 鎖の変形と環平面の傾きの影響が, A では相殺す
る一方, B では共に二面角を小さくするように働いているためである. 細部に関しては精度改善の余地があるが,分子全体の概形は図 3 と良く一致している.

4.2

表 1 安定構造における回転軸周りの二面角（deg.）.
A

B

実験

−76.9

101.0

SCC DFTB

−66.3

73.0

内部回転の有効ポテンシャルエネルギー曲線

上述のように, RotF の単位格子は 780 個の原子を内包する超多次元系である. 回転の活性化エ
ネルギーを見積もるためには, その中から有効な反応座標を抽出する必要がある. そこで, 4 分子が同じ安定構造にある単位格子を起点として, 任意の 1 分子の二面角を変数として 1°ずつ変化さ
せ, それぞれの角度において他の自由度を最適化した. 360 回の構造最適化には時間を要するので,この計算には Non-SCC DFTB 法を用いた（極小点と遷移状態のみ SCC DFTB 法でも求めた）.

— 7 —

結晶性分子ジャイロスコープの構造と回転動力学の密度汎関数強束縛法によるシミュレーション

Potential energy (kcal/ml)

図 7 は二面角の関数として得られた内部回転の有効ポ
BA
テンシャルエネルギー曲線である. 80°の極小点では 4 分子が全て安定構造 B にあり, −65°の極小点は 1 分子のみ安定構造 A に移ったことに対応する. 環が伸びた鎖に重なる 0°付近と±180°付近に遷移状態が存在し, 活性化エネル
ギーの値はそれぞれ 5.6 kcal/mol と 4.5 kcal/mol（SCC DFTB
法では 7.0 kcal/mol と 5.8 kcal/mol）である。フッ素で置換されたために RotH の値（0.7 kcal/mol）と比べて増加したが, それでも図 1 のかご状骨格を持たない分子ジャイロス
コープの中で最も活性化エネルギーが低いものと同程度
Dihedral angle (deg.)である.
図 7 内部回転の有効ポテンシャ
回転に伴う固定子の構造変化を詳しく調べたところ,
ルエネルギー曲線. 回転を起こす
遷移状態の周辺では環との立体反発を避けるように鎖が分子が 2 つの極小点で成す構造を変形する. 特に, フッ素と鎖が接近する 0°付近では大きく明記した. 他の分子の構造は B.
歪む. RotF の有効反応座標は環が結合軸周りに回転するだけの単純なものではなく, 回転子と固定子が互いに相関した複雑な運動である. 一方, 2 つの極小
点のエネルギーがほぼ等しい（SCC DFTB 法でも差は 0.4 kcal/mol）ことから, 分子間の双極子−
双極子相互作用は小さいと示唆される. これは X 線結晶構造解析において回転子の配向秩序が現れない（ランダムに A または B の構造にある）事実と合致する. しかし, 反応座標の解析からは, 回転に際して固定子の変形が隣接する分子の構造にわずかながら影響することが分かった. 図 7 のポテンシャルエネルギー曲線が回転方向に対して非対称であることも結晶中の分子間相関を反映
している.

4.3

熱運動の分子動力学シミュレーション

有限温度における熱運動が誘起する内部回転の挙動を解明するために, 全エネルギーが一定（小正準集団）の条件下で Non-SCC DFTB 法を用いた古典トラジェクトリ計算を実行した. 始時
刻の単位格子は 4 分子が全て安定構造 B にあると設定した. 運動温度の時間平均として系の温度を定義し, 300 K から 600 K の範囲を調査した.
300 K の場合, ピコ秒の時間領域では回転に至らなかった. 一方, 400 K では環がどちらの遷移状態も越えて回転することが確認された. 図 8(a)は 400 K の典型的なトラジェクトリにおける回転
軸周りの二面角の時間変化をプロットしたものである. 環が 200 ps の間に 2 つの安定構造間を移動しながら回転することが読み取れる. また, 古典トラジェクトリ計算においても, 環が遷移状態を越える瞬間に固定子が変形して回転を促す様子が見られた.

図 8 古典トラジェクトリ計算の結果. (a) 温度 400 K の典型的なトラジェクトリにおける回転軸
周りの二面角の時間変化. 有効ポテンシャルエネルギー曲線の極小点の位置を縦線で示した. (b)±180°付近の低い遷移状態を越える移動に関するアレニウスプロット.

— 8 —

SENAC Vol. 48, No. 4（2015. 10）環が安定構造に滞在する時間と回転に要する所要時間を比較すると, 前者の方が十分に長いとみなせる. そこで, 安定構造間の移動の速度定数 k を平均滞在時間の逆数として定義し, アレニウ
スプロットを作成した. 図 8(b)は±180°付近の低い遷移状態を越える移動に関するアレニウスプロ
ットである. 直線でフィッティングして得られた活性化エネルギーは約 4 kcal/mol であり, 図 7 の
有効ポテンシャルエネルギー曲線から算出した値に近い. この直線を外挿することにより, 室温（298 K）における速度定数を k = 4×108 s−1 程度と予測した. これは約 3 ns に 1 回の頻度で低い遷移状態を越えて 2 つの安定構造間を移動することを意味する. 実験では, RotF の内部回転は室温
で NMR の時間スケールよりはるかに速いことが分かっており, 回転子−固定子間の立体反発が著しく大きい塩素置換体 RotCl では 83 μs に 1 回の頻度と見積もられている[7]. 我々の計算はこれらの実験結果と矛盾しないものである.

5. まとめと展望
次世代の光機能性材料として有望視される結晶性分子ジャイロスコープ RotF の回転動力学を解明するために, DFTB 法によるシミュレーションを行った. 計算結果の詳細な解析から得られた成果を以下にまとめる.
(ⅰ) 周期境界条件を課した DFTB 法は RotF の結晶構造を半定量的に再現できた.
(ⅱ) 内部回転の有効ポテンシャルエネルギー曲線の探索から, 2 つの遷移状態（活性化エネルギ
ーは 5.6 kcal/mol と 4.5 kcal/mol）の存在が明らかになった.
(ⅲ) 環の回転に際しては固定子も相関して動く.
(ⅳ) アレニウスプロットから室温での回転頻度を算出し, 約 3 ns に 1 回の頻度で低い遷移状態を超えると予測した.
今後は, テラヘルツ波照射による RotF の回転駆動シミュレーションを行う計画である. RotH と比較することで, テラヘルツ波から回転子へのエネルギー供給に対するフッ素置換の効果を調べていく. 分子ジャイロスコープは分子機械としても位置づけることができ[17]，今後のさらなる展開が期待されている. このような大型計算機を活用した実在系のシミュレーション研究から, 高い機能性を有する分子の設計に繋げていきたい.

謝辞
本研究の計算の一部は, 東北大学サイバーサイエンスセンターの並列コンピュータシステム
NEC LX 406Re-2 を利用することで実現できた. なお, 本研究の一部は JSPS 科研費 15KT0138 および文部科学省 HPCI 戦略プログラム分野 2「新物質・エネルギー創成」の助成を受けて行われた.

参考文献[1] Z. Dominguez, H. Dang, M. J. Strouse, and M. A. Garcia-Garibay, Molecular “Compasses” and“Gyroscopes.” III. Dynamics of a Phenylene Rotor and Clathrated Benzene in a Slipping-Gear Crystal
Lattice, J. Am. Chem. Soc. 124, 7719 (2002).
[2] T.-A. V. Khuong, G. Zepeda, R. Ruiz, S. I. Khan, and M. A. Garcia-Garibay, Molecular Compasses
and Gyroscopes: Engineering Molecular Crystals with Fast Internal Rotation, Cryst. Growth Des. 4,
15 (2004).
[3] Z. J. O’Brien, A. Natarajan, S. I. Khan, and M. A. Garcia-Garibay, Synthesis and Solid-State
Rotational Dynamics of Molecular Gyroscopes with a Robust and Low Density Structure Built with a
Phenylene Rotator and a Tri(meta-terphenyl)methyl Stator, Cryst. Growth Des. 11, 2654 (2010).
[4] W. Setaka, S. Ohmizu, C. Kabuto, and M. Kira, A Molecular Gyroscope Having Phenylene Rotator
Encased in Three-spoke Silicon-based Stator, Chem. Lett. 36, 1076 (2007).

結晶性分子ジャイロスコープの構造と回転動力学の密度汎関数強束縛法によるシミュレーション— 9 —[5] A. B. Marahatta, M. Kanno, K. Hoki, W. Setaka, S. Irle, and H. Kono, Theretical Investigation of theStructures and Dynamics of Crystalline Molecular Gyroscopes, J. Phys. Chem. C 116, 24845 (2012).
[6] W. Setaka and K. Yamaguchi, Thermal modulation of birefringence observed in a crystalline
molecular gyrotop, Proc. Natl. Acad. Sci. USA, 109, 9271 (2012).
[7] W. Setaka, S. Ohmizu, and M. Kira, Molecular Gyroscope Having a Halogen-substituted p-Phenylene
Rotator and Silaalkane Chain Stators, Chem. Lett. 39, 468 (2010).
[8] D. Porezag, Th. Frauenheim, Th. Köhler, G. Seifert, and R. Kaschner, Construction oftight-binding-like potentials on the basis of density-functional theory: Application to carbon, Phys.
Rev. B 51, 12947 (1995).
[9] G. Seifert, D. Porezag, and Th. Frauenheim, Calculations of Molecules, Clusters, and Solids with a
Simplified LCAO-DFT-LDA Scheme, Int. J. Quantum Chem. 58, 185 (1996).
[10] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, Th. Frauenheim, S. Suhai, and G. Seifert,
Self-consistent-charge density-functional tight-binding method for simulations of complex materials
properties, Phys. Rev. B 58, 7260 (1998).
[11] M. Elstner and G. Seifert, Density functional tight binding, Phil. Trans. R. Soc. A 372, 20120483(2014).
[12] B. Aradi, B. Hourahine, and Th. Frauenheim, DFTB+, a Sparse Matrix-Based Implementation of the
DFTB Method, J. Phys. Chem. A 111, 5678 (2007).
[13] R. S. Mulliken, Electronic Population Analysis on LCAO–MO Molecular Wave Functions. I, J. Chem.
Phys. 23, 1833 (1955).
[14] J. Hubbard, Electron Correlations in Narrow Energy Bands, Proc. R. Soc. A 276, 238 (1963).
[15] A. Knüpfer, C. Rössel, D. an Mey, S. Biersdorff, K. Diethelm, D. Eschweiler, M. Geimer, M. Gerndt,
D. Lorenz, A. Malony, W. E. Nagel, Y. Oleynik, P. Philippen, P. Saviankou, D. Schmidl, S. Shende, R.
Tschüter, M. Wagner, B. Wesarg, and F. Wolf, Score-P: A Joint Performance Measurement Run-Time
Infrastructure for Periscope, Scalasca, TAU, and Vampir, in Tools for High Performance Computing
2011, 79 (Springer, 2012).
[16] M. Geimer, F. Wolf, B. J. N. Wylie, E. Ábrahám, D. Becker, and B. Mohr, The Scalasca performancetoolset architecture, Concurrency Computat.: Pract. Exper. 22, 702 (2010).
[17] S. Erbas-Cakmak, D. A. Leigh, C. T. McTernan, and A. L. Nussbaumer, Artificial Molecular Machines,
Chem. Rev. 115, 10081 (2015).

— 10 —

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

グリーンランド氷床モデルを用いた地震波伝播シミュレーション

豊国源知 1・竹中博士 2・金尾政紀 3・坪井誠司 4・東野陽子 5
1 東北大学大学院 理学研究科 地震・噴火予知研究観測センター
2 岡山大学大学院 自然科学研究科 地球生命物質科学専攻
3 国立極地研究所
4 海洋研究開発機構
5 文部科学省
平成 26 年度の共同研究課題では，局地的な地震波伝播を効率良くモデリングできる数値計算手法を開発した．また本手法による計算例として，北極グリーンランド氷床を対象とした地震波
伝播シミュレーションを行った．グリーンランドでは 11 ヶ国の国際協力で，2011 年から本格的な地震観測が進められてきた．特に氷床上の観測点は，日米共同観測隊が毎年現地で行うメンテ
ナンスによって維持されており，得られたデータの解析が鋭意進められている．解析の際には，観測地震波形に与える氷床の影響を正確に考慮する必要があるが，これまで研究はほとんど行われておらず，本研究がその端緒を開くと期待される．本稿では，我々の計算手法とシミュレーシ
ョン成果の概要を紹介する．

1.

はじめに

グリーンランドは北極圏に位置する世界最大の島で，表面の 80%以上は氷床に覆われている．
氷床とは，地表を覆う全面積が 5 万 km2 以上の氷塊のことで，グリーンランドと南極にしか存在しない．グリーンランド氷床の厚さは平均で 2 km，最厚部では 3 km にも達し，すべてが融解した場合，海水準が 7 m 以上も上昇すると予測されている．近年は地球温暖化の影響により，氷床の融解が指摘されていることから，融解の状況をモニタリングする方法として，地震観測に期待が寄せられることとなった．融解に伴い，氷が割れたり滑ったりする際に生じる振動は，地震計で即時に検出できるからである．しかしグリーンランドの内陸部は，人を寄せ付けない過酷な環境にあるため，2008 年の時点では，周辺の島々を合わせても地震観測点は 16 点，うち氷床上はわずか 1 点という状況で，観測網の拡充が望まれていた．
こうした経緯で，2009 年に「グリーンランド氷床モニタリング観測網（GreenLand Ice Sheet
グ リ ス ン
monitoring Network，略称GLISN）」と呼ばれる国際プロジェクトが発足し，観測点の増強が進められてきた．現在は 11 ヶ国が参加し，氷床上 4 点を含む計 33 点の地震観測点が運用されている（図 1）．日本はプロジェクト発足時からの参加国で，2011 年以降は毎年アメリカ隊との共同観測隊を派遣して，氷床上 3 観測点と沿岸域 3 観測点の設置やメンテナンスに携わってきた．氷床上の観測点は維持に多大な困難が伴うが，良質な観測網の構築には均一な観測点分布が不可欠であり，日米共同観測隊の貢献はプロジェクト全体でも極めて大きいといえる．2014 年には衛星通信を用いた氷床上からの即時データ転送にも成功し，極地で観測された地震波形を，世界中ですぐさまダウンロードできるようになった[1-5]．
GLISN 観測網では，広帯域 3 成分地震計を用いた連続観測を行っているため，得られるデー
タの情報量は極めて多く，氷床融解の検出以外の目的にも利用することができる．例えば，グリーンランドは詳細な地下構造が明らかされておらず，GLISN データの解析による研究が待たれている．データ解析の際には，氷床が観測波形に与える影響がノイズとなることが予想されるため，解析に先立って氷床の影響を正しく見積もっておく必要がある．しかし氷床中を伝播する地震波に関する先行研究は，1960 年代に 1 篇が見られるのみで[6]，数値計算等で理論的に取り扱った研究はいまだ行われていない．

グリーンランド氷床モデルを用いた地震波伝播シミュレーション— 11 —

図 1：（左）GLISN 地震観測点の分布．
赤色は日本隊が設置・メンテナンスを行った点．黄色はその他の点．（右）日本隊が設置を行った ICESG 観測点の外観（2011 年）．

平成 26 年度の共同研究課題では，局地的な地震波伝播を効率よくシミュレーションすることができる数値計算手法を開発した．またこの手法に氷床を含むグリーンランドの地下構造モデルを入力してシミュレーションを行い，観測地震波形に氷床が与える影響を調べた．以降では本年出版した論文[7]の流れに沿いながら，手法と結果について述べていく．

2.

円筒座標系 2.5 次元モデリング

地震波伝播シミュレーションでは，ある地下構造モデルと震源モデルを与えて，地震波伝播の支配方程式を数値的に解くことで，任意の地点で観測される地震波形を理論的に予測する．モデリング手法には様々なものが存在するが，計算精度と効率の良さを併せ持っている手法は少ない．
本誌 46 巻 3 号で紹介した「軸対称 2.5 次元モデリング」は，その代表的な手法である[8]．この手法では，震源を通る鉛直軸の周りに構造の軸対称性を仮定することで，3 次元的な空間の広がりの効果を 2 次元断面に落とし込んだのち，その断面内において，面内運動 2 成分と，面外運動
1 成分の，3 次元の波動伝播を計算する．
「波動場 3 次元，構造 2 次元」という意味で，「2.5 次元」という呼称が用いられている．計算は断面上だけで行うので計算資源を節約でき，計算時間やメ
モリは 3 次元でフルに計算する場合の数千～数万分の 1 で済む．また 3 次元の波動場が得られる
ので，観測地震波形と直接比較できるというメリットを持つ．計算では最初に回転対称軸を与える必要があることから，座標系としては球座標系と円筒座標系が利用される．
我々はこれまで，地球や他の天体の全球地震波伝播をシミュレーションするため，球座標系を使ったプログラム開発を行ってきた[8-13]．しかしグリーンランドのような局地的な地震波伝播をシミュレーションする場合には，地球の中心まで含む計算を行うとメモリや計算時間が余分にかかってしまうことから，地球表層に特化したプログラム開発が必要であった．
円筒座標系��� �� ��において，�軸を鉛直下方，�軸を水平方向にとって，方位角を�で表すと，� � 0を対称軸とした��断面上で 2.5 次元計算が可能となる．��座標は，�軸が� � 0で反転することを除けば 2 次元デカルト座標系と同じ扱いができるので，地表面を� � 0とおいた矩形断面上で計算が行える．矩形断面では地表の曲率を考慮できないので，シミュレーションの対象としては，地表を平面として取り扱える，局地的（伝播距離≲ 1000 km）な地震波伝播が適している．円筒座標系を用いた 2.5 次元計算は，1960 年代から局地的な地震波伝播シミュレーションに活用されてきた[14-17]．
しかし上述の�軸が反転する問題は，従来のシミュレーションに深刻な影響を与えてきた．すな
わち，ある構造断面を��平面上で表そうとすると，� � 0の軸を挟んで，� � 0と� � �における 2

— 12 —

SENAC Vol. 48, No. 4（2015. 10）

図 2：（a）通常の円筒座標領域と（b）準円筒
座標領域[18]．星印は震源．

枚の矩形を用いなければならない．このため任意の構造モデルを߶ ൌ Ͳの面に入力すると，軸対称性によって߶ ൌ ߨの面における構造は߶ ൌ Ͳのものと対称にならざるを得ず，現実的な地下構造を取り扱うことができない（図 2a）．また対称軸である‫ ݎ‬ൌ Ͳが，剛体の壁と同じ働きをして，そこで人工的な反射波を生じさせる問題もあった．これらのデメリットは，円筒座標領域が通常（െλ ൏ ‫ ݖ‬൏ λǡ Ͳ ൑ ‫ ݎ‬൏ λǡ െߨ ൑ ߶ ൑ ߨ）という範囲を用いていたことに起因している．
これに対して文献[18]では，準円筒座標領域（െλ ൏ ‫ ݖ‬൏ λǡ െλ ൏ ‫ ݎ‬൏ λǡ െߨȀʹ ൑ ߶ ൑ ߨȀʹ）が提案された．このような範囲を用いると，構造断面は߶ ൌ Ͳにおける 1 枚の矩形のみで表され，この面には任意の不均質構造モデルが入力できる（図 2b）．座標軸の取り方を変えるだけなので，従来の軸対称モデリングのメリットは保存したまま，非対称な構造モデルが取り扱える．文献[18]は爆破地震探査を想定した研究であったことから，計算では軸対称な震源のみを扱っていた．我々はこの手法をさらに拡張し，断層を含む任意の震源メカニズムと任意形状の地形の導入も行った[7]．これにより，現在のプログラムでは極めて現実に近い地震波伝播をシミュレーションできる．
手法の詳細については，文献[19]も参照されたい．
円筒座標系を使った地震波伝播シミュレーションで基本となる方程式は，3 次元の弾性波の支配方程式である．ここで‫ݐ‬は時間，ߩは密度，ߨとߤは弾性定数，‫ݒ‬௜ ，݂௜ ，ߪ௜௝ ，‫ܯ‬௜௝ はそれぞれ粒子速度，体積力，応力テンソル，モーメントテンソルの成分である（݅ǡ ݆ ‫ א‬ሼ‫ݖ‬ǡ ‫ݎ‬ǡ ߶ሽ）．変数の上についたドットは時間微分を表す．弾性定数としてはラメ定数ߣǡ ߤを用いるのが一般的であるが，ߨ ‫ ߣ ؠ‬൅ ʹߤを導入すると，P 波の性質を単独の弾性定数で表せるメリットがある．
我々の開発したプログラムでは，この方程式に現れるすべての物理量を߶方向にフーリエ級数
展開したのち，差分法で離散化して計算を行う．フーリエ級数展開によって，震源メカニズムの߶— 13 —

グリーンランド氷床モデルを用いた地震波伝播シミュレーション

図 3：計算で用いた差分格子の模式
図．黒が粒子速度，
白が応力の格子を表す．星印は震源．

依存性が考慮でき，任意のメカニズムによる波動場の励起が可能となる．図 3 に差分法で用いる格子配置の模式図を示した．気体・液体を含む媒質でも安定に計算できるよう，粒子速度と応力の別々の格子に配置される「スタガード格子」を採用している．震源は必ず� = �の軸上に置かなければならないので，水平方向の震源位置を変えると�座標も変化する．

3.

グリーンランド氷床モデルによるシミュレーション

本節では円筒座標系 2.5 次元モデリングによるグリーンランドの地震波伝播シミュレーションの成果について述べる．
図 4 はグリーンランドの南端，深さ 10 km で発生した地震（実体波マグニチュード�� = 4.8）の上下動記録である．GLISN 観測網の露岩上 1 観測点（NRS），氷床上 2 観測点（DY2G ，ICESG）の波形を示している．氷床上の観測波形には，S 波の後続波部分に群速度 3.1‐2.6 km/s の特徴
的な振幅増幅が見られることがわかる．氷床に関連すると思われるこのような特徴的な波形の成因について調べるため，今回は表 1 と図 5 で示す 4 つのシミュレーションを行った．
図 4：2013 年 2 月 19 日にグリーンランド南端で発生した地震（深さ 10 km，�� = 4.8）の，3 観測点における上下動記
録 [7] ． 左 に 観 測 点 名 と 氷 床上・露岩上の別，震央距離（Δ）
を示した．観測点の位置は図 1
参照．すべての波形に 0.5‐833 s のバンドパスフィルタ
ー を 適 用 ． 群 速 度 3.1 ‐ 2.6
km/s の波群を太線で示す．
表 1：4 つのシミュレーションの諸元
識別名

氷床モデル

FLAT2.0-0

水平均一厚（厚さ 2.0 km）

氷床表面からの震源深さ（km） 氷床基盤からの震源深さ（km）
2.2

0.2

FLAT0.5-0

水平均一厚（厚さ 0.5 km）

0.7

0.2

TOPO-0

図 4 の A-A’断面

0.6

0.2

TOPO-5図 4 の A-A’断面

5.6

5.2

— 14 —

SENAC Vol. 48, No. 4（2015. 10）

図 5：4 つのシミュレーションで用いた構造モデルと震源位置の違い（文献[7]の図を修正）．(a)～(c)は計算に用いた構造断面．黄色が空気層，青色が氷床，赤色は岩石層．地表付近を鉛直方向に拡大して表示．星印は震源．(a) 水平成層，厚さ 2 km の氷床モデルの直下に震源を入力（FLAT2.0-0），(b) 水平成層，厚さ 0.5 km の氷床モデルの直下に震源を入力（FLAT0.5-0），(c) 現実的な地形の起伏と氷床厚分布を使用し，震源は氷床直下（TOPO-0）または氷床下 5 km
（TOPO-5）に置いた．(d) 構造断面の A-A’測線の位置を示す平面図．背景色は氷床厚．
計算のメインは，実際の地形の起伏や氷床厚分布を考慮した TOPO-0 および TOPO-5 であるが，比較対象として，氷床を厚さ一定の水平成層構造としてモデル化した FLAT2.0-0，FLAT0.5-0
でもシミュレーションを行った．氷床の影響を調べるため，氷床部分を岩石で置き換えた計算も行ったので，トータルの計算回数は 8 回である．構造モデル作成の際，氷床以外の部分は 1 次元
地球モデルである PREM[20]を用いた．TOPO-0，TOPO-5 で用いた地形や氷床厚分布は，文献，S 波速度（�� ）は，先行研究[23-25]等を参[21,22]に基づく．氷床の密度（�），P 波速度（�� ）�
考に(�, �� , �� ) = (0.914g/cm , 4.0 km/s, 2.0 km/s)を用いた．また空気層は(�, �� , �� ) = (0.00129g/
cm� , 0.0 km/s, 0.0 km/s)とした．空気層を導入すると，地震波から音波への変換波もシミュレー
ションできるが，音速は岩石中の地震波速度に比べて極端に遅く，精度良くシミュレーションするためには格子間隔を細かく切らなければならない．今回の計算は音波を対象としたものではな
いので，空気層の地震波速度をゼロとすることで，音波の伝播を抑制した[26-28]．
計算は SX-9 のジョブクラス p8 による自動並列演算で，水平（�）700 km，鉛直（�）150 kmの構造断面を，14000×3000 の空間差分格子に分割して行った．格子間隔は水平・鉛直ともに0.05 km である．時間刻みは 0.0025 s とし，励起後 150 s まで計算した．震源にはピュアな縦ず
れ断層型のメカニズムを用い，震源時間関数として幅 2 Hz のベル型パルスを入力した．図 6 は，
4 つのシミュレーションで得られた理論波形を，3 つの震央距離（Δ =50，100，150 km）について並べて表示した図である．また図 7～11 では，励起後 10，30，50，70，90 s における構造断面上の地震波伝播のスナップショットを並べている．赤色が P 波，緑色が S 波で，色が濃いほど振幅が大きい．以下ではこれらの理論波形とスナップショットを用いて，4 つのシミュレーショ
ンの結果について述べる．

グリーンランド氷床モデルを用いた地震波伝播シミュレーション— 15 —

図 6：4 つのシミュレーションによる，震央距離Δ =50，100，150 km における理論波形[7]．
（左）上下動，（右）水平動．震央距離とシミュレーションの識別名はパネルの左に表示した．
赤線は氷床あり，黒線は氷床なし．主要なフェーズ（Pg，Sg など）の到着時刻を青点線で示した．太線は図 4 に対応する群速度 3.1‐2.6 km/s の波群．

3.1 FLAT2.0-0
氷床が厚さ 2 km 一定の水平成層構造と仮定し，震源を氷床直下に置いた場合，極めて継続時
間の長い S 波後続波が現れた．氷床を入れない場合（図 6 黒線），地殻の中を通過した S 波（Sg）と表面波による単純なパルスの到着後，振幅はほぼゼロになっている．一方，氷床を入れた場合（図 6 赤線），Sg 到着後，数 10 秒以上大振幅が継続する．大振幅の継続時間は震央距離に比例し，Δ = 150 kmの地点では 80 s にも達している．
図 7 は氷床を入れない場合，図 8 は入れた場合の地震波伝播のスナップショットである．図 7と図 8 の比較から示されるように，氷床を入れた場合の継続時間の長い S 波後続波は，氷床内部に地震波のエネルギーが強くトラップされることで起きる現象である．氷床直下に与えられた震源から放出された地震波は，エネルギーの大部分が氷床に入射し，薄く低地震波速度の氷床内部を多重反射しながら伝播する．氷床と基盤岩との境界が水平な場合，全反射が卓越することで，地震波のエネルギーは氷床から漏れ出しにくく，継続時間の長い S 波後続波が形成される．
低地震波速度の層内で地震波がトラップされる現象としては，地殻内部を多重反射する S 波である「Lg 波」がよく知られている．一方，氷床内部の多重反射は，P 波の多重反射の観測事例が文献[6]で示されているものの，理論的研究は本研究が最初である．本研究では，氷床を意味するドイツ語「Eisdecke」の頭文字から，氷床内を多重反射する S 波を「Le 波」と命名した．同様に，地殻表面を伝播する表面波（レイリー波）は「Rg 波」と呼ばれることから，氷床表面を伝播する表面波を「Re 波」と命名した[7]．

— 16 —

SENAC Vol. 48, No. 4（2015. 10）

図 7：FLAT2.0-0 で氷床がない場合の地震波伝播の様子．図 5(a)の氷床モデルを直下の地殻構
造で置き換えて計算した．上から励起後 10，30，50，70，90 s のスナップショット．赤色は P
波，緑色は S 波を表す．色が濃いほど振幅が大きい．構造断面のうち，深さ 30 km 以浅の部分を拡大して表示した．震源の位置は星印．Pg，Sg は地殻中を伝播する P 波と S 波，Pn，Snはマントル中を伝播する P 波と S 波，Lg は地殻内部にトラップされた S 波，Rg は地殻表面を伝播する表面波．

グリーンランド氷床モデルを用いた地震波伝播シミュレーション— 17 —

図 8：FLAT2.0-0 による地震波伝播の様子（文献[7]の図を修正）．上から励起後 10，30，50，
70，90 s のスナップショット．赤色は P 波，緑色は S 波を表す．色が濃いほど振幅が大きい．
構造断面のうち，深さ 30 km 以浅の部分を拡大して表示した．震源の位置は星印．Pg，Sg は地殻中を伝播する P 波と S 波，Pn はマントル中を伝播する P 波，Lg は地殻内部にトラップされた S 波，Le は氷床内部にトラップされた S 波，Re は氷床表面を伝播する表面波．図 7 と比
較すると，トラップ波によって形成された地表付近の大振幅が，長い水平距離にわたって続いていることがわかる．

— 18 —

SENAC Vol. 48, No. 4（2015. 10）

図 9：FLAT0.5-0 による地震波伝播の様子（文献[7]の図を修正）．上から励起後 10，30，50，
70，90 s のスナップショット．赤色は P 波，緑色は S 波を表す．色が濃いほど振幅が大きい．
構造断面のうち，深さ 30 km 以浅の部分を拡大して表示した．震源の位置は星印．Pg，Sg は地殻中を伝播する P 波と S 波，Pn，Sn はマントル中を伝播する P 波と S 波，Lg は地殻内部にトラップされた S 波，Le は氷床内部にトラップされた S 波，Re は氷床表面を伝播する表面波．氷床が薄いため，トラップされた波が十分成長できず，図 8 に比べて細かな振動が卓越している．大振幅の部分の水平距離は，図 8 のものとほぼ同じである．

グリーンランド氷床モデルを用いた地震波伝播シミュレーション— 19 —

図 10：TOPO-0 による地震波伝播の様子（文献[7]の図を修正）．上から励起後 10，30，50，
70，90 s のスナップショット．赤色は P 波，緑色は S 波を表す．色が濃いほど振幅が大きい．
構造断面のうち，深さ 30 km 以浅の部分を拡大して表示した．震源の位置は星印．Pg，Sg は地殻中を伝播する P 波と S 波，Pn はマントル中を伝播する P 波，Lg は地殻内部にトラップされた S 波，Le は氷床内部にトラップされた S 波，Re は氷床表面を伝播する表面波．地形や氷床基盤の起伏によってトラップ波のエネルギーが漏出するため，図 8，9 と比べて地殻内部に緑色部分が卓越している．

— 20 —

SENAC Vol. 48, No. 4（2015. 10）

図 11：TOPO-5 による地震波伝播の様子（文献[7]の図を修正）．上から励起後 10，30，50，
70，90 s のスナップショット．赤色は P 波，緑色は S 波を表す．色が濃いほど振幅が大きい．
構造断面のうち，深さ 30 km 以浅の部分を拡大して表示した．震源の位置は星印．Pg，Sg は地殻中を伝播する P 波と S 波，Pn はマントル中を伝播する P 波，Lg は地殻内部にトラップされた S 波．氷床に入射するエネルギーが少なく，トラップ機構も弱いため，氷床内の振幅は図8～10 と比べて小さく，氷床がない場合の波動場（図 7）に近い．

グリーンランド氷床モデルを用いた地震波伝播シミュレーション— 21 —

3.2 FLAT0.5-0
氷床が厚さ 0.5 km 一定の水平成層構造と仮定し，震源を氷床直下に置いた場合，前節の結果と似た継続時間の長い S 波後続波が現れた．図 6 の赤線が示すように，S 波後続波の継続時間は
FLAT2.0-0 とほとんど変わらず，氷床の厚さには依存しないと考えられる．一方，周波数は前節の場合よりも高周波が卓越している．図 9 からは，波がトラップされる過程は前節の結果と同様
であるが，氷床が薄いことで波の成長が阻害され，高周波成分が卓越したことがわかる．
3.3 TOPO-0
現実的な地形・氷床厚分布を用い，震源を氷床直下に置いた場合は，水平成層の場合と大きく異なる波形が得られた．S 波後続波の大振幅は Sg の到着からかなり遅れて現れるうえ，継続時
間は短く，波群がコンパクトにまとまっている（図 6 赤線）．このケースも震源が浅いため，地震波のエネルギーの大部分が氷床内に入射する点は前節までと同じである．しかし地形や基盤岩の凹凸によって P 波から S 波，S 波から P 波への変換が起こり，エネルギーが分散されることに加え，波の全反射が崩れてエネルギーは次々と地殻内に漏れ出していく．結果として，氷床内をほぼ水平に伝播する波が選択的に残され，コンパクトな Le 波の波群が現れる．氷床から S 波エ
ネルギーが漏れ出していることは，図 10 のスナップショットで，地殻内に S 波を示す緑色が卓越していることから明瞭に確認できる．
3.4 TOPO-5
現実的な地形・氷床厚分布を用い，震源を氷床の下 5 km に置いた場合は，明瞭な Le 波が現れなかった．このケースでは，震源が深いことで，そもそも氷床内部に入射するエネルギーが少ない．また地形・基盤岩の凹凸によって氷床内へのエネルギーのトラップも弱められている．従って地殻内を伝播する波が卓越し，氷床内を伝播する波は不明瞭となる（図 6 赤線，図 11）．しかし氷床の影響は，氷床がない場合（図 6 黒線）に比べた Sg 波の増幅や後続波の継続時間の増大，といった現象に現れている．これは柔らかい地盤で地震動が増幅される現象と同じメカニズムである．図 6 の太線で示した部分は，群速度 3.1‐2.6 km/s の波群で，図 4 の観測波形の太線部分と対応している．今回の計算例は構造断面が観測のものと異なるため波形の直接比較はできないが，震源が氷床直下にない場合の波形の特徴は理論波形でよく再現できたといえる．

4.

まとめ

本稿では平成 26 年度の共同研究の成果として，局地的な地震波伝播を精度と効率よくモデリングする手法の開発と，グリーンランド氷床の現実的な地形・氷床厚分布を考慮して実行した地
震波伝播シミュレーションの成果を報告した．南極では以前から氷床上での地震観測が続けられているうえ，グリーンランドでも 2011 年から本格的な観測がスタートした．航空機によるメンテナンスや，衛星通信によるデータ転送によって，氷床上で得られた観測地震波形は急速に身近になりつつある．しかし氷床が観測地震波形に与える影響についての本格的な研究はほとんどな
く，試算が行われる場合も水平成層構造が用いられてきた．本研究では，水平成層構造モデルと現実的な地形・氷床厚モデルによるシミュレーションを比較して，構造や震源位置のわずかな違
いが，大きく異なる波形を生み出すことを明らかにした．シミュレーションの結果は次のようにまとめられる：(1) 震源が氷床直下にある場合，氷床の影響が波形に強く表れる；(2) 氷床内部での多重反射による地震波エネルギーのトラップ機構と，地形・基盤岩の起伏によるトラップ崩壊
機構の複合度合いによって，特徴の異なる波動場が形成される．さらに本研究では，氷床内にトラップされた S 波を「Le 波」，氷床表面を伝播する表面波を「Re 波」と命名した．氷床上で得られたデータの普及に伴い，氷床そのものの構造解析や氷床の影響の除去を目的として，観測波

— 22 —

SENAC Vol. 48, No. 4（2015. 10）形に氷床が与える影響の研究は重要な課題となりつつある．本研究は，世界に先駆けてその端緒を開く成果を挙げた．
なお本研究の計算は弾性波のみを扱っているが，実際の氷床や岩石は完全弾性体からずれた性質を持つため，波動場を減衰させる（＝非弾性減衰）．我々は本年，円筒座標系 2.5 次元差分法のプログラムに，非弾性減衰を導入した[29]．今後は氷床の非弾性減衰も考慮し，さらに多くの構造断面について大規模な計算を実行し，観測波形と理論波形の直接比較を行っていきたい．

謝辞
サイバーサイエンスセンターの小林広明センター長，江川隆輔准教授，小松一彦助教，大泉健治氏，山下毅氏には，本共同研究の機会をお与えいただいたうえ，プログラム開発やシステム利用について様々なご助言・ご助力をいただいた．GLISN 日本隊の活動は，科研費（課
題番号 24403006）により運営されている．記して感謝する．
参考文献[1] 豊国源知, グリーンランド氷床上での地震観測点設置, 極地, 94, 30-39, 2012.
[2] 豊国源知, GLISN 日本隊：2012 年の活動紹介, 極地, 96, 16-26, 2013.
[3] Toyokuni, G., Kanao, M., Tono, Y., Himeno, T., Tsuboi, S., Childs, D., Anderson, K.,
Takenaka, H., Monitoring of the Greenland ice sheet using a broadband seismometer
network: the GLISN project, Antarctic Record, 58(1), 1-18, 2014.
[4] 豊国源知・金尾政紀・東野陽子・姫野哲人・坪井誠司, GLISN 計画による日米共同地震観測隊の活動報告(2011～2014 年), 月刊地球, 37(3), 83-91, 2015.
[5] 豊国源知, 地震で調べる北極グリーンランドの環境変動, 科学フォーラム, 370, 28-33, 2015.
[6] Robinson, E.S., Seismic wave propagation on a heterogeneous polar ice sheet, J. Geophys.
Res., 73(2), 739-753, 1968.
[7] Toyokuni, G., Takenaka, H., Kanao, M., Tsuboi, S., Tono, Y., Numerical modeling ofseismic waves for estimating the influence of the Greenland ice sheet on observed
seismograms, Polar Sci., 9(1), 80-93, doi:10.1016/j.polar.2014.12.001, 2015.
[8] 豊国源知・竹中博士・趙大鵬・石原吉明, 月と火星の全球地震波伝播シミュレーション,
SENAC, 46(3), 1-12, 2013.
[9] Toyokuni, G., Takenaka, H., Wang, Y., Kennett, B.L.N., Quasi-spherical approach forseismic wave modeling in a 2-D slice of a global Earth model with lateral heterogeneity,
Geophys. Res. Lett., 32, L09305, doi:10.1029/2004GL022180, 2005.
[10] Toyokuni, G., Takenaka, H., FDM computation of seismic wavefield for an axisymmetricearth with a moment tensor point source, Earth Planets Space, 58(8), e29–e32, 2006.
[11] Toyokuni, G., Takenaka, H., ACE–A FORTRAN subroutine for analytical computation ofeffective grid parameters for finite-difference seismic waveform modeling with standard
Earth models, Comput. Geosci., 35(3), 635–643, doi:10.1016/j.cageo.2008.05.005, 2009.
[12] Toyokuni, G., Takenaka, H., Accurate and efficient modeling of global seismic wave
propagation for an attenuative Earth model including the center, Phys. Earth Planet. Int.,
200-201, 45-55, doi:10.1016/j.pepi.2012.03.010, 2012.
[13] Toyokuni, G., Takenaka, H., Kanao, M., Wiens, D., Nyblade, A., Comparison of global
synthetic seismograms calculated by the spherical 2.5-D finite-difference method withobserved long-period waveforms including data from intra-Antarctic region, Polar Sci., 6,155-164, doi:10.1016/j.polar.2012.06.001, 2012.
[14] Alterman, Z., Karal, F.C., Propagation of elastic waves in layered media by finite
difference methods, Bull. Seism. Soc. Am., 58(1), 367–398, 1968.

グリーンランド氷床モデルを用いた地震波伝播シミュレーション— 23 —[15] Stephen, R.A., A comparison of finite difference and reflectivity seismograms for marine
models, Geophys. J. R. Astr. Soc., 72(1), 39–57, 1983.
[16] Stephen, R.A., A review of finite difference methods for seismo-acoustics problems at theseafloor, Rev. Geophys., 26(3), 445–458, 1988.
[17] Igel, H., Djikpéssé, H., Tarantola, A., Waveform inversion of marine reflection
seismograms for P impedance and Poisson's ratio, Geophys. J. Int., 124(2), 363–371, 1996.
[18] Takenaka, H., Tanaka, H., Okamoto, T., Kennett, B.L.N., Quasi-cylindrical 2.5D wave
modeling for large-scale seismic surveys, Geophys. Res. Lett., 30(21), 2086,
doi:10.1029/2003GL018068, 2003.
[19] Toyokuni, G., Takenaka, H., Kanao, M., Quasi-axisymmetric finite-difference method forrealistic modeling of regional and global seismic wavefield-review and application. In: Kanao,
M. (Ed.), Seismic Waves, Research and Analysis. InTech, ISBN 978-953-307-944-8, 85-112,
2012.
[20] Dziewonski, A.M., Anderson, D.L., Preliminary reference Earth model, Phys.Earth
Planet. Int., 25(4), 297–356, 1981.
[21] Bamber, J.L., Ekholm, S., Krabill, W.B., A new, high-resolution digital elevation model ofGreenland fully validated with airborne laser altimeter data, J. Geophys. Res., 106(B4),6733–6745, 2001.
[22] Bamber, J.L., Layberry, R.L., Gogineni, S.O., A new ice thickness and 400 bed data set
for the Greenland ice sheet, 1. Measurement, data reduction, and errors, J. geophys. Res.,
106(D24), 33773–33780, 2001.
[23] Livingston, C.W., Explosions in ice, Technical Report 75, U.S. Army Snow, Ice, andPermafrost Research Establishment, Engineer Research and Development Center, pp 100,
1960.
[24] Robinson, E.S., Seismic wave propagation on a heterogeneous polar ice sheet, J. Geophys.
Res., 73(2), 739–753, 1968.
[25] Ishizawa, K., The measurement of the velocities of P and S waves propagating in thesurface layer of ice sheet at Mizuho Station, East Antarctica, Antarctic Record, 73, 147–160,
1981.
[26] Takenaka, H., Nakamura, T., Okamoto, T., and Kaneda, Y., A unified approach
implementing land and ocean-bottom topographies in the staggered-grid finite-difference
method for seismic wave modeling, Proc. 460 9th SEGJ Int. Symp., Sapporo,doi:10.1190/SEGJ092009-001.13, 2009.
[27] Nakamura, T., Takenaka, H., Okamoto, T., and Kaneda, Y., FDM simulation ofseismic-wave propagation for an aftershock of the 2009 Suruga Bay earthquake: Effects 440
of ocean-bottom topography and seawater layer, Bull. Seism. Soc. Am., 102(6), 2420–2435,
doi:10.1785/0120110356, 2012.
[28] Nakamura, T., Takenaka, H., Okamoto, T., and Kaneda, Y., Seismic wavefields in the
deep seafloor area from a submarine landslide source, Pure Appl. Geophys., 171, 1153–1167,
doi:10.1007/s00024-013-0717-3, 2014.
[29] Toyokuni, G., Takenaka, H., Zhao, D., Quasi-cylindrical 2.5-D wave modeling with ageneral point source and the anelastic attenuation, Phys. Earth Planet. Int., submitted, 2015.

— 24 —

SENAC Vol. 48, No. 4（2015. 10）[研究成果]［共同研究成果］
CFD による血管壁拍動挙動における Vasomotion の再現解析
工藤和貴、横堀壽光、西川１．

優、大見敏仁、門間良平、井上恭介(東北大学大学院工学研究科)

緒言

拍動流によって誘起される血管壁の動的挙動は、血管疾患に伴う壁の力学的性質変化により変化すると考えられ、その解明には流体解析と流体圧力により生じる血管壁の変形挙動の構造解析が必要である。
流体解析には、2 次元あるいは軸対称流れにおいては、渦度―流れ関数法、3 次元流れに対しては、MAC 法が一般的には用いられる。しかしながら、前者は連続の式を厳密に満足するが、血管壁の動的挙動の解明に必要な圧力を求める解析への適用が困難である[1]。また、MAC 法は、連続の式が近似的にしか満足されていないことが指摘されている[1]。
このような問題を解決するため、著者らは、渦度―流れ関数法と圧力のポアソン方程式を連成して解く方法により、血管壁の拍動挙動を解析している[6-8]。血管疾患に影響する血管拍動を構成する振動成分の特定は、より高精度に血管疾患を診断するために必要である。この手法により、拍動下での血管壁粘弾性発現度を定量化して動脈硬化の進行度を非侵襲で診断することが可能となる。
本稿では、著者のグループで検討・発展させてきた、渦度―流れ関数法と圧力のポアソン方程式の連成による血管壁の拍動挙動解析手法について述べる。
２．

解析方法

2.1 解析モデル
本研究で用いた、直管を模した解析モデルを図 1 に示す。流体は均質な非圧縮性 Newton 流体、流れは軸対称流とし、流れ方向および半径方向の二次元で解析を行った。

図 1 解析モデル
2.2 基礎式
二次元円筒座標系における連続の式と Navier-Stokes 方程式をそれぞれ（１）及び（２）式に示す。
���
��

���
��

+ ��
�

���
��

+ ��

+

� �

� ��

���
��

(�� � ) = 0

=−

� ��

� ��

+��

（１）

� � ��
��

+
�

� ���
� ��

+

�� ��
�� �

�

— 25 —CFD による血管壁拍動挙動における Vasomotion の再現解析

�� �
��

+ ��
�

�� �
��

+ ��

�� �

=−

��

� ��

+��

� ��

�� � ���+�� �� �
� ��

u+ : 軸方向速度m/s

v+ : 半径方向速度m/s

ft+ : 時間

s

p:

圧力

Pa=N/m2

ρ:

密度

kg/m3

ν:

動粘性係数

m2/s

+

�� � �
��

��

− ��
�
�

（２）

なお、添え字の+は有次元量であることを表す。管内流れは（１）式を満たすことから、（３）及び（４）式を満たす Stokes の流れ関数�が存在する。

�� = −� ��

�� = −

� ��

�=

−

また、渦度�は（５）式で定義される。

（４）

� ��

�� �
��

（３）

� ��

���

（５）

��

本研究では、流れ関数渦度法により流れ関数�と渦度�を計算することで流れ場の解析を行った。また、初期半径�� 、代表流速�� を用いて無次元化を行った。ここで、�� は（６）式に示す Reynolds 数の定義から、Re の値を定めることにより�� を決定した。

�� =

�� ��

（６）

�

以上により無次元変数及び無次元流れ関数�、無次元渦度�を導入することにより、（１）、（２）式
を連立することで以下の Poisson 方程式及び渦度輸送方程式を以下の（７）及び（８）式のように導
出できる。

��
��

+�

��
��

ΩR =

+�

��
��

�� �
��

−
�

�

:

管出口

:

管中心軸上

:

壁面上

:

��
��

��
��

��

��

= 0,
= 0,

��
��

��
��

� ��

− Ω=
�

また、境界条件は以下のように定義される。
管入口

� ��

= 0,

�� �

�� �

�

+

�

�� �
�� �

�� �

�� ��

（７）

+
�

�� �
�� �

+

�

�

� ��

�� �

（８）

=0

= 0, � = �00����= 0,� = 0� � = 0,� = 0,

��

��

=0

� = 0� � = 0� � = −��0�定常流),� = −��0�� − ����)�拍動流),� =� �� � ��,=0� �� � ��

ここで、壁面上の流れ関数を余弦関数により振動させることで、拍動を再現している。

解析において、血流の拍動に伴う血管の拡張収縮により、血管壁は直線ではなく複雑な曲線形状を

— 26 —

SENAC Vol. 48, No. 4（2015. 10）持つ。そのような非直線形状を持った領域における流れを取り扱う場合、長方形格子を用いると曲線
境界が階段形状に近似されるため、境界条件を正確に課すことは困難である。一方、境界に沿った格子を用いれば正確な境界条件を容易に課すことができる。そこで、座標変換により基礎方程式の変換を行い、解析領域を壁面形状で数値的に無次元化することにより物理領域である(R,Z)座標系から計算領域である(η,Z)座標系への変換による写像を施す。(R,Z)座標系における壁面形状を S(Z)で表すと、以下の（９）式で示される。

�=

�

�(�)

�� = �

（９）

これらを基礎方程式に適用すると、無次元基礎方程式である（７）及び（８）式は以下の（１０）及び（１１）式のように表される。

��
��

=

�

��

���

�� �
�� �

��

� ��

�� �
��

� ��
�

�� �
�� �

各係数は以下の通りである。

� ��

��
��

���

����

�� =

�� =

�� �

� ��

� ��

����

��
��

� ��

�� �
�� �

= �� Ω� ��� Ω� − ���1��(�) �
�
�1
�
�
�
� �
�(�)�
��

��
��

（１０）

− ���

��
��

− ��� Ω （１１）

�

1
1
��(�)
� � �(�)
�−
�
2�
�
�
−
��(�)
�
�� �
�(�)�
�
��
�� = −

2� ��(�)�(�) ��

�� = 1

�� = ��(�)1
��(�) ���� =
�1 � � �
� �
�(�)�
��
�� = 1

�� = −

2� ��(�)�(�) ��
�

1
1
��(�)
� � �(�)
�� =
�
�
2�
�
�
−
��(�)
�
�� �
�(�)� ������ = − �
��� =

��� =

1 �
�
��(�)

2 ��
��(�)� ��

��� = −2 ��
��(�)� ��2� ��(�) �� ���−��
��(�) �(�) �� �� ��（１０）、(１１)式を連立して�、�を計算することにより、Navier-Stokes 方程式の解を求める。なお、
�は SOR 法、�は陽解法により計算した。得られた�から、軸方向速度 u と半径方向速度 v は次のよ
うに求められる。

�=−

�

��

��(�)� ��（１２）

— 27 —CFD による血管壁拍動挙動における Vasomotion の再現解析

�=

�

�−�

��(�)

� ��(�) ��

�(�) ��

��

+

��
��

�

（１３）

次に、圧力の定式化について述べる。（２）式において圧力の半径及び軸方向の一次微分項を作り、（１４）式に示すポアソン方程式を求め圧力を計算した。
�
�
�
�
���
�� �
��� �� ��������� = − � − �� � + 2+� � �− �
��
��
�� ��
��
�
�
∆は軸対称流における円筒座標系の微分演算子で、（１５）式で表される。

∆=

��

�� �

+

� �

� ��

+

��

�� �

（１４）

（１５）

（１４）式における D+は、連続の式の値として以下のように定義する。

�� � � � � =

���
��

+

��
�

+

�� �
��

（１６）

本来、流れ関数渦度法では連続の式が満たされるが、圧力式を差分化して解く際に離散化による誤差が発生する。より高精度な計算のため、D を含む項のうち計算に影響を及ぼすと考えられる時間微分項を式中に残してある[1]。（１４）式を無次元化及び写像、差分化した最終的な計算式を（１７）式に示す。

���� =

0��� � ℎ�

1 + � � + �� �

��(�)
�
��

�

�

�

� � ��(�) ������ + ������
2� ��(�)� � � �(�) ������ − ������
�
�
�
�
+
�
�
�
−
�
� � ��ℎ������ �� �2ℎ����2� ��(�) �������� − �������� − �������� + �������� ������ + ������
�
�
−
+
4ℎ�
ℎ�
� ��
�
�
�
�
−
�
�
+
�
1
1
�
�
�
�����
�����
�����
�����
+ �
+ �
− ��
�
�
��
2ℎ
�
ℎ�
��
�
�
− ������ � 2 ������ − ������
1 �
� ��(�) ������ − ������ ������ − ������ ��� � � �����
�
� +
�−
+
�
��
��
2ℎ
2ℎ
2ℎ
2ℎ
�
� ��
�
�+��
�
�
���� �� ��(�) ������ − ������ ������ − ������
�
��
�
+ �−
+
� +� �
�
�
� ��
2ℎ
2ℎ
��
�
�
��
（１７）

ここで、D の差分化において、t について離散化し、��

れる圧力の値は、� ≠ 0であっても�

���

���

= 0と置いている。このようにして求めら= 0となるように定めているため、境界条件や時間進行による誤差集積に対して D を常に小さな値に抑えておくことが出来る[1]。この手法は MAC 法で用いられるものであるが、本研究における解析手法は、流れ関数渦度法から得られる流れ関数から流速を計算し、時間発展させているという点で MAC 法とは異なっている。
先に述べたように、本解析手法では、流れ関数の壁面境界条件により流量及び流速に変動が生じ、それに伴い管内圧力が変動する。そこで、血管壁が管内圧力に応じて変形する効果を考える。血管壁の剛性はひずみに対して非線形性を示すことが知られており、血管壁の剛性を（１８）式に示す、fungによる理論式[9]より求めた。

�=

�∗

�
�∗ � ∗������� − �� �����(� − �∗ )�
�:血管壁に作用する応力[g/mm2]

:

ひずみ

（１８）

— 28 —

SENAC Vol. 48, No. 4（2015. 10）（１８）式をひずみで微分した（１９）式を計算することで、剛性 E を算出した。
��
�∗
�
�
= ∗ � �� � � � � �� � ��� �����(� � �∗ )�
（１９）
��
�
�
� � ∗�
�

ここで a、� ∗ 、�∗ は定数であり、体重 15~25kg のイヌの胸部大動脈を用いた引張強度試験による実験結果[10]から、以下の値を用いた。

a

�

∗

�∗

:

1.9

:

25.83g/mm2:

1.6

内圧 p、壁厚 h、半径 R の薄肉円管モデルにおける円周応力は以下のようになる。

�=

��
�

（２０）

また、円周方向のひずみ�と応力�の関係は以下のようになる。

�� = ���

（２１）
以上より、半径は（２２）式で表される[11]。なお、以下のように漸化式の形とすることで、圧力変動に対し各ステップの半径を計算することができる。

��ℎ���� =

�� = ���� � ��

��

����
�� �

:

��

:

管壁の厚さ m:

圧力の増分 Pa=N/m2（２２）

管壁の剛性増分 Pa=N/m2

壁面形状関数として定義された S(Z)は（２２）式により計算される。これにより、各軸方向位置における半径を用いて写像を行うことで、計算領域を長方形格子に変換している。
次に、解析のフローチャートを図 2 に示す。第一節で述べたように、本解析では渦度流れ関数法と圧力のポアソン方程式を連成して解く手法を取っている。一般的に、レギュラー格子を用いた差分法
による流体解析では、求めた圧力から流速を時間発展させることにより連続の式が完全には満足されないという問題がある。しかし、本解析手法では流れ関数渦度法により時間発展させるため、レギュ
ラー格子を用いながらも連続の式を十分に満たし、かつ圧力を計算することが可能となる。

— 29 —CFD による血管壁拍動挙動における Vasomotion の再現解析

図 2 フローチャート
2.3 解析条件
2.3.1

Vasomotion のモデル化

血管運動には、心拍動の周期とは異なる、自律神経活動に起因する自然発生的で律動的な長周期的収縮拡張運動(Vasomotion)が存在する[12]。一般に、Vasomotion の生理学的意義については多くの研究がなされている[13]。しかし、Vasomotion が血管壁拍動挙動に及ぼす影響についてはまだ解明されていない。そこで、Vasomotion が関与する血管壁運動を数値解析において再現するため、拍動の周波数
�� より低い周波数�� で振動する関数を用いて、（２３）式で示す拍動条件を考えた。
��Ψ���� = −1.0 × (� − ����� �) × �

��

−

�

��

����� ������ = �� �
�（２３）

この拍動条件を用いて、Vasomotion が関与する血管壁の拍動挙動に与える影響について検証した。
2.3.2 粘弾性モデル
血管壁は弾性的性質が劣化すると、その因子であるコラーゲン等が有する粘弾性的性質により、ク
リープや応力緩和に起因する粘弾性的特性が発現する。この時、血管壁がどのような拍動挙動を示す
かを明らかにすることは血管疾患の診断において重要である。本解析を行うために、血管壁の応力ひ— 30 —

SENAC Vol. 48, No. 4（2015. 10）

ずみ構成則として図 3 に示すような粘弾性ばねモデルを構築した。このモデルは、粘弾性 Voigt モデ
ルと Maxwell モデルを組み合わせたものであり、モデル血管における血管壁のクリープ特性と応力緩和特性を同時に再現している[5]。この粘弾性モデルに引張応力を負荷した場合の構成式は、（２４）、（２５）式で表される。
�σ = �� � � �� �� � ��� �� � ���
E=

��
��

�

��

= �� � �� ��� �� � ���� � ��
�

�

��

:

弾性係数

:

粘性係数

:

ひずみ速度

（２４）
（２５）

図 3 粘弾性ばねモデル
ところで、粘弾性特性には、応力負荷時と除荷時で応力ひずみ曲線が異なる経路をたどる、ヒステリ
シス特性が存在する。より詳細に血管壁粘弾性をシミュレートするため、ヒステリシス特性を再現する条件を導入した。除荷条件として、最大ひずみ���� からの減少量���� を導入し、（２４）式をΔσとして最大負荷応力���� からの減少量と定めた。これにより、負荷曲線を 180°回転した形の除荷曲線が決定し、ヒステリシス特性が再現される。
2.4 カオス理論に基づく不安定挙動解析

過去の研究において、血管壁粘弾性の発現により血管挙動の不安定性が増加することが臨床において示されている[4]。著者の一人が提案してきた DPC 法[2-5]では、測定した血管壁ひずみ速度から動脈硬化進行度を定量的に評価している。さらに、血管壁挙動の振動成分をより詳細に解析するために、
アトラクタ解析による軌道不安定性の定量化が行われている[4]。図 4～6 に示すように、DPC 出力波形に対しアトラクタ解析を行うことで、血管壁挙動の拍動軌跡が得られる。また、以下の式で示されるエントロピーS を計算することで、軌道の不安定性を定量化することが出来る。

— 31 —CFD による血管壁拍動挙動における Vasomotion の再現解析(a) 健常例(b)動脈硬化症例

図 4 超音波ドップラー法で検出された血管壁拡張時の血管壁変形速度非侵襲測定例[2]

(a)

(b)

V(t+⊿t)V(t+⊿t)

図 5 拍動下の(a)血管壁速度非侵襲測定結果と(b)拍動軌跡特性[4]

V(t)

V(t)

(a) 健常例(b)動脈硬化症例
図 6 血管壁の拍動軌跡[4]
�

�� :

� :

� � � � �� �����
���

i 番目の領域に軌道が落ちる確率
領域分割数

（２６）

— 32 —

SENAC Vol. 48, No. 4（2015. 10）

健常者において、血管壁の拍動挙動には図 6(a)に示されるように、サーキットループ(リミットサイクル)の軌跡に長周期的変動が重畳されていることが分かる。動脈硬化症の血管においては、このループが半楕円状になり、また拍動軌道に乱れが生じていることが分かる。また、血管壁拍動軌道の乱れがエントロピーの増加として発現している [4] 。本稿では、数値解析結果に対し血管壁拍動軌跡を解析し、臨床結果との比較を行った。
３．

解析結果

図 7～9 に弾性変形モデルと Vasomotion を考慮した弾性変形モデル、粘弾性モデルのそれぞれについて数値解析を行い、得られた血管壁ひずみ速度波形及び拍動軌跡を示す。また、これらの各解析条
件における、血管壁ひずみ速度波形から求めたエントロピーS を表１に示す。
弾性変形モデルでは、血管壁ひずみ速度波形は健常例のような正弦波の形を示し、また波形の経時的変化も見られなかった。このような定拍動振幅条件においては、拍動軌道は図 7(b)のようにただ一つの軌道を示すことが解析により示された。
Vasomotion を考慮した弾性変形モデルでは、通常の弾性変形モデルと比べ血管壁ひずみ速度波形の形状的変化は見られなかったが、時間進行によりその振幅が異なる波形が得られた。この速度波形に対する拍動軌跡には、図 8(b)に示されるような長周期的なゆらぎが確認された。このように、
Vasomotion の導入により拍動軌跡に長周期的軌道が重畳され、図 6(a)で示すような健常者の臨床例に一致することが分かった。すなわち、健常血管においても血管壁拍動挙動に Vasomotion が存在することが CFD により検証された。
次に、粘弾性効果を導入したモデルでは、図 9 に示すように時間進行に伴う血管壁ひずみ速度波形の拍動振幅に変動がみられ、また拍動軌道のゆらぎが確認された。このモデルのように、定拍動振幅
条件においても粘弾性の効果により粘弾性特有の Vasomotion が発現することが解析により示された。
さらに、Vasomotion を考慮した粘弾性モデルでの解析結果を図 10 に示す。この解析条件では、通常の Vasomotion によるゆらぎに対し粘弾性に起因する Vasomotion が重畳され、より顕著な軌道の乱れが確認された。また表 1 から、粘弾性の発現による軌道の乱れがエントロピーの増加として反映されていることが示された。
４．

結言

渦度－流れ関数法と圧力のポアソン方程式を連成し血管壁の拍動挙動を解析する手法により、健常な血管壁拍動における長周期的振動の存在と、壁粘弾性発現による特有の長周期的振動という二つのVasomotion の存在を再現できた。今後、さらに本手法を拡張発展させることで、詳細な血管疾患との対応が可能となるだろう。

(10-5)V(t+⊿t)Wallspeed

3

2

1

0
500000Timestep

600000(a) 血管壁ひずみ速度波形
図 7 弾性変形モデル

V(t)(b) 拍動軌跡

— 33 —CFD による血管壁拍動挙動における Vasomotion の再現解析(10-5)4

V(t+⊿t)Wallspeed

3

2

1

0
500000600000TimestepV(t)(a) 血管壁ひずみ速度波形
図8(b) 拍動軌跡

Vasomotion を考慮した弾性変形モデル(10-5)2

V(t+⊿t)Wallspeed310500000600000TimestepV(t)(a) 血管壁ひずみ速度波形(b) 拍動軌跡

図 9 粘弾性モデル(10-5)V(t+⊿t)Wallspeed

3

2

1

0
500000Timestep

600000V(t)(a) 血管壁ひずみ速度波形
図 10(b) 拍動軌跡

Vasomotion を考慮した粘弾性モデル
表 1 エントロピー
解析モデル

エントロピー

弾性変形モデル

4.417

Vasomotion を考慮した弾性変形モデル

4.942

粘弾性モデル

4.817

Vasomotion を考慮した粘弾性モデル

5.186

謝辞
本解析手法の開発にあたり、東北大学サイバーサイエンスセンターの SX シリーズを使用させて頂いたことを感謝する。

— 34 —

SENAC Vol. 48, No. 4（2015. 10）

参考文献[1] 河村哲也

流体解析１、1996,朝倉書店[2] A.T.Yokobori,Jr., T.Ohkuma, S.Sasak, H.Yoshinari, T.Yokobori, H.Ohuchi and S.Mori,
Bio-Medical Materials and Engineering, 4, (1994), pp.87−96.
[3] A.T.Yokobori,Jr., M.ichiki, H.Ohuchi, T.Kobayashi, T.Satoh and Y.Kinoshita, Bio-Medical
Materials and Engineering, 14,3 (2006), pp.241−249.
[4] A.T. Yokobori,Jr., M. Owa, M. Ichiki, T. Satoh, Y. Ohtomo, Y. Satoh, S. Ohogoshi, Y. Kinoshita
and S. Karino, Journal of Atherosclerosis and Thrombosis, 13,4, (2006),pp163-174.
[5] A.T.Yokobori,Jr., T.Ohmi, R.Monma, Y.Tomono,K.Inoue, M.Owa, M.Ichiki, N.Mochiduki and
H.Yamashita, Bio-Medical Materials and Engineering, 23, (2013), pp.75−91.
[6] A.T.Yokobori.Jr., Y.Nishikawa, The analysis of pulsatile behavior of blood vessel wall withaneurysm based on hybrid analysis of computer fluid dynamics (CFD) and solid mechanics,
Proc. of the Ann Meeting of Atherosclerosis and Thrombosis, (2013)[7] A.T.Yokobori.Jr., Y.Nishikawa, The analyses of mechanical behaviors of blood vessel wall withaneurysm at the progressive stage based on simulate experiment and computer fluid
dynamics(CFD), Proc. of the Ann Meeting of Atherosclerosis and Thrombosis, (2014)[8] A.T.Yokobori.Jr., K.Kudo, Long Periodic Fluctuation of Blood Vessel Wall Caused by the
Viscoelasticity Under Pulsatile pressure Flow, Proc. of the Ann Meeting of Atherosclerosis
and Thrombosis, (2015)[9] Y.C.Fung, Appl. Mech. Rev., 1968,21(1),pp1-20.
[10] 横堀壽光，大熊恒郎，前山俊秀，横堀武夫，大内博，奈良茂樹，佐々木久雄，葛西森夫，生体材
料，1986,4,pp27-32.
[11] 岡小天，ながれ，1993，12，pp249-254.
[12] 南山求，中野厚史，日本レーザー医学会誌，1994,9,pp27-38.
[13] T.Sakurai, N,Terui, Am J Physiol Heart Circ Physiol 291, (2006), pp.1761-1767

SENAC Vol. 48, No. 4（2015. 10）— 35 —［共同研究成果］[共同研究成果]

大規模マルチビームリフレクトアレーの設計法とベクトル型スーパーコンピュータによる高速化の研究
今野 佳祐 † , 陳 強 †東北大学大学院工学研究科 通信工学専攻 †

まえがき

1

高利得・広帯域なアンテナの 1 つとして, パラボラアンテナがこれまで用いられてきたが, 大型でかさばるという欠点があった. そこで近年, パラボラアンテナに代わる新たなアンテナとして,
平面型リフレクトアレーが大きな注目を集めている [1]. 平面型リフレクトアレーは, 電磁波を送受信する一次放射器と, 多数の素子から構成される反射面から成る [2]. 平面型リフレクトアレーの反射面上の素子は, 一次放射器から入射した電界の位相を回転させ, 所望の方向で平面波を形成するように設計される. 入射電界の位相は素子の位置によって変わることから, 反射面上にあるリフレ
クトアレー素子は 1 つ 1 つ異なる大きさにしなければならない. 従って, 反射面状の素子サイズを決定する前に, リフレクトアレーの素子の大きさに対する反射電界の位相変化量, すなわち反射係
数の位相の情報が必要となる. これまで筆者らは，反射係数の位相を求める手法として起電力法を用いた手法を提案し，その有効性を明らかにしてきた [3]-[7]. 筆者らの提案手法は，スペクトル
領域モーメント法 [8]-[10] や有限要素法 [11]-[14] と比べ，一次放射器を厳密に考慮できるという点で優れている．
一方, 衛星通信やレーダー等への応用が期待されることから, マルチビームリフレクトアレーに関する研究が盛んに行われている. これまで, 様々なマルチビームリフレクトアレーの設計法が提案されてきた [15] が, リフレクトアレー素子と一次放射器の相互結合を厳密に考慮した設計法は提案されていない. また, 素子間相互結合を考慮してマルチビームリフレクトアレーを設計するのは非常に長い計算時間がかかるため, 設計法の高速化が欠かせない.
本報告では, 文献 [3]-[7] に示されたシングルビームリフレクトアレーの設計法をベースにした，マルチビームリフレクトアレーの設計法を提案する．また，提案法をベクトル型スーパーコンピュータを用いて高速化するためのチューニング法を明らかにする．最後に，大規模マルチビームリフ
レクトアレーの設計を行い, その有効性を数値的に明らかにする.

2起電力法を用いたリフレクトアレー設計法

2.1

線状素子シングルビームリフレクトアレーの設計法

文献 [3]-[7] の設計法によって, 線状素子のシングルビームリフレクトアレーは以下のように設計される.

— 36 —

SENAC Vol. 48, No. 4（2015. 10）

1. リフレクトアレー素子数, 素子間距離, 反射板と素子間の距離や一次放射器の位置など, リフ
レクトアレーの基本的な構造を決定する.

2. 起電力法を用い, 各位置にあるリフレクトアレー素子の電流を素子長に対して数値計算する.
このとき, 反射板の大きさは無限とし, 鏡像法を用いる. また, 素子の電流を計算するとき,
リフレクトアレー素子間の相互結合は無視する.

3. 所望の散乱方向 (θd , ϕd ) を決定する. そして, (2) で数値計算した電流を用い, 各位置にある素子の (θd , ϕd ) 方向の散乱電界の位相と素子長との関係を求める.

4. 任意の位置にある素子の長さを適当に固定し, その素子が (θd , ϕd ) 方向へ散乱する電界の位
相 P を位相の基準とする. その他の位置にあるリフレクトアレー素子の長さは, (θd , ϕd ) 方
向へ散乱する電界の位相が P となるように決める.
提案設計法の利点は 2 つ挙げられる. 1 つは, 一次放射器からリフレクトアレー素子への入射電界の角度や距離が, リフレクトアレー素子の散乱電界の位相に与える影響を厳密に考慮できることである. もう 1 つは, リフレクトアレー素子を線状素子とし, かつリフレクトアレー素子間の相互結合を無視しているため, 設計に要する時間が短いことである.

2.2

線状素子マルチビームリフレクトアレーの設計法

線状素子マルチビームリフレクトアレーの設計法を図 1 に示す. 提案設計法は，起電力法によって設計したシングルビームリフレクトアレーを初期構造とし，素子長を変えながらマルチビーム
リフレクトアレーを設計していくアルゴリズムになっている．提案設計法では，まず初期構造を設計した後に 2 つ目のメインビーム方向を決める．次に，ある素子の素子長をわずかに増減させ，そのときの 2 つのビーム方向の利得を計算する．このとき，素子間相互結合の影響を厳密に考慮するため，リフレクトアレーの電流は前処理付き共役勾配法で求める．そして，素子長を変える
前後での 2 つ目のビーム方向の利得をそれぞれ比較する．素子長を増やしたときに利得が大きくなった場合は，利得が増え続けるか素子長が 0.5λ に達するまで素子長を増やして同様の計算を行う．一方，素子長を減らしたときに利得が大きくなった場合は，利得が増え続けるか素子長が 0.1λに達するまで素子長を減らして同様の計算を行う．これを全ての素子について繰り返し，2 つのビームの利得が同じ値になったら終了する．
リフレクトアレーの数値シミュレーションが悪条件問題であることを考慮し，1 列分の素子長が更新し終わった段階で前処理行列を求め直すこととする．これにより，リフレクトアレーの素子長が変わっても良好な収束性を維持することができる．

3

ベクトル型スーパーコンピュータによる高速化
ベクトル型スーパーコンピュータは, プログラムのループ中で繰り返し処理される配列データの演算を一括実行するというベクトル演算機能を持っている [16]. 加えて, プログラムのループ中の

— 37 —

大規模マルチビームリフレクトアレーの設計法とベクトル型スーパーコンピュータによる高速化の研究

1. シングルビームリフレクトアレーの設計.
シングルビームリフレクトアレーの設計(第一・第二メインビーム方向利得G1(lmi), G2(lmi)を計算. ただしG1>>G2)
2.前処理行列の計算
前処理行列の計算.
前処理行列の計算(悪条件問題のため, Full-sizeの前処理行列を利用. )3. 第i素子長
素子長lだけ変化させ,素子長 miを∆lだけ変化させだけ変化させii
G2(lm +∆l)およびおよびGを前処理付き共役勾配法で計算.
および 2(lm -∆l)を前処理付き共役勾配法で計算を前処理付き共役勾配法で計算(y軸に対して線対称な位置の素子長も同時に変化. )If G2(lmi+∆l) >G2(lmi)
If G2(lmi-∆l) >G2(lmi)iiand G2(lm +∆l) >G2(lm -∆l)and G2(lmi-∆l) >G2(lmi+∆l)
Increment of n until G2 =G1 ,or until G2(lmi+n∆l)→max.
where lmi+n∆l<0.5λ

Increment of n until G2 =G1 ,or until G2(lmi-n∆l)→max.
where lmi-n∆l>0.1λOtherwise

No increment
If G2 =G1 or i=MxMy/2
If G2 =G1
ori=MxMy/2

Elseif i=My, 2My,…,

Elseif i=My, 2My,…,4. 1列分の素子長更新が列分の素子長更新が終了する度に前処理行列の更新.
終了する度に前処理行列の更新(収束性向上のため. )4. 設計したリフレクトアレーの数値シミュレーション.
設計したリフレクトアレーの数値シミュレーション

図 1: マルチビームリフレクトアレーの設計法.
繰り返し処理を並列実行する並列演算機能も備えている. ベクトル演算は最も内側のループに, 並
列演算は外側のループに適用されるため, 多重ループ内で最も内側のループ長が長くなるように
チューニングを行えば, ベクトル演算・並列演算によって計算時間の大幅な短縮が可能となる. ベ
クトル演算による高速化のためのチューニングの概念を図 2 に示す.
マルチビームリフレクトアレーの設計法では, 素子長を変えた時のインピーダンス行列要素の更新, 前処理付き共役勾配法による行列方程式の求解, リフレクトアレー遠方界の演算といった 3 つの処理を繰り返し行う. 前処理付き共役勾配法で最も大きな計算時間を要するのは行列-ベクトル積の演算であり, 行列-ベクトル積は, ループ長が行列サイズと同じ長さの 2 重ループで演算される. したがって, 行列方程式が大規模であれば, 自動的に最も内側のループはベクトル演算で, 外側のループは並列演算によってそれぞれ処理され, 行列方程式の求解は自動的に高速化される. その一方で, インピーダンス行列要素の更新と, リフレクトアレー遠方界の演算を行う部分はそれぞれループ長の異なる多重ループで処理されるため, ベクトル演算・並列演算による高速化を実現する
ためにはチューニングが必要となる. 大規模リフレクトアレーでは素子数が非常に多くなること

— 38 —

SENAC Vol. 48, No. 4（2015. 10）

Low vector operation ratio
Do I=1, N1, 1
Do I=1, N2, 1

ENDDO
ENDDO

TuningHigh vector operation ratio
Do I=1, N2, 1
Do I=1, N1, 1

ENDDO
ENDDO

where N1>N2

where N1>N2

図 2: ベクトル演算による高速化のためのチューニングの概念.

y
Reflectarray element

x

d

lm or

lm

lp

2ap
2am
Dipoledx

2am
Parasitic dipolezdy

Dipoleantenna

Mx

lt
t

My

h

2at

Primary source
Infinite ground plane

図 3: 線状素子リフレクトアレーの構造.
から, 本報告では素子数に対応するループが最も内側になるようにチューニングを行い, インピー
ダンス行列要素の更新およびリフレクトアレー遠方界の演算を行う部分をベクトル演算・並列演
算によって高速化する.

4

数値シミュレーション結果
提案設計手法により, 大規模なシングルビームおよびマルチビームリフレクトアレーを設計した. リフレクトアレーの構造を図 3 に示す. リフレクトアレーは, ダイポール素子および寄生素子付きダイポール素子から成り, 背面には反射板を設けている. 本報告の数値シミュレーションは東北大学サイバーサイエンスセンターのベクトル型スーパーコンピュータ SX-9 を利用して 8 並列で行われ, ベクトル演算率は 99.5%であった. なお, 本報告では E 面でビーム走査を行うリフレクト
アレーを設計したため, リフレクトアレーは H 面について対称な構造を持ち, 設計にも対称性を利用している.
設計した 50 × 50 素子シングルビームリフレクトアレーの指向性を図 4 に，開口面の位相分布

— 39 —

大規模マルチビームリフレクトアレーの設計法とベクトル型スーパーコンピュータによる高速化の研究

50

Directivity [dBi]f=20 GHz, h=364 mm, t=1.7 mm, am=ap=0.45 mm, lp=0.97lm
d=1.9 mm, Mx=My=50, dx=dy=10.5 mm, lt=6.6 mm
40 ap=1.075 mm, (θd, φd)=(20o, 90o ), Eθ on yz plane
30

20

10

0

-10
30

20

10

0

θ [deg.]

10

20

30

図 4: シングルビームリフレクトアレーの指向性.

図 5: シングルビームリフレクトアレーの開口面位相分布.
を図 5 にそれぞれ示す. 素子ごとに長さを異なる値にすることで，図 5 に示すような不均質な開口面位相分布を実現し，その結果所望の (20◦ , 90◦ ) 方向に強いビームが出ていることが分かる．実
効開口効率は 36.2%であり，設計と数値シミュレーションに要した時間は 375 秒であった．なお，

— 40 —

SENAC Vol. 48, No. 4（2015. 10）

50

Directivity [dBi]f=20 GHz, h=364 mm, t=1.7 mm, am=ap=0.45 mm, lp=0.97lm
d=1.9 mm, Mx=My=50, dx=dy=10.5 mm, lt=6.6 mm
40 ap=1.075 mm, (θd , φd)=(20o, 90o), (10o , 270o), Eθ on yz plane
Single-beamMulti-beam

30

20

10

0

-10
30

20

10

0

θ [deg.]

10

20

30

図 6: マルチビームリフレクトアレーの指向性.

図 7: マルチビームリフレクトアレーの開口面位相分布.
同様のリフレクトアレーを Core i7 搭載の PC で設計し，その数値シミュレーションをするのに要した時間は 60,720 秒であり，CPU 比 162 倍の高速化を実現した．
次に，設計した 50 × 50 素子マルチビームリフレクトアレーの指向性を図 6 に，開口面の位相大規模マルチビームリフレクトアレーの設計法とベクトル型スーパーコンピュータによる高速化の研究— 41 —

分布を図 7 にそれぞれ示す. 図 6 から, 所望の方向 (20◦ , 90◦ ) および (10◦ , 270◦ ) に 2 つのビームがそれぞれ形成されていることが分かる. シングルビームリフレクトアレーとは各素子の長さが異なり，結果的に開口面の位相分布が変わっていることが図 7 から分かる．なお，2 方向にビー
ムを向けているため, 1 方向にビームを向けたリフレクトアレーに比べて実効開口効率が 12.5%に下がり, 指向性利得はその結果 4.8 dB 低下している. 設計に要した時間は 20,097 秒であり，設計
中にリフレクトアレー構造全体の数値シミュレーションを 5,092 回行ったことから，1 回あたりの数値シミュレーション時間は約 4 秒である．リフレクトアレー構造全体を解析するとき，未知数

N = 10, 610 の行列方程式を解く必要があるが，PC でこの規模の行列方程式を解くと最低でも数時間はかかる．したがって，スーパーコンピュータによる高速化が，マルチビームリフレクトア
レーの設計には不可欠である．

5

むすび
本報告では, 大規模なマルチビームリフレクトアレーの設計法を提案し，ベクトル型スーパーコ

ンピュータによる高速化チューニングを行った．数値シミュレーションによって大規模なマルチビームリフレクトアレーを設計し, 提案法の有効性を数値的に明らかにした.

謝辞
本研究は, 東北大学サイバーサイエンスセンターのスーパーコンピュータ SX-9 を利用して行われたものである. スタッフを始め, 関係各位に深く感謝する.

参考文献[1] J. Huang, “Analysis of a microstrip reﬂectarray antenna for microspacecraft applications,”
TDA Progress Report 42-120, pp. 153-173, Feb. 1995.
[2] J. Huang and J.A. Encinar, Reﬂectarray Antennas, John Wiley & Sons, 2008.
[3] 今野 佳祐, 陳 強, 澤谷 邦男, 亀田 卓, 末松 憲治, “起電力法によるリフレクトアレー設計法の一検討,” 信学技報, vol. 112, no. 491, AP2012-164, pp. 1-5, 2013 年 3 月.

[4] 今野 佳祐, 陳 強, 澤谷 邦男, 亀田 卓, 末松 憲治, “起電力法を用いた線状素子リフレクトア
レーの設計,” 信学総大, B-1-138, p. 138, 2013 年 3 月.

[5] K. Konno, Q. Chen, K. Sawaya, S. Kameda and N. Suematsu, “Novel design method forreﬂectarray by induced electromotive force method,” Proc. IEEE AP-S Int. Symp., 429.3,
pp. 1342-1343, July 2013.

— 42 —

SENAC Vol. 48, No. 4（2015. 10）[6] K. Konno, Q. Chen, S. Kameda and N. Suematsu, “Design of ﬁnite FSS-backed reﬂectarrayby using BDP-CG method,” Proc. iWAT2014, pp. 200-202, March 2014.
[7] 今野 佳祐, 陳 強, “アレー給電による大規模な低姿勢リフレクトアレーの高利得化,” 信学技報, vol. 114, no. 294, AP2014-138, pp. 55-59, 2014 年 11 月.

[8] D.M. Pozar and T.A. Metzler, “Analysis of a reﬂectarray antenna using microstrip patchesof variable size,” Electr. Lett., vol. 29, no. 8, pp. 657-658, April 1993.
[9] C. Wan and J.A. Encinar, “Eﬃcient computation of generalized scattering matrix for analyzing multilayered periodic structures,” IEEE Trans. Antennas Propag., vol. 43, no. 11,
pp. 1233-1242, Nov. 1995.
[10] M. Zhou, E. Jørgensen, O.S. Kim, S.B. Sørensen, P. Meincke, and O. Breinbjerg, “Accurate
and eﬃcient analysis of printed reﬂectarrays with arbitrary elements using higher-orderhierarchical Legendre basis functions,” IEEE Antennas Wireless Propag. Lett., vol. 11, pp.
814-817, 2012.
[11] R. Remski, “Analysis of PBG surfaces using Ansoft HFSS,” Microw. J., vol. 43, no. 9, pp.
190-198, Sept. 2000.
[12] I. Bardi, R. Remski, D. Perry, Z. Cendes, “Plane wave scattering from frequency-selective
surfaces by the ﬁnite-element method,” IEEE Trans. Magn., vol. 38, no. 2, pp. 641-644,
March 2002.
[13] L. Li, Q. Chen, Q. Yuan, K. Sawaya, T. Maruyama, T. Furuno, and S. Uebayashi, “Novel
broadband planar reﬂectarray with parasitic dipoles for wireless communication applications,” IEEE Antennas Wireless Propag. Lett., vol. 8, pp. 881-885, 2009.
[14] Q.-Y. Li. Y.-C. Jiao, G. Zhao, “A novel microstrip rectangular-patch/ring-combination
reﬂectarray element and its application,” IEEE Antennas Wireless Propag. Lett., vol. 8,
pp. 1119-1122, 2009.
[15] P. Nayeri, F. Yang, and A. Z. Elsherbeni, “Design of single-feed reﬂectarray antennas withasymmetric multiple beams using the particle swarm optimization method,” IEEE Trans.
Antennas Propag., vol. 61, no. 9, pp. 4598-4605, Sept., 2013.
[16] サイバーサイエンスセンター スーパーコンピューティング研究部, “スーパーコンピュータシ
ステム SX-9 利用ガイド”, p. 2, 2014 年 4 月.

— 43 —

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

アンテナの放射効率低下メカニズムに関する一検討

東京農工大学

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

アンテナは，電磁波と通信信号のインターフェイスを果たす重要なデバイスである．一方，良導体のみで構成したアンテナにおいてもその条件により，アンテナ放射効率の低下が起きる事が知られている．このメカニズムは解明されていない．本研究では，このメカニズム解明のために，
大型計算機を用いた大規模シミュレーションにより解明するものである．本報告では，アンテナの放射効率の要因の一つであると考えられる，良導体の表面粗さを考慮した解析を行う．そして，導体の表面粗さが，その導電率に及ぼす影響について検討した．導体の表面粗さは実際に測定した値を用い，放射効率に関しては，FDTD 法に表面粗さを組み込んでシミュレーションを行い評価
している．

1.

はじめに

スマートフォンなどの携帯情報端末をはじめとする，電磁波を用いた機器の増加により，アンテナの重要性はますます高まっている．アンテナの放射効率は，入力電力に対してどれだけ電磁波を放射したかの比で定義されるが，当然その効率はより高いものが望まれる．そのため，アンテナは高い放射効率を得る目的で良導体を用いて構成される．一方，これまでの研究により，良導体のみで構成したアンテナにおいてもその条件により，アンテナ放射効率の低下が起きる事が知られている[1,2]．このメカニズムは十分に解明されているとは言えない[3,4]．一方，アンテナを構成する導体板の表面の粗さがその導電率に及ぼす影響について議論されている．すなわち導体自体の導電率は高いのであるが，その表面の粗さによって導電率が低下することが指摘されている．アンテナにおいては，放射において表面上に電流が集中することからその影響は大きいと考えられる．文献[5]においては，導体の表面粗さを等価的な式で導電率の変化を表している．
本研究では，文献[5]と同様にアンテナを構成する導体板の表面粗さを，導電率の変化に置換え，その放射効率に対する影響を検討した．導体板の表面粗さはいくつかのサンプルを用いて実際に測定した値を用いた．また，放射効率については，表面粗さを導電率の変化に置換えた値を用いて，FDTD 法[6]を用いたシミュレーションにより評価した．FDTD 法の解析時間の高速化について
は，東北大学サイバーサイエンスセンターの協力を得た．

２.

導体板の表面粗さの測定結果

本節では，実際の導体板の表面粗さを調べるために，レーザマイクロスコープと呼ばれる光学機器を用いて導体板表面の粗さを測定し，二乗平均粗さ(Root Mean Square: RMS)を算出する．RMS
は次式で定義される．
��

RMS = � �� ����� ���

(1)

ここで f(x) は基準としたところからの高さであり，この変化が粗さに相当する．レーザマイク
ロスコープは KEYENCE 製 VK-X100 を用い，RMS の算出には KEYENCE 製の VK Analyzer というソフ

— 44 —

SENAC Vol. 48, No. 4（2015. 10）

トウェアを用いた．本研究では，誘電体基板上に作製されるアンテナを想定しているために，測定は新品の両面銅箔誘電体基板（厚み 1.6 mm）を 4 つのサンプルで測定した．サンプルはそれぞれ，サンプル A～サンプル D と区別する．未使用の状態のサンプル A，サンプル B，数カ月屋外に放置したサンプル C，サンプル D について表面粗さを測定した．その結果，サンプル A が 0.309 m，
サンプル B が 0.337 m，数カ月屋外に放置したサンプル C，サンプル D は共に 0.396 m となった．銅箔基板の表面粗さは，雨などにさらされることにより大きくなることが分かった．これは錆や誘電体基板の特性変化によると考えられる．各サンプルについて 234 m×175 m の範囲で表面粗さを測定し，位置が高い部分を赤く，低い部分を青く表示した結果をそれぞれ図 1 に示す．
図 1 より，周期的に山と谷が分布しているのが確認でき，およそ 1m 周期で山と谷があるのが確認できる．最大値と最小値の差は，およそ 3m であり，1m の周期で表面粗さが発生していることが分かる．これは，銅箔基板が作製される際の工程によるものと推測される．また，図より全体的な周期の中に複雑な粗さが混在していることが分かる．

図 1： 導体板表面の粗さ測定結果

３.

表面粗さと損失の関係

本研究では文献[5]にならって，導体の表面粗さによる導電率の低下としてモデル化する．文献[5]によると，粗さによる損失の増加係数αを表面が滑らかな場合の抵抗率α0 で表すと，次式のようになる．
α = α� K ��

(2)

�

∆

K �� = � � arctan ���� � ����� � (3)��

ここで，Ksr は Hammerstad と Jensen によりマイクロストリップ線路の計測から得られた式[5]，∆は導体表面の粗さの RMS(二乗平均平方根)，は表皮深さを表す． は周波数と導電率などによって変化し  

2 / ()の関係がある． (2)，(3) 式より，∆が大きく，が小さくなるほど導体の損失係数が増加することがわかる．この変化を図 2 に示す．

— 45 —

アンテナの放射効率低下メカニズムに関する一検討

7 mm
75 mm

75 mm

図 2：

Ksr の変化V(t)図 3 ループアンテナモデル(1GHz)

図 2 では∆を 0.2m～1.6m まで変化させた際の Ksr を示している．図より∆と周波数が大きくなるにつれて，Ksr がおよそ 2 に収束しているのがわかる．これは，周波数が十分に大きくなることで完全に導体表面の粗さに沿って電流が流れるようになり，これ以上周波数が増加しても電流の移動距離が増加しなくなるからであると考えられる．この式を用いることで，使用する導体の透磁率，導電率，表面粗さの RMS をパラメータとして入力すれば，導体の周波数による導電率の変化を求めることが可能となる．

４.

ループアンテナの放射効率と表面粗さの関係

本節では，これまでの検討より求まった表面粗さと導電率の変化を実際にアンテナ問題に適応しその放射効率の変化を FDTD 法シミュレーションにより求める．FDTD 法の高速化に関しては東北大学サイバーサイエンスセンターの協力を得て実施した．ループアンテナのモデルを図 3 に示
す．1 GHz の解析においては図 3 に示すように，アンテナの幅を 7 mm，一辺の長さを 75 mmとしている．解析は 1 GHz，5 GHz，10 GHz としているので，すべての周波数においてアンテナの全長が 1 波長となるように，周波数によってアンテナの大きさを変化させている．この正方形ループアンテナに正弦波を入力し，FDTD 法により放射効率を計算した．FDTD 法における放射効率の求め方は，入力電力をアンテナの給電電流 I(t)と電圧 V0(t)より下記(4)式で，また放射電力 Prad
を(5)式で，そして放射効率を(6)式で求めている．

Pin 

T

1
V ( t )I( t )dt (4)T 0 0Prad 2

Prad
1 T 2  D( ,  )(6)sin dddt (5)  
Pin
T 0 0 0Z0 

ここで，T は励振波の周期，Prad の積分範囲はアンテナを囲む閉曲面である．また，FDTD 法のセルサイズは 1 GHz の解析では 0.1 mm と細かく設定している．
アンテナ導体の表面が滑らかな場合と粗さの RMS が 0.3 m，1.0 m の場合について，周波数を 1 GHz，5 GHz，10 GHz としたときの放射効率と導電率の関係をそれぞれ図 4 に示す．ここでは，アンテナの材料として銅およびその他の金属を想定している．銅の一般的な導電率は，5.8
×107 S/m 程度であるが，鉄などの導電率が若干下がる媒質ではその数百分の一程度であることから，ここでは，導電率を 5.8×107 S/m から 5.8×104 S/m 程度まで計算を行った．表面粗さ RMSを，測定地に近い 0.3 m，粗さの無い 0 m，および測定地の 3 倍程度粗い 1 m として解析を行った．図 4 に示す放射効率の結果では，表面が粗い場合の効率は，滑らかな場合に比べて最大で 2.4 %ほど減少していることがわかった．

— 46 —

SENAC Vol. 48, No. 4（2015. 10）

図 4： ループアンテナ放射効率と表面粗さ（左：1GHz，中：5GHz，右：10GHz）
続いてアンテナが誘電体上に構成された際のシミュレーションを行った．誘電体の比誘電率は4.3 とした．その解析結果を図 5 に示す．

図 5： 誘電体上ループアンテナ放射効率と表面粗さ（左：1GHz，中：5GHz，右：10GHz）
このモデルにおいても先程と同様に，導電率を 5.8×107 S/m から 5.8×104 S/m 程度まで変化させながら計算を行った．表面粗さについてもさきほどと同様に RMS を測定地に近い 0.3 m，粗さの無い 0 m，および測定地の 3 倍程度粗い 1 m として解析を行った．図 5 に示す放射効率の結果
では，表面が粗い場合の効率は，滑らかな場合に比べて最大で 2.1%ほど減少していることがわかった．
以上の検討より，良導体を用いたアンテナにおいても，その表面粗さにより放射効率が低下す
ることがわかった．しかし，その低下量は数%にとどまっており，今後は表面粗さを実際にモデル化した解析等が必要であると考えられる．

５.

まとめ

本研究では，アンテナに実際に用いられる両面銅箔基板の表面粗さを実測し，その粗さを基に，表面粗さを導電率の変化として FDTD シミュレーションに組み込み，アンテナ効率の検討を行った．
検討には FDTD 法を用いたシミュレーションにより実施した．解析モデルはループアンテナ単体と誘電体基板上のループアンテナとした．しかし，その低下量は数%にとどまっており今後さらなる検討が必要であると考えられる．今後は，実験や数マイクロメートル以下の表面粗さを実際にモデル化した解析等による検討が必要と考えられる．

アンテナの放射効率低下メカニズムに関する一検討— 47 —

謝辞
本研究は，東北大学サイバーサイエンスセンターのスーパーコンピュータを利用することで実現することができた．また，研究にあたっては同センター関係各位に有益なご指導とご協力をいただいた．
参考文献[1] 山本雄太, 宇野 亨, 有馬卓司, 何 一偉 "FDTD 法におけるアンテナ放射効率計算法に関する検討", 電子情報通信学会技術研究報告，AP2011-71，pp.37-42[2] 山本雄太，宇野 亨，有馬卓司” FDTD 法における遠方界を用いたアンテナ放射効率計算法に関する研究 “，2012 電子情報通信学会全国大会論文集，B-1-175, p.175[3] Naoko Yamada, Kentaro Asami, Takuji Arima, and Toru Uno, “Comparison of MoM and
FDTD Calculations of Impedance and Radiation Efficiency for Lossy Plannar Antenna，”
2005 International Symposium on Antennas and Propagation.
[4] 浅見 賢太郎，宇野 亨，有馬 卓司，"FDTD 法による板状アンテナ効率計算の問題点とその一解決法，" 電子情報通信学会総合大会，2005.3[5] Hammerstad and Jensen, “Accurate Models for Microstrip Computer-Aided Design”,
IEEE MTT-S Digest, pp407-409, May 1980
[6] 宇野亨，FDTD 法による電磁界およびアンテナ解析，コロナ社，1998.

— 48 —

SENAC Vol. 48, No. 4（2015. 10）［研究成果］

ベクトルスーパーコンピュータ SX-ACE
による暑熱環境下体温上昇の高速解析
西尾 渉
浅野 陽平 平田 晃正
名古屋工業大学 大学院工学研究科

あらまし 熱ストレスの人体へのリスク評価には，体温上昇，発汗の時間変化を把握することが重要となり，それには熱解析と温熱生理反応を融合した定式化を用いる必要がある．しかしながら，これまで名古屋工業大学の研究グループが開発してきた体温上昇解析コードを用いて， 2
時間の暑熱環境におけるばく露状況を解析するには，通常のワークステーションを用いた場合 7
時間程度必要となる．上記のばく露解析を 20 秒程度で完了することが可能となれば，個体差によるばらつきなど様々な環境を模擬した解析を即時的に処理，普及啓発などへの応用展開が期待できる．本稿では，ベクトルスーパーコンピュータ SX-ACE に上記解析コードを実装し，解析時間の高速化について検討するとともに，解析例を示す．

1. まえがき
暑熱環境などの熱ストレスに対する人体安全性を評価するには，体内温度上昇，発汗の時間変化を解析，検討することが重要となる．名古屋工業大学の研究グループでは，部位ごとの血液温度差を考慮し，かつ温熱生理を考慮した混成熱解析手法を開発してきた[1]．この手法を精巧な数
値人体モデルに適用し，暑熱環境に配置した場合の体内温度上昇，発汗量の解析を行い，測定値との比較によりその有効性を確認している[2]．しかしながら，物理計算に加え，システムバイオ
ロジー，具体的には発汗など体温調節反応を加味した複合的な解析を実施する場合，ベクトル化には適しておらず解析負荷も小さくない．例えば，2 時間の暑熱環境におけるばく露状況を解析するには，通常のワークステーション（CPU ：Intel® Xeon® W5590 @3.33GHz , 4 コア×2）において 7 時間程度必要である．たとえば，このばく露状況の解析を 20 秒程度で実行可能となれば，現実的なばく露環境下において，年齢や性別など複数の人体モデル，さらには熱調整機能の相違などの個人差を考慮した，様々なパターンのデータを即時的に処理することができる．得られたデータを総合的に評価することで，熱ストレスに関する準リアルタイムでのリスク評価が可能となる．本稿では，4 コアで構成されるノードを多数接続したベクトルスーパーコンピュータ
SX-ACE[3]へのコードの実装，および並列実行環境による計算時間の加速率について検討した結果を紹介するとともに，解析例を紹介する．

2. 解析手法
2.1 解剖学的数値人体モデル
解 析 に は ， 情 報 通 信 研 究 機 構 (NICT: National Institute of Information and Communication
Technology )開発の成人男性に相当する日本人の解剖学的数値人体モデルを使用した[4]．これらのモデルは皮膚，筋肉，脂肪，骨，骨髄液，心臓，レンズ，血管など 51 種類の組織から構成されており， 2mm の分解能を有している．また，ハイリスク群である 3 歳児相当のモデルについては[4]に基づき開発された[5]を用いた．
本研究では従来の血液温度計算を改良し，部位ごとに血液温度を考慮出来るようにするため，それぞれのモデルを部位ごとに分割して用いることにした．なお，頭部と胴体部における血流量ベクトルスーパーコンピュータ SX-ACE による暑熱環境下体温上昇の高速解析— 49 —

は十分に大きいため，2 つの部位の血液温度は同一なものと仮定して扱う．数値人体モデルの概観と部位分割の定義を図 1 の(a)，(b)にそれぞれ示す．

(a)

(b)

図 1. 数値人体モデルの外観(a)と部位分割(b)の定義図

2.2 温度上昇解析手法
2.2.1 生体熱輸送方程式
生体熱輸送方程式を計算することにより，温度上昇の数値計算を行う．一般に，生体内の温度上昇計算は，組織内および組織間の熱伝導による熱移動，血流による熱冷却，皮膚表面から外気への熱伝達を考慮することで導出され，体内座標を表す位置ベクトルを r とすると，次式で表される[6]．
T (r, t )C(r )   (r )   (K (r )T )  A(r, t )  B(r, t )  (T (r, t )  Tb,m (t ))(1)t
ここで，C は組織の比熱[J/kg·oC]， は密度[kg/m3]，T は組織の温度[oC]，K は熱伝導率[W/m·oC]，
A は代謝熱[W/m3]， B は血流定数[W/m3·oC]， Tb は血液温度[oC]である．また，境界条件はモデルの体表面および肺内部に対して，次式が適用される．
T (r , t ) K (r ) H (r )  (Ts (r , t )  Te (t ))  40.6  ( SW (r , Ts (r , t ))  PI ) / Sn(2)

ここで， H は熱伝達率[W/m·oC]， TS は外気に接する組織の温度[oC]， Ta は外気の温度[oC]， SW は汗に関する係数[W/m2]，PI は不感蒸散量，n はその表面の法線方向の座標，S は人体の表面積[m2]である．また，40.6 は変換係数である．各種熱定数は，文献[1]に示したものを用いた．式(1)中の血液温度 Tb,m に関しては文献[7]を拡張し，Tb ,m (t )  Tb 0  

Qbtot (t )  Qb ,n (t )C b   b  Vb

dt

(3)

Qbtot (t )   B (t )(Tb (t )  T ( r , t )) dVにより部位間の血液による熱交換を模擬し，部位ごとに血液温度の時間変化を計算することとした．ここで，m(=1,…,5)は分割した部位，n(=1,…,4)は部位間の隣接部分を表している．また，Tb0
は血液の初期温度[oC]，Qbtot は当該部位で血液が獲得する熱量，Cb は血液の比熱[W/ oC]，  b は血液の密度[kg/m3]，Vb は血液の体積[m3]である．分割した部位間で血液により交換される熱量は文— 50 —

SENAC Vol. 48, No. 4（2015. 10）

献[8]より，Qb,n (t )   hx (Tb _ trunk (t )  Tb_limb(t ))dV

(4)で与えられる．ここで，hx は熱交換係数[W/ oC]である．なお，前述の通り，頭部と胴体は血流量が多いため両部位の血液温度は同一であると仮定した[8]．
2.2.2 発汗のモデル化
汗に関する係数 SW は，文献[8]に基づき，SW ( r , t )   (r )[{11 tanh( 11ΔTS  10 )  10 }ΔTS
 { 21 tanh(  21ΔTH   20 )   20 }ΔTH ](5)で与える．式(4)の  TS ，TH は，それぞれ皮膚と視床下部の平均温度上昇である． および  は，
文献[7][8]で示したように，10 =1.20，11 =0.80，10 =0.19，11 =0.59， 20 =6.30， 21 =5.70， 20 =1.03，
 21 =1.98 とした．また，χ は，部位に依存する発汗係数である．
一般に，皮膚表面の蒸発熱量 E は次式で表される[8]．

E  2 .2  hc wf pcl ( Ps   A PA )(6)

ここで，hc は対流性熱伝達係数，w は皮膚の濡れ率，fpcl は被服の透湿率，PS,PA はそれぞれ皮膚
表面および外気の飽和水蒸気圧[mmHg]， A は相対湿度である．皮膚の濡れ率 w は，
“皮膚全体が汗で均一に覆われて，すべて蒸発すると仮定した場合の蒸発熱量（w=1 における E）”に対する“実際の人体表面からの蒸発熱量（式(4) に基づく発汗量）”の比で定義されるものである[8]．発汗が皮膚全面を完全に濡らすほどに達した場合は，それ以上の汗は熱冷却には有効に作用しない．このことから本解析では，式 (5) において，w=1 の値をとるときの蒸発熱量を，熱冷却として働く
発汗量の上限とし，解析に組み込み，実装した．
高齢者のモデル化において，脚部においてのみ[9][10]（一部の文献では背中[10]も含む），最大発汗率の減少が報告されている．そこで，高齢者における発汗率のモデル化に際し，式(4)に対して，以下の変更を加えることとした． i)脚部の発汗率の減少を考慮するために，係数  r  を導入
する．ii)加齢に伴い，熱感度が低下するとの報告に基づき，発汗反応を誘発する閾値として  TS,dec
および  TH,dec を導入した．これらを加味した上で得られる式は以下の通りとなる．

SW (r , t )   (r )  ( r )[{11 tanh( 10 )  10 }ΔTS
 { 21 tanh(  21ΔTH   20 )   20 }ΔTH ]

(7)

ここで，

TS  TS ,0  ΔTS ,dec0ΔTS  TS  (TS ,0  ΔTS ,dec ),TS  TS ,0  ΔTS ,dec0ΔTH  TH  (TH , 0  ΔTH ,dec ),TH  TH , 0  ΔTH ,dec
TH  TH ,0  ΔTH ,dec

(8)

(9)

である．なお，ここでは示さないが，上記の発汗率のモデル化を導入したうえで，暑熱環境下における高齢者の体内温度上昇および発汗量を解析した結果，平均皮膚温度および体内深部温度上昇において文献[10]での実測値とよい一致が得られている．
2.2.3 血流量のモデル化
成人の皮膚における血流定数 B は，

ベクトルスーパーコンピュータ SX-ACE による暑熱環境下体温上昇の高速解析

B(r )  [ B0 (r )  FHB (TH  TH 0 )  FSB  TS ]  2— 51 —(T ( r )T0 ( r ))6(10)で表される[11]．FHB , FSB はそれぞれ 17,500 W/(m3·oC2)， 1100 W/(m3·oC2)という値を有する．また，T (r ) ，T0 (r ) は，皮膚の微小組織における温度と基準温度である．さらに， TS は皮膚の平均温度上昇である．なお，本研究では FDTD 法の時間ステップ Δt および空間差分⊿を，それぞれ 2.0[s]，2.0[mm](8 ヶ月幼児では空間差分 Δ を 1.7[mm])とした．また，皮膚以外の血流定数 B は，

B((r ),T (r ))  B0 (r )T (r )  39o C(11)で表される[6]．式中の T は組織の温度であり，B0 は組織の基準となる血流量であるため，皮膚以外の血流定数 B は一定であると取り扱う．

2.3 解析のフローチャート
上記 2.2 の計算を含む体温上昇解析のフローチャートを図 2 に示す．図より，解析を実行すると，始めに 28 oC 環境下における熱定常状態の初期温度分布を読み込む．次に，人体モデルや発汗モデルの設定，また外気温などの環境設定を入力する．それから，2.2.1 で示した生体熱輸送方
程式による温度上昇を計算，またその計算を行う際，2.2.2，2.2.3 で示した発汗率の変化や血流量の変化など，熱調整機能を組み込んだ計算を行う．このフローチャートに従って，解析を行った．
Start the program

Initialize calculation to obtainthermoneutral state at 28C

Input parameters
Input model data

Loop through Time

Calculate temperature, Eqs. (1), (2)TB,m Eqs. (3), (4)Update thermoregulatory response• Sweating rate, Eqs. (5) – (11)• Blood perfusion, Eq. (12)Thermal conductivity, Specific Heat
Blood perfusion rate

Input sweating model
High, Normal, Low

Setup ambient condition
If Time= Simulation timeNo

Temperature, Humidity

Yes
End simulation (Store results)
temperature distribution, total sweating

図 2. 体温上昇解析のフローチャート
2.4 解析プログラムの最適化
2.4.1 ベクトル化
今回解析に用いたスーパーコンピュータ SX-ACE はベクトル処理を活用すると，高い性能を得ることができる．このことから，ベクトル計算に適したコードを作成する必要がある．ベクトル計算に適したループの例を図 3 に示す．言語は fortran90 にて記載する．図に示すように，最内側

— 52 —

SENAC Vol. 48, No. 4（2015. 10）のループ長を最長にすることで，ノード当たりの計算要素数を増加させることが可能となる．
このように，さまざまなループレベル最適化[12]を施すことで，ベクトル化率とノードあたりの性能を向上させた．

MODELX=160
MODELY=320
MODELZ=866

MODELX=866
MODELY=320
MODELZ=160

DO K = 1,MODELZDO J = 1, MODELY
DO I = 1, MODELXTEMP(I,J,K)
END DOEND DOEND DODO I = 1,MODELZDO J = 1, MODELYDO K = 1, MODELXTEMP(K,J,I)
END DOEND DOEND DO

図 3. ベクトル化に適したループ
2.4.2 並列化
コードを並列化する際，MPI(Message Passing Interface)を用いて並列化を行った．MPI とは，分散メモリ間のデータ処理分割やデータ転送を担うメッセージパッシングの標準規格であり，コードから呼び出すサブプログラムのライブラリである．MPI を用いて並列化を行った理由として，並列化手法には，タスク並列用の自動並列化や OpenMP(Open Multi-Processing)，データ並列用のMPI などがある．今回，解析コードを並列化する上で，タスク並列と MPI を併用したハイブリッ
ド MPI，MPI のみを用いたフラット MPI をそれぞれ実装して，実行時間を比較した．その結果，フラット MPI の方が効率よく計算を高速化できたため，フラット MPI を用いた．
具体的に行った並列化処理としては，人体モデルを Z 軸方向各平面に分割し，それぞれの 1 平面
を 1 プロセスとして割り当てた．

3. 解析結果異なる並列プロセス数における計算時間を表 1 に示す．表 1 より算出したプロセス数の増加に対する計算時間の加速率を図 4 に示す．図より，並列プロセス数の増加に伴って並列処理効率が顕著に低下することはなく，加速率は上昇している．しかしながら，理論値 64 倍に対して，並列
プロセス数が 64 の場合に得られた加速率は，34.6 倍となった．これは主に演算量のインバランス
による並列化率の低下や，ノード間通信時のデータ転送量が増加する際に発生する遅延によるものと考えられる．

表 1. 異なる並列プロセス数における計算時間
並列プロセス数

1

4

8

16

32

64

計算時間(sec)

2329

721

391

213

124

67

— 53 —

ベクトルスーパーコンピュータ SX-ACE による暑熱環境下体温上昇の高速解析

70
Speed Up

60

Theoretical Value

Speed Up

50

Approxiamte Curve

40
30
20
10
0

0

8

16

24
32
40
Number of Process

48

56

64

図 4. 並列プロセス数の増加に対する計算時間の加速率

ベクトル化，並列化を行ったコードを SX-ACE に実装し，暑熱ばく露における体温上昇を解析した例を図 5 に示す．図は 3 時間の暑熱ばく露(37 oC)を行った際の体表面温度を表している．左から 3 歳児，22 歳相当の成人男性，65 歳相当の成人男性，75 歳相当の成人男性，成人女性となっている．
図より，他のモデルに比べて，75 歳相当の成人男性の体表面温度が最も大きく上昇していることが確認できる．これは，発汗機能の相違によるもので，高齢者は若年者よりも発汗量が少なく，発汗の開始が遅れるためである．

図 5. 暑熱ばく露における表面温度上昇の解析結果

— 54 —

SENAC Vol. 48, No. 4（2015. 10）

4. まとめ本研究では，SX-ACE による体温上昇解析の高速化を目標とし，コードのベクトル化，並列化を行い，異なる並列プロセス数における計算時間を算出した．算出した結果，プロセス数における加速率が理論値 64 倍に対して，計算時間の加速率は 34.6 倍となった．これは，逐次処理による並列化率の低下や，データ転送時の遅延によるものと考えられる．また，解析例として暑熱ばく露における年齢や性別など 5 体の人体モデルに対して体表面温度上昇解析を行った．他のモデルに比べて，75 歳相当の成人男性が最も大きく上昇するという結果が得られた．これは，高齢者の発汗反応の遅延および発汗量に起因するものである．
現在，図 5 に示した体表面温度上昇の解析例などの体温上昇解析に必要な計算時間は，ばく露
時間によって差はあるものの，目標である 20 秒に達していない．目標を達成するためには更なる並列化率の向上が今後の課題となる．

謝辞
本研究は，東北大学サイバーサイエンスセンターのスーパーコンピュータ SX-ACE を利用することで実現でき，計算コードの高速化および MPI 並列化にあたっては，同センターの佐々木技術
職員ならびに山下技術専門職員に多大なるご協力をいただきました．この場を借りて深謝致します．
参考文献[1] A. Hirata, T. Asano, O. Fujiwara : “FDTD analysis of body-core temperature elevation in children
and adults for whole-body exposure”, Phys. Med. Biol., vol.53. Pp.5223-5238, (2008).
[2] A. Hirata, T. Nomura and I. Laakso : “Computational estimation of decline in sweating in the elderly
from measured body temperatures and sweating for passive heat exposure,” Physiol. Meas., vol. 33,
pp.N51-60, (2012).
[3] R. Egawa, S. Momose, K. Komatsu, Y. Isobe, H. Takizawa, A. Musa, H. Kobayashi : “Early
Evaluation of the SX-ACE Processor”, Proc. Int’l Conf. High Pef. Comput. Storage & Anal., no.14,(2014).
[4] T. Nagaoka, S. Watanabe, K. Sakurai, E. Kunieda, M. Taki, Y. Yamanaka : “Development of realistic
high-resolution whole-body voxel models of Japanese adults male and female of average high andweight, and application of models to radio-frequency electromagnetic-field dosimetry”, Phys. Med.
Biol., vol.49, pp.1-15, (2004).
[5] T. Nagaoka, E. Kunieda, and S. Watanabe : “Proportion-corrected scaled voxel models for Japanesechildren and their application to the numerical dosimetry of specific absorption rate for frequencies
from 30 MHz to 3 GHz”, Phys. Med. Biol., vol.53, no.23, pp.6695-6712, (2008).
[6] H. H. Pennes : “Analysis of tissue and arterial blood temperatures in resting forearm”, J.Appl.Physiol.,vol.1, pp.93-122, (1948).
[7] A. Hirata, T. Asano and O. Fujiwara : “FDTD analysis of human body-core temperature elevation due
to RF far-field energy prescribed in the ICNIRP guidelines,” Phys. Med. Biol., vol.52, pp.5013-23,(2007).
[8] D. Fiala, K. J. Lomas, and M. Stohrer : “Computer prediction of human thermoregulation andtemperature responses to a wide range of environmental conditions,” Int J Biometeorol, vol.45,
pp.143-159, (2001).
[9] A. Dufour and V. Candas : “Ageing and thermal responses during passive heat exposure: sweating
and sensory aspects,” Eur. J. Applied Physiol., vol. 100, pp.19-26, (2007).
[10] Y. Inoue and M. Shibasaki : “Regional differences in age-related decrements of the cutaneous
vascular and sweating responses to passive heating,” Eur. J. Appl. Physiol. & Occup. Physiol., vol.74,
pp.78-84, (1996).
[11] P. Bernardi, M. Cavagnaro, S. Pisa and E. Piuzzi : “Specific absorption rate and temperature elevation
in a subject exposed in the far-field of radio-frequency sources operating in the 10-900-MHz range,”
IEEE Trans. Biomed. Eng., vol.50, pp.295-304, (2003).
[12] 王 建青, 岡部 公起, 曽根 秀昭, 藤原 修 : “ベクトル化・並列化に適した FDTD プログラミ
ングの一手法” ,信学技報.,Vol.103, No.223, pp31-35, (2003).

— 55 —

SENAC Vol. 48, No. 4（2015. 10）

[報 告]
［報 告］

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

— 56 —

SENAC Vol. 48, No. 4（2015. 10）

[報 告]
［報 告］

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

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

— 57 —

SENAC Vol. 48, No. 4（2015. 10）

［報告]
告］
[報
＜プレスリリース報告＞

準リアルタイムで熱中症リスクを評価する技術の開発～個人に応じた熱中症対策を提案する、新たな情報～

平成 27 年 7 月 21 日、以下の内容でプレスリリースをいたしましたのでご報告いたします。
【概要】
東北大学サイバーサイエンスセンター、名古屋工業大学、一般財団法人日本気象協会の共同研
究 グループは、1）乳幼児や高齢者などの個人特性を考慮した熱中症リスク評価のための複合物
理・システムバイオロジー統合シミュレーション技術をスーパーコンピュータに効率的に実装、高速化、2）気象予報データと融合させることによる、個人特性を考慮した 3 時間後の熱中症のリ
スクを 10 分で評価する技術の開発に成功しました。今後予定されている大規模なスポーツ大会や屋外イベント等において、個人属性を考慮した適切な熱中症リスク評価技術を活用することで、今まで以上に場面に応じた発症数の低減に貢献することが期待できます。
【これまでの研究・開発の経緯】(1) 名古屋工業大学の研究は、1 辺 2 ミリメートルの微小立法体から構成され（総計約 800 万点）、50 以上の体内組織を考慮した温熱人体モデルを構築し、熱中症のリスク要因である体内深部温度（以下、体温）上昇および発汗量を解析してきました。具体的には、①太陽光の吸収エネルギ
ー、②環境温度というマルチフィジックス現象による体温上昇、さらにはそれに伴うシステムバ
イオロジー（発汗、血流の変化などの熱調整機能）を統合的に取り扱う技術を有しており、サウナなど特殊な温熱環境まで有効性を確認した唯一の技術を有しています。3 時間後のリスク評価を行うには、通常ワークステーションを用いた場合、解析には 12 時間を要していました。
(2) 東北大学サイバーサイエンスセンターは全国共同利用施設として、高メモリバンド幅と高い
コ ア性能を有するベクトル並列型スーパーコンピュータを整備・運用し、ベクトル処理および並列処理を高度に融合させたメモリ性能制約型アプリケーションの高度化、高速化に関する計算 科
学者（ユーザ）と計算機科学者の学際的共同研究を積極的に推進しています。今回は、この取り組みを通して当該コードにベクトル化、並列化を施すことでスーパーコンピュータ SX-ACE の高い演算性能を引き出すことに成功しています。
(3) 日本気象協会では、これまでにもさまざまなメディアや気象情報提供を通して、広く一般の方に向けた熱中症予防に関する啓発活動を行ってきました。2012 年から推進している熱中症予防啓発プロジェクト『熱中症ゼロへ®』では普段は積極的に予防に関する行動をとらない層に向けても、イベントや協賛企業、各種メディアなど日々多くの人が接する媒体を通して、熱中症予防に呼びかけてきました。また、全国のリアルタイムの気象情報や天気予報を提供する日本気 象協会の天気予報専門サイト「tenki.jp」では、気象予測情報を活用した地点別の熱中症情報の予測や、

— 58 —

SENAC Vol. 48, No. 4（2015. 10）

熱中症予防に関するコラムなどを公開しています。そこで、共同研究グループで得られた知見を
このような実績のある媒体などを通して広く一般の方に発信することで、開発された技術を有効に活用してもらうプラットホームを提供します。さらに、個人属性を考慮した熱中症リスク評価を提供情報として活用することで、より「他人事にしない」熱中症予防啓発 につなげることを目
指します。
【新規性】
今回の共同研究では、名古屋工業大学が開発してきたマルチフィジックスとシステムバイオロ
ジー融合した解析コードを、東北大学が中心となりベクトル化、並列化を同時に施し、スーパーコ ンピュータ SX-ACE に効率的に実装しました。なお、これまでの物理解析とは異なり、ベク
トル化には適していないシステムバイオロジーの部分を効率的に実装するよう工夫しています。
スーパーコンピュータにおいて 64 コアを同時に用いた結果、3 時間後のリスク評価に要する解析
時間は、成人男性で 45 秒（3 歳児は 15 秒）となりました。
この解析手法に、日本気象協会が提供する気象予測データを入力情報として組み込むことにより、例えば、大規模イベントなどにおける熱中症リスク評価シミュレーションが準リアルタイム
で実施できるようになりました。年齢や性別など複数の人体モデル（5 体）、さらには熱調整機能のばらつき（各 3 種類）などの複数のパターンのデータを並列的に処理（64Core で単一のシミュ
レーション、同時に 2 つの解析を実施）した場合、10 分以内に大半の人口をカバーした評価が可能となります（下 記図は、2014 年 7 月 25 日東京での評価例）。東北大学サイバーサイエンスセ
ンターのスーパーコンピュータでは、最大 2048 コアの利用が可能であり、さらなる高速化も可能です。

【今後の展開】 気象予測情報を組み合わせた新しい熱中症に関する情報提供や熱中症予防を呼びかけるコラムなど、日本気象協会が提供する各種情報として活用いたします。

— 59 —

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

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

SX-ACE 1,024 ノード実行の提供開始について（No.203）
スーパーコンピュータ SX-ACE での最大実行ノード数が、これまでの 512 ノードから 1,024 ノードに変更になります。
1,024 ノードを利用した場合、理論演算性能 283TFLOPS、（最大ベクトル性能 262TFLOPS）、最大メ
モリ 60TB での実行が可能です。キューの指定方法等に変更はありません。-b オプションで利用ノード数を指定してください。なお、2 ノード以上を利用した実行の場合、MPI ライブラリの利用が必要です。
257 ノード以上を利用するリクエストの実行時間（経過時間）の制限値は、規定値・最大値ともに
1 ヶ月間（744 時間）です。計算機資源の効率的な利用のため、-l elapstim_req オプションにより必要十分な実行時間の指定をお願いいたします。
最大実行ノード数変更日 ：（バッチリクエストファイル例）#!/bin/csh#PBS -q sx -b 1024
#PBS -l elapstim_req 01:00:00
mpirun -np 4096 -nn 1024 ./a.out

2015 年 9 月 18 日（金）# SX-ACE を 1024 ノード利用する# 実行時間を 1 時間に制限する# 4,096 プロセスで実行する（Flat MPI）

SX-ACE の利用方法については以下のページをご覧ください。
http://www.ss.cc.tohoku.ac.jp/super/guide.html
利用負担金については以下のページをご覧ください。
http://www.ss.cc.tohoku.ac.jp/utilize/academic.html
http://www.ss.cc.tohoku.ac.jp/utilize/business.html（共同利用支援係）

— 60 —

SENAC Vol. 48, No. 4（2015. 10）

非線形構造解析汎用プログラム MSC.Marc/Mentat のバージョンアップについて（No.203）
非線形構造解析汎用プログラム「MSC.Marc/Mentat」のバージョンアップを行いますのでお知らせ
いたします。
MSC.Marc は有限要素法による非線形汎用構造解析プログラムです。世界中で広く利用され最も評価を受けているプログラムの一つです。
Mentat2014 では新 GUI を採用しています。従来の Classic GUI は mentat.classic というコマンド
でご利用できます。
バージョン名：バージョンアップ日 ：サービスホスト：起動コマンド：2014.22015 年 9 月 18 日（金）front.cc.tohoku.ac.jp (並列コンピュータ)
run_marc，mentat, mentat.classic

新機能の概要、機能の詳細、実行方法については開発元 HP、およびセンターHP をご覧ください。
開発元 HP
http://www.mscsoftware.co.jp/products/marc/index.php
https://simcompanion.mscsoftware.com/infocenter/index?page=home&locale=ja_JP
センターHP
http://www.ss.cc.tohoku.ac.jp/application/marc_mentat.html(共同利用支援係)

数式処理プログラム Mathematica のバージョンアップについて（No.203）
数式処理プログラム「Mathematica」のバージョンアップを行いますのでお知らせいたします。
Mathematica は Stephen Wolfram によって作られた、プログラミング言語を備えた数式処理システ
ムです。Mathematica の機能は、数値計算、記号計算、グラフィックスという 3 つに大別でき、この3 つが一体となって使いやすいインターフェイスを提供しています。
バージョン名:バージョンアップ日 :サービスホスト:起動コマンド::Mathematica 10.2
2015 年 9 月 18 日（金）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（共同利用支援係）

Web 版大規模科学計算システムニュースより— 61 —

科学技術計算言語 MATLAB のバージョンアップについて（No.203）
科学技術計算言語「MATLAB」のバージョンアップを行いますのでお知らせいたします。
MATLAB は、アルゴリズム開発、データの可視化、数値計算を行うための高レベルなテクニカルコ
ンピューティング言語と対話型環境です。 MATLAB® を利用することにより、C、C++、Fortran といった伝統的なプログラミング言語よりも短時間で科学技術計算の問題を解決することが可能です。
バージョン名：バージョンアップ日：
サービスホスト：起動コマンド：：

MATLAB R2015a (Ver.8.5)
2015 年 9 月 18 日（金）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（共同利用支援係）

— 62 —

SENAC Vol. 48, No. 4（2015. 10）

反応経路自動探索プログラム GRRM14 のサービス開始について（No.203）
GRRM は、2002 年に東北大学（教授：大野公一、修士１年：前田理、当時）で制作が開始され、その後開発が進められて、2011 年に GRRM11、2014 年に GRRM14 が発表され、広く利用されるようになりました。
GRRM には、以下のような特長があります。






Gaussian プログラム（g09、g03）などの非経験的量子化学計算に基づいて、各化学式で表される構造や反応経路を自動的に探索します。
平衡構造から出発し、その周囲に存在する反応経路を、ポテンシャルの非調和下方歪みを検出して、系統的に調べ上げる超球面探索アルゴリズムが搭載されており、反応経路自動探索を行うことができます。
励起状態のポテンシャル交差を自動的に調べることができます。
解離した状態から、人工力誘起反応法で、反応経路を効率的に調べることができます。

バージョン名：バージョンアップ日：
サービスホスト：

14.01
2015 年 9 月 18 日（金）front.cc.tohoku.ac.jp（並列コンピュータ）

利用方法についての詳細は NPO 法人 量子化学探索研究所の HP、および当センターのGRRM14 紹介 Web ページをご参照下さい。
NPO 法人 量子化学探索研究所
http://iqce.jp/
センターGRRM14 紹介 Web ページ
http://www.ss.cc.tohoku.ac.jp/application/grrm14.html（共同利用支援係）

— 63 —

SENAC Vol. 48, No. 4（2015. 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)投稿予定の原稿が 15 ページを超す場合は共同利用支援係まで前もってご連絡ください。
(3)初回の校正は、執筆者が行って、誤植の防止をはかるものとします。
(4)原稿の提出先は次のとおりです。
東北大学サイバーサイエンスセンター内 情報部情報基盤課共同利用支援係
e-mail uketuke@cc.tohoku.ac.jp
TEL022-795-3406

— 64 —

SENAC Vol. 48, No. 4（2015. 10）

スタッフ便りスーパーコンピュータ．
．．
．皆さん，良くご存知のように，それはその時々の最速のコンピュータのことを云います．世界最初の電子計算機と云われる ENIAC（1946 年）は加減算が 1 秒間に 5000 回しか出来ない程度の性能でしたが，現在のスーパーコンピュータは 1 秒間に 1017 回の演算ができる性能に迫る勢いです．70 年間で約 10 兆倍の性能をもつまでに急成長を遂げた機械（装置）は他に見当たりません．我々は，この格段に性能向上を達成した稀有な装置を十分に活かしきれているでしょうか？
昨年設置された「高性能計算技術開発（NEC）共同研究部門」は，まさしく高性能計算技術を研究する部門です．昨年 10 月からこの研究部門に客員として着任した私は，これまで国の 2 つのスーパーコ
ンピュータ開発プロジェクトに携わってきました．高度なアプリケーションには，ハードウェア動作や基本ソフトウェアの癖を知ることが重要であるけれども，それをすべて理解することの難しさを実感してきました．
この研究部門での活動を通して，多くの困難な課題が解決できるようなアプリケーションの開発に適ったスーパーコンピュータ利用技術の道筋が出来ればいいと思っています．(M.Y）

7 月に異動してきてから，早いもので 3 ヶ月が過ぎました．センター内で飛び交う聞き慣れない用語に，頭の中はまだまだ混乱していますが，新しい環境に刺激を受け，新しい業務内容に悪戦苦闘しながら充実した日々を過ごしています．
さて，地下鉄東西線開業まで約 2 ヶ月となりました．青葉山駅の見学会がおこなわれたり，青葉山
連絡バスの試運転が始まったりと，開業に向けて盛り上がりを感じます．開業後は，いつもは静かなセンター周辺も少しにぎやかになるのかなと，少々楽しみになってきました．
話は変わりますが，2 歳半を過ぎた娘はおしゃべりが急激に上手になってきました．ずっと止まらな
いおしゃべりを聞いていると，こちらの想像以上に様々なことを考え，理解しているよう…「この小さな頭の中はどうなっているのだろう」と覗いてみたくなります．と同時に，あれもこれも知りたい！
と思う意欲やどんどん新しいことを吸収できる力をとても羨ましく感じています．
そんな彼女に負けないよう，私もセンターや担当する業務に関する知識や技術をどんどん習得し，少しでも皆様のお役にたてるよう精一杯努めたいと思いますので，どうぞよろしくお願いします．（M.K）

サイバーサイエンスセンター・情報部情報基盤課スタッフ転入・転出者のお知らせ[転入]2015.7. 1
2015.7. 1[転出]
2015.7. 1
2015.7. 1

栗原尚志（情報基盤課長補佐）
加藤美久（情報基盤課総務係）

病院施設企画室室長補佐から理学部総務課人事係から田口秀樹（情報基盤課長補佐）
加藤紀恵（情報基盤課総務係）

宮城教育大学財務課副課長へ
病院総務課人事係へ

SENAC 編集部会
小林広明 曽根秀昭 水木敬明 後藤英昭
江川隆輔 佐藤恵美子 高杉佳奈 大泉健治
小野 敏 斉藤くみ子

編集・発行

印

整備中の青葉山新キャンパス

刷

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

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

/;5H

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

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

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

ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢱ

㐃

⥆

㐠

㌿

୪ิࢥࣥࣆ࣮ࣗࢱ

㐃

⥆

㐠

㌿

ྍど໬ᶵჾᐊ

ᖹ᪥ 㹼

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

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

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

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

ࠥ ᫬㛫

ࠥ ᫬㛫

㏻ᖖ

↓ᩱ
ࢹࣂࢵࢢ

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

*%ࣀ࣮ࢻᩘ

฼⏝ࣀ࣮ࢻᩘ

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

฼⏝ࣀ࣮ࢻᩘ

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

㏻ᖖ

ࠥ

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

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

࡞ࡋ

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

*%

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