1. No category

Simulation Engineering Solution’s
No . １ ０
January 2014
【
【 第
】
拶 】
挨拶
ご挨
のご
刊の
発刊
号発
第 1100 号
“
”
て”
って
返っ
り返
振り
を振
号を
∼99 号
号∼
１号
“１
【
】
文 】
論文
術論
技術
【 技
「
１
た
した
携し
連携
を連
ムを
ラム
グラ
ログ
プロ
種プ
異種
る異
ける
おけ
にお
野に
分野
造分
構造
鋼構
築鋼
建築
「建
．
１．
統
」
て」
いて
つい
につ
法に
手法
析手
解析
化解
合化
統合
「
２
」
て」
いて
つい
につ
析に
解析
成解
連成
M連
た CCFFDD--DDEEM
いた
用い
を用
Mを
ス CCFFDDEEM
ース
ソー
ンソ
プン
ープ
オー
「オ
．
２．
３
「
３．
．
「自
自動
動メ
メッ
ッシ
シュ
ュ生
生成
成機
機能
」
ム」
能を
テム
ステ
シス
析シ
解析
展解
進展
裂進
亀裂
元亀
次元
る 33 次
する
有す
を有
ＮＳプラント設計株式会社
NS Plant Designing Corporation
Contents
○巻頭言
「ＮＰＤ数値解析技報第１０号発刊のご挨拶」
・・・・・・・・・・・・・・・・・・・・・・1
○技術論文
「 建 築 鋼 構 造 分 野 に お け る 異 種 プ ロ グ ラ ム を 連 携 し た 統 合 化 解 析 手 法 に つ い て 」・ ・ ・ ・ 5
「 オ ー プ ン ソ ー ス C F D E M を 用 い た C F D - D E M 連 成 解 析 に つ い て 」・ ・ ・ ・ ・ ・ ・ ・ 9
「 自 動 メ ッ シ ュ 生 成 機 能 を 有 す る 3 次 元 き 裂 進 展 解 析 シ ス テ ム 」・ ・ ・ ・ ・ ・ ・ ・ 1 2
○技術レポート
「配管網 脈動圧シミュレーション技術の紹介」・・・・・・・・・・・・・・・・・・・・17
「 乱 流 モ デ ル と は な に か 」・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ 1 9
「SPH 法のプログラム開発と解析事例」 ・・・・・・・・・・・・・・・・・・・・・・・・21
「ホワイトボディ振動解析」・・・・・・・・・・・・・・・・・・・・・・・・・・・・・23
「地盤専用ＦＥＭプログラム「PLAXIS」解析事例紹介」・・・・・・・・・・・・・・・・・25
「 E XCE L を 使 用 し た 操 業 シ ミ ュ レ ー タ 」
・・・・・・・・・・・・・・・・・・・・・・・・27
○製品紹介
汎用有限要素解析ソフト Marc2013.1 ・・・・・・・・・・・・・・・・・・・・・・・・・29
統合流体構造連成解析システム「F.S.Comm」・・・・・・・・・・・・・・・・・・・・・・30
FLUENT 用流体構造連成解析モジュール「F.S.Addin」・・・・・・・・・・・・・・・・・・32
光ファイバセンシングの概要・・・・・・・・・・・・・・・・・・・・・・・・・・・・・34
○トピックス
・・・・・・・・・・・・・・・・・・・・・・・・・・36
「FOCUS スパコンの利用について」
「九州デジタルエンジニアリング研究会（KDK）の活動紹介」
・・・・・・・・・・・・・・・37
「日本機械学会
第２６回計算力学講演会 出席報告」 ・・・・・・・・・・・・・・・・38
技報
第 10 号発刊のご挨拶「１号∼9 号を振り返って」
2014 年 1 月
ｼﾐｭﾚｰｼｮﾝｴﾝｼﾞﾆｱﾘﾝｸﾞ･ｿﾘｭｰｼｮﾝ部
日本機械学会
計算力学技術者認定
上級ｱﾅﾘｽﾄ
ﾘｰﾀﾞｰ
博士（工学）上級ｱﾅﾘｽﾄ
湯本 淳
大神
勝城
NS プラント設計株式会社 シミュレーションエンジニアリング・ソリューション部は
2002 年 か ら NPD 数 値 解 析 技 報 を 発 刊 し 、 当 社 ホ ー ム ペ ー ジ （ URL ：
http://www.npd.eng.nssmc.com）にて掲載しております。今回の技報 10 号は、ここ 10 年
間の区切りの号となりましたので、以下に これまで発刊した技術内容を振り返り、当社
におけるシミュレーション技術の推移について述べてみたいと思います。現在 NPD が保有
する高度解析技術の“魁”になったものについてはコメント（＊印）を記入しております。
少々、挨拶が長くなりますが 最後まで通読していただければ幸いです。
１．これまで NPD 数値解析技報の発刊内容
●第 1 号（創刊号）2002 年 7 月
技術論文: 「有限要素法における座標変換則について」
トピックス：Marc_DDM、肉盛溶接解析(by Marc)
事例：加熱炉ガス流動(by Fluent)、スティール缶落下(by LS-dyna)、
誘導加熱(by Fledy)
＊陽解法ソルバー LS-dyna をエンジニアリングへ適用開始。
解析対象を衝撃問題へと広げた。
＊当時、溶接解析に関する標準機能がなかったため、
Marc ユーザサブルーチンを使って構築した。
●第２号 2003 年 4 月
技術論文: 「有限要素法をつくる」
トピックス：流体構造連成問題（Marc⇔Fluent)、
領域結合システム開発（NPD 自作：拘束条件自動設定ツール）、
Marc_DDM による大規模解析（対象：機械継手）、
ドラム缶巻き締め解析(by Marc)
事例：円筒の衝撃解析(by Dytran)、
容器の揺動解析(by LS-dyna：SPH 法)、
液滴解析(by Fluent：混相流、自由表面問題)
＊流体構造連成解析を可能とする自作インターフェース（FSI）
の開発開始。＝＞現在の NPD のツール商品に仕上げた。
＊同様に、領域結合システム開発により、作業効率を向上させるとともにツールとて商品化した。
●第３号 2003 年 4 月
技術論文: 「シェル要素用いたい複合非線形有限要素法について（その１）」
トピックス：拡散燃焼問題(by Fluent：反応)、
耐火材熱衝撃問題(by Marc)、解析システム開発
（NPD 自作：エクセルと Marc の連携）、
オープンソース Adventure＿Pro.による解析
事例：ボルト締結力(by Marc)、土石流問題＆風圧力による板の振舞
(by LS-dyna：SPH 法)、
ブレーキ・ディスク冷却解析(by Fluent：Moving_Mesh の適用）
数値解析技報 第 10 号，2014
1
＊エクセルを GUI とした解析ツールは その後も多数 受託開発を行い、NPD の保有技術の１つとなった。
＊Fluent に MovingMesh（境界を移動回転させる機能）が取り入れされ、軸の回転など機構を有する流動
現象へのシミュレーションが可能となった。
●第４号 2004 年 10 月
技術論文: 「シェル要素用いたい複合非線形有限要素法について（その２）」
トピックス：インターネット解析(NPD と阪大の共同研究)、
流体・構造連成解析(by LS-dyna：ALE 法、オイラー法、SPH 法)
事例：低減積分シェル要素を使用した異方性深絞り解析(by Marc)、
球形永久磁石の静磁場解析(by Marc)、3 次元音響解析(by Nastran)、
造波解析(by Fluent：２相流、自由表面）、
コンクリート・クラック問題(by LS-dyna)
＊LS-dyna の流体・構造連成解析機能に注目し、
特に粒子法の取り扱いを習得し、その後の具体的テーマに適用した。
＊音響（騒音）解析の需要が住宅関連、環境騒音で予想されたため、
まずは解析技術の習得を行った。
●第５号 2005 年 9 月
技術論文: 「3 次元静磁場辺要素について」
技術レポート：流体−構造連成解析システムの開発(NPD 自作の FSI)、
ダイナミックダンパーによる防振効果の確認(by Marc)、ケーブル要素の利用に関
する検討(by Marc)、加振テストに関する 3 次元構造-音場連成解析(by
Nastran)
事例： 溶接シミュレーション機能(by Marc)、動的応答解析手法を用いた座屈
解析(by Marc)、薄板のしわシミュレーション(by Marc)、εーN 寿命予測解析
(MSC.Fatigue)、地震発生時のタンク内液の挙動解析(by Fluent：２相流、自由
表面 VOF）
＊自作 FSI の機能向上、音場解析レベルの向上
＊Fluent による混相流の解析へ注力しだした＝＞その後、気―液、気―液―固へ課題へ発展させた
●第６号 2006 年 9 月
技術論文: 3 次元静磁場辺要素について（続編）
技術レポート：バウシンガー効果を考慮した交番塑性構成則について(NPD と阪
大との共同研究： Marc ユーザーサブルーチンの作成)、水中における板の振動解
析 （ 連 成 問 題） （ by FSI:Marc ⇔ Fluent) 、 砂 防 ダム 緩 衝 材 の 衝 撃 解 析 (by
LS-dyna)、Marc をステーション化した 4 層骨組の地震応答統合解析(NPD と
阪大との共同研究)、群杭を対象とした地盤解析事例(by Marc)
製品紹介：FSI、光ファイバーセンシング、座屈問題における人工減衰他の新機
能紹介 (by Marc)、並列解析機能(by Fluent）
＊自作 FSI の機能向上（水中の板振動）
＊阪大との共同研究による統合化解析手法の構築
＊光ファイバーセンシング（商品化に成功）＝＞その後、環境プラント設備の温度計測などに適用
2
数値解析技報 第 10 号，2014
●第７号 2008 年 1 月 特集：「ものづくりを支えるシミュレーション技術」
（当時の新日鉄エンジニアリング㈱ 技術部長の意見を掲載）
技術論文: 鉄鋼プロセスにおける電磁誘導加熱の高品位化対策に関する研究
技術レポート：FLUENT 混相流解析モデルについて、光ファイバセンシングの研究
動向、音響解析の基礎的検討(by Nastran)、Marc/HYPELA2 による材料構成
則の組み込み、落石防止柵の間伐材に対する LS-DYNA の適用事例、水中に
おける円柱の自由振動シミュレーションの検証(NPD 自作の FSI 適用)
機能紹介：亀裂進展（VCCT 法：Marc2007r1 機能）、ウォーターハンマー現象の
表現(by Fluent）
＊LS-dyna により、衝撃によるエネルギー吸収を検討する事例が増加。
＊Marc で亀裂伝播解析機能が使用できるようになった。（ただし、２次元ベース） MSC にて開発継続中。
●第８号 2010 年 8 月 技術論文 1: 合成床版の輪荷重走行試験に対するシ
ミュレーション手法の一提案(by Marc：コンクリートのひび割れ問題)、
技術論文 2: プロセス配管系流体騒音問題に対するシミュレーション手法の
一提案(by Marc＆ Fluent）
技術レポート：マイクロ AVS によるカスタマイズ例（可視化事例）、オープンソース
OpenFOAM の紹介、統合型ソルバー非線形機能調査(by SolidWorks)、燃焼
モデルの説明(by Fluent）、電磁界における電磁力表現(by Marc)、輻射形態係
数の計算精度(by Marc)、延性破壊解析(by Marc)、３次元疲労亀裂進展解
析(by Marc：NPD と東京理科大の共同研究)、
製品紹介：Marc2010 新機能紹介
＊オープンソース OpenFOAM の利用開始。
＊応力解析を Marc で行う亀裂伝播解析システムを大学との共同研究により開発開始。
（技報 10 号の技術論文 参照）
●第９号 2011 年 12 月
特集 ：「大規模解析の実現に向けて」
技術論文: OpenFOAMを用いたアプリケーション・ソフトウェアの開発
技術レポート：LS-DYNA による粉体圧縮挙動シミュレーション、実験計画法と
FEM による最適構造設計(by Marc)、非線形座屈解析における問題と対策
(by Marc)
製品紹介：Marc 2011 新機能紹介
＊オープンソース OpenFOAM の利用技術を向上させた
＝＞ 現在、さらに高度化を狙い継続中
＊LS-dyna による粒子法（SPH 法）難度の高い粉体移送のテーマを
対象にして検討（構成則の精度を向上させることが課題）
数値解析技報 第 10 号，2014
3
２．10 年間を振り返って
以上を俯瞰すると次のように整理することができます。
（１）ここ 10 年間で、当社の解析技術分野が広がり、かつ新規技術の獲得と既存技術の
高度化を併せて達成してきた。
（１．参照）
（２）具体的には、構造解析分野では、静的∼振動∼衝撃にいたる時間軸の解析、コン
クリートひび割れ問題、地盤解析、粒子法（SPH)の適用、音響問題、さらにエクセル
を用いた解析システム開発まで技術の幅が広がった。一方、非線形弾塑性構成則の構
築や、統合化解析手法、および亀裂伝播問題に関しては、大学との共同研究による取
り組みにより新規技術の獲得と高度化を達成した。
（３）流体解析分野では、早い時期から混相流の分野を開拓するとともに、燃焼問題（反
応）
、流体構造連成問題に取り組んでおり順調に解析技術の幅が広がってきた。後半
期では、オープンソースによる解析技術の獲得を行い、流動解析、DEM（粒子法）開
発、さらに DEM と CFD の連成問題まで新規技術の獲得と高度化を行ってきた。
（その成果の一部を 10 号の技術論文に掲載）
（４）電磁場解析分野では、プラント設備の誘導加熱装置の性能検証のため、継続して
解析を実施しており、現在 数値解析分野の１つの柱に成長してきた。
（５）また、MSC ソフトウエアーの販売代理店および Marc ユーザーの立場から、Marc 新
機能と Marc の解析事例を一般ユーザーに対して継続して発信した。また、別途 Marc
ユーザーに対しては、１年に一度のペースで北九州において定期的に CAE 活用事例セ
ミナーを開催し、今年で 16 回目となった。
（６）後半期では、可視化技術に関するシステム開発の例を技報に掲載した。特に光フ
ァイバーセンシング技術は、プラント設備の操業監視（温度、ひずみ計測）に適用し
ており、シミュレーションから現場での計測まで含めた総合シミュレーションエンジ
ニアリングの一端を担っている。
（７）１．では特に記載しなかったが、対外活動として前述の CAE 活用事例セミナーの
開催、および KDK（九州デジタルエンジニアリング研究会）の副会長と事務局を歴任
してきた。また、機械学会 計算力学技術者認定 付帯講習会 講師を担う等 九州
地区の地場関係者へ CAE 技術の浸透に貢献したとのことで、2012.3 月 日本機械学会
から九州支部賞を授与された。
３．今後の 10 年間に向けて
今後とも、
「シミュレーション技術による社会へ貢献」を組織の使命として努力してまい
ります。その一つの取り組みとして、NPD 数値解析技報による一般ユーザーの皆様への情報
発信を継続することが、その一端を担うと考えております。
先の 2011 年 3 月に発生した東北地方太平洋沖地震による被害状況を考えると、被害の予
測と事前の対策がいかに重要であるかを考えさせられました。今後、津波荷重の予測や防
潮堤や防波堤の基礎地盤の安全性を粒子法（DEM)を用いて検討することが国レベルで進ん
でいくと思われます。私共もエンジニアリング会社の組織の一員として、技術的に貢献で
きるように切磋琢磨していこうと考えています。
その他、国の懸案となっている環境問題、エネルギー資源開発、製造プロセスの効率化
など枚挙にいとまがないのですが、技術課題の解析難度が高くなることが予想されます。
それらの分野では連成場（流体∼構造∼化学反応∼電磁∼伝熱）の検討が必要になってき
ます。私共は これらの分野で一層力を発揮できるように努力してまいりますので、今後
とも よろしく お願いいたします。
敬具
4
数値解析技報 第 10 号，2014
技術論文１．
「建築鋼構造分野における異種プログラムを連携した統合化解析手法について」
統合化解析手法とは、大阪大学大学院 多田元英教授が提唱した鉄骨構造の非線形挙動を並列処理で解
析する手法の呼称である。当社は、多田教授が 2002 年に立ち上げたインターネット解析研究会に参画し現
在に至っている。統合化解析手法の目的は、研究者が開発した単機能プログラム群を、相互に統合するこ
とで複合機能を有する数値解析を達成することにある。これを多田教授は、
「異種プログラム間のコラボレ
ーション」と称している。
大学における研究対象の１つに、建物全体の倒壊挙動解析の実現があり、2009 年度には、E－DEFENCE
で実施された実大 4 層鉄骨造建物の神戸地震動を想定した立体倒壊実験の実験シミュレーションを行い、
倒壊挙動を精度よく再現しうることを示した。
なお、以降に記載する内容は、機械学会第２６回計算力学講演会にて主著者 多田教授、 共著者
大神にて論文発表したものである。
博士（工学） 大神 勝城
1.緒言
ここでは，民間で使用されている汎用非線形構造解析ソルバーMarc と大学で開発された骨組構造解析ソル
バー等を連携させ，より多機能化した統合化解析手法について説明する．さらに，E-ディフェンスにて震
動台実験が実施された実大 4 層鋼構造建物の完全倒壊実験に対し，この解析手法を適用した事例を述べる．
それらのソルバーは，ネットワーク回線を介して分散処理を行う．勿論，解析対象によっては，プログラ
ム構成を変更することは可能である．今回は，骨組構造解析ソルバーNETLYS1）をホストとし，部材に発生
する局部座屈を考慮できるようにシェル要素を用いた解析ソルバーをステーションとして用いた．
2. 統合化構造解析システム
2-1 解析法の目的と概要
骨組の倒壊挙動は局部座屈をはじめとした複数の不
安定現象が互いに連成しており，さまざまな研究者が
開発した詳細解析ツールをネットワーク回線により
統合すれば，従来になく高い解析精度で解析すること
が可能である．文献 1 では複数の研究者が開発した詳
細解析ツールを統合して建物全体を高精度に解析し
ようという統合化構造解析システムを発表した．複数
の既存プログラムを図 1 に示すように統合することで，
各構造部分の挙動を詳細に解析しながら互いの連成
効果を考慮しながら建物全体の挙動を解明するのが
本システムである．
図1．統合化解析システムの構成イメージ
2-2 局部座屈解析ソルバーのステーション化方法
統合化解析のデータ通信の流れを図 2 に示すが，ホストからステーションへ入力条件として通信するのは，
ステーションが担当する部分構造の境界自由度に対応した変位増分②であり，一方 ステーションからホ
ストへ通信する解析結果はそれらの復元力成分④と境界自由度に縮約した剛性行列（①、③）である．こ
のとき，部分構造の剛性については、弾塑性や局部座屈の影響により刻々変化するため，ステップ毎に算
出する．具体的には，図 3 に示すように，部分構造全体の接線剛性行列を境界自由度に対応する成分[K*11]
に縮約し，それをホストに送信する.
②境界変位増分3成分：Δu,Δθi,Δθj
骨組解析
（ホスト）
i
j
u
・・・境界自由度成分
局部座屈解析
（ステーション）
・・・内部の自由度成分
局部座屈解析の境界変位
④境界力3成分(N,Mi,Mj )
①、③：３ｘ３剛性行列（Kij)
図２． 統合解析の通信の模式図
図３．ステーションにおける剛性行列の縮約について
数値解析技報 第 10 号，2014
5
これらの手順については，汎用ソルバーMarc の標準
機能では対応できないため，ユーザーサブルーチンに
必要な機能を組み込んだ．これら 4 種類のユーザーサ
ブルーチンについて表 1 に整理している．なお，ステ
ップ毎の接線剛性行列を抽出してホストに送信する
ことは汎用ソルバーでは難しい．そこで，運動方程式
を数値積分する際に，接線剛性行列の情報が不要で部
分構造の復元力だけで求解可能な OS 法 2)を第 3 節の
統合解析では用いている．この OS 法は，数値解析と
構造実験を組み合せた「オンライン実験」または「サ
ブストラクチャ仮動的実験」3)などで従来利用されて
いる方法である．
表１． ステーションで使用したユーザーサブルーチンの説明
NO
ユーザーサブルーチン名
1
SUBROUTINE
UBGINC(INC,INCSUB)
2
SUBROUTINE
FORCDT(U,V,A,・・・）
3
SUBROUTINE
IMPD (LNODE,DD,・・・)
4
SUBROUTINE UEDINC
(INC,INCSUB)
用途
入力
ホストからの
通信情報
変数の初期化機能
荷重設定機能。ホストから
の通信情報（変位増分）を
境界節点の強制変位成分と
して入力。
復元力抽出機能。インクリメ
ント毎に境界節点の強制変
位成分に対応する復元力
（反力）を出力。
最終INCで計算をSTOPする
機能。およびインクリメント
毎に解が収束すれば、復元
力その他をホストへ通信。
出力
Marc側INC番号
通信情報、
節点番号
境界節点番号
自由度、
強制変位増分
通信情報、
状態変数情報
復元力
通信情報、
状態変数情報
節点番号
復元力、
通信情報
2-3 局部座屈の影響を考慮した門型フレームの解析事例
ここでは，局部座屈解析ソルバーに NASP4）を使用した例を示す．なお，Marc を適用した事例については，
第 3 節で説明する．
(1) 解析対象の説明
H-400×200×8×13
図 4 に対象とする門型フレームを示す．局部
① N= 0.3×Nｙ
座屈解析の対象は，左右の柱基部に材軸方向
長さが 700ｍｍの範囲，および梁の左右端に
②
②δ
′-350×12
800ｍｍの範囲の計 4 箇所とした．なお，前述
δ=200 mm
N=1581.9ｋN
4000
の材軸方向長さは部材せいの 2 倍である．材
料特性は，ヤング率 20.58GPa，降伏点 325 MPa,
：局部座屈検討箇所
および歪硬化係数を極小とした．
単位：ｍｍ 8000
(2) 荷重条件および拘束条件
図５． 荷重条件
図４．モデル寸法
図 5 に示すように，骨組の境界条件は柱端部
を固定し，柱に降伏軸力の 30％の軸力を掛けた後に水平方向に回転角 1/20 まで強制変位を 200 回の増分
で載荷した．
せん断力（ｋN)
600
δ
500
θ
400
P
θ
175ｍｍ
座屈振幅 約30ｍｍ
H
300
統合解析
骨組解析
200
100
部材角：θ=δ/H
せん断力 ：P
部材角（rad)
回転角
変形倍率：２倍
0
0
0.01
0.02
0.03
0.04
0.05
図６．門型フレームの荷重-変位関係
下端側
断面の座屈形状
図７．右柱基部の局部座屈（回転角θ＝0.05のとき）
(3) 骨組解析結果と統合化解析結果の比較と考察
NETLYS の骨組解析では軸力の影響を考慮した一般化硬化塑性ヒンジ法を適用した．NASP はシェルの高次要
素を用いており，塑性構成則にはミーゼスの降伏条件を適用した．図 6 に横軸を回転角，縦軸をせん断力
の解析結果，および図 7 に局部座屈モードを示す．骨組解析の結果と今回の統合解析結果を比較すると，
後者は 1)初期剛性がほとんど一致，2）最大耐力が 10％ほど大きい，3）左右の柱基部に発生した局部座
屈により劣化勾配が鋭い，などの特徴が出ている．
なお，2）は統合解析において局部座屈解析の
多軸方向応力場の塑性構成則の影響が大きく出たと考えられる．
3. 実大 4 層鉄骨造建物の立体倒壊解析 5)
3-1 解析対象の試験体
解析対象の試験体の模式図と立体解析法の説明を図 8 に示す．各平面架構の柱梁節点の水平変位を，各層
床の平面重心位置における X または Y 方向水平変位と床のねじれ回転に伴う水平変位との和に等置するこ
とで立体挙動を表現する．
6
数値解析技報 第 10 号，2014
表 2 に素材試験結果を示す．柱は冷間成形角形鋼管で材質が BCR295 であり，梁は H 形鋼で材質が SN490B
である．柱梁接合部は通しダイヤフラム形式で接合している．屋根は厚さが 150 mm の鉄筋コンクリート
スラブであり，2∼ 4 階床は高さが 75 mm のデッキプレートと厚さが 100 mm の等厚部からなる合成床版
である．柱脚
表2．素材試験結果
は 50 mm のベ
ースプレート
と M36 の ア
Z
ンカーボルト
8 本からなる
Y
露出形式であ
り，震動台上
X
の鋼製基礎に
図８．立体解析法
アンカーボル
トで緊結され
ている．
3-2 統合解析のプログラム編成
柱は，2 軸曲げ，軸力，せん断力を考慮し
た立体解析を行い，梁と柱梁接合部パネル
と柱脚は平面解析を行う．本統合化構造解
析では，全体制御を司るホストプログラム
に NETLYS を用いる．図 9 に本統合化構造
解析で使用するプログラムの編成を示す．
ステーションプログラムに合成梁解析プ
ログラム COMPO と，汎用ソルバーMARC を
用いる．後者は，1 階の層崩壊の起点とな
る角形鋼管柱の端部の局部座屈を表現す
る．図 10 に柱の解析モデルを示す．
図９．統合化解析のプログラム編成
3-3 地震応答解析結果
(1) 入力地震動と解析条件
入力地震動に震動実験時に入力した JR Takatori 波を
用いる．X 方向に EW 成分，Y 方向に NS 成分，Z 方向に
UD 成分を与え， 100% 波について地震応答解析をする．
60% 波と 100% 波の震動台上で観測した最大速度を表 3
に示す． 100% 波では 2 次設計レベルの約 2.5 倍の最
大速度を記録し，設計での想定を大きく上回った地震動
が入力されたことがわかる．試験体の固有周期の解析結
果を表 4 に示す．減衰特性は，X 方向の 1 次(0.81s) と
2 次(0.26s) の固有周期に対する減衰定数が 0.02 の
レーリー型とする．
(2) 統合解析結果と実験結果の比較
以降に解析結果と実験結果の比較を示す．図 11，12 の
1 階層間変形角の時刻歴に示すように，1 階に
おいて X 方向の負側に，Y 方向の正側に層間
変形角が大きくなり，倒壊に至ったことが，解
析で追跡できていることがわかる．
図 13 に示す 1 階の層せん断力−層間変形角関
係において，局部座屈による耐力劣化が解析結
果と実験結果で概ね一致していることがわか
る．また、この 1 階の柱の局部座屈による層
崩壊の様子が，図 14 に示す解析終了時の塑性
化部位の分布の図からもわかる．
図10．１階柱の構造モデル
図１１．１階の層間変形角の時刻歴（0.6×Takatori)
数値解析技報 第 10 号，2014
7
さらに、図 15 に示すように、局部座屈モードについても実験結果とほぼ一致していることがわかる。
図１2．１階の層間変形角の時刻歴（1.0×Takatori)
図１3．１階の層せん断力-層間変形角関係（1.0×Takatori)
震動台実験
図１4．解析終了時の塑性化部位の分布（1.0×Takatori)
倒壊解析結果（by Marc)
図１5．柱端部の局部座屈モード
4. 結言
本報では，異種プログラムを統合した詳細な解析システムについて説明し，その中で局部座屈解析ソル
バーを統合する手法について述べた．実用例として E- ディフェンスにて震動台実験が実施された実大 4
層鉄骨造建物に対して，統合化構造解析システムを用いた立体解析について述べた．
これらの成果を整理すると以下のとおりである．
(1) 異種プログラムを統合することで，局部座屈の表現など多機能かつ詳細な解析システムを構築
できることを示した．
(2) 局部座屈解析ソルバーの統合にあたっては，担当する部分構造の境界変位成分に対応した復元力と
剛性行列を計算し，ホスト側へ通信する仕組みを構築することで実現することを示した．
(3) また，実大 4 層鉄骨造建物の倒壊実験に関して，JR Takatori 100% 波に対する応答にて，
損傷が 1 階に集中して柱の両端に局部座屈が起こり，1 階の層崩壊によって建物が倒壊する
現象を精度よく表現できた．
以上
＜参考文献＞
1) 多田元英，桑原 進：インターネットで異種プログラムを統合した構造解析システムの基本考察，日本建築学会構造
系論文集，No.580，pp.113-120，2004.6
2) 多田元英，潘 鵬：統合化構造解析へのオペレータ・スプリッティング法の適用，日本建築学会構造系論文集，第 604
号，pp.167-174，2006.6
3) 中島正愛，石田雅利，安藤和博：サブストラクチャ仮動的実験のための数値積分法(サブストラクチャ法を用いた仮
動的実験の開発)，日本建築学会構造系論文集，第 417 号，pp.107-117，1990.11
4) 大神勝城，多田元英：局部座屈解析と骨組解析を統合した門型フレームの複合非線形解析，日本建築学会大会学術講
演梗慨集，C-1，pp.991-992，2004 年 8 月
5) 堀本明伸，多田元英，玉井宏章，大神勝城，桑原 進，三谷 淳：実大４層鉄骨造建物の統合化構造解析システムに
よる立体倒壊解析，構造工学論文集，Vol. 55B，pp.277-283，2009.3.
8
数値解析技報 第 10 号，2014
技術論文２．
「オープンソース
オープンソース CFDEM を用いた CFD-DEM 連成解析について」
連成解析について
コンピュータ性能の格段の向上により
格段の向上により、プラント設備、海洋・土木分野などにおいて粒子法を用い
、プラント設備、海洋・土木分野などにおいて粒子法を用い
たシミュレーションが活発になってきている。また、粒子的挙動と連続体力学を前提にした流体解析
（CFD）
、構造解析（FEM）などを
を組み合わせることで製造プロセスから、防災に関わる海底地盤の挙
製造プロセスから、防災に関わる海底地盤の挙
動検討にいたる幅広い分野で適用が可能である
動検討にいたる幅広い分野で適用が可能であると考えられる。
。ここでは、当社が開発中の粒子法
当社が開発中の
（DEM）
について説明する。なお、以降に記載する内容は、機械学会第２６回計算力学講演会にて
以降に記載する内容は、機械学会第２６回計算力学講演会にて
著者 首藤が論文発表したものである。
論文発表したものである。
水産学博士
1.はじめに
当社ではクラウドコンピューティングサービスや自社の計
算機環境を有効活用するために，解析業務へのオープンソー
スの利用を推進している．その取り組みの中から離散要素法
を用いた流体−粒子連成解析例を紹介する．以下では，まず
流体−粒子連成解析に用いたオープンソース CFDEM1）の概要
を説明し，次に解析例として二成分粒子の気系流動層におけ
る粒子混合と偏析の解析，および伝熱や反応を伴う流動層の
解析を紹介する．
首藤 史
CFDEMソルバ
LIGGGHTS
スタート
DEM粒子データ受信
CFDにデータ転送
セルに粒子を配置
粒子運動の計算
体積率と運動量交換量の計算
CFDデータの受信
2. オープンソース CFDEM の概要
DEMにデータを転送(MPI)
CFDEM1）は CFD−DEM 連成解析を可能としたオープンソース
2）
であり，OpenFOAM をベースとした CFD ソルバーと DEM ソフ
流れの計算
トウェア LIGGGHTS3）で構成されている 4）5）．
CFDEM では CFD-DEM
連成解析手法として Euler-Lagrange
Lagrange 法と埋め込み境界法の
２種類のソルバーが提供されているが，ここでは本報の解析
図1．CFDEM
CFDEMの計算手順
に使用した Euler-Lagrange 法のソルバーについてのみ説明
する．
図 1 に Euler-Lagrange 法の CFDEM の計算手順を示す．まず CFD 側では DEM 粒子データ (粒子の中
心座標値と速度) の受信後，流体計算セルに粒子を配置する．次に各流体計算セルにおける流体体積
率αf (= 1 - 粒子体積率)と粒子−流体間の運動量交換量
と粒子−流体間の運動量交換量 (流体力) が計算され，流体力は DEM 側に
転送され，流体体積率と流体力は流体解析に用いられる．
転送され，流体体積率と流体力は流体解析に用いられる．DEM 側では，流体力の受信後，粒子運動を
計算し，粒子の中心座標値と速度を CFD 側へ転送する．ここで，CFD と DEM 間のデータ送受信には MPI
が用いられている．また，流体力 (抗力，浮力，仮想質量力，揚力)のモデルは既に組み込まれてい
のモデルは既に組み込まれてい
るモデルを参考にすれば，新たなモデルを組み込むことも容易である．
ば，新たなモデルを組み込むことも容易である．
3. オープンソース CFDEM を用いた CFD−DEM 連成解析事例
3-1 二成分粒子の流動層における粒子混合と偏析の解析
粒子の流動層における粒子混合と偏析の解析
ここでは文献 6 を参考にして実施した，矩形流動層 (幅 ＝ 0.065 m, 高さ ＝ 0.26 m, 奥行 ＝
0.0081 m) における二種混合粒子（小粒子径 ＝ 0.001 m および 大粒子径 ＝ 0.002 m）の流動化シ
ミュレーションについて紹介する．総粒子数は 25000 で，各粒径の割合は重量で 50%ずつである．ガ
スは底面から一様に入れている．
一様に入れている．シミュレーションに用いたパラメータおよび抗力係数モデルは文献
6 と同じである．このシミュレーションは CFDEM が提供しているソルバーをそのまま使用すれば実施
できる．
図 2 に流入ガス流速 1.4 および 1.6 m/s における粒子流動状態の時間変化を示す．まず，ガス流速
1.4 m/s では小粒子が層上部に，大粒子が層下部に偏在していることが確認できる．一方，ガス流速
1.6 m/s になると層全体が流動化し，層底部を除くと大小粒子が良く混合していることがわかる．
図 3 に実験 6)と解析における粒子流動状態の比較を示す．実験の方は暗い色が小粒子，明るい色が
大粒子である．この図から流入ガス流速による大粒子の層高さの違いが実験と解析で良く一致してい
ることがわかる．
数値解析技
数値解析技報
第 10 号，2014
9
ガス流速＝1.4 m/s
0s
2s
5s
10s
15s
実験
20s
解析
実験
ガス流速＝1.4 m/s
ガス流速＝1.6 m/s
解析
ガス流速＝1.6 m/s
図3．実験6) と解析における粒子流動状態の比較
[K]
0s
2s
5s
10s
15s
20s
図2．流入ガス流速1.4および1.6 m/sにおける
粒子流動状態の時間変化 (青：小粒子，赤：大粒子)
3-2 伝熱を伴う流動層の解析
伝熱を伴う流動層を扱う場合，粒子間およ
び粒子−壁間の接触および輻射伝熱や粒子−
流体間の対流伝熱の考慮が必要となる．この
内，輻射伝熱以外は CFDEM に実装されている
ため取り扱うことができる．また，CFDEM は
非圧縮性流体ソルバーのみ提供しているため，
圧縮性流体およびエネルギー方程式のソルバ
ーは新たに作成する必要があるが，これらは
OpenFOAM を用いて比較的容易に作成できる．
図 4 に伝熱を伴う流動層の解析例を示す．
初期状態は粒子温度＝400 K,ガス温度＝292 K
であり，底部からの流入ガス温度は 292 K で
ある．この図から時間と共にガスに熱を奪わ
れ粒子温度が低下して行く様子がわかる．
10
数値解析技報 第 10 号，2014
0.2s
0.5s
1s
2s
3s
2s
3s
ガス温度
0.2s
0.5s
1s
粒子温度
図4．伝熱を伴う流動層におけるガス温度と粒子温度の時間変化
3-3 反応による粒径変化を伴う流動層の解析
反応を伴う流動層を扱う場合，反応による
粒子−流体間の物質移動やそれに伴う粒子密
度や粒径の変化を考慮する必要がある．これ
らを考慮するには LIGGGHTS のカスタマイズ
が必要になる．また，化学種輸送方程式のソ
ルバーも新たに作成する必要がある．
図 5 に反応による粒径変化を伴う流動層の
解析例を示す．これは粒子の揮発分がガスに
放出され，時間と共に粒径が縮小して行く様
子を解析したものである．
粒径
[m]
0s
0.5s
1s
1.5s
2s
図5．反応による粒径変化を伴う粒子の流動状態の時間変化
4. まとめ
以上では，オープンソース CFDEM の概要を説明し，その適用例として二成分粒子の流動層における
粒子混合と偏析の解析，および伝熱や反応を伴う流動層の解析を紹介した．二成分粒子の流動層の解
析については実験と解析の粒子流動状態を比較した結果，CFDEM により粒子の混合や偏析を精度良く
予測できることを確認した．したがって，CFDEM は流体−粒子連成解析の業務に利用できるレベルに
あると考えている．また，伝熱や反応を伴う流動層の解析例では，より複雑な現象を考慮した CFD-DEM
連成解析が CFDEM のカスタマイズにより実現できることを示した．
以上
＜参考文献＞
1) CFDEM, Open Source CFD, DEM and CFD-DEM, URL: www.cfdem.com
2) OpenFOAM, The Open Source CFD-toolbox, URL: www.openfoam.com
3) LIGGGHTS, LAMMPS improved for general granular and granular heat transfer simulations, URL:
www.liggghts.com
4) Kloss C., Goniva C, Aichinger and G., Pirker S., Comprehensive DEM-DPM-CFD simulations - model synthesis,
experimental validation and scalability, Seventh International Conference on CFD in the Minerals and
Process Industries, 2009.
5) Goniva C., Kloss C., Hager A. and Pirker S., An Open Source CFD-DEM Perspective, 5th OpenFOAM Workshop,
2010.
6) Feng, Y. Q., Yu, A. B., Assessment of model formulations in the discrete particle Simulation of gas-solid
flow. Industrial and Engineering Chemistry Research 43, pp. 8378-8390, 2004.
数値解析技報 第 10 号，2014
11
技術論文３．
「自動メッシュ生成機能を有する
自動メッシュ生成機能を有する 3 次元き裂進展解析システム
解析システム」
鋼構造物の疲労問題は、プラント設備の保全
は、プラント設備の保全において重要な課題の１つである。特に、
重要な課題の１つである。特に、溶接継手部
の疲労き裂の発生と進展について精度良く予測できれば、
の発生と進展について精度良く予測できれば、保全対策に有効であることは言うまでもな
い。き裂進展解析では、欠陥が内在すること
欠陥が内在することを前提として、荷重履歴によるき裂
き裂の広がりを数値解析
で自動的に追跡することが必要であるが、種々の解析上の課題があり、まだ
することが必要であるが、種々の解析上の課題があり、まだ研究途上
研究途上にある。当社で
は、東京理科大 岡田教授のご指導をうけ、
岡田教授のご指導をうけ、Marc をソルバーとしたき裂進展解析
解析システムを構築し
た。さらに、基本的な課題による
基本的な課題による解析精度の検証を行い、実用的な課題にも取り組んだので本論文に
実用的な課題にも取り組んだので本論文に
て説明する。
博士(工学)
大神 勝城
計算（固体）力学技術者 １級認定 須藤 浩章
１．はじめに
米国の戦時標準船（リバティー船）の溶接部に脆性破壊
脆性破壊が生じ、船体
第 2 次世界大戦の戦中戦後、米国の戦時標準船
が破断した事故が多発した。それを
が多発した。それを切掛けにして、20 世紀半ばから破壊力学に関する
に関する研究が始まっ
た。それらの技術成果が線形破壊力学として確立され、き裂先端の応力状態を応力拡大係数という指
先端の応力状態を応力拡大係数という指
標で表現し、き裂進展則に取り入れることで実験結果を
進展則に取り入れることで実験結果を定性的に説明できるようになった。
説明できるようになった。
ただし、実構造物のき裂進展量を精度良く予測することは一般的に難しく、
進展量を精度良く予測することは一般的に難しく、FEM
FEM と破壊力学を統合
する手法が現在も研究されている。
研究されている。その重要な技術課題として、①き裂先端の応力拡大係数を精度よ
先端の応力拡大係数を精度よ
く算出するメッシュ作成法、および
く算出するメッシュ作成法、および②き裂進展方向の自動化が挙げられる。ここでは、
ここでは、当社で開発し
た解析システムについて概要を説明し、基本的な問題、実用問題の２例について
実用問題の２例についてき裂進展解析結果を
説明する。
２．き裂進展解析手法について
２−１．本手法の特徴
面体
一般的に、任意の三次元形状で応力拡大係数を精度よく計算するためには 6 面体要素のメッシュ作
成をする必要があるが、き裂前縁部のような複雑な形状に対しては非常に難しく
前縁部のような複雑な形状に対しては非常に難しく 多大なコストがか
前縁部のような複雑な形状に対しては非常に難しく、
かると考えられる。一方、岡田 等の研究 1）、2）、3）により、4 面体要素であってもメッシュ分割を適
あってもメッシュ分割を適
正に行うことで、応力拡大係数を
応力拡大係数を精度よく計算する手法が構築された。あえて、
あえて、4 面体要素を使うこ
とによって、6 面体要素に比べて
要素に比べて より簡潔にき裂前縁部周辺をメッシュ分割することを実現した。
メッシュ分割することを実現した。
この手法は、VCCM 法（Virtual
Virtual Crack Closure
Closure-Integral Method）を用いた三次元き裂進展解析手
を用いた三次元き裂進展解析手
法であり、仮想き裂閉口積分法を用いて
を用いて混合モード環境下においても応力拡大係数の各モード成分を
混合モード環境下においても応力拡大係数の各モード成分を
個別に算出することが容易である
である 2）。以下は、粗密を有する空間の節点群を 4 面体メッシュに変換す
る様子をイメージ化したものである。
tetrahedralize
図 1．メッシュイメージ
12
数値解析技報 第 10 号, 2014
２−２．手順説明
き裂進展解析のフローを以下に示す。き裂の自動進展解析は、岡田等の研究 4),5)により実現された。
以下に手順説明を簡単に行う。
（１）通常のき裂なしの構造解析モデルを作成する。別途、初期クラックの位置とき裂形状を設
定して、図 2 の要領に従って初期き ①初期健全解析モデル作成
②初期き裂位置、き裂大きさ
（き裂なしmarcモデル）
データ情報入力
裂形状を考慮した該当領域のみを自
動的にメッシュ分割する。（フェー
③初期き裂有りデータの作成（き裂周辺のみ）
ズ： ①→②→③）
（２）初期き裂部分のメッシュと そ
④き裂部分とその他の部分を合体し解析データ作成
の他の領域を、Marc の機能であるグ
ル―オプション（接着機能）を用い
⑤Marcにて応力解析
て合体する。
（フェーズ：③→④）
（３）通常の構造解析を実行して、変
⑥応力拡大係数,き裂進展方向、き裂進展量の計算
位・応力・歪を算出し、岡田等が開
発した手法に則り応力拡大係数、き
⑦き裂進展後のメッシュデータ作成
裂進展方向を計算し、予め解析者が
繰り返し計算：yes
繰り返し計算：No
設定した「パリスのき裂進展則」に
⑧解析評価
従って、き裂進展量の増分値を算出
する。
（フェーズ：④→⑤→⑥）
図 2．解析手法のフロー
（４）き裂進展後のメッシュ分割を図 1 のフローに従い作成する。
（フェーズ：⑥→⑦）
（５）以降、荷重条件の最後まで増分解析を行う毎に④→⑦の手続きを行い、き裂進展量を予測
する。
２−２．応力拡大係数とき裂進展方向の算出法について
図 3 に示すように、応力拡大係数の各モード成分を個別に算出する。式中の E'はヤング率 E と ポ
アソン比νを用いて,E'=E(平面応力状態), E'=E/(1-ν2)(平面ひずみ状態)である。また、G は各モー
ド毎のエネルギー解放率である。エネルギー解放率は、き裂前縁の節点力や節点開口変位による計算
されることが文献 1）、2）に示されている。また、図 4 に示すように、進展方向については各モードの
応力拡大係数を用いて算出される 5）。
き裂進展方向φ
応力拡大係数

 2 KⅡ
  2 tan 1 
 KⅠeq 
