巻頭言

Vol.19, No.1 ߸ͷൃ‫ʹץ‬౰ͨͬͯ‫౎ژ‬େֶֶज़৘ใϝσΟΞηϯλʔɹ‫ڇ‬ౡ লຊ߸Ͱ͸ɼ‫౎ژ‬େֶֶज़৘ใϝσΟΞηϯλʔʹ͓͍࣮ͯࢪ͞Εͨ 2019 ೥౓ͷ‫ڞ‬ಉ‫ڀݚ‬ใࠂʹ͍ͭͯಛू͍ͨ͠·͢ɽ2019 ೥౓͸ɼ13 ݅ͷʮएखɾঁੑ‫ऀڀݚ‬঑ྭ࿮ʯɼ1 ݅ͷʮେ‫ن‬໛‫ࢧࢉܭ‬ԉ࿮ʯɼ
ͦͯ͠ 3 ݅ͷʮϓϩάϥϜߴ౓Խ‫ڞ‬ಉ‫ڀݚ‬ʯʹΑΔ‫ڞ‬ಉ‫͞ࢪ࣮͕ڀݚ‬Ε·ͨ͠ɽֶज़৘ใϝσΟΞηϯλʔͰ͸ɼεʔύʔίϯϐϡʔλ‫ڞ‬ಉ‫ͯ͠ͱ౓੍ڀݚ‬ɼʮएखɾঁੑ‫ऀڀݚ‬঑ྭ࿮ʯɼʮେ‫ن‬໛‫ࢧࢉܭ‬ԉ࿮ʯɼʮϓϩάϥϜߴ౓Խ‫ڞ‬ಉ‫ڀݚ‬ʯͱ͍͏̏छྨͷ‫ࢧڀݚ‬ԉࣄ‫ۀ‬Λ࣮ࢪ͍ͯ͠·͢ɽ՝୊ืू͸೥౓։࢝લ‫ߦʹޙ‬ΘΕ·͢ɽ·ͨɼ೥౓్தͰ௥Ճืू͕࣮ࢪ͞ΕΔ৔߹͕͋Γ·͢ɽԠื͞Εͨ՝୊ͷ࠾൱͸ɼεʔύʔίϯϐϡʔλ‫ڞ‬ಉ‫اڀݚ‬ըҕһձʹͯ৹ࠪ͞Ε·͢ɽຊ‫ڞ‬ಉ‫΁౓੍ڀݚ‬ͷԠืํ๏΍ਃࠐΈ࣌‫ظ‬౳ͷৄࡉʹ͖ͭ·ͯ͠͸ɼhttp://www.iimc.kyoto-u.ac.jp/ja/
services/comp/support/kyodo/ Λ͝ࢀরԼ͍͞ɽ
্‫ه‬ͷ‫ڞ‬ಉ‫౓੍ڀݚ‬ͷ͏ͪɼʮएखɾঁੑ‫ऀڀݚ‬঑ྭ࿮ʯͷԠืࢿ֨͸ɼ(1) 40 ࡀະຬͷएख‫ऀڀݚ‬(ֶੜΛ‫ؚ‬Ή, ੑผ͸໰Θͳ͍)ɼ͋Δ͍͸ɼ(2) ঁੑ‫( ऀڀݚ‬೥ྸ͸໰Θͳ͍)ɼͱ͞Ε͍ͯ·͢ɽ͜ͷʮएखɾঁੑ‫ऀڀݚ‬঑ྭ࿮ʯͰ͸ɼεʔύʔίϯϐϡʔλΛར༻͢Δ͜ͱͰֶज़తʹΠϯύΫτ͕͋Δ
੒ՌΛ૑ग़Ͱ͖Δ՝୊ʹରͯ͠ɼ‫ػࢉܭ‬ར༻ෛ୲ۚͷશֹ·ͨ͸Ұ෦Λຊηϯλʔ͕ෛ୲͍ͯ͠·͢ɽ͞Βʹɼʮएखɾঁੑ‫ऀڀݚ‬঑ྭ࿮ʯ͸ɼJHPCN (ֶࡍେ‫ن‬໛৘ใ‫ج‬൫‫ڞ‬ಉར༻ɾ‫ڞ‬ಉ‫ )఺ڌڀݚ‬ͷ‫׆‬ಈͷҰ‫ͱ؀‬Ґஔ͚ͮΒΕ͓ͯΓɼ༏Εͨ՝୊͸ɼJHPCN ͷ๖ժ‫ڞܕ‬ಉ‫΁ڀݚ‬ਪન͞Ε·͢ɽਪન͞
Εͨ՝୊͸ɼJHPCN γϯϙδ΢ϜʹͯϙελʔൃදΛߦ͏͜ͱ͕ՄೳͱͳΔͳͲͷಛయ͕͋Γ·͢ͷͰɼੋඇ͜ͷ‫ڞ‬ಉ‫౓੍ڀݚ‬Λ͝‫͍ͩ͘͞༻׆‬ɽຊ߸ʹ‫͞ࡌܝ‬Εͨ 13 ݅ͷʮएखɾঁੑ‫ऀڀݚ‬঑ྭ࿮ʯͷ‫ڀݚ‬੒ՌใࠂͰ͸ɼෳࡶྲྀମ͓Αͼྲྀମߏ଄࿈੒໰୊ɼϦϞʔτηϯγϯάʹΑΔ‫ڥ؀‬ϞχλϦϯάɼ΢Ϛͷ‫܈‬Εͷ‫ݸ‬ମؒͷྗֶϞσϧ‫ࢉܭ‬ɼίϯΫϦʔτͷഁյྗֶɼྔࢠԽֶɼ֤छͷ෼ࢠಈྗֶ‫ࢉܭ‬ɼ݁থੜ੒ͳͲʹؔ͢Δଟ༷ͳ‫ߦ͕ڀݚ‬ΘΕɼશࠃ 7 େֶͱࠃཱߴ౳ઐ໳ֶߍ 1 ߍͷएखɾঁੑ‫ऀڀݚ‬ͷ‫ڀݚ‬Λαϙʔτ͍ͨ͠·ͨ͠ɽͳ͓ɼ
JHPCN ͷ๖ժ‫ڞܕ‬ಉ‫ʹڀݚ‬͸ 10 ݅ͷ՝୊͕࠾୒͞Ε·ͨ͠ɽʮେ‫ن‬໛‫ࢧࢉܭ‬ԉ࿮ʯ͸ɼλΠϓ A Ͱ͸ 128 ϊʔυ × 2 िؒɼλΠϓ B Ͱ͸ 64 ϊʔυ × 2 िؒ·Ͱͷ‫ݯࢿࢉܭ‬ɾ‫ؒظ‬Λ͝ར༻௖͘‫ڞ‬ಉ‫Ͱ౓੍ڀݚ‬ɼར༻ෛ୲ۚ͸ຊηϯλʔ͕શֹෛ୲͢Δ΋ͷͰ
͢ɽ2019 ೥౓ͷʮେ‫ن‬໛‫ࢧࢉܭ‬ԉ࿮ʯͰ͸ɼಁաੑΛ༗͢Δฏߦฏ൘ؒͷཚྲྀ‫ߦ͕ࢉܭ‬ΘΕ·ͨ͠ɽ໿ 4 ԯηϧͷେ‫ن‬໛ྲྀମ‫ʹࢉܭ‬ΑΓɼཚΕ‫ͱ౓ڧ‬ฏ൘ͷಁաੑͷؔ܎ͳͲ͕ৄࡉʹ‫ݕ‬౼͞Εɼ༗༻ͳ੒Ռ͕ಘΒΕ·ͨ͠ɽ·ͨɼʮϓϩάϥϜߴ౓Խ‫ڞ‬ಉ‫ڀݚ‬ʯͰ͸ɼ‫౎ژ‬େֶͷεʔύʔίϯϐϡʔλΛάϧʔϓίʔε·ͨ͸ઐ༻ΫϥελίʔεͰ͝ར༻ͷ‫( ऀڀݚ‬JHPCN ͓Αͼ HPCI ʹͯ‫౎ژ‬େֶͷεʔύʔίϯϐϡʔλΛར༻͢Δ‫ऀڀݚ‬΋‫ؚ‬Ή) Λର৅ʹɼϓϩάϥϜͷ੍‫ߏޚ‬଄ɾσʔλߏ଄ͷվྑʹΑΔੑೳνϡʔχϯά΍ฒྻԽख๏ͷվྑɼ໰୊෼ׂɾෛՙ෼ࢄํࣜͳͲͷվྑʹΑΔେ‫ن‬໛‫ࢉܭ‬ϓϩάϥϜͷߴ౓ԽɾߴੑೳԽͳͲΛࢧԉ͍ͯ͠·͢ɽ2019 ೥౓͸ɼ(1) 3 ࣍‫߶ݩ‬ମόωϞσϧʹΑΔίϯΫϦʔτߏ଄෺ͷม‫ܗ‬ɾഁյγϛϡϨʔγϣϯɼ(2) ਫΛਪਐࡎͱ͢Δ௒খ‫ܕ‬ϚΠΫϩ೾์ిࣜϓϥζϚ‫ݯ‬ͷ 3 ࣍‫ݩ‬ཻࢠ‫ࢉܭ‬ɼ(3) ௡೾ʹΑΔ໦଄Ո԰ͷ‫נ‬᛽ੜ੒γϛϡϨʔγϣϯͷߴ౓Խɼͱ͍͏ 3 ݅ͷϓϩάϥϜ
ߴ౓ԽࢧԉΛߦ͍·ͨ͠ɽ͜ΕΒͷ‫ڞ‬ಉ‫Ͱڀݚ‬͸ɼεϨουฒྻԽ΍ೖग़ྗ࣌ͷόοϑΝϦϯάɼϕΫτϧԽ࠷దख๏ͷಋೖɼ·ͨɼྖҬ෼ׂʹΑΔཻࢠιʔςΟϯάͷ‫ݕ‬౼ͳͲ͕࣮ࢪ͞Εɼ‫ࢉܭ‬଎౓Λߴ଎Խ͢Δͱͱ΋ʹɼϓϩάϥϜͷվྑҊͳͲ͕ࣔ͞Ε·ͨ͠ɽࠓ‫ޙ‬΋օ༷ͷ‫ڀݚ‬ɼ‫ڭ‬ҭʹ͝‫͚͍ͩͨ༻׆‬ΔΑ͏ʹηϯλʔ‫ڭ‬৬һ΋ਚྗ͍͖ͯ͠·͢ͷͰɼ͝ར༻ɾ͝ࢧԉͷ΄ͲɺΑΖ͓͘͠‫͍ئ‬ਃ্͛͠·͢ɽ

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 1

スーパーコンピュータ共同研究制度（若手・女性研究者奨励枠）研究報告

液体攪拌操作における翼端渦と多相分散の関係性の解明
山本 卓也
東北大学大学院環境科学研究科先端環境創成学専攻、工学研究科金属フロンティア工学専攻

1 緒言

ケードによって翼端渦から変化した渦が気泡を微細化するかに関しては理解されていない。気泡と液体の攪拌操作は化学、金属、食品工場等各種

同程度のサイズの乱流渦が分裂に寄与するのであ産業で用いられるため、その原理を理解することれば、エネルギーカスケードによって生じた渦には応用面で重要である。この液体をかき混ぜるとよって分裂される可能性も否定できない。

いう操作は非常に単純な操作ではあるものの、そ本研究では、上記メカニズムを解明する第一段の混合中の液体中で発生する現象は非常に複雑で階として、翼周辺の翼端渦構造を数値解析によっある。また、液体だけではなく、気体や相分離すて解明した。

る液体が存在している条件では連続相内の流動に加えて分散相との関係性も重要となる。

2 数値解析手法

上記の攪拌操作に対して、これまでも数多くの研究が行われてきた。例えば、気泡吹き込みを伴

内径 192 mm, 高さ 200 mm の円筒容器内にアルう攪拌操作では、Hinze らの気泡径に対する相関

ミニウム産業界で利用される 6 枚軸流翼である攪

式[1]に始まり、
Zhou and Kresta [2]や Paul ら[3]によ拌翼[8]を挿入し、300, 400 rpm で回転させた。こって review されているように、気泡径に対する大の際の攪拌 Reynolds 数は 18000, 24000 である。

量の相関式が実験的に作成された。Zhou and

支配方程式は Navier-Stokes 式、連続式であるとKresta [2]は平均エネルギー散逸率ではなく、局所した。攪拌翼の回転は Sliding grid 法を利用するこ最大エネルギー散逸率でスケーリングするとよりとで直接回転させ、sliding 界面では cyclic AMI 境

正確に示すことができると示した。一方で、界条件を利用した。また、容器上面の境界を除いAndersson and Andersson [4]のように、乱流中の気て壁面条件を課し、上面では圧力勾配に応じた自泡や液滴の分裂モデルは気泡と乱流渦の関係性か

由流入出条件を利用した。LES によって乱流をモらモデル化されてきた。気泡や液滴と同程度のサデル化し、一方程式型サブグリッドスケールモデ

イズの乱流渦が分裂に寄与するとされ、このようルを利用した。

な分裂モデルを元にポピュレーションバランスモ

空間の離散化手法として 2 次精度線形補間法、デルが攪拌槽に拡張されている[5]。しかしながら、時間の離散化手法として 2 次精度 backward 法、速攪拌槽内での気泡分裂と乱流渦の関係性を実験、度圧力連成手法として PIMPLE アルゴリズムを利

もしくは数値解析を利用して直接的に捉えたもの用した。計算格子は翼付近で細分化し、総計算格は少ない。著者らは気泡分裂と渦の関係性を明確

子数は約 1400 万であった。
本計算はオープンソーにするため一連の研究[6,7]を行い、攪拌翼後方かスである OpenFOAM を利用した。数値解析コーら発生する翼端渦対間のせん断流によって気泡がドの検証は既往の研究[8]で行い、攪拌翼からの吐分裂することを見出した。この研究では翼端渦対出流を正確に計算できた。

のせん断流によって段階的に小さい気泡へと分裂される現象（カスケード分裂）も確認されたが、3 結果と考察

翼端渦によって分裂した後にもう一段階分裂するメカニズムは定かでない。特に、エネルギーカス

2 − Academic Center for Computing and Media Studies, Kyoto University

流れ場を発達させた後の翼端渦形状を Fig. 1 に示す。翼端渦は速度勾配テンソルの第二不変量のとが難しくなり、気泡分裂に寄与するのが難しく

等値面から渦を描画した。周方向渦度を利用してなると考えられる。

回転方向により渦を色分けした。翼端渦は翼後面で渦対として生成し、吐出流に乗って周方向へと4 結論移動する。攪拌槽側壁と軸の中央付近で翼端渦は散逸、構造が変化し、最終的には壁面に衝突する

攪拌操作中の翼端渦と多相分散の関係性についことでその向きを変え、上下縦に伸びた渦へと変て解明するために、攪拌操作時の翼端渦構造を調化し、最終的に散逸する。また、回転数を増加さ査した。攪拌操作中では攪拌翼後面で翼端渦が形せることでより複雑な渦構造へ変化し、翼端渦の成し、それが吐出流に乗って径方向へと広がる。

長さが伸長する。

その広がる過程で翼端渦は構造が変化し、散逸するのに伴い、螺旋状の渦へと変化することが分かった。

Fig. 2 Separated trailing vortices at 400 rpm.

参考文献
[1] J. O. Hinze, AIChE J. 1, 289-295 (1955).
Fig. 1 Trailing vortices during mechanical stirring: (a)at 300 rpm and (b) 400 rpm.
翼端渦の特徴を評価するために、400 rpm の条
件における回転方向によって渦を分離し表示した
結果を Fig. 2 に示す。この際計算結果を見やすくするため、半径 35-90 mm の領域に存在する渦のみを示す。特徴的な構造として、翼端渦は翼付近から翼後方へ向かって真っ直ぐの渦が生成するが、距離が翼から離れるに従い、散逸するのに伴い螺旋状の渦が生成する。また、回転方向の異なる渦が交互に存在する。これは、翼付近で発生した翼端渦が交互に渦対として発生し、その後翼端渦が散逸することで、回転方向の異なる渦群が交互に存在するようになるからである。
先行研究によると、攪拌操作中に吹き込んだ気泡は主に翼から発生した翼端渦対の渦対間のせん[2] G. Zhou and S. M. Kresta, Chem. Eng. Sci. 53,
2063-2079 (1998).
[3] E. L. Paul et al., Handbook of industrial mixing
science and practice (Wiley, Hoboken, 2004).
[4] R. Andersson and B. Andersson, AIChE J. 52,
2031-2038 (2006).
[5] S. Solsvik et al., Ind. Eng. Chem. Res. 55,
2872-2882 (2016).
[6] T. Yamamoto et al., Chem. Eng. Sci. 197, 26-36(2019).
[7] T. Yamamoto and S. V. Komarov, Chem. Eng. Sci.
207, 1007-1016 (2019).
[8] T. Yamamoto et al., J. Mater. Process. Technol. 261,
164-172 (2018).

謝辞

断によって分裂する[7]。一方で、今回解明したよ
うに、攪拌翼から離れるにつれ渦構造が変化する
が、渦の回転方向の揃った領域が発生する。このため、渦対間のような強いせん断力が発生するこ本研究制度（若手・女性研究者奨励枠）を活用させて頂いたことを、この場を借りて厚く御礼申し上げる。

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 3ք໘ͷຎࡲ઀৮Λߟྀͨ͠ଛইϞσϧʹΑΔ
మ‫ے‬ίϯΫϦʔτͷ 3 ࣍‫ݩ‬ഁյγϛϡϨʔγϣϯ૬അ ༔ਓ ∗∗

ҵ৓େֶେֶӃཧ޻ֶ‫ڀݚ‬ՊࣾձΠϯϑϥγεςϜՊֶઐ߈1 ͸͡Ίʹ

450400

౔໦ߏ଄෺Λߏ੒͢Δమ‫ے‬ίϯΫϦʔτ͸ɼࡐྉ450ಛੑͷҟͳΔమ‫ͱے‬ίϯΫϦʔτΛ૊Έ߹Θͤͨߏ଄Ͱ͋Γɼమ‫ے‬ίϯΫϦʔτͷྗֶੑೳ͸ɼమ‫ͱے‬ί

200

y
1200

ϯΫϦʔτؒͷ෇ண‫ڍ‬ಈʹେ͖͘Ө‫͞ڹ‬ΕΔɽమ‫ے‬ͱίϯΫϦʔτؒͷ෇ண͸ɼओʹք໘ͷ೪ணྗɼຎ

x
z

150

ࡲ఍߅ྗɼమ‫ے‬ͷઅʹΑΔࢧѹ఍߅ྗʹΑͬͯ੒Γ
Interfaceཱ͍ͬͯΔɽమ‫ے‬ίϯΫϦʔτͷ਺஋ղੳͰ͸ɼ͜ΕΒͷ෇ணཁҼΛؒ઀తʹϞσϧԽ͍ͯ͠Δํ๏͕΄ͱΜͲͰ͋Γɼ௚઀ղੳʹ൓өͤͨ͞ํ๏͸΄ͱΜͲͳ͘ɼమ‫ے‬ͷ‫ܗ‬ঢ়ͷϞσϧԽͷΈͰమ‫ے‬ίϯΫϦʔτͷ෇ணੑೳΛධՁͰ͖Δํ๏͸΄ͱΜͲͳ͍ɽ

Round bar ( 16)

Deformed bar (D16)

Cross sectionͦ͜Ͱຊ‫Ͱڀݚ‬͸ɼ෇ணཁҼΛ௚઀൓өͤͨ͞਺஋ղੳΛ࣮ࢪ͠ɼ෇ணੑೳͷҟͳΔ RC ͸Γͷഁյਤ 1: RC ͸Γͷ༗‫ݶ‬ཁૉϞσϧ‫ڍ‬ಈͷ࠶‫ݱ‬ੑΛ‫͢ূݕ‬Δɽ
ͯ͠๏ઢํ޲ͷͻͣΈ͕ ε′1 ≤ 0 ͷͱ͖ʹ઀৮ͱ൑ఆ

͠ɼ๏ઢํ޲ͷ߶ੑΛอ࣋ͤ͞Δ͜ͱͰ઀৮ྗΛද

2 ਺஋ղੳख๏‫͢ݱ‬Δɽ๏ઢํ޲ͷԠྗ σ′1 ͸࣍ࣜͰද͞ΕΔɽమ‫ʹے‬͸ɼඇઢ‫ܗ‬౳ํߗԽଇʹ‫ ͮ͘ج‬von-Mises
઼ੑϞσϧɼίϯΫϦʔτ͓Αͼք໘ʹ͸ɼഁյ໘σ′1 = C11 ε′1 + (1 − D)C12 ε′2 + (1 − D)C13 ε′3(2)ͷຎࡲ઀৮Λߟྀͨ͠ଛইϞσϧ [1] Λద༻͢Δɽຊ
‫Ͱڀݚ‬ద༻͢ΔଛইϞσϧͰ͸ɼഁյ໘ͷ։‫ڍޱ‬ಈͷຎࡲଇΛಋೖ͢ΔɽຎࡲԠྗ τf ͸࣍ࣜͰද͞ΕΔɽ

͸࣍ࣜͰද͞ΕΔɽ

σ = (1 − D)Cεഁյ໘ͷ‫ݻ‬ண΍‫׈‬ΓΛఆࣜԽ͢ΔͨΊɼΫʔϩϯ(1)͜͜Ͱɼσ ͱ ε ͸ϑΥʔΫτද‫ʹه‬ΑΔԠྗϕΫτϧͱͻͣΈϕΫτϧɼC ͸஄ੑ܎਺ϚτϦοΫεɼD͸ 0 ∼ 1 ͷ஋ΛऔΔଛইม਺Ͱ͋Δɽഁյ໘ͷຎࡲ઀৮ΛఆࣜԽ͢ΔͨΊɼഁյ໘ʹ௚

