2012年8月6日
有限要素法による２次元応力解析
— No-tension
材料の考慮
—
目次
1. 2次元弾性問題の剛性方程式 1 1.1 仮想仕事の原理による釣合式の表示 . . . 1 1.2 要素剛性方程式 . . . 2 2. ４節点パラメトリック要素 5 2.1 ひずみ-変位関係行列[B]の導入. . . 5 2.2 要素剛性行列および節点力ベクトル . . . 6 2.3 要素応力の算定 . . . 7 2.4 物体力による節点力 . . . 7 3. 2次元弾性応力解析のNo-tension解析への拡張 8 3.1 解析の基本手順 . . . 8 3.2 引張応力成分除去の方法 . . . 11 3.3 収束の考え方. . . 12 4. プログラミング 13 4.1 ２次元応力解析用subroutine . . . 134.2 FEMプログラム共通のsubroutineおよびfunction . . . 14
付録A 三角形定ひずみ要素 15 A.1 変位関数 . . . 15
A.2 ひずみ-変位関係行列. . . 16
A.3 要素剛性行列および節点力ベクトル . . . 16
1.
2
次元弾性問題の剛性方程式
1.1 仮想仕事の原理による釣合式の表示 座標系として，水平方向にx座標，水平方向に直交する上下方向にy座標を設定する． 2次元弾性体が釣合状態にある時，直交座標系x-yにおいて，以下の関係が成立している． 物体内部 ∂σx ∂x + ∂τxy ∂y + XB = 0 ∂σy ∂y + ∂τxy ∂x + YB= 0 物体表面 Sx= σx· ` + τxy· m Sy= σy· m + τxy· ` ここに， σx，σy，τxy ：内部応力の成分 XB，YB ：物体力 Sx，Sy ：表面力 `，m ：境界上法線の方向余弦 釣り合い状態にある場合は任意の仮想変位δuおよびδvに対し，次の仮想仕事の式が成り立つ． ∫ A [ δu ( ∂σx ∂x + ∂τxy ∂y + XB ) + δv ( ∂σy ∂y + ∂τxy ∂x + YB )] dA = 0 (1) 式(1)は，部分積分公式f · g0= (f· g)0− f0· gを用いて，式(2)の通り変形できる． 0 = ∫ A [ ∂(δu· σx) ∂x + ∂(δu· τxy) ∂y + ∂(δv· σy) ∂y + ∂(δv· τxy) ∂x ] dA − ∫ A [ ∂(δu) ∂x σx+ ∂(δu) ∂y τxy+ ∂(δv) ∂y σy+ ∂(δv) ∂x τxy ] dA + ∫ A (δu· XB+ δv· YB) dA (2) 式(2)右辺第1項は，Greenの公式 ∫∫ D ( ∂P ∂x + ∂Q ∂y ) dA = ∫ C ( P∂x ∂n+ Q ∂y ∂n ) ds を用いて，以下の通り変形できる． ∫ A [ ∂(δu· σx) ∂x + ∂(δu· τxy) ∂y + ∂(δv· σy) ∂y + ∂(δv· τxy) ∂x ] dA = ∫ s ( δu· σx· ∂x ∂n+ δu· τxy· ∂y ∂n+ δv· σy· ∂y ∂n+ δv· τxy· ∂x ∂n ) ds = ∫ s [δu(σx· ` + τxy· m) + δv(σy· m + τxy· `)]ds = ∫ s (δu· Sx+ δv· Sy)ds (3) また，ひずみ-変位関係 ∂u ∂x = x ∂v ∂y = y ∂v ∂x + ∂u ∂y = τxy より，式(2)右辺第2項は， ∫ A [ ∂(δu) ∂x σx+ ∂(δu) ∂y τxy+ ∂(δv) ∂y σy+ ∂(δv) ∂x τxy ] dA = ∫ A (δxσx+ δyσy+ δγxyτxy)dA (4) よって，仮想仕事の式は最終的に以下の形となる． ∫ A (δxσx+ δyσy+ δγxyτxy)dA = ∫ s (δuSx+ δvSy)ds + ∫ A式(5)をベクトル表示すると以下の通り． ∫ A {δ}T{σ}dA = ∫ s {δu}T{S}ds + ∫ A {δu}T{F }dA (6) 一般に３次元問題であっても，釣り合い状態にあるなら式(6)の面積積分を体積積分（A→ V）に，線積分 を面積積分（s→ A）とすることにより以下が成り立つ． ∫ V {δ}T{σ}dV = ∫ A {δu}T{S}dA + ∫ V {δu}T{F }dV (7) 板厚を考える場合，考える領域の板厚が一定値（= t）であれば，その方向をz方向として， ∫ V {δ}T{σ}dV = ∫ t 0 (∫ A {δ}T{σ}dA ) dz = t· ∫ A {δ}T{σ}dA となり，式(6)の両辺に板厚tを乗じればよいことになる． 1.2 要素剛性方程式 ■板厚について 実際の解析では平面応力状態で板厚の異なる要素を組み合わせた場合も考えられるため，以降各要素は一定の 板厚tを持つものとして式展開を進める． (1) 一般式 ここで，内挿関数[N ]を適切に選定し，要素内任意点の変位{u}および要素任意点のひずみ{}が，節点 変位{und}を用いて以下のように表現できるものとする． {u} =[N]{und} (8) {} =[B]{und} (9) 上記関係を用い，仮想仕事の式(6)を節点変位{und}を用いて表記すると {δund} ∫ A [B]T{σ}dA = {δund} ∫ s [N ]T{S}ds + {δund} ∫ A [N ]T{F }dA (10) よって釣合式は以下の通りとなる． ∫ A [B]T{σ}dA = ∫ s [N ]T{S}ds + ∫ A [N ]T{F }dA (11) ここで，応力-ひずみ関係行列[De]を導入し，要素内板厚を一定値tとすることにより，一般に知られてい る以下の要素剛性方程式が得られる． {fnd} =[k]{und} (12) [k] = t· ∫ A [B]T[De][B]dA {fnd} = t · ∫ s [N ]T{S}ds + t · ∫ A [N ]T{F }dA {} = [B]{und} {σ} = [De]{} = [De][B]{und} ここに， {fnd} ：節点外力ベクトル [k] ：要素剛性行列 {und} ：節点変位ベクトル [B] ：節点変位-要素任意点ひずみ関係行列 A ：要素面積 t ：要素板厚 {} ：要素任意点ひずみ {σ} ：要素任意点応力 [De] ：応力-ひずみ関係行列 [N ] ：内挿関数 {S} {F }
■要素表面力による節点外力ベクトル 要素任意表面の表面力{S}による節点外力ベクトル(t·∫C[N ]T{S}ds)の成分は，等価節点外力として直接 入力できる． ■物体力ベクトル 要素内任意点に作用する物体力（ここでは慣性力とする）{F } = {FBx, FBy}T は物体力に関わる節点力ベク トル{w}と内挿関数[N ]を用い以下の通り表せる． {F } = [N]{w} (13) よって物体力（慣性力）による等価節点外力ベクトル{FB}は， {FB} = t · ∫ A [N ]T{F }dA = t · ∫ A [N ]T[N ]dA· {w} (14) として算定される． (2) 初期ひずみを考慮する場合 温度などによる初期ひずみ{0}が存在する場合，要素応力は以下の通りとなる． {σ} = [De]{ − 0} = [De][B]{und} − [De]{0} (15) よって，初期ひずみによる節点外力ベクトルに相当する節点力{Tnd}は， {Tnd} = t · ∫ A [B]T[De]{0}dA (16) となり，初期ひずみを考慮した剛性方程式は，以下の通りとなる． {fnd} = [k]{und} − {Tnd} (17) 上式より節点変位は， {und} = [k]−1{fnd+ Tnd} (18) として求まるが，この節点変位から算定した応力は，初期ひずみを外力として加えた影響を含んでいる．この ため，構造系の応力は，上記節点変位より算定した応力から初期ひずみ相当応力を差し引いたもの，すなわち 式(15)より算定したものとする必要がある． (3) 応力-ひずみ関係 一般に３次元問題における等方均質弾性体の応力-ひずみ関係は，弾性係数をE，ポアソン比をν，材料の 線膨張係数をα，温度変化量をTとして以下の通りである（Hookeの法則）． x− αT = 1 E[σx− ν(σy+ σz)] y− αT = 1 E[σy− ν(σz+ σx)] z− αT = 1 E[σz− ν(σx+ σy)] γxy= 2(1 + ν) E τxy γyz= 2(1 + ν) E τyz γzx= 2(1 + ν) E τzx ここでx− y平面上を考える．平面応力状態では， σz= 0 τyz = 0 τzx= 0 であるため， x− αT = 1 E(σx− νσy) y− αT = 1 E(σy− νσx) γxy= 2(1 + ν) E τxy =⇒ σx σy τxy = E 1− ν2 1 ν 0 ν 1 0 0 0 1− ν 2 x− αT y− αT γxy
平面ひずみ状態では， z= 0→ σz= ν(σx+ σy)− EαT τyz= 0 τzx= 0 であるため， x− αT = 1 E[(1− ν 2)σ x− ν(1 + ν)σy] y− αT = 1 E[(1− ν 2)σ y− ν(1 + ν)σx] γxy= 2(1 + ν) E τxy =⇒ σx σy τxy = E (1 + ν)(1− 2ν) 1− ν ν 0 ν 1− ν 0 0 0 1− 2ν 2 x− (1 + ν)αT y− (1 + ν)αT γxy となる．これらの関係を再整理すると以下の通りとなる． ■２次元問題での応力とひずみの成分 応力成分：{σ} = σx σy τxy ひずみ成分：{} = x y γxy 平面応力温度ひずみ成分：{0} = αT αT 0 平面ひずみ温度ひずみ成分：{0} = (1 + ν)αT (1 + ν)αT 0 応力・ひずみの符号は，引張を正とし，温度の符号は温度上昇を正とする． ■２次元問題での応力-ひずみ関係式 平面応力：[De] = E 1− ν2 1 ν 0 ν 1 0 0 0 1− ν 2 Eν ：ポアソン比：弾性係数 平面ひずみ：[De] = E (1 + ν)(1− 2ν) 1− ν ν 0 ν 1− ν 0 0 0 1− 2ν 2
2.
４節点パラメトリック要素
2.1 ひずみ-変位関係行列[B]の導入四角形要素内任意点の変位u,vを以下のように仮定する．ここに座標(a,b)は[−1,1]の範囲を有する正規 化パラメータ座標，ui,j,k,lおよびvi,j,k,lは要素を構成する節点の変位を示す．
u =N1(a, b)· ui+ N2(a, b)· uj+ N3(a, b)· uk+ N4(a, b)· ul v =N1(a, b)· vi+ N2(a, b)· vj+ N3(a, b)· vk+ N4(a, b)· vl
要素内任意点の変位を{u}，節点変位を{und}とすれば， {u} = { u v } = [ N1 0 N2 0 N3 0 N4 0 0 N1 0 N2 0 N3 0 N4 ] ui vi uj vj uk vk ul vl = [N ]{und} (19) N1= 1 4(1− a)(1 − b) N2= 1 4(1 + a)(1− b) N3= 1 4(1 + a)(1 + b) N4= 1 4(1− a)(1 + b) (20) ∂N1 ∂a =− 1 4(1− b) ∂N2 ∂a = + 1 4(1− b) ∂N3 ∂a = + 1 4(1 + b) ∂N4 ∂a =− 1 4(1 + b) ∂N1 ∂b =− 1 4(1− a) ∂N2 ∂b =− 1 4(1 + a) ∂N3 ∂b = + 1 4(1 + a) ∂N4 ∂b = + 1 4(1− a) (21) 要素内任意点のひずみ{}は，節点変位{und}を用いて， {} = ∂u ∂x ∂v ∂y ∂u ∂y+ ∂v ∂x = ∂N1 ∂x 0 ∂N2 ∂x 0 ∂N3 ∂x 0 ∂N4 ∂x 0 0 ∂N1 ∂y 0 ∂N2 ∂y 0 ∂N3 ∂y 0 ∂N4 ∂y ∂N1 ∂y ∂N1 ∂x ∂N2 ∂y ∂N2 ∂x ∂N3 ∂y ∂N3 ∂x ∂N4 ∂y ∂N4 ∂x ui vi uj vj uk vk ul vl = [B]{und} (22) ここで，内挿関数Niに関する微分は，Jacobi行列[J ]を用いて以下の通り表せる． ∂Ni ∂a ∂Ni ∂b = [J ] ∂Ni ∂x ∂Ni ∂y ∂Ni ∂x ∂Ni ∂y = [J ]−1 ∂Ni ∂a ∂Ni ∂b (23) [J ] = ∂x ∂a ∂y ∂a ∂x ∂b ∂y ∂b = 4 ∑ i=1 ( ∂Ni ∂axi ) 4 ∑ i=1 ( ∂Ni ∂ayi ) 4 ∑ i=1 ( ∂Ni ∂bxi ) 4 ∑ i=1 ( ∂Ni ∂byi ) = [ J11 J12 J21 J22 ] (24) [J ]−1= 1 det(J ) [ J22 −J12 −J21 J11 ] det(J ) = J11· J22− J12· J21 (25)
式(22)(24)より， ∂Ni ∂x = 1 det(J ) { J22 ∂Ni ∂a − J12 ∂Ni ∂b } ∂Ni ∂y = 1 det(J ) { −J21 ∂Ni ∂a + J11 ∂Ni ∂b } (26) 式(24)の(xi,yi)を要素節点の座標とすれば，式(21)(24)より[J ]の要素が定まる．また，[J ]が定まれば， 式(21)(25)(26)より，式(22)で示されるひずみ-変位関係行列[B]が定まる． なお，この段階では[B]は変数aおよびbの関数として定義されていることに注意する． 以下に[J ]の要素を示しておく．ここにxi,j,k,lおよびyi,j,k,lは要素を構成する節点座標である． J11= ∂N1 ∂a xi+ ∂N2 ∂a xj+ ∂N3 ∂a xk+ ∂N4 ∂a xl J12= ∂N1 ∂a yi+ ∂N2 ∂a yj+ ∂N3 ∂a yk+ ∂N4 ∂a yl J21= ∂N1 ∂b xi+ ∂N2 ∂b xj+ ∂N3 ∂b xk+ ∂N4 ∂b xl J22= ∂N1 ∂b yi+ ∂N2 ∂b yj+ ∂N3 ∂b yk+ ∂N4 ∂b yl 2.2 要素剛性行列および節点力ベクトル 4節点パラメトリック要素の剛性行列[k]は，要素の板厚を一定値tとし，応力-ひずみ関係行列を[De]と すれば以下の通りとなる． [k] = t ∫ 1 −1 ∫ 1 −1 [B]T[De][B]· det(J) · da · db (27) また，初期ひずみ{0}による節点力ベクトル{Tnd}および要素応力{σ}による内力ベクトル{Rnd}は 以下の通りとなる． {Tnd} = t ∫ 1 −1 ∫ 1 −1 [B]T[De]{0} · det(J) · da · db (28) {Rnd} = t ∫ 1 −1 ∫ 1 −1 [B]T{σ} · det(J) · da · db (29) ■温度変化による初期ひずみ 平面応力：{0} = αTp αTp 0 = α· (N1· φi+ N2· φj+ N3· φk+ N4· φl) α· (N1· φi+ N2· φj+ N3· φk+ N4· φl) 0 平面ひずみ：{0} = (1 + ν)αTp (1 + ν)αTp 0 = (1 + ν)· α · (N1· φi+ N2· φj+ N3· φk+ N4· φl) (1 + ν)· α · (N1· φi+ N2· φj+ N3· φk+ N4· φl) 0 ここに，αは線膨張係数，νはポアソン比，Tpは要素任意点の温度変化量，{φi φj φk φl}Tは節点温度変化量， [N1 N2 N3 N4]は[N ]の要素である．また，温度変化量の符号は，温度上昇を正とする． 積分は以下に示すGauss積分により行う． ∫ 1 −1 ∫ 1 −1 f (a, b)· da · db + n ∑ i=1 n ∑ j=1 Hi· Hj· f(ai, bj) (30)
上式より [k]+ t · n ∑ i=1 n ∑ j=1
Hi· Hj· [B(ai, bj)]T[De][B(ai, bj)]· det(J(ai, bj)) (31)
{Tnd} + t · n ∑ i=1 n ∑ j=1 Hi· Hj· [B(ai, bj)]T[De]{0} · det(J(ai, bj)) (32) {Rnd} + t · n ∑ i=1 n ∑ j=1 Hi· Hj· [B(ai, bj)]T{σ} · det(J(ai, bj)) (33) ここで，積分点を4点(n = 2)とすれば，a，bおよび重みH の値は下表の通りであり，a，bを変化させた 4回の値の総和が積分の近似値となる． i j a b H 1 1 −0.5773502692 −0.5773502692 1.0000000000 1 2 +0.5773502692 −0.5773502692 1.0000000000 2 1 +0.5773502692 +0.5773502692 1.0000000000 2 2 −0.5773502692 +0.5773502692 1.0000000000 2.3 要素応力の算定 要素の応力は{σ}は，節点変位{u}を用いて以下の通り表される． {σ(ai, bj)} = [De][B(ai, bj)]{und} (34) ここで，ひずみ-変位関係行列[B]はaおよびbの関数であるため，{σ}もaおよびbの関数として表現さ れている．このため，4積分点要素の場合，4つの積分点での応力が定まる． 2.4 物体力による節点力 要素内任意点に作用する物体力（ここでは慣性力とする）{F } = {FBx, FBy}T は物体力に関わる節点力ベ クトル{w}と内挿関数[N ]を用い以下の通り表せる． {F } = [N]{w} (35) [N ] = [ N1 0 N2 0 N3 0 N4 0 0 N1 0 N2 0 N3 0 N4 ] (36) {w} = γ ·{kh kv kh kv kh kv kh kv }T (37) ここにγは要素の単位体積重量（質量×重力加速度），khおよびkvは水平方向・鉛直方向の慣性力の係数 （重力加速度に対する比率）である． 以上より，物体力（慣性力）による等価節点外力ベクトル{FB}は， {FB} = t · ∫ A [N ]T{F }dA = t · γ · ∫ A [N ]T[N ]dA· {w0} (38) として算定する．ここでベクトル{w}についてはγをくくり出すととにより， {w0} ={k h kv kh kv kh kv kh kv }T (39) とおいた．
3.
2
次元弾性応力解析の
No-tension
解析への拡張
3.1 解析の基本手順 引張に抵抗しない材料を含む構造の解析手法として，Zienkiewiczらが提案した応力遷移法を用いる．この 手法は，当初作成した線形弾性剛性行列および応力の座標変換のみでNo-tensionの釣り合い状態を得るもの である．また，No-tension状態となった要素では，ポアソン比を0としている． No-tension材料のイメージは，発生応力が材料の引張強度を超過した場合，引張応力は分担しないとする ものであり，例えばコンクリートの発生応力が引張強度以下の場合は引張応力を受け持つが，発生応力が引張 強度を超過しひび割れが発生した後は引張応力を受け持たなくなる挙動を想定している． 基本的な解析手順は，以下の通り．この手順では，No-tension状態が発生しない場合でも，最低2回の計 算が行われる．これは，温度変化および変位指定を与えた場合の処理を，荷重として考慮するのではなく，内 力計算過程に組み込んだためである．すなわち，温度変化および指定変位を用いて計算した内力は，これらを 考慮しない外力とは釣り合わず，1回目の繰り返しでは必ず不平衡力が生じ，2回目以降の繰り返し計算で釣 り合い状態が得られることになるためである． No-tension解析の基本手順 1 剛性行列[k]の組立． 2 初期の不平衡力{∆f}に，外力{f}および物体力ベクトル{fm}をセットする．また，全変位 {u}を零セットする． {∆f} = {f} + {fm} {u} = {0} 3 与えられた不平衡力ベクトル{∆f}に対し弾性剛性行列[k]を用い変位増分{∆u}を算定する． {∆u} = [k]−1{∆f} 4 算定された変位増分よりそのステップでの全変位，およびこれを用いた全ひずみを算定する．こ の際，既知節点変位（強制変位）があれば全変位はこれを考慮したものに修正する． {u} = {u} + {∆u} {} = [B]{u}5 算定された全ひずみを用いて応力を算定する．この時温度変化によるひずみを除去する． {σ} = [De]{ − 0} 6 算定された応力（主応力：σ1，σ2）が材料の引張強度σtsを超過した場合，算定された応力{σ} から引張成分{σt}を除去する． {σ} = {σ} − {σt} 7 引張成分を除去した応力から内力ベクトル{R}を算定する．内力ベクトルは，既知節点変位（強 制変位）および温度ひずみを考慮して算定しているため，釣り合い条件を満足しなければ，自動的 に不平衡力が発生する． {R} = t · ∫ A [B]T{σ}dA 8 不平衡力{∆f}に，初期の荷重（外力および物体力）と内力との差分をセットする． {∆f} = {f} + {fm} − {R} 9 不平衡力{∆f}あるいは変位増分{∆u}が十分小さければ計算終了とし，そうえなければ 3 に 戻り収束計算を行う．本プログラムでは，実際には不平衡力ではなく，変位増分が十分小さくなっ た時点で収束と判定している．
Setting of [k] ? {∆f} = {f} + {fm} ? {u} = {0} ? i=0 ? i=i+1 ? {∆u} = [k]−1{∆f} ? {u} = {u} + {∆u}
? {} = [B]{u} ? {σ} = [D]{ − 0} ? σ1, σ2< σts No Yes ? {R} = t ∫ A [B]T{σ}dA ? {∆f} = {f} + {fm} − {R} ? |∆f| < δ No Yes ? End -? {σ} = {σ} − {σt}
■収束計算手順 i − 1番目の不平衡点における応力をσi−1，不平衡力を∆fi−1= f− Ri−1とする． 不平衡力∆fi−1より変位増分∆uiを算定する． ∆uiより全変位ui，全ひずみiを求め，応力σiを算定する． 温度ひずみ0はこの応力算定時に考慮する． i番目の不平衡点における応力をσi= σi− σtとして算定する． ここにσtはσiの引張成分である． σiより内力ベクトルRiを算定する． ∆fi= f− Riとしてi番目の不平衡力を算定し，繰り返し計算を行う． 繰り返し過程における剛性行列[K]は常に一定である． -6 変位 外力・内力 [K ] [K ] [K ] 6 ? { R i− 1 } 6 ? { ∆ fi− 1 } 6 ? {ui} { R i } 6 ? { ∆ fi } - {∆ui} {f} 収束点 {σi} = [D]{i− 0} {σi} = {σi} − {σt} 図1 収束計算概念図
3.2 引張応力成分除去の方法 No-tension解析では内力ベクトルを求める際，引張応力成分を除去する必要がある．この方法としては， 算定した全応力{σx, σy, τxy}T に対し主応力およびその方向を計算し，主応力が引張状態であれば，これを {σx, σy, τxy}Tの成分に座標変換し，差し引くことにより引張応力成分を除去する．具体的手順を以下に示す． (1) 主応力の計算 水平方向をx軸，鉛直方向をy軸とする座標において，引張を正とするある応力状態が存在するとする．こ の状態に対し，反時計回りにθだけ回転させたところに主軸(I軸)があるとすると，主応力は以下の通り算定 できる． σ1= σx+ σy 2 + 1 2 √ (σx− σy)2+ 4τxy2 (40) σ2= σx+ σy 2 − 1 2 √ (σx− σy)2+ 4τxy2 (41) θ =1 2tan −1 ( 2τxy σx− σy ) (42) ここに， σx> σyかつτxy= 0なら θ = θ σx> σyかつτxy< 0なら θ = θ + 180◦ σx< σyなら θ = θ + 90◦ σx= σyかつτxy> 0なら θ = 45◦ σx= σyかつτxy< 0なら θ = 135◦ σx= σyかつτxy= 0なら 不定(θ = 0) (2) 引張応力成分の除去 水平方向をx軸，鉛直方向をy軸とする座標において，引張を正とするある応力状態が存在するとする．こ の座標を反時計回りにφだけ回転させた場合の応力を「0」をつけて表現すると以下の通りとなる． σ0x= σx+ σy 2 + σx− σy 2 cos 2φ + τxysin 2φ (43) σy0 =σx+ σy 2 − σx− σy 2 cos 2φ− τxysin 2φ (44) τxy0 =−σx− σy 2 sin 2φ + τxycos 2φ (45) ここで，回転前のx-y座標系が主軸とすれば， σx0 =σ1+ σ2 2 + σ1− σ2 2 cos 2φ (46) σ0y=σ1+ σ2 2 − σ1− σ2 2 cos 2φ (47) τxy0 =−σ1− σ2 2 sin 2φ (48) となり，座標回転後におけるσ1およびσ2の成分は以下の通りとなる． σ1の成分 σx0 = σ1 2(1 + cos 2φ) σy0 =σ1 2(1− cos 2φ) τxy0 =−σ1 2 sin 2φ (49) σ2の成分 σx0 = σ2 2(1− cos 2φ) σy0 =σ2 2(1 + cos 2φ) τxy0 = σ2 2 sin 2φ (50)
No-tension解析を行う場合，計算された応力{σx,σy,τxy}から引張成分を除去する必要がある．このため， σ1> 0(引張)の場合，σ1による引張成分の除去は以下の通り行う． σ1による引張成分の除去 σxc1= σx− σ1 2(1 + cos 2φ) σyc1= σy− σ1 2(1− cos 2φ) τxyc1= τxy+ σ1 2 sin 2φ (51) さらに，σ2が引張の場合， σ2による引張成分の除去 σc2 x = σc1x − σ2 2(1− cos 2φ) σc2 y = σc1x − σ2 2(1 + cos 2φ) τc2 xy= τxyc1− σ2 2 sin 2φ (52) ここでφ =−θである．また，σ1を最大主応力としているため，計算された応力から差し引く順序はσ1の 成分からとしている． 不平衡力を求めるための内力ベクトル{R}は，ここで計算した{σc1 x , σc1y , τxyc1}T あるいは{σxc2, σc2y , τxyc2}T を用いて計算する． 3.3 収束の考え方 収束判定は，変位と応力の双方で行うことが好ましいとの考えもあるが，本プログラムでは単純化を図るた め，変位増分が十分小さくなった場合に収束と判定することとした．具体的には，以下に示すように，変位指 定していない自由度箇所の収束過程における変位増分と，それまでの累計変位の比を検査している． 収束判定条件 ∑(∆ui)n {(∆ui)n} < 1 × 10−6 ここに，∆uiは変位増分， ∑ {∆ui}は収束計算過程での全累計変位（その時点での変位の解）であり，iは自 由度，nは繰り返し回数のインクリメントである． 実際の構造解析では，構造物の寸法は，m単位あるいはmm単位で入力することを考えると，10−3mm程 度の変位の精度があれば差し支えないという考えに基づいてこの値を定めた．
4.
プログラミング
4.1 ２次元応力解析用subroutine (1) subroutine XGYG3M(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod) 三角形要素における要素重心計算．要素応力出力時に要素重心座標を併記している． (2) subroutine XGYG4M(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod) 四角形要素における要素重心計算．要素応力出力時に要素重心座標を併記している． （ただしIPR=1：四角形要素の計算で応力の出力を要素平均に指定した場合） (3) subroutine XGYG4I(ne,kk,kakom,x,y,xg,yg,NELT,NODT,nod) Gauss積分点指定．要素応力出力時にGauss積分点座標を併記している．(4) subroutine DMAT PL(NSTRES,Eme,poe,d)
三角形・四角形要素における応力ーひずみ関係行列
(5) real(8) function CALA PL3(ne,kakom,x,y,NELT,NODT)
三角形要素における要素面積計算
(6) subroutine BMAT PL3(ne,kakom,x,y,b,NELT,NODT)
三角形要素におけるひずみ-変位関係行列作成
(7) subroutine SM PL3(ne,kakom,x,y,NSTRES,Em,po,t,sm,NELT,NODT)
三角形要素における要素剛性行列計算
(8) subroutine MV PL3(ne,kakom,x,y,t,gamma,gkh,gkv,fevec,NELT,NODT)
三角形要素における要素慣性力ベクトル計算
(9) subroutine CALST PL3TS(ne,kakom,x,y,deltaT,dist,NSTRES,Em,po,t,alpha,ts, strsave,fevec,noten,NELT,NODT,nt)
三角形要素における要素内力計算
(10) subroutine BMAT PL4(ne,kakom,x,y,bm,detJ,a,b,NELT,NODT)
四角形要素におけるひずみ-変位関係行列作成 (11) subroutine IPOINT4(kk,a,b) 四角形要素の平面積分におけるGauss積分点指定 (12) subroutine SM PL4(ne,kakom,x,y,NSTRES,Em,po,t,sm,NELT,NODT) 四角形要素における要素剛性行列計算 (13) subroutine MV PL4(ne,kakom,x,y,t,gamma,gkh,gkv,fevec,NELT,NODT) 四角形要素における要素慣性力ベクトル計算
(14) subroutine CALST PL4TS(ne,kakom,x,y,deltaT,dist,NSTRES,Em,po,t,alpha,ts, strsave,fevec,noten,NELT,NODT,nt)
(15) subroutine PST PL(sigx,sigy,tauxy,ps1,ps2,ang)
主応力計算
4.2 FEMプログラム共通のsubroutineおよびfunction
(1) subroutine ASVEC(ne,nod,nhen,kakom,fevec,ftvec,NELT,NODT) 要素節点ベクトル（vector）の全体節点ベクトルへの組み込み. 浸透流解析プログラムでは使用しない． (2) subroutine ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT) 要素行列（matrix）の全体行列への組み込み (3) subroutine CHBAND10(n1,m1,array,nt) Cholesky法三角行列作成．バンドマトリックスを用いたCholesky法によるルーチンを２分割し，時刻や 収束計算過程により変化しない三角行列作成部分は１回のみの呼び出しとし，以降保存利用する． (4) subroutine CHBAND11(n1,m1,array,vector,nt) Cholesky法前進消去・後退代入．時刻や収束計算過程で条件変化する部分を各過程毎に更新しながら解を 求める．
(5) integer function IBAND(mm,array,fact0,nt)
Cholesky法適用のため，バンド幅（格納変数：ib）の計算を行う整数型function．入力は，拘束節点処理 後の縮小したmm×mmの正方行列であり，行毎に最終列から連続して0となっている列数を数え，バンド幅 を算定している．数値が0であるか否かの判定は，各行列要素数値の絶対値が1.0×10−30であるか否かを，0 判定数値格納変数fact0との比較で判定している． (6) subroutine BAND(mm,ib,array,nt) Cholesky法適用のための，フルマトリックスのバンド化記憶を行う．入力は，拘束節点処理後の縮小した mm×mmの正方行列であり，出力はCholesky法を適用するmm×ibのサイズの行列となる．なお計算に用 いない要素には0を格納している．
(7) subroutine del spaces(s)
付録
A
三角形定ひずみ要素
A.1 変位関数 三角形定ひずみ要素では，１節点２自由度３節点であるため，１要素あたり６自由度を有する．このため， ６個の未定係数α1∼6を用い，以下の変位を仮定する． u =α1+ α2· x + α3· y (53) v =α4+ α5· x + α6· y (54) ここで，三角形の３つの頂点をi，j，kと反時計回りに設定し，それぞれの節点変位(u, v)を節点座標(x, y) とαで表せば以下の通りとなる． ui vi uj vj uk vk = 1 xi yi 0 0 0 0 0 0 1 xi yi 1 xj yj 0 0 0 0 0 0 1 xj yj 1 xk yk 0 0 0 0 0 0 1 xk yk α1 α2 α3 α4 α5 α6 (55) これを{α}について解けば {α} = [C]{und} (56) ここに， [C] = 1 2∆ ai 0 aj 0 ak 0 bi 0 bj 0 bk 0 ci 0 cj 0 ck 0 0 ai 0 aj 0 ak 0 bi 0 bj 0 bk 0 ci 0 cj 0 ck 2∆ = 1 xi yi 1 xj yj 1 xk yk (57) ai= xjyk− xkyj bi= yj− yk ci= xk− xj aj = xkyi− xiyk bj = yk− yi cj= xi− xk ak= xiyj− xjyi bk = yi− yj ck= xj− xi (58) ここに，∆は三角形要素の面積，{und}は節点変位ベクトル，{α}は一般化変位ベクトルである．また， 2∆を直接の計算式で示すと以下の通りである． 2∆ = (xk− xj)yi+ (xi− xk)yj+ (xj− xi)yk (59) ■参考：Sarrusの法則（3×3の行列式の値） a11 a12 a13 a21 a22 a23 a31 a32 a33 = a11a22a33+ a21a32a13+ a31a12a23− a13a22a31− a23a32a11− a33a12a21 以上により，要素内任意点の変位は以下の通り表される． { u v } = [ 1 x y 0 0 0 0 0 0 1 x y ] {α} = [ 1 x y 0 0 0 0 0 0 1 x y ] [C]{und} = [N]{und} (60) ここに[N ]は内挿関数であり，要素は以下の通り． [N ] = [ 1 x y 0 0 0 0 0 0 1 x y ] [C] = [ Ni 0 Nj 0 Nk 0 0 Ni 0 Nj 0 Nk ] (61) Ni= ai+ bix + ciy 2∆ Ni= aj+ bjx + cjy 2∆ Ni= ak+ bkx + cky 2∆ (62)A.2 ひずみ-変位関係行列 要素ひずみ{}は変位(u, v)を微分することにより以下の通り表される． {} = ∂u ∂x ∂v ∂y ∂u ∂y+ ∂v ∂x = ∂Ni ∂x 0 ∂Nj ∂x 0 ∂Nk ∂x 0 0 ∂Ni ∂y 0 ∂Nj ∂y 0 ∂Nk ∂y ∂Ni ∂y ∂Ni ∂x ∂Nj ∂y ∂Nj ∂x ∂Nk ∂y ∂Nk ∂x {und} = [B]{und} (63) よってひずみ-変位関係は以下の通りとなる． [B] = 1 2∆ b0i c0i b0j c0j b0k c0k ci bi cj bj ck bk = 1 2∆ yj− y0 k xk− x0 j yk− y0 i xi− x0 k yi− y0 j xj− x0 i xk− xj yj− yk xi− xk yk− yi xj− xi yi− yj (64) ∆ = 1/2· [(xk− xj)yi+ (xi− xk)yj+ (xj− xi)yk] (65) A.3 要素剛性行列および節点力ベクトル 要素剛性行列[k]は，ひずみ-変位関係行列[B]および応力-ひずみ関係行列[De]を用い，以下の通りとなる [k] = t· ∫ A [B]T[De][B]dA = t· ∆ · [B]T[De][B] (66) また，初期ひずみ{0}による節点力ベクトルおよび要素応力{σ}による内力ベクトル{Rnd}は以下の通 りとなる． {Tnd} = t · ∫ A [B]T[De]{0}dA = t · ∆ · [B]T[De]{0} (67) {Rnd} = t · ∫ A [B]T{σ}dA = t · ∆ · [B]T{σ} (68) ■温度変化による初期ひずみ 平面応力：{0} = αTm αTm 0 = α· (φi+ φj+ φk)/3 α· (φi+ φj+ φk)/3 0 平面ひずみ：{0} = (1 + ν)αTm (1 + ν)αTm 0 = (1 + ν)· α · (φi+ φj+ φk)/3 (1 + ν)· α · (φi+ φj+ φk)/3 0 ここに，αは線膨張係数，νはポアソン比，Tmは要素を取り囲む節点の平均温度変化量，{φi, φj, φk}Tは節点温 度変化量である．また，温度変化量の符号は，温度上昇を正とする． A.4 物体力による節点力 (1) 面積座標 慣性力など物体力による節点力ベクトルを求める場合，内挿関数の積分が必要となる．このため，内挿関数 の書き換えを行う． [ Ni 0 Nj 0 Nk 0 ] [ ζi 0 ζj 0 ζk 0 ]
ζi= ∆i ∆ ζj = ∆j ∆ ζk = ∆k ∆ (70) 2∆i= ai+ bix + ciy = 1 x y 1 xj yj 1 xk yk (71) 2∆j = aj+ bjx + cjy = 1 x y 1 xk yk 1 xi yi (72) 2∆k = ak+ bkx + cky = 1 x y 1 xi yi 1 xj yj (73) この座標(ζi, ζj, ζk)は面積座標を表し，明らかに ζi+ ζj+ ζk= 1 05 ζr5 1 (r = i, j, k) (74) である． (2) 面積座標で表された内挿関数の積分 物体力を求める場合，∫A[N ]T[N ]dAの形の積分が必要となるため，まず[N ]T[N ]を求める． [N ]T[N ] = ζi 0 0 ζi ζj 0 0 ζj ζk 0 0 ζk [ ζi 0 ζj 0 ζk 0 0 ζi 0 ζj 0 ζk ] = ζi2 0 ζiζj 0 ζkζi 0 0 ζi2 0 ζiζj 0 ζkζi ζiζj 0 ζj2 0 ζjζk 0 0 ζiζj 0 ζj2 0 ζjζk ζkζi 0 ζjζk 0 ζk2 0 0 ζkζi 0 ζjζk 0 ζk2 (75) 積分変数の変換公式より ∫∫ A f (x, y)dxdy = ∫ 1 0 ∫ 1 0 g(ζi, ζj, ζk) ∂(x, y) ∂(ζi, ζj) dζidζj (76) 要素内任意点を内挿関数で表し，ζi+ ζj+ ζk= 1の関係を用いると x =ζixi+ ζjxj+ ζkxk → x = ζi(xi− xk) + ζj(xj− xk) + xk y =ζiyi+ ζjyj+ ζkyk → y = ζi(yi− yk) + ζj(yj− yk) + yk (77) ヤコビアンは ∂(x, y) ∂(ζi, ζj) = ∂x ∂ζi ∂x ∂ζj ∂y ∂ζi ∂y ∂ζj = xi− xk xj− xk yi− yk yj− yk =(xk− xj)yi+ (xi− xk)yj+ (xj− xi)yk = 2∆ (78) よって ∫∫ A f (x, y)dxdy = 2∆ ∫ 1 0 ∫ 1 0 g(ζi, ζj, ζk)dζidζj (79) また面積座標(ζi, ζj, ζk)に対し，以下の公式が成り立つ． ∫ 1 0 ∫ 1 0 ζi`· ζjm· ζkn· dζidζj= `!· m! · n! (` + m + n + 2)! （`, m, nは負でない整数） (80) たとえば ` = 2, m = n = 0 → ∫ 1 0 ∫ 1 0 ζi`· ζjm· ζkn· dζidζj= 1/12 ` = 1, m = 1, n = 0 → ∫ 1 0 ∫ 1 0 ζi`· ζjm· ζkn· dζidζj= 1/24
したがって（ζでの積分ではなくAでの積分であることに注意） ∫ A ζi2dA = ∫ A ζj2dA = ∫ A ζk2dA = 1/12· 2∆ = 1/6 · ∆ ∫ A ζiζjdA = ∫ A ζjζkdA = ∫ A ζkζidA = 1/24· 2∆ = 1/12 · ∆ 以上により ∫ A [N ]T[N ]dA = ∆ 1/6 0 1/12 0 1/12 0 0 1/6 0 1/12 0 1/12 1/12 0 1/6 0 1/12 0 0 1/12 0 1/6 0 1/12 1/12 0 1/12 0 1/6 0 0 1/12 0 1/12 0 1/6 (81) (3) 物体力による節点力ベクトル 物体力（慣性力）による節点力ベクトル{FB}は， {FB} = t · ∫ A [N ]T{F }dA = t · ∫ A [N ]T[N ]dA· {w} (82) と表せる． ここで，ベクトル{w}を以下の通りおく． {w} = γ ·{kh kv kh kv kh kv }T (83) ここにkhは水平方向加速度の重力加速度に対する比率（右向き：正），kvは鉛直方向加速度の重力加速度に 対する比率（上向き：正）であり，γは要素の単位体積重量（質量×重力加速度）である． これにより物体力（慣性力）による節点力ベクトル{FB}は以下の通り定義される． {FB} = γ· ∆ · t 3 1/2 0 1/4 0 1/4 0 0 1/2 0 1/4 0 1/4 1/4 0 1/2 0 1/4 0 0 1/4 0 1/2 0 1/4 1/4 0 1/4 0 1/2 0 0 1/4 0 1/4 0 1/2 kh kv kh kv kh kv (84) ここで，∆は要素面積，tは板厚であり，要素の単位体積重量γは要素面積・板厚と共にくくりだしている．
参考文献
[1] 三本木茂夫・吉村信敏：コンピュータによる構造工学講座I-1-B 日本鋼構造協会編 有限要素法による
構造解析プログラム 考え方と解説，培風館，昭和54年5月
[2] 山田嘉昭：有限要素法の基礎と応用シリーズ２ マトリックス法材料力学，培風館，昭和55年5月
[3] S.P. Timoshenko, J.N. Goodier : Theory of Elasticity Third edition, McGRAW-HILL BOOK CAM-PANY, 21st printing 1985 [4] 宮沢健二：パソコン構造力学–力学挙動の可視化–，啓学出版株式会社，1984年6月 [5] H.C. Martin / G.F. Carey共著，鷲津久一郎/山本善之共訳：有限要素法の基礎と応用，培風館，昭和 56年10月 [6] O.C.Zienkiewicz著，吉識雅夫/山田嘉昭監訳：基礎工学におけるマトリックス有限要素法，培風館，昭 和50年12月 [7] 竹内則雄：コンピュータによる極限解析法シリーズ４ 地盤力学における離散化極限解析，培風館， 1991.7 [8] 小堀為雄・吉田博,有限要素法による構造解析プログラム,丸善株式会社,昭和55年12月