KⅠeq 2


2 
 8( KⅡ ) 

⊿K Ieq  ⊿K I  B ⊿K III
図 4．進展方向の計算式
K I  E ' GIVCCM
GⅡVCCM
KⅡ  VCCM
GⅡ
E ' | GⅡVCCM
GⅢVCCM
KⅢ  VCCM
GⅡ
2μ| GⅢVCCM
図 3．各モードの応力拡大係数
K I
2
1
2
2

 K 2I  4 1  K II   4 2  K III 
2
 1  1 .155 ,  2  1
 K eq 
・・・（1）

da
 C ⊿K eq
dn
m
・・・（2）
数値解析技報 第 10 号, 2014
13
２−３．パリスのき裂進展側の導入
ここでは、き裂前縁にある節点群の混合モード（3 成分）の応力拡大係数を算出して(1)式により
相当応力拡大係数 Keq を算出した。さらに、(2)式のように一般的なパリスのき裂進展則を構成則と
して Marc のユーザーサブルーチンに組み込んだ。ただし、式(2)のｎは繰り返し数、a はき裂長さで
あり、左辺はき裂進展速度を示す。一方、右辺に含まれる C、m は実験で得られた材料定数である。
３．解析手法の検証 ： ﾌﾞﾛｯｸ表面の半円き裂問題
本問題については、解析モデルが簡潔であり、荷重の作用もモード１主体であるため、これまでも
対象になり計算されてきた。ここでは、Raju と J.C.Newman が解析した結果(たとえば、NASA レポー
ト参照）と比較する。
計算モデルを図 5 に示す。前述した手法により、き裂近傍を抜き出してメッシュ分割を行った。ブ
ロックのメッシュの粗密の境界は、Marc のグル―機能により両者を接着している。関係する材料特
性は以下のとおりである。
・ヤング率 E=2.06E+5[MPa],ν=0.3
・パリス則 C 値：7.11×10-13 , m 値：3.65（一般鋼）
荷重条件は 0∼200MPa 面圧を繰り返し端部に与えた。境界条件は、下端の３軸方向の並進変位成分
を拘束した。ここでは、計算時間の短縮を図るため、１回の計算で荷重の繰り返し数を数万程度にま
とめてき裂進展量を計算したところ図 6 に示すようにスムーズに進展した。図 7 に示すように、参照
解と比較して 10 万回の繰り返し後、約 10％の誤差があったが、概ね良好な結果を得た。
解析モデル
解析対象物
200 [MPa]
50
き裂進展履歴図
10
58,293サイクル後
20,703サイクル後
73,826サイクル後
き裂周辺メッシュ図
3
6
初期き裂
初期亀裂
（半円）
図 6．き裂進展の様子
図 5．解析モデル
亀裂進展長さ履歴
12.0
き裂長さ [mm]
10.0
8.0
a
参照解：a
2c
参照解：2c
6.0
4.0
2.0
0.0
0.00E+00
2.00E+04
4.00E+04
6.00E+04
繰り返し数
8.00E+04
図 7．き裂進展長さの比較
14
数値解析技報 第 10 号, 2014
1.00E+05
４．実用的課題への適用：鋼管継手部からのき裂進展問題（文献 6 参照）
（実験結果との比較検証）
ここでは、鋼管に鞘管が接合された溶接継手を対象に検討する。この継手部のき裂進展実験を参照
して検討を進める。図 8 に鋼管継手の模式図、図 9 に解析モデルを示す。解析モデルはき裂周辺を密
なメッシュサイズとし、接続する区画のメッシュとの間をグル―で接着した。き裂進展解析の初期条
件として、初期き裂を実験で確認されたき裂位置に埋め込み、その形状は深さ 1mm、周長 4mm の半楕
円形状の表面き裂とした。パリスのき裂進展則は３．で適用したものを用いた。荷重の載荷条件は
235kN∼25kN の繰り返し荷重である。
解析結果としてき裂進展の様子を図 10 に示す。また、図 11 に周方向と深さ方向のき裂進展長さの
履歴を示す。図によると周方向のき裂進展長さは、実験と概ね一致している。なお、実測値はき裂深
さ 2ｍｍ（周方向の進展長さ 30ｍｍ）以降の値のみであるため、き裂長さ 30ｍｍ未満については、
実測値に仮の予測ラインを接続して点線で示している。両者のグラフ間に間隙がある理由は、実測値
は初期き裂が微小（原点位置）であったと仮定したためである。ちなみに、図 12 において解析結果
の周長き裂長さを単に零点補正してみた。この試みは、き裂進展速度が初期き裂形状に依存しないと
の仮定に基づいている。図によると、周長方向については、あらためてき裂進展速度が実験と概ね一
致していることがわかる。
WELD
WELD
図 8．繰り返し載荷実験モデル
き裂部詳細
き裂位置
周方向→
40
35
30
25
20
15
10
5
0
5
10
15
20
25
30
35
40
5
着色部：細メッシュ
10
深さ方向↓
15
図 9．詳細メッシュの位置
図 10．き裂進展の様子
18
100
解析（深さ）
16
12
き裂長さ [mm]
実験（深さ）
14
き裂長さ [mm]
亀裂進展長さ（周方向）の履歴
亀裂進展長さ（深さ方向）の履歴
20
鋼 C=7.11×10-13 , m=3.65
10
8
6
4
2
90
解析（周長）
80
実験（周長）
70
実験／外挿（周長）
60
50
40
30
20
10
0
0
0
50000
100000
150000
200000
250000
300000
350000
0
繰り返し数
50000
100000
150000
200000
250000
300000
350000
繰り返し数
図 11．き裂進展長さの比較
数値解析技報 第 10 号, 2014
15
繰り返し数−き裂周長の関係（その２）
き裂長さ [mm]
100
90
解析（零点補正）
80
実験（周長）
70
実験／外挿（周長）
60
50
40
30
20
10
0
0
50000
100000
150000
200000
250000
300000
350000
繰り返し数
図 12．き裂進展長さの比較（その２）
以上 述べてきたように、図 11 から周長方向については き裂長さ 30ｍｍを超える場合は、実測
値と解析結果がオーダー的には合っていると考えられる。一方、深さ方向のき裂進展速度については、
解析結果は実測値の数倍になっており、表現できていない。その理由として、深さ方向については残
留応力の影響等が強く出ている可能性があり、今後の課題と考えている。
５．まとめ
本論文では、Marc を応力解析ソルバーとしたき裂進展解析手法について説明した。今後の課題と
して、解の精度を上げるには計算時間がかかること、また残留応力等の影響を評価する必要があるこ
とが挙げられる。これらの課題に対して、本システムを疲労き裂進展量の予測ツールとして活用しな
がら改善していきたいと考えている。
以上
参考文献
１）岡田裕 他,大規模破壊力学解析のための応力拡大係数計算手法(四面体要素を用いた仮想き裂閉口積分法
(VCCM)),日本機械学会論文集(A 編)73 巻 728 号 (2007-4),PP.498−504
２）岡田裕他,大規模破壊力学解析のための応力拡大係数計算手法(混合モード/複雑形状き裂解析のための四面体
要素用仮想き裂閉口積分法(VCCM))，日本機械学会論文集(A 編)73 巻 733 号(2007-9),PP.997-1004
３）河合浩志 他, 四面体要素向け仮想き裂閉口積分法(VCCM)を用いた
三次元き裂進展解析のためのメッシュ生成技術,日本機械学会論文集(A 編)74 巻 742 号(2008-6),PP.819-826
４ ） 徳 田 貴 志 他 ， 三 次 元 き 裂 進 展 自 動 解 析 シ ス テ ム の 構 築 , 日 本 機 械 学 会 論 文 集 (A 編 )76 巻 770 号
(2010-10),PP.1255-1262
５）岡田 裕 他，三次元き裂進展自動解析システムの構築（第 2 報，き裂進展自動解析システムとき裂進展用有限
要素モデル生成）,日本機械学会論文集(A 編)76 巻 772 号(2010-12),PP.1681-1688
６）Xudong Qian etc.,"FATIGUE BEHAVIOR OF CHS JOINTS FABRICATED USING ENHANCED PARTIAL JOINT
PENETRATION WELDS",Proceedings of the ASME 2011 30th International Conference on Ocean, Offshore and
Arctic Engineering OMAE2011,June 19-24, 2011, Rotterdam, The Netherlands
16
数値解析技報 第 10 号, 2014
技術レポート No. 1
配管網 脈動圧シミュレーション技術の紹介
ソリューション第二チーム
日本機械学会認定
上級アナリスト
∼脈動圧計算プログラムのプロトタイプ構築∼
西村 誠一
概要
配管網が存在するプラント設備では、往復式圧縮機により配管内の流体を移送させることが、
ごく一般に行われる。その際に、圧縮機の周期的外力により配管網の管内に大きな脈動圧（気柱
振動と共振）が発生し、配管が大きく振動する場合があることが知られている。この為、API618
規格(米国石油協会)などを用いて、脈動圧が許容値内であることを検証する必要がある。当社で
はこれまで、この種の検討ツールが無かった為、今回、脈動圧計算の最も一般的な手法である伝
達マトリックス法[1]を用いて、脈動圧計算プログラムのプロトタイプを構築した。そして、機械
学会論文集に記載があった事例を対象に、往復式圧縮機を有する配管系の気柱共鳴を考慮した脈
動圧を計算した所、ほぼ同等な解を得ることができたので、ここに報告する。なお、プログラム
言語としては、エクセル VBA を用いた。
１．伝達マトリックス法
伝達マトリックス法は、昭和４０年代頃からメ
ーカー、大学で開発が進められたものである。こ
の手法は、まずは配管系を管・容量・分岐点など
の要素に分解した上で、各要素を未知数を含む伝
達式で表現する。次に、入口の体積速度脈動（圧
力脈動は未知数）および、出口の圧力脈動、体積
速度脈動を境界条件として、各伝達式で得られた
系を入口から出口までつないでいき、出口の既知
条件から、系の未知数を解くことで、配管系各部
の脈動圧が得られるものである。この基礎式は、
非圧縮性流体理論に基づいている。その為、この
手法は、脈動圧が平均圧に比べて小さく、また平
均流が音響理論で扱える程度の大きさのものに
使用が限定される。
２．圧力脈動の基礎式 （管に対して）
図１に示す一様断面をもつ長さ L の管の基礎
式には、管内流体の密度ρと音速ａ一定値とし
た、式（１）に示す波動方程式が用いられる。
 2 P*
R P *
1  2 P*
 2
 2
2

a ρ t
a t 2
Pi
（１）
Pi1
・
・
Xi
l
ξ
X i 1
また、式（１）の波動方程式は、Ａ、Ｂを任意
定数とした一般式（２）
、
（３）によって表現され
る。


P *  A ' e (α jk )ξ  B ' e (α jk )ξ e jωt
 A cosh(  jk )＋B sinh(  jk ) e jωt
・＊
X
（２）
S

A sinh(  jk )  B cosh(  jk ) e jωt （３）
ρa
・
ここで、P は圧力脈動の振幅、X は体積速度脈
動の振幅、P*は脈動圧、ξは管の座標、R は管の
摩擦係数、t は時間、αは減衰係数、ｋはω/a、
・*
ωは角周波数、 X は体積速度、S は管の断面積
である。
３．各要素の伝達マトリックス式
管路、容量、オリフィスの伝達マトリックス式
は、それぞれ式（４）
、
（５）
、
（６）で表される。
管路
ρa


 Pi 1   cosh(α  jk )l − S sinh(α  jk )l   Pi 
  ・  （４）
・   
 X i 1    S sinh(α  jk )l cosh(α  jk )l  Ｘi 

 ρ a
容量（タンク等）
1
 Pi 1  

ωV
・   
 X i 1    jρ a 2

0  P 
 i 
・
1  Ｘ 
  i 
（５）
オリフィス（バルブ含む）
ρ, a
図１ 管
 Pi 1  
・   
 X i 1  
1
0
ξ  Pi 
・
1  Ｘ 
 i
（６）
ここで、ξは、集中抵抗係数[2]である。
数値解析技報 第 10 号, 2014
17
４．気柱共鳴の確認
図２に示すような閉管（一端開放、一端閉塞）
を対象にして、本プログラムを用いて脈動圧計算
を行った。結果、管の入側圧力の周波数を漸増さ
せると、図３に示すように、気柱共鳴に呼応して
脈動圧が増加し、共鳴振動数（理論値）を明快に
捉えていることが確認出来る。また、図４に示す
ように、管の長さ方向に対する脈動圧分布は、妥
当性のある分布になっていることも確認出来た。
５．適用例
文献[1]に記載されている脈動圧計算をトレー
スした。計算対象は、図５に示すようなアンモニ
ア合成設備の配管系であり、論文中に記載されて
いない条件などは、独自に設備 [3]を調査して仮
定した。
結果、図６に示すように、本プログラムを用い
て、論文の各計測点の脈動圧（振幅）と比較検討
した結果、ほぼ一致していることを確認した。
Ｌ
周波数
を漸増
混合域
窒素
1
2
3
4 5
6
13
16
17
18 21
22 23
20
26
28
29
30
25
24
27
開
1
2
3
4
5
6
7
8
9 閉
12
7
8
9
水素
Point-5の脈動圧
f4=
700Hz
脈動圧（本Progと論文結果の比較）
圧力
0.3
2
|P| kg/cm
600
周波数(Hz)
900
圧力
1次ﾓｰﾄﾞの脈動圧分布
3
1次ﾓｰﾄﾞの定在波形
4
5
6
節点位置
7
8
9
同様
1次ﾓｰﾄﾞの圧力分布
(定在波の節が最大圧力)
図４ 閉管の 1 次モードの脈動圧分布
18 数値解析技報 第 10 号，2014
0.1
0
図３ 閉管の周波数-脈動圧関係
2
0.2
0.15
0.05
気柱共鳴の固有振動数(閉管の場合)
理論値：fn＝(2n-1)v/4L
1次ﾓｰﾄﾞ f1＝100Hz
2次ﾓｰﾄﾞ f2＝300Hz
3次ﾓｰﾄﾞ f3＝500Hz
ここで、v=400m/s、L=1m(管長)
1
本Progの計算結果
論文結果
0.25
f1=
100Hz
300
アンモニア
図５ アンモニア合成工程
f5=
900Hz
f3=
f2=
300Hz 500Hz
0
10 11
19
図２ 閉管
理論値と一致
14 15
15
17
21
23
節点位置
27
30
図６ 各計測点の脈動圧（振幅）
４．おわりに
本プログラムは往復式圧縮機を有する配管系
の脈動圧を求めるものである。現在、そのプロト
タイプが構築できており、今後とも実設備に適用
しながら機能向上に努めたい。
参考文献
[1] 阿部 他、神戸製鋼所：圧縮機配管系の脈動
圧計算法、機械学会論文集、1968.4.30
[2] 毛利 他、東京大学：管路内圧力脈動の共振
振幅、機械学会論文集、1979.8.30
[3] 第 6 版 化学便覧 応用化学編Ⅰ 日本化
学会編 出版社：丸善株式会社
技術レポート No.2
乱流モデルとはなにか
∼乱流モデルの種類とモデルの選択∼
東京分室
春日 悠
概要
乱流モデルは、流体解析において乱流の挙動を決める重要なモデルである。乱流はミクロな現
象を含むため、まじめに計算すると非常に時間がかかり、実用的ではない。そのため、流れ場ご
とに適した乱流モデルが提案されている。本記事では、乱流モデルとはなにか、またどのような
種類がありどう選択するかについて概説する。
1．層流と乱流
流体の流れには層流と乱流がある。層流か乱流
かは、レイノルズ数 Re で判断することができる
層流の速度をどんどん大きくしていくと、どこか
で乱流になる。乱流になるレイノルズ数の値は流
れ場によるが、オーダーで 1000 以上であれば乱
流とみなせる。たとえば、径 1 cm の水道管を水
が 1 m/s で流れる場合、 Re は約 10000 である。
乱流は特別なものではなく、日常的なものであ
る。
2． 乱流モデルの必要性
乱流は、渦が生じては消えていく現象である。
大きな渦は小さな渦にわかれていき、最後には流
体の粘性によって消えてしまう。この過程を厳密
に数値計算することは可能だろうか。
流体解析では領域を格子に分割する必要があ
る。消滅する渦の解像に必要な分割数は、理論的
にレイノルズ数の 9/4 乗のオーダーと見積もら
れる。上で示した水道の例の場合、分割の数は
10 の 9 乗、つまり 1 億である。比較的低速な流
れでもこの程度である。日常的なエンジニアリン
グにおける格子数の規模がせいぜい数百万、がん
ばっても数千万程度であることを考えると、乱流
をまともに数値計算することは現実的ではない。
そこで、乱流のモデル化が必要ということにな
る。
3． 乱流のモデル化
乱流のモデル化には、ナビエ・ストークス方程
式 (NS 方程式) の平均化が用いられる。平均化
にはレイノルズ平均と空間平均がある。
レイノルズ平均は、時間平均やアンサンブル平
均 (個数平均) に相当するものである。流れ場が
時間的になまるため非定常性の再現に弱いが、定
常解析や 2 次元解析が可能なため、比較的計算
リソースが少なくて済む。日常的なエンジニアリ
ングで用いられるのは主にこのタイプのモデル
である。ただし、この方法では普遍的なモデルが
作りにくく、流れ場ごとにモデルを用意する必要
がある。
空間平均を用いた方法としては、Large Eddy
Simulation (LES) がある。フィルターを用いて
流れ場を格子以上のスケールと格子以下のスケ
ールにわけ、格子以上のスケールは直接解き、格
子以下のスケールのみモデル化する。普遍的な性
質が成り立つとされるスケール以下をモデル化
するため、レイノルズ平均モデルよりも普遍性の
高いモデルになる。つまり、高い精度を期待でき
る。ただし、理論的前提により格子をかなり小さ
くしなければならないこと、必ず非定常 3 次元
解析になることにより、かなりの計算リソースが
要求される。日常的なエンジニアリングでの利用
はあまり多くない。
4．RANS モデル
レイノルズ平均を施した NS 方程式をレイノ
ル ズ 平 均 ナ ビ エ ・ ス ト ー ク ス (Reynolds
Averaged Navier-Stokes) 方程式 (RANS 方程
式) と呼び、これをベースにした乱流モデルを
RANS モデルと呼ぶ。NS 方程式と RANS 方程
式の違いは、速度が平均化されているのみなら
ず、レイノルズ応力と呼ばれる項が生じることで
ある。
数値解析技報 第 10 号, 2014 19
流体のせん断応力が分子粘性係数と変形速度
により表されることから、レイノルズ応力を「乱
流粘性係数」と変形速度により表現するのがひと
つの方法である。そうすると、NS 方程式の
RANS 方程式の違いは主に粘性のみとなり、乱
流の効果は渦による「粘っこさ」として現れる。
この方法では、乱れが等方的であると仮定するこ
とになり、乱れの非等方性が現れる旋回流などは
表現しにくい。
輸送モデル (Reynolds Stress Model : RSM) が
ある。本モデルでは乱れの非等方性を考慮できる
ため、 k -  モデルなどよりも旋回流などの再現
性が高い。
ただし、レイノルズ応力は対称テンソルとして
6 つの成分を求める必要がある。さらに  の方程
式も解くため、RANS 方程式に加えて 7 つの方
程式を解くことになり、計算に k -  モデルの 3
倍以上の時間がかかる。
5． k -  モデル
乱流粘性係数を計算する方法として最も広く
使われている方法は、乱流エネルギー k とエネル
ギー散逸率  を用いるものである。RANS 方程
式と k 、 の方程式を解いて乱流を表現するモデ
ルを、一般に k -  モデルという。
標準の k -  モデルは、エンジニアリングにお
いて広く実績のあるモデルであるが、曲率や旋
回、はく離などが苦手である。そのためいくつも
の改良モデルが提案されている。たとえば、流体
解析ソフト ANSYS FLUENT では、以下のよう
なモデルが用意されている。
7．乱流モデルの選択手順
乱流モデルの選択手順は、以下のようになる。
・RNG k -  モデル
繰り込み群 (renormalization group, RNG)
理論を用いたもので、平均ひずみ効果の補正が加
えられており、平均ひずみの大きな流れに有効で
ある。
・Realizable k -  モデル
k などが負の値を取らないように物理的な実
現性 (realizability) の制限を課したモデルであ
る。曲率や旋回がある流れなどに有効とされる。
 の代わりに比散逸率  を用いる k -  モデル
や、それを改良した SST k -  モデルもある。
SST k -  モデルは、比較的旋回流れに強いと言
われている。
また、境界層流れの扱い方にもいくつかモデル
があり、はく離、層流-乱流遷移がある場合や熱
伝達率が問題となる場合などには注意が必要で
ある。
6．レイノルズ応力輸送モデル
レイノルズ応力を乱流粘性係数で表さずに、方
程式を解いて求める方法としてレイノルズ応力
20 数値解析技報 第 10 号，2014
1. 実測値を用意できる場合は、まずそれに合う
ような乱流モデルを計算により探すところから
始める。実測値を用意できない場合は、似たよう
な問題のデータを文献などから入手する。
2. 実測値も文献値も用意できない場合、あるい
はモデルを検討する手続きをとれない場合は、文
献において同様の問題で採用されているモデル
を使うか、過去の実績を参考にする。
3. 以上の手順を取れない場合は、それぞれの乱
流モデルについて一般的に言われている特徴を
参考にしてモデルを選択する。ただし、結果は疑
ってかかる必要がある。
８．おわりに
最後の項目に対する注意として、以下のことが
挙げられる。
・ソフトに実装されているモデルが、一般的にそ
の名前で示されるモデルと同じものとは限らな
い。この場合、一般的な特徴は当てはまらない可
能性がある。
・モデルの一般的な特徴がまったく当てはまらな
い問題もしばしば存在する。旋回流に強いとされ
るモデルが、旋回流において精度が悪いなど。
要するに、やってみなければわからないところ
が乱流モデル (特に RANS モデル) にはある。
参考文献
H.K.Versteeg, W.Malalasekera: 数 値 流 体 力 学 ,
森北出版 (2011)
技術レポート No.3
SPH 法のプログラム開発と解析事例
ソリューション第二チーム
日本機械学会認定
∼SPH 法の基礎理論と２次元平面問題の事例∼
計算力学技術者１級
今道 功治
概要
近年、構造物の数値解析にメッシュに依存しない粒子法が適用される事例が多くなってきた。
当部でも粒子法の技術習得に取り組んでおり、その一環で粒子法（SPH 法）のプログラム開発を
行っている。そこで SPH 法の理論とプログラム開発による解析事例を紹介する。
１．SPH 法の基礎理論
SPH(Smoothed Particle Hydrodynamics)法
は、物体を有限個の粒子（評価点）と呼ば
れる物体素片に分割し、その任意の粒子の
集合を用いた近似化アルゴリズムを持つ。
この手法は，計算格子を必要としないメッ
シュフリー法であるため、大変形領域にお
いても解析対象の変形を表現することがで
きる。
SPH 法と個別要素法（以下、DEM）の違い
は、DEMのように粒子個々に着目して運動方
程式を解くのではなく、連続体内部にラン
ダムに分布した粒子から近似関数の構築を
行うところにある。
図1に連続体を粒子場へ離散化したモデ
ルを示す。粒子i は、ある長さh 内に存在
する他の粒子による相互作用を受ける。こ
の長さhをSmoothing長さと呼ぶ。粒子は物
理量を評価する一種の評価点であり、有限
要素法における節点と積分点の双方の役割
を持つ存在に近い。
図2 物理量の重ね合わせ
この操作により、粒子のもつ物理量は連続
体の物理量に置き換えられる。言い換えれ
ば粒子という物質の孤立したモデルから連
続体モデルへとSmoothing 化される。
(1) Kernel 近似
SPH法の基礎概念は、Kernel近似式にある。
f ( x )   f ( x ) ( x  x ) dx 
ここで、f は３次元位置ベクトルx の関
数であり、dx’は体積、δ（x−x’）はデ
ルタ関数を表わす。デルタ関数をある長さh
をもった適当な重み関数W(x−x’,h)に置
き換えると
f ( x )   f ( x )W ( x  x , h ) dx 
f(x)を   f (x ) と置き、部分積分とGauss
の定理を用いて以下のように表すことが出
来る。
  f ( x )    f ( x )   W ( x  x , h ) dx 
図1 SPH法のイメージ図
ある物理量を重み関数（kernel関数）を
用いて重ね合わせ、ある広がりを持った領
域に滑らかに分布させる概念を図2に示す。
すなわち１階の偏微分項は、重み関数の空
間微分のみで表すことができ、物理量の偏
微分値は微分操作を経ることなく求められ
る。右辺の積分の離散化手法については、
Smoothing長さの内部にある粒子の内挿点
数値解析技報 第 10 号，2014
21
合計により近似する。ここで、重み関数W(x
−x’,h)は、以下のような性質を有する必
要がある。
1) 重み関数の極限がデルタ関数である
lim W ( x  x , h )   ( x  x )
h0
2) 正規化されている
 W ( x, h ) dx  1
