平成１３年度
卒業論文

特異要素を用いたき裂先端応力場の有限要素解析

高知工科大学工学部
知能機械システム工学科
材料強度学研究室
藤原 誠敬

目次
第１章 序論……………………………………………… １

第1節 はじめに第2節 本研究の目的と概要
第3節 有限要素法の特徴
第２章 応力拡大係数……………………………………… ２

第1節

き裂の変形様式

第2節

節応力分布

第3節

次元解析

第4節

有限要素法による K 値の解析

第３章 有限要素法による K 値の解析
第1節

有限要素法の概要

第2節

K の解析

第3節

Barsoum の特異要素……………………… ５

第４章 結果と考察 ……………………………………………… ８
第1節

対象領域

第2節

境界条件

第3節

要素分割

第4節

解析結果及び考察第５章 緒言 ……………………………………………………… ２４第１章 緒言
第１節 はじめに構造物や機械の設計・使用に際して強度安全性を推測する方法として従来は材料力学が重要な役割を果たしてきた。
しかし、今日における構造物や機械類の大型化高速化に加え過酷な使用条件が強いられるようになり、従来の材料力学的手法では解決仕切れない問題が生じるようになった。材料力学的手法では応力集中率や切欠き係数と言った概念から想像できるように、切欠き近傍の応力集中を基本に考えられているため鋭いき裂状の欠陥に対して何ら考慮がなされていない。この重要な問題を解決するために発展したものが破壊力学である。
すなわち、破壊力学はき裂状の欠陥を含む構造物あるいは機械を安全に使用するための１つの理論体系である。
破壊力学は、破壊機構の解明を主眼とし、き裂を含む部材の応力解析と破壊じん性試験との２つから構成されている。き裂先端の弾性応力場を特性づけるパラメータ応力拡大係数 K を用いて、ぜい性破壊を防止する条件を考慮しようとしたものである。この方法は弾性解析に基づくパラメータ K により破壊条件が記述されるので、き裂先端近傍に生じる塑性域がき裂や他の部材寸法に比べて小さい小規模降伏の範囲で有効である。破壊力学は、設計、材料の選択あるいは破壊事故の診断等に際して広く使用されるようになった。一方でこの破壊力学は、有限要素法の研究者たちが活発に応用を試みた分野の１つで、近年のその技術の進歩には著しいものがある。現在では、応力解析と材料試験からなる破壊力学において、有限要素法は数値解析を行うための重要な手段となっている。
第２節 本研究の目的と概要
前節で述べたように、破壊力学に応用でき、複雑なき裂先端近傍の応力解析に適している有限要素法を理解し、特に特異要素を適応したき裂先端応力解析により破壊において重要なパラメータである応力拡大係数を求めき裂先端近傍応力場の解析を行うことが目的である。本研究では、八節点四角形アイソパラメトリック要素を用い、き裂先端に Barsoum の特異要素を適応し、中央切欠き試験片における K の算出を行う。さらに、三角形要素と比較し特異要素の有用性の評価を行う。
第３節 有限要素法の特徴
有限要素法とは、与えられた偏微分方程式の境界条件問題を、有限要素という離散化近似を行うことにより、高次元の連立１次方程式の形に変換して解く一種の近似法である。偏微分方程式の境界条件問題として記述される物理現象は極めて多用なため、使用範囲も広範にわたっている。特にき裂の場合、き裂を含む試験片や構造部材は一般に極めて複雑な形状をしており、これを理論的に解くことは困難である一方、有限要素法はこのような複雑な境界条件においても容易に解を求めることができる。き裂先端近傍の特異性に対して取り扱いの困難さが従来の問題であったが研究の進歩より特異性を含む解析を行うことが可能である。弾性問題だけでなく弾塑性問題まで拡張でき、理論的に不可能といってもよい非線型問題への応用も可能である。

1

第２章 応力拡大係数
第１節 き裂の変形様式
き裂先端の応力場は基本的な３つの変形様式に分類することが出来る。すなわち、模式的に図 2-1.1 に示すように開口型（モードⅠ）、面内せん断型（モードⅡ）および面外せん断型（モードⅢ）の３つである。これらの３つの変形様式を適当に重ね合わせることにより、き裂先端近傍の任意の変位場および応力場を記述することが出来る。
vyyyxxxuzzz
モードⅢモードⅠモードⅡw

