特表 2024-513620 仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法及びシステム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公表特許公報(A)

(11)【公表番号】

(43)【公表日】2024-03-27

(54)【発明の名称】仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法及びシステム

(51)【国際特許分類】

   G06F 30/23 20200101AFI20240319BHJP

   G01N 3/00 20060101ALI20240319BHJP

   G06F 119/04 20200101ALN20240319BHJP

【ＦＩ】

G06F30/23

G01N3/00 Z

G06F119:04

【審査請求】有

【予備審査請求】未請求

(21)【出願番号】P 2022577653

(86)(22)【出願日】2022-07-22

(85)【翻訳文提出日】2022-12-16

(86)【国際出願番号】 CN2022107269

(87)【国際公開番号】W WO2023151233

(87)【国際公開日】2023-08-17

(31)【優先権主張番号】202210120848.2

(32)【優先日】2022-02-09

(33)【優先権主張国・地域又は機関】CN

(81)【指定国・地域】

【公序良俗違反の表示】

（特許庁注：以下のものは登録商標）

１．ＭＡＴＬＡＢ

(71)【出願人】

【識別番号】519291342

【氏名又は名称】南京航空航天大学

(74)【代理人】

【識別番号】100095407

【弁理士】

【氏名又は名称】木村 満

(74)【代理人】

【識別番号】100132883

【弁理士】

【氏名又は名称】森川 泰司

(74)【代理人】

【識別番号】100148633

【弁理士】

【氏名又は名称】桜田 圭

(74)【代理人】

【識別番号】100147924

【弁理士】

【氏名又は名称】美恵 英樹

(72)【発明者】

【氏名】王 英玉

(72)【発明者】

【氏名】王 文軒

(72)【発明者】

【氏名】▲ゴン▼ 帥

(72)【発明者】

【氏名】張 暁凡

(72)【発明者】

【氏名】姚 衛星

【テーマコード（参考）】

2G061

5B146

【Ｆターム（参考）】

2G061AA08

2G061AB05

2G061BA15

2G061DA11

2G061DA12

2G061EC02

2G061EC04

5B146DJ02

5B146DJ07

5B146DJ15

(57)【要約】

本発明は、仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法及びシステムを開示する。方法は、以下のステップを含む。金属部材に対して有限要素分析を行い、危険点を確定し、危険点における応力歪み荷重時間の経過を求める。危険点を通過した各平面上のせん断歪みエネルギーを計算し、せん断歪みエネルギーが最大の平面を臨界面とし、臨界面上のせん断歪みエネルギー及び多軸応力比を得る。危険点における仮想歪みエネルギーを計算し、ねじり荷重のせん断歪みエネルギー寿命曲線を用いて寿命予測を行う。本発明における仮想歪みエネルギーは、単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線におけるせん断歪みエネルギーを正規化処理した後の歪みエネルギーである。本発明が提案した仮想歪みエネルギーに基づく方法は、同時に荷重経路が損傷パラメータ係数と指数に与える影響を考慮し、パラメータの確定が簡単であり、同時に低サイクル、高サイクルの異なる荷重経路下の疲労寿命分析を適用することができる。
【選択図】図１

【特許請求の範囲】

【請求項1】

以下のステップを含むことを特徴とする仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法。
（１）金属部材に対して有限要素分析を行い、危険点を確定し、危険点における応力歪み荷重時間の経過を求める。
（２）危険点を通過した各平面上のせん断歪みエネルギーを計算し、せん断歪みエネルギーが最大の平面を臨界面とし、臨界面上のせん断歪みエネルギー及び多軸応力比を得る。ここで、せん断歪みエネルギーは、せん断応力変化幅とせん断歪み変化幅の積として定義され、多軸応力比は、臨界面上の正応力最大値とせん断応力振幅の比として定義される。
（３）危険点における仮想歪みエネルギーを計算し、ねじり荷重のせん断歪みエネルギー寿命曲線を用いて寿命予測を行う。前記仮想歪みエネルギーは、単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線におけるせん断歪みエネルギーを正規化処理した後の歪みエネルギーである。

【請求項2】

前記ステップ（１）は、
所定の外荷重に対しては、３Ｄモデリングソフトウェアを使用して部材の幾何学モデルを描画し、メッシュ分割を行い、疲労危険点付近をメッシュ精密化することと、
材料属性を有限要素モデルに付与し、有限要素モデルに境界条件を追加し、実際の環境における部材の荷重状況を模擬することと、
有限要素解析により、危険点における応力歪み荷重時間の経過を得ることと、を含むことを特徴とする、
請求項１に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法。

【請求項3】

前記ステップ（２）における危険点を通過した各平面上のせん断歪みエネルギー計算方法は、以下を含むことを特徴とする、
請求項１に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法。
金属部材が時間［０，Ｔ］内に荷重を受けるとし、疲労危険点Ｏに局所座標系Ｏｘｙｚを構築し、危険点Ｏを通過する任意の平面Δは、その平面上の２つの方向ベクトルと１つの平面の法ベクトルで記述される。具体的には、ベクトルａは、平面Δと座標系Ｏｘｙｚにおける平面Ｏｘｙとの交線であり、ベクトルｂは、平面Δ上のａ軸に垂直なベクトルであり、ｎは、平面Δの法ベクトルである。ベクトルｎ，ａ，ｂと座標系Ｏｘｙｚとの位置関係は、３つの角度（φ，θ，α）で表される。ここで、φは、平面Δの方向法ベクトルｎとｘ軸の角度であり、θは、平面Δの方向法ベクトルｎとｚ軸の角度であり、αは、平面Δ上の任意のベクトルｑとベクトルａの角度である。危険点Ｏにおける応力歪み荷重時間の経過は、次の行列で表される。

【数1】

