マルチスケール法解説(2)

VOXELCONによる解析事例
不確かさのモデリング・シミュレーション法に関する研究会（第2回）
2015.9.10 @慶應義塾大学

株式会社くいんと 月野 誠

会社概要


株式会社くいんと 設立：1985年3月22日
 本社：東京都府中市

 代表取締役会長：石井 惠三、代表取締役社長：月野 誠
 業務内容：計算力学関連ソフトウェアの開発、販売、コンサルティン

グ～夢のあるCAEを日本から～圧倒的に海外製品が優勢なCAEソフトウェア分野において、大学のユニークな研究をもとに、独創的なオリジナル製品で製造業界・社会への貢献を目指しています。
2

くいんとプロダクト
OPTISHAPE-TS
教育用位相最適化ソフトウェア

OPTISHAPE-ES

Surface Generator

構造最適設計ソフトウェア

OPTISHAPE-TS
最適な「かたち」を提案

パラメータ最適化支援ソフトウェア

AMDESS
汎用パラメータ最適化
最適な「値」を提案

イメージベース構造解析ソフトウェア

VOXELCON
現物イメージから計測・解析

3

イメージベース構造解析ソフトウェア

VOXELCON（since 1996）

鋳造欠陥検証

形状比較

構造解析

データ提供：株式会社富士通長野システムエンジニアリング様

静応力解析
定常熱伝導解析
位相最適化

VOXELCON

等値面STLモデル作成
データ提供：群馬産業技術センター様

マルチスケール解析
粒度分布計測、塊・穴のサイズ計測
データ提供：ファインセラミックス技術研究組合/財団法人ファインセラミックスセンター様

均質化法
重合メッシュ法4

VOXELCONの取扱いデータ


データの流れＶＯＸＥＬＣＯＮ

ＣＴスキャン

現物

CT画像
2値化

等値面設計・CAD
3次元計測

ボクセルモデル

解析
分析

ボクセル化
サーフェスモデル（STL）
5

イメージベース解析のはじまり
こんなモデルを作れと言われたら？

1990年代半ば～

6
提供：ミシガン大学 菊池昇 教授

イメージベース解析のはじまり

CT断層画像応力解析

ボクセルモデル（約300万ボクセル）

7
提供：ミシガン大学 菊池昇 教授

イメージベース解析例
デンタルインプラントのイメージベース静応力解析
 高分解能CT画像（15μm）

ボクセル数：8,920,120
CPU : Intel Core i7(3.2GHz)（4並列）
解析時間：2h21m使用メモリ：3,043MB
8
提供：東京歯科大学解剖学講座教授 井出吉信様

イメージベース解析例
恐竜頭骨の解析事例（CT画像から）

提供：東京大学生産技術研究所吉川研究室 大橋智之 様 協力：国立科学博物館 様9CADデータからの解析例


複雑な形状を有するモデルへのボクセル解析の適用
 通常のFEMで解析する際の問題点

そのまま要素分割すると・・・

・要素数が膨大・ハンドリングが困難
・計算コスト大・場合によっては、自動要素分割できない簡略モデルを作成しようとすると・・・
・簡略化の手間がかかる・簡略化が解析精度へ影響

10

CADデータからの解析例
 ボクセル解析法では、誰でも容易に同じボクセルモデルを作成することが可能。

 Robust
 Easy High-speed
CADからエクスポート
VOXELCONに読み込み、ボクセル分割

CAD

STLファイル
ボクセル分割：10秒以下

ボクセルモデル11

CADデータからの解析例
 境界条件を設定して解析した結果（ミーゼス応力）

総ボクセル数：2,313,105
使用メモリ量：777[MB]
計算時間：4時間12分
Xeon E5450（3.0GHz） （4Core使用）

12CT画像と設計データの組み合わせによるモデリング

CT画像とSTL（CAD）を組み合わせたボクセルモデル生成例

大腿骨（CT画像）

ステム（STL）

位置合わせ一体のボクセルモデル生成

13
提供：北里大学病院様ボクセル解析ソルバー
Element by Element PCG法ソルバー
 全体剛性マトリクスを構築しないので、大規模な問題でも省メモリで解析可能。
 ボクセル解析では、すべての要素剛性マトリクスが同じであるため、Element by Elementによる処理に適している。
 反面、全体剛性マトリクスを作成しないため、単純な前処理しか出来ない。
そのため、収束するまでの反復回数が多くなってしまう。

直接法

反復法

60万自由度14

マルチスケール解析～均質化法～
均質化法
 微視的に周期構造を持つ材料の平均的な物性を評価する理論
 マルチスケール解析の一手法

材料定数の均質化

ミクロモデル（ユニットセル）

着目点のひずみ

マクロモデルの解析ミクロモデルの応力解析

 得られる均質化物性値

均質化弾性定数 均質化線膨張率
 均質化熱伝導率
 均質化浸透係数