図 2-1.1 き裂先端様式
第２節 き裂先端近傍応力分布
先に述べたモードⅠとモードⅡは弾性問題の平面問題として解析することができ、き裂面に関してそれぞれ対称部分とに分けられる。また、モードⅢは純粋せん断あるいはねじり問題と考えることができる。座標及び応力成分に関して図 2-2.1 の記号を用いると、各モードにおける応力は以下のようにあらわされる。

図 2-2.1 き裂先端の局所座標と応力成分2

モードⅠKⅠθθ3θ cos 1 − sin sin 2
222πrθθK3θ σ y = Ⅰ cos 1 + sin sin 222 2πrθθ3θKτ xy = Ⅰ sin cos cos
2222πrσ z = ν (σ x + σ y ), τ xy = τ yz = 0

σx =(2-2.1)
モードⅡKⅢθsin
22πrKθτ yz = Ⅲ cos
22πr
σ x = σ y = σ z = τ xy = 0

τ xz = −(2-2.2)
モードⅢKⅡθθ
3θ sin  2 + cos cos 2
22
2πr rKθθ3θ
σ y = Ⅱ sin cos cos
2222πrσx = −θθ
3θ KⅡcos 1 − sin sin 222
2πrr
σ z = ν (σ x + σ y ), τ xy = τ yz = 0

τ xy =(2-2.3)
この式中の KⅠ, KⅡ , K Ⅲ はそれぞれのモードにおける応力場に対応する応力拡大係数である。応力拡大係数は
局座標 r ,θ には依存しない。したがってそれぞれの変形様式に対して応力場の強さに影響するが、分布の形には影響を与えない。弾性体の場合、応力拡大係数は外力の大きさに比例し、き裂を含む部材の形状に依存することがわかる。応力拡大係数は、部材にき裂を導入したことによる応力の再分配、き裂の変形様式とき裂先端部を通過する力の大きさとをあらわすパラメータであることがわかる。

3第３節 次元解析
通常のき裂近傍の任意の応力場を記述する場合３つの基本的なき裂様式を適当に重ね合わせることにより求めることができるが、実験レベルの場合は単一のモードに限った実験が一般的である。
σ2aσ

図 2-3.1 垂直応力を受ける無限平板中のき裂

切欠の場合、寸法や応力が異なる時、切欠底の最大応力あるいは応力集中係数が切欠先端の力学的厳しさの尺度の一つとなるが、き裂の場合には最大応力は、遠方の応力、き裂長さに関係無く無限大となるので、厳しさの尺度として最大応力あるいは応力集中係数を用いることは適切ではなく、応力拡大係数のような別の尺度を必要とする。図 2-3.1 に示すようように長さ 2a のき裂を持つ無限板の遠方に、き裂に垂直な方向に一様な応力σが作用している場合を考える。き裂の形状は面に関して対称であるため、モードⅠのき裂様式のみが適応される。したがって式（2-2.1）に関する簡単な次元解析により

KⅠ = C1σ a , KⅡ = K Ⅲ = 0(2-3.1)が求められる唯一の結果となる。
対称性の考察あるいは次元解析は応力拡大係数を決定する手段となっている。
このような考察のみからは、C1 を決定することはできないが、 π であることがわかっている。しかしたとえ
C1 が決まらなくても、この形状に関する破壊の寸法効果は予測することができる。応力拡大係数 K は、（ 応力× き裂長さ ）で単位は MPa m なる。

4

第 3 章 有限要素法による K 値の解析
第１節 有限要素法の概要
2 次元の領域 D、境界条件 C で定義される物体が空間に固定されている、これに外力が作用した場合、弾性の基礎方程は式(3-1.1)〜(3-1.5)となる。

σy + ∂σy dy∂yyτxy + ∂τxy dy∂yν(cosα,sinα)CσYσyD

τxyσxσxdy
τxydx

Cu

σx + ∂σx dx∂xXτxy + ∂τxy dx
∂x

τxyx

Oσy(a)(b)
図 3-1.1 弾性問題における座標系つり合いの方程式∂σ x ∂τ xy+ X = 0,+
∂y∂x

∂τ xy
∂y+

∂σ y
∂y+Y = 0(3-1.1)( X , Y : 体積力)

変位とひずみの関係

εx =∂v ∂u∂x∂u, ε y = , τ xy =+∂x ∂y∂y∂x(3-1.2)構成方程式（Hooke の法則）

