任意の構成関係を有する材料の応力解析に関する研究

(1)

任意の構成関係を有する材料の応力解析に関する研究

2015 年 12 月

長崎大学大学院工学研究科

伊野拓一郎

(2)

主要記号

各記号の下付き添字，i, j, k, l, m, n, o, p, ξ, η はテンソル記号のダミーインデックスであり，x, y, z, r, θ などの成分を表す．

P : 注目点

Q : 着力点

Γ : 境界

ΩP : 塑性領域

Ω_I : 介在物領域

σij : 応力

εij : ひずみ

ε^p_ij : 塑性ひずみ dε^p_ij : 塑性ひずみ増分

u, v : 変位

Sij : 偏差応力

Γ_ij : 力対の大きさ

σ^∞ : 無限遠方から作用する外力 σ₀ : 有限板の縁に作用する外力 σij(P) : 注目点P の応力

dσij(P) : 注目点P の応力増分

σ^∗_ij(P, Q_η) : 着力点Qにη方向の単位大きさの集中力が作用した時に注目点P に生じる応力

σ_ij^†(P, Q_ξη) : 着力点Qにξη方向の単位大きさの集中力対が作用した時に注目点P に生じる応力

φ_η(Q) : 着力点Qのη 方向の体積力の密度 Tij(Q) : 着力点Qに作用する体積力対 dTij(Q) : 着力点Qに作用する体積力対増分

D_ijkl : 弾性定数テンソル

Cijkl : 弾性コンプライアンステンソル E_ijkl : 単位行列

(5)

z(=x+iy) : 複素平面上での注目点の座標 z0(=ξ+iη) : 複素平面上での着力点の座標

Ω(z), ω(z) : Goursat の複素応力関数 F(=Fx+iFy) : 集中体積力の複素数表記

ρ(=ρ_x+iρ_y) : 体積力の密度の複素数表記

ν : Poisson 比

κ : Kolosov の定数

E : 縦弾性係数，Young 率 G : 横弾性係数

σ_Y : 降伏応力

(6)

第 1 章緒言

1.1 研究背景

有限要素法をはじめとして，現在多くの離散化解析手法が開発されている．離散化解析手法は大きく分けると、領域型と境界型に分けられる．領域型解法の代表例として，有限要素法，有限体積法，有限差分法などがある．これらの解析手法の特徴としては，領域全体を離散化し各離散化要素内を近似関数により近似して，支配微分方程式を解く方法である．特に有限要素法は，1956年にM.J.Turner らが梁の問題に用いてから，さまざまな改良や発展が行われてきた．¹⁾ しかし，領域型解法には本質的な問題として，領域全体を離散化するため節点番号や節点座標などの入力データを作成する作業が必要となり計算コストを増大させている．一方，境界型解法は，境界条件を基に積分方程式を作成し，境界の離散化解析により解を求めるため領域全体の離散化が必要ないという特徴があり，必然的に領域型解法と比較して計算コストが抑えられる．この境界型解析方法には，境界要素法や体積力法などがある．境界要素法は Betti の相反定理に基づいて積分方程式を構成する．²⁾ 一方，体積力法³⁾^∼²⁶⁾は重ね合わせの原理に基づいて積分方程式を構成する．つまり，基本的な応力場の重ね合わせにより解を得るため，直感的に解を得ることができ，また種々の直感的な工夫が導入し易いため，高精度の解を得る事が出来る．

この方法は1967年に西谷弘信が開発した解法である．本章では，本質的な理解の為に，体積力法の発表の4年前に公開された，だ円孔列を有する無限板の引張りの研究について注目する．²⁷⁾^∼²⁸⁾ この論文では，Fig.1に示すような円孔列を有する無限板の問題を取り扱っている．この問題を西谷は，各円孔がそれぞれ，ある大きさの引張を受ける問題に置き換えて解を得ている．つまり，この問題の境界条件は各円孔の１点に生じる表面力を境界条件とし，各円孔に作用させる荷重を未知数としている．

(7)

Infinite Plate

Unknown

Reference Point

Fig.1 円孔列を有する無限板の引張り

この考え方は体積力法の考え方の基礎になっている．体積力法では，仮想境界を分割し，適当な大きさの体積力を境界に分布させる事で自由境界を表現している．例え

ば，Fig.2に示すような円孔を有する無限板の問題の場合，円孔仮想境界を分割し，各

境界に体積力を分布させて，各境界の参照点の表面力が境界条件を満たすように仮想境界に分布させる体積力の大きさを決定する．

Reference Point

Unknown Infinite Plate

Fig.2 円孔を有する無限板の引張り

どちらも，簡単な応力場を重ね合わせる事により目的の解を得ており，解を得るまでの過程が非常に直感的である．この問題以外にも，半円切欠きを有する丸棒の引張りの問題²⁹⁾ でも同様の考え方に基づいて数値解析を行っている．このような考え方が体積力法の原点になったと，その後，西谷は述懐している．³⁰⁾^∼³¹⁾

ここまでは，体積力法における長所を紹介したが，次に，体積力法の短所について述べる．体積力法に用いる基本解は弾性解であり，各注目点に生じる応力やひずみと外力は比例関係にある事をうまく利用して目的の解を得ている．そのため，応力-ひずみ関係がHooke の法則のみに依らないような問題を体積力法は不得意としている．

(8)

方法 ⁴³⁾^∼⁴⁶⁾ がある．力対とは等しい大きさを持ち反対方向の集中力の対で，き裂問

題^47)∼48)を解く際などによく用いられてきた．この力対によって微小領域に外力と

