ISSN 2436-0066

Vol.55 No.2

2022-4኱つᶍ⛉Ꮫィ⟬ࢩࢫࢸ࣒㛵㐃᱌ෆ㸺኱つᶍ⛉Ꮫィ⟬ࢩࢫࢸ࣒㛵㐃ᴗົࡣࠊࢧ࢖ࣂ࣮ࢧ࢖࢚ࣥࢫࢭࣥࢱ࣮ᮏ㤋ෆࡢ᝟ሗ㒊᝟ሗᇶ┙ㄢࡀᢸᙜࡋ࡚࠸ࡲࡍࠋ㸼https://www.ss.cc.tohoku.ac.jp/
*

㝵

ಀ䞉ᐊྡ

฼⏝┦ㄯᐊ

㟁ヰ␒ྕ䠄ෆ⥺䠅e-mail022-795-6153(6153)
sodan@cc.tohoku.ac.jp┦ㄯဨ୙ᅾ᫬

୺䛺䝃䞊䝡䝇ෆᐜィ⟬ᶵ฼⏝඲⯡䛻㛵䛩䜛┦ㄯ䝃䞊䝡䝇᫬㛫
ᖹ ᪥
8:30䡚17:15

኱ุ䝥䝸䞁䝍䚸฼⏝⪅➃ᮎ➼䛾฼⏝

8:30䡚21:00

022-795-3406(3406)
฼⏝⪅ㄯヰᐊ

(3444)⮬㈍ᶵ

8:30䡚21:00

ᒎ ♧ ᐊ* (ศᩓ*ぢᏛᕼᮃ䛾᪉䛿ඹྠ฼⏝ᨭ᥼ಀ䜎䛷䛤䡶䢙䢇䢛䡩䡬䡼༤≀㤋)*㐃⤡䛟䛰䛥䛔䚹

Ṕ௦䛾኱ᆺィ⟬ᶵ➼䛾ᒎ♧

9䠖00䡚17䠖00⥲ົಀ022-795-3407(3407)
cc-som@grp.tohoku.ac.jp⥲ົ䛻㛵䛩䜛䛣䛸

8:30䡚17:15

఍ィಀ

022-795-3405(3405)
cc-kaikei䠜䡃rp.tohoku.ac.jp఍ィ 䛻㛵䛩䜛䛣䛸䚸㈇ ᢸ㔠 䛾ㄳồ䛻㛵䛩䜛䛣䛸

8:30䡚17:15

ඹྠ฼⏝ᨭ᥼ಀ(ཷ ௜)022-795-3406䠄3406䠅
022-795-6251(6251)
uketuke@cc.tohoku.ac.jp฼⏝ᡭ⥆䛝䚸฼⏝┦ㄯ䚸ㅮ⩦఍䚸䝷䜲䝤䝷䝸䚸ぢᏛ䚸䜰䝥䝸䜿䞊䝅䝵䞁䛻㛵䛩䜛䛣䛸

8:30䡚17:15

ඹྠ◊✲ᨭ᥼ಀ

022-795-6252(6252)
rs-sec@cc.tohoku.ac.jpඹྠ◊✲䚸ィ⟬ᶵ䝅䝇䝔䝮䛻㛵䛩䜛䛣䛸

8:30䡚17:15

䝛䝑䝖䝽䞊䜽ಀ022-795-6253(6253)
i-network@grp.tohoku.ac.jp䝛䝑䝖䝽䞊䜽䛻㛵䛩䜛䛣䛸

8:30䡚17:15

ᅄ
㝵

◊✲㛤Ⓨ㒊

022-795-6095(6095)

஬
㝵

➃ᮎᶵᐊ

(3445)

PC ➃ᮎᶵ䠄X ➃ᮎ䠅

8:30䡚17:15*㸦㸧ෆࡣᮾ໭኱Ꮫෆࡢࡳࡢෆ⥺㟁ヰ␒ྕ࡛ࡍࠋ㟷ⴥᒣ࣭ᕝෆᆅ༊௨እ࠿ࡽࡣ㢌࡟ 92 ࢆຍ࠼ࡲࡍࠋ

୍
㝵

୕
㝵

ᮏㄅ䛾ྡ๓䛂SENAC䛃䛾⏤᮶᫛࿴ 33 ᖺ䛻ᮾ໭ᆅ༊䛾᭱ึ䛾㟁Ꮚィ⟬ᶵ䛸䛧䛶䚸ᮾ໭኱Ꮫ㟁Ẽ㏻ಙ◊✲ᡤ䛻䛚䛔䛶᏶ᡂ䛥䜜䛯䝟䝷䝯䝖䝻䞁ᘧィ⟬ᶵ䛾ྡ๓䛷 SENAC-1(SENdai Automatic Computer-1)䛛䜙䛸䛳䛶࿨ྡ䛥䜜䛯䚹SENAC Vol. 55, No. 2（2022. 4）― 1 ―［共同研究成果］

更新され たベクト ル機での 日本域を 対象とし た長期領 域再解析 の実行
福井真, 曽我大輝, 小原涼太, 臼井健, 瀧口海人, 小野佳祐, 廣瀬大河, 松島沙苗,
白川栄一, 伊藤純至, 岩崎俊樹, 山崎剛(東北大学大学院理学研究科),
斉藤和雄(気象業務支援センター/東京大学大気海洋研究所/気象研究所), 瀬古弘(気象研究所)
日本域を対象とした過去の気象場を再現するアンサンブルカルマンフィルタと非静力学数値気象モデルによる領域再解析を更新されたベクトル機(SX-Aurora TSUBASA)上で実行した。2018
年 7 月の 1 か月間を対象として、旧機(SX-ACE)による実行結果と比較したところ、正常に動作することが確認できた。旧機と比べ、ノード数 1/32 で 12％程度高速に実行できた。今後 2 年程
度で約 60 年の長期間の領域再解析が実行可能な環境が整った。今後は、近年の日本における地域
気候の変動の把握を目指し、長期領域再解析の実行を進めていく。また、更なる安定した領域再解析のために、サンプリングエラー軽減が課題の一つである。

1.

はじめに

我々は、日本域を対象とした気象場を長期間にわたり高解像度且つ精度良く再現することを目標に、長期領域再解析に取り組んでいる。長期領域再解析を通じて、日本域における近年の気候変動の実態の把握、過去のメソスケール極端気象の分析、防災や農業、エネルギー、水資源管理などでの気象/気候情報の高度利用に向けた基礎資料の提供、数値気象予報/気候予測システムの改良に資する知見の獲得などを目指す。
地上観測に加え、ラジオゾンデ高層観測網が全球的に整備されてから 60 年以上が経過した。その間に、数値気象予報モデルやデータ同化手法が、計算機の発展によって着実に向上してきた[1]。
こうした背景によって、蓄積された観測と高度化された数値予報技術を組み合わせ、観測できな
い量も含めた全球の大気の尤もらしい状態を長期にわたって推定する長期全球再解析が、世界の主要気象機関によって実施されるようになり、研究から現業利用まで、気象・気候に関わる様々な用途で活用されている。
しかし、複雑な地形の影響を受ける地域気候やメソスケール気象を対象とするには、全球再解
析では最新のもの（e.g. ERA5[2]）でも水平解像度が 30–40km 格子相当であるため不十分である。
これまで、解像度の問題に対応するために、高解像度の領域気候モデルを用いた力学的ダウンス
ケールが広く実施されてきた。しかし、積分時間が短いとスピンアップの問題が生じ、逆に積分時間が長いとモデル領域内部の予測値が時間と共に境界値として与えた場から乖離してしまう問
題を抱えている。モデルによる時間積分だけではなく領域内部の観測の同化も行う領域再解析によって、上述の問題が軽減できることが期待される。
これまでに、我々は日本域を対象とした水平格子間隔 5km の長期領域再解析のためのシステムを構築し、領域モデルのみによる力学的ダウンスケールや全球再解析に対する有効性を示してきた[3]。しかし、計算機の性能の限界によって、地域の気候変動を調査するために必要な長期間（数十年以上）を対象とした領域再解析の実行は困難であった。2021 年 10 月に東北大学サイバーサ
イエンスセンター大規模科学計算システムが更新され性能が向上した。この更新に伴い、長期領域再解析の実行が可能になることが期待できる。そこで、本稿では、これまで構築してきた長期領域再解析システムの概要を紹介し、2018 年 7 月を対象とした領域再解析を新機で実行した結果を旧機で実行された結果と比較しながら記述する。

2.

領域再解析システム

用いた領域再解析システムは、全球再解析 JRA-55[4]を境界値に与えながら気象庁非静力学モデル(NHM[5])によって時間発展を求め、局所アンサンブル変換カルマンフィルタ(LETKF[6])によって観測データを同化する解析予報サイクルを実行する。水平格子間隔は、 NHM-LETKF を一

― 2 ―

SENAC Vol. 55, No. （2 2022. 4）

方向ネストさせることによって、25km で東アジア域を、5km で日本域を対象とするよう設定した。同化する観測は、長期間入手可能な地上気圧観測、ラジオゾンデ高層観測、及び海上にある
台風の中心位置に限定し、近年になって観測網が整備されたレーダや衛星による観測は利用しない。これは、気候変動シグナルを抽出するには、数値モデル及びデータ同化手法に加え同化する
観測も再解析の期間で一貫させることが有効であるからである[7]。本システムでは、時間発展を解く NHM では 1 本の標準ランと 30 本の摂動ランを実行し、LETKF では標準ランの結果を第一推定値として用い、その誤差共分散を摂動ランから評価するようにした。第一推定値として、
LETKF で一般に用いられる摂動ランのアンサンブル平均は利用しなかったのは、本システムでは解析場が平滑化され現実的な変動が再現できなくなるのを回避するためである。側面境界摂動について、25km 格子の NHM-LETKF の摂動ランには、JRA-55 を EOF 解析して得られた上位 30 モードをランダムに与えるようにした。実行は、６時間を１サイクルとして、開始時刻を 25km 格
子の NHM-LETKF は 6 月 20 日 12UTC、5km 格子の NHM-LETKF は 6 月 29 日 12UTC からとした。
6 月中はスピンアップ期間として、比較には 7 月の再解析結果を用いた。

3.

結果

新機と旧機を用いて計算した領域再解析の再現結果を比較する。図１は、2018 年 7 月の海面較
正気圧について、JRA-55、新機及び旧機の領域再解析間の RMSD を示している。JRA-55 との RMSD
は、新機と旧機の領域再解析で共に期間の平均が 0.7hPa で同程度だったのに対して、新機と旧機の領域再解析間の違いは、平均 0.4hPa と JRA-55 の RMSD に比べ半分程度に収まった。新機においても概ね旧機での計算と同等の総観場を再現できており、領域再解析システムが新機においても旧機と同様の動作をすることが確認できた。なお、7 月 22 日周辺で RMSD が大きくなった期間があったが、この期間は領域内に解析された台風の位置にずれがあったことに起因していた。
台風は、それ自体がシャープな構造を有するために位置ずれが RMSD に反映されやすい。さらに、観測の疎らな日本の南方海上から接近するため場の不確実性が大きいことを示唆する。また、本システムのアンサンブルメンバー数は 30 と限られており、サンプリングエラーの影響を受けやすい。本システムでは、境界摂動をランダムに与えているので、新機の実験と旧機の実験とで摂動が異なり、再解析の結果に違いをもたらしうる。このことは、特にアンサンブル数が少ない場合において、成長する誤差モードを効率的に捉え領域再解析のより安定した再現性を確保するために、側面境界摂動の与え方が重要である可能性を示唆する。
極端豪雨の事例として、平成 30 年 7 月豪雨の新機における再現結果を示す。図 2 は、2018 年 7
月 5 日からの 72 時間積算雨量を示したものである。地上レーダによる降水強度の観測を地上での降水量の直接観測で補正した解析雨量によると、72 時間積算雨量 200mm を超える強雨域が九州から中部山岳の東側にかけて伸びており、九州北部や四国の南側、岐阜県北部などでは 500mm
を超える地域が局所的に存在した(図 2a)。水平格子間隔約 50km 相当の解像度の全球再解析 JRA-55
では、西日本に強雨域が広がっていた様子は再現できているが、局地的な分布は十分捉えることができていない(図 2b)。一方、水平格子間隔 5km の領域再解析では、必ずしも観測と一致するわけではないが、局地的な強雨域も含めて概ね再現することができていた。領域再解析では、高解像度モデルを用いることで、細かい地形の影響やメソスケールの降水系に伴う強い降水の再現性向上につながったと考えられる。なお、これらの特徴は、旧機での本事例の再解析結果[8]と整合的であった。更新されたベクトル機上でも、本領域再解析システムが想定の動作することが確認できた。
計算に要した資源量について、計算機の更新前後で比較する。1 サイクル (6 時間分の再解析) の実行に要した時間は旧機では 32 ノードで平均約 1.9 時間に対して、新機では 1 ノード(8VE)で平均約 1.6 時間となった。計算機の更新に伴い、1/32 の資源量で 12%程度高速に実行可能となった。
これは、1 ノード当たりの演算性能の向上が主要因であるが、1 ノードに搭載されるメモリが増加したことでノード間通信をせずに計算が可能になったことも一部寄与している。1 サイクルの時

更新されたベクトル機での日本域を対象とした長期領域再解析の実行― 3 ―間積分及び解析が少ない資源で且つ高速に実行可能となったことで、複数の計算ストリームを並行した実行も現実的となった。これによって、近年の気候変動を捉えるための長期領域再解析実行が今後数年間で可能な環境が整った。

4.

まとめ

日本域を対象に 2018 年 7 月の 1 か月分の領域再解析を更新されたベクトル機上で実行し、旧機での実行結果と比較し、正常に動作することが確認できた。旧機と比べノード数 1/32 で 12％程
度高速に実行できた。今後 2 年程度で約 60 年の長期間の領域再解析が実行可能な環境が整った。
近年の日本における地域気候の変動の把握を目指し、長期領域再解析の実行を進めることが今後の課題である。
また、より安定した再現性に向けた領域再解析システムの開発も重要な課題である。一つとして、誤差共分散の推定の際のサンプリングエラーを軽減がある。そのために、アンサンブル数を増やすこと、場に依存した側面境界摂動の導入が挙げられる。サンプリングエラーによる影響は抑えられる一方、力学的整合性を崩してしまう局所化も緩和できれば、相乗効果として再現性向上が期待できる。これらは、更なる計算資源を必要となるため、計算機の一層の性能向上も併せて望まれる。

謝辞
東北大学サイバーサイエンスセンターのスーパーコンピュータの令和 2 年度（第２期）共同研
究の一般課題として本研究は実施した。同センター関係各位に有益なご指導とご協力をいただいた。また、本研究は、東北大学と気象庁気象研究所の共同研究の一環として行った。東北大学（指定国立大学）における災害科学世界トップレベル研究拠点による支援を受けた。

参考文献[1] Bauer, P., A. Thorpe, and G. Brunet, The quiet revolution of numerical weather prediction. Nature, 525,
47–55, 2015.
[2] Hersbach, H., B. Bell, P. Berrisford, S. Hirahara, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey,
R. Radu, D. Schepers, A. Simmons, C. Soci, S. Abdalla, X. Abellan, G. Balsamo, P. Bechtold, G. Biavati, J.
Bidlot, M. Bonavita, G. De Chiara, P. Dahlgren, D. Dee, M. Diamantakis, R. Dragani, J. Flemming, R.
Forbes, M. Fuentes, A. Geer, L. Haimberger, S. Healy, R. J. Hogan, E. Hólm, M. Janisková, S. Keeley, P.
Laloyaux, P. Lopez, C. Lupu, G. Radnoti, P. de Rosnay, I. Rozum, F. Vamborg, S. Villaume, and J.-N.
Thépaut, The ERA5 global reanalysis. Quart. J. Roy. Meteor. Soc., 146, 1999–2049, 2020.
[3] Fukui, S., T. Iwasaki, K. Saito, H. Seko, and M. Kunii, A feasibility study on the high-resolution
regional reanalysis over Japan assimilating only conventional observations as an alternative to the
dynamical downscaling. J. Meteor. Soc. Japan, 96, 565–585, 2018.
[4] Kobayashi, S., Y. Ota, Y. Harada, A. Ebita, M. Moriya, H. Onoda, K. Onogi, H. Kamahori, C.
Kobayashi, H. Endo, K. Miyaoka, and K. Takahashi, The JRA-55 Reanalysis: General specifications andbasic characteristics. J. Meteor. Soc. Japan, 93, 5–48, 2015.
[5] Saito, K., J. Ishida, K. Aranami, T. Hara, T. Segawa, M. Narita, and Y. Honda, Nonhydrostatic
atmospheric models and operational development at JMA. J. Meteor. Soc. Japan, 85B, 271–304, 2007.
[6] Hunt, B. R., E. J. Kostelich, and I. Szunyogh, Efficient data assimilation for spatiotemporal chaos: A
local ensemble transform Kalman filter. Physica D, 230, 112–126, 2007.
[7] Kobayashi, C., H. Endo, Y. Ota, S. Kobayashi, H. Onoda, Y. Harada, K. Onogi, and H. Kamahori,
Preliminary results of the JRA-55C, an atmospheric reanalysis assimilating conventional observations only.
SOLA, 10, 78-82, 2014.
[8] 曽我大輝, 日本域高解像度領域再解析の夏季における再現性に関する研究, 東北大学大学院
理学研究科地球物理学専攻修士論文, 2019.

― 4 ―

SENAC Vol. 55, No. （2 2022. 4）

図１2018 年 7 月の海面較正気圧場についての根平均二乗差。対象領域は図２(c)の領域。

図２

2018 年 7 月 5–7 日における 72 時間積算雨量。(a)解析雨量、(b)JRA-55、(c)領域再解析。
(a)及び(c)の灰色域は各データの領域外を示す。

SENAC Vol. 55, No. 2（2022. 4）― 5 ―［共同研究成果］

熱 損 失 を 考 慮 し た F l am el e t a pp ro a ch に 基 づ く
燃焼シミ ュレーシ ョンと N O 生成の予測
松川 嘉也：東北大学大学院工学研究科
小澤 龍磨：東北大学大学院工学研究科
赤尾津翔大：東北大学大学院工学研究科
松下 洋介：弘前大学大学院理工学研究科
青木 秀之：東北大学大学院工学研究科
Weeratunge Malalasekera：Loughborough University

本研究では，乱流拡散火炎を対象にふく射伝熱による熱損失を考慮した Flamelet approach に基づく燃焼シミュレーションを実施し，一酸化窒素 NO 生成も予測した．その結果，燃焼の Large Eddy
Simulation により主要な化学種と温度の時間平均値を概ね予測することができた．また，Enthalpy
defect と呼ばれる手法を導入することで，ふく射伝熱による熱損失に伴う温度低下を表現することができた．しかしながら，NO 生成の予測精度は Enthalpy defect よりもむしろ NO 生成の予測手
法自体の感度の方が大きかった．

1.

緒言

前報[1]ではメタンのパイロット火炎である Sandia Flame D [2-4]を対象に Flamelet/Progress-Variable approach (FPV) [5]に基づく燃焼の Large Eddy Simulation (LES)を実施した．その結果，他の研究[6]と同様に流れと混合分率の時間平均値だけでなくその分散もほぼ完全に再現することができた．また，主要な化学種に加えてラジカル種などのマイナー化学種の濃度の時間平均値についてもほぼ完全に予測することができた．そのため，FPV に基づく LES の有用性を確認することができたと言える．
現在，Flamelet approach として Flamelet model [7,8]，Flamelet-Generated Manifolds method (FGM)[9]と FPV [5]が広く用いられている．そもそも拡散燃焼を対象としている Flamelet model では
Flamelet library のパラメータは混合分率とスカラー消散率，FPV では混合分率と進行変数であり，いずれのモデルにおいてもエンタルピーは Flamelet library のパラメータではない．そのため，Flamelet model [7,8]や FPV [5]において熱損失を考慮するためには，Flamelet library のパラメータとしてエンタルピーを追加する必要がある．一方，予混合燃焼を対象としている FGM [9]ではFlamelet library のパラメータは進行変数とエンタルピーであり，エンタルピーはすでに Flamelet
library のパラメータである．そのため，FGM [9]は熱損失を表現可能である．しかしながら，FGM[9]を拡散燃焼に適用する場合，Flamelet library のパラメータに混合分率を追加する必要がある．
以上より，いずれの Flamelet approach もそのまま用いても熱損失を考慮した拡散燃焼シミュレー
ションを実施することはできない．なお，FPV [5]の Flamelet library のパラメータにエンタルピーを追加したモデルと FGM [9]の Flamelet library のパラメータに混合分率を追加したモデルはFPV [5]の Flamelet library を対向流拡散火炎から，FGM [9]の Flamelet library を予混合火炎からそれぞれ構築すること以外，基礎式や Flamelet library のパラメータなどは同一となる．
本研究では，水素の噴流の乱流拡散火炎である H5 Flame [10]を対象に Flamelet approach に基づ
く燃焼の Large Eddy Simulation あるいは Reynolds-Averaged Navier-Stokes を実施する．Flamelet
approach には Flamelet-Generated Manifolds method [9]あるいは前報と同様に Flamelet/Progress-Variable approach [5]を用いる．本研究では，さらに Enthalpy defect [11,12]と呼ばれる熱損失を考慮可能な手法を新たに追加して，ふく射伝熱による熱損失を考慮した水素の噴流の乱流拡散火炎

― 6 ―

SENAC Vol. 55, No. （2 2022. 4）を対象とした燃焼シミュレーションを実施する．また，種々の手法を用いて乱流拡散火炎における環境汚染物質の 1 つである一酸化窒素 NO 生成の予測も試みる．

2.

解析手法

2.1 解析対象と解析条件
解析対象は水素の噴流の乱流拡散火炎である H5 Flame である．内径𝐷𝐷𝐷𝐷 = 8 mm のパイプから 295
K の窒素で希釈した水素を 21.7 m/s で，周囲から co-flow として 295 K の空気を 0.3 m/s でそれぞれ供給する．解析領域は図 1 に示すバーナ上部の開放空間12.5𝐷𝐷𝐷𝐷 × 12.5𝐷𝐷𝐷𝐷 × 70𝐷𝐷𝐷𝐷であり，燃料ジェット基準のレイノルズ数は 6,200 である．

図1

2 .2

バーナの概念図 [10]Fla mele t lib rar y の 構 築

