板材料創製・成形シミュレーションのための動的陽解法・結晶均質化弾塑性有限要素解析プログラムの開発
大阪工業大学 工学部
機械工学科 教授 仲町英治
経営工学科 助教授 倉前宏行
(平成 15 年度研究開発助成 AF-2003002)
キーワード : 並列処理, 有限要素法, マルチスケール解析,結晶均質化法, 動的陽解法

1. 緒言
金属材料の微視結晶体をミクロン単位で観察した場
合，結晶形態と呼ばれる結晶方位，粒サイズ，粒形状，粒界構造，格子構造，析出物分布など不均質な特性を持つことがわかる．しかし，巨視連続体としての構造体のサイズと微視結晶体のサイズは大きく異なるため，不均質性を正確に反映した微視結晶体から直接巨視連続体の特性を記述する解析手法を展開することは，現在のコンピュータシミュレーションの能力を遥かに超えている．そこで，ユニットセル，あるいは代表体積要素 (Representative Volume Element) と呼ばれる，構造周期性を考慮することで，寸法の階層を上げる，つまりスケールアップの手法として均質化法が提案された．さらに，この均質化法を巨視連続体力学に基づく有限要素法に導入する試みが広く行われるようになってきた1–4) ．均質化法を導入した有限要素解析は，微視結晶体と巨視連続体の連成を満足する．つまり，適合条件および平衡条件を同時に満足できる手法であり，微視結晶体の不均質性を正確に反映できる解析手法といえる．ここで，巨視連続体としては，自動車等の一般機械構造物用の板材を検討対象とする．そこで，板材の強度および成形性を，微視結晶体の初期構造および塑性変形誘起の硬化発展・集合組織発展から予測する結晶均質化有限要素法定式化を行う．
成形性に大きく影響する絞り特性や硬化指数は，初期結晶方位分布，塑性変形誘起の結晶方位回転，つまり集合組織発展，および結晶すべり系で定義される硬化発展と大きく関係する．初期弾性異方性，および塑性域での異方硬化と呼ばれる巨視的特性の異方性と成形性は強く関連している．さらに，塑性変形の進展に伴い，すべり系での転位の運動に差異が生じることでバウシンガー効果と呼ばれる再降伏応力の低下および応力空間内で定義される降伏曲面の移動および変形が起こる5–14) ．そこで，本研究では，従来広く採用されてきた等方硬化則および異方硬化のひとつである移動硬化則を採用した結晶塑性理論に基づく構成式を採用する．
本研究では，結晶均質化・動的陽解法による弾/結晶
塑性有限要素解析を高速化するため，巨視連続体の有限要素メッシュを領域分割することによって並列化を行なう方法を提案する．解析中，弾性域から塑性域へ推移する際，計算コストが大きく変化するため，並列計算を行なう場合には，計算ノード における負荷の均等化を保持しなければならない．そこで，解析中の計算負荷変化に応じて自動的に領域サイズを調整する適応的負荷分散法を導入する．並列性能を評価するため，通信時間を含めた並列解析時間を問題の規模や並列数といった基本量をもとに定量化する．これに基づき，問題の規模や並列数に対する並列性能の依存性について議論するとともに，実際の並列シミュレーションにおいて並列化効率を検討する．

2. 結晶均質化有限要素法
結晶塑性論は，多結晶体の巨視的変形挙動をその構成要素である単結晶体の塑性変形挙動を用いて記述することが可能な理論である．単結晶の塑性変形は，ある特定面上の特定方向に相対的なずれを生じることが知られている．このずれは結晶内に存在する線欠陥である転位の運動に起因し，ずれを生じる面をすべり面，方向をすべり方向と呼び，それらの組合せをすべり系と呼ぶ．結晶
塑性理論は，これらすべり系のすべりを塑性変形の素過程として，単結晶の塑性変形を記述し，さらには多結晶体の巨視的な塑性変形を記述するものである．塑性変形の素過程として結晶のすべりのみを考慮し，弾性変形は原子間距離の変化によって生じると考える．
図 1に示す周期性を満足する微視結晶体を考える．巨視連続体の座標 xi ，微視結晶体の座標 yi を用い，スケール比 λ をxi(1)λ=yiと定義する．

