Suuchi Kaiseki Gihou 12 NPD 2017

★技術論文



オープンソース CFD-DEM を用いた球形および

非球形粒子挙動シミュレーション

★技術レポート



進化的計算手法と確率的戦略手法の最適化技術の紹介∼DE・PSO・SA の最適化手法∼



パラメトリック流体形状最適化技術の紹介



粒子を含む流体解析手法について



摩擦熱の解析適用∼SSB（すべり支承）を例として∼



解析によるJ積分の算出、評価について



杭の中抜け現象へのＳＰＨ法適用例題 ∼地盤解析への取り組み（その３）∼



MATLAB Simulink によるプロセス制御シミュレーション技術紹介



製品紹介 : 汎用有限要素解析ソフト Marc2016 新機能紹介



社外情報発信活動

目次

数 値 解 析 技 報 第 １ ２ 号 発 刊 の ご 挨 拶 ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ 1

【技術論文】

○

オープンソース CFD-DEM を用いた球形および

非球形粒子挙動シミュレーション

・・・・・・・・・・・・・・・・

2

【技術レポート】

1.

進化的計算手法と確率的戦略手法の最適化技術の紹介

∼ D E ・ P S O ・ S A

の 最 適 化 手 法 ∼ ・ ・ ・ ・ 6

2.

パ ラ メ ト リ ッ ク 流 体 形 状 最 適 化 技 術 の 紹 介

・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ 8

3.

粒 子 を 含 む 流 体 解 析 手 法 に つ い て ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ 1 0

4.

摩擦熱の解析適用―SSB（すべり支承）を例として・・・・・・・・・・・・・・12

5.

解 析 によ る J積 分の 算出 、 評価 に つい て・・ ・ ・・ ・ ・・ ・・・ ・ ・・ ・ ・・ ・1 4

6.

杭の中抜け現象へのSPH法適用例∼地盤解析への取り組み(その3）・・・・・・・16

7.

M A T L A B S i m u l i n k

に よ る プ ロ セ ス 制 御 シ ミ ュ レ ー シ ョ ン 技 術 紹 介 ・ ・ 1 8

【製品紹介、情報発信】

1.

汎用有限要素解析ソフト

Marc2016

新機能トピックス・・・・・・・・・・・・20

1 数値解析技報 第 12 号 2017

数値解析技報第 12 号発刊のご挨拶

2017 年 01 月

数値解析技報１２号の発刊にあたり、ご挨拶を申し上げます。

この技報では、技術論文として「オープンソースＣＦＤ-ＤＥＭを用いた球形および非球形粒子挙動 シミ ュレ ーショ ン」を 掲載し ており ます。 最近 のトレンドで あり ますオ ープンソース 、連成 を主体 とし、 さらに球形に限定されていたＤＥＭを非球形に開発、解析適用可能なところまで構築しております。

詳細は本文をご覧ください。

技術トピックスとして流体最適化技術の紹介、粒子を考慮した流体解析、摩擦熱の解析適用、非線形J 積分の算出、ＳＰＨ法の地盤への適用、制御シミュレーション技術等を掲載しております。皆様の技術 の一助となればと考えております。

さ て 、 Ｎ Ｓ プ ラ ン ト 設 計 ㈱ シ ミ ュ レ ー シ ョ ン エ ン ジ ニ ア リ ン グ ・ ソ リ ュ ー シ ョ ン 部 （ 以 降 、 SiES’S 部と記す）は、新日鉄住金エンジニアリング株式会社殿および新日鐵住金株式会社殿向け の研究開発、設計支援を中心に社外案件も含め３０年間で約 1 万件超のシミュレーションの実績を 積み上げてきました。構造解析、流体解析、電磁場解析、システム開発の４本柱とともに、互いの 場 を 連 成 さ せた 連 成 解析も 行 っ て き てお り ま す。ま た 、 シ ス テム 開 発 分野で は プ ロ グ ラム 開 発 、 操業データ可視化、データ分析、最適化技術などにも取り組んでおります。組織の陣容は約５０名 を擁し、博士取得者３名、日本機械学会計算力学技術者 認定資格 上級アナリスト５名、固体分野 1 級、熱流体分野 1 級、振動分野 1 級の認定者合わせて１６名を抱え、部員一同 皆様のご愛顧に 応 え る た め に、 数 値 解析、 シ ス テ ム 開発 に 関 する技 術 を 更 に 高度 化 し ており ま す 。 ま た、 東 京 の 新日鉄住金エンジニアリング本社１５階に東京分室を設け、関東、東海地区の営業を担当しており ます。

今後とも、微力ながら皆様のお役に立ちたいと考えております。

ご一読いただき、ご質問・ご意見などありましたら、以下にメールいただければ幸甚です。

メールアドレス：[email protected] シミュレーションエンジニアリング・ソリューション部

オープンソース CFD-DEM を用いた

球形および非球形粒子挙動シミュレーション

当社では流体解析（CFD）と離散要素法（DEM）の連成解析ができるオープンソース CFDEM を公開当初

から利用し，解析時間や精度の面から CFDEM が粉体解析業務に使えるレベルにあることを確認してきた。

本稿ではそのCFDEM を用いた球形および非球形粒子の CFD-DEM 連成解析事例を紹介する。なお、以降に記 載する内容は、設計工学会 2016 年度秋季大会にて著者 首藤が論文発表したものである。

博士（水産科学） 首藤 史

１．はじめに

コンピュータ性能の格段の向上により，プラント設 備 ， 海 洋 ・ 土 木 分 野 な ど に お い て 離 散 要 素 法 （ DEM: Discrete Element Method ） を用 いた シミ ュ レ ーシ ョ ンが活発になってきている．DEM は粒子間の衝突を考 慮した粒子挙動をシミュレーションする手法であり、 連続体シミュレーション（Euler-Euler 法）では表現 が困難な粒状体の自由表面形状や応力ネットワーク構 造を算出できるという特徴をもつ。また，粒子的挙動 と連続体力学を前提にした流体解析（CFD），構造解 析（ FEM）などを組み合わせることで，製造プロセス から防災に関わる海底地盤の挙動検討に至る幅広い分 野での適用が可能となってきている

1） ．

しかしながら，実現象においては，扱う粒子数が億 から兆のオーダーとなる事が多々ある．現在の計算機 でその粒子数のDEM解析を現実的な時間で実行するの は困難であるため，粒径を実際のものより大きくする などして粒子数を減らすことで対応している．それで もDEMの計算には大規模高速計算機が必要となること から，ソルバーの並列化は不可欠と言える．この並列 化への対応としては，各種解析ソフトの並列化ソルバ ー導入が最も手早いと考えられるものの，ソフト導入 費用が高騰する問題があることから，これを解決する ためにオープンソースを用いてソルバーを開発する流 れが盛んになりつつある．

また，DEM では球形粒子を用いることが多いが，実 際の粒子形状は非球形であるものが大半であるため， これら粒子挙動を高精度に把握するには，非球形粒子 での検討が必要になる．

本稿ではまず流体-粒子連成解析に用いたオープン ソース CFDEM

2)

について概説した後，解析例として二 成分粒子の気系流動層における粒子混合と偏析の解析 および，伝熱や反応を伴う流動層の解析を紹介する． 続いて，非球形粒子のDEMについて概要を述べた後， 具体的な適用事例について紹介する．

２．オープンソース CFDEM の概要

CFDEM は CFD-DEM 連成解析を可能としたオープン ソースであり，OpenFOAM

3）

をベースとした CFD ソルバ ーと DEM ソフトウェア LIGGGHTS

4）

で構成されている． CFDEM では CFD-DEM 連成解析手法として

Euler-Lagrange 法と埋め込み境界法の２種類のソルバーが 提供されているが，本章では Euler-Lagrange 法のソ ルバーについて説明する．

図 1に Euler-Lagrange 法の CFDEM の計算手順を示 す．まず CFD 側では DEM 粒子データ (粒子の中心座標 値と速度) の受信後，流体計算セルに粒子を配置する． 次に各流体計算セルにおける流体体積率α

