長岡工業高等専門学校研究紀要 第 52 巻（2016）論文

有限要素解析に用いる三角形メッシュの自動生成法－き裂・切り欠き先端近傍の応力集中を考慮したメッシュ生成－

倉橋 貴彦1・鋤柄 あかね2・山際 健吾3・丸岡 宏太郎4・井山 徹郎5
1 長岡技術科学大学大学院（Department of Mechanical Engineering, Nagaoka University of Technology）
2 長岡技術科学大学（Department of Mechanical Engineering, Nagaoka University of Technology）
3 長岡技術科学大学大学院（Department of Mechanical Engineering, Nagaoka University of Technology）
4 長岡技術科学大学大学院（Department of Mechanical Engineering, Nagaoka University of Technology）
5 機械工学科（Department of Mechanical Engineering, National Institute of Technology, Nagaoka College）

Automatic Mesh Generation of Triangular Meshes
in Finite Element Analysis-Mesh Generation Considering Stress Concentration near Crack and Notch TipsTakahiko KURAHASHI1, Akane SUKIGARA2, Kengo YAMAGIWA3 ,Kotaro MARUOKA4 and Tetsuro IYAMA5Abstract
In this study, we present triangular mesh generation method considering stress concentration near internalcrack tip. The total number of nodes and elements affect the computational time, and it is desired that thelocation of nodes is distributed such that distribution of physical variables is appropriately expressed. In
this study, we focus on the distribution of stress components near crack tip, and propose the automatic meshgeneration method considering the stress concentration. In addition, the stress analysis using the finite
element meshes obtained by the present method is carried out, and the stress distribution near crack tip iscompared with that obtained by using the other meshes. Furthermore, the stress analysis for V-notchedplate is carried out using the finite element meshes obtained by the present method.
Key Words : Automatic mesh generation, Triangular mesh, Finite element method, Stress analysis, Stressconcentration

１．序論
従来より，構造部材のき裂先端や接合部材の接合界面端近傍における応力集中部の解析が行われている 1), 2)．応力集中部では応力成分の値が急激に上昇をするため，応力解析のために部材に対するメッシ- 11 -

倉橋 貴彦・鋤柄 あかね・山際 健吾・丸岡 宏太郎・井山 徹郎

ュ分割をする際にも応力集中部近傍では特にメッシュを細かく分割する等の配慮が必要となる．また，
応力集中部先端の点からの距離をr，応力集中部における応力の増幅率を表す特異性のオーダをλとす
ると，応力集中部近傍における応力分布はσij∝r-λにより表される 1), 2)．き裂の問題の場合はλの値は0.5となり，異なる材料の接合部材においては，材料の特性（ヤング率・ポアソン比）および接合部の形状によってこの値は変わることになる．ここに，記号∝は比例関係を示す．このため，この関係を考慮したメッシュ分割を行うことにより，応力集中部近傍における応力分布を適切に解析できるものと考えられる．そこで本研究では，有限要素法に基づく応力解析を考え，有限要素メッシュの分割法3), 4)に着目し検討を行う．

２．メッシュ生成アルゴリズム
2.1 応力集中部の応力分布を考慮したメッシュ間隔
応力集中部においては，応力の最も高い点（特異点）からの距離をr，特異点周りの角度をθ，特異性のオーダをλとすると応力場は式(1)のように表すことができる．

σ ij (r, θ ) = K ij f (θ )r −λ + (Other terms)(1)

ここに，Kijは特異応力場の強さ，f(θ)は角度依存関数，(Other terms)は式(1)の左辺第一項のみによって表すことができない項の集まりを示す．同様にひずみ場も式(2)のように書くことができる．

ε ij (r, θ ) = Aij g (θ )r −λ + (Other terms)(2)

立てると，応力集中点近傍の応力成分に特異性の生じる特異性のオーダλの範囲は0<λ<1であることから，式(3)のように特異点から離れるに従い上昇する関数に変えることができる．