図 1 巨視連続体と微視結晶体
巨視連続体および微視結晶体の速度 U̇i および u̇i をU̇i (x, y) = u̇0i (x) + λu̇1i (y)(2)∂ u̇0i (x)yj + u̇1i (y)u̇i (x, y) =∂xj(3)と定義すると，巨視連続体および微視結晶体に対する仮想仕事の原理式は

ρÜi (x, y)δ U̇i dΩ +Ων U̇i (x, y)δ U̇i dΩΩf¯i δ U̇i (x, y)dΩ +=ΩT̄i δ U̇i (x, y)dΓΓσ−H
σij(x, y)Ω∂δ U̇i (x, y)dΩ∂xj(4)ν u̇i δ u̇i dY = −ρüi δ u̇i dY +YY

σij δ u̇i,j dY(5)YHと表される．ここで，σijは均質化応力であり，微視結晶体で算出されるコーシー応力の体積平均Neσee=1

σH = Ne

(6)

クスとなることから，式 (13) 右辺は，ベクトルの演算のみで算出できる．
結晶均質化弾塑性有限要素解析においては，節点内力
F の算出に必要となる均質化応力算出のため，図 2 に示すように，巨視連続体の要素積分点ごとに微視結晶体の有限要素解析を行なう必要がある．よって本解析は，F の算出時間が支配的となる．さらに，増分時間間隔∆t は，応力波の伝播を厳密に追跡できるように，有限要素の固有周期よりも小さく設定しなければならない．
したがって，増分ステップ数が増大し多大な計算時間を要する．

|Je |e=1により定義される．
微視結晶体に対しては，次に示す速度形構成式を採用した．
eDkl −σ̂ij = Cijkl

(α)

Rij(α)f˙(α) + Fkl Dkl(a) Macro-meshxxxx(b) Gaussian integration points
on a macro-element(α)

V= Cijkl
Dkl −

(α)
Rij f˙(α)(7)(c) Micro-mesh(α)

弾/結晶粘塑性構成式として，Pan-Rice らの提案による結晶粘塑性せん断ひずみ速度 γ̇ (α) と分解せん断応力
τ (α) の関係を採用した．
(α)γ̇ (α) = γ̇0 (α)   (α)  m1 −1τ τ(α)
 g(α) g(8)(α)

ここで，γ̇0 は基準せん断ひずみ速度であり，m はひずみ速度感度指数である．
基準せん断応力 g (α) の硬化発展式を以下に示す．

hαβ γ̇ (β) (9)hαβ = qh(γ) + (1 − q)h(γ)δαβ(10)ġ(α) =

図 2 マルチスケール解析の階層構造
巨視連続体の積分点ごとに行なう微視結晶体の有限要素解析は，完全に独立に実行可能である．したがって，
図 3 に示すように，巨視連続体を領域分割し計算ノード(スレーブ) に割り当て，微視結晶体の有限要素解析による均質化応力の算出を並列化する．マスタから全スレーブに対して変位増分 ∆u をブロード キャスト送信し，領域ごとに算出された均質化応力増分 ∆σ H を受信する．
これをもとに節点内力 F を算出し，工具接触などによる節点外力 P を求めて再び巨視連続体の変位 ut+∆t を更新する．

(β)h(γ) = h0 sech2h0 γ
τs − τ0Mastergenerate F
evaluate P+ ∆tupdate ut

(11)

∆u

ここで，τ0 , τs , h0 , hαβ はそれぞれ，臨界分解せん断応力，飽和分解せん断応力，初期硬化係数，硬化係数を表す．