f

(= 1 - 粒子体積率)と粒子−流体間の運動量交換量 (流体力) が計算され，流体力は DEM 側に転送され，流体体積率 と流体力は流体解析に用いられる．DEM 側では，流体 力の受信後，粒子運動を計算し，粒子の中心座標値と 速度を CFD 側へ転送する．ここで，CFD と DEM 間のデ ータ送受信には MPI が用いられている．また，流体力 (抗力，浮力，仮想質量力，揚力)のモデルは既に組み 込まれているモデルを参考にすれば，新たなモデルを 組み込むことも容易である．

図 1 CFDEM の計算手順 CFDEMソルバ

DEM粒子データ取得

セルに粒子を配置

体積率と運動量交換量の計算

DEMにデータを転送(MPI)

流れの計算

LIGGGHTS

CFDにデータ転送

粒子運動の計算

CFDデータの取得 ス タート

３．CFDEM を用いた CFD-DEM 連成解析事例

3.1 二成分粒子流動層の粒子混合と偏析の解析 ここでは文献5)を参考にして実施した，矩形流動層 における二種混合粒子の流動化シミュレーションにつ い て 紹 介 す る ． 槽 の サ イ ズ は 幅 が 0.065 m, 高 さ が 0.26 m, 奥行が 0.0081 m である。総粒子数は 25000 で，粒径 1mm（小粒子）、粒径 2mm（大粒子）がそれ ぞれ重量比 50%ずつの割合で混ざった状態で存在して いる．ガスは底面から一様に送入される系である．シ ミュレーションに用いたパラメータおよび抗力係数モ デルは文献5)と同じである．このシミュレーションは CFDEM が提供しているソルバーをそのまま使用すれば 実施できる．

図2に流入ガス流速1.4と 1.6 m/sの時の粒子流動 状態の時刻歴を示す．まず，ガス流速1.4 m/s では小 粒子が層上部に，大粒子が層下部に偏在していること が確認できる．一方，ガス流速 1.6 m/sになると層全 体が流動化し，層底部を除くと大小粒子が良く混合し ていることがわかる．

図3に実験と解析における粒子流動状態の比較を示 す．実験の方は暗い色が小粒子，明るい色が大粒子で ある．この図から流入ガス流速による大粒子の層高さ の違いが実験と解析で良く一致していることが確認で きる．

図 2 流入ガス流速(a) 1.4 m/s と(b) 1.6 m/s の時の粒子流動状態の時刻歴

図 3 実験と解析の粒子流動状態の比較 （実験の写真は文献5)から転載）

3.2 伝熱を伴う流動層の解析

伝熱を伴う流動層を扱う場合，粒子間と粒子−壁間 の接触および輻射伝熱や粒子−流体間の対流伝熱の考 慮が必要となる．この内，輻射伝熱以外は CFDEMに実 装 さ れ て い る た め 取 り 扱 う こ と が で き る ． ま た ， CFDEM は非圧縮性流体ソルバーのみ提供しているため， 圧縮性流体およびエネルギー方程式のソルバーは新た に作成する必要があるが，これらは OpenFOAM を用い て比較的容易に作成できる．解析例を図 4に示す．

図 4 伝熱を伴う流動層の (a)ガス温度と(b)粒子温度の時刻歴

(a) ガス温度

0.2s 0.5s 1 s 2 s 3 s

(b) 粒子温度

0.2s 0.5s 1 s 2 s 3 s (K)

(a) ガス流速 1.4 ｍ/s 実験 解析

(b) ガス流速 1.6 ｍ/s 実験 解析

(a) ガス流速 1.4 m/s

0 s 2s 5s 10 s 15s 20 s

(b) ガス流速 1.6 m/s

初期状態は粒子温度 400 K,ガス温度 292 K であり， 底部からの流入ガス温度は 292 K である．この図から 時間と共にガスに熱を奪われ粒子温度が低下して行く 様子がわかる．

3.3 反応による粒径変化を伴う流動層の解析 反応を伴う流動層を扱う場合，反応による粒子−流 体間の物質移動やそれに伴う粒子密度や粒径の変化を 考 慮 す る 必 要 が あ る ． こ れ ら を 考 慮 す る に は LIGGGHTS のカスタマイズが 必要にな る．また ，化学 種輸送方程式のソルバーも新たに作成する必要がある．

図 5に反応による粒径変化を伴う流動層の解析例を 示す．これは粒子の揮発分がガスに放出され，時間と 共に粒径が縮小して行く様子を解析したものである．

４．非球形粒子の DEM

4.1 非球形粒子の DEM の概要

非球形粒子の DEM 解析では、個々の粒子の姿勢(物 体固定座標系)を時々刻々更新する必要があるため、 球形粒子の場合と比べて計算量が多くなる。そのため， 球形粒子では数百万粒子が実用的な時間で解析できる ようになってきたが，非球形粒子では未だ困難である．

非球形粒子の形状表現方法には，超二次関数 6)

や超 二次曲面

7)

を用いる方法，多面体で表現する方法 8)

， 球の集合体で表現する方法

9)

がある．球の集合体で表 現する方法は多面体で表現する方法より複雑な形状を 表現し難いことが欠点だが，接触力計算に球形粒子と 同じ方法が使えるため計算量が少ないことが利点であ る．LIGGGHTS には非球形粒子を球の集合体で表現す る方法が実装されているため，その解析例を以下に紹 介する．

4.2 安息角測定試験の解析

ここでは文献10)を参考にして実施した，安息角測 定試験の解析を紹介する．要素モデルは図 5 に示す 6 個の等球形要素を対称配列したものである．球要素の 重複程度によって凹凸度合いを表現するパラメータと して，次式で表す凸度bを用いる

11) ．

Dn = bD (1)

ここで，D_n：非重複領域の長さ，b：凸度を表すパラ メータ，D：球形要素の直径である．安息角測定装置 は文献10)と同様のものとしたが，要素の径は 30 mm （文献10)では 5∼10 mm）としたため，それに合わせ て装置も拡大した．

図 7に凸度b = 0.2と0.5の場合の粒子挙動の時刻歴 を比較して示す．また，図 8に安息角と凸度の関係を 示す．これらの図から凸度bが大きいほど安息角が大 きくなることがわかる．したがって，粒子の凸度bの 調整で安息角を任意の値に調整できるといえる．

図 5 反応による粒径変化を伴う 粒子の流動状態の時刻歴

図 6 六等球形対称配列集合体要素 （文献11)から転載）

図 7 凸度b = 0.2と0.5の粒子挙動の時刻歴比較 （左：凸度b = 0.2，右：凸度b = 0.5） 粒径

(K)

0 s 0.5s 1 s 1.5 s 2s

(a) 基本形状 (b)重複程度

4 s

8 s

12 s

図 8 安息角と凸度bの関係

4.3 埋め込み境界法を用いた CFD−DEM 連成解析 流体と粒子の連成解析手法に埋め込み境界法を用い， 水中に連続的に落下する物体の様子をシミュレーショ ンした結果を図 9に示す。物体は 27 個の球で立方体 を近似したものである。気液界面追跡法には界面捕獲 法（Volume of Fluid 法）を用いている．埋め込み境 界法では，メッシュを粒子サイズより小さくするため， 物体が水中に落下した時に巻き込む気泡まで捉えられ ていることがわかる．

５．おわりに

以上では，オープンソース CFDEM の概要を説明し， 球形粒子と非球形粒子を用いた解析例を紹介した．二 成分粒子の流動層の解析については実験と解析の粒子 流動状態を比較した結果，CFDEM で粒子の混合や偏析 を精度良く予測できることを確認した．したがって， CFDEM は流体−粒子連成解析の業務に利用できるレベ ルにあると考えている．また，非球形粒子を用い，安 息角測定試験を模擬した解析では，凸度bを変えるこ とで粉体の安息角を任意の値に調整できることを確認 した．尚，CFDEM で伝熱や反応を伴う解析や埋め込み 境界法を用いた解析は実行可能だが，今後，実験との 比較等の妥当性検証は必要である．

