Progresses in Analytical Method of Metal Forming and Material Modeling

Synopsis:

The theory of plasticity was first systematized in a textbook by R. Hill in 1950. This theory was further generalized in the framework of rational mechanics and is nowadays widely used in industrial as well as academic fields by means of the finite-element method. Because of the highly nonlinear process in metal forming, a finite-element simulation of metal forming processes suffers from various theoretical and technical problems. One of them is material modeling, i.e. a plastic constitutive model. Compared to the rapid progress in the analytical method of metal forming, a progress in the material modeling is not significant and a lot of problems still remain. Among them, a modeling of anisotropy is one of the problems of great importance. Recently a crystal plasticity model receives attention again as an alternative method to the conventional phenomenological constitutive model because this model can represent various phenomena including the anisotropy without giving specific assumptions. In this paper, progresses in the analytical method of metal forming and the material modeling during the period of the latter half of the twentieth century until now are comprehensively reviewed. Remaining problems to be solved and future perspectives are also discussed.

1. はじめに

本稿は金属塑性加工の成形解析とそこに適用された塑性構成式の研究に関する20世紀後半から現在までの展望である。

金属成形加工(塑性加工学)をメインフィールドする日本塑性加工学会が発足以来半世紀を越える年月が経過した。設立当初に中心的な役割を果たされた諸先生，諸先輩は日本における塑性加工学の草分け的存在である故福井伸二先生から直接薫陶を受けられた世代である。この方々を第1世代として，筆頭著者は第1世代に属する先生方からご指導を賜った第2世代に属し，その私も10年前に長年勤めた大学を定年退官したので，今や塑性加工学会は第3世代が中心的役割を果たす時代に突入したことになる。

私が塑性力学の研究に関わり始めたのは塑性加工学会が発足した1960年代である。その頃の日本経済は高度成長期にあり，製造技術の中核としての塑性加工への期待も大きかった。また，塑性加工学の土台となる塑性力学も黎明期にあり，産業界の期待を背に新しい学問に携わることの喜びを抱いた多くの研究の徒が参集し，小さいながらも活気あふれる学会であった。やがて高度成長の終焉に相前後して，塑性加工の学問技術も成熟期を迎え，次世代へ向けてのイノベーションを模索しつつ第3世代に突入したのである。

本記事は「鉄と鋼」誌の100巻記念特集号の企画であるが，鉄鋼協会は塑性加工学会の倍の学術活動の歴史を有していることになる。そのような伝統ある学協会誌への展望記事を執筆する機会を与えられたことは筆者二人にとって大変光栄なことと感ずる次第である。

筆頭著者の活動期は主として20世紀後半であり，若手研究者としての共著者は21世紀以降に活躍している。そこでシニア伊藤は主として20世紀後半の塑性力学研究の経過を展望し，世紀末における現状と課題を独断的に総括した上で，共同執筆者には主として21世紀の展望を期待してバトンタッチをする。ただし，それぞれが展望する年代の境界はシームレスとせず重複ないしは片方のみの引用項目や重複項目についての評価が互いに異なることがあり得ることを敢えて認めている。これは時代の変遷に伴う研究トレンドの変化や価値観の変化を反映したものと考え，そのまま読者に提示させていただく。

2. 金属成形解析手法の発展

2・1　20世紀後半の概括的展望

塑性加工過程の解析は塑性力学の問題である。その塑性力学は1950年に発刊されたHill¹⁾の著書によりはじめて体系化された。この書物は第1世代の先生方で翻訳²⁾され，当時は他に類書が無く正に塑性力学のバイブルであったが，その地位は現在までも維持し続けている名著である。改めて読み直してみても，現代流行のマルチスケール的な視点を持ちながら連続体力学の上で見事なまでに体系化されていることが感じ取れる。しかし，この本では350ページに及ぶ訳本の中で一般論は冒頭のわずか70ページ足らずを占めるに過ぎない。その他は多くの塑性加工問題を含む個別問題の解析手法の紹介に費やされている。

それから約10年後，Truesdellらの著書^3,4)を皮切りに有理力学(Rational Mechanics)なる新たな分野が展開され，塑性力学を含む連続体力学の一般化・体系化が盛んに行われた。Archive of Rational Mechanicsというジャーナルも定期刊行されて，関連する研究者の関心を大いに集めたものである。しかし有理力学の指向するところは概念の抽象化であったがために，現実の塑性加工問題に関心のある研究者との間の乖離がむしろ広がった。客観性のある応力速度テンソルや速度型構成式の亜弾性としての性質などが盛んに議論された時期があったが，その熱が冷めた今改めて見直してみると，塑性加工問題の解析に利用されている応力速度の定義も構成式もHill の一般論で使われているものと本質は何も変わっていないことに気が付く。つまるところ，有理力学の恩恵を塑性加工学はいまのところをあまり受けていないようにも思われるが，このような抽象的な一般論が実用の場においても重要視されるようになったのはZienkiewickの著書⁹⁾により，有限要素法が画期的な数値解析手法として登場してからである。一つの解法で関連する問題の殆どを解析できるところの解法の一般化を実現したもので，これにより塑性加工問題の数値解析を行う研究者・技術者は課題ごとに解法を考えてコーディングをする煩わしさから解放されると共に，一つの開発されたソフトが多方面に使えるといったソフト開発の効率の良さをもたらし，コンピューターの驚異的な進歩と相まって急速に普及したのは周知の通りである。有限要素法のもたらした数値解法の一般化・共通化は塑性加工技術に対して功罪半ばであると思われる。それまでは解析することが不可能であった複雑な実加工問題を，殆ど簡略化することなく解析できるようにしたことは大きな進歩である。一方では有限要素はモデルの領域内を基本的に平等に扱うので，ひずみの局所化などの細かい現象の解析精度が低下し，場合によってはそれを予測できないこともある。また，良くできたGUI環境のせいもあって，解析対象固有の性質をよく考察をせずにソフトの丸投げをした結果，有益な情報をそこから引き出せない例もよく見られる。個別問題の特徴を見極めたうえで，適切な要素の選択や解析条件を設定したモデルの解析結果を対象に即して解釈する努力が求められると言う意味で個別問題への回帰も要請される。本来一般性の確保と個別の課題に対する精度保証とは相反する宿命にある。

2・2　成形解析問題の定義と近似解法

成形解析は連続体力学としての塑性力学の問題であり，次のような微分方程式の境界値問題として定義される。対象とする物体の表面Sの一部S_F上で表面力t，残りの表面S_Uで変位uが指定されているとき   

$$\begin{array}{c}力の境界条件:\sigma \cdot n=t\text{on}{\text{S}}_{\text{F}}\\ 変位拘束条件:u=U\text{on}{\text{S}}_{\text{U}}\end{array}\}$$ (1)

なる境界条件を満たす次式の解を求めるのが成形解析問題である。   

$$\begin{array}{c}平衡方程式:\nabla \cdot \sigma =0\\ 構成式:d\epsilon =Cd\sigma \end{array}\}$$ (2)

ここにnは物体表面外向き法線ベクトルであり，Uは物体表面で指定された変位の値である。この境界値問題で中心となる場の方程式(平衡方程式)は偏微分方程式であり，一般に正解を求めることが困難である。そこで何らかの近似解で妥協せざるを得ない。成形解析の発展の歴史は近似解法の発展の歴史でもある。

Hillの著書に掲載されている塑性加工問題の解析事例はコンピュータが実用になる以前のものであり，殆どが閉じた式表現による解析解である。これらは解表現を得るために対象となる現象を大胆に単純化している。その意味で解析解の殆どは近似解である。ここで，境界値問題の近似解の分類とそれらの性格について整理をしておく。

式(1)，(2)の全てを満足する解が正解である。それを求めることが困難であるとき，これらの式のある部分を緩めるような妥協をした結果の解が近似解である。多くの近似解法は次の2つのカテゴリーに分類される。

・静力学的可容解

平衡方程式と力の境界条件だけを満たす応力を求め，得られた応力から構成式によりひずみを評価する解法であり，ひずみ−変位の関係と変位の境界条件を満足する保証は無い。これは材料内部の物質点が変形後も連続となる保証はないことを意味している。多結晶モデルとしてのSacksモデル⁵⁾がこれに該当する。

・運動学的可容解

変位の境界条件を満足する連続な変位場を求め，ついでひずみの変位による定義式から求まるひずみに対応する応力を構成式により求める解法であり，得られた応力場が平衡条件および力の境界条件を満たす保証は無い。多結晶モデルとしてのTaylorモデル⁶⁾がこれに該当する。Hillは教科書出版後の1950年代から1960年代にかけて弾塑性体，剛塑性体に関する境界値問題の解の唯一性と平衡の安定性に関する系統的な理論を展開しており，その一連の論文で剛塑性体に対する2種類の近似解すなわち静力学的可容応力場と運動学的可容速度場に関連する極値定理(変分原理)⁷⁾を提示している。

・第1変分原理(下界定理)

静力学的可容応力場による物体表面外力をF*として   

$$\int F^{*}\cdot vdS\le \int F\cdot vdS$$ (3)

・第2変分原理(上界定理)

運動学的可容速度場をv*とすると   

$$\int F\cdot v^{*}dS\ge \int F\cdot vdS$$ (4)

この変分原理は塑性仕事速度について静力学的可容応力場はその下限値を与え，運動学的可容速度場はその上限値を与えることを示している。そこで，塑性仕事速度が極大となるような静力学的可容応力場，あるいは極小となるような運動学的可容速度場を求めることで精度の高い近似解を得ることが出来る。そのためには応力場ないし速度場の数式表現が必要となるが，静力学的可容応力場は偏微分方程式である平衡方程式の一般解であり，その表現を得るのは困難であるが，後者は多項式などの連続関数の組み合わせで対応出来る。そのため多くの近似解法は上界定理に基づいている。圧延，押し出し，引き抜きなどの定常解に対する近似解法として展開された上界接近法の離散化解法(UBET)⁸⁾もその範疇の数値解法である。古来より変分法の直接法による近似解では解析対象の全領域で定義される関数表現を基本としていたが，上界接近法では区分的な領域のみで定義される解の集合として近似解を構成する方法もとっており，後で述べる変位法による有限要素法のコンセプトに近い。

2・3　成形解析における有限要素法の現状と課題

