科学館員の独り言

また！平板ブーメラン（１５）

未来人 — 2008-07-26T16:14:20+09:00

　エクセルのマクロで今回のシミュレーション計算が簡単にできます。添付ファイルはエクセルのワークシート*と、ビジュアル・ベーシックのソース･ファイル**に分けてあります。ソース･ファイルはエディターで閲覧・編集ができます。

　使うときは、まず、［ファイル］→［名前の変更］で、ソース･ファイルの拡張子をtxtからbasに変更します。次に、エクセルを開いて、［ツール］→［マクロ］→［Visual Basic Editor］としてビジュアル・ベーシック・エディターを開き、［ファイル］→［ファイルのインポート］でソース･ファイルを選択します。その後は、ワークシートの初期条件を入力して、［ツール］→［マクロ］→［マクロ］→［実行］で計算します。

　初速を大きくしたり、迎え角の初期値を大きくしたり、逆にマイナスにすると、過大な初期条件となって発散します。また、時間ステップも大きくすると発散します。

　入力が必要なのは、迎え角の増大係数CA、揚力係数CL（3次式まで可能）、抗力係数CD（2次式まで可能）、時間ステップT、ブーメランの重量W、実翼面積S、重力加速度G、揚力倍率Fl、抗力倍率Fd、翼面積倍率でFs、初期値としては、迎え角αと面の法線角φのほか、表中の黒枠部分に、位置（x,y,z）と初速度（vx,vy,vz）を入力します。揚力倍率、抗力倍率、翼面積倍率は、通常は1のままで構いませんが、係数を変えて結果がおかしいときなどに、倍率を変えて様子を見ます。

　ワークシートには、平板ブーメラン用のflatと、回転揚力のあるブーメラン用のflapの2種類があります。flapの場合は、回転揚力Lfを入力する必要があります。

　計算結果は、時間ステップごとに表として出力されるほか、x-y断面とx-z断面、y-z断面、xn-yn断面とxn-zn断面で、軌跡が作図されます。表の出力は、時間t、位置成分（x,y,z）、速度成分（vx,vy,vz）、速度v、法線ベクトル成分（xn,yn,zn）、法線角φ、迎え角α、揚力係数CL、揚力L、抗力係数CD、抗力D、揚力成分（Lx,Ly,Lz）、抗力成分（Dx,Dy,Dz）、上限速度v0、速度比ce、座標変換角（θ、θ'、θ''）、位置増分（⊿x,⊿y,⊿z）、進行ベクトル原成分（x3,y3,z3）、法線ベクトル原成分（xn3,yn3,zn3）となっています。

　ソース･ファイルを書き換えれば、もっといろいろなことができるようになります。あるいは、自分の考えで、改造することも可能です。これをヒントにして、やりたいようにやってみてください。完成したプログラムではないので、いろいろと問題がでてきたり、バグがあったりするかもしれませんが、ご容赦ください。

添付ファイル：
*simulation3s6.xls「simulation3s6.xls」をダウンロード
**boo_sim3s6.txt「boo_sim3s6.txt」をダウンロード

また！平板ブーメラン（１４）

未来人 — 2008-07-25T20:13:29+09:00

　2005年12月から約2年半をかけて、一応、満足できる平板ブーメランの飛行シミュレーションが可能になりました。
　この、こだわりのきっかけは、ブーメランの飛行原理の通説に疑問を持ったからです。あちこちの科学館でも、「ブーメランはジャイロ効果で旋回する」と解説しています。ホームページでも、いろいろと解説されていましたが、難解との指摘があることも分かりました。一方、平板ブーメランは飛行原理そのものが不明だったので、新たに原理を構築する必要があります。平板ブーメランは、シンプルなので理論化が容易になるとも考えました。

　平板ブーメランの飛行原理を検討しているうちに、「平板翼では、気流の中では、面が気流に正対する性質がある」ことに気が付き、実験的に定式化して、飛行のシミュレーションができるのではないかと考えたのですが、シミュレーションの開発過程で、ジャイロ効果を表す式が、シミュレーションにはまったく役に立たないことが分かりました。飛行そのものは重心移動なので、ジャイロ効果が示す回転運動だけで旋回飛行の運動方程式は解けません。むしろ、ジャイロ効果を考慮しなくても、「面が気流に正対する」ことを定式化することで、シミュレーションが可能になりました。