参考文献

1) Zhu, H.P., Zhou, Z.Y., Yang, R.Y., Yu, A.B.： Discrete particle simulation of particulate systems: A review of major applications and findings, Chemical Engineering Science, 63, (2008), 5728.

2) Goniva, C., Kloss, C., Deen, N.G., Kuipers, J.A.M., Pirker, S.：Influence of rolling friction modelling on single spout fluidized bed simulations”, Particuology, 10, 5, (2012).

図 9 水中に落下する物体のシミュレーション

3) Weller, H. G., Tabor, G., Jasak, H., Fureby, C.：A tensorial approach to computational continuum mechanics using object-oriented techniques, Computers in Physics, 12(6), (1998), 620.

4) Kloss, C., Goniva, C., Hager, A., Amberger, S., Stefan Pirker, S.：Models, algorithms and validation for opensource DEM and CFD-DEM, Progress in Computational Fluid Dynamics, An Int. J. 12, (2012), 140.

5) 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, (2004), 8378.

6) Cleary, P.W.：Industrial particle flow modeling using discrete element method, Eng. Comput. , 26, 6 (2009), 698.

7) Lillie, C., Wriggers, P.：Three-dimensional modeling of discrete particles by

superellipsoids, Proc. Appl. Math. Mech. , 6 (2006), 101.

8) Zhao, D., Nezami, E. G., Hashash, Y. M. A.： Three-dimensional discrete element

simulation for granular materials., Eng. Comput. Int. J. Comput. Aided Eng. Softw., 23, 7 (2006), 749.

9) Bierwisch, C., Kübler, R., Kleer, G., Moseler, M.：Modelling of contact regimes in wire sawing with dissipative particle dynamics, Phil. Trans. R. Soc. A, 369 (2011), 2422. 10) 堀口俊行，澁谷一，香月智：個別要素法による

礫材の安息角解析，構造物の衝撃問題に関する シンポジウム論文集，10 (2011), 43.

11) 堀口俊行，香月智，田附正文：集合体要素を用

いた個別要素法による柔性鋼製枠堰堤の抵抗力 推定法，構造工学論文，59A (2013), 208.

0 s 0.6 s 1.2 s 1.8 s 0

5 10 15 20 25 30 35

0 0.1 0.2 0.3 0.4 0.5 0.6

安息

角

(°

)

6 数値解析技報 第12号,2017 概要

人工知能の一分野である組合せ最適化技術を紹介する。本稿では特に、進化的計算手法である

