局所座標系のはなし

STEP3:弾性方程式の離散化(後編)

弾性方程式の離散化

前編で述べたようにある構造物に仮想変位が生じているとしてもその構造物がつりあいの状態にあるとみなすことができるので、 (1)δW+δU=0 が成り立ちます。このとき外力Tiによって仮想変位uを生じたとすると、この外力の仮想仕事量は、 (2)δW=(δu)TTdΓ ですから、両式を式(1)に代入すると、 Se(δu)TTdΓ+Ve{δεTσ(δu)TF}dΩ=0Se(Nδue)TTdΓ+Ve{(Bδue)Tσ(Nδue)TF}dΩ=0(δue)TSeNTTdΓ+(δue)TVe(BTσNTF)dΩ=0(3)(δue)T{SeNTTdΓ+Ve(BTσNTF)dΩ}=0 となります。これは任意の仮想変位についてなり立たなければならないので、 (4)SeNTTdΓ+Ve(BTσNTF)dΩ=0 が得られます。応力を節点変位で表した式 (5)ε=[B1B2B20]ue=Bue を代入すると式(4)は、 (6)SeNTTdΓ+Ve(BTDBueNTF)dΩ=0(7)VeBTDBdΩue=VeNTFdΩSeNTTdΓ となります。したがって (8)Ke=VeBTDBdΩ,Fe=VeNTFdΩ,Te=SeNTTdΓ と置くことにより、式(7)は、 (9)Keue=Fe+Te と表され、これが要素eについて弾性方程式を離散化した式となっています。

数値積分法

各要素の節点変位は式(9)で定義されましたが、この積分方程式をコンピュータで算出するために、Gaussの数値積分を活用します。式(9)は直行座標系の式であるため、まずは前編で定義したBマトリクスの成分を局所座標系で表現する工夫が必要となります。これに伴って、各要素の体積(あるいは面積)積分の被積分領域も局所座標系に変換します。

局所座標への変換

節点iの形状関数Niは、 (10)Ni=Ni(x,y,z),x=x(ξ,η,ζ),y=y(ξ,η,ζ),z=z(ξ,η,ζ) であるため、1次偏導関数は、 Niξ=Nixxξ+Niyyξ+Nizzξ(11)Niη=Nixxη+Niyyη+NizzηNiζ=Nixxζ+Niyyζ+Nizzζ と表わされます。したがってこれをマトリクス表示すると、 [NiξNiηNiζ]=[xξyξzξxηyηzηxζyζzζ][NixNiyNiz]=[J11J12J13J21J22J23J31J32J33][NixNiyNiz](12)=J[NixNiyNiz] と表すことができます。上式で(x,y,z)直交座標座標を形状関数を用いて(ξ,η,ζ)局所座標で表現すると、(x,y,z)は形状関数を用いて式(???)によって定義されるため、式(12)のJacobianマトリクスの成分は、 J11=i=120NiξxiJ12=i=120NiξyiJ13=i=120Niξzi(13)J21=i=120NiηxiJ22=i=120NiηyiJ23=i=120NiηziJ31=i=120NiζxiJ32=i=120NiζyiJ33=i=120Niζzi となります。したがってJacobianの逆行列が定まれば式(12)からすぐさま、Bマトリクスの成分が求まることになります。そこで逆行列と行列式を求めると、 J1=1|J|[J22J33J23J32(J21J33J23J31)J21J32J22J31(J12J33J13J32)J11J33J13J31(J11J32J12J31)J12J23J13J22(J11J23J13J21)J11J22J12J21]T(14)=1|J|[J22J33J23J32J13J32J12J33J12J23J13J22J23J31J21J33J11J33J13J31J13J21J11J23J21J32J22J31J12J31J11J32J11J22J12J21] (15)|J|=J11J22J33J11J23J32+J12J23J31J12J21J33+J13J21J32J13J22J31 となるので、これらを次式 (16)[NixNiyNiz]=J1[NiξNiηNiζ] に代入することでBマトリクスの成分は、 Nix=1|J|{(J22J33J23J32)Niξ+(J13J32J12J33)Niη+(J12J23J13J22)Niζ}(17)Niy=1|J|{(J23J31J21J33)Niξ+(J11J33J13J31)Niη+(J13J21J11J23)Niζ}Niz=1|J|{(J21J32J22J31)Niξ+(J12J31J11J32)Niη+(J11J22J12J21)Niζ} と定まり、これが局所座標系で表わされた式となっています。 次に、被積分領域を局所座標へ変換することについて論じます。まず、被積分領域が体積である場合、直行座標系で任意の座標を、 (18)r=xi+yj+zk,x=x(ξ,η,ζ),y=y(ξ,η,ζ),z=z(ξ,η,ζ) とすると局所座標の任意の座標(drξ,drη,drζ)は、 drξ=rξdξ=(xξi+yξj+zξk)dξ(19)drη=rηdη=(xηi+yηj+zηk)dηdrζ=rζdζ=(xζi+yζj+zζk)dζ ですので、直行座標系の体積dΩを局所座標系に変換したものは、(x,y,z)座標において、 dΩ=drξ×drηdrζ(20)=|ijkxξyξzξxηyηzη|[xζyζzζ]=|xξyξzξxηyηzηxζyζzζ|dξdηdζ=|J|dξdηdζ と表わされます。 次に被積分領域が面積である場合、直行座標系で r=xi+yj+zk,x=x(ξ,η,1),y=y(ξ,η,1),z=z(ξ,η,1) と表わされるベクトルは、(ξ,η,ζ)空間における局所座標平面ζ=1を表します。そこで表面荷重のかかる要素面をこれにとると、平面ζ=1上の点では、N9N20=0ですから、この面に応じた(x,y,z)座標は、 x=i=120Nixi=i=18Nixi,y=i=120Niyi=i=18Niyi,z=i=120Nizi=i=18Nizi と表すことができます。このとき、 dr1=rξdξ=(xξi+yξj+zξk)dξ(21)dr2=rηdη=(xηi+yηj+zηk)dη と定義したdr1dr2(ξ,η,ζ)空間における局所座標平面ζ=1上のξ方向のベクトルと、η方向のベクトルとなりますから、 dr1dr2で張られる面の法線方向の成分dΓp(x,y,z)空間において、 dΓp=dr1×dr2=|ijkxξyξzξxηyηzη|dξdη=|ijki=18Niξxii=18Niξyii=18Niξzii=18Niηxii=18Niηyii=18Niηzi|dξdηdζ(22)=|ijkJ11J12J13J21J22J23|dξdη で表されます(ただし方向はdr1からdr1に右ねじを回して進む向きです。そこでこの方向をp方向としています)。これをベクトル表示したものを、 (23)dΓp=[J12J23J13J22J13J21J11J23J11J22J12J21]dξdη=[VxpVypVzp]dξdη と定義しておきます。 このときp方向、ξ方向、η方向それぞれの単位ベクトルは(x,y,z)座標において (24)np=1(Vxp)2+(Vyp)2+(Vzp)2[VxpVypVzp]=[npxnpynpz](25)nξ=1J112+J122+J132[J11J12J13]=[nξxnξynξz](26)nη=1J212+J222+J232[J21J22J23]=[nηxnηynηz] です。また、 (27)dΓξ=nη×dΓp,dΓη=dΓp×nξ ですから、式(23)(26)を代入すると、 (28)dΓξ=[nηxnηynηz]×[VxpVypVzp]dξdη=[nηyVzpnηzVypnηzVxpnηxVzpnηxVypnηyVxp]dξdη=[VxξVyξVzξ]dξdη(29)dΓη=[VxpVypVzp]dξdη×[nξxnξynξz]=[nξzVypnξyVzpnξxVzpnξzVxpnξyVxpnξxVyp]dξdη=[VxηVyηVzη]dξdη となり、これをξη平面に係る表面荷重の演算に用います。

Gaussの数値積分

Gaussの数値積分は多項式f(x)に対して、 (30)Δ=11f(x)dxi=1nf(xi)wiΔを最小にするように幾何対称を利用して表2.1のようにGauss point xiとその重みwiを定めるため、多項式の積分はGaussの数値積分によって厳密な解を得ることができます。 3次元問題においてはn=3で十分であり、次の積分: (31)I=111111f(ξ,η,ζ)dξdηdζ を求めるにあたり、まず、ηζを固定したξのみの関数とみて重みwiを考慮して、次に示すように内側の積分を行います。 (32)11f(ξ,η,ζ)dξ=inf(ξi,η,ζ)wi=ψ(η,ζ) すると、Iは、 (33)I=1111ψ(η,ζ)dηdζ となります。同様に、ζを固定してηのみの関数とみると、 (34)11ψ(η,ζ)dη=jnψ(ηj,ζ)wj=ϕ(ζ) となります。この式に式(33)(34)を代入すると (35)I=knjninf(ξi,ηj,ζk)wiwjwk としてGaussの数値積分を表す式が得られます。表1にGauss Pointとその重みを示します。

表1 Gauss Pointとその重み
n ±xi wi
1 0.00000 00000 00000 2.00000 00000 00000
2 0.57735 02691 89626 1.00000 00000 00000
3 0.77459 66692 41483 0.55555 55555 55555
0.00000 00000 00000 0.88888 88888 88889
4 0.86113 63115 94053 0.34785 48451 37454
0.33998 10435 84856 0.65214 51548 62546
5 0.90617 98459 38664 0.23692 68850 56189
0.53846 93101 05683 0.47862 86704 99366
0.00000 00000 00000 0.56888 88888 88889

要素剛性マトリクスと荷重ベクトルのGauss積分表示

前節では離散化された弾性方程式を式(9)として示しました。この要素剛性マトリクスと荷重ベクトルを局座標変換してGauss積分にて表示し、コンピュータプログラムにセットできる状態まで、より具体的な展開を行います。

要素剛性マトリクス

要素剛性マトリクスについては被積分領域が体積ですから、式(20)の変換により、 Ke=VeBTDBdΩ=111111BTDB|J|dξdηdζ(36)=k=13j=13i=13{BTDB}ijk|J|wiwjwk\label{64} と表わされます。Bマトリクスが6×60Dマトリクスが6×6の行列ですから、上式(36)の{ }の中身は60×60の行列となっています。この行列を (37)Ke=[Kmne](m=1,60n=1,60) と表すと、その成分は、 (38)Kmne=k=13j=13i=13{p=16(BmpTq=16DpqBqn)}ijk|J|wiwjwk で与えられます。このとき材料マトリクスD=[Dij]とは、 (39){Dij=Eν(1+ν)(12ν)(i=13j=13ij)Dij=E(1ν)(1+ν)(12ν)(i=13i=j)Dij=E2(1+ν)(i=46i=j)Dij=0(others) と表されるものです。

荷重ベクトル

荷重ベクトルも被積分領域が体積ですから、式(20)の変換により、 Fe=VeNTFdΩ=111111NT|J|dξdηdζF(40)=[k=13j=13i=13NijkT|J|wiwjwk]F と表わされます。この行列を (41)Fe=[Fm1e](m=1,60) と表すと、その成分は、 (42)Fm1e=[n=13{k=13j=13i=13(NmnT)ijk|J|wiwjwk}]Fn1 で与えられます。なお、この荷重ベクトルは重力などの体積力によるものです。

表面荷重ベクトル

表面荷重ベクトルは、 Te=SeNTTdΓ(43)=1111NTNdΓξTξ1111NTNdΓηTη1111NTNdΓpTp のように局座標で表現できます。この被積分領域は面積ですから、式(23)(28)(29)の変換により、ξη平面に係る表面荷重は、 (44)Te=Se[N1IN2IN8I][VxξVyξVzξ][N1N2N8]dξdη[T1ξT2ξT8ξ]Se[N1IN2IN8I][VxηVyηVzη][N1N2N8]dξdη[T1ηT2ηT8η]Se[N1IN2IN8I][VxpVypVzp][N1N2N8]dξdη[T1pT2pT8p] で与えられます。