Title
Mn基合金の内部摩擦を考慮した振動シミュレータの検討
Author(s)
鞆田 顕章
Citation
福岡工業大学総合研究機構研究所所報 第1巻 P13-P17
Issue Date
2018-12
URI
http://hdl.handle.net/11478/1215
Right
Type
Departmental Bulletin Paper
Textversion
Publisher
福岡工業大学 機関リポジトリ
FITREPO
Mn 基合金の内部摩擦を考慮した振動シミュレータの検討
鞆田 顕章（工学部知能機械工学科）
Study on the Vibration Simulator Considering Internal Friction in Manganese-Base Alloys
Akinori TOMODA （Department of Intelligent Mechanical Engineering, Faculty of Engineering）
Abstract
Accurate estimation method of the vibration response for a mechanical system using Molecular dynamics method (MD) is considered. This study focuses on the MD modeling for estimating energy dissipation in Manganese-based alloys under vibration. The phenomenon of energy dissipation in metal is caused by internal friction. In manganese-based alloys such as M2052 alloy, internal friction occurs by the formation of twins in the material under external force. In this paper, the result of experimental modal analysis for obtaining values of natural frequency and damping ratio of M2052 plate is reported. In addition, the development of a large-scale MD simulator for estimating of displacement and velocity of atoms in pure Manganese single crystals under a very low temperature is investigated.
Keywords：Damping Alloy, Vibration, Twin, Molecular dynamics simulation, Modified embedded atom method, Experimental modal analysis
1. 緒言
機械システムにおいては，重要部品の損傷や騒音発生を 防止する観点から，振動入力時に系の応答が著しく増大し ないように設計する必要がある．しかし，自動車やハード ディスクドライブ，原子間力顕微鏡等，多数の部品で構成 される機械システムの設計開発においては，その仕様上， やむを得ず振動特性の異なる金属材料を用いた部品を多数 使用しなければならず，振動入力時において系の応答が増 大する可能性を有する．このような系の振動特性の推定に は，FEM（有限要素法）が用いられることが多い．FEM に よる振動解析では，系の特性行列（質量，剛性，減衰行列） を予め求めておく必要がある．特に，減衰行列は振動入力 時における系の応答値に大きく影響を及ぼすため，正確に 決定するための手法が求められる． FEM における系の減衰行列の決定手法に関する過去の研 究では，対象構造物の周波数応答の実験結果から同定する 手法(1)，系の1 周期あたりのエネルギー散逸量がモード減衰 比で表現されるものと同値になるように同定する手法(2)が 提案されている．しかし，上述の手法は金属材料における 振動エネルギーの散逸機構を考慮したものではないため， FEM による複雑な機械システムの振動応答の高精度推定は 困難であり，CAE（Computer Aided Engineering）よる最適設 計の障壁となっている． 近年，機械構造物の挙動を推定するためのシミュレーシ ョン手法として，MD（分子動力学法）が注目されている． MD は，材料内部に存在する原子，分子の挙動についてニュ ートンの運動方程式を用いて計算する手法であり，新材料 開発や超精密加工シミュレーション等での実績を有する (3-6)．ここで，金属材料の振動減衰能の発現機構は，材料内 部における転位や双晶などの内部摩擦に起因する．金属材 料の結晶構造をモデリングすることが可能なMD を用いる ことで，その振動減衰能を精度良く推定することができる ようになると考えられる．最近では，転位生成および粒界 での転位の移動が超微細粒金属の機械特性に与える影響を 明らかにするため，約 1 億個の原子を用いたアルミニウム 多結晶のMD モデルを構築し，塑性変形時の材料内部の状 態の変化を考察できたことが報告されている(7)．MD のモデ リング手法に関する研究が進めば，金属材料の内部構造に 基づく新たな振動解析法を確立することが可能となり， CAE による複雑な機械システムにおける最適設計の高度化 が期待できる． そこで，本研究では，MD を用いて金属材料内部の転位や 双晶などの内部摩擦を高精度に再現することを目的とし， 金属材料の振動減衰能を高精度に推定することができる振 動シミュレータを開発する．本報告では，マンガン基制振 合金であるM2052 合金を解析対象とし，制振合金板の振動 特性の取得および大規模 MD 解析プログラムの構築につい て述べる．2. 制振合金板の振動特性
〈2･1〉 M2052 合金の振動減衰能発現機構 制振合金における振動減衰能の発現機構においては，複 合型，強磁性型，転位型，双晶型の4 種類に分類される(8)．鞆田 顕章 本研究では，マンガン基制振合金であるM2052 合金（Mn： 73%，Cu：20%，Ni：5%，Fe：2％）を解析対象とし，MD モデルを構築する．この合金は双晶型に属しており，マル テンサイト変態により合金内部に双晶が発生する．本合金 は，外力の作用時において材料内部に双晶が発生し，除荷 時に発生した双晶は消失する．このとき，振動エネルギー の一部が熱エネルギーに変換され，材料の振動低減がなさ れる．また，本合金の内部に発生する双晶は大小様々なサ イズであり，幅広い入力振動数に対して振動減衰能を発生 させることが可能である．本合金の減衰能の温度依存性お よび周波数依存性については，川原らにより実験的に明ら かにされている(9)．特に，周波数依存性については，入力周 波数によりわずかに変化が見られるものの，概ね少ないと いう報告がなされている．本報告では，M2052 合金製薄板 である WAVLES1（リックス（株））を用いて，打撃試験に より薄板の固有振動数および減衰比を求めることで本合金 の振動減衰能を取得する．取得した実験データは，開発す る振動シミュレータの妥当性検証時に利用する． 〈2･2〉 M2052 合金製薄板の振動減衰能取得試験 前述のM2052 合金製の薄板の振動減衰能を取得するため の実験装置を図1 に示す．薄板は幅 300mm，高さ 500mm， 厚さ0.1mm であり，支持構造物の上部からナイロン製の釣 り糸により薄板を吊り下げ，自由支持となるようにした． 薄板への振動入力は，実験モード解析用インパクトハンマ （PCB Piezotronics 製 086C03）を用いて，予め設定した薄板 表面の打撃点44 点（図 2）について各 5 回ずつ打撃した． 薄板の振動応答については，薄板の下方に取り付けた加速 度ピックアップ（PCB Piezotronics 製 352A24）により取得し た．インパクトハンマおよび加速度ピックアップをFFT ア ナライザ（小野測器製DS-3000）に接続し，薄板の伝達関数 を取得した．各打撃点における伝達関数を用いて曲線適合 を行い，薄板の固有振動数，減衰比を算出した．打撃試験 により得られた薄板の各次の固有振動数fn，減衰比n を表 1 に示す．本試験においては，1～10 次モードにおける各モー ドの固有振動数，減衰比を取得することができた．各次の 減衰比は概ね1.2～2.5％程度であることが分かった． 表1 M2052 合金製薄板の固有振動数および減衰比 Modes Natural frequency
fn [Hz] Damping ratio n [%] 1 16.34 2.54 2 21.24 2.04 3 35.70 2.43 4 56.16 1.42 5 63.41 1.81 6 74.93 1.69 7 90.83 1.48 8 107.9 1.26 9 126.9 1.56 10 156.7 1.73 図1 振動減衰能取得試験 図2 M2052 合金製薄板および打撃点 次に，薄板の1 次の減衰比1を用いて対数減衰率を求 めた．1との関係は式(1)で表される(10)． 1 1 2 1 2 1 ... (1) 上式および表1 より= 0.160 が得られた．本研究で使用し たM2052 合金製薄板の対数減衰率は参考文献(9)に記載され ている値（0.1～0.2 程度）とほぼ同値であることが分かった． ただし，本研究で用いたインパクトハンマでは，系に大き な振動入力を与えることが困難であるため，本報で示した 各振動モードにおける減衰比はいずれも小ひずみ領域にお ける値であるということに注意する必要がある．薄板の振 動減衰能の周波数依存性およびひずみ依存性を詳細に調査 Accelerometer (PCB 352A24) Impact hammer (PCB 086C03) String M2052 plate 11 12 13 14 15 1 2 3 4 5 6 7 8 9 10 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 Accelerometer (PCB 352A24)
するためには，電磁型加振器のように加振力および入力波 の周波数特性の調節が可能な装置を用いて薄板を加振する 必要がある． 以上より，小ひずみ領域におけるM2052 合金の薄板の振 動減衰能を取得することが出来た．
3. 分子動力学法による解析モデルの構築
〈3･1〉 分子動力学法について MD は，個々の原子に作用する力に基づき，ニュートンの 運動の第 2 法則を適用し原子挙動を追跡する方法である (11-12)．解析対象における系内の各原子の運動方程式は，以 下のように表される． 2 2 i id i m dt R f ... (2) ここで，miおよびRiはi 番目の原子の質量と位置ベクトル， fiは原子i に作用する力であり，系のポテンシャルエネルギ ーEtotを用いて，以下のように表すことができる． tot i i tot i E E f R ... (3) ここで，iは勾配を与える原子位置Ri (xi, yi, zi)に関する微 分演算子であり，直交座標系においてx, y, z 軸の正方向の単 位ベクトルをそれぞれex，ey，ezとして，以下のように表す ことができる． i x y z i i i i x y z e e e R ... (4) 系のポテンシャルエネルギーを求めるための原子間ポテン シャルは，解析対象の原子および結合の種類に対応したも のが多数開発されている．本研究では，Mn，Cu，Ni，Fe および多元系合金を扱うことができる MEAM（Modified Embedded Atom Method）ポテンシャルを用いて解析モデル の構築を行った．MEAM ポテンシャルは Baskes らによって 開発され，金属結合に関わる原子近傍の電子に起因するエ ネルギー多体項と，二体相互作用の項との和で表される(13)．MEAM ポテンシャルでは，系のポテンシャルエネルギーEtot
は次式により与えられる．
1 2 EAM tot i i ij ij i j i E F R
... (5) ここで，
iは原子i の位置の電子密度，FiEAMは電子密度 が
iの位置に原子i を埋め込むのに必要なエネルギーであ る．ijは距離がRijの2 原子間に作用する相互作用エネルギ ーである．式(2)～(5)を用いて各原子に作用する力を計算す ることで，原子の挙動を追跡することが可能となる．なお， Cu，Ni，Fe の 3 元素のポテンシャルパラメータについては Baskes らの報告に記載されている値(13)を，Mn については Kim らの報告に記載されている値(14)を用いることにした． ここで，本報の著者は，過去の研究において銅原子 500 個 程度の小規模 MD 解析モデルを用いて系の表面エネルギー の導出を行った(15)．本研究では，上述のMD 解析モデルを 構築するための数値解析用プログラムを基に，数万～数千 万原子の挙動が計算可能となるよう，プログラムの再構築 を行う． 〈3･2〉 大規模 MD 解析プログラムの構築 大規模MD 解析プログラムについては，Fortran 言語によ り記述した．実行ファイル生成のためのコンパイラについ ては，全国共同利用が可能なスーパーコンピュータ（九州 大学，Fujitsu PRIMERGY CX2550/CX2560 M4，Intel Xeon Gold 6154，各ノードに 2CPU 搭載）に対応，かつスレッド 並列およびプロセス並列用実行ファイルを容易に生成でき るIntel Fortran Compiler を採用した．Mn 原子を解析対象とし，結晶構造の構築および物体内部 の原子挙動の計算を行うための解析プログラムを構築し た．ここで，Mn 原子は，材料の温度によって結晶構造が複 雑に変化することが分かっている．本研究では，系の温度T が0.01K と極低温であることを仮定して-Mn に近い結晶構 造の構築が可能か検討を行った．-Mn の結晶構造は，単位 胞中に58 個の原子を含む体心立方格子であり，複雑な構造 を有する(16)．本研究においては，平衡状態における Mn の 各原子の配置が不明であったため，単純な BCC（体心立方 格子）構造をコンピュータ内で再現することが可能か検討 した．計算開始時における各原子の初期速度については， 乱数および系の温度を考慮して決定した．乱数については， Intel Math Kernel Library（数値計算・統計関数ライブラリー） に含まれるMersenne twister 疑似乱数列生成器（MT19937） により生成される乱数を用いた．なお，本研究で構築した 大規模MD 解析用プログラムにおいては，計算開始時に疑 似乱数を生成するための初期値を明示的に入力するように しており，この初期値を含めた解析条件が全て同じであれ ば，再計算を行った場合においても全く同じ結果が得られ るようにしている．計算開始直後において系は平衡状態で はないため，次式で示される系の運動エネルギーと温度の 関係を用いて各原子の速度を制御していき，平衡状態にな るようにしている． 2 3 2 2 N i i B i m v Nk T
... (6) ここで，N は系に存在する原子の個数，kBはボルツマン定数， T は系の温度，mi，viはそれぞれ原子i の質量および速度で ある．初期速度を与えた後，各原子に作用する力を式(3)， (5)により算出し，運動方程式から各原子の変位および速度 を更新する．運動方程式の数値積分は速度ベルレー法を用 いた(17)．速度ベルレー法は，Newton の運動方程式を満たす ように原子挙動を追跡できる数値積分法である．微小時間 をt としたとき，ある時刻 t からt だけ経過した後の原子 の位置と速度をテイラー展開により求めると，式(7)および (8)のようになる．鞆田 顕章
2 2 2 3 3 4 3 2! 3! i i i i i dr t t d r t r t t r t t dt dt t d r t O t dt ... (7)
2 2 2 3 3 4 3 2! 3! i i i i i dv t t d v t v t t v t t dt dt t d v t O t dt ... (8) 高次項を無視して整理すると，
2 2 2i
2 i i i t d r t r t t r t t v t dt ... (9)
2 2 2 2 2 2 i i i i d r t v t t v t t dt t d d r t dt dt ... (10) となる．ここで，原子i における運動方程式
2 2 i i i d r t f t m dt ... (11) を式(9)，(10) に代入すると，
2 2
i i i i i t r t t r t t v t f t m ... (12)
2 2
2
i
i i i i t df t v t t v t f t m t dt ... (13) となる．ここで，式(13)の右辺第 2 項に注目し，前進差分を 用いると，
i i i df t f t t f t dt t ... (14) となる．これを式(13)に代入すると，
2
i i i i i t v t t v t f t t f t m ... (15) となる．式(12)，(15)より，微小時間t 経過後の原子の位置 と速度を求めることが可能となる． MD 解析では，式(3)，(5)を計算する際に系内の原子全て を対象とすると，N (N - 1) 回の計算が必要となる．系内に 存在する原子の個数が多い場合，1 ステップ分の速度および 位置の更新に必要な計算回数が膨大となるため，解析結果 を得るまでに相当な時間を要することになる．ここで，金 属における結晶中の原子は，第二近接原子間距離よりも離 れた原子同士は互いにほとんど影響を及ぼさないことが分 かっている．この距離はカットオフ距離と呼ばれている． MD 解析においては，解析開始前に予めカットオフ距離を設 定し，カットオフ距離よりも離れた原子同士については相 互作用の計算を一切行わないようにすることで，計算の高 速化を図る．本研究では，複数のCPU および CPU 内の複数 のコアを同時に用いる並列計算を行うことを予定している ため，一辺がカットオフ距離の 2 倍程度の長さの立方体の ブロックを設定し，解析対象系を複数のブロックに分割し た．各ブロックには予め番号を設定しておき，番号の小さ いブロックから順番にブロック内の原子における相互作用 を計算した．上述の方法は領域分割法と呼ばれており，大 規模FEM，MD 等による数値解析モデルを構築する際に用 いられる．ただし，計算対象のブロック内に存在する原子 の一部は，隣接するブロックの原子の一部との距離がカッ トオフ距離以内であることが想定される．そこで，解析時 においては，計算対象ブロックとその周囲26 ブロック内に ある原子の間の距離を全て計算し，カットオフ距離以下の 原子ペアについては相互作用の計算を行うことにした．上 述の方法の導入により，解析時に必要なメインメモリの容 量を削減し，多数のCPU を有するスーパーコンピュータに おいて大規模並列計算を可能にした．なお，次節で示す解 析結果については，九州大学スーパーコンピュータ（Fujitsu PRIMERGY CX2550/CX2560 M4）を 1 ノード（2CPU 搭載， 18core/CPU）分利用し計算したものである． 〈3･3〉 解析結果 前節で述べた MD 解析用プログラムを用いて，単結晶の Mn 原子の挙動を計算した．数値解析モデルとして，Mn 原 子の平衡原子間距離R0 = 2.53Å および格子定数 2.92Å を考 慮し，平衡状態におけるBCC 構造の各原子の位置に Mn 原 子を配置した．図3 に示すように，解析対象の原子数 N = 128000，解析対象セルを一辺が約 11.7 nm の立方体とした． 本研究では，解析対象セルの外側に解析対象セルの原子の 情報をコピーすることで，計算領域外においても原子が無 数に存在しているという条件（周期境界条件）を適用した． これにより，解析対象セルの境界面付近に存在する原子に 作用する力を適切に計算することができる．数値積分にお ける時間刻みt は原子の振動周期を考慮し 1fs，総ステップ 数を10000（10ps），系の目標温度 TRは0.01K とし，系の原 子数N，体積 V を一定として温度が TRとなるように制御し ながら解析を進めた．解析結果を図4 に示す．解析開始後， 0.1ps まではポテンシャルエネルギーの値に変化がみられる が，その後すぐに一定値となり，系は安定した状態を維持 していた．さらに，各時刻における原子の位置に大きな変 化は見られず，解析中において計算結果が不安定になる等 の不具合は認められなかった．また，Mn 原子 128000 個の 挙動を10000 ステップ分計算した際に要した時間は約 17.4 時間であった．著者らが前年度に構築した小規模MD 解析 用プログラム（MATLAB 使用）では，Cu 原子 500 個におけ る10000 ステップ分の計算に要した時間は約 0.8 時間であっ た．解析対象の原子の個数を考慮した上で両者を比較する と，本研究で構築した大規模MD 解析用プログラム（Fortran） は，前年度までの研究で構築した小規模MD 解析用プログ ラム（MATLAB）よりも計算時間が約 1/12 程度となった． 大規模MD 解析用プログラムを用いることにより，系内に 原子が大量に存在する場合においても高速に計算を行うこ とが可能となった．本報告で紹介した大規模MD 解析用プ ログラムは単一の計算ノードのみを利用して解析を行う仕 様としたが，今後は，MPI（メッセージパッシングインター フェース）の利用により複数ノードを同時に使用すること が可能なMD 解析用プログラムを構築することを予定している． 以上より，本研究で作成した MD 解析用プログラムを用 いたMn 原子 128000 個の極低温状態における結晶構造の構 築が可能であることを示した． 図3 MD 解析モデル（Mn 原子，N = 128000） Potential Energy [J] 図4 系のポテンシャルエネルギー
4. 結言
本報では，金属材料の内部摩擦を考慮した振動シミュレ ータ開発のため，Mn 基制振合金板の振動減衰能取得および 分子動力学モデルの検討を行い，以下の結果を得た． (1) M2052 合金製薄板の振動減衰能取得試験を行い，1～10 次モードにおける各モードの固有振動数，減衰比を取得 することができた．各次の減衰比は概ね1.2～2.5％程度 であった． (2) Mn 原子を解析対象とし，分子動力学法を用いた大規模 解析プログラムを構築した．極低温時における-Mn の 原子挙動を解析するため，128000 個の Mn 原子を体心 立方格子中に配置した．その結果，原子の位置に大きな 変化は見られず，安定した構造となることが示された． ただし，今回の解析においては，実際の-Mn とは異な る結晶構造を構築した．今後は，-Mn の実際の原子配 置を調査するとともに，多元系合金の解析を行うための MD 解析用プログラムの構築を目指す．謝 辞
本研究は，福岡工業大学平成29 年度新任教員スタートア ップ支援の助成を受けた．また，本研究の一部は，JSPS 科 研費若手研究B（課題番号 JP16K18037）の助成を受けた． 本研究における数値解析は，九州大学情報基盤研究開発セ ンターの研究用計算機システムを利用して実施した．ここ に記して謝意を表す． （平成３０年９月２０日受付）文 献
(１) 大熊正明，長松昭男：「特性行列の実験的決定による系の同定（第 1 報，方法の提案と基礎的検討）」，日本機械学会論文集C 編，Vol.51， No.464，pp.719-728 (1985) (２) 中野陽介，吉村卓也：「モード特性を用いた減衰行列の同定」，日本 機械学会機械力学・計測制御講演論文集，pp.(323-1)-(323-9) (2013) (３) 稲村豊四郎，武澤伸浩，熊木保浩：「原子系/連続体変換とナノスケ ー ル 切 削 の 応 力 ひ ず み 解 析 」， 精 密 工 学 会 誌 編 ，Vol.60 ， No.4, pp.570-575 (1994) (４) 北村隆行，屋代如月，大谷隆一：「ニッケルナノ単結晶の引張変形と 破壊の原子シミュレーション」，日本機械学会論文集A 編，Vol.62， No.600，pp.1813-1819 (1996) (５) 中谷彰宏，北川浩，下川智嗣：「分子動力学法によるナノ多結晶体の 力学特性と引張荷重下における非弾性変形機構の検討」，日本機械学 会論文集A 編，Vol.66，No.643，pp. 435-441 (2000)(６) Hara, S., Kumagai, T., Izumi, S., Sakai, S.：“Multiscale analysis on the onset of nanoindentation-induced delamination: Effect of high-modulus Ru overlayer”, Acta MATERIALIA, Vol.57, pp.4209-4216 (2009)
(７) Tsuru, T., Aoyagi, Y., Kaji, Y., Shimokawa, T.：“Influence of Competition between Intergranular Dislocation Nucleation and Intergranular Slip Transfer on Mechanical Properties of Ultrafine-Grained Metals”, Materials Transactions, Vol.54, No.9, pp.1580-1586 (2013)
(８) 杉本孝一：「防振合金の最近の進歩」，鉄と鋼，Vol.60，pp.2203-2220 (1974) (９) 川原浩司，殷福星：「新しいタイプの振動吸収材料，M2052 制振合 金」，真空，Vol.42，No.1，pp.11-17 (1999) (10) 入江敏博：「機械振動学通論 第 2 版」，朝倉書店，pp.24-28 (2001) (11) 日本機械学会編：「コンピュータアナリシスシリーズ 7 原子・分子 モデルを用いる数値シミュレーション」，コロナ社，pp.12-21 (1996) (12) 上田顕：「分子シミュレーション －古典系から量子系手法まで」，裳 華房，pp.64-66 (2003)
(13) Baskes, M. I.：“Modified embedded-atom potentials for cubic materials and impurities”, Physical Review B, Vol.46, No.5, pp.2727-2742 (1992) (14) Kim, H.-K., Jung., W.-S., Lee, B.-J.：“Modified embedded-atom method
interatomic potentials for pure Mn and the Fe–Mn system”, Acta MATERIALIA, Vol.57, pp. 474-482 (2009)
(15) 鞆田顕章：「接触を伴う振動系の分子動力学解析」，日本機械学会 Dynamics and Design Conference 2018 講演論文集，No.158 (16) Bradley, A. J., Thewlis, J.：”The crystal structure of α-manganese”,
Proceedings of the Royal Society of London. Series A, Vol.115, Issue 771, pp.456-471 (1927)
(17) Andersen, H. C., Berens, P. H., Wilson, K. R.：“A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters”, The Journal of Chemical Physics, Vol.76, Issue 1, pp.637-649(1982)