差分進化(Differential Evolution、以下 DE と略す）と粒子群最適化(Particle Swarm Optimization、

以下 PSOと略す）、確率的戦略手法である焼きなまし法(Simulated Annealing、以下SA と略す）

について簡単に説明し、離散型(Discrete)組合せ最適化の例を紹介する。

1． はじめに

近年、コンピューターの性能が向上したことによ り、人工知能の技術も飛躍的に向上している。その 人 工 知 能 の 一 分 野で あ る最 適 化 技 術を エ ン ジ ニ ア リングに応用している事例も多い。

そこで、当部で最適化技術の獲得に向けた取り組 みを開始し、今後増加することが見込まれる最適化 問題への対応力強化を進めている。

2． 進化的計算手法と確率的戦略手法

進化的計算手法とは、生物の進化や行動を元にコ ンピューター上でシミュレーションを行い、最適解 探 索 を 行 う こ と を特 徴 とし て い る 手法 の こ と で あ る。例として、遺伝的アルゴリズム・差分進化・粒 子群最適化があげられる。

確率的戦略手法とは、最適解探索に確率を用いて 乱択するアルゴリズムのことである。例として、焼 きなまし法がある。

多くの進化的計算手法でも確率を用いており、こ れらの類似性は非常に高いと言える。

3． DE の概要

DE は、突然変異・交叉・適者生存を繰り返しな が ら 大 域 的 最 適 解ま た は準 最 適 解 を求 め る も の で ある。

DEのアルゴリズム (1)探索点の初期化。

探索点はランダムに与える。

(2)全探索点に対して以下の操作を繰り返す。 (2-1) d 番の探索点

d

x

がランダムに3点選ぶ。

(2-2) 突然変異によって、新しい探索点

d

v

を

生成する。

)

(

_{2 3}

1

F

x

x

x

v

_d







Fは突然変異確率を表す。 (2-3) 探索点

d

x

と

d

v

が交叉確率 Cr によって

交叉し、新しい点

d

u

を生成する。

CrossoverPointはランダムに選択。 rは 0から 1の一様乱数を表す。

図1 交叉の例

(2-4) 目的関数の評価をし、探索点を更新。 (3)探索回数を更新。

(4)終了判定を満たしていれば終了し、そうでな ければ(2)に戻る。

(最大世代に到達、解が一定値に収束)

4．PSO の概要

PSOは、生物において一匹がよさそうな場所を発 見すると、群れの残りもそこに集まって周辺を探索 する性質を利用したものである。

PSOのアルゴリズム (1)粒子群の初期化。

各 粒 子 の 位 置 お よ び 速 度 を ラ ン ダ ム に 与 え る。

(2)全粒子に対して以下の操作を繰り返す。 (2-1) 全粒子の評価値を計算する。 (2-2) 個々の粒子の最良位置を保存する。

技術レポート No.１

進化的計算手法と確率的戦略手法の最適化技術の紹介

∼DE・PSO・SA の最適化手法∼

ソリューション第二ﾁｰﾑ

日本機械学会認定計算力学技術者 2級

7 数値解析技報 第12号,2017 これを

i

pBest

とする。

(2-3) 全粒子から最良の位置を保存する。 (評価値の一番いいもの)

これを

gBest

とする。 (2-4) 各粒子の速度を更新する。

)

(

)

(

2 2 1 1 1 i i i i i

x

gBest

c

r

x

pBest

c

r

wv

v













wは慣性力、r

1,r2は0 から 1の一様乱 数、c

1,c2は定数を表す。 (2-5) 各粒子の位置を更新する。

1 1 





i



i i

x

v

x

(3)終了判定を満たしていれば終了し、そうでな ければ(2)に戻る。

(最大反復回数に到達、解が一定値に収束)

5．SAの概要

SA は、金属工学における焼きなましを模したも ので、熱が高い時は大域的に解を探索し、熱が冷め て く る と 大 域 探 索で 得 られ た 最 適 解近 傍 を 探 索 す るものである。

SA のアルゴリズム (1)初期値の設定。

初期温度T、冷却率ε、初期値 x を設定。 (2)以下の操作を繰り返す。

(2-1) 評価値の計算をする。 (2-2) 変更する変数

i

x

を1 つ決定する。

(2-3) 変数の値の増減を決定。



















₀

_.

₅

5

.

0

1

r

x

r

x

x

i i i





αは任意定数、rは0から1の一様乱 数を表す。

(2-3) 前の解と比較し、良質な解の場合は採 用。ただし、一定確率で不良な解も採 用。











 

_



T

E

A

exp

A は解の採用確率、⊿E は今回と前回 の解の差を表す。

(2-4) 温度を冷却する。







 i i

T

T

₁

(3)終了判定を満たしていれば終了し、そうでな ければ(2)に戻る。

(温度が設定温度以下、解が一定値に収束)

6．適用例

歯車の歯数列の設計問題を考える 1)

。下図に示す 4 つ の 歯 車 D,B,A,F が あ り 、 そ れ ら の 歯 数 を x

1,x2,x3,x4とする。

図2 歯車

所 望 の ギ ア 比 と現 状 のギ ア 比 の 差を 最 小 化 す る 問題を考えると、次のように定式化できる。

4

,

3

,

2

,

1

60

12

min

931

.

6

1

)

(

2 4 3 2 1

















_



i

x

x

x

x

x

x

f

i

大域的最適解は総当たりにより得られており、こ の問題では4 つの大域的最適解が存在している。



 













19

,

16

,

43

,

49



43

,

49

,

19

,

16

49

,

43

,

16

,

19

49

,

43

,

19

,

16

,

,

,

_{2 3 4}

1









x

x

x

x

各アルゴリズムで探索した最適解を示す。 表1 探索結果

歯数 DE PSO SA

x

1 16 19 16

x

2 19 16 19

x 3

43 43 43

x

4 49 49 49

f(x) 2.7e-12 2.7e-12 2.7e-12 すべてのアルゴリズムで最適解を得られた。

7．おわりに

最 適 解 探 索 は 必ず 最 適解 が 見 つ かる と い う わ け ではない。しかし、準最適解は見つけることが出来 る。今後もさらに最適化技術を高めていく予定であ る。

参考文献

1) Sandgren, E., Nonliear Integer and Discrete Programming in Mechanical Design Optimization,

ASME/J. of Mechanical

8 数値解析技報 第12号 概要

近年、最適な設計寸法や形状を予測するシミュレーション技術（最適化技術）が注目されてい

る。その要因として、これまでの形状最適化検討は設計者が手動でケーススタディを実施し最適

形状を求めていたが、最適化ソフトの性能向上により、半自動的に検討が可能になったことが挙

げられる。そこで、本稿では熱流体解析ソルバーFluent（ANSYS社）と最適化ソフトHyperStudy

（Altair社）を連成させたパラメトリック流体形状最適化手法について紹介する。

1．はじめに

一 般 に 最 適 化 手 法 は ① ノ ン パ ラ メ ト リ ッ ク 最 適

化、②パラメトリック最適化の2つに大別される。

①は制約条件を設けないため、設計者が予想しない

ような最適形状を得られることがある。そのため、

設 計 初 期 の ア イ デ ア 出 し の 場 面 で 用 い ら れ る こ と

が多いが、得られた最適形状は複雑なため、実際に

は製作できない場合が多々ある。流体形状最適化の

場合、Fluentのみでこの検討が可能である。

一方、②のパラメトリック最適化は、設計者が決

めた制約条件の中で最適形状を見出す手法であり、

詳細設計の段階で用いられることが多い。流体形状

最適化の場合、Fluent および HyperStudy を連成

させて、最適形状を求める必要がある。

今 回 は あ る 制 約 条 件 の 中 で ダ ク ト の 圧 力 損 失 を

最小化することを目的として、角型ダクトのコーナ

ーベー ン形状の最 適化検討 を Fluent と HyperStudyを用いて実施した。

2．連成シミュレーション手法

今 回 の パ ラ メ ト リ ッ ク 流 体 形 状 最 適 化 検 討 の 処

理フローを図1に示す。まず、通常の流体解析の時

と同様にGambitにて形状を定義し、流体解析の設

定を終えた Fluent のケースファイルを作成する。

作成したケースファイルをHyperMeshに読み込ま

せ、変形領域および設計変数の定義を行い、最後に HyperStudyにて、形状最適化手法の設定を行うと

いう流れである。

図1. 処理フロー

3．解析モデル

解析対象を図 2 に示す。角型ダクト内のコーナー

部にベーンと呼ばれる整流板を 3 枚等間隔に配置し

たモデルである。HyperMeshでは2 次元の対象物

を設計変数とすることができないため、解析モデル

は3次元軸対称モデルとした。解析条件を表1 に示

す。

表1. 解析条件

流体解析条件

入口 10m/sの空気が流入

出口 大気圧固定

形状最適化

検討の条件

設計変数

ベーン直線部

（直線部1000mmを初期

状態 とし 、3 枚 の直線

部は同様に変形。）

制約条件

垂 直 方 向 に ±500mm

の範囲内で変形

目的関数 圧力損失の最小化

技術レポート No.2

パラメトリック流体形状最適化技術の紹介

∼Fluent と HyperStudy の連成解析∼

9 数値解析技報 第12号

図2. 解析モデル

4．解析結果

4.1 実験計画法

最適形状の傾向を掴むため、HyperStudyにて実

験計画法を用いた検討（全9水準）を実施した。ベ

ーン直線部は図3に示す通り、1000mmを基準とし

て±500mmの範囲内で、125mmずつ変形する条件

である。実験計画法から得た結果を最小二乗法で近

似し た曲 線を図 4 に示 すが 、ベー ン直 線部長 さが 1000mm 付 近 に ダ ク トの圧 力損 失が 最 小と なる 最

適形状が存在することが確認できる。

図3. 変形領域（ベーン直線部）

図4. ベーン直線部長さと圧力損失の関係（実験計画法）

4.2 ARSM を用いた最適化検討

続いて、前述の実験計画法で得られた結果を用い

て、最適形状の絞り込みを行った。最適形状を特定

す る た め に 今 回 、“ARSM（Adaptive Response Surface Method：逐次更新型応答局面法）”と呼ば

れる最適化手法を採用した。ARSMは最小二乗回帰

法をベースとした最適化手法であり、図4の近似式

か ら ダ ク ト の 圧 力 損 失 が 最 小 と な る 設 計 変 数 の 値

を算出するために用いられる。今回の場合は設計変

数が1つであるため、最小値を見つけることは容易

であるが、設計変数が複数になると、人間の手で求

めることは困難である。そのため、このような場合

は ARSM などの最適化手法を適用することが有益

である。得られた結果を図5に示すが、ARSMを適

用すると、ベーン直線部長さが約 998.6mmの時に

ダクトの圧力損失が最小（12.188Pa）となることが

確認できた。

図5. ベーン直線部長さと圧力損失の関係（ARSM）

5．おわりに

FluentおよびHyperStudyを連成させたパラメト

リック流体形状最適化検討行い、制約条件内での最

適形状を見出すことができた。この手法を実際の設

備に適用することで、様々な流体問題における形状

最適化のニーズに対応できると考えられる。

参考文献

Altair HyperStudyチュートリアル

ガス流れ

直線部（赤色）

が垂直方向に

±500mmの範

囲内で125mm

1 0 数値解析技報 第12号 概要

弊社の主たる解析対象である製鉄プラントや溶融炉およびそれらの付帯設備では、鉄鉱石や石炭、コークス、

微粉炭、ダストといった粒子状物質が多く扱われており、これら粒子の輸送や反応など、粒子を含む流体解析が

多く利用されている。本稿では、代表的な粒子の計算手法および流体解析との相互作用について概説する。

１．はじめに

粒子の計算方法は、流体計算と同様ラグランジュ

型とオイラー型に大別され、問題に応じて適切な計

算手法を選定する必要がある。前者は運動方程式に

も と づ い て 個 々 の粒 子 位置 や 速 度 を直 接 計 算 す る

方法であり、粒子が離散的に存在する場合や粒子同

士 の 衝 突 を 厳 密 に計 算 する 場 合 に 用い ら れ る こ と

が多い。一方、後者は粒子の集まりを流体と同じよ

うに固定した計算格子上で計算する方法であり、流

動 層 の よ う に 粒 子数 が 多く 直 接 計 算す る こ と が 困

難な場合などに適用される。

粒子を含む流体解析では、これら粒子の計算と流

体の計算を連成させて計算する。

２．粒子の計算手法

（１）粒子追跡法

粒子追跡法は、流体解析の入口境界、あるいは内

部領域に通常複数粒子投入点を定め、運動方程式を

用 い て 投 入 点 か らの 各 粒子 位 置 を 追跡 す る 方 法 で

ある。本手法は一般に粒子同士の衝突が無視するこ

とが出来る希薄な状態を想定している。

流 体 解 析 で 得 ら れ た流 れ場 を 固 定 し て 粒 子 の 計

算を行う1-wayの場合には、流体の計算と粒子の計

算を分離することが可能であり、流体解析のいわゆ

るポスト処理の感覚で計算することが出来る。図 1

は、ダストを含むガス配管の計算例を示したもので

あり、粒子追跡法で得られたダストの衝突頻度や速

度、角度情報からNeilsonとGilchristの式 3)4)

