NAZCA-TR-23-008A

ゴム動力飛行機の３次元運動解析

滝敏美

2023 年 5 月 18 日
A 改訂：2023 年 6 月 16 日i

ii

ゴム動力飛行機の３次元運動解析
目次
1

研究の目的 .............................................................................................................................................................. 11.1

ゴム動力機の概要 .................................................................................................................................................. 1

1.2

解析で考慮する項目 .............................................................................................................................................. 22

模型飛行機の３次元運動方程式........................................................................................................................... 32.1

機体の姿勢と速度 .................................................................................................................................................. 32.2

運動方程式 .............................................................................................................................................................. 4

2.3

機体に作用する空気力 .......................................................................................................................................... 62.4

空力微係数の計算方法 .......................................................................................................................................... 83

動力ゴム .................................................................................................................................................................. 9

3.1

使用するゴム .......................................................................................................................................................... 93.2

巻き数とトルクの関係 .......................................................................................................................................... 93.3

巻き数とトルクのモデル .................................................................................................................................... 114プロペラ ................................................................................................................................................................ 134.1

プロペラの形状と寸法の定義 ............................................................................................................................ 134.2

機体固定座標系とプロペラ座標系の関係 ......................................................................................................... 144.3プロペラの空気力 ................................................................................................................................................ 154.3.1

プロペラの性能曲線の計算法 ........................................................................................................................ 15

4.3.1.1プロペラの性能曲線 .............................................................................................................................. 15

4.3.1.2

Larrabee の方法 ....................................................................................................................................... 16

4.3.2

気流がプロペラの面に正対していない場合に作用する力 ........................................................................ 22

4.3.2.1

Ribner の理論 .......................................................................................................................................... 22

4.3.2.2

McCormick の理論 .................................................................................................................................. 25

4.3.2.3

機体固定座標系から見たプロペラの横力，モーメントの方向 ....................................................... 26

4.3.3
4.4

プロペラの空転................................................................................................................................................ 27
プロペラの回転慣性の影響 ................................................................................................................................ 284.4.1

プロペラのトルクの釣り合い ........................................................................................................................ 284.4.2

プロペラのジャイロ効果 ................................................................................................................................ 31

5

ゴム動力機の３次元運動方程式とその解法 ..................................................................................................... 32

5.1

運動方程式のまとめ ............................................................................................................................................ 32

5.2

運動方程式の解法 ................................................................................................................................................ 32

6

ゴム動力機の３次元運動解析の例..................................................................................................................... 33

6.1

解析対象 ................................................................................................................................................................ 33

6.2

機体の空力特性 .................................................................................................................................................... 39

6.3

慣性データ ............................................................................................................................................................ 446.4

動力ゴムの特性 .................................................................................................................................................... 45

6.5

プロペラ特性 ........................................................................................................................................................ 48

6.5.1

プロペラ仕様 ................................................................................................................................................... 48

6.5.2

プロペラ性能曲線............................................................................................................................................ 51

6.5.3

プロペラに見流が斜めからあたる影響 ........................................................................................................ 53

6.5.4
6.6

プロペラ慣性能率............................................................................................................................................ 54
３次元運動解析結果 ............................................................................................................................................ 55

iii

6.6.1

機体固定座標系................................................................................................................................................ 55

6.6.2

解析ツール ....................................................................................................................................................... 56

6.6.3

解析ケースと解析結果 .................................................................................................................................... 56

7

結論 ........................................................................................................................................................................ 818

参考文献 ................................................................................................................................................................ 82

付録１ ゴム動力機の３次元運動方程式を解く Octave スクリプト .......................................................................... 83
付録２ 機体の図を描くスクリプトの解説記事 ........................................................................................................... 96iv

図のリスト
図 1.1-1

ゴム動力機の実例........................................................................................................................................... 1

図 2.1-1

機体の姿勢と速度........................................................................................................................................... 3

図 2.1-2

機体に作用する力........................................................................................................................................... 3

図 2.1-3

地面固定座標系と機体固定座標系の関係，オイラー角の定義................................................................ 4

図 2.3-1

垂直尾翼の取付角........................................................................................................................................... 7

図 2.3-2

水平尾翼の傾きとエルロン舵角の符号 ....................................................................................................... 7

図 2.3-3

主翼単独の無揚力線と水平尾翼取付角の関係 ........................................................................................... 7

図 3.1-1

結び目を作らない輪ゴムのつなぎ方 ........................................................................................................... 9

図 3.2-1

動力ゴムのループ長と条数の定義 ............................................................................................................. 10

図 3.2-2ゴムの巻き数とトルクの関係 ..................................................................................................................... 10

図 3.3-1

巻き数-巻き戻し時トルクのグラフの例 .................................................................................................... 11

図 3.3-2

巻き込み数と最大トルクの関係の例 ......................................................................................................... 12

図 3.3-3

巻き数とトルクのモデル ............................................................................................................................. 12

図 4.1-1

プロペラの形状と寸法 ................................................................................................................................. 13

図 4.2-1

機体固定座標系とプロペラ座標系 ............................................................................................................. 14

図 4.3-1

プロペラの運動............................................................................................................................................. 15

図 4.3-2プロペラ性能曲線の例 ................................................................................................................................. 16

図 4.3-3

プロペラブレードの断面 ............................................................................................................................. 17

図 4.3-4

プロペラ後流................................................................................................................................................. 17

図 4.3-5ブレードの翼素に作用する力 ................................................................................................................... 19

図 4.3-6

失速後の２次元翼空力データの推定例（参考文献[9]より） ................................................................. 21

図 4.3-7プロペラと気流の角度 ................................................................................................................................. 22

図 4.3-8ブレード翼素と空気流の関係 ..................................................................................................................... 25

図 4.3-9

横力の方向 .................................................................................................................................................... 26

図 4.3-10プロペラの性能曲線の例 ........................................................................................................................... 27

図 4.4-1

プロペラの慣性能率..................................................................................................................................... 28

図 4.4-2プロペラのトルクの釣り合い ..................................................................................................................... 29

図 4.4-3プロペラのモーメントの釣り合い ............................................................................................................. 29

図 4.4-4

ゴムとプロペラ軸の成す角 ......................................................................................................................... 30

図 4.4-5プロペラのジャイロ効果 ............................................................................................................................. 31

図 6.1-1

解析対象の機体............................................................................................................................................. 33

図 6.1-2

解析対象の機体の図面 ................................................................................................................................. 34

図 6.1-3

水平尾翼容積と重心位置の関係 ................................................................................................................. 38

図 6.2-1

AVL モデルの座標系 .................................................................................................................................... 39

図 6.2-2

AVL モデル .................................................................................................................................................... 40

図 6.2-3

主翼と水平尾翼の揚力係数 ......................................................................................................................... 42

図 6.2-4

失速後の全機揚力係数と全機抗力係数 ..................................................................................................... 43

図 6.4-1

使用する動力ゴムの巻き数-トルク関係 .................................................................................................... 45

図 6.4-2

巻き数とトルクの関係 ................................................................................................................................. 47

図 6.5-1

ピッチ比 1.4 のプロペラ .............................................................................................................................. 48

図 6.5-2翼弦長分布 .................................................................................................................................................... 49

図 6.5-3

２次元翼の揚力係数と抗力係数データ – 薄肉円弧翼（キャンバー5%），レイノルズ数 20,000[23].................................................................................................................................................................................. 49v

図 6.5-4プロペラ性能曲線......................................................................................................................................... 51

図 6.5-5プロペラ性能曲線の近似式に使ったニューラルネットワーク.............................................................. 52

図 6.5-6

Ribner の理論と McCormic の理論の比較 .................................................................................................. 53

図 6.6-1

機体固定座標系............................................................................................................................................. 55

図 6.6-2３次元運動解析結果 – ケース 1（ベースライン） ................................................................................ 60

図 6.6-3３次元運動解析結果 – ケース 1a，プロペラに気流が斜めにあたる影響を無視 ............................... 61

図 6.6-4３次元運動解析結果 – ケース 1b，プロペラのジャイロ効果を無視 ................................................... 62

図 6.6-5３次元運動解析結果 – ケース 1c，プロペラの回転慣性を無視 ........................................................... 63

図 6.6-6３次元運動解析結果 – ケース 2，プロペラのアップスラスト 1.0 度 .................................................. 64

図 6.6-7３次元運動解析結果 – ケース 3，プロペラのダウンスラスト 1.0 度 .................................................. 65

図 6.6-8３次元運動解析結果 – ケース 4，プロペラのダウンスラスト 5.0 度 .................................................. 66

図 6.6-9３次元運動解析結果 – ケース 5，右ラダー舵角有り，サイドスラスト左 0.8 度 .............................. 67

図 6.6-10３次元運動解析結果 – ケース 6，ラダー舵角無し，サイドスラスト無し ....................................... 68

図 6.6-11３次元運動解析結果 – ケース 7，左ラダー舵角有り，サイドスラスト左 0.5 度 ............................ 69

図 6.6-12３次元運動解析結果 – ケース 8，左ラダー舵角有り，サイドスラスト右 1.0 度 ............................ 70

図 6.6-13３次元運動解析結果 – ケース 9，発進角度 25 度 ............................................................................... 71

図 6.6-14３次元運動解析結果 – ３次元運動解析結果 – ケース 10，発進角度 35 度 .................................... 72

図 6.6-15

３次元運動解析結果 – ３次元運動解析結果 – ケース 11，発進角度 55 度 .................................... 73

図 6.6-16３次元運動解析結果 – ３次元運動解析結果 – ケース 12，発進角度 65 度 .................................... 74

図 6.6-17３次元運動解析結果 – ３次元運動解析結果 – ケース 13，FAI TAN ゴム，発進角度 45 度 ........ 75

図 6.6-18

３次元運動解析結果 – ３次元運動解析結果 – ケース 14，FAI TAN ゴム，発進角度 50 度 ........ 76

図 6.6-19３次元運動解析結果 – ３次元運動解析結果 – ケース 15，FAI TAN ゴム，発進角度 60 度 ........ 77

図 6.6-20

３次元運動解析結果 – ３次元運動解析結果 – ケース 16，FAI TAN ゴム，発進角度 70 度 ........ 78

図 6.6-21

３次元運動解析結果 – ３次元運動解析結果 – ケース 17，FAI TAN ゴム，発進角度 75 度 ........ 79

図 6.6-22

プロペラのパワー効果 ............................................................................................................................... 80vi

表のリスト
表 6.1-1

解析対象の機体の諸元 ................................................................................................................................. 35

表 6.1-2

主翼面積，平均空力翼弦の計算 ................................................................................................................. 36

表 6.1-3

水平尾翼容積，垂直尾翼容積の計算 ......................................................................................................... 37

表 6.2-1

全機 AVL 解析結果 ....................................................................................................................................... 41

表 6.2-3

運動解析に使用する空力微係数 ................................................................................................................. 41

表 6.3-1

機体の慣性データの計算結果 ..................................................................................................................... 44

表 6.4-1

巻き数-トルク関係式 ................................................................................................................................... 46

表 6.5-1

プロペラ性能計算に使用する２次元翼の空力データ ............................................................................. 50

表 6.5-2プロペラの Solidity の計算 .......................................................................................................................... 50

表 6.5-3プロペラ性能曲線の近似式のパラメータ ................................................................................................. 51

表 6.5-4

プロペラ慣性能率の計算 – プロペラ直径 170 mm の場合 .................................................................... 54

表 6.6-1

プロペラの効果............................................................................................................................................. 57

表 6.6-2

ダウンスラスト角の影響 ............................................................................................................................. 57

表 6.6-3

サイドスラスト角の影響 ............................................................................................................................. 58

表 6.6-4

発進角度の影響............................................................................................................................................. 58

表 6.6-5

発進角度の影響，FAI TAN Super Sport ゴム，3/32 インチ幅を使用..................................................... 59

vii

1

研究の目的
フリーフライト模型飛行機の一種であるゴム動力模型飛行機（以下，ゴム動力機という）は模型飛行機のうちで最も簡単な動力付き飛行機であり，模型飛行機の初心者が製作して飛ばすことが多い．しかし，ゴム動力機をうまく飛ばすのは非常に難しい．うまく飛ばないために，模型飛行機に対する興味を失うケースも多い．難しい理由は，調整箇所が多く，しかも試行錯誤による調整が必要なことにある．機体のどこを，どういじったら，どういう飛び方をするかがわからないと調整ができない．
残念なことに，ゴム動力機の３次元運動解析はこれまで行われておらず，飛行調整に関する理論的な知見は得られていない．筆者は，ハイスタート発進の模型グライダーについて，３次元運動解析を行ってその調整方法に関する知見を明らかにした．その解析法を拡張すれば，ゴム動力機の３次元運動解析も可能であると考えた．
本研究の目的は，ゴム動力機の３次元運動解析法を開発し，その解析法を使ってゴム動力機の合理的な飛行調整法を見出すことである．

1.1 ゴム動力機の概要
最も単純な構造のゴム動力機の実例を図 1.1-1 に示す．胴体は１本の棒でできており，前方には主翼が固定され，尾部には水平尾翼と垂直尾翼が固定されている．胴体の先端にプロペラ軸を保持するプロペラ・ハンガーが取り付
けられている．プロペラ・ハンガーの軸受けには，プロペラ軸が通っており，プロペラが取り付けられている．プ
ロペラ軸の後方はフック状になっており，ループ状にした動力ゴムの前端を取り付ける．胴体の後部下側にゴム・フックが固定されており，動力ゴムの後端を取り付ける．
プロペラは機体の後方から見て，右回りに回転して推力を出すようになっている．プロペラ・ハブの前側には，空転ラッチ機構がついており，ゴムが完全に巻き戻ったあとはプロペラが空転して空気抵抗を減らすようになっている．
プロペラ・ハンガーは太めの針金とそれに固定された軸受けから成り，針金を曲げることにより，プロペラ軸の角度を微調整する．このプロペラ軸の角度はゴム動力機の調整に非常に重要である．
ゴム動力機には，室内で飛ばす特殊な機種（重量数グラム程度）や屋外での国際競技用の大型で高性能な機種もあるが，本研究では製作のしやすさを考慮して，図 1.1-1 に示すような翼幅 400mm 程度の初心者向け小型機体を研究対象とする．

図 1.1-1

ゴム動力機の実例

1

1.2 解析で考慮する項目
本研究の解析で考慮する項目は以下のとおりである．


ゴムの巻き数-トルク曲線



プロペラの影響

プロペラの推力とトルク



プロペラに斜めに気流があたる影響

プロペラ軸の位置，偏心，向き（ダウンスラスト，サイドスラスト）



動力飛行後のプロペラ空転による抗力



プロペラの回転慣性の影響



プロペラのジャイロ効果

垂直尾翼の傾き（垂直尾翼後縁に貼るガーニーフラップの効果を模擬）

エルロン舵角（主翼外舷部の後縁に貼るガーニーフラップの効果を模擬）

水平尾翼の左右方向の傾き（いわゆるスタブ・ティルト）

考慮していないのは次の項目である．


プロペラ後流が主翼，尾翼にあたる影響

2

2

模型飛行機の３次元運動方程式
参考文献[2] を参考にしてゴム動力模型飛行機の３次元運動方程式を示す．本レポートでは微小擾乱の仮定は用いていない．

2.1 機体の姿勢と速度
飛行中の機体の姿勢と速度を図 2.1-1 に示す．地面固定座標系を XE-YE-ZE とし，機体固定座標系（重心を原点とする）を X-Y-Z とする．機体速度 VC の機体固定座標系に関する成分を U，V，W とする．
機体の重心には，重力，空気力等の外力 FX，FY，FZ とモーメント𝐿，M，N が作用する．
（図 2.1-2）．
機体固定座標系と地面固定座標系の関係を表すオイラー角の定義を図 2.1-3 に示す．

図 2.1-1

機体の姿勢と速度

図 2.1-2

機体に作用する力3図 2.1-3

地面固定座標系と機体固定座標系の関係，オイラー角の定義

オイラー角を使うと，機体固定座標系と地面固定座標系の関係は次のように表すことができる．
𝑋𝑋
𝑌 = 𝐸(𝛷, 𝛩, 𝛹) 𝑌𝑍𝑍
cos 𝛩 cos 𝛹= sin 𝛷 sin 𝛩 cos 𝛹 − cos 𝛷 sin 𝛹
cos 𝛷 sin 𝛩 cos 𝛹 + sin 𝛷 sin 𝛹

cos 𝛩 sin 𝛹
sin 𝛷 sin 𝛩 sin 𝛹 + cos 𝛷 cos 𝛹
cos 𝛷 sin 𝛩 sin 𝛹 − sin 𝛷 cos 𝛹− sin 𝛩sin 𝛷 cos 𝛩
cos 𝛷 cos 𝛩𝑋𝑌𝑍(2.1 − 1)

2.2 運動方程式
機体の重心まわりの機体固定座標系における運動方程式は次のように表される．
𝑚 𝑈̇ + 𝑄𝑊 − 𝑅𝑉 = 𝐹(2.2 − 1)𝑚 𝑉̇ + 𝑅𝑈 − 𝑃𝑊 = 𝐹(2.2 − 2)𝑚 𝑊̇ + 𝑃𝑉 − 𝑄𝑈 = 𝐹(2.2 − 3)𝐼 𝑃̇ − 𝐼 𝑅̇ + (𝐼𝐼 𝑄̇ + (𝐼− 𝐼 )𝑄𝑅 − 𝐼 𝑃𝑄 = 𝐺(2.2 − 4)− 𝐼 )𝑅𝑃 + 𝐼 (𝑃 − 𝑅 ) = 𝐺(2.2 − 5)−𝐼 𝑃̇ + 𝐼 𝑅̇ + (𝐼− 𝐼 )𝑃𝑄 + 𝐼 𝑄𝑅 = 𝐺(2.2 − 6)(2.2-4)～(2.2-6)式を P，Q，R に
̇ ついて解いて，

4

𝑃̇ =

𝑄̇ =
𝑅̇ =

1
𝐼 𝐼

−𝐼

1
[𝐺 − (𝐼
𝐼
1
𝐼 𝐼

−𝐼

−𝐼

+𝐼 𝐼

−𝐼

𝑄𝑅 + 𝐼 (−𝐼+𝐼+ 𝐼 )𝑃𝑄 + 𝐼 𝐺+𝐼 𝐺(2.2 − 4𝑎)− 𝐼 )𝑅𝑃 − 𝐼 (𝑃 − 𝑅 )](2.2 − 5𝑎)

𝐼 (−𝐼

+𝐼

− 𝐼 )𝑄𝑅 + −𝐼 𝐼

+𝐼

+𝐼

𝑃𝑄 + 𝐼 𝐺+𝐼 𝐺(2.2 − 6𝑎)

ここで，FX, FY, FZ, GX, GY, GZ は重心まわりの外力とモーメント，
U は機軸方向（X 方向）の速度，V は横方向（Y 方向）の速度，W は機軸に垂直方向下向き（Z 方向）の速度，P はローリング角速度，Q はピッチング角速度，R はヨーイング角速度，m は機体質量，IXX は重
心まわりのローリング慣性能率，IYY は機体の重心まわりのピッチング慣性能率，IZZ は重心まわりのヨー
イング慣性能率，IXZ は重心まわりの慣性乗積である．
記号の上のドット「・」は時間微分を表す．
オイラー角と機体の角速度の関係は次のように表される．
𝑃 = 𝛷̇ − 𝛹̇ sin 𝛩(2.2 − 7)

𝑄 = 𝛩̇ cos 𝛷 + 𝛹̇ sin 𝛷 cos 𝛩(2.2 − 8)

𝑅 = −𝛩̇ sin 𝛷 + 𝛹̇ cos 𝛷 cos 𝛩(2.2 − 9)オイラー角の時間微分について解くと，
𝛷̇ = 𝑃 + 𝑄 sin 𝛷 tan 𝛩 + 𝑅 cos 𝛷 tan 𝛩(2.2 − 7𝑎)𝛩̇ = 𝑄 cos 𝛷 − 𝑅 sin 𝛷(2.2 − 8𝑎)𝛹̇ = 𝑄 sin 𝛷 sec 𝛩 + 𝑅 cos 𝛷 sec 𝛩(2.2 − 9𝑎)

速度ベクトル成分と迎角，横滑り角 の関係は次のようになる（図 2.1-2）．
tan 𝛼 =

𝑊
𝑈

(2.2 − 10)tan 𝛽 =

𝑉𝑈

(2.2 − 11)

地上固定座標系における機体の重心位置を（xe, ye, ze）とすると，重心位置の速度は次のようになる．
𝑑𝑥= 𝑈 cos 𝛩 cos 𝛹 + 𝑉(sin 𝛷 sin 𝛩 cos 𝛹 − cos 𝛷 sin 𝛹)𝑑𝑡+ 𝑊(cos 𝛷 sin 𝛩 cos 𝛹 + sin 𝛷 sin 𝛹)(2.2 − 12)𝑑𝑦= 𝑈 cos 𝛩 sin 𝛹 + 𝑉(sin 𝛷 sin 𝛩 sin 𝛹 − cos 𝛷 cos 𝛹)𝑑𝑡+ 𝑊(cos 𝛷 sin 𝛩 sin 𝛹 − sin 𝛷 cos 𝛹)𝑑𝑧= −𝑈 sin 𝛩 + 𝑉 sin 𝛷 cos 𝛩 + 𝑊 cos 𝛷 cos 𝛩𝑑𝑡(2.2 − 13)(2.2 − 14)5

2.3 機体に作用する空気力
機体の重心に作用する機体固定座標系方向の空気力 Xa，Ya，Za とモーメント𝐿，M，N は次のように表すことができる．
𝑋 = 𝐿 sin 𝛼 − 𝐷 cos 𝛼(2.3 − 1)

1
𝑌 = 𝜌(𝑈 + 𝑉 + 𝑊 )𝑆 𝐶
2

(2.3 − 2)

𝑍 = −𝐿 cos 𝛼 − 𝐷 sin 𝛼(2.3 − 3)

1
𝐿 = 𝜌(𝑈 + 𝑉 + 𝑊 )𝑏𝑆 𝐶
2(2.3 − 4)1
𝑀 = 𝜌(𝑈 + 𝑉 + 𝑊 )𝑐̅𝑆 𝐶
2(2.3 − 5)1
𝑁 = 𝜌(𝑈 + 𝑉 + 𝑊 )𝑏𝑆 𝐶
2(2.3 − 6)

ここで，L，D は全機の揚力と抗力で次のように表される．
1
𝐿 = 𝜌(𝑈 + 𝑉 + 𝑊 )𝑆 𝐶
2

(2.3 − 7)

1
𝐷 = 𝜌(𝑈 + 𝑉 + 𝑊 )𝑆 𝐶
2

(2.3 − 8)

：空気密度，Sw：主翼面積，𝑐̅：主翼の平均空力翼弦，b：翼幅，CL：全機揚力係数，
CD：全機抗力係数，Cy：全機横力係数，Cl：全機ローリングモーメント係数，
Cm：全機ピッチングモーメント係数，Cn：全機ヨーイングモーメント係数
全機の空気力係数は空力微係数を使って次のように近似できる（参考文献[3]）．垂直尾翼の傾き，水平尾翼の傾き，
エルロンの影響を含めた．
𝐶 = 𝐶 (𝛼 − 𝛼 ) − 𝐶 𝑞(2.3 − 9)

𝐶 =𝐶(2.3 − 10)+ 𝑘𝐶

𝐶 = 𝐶 (𝛽 − 𝑖 ) + 𝐶 𝑝 + 𝐶 𝑟 + 𝑎
𝐶 = 𝐶 𝛽+𝐶 𝑝+𝐶 𝑟+𝐶
𝐶 =𝐶

𝛼−𝛼

+𝐶

𝑆
𝑆

𝛼 −𝜀 −

𝜕𝜀𝛼 − −𝑖 𝜕𝛼𝛿(2.3 − 12)𝑞

𝐶 = 𝐶 𝛽+𝐶 𝑝+𝐶 𝑟+(2.3 − 11)(2.3 − 13)𝑙𝑙
𝑆
𝐶 𝑖 − 𝑎𝑏𝑏 𝑆

𝛼 −𝜀 −

𝜕𝜀𝛼 −𝑖 
𝜕𝛼

(2.3 − 14)

ここで，lf：機体重心から垂直尾翼空力中心までの距離，

a：エルロン舵角，iv：垂直尾翼取付角（図 2.3-1），th：水平尾翼の傾き（図 2.3-2，右翼下げが正），
lt：機体重心から水平尾翼空力中心までの距離，
it：水平尾翼取付角＝主翼無揚力線からの水平尾翼の角度（図 2.3-3，後縁上げが正），St：水平尾翼面積，at：水平尾翼揚力傾斜，w：主翼無揚力線からの主翼の迎角，

：吹き下ろし角，0：主翼揚力がゼロのときの吹き下ろし角，
0：ゼロ揚力迎角，eq：釣り合い迎角，CDP：有害抵抗係数，k：誘導抵抗の係数，
p，q，r：無次元角速度で，次のように定義している．
𝑝=

𝑏𝑃
2𝑈

(2.3 − 15)

6

𝑞=

𝑐̅𝑄
2𝑈

(2.3 − 16)

𝑟=

𝑏𝑅
2𝑈

(2.3 − 17)

ガーニーフラップを垂直尾翼後縁左側につけたときが＋
図 2.3-1

垂直尾翼の取付角

ガーニーフラップを主翼左舷後縁下側につけたときが＋
図 2.3-2

図 2.3-3

水平尾翼の傾きとエルロン舵角の符号

主翼単独の無揚力線と水平尾翼取付角の関係

7

(2.3-11)式と(2.3-14)式に出てくる
𝑎 𝛼 −𝜀 −𝜕𝜀𝛼 − −𝑖𝜕𝛼は，水平尾翼荷重係数であり，次のように表すことができる．
(𝐶 )=𝑎 𝛼 −𝜀 −ここで， 𝐶𝜕𝜀𝛼 − −𝑖
𝜕𝛼

= 𝐶

𝛼−𝛼 ,

：水平尾翼の揚力傾斜，𝛼 ,(2.3 − 18)：水平尾翼のゼロ揚力迎角(2.3-11)式と(2.3-14)式は次のように書き換えることができる．
𝐶 = 𝐶 (𝛽 − 𝑖 ) + 𝐶 𝑝 + 𝐶 𝑟 +𝐶 = 𝐶 𝛽+𝐶 𝑝+𝐶 𝑟+

𝑆
𝑆

𝐶

𝑙
𝑙 𝑆
𝐶 𝑖 −
𝑏
𝑏𝑆



𝛼−𝛼 ,
𝐶

𝛼−𝛼 ,

(2.3 − 11𝑎)

(2.3 − 14𝑎)

2.4 空力微係数の計算方法
前項で出てきた機体の空力微係数を計算するには，渦格子法プログラム AVL [19], [20] を使用する．

8

3

動力ゴム

3.1 使用するゴム
動力ゴムにはゴム動力飛行機専用の高性能ゴムがある．一般的には専用のゴムを使用するが，本研究では入手が容易で，安価な輪ゴムも使用する．輪ゴムは専用のゴムに比べて伸びが半分程度で性能は落ちるが，大量生産される工業製品であるため，物理特性が安定していると考えられる．
輪ゴムを使う場合に気を付けなければならない点は，


ワインダーを使って伸ばして巻くとトルクが低下する．



つなぐために結び目を作ると「こぶ」が均一にできず，振動が起きやすい．

そこで，以下の対策をとる．


手巻きとする．または，ワインダーを使っても伸ばさずに巻く．



結び目を作らないでつなぐ．結び目を作らないつなぎ方を図 3.1-1 に示す．

図 3.1-1

結び目を作らない輪ゴムのつなぎ方

3.2 巻き数とトルクの関係
動力ゴムのループ長と条数の定義を図 3.2-1 に示す．動力ゴムの巻き数とトルクの関係を測定した例を図 3.2-2に示す．このグラフからわかるように，巻き込んでいくときのトルクと巻き戻すときのトルクは大きく異なり，ヒステリシス現象がある．ゴム動力機を飛ばすときにプロペラに伝達されるトルクは巻き戻しのときのトルクである．巻き戻し初めの４０回程度で急激にトルクが低下し，その後はトルクが緩やかに低下していく．急激にトルクが減少する領域を「バースト領域」という．その後のほぼ一定トルクの領域を「クルーズ領域」という．

9

ループ長 ＝ ゴムの全長÷条数
輪ゴムの場合は，輪ゴムを切断して全本数をつなぎ合わせた長さを条数で割った値をループ長とする．
図 3.2-1

動力ゴムのループ長と条数の定義

0.025#18輪ゴム，結び目無し，潤滑あり，手巻き
輪ゴム本数：18本
条数：8本
ループ⾧：315 mm
フック間距離：278 mm
平均破断巻き数：525回

0.02

トルク
(N-m)
0.015

巻き込み
0.01

巻き戻し

0.005

バースト領域
クルーズ領域
0
0

50

100

150

200

250

300

巻き数（手巻き）

図 3.2-2ゴムの巻き数とトルクの関係

10

350

400

450

3.3 巻き数とトルクのモデル
ゴム動力機用ゴム FAI Tan Super Sport のトルクデータの例を図 3.3-1 に示す．トルクはクルーズトルクと付加トルクの和であり，クルーズトルクは巻き込み数に関わらず同じである．このグラフから得た巻き込み数と最大トルクの関係を図 3.3-2 に示す．
これらのグラフから，巻き数とトルクの関係は図 3.3-3 に示すモデルで表すことができる．モデルに使用する関数はゴムの種類・条数ごとの計測データ（例：図 3.3-1）を使って設定する．
ゴムの巻き数を次のように表す．
𝑚

=𝑚 −𝑚

(3.3 − 1)

ここで，mrubber は巻き数，ms は巻き込み数，mprop はプロペラ回転数
プロペラの単位時間あたりの回転数 n は，
𝑛=

𝑑𝑚
𝑑𝑡

(3.3 − 2)

ゴムのトルク Qr の具体的な式は 6.4 項に示す．

25

FAI TAN Super Sport
3/32 inch幅，4条

20

最大トルク
トルク
(N-mm)

15

10

付加トルク
5

クルーズトルク
0
0

0.5

1

1.5

2

単位⾧さあたり巻き数 （回/mm）

図 3.3-1

巻き数-巻き戻し時トルクのグラフの例

11

2.5

3

25

FAI TAN Super Sport
3/32 inch幅，4条

20

最大トルク 15(N-mm)

10

5

0
0

0.5

1

1.5

2

単位⾧さあたり巻き込み数 （回/mm）

図 3.3-2

巻き込み数と最大トルクの関係の例

図 3.3-3

巻き数とトルクのモデル

12

2.5

3

4

プロペラ

4.1 プロペラの形状と寸法の定義
図 4.1-1 に示すように，プロペラの直径を Dp，プロペラ軸からの距離（半径）r の位置の翼素の翼弦長を c，翼
角を とする．半径 r の位置での幾何学的ピッチ P は翼角を使って次のように表される．
𝑃 = 2𝜋𝑟 tan 𝛽(4.1 − 1)

図 4.1-1

プロペラの形状と寸法134.2 機体固定座標系とプロペラ座標系の関係
プロペラ座標系と機体座標系の関係を図 4.2-1 に示す．プロペラ軸方向の単位ベクトル iprop の成分（機体固定座標系）は，
𝐢

𝑖
= 𝑖
𝑖

=

cos 𝜃 cos 𝜃
cos 𝜃 sin 𝜃− sin 𝜃(4.2 − 1)

ここで，up はアップスラスト角，side はサイドスラスト角

プロペラの中心（Xprop, Yprop, Zprop）
アップスラスト角度
機体固定座標系

プロペラ軸方向の単位ベクトル

右サイドスラスト角度

プロペラ座標系

図 4.2-1

プロペラ軸方向の単位ベクトル（機体固定座標系）

機体固定座標系とプロペラ座標系

144.3 プロペラの空気力
4.3.1

プロペラの性能曲線の計算法

Larrabee の方法を使ってプロペラの性能曲線を計算する方法を示す．
4.3.1.1 プロペラの性能曲線
図 4.3-1 の上の図に示すような直径 Dp = 2R のプロペラを考える．ブレード枚数を B とする．このプロペラが回転角速度  （rad/s）で回転しており，プロペラ（機体）の前進速度が V であるとする．
プロペラの回転速度 n（回/s）と回転角速度の関係は，
𝛺 = 2𝜋𝑛(4.3 − 1)𝛺
2𝜋(4.3 − 2)𝑛=前進係数を次のように定義する．
𝐽=

𝑉
2𝜋𝑉
=
𝑛𝐷
𝛺𝐷

(4.3 − 3)プロペラが発生する推力 T とトルク Q は無次元係数 CT と CQ を使って次のように表すことができる．
𝑇 = 𝐶 𝜌𝑛 𝐷

(4.3 − 4)

𝑄 = 𝐶 𝜌𝑛 𝐷

(4.3 − 5)

ここで， は空気密度である．
プロペラ効率p は次のように表される．
𝜂 =

𝑇𝑉
𝑄𝛺

(4.3 − 6)

図 4.3-1

プロペラの運動15

無次元係数 CT と CQ は前進係数 J の関数であり，図 4.3-2 に示すようなグラフで表示される．これをプロペラの性能曲線という．レイノルズ数，マッハ数の影響が無いと仮定すれば，ひとつのプロペラの性能曲線があれば，相似形のプロペラの推力とトルクはこの無次元係数を使って容易に算出することができる．

0.8

p

0.7

0.6

CT
0.5

10×CQ

p

0.4

10×CQ

0.3

0.2

CT

0.1

0

-0.1
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

1.5

前進係数 J
図 4.3-2プロペラ性能曲線の例

4.3.1.2 Larrabee の方法
参考文献[5]，[6]にしたがって，プロペラの性能を計算する方法を示す．
（１）ブレードの翼素の空気流の速度
半径 r の位置にあるブレード断面を考える．図 4.3-3 に示すように，弦長 c の翼素（翼角）を考え，この翼素にあたる空気の速度ベクトル w の成分を次のように表す．
𝑤 = 𝑉(1 + 𝑎)(4.3 − 7)𝑤 = 𝛺𝑟(1 + 𝑎′)(4.3 − 8)𝑤=(4.3 − 9)𝑤tan  =+𝑤

𝑤
𝑉(1 + 𝑎)=𝑤
𝛺𝑟(1 + 𝑎′)(4.3 − 10)

ここで，a は軸方向誘導速度係数，a’ は接線方向誘導速度係数である．

16

図 4.3-3

プロペラブレードの断面（２）プロペラ後流の速度
プロペラ後流の軸方向の速度の変化を v’ とする．軸方向速度変化率  を次のように定義する．
𝜁=

𝑣′
𝑉

(4.3 − 11)

プロペラ後流の速度変化と翼素の接線方向速度との関係は図 4.3-4 より
𝑤 = 𝑤 sin  = 𝑣′ sin  cos (4.3 − 12)

半径 r の単位幅円環のブレード１枚あたりの循環を とすると，プロペラの円環の循環は
𝐵𝛤 = 2𝜋𝑟𝑤 𝐹 = 2𝜋𝑟𝑣′ sin  cos  𝐹(4.3 − 13)

ここで，F はプラントルの翼端損失係数と呼び，翼端渦の影響を考慮する係数であり，次のように定義され
る．
𝐹=

2
cos
𝜋

𝑒

(4.3 − 14)

ここで，f は V，，r で決まる係数，e は自然対数の底

図 4.3-4

プロペラ後流

この循環に現れる角度 を近似的に表す（図 4.3-4 で a と a’ を小さいとして無視）と次のようになる．
sin  =

1

(4.2 − 15)√1 + 𝑥

ここで，𝑥 = 𝛺𝑟⁄𝑉

17

これらの式を使って循環の式を書き直すと，
𝐵𝛤𝛺
𝑥
=
𝐹=𝐺
2𝜋𝑣′ 1 + 𝑥(4.2 − 16)

f は次のように表される．
𝑓=𝐵 √𝜆 + 1(1 − 𝜉)2 𝜆(4.2 − 17)

ここで，𝜆 = 𝑉⁄𝛺𝑅 , 𝜉 = 𝑟⁄𝑅
プロペラの円環の循環をさらに書き換えると，
𝛤=

2𝜋𝑉𝑣′𝐺 2𝜋𝐺
=
𝑉 𝜁
𝐵𝛺
𝐵𝛺

(4.2 − 18)（３）ブレードに作用する力
ブレードの翼素に作用する力を図 4.3-5 に示す．
ブレード単位幅に作用する推力 T’ とトルク Q’ は次のように表される．
𝑑𝑇= 𝐿 cos  − 𝑇 sin  = 𝐿 cos  (1 − 𝜖 tan )𝑑𝑟𝑇 =(4.3 − 19)

𝑄
𝑑𝑄
𝜖
=
= 𝐿 sin  1 +
𝑟
𝑟𝑑𝑟
tan 

(4.3 − 20)

ここで，ϵ は２次元翼の揚抗比の逆数である．
𝜖=𝑐𝑐(4.3 − 21)

循環と揚力の関係（クッター・ジューコフスキーの定理）を使うと，上の式は次のように書き換えられる．
𝑇 = 𝐵𝜌𝑤𝛤 cos  (1 − 𝜖 tan ) = 𝐵𝜌𝛤𝛺𝑟(1 − 𝑎′)(1 − 𝜖 tan )(4.3 − 22)𝑄′𝜖𝜖= 𝐵𝜌𝑤𝛤 sin  1 += 𝐵𝜌𝛤𝑉(1 + 𝑎) 1 +
𝑟
tan 
tan 

(4.3 − 23)

接線方向誘導速度係数 a’ を後流 wt の 1/2 であると仮定すると，
𝑎 =1𝑤
1 𝑣′ sin  cos  1
1== 𝜁
2 𝛺𝑟 2𝛺𝑟
2 1+𝑥(4.3 − 24)

この式を推力の式に代入して
1
1
𝑇 = 𝐵𝜌𝛤𝛺𝑟 1 − 𝜁
2 1+𝑥(1 − 𝜖 tan )(4.3 − 25)

さらに書き換えると
1
1
𝑇 = 2𝜋𝑟𝜌𝐺𝑉 𝜁 1 − 𝜁
2 1+𝑥

1−

𝜖
𝑥

(4.3 − 26)

同様に，
1
1
𝑥
𝑎 = 𝑣 (cos ) = 𝜁2
2 1+𝑥(4.3 − 27)

この式をトルクの式に代入して
𝑄′1𝑥= 𝐵𝜌𝛤𝑉 1 + 𝜁
𝑟
2 1+𝑥

1+

𝜖
tan 

(4.3 − 28)

さらに書き換えると𝑄′ 2𝜋𝜌𝐺𝑉 𝜁
1
𝑥
=
1+ 𝜁
𝑟
𝛺
2 1+𝑥

1+

𝜖
tan 

(4.3 − 29)18図 4.3-5ブレードの翼素に作用する力（４）推力と後流の軸方向速度との関係
推力 T から後流の軸方向速度 v’ を求めるには，推力をプロペラ半径で微分して
𝑇 =

𝑑𝑇 1
𝑑𝑇
= 𝜌𝑉 𝜋𝑅𝑑𝑟 2𝑑𝜉(4.3 − 30)

ここで，Tc は推力係数で
𝑇 =

2𝑇
𝜌𝑉 𝜋𝑅

(4.3 − 31)

推力の式を代入して次の式が得られる．
𝑑𝑇𝜖= 4𝜉𝐺𝜁 1 − − 2𝜉𝐺𝜁
𝑑𝜉
𝑥

1−𝜖
𝑥 +1

(4.3 − 32)

この式をプロペラの半径方向に積分して， について解くと次の式が得られる．
𝜁=

𝐼
2𝐼

1− 1−

4𝐼 𝑇

(4.3 − 33)

𝐼

ここで，I1，I2 は積分で，次のように表されるので数値積分することができる．
𝜖
𝑑𝜉
𝑥

𝐼 =

4𝜉𝐺 1 −

𝐼 =

𝜖
𝑥
2𝜉𝐺
𝑑𝜉
𝑥 +1

(4.3 − 34)

1−(4.3 − 35) が求められると，a，a’ を計算できるので，流れの角度 を次の式で計算でき，迎角もわかる．
 = tan𝑉(1 + 𝑎)𝛺𝑟(1 − 𝑎′)(4.3 − 36)

19

（５）プロペラの性能計算手順
プロペラの翼角分布，翼弦長分布が与えられたときにプロペラの性能を計算する手順は以下のとおり．
① プロペラの寸法（R，翼角分布，翼弦長分布 c）を与える．
② 前進速度 V と回転角速度 を与える．
③ 流れの角度を仮定する． = tan𝑉⁄𝛺𝑟 を使えばよい．

④ 迎角 を計算する．𝛼 = 𝛽 − ⑤ ２次元翼の揚力係数 cl ，抗力係数 cd を与える．
⑥ プラントルの翼端損失係数を計算する．
⑦ 揚力係数と抗力係数を推力係数とトルク係数に座標変換する．
⑧ a，a’ を計算し，速度 w を計算する．
⑨ 新しい流れ角度new を計算する．
⑩ 流れ角度を更新して，④に戻って再計算を行う．
⑪ 流れ角度が収束したら，推力，トルク，パワー，効率を計算する．
必要な計算式を以下に示す．
𝑐 = 𝑐 (cos  − 𝜖 sin )(4.3 − 37)𝑐 = 𝑐 (sin  + 𝜖 cos )(4.3 − 38)

1
𝑇 = 𝜌𝑤 𝐵𝑐𝑐
2(4.3 − 39)𝑄′ 1= 𝜌𝑤 𝐵𝑐𝑐𝑟2(4.3 − 40)𝜎𝑐4𝐹(sin ) − 𝜎𝑐(4.3 − 41)𝜎𝑐
4𝐹 sin  cos  + 𝜎𝑐

(4.3 − 42)

ここで，𝜎 = 𝐵𝑐 ⁄(2𝜋𝑟)(4.3 − 43)

𝑎=

𝑎 =

前進速度と回転角速度を変化させて推力とトルクを計算することにより，プロペラの性能曲線を得ることができる．

20

（６）２次元翼の空力データ
プロペラの性能計算の入力データとして，失速後の大きな迎角までの２次元翼の揚力係数と抗力係数を必要とする．ゴム動力機用プロペラには，低レイノルズ数における２次元翼の風洞試験データを使用する．失速後のデー
タは発表されていないので，推定値を使用する．失速後の迎角における２次元翼の空力データを推定する方法としては Viterna and Corrigan の方法（参考文献[8]）がある．
失速前𝑐 = 𝑐 (𝛼 − 𝛼 )(4.3 − 44)𝑐 =𝑐(4.3 − 45)



+ 𝑘𝑐

失速後
(cos 𝛼)sin 𝛼(4.3 − 46)𝑐 = 𝐵 (sin 𝛼) + 𝐵 cos 𝛼(4.3 − 47)𝑐 = 𝐴 sin 2𝛼 + 𝐴𝐴 =𝑐(4.3 − 48)

2

𝐴 = 𝑐

−𝑐

sin 𝛼

cos 𝛼

sin 𝛼
(cos 𝛼(4.3 − 49)

)

𝐵 =𝑐
𝐵 =
𝑐

𝑐

(4.3 − 50)(sin 𝛼−𝑐cos 𝛼)(4.3 − 51)= 1.11 + 0.018𝐴𝑅(4.3 − 52)

ここで，AR はブレードのアスペクト比
ただし，低レイノルズ数の場合は(4.3-52)式が成り立たないので，失速後の空力データとのマッチングを行って cdmax を決める必要がある．

図 4.3-6

失速後の２次元翼空力データの推定例（参考文献[9]より）

21

4.3.2気流がプロペラの面に正対していない場合に作用する力

プロペラが回転角速度  で回転しており（後ろから見て時計回り方向の回転とする），空気流の速度が V であるとする．気流はプロペラの面に対して角度P だけ傾いているとする（図 4.3-7）．

図 4.3-7プロペラと気流の角度

このとき，プロペラには図 4.3-7 の向きに以下の力が作用する．


推力 T：プロペラ面に垂直方向

トルク Q：プロペラ軸まわり

プロペラ面内に作用する横力 PN：プロペラ面内にあり，プロペラ軸と 90 度

モーメント NP：P-Factor

これらの力を計算する方法には，Ribner の理論[10]，[11] と McCormick の理論[14]，[15] がある．どちらの理論を使ってもほぼ同じ結果が得られるが，McCormick の理論のほうが計算の手間がかからない．
4.3.2.1 Ribner の理論
参考文献 [10]，[11]，[12] によると，プロペラが気流に対して傾いているときのプロペラに作用する力とモーメ
ントの計算式は以下のようになる．傾きP は小さい（20 度以下）と仮定している．
プロペラが気流に正対している場合の推力とトルクは Larrabee の方法で計算してあるとする．
前進係数を次のように表す．
𝐽=

𝑉
𝑛𝐷

(4.3 − 53)

ここで，n は単位時間のプロペラ回転数（１）推力とトルクプロペラが気流に正対しているときの推力係数 CT とパワー係数 CP は次のように定義される．
𝐶 =
𝐶 =

4𝜋 𝑇
𝜌𝛺 𝐷
8𝜋 𝑄
𝜌𝛺 𝐷

=
=

𝑇

(4.3 − 54)

𝜌𝑛 𝐷
2𝜋𝑄

(4.3 − 55)

𝜌𝑛 𝐷

22

ここで， は回転角速度で，𝛺 = 2𝜋𝑛， は空気の密度，T と Q はプロペラが気流に正対しているときの推力とトルク．
プロペラが気流に対して傾いている場合の推力係数とパワー係数を次のように表す．
𝐶 (𝛼 , 𝐽) = 𝜂 𝐶 (0, 𝐽 cos 𝛼 )(4.3 − 56)

𝐶 (𝛼 , 𝐽) = 𝜂 𝐶 (0, 𝐽 cos 𝛼 )(4.3 − 57)

ここで，CT (P, J) はプロペラが気流に対して傾いているときの推力係数，
CT (0, J cos P ) はプロペラが気流に対して正対しているとし，プロペラ軸方向の気流速度がVcos P であるときの推力係数，

T はプロペラの傾きの影響を表す係数，
CP (P, J) はプロペラが気流に対して傾いているときのパワー係数，
CP (0, J cos P ) はプロペラが気流に対して正対しているとし，プロペラ軸方向の気流速度がVcos Pであるときのパワー係数，

P はプロペラの傾きの影響を表す係数
プロペラの傾きの影響を表す係数は次の計算式で表される．
𝜂 = 1+

𝐽 sin 𝛼0.75𝜋
𝐽 cos 𝛼
2 1−
𝐽

𝛿(𝛼 )

(4.3 − 58)𝜂 =1+

𝐽 sin 𝛼0.75𝜋
𝐽 cos 𝛼
2 1−
𝐽

𝛿(𝛼 )

(4.3 − 59)

ここで，J0T は推力がゼロになる推力係数の値，J0P はトルクがゼロになる推力係数の値関数 (P) は次の式で表される．
3
𝛿(𝛼 ) = cos 𝛽 , .
2
ここで，𝛽 , .

1+𝜎.
tan 𝛽 , .

1+ 1+

2 tan 𝛽 , .
𝜎.

(1 − cos 𝛼 )(4.3 − 60)

は 0.75R 位置でのブレード翼素のゼロ揚力線から測った翼角，𝜎.

は 0.75R 位置での Solidity で

𝜎.

=

4𝐵𝑐 .
3𝜋𝐷

(4.3 − 61)（２）横力とモーメント
参考文献 [10]，[11] による計算式を以下に示す．
横力係数 CNP とモーメント係数 CnP は次のように定義される（参考文献 10]）
．この定義は推力係数とトルク係数のときとは形が異なることに注意すること．
𝐶

=

𝐶

=

𝑃
𝜋𝐷
1
𝜌𝑉
2
4
𝑁
𝜋𝐷
1
𝜌𝑉
2
4

=

=

8𝑃

(4.3 − 62)

𝜋𝜌𝑉 𝐷

8𝑁

(4.3 − 63)𝜋𝜌𝑉 𝐷

プロペラの傾きに関する横力係数の微係数，モーメント係数の微係数は，
𝜕𝐶
𝜕𝛼

=

𝑘 𝑓(𝑎 )𝜎 . 𝐼𝐼+𝑘 𝜎 . 𝐼
𝐼 −∆(4.3 − 64)

23

𝜕𝐶
𝜕𝛼

=

𝑘 𝑓(𝑎 )𝑚
1 + 𝑘 𝜎 . (𝐼 − ∆)(4.3 − 65)

ここで，
ks = 1.14 はスピナー係数
ka = 0.4 はサイドウォッシュ係数
ai はインフロー係数で，
1+
𝑎 =

8𝑇
−1
𝜋
=
2

1+

8𝐶
−1
𝜋𝐽
2

(4.3 − 66)

関数 f(ai) は次のように表される．
𝑓(𝑎 ) =

(1 + 𝑎 )[(1 + 𝑎 ) + (1 + 2𝑎 ) ]1 + (1 + 2𝑎 )(4.3 − 67)

 は
∆=

𝜎.

2𝑎 𝐽
𝜎.
𝜋
𝜎 . (1 + 𝜎 .

𝐼 −

𝐼 +

4𝑎 𝐽
𝜋

(4.3 − 68)

𝐼 )I1 は側面積インデックスで
3
𝐼 = 𝐶
4

𝑐
𝑐,

sin 𝛽 𝑑𝜉(4.3 − 69)

ここで，𝜉 = 𝑟⁄𝑅 ，𝜉 はブレードの有効な部分の最小半径比c はブレードの翼弦長，𝐶 はブレードの２次元翼揚力係数
𝛽 はブレード翼素のゼロ揚力線から測った翼角
I2，I3 は
3
𝐼 = 𝐶
4

𝑐
𝑐,

3
𝐼 = 𝐶
4

𝑐
𝑐,

cos 𝛽 𝜉𝑑𝜉(4.3 − 70)(cos 𝜙)𝜉 𝑑𝜉
sin 𝜙(4.3 − 71) は effective helical angle で𝑉(1 + 𝑎 )𝛺𝑟 − 𝑆𝑙𝑖𝑝𝑠𝑡𝑟𝑒𝑎𝑚 𝑟𝑜𝑡𝑎𝑖𝑜𝑛𝑎𝑙 𝑣𝑒𝑙𝑜𝑐𝑖𝑡𝑦
4𝐽𝑎
𝜎. 𝐼 +𝜋𝑚=2(1 + 𝜎 . 𝐼 )𝜙 = tan24(4.3 − 72)(4.3 − 73)

図 4.3-8ブレード翼素と空気流の関係

4.3.2.2 McCormick の理論（１）推力とトルクプロペラが気流に対して傾いているときの推力 T とトルク Q は，次の式の前進係数 J を使って，プロペラ性能
曲線から推力係数 CT とトルク係数 CQ を読み取り計算する（参考文献[13]）．
𝐽=𝑉
𝑉 cos 𝛼=𝑛𝐷
𝑛𝐷

(4.3 − 74)

ここで，VN は気流の速度のプロペラ面に垂直な成分
推力とトルクの計算式は，
𝑇 = 𝜌𝑛 𝐷 𝐶(4.3 − 75)

𝑄 = 𝜌𝑛 𝐷 𝐶(4.3 − 76)（２）横力とモーメント
参考文献[14]の計算式を以下に示す．
パラメータ を次の式で表す．推力 T は前項で計算した値を使う．
𝜏=

1
𝑇
1+ 1+
2
𝑞𝐴

(4.3 − 77)

ここで，A はプロペラのディスク面積で，𝐴 = 𝜋𝑅q は動圧で，𝑞 = 𝜌𝑉 ⁄2プロペラの solidity  は𝜎=𝐵𝑙𝑎𝑑𝑒 𝐴𝑟𝑒𝑎 𝐵 ∫ 𝑐𝑑𝑟=
𝐷𝑖𝑠𝑐 𝐴𝑟𝑒𝑎
𝜋𝑅(4.3 − 78)

ブレードの平均揚力係数 Cl,av と抗力係数 Cd は，
𝐶,

=

3𝑇
𝜋
𝐴𝑞𝜎
𝐽

(4.3 − 79)

25

𝑄

𝐶 =8

−

𝜋
𝐽

𝐴𝑞𝜎𝐷

𝐶, 𝐽
6𝜋

(4.3 − 80)

トルク Q と前進係数 J は前項で計算した値を使う．
横力 PN とモーメント NP は，
𝑃 =

𝜎𝑞𝐴
𝜏 𝐶,
2

𝑁 =

𝜎𝑞𝐴𝑅 2𝜋𝐶 ,
2
3𝐽

+

𝐶 𝐽
𝜋
τ ∙ ln 1 +
2𝜋
𝐽𝜏
+

𝐶
2
1− 𝜏
2
3

𝐽
𝜋

+

𝜋
𝐶
𝐽𝜏

ln 1 +

sin 𝛼𝜋𝐽𝜏(4.3 − 81)
𝜋
− 𝐶
𝐽

sin 𝛼

(4.3 − 82)

ここで，Cl は２次元翼の揚力傾斜で 2
McCormick の理論では，プロペラの性能曲線があれば，横力とモーメントを簡単に計算できる．
4.3.2.3 機体固定座標系から見たプロペラの横力，モーメントの方向
機体軸とプロペラの向きの関係は図 4.2-1 のように表すことができる．
プロペラの向きと気流との角度P は次の式で計算できる．
cos 𝛼 =

1
𝑉

𝑖
𝑈
𝑉 ∙ 𝑖
𝑖
𝑊

=

1
𝑖
𝑉

𝑈+𝑖

𝑉+𝑖

𝑊

(4.3 − 83)

プロペラに作用する横力とモーメントの方向の単位ベクトル j は図 4.3-9 から求める．
まず，速度ベクトル Vc とプロペラ軸の単位ベクトルの作る面に垂直な単位ベクトル k を次の式で求める．
𝐤=

1

𝐕 ×𝐢

𝐕 ×𝐢

𝑉𝑖
𝑊𝑖
𝑈𝑖

1

=
𝑉𝑖

− 𝑊𝑖

+ 𝑊𝑖

− 𝑈𝑖

+ 𝑈𝑖

− 𝑉𝑖

− 𝑊𝑖
− 𝑈𝑖
− 𝑉𝑖

(4.3 − 84)

単位ベクトル j は次の式で計算できる．
𝐣=𝐤×𝐢

𝑘 𝑖
= 𝑘 𝑖
𝑘 𝑖

−𝑘 𝑖
−𝑘 𝑖
−𝑘 𝑖

(4.3 − 85)

横力の方向の単位ベクトル
プロペラ軸方向の単位ベクトル

単位ベクトル

速度ベクトル

⾧さ Vc

図 4.3-9横力の方向

26

4.3.3プロペラの空転

プロペラが空転する場合にプロペラにどういう力が発生するかを考える．図 4.3-10 にプロペラの性能曲線の例を示す．プロペラの空転時にはプロペラとプロペラシャフト間の摩擦力が無いと仮定すると，ゴムが完全に巻き戻
り，空転しだすとゴムからプロペラに伝達されるトルクはゼロになる．図 4.3-10 のプロペラ性能曲線ではトルク係数 CQ がゼロで前進係数 J0P の状態となる．このときの推力係数 CT は負の値となり，プロペラに空気抵抗が発生することがわかる．

0.8

p

0.7

0.6

CT
10×CQ

p

0.5

0.4

10×CQ

0.3

0.2

CT

0.1

0

-0.1
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

前進係数 J
図 4.3-10プロペラの性能曲線の例

27

1.2

1.3

1.4

1.5

1.6

1.7

1.8

4.4 プロペラの回転慣性の影響
プロペラが回転すると，プロペラが角運動量を持つ．プロペラ角運動量は機体の運動に影響をおよぼす．
プロペラの慣性能率 Ip の計算式を以下に示す．

𝐼 =𝐵

𝑟 𝜇𝑑𝑟 = 𝐵

𝑟

+𝑟
2

𝑚

(4.4 − 1)

ここで，B はブレードの枚数，Dp はプロペラの直径，m は単位長さあたりの質量，mi は区間 ri-1～ri の質量

図 4.4-14.4.1

プロペラの慣性能率

プロペラのトルクの釣り合い

プロペラのトルクのフリーボディ・ダイヤグラムを図 4.4-2 に示す．プロペラだけを取り出してモーメントの釣り合いを示したのが図 4.4-3 である．この図は，ゴムとプロペラ軸が作る平面に垂直に見た図である．ゴムのトル
クと軸受けからの支持モーメントの合計とプロペラの空気力によるトルクの差がプロペラの回転角加速度を発生する．釣り合い式は次のように表される．
𝑄
−𝑄
cos 𝛾

=𝐼

𝑑𝛺
2𝜋𝑑𝑛
=𝐼
𝑑𝑡
𝑑𝑡

(4.4 − 2)

ここで，Qr はゴムのトルク，Qprop はプロペラの発生トルク， はゴムとプロペラ軸の成す角， はプロ
ペラの回転角速度ゴムとプロペラ軸の成す角 の計算式を以下に示す（図 4.4-4）．
ゴムの方向を表す単位ベクトルの機体固定座標系の成分が次のようになっているとする．
𝐢

𝑖
= 𝑖
𝑖

(4.4 − 3)

ゴムとプロペラ軸の成す角 は，ゴムの方向ベクトルとプロペラ軸方向ベクトルの内積を使って求めることができる．
cos 𝛾 = 𝐢∙𝐢(4.4 − 4)

機体には，プロペラ軸方向後ろ向きに𝑄 ⁄cos 𝛾 の大きさのモーメント Mr が作用する．

28

𝑖
𝐌 =− 𝑖𝑖𝑄
cos 𝛾(4.4 − 5)

図 4.4-2

図 4.4-3プロペラのトルクの釣り合い

プロペラのモーメントの釣り合い29

図 4.4-4

ゴムとプロペラ軸の成す角30

4.4.2プロペラのジャイロ効果

図 4.4-5 に示すように，プロペラが機体後方から見て右回転しているとする．機首上げしようとすると（ピッチング角速度が正），右旋回しようとするモーメント（ヨーイングモーメントが正）が機体に作用する．このように回転体に角速度を加えると，この回転体にモーメントが発生するのがジャイロ効果である．

図 4.4-5プロペラのジャイロ効果

参考文献[16] により，機体重心に作用するプロペラのジャイロ効果の力の計算式を以下に示す．
機体の運動は，機体重心の速度と機体重心まわりの角速度で表される．速度ベクトルと角速度ベクトルを機体固定座標系の成分で表すと，
𝐕 = 𝑈𝐢 + 𝑉𝐣 + 𝑊𝐤

(4.4 − 6)

𝛚 = 𝑃𝐢 + 𝑄𝐣 + 𝑟𝐤

(4.4 − 7)

参考文献[15] によると，ジャイロ連成によるモーメントは次のように表される（機体固定座標系）
．
𝐌 =

0
sin 𝜃
sin 𝜃
cos 𝜃

− sin 𝜃0− cos 𝜃cos 𝜃

− sin 𝜃
cos 𝜃

cos 𝜃
cos 𝜃
0

機体重心には －MG のモーメントが作用する．

31

𝑃
𝑄 𝐼 𝛺
𝑅

(4.4 − 8)

5

ゴム動力機の３次元運動方程式とその解法

5.1 運動方程式のまとめ
微分方程式をまとめると，𝑚 𝑈̇ + 𝑄𝑊 − 𝑅𝑉 = 𝐹(5.1 − 1)𝑚 𝑉̇ + 𝑅𝑈 − 𝑃𝑊 = 𝐹(5.1 − 2)𝑚 𝑊̇ + 𝑃𝑉 − 𝑄𝑈 = 𝐹(5.1 − 3)

𝑃̇ =

𝑄̇ =
𝑅̇ =

1
𝐼 𝐼

−𝐼

1
[𝐺 − (𝐼
𝐼
1
𝐼 𝐼

−𝐼

−𝐼

+𝐼 𝐼

−𝐼

𝑄𝑅 + 𝐼 (−𝐼+𝐼+ 𝐼 )𝑃𝑄 + 𝐼 𝐺+𝐼 𝐺(5.1 − 4𝑎)− 𝐼 )𝑅𝑃 − 𝐼 (𝑃 − 𝑅 )](5.1 − 5𝑎)

𝐼 (−𝐼

+𝐼

− 𝐼 )𝑄𝑅 + −𝐼 𝐼

+𝐼

+𝐼 𝐺

+𝐼

𝑃𝑄 + 𝐼 𝐺(5.1 − 6𝑎)𝛷̇ = 𝑃 + 𝑄 sin 𝛷 tan 𝛩 + 𝑅 cos 𝛷 tan 𝛩(5.1 − 7𝑎)𝛩̇ = 𝑄 cos 𝛷 − 𝑅 sin 𝛷(5.1 − 8𝑎)𝛹̇ = 𝑄 sin 𝛷 sec 𝛩 + 𝑅 cos 𝛷 sec 𝛩(5.1 − 9𝑎)𝑑𝑥= 𝑈 cos 𝛩 cos 𝛹 + 𝑉(sin 𝛷 sin 𝛩 cos 𝛹 − cos 𝛷 sin 𝛹)𝑑𝑡+ 𝑊(cos 𝛷 sin 𝛩 cos 𝛹 + sin 𝛷 sin 𝛹)(5.1 − 12)𝑑𝑦= 𝑈 cos 𝛩 sin 𝛹 + 𝑉(sin 𝛷 sin 𝛩 sin 𝛹 − cos 𝛷 cos 𝛹)𝑑𝑡+ 𝑊(cos 𝛷 sin 𝛩 sin 𝛹 − sin 𝛷 cos 𝛹)𝑑𝑧= −𝑈 sin 𝛩 + 𝑉 sin 𝛷 cos 𝛩 + 𝑊 cos 𝛷 cos 𝛩𝑑𝑡𝑛=(5.1 − 14)

𝑑𝑚
𝑑𝑡

𝑄
−𝑄
cos 𝛾

(5.1 − 13)(3.3 − 2)
=𝐼

2𝜋𝑑𝑛
𝑑𝑡

変数が，𝑈, 𝑉, 𝑊, 𝑃, 𝑄, 𝑅, , , , 𝑥 , 𝑦 , 𝑧 , 𝑚(4.4 − 2), 𝑛 の 14 個で，微分方程式が上の 14 個で，この微分方程式は解くことができる．

5.2 運動方程式の解法
MATLAB 互換の高機能数値計算ツール Octave を使って微分方程式を解く（参考文献[21]，[22]）．
解析結果は変数のグラフ，飛行軌跡/機体姿勢の図で示す．

32

6

ゴム動力機の３次元運動解析の例

6.1 解析対象
解析対象の機体の写真を図 6.1-1 に示す．図面を図 6.1-2 に示す．機体諸元を表 6.1-1 に示す．主翼面積，平均空力翼弦の計算を表 6.1-2 に示す．水平尾翼容積，垂直尾翼容積の計算を表 6.1-3 に示す．水平尾強容積と重心位置の関係を図 6.1-3 に示す．

図 6.1-1

解析対象の機体33

図 6.1-2

解析対象の機体の図面

34

表 6.1-1

解析対象の機体の諸元

機体諸元
重量

プロペラ

空虚重量（実測）

11.75 gf

ゴム重量（実測）

2.82 gf

全備重量（実測）

14.57 gf

全備重量（運動解析）

15 gf

直径

170 mm

ピッチ比

1.4

ブレード形状

一定翼弦⾧

重心位置
主翼

400 mm

面積29100 mm^2

アスペクト比5.50

平均空力翼弦

73.3 mm
6%キャンバー折り曲げ翼

翼幅160 mm

面積8400 mm^2

アスペクト比3.05

平均空力翼弦

53.5 mm

翼型

垂直尾翼

Peckpolymer社プロペラと同じ

30 %MAC
翼幅

翼型
水平尾翼

#18輪ゴム15本＋S環

平板翼

アーム⾧

162.9 mm

水平尾翼容積

0.641

翼⾧

55 mm

面積

2750 mm^2

平均空力翼弦翼型

51.5 mm
平板翼

アーム⾧168.1 mm

垂直尾翼容積0.039735

表 6.1-2

主翼面積，平均空力翼弦の計算

主翼
0

翼幅 b

-10

400 mm

C.G.

-20
-30

mm

S

-40

x_bar

MAC

-50

LE 95.32901

29100 mm^2

-60

TE 95.32901 -73.34085052 -73.3409

y'
0

0

-70

AR

5.498282

MAC

73.34085 mm

x_bar

95.32901 mm

e_bar

0 mm

断面 i

y

-80
-90
0

BL

前縁

後縁

x

x

y_LE

y_TE

(mm)

(mm)

(mm)

(mm)

c

50

100

150

200

250

delta_xc_ave

2*delta_Sc^2c^2_avedelta_c_bar

y_ave

delta_y_bar

e

6084

e_ave

delta_e_bar

y'_LE

0

y'_TE

1

0

0

0

-78

78

0

-78

2

23.75

23.75

0

-78

78

23.75

78

3705

6084

6084

9.930927835

11.8751.511920103

0

0

0

0

-78

3

47.5

47.5

0

-78

78

23.75

78

3705

6084

6084

9.930927835

35.6254.535760309

0

0

0

0

-78

4

71.25

71.25

0

-78

78

23.75

78

3705

6084

6084

9.930927835

59.3757.559600515

0

0

0

0

-78

5

95

95

0

-78

78

23.75

78

3705

6084

6084

9.930927835

83.125

10.58344072

0

0

0

0

-78

108.125 108.125

0

-75.5

75.5

13.12576.75

2014.688 5700.25 5892.125 5.315061211

101.5625 7.031501691

0

0

0

0

-75.5

7

121.25121.250-73

73

13.12574.25

1949.063

5514.625 4.974532861

114.6875 7.681550016

0

0

0

0

-73

8

134.375 134.375

0

-70.5

70.5

13.12571.75

1883.438 4970.25 5149.625 4.645280284

127.8125 8.272400532

0

0

0

0

-70.5

6

5329

9

147.5

147.5

0

-68

68

13.125

69.25

1817.813

4797.125 4.327303479

140.9375 8.804053238

0

0

0

0

-68

10

160.625 160.625

0

-65.5

65.5

13.12566.75

1752.188 4290.25 4457.125 4.020602448

154.0625 9.276508135

0

0

0

0

-65.5

11

173.75173.750-63

63

13.12564.25

1686.563

167.1875 9.689765222

0

0

0

0

-63

12

186.875 186.875

0

-60.5

60.5

13.12561.75

1620.938 3660.25 3814.625 3.441027706

180.312510.0438245

0

0

0

0

-60.5

0

-58

58

13.12559.25

1555.313

193.4375 10.33868597

0

0

0

0

-58

13

200

200

4624
3969
3364

4129.625 3.725177191
3512.125 3.168153995

36表 6.1-3

機体名

水平尾翼容積，垂直尾翼容積の計算

 一般の翼の場合
30 %MAC

重心位置

：翼面積

主翼
翼幅

b

400 mm

半翼幅b/2

200 mm

翼弦長，翼根c_r

78 mm

翼弦長，翼端c_t

58 mm

翼端の後退長e_t

テーパー比

lambda

翼面積

S

アスペクト比AR

平均空力翼弦c_bar

73.3 mm

平均空力翼弦の位置

y_bar

95.3 mm

e_bar

0 mm

0 mm
0.74359

1/4 c

 テーパー翼の場合：テーパー比29100 mm^25.498282

1/2 c：後退角 L

水平尾翼
翼幅b

160 mm

半翼幅

b/2

80 mm

翼弦長，翼根c_r

65 mm

翼弦長，翼端c_t

40 mm

翼端の後退長e_t

テーパー比

lambda

翼面積S_h

アスペクト比AR3.047619

平均空力翼弦c_bar

53.49206 mm

平均空力翼弦の位置

y_bar

36.8254 mm

e_bar

11.50794 mm

LAMBDA

0.23022 radian

後退角

25 mm
0.615385
8400 mm^2

 水平尾翼容積：主翼の平均空力翼弦：主翼面積

 垂直尾翼容積：主翼の翼弦⾧

LAMBDA 13.19061 degree
垂直尾翼
翼長

h

55 mm

翼弦長，翼根c_r

65 mm

翼弦長，翼端c_t

35 mm

翼端の後退長e_t

テーパー比

lambda

翼面積S_v

平均空力翼弦c_bar

51.5 mm

平均空力翼弦の位置

y_bar

24.75 mme_bar

13.5 mm

後退角：平均空力翼弦30 mm
0.538462
2750 mm^2

LAMBDA 0.388319 radian
LAMBDA 22.24902 degree

主翼前縁と水平尾翼前縁の距離

dist_w-h

主翼前縁と垂直尾翼前縁の距離

dist_w-v

160 mm水平尾翼モーメントアーム長

l_h162.891 mm

垂直尾翼モーメントアーム長

l_v

168.05 mm

水平尾翼容積V_h

0.641

垂直尾翼容積

V_v

0.039703

160 mm37

図 6.1-3

水平尾翼容積と重心位置の関係

38

6.2 機体の空力特性
機体の空力特性を計算するのに使用した AVL モデルの座標系を図 6.2-1 に示す．
AVL モデルを図 6.2-2 に示す．
モデルの座標系の原点を胴体の棒の先端上面にとり，空力係数の基準点を重心位置（xcg = 92，ycg = -1.7，zcg = 0.0）としている．主翼の取付角は 2 度，水平尾翼の取付角は定常滑空時の全機 CL が 0.8 となる角度とした．有害抵抗係数 CDP = 0.06 とした．

図 6.2-1

AVL モデルの座標系

39図 6.2-2

AVL モデルAVL 解析結果のまとめを表 6.2-1 に示す．定常滑空時の迎角は 5.75 度で，全機 CL は 0.8，全機 CD は 0.0977 である．全機 AVL モデルを使って主翼と水平尾翼の揚力係数も計算した．
運動解析に使用する空力微係数のまとめを表 6.2-2 に示す．空力微係数のうち，CL，Cm，Clq，Cmq，Cy，Cl，
Clp，Cyr，Cnr は迎角の影響はほとんどないので，定常滑空時の値を使う．CD は CL の２次関数，Cnr は迎角の 2 次の関数，CL，Cyp，Cnp，Clr は迎角の 1 次の関数であるとして近似式を導出した．
全機失速迎角をstall = 9.43 度，そのときの最大揚力係数を CLmax = 1.10 とした．失速後の全機揚力係数と全機抗力係数は Viterna and Corrigan の方法（参考文献[8]）を使って図 6.2-4 のように推定した．

40

機体軸 変更= +表 6.2-1

全機 AVL 解析結果

CDtot

Cma

0 deg

alpha

alpha

alpha

(degree)(radian) (degree)

Cmtot

CLtot

CLa

CYb

Clb

Cnb

CYp

Clp

Cnp

CLq

Cmq

CYr

Clr

Cnr

-3.94221-0.0688-3.94

0.1596

00.06454.72 -0.780 -0.431 -0.175 0.0877 -0.370-0.471 -0.0087

8.46

-10.1

0.275

0.060

-0.0976-1.52283-0.02658-1.52

0.12490.2

0.06564.75 -0.867 -0.442 -0.192 0.0868 -0.321-0.470 -0.0242

8.38

-10.1

0.318

0.103

-0.1030

0.89152

0.01556

0.89

0.08660.4

0.07154.74 -0.948 -0.451 -0.210 0.0874 -0.269-0.467 -0.0400

8.29

-10.1

0.358

0.146

-0.1091

3.31250.057814

3.31

0.04500.6

0.08214.72 -1.022 -0.460 -0.227 0.0895 -0.214-0.462 -0.0563

8.19

-10.1

0.393

0.188

-0.11614.52920.079049

4.53

0.02290.7

0.08934.70 -1.056 -0.464 -0.236 0.0911 -0.185-0.458 -0.0646

8.13

-10.1

0.409

0.208

-0.1198

5.75210.100393

5.75

0.00000.8

0.0977

4.67 -1.089 -0.468 -0.245 0.0931 -0.156-0.455 -0.0732

8.07

-10.1

0.423

0.229

-0.12376.98283

0.121873

6.98-0.02370.9

0.10744.64 -1.120 -0.472 -0.253 0.0955 -0.125-0.451 -0.0820

8.00

-10.1

0.437

0.249

-0.1278

8.223120.143521

8.22-0.048310.11834.60 -1.150 -0.476 -0.262 0.0983 -0.094-0.446 -0.0910

7.93

-10.1

0.449

0.270

-0.1321alpha

CL_wing

CL_hs*St/Sw

(deg)

(rad)

CLa_hs*St/Sw =0.008658

0.49608

alpha_hs0 =8.727995

0.15233

CL_hs/Sw

誤差

誤差^2(degree)
-3.94

0.109

-0.109-0.10970

0.00

0.00000

-1.52

0.289

-0.089-0.08875

0.00

0.00000

0.89

0.468

-0.068-0.06785

0.00

0.00000

3.31

0.647

-0.047-0.04689

0.00

0.00000

4.53

0.736

-0.036-0.03635

0.00

0.00000

5.75

0.826

-0.026-0.02577

0.00

0.00000

6.98

0.915

-0.015-0.01511

0.00

0.00000

8.22

1.004

-0.004-0.00437

0.00

0.00000

合計

0.00000

表 6.2-2

運動解析に使用する空力微係数

rad

deg

alpha_0-0.0709-4.0600

失速迎角alpha_stall

0.16469.432

1.1

定常滑空迎角alpha_eq

0.10045.752

0.8

失速時抵抗係数

CD_stall

0.13023

ゼロ揚力迎角CYb

CL_total

CLa

Cma

Clb

Cnb

CYp

Clp

Cnp

CLq

Cmq

CYr

Clr

Cnr

4.671

-1.089 -0.468 -0.245

0.093

１次式

-0.455１次式

8.067-10.101

0.423

１次式

２次式

近似式

a

1.2888-0.3858

0.9928-0.2680の係数b-0.2829-0.0349

0.1286-0.1225

c

--

--

--

-0.1048

１次式：𝐶= 𝑎𝛼 + 𝑏

２次式：𝐶= 𝑎𝛼 + 𝑏𝛼 + 𝑐

ここで，迎角の単位は rad である．

41

1.1
1
0.9

揚力係数
CL

0.8

主翼
0.7
CL,hs*St/Sw 近似式

0.6

CL,wing AVL結果
0.5

CL,hs*St/Sw AVL結果
全機CL AVL結果

全機

0.4
0.3
0.2

水平尾翼

0.1
0
-4

-3

-2

-1

0

1

2

3

-0.1

4

5

6

7

8

9

機体迎角 （度）-0.2

図 6.2-3

主翼と水平尾翼の揚力係数

全機揚力係数と全機抗力係数は次の式で近似する．係数は図 6.2-4 を参照のこと．ここで，迎角の単位は rad で
ある．




失速前
𝐶 = 𝐶 (𝛼 − 𝛼 )(6.2 − 1)

𝐶 =𝐶(6.2 − 2)

+𝑘

𝐶

失速後
(cos 𝛼)sin 𝛼(6.2 − 3)

𝐶 = 𝐵 (sin 𝛼) + 𝐵 cos 𝛼(6.2 − 4)

𝐶 = 𝐴 sin 2𝛼 + 𝐴

𝐴 =𝐶(6.2 − 5)

2

𝐴 = 𝐶

−𝐶

sin 𝛼

cos 𝛼

sin 𝛼
(cos 𝛼(6.2 − 6)

)

𝐵 =𝐶
𝐵 =

𝐶

(6.2 − 7)(sin 𝛼−𝐶cos 𝛼)(6.2 − 8)42Viterna and Corrigan の方法
CLa(1/rad)

4.67

a0

(deg)

-4.06023(rad)-0.07086

CLstall

1.1

a_stall(deg)

9.431722(rad)0.164615

CD0

0.06

CDstall

0.130226

CDmax

1.5

k_CD

0.058038

A1

0.75

A2

0.144401

B1

1.5

B2

0.091178

1.6
CL, Viterna and Corrigan の方法

1.5

全機 CL，CD

1.4

CD, Viterna and Corrigan の方法

1.3

CL, AVL

1.2

CD, AVL

1.1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
-20

-10

0

10

20

30

40

50

60

機体迎角 （度）図 6.2-4

失速後の全機揚力係数と全機抗力係数

43

70

80

90

100

6.3 慣性データ
慣性データの基準座標系は AVL モデルの座標系（図 6.2-1）と同じである．図面に基づいて機体の慣性データを計算した結果を表 6.3-1 に示す．実際の模型の全備重量はゴム（2.6 gf）を含めて 14.3 gf であった．実機のほうが11%重いが，重心位置は計算とほぼ一致していた．運動解析には回転半径の計算値を使い，全機質量を 15 g とする．x 方向の重心位置は 30%MAC（x = 92 mm）を使い，z 方向の重心位置は計算値の -1.7 mm を使う．

表 6.3-1
慣性能率，慣性乗積の計算

機体の慣性データの計算結果

解析用ゴム動力機，スパン400mm

作成：滝

2022/9/26

機体重心位置
X_cg

Y_cg

Z_c

(mm)

(mm)

(mm)

95.3
質量

0.0

-1.7 <== 下の計算値を入力

重心位置 - 機体座標系

機体重心まわりの慣性 - 機体軸方向

M

X_cg

Y_cg

Z_cg

(g)

(mm)

(mm)

(mm)

IX_cg

IY_cg

IZ_cg

IXY_cg

IYZ_cg(g-mm^2) (g-mm^2) (g-mm^2) (g-mm^2) (g-mm^2)

IZX_cg

M*X_cg

M*Y_cg

M*Z_cg(g-mm^2)(g-mm)(g-mm)(g-mm)

主翼 LH+RH

1.5

103.8

0.0

25.3

34535

1937

34058

0

0

322

158

0

39

主翼中央部

1.9

109.0

0.0

12.8

8448

1694

9351

0

0

371

204

0

24

水平尾翼

0.8

268.4

0.0

1.0

835

23609

24433

0

0

364

210

0

1

垂直尾翼

0.3

269.4

0.0

27.4

280

8098

7818

0

0

1310

69

0

7

胴体

8.3

69.1

0.0

-11.0

1013

64372

63365

0

0

4294

575

0

-92

45112

99711

139026

0

0

6661

1215

0

-22

合計
質量

合計

重心位置 - 機体座標系

機体重心まわりの慣性 - 機体軸方向

M

X_cg

Y_cg

Z_cg

(g)

(mm)

(mm)

(mm)

12.8

95.3

0.0

-1.7

IX_cg

IY_cg

IZ_cg

IXY_cg

回転半径

IYZ_cg(g-mm^2) (g-mm^2) (g-mm^2) (g-mm^2) (g-mm^2)
45112

99711

44

13902600

IZX_cgrX

rY

rZ(g-mm^2)

(mm)

(mm)

(mm)

6661

59.5

88.4

104.4

rZX
22.9

6.4 動力ゴムの特性#18 輪ゴムの 8 条と FAI TAN Super Sport 3/32 インチ幅 4 条の２種類を使用する．図 6.4-1 に巻き数-トルク関係の測定結果を示す．このグラフは巻き戻し時のトルクを表している．輪ゴムは手巻きの場合を，FAI ゴムは伸ばしてワインダーで巻いた場合のデータである．
0.025

実線：#18輪ゴム，8条，手巻き
破線：FAI TAN Super Sport，3/32”幅，4条
0.02

トルク
(N-m)

FAI TAN Super Sport，3/32”幅，4条
実用最大巻き数：2.3回/mm#18輪ゴム
破断巻き数：1.69回/mm#18輪ゴム，8条，手巻き実用最大巻き数：1.3回/mm

0.015

0.01

0.005

0
0

0.5

1

1.5

2

2.5

3

単位⾧さあたり巻き数 （回/mm）
図 6.4-1

使用する動力ゴムの巻き数-トルク関係

これらのデータを使って解析用にモデル化した巻き数-トルク関係式を以下に示す．
動力ゴムの巻き数-トルクの関係を表 6.4-1 に示すように数式で表現する．数式と必要なパラメータは以下のとおり．
𝑁 ≤ 𝑁 − ∆𝑁 のとき，𝑄 = 𝑄
𝑁 > 𝑁 − ∆𝑁 のとき，𝑄 = 𝑄+ ∆𝑄 = 𝑄+ 𝑏[𝑁 − (𝑁 − ∆𝑁 )]ここで，Qr：ゴムのトルク，N：巻き戻し時の単位長さあたりの巻き数，
NS：巻き完了時の単位長さあたりの巻き数，
𝑄= 𝑓 (𝑁 )：巻き完了時のトルク，4 次までの多項式で表す．

𝑓 (𝑁 ) = 𝑎 𝑁 + 𝑎 𝑁
𝑄

+𝑎 𝑁

+𝑎 𝑁

= 𝑓 (𝑁)：ヒステリシスの無いトルクカーブ，0.5 次の項を含む多項式で表す．

𝑓 (𝑁) = 𝑐 . 𝑁 . + 𝑐 𝑁 + 𝑐 𝑁∆𝑇 = 𝑏[𝑁 − (𝑁 − ∆𝑁 )] ：トルクのヒステリシス分𝑏=𝑄− 𝑓 (𝑁 )∆𝑁=𝑓 (𝑁 ) − 𝑓 (𝑁 )∆𝑁

NB：ヒステリシスがある巻き数範囲
機体のゴムの長さで実用最大巻き数にした場合の，このモデルを使って求めた巻き数-トルク関係を図 6.4-2 に示す．

45

表 6.4-1

巻き数-トルク関係式

ゴム種類
幅

FAI TAN Super Sport#18輪ゴム(inch)

0.09375

(mm)

2.38

(inch)

0.042

(mm)

1.07

1.1

条数

本

4

8

断面積

(mm^2)

10.16

9.68

質量

(g/mm)

0.00944

0.01

厚さ

巻き方
実用最大巻き数

ワインダー

-1.1
--

手巻き

(回/mm)

2.3

1.3

a1

10.71

34.62

a2

-6.226-70.789

a3

1.608

60.897

トルクカーブのa4

0.0446-16.056

パラメータ

c0.5

4.487

6.470

c1

-0.957-2.649c2-0.07080.24920.8

0.4

delta NB

注１：使用した単位は，N，mm である．
注２：#18 輪ゴム１本の重量は 0.175 gf，1 本を切断した長さは 140 mm である．
注３：FAI ゴムの比重は 1.0 である．

46

15

巻き数：340回
巻き数：640回

トルク(N-mm)#18輪ゴム，8条，手巻き
ループ⾧ 263 mm

105

FAI TAN Super Sport，3/32”幅，4条
ループ⾧ 280 mm

0
0

100

200

300

400

巻き数 （回）図 6.4-2

巻き数とトルクの関係

47

500

600

700

6.5 プロペラ特性
6.5.1

プロペラ仕様プロペラを図 6.5-1 に示す．円筒巻き付け法で製作した．プロペラの仕様は以下のとおり．


プロペラ仕様

ピッチ比 P/D：1.4 一定



翼弦長分布：Peck Polymer 社のプロペラの翼弦長分布を採用（図 6.5-2）

翼型：薄肉 5%キャンバー翼，Reynolds 数 20,000 の空力データ（図 6.5-3，参考文献[23]）を使ってViterna and Corrigan の方法で失速後の特性を推定（表 6.5-1）．

Solidity：0.124（計算を表 6.5-2 に示す）

図 6.5-1

ピッチ比 1.4 のプロペラ

48

0.3

0.25

翼弦⾧分布
c/Rp

0.2

0.15

0.1

0.05

0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

半径 r/Rp
図 6.5-2

図 6.5-3翼弦長分布

２次元翼の揚力係数と抗力係数データ – 薄肉円弧翼（キャンバー5%），レイノルズ数 20,000

49

[23]

表 6.5-1

プロペラ性能計算に使用する２次元翼の空力データ

Viterna and Corrigan の方法

cl_a

2.5

5.79

a0

(deg)

-2.4

(rad)

-0.04189(deg)

8.287296(rad)0.144641

cl_stall2

1.08

a_stall
cd_0

0.038

cd_stall

0.078

AR

6

cd_max

2

k

1.5
1
0.5
0
-20

0.034294

A1

1

A2

0.116981

B1

2

B2

0.036834

0

20

40

60

80

100

120

-0.5

表 6.5-2プロペラの Solidity の計算

プロペラのSolidity -- Peckpolymer社製プロペラ展開形
0.3

プロペラ直径

D

(mm)

100

(m)

0.1

プロペラ半径

R

(mm)

50

(m)

0.05

ブレード枚数B

Solidity

sigma

0.25

𝑅

𝜎=

𝐵𝑙𝑎𝑑𝑒 𝐴𝑟𝑒𝑎 𝐵 ∫0 𝑐𝑑𝑟=
𝐷𝑖𝑠𝑐 𝐴𝑟𝑒𝑎
𝜋𝑅2

0.2

2

0.15

0.12409

合計 =0.0004873
0.1

xi

r

c

c/R

(m)

(m)

xi

r

delta_rcc*dleta_r

(m)

(m)

(m)

0.05
0

0.125

0.00625

0.00250.05

0.15

0.00750.0025 0.003675

9.188E-06

0.175

0.00875 0.00485

0.097

0.2

0.01

0.0025 0.005925

1.481E-05

0.225

0.01125

0.007

0.14

0.25

0.01250.00250.008

0.000020.275

0.01375

0.009

0.18

0.3

0.015

0.00250.00965

2.413E-05

0.325

0.016250.0103

0.206

0.35

0.01750.0025 0.010825

2.706E-05

0.375

0.01875 0.01135

0.227

0.4

0.02

0.0025 0.011725

2.931E-05

0.425

0.02125

0.0121

0.242

0.45

0.02250.0025 0.012325

3.081E-05

0.475

0.02375 0.01255

0.251

0.5

0.025

0.00250.01265

3.163E-05

0.525

0.02625 0.01275

0.255

0.55

0.02750.00250.0128

3.2E-05

0.575

0.02875 0.01285

0.257

0.6

0.03

0.00250.01293.225E-05

0.625

0.03125 0.01295

0.259

0.65

0.03250.0025 0.012925

3.231E-05

0.675

0.03375

0.0129

0.258

0.7

0.035

0.0025 0.012875

3.219E-05

0.725

0.03625 0.01285

0.257

0.75

0.03750.0025 0.012775

3.194E-05

0.775

0.03875

0.0127

0.254

0.8

0.04

0.825

0.04125 0.01255

0.251

0.85

0.04250.875

0.04375 0.01245

0.249

0.9

0.925

0.04625 0.01235

0.247

0.0025 0.012625

3.156E-05

0.00250.01253.125E-05

0.045

0.00250.01240.000031

0.9375 0.046875

0.00125

0.0123

1.537E-05

0.95

0.0475 0.01225

0.245

0.9625 0.048125

0.00125 0.012225

1.528E-05

0.975

0.04875

0.01220.244

0.9875 0.049375

0.00125

1.519E-05

1

0.05

0.0121

0.242

0.01215

50

0

0.2

0.4

0.6

0.8

1

6.5.2プロペラ性能曲線

Larrabee の方法を使って計算したプロペラの性能曲線を図 6.5-4 に示す．このグラフを人工ニューラルネットワークの近似式で表したものが表 6.5-3，図 6.5-5 である．トルクがゼロとなる前進係数 J0P とその時の推力係数 CTの値は，𝐽= 1.586, 𝐶 = −0.020 である．

0.8

p

0.7

0.6

CT
10×CQ

p

0.5

0.4

10×CQ

0.3

0.2

CT

0.1

0

-0.1
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

前進係数 J
図 6.5-4

表 6.5-3プロペラ性能曲線

プロペラ性能曲線の近似式のパラメータ

図 6.5-5 に示すニューラルネットワークの正規化係数，重み，バイアスの値を示す．

正規化係数
J

CT
1

重み

CQ
0.2

バイアス

0.130008 0.016619
0.577394

0.01457-0.00011 -1.00769

0.05

重み

バイアス

重み

バイアス

-3.20698 0.726817 1.271703 0.055661 1.115162

7.210628 -3.76679 4.284665 0.877135-3.05972 6.920664 -0.11869 0.627164 -3.49764

7.239478 0.144716 -0.48386 -0.05601

0.528821 0.348424 0.751921

1.22375 -0.40249

3.250049 -2.62554

51

入力

入力層
（4ニューロン）
tanh 関数隠れ層（4ニューロン）
tanh 関数

S

S
tanh関数S

tanh関数S

S

線形関数

tanh関数S

出力

出力層（2ニューロン）
線形関数

tanh関数S

S

線形関数

tanh関数tanh関数S
tanh関数

入力

和（合計）
を表す

多入力ニューロン

伝達関数
出力

S

バイアス
バイアスに関する入力（常に１）𝑛 = 𝑤 , 𝑝 + 𝑤 , 𝑝 + 𝑤 , 𝑝 + ⋯+ 𝑤 , 𝑝 + 𝑏
重み係数

図 6.5-5プロペラ性能曲線の近似式に使ったニューラルネットワーク

52

6.5.3プロペラに見流が斜めからあたる影響

横力の係数とモーメントの係数を Ribner の理論で計算した場合と，McCormic の理論で計算した場合を比較したのが図 6.5-6 である．実際に飛行する領域（J > 0.5）では両者はほぼ等しい．運動解析では McCormic の理論を使う．

0.5
McCormicの理論

0.45

Ribnerの理論
0.4

0.35

0.3

0.25

0.2

0.15

0.1

0.05

0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

前進係数 J
図 6.5-6

Ribner の理論と McCormic の理論の比較

53

1.2

1.3

1.4

1.5

1.6

1.7

6.5.4プロペラ慣性能率

プロペラ直径 170 mm（図 6.5-1）の場合のプロペラ軸まわりの慣性能率の計算結果を表 6.5-4 に示す．プロペラブレードの材料は 1 mm 厚のプリスチレン板である．質量 3.40g，回転半径 50.0 mm となった．プロペラ直径が変わる場合は，質量と回転半径がプロペラ直径に比例すると仮定する．

表 6.5-4

プロペラ慣性能率の計算 – プロペラ直径 170 mm の場合

慣性能率 -- プロペラ
Dp=

170 mmブレード比重 =

半径

半径比

翼弦⾧比

r

r/Rp

c/Rp

半径
c/Rp

(mm)

1

翼弦⾧

厚さ

体積

質量

r

delta_xctvolumem

慣性モーメント
m*r^2

(mm)

(mm)

(mm)

(mm)

(mm^3)(g)(g-mm^2)

10

0.118

0.045

20

0.235

0.148

0.0965

15

10

8.2

1

82.025 0.082025 18.45563

30

0.353

0.219

0.1835

25

10

15.6

1

155.975 0.155975 97.48438

40

0.471

0.25

0.2345

35

10

19.9

1

199.325 0.199325 244.1731

50

0.588

0.257

0.2535

45

10

21.5

1

215.475 0.215475 436.3369

60

0.706

0.258

0.2575

55

10

21.9

1

218.875 0.218875 662.0969

70

0.824

0.251

0.2545

65

10

21.6

1

216.325 0.216325 913.9731

80

0.941

0.245

0.248

75

10

21.1

1

85

1.000

0.242

0.2435

82.5

5

20.7

1

210.8

0.21081185.75

103.4875 0.103488 704.3618

合計（１ブレード） 1.402288 4262.632
２ブレード 2.804575 8525.264
ハブ質量
合計

0.6

回転半径

54

0

3.40 8525.264
50.0 mm

6.6 ３次元運動解析結果
6.6.1

機体固定座標系

３次元運動解析では機体固定座標系を図 6.6-1 のようにとる．機体重心が機体固定座標系の原点であり，X 軸を胴体軸と平行にとる．ゴムの方向は X 軸に平行であるとする．プロペラ中心の座標は，(118, 0, 22) である．プロ
ペラ中心と重心位置を結んだ線は X 軸に対して 10.56 度下を向いている．

図 6.6-1

機体固定座標系

55

6.6.2解析ツール

３次元運動方程式を解く Octave のスクリプトを付録に示す．
6.6.3

解析ケースと解析結果

解析ケースを表 6.6-1～表 6.6-5 に示す．#18 輪ゴム×8 条を動力ゴムとして使用し，340 回巻くとした．プロペ
ラはゴムの断面積に適合するように直径 170 mm とした．このゴム断面積は急上昇をねらったものである．左旋回
直径 12 m の滑空をするように，ラダー舵角を決め，きれいに上昇するようにダウンスラスト角 3 度，右サイドス
ラスト角 0.5 度としたものをベースライン（ケース 1）とした．発進条件はすべてのケースで，1.5 m の高さから，
速度 3 m/s で発進するとした．ベースラインでは 27.5 m まで上昇する．
解析結果を図 6.6-2～図 6.6-21 に示す．解析結果からわかることは以下のとおり．


プロペラの効果
ケース 1．ケース 1a，ケース 1b，ケース 1c から，プロペラがゴム動力機の運動に及ぼす影響がわかる（表6.6-1，図 6.6-3～図 6.6-5）．
最も大きいのは，プロペラに気流が斜めからあたることによる横力とモーメントの影響である（パワー効果，
ケース 1a）．この機体の場合，プロペラの軸が機体軸よりかなり下を向いているので，動力飛行では気流がプロペラの上のほうからあたる．したがって，図 6.6-22 に示すように，プロペラに下向きの力と右向きのヨー
イングモーメントが作用する．このヨーイングモーメントは，ゴムによる左翼下げのローリングモーメントによる左旋回の傾向を打ち消す効果がある．この効果が無いと，ケース 1a のように左旋回がきつくなり，機体が上昇しない．機体が上昇しないので，ゴムのエネルギが主に速度増加に使われる．
次に影響が大きいのは，プロペラのジャイロ効果である（ケース 1b）．左旋回するので，機首を左に向ける
ヨーイング角速度がある．このヨーイング加速度により，プロペラのジャイロ効果で頭上げのモーメントが発生する．ジャイロ効果がないと頭上げが生じないので，図 6.6-4 のように上昇がゆるくなる．
プロペラの回転慣性の影響は小さい（ケース 1c，図 6.6-5）．



ダウンスラストの影響
プロペラ中心と機体重心を結んだ線を基準としたプロペラ軸のダウンスラスト角の影響をケース 2～ケー
ス 4（表 6.6-2，図 6.6-6～図 6.6-8）に示す．ダウンスラスト 1 度～5 度では最大高度はほとんど変化しない
が，アップスラストにすると，最大高度は急に減少する．



ラダー舵角とサイドスラストの影響
ラダー舵角とサイドスラストの組み合わせを変えて計算した結果をケース 5～ケース 8（表 6.6-3，図 6.6-9～図 6.6-12）に示す．ラダー舵角（滑空の旋回の向きが決まる）とサイドスラストをうまく組み合わせることによって，高度を大きくすることができる．左旋回滑空も右旋回滑空も可能である．



発進角度の影響
発進時の上向き角度によって，飛行パターンが変わる（表 6.6-4，図 6.6-13～図 6.6-16）．35 度～55 度の間が適切である．角度が小さすぎると速度が大きいままで，上昇しない（ケース 9，図 6.6-13）．



FAI TAN Super Sport ゴムの場合
FAI TAN Super Sport ゴムの 3/32 インチ幅×4 条を使った場合の計算結果を表 6.6-5 と図 6.6-17～図 6.6-21に示す．機体の舵角，プロペラ軸の角度の設定は輪ゴムの場合と同じでよい．発進角度は輪ゴムの場合より大きくて，50 度～70 度が適正である．バースト領域のゴムのエネルギが大きいためである．最大高度は輪ゴムの場合の 2 倍弱の 49 m に達する．

56

表 6.6-1プロペラの効果ケース１（ベースライン）

ケース１a

ケース１b

ケース１c
〇

プロペラに斜めから気流があたる影響

〇

×

〇

プロペラのジャイロ効果

〇

〇

×

〇

プロペラの慣性の影響

〇

〇

〇

×

全機質量

(g)

15

15

15

15

プロペラ 直径

(mm)

170

170

170

170
1.4

1.4

1.4

1.4

プロペラ中心-重心からのダウンスラスト

プロペラ

ピッチ比
(deg)

3

3

3

3

サイドスラスト

(deg)

0.5

0.5

0.5

0.5
左旋回 直径12m

ラダー舵角

左旋回 直径12m

左旋回 直径12m

左旋回 直径12m

エルロン舵角

舵角無し

舵角無し

舵角無し

舵角無し

使用ゴム

#18輪ゴム×15本#18輪ゴム×15本#18輪ゴム×15本#18輪ゴム×15本

8条

8条

8条

8条

巻き数

340

340

340

340
0

巻き戻し数

0

0

0

発進速度

(m/s)

3

3

3

3

発進角度

(deg)

45

45

45

45

最大高度

(m)

27.5

15.4

23.7

27.8

表 6.6-2

ダウンスラスト角の影響
ケース１

ケース２

ケース3

プロペラに斜めから気流があたる影響

〇

〇

〇

〇

プロペラのジャイロ効果

〇

〇

〇

〇

プロペラの慣性の影響

〇

〇

〇

〇

(g)

15

15

15

15

(mm)

170

170

170

170

1.4

1.4

1.4

1.4

全機質量
プロペラ
プロペラ

直径

ピッチ比

（ベースライン）ケース4

プロペラ中心-重心からのダウンスラスト

(deg)

-1

1

3

5

サイドスラスト

(deg)

0.5

0.5

0.5

0.5
左旋回 直径12m

ラダー舵角

左旋回 直径12m

左旋回 直径12m

左旋回 直径12m

エルロン舵角

舵角無し

舵角無し

舵角無し

舵角無し

使用ゴム

#18輪ゴム×15本#18輪ゴム×15本#18輪ゴム×15本#18輪ゴム×15本

8条

8条

8条

8条

巻き数

340

340

340

340

0

0

0

0

3

3

3

3

巻き戻し数
発進速度

(m/s)

発進角度

(deg)

45

45

45

45

最大高度

(m)

16.9

25.2

27.5

24.5

57

表 6.6-3

サイドスラスト角の影響
ケース１

ケース5

ケース6

ケース7

プロペラに斜めから気流があたる影響

〇

〇

〇

〇

プロペラのジャイロ効果

〇

〇

〇

〇

〇

プロペラの慣性の影響

〇

〇

〇

〇

〇

（ベースライン）

ケース8
〇

全機質量

(g)

15

15

15

15

15

プロペラ 直径

(mm)

170

170

170

170

170

1.4

1.4

1.4

1.4

1.4

プロペラ中心-重心からのダウンスラスト

プロペラ ピッチ比
(deg)

3

3

3

3

3

サイドスラスト

(deg)

-0.8

0

-0.5

0.5

1

右旋回 直径12m

舵角無し

左旋回 直径12m

左旋回 直径12m

左旋回 直径12m

ラダー舵角
エルロン舵角

舵角無し

舵角無し

舵角無し

舵角無し

舵角無し

使用ゴム

#18輪ゴム×15本#18輪ゴム×15本#18輪ゴム×15本#18輪ゴム×15本#18輪ゴム×15本

8条

8条

8条

8条

8条

巻き数

340

340

340

340

340
0

巻き戻し数

0

0

0

0

発進速度

(m/s)

3

3

3

3

3

発進角度

(deg)

45

45

45

45

45

最大高度

(m)

27.5

27.8

23.7

27.5

26.1

ケース11

ケース12

表 6.6-4

発進角度の影響

ケース9

ケース10

ケース１
（ベースライン）プロペラに斜めから気流があたる影響

〇

〇

〇

〇

〇

プロペラのジャイロ効果

〇

〇

〇

〇

〇

プロペラの慣性の影響

〇

〇

〇

〇

〇

全機質量

(g)

15

15

15

15

15

プロペラ 直径

(mm)

170

170

170

170

170

1.4

1.4

1.4

1.4

1.4

プロペラ ピッチ比
プロペラ中心-重心からのダウンスラスト

(deg)

3

3

3

3

3

サイドスラスト

(deg)

0.5

0.5

0.5

0.5

0.5
左旋回 直径12m

ラダー舵角

左旋回 直径12m

左旋回 直径12m

左旋回 直径12m

左旋回 直径12m

エルロン舵角

舵角無し

舵角無し

舵角無し

舵角無し

舵角無し

使用ゴム

#18輪ゴム×15本#18輪ゴム×15本#18輪ゴム×15本#18輪ゴム×15本#18輪ゴム×15本

8条

8条

8条

8条

8条

巻き数

340

340

340

340

340
0

巻き戻し数

0

0

0

0

発進速度

(m/s)

3

3

3

3

3

発進角度

(deg)

25

35

45

55

65

最大高度

(m)

12.5

27.3

27.5

26.3

24.6

58

表 6.6-5発進角度の影響，FAI TAN Super Sport ゴム，3/32 インチ幅を使用
ケース１（ベースライン）

ケース13

ケース14

ケース15

ケース16

ケース17
〇

プロペラに斜めから気流があたる影響

〇

〇

〇

〇

〇

プロペラのジャイロ効果

〇

〇

〇

〇

〇

〇

プロペラの慣性の影響

〇

〇

〇

〇

〇

〇

全機質量
プロペラ

直径

(g)

15

15

15

15

15

15

(mm)

170

170

170

170

170

170
1.4

1.4

1.4

1.4

1.4

1.4

プロペラ中心-重心からのダウンスラスト

プロペラ ピッチ比
(deg)

3

3

3

3

3

3

サイドスラスト

(deg)

0.5

0.5

0.5

0.5

0.5

0.5
左旋回 直径12m

ラダー舵角

左旋回 直径12m

左旋回 直径12m

左旋回 直径12m

左旋回 直径12m

左旋回 直径12m

エルロン舵角

舵角無し

舵角無し

舵角無し

舵角無し

舵角無し

舵角無し

使用ゴム

#18輪ゴム×15本

FAI TAN 3/32"幅FAI TAN 3/32"幅FAI TAN 3/32"幅FAI TAN 3/32"幅FAI TAN 3/32"幅

8条

4条

4条

4条

4条

4条

巻き数

340

640

640

640

640

640

巻き戻し数

0

0

0

0

0

0

3

3

3

3

3

3

発進速度

(m/s)

発進角度

(deg)

45

45

50

60

70

75

最大高度

(m)

27.5

12.8

46.8

47.8

49.0

20.5

59

図 6.6-2３次元運動解析結果 – ケース 1（ベースライン）60

図 6.6-3３次元運動解析結果 – ケース 1a，プロペラに気流が斜めにあたる影響を無視

61

図 6.6-4３次元運動解析結果 – ケース 1b，プロペラのジャイロ効果を無視

62

図 6.6-5３次元運動解析結果 – ケース 1c，プロペラの回転慣性を無視

63

図 6.6-6３次元運動解析結果 – ケース 2，プロペラのアップスラスト 1.0 度

64

図 6.6-7３次元運動解析結果 – ケース 3，プロペラのダウンスラスト 1.0 度

65

図 6.6-8

３次元運動解析結果 – ケース 4，プロペラのダウンスラスト 5.0 度

66

図 6.6-9３次元運動解析結果 – ケース 5，右ラダー舵角有り，サイドスラスト左 0.8 度

67

図 6.6-10

３次元運動解析結果 – ケース 6，ラダー舵角無し，サイドスラスト無し

68

図 6.6-11

３次元運動解析結果 – ケース 7，左ラダー舵角有り，サイドスラスト左 0.5 度

69

図 6.6-12

３次元運動解析結果 – ケース 8，左ラダー舵角有り，サイドスラスト右 1.0 度

70

図 6.6-13

３次元運動解析結果 – ケース 9，発進角度 25 度

71

図 6.6-14

３次元運動解析結果 – ３次元運動解析結果 – ケース 10，発進角度 35 度

72

図 6.6-15

３次元運動解析結果 – ３次元運動解析結果 – ケース 11，発進角度 55 度

73

図 6.6-16

３次元運動解析結果 – ３次元運動解析結果 – ケース 12，発進角度 65 度

74

図 6.6-17３次元運動解析結果 – ３次元運動解析結果 – ケース 13，FAI TAN ゴム，発進角度 45 度

75

図 6.6-18

３次元運動解析結果 – ３次元運動解析結果 – ケース 14，FAI TAN ゴム，発進角度 50 度

76

図 6.6-19

３次元運動解析結果 – ３次元運動解析結果 – ケース 15，FAI TAN ゴム，発進角度 60 度

77

図 6.6-20

３次元運動解析結果 – ３次元運動解析結果 – ケース 16，FAI TAN ゴム，発進角度 70 度

78

図 6.6-21

３次元運動解析結果 – ３次元運動解析結果 – ケース 17，FAI TAN ゴム，発進角度 75 度

79

図 6.6-22プロペラのパワー効果

80

7

結論
ゴム動力機の３次元運動解析を行うための詳細な計算式を示し，Octave による解析ツールを作成した．
実際のゴム動力機を例にして，３次元運動解析を行った．プロペラのパワー効果とジャイロ効果の影響が大きい

ことがわかった．ダウンスラスト角を決めておき，ラダー舵角とサイドスラスト角を調整することによって，飛行パターンを調整することができる．
解析結果は定性的には実際の飛行と整合しているが，解析による最大高度は実際より大きいようである．これは解析に使用した有害抵抗係数 CDP = 0.06 が小さすぎるためかもしれない．
今後は低翼機の解析を行って，本解析の妥当性を検証したい．

81

8

参考文献

[1] 滝敏美，「模型グライダーのハイスタート発進の３次元運動解析」，第 60 回飛行機シンポジウム講演集，2022
年 10 月 11 日(火)～10 月 13 日(木)，JSASS-2022-5039．
[2] 加藤寛一郎，大屋昭男，柄沢研治，「航空機力学入門」，東京大学出版会，1982 年．
[3] 金田弥奈，加藤寛一郎，「模型滑空機の最長時間滞空について」，日本航空宇宙学会誌，第 39 巻，第 453 号，1991 年 10 月，pp.539-544．
[4] C. D. Perkins and R. E. Hage, “Airplane Performance, Stability and Control,” John Wiley & Sons, Inc., 1949.
[5] A. Capitao Patrao, “Implementation 0f Blaed Element Momentum/Vortex Methods for the Design of Aero Engine Propellers,”
Chalmers University of Technology, 2017. Downloaded from: http://publications.lib.chalmers.se/publication/253423.
[6] E. E. Larrabee, “Design of Propellers for Motorsoarers,” NASA. Langley Res. Center The Sci. and Technol. of Low Speedand Motorless Flight, Pt. 1, pp.285-303, 1979.
[7] C. N. Adkins and R. H. Liebeck, “Design of Optimum Propellers,” J. Propulsion and Power, Vol.10, No.5, Sept-Oct., 1994,
pp.676-682.
[8] L. A. Viterna and R. D. Corrigan, “Fixed Pitch Rotor Performance of Large Horizontal Axis Wind Turbines,” Technical
Report, NASA Lewis Research Center, 1982.
[9] R. MacNeill and D. Verstraete, “Blade Element Momentum Theory Extended to Model Low Reynolds Number Propeller
Performance,” 20th Australasian Fluid Mechanics Conference, Perth, Australia, 5-8 December 2016.
[10] H. S. Ribner, “Propellers in Yaw,” NACA Wartime Report, ARR No.3L09, 1943.
[11] H. S. Ribner, “Formulas for Propellers in Yaw and Charts of the Side-Force Derivative,” Technical Report NACA-TR-819,
1945.
[12] Y. Leng, J.-M. Moschetta, T. Jardin, and M. Bronz, “An Analytical Model for Propeller Aerodynamic Efforts at High
Incidence,” (2019) In: 54th 3AF International Conference on Applied Aerodynamics, 25 March 2019 - 27 March 2019
(Paris, France).
[13] M. S. Selig, “Modeling Propeller Aerodynamics and Slipstream Eﬀects on Small UAVs in Realtime,” AIAA AtmosphericFlight Mechanics 2010 Conference, 2-5 Aug 2010, Toronto, Ontario, Canada, No.2010-7938.
[14] O. Gur, “Propeller Effect on Air-Vehicle,” 出典不明.
[15] B. W. McCormick, “Aerodynamics, Aeronautics, and Flight Mechanics,” John Wiley & Sons, New York, 2nd ed., 1995.
[16] 嶋田有三，佐々修一，
「増補版

飛行力学」
，森北出版，2021 年．

[17] 藤原邦男，「物理学序論としての力学」，東京大学出版会，1984 年．
[18] T. v. カルマン，M. A. ビオ，「工学における数学的方法（上巻）」，法政大学出版局，1954 年．
[19] M. Drela and H. Youngren, “AVL 3.36 User Primer,” 2012.
[20] K. Budziak, “Aerodynamic Analysis with Athena Vortex Lattice (AVL),” Hamburg University of Applied Sciences, 2015.
[21]
[22]

松田七美男，「Octave の精義（第２版）」，カットシステム，2019 年．
松岡辰郎，「平成２３年度
資料コンピュータアルゴリズム（松岡担当分

Octave による数値計算入門）講義

公開版」，名古屋大学大学院工学研究科物質制御工学専攻，2011 年．

[23] J. McArthur, “Aerodynamics of Wings at Low Reynolds Numbers: Boundary Layer Separation and Reattachment,” PhD
Thesis, University of Southern California, 2008.

82

付録１

ゴム動力機の３次元運動方程式を解く Octave スクリプト

clear%機体姿勢の描画，2023/6/15%レポート A 改訂，2023/5/22%全面的に書き換え，2023/5/16%%プロペラ横力の影響（McCormic の理論）を組み込み，2023/5/16%%動力ゴムのモデルを組み込み，2023/4/6%%スタブ・ティルトを組み込み，2023/3/9%%ジャイロ効果を組み込み，スラスト角の定義を変更，2023/3/30%%プロペラ回転慣性を考慮，2023/4/24
global rho Dprop mm_start thr_up_deg thr_side_deg U0 iv delta_a Xprop Yprop Zprop%初期値%初速
U0 = 3.0;%姿勢（deg）
Phi0 = 0.0;
Theta0 = 45.0;%プロペラ回転数（回/s）n0 = 15.0;%時間刻み
tmax =30.0;%飛行時間[sec]step = 10;%1 秒あたりのステップ数
nt = step*tmax+1;%
function dx = equilibrium(x,t)
global rho Dprop mm_start thr_up_deg thr_side_deg U0 iv delta_a Xprop Yprop Zprop%プロペラ効果で考慮するもの%プロペラに斜めから気流があたる影響：switch_incline = 1%プロペラのジャイロ効果：switch_jyro = 1%プロペラ慣性：switch_inertia = 1%プロペラ空転抵抗：switch_free_wheel = 1
switch_incline = 1;
switch_jyro = 1;
switch_inertia = 1;
switch_free_wheel = 1;
%定数g = 9.80665;

83

rho = 1.225;%% ゴムのデータ，単位は N, mm%#18 輪ゴム×8 条×15 本
L_rubber = 262.5;
ns = 340;
ns_minus = 0;
mm_start = ns - ns_minus;
NNS = ns/L_rubber;delta_NNB = 0.4;
a1_rubber = 34.62;
a2_rubber = -70.789;
a3_rubber = 60.897;
a4_rubber = -16.056;
TS = a1_rubber*NNS+a2_rubber*NNS^2+a3_rubber*NNS^3+a4_rubber*NNS^4;
c05_rubber = 6.470;
c1_rubber = -2.649;
c2_rubber = 0.2492;
b_rubber = (TS-c05_rubber*NNS^0.5-c1_rubber*NNS-c2_rubber*NNS^2)/delta_NNB^4;%%FAI TAN Super Sport, 3/32"幅，4 条%L_rubber = 280.0;%ns = 640;
%ns_minus = 0;
%mm_start = ns - ns_minus;%NNS = ns/L_rubber;%delta_NNB = 0.8;
%a1_rubber = 10.71;%a2_rubber = -6.226;%a3_rubber = 1.608;%a4_rubber = 0.0446;%TS = a1_rubber*NNS+a2_rubber*NNS^2+a3_rubber*NNS^3+a4_rubber*NNS^4;%c05_rubber = 4.487;%c1_rubber = -0.957;%c2_rubber = -0.0708;%b_rubber = (TS-c05_rubber*NNS^0.5-c1_rubber*NNS-c2_rubber*NNS^2)/delta_NNB^4;% ゴムの方向
iXrubber = 1.00;
iYrubber = 0.00;
iZrubber = 0.00;%機体質量（kg），慣性データm = 0.015;
IXX = 0.0595^2*m;84

IYY = 0.0884^2*m;
IZZ = 0.1044^2*m;
IXZ = 0.0229^2*m;%垂直尾翼取付角，(-)が右旋回
iv = 0.01*1.0;%エルロン舵角，(+)が右旋回
delta_a = 0.1*0.0;%スタブ・ティルト（deg），右翼下げが(+)
stab_tilt_deg = 0.0;
stab_tilt_rad = stab_tilt_deg/180*pi;%プロペラ・データ
Dprop = 0.170;
Iprop = 0.0034/0.170*Dprop*(0.050/0.170*Dprop)^2;
Xprop = 0.118;
Yprop = -0.000;
Zprop = 0.022;
thr_up_deg = -10.56 - 3.0;
thr_side_deg = 0.5;

% プロペラ性能曲線近似式 -- 5%キャンバー円弧翼型，Peckpolymers タイプ展開形，ピッチ比 1.4
WT1 =[0.130007914;
0.577394189;-0.000106511;
3.250049106];
BT1 =[0.016618866;
0.014569938;-1.007691246;-2.625542778];
WT2 =[-3.206976146 0.726816739 1.27170327 0.055661233;-3.059721399 6.920664457 -0.118687632 0.627164185;
0.528820649 0.348423837 0.751921175 1.223750341];
BT2 =[1.115162042;-3.497636093;-0.402494442];

85

WT3 =[7.210628493 -3.766789952 4.2846653;
7.239477546 0.144715966 -0.48386372];
BT3 =[0.87713492;-0.056013793];
Solidity = 0.124;%プロペラ空転
J0 = 1.586;
CT0 = -0.020;

iXprop = cos(thr_up_deg/180*pi)*cos(thr_side_deg/180*pi);
iYprop = cos(thr_up_deg/180*pi)*sin(thr_side_deg/180*pi);iZprop = -sin(thr_up_deg/180*pi);
gamma = acos(iXprop*iXrubber + iYprop*iYrubber + iZprop*iZrubber);%空力データSW = 0.029100;
b = 0.400;
cbar = 0.0733;
lt = 0.1629;
CLa = 4.671;
alpha_eq = 0.1004;
CLeq = 0.8;
alpha_0 = -0.0709;
alpha_stall = 0.1646;
CLmax = 1.1;
CD0 = 0.06;
CD_stall = 0.13023;
CD_max = 1.5;
k_cd = (CD_stall-CD0)/CLmax^2;
A1 = CD_max/2;
A2 = (CLmax-CD_max*sin(alpha_stall)*cos(alpha_stall))*sin(alpha_stall)/(cos(alpha_stall))^2;
B1 = CD_max;
B2 = (CD_stall-CD_max*(sin(alpha_stall))^2)/cos(alpha_stall);Cma = -1.089;
CLq = 8.067;
Cmq = -10.101;
Cyb = -0.468;

86

Cyr = 0.423;
Clb = -0.245;
Clp = -0.455;
Cldelta_a = 0.18;
Cnb = 0.093;
k1_Cyp = 1.2888;
k2_Cyp = -0.2829;
k1_Cnp = -0.3858;
k2_Cnp = -0.0349;
k1_Clr = 0.9928;
k2_Clr = 0.1286;
k1_Cnr = -0.2680;
k2_Cnr = -0.1225;
k3_Cnr = -0.1048;% Stab Tilt，H.S.揚力係数
CLa_hs_St_Sw = 0.49608;
alpha_hs0 = 0.15233;%
Vc2 = x(1)^2+x(2)^2+x(3)^2;
Vc = sqrt(Vc2);
alpha =atan2(x(3),x(1));
beta =atan2(x(2),x(1));
p = b*x(4)/2/x(1);q = cbar*x(5)/2/x(1);r = b*x(6)/2/x(1);Cyp = k1_Cyp*alpha+k2_Cyp;
Cy = Cyb*(beta-iv)+Cyp*p+Cyr*r+CLa_hs_St_Sw*(alpha-alpha_hs0)*stab_tilt_rad;
Clr = k1_Clr*alpha+k2_Clr;
Cl = Clb*beta+Clp*p+Clr*r+Cldelta_a*delta_a;
Cm = Cma*(alpha-alpha_eq)+Cmq*q;
Cnp = k1_Cnp*alpha+k2_Cnp;
Cnr = k1_Cnr*alpha^2+k2_Cnr*alpha+k3_Cnr;
Cn = Cnb*beta+Cnp*p+Cnr*r+0.328*Cyb*iv-lt/b*CLa_hs_St_Sw*(alpha-alpha_hs0)*stab_tilt_rad;qq = 0.5*rho*Vc2;87

CL = CLa*(alpha-alpha_0)+CLq*q;
CD = CD0+k_cd*CL^2;if(alpha>alpha_stall)
CL = A1*sin(2*alpha)+A2*(cos(alpha))^2/sin(alpha);
CD = B1*(sin(alpha))^2+B2*cos(alpha);endif

L = qq*CL*SW;
D = qq*CD*SW;
Xa = L*sin(alpha)-D*cos(alpha);
Za = -L*cos(alpha)-D*sin(alpha);
Ya = qq*SW*Cy;
M = qq*cbar*SW*Cm;
L_moment = qq*b*SW*Cl;
N = qq*b*SW*Cn;%動力ゴムのモデル：新トルクカーブ
mm = mm_start-x(13);%ゴムの巻き数がある場合
if (mm>=0)
NN = mm/L_rubber;
Qr = c05_rubber*NN^0.5+c1_rubber*NN+c2_rubber*NN^2;
if (NN > NNS-delta_NNB)
Qr = Qr + b_rubber*(NN-(NNS-delta_NNB))^4;
endif
Qr = Qr/1000;%動力ゴムのモデル

ここまで

alpha_P = acos((iXprop*x(1)+iYprop*x(2)+iZprop*x(3))/Vc);
J = Vc*cos(alpha_P)/x(14)/Dprop;%プロペラ荷重計算
nT1 = WT1*J+BT1;
aT1 = tanh(nT1);
nT2 = WT2*aT1+BT2;
aT2 = tanh(nT2);
nT3 = WT3*aT2+BT3;
CT = nT3(1)*0.2;
CQ = nT3(2)*0.05;

88T = CT*rho*x(14)^2*Dprop^4;
Q = CQ*rho*x(14)^2*Dprop^5;% プロペラに作用する横力の計算 -- McCormick の理論
A_prop = pi*Dprop^2/4;
tau = 1/2*(1+sqrt(1+T/qq/A_prop));cl_av = 3*T/A_prop/qq/Solidity/(pi/J)^2;
cd = 8*(Q/A_prop/qq/Solidity/Dprop/(pi/J)^2-cl_av*J/6/pi);
PN = Solidity*qq*A_prop/2*tau*( cl_av + 2*pi*J/2/pi*tau* log(1+(pi/J/tau)^2) + pi/J/tau*cd )*sin(alpha_P);
NPP=

Solidity*qq*A_prop*Dprop/2/2*(

2*pi*cl_av/3/J+

2*pi/2*(1-2/3*tau^3*(J/pi)^2*log(1+(pi/J/tau)^2))-

pi/J*cd )*sin(alpha_P);
if (switch_incline == 0)
PN = 0.0;
NPP = 0.0;
endif%%
%プロペラのジャイロ効果
omega = 2.0*pi*x(14);MGx = (-sin(thr_up_deg/180*pi)*x(5) - sin(thr_side_deg/180*pi)*cos(thr_up_deg/180*pi)*x(6))*Iprop*omega;
MGy = (sin(thr_up_deg/180*pi)*x(4) + cos(thr_side_deg/180*pi)*cos(thr_up_deg/180*pi)*x(6))*Iprop*omega;MGz=(sin(thr_side_deg/180*pi)*cos(thr_up_deg/180*pi)*x(4)cos(thr_side_deg/180*pi)*cos(thr_up_deg/180*pi)*x(5))*Iprop*omega;
if (switch_jyro == 0)
MGx = 0.0;
MGy = 0.0;
MGz = 0.0;
endif%
%プロペラのジャイロ効果おわり
endif%ゴムの巻き数がある場合の計算おわり%ゴムの巻き数がなくなった%空転プロペラの抵抗だけがプロペラの影響
if (mm < 0)
Qr = 0.0;
alpha_P =0.0;
T = CT0*rho*Vc^2/J0^2*Dprop^2;
Q = 0.0;89-PN = 0.0;
NPP = 0.0;
MGx = 0.0;
MGy = 0.0;
MGz = 0.0;
if (switch_free_wheel == 0)
T = 0.0;
endif
%
endif
%ゴムの巻き数なしの計算おわり%機体に作用する力 -- プロペラの影響のまとめ
kX = x(2)*iZprop - x(3)*iYprop;
kY = x(3)*iXprop - x(1)*iZprop;
kZ = x(1)*iYprop - x(2)*iXprop;
kkk = sqrt(kX^2+kY^2+kZ^2);
kX = kX/kkk;
kY = kY/kkk;
kZ = kZ/kkk;
jX = kY*iZprop - kZ*iYprop;
jY = kZ*iXprop - kX*iZprop;
jZ = kX*iYprop - kY*iXprop;
Xp = T*iXprop+PN*jX;
Yp = T*iYprop+PN*jY;
Zp = T*iZprop+PN*jZ;
Lp = -Yp*Zprop+Zp*Yprop -Qr*iXprop/cos(gamma) +NPP*jX*(1) -MGx;
Mp =Xp*Zprop-Zp*Xprop -Qr*iYprop/cos(gamma) +NPP*jY*(1) -MGy;

Np = -Xp*Yprop+Yp*Xprop -Qr*iZprop/cos(gamma) +NPP*jZ*(1) -MGz;%機体に作用する力 -- プロペラの影響おわり

dx(1) = -g*sin(x(8))+(Xa+Xp)/m-x(5)*x(3)+x(6)*x(2);dx(2) = g*cos(x(8))*sin(x(7))+(Ya+Yp)/m-x(6)*x(1)+x(4)*x(3);dx(3) = g*cos(x(8))*cos(x(7))+(Za+Zp)/m-x(4)*x(2)+x(5)*x(1);
dx(4)

=

1/(IXX*IZZ-IXZ^2)*((-IZZ^2-IXZ^2+IYY*IZZ)*x(5)*x(6)+IXZ*(IXX-

IYY+IZZ)*x(4)*x(5)+IZZ*(L_moment+Lp)+IXZ*(N+Np) );dx(5) = 1/IYY*( -(IXX-IZZ)*x(6)*x(4)-IXZ*(x(4)^2-x(6)^2)+(M+Mp) );
dx(6)

=

1/(IXX*IZZ-IXZ^2)*(

IXX*IYY+IXZ^2)*x(4)*x(5)+IXZ*(L_moment+Mp)+IXX*(N+Np) );dx(7) = x(4)+x(5)*sin(x(7))*tan(x(8))+x(6)*cos(x(7))*tan(x(8));

90

IXZ*(-IXX+IYY-IZZ)*x(5)*x(6)+(IXX^2-

dx(8) = x(5)*cos(x(7))-x(6)*sin(x(7));dx(9) = x(5)*sin(x(7))*sec(x(8))+x(6)*cos(x(7))*sec(x(8));
dx(10)=

x(1)*cos(x(8))*cos(x(9))+x(2)*(sin(x(7))*sin(x(8))*cos(x(9))-

cos(x(7))*sin(x(9)))+x(3)*(cos(x(7))*sin(x(8))*cos(x(9))+sin(x(7))*sin(x(9)));dx(11)=

x(1)*cos(x(8))*sin(x(9))+x(2)*(sin(x(7))*sin(x(8))*sin(x(9))-

cos(x(7))*cos(x(9)))+x(3)*(cos(x(7))*sin(x(8))*sin(x(9))-sin(x(7))*cos(x(9)));dx(12) = -x(1)*sin(x(8))+x(2)*sin(x(7))*cos(x(8))+x(3)*cos(x(7))*cos(x(8));
dx(13) = x(14);dx(14) = (Qr/cos(gamma)-Q)/2/pi/Iprop;
if (switch_inertia == 0)dx(14) = (Qr/cos(gamma)-Q)/2/pi/(Iprop*0.001);
endif
endfunction

x0 = [U0, 0.0, 0.0, 0.0, 0.0, 0.0, Phi0/180*pi, Theta0/180*pi, 0.0, 0.0, 0.0, -1.5, 0.0, n0];t = linspace(0.0,tmax,nt);
x = lsode(@equilibrium,x0,t)% グラフ作図
subplot(2,4,1)
plot(x(:,10),-x(:,12))
grid on;subplot(2,4,2)
plot(x(:,11),-x(:,12))
grid on;subplot(2,4,3)
plot(t,-x(:,12))
grid on;subplot(2,4,4)
plot(t,x(:,1))
grid on;subplot(2,4,5)
plot(x(:,10),-x(:,11))grid on;subplot(2,4,6)
plot3(x(:,10),-x(:,11),-x(:,12))
grid on;subplot(2,4,7)

91

plot(t,mm_start-x(:,13))
grid on;subplot(2,4,8)
plot(t,x(:,14))
grid on;% グラフの終了%飛行軌跡と機体姿勢の図示
figure(); hold on; box('off')
axis equal;set(gca,'YDir','rev','ZDir','rev'); %飛行機でよく使われる右手系にする
xlabel('x');ylabel('y');zlabel('z');
grid on;view (3)%見やすくなるように適宜調整するinterval = step*1; %飛行機の表示間隔
scale = 0.4; %飛行機のサイズ
for i=1:interval:tmax*step
aircraft_figure_a(x(i,10),x(i,11),x(i,12),x(i,9),x(i,8),x(i,7),scale);
stem3(x(1:interval:i,10),x(1:interval:i,11),x(1:interval:i,12),'.:m');plot3(x(1:1:i,10),x(1:1:i,11),x(1:1:i,12));refresh;end===========================================
機体の図を描くスクリプト：aircraft_figure_a.m===========================================% 機体の図function [] = aircraft_figure(x,y,z,Psi,Theta,Phi,scale)%飛行機の 3D グラフィックを描画する関数%引数は重心座標，ロール角，ピッチ角，ヨー角，飛行機の寸法%ポリゴンの頂点を定義して，line オブジェクトで繋ぐ%回転行列を定義したスクリプトを読み込む
source('Rotation_Matrix.m');

p1 = [0+1.8, 0, 0];p2 = [-1.3+1.8, 0, 0];92p3 = [-2.3+1.8, 0, 0];
p4 = [-4+1.8, 0, 0];
p5 = [-4-0.8+1.8, 0, 0];
p6 = p2 + [0, 0, -0.2];
p7 = p3 + [0, 0, -0.1];
p8 = p4 + [-0.2, 0, -1];
p9 = p5 + [0, 0, -1];
p10 = p6 + [0, 3, -0.5];
p11 = p7 + [0, 3, -0.5];
p12 = p6 + [0, -3, -0.5];
p13 = p7 + [0, -3, -0.5];
p14 = p4 + [0, 1.5, 0];
p15 = p5 + [0, 1.5, 0];
p16 = p4 + [0, -1.5, 0];
p17 = p5 + [0, -1.5, 0];
p21 = [1.8, 1.2*cos(pi*0/6), -1.2*sin(pi*0/6)];
p22 = [1.8, 1.2*cos(pi*1/6), -1.2*sin(pi*1/6)];
p23 = [1.8, 1.2*cos(pi*2/6), -1.2*sin(pi*2/6)];
p24 = [1.8, 1.2*cos(pi*3/6), -1.2*sin(pi*3/6)];
p25 = [1.8, 1.2*cos(pi*4/6), -1.2*sin(pi*4/6)];
p26 = [1.8, 1.2*cos(pi*5/6), -1.2*sin(pi*5/6)];
p27 = [1.8, 1.2*cos(pi*6/6), -1.2*sin(pi*6/6)];
p28 = [1.8, 1.2*cos(pi*7/6), -1.2*sin(pi*7/6)];
p29 = [1.8, 1.2*cos(pi*8/6), -1.2*sin(pi*8/6)];
p30 = [1.8, 1.2*cos(pi*9/6), -1.2*sin(pi*9/6)];
p31 = [1.8, 1.2*cos(pi*10/6), -1.2*sin(pi*10/6)];
p32 = [1.8, 1.2*cos(pi*11/6), -1.2*sin(pi*11/6)];
p33 = [1.8, 1.2*cos(pi*12/6), -1.2*sin(pi*12/6)];
w1 = vertcat(p1,p2,p3,p4,p5);
w2 = vertcat(p2,p6,p7,p3);
w3 = vertcat(p6,p12,p13,p7,p11,p10,p6);
w4 = vertcat(p4,p16,p17,p5,p15,p14,p4);
w5 = vertcat(p4,p8,p9,p5);
w6 = vertcat(p21,p22,p23,p24,p25,p26,p27,p28,p29,p30,p31,p32,p33);%オブジェクトのスケールを変更する
w1 = w1 * scale;w2 = w2 * scale;
w3 = w3 * scale;w4 = w4 * scale;93w5 = w5 * scale;
w6 = w6 * scale;%オブジェクトを回転させる
DCM = Rotation_X(-Psi) * Rotation_Y(-Theta) * Rotation_Z(-Phi);
w1 = (DCM * w1')';
w2 = (DCM * w2')';
w3 = (DCM * w3')';w4 = (DCM * w4')';w5 = (DCM * w5')';w6 = (DCM * w6')';%オブジェクトを平行移動させる
w1= translational_shift(w1,x,y,z);

w2= translational_shift(w2,x,y,z);

w3= translational_shift(w3,x,y,z);

w4= translational_shift(w4,x,y,z);

w5= translational_shift(w5,x,y,z);

w6= translational_shift(w6,x,y,z);%描画する
line(w1(:,1),w1(:,2),w1(:,3));line(w2(:,1),w2(:,2),w2(:,3));line(w3(:,1),w3(:,2),w3(:,3));line(w4(:,1),w4(:,2),w4(:,3));line(w5(:,1),w5(:,2),w5(:,3));line(w6(:,1),w6(:,2),w6(:,3));end
function result = translational_shift(object,x,y,z)%頂点群 object を[x,y,z]だけ平行移動させる
shift_x = x*ones(rows(object),1);
shift_y = y*ones(rows(object),1);
shift_z = z*ones(rows(object),1);
shift= horzcat(shift_x,shift_y,shift_z);

result = object + shift;end=======================================================
回転マトリックスを定義するスクリプト：Rotation_Matrix.m=======================================================
1; %各軸周りの回転行列の定義

94

%x 軸回転function R_x = Rotation_X(Psi)
R_x = [ cos(Psi), sin(Psi), 0;-sin(Psi), cos(Psi), 0;0,0, 1 ];

end%y 軸回転function R_y = Rotation_Y(Theta)
R_y = [ cos(Theta), 0, -sin(Theta);
0, 1,

0;

sin(Theta), 0,cos(Theta) ];end%z 軸回転function R_z = Rotation_Z(Phi)R_z = [1,

0,
0,

0;
cos(Phi), sin(Phi);

0, -sin(Phi), cos(Phi) ];
end

95

付録２

機体の図を描くスクリプトの解説記事

96

97

98

99

100

101

102

103

104