「テーマ」

　フラッタ解析に遠心力を考慮するためまずはScilabを用いて1次元の弦の解析にトライしてみます。

「方針」

　1次元梁のモデルに張力を設定して弦の運動を求めてみます。

---

1．解析モデル

１．１　1次元梁モデル

　梁の曲げ方程式は以下となります(参考文献[1]3.1章より)。

  　　　EId⁴ω/dｘ⁴ = 0

　　　　E : 梁要素の縦弾性係数[Pa]

　　　　I：梁の断面2次モーメント[m⁴]

　　　　L：梁のx方向の微小長さ[m]

　　　　ω：梁のz方向の変位[m]

　　　　β：梁の傾き-ｄω/dx(時計回りを＋)[rad]

　ここでは細部を割愛しますが、板の解析で示した方法と同様にして上図の一つの要素の力の釣り合い式は以下となります。

一般化すると

　　{F}={K} x {δ}

　　{F}={F₁、M₁ 、F₂、M₂}^T　　：境界でのせん断及び曲げモーメント

　　{δ}={ω₁、β₁ 、ω₂、β₂}^T   　： 境界条件

　　{K}  　　　　　　　　　　：剛性行列と呼びます

　また質量に対応する整合質量行列は以下となります

　(整合質量行列の詳細は[Detail1:変位法による薄板FEM解析詳細]に示します)。

　　　　ρ: 梁の密度[kg/m³]

　　　　A：梁の断面積[m²]

　　　　L：梁のx方向の微小長さ[m]

１．２　1次元弦モデル

　弦の運動方程式は以下となります(参考文献[2]7.1章より)。

  　　　ρd²ω/dt² = Td²ω/dx²

　　　　ρ: 弦の線密度[kg/m]

　　　　T：弦への張力[N]

　　　　ω：弦のz方向の変位[m]

　張力Tの分力が弦要素のｚ方向運動を担います。

　図示した張力Tの傾きにより発生するz方向の分力Tβを梁の運動方程式に代入してみます(当事業所オリジナルアイデアのため検証が重要となります)。

１．３　梁の弦のモデル統合

　例として3点(下図)で具体的計算を実施してみます。

　まずは梁の運動方程式を羅列します。

　この2式をまとめると以下となります。

　同様に整合質量行列は以下となります。

　最後に張力を表現する行列は以下となります。

張力のｚ方向分力は上記のような足し算をすると打ち消してしまうので以下のように半分ずつに分けてみます。

　以上の行列をScialbに用いることで梁と弦の結合モデル(下図)が出来上がりました。

　尚両端ω₁、ω₂を拘束する初期条件とした場合、対応する列と行(黄色部)を除外する必要があります。

（[MI]は[M]の逆行列です）

---

２．計算その１

２．１　1次元計算モデル(両端支持)

　具体的な計算を行い、理論式と比較することでモデルの妥当性を確認してみます。

　以下のモデルを考えて行きます。

・材料：鉄(E=2.1e11[Pa]、ρ=7800[kg/m³])

・寸法：直径1mm×50cm

・分割数：8要素9節点

・両端支持：ω₁=ω₂=0

・張力：3[N]

　Scilabモデルは以下となります。

　　・積分方式：ルンゲクッタ

　　・積分幅：1e-5[s]

　　・積分時間：3[s] データ抽出はラスト2[s]

　　・初期位置：非対称矩形(適当)を設定

　　尚モデル安定のため以下を追加。

　　・構造減衰模擬行列[KMsqrt]（減衰率ξ=0.001）：ノイズカット

　　　構造減衰模擬行列については [Detail2:変位法による薄板FEM動的解析詳細]を参照方　　

　　・張力行列[Nstr]を1秒間フェードインするためにランプブロック及びリミッタブロック

２．２　計算結果

　代表として節点2～5の動きを計算結果を以下に示します。

減衰力のおかげて不必要な高周波成分がカットされていきます。

　更に本計算結果のうちラスト2秒のFFT結果を示します。

　FFT結果から23Hzと53Hz、94Hzに大きなピークが出来ています。

　エクセルのFFT関数を使い、それぞれの周波数成分を抽出し、2.98秒から3秒までの20msecを表示してみました。

　結果、それぞれが波の１次、２次、３次の固有形状に対応していることが確認出来ました。

　ここで本結果を理論値と比較してみます。

・梁の固有振動数理論値(両端支持条件)

　 fi=1/2π×(λi/ℓ)²×√(EI/(ρA)) 　[関連文書[2]7.3章]

　　　　E : Feヤング率 2.1e11[Pa]