3. 並列解析の方法∆σ H

subdomain 1

3.1 PC クラスタによる並列計算法
限要素運動方程式は

巨視連続体の有Mü + Cu̇ + F = P(12)のように表せる．ここで，M は質量マト リックス，C は減衰マト リックス，P は節点外力ベクトル，F は節点内力ベクトルである．
中央差分を用いると，式 (12) の運動方程式はdivide into subdomains

xxxx

Slave 1

subdomain n

subdomain 2xxxx

Slave 2xxxx....

Slave n

図 3 PC クラスタを用いた並列解析システム(13)

本研究においては，並列計算環境として，Gigabit Ethernet (1000Base-T, 1000 Mbps) によって接続された 64bit PC (AMD Opteron 1.8 GHz ×2, 2 GB RAM, SuSELinux OS, PGI Workstation 5.1 Fortran 95 Comliler) 9 台(計 18 CPU) を用いる．マスタとスレーブ間におけるデーの通りとなる．この式を解くことにより，時刻 t + ∆t における変位 ut+∆t を求めることができる．ここで，M, Cの両マト リックスは，集中化の操作により対角マト リッ

タ通信 (メッセージパッシング) には，MPICH 1.2.515)(ch p4 デバ イス) を用いる．
3.2 適応動的負荷分散法の導入 弾/塑性解析においては，応力状態が弾性域から塑性域へと変化する際，応力−111CM+ut+∆t =∆t22∆t+M1∆t2





P−F

2ut − ut−∆t + C

1 t−∆tu2∆t評価に要する計算時間が大幅に増加する．このため，並列計算のための領域分割の際，要素数を等しく設定しても，一部の部分領域が集中して塑性しはじめると，計算ノード 間の計算負荷が不均等となり並列性能が大きく低下する．
そこで，各領域の計算負荷の変化に応じて自動的に領域サイズを再設定することにより，計算ノード 間の負荷の均等化を保つようにする．領域サイズを等しく初期設定しておき，各計算ノード における計算時間を測定しながら解析を進めていく．図 4 に示すように，増分計算の反復一定ステップ間隔 (例えば，数百ステップごと) において，それまでの計算時間に応じて領域サイズを再設定する．このとき，巨視連続体の有限要素がスレーブ間を移動するため，対応する微視結晶体の結晶方位などの計算履歴を転送する必要がある．

デルに 8 節点アイソパラメト リック要素の積分点に対応する計測結晶方位を内挿して与えた．図 5 に，8 および1000 有限要素分割および方位カラーマップ 図を比較して示す．

(a) 有限要素モデル(b) 8 要素(c) 1000 要素

図 5 代表体積要素 (RVE)
44 elements

Elasticdeformation

54 elements

79 elements

Elastic + Plastic
deformation

44 elementsPlastic
deformation

Initial Domain
8 subdomains,
8 slaves

44 elements

43 elements

24 elements

44 elements

図 4 負荷均等化のための領域サイズの適応化

初期状態の結晶方位頻度データを用い，RVE の有限要素モデル化の検討を行った．[111] 方向に対する RD お
よび ND 方向との角度差の結晶方位頻度分布を用いた．
9 種類の有限要素モデルの結晶方位頻度分布について，それぞれ母集団と比較したときの誤差平方和を求めた．
有限要素モデルの積分点数と誤差平方和との関係を図 6に示す．これより，1000 積分点以上の有限要素モデルにおいて誤差平方和が収束しており，RVE 有限要素モデルとして成立することがわかった．このときの有限要素数は，125 (= 5 × 5 × 5) 要素である．