2 .2 .1 Fla mel et /Pr ogr ess -V ari abl e a pp roa ch (FP V)
Flamelet/Progress-Variable approach (FPV) [5]は Flamelet model の安定した燃焼の状態に加えてスカラー消散率では一意に決定できない不安定な燃焼の状態まで表現可能であることが特徴である．FPV では，Flamelet model と同様に対向流拡散火炎を対象に Flamelet library を作成するものの，燃焼を特徴づけるパラメータがスカラー消散率ではなく Progress Variable (PV)と呼ばれる進行変数とすることが Flamelet model と異なる．なお，本研究では PV を H2 O の質量分率と定義する．具体的には，Flamelet model と同様に種々のスカラー消散率に対して対
向流拡散火炎を対象に燃焼シミュレーションを実施し，混合分率とスカラー消散率に対してPV の正味の生成速度を含む変数を保存した後にスカラー消散率を PV に変換することで，混合分率と PV に対して変数を Flamelet library に保存する．
2 .2 .2 Fla mel et -Ge ner ate d Man ifo lds ( FGM )Flamelet-Generated Manifolds Method (FGM) [9]は，一次元の対向流拡散火炎ではなく予混合
火炎を対象に Flamelet library を作成する点で Flamelet model や FPV と異なる．ここでは，拡散燃焼を対象とした燃焼シミュレーションに FGM を適用することを考える．具体的には，目的の燃焼シミュレーションに先立ち，種々の等量比に対して一次元の予混合火炎の燃焼シ
ミュレーションを実施し，距離を Controlling Variable (CV)と呼ばれる進行変数に，等量比を混合分率にそれぞれ変換し，混合分率と CV に対して CV の正味の生成速度を含む変数をFlamelet library を保存する．なお，CV は PV と名称が異なるものの，本質的には同じパラメ

熱損失を考慮した Flamelet approach に基づく
燃焼シミュレーションと NO 生成の予測― 7 ―ータである [11]．
2 .2 .3 Ent hal py de fec t
Flamelet approach においてふく射伝熱などによる熱損失を考慮する場合，一般に Enthalpy defectと呼ばれる手法が用いられる．Enthalpy defect では，これまで説明した Flamelet library を構築する
条件に加えて，あらかじめ熱損失を与えた一次元の燃焼シミュレーションを実施し，Flamelet
library を構築する．著者らの知る限り，熱損失を与える手法には，系全体に対して一様に熱損失を与える方法(scheme 1) [12]，正味の熱発生速度(Heat release rate)に 1 より小さい値を乗じて熱損失を与える方法(scheme 2) [13]，人工的なふく射伝熱を考えることで熱損失を与える方法(scheme3) [13]がある．系全体に対して一様に熱損失を与える方法(scheme 1)は直感的で Flamelet library も構築しやすいものの，考慮することができる熱損失量に限界があり，大きな熱損失を与えることができない．正味の熱発生速度に 1 より小さい値を乗じて熱損失を与える方法(scheme 2)，人工的
なふく射伝熱を考えることで熱損失を与える方法(scheme 3)では大きな熱損失を与えることができるものの，熱損失が不均一となるため，Flamelet library の構築が煩雑になってしまう．本研究では，Flamelet library の構築が容易な scheme 1 を用いて熱損失を考慮した Flamelet library を構築する．
2.3 三次元の燃焼シミュレーションの解析方法
水素の噴流の拡散火炎(H5 Flame) [10]を対象とした乱流流れと Flamelet-Generated Manifolds
method (FGM)と Flamelet/Progress-Variable approach (FPV)の Large Eddy Simulation の基礎式は 3 次
元のファーブル平均を施した連続の式，運動量保存式，混合分率𝑍𝑍𝑍𝑍�と進行変数𝐶𝐶𝐶𝐶̃ の保存式である．
𝜕𝜕𝜕𝜕𝜌𝜌𝜌𝜌̅
𝜕𝜕𝜕𝜕+�𝜌𝜌𝜌𝜌𝜌𝜌𝜌𝜌𝜌�𝚥𝚥𝚥𝚥 � = 0𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗

𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕𝑝𝑝𝑝𝑝
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢�𝑖𝑖𝑖𝑖 𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢�𝑗𝑗𝑗𝑗 2 𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢�𝑘𝑘𝑘𝑘(𝜌𝜌𝜌𝜌𝑢𝑢𝑢𝑢𝑢𝑖𝑖𝑖𝑖 ) +�𝜌𝜌𝜌𝜌𝑢𝑢𝑢𝑢𝑢𝑖𝑖𝑖𝑖 𝑢𝑢𝑢𝑢𝑢𝑗𝑗𝑗𝑗 � = −+�(𝜇𝜇𝜇𝜇 + 𝜇𝜇𝜇𝜇𝑡𝑡𝑡𝑡 ) �+− 𝛿𝛿𝛿𝛿
��𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖 3 𝑖𝑖𝑖𝑖𝑗𝑗𝑗𝑗 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑘𝑘𝑘𝑘
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜇𝜇𝜇𝜇
𝜇𝜇𝜇𝜇𝑡𝑡𝑡𝑡 𝜕𝜕𝜕𝜕𝑍𝑍𝑍𝑍��𝜌𝜌𝜌𝜌𝑍𝑍𝑍𝑍�� +�𝜌𝜌𝜌𝜌𝑢𝑢𝑢𝑢𝑢𝑗𝑗𝑗𝑗 𝑍𝑍𝑍𝑍�� =��+��𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑍𝑍𝑍𝑍 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑍𝑍𝑍𝑍,𝑡𝑡𝑡𝑡 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗

∂
∂
∂
𝜇𝜇𝜇𝜇
𝜇𝜇𝜇𝜇𝑡𝑡𝑡𝑡 ∂𝐶𝐶𝐶𝐶̃��+�� + 𝜔𝜔𝜔𝜔𝜔𝐶𝐶𝐶𝐶
�𝜌𝜌𝜌𝜌𝐶𝐶𝐶𝐶̃ � +�𝜌𝜌𝜌𝜌𝑢𝑢𝑢𝑢𝑢𝑗𝑗𝑗𝑗 𝐶𝐶𝐶𝐶̃ � =∂𝜕𝜕𝜕𝜕
∂𝑥𝑥𝑥𝑥𝑗𝑗𝑗𝑗∂𝑥𝑥𝑥𝑥𝑗𝑗𝑗𝑗 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝐶𝐶𝐶𝐶 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝐶𝐶𝐶𝐶,𝑡𝑡𝑡𝑡 ∂𝑥𝑥𝑥𝑥𝑗𝑗𝑗𝑗

(1)

(2)
(3)
(4)

ここで，進行変数𝐶𝐶𝐶𝐶̃ は H2O の質量分率としている．
本研究では，ふく射伝熱による熱損失を表現可能ないわゆる Enthalpy defect を考慮するため，
ファーブル平均を施した混合分率𝑍𝑍𝑍𝑍�と進行変数𝐶𝐶𝐶𝐶̃ に加えてエンタルピーℎ�の保存式を解く．
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜇𝜇𝜇𝜇
𝜇𝜇𝜇𝜇𝑡𝑡𝑡𝑡 𝜕𝜕𝜕𝜕𝜕�
(5)�𝜌𝜌𝜌𝜌ℎ�� +�𝜌𝜌𝜌𝜌𝑢𝑢𝑢𝑢𝑢𝑗𝑗𝑗𝑗 ℎ�� =�� +�� − 𝛻𝛻𝛻𝛻 𝛻 q𝑟𝑟𝑟𝑟
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗 𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃 𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑡𝑡𝑡𝑡 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗
ここで，q𝑟𝑟𝑟𝑟 はふく射伝熱による熱損失を表すふく射熱流束である．
(6)−𝛻𝛻𝛻𝛻 𝛻 q 𝑟𝑟𝑟𝑟 = 𝑎𝑎𝑎𝑎(𝐺𝐺𝐺𝐺 𝐺 4𝜎𝜎𝜎𝜎𝜎𝜎𝜎𝜎 4 )𝐺𝐺𝐺𝐺は入射熱流束であり，P1 approximation method に基づき求める．
1 𝜕𝜕𝜕𝜕𝐺𝐺𝐺𝐺
q𝑟𝑟𝑟𝑟 = −(7)
3𝑎𝑎𝑎𝑎 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗
a はガスの吸収係数であり，Weighted-sum-of-gray-gasses model (WSGGM)を用いて推算する．
さらに，乱流拡散燃焼における NO 生成の予測を試みる．NO の濃度は他の化学種と同様 Flamelet
library を参照することで求めることができる．しかしながら，NO の生成の時間スケールは燃焼に関与する化学種の中では極めて大きく，流体の運動に対して時間スケールが十分に小さく，定常状態を仮定する Flamelet approach の仮定は成り立たないことが懸念される．そこで本研究では，
Flamelet library を参照することで求めるのではなく，NO の化学種の保存式を解くことで NO の質
量分率を求める．

― 8 ―

SENAC Vol. 55, No. （2 2022. 4）

𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑡𝑡𝑡𝑡

𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
𝜇𝜇𝜇𝜇
�𝜌𝜌𝜌𝜌𝑌𝑌𝑌𝑌�𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 � +�𝜌𝜌𝜌𝜌𝑢𝑢𝑢𝑢𝑢𝑗𝑗𝑗𝑗 𝑌𝑌𝑌𝑌�𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 � =��+𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗

𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗

𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁

𝜇𝜇𝜇𝜇𝑡𝑡𝑡𝑡

𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁,𝑡𝑡𝑡𝑡

𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕�� 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 � + 𝜔𝜔𝜔𝜔𝜔𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁

(8)

𝑗𝑗𝑗𝑗

ここで，𝜔𝜔𝜔𝜔𝜔𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 は NO の正味の生成速度であり，Flamelet library に保存した NO の正味の生成速度𝜔𝜔𝜔𝜔𝜔𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 を用+いる方法(case 1)およびFlamelet library に保存した正味の NO の生成速度𝜔𝜔𝜔𝜔𝜔𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 を生成速度𝜔𝜔𝜔𝜔𝜔𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 および消費速−𝜔𝜔𝜔𝜔̇𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁
𝑙𝑙𝑙𝑙𝑖𝑖𝑖𝑖𝑙𝑙𝑙𝑙
度を NO の質量分率𝑌𝑌𝑌𝑌�𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁で除した値 � 𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿
𝐿𝐿𝐿𝐿 に分離して評価する方法(case 2)を検討する．
𝑌𝑌𝑌𝑌𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁

2 .4

離散化スキームと数値解法−
𝜔𝜔𝜔𝜔𝜔𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁+𝜔𝜔𝜔𝜔𝜔𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 = 𝜔𝜔𝜔𝜔𝜔𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 + 𝑌𝑌𝑌𝑌�𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 � 𝐿𝐿𝐿𝐿𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 �
𝑌𝑌𝑌𝑌�𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁(9)非構造格子の有限体積法に基づき，二次中心差分法を 95%と一次風上差分法を 5%混合したいわゆるブレンディング・スキームと二次中心差分法を用いて運動量保存式の対流項と拡散項をそれぞれ離散化し，3 次の Adams-Bashforth 法を用いて時間進展させる．また，Hybrid法を用いて混合分率𝑍𝑍𝑍𝑍�，進行変数𝐶𝐶𝐶𝐶̃ ，エンタルピーℎ�と NO の質量分率𝑌𝑌𝑌𝑌𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁 の保存式の対流項と拡散項をそれぞれ離散化し，陰解法を用いて時間進展させる．圧力の解法には Simplified
Marker And Cell (SMAC) [14]を用いて非定常計算を実施し，時間平均値を求める．SMAC の圧力補正値の連立方程式の解法には AMGS [15]を，混合分率と進行変数の連立方程式の解法には多項式前処理付き安定化双共役勾配法[16]をそれぞれ用いた．求めた混合分率𝑍𝑍𝑍𝑍�，進行変数′′2 を用いて Flamelet library 𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓′′2 , 𝐶𝐶𝐶𝐶̃ , ℎ� �を参照�, 𝑍𝑍𝑍𝑍�𝐶𝐶𝐶𝐶̃ とエンタルピーℎ�に加えて混合分率の分散𝑍𝑍𝑍𝑍�

することで，目的の燃焼シミュレーションに必要な混合ガスの物性値である密度𝜌𝜌𝜌𝜌や粘度𝜇𝜇𝜇𝜇，
反応進行変数𝐶𝐶𝐶𝐶̃ の正味の反応速度𝜔𝜔𝜔𝜔𝜔𝐶𝐶𝐶𝐶 に加えて化学種の濃度や温度など興味あるあらゆる変数を求める．METIS-5.1.0 [17]を用いて解析領域を分割し，Message Passing Interface (MPI)を用
いて領域分割に基づく並列計算を実施した．また，通信の多い領域間から可能な限り同時に通信する独自のアルゴリズム[18]を採用し，通信に要する時間を最小化し，オーバヘッドを最小化することで並列化効率を向上した．なお，東北大学サイバーサイエンスセンター 所有の並列コンピュータ LX 406Re-2 を用い，4 ノード 48 コアを用いた並列計算を実施した．

3.

結果と考察

3 .1
RANS と LES の数値解および測定値との比較
まず，乱流流れの解析手法が乱流拡散燃焼の数値解に及ぼす影響を検討する．図 2 に Enthalpy
defect を考慮しない場合の，Reynolds-Averaged Navier-Stokes (RANS)あるいは Large Eddy Simulation(LES)に基づき，燃焼モデルには Flamelet/Progress-Variable approach (FPV)を用いて求めた中心軸上における混合分率，その Root-mean-square (RSM)，温度および H2 と H2 O の質量分率の時間平均値の分布を測定値とともに示す．RANS の数値解は測定値と比較して混合分率𝑍𝑍𝑍𝑍�の減衰を過大に見積もった結果，H2 の低下と H2 O の増加を過大に見積もり，温度の上昇を過大に見積もっている．
これは RANS では噴流の拡散を過大に評価するためであると考えられる．一方，LES の数値解は混合分率の測定値をほぼ完全に再現した結果，H2 の低下と H2 O の増加もほぼ完全に表現している．LES や RANS の数値解は下流において温度の測定値を若干過大に見積もっている．これは，燃焼により生成した吸収性ガスである H2 O が広範囲に分布する領域であり，ふく射伝熱による熱熱損失を考慮した Flamelet approach に基づく
燃焼シミュレーションと NO 生成の予測― 9 ―損失の影響であると考える．

(a) 混合分率とその変動(b) 温度(c) H2 の質量分率
H2O の質量分率
図 2 中心軸上における数値解と測定値
3 .2

Ent halp y def ect

次に，Enthalpy defect が乱流拡散燃焼の数値解に及ぼす影響を検討する．図 3 に Enthalpy
defect を考慮した場合としない場合の，Reynolds-Averaged Navier-Stokes (RANS)に基づき，燃焼モデルに Flamelet-Generated Manifolds method (FGM)あるいは Flamelet/Progress-Variableapproach (FPV)を用いて求めた中心軸上における温度の分布を測定値とともに示す．なお，
Enthalpy defect を考慮せず，燃焼モデルに FPV を用いて求めた数値解は前節と同一の結果である．用いた燃焼モデルによらず，Enthalpy defect を考慮することで，特に下流に向かって温度が低下している．これは，燃焼により生成した吸収性ガスである H2 O が生成し，半径方向に分布することでふく射伝熱による熱損失の影響が大きくなったためであると考える．
FGM の数値解は FPV と比較して温度の低下が若干大きく，熱損失の影響が若干大きい．これは，FGM の数値解が FPV と比較して燃焼反応の進行を若干早く見積もっているためである．以上より，FGM と FPV のパラメータに混合分率，その分散と進行変数にエンタルピー
を加え，Flamelet library を 4 次元のデータベースとすることで，熱損失を考慮することが可能であった．今後，LES に Enthalpy defect の影響を考慮することで，解析精度が向上することを示したいと考える．

― 10 ―

SENAC Vol. 55, No. （2 2022. 4）

図3

図4

中心軸上の温度分布の数値解と測定値中心軸上における NO の質量分率の数値解と測定値

3.3 NO 生成の予測
最後に，Enthalpy defect と NO 生成の予測方法が NO 生成の予測精度に及ぼす影響を検討する．
結果の一例として，図 4 に Reynolds-Averaged Navier-Stokes (RANS)に基づき，燃焼モデルに
Flamelet-Generated Manifolds method (FGM)を用いて求めた中心軸上における NO の質量分率の分布を測定値とともに示す．Enthalpy defect の考慮および NO 生成の予測方法によらず，NO の質量熱損失を考慮した Flamelet approach に基づく
燃焼シミュレーションと NO 生成の予測― 11 ―

分率は測定値と同様に上流から下流に向かって増加し，その後減少した．これは温度の上昇に伴い Thermal NO が生成し，その後半径方向に拡散したためであると考える．Enthalpy defect の考慮によらず，NO の質量分率は Flamelet library を参照することで求める方法では測定値を過小評価
し，生成項に正味に生成速度を有する NO の保存式を解く方法(case 1)では測定値を過大評価し，生成項に生成と消費に分離した NO の生成速度を有する NO の保存式を解く方法(case 2)では比較的測定値を再現した．また，いずれの方法においても Enthalpy defect を考慮することで NO の濃度が低下した．これは，H5 Flame では NO は Thermal NO により生成するが，ふく射伝熱による熱損失を考慮することで火炎の温度が低下し，Thermal NO の生成が抑制されたためである．ただし，本結果からいずれの方法が NO 生成の予測に最も適しているか結論づけるのは難しく，今後さらに検討を進める必要があると考える．

4.

結言

本研究では，乱流拡散火炎を対象に Enthalpy defect によりふく射伝熱による熱損失を考慮した
Flamelet/Progress-Variable approach を用いた燃焼の Large Eddy Simulation (LES)あるいは Reynolds-Averaged Navier-Stokes (RANS)を実施し，NO 生成の予測も試みた．その結果，燃焼の LES により主要な化学種と温度の時間平均値を概ね予測することができた．また，Enthalpy defect を導入することで，ふく射伝熱による熱損失に伴う温度低下を表現することができた．NO 生成の予測
精度は Enthalpy defect よりもむしろ NO 生成の予測手法自体の感度の方が大きかった．

謝辞
本研究は，東北大学サイバーサイエンスセンターの並列コンピュータ LX 406Re-2 を利用することで実現することができた．また，本研究の一部は特別研究員奨励費(18J11135)，JSPS 科研費(JP18K03964)および東北大学若手リーダー研究者海外派遣プログラムの助成を受けたものであり，
ここに謝意を表す．
参考文献[1] 赤尾津翔大ら, Flamelet approach に基づくガス燃料を対象とした燃焼シミュレーション,
SENAC, 52(3), 1–7 (2019)[2] Barlow, R. S. and Frank, J. H., Effects of turbulence on species mass fractions in methane/air jet
flames, Proc. Combust. Inst. 27, 1087–1095 (1998)[3] Barlow, R. S. et al., Piloted methane/air jet flames: Scalar structure and transport effects, Combust.
Flame, 143, 433–449 (2005)[4] Schneider, Ch. et al., Flow field measurements of stable and locally extinguishinghydrocarbon-fuelled jet flames, Combust. Flame, 135, 185–190 (2003)[5] Pierce, C. and Moin, P., Progress-variable approach for large-eddy simulation of non-premixed
turbulent combustion, J. Fluid Mech., 504 73–97 (2004)[6] Raman, V. and Pitsch, H, A consistent LES/filtered-density function formulation for the simulation of
turbulent flames with detailed chemistry, Proc. Combust. Inst., 31(2), 1711–1719 (2007)[7] Peters, N., Laminar diffusion flamelet models in non-premixed turbulent combustion, Prog. Energy
Combust. Sci., 10(3), 319–339 (1984)[8] Peters, N., Laminar flamelet concepts in turbulent combustion, Symp. (Int.) on Combust., 21, 1231–
1250 (1998)[9] van Oijen, J. A. and de Goey, L. P. H., Modelling of Premixed Laminar Flames using Flamelet-Generated Manifolds, Combustion Science and Technology. 161, 113–137 (2000)[10] W. Meier et al., Characterization of Turbulent hytVAir Jet Diffusion Flames by Single-Pulse
Spontaneous Raman Scattering, Combust. Sci. Technol. 118, 293–312 (1996)

― 12 ―

SENAC Vol. 55, No. （2 2022. 4）[11] 松下洋介ら, Flamelet Model に基づく乱流燃焼シミュレーション, 金属, 85(11), 915–921 (2015)[12] Hossain, M. et al., Modelling of a Bluff-Body Nonpremixed Flame using a Coupled
Radiation/Flamelet Combustion Model, Flow, Turbul. Combust., 67(3), 217–234 (2001)[13] Wollny, P. et al., Modelling heat loss effects in high temperature oxy-fuel flames with an efficientand robust non-premixed flamelet approach, Fuel, 216 44–52 (2018)[14] Harlow, F. H. and Welch, J. E., Numerical calculation of time-dependent viscous incompressible flowof fluid with a free surface, Phys. Fluids, 8 2182–2189 (1965)[15] 藤井昭宏ら, 領域分割による並列 AMG アルゴリズム, 情報処理学会論文誌コンピューティ
ングシステム（ACS）, 44, SIG06(ACS1), 9–17 (2003)[16] Van der Vorst, H. A., Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the
Solution of Nonsymmetric Linear Systems, SIAM J. Sci. Stat. Comput., 13(2), 631–644 (1992)[17] Karypis, G. and Kumar, V., A Fast and Highly Quality Multilevel Scheme for Partitioning Irregular
Graphs, SIAM J. Sci. Comput., 20(1), 359–392, (1999)[18] Matsushita, Y. et al., Efficient communication strategy in parallel computation based on domain
partitioning, J. Chem. Eng. Jpn., 51(1), 79–82 (2018)

SENAC Vol. 55, No. 2（2022. 4）― 13 ―［共同研究成果］

Computational analysis of the body-force models forplasma actuator-assisted flow control simulations
Di Chen*: Dept. Information and Comp. Tech., Tokyo University of Science
Kengo Asada: Dept. Information and Comp. Tech., Tokyo University of Science
Satoshi Sekimoto**: Dept. Information and Comp. Tech., Tokyo University of Science
Kozo Fujii: Dept. Information and Comp. Tech., Tokyo University of Science
Hiroyuki Nishida: Tokyo University of Agriculture and Technology*

Currently, Institute of Industrial Science, University of Tokyo
Currently, Tokyo University of Agriculture and Technology**

1．
． Introduction
The well-known “ion wind” induced by a dielectric barrier discharge plasma actuator (DBDPA) has been extensively used as an active flow control device to avoid flow separation thanks to thehigh availability and feasibility. DBD-PA only consists of two electrodes with a dielectric layer
between them, see Fig. 1(a). Due to the thin and light structure, DBD-PA can be attached on any flator curved surfaces, corner, or edge, where flow separation control1, 2is considered without changing
the original shape of the body surface. A review conducted by Wang et al.3 shows a number of flow
control applications of DBD-PA till the year 2013. The applications include not only on separation
control shown in Fig. 1(b,c)4,5 but also on noise reduction,6 skin-friction reduction,7and else.
DBD-PA for flow separation control is typically installed on the suction surface near the leadingedge1, 2 It can generate a wall-jet flows with the maximum velocity up to around 10m/s, when thepeak-to-peak voltage (𝑉𝑉𝑝𝑝𝑝𝑝 ) and the base frequency of operating alternating current (AC) are 5-20kV
and 1-10kHz, respectively. The ionization effect is largely determined by the applied voltage as wellas electric permittivity and thickness of the dielectric layer. The mechanism of plasma-assisting flow
control is described as the interaction of the ionized gas and the neutral air, which results in anelectrohydrodynamic body-force field coupling with the momentum transfer in the external flow
downstream of the exposed electrode.8-11
Massines et al.12 and Roth et al.,13 proposed the early model of body force field generated by asingle DBD-PA, which was one-dimensional (1-D) based on static formulation and neglected thepresence of the charged particles, therefore it barely fitted two-dimensional (2-D) or threedimensional (3-D) applications. Semi-empirical models of 2-D plasma flow modelling, using linear,14
exponential functions,10, 11 and Gaussian distribution2 of the spatial decay for the 2-D body force
component were reviewed by Corke et al.15 Numerous studies have made great efforts on the bodyforce modelling.2, 12-14