　このメカニズムが、平板ブーメランだけでなく、一般のブーメランに適用できることは自明です。ブーメランが旋回飛行をするのは、（自転によって等方化された）自転面に作用する揚力が作るトルクが支配的で、ジャイロ効果の影響は小さい、というのが今回の結論です。ブーメランが旋回するには、自転が不可欠なのはいうまでもありませんが、自転の効用は、主としてコマと同様な自転軸の維持というジャイロ効果にあります。つまり、ジャイロ効果は、ブーメランの飛行安定のための復元力として作用しているといえます。

また！平板ブーメラン（１３）

未来人 — 2008-07-24T15:18:13+09:00

　さて、このシリーズの（１）回で述べた、面白い現象とは、次のようなものです。ブーメランの自転落下実験を側面からではなく、正面から観測すると、ブーメランが水平方向にわずかに移動します。自転落下実験では、回転方向はすべて、正面から見て反時計回りに回転させています。

　自転によって進行方向が曲がる現象では、「マグヌス効果」が知られています。例として挙げられる、野球のボールやゴルフボールの軌道の曲がりなどは、「マグヌス効果」では、揚力の方向は球の回転に一致する方向であるとして説明されています。つまり、上から見てボールが時計回りに回転すればシュートに、反時計回りに回転すれば（昔風の）カーブになります。ドライバーで打ったゴルフボールがホップするのは、縦に回転し、右側面から見た場合に反時計回りに回転するからです。

　しかし、観測されたブーメランの曲がり方は、「マグヌス効果」とは違って自転とは逆方向になっています。この解釈については、プリミティブには、落下の際に進行方向となる側での翼の「漕ぎ方」の方が、後方の逆向きの「漕ぎ方」よりも効果が大きくなるため、といった「後流効果」でも説明できそうですが、次のような解釈の方が、もっともらしいと考えています。
　翼は回転しているので、進行方向の速度に加えて、右の翼では自転による回転が速度として加わり、左の翼では、逆に、速度が減じます。迎え角は同じで、翼の断面形状も同じなので、翼の抗力係数が同じと仮定すれば、抗力は速度の二乗に比例することから、反時計回りの回転の場合は、正面から向かって左側の翼に働く抗力の方が右側よりも大きくなります。そのため、左右の翼に抗力の差が生じ、結果として向かって左方向に向かう力が働く、という解釈です。
　今回のシミュレーション改良版では採用していませんが、この効果を取り入れると、飛行初期の落下が抑えられるので、より現実に近づく可能性もあります。より現実に近いシミュレーションには、やはり、本格的な風洞実験が必要です。今回のような、簡易な自転落下実験には限界があります。

また！平板ブーメラン（１２）

未来人 — 2008-07-23T17:08:20+09:00

　さて、このように、手頃なツールができたので、いろいろと遊ぶことができます。そこで、前回と同様に、宇宙ステーションでブーメランを飛ばしたらどうなるか、を計算してみました。そのために、普通のブーメランと平板ブーメランとの違いは、面の法線方向に自転による揚力があるいかないかである、と考えてみます。普通の紙製ブーメランを模擬した「フラップ付き」の自転落下実験もやってみましたが、迎え角の増大も、揚力係数も、抗力係数も、平板よりは大きめに出てきます。傾向はほとんど同じなので、実験式を入れて計算してみると、なかなかうまく行きません。発散することが多いようです。そこで、実験式は平板ブーメランと同じにして、普通のブーメランの場合は、面の法線方向に一定の揚力を加えてみました。

　「フラップ付き」の計算例として、（１１）回の旋回飛行のケースに準じて、高さ2mから、初速度をx軸の正方向に8m/s、迎え角の初期値を0度、面の法線の傾きを70度とし、面の法線方向に0.08Nの揚力を加えてみました。この0.08Nという値は自転落下実験のデータから見ても、無理な仮定ではありません。

　これは、（１１）回で示したフラップなしの結果と比較して、一周して戻る傾向が強く出るようになりました。

　次に、面の法線方向に0.08Nの揚力を加えたうえに、面の法線の傾きを80度とし、宇宙ステーション内を模擬して重力加速度を0にしてみました。

　きれいに一周する結果となりました。前回のらせん状に上昇するシミュレーションとは異なる結果となりましたが、揚力の値を小さくすると、らせん状に上昇する結果にもなります。投げ方で変るということでしょうか。

　いろいろと条件を変えて計算することができますが、条件によっては発散したりします。また、時間ステップを変えると結果が微妙に変りますが、原因は不明です。
前回のシミュレーションでは、平板ブーメランでは重力による落下で迎え角を生ずるとしていましたが、今回は、迎え角の初期値を入力できるようにしたので、逆に迎え角の初期値を0とすると、なかなか旋回しない計算結果もあります。シミュレーションによれば、平板ブーメランでも迎え角を付して投げると、旋回性が良くなるようです。実際には、気付かずに迎え角を付して投げているのかもしれません。

また！平板ブーメラン（１１）

未来人 — 2008-07-22T16:13:07+09:00

　さて、平板ブーメランの飛行シミュレーション改良版の検証です。まず、典型的な例として、高さ2mから、水平に、初速度のx成分を8m/s、y成分を-1m/s、z成分を0.5m/sとして投げ、面の法線の傾きが70度、迎え角の初期値が4度の場合の、x-z断面とx-y断面の計算結果を示します。

　旋回直径が実態より小さめですが、飛行軌跡がほぼ再現できています。同じ条件での、y-z断面の計算結果では、

中間部でやや降下が早いようですが、まずまずの結果と思います。

　次に、今回の目標の一つである、「ブーメランを高いところから傾けて落下すると、スロープを描いて軸が立ち上がるような軌道」の例としては、高さ3mから、面の傾きを65度（初期迎え角としては、25度に相当）として、x軸の正方向に初速度0でリリースした場合の計算結果を示します。

　これも、落下地点までの距離が短いようですが、傾向は再現されています。

もう一つの目標である、「ブーメランを水平に投げると、上昇した後、突然、反転して戻ってくるような軌道」については、かなり微妙な結果になりました。現実もそうだ、ということの反映ではないかと思いますが、具体的には、ブーメランを旋回させるときよりも初速度を大きくする必要があり、わずかな条件の違いで、山を登ってから向こうに下りていくようなケースと、反転して戻るケースと、宙返りをするケースが出てきました。

　反転して戻る例では、初速度14m/sで、面の傾きをほぼ0度にして、迎え角の初期値を5.5度とした場合の計算結果を示します。

　宙返りをする例としては、初速度12m/sで、やはり面の傾きをほぼ0度にして、迎え角の初期値を4.5度とした場合の計算結果を示します。

　このほか、初速度が小さい場合は、山を登りきってから向こうに下るケースになります。

　同じ実験式のセットを適用して、目標とした3パターンのシミュレーションが良好に検証できました。まだ、満点ではありませんが、一応、合格点と考えています。

　ここで重要なことは、これらのシミュレーションでは、ブーメランの自転に、揚力などの外力がトルクとして加わった場合の、際差運動、いわゆるジャイロ効果をアプリオリには考慮していないということです。つまり、際差運動を考慮しなくても、ブーメランの旋回がシミュレーションできることから、ブーメランの旋回運動にとってジャイロ効果は主要な因子ではない、ということになります。

また！平板ブーメラン（１０）

未来人 — 2008-07-21T14:49:12+09:00

　エクセルのマクロ（ＶＢＡ）でプログラムを組むときに注意しなければいけないのは、角度の作り方です。例えば、xyz座標系からx'y'z'座標系に変換するための角度θは、tanθ=⊿y/⊿xの関係があるので、プログラム上は⊿y/⊿xのアークタンジェント（逆正接）として算出します。ところが、ＶＢＡのアークタンジェントで算出した角度は、-π/2からπ/2の間で定義されているので、今回のように、-πからπまで変化させるには、⊿yと⊿xの組み合わせを考えて、物理的に意味のある変換をしないと、正弦や余弦の値（とくに符号）がおかしくなります。突然に、プラスからマイナスに飛ぶような変化や、その逆の変化が現れた場合は、その可能性を疑う必要があります。

　今回は、角度θに関しては、次のような変換をしています。
If x = 0 Then th = 0 Else th = Atn(y / x)
If x < 0 And y > 0 Then th = th + Pi
If x < 0 And y < 0 Then th = th - Pi
ここで、thはθ、xは⊿x、yは⊿y、に相当します。なお、Piは、Pi = 4 * Atn(1) として定義した定数です。

　同様に、他の角度についても、それぞれの考え方で変換を行っています。x'y'z'座標系からx''y''z''座標系に変換するための角度θ'については、
If x1 = 0 Then th1 = Pi / 2 Else th1 = Atn((-z1) / x1)
If x1 < 0 And -z1 > 0 Then th1 = th1 + Pi
If x1 < 0 And -z1 < 0 Then th1 = th1 - Pi
としました。ここで、th1はθ'、x1は⊿x'、z1は⊿z'、に相当します。

　また、x''y''z''座標系からx'''y'''z'''座標系に変換するための角度θ''についても、
If zn2 = 0 Then th2 = 0 Else th2 = -Atn(yn2 / zn2)
If zn2 < 0 And yn2 < 0 Then th2 = th2 + Pi
If zn2 < 0 And yn2 > 0 Then th2 = th2 - Pi
となっています。ここで、th2はθ''、yn2は⊿yn''、zn2は⊿zn''、に相当します。角度θ''だけは、通常、負なので、他とは扱いが異なります。

　また、本来は、抗力は、推進力と揚力の結果として現れる力のはずですが、シミュレーションでは抗力係数と速度から自動的に算出されるので、過大になる可能性があります。そのため、ステップごとの速度には、前のステップの運動エネルギーと、落下によって加わった位置エネルギーの和を超えないというエネルギー上の制限を設けています。

また！平板ブーメラン（９）

未来人 — 2008-07-20T18:42:35+09:00

　ブーメランに作用する力は、揚力L、抗力D、重力-Wgだけで、与えられた初速によって推進します。運動方程式から導かれる差分は、

図の定義に基づいてxyz座標系に変換した後、

などのようになります。これで、迎え角αに対する揚力L、抗力Dから、位置の差分⊿x、⊿y、⊿zと、速度の差分が求められます。

　実験式は、このシリーズの（６）、（７）で述べたように、自転落下実験のデータから、エクセルの最小二乗近似式で求めています。今回、採用した実験式は以下のとおりです。簡単な直線近似から始めてたどりついた、現状でベストと思われるセットです。

　揚力Lと抗力Dの和Fと時間ステップΔtから、迎え角の増分Δαを算出し、次のステップの迎え角αを決めます。迎え角αが決まれば、実験式から揚力係数と抗力係数が決まるので、既に求められている速度vを前出の式に代入して、次のステップの揚力Lと抗力Dが求められます。その後は、これらの操作を繰り返します。

また！平板ブーメラン（８）

未来人 — 2008-07-19T17:47:04+09:00

　次は、計算シミュレーションの方法です。前回と同様に、エクセルのワークシートとマクロ（ＶＢＡ）を組み合わせました。運動方程式の解法については、全面的に改めて、素直に座標変換を行いました。図のように、1軸の回りの座標変換を3回繰り返します。

　まず、xyz座標系をブーメランの進行方向に合わせ、z軸の回りに回転させて、x'y'z'座標系に変換します。次は、このx'y'z'座標系をブーメランの進行方向にx'軸が一致するようにy'軸の回りに回転させて、x''y''z''座標系に変換します。最後は、x''y''z''座標系を、ブーメランの面の法線方向を含むようにy''-z''面に一致させて、x'''y'''z'''座標系に変換します。その結果、x'''y'''z'''座標系では、ブーメランを側面から見た形となるので、その面とx'''軸がなす角度が迎え角となり、z'''軸の正方向に揚力が、x'''軸の負方向に抗力が作用することになります。また、面の法線ベクトルを単位ベクトルとすれば、法線ベクトルはx'''-z'''面内で、z'''軸から迎え角だけ傾き、xn'''は-sinαと、zn'''はcosαと表現できます。

　ブーメランの進行ベクトル、面の単位法線ベクトルの座標を、xyz座標系から、x'y'z'座標系、x''y''z''座標系、x'''y'''z'''座標系と変換し、αの増分によって揚力L'''と抗力D'''を決定し、点N'''を決定した後に、逆の手順でxyz座標系に戻して表現します。
　一例として、xyz座標系をz軸の回りに角度θだけ回転した場合の座標変換は、
　　x'=xcosθ+ysinθ
　　y'=-xsinθ+ycosθ
　　z'=z
で変換でき、逆には、
　　x=x'cosθ-y'sinθ
　　y=x'sinθ+y'cosθ
　　z=z'
で変換できます。このルールを、固定軸に一致させて、回転軸の対応を都度決めれば、同じ形式で全部の変換が可能です。ただし、座標系は右手系で考えます。また、式の定義では角度θの正負は問いませんが、座標から角度を決める際には、角度の正負を物理的に考慮する必要があります。

　なお、このような座標変換では、基本的にブーメランの進行方向をx'''軸に一致させています。したがって、ブーメランの進行方向が上下、つまりz軸方向を向く場合は、そのままでは無理があります。今回は⊿z=0となる場合は、便宜的に、ケースを仮定して角度を読み替える解決をしています。そのケースとは、今回の目標としている「ブーメランを高いところから傾けて落下すると、スロープを描いて軸が立ち上がるような軌道」と、「ブーメランを水平に投げると、上昇した後、突然、反転して戻ってくるような軌道」の場合です。

また！平板ブーメラン（７）

未来人 — 2008-07-18T20:35:58+09:00

　揚力係数は、検討の結果、単純な2次式で近似することにしました。

　上の近似式はデータの最小二乗式ですが、αが0とπ/2で揚力係数が0となることを仮定して、最大値を合わせた2次式に換算してあります。この式ではαがπ/4で最大値となりますが、揚力係数の一般的な傾向からは、実際にはαが0.5以下で最大になると考えられます。その意味では、この近似式を高角で使用するには、注意が必要です。なお、αが負の場合は負号で表現し、αがπ/2以上の場合はπ-αをαと読み替えて近似式を適用します。

　また、抗力係数も2次式で近似することにしました。

　αが負の場合は、-αをαと読み替え、αがπ/2以上の場合は、π-αをαと読み替えて使用します。

　なお、これらは、回転速度が毎秒約2回転のデータによるものです。実際には揚力係数も抗力係数も回転速度依存があるので、ずれが生じてきます。元来、ばらつきの大きいデータなので、その中からある程度の取捨選択をした結果ですが、結果的には妥当な線と思われます。

また！平板ブーメラン（６）

未来人 — 2008-07-17T16:36:29+09:00

　自転落下実験は、約1mの高さのほかに、約2mの高さでも実施しました。高さ1mでは最速3m/s、高さ2mでは最速7m/s程度になります。また、側面以外に、正面から撮影したケースもあります。側面の場合の傾斜角度は、0度から70度までの10度刻みと、90度で実施しています。回転速度は、毎秒約2回転と4回転を目標にしましたが、電池駆動のため、回転数の測定中に回転数が低下するなどの原因で、多少、ずれがありました。
　この2月と8,9月にかけて、トータルで16シリーズ264回の実験を行っていますが、途中の装置改良やデータ処理の見直しもあって、最終的に採用したデータは、4シリーズ36回分のデータに止まりました。しかし、1回の落下で動画のコマ数は10ないし20なので、それでも結構な数になります。その結論は以下のとおりです。

　まず、シミュレーションの際に、揚力係数や抗力係数を算出するには、迎え角を決める必要があります。迎え角の増大は、揚力が翼の揚力中心に作用し、一般に揚力中心は重心よりも翼断面の前方にあるために、トルクを生じて回転する結果、生ずる現象です。厳密には回転の運動方程式を解く必要がありますが、回転の運動方程式は回転角の2階微分が外力によるモーメントに等しいとするので、自転落下実験において、迎え角の増分⊿αを微小時間⊿tの二乗で除した値が、ブーメランに作用する揚力と抗力の和に比例すると考えました。揚力と抗力は直交するので、それらの和は揚力と抗力がつくる矩形の対角線に相当します。その結果が次の図です。

図のFが揚力Lと抗力Dの和です。