ަ͢Δ‫࠲ॴہ‬ඪ‫ܥ‬Λߟ͑ΔɽҎԼͰ͸ɼμογϡ෇͖ද‫)•( ه‬′ Λ‫࠲ॴہ‬ඪ‫ܥ‬ͷม਺ͱ͢Δɽഁյ໘ʹର

4 − Academic Center for Computing and Media Studies, Kyoto University

τf = µσn = µ|σ′1 |(3)͜͜Ͱɼµ ͸ຎࡲ܎਺ɼσn ͸ਨ௚ԠྗͰ͋Δɽ
ຎࡲԠྗͱഁյ໘ʹฏߦʹ࡞༻͢ΔͤΜஅྗͷେ͖͞ τ′s Λൺֱ͢Δ͜ͱͰɼഁյ໘ͷ‫ݻ‬ணͱ‫׈‬ΓΛ൑ఆ͢Δɽτf ≥ |τ′s | Ͱ͸ɼͤΜஅԠྗ͸࣍ࣜͰද͞ΕΔɽ
′τ′12 = Gγ12(4)′τ′31 = Gγ31

(5)

150

E



 y0

Q

b

200 GPa

0.3

380 MPa

200 MPa

40

E



k

Gf

0



Concrete

30 GPa0.2

20

0.1 N/mm

0.0001

Interface

200 MPa

0.2

20

0.4 N/mm

0.004

Steel

Load (kN)ද 1: ࡐྉύϥϝʔλ100

50Ana. (Round bar)Ana. (Deformed bar)

0.2
0

0

5
Displacement (mm)10

ਤ 2: ‫ؙ͓߯‬Αͼҟ‫ܗ‬మ‫ے‬Ϟσϧʹ͓͚Δՙॏ − มҐؔ܎

͜͜ͰɼG ͸ͤΜஅ߶ੑͰ͋Δɽ
ରͯ͠ɼτf < |τ′s | Ͱ͸ɼͤΜஅԠྗ͸࣍ࣜͰද͞

ΕΔɽ

γ′
′
τ′12 = (1 − D)Gγ12+ Dτf √ 12
′2
′2
γ12
+ γ31

(6)

γ′
′
+ Dτf √ 31
τ′31 = (1 − D)Gγ31
′2
′2
γ12
+ γ31

(7)
(a) RC beam with round bars: 40/300 (49.1 kN)

3 ෇ணੑೳͷҟͳΔ RC ͸Γͷഁյ‫ڍ‬ಈͷ࠶‫ݱ‬ੑ
3.1 ղੳ৚݅
ਤ 1 ʹ෇ணੑೳͷҟͳΔ‫ͱؙ߯‬ҟ‫ܗ‬మ‫ے‬Λ࢖༻͠(b) RC beam with deformed bars: 220/300 (120.0 kN)ͨ‫͛ۂ‬ഁյ‫ܕ‬ͷ RC ͸Γͷ༗‫ݶ‬ཁૉϞσϧΛࣔ͢ɽରশੑΛߟྀ͠ɼղੳൣғΛ 1/4 ͱͨ͠ɽཁૉʹ͸ɼ࢛
໘ମ 1 ࣍ཁૉΛ࢖༻͢Δɽమ‫ے‬ͷද໘‫ܗ‬ঢ়ΛৄࡉʹϞσϧԽ͠ɼమ‫ے‬ද໘ʹް͞ 0.5 mm ͷք໘૬Λઃ

1.0

0.0
Internal damage variableਤ 3: ‫ؙ͓߯‬Αͼҟ‫ܗ‬మ‫ے‬Ϟσϧʹ͓͚Δ಺෦ͷଛই෼෍͚ͨɽ֤Ϟσϧͷઅ఺਺͸໿ 15 ສɼཁૉ਺͸໿ 80
ສͰ͋Δɽද 1 ʹࡐྉύϥϝʔλΛࣔ͢ɽ‫੍ڧ‬มҐ

9 mm Λ 300 step Ͱ༩͑ͨɽమ‫ے‬ͷද໘‫ܗ‬ঢ়ͷϞσϧԽͷΈΛม͑ɼ2 έʔεͷղੳΛ࣮ࢪͨ͠ɽ

4͓ΘΓʹຊ‫Ͱڀݚ‬͸ɼ෇ணੑೳͷҟͳΔ‫ͱؙ߯‬ҟ‫ܗ‬మ‫ے‬Λ࢖༻ͨ͠‫͛ۂ‬ഁյ‫ܕ‬ͷ RC ͸ΓΛର৅ͱͨ͠਺஋ղ
ੳΛ࣮ࢪ͠ɼమ‫ے‬ͷද໘‫ܗ‬ঢ়Λม͑Δ͚ͩͰɼ෇ண

3.2 ղੳ݁Ռ

ੑೳͷҟͳΔ RC ͸Γͷഁյ‫ڍ‬ಈΛ࠶‫͖Ͱݱ‬Δ͜ͱਤ 2, 3 ʹՙॏ − มҐؔ܎ͱ಺෦ͷଛই෼෍Λࣔ͢ɽΛࣔͨ͠ɽҟ‫ܗ‬మ‫ے‬Λ༻͍ͨ RC ͸ΓͰ͸ɼయ‫ܕ‬తͳ‫͛ۂ‬ഁյͷྗֶ‫ڍ‬ಈΛࣔ͠ɼͻͼׂΕ͕෼ࢄͨ݁͠Ռͱͳͬ
ͨɽରͯ͠ɼ‫ؙ߯‬Λ༻͍ͨ RC ͸ΓͰ͸ɼҟ‫ܗ‬మ‫ے‬ͷઅʹΑΔࢧѹ఍߅ྗ͕ͳ͍ͨΊɼ࠷େՙॏ͕௿͘ɼ
͞Βʹ௿͍ՙॏϨϕϧͰ෇ண͕ফࣦ͢ΔͨΊɼ‫͛ۂ‬ͻͼׂΕͷΈ͕੒௕͢Δҧ͍͕‫ݟ‬ΒΕͨɽࢀߟจ‫ݙ‬[1] ૬അ ༔ਓɼं୩ ຑॹɿք໘ͷຎࡲ઀৮Λߟྀͨ͠ଛই
ϞσϧʹΑΔమ‫ے‬ίϯΫϦʔτͷഁյγϛϡϨʔγϣϯɼ౔໦ֶձ࿦จू A2(Ԡ༻ྗֶ)ɼVol.75, pp.I 165I 173, 2019.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 5

高効率有機系太陽電池の実現に向けた光機能性分子の構造と電子物性の
相関解明
東野 智洋
京都大学大学院工学研究科 分子工学専攻

1

緒言を向上させることができると考えた。そこで本研
究では、カテコールに電子求引性のチアゾール基

ポルフィリンは400–450 nmにSoret帯と呼ばれる強を縮環させた、チアゾロカテコール基を吸着基とい吸収と550–600 nmにQ帯と呼ばれる中程度の吸収して用いたポルフィリン色素 ZnPTC を設計・合をもち、増感色素として有望である。そして実際に、成した（Figure 1）[3]。

様々なプッシュ−プル型ポルフィリン色素が開発され、高い光電変換効率が実現されてきた。また、色

2

結果と考察

素を酸化チタン上に固定する吸着基はエネルギー
変換効率だけでなく、セルの長期耐久性においても重要な役割を果たす。ポルフィリン色素の多くではカルボン酸吸着基が用いられてきたが、カルボン酸
は酸化チタンから脱着しやすいことが太陽電池素子の耐久性を低下させる一因となっているため、より強く酸化チタンに吸着する吸着基の開発が望まれる。一方、カテコールは酸化チタンに強く吸着することが知られているが、カテコールの電子供与性のために、カテコールを用いた色素では光電変換効率が低い値に留まっていた[1,2]。そのため、カテコー
ルに電子求引性の置換基を導入することにより、カテコールの強い吸着能を活用しつつ、太陽電池性能
C6H13

C6H13

C6H13

C6H13

N
C8H17 O

C8H17 O

N
N

Zn

N
N
N

O C8H17 C8H17 O

O C8H17 C8H17 O

N
N

Zn

2.1 新規ポルフィリン色素の物性と太陽電
池性能評価
参照色素 YD2-o-C8 を用いて色素増感太陽電池
セル作製条件の最適化を行った結果、ZnPTC を用いたセルでの光電変換効率は 4.87%となった。
カルボン酸吸着基をもつ YD2-o-C8 を用いたセルでの光電変換効率は 9.93%であり、チアゾロカテ
コール基を用いた色素では変換効率が低下することがわかった。しかしながら、カテコール吸着基を用いた色素で光電変換効率は 2%を超えるようなものはなく[1]、特にポルフィリンとカテコール
を組み合わせた色素では 1%にも満たないため[2]、ZnPTC はカテコールを吸着基として用いた色素の太陽電池性能の最高値を更新した。これは、電子求引性のチアゾール基を縮環させたことにより、

N
N

O C8H17

電子注入効率が向上したためであると考えられる。
また、酸性あるいは塩基性溶液中における色素のO C8H17

酸化チタンからの脱着挙動を調べた結果、ZnPTC
はカルボン酸吸着基をもつ色素 YD2-o-C8 と比較

Sして脱着しにくく、強い吸着能を持つことを明ら

N

HO
OH
ZnPTC

O

OH

YD2-o-C8

Figure 1. Molecular structures of porphyrin dyes.

6 − Academic Center for Computing and Media Studies, Kyoto University

かにした（Figure 2）。このことから、カテコール
基に電子求引性置換基を導入した吸着基を適切に設計することにより、高い色素増感太陽電池性能と長期耐久性の実現が期待できると言える。

効率は 4.87%と中程度に留まったと考えられる。
このように、本計算結果は実験結果の理論的解釈の一助となった点で意義がある。
(a)

(b)

LUMO

Figure 2. Plots of the amount of porphyrin dyes on the
TiO2 surface as a function of immersion time for
ZnPTC (gray) and YD2-o-C8 (black) under THF/H2O(1:1 v/v) solution containing 28 mM acetic acid (solidlines) and 0.1 mM NaOH (dashed lines).
HOMO

2.2 理論計算
ポルフィリン色素の最安定化構造およびそのフロンティア軌道の電子構造について知見を得るた
めに、Gaussian09 プログラムを用いて密度汎関数
法 （ DFT ） に よ る 理 論 計 算 を 行 っ た（B3LYP/6-31G(d)）
。ポルフィリン色素 ZnPTC の構造は、カルボン酸吸着基をもつ色素 YD2-o-C8とほぼ同じであり、カテコール吸着基は分子構造にほとんど影響を与えないことがわかった。
また、
Figure 3 にポルフィリン色素 ZnPTC および
YD2-o-C8 の最高被占軌道（HOMO）と最低空軌
道（LUMO）における軌道分布を示す。ZnPTCのカテコール部分における LUMO の軌道分布は、
YD2-o-C8 のカルボキシル基部分における LUMOの軌道分布よりも小さくなっていた。吸着基部分の LUMO の軌道分布が小さいと、励起状態の色
素からの電子注入の効率が低くなってしまうことが知られているため、電子注入効率はカルボキシル基の場合よりも低下してしまうと考えられた。
加えて、ZnPTC のカテコール部分におけるFigure 3. Selected Kohn−Sham orbitals for (a) ZnPTC
and (b) YD2-o-C8 obtained by DFT calculations with the
B3LYP/6-31G(d) level. To simplify the calculations,
hexyl groups on the diarylamino groups and octyloxygroups on the meso-aryl groups were replaced withmethyl and methoxy groups, respectively.

3

参考文献

[1] E. L. Tae, S. H. Lee, J. K. Lee, S. S. Yoo, E. J. Kang,
K. B. Yoon, J. Phys. Chem. B 2005, 109, 22513–
22522.
[2] M. Adineh, P. Tahay, M. Ameri, N. Safari, E.
Mohajerani, RSC Adv. 2016, 6, 14512–14521.
[3] T. Higashino, H. Iiyama, Y. Kurumisawa, H.
Imahori, ChemPhysChem 2019, 20, 2689–2695.

HOMO の軌道分布は、YD2-o-C8 のカルボキシル
基部分における HOMO の軌道分布よりも大きくなっていた。吸着基部分の HOMO の軌道分布が大きいと、酸化チタンからの逆電子移動が起こりやすくなり、太陽電池性能が低下してしまうことが知られている。これらの要因により、光電変換

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 7

時系列リモートセンシングデータによる全球陸域統合環境モニタリング
堤田 成政
京都大学地球環境学堂

1 はじめに人工衛星をもちいた地球観測が 1980 年代から

デ ー タ に 着 目 し 、 ９ つ の 分 類 ク ラ ス （ 1.
Evergreen/Deciduous

Needleleaf

Trees;

2.

Evergreen/Deciduous

Broadleaf

Trees;

3.

活発化し、これまでに膨大な衛星観測データが蓄

Mixed/Other Trees; 4. Shrubs; 5. Herbaceous

積されてきた。それに伴い、断片的に撮影されてVegetation;

きたデータを経時的に、かつ全球的にくまなく利

Vegetation; 7. Non (Sparse) Vegetation; 8.

用可能とする観測データの整備が進められてきた。

Snow/Ice, and 9. Water-related ） に お い て本研究では、そのような大規模な時系列リモート

AVHRR CRD データと同様の空間解像度に再集

センシングデータをもちい、気候変動や人為的改計し、均一な土地被覆を示したピクセルを地上参変の陸域への影響を明らかにし、気候変動予測や照データとした。AVHRR CDR は、観測バンドを緩和策のための科学的知見を提供することを目的６つ有した全球をカバーする日次データであり、とする。具体的には、過去 35 年における森林増これを一年毎の平均、標準誤差を算出しまとめた。

加／減少、農地拡大、砂漠化といったグローバル

これらのデータを一年ごとに Random Forestに起こる地理的現象を統合的にモデル化する。全

Classification（Breiman, 2001）に適用した。一球レベルで過去の陸域観測から蓄積された

年ごとに独立した土地被覆分類図を比較すると、

AVHRR CDR （Vermote et al., 2019）の日次デ

分類誤差が伝搬してしまい（Post classificationータ（空間解像度約 5km）をもちいて、1982-2016

comparison problem）、正しく土地被覆変動をモにおける土地被覆変動のダイナミクスをモデル化

デル化できない。そのため、一年ごとの Random

する。使用データは米国の NOAA より無償提供

Forest Classification の分類確率に着目し、これされているプロダクトを使用するが、データ量がを時系列データとしてピクセルごとに多項ロジス

膨大であり通常のワークステーションでの処理にティック回帰を適用し、土地被覆変動をモデル化は限界がある。そのためスーパーコンピュータを

した。

6.

CultivatedandManaged利用し研究を実施する。

2 手法

3 結果
１年毎の土地被覆分類図を作成した。これらの過去 35 年間の土地被覆変動をモデル化するに分類精度は平均で 96%となった。対象とした 35

あたり、土地被覆の地上参照データを収集するこ

年 間 で は 「 Evergreen/Deciduous Needleleafとが必須である。ここでは、DISCover（1992-1993、Trees」が約 1,545,000km2 、「Cultivated and

空間解像度 1km）（Loveland et al., 2000）,

Managed Vegetation」が約 99,140,000km2 増加

GLC2000（2000、
空間解像度 1km）
（Bartholome

した一方で、その他の土地被覆は減少傾向を示し

et al., 2005）, GlobCover （2005、空間解像度

た。地域レベルでみると、サブサハラの地域で300m）（Arinoet et al., 2007）の過去の土地被覆「Non (Sparse) Vegetation」から「Herbaceous

8 − Academic Center for Computing and Media Studies, Kyoto University

Vegetation 」 へ 、 ま た ブ ラ ジ ル 南 東 部 で はの共同研究制度（若手奨励枠）を活用させていた「Evergreen/Deciduous Broadleaf Trees」からだきました。

「Cultivated and Managed Vegetation」への土地被覆変化が顕著であった。

6 参考文献[1] Arino, O., Gross, D., Ranera, F., Bourg, L.,
Leroy, M., Bicheron, P., … Weber, J.-L. (2007).
GlobCover: ESA service for Global Land Coverfrom MERIS. In I. 2007 (Ed.), Geoscience andremote sensing symposium, 2007 (pp. 2412–
2415). IEEE.
[2] Bartholome, E., & Belward, A. S. (2005).
GLC2000: A new approach to global land cover
mapping

from

Earth

observation

data.

International Journal of Remote Sensing, 26(9),1959−1977.
[3] Breiman, L. (2001). Random forests.
Machine Learning, 45(1), 5–32.
[4] Eric Vermote, Chris Justice, Ivan Csiszar,
Jeff Eidenshink, Ranga Myneni, Frederic Baret,
Ed Masuoka, Robert Wolfe, Martin Claverie
and NOAA CDR Program (2019): NOAA
Climate Data Record (CDR) of AVHRR Surface
Reflectance, Version 5. NOAA National Centers
for Environmental Information.
[5] Hansen, M. C., DeFries, R. S., Townshend, J.
図１：土地被覆図の例（1982, 1995, 2005, 2016

R., & Sohlberg, R. (2000). Global land coverclassifcation at 1 km spatial resolution using a年）

classifcation

4 おわりにtree

approach.

International

Journal of Remote Sensing, 21(6-7), 1331–1364.
[6] Loveland, T. R., Reed, B. C., Brown, J. F.,

本研究では、1982〜2016 年の土地被覆変動を

Ohlen, D. O., Zhu, Z., Yang, L., & Merchant, J.

地球観測データをもとにしてモデル化を試みた。

W. (2000). Development of a global land cover

全体的な傾向の把握には成功した一方で、細かなcharacteristics database and IGBP DISCover

課題として、(i) 時空間土地被覆参照データの収集

from 1 km AVHRR data. International Journal

方法の再検討、(ii) 本事例に適した機械学習モデ

of Remote Sensing, 21(6-7), 1303–1330.

ルの選択・構築、(iii) 類似の土地被覆分類プロダ
クトとの比較、があげられる。今後はこれらの課題を中心に引き続き研究をすすめ、論文化をすす
めたい。

5 謝辞
本研究では、京都大学学術情報メディアセンタ ー京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 9界面分光の分子動力学シミュレーション
城塚 達也
茨城大学工学部物質科学工学科

1 はじめに

面選択性（Domen et al., Nature 581, 411, 2020など）に注目しており、Low-pressure chemical

固体と液体の界面に代表される埋もれた界面はvapor deposition による実験との共同研究（論文

学術的・工業的に重要であるが、その分子レベル投稿中）により、（112）面が特に高い反応性を示での反応メカニズムの解析は実験的・理論的に難すことが実験的だけでなく理論的にも分かったた

しい。そこで本研究では、分子動力学シミュレー

め、この界面でのホールの安定性を調べる。シミションにより界面構造・ダイナミクスを解析し、ュレーションには anatase (101), (001), (112),

近年発展が著しい和周波発生（SFG）分光法を申(100), (110)面を用いる。Generalized gradient

請者がこれまで開発してきた手法により解析する。

approximation (GGA) に基づく汎関数では局在

特に、本研究では酸化チタン・水界面などにおけしたホールは安定化しないため、本研究ではる構造やダイナミクスを調べる。第一原理分子動

CDFT 法と Hartree-Fock 交換エネルギーを力学（MD）シミュレーションを用いて固液界面

10.5%混合した hybrid 汎関数を用いる。

をシミュレーションし、最終的に界面分光の１つである SFG 分光法を解析する。

2 手法

3 結果
ケイ酸の脱プロトン化反応をシミュレーショ
ン (Fig. 1 top) では、自由エネルギー摂動法を用

まず、第一原理（DFT）計算を用いて界面でのいて計算されたケイ酸の脱プロトン化の自由エプロトン移動をシミュレーションするための一般

ネルギー変化は 13.9 kcal/mol となり実験値の的な手法を開発した。[1]これまでのシミュレーシ

13.5 kcal/mol と良い一致を示した。一方で、3

ョンでは化学反応を記述する反応座標は結合長な点近似やガウス近似を用いると自由エネルギー

どの構造を用いることが多かった。しかし、この変化はそれぞれ 12.4 kcal/mol と 40.3 kcal/mol

ような座標は固液界面などの不均一な系ではレアとなりより正確な自由エネルギー摂動法によるイベントをサンプリングするために高コストなシ結果と異なる。
よって、
化学精度（約 1 kcal/mol）

ミュレーションが必要であった。そこで、本研究が求められる計算の場合、近似的手法ではなく自で は Constrained density functional theory

由エネルギー摂動法を用いる必要があるとわか(CDFT)法を用いてドナーとアクセプターのエネった。また、自由エネルギー変化は約 10 ps 程

ルギー差を反応座標に採用する。プログラムパッ度の分子動力学シミュレーションで収束したこ

ケージは cp2k を用いて、ケイ酸の脱プロトン化とから、CDFT 法により高精度に反応自由エネ

反応をシミュレーションした。自由エネルギー曲

ルギー曲面を計算できることがわかった。

面の計算では自由エネルギー摂動法を使用した。

アナターゼ型酸化チタン界面での計算 (Fig. 1

その後、アナターゼ型酸化チタン界面での実験

bottom) では、ヒドロキシ基がない場合ではアナとの共同研究を実施した。申請者は光触媒性能のターゼ型酸化チタン表面でホールは安定に局在化

10 − Academic Center for Computing and Media Studies, Kyoto University

しないのに対し、ヒドロキシ基がある場合ではヒドロキシ基にホールが安定に局在することがわかった。これにより。アナターゼ型酸化チタン表面

B, in press.
[2] T. Joutsuka and S. Yamauchi, Chem. Phys.
Lett., 741, 137108, 2020.

ではヒドロキシ基にホールが流れ、光触媒による[3] T. Nishikawa, K. Horiuchi, T. Joutsuka

酸化反応が進行すると考えられる。加えて、アナand S. Yamauchi, J. Cryst. Growth, doi:

ターゼ型酸化チタンでは(101)面が熱力学的に安10.1016/j.jcrysgro.2020.125849.

定であるが、(112)や(001)表面の方がより安定に[4] S. K. Reddy, R. Thiraux, B. A. W. Rudd, L.

ホールをトラップすることが判明した。この結果

Lin, T. Adel, T. Joutsuka, F. M. Geiger, H. C.

は(112)面における実験で見られた高い反応性を

Allen, A. Morita, and F. Paesani, Chem, 4 (7),支持している。 (001)表面での反応性はよくわか

1629–1644, 2018.

っていたが、(112)面の報告例は少ないため、今後[5] T. Joutsuka, T. Hirano, M. Sprik, and A.

の酸化チタンを使った光触媒材料開発にとって有Morita, Phys. Chem. Chem. Phys., 20 (5),

益な知見である。

3040–3053, 2018.
[6] T. Joutsuka and A. Morita, J. Phys. Chem.

4 おわりにC, 122 (21), 11407–11413, 2018.

界面での脱プロトン化をシミュレーションする
ための CDFT 分子動力学シミュレーションを応用し、自由エネルギー摂動法により高精度な自由エネルギー曲面を計算する手法を構築した。現在はこの手法を改良し通常の DFT 計算と同程度の計算コストで自由エネルギー計算する手法を開発
し、アモルファスシリカ・水界面に応用し界面構造を明らかにする予定である。
近年では電極界面での界面分光測定も可能になりつつあるため、金属や半導体界面でのシミュレ
ーションを行う。銅の CVD 法による堆積のシミュレーションを行い実験的に観測された律速段階を明らかにすることができた。[2-3]今後は、酸化チタン・水界面での水分子の役割や異なる界面における反応性の違いなどを解析する予定である。
以上で計算した界面について、現在界面分光ス
ペクトルの計算法の構築をしている。申請者の先行研究[4-6]をさらに拡張し、分光スペクトルからわかる界面構造と光触媒反応メカニズムを議論する予定である。
謝辞
本研究では、京都大学学術情報メディアセン
ターの共同研究制度（若手奨励枠）を活用させていただきました。
Figure 1. (Top) Aqueous silicic acid [1] (Bottom)

参考文献[1] T. Joutsuka and K. Ando, J. Phys. Chem.

Photocatalysis at anatase TiO2(112).

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 11

水平基板上を濡れ拡がる液膜と単一球形粒子との相互作用による接触界線近傍流体の加速現象に関する数値解析中村 颯
東京理科大学理工学研究科機械工学専攻

1 緒言

2 計算手法濡れとは一般的に液体が固体表面上を拡がる現本研究では，オープンソース数値流体解析ソフ象である. 濡れは自然界に偏在し, 自然は濡れをトウェアの OpenFOAM (Ver.5.0)を用いて数値解有効に活用している.濡れの特性は工業製品にも析を行った.気液二相流を扱う計算手法として数多く応用されている. 特に, 熱の移動効率を高Volume of fluid (VOF) 法(3)を用いた.本研究で扱めるヒートパイプや, lab on a chip デバイスなど,う支配方程式は以下の三式となる.
(1)

液体の効率輸送を目的とした工業機器の向上のためには, 固体表面上の液体の濡れ拡がりにおける(2)

物理的な理解とその効率的な制御が必要不可欠である. 最近, 水平基板上に静置された球状粒子とM-CL 近傍流体との相互作用の微視的な動力学が(3)示された. Mu et al.(1,2)は, 水平基板上を濡れ拡が式 (1)-(3) はそれぞれ, 非圧縮性流体の連続の式,

る液膜と基板上に静置された単一球状粒子との相ナビエ・ストークス方程式, 含有率αの移流方程互作用が, M-CL 近傍流体の急激な加速現象を誘式である.αは一つのセルにおける液相の体積分起することを示した.また, 連続的に微小構造体率を表す. 本研究では Continuum surface forceを配列することで, 液体の輸送方向の制御が可能(CSF) 法(4)により界面張力を FSV = σκnn =σκになることを示した.これらの加速現象は，球状粒∇αと表す. ここで，σ,κ, n はそれぞれ界面張力

子周りに形成される曲率を持った液体（メニスカ

係数, 界面曲率, 界面の気体方向の法線ベクトルス）に起因すると考えられている．実験では球形を表す.

粒子及び円柱構造物と液膜の相互作用時のメニス

本数値解析での解析対象を Fig.1 に示す．計算カス形成過程を詳細に観察し, 構造物周りに形成領域の寸法と格子数はそれぞれ (Lx, Ly, Lz) =されるメニスカス形状が加速現象に及ぼす影響を(300+L, 80, 150) µm, (Nx, Ny, Nz) = (Lx / 2, Ly / 2,示している.しかし、微小系のため液体の対流場のLz / 2) とした. ただし L は 2 粒子の中心間距離で観察は困難であり、加えて加速現象に重要とされあり, 本数値解析では 3 ケース (L = 75, 100,る圧力場の議論は欠けている．障害物形状の違い

125)の計算を行った.基板表面は, no-slip 条件を課が液膜に及ぼす影響を十分に理解するためには,している. 上面(x, y = Ly, z)は開放系,側面 ((x, y,

数値計算を用いた解析が必須である.

z = 0); (x, y, z = Lz))は対称系である. 2 つの粒子半本研究では, 数値解析を用いて粒子周りに形成

径は Rp = 25 µm であり，粒子，基板の接触角はされる圧力場に着目し, 動力学的な観点から,連をそれぞれθs = 5°，θp = 20°とした．液体は 2

続的に配列された微小構造物と液膜間の相互作用

cSt シリコーンオイル，気体は空気の物性を模擬による M-CL 近傍流体加速現象を描像する.本数して用いる. 液体の表面張力係数, 密度, 動粘度

値解析では粒子間距離をパラメータとし, その影はそれぞれσ = 1.83 × 10−2 N/m, ρl = 873 kg/m3,

響を考察した.

νl = 2.0 × 10−6 m2/s である. 加えて, 気体の密度,
動粘度はそれぞれρg = 1 kg/m3, νg = 1.48 × 10−5
m2/s である．重力加速度は g = 9.81 m/s2 である.

12 − Academic Center for Computing and Media Studies, Kyoto University

3 結果と考察

Computational Physics, 201 (1981) 39-1.
[4] J. U. Brackbill, D. B. Kothe and C. Zemach,

数値解析によって同現象の再現を試みた．そのJournal of Computational Physics, 335 (1992) 100-2.

結果、粒子周りにメニスカスが形成され、実験を忠実に再現することができた（Fig.2 挿入図）．

謝辞：本研究を進めるにあたり，京都大学学術情

M-CL 速度 UCL は一つ目の粒子と接触後すぐに

報メディアセンターの方々には多大なご協力を賜

増加し、最大速度 max(U1)に達した後減少した

りました．ここに記し，感謝の意を表します．

（Fig. 2）
．二つ目の粒子と接触後も M-CL は同様の加速及び減速を示した．異なる粒子間距離で比較をすると、粒子間距離 L が大きい程 2 つ目の粒子と接触後の最大速度 max(U2)が小さくなった．
粒子周りの圧力分布に着目すると、液膜が粒子と接触後、メニスカス内部は低圧になった（Fig. 3（a））．また、粒子の上流及び下流側のメニスカ
ス内部に圧力差Δp が生じた．圧力差が最大になる時刻とM-CLが加速する時刻は殆ど一致しており、粒子周りの圧力差が液膜の駆動源であることを明らかにした（Fig. 3（b））
．加えて、粒子間距

Fig. 1 Computational domain of spreading

離 L が大きい程 2 つ目の粒子周りの最大圧力差

liquid film on horizontal substrate with twomax(Δp2)が小さくなることも本研究より明らか

particles.

になった．

4 結言
以上より、粒子間距離の違いがメニスカス内部の圧力差に影響を及ぼし、その結果液膜挙動に違いが生じることを数値解析によって明らかにした．
本研究では、簡易的に静置された基板上の粒子がFig. 2 Temporal variations of the CL velocity

省エネルギで液膜を駆動、すなわちポンプの役割

UCL on substrates in the case of L = 75, 100, 125.

を果たしており、粒子間距離の違いによってポン

The CL position is derived on the line (x, y = 0, z

プの性能が変化することを示している．本研究か= zp). Note that the CL position is not wellら得られた知見は、高効率液体輸送デバイスにおdefined in the gray area, in which the CL is

ける、構造の最適化に密接に関わってくる．

traveling around the particle foot.

5 参考文献[1] L. Mu, D. Kondo, M. Inoue, T. Kaneko, H. N.
Yoshikawa, F. Zoueshtiagh and I. Ueno:
Journal of Fluid Mech., 830 (2017) R1.
[2] L. Mu, H. N. Yoshikawa, D. Kondo, T.

Fig. 3 (a) Temporal variation of pressure

Ogawa, M. Kiriki, F. Zoueshtiagh, M. Motosuke,around particles. (b) Temporal variations of

T. Kaneko and I. Ueno: Colloids and Surf. A,

pressure difference between upstream and555 (2018) 615–620

downstream sides in meniscus, and CL velocity[3] C. W. Hirt and B. D. Nichols, Journal of

UCL.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 13

Ｎ結合型糖鎖修飾によるタンパク質の機能制御の関連性
リントゥルオト正美
京都府立大学大学院生命環境科学研究科応用生命科学専攻

1 緒言

影響を受けることがわかった。
これまでの計算では複雑なエネルギー曲面上の細胞外のタンパク質や膜たんぱく質の多くが糖比較的似通った局所的な構造のサンプリングのみ

鎖修飾を受け、タンパク質の構造や機能に大きな影を行っている可能性はぬぐえない。今回はさらにサ響を与えている。糖鎖修飾は新規医薬品や医療材料

ンプリング空間を広げた計算を行うことを目的ととして注目され研究が多くなされてきたが、糖鎖はして、拡張アンサンブル法の一つであるレプリカ交

そのフレキシビリティーの高さから固体結晶構造換法を用いた計算（REMD）を用い、N 結合型糖鎖は得られておらず、機能との関連性は完全に解き明修飾のタンパク質の構造に与える影響と反応制御

かされてはいない。

の関連性について調べた。

多様な組織で発現し様々な疾病との関連性が示
唆されている kallikrein（KLK）ファミリーは異な

2 方法

る基質特異性を有するセリンプロテアーゼでこれまでに 15 のタンパク質が報告されている。KLK フ

N結合型糖鎖修飾モデルとして、
ヒトKLK8(PDB

ァミリーの反応サイトに共通して存在する Loop 99

ID:1npm)の Loop 99 内 Asn95 にアセチルグルコサ

は酵素によってその長さは異なるが糖鎖結合部位

ミン(GlcNAc)、マンノース(Man)、ガラクトースを有しセリンプロテアーゼの基質特異性に関係し(Gal)から成る糖鎖を結合したモデル、糖鎖なしのていると考えられている 1)。