均質化法
漸近展開法に基づく数学的均質化法を理論背景とし、1980年代半ばより応用数学に通じた工学者の間から、複合材料をはじめとした工学的な問題への応用が始まった。近年では、複合材料や多孔質材料などの工学分野で幅広く用いられている。

15

マルチスケール解析～重合メッシュ法～
重合メッシュ法
 解析領域全体（グローバルメッシュ）に，ローカルメッシュ（詳細に解析したい領域）を重ね合わせて解析 グローバルとローカルが連成

従来のズーミング解析のような、グローバル→ローカルの一方通行ではない。

 モデリングが容易

メッシュの整合性を考えなくてよい。

詳細に解析したい領域（ローカルメッシュ）

重合メッシュ法の発展
T.Belytschko, J.Fishら（1990） ： メッシュ重ね合わせの概念の提案
高野（慶應大）ら（1999～） ： 複合材料構造物への階層モデリングをはじめとしたマルチスケール解析への適用
現在も、研究レベルで破壊力学、非線形解析、振動問題、バイオメカニクスなどへの適用が進められている。

16

重合メッシュ法によるマルチスケール解析例


適用例：鋳造欠陥まわりの応力解析

ボクセルモデル化

グローバルメッシュ（粗）※鋳造欠陥はモデル化されていない

鋳造欠陥※鋳造欠陥が精度よくモデル化されている

ローカルメッシュ（細）17重合メッシュ法によるマルチスケール解析例

解析結果（ミーゼス応力）

グローバルメッシュ：約19万ボクセル
ローカルメッシュ：約58万ボクセル

約2,996万ボクセル

ローカルメッシュの解像度は同程度
計算時間：1時間20分
使用メモリ：3,341[MB] （ICCG法）

計算時間：43時間30分
使用メモリ：6,245[MB] （EBE-PCG法）※Pentium D 3.2GHz

重合メッシュ解析※Core2Quad 2.66GHz（1Core）

計算コストの大幅削減

通常のボクセル解析18

重合メッシュ法と均質化法によるポーラスチタンの強度予測

慶應義塾大学理工学部 高野 直樹 教授のご研究のご紹介

19ポーラスチタン部材の破壊試験
マクロ構造

ミクロ構造

10mm
4mm

気孔径180mm
20mm

マクロ荷重mCT撮像分解能：3.6mm=0.0036mm
気孔径90mm(t=2mm)4mm

60
気孔径180μ m，気孔率70%

マクロ荷重 (N)50

気孔径90μ m，気孔率63%4030
2010000.5

11.5
変位 (mm)2 20

2.5

マイクロCTイメージベースモデリング
1.マイクロCT撮像断面画像

2.三次元構造再現
1504.8mm

二値化処理
ボクセル生成2012.4mm
Y
3.6×599=2156.4mm

z2156.4mm

500 1000 1500 2000
原点からの位置(mm)9590858075706560
55
075
73716967656361595755

500 1000 1500 2000
原点からの位置(mm)xxy断面yz平面zx平面

気孔率(%)

気孔率(%)

気孔率(%)

yz平面

xy断面

0yX

2012.4mm

均質化ユニットセル選択
8078767472706866646260

2012.4mm0500
1000
150021 2000
原点からの距離(mm)均質化法の妥当性確認（１）・材料1
気孔径：90μm
気孔率：63%
解像度：3.6μm

均質化計算により得られたマクロ特性・材料2
気孔径：180μm
気孔率：70%
解像度：3.6μm

ExEy

気孔径90μ m荷重-変位曲線

60Ez

ポアソン比

せん断係数

 zx

Gxy G yz Gzx

 xy

 yz

材料１ 9.41 9.02 8.09 0.224 0.236 0.201 3.72 3.49 3.54
材料２ 6.67 6.21 6.63 0.252 0.233 0.232 2.56 2.56 2.73

気孔径90μ m荷重-変位曲線60

気孔径180μ m荷重-変位曲線
m荷重-変位曲線
気孔径180μ606050505040404040

荷重(N)
荷重 (N)

50303020203030202010101010

気孔径180μ m荷重-変位曲線00600.50.55000)

荷重
荷重(N)
(N)

ヤング率(E,G : GPa)40

111.51.5
変位
変位 (mm)(mm)

222.52.500
00

0.50.5

実験データ
マクロ物性値を用いたFEAデータ

111.51.5
変位
変位 (mm)(mm)

22.522

均質化法の妥当性確認（２）
ひずみ

1.3%

マクロ荷重
Load-Displacement curve (90μ m)

0%

6050

Load (N)4030201000ARAMIS0.5

11.5
2Displacement Y (mm)2.53Experimental data
Elastic finite elemental analtsis data

表面ひずみ比較結果

Y

A1 A2 A3 A4

座標 ARAMIS Simulation Error (%)

B1 B2 B3
C1 C2 C3
D1 D2 D3A4B3C2D1

番号の間隔 0.25mm

X

0.429
0.860
0.587
0.441