Zienkiewick⁹⁾は有限要素法を「無限の自由度もつ連続体を，有限個の未知量を含んだ部分領域としての要素の集合体で近似する解法」と定義している。換言すれば連続体を離散化集合体で近似する手法であるが，その際

1.隣接する要素間の変位の連続と力の釣り合いは有限個の節点においてのみ行う。

2.要素内の未知量としての応力，ひずみ，および変位は要素毎に独立に定義する。

なる方針でモデリングすることで，いかなる複雑な対象であっても，要素単位の単純な操作の多数回の繰り返しで処理できることになる。有限要素法では基本未知量を変位とする変位法の他に応力を基本未知量とする応力法，変位の他にひずみも未知量とする仮想ひずみ法およびこれらの組み合わせによるハイブリッド法などが提案されているが，最も多く利用されているのは変位法である。変位法では要素毎に変位場を定義し，これよりひずみ変位関係および構成式を通してひずみおよび応力を評価するので運動学的可容解に属する近似解法となる。また，要素内の応力場の釣り合いの微分方程式(2)に代えて   

$$\int \sigma :\delta \epsilon d{V}_e=\int t\cdot \delta vd{S}_e$$ (5)

なる仮想仕事等式を適用することで微分の概念を放棄して，すべての関係式を代数式で表現可能としている。ここに左辺の積分領域は要素であり，右辺の積分はその要素境界面であるから，要素内の応力のなす仮想仕事と等しい仕事をなす節点力tを定義することで応力の釣り合い条件としていることになる。

有限要素法出現以前の主な数値解法は微分方程式の差分解法であった。差分法では支配方程式のタイプと境界条件に応じてスキームを変える必要があり個別対応であったが，有限要素法は全体の解を局所的な関数の集合で表現することで複雑な境界条件を普遍的な形式で扱えるなどの汎用性を確保できた。そのため，工学分野での数値解法として急速に普及し，現在では成形解析分野でも数値シミュレーションの代名詞となっている。

一方，有限要素法は，Turnerら¹⁰⁾が航空機の骨組み構造解析手法として開発したのがそのルーツであり，同じ固体力学分野であっても成形解析とは目的が異なり，有限要素法に期待する出力情報も異なる筈である。たとえば破壊の評価に対して，構造設計では破壊に耐えうる荷重レベルの予測が主な命題であり変形は特別な場合を除いて殆ど関心がない。一方，成形解析では破壊に至る迄の変形量を的確に予測することが求められる。変形状態の予測精度を向上させるためには変形後の形状に基づいた釣り合条件式を解かねばならないが。変形後の形状は未知量であるから，いわゆる幾何学的非線形性を有する式を解かねばならない。

幾何学的非線形を有する問題の解を線形近似が許容されるような微少増分解の積み重ねで構成する手法は古くから知られているが，この解法に基づく有限要素法は，Hill¹¹⁾が提示した速度型境界値問題の変分原理を用いてMacMeeking and Rice¹²⁾が定式化したのが最初であり，それは次の速度型仮想仕事の原理を基礎としている。   

$$\int {\dot{\Pi }}^T:\left(\delta v\otimes \nabla \right)dV=\int \dot{t}\cdot \delta vd{S}_F$$ (6)

ここに， $\dot{\Pi }$ はいわゆる公称応力速度テンソルであり，Cauchy応力速度テンソルとは次の関係にある。   

$$\dot{\Pi }=\stackrel{\nabla }{\sigma }+\omega \cdot \sigma -\sigma \cdot \omega +\text{tr}\left(L\right)\sigma -L\cdot \sigma$$ (7)

右辺第1項の $\stackrel{\nabla }{\sigma }$ はCauchy応力のJaumann速度と呼ばれ，客観性のある応力速度テンソルの一つである。ここで，基準配置を変形前(t₀=0)の配置に固定するtotal Lagrange形式と現配置(t₀=t)に更新しながら求解を進めるupdated Lagrange形式とがあり，構成式が変形前を基準に表現される超弾性では前者が，現配置で構成式を記述する弾塑性では後者の形式が都合が良い。これは弾性変形はその大小に拘わらず無応力状態では無変形状態に復帰するという可逆性を有するが，塑性変形は非可逆であり変形前の配置は変形後に参照不可能となるからである。そのため，この塑性変形特有の材料非線形(非可逆性)の速度型構成式を積分する際には付帯条件としての降伏条件を満たすようになされなければならない。これを解決する目的で提案された積分法のなかで代表的なものはKrieg and Krieg¹³⁾によるRadial return法であろう。

さて，有限要素法では全体剛性方程式と呼ばれる総節点数を次元とするベクトル方程式を解かねばならなない。成形解析ではそれに要する長大な計算時間と要求メモリー数の多さが実モデルの解析の障害となている。この問題を解決すべく登場したのが動的陽解法である。これは本来クラッシュなどの動的な問題解析のための有限要素法であり，解くべき方程式は次式で表される。   

$$M\ddot{u}+C\dot{u}+Ku=F$$ (8)

第1項の慣性項，第2項の減衰項は動的問題に固有な項であり静的釣り合いの下での変形を求める通常の成形解析では存在しない。動的陽解法による成形解析では，あえてこの項を加えた上で要素の質量を節点のみに配分する集中モデルを採用して質量マトリックスMを対角マトリックス化し，その他の項を現在配置で評価することで増分形式をスカラー方程式化して大幅な計算時間の短縮と要求メモリー量の削減を図った手法であり，主として板材成形分野で広く利用されている。ただし，現実的な解析条件を得るために恣意的な質量や減衰係数を設定することや，応力の振動などの影響で応力値の評価精度は静的解法には及ばない。

固体の変形解析では物質点に着目して変形履歴を追跡するので要素を材料に密着させているので変形に伴って要素形状が初期形状から変化し，やがては要素体積が0となる要素のつぶれに至ることが少なくない。この問題の対処法として，要素のゆがみが著しい領域で要素の再分割を行う方法がよく行われているが，この方法ではそれまでの要素積分点情報を再分割後新たに定義された要素積分点に再配分するための原理，つまり新旧要素間での遺産相続のための物理法則がないのが欠点である。これに対してALE(Arbitrary Lagrangean and Eulerlian)法¹⁴⁾は合理的な手法であるが，自由表面が多く材料の空間移動の大きい板材成形問題への適用が難しそうである。

20世紀後半には成形解析手法は目覚ましい発展を遂げきた。まだ解決しなければ課題が多数残されているにせよ，大部分の実成形問題の解析が可能で多くの企業の設計生産現場で有効に利用されている現実を見ると正に隔世の感を禁じ得ない。現在のところ成形解析の手法の中心は有限要素法であるが，それによる成形解析に期待されるものの多くは成形過程に発生する不良現象，例えば板材成形ではわれ，しわ，スプリングバックなどの事前予測とその対策の指針となるを情報を提供することである。これらの不良現象の要因の殆どは成形過程で材料内部に発生する応力(ひずみではない)であるが，内部応力の評価は特別な場合を除き実測不可能であり，それが可能な有限要素解析は他に代え難い貴重なツールである。しかしながら，有限要素法はこれらの不良現象を直接予測・評価し，かつその対策の指針を与える機能はないので，有限要素解析がもたらす応力データに基づいて，不良現象を解析する理論も必要である。これらについても多くの提案がなされているが紙数の都合で具体的な紹介は割愛する。

2・4　成形解析手法としての弾塑性有限要素法

2・4・1　成形解析で考慮すべき非線形性

本節では，成形解析手法としての有限要素法の考え方を概観する。金属材料の成形加工プロセスを有限要素法により解析するためには，主に幾何学非線形，材料非線形，そして接触非線形を精度良くかつ合理的に取り扱うことが不可欠である。材料非線形については第3章にて詳述するため，本節では幾何学非線形と接触非線形について概説する。なお以降では解析手法として弾塑性有限要素法を用いることを前提とする。

(1)有限変形に起因する非線形(幾何学非線形)

変形やひずみが十分小さいと仮定した理論を微小変形理論，またこの仮定を排除した理論を有限変形理論などと呼ぶ。加工成形では大きな変形あるいはひずみを伴うため，その解析には有限変形理論に基づく定式化を行うことが不可欠である。有限変形理論では通常，次式の変形勾配テンソルFとその(右)極分解が冒頭で定義される。   

$$F=\frac{\partial x}{\partial X},F=RU$$ (9)

ただしXおよびxはそれぞれ変形前後の配置に関する位置ベクトル，またRおよびUは直交テンソルおよび正値対称テンソルである。第一式は変形前後で配置を厳密に区別することを，また第二式は物体の変形はRによる剛体回転とUによる伸縮に分解できることを示す。このように配置と剛体回転の考え方が厳密に定義されているか否かが，微小変形理論と有限変形理論の大きな違いである。

例として変形の評価に及ぼす剛体回転の影響を考える。Fig.1のように，時刻tにおいてx方向に引張荷重を受けた棒が，微小時間増分Δt後にはその状態を保ったまま90°回転したとする。空間固定されたxy座標系を参照して変形前後の応力を成分表示すると，以下のように回転の影響で成分は大きく変化する。   

$${\left[\begin{array}{cc}{\sigma }_{xx}& {\sigma }_{xy}\\ {\sigma }_{yx}& {\sigma }_{yy}\end{array}\right]}^t=\left[\begin{array}{cc}\sigma & 0\\ 0& 0\end{array}\right],{\left[\begin{array}{cc}{\sigma }_{xx}& {\sigma }_{xy}\\ {\sigma }_{yx}& {\sigma }_{yy}\end{array}\right]}^{t+\Delta t}=\left[\begin{array}{cc}0& 0\\ 0& \sigma \end{array}\right]$$ (10)

Fig. 1.

 A rod subjected to rigid body rotation.

一方棒に作用する引張荷重は回転前後で変化しないため，ひずみ増分はゼロである。したがって材料構成式から計算される応力増分もゼロになるはずであり，式(10)に示す結果とは大きく矛盾する。この二つの考え方で見られる不整合は，応力を観測する基準枠の違いに起因する。例えばFig.1のx’y’座標系のように棒とともに回転する共回転座標系を参照すると，次式のように応力成分は回転前後で変わらない。   