を

用いて壁面摩耗速度を評価している。

粒 子 の 濃 度 が 高 く 流れ が粒 子 の 影 響 を 受 け る 場

合 や 微 粉 炭 燃 焼 にお け るガ ス の 揮 発や 粒 子 表 面 反

応、液滴噴霧における液体の揮発・蒸発のように、

粒子と流体との運動量や物質、熱の授受を考慮する

必要がある場合には双方向（2-way）の連成計算

が必要となり、1-way の場合と比べて計算の安定性

や計算負荷の面でハードルが高くなる。

(a)ダストの飛跡 (b)摩耗速度分布

図1 ダストによる摩耗速度の評価例

（２）DEM/離散要素法

DEM も 粒 子 追 跡 手 法 と 同 様 ラ グ ラ ン ジ ュ 型 の 計

算方法であるが、個々の粒子同士の相互作用を計算

する点で大きく異なる。DEMでは互いに接触した粒

子 間 に 働 く 弾 性 反発 力 や摩 擦 な ど の接 触 力 を 、 ば

ね、ダッシュポットおよびスライダーから構成され

るモデルを用いて直接的に計算する。摩擦係数など

の物性値を適切に与えることで、粒子の流動化速度

を直接計算することも可能である。DEM の詳細につ

い て は 本 数 値 解 析技 法 に掲 載 の 技 術論 文 を 参 照 さ

れたい。

（３）二流体モデル

無 数 の 粒 子 が 対 象 であ り且 つ 局 所 的 に 均 質 な 挙

動を仮定できる場合には、計算速度が粒子数に依存

しない二流体モデルの適用が有効である。二流体モ

デルはオイラー型混相流モデルの一種であり、流体

の 計 算 に も 粒 子 の計 算 にも オ イ ラ ー型 の 計 算 手 法

が適用される。

粒 子 が 固体 の場 合 に は 固 体 相 の 内部 摩擦 や 圧 力

を考慮することで、堆積など固体粒子特有の挙動を

技術レポート No. 3

粒子を含む流体解析手法について

ソリューション第二ﾁｰﾑ

日本機械学会認定計算力学技術者１級

塩月 浄志

ソリューション第二ﾁｰﾑ 博士（環境科学）

山田 貴啓

滞留時間 摩耗速度

11 数値解析技報 第 12 号 適切に計算することが可能である。図2 は沈殿槽に

堆積する固体粒子の挙動を計算した例である。(a)

の 固 体 粒 子 間 の 作用 を 考慮 し な い ケー ス で は 体 積

濃度が 100%まで増加しているのに対して、相互作

用を考慮している(b)では、粒子間の空隙を考慮し

た妥当な体積濃度が得られている。

(a)固体の相互作用を考慮しない場合

(b)固体の相互作用を考慮した場合

図2 沈殿槽の計算例

３．流体との相互作用

粒 子 を 含 む 流 体 解 析で は、 粒 子 と 流 体 の 相 互 作

用、すなわち運動量や熱、物質の移動を適切にモデ

ル化する必要がある。物質の移動とは、例えば粒子

から流体側への揮発・溶出、あるいは逆に流体から

粒子への収着である。微粉炭におけるチャー燃焼で

は、流体側の酸素と粒子側の炭素が粒子表面で反応

す る こ と で 二 酸 化炭 素 が生 成 す る とと も に 反 応 熱

も生じることから、両者が複雑に作用している。

粒子と流体の間で交換する代表的な運動量には、

抗力や揚力がある。この他に、乱流場では乱流分散

力、流体や粒子が加速度を持つ場合には仮想質量力

やバセット力といった力が作用する。粒子に作用す

る 流 体 力 を 直 接 流体 計 算で 得 ら れ ない 場 合 に は 何

ら か の モ デ ル 化 によ り これ ら の 力 を計 算 す る 必 要

がある。以下に抗力と揚力について簡単に記す。

（１）抗力

抗 力 を 動 圧 と 粒 子 の投 影面 積 で 無 次 元 化 し て 得

られる抗力係数 C

Dは、よく知られているように粒 子レイノルズ数の関数として整理される。粒子レイ

ノルズ数は、代表寸法に粒子径を、代表速度に粒子

と流体の相対速度を適用して算出する。粒子が球体

の場合には以下に示すSchillerとNaumannの経験

式などを適用して算出する。











44

.

0

/

)

15

.

0

1

(

24

Re

0.687

Re

C

D

1000

1000





Re

Re

実 際 に は 対 象 と す る粒 子が 完 全 な 球 体 で あ る ケ

ースは少なく、粒子の性状に応じて抗力を補正する

必要がある。例えばHaiderとLevenspielの式 5)6)

では、粒子の実際の表面積と、粒子と同じ体積を有

す る 球 の 表 面 積 との 比 をパ ラ メ ー タと し て 非 球 形

粒子の抗力係数を計算する。また、液体粒子の場合

に は 流 動 様 式 に 応じ て 抗力 が 変 化 する こ と を 考 慮

しなければならない。

（２）揚力

揚力は流れに直行する方向に作用する力であり、

流体の速度勾配や粒子の回転によって生じる。野球

で投手が投げる変化球の多くは、後者の力すなわち

ボ ー ル に 回 転 を 与え て 生じ る 揚 力 によ り 球 筋 を 変

化させている。揚力は抗力と比べて小さいため無視

される場合が多いが、非球形粒子などでは揚力の影

響を慎重に評価する必要がある。

４．おわりに

粒子を含む流体解析について概説した。計算機性

能が向上した現在においても、実設備を対象として

粒 子 個 々 の 挙 動 や粒 子 周り の 流 れ の構 造 を 計 算 で

求めることは計算負荷の点で困難な課題である。汎

用的な計算方法は無いため、各種計算手法の中から

そ れ ぞ れ の 前 提 条件 や 限界 を 把 握 した 上 で 計 算 対

象に応じて適切な計算方法を適用する必要がある。

参考文献

1) 日本混相流学会 編集,混相流ハンドブック,朝

倉書店

2) 太田光浩ほか,混相流の数値シミュレーション,

丸善出版

3) NEILSON, J. H.; GILCHRIST, A. Erosion by a

stream of solid particles. Wear, 1968, 11.2:

111-122.

4) NEILSON, J. H.; GILCHRIST, A. An experimental

investigation into aspects of erosion in

rocket motor tail nozzles. Wear, 1968, 11.2:

123-143.

5) ANSYS, ANSYS17.2 FLUENT Theory Guide.

6) HAIDER, A.; LEVENSPIEL, O. Drag coefficient

and terminal velocity of spherical and

nonspherical particles. Powder technology,

1989, 58.1: 63-70.

12 数値解析技報 第12号

概要

鋼製の振り子型免震装置の球面すべり支承における摩擦係数は、免震建物の地震応答に影響を及ぼすことが知ら

れている。摩擦係数は温度と滑り速度の関数 1)2)3)

であることが示されており、また、すべり材をスライダー表面に

接着しているため、温度上昇が接着強度に影響する。従って、地震時の繰り返しにより摩擦熱で生じる温度上昇を

把握することは重要である。摩擦熱を考慮した温度解析は、熱と構造の連成解析で行う必要があるが、構造解析ル

ープで計算時間がかかることが課題である。そこで温度解析のみで温度予測を行う手法を試み、実用的な手法を開

発した。これによりパラメータスタディを効率的に行うことが可能になった。

1． SSB について

SSB装置3)を図１および図2 に示す（建築学会大

会概要より引用）。上コンケーブが建造物に、下コ ンケーブが基礎に固定される。上下コンケーブ間に スライダーは設置され、両者間に滑りが生じて復元 力が発生する。この滑りで摩擦熱が生じる。

２．温度予測 FEM モデルについて

２−１熱構造連成モデル