― 14 ―

SENAC Vol. 55, No. （2 2022. 4）

Figure 1. (a) Sketch of DBD-PA. (b) Wind tunnel test on airfoil4. (c) DBD-PA attached on theleading edge of a small model plane (wingspan of 3m) in a flight test5.
The phenomenological such as Suzen-Huang (S-H) model2 and Shyy model14 were largely
employed because of the simplicity in the computational flow simulations of Asada et al.16,17 andVisbal et al.,1 respectively. The well-known analytical S-H model is an electrostatic (semi-empirical)
model derived from Enloe et al.’s work,10 the body force vector can be expressed as,(1)
𝒇𝒇𝑖𝑖 = 𝑄𝑄𝑐𝑐 𝑬𝑬𝑖𝑖 = 𝑄𝑄𝑐𝑐 (−∇𝜙𝜙),

where 𝑬𝑬𝑖𝑖 is the electric field vector. The force is contributed by two different parts: the external
electric potential 𝜙𝜙 and the electric field created by the net charged density 𝑄𝑄𝑐𝑐 . The net body force
obtained by the closed-form solution of the S-H model was well validated in the experimental study,21
by calibrating the nondimensional plasma parameter 𝐷𝐷𝑐𝑐 representing the scaling of the electrical toinertial forces, which also describes the strength and the scale of body force from anotherperspective.16,17, 21,22 These studies indicated that the flow control authority of the plasma actuators
are well realized by the simulations using these models. Airfoil stall characteristics are well simulated,
and the computational and experimental pressure distributions over an airfoil agree pretty well.
Despite the low cost of using the analytical model, it is true that the S-H model neglects thecomplex plasma chemistry which usually causes the highly unsteady forcing on the plasma flow. On
the other hand, the charged-particle models associated with the fully-coupled approaches consider theordinary force diffusion, drift motion, Coulomb acceleration of electrons, and positive and negativeions, respectively. Drift-Diffusion (D-D) model was first developed from the physical modelling,23-25
which focused on the electric-field effects on the charged particles. More recently the D-D model wasused to simulate the discharge plasma evolution and the DBD-PA-induced body-force field.9,26
However, due to the time-consuming computation, few studies applied the D-D model in flow control
simulations. Gaitonde et al.27 conducted plasma-based stall control simulations with coupledapproaches that largely reduced the complexity of the broad-spectrum problem, nonetheless, theinduced flow field of high temporal resolution during a single discharge cycle still remains unclear.
In the present study, the body-force field was obtained by the S-H model, as well as by the D-Dcomputation of high temporal resolution. A high-fidelity Navier-Stokes solver is employed to simulate
the DBD-PA induced flow with the incorporating body force field. Based on the published results by
Chen et al.,28 this report summarizes the nature of the S-H model and the D-D model, clarifies thedifference between the two models from the viewpoint of the DBD-PA flow control authority, andfinally provides several suggestions in using the body-force models for plasma actuator-assisted flow
control simulations.

Computational analysis of the body-force models forplasma actuator-assisted flow control simulations― 15 ―

2．
． Methodology
2.1 Body force modelling
2.1.1 Suzen-Huang model
As we introduced in Eqn. (1), the body force vector is computed by multiplying the charge
density 𝑄𝑄𝑐𝑐 and electric field vector 𝑬𝑬𝑖𝑖 , which are solved in the Maxwell’s equations of the external
electric potential 𝜙𝜙 and the charged particle potential, respectively, as follows∇ ∙ (𝜀𝜀𝑟𝑟 ∇𝜙𝜙) = 0,(2)
𝑄𝑄𝑐𝑐
∇ ∙ (𝜀𝜀𝑟𝑟 ∇𝑄𝑄𝑐𝑐 ) = 2 ,(3)𝜆𝜆𝑑𝑑
where 𝜀𝜀𝑟𝑟 denotes the relative permittivity of the dielectric layer, and 𝜆𝜆𝑑𝑑 denotes the Debye length.
In Fig. 2(a), the boundary conditions for solving Eqn. (2) and (3) on the exposed electrode and thewall above the covered electrode can be written as𝜙𝜙(𝑡𝑡) = 𝜙𝜙 𝑚𝑚𝑚𝑚𝑚𝑚 𝑓𝑓(𝑡𝑡),(4)(5)
𝑄𝑄𝑐𝑐,𝑤𝑤 (𝑥𝑥, 𝑡𝑡) = 𝑄𝑄𝑐𝑐𝑚𝑚𝑚𝑚𝑚𝑚 𝐺𝐺(𝑥𝑥)𝑓𝑓(𝑡𝑡),respectively, where 𝜙𝜙 𝑚𝑚𝑚𝑚𝑚𝑚 and 𝑄𝑄𝑐𝑐𝑚𝑚𝑚𝑚𝑚𝑚 are the maximum values of the external electric potential andthe charge density, respectively. Consequently, the time variation of the non-dimensional S-H bodyforce is given as follows(6)
𝑺𝑺𝑖𝑖 (𝑥𝑥, 𝑦𝑦, 𝑧𝑧, 𝑡𝑡) = 𝐷𝐷𝑐𝑐 sin2 (2𝜋𝜋𝐹𝐹𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 𝑡𝑡) 𝑺𝑺𝑆𝑆𝑆𝑆 (𝑥𝑥, 𝑦𝑦, 𝑧𝑧),
where 𝐹𝐹𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 is the nondimensional base frequency of the AC power, 𝑺𝑺𝑆𝑆𝑆𝑆 is the non-dimensional
body-force vector of S-H model. 𝑺𝑺𝑆𝑆𝑆𝑆 follows a half Gaussian function 𝐺𝐺(𝑥𝑥) in two-dimensional
spatial distribution described in Fig. 2(b).2 The magnitude of the S-H body force is manipulated by
the nondimensional plasma parameter 𝐷𝐷𝑐𝑐 , which is given empirically as
𝑄𝑄𝑐𝑐,𝑟𝑟𝑟𝑟𝑟𝑟 𝐸𝐸𝑟𝑟𝑟𝑟𝑟𝑟 𝐿𝐿𝑟𝑟𝑟𝑟𝑟𝑟
𝐷𝐷𝑐𝑐 =.
(7)𝜌𝜌𝑟𝑟𝑟𝑟𝑟𝑟 𝑈𝑈𝑟𝑟𝑟𝑟𝑟𝑟 2
However, it needs further calibration by experiments and numerical simulations according to the
induced flow field.21
2.1.2 Drift-Diffusion modelNotwithstanding the capability of the D-D model in the 3-D body force simulation by Nishida
et al.,26 the model comparison is conducted in 2-D simulations for convenience and simplicity in thisstudy. The electron, the positive ion, and the negative ion are considered with basic plasma chemistry
including electron impact ionization, attachment, and recombination. The ionization and attachment
coefficients and electron mobility are calculated by BOLSIG+29 a free electron Boltzmann equation
solver, assuming the ambient gas is air (N2:O2 = 0.8:0.2).
Here the governing equations to be solved are explained in detail in the previous studies,9, 24-26,
28, 30the computational electrohydrodynamic (EHD) force regarded as the two-dimensional bodyforce is expressed as the rate of momentum transfer per unit volume due to collisions,
𝐷𝐷𝑝𝑝
𝐷𝐷𝑒𝑒
𝐷𝐷𝑛𝑛
𝒇𝒇 = 𝑒𝑒𝑐𝑐 (𝑛𝑛𝑝𝑝 − 𝑛𝑛𝑒𝑒 − 𝑛𝑛𝑛𝑛 )𝑬𝑬 − [ ∇𝑛𝑛𝑝𝑝 + ∇𝑛𝑛𝑒𝑒 +∇𝑛𝑛 ],
(8)
𝜇𝜇𝑝𝑝
𝜇𝜇𝑒𝑒
𝜇𝜇𝑛𝑛 𝑛𝑛where 𝑛𝑛 is the plasma density, 𝜇𝜇 is the charged particle mobility, and their subscript e, p, n denotethe electron, the positive ion, and the negative ion, D is the coefficients of diffusion. The first term on
the right side of Eqn. (8) plays a dominant role and corresponds to the expression used in S-H modelin Eqn. (1), while the second term represents the diffusion effect. The obtained 2-D body force in Eqn.
(8) is first averaged in phases, then input into the N-S solver for the 2-D simulation of the induced
flow, however the spanwise uniformity of body-force distribution is assumed in the 3-D simulation
for simplicity.

― 16 ―

SENAC Vol. 55, No. （2 2022. 4）

Figure 2. (a) Boundary condition of charged particle in S-H model, following a half
Gaussian distribution.2 (b) Spatial and temporal distribution of body force in S-Hmodel.
2.2 Computational fluid dynamics (CFD)
Two typical flows are considered. First, the flow field is assumed to be globally quiescent andlaminar over a flat plate, then the strength of the induced flow is investigated. Second, the body-force
models are applied and examined in practical flow control simulations, where the separated flow overa NACA0015 airfoil with a light-stall angle of attack of 12° in turbulence transition state (Reynolds
number is 63,000) is investigated. The Mach number is 0.2, the specific heat ratio (𝛾𝛾) is 1.4, and the
Prandtl number (Pr) is 0.72, those are keeping the same as the previous simulation16 and theexperiment setup.21 On the surface of the flat plate and the airfoil, no-slip and adiabatic conditions
are imposed.
The flow field is described by the 3-D compressible Navier-Stokes equations as below,
augmented by the term 𝑺𝑺𝑖𝑖 representing the local forcing on the ionized region by DBD-PA. The nondimensional forms of the continuity, momentum, energy equations, as well as ideal gas equation are
written as follows:𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕𝒖𝒖𝑘𝑘
+
= 0,
(9)
𝜕𝜕𝜕𝜕
𝜕𝜕𝑥𝑥𝑘𝑘𝜕𝜕𝜕𝜕𝒖𝒖𝑖𝑖 𝜕𝜕(𝜌𝜌𝒖𝒖𝑖𝑖 𝒖𝒖𝑘𝑘 + 𝑝𝑝𝜹𝜹𝑖𝑖𝑖𝑖 )1 𝜕𝜕𝝉𝝉𝑖𝑖𝑖𝑖
+=+ 𝑺𝑺𝑖𝑖 ,(10)
𝜕𝜕𝒙𝒙𝑘𝑘𝑅𝑅𝑅𝑅 𝜕𝜕𝒙𝒙𝑘𝑘
𝜕𝜕𝜕𝜕
1 𝜕𝜕𝒖𝒖𝑙𝑙 𝝉𝝉𝑘𝑘𝑘𝑘
1
𝜕𝜕𝒒𝒒𝑘𝑘𝜕𝜕𝜕𝜕 𝜕𝜕((𝑒𝑒 + 𝑝𝑝)𝒖𝒖𝑘𝑘 )+=(−) + 𝑺𝑺𝑘𝑘 𝒖𝒖𝑘𝑘
(11)2(𝛾𝛾 − 1)𝑃𝑃𝑃𝑃𝑀𝑀∞ 𝜕𝜕𝒙𝒙𝑘𝑘
𝑅𝑅𝑅𝑅 𝜕𝜕𝒙𝒙𝑘𝑘
𝜕𝜕𝒙𝒙𝑘𝑘𝜕𝜕𝜕𝜕
𝑝𝑝 = (𝛾𝛾 − 1)(𝑒𝑒 − 0.5𝜌𝜌𝒖𝒖𝑘𝑘 𝒖𝒖𝑘𝑘 )(12)
where 𝒙𝒙𝑖𝑖 is the position vector, 𝒖𝒖𝑖𝑖 is the velocity vector, 𝒒𝒒𝑘𝑘 is the heat flux vector, 𝜌𝜌 is thedensity, 𝑝𝑝 is the static pressure, 𝑒𝑒 is the total energy per unit volume, 𝝉𝝉𝑖𝑖𝑖𝑖 is the viscous stresstensor, 𝜹𝜹𝑖𝑖𝑖𝑖 is the Kronecker delta, 𝑺𝑺𝑖𝑖 is the body force vector, 𝛾𝛾 is the ratio of specific heats, 𝑡𝑡 is
the time. In S-H model, 𝑺𝑺𝑖𝑖 is directly given by Eqn. (6) with a proper 𝐷𝐷𝑐𝑐 value, while 𝒇𝒇𝑖𝑖 in D-D
model (Eqn. (8)) is normalized as
1
𝒇𝒇 .
(13)
𝑺𝑺𝑖𝑖 =2 𝑖𝑖
𝜌𝜌∞ 𝑈𝑈∞
2.3 Numerical schemes
To solve the governing equations in Eqn. (9) – (12), this research employs a compressible NavierStokes flow solver LANS3D,31,32 which is extensively used and well validated in recent 30 years. The
implicit large-eddy simulations (LES) are conducted for solving the more complex flow field above
the airfoil in turbulent transition state.16 All the spatial derivatives are obtained with a sixth-order

Computational analysis of the body-force models forplasma actuator-assisted flow control simulations― 17 ―

compact difference scheme.33 Lower-upper symmetric alternating direction implicit and symmetric
Gauss-Seidel (ADI-SGS) method is utilized for time integration. 10th order filtering33 is applied witha filtering coefficient of 0.42.
The current numerical schemes in D-D force computation largely follow Nishida et al.’s previous
studies.9, 26, 30 given by the formula in Eqn. (6). The input time interval of the S-H and the D-D force
match the time step of CFD simulation in real scale.
2.4 Computational grids
The information of the computational grid systems for the body force computation and the CFD
are listed in Table 1. In the 2-D simulation of the quiescent air, the 2-D computational domains andgrids are shown in Fig. 3(a), the D-D force (red zone) is computed and interpolated into the CFD grids(white background). In the 3-D simulations of the separated flow above the airfoil, the C-type grid isdesigned, the near and far field are shown in Fig. 3(b) and (c), respectively. The spanwise length ofthe computational region is set to 0.2 times the chord length. The D-D force is interpolated into theactuator zone (highlighted in green). A message passing interface (MPI) between the fluid zones is
employed for the parallel computing in CFD simulations.

Figure 3. Domain incorporation of the D-D body force (in red) into CFD (in white) for the flat
plate (a), and the NACA0015 airfoil (b), the green grids denote the overset actuator zone, theexposed electrode is highlighted in the zoom-in snip; (c) Computational domain of the airfoil andthe schematic of inflow, 𝐿𝐿𝑟𝑟𝑟𝑟𝑟𝑟 = 𝑐𝑐 = 0.1m.

Table 1. Grid density of the body-force computation and CFD
Computation domain𝜉𝜉𝜂𝜂𝜁𝜁D-D body force (red in Fig. 3a, b)
600
3
250
Flat plate (white in Fig. 3a)
873
3
416
Airfoil surface (white in Fig. 3b)
759
134
179

Total points
450,000
1,089,504
18,205,374

― 18 ―

SENAC Vol. 55, No. （2 2022. 4）

Airfoil actuator (green in Fig. 3b)
149
134
111
2,216,226* In curvilinear coordinates (𝜉𝜉, 𝜂𝜂, 𝜁𝜁), 𝜉𝜉 and 𝜁𝜁 denote wall-parallel and wall-normal directions,
respectively, 𝜂𝜂 denotes the spanwise direction in the airfoil case.

2.5 DBD-PA setup
Fig. 4(a) shows the configuration of DBD-PA in the D-D computational domain, which isconsistent with Sekimoto et al.’s experiment.34 The physical parameters as shown in Fig. 4(a) areprecalculated in the previous numerical studies using D-D model.26, 29 The plasma actuator consistsof a 160𝜇𝜇𝜇𝜇 thick Kapton dielectric layer and two 50𝜇𝜇𝜇𝜇 thick copper electrodes. The thickness ofthe exposed electrode can be neglected in CFD simulations. A sinusoidal form 0.5 𝑉𝑉𝑝𝑝𝑝𝑝 sin (2𝜋𝜋𝑓𝑓𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 𝑡𝑡)
of AC power is applied on the electrodes, the peak-to-peak voltage (𝑉𝑉𝑝𝑝𝑝𝑝 ) amplitude is from 7kV to
20kV and the base frequency (𝑓𝑓𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 ) is set to 10kHz. In the quiescent flow simulation, the continuousmode is applied, therefore the flow is permanently driven by the DBD-PA.
In the separated flow over the airfoil, the exposed electrode of DBD-PA is placed at 5% chord
length from the leading edge, the DBD-PA is periodically activated in the duty cycle as shown in Fig.
4(b). It is called burst mode actuation compared to the continuous mode actuation generally used. The
burst frequency 𝑓𝑓 + = 500Hz and the burst ratio BR = 0.1. The burst mode with a proper 𝐹𝐹 + is wellknown to have better performance in flow separation control than the continuous mode.15,16, 19, 22, 34
𝐹𝐹 + and 𝐹𝐹𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 are the normalized values of burst and base frequency, respectively.

Figure 4. (a) PA configuration and physical parameters in the D-D computation, including
electrodes and dielectric layer highlighted in yellow and grey, respectively (𝐿𝐿𝑟𝑟𝑟𝑟𝑟𝑟 = 0.1𝑚𝑚);(b) Schematic diagram of a bursting wave in burst actuation mode.
2.6 Computational resources
To obtain the time-varying D-D force fields which are necessary as an input data for the flow
field simulations, an in-house computation server is used. It takes around 90 hours using one node
per case, which computes 4.5 discharge cycles. For the CFD simulations, the supercomputer
subsystem named AOBA-A of the Cyberscience Center in Tohoku University is used. 8 nodes areused for computing the flow fields in the divided subdomains with MPI techniques, it takes around
26 hours per case (non-dimensional time 𝑡𝑡𝑈𝑈𝑟𝑟𝑟𝑟𝑟𝑟 /𝐿𝐿𝑟𝑟𝑟𝑟𝑟𝑟 = 80) of the flat-plate flow simulation, andaround 70 hours per case (non-dimensional time 𝑡𝑡𝑈𝑈∞ /𝑐𝑐 = 9.6) of the airfoil flow simulation.
． Results
3．

3.1 Body force fields of DBD-PA
The body force terms on the right side of Eqn. (10) and (11) are computed by the D-D modeland the S-H model. Fig. 5 show the time history of the wall-parallel EHD net force within 3 cycles,Computational analysis of the body-force models forplasma actuator-assisted flow control simulations― 19 ―

the D-D force computations are conducted at the 𝑉𝑉𝑝𝑝𝑝𝑝 from 7kV to 20kV. In the case of S-H model,
the nondimensional plasma parameter 𝐷𝐷𝑐𝑐 is set to 0.0117 referred to the previous parametricstudies,21 in which the near-field induced flow agrees well with the experiment of 𝑉𝑉𝑝𝑝𝑝𝑝 = 7𝑘𝑘𝑘𝑘. Dueto the lack of strong evidence relating the controlling parameter 𝐷𝐷𝑐𝑐 to the applied voltage, we only
plot the results of the S-H model of 𝐷𝐷𝑐𝑐 = 0.0117 (SH-Dc0.0117) at the corresponding voltage of
7kV for all the following comparisons.
As shown in Fig. 5 the time variation of body force in S-H model which was developed by Asadaet al.16, 17 assumes to follow the function of sin(2𝜋𝜋𝐹𝐹𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 𝑡𝑡) in Eqn. (14). The net force is positive all
the time and symmetrical between negative- and positive-going phases, that probably overestimatesthe input power when the body force peak is approximately same with that in the D-D case of 7kV.
The time-averaged value of the net body force in Fig. 5(c) are 10.21mN for D-D model at 7kV,
18.45mN for 10kV, 21.07mN for S-H model with 𝐷𝐷𝑐𝑐 = 0.0117. The power law of 𝑉𝑉𝑝𝑝𝑝𝑝1.49 agrees
with Nishida et al.’s 2-D simulation results26 in a range of 𝑉𝑉𝑝𝑝𝑝𝑝 from 10kV to 24kV when the relative
permittivity of the dielectric layer (glass epoxy) is set to 5.0.

Figure 5. Time history (three periods) of the dimensional EHD net force computed in the entire
body-force domain using the D-D model at 7kV, 8kV, 9kV, 10kV (a), and 12kV, 15kV, 20kV (b),as well as the S-H model with 𝐷𝐷𝑐𝑐 = 0.0117; (c) Time-averaged D-D body force .vs. peak-topeak voltage and S-H body force.

3.2 Induced flow in quiescent field

When the body force fields are obtained by D-D or S-H model, it is quite straightforward to
relate the force production to the induced flow structure in quiescent air over a flat plate. The timevarying body force terms are incorporated into the right-hand side of the N-S equations, the induced
flow fields are then computed by the CFD solver. The time-averaged wall-parallel velocity contours
using the D-D model from 7kV to 9kV are compared with the corresponding experimental results,34

― 20 ―

SENAC Vol. 55, No. （2 2022. 4）

as well as the S-H model with 𝐷𝐷𝑐𝑐 = 0.0117 in Fig. 6(a), which is determined based on the previous
study, where Aono et al. conducted showed the guideline of Dc values by the comparison ofexperimental and computational result with S-H model.21 Fig. 6(b) shows the maximum velocity
measured at 20mm downstream (position A), the far-field numerical and experimental results arecompared. The current results are obtained in the flow field averaged from 0.4s to 0.8s.

Figure 6. (a) Wall-parallel velocity fields in the present CFD simulations (left) and our previous
experiments34 (right); (b) Far-field velocity (𝑢𝑢𝐴𝐴 ) at 20mm downstream.

The induced flow structures and the downstream velocity (𝑢𝑢𝐴𝐴 ) computed by the D-D model show
very good agreement with the experiment. Compared to 8kV, the velocity decline at 9kV is probably
caused by the less body force production in Fig. 5(c) due to the unstable EHD force pushing fluid
towards downstream in positive-going phase. However, the induced flow of S-H model is muchstronger than those of D-D models from 7kV to 9kV, the local velocity at 20mm downstream is evenhigher than the D-D case of 12kV in Fig. 6(b).
Fig. 7 further shows the near-field velocity profiles in a range of 𝑥𝑥/𝐿𝐿𝑟𝑟𝑟𝑟𝑟𝑟 = 0.02 to 0.08 (2mmto 8mm) in the downstream of the exposed electrode, as shown in the top of Fig. 7(b). The experimentof 9kV shows the remarkable unsteadiness, which is the potential cause of the reversed induced flowat 𝑥𝑥/𝐿𝐿𝑟𝑟𝑟𝑟𝑟𝑟 = 0.02, and the exceeding velocity compared to the D-D result at the other downstream
positions. On the other hand, the near-field velocity profiles of the S-H model show the much smaller
gradient of the velocity field, compared to the D-D results.