εx =

τ xy
11(σ x − ν 'σ y ), ε y = (σ y − ν 'σ x ), τ xy =E'E'
G(3-1.3)

境界条件(i) 幾何学的境界条件

C u 上で u = u , v = v(3-1.4)(ii)力学的条件

Cσ 上で X ν = X ν , Yν = Yν

(3-1.5)式(3-1.1)および式(3-1.5)を付帯条件式(3-1.2)および式(3-1.4)の下に書き直すと、いわゆる仮想仕事の原理が得られる。

∫∫ (σ δε + σ δε + τ δγ ) − tdxdx − ∫∫ ( Xδu + Y δv)tdxdy − ∫ σ ( X ν δu + Yν δv)tds = 0
DxxyyxyxyDC(3-1.6)5

ここにおいて t は板厚、ds は境界を指す。また付加される仮想変位δu , δv は式(3-1.2)を満足し、境界条件 Cu上でδu = 0, δv = 0 という条件を満足する許容変位である。式(3-1.6)を導くに当たって構成方程式は持ちられていないので、この式は弾塑性体に対しても適用可能である。
有限要素法は式(3-1.1)および(3-1.5)の変わりに式(3-1.6)を用いて定式かされる。この定式かは以下に述べる
ように、節点変位を未知パラメータとして扱っているため変位法と呼ばれている。
次に有限要素法において最も重要な仮定、離散化近似を行う。領域を有限個の要素に分割し、個々の要素内における変位のモードを仮定する。これは、一般には次のように表す。

{ f } = u, v T = [N ]{δ }m , {δ }Tm = u i , vi , u j , v j , u k , v k (3-1.7)

ここに{δ }m は要素 m の節点変位ベクトル、[N]は形状関数あるいは変位関数である。要素および形状関数の選択が解の精度に大きな影響を与える。
式(3-1.7)を変位とひずみの関係式(3-1.2)に代入して整理すると{ε } = ε x , ε y , γ xy T = [B]{δ }m(3-1.8)

ここに[B]の成分は形状関数の座標微分により求められる。
次に応力とひずみの関係は弾性体に対してはHookeの法則、式(3-1.3)をマトリックスで表して、

σ x 1 ν'0  ε x − ε x 0 E' ν ' 10  ε y − ε y 0 
σ y  =2τ  1 − ν '  0 0 1 − ν '  ε − ε z0  z xy 0 


(3-1.9)あるいは、簡単に{σ } = [D ]({ε } − {ε 0 })ここに{ε 0 }は初期ひずみを表す。式(3-1.3)および(3-1.4)を式(3-1.6)に代入して整理すると{F }m = [K ]m {δ }m + {Fε 0 }m + {Fp }m(3-1.10)
ここに[K ]m = ∫∫ [B]T [D][B ]tdxdy(3-1.11){Fε 0 }m = − ∫∫ [B]T [D]{ε 0 }tdxdy(3-1.12){F } = − ∫∫ [N ] {p}tdxdy, {p} = X , Y TTp m(3-1.13){ }また F p m は要素の等価節点
第 2 節 K の解析
有限要素解析で得られたき裂先端近傍の応力を直接式(2-2.1)〜(2-2.3)に代入してK 値を求める方法が直接応力法である。この時、代入する値は変形モードに対して支配的な成分を選択する。モードⅠ型のき裂に対しては、x 軸上の y 方向の垂直成分応力σ y を代入する。有限要素法で節点変位を未知パラメータとする変位法が6応力法に比べると精度の点で信頼度が高い。したがって、これらの値を代入して得られるＫの値に関しても同様のことがいえる。
直接法によりＫ値を求める場合、有限要素法で得られた解がき裂先端近傍において十分に精度のよいものでなければならない。一般に使用されている要素を用いた場合、式(2-2.1)〜(2-2.3)で表わされるような r−1 / 2の特異性を表現することはできない。そのため、き裂先端近傍における応力の精度はそれほど期待できない。したがって、これらよりより求められるＫ値の精度もあまりよいとはいえない。
以上より直接 K 値を求める場合解の精度を改善するために一般的に以下に示すような手法がとられている。
1. き裂先端近傍での要素分割をできるだけ細かくする。
2. き裂先端から距離 r の異なる点で式(2-1.1)〜(2-1.3)より K 値を計算し、K と r の関係をプロットして r→0に外挿した点をもってこのき裂の K 値とする。
（外挿法）
3. き裂先端で r−1 / 2の応力の特異性を表現できる特異要素を使用する。