式中、σ_ｘ（ｔ），σ_ｙ（ｔ），σ_ｚ（ｔ）は、正応力成分であり、τ_ｘｙ（ｔ），τ_ｙｚ（ｔ），τ_ｘｚ（ｔ）は、せん断応力成分であり、ε_ｘ（ｔ），ε_ｙ（ｔ），ε_ｚ（ｔ）は、正歪み成分であり、γ_ｘｙ（ｔ），γ_ｙｚ（ｔ），γ_ｘｚ（ｔ）は、せん断歪み成分であり、ｔは、［０，Ｔ］区間内のいずれかの時刻である。
平面Δを記述する３つの単位ベクトルｎ、ａ、ｂは、座標変換によって次のように表される。

【数2】

平面Δ上の任意の単位ベクトルｑは、次のように表される。

【数3】

任意角度の正応力σ_ｎ（ｔ）、せん断応力τ_ｑ（ｔ）、せん断歪みγ_ｑ（ｔ）は、次のように表される。

【数4】

任意平面上の正応力平均値σ_ｎ，ｍを［０，Ｔ］区間における正応力積分の平均値に定義し、正応力振幅σ_ｎ，αは、正応力の分散によって定義される。

【数5】

式中、

【数6】

せん断応力、せん断歪みの平均値及び振幅τ_ｍ，τ_α，γ_ｍ，γ_αを定義する。

【数7】

式中、

【数8】

せん断歪みエネルギーΔτΔγは、部材の安定サイクル下でヒステリシスループを囲む最小矩形、すなわちせん断応力変化幅Δτとせん断歪み変化幅Δγとの積として定義し、変化幅Δτ＝２τ_α、Δγ＝２γ_α、そこで、ΔτΔγ＝４τ_α×γ_α。

【請求項4】

多軸応力比ρは、次のように表されることを特徴とする、
請求項３に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法。

【数9】

ここで、σ_{ｎ，ｍａｘ}は、臨界面における正応力の最大値であり、σ_ｎ，ｍは、正応力平均値であり、σ_ｎ，αは、正応力振幅であり、τ_αは、せん断応力振幅である。

【請求項5】

前記ステップ（３）における単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線は、次の方法に基づいて決定されることを特徴とする、
請求項１に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法。
単軸引張圧縮、ねじり荷重下の各疲労データ点臨界面上のせん断歪みエネルギーをそれぞれ計算し、首尾連係してせん断歪みエネルギー－寿命曲線を描画し、べき関数で２種類の荷重状況下でせん断歪みエネルギー－寿命曲線の各パラメータＡ_１，Ｂ_１，Ａ_３，Ｂ_３をそれぞれフィッティングして解く。
２つの荷重におけるせん断歪みエネルギー寿命曲線関数は、次のように表されることを特徴とする。

【数10】

ここで、ΔτΔγは、せん断歪みエネルギーであり、Ｎ_ｆは、寿命である。

【請求項6】

前記ステップ（３）において正規化された仮想歪みエネルギーＷ_Ｖは、次のように表されることを特徴とする、
請求項５に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法。

【数11】

ここで、

【数12】

は、正規化係数であり、ρは、多軸応力比である。

【請求項7】

金属部材の幾何学モデル及び外荷重を入力するための入力モジュールと、
処理モジュールと、
所定の外荷重における金属部材の予測寿命を表示するための出力モジュールを含み、
前記処理モジュールは、
金属部材に対して有限要素分析を行い、危険点を確定し、危険点における応力歪み荷重時間の経過を求めるための応力歪み計算ユニットと、
危険点を通過した各平面上のせん断歪みエネルギーを計算し、せん断歪みエネルギーが最大の平面を臨界面とし、臨界面上のせん断歪みエネルギー及び多軸応力比を得るための臨界面計算ユニットであって、ここで、せん断歪みエネルギーは、せん断応力変化幅とせん断歪み変化幅の積として定義され、多軸応力比は、臨界面上の正応力最大値とせん断応力振幅の比として定義される臨界面計算ユニットと、
危険点における仮想歪みエネルギーを計算し、ねじり荷重のせん断歪みエネルギー寿命曲線を用いて寿命予測を行うための予測ユニットであって、ここで、前記仮想歪みエネルギーは、単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線におけるせん断歪みエネルギーを正規化処理した後の歪みエネルギーである予測ユニットと、を含むことを特徴とする、
仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測システム。

【請求項8】

前記単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線は、べき関数フィッティングユニットによって生成され、
前記べき関数フィッティングユニットは、単軸引張圧縮、ねじり荷重下の各疲労データ点臨界面上のせん断歪みエネルギーをそれぞれ計算し、首尾連係してせん断歪みエネルギー－寿命曲線を描画し、べき関数で２種類の荷重状況下でせん断歪みエネルギー－寿命曲線の各パラメータＡ_１，Ｂ_１，Ａ_３，Ｂ_３をそれぞれフィッティングして解き、
２つの荷重におけるせん断歪みエネルギー寿命曲線関数は、次のように表されることを特徴とする、
請求項７に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測システム。

【数13】

ここで、ΔτΔγは、せん断歪みエネルギーであり、Ｎ_ｆは、寿命である。

【請求項9】

前記予測ユニットで正規化された仮想歪みエネルギーＷ_Ｖは、次のように表されることを特徴とする、
請求項８に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測システム。

【数14】

ここで、

【数15】

は、正規化係数であり、ρは、多軸応力比である。

【請求項10】

メモリ、プロセッサ、及びメモリ上に記憶されプロセッサ上で動作可能なコンピュータプログラムを含むコンピュータシステムであって、
前記コンピュータプログラムがプロセッサにロードされると、請求項１～６のいずれか１項に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法を実現することを特徴とする、
コンピュータシステム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法及びシステムに関し、増材製造金属材料の寿命予測分野に属し、航空システムの技術分野に用いることができる。

【背景技術】

【0002】