WT モデルを用いた。これらのモデル周辺 20 Å に

KLK ファミリーに属する KLK8（neuropsin）は水を、生理条件になるように Na+と Cl-を配置した。

中枢神経系で発現しており、このタンパク質の機能

MD 計算には Gromacs3)を用い、アミノ酸残基に障害により、記憶障害や精神疾患との関連性が示唆

は amber ff14SB 力場 4)、
糖鎖には GLYCAM06-j5)、

されてきた。本研究に先立って KLK8 とその基質で水分子には TIP3P6)を用いた。最急下法によるエネある neuregulin-1 の複合体のモデルを用いた分子

ルギー最小化に続いてモデルに位置拘束を課した

動力学計算を行った。その結果、KLK8 の基質特異状態で 300 K まで昇温した。次に 1 ns かけて徐々性は基質のアミノ酸配列に左右され、基質が Loopに拘束を解除しながら、NPT 計算をおこなった。

99 と効率的に相互作用を行うことが重要であるこ

300 K から 450 K までの間の 70 個のレプリカを用とがわかった 2)。さらに Loop 99 の N 結合型糖鎖修意しそれぞれ 2 ns の平衡化を行い平衡構造とした。

飾と機能制御の関連性について分子動力学法（MD）

この際に Loop99 および糖鎖を除くたんぱく質の位法を用いて研究を進めた。糖鎖修飾がない場合に比置拘束を行い、Loop99 および糖鎖のみの構造サン

べ糖鎖修飾したモデルでは結合部位周辺のループ

プリングを行った。70 個のレプリカを用いて、の揺らぎに変化が表れ、特に Loop 99 の揺らぎが変

REMD で 100 ns の計算を行った。

化した。また、糖鎖修飾によって結合部位の構造が14 − Academic Center for Computing and Media Studies, Kyoto University

Fig.１ Principle component analysis for 100 ns trajectory of N-linked glycosylated model. The sticks
represent the sugar chain.

3 結果と考察なると考える。
これらの構造では Loop99 と糖鎖以外のタンパク

100 ns の trajectory を用いて主成分分析を行った

質構造の位置拘束を行っていたため、糖鎖の揺らぎ

結果、N 結合型糖鎖モデルにおいて糖鎖は（Fig. 1による基質ポケットの構造変化についての情報が右上）主に初期構造の周り、すなわち主に溶媒中を得られない。これらの糖鎖の振る舞いによる基質ポ

揺らいでいた。しかし、Fig. 1、２に示すように糖鎖ケットの構造変化に関するシミュレーション及びがタンパク質表面上に沿う形の構造が現れており、解析をさらに行うことによって触媒活性との関連

これらの構造はレプリカアンサンブル法を使う前性が明らかになるであろう。また、シアル酸で延長では見られなかった。また、タンパク質表面上に静した糖鎖は他タンパク質においてタンパク質機能

電ポテンシャルをマッピングすると、糖鎖は正電荷との関連性も指摘されていることから検討が必要を有するタンパク質表面上に沿うように存在してであると考えている。

いた。

4 参照文献1)

Skala, Wolfgang et al. J. Biol. Chem. 289,

34267, (2014)2)

リントゥルオト正美他、日本コンピュータ化学会 2017 秋季年会精選論文特集号 16 巻（5号）、160-162
Fig. 2

The fluctuation of sugar chain on theprotein surface of the glycosylated model.

Thesticks represent the sugar chain.

3)

S. Pronk, et al., Bioinformatics, 29, 845,

(2013).
4)

J. A. Maier, et al., J. Chem. Theory

Comput., 11, 3696, (2015).
基質ポケット内は負電荷を有し、基質の切断点である Arg が接近しやすくなっている。一方で基質ポケット周囲は正電荷を有しておりこれをふさぐ形
で糖鎖が存在していることから、糖鎖のこのような

5)

K. N. Kirschner, et al., J. Comput. Chem.,

29, 622, (2008).
6)

W. L. Jorgensen, et al., J. Chem. Phys.,

79, 926, (1983).

conformation は基質の接近に影響を与えることに京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 15

Numerical simulation of deepwater oil blowout
Crossflow and droplet size distribution
Daniel Cardoso Cordeiro

大阪大学大学院基礎工学研究科物質創成専攻化学工学領域

1 Introduction

2 Numerical methods

In 2010, one of the largest offshore oil spills of

To investigate the relationship between oil

all time happened in the Gulf of Mexico, USA.

plumes and the aforementioned dimensionless

For the first time, the sub-sea injection of

numbers, numerical simulations of oil blowout

chemical dispersants was used to treat

in a cross flow were executed using the

deepwater oil spills. However, with only few

schematics shown in Fig. 1.

studies prior to its application, the overall

The governing equations solved are theeffectiveness of this method is still questioned

Continuity (Eq. 1) and Navier-Stokes (Eq. 2)

as appropriate measures of the oil droplets

equations:were not performed in situ.

(1)

In order to investigate the impact of thewater current crossflow in an oil blowout, as
real-scale experiments are clearly unfeasible,(2)

with the aid of numerical simulations it ispossible to get some information that might
help

us

improve

ourunderstandingof

deepwater blowouts in a cross flow. The jet
trajectories can be analyzed [1]. Specially, byusing dimensionless numbers, it is possible toassess the effect of several parameters in the oil
jet behavior. In this research, three parametersthat influence the jet under crossflow areconsidered: the Weber number (Wew = ρwU2D/σ),
the Reynolds number (Rew = ρwUDμ), and theinlet velocity ratio (r = Uo/Uw) [2], where ρ is thedensity, D is the oil inlet diameter, U is the inlet
velocity, μ is the water viscosity, and thesubscripts w and o indicates water and oil,respectively. We investigate the effect of theseparameters on the behavior of the oil jet andthe velocity field.

16 − Academic Center for Computing and Media Studies, Kyoto University

where u isthevelocity, isthevolumefraction, t is time,  is the density, p is thepressure, g is the gravity acceleration, FS,k isthe surface tension force, FD,k is the drag force,

Fvm,k is the virtual mass force and thesubscript k indicatesthe fluid

phase. The

Euler-Euler method was used to model themultiphaseflow

Schiller-Naumann

coupleddragwithmodel.

aThe

open-source OpenFOAM v4.0 was used for the
calculations.

3 Results and Discussion
Averaging the velocity and volume fraction
for 30 seconds (t = 10 s – t = 40s ). In Fig. 2, theeffect of varying the velocity ratio (r) on the oil

areDetailed simulations of primary breakup of

non-dimensionalized by the oil inlet diameter,

turbulent liquid jets in crossflow. CTR Annual

D. With the increase of Rewater, Wewater, and r,

Research Briefs.

plume

trajectory.

The

distances

the trajectories of the oil plumes centerlines
were ascendant, while with lower r the plumes
quickly became horizontal. These effects arehighly dependent on the value of r, which
seems to change the behavior of the plume in
the interval 1 ≤ r ≤ 2.
In Fig. 3, the oil mean velocity profiles for r =1 were shown. As the oil plume developed, theoil plume was caught in the effect of thevortexes near the inlet, which did not happen sostrongly in the values with r above 2. Thismechanism

could

explainthestrongdependence of the plume centerline tendency

Figure 1: Schematics of oil jet in water crossflow
including boundary conditions.

according to the r ratio.

4 Conclusion
The effects of the oil jet in a crossflow were
investigated using an Euler-Euler multiphase
numerical method.
The trajectory of the oil plume centerline wasstrongly affected by the ratio r between the oil
inlet velocity and the water inlet velocity.
Vortexes near the inlet could be a mechanism
in favoring the faster horizontal trend of theplume trajectory of the instead of a fasterFigure 2: Oil jet centerline trajectory for varying

Rewater, Wewater, and r

plume ascendance.
Techniques that could generate such a vortex
near the well inlet could help in extending theresidence time of oil in the deeper ocean which
could enhance biodegradation.
Further research is needed to investigate thejet in the crossflow phenomenon.

References[1] New, T. H. T. Lim, T & C. Luo, S. (2006).
Effects of jet velocity profiles on a round jet in

Figure 3: Oil mean velocity field and streamline forthe case with r = 1.

cross-flow. Experiments in Fluids.
[2] M.G. Pai, O. Desjardins, H. Pitsch. (2008).

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 17ຌ੄ʍφഷॐߢ‫ػ‬໑˙ƪˑʊ੆ɸʪ
Ցҗӌࡌʊʧʪ˴˙˽ʍ‫ۥ‬२઺π٩‫∗ ڭ‬∗୼֮Ҙ๎੝ӌ ӌࡱٰ֖ϙ ຌ૾࣮൙‫ۑ‬ӌ೼ฉච‫ۮ‬ʎঀ෼ ոຒ߅ ( φׄ੝ӌ‫ؼ‬ϼԛ๽ٰ֖Ѡ ) ʇʍִடٰ֖ [5] ʊՂʄɮʡʍʆɡʪ.

1

ʎɷʠʊ.
Ցҗӌࡌʎഐ๽ӌʍํƧʉഒฐʆુ෾ɴʫʅɣʪ.

‫ׯ‬௻, Ցҗӌࡌʍφ࠱ʆ, ˼ʽƪ˾̅˚˟˕˚́ƪˁʱʡʂɾ˼ˈƪˢƪ˅̅˦˷ƪ˘ʵ̅˂ [3, 6] ɫߢ‫ػ‬໑ʣ˼˵˩ˠ˧޼ॐʉʈʍ฽ਢʊɩɣʅอ‫ۄ‬ʆɡʪɲʇɫ൙‫ܘ‬ɴʫʅɣʪ. ѼƧʎɲʍՑҗӌࡌࠬ൥ʱ๑ɣʅຌ੄ʍʺ˟˽ˀƪഷॐʍ˴˙˽ʱ‫ۥ‬२ɶ, ฽৳ʱ२‫ۃ‬ɴɺʅɣʪ [4].
φൣʆ, ຌ੄ʱ޳ʠʇɶɾ‫ٿߣ܊‬ʍອӌ‫ػ‬ʍߢ‫ػ‬໑˙ƪˑʱॴ୩ວɮ˴˙˼̅˂ɸʪɾʠʊʎ‫ۄ‬ຆʧɮӌࡌɸʪɲʇɫ֑ʠʨʫʪ. ɼɲʆ, ߢԨભʫ݈ೀ
‫ػ‬ʱެɥɲʇʆ฽৳ʍॴ୩ʱࣣɱɾ. ʝɾ, ວɣߢԨભʫ݈ೀ‫ػ‬ʍʇʩൣʱ๑ɣʪɲʇʆຌ੄ʍˋʽ˻ƪ˰ˁ˿ഷॐʍߢԨౙୈ˴˙˽ʍ‫ۥ‬२ʱɸʪ.

2˼ˈƪˢƪ˅̅˦˷ƪ˘ʵ̅˂.

dφ/dt = f (φ) ʆ೅ɴʫʪອӌ‫ػ‬ʍഷॐ, u = h1 (φ)∈ RM ʇ s = h2 (φ) ∈ RP ʊʃɣʅ۵ɧʪ. ɡʪߢ‫ܗ‬ʝʆ (ӌࡌߢԨʇ‫ڐ‬ʕ ) ʍ௬ອഷॐ u, ࡰອഷॐ sʍߢ‫ػ‬໑ʎՍઢʇɸʪ. s ʍߢ‫ػ‬໑ʱ฽ਢɸʪߢ‫ܗ‬ʊɩɣʅ, ߢ‫ػ‬໑˙ƪˑ u ʎෆઢʇɸʪ. ഐ๽଺ʉઢ‫ٵ‬ʱ๑ɣɹʊӌࡌߢԨʍ˙ƪˑˍ˕˚ {u, s} ʍʞɪʨ˴˙˽ʱ‫ۥ‬યɸʪ.
˼ˈƪˢƪ˅̅˦˷ƪ˘ʵ̅˂ʍற૙ʍφʃʎ˝˷ƪ˻˽˟˕˚́ƪˁʍ‫ۥ‬ਚ߭੄ʱӌࡌɶʉɣɾʠ
Ώ୭଺ʊ‫ޟك‬ສɫࢭʉɮʅ‫ݗ‬ʞ, ਜ਼ʮʩʊ˝˷ƪ˻˽˟˕˚́ƪˁʍߣ‫ٿ‬ʱ੝ɬɮɶʅɣʪ୐ʆɡʪ (ࣚ‫ݟ‬ʎ [5]). ɲʍɲʇʊʧʩ˒ʶ˜˱ˁˋɫഉ‫ތ‬ʉʔ

18 − Academic Center for Computing and Media Studies, Kyoto Universityʪʝɣʱߪɸʇɶʅʡ, ɼʫɫٔଜ໼଺ʆɡʪ࣪܏
ʊʎɲʍ࠱ʍӌࡌൣ൥ɫอ‫ۄ‬ʆɡʪ.

3

˴˙˽ʍঙଜʇٗѢ.

3.1

௬ࡰອഷॐʍঙଜ.

˘ƪ˻ƪ˰ʶˁ˿ˋ˃ƪ˽˾ʶˠ˽ˌॐ. ࡀՎֻҚࣰٛʍђʆ 3 ߣ‫ٿ‬ಝΏ࡬ Navier–Stokes ൣପ߲ʍ૰খॐડ‫[ ޟك‬2] ʊʧʂʅமʨʫɾߢ‫ػ‬໑˙ƪˑʱ˼ˈƪˢƪ˅̅˦˷ƪ˘ʵ̅˂ʍӌࡌ˙ƪˑʇɶʅ
๑ɣʪ. ɲɲʆʎ, Һߢ‫ܗ‬ʆʍຌʫʍഉ‫ތ‬ɴʱ೅ɸ˘ƪ˻ƪ˰ʶˁ˿ˋ˃ƪ˽˾ʶˠ˽ˌॐ R̃λ (t) [1] ʊુ෾ɸʪ.
ߢԨભʫ݈ೀ‫ػ‬. M ߣ‫ٿ‬ʍ ∆τ ߢԨભʫ݈ೀ‫ػ‬ʱ௬ອഷॐ u(t) = (u1 (t), u2 (t), · · · , uM (t))T , ࡰອഷॐ s(t) = (s1 (t), s2 (t), · · · , sM (t))T ʇɶʅ஡௬ɸʪ.
ɾɿɶ, ui (t) = si (t) = R̃λ (t − (i − 1)∆τ ) ʇɸʪ.

3.2

ߢԨભʫ݈ೀʇભʫߢԨʍٔʠൣ.

߭‫ਂڔ‬ԪԪॐ. ˾ʶˠ˽ˌॐ {Rλ (t)} ʍ߭‫ਂڔ‬ԪԪॐʱ C(x), ɼʍൗ๭সʱ Ce (x) ʇɸʪ. ऺ 1
ʊ, C(x) ʇ Ce (x) ʱ࢑ɬɿɶɾ. ऺ 1 (ܾऺ) ʧʩ,C(x) = 0.7 ɪʨ 0.3 ʍડʱʇʪʍʎ, x ≈ 3.0 ɪʨ
5.0 ʍߢʆɡʪ. ɲʍɲʇɪʨ ∆τ < 3.0, ʉʨʏ,Rλ (t) ʇ Rλ (t − ∆τ ) ɫߛɾஞɬʱɸʪɲʇʊʉʩ,φൣʆ, ∆τ > 5.0 ʉʨʏ, ຜࠖʍਂԪɫࢬɴɮʉʪ.
ɲʍɲʇɪʨ, ߢԨભʫ݈ೀ‫ػ‬ʍભʫߢԨ ∆τ ʇɶ
ʅ 3.0-5.0 ɪʨূʕʇວɣʇीਢɴʫʪ.
ࠄ‫ݣ‬, ∆τ = 3.0-5.0 ʊʃɣʅ, ‫଼ݍ‬ʉભʫߢԨ∆τ ʇ ભʫߣ‫ ٿ‬M ʍৠʞ܏ʮɺʊʃɣʅ૦݀ɶɾ(೅ 1). Һ (∆τ, M ) ʊʃɣʅ, ҺƧ8160 ˣˑƪ̅ʍ˝˷ƪ˻˽˟˕˚́ƪˁʊʃɣʅӌࡌʱ‫ۼ‬ɣ, ฽৳

1

0

s1

C(x)

C(x)

0.5
0.3

0

0

2

4

x

6

8

10

-1

2

reservoir output
actual

-Ce
0

20

40

x

60

80

100

reservoir output
actual

1.5

1

1

0.5

0.5

0
-0.5

-0.3
-0.5

-0.5
-1

2
1.5

Ce

s1

1
0.7
0.5
0.3

0
-0.5

-1

-1

-1.5

-1.5

-2
1000

1050

1100

1150

1200

1250

-2
5000

t'

5050

5100

5150

5200

5250

t'

ऺ 1: ߢԨભʫ x ʍߢԨʊ੆ɸʪ Řλ (t) ʍ߭‫ਂڔ‬Ԫऺ 2: ӌࡌʊʧʩமɾ˴˙˽ʍφʃʱ๑ɣʅํƧʉԪॐ C(x) ʇൗ๭স Ce (x) (Ϟऺ) ʇɼʍҼ੝ऺ (ܾߢԨ Tout < t′ < Tout + 250 (Tout = 1000, 5000) ʆऺ).

ʍߢ‫ػ‬໑˙ƪˑ฽ਢʱ‫ۼ‬ʂɾ.

ɶɾഷॐ ŝ1 (t) ʇ s1 (t) ʍߢ‫ ܗ‬t ∈ [0, 60](t ∈ [0, 90])ߢԨભʫ݈ೀ‫ػ‬ʍભʫߢԨ ∆τ ʇ ભʫߣ‫ ٿ‬M ʊʃʍܿʍ‫ݍ‬੝ડɫ e60 = 0.14(e90 = 0.30) Τђʊʉʪɣʅ, ߭‫ਂڔ‬Ԫ‫ؤ‬ॐ C(∆τ ) ɫ 0.45-0.55 ɪʃ, ൗ๭স Ce (M ∆τ ) ɫ 0.35-0.40 ʱෂɾɸʧɥʊ (∆τ, M )ʡʍʍ‫ڎ‬ॐʱॐɧɾ.
∆τ \ M
3.0
3.5
4.0
4.5
5.0

(e60 , e90 ) = (0.14, 0.30)
11
12
13
14
0
0
0
0
0
0
11
20
3
18
43
107
14
43
54
21
24
26
19
9

ʱٔʠʪʇ‫଼ݍ‬ʆɡʪɲʇɫʮɪʂɾ. ʝɾ, ɲʍʧ
15
1
28
59
15
1

16
19
57
21
8
1

17
24
48
14
1
1

18
43
21
2
1
0

೅ 1: ߢԨભʫ݈ೀ‫ػ‬ʍҺભʫߢԨ ∆τ ʇߣ‫ ٿ‬M ʆʍ฽ਢʍ२‫҉ۃ‬ॐɥʊɶʅٔʠɾભʫߢԨ݈ೀ‫ػ‬ʱ‫ٿ‬ʊӌࡌɶɾ˴˙˽ʎ, ߢ‫ػ‬໑˙ƪˑʍ฽৳ɫʆɬ, ਤ୩࣪ɫѐʣɪʉຌʫɪʨഉ‫ތ‬ʉຌʫʊߢԨԨٓ଺ʊৃζɸʪํ޶ʡ
‫ڇ݌‬ɸʪɲʇɫӂ௳ʆɬɾ.