第 3 節 Barsoum の特異要素
２次元問題において、８節点アイソパラメトリック要素では図 3-1.1(a)で示すように辺の中点に節点がある。
この要素において節点 1、8、4 を１つの点とするとともに節点 5、7 をそれぞれ辺の４分の１の点へずらすことにより、r−1 / 2, r −1 のひずみの特異性を同時に持つ要素を作ることができる。この要素は弾性体および加工 k

硬化のない完全弾性体におけるき裂を有する解析を行うことができる。このアイソパラメトリック要素は、プログラム段階で特別な工程は必要とせず、データを準備する段階において図 3-1.1(b)に示すような節点と座標を選ぶだけである。

37437868L24L1/4161253L1/45L1(a)(b)
図 2-3.1 Barsoum の特異要素

72

第４章 結果および考察
第１節 対象領域
本研究では、疲労き裂で伝ぱ試験に使用されている中央切欠き試験片について、Ｃ言語により作成したプログラムを用いて応力場および応力拡大係数を解析した。図 4-1.1 に寸法を示す。板厚は 5(mm)である。試験片中央には長さ 2a のき裂が存在する。

2a34

図 4-1.1 試験片の形状 (mm)
図 4-1.1 に示す対象領域は試験片に対して一点鎖線で区切られら領域はすべて対象であるため、実際に計算する解析対象領域は４分の１のみとする。従って解析領域は図 4-1-2 となる。

17

図 4-1.2 解析領域 (mm)
また、試験片の材料特性に関する諸定数には、ヤング率 E=206(GPa),ポアソン比ν=0.3 を用いた。
8

第２節 境界条件23.5(MPa)の引張り負荷荷重をあたえた。境界条件としては、解析領域が 4 分の１であることを考慮し、図4-2.1 に示すように、Y 軸対称軸上の接点は x 方向に、き裂面対称軸上の接点は y 軸に対しての拘束する。
σy : 23.5(MPa)

き裂先端

4-2.1 境界条件
第３節 要素分割
本解析に用いた要素分割を説明する。要素分割には八節点アイソパラメトリック要素を用い、き裂先端には特異要素を適用する。例としてき裂長さ 5(mm)の場合の要素分割図 4-3.1〜4.3.3 を示す。この時点の節点数は722 固、要素数は 217 個である。

図4‑3.2

図 4-3.1 メッシュ分割全体9

図4‑3.3

図 4-3.2 メッシュ分割き裂先端付近拡大

図 4-3.3 メッシュ分割き裂先端

図 4-3.2 は図 4-3.1 からき裂先端付近の要素分割を拡大して示したものである。また、図 3-3.3 はさらに拡大したものである。
図 4-3.3 のき裂先端には Barsoum の特異要素を用いて要素分割を行っている。本解析では、き裂の長さを5,6,8,10(mm)と変えて解析を行った。き裂長さを変える場合、図 4-3.2 の要素分割の部分をスライドし、空いた領域に新しい要素を作ることにより対応した。
10

比較のために三角形要素を用いての有限要素解析も併せて行った。この時の、要素分割を図 4-3.4 に示す。三角形用は特異性を表すことが出来ないため、き裂先端の要素分割を全体に比べると細かくした。き裂長さ5(mm)において,節点 697 個、要素数は 1274 個である。

図 4-3.4 三角形要素分割11

第４節 解析結果及び考察
き裂長さ 5・6・8・10(mm)の場合における特異要素および三角形要素を用いた解析の結果をそれぞれ図 4-4.1〜4-4.16に示す。
図4-4.1〜4-4.4 はa=5(mm)、
図4-4.5〜4-4.8 はa=６(mm)、
図4-4.9〜4-4.12 はa=8(mm)、
図 4-4.13〜4-4.16 は a=10(mm)における結果を示している。それぞれの要素数及び節点数を表 4-4.1 及び4-4.2 に示す。

表 4-1.1 特異要素
き裂長さ (mm)

節点数 （個）56810

要素数 （個）722751751722217
226226217

表 4-1.2 三角形要素
き裂長さ (mm)

節点数 （個）56810

要素数 （個）69751753863912749319771164