航空機部品は、服役中に循環荷重を受ける。疲労損傷は、航空機構造の故障の主な原因の一つである。航空会社は、運営コストを節約するために、軽量化と長寿命化の設計案に対する要求がますます切実になっている。増材製造（ＡＭ）技術は、複雑な構造のニアネットシェイプに適用され、構造トポロジー最適化の要求を満たすことができるとともに、快速成形と材料利用率が高い利点があり、航空工業の広範な注目を集めている。

【0003】

研究により、増材製造加工の部品の静力性能は、通常の加工方式で加工された部品と大差ないが、疲労性能の差が大きく、疲労寿命が低く、疲労分散性が大きいことが示された。そのため、増材製造技術は、航空機の主支持力構造の製造には応用されていない。増材製造加工技術の特徴により、加工過程において大量の欠陥が避けられず、疲労亀裂の発生に良い環境を提供した。現在、科学界は、増材製造における疲労挙動に対して大量の研究を行い、ミクロとマクロの２つの研究方向に分けることができる。例えばミクロの角度から分析すると、ＣＴスキャン、顕微鏡撮影などの手段を通じて、欠陥、表面粗さ、亀裂拡張速度などの方面を通じて分析を行うことが多い。マクロの分析方法は、増材製造材料に対する応力寿命曲線を確立するか、大域パラメータの代わりに穴付近の局所パラメータを採用するなどの方法に分けられる。ミクロ分析は、疲労破壊の原理を深く理解するのに役立ち、マクロ分析は、工学応用においてより便利である。総じて言えば、科学界は、増材製造の研究方向に対して、主に単軸状況の予測分析に集中しており、多軸荷重下の疲労問題がまだよく解決されていない。

【0004】

飛行機の構造は、服役中に複雑な多軸荷重を受ける。例えば、翼全体は、せん断力、曲げモーメント、トルクの作用を受ける。翼のスキンは、正応力とせん断応力を同時に受ける。飛行機の着陸機構造は、離陸と着陸の過程で複数の方向からの荷重を受ける。多軸疲労故障は、通常、これらの構造が故障する原因である。現在、多軸疲労寿命の予測は、主に多軸応力状態を単軸応力状態に等価化し、更に単軸応力寿命曲線に基づいて疲労寿命を予測する。これらの方法は、応力基準、歪み基準、エネルギー基準に分けることができる。エネルギー基準によると、疲労危険点付近のエネルギーの蓄積が疲労破壊の原因であり、疲労破壊の内在メカニズムに比較的符合する。同時にエネルギーは、スカラーであり、計算過程も簡単である。しかし、各種の方法は、特定の材料又は荷重条件に適用され、現在まだ公認された工程構造設計要求に符合する増材製造金属材料の多軸疲労寿命予測方法がない。

【発明の概要】

【発明が解決しようとする課題】

【0005】

本発明の目的は、上記従来技術の不足に対して、多軸仮想歪みエネルギー処理方法を提供し、増材製造金属材料の多軸疲労寿命予測問題を解決することにある。

【課題を解決するための手段】

【0006】

上述した目的を達成するために、本発明は、以下の技術手段を採用する。

【0007】

仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法であって、以下のステップを含む。
（１）金属部材に対して有限要素分析を行い、危険点を確定し、危険点における応力歪み荷重時間の経過を求める。
（２）危険点を通過した各平面上のせん断歪みエネルギーを計算し、せん断歪みエネルギーが最大の平面を臨界面とし、臨界面上のせん断歪みエネルギー及び多軸応力比を得る。ここで、せん断歪みエネルギーは、せん断応力変化幅とせん断歪み変化幅の積として定義され、多軸応力比は、臨界面上の正応力最大値とせん断応力振幅の比として定義される。
（３）危険点における仮想歪みエネルギーを計算し、ねじり荷重のせん断歪みエネルギー寿命曲線を用いて寿命予測を行う。前記仮想歪みエネルギーは、単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線におけるせん断歪みエネルギーを正規化処理した後の歪みエネルギーである。

【0008】

好ましくは、前記ステップ（１）は、所定の外荷重に対しては、３Ｄモデリングソフトウェアを使用して部材の幾何学モデルを描画し、メッシュ分割を行い、疲労危険点付近をメッシュ精密化することと、材料属性を有限要素モデルに付与し、有限要素モデルに境界条件を追加し、実際の環境における部材の荷重状況を模擬することと、有限要素解析により、危険点における応力歪み荷重時間の経過を得ることとを含む。

【0009】

好ましくは、前記ステップ（２）における危険点を通過した各平面上のせん断歪みエネルギー計算方法は、以下を含む。
金属部材が時間［０，Ｔ］内に荷重を受けるとし、疲労危険点Ｏに局所デカルト座標系Ｏｘｙｚを構築し、危険点Ｏを通過する任意の平面Δは、その平面上の２つの方向ベクトルと１つの平面の法ベクトルで記述される。具体的には、ベクトルａは、平面Δと座標系Ｏｘｙｚにおける平面Ｏｘｙとの交線であり、ベクトルｂは、平面Δ上のａ軸に垂直なベクトルであり、ｎは、平面Δの法ベクトルである。ベクトルｎ，ａ，ｂと座標系Ｏｘｙｚとの位置関係は、３つの角度（φ，θ，α）で表される。ここで、φは、平面Δの方向法ベクトルｎとｘ軸の角度であり、θは、平面Δの方向法ベクトルｎとｚ軸の角度であり、αは、平面Δ上の任意のベクトルｑとベクトルａの角度である。危険点Ｏにおける応力歪み荷重時間の経過は、次の行列で表される。

【0010】

【数1】

【0011】

式中、σ_ｘ（ｔ），σ_ｙ（ｔ），σ_ｚ（ｔ）は、正応力成分であり、τ_ｘｙ（ｔ），τ_ｙｚ（ｔ），τ_ｘｚ（ｔ）は、せん断応力成分であり、ε_ｘ（ｔ），ε_ｙ（ｔ），ε_ｚ（ｔ）は、正歪み成分であり、γ_ｘｙ（ｔ），γ_ｙｚ（ｔ），γ_ｘｚ（ｔ）は、せん断歪み成分であり、ｔは、［０，Ｔ］区間内のいずれかの時刻である。
平面Δを記述する３つの単位ベクトルｎ、ａ、ｂは、座標変換によって次のように表される。