5ࠗ߯઺πʎ, චٰ֖ʎ JSPS Ѡٰಛ JP19J12482 ʍ࢕∆τ ≈ 4.0 ɪʃ M ∆τ ≈ 55-60 ʇʉʪʧɥ (∆τ, M )२ʱ࠷ɰʅɣʪ. ʝɾ, චٰ֖ʆɩɲʉʂɾ‫ޟك‬ʍφʱଜʠʪʇ, ฽৳ʍʺ˻ƪɫɶɬɣડΤђʊʉʪ‫ڎ‬ॐ೼ʎ֮୤੝ӌʍˋƪˣƪ˅̅˦˷ƪˑִடٰ֖ॣ୩

ɫ‫ݍ‬ʡਵɮʉʪɲʇɫ‫ٵ‬ʅࠪʫʪ. ɲʫʨʍડʎɼ( ࠥࠬƔࢗ१ࢧໃ༄ ) ʊՂʄɮ. ɲɲʊԈࠗʍίʱ೅ɸ.

ʫɽʫ, C(∆τ ) ≈ 0.45-0.55 ʇൗ๭স Ce (M ∆τ ) ≈0.35-0.40 ʊਂஆɸʪ (ऺ 1 (Ϟऺ)).

3.3‫ޖ‬۵ഞ‫ٯ‬ӌࡌ˴˙˽ʊʧʪߢ‫ػ‬໑˙ƪˑ฽ਢ.

(∆τ, M ) = (4.0, 14) ʍʡʇʆߢԨભʫ݈ೀ‫ػ‬ʱ۵ɧ, ˼ˈƪˢƪ˅̅˦˷ƪ˘ʵ̅˂ʊʧʩߢԨౙୈ˴
˙˽ʱ 1 ʃமɾ. ɲʍ˴˙˽ʱ๑ɣʅํƧʉࢊɪʨʍ
ߢԨౙୈ฽ਢʱ‫ۼ‬ʂɾ. ˾ʶˠ˽ˌॐ s1 = R̃λ ʍߢ[1] T. Ishihara and Y. Kaneda, In Statistical Theories and Computational Approaches to Turbulence(2003), pp. 177–188.
[2] K. Ishioka, GFD Dennou Club, (1999), http://
www.gfd-dennou.org/arch/ispack/.

‫ػ‬໑˙ƪˑʍ฽৳ɫɥʝɮ‫ۼ‬ʂʅɣʪɲʇʱӂ௳ɶɾ

[3] H. Jaeger, and H. Haas, Scince, 304, (2004), pp.
78-80.

(ऺ 2). ऺʆʎ, ಐӇʍɾʠ Navier–Stokes ൣପ߲ʍ૰খॐડ‫ޟك‬ɪʨமʨʫɾ९҈ʍߢ‫ػ‬໑ s1 = R̃λ (t)[4] K. Nakai, and Y. Saiki, Physical Review E 98,(2018), 023111:1-6.

ʇɡʮɺʅ࢑ɬɿɶɾ.

4[5] K. Nakai, and Y. Saiki, Discrete and Continuous Dynamical Systems Series S, (2020), doi:
10.3934/dcdss.2020352.

ʝʇʠຌ੄ʍ˰ˁ˿ഷॐʇɶʅφഷॐߢ‫ػ‬໑˙ƪˑʆɡʪ˾ʶˠ˽ˌॐʍ֣ஞʱ೅‫ڇ‬ɸʪ˼ˈƪˢƪ˴˙˽[6] Z. Lu, J. Pathak, B. Hunt, M. Girvan, R. Brockett,and E. Ott, Chaos, 27, (2017), pp. 041102. pp.
391-403.

ʱߢԨભʫ݈ೀ‫ػ‬ʱެɥɲʇʆ‫ۥ‬२ɶɾ. ʝɾ, ɲʍ

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 19

Numerical simulation of the transport phenomena occurring in the
Bottom Seeded Solution Growth (BSSG) process of SiC
ZHANG Jiangao
Graduate School of Engineering Science, Osaka University

1 Introduction

frequency coil on the crystal growth rate and its uniformity
will be comprehensively considered.

Silicon carbide (SiC) is a promising wide-bandgap
semiconductor material for power devices due to its high

2 Numerical method

electrical breakdown field and high thermal conductivity. The
The BSSG furnace configuration is schematicallyshown in Fig. 1. The radio coil is set outside the insulator, and--order to avoid these disadvantages, solution growth processes
have been developed to produce high-quality crystals sincethese processes are performed at the condition close tothermodynamic equilibrium. One example of such solution
growth processes is the Top-Seeded Solution Growth(TSSG) method. However, it was found that, in this growth
process, the Marangoni convection developing along the melt

the Si-C melt is heated by the crucible, which is heated by the
induction coil. The simulation was based on the integrated
process model [5, 6] and carried out separately for the threestages: (i) the electromagnetic field, (ii) the temperature field,and (iii) the fluid flow in the melt. We first get the Lorentz
force, and the heat generation densities from step (i). After
that, we compute the temperature field (step (ii)) by using thesurface is strong and affects the crystal homogeneity. To

heat generation densities from step (i). In the end, we carry out

overcome this issue, we numerically study the Bottom

the melt flow in the step (iii) by using the heat generation

Seeded Solution Growth (BSSG) method for the growth ofdensities and Lorentz force from step (i) and the boundary

SiC, since the Marangoni convection is not significant

condition of temperature from the step (ii).

because the seed is at the bottom of the crucible. Also, thistechnique had shown a reduction in dislocation density in thegrowth of GaAs [1].
Most numerical simulations were conducted in the 2D
axisymmetric model [2, 3]. In this study, we carried out a 3D
numerical simulation to investigate the transport phenomena
occurring during the BSSG process of SiC. The coilfrequencies have a significant effect on electromagnetic andbuoyancy flows, and the Marangoni, electromagnetic, andbuoyancy forces are the essential parameters of the growth of
SiC crystal. Yamamoto et al. [4] examined, in the TSSG
process, the effect of coil frequency on flow patterns in adifferent furnace configurator, but the crystal growth rate andits uniformity were not considered. In the present BSSG
method, the effect of the magnetic field induced by the radio

20 − Academic Center for Computing and Media Studies, Kyoto University

Figure 1 Schematic view of the growth setup

3 Results and discussion
The computed maximum values of the Lorentz force
density in the melt and the Joule heat density in the wholefurnace are plotted in Fig. 2 and Fig. 3. The maximum Jouletemperature value when the coil frequency is lower than 30

heat density increases with the coil frequency. On the other

kHz. It happens because the stronger Lorentz force will behand, the maximum Lorentz force in the melt increases

generated in the melt at the low frequencies, and thequickly at the low coil frequency. When the frequency

electromagnetic flow is dominant. This developed flow thinsexceeds a certain value (5kHz), the Lorentz force decreases

the concentration boundary layer along the crystal growthwith the increase of the coil frequency.

interface and consequently makes the growth rate higher. The
flow structures in the melt are almost identical, and this9

indicates that the frequency has no significant effect on the8

uniformity of crystal growth. When the coil frequency ishigher than 30 kHz, the growth rates almost keep the same

7
6
6

6

because Marangoni flow is dominant.

2
2

6

2

800

5

1kHz
10kHz
15kHz
30kHz
80kHz

600

Figure 2 The maximum value of Joule heat density in the
furnace at various frequencies and current densities.

0
0
0

6

200
0

2
6

6

400

1600
2

2

1800

2000

2200

2400

2600

Max temperature in the melt (K)2800

Figure 4 The time-averaged growth rate at various
requencies

4 Conclusion
The influence of coil frequencies on the crystal
Figure 3 The maximum value of Lorentz force density in themelt at various frequencies and current densities.

growth rate and its uniformity during the BSSG
process of SiC is numerically investigated. The result showsthat the coil frequencies strongly affect the growth rate andWe can also find that the maximum Lorentz force

have no significant effect on the uniformity of crystal growth.

value is located at the center edges of the melt for low

The growth rate can be enhanced at low frequency because of

frequencies while it is located at the top and bottom corner ofthe more potent Lorentz force. For low-frequency coil, thethe melt for high coil frequencies due to the skin effect, whichelectromagnetic flow is dominant. However, Marangoni flow

may affect the flow structure in the melt.

is dominant when the coil frequency exceeds a specificFigure 4 shows the effect of the coil frequencies on thegrowth rate. For all the cases investigated, we only enhance
the current intensity to increase the temperature of the melt.
As seen, the growth rate increases with the increase of thetemperature. More carbon will dissolve from the cruciblewalls, and the concentration of carbon in the Si-melt will behigh at a higher temperature. As a result, carbon flux at thegrowing interface will increase.
Besides, we can also find that the growth rate will
increase with the decrease of the coil frequency at the same

frequency.

5 Reference[1] K. Hoshikawa, et al., J. Crys. Growth, 94 (1989) 643.
[2] K. Ariyawong, et al., Mater. Sci. Forum, 71 (2014) 778.
[3] T. Umezaki, et al., Mater. Sci. Forum, 31 (2015) 821.
[4] T. Yamamoto, et al., J. Cryst. Growth, 474 (2017) 50.
[5] P.M. Gresho, et al., J. Cryst. Growth, 85 (1987) 40.
[6] J.J. Derby, et al., J. Cryst. Growth, 97 (1989) 792.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 21

高プラントル数液柱温度差マランゴニ対流における自由表面を介した熱損失の影響に関する数値解析
藤本 空知
東京理科大学大学院理工学研究科機械工学専攻

1 緒言

目的とする。また、今後行われる予定の流体実験である JEREMI(Japanese European Research

シリコン単結晶生成の一つにフローティング
ゾーン(FZ)法がある。多結晶ロッドの一部をリンExperiment on Marangoni Instabilities)に向けて知見を蓄積することも目的である。

グ状のヒーターで溶解させる方法で、坩堝など他∗の物質に接触することがないため高品質な結晶

T = THを生成できる。溶解部は表面張力により保持されるため地上では小径でなければならないが、軌道上、つまり微小重力下では大径化が可能である。

H

さて、溶解部を保持する表面張力だが、温度差に∗z∗

T = TC

より張力差が生じ、溶解部内にはマランゴニ対流が生じる。微小重力化では浮力対流よりもマラン

∗

∗

r

∗

R

ゴニ対流が支配的となり、単結晶表面に欠陥縞を作ることが明らかになった[1]。以上を背景に、マθ

Fig.1 Schematic diagram of Half-Zone model

ランゴニ対流の不安定性の現象解明、制御等を目

2 計算手法的に研究が行われるようになった。
本研究では FZ 法を簡略化させ、図 1 に示すハーフゾーン液柱モデルを用いる。高温の上部ディ

図１のように液柱高さを H、半径を R とする。

スク、低温の下部ディスクの間に液体が保持され� � �/�で定義されるアスペクト比は 2.0 とした。

ており、自由表面上の表面張力差によりマランゴ
ニ対流が発生する。温度差がある閾値を超えると二次元対称流が不安定性を起こす。この不安定性に様々なパラメータが影響を与えることが明らかにされている。プラントル数（Pr � 𝜈𝜈/𝜅𝜅, 𝜈𝜈 :動粘性係数, 𝜅𝜅 : 熱拡散率）を例に挙げると、低Pr 数流体の場合は不安定後三次元定常流となり

その後三次元振動流となるが、高 Pr 流体の場合は[2][3][4]不安定後直ちに三次元振動流へ遷移する。

そして近年、液体が周囲の気体に熱を損失する
度合いも不安定性に影響を与えることが確認された[5][6]が、その原因の解明は行われていない。
本研究では熱伝達の度合いをビオ数（B� � ��/𝜆𝜆,h : 周囲気体の熱伝達率, 𝜆𝜆 : 液体の熱伝導率）で表し数値計算及び解析を行い上記原因の解明を

22 − Academic Center for Computing and Media Studies, Kyoto University

上下壁面の温度をそれぞれ𝑇𝑇� , 𝑇𝑇� とし、温度差は∆𝑇𝑇 � 𝑇𝑇� � 𝑇𝑇� と定義する。非圧縮性ニュートン

流体を仮定し、表面張力を除く物性値の温度依存性を無視した。円柱座標系�𝑟𝑟𝑟 𝑟𝑟𝑟 𝑟𝑟�を用いて、無次元化した支配方程式（連続の式、Navier-Stokes
式、エネルギー方程式）を以下に示す。
𝛻𝛻 𝛻 𝛻𝛻�⃗∗ � �

�⃗∗
𝜕𝜕𝜕𝜕
1
Bd ∗
� �𝑢𝑢
�⃗∗ ∙ 𝛻𝛻�𝑢𝑢�⃗∗ � �𝛻𝛻𝑝𝑝∗ � 𝛻𝛻 � 𝑢𝑢�⃗∗ �𝑇𝑇 ���⃗
𝑒𝑒�
𝜕𝜕𝜕𝜕
Re
Re
𝜕𝜕𝜕𝜕 ∗1� �𝑢𝑢�⃗∗ ∙ 𝛻𝛻�𝑇𝑇 ∗ �𝛻𝛻 � 𝑇𝑇 ∗
𝜕𝜕𝜕𝜕
RePr

�⃗∗ は無次元速度ベクトル、𝑝𝑝∗ は無次元圧
ここで𝑢𝑢

力、𝑇𝑇 ∗ � �𝑇𝑇 � 𝑇𝑇� �/∆𝑇𝑇は無次元温度である。以下では無次元を表す∗は省略する。Pr は 2cSt シリ
コンオイルを模擬して 28.6 とする。Re は対流の

強さを表すパラメーターであり、以下のように定

3 結果

義する。
�� �

|𝜎𝜎� |𝐻𝐻𝐻𝐻𝐻𝜌𝜌𝜌𝜌 �

LSA による中立安定曲線を図 2 に示す。今回のここで𝜎𝜎� は表面張力温度系数、𝜌𝜌は密度である。

計算パラメータでは周方向波数 m = 1 が臨界とな

Bd は動的ボンド数であり以下のように定義されった。Bi 数によって臨界値での固有値（ここではる。

Hydrothermal wave, HTW と呼ぶこととする）が変�� �𝜌𝜌𝜌𝜌𝜌𝜌𝜌𝜌
|𝜎𝜎� |

�

化することが確認できた。Bi の低いものから
HTW(b),(a’),(a)とする。固有関数が変化する Bi 域

ここで g は重力加速度、𝛽𝛽は体積膨張率である。

今回は地上実験を模擬して 0.34 とした。

では急激に対流が安定化されることが明らかになった。

上下ディスクとの速度の境界条件はノンスリ

ップ、温度はそれぞれ TH, TC で一定とした。自由表面における速度の境界条件は表面張力と剪断応力の釣り合いから導出される、いわゆるマラン
後に境界条件を採った。さらに熱の授受に関する境界条件は Bi でコントロールするよう以下に示す式で定義した。
𝜕𝜕𝜕𝜕�� ����𝐻𝐻 � 𝐻𝐻��� �𝜕𝜕𝜕𝜕 �����

𝑇𝑇��� は周囲気体温度を示す。今回は𝑇𝑇��� � � 𝐻𝐻�

としている。先行研究により Bi が正の領域、すなわち液柱が熱を損失する場合に臨界値が Bi に

Fig.2 Neutral curves as a function of Bi

敏感であることが明らかにされたため、本研究でも Bi 数を正の領域で変化させる。

固有値の変化や安定化の原因を調査するため、線形安定性解析（Linear Stability Analysis, LSA）

エネルギー解析を行った。計算手法は Wanschuraと直接数値計算（Direct Numerical Simulation, DNS）

et al.(1995)によるものを参考にした。熱的不安定により計算する。LSA 手法は Fujimoto et al.(2019)に注目するため、速度変動による基本温度場の熱と同様である[7]。

エネルギー対流輸送から生じる It1（半径方向エネ

パラメータを表 1 にまとめる。

ルギー生成）、It2（高さ方向エネルギー生成）を出力する。HTW(a’)と(a)の識別は液柱内部の温度

Table1 Numerical conditions

Prandtl number
Aspect ratio
Ambient gas temperature
Dynamic Bond number
Biot number
Reynolds number
Re
Grid size

変動分布が指標となるため、三次元構造を観察し

Pr
28.6
Γ
2.0
Tamb Tc (= 0)Bd 0.34
Bi
0.0 ~ 4.0
LSA : FreeDNS : 900 ~ 1450
�𝑁𝑁� , 𝑁𝑁� , 𝑁𝑁� � LSA :(67, 61, 120)
DNS :(50,48,100)

た。各固有関数を観察できる Bi = 0.0, 0.5, 2.0 における臨界直上の温度変動、It1、It2 の三次元分布を図 3 に示す。HTW(b)に関しては、中央の大きな構造は It1、上部ディスク付近の小さな構造は
It2 による寄与があ。HTW(a’)、(a)に関しては中央の細くつながる構造は It1、上下ディスク付近の大きな構造は It2 による寄与があった。It2 は発現した全固有関数でそれぞれ違う分布を示したため、It2 の方が強く固有値の分岐を引き起こしていると考察できる。
宇宙実験や地上実験で得られた粒子集合現象(Particle Accumulation Structure, PAS)について、重京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 23

力の有無により粒子不在領域に違いが確認された。Bd = 0.0 とした計算結果も併せて実験結果と比較することで温度変動分布との相関が確認されたが、原因解明には至らなかった。更なる解析を進め原因の調査を行うことが今後の方針であ
る。

𝑇𝑇�
It1

𝑇𝑇�

It2
Re

It1

1050

1250

1450

Fig.4 Bird-eye view of the isosurfaces of𝑇𝑇� � � � �����, It1 = +0.2 and It2 = +0.2
(Pr = 28.6, Γ = 2.0, Tamb = Tc, Bd = 0.34, DNS)

It2

4 結言

HTW

(b)

(a’)

(a)

高 Pr 数流体長液柱について地上実験を模擬し

Bi

0.0

0.5

2.0

た線形安定性解析及びエネルギー解析を行った。

Fig.3 Bird-eye view of the isosurfaces of𝑇𝑇� � � � ����, It1 = +5 and It2 = +5
(Pr = 28.6, Γ = 2.0, Tamb = Tc, Bd = 0.34, LSA)

エネルギー解析結果の三次元構造可視化から固

直接数値計算においても同様のエネルギー解PAS に関しても実験との相関を確認することが析を行った。直接数値計算は計算負荷が高いができ、宇宙実験への知見蓄積に貢献できた。DNS

Re 数を指定して対流を再現できる。Bi = 0.5 とし、によりより詳しく対流を調査でき、二次不安定現

Re = 1050, 1250, 1450 で一次不安定後定常状態と象解明の端緒を開くことができた。

有値分岐の原因調査が進み、ねじれ成分を生み出しているのは It1 であり、さらに HTW(a’)と(a)の分岐には It2 が影響していることを示した。またなった時間での温度変動、It1、It2 の三次元分布を図 4 に示す。It1,It2 のカラーリングは温度変動

5 参考文献を表す。LSA と無次元化の方法が異なるため絶対値は異なるが、構造は似たものを再現できた（図[1] A Cröll et al., J. Cryst. Growth, 191 (1998) 365.

3 の Bi = 0.5 参照）。Re = 1050, 1250 には大きな変[2] Wanschura et al., Phys. Fluids, 7 (1995) 912.

化は確認できないが、Re = 1450 まで増加させる[3] Preisser et al., J. Fluid Mech., 126 (1983) 545.

と液柱中心付近で It1 が弱まり、It2 が強まること[4] Leypoldt et al., J. Fluid Mech., 414 (2000) 285.

が分かった。また、時間変化による構造変化も確[5] B. Xun et al., Int. J. Heat Mass Transf., 54 (2011)

認されたことから、Re = 1450 では二次不安定が起こっている可能性も示唆された。エネルギー生
成の成分のうちどの項が支配的であるか、またエ
ネルギー生成の時間変化がどのように二次不安定に影響を与えているかの調査と解明を今後の課題とする。

24 − Academic Center for Computing and Media Studies, Kyoto University

1698.
[6] Wang et al., Int. J. Heat Mass Transf., 50 (2007)
4195.
[7] Fujimoto et al., J. Microg. Sci. Appl., 36 (2019)
360204.

数値シミュレーションを用いたウマの個体間に働く力の推定
井上漱太 1、深沢圭一郎 2、平田聡 1
1 : 京都大学野生動物研究センター
2 : 京都大学学術情報メディアセンター

1 はじめに互作用を、距離を変数にもつ関数によってモデル化
し、シミュレーションをおこない、野外観察によっ動物が群れをなし、動きを同調させる現象は古くて得られたデータの再現を目指した。

から研究者の興味をひいてきた。近年、計測・計算技術の急激な発展によって、一個体レベルでの行動

2 野外観察ルールの解明を目指す研究が高まりを見せている。
Couzin らは Attration、Alignment、Repulsion と野外観察におけるデータ収集は 2016 年 6 月およいう３種類の力の相互作用により、動物の群れにおび 2017 年 5-7 月にポルトガルのアルガ山でおこな

けるさまざまな振る舞いが説明できると提唱した

われた。合計 3 群を終日追跡し、ドローン(Phantom[1]。これにより、巨大な群れにおいても、近傍の個

3、Mavic Pro)により画像データを収集した。それぞ体との局所的なインタラクションが、群れの全体のれの群れの個体数は 7-8 頭だった。画像は 30 分ご維持機構であることが理論的に提示された。その後、とに高度 25-80m 程度で撮影され、一枚の画像が群

魚類の群れにおいて、個体の加速度を定量的に計測れの構成個体全てを含むように高度を調節した。ウ

することにより、この理論の有効性がすでに検証さマの群れは採食、移動、休息という 3 種類の状態にれている[2]。

大別することができ、本研究では採食場面のみを解

この分野の根底にある一つの大きな疑問は、種間析した。それぞれの群れに関して、合計 60 枚程度

で共通する個体間のインタラクションルールの存の画像を解析に使用した。個体の首の付け根と尾の在である。現状では、多くの研究が魚類や鳥類を対基部の中点に相当する点を個体の座標として、個体象にしており、陸上性の哺乳類の群れにおける研究

間距離の計測に使用した。

は少ない[3]。特に 体サイズの大きな哺乳類の群れを対象にして個体の動きや群れ内での相対位置を定量的に捉えた研究は非常に少ない。そこで、本研究

3 ウマの個体間に働く力の数値シミュ
レーション

ではポルトガルのアルガ山に生息する野生のウマの群れを対象に個体間距離を調節している力の推力の数値シミュレーションのために、一個体の周定を目指した。

囲に重力のように”場”を導入し、距離が遠いと一般的にウマは極めて社会性の高い動物と考えられている。野生環境下において、ウマの個体群はAttraction として、距離が近いと Repulsion として働く力を以下のように定めた。

