数値材料試験による異方性非線形材料特性の算出 〜 その2：繊維強化プラスチックのマルチスケール解析 〜

東北大学
寺田賢二郎
日東紡績株式会社
平山紀夫

1．はじめに

前稿[1]では、繊維強化プラスチック（以下、FRP）の微視的（ミクロ）構造と巨視的（マクロ）構造の力学的相互作用の評価を意図して、マクロ材料特性を把握するための「数値材料試験」、すなわちミクロ構造に対する“均質化解析”に用いるプラスチック（エポキシ樹脂）の材料試験データの計測手順と、データ処理・分析について説明しました。本稿では、このプラスチックを母材に用いたFRPのマクロ的な材料挙動を表現可能な異方性粘弾性構成則を提示し、その材料パラメーターを数値材料試験により同定する方法を解説します。また、そのFRPの供試体をマクロ構造物とするクリープ試験の“マクロ解析”を行い、実験結果と比較することで提示した構成則を検証するとともに、その変形履歴を反映したミクロ構造の応答を再現する“局所化解析”の一例を紹介します。

この均質化解析−マクロ解析−局所化解析という一連の解析を（均質化理論に基づく）「マルチスケール解析」と呼んでいますが、その理論の詳細については文献[2]等を、数値材料試験の考え方や方法については文献[3]〜[5]等を参照して下さい。

2．異方性粘弾性構成則

前稿[1]では、一般に広く利用されている等方性一般化Maxwellモデルを、FRPミクロ構造内の母材（プラスチック）の粘弾性特性として用いて数値材料試験を行いました。その数値材料試験で得られたマクロ応力緩和応答が異方性を示したことをから、本稿では等方性一般化Maxwellモデルを異方性粘弾性体に拡張したモデル[6]を紹介し、マクロ構成則として採用することにします。

2.1　“等方性”一般化Maxwellモデル

Ansys等の汎用FEMコードでは一般に、一般化Maxwellモデルが等方的な粘弾性モデルとして用意されています。このモデルは図1（a）に示すような“ばね−ダッシュポット”系を想定して、時間 t に関する緩和（せん断）弾性率を

と与えます。ここで、 は平衡応力に関するせん断弾性係数、 は初期せん断弾性率、 と はそれぞれMaxwell要素αのせん断弾性係数と粘弾性係数です。また、緩和時間 は次式で定義されています。

このモデルは、前稿1）でプラスチックの材料モデルとして採用したものであり、同定された材料パラメーターは数値材料試験に際しての入力データとなります（文献[1]の図5を参照）。
通常は、応力の静水圧成分は粘弾性挙動を示さないと仮定するので、緩和試験を想定して一定ひずみ を与えると、応力 は、

と計算されます。ここで、添え字“T”は行列・ベクトルの転置を表し、Kは（時間に依存しない）体積弾性係数です。また、 i やI₀ は次式で定義されています。

この式（3）に緩和せん断弾性率（1）（第1式）を代入すると、次式のように応力の平衡部分 と非平衡部分 （時間とともに緩和する成分）に分解できます。

図1　一般化Maxwellモデル

2.2　“異方性”一般化Maxwellモデルへの拡張

前稿[1]では、FRPの母材であるプラスチックには前節で示した等方性一般化Maxwellモデルを用いて、繊維は粘弾性挙動を示さない弾性体とみなしました。そして、そのFRPミクロ構造に対して数値材料試験を行うと、マクロ材料挙動としては異方的な粘弾性応答が得られることを、前稿[1]の図7で示しました（図2）。本節では、そのような異方性粘弾性挙動を表現できる構成則として、前節の一般化Maxwellモデルを拡張したモデルを紹介します。

拡張に際して、まず「FRPの異方性粘弾性挙動で観察される緩和時間は、母材の緩和時間と同一である」という仮説を導入することにします。この仮定は、一見大胆に思えますが、「応力は変形と弾性率のみで決まる」という原則に照らせば妥当な仮説だといえます。すなわち、緩和せん断弾性率（1）から、初期と平衡時の応力は弾性率のみで決定されることが分かり、緩和弾性率は非平衡時の応力緩和挙動のみを司ります。そして、緩和を特徴づける緩和時間を有する材料は母材のプラスチックのみであり、繊維は緩和しません。このことは、マクロな緩和挙動は繊維の有無に無関係であり、母材のそれを継承することを意味します。

次に、図1（b）に示すように等方性一般化Maxwellモデルのスカラー値の弾性係数・粘性係数を、異方性を表す行列形式の係数に置き換えると、異方性緩和弾性率が、