【0012】

【数2】

【0013】

平面Δ上の任意の単位ベクトルｑは、次のように表される。

【0014】

【数3】

【0015】

任意角度の正応力σ_ｎ（ｔ）、せん断応力τ_ｑ（ｔ）、せん断歪みγ_ｑ（ｔ）は、次のように表される。

【0016】

【数4】

【0017】

任意平面上の正応力平均値σ_ｎ，ｍを［０，Ｔ］区間における正応力の積分の平均値に定義し、正応力振幅σ_ｎ，αは、正応力の分散によって定義される。

【0018】

【数5】

式中、

【数6】

【0019】

同様に、せん断応力、せん断歪みの平均値及び振幅τ_ｍ，τ_α，γ_ｍ，γ_αを定義する。

【0020】

【数7】

式中、

【数8】

【0021】

せん断歪みエネルギーΔτΔγは、部材の安定サイクル下でヒステリシスループを囲む最小矩形、すなわちせん断応力変化幅Δτとせん断歪み変化幅Δγとの積として定義し、変化幅Δτ＝２τ_α、Δγ＝２γ_α、そこで、ΔτΔγ＝４τ_α×γ_α。多軸応力比ρは、正応力の最大値σ_{ｎ，ｍａｘ}とせん断応力振幅τ_αとの比として定義される。

【0022】

【数9】

【0023】

好ましくは、前記ステップ（３）における単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線は、次の方法に基づいて決定される。
単軸引張圧縮、ねじり荷重下の各疲労データ点臨界面上のせん断歪みエネルギーをそれぞれ計算し、首尾連係してせん断歪みエネルギー－寿命曲線を描画し、べき関数で２種類の荷重状況下でせん断歪みエネルギー－寿命曲線の各パラメータＡ_１，Ｂ_１，Ａ_３，Ｂ_３をそれぞれフィッティングして解く。

【0024】

【数10】

【0025】

ここで、ΔτΔγは、せん断歪みエネルギーであり、Ｎ_ｆは、寿命である。

【0026】

好ましくは、前記ステップ（３）において正規化された仮想歪みエネルギーＷ_Ｖは、次のように表される。

【0027】

【数11】

ここで、

【数12】

は、正規化係数であり、ρは、多軸応力比である。

【0028】

同じ発明思想に基づいて、本発明による仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測システムは、金属部材の幾何学モデル及び外荷重を入力するための入力モジュールと、処理モジュールと、所定の外荷重における金属部材の予測寿命を表示するための出力モジュールを含む。前記処理モジュールは、金属部材に対して有限要素分析を行い、危険点を確定し、危険点における応力歪み荷重時間の経過を求めるための応力歪み計算ユニットと、危険点を通過した各平面上のせん断歪みエネルギーを計算し、せん断歪みエネルギーが最大の平面を臨界面とし、臨界面上のせん断歪みエネルギー及び多軸応力比を得るための臨界面計算ユニットであって、ここで、せん断歪みエネルギーは、せん断応力変化幅とせん断歪み変化幅の積として定義され、多軸応力比は、臨界面上の正応力最大値とせん断応力振幅の比として定義される臨界面計算ユニットと、危険点における仮想歪みエネルギーを計算し、ねじり荷重のせん断歪みエネルギー寿命曲線を用いて寿命予測を行うための予測ユニットであって、ここで、前記仮想歪みエネルギーは、単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線におけるせん断歪みエネルギーを正規化処理した後の歪みエネルギーである予測ユニットとを含む。

【0029】

好ましくは、前記単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線は、べき関数フィッティングユニットによって生成され、前記べき関数フィッティングユニットは、単軸引張圧縮、ねじり荷重下の各疲労データ点臨界面上のせん断歪みエネルギーをそれぞれ計算し、首尾連係してせん断歪みエネルギー－寿命曲線を描画し、べき関数で２種類の荷重状況下でせん断歪みエネルギー－寿命曲線の各パラメータＡ_１，Ｂ_１，Ａ_３，Ｂ_３をそれぞれフィッティングして解く。

【0030】

好ましくは、前記予測ユニットで正規化された仮想歪みエネルギーＷ_Ｖは、次のように表される。

【0031】

【数13】

【0032】

同じ発明思想に基づいて、本発明によるコンピュータシステムは、メモリ、プロセッサ、及びメモリ上に記憶されプロセッサ上で動作可能なコンピュータプログラムを含み、前記コンピュータプログラムがプロセッサにロードされると、上記の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法を実現する。

【発明の効果】

【0033】

本発明は、典型的な靭性金属材料の単軸引張圧縮、ねじり状態におけるせん断歪みエネルギー寿命曲線を分析し、２種類の荷重状況におけるせん断歪みエネルギーと寿命とは比較的に良い二対数線形関係を持つことを発見し、以下の２点に基づいてせん断歪みエネルギーを正規化処理して仮想歪みエネルギーを得、ねじり荷重のせん断歪みエネルギー寿命曲線を用いて寿命予測を行うことによって、異なる荷重経路の寿命分析問題に用いることができる。
１）靭性金属材料の疲労亀裂挙動がせん断破壊として現れ、ねじり荷重のせん断歪みエネルギー寿命曲線を用いたほうががより直感的である。
２）ねじり荷重のせん断歪みエネルギー寿命曲線は、すべての曲線の一番上に位置することが多い。この曲線の下にフィットするのは、より直感的であり、正規化係数形式が比較的簡単である。
また、従来の応力基準は、塑性歪みの小さい高サイクル疲労に用いられることが多く、歪み基準は、塑性歪みの大きい低サイクル疲労に用いられることが多い。本発明は、応力と歪みを掛け合わせたエネルギー基準を用いて低サイクル疲労、高サイクル疲労を両立することができる。