ハレム群とバチェラー群に分けられ、ハレム群は有蹄類としては珍しく恒常的な群れ構成を維持する。
群れのまとまりが安定して観察できることから個体間に働く力を推定するには、非常に適した動物だと言える。本研究では Attraction と Repulsion の相𝐹𝐹 = −𝑎𝑎(𝑟𝑟 − 𝑟𝑟( )− 𝑏𝑏(𝑟𝑟 − 𝑟𝑟( )*+ 𝑑𝑑𝑐𝑐

この関数において r は個体間の距離をあらわし、

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 25

ra は Attraction と Repulsion の境界値である。本研

topologicalmetric

究では、各群れの個体数と同数の点を 2 次元平面上

0.15にランダムに配置した状態から、F を個体間に作用

16-Kyoto

0.10

0.05

させ、一定ステップ後に点間の距離を測定した。こ0.00

れを観察で得られた画像の枚数分繰り返した。そし

0.15

cRMSE

0.05

(RMSE)をもとに比較し、a、b、c、d、ra の最適パラ

group17hyogoり得られた分布とピーク位置および残差平方和

0.10

17-Hyogoて、測定された個体間距離の分布を、野外観察によ16-Kyoto

17-Hyogo

17-Kyoto

0.00
0.15

メータセットを探索した。合計で各群れにつき

17-Kyoto

0.101171800 通りのモデルを試した。

0.05

0.00

4 相互作用モデル

6

7

8

distance D

9

10

2

3

neighbor K

4ある個体の力の作用範囲を定めるモデルはいく

Figure 1 インタラクションタイプおよびつかあるが、本研究ではメトリック距離およびトポ

その範囲におけるシミュレーション結果とロジカル距離を使用した。メトリック距離において実測値との誤差

は、個体は周囲 Rm 内に存在する他個体全てとイン
group

タラクションを持つ。トポロジカル距離においては、距離の近い順に Nt 個体とインタラクションを持つ。

16-Kyoto

17-Hyogo

17-Kyoto

****

4.50

****
n.s.

distance (BL)

5 結果

4.25

4.00

ウマの群れにおける任意の 2 個体の距離を最もよく説明するモデルは、メトリック距離を 7 body

3.75

length (BL) 以内に作用させるものだった (Figure1)。また、Attraction と Repulsion は 3.5-4 BL 以
内で逆転すると推定された (Figure 2)。

6 考察
これまで大きな体サイズ動物の群れにおける局所的行動ルールの解明に試みた研究は限られていた。本研究では、ドローン技術を用いた高精度計

3.50
16-Kyoto

17-Hyogogroup

17-Kyoto

Figure 2 推 定 さ れ た Attraction とRepulsion の逆転距離

6 参考文献[1] ID. Couzin, et al. "Collective memory andspatial sorting in animal groups." Journal of

測とスーパーコンピュータを用いたシミュレーシ

theoretical biology 218.1 (2002): 1-11.

ョンによって、ウマの個体間の距離を調節している[2] JE. Herbert-Read et al. "Inferring the rules of

インタラクションの推定をおこなった。その結果、interaction of shoaling fish." PNAS 108.46 (2011):

ウマにおいても、近傍の個体との局所的な相互作用

18726-18731.

によって距離を調節していることが示唆された。ま[3] S. Inoue, et al. "Spatial positioning of

た、この結果は、群れサイズにかかわらず局所的なindividuals in a group of feral horses: a case study

相互作用で群れのふるまいが説明できる可能性をusing drone technology." Mammal Research 64.2

提示しており、群れサイズの小さな動物種という観(2019): 249-259.

点からも興味深い知見を提供している。

26 − Academic Center for Computing and Media Studies, Kyoto University リモートセンシングなどの大規模環境データ解析による環境社会問題への応用、社会実装実験佐藤 建
国立高等専門学校機構一関工業高等専門学校未来創造工学科

1 宮城県伊豆沼の環境

湖沼でも冬季には湖面が部分的には凍結して鳥の行動も制限される。

宮城県伊豆沼内沼自然公園は宮城県北部に位置

広域の情報を取得できる UAV や IoT、衛星デする自然湖沼である。この伊豆沼周辺では毎年 10ータによるセンシングの応用は伊豆沼のような広万羽以上の渡り鳥が飛来して越冬しており、湿地い環境を調査するために有効な手段となる。湖沼保全に関するラムサール条約にも登録される国際内および湖沼周辺の植生の分布とその季節的な変的にも重要な湿地である。かつてはアサザやガガ化を把握する手法として IoT による湖面撮影デー

ブタ、ヒシなどの豊かな水性植物群落が広がって

タ、UAV への搭載が可能な農業用の簡易マルチス

いたが最近の 30 年程度で大きく環境が変化してペクトルカメラを応用した手法を提案し、伊豆沼いる。

を対象としたモニタリング観測を行った。これら

伊豆沼では最多時期には水面の 85%をハスが占めるなど環境の生物多様性が失われつつあるこのデータを機械学習、統計、画像処理手法を用いて解析をおこなった。

とが懸念されている。宮城県ではヒメビシ、ガガ
ブタは県の絶滅危惧種 I 類、アサザはⅡ類となっ

2 環境データ取得ており、環境省ではヒメビシがⅡ類、ガガブタ、アサザは準絶滅危惧種とされている（宮城県レッ

寒冷地域では低温環境での行動のリスク、積雪

ドリスト、2016; 環境省レッドリスト、2020）。

環境での広範囲移動の困難、交通・行動範囲の減

伊豆沼は春季から夏季にかけて植生が成長し、少から人間活動が制限される。このため UAV や

現在では上記のように夏季にハスによって湖面がIoT データによるセンシングの応用は有効な手段

ほぼ覆われる。ハスの繁茂は水質の悪化や浅底化となる。本研究では伊豆沼での計測実験のためににも寄与していることが分かっている。また伊豆遠隔カメラやUAVによる分光画像計測を行った。

沼の水質悪化については COD 値が観測されてい伊豆沼では環境省によって伊豆沼・内沼サンク

る湖沼の中で最も高いような状態となっている

チュアリセンターに設置されている遠隔カメラ(環境省、平成３０年度公共用水域水質測定結果)。

(Web カメラ)によるタイムラプス撮影が行われ湖沼環境の人為的復元およびそのための研究開発ている。伊豆沼の冬季環境データとこのタイムラが試みられており、ハスの刈り取りのための機器

プス画像を利用して伊豆沼の冬季環境情報の解析の開発と運用を行うほか、マコモなどの植栽が行を行った。撮影間隔は 1 時間に 1 回であり途中配われている。

置やカメラの変更があるもが 2002 年から現在ま冬季は渡り鳥が飛来する時期である。伊豆沼周で継続して撮影されている。

辺に数多く訪れる渡り鳥は湖面から潜水してハスUAV による撮影は 2019 年 5 月から 11 月にかの地下茎（レンコン）などの植物を食べている。

けて週 1 回程度で行われた。分類のために室内実

これらの渡り鳥は湖が完全に凍結するような場所

験により各植物の分光特性の計測実験も行っていをさけて来訪しているものの、伊豆沼など周辺のる。UAV にはマルチスペクトルカメラが搭載され、京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 27

撮影により可視から近赤外までの波長ごとの反射に使用したデータとは別のシーンから作成した検強度が得られる。

証用データセットを用いて判別モデルの精度検証

3 伊豆沼の凍結情報解析を行った。
精度検証の結果イネ科植物と浮葉植物の 2 群判

伊豆沼・内沼サンクチュアリセンターに設置さ別の正答率が 90%以上となり、2 群判別が十分可

れている遠隔カメラには 2002 年から現在まで大能であると言える。各植物の分類を行う 6 群判別

量の撮影画像が記録されており、これを管理するもトップ 1 分類結果が 60％となり詳細分類につ環境省自然環境局生物多様性センターから画像デいてもある程度の分類が可能であると考えられるータを提供いただいた。これに対し Superpixel

が、イネ科植生間の分類に困難があった。

によるセグメンテーションを行ったところ湖沼が凍結していないときは矩形での協会分割が行われ5 総括るのに対し、湖沼が凍結しているときは氷縁に沿って境界分割された。

宮城県伊豆沼において UAV やリモートカメラ

こうした画像特徴を利用し、さらにタイムラプ

画像等を用いた湖面環境情報の取得を行った、機ス画像を動画化して湖氷変動の確認もおこない湖械学習や画像処理などのデータ処理によるタイム

面凍結の日数の変化を調べた。集計を行ったとこ

ラプスカメラ画像解析から冬季凍結期間の定量化

ろ 2010 年代前半でも冬季(10～12 月)の凍結日数を行いその減少を示した。分光画像情報を利用しが 2010 年代後半では 40 日以上あったのに対し

た UAV による植生撮影画像の機械学習により夏

20 日程度と減少していることがわかった。

季の植生分類モデルを作成した。

冬季期間の湖面は 3 つの期間に分類される（i）.

 今後は冬季の湖面環境の分類の自動化のための初頭の時期。
ハスなどの植生の一部が湖面に残る。

統計・機械学習解析を行う。植生は成長段階によ

凍結はほとんど行らない。
(ii). 湖面凍結が起こり、り期間毎に植生の情報が変動することも考えられ時に数日以上の凍結が起きる。湖面の植生はなくるためそれらを考慮した分類を検証することなどなっている。(iii). 晩冬の時期となると湖面の凍結の課題があげられる。

はおこらず、湖面に植生が確認できない。

 伊豆沼の湖面環境の画像データは 2002 年から

以前は凍結後に数日以上継続して凍結している

一貫して計測される長期的なデータとなっており、ことが多かったが近年は凍結したその日や翌日に

年代ごとの湖氷変化を示すことが可能となった。

は融解していることも見受けられる。画像解析か

過去の画像には大型台風による水位の変動なども

らは凍結が起こりやすい領域も確認された。

記録されている。自然公園の管理や防災データ取得などのため継続的なデータの取得と保全が必要

4 UAV 撮影画像解析による植生分類であり、そうしたデータの解析により新しい知見を得ることが可能になる。

UAV による分光計測画像から植生判別モデルを作成することで、植生分布をマッピングする。

謝辞：本研究の一部はスーパーコンピュータ共同

各バンドの画像情報を用い画像の選択領域内につ研究制度(若手・女性奨励枠)によるものである。

いて特徴量を説明変数として機械学習を行う。計

伊豆沼における遠隔カメラ画像は自然環境局環境測画像群のマルチスペクトル情報からイネ科（ヨ省生物多様性センターから取得した。

シ、マコモなど）と浮葉植物（ハス、アサザ、ガガブタ、ヒシなど）の植生判別を行うための判別
モデルの作成を行った。
計測を行った測定日の日別で植物種ごとに 1 シーンを選択し、分光情報を利用した特徴量を判別
モデル学習用のデータセットとした。機械学習の手法はランダムフォレストである。モデルの作成

28 − Academic Center for Computing and Media Studies, Kyoto University

スーパーコンピュータ共同研究制度（大規模計算支援枠）研究報告

ϝογϡɾνϟωϧཚྲྀͷ௒ߴղ૾౓ DNS
ԣౢ ఩ ∗ ɼੁ໺ խ໵ †∗†੩Ԭେֶ ֶज़Ӄ޻ֶྖҬ ਺ཧγεςϜ޻ֶ‫ྻܥ‬੩Ԭେֶ େֶӃ૯߹Պֶٕज़‫ڀݚ‬Պ ޻ֶઐ߈ ਺ཧγεςϜ޻ֶίʔε

͸͡Ίʹ

1

periodicity imposed in x2 as well as in x1 & x3
x2

ಁաੑΛ༗͢Δ‫ڥ‬ք͸ɼਁಁচ΍थྛଳʹ୅ද͞
ΕΔΑ͏ʹɼಛʹࣗવքʹวࡏ͢Δ΋ͷͷɼͦͷӨ

x1

L3

x3

L2

‫ڹ‬ͷཧղ͸ෆಁաͳ‫ڥ‬քͱൺ΂ͯෆे෼Ͱ͋ΔɽஶऀΒ͸͜Ε·Ͱɼ‫׬‬શൃୡͨ͠νϟωϧཚྲྀͷ‫ڥ‬քFlowʹγϯϓϧͳಁաੑΛ෇Ճͨ͠‫ܥ‬Λର৅ͱͨ͠਺஋
࣮‫ݧ‬Λߦ͍ɼͦͷ఍߅ಛੑ΍ཚྲྀมௐʹ͍ͭͯ‫ݕ‬౼D͖ͯͨ͠ [1, 2, 3]ɽಛʹจ‫[ ݙ‬2, 3] Ͱ͸ɼνϟωϧ

W

S
L1

ཚྲྀͷ্ԼͷนΛɼϝογϡঢ়ͷಁաੑ‫ڥ‬քͰஔ͖(a)‫͑׵‬ɼ‫ڥ‬ք๏ઢํ޲ʹ΋प‫ظ‬ੑΛ՝ͨ͠‫ܥ‬ʢҎԼɼ“ϝογϡɾνϟωϧྲྀ” ͱ‫ݺ‬শʣΛର৅ͱͯ͠ɼͦ
ͷۭܺಛੑΛ‫ܥ‬౷తʹมԽͤͨ͞਺஋࣮‫ݧ‬Λߦ͖ͬͯͨɽຊใ͸ɼจ‫[ ݙ‬3] Ͱใࠂͨ͠਺஋࣮‫݁ݧ‬Ռʹؔ͠

P

P+∆P

ͯɼಛʹۭؒղ૾౓ґଘੑʹ͍ͭͯ‫ݕ‬౼Λߦ͏΋ͷ
Ͱ͋Δɽ
(b)

2

ର৅ྲྀΕɼ਺஋‫ࢉܭ‬๏ɼઌߦ‫ݚ‬ਤ 1: (a) ϝογϡɾνϟωϧཚྲྀͷུ֓ɽશͯͷํ‫Ͱڀ‬ಘΒΕͨओཁͳ஌‫ͱݟ‬ຊ‫ݚ‬ϝογϡঢ়‫ڥ‬ք͸ಉҰ‫ڥ‬քΛࢦ͢ɽ(b) (a) ͸͜ͷΑ޲ʹप‫ڥظ‬ք৚݅Λ՝͓ͯ͠Γɼ্Լʹ഑ஔ͞Εͨ͏ʹଟ਺ͷϝογϡঢ়ฏ൘‫ڥ‬ք͕౳ִؒʹฒΜͩ‫ܥ‬‫ڀ‬ͷҐஔ෇͚

ͷ‫ز‬Կֶత࠷খϢχοτʢ͢ͳΘͪ 1 ૚෼ʣΛ‫ࢉܭ‬ຊষͰ͸ɼຊ‫ڀݚ‬ͷઌߦ‫ͨ͋ʹڀݚ‬Δจ‫[ ݙ‬2, 3]ର৅ͱͨ͠ͱ΋ղऍͰ͖Δɽ

ͷུ֓ͱɼͦΕΒʹର͢Δຊ‫ڀݚ‬ͷҐஔ෇͚Λࣔ͢ɽ૚෼ʣͱղऍͰ͖Δɽ(a) Ͱ͸શͯͷํ޲ʹͦΕͧΕ

2.1

ର৅ྲྀΕɼ਺஋‫ࢉܭ‬๏ɼઌߦ‫Ͱڀݚ‬ಘΒΕͨओཁͳ஌‫ݟ‬प‫ظ‬௕ L1 , L2 , L3 Ͱप‫ڥظ‬ք৚݅Λ՝͓ͯ͠Γɼ(a)ͷ্Լ 2 ຕͷϝογϡঢ়ฏ൘͸࣮ࡍʹ͸ಉҰ‫ڥ‬քΛࢦ͢ɽฏ൘ؒ‫཭ڑ‬Λ 2δ (= L2 ) ͱͯ͠ɼओྲྀํ޲ʹ‫ݕ‬౼ର৅Ͱ͋Δϝογϡɾνϟωϧྲྀͷུ֓Λਤ

L1 = 3.84δ ɼԣஅํ޲ʹ L3 = 1.92δ ͷ‫ྖࢉܭ‬ҬΛ֬

1 ʹࣔ͢ɽ(a) ͸࣮ࡍʹ‫ݕ‬౼ର৅ͱͨ͠‫ܥ‬ͷུ֓Λද͢ɽ(a) ͸ɼ(b) ͷΑ͏ʹଟ਺ͷϝογϡঢ়ฏ൘‫ڥ‬քอͨ͠ɽϝογϡঢ়‫ڥ‬քͷ‫ܗ‬ঢ়͸ϝογϡૉࡐ෯ W ɼ͕౳ִؒʹઃஔ͞Εͨ‫ܥ‬ͷ‫ز‬Կֶత࠷খϢχοτʢ1Ծఆʣɼฏ൘ް D Ͱ‫ن‬ఆ͞ΕΔɽຊ‫Ͱڀݚ‬͸ฏ൘໘্ۭܺαΠζ S ʢۭܺ͸ x1 –x3 ฏ൘໘಺Ͱਖ਼ํ‫ܗ‬ঢ়ͱ

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 29ද 1: จ‫[ ݙ‬3] Ͱ‫ݕ‬౼ର৅ͱͨ͠ϝογϡঢ়ฏ൘‫ڥ‬քͷ‫ز‬ԿֶಛੑͱຎࡲϨΠϊϧζ਺ Reτ ͷҰཡɽM00-00 ͸ۭܺαΠζ S = 0ɼ͢ͳΘͪ‫ݹ‬యతͳෆಁաνϟωϧཚྲྀΛද͢ɽ

Case

ϕ

S/δ

W/δ
∞

S+
0

W+

S/∆x

W/∆x

∞

0

∞

M00-0000

M25-02
M25-04

0.25
0.25

0.02
0.04

0.02
0.04

3.0
7.7

3.0
7.7

3
6

M25-08

0.25

0.08

0.08

16.2

16.2

M25-16
M25-24

0.25
0.25

0.16
0.24

0.16
0.24

28.8
40.6

M39-05
M39-10

0.39
0.39

0.03
0.07

0.02
0.04

8.9
21.8

L1 /(S + W )L3 /(S + W )

Reτ

0

0

152.2

3
6

96
48

48
24

152.4
192.0

12

12

24

12

202.1

28.8
40.6

24
36

24
36

12
8

6
4

179.8
169.1

5.3
13.1

5
10

3
6

72
36

36
18

266.9
327.7

M39-15

0.39

0.10

0.06

28.4

17.1

15

9

24

12

284.3

M39-20

0.39

0.13

0.08

33.5

20.1

20

12

18

9

251.4

M44-05
M44-11

0.44
0.44

0.05
0.11

0.03
0.05

19.6
33.0

9.8
16.5

8
16

4
8

48
24

24
12

366.7
309.7

M44-21
M44-32

0.44
0.44

0.21
0.32

0.11
0.16

47.9
61.7

23.9
30.9

32
48

16
24

12
8

6
4

224.3
192.9

M56-06

0.56

0.06

0.02

25.1

8.4

9

3

48

24

418.8

M56-12
M56-24

0.56
0.56

0.12
0.24

0.04
0.08

39.2
59.0

13.1
19.7

18
36

6
12

24
12

12
6

326.5
245.7

M56-36

0.56

0.36

0.12

76.5

25.5

54

18

8

4

212.6

M69-10
M69-20

0.69
0.69

0.10
0.20

0.02
0.04

39.2
56.6

7.8
11.3

15
30

3
6

32
16

16
8

392.1
282.8

M77-14

0.77

0.14

0.02

49.2

7.0

21

3

24

12

351.3

M77-28

0.77

0.28

0.04

72.2

10.3

42

6

12

6

257.8

Λۭ͕ܺ઎ΊΔׂ߹Λදۭܺ͢཰ ϕ ≡ S 2 /(S + W )2ͱۭܺαΠζ S Λύϥϝʔλͱͯ͠‫ڥ‬քͷۭܺಛੑจ‫[ ݙ‬3] Ͱ͸ɼจ‫[ ݙ‬2] ͱൺ΂ۭͯؒղ૾౓Λશํ޲ʹ 2 ഒʢ͢ͳΘͪ 576 × 300 × 288 ͷཱํମ‫ܭ‬Λίϯτϩʔϧͨ͠ɽ͜Εͱ͸ಠཱʹްΈ D ʹΑͬࢉηϧͰ཭ࢄԽɼ࣌ؒ෼ղೳ΋ 2 ഒʹʣͱͨ͠ɽจ

ͯ‫ڥ‬քૈ౓ͷӨ‫ڹ‬Λ੍‫͖Ͱޚ‬Δ΋ͷͷɼຊ‫Ͱڀݚ‬͸‫[ ݙ‬3] ʹ͓͍ͯ΋ɼൺֱ‫ݕ‬౼ͷϕʔεͱͳΔྲྀΕ͸
Reb = 2300ɾReτ = 150 ͷ௨ৗͷνϟωϧཚྲྀͰ͋‫ڥ‬քͷಁաੑͷӨ‫ڹ‬ͷཧղʹٞ࿦Λूத͢ΔͨΊɼ

D = 0 ͱͨ͠ɽจ‫[ ݙ‬2] Ͱ͸ɼਤ 1(a) ͷ௚ํମྖҬʢ3.84δ × 2δ ×Δ΋ͷͷɼ͜͜Ͱ͸શͯͷ‫ݕ‬౼έʔεͰ Reb = 2300
Λॠ࣌ϨϕϧͰҰఆʹอͭΑ͏ʹɼओྲྀํ޲ͷྲྀΕ

1.92δ ʣΛ 288 × 150 × 144 ͷཱํମ‫ࢉܭ‬ηϧͰ཭ࢄԽͨ͠ɽஅ໘ฏ‫ྲྀۉ‬଎ Ub ͱνϟωϧ൒෯ δ ʹ‫ͮ͘ج‬όϧΫɾϨΠϊϧζ਺ Reb = 2300ɾຎࡲϨΠϊϧζ

‫ۦ‬ಈྗΛ੍‫ͨ͠ޚ‬ɽ͜ͷ৔߹ɼ‫ࢠ֨ࢉܭ‬෯͸‫ڥ‬ք͕਺ Reτ = 150 ͷ௨ৗͷνϟωϧཚྲྀʢ೪ணɾෆಁա

ͨέʔεʹ͓͚Δน࠲ඪͰͷ‫ࢠ֨ࢉܭ‬ղ૾౓͸‫ࢉܭ‬ͳ‫׈‬໘ฏ൘‫ڥ‬քͰߏ੒͞ΕΔʣΛൺֱ‫ݕ‬౼ͷϕʔε݁Ռʹґଘ͢Δʢ‫ޙ‬ड़ͷද 1 தͷ Reτ ʣɽෆಁաͳ৔߹ʢද 1 தͷέʔε M00-00ʣʹ͸શͯͷํ޲ʹҰ༷ʹ 1 ೪ੑ௕ͱͳΔ΋ͷͷɼಁաੑΛ෇͠ͱͯ͠ɼશͯͷ‫ݕ‬౼έʔεͰྲྀΕͷ‫ۦ‬ಈྗΛಉҰͷϝογϡঢ়ฏ൘‫ڥ‬ք͕ྲྀΕʹ‫͢΅ٴ‬Ө‫ڹ‬͸ Uhlmann஋ʹ‫ݻ‬ఆͨ͠ʢ͢ͳΘͪɼReτ ͸શέʔεͰ 150ʣɽ

ͷఏҊ͢ΔຒΊࠐΈ‫ڥ‬ք๏ [4] Λ༻͍ͯද‫ͨ͠ݱ‬ɽ਺͜ͷ৔߹ɼ‫ࢠ֨ࢉܭ‬෯͸શͯͷํ޲ʹҰ༷ʹ 2 ೪ੑ஋‫ࢉܭ‬๏ͷৄࡉʹ͍ͭͯ͸จ‫[ ݙ‬5, 6] Λࢀর͞Ε͍ͨɽ௕ͱͳΔɽ

30 − Academic Center for Computing and Media Studies, Kyoto Universityจ‫[ ݙ‬3] Ͱ‫ݕ‬౼ͨ͠ϝογϡ‫ڥ‬քͷ‫ز‬ԿֶಛੑΛ

0.08

͸ ϕ = 0 ͱද‫ه‬ʣ͔Β‫ͨͬ·ٻ‬ຎࡲ܎਺ͷ࣌ؒฏ‫ۉ‬

0.07

஋͸ C̄f = 8.61 × 10−3 Ͱ͋ΓɼIwamoto et al.[7] ʹ

0.06

ΑΔεϖΫτϧ๏ʹ‫ͮ͘ج‬νϟωϧཚྲྀͷ DNS ݁

0.05

Ռ (Reτ = 150) Ͱ͋Δ C̄f = 8.76 × 10−3 ͱ‫ ࠩޡ‬1%

Cf

0.04

ఔ౓ͰҰகͨ͠ɽਤ 2 ͔Β͸ɼྲྀΕ఍߅ͱ‫ڥ‬ք‫ܗ‬ঢ়

0.03

ͷؔ܎ʹ͍ͭͯɼҎԼͷ 2 ͭͷ܏޲͕ಡΈऔΕΔɽ(i) ۭܺ཰ ϕ ͕ߴ͍΄Ͳ఍߅͸૿Ճ͢Δɽ(ii)ϕ ͷ஋

0.02
0.01
0

͕ಉ͡Ͱ͋Ε͹ɼۭܺαΠζ S ͕খ͍͞΄Ͳ఍߅͕& ?=0/ ? = 0:25
5 ? = 0:39
0

/ ? = 0:44
. ? = 0:56
4 ? = 0:69

0.1

0.2

S=/

$ ? = 0:77‫·ڧ‬Δɽ͔͠͠ɼϕ = 0.25 ͱ ϕ = 0.39 ͷ৔߹ʹ্

0.3‫ه‬ͷ܏޲͔Βဃ཭͢Δσʔλ΋‫ݟ‬ΒΕͨɽ͢ͳΘͪɼ
ϕ = 0.25 ͷ৔߹ʹ͸ S/δ = 0.08ɼϕ = 0.39 ͷ৔߹ʹ͸ S/δ = 0.07 ͷ৔߹ʹ఍߅͕࠷΋‫͘ڧ‬ɼS Λ͞Β

0.4

ਤ 2: ຎࡲ܎਺ Cf ͷ࣌ؒฏ‫ۉ‬஋ʢϓϩοτʣͱඪ४ʹখ͘͢͞Δͱɼ૒ํͷέʔεͱ΋ʹ‫ݹ‬యతνϟωภࠩʢΤϥʔόʔʣͷ‫ڥ‬ք‫ز‬Կֶಛੑґଘੑɽͨͩ

ϧཚྲྀʢϕ = 0ʣͷσʔλʹ‫ؼ‬ண͢Δ܏޲͕ೝΊΒ͠ඪ४ภࠩ͸ 2 ഒʹ֦େͯ͠දࣔɽ

ΕΔɽ࣍ʹɼ֤έʔεͰ‫ͨͬ·ٻ‬ཚΕ‫౓ڧ‬ͷ‫ڥ‬ք๏ઢํͷ‫ڥ‬ք্Ͱͷ஋ʢ‫ڥ‬քۭܺͰͷྲྀମͷ޲੒෼ urms2

1.0

ਧ͖ग़͠/‫ࠐ͍ٵ‬Έͷ‫౓ڧ‬Λද‫ݱ‬ʣͷϝογϡ‫ز‬Կ

urms;+
2;meshֶಛੑʹର͢ΔґଘੑΛਤ 3 ʹࣔ͢ɽਤ 2 ͱ܏޲͸Ұக͓ͯ͠Γɼ͢ͳΘͪɼϝογϡۭܺͰͷਧ͖ग़͠/‫ࠐ͍ٵ‬Έ͕‫·ڧ‬Δͱओྲྀͷૄ௨ೳྗ͕௿Լ͢Δ͜

0.5

ͱ͕ಡΈऔΕΔɽ·ͨɼਤ 2 Ͱۭܺ཰ ϕ Λ‫ݻ‬ఆͨ͠& ?=0৔߹ʹ Cf ͕ϐʔΫ஋ΛͱΔۭܺαΠζɼ͢ͳΘͪ/ ? = 0:44 $ ? = 0:77

/ ? = 0:25 . ? = 0:56
5 ? = 0:39

0
0

0.1

ϕ = 0.25 ʹ͓͚Δ S = 0.08δ ͓Αͼ ϕ = 0.39 ʹ͓

4 ? = 0:69

͚Δ S = 0.07δ ͸ɼ͜ΕҎ্খ͞ͳۭܺ͸ྲྀମ͕ࣗ

0.2

༝ʹԟདྷͰ͖ͳ͍ྟք஋Λද͢ͱߟ͑ΒΕɼ͜ͷྟ

S=/

0.3

0.4

ਤ 3: ϝογϡঢ়‫ڥ‬քۭܺͰͷਧ͖ग़͠/‫ࠐ͍ٵ‬Έ‫ڧ‬౓ʢϝογϡঢ়ฏ൘‫ڥ‬ք্Ͱͷ‫ڥ‬ք๏ઢํ޲ཚΕ‫ڧ‬քۭܺαΠζ͕ԿͰ‫ن‬ఆ͞ΕΔ͔͸ඇৗʹ‫ڵ‬ຯਂ͘ɼࠓ‫ޙ‬ͷղ໌͕‫ظ‬଴͞ΕΔɽͷ‫ڥ‬ք‫ز‬Կֶಛੑґଘੑɽ౓ʣurms22.2
ද 1 ʹ·ͱΊͨɽۭܺ཰ ϕ Λ 6 ஈ֊ʢ0.25, 0.39, 0.44,

0.56, 0.69, 0.77ʣʹมԽͤ͞ɼͦΕͧΕͷۭܺ཰ʹରͯ͠ϝογϡִؒ S ʢۭܺαΠζʣΛ‫ܥ‬౷తʹมԽͤ͞ɼ‫ ܭ‬21 छྨͷϝογϡঢ়‫ڥ‬քΛର৅ͱͨ͠ɽจ‫[ ݙ‬3] Ͱ͸ۭؒղ૾౓Λ্͛ͨ͜ͱʹΑͬͯɼจ‫ݙ‬[2]ʢۭܺ཰͸ 0.25, 0.44, 0.56 ͷ 3 ஈ֊ɼ‫ ܭ‬9 छྨʣͱൺ΂ͯɼ‫ڥ‬քͷ‫ز‬ԿֶಛੑΛΑΓ޿ൣʹมԽͤ͞Δ͜ͱ͕ՄೳͱͳͬͨɽM00-00 ͸ S = 0 ͷ‫ݹ‬యతνϟωϧཚྲྀΛҙຯ͢Δɽ
ਤ 2 ͸ɼ֤έʔεͰಘΒΕͨຎࡲ܎਺ Cf ≡ 2u2τ /Ub2

ʢuτ ͸น໘ຎࡲ଎౓Λද͢ʣͷ࣌ؒฏ‫ۉ‬஋ͱඪ४ภ
ࠩͷϝογϡ‫ز‬Կֶಛੑʹର͢ΔґଘੑΛද͢ɽຊ਺஋࣮‫ݧ‬ͷ M00-00ʢ‫ݹ‬యతνϟωϧཚྲྀɼਤ 2 Ͱຊ‫ڀݚ‬ͷҐஔ෇͚

ද 1 ͔Β΋ಡΈऔΕΔΑ͏ʹɼۭܺ཰ ϕ Λ͋Δ஋ʹ‫ݻ‬ఆͨ͠৔߹ʹ࠷΋ۭܺαΠζ S ͕খ͍͞έʔεʢ͢ͳΘͪɼM25-02, M39-05, M44-05, M56-06, M69-

10, M77-14ʣʹ͓͍ͯ͸ɼۭܺαΠζ S ·ͨ͸ϝογϡૉࡐ෯ W ͸ߴʑ‫ࢉܭ‬ηϧ 3 ·ͨ͸ 4 ͭ෼ʹա͗ͣʢද 1 ͷ S/∆x, W/∆x ͕ͦΕΒʹରԠʣɼ‫ڥ‬քۙ๣Ͱ͸ྲྀ଎͸஗͍͜ͱΛߟྀͯ͠΋ɼ͜ͷۭؒ෼ղೳʹΑΔϝογϡۭܺ಺΍ૉࡐ෯पΓͷྲྀಈͷ࠶‫ݱ‬ੑʹ͸‫ݒ‬೦͕࢒Δɽ͜ͷ఺ʹ͍ͭͯ‫ݕ‬౼͢ΔͨΊɼຊ‫Ͱڀݚ‬͸ಛʹɼ࠷΋ຎࡲ܎਺ Cf ͕େ͖ͳ M56-06ɼ͓Αͼ࠷΋ຎࡲ܎਺͕খ͍͞ʢ‫ݹ‬యతͳνϟωϧཚྲྀʹ࠷΋͍ۙʣM25-02 ͷ 2 έʔεʹ͍ͭͯɼۭ࣌
ؒղ૾౓Λจ‫[ ݙ‬3] ʹରͯ͠͞Βʹ 2 ഒʢ͢ͳΘͪɼ

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 31

ਤ 1(a) ͷ௚ํମྖҬΛ 1152 × 600 × 576 ͷཱํମ‫ܭ‬

0.08ࢉηϧͰ཭ࢄԽɼ࣌ؒ෼ղೳ΋ 2 ഒʹʣͱͨ͠਺஋

0.07

࣮‫ݧ‬Λߦ͏ɽ

0.06
0.05

ωϧཚྲྀʢ‫ڥ‬ք͕ෆಁաʣʹ͓͍ͯ͸͍ΘΏΔϛχ

0.04

ϚϧɾϢχοτ [9] ʹ͍ۙ৚݅Ͱ͸͋Δ΋ͷͷɼද 1

0.03

Cf

ͳ͓ɼ3.84δ ×2δ ×1.92δ ͷ‫ྖࢉܭ‬Ҭ͸ɼ௨ৗͷνϟதͷҰ෦ͷέʔεΛର৅ͱͯ͠ɼओྲྀํ޲ͷ‫ྖࢉܭ‬0.02

ҬΛ 2 ഒʢL1 = 7.68δ ʣɼ͋Δ͍͸ԣஅํ޲ͷ‫ࢉܭ‬0.01

ྖҬΛ 2 ഒʢL3 = 3.84δ ʣͱͨ͠‫ࢉܭ‬Λ࣮ࢪ͠ɼ͜

0/ ? = 0:25

. ? = 0:56

ΕΒਫฏํ޲ͷ‫ྖࢉܭ‬ҬαΠζ͸ 2.1 અͷٞ࿦ʹ͋·ΓӨ‫ͱ͍͜ͳ͠ڹ‬Λ֬ೝͨ͠ [8]ɽ·ͨɼ‫ڥ‬ք๏ઢ
ʢx2 ʣํ޲ʹରͯ͠ɼਤ 1(b) ͷ‫ز‬Կֶత࠷খϢχο

0

0.1

0.2

S=/

0.3

0.4

τʢ1 ૚෼ʣΛର৅ͱͯ͠प‫ڥظ‬ք৚݅Λ՝͢ʢਤਤ 4: ຎࡲ܎਺ Cf ͷ‫ڥ‬ք‫ز‬Կֶಛੑґଘੑ: ੺,1(a)ʣ͜ͱʹ͍ͭͯ΋ɼ2 ૚෼ʹରͯ͠प‫ڥظ‬ք৚݅

576 × 300 × 288 ͰಘΒΕͨ݁Ռʢจ‫[ ݙ‬3]ʣ; ࠇɼ1152 × 600 × 576 ͰಘΒΕͨ݁Ռʢຊ‫ڀݚ‬ʣɽΛ՝ͨ͠‫ࢉܭ‬Λද 1 தͷҰ෦ͷέʔεΛର৅ͱͯ͠
࣮ࢪ͠ɼ2.1 અͷٞ࿦ʹ΄ͱΜͲӨ‫ͱ͍͜ͳ͠ڹ‬Λ֬ೝࡁΈͰ͋Δ [8]ɽ਺஋࣮‫݁ݧ‬Ռͱߟ࡯ਤ 4 ͓Αͼਤ 5 ͸ͦΕͧΕਤ 2 ͓Αͼਤ 3 ʹରԠurms;+
2;mesh

3

1.0

0.5

͢Δ΋ͷͰɼࠓճ৽ͨʹ 1152 × 600 × 576 ͷ‫֨ࢉܭ‬ࢠΛ༻͍ͯ਺஋࣮‫ݧ‬Λߦͬͨ M56-06 ͱ M25-02 ͷ

2 έʔεʹ͍ͭͯɼจ‫[ ݙ‬3] ͷ‫݁ࢉܭ‬ՌͱൺֱΛߦͬͨ΋ͷͰ͋Δɽۭؒղ૾౓Λ্͛ͨ݁ՌɼM56-06 Ͱ͸ຎࡲ఍߅܎਺ Cf ͕௿ԼɼM25-02 Ͱ͸΍΍૿Ճͨ͠ɽ·ͨɼM56-06 Ͱ͸ਧ͖ग़͠/‫ࠐ͍ٵ‬Έ‫ ౓ڧ‬urms
2͸௿ԼɼM25-02 Ͱ͸૿Ճͨ͠ɽ͜ΕΒͷ݁Ռ͸ɼจ‫[ ݙ‬3] ͷۭؒղ૾౓Ͱ͸ɼগͳ͘ͱ΋ࠓճ‫ݕ‬౼ର৅ͱͨ͠ 2 έʔεʹؔͯ͠ɼఆྔతͳٞ࿦͕Մೳͳσʔλ͸ಘΒΕ͍ͯͳ͍͜ͱΛҙຯ͢ΔɽଞํͰɼจ‫ݙ‬/ ? = 0:25 . ? = 0:56
0
0

0.1

0.2

S=/

0.3

0.4

ਤ 5: ϝογϡঢ়‫ڥ‬քۭܺͰͷਧ͖ग़͠/‫ࠐ͍ٵ‬Έ‫ڧ‬ͷཚΕ‫౓ڧ‬ͷ‫ڥ‬ք‫ز‬Կֶಛੑґଘੑɽ: ੺,౓ urms2

576 × 300 × 288 ͰಘΒΕͨ݁Ռʢจ‫[ ݙ‬3]ʣ; ࠇɼ1152 × 600 × 576 ͰಘΒΕͨ݁Ռʢຊ‫ڀݚ‬ʣɽ[3] ͓Αͼจ‫[ ݙ‬2] ͰಘΒΕͨɼຎࡲ܎਺ Cf ͱϝοͷγϡঢ়‫ڥ‬քۭܺͰͷਧ͖ग़͠/‫ࠐ͍ٵ‬Έ‫ ౓ڧ‬urms
2ࢀߟจ‫ݙ‬

૬ؔؔ܎ɼ͢ͳΘͪɼϝογϡۭܺͰͷਧ͖ग़͠/‫ٵ‬[1] Yokojima, S., Eﬀects of wall permeability on
wall-bounded turbulent ﬂows, J. Phys. Soc. Jpn.,
Vol.80, 033401, (2011).

͍ࠐΈ͕‫·ڧ‬Δͱओྲྀͷૄ௨ೳྗ͕௿Լ͢Δɼͱ͍͏஌‫ݟ‬͸ຊ਺஋࣮‫݁ݧ‬ՌΛՃ͑ͯ΋ͳ͓੒ཱ͓ͯ͠
Γɼਤ 4 ͔Β‫ݟ‬ग़͞Εͨ 2 ͭͷ܏޲ʢ(i) ۭܺ཰ ϕ ͕ߴ͍΄Ͳ఍߅͸૿Ճ͢Δɽ(ii)ϕ ͷ஋͕ಉ͡Ͱ͋Ε͹ɼʣ͓ΑͼɼۭܺαΠζ S ͕খ͍͞΄Ͳ఍߅͕‫·ڧ‬Δɽ

S/δ ≈ 0.07 − 0.08 ۙลʹɼྲྀମ͕ࣗ༝ʹߦ͖དྷ͕ՄೳͳۭܺαΠζͷྟք஋͕ଘࡏ͢Δͱ͍͏஌‫ݟ‬͸ͦ
ͷ··੒Γཱͭ͜ͱ΋ߟ͑ΒΕɼࠓ‫ޙ‬ͷ‫ݕ‬౼՝୊Ͱ͋Δɽ

32 − Academic Center for Computing and Media Studies, Kyoto University[2] ԣౢ ఩, ཚྲྀதͷಁաੑ‫ڥ‬քͷ఍߅ಛੑ, ͳ͕Ε,
Vol.35, No.4, pp.307–312, (2016).
[3] ੁ໺ խ໵, ԣౢ ఩, ಁաੑϝογϡঢ়‫ڥ‬քΛ༻͍ͨνϟωϧཚྲྀͷ DNS, ೔ຊ‫ػ‬ցֶձ ୈ 97 ‫ྲྀظ‬ମ޻ֶ෦໳ߨԋձ ߨԋ࿦จू, OS 8-06, (2019).
[4] Uhlmann, M., An immersed boundary methodwith direct forcing for the simulation of particulate ﬂows, J. Comput. Phys., Vol.209, pp.448–476,(2005).

[5] ԣౢ ఩, ʮϝογϡঢ়ಁաੑฏߦฏ൘ؒཚྲྀͷ௚઀਺஋γϛϡϨʔγϣϯʯʹؔ͢ΔϓϩάϥϜߴ౓Խ,‫౎ژ‬େֶ ֶज़৘ใϝσΟΞηϯλʔ શࠃ‫ڞ‬ಉར༻൛޿ใ, Vol.14, No.2, pp.48–51, (2016).
[6] ԣౢ ఩, ʮϝογϡঢ়ಁաੑฏߦฏ൘ؒཚྲྀͷ௚઀਺஋γϛϡϨʔγϣϯʯʹؔ͢Δେ‫ن‬໛਺஋‫ࢉܭ‬, ‫ژ‬౎େֶ ֶज़৘ใϝσΟΞηϯλʔ શࠃ‫ڞ‬ಉར༻൛޿ใ, Vol.14, No.2, pp.60–63, (2016).
[7] Iwamoto, K., Suzuki, Y. & Kasagi, N., Reynolds
number eﬀect on wall turbulence: toward eﬀective
feedback control, Int. J. Heat Fluid Flow, Vol.23,
pp.678–689, (2002).
[8] ੁ໺ խ໵, ԣౢ ఩, ϝογϡঢ়ಁաੑ‫ڥ‬քΛ༗͢Δ
νϟωϧཚྲྀ DNS ʹ͓͚Δ‫ྖࢉܭ‬ҬαΠζͱۭؒղ૾౓ґଘੑ, ೔ຊ‫ػ‬ցֶձ ౦ւࢧ෦ ୈ 69 ‫ظ‬૯ձɾߨԋձ, (2020).
[9] Jimenez, J. & Moin, P., The minimal ﬂow unit
in near-wall turbulence, J. Fluid Mech., Vol.225,
pp.213–240, (1991).

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 33

プログラム高度化共同研究報告

水・熱移動を考慮した 3 次元剛体バネモデルによるコンクリート構造物の長期変形・ひび割れ進展挙動および動的破壊挙動の高解像度シミュレーション
山本 佳士
法政大学 デザイン工学部 都市環境デザイン工学科

1 はじめに手法は、地震動等の作用を受ける RC 構造物の耐力やマクロな変形性能を精度良く予測することには

近年、インフラ構造物の適切な維持管理、あるい

成功しているものの、前述のように複数のひび割れ

は長寿命化という観点から、様々な要因により劣化を平均化して扱っているため、個々のひび割れの位したコンクリート構造物の残存性能評価が重要な

置、幅、角度等の詳細情報を評価することができな

課題になってきている。特に、原子力発電施設等のい。したがって、様々な要因で発生、進展するひび

重要構造物では、精度の高い数値解析モデルにより割れに起因して局所的に変化する透気性、透水性を長期挙動予測および性能評価を行うことが望まれ評価することは難しい。その他、高解像度メッシュている。このような背景から、近年、コンクリート

を用いた手法 2)や、X-FEM の適用 3)などにより、構造物の各種応答予測を対象としたベンチマーク