は独立した荷重を作用させる事が可能になり，それにより生じるひずみを非弾性ひずみとして取り扱う方法である．これにより取り扱う事が可能となる問題としては，以下の問題が挙げられる．

• ^{弾塑性問題}

• 熱応力問題

• ^{粘弾性問題}

• ^{傾斜機能材料の問題}

これらの問題は，物理的現象は全く異なるが非弾性ひずみの大きさを決定する構成式が得られれば，陳らの方法により同じ方法で解を得ることができる．本研究ではこの事に注目して，上記の問題のうち弾塑性問題と傾斜機能材料の問題に体積力法を適用し，その高度化を行った．

(9)

1.2 研究目的

1.2.1 傾斜機能材料の問題

傾斜機能材料はセラミックから金属へと連続的に組成を変化させる事によって熱応力を緩和させる事を目的として開発された材料である．開発当初の仕様用途としては，宇宙往還機(スペースプレーン)を想定していたが⁴⁹⁾，現在では，核融合炉やエンジンまた，切削工具⁵⁰⁾の材料として用いられている．作成方法としては，ガス⁵¹⁾ や放電プラズマ⁵²⁾ による焼結法や，遠心力を利用した方法⁵³⁾ などがある．このように，さまざまな製造方法が提案されており，仕様用途に合わせて自由に材料の組成設計を行うことが出来るため，傾斜機能材料の設計式を作成する事は不可能である．

そこで，対象の材料組成に合わせて数値解析を行うプログラムの開発が求められる．

先行研究を挙げると，八田正俊らが介在物の問題として無限板中に一様な材料定数の介在物がある問題を取り扱っている．⁵⁴⁾ 彼らは，介在物と母材の間に体積力を分布させる事により解を得ているため，少ない方程式系で目的の解を得る事が出来る反面，本研究で取り扱うような材料定数が空間座標の関数として変化するような材料に関しては同様の解析方法を適用出来ない．そこで，介在物領域の材料定数と母材の材料定数から力対の大きさを表す式を構成し，陳らの方法を適応させた．本式の検証の為に，無限板中に円形の介在物がある問題を解くことで検証を行った．この検証の結果，無限板中の材料定数が任意に変化するような問題についても体積力法を用いて応力解析を行うことが可能になった．

1.2.2 弾塑性問題

応力集中部に弾性限度内の応力が生じる場合や，応力集中部に生じる塑性域が構造物の寸法に対して小さい小規模降伏条件にあるときは，弾性解析により，比較的容易にかつ高精度に現象を予測する事が出来る．そのため，機械構造物の設計を行う際に，弾性解析を行い，応力集中部の弾性応力状態について検討する事はもちろん重要であるが，応力集中部の応力が弾性限度を越えた状態についても，構造の強度を検討する塑性設計を行う事も重要である．その為，弾塑性問題は古くから多くの研究がなされている例えば，山田嘉昭は日本における有限要素法の第一人者であり，特に非

(10)

し有限要素法の拡張を行った．³²⁾^∼³³⁾ また，粘弾性問題についても有限要素法で取り扱う事を可能とした．³⁴⁾ そしてこの方法を応用して，塑性領域を有する切欠き底近傍の応力集中問題を松井が取り扱った．³⁵⁾ 有限要素法以外の解析的手法としては，

弾塑性破壊力学パラメータであるRice のJ積分を用いる方法を久保司郎が研究をした．³⁶⁾^∼⁴¹⁾ 実験的手法では，河田によって皮膜法による光弾塑性解析が提案された．

42) 体積力法においては，陳らが力対の分布により非弾性ひずみを表現する前節の方法を繰り返し行い，塑性ひずみを表現する．この点が傾斜機能材料の問題との大きく異なる点である．傾斜機能材料の問題の場合，局所的に材料定数が異なる線形問題であるが，弾塑性問題の場合，塑性ひずみが荷重履歴に依存するため，応力-ひずみ関係が非線形である．このことから，弾塑性問題は次に示すような特徴がある．

(i) 塑性領域が降伏条件によって決定される (ii) 塑性ひずみの大きさが荷重履歴に依存する

体積力法では，この特徴を表現するにあたり，塑性領域と思われる領域を離散化し塑性ひずみに相当する力対を内挿して解析を行う．そのため，解析対象の全領域を離散化する必要がなく，力対が埋め込まれた弾性問題として解析を行う事が出来る．これらの長所を有する為，体積力法による弾塑性解析を研究する事は有意義である．そこで，陳らの示した弾塑性解析の手法と本研究で行った手法の相違点及び改善点を示す．

(i) 塑性ひずみ増分をPrandtl-Reuss の式により決定した．

(ii) 力対を内挿するための離散化をDelaunay 分割により行った．

(iii) 力対を決定する際の解法の効率化を行った．

この (i), (ii), (iii)は一見別々の課題のようであるが，密接に関係している．まず，

“(i) 塑性ひずみ増分をPrandtl-Reuss の式により決定した．”についてだが，陳らは加工硬化指数H^′ に基づいて，塑性マトリックスS_ijkl を構成し応力増分から塑性ひずみ増分を決定し弾塑性解析を行った ⁵⁹⁾^∼⁶¹⁾．それに対して，本手法は塑性ひずみ増分を塑性構成式と降伏条件に基づいて決定し解析を行った．この塑性構成式は，塑性ひずみ増分が偏差応力と比例関係にあることを表しており，加工硬化指数に依らず任意の弾塑性体について解析を行う事が出来る反面，降伏条件が非線形方程式であるため非線形連立方程式により各要素の塑性ひずみの大きさを決定しなければならない．しかも，“(ii)力対を内挿するための離散化をDelaunay 分割により行った．”に示すように本研究ではDelaunay 自動分割により予め入力データとして要素データを用意する事無く多くの要素を生成する事が可能となったため，この非線形連立方程式