Computational analysis of the body-force models forplasma actuator-assisted flow control simulations― 21 ―

Figure 7. (a) Wall-parallel velocity profiles in the near field; (b) Thrust at 6mm downstream.
3.3 Flow separation control over NACA0015 airfoil
The PA-assisted flow control have gained great interest in a variety of separated flow fields overthe Reynolds number from 103 to 106 ,36and the angle of attack of pre-to-post stall from 10 degrees
to 16 degrees. In this report, a comparative study among the D-D model, S-H model and experimentis conducted in the separated flow over a NACA0015 airfoil at a light-stall angle of attack of 12°
and a transitional Reynolds number of 63,000. Burst actuation (𝐹𝐹 + = 5, BR = 0.1) is set to maximize
the separation control effect.35, 37 In the following results, “EXP-7kV” denotes the experimental resultwith the same setups of flow filed and DBD-PA in Asada et al.’s work,37 “NOACT” is the case without
control.
Fig. 8 shows that the flows are almost completely attached in all the three controlled cases, italso shows little difference between the S-H and the D-D cases either in velocity field or in TKE(turbulent kinetic energy) field. The TKE has overshoot in the transition position of the shear layer,
and the overshoot becomes larger when the transitional shear layer is thick.36 The chordwise locationsof 𝑇𝑇𝑇𝑇𝑇𝑇𝑚𝑚𝑚𝑚𝑚𝑚 peaks in Fig. 8(c) show very good agreement with the reattachment positions in Fig.
9(b). However, the S-H case shows the relatively larger TKE overshoot compared to the D-D cases,
due to the delay of turbulent transition. On the contrary, the early and smooth transition process in theD-D cases results in the smaller TKE overshoots.

― 22 ―

SENAC Vol. 55, No. （2 2022. 4）

Figure 8. Distribution of chordwise velocity (a), turbulence kinetic energy (b); Spatial
development of maximum turbulent kinetic energy (c) obtained from (b).
Distributions of the aerodynamic coefficients in the chordwise direction are further shown in Fig.
9. In the cases with PA control, the time- and span-averaged pressure coefficient (𝐶𝐶𝑝𝑝 ) shows the
suction peaks as well as the plateau distributions which corresponds to the so-called laminar
separation bubble. The 𝐶𝐶𝑝𝑝 curves indicate the lift increase in the cases with PA control, while theskin-friction coefficient (𝐶𝐶𝑓𝑓 ) shows the reduction of pressure drag, as well as the flow-separated
regions shown in the bottom of Fig. 9(b). Better agreement of 𝐶𝐶𝑝𝑝 distribution can be observed
between the D-D cases and the experiment at 7kV, compared with the S-H case. As to the 𝐶𝐶𝑓𝑓 curves
in Fig. 9(b), it is found that the reattachment points in the D-D cases move closer towards the leadingedge compared to the S-H case. The discrepancy of the reattachment position is partly caused by thedifferent location of the maximum induced velocity in the quiescent flow simulation, which is
discussed in detail by Chen et al.28

Computational analysis of the body-force models forplasma actuator-assisted flow control simulations― 23 ―

Figure 9. Pressure coefficient (a), skin-friction coefficient (b) along the airfoil chord.
4．
． Summary and suggestions
In computational studies of PA-assisted flow control, S-H model can be quickly implemented
because of the low cost, it has been widely used and also well validated by experiments. On the other
hand, D-D model as a high-fidelity body force model has never been employed before in practicalflow control simulations. The study shown in Ref. 28 and summarized in this report by the present
authors first shows the usage of the two models, as well as the result comparison in the quiescent flow
and airfoil flow simulation. Using the different body force model, the effect of separation control isalmost the same in the separated airfoil flow, despite the large discrepancy of the induced flow in thequiescent air.
When we consider the cost reduction of the body force computation, for the current flow
condition in light separated and turbulent transition state, S-H model is acceptable within thereasonable range of computational accuracy so far as discussions are qualitative. However, it isprobably necessary to use the D-D model of high temporal fidelity in fully turbulent flows, as the
body-force unsteadiness within the duty cycle may not be neglected.
Acknowledge
The flow field computations presented in this report is performed on the “SX-Aurora TSUBASA”
in the Cyberscience Center, Tohoku University. Supports by the center staff are very muchacknowledged. This study is also partly supported by the Grants-in-Aid for Scientific Research A(JSPS 18H03816).

References
1
2
3
4
5
6
7

M. Visbal and D. Gaitonde, AIAA paper 2006-3230 (2006).
Y. Suzen and G. Huang, AIAA paper 2005-4633 (2005).
J.-J. Wang, K.-S. Choi, L.-H. Feng, T.N. Jukes, and R.D. Whalley, Prog. Aerospace Sci. 62, 52(2013).
S. Sekimoto, T. Nonomura, and K. Fujii, AIAA J. 55, 1385 (2017).
S. Sekimoto, K. Fujii, H. Yoneda, AIAA paper 2020-0825 (2019).
E. Peers, X. Huang, and X. Luo, IEEE Trans. Plasma Sci. 37, 2250 (2009).
S. Roy, P. Zhao, A. DasGupta, and J. Soni, AIP Adv. 6, (2016).

― 24 ―

8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

32
33
34
35
36
37

SENAC Vol. 55, No. （2 2022. 4）

J.P. Boeuf and L.C. Pitchford, J. Appl. Phys. 97, 103307 (2005).
H. Nishida and T. Abe, J. Appl. Phys. 110, 013302 (2011).
C.L. Enloe, T.E. McLaughlin, R.D. VanDyken, K.D. Kachner, E.J. Jumper, T.C. Corke, M. Post,
and O. Haddad, AIAA J. 42, 595 (2004).
D.M. Orlov, Modelling and Simulation of Single Dielectric Barrier Discharge Plasma
Actuators, Thesis (Ph.D.), University of Notre Dame (2006).
F. Massines, A. Rabehi, P. Decomps, R. Ben Gadri, P. Ségur, and C. Mayoux, J. Appl. Phys. 83,
2950 (1998).
J.R. Roth, D.M. Sherman, and S.P. Wilkinson, AIAA J. 38, 1166 (2000).
W. Shyy, B. Jayaraman, and A. Andersson, J. Appl. Phys. 92, 6434 (2002).
T.C. Corke, C.L. Enloe, and S.P. Wilkinson, Annal Rev. Fluid Mech. 42, 505 (2009).
K. Asada and K. Fujii, in AIAA 5th Flow Control Conf. Reston, Virginia (2010), pp. 1–12.
K. Asada, T. Nonomura, H. Aono, M. Sato, K. Okada, and K. Fujii, Int. J. Comp. Fluid Dyn. 29,
215 (2015).
S. Kawai, T. Bouwhuis, Y. Abe, A. Yakeno, T. Nonomura, H. Aono, A. Oyama, H.W.M.M.
Hoeijmakers, and K. Fujii, Theor. Comput. Fluid Dyn. 32, 805 (2018).
S. Sekimoto, T. Nonomura, and K. Fujii, AIAA J. 55, 1385 (2017).
M. Kotsonis, S. Ghaemi, L. Veldhuis, and F. Scarano, J. Phys. D. Appl. Phys. 44, 045204 (2011).
H. Aono, S. Sekimoto, M. Sato, A. Yakeno, T. Nonomura, and K. Fujii, Mech. Eng. J. 2, 15 (2015).
K. Fujii, Appl. Sci. 8, 546 (2018).
V.I. Gibalov and G.J. Pietsch, J. Phys. D. Appl. Phys. 37, 2082 (2004).
J.P. Boeuf and L.C. Pitchford, J. Appl. Phys. 97, 103307 (2005).
A. V Likhanskii, M.N. Shneider, S.O. Macheret, R.B. Miles, A. V Likhanskii, M.N. Shneider, S.O.
Macheret, and R.B. Miles, J. Appl. Phys. 103, 13 (2008).
H. Nishida, T. Nonomura, and T. Abe, AIAA J. 54, 659 (2016).
D.V. Gaitonde, M.R. Visbal, and S. Rouy, in Proc. ASME/FEDSM2006 ( 2006), pp. 1339–1356.
D. Chen, K. Asada, S. Sekimoto, K. Fujii, and H. Nishida, Phys. Fluids 33, 037115 (2021).
URL: http://www.bolsig.laplace.univ-tlse.fr/ for more information about the BOLSIG+ database.
K. Nakai, A. Nakano, and H. Nishida, AIAA J. 1 (2021).
K. Fujii, H. Endo and M. Yasuhara, Activities of computational fluid dynamics in Japan:
compressible flow simulations, In High performance computing, 139, John Wiley and Sons Ltd.
(1992).
K. Fujii and S. Obayashi, J. Aircraft. 26, 1123 (1989).
S.K. Lele, J. Comp. Phys. 103, 16 (1992).
S. Sekimoto, K. Fujii, S. Hosokawa, and H. Akamatsu, Sensors Actuators, A Phys. 306, 111951(2020).
H. Aono, S. Kawai, T. Nonomura, M. Sato, K. Fujii, and K. Okada, AIAA J. 55, 3789 (2017).
M. Sato, K. Okada, K. Asada, H. Aono, T. Nonomura, and K. Fujii, Phys. Fluids 32, 025102 (2020).
K. Asada, Y. Ninomiya, K. Fujii, and A. Oyama, in 47th AIAA Aerospace. Sci. Meet. Reston,
Virginia, (2009), pp. 1–14.

SENAC Vol. 55, No. 2（2022. 4）― 25 ―［共同研究成果］

直交格子 法による 移動境界 問題の解 法に関す る研究高橋俊：東海大学工学部動力機械工学科
佐々木大輔：金沢工業大学工学部航空システム工学科
山田剛治 ：東海大学工学部機械工学科

Euler 方程式と断熱壁面，滑り壁面境界条件を用いた埋め込み境界法により，広範囲の速度域を安定に解析できる解析手法を構築した．極超音速流中における円柱と球の解析を実施して，先行研究の数値解析結果や，実験結果との定性的な一致を確認した．また複雑形状の例として２球から生じる衝撃波が干渉する流れ場の解析を行って，本手法の適用性について確認した．今後は移動境界問題に適用すべく本手法の拡張を行い，最終的には多数の球と衝撃波が干渉した際の運動量，エネルギーの変換や，流れ中において様々な相対移動を行う物体周りの解析と現象解明に取
り組む．

1.

はじめに

流体中を伝搬する衝撃波が物体と干渉する現象は以前から幅広く研究されており，爆発や衝撃による物体の破壊から，衝撃波が誘起する相変化[1]や，体内の結石の破砕[2, 3]など様々な問題設定の研究が行われてきた．衝撃波が物体表面を通過する時に物体に作用する瞬間的な非定常抗力
は衝撃波と物体の干渉現象において特徴的な現象の一つとして挙げられる．Tanno[4]や Sun[5]らは単一の球に衝撃波を負荷した時に生じるこの非定常抗力について実験と数値計算によって計測，調査し，非定常抗力係数は定常時の抗力係数の数倍に達することを明らかにした．またこれらの研究では非定常抗力が生じる際の衝撃波の可視化から，衝撃波が球体を通過して後流側よどみ点で収束した際に抗力が大きく低下する現象や，Re 数によっては瞬間的に抵抗が負になる現象を示した．また数値解析では Re 数依存性の議論から Re 数が大きい場合には非粘性解析でも十分に非定常抗力が議論できると結論付けられた．彼らの研究はこれまで多くの非定常抗力の研究者から引用され，現在まで世界中で詳細なメカニズムが研究されている．例としては衝撃波を微粒子に当てた際の軌道解析[6]，衝撃波管の境界層の影響なども考慮して計測を高精度化した実験的な研究[7]や，複数の粒子の流体力学的相互作用に関する CFD 研究[8, 9]などが挙げられる．これらの研究の応用の一つが粒子幕（particle curtain，particle cloud）のモデル化と現象解明である．質点モ
デルを CFD に組み入れた研究[10-13]では，質点モデルの高度化とそれが運動量交換，エネルギー
交換に及ぼす影響について検討されている．実験での粒子幕の研究には[14, 15]などが存在し，CFDによる粒子幕の研究では埋め込み境界法（Immersed Boundary Method: IBM）を用いて，乱流エネ
ルギー生成と散逸，運動量交換，エネルギー交換などが詳細に議論されてきた[16-18]．実際に非定常抗力や定常抗力を受けた微粒子は流れの中で加速度運動して流れ場も同時に変化する．IBMを用いた非粘性 CFD 解析で粒子の移動を調査した研究例も存在する[8, 19]が，CFD で粒子の運動の影響を議論した研究例は少ない．これまでは単一球周りの高 Ma 数，低 Re 数の流れ場の詳細な数値解析によるデータベース構築の研究[20-23]や，実験的研究[24]，移動物体周りの CFD による研究[25]，CFD を用いた抵抗モデルの高精度化[26]と，継続的に傾向を調査してきた．本研究ではさらなる粒子幕に関する研究の深化に向けて，1 つもしくは 2 つの粒子周りの流れ場を IBM により解析する計算手法の構築と検証計算を行った．

― 26 ―

2.

SENAC Vol. 55, No. （2 2022. 4）

計算方法

2.1 支配方程式
支配方程式には無次元化された三次元圧縮性 Euler 方程式を用いる．(1)，(2)，(3)はそれぞれ連続の式，運動量保存式，エネルギー保存式である．
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕(𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢𝑖𝑖𝑖𝑖 )+=0𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕𝑥𝑥𝑥𝑥𝑖𝑖𝑖𝑖
𝜕𝜕𝜕𝜕(𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢𝑖𝑖𝑖𝑖 ) 𝜕𝜕𝜕𝜕�𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢𝑖𝑖𝑖𝑖 𝑢𝑢𝑢𝑢𝑗𝑗𝑗𝑗 �𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕
+=−𝜕𝜕𝜕𝜕𝑥𝑥𝑥𝑥𝑗𝑗𝑗𝑗𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕𝑥𝑥𝑥𝑥𝑖𝑖𝑖𝑖
𝜕𝜕𝜕𝜕(𝜕𝜕𝜕𝜕𝜌𝜌𝜌𝜌)𝜕𝜕𝜕𝜕+
�(𝜕𝜕𝜕𝜕𝜌𝜌𝜌𝜌 + 𝜕𝜕𝜕𝜕)𝑢𝑢𝑢𝑢𝑗𝑗𝑗𝑗 � = 0𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕𝑥𝑥𝑥𝑥𝑗𝑗𝑗𝑗

(1)

(2)

(3)

ここで位置 xi は球形粒子の直径 d，流速 ui は基準音速 c0，時間 t は無次元時間 d/c0，総エネル
ギーE は基準音速の二乗 c02，密度 ρ は基準密度 ρ0，圧力 p は ρ0c02 でそれぞれ無次元化されている．また Euler 方程式内で用いる諸量を求めるための無次元化された状態方程式は以下の式で表される．
1𝜕𝜕𝜕𝜕 = (𝛾𝛾𝛾𝛾 − 1) �𝜕𝜕𝜕𝜕𝜌𝜌𝜌𝜌 − 𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢𝑖𝑖𝑖𝑖 𝑢𝑢𝑢𝑢𝑖𝑖𝑖𝑖 �2(4)

計算格子には等間隔直交格子を用いて，非粘性流束計算には SLAU スキーム[31]を用い，時間積分は 3 次精度の Total Variation Diminishing (TVD) Runge-Kutta 法を用いて行う．物体境界はレベルセット関数と埋め込み境界法を用いて表現する．これらの手法の計算精度は著者らの先行研究で既に実施済みである[27, 28]．また本研究では粒子の運動方程式との連成解析を行うべく進めている．その際，流体力の算出には先行研究により提案されたポリゴン形状を不要とする簡易な方法を用いる[29]．Navier-Stokes 方程式を用いた際の，この手法を応用した運動方程式との連成解析の精度については既に検証済みである[25, 30]．

2.2 埋め込み境界法
これまでの研究では滑り無し壁面境界条件を用いた検証結果を報告していたが，本研究では高Re 数の流れ場を模擬した解析を高効率に実施するために滑り壁面境界条件の埋め込み境界法を適用する．図 1 が滑り壁面用の埋め込み境界法の模式図で，黒，黒縁灰色．灰色，白抜き緑，黄色の丸印がそれぞれ計算セル，ゴーストセル，物体セルの中心位置，壁面上の点，イメージポイントを示す．黄色，緑，紫，赤，青の線と矢印がそれぞれ壁面からの距離，物体境界線，速度ベク
トル，速度ベクトルの壁面法線方向成分，接線方向成分である．φg，φi はそれぞれ着目しているゴーストセルと壁面上の点までの最小距離と，イメージポイントと壁面上の点までの最小距離（プローブ長さ）で，今回はイメージポイントにおける再帰参照を防止するため φi=1.75Δx（Δx は格子幅）の固定値とした．イメージポイントには trilinear 内挿で周囲の 8 点から物理量を内挿する．

直交格子法による移動境界問題の解法に関する研究

図1

𝒗𝒗𝒗𝒗𝑔𝑔𝑔𝑔,𝑛𝑛𝑛𝑛 = −― 27 ―

滑り壁のゴーストセルの速度の定義方法

𝜙𝜙𝜙𝜙𝑔𝑔𝑔𝑔
�𝒗𝒗𝒗𝒗 ∙ 𝒏𝒏𝒏𝒏�𝒏𝒏𝒏𝒏
𝜙𝜙𝜙𝜙𝑖𝑖𝑖𝑖 𝑖𝑖𝑖𝑖,𝑛𝑛𝑛𝑛(5)

𝒗𝒗𝒗𝒗𝑔𝑔𝑔𝑔,𝑡𝑡𝑡𝑡 = 𝒗𝒗𝒗𝒗𝑖𝑖𝑖𝑖,𝑛𝑛𝑛𝑛 − �𝒗𝒗𝒗𝒗𝑖𝑖𝑖𝑖,𝑛𝑛𝑛𝑛 ∙ 𝒏𝒏𝒏𝒏�𝒏𝒏𝒏𝒏(6)

𝜕𝜕𝜕𝜕𝑔𝑔𝑔𝑔 = 𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖(8)𝜌𝜌𝜌𝜌𝑔𝑔𝑔𝑔 =𝒗𝒗𝒗𝒗𝑔𝑔𝑔𝑔 = 𝒗𝒗𝒗𝒗𝑔𝑔𝑔𝑔,𝑛𝑛𝑛𝑛 + 𝒗𝒗𝒗𝒗𝑔𝑔𝑔𝑔,𝑡𝑡𝑡𝑡(7)

𝜕𝜕𝜕𝜕𝑔𝑔𝑔𝑔 = 𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖

(9)

1 𝜕𝜕𝜕𝜕𝑔𝑔𝑔𝑔 1+ �𝒗𝒗𝒗𝒗 ∙ 𝒗𝒗𝒗𝒗 �𝛾𝛾𝛾𝛾 𝛾 1 𝜕𝜕𝜕𝜕𝑔𝑔𝑔𝑔 2 𝑔𝑔𝑔𝑔 𝑔𝑔𝑔𝑔

(10)

下添え字 i，g，n，t はそれぞれイメージポイント，ゴーストセル，壁面法線方向成分，壁面接
線方向成分の諸量を指す．イメージポイントに内挿した速度 vi，密度 ρi，圧力 pi を用いて，式(5)(10)を順に計算して断熱条件と滑り壁条件を満足するようにゴーストセルの諸量である速度 vg，
密度 ρg，圧力 pg，エネルギーEg を求める．本研究では物体境界をレベルセット関数によって表現しているため，壁面との単位法線ベクトル n はレベルセット関数の勾配から容易に求められる．
これまでの研究では滑り無し壁面境界条件を適用していたが，本研究では高 Re 数の流れ場を模擬した解析を高効率に実施するために，滑り壁面境界条件の埋め込み境界法を適用した．また今
回は物体境界からプローブ長さの 2 倍程度まで物体の内部に位置しているセルをゴーストセルとした．

― 28 ―

3.

SENAC Vol. 55, No. （2 2022. 4）

結果と考察

3.1 円柱周りの解析（M∞=8.1）
構築した手法の検証として一様流マッハ数 M∞=8.1 の中に置かれた円柱周りの解析を行った．
等間隔直交格子を用いて主流方向と主流垂直方向にそれぞれ 100，200 セル，合計 20000 セルを用いて，円柱の直径を 80 セルで分割した．図 2 は用いた格子と解析により得られた等圧力線図を示しており，この結果が構造格子によって得られた先行研究[31]と良好な一致を示した．このことか
ら本解析手法の妥当性と有効性を確認した．

図2

M∞=8.1 における格子と円柱周りの等圧力線図

3.2 球周りの解析（M∞=6.7）
一様流マッハ数 M∞=6.7 の流れ中の球周りの流れ場を，直径を 40 セルと 80 セルで分割して解析を行った．図 3，図 4 はそれぞれの解析から得られた中央断面における等圧力線図と密度勾配の絶対値の分布図である．両者とも離脱衝撃波が安定して捉えられ，80 セルの方が等圧力線が滑らかで，また密度勾配の絶対値の濃淡が明瞭であった．衝撃風洞で本解析と同様の条件で可視化実験を実施した結果，衝撃波離脱距離が両者とも可視化結果とほぼ一致することを確認しており，これらの解析結果が妥当であることを確認した．

直交格子法による移動境界問題の解法に関する研究

図3

図4― 29 ―

M∞=6.7 における球周りの等圧力線図（直径 40 分割，80 分割）

M∞=6.7 における球周りの密度勾配の絶対値分布図（直径 40 分割，80 分割）

3.3 ２球周りの解析（M∞=6.7）
一様流マッハ数 M∞=6.7 の流れ中の２球周りの流れ場を，直径を 40 セルで分割して解析を行った．図 5，図 6 はそれぞれ中央断面における等圧力線図と，密度勾配の絶対値の分布図である．前方の球周りで生じた離脱衝撃波が後方の球に当たり流れ場中で最大の圧力を生じた．図 5 の等圧力分布図はこの最大圧力で無次元化して示してある．球の後流では最大圧力の 1%未満の領域も確認される広範囲の速度と圧力が生じる運動量変化が大きな流れ場であったが，解析は安定に実施できたことから本手法の堅牢性が確認できた．図 6 に示した密度勾配の絶対値の分布には２球から生じる衝撃波の干渉と衝撃波の三次元的な分布が見られた．本解析で採用した埋め込み境界法
はこのような三次元任意形状の複数物体周りの計算格子作成が容易に行えることから，これらの三次元的な流れ場構造の把握に対する適用性が確認された．

― 30 ―

SENAC Vol. 55, No. （2 2022. 4）

図5

図6

4.

M∞=6.7 における２球周りの無次元圧力分布と無次元等圧力線図