FEM を用いて RC 部材の不連続挙動を直接扱う手実験およびブラインド解析が国内外で精力的に実法が提案されているものの、まだ検討段階である。

施されている。例えば、フランス電力（EDF）は、

一方、著者らは、現在の数値解析技術では評価が実規模の原子炉格納容器のモックアップを製作し、難しい、前述のひび割れ詳細情報および RC 構造物

施工段階から経時変形・ひび割れ進展挙動、透気性の破壊・倒壊挙動を精度良く再現することを目的と（遮蔽性）を計測するとともに、それらを対象としして、剛体バネモデル 4)（以下、RBSM）を基本とたオープン参加のブラインド解析を実施し、現在のした解析手法・モデルの開発を行ってきた。RBSM

数値解析技術の予測性能を評価する試みを行っては対象を剛体要素と要素間のバネで離散化する手いる（VeRCoRs プロジェクト 1)）。なお、同プロジ法であり、開発当初は規則配列要素が用いられていェクトでは、最終的にシビアアクシデント下におけたが、Voronoi 分割を用いたランダム形状要素（図-る予測精度の確認を行うために、前述の実規模モッ1）の適用 5)により、あらかじめ破壊面を想定する必クアップに対して火災を想定した温度変化作用下

要がなくなり、コンクリートのひび割れ進展予測手の応答を計測する予定になっている。

法として有用なツールになってきている。著者らは、コンクリート構造物の代表的な数値解析手法と境界面頂点

しては非線形 FEM が挙げられ、これまでに多くの

積分点

境界面重心

研究開発と適用が行われてきた。コンクリート材
料・構造の力学挙動の特徴としては、ひび割れ、ひび割れ面のせん断伝達、鉄筋等、異材料間の付着す
べり等の不連続挙動が挙げられる。連続体力学を基本とした FEM で上記のような不連続挙動を表現するためには工夫が必要であり、現在では、複数のひび割れおよび鉄筋を含んだ鉄筋コンクリート（RC）
要素の力学応答を空間平均化して評価する、分散ひび割れ・分散鉄筋モデルの適用が主流である。この34 − Academic Center for Computing and Media Studies, Kyoto University(a) Voronoi 分割(b) 2 要素力学モデル

図-1 RBSM の概要

これまでに各種応力下におけるコンクリートのひ用いることで、各種応力下におけるコンクリートのび割れ詳細情報、ひずみ軟化・局所化挙動等を定量

ひずみ軟化・局所化挙動を再現することができる。

的に再現できる構成モデルの開発、付着すべり挙動

鉄筋は梁要素でモデル化し、鉄筋軸方向および鉄筋を考慮した鉄筋モデルの導入、動的解析への拡張等

軸直交方向のバネからなるリンク要素を介してを行ってきた 6) 7) 8)。さらに、水・熱移動解析と前述

RBSM 要素に結合している 7)。鉄筋軸方向のバネにの RBSM とを連成させた手法を開発し、コンクリー

非線形性を導入することで鉄筋－コンクリート間

ト構造物の、施工段階から供用期間中の任意の時点の付着すべり挙動を再現することができる。

までの、1)水和反応、温度応力、乾燥収縮、クリープ

水分移動および熱移動は、拡散方程式を離散化し等に伴う変形・ひび割れ進展性状、および任意の時

た 2 節点 1 次元 FEM 要素により表現する。同要素

点における 2)透気性（遮蔽性能）、3)耐震性能、4)耐は、Voronoi 分割の母点および Voronoi 面重心、火性能を、再現可能な解析システムを開発している

Voronoi 辺上の点を結ぶネットワーク上に配置され

9) 10) 11) 12)。なお、RBSM を用いて、コンクリートのる（図-2）。RBSM では、ひび割れの幅を要素間の相軟化・変形局所化挙動を再現するために、数 mm か

対変位を用いて定義しているが、このひび割れ幅のら数 10mm スケールのメッシュを使用する必要が拡大に応じて、水分・熱移動を表現する 1 次元 FEMある。すなわち、手法の性能を維持した上で実規模

要素の面積を増加させることにより、ひび割れの増の構造物の性能評価、あるいは前述のベンチマーク

加に応じて、局所的かつ不連続に見かけの物質拡散実験（数 m～数 10m スケール）の数値解析を効率

係数が変化する挙動を表現できる。

的に行うためには、解析プログラムの高速化が必須

本解析プログラムでは、時間ステップ毎に、まず、である。このような背景から著者は本プログラム高水分移動および熱移動を解き、それぞれの解析から度化事業に応募した。以下に、提案する解析手法の

各要素の水分量の変化（すなわち乾燥の進展）およ

概要とその適用例、プログラムの高度化の成果、おび温度の変化を算定する。水分の変化から乾燥収縮

よび今後の課題について述べる。

ひずみを、温度変化から熱膨張（収縮）ひずみを評価し、それらと等価なバネひずみを初期ひずみとし

2 解析手法の概要て RBSM に与えて全体の変形挙動を解いている。
時間ステップ毎に、上記の計算を繰り返していくこ

RBSM は前述のように対象を剛体要素と要素間とで、弱連成の枠組みで簡単に温度変化や乾燥収縮のバネで離散化する手法である。剛体要素は要素内に伴う変形・ひび割れ進展解析を行っている。また、の自由度評価点に、並進 3 方向、回転 3 方向の自由

任意の時刻から、単調載荷解析あるいは地震応答解度を持ち、剛体変位場を仮定して 2 つの要素が共有

析を実施することも可能である。

する境界面上の任意点のバネのひずみを評価する（図-1）。バネに引張軟化、すべり軟化等の非線形構

3 解析手法の適用例

成モデルを適用して、仮想仕事の原理を用いて増分形の線形化剛性方程式を導出し、ニュートンラプソ

前章で示したように、本解析プログラムは、地震・ン法等を用いて非線形問題を解いている。定式化お衝撃作用下の動的応答解析や水分・熱移動に伴う乾よび提案する構成モデル、材料パラメータの決定方

燥収縮、高温作用下の応答、各種劣化進展、長期変法等の詳細は参考文献 6)を参照いただきたい。前述

形・ひび割れ進展解析など、広範な問題の数値解析のように数 mm から数 10mm スケールのメッシュが可能である。本章では本解析プログラムの適用例

を使用し、提案する構成モデルと材料パラメータをの概略を示す。

3.1 地震・衝撃・火災作用を受けるコンクリー
ト構造物の応答シミュレーション
図-3 および図-4 に、正負交番繰返し載荷を受ける
図-2 Voronoi 分割上の FEM 要素配置

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 35

Unit : mm
200

200
100

Heating area

0.1mm

境界条件

0.5mm

図-3 繰返し正負交番載荷を受ける RC 壁のひび割れ詳細情報評価（ConCrack21) 7)）

Symmetric
boundary

Outer surface
boundary

温度分布

蒸気圧分布

図-4

RC 建屋の模型供試体（1/4 スケール）の振動実験の解析（SMART201313)）

爆裂開始
爆裂開始
13.2 msec 後

爆裂開始
28.8 msec 後

変形性状
図-5 高速衝突を受ける RC 板の局部破壊解析 8)

図-6 高温負荷を受けるコンクリートの爆裂挙動シミュレーション 12)

RC 壁の応答シミュレーション、および RC 建屋の内の蒸気圧が急激に増加し、蒸気圧の作用を受けて地震応答シミュレーションの例を示す。本解析は、コンクリートが爆裂している。

