有限要素法は分かり難い理論です。でも、ちゃんと理解していないと、現実の問題にに適用したときに、公式だけを暗記している高校生のように頓珍漢な答えになってしまう危険性を内包しています。だれが使っても正しく動作する fool proof な有限要素法のソフトは作れないでしょう。
有限要素法を理解するためには、単純な問題を理論に従った自分で計算してみることが良い方法です。納得できるまで、ああだこうだと数値実験をしてみるべきです。でも電卓では有限要素法の計算するには手間がかかりすぎます。コンピュータのプログラムを作れば計算の手間が省けますが、デバッグの手間がそれ以上にかかってしまいます。そこで sf です。sf を使えば、電卓のような計算の手間もかかけずに、プログラム開発でのデバッグをすることなく、有限要素法の数値実験が可能です。
以下、sf による有限要素法の数値実験を truss に適用した例を示していきます。
最初に有限要素法の勉強を始めたとき「トラス要素の誘導」を読みました。解かりやすく書くことに力を注いでいます。有限要素法が初めての方は、私の説明では解かり難いときに、こちらの説明も読んでみてください。
「トラス要素の誘導」を参照するとき楽なように、「例題 1」についての問題については、形状寸法、材料定数を意識して同じにしておきました。記号もなるべく一致させておきました。
平面に置かれた棒の両端に力が、どのように変移するかを考えます。棒の両端の位置と、両端の変移と Young 率と両端に加わる力の関係を求めます。これを求めると、棒の組み合わせで作られた truss では、個別の棒要素の位置と変移と力の関係が線形に重ね合わせたものとして、truss として組み合わさった条件での関係を導けます
y 軸
↑
│
│ <x0,y0> <x1,y1> 位置
│ ← <Fx0,Fy0> → <Fx1,Fy1> 力
│ ━━━━━━━
│ <Ux0,Uy0> <Ux1,Uy1> 変移
│
─┼──────────→ x 軸
│
│
棒の断面積を A
長さを L
Young 率を E
とします。
<Fx0,Fy0, Fx1,Fy1> = E A/L Mat <Ux0,Uy0, Ux1,Uy1> --- (4 式)
Mat = < 1, 0,-1, 0>
< 0, 0, 0, 0>
<-1, 0, 1, 0>
< 0, 0, 0, 0>
水平に置かれていること、釣り合い状態にあることから
y0 == y1
Uy0 == Uy1 == 0
Fy0 == Fy2 == 0
が成り立ちます。
Trus として組み合わさった状態でも、変移が小さいときは 4 式が成り立っています。逆に変移が大きくなると 4 式が成り立ちません。バネのように何倍にも伸びる変形は 4 式がなりたたなくなります。有限要素法ではでは(4 式)が成り立つとき小さな変移を対象とします。
<x0,y0> <x1,y1> 位置
← <Fx0,Fy0> → <Fx1,Fy1> 力
━━━━━━━
y'軸 <Ux0,Uy0> <Ux1,Uy1> 変移
／
／
＼
＼ x'軸
力と変移の関係式下のようになります。
Rotate = < cos(φ), -sin(φ), 0, 0>
< sin(φ), cos(φ), 0, 0>
<0, 0, cos(φ), -sin(φ)>
<0, 0, sin(φ), cos(φ)>
# @Rotate(φ)^-1 を変移に作用させることで、棒の軸方向が x 軸になるような座標系での変移を表します
# @Rotate(φ)^-1 を力 に作用させることで、棒の軸方向が x 軸になるような座標系での力 を表します
Rotate^1 <Fx0,Fy0, Fx1,Fy1> = E A/L Mat Rotate^1 <Ux0,Uy0, Ux1,Uy1>
Mat = < 1, 0,-1, 0>
< 0, 0, 0, 0>
<-1, 0, 1, 0>
< 0, 0, 0, 0>
<Fx0,Fy0, Fx1,Fy1> = E A/L Rotate Mat Rotate^-1 <Ux0,Uy0, Ux1,Uy1> --- (7 式)
Mat = < 1, 0,-1, 0>
< 0, 0, 0, 0>
<-1, 0, 1, 0>
< 0, 0, 0, 0>
E A/L Rotate^-1 Mat Rotate のことを要素剛性行列と言います。Truss 全体での剛性行列は、この要素剛性行列を線形に足し合わせてやることで計算できます。要素ごとの変移と力の関係を線形に重ね合わせた状態が truss で発生しているはすだからです。
2
5000mm／|→ Fx3 == 2000kgf
③／ |②
／ | 4000mm
0-------+1
△ ① ◎
3000mm
断面積: A=100 mm^2, Young 率: E=21000 Kg/mm^2 とします
剛性行列 Ke は節点 0,1,2 の変移に対して、節点 0,1,2 に加わる力の関係式を定めます。
Ke は 6 x 6 正方行列
力 = Ke * 変移
<Fx0,Fy0, Fx1,Fy1, Fx2,Fy2> = Ke <Ux0,Uy0, Ux1,Uy1, Ux2,Uy2>
Ke を求めるため、① ② ③ 三つの棒の要素剛性行列 Ke1, Ke2, Ke3 を求めます。まず、水平な棒の変移の差分を 計算する Mat を作ります。同時に、棒の断面積、Young 率の sf ファイル変数: A.val, E.val も宣言しておきます。
//@@
Mat =[[4]]
Mat[0,*] = < 1, 0, -1, 0 >
Mat[1,*] = < 0, 0, 0, 0 >
Mat[2,*] = < -1, 0, 1, 0 >
Mat[3,*] = < 0, 0, 0, 0 >
A=100 # mm^2 断面積
E=21000 # Kg/mm^2 Young 率
//@@@
この Mat, E, A より棒 ① の要素剛性行列 Ke1 は下の様に求められます。
Ke1 の計算;; L@=3000, Ke1= E A/L Mat
< 700, 0, -700, 0 >
< 0, 0, 0, 0 >
< -700, 0, 700, 0 >
< 0, 0, 0, 0 >
6 x 6 の剛性行列 Ke の節点 0, 1 に対応する部分に Ke1 要素を代入します。
//@@
Ke = [[6]]
##要素 ① の変移(節店 0,1)と力の関係を Ke に設定する
#要素 ① による力 Fx0 と変移 Ux0,Uy0, Ux1,Uy1 の関係 Ke1 から Ke に移す
Ke[0,0] = Ke1[0,0]
Ke[0,1] = Ke1[0,1]
Ke[0,2] = Ke1[0,2]
Ke[0,3] = Ke1[0,3]
#要素 ① による力 Fy0 と変移 Ux0,Uy0, Ux1,Uy1 の関係 Ke1 から Ke に移す
Ke[1,0] = Ke1[1,0]
Ke[1,1] = Ke1[1,1]
Ke[1,2] = Ke1[1,2]
Ke[1,3] = Ke1[1,3]
#要素 ① による力 Fx1 と変移 Ux0,Uy0, Ux1,Uy1 の関係 Ke1 から Ke に移す
Ke[2,0] = Ke1[2,0]
Ke[2,1] = Ke1[2,1]
Ke[2,2] = Ke1[2,2]
Ke[2,3] = Ke1[2,3]
#要素 ① による力 Fy1 と変移 Ux0,Uy0, Ux1,Uy1 の関係 Ke1 から Ke に移す
Ke[3,0] = Ke1[3,0]
Ke[3,1] = Ke1[3,1]
Ke[3,2] = Ke1[3,2]
Ke[3,3] = Ke1[3,3]
//@@@
< 700, 0, -700, 0, 0, 0 >
< 0, 0, 0, 0, 0, 0 >
< -700, 0, 700, 0, 0, 0 >
< 0, 0, 0, 0, 0, 0 >
< 0, 0, 0, 0, 0, 0 >
< 0, 0, 0, 0, 0, 0 >
Ke2, Ke3 を作るために回転を行わせる sf サブルーチン・ファイルを Rotate.sf を作ります。
//@@
@@subFile Rotate # _prm1 is angular of rotation
φ @= _prm{1}
_rt = [[4]]
_rt[0,*] =
_rt[1,*] =
_rt[2,*] = <0, 0, !cos(φ), -!sin(φ)>
_rt[3,*] = <0, 0, !sin(φ), !cos(φ)>
@@endSubFile
//@@@
この Rotate 関数を使って、下の様に ② と ③ の要素剛性行列 KKe2, Ke3 を作ります。
Ke2 の計算;; L@=4000, φ=`π/2, Ke2= @Rotate(φ) E A/L Mat @Rotate(φ)^-1
< 7.6955e-029, -2.01001e-013, -7.6955e-029, 2.01001e-013 >
< -2.01001e-013, 525, 2.01001e-013, -525 >
< -7.6955e-029, 2.01001e-013, 7.6955e-029, -2.01001e-013 >
< 2.01001e-013, -525, -2.01001e-013, 525 >
Ke3 の計算;; L@=5000, φ@=!acos(3/5), Ke3= @Rotate(φ) E A/L Mat @Rotate(φ)^-1
< 151.2, 201.6, -151.2, -201.6 >
< 201.6, 268.8, -201.6, -268.8 >
< -151.2, -201.6, 151.2, 201.6 >
< -201.6, -268.8, 201.6, 268.8 >
Ke1,Ke2 を Ke の節店 2,3 に対応する個所、節店 0,3 に対応する個所に足し合わせます。
//@@
##要素 ② の変移(節店 1,2)と力の関係を Ke に足し合わせる
Ke[2,2] = Ke[2,2] + Ke2[0,0]
Ke[2,3] = Ke[2,3] + Ke2[0,1]
Ke[2,4] = Ke[2,4] + Ke2[0,2]
Ke[2,5] = Ke[2,5] + Ke2[0,3]
Ke[3,2] = Ke[3,2] + Ke2[1,0]
Ke[3,3] = Ke[3,3] + Ke2[1,1]
Ke[3,4] = Ke[3,4] + Ke2[1,2]
Ke[3,5] = Ke[3,5] + Ke2[1,3]
Ke[4,2] = Ke[4,2] + Ke2[2,0]
Ke[4,3] = Ke[4,3] + Ke2[2,1]
Ke[4,4] = Ke[4,4] + Ke2[2,2]
Ke[4,5] = Ke[4,5] + Ke2[2,3]
Ke[5,2] = Ke[5,2] + Ke2[3,0]
Ke[5,3] = Ke[5,3] + Ke2[3,1]
Ke[5,4] = Ke[5,4] + Ke2[3,2]
Ke[5,5] = Ke[5,5] + Ke2[3,3]
//@@@
< 700, 0, -700, 0, 0, 0 >
< 0, 0, 0, 0, 0, 0 >
< -700, 0, 700, -2.01001e-013, -7.6955e-029, 2.01001e-013 >
< 0, 0, -2.01001e-013, 525, 2.01001e-013, -525 >
< 0, 0, -7.6955e-029, 2.01001e-013, 7.6955e-029, -2.01001e-013 >
< 0, 0, 2.01001e-013, -525, -2.01001e-013, 525 >
//@@
##要素 ③ の変移(節店 0,2)と力の関係を Ke に足し合わせる
Ke[0,0] = Ke[0,0] + Ke3[0,0]
Ke[0,1] = Ke[0,1] + Ke3[0,1]
Ke[0,4] = Ke[0,4] + Ke3[0,2]
Ke[0,5] = Ke[0,5] + Ke3[0,3]
Ke[1,0] = Ke[1,0] + Ke3[1,0]
Ke[1,1] = Ke[1,1] + Ke3[1,1]
Ke[1,4] = Ke[1,4] + Ke3[1,2]
Ke[1,5] = Ke[1,5] + Ke3[1,3]
Ke[4,0] = Ke[4,0] + Ke3[2,0]
Ke[4,1] = Ke[4,1] + Ke3[2,1]
Ke[4,4] = Ke[4,4] + Ke3[2,2]
Ke[4,5] = Ke[4,5] + Ke3[2,3]
Ke[5,0] = Ke[5,0] + Ke3[3,0]
Ke[5,1] = Ke[5,1] + Ke3[3,1]
Ke[5,4] = Ke[5,4] + Ke3[3,2]
Ke[5,5] = Ke[5,5] + Ke3[3,3]
//@@@
< 851.2, 201.6, -700, 0, -151.2, -201.6 >
< 201.6, 268.8, 0, 0, -201.6, -268.8 >
< -700, 0, 700, -2.01001e-013, -7.6955e-029, 2.01001e-013 >
< 0, 0, -2.01001e-013, 525, 2.01001e-013, -525 >
< -151.2, -201.6, -7.6955e-029, 2.01001e-013, 151.2, 201.6 >
< -201.6, -268.8, 2.01001e-013, -525, 201.6, 793.8 >
上の操作で ① ② ③ 三つの棒の要素剛性行列を全て足し合わせたので、全体の剛性行列 Ke が求まりました。上の行列要素の値を見れば解かるように Ke は対称行列です。Ke の計算に誤りがないか確認の意味で Ke の固有値を Jacobi.exe で計算させます
# コンソール・コマンド, kShell を使っていれば ctrl O + S
Jacobi Ke
< -1.53392e-013, -8.90445e-014, 1.13341e-014, 478.317, 1205.28, 1606.41 >
三個の固有値が実質的に 0 です。これは、二次元の平行移動と回転の三自由度分だけ変移が不定になる物理的な状況を反映しています。途中で操作ミスが紛れ込むと、剛性行列の個誘値が 0 でなくなります。Jacobi.exe を使って剛性行列の固有値を求めることは、計算ミスが入りこんでいないことを確認する簡単で有効な方法です。
Ke が求められたことで下の式を使って変移から力が求められます。
力 ベクタ = Ke * 変移ベクタ
でも、与えられた問題は、力ベクタが与えられているときの変移を求める問題です。しかも Ke には逆行列が存在しません。このままでは、力ベクタから変移ベクタを求められません。
二次元の平行移動と回転の三自由度分が Ke の行列式を 0 にしています。ならば、三つの自由度をなくしてやれば、固有値から 0 がなくなるはずです。
節点 0 を固定し、節点 1 をローラー固定します。そのときの変移と力をさだめる Ke は次のような K に変形されます。
//@@
K=Ke
e1 = <<6>>
e1[0]=1
# x0 の固定
K[0,*] = e1
K[*,0] = e1
# y0 の固定
K[1,*] = ~shift(e1,1)
K[*,1] = ~shift(e1,1)
# y1 の固定
K[3,*] = ~shift(e1,3)
K[*,3] = ~shift(e1,3)
//@@@
< 1, 0, 0, 0, 0, 0 >
< 0, 1, 0, 0, 0, 0 >
< 0, 0, 700, 0, -7.6955e-029, 2.01001e-013 >
< 0, 0, 0, 1, 0, 0 >
< 0, 0, -7.6955e-029, 0, 151.2, 201.6 >
< 0, 0, 2.01001e-013, 0, 201.6, 793.8 >
!det(K)
< 5.5566e+007 >
上のように逆行列の存在する剛性行列 K を作れました。K を使って節点 2 に x 軸方向に 2000 kgf の力を加えたときの変移を求めます
displacement = K^-1 <0,0, 0,0, 2000,0>
< 0, 0, 1.45851e-015, 0, 20, -5.07937 >
上の計算結果より、2000 kg重の力が働いたことで、 節点 2 が x 方向に 20mm, y 方向に -50 mm 変移していることが解かります。節点 1 が x 方向に 1.5e-15 だけ変移していますが、これは浮動小数点計算の誤差です。1 に比べて 15 桁も小さい値であり、無視できます。
逆に変移 displacement から節店 0,1,2 に加わる力を求められます
force = Ke displacement
< -2000, -2666.67, -2.01948e-028, 2666.67, 2000, 9.09495e-013 >
->
①・②・③ 要素剛性行列 Ke1, Ke2, Ke3 と 0・1・2 節点の変移 displacement sf 変数として残っています。これらを使って、棒 ③ の応力を Ke3 と 0 と 2 節点の変移から計算できます。
Ke3 < displacement[0], displacement[1], displacement[4],displacement[5]>
< -2000, -2666.67, 2000, 2666.67 >
->
他の計算も試してください。加える力や固定点を変えてみて、納得できるまで応力や変移の数値実験を行ってみてください。sf を使えば容易ですから。
ここまでの議論は、より多い数の棒の組み合わせに拡張できます。、三つの棒でできた Trus に限定されません。
L34③ 2
3 ／|→ Fx3 == 2000kgf
L14 ④＼ |②
／ ＼| 4000mm
0-------+1
△ ① ◎
3000mm
L14 = 3000 !cos(!acos(3000/5000))
< 1800 >
L34 = 5000-L14
< 3200 >
//@@
Node = [[4,2]]
Node[<0,4*2,1>] = <\
0 , 0, # 節点 0
3000 , 0, # 節点 1
3000 , 4000, # 節点 2
1800 !cos(!acos(3000/5000)), 1800 !sin(!acos(3000/5000))> # 節点 3
Rod = [[5,2]]
Rod[<0,5*2,1>] = <\
0, 1, # 節点 0 と 1 を結ぶ棒
1, 2, # 節点 1 と 2 を結ぶ棒
2, 3, # 節点 2 と 3 を結ぶ棒
3, 1, # 節点 3 と 1 を結ぶ棒
3, 0> # 節点 3 と 0 を結ぶ棒
//@@@
//@@
@@subFile Rotate # r1,r2,L より決まる
_rt = [[4]]
_rt[<0*4+0,2,1>] = <(r2[0]- r1[0])/L, -(r2[1]- r1[1])/L>
_rt[<1*4+0,2,1>] = <(r2[1]- r1[1])/L, (r2[0]- r1[0])/L>
_rt[<2*4+2,2,1>] = <(r2[0]- r1[0])/L, -(r2[1]- r1[1])/L>
_rt[<3*4+2,2,1>] = <(r2[1]- r1[1])/L, (r2[0]- r1[0])/L>
@@endSubFile
@@subFile SetKatKe # p,q より決まる
Ke[< 2 p N + 2p, 2, 1>] = Ke[< 2 p N + 2p, 2, 1>] +K[<0*4+0, 2,1>]
Ke[< 2 p N + 2q, 2, 1>] = Ke[< 2 p N + 2q, 2, 1>] +K[<0*4+2, 2,1>]
Ke[<(2 p+1) N + 2p, 2, 1>] = Ke[<(2 p+1) N + 2p, 2, 1>]+K[<1*4+0, 2,1>]
Ke[<(2 p+1) N + 2q, 2, 1>] = Ke[<(2 p+1) N + 2q, 2, 1>]+K[<1*4+2, 2,1>]
Ke[< 2 q N + 2p, 2, 1>] = Ke[< 2 q N + 2p, 2, 1>] +K[<2*4+0, 2,1>]
Ke[< 2 q N + 2q, 2, 1>] = Ke[< 2 q N + 2q, 2, 1>] +K[<2*4+2, 2,1>]
Ke[<(2 q+1) N + 2p, 2, 1>] = Ke[<(2 q+1) N + 2p, 2, 1>]+K[<3*4+0, 2,1>]
Ke[<(2 q+1) N + 2q, 2, 1>] = Ke[<(2 q+1) N + 2q, 2, 1>]+K[<3*4+2, 2,1>]
# Dummy return value. Sub file demands _rt value
_rt=0
@@endSubFile
//@@@
//@@
NN@=4
N @= NN * 2 # 節点数 x 節点次元数
Ke @= [[N]] # temporary 変数にしておかないと何時間もの計算時間になってしまう
<<0,5,1@j|
p = Rod[j,0]
q = Rod[j,1]
r1 @= Node[p,*]
r2 @= Node[q,*]
L @= !norm(r2-r1)
K = E A/L @Rotate() Mat @Rotate()^-1
@SetKatKe()
>>
Ke{8} = Ke
//@@@
< 1120, 560, -700, 0, 0, 0, -420, -560 >
< 560, 746.667, 0, 0, 0, 0, -560, -746.667 >
< -700, 0, 1260, -420, 0, 0, -560, 420 >
< 0, 0, -420, 840, 0, -525, 420, -315 >
< 0, 0, 0, 0, 236.25, 315, -236.25, -315 >
< 0, 0, 0, -525, 315, 945, -315, -420 >
< -420, -560, -560, 420, -236.25, -315, 1216.25, 455 >
< -560, -746.667, 420, -315, -315, -420, 455, 1481.67 >
# コンソール・コマンド, kShell を使っていれば ctrl O + A
Jacobi Ke{8}
< -1.60424e-013, 3.85179e-014, 1.49755e-013, 328.844, 1004.08, 1185.44, 2306.98, 3020.49 >
今度も、全体の剛性行列 Ke{8} の固有値は 3 つが 0 になっています。誤りはなさそうです。
節店 0 を固定し、節点 1 をローラー固定します。そのときの剛性行列 K を下の様に計算します。
//@@
K=Ke{8}
e1 = <<8>>
e1[0]=1
# x0 の固定
K[0,*] = e1
K[*,0] = e1
# y0 の固定
K[1,*] = ~shift(e1,1)
K[*,1] = ~shift(e1,1)
# y1 の固定
K[3,*] = ~shift(e1,3)
K[*,3] = ~shift(e1,3)
!det(K)
< 5.5566e+007 >
//@@@
節点 2 に x 軸方向に 2000 kgf の力を加えたときの変移を求めます
displacement = K^-1 <0,0, 0,0, 2000,0, 0,0>
< 0, 0, 2.64339e-016, 0, 20,
-5.07937, 1.71429, 2.28571 >
新たに追加した節点 3 にかかる力は 0 です。節点 2 の位置の変移は計算例 1 と同じです。斜めの棒を追加しても、この場合は剛性には影響しないことになりました。
例題 2 で行った Node と Rod を使った計算のやりかたは汎用性を持ちます。ここまで出来上がっていれば 一般的な truss 有限要素法のプログラムを C などで作ることも容易でしょう。 sf には行列の操作や LU 分解ルーチンやを template library として実装してあります。これらを利用して是非と汎用的な trust 有限要素法プログラムの作成にチャレンジしてみてください。
自転車のホイールを例題とする、もう少し本格的な有限要素法の例題に挑戦します。sf が十分に実用的であると納得してもらえるはずです。
下の様に 32 hole のホイールを考えます。平面に単純化して考えます。
//@@
Node = [[33,2]]
Node[0,*] = <0,0>
r@=622/2, <<0,32,1@j| θ@= 2`π/32, Node[<(j+1)*2,2,1>] = >>
//@@@
また Node 間の接続を既定する Rod.val ファイル変数を下のように作ります。
//@@
Rod = [[32 * 2,2]]
Rod[<0,32*2 * 2,1>] = <\
0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 8, 0, 9,
0,10, 0,11, 0,12, 0,13, 0,14, 0,15, 0,16, 0,17, 0,18, 0,19,
0,20, 0,21, 0,22, 0,23, 0,24, 0,25, 0,26, 0,27, 0,28, 0,29,
0,30, 0,31, 0,32,
1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9,10,
10,11, 11,12, 12,13, 13,14, 14,15, 15,16, 16,17, 17,18, 18,19, 19,20,
20,21, 21,22, 22,23, 23,24, 24,25, 25,26, 26,27, 27,28, 28,29, 29,30,
30,31, 31,32, 32, 1>
//@@@
//@@
/s # 出力を抑制します。計算を早くするためです。
NN@=33
N @= NN * 2 # 節点数 x 節点次元数
Ke @= [[N]] # temporary 変数にしておかないと何時間もの計算時間になってしまう
# スポーク部分の剛性行列の足し合わせ
# `gH は < 9.80665 > m/sec^2 重力加速度です
E = 2e+5/`gH # 2e+5 N/mm^2 ステンレス ついでに引張り強さ 540N/mm^2
A = 2 # mm^2 14 番線
<<0,32,1@j|
p = Rod[j,0]
q = Rod[j,1]
r1 @= Node[p,*]
r2 @= Node[q,*]
L @= !norm(r2-r1)
K = E A/L @Rotate() Mat @Rotate()^-1
@SetKatKe()
>>
# リム部分の剛性行列の足し合わせ
# リムの重量 450g 直径 622mm, アルミの比重 0.0027g/mm^3 より断面積を計算する
r@=622/2, A = 450 /(2`π r 0.0027)
E = 7e+4/`gH # 7e+4 N/mm^2 アルミ
<<32,32,1@j|
p = Rod[j,0]
q = Rod[j,1]
r1 @= Node[p,*]
r2 @= Node[q,*]
L @= !norm(r2-r1)
K = E A/L @Rotate() Mat @Rotate()^-1
@SetKatKe()
>>
Ke{33} = Ke
//@@@
# コンソール・コマンド, kShell を使っていれば ctrl O + A
Jacobi Ke{33}
< -1.32508e-012, -2.98504e-013, 2.65126e-012, 103.816, 103.816, 118.231, 118.231,
123.99, 123.99, 126.785, 126.785, 128.335, 128.335, 129.276,
129.276, 129.885, 129.885, 130.297, 130.297, 130.585, 130.585,
130.789, 130.789, 130.934, 130.934, 131.036, 131.036, 131.103,
131.103, 131.141, 131.141, 131.153, 514.911, 726.969, 726.969,
1902.17, 1902.17, 2262.78, 2262.78, 3697.91, 3697.91, 6128.21,
6128.21, 9093.76, 9093.76, 12478.8, 12478.8, 16152.5, 16152.5,
19973.4, 19973.4, 23794.4, 23794.4, 27468.8, 27468.8, 30855.1,
30855.1, 33823.3, 33823.3, 36259.3, 36259.3, 38069.4, 38069.4,
39184, 39184, 39560.4 >
節点番号 25 を固定点とします。節店番号 0 即ちホイールの中心 x 軸を固定します。ホイールが床と固定されて、ハブが上下にのみ動ける状態です。
//@@
K=Ke{33}
e1 = <<33*2>>
e1[0]=1
# 節店番号 25 の固定
K[25*2,*] = ~shift(e1,25*2)
K[*,25*2] = ~shift(e1,25*2)
K[25*2+1,*] = ~shift(e1,25*2+1)
K[*,25*2+1] = ~shift(e1,25*2+1)
# 節店番号 0 の x 軸固定
K[0*2,*] = ~shift(e1,0*2)
K[*,0*2] = ~shift(e1,0*2)
//@@@
自転車のホイールはニップルを回してスポークとリムの間に張力を与えます。その張力を初期応力 iniTension として下のように設定します。一本のスポークに 20 kgf の張力を与えた状態です。
iniTension = <<33*2>>
iniTension[<2,32*2,1>] = (f@=20, <<0,32,2`π/32@θ| = <-f!cos(θ), -f!sin(θ)> >>)[<0,32*2,1>]
//@@
iniTension = <<33*2>>
f@=20, <<0,32,1@j| θ@= 2`π/32, iniTension[<(j+1)*2,2,1>] = <-f!cos(j θ), -f!sin(j θ)> >>
//@@@
このときの変移 iniDisplacement は下のように計算されます
iniDisplacement = K^-1 iniTension
< 0, -0.0388416, -0.0388416, -0.0388416, -0.0380953, -0.0464193,
-0.035885, -0.0537057, -0.0322956, -0.0604209, -0.0274652, -0.0663068,
-0.0215793, -0.0711373, -0.0148641, -0.0747266, -0.00757763, -0.076937,
-2.23779e-015, -0.0776833, 0.00757763, -0.076937, 0.0148641, -0.0747266,
0.0215793, -0.0711373, 0.0274652, -0.0663068, 0.0322956, -0.0604209,
0.035885, -0.0537057, 0.0380953, -0.0464193, 0.0388416, -0.0388416,
0.0380953, -0.031264, 0.035885, -0.0239776, 0.0322956, -0.0172624,
0.0274652, -0.0113765, 0.0215793, -0.006546, 0.0148641, -0.00295664,
0.00757763, -0.000746331, 3.67382e-015, 20, -0.00757763, -0.000746331,
-0.0148641, -0.00295664, -0.0215793, -0.006546, -0.0274652, -0.0113765,
-0.0322956, -0.0172624, -0.035885, -0.0239776, -0.0380953, -0.031264 >
# Node 25 の y 方向の張力の値 20 自体が、固定点の変移として出てきてしまうので 0 にします
iniDisplacement[25*2+1] = 0
80 kg の体重の半分がホイールにかかったときの変移 loadDisplacement は下の様に計算されます
//@@
load @= <<33*2>>
load[1] = -40
loadDisplacement = K^-1 load
//@@@
< 0, -0.297884,
0.00710324, -0.282864, 0.0044029, -0.283608, 0.00225181, -0.284758, 0.000690955, -0.286132,
-0.000287375, -0.287551, -0.000736685, -0.288848, -0.000750713, -0.289884, -0.000455565, -0.290551,
-7.66464e-016, -0.29078, 0.000455565, -0.290551, 0.000750713, -0.289884, 0.000736685, -0.288848,
0.000287375, -0.287551, -0.000690955, -0.286132, -0.00225181, -0.284758, -0.0044029, -0.283608,
-0.00710324, -0.282864, -0.0102631, -0.282697, -0.0137471, -0.283257, -0.0173795, -0.284659,
-0.0209531, -0.286976, -0.0242395, -0.29023, -0.0270013, -0.294388, -0.0290058, -0.299356,
0, 0, 0.0290058, -0.299356, 0.0270013, -0.294388, 0.0242395, -0.29023,
0.0209531, -0.286976, 0.0173795, -0.284659, 0.0137471, -0.283257, 0.0102631, -0.282697 >
ノードの配置データと変移データである sf 変数、 Node : loadDisplacement と iniDisplacement の二つの変移を gnuplot で可視化します。gnuplot で二つのデータを重ね合わせて表示させるために、エディタを使った手作業が必要となります。
まず最初に iniDisplacement 変移分だけ変化したノードの cdspl を求めます。下の sf 式で計算します。変移量は小さくそのまま表示しても判別できないので、30 倍の大きさに拡大して表示します。
//@@
cdspl = <<33>>
r@=30, <<0,33,1,@I| cdspl[I] = Node[I,0] + r iniDisplacement[2I] + (Node[I,1]+r iniDisplacement[2I+1])i>>
//@@@
コンソールで cdspl データの x 成分のみを表示させます。
gdsp cdspl
カレント ディレクトリの gnuplot.ini ファイルには下のデータ部分が書かれています。y 成分も書かれています。
# gnuplot.ini にできているデータ
0 0 -1.1652
1 309.83 -1.1652
2 303.88 59.281
3 286.25 117.4
4 257.62 170.97
5 219.09 217.92
6 172.13 256.45
7 118.57 285.08
8 60.446 302.72
9 -4.8091e-014 308.67
10 -60.446 302.72
11 -118.57 285.08
12 -172.13 256.45
13 -219.09 217.92
14 -257.62 170.97
15 -286.25 117.4
16 -303.88 59.281
17 -309.83 -1.1652
18 -303.88 -61.611
19 -286.25 -119.73
20 -257.62 -173.3
21 -219.09 -220.25
22 -172.13 -258.78
23 -118.57 -287.42
24 -60.446 -305.05
25 5.3087e-014 -311
26 60.446 -305.05
27 118.57 -287.42
28 172.13 -258.78
29 219.09 -220.25
30 257.62 -173.3
31 286.25 -119.73
32 303.88 -61.611
このコンソー gnplot.ini 表示データの "0 0 -1.1652" 以下の行をコピーして wheel.data にペーストします。
同様に 80 Kg の負荷がかかったときの変移 loadDisplacement の可視化データを作成します。下の sf 式で計算します。ここでも変移量を、30 倍に拡大しています。
//@@
cdspl = <<33>>
r@= 30, <<0,33,1,@I| cdspl[I] = Node[I,0] + r (loadDisplacement[2I]+iniDisplacement[2I]) + (Node[I,1]+r (loadDisplacement[2I+1]+iniDisplacement[2I+1]) )i>>
//@@@
ここでも gdsp cdspl によって gnuplot.ini に出力された数値データを wheel.dat にコピー・ペーストします。ただし、二つのプロックの間に空行を挿入します。gnuplot で二つのデータを重ねて表示するためです。
gdsp cdsp
iniDisplacement と loadDisplacement 二つの変移を重ね合わせた Node の可視化データ wheel.dat を下のように作ります。ここで二つのデータ・ブロックの間に二行以上の改行が必要です。一個の改行だけではエラーになります。ご注意ください
# 最終的な wheel.data
type wheel.data
0 0 -1.1652
1 309.83 -1.1652
2 303.88 59.281
3 286.25 117.4
4 257.62 170.97
5 219.09 217.92
6 172.13 256.45
7 118.57 285.08
8 60.446 302.72
9 -4.8091e-014 308.67
10 -60.446 302.72
11 -118.57 285.08
12 -172.13 256.45
13 -219.09 217.92
14 -257.62 170.97
15 -286.25 117.4
16 -303.88 59.281
17 -309.83 -1.1652
18 -303.88 -61.611
19 -286.25 -119.73
20 -257.62 -173.3
21 -219.09 -220.25
22 -172.13 -258.78
23 -118.57 -287.42
24 -60.446 -305.05
25 5.3087e-014 -311
26 60.446 -305.05
27 118.57 -287.42
28 172.13 -258.78
29 219.09 -220.25
30 257.62 -173.3
31 286.25 -119.73
32 303.88 -61.611
0 0 -10.102
1 310.05 -9.6512
2 304.01 50.772
3 286.32 108.86
4 257.64 162.39
5 219.08 209.29
6 172.11 247.79
7 118.55 276.39
8 60.432 294
9 -7.1085e-014 299.95
10 -60.432 294
11 -118.55 276.39
12 -172.11 247.79
13 -219.08 209.29
14 -257.64 162.39
15 -286.32 108.86
16 -304.01 50.772
17 -310.05 -9.6512
18 -304.19 -70.092
19 -286.66 -128.23
20 -258.14 -181.84
21 -219.71 -228.86
22 -172.86 -267.49
23 -119.38 -296.25
24 -61.316 -314.03
25 5.3087e-014 -311
26 61.316 -314.03
27 119.38 -296.25
28 172.86 -267.49
29 219.71 -228.86
30 258.14 -181.84
31 286.66 -128.23
32 304.19 -70.092
最後に、重ねた表示を行わせるためにコマンド下のような gnuplot.ini ファイルを作り、 gnuplot を起動します。
type gnuplot.ini
plot "wheel.data" index 0 using 2:3,\
"wheel.data" index 1 using 2:3
start wgnuplot
問題がおきるのはホイールの中心を固定点としたときです。即ち節点番号 0 を固定点とします。そして節点番号 1 の x 軸を固定します。ハブ軸が固定されて、回転しない状態です。
//@@
K=Ke{33}
e1 = <<33*2>>
e1[0]=1
# 節店番号 0 の固定
K[ 0*2,*] = ~shift(e1, 0*2)
K[*, 0*2] = ~shift(e1, 0*2)
K[ 0*2+1,*] = ~shift(e1, 0*2+1)
K[*, 0*2+1] = ~shift(e1, 0*2+1)
# 節店番号 1 の x 軸固定
K[1*2,*] = ~shift(e1,1*2)
K[*,1*2] = ~shift(e1,1*2)
//@@@
でも、このときの剛性行列 K の固有値には、上のように 0 が入ってきます。逆行列が作れません。
Jacobi K
< 2.37243e-013, 1, 1, 1, 59.9609, 61.8562, 103.816, 105.18,
118.231, 118.993, 123.99, 124.453, 126.785, 127.092, 128.335, 128.552,
129.276, 129.435, 129.885, 130.006, 130.297, 130.391, 130.585, 130.658,
130.789, 130.846, 130.934, 130.978, 131.036, 131.068, 131.103, 131.123,
131.141, 131.15, 501.212, 805.571, 831.335, 1877.91, 1902.17, 3675.4,
3697.91, 6107.5, 6128.21, 9075.01, 9093.76, 12462.2, 12478.8, 16138.2,
16152.5, 19961.4, 19973.4, 23784.8, 23794.4, 27461.4, 27468.8, 30849.8,
30855.1, 33819.9, 33823.3, 36257.3, 36259.3, 38068.5, 38069.4, 39183.8,
39184, 39560.4 >
ちなみに、このぐらいのサイズになると逆行列の有無の判定に行列式の値は使えません。行列式の値は下のように大きな値になります。一つの固有値が 1e-13 の小ささでも、他の要素が 10^3 のオーダーならば、たくさんの要素の積を計算したときには大きな値になってしまうからです。
!det(K)
< -1.72551e+181 >
変です。三つの自由度を拘束したのですから、K の固有値に 0 が入り込むのも不思議なことです。－－－－－－－ 誤りは "節点番号 1 の x 軸を固定します。ハブ軸が固定されて、回転しない状態です" とした所にありました。"節点番号 1 の y 軸を固定します。ハブ軸が固定されて、回転しない状態です" とするべきだったのです。実際、このときの剛性行列の固有値に 0 は入ってきません。
//@@
K=Ke{33}
e1 = <<33*2>>
e1[0]=1
# 節店番号 0 の固定
K[ 0*2,*] = ~shift(e1, 0*2)
K[*, 0*2] = ~shift(e1, 0*2)
K[ 0*2+1,*] = ~shift(e1, 0*2+1)
K[*, 0*2+1] = ~shift(e1, 0*2+1)
# 節店番号 1 の x 軸固定
K[1*2+1,*] = ~shift(e1,1*2+1)
K[*,1*2+1] = ~shift(e1,1*2+1)
//@@@
# コンソール実行
Jacobi K
< 1, 1, 1, 22.1727, 59.9609, 89.7092, 103.816, 113.837, 118.231, 122.264, 123.99, 125.984,
126.785, 127.919, 128.335, 129.044, 129.276, 129.75, 129.885, 130.217, 130.297, 130.537, 130.585, 130.761,
130.789, 130.918, 130.934, 131.028, 131.036, 131.099, 131.103, 131.14, 131.141, 131.153, 514.911, 561.285,
831.335, 1245.31, 1902.17, 2691.68, 3697.91, 4817.95, 6128.21, 7531.17, 9093.76, 10724, 12478.8, 14272.7,
16152.5, 18040.5, 19973.4, 21882.3, 23794.4, 25650.5, 27468.8, 29200.1, 30855.1, 32394.7, 33823.3, 35111.5,
36259.3, 37246.1, 38069.4, 38716.5, 39184, 39466 >
考え直してみれば節点 1 の x 軸方向の自由度だけを拘束しても、まだホイールには回転の自由度が残っていました。ですから三つの拘束条件を与えただけでは剛性行列 K の固有値は 0 になる場合があることになります。やっぱり物理的直感も働かせてやらないと失敗します。
このことは、下手に三つの拘束条件を与えるだけでは、剛性行列の逆行列が信用できないことがあることを意味します。ホイールのような規則的で単純なときは、逆行列が計算できなくなって変なことが判りました。でも、なまじ逆行列が計算できてしまうような規則的ではない配置のときは見逃してしまったはずです。拘束条件が本当に物理的に拘束できているか判定しなければならない例として書いておきます。