(11)

の規模は莫大となった．そこで，“(iii) 力対を決定する際の解法の効率化を行った．” の研究に取り組んだ．非線形方程式は各要素の代表点の応力により構成されるため，

各要素に単位大きさの力対が内挿され，また，自由表面の境界条件を満たした各要素の応力成分を予め計算した．この事により，非線形方程式を構成する代表点の応力を線形和により構成する事が可能になった．また，この事を様々な形状の応力集中問題について検証を行った．

(12)

1.3 論文の構成

本章では，体積力法の歴史および特徴また，本研究で取り扱う問題についての目的について述べた．第 2章では弾性解析について取り扱った．本研究で取り扱う任意の応力-ひずみ関係を有するような材料でも，すべて力対を埋め込んだ弾性解析に置き換えて解析を行う為，弾性解析は本研究において非常に重要である．その為，他の章に比べより詳細に述べた．第3章については，力対の埋め込みによる方法を用いた傾斜機能材料の問題の解法について述べた．第4章では力対の埋め込みによる方法を繰り返し用いた弾塑性問題の解法について述べた．第5章は本研究を総括し今後の展開について述べた．以下に各章各節の構成を箇条書きで示した．

1.3.1 第 1 章

第1節

• ^{体積力法の位置づけ}

• 体積力法が開発される前に考えだされた解析手法と体積力法の類似点の考察

• 非線形問題の数値解析に関する取り組み第2節

• 傾斜機能材料を数値解析する事の目的と研究方針

• 陳らが定義した弾塑性問題の解法を汎用化させる事の目的と研究方針第3節

• ^{論文の構成を記載}

1.3.2 第 2 章

第1節

• ^第2章の要約第2節

• 体積力法で用いられる各記号について

• 弾性解析の支配積分方程式第3節

(13)

• Kelvin の解について

• Goursat の複素応力関数について第4節

• ^{体積力の密度について}

• 座標変換を用いた任意直線境界に分布する体積力による応力場を表す複素応力関数

• 引張りを受ける正方形矩形板の弾性解析第5節

• 任意の円弧境界に分布する体積力による応力場を表す複素応力関数

• 引張りを受ける円孔を有する無限板の解析

• 内圧を受ける円孔を有する無限板の解析第6節

• 有限板の解析を行う際に生じる剛体変位の除去について第7節

• 引張りを受ける円孔を有する有限板の解析

• ^{切欠きを有する帯板の}3点曲げの解析第8節

• ^第2章のまとめ

1.3.3 第 3 章

第1節

• 介在物領域内部と介在物領域外部の支配積分方程式

• 引張りを受ける円形介在物を有する無限板の解析

• 引張りを受ける円形介在物と円孔を有する無限板の解析

• 引張りを受ける円形介在物とき裂を有する無限板の解析第3節

• ^第3章のまとめ

(14)

1.3.4 第 4 章

第1節

• 力対埋め込みモデルの概念について

• 力対の定義および，集中力対の基本解

• 弾塑性解析の支配積分方程式第3節

• Delaunay 分割を援用した塑性領域の離散化について

• 塑性領域の予測について第4節

• 塑性ひずみの大きさの効率的な決定方法について第5節

• Nadai の弾塑性解

• 内圧を受ける円孔を有する無限板の弾塑性解析第6節

• 線形切欠き力学について

• ^{切欠きを有する帯板の}3点曲げ解析第7節

• ^第4節のまとめ

1.3.5 第 5 章

• 本研究により得られた結論

• ^{今後の本研究の展望}

(15)

第 2 章体積力法による弾性解析

⁶⁸⁾^∼⁷⁴⁾

2.1 序論

本章では，体積力法による弾性解析³⁾^∼⁴⁾の方法を紹介する．まずは，弾性解析における基礎積分方程式を示す．本論文では積分方程式を解く際に解がない場合を除いて数値積分を用いず解析解を用いているので，直線や円境界の微小区間に体積力が一様に分布した時の解を組み合わせて様々な問題を解いている．そこで，直線境界や円境界に体積力が分布した時の基本解を示し，その基本解を用いて，厳密解がすでに得られている問題を解いて解析精度や分割数による影響を議論する．

2.2 弾性解析における基礎方程式

体積力法は重ね合わせの原理に基づいた境界型の弾性応力解析法であるため，注目点の応力は積分項を含む和で表される．例えば，Fig.3のような円孔を有する無限板を引張る問題を考える．

Infinite Plate

Imaginaly Boundary

σ

^∞

σ

^∞

Γ Q

P

σ

_ij

(P )

Fig.3 引張を受ける円孔を有する無限板

(16)

この時，実際に円孔を有する無限板の解を用いるのではなく，円孔と同じ大きさの仮想境界に体積力を分布させて解を得る．体積力法では，境界の単位長さに作用する体積力の大きさを体積力の密度という．この密度関数をφη(Q)と表す．Qは着力点を表し，円孔を有する無限板の問題では円境界上の点である．η は分布する体積力の方向である．基本解としては集中力の解を用いる．集中力の解^*1はσ^∗_ij(P, Q_η)と表される．着力点Qにη方向の集中力が作用した時に注目点P に作用するij方向の応力を意味している．そして，外部荷重をσ^∞，注目点P のij 方向の応力をσ_ij(P)，円孔となるべき仮想境界をΓと定義する．この時，基礎方程式は次のようになる．