∫ ε dr = ∫ A g(θ )r dr−λrrij= ∫ α (θ )r −λ dr =α (θ ) −λ +1r+C− λ +1(3)仮に，式(3)により刻み幅Δrが与えられるとすると，式(4)のように表される．

Δr =α (θ ) −λ +1r+C− λ +1(4)

式(4)について，き裂を有する問題を考えると，き裂先端の特異点近傍における特異性のオーダλは0.5であるため，刻み間隔Δrは式(5)のように与えられる．

Δr =α (θ )− 0.5 + 1

r −0.5+1 + C = 2α (θ ) r + C(5)

積分定数Cはメッシュ生成時に仮定するものとし，き裂を有する材料モデルに対しては式(5)に示す評
価式をもとに，また異材接合体や切り欠きを有する材料モデルに対しては式(4)に示す評価式をもとにメッシュ分割を行う方法について検討を行う．
2.2 メッシュ生成のアルゴリズム
本論文におけるメッシュ生成の流れは大きく7段階に分けられており，メッシュ生成の計算は以下の手順により行われる．
Step1：対象領域の外形点の座標値（節点）を入力

特異点近傍において，式(2)により表されるひずみεrrの分布は，特異点に向かうに従い累乗関数的に上昇する関数として表すことができる．また，特異点近傍における応力解析の精度を上げるため，特異点近傍の刻み幅Δrを小さく，特異点から離れるに従いΔrが大きく与えられるように，特異点からの距離rの関数として式(2)を適用することを考えると，式(2)は特異点に向かうに従い累乗関数的に上昇する関数であることから，Δrの計算式として，直接的な適用は困難なものとなる．そこで，一つの考えとして，ひずみ成分εrrの半径r方向に関する積分式をする．（図-1）
Step2：ステップ1において入力した節点により粗い三角形要素を生成する．（図-2）
Step3：式(4)（き裂を有する問題では式(5)）をもとに，境界線上に細かい節点を生成する．
Step4：ステップ3と同様に領域内の節点を生成する．
（図-3）（※節点の生成方法の詳細は本アルゴリズ
ムの後に説明する．）
Step5：ステップ3，4において生成された節点をもとに細かい三角形要素を生成する．（図-4）
Step6：ステップ5において生成された三角形要素に- 12 -

有限要素解析に用いる三角形メッシュの自動生成法

おいて，各節点に対して，節点の再配置をN回繰り返す．（図-5）（※節点の再配置の詳細は本アルゴリズムの後に説明する．）
Step7：Step6において出力されたメッシュをメッシュの再集形とする．（図-6）

図-3

図-1

図-2

Step3,4 における節点の生成

Step1 で入力された境界の節点と内部節点

図-4

Step5 において作成される細かい三角形メッシュ

図-5

Step6 の節点再配置の過程(き裂先端部拡大図)

Step2 で作成される粗い三角形メッシュ- 13 -

倉橋 貴彦・鋤柄 あかね・山際 健吾・丸岡 宏太郎・井山 徹郎

C(x3,y3)δΔr(x,y)δ
B(x2,y2)

P(x,y)

図-6

A(x1,y1)

提案手法により生成された三角形メッシュ： New interior nodes： Nodes of rough mesh

ここでStep4およびStep6の詳細について説明する．
Step4における粗い三角形メッシュ内における細か図-7い節点を生成するための説明図を図-7に示す．節点の生成の仕方としては，まず仮の細かい節点を間隔δごとに等間隔に配置する．（本検討ではδ=0.1としている．）次に，各々の仮の細かい節点上において，特異点から半径方向に対するひずみ成分εrrの変化量が一定となるように計算したΔr（※Δrの計算には，異材接合体および切り欠きを有する材料モデルの問題に対しては式(4)を適用し，き裂を有する材料モデルの問題に対しては式(5)を適用する．）を計算する．計算の方法としては，図-7の菱形の点上において与えられているΔr（粗いメッシュの各節点上におけるΔr）の値を用いて式(6)のように線形補間する．