図3にモデルを示す。コンケーブ・スライダーと もにメッシュを作成する必要がある。球面の R は

4000 ㎜程度と大きく、コンケーブの相対変位に対

し て 生 じ る 復 元 力を 算 出す る 必 要 が無 い の で 平 面 形状としている。コンケーブにGLUE結合した剛 体に強制変位を与えることで相対運動を与える。

接触面の節点に作用する垂直力 F に摩擦係数μ をかけて速度Vをかけると発熱量Qが計算される。

Q=μFV [Nm/s] =μFV [Watt]

この発熱を節点 FLUXとして与えて温度解析を行 う事になる。MARC ではこの換算を自動的に行っ て熱構造連成を実施することが出来る（図4参照）。

２−２ 熱解析モデル

図 5 に示す通り、熱モデルでは、構造解析と同じ メッシュを用いるが、コンケーブの要素タイプは剛 体メッシュとする。構造解析のように別途剛体を設 ける必要はない。

摩 擦 に よ る 発 熱を 境 界条 件 と し て温 度 解 析 モ デ ル に 与 え る こ と で構 造 解析 を 省 く こと が 可 能 に な る。

２−３ 熱解析モデルの発熱の考慮方法

摩擦熱は、両方の接触面に生じるが、接触面が移 動するために、コンケーブ側に発熱を与える領域を

技術レポート No.4

摩擦熱の解析適用

∼SSB（すべり支承）を例として∼

東京分室シニアスタッフ

西 俊二

F V

図3 連成解析モデル 図4 摩擦熱Q

13 数値解析技報 第12号

認識することが困難である。そこで発熱面が一定の スライダー面に発熱 Q を与える。接触面に適当な 熱 伝 達 係 数 を 設 定す る こと で コ ン ケー ブ 側 に も 熱 を与える（図6参照）。

ス ラ イ ダ ー 接 触 面 の垂 直力 分 布 は 一 定 で は な い が、負荷される垂直荷重を面積で割った平均値とし て考える。従って、スライダー面に発生する摩擦を 一定とすることよる誤差が生じることになる。

日本機械学会発行の「伝熱資料」によると、図７ に示す通り接触熱伝達係数Hは接触面圧pにより 変化し、H=１／RH＝0.01[W/(㎜2K)]程度と推定で き る 。 こ の 値 を 参考 に ケー ス ス タ ディ を 行 う 事 と し、本モデルではH=1.0[W/(㎜2K)]とした。

３．熱構造連成モデルと熱モデルの比較

３−１ 要素タイプについて

先に示した様にメッシュ形状は同じであるが、要 素タイプは下表の通りとする。

何 れ も 上 下 コ ンケ ー ブに 反 対 方 向に 強 制 変 位 を 与え運動を再現する。

３−２．計算結果

図 8 にそれぞれのモデルを用いて同一条件で算 出 し た ス ラ イ ダ ーの す べり 面 中 心 と側 面 の ２ 点 で

の温度履歴を比較している。接触面から遠い側面の 結果はほぼ一致していることがわかる。

接触面の結果は、熱解析のほうが少し高めに算出 され安全サイドの評価になっている。

図 9 は最終時刻での温度分布で左が連成解析、右 が熱解析である。

３−３ 計算時間

計算時間は、大凡下表の通りであり、熱構造連成 モデルは熱モデルの8倍程度の時間が必要である。

４．おわりに

熱 構 造 連 成 モ デル を 代用 す る 熱 解析 モ デ ル を 開 発した。スライダー面圧を一定とする仮定を設けて お り 、 誤 差 が 生 じる こ とを 認 識 し た上 で 使 用 す れ ば、パラメータスタディを効率的に実施することが できる。

参考文献；

1) 西本晃治他：新日鉄住金技報,第405号（2016）

p124-128,球面すべり支承NS-SSBの開発

2) 中 村 秀 司 他 ：新 日 鉄 住金 ｴ ﾝ ｼ ﾞ ﾆ ｱ ﾘ ﾝ ｸ ﾞ 技報,Vol.6

（2015）,p28-35,球面すべり支承NS-SSBの開発

3) 中 村 秀 司 他 ： 日 本 建 築 学 会 大 会 学 術 講 演 梗 概

集,2015,p465-466,球 面 す べ り 支 承(SSB)の 繰 り 返しによる温度上昇と摩擦係数の予測

4) 日 本 機 械 学 会 ： 伝 熱 工 学 資 料 （ 改 訂 第 ５

版）,2009,p17

Q =μFV Watt をスライダー側に

図６ 摩擦熱の与え方

図７ 接触熱抵抗RH（熱伝達係数Hの逆数）4)

図８ 温度履歴

14 数値解析技報 第12号

概要

一般的に、構造解析における強度評価指標としては応力が用いられる場合が多い。しかし亀裂先端や 溶接不溶着部のような鋭い切り欠き部においては、理論上、応力集中係数が無限大となってしまうため、 応力による評価は困難である。

破壊力学はこのような問題を取り扱う分野であり、応力に代わる強度指標として、応力拡大係数、エ ネルギー解放率、J 積分などが用いられる。

応力拡大係数とエネルギー解放率は線形弾性材料を前提としているため、金属のような延性材料が広 範囲に塑性降伏する問題には適用できない。これに対し J 積分は非線形材料まで適用することが可能で あり、弾塑性材の破壊力学的評価指標として良く用いられている。本稿では汎用ソルバーMarc による J 積分の算出、評価方法等について紹介する。

1．J 積分の定義

J積分は、き裂先端を含む経路Γ(図 1)に沿った 線積分として下式(1)のように表すことができる。













(

ds

)

x

u

T

Wdy

J

･･･(1)

Γ：き裂先端を囲む積分経路（図 1） W：ひずみエネルギー密度、

T：経路 Γ に沿う表面力ベクトル、 u：経路Γに沿う変位ベクトル、 ds：経路 Γに沿う微小要素

MARC では、(1)式を面積積分に置き換えた(2)式に てJ 値を計算している

1)

。

dA

x

q

W

x

u

J

i

A i

j ij









1

1 1

)

(











･･･(2)

J積分は弾性材料に対して適用されるエネルギー 解放率の考え方を、非線形材料にまで拡張したもの と解釈することができる（よって、弾性材料に対す るJ 積分はエネルギー解放率と同じとなる）。

2．解析によるＪ積分の算出

MARC によるJ 積分算出の例として、破壊靭性評 価の為のASTM規格による3 点曲げ試験

2)

を解析に て再現した(図 2-1)。

① J値の積分経路依存性について

MARC でのJ 積分値は、き裂先端を囲う積分経路(図 2-2)毎の値の平均値として出力される。

荷重点変位 0.4mm 時の塑性域の分布と積分経路毎 のJ 値を図2-3、図2-4 に示す。

技術レポート No.5

解析によるJ積分の算出、評価について

ソリューション第一ﾁｰﾑ

日本機械学会認定計算力学技術者１級

坂田 健

図1 J積分概念図

1 5 数値解析技報 第12号 J値は、き裂先端において若干小さな値をとるが、

おおむね一定しており、J 積分値は経路依存性が無 い事を示している。

②厚さ方向のJ値分布(3 次元効果の影響) 試験片の厚さ方向のJ値の分布を図2-5に示す。J 値は、厚さ方向の中央と端部で差が有り、変位が大 きくなるとその差が拡大している。

3．実験式との比較

J 積 分 は 変 形 に 要 する ポテ ン シ ャ ル エ ネ ル ギ ー の変化を表すものであり、3 点曲げにおける荷重-変位グラフ(図 3-1)より得られる面積 Apl から次の 式(3)で表わす事が出来る。

･･･(3)

K：応力拡大係数、 ν：ポアソン比

E：ヤング率、 B、b：部材寸法

η：形状係数(3点曲げの場合η=1.9)

解析結果より得られるJ値と、荷重-変位グラフお よび(3)式から算出したJ値を比較した。厚さ方向 のJ 値については、平均値、および最大値での比較 を行ったところ、最大値での比較の方がより算出値 に近い値が得られた(図3-2)。