σ_ij(P) =

∫

Γ

σ_ij^∗(P, Q_η)φ_η(Q)dΓ +σ^∞ (1) 本章で行った弾性解析はこの積分方程式に基づいて解析を行った．

*1集中力の解の詳細については次節にて示す．

(17)

2.3 基本解

無限板の一点に集中力が作用する際に生じる応力場の解をKelvin の解³⁾^∼^4), ⁶⁴⁾^∼

78) という．Fig.4は注目点と着力点の関係及び，記号の意味を示している．

Infinite Plate

Reference Point

Source Point

x, ξ y, η

z0 = ξ+iη

z = x+iy

F x F y

Fig.4 着力点と注目点の定義

Fx もしくは，Fy が着力点z0 に作用した時の，注目点 z の変位u, v および応力 σ_xx, σ_yy, τ_xy を以下に示す．

u^F^x(x, y) = 1 2πG(κ+ 1)

[ κln 1

r₁ + A² A²+m²

]

Fx (2)

v^F^x(x, y) = 1 2πG(κ+ 1)

[ Am A²+m²

]

Fx (3)

σ_xx^F^x(x, y) = A

2π(κ+ 1)y(A²+m²)²

[κ(A²+m²) + (3A²−m²)]

Fx (4) σ_yy^F^x(x, y) = −A

2π(κ+ 1)y(A²+m²)²

[κ(A²+m²) + (−A²−5m²)]

Fx (5) τ_xy^F^x(x, y) = m

2π(κ+ 1)y(A²+m²)²

[κ(A²+m²) + (3A²−m²)]

Fx (6) u^F^y(x, y) = 1

2πG(κ+ 1)

[ Am A²+m²

]

Fy (7)

(18)

v^F^y(x, y) = 1 2πG(κ+ 1)

[ κln 1

r1

+ m² A²+m²

]

Fy (8)

σ^F_xx^y(x, y) = −m

2π(κ+ 1)y(A²+m²)²

[κ(A²+m²) + (−5A²−m²)]

Fy (9) σ^F_yy^y(x, y) = m

2π(κ+ 1)y(A²+m²)²

[κ(A²+m²) + (−A²+ 3m²)]

F_y (10) τ_xy^F^y(x, y) = A

2π(κ+ 1)y(A²+m²)²

[κ(A²+m²) + (−A²+ 3m²)]

Fy (11) 上付き添字がFx のものは着力点にFx が作用した場合，Fy のものは着力点にFy が作用した場合をあらわしている．また，A, m, r₁ は以下のように定義する．

A= ξ−x

y (12)

m= η−y

y (13)

r1 =√

(ξ−x)²+ (η−y)² (14)

κはKolosov の定数で，各応力状態により以下のように定義する．

κ= 3−ν

1 +ν (平面応力状態) (15) κ= 3−4ν (平面ひずみ状態) (16) ν はポアソン比である．

式 (2)∼^式(11) は各変位成分，各応力成分ごとに表しているが，二次元問題では，

Goursat の複素応力関数^64)∼^{78) *2}を用いる事で取り扱う式が大幅に少なくなる．

Kelvin の解のGoursat の複素応力関数表示は以下のようになる．

Ω(z) =− F

2π(κ+ 1)ln(z−z0) (17)

ω(z) = κF

2π(κ+ 1)ln(z−z0) + F 2π(κ+ 1)

z0

z−z0

(18) F は体積力の成分を複素数で表記したものでF =Fx+iFy のように定義する．複素応力関数と各応力成分と変位成分は以下の関係式で示される．

σxx+σyy = 2 (dΩ

dz + dΩ dz

)

(19) σyy−σxx+ 2iσxy = 2

( zd²Ω

dz² + dω dz

)

(20)

*2Goursatの複素応力関数は付録にて詳細を示す．

(19)

2G(u+iv) =κΩ−zdΩ

dz −ω (21)

(20)

2.4 直線境界に体積力が分布するときの基本解

体積力法は境界型の数値解析法であるため，境界を離散化し解析を行う．その時，

各境界作用させる体積力の単位長さあたり大きさは一定もしくは重み関数の係数を一定として解析を行う．この各区間の単位長さあたりの体積力の大きさを体積力の密度 ρという．本論文上では体積力の密度も複素変数として取り扱うため，x方向の体積力の密度ρ_x，y方向の体積力の密度ρ_y としたとき体積力の密度ρは

ρ=ρ_x+iρ_y (22)

と表し，微小区間dsに作用する体積力の大きさF(=F_x+iF_y)は

F =ρ ds (23)

となる．

本節では，直線境界を表現する際に用いる，直線境界に体積力が一様に分布した時の解を示す．任意の直線境界を積分するよりも，x軸上に分布する体積力による応力を計算し，応力の座標変換を行う方が計算が容易であるのでFig.5のようなローカル座標系(x^′, ξ^′-y^′, η^′ 座標)で積分⁷⁹⁾を行う．グローバル座標系(x, ξ-y, η座標)での着力区間を[z01(=ξ1+iη1) :z02(=ξ2+iη2)]，着力区間の長さの半長を2a，着力区間の中心をz_c，グローバル座標系とローカル座標系のなす角度をθ とする．

Table.1 座標変換

Before coordinate transformation

After coordinate transformation

Reference point z z^′ = (z−zc)e⁻^iθ