Δr ( x, y ) = N A ΔrA + N B ΔrB + N C ΔrC

新しい節点の生成方法の説明図

また，Step6における節点の再配置（節点座標値の修正）については，図-8に示すように，ある対象としている節点Pを取り囲む要素の各重心点座標を用いて，節点を取り囲む要素の面積を使用した加重平均により節点座標値の修正を順次行う．全節点における座標の再計算を1回とし，この操作をN回繰り返す．

(x4,y4)(x3,y3)(6)(x1,y1)

ここに，NA，NB，NCは三角形一次要素における形状関数を示す6)．ここで，図-7に示すように，ある対象とする点Pにおいて，Δrの円を描き，その中に節点が新しく生成した節点が含まれていなければ新しい節点とする．この操作を順次繰り返すことにより，粗いメッシュ内において新しい節点が生成されることになる．

P(xP,yP)(x2,y2)図-8

節点の再配置の説明図

本論文におけるメッシュ生成では，式(4) および式(5)を要素サイズ関数とし，パラメータα(θ)の値は定- 14 -有限要素解析に用いる三角形メッシュの自動生成法

数として0.1，積分定数Cの値は0.1と設定する．また，Step6における節点再配置の回数Nは3とする．

３．提案メッシュの適切性に関する検証
図-9にき裂を有する材料の1/4モデル図を示す．
サイズおよび境界条件は図-9に示すとおりであり，x=10mm,y=0mmの点（特異点）から，θ=0°の半径r方向に対する応力成分σyyの分布について，メッシュを変えたモデルに対して比較を行う．材料は，鋼を想定し，表-1に示す材料定数を用いる．
図-6に示す非構造メッシュと比べて，節点数はおおむね等しく且つ一様分割した構造メッシュ（Case1）と，特異点近傍も細かくした構造メッシ
ュ（Case2）を用意する（図-10，図-11）．図-6に示すメッシュをCase3とし，節点数，要素数，最小メッシュを表-2に整理する．
図-9に示す計算条件のもとCase1，Case2に対して応力解析をMarc（MSC Software）により行った．応図-9計算モデル図表-1

材料定数2000.3

ヤング率, GPa
ポアソン比

力成分σyyの分布を図-12，図-13に示す．特異点近傍において，応力成分σyyの値は急激に高くなる傾向にあり，メッシュの細かさによっても，得られる値の違いがあることがわかる．
また，図-6に示すメッシュをCase3とし，Fortranにより作成した線形弾性体の応力解析のプログラムにより計算した結果との比較を行う．特異点を原点とし，θ=0°の半径r方向に対する応力成分σyyの分布
図の比較を図-14に示す．また，各ケースに対して特異点近傍における応力成分σyyの分布をσyy=Kyyr-λによる累乗近似を行い，応力拡大係数Kyy および特異性のオーダλを整理したものを表-3に示す． 結果として，き裂先端近傍における特異性のオーダの理論値である0.5に最も近い結果はCase2の値であるが，本論文において提案するメッシュを用いて計算したCase3の結果もCase2に非常に近い値を示していることがわかる．また，Case3は，他のケースに比べて最も節点数が少ないメッシュであるため，最も計算時間を必要とせず応力成分σyyが得られている．このような点から，本論文において紹介したメッシュ生成法は，応力集中を有する問題に対する解析においては非常に有効であることがわかる．

図-10計算メッシュ図 (Case1)

図-11

計算メッシュ図 (Case2)
表-2

節点数
要素数
最小メッシ
ュ寸法,mm- 15 -

計算条件

Case1

Case2
（参照解）

1013
1936
1.286

3671
7200
0.189

Case3
（提案メッシュ）870
16270.131

倉橋 貴彦・鋤柄 あかね・山際 健吾・丸岡 宏太郎・井山 徹郎