3) 影響半径の範囲以外で0 になるコンパ
クト化が可能である
W ( x, h )  0
x h
応力の空間１階偏微分が加速度を表すこと
を意味し、SPH法の基礎式をそのまま適用で
きる。
また、弾塑性解析アルゴリズムは、有限
要素解析で開発された手法を利用すること
が可能である。ただし有限要素法では降伏
判定や弾塑性ひずみマトリクス作成は要素
の積分点で行うのに対し、SPH法では粒子で
行う。
３．プログラム開発による解析事例
SPH 法のプログラム開発において、プログ
ラムフローと解析事例を以下に示す。
(2) 重み関数（kernel関数）
SPH 法では重み関数（kernel関数）を用
いて、あるSmoothing長さｈの内部にある複
数の粒子位置の物理量を中心位置の粒子に
内挿近似する。そのため、この近似精度は、
重み関数の精度に大きく依存される。重み
関数には多数の関数が提案されているが、
よく利用されるのは次の3 次Spline 関数
である。
 C  3 2 3 3
z 1
 h D 1  2 z  4 z 



 C
1 z  2
W 
2  z 3
D
 4h
0
Z 2


ここで、D は問題の次元を表し、Cは定数、
z は次のような代表長さである。
Z
x  x
h
２．運動方程式と弾塑性解析
SPH 法により構造解析を行う場合、個々
の粒子における運動方程式を陽的に解い
て、粒子の運動を求める。この時、粒子の
位置を記述する座標系はLagrange記述とす
る。運動方程式は応力勾配を用いて以下の
ように表される。
a
1