Initial point of source section z01 z₀₁^′ =−a End point of source section z₀₂ z₀₂^′ =a

Body Force F F^′ =F e⁻^iθ

∫ a

−a

dΩ =

∫ a

−a

− F^′

2π(κ+ 1)ln(z^′−z₀^′)dξ^′

=

∫ a

−a− F^′

2π(κ+ 1)ln (z^′−ξ^′)dξ^′

(21)

Infinite Plate

Global coordinate

Local coordinate

2a a

z01

z02

θ

zc

y, η

x, ξ y^′, η^′

x^′, ξ^′ z0(x, y)→z^′0(x^′, y^′)

Fig.5 任意直線上に分布する体積力

= F^′

2π(κ+ 1) {

(z^′−a) ln (z^′−a)−(z^′+a) ln (z^′+a) + 2a }

(24) d

dz^′

∫ a

−a

dΩ = F^′ 2π(κ+ 1)ln

(z^′−a z^′+a

)

(25) d²

dz^′²

∫ a

−a

dΩ = F^′ 2π(κ+ 1)

2a

z^′²−a² (26)

∫ a

−a

dω =

∫ a

−a

κF^′

2π(κ+ 1)ln(z^′−z₀^′)dξ^′+

∫ a

−a

F^′ 2π(κ+ 1)

z^′₀ z^′−z₀^′

=

∫ a

−a

κF^′

2π(κ+ 1)ln (z^′−ξ^′)dξ^′+

∫ a

−a

F^′ 2π(κ+ 1)

ξ^′ z^′−ξ^′

=− κF^′ 2π(κ+ 1)

{

(z^′−a) ln (z^′−a)−(z^′+a) ln (z^′+a) + 2a }

− F^′ 2π(κ+ 1)

{ z^′ln

(z^′−a z^′+a

) + 2a

}

(27) d

dz

∫ a

−a

dω =− κF^′ 2π(κ+ 1)ln

(z^′−a z^′+a

)

− F^′ 2π(κ+ 1)

{ 2az^′

z^′²−a² + ln

(z^′−a z^′+a

)}

(28)

(22)

目点 zの応力を求める事が出来る．この事を利用して，Fig.6のような単純な正方形板を引張る問題を体積力法によって解く．境界条件は，要素の中心の表面力で満たすように体積力の大きさを決定する．また，Poisson 比は0として計算を行った．

Plane plate

parts Reference point

(center of this plate)

Divided into

σ

⁰

σ

⁰

N

Fig.6 一様引張をうける正方形板

注目点は正方形板の中心としy方向応力σy は，当然，荷重σ0と等しくなければならない．そこで，σ_y/σ₀ と分割数について，Fig.7に示す．Table.2は解析結果を数値で示している．横軸を1/N としており，横軸が0の時の値にσy/σ0が1に限りなく近づいているため，直線境界に体積力を分布させる事により正しく有限板を表している事を示している．また，Table.2より，32分割(1辺あたり8分割)程度で，解析誤差は1%未満になる事が分かった．

(23)

0.975 0.980 0.985 0.990 0.995 1.000 1.005 1.010 1.015 1.020 1.025

0.000 0.050 0.100 0.150 0.200 0.250 Exact solution

σy/σ⁰

1/N

Fig.7 正方形板中央のy方向応力と分割数N

Table.2 正方形板中央のy方向応力と分割数N(Fig.7の参照)

N σ_y/σ₀ N σ_y/σ₀ N σ_y/σ₀ N σ_y/σ₀ 4 0.976473 32 1.009765 240 1.002839 1600 1.000974 8 1.020028 36 1.009054 280 1.002598 2000 1.000859 12 1.017171 40 1.008461 320 1.002406 2400 1.000776 16 1.014933 80 1.005465 360 1.002250 2800 1.000712 20 1.013124 120 1.004270 400 1.002119 3200 1.000661 24 1.011728 160 1.003598 800 1.001434 3600 1.000619 28 1.010637 200 1.003156 1200 1.001143 4000 1.000583

(24)

2.5 円弧境界に体積力が分布するときの基本解

体積力法で切欠きなどを表現する際に，円境界にそって体積力を分布させる事により表現する．前節では，要素内に分布する体積力は一定として解析を行ったが，問題の特徴を取り入れた重み関数を取り入れる事で解析精度を向上させる事が出来る．この重み関数を体積力法では基本密度関数φi という．基本密度関数を用いる例として，

円孔を有する無限板の問題である．この問題の場合，体積力の分布が境界の投影長さに比例するという特徴がある．投影長さとは，Fig.8のdy, −dxに相当し，x方向の体積力はdyにy方向の体積力には−dxに対応する．

Infinity Plate

Projection length a

σ^∞ σ^∞

−dx dy

θ

Fig.8 一様引張りを受ける円孔を有する無限板

そこで，体積力と体積力の密度は次のように表される．

dFx =ρxdy (29)

dFy =ρy(−dx) (30)

(25)

また，dy, −dxはそれぞれ次のように表される．

x=acosθ ⇒ dx=−asinθdθ (31) y =asinθ ⇒ dy =acosθdθ (32) 式(29), (30)に式(31), (32)を代入すると，

dFx =ρxacosθdθ=aρxφxdθ (33) dF_y =ρ_yasinθdθ=aρ_yφ_ydθ (34) よって円孔を有する無限板の引張の問題の基本密度関数φ_iが導かれた．この基本密度関数を考慮した基本解をGoursaの複素応力関数表示で示す．体積力が分布する範囲l は，半径r の境界上の角度αからβまでの範囲である．

∫