４．切り欠きを有する材料に対する応力解析
次に，切り欠きを有する材料に対するメッシュ分割 に 関 す る 検 討 結 果 を 示 す ． 図 -9に お い て は ，θ=180°のき裂モデルであったが，図-15に示すようにθ=170°，160°，150°，140°の切り欠きを有するモデルに対して特異性のオーダλを考慮したメッシュを作成し，応力解析結果の妥当性について検証を行う．なお，図-9のように1/4領域の対称モデルとし，矩形領域の中心に菱形の空洞を有する解析モデルを想定している．特異性のオーダλは，Bogyの特性方程式5)の特性根pを求め，特異性のオーダλと特性根pがλ=1-pとなる関係から求めたものである．
切り欠きを有するモデルではき裂を有するモデルと特異性のオーダλの値が異なるため，式(4)に示す評
価式のパラメータα(θ)の値は定数として0.1，積分定数Cの値は0.1と設定し，各切り欠きにおける特異性のオーダλの値を変えることで，2.2節に基づく方法
によりメッシュ分割を行った（式(7)）．

応力成分 σyy の分布図 (Case1)

図-12Δr =(7)0.1 −λ +1=r+ 0.1− λ +1

応力成分 σyy の分布図 (Case2)

図-13α (θ ) −λ +1r+C− λ +1

1.00E+01
Case1:粗い構造メッシュ

Case3
σyy = 3.0529r‐0.438

Case2:メッシュの粗密も考慮した構造メッシュ

σyy, MPaCase2
σyy = 3.0827r‐0.446

Case3:応力分布を考慮した非構造メッシュ

Case1
σyy = 2.8971r‐0.416

累乗 (Case1:粗い構造メッシュ)
累乗 (Case2:メッシュの粗密も考慮した構造メッシュ)

1.00E‐01

累乗 (Case3:応力分布を考慮した非構造メッシュ)

1.00E+00
1.00E+00

1.00E+01

distance from crack edge "r", mm

図-14Case A : θ=170°，λ= 0.4996

θ=0°の半径 r 方向に対する応力成分 σyy の分布図

表-3

（総節点数：839，総要素数：1567）

各メッシュ条件に対する累乗近似により求まった特異性のオーダλと特異応力場の強さKyyの比較
特異性のオーダλ
応力拡大係数Kyy

Case1
0.416
2.8971

Case2
0.446
3.0827（図-15-1）

Case30.438
3.0529- 16 -有限要素解析に用いる三角形メッシュの自動生成法

応力解析における検討では，図-9に示す境界条件（荷重条件・対称面を考慮した変位に対する固定条件）および表-1に示す材料定数を用いることにする．
CaseA～Dのモデルにおいて，特異点を原点とし，θ=0°の半径r方向に対する応力成分σyyの分布図の比較を図-16に示す．図は両対数のグラフにより示しており，十字点で表した応力解析結果および図-15に示す特異性のオーダλの値を用いたσyy=Kyyr-λ による最小二乗近似曲線を示す．全てのケースにおいて，応力解析結果と最小二乗近似曲線は良い一致を示しているもの（精度の良い応力解析結果が得られている）と考えられる．
Case B : θ=160°，λ= 0.4965（図-15-2）（総節点数：804，総要素数：1498）

Case C : θ=150°，λ= 0.4878（図-15-3）（総節点数：762，総要素数：1416）

Case A : θ=170°，λ= 0.4996，Kyy= 2.8575 MPamm0.4996（図-16-1）

Case D : θ=140°，λ= 0.4696（図-15-4）（総節点数：710，総要素数：1313）
図-15

切り欠き先端の特異性のオーダ λ を考慮し作成した有限要素メッシュ- 17 -倉橋 貴彦・鋤柄 あかね・山際 健吾・丸岡 宏太郎・井山 徹郎