$${\left[\begin{array}{cc}{\sigma }_{x^{\prime }{x}^{\prime }}& {\sigma }_{x^{\prime }{y}^{\prime }}\\ {\sigma }_{y^{\prime }{x}^{\prime }}& {\sigma }_{y^{\prime }{y}^{\prime }}\end{array}\right]}^t={\left[\begin{array}{cc}{\sigma }_{x^{\prime }{x}^{\prime }}& {\sigma }_{x^{\prime }{y}^{\prime }}\\ {\sigma }_{y^{\prime }{x}^{\prime }}& {\sigma }_{y^{\prime }{y}^{\prime }}\end{array}\right]}^{t+\Delta t}=\left[\begin{array}{cc}\sigma & 0\\ 0& 0\end{array}\right]$$ (11)

この結果は後者の考察と整合する。このように，共回転する基準枠から観測した応力を共回転応力と呼ぶ。共回転応力速度 $\stackrel{\nabla }{\sigma }$ の具体形として，次式が用いられることが多い。   

$$\stackrel{\nabla }{\sigma }=\dot{\sigma }-W\cdot \sigma +\sigma \cdot W$$ (12)

ただし $\dot{\sigma }$ は空間固定された基準枠から観測した応力速度，Wは(連続体)スピンテンソルである。式(12)をJaumann応力速度と呼ぶ。共回転応力速度 $\stackrel{\nabla }{\sigma }$ は客観性を満たすことから，有限変形理論における弾塑性構成式では応力速度として $\dot{\sigma }$ の代わりに $\stackrel{\nabla }{\sigma }$ を用いることで構成式の客観性を保証している。

なお式(12)は，物質点はスピンテンソルWとともに回転することを意味している。1980年代頃までは物体の瞬間的な回転(スピン)を表すテンソルとして何を用いるべきか活発に議論され^15,16)，様々な客観性のある応力速度が提案された。現在では実用上大きな問題がないとしてJaumann応力速度が広く用いられている。

(2)接触非線形

有限要素法では，成形中に生じる素材と工具の接触は，素材側の節点(材料節点)が工具面を構成する要素(工具要素)に幾何学的に接触したかどうかを判断基準とする場合が多い^16,17)。幾何学的に工具要素に接触したと判定された材料節点は，それ以後工具内部へ潜り込まないように変位の工具面法線方向成分が拘束される。この変位拘束によって生じる節点力は工具反力に相当し，工具反力が生じている間は接触状態が持続される。その後工具反力がゼロになると材料節点は工具から離脱したと判定され，それ以後変位拘束は与えられない。以上のように，素材と工具の接触状況が変化すると境界条件や摩擦状態が変化するため，結果として強い非線形性がもたらされる。

材料節点と工具要素の幾何学的関係から接触判定を行うアルゴリズムを，接触探索アルゴリズムと呼ぶ。材料節点が工具要素上に厳密に位置することは数値計算上ほぼありえないため，解析上で“接触”をどう定義するかが解析精度上重要である。広く用いられている手法では，材料節点から工具要素へ垂線を下した時の距離が許容値以内のとき接触と判定される。しかしこの方法では，垂線の定義によってはデッドゾーンなどの問題により適切に接触判定できない場合があり，その対策が必要となる¹⁶⁾。また，それぞれの材料節点に対して数多くある工具要素の中から接触関係にある要素を見つけ出す作業は，計算時間に大きく影響する。そこで通常は接触探索を，接触する可能性がある工具要素を絞り込む粗探索と詳細な計算に基づいて接触判定を行う詳細探索の2段階に分ける場合が多い¹⁶⁾。

接触した材料節点に与える変位拘束には，ラグランジュ未定乗数法やペナルティ法が用いられる場合が多い。ラグランジュ未定乗数法は変位拘束を付帯条件として方程式の極値を求める数学的手法である。この方法では，計算の過程で剛性マトリクス対角項にゼロが入るため計算処理が少し煩雑になるが，境界条件を厳密に満たすことができる。それに対してペナルティ法は，工具内部へ潜り込もうとする材料節点に対して非常に大きな抵抗力を発生させる方法である。抵抗力の大きさは予め設定されたペナルティ数によって決定される。ペナルティ法では，対応する剛性マトリクス対角項にペナルティ数を加えるなどの処理により変位を拘束することができるため，計算上の処理は簡便である。その一方でペナルティ数が解析結果に影響を及ぼすリスクがある。基本的にはペナルティ数が大きいほど高精度な変位拘束が実現できるが，Fig.2のように工具面法線方向が空間固定された座標軸に対して斜めを向く部位に材料節点が接触した場合は非対角項にもペナルティ数が入るため，ペナルティ数が過度に大きいと数値計算上の丸め誤差により精度が低下する場合がある。ただし，接触した節点の成分を工具面法線方向を一つの軸とした局所座標系(Fig.2におけるx’y’座標)で表示すれば，ペナルティ数は必ず対角項にしか入らないためこの問題を回避することができる。

Fig. 2.

 A work-piece node sliding on a curved part of tool surface. □ and ● denote the tool and work-piece nodes, respectively.

有限変形問題において剛体回転の取り扱いが極めて重要なことは前述のとおりであるが，同様の問題が工具反力でも生じる。例えば工具の曲率部では，部位によって工具面法線ベクトルの向きが異なる。そのため接触した材料節点が工具の曲率部上を移動する場合，Fig.2のように移動に伴って変位拘束する方向を変化させる必要がある。またさらに工具反力(節点力)ベクトルも工具面法線ベクトルと共回転させる必要がある¹⁸⁾。これは言い換えれば，工具反力ベクトルfの速度を基底ベクトルの変化率 ${\dot{e}}_i$ を考慮した次式で与えることに相当する。   

$$\dot{f}={\left(f_i{e}_i\right)}^{\cdot }={\dot{f}}_i{e}_i+f_i{\dot{e}}_i$$ (13)

このとき， ${\dot{e}}_i$ は工具形状と材料節点の速度から決定される。したがって式(13)第二項は工具の幾何情報に基づいて離散化し，材料節点速度の関数として表す必要がある¹⁹⁾。

2・4・2　弾塑性有限要素法の各種解法

成形解析のための弾塑性有限要素法としては，静的陰解法，動的陽解法，静的陽解法などの手法が良く用いられている。静的/動的とは静解析と動解析のどちらを採用しているかを表し，具体的には力の釣り合い条件式に慣性項を考慮するか否かが大きな違いとなる。動的陽解法の考え方は2・3節で述べられた通りである。一方陰解法/陽解法は数値計算上のアルゴリズムの違いを表す。弾塑性有限要素法では通常プロセスを多くの微小区間(ステップあるいはインクリメントと呼ばれる)に分割して解析を進める増分解析が用いられる。静的陽解法では各ステップにおける非線形性を無視して区分線形近似を行う^16,18)。これは，ステップ内では接線剛性マトリクスおよび境界条件が変わらないと仮定することに相当する。しかしながら2・4・1節で述べた各種非線形性の影響により，実際には区分線形近似が成り立つ保証はない。そこで通常は，r_min法^16,20)により非線形性の発生要因を監視して区分線形近似が成り立つようにステップ幅を順応的に調整する。そのため静的陽解法ではステップ数が膨大になる傾向にあり，計算コスト上大きなデメリットである。また以上の計算上の工夫を行ったとしても，内力と外力の釣り合いが満足される保証はない。そこで静的陽解法では，速度型仮想仕事の原理式を次式のように修正して用いる場合が多い。   

$$\int \Delta {\Pi }^T:\left(\delta v\otimes \nabla \right)dV=\int \Delta t\cdot \delta vd{S}_F+\int t\cdot \delta vd{S}_F-\int {\sigma }^T:\left(\delta v\otimes \nabla \right)dV$$ (14)

右辺第2項と第3項の差が不釣り合い力であり，発生した不釣り合い力を次ステップ以降で補正する役割を担う^16,19)。またこの他に，これらの項を加える代わりに陽的に不釣り合い力補正を行うステップを随時導入する手法(ALGONEQ)²¹⁾などが提案されている。

一方静的陰解法では，各ステップにおける非線形性を認めてニュートン・ラフソン法等により力の釣り合いを満たしながら解析を進める。このとき，幾何学非線形および材料非線形(応力の時間積分)によって生じる不釣り合い力は同一の反復計算中で処理される。一方反復計算中は境界条件を更新しないため，反復計算収束後に接触状態を更新すると新たに不釣り合い力が発生する可能性がある。そこで厳密に釣り合いを満たすためには，接触状態の更新により生じた不釣り合い力を改めて反復計算によって補正し，最終的に不釣り合いが完全に補正できるまでこの作業を繰り返す必要がある。そのため，複雑な成形プロセスではしばしば計算が収束せず解が得られない場合が生じる。以上のような技術的課題から，実際には接触状態の更新に伴う不釣り合い力を無視する準陰解法が採用される場合も多いようである。この問題に対して，接触による変位拘束にaugmented Lagrangian法を用いることで反復計算の多重構造を回避する方法も提案されている²²⁾。

以上のように，解法によって前提とする考え方や計算の進め方に大きな違いがある。そのため解法が解析結果に及ぼす影響についても様々検討されている^23,24)。一方で，静的陽解法における解析精度の問題と静的陰解法における計算の収束性の問題から，両者は徐々に接近しているように感じられる。静的陽解法/陰解法という分類自体の意味がなくなる日も近いかもしれない。

2・5　今後解決すべき課題と最近の展開

計算機と解析技術の発達に伴い，最近では有限要素法によりあらゆる加工法の解析が可能になりつつある。しかしながら，“解析できる”ことと“高精度に予測する”ことは必ずしも等価ではなく，解析の要素技術としては未だ発展途上の課題も多い。本節では加工成形解析の現状と今後解決すべき課題を展望する。

(1)有限要素に関する問題。

標準的な要素としてはソリッド要素(連続体要素)が良く用いられるが，板材成形ではシェル要素(構造要素)が用いられることが多い。またこれらの要素において次数や体積積分法を改良した多くの亜種が提案されている。一方で，要素タイプによって異なる解析結果が得られる場合が多いこともよく知られている。一例として，初期張力を与えたハット曲げ成形におけるスプリングバック解析をシェル要素およびソリッド要素を用いて行ったときの，初期張力と縦壁部曲率の関係をFig.3に示す²⁵⁾。要素タイプ以外の解析条件は全く同一であり，ダイとブランクホルダの間隔は定隙間としている。要素タイプによって解析結果が大きく異なることが明らかである。しかしながらこの原因は十分解明されておらず，解析対象に合わせた要素の使い分けや要素に合わせた解析パラメータの微調整などが不可欠なのが現状である。