実問題への適用法としては、3点曲げ試験を実施し て破断時のJ 値を算出し、実機モデル解析結果のJ 値と比較する事で亀裂進展有無を判定する、等が考 えられる。

4．おわりに

弾塑性材に亀裂が生じた場合の亀裂進展の評価 手法として、解析によるJ積分算出を紹介した。 今後は業務への適用に向けて、実測データの収集、 解析との比較等、知見の習得に取り組む予定であ る。

参考文献

1)Marc2014 Volume A: Theory and User Information CHAPTER 5 Structural Procedure Library Numerical Evaluation of the J-integral pp.149-151

2) 山 本 琢 也 : 破 壊 靭 性 試 験 法 と デ ー タ 解 析 の 実 例 ( 講 座 : 核 融 合 構 造 材 料に お け る 機 械的 特 性の 評 価手法とデータ解析、2015)

図2-2 MARCの積分経路

図2-3 塑性ひずみ分布

図2-4 積分経路毎の J値分布

図2-5 厚さ方向のJ値分布

図3-1 3 点曲げ解析 荷重-変位グラフ

16 数値解析技報 第12号 概要

ここ数年、本数値解析技報では地盤解析技術の向上に向けた取り組みを報告してきた。今回は

（その３）として杭の中抜け現象などに適用できると考えられるＳＰＨという解析手法への取り組

み例を報告する。

1．ＳＰＨとは

（１）ＳＰＨとは 1)

ＳＰＨ法はメッシュレス（フリーメッシュ）法の

一つに分類される数値解析手法である。ＦＥＭ解析

で定義するような節点と要素の代わりに、「（疑似）

粒子」の集合のみで与えられた物体を表現する。

ただしＳＰＨは、ＤＥＭなどの離散要素法のよう

に互いに衝突や粘着挙動を示す粒子ではなく、偏微

分 方 程 式 の 離 散 化方 法 とし て 粒 子 を用 い て い る だ

けで、基本的にはＦＥＭに近い手法と言える。

（２）SPH モデルの作成

今回用いた解析ソフトでは通常のＦＥＭ（連続体

要素）と同様にモデル作成したのち、計算開始時（計

算途中でも可）にボタン一つでＳＰＨ粒子への変換

が可能であり、モデル作成作業はＤＥＭなどに比べ

ると容易であった。

2．解析対象と解析要領

（１）対象モデル

モデル図を図１に示す。角度20°層厚2mの斜面

を想定し、杭１間隔分のみをモデル化した。すべり

面および杭位置の対称面および「杭」は剛板で表現

した。なお、杭間隔は２ｍ∼４ｍとした。

図１ モデル図

（２） 解析条件

今回想定した斜面（移動層）の材料パラメータを

表１に示す。また、移動層と各剛板には接触条件を

付与し、摩擦係数は表２とした。

作 用 荷 重 は 移 動 層 の自 重の み で 解 析 開 始 直 後 に

作用させたのち、３秒後まで動的解析を行った。

表１ 材料パラメータ

移動層

単位重量 γ 16 kN/m 3

ヤング率 Ｅ 20 x 10 3

kN/m 2

ポアソン比 ν 0.40

内部摩擦角 φ 20°

粘着力 ｃ 10∼50 kN/m 2

表２ 境界部の摩擦係数

接触面 摩擦係数μ

移動層−すべり面 0.01

移動層−杭 0.3

移動層−対称面 0

3． 解析結果

（１）粒子数の影響

本 検 討 で は図 １に 示 す 通 常 の メ ッ シ ュ 作 成 を 行

い、そのメッシュを粒子に変換する。ここでは１メ

ッシュを１次（粒子１個）および２次（粒子４個）

で 変 換 し た 際 の 差異 を 確認 し た 。 対象 は 杭 間 隔 ２

ｍ、移動層の粘着力c=50kN/m2 とした。

解析結果から杭間方向応力（３秒後）を図２に示

す。１次、２次とも杭間や杭上流側で圧縮応力が高

まる状況を捉えているが、２次の方がより明確なコ

ンタでかつ最大値も大きい。当然ながら２次の方が

精度は良いと推定される。

一方、２次の計算時間は１次の約１８倍となり解

析 規 模 に よ り 精 度と 計 算時 間 の バ ラン ス が 課 題 と

技術レポート No.6

杭の中抜け現象へのＳＰＨ法適用例

∼地盤解析への取り組み（その３）∼

ソリューション第一ﾁｰﾑ 日本機械学会認定計算力学技術者 2級

岡本 有造

すべり面（剛板） 対称面（剛板）

※両側に設置 移動層（厚さ2m）

17 数値解析技報 第12号 なることも分かった。なお、今回の以降の検討はす

べて１次変換としている。

図２ 杭間方向応力（左：１次、右：２次）

（２）杭間隔の影響

杭 間 隔 を 変 更 した ケ ース と 杭 な しケ ー ス の Ｘ 方

向（モデル長手方向）の変位コンタを図３に示す。

また、杭間中央点（杭から上流側約１ｍ）の粒子の

速度履歴を図４に示す。

図３ Ｘ方向変位 杭間隔変更（３秒時）

図４ 杭間中央点の粒子の速度履歴

こ れ ら の 結 果 につ い て定 量 的 な 妥当 性 評 価 は 難

しいが、定性的には杭の中抜け現象をある程度表現

できていると考えている。

（３）物性値の影響

移 動 層 の 粘 着 力 c を 変 化 さ せ た 結 果 か ら

c=20kN/m2 のケースのＸ方向変位を図５に示す。

図 ３の 杭 間 隔 ２ ｍ と 比 較 し て 変 位 が 増 大 し て お

り、また、杭近傍と杭間中央部での変位差も大きく

なっていることが確認できた。

図５ Ｘ方向変位 粘着力 c=20kN/m2（３秒時）

4．今後の展開

今回、杭の中抜け現象に対して定性的な評価では

あるが、ＳＰＨの適用可能性を確認できた。今回の

モデルでは杭を剛板としてモデル化しているが、今

後は変形体とすることで、杭の変形や応力状態も把

握できるよう検討を進めていきたい。

5．おわりに

本稿では杭の中抜けに着目した検討を行ったが、

ＳＰＨの技術は汎用性があり、特に変形の大きな領

域 を 対 象 と す る 地盤 解 析に は 有 用 であ る と 考 え て

いる。ＦＥＭに限界を感じた際など積極的にＳＰＨ

を活用して、その可能性を広げていく所存である。

参考文献

1) Abaqus 6.14 マニュアル

0 100 200 300 400 500 600 700 800 900 1000

0 0.5 1 1.5 2 2.5 3

Ｘ

方

向

速

度

（

m

m

/s

e

c

）

時間（sec） 杭間粒子の速度履歴

杭間隔２ｍ 杭間隔３ｍ 杭間隔４ｍ 杭なし

変位量：約0.3m 杭なし

杭間隔２ｍ

杭間隔３ｍ

変位量：約0.9m 変位量：約6.6m

変位量（杭間中央）

：約2.3m 変位量（杭近傍）

18 数値解析技報 第12号 概要

プロセス制御シミュレーションを実現するにあたり、動特性をもつプロセスモデルを MATLAB Simulink上で構築する方法について紹介する。

1． はじめに

プ ロセ ス制 御シ ミュ レー シ ョ ン を 実 現 す る た め

には、図-１に示 すよう にプロセスの入力 （バルブ

操作 etc）に対する出力（温度、圧力 etc）がどの

よう に変 化する かを表 現するプロセスモデルとプ

ロセ スモデル が目標と する出力を 得るため の制 御モデルの２つが必要となる。

図-１ 制御モデル図

プロセスモデルでは、コントローラからの入力を受

け取り、演算結果を出力する。制御モデルは演算値

と目標値を比較し、出力が目標値となるようにコン

トローラを操作する。以下よりプロセスモデルのモ

