「テーマ」
・フラッタ解析に遠心力を考慮するため、1次元の弦に続いて2次元膜の解析にトライしてみたいと思います。
「方針」
２次元薄板モデルに張力を設定して膜の振動を解いてみます。
以降張力を考慮した薄板モデルを膜モデルと記述し、張力なしの薄板モデルと区別します。
1．解析モデル
１．１ 2次元薄板モデル
薄板の曲げ方程式は以下となります(参考文献[1]4.1章より)。
D(∂4ω/∂x4+2∂4ω/∂x2∂y2+∂4ω/∂y4)=0
D=Et3/12/(1-ν2)
E : 板要素の縦弾性係数
ｔ：板厚
ν：板要素のポアソン比
ω：板のz方向の変位
上記方程式より１次元梁と同様に剛性行列[K]及び整合質量行列[M]を求めることが出来ます。
（詳細は「Detail1」変位法による薄板FEM解析詳細に示しました。）
更に張力の分力を表現する行列[Nstr]を作成することで以下のscilabモデルにより膜モデルの振動を解析してみます。
尚薄板の運動に下図のような張力を入れ込む解析方法は当事業所オリジナルアイデアのため検証が重要となります。
１．２ 2次元膜モデル
膜の運動方程式は以下となります(参考文献[2]7.1章より)。
ρd2ω/dt2 = T(d2ω/dx2+d2ω/dy2)
ρ: 膜の面密度[kg/m]
F：膜への外力[N]
T：各辺に加わる張力(Fの1/2とします)
ω：膜のz方向の変位[m]
下図のように張力Tの分力が膜要素のｚ方向運動を担います。
図中Nijを行列で表現すると
N12の半分が点1と点2に2等分されるとします。
これを行列で表現すると
Scilabに用いる張力の分力を表現する行列[Nstr]は以下となります。
尚上述膜の運動方程式にｘとｙのカップリング項がないため、x軸とy軸の張力を独立に扱う事が出来ると考えます。
１．３ 膜のモデル計算例
例として6点(下図)で具体的に張力の行列を計算してみます。
四角①②③④に関しては前項より
四角③④⑤⑥も同様に
2つを合わせること(黄色部分)で6点モデルでの張力を表現する行列が出来上がります。
２．計算その１
２．１ 2次元膜モデル(全周固定)
具体的な計算を行い、理論式と比較することでモデルの妥当性を確認してみます。
以下のモデルを考えて行きます。
・材料：AL(E=7e10[Pa]、ρ=2800[kg/m3]、ポアソン比 0.3)
・寸法：500×500×0.2ｍｍ
・分割数：11×11節点10×10要素
・全周固定：ωi=ωxi=ωyi=0
・張力：面圧1e6[Pa]
Scilabモデルは以下となります。
・積分方式：ルンゲクッタ
・積分幅：1e-5[s]
・積分時間：3[s] データ抽出はラスト2[s]
・初期位置：中心点61を3mm(適当)に設定
・初期荷重：IniF＝[0]
尚モデル安定のため以下を追加。
・構造減衰模擬行列[KMsqrt]（減衰率ξ=0.001）：ノイズカット
・張力行列[Nstr]を1秒間フェードインするためにランプブロック及びリミッタブロック
２．２ 計算結果
計算結果の代表として節点59、61、63の動きを以下に示します。減衰力のおかげて不必要な高周波成分がカットされていきます。
更に膜モデル計算結果のうち中心点61のラスト2秒のFFT結果(500Hzサンプリング)を示します。
FFT結果から27.3Hzと61.0Hz、85.9Hz、98.6Hzに大きなピークが出来ています。
エクセルのFFT関数を使い、それぞれの周波数成分を抽出し、2秒付近の形状を表示してみました。
結果、膜モデルの1次、2次及び3次は薄板固有モード1次、5次及び7次に対応していますが、4次以上は薄板固有モードとは一致していないようです(モデルの粗さが原因と推定されます)。
(コメント：FFTを用いた周波数成分抽出では復元されたデータの終盤時間(ここでは3秒直前)では波形が不安定になるようです)
(フィルター25-30Hz)
(フィルター55-65Hz)
(フィルター80-90Hz)
(フィルター95-105Hz)
ここで本結果を理論値と比較してみます。
・板の固有振動数理論値(全周固定条件)[固定：ω=ωx=ωy=0]
fi=1/2π×(λi2/a2)×√(D/(ρt)) [関連文書[2]7.4章]
D=Et3/12(1-ν2)
E : ALヤング率 7e10[Pa]
ρ：Fe密度 2800[kg/m3]
a：正方形１辺の長さ 0.5[m]
t：板の厚さ 0.2[mm]
ℓ：0.5[m]
ν：ポアソン比 0.3
λ12=35.99、λ22=73.41、λ32=108.3、λ42=131.6、λ52=132.2
∴f1=6.93Hz、f2=14.14Hz、f3=20.86Hz、f4=25.35Hz、f5=25.47Hz
(f6=30.69Hz、f7=40.20Hz、f8=40.21Hz、f9=45.01Hz、f10=45.30Hz、f11=53.73Hz)※1
※1：括弧内は理論値のλ値が見つけられなかったためscilab計算値を示しています
・膜の固有振動数理論値(全周支持条件)[支持：ω=0] (ωxとωyは拘束されていません)
fi=1/2×√(i2/a2+j2/a2)×√(P/ρ) [関連文書[2]7.4章]
P : 張力の元となる面圧 1e6[Pa]
ρ：AL密度 2800[kg/m3]
a：0.5[m]
i、j=1、2、3
∴f1=26.73Hz(i=1、j=1)、f2=42.26Hz(i=2、j=1)、f3=53.45Hz(i=2、j=2)、f4=59.76Hz(i=3、j=1)、 f5=68.14Hz(i=3、j=2)、f6=77.92Hz(i=4、j=1)、f7=80.18Hz(i=3、j=3)、f8=84.52Hz(i=4、j=2)、f9=94.49Hz(i=4、j=3)
・Scailab計算結果 f1=27.3Hz、f2=61.0Hz、f3=85.9Hz、f4=98.6Hz
周波数から推定されるカップリングの様子を図示してみます。
結果、薄板の固有モードと張力による膜のモードがカップリングして新しい固有振動数のモードになる様子が見て取れます。なおこの結果は中心点61でのカップリングであり、他の点では別のカップリングが起きていることを確認しました。非常に複雑に絡み合っています。
３．計算その２
３．１ 2次元計算モデル(1辺固定)
続いて以下のような1辺固定の場合を考えてみます。
・材料：AL(E=7e10[Pa]、ρ=2800[kg/m3]、ポアソン比 0.3)
・寸法：500×500×0.2mm
・分割数：11×11節点10×10要素
・1辺固定：①～⑪ ωi=ωxi=ωyi=0
・張力：面圧1e4[Pa]・初期位置：中心点61を2mm(適当)に設定
３．２ 計算結果
計算方法は2章と同じとしました。12秒の計算のうち、中心点61のラスト10秒のFFT結果(200Hzサンプリング)を示します。
FFT結果から0.49Hzと4.59Hz、5.27Hz、10.64Hzに大きなピークが出来ています。
エクセルのFFT関数を使い、それぞれの周波数成分を抽出し、10秒付近の形状を表示してみました。
結果、膜モデルの1、2、3、4次モードは薄板モデルの1次、3次、4次、6次モードに対応していました。
その形状の様子を示します(上段：膜モデル、下段：薄板モデル)。
ここで本結果を理論値と比較してみます。
・薄板の固有振動数理論値(片端固定条件)[固定：ω=ωx=ωy=0]
fi=1/2π×(λi2/a2)×√(D/(ρt)) [関連文書[2]7.4章]
D=Et3/12(1-ν2)
λ12=3.494、λ22=8.547、λ32=21.44、λ42=27.46、λ52=31.17
∴f1=0.67Hz、f2=1.65Hz、f3=4.13Hz、f4=5.29Hz、f5=6.00Hz
(f6=10.39Hz、f7=11.83Hz、f8=12.37Hz、f9=13.68Hz、f10=17.72Hz、f11=18.53Hz)※1
※1：括弧内は理論値のλ値が見つけられなかったためscilab計算値を示しています
・膜の固有振動数理論値(両端支持条件:片端固定の式が見あたないため、両端固定を使用)
fi=i/(2*ℓ)×√(P/(ρ)) [関連文書[2]7.1章]
P : 張力の元となる面圧 1e4[Pa]
ρ：AL密度 2800[kg/m3]
ℓ：0.5[m]
i、j=1、２、３
∴f1=1.89Hz、f2=3.78Hz、f3=5.67Hz、f4=7.56Hz、f5=9.45Hz
・Scailab計算結果
f1=0.49Hz、f2=4.59Hz、f3=5.27Hz、f4=10.64Hz
周波数から推定されるカップリングの様子を図示してみます。
薄板モードの周波数が変わっていくという様子が伺えます(厳密には新しいモードが作成されているのですが、おおよそ薄板のモードを引き継いでいるようです)。
また薄板で存在していた低周波数のモードが更に低い周波数に移動しています。尚初期位置を変更することで以下のカップリングの様子も変わってきます。
４．まとめ
両端支持(2章)の場合、張力を上げていくことで全体の剛性が上がり、基本的に膜の周波数に梁の剛性が補正される様子を呈していました。
しかし片端固定(3章)の場合、張力を上げていっても取り残される薄板モードが存在し、しかもその周波数は低下しています(工学上、注意が必要)。
この様子は1次元の解析結果と一致しています。
以上より薄板運動に張力を加えることで薄板固有モードの振動数が上がることを確認出来ました。
そうなると遠心力を考慮したフラッタ運動を解析した場合、遠心力が張力として働きフラッタ運動が消えるという解析結果が推定されます。
「注意」 本計算は当事業所オリジナルアイデアであり、残念ながら試験値との比較が出来ていません。
「補足」
・張力のない場合、動的挙動を計算する際の時間積分幅はその系の最大固有振動数をカバーする程度にすれば発散を防げます。具体例として系の最大固有値振動数が1000Hzの場合(2.2項の計算)、おおよそ5倍以上の周波数つまり5000Hz以上の積分幅にすると発散を防ぐことを経験しています(1波形を5点で表現している状態です)。
・弱い張力ならば上記5000Hzでも発散しません(2.2項の計算)。
・ところが張力を4倍にすると積分幅を200倍の1000000Hzにしても発散してしまいます。つまり張力による発散は系の最大固有振動数には関係していないようです。
・そこで以下のような推定を立てました。
① 波数
最小要素の長さをdL(2.2項では5cm)とすると最小の波は2*dLで表現出来ます。すると１辺Lでの波数はk=L/(2*dL)となります(2.2項では50cm/(2*5cm)=5個)。
② 膜の運動方程式からの最大周波数
薄板における膜の周波数は以下の式から求められます。
fi=1/2×√(i2/a2+j2/a2)×√(P/ρ) [関連文書[2]7.4章]
i、j=1、2、3：i、jは半波長を示している
①で示した波数kはi、jと以下の関係となります。
i、ｊ=2*k=L/dL(2.2項では10)
よって2.2項ではfmax=267Hz
③ 薄板固有振動数の最大値と弦の運動の最大周波数
当事業所で設定したscilabモデルはベースが薄板曲げモデルであり、そこに弦の運動方程式を補間したような形となっています(本来は曲げ運動と弦運動を同時に満たす基底を求めるべきです)。
そのような摂動モデルにおいてベースとなる曲げモデルの最大周波数が2.2項では1000Hz。対して弦の運動は267㎐で約4倍。この程度差があることで弦の運動モードを曲げモデルモードが表現出来ているのではないかと考えます。
この状態で張力を4倍にすると弦の最大周波数は2倍の535㎐になってしまい、曲げモデル周波数1000Hzでは表現出来なくなり、scilab計算で発散してしまうと推定します。
尚１次元の弦解析では梁の曲げ最大周波数は弦運動の最大周波数の3倍程度となりました(上述4倍に対応)。
参考文献
[１]機械構造 振動学 MATALBによる有限要素法と応答解析 小松敬冶 著 森北出版株式会社
[２]機械振動学 岩田佳雄・佐伯暢人・小松崎俊彦 共著 数理工学社
HA-technodefy 