と表されます。 は、それぞれ平衡応力に関する異方性弾性係数行列、初期の異方性弾性係数行列、Maxwell要素の異方性弾性係数行列であり、 6×6次の対称行列です。例えば直交異方性を仮定すると、緩和弾性係数行列の成分は

で与えられます。そして、一定ひずみ を与えたときの応力の時間変化は、

となり、この応答が前稿（1）図7の緩和曲線に対応します（図2）。なお、緩和時間は は、上記の仮説により母材の緩和時間（式（2））と同じであると仮定したスカラー値ですが、Maxwell要素αの異方性粘弾性数行列を と表すと のように定義したとものと解釈できます。ここで、 I は6×6の単位行列です。

式（5）と（8）と比較すれば分かりますが、本稿で提示した異方性一般化Maxwellモデルは、等方的な一般化Maxwellモデルにおける 1＋n 個のスカラー値のせん断弾性率を、1＋n 個の異方性弾性行列に代えただけのものです。しかし、残念ながらAnsysを含む多くの汎用FEMソフトには異方性粘弾性モデルが具備されていないようですので、今回はこのモデルをユーザーサブルーティンプログラムとしてAnsysに実装し、FRPのマクロ構成則として用いることにします。

図2　 一方向強化複合材料の数値応力緩和試験の結果
（「CAEのあるものづくり Vol.10」 p12 図7に掲載。ZZ方向に繊維が配向している。）

3．異方性粘弾性構成則のパラメーター同定

3.1　数値材料試験

数値材料試験[3]〜[5]とは、均質化法に基づくマルチスケール解析に際して、マクロ構成則のパラメーターを同定するために行う、複合材料の単位周期ミクロ構造（以下、ユニットセル）のFEMモデルを“数値”供試体とみたてた数値実験のことです。ですから、FEMの計算結果として得ようとするものは、ミクロの変形や応力状態ではなく、マクロひずみおよび応力の応答データです。均質化法では、この数値実験プロセスのことを「均質化」、対応するユニットセルに対する解析を「均質化解析」と呼んでいます。

図3に示すような1方向繊維強化複合材料（FRP）について、そのミクロ構造（ユニットセル）内のエポキシ樹脂（プラスチック）に、前稿で実験的に求めた粘弾性材料パラメーターを、繊維には適当な弾性係数を入力データとして与え、周期境界条件の下で各軸垂直方向にマクロ的な一軸引張、各面内に純せん断を与える数値材料試験（全6ケース）を実施しました。得られたマクロ応力緩和曲線を図4に示します。
ここで、前節で参照した座標系とユニットセルに設定した座標系の関係は y₁= x₁, y₂= x₂, y₃= x₃ となっています。この図から、等方性粘弾性体の母材と等方性弾性体の繊維からなるFRPの材料挙動が、繊維の配向を反映して異方的な緩和挙動を示すことが分かります。なお、繊維軸（y₃ = x₃軸）方向はほとんど緩和しないので は示していません。

図3　 1方向繊維強化複合材料（FRP）のマクロ構造と単位周期ミクロ構造
（ユニットセル）

図4　数値材料試験結果

3.2　パラメーター同定

本項では、前項で求めた数値材料試験結果を用いて、前節で示した異方性粘弾性構成則のパラメーターを同定します。
まず、前節の式（6）で与えた異方性緩和弾性係数行列 D_r の成分を改めて

と書くことにします。ここで、 と定義し、

とおいています。前項のユニットセルのように、マクロ材料挙動が直交異方性の場合、同定しなければならないDrの独立な成分は、各  について9個 あります。したがって、同定すべきパラメーター数は9×（n＋1）個です。

前項で行った各マクロひずみパターンに対する数値材料試験では、マクロ緩和弾性係数行列（7）の各列が得られていますから、その全成分がミクロ解析の時間ステップごとに“計測”されていることになります。いま、マクロひずみパターンm、時刻tpで計測された緩和弾性率の成分を 、対応して式（9）から算出される異方性緩和弾性係数行列の成分を と書くことにします。そして、これら と から作る二乗誤差関数

を定義して、この停留条件、すなわちパラメーター についての微分がゼロ を求めると、 に関するn＋1元連立一次方程式

が得られます。この方程式を満たす は、誤差関数の最小値（理想的にはゼロ）を与えるもので、各成分IJ（直交異方性の場合は先に挙げた9成分）について同様の方程式を解けば、所望の材料パラメーターのすべてが同定されることになります。

このようにして同定した材料パラメーターを用いて緩和弾性係数行列を式（9）から算出し、前項で数値材料試験により“計測”したマクロな応力緩和応答と重ねて緩和曲線を描くと図5のようになります。前節で示した異方性粘弾性構成則が、同定したパラメーターを用いることで、数値材料試験で得られた材料挙動を適切に再現していることが分かります。