Fig. 3.

 Relationship between initial tensile stress and sidewall curvature²⁵⁾.

またバルク加工のように数十%のひずみが生じる問題では，変形中の過度の要素ゆがみにより解析精度の低下や解析の停止を招くおそれがある。そのため適宜リメッシュ(メッシュの再分割)やそれに伴うリマップ(物理量の再分配)を行うことが不可欠である。最近の塑性加工解析用ソフトウェアでもこれらは標準的な機能として搭載されている。しかしながら2・3節でも述べられているように，物理的整合性を保ったままリマップを行うことは困難である。また六面体要素を対象としたリメッシング技術も発展途上であり，今後さらなる研究開発が求められている。

これらの問題に対して，ラグランジュ記述とオイラー記述の特徴を兼ね備えたALE法を押し出し加工²⁶⁾やリングローリング²⁷⁾などの大ひずみを伴う成形加工問題へ適用する試みが数多く報告されている。またオイラー記述による成形加工解析に関する研究も進められており，最近でもロール成形プロセス²⁸⁾などへの適用事例が報告されている。一方，2・3節でも述べられているようにオイラー記述では自由表面の取り扱いに課題があり，加工成形解析上で大きな欠点となっている。そこで最近注目されているのがメッシュフリー解析法の加工成形プロセスへの適用である。メッシュフリー法は文字通りメッシュを用いずに連続体の解析を行う手法であり，有限要素法を拡張した手法から連続体を粒子の集合としてモデル化する手法まで様々提案されている。メッシュフリー法ではメッシュに起因する種々の問題に悩まされる心配がないため，最近ではバルク加工を対象に解析事例^{例えば29,30,31)}が増えている。実用的な利用には解決すべき課題も多いが，今後の発展が期待される。

(2)素材と工具の接触に関する問題。

加工成形の有限要素法解析では，工具についてもその表面を三角形や四角形などの面要素(パッチ)により離散化する場合が多い。しかしながらこの方法は工具面形状を区分線形近似することとほぼ同義であり，近似による形状誤差の発生が不可避である。この形状誤差は工具に接触した材料節点の移動経路に影響を及ぼし，解析精度を低下させる原因になりうる¹⁹⁾。したがって解析では高精度な工具モデルを用いることが重要である。

最も形状表現精度が高いのは，CADデータを直接解析で利用する方法であろう^16,22)。しかしながらCADデータの互換性や複雑さなどの問題からあまり普及していない。一方CAD表現でも用いられる非有理関数を用いて，一度離散化された表面を再び平滑化する方法も提案されている³²⁾。この手法では，表面は平滑化されるが形状精度は全く保証されない問題がある。それに対して，Nagataパッチ³³⁾を用いた表現手法が提案されている¹⁹⁾。Nagataパッチは一度離散化された表面を離散化前の元形状へ高精度に復元できるのが特長である。この手法を用いることで表面形状精度を大幅に向上させ，その結果工具モデルに影響されにくい接触解析が遂行できることが報告されている。

続いて被加工材の要素タイプが接触解析へ及ぼす影響を考える。板材成形解析で用いられるシェル要素は，板厚中央面にしか節点を持たず板厚は仮想的に与えられるのみである。これは板材成形では平面応力状態が成り立つという仮定に基づくモデル化だが，このモデル化は接触解析上次のような欠点をもたらす。まず，板厚方向には一つの節点しか持たないため板両面からの接触を表現できない。そのため，フランジ部におけるダイとブランクホルダによる両面接触やパンチ底部におけるパンチとカウンターパンチによる両面接触は原理的に扱えない。また，シェル要素は応力の板厚方向成分を無視しているため，ブランクホルダ力や決め押しを与えてもその効果は原理的に再現できない。液圧成形において板両面から工具反力と液圧が同時に作用する場合にも同様の問題が生じる³⁴⁾。そこで最近では，シェル要素による板材成形解析精度を向上させるため，シェル要素とソリッド要素の特徴を兼ね備えたソリッド−シェル要素^例えば35)の研究開発が盛んである。適用事例も少なく未だ研究段階だが，今後の発展に期待したい。

(3)最適化技術。

最近では，加工成形解析と最適化技術を組み合わせることで加工条件の最適化を行う取り組みが盛んである。初期ブランク形状の最適化^例えば36)や各種加工プロセスの最適化^例えば37)など多くの適用事例が報告されている。近年のこの流れは，産業界における有限要素法解析に対するニーズ拡大とともに解析結果は実用に十分資するとの認識が定着した証であろう。今後も，解析精度の向上を目指した基礎技術の開発とともに最適化をはじめとする応用技術の開発がますます加速すると考えられる。

3. 金属成形解析のための塑性構成式の発展

3・1　塑性構成式の枠組みと概括的展望

20世紀後半において著しい進展を遂げた成形解析手法の理論と技術に較べて，塑性構成式論はそれほど大きな進展はなく，Hillの著書¹⁾で確立された基本の枠組みを超えるような新たな提案はそう多くはない。

塑性構成式の枠組みは

1.降伏条件:f(σ)=σ_Y

2.加工硬化則:σ= σ=F(∫σ:dε^P)，σ，dε^P:相当応力と相当塑性ひずみ増分

3.流れ則:dε^P=C^P(σ，κ)dσ

の3個の基礎的関係で構成されている。

1.降伏条件と2.加工硬化則は応力とひずみのテンソル成分で表わされているがスカラー式であり，最後の流れ則のみが応力増分テンソルと塑性ひずみ増分テンソルの間のテンソル的な関係式である。一般的に応力テンソルは架空の量であり，いかなる実験でも直接測定することはできない^*1。材料変形に関して実験で直接認知できるのは，力と変位であるから材料試験では測定された力と変位から直接応力とひずみのテンソルを評価できなければならない。それにはごく限られた手段しかなく，そこで得られた情報は一般の変形状態に対して，限られた状況での情報しか得られない。そこで，これを一般化するためには物理的な背景に基づいた仮説が必要となる。

*1　応力テンソルはそれ自身が直接内力を表すものではなく，Cauchyの公式t=σ^T・nを通して単位法線ベクトルnなる物体内の単位面要素に作用する内力ベクトルtを評価する。

この観点からは，最初の二つはまさしく，単軸引張り試験で得られた応力−ひずみ曲線の一般化であるが，唯一のテンソル的な関係式である「流れ則」の具体的な表現を得るための情報は単軸引張り試験からは得られないので，新たな仮説が必要となる。

Hillが提示した塑性ポテンシャル論では，次の法線則(Normality rule)   

$$d{\epsilon }^P=\frac{\partial g\left(\sigma \right)}{\partial \sigma }d\lambda$$ (15)

がそのための仮説であり，塑性ひずみ増分は塑性ポテンシャル関数g(σ)の現在の応力点の外向き法線方向にのみ生ずることを主張している。法線則はさらに降伏関数f(σ)を塑性ポテンシャルとする関連流れ則(Associated flow rule)と降伏関数とは異なる関数を塑性ポテンシャルとする非関連流れ則(Non- associated flow rule)に分類されるが，成形解析に用いられる塑性構成式の殆どが関連流れ則の枠組みを踏襲している。従って20世紀後半の塑性構成式についての展望は，上述の法線則の特徴に由来する流れ則の欠陥を解消する提案の他は関連流れ則の枠組みを構成する基本要素である降伏関数，加工硬化則の表現の改良にとどまるものが殆どである。

3・2　法線則による流れ則の問題点とその解決策

法線則では塑性ひずみ増分dε^Pの方向は応力増分dσの方向とは無関係に塑性ポテンシャル面の現在の応力点での法線方向に唯一に定まる。つまりdσの方向変化にdε^Pの方向が追従することを完全に否定するので，しわや局所くびれなどのひずみ増分方向の急激な変化を伴う分岐現象に適用できない。Hillの局所くびれ理論が張り出し領域に適用できないことも法線則のこの欠陥に由来している。そもそも，この問題は1940年代に弾塑性平板の座屈解析でJ₂流れ理論^*2による解よりも，塑性力学的には不合理な関係とされるJ₂変形論の増分形式による解の方が現実的な値を与えるというパラドックスが露呈したことが発端となっている。これはJ₂流れ理論では塑性ひずみ増分方向が応力増分方向に依存しないことでdε^Pとdσの1対1対応が喪失し，流れ則dε^P=C^Pdσの係数テンソルはdetC^P=0なる特異性を示すことが原因であることが明かにされた³⁸⁾。

*2　降伏関数を偏差応力テンソルの2次の不変量J₂≡σ':σ'りする関連流れ則，歴史的には塑性ひずみ増分が偏差応力に比例するとしたLevy-Misesの流れ則と同一

J₂変形論の増分形式とはHenckyの変形論   

$${\epsilon }^P=\frac{3{\overline{\epsilon }}^P}{2\overline{\sigma }}{\sigma }^{\prime }$$ (16)

を増分形にした次式である   

$$d{\epsilon }^P=\frac{3}{2}\left\{\frac{1}{E_S}d{\sigma }^{\prime }+\left(\frac{1}{H^{\prime }}-\frac{1}{E_S}\right)\frac{d\overline{\sigma }}{\overline{\sigma }}{\sigma }^{\prime }\right\}$$ (17)

ここにE_S，H'は加工硬化曲線のセカント係数および接線係数であり，右辺第1項が塑性ひずみ増分方向への応力増分方向依存性を表している。後にBudianskyによって降伏曲面上に角点の存在を前提に，ある限られた負荷方向に対してはHenckyの変形論の正当性が評価された³⁹⁾。いわゆるコーナーモデルの草分けで，板材のくびれ問題への適用例⁴⁰⁾もある。その後Christoffersen and Hutchinson⁴¹⁾，Gotoh⁴²⁾などの提案があり，座屈や局所くびれ問題に適用され概ね良好な結果を与えている。降伏曲面の角点は曲面の特異点であるため，その法線方向が一意に定まらないことを利用して塑性ひずみ増分方向の応力増分依存性を導入したもので，法線則の枠組みにおける解決策である。ただ，角点の頂角の実験的な同定が困難なことや，降伏曲面の至る所に角点が存在する曲面の表現が不可能等の問題がある。そこで，Itoら⁴³⁾は，法線則を放棄して次に示すような塑性ひずみ増分の応力増分依性構成式を提案し3次元局所分岐解析に適用していくつかの有用な成果を得ている。   