【図面の簡単な説明】

【0034】

【図1】本発明の実施例の全体フローチャートである。

【図2】本発明の実施例に係る危険点における局所座標系の定義概略図である。

【図3】本発明の実施例における疲労寿命の詳細な計算フローチャートである。

【図4】本発明の実験例の３１６Ｌ鋼の単軸引張圧縮、ねじり荷重下のせん断歪みエネルギーの寿命変化曲線図である。

【図5】本発明の実験例の有限要素モデリング・メッシュ分割の概略図である。

【図6】本発明の実験例の予測結果と試験結果の比較図である。

【図7】本発明の実施例のモジュール構成図である。

【発明を実施するための形態】

【0035】

当業者に分かりやすくするために、以下に実施例と図面を参照して本発明を更に説明するが、実施形態で言及される内容は、本発明を限定するものではない。

【0036】

図１に示すように、本発明の実施例によって開示される仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法は、主なステップが以下である。

【0037】

（１）金属部材に対して有限要素分析を行い、危険点を確定し、危険点における応力歪み荷重時間の経過を求める。主に以下を含む。
所定の外荷重に対して、３Ｄモデリングソフトウェアを用いて部材の幾何学モデルを構築し、有限要素ソフトウェアを導入する。
幾何学モデルをメッシュ分割し、計算結果が収束するまで疲労危険点付近をメッシュ精密化する。
材料属性を有限要素モデルに付与し、材料の疲労試験過程における塑性歪みの有無に基づいて、線弾性又は弾塑性の構成モデルを選択する。
有限要素モデルに境界条件を追加し、実際の環境における部材の荷重状況を模擬する。
最後、有限要素解析により、危険点Ｏにおける応力歪み荷重時間の経過を得て、後続の計算に備える。

【0038】

（２）危険点を通過した各平面上のせん断歪みエネルギーを計算し、せん断歪みエネルギーが最大の平面を臨界面とし、臨界面上のせん断歪みエネルギー及び多軸応力比を得る。主に以下を含む。
金属部材が時間［０，Ｔ］内に荷重を受けるとし、疲労危険点Ｏに局所デカルト座標系Ｏｘｙｚを構築し、危険点Ｏを通過する任意の平面Δは、その平面上の２つの方向ベクトルと１つの平面の法ベクトルで記述される。具体的には、ベクトルａは、平面Δと座標系Ｏｘｙｚにおける平面Ｏｘｙとの交線であり、ベクトルｂは、平面Δ上のａ軸に垂直なベクトルであり、ｎは、平面Δの法ベクトルである。ベクトルｎ，ａ，ｂと座標系Ｏｘｙｚとの位置関係は、３つの角度（φ，θ，α）で表される。ここで、φは、平面Δの方向法ベクトルｎとｘ軸の角度であり、θは、平面Δの方向法ベクトルｎとｚ軸の角度であり、αは、平面Δ上の任意のベクトルｑとベクトルａの角度である。危険点Ｏにおける応力歪み荷重時間の経過は、次の行列で表される。

【0039】

【数14】

【0040】

式中σ_ｘ（ｔ），σ_ｙ（ｔ），σ_ｚ（ｔ）は、正応力成分であり、τ_ｘｙ（ｔ），τ_ｙｚ（ｔ），τ_ｘｚ（ｔ）は、せん断応力成分であり、ε_ｘ（ｔ），ε_ｙ（ｔ），ε_ｚ（ｔ）は、正歪み成分であり、γ_ｘｙ（ｔ），γ_ｙｚ（ｔ），γ_ｘｚ（ｔ）は、せん断歪み成分であり、ｔは、［０，Ｔ］区間内のいずれかの時刻である。
平面Δを記述する３つの単位ベクトルｎ、ａ、ｂは、座標変換によって次のように表される。

【0041】

【数15】

【0042】

平面Δ上の任意の単位ベクトルｑは、次のように表される。

【0043】

【数16】

【0044】

任意角度の正応力σ_ｎ（ｔ）、せん断応力τ_ｑ（ｔ）、せん断歪みγ_ｑ（ｔ）は、次のように表される。

【0045】

【数17】

【0046】

任意平面上の正応力平均値σ_ｎ，ｍを［０，Ｔ］区間における正応力の積分の平均値に定義し、正応力振幅σ_ｎ，αは、正応力の分散によって定義される。

【0047】

【数18】

式中、

【数19】

【0048】

同様に、せん断応力、せん断歪みの平均値及び振幅τ_ｍ，τ_α，γ_ｍ，γ_αを定義する。

【0049】

【数20】

式中、

【数21】

【0050】

せん断歪みエネルギーΔτΔγは、部材の安定サイクル下でヒステリシスループを囲む最小矩形、すなわちせん断応力変化幅Δτとせん断歪み変化幅Δγとの積として定義し、変化幅Δτ＝２τ_α、Δγ＝２γ_α、そこで、ΔτΔγ＝４τ_α×γ_α。

【0051】

多軸応力比ρは、正応力最大値σ_{ｎ，ｍａｘ}とせん断応力振幅τ_αの比として定義される。

【0052】

【数22】

【0053】

（３）危険点における仮想歪みエネルギーを計算し、ねじり荷重のせん断歪みエネルギー寿命曲線を用いて寿命予測を行う。ここでの仮想歪みエネルギーは、単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線におけるせん断歪みエネルギーを正規化処理した後の歪みエネルギーである。単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線は、予め、材料の単軸引張圧縮、ねじり疲労試験データに基づいて決定される。具体的には、単軸引張圧縮、ねじり荷重下の各疲労データ点臨界面上のせん断歪みエネルギーをそれぞれ計算し、首尾連係してせん断歪みエネルギー－寿命曲線を描画し、２種類の荷重状況下でせん断歪みエネルギー－寿命曲線の各パラメータＡ_１，Ｂ_１，Ａ_３，Ｂ_３をそれぞれフィッティングして解く。
典型的な靭性金属材料の単軸引張圧縮、ねじり状態におけるせん断歪みエネルギー?寿命曲線を分析した。２つの荷重状況におけるせん断歪みエネルギーと寿命とは、比較的に良い二対数線形関係を持つことを発見し、それぞれ２つの曲線にべき関数フィッティングを行うことで、寿命に対応するせん断歪みエネルギーを記述するために使用することができる。