ひび割れ詳細情報、3 次元地震応答等に対する予測
技術の現状を把握することを目的として実施された 、 ベ ン チ マ ー ク 解 析 (ConCrack21) お よ び
SMART201313))で実施したものである。ベンチマー
ク解析の結果およびその評価は、参考文献 1)7)において紹介されているが、提案手法は、マクロな荷重
変位応答に加え、ひび割れの幅、角度、間隔などのひび割れ詳細情報、破壊モードを良く再現していた。
図-5 には、剛飛翔体の高速衝突(飛翔体：50g の鋼球、衝突速度：300～500m/s)を受ける RC 板の解析の例8)を示す。参考文献に示すように提案モデルは、破壊限界速度、破壊モードに加え、高速カメラによって観察された、ひび割れ発生のタイミング、ひび割れの進展速度までを再現している。図-6 には、コンク
リート供試体表面に高温を与えた解析の一例を示す 12)。解析では、蒸気圧および熱の移動をシミュレ

3.2 コンクリート構造物の長期変形・劣化挙動
シミュレーション
図-7 に乾燥収縮によりひび割れが生じた RC 壁部
材の正負交番載荷実験の解析結果を示す 9)。解析は、
乾燥収縮ひび割れが生じた RC 壁の剛性を過小評価する傾向にあるものの、乾燥収縮ひび割れの有無によって変化する載荷時のひび割れ分散性状や、破壊モード、最大荷重を妥当に再現していた。図-8 にはアルカリ骨材反応（ASR）により劣化した RC 壁部
材の正負交番載荷実験の解析結果を示す 10)。なお、
本解析は、経済協力開発機構／原子力機関(OECD/NEA)の原子力施設安全委員会(CSNI)の下で実施されたベンチマーク解析(ASCET phase II)
で実施したものである。提案モデルは、ASR 劣化にートする２層の FEM 要素と RBSM を連成させて伴うひび割れが、載荷時に発生・進展するひび割れいる。高温負荷により、表面からコンクリートの温健全供試体と ASR 劣化供試体の間で、終局変位、度が上昇し、さらに温度上昇に応じてコンクリート

36 − Academic Center for Computing and Media Studies, Kyoto Universityの分散性や角度に与える影響を評価し、結果として破壊モード、局所化位置が変化する様子を参加チー(a) 相対含水率(c)荷重作用下における終局時ひび割れ（湿潤養生供試体）(b) 乾燥収縮ひび割れ(d)荷重作用下における終局時ひび割れ（乾燥供試体）

図-7 乾燥収縮によりひび割れが生じた RC 壁の正負交番載荷解析 9)

ムの中で唯一再現していた。

4 プログラム高度化の成果と今後の課題
本プログラム高度化事業による、解析プログラムの計算速度向上効果を評価するために、図-7 に示したシミュレーションを、改善前のプログラムおよび改善後のプログラムを使用して実施した。なお、プ
ログラム改善のための主な実装の内容は、プログラ
ム内のスレッド並列化領域の拡大である。RBSM の要素数を約 10000 および約 30000 とした、2 種類の計算規模で検討を行ったところ、いずれの計算規模でも計算速度は約 2 倍程度向上した。
解析プログラムの改善により、ここで検討した範囲内では約 2 倍程度の速度の向上が確認されたが、前述のような実規模構造レベルのシミュレーショ
ンのためにはさらなる高速化、効率化が必要である。
今後は、詳細なひび割れ・損傷情報を評価したい領域では RBSM を用い、その他の領域では FEM を用いて粗い解像度で効率的に解析を行う、RBSM と
FEM とを結合した解析手法の開発とその高度化を行っていく予定である。

参考文献
1) J. Mazars et al.：Behavior and Assessment of Massive

Structures : An Overview of the French Research Programs
CEOS.fr and VeRCoRs, Proceeding of FraMCoS-9, (CDROM), 2016.
2) M. Kurumatani et al. ：Damage model for simulating(a) 健全供試体(b) ASR 劣化供試体

図-8 ASR 劣化が生じた RC 壁の正負交番載荷解析 10)（ASCET phase-II）chloride concentration in reinforcedconcrete with internalcracks, Cement and Concrete Composites, Vol.84, pp.62-73,
2017.
3) C. Chung et al. ：Diagonal Crack Propagation Analysis ofReinforced Concrete Beams Using XFEM, Proceeding of
ASEM 17, 2017.
4) Kawai, T. : New Element Models in Discrete Structure
Analysis, 日本船舶学会論文集，第 141 号，1977.
5) 竹内則雄ら：ボロノイ分割を用いた RBSM による無筋コンク
リート梁の寸法効果解析，構造工学論文集，40A，519-527，
1994.
6) 山本佳士ら：3 次元剛体バネモデルによるコンクリート供試体の圧縮破壊解析，土木学会論文集 E，Vol.64 No. 4，pp. 612630，2008.
7) Y. Yamamoto et al. : Crack propagation analysis ofreinforced concrete wall under cyclic loading using RBSM,
European Journal of Environmental and Civil Engineering,18(7), 780–792, 2014.
8) 山本佳士ら：RBSM による高速衝突を受けるコンクリート板の局部破壊挙動評価に関する基礎的研究，構造工学論文集
Vol.61A，2015.
9) H. Sasano et al. : Impact of Drying on Structural
Performance of Reinforced Concrete Shear Walls, Journal ofAdvanced Concrete Technology, Vol.16, No.5, pp.210-232,
2018.
10) 杉本啓太：ひび割れ詳細情報の再現に着目したASR-RBSMの開発および ASR 劣化した RC 部材の破壊機構評価，名古屋大
学修士論文，2018.
11) P. SRIMOOK : Development of Simulation Method for Air
Leakage Through Concrete Containment Structure underService and Accidental Condition using RBSM 名古屋大学
修士論文，2020.
12) H. Nakamura et al. : Simulation about Effect on Explosion
Spalling of Thermal Stress and Vapor Pressure, 2ndInternational Workshop on Concrete Spalling due to Fire
Exposure, pp.213-219, 2011.
13) B. Richard et al. : SMART 2013 Lessons learned from theinternational benchmark about the seismic margin
assessment of nuclear RC buildings, Engineering
Structures, 161, 207-222, 2018.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 37

水を推進剤とした超小型マイクロ波放電式プラズマ源を対象とした

3 次元プラズマ粒子計算
鷹尾 祥典**横浜国立大学 大学院工学研究院 システムの創生部門

1 はじめに要があり、性能低下を招いている。この性能改善ため中和器から効率的に電子を引き出すことが必近年、大学や中小企業など比較的小規模な組織要であるが、小さい放電室のため実験だけですべでも運用が可能な超小型衛星の研究開発が盛んにての現象を把握することは難しい。そのため、数行われている。超小型衛星は短期間・低コストで値計算を用いたプラズマ放電ならびに電子引き出製作できるため、教育や技術実証、商用利用など

し機構の解析が求められている。

幅広いミッションに用いられている。既に年 300

これまでに、プラズマ生成から電子引き出しま基が打ち上げられ、今後もさらに伸びると予測さでを統一的に解析できる 3 次元プラズマ粒子計算

れている [1]。

モデルを構築し [3]、まずは Xe を対象に中和器か

これまでは技術的な制約から推進機を搭載しならの電子引き出し機構の解析を行ってきた [4, 5]。

い小型衛星がほとんどであったが、最近は小型推当時、１つの条件当たり十分な時間平均を取って進機の開発も進み、搭載する例が増えてきている。

統計的に信頼できる解析結果を得るためには、最

例えば、東京大学と NExt generation Space system

低でも 1 ヵ月以上（～数ヵ月）計算を続ける必要

Technology Research Association (NESTRA) によっがあり、
平成 29 年度プログラム高度化支援事業第て 開 発 さ れ た 50 kg 級 の 超 小 型 衛 星

2 期を利用するに至った。

HODOYOSHI-4 には小型イオン推進システム

前回の高度化支援により、システムB において、(MIPS: Miniature Ion Propulsion System) が搭載さ

72 スレッドで 23 倍のスケーラビリティを得られれ、2014 年10 月に宇宙空間での運用に成功した。

た。今回、分子性の水を利用することで、考慮す

この MIPS にはイオン源と中和器にマイクロ波アる反応および粒子種が大幅に増え [6]、また、よンテナとリング状磁石を用いた ECR（Electron

り細かい現象を把握するためセル幅を半減させた

Cyclotron Resonance）放電を採用している [2]。

都合、更なる高度化が必要となった。

MIPS は従来から利用されているXe ガスを推進剤に利用する都合、
これ以上の小型が難しい一方、

2 数値計算手法

今でも大半の小型衛星は 10 kg 以下であり、更なる小型化が求められている。そこで、高圧ガス系

推進剤が Xe から水に変更となっても、計算対が不要な液体推進剤として、安価、かつ、完全無象の構造は同じである。また、解析手法について

毒で安全な水を新たに推進剤として利用することも同様であるが、以下に簡単に示す。

で推進システムの小型化を目指している。

図 1 に本研究が対象とする小型マイクロ波放電

一方、水分子は Xe と比較して放電を維持する

式中和器の計算領域を示す [7]。金属面とプラズことが難しく、特にイオンビームを電気的に中和

マ源の境界においてその中心を原点とする直交座する中和器には Xe より多くの流量を投入する必標系を設定した。プラズマ源を 20×20×4 mm3 の

38 − Academic Center for Computing and Media Studies, Kyoto University

2 mm 4 mm
Antennayx

Plasma4.2-GHzMicrowave
Excitation plane(a) z-y/x plane (x/y = 0 mm)zBN

Orificeplate(Metal)

Metal

10 mm
Magnet(Metal)

図 2: 計算フローチャート(b) x-y plane (z = 1.0 mm) : antenna

粒子は電子および H2O+, OH+, H+, H−, O−, OH−とす

Xenon Plasma

る。(ii) 中性粒子は300 K のマクスウェル分布で、Ring-shaped
Antenna

時間的・空間的に一様とする。(iii) 荷電粒子と中性粒子の衝突は、弾性・励起・電離・解離・電子yx

付着の各衝突および電荷交換衝突を考慮する。

20 mm

z

(iv) マイクロ波による磁場は、永久磁石による静磁場と比べて小さいため無視する。(v) 本研究ではマイクロ波の電力が小さいためにプラズマ電流

20 mmの影響を無視する [8]。
図 2 に本計算のフローチャートを示す。はじめ(c) x-y plane (z = 4.0 mm) : orifice plateに初期条件を与え、定常状態におけるマイクロ波
y𝝓𝝓 2.2 mmzxに よ る 電 磁 場 を 得 る た め に Finite-difference
Time-domain (FDTD) 法によりマクスウェル方程式を解く。このとき、時間ステップは∆tEM =1.49×10−13 s（4.2 GHz マイクロ波の 1/1600 周期）である。次に、PIC/MCC 法により粒子の運動を、
時間ステップ∆te = 5.95×10−12 s（マイクロ波の 1/40

Orifice forElectron Extraction

周期）の下で解く。
マイクロ波電場の時間変化は、

図 1: 解析対象とする電子源の構造 [7](JSASS-2016-4687 Ⓒ日本航空宇宙学会)
領域とし、その内部（z = 1.0 mm）にはリング状ア

先ほど求めた定常状態のマイクロ波電場 E0 を用
いて、EEM = E0 cos(ωt) で与えられる。また、永久
10ンテナが配置されている。また、プラズマ源の下|B|(T)
5

流側（z = 4.0 mm）には 4 つの穴を持つ厚さ 0.6 mm
mm までを電子引き出しの解析領域とした。オリ
フィス面を含む金属面のポテンシャルはゼロとし、電子引き出しのために z = 10 mm でのポテンシャ
ルは 20 V とした。

y(mm)のオリフィス面を備えており、この面から z = 10

0.6
0.55
0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05

ECR

0

layer

-5

本研究で用いる粒子計算法は Particle-in-Cell(PIC) 法であり、荷電粒子と中性粒子との衝突をMonte Carlo collision (MCC) 法により考慮している。計算セルは各方向とも 0.1 mm の等間隔格子であり、計算は以下の仮定の下で行った。(i) 荷電

-10
0

5

z(mm)

10

図 3: 磁束密度分布 [7](JSASS-2016-4687 Ⓒ日本航空宇宙学会)

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 39

磁石による静磁場Bst はANSYS を用いて解析して行う必要が生じるが、複雑な境界条件を判定する

いる（図 3）。マイクロ波周波数 4.2 GHz はイオン

必要なく並列化することが可能になる。

プラズマ周波数と比べて非常に大きいため、計算

次に、ループ内の依存関係を回避するために

時間短縮の都合イオンの運動では EEM を無視し、RedBlack 法を採用することにした。また、現状扱

静電場 EES はマイクロ波 1 周期で平均化した値をっている粒子数から 1000 倍以上粒子を増やさな用いて、時間ステップ∆ti = 2.38×10s（マイクロ

い限り、システム B の 1 ノードで実行可能なレベ

波の 1 周期）で更新している。本研究ではプラズルであり、かつ、ポアソン方程式を解くサブルー

マ吸収電力 P0 を入力パラメータとしており、イオチン内でのメモリ使用量も 100 MB 程度と小さく、ンと電子の運動方程式を解く前後での粒子の運動

格子間隔や粒子数に依存しないことからプロセス

エネルギーの差から得られる電力Pabs がP0 を満た

並列ではなくスレッド並列を採用した。

−10

すように E0 を調節している。上記の計算を定常解が得られるまで繰り返し行う。

これらの対策の結果、ポアソン方程式を解くサ
ブルーチンにおいて、オリジナル（1n1p1t）に対計算条件は、マイクロ波周波数 f = 4.2 GHz、吸してスレッド並列化（1n1p36t）の適用により、43.7

収電力 P0 = 2.0 W である。電子とイオンの初期温倍の高速化が実現できた。その他、細かいチュー

度をそれぞれ 2.0 eV、0.05 eV とし、いずれも初期ニングを施すことで最終的には約 53 倍の高速化

密度 1×10 m で一様に分布させた。その他、モを達成している。なお、SOR 法とは演算順序が異デルの詳細は過去の文献に委ねる [3-6]。

なるため、計算結果の完全一致は得られないが、

16−3

両解法における計算結果の差異は実用上問題ない

3 プログラム高度化内容と性能評価レベルであった。
ポアソン方程式を解くサブルーチンの高速化が

図 2 のフローチャートが示すように、本プログ実現出来た結果、全計算時間に占める割合は 50％

ラムで大部分の計算時間を占めるのは個々のプラ程度まで削減できた。その結果、計算格子内のプズマ粒子の運動変化を求めるループとなる。

ラズマ粒子の位置から計算格子の各コーナーのグ

平成 29 年度におけるプロファイル計測の結果

リッドへ電荷などを割り当てる allocation のサブによると、中でもポアソン方程式を解いてプラズ

ルーチンが計算時間全体に占める割合が大幅に増

マポテンシャルを求めるサブルーチンが実に全体加し 30%程度を占めるに至った。この allocationの 90%を占めている事が明らかになった。これはの計算時間削減のため、粒子番号によるデータ並

高度化前のプログラムにおいて、図 1 の計算領域

列化、omp atomic によるループ内排他処理、デーが示すようにプラズマが分布する領域にアンテナ

タ移動量削減、スレッド同期削減など施すことにやオリフィスなどが存在するため、遅くても最もより 15%まで低下した。

単純でバグが潜みにくいことを優先し SOR 法を

その他のサブルーチンについても細かいチュー

採用した事に起因している。また、最内ループにニングを行った結果、プログラム全体としては、

IF 文を用いて計算領域内のプラズマの有無を判オリジナルのコードに対してスレッド並列化

断している点も要因として挙げられる。

（1n1p36t）
で約 28.6 倍まで高速化することができ

ポアソン方程式の高速化において、SOR 法ではた。

配列の FLOW 依存（前方依存）、ANTI 依存（後この状況において、今回、更なる高度化を図る

方依存）が存在するためコンパイラによるベクト

ため、
領域分割法と粒子ソート法を組み合わせた。

ル化およびスレッド並列が効かない状況である。

具体的には、前述の allocation のサブルーチンにおこの状況を改善するため、まずは最内ループ内にいて、粒子数による並列化から、格子空間領域の存在する IF 文において TRUE と FALSE のいずれ

分割による並列化への対応となる。そのため、分の場合でも演算式を共通化するマスク方式を採用

割した各領域に含まれる粒子数と粒子番号の配列

して IF 文を除去した。これにより、本来ポアソンを生成する必要があり、新たなサブルーチン div1d方程式を解かずに済むアンテナやオリフィスプレを作成し、allocation 呼び出しの前に div1d を呼びート内部についてもポアソン方程式を解く演算を出す方式とした。その結果、allocation そのものの40 − Academic Center for Computing and Media Studies, Kyoto University図 4: 負イオン密度および電子電流密度の周方向分布 [6] (ⒸAIP Publishing)
速度はメモリアクセスの改善から 3 倍程度の速度

K.: Engineering Model of the Minature Ion Propulsion

向上が得られた一方、前処理の div1d における経System for the Nano-satellite: HODOYOSHI-4, Trans.

過時間が加算されるため、全体としては約 1 割程

JSASS Aerospace Tech. Japan, 12 (2014), pp. Tb_19–

度の性能低下に至った。

Tb_24.

次に、上記に加えて、粒子ソートを行うことに[3] Takao, Y., Koizumi, H., Komurasaki, K., Eriguchi, K., andより、div1d におけるデータ局所性による性能向Ono, K.: Three-dimensional particle-in-cell simulation of a

上を図った。結果として、div1d は 1.5 倍程度の速

miniature plasma source for a microwave discharge ion

度向上となったが、全体としては高度化前のコー

thruster, Plasma Sources Sci. Technol., 23 (2014), pp.

ドと同等程度の性能になった。

064004-1–11.

結果として、今回の高度化においては、平成 29[4] Takao, Y., Hiramoto, K., Nakagawa, Y., Kasagi, Y.,

年度の高度化を凌駕する性能は得られなかった。

Koizumi, H., and Komurasaki, K.: Electron extraction

しかしながら、
平成 29 年度の高度化により得られmechanisms of a micro-ECR neutralizer, Jpn. J. Appl.

たコードを水に適用することにより、負イオンを

Phys., 55 (2016), pp. 07LD09-1–5.

含む放電に特有な不安定性などが分かっており[5] Hiramoto, K., Nakagawa, Y., Koizumi, H., and Takao, Y.:（図 4）[6]、今後このコードを利用して、電子引

Effects of E × B drift on electron transport across the

き出し効率の向上を目指した解析を進めていく予

magnetic field in a miniature microwave discharge

定である。

neutralizer, Phys. Plasmas, 24 (2017), pp. 064504-1–5.
[6] Nakamura, K., Koizumi, H., Nakano, M., and Takao, Y.:
Effects of negative ions on discharge characteristics of

謝辞プログラム高度化の実施にあたって、HPE のwater plasma source for a miniature microwave discharge

伊藤 泰善 様には大変ご尽力いただきました。ま

ion thruster, Phys. Plasmas, 26 (2019), pp. 043508-1–11.

た本研究は、JSPS 科研費 JP16H06370 の助成[7] 平本謙太, 中川悠一, 小泉宏之, 小紫公也, 鷹尾祥を受けて実施されました。ここに感謝の意を表し

典：小型マイクロ波放電式中和器のプラズマ源におます。

ける E×B ドリフトの効果、第 58 回宇宙科学技術連
合講演会講演集, 2016, JSASS-2016-4687。
[8] Takao, Y., Eriguchi, K., Ono, K., Sugita, Y., Koizumi, H.,

参考文献

and Komurasaki, K.: A Validation Study of a 3D PIC[1] Nano/Microsatellite Market Forecast 2020 (SpaceWorks

Model for a Miniature Microwave Discharge Ion Thruster,
50th AIAA/ASME/SAE/ASEE Joint Propulsion Conf.,

Enterprises, Inc., 2020).
[2] Koizumi, H., Komurasaki, K., Aoyama, J., and Yamaguchi,

2014, AIAA-2014-3829.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 41

巨大津波遡上時の木造家屋の瓦礫生成過程シミュレーションの高度化
浅井光輝
九州大学大学院工学研究院社会基盤部門

1 はじめに線形チモシェンコはり要素だけで分割し，数値積分
点をシフトすることではりの崩壊現象を適切に表2011 年東北地方太平洋沖地震で生じた大規模な現する手法であり，Updated Lagrangian の定式化

津波により，沿岸域の木造家屋の多くは瓦礫となっにより大変形までを解析可能な有限変形弾塑性解た．津波と瓦礫が混在して遡上することで，人的・析方法である．また破断後は，要素の断面力を解放

経済的損失が拡大した．瓦礫の総量は約 2000 万ト

することで表現している．降伏関数として次式を採ンにも及び，災害復旧時にはその膨大な瓦礫の処理用した．本研究では Mx，My，N はそれぞれ x，yに時間を有し，復旧・復興に遅れが生じた．今後，軸回りの曲げモーメント，軸力とし，次式に示す降南海トラフなど，同規模の地震・津波被害が生じる

伏関数を用いた．ここで，右下添え字の“0”は全断危険性が高いことが指摘されている．この規模の津面塑性値であることを示す．

波からの構造物等の被害をゼロにすることは非現実的であることから，人的被害をゼロにするために

𝑀𝑀

2

𝑀𝑀

2

𝑁𝑁 2

𝑓𝑓𝑦𝑦 = (𝑀𝑀 𝑥𝑥 ) + (𝑀𝑀 𝑦𝑦 ) + ( ) = 1   (1)
𝑁𝑁
𝑥𝑥0

𝑦𝑦0

0

津波遡上現象の予測精度を向上させると同時に，瓦

従来の ASI-Gauss 法は，並列計算に対応してお礫の総量とその拡散状況を事前評価することで災らず，また要素数が数千万オーダーのモデル解析を害後の復旧活動に備えておく必要がある．

想定していないため，計算効率の低いアルゴリズム

以上の背景から地震および津波に対する木造家が存在していた．都市全域の大規模地震応答解析を屋倒壊解析を都市全体規模で実施するためのツーするにあたって，ASI-Gauss 法の並列化およびアルルの開発を行った．この解析のためには，GIS デー

ゴリズムの改良が必要となった．

タ等から都市全域の建物のモデル化を自動化し，まASI-Gauss 法を計算並列化するため，ソースコー

た広域な領域に含まれる多数の家屋の非線形解析

ドのスレッド並列化（以下，OpenMP と記述）を段

をロバストに実施できるツールが望まれる．そこで階的に行った．計算コストが高く非効率な場所を計

本研究では建物の倒壊解析に実績のある ASI-

測し，連立一次方程式ソルバーと接触探査アルゴリ

Gauss 法を解析手法として選択し，大規模解析を実ズム，全体剛性行列構築アルゴリズムの改良を行っ施するための改良を行った．また木造家屋のモデル

た．ソルバーについては剛性行列の格納法として化には建築分野の設計基準に準拠した強度を付与

CRS 法（Compressed Row Storage）を採用し，スすることで建築年代に整合した倒壊解析を実施し，レッド並列効率の高い共役勾配法ライブリーを移地震動による都市全体に存在する木造家屋の崩壊

植した．また接触探査には，特に粒子系の解析にお予測を可能とした．

ける近傍粒子探索法で用いられているバックグラ
ウンドセルを使ったセル検索方法を応用すること

2 ASI-Gauss 法とその高速化で，接触探索すべき要素を限定することで大幅な高速化を実現した．

本研究では構造解析手法として ASI-Gauss 法を

前述のソルバーOpenMP 化と，接触探査アルゴリ

採用した．ASI-Gauss 法は筑波大学磯部らによってズムを改良した結果，十分な高速化を達成した．性

開発された FEM 解析手法である．1 部材を 2 つの能評価に使用したモデルは，図－1 に示す要素数：

42 − Academic Center for Computing and Media Studies, Kyoto University 
図－1 モデル小            図－2 モデル大
（要素数：43934，節点数：30370）   （要素数： 217 万，節点数：155 万）
43934，節点数：30370 のモデルである．このモデしていることが問題であった．そこで上記アルゴリルを用いて，
1000 ステップの計算に要する時間を測

ズムの見直しを行い，約 9 秒で計算できるように改定した．OpenMP 化コードの解析環境は，京都大学良した（約 1830 倍の高速化）．

「Camphor 2」，1 ノード，64 スレッドである（京以上，ソースコードの OpenMP 化，接触探査・全コンピュータ停止後、別途自身で計算機環境は調達

体剛性行列構築アルゴリズムの改良を行い，ASI-した）．その結果，従来のコードでは 33 時間 32 分

Gauss 法は大規模崩壊解析が実施可能となった．

要した計算時間が，OpenMP 化により 13 時間 37分になり，約 2.5 倍の高速化を達成した．特にソルバーのサブルーチン単体に着目すると，35 倍高速化

3 木造家屋モデル化指針（並列化効率：57％）した．この結果から，ソルバ