$$d{\epsilon }^P=\Lambda \left[K_{C}d{\sigma }^{\prime }+\left(1-K_C\right)\left(n_F:d{\sigma }^{\prime }\right)n_N\right]$$ (18)

  

$$n_F=\frac{\partial \overline{\sigma }}{\partial \sigma }/\left|\frac{\partial \overline{\sigma }}{\partial \sigma }\right|,n_N=\frac{\partial {\overline{\sigma }}^{*}}{\partial \sigma }/\left|\frac{\partial {\overline{\sigma }}^{*}}{\partial \sigma }\right|$$ (19)

ここに，0<K_C<1は応力増分依存性の程度を示す材料パラメータであり，K_C=0ではn_Nを塑性ポテンシャル曲面の単位法線ベクトルとする非関連流れ則に形式が一致する。

3・3　異方性の表現に関する展望

初期異方性にはR値で通称される変形異方性と降伏応力異方性の2種類がある。これらは互いに連成することはあっても基本的には互いに異なる異方性である。しかるに関連流れ則では降伏関数はひずみ比を決定するための塑性ポテンシャルの役割も兼務させられるので異方性構成式への拡張の際に過大な任務を負わされる。

異方性降伏関数は等方性降伏関数の異方性への拡張であるから，すべての異方性降伏関数は基準となる等方性降伏関数に帰着可能でなければならない。そこで読者には周知のことではあるが最初に基準となる等方性降伏関数を提示しておく。

特別な場合を除いて金属の塑性変形は体積一定で静水圧不依存であることおよび材料要素に対する座標変換に対して不変な形式であるべき条件を満たす最も単純な形式がMisesの降伏条件として知られている次式で等方性降伏関数のスタンダードである。   

$$f\left(\sigma \right)=\overline{\sigma }=\sqrt{\frac{3}{2}{{\sigma }^{\prime }}_{ij}{{\sigma }^{\prime }}_{ij}}=\sqrt{\frac{3}{2}{\sigma }^{\prime }:I:{\sigma }^{\prime }}$$ (20)

Hill⁴⁴⁾は式(20)によるMIsesの降伏関数の表現形式(偏差応力テンソルの2次形式)を踏襲して恒等テンソルIを異方性テンソルAに置き換えた   

$$f\left(\sigma \right)\sqrt{\frac{3}{2}{\sigma }^{\prime }:A:{\sigma }^{\prime }}$$ (21)

を初期異方性降伏関数の表現として提示をした。これはHill48と通称される異方性降伏関数のはしりで，平面応力状態では次のように表される。   

$$f=\left(G+H\right){\sigma }_x^2+\left(H+F\right){\sigma }_y^2-2H{\sigma }_x{\sigma }_y+2N{\tau }_{xy}^2=0$$ (22)

ここで，F，G，H，Nは異方性テンソルAの非ゼロ成分である。関連流れ則を前提として，3つの独立な方向の単軸引張試験によるR値の測定から，これらの異方性パラメータ決定する方法が一般的であるが，このように決定された異方性パラメータを用いた降伏関数から評価される降伏応力の方位分布が現実に適合しない場合が少なくない。本来降伏関数は降伏応力の一般化であるから逆に降伏応力の方位分布から異方性パラメータを決めるべきであるが，その場合はR値の方位分布に適合しないことは自明である。これは関連流れ則における異方性降伏関数としての2次形式の表現能力の限界を越えたと考えるべきである。

また，板材成形の分野では面内は等方とし，板厚異方性のみを考慮した取り扱いも良く行われる。その場合平均的なR値(R)のみが異方性パラメータとなるが，Hillの異方性降伏関数ではX値(等2軸引張り降伏応力/単軸引張り降伏応力)   

$$X=\sqrt{1+\overline{R}/2}$$ (23)

はR値とともに増加すことになる。しかしながら，アルミ合金などはR<1でX>1であることが知られたのを契機として，異方性降伏関数の改善の提案が種々なされてきた。大別して非多項式を採用したものと高次の多項式による表現である。非多項式型の異方性降伏関数の表現は異方性主軸と応力主軸が一致する場合についてHill⁴⁵⁾，Logan and Hosford⁴⁶⁾の提案があり，これをBarlat and Lian⁴⁷⁾が異方性主軸と応力主軸が一致しない場合へ次のように拡張した。   

$$f=a{\left|K_1-K_2\right|}^M+a{\left|K_1+K_2\right|}^M+\left(2-a\right){\left|2k_2\right|}^M$$ (24)

  

$$K_1=\frac{{\sigma }_x+h{\sigma }_y}{2},k_2=\sqrt{{\left(\frac{{\sigma }_x-h{\sigma }_y}{2}\right)}^2+p^2{\tau }_{xy}^2}$$ (25)

上式はBarlatモデルのプロトタイプで異方性パラメーターはa，h，pの3個であるが指数Mが材料パラメータに加わる。むろんM=2と選べば2次降伏関数に一致する。指数Mの値が大きくなるほど，降伏曲面の形状はトレスカ型に近づく傾向を持っており，Mを適当に選ぶことでアルミ板のようなFCC金属の降伏曲面に適する。Barlatらは結晶方位分布関数に基づく重みをつけたTaylor modelにより指数Mを決定している。しかしながら，Hill48の問題点は，そのまま上式のBarlat Modelにもあてはまる。すなわち，応力の異方性と変形の異方性を両立させられるかどうかということである。ただ，Hill48よりパラメータが一つ多い分適合精度が上がる可能性はあるが，それでも不十分ということで，さらにパラメーターの数を増やしたモデルに年毎にバージョンアップし，2003年モデル(Yld2000)はパラメーター数8個の極めて複雑な式となっている⁴⁸⁾。一方，多項式の形式で，パラメータを増やす目的での改善策としてGotoh²⁶⁾は4次降伏関数を提案している⁴⁹⁾。これも異方性パラメーターは8個である。いずれもHill48 の3個に較べてはるかに多くのパラメータを導入しているので，それだけ精度は高いのは当然であるが，その分要求される材料試験の負荷が高いことと複雑な数式故に，それぞれの相当応力に対して塑性仕事共役な相当塑性増分の表現を見出すのが困難であるのが欠点である。その根本の原因は関連流れ則の枠組みが異方性表現の役割を降伏関数のみに押し付けていることに他ならない。異方性降伏関数の表現はいかに複雑で技巧的な数式表現であっても，応力ひずみ曲線におけるSwift 則同様，特別な物理的な意味は無く所詮は曲線のあてはめの近似式に過ぎない。この問題の解決には関連流れ則の枠組み自体の見直しから始めるのが王道と思われる。

さて，これまで展望して異方性降伏関数は基本的に初期異方性の表現に関するものである。これとは別に硬化発展則の異方性表現の問題がある。硬化発展則とは塑性変形の進行に伴う降伏応力(塑性流動応力)の変化を記述した式である，一方塑性負荷が継続している場合の塑性流動応力は降伏曲面上にあるので，硬化発展則は降伏曲面の塑性変形に伴う形状，大きさおよび原点に位置の変化を記述するとも考えられる。等方硬化(Isotropic hardening)とはこの硬化発展則が負荷方向および負荷経路履歴に依存しない硬化特性であり，降伏曲面は原点の位置と形状を変えずにそのサイズだけが塑性変形とともに膨張するモデルで表現されている。これに対して，異方硬化モデルは降伏曲面の位置と形状も変化させる硬化モデルであるが，後者はあまりにも複雑なモデルになるため，少なくとも成形解析には殆ど登場していない。前者のモデルはPragar⁵⁰⁾をによる移動硬化モデルおよびその発展形^51,52)の提案を経て，主として曲げ曲げ戻しを受ける板成形のスプリングバック解析に適用されているものもある⁵³⁾。これらはバウシンガー効果の表現を対象としているが，これと同種の異方性硬化現象である交差硬化を対象としたモデルは殆どみられない。

3・4　最近の研究動向

最近の加工成形解析に関する研究動向を調べると，材料モデルに関連した研究開発が例年大きな割合を占めている。特に板材成形分野では解析における材料モデルの重要性が認識され，弾塑性構成式の精緻化を目指して精力的な研究が展開されている。弾塑性構成式の高精度化を達成するには，バウシンガー効果と塑性異方性を精度良く再現することが重要な鍵となる。塑性異方性を高精度に表現する手段として，異方性降伏関数の研究が活発に展開されている。3・3節で詳述されているように，Hill^44,45)やGotoh⁴⁹⁾をはじめとしてこれまで多くのモデルが提案されてきた。成形解析における活用度と解析精度の観点からみると，Barlatらによる一連の降伏関数はその中でも特筆すべき成果であろう。このモデルは，当初はアルミニウム合金板を対象として発展したが⁴⁷⁾，最近のモデル(Yld2000-2d)⁴⁸⁾ではパラメータを適切に選択することで種々の金属材料の特性を再現できることが報告されている^67,68)。高い解析精度に加えてこのような材料を選ばない高い汎用性が，成形解析において幅広く利用される理由であろう。またその他，Yoshidaら⁶⁹⁾による各種鋼板の異方性表現を目的とした6次の降伏関数や，Cazacuら⁷⁰⁾によるマグネシウムやチタンのように強い引張−圧縮非対称性を示す材料に適用可能なモデルなどが提案され，その検証と加工成形解析への適用が進められている。

以上のモデルでは，初期異方性が変形中も持続されること，また等方硬化であることが仮定されている。しかしながら実際には，3・3節でも述べられているように多くの板材は変形に伴って塑性異方性が変化し異方硬化挙動を示す。そこで最近では，塑性異方性の発展を考慮したモデルの検討が進められている。例えばYanagaら⁷¹⁾は，Yld2000-2dでは諸パラメータを変化させることで降伏曲面の形状を変化させることができる特徴を利用して，Yld2000-2dにおける諸パラメータを塑性ひずみの関数として表すことで異方硬化を考慮した解析技術を提案している。