4. 微視結晶体
4.1 代表体積要素の検討 結晶方位に関して周期性を満足する 3 次元多結晶代表体積要素 (RVE) を，SEM・EBSD 計測データを用いて導出する．純鉄多結晶板材の試験片を研磨と SEM・EBSD 計測を繰り返しながら，研磨層ごとの結晶方位データを得た．計測解像度は板面内の RD および TD 方向では，3.8 µm，板厚 ND 方向には，5 µm であった．計測範囲は RD 209.55 µm，TD 156.21µm，ND 200 µm であり，この範囲のデータを母集団とした．母集団の方位分布を基準とし方位の周期性を満たす立方体のサイズを誤差平方和を基準に求めた．収束の判定条件として，誤差平方和が 1.0×10−6 を下回ることとすると，1 辺が 100 µm の立方体が RVE として妥当であることがわかった．母集団の平均結晶粒径は 27 µm であることから，この立方体には結晶粒が 50 個存在することがわかった．
実際には塑性変形に伴う転位の蓄積および結晶方位回転によりサブグレインの発生が起こる．そこで，初期および，30 % 引張りひずみ後の結晶方位分布の計測結果と，有限要素解析結果を比較し，必要な微視結晶体有限要素分割数の検討を行った．有限要素分割数を 2×2×2(=8) から 10 × 10 × 10(= 1000) まで変え，それぞれのモ

図 6 初期状態における誤差平方和
4.2 単軸引張り有限要素解析結果 9 種類の有限要素モデルを用いて，単軸引張り解析を行った．純鉄多結晶の材料特性 (ヤング率 210 GPa, ポアソン比 0.3) を用い，式 (9)〜(11) に示す結晶塑性特性値を τ0 = 50.0 MPa,h0 = 80.0 MPa, τs = 80.0 MPa とした．図 7 に要素数と30 % 引張り解析の真応力・真ひずみ線図を示す．
30 % 引張り解析の 1000 要素 RVE の方位カラーマッ
プ 図を図 8 に示す．結晶方位が [101] 方向へ回転することが顕著であることがわかった．このときの結晶方位の図 7 真応力・真ひずみ線図
図 10 30 % ひずみ後の誤差平方和

集積状態を，逆極点図により図 9 に示す．

5. 並列性能の評価

図 8 30 % ひずみ後の代表体積要素

5.1 弾性解析における並列性能の定量化と評価 動的
陽解法の並列計算は，(1) マスタから全スレーブに対して，変位増分ベクトルを送信する，(2) スレーブにおいて担当領域の応力を算出する，(3) スレーブからマスタへ応力を返送する，(4) マスタにおいて外力を算出し変位を更新する，の処理を繰り返し実行する．これらの処理の流れを図 2 に示す．

1st time step
TmMasterSlave 1

2nd time step

Ts Tc

Tb

Slave 2
Slave 3time
図 11 並列計算のタイムチャート

図 9 30 % ひずみ後の逆極点図
4.3 30 % 変形後の代表体積要素 (RVE) の検討 30 %ひずみを与えた試験片の計測解像度は板面内の RD お
よび TD 方向では，3.8 µm，板厚 ND 方向には，4.8 µmであった．計測範囲は，RD 209.55 µm，TD 156.21 µm，ND 177.6 µm であり，この範囲を 30 % ひずみの結晶方位母集団とした．計測結果より平均粒径は 16 µm であり，初期状態より小さいことが確認できた．解析した 9種類の有限要素分割モデルの 30 % ひずみ母集団に対する誤差平方和を求めた．積分点と誤差平方和の関係を図10 に示す．収束判定を行うことにより，2744 積分点以上の分割数により，RVE 有限要素モデルが成立することがわかり，解析に用いる要素数を 343 要素とした．

同図において，方位の頻度分布の相違が大きい部分がみられた．図 9 に示す実験結果の ND 逆極点図においても，[001] 方向に結晶方位の集積が確認される．これらは，解析における単軸引張り境界拘束条件が，計測領域境界と完全に一致しないためであると考えられ，境界近傍の結晶による近距離・遠距離相互作用により，結晶方位回転が影響を受けたためと考える．