デル化について述べる。

2．状態方程式のモデル化

動 特 性 を も つ プ ロ セ ス モ デ ル は プ ロ セ ス の状 態

を 示 す 変 数 と そ の時 間 変化 を 表 す 変数 を 用 い た 状

態方程式として表せる。例えばプロセスモデルが一

階微分方程式で表せる場合は、

x’=u ･･･プロセスの時間変化

y=x+C ･･･プロセスの状態

となる。ここで、例題としてディーゼル発電機を冷

却水で冷却するシステムを考える。冷却水は循環し

ており、発電機の熱で温められた冷却水は放熱器で

冷却され発電機へ戻る。放熱器は電動ファンが付い

ており、冷却水温度が一定になるように送風量を制

御する。

図-２ 冷却システムモデル図

表-１ 変数表

パラメータ変数 設定値

Qin 発電機発熱量[W] 5000

Tgoal 冷却水目標温度[℃] 60

Text 外気温度[℃] 20

rho 冷却水密度[kg/L] 1

Cp 冷却水比熱[J/kg・K] 4217

V 冷却水容量[L] 5

A 放熱面積[m

2_{] 2}

T0 冷却水初期温度[℃] 20

計算式

s ファンコントローラ指令値

Alpha(s) 系外への熱伝達率[W/m

2

・K]

=-0.0875・s

2 _{+ 7.875}

・s + 23.75

Qout 放熱量[W]

=A・Alpha(s)・(T-Text)

dQ プロセスの時間変化

=Qin-Qout

Q プロセスの状態（熱量[Ｗ]）

=òdQ・dt+T0・rho・Cp・V

T 出力変数（冷却水温度[℃]）

=Q/(rho・Cp・V)

技術レポート No.7

MATLAB

Simulink による

プロセス制御シミュレーション技術紹介

ソリューション第三チーム

石田 崇

目標値

プロセス

モデル コント

ローラ

入力

出力 −

＋

放

熱

器

発

電

機

Qin T Qout Text

ファン コントローラ

冷却水

Tgoal

−

＋

19 数値解析技報 第12号 3． MATLAB Simulink によるモデリング

図-2に示す冷却システムをSimulinkで構築した

ものが図-３とな る。こ の中で制御モデル に相当す

るものがオレンジで示されたブロックで、プロセス

モ デ ル に 相 当 す る も の が 赤 丸 で マ ー キ ン グ さ れ

た’Heat Model’ブロックである。

図-３ Simulink モデル図

‘Heat Model’は’Level2 MATLAB S-Function’と呼

ばれるAPI*1を通じてSimulink側から呼び出され

るユーザー定義のMATLAB関数である。この機構

を利用することで、制御モデルをSimulink、プロ

セスモデルを MATLAB でそれぞれ構 築する こと

が可能となる。MATLABでプロセスモデルを構築

するメリットは MATLAB の豊富な数学 ライブ ラ

リやソルバを利用できるため、有限差分法や有限要

素 法 な ど を 用 い たよ り 厳密 な プ ロ セス モ デ ル を も

構築可能ということである。

*1 Application Program Interface:相互利用関数のイン

ターフェイス仕様

4．Level2 MATLAB S-Function によるコーディング

’Level2 MATLAB S-Function’では４つの関数を

記述する必要があり、それぞれSimulinkが「ブロ

ック初期化」「状態変数初期化」「微分係数計算」「出

力」を行うタイミングで呼び出される。ユーザーは

呼び出されたタイミングで引数 block に適切な値

を セ ッ ト す る 必 要 が あ る 。 図 − ４ に’Level2

MATLAB S-Function’の コ ー デ ィ ン グ 例 を 示 す 。

(変数及び計算式は表-１参照)。

（つづき）

図-４ ’Level2 MATLAB S-Function’ コーディング例

5．シミュレーション結果

例題のシミュレーション結果を図-５に示す。

シミュレーション開始と共に冷却水温度が上昇し、

180秒から電動ファンが稼動、300秒で目標温度に

到達し、以後温度、指令値共に安定状態になった。

図-５ シミュレーション結果

6. おわりに

今回、プロセスモデリングについて

’Level2 MATLAB S-Function’’を用いた方法を紹

介 し た 。 今 回 の 例 題 規 模 で あ れ ば Simulink の

State-Space ブ ロ ッ ク を 用 い る 方 法 で も 十 分 構 築

可能ではあるが、今後のプロセスシミュレータを作

る 基 盤 技 術 と し て積 極 的に 活 用 し てい く 予 定 で あ

る。

参考文献

Simulink ドキュメント:Level-2 MATLAB S-Function

https://jp.mathworks.com/help/simulink/slref/level

2matlabsfunction.html

%ブロック初期化

functionHeatModel( block )

block.RegBlockMethod('InitializeConditions',

@InitConditions);%状態変数初期化関数の登録 block.RegBlockMethod('Outputs',

@Output);%出力関数の登録

block.RegBlockMethod('Derivatives',

@Derivative);%微分係数計算関数の登録

end

%状態変数初期化関数

function InitConditions(block)

%プロセスの状態変数Qの初期値 Q=T0*rho*Cp*V;

block.ContStates.Data=Q;

end

%微分係数計算関数

function Derivative (block)

%プロセスの時間変化dQ

Alpha=@(s) -0.0875*s^2 + 7.875*s + 23.75; Qin=block.InputPort(2).Data;

s=block.InputPort(1).Data; Qout=A*Alpha(s)*(T-Text); dQ=Qin-Qout;

block.Derivatives.Data=dQ;

end

%出力関数

function Output (block)

%プロセスの状態変数Qと出力変数T Q=block.ContStates.Data

T=Q/(rho*Cp*V);

block.OutputPort(1).Data =T;

end

冷却水温度[℃]

数値解析技報 第12 号, 2017

製品紹介 :

汎

_汎

_汎

用

_用

_用

有

_有

_有

限

_限

_限

要

_要

_要

素

_素

_素

解

_解

_解

析

_析

_析

ソ

_ソ

_ソ

フ

_フ

_フ

ト

_ト

_ト

M

_M

_M

a

_a

_a

r

_r

_r

c

_c

_c

2

₂

₂

0

₀

₀

1

₁

₁

6

₆

₆

Marc2016 の新機能 トピックス

Marcは前バージョンのMarc2015で、従来のGUIであるClassic版が廃止されました。

今後リリースされていくバージョンではすべて新しいGUIとなります。これまでClassic版をお使いのユーザの方で Marcの新しい機能を利用するためには新しいGUIへの切替をお願いします。

◆ 新しくなった GUI

Marc 用プリポストソフトウェア Mentat のメニュー表記が日本語化されました。これから Mentat の利用を始められる方にもわかりや すいインタフェースとなっています。

また解析モデルの設定にダイレクトにアクセスできる「モデルブラウザによるツリー表示」や「操作アイコン」などこれまでに比較して格段に 操作性が向上しています。

Classic 版をお使いのユーザの方は、この機会にぜひ新しい GUI の利用をご検討ください。

モデルブラウザによる

ツリー表示

荷重ケースの設定が Drag&Dropで対応可能に

◆ 相変態機能の追加

MSC ソフトウェアが販売している Simufact 製品からの技術フィードバックとして、Marc の材料データに相変態が設定できるように なりました。

金属は融解温度まで熱せられると相変化と呼ばれる冶金プロセスを受けます。鋼鉄であれば、オーステナイト、フェライト、パーライ ト、ベイナイト、マルテンサイトの分布状態を再現できます。

Marc では、TTT 線図(Time-Temperature-Transformation：時間-温度-変態)、もしくは CCT 線図(Continuous Cooling Transformation：連続冷却変態)を材料データに組み込むことで、材料の硬度と強度の予測が可能になります。

表示操作のアイコン化により、

操作性が向上しました。

図1．連成シミュレーションのフロー

図2．計算例−オーステナイト体積率の表示

20

熱伝導解析 機械的解析