バウシンガー効果は，移動硬化あるいは等方硬化と移動硬化を組み合わせた複合硬化により表現する方法が一般的である。最近の特筆すべき成果として，Yoshida and Uemori⁷²⁾によるモデルが挙げられる。このモデルでは，反転負荷後の応力レベルが単調負荷時よりも低下する傾向や急速な加工硬化係数の減少が適切に再現でき，特に鋼板の繰り返し塑性変形を精度良く記述できることが報告されている。最近の高張力鋼板の加工成形解析では最も良く用いられるモデルの一つである。また，Teodosiu and Hu⁷³⁾は材料の微視的な構造変化を内部変数として巧みに取り込んだモデルを提案している。本モデルは材料物理的に明瞭という特徴がある一方で，膨大な数のパラメータを同定する必要がある点がネックとなり実用化は進んでいないようである。

一方最近Barlatらは，移動硬化を用いずに降伏曲面の変形によってバウシンガー効果や異方硬化，交差硬化を表現するモデルを提案している⁷⁴⁾。本モデルにより，低炭素鋼板や二相鋼板など様々な鋼種において各種負荷経路下での応力−ひずみ関係を精度良く再現できることが報告されている。

なお，このBarlatらの新しいモデルをはじめとした上述の最新の研究成果は必ずしも完成されたモデルではなく，現在でも引き続き研究されている。それゆえ今後成形解析において幅広く適用され，3・3節で取り上げられたような歴史の一部となりえるかはこれからの研究に委ねられよう。しかしながら，いずれもそのポテンシャルは十分に備えていると考えられ，今後の研究に注目したい。

3・5　塑性構成式と有限要素法解析技術

材料モデルから予測される応力挙動を解析で適切に再現するためには，空間的および時間的離散化後も材料モデルの理論的整合性を保持できるように対策を講じる必要がある。時間的離散化に伴う問題として，速度型弾塑性構成式における無限小時間増分の仮定と増分解析における有限時間増分の仮定の不整合が挙げられる。例として応力の時間積分を考える。通常弾塑性解析では，時刻t+Δtにおける応力テンソルσ^t+Δtは時刻tにおける応力テンソルσ^tとそのステップで生じた応力増分テンソルΔσを用いてσ^t+Δt=σ^t+Δσのように計算される。ただしΔtは解析における有限時間増分である。このとき材料非線形に起因して，式(12)から計算される応力速度 $\dot{\sigma }$ に単に時間増分Δtを乗じた結果をΔσとすると，時刻t+Δtにおいて降伏条件式を満足する保証はない。そこで2・3節でも述べられたように，Radial return法やReturn mapping法，r_min法などを用いることで時刻t+Δtにおいても(近似的に)降伏条件式を満足しうる適切な応力増分を求める必要がある。また，Radial return法やReturn mapping法のように反復計算を必要とする手法の場合は，反復計算の収束性を向上させるために弾塑性構成テンソルC^epを応力積分法に整合した形(整合した接線剛性)で用いる必要がある⁷⁵⁾。

さらに，剛体回転の影響も無視できない。式(12)において客観性が満足されるのはあくまでも無限小時間増分の仮定が成り立つ場合である。そこで有限時間増分を仮定した増分解析でも客観性が破綻しないように，r_min法などによりステップ幅を極微小にとどめる方法やHughes and Wingetによる有限時間増分でも客観性が保たれる手法⁷⁶⁾を用いるなど，アルゴリズム上の工夫が不可欠である。

一方空間的離散化については，有限要素と材料モデルの整合性が重要である。例えばシェル要素を用いて板材成形解析を行う場合は，材料モデルも平面応力の形式に変形する必要がある。一方近年提案されている材料モデルの中には予め平面応力を仮定しているものも多く⁴⁸⁾，この場合は逆にソリッド要素への適用が困難である。したがって，材料モデルの選択には精度だけでなく解析モデルとの整合性に十分留意する必要がある。

最後に座標系の取り扱いについて述べておく。等方性材料では構成テンソルはどの座標系から記述しても同一の成分表示に帰着するので，構成式の記述にも空間固定された全体座標系を用いることができる。一方異方性を有する材料の場合，構成テンソルは異方性主軸に基づいて構築する必要がある。このとき，異方性主軸は通常物体の変形に伴って回転するとモデル化されるため，解析では異方性主軸とともに回転する局所座標系(物質座標系)を積分点ごとに定義し，全体座標系との変換則を随時更新していく必要がある。またシェル要素の場合は平面応力状態を保持する必要がある。そこで，物質座標系が異方性主軸と平面応力状態を同時に表現できるようにモデル化するか⁷⁷⁾，あるいは平面応力状態を表す局所座標系(シェル座標系)を別途定義し，全体座標系−シェル座標系−物質座標系の変換則を随時更新する必要がある⁷⁸⁾。なお異方性主軸の回転は2・4節で述べたスピンの選択と密接に関連している。しかしながらこの点については未解明な点が多く，現在でも研究が続けられている⁷⁹⁾。このように，加工成形解析において材料の異方性を再現するためには，材料モデルそのものの高度化だけでなく計算アルゴリズムの精緻化も不可欠であるといえる。

4. 結晶塑性解析

4・1　Taylorモデルから多結晶塑性有限要素法に至る展望

最後に結晶塑性学的な観点での解析に関して簡単に触れておく。成形用の金属材料の多くは多結晶体である。ここではSchmid則に基づく単結晶の塑性変形特性を前提とした多結晶体の塑性変形挙動の解析を結晶塑性解析と呼ぶことにする。いわゆる多結晶モデルであり，先駆的なモデルはいわずと知れたTaylorモデル⁵⁴⁾であろう。Taylorモデルは「ひずみ一定」と「最小辷りに原理」を基本的な仮説としている。いずれも証明なしの直感的な提案である。前者の仮説は多結晶を構成する結晶粒の塑性ひずみを等しくすることで結晶粒間の変位の連続を保証するもので運動学的可容解のひとつに他ならない。また，後者の仮説「最小辷りの原理」もBishop and Hillの最大塑性仕事の原理⁵⁵⁾を通して理論的な根拠が与えられており，近年に至るまで多結晶金属の様々な塑性変形挙動の解析に採用されている。

さて，多結晶モデルは結晶レベルの微視的な塑性変形挙動と成形性に影響を及ぼす巨視的な塑性変形挙動と結晶レベルの微視的挙動を結び付ける役割が期待される。そのためには一般に結晶粒毎に異なる微視的な応力増分dσ，塑性ひずみ増分dε^Pと巨視的な応力速度dS，ひずみ速度dE^Pを関連づけなければならない。いわゆるセルフコンシステント(self-consistent)多結晶モデルにおける   

$$\overline{d\sigma }=dS,{\overline{d\epsilon }}^P=d{E}^P$$ (26)

はその為の付帯条件である。ここにdσ，dε^Pは微視的応力増分，塑性ひずみ増分の平均値である。Eshelbyの変態問題の解⁵⁶⁾に基礎をおくKröner⁵⁷⁾とBudiansky and Wu⁵⁸⁾のモデル(KBWモデル)や熱弾性問題におけるDuhamelのアナロジー法によるLin モデル⁵⁹⁾などがその系統に属する。筆者らもこのモデルにより多結晶体の応力増分依存性の解析を行い^60,61)，得られた知見に基づいて前出の応力増分方向依存性構成式の定式化した⁴³⁾。この多結晶モデルは巨視応力も不均一であることを無視出来ない場合には適用出来ない。この問題のブレークスルーを期待して登場したのが有限要素多結晶モデルである。国内でのルーツはMiyamotoら⁶²⁾，Jimma⁶³⁾らに遡るが，基本的には1要素を1結晶粒とし，構成式としてSchmid則に基づく単結晶の構成式を導入したものであるが，塑性ひずみ増分の決定のプロセスが連続体モデルとは異なる。結晶塑性では塑性ひずみ増分は実活動辷り系の辷り増分のSchmid テンソルの重み付き総和で定義されるが，実活動辷り系の組み合わせも未知であるので，現応力に対応する塑性ひずみ増分を事前に決定できないという厄介な問題がある。この解決法としてはPierceらによる速度依存型解法解法⁶⁴⁾とTakahashiによる逐次累積法⁶⁵⁾がある。後者は物理的に合理的な解法であるが陰解法であり，計算時間がかかる。一方，Pierceらの方法は単結晶の硬化則に巨視的クリープ則で広く採用されているNorton則によるひずみ速度依存性硬化則を採用し，すべての辷り系が常に活動しているとの前提で   

$$\left|{\dot{\gamma }}^{\left(i\right)}\right|={\dot{\gamma }}_0{\left|\frac{{\tau }^{\left(i\right)}}{{\tau }_c^{\left(i\right)}}\right|}^{1/m}$$ (27)

で評価している。ここにγ⁽ⁱ⁾，τ⁽ⁱ⁾，τ_c⁽ⁱ⁾はi番目の辷り系の辷り速度，分解せん断応力，せん断降伏応力，γ₀は代表辷り速度である。ここでクリープ硬化指数mを小さな値に設定すると，分解せん断応力が降伏せん断応力以下の辷り系の辷り速度は急激に小さくなり，事実上単結晶の降伏条件を満す解が得られる。これは活動辷り系の事前の選択を回避した陽解法であり，計算効率が良く簡便であるため，多くの多結晶有限要素解析で採用されている。しかしながら，クリープ硬化係数mの値は対象とする単結晶のひずみ速度依存性の実測値から評価したのではく，計算上の都合で恣意的に設定され物理的な意味はない。また，巨視応力の方向が急変する場合には活動辷り系の組み合わせがm に依存して変化するなどの問題がある。

ここに紹介した多結晶モデルをダイレクトに有限要素法などに導入した成形解析は材料特性に関する巨視的モデルの物理的な曖昧さを排除でき理想的であるが，現在の環境では現実的ではないと思われる。仮に1 mm×10 cm×10 cmの板材を多結晶有限要素で解析する場合，結晶粒径を100 μとして約3×10⁷元の連立1次方程式を解く必要があり，通常の計算機環境では不可能に近い。また，異方性材料に対しては集合組織情報も要求されるが材料のユーザー側では普通には入手困難である。それ故，現時点では結晶塑性学的視点から成形性に影響する巨視的な材料特性の物理的な解釈および巨視的構成式への導入法の道筋を示すことをマルチスケールモデルである多結晶モデルに期待したい。そんな示唆に富む高橋の著書⁶⁶⁾を紹介して20世紀後半の展望を終える。

4・2　結晶塑性有限要素法における最近の展開