巨視連続体の有限要素数，節点数をそれぞれ ne ，npとするとき，ns 台のスレーブを用いると，各処理に要する時間は，それぞれ以下の通り表すことができる．
• 変位増分ベクトルの送信時間
Tb = αb np + βb
• 均質化応力の算出時間 (並列計算部分)
neT s = αsns(14)(15)• 各領域の応力返送時間
Tc = αcne+ βbns(16)• 変位更新時間
Tm = αm np(17)ここで，αb , αs , αc , αm , βb は，演算や通信の速度を表す係数である．これらは簡単な例題を実際に解析し，その際の各処理時間を計測して求めることができる．
以上の処理時間を図 11 に示したタイムチャートに基づいて積算すると，並列解析時間はT = (Tb + Ts + Tc ns + Tm ) It=αsne+ βb ns + (αb np + βb + αc ne + αm np ) Itns

(18)と表すことができる．ここで，It は時間積分のステップ数である．式 (18) の { } 内第 1 項目は，スレーブ数 ns が増すと並列処理の効果によって反比例して減少するが，第 2 項目は比例して増加，第 3 項目以降の ( ) 内 は nsに対して一定値となる．
逐次計算を行なう場合には，通信時間を Tb = Tc = 0とし，スレーブ台数 ns = 1 とおくことにより，解析時間を
T1 = (Ts + Tm ) It= (αs ne + αm np ) It(19)と表すことができる．
したがって，スピード アップ (加速率) Sp や並列効率
Ef といった並列性能は，式 (18), (19) に基づいて
図 13 微視結晶体モデルの規模と並列効率との関係

T1T
T1
Ef =T ns
Sp =(20)(21)Dieと算出することができる．
以上，式 (18) に基づいて 1,000 ステップ分の並列解析(弾性計算) に要する時間を評価した結果を，実際の解析時間と比較して図 12 に示す．評価時間と実際の解析時間は良く一致しており，並列性能の定量化方法は妥当であると言える．

HolderPunch

φ 80(a) 工具設定(b) 巨視連続体の有限要素 メッシュ
図 14 並列性能評価のための解析例題

3000

図 15 に増分ステップ数 10,000 まで解析した際の並列解析時間を示す．ここで，領域サイズの適応化は 200 ス
テップ ごとに計 50 回行なった．本解析は，4,000 ステップ以降，塑性変形が生じはじめるが，適応的負荷分散の効果により，解析が高速化されている．これより，領域サ イズの再設定のための計算および通信コスト はごくわずかであり，並列性能に影響を及ぼさないことが確認される．

2000

1000

0010
Number of Slaves ns20

図 12 評価時間と実行時間との比較
微視結晶体の有限要素モデルを 8 (= 2 × 2 × 2) 要素から 216 (= 6 × 6 × 6) 要素へ変化させたときの並列性能を図 13 に示す．マスター，スレーブ間の通信は，巨視連続体の変位と，微視結晶体において均質化された (体積平均後の) 応力であることから，微視結晶体の有限要素モデルの規模によらず一定である．したがって，有限要素モデルの規模が大きくなるにしたがい，並列計算される計算量のみが大幅に増加するため，並列効率は向上する．
5.2 弾塑性解析例における並列性能評価 並列性能を評価するための例題として，図 14 に示す円筒パンチによる深絞り問題を解析する．巨視連続体，微視結晶体ともに 3 次元 8 節点アイソパラメト リック要素を用いて要素分割した．このときの巨視連続体の有限要素数は 352，微視結晶体は立方体を 2 × 2 × 2 要素分割し，各積分点に結晶方位を与えた．

Parallel Processing Time [h]Time [sec]Estimation
Execution

Adaptive workload allocation
Static workload allocation

1050352 macro−elements,
8 micro−elements

4
Number of Slaves8

図 15 並列解析時間
図 16 に，8 台のスレーブを用いて並列解析を行なった際の，領域サイズの適応履歴を示す．弾塑性計算となる 4,000〜22,000 ステップ の範囲において，計算負荷に応じて領域サイズが適応化されている．また，弾性計算