図 4-4.1 および図 4-4.2 は、き裂長さ a=5(mm)における応力分布図である。垂直応力成分σ y はき裂先端近傍において急激に上昇していることがわかる。このσ y から K を算出する方法は以下のとおりである。
(2-1.1)式の第 2 式でθ = 0 とする。

σy =KⅠ(4-4.1)

2πr

よって、

KⅠ = σ y 2πr(4-4.2)

解析結果でσ y を求めて位置 r とσ y より KⅠを(4-4.2)式より計算する。より精度を高めるため K-r をプロットし、r→0 に外挿した点の K 値を対象き裂の K とする。図 4-4.1 の結果から K-r 関係を示したものが図 4-4.3である。図 4-4.3 および図 4-4.4 は K-r 関係を示している。特異要素において、K-r 関係はほぼ線形である。しかし、三角形要素については、若干のバラツキが見られる。この傾向はき裂長さが長くなるほど顕著に表れている。

126000

σy (MPa)

5000
4000
3000
2000
1000
000.002

0.004

0.006

0.008

0.01

き裂先端からの距離 (mm)
図 4-4.1

き裂長さ 5（ｍｍ）における応力分布

140

σy (MPa)

120
100
80
60
402002468

き裂先端からの距離 （ｍｍ）
図 4-4.2 三角形要素におけるき裂長さ 5（ｍｍ）の応力分布
131012

3.463.45

K (MPa・m 1/2)

3.44
3.43
3.42
3.41
3.4
3.39

0.002

0.004

0.006

0.008

0.01

き裂先端からの距離 (mm)
図 4-4.3 特異要素によるき裂長さ 5(mm)における応力拡大係数

8K (MPa・m1/2)

7654320246810

き裂先端からの距離（ｍｍ）
図 4-4.4 三角形要素によるき裂長さ 5(mm)における応力拡大係数

1412

6000

4000
3000

σy (MPa)

5000

2000
1000
000.002

0.004

0.006

0.008

0.01

き裂先端からの距離 (mm)
図 4-4.5 特異要素によるき裂長さ 6(mm)における応力分布100

σ y (MPa)

8060

40

2002468

き裂先端からの距離 （ｍｍ）
図 4-4.6 三角形要素によるき裂長さ 6(mm)における応力分布
151012

3.82

3.8
3.79
3.78
3.77
3.76

0.002

0.004

0.006

0.008

0.01

き裂先端からの距離 (mm)
図 4-4.7 特異要素によるき裂長さ 6(mm)における応力拡大係数

76K (MPa・m1/2 )K (MPa・m1/2 )3.81

5430246810

き裂先端からの距離 (mm)
図 4-4.8 三角形要素におけるき裂長さ 6(mm)における応力拡大係数
1612

7000
6000

σy (MPa)

5000
4000
3000
2000
1000
000.002

0.004

0.006

0.008

0.01

き裂先端からの距離 (mm)
図 4-4.9 特異要素によるき裂先端 8(mm)における応力分布9080

σy (MPa)

70
60
50
4030
2002468

き裂先端からの距離 (mm)
図 4-4.10 三角形要素によるき裂長さ 8(mm)における応力分布
1710

K (MPa・m1/2 )

4.47
4.46

4.45
4.44
4.43

0.002

0.004

0.006

0.008

0.01

き裂先端からの距離 (mm)
図 4-4.11 特異要素によるき裂長さ 8(mm)おける応力拡大係数

7K (MPa・m1/2)

6.5
65.554.5
402468

き裂先端からの距離 (mm)
図 4-4.12 三角形要素によるき裂長さ 8(mm)における応力拡大係数
1810

8000

4000

σy (MPa)

6000

2000

000.002

0.004

0.006

0.008

0.01

き裂先端から距離 (mm)
図 4-4.13 特異要素によるき裂長さ 10(mm)における応力分布
図 4-4.14 三角形要素によるき裂長さ 10(mm)における応力分布

250

σy (MPa)

200
150
100
50
0

12345

き裂先端からの距離 (mm)19675.32

K (MPa・m1/2 )

5.31
5.3
5.29
5.28
5.27
5.26
5.25

0.002

0.004

0.006

0.0080.01

き裂先端からの距離 (mm)
図 4-4.15 特異要素によるき裂長さ 10(mm)における応力拡大係数

6.6K (MPa・m 1/2)

6.46.265.85.65.45.2012345