M∞=6.7 における２球周りの密度勾配の絶対値分布（中央断面，前方から見た等値面）

結論

Euler 方程式と断熱壁面，滑り壁面境界条件を用いた埋め込み境界法により，広範囲の速度域を安定に解析できる解析手法を構築した．極超音速流中における円柱と球の解析を実施して，先行研究の数値解析結果や，実験結果との定性的な一致を確認した．また複雑形状の例として２球から生じる衝撃波が干渉する流れ場の解析を行って，本手法の適用性について確認した．今後は移動境界問題に適用すべく本手法の拡張を行い，最終的には多数の球と衝撃波が干渉した際の運動量，エネルギーの変換や，流れ中において様々な相対移動を行う物体周りの解析と現象解明に取り組む．

直交格子法による移動境界問題の解法に関する研究― 31 ―

謝辞
本研究は東北大学サイバーサイエンスセンターのスーパーコンピュータを利用することで実現することができた．また研究にあたって同センター関係各位に有益なご指導とご協力をいただいた．この場をお借りして謝意を表する．
また本解析手法の開発において，東京大学大学院工学系研究科航空宇宙工学専攻航空宇宙
システム学講座 今村太郎 准教授には貴重な意見をいただいた．ここに謝意を表する．

参考文献[1] Duvall, G. E., et al. “Phase transitions under shock-wave loading.” Reviews of Modern Physics 49.3
(1977): 523.
[2] Chaussy, C., et al. “First clinical experience with extracorporeally induced destruction of kidney stones
by shock waves.” The Journal of urology 127.3 (1982): 417-420.
[3] Takayama, K., et al. “Applications of shock wave research to medicine.” WIT Transactions on Modellingand Simulation 41 (2005).
[4] Tanno, H., et al. “Interaction of a shock with a sphere suspended in a vertical shock tube.” Shock Waves
13.3 (2003): 191-200.
[5] Sun, M., et al. “Unsteady drag on a sphere by shock wave loading.” Shock waves 14.1 (2005): 3-9.
[6] Jourdan, G., et al. “Drag coefficient of a sphere in a non-stationary flow: new results.” Proceedings of theRoyal Society A: Mathematical, Physical and Engineering Sciences 463.2088 (2007): 3323-3345.
[7] Wagner, Justin L., et al. “Shock tube investigation of quasi-steady drag in shock-particle interactions.”
Physics of Fluids 24.12 (2012): 123301.
[8] Sridharan, P., et al. “Shock interaction with one-dimensional array of particles in air.” Journal of Applied
Physics 117.7 (2015): 075902.
[9] Mehta, Y., et al. “Numerical investigation of shock interaction with one-dimensional transverse array ofparticles in air.” Journal of Applied Physics 119.10 (2016): 104901.
[10] Ling, Y., et al. “Interaction of a planar shock wave with a dense particle curtain: Modeling andexperiments.” Physics of Fluids 24.11 (2012): 113301.
[11] Regele, J. D., et al. “Unsteady effects in dense, high speed, particle laden flows.” International Journalof Multiphase Flow 61 (2014): 1-13.
[12] Jiang, L.-J., et al. “Particle-resolved simulations of shock-induced inviscid flow through particle-curtain
at initial stage.” Computers & Fluids (2021): 105196.
[13] Shallcross, G. S., et al. “A volume-filtered description of compressible particle-laden flows.”
International Journal of Multiphase Flow 122 (2020): 103138.
[14] Bordoloi, A. D., et al. “Relaxation drag history of shock accelerated microparticles.” Journal of Fluid
Mechanics 823 (2017).
[15] DeMauro, E. P., et al. “Unsteady drag following shock wave impingement on a dense particle curtain
measured using pulse-burst PIV.” Physical review fluids 2.6 (2017): 064301.
[16] Das, P., et al. “A sharp interface Cartesian grid method for viscous simulation of shocked particle-laden
flows.” International Journal of Computational Fluid Dynamics 31.6-8 (2017): 269-291.
[17] Hosseinzadeh-Nik, Z., et al. “Investigation and quantification of flow unsteadiness in shock-particle
cloud interaction.” International Journal of Multiphase Flow 101 (2018): 186-201.
[18] Osnes, A. N., et al. “Computational analysis of shock-induced flow through stationary particle clouds.”
International Journal of Multiphase Flow 114 (2019): 268-286.
[19] Luo, K., et al. “Studies on shock interactions with moving cylinders using immersed boundary method.”
Physical Review Fluids 2.6 (2017): 064302.

― 32 ―

SENAC Vol. 55, No. （2 2022. 4）[20] Nagata, T., et al. “Investigation on subsonic to supersonic flow around a sphere at low Reynolds numberof between 50 and 300 by direct numerical simulation.” Physics of Fluids 28.5 (2016): 056101.
[21] Nagata, T., et al. "Direct numerical simulation of flow past a transversely rotating sphere up to a
Reynolds number of 300 in compressible flow." Journal of Fluid Mechanics 857 (2018): 878-906.
[22] Nagata, T., et al. “Direct numerical simulation of subsonic, transonic and supersonic flow over an
isolated sphere up to a Reynolds number of 1000.” Journal of Fluid Mechanics 904 (2020).
[23] Nagata, T., et al. “Direct numerical simulation of flow around a heated/cooled isolated sphere up to a
Reynolds number of 300 under subsonic to supersonic conditions.” International Journal of Heat andMass Transfer 120 (2018): 284-299.
[24] Nagata, T., et al. “Experimental investigation of transonic and supersonic flow over a sphere for
Reynolds numbers of 103–105 by free-flight tests with schlieren visualization.” Shock Waves 30 (2020):
139-151.
[25] Mizuno, Y., et al. “Direct numerical simulation of gas–particle flows with particle–wall collisions using
the immersed boundary method.” Applied Sciences 8.12 (2018): 2387.
[26] Loth, E., et al. “Supersonic and Hypersonic Drag Coefficients for a Sphere.” AIAA Journal 59.8 (2021):
3261-3274.
[27] Takahashi, S., et al. “A numerical scheme based on an immersed boundary method for compressible
turbulent flows with shocks: application to two-dimensional flows around cylinders.” Journal of Applied
Mathematics 2014 (2014): 252478.
[28] Mizuno, Y., et al. “A simple immersed boundary method for compressible flow simulation around a
stationary and moving sphere.” Mathematical Problems in Engineering 2015 (2015): 438086.
[29] Nonomura, T., et al. “A comparative study on evaluation methods of fluid forces on Cartesian grids.”
Mathematical Problems in Engineering 2017 (2017): 8314615.
[30] Nagata, T., et al. “A simple collision algorithm for arbitrarily shaped objects in particle‐resolved flow
simulation using an immersed boundary method.” International Journal for Numerical Methods in Fluids
92.10 (2020): 1256-1273.
[31] Kitamura, K., et al. “Towards shock-stable and accurate hypersonic heating computations: A newpressure flux for AUSM-family schemes.” Journal of Computational Physics 245 (2013): 62-83.

SENAC Vol. 55, No. 2（2022. 4）― 33 ―[お知らせ]

コンパイ ラのバー ジョンア ップにつ いて
2022 年 4 月 11 日に AOBA のコンパイラをバージョンアップいたしました。
システム

コンパイラ名

旧バージョン

新バージョ

ドキュメントン
AOBA-A

AOBA-BSX-AT マニュアル

Fortran Compiler

3.3.0

3.4.0

C/C++ Compiler

3.3.0

3.4.0

MPI※1

2.18.03.0.0

AOCC Compiler

3.1

3.2

リリースノート＆マニュアル

Intel Compiler※2

oneAPI 2021.3

oneAPI 2022

oneAPI マニュアル関連※1 MPI を利用するプログラムは再コンパイルが必要※2 Intel oneAPI 2022 の環境変数設定ファイルは、bash 向けのみの提供（共同利用支援係，共同研究支援係）【重要】AOBA-A で MPI ライブラリを利用するプログラムについての注意点
2022 年 4 月のサービス再開後から AOBA-A の Fortran コンパイラ、C/C++コンパイラおよび MPI ライ
ブラリをバージョンアップします。今回のバージョンアップでは、MPI ライブラリを利用したプログラ
ムは再コンパイルと MPI ライブラリのリンクが必須となります。バージョンアップ前にコンパイルされた MPI ライブラリを利用したプログラムは、正常に実行できませんので必ず再コンパイル・MPI ライブ
ラリのリンクを行って下さい。なお、コンパイラコマンド名、コンパイルオプション等に変更はありません。
再コンパイルが必要なプログラム
サブシステム AOBA-A（SX-Aurora TSUBASA）で実行する MPI ライブラリを利用したプログラム
コマンド名 mpinfort でコンパイルする Fortran プログラム および
コマンド名 mpincc/mpinc++でコンパイルする C/C++プログラムの全てが対象になります。
（共同利用支援係，共同研究支援係）

― 34 ―

SENAC Vol. 55, No. （2 2022. 4）[お知らせ]

令和 4 年度 サイバーサ イエンス センター 講習会計 画
東北大学サイバーサイエンスセンターでは、以下の講習会の開催を予定しています（zoom によるオンライ開催）。多くの皆様のご参加をお持ちしております。申し込み方法等、詳しくは以下をご覧ください。
https://www.ss.cc.tohoku.ac.jp/lectures/
No.

講習会名

開催日時

募集
講師
人数

1

はじめての Linux

5 月 24 日(火)13:30-15:30

15

2

はじめてのスパコン

5 月 26 日(木)13:30-15:30

15

3

はじめての並列化

6 月 6 日(月)13:30-15:30

15

6 月 24 日(金)13:30-16:30

15

ネットワークとセキ
ュリティ入門

8 月 3 日(水)13:30-15:30

上限

水木

なし

(サイバーサイエンスセンター)

Gaussian 入門

8 月 25 日(木)13:30-16:30

15

岸本

8 月 30 日(火)13:30-16:30

15

9 月 9 日(金)13:30-16:30

15

10 月 18 日(火)13:30-15:30

15

10 月 20 日(木)13:30-15:30

15

内容

小野

・Linux システムの基本的な(情報部情報基盤
使い方・エディタの使い方課)
山下・スーパーコンピュータの紹介と(情報部情報基盤
利用法入門
課)

4

5

6

7

8

9

10

MATLAB 入門

Mathematica 入門

Fortran 入門

はじめての Linux

はじめてのスパコン

小松・並列プログラミングの概要(サイバーサイエンスセンター)

陳・MATLAB の基本的な使い方(秋田県立大学)

・ネットワークの基本的な仕組み
・ネットワークの危険性と安全対策
・Gaussian の基本的な使い方(理学研究科)横井・Mathematica の基本的な使い方(尚絅学院大学)

田口・Fortran の入門編

(日本原子力研
究開発機構)

大泉・Linux システムの基本的な使い方・エディタの使い方
佐々木・スーパーコンピュータの紹介と(情報部情報基盤 利用法入門(情報部情報基盤課)

課)

11

12

13

SX-Aurora TSUBASA の性能分析・高速化

10 月 25 日(火)13:30-16:30

15並列プログラミング
入門Ⅰ（OpenMP）

10 月 27 日(木)13:30-16:30

15並列プログラミング
入門Ⅱ（MPI）

10 月 31 日(月)13:30-16:30

15江川(サイバーサイエンスセンター)小松(サイバーサイエンスセンター)下村(サイバーサイエンスセンター)・スーパーコンピュータでの性能解析から最適化まで・並列プログラミングの概要
・OpenMP による並列プログラミングの基礎
・利用法
・MPI による並列プログラミングの
基礎
・利用法（共同利用支援係）

SENAC Vol. 55, No. 2（2022. 4）― 35 ―［大規模科学計算システム］

鍵ペアの作成とログイン方法
情報部情報基盤課 共同研究支援係 共同利用支援係

1. はじめに本センターのシステムは、セキュリティ対策として、公開鍵認証方式による SSH 接続を採用しています。また、フロントエンドサーバは、ログインサーバを経由しなければログインできない構成としています。
本稿では、公開鍵認証方式で使用する鍵ペアの作成と各サーバのログイン方法についてご紹介します。解説では以下のターミナルソフトを使用する例をご紹介します。
（Windows の場合）

Windows PowerShell（macOS／Linux の場合）

ターミナル

本センターのシステムをはじめて利用する方は、以下の手続きが必要です。
(1) 利用者番号の取得（ 利用申請：https://www.ss.cc.tohoku.ac.jp/apply-for-use/ ）(2) 鍵ペアの作成（4 章）
以前のシステムを利用していた方は、(1)(2)の手続きは不要です。以前使用していた利用者番号および鍵ペアをそのままご利用いただけます。5 章からお読みください。

2. ログイン認証方式
表 1 に、各サーバのログイン認証方式を示します。
表 1 各サーバのログイン認証方式
サーバ名
ログインサーバ

用途
フロントエンドサーバの入口（踏み台サーバ）

フロントエンド

計算機の利用

サーバ（コンパイル、ジョブ実行、等）

データ転送サーストレージシステムとの大容量

バのデータ転送HPCI 用ログイ

HPCI、HPCI-JHPCN ユーザ専

ンノード

用ログインノード-

センター内施設の利用（※2）

ログインホスト名

認証方式

login.cc.tohoku.ac.jp

公開鍵

（※1）

公開鍵またはパスワード

file.cc.tohoku.ac.jp

公開鍵

hpcif.cc.tohoku.ac.jp

公開鍵

-

パスワード

― 36 ―

SENAC Vol. 55, No. （2 2022. 4）（※1） フロントエンドサーバは、ログインサーバからしかログインできません。本稿では
多段 SSH による接続方法を解説します。
（※2） 本センター内の施設（大判カラープリンタ、利用者端末、講習会端末）はパスワー
ド認証でご利用いただけます。利用にあたり、秘密鍵を持参する必要はありません。

3. 鍵ペアの作成からログインまでの流れ
図 1 に、鍵ペア作成からログインまでの流れを示します。①は 4 章、②③は 5 章で詳し
く解説します。

図 1 鍵ペア作成からログインまでの流れ① 鍵ペアの作成（初回ログイン時、および、ログイン端末を追加する場合）
利用者ポータルで鍵ペアを作成します。作成された秘密鍵は、利用者のローカル PC に保存します。公開鍵は、スーパーコンピュータ AOBA のホームディレクトリ上に自動で保存されます。

鍵ペアの作成とログイン方法②― 37 ―

ターミナルソフトの設定（初回ログイン時）
各サーバにログインするための設定を行います。使用するターミナルソフトによって必要な設定は異なります。
③

ログイン
利用者のローカル PC に保存した秘密鍵を使ってログインします。フロントエンドサーバ

は、ログインサーバを経由して多段 SSH でログインします。

4. 公開鍵認証方式で使用する鍵ペアの作成
4.1. 公開鍵認証方式を使用する上での注意事項
以下の注意事項を必ず守ってください。守らない場合、不正アクセス（不正ログイン、
クライアントのなりすまし、暗号化された通信の暴露、他サーバへの攻撃、等）のリスクが非常に高まり、大変危険です。ご注意願います。
・ パスフレーズなしの秘密鍵を使用しないこと
・ 秘密鍵、パスフレーズを使いまわさないこと
・ 秘密鍵を持ち出さないこと（メールに添付しない、USB メモリ等に保存しない）
・ 秘密鍵をスーパーコンピュータ AOBA のホームディレクトリに保存しないこと
・ 公開鍵と秘密鍵の鍵ペアを同一ノード上に保存しないこと

4.2. 鍵ペアの作成（初回ログイン時、および、ログイン端末を追加する場合）〇初回ログイン時
鍵ペアの作成は、利用者ポータルで行います。
(1) 以下の URL 先から利用者ポータルを開きます。
利用者ポータルには、利用者番号とパスワード（※）でログインします。
利用者ポータル：https://www.ss.cc.tohoku.ac.jp/portal/(2) 「SSH 公開鍵登録」ボタンをクリックします。
(3) 利用者ポータルの画面の説明に従い、鍵ペアを作成します。
（パスフレーズを設定し、鍵生成・登録ボタンをクリック）(4) 作成された秘密鍵を利用者のローカル PC に保存します。保存先は以下を推奨します。
フォルダがない場合は新規作成します。
（Windows の場合）

C:¥Users¥ユーザ¥.ssh（macOS／Linux の場合） ~/.ssh
公開鍵は、ホームディレクトリ（~/.ssh/authorized_keys）に自動で保存されます。
※利用者ポータルで使用するパスワードの変更方法は、6 章を参照してください。

― 38 ―

SENAC Vol. 55, No. （2 2022. 4）〇別の PC からログインする場合（ログイン端末を追加する場合）
既存の秘密鍵を使いまわすのではなく、ログイン端末ごとに鍵ペアを作成してください。
初回ログイン時と同じ手順で、新しい鍵ペアを追加します。

5. 各サーバのログイン方法
5.1. ターミナルソフトの設定（初回ログイン時）
利用者のローカル PC 上で、ターミナルソフトの設定を行います。
以降の解説は、次のフォルダを「.ssh フォルダ」と呼び、秘密鍵を「id_rsa_cc」というファイル名で.ssh フォルダに保存した場合とします。
（Windows の場合）

C:¥Users¥ユーザ¥.ssh（macOS／Linux の場合）~/.ssh

各ログインホストのホスト名は、次の文字列で設定するものとして解説します。ホスト
名には任意の文字列を設定することができます（他の設定との重複は不可）。他の文字列を設定した場合は、以降の解説におけるホスト名を読み替えてください。
（ログインサーバ）

login（フロントエンドサーバ）

front（データ転送サーバ）

file（HPCI 用ログインノード） hpcif(1) macOS／Linux の場合は、秘密鍵のパーミッションの変更（600 に設定）が必要です。
ターミナルソフトを起動し、以下のコマンドを実行します。
$ chmod 600 ~/.ssh/id_rsa_cc
以降は Windows、macOS／Linux 共通です。
(2) .ssh フォルダの「config」というファイルをテキストエディタで開きます。ファイルがない場合は新規作成します。拡張子はつけません。
（フォルダの設定を「拡張子を表示しない」にしている場合、意識せずに拡張子つきのファイルを作成している可能性があります。config ファイルに拡張子がついていると、ログインできません。ご注意ください）(3) config ファイルに以下の設定を記述します。太
太字下線の部分は、ご自身の環境に合わせて読み替えてください。

鍵ペアの作成とログイン方法― 39 ―〇フロントエンドサーバを利用するための設定（※）# ログインサーバの設定（ホスト名を”login”とする場合）
Host login# ホスト名を指定

HostName login.cc.tohoku.ac.jp# ログインホスト名を指定

User 利用者番号# 利用者番号を指定

IdentityFile ~/.ssh/id_rsa_cc# 秘密鍵の保存場所とファイル名を指定# フロントエンドサーバの設定（ホスト名を”front”とする場合）
Host front
HostName front.cc.tohoku.ac.jp
User 利用者番号
ProxyCommand ssh -CW %h:%p login# login 経由で多段 SSH する設定

IdentityFile ~/.ssh/id_rsa_cc〇データ転送サーバを利用するための設定# データ転送サーバの設定（ホスト名を”file”とする場合）
Host file
HostName file.cc.tohoku.ac.jp
User 利用者番号
IdentityFile ~/.ssh/id_rsa_cc〇HPCI 用ログインノードを利用するための設定# HPCI 用ログインノードの設定（ホスト名を”hpcif”とする場合）
Host hpcif
HostName hpcif.cc.tohoku.ac.jp
User 利用者番号
IdentityFile ~/.ssh/id_rsa_cc（※）Windows の場合、フロントエンドサーバへのログイン時に以下のようなエラーが出る場合があります。
$ ssh front
CreateProcessW failed error:2
posix_spawn: No such file or directory

― 40 ―

SENAC Vol. 55, No. （2 2022. 4）

エラーが出た場合は次の要領で config ファイルを書き換えてください。
[1] ターミナルソフトを起動し、以下のコマンドで ssh の絶対パスを調べる。
$ gcm ssh
CommandType

Name

Version Source----------------------------

Application

ssh.exe

x.x.x

C:¥WINDOWS¥System32¥OpenSSH¥ssh.exe[2] config ファイルの「ProxyCommand ssh …」の行の「ssh」の部分を、絶対パス（[1]で「Source」に表示された文字列）に書き換える。
# 修正前
ProxyCommand ssh -CW %h:%p login# 修正後ProxyCommand C:¥WINDOWS¥System32¥OpenSSH¥ssh.exe -CW %h:%p login

5.2. フロントエンドサーバのログイン方法
ターミナルソフトを起動し、以下のコマンドを実行するとログインします。ホスト名を別の文字列で設定している場合は「front」の部分を読み替えてください。
$ ssh front
フロントエンドサーバは冗長構成になっており、自動的に front1 または front2 が選択されます。どちらにログインしても、動作は変わりません。
なお、フロントエンドサーバでは一定時間以上のプロセスは実行できません。また、大容量のデータ転送はシステムに高い負荷がかかります。大容量のデータ転送を行う場合は、データ転送サーバをご利用ください。

5.3. データ転送サーバの利用方法
データ転送サーバは、ログインして利用するのではなく、利用者のローカル PC 上から
scp コマンドや sftp コマンドで利用します。詳しくは以下をご参照ください。
データ転送（ストレージ）：https://www.ss.cc.tohoku.ac.jp/storage/

5.4. HPCI 用ログインノードのログイン方法
ターミナルソフトを起動し、以下のコマンドを実行するとログインします。ホスト名を別の文字列で設定した場合は「hpcif」の部分を読み替えてください。
$ ssh hpcif

鍵ペアの作成とログイン方法― 41 ―

5.5. ログインシェルの確認と変更
ログインシェルは、デフォルトでは csh が設定されています。設定の確認および変更は以下の手順で行います。ログインシェルの変更がシステム全体に反映されるまで、15 分程度かかります。
(1) フロントエンドサーバにログインする。
(2) 以下のコマンドを実行する。
〇ログインシェルの確認
front1 $ fchsh （ログインシェルの確認）
Enter Password: （パスワードを入力）
loginShell: /bin/tcsh（現在のログインシェルが表示される）〇ログインシェルの変更
front1 $ fchsh /bin/bash （ログインシェルを/bin/bash に変更）
Enter Password: （パスワードを入力）
Changed loginShell to /bin/bash（ログインシェルが変更された）

6. パスワードの変更
利用者ポータルなどで使用するパスワードの変更は、以下の手順で行います。
(1) 以下の URL 先から利用者ポータルを開きます。
利用者ポータルには、利用者番号とパスワードでログインします。
利用者ポータル：https://www.ss.cc.tohoku.ac.jp/portal/(2) 「パスワード変更」ボタンをクリックします。
(3) 利用者ポータルの画面の説明に従い、新しいパスワードを設定します。
(4) 以下で使用するパスワードが変更されます。
・利用者ポータルへのログインパスワード
・大判カラープリンタのプリンタサーバへのログイン
・ログインシェルの変更時のパスワード

7. おわりに
本稿では、鍵ペアの作成とログイン方法についてご紹介しました。センターのシステム
を安全にご利用いただければ幸いです。ご不明な点、ご質問等ございましたら、お気軽にセンター（利用相談）までお問い合わせください。
利用相談：https://www.ss.cc.tohoku.ac.jp/consultation/
また、センターからのお知らせは、ウェブサイトにてご確認ください。
センターウェブサイト：https://www.ss.cc.tohoku.ac.jp/

― 42 ―

SENAC Vol. 55, No. （2 2022. 4）［大規模科学計算システム］

ストレージシステムの利用法
ストレージシステムの利用法
ストレージシステムの利用法

情報部情報基盤課 共同研究支援係 共同利用支援係
情報部情報基盤課 共同研究支援係 共同利用支援係
2020 年 10 月 1 日 初版2020 年 10 月 1 日 初版

1 はじめに1 はじめに2022 年 5 月 13 日 第 3 版
2022 年 5 月 13 日 第 3 版

本稿では，スーパーコンピュータ AOBA のストレージシステムの利用法について紹介します．スト
本稿では，スーパーコンピュータ AOBA のストレージシステムの利用法について紹介します．スト
レージ環境にあるホームディレクトリと課題領域の容量確認方法と実行した結果ファイル等のデータ
レージ環境にあるホームディレクトリと課題領域の容量確認方法と実行した結果ファイル等のデータ(ストレージ環境) をローカル PC へ転送する方法およびローカル PC からストレージ環境へ転送する方(法について説明します．
ストレージ環境) をローカル PC へ転送する方法およびローカル PC からストレージ環境へ転送する方法について説明します．

2 ストレージ環境
2 ストレージ環境
2.1
2.1

ホームディレクトリ (uhome)
ホームディレクトリ (uhome)

プログラムファイル等を置く自分専用のホームディレクトリになり，ScaTeFS マウントし，AOBA-A
ScaTeFS マウントし，AOBA-Aとプログラムファイル等を置く自分専用のホームディレクトリになり，
AOBA-B の両方に共有しています．
と AOBA-B の両方に共有しています．

ディレクトリ名： /uhome/利用者番号/uhome/
利用者番号
ディレクトリ名：
クォータ（容量）制限：
5TB
クォータ（容量）制限： 5TB

クォータ制限を超過した場合，新規の書き込みができなくなりますのでご注意ください． クォータ制
クォータ制限を超過した場合，新規の書き込みができなくなりますのでご注意ください． クォータ制
限を下回るように容量を削減すれば再度書き込みが可能になります．
限を下回るように容量を削減すれば再度書き込みが可能になります．
• ホームディレクトリの容量確認コマンド• ホームディレクトリの容量確認コマンド
front$ uquota

front$ uquota
表示例

表示例
Disk quotas for user 利用者番号

Diskquotas for user 利用者番号
Filesystemused(KB)
Filesystemused(KB)/uhome/利用者番号4/uhome/利用者番号

2.2
2.2

4

quota(TB)
quota(TB)
5
5

ホームディレクトリの容量追加申請については，2.3 章をご参照ください．
ホームディレクトリの容量追加申請については，2.3 章をご参照ください．

課題領域（/short/プロジェクトコード)
課題領域（/short/プロジェクトコード)

課題領域は事前申請となり，同一プロジェクトコードの利用者間で大規模なデータ容量を利用され課題領域は事前申請となり，同一プロジェクトコードの利用者間で大規模なデータ容量を利用される領域になります．申請を希望される際は，ストレージ資源の兼ね合いもありますので，「利用相談」「利用相談」（る領域になります．申請を希望される際は，ストレージ資源の兼ね合いもありますので，https://www.ss.cc.tohoku.ac.jp/consultation/）から事前にご連絡をお願いします．
（https://www.ss.cc.tohoku.ac.jp/consultation/
容量については，申請されたディスク容量）から事前にご連絡をお願いします．
(quota 値)になります．ホームディレクトリ同様に AOBA-(quota 値)になります．ホームディレクトリ同様に AOBAと AOBA-B の両方に共有しています．
A 容量については，申請されたディスク容量
A と AOBA-B の両方に共有しています．

1
1

ストレージシステムの利用法― 43 ―

課題領域の容量については，以下のコマンドの quota(TB) 部分をご確認お願いします．

• 課題領域の容量確認コマンド

front$ uquota -A プロジェクトコード
表示例

Disk quotas for project プロジェクトコード
Filesystem/short/プロジェクトコード

used(KB)

quota(TB)

1020• 利用方法– 同一プロジェクトコードの利用者間でデータを共有する．

– 対象ディレクトリ内で利用者がそれぞれの専用サブディレクトリ (パーミッション：700) を作成し使用する．

利用方法の一例

プロジェクトコード：xx200001 の場合1 ）プロジェクトコードの利用者間で share を作成しデータを共有
コマンド例

front$ cp ホームディレクトリデータ /short/xx200001/share/2 ） 利用者の専用サブディレクトリを作成し使用
コマンド例

front$ mkdir /short/xx200001/利用者番号

front$ chmod 700 /short/xx200001/利用者番号

front$ cp ホームディレクトリデータ /short/xx200001/利用者番号/【留意事項】– ファイル同期コマンド（rsync コマンド），コピーコマンド (cp コマンド) を使用する際は，グ

ループ権限を保持するオプションは設定せずにご利用ください．オプションを付けた場合，同一プロジェクトコード間のグループによる容量制限で正しく管理できなくなる恐れがあります．また，移動コマンド（mv コマンド）によるファイルの移動を行った場合，元のファ

イルのグループ権限が保持されてしまいますので，rsync コマンド，cp コマンドを利用するようにしてください．

– 課題領域を当年度までのご利用の際，翌年度はデータ保管を行っていません．猶予期間後，

対象課題領域を削除しますので，ローカル PC のディスクへ移行を速やかに進めてください．

2.3 ストレージ申請