l

dΩ(z) =

∫

l

− dF

2π(κ+ 1)ln (z−z₀)

=

∫ β α

−a(ρxcosθ+iρysinθ)

2π(κ+ 1) ln(z−ae^iθ)dθ

=−a(ρx+ρy) 4π(κ+ 1)

[i a

(z−ae^iθ) { ln(

z −ae^iθ)

−1}]^β

α

−a(ρx−ρy) 4π(κ+ 1)

[ iln(

z−ae^iθ) e^iθ + ia

z {

ln

( e^iθ z−ae^iθ

)}]^β

α

(35) d

dz

∫

l

dΩ(z) =−a(ρ_x+ρ_y) 4π(κ+ 1)

[i aln(

z−ae^iθ)]^β

α

−a(ρ_x−ρ_y) 4π(κ+ 1)

[ i

ze^iθ − ia z² ln

)]β α

(36) d²

dz²

∫

l

dΩ(z) =−a(ρx+ρy) 4π(κ+ 1)

[ i a(z−ae^iθ)

]β

α

−a(ρx−ρy) 4π(κ+ 1)

[ i

(a z − 1

e^iθ

) 1 (z−ae^iθ)² +2ia

z³ ln e^iθ

z−ae^iθ

+ 2ia z²(z−ae^iθ)

]^β

α

(37)

∫

l

dω(z) =

∫

l

κdF

2π(κ+ 1)ln(z−z0) +

∫

l

dF 2π(κ+ 1)

z0

z−z₀

=

∫ β α

aκ(ρ_xcosθ−iρ_ysinθ)

2π(κ+ 1) ln(z−ae^iθ)dθ +

∫ β α

a(ρxcosθ+iρysinθ) 2π(κ+ 1)

ae⁻^iθ z−ae^iθdθ

(26)

= aκ(ρ_x−ρ_y) 4π(κ+ 1)

[i a

(z−ae^iθ) { ln(

z−ae^iθ)

−1}]^β

α

+ aκ(ρx+ρy) 4π(κ+ 1)

[ iln(

z −ae^iθ) e^iθ + a

z ln

)]^β

α

+ a(ρx+ρy) 4π(κ+ 1)

[ia z ln

(z−ae^iθ e^iθ

)]β α

+ a(ρx−ρy) 4π(κ+ 1)

[ ia

2ze^2iθ + ia²

z²e^iθ + ia³ z³ ln

(z−ae^iθ e^iθ

)]β α

(38)

∂

∂z

∫

l

dω(z) = aκ(ρx−ρy) 4π(κ+ 1)

[i a ln(

z−ae^iθ)]

+ aκ(ρx+ρy) 4π(κ+ 1)

[ i

ze^iθ − ia z² ln

)]β α

+ a(ρ_x+ρ_y) 4π(κ+ 1)

[ a iz² ln

(z−ae^iθ e^iθ

)

− a

iz(z−ae^iθ) ]^β

α

+ a(ρx−ρy) 4π(κ+ 1)

[ a

2ie^2iθz² + 2a²

ie^iθz³ + 3a³ iz⁴ ln

(z−ae^iθ e^iθ

)

− a³ iz³(z−ae^iθ)

]β α

(39) この基本解を用いて円孔を有する無限板の引張り問題の解析を行った．解析条件は以下の通りである．

Table.3 解析条件

Load (σ^∞) 1.0

Poisson’s ratio (ν) 0.0

Number of division (N) 4

Radius of circular hole (a) 1.0

Reference point (z) (1.0, 0.0), (1.1, 0.0), (1.2, 0.0), · · · ,(2.0, 0.0) (0.0, 1.0), (0.0, 1.1), (0.0, 1.2), · · · ,(0.0, 2.0)

解析結果の検討にはG.Kirsh の解を参照した．⁸⁰⁾ 半径方向応力 σrr，周方向応力 σ_θθ，せん断応力τ_rθは以下の通りである．

σ_rr = σ^∞ 2

( 1− a²

r² )

− σ^∞ 2

(

1 + 3a⁴

r⁴ − 4a² r²

)

cos 2θ (40)

σθθ = σ 2

( 1 + a²

r² )

+ σ^∞ 2

(

1 + 3a⁴ r⁴

)

cos 2θ (41)

(27)

τrθ =−σ^∞ 2

(

1− 3a⁴

r⁴ + 2a² r²

)

sin 2θ (42)

この解と解析結果の比較をしたものを Fig.9∼10に示す．Fig.9は，x 軸上のy 方向応力を表しており，参照点での応力は厳密解と非常によく一致しており，有効数字6 桁まで一致している．Fig.10は，y軸上のx方向応力を表しており，Fig.9と同様に厳密解と非常よく一致している．

1.0 1.5 2.0 2.5 3.0

1.0 1.2 1.4 1.6 1.8 2.0

’s solution G.Kirsh

Analytic solution

Infinite plate

x/a σy/σ∞

σy

y

x

Fig.9 x軸上のy方向応力σy

(28)

-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2

1.0 1.2 1.4 1.6 1.8 2.0

’s solution G.Kirsh

Analytic solution

Infinite plate

x σx/σ∞

σx

y

y/a

Fig.10 y軸上のx方向応力σ_x

Table.4 x軸上のy方向応力(Fig.9を参照)

x/a Analytic solution G.Kirsh’s solution

1.0 2.999999 3.000000

1.1 2.437743 2.437743

1.2 2.070601 2.070602

1.3 1.821049 1.821050

1.4 1.645564 1.645564

1.5 1.518518 1.518519