および全領域が塑性計算となる範囲においては，領域サイズは自動的に均等化されている．
350
Slave 8300

Slave 6200

Slave 5
Slave 4150

Slave 3100

Slave 2
Slave 1500110001

20001
Time Steps

図 18 解析結果 (板厚ひずみ)30001

図 16 スレーブにおける計算負荷量の推移
10,000 ステップまでの解析における並列性能 (加速率)を図 17 に示す．適応的負荷分散を行なっても，スレーブ数が 8 を超えると並列性能が飽和状態となっている．
これは，並列台数の増加に伴って通信コストが増すためであると考えられる．

コード を開発した．解析中の計算ノード 間における計算負荷の均等化を保つため，計算負荷に応じて領域サイズを自動的に調整変更する適応動的負荷分散法を導入し，その効果を確認した．並列性能の定量的な評価，ならびに実際の解析例により，開発した並列化コード の有効性を確認した．純鉄多結晶板材に対し，SEM・EBSD 計測によって得られた結晶方位データをもとに，微視結晶体(代表体積要素) の解像度について検討を行なった．

謝 辞

Ideal Speed−up

16

本解析の一部は，大阪工業大学大学院工学研究科博士前期課程 岡田健二君，山本教範君らにご協力いただいた．ここに，感謝の意を表わす．

Adaptive Workload Allocation
Static Workload Allocation

12
Speed−up

Number of Elements

Slvae 7250

参考文献1) 飯星真治・仲町英治・陳胎平・上田整・上辻靖智: 機論 A, 70-690 (2004), 191–197.

82) 比嘉吉一・澤田幸秀・冨田佳宏: 機論 A, 69-697 (2003),523–529.
403) 松井和己・寺田賢二郎・京谷孝史・岩熊哲夫: 土木工学会論文集, 682/I-56 (2001), 335–346.
4812
Number of Slaves

16図 17 並列性能 (10,000 ステップ まで)

4) Miehe, C., Schotte. J., Lambercht. M.I.: Journal of
the Mechanics and Physics of Solids, 50 (2002), 2123–2167.
5) 池上皓三: 材料, 24 (1975), 491–504.
6) 池上皓三: 材料, 24 (1975), 709–719.

しかし，同じ解析を最終計算 (35,000 ステップ) まで行なうと，塑性計算が支配的となる．したがって，表 1 に示すように，16 スレーブにおいても 80 % と高い並列性能が得られた．このときのシミュレーション結果 (板厚ひずみ) を図 18 に示す．
表 1 並列性能 (最終 35,000 ステップ まで)

7) 五弓勇雄・岸輝雄: 塑性と加工, 107 (1969),863–874.
8) Prager, W.: ASME J. App. Mech. Trans., 78 (1956),
493-496.
9) Ziegler, H.: Quart. Appl. Math., 7 (1959), 55–65.
10) Mroz, Z.: J. Mech. Phys. Solids., 15 (1967),163–175.
11) Harder, J.: Int. J. Plast., 15 (1999), 625–624.

Number ofslaves

Processing
time [h]Speed-up[–]

Paralleleﬃciency [% ]12) Meissonnier, E.P., Busso, E. P. and O Dowd, N.P.:Int. J. Plast., 17 (2001) 601–640.

1816435.12
54.39
33.99

1.0
8.0
12.8

1001008013) Xie, C. L., Groeber, M., Butler, R., Ghosh, S.: Spec.
Publ. Soc. Automot. Eng., SP-1736 (2003), 97–108.

6. 結言
結晶均質化マルチスケール弾/塑性有限要素解析を高速化するため，巨視連続体の領域分割に基づいた並列化

14) Asai, M., Ghosh, S.: Proc. of IWACOM, Japan,(2004), 71.
15) http://www.mcs.anl.gov/mpi/mpich/