本研究で対象とする木造家屋の軒数は非常に多ーは大幅に高速化したものの，他の計算効率が低い

く，建物ごとの詳細なモデル化は不可能である．そサブルーチンが，解析時間の多くを占めていることのため一軒一軒の詳細なモデル化はせず，簡易的にが 2.5 倍の高速化に留まっている要因であった．

モデル化し，建築年代ごとに適用された建築基準法

次に接触探査アルゴリズムの改良を実施した．従に基づいた耐力を保有させる方針とした．図－3 に来のコードでは，一要素が，その他すべての要素に示すように，木造家屋をはり要素によりモデル化を

対し，接触探査のための計算を行っていた．無駄な行った．厳密なモデル化はせず，層せん断力と層間計算を削減するために，解析領域をセル単位で分割

変形角の関係から，構造物全体としての耐力を設定

し，接触探査の候補を自身のセルと近傍のセルのみした．木造家屋を柱，梁，壁，屋根（床）の 4 つのに限定することで接触探査の高速化を実現した．そ

要素に分類し，壁と屋根（床）要素はブレースに置の結果，図－1 と同様のモデル・解析条件において，換，屋根（床）は剛体床と仮定した．モデルの重量

接触探査アルゴリズム改良後のコードでは，解析時

は，「2012 年改訂版 木造住宅の耐震診断と補強方間が 13 時間 37 分から 2 時間 11 分になり，さらに法」より代表的な建築年代ごとの床面積当たりの重6 倍近く高速化することができた．
（ソルバーと併せ量を参照することで決定した．屋根（床）の重量をると約 15 倍の高速化を実現．
）

梁と屋根要素に，外壁，内壁の重量を柱と壁要素に以上の高速化は，非線形解析までに留まっており，分配した．部材の降伏は，層間変形角がある閾値に

ASI-Gauss 法一番の特徴である崩壊後の解析まで達したとき，降伏と判断した．例えば，1959 年基準

実施すると，急激に解析速度の低下を招いていた．

に準拠する図－4 の場合，壁を構成する要素は層間

図－1 と同様のモデルでは，全体剛性行列構築に 7変形角が 1/120 rad，柱を構成する要素は 1/30 rad

秒要したのに対し，図－2 に示す要素数：217 万，に達した時に降伏させ，要素端部に塑性ヒンジを形

節点数 155 万オーダーのモデルでは約 5 時間要して成させた．スギの材料特性を用い，降伏後は完全弾

いた．接触計算では，ギャップ要素を追加し，その塑性体と仮定した．以上のように，木造家屋をはり要素を介して接触要素間の力を伝えている．そして要素にてモデル化し，建築基準法に準拠した耐力を

全体剛性行列を再び構築している．そのため接触判持たせる方針とした．

定・接触解除判定されるたびに，ギャップ要素が増減し，約 5 時間かけて全体剛性行列を一から作り直京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 43

1,514,782（絶対偏差：1.1 %）であった．領域分割の結果，各領域の平均解析時間は 28 時間 28 分
（Camphor 2：20 ノード，68 スレッド，1000 ステップ）に抑えることができ，都市全体の弾塑性地震
応答解析（変位入力）が数日の内に完了できるようになった．図－7 に 2 区画を可視化した結果を示す．
カラーコンターは水平変位を意味する．背の高い建物ほど大きく揺れていることがわかる．なお，解析図－3 木造家屋モデル

では要素の破断，接触はなく，都市全域モデルには，



個々のモデル耐力が未設定，壁ブレースの未設定であること留意しておく．

層せん断力係数

0.25

1/30

0.207 0.2
1/120

0.15

0.115

0.1
0.05

0

0

0.02

0.04

0.06

0.08

0.1

0.12

層間変形角(rad)図－4 層せん断力と層間変形角の関係（ 年基準）
図－5 高知市モデル

4 都市全域モデルの自動化処理
都市全域のモデル化には，地理情報システム（GIS）
を用いた．GIS より得られる建物の立体位置情報と
数値標高データから，都市全域の建物モデルを構築できる．また同データから，津波遡上解析に必要な
地表面モデルの作成も可能である．本研究で対象とする都市は，高知県高知市とした．ASI-Gauss 法を用いた都市全域のモデル化及び地震応答解析を達成後は，SPH 法と連成解析することで，地震・津波
に対する都市全域の木造家屋倒壊解析を実施する方針である．
高知市モデルの概要を次に示す．領域：7 km×10

図－6 高知市モデルの領域分割

km，建物数：82,916，要素数：43,411,162，節点数：
30,295,639，データサイズ：9.5 GB（図－5）．単一の計算機でも解析できないことはないが，数日の間に解析を終了させるためには，複数の計算機による同時処理を行うことにした．このため，複数の計算機に計算負荷が一致するように領域分割を行った．
その例として，高知市全域モデルを 20 区画に分割した例を示す（図－6）．このモデルの平均要素数は
2,170,558 （絶対偏差：0.3 % ），平均節点数は44 − Academic Center for Computing and Media Studies, Kyoto University

図7 地震応答解析結果（2 区画分のみ）

5 都市全域モデルによる木造家屋倒壊解析例

6 おわりに地震動による都市全域の木造家屋倒壊解析を実

都市全域の木造家屋倒壊解析の実現可能性を確認するため，1 区画のモデルを用いて，1000 ステッ

施するために，ASI-Gauss 法を高速化させ，また建築基準に整合した木造家屋のモデル化方法を整備しプ，実時間 20 秒の計算を実施した．木造家屋相当

た．都市全体の木造家屋を効率よくモデル化するたの節点数 1000 未満のモデルに対して，弱い強度をめ，GIS で登録された情報から建物モデルを自動生設定し，地震動は変位入力した．要素の破断，接触，成した．また構造解析手法の OpenMP 化，アルゴリ

接触解除，地面との接地を考慮した．その結果，56.6

ズム改良をすることで，上記の解析を現実的な時間時間で計算を終えることを確認した（結果出力を除内で実施可能にした．今後は，まず過去の地震被災

けば 42.2 時間）．図－8 に 1 区画モデルを拡大した事例との比較による木造家屋倒壊解析の妥当性を結果を示す（0 秒，7 秒）．カラーコンターは鉛直変確認し，さらには ASI-Gauss 法による家屋の解析と位を意味する．要素の破断が 210809 回，要素間の接触が 5524 回発生しており，破断後に接地したと判定された回数は 29311 回あり，全体剛性行列の形

津波遡上解析法を連成することで，地震だけでなく津波に対する都市全域の木造家屋倒壊予測技術へ状が 37854 回変わる計算であった．全体剛性行列構の発展が望まれる．
 本研究は京都大学学術情報メディアセンター

築アルゴリズムの改良により，現実的な時間内で計2019 年度スーパーコンピュータ共同研究制度・プロ

算でき都市規模の崩壊解析が可能になった． 以上，

グラム高度化に採択いただき，特に連立一次方程式

都市全域の木造家屋倒壊解析の実現可能性検討を

ソルバーの高速化を実施していただいた．ここに記行った．各種アルゴリズムの改良により，大規模崩して感謝の意を表す．

壊解析が数日のうちに終了することを確認できた．
今後，壁ブレースを追加した都市モデルに木造家屋の物性値・断面形状を適用し，年代ごとの倒壊過程の評価へと発展させる予定である．

図－8 木造家屋倒壊解析（初期モデルと地震動流力後 7 秒後）

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 45

サービスの記録・報告システム A 運転状況 ( 2019 年 10 月 ～ 2020 年 3 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時

サービス再開日時

システムダウン障害発生状況
障害発生日時

サービス再開日時

2019/11/11

2019/11/11

保守時間[h]2019/10/14

9:00

2019/10/17

16:55

81.00

2019/12/05

9:00

2019/12/05

11:30

6.00

2020/03/30

9:00

2020/03/31

24:00

39.00

14:20

ダウン時間[h]20:08２）サービス状況
サービス
時間[h]

ジョブ
処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

ノード
利用率

10 月

662.4036,918

182,30257,866,786

46,713,480

1745.7

76.0

%

11 月

715.2036,386

222,566

69,326,480

57,355,983

1798.2

81.6

%

12 月

744.0050,889285,720

73,444,981

57,553,181

1795.6

81.7

%

1月

744.0073,141

284,384

76,523,814

54,088,568

1795.5

83.0

%

2月

696.0098,596230,00767,725,077

49,308,813

1792.6

79.3

%

3月

705.6050,476250,096

71,487,519

56,244,208

1798.780.6%

計4,267.20

346,406

1,455,075

416,374,657

321,264,233

1787.780.4%

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

10月

11月

12月

1月

2月

3月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

46 − Academic Center for Computing and Media Studies, Kyoto University

5.80

システム B 運転状況 ( 2019 年 10 月 ～ 2020 年 3 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時

サービス再開日時

システムダウン障害発生状況
障害発生日時

サービス再開日時

2019/11/11

2019/11/11

保守時間[h]2019/10/14

9:00

2019/10/17

16:55

81.00

2019/12/05

9:00

2019/12/05

11:30

6.00

2020/03/30

9:00

2020/03/31

24:00

39.00

14:20

ダウン時間[h]

20:08

5.80

２）サービス状況
サービス
時間[h]

ジョブ
処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

ノード
利用率

10 月

662.40119,090

657,152

11,986,043

9,221,575

779.0

70.0

%

11 月

715.20138,853765,685

14,054,172

11,550,943

807.9

73.2

%

12 月

744.00153,638

791,775

14,589,627

12,012,627

807.6

72.3

%

1月

744.00336,625

893,723

14,805,263

12,319,297

808.0

72.2

%

2月

696.00165,011753,756

13,723,839

11,204,145

808.0

70.8

%

3月

705.60113,642624,961

13,506,733

10,670,110

808.0

69.1

%

計

4,267.20

1,026,859

4,487,052

82,665,677

66,978,697

803.1

71.3

%

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

10月

11月

12月

1月

2月

3月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 47

システム C 運転状況 ( 2019 年 10 月 ～ 2020 年 3 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時

サービス再開日時

システムダウン障害発生状況
障害発生日時

サービス再開日時

2019/11/11

2019/11/11

保守時間[h]2019/10/14

9:00

2019/10/17

16:55

81.00

2019/12/05

9:00

2019/12/05

11:30

6.00

2020/03/30

9:00

2020/03/31

24:00

39.00

14:20

ダウン時間[h]20:08２）サービス状況
サービス
時間[h]

ジョブ
処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

ノード
利用率

10 月

662.403,705

32,427

376,447

256,345

15.7

72.0

%

11 月

715.2020,80825,688508,109478,321

16.0

83.3

%

12 月

744.004,314

31,717526,731550,943

15.9

84.7

%

1月

744.004,613

37,859516,085573,616

16.0

72.6

%

2月

696.003,134

22,247

378,711

251,428

16.0

70.6

%

3月

705.602,300

28,396

304,021236,200

15.9

46.8

%

計

4,267.20

38,874

178,334

2,610,104

2,346,853

15.9

71.7

%

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

10月

11月

12月

1月

2月

3月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

48 − Academic Center for Computing and Media Studies, Kyoto University

5.80京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 49

資料

別表１ スーパーコンピュータシステム
提供サービス
利用負担額
コース

タイプ

セット

エントリ

タイプA
タイプB
タイプC

基本

パーソナル

タイプA1
タイプA2
タイプA3
タイプB1
グループ
タイプB2
タイプB3
タイプC1
タイプC2
タイプA
大規模ジョブ タイプB
タイプC
専用クラスタ

-

基本
基本
基本
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位

ストレージ容量追加
ライセンスサービス

システム バッチ

12,600 円/年
100,000 円/年
100,000 円/年
100,000 円/年
200,000 円/年
100,000 円/年
240,000 円/年
60,000 円/年
600,000 円/年
300,000 円/年
210,000 円/年
105,000 円/年
252,000 円/年
63,000 円/年
630,000 円/年
315,000 円/年
130,000 円/年
130,000 円/年
156,000 円/年
78,000 円/年
20,000 円/週(7日）
10,000 円/週(7日）
21,000 円/週(7日）
10,500 円/週(7日）
13,000 円/週(7日）
6,500 円/週(7日）
630,000 円/年
315,000 円/年
10,000 円/年
20,000 円/年

システム資源

経過時間 ストレージ
無料
（時間）
（TB)
利用者数

最大1ノード相当((36コア、128GBメモリ)×1)
1
最大4ノード相当((68コア、16+96GBメモリ)×4)
最大4ノード相当((36コア、128GBメモリ)×4)
168
B
最大1ノード相当((72コア、3072GBメモリ)×1)
C
4ノード((68コア、16+96GBメモリ)×4)
優先
2ノード((68コア、16+96GBメモリ)×2)
8ノード((68コア、16+96GBメモリ)×8)A
準優先
2ノード((68コア、16+96GBメモリ)×2)
8ノード((68コア、16+96GBメモリ)×8)
占有
4ノード((68コア、16+96GBメモリ)×4)
4ノード((36コア、128GBメモリ)×4)
優先
2ノード((36コア、128GBメモリ)×2)
336
8ノード((36コア、128GBメモリ)×8)
B
準優先
2ノード((36コア、128GBメモリ)×2)
8ノード((36コア、128GBメモリ)×8)
占有
4ノード((36コア、128GBメモリ)×4)
1ノード((72コア、3072GBメモリ)×1)
優先
1ノード((72コア、3072GBメモリ)×1)
C
2ノード((72コア、3072GBメモリ)×2)
準優先
1ノード((72コア、3072GBメモリ)×1)
8ノード((68コア、16+96GBメモリ)×8)A
占有
4ノード((68コア、16+96GBメモリ)×4)
8ノード((36コア、128GBメモリ)×8)B
占有
168
4ノード((36コア、128GBメモリ)×4)
2ノード((72コア、3072GBメモリ)×2)
C
占有
1ノード((72コア、3072GBメモリ)×1)
8ノード((36コア、128GBメモリ)×8)
B
4ノード((36コア、128GBメモリ)×4)
ストレージ容量10TBの追加につき可視化ソフト（AVS,ENVI/IDL)およびプリポストウェアの1ライセンスにつき
B
A

共有
共有
共有
共有

0.2
3.0
3.0
3.0
24.0
12.0
28.8
7.2
48.0
24.0
24.0
12.0
28.8
7.2
48.0
24.0
24.0
24.0
28.8
14.4

-

-

-

8
4
16
4
16
8
8
4
16
4
16
8
8
8
16
8

48.0
24.0

16
8

備考
1． 利用負担額は、年度単位で算定している。また、総額表示である。パーソナルコース、グループコース又は専用クラスタコースを、
年度途中から利用を開始する場合及び年度途中で利用を終了する場合の利用負担額は、上記表中の利用負担額を12で除した後、利用月数を乗じて算出するものとし、100円未満に端数が出た場合は、10円単位を四捨五入するものとする。
なお、月途中から利用を開始する場合及び月途中で利用を終了する場合は、それぞれ1月の利用とする。
2． 大型計算機システムの全ての利用者は、上記表のサービスの他、次のサービスを受けることができる。
1） 大判プリンタサービス
2） その他、大型計算機システムが提供するサービス、機器の利用
3． 上記表の大規模ジョブコース、ストレージ容量追加、ライセンスサービスの申請には、スーパーコンピュータシステムの利用者であることが必要である。
4． 「共有」： 当該カテゴリのユーザ間で一定の計算資源を共有するベストエフォートのスケジューリングを行う。
「準優先」： 定常稼働状況において記載値（以上）の計算資源が確保されるように優先スケジューリングを行う。
また、稼働状況によらず記載値の1/4の計算資源が確保されることを保証する。
「優先」： 定常稼働状況において記載値（以上）の計算資源が確保されるように優先スケジューリングを行う。
また、稼働状況によらず記載値の1/2の計算資源が確保されることを保証する。
「占有」： 稼働状況によらず記載値の計算資源が確保されることを保証する。
5． ストレージ容量はバックアップ領域（最大で総容量の1/2）を含む。
6． グループコース及び専用クラスタコースの利用者番号は利用者あたり年額5,000円を負担することで追加できる。
7 ． 機関・部局定額制度
他機関又は学内における部局（『国立大学法人京都大学の組織に関する規程』第3章第2節から第11節で定める組織をいう。）の組織が、その組織単位でグループコースサービスを利用申請する場合の利用負担額は、別表１に規定する
1.5倍の額とする。なお、利用負担額が年額150万円未満の場合は100人、年額150万円を超える場合は、150万円毎に100人までの利用者を認める。ストレージは、1.5倍の容量とする。
8 ． スパコン連携サービス
学術情報メディアセンターのスーパーコンピュータシステムと密な連携により、学内における部局の組織が計算サーバ等を設置する
場合、下記の負担額を支払うものとする。
冷却方式
水冷
空冷

利用負担額
9,800 円/月
11,500 円/月

利用負担額算定単位
水冷冷却方式の計算サーバ等の定格電力 1kWにつき
空冷冷却方式の計算サーバ等の定格電力 1kWにつき50 − Academic Center for Computing and Media Studies, Kyoto University

別表２ 汎用コンピュータシステム
区

分

仮想サーバ
ホスティングサービス

利用負担額

単

位

36,000円／年 1仮想サーバにつき

備考
1． 利用負担額は、総額表示である。
2． 上記表の仮想サーバホスティングサービスを利用するには、
スーパーコンピュータシステムの利用者であること。
3． 1仮想サーバに割当てるシステム資源は、CPU：2コア、メモリ：4GB、ディスク：100GBである。
4． 仮想サーバホスティングサービスにおいて、下記の負担額を支払うことによりCPU、メモリ、ディスクを増量することができる。
区

分

利用負担額

単

位

CPU増量

3,000円／年 2コアにつき(最大8コアまで)

メモリ増量

3,000円／年 4GBにつき(最大64GBまで)

ディスク増量6,000円／年 100GBにつき(最大1,000GBまで)

5． 利用負担額は、当該年度(4月から翌年3月まで)の利用に対して年額として算定するが、年度途中から利用を開始する場合には月数に応じて減額する。

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 51

別表３ スーパーコンピュータシステム
システム

A

B

C

システム資源

経過時間 ストレージ
無料
（時間）
（TB)
利用者数

利用負担額

8ノード(68コア、16+96GBメモリ)×8)

336

28.8

16

960,000 円/年

12ノード(68コア、16+96GBメモリ)×12)

336

43.2

24

1,440,000 円/年

16ノード(68コア、16+96GBメモリ)×16)

336

57.6

32

1,920,000 円/年

8ノード(36コア、128GBメモリ)×8)

336

28.8

16

1,008,000 円/年

12ノード(36コア、128GBメモリ)×12)

336

43.2

24

1,512,000 円/年

16ノード(36コア、128GBメモリ)×16)

336

57.6

32

2,016,000 円/年

2ノード(72コア、3072GBメモリ)×2)

336

28.8

16

624,000 円/年

3ノード(72コア、3072GBメモリ)×3)

336

43.2

24

936,000 円/年

4ノード(72コア、3072GBメモリ)×4)

336

57.6

32

1,248,000 円/年

備考
1． 利用負担額は、年度単位で算定している。また、総額表示である。パーソナルコース、グループコース又は専用クラスタコースを、
年度途中から利用を開始する場合及び年度途中で利用を終了する場合の利用負担額は、上記表中の利用負担額を12で除した後、利用月数を乗じて算出するものとし、100円未満に端数が出た場合は、10円単位を四捨五入するものとする。
なお、月途中から利用を開始する場合及び月途中で利用を終了する場合は、それぞれ1月の利用とする。
2． ストレージ容量はバックアップ領域（最大で総容量の1/2）を含む。

52 − Academic Center for Computing and Media Studies, Kyoto University

全国共同利用版広報・Vol.18(2019)総目次［巻頭言］
Vol.18，No.1 号の発刊にあたって ................................................................................................................1-1
Vol.18，No.2 号の発刊にあたって ................................................................................................................2-1［スーパーコンピュータ共同研究制度（若手・女性研究者奨励枠）研究報告］
機械攪拌操作時気泡巻き込みに対する数値解析 ..........................................................................................1-2
大規模溶解炉における溶融アルミニウム機械撹拌時の流動解明 ...............................................................1-4
界面分光の分子動力学シミュレーション......................................................................................................1-6
高効率有機系太陽電池の実現に向けた光機能性分子の構造と電子物性の相関解明 ................................1-8
無重力下での高プラントル数流体における温度差マランゴニ効果に起因する液柱内対流場の二次不安定性 .............................................................1-10
Adjoint sensitivity 解析を用いたパワーデバイス用半導体製造装置の最適設計のための数値解析 ....1-13
ウマの個体間に作用する力の解明に向けた数値シミュレーション .........................................................1-15
界面の摩擦接触を考慮した損傷モデルによる鉄筋コンクリートの 3 次元破壊シミュレーション.......1-17
N 結合型糖鎖修飾によるタンパク質の機能制御の関連性 ........................................................................1-19
都市構造物の幾何的特徴がもたらす大気乱流の空間スケールへの影響..................................................1-21
Numerical simulation of InGaSb crystal growthunder micro-gravity onboard the International Space Station .............................................................1-23
3 次元流体変数の予測 ...................................................................................................................................1-25
Numerical simulation of deepwater oil blowout--Turbulent jets and droplet size distribution .......1-27
分子動力学計算による膜貫通型ペプチドとリン脂質二重膜の相互作用ダイナミクス ..........................1-29［プログラム高度化共同研究報告］
分子シミュレーションによるヌクレオソーム構造変化の網羅的探索 .....................................................1-31
巨大津波遡上時の木造家屋の瓦礫生成過程シミュレーション.................................................................1-33
異方性弾性波動問題に対する演算子積分時間領域境界要素法の高性能化..............................................1-37
飽和土の大規模変形・流動計算を目的とした固液混合 MPM の開発.....................................................1-41
強い音響異方性を有する CFRP に対する開口合成法の高速実行............................................................1-45［ARM プロセッサの性能評価］
Marvell ThunderX2 評価機の紹介 ...............................................................................................................2-2
ThunderX2 評価機を利用した MHD シミュレーションコードの性能評価.............................................2-4［サービスの記録・報告］
スーパーコンピュータシステムの稼働状況.................................................................................... 1-50，2-10
センター利用による研究成果（平成 30 年度）..........................................................................................2-13［資料］
大型計算機システム利用負担金 別表...............................................................................................1-54,2-18
全国共同利用版広報・Vol.17(2018)総目次 .................................................................................................1-57

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 53

サービス利用のための資料一覧 ..........................................................................................................1-59,2-21［奥付］
奥付 ........................................................................................................................................................1-60,2-22

54 − Academic Center for Computing and Media Studies, Kyoto University

―1.

サービス利用のための資料一覧―

スーパーコンピュータシステム・ホスト一覧
システム A：camphor.kudpc.kyoto-u.ac.jp



システム B・C：laurel. kudpc.kyoto-u.ac.jp


システム B（SAS 利用時）
：sas.kudpc.kyoto-u.ac.jp※ ホストへの接続は SSH(Secure SHell) 鍵認証のみ、パスワード認証は不可

2.

問い合わせ先 & リンク集


情報環境機構のホームページ
http://www.iimc.kyoto-u.ac.jp/



学術情報メディアセンターのホームページ
http://www.media.kyoto-u.ac.jp/



スーパーコンピュータシステムに関する問い合わせ先

利用申請などに関する問い合わせ先【情報環境支援センター】
E-mail : zenkoku-kyo@media.kyoto-u.ac.jp / Tel : 075-753-7424
URL: http://www.iimc.kyoto-u.ac.jp/ja/services/comp/

システムの利用など技術的な問い合わせ先【スーパーコンピューティング掛】
E-mail : consult@kudpc.kyoto-u.ac.jp / Tel : 075-753-7426
URL: http://www.iimc.kyoto-u.ac.jp/ja/services/comp/contact.html

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.19 No.1（ Oct. 2020）− 55

奥付京都大学学術情報メディアセンター全国共同利用版広報 Vol．19，No．1
２
０２０年
３０日xx
発行
2020
年 １０月
10 月
日発行
広報編集部会
編集者 京都大学学術情報メディアセンター

深沢 圭一郎（部会長）

全国共同利用版広報編集部会

平石 拓（副部会長）

発行者 〒606-8501 京都市左京区吉田本町
京都大学学術情報メディアセンター

熊谷 真由美尾形 幸亮

Academic Center for Computing and Media Studies
Kyoto University
Tel. 075-753-7407
http://www.media.kyoto-u.ac.jp/
印刷所 〒616-8102 京都市右京区太秦森ヶ東町 21-10
株式会社エヌジーピー

56 − Academic Center for Computing and Media Studies, Kyoto University

表紙デザイン：中山 豊（中山商店）