４．おわりに
本プログラムは、2 次元平面弾塑性問題ま
でを計算可能としている。３次元問題、大
変形問題、接触問題などについては、プロ
グラム構築できたら、またここで紹介した
い。
参考文献
「SPH 粒子法の基礎と応用」；横浜国立大学
酒井譲
「粒子法による大変形弾塑性解析」
；横浜国立大学酒井 譲
 
22 数値解析技報 第 10 号，2014
「SPH 法による実大規模斜面掘削実験の二次元再現解
析」；名古屋大学大学院
野々山栄人ら
技術レポート No.4
ホワイトボディ振動解析
ソリューション第二
二ﾁｰﾑ
∼自動車骨格の剛性∼
出口 泰明
概要
MSC.Nastran の振動解析機能を用いて
の振動解析機能を用いて、自動車の骨組（ホワイトボディ）の固有値解析、
の固有値解析
および上下振動を想定した周波数応答解析を行ったので紹介する。また、現在の当社の解
析環境の大規模問題への対応力
析環境の大規模問題への対応力を調査するために、解析モデルの規模を 30 万、60 万、120
万要素の 3 ケースについて
ケースについて固有値解析を行い、計算時間の比較を行ったので併せて紹介す
たので併せて紹介す
る。(解析環境：HP
HP xw8400 Workstation Intel(R) Xeon(R) CPU3.0GHz、MEMORY32
MEMORY32GB)
表１ 解析条件
1．はじめに
近年では、製品開発の場で
製品開発の場で数値シミュレ
解析条件
車体モデル
フォード・トーラス 1)
ーションを活用することにより
ーションを活用することにより開発期間や
要素数
約 30 万要素
コストの削減を行っている。
。自動車メーカ
解析法
Lanczos 法
ーにおいても自動車の低騒音
騒音性を実現する
材料モデル
等方性
ため、開発段階で振動解析により振動特性
により振動特性
重量
約 3000kN
を改善する試みが盛んにおこなわれている
を改善する試みが盛んにおこなわれている。
ここでは、ホワイトボディ
ホワイトボディに対して、以
下の解析を実施した。
表 2 材料条件
骨格
骨格部
材料条件
1) 固有値解析
板厚、スポット径 mm
2) 周波数応答解析
ヤング率 MPa
3
質量密度 ton/mm
2．固有値解析
スポット部
1.2
5
200
200,000
50,000
7.89E
7.89E-09
7.80E-09
0.3
0.3
ポアソン比
先ず、振動的な特性を求めるために
振動的な特性を求めるために実固
有値解析を行った。実固有値解析によって
解析を行った。実固有値解析によって
算出された代表的な固有振動数
固有振動数・固有モー
ド図を図 1 に示す。また、解析
に示す。また、解析条件、材料
条件をそれぞれ表 1、2 に記載する。
記載する。
一般に、ボディー骨格全体
全体で見た振動の
曲げおよび捩り振動モードは
振動モードは 10∼50Hz2)
19.18Hz
25.32Hz
にある。本解析結果でもその
その周波数領域に
あることを確認した。また、支配的なモー
また、支配的なモー
ドは特定できないが、捩れモードと
捩れモードと曲げモ
ードが連成していることがわかる。
28.86Hz
図 1 代表的な固有振動数・モード図
数値解析技報 第 10 号，2014
23
4．当部の解析設備能力
3．周波数応答解析
走行時の上下振動(シェイク)を想定し、図
現モデルは約 30 万要素で構成されてお
2 に示す加振点（赤矢位置）に 周期的強制
り、大半を占めるのは CQUAD4 要素と
変位を与え、モーダル法により減衰を無視
CTRIA3 要素である。そこで、上記の 2 種
した周波数応答解析を行った。観測点は図
類の要素を同要素で分割し約 60 万要素、約
中の青点であり、解析結果を図 3 に示す。
120 万要素、約 240 万要素に規模増とした
図によると図 1 の固有振動数近くで共振
場合の計算時間の比較を行った。表による
が起き、応答倍率が非常に大きくなってい
と Lanczos 法を適用すれば、約 240 万要素
る。また、上記の固有振動数と周波数応答
では約 4 時間かかるとなっている。
解析で求めた共振を起こす周波数では僅か
表4
な差がある。これは、加振点に縦振動を与
解析時間
える際に縦方向には強制変位を与えるとし
要素数
解析時間
て他の 5 つの自由度は全て固定したことが
30 万要素
約 30 分
原因だと考えた。そこで、加振点 4 点の 6
60 万要素
約 1 時間
自由度を固定し実固有値解析を試みたとこ
120 万要素
約 2 時間
ろ、ほぼ共振点に一致する固有周波数を得
240 万要素
3 時間 50 分
ることができた。
5．まとめ
実固有値解析では過去の文献と同じ固有
観測点
周波数の範囲にあるということがわかった。
周波数応答解析では、4 点固定すると固有
周波数が僅かにずれるが、自由モードで得
られた固有周波数と近いということがわか
加振点
図2
った。当部の解析環境で 240 万要素には約
解析モデル
4 時間の計算時間がかかることが確認でき
た。したがって、240 万要素程度の規模で
応答倍率
30
25
あれば、当社の環境で十分対応可能である
20
ことがわかった。
15
参考文献
10
1)http://www.ncac.gwu.edu/vml/models.h
5
tml、Ford Taurus, Modified model
0
18
20
22
24
26
周波数[Hz]
図3
周波数応答
24 数値解析技報 第 10 号，2014
28
30
2)社団法人
自動車技術会編、新編自動車
工学ハンドブック、図書出版社
技術レポート No.５
地盤専用ＦＥＭプログラム「PLAXIS」解析事例紹介
ソリューション第一ﾁｰﾑ
∼地盤解析への取り組み∼
岡本 有造
概要
近年、異常気象による自然災害が多発しており、例えば豪雨による斜面崩壊の対策工などのイン
フラ整備が急務であると痛感している。そこで、それら対策工の開発や実行を数値解析でサポート
するため、当部は地盤解析の技術習得に取り組んでおり、本稿ではその一例として、地盤専用 FEM
解析プログラムである PLAXIS を使った解析事例を紹介する。
１．PLAXIS の概要
（１） PLAXIS とは
PLAXIS はオランダの PLAXIS 社が開発した
地盤解析専用の FEM 解析ソフトウェアで、現在
までに全世界で約 15,000 ライセンスの販売実績
がある。海外、特に東南アジアでは実績および信
頼があり、国内でも PLAXIS のユーザーは海外
事業を展開している企業が多い。
PLAXIS には 2 次元版（2D）と 3 次元版（3D）
のパッケージがあるが、当部では非定常の飽和−
不飽和浸透解析が可能な PLAXIS-2D を新規導
入した。本稿での紹介はこの PLAXIS-2D のみと
なる（次バージョンでは 3D も浸透流解析可能）
。
（２） PLAXIS の特徴
PLAXIS で実施できる計算タイプを表１に示
す。また、PLAXIS に搭載されている地盤の構
成則の一覧を表２に示す。
表１ PLAXIS の計算タイプ
計算タイプ
対象となる事象例
塑性（変形）解析
掘削、斜面崩壊
圧密解析（土−水連成）
盛土、地盤沈下
安全性解析（ｃ-φ低減法） 斜面安定性照査
動的／自由振動解析
構造物の地震時挙動
定常／非定常浸透流解析
堤防、地下水位変動
PLAXIS の特徴に、様々な計算タイプ（表１）
が１パッケージとなっていること、豊富な地盤の
構成則（表２）が搭載されていることがある。
また、図１に示すような段階施工解析の入出力が
可能で、現場での各施工ステージにおける安全性
照査などにも適している。
図１
段階施工解析（地盤掘削の例）
２．解析事例
（１） 斜面の安定解析
簡単な例題として図２に示す斜面の安定性解
析を実施した。地盤は均一とし、各パラメータは
表３に示すとおりとした。
表３ 地盤パラメータ
構成則
表２ 地盤の構成則一覧
線形弾性
モール・クーロン
モール・クーロン
単位重量
kN/m3
18.0
kN/m2
2.0 x104
Hardening soil
HS small
弾性係数
Soft soil
Soft soil creep
ポアソン比
Jointed rock
修正カムクレイ
内部摩擦角
NGI-ADP
Hoek-Brown
粘着力
関口・太田（弾塑性）
関口・太田（弾粘塑性）
ダイレタンシー角
0.3
°
20
kN/m2
10.0
°
0
数値解析技報 第 10 号，2014
25
今回の解析ではまず初期状態として自重解析
を実施した後、せん断強度（ｃ−φ）低減法によ
る安定性解析を行い、斜面の円弧すべりに対する
安全率を算出した。せん断強度低減法とは FEM
における斜面安定評価手法の一つで、地盤のせん
断強度を徐々に低減させることで、破壊を生じさ
せ、その時点の強度低減率の逆数を安全率とする
手法である。一般的に円弧すべりの計算に用いる
極限平衡法ではすべり面を仮定する必要がある
が、せん断強度低減法では、地盤の潜在的な弱点
から破壊が生じ、すべり面を計算上で算定できる
という優位性がある。
今回の例題での計算結果として地盤の変形図
および地すべり対策工として杭を設置したケー
スと無対策のケースの安全率の比較を図３示す。
15m
表４ 地盤パラメータ（フローパラメータ）
サクション水頭 ψ
径 1m、長さ 10m
m
1.5x10-3
図４
有効飽和度 Sr
図４
比透水係数 kr
図４
図４
不飽和水理特性（左：kr-ψ、右：Sr-ψ）
計算結果として、降雨の浸透による地下水位の
変動が分かる間隙水圧分布を図５に、降雨の浸透
による斜面の安全率の変化を図６に示す。
降雨 24h
初期状態
10m
cm/s
飽和透水係数
地下水位
の鋼管杭（3m 間隔）
地下水位はモデル下端
図２
モデル図（杭あり）
図５
円弧すべり
無対策
杭あり
図３
地下水位上昇
飽和状態
降雨 48h
を想定
間隙水圧コンタ
安全率
1.15
1.39
解析結果（無対策の変形図、安全率）
無対策のケースについては、手計算による安全
率とよく一致しており、解析の妥当性が確認でき
た。また、今回のケースでは対策工による安全性
の向上効果は 21%と試算できた。
（２） 浸透流解析との連成（土−水連成）解析
次に、降雨（10mm／h）の地盤への浸透が斜
面の安定性に与える影響を計算した。
対象とした斜面形状は図２と同様（ただし解析
領域は変更）とし、地盤のフローパラメータは表
４とした。今回は不飽和土への降雨の浸透を考慮
する不飽和浸透流解析となるため、一般的な（飽
和）透水係数のほかに、不飽和土の水理特性をモ
デル化する必要がある。PLAXIS には代表的な
地盤について水理特性を規定するフローパラメ
ータが数種類用意されているが、今回は他ソフト
に よ る あ る 例 題 を 再 現 す る た め VanGenuchten モデルによるユーザー設定を行っ
た。設定したフローパラメータを図４に示す。
26 数値解析技報 第 10 号，2014
斜面下端に水平排水孔および抑止杭（水平）を追加
図６
降雨時間と斜面安全率の変化
本解析の妥当性検証は難しいが、例題として参
照した他ソフトによる計算結果とは比較的良い
一致をみており、降雨による斜面安定性の低下や
排水効果を表現できたと考えている。
３．今後の展開
今回報告できていない PLAXIS の動的解析の
機能を活用し、対応できる地盤解析の幅を広げ
る。また、地盤は材料特性が非常に複雑で、材料
パラメータを如何に設定するかが解析の良し悪
しを左右するため、経験を積み、地盤解析技術の
底上げを継続的に図っていく。
４．おわりに
PLAXIS というツールを活用しながら地盤解
析技術を担保し始めており、地盤関連の商品開発
や施工計画のサポートに尽力する所存である。
技術レポート No.6
EXCEL を使用した操業シミュレータ
ソリューション第一ﾁｰﾑ
∼操業状態の可視化∼
長嶋 満也
概要
物流や操業のシミュレータとして世の中に出回っている市販ソフトは、高機能ではあるが
比較的簡単なシミュレーションを行うには高価なものが多い。
そこで、今回身近なツール
である EXCEL を使用して、リサイクル施設の操業シミュレータを開発したので紹介する。
１．シミュレーションの前提条件
本シミュレータはリサイクル施設において車
両にて搬入されてくる不燃物に対して、荷下ろし
や仕分け作業を実施する作業員の負荷バランス
が、どのようになるかを検討するものである。
シミュレーション対象となる物には
① 搬入車両
② 搬出クレーン
③ 場内作業車
④ 置き場（ヤード或いはピット）
⑤ 作業員
の５つが存在し、それぞれ以下のようなパラメ
タが設定可能である。
（１）搬入車両
・時間当たりの搬入台数及び積載量
・作業時間（荷下ろし）及び場内移動時間
（２）搬出クレーン
・搬出時間間隔及び１回の搬出量
（３）場内作業車
・すくい、おろし作業時間
・場内移動時間
（４）置き場（ヤード或いはピット）
・貯留容量、待ち車両台数
（５）作業員
・人数（作業量に反映）
・作業時間（仕分け、検査、荷下ろし、誘導）
・時間当たりの作業量
・場内移動時間
又、共通パラメタとして各々稼働時間を持たせ
ている。
これらのパラメタはマスタテーブルとして保
持し、施設の仕様によって変更可能としている。
搬入車両の台数は実績値を使用するケースと実
績値より求めた日別・時間帯別の平均データを使
用するケースの２通りでシミュレーション可能
としている。
２．プログラムの特徴
本プログラムの最大の特徴は、操業状態を図形
を使用して、物の動きまでを可視化しているとこ
ろにある。 これにより繁忙期の操業状態も視覚
的に捉えることができる。 例えば、場内に車両
が滞留していて侵入できない場合、図１に示すよ
うに場外の待機場で車両が待ち状態となる様子
が確認できる。
図１
車両の待ち状況
さらに、本プログラムの目的である作業員の負
荷状況について、シミュレーション実行中にその
時の稼働率を図２に示すようにグラフ化して表
示しているため、リアルタイムに確認できる。
図２
作業員のﾘｱｲﾙﾀｲﾑ稼働率
数値解析技報 第 10 号，2014
27
この他、ロジック的な工夫も施してある。
以降にその代表的部分を記述する。
（１）場内移動ルートの可変化
車両や作業車が場内を移動する際の移動ルー
トは EXCEL の特徴を活かして図３に示すように、
移動物毎にセルの背景色を関連付けし、その同一
背景色の帯（道）で目的地間を結ぶ事で、移動ル
ートを示している。
従って、この道を施設レイアウトに合わせて変更
する事により、プログラムを変更する事なく簡単
に移動ルート（動き）の変更ができる。
合、アニメーション使用時で、１時間のシミレー
ションには約３分の時間がかかる。
よって、長期間のシミュレーションを実施した場
合に、ある特定時間のものを後から再度確認でき
るように時間単位で状態を記録しておき、これを
再現する機能を有している。
（４）必要パラメタのテーブル（マスタ）化
前提条件で述べたパラメタも含め、シミュレーシ
ョン上で必要とするパラメタの殆どはマスタ値と
してテーブル化されている。
主なテーブルとしては
① 固定物情報
⇒ヤード・ピットの情報
② 移動物情報
⇒車両、作業員の情報
③ ルート情報
⇒移動ルート情報
④ 動作情報
⇒動作に関する情報
⑤ 動作パターン定義情報
⑤は、他のテーブルを参照して、各移動物の動作
パターンを定義しているものである。 動作パター
ンとは“停止”
、
“移動”、
“作業”の基本動作を組み
合わせて一連の動作を定義したものであり、本プロ
図３
背景色で帯状にしたルート
移動物はここで定義された帯上をセル単位で
移動していく事で動きを表現する。
グラムでは非常に重要なテーブルで、例えば施設の
プロフィール変更や動作順序の変更等にも、よほど
特殊な条件でない限りプログラムの変更無しで対
応できる柔軟性を有している。
（２）動作毎に図形を定義可能
本プログラムでは、物の動きを図形により表現
しており、各動作に応じて図を変えられるように
している。 例えば作業員の動きとしては“待
ち”
、
“移動”
、
“作業”の３つがあるがこれを表現
する図形を図４に示す。
３．シミュレーション結果
本プログラムで実行したシミュレーション結
果については、図５の作業員稼働率グラフとか、
図６のタイムチャートの形で確認できる。
図５
待ち
図４
移動
作業員稼働率グラフ
作業
作業員の動作図
この動作毎に図形を定義出来ることにより、さ
らに視覚的にわかりやすいシミュレーションが
表現可能となっている。
なお、この図形についてはユーザが独自で描画・
編集したものを使用できる。
（３）過去のシミュレーション結果の再現機能
本プログラムで、シミュレーションをした場
28 数値解析技報 第 10 号, 2014
図６
タイムチャート
４．おわりに
本プログラムは平面的な動きで表現可能な操
業或いは物流のシミュレーションに幅広く応用
可能なものであるため、同様のニーズを持ってい
る方がおられたら、是非相談して頂きたい。
製品紹介 :
汎用有限要素解析ソフト Marc2013.1
Marc2013.1 の新機能 トピックス
Marcの新バージョン2013.1が2013年10月にリリースされ、より使いやすく、処理速度も速くなりました。
強化された項目の一部を下記にご紹介します。
項目
内容
Mentat モデルナビゲーター機能
モデル内のアイテムをツリー形式で表示
視覚的に使え、作業性も向上。
並列計算機能の拡充
マルチコア環境へ拡充された並列計算機能。
グローバルリメッシングとの併用可能。
新しい材料モデル
異方性ゴム材料モデル。
繊維等を含む非圧縮性材料解析に適用可能。
き裂進展とグローバルリメッシング
グローバルリメッシング機能との併用により
き裂先端の進展と成長の表現を解析可能
◆ Mentat2013.1 モデルナビゲーター機能
解析条件設定にマウス操作でダイレクトにアクセスできるツリー表示機能。
従来よりも作業効率の大幅な向上が期待できます。
◆ 並列計算機能の拡充
従来のマトリクス求解の並列化に加え、
剛性マトリクス作成と応力計算プロセスの並列化が可能となりました。
8 並列で
約 4 倍に
速度アップ
弊 社取 扱 製 品
弊社では構造解析、流体解析などの受託解析に加え、エムエスシーソフトウェア株式会社の代理店として、
Marc 等 CAE ソフトウェアの販売を行っております。
＜お問合せ先＞ソフトウェア契約担当窓口 TEL: 093-588-7233 ／ E-mail：[email protected]
※ 記載の製品名等の固有名詞は各社の商標または登録商標です。
数値解析技報 第 10 号, 2014
29
ＮＳプラント設計製品紹介 :
30 数値解析技報 第 10 号，2014
F.S.Comm
数値解析技報 第 10 号, 2014 31
ＮＳプラント設計製品紹介:
32
数値解析技報 第 10 号，2014
F.S.Addin
数値解析技報 第 10 号, 2014
33
システム化技術紹介 :
34 数値解析技報 第 10 号，2014
光ファイバセンシング
数値解析技報 第 10 号, 2014 35
トピックス 1
【
FOCUS スパコンの利用について
】
熱流体解析チーム 大木 勉、首藤 史
大規模な気泡塔の流動計算を迅速に実施するために、公益財団法人計算科学振興財団（以
下 FOCUS）のスパコンを利用しています。当初は、当部が保有する 120 コアのクラスタ
ーコンピュータで計算を行ってきましたが、検討する項目が多くなったため、今年の夏か
ら FOCUS スパコンを利用しています。
FOCUS スパコンはスーパーコンピュータの産業利用促進を図ることと技術高度化支援、
ならびに企業技術者の育成を推進するために国の補助金によって FOCUS に整備された産
業界専用のスーパ―コンピュータです。独立行政法人理化学研究所が管理する有名なスパ
コン「京」と同系統のコンピュータ仕様になっており、
「京」の建屋に隣接して FOCUS ス
パコンの建屋が設置されています。無料講習会が毎月行われており、スパコン利用拡大に
向けた取り組みがなされています。利用に際しては、利用契約手続きと利用者の登録が必
要ですが、特に手間がかかるようなことはなく、申請から 1 ヶ月弱で利用を開始すること
ができました。利用する CPU 数を増やしたいと思った時には、1 ヶ月の短期利用であれば
Web 上での申請が可能なため、利用者の使い勝手を考えた運用になっています。FOCUS
スパコンの利用は 2013 年 12 月 2 日現在で 118 法人です。NPD の名前も Web 上で公開さ
れております。
当部では FOCUS スパコンの利用に際して先ず、オープンソースの流体解析ソフト
OpemFOAM をインストールして環境設定を行いました。これは Linux に慣れていれば 1
日程度の作業です。計算の実行やデータ処理はインターネットを介して全て戸畑の PC から
行っています。現在、複数の非定常大規模解析を計算しており、数多くの検討を迅速に進
めるうえで FOCUS スパコン大変役立っています。今後も当部のクラスターコンピュータ
だけでは計算機環境が不足するような場合には、FOCUS スパコンを活用していきたいと考
えています。
36 数値解析技報 第 10 号，2014
トピックス 2
【九州デジタルエンジニアリング研究会（KDK）の活動紹介】
九 州デ ジ タル エンジニ アリ ン グ研 究会 （KDK）は九州 地 域近郊の 産学官 連携 によ り、
CAE/CAD/CAM/ＣＧ等
コンピュータ援用技術の向上と、デジタルエンジニアリング技術を駆
使した２１世紀型モノづくり支援への活用を進めるオープンな情報交換の場、新技術習得の場と
して 2004 年に発足しました。弊社は KDK 設立時より主要幹事として会の企画運営に深く関わ
っております。KDK 会員は九州内外の民間企業や工業技術センター等の公的機関、大学・高専
等の教育機でデジタルエンジニアリングに関わられている方々が大部分で、メーリングリスト登
録者数は 300 名を超えています。KDK の会合では、下表に示すような CAE 関わる興味深い講
演、工業技術センターによる中小企業の支援事例、民間企業での適用成果の紹介が行われており
活発な質問や意見が交わされています。
詳細は KDK ホームページをご覧ください。
KDK ホームページ
http://kdk.firebird.jp/
表１
開催日
2011 年
6月
項目
熊本例会
2011 年 10 月
KDK セミナー2011
2012 年 6 月
鹿児島例会
2012 年 10 月
KDK セミナー2012
2013 年 5 月
宮崎例会
2013 年 10 月
KDK セミナー2013
最近の KDK の活動
主な講演
開催場所
「脳とコンピュータを繋ぐ」
熊本大学大学院 自然科学研究科 西山教授
「オープンソース CAE ソフト ADVENTURE の開発秘話
と世界最速のペタスケールコンピュータ"京"のその先」
東洋大学 総合情報学部 塩谷教授
「衛星の恵み “うれしの茶”の誕生秘話」
佐賀大学 知能情報システム学科 新井教授
「北九州市における水素エネルギーに関する取組」
北九州市環境局 平石課長
「車両（2 輪/４輪）開発と部品開発「盗難」と向き合った
部品開発事例」
㈱ホンダロック OB 野田様
「グロ−バル化時代の自動車の開発革新を目指した、考え
方と計算工学の適用」
呉工業高等専門学校 機械工学分野 尾川教授
熊本
数値解析技報 第 10 号, 2014
北九州
鹿児島
北九州
宮崎
北九州
37
トピックス 3
「日本機械学会
第２６回計算力学講演会（佐賀大学にて開催） 出席報告」
- 2013/11/2∼11/4 ‒
フォーラムコーディネータ/座長/論文発表・・・大神勝城
フォーラム パネラー・・・湯本 淳
論文発表・・・首藤 史
【NPD 関係者の講演会での活動】
今回、九州地区での開催となり、企画段階から実行委員として準備を進めました。
NPD から論文発表２件、セッションの座長、フォーラム（
「企業と大学の研究者・実務者が語る V&V
の本音」
）の企画と当日の司会進行、パネラーとして多岐にわたって活動致しました。
その概要は以下の通りです。
（１）論文発表
論文講演者：首藤 史
テーマ：オープンソース CFDEM を用いた CFD-DEM 連成解析について
論文発表者：大神勝城
テーマ：建築鋼構造分野における異種プログラムを連携した統合化解析手法について
（２）セッション 座長
座長者：大神勝城
セッション名：OS27-2.企業での CAE 事例と産官学連携(2)
（３）フォーラム/コーディネーター 、パネラー
フォーラム企画・司会進行 ： 大神勝城
パネラー ： 湯本 淳
テーマ：[ F1. 企業と大学の研究者・実務者が語る V&V の本音]
フォーラムの企画は、吉村忍教授（東大）
、萩原世也教授（佐賀大）と大神（NPD）の３名で企画しま
した。開催に至るまでには、学会関係者、民間技術者の方々からの協力を仰ぐなど慣れない作業を行い
ましたが、無事開催に至ることができました。特に、V&V（Verification &Validation)について、本音を
語らせるのがコーディネーターの役割でしたが、開催時間３時間の長丁場の中で、終始
でて一安心いたしました。このフォーラムの状況については、別途
活発な意見が
機械学会のホームページに掲載さ
れます。
湯本
首藤
（司会：大神）
（備考： 写真は東大 吉村教授 撮影分から一部編集したもの）
38 数値解析技報 第 10 号，2014
【システム開発】
MicroAVSによる可視化
設計・研究開発等の種々の課題に対して、
Marcによる切削解析
数値シミュレーションを軸としつつ、システム開発および
【ソフトウェア販売】
数値解析ソフトウェア販売の連携により、
鋼管の塑性座屈解析
総合シミュレーションエンジニアリングを提供します。
【構造解析】
【熱流体解析】
円形容器のスロッシング解析
【電磁場解析】
誘導加熱
【光ファイバセンシング】
LNGタンクの温度監視システム
MSC.Software 社製品 取扱品目
Marc 、MSC Nastran 、Patran、
University 他
お問合せ窓口
契約担当窓口：
[email protected]
サポート窓口：
[email protected]
※お気軽にご相談ください。
MSC.Software 正規代理店
シミュレーションエンジニアリング・ソリューション部
http://www.npd.eng.nssmc.com/
＜本社＞ 担当 湯本
福岡県北九州市戸畑区中原46-59
TEL：093-588-7233 FAX：093-882-7655
＜東京分室＞ 担当 加来
東京都品川区大崎1-5-1 大崎センタービル
TEL：03-6665-4606 FAX：03-3492-2023
Copyright© 2013 NS Plant Designing Corporation. All rights reserved
NPD 数値解析技報
発
本
行
第 10 号
2014 年 1 月 発行
ＮＳプラント設計株式会社
社 ： 〒804-0002 福岡県北九州市戸畑区中原 46 番地の 59
東京分室 ： 〒141-0032 東京都品川区大崎 1 丁目 5 番 1 号 大崎センタービル 12 階
発行責任者
シミュレーションエンジニアリング・ソリューション部 【SiES’S 部】 湯本 淳
http://www.npd.nsc-eng.co.jp/
TEL： 093-588-7236
FAX：093-882-7655
※本誌掲載の写真・図版・記事などの無断使用・転載・複製を禁じます