き裂先端からの距離 (mm)
図 4-4.16 三角形要素によりき裂長さ 10(mm)における応力拡大係数206図 4-4.1〜4-4.16 に示されているように応力分布および応力拡大係数は特異要素と三角形要素では若干の違いが見られる。特異要素の応力分布はき裂先端に対してごく微小な領域であるためき裂先端近傍において応力集中を顕著に反映している。逆に三角形要素はき裂からの距離が特異要素に比べ大きいため、き裂先端での応力値は特異要素に比べ低い値となっている。また、値のばらつきの原因になっていると思われる。以上のことより、応力値より求められる応力拡大係数についても同様のことが考えられる。特異要素は三角形要素にくらべ精度がよくき裂先端近傍の応力場解析に適していることが分かる。
応力場解析より選られた特異要素および三角形要素の応力拡大係数と一般に使われている式(4-4.1)の近似式 Secant 公式の解との比較を図 4-4.17 で行った。

KⅠ = σ πa secπa(4-4.1)

W

ここで、a はき裂長さ、W は試験片幅である。この図から、3 つの方法による応力拡大係数の値がほぼ同等であることが分かる。若干き裂長さが短い時、特異要素については K 値が高くなる傾向にある。しかしながら、三角形要素については K 値を求める場合、K-r 関係にばらつきが見られた。これに対して、特異要素を用いた場合、き裂先端における特異性が保証されているため K-r 関係のばらつきは小さい。したがって、特異要素の場合、三角形要素に比べ、要素分割数が少なくても同等以上の結果を得ること出来ると思われる。

1614Secant公式
特異要素
三角形要素

K (MPa・m 1/2)

1210864200246810121416

き裂長さ (mm)
図 4-4.17 応力拡大係数の比較
図 4-4.18・図 4-4.19 は 5(mm)でのそれぞれの要素、また図 4-4.20 は式(2-2.1)より得られたき裂先端における応力場を示した。図 4-4.18 に示されている特異要素による結果では、数値積分点でのσ y の値から等高線を求めた結果である。三角形要素は三角形の面積全体に一様な応力として計算されているため、図 4-4.19 に示されている三角形要素は要素の重心を点として応力場が示されている。いずれの場合もσ y の広がりはほぼ 60 度であることが分かる。

21530353022030012D08)mm(5 53 3pit05 6k0c5ar0c206E41mor4fecnatsi147 72 2)
aPM(
5331313302200102D8)
mm
(pit
336kc
3 7a5 4r0c046mor
4f7.
e6cnatsi

1

222303022 5
3 30 03 33

図 4-4.19 三角形要素 5(mm)の応力分布

0323)
aP
M
(4)
mm
(ecafruskcarcmor
fecnatsi
D025)
mm(ecafruskcarcmor
fecnatsiD

図 4-4.18 特異要素き裂長さ 5(mm)の応力分布5)
aPM
(
407.
0301.
0001

図 4-4.20 Sceont 公式による近似応力分布8)mm(02Dpit36k.
1c9a6.
r.
11cmor4fecnatsi

12300.
104.
02)
mm(ecafruskcarcmor
fecnatsiD

第６章

結言

本研究では、特異要素を使用した八節点アイソパラメトリック要素を用いた中央切欠き試験片におけるき裂先端応力場の解析を行った。
・Secant 公式と解析結果を比較しても有用性が証明された。
・き裂先端においての r−1 / 2の特異性を示すことが出来るため、微小なき裂先端近傍の応力集中を反映する ことができる。
・特異性を考慮しない三角形要素を用いた場合に比べ、要素数が少なくても精度よく解析を行うことができる。
以上の点より特異要素を用いたアイソパラメトリック要素による有限要素法はき裂先端近傍における応力場を解析するのに適していることが証明できた。

参考文献（１）材料強度学入門
村上 理一・高尾 健一・萩山 博之(1995)（２）数値破壊力学白鳥 正樹・三好 俊郎・松下 久雄(1980)（３）計算力学ハンドブック 有限要素法・構造編日本機械学会（４）C 言語による有限要素法入門
森 博嗣 (1989)
森北出版（５）http://augustine.mit.edu/Christophe.Prudhomme/gfem.php

24

西日本法規出版
実教出版

謝辞
本研究を遂行するにあたり、終始懇切丁寧な指導を賜りました楠川量啓助教授に対し、深く感謝いたします。
また、高知工科大学材料強度研究室の皆さんからは多大なるご援助いただき、あわせて感謝いたします25