1.6 1.424194 1.424194

1.7 1.352605 1.352606

1.8 1.297210 1.297211

1.9 1.253604 1.253605

2.0 1.218750 1.218750

(29)

Table.5 y軸上のx方向応力(Fig.10を参照)

y/a Analytic solution G.Kirsh’s solution

1.0 -1.000000 -1.000000

1.1 -0.611297 -0.611297

1.2 -0.376157 -0.376157

1.3 -0.229333 -0.229334

1.4 -0.135360 -0.135360

1.5 -0.074074 -0.074074

1.6 -0.033569 -0.033569

1.7 -0.006585 -0.006585

1.8 0.011431 0.011431

1.9 0.023403 0.023404

2.0 0.031250 0.031250

(30)

External load

Circular boundary

Distributing body force

External load

Fig.11 円境界に作用させる体積力の分布

この時の自由境界を表現するために円境界に分布させる体積力の分布をFig.11に示す．つまり，無限板中の仮想円境界の内側に Fig.11と同じ大きさ同じ方向の体積力を作用させることで，円境界外部に生じる応力は円弧がある場合と同じ応力を生じる．しかし，式(35)∼(39)は基本密度関数を含んでいるため，遠方から引張を受けるような場合には有効であるが，それ以外の場合には用いる事ができない．そこで，基本密度関数を含まない場合についても示す．この時，積分範囲等は式(35)∼(39)の場合と同じである．∫

ldΩ(z)と∫

ldω(z)は積分の性質上解析解を得る事ができない．^*3

*3 ∫ _β

α

ln(z−z₀)dθ=

∫ _β

α

ln(z−ae^iθ)dθ

(31)

しかし，式(19)∼(21)より，応力を計算する際にはΩの１階微分，と２階微分そして，ω の１階微分しか必要がないため，積分と微分の順序を入れ替えて計算を行う．

変位を計算する場合には，解析解が存在しない項のみを数値積分を行う事により解決する．

∫

l

dΩ(z) =

∫

l

− dF

2πi(κ+ 1)ln (z−z0)

=

∫ β α

−a(ρ_x+iρ_y)

2πi(κ+ 1) ln(z−ae^iθ)dθ (43)

d dz

∫

l

dΩ(z) =

∫ β α

−a(ρx+iρy) 2πi(κ+ 1)

1 z−ae^iθdθ

= a(ρ_x+iρ_y) 2πi(κ+ 1)

1 z

{ ln

(z−ae^iβ z−ae^iα

)

−i(β−α) }

(44) d²

dz²

∫

l

dΩ(z) = a(ρx+iρy) 2πi(κ+ 1)

1 z²

{ z

z −ae^iβ − z

z−aeîα + z(eîβ −eîα) eî(β+α)

−2aln

(z−ae^iβ z −ae^iα

)

+ 2ia(β−α) }

(45)

∫

l

dω(z) =

∫

l

κdF

2πi(κ+ 1)ln(z−z0) +

∫

l

dF 2πi(κ+ 1)

z0

z−z₀

=

∫ β α

aκ(ρ_x−iρ_y)

2πi(κ+ 1) ln(z−ae^iθ)dθ+

∫ β α

a(ρ_x+iρ_y) 2πi(κ+ 1)

ae^−iθ z−ae^iθdθ

=

∫ β α

aκ(ρx−iρy)

2πi(κ+ 1) ln(z−ae^iθ)dθ + a(ρx+iρy)

2πi(κ+ 1) a iz³

[ az

z−ae^iθ − z

e^iθ −2aln

(z−ae^iθ e^iθ

)]β α

(46) d

dz

∫

l

dω(z) =−aκ(ρx−ρy) 2πi(κ+ 1)

1 z

{ ln

(z−ae^iβ z−ae^iα

)

−i(β−α) }

− a²ρ 2πi(κ+ 1)z³

[ a²z(e^iβ −e^iα)

(z−aeîβ)(z−aeîα) + z(eîβ −eîα) eî(β+α)

−2aln

(z−re^iβ z −re^iα

)

+ 2ia(β−α) ]

(47)

t=e^iθとするとdt/dθ=ie^iθとなり置換積分を行うと

∫ β α

ln(z−z0)dθ=

∫ e^iβ e^iα

ln(z−at)dt it

= 1 i

∫ e^iβ e^iα

ln(z−at)

t dt

= ln(z−at) ln(−at

) +z−at

+(z−at)²

2 2 +(z−at)³

2 3 · · ·

(32)

式(43)∼(47)の基本解を用いる問題は多く考えられるが，最も容易に厳密解の比較が出来る問題として，Fig.12のような内圧をうける円孔の問題がある．

Inner Pressure Infinite Plate

σ^r(r), σ^θ(r), τ^rθ(r)

r

a p

Fig.12 内圧を受ける円孔を有する無限板

この問題の厳密解は次の円筒問題のAiryの応力関数⁸⁰⁾を解く事で得られる．

σrr(r) =C1+ C2

r² (48)

σθθ(r) =C1− C2

r² (49)

τrθ(r) = 0 (50)

この時，円孔半径はa，内圧はpであるため，円孔中心からr 離れた点の半径方向応力σr(r)，は境界条件(σr(a) =−p, σr(∞) = 0)から以下のようになる．

σ_rr(r) =−a²p

r² (51)

σ_θθ(r) = a²p

r² (52)

同様の問題を以下の解析条件に基づいて体積力法で解析をした．

(33)

Table.6 解析条件

Poisson’s ratio (ν) 0.0

Number of division 4