図5　数値材料試験結果

4．マルチスケール解析 〜マクロ解析とミクロ解析〜

第2節で示した異方性粘弾性構成則に前節で同定した材料パラメーターを用いて、実際に実験を行った3点曲げクリープ試験を再現するための「マクロ解析」を行い、次いで「局所化解析」により対応するミクロ構造の応答を評価してみます。

4.1　マクロ解析

長手方向に対して炭素繊維が15°および45°に配向されているFRPの試験片を2種類（15°補強材および45℃補強材）用意して、クリープ試験をJIS K 7166の試験方法に準じて行いました。対応するマクロ解析には、図6に示すようなモデル・負荷／支持条件を設定し、マクロ要素の局所座標をユニットセルの座標系に一致させました。ここで、繊維方向は、図4に示したユニットセルの座標系のy3軸の、マクロ構造の座標系のx1軸に対する角度です。また、解析温度は40℃と60℃の2水準を設定し、温度ー時間換算則にはWLF式を採用しました。

マクロ解析結果として、クリープ解析の初期ステップにおける相当応力のコンター図を図7に示します。母材であるエポキシ樹脂より剛な材料である炭素繊維の配向を反映した応力分布になっており、ユーザーサブルーティンとして実装した異方性粘弾性構成則が適切に機能していることが分かります。また、実験および解析で得られた各解析温度におけるクリープひずみの時刻歴応答を図8に示します。実験条件と解析条件の不一致に起因して多少ずれていますが、マクロ解析はおおむね実験を再現しているといえます。

図6　FRPのクリープ試験のための供試体の有限要素モデル（マクロ構造）
※ 力Fは15°モデルには320[N]、45°モデルには70[N]を負荷

図7　負荷直後のマクロ応答

4.2　局所化解析

次に、マクロ解析の結果として得られるマクロ構造内の任意の点におけるマクロ変形履歴を、ユニットセルに与えて解析することで、マクロな応答に対応するミクロ構造の力学的挙動を評価します。均質化理論の枠組みでは、このミクロ構造に対する解析のことを「局所化解析」と呼んでいます。

15°補強材を用いた試験片に対するマクロ解析結果（設定温度60℃）について、最も大きな軸方向引張応力が作用した要素の積分点を選定し、そこでのマクロ変形履歴（ただし局所座標系を参照した成分）をマクロ解析の出力データから抽出します。これを文献5）に詳しく解説されている、外部制御節点を用いる方法でユニットセルに与えて、ミクロ解析を行いました。解析結果として、ミクロ領域におけるvon Mises応力分布をマクロ解析結果と併せて図9示します。このように、均質化法に基づくマルチスケール解析によれば、実際の使用環境下で粘弾性挙動を示しているFRP製品の任意の場所で、材料内部の力学現象を評価できるので、FRPのミクロ構造の設計支援も可能となります。

5．おわりに

本稿では、前稿における樹脂の粘弾性特性の計測方法、パラメーター設定方法についての解説、ならびに数値材料試験結果を受けて、これを母材とするFRPのマルチスケール解析について解説しました。
ここで紹介した一連の解析は、原理的にはすべて汎用FEMコードで実行可能であり、特にAnsysにはこの機能を追加して利用できる専用ツール7）も用意されています。

本稿で提案した解析方法では、汎用FEMソフトに具備されている等方性の一般化Maxwellモデルを異方性モデルに拡張した粘弾性構成則を使用しています。したがって、汎用FEMソフトで普段から粘弾性解析を行っているユーザーには、比較的なじみやすい手法だと思います。実際に解析を行ってみると、各材料軸方向で大きく異なる材料挙動が表現でき、実用的にも有用であることが実感できます。しかし、その材料パラメーターを実際の材料試験から計測するのでは、試験片の準備の大変さ、そして計測方法の困難さが故に、せっかくの異方性粘弾性構成則も利用し易いとは言えません。そのような問題に直面したときには、それを解決する方法を紹介した前稿と本稿を、是非いま一度、読み返して見てください。

また、本稿でも利用したユーザーサブルーティンという汎用FEMソフトの拡張機能は、“CAEのあるものづくり”に新たな付加価値を与えてくれます。本稿では紙面の関係で説明を省略しましたが、機会があれば解説してみたいと思います。

図8　マクロクリープ曲線の比較

図9　局所化解析結果

謝辞

本稿は、東北大学大学院の濱名康彰氏（現：JFEエンジニアリング）の研究成果をもとに著しました。ここに深く感謝の意を表します。