0.405
0.833
0.608
0.450

-5.5-3.2
3.5
2.123

重合メッシュ法による強度予測結果
キャリブレーション(気孔径90mm，非線形挙動開始点)

25MPa

Volumefrequency
体積頻度

0.06

200MPa0.05
0.04c

0.03

V 200MPa
microV

micro  0.375

0.02YZ0.01

0200400600

800

Mises stress (MPa)ミクロ相当応力(MPa)マクロ荷重37.6Nにおいてc=0.375を満たす．

40

0.040.030.02

荷重(N)
荷重 (N)

50

X

気孔径180μ m荷重-変位曲線
m荷重-変位曲線
気孔径180μ60

適用結果(気孔径180mm，c=0.375)
気孔径180μ mミクロ応力ヒストグラム

0.05

体積頻度

XY0MPa
重合メッシュ法による解析

00.06

0MPa

強度予測結果30
20
100.01000200400600
ミクロ相当応力 (MPa)

800

00.5

11.5
変位 (mm)(mm)

2242.5

均質化解析事例
～ファインセラミックス～
ミクロモデル（球状孔アルミナ）孔

 解析解像度：4μm
 要素数：2,457,079324mm

 気孔率：7.4%
均質化解析結果（均質化弾性マトリクス）
4.04E+02 1.24E+02 1.23E+02 -5.89E-01 1.24E-01 6.83E-01
4.04E+02 1.23E+02 -6.77E-01 -3.44E-01 7.16E-02
3.95E+02 6.04E-03 -2.35E-01 7.07E-01
1.40E+02 2.44E-01 7.42E-02
1.38E+02 -8.74E-02
1.38E+02

720mm

728mm

CT画像からモデリング

直交異方性を仮定すると、ヤング率(GPa) (Ex,Ey,Ez)=

ポワソン比 (Vxy,Vyz,Vzx)= 0.23537E+00 0.23710E+00 0.23234E+00

せん断弾性率(GPa) (Gxy,Gyz,Gzx) =等価なヤング率は母相の約84%0.34604E+03 0.34556E+03 0.33832E+03

0.14015E+03 0.13781E+03 0.13817E+03

25
提供：ファインセラミックス技術研究組合/財団法人ファインセラミックスセンター様均質化解析事例
～ファインセラミックス～

均質化法による計算結果と実験値との比較

ヤング率(GPa)

410

均質化法

400

超音波パルス法

390

圧縮試験
Hashin-Shtrikman bound

380
370
360
350
340
0246

気孔率(%)
木村圭一、高野直樹、久保弘、小川秋水、河本洋、座古勝：
日本セラミックス協会学術論文集、第110巻、通巻1282号(2002年6月)、pp.567-575「セラミックス多孔体のイメージベースモデリングと均質化法による弾性解析」

8

均質化法で得られる等価物性値が、実験結果ともよく一致している
26

均質化解析事例
～ハニカムサンドイッチパネル～
均質化法を用いた等価材料定数の算定例

 ハニカムサンドイッチパネルによる事例

t方向

寸法

単位mm

Tc=0.0508

b=6.3

Hc=9.8

Tf=0.1

Ls=77.0

Lt=63.0s方向
27均質化解析事例
～ハニカムサンドイッチパネル～

ハニカム部の等価材料定数算出（均質化計算）
等価材料定数ヤング率（E11, E22, E33）
0.30032E-01 0.29434E-01 0.14448E+03

ポアソン比（ν12, ν23, ν31）
0.10096E+01 0.60602E-04 0.30251E+00

せん断弾性率（G12, G23, G31）
0.19717E-01 0.29860E+02 0.19475E+02

ユニットセル28

均質化解析事例
～ハニカムサンドイッチパネル～


解析精度検証簡略モデル
ソリッド要素数：9,600（均質化法を用いない）詳細モデル

等方性材料

シェル要素数：10,080

均質化計算で求めた直交異方性材料●2方向の曲げ解析結果
詳細モデル

簡略モデル

s方向の最大変位

0.638

0.672

t方向の最大変位

1.271

1.333

29

有限被覆法（FCM）への取り組み* 開発中


有限被覆法（Finite Cover Method, FCM）
解析に用いるボクセル
物体（積分領域）：STL
解析のためのボクセルと数値積分の領域を分離

・応力分布の波打ちの大幅な改善
・形状再現性UPによる精度向上
メッシュ分割の簡便性・ロバスト性は維持

30

有限被覆法（FCM）への取り組み
解析例

３次元の円孔モデル（STL）
擬似的な平面応力解析

31

有限被覆法（FCM）への取り組み

従来のボクセル解析

有限被覆法（FCM）

32

まとめ
VOXELCONによるイメージベース解析およびCAD（STL）
データからメッシュ分割について
 均質化法および重合メッシュ法によるマルチスケール解析について、およびその事例 新しい取り組み 有限被覆法（FCM）による従来のボクセル解析の問題点の克服

33