Stress state Plane stress

応力境界値問題としてこの問題を解き，表面力を境界条件とした．その時の半径方向応力 σr と周方向応力σθ をFig.13に示す．そして，Table 7に解析結果と厳密解の数値の比較を行った．

-1.0 -0.5 0.0 0.5 1.0

1.0 1.5 2.0 2.5 3.0

Numerical solution

Numerical solution Exact solution

Exact solution

r/a σr/p,σθ/p

σ^r

σ^θ

Fig.13 内圧を受ける円孔を有する無限板に生じる半径方向応力σ_rと周方向応力σ_θ

(34)

Table.7 内圧を受ける円孔を有する無限板に生じる半径方向応力σ_r と周方向応力σ_θ(Fig.13の参照)

r σ_r/p (Numeri- cal solution)

σ_r/p (Exact solution)

σ_θ/p (Numeri- cal solution)

σ_θ/p (Exact solution) 1.0 -1.000000 -1.000000 1.000000 1.000000 1.1 -0.826446 -0.826446 0.826446 0.826446 1.2 -0.694444 -0.694444 0.694444 0.694444 1.3 -0.591716 -0.591716 0.591716 0.591716 1.4 -0.510204 -0.510204 0.510204 0.510204 1.5 -0.444444 -0.444444 0.444444 0.444444 1.6 -0.390625 -0.390625 0.390625 0.390625 1.7 -0.346020 -0.346021 0.346020 0.346021 1.8 -0.308642 -0.308642 0.308642 0.308642 1.9 -0.277008 -0.277008 0.277008 0.277008 2.0 -0.250000 -0.250000 0.250000 0.250000 2.1 -0.226757 -0.226757 0.226757 0.226757 2.2 -0.206611 -0.206612 0.206611 0.206612 2.3 -0.189035 -0.189036 0.189035 0.189036 2.4 -0.173611 -0.173611 0.173611 0.173611 2.5 -0.160000 -0.160000 0.160000 0.160000 2.6 -0.147929 -0.147929 0.147929 0.147929 2.7 -0.137174 -0.137174 0.137174 0.137174 2.8 -0.127551 -0.127551 0.127551 0.127551 2.9 -0.118906 -0.118906 0.118906 0.118906 3.0 -0.111111 -0.111111 0.111111 0.111111

この場合でも，厳密解と解析解は有効数字５桁まで一致している事が分かる．式 (43)∼(47)を用いて，内圧を受ける円孔を有する無限板の問題を解くことが出来る事が示された．

(35)

2.6 剛体変位

本章では，有限板の問題についての取り扱いについて言及する．有限要素法において，解析条件として，拘束条件を正しく設定する事は非常に重要である．正しい拘束条件を指定していない場合，解析対象を正しく表現することができなくなり，場合によっては，剛体変位により連立方程式を解く事ができなくなる場合もある．この剛体変位に関しては，有限板の問題を取り扱う場合のみ，体積力法でも問題となる．境界条件として，変位を設定する場合は剛体変位を考慮しなくても良いが，応力や表面力を境界条件として設定する場合は剛体変位を考慮しなくてはならない．つまり，内部の応力状態が境界条件をみたしているなら，解析対象がどこに剛体移動していても問題は無いからである．図示すると次のようになる．

Finite

plate Finite

plate

Finite plate Same

shape

Same load Same stress field

Fig.14 同一荷重を受ける同一形状の３つの有限板

３つの同一荷重を受ける同一形状の有限板は内部に生じている応力は同じであるが，応力境界条件だけでは剛体変形に関する条件を満たす事ができない．では，なぜ

(36)

限遠方の変位が正則でなければならないという条件が含まれているからである．

そこで，本研究では有限板の問題を応力境界条件で解く場合について示す．応力境界条件には，変位に関する条件がないため，以下の条件を追加する．

(i) 全境界に作用する体積力の釣り合い

(ii) 任意の１点を基準とした時の，全境界に作用する体積力により生じるモーメントの釣り合い

境界の分割数n，i番目の境界の合力 F_xⁱ, F_yⁱ, 任意の基準点から各境界の中心までの距離のベクトルを⃗lⁱ(lⁱ_x, lⁱ_y)とする時，条件式は以下のようになる．

∑n i=1

F_xⁱ = 0 (53)

∑n i=1

F_yⁱ = 0 (54)

∑n i=1

lⁱ_xF_yⁱ−lⁱ_yF_xⁱ = 0 (55) 式(53)∼(54)は条件(i)に相当し，式(55)は条件(ii)に相当する．つまり，体積力の密度を決定する為に解く行列式の大きさは((総分割数)×2+3)となる．2.4節の正方形板の問題についても上記の条件を用いて解いた．

(37)

2.7 円孔を有する有限板の問題

本節では，さらにこの条件式を用いて円境界の解と直線境界の解の両方を用いて解く問題を取り扱う．まず，円孔を有する有限板の問題を取り扱う．解析モデルを Fig.15に示す．

Square plate

Circular hole

Load Load

a b

b

s

t

σ

^∞

σ

^∞

Fig.15 円孔を有する引張りを受ける有限板

この時，s=2a,t=2a,b=4aとすると，円孔縁に生じる応力集中係数K_t は6.389⁶⁵⁾ となる．直線境界と円弧境界の基本解は式 (24)∼(28)，式 (35)∼(39)を用いた．ま

た，Table.8の解析条件に基づいて解析を行った場合の解析結果を厳密解と境界の分

割数に注目して比較する．

任意の構成関係を有する材料の応力解析に関する 研究