【0054】

【数23】

【0055】

仮想歪みエネルギーは、異なる荷重経路のせん断歪みエネルギーを正規化することにより、疲労寿命を計算する。

【0056】

【数24】

【0057】

式中、多軸応力比ρは、異なる荷重経路の影響を反映する。

【0058】

仮想歪みエネルギーのねじり荷重下の寿命曲線に基づいて疲労寿命を計算する詳細な過程を図３に示す。

【0059】

以下、具体的な実験例を用いて本発明の提案を更に検証する。本例で用いた材料は、領域選択的レーザー溶融により生産された３１６Ｌステンレス鋼であり、多種類の荷重経路の疲労試験、すなわち単軸引張圧縮、ねじり、比例、９０°非比例多軸荷重を行った。まず試験片に対して有限要素分析を行い、危険点を確定し、危険点における応力歪み荷重時間の経過［σ（ｔ）］、［ε（ｔ）］を求めた。そして座標変換により、臨界面におけるせん断歪みエネルギーΔτΔγと多軸応力比ρを解く。材料の単軸引張圧縮、ねじり疲労試験データに基づいてモデルのパラメータＡ_１，Ｂ_１，Ａ_３，Ｂ_３を確定する。最後に、求められる荷重下の危険点における仮想歪みエネルギーＷ_Ｖを求め、寿命Ｎ_ｆを解く。詳細な手順は次のとおりです。

【0060】

Ｓ１：本例で用いた試験片材料は、増材製造３１６Ｌ鋼であり、ゲージ長が３８ｍｍ、直径が１２ｍｍの中実丸棒である。ＭＴＳ８０９引張ねじり試験機を用いて、平均応力を含まない単軸引張圧縮、ねじり、多軸比例荷重、９０°非比例荷重及び平均応力を含む単軸引張圧縮、多軸比例荷重、９０°非比例荷重試験を行った。単軸引張圧縮、ねじり試験片の疲労データは、本発明のモデルのパラメータを決定するために用いられ、他の荷重経路の疲労データは、本発明のモデルの正確性を検証するために用いられる。具体的な試験データは、表１に示されている。ここで、荷重比Ｒは、応力最小値と最大値の比であり、平均応力の大きさを表す。位相角Φは、正応力とせん断応力の位相差であり、非比例の程度を説明するために使用される。

【0061】

【表1】

【0062】

本例では、有限要素ソフトウェアＰａｔｒａｎ＆Ｎａｓｔｒａｎを使用して有限要素解析を行う。試験片の重心を原点とし、試験片の軸線方向をｘ軸とし、全体座標系を構築する。平滑化試験片のゲージ表面の各点の応力歪み状態は同じであるため、ゲージ表面の任意の位置から亀裂が発生する可能性があり、疲労危険点を座標（０、６、０）の点と仮定してもよい。

【0063】

メッシュ分割は、まず２次元シェルユニットを確立し、軸線に沿って回転スイープして３次元エンティティユニットを生成する方法を採用し、疲労危険点を局所的に暗号化し、１ｍｍに２０ユニットを分割し、有限要素モデリング結果を図５に示す。

【0064】

疲労試験の過程で、塑性歪みが発見されず、又は、塑性歪みが無視されてもよい。したがって、線形弾性構成モデルを採用し、弾性係数Ｅ＝１９０．８ＧＰａ、ポアソン比ｖ＝０．３は、静力試験によって測定された。

【0065】

試験機における試験片の荷重状況を模擬するために、試験片の一端を固定し、他端にＭＰＣを構築し、荷重を加えた。

【0066】

最後に有限要素計算を行い、応力歪み荷重時間の経過［σ（ｔ）］、［ε（ｔ）］を導出する。

【0067】

Ｓ２：ｍａｔｌａｂソフトウェアを用いてプログラミングを行い、３つの角度（φ，θ，α）に対して０°～１８０°の範囲内にステップサイズ１で各平面Δのせん断歪みエネルギーを求め、その中でせん断歪みエネルギーが最も大きい平面を探し出して臨界面とし、そして本モデルが必要とするパラメータを解いて、表２の３－７列に列挙した。

【0068】

【表2】

【0069】

Ｓ３：単軸引張圧縮（番号Ｐ－０３～Ｐ－１２）、ねじり（番号Ｐ－２１～Ｐ－２８）の場合のせん断歪みエネルギーΔτΔγ－寿命Ｎ_ｆ曲線を描画し、図４に示すように、ｍａｔｌａｂフィッティングプログラム又はＥｘｃｅｌフィッティング機能を用いて曲線に対してべき関数フィッティングを行い、Ａ_１＝５．４２、Ｂ_１＝－０．１３、Ａ_３＝３．０３、Ｂ_３＝－０．０６を求めた。

【0070】

Ｓ４：モデル：

【数25】

によって、各荷重状況における仮想歪みエネルギー値Ｗ_Ｖを計算し、寿命Ｎ_ｆ，ｅを計算し、それぞれ表２の８、９列目に列挙した。結果を比較しやすいように、計算結果と試験結果を図６に示す。図から分かるように、７９％の寿命予測点は、２倍誤差帯内にあり、９３％の寿命予測点は、３倍誤差帯内にあり、各種荷重経路の疲労寿命予測に対して良い予測効果を得た。

【0071】