Case B : θ=160°，λ= 0.4965，Kyy= 3.0558 MPamm0.4965

Case D : θ=140°，λ= 0.4696 ，Kyy= 3.8379MPamm0.4696（図-16-2）（図-16-4）
図-16

各ケースに対する応力解析結果

５．結論
本論文では，応力集中時における応力分布を考慮した三角形メッシュの自動生成法について検討を行った．構造のき裂先端近傍における応力場・ひずみ場の計算式より，ひずみ場を特異点からの距離で積分した値により，三角形メッシュの代表長さを決定する方法について検討を行った．結果として，き裂を有する材料モデルに対する応力解析では，節点数が少ないにも関わらず，精度の良い応力分布が得られることがわかった．また，切り欠きを有する材料モデルに対する応力解析では，切り欠きの角度を変えた全てのケースにおいて，応力解析結果と最小二乗近似曲線は良い一致を示すことを確認できた．
Case C : θ=150°，λ= 0.4878，Kyy= 3.2614 MPamm0.4878（図-16-3）

付録

クラック付き材料に対する引張り試験

クラック付きアクリル板および塩化ビニル板に対する引張り破壊試験を行い，破壊荷重に関する比較を行った．図-17に示すように，各種材料に対して100mm×25mmの試験の両端に5mmのクラックを入- 18 -

有限要素解析に用いる三角形メッシュの自動生成法れた表面き裂の試験片を用意し，引張り試験を行った．実験結果を図-18，図-19に整理する．
結果より，アクリル板の場合，破断時に荷重は瞬時に零になっていることに対して，塩化ビニル板の場合はなだらかに荷重が下がっていることがわかる．
図-19の結果は，塑性変形をしながら破断に至っているものと考えられ，将来的には，このような塑性変形も考慮した応力解析および応力拡大係数に関する検討が必要になると考えられる．
2回目
謝辞：本論文の解析結果は，九州大学情報基盤研究開発センターの高性能アプリケーションサーバ
SR16000を使用し計算を行ったものである．計算機（図-18-2）

図-18 アクリル板（表面き裂の場合）に対する実験結果の使用について，センター関係者の方へ謝意を表す．

図-17

1回目

（図-19-1）

2回目（図-19-2）

図-19 塩化ビニル板（表面き裂の場合）に対する実験結果

試験片および表面き裂の寸法図

参考文献1)

倉橋貴彦，中島正人，石川晃広，星和久，古口日出男， 三次元異材接合体界面端における臨界特異応力場の強さ（接着形状が矩形の場合）, 日本機械学会論文集 A 編, 78 巻, 794 号, pp.1382-1399, 2012.

2)

倉橋貴彦，渡辺裕太郎，近藤俊美，古口日出男， 三
次元特異性オーダを用いた Akin 特異要素による異材接合体の界面端角部近傍における特異応力場の評価

1回目（四面体要素と特異要素使用時の結果の比較），日（図-18-1）- 19 -倉橋 貴彦・鋤柄 あかね・山際 健吾・丸岡 宏太郎・井山 徹郎

本機械学会論文集, 第 80 巻，第 809 号 pp.1-18, 2014.
3)

S.H.Lo, A New Mesh Generation Scheme for Arbitorary
Planar Domains, International Journal for Numerical
Methods in Engineering, Vol.21, pp.1403-1426, 1985.

4)

K.Kashiyama and T.Okada, Automatic Mesh Generation
Method for Shallow Water Flow Analysis, International
Journal for Numerical Methods in Fluids, Vol.15, pp.10371057, 1992.

5)

D.B.Bogy, Two Edge-Bonded Elastic Wedges of Different
Materials and Wedge Angles Under Surface Tractions，
Journal of Applied Mechanics, Vol.38, pp.377-386, 1971.

6)

竹内則雄，寺田賢二郎，樫山和男，計算力学 有限要素法の基礎，森北出版(株)，2003．
(2016. 9. 16 受付)- 20 -