第3章で概説した現象論モデルは，実験的に測定された巨視的な加工硬化挙動の数式化を目的としている。そのためそれまで想定されなかった加工硬化挙動を示す金属材料が開発されると，それに応じてモデルを再構築する必要がある。また前述のように異方性の発展を考慮することが難しく，広範囲な変形モードを正しく記述するには至っていない。例えば近年軽量化材料として注目されるマグネシウム合金は，引張時と圧縮時で加工硬化挙動が大きく異なり，また異方性の発展も著しい⁸⁰⁾。この特異な変形挙動を既存の現象論モデルで再現することは困難であり，またマグネシウム合金用に新しく提案されたモデル^81,82)も複雑かつパラメータ同定が困難などの問題がある。

一方最近，有限要素法および計算機の発展に伴って結晶塑性解析の研究が再び活発化している。4・1節でも述べられているように，結晶塑性解析技術は1980年代にはその基盤が既に確立されている⁶⁴⁾。歴史的には変形に伴う集合組織の発展予測に用いられることが多かったが，最近ではすべり系レベルの微視的な加工硬化挙動を適切にモデル化することで巨視的な加工硬化挙動の予測も可能になりつつある。結晶塑性モデルが現象論モデルと大きく異なる点は，加工硬化挙動や集合組織の発展とそれに伴う異方性の発展など，現象論モデルでは“入力値”として用いられてきた諸物理量が結晶塑性モデルでは予測可能な点である。また理論的には，塑性スピンのように現象論モデルでは取り扱いの困難な物理量が陽に表現できるのも利点の一つである。加えて，すべり系レベルの微視的な加工硬化挙動を転位密度に基づいてモデル化することも可能であり，転位論の知見が巨視的な変形予測に直接活用できるのも大きな特徴である。4・1節で詳述されているように，歴史的にはTaylorモデル⁵⁴⁾やセルフコンシステントモデル^56,57,58,59)が主流であった。最近では弾塑性構成式に代わって結晶塑性モデルを有限要素法へ導入した結晶塑性有限要素法が幅広く用いられており，加工成形解析との親和性も高まってきている。

歴史的には，活動しうるすべり系の種類が比較的少ない面心立方(FCC)金属を対象とした解析が主流であったが，最近では六方最密(HCP)金属や体心立方(BCC)金属を対象とした研究が活発である。鋼を対象とした研究事例を以下に紹介する。Kitayamaら⁸³⁾はIF鋼の弾塑性変形挙動を結晶塑性モデルにより解析し，定性的な傾向は精度良く予測できるもののバウシンガー効果やr値などの定量的な予測精度は未だ不十分であると結論づけている。Peetersら^84,85)はBCC金属における転位下部組織と集合組織の発展を考慮した転位密度ベースの結晶塑性モデルを提案し，IF鋼におけるバウシンガー効果や交差硬化，また転位下部組織の発展を再現できることを示した。Watanabeら⁸⁶⁾はフェライトーパーライト鋼の力学挙動予測を目的として，結晶塑性FEMをベースとした解析システムの構築を試みている。

一方HCP金属については，軽量化材料として注目されるマグネシウム合金への適用が盛んである。マグネシウム合金では，塑性変形に対して変形双晶が重要な役割を果たすことが知られている。変形双晶はその活動に極性があり，また一旦双晶回転した領域が反転負荷を受けることで回復する(元の方位へ戻る)ことも明らかにされている⁸⁰⁾。そのため結晶塑性解析の枠組みで様々な変形双晶モデルが提案されており^87,88)，それにより各種負荷経路における応力挙動やそれに伴う集合組織の発展を高精度に予測できることが報告されている。

以上のように結晶塑性モデルは，結晶レベルの微視的変形を適切にモデル化することで特別な仮定を設けることなく比較的シンプルに様々な力学特性を再現できる利点を持つ。そのため現象論モデルによる弾塑性構成式に代わって加工成形解析で用いることができれば，これまで以上の解析精度の向上と幅広い活用につながることが期待される。しかしながら現状ではいくつかの技術的な課題のためあまり普及していないのが現状である。第一の問題として，計算モデルのサイズが挙げられる。4・1節でも指摘されているように1要素1結晶粒としてモデリングすれば非現実的な要素数が必要になる。またTaylorモデルを各積分点に導入することで現実的な要素数で解析する方法も提案されている⁸⁹⁾が，有限要素法の特性が十分活用できているとは言い難い。これに対して最近では，均質化法を用いることで有限要素法の特性を活かしたまま各積分点へ多結晶塑性モデルを導入する手法が用いられている⁹⁰⁾。この手法では，各積分点に導入されるミクロスケール有限要素モデルも毎ステップ解く必要があるため，現実的な加工プロセスを解くにはやはり膨大な時間がかかる。しかしながら解析精度と計算時間の両面から最も実用に近い手法の一つと思われる。

また第二に，材料パラメータ同定の問題がある。例えばマグネシウム合金のようにすべり系によって臨界分解せん断応力や加工硬化が異なる材料では，膨大な数のパラメータを同定する必要がある。本来であればこれらのパラメータを単結晶材の変形に基づいて同定するのが理想的だが，現実問題として困難である。したがって通常は多結晶材の一軸引張変形で得られた応力−ひずみ曲線を参照解として同定する場合が多い。しかしながらこれだけでは同定されたパラメータの一意性が保証できず，その他の変形モードでは実験と全く異なる解が得られる場合がある⁹¹⁾。結晶塑性解析をより実用化に近づけるには，比較的容易かつ高精度にパラメータ同定できる手段を構築することが不可欠である。最近共著者らにより，マグネシウム合金における高精度なパラメータ同定法が提案された⁹¹⁾が，同様の枠組みを各種材料で確立することが重要であろう。

以上のように結晶塑性モデルは計算時間や使い勝手の面では現象論モデルには及ばず，直接加工成形解析へ活用することは現実的ではないのが現状である。一方，結晶塑性モデルを現象論モデルのパラメータ同定に利用するのは一つの現実的な活用法であろう。上述のように現象論モデルでは複雑化に伴って同定すべきパラメータ数も増えている。通常これらのパラメータは各種の材料実験により同定されるが，この材料実験には専用の装置を必要とするものも多い。そこでこの材料実験を結晶塑性モデルにより数値的に行うことができれば，計算機により簡便にパラメータを同定することができ，これまで以上に各種現象論モデルを有効活用することができる。このシステムの確立にも解決すべき技術課題が多いが，近い将来の目標として現象論モデルと結晶塑性モデルの両方を有効活用できる一つの方法論となりうるのではないかと考えている。

5. おわりに

本稿では20世紀後半から現在に至るまでの成形解析手法と材料モデリングの発展について，筆者らの私見を交えながら概説した。現在加工成形解析で活用されている基礎理論の大半は1980年代頃までには完成しており，以降はそれを土台にして応用が積み重ねられてきた。また数値解析技術についても計算機の発展と同期して2000年代に急成長し，現在は成熟の域に入っていると考えてよいだろう。一方，本文でも述べたように解決すべき課題は数多く残されており，未だ大きく発展する余地がある。特に，材料モデリングについては結晶塑性モデルをはじめとする従来の現象論モデルとは異なるアプローチが急発展しており，今後大きなブレークスルーを期待したい。

最後に，有限要素法のユーザーについて少し言及したい。現在では有限要素法ソフトウェアの高度化に伴って，かつてMakinouchi⁹²⁾が指摘したようにソフトウェアの開発者とユーザーが完全に分断されている状況である。その結果，同一ソフトウェアを使用して同一の問題を解いているつもりでも，ユーザーによって解析結果が大きく異なる事態が日常茶飯事となっている。高度化した数値解析技術を最大限に有効活用するためにはユーザーにも最低限の知識が求められるが，必ずしもその状況が整っていないのが現状のようである。その意味において，本稿が専門の研究者だけでなく一般のユーザーにとっても参考となれば望外の喜びである。なお，本稿では紙面の都合上参考文献には主要なもののみを厳選して掲載した。本稿で説明しきれなかった詳細についてはこれら原著文献を参照いただきたい。

文献