ファイル容量の追加は 1TB 単位から申請可能です．ホームディレクトリ，課題領域ともに利用負担金が発生しますので，詳しくは「利用負担金」(https://www.ss.cc.tohoku.ac.jp/charge/) をご参照ください．
申請用紙は，「ストレージ容量申請書」(https://www.ss.cc.tohoku.ac.jp/application-form/) を使い申請をお願いします．

2

― 44 ―

SENAC Vol. 55, No. （2 2022. 4）

3 データ転送方法

デ ー タ 転 送 サ ー バ へ ロ グ イ ン し ，SSH に よ る 暗 号 化 を 行 う scp(Secure CoPy)，SFTP(SshFile Transfer Protocol) を 利 用 し ま す ．接 続 方 法 に つ い て は「 利 用 申 請 か ら ロ グ イ ン ま で 」(https://www.ss.cc.tohoku.ac.jp/ﬁrst-use/) をご参照ください．

図 1 アクセスイメージ

3.1 Powershell(Windows)・MAC・Linux

標準で SSH クライアントがインストールされています．インストールされている各 Terminal ソフトでデータ転送サーバへログインします．

• scp コマンド

SSH 利用し，ネットワーク・ホスト間でファイルを安全にコピーするためのコマンドです．

ローカル PC からリモート (ストレージ環境) に転送$ scp -i 秘密鍵ファイル ローカル PC のファイル名 利用者番号@ﬁle.cc.tohoku.ac.jp:リモートの保存先パス

scp コマンドの実行例（ローカル PC 上にある sample.txt ファイルをホームディレクトリへ転送)$ scp -i ~/.ssh/id_rsa_cc sample.txt 利用者番号@ﬁle.cc.tohoku.ac.jp:sample.txt
パスフレーズを聞かれますので入力します．

リモート (ストレージ環境) からローカル PC に転送$ scp -i 秘密鍵ファイル 利用者番号@ﬁle.cc.tohoku.ac.jp:ファイル名 ローカル PC の保存先パス
scp コマンドの実行例（ホームディレクトリ上にある sample.txt をローカル PC へ転送）

$ scp -i ~/.ssh/id_rsa_cc 利用者番号@ﬁle.cc.tohoku.ac.jp:sample.txt ./詳しい用例については man コマンド等を利用し，scp コマンドのマニュアル閲覧をお願いします．

$ man scp3

ストレージシステムの利用法― 45 ―• sftp コマンド

SSH 利用し，対話的なファイル転送を行うことができるコマンドです．

ローカル PC からリモート (ストレージ環境) に転送$ sftp -i 秘密鍵ファイル 利用者番号@ﬁle.cc.tohoku.ac.jp
パスフレーズを聞かれますので入力します．

sftp> と表示されたら成功です．

sftp> put ファイル名 保存先フォルダ
sftp コマンドの実行例（ローカル PC 上にある sample.txt ファイルをホームディレクトリへ転送)

sftp> put パス名/sample.txt ./リモート (ストレージ環境) からローカル PC に転送$ sftp -i 秘密鍵ファイル 利用者番号@ﬁle.cc.tohoku.ac.jp
パスフレーズを聞かれますので入力します．

sftp> と表示されたら成功です．

sftp> get ファイル名 保存先フォルダ
sftp コマンドの実行例（ホームディレクトリ上にある sample.txt をローカル PC へ転送）

sftp> get ./sample.txt ./詳しい用例については man コマンド等を利用し，sftp コマンドのマニュアル閲覧をお願いします．

$ man sftp

3.2 WinSCP(Windows ソフト)

標準で scp，sftp に対応したソフトウェアがインストールされていないため，はじめにインストールする必要があります．

ここでは，代表的なソフトウェアである WinSCP を利用したファイル転送方法を説明します．

図 2 WinSCP 設定画面

4

― 46 ―

SENAC Vol. 55, No. （2 2022. 4）

1. WinSCP を起動します．

2. ホスト名 （上図（1））に ﬁle.cc.tohoku.ac.jp と入力します．

3. ユーザ名 （上図（2））に利用者番号を入力します．

4. 設定 （上図（3））のプルダウンメニューから設定 （上図（4））をクリックします．

図 3 WinSCP 鍵認証設定画面

図 4 WinSCP 鍵認証画面

5. 左側ナビゲーションメニューの認証（上図（5））を選択します．

6. 秘密鍵のプルダウンメニュー （上図（6））をクリックし，ログインに使用する秘密鍵を指定します．秘密鍵 (.ppk ファイル) を未生成の場合，3.2.1 章をご参照ください．

7. OK （上図（7））をクリックし，鍵の設定を保存します．

8. 元の画面に遷移しますので，ログインをクリックしてください．

9. パスフレーズの入力画面が出ますので，パスフレーズを入力（上図（8））した上で， OK （上図（9））をクリックしてください．成功しますと WinSCP の画面が表示されファイル転送が可能になります．

3.2.1 WinSCP 用の鍵生成手順

1. WinSCP を起動した後，「ツール」をクリックし，「PuTTYgen を実行」を選択します．

2. PuTTYgen を起動すると，「PuTTY Key Generator」ダイアログボックスが表示されますので，「Load」をクリックします．

図 5 WinSCP 初期設定画面

図 6 PuTTY Key Generator 画面

5

ストレージシステムの利用法― 47 ―

3. ファイルの選択画面が表示されますので，ポータルサイトで作成した秘密鍵「id_rsa_cc」を選択
すると，パスフレーズの入力を求められます．

4. パスフレーズの内容が一致すると Notice(情報) メッセージが表示されるので，OK をクリックします．

5.「Save private key」をクリックし，ファイル名を設定します．

6. 設定が完了しましたら，「PuTTY Key Generator」の画面は閉じてください．

4 おわりに本 稿 で は ，ス ト レ ー ジ シ ス テ ム の 利 用 法 を 紹 介 し ま し た ．ご 不 明 な 点 ，ご 質 問 等 ご ざ い ま

し た ら ，お 気 軽 に セ ン タ ー ま で お 問 い 合 わ せ く だ さ い ．問 い 合 わ せ 先 に つ い て は「 利 用 相 談 」（https://www.ss.cc.tohoku.ac.jp/consultation/）をご参照下さい．

6

― 48 ―

SENAC Vol. 55, No. （2 2022. 4）［大規模科学計算システム］

アプリケーションサービスの紹介
アプリケーションサービスの紹介
アプリケーションサービスの紹介

1 はじめに1 はじめに情報部情報基盤課 共同利用支援係
情報部情報基盤課 共同利用支援係

本センター大規模科学計算システムでは、分子軌道計算、数式処理、データ処理等の各アプリケー
本センター大規模科学計算システムでは、分子軌道計算、数式処理、データ処理等の各アプリケー
ションソフトウェアをサービスしています。本稿では、スーパーコンピューター
AOBA でサービスを

ションソフトウェアをサービスしています。本稿では、スーパーコンピューター
AOBA でサービスを行っているアプリケーションソフトウェアの紹介をします。
行っているアプリケーションソフトウェアの紹介をします。
アプリケーションソフトウェアの紹介は、以下 URL の本センター大規模科学計算システム Web サイ

アプリケーションソフトウェアの紹介は、以下 URL の本センター大規模科学計算システム Web サイトにも掲載しています。
トにも掲載しています。
■アプリケーションサービス https://www.ss.cc.tohoku.ac.jp/software-service/■アプリケーションサービス https://www.ss.cc.tohoku.ac.jp/software-service/
表 1 アプリケーションソフトウェアとサービスホスト
表 1 アプリケーションソフトウェアとサービスホスト
アプリケーションソフトウェア
アプリケーションソフトウェア
非経験的分子軌道計算プログラム
Gaussian16
Gaussian16
非経験的分子軌道計算プログラム
反応経路自動探索プログラム
GRRM17 ※ 1

反応経路自動探索プログラム
統合型数値計算ソフトウェア
統合型数値計算ソフトウェア
数値解析ソフトウェア
数値解析ソフトウェア
熱流体解析ソフトウェア
熱流体解析ソフトウェア
第一原理計算統合パッケージ

第一原理計算統合パッケージ
分子動力学シミュレーション
分子動力学シミュレーション
古典分子動力学計算プログラム
古典分子動力学計算プログラム
コンテナプラットフォーム
コンテナプラットフォーム
プログラム対話型実行環境

GRRM17※ 1※ 2Mathematica
Mathematica
MATLAB※※
1 2
MATLAB ※ 1
OpenFOAM
OpenFOAM
Quantum Espresso
Quantum
Espresso
GROMACSGROMACSLAMMPSLAMMPS
Singularity
Singularity
Jupyter Notebook
Jupyter Notebook

プログラム対話型実行環境
RIST が整備するアプリケーションソフトウェア
RIST が整備するアプリケーションソフトウェア※ 1 学術目的での利用（非商用利用）のみ可能です。
※
2 学術目的での利用（非商用利用）のみ可能です。
東北大学の構成員のみ利用可能です。
※1

サービスホスト
サービスホスト

AOBA-B
AOBA-B
AOBA-B
フロントエンドサーバ

フロントエンドサーバAOBA-B
フロントエンドサーバおよびフロントエンドサーバおよびAOBA-B
AOBA-B

AOBA-AAOBA-B
および AOBA-B
AOBA-AAOBA-B
および AOBA-B
AOBA-B
AOBA-B
AOBA-B
フロントエンドサーバ

フロントエンドサーバ
AOBA-A
および AOBA-B
AOBA-A および AOBA-B※ 2 東北大学の構成員のみ利用可能です。

2 ご利用の前に
2 ご利用の前に

2.1 リモートログイン
2.1 リモートログイン

アプリケーションを利用するためには、フロントエンドサーバに SSH（Secure SHell）を用いてリ
アプリケーションを利用するためには、フロントエンドサーバに SSH（Secure SHell）を用いてリモートログインする必要があります。フロントエンドサーバにリモートログインする方法については、

モートログインする必要があります。フロントエンドサーバにリモートログインする方法については、
本誌 35 ページ「鍵作成とログイン方法」をご参照下さい。
本誌 35 ページ「鍵作成とログイン方法」をご参照下さい。
■利用申請からログインまで https://www.ss.cc.tohoku.ac.jp/ﬁrst-use/■利用申請からログインまで https://www.ss.cc.tohoku.ac.jp/ﬁrst-use/

2.2 ファイル転送
2.2 ファイル転送

ローカル PC で作成したインプットファイルをストレージシステムに転送したり、アプリケーション
ローカル PC で作成したインプットファイルをストレージシステムに転送したり、アプリケーションで実行した結果ファイルをローカル
PC に転送する場合は、ファイル転送が必要になります。ファイルで実行した結果ファイルをローカル
PC に転送する場合は、ファイル転送が必要になります。ファイル
転送については本誌 42 ページ「ストレージシステムの利用法」をご参照下さい。
転送については本誌 42 ページ「ストレージシステムの利用法」をご参照下さい。
■データ転送（ストレージ） https://www.ss.cc.tohoku.ac.jp/storage/■データ転送（ストレージ） https://www.ss.cc.tohoku.ac.jp/storage/

アプリケーションサービスの紹介― 49 ―

2.3 バッチジョブ形式での実行方法

バッチジョブ形式で実行するアプリケーションは、並列数の指定方法などが計算機やコンパイルに使用したコンパイラによって異なります。詳しくは各計算機の利用方マニュアルをご参照ください。
■ジョブの実行方法

https://www.ss.cc.tohoku.ac.jp/nqs/■AOBA-A の概要と利用方法■AOBA-B の概要と利用方法

https://www.ss.cc.tohoku.ac.jp/sx/
https://www.ss.cc.tohoku.ac.jp/lx/

2.4 GUI アプリケーションの利用

GUI アプリケーション（Mathematica および MATLAB）を利用する場合は、ローカル PC に X Window

System 環境の設定が必要です。
2.4.1 Linux からの利用

標準で X Window System がインストールされています。ローカル PC からリスト 1 のように-X オプション（または-Y オプション）を付けてログインしてください。

リスト 1 GUI アプリケーションを利用する場合のログイン方法（ 秘 密 鍵 の フ ァ イ ル 名 を i d _ r s a _ c c と し て ロ ー カ ル P C の ˜/. s s h 以 下 に 作 成 し た 場 合 ）【 ロ ー カ ル PCからログインサーバにSSH接続】
localhost$ ssh -i ˜/. ssh/ id_rsa_cc -X 利 用 者 番 号 @login .cc. tohoku .ac.jp【 ロ グ イ ン サ ー バ か ら に フ ロ ン ト エ ン ド サ ー バ に SSH接続】

login$ ssh -X front

X Forwarding によりローカル PC の画面にフロントエンドサーバで動作する GUI アプリケーション

画面が表示されます。

2.4.2 macOS からの利用

macOS では X Window System 環境の「XQuartz」をインストールして下さい。Linux からの利用と同様に利用可能です。接続方法はリスト 1 と同様です。

• XQuartz (https://www.xquartz.org/)
2.4.3 Windows からの利用■商用のアプリケーションを利用する場合のメーカから販売されています。

Windows 用 X サーバは、X サーバソフトとしていくつか• ASTEC-X (http://www.astec-x.com)• OpenText Exceed (https://www.macnica.net/opentext/exceed.html/)■無償のアプリケーションを利用する場合 VcXsrv Windows X Server が無料で利用できます。VcXsrvのインストールには WSL2 が必要になります。

• XcXsrv (https://sourceforge.net/projects/vcxsrv/)• WSL2 (https://docs.microsoft.com/ja-jp/windows/wsl/)

― 50 ―

SENAC Vol. 55, No. （2 2022. 4）

3 非経験的分子軌道計算プログラム Gaussian16

Gaussian は、Carnegie-Mellon 大学の Pople を中心として開発された分子軌道計算プログラムパッ

ケージです。広範囲にわたる非経験的モデルおよび半経験的モデルをサポートしています。
本センター AOBA-B での Gaussian の実行には、以下のような特長があります。

• 最大 128 並列までの並列処理が行え、実行時間の短縮が可能です。

• 256GB の搭載メモリを用いて、大規模なモデルの解析が可能です。

• スクラッチファイル (テンポラリファイル) を高速な SSD ディスクに置くことにより、ファイル
入出力時間が短縮されます。

■サービスホスト AOBA-B■バージョン Gaussian16 C.01

3.1 利用方法

以下は Gaussian 利用方法の概要です。

3.1.1 実行コマンド

Gaussian のインプットファイルは、拡張子を .com として作成します。 (リスト 2：e2-01.com )インプットファイルを Windows のエディタで作成した場合、拡張子.com のファイルは Windows では実行ファイルと認識されるため、 誤ってダブルクリックなどでインプットファイルを実行しないようご注意ください。また、ファイル転送ソフトでストレージシステムにインプットファイルを転送する際には、アスキーモードを指定して転送してください。

フロントエンドサーバにログイン後、 subg16 コマンドにキュー名 (lx) と利用ノード数 (1)、および入力プログラム名を指定することにより、AOBA-B のバッチリクエストとして実行されます。（リスト 3）
バッチリクエストについては、「ジョブの実行方法」をご参照下さい。

■ジョブの実行方法

https://www.ss.cc.tohoku.ac.jp/nqs/%NProc =128
#並 列 数%Mem =128 Gb#確 保 メ モ リ 量# RHF /6 -31G(d) Pop=Full Test

リスト 2 サンプルファイル e2-01.com

Formaldehyde Single Point
0 1
C
0.
O
0.
H .94
H -.94

0.
1.22
-.54
-.54

0.
0.
0.
0.

（最終行にも空行が必要）

リスト 3 e2-01.com を解析するジョブ投入コマンド[ front1 ˜]$ subg16 -q lx -b 1 e2 -01 ( s u b g 1 6 コ マ ン ド に 入 力 フ ァ イ ル を 指 定 す る 際 は 拡 張 子 . c o m を 省 き ま す )プ ロ ジ ェ ク ト コ ー ド ： un0000にリクエストを投入します
Request 12345. job submitted to queue: lx.

アプリケーションサービスの紹介― 51 ―

3.1.2 実行結果の確認

計算が終了すると、インプットファイル名に拡張子.log がつけられた結果ファイル (例： e2-01.log )が作成されます。計算結果をはじめ、CPU 時間などの計算機使用量に関する情報もここに含まれます。
正常終了ならば、このファイルの末尾に 「Normal termination of Gaussian 16.」というメッセージが出力されます。ファイルの末尾を表示する tail コマンドで確認できます。
（リスト 4）
リスト 4 実行結果の確認[ front1 ˜]$ tail e2 -01. log:
Job cpu time:
0 days 0 hours 0 minutes 10.9 seconds .
Elapsed time:
0 days 0 hours 0 minutes 0.7 seconds .
File lengths ( MBytes ): RWF=
72 Int=0 D2E=
0 Chk=
Normal termination of Gaussian 16 at Fri Apr 1 12:00:00 2021.

8 Scr=

8

結果ファイルの詳細な見方は、マニュアル等をご参照ください。

3.1.3 ユーティリティプログラム

formchk などのユーティリティプログラムは、/mnt/stfs/ap/g16 以下に配置されています。実行パスの設定を行うか、絶対パスでユーティリティプログラムを実行します。

3.2 マニュアル

本センター本館 1 階 利用相談室に以下の資料を備えてあります。

• 電子構造論による化学の探求 第 3 版, ガウシアン社, 2017

• Gaussian 16 Online Manual, https://gaussian.com/techsupport/• Gaussian プログラムによる量子化学計算マニュアル：堀憲次、丸善出版• すぐできる量子化学計算ビギナーズマニュアル：武次鉄也、講談社• すぐできる分子シミュレーションビギナーズマニュアル：長岡正隆、講談社• Gaussian プログラムで学ぶ情報化学・計算化学実験：堀憲次、丸善出版

3.3 利用方法についての問合せ

AOBA-B での利用方法についてのご質問は、サイバーサイエンスセンターの利用相談フォームをご利

用下さい。
■利用相談

https://www.ss.cc.tohoku.ac.jp/consultation/

3.4 Gaussian、 GaussView サイトライセンス

東北大学では Gaussian および GaussView のサイトライセンスを取得しており、東北大学内のコンピュータにインストールすることができます。

東北大学所属の教職員または学生の利用が可能です。学外で利用するコンピュータにインストール

することはできません。東北大学内でインストールしたコンピュータを学外に持ち出すこともできません。

3.4.1 申込み方法

インストールに必要な情報をお知らせしますので、利用を希望される方はサイトライセンス利用申請書作成フォーム（東北大学内からのみアクセス可能）

― 52 ―

SENAC Vol. 55, No. （2 2022. 4）https://www.ss.cc.tohoku.ac.jp/sscc/Gaussian-media/gaussian_license.html で PDF

ファイルを作成し、押印（またはサイン）の上、共同利用支援係までメールでご送付ください。
なお、申請には東北大メール（Gmail）アドレスが必要です。

4 反応経路自動探索プログラム GRRM17

GRRM は，NPO 法人 量子化学探索研究所 （https://iqce.jp/）で開発される反応経路自動探索プログ

ラムです。利用の方法，最新情報については https://iqce.jp/GRRM/ をご参照ください。

■サービスホスト

AOBA-B

4.1 利用方法

フロントエンドサーバにログイン後、ジョブスクリプトファイルと入力データファイル（.com）を作

成またはアップロードし、AOBA-B のバッチリクエストとして実行します。
バッチリクエストについては、「ジョブの実行方法」をご参照下さい。

■ジョブの実行方法

https://www.ss.cc.tohoku.ac.jp/nqs/

4.1.1 ジョブスクリプトファイルの作成

GRRM プログラムを実行するためのスクリプトファイルの作成例です。リスト 5 の例ではファイル

名を run.csh としてテキストファイルで保存します。

リスト 5 ジョブスクリプトファイル run.sh#!/ bin/sh#PBS -q lx -b 1
#PBS -l elapstim_req =1:00:00#PBS -N grrm -testsource /usr/ap/etc/ GRRM17 / config .sh
cd $PBS_O_WORKDIR
GRRMp test -p16 -h1#投 入 す る キ ュ ー 名 と ノ ー ド 数 を 指 定 （ 固 定 ）#実 行 す る 時 間 を 指 定#リ ク エ ス ト 名 を 指 定 。 指 定 し な い 場 合 は リ ク エ ス ト 番 号 が 指 定 さ れ る#GRRM17の環境変数ファイルを読み込む#リ ク エ ス ト を 投 入 し た デ ィ レ ク ト リ に 移 動#GRRMを起動# こ の 例 で は test. c o m フ ァ イ ル を 実 行 し 、 G R R M の J O B 並 列 数 が 1 6 、 計 算 時 間 の 上 限 値 が 1 時 間 に な る

4.1.2 入力データファイル (.com) の作成、またはファイルの転送

入力データファイルの拡張子は.com にします。リスト 6 では test.com ファイルを例に解説します。

test.com ファイルはフロントエンドサーバにログインしてテキストエディタで作成するか、ローカル PCの GaussView などで作成したものを、ストレージシステムにテキストモードでアップロードします。
# GRRM/MP2 /6 -31G

リスト 6 入力データファイル test.com

0 1
C -0.000000000000 -0.000000000000 -0.549482561269
O 0.000000000000 0.000000000000 0.708343639882
H 0.000000000000 0.934113144104 -1.131025039307
H -0.000000000000 -0.934113144104 -1.131025039307
Options
GauProc =2（最終行にも空行が必要）

この入力データファイルの場合、Gaussian の並列度は「2」になります。GRRM の並列 JOB 数との積が AOBA-B1 ノードの利用可能並列数 (128) を越えないように注意して下さい。test.com を「ジョブスクリプトファイルの作成」で作成した run.sh でジョブを投入した場合、並列数は以下のようになります。

(Gaussian 並列度 2) × (GRRM 並列 JOB 数 16) ＝ (並列数 32)

アプリケーションサービスの紹介― 53 ―

4.1.3 ジョブの投入、状況確認、取り消し・削除

front のコマンドプロンプト上で、以下のコマンドでリクエストを投入します。（リスト 7）
リスト 7 ジョブの投入[ front1 ˜]$ qsub run.sh
プ ロ ジ ェ ク ト コ ー ド ： un0000にリクエストを投入します
Request 12345. job submitted to queue: lx.

ジョブが投入され、待ち順に従ってリクエストが実行されます。リクエストを投入したディレクトリ

に test.com が保存されている必要があります。また、計算結果のファイルも同じディレクトリに作成されます。

投入したリクエストの状況確認と、取り消しは以下のコマンドで行います。
（リスト 8）[ front1 ˜]$ reqstat

リスト 8 リクエストの状況確認と、取り消し( リ ク エ ス ト の 状 況 確 認 、 リ ク エ ス ト 番 号 の 確 認)

[ front1 ˜]$ qdel 12345. job ( リ ク エ ス ト の 取 り 消 し ・ 削 除 )4.1.4 実行結果の確認、ファイルの転送

reqstat コマンドで実行したリクエストが表示されなければ実行が終了しています。実行結果のファイ

ルは リクエストを投入したディレクトリに作成されます。テキストエディタでファイルの内容を確認す
るか、ローカル PC にダウンロードして下さい。

4.1.5 実行についての注意事項

作業ディレクトリに「test.log」が存在する状態で、その GRRM JOB を開始する際に投入したものと同一のシェルスクリプトを再投入すると、再び GRRMp が起動され、GRRM JOB が再開されます。

GRRM のジョブを最初からやり直すときは、作業ディレクトリから識別子が.log や.rrm の全ファイルを消去してか ら (あるいは新しいディレクトリを作成して) 行う必要があります。

― 54 ―

SENAC Vol. 55, No. （2 2022. 4）5 統合型数値計算ソフトウェア Mathematica

Mathematica は Stephen Wolfram によって作られた、プログラミング言語を備えた数式処理システムです。 Mathematica の機能は、数値計算、記号計算、グラフィックスという 3 つに大別でき、この 3 つが一体となって使いやすいインタフェースを提供しています。
■サービスホスト

フロントエンドサーバ■バージョン version 13.0

5.1 利用方法

5.1.1 Mathematica の起動■GUI 版

GUI 版の Mathematica の起動には、フロントエンドサーバに接続する際に X forwarding の設定を行う必要があります。（2.4 項を参照）[ front1 ˜]$ mathematica

リスト 9 GUI 版の起動方法

GUI 版の Mathematica が起動します。（図 1）

図 1 GUI 版 Mathematica の起動画面■CUI 版

GUI を使用せず、コマンドライン上で起動することもできます。（リスト 10）
リスト 10

テキスト版の起動方法[ front1 ˜]$ math
Mathematica 12.2.0 Kernel for Linux x86 (64- bit)
Copyright 1988 -2021 Wolfram Research , Inc.
In [1]:=

5.2 マニュアル・参考資料

Mathematica の使い方は、 マニュアル・参考資料 や、Web などをご参照ください。
• Wolfram 社 https://www.wolfram.com/

本センター本館１階 利用相談室に、以下の資料を備えてあります。

アプリケーションサービスの紹介― 55 ―• スティーブンウルフラム Mathematica ブック (日本語版)：トッパン• Mathematica 方法と応用：J.W. グレイ、サイエンティスト社• Mathematica プログラミング技法：R. メーダー、トッパン• 入門 Mathematica：日本 Mathematica ユーザー会、東京電機大学出版局• はやわかり Mathematica：榊原進、共立出版• もっと Mathematica で数学を：吉田孝之、培風館

6 数値解析ソフトウェア MATLAB

MATLAB は高機能な数値計算機能と多彩な可視化機能を備えた技術計算ソフトウェアです。科学的、工学的分野の様々な数値計算 (特に行列演算)、データ解析、シミュレーション、およびビジュアラ
イゼーションのための統合環境を提供しています。提供される全ての Toolbox が利用可能です。

■サービスホスト

フロントエンドサーバ、AOBA-B■バージョン R2022a

6.1 利用方法

6.1.1 MATLAB の起動■GUI 版

GUI 版 MATLAB の起動には、フロントエンドサーバに接続する際に X forwarding の設定を行う必要があります。（2.4 項を参照）[ front1 ˜]$ matlab

リスト 11 GUI 版 MATLAB の起動方法

GUI 版の MATLAB が起動します。（図 2）

図 2 GUI 版 MATLAB の起動画面

― 56 ―■CUI 版

SENAC Vol. 55, No. （2 2022. 4）

GUI を使用せず、コマンドライン上で起動することもできます。（リスト 12）
リスト 12 CUI 版 MATLAB の起動方法[ front1 ˜]$ matlab -nojvm -nosplash -nodesktop -nodisplay< M A T L A B (R) >Copyright 1984 -2022 The MathWorks , Inc.
R2022a (9.12.0.1884302) 64- bit ( glnxa64 )February 16, 2022
For online documentation , see https: // www. mathworks .com/ support
For product information , visit www. mathworks .com.
>>

6.2 バッチジョブ実行

MATLAB の並列処理機能を使用し、AOBA-B の 1 ノードを利用して 128 並列までの処理が可能です。最大メモリも 256GB まで利用可能です。大規模な計算にご利用ください。

バッチジョブ実行ではグラフ描画など画面出力のあるプログラムや、対話的な処理は行えません。
バッチリクエストについては、「ジョブの実行方法」をご参照下さい。

■ジョブの実行方法 https://www.ss.cc.tohoku.ac.jp/nqs/

function として作成した test（リスト 13）を実行するためにはジョブスクリプトファイル（リスト 14）を作成します。

リスト 13 サンプルファイル test.mfunction test

p= parpool (’local ’ ,128) % 並 列 プ ー ル を 作 成 （ 128 個 ）n = 1024;
A = 500;
a = zeros (1,n);parfor i = 1:na(i) = max(abs(eig(rand(A ))));
end
poolobj = gcp(’nocreate ’);delete ( poolobj ) % 並 列 プ ー ル を 削 除#!/ bin/sh#PBS -q lx -b 1
#PBS -l elapstim_req =1:00:00
cd $PBS_O_WORKDIR
matlab -batch test

リスト 14

バッチリクエストファイル run.sh#投 入 す る キ ュ ー 名 と ノ ー ド 数 を 指 定 （ 固 定 ）#実 行 す る 時 間 を 指 定#リ ク エ ス ト を 投 入 し た デ ィ レ ク ト リ に 移 動#バ ッ チ ジ ョ ブ 形 式 で 実 行

以下のコマンドでリクエストを投入します。

リスト 15 バッチリクエストの投入方法[ front1 ˜]$ qsub run.sh
プ ロ ジ ェ ク ト コ ー ド ： un0000にリクエストを投入します
Request 12345. job submitted to queue: lx.

6.3 Parallel Server による複数ノード実行

Parallel Server の機能により、AOBA-B の複数ノードを利用した実行が可能です。利用方法については以下のマニュアルをご参照下さい。

■AOBA-B での MATLAB Parallel Server の使い方

https://www.ss.cc.tohoku.ac.jp/pdf/AOBA-B-Parallel-Server.pdf

アプリケーションサービスの紹介― 57 ―

6.4 サンプルプログラム

MATLAB には豊富なデモがありますので、ご利用ください。MATLAB 上で、demo コマンドを実行

すると、デモ画面が開きます。

6.5 マニュアル・参考資料

MATLAB の使い方は、 マニュアル・参考資料などをご参照ください。

■マニュアル日本語オンラインマニュアルが公開されています。以下のページをご参照ください。

• https://www.mathworks.co.jp/help/ja JP/techdoc/index.html■参考資料

本センター本館１階 利用相談室に、以下の資料を備えてあります。

• MATLAB による制御理論の基礎：野波健蔵、東京電機大学出版局• MATLAB による制御のためのシステム同定：足立修一、東京電機大学出版局• だれでもわかる MATLAB：池原雅章、培風館• はやわかり MATLAB 第 2 版：芦野隆一、共立出版• 最新 MATLAB ハンドブック第 3 版：小林一行、秀和システム• MATLAB グラフィックス集：小国 力、朝倉書店• MATLAB と利用の実際：小国 力、サイエンス社• MATLAB の総合応用：高谷邦夫、森北出版• 最新使える！ MATLAB：青山貴伸、講談社• 使える！ MATLAB/Simulink プログラミング：青山貴伸、講談社• MATLAB による画像＆映像信号処理：村松正吾、CQ 出版• Matlab によるグラフ描画：西村竜一(広報誌 SENAC Vol.37 No.1 (2004-1))• 高機能数値計算・可視化機能ソフト MATLAB の基本的な使い方：陳国曜 他
Vol.46 No.3 (2013-7))

6.6 利用方法についての問合せ(広報誌 SENAC

利用方法についてのご質問は、MathWorks サポート窓口、コミュニティ Q&A サイトをご利用下さい。

■MathWorks サポート窓口 https://jp.mathworks.com/support/contact us.html■コミュニティ Q&A サイト https://jp.mathworks.com/matlabcentral/answers/

― 58 ―

SENAC Vol. 55, No. （2 2022. 4）

7 熱流体解析ソフトウェア OpenFOAM

OpenFOAM（Open source Field Operation And Manipulation の略称) は数値解析開発、及び数値流体力学を含む連続体力学の前後処理用の C++ 製ツールボックスです。

使い方についての問合せは開発元 Web サイトやユーザコミュニティをご利用下さい。

7.1 OpenCFD 版■開発元 Web サイト https://www.openfoam.com■サービスホスト

AOBA-B■使用コンパイラ

Intel OneAPI■バージョン v2112

OpenFOAM の実行文を記述したジョブスクリプトファイル（リスト 16）を作成して qsub コマンドでリクエストを投入します。

リスト 16 ジョブスクリプトファイル例 job.sh#!/ bin/sh#PBS -T intmpi#PBS -l elapstim_req =24:00:00#PBS -q lx -b 1
source /usr/ap/etc/ OpenFOAM / v2112 / config .sh
cd $PBS_O_WORKDIR・ ・ ・ OpenFOAM実行文・・・

以下のコマンドでリクエストを投入します。

リスト 17 バッチリクエストの投入方法[ front1 ˜]$ qsub job.sh
プ ロ ジ ェ ク ト コ ー ド ： un0000にリクエストを投入します
Request 12345. job submitted to queue: lx.

7.2 The OpenFOAM Foundation 版■開発元 Web サイト https://openfoam.org■サービスホスト

AOBA-B■使用コンパイラ

Intel OneAPI■バージョン 9

OpenFOAM の実行文を記述したジョブスクリプトファイル（リスト 18）を作成して qsub コマンドでリクエストを投入します。

リスト 18 ジョブスクリプトファイル例 job.sh#!/ bin/sh#PBS -T intmpi#PBS -l elapstim_req =24:00:00#PBS -q lx -b 1
source /usr/ap/etc/ OpenFOAM /9/ config .sh
cd $PBS_O_WORKDIR・ ・ ・ OpenFOAM実行文・・・

以下のコマンドでリクエストを投入します。

アプリケーションサービスの紹介― 59 ―

リスト 19 バッチリクエストの投入方法[ front1 ˜]$ qsub job.sh
プ ロ ジ ェ ク ト コ ー ド ： un0000にリクエストを投入します
Request 12345. job submitted to queue: lx.

バッチリクエストについては、「ジョブの実行方法」をご参照下さい。

■ジョブの実行方法

https://www.ss.cc.tohoku.ac.jp/nqs/

8 第一原理計算統合パッケージ Quantum Espresso

Quantum ESPRESSO は、第一原理電子構造計算と材料モデリングのためのオープンソース統合パッケージです。

使い方についての問合せは開発元 Web サイトやユーザコミュニティをご利用下さい。

■開発元 Web サイト https://www.quantum-espresso.org■サービスホスト

AOBA-A および AOBA-B

8.1 AOBA-A での実行

■バージョン 6.3（pw.x のみの提供）

Quantum Espresso の実行文を記述したジョブスクリプトファイル（リスト 20）を作成して qsub コマ

ンドでリクエストを投入します。

リスト 20 ジョブスクリプトファイル例 job.sh#!/ bin/sh#PBS -l elapstim_req =24:00:00#PBS -q sx --venode =1
source /usr/ap/etc/QE/SX /6.3/ config .sh
cd $PBS_O_WORKDIR・ ・ ・ Quantum E s p r e s s o 実 行 文 ・ ・ ・

以下のコマンドでリクエストを投入します。

リスト 21 バッチリクエストの投入方法[ front1 ˜]$ qsub job.sh
プ ロ ジ ェ ク ト コ ー ド ： un0000にリクエストを投入します
Request 12345. job submitted to queue: sx.

8.2 AOBA-B での実行
■バージョン 7.0■使用コンパイラ

Intel OneAPI

Quantum Espresso の実行文を記述したジョブスクリプトファイル（リスト 22）を作成して qsub コマ

ンドでリクエストを投入します。

リスト 22 ジョブスクリプトファイル例 job.sh#!/ bin/sh#PBS -T intmpi#PBS -l elapstim_req =24:00:00#PBS -q lx -b 1
source /usr/ap/etc/QE/LX /7.0/ config .sh
cd $PBS_O_WORKDIR・ ・ ・ Quantum E s p r e s s o 実 行 文 ・ ・ ・

以下のコマンドでリクエストを投入します。

― 60 ―

SENAC Vol. 55, No. （2 2022. 4）

リスト 23 バッチリクエストの投入方法[ front1 ˜]$ qsub job.sh
プ ロ ジ ェ ク ト コ ー ド ： un0000にリクエストを投入します
Request 12345. job submitted to queue: lx.

バッチリクエストについては、「ジョブの実行方法」をご参照下さい。

■ジョブの実行方法

https://www.ss.cc.tohoku.ac.jp/nqs/

9 分子動力学シミュレーション GROMACSGROMACS は、分子動力学 (MD) 計算プログラムです。MPI あるいはスレッドでの効率のよい並列計算を用いて高速処理が可能です。

使い方についての問合せは開発元 Web サイトやユーザコミュニティをご利用下さい。

■開発元 Web サイト https://www.gromacs.org/■サービスホスト

AOBA-B■使用コンパイラ

Intel OneAPI■バージョン 2021.2

GROMACS の実行文を記述したジョブスクリプトファイル（リスト 24）を作成して qsub コマンドでリクエストを投入します。

リスト 24 ジョブスクリプトファイル例 job.sh#!/ bin/sh#PBS -T intmpi#PBS -l elapstim_req =24:00:00#PBS -q lx -b 1
source /usr/ap/etc/ GROMACS /2021.2/ config .sh
cd $PBS_O_WORKDIR・ ・ ・ GROMACS実行文・・・

以下のコマンドでリクエストを投入します。

リスト 25 バッチリクエストの投入方法[ front1 ˜]$ qsub job.sh
プ ロ ジ ェ ク ト コ ー ド ： un0000にリクエストを投入します
Request 12345. job submitted to queue: lx.

バッチリクエストについては、「ジョブの実行方法」をご参照下さい。

■ジョブの実行方法

https://www.ss.cc.tohoku.ac.jp/nqs/

10 古典分子動力学計算プログラム LAMMPSLAMMPS は古典分子動力学 (MD) 計算プログラムです。金属や半導体などの固体、生体分子やポリマーなどのソフトマター、メゾスコピック系物質などに対応した多彩なポテンシャルが用意されています。

使い方についての問合せは開発元 Web サイトやユーザコミュニティをご利用下さい。

■開発元 Web サイト https://www.lammps.org/■サービスホスト

AOBA-B■バージョン 29Sep2021

アプリケーションサービスの紹介■使用コンパイラ― 61 ―

Intel OneAPILAMMPS の実行文を記述したジョブスクリプトファイル（リスト 26）を作成して qsub コマンドでリクエストを投入します。

リスト 26 ジョブスクリプトファイル例 job.sh#!/ bin/sh#PBS -T intmpi#PBS -l elapstim_req =24:00:00#PBS -q lx -b 1
source /usr/ap/etc/ LAMMPS /29 Sep2021 / config .sh
cd $PBS_O_WORKDIR・ ・ ・ LAMMPS実行文・・・

以下のコマンドでリクエストを投入します。

リスト 27 バッチリクエストの投入方法[ front1 ˜]$ qsub job.sh
プ ロ ジ ェ ク ト コ ー ド ： un0000にリクエストを投入します
Request 12345. job submitted to queue: lx.

バッチリクエストについては、「ジョブの実行方法」をご参照下さい。

■ジョブの実行方法

https://www.ss.cc.tohoku.ac.jp/nqs/

11 コンテナプラットフォーム Singularity

Singularity は HPC (High Performance Computing) 環境向けに設計されたオープンソースのコンテナ型プラットフォームです。様々なアプリケーションのイメージファイルがインターネット上で公開されており、インストールが複雑なアプリケーションでも容易に実行環境を構築することができます。
■サービスホスト

AOBA-B

Singularity の利用方法については、以下のマニュアルをご参照下さい。

■Singularity の利用方法

https://www.ss.cc.tohoku.ac.jp/pdf/AOBA-B-Singularity.pdf

12 プログラム対話型実行環境 Jupyter Notebook

Jupyter Notebook（ジュピター・ノートブック）は、ブラウザ上で動作するプログラムの対話型実行環境です。

ノートブックと呼ばれるファイルに python などのプログラムを記述し、実行結果を逐次確認しながら、データ分析を進めることができます。Jupyter Notebook 上では、様々なライブラリを使うことができ、TensorFlow（テンソルフロー）などの機械学習やディープラーニング用のライブラリも動作させることができます。
■サービスホスト

フロントエンドサーバ

Jupyter Notebook の利用方法については、以下のマニュアルをご参照下さい。

■Jupyter Notebook の利用方法

https://www.ss.cc.tohoku.ac.jp/pdf/AOBA-jupyterNotebook.pdf

― 62 ―

SENAC Vol. 55, No. （2 2022. 4）

13 RIST が整備するアプリケーションソフトウェア

RIST では「富岳」を中核とした HPCI において利用が多い OSS（オープンソースソフトウェア）

や、ポスト「京」重点課題アプリケーション開発に関するプロジェクトなどより推薦された、利用が多いまたは今後期待される日本で開発されたアプリケーションソフトウェアを対象として、プリインス
トールや利用者に有益な情報を提供し、容易に利用いただけるように整備しています。

サイバーサイエンスセンターではその一部のアプリケーションソフトウェアについて利用環境を提供しています。利用方法については、以下のリンク先の情報をご参照ください。
■RIST が整備するアプリケーションソフトウェア■サービスホスト

AOBA-A■サービスホスト

AOBA-BHΦ

材料バイオ分野 ABINIT-MP
材料バイオ分野 ABINIT-MPHΦ

OpenMXPHASE/0

SMASH

CAE 分野 FrontFlow/bluehttps://www.hpci-oﬃce.jp/pages/appli software

FrontISTR

GENESISMODYLAS

NTChemSALMONSENAC Vol. 55, No. 2（2022. 4）― 63 ―[報 告]

東北大学と NEC、マテリアルズインテグレーションシステムによる航空機用複合材料開発シミュレーションのサービスを提供開始～スーパーコンピュータを活用し、高機能材料開発に貢献～
サイバーサイエンスセンターのスーパーコンピュータ AOBA を活用したシミュレーションサー
ビスの提供が開始されました。
詳細は、東北大学 HP よりご確認ください。
【東北大学 HP】
https://www.tohoku.ac.jp/japanese/2022/03/press20220318-02-service.html

― 64 ―

SENAC Vol. 55, No. （2 2022. 4）[スーパーコンピュータ AOBA のお知らせより]
東北大学サイバーサイエンスセンター大規模科学計算システムウェブサイトに掲載されたお知らせの一部を転載しています。

https://www.ss.cc.tohoku.ac.jp/information/

令和 4 年度利用負担金について
令和 4 年度の利用負担金についてお知らせします。詳細は以下をご覧ください。
別表１ 基本利用負担金【大学･学術利用】
区 分
項 目
利用
負担額及び課金対象時間
形態
スーパー
共有
利用 VE 数 1(実行数、実行時間の制限有)
無料
演 算コンピュータ(無料)
負担経費
課金対象時間 ＝共有（利用 VE 数÷8 を切上げた数）× 経過時間(秒)(従量)
課金対象時間 1 時間につき
125 円

並列
コンピュータ

共有
(定額)
占有
共有

負担額 10 万円につき課金対象時間 800 時間分使用可能
利用 VE 数 8
利用期間 3 ヶ月につき270,000 円
課金対象時間 ＝ 利用ノード数× 経過時間(秒)
課金対象時間 1 時間につき22 円
負担額 10 万円につき課金対象時間 4,600 時間分使用可能
利用ノード数 1
利用期間 3 ヶ月につき
5TB まで無料、追加容量 1TB につき年額

47,000 円
3,000 円

占有
10TB まで無料、追加容量 1TB につき年額
大判プリンタによるカラープリント
フォト光沢用紙 1 枚につきクロス 1 枚につき3,000 円
600 円
1,200 円

ファイル負担経費
出力
負担経費

共有
(定額)
占有
共有
(従量)

備考
１ 負担額が無料となるのは専用のキューで実行されたものとし、制限時間を超えた場合は強制終了する。
２ 演算負担経費の課金対象時間については半期毎(4 月から 9 月及び 10 月から 3 月)に合計し、1 時間
未満を切上げて負担金を請求する。
３ 演算負担経費について定額制を選択した場合はスーパーコンピュータ及び並列コンピュータを課
金対象時間の範囲内で共用できる。
４ 占有利用期間は年度を超えないものとし、期間中に障害、メンテナンス作業が発生した場合においても、原則利用期間の延長はしない。
５ ファイル負担経費については申請日から当該年度末までの料金とする。運用期間が 1 年に満たない
場合は、月割りをもって計算した額とする。

スーパーコンピュータ AOBA のお知らせより― 65 ―

別表２ 基本利用負担金【民間機関利用（成果公開型）】
区 分
項 目
利用
負担額及び課金対象時間
形態
スーパー
共有
利用 VE 数 1(実行数、実行時間の制限有)
無料
演 算コンピュータ(無料)
負担経費
課金対象時間 ＝共有（利用 VE 数÷8 を切上げた数）× 経過時間(秒)(従量)
課金対象時間 1 時間につき
250 円

並列
コンピュータ

共有
(定額)
占有
共有

負担額 20 万円につき課金対象時間 800 時間分使用可能
利用 VE 数 8
利用期間 3 ヶ月につき
540,000 円
課金対象時間 ＝ 利用ノード数× 経過時間(秒)
課金対象時間 1 時間につき44 円

負担額 20 万円につき課金対象時間 4,600 時間分使用可能
利用ノード数 1
利用期間 3 ヶ月につき
5TB まで無料、追加容量 1TB につき年額

94,000 円
6,000 円

占有
10TB まで無料、追加容量 1TB につき年額
大判プリンタによるカラープリント
フォト光沢用紙 1 枚につき6,000 円
1,200 円

ファイル負担経費
出力
負担経費

共有
(定額)
占有
共有
(従量)

クロス 1 枚につき2,400 円

備考
１ 負担額が無料となるのは専用のキューで実行されたものとし、制限時間を超えた場合は強制終了する。
２ 演算負担経費の課金対象時間については半期毎(4 月から 9 月及び 10 月から 3 月)に合計し、1 時間
未満を切上げて負担金を請求する。
３ 演算負担経費について定額制を選択した場合はスーパーコンピュータ及び並列コンピュータを課
金対象時間の範囲内で共用できる。
４ 占有利用期間は年度を超えないものとし、期間中に障害、メンテナンス作業が発生した場合においても、原則利用期間の延長はしない。
５ ファイル負担経費については申請日から当該年度末までの料金とする。運用期間が 1 年に満たない
場合は、月割りをもって計算した額とする。

― 66 ―

SENAC Vol. 55, No. （2 2022. 4）

別表３ 基本利用負担金【民間機関利用（成果非公開型）
】
区 分
項 目
利用
負担額及び課金対象時間
形態
スーパー
共有
利用 VE 数 1(実行数、実行時間の制限有)
無料
演 算コンピュータ(無料)
負担経費
課金対象時間 ＝（利用 VE 数÷8 を切上げた数）× 経過時間(秒)共有(従量)
課金対象時間 1 時間につき
500 円

並列
コンピュータ

ファイル負担経費
出力
負担経費

共有
(定額)
占有
共有
(従量)
共有
(定額)
占有
共有

負担額 40 万円につき課金対象時間 800 時間分使用可能
利用 VE 数 8
利用期間 3 ヶ月につき1,080,000 円
課金対象時間 ＝ 利用ノード数× 経過時間(秒)
課金対象時間 1 時間につき88 円
負担額 40 万円につき課金対象時間 4,600 時間分使用可能
利用ノード数 1
利用期間 3 ヶ月につき
5TB まで無料、追加容量 1TB につき年額

占有10TB まで無料、追加容量 1TB につき年額
大判プリンタによるカラープリント
フォト光沢用紙 1 枚につきクロス 1 枚につき188,000 円
12,000 円
12,000 円
2,400 円
4,800 円

備考
１ 負担額が無料となるのは専用のキューで実行されたものとし、制限時間を超えた場合は強制終了する。
２ 演算負担経費の課金対象時間については半期毎(4 月から 9 月及び 10 月から 3 月)に合計し、1 時間
未満を切上げて負担金を請求する。
３ 演算負担経費について定額制を選択した場合はスーパーコンピュータ及び並列コンピュータを課
金対象時間の範囲内で共用できる。
４ 占有利用期間は年度を超えないものとし、期間中に障害、メンテナンス作業が発生した場合においても、原則利用期間の延長はしない。
５ ファイル負担経費については申請日から当該年度末までの料金とする。運用期間が 1 年に満たない
場合は、月割りをもって計算した額とする。

（共同利用支援係）

スーパーコンピュータ AOBA のお知らせより― 67 ―

令和 4 年度共同研究について

本センターでは、大規模科学計算システムの利用者と共同でプログラムやアルゴリズムを開発する
共同研究を行っています。令和 4 年度の募集に応募されたものについて共同研究専門部会で審査の結果、14 件が採択されましたのでお知らせします。
[A] 若手・女性研究者支援課題
No. 申請者
所属
A-1 松川 嘉也 東北大学
大学院工学研究科
A-2 PHAM Ngoc
大阪大学
大学院工学研究科
Thanh
A-3

WANG
YUELIN[B] 萌芽型課題
No. 申請者B-1 有馬 卓司
B-2

春日 貴志
林 優一

B-3
B-4

佐々木大輔
高橋 俊
松岡 浩

B-5

山本

義暢

[C] 一般課題
No. 申請者
C-1 伊藤 純至
C-2

越村

俊一

C-3

茂田

正哉

C-4

塚原

隆裕

C-5

藤井

孝藏

C-6

藤崎

貴也

研究課題
熱分解反応場における温度の変動が化学反応速度に及ぼす影響の解明
First-principles multiscale simulation ofsintering process of perovskite-supported
metal nanoclusters
貴金属フリー炭素系材料における酸素還元反応の大規模第一原理電子状態計算による研究

大阪大学
大学院工学研究科
所属
東京農工大学
大学院工学研究院
長野工業高等専門学校
奈良先端科学技術大学院大学
金沢工業大学 工学部
東海大学 工学部
技術士事務所 AI コンピュー
ティングラボ
山梨大学
大学院総合研究部

研究課題
大規模電磁界解析を可能とする複数領域 FDTD 法の開発に関する研究
通信用コネクタ構造を組み込んだプリント基板
解析の高速化所属
東北大学
大学院理学研究科
東北大学
災害科学国際研究所
東北大学
大学院工学研究科

研究課題
アンサンブルカルマンフィルタと非静力学数値モデルを用いた日本領域の気象の長期再解析
SX-Aurora TSUBASA の津波のリアルタイム予測計
算の性能評価高エンタルピープラズマ－非電離気体共存系乱流場における物質変換プロセスの大規模シミュ
レーション
表面粗さ・主流乱れが影響した後退平板境界層における乱流遷移の直接数値解析
複雑流れ現象の理解とその効果的制御法に関する研究− １クラスター展開モンテカルロ法による有限温度でのイオン伝導体の原子配置の決定

東京理科大学
理工学部
東京理科大学
工学部
東北大学
多元物質科学研究所

直交格子法による移動境界問題の解法に関する研究
リカレント型ビット演算による縦渦挙動のマル
チスケール創発解析
HPC 及びデータサイエンス融合による高レイノル
ズ数乱流解析の新展開（スーパーコンピューティング研究部、共同研究支援係）

― 68 ―

SENAC Vol. 55, No. （2 2022. 4）

商用アプリケーションのバージョンアップについて

数式処理システム「Mathematica」および、数値解析ソフトウェア「MATLAB」のバージョンアッ
プを行いましたのでお知らせいたします。
新機能の概要、機能の詳細については開発元 Web サイトをご参照ください。
Mat hemati ca
・バージョン：13.0
・バージョンアップ日：2022 年 4 月 11 日・サービスホスト：フロントエンドサーバ
・起動コマンド：（GUI 版）mathematica （テキスト版）math・開発元 Web サイト：https://www.wolfram.com/mathematica/new-in-13/
MAT LAB・バージョン：2022a・バージョンアップ日：2022 年 4 月 11 日・サービスホスト：フロントエンドサーバ、AOBA-B
・起動コマンド：（GUI 版）matlab （テキスト版）matlab –nojvm –nosplash –nodesktop -nodisplay
・開発元 Web サイト：
https://jp.mathworks.com/products/new_products/latest_features.html?s_tid=hp_release_2022a(共同利用支援係）

スーパーコンピュータ AOBA のお知らせより― 69 ―

計算科学・計算機科学人材育成のためのスーパーコンピュータ無償提供制度について

東北大学サイバーサイエンスセンターでは、計算科学・計算機科学分野での教育貢献・人材育成を目的として、無料で大規模科学計算システムを利用できる制度を用意しております。提供の対象
は、大学院・学部での講義実習等の教育目的(卒業論文、修士論文、博士論文での利用を除く)に限ります。利用を希望される場合は以下の情報を添えて、講義開始の 2 週間前までに edu-prog@cc.tohoku.ac.jp
宛お申し込みください。
・講義担当者氏名・同所属・同連絡先（住所，電話，電子メール）・講義名・講義実施日時（１セメスターの中で実習を予定している回数）
・センター端末機室等での実習利用希望の有無（必要であれば予定日時）
・講師派遣の希望の有無
・講義シラバス
・講義ウェブ（もし用意されていれば）・受講者数（予定）・必要とする理由（利用目的：例えば、数値シミュレーションの研修を行うなど）・期待できる教育効果
・居住性チェックリストの提出（受講者に外国人が居る場合）
参照：https://www.ss.cc.tohoku.ac.jp/apply-for-use/#toc3
・その他（センターへの要望等）
なお、講義終了後、報告書（広報誌 SENAC へ掲載）の提出をお願いいたします。
たくさんのお申し込みをお待ちしております。不明な点は、edu-prog@cc.tohoku.ac.jp までお問い合
わせください。
(スーパーコンピューティング研究部，共同利用支援係)

― 70 ―

SENAC Vol. 55, No. （2 2022. 4）

民間企業利用サービスについて
東北大学サイバーサイエンスセンターでは、社会貢献の一環として大学で開発された応用ソフト
ウェアとスーパーコンピュータを、民間企業の方が無償または有償にてご利用頂ける制度を用意しております。本サービスにおける利用課題区分は以下の２つとなります。
・大規模計算利用(有償利用)
・トライアルユース(無償利用)
詳細については以下を参照し、利用を希望される場合は共同利用支援係までお申し込みください。
https://www.ss.cc.tohoku.ac.jp/business/【問い合わせ先】
共同利用支援係（uketuke@cc.tohoku.ac.jp）(共同利用支援係)

大規模科学計算システムの機関（部局）単位での利用について

東北大学サイバーサイエンスセンターでは、大規模科学計算システムをご利用いただくにあたり、利用負担金を利用者単位のほか、機関（部局）単位で年間定額をお支払いいただくことで利用できるサービスも提供しております。このサービスは、機関（部局）単位でお申し込みいただくことにより、その構成員であれば、各研究室が個別に利用負担金を支払うことなく、下記システムを利用できる仕組みとなっております。
これまで計算機を利用する機会がなかった研究者による新たなニーズへの対応や研究室の計算機では実行できなかった大規模シミュレーションが実行可能であり、また自前で計算機を導入するためのコストや運用コストも削減可能です。すでにご利用いただいている機関（部局）からは、当初の予想を上回るご利用をいただき、ご好評をいただいております。
占有利用・共有利用については必要に応じて取り混ぜながら、ご予算に合わせて、年間定額により利用することが可能となっておりますので、ぜひご相談ください。
記【利用可能なシステム】
・サブシステム AOBA-A・サブシステム AOBA-B
・ストレージシステム
・大判カラープリンター（光沢紙、ソフトクロス紙）【問い合わせ先】
共同利用支援係(uketuke@cc.tohoku.ac.jp)（共同研究支援係，共同利用支援係）

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

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

― 72 ―

SENAC Vol. 55, No. （2 2022. 4）スタッフ便りスタッフ便りを書く順番が回ってきました。←四年前と二年前の書き出しがこれでした。四年前はその日の出来事を日記風に記述するという手法でスペースを埋め、二年前は大学での業務スタイルがオンラインに移行して一変したという誰もが書きそうなことを書いていますね。それより前の 2016
年は例年 8 月に講習会を担当していることを、2013 年はオープンキャンパスでカード組を用いた秘密計算のデモをした様子を書いていたようです。カード組を用いた秘密計算は私の研究テーマのひとつであり、当時は細々と（？）研究を進めていましたが、その後、この「カードベース暗号」の研究
分野は急速に発展し、今ではサイバーサイエンスセンターの研究力（トップ 10％論文、FWCI）の中心となるまでに成長しました。今後も、この東北大学発の特徴ある研究分野の発展に貢献したいと思い
ます。
（T.M）
みなさんはラジオを聴く機会はありますでしょうか？防災グッズの中に災害時や、停電時用に携帯ラジオが入ってるけど、普段は車の中でも聴かなくなった、という方も多いかと思います。
私は小学生の頃からラジオを聴くのが大好きで、中学生の頃は受験勉強をしながら深夜ラジオを、
高校生の頃はカセットテープに録音してウォークマンで聴きながら電車通学を、大学生になって来仙してからは関東のラジオ番組を雑音混じりでもなんとか聴こうと試行錯誤、そんな思い出がはっきりと記憶に残っています。
今では「radiko（ラジコ）」で全国各地のラジオ番組をインターネットで簡単に聴ける、夢のような時代になりました。北海道も広島も沖縄のラジオ番組も、雑音無しで楽しめてしまうんです。数えてみたら全国各地のラジオ番組を毎週 50 時間分も聴いていました！「〜〜しながら」聴けるラジオだ
からこその楽しみ方ですね。（Y.T）

【サイバーサイエンスセンター・情報部情報基盤課スタッフ退職者、異動、新規採用のお知らせ】
2022.3.31 付け
[退職]
八巻 俊輔
助教（辞職・東北工業大学准教授へ）
[転出]
伊藤 昭彦
総務係長（工学研究科総務課専門職員へ）
佐々木徳仁
会計係主任（病院医事課収入管理係主任へ）
2022.4.1 付け
[転入]
杉田 典大
教授（工学研究科准教授から）ギリエ ハルハ ルイス アルベルト
助教（電気通信研究所特任助教から）
今野 義則
総務係長（宮城教育大学学生課専門職員から）[新規採用]
小林 一斗
会計係

SENAC 編集部会
滝沢寛之 水木敬明 後藤英昭
今野義則 早坂和勝 大泉健治
斉藤くみ子

編集・発行

PDF 作成高橋慧智
小野 敏

令和 4 年 4 月発行
東北大学
サイバーサイエンスセンター
仙台市青葉区荒巻字青葉 6-3
郵便番号 980-8578
株式会社 東誠社

ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢱ AOBA ࢩࢫࢸ୍࣒ぴィ⟬ᶵࢩࢫࢸ࣒

ᶵ✀ࢧࣈࢩࢫࢸ࣒ AOBA-A

SX-Aurora TSUBASAࢧࣈࢩࢫࢸ࣒ AOBA-B

LX 406Rz-2

ࢧ࣮ࣂ࡜࣍ࢫࢺྡ
ࣟࢢ࢖ࣥࢧ࣮ࣂ

login.cc.tohoku.ac.jp

ࢹ࣮ࢱ㌿㏦ࢧ࣮ࣂ

file.cc.tohoku.ac.jpࢧ࣮ࣅࢫ᫬㛫
฼⏝ࢩࢫࢸ࣒ྡ➼฼⏝᫬㛫ᖏࢧࣈࢩࢫࢸ࣒ AOBA-A

㐃

⥆

㐠

㌿

ࢧࣈࢩࢫࢸ࣒ AOBA-B

㐃

⥆

㐠

㌿

ࢧ࣮ࣂ

㐃

⥆

㐠

㌿

㤋ෆ฼⏝

ᖹ᪥ 8:30㹼21:00

ࢧࣈࢩࢫࢸ࣒ AOBA-A ࡢ฼⏝ᙧែ࡜ไ㝈್
฼⏝ᙧែ

࣮࢟ࣗྡ

↓ᩱ

sxf

VE ᩘͤ
1

1VE

1

1VE

sx
2ࠥ256

ඹ᭷
sxmix

ᐇ⾜ᙧែ

2ࠥ8

༨᭷

᭱኱⤒㐣᫬㛫
᪤ᐃ್/᭱኱್

࣓ࣔࣜࢧ࢖ࢬ

1 ᫬㛫/1 ᫬㛫

8VE ༢఩࡛☜ಖ㸦VH ࢆඹ⏝ࡋ࡞࠸㸧 72 ᫬㛫/720 ᫬㛫

48GB×VE ᩘ

1VE ༢఩࡛☜ಖ㸦VH ࢆඹ⏝ࡍࡿ㸧ಶูタᐃ

ͤ 2VE௨ୖࢆ฼⏝ࡋࡓ୪ิᐇ⾜࡟ࡣMPIࡢ฼⏝ࡀᚲ⏝ࢧࣈࢩࢫࢸ࣒ AOBA-B ࡢ฼⏝ᙧែ࡜ไ㝈್
฼⏝ᙧែ

࣮࢟ࣗྡ

ඹ᭷

lx

ࣀ࣮ࢻᩘͤ
1ࠥ16

᭱኱⤒㐣᫬㛫
᪤ᐃ್/᭱኱್

࣓ࣔࣜࢧ࢖ࢬ

72 ᫬㛫/720 ᫬㛫
256GBࣀ࣮ࢻᩘ

༨᭷

ಶูタᐃ

ͤ 2ࣀ࣮ࢻ௨ୖࢆ฼⏝ࡋࡓ୪ิᐇ⾜࡟ࡣMPIࡢ฼⏝ࡀᚲ⏝

Vol.55 No.2 2022-4［共同研究成果］
更新されたベクトル機での日本域を対象とした長期領域再解析の実行…………………………福井真・曽我 大輝・小原 涼太・臼井健
瀧口 海人・小野 佳祐・廣瀬 大河・松島 沙苗
白川 栄一・伊藤 純至・岩崎 俊樹・山崎
剛
斉藤 和雄・瀬古弘
熱損失を考慮した Flamelet approach に基づく燃焼シミュレーションとNO 生成の予測 ……………………松川 嘉也・小澤 龍磨・赤尾津翔大・松下
青木 秀之・Weeratunge Malalasekera

洋介

Computational analysis of the body-force models forplasma actuator-assisted ﬂow control simulations …… Di Chen・Kengo Asada
Satoshi Sekimoto・Kozo Fujii
Hiroyuki Nishida
直交格子法による移動境界問題の解法に関する研究………高橋俊・佐々木大輔
山田 剛治

1

5

13

25

［お知らせ］
コンパイルのアップデートについて ………………………………………………………
【重要】AOBA-A で MPI ライブラリを利用するプログラムについての注意点 ………
令和 4 年度サイバーサイエンスセンター講習会のご案内 ………………………………

33
33
34

［大規模科学計算システム］
鍵ペアの作成とログイン方法 ………………………………………………………………
ストレージシステムの利用法 ………………………………………………………………
アプリケーションサービスの紹介 …………………………………………………………

35
42
48

［報 告］
東北大学と NEC、マテリアルズインテグレーションシステムによる航空機用複合材料開発
シミュレーションのサービスを提供開始〜スーパーコンピュータを活用し、高機能材料開発に貢献〜 …………………………

63［スーパーコンピュータ AOBA のお知らせより］
令和 4 年度利用負担金について ……………………………………………………………
令和 4 年度年度共同研究について …………………………………………………………
商用アプリケーションのバージョンアップについて ……………………………………
計算科学・計算機科学人材育成のためのスーパーコンピュータ無償提供制度について
民間企業利用サービスについて ……………………………………………………………
大規模科学計算システムの機関（部局）単位での利用について ………………………

64
67
68
69
70
70

執筆要項 …………………………………………………………………………………………

71

スタッフ便り ……………………………………………………………………………………

72