「テーマ」
・遠心力を考慮したブレードのフラッタ解析に挑戦してみたいと考えます。
私自身ブレード設計の経験がなく、また設計データや試験データを入手出来なかったためあくまで試算となります。
「方針」
すでにフラッタ計算はScilabを用いた薄板フラッター解析にて検討しており、また張力の影響はScilabを用いた膜の解析にて検討しています。そこで張力を遠心力に置き換え、フラッタ解析に含めることでブレードのフラッタ解析を挑戦してみます。
注意：上述の通り、あくまで試算のため参考程度にご覧ください。
１．解析モデル
１．１ モデル方針
Scilabを用いた薄板フラッター解析にて用いたScilabモデルは以下となります。
特徴は薄板構造モデルに1次曲げ空気力と1次捩り空気力を加えた点となります(薄板構造モデルの詳細は「Detail1」変位法による薄板FEM解析詳細及び「Detail2」変位法による薄板FEM動的解析詳細をご参照下さい)。
このモデルによりフラッタの中でも基本的な1次曲げモードと1次捩りモードによる振動発散を表現しています。
一方Scilabを用いた膜の解析にて用いたScilabモデルは以下となります。
特徴は薄板構造モデルに張力を加えた点となります。
このモデルにより張力による膜振動を表現しています。 今回はこの張力を遠心力とすることで遠心力を考慮したブレードのフラッタ解析を試みます。
１．２ 解析モデル
ブレードとして以下の寸法を想定しました。フラッタが起きやすいようにあえて付け根の捩れ剛性が低くなるように細くしています。
この形状を以下のモデルで解析します。
・300×800×2mmの羽子板形状
・71節点の三角マスで分割・材質：AL材
１．３ Scilabモデル
フラッタ解析に使用するScilabモデルを以下に示します。見た目は非常に簡単なモデルになっていますが、豊かな解析結果を示してくれます。
モデルの構成は以下となります(詳細を3項に示します)。
①薄板構造部
薄板の剛性行列[Kf]と質量行列の逆行列[MIf]により構成されており、本モデルの積分部分となります。
②通常外力(空気力)[Ini_forceRPM]
定常状態(迎角2°)での空気力[Ini_forceRPM]を印加します。
③1次捩り空気力[Air_Mode201]
捩り運動によるばね要素をモデル化します。具体的には1次捩りの形状により発生する空気力をOpenFoamにて事前計算し、[Air_Mode201]として印加します。
④1次曲げダンピング力[Air200_AOA]
曲げ運動による局所速度変化は周囲速度との合成により迎角を変化させます。
⑤遠心力[Nstr]
遠心力により発生する面直方向の荷重を[Nstr]として印加します。 今回はブレードの発生する推力は省略しましたが、この[Nstr］に含めることが可能です。
２．計算結果
Scilab計算時のパラメータは以下となります。
・RPM：回転数[rpm]
・ξ：構造減衰係数 一般的な構造で1％位
・かけ数：[Air_Mode201]及び[Air200_AOA]に作用させる１次遅れ時定数taoの割合でここでは0.5としました(詳しくは3.1項)
まずは遠心力を考慮しない場合の計算結果を示します。図中Point69は節点69の面直方向の変位、Mode200及びMode201はscialbモデル上でのモード存在量を示しています。
〇RPM=380[rpm](遠心力なし)
この回転では発散まではいかないまでも
わずかな持続振動が見られます。
〇RPM=400[rpm](遠心力なし)
急激な拡散振動が見て取れます。
フラッタを表現していると考えます。
回転数400[rpm]でのフラッタの様子を計算時間4秒から0.1秒ごとのコマ送りで示します(左⇒右の順)。１次曲げモードと1次捩れモードが混ざりあっている様子が確認できます。
続いて遠心力を考慮した場合の計算結果を示します。
〇RPM=450[rpm](遠心力あり)
遠心力なしの場合と同じような
拡散振動が見て取れます。
〇RPM=600[rpm](遠心力あり)
振動発散は完全になくなっています。
結果として遠心力を考慮しない場合、400[rpm]でフラッタが発生するのですが、遠心力を考慮すると600[rm]では完全にフラッタが消失しています。しかし遠心力を考慮しても450[rpm]当たりまでは振動発散してしまいます。この辺りを次項にて確認します。
３．計算詳細
３．１ 遠心力なしの場合
３．１．１ Scilab基本モデル
以下の薄板のScilab変位法を用いたMode解析結果を示します。
・300x800x2mmアルミ板
・ヤング率70GPa / 密度2800kg/m3 / ポアソン比0.3
・71節点(1要素は三角形状)で1点辺りの自由度3(ω、ωx、ωy)※１
・端部固定(①②③)により自由度は3点減って、全自由度は68×3=204。よってMode数も204。
Mode解析の結果は以下となりました。
1次曲げモード
1次捩りモード
2次曲げモード
2次捩りモード
このModeのうち、最も低次である1次曲げModeと１次捩りModeを以後のScilabモデルの対象とします。尚ここで用いたMode番号(200～203)はPCを変えただけで変化する不特定な番号なので、あくまで今回の計算に限定された指定番号となります。
※１：ωは面直方向変位(Z軸)、ωx、ωyはx、y方向の傾き(ωx=∂ω/∂x、ωy=∂ω/∂y)
３．１．２ 動的解析モデル
前項を受けて動的解析用に以下のscilabモデルを作成しました。
（１）薄板構造部
3.1項で設定した節点モデルを元に以下の行列を作成します(詳細は「Detail1」変位法による薄板FEM解析詳細及び「Detail2」変位法による薄板FEM動的解析詳細をご参照下さい)。
・剛性行列[Kf]
・質量行列[Mf]の逆行列[MIf]
・減衰係数1％(構造減衰率2％)とした減衰項2*.01*[KMsqrt](詳細は「Detail2」変位法による薄板FEM動的解析詳細をご参照下さい) この減衰係数は薄板内部構造に依存しており、試験値及び経験値との比較が重要となります。
（２）定常外力
初期状態の定常外力の条件を以下としました。
・Max300xMax800x6mm板(メッシュ粗さから2mmを6mmに変更)
・迎角2°
・2[m/s]
この条件のもと、OpenFoamを用いて薄板の上下面に働く圧力を計算しました。
OpenFoam設定
・空間2.4(長さ) x 1(高さ) x1.2(幅)[m]を48x20x24の合計23040マスに分割(1マス50[mm]立方)
・薄板近辺のメッシュ粗さは上記１マスを1/16に緻密化(1マス3mm立方)
・入口速度2[m/s]、出口は内部計算値、薄板表面及び薄板を取付た壁は表面速度0、それ以外の面はslip条件
・出口圧力0、その他はzeroGradien条件
・空気動粘性：2.0e-5 [m2/s]
・解析ソルバーはpisoFoam(非圧縮・非定常乱流解析ソルバー)
・乱流モデルはRAS(レイノルズ平均に基づく乱流モデル)でkEpsilonモデル(標準ｋ-εモデル)
・epsilon、kは初期値0.1、nutの初期値は0（適当）。
上記以外のパラメータは以下のOpenFoamチュートリアルの物をそのまま使用しました。
tutorials\incompressible\pisoFoam\RAS\cavity
計算結果をparaviewのsave dataコマンドでCSVに吐き出し、71点に対応した点の表面と裏面の圧力差分を行列[Ini_force]とします。この計算では周囲速度を2[m/s]としました。一方今回は回転数RPMをパラメータとしますので速度U=回転数RPM×回転半径Rとして各節点での圧力差行列を補正した[Ini_forceRPM0]を生成します(基準回転数をRPM0=20[rpm]とします)。
例えば節点56ではX=.675m。これにスピナー半径.15[ｍ]を加えるとU56=20/60*2π*(.675+.15)=1.73[m/s]。そこで節点56の圧力差は(1.73/2)2の補正を加えました(回転運動では半径毎に速度が変化して半径別の流管間でせん断力が発生するため、平行気流とは流れが違いますが、ここでは無視します)。
さらに回転数RPMの時は行列[ini_forceRPM]=[Ini_forceRPM0]x(RPM/RPM0)2の補正を行いました。
尚薄板の縁は圧力変化が急なので6[mm]内側にオフセットした圧力を読み取ることとしました。
またScilab計算上、急激な入力による計算発散を防ぐため200[msec]のランプ入力としています(モデル上ランプ・ブロック)。
[Ini_forceRPM0]の様子(RPM=20[rpm])を以下に示します(縦横比適当)。揚力のほとんどを前縁が担っています。
ここで定義した定常外力に対する変形量を上述Scilabモデルで計算します(減衰係数10％で5 秒間)。
結果は以下で、端部で0.53[mm]浮き上がっています(下図は変位を拡大して表示しています)。
この計算結果の行列を[Ini_deflectionRPM0]としておきます。こちらにも速度補正を加えておきます。
[Ini_deflectionRPM]=[Ini_deflectionRPM0]x(RPM/RPM0)2
この変形量は次の項で使用します。
（３）１次捩り空気力
(2)項で定義した定常外力に対し、追加で微小な１次捩りMode変形が起きた際の空気力増減をOpenFoamで計算します。
・Scilabで計算したMode201の薄板形状を6[mm]厚にSolid化。
(今回は捩り形状が大きかったので捩りを半分とし、計算結果を2倍としました)
・迎角0°
・2[m/s]
・OpenFoam設定は(2)項と同じ
計算結果をparaviewのsave_dataコマンドでCSVに吐き出し、先ほどと同様に表面と裏面の圧力差分を求め、行列[Air_Mode201]とします。この行列は迎角0°の薄板がMode201(1次捩り)の変形をした場合の空気力の増加分を意味します。前項と同様な回転数補正を実施した後の[Air_Mode201]の様子を以下に示します(縦横比適当)。外舷側前縁が最も空気力を発生している様子が分かります。
この空気力増加分[Air_Mode201]をScilabモデルに取り込みます。取り込み方法はScilabを用いた薄板フラッター解析と同じためここでは省略します。
（４）1次曲げダンピング力
(2)項で定義した定常外力に対し、１次曲げモードの速度変形が起きた際の空気力増減を下記方法で求めます。
・[Ini_force]を補正することで作成。節点68個の点毎に実施
例として節点63の計算を示します。
・Mode200の形状を単位時間(1秒間)に変形した量と仮定。節点63では-73.6[mm]
・基準回転数RPM0(20[rpm])における局所速度は20/60*2π*(.15+.75)=1.88[m/s]
と上記薄板変形速度-73.6[mm/s]により局所的に迎角が変化。
・局所速度と変形速度による迎角変化量tan-1(-73.6/1880)=-2.2°
・[Ini_force]は迎角2°の揚力なので、[Ini_force]の節点63の値を2°で割り、-2.2倍したものを求めるべきダンピング力[Air120_AOA]とする。
尚薄板は気流に対して向かっていく方向に変形しているのでダンピング力[Air120_AOA]の極性が「+」になるようにScilabモデルのSUMブロックで「-」を使います。
取り込み方法はScilabを用いた薄板フラッター解析と同じためここでは省略します。
上記計算を各列毎に実施したダンピング行列[Air120_AOA]が以下となります(縦横比適当)。
（５）1次遅れ
本来薄板のフラッタ現象では構造変形に対し空気力変化が一定の遅れを持つことが重要と言われています。この点を時定数taoの１次遅れブロックによりの簡易に模擬しました。
方法はScilabを用いた薄板フラッター解析と同じとなりますが、回転しているため以下を変更する必要があります。
回転がない場合の時定数taoは以下となります。
tao=薄板幅(0.3[m])÷周囲速度(2[m/s]) x 0.5 (掛け数：適当)= 75[msec]
回転がある場合は周囲速度の代表を節点63とし、以下の式で時定数taoを規定しました。
tao=薄板幅(0.3[m])÷局所速度(RPM/60*2π*(.15+.75)[m/s])x 0.5 (掛け数：適当)
回転数RPM=20[rpm]の場合、tao=80[msec]となります。
これはかなり簡易な模擬なので実際の使用では掛け数0.5を0～1位で変化させ、解析結果がどれくらい振れるかを確認すべきと考えます。
以下に回転数300、380、400[rpm]の計算結果を示します。
中回転(380[rpm])では振動発散していないのですが、僅かに回転数の上がった400[rpm]で急激な振動発散を起こしています。この発散の際、Mode200(1次曲げモード 1.7[Hz])とMode201(1次捩りモード 10.1[Hz])の周波数は共に中間値5.9[Hz]近傍の6.6[Hz]に変移しており、古典的フラッタの様子を示していると考えます。
３．２ 遠心力のある場合
３．２．１ Scilab基本モデル
基本的なScilabモデルは3.1項と同じ以下となります。
・300x800x2mmアルミ板
・ヤング率70GPa / 密度2800kg/m3 / ポアソン比0.3
・71節点(1要素は三角形状)で1点辺りの自由度3(ω、ωx、ωy)
・端部固定(①②③)により自由度は3点減って、全自由度は68×3=204。よってMode数も204。生成されるModeは3.1項と同じとなります。
３．２．２ 動的解析モデル
動的モデルは以下となります。3.1項のモデルに遠心力を表現する行列[Nstr]を追加しています。
（１）遠心力部
Scilabを用いた膜の解析にて膜に張力を与えた場合のモデルを作成しました。ここではその張力を遠心力に差し替えることを考えます。
遠心力の様子を下図に示します。
①-④間要素を注目します。
①に働く遠心力CF1のωx1分と④に働く遠心力CF4のωx4分がN14として要素①-④に働きます。
N14が半分ずつ節点①と節点④に働くとすると式は以下となります。
F1＝N14/2＝(-CF1*ωx1+CF4*ωx4)/2： (遠心力CFi＞０として符号を合わせます)
F4＝N14/2＝(-CF1*ωx1+CF4*ωx4)/2
同様にしてF2、F3、F4・・・と式を作り、行列としてまとめると
この行列にωとωyの行、列を「0」として入れ込むことで[Nstr]が完成します。
尚CF1、CF2、CF3、・・・については以下の図のように区切りをつけて計算しました。 例えば区分9の質量が重心に集中したとして遠心力mrω2を計算します。同じように区分別の遠心力を計算します。節点⑧に印加するCF8は3～12区分の遠心力を合計し、それを3点分ということで3で割ることで計算します(尚[Nstr]では更にこの遠心力を２で割ります)。
以下に回転数300、380、600[rpm]の計算結果を示します。
低回転(300[rpm])及び高回転(600[rpm])でフラッタが起きていません(遠心力なしでは400[rpm]でフラッタに移行しています)。
但し中回転(380[rpm])では振動発散しています(遠心力なしでは380[rpm]ではフラッタに移行していません)。
ここでもう少し中回転数と高回転数の計算結果を確認してみます。
まずは高回転数(600[rpm])の様子を再確認します。
振動が残っているMode200の信号にFFTをかけてみました。結果、推定理論値※1である9.1[Hz]と27.2[Hz]に十分近い11.5[Hz]と31.3[Hz]の周波数ピークを確認できました。すなわち回転数が600[rpm]まで来ると弦の理論式と一致しだすようです。
更にピーク周波数では薄板１次曲げモードと１次曲げモードが混ざりあっている様子が確認できました。 今まで調べてきた古典的フラッタは１次捩りモードの周波数が落ちていき、１次曲げモードと連成していました。しかしここでは既に１次曲げモードと１次捩りモードが混ざりあい、且つ周波数の変移が起きていないため振動発散にならないようです。
※１：膜の固有振動数理論値(片端固定条件：1次元の弦とした)
fi=(2i-1)/(4*ℓ)×√(P/(ρ)) =9.1[Hz]、27.2[Hz]
P : 張力の元となる面圧 2.36e6[Pa]
ρ：AL密度 2800[kg/m3]
ℓ：0.8[m]
i=1、2・・・
Ｐ=CF/S/2 遠心力は傾斜配分なので1/2とした
CF：遠心力。代表点として節点㊷を選択
CF=0.19m2(面積)×0.002m(厚み)×2800(密度)×0.675[m](半径)×6002[rpm2]=2835N
S：断面積。ブレード本体の幅0.3[m]×0.002[m]=0.0006m2
続いて中回転数(380[rpm])の様子を再確認します。
遠心力無しの場合、380[rpm]では振動発散していないのに、遠心力を考慮すると振動発散しています。
詳しく見てみるとMode200及びMode201共に8.8[Hz]で、古典的フラッタにみられる1次曲げ1.7[Hz]と1次捩り10.1[Hz]の平均である5.9[Hz]に近づいていません。つまりこの振動発散は古典的フラッタの動きではなく、別の原因で振動発振してるようです。
試しに1次曲げによる空気力をなくし、且つ1次捩れモードの1次遅れを無くしてみました。またMode200の抽出点を2階積分後に移動してみました。
結果が以下で、1次曲げ空気力[Air200_AOA]がなくても振動発散しています。図中からわかることとして1次曲げMode200と遠心力[Nstr]の波形位相が約30°ずれており、これが振動発散の原因と考えられます。遠心力[Nstr]と[Point69]、[Mode200]は位相が一致しており、これだけならば発振しないのですが、そこに1次捩り空気力[Air_Mode201]が位相をずらして加振してくるため振動発散していると考えられます。1次捩り空気力[Air_Mode201]の中には複数のモードが含まれているのですが、ちょうど1次曲げMode200は成分が「-」で発散側となっていました。
そこで[Air_Mode201]成分からMode200の成分を抜き取る※2と振動が収束していくことを確認しました。
振動発散を考える時、Mode間の位相も注意深く観察する必要があることを認識しました。
※2：1次捩り空気力[Air_Mode201]からMode200の成分を抽出
・基底成分を抽出：[A]＝[FAII]×[Air_Mode201] (FAII：固有ベクトル逆行列)
・Mode200の基底成分を「0」に変更
・基底を一般座標に戻す：[B]=[FAI]×A (FAI：固有ベクトル行列)
・この時、行列Bのωxとωyに変換によるカスがたまるので手動で「0」に修正
・最後に[Air_Mode201]=[B]
「おわりに」
ここまで示した当事業所オリジナルのブレード・フラッタ解析は試験値との比較を行っていないため、あくまで参考計算となります。
その上でscilabモデルから分かったフラッタの特徴を下記に示しておきます。
1. １次捩りモードのみのモデル(１次曲げモードなし)
１次捩りモードの固有振動数は低回転では薄板構造剛性[K]に従ったモード振動数(本稿では10.1[Hz])を示します。
一方回転を上げていくと1次捩り空気力[Air_Mode201]は回転の二乗に比例して大きくなり、薄板構造剛性[K]を打ち消していきます。このため１次ねじりモードの固有振動数はひたすら小さくなっていきます。
そしてある速度で薄板構造剛性[K]を完全に打ち消す回転数に到達します。この段階でブレードは振動を伴わない発散(ダイバージェンス)を起こすことを確認しました。
つまりフラッタを起こさないブレードの限界はこのダイバージェンスと考えられます。
2. 1次曲げモードと1次捩りモードのモデルの空気力の時間遅れあり/なし
1次曲げ運動と1次捩り運動は空気力[Air200_AOA][Air_Mode201]を通してカップリングしています。このためかなり低い回転からでも時間をかけて振動発散に移行していきます。
結果、このモデルでは極低回転からフラッタに入ってしまうモデルになっています。
一方空気力の時間遅れモデルは1次フィルタとして働き、特定の回転数までは1次モードの振動を吸収していると考えられます。
しかし特定の回転数を超えると一気に振動発散に移行します。これがフラッタを表現しています。この非線形性は1次フィルタのフラット部と斜辺部により生成されていると考えます。
3. 遠心力の効果
遠心力[Nstr]は空気力[Air_Mode201]と反対の極性を持っています。このため遠心力が優っている場合、フラッタ・モデルの特徴は消失していきます。
遠心力が十分に大きくなると、膜の１次曲げ、1次捩りモードがそれぞれ薄板1次曲げ、1次捩りモードのまじりあった状態になり、且つ1次捩りモードの周波数が低下する現象もなくなるためここで取り上げてきた古典的フラッタはもう起きない状況となりました。
但し中回転では1次捩りモードと遠心力モードの過渡期ではMode間の位相ずれから振動発散が見られました。
HA-technodefy 