• 1)   R.  Hill: Mathematical Theory of Plasticity, Oxford at the Cralenton Press, (1950).
• 2)  鷲津久一郎，山田嘉昭，工藤英明 共訳:塑性学，培風館，東京，(1954).
• 3)   C.  Truesdell and  R.  Toupin: The Non-Linear Field Theories of Mechanics, Springer-Verlag, (1960).
• 4)   C.  Truesdell and  W.  Noll: The Non-Linear Field Theories of Mechanics, Springer-Verlag, (1964).
• 5)   G.  Sacks: Zeit. Ver. deut. Inge., 72(1928), 734.
• 6)   G.I.  Taylor: J. Inst. Met., 62(1938), 307.
• 7)   R.  Hill: J. Mech. Phys. Solids, 4(1956), 245.
• 8)  材料加工の計算力学，日本塑性加工学会編，コロナ社，東京，(1992), 80.
• 9)   O.C.  Zienkiewicz: The Finite Element Method in Structural and Continuum Mechanics, McGraw -Hill, New York, (1967).
• 10)   M.J.  Turner,  R.W.  Clough,  H.C.  Martin and  L.J.  Topp: J. Aero. Sci., 3(1956), 805.
• 11)   R.  Hill: J. Mech. Phys. Solids, 7(1959), 209.
• 12)   R.M.  MacMeeking and  J.R.  Rice: Int. J. Solids Struct., 6(1975), 11.
• 13)   R.D.  Krieg and  D.R.  Krieg: Trans. ASME, J. Press. Vessel Tech., 99(1977), 51.
• 14)   S.  Gosh and  N.  Kikuchi: Int. J. Eng. Sci., 26(1988), 143.
• 15)   M.  Gotoh: J. Jpn. Soc. Technol. Plast., 27(1986), 25.
• 16)  静的解法FEM−板成形，日本塑性加工学会編，コロナ社，東京，(2004).
• 17)   Y.T.  Keum,  E.  Nakamachi,  R.H.  Wagoner and  J.K.  Lee: Int. J. Numer. Methods Engrg, 30(1990), 1471.
• 18)   M.  Kawka and  A.  Makinouchi: J. Mater. Process. Technol., 50(1995), 105.
• 19)   T.  Hama,  T.  Nagata,  C.  Teodosiu,  A.  Makinouchi and  H.  Takuda: Int. J. Mech. Sci., 50(2008), 175.
• 20)   Y.  Yamada,  N.  Yoshimura, and  T.  Sakurai: Int. J. Mech. Sci., 10(1968), 343.
• 21)   N.  Yamamura,  T.  Kuwabara,  A.  Makinouchi and  C.  Teodosiu: J. Jpn. Soc. Technol. Plast., 43(2002), 432.
• 22)   M.C.  Oliveira,  J.L.  Alves and  L.F.  Menezes: Int. J. Numer. Methods Engrg, 58(2003), 2083.
• 23)   T.  Hama,  K.  Iguchi,  H.  Hishida and  H.  Takuda: Metal Forming Conference 2008, Steel Research International 79, Special Edition, 2(2008), 587.
• 24)   D.Y.  Yang,  D.W.  Jung,  I.S.  Song,  D.J.  Yoo and  J.H.  Lee: J. Mater. Proc. Technol., 50(1995), 39.
• 25)   T.  Hama: J. Jpn. Soc. Technol. Plast., 51(2010), 1017.
• 26)   K.S.  Al-Athel and  M.S.  Gadala: Comput. Methods Appl. Mech. Engrg, 200(2011), 2145.
• 27)   M  Nakazaki,  I.  Takasu and  H.  Utsunomiya: J. Jpn. Soc. Technol. Plast., 50(2009), 349.
• 28)   R.  Bonan and  J.P.  Ponthot: Int. J. Numer. Methods Engrg, 92(2012), 857.
• 29)   P.W.  Cleary,  M.  Prakash and  J.  Ha: J. Mater. Process. Technol., 177(2006), 41.
• 30)   I.  Alfaro,  J.  Yvonnet,  E.  Cueto,  F.  Chinesta and  M.  Doblare: Comput. Methods Appl. Mech. Engrg., 195(2006), 6661.
• 31)   S.  Xiong and  P.A.F.  Martins: J. Mater. Process. Technol., 177(2006), 49.
• 32)   M.A.  Puso and  T.A.  Laursen: Int. J. Numer. Methods Engrg, 54(2002), 1161.
• 33)   T.  Nagata: Comput. Aided Geom. Design, 22(2005), 327.
• 34)   T.  Hama,  M.  Asakawa,  S.  Fuchizawa and  A.  Makinouchi: Mater.Trans., 44(2003), 940.
• 35)   R.P.R.  Cardoso,  J.W.  Yoon,  M.  Mahardika,  S.  Choudhry,  R.J.  Alves de Sousa and  R.A.F.  Valente: Int. J. Numer. Meth. Engrg, 75(2008), 156.
• 36)   R.  de-Carvalho,  S.  Silva,  R.A.F.  Vanente and  A.  Andre-Campos: Finite Elem. Anal. Design, 61(2012), 75.
• 37)   H.  Ou,  P.  Wang,  B.  Lu and  H.  Long: Comput. Struct., 90-91(2012), 13.
• 38)   K.  Ito: J. Jpn. Soc. Technol. Plast., 25(1984), 577.
• 39)   B.  Budiansky: Trans. ASME, J. Ser. E, 26(1959), 259.
• 40)   S.  Storen and  J.R.  Rice: J. Mech. Phys. Solids, 23(1975), 421.
• 41)   J.  Christoffersen and  J.W.  Hutchinson: J. Mech. Phys. Solids, 27(1979), 465.
• 42)   M.  Gotoh: Trans. JSME, Ser.A, 47(1981), 1389.
• 43)   K.  Ito,  K.  Satoh,  M.  Goya and  T.  Toshida: Int. J. Mech. Sci., 42(2000), 2234.
• 44)   R.  Hill: Proc. R. Soc. London, A, 193(1948), 281.
• 45)   R.  Hill: Math. Proc. Camb. Phil. Soc., 85(1979), 179.
• 46)   R.W.  Logan and  W.F.  Hosford: Int. J. Mech. Sci., 22(1980), 419.
• 47)   F.  Barlat and  J.  Lian: Int. J. Plast., 5(1989), 51.
• 48)   F.  Barlat,  J.C.  Brem,  J.W.  Yoon,  K.  Chung,  R.E.  Dick,  D.J.  Lege,  F.  Pourboghrat,  S.-H.  Choi and  E.  Chu: Int. J. Plast., 19(2003), 1297.
• 49)   M.  Gotoh: Trans. JSME, Ser.A, 48(1982), 498.
• 50)   W.  Pragar: Proc. Inst. Mech. Eng., 169(1955), 41.
• 51)   H.  Zieger: Quart. J. Appl. Math., 17(1959), 55.
• 52)   Z.  Mroz: J. Mech. Phys. Solids, 15(1967), 163.
• 53)   F.  Yoshida and  T.  Uemori: Int. J. Mech. Sci., 45(2003),1687.
• 54)   G.I.  Taylor: J. Inst. Met., 62(1938), 307.
• 55)   J.F.W.  Bishop and  R.  Hill: Philos. Mag., 42(1951), 414.
• 56)   J.D.  Eshelby: Proc. R. Soc. London, A, 241(1957), 376.
• 57)   E.  Kröner: Acta. Metall., 9(1961), 155.
• 58)   B.  Budiansky and  T.T.  Wu: Proc. 4th Congr. Appl. Mech., 2(1962), 1175.
• 59)   T.H.  Lin: Trans. ASME, J. Appl. Mech., 34(1961), 606.
• 60)   K.  Ito and  H.  Takeyama: J. Jpn. Soc. Technol. Plast., 23(1982), 410.
• 61)   H.  Takahashi,  K.  Ito and  M.  Goya: Int. J. Plast., 6(1990), 615.
• 62)   H.  Miyamoto,  Y.  Ishijima and  T.  Miyoshi: J. Jpn. Soc. Mech. Eng., 75(1972), 602.
• 63)   T.  Jimma,  T.  Murota and  T.  Ichiyanagi: J. Jpn. Soc. Mech. Eng., 75(1972), 575.
• 64)   D.  Peirce,  R.J.  Asaro and  A.  Needleman: Acta. Metall., 31(1983), 1951.
• 65)   H.  Takahashi,  H.  Motohashi,  M.  Tokuda and  T.  Abe: Int. J. Plast., 10(1994), 63.
• 66)  高橋寛:多結晶塑性論，コロナ社，東京，(1999).
• 67)   T.  Kuwabara,  R.  Enatsu,  S.  Yamagishi and  F.  Sugawara: Tetsu-to-Hagané, 98(2012), 61.
• 68)   K.  Hashimoto,  T.  Kuwabara,  E.  Iizuka and  J.W.  Yoon: Tetsu-to-Hagané, 96(2010), 27.
• 69)   F.  Yoshida,  H.  Hamasaki and  T.  Uemori: Int. J. Plast., 45(2013), 119.
• 70)   O.  Cazacu,  B.  Plunkett and  F.  Barlat: Int. J. Plast., 22(2006), 1171.
• 71)   D.  Yanaga,  H.  Takizawa and  T.  Kuwabara: J. Jpn. Soc. Technol. Plast., 55(2014), 55.
• 72)   F.  Yoshida and  T.  Uemori: Int. J. Plast., 18(2002), 661.
• 73)   C.  Teodosiu and  Z.  Hu: Proc. the NUMIFORM’95, ed. by S.F. Shen and P.R. Dawson, Balkema, Rotterdam, (1995), 173.
• 74)   F.  Barlat,  J.J.  Gracio,  M.G.  Lee,  E.F.  Rauch and  G.  Vincze: Int. J. Plast., 27(2011), 1309.
• 75)  久田俊明，野口裕久:非線形有限要素法の基礎と応用，丸善，東京，(1996).
• 76)   T.J.R.  Hughes and  J.  Winget: Int. J. Numer. Methods Engrg, 15(1980), 1862.
• 77)   J.W.  Yoon,  D.Y.  Yang and  K.  Chung: Comput. Methods Appl. Mech. Engrg, 174(1999), 23.
• 78)   M.  Kawka and  A.  Makinouchi: J. Mater. Process. Technol., 60(1996), 239.
• 79)   L.  Duchêne,  T.  Lelotte,  P.  Flores,  S.  Bouvier and  A.-M.  Habraken: Int. J. Plast., 24(2008), 397.
• 80)   X.Y.  Lou,  M.  Li,  R.K.  Boger,  S.R.  Agnew and  R.H.  Wagoner: Int. J. Plast., 23(2007), 44.
• 81)   M.G.  Lee,  S.J.  Kim,  R.H.  Wagoner,  K.  Chung and  H.Y.  Kim: Int. J. Plast., 25(2009), 70.
• 82)   M.  Li,  X.Y.  Lou,  J.H.  Kim and  R.H.  Wagoner: Int. J. Plast., 26(2010), 820.
• 83)   K.  Kitayama,  T.  Kobayashi,  T.  Uemori and  F.  Yoshida: Tetsu-to-Hagané, 97(2011), 57.
• 84)   B.  Peeters,  M.  Seefeldt,  C.  Teodosiu,  S.R.  Kalidindi,  P.  van Houtte and  E.  Aernoudt: Acta Mater., 49(2001), 1607.
• 85)   B.  Peeters,  B.  Bacroix,  C.  Teodosiu,  S.R.  Kalidindi,  P.  van Houtte and  E.  Aernoudt: Acta Mater., 49(2001), 1621.
• 86)   I.  Watanabe,  D.  Setoyama,  N.  Nagasako,  N.  Iwata and  K.  Nakanishi: Int. J. Numer. Methods. Engrg, 89(2012), 829.
• 87)   H.  Wang,  P.D.  Wu,  C.N.  Tome and  J.  Wang: Mater. Sci. Eng. A, 555(2012), 93.
• 88)   T.  Hama,  N.  Kitamura and  H.  Takuda: Mater. Sci. Eng. A, 583(2013), 232.
• 89)   T.  Walde and  H.  Riedel: Acta Mater., 55(2007), 867.
• 90)   E.  Nakamachi,  N.N.  Tam and  H.  Morimoto: Int. J. Plast., 23(2007), 450.
• 91)   T.  Hama,  N.  Hosokawa and  H.  Takuda: Proc. NUMISHEET2014, Eds. J.W.Yoon, B.Rolfe, J.H.Beynon and P.Hodgson, AIP Publishing, New York, (2014), 692.
• 92)   A.  Makinouchi: J. Jpn. Soc. Technol. Plast., 40(1999), 414.