　　　　ρ：Fe密度 7800[kg/m³]

　　　　A：断面積 π*0.0005²[m²]

　　　　I：断面２次モーメント π*0.0005⁴/4[m⁴]

　　　　ℓ：0.5[m]

　 　　 λ₁=π、λ₂=2π、λ₃=3π

　 ∴f₁=8.15Hz、f₂=32.60Hz、f₃=73.35Hz

・弦の固有振動数理論値(両端支持条件)

　 f_i=i/(2*ℓ)×√(T/(ρA)) 　[関連文書[2]7.1章]

　　　　T : 張力3[N]

　　　　ρ：Fe密度 7800[kg/m³]

　　　　A：断面積 π*0.0005²[m²]

　　　　ℓ：0.5[m]

　 　　 i=1、2、3

　 ∴f₁=22.13Hz、f₂=44.26Hz、f₃=66.39Hz

・Scailab計算結果

　 　f₁=23Hz、f₂=53Hz、f₃=94Hz

図示してみます。

結果、梁の固有モードと張力による弦のモードがカップリングして新しい固有振動数のモードになる様子が見て取れます。

今回の計算では前半は弦側の運動が支配的になっているようです。

---

３．計算その２

３．１　1次元計算モデル(片端固定)

　続いて以下のような片端固定の場合を考えてみます。

・材料：鉄(E=2.1e11[Pa]、ρ=7800[kg/m³])

・寸法：直径1mm×50cm

・分割数：8要素9節点

・片端固定：ω₁=β₁=0

・張力：3[N]

３．２　計算結果

　計算方法は2章と同じとしました。3秒の計算のうち、後半の2秒のFFT結果を示します。

　FFT結果から0.5Hzと29Hz、66Hz、115Hzに大きなピークが出来ています。

　エクセルのFFT関数を使い、それぞれの周波数成分を抽出し、2.96秒から3秒までの40msecを表示してみました。

　結果、それぞれが梁の1次、2次、3次、4次の固有形状に対応していることが確認出来ました。

ここで本結果を理論値と比較してみます。

・梁の固有振動数理論値(片端固定条件)

　 f_i=1/2π×(λi/ℓ)²×√(EI/(ρA)) 　[関連文書[2]7.3章]

　 　　 λ₁=1.875、λ₂=4.694、λ₃=7.855、λ₄=10.996

　 ∴f₁=2.90Hz、f₂=18.20Hz、f₃=50.95Hz、f₄=99.85Hz

・弦の固有振動数理論値(両端支持条件:片端固定の弦の公式が見当たりません)

　 f_i=i/(2*ℓ)×√(T/(ρA)) 　[関連文書[2]7.1章]

　 ∴f₁=22.13Hz、f₂=44.26Hz、f₃=66.39Hz

・Scailab計算結果

　 　f₁=0.5Hz、f₂=29Hz、f₃=66Hz、f₄=115Hz

　図示してみます。

　結果、両端支持(2章)の時と違って、カップリングに次数のずれが起きています。カップリングは次数ではなく周波数と波形により選別されると考えられます。

---

４．まとめ

　両端支持(2章)の場合、張力を上げていくことで全体の剛性が上がり、基本的に弦の周波数に梁の剛性が補正される様子を呈していました。

　しかし片端固定(3章)の場合、張力を上げていっても取り残される梁の剛性モードが存在し、しかもその剛性は低下しています。これは工学上、注意が必要と考えます(低下するモードでは弦による荷重と梁による荷重の位相がおよそ反転していることをScialb計算で確認しています)。

「注意」

　本計算は当事業所オリジナルアイデアであり、残念ながら試験値との比較が出来ていません。

　１次モードは問題ないと考えますが2次、3次の挙動が未確認となります(ネット上で対応する試験データを見つけられませんでした)。

「所感」

・FreeCADで同じ計算を試みたのですが、メッシュ形状とソルバーによりモードが大きく変化し、Scilab計算よりも理論値との乖離が顕著に出てしまい、なにが真に近いのか分からなくなってしまいました。むしろ当事業所の実施している簡易モデルの方が現実を理解し易いと感じました。

・梁の固有振動数を議論する時は空間次元(x軸)の要素数はそれほど重要ではなかったのですが、弦の解析では5要素で計算が発散してしまいました。進行波を表現するには空間次元(x軸)もある程度の要素数が必要のようです。今回は8要素で計算しました。多いほど発散しないのですが、行列作成は大変となります。

・張力を上げていくと計算が発散してしまいます。上述と重なりますが要素数を適宜上げる必要があります。

---