同じ発明思想に基づいて、本発明の実施例による仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測システムは、図７に示すように、金属部材の幾何学モデル及び外荷重を入力するための入力モジュールと、処理モジュールと、所定の外荷重における金属部材の予測寿命を表示するための出力モジュールを含み、前記処理モジュールは、金属部材に対して有限要素分析を行い、危険点を確定し、危険点における応力歪み荷重時間の経過を求めるための応力歪み計算ユニットと、危険点を通過した各平面上のせん断歪みエネルギーを計算し、せん断歪みエネルギーが最大の平面を臨界面とし、臨界面上のせん断歪みエネルギー及び多軸応力比を得るための臨界面計算ユニットであって、ここで、せん断歪みエネルギーは、せん断応力変化幅とせん断歪み変化幅の積として定義され、多軸応力比は、臨界面上の正応力最大値とせん断応力振幅の比として定義される臨界面計算ユニットと、危険点における仮想歪みエネルギーを計算し、ねじり荷重のせん断歪みエネルギー寿命曲線を用いて寿命予測を行うための予測ユニットであって、ここで、前記仮想歪みエネルギーは、単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線におけるせん断歪みエネルギーを正規化処理した後の歪みエネルギーである予測ユニットと、を含む。単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線は、べき関数フィッティングユニットによって生成される。前記べき関数フィッティングユニットは、単軸引張圧縮、ねじり荷重下の各疲労データ点臨界面上のせん断歪みエネルギーをそれぞれ計算し、首尾連係してせん断歪みエネルギー－寿命曲線を描画し、べき関数で２種類の荷重状況下でせん断歪みエネルギー－寿命曲線の各パラメータをそれぞれフィッティングして解く。

【0072】

当業者であれば、上述した各ユニットの具体的な動作手順は、説明の容易さと簡潔さのために、上述した方法の実施例における対応する手順を参照することができ、ここではこれ以上説明しない。前記ユニットの分割は、１つの論理機能分割にすぎず、実際に実装される場合には、複数のユニットが結合されてもよく、又は別のシステムに統合されてもよいなど、別の分割方法があってもよい。

【0073】

同じ発明思想に基づいて、本発明の実施例によるコンピュータシステムは、メモリ、プロセッサ、及びメモリに記憶されプロセッサ上で実行可能なコンピュータプログラムを含み、前記コンピュータプログラムがプロセッサにロードされると、前記仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法を実現する。

【0074】

本発明の技術的態様の本質的又は先行技術に寄与する部分は、ソフトウェア製品の形態で具現化することができ、これは当業者にとって明らかである。該コンピュータソフトウェア製品は、本発明の実施例に記載された方法のすべて又は一部のステップを１台のコンピュータ装置（パーソナルコンピュータ、サーバ、又はネットワーク装置など）に実行させるためのいくつかの命令を含む１つの記憶媒体に記憶される。記憶媒体は、Ｕディスク、リムーバブルハードディスク、読み取り専用メモリＲＯＭ、ランダムアクセスメモリＲＡＭ、磁気ディスク、又は光ディスクなど、コンピュータプログラムを記憶することができる様々な媒体を含む。

【0075】

なお、上記は、本発明の好適な実施形態にすぎず、当業者にとっては、本発明の原理を逸脱することなく、いくつかの改良を行うことができ、これらの改良も本発明の保護範囲と見なすべきである。

【0076】

（付記）
（付記１）
以下のステップを含むことを特徴とする仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法。
（１）金属部材に対して有限要素分析を行い、危険点を確定し、危険点における応力歪み荷重時間の経過を求める。
（２）危険点を通過した各平面上のせん断歪みエネルギーを計算し、せん断歪みエネルギーが最大の平面を臨界面とし、臨界面上のせん断歪みエネルギー及び多軸応力比を得る。ここで、せん断歪みエネルギーは、せん断応力変化幅とせん断歪み変化幅の積として定義され、多軸応力比は、臨界面上の正応力最大値とせん断応力振幅の比として定義される。
（３）危険点における仮想歪みエネルギーを計算し、ねじり荷重のせん断歪みエネルギー寿命曲線を用いて寿命予測を行う。前記仮想歪みエネルギーは、単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線におけるせん断歪みエネルギーを正規化処理した後の歪みエネルギーである。

【0077】

（付記２）
前記ステップ（１）は、
所定の外荷重に対しては、３Ｄモデリングソフトウェアを使用して部材の幾何学モデルを描画し、メッシュ分割を行い、疲労危険点付近をメッシュ精密化することと、
材料属性を有限要素モデルに付与し、有限要素モデルに境界条件を追加し、実際の環境における部材の荷重状況を模擬することと、
有限要素解析により、危険点における応力歪み荷重時間の経過を得ることと、を含むことを特徴とする、
付記１に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法。

【0078】

（付記３）
前記ステップ（２）における危険点を通過した各平面上のせん断歪みエネルギー計算方法は、以下を含むことを特徴とする、
付記１に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法。
金属部材が時間［０，Ｔ］内に荷重を受けるとし、疲労危険点Ｏに局所座標系Ｏｘｙｚを構築し、危険点Ｏを通過する任意の平面Δは、その平面上の２つの方向ベクトルと１つの平面の法ベクトルで記述される。具体的には、ベクトルａは、平面Δと座標系Ｏｘｙｚにおける平面Ｏｘｙとの交線であり、ベクトルｂは、平面Δ上のａ軸に垂直なベクトルであり、ｎは、平面Δの法ベクトルである。ベクトルｎ，ａ，ｂと座標系Ｏｘｙｚとの位置関係は、３つの角度（φ，θ，α）で表される。ここで、φは、平面Δの方向法ベクトルｎとｘ軸の角度であり、θは、平面Δの方向法ベクトルｎとｚ軸の角度であり、αは、平面Δ上の任意のベクトルｑとベクトルａの角度である。危険点Ｏにおける応力歪み荷重時間の経過は、次の行列で表される。

【数26】

式中、σ_ｘ（ｔ），σ_ｙ（ｔ），σ_ｚ（ｔ）は、正応力成分であり、τ_ｘｙ（ｔ），τ_ｙｚ（ｔ），τ_ｘｚ（ｔ）は、せん断応力成分であり、ε_ｘ（ｔ），ε_ｙ（ｔ），ε_ｚ（ｔ）は、正歪み成分であり、γ_ｘｙ（ｔ），γ_ｙｚ（ｔ），γ_ｘｚ（ｔ）は、せん断歪み成分であり、ｔは、［０，Ｔ］区間内のいずれかの時刻である。
平面Δを記述する３つの単位ベクトルｎ、ａ、ｂは、座標変換によって次のように表される。

【数27】

平面Δ上の任意の単位ベクトルｑは、次のように表される。

【数28】

任意角度の正応力σ_ｎ（ｔ）、せん断応力τ_ｑ（ｔ）、せん断歪みγ_ｑ（ｔ）は、次のように表される。

【数29】

任意平面上の正応力平均値σ_ｎ，ｍを［０，Ｔ］区間における正応力積分の平均値に定義し、正応力振幅σ_ｎ，αは、正応力の分散によって定義される。

【数30】

式中、

【数31】

せん断応力、せん断歪みの平均値及び振幅τ_ｍ，τ_α，γ_ｍ，γ_αを定義する。

【数32】

式中、

【数33】

せん断歪みエネルギーΔτΔγは、部材の安定サイクル下でヒステリシスループを囲む最小矩形、すなわちせん断応力変化幅Δτとせん断歪み変化幅Δγとの積として定義し、変化幅Δτ＝２τ_α、Δγ＝２γ_α、そこで、ΔτΔγ＝４τ_α×γ_α。

【0079】

（付記４）
多軸応力比ρは、次のように表されることを特徴とする、
付記３に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法。

【数34】

ここで、σ_{ｎ，ｍａｘ}は、臨界面における正応力の最大値であり、σ_ｎ，ｍは、正応力平均値であり、σ_ｎ，αは、正応力振幅であり、τ_αは、せん断応力振幅である。

【0080】

（付記５）
前記ステップ（３）における単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線は、次の方法に基づいて決定されることを特徴とする、
付記１に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法。
単軸引張圧縮、ねじり荷重下の各疲労データ点臨界面上のせん断歪みエネルギーをそれぞれ計算し、首尾連係してせん断歪みエネルギー－寿命曲線を描画し、べき関数で２種類の荷重状況下でせん断歪みエネルギー－寿命曲線の各パラメータＡ_１，Ｂ_１，Ａ_３，Ｂ_３をそれぞれフィッティングして解く。
２つの荷重におけるせん断歪みエネルギー寿命曲線関数は、次のように表されることを特徴とする。

【数35】

ここで、ΔτΔγは、せん断歪みエネルギーであり、Ｎ_ｆは、寿命である。

【0081】

（付記６）
前記ステップ（３）において正規化された仮想歪みエネルギーＷ_Ｖは、次のように表されることを特徴とする、
付記５に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法。

【数36】

ここで、

【数37】

は、正規化係数であり、ρは、多軸応力比である。

【0082】

（付記７）
金属部材の幾何学モデル及び外荷重を入力するための入力モジュールと、
処理モジュールと、
所定の外荷重における金属部材の予測寿命を表示するための出力モジュールを含み、
前記処理モジュールは、
金属部材に対して有限要素分析を行い、危険点を確定し、危険点における応力歪み荷重時間の経過を求めるための応力歪み計算ユニットと、
危険点を通過した各平面上のせん断歪みエネルギーを計算し、せん断歪みエネルギーが最大の平面を臨界面とし、臨界面上のせん断歪みエネルギー及び多軸応力比を得るための臨界面計算ユニットであって、ここで、せん断歪みエネルギーは、せん断応力変化幅とせん断歪み変化幅の積として定義され、多軸応力比は、臨界面上の正応力最大値とせん断応力振幅の比として定義される臨界面計算ユニットと、
危険点における仮想歪みエネルギーを計算し、ねじり荷重のせん断歪みエネルギー寿命曲線を用いて寿命予測を行うための予測ユニットであって、ここで、前記仮想歪みエネルギーは、単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線におけるせん断歪みエネルギーを正規化処理した後の歪みエネルギーである予測ユニットと、を含むことを特徴とする、
仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測システム。

【0083】

（付記８）
前記単軸引張圧縮荷重、ねじり荷重下のせん断歪みエネルギー寿命曲線は、べき関数フィッティングユニットによって生成され、
前記べき関数フィッティングユニットは、単軸引張圧縮、ねじり荷重下の各疲労データ点臨界面上のせん断歪みエネルギーをそれぞれ計算し、首尾連係してせん断歪みエネルギー－寿命曲線を描画し、べき関数で２種類の荷重状況下でせん断歪みエネルギー－寿命曲線の各パラメータＡ_１，Ｂ_１，Ａ_３，Ｂ_３をそれぞれフィッティングして解き、
２つの荷重におけるせん断歪みエネルギー寿命曲線関数は、次のように表されることを特徴とする、
付記７に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測システム。

【数38】

ここで、ΔτΔγは、せん断歪みエネルギーであり、Ｎ_ｆは、寿命である。

【0084】

（付記９）
前記予測ユニットで正規化された仮想歪みエネルギーＷ_Ｖは、次のように表されることを特徴とする、
付記８に記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測システム。

【数39】

ここで、

【数40】

は、正規化係数であり、ρは、多軸応力比である。

【0085】

（付記１０）
メモリ、プロセッサ、及びメモリ上に記憶されプロセッサ上で動作可能なコンピュータプログラムを含むコンピュータシステムであって、
前記コンピュータプログラムがプロセッサにロードされると、付記１～６のいずれか１つに記載の仮想歪みエネルギーに基づく金属材料の多軸疲労寿命予測方法を実現することを特徴とする、
コンピュータシステム。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【国際調査報告】