学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

jh210031-NAH

氷海船舶に作用する氷荷重のシミュレーション
渡辺 勢也（九州大学）
概要
北極海など氷で海面が覆われた海域を航海する氷海船舶は氷と衝突しながら進むため，船舶に作用する氷荷重の評価が極めて重要であり，航海の安全性に直結する．氷を張った水路で小型の模型船を移動させて氷荷重を評価する氷海水槽試験が有効であり，それを CFD で模擬する数値水槽試験の研究が盛んに行われている．本研究課題では，格子ボルツマン法による流体計算，フェーズフィールド法による自由表面計算，個別要素法による剛体計算を組み合わせた手法を開発し，GPU スパコンを用いた大規模計算により流体と氷の直接相互作用に基づく氷海船舶の高精度な解析の実現を目的とする．局所的にメッシュ解像度を上げる AMR（Adaptive Mesh Refinement）法を導入することで計算を効率化し，氷の形状や密接度，船舶のモデルや航行速度を変えた様々な条件のシミュレーションを実施した．氷荷重に対して先行研究の数値解析やモデル式と比較し妥当な計算結果であることを確認し，また，船舶の速度や氷海の密接度の増加に伴う氷荷重の増加を再現することができた．

1. 共同研究に関する情報

する氷荷重の評価が極めて重要であり，航海の安全(1) 共同研究を実施した拠点名性に直結する．CFD による数値水槽試験が有効であ

東京工業大学るが，一般的な船舶 CFD と異なり，流れの計算だけでなく船舶への氷の衝突，氷の運動，氷と氷の接触(2) 共同研究分野などの様々な現象を扱う必要がある．特に，流体と超大規模数値計算系応用分野

氷の相互作用を直接計算するには，氷の近傍で高解像度メッシュを用いた流体計算が要求され，必然的(3) 参加研究者の役割分担に大規模な流体計算となる．

渡辺勢也（九大）：コード開発，ジョブ実行

本課題では，大規模流体解析に適した完全陽解法

青木尊之（東工大）：GPU，大規模計算の助言である格子ボルツマン法に高解像度格子を局所的に胡長洪（九大）：解析条件の選定，結果の考察割り当てる AMR（Adaptive Mesh Refinement）法を導Mohamed M. Kamra（九大）：乱流モデルの検討

入し，複数 GPU 計算を実行する自由界面流れのコー

Jabir Al Salami（九大）：GPU 計算高速化

ドを開発する．それを用いて，氷海中を航行する船楼銘濤（九大）：界面捕獲手法の検討

舶の数値水槽試験を実施し，氷海船舶に作用する氷後藤優弥（九大）：ポスト処理

荷重の高精度な評価を実現することを目的とする．
氷のサイズや氷の量（密接度）を指定できる流氷モ

2.

研究の目的と意義

デルを導入することで，形状や密接度，船舶の速度地球温暖化による北極海の氷の減少に伴い，北極が氷荷重に与える影響を明らかにする．

圏の研究観測を行う観測船だけでなく，北極航路を3.

利用する外航商船が増加している．これらの船は流当拠点公募型研究として実施した意義

氷を割りながら進む砕氷船とは異なり，砕氷能力が本研究の流体シミュレーションに用いる格子ボルなくとも航海が可能な数メートルの小さい浮氷の集ツマン法は GPU と相性が良いため東京工業大学のまった海域を選択して航海する．船舶は氷と衝突しTSUBAME3.0 の利用が有効である．研究副代表者のながら進むため，氷海船舶の設計では，船舶に作用

東工大・青木教授は GPU コンピューティングの第一1

学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

人者であり，かつ AMR 法や負荷分散など HPC 分野と求める．

での最先端の研究を進められている．JHPCN の枠組

船舶や氷の固体表面の境界条件として，格子点のみを使い，HPC 分野の研究者からスーパーコンピュ間に位置する壁を 2 次精度で扱える Interpolatedーティングの助言をいただくことで，初めて高精度bounce-back を用い，滑りなし境界条件を与える．

な氷海船舶のシミュレーションを実現できると考え格子点から壁までの距離𝑞に応じて以下の式から壁る．

で跳ね返る分布関数を以下のように計算する．
𝑓𝑖𝑗𝑘 (𝒙, 𝑡 + Δ𝑡)

4. 前年度までに得られた研究成果の概要

2𝑞𝑓̅̅̅̅
𝑖𝑗𝑘 (𝒙, 𝑡) + (1 − 2𝑞)𝑓̅̅̅̅𝑖𝑗𝑘 (𝒙 + 𝝃𝑖𝑗𝑘 Δ𝑡, 𝑡) +=新規課題であり該当しない．

5.

6𝜔𝑖𝑗𝑘 𝜌(𝝃𝑖𝑗𝑘 ∙ 𝒖wall ),𝑐2

6𝜔𝑖𝑗𝑘 𝜌(𝝃𝑖𝑗𝑘 ∙ 𝒖wall )1
2𝑞 − 1[𝑓̅̅̅̅]+𝑓 (𝒙 + 𝝃𝑖𝑗𝑘 Δ𝑡, 𝑡 + Δ𝑡),𝑖𝑗𝑘 (𝒙, 𝑡) +𝑐2
2𝑞 𝑖𝑗𝑘
{2𝑞0 ≤ 𝑞 < 0.5
0.5 ≤ 𝑞 < 1.0̅̅̅̅
ここで，𝒖wall は物体表面の移動速度であり，添字𝑖𝑗𝑘

今年度の研究成果の詳細

5.1 流体・船舶・海氷の連成解析手法の開発は反対向きの速度分布関数の方向を示す．物体に作本研究課題では，流体計算に格子ボルツマン法，自用する流体力は，物体境界における速度分布関数と由表面形状の計算にフェーズフィールド法，氷の計

物体の運動量交換に基づいて評価される．ある格子算に粒子ベースの剛体計算法を用い，これらを組み点における速度分布関数が物体に与える運動量は，合わせることで氷海船舶の計算手法を開発した．

物体表面で跳ね返る前後の速度分布関数から

格子ボルツマン法は圧力ポアソン方程式を解かない完全陽解法であり，計算アルゴリズムがシンプル

𝑭𝑖𝑗𝑘 (𝒙, 𝑡) = 𝝃𝑖𝑗𝑘 𝑓𝑖𝑗𝑘 (𝒙, 𝑡) − 𝝃̅̅̅̅
𝑖𝑗𝑘 𝑓̅̅̅̅
𝑖𝑗𝑘 (𝒙, 𝑡 + Δ𝑡)

でメモリ参照が局所的という特徴を持つため，大規模な流体シミュレーションに適した手法である．速と計算され，物体表面上で𝑭𝑖𝑗𝑘 を積分することで物体度空間で離散化したボルツマン方程式を用い，速度に作用する流体力𝑭fluid とモーメント𝑴fluid を求める．

分布関数𝑓𝑖𝑗𝑘 の時間発展方程式は𝑭fluid =𝑓𝑖𝑗𝑘 (𝒙 + 𝝃𝑖𝑗𝑘 Δ𝑡, 𝑡 + Δ𝑡) = 𝑓𝑖𝑗𝑘 (𝒙, 𝑡) + 𝛺𝑖𝑗𝑘 (𝒙, 𝑡)∑∑ 𝑭𝑖𝑗𝑘 (𝒙, 𝑡)𝒙∈boundary 𝑖𝑗𝑘である．ここで，𝑖𝑗𝑘は分布関数の速度方向で(𝑖, 𝑗, 𝑘) ∈𝑴fluid =∑𝒙∈boundary{−1, 0, 1}3 ，𝝃は分布関数の速度で，格子幅𝛥𝑥と時間[(𝒙 − 𝒙object ) × ∑ 𝑭𝑖𝑗𝑘 (𝒙, 𝑡)]𝑖𝑗𝑘

刻 み 𝛥𝑡 か ら 定 義 さ れ る 𝑐 = 𝛥𝑥/𝛥𝑡 を 用 い て 𝝃𝑖𝑗𝑘 =(𝑖𝑐, 𝑗𝑐, 𝑘𝑐)𝑇 である．衝突項𝛺𝑖𝑗𝑘 (𝒙, 𝑡)に安定性の高いキ

自由表面の形状はフェーズフィールド法を用いてュムラントモデルを用いることで，高レイノルズ数陰的に表現される．液相領域で𝜙 = 1，気相領域での流れに適用する．キュムラントモデルは速度分布

𝜙 = 0とし，界面では0 < 𝜙 < 1の値を取る．時間発展関数を統計量であるキュムラントに非線形変換し衝方程式には以下の保存型 Allen-Cahn 方程式を用いる．

突項を計算する．流体の密度と速度は，速度分布関数から

𝜕𝜙
1 − 4(𝜙 − 0.5)2+ 𝛻 ∙ (𝒖𝜙) = 𝛻 ∙ [𝑀 (𝛻𝜙 −𝒏)]𝜕𝑡𝑊𝜌(𝒙, 𝑡) = ∑ 𝑓𝑖𝑗𝑘 (𝒙, 𝑡)
ここで，𝑊は界面厚さ，𝑀はモビリティである．𝒏は𝑖𝑗𝑘

𝒖(𝒙, 𝑡) =1∑ 𝝃𝑖𝑗𝑘 𝑓𝑖𝑗𝑘 (𝒙, 𝑡)𝜌(𝒙, 𝑡)界面の単位法線ベクトルであり，界面の符号付き距離関数から計算する．本研究の自由表面流れの計算

𝑖𝑗𝑘

では，液相領域のみを格子ボルツマン法で計算し，2

学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

自由界面で境界条件を課す．保存型 Allen-Cahn 方程の重心の位置，𝑛𝑘 は剛体を形成している球形粒子の式の計算に必要な気相領域における速度場は数である．これらの式を 1 次オイラー法で時間積分Velocity Extension 法により補外している．また，界面することで剛体の並進運動と回転運動を計算する．

近傍の格子点において，気相領域から移流してくる本研究では計算高速化のため，図 1 に示す Octree

速度分布関数をデータ構造を用いて氷片や船舶，海面，渦の近傍に高解像度格子を適合させた細分化格子を生成する．

eq𝑓𝑖𝑗𝑘̅̅̅̅ (𝒙, 𝑡 + Δ𝑡) = −𝑓𝑖𝑗𝑘 (𝒙, 𝑡) + 𝑓𝑖𝑗𝑘 (𝜌b , 𝒖b )eq+ 𝑓ijk̅̅̅̅ (𝜌b , 𝒖b )GPU でのスレッド並列化を容易に行うために，octree
構造で細分化された小領域に対して8 × 8 × 8の均一格子を割り当てるブロック構造 AMR 法を用いる．と設定する．ここで，添字eqは局所平衡状態の分布大規模計算を行うために，空間充填曲線のひとつで関数，𝑓𝑖𝑗𝑘
̅̅̅̅ は𝑓𝑖𝑗𝑘 の反対向きの分布関数，𝜌b は基準とあるモートン曲線による領域分割法を用いた複数なる密度（大気圧），𝒖b は界面の速度を示す．

GPU 実装を行う．1 つの MPI プロセスに 1 台の GPU

複雑形状の物体同士の衝突判定と接触力の計算にを割り当てる Flat-MPI で並列化を行う．GPU 間通信は，粒子ベースの剛体シミュレーション手法を用いのチューニングとして計算と通信のオーバーラップる．船舶や氷片の形状は物体表面に配置された複数法を導入する．図 2 に示すように，各小領域の計算の球形粒子で表現される．固体表面に配置された球ブロックを隣接領域のデータが必要な外側ブロック

形粒子は，他の物体との衝突検出と接触力の計算にと通信が必要ない内側ブロックに分ける．GPU 間通利用される．計算安定性向上のため，物体同士の衝信と同時に内側ブロックの計算を行うことで，通信突時に物体間に LBM の数メッシュ分の隙間が空くよ

オーバーヘッドを隠蔽し性能を向上させる．

うに，球形粒子にはある程度の大きさを持たせる．
接触力の計算には，個別要素法のバネとダッシュポットのモデルを用いる．バネは粒子の食い込み深さに比例した反発力を，ダッシュポットは粒子の相対速度に応じた減衰力を粒子に与える．粒子間の摩擦

図 1 木データ構造に基づくブロック構造 AMR 格子を表現する接線方向のモデルには摩擦スライダーを導入し，クーロンの摩擦法則に従い接線方向の力に上限値を設定する．各球形粒子𝑖に作用する接触力𝑭𝑖の総和を計算し，剛体に加わる合力と重心周りのトルクを求める．
剛体𝑘の並進と回転の運動方程式は，剛体の質量を𝑚𝑘 ，並進速度を𝒗𝑘 ，角速度を𝝎𝑘 ，角運動量を𝑳𝑘 ，慣図 2 AMR 格子での計算と通信のオーバーラップ性テンソルを𝐼𝑘 とすると，5.2 TSUBAME3.0 での実行性能測定
𝑛𝑘

東京工業大学の TSUBAME3.0 において弱スケー

d𝒗𝑘
𝑚𝑘= ∑ 𝑭𝑖 + 𝑚𝑘 𝒈 + 𝑭fluidd𝑡

リングと強スケーリングを測定した．GPU には Tesla

𝑛𝑘P100 を用い，1 ノードあたり 4MPI プロセスと 4GPU

d𝑳𝑘= ∑((𝒙𝑖 − 𝒙𝑘 ) × 𝑭𝑖 ) + 𝑴fluidd𝑡を割り当てた．スケーリングは 16GPU から 256GPU

𝝎𝑘 = 𝐼𝑘−1 𝑳𝑘で評価した．オーバーラップ手法によるスケーリン
グ向上の効果を確認するため，オーバーラップ手法となる．ここで，𝒙𝑖 は剛体を形成している球形粒子𝑖ありとなしの両方の場合で計測した．

3

学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

弱スケーリングの結果を図 3 に示す．1 台の GPU

強スケーリングの結果を図 5 に示す．格子点数はあたり 25,225 ブロック，12,915,200 格子点を割り当206,643,200 であり，弱スケーリングと同様に最細格て，最も細い格子で 256 ステップの計算に要する時子で 256 ステップの計算に要した時間を示している．

間を測定した．AMR 法の格子生成はでの 128 ステッ

16GPU から 64GPU までは並列数の増加により実行プに一度の頻度で行った．GPU 数の増加に伴い実行時間を短縮できているが，それ以上の並列数では実時間が長くなっており，理想的な弱スケーリングは行時間はほぼ横ばいとなっている．計算時間と通信達成できていない．格子ボルツマン法のステンシル時間の内訳を見ると，並列数の最も少ない 16GPU の計算の時間は GPU 数によらず一定であるが，AMR条件で通信時間は計算時間とほぼ同じである．GPU法のメッシュ生成時間は並列数の増加に伴い長くな数の増加に伴い通信時間の割合は増え，256GPU ではっている．これは，格子生成の処理に関しては MPI全体の実行時間の約 7 割が通信時間となっている．

による分散並列が行えていない箇所があるためであ

このことから，本研究課題で開発したコードは通信り，問題規模の増加に伴い格子生成時間が長くなっがボトルネックになっており，効率的な計算を行うてしまうためである．また，どの条件でもオーバーには利用する GPU 数を少なくする必要があること

ラップ手法を導入したほうが計算時間は短く，例えが確認できた．また，強スケーリングにおいても，オば，16GPU ではオーバーラップ手法を用いた場合ーバーラップ手法により実行性能を改善できた．

1.29 倍の高速化となった．図 4 は 16GPU における計算時間と通信時間の内訳をサブルーチンごとに示している．オーバーラップ手法が有効なサブルーチンと効果のないサブルーチンがあることがわかる．オーバーラップ手法を用いた場合，キュムラント型の格子ボルツマン法の計算時間は 3.6 秒であり，計算全体の 30.9 秒に対して約 11%である．自由表面やVelocity Extension 法の処理に時間がかかっていることがわかる．

図 5 強スケーリングの性能測定結果

5.3 模型船に作用する氷荷重の評価開発コードの検証を兼ね，砕氷船の模型船舶を用いて氷荷重の評価を行った．模型船には，砕氷実験用の単純形状模型を用い，全量が 3.5 m，船幅が 0.7m の六角錐台形形状で，実船の 25 分の 1 スケールである．船体両端の船首の両方が実験に利用できる形状であり，片側は船首傾斜が 30 度（船首 A），も

図 3 弱スケーリングの性能測定結果

う一方は 45 度（船首 B）である．計算領域を 1.39 m× 11.15 m ×1.39 m とし，水深を 1.1 m と設定した．
FSICR（Finnish-Swedish Ice Class Rules）のガイドラインに基づき，小氷片が密集する brash ice channel を模擬して水面には一辺の大きさが 41 mm，密度が 950kg/m3，摩擦係数 0.1 の立方体形状の氷を等間隔に配置する．
水面に配置する氷の個数を 2034 個，3136 個，4096 個の 3 パターンとし，それぞれ氷の密接度は

図 4 弱スケーリングの計算時間内訳（16GPU）

4

学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

50%，68%，88%である．また，船舶の速度を 0.514

船体下部に潜り込む氷片が多くなり，より深くまでm/s，0.771 m/s，1.028 m/s の 3 種類の条件で設定し，沈むことがわかった．また，船首 A（30 度）と船首船舶の速度や氷の密接度が氷荷重に与える影響を調B（45 度）の結果を比較すると，船首傾斜の浅い船首査する．最細格子の格子幅は 2.72 mm であり，氷のA のほうが氷片が船体下部に潜り込む領域が広いこ

一辺の長さ 41 mm に対して約 15 メッシュ，計算領とがわかる．

域に対しては 512 × 4096 × 512 メッシュ相当である．
最細格子は自由表面，船舶，氷，渦の近傍に配置される．計算格子点数は計算初期で約 1.6 億であり，計算には NVIDIA Tesla P100 を 32 台用いた．計算はどのケースも 24 時間以内に完了している．
船舶の航行速度を 0.514 m/s，密接度を 68%，船首(a) 0.514 m/s(b) 1.028 m/s

A の条件における計算結果のスナップショットを図

図 7 航行速度の氷片挙動への影響（密接度 68%，船6 に示す．上のパネルが計算領域の上面図であり，自首 A）

由表面と氷，船舶を示している．下のパネルは側面図であり，氷と船舶，計算領域の中央における速度プロファイルを表示している．船舶と衝突した氷片は船体下部に潜り込み，その後，船舶の横側に浮き上がっていく様子が確認できる．それにより船舶が(a) 50%(b)88%移動した後方にはチャネルが形成されている．船舶

図 8 密接度の氷片挙動への影響（航行速度 0.514 m/s，の移動により水面に波が発生し，氷片も水面の揺れ船首 A）に合わせて上下運動する様子が確認された．

(a) 船首 A(b)船首 B

図 9 船首傾斜の氷片挙動への影響（航行速度 0.514m/s，密接度 68%）

模型船舶に加わる氷荷重の時刻歴の一例として，
図 6 航行速度 0.514 m/s，密接度 68%，船首 A での航行速度 0.514 m/s，密接度 68%，船首 A における結計算結果

果を図 10 に示す．横軸は船舶の移動距離，縦軸は氷荷重をログスケールで示す．赤い線は時間平均した氷片の運動に対する船舶の速度，氷の密接度，船氷荷重である．接触荷重は氷の衝突個数に大きく影

首傾斜の違いの影響を順に図 7 から図 9 に示す．上

響され，氷荷重は 10-3 N から 102 N のオーダーで大図が横から，下図が船底側からの可視化図である．

きく変動することが確認できた．

航行速度が異なる場合でも氷片の船体下部への潜り航行速度が氷荷重の時間平均値に与える影響を図

込みの様子はあまり変わらないことがわかる．航行

11 の上図に示す．航行速度の増加に伴い氷荷重も大速度が大きいほど船舶の移動により発生する波が高きくなる傾向が確認できる．一般的に接触力は速度くなるため，氷片も波に合わせてより大きく上下運の 2 乗に比例するが，本研究の結果では，氷荷重は動する様子が確認できた．氷の密接度が増加すると速度の 2 乗よりも低い増加傾向となった．これは，

5

学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

船舶の速度が増加するほど氷片の運動も大きくなり，接触回数が少なくなったためと考えられる．図 11 の下図は密接度と氷荷重の関係を示し，密接度の増加に伴い氷荷重が増加することがわかる．また，どの条件においても船首傾斜の大きい船首 B のほうが船首 A よりも抵抗が大きい結果となった．
スケールの 3 乗則にもとづいて実船スケールでの氷抵抗を推定し，FSICR の抵抗評価式に基づく抵抗と比較する．FSICR のアイスクラス 1A で航行速度を 5 ノットとすると，brash ice channel で船首 A では300 kN，船首 B では 362 kN の氷荷重となる．本研究の航行速度 0.514 m/s，密接度 88%の計算結果から実スケールの荷重を推定すると，船首 A では 188 kN，船首 B では 330 kN となる．抵抗評価式とオーダーは一致しているが，船首 A では抵抗を過小評価してい図 11 氷荷重の評価（上：航行速度の影響，下：密接る．船首 A と B の抵抗の比を比べると，FSICR の評度の影響）

価式では 1.21 に対して計算結果では 1.76 と異なる結果となった．また，先行研究の工学院大学の金野ら

5.4 Pack ice channel の氷荷重が行った氷片モデルのみの数値解析では，FSICR 評氷形状の船体抵抗への影響を調べるため，板氷で価式の 2 倍以上の氷抵抗を見積もっており，本研究覆われたチャネル（Pack ice channel）での計算を行い，のシミュレーションでより妥当な結果が得られた．

5.3 節の brash ice channel の結果と比較する．板氷モ
デルは厚さを 41 mm，縦と横の大きさを 150 mm とし，実スケールでは厚さ 1 m，直径が 3.75 m の流氷に対応する．水面には合計で 256 個の板氷モデルを配置し，その他の条件は 5.3 節と同じである．
航行速度 0.514 m/s における計算結果を図 12 に示す．船首の傾斜に沿って板氷が回転しながら船体下部に潜り込んでいく様子が確認できる．船首傾斜が大きい船首 B ほど，船底まで板氷が沈み込むのが早いことがわかる．
航行速度に対する氷荷重の平均値のグラフを図 13

図 10 氷荷重の時刻歴に示す．Brash ice channel と同様に，航行速度の増加に伴い荷重は大きくなり，船首傾斜が大きい船首 Bは抵抗が大きくなる傾向となった．また，Brash icechannel の場合よりも pack ice のほうが氷荷重は大きいことがわかった．

6学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

いほど氷抵抗が大きくなり，妥当な結果を得られている．模型船舶と比較すると，KCS コンテナ船型の氷荷重は半分以下となっている．これは船体下部に潜り込む氷片が少なく，船舶の両側に押しのけられる氷片が多くなったためと考えられる．

(a) 船首 A

図 14 Brash ice channel での KCS 船舶の計算結果(b) 船首 B
図 12 Pack ice channel の計算結果

図 15 Brash ice channel での KCS 船舶の氷荷重
図 13 Pack ice channel での氷荷重

6. 今年度の進捗状況と今後の展望
5.5 商船の船体データを用いた解析

本研究課題では，格子ボルツマン法による流体計実際の船舶形状への適用として，韓国 KRISO が公算，フェーズフィールド法による自由表面計算，個開している KCS（KRISO Container Ship）コンテナ船別要素法による剛体計算を組み合わせ，GPU スパコ型の船舶形状データを用いて，5.3 節の brash iceンを用いて流体と氷の直接相互作用に基づく氷海船channel での数値抵抗試験を実施した．計算条件は単舶の大規模シミュレーションを実現した．先行研究

純形状模型とほぼ同じであり，KCS は全長が単純模の数値解析や氷抵抗のモデル式と本計算結果を比較型よりも長いため計算領域を航行方向に拡大した．

し，氷荷重や氷の挙動に対し妥当な計算が行えていKCS コンテナ船型の 3DCAD データを，STL データることを確認した．氷の形状や密接度，船舶の航行を介して符号付き距離関数場に変換し，LBM の境界速度，船首形状を変えた多数のシミュレーションを条件や個別要素法の球形粒子の初期配置生成に利用行い，それらが氷荷重に与える影響を調査し，研究した．

目的は概ね達成できた．

KCS コンテナ船型での計算結果を図 14 に示す．

課題申請時の研究計画では，水槽試験で用いられKCS 船型は砕氷船模型モデルよりも喫水が浅いため， る割れた氷で覆われた氷海水路（broken ice channel）
氷片の船体下部への潜り込みが浅くなり，より後方も再現する予定であったが，複雑な形状の氷を水面まで氷片が水中を移動することがわかったに密に充填する初期配置作成が課題となり，broken

氷荷重の時間平均値を図 15 に示す．航行速度が速

ice channel でのシミュレーションは実施できなかっ7

学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書た．また，計算の安定性に関しても課題があり，計算なしに用いたよりも小さいサイズの氷片や 90%以上の密接度の条件では，氷片同士のめり込みなどの非物理(6) その他（特許，プレスリリース，著書等）

的な挙動が発生し計算は行えなかった．

[5]【受賞】渡辺勢也，若手優秀講演表彰，日本流体力割り当てられた計算資源の約 92%を使用した．

学会 第 35 回数値流体力学シンポジウム，2021 年

今後の展望としては，氷モデルや計算安定性を改12 月 16 日

善し氷海水槽試験を高精度で再現することで，実験では計測が困難な氷荷重の分布や船舶周囲の流れ場の氷挙動への影響の調査が期待できる．

7. 研究業績一覧（発表予定も含む．投稿中・投稿予定は含まない）(1)

学術論文 （査読あり）[1] Seiya Watanabe, Takayuki Aoki: Large-scale flow
simulations using lattice Boltzmann method with AMR
following free-surface on multiple GPUs, Computer
Physics Communications, Volume 264, 107871, 2021
[2] Seiya Watanabe, Shintaro Fujisaki, Changhong Hu,
Numerical simulation of dam break flow impact onvertical cylinder by cumulant lattice Boltzmann
method, Journal of Hydrodynamics, Volume 33, pp
185-194, 2021

(2) 国際会議プロシーディングス （査読あり）
なし(3) 国際会議発表（査読なし）

なし(4) 国内会議発表（査読なし）[3]渡辺勢也，胡長洪，青木尊之，AMR 法を導入した
格子ボルツマン法による船舶-氷相互作用のシミュレーション，第 35 回数値流体力学シンポジウ
ム，2021 年 12 月
[4]渡辺勢也，自由表面流れに対する AMR 法を導入した LBM コードの開発と海洋工学への適用，第33 回 LBM 研究会，オンライン，2021 年 6 月(5) 公開したライブラリなど
8