人力飛行機の主翼設計を銀本より良くするたった一つの方法

タイトルは流行りの煽り文句です.9割嘘です.自分で書いてイライラしますね.

下にも同じリンク貼ってるけど,今回作ったプログラム

プログラム:https://gist.github.com/ina111/5053876

解説:https://gist.github.com/ina111/5053903

(追記)pythonで書き直しました。こちらの方が現代では使いやすいと思います。

https://gist.github.com/ina111/f0cedfb1c9e9055af72f8af26f31daca

(追記終わり)

鳥人間関係で面白いネタを教えてもらったので,議論するのもなんなのでサクッと作ってみました.

人力飛行機など(効率重視の飛行機ならなんでも)の主翼の設計をする時の話です.

まえがき

飛行機を低出力で飛ばすためには揚抗比というのが大事になります.揚力の大きさに対しての抗力を小さくしろってことです.人力飛行機では抗力の中でも誘導抗力というものが全体の1/3もあります.これは翼の上下面で圧力差が出来ることから翼端で渦発生し,その生成エネルギーで消費されるものです.

この誘導抗力を最小にする主翼の平面形というのは航空力学の教科書に載っています.それが楕円循環分布と呼ばれる.揚力の分布が楕円状になるものです.昔の航空機で主翼の平面形が楕円形になっているものがあったりします.循環っていうのは揚力や誘導抗力などを上手く数式で表せる便利なやつです.

主翼の設計方法

十分な揚力が得られる翼面積を保ちつつ,誘導抗力が小さくなるように楕円循環分布に従うように主翼の平面形を決めていきます.

教科書に書いてるこの方法が,イケテナイです.

楕円循環分布は揚力とスパン幅が一定の場合の誘導抗力が最小の循環分布です.つまり,空力面だけでの最適な形です.

空力だけではなく構造も含めて最適な主翼が設計したいなら曲げモーメントと誘導抗力を同時に制限して最適問題を解かないといけません.

そのことを書いた論文がR.T.JonesのNACA(NASAの前身)のテクニカルノート(TN-2249,”The spanwise disribution of lift for minimum induced drag of wings having a given lift and a given bending moment”)です.これは平面翼について解析的に書いてあります.

これを非平面翼で数値計算で解いて日本語で書いてあるのが浅井先生のNAL(JAXAの前身)のテクニカルレポート(TR-797,”非平面翼の最適設計-揚力と翼根曲げモーメントを与えた時の最小誘導抵抗-“)です.

プログラム

このTR-797をMatlab(Octave)で実装してみました.200行程度だったのでGithubのサービスの中のGistを使ってみました.詳しくは論文と下の解説を読んで下さい.論文を読んで下さい(大事なので2回)

プログラム:https://gist.github.com/ina111/5053876

解説:https://gist.github.com/ina111/5053903

プログラム中のbetaの値を0.9にしています.これは翼根での曲げモーメントを楕円循環分布での値から0.9倍にした制限をかけた時の誘導抗力が最小になる揚力の分布が以下になります.

翼根のあたりでより多くの揚力がが必要な分布になっています.

この条件(Gistに載せているデフォルトの条件)のときはbeta=0.9で飛行効率0.92となりました.だいたい1[N]ぐらいの違いでしょうか.この違いが構造重量で吸収できるようなら,楕円循環分布を捨てて,こちらの循環分布を用いるべきなことがわかります.

実際は構造を上手く攻めたり,できるかどうか,あと複数パラメータ振ってみての比較をするかどうかで,採用できるかどうか変わってくるかと思います.

感想

揚力線理論なプログラムは組んだこと合ったけど,渦格子法っぽいプログラムは初めて組んだので,これで人力飛行機なプログラムはかなり満足しました.やり残したことキチンと無くして社会人になりたいですね.

翼型とかの図面出力について、C#で作ってみた

まえがき

前の記事の続き

プロペラの渦法のGUIソフトを引き続き作っています。空いた時間にやってるのでのんびりのんびり・・・

現在のところ渦法とLarabeeの設計法、Adkins & Liebeckの設計法(翼素運動量理論Blade Element Momentum theory:BEMとかとも呼ばれたりする)と設計した翼の性能解析を実装しました。基本的には以前に作ったものの移植です。
機能は整ってきたので公開しようとおもうのですが、クラス設計を見直したいのでソースが整ってきたら公開します。

そんな中

製作段階になると、設計できたプロペラの翼型を印刷したいという要求があります。作っている人達は大抵、翼型データベースから持ってきた翼型datファイルをCADに持ってきて印刷したりしてるかと思います。
設計結果からExcelなどで翼型datファイルを拡大、回転させてCADにインポートして印刷してたりするんでしょうか。
業務としてやるなら、CADについてるスクリプトで書くのが普通でしょうか。
そこで、いま作っているプロペラ設計ソフトに翼型印刷機能を付けようと計画してコード書いていました。
しかし、ソフトに組み込むような機能では無いなと思ったので、ここで翼型印刷機能だけ別ソフトにして公開しておきます。
他記事
鳥人間の現役の人が翼型印刷について頑張って書いた記事は以下のようなものがあります。
設計の為の言語としてOctaveを勧めてしまった経緯もあって細かい技術的な工夫をして作っているのがわかります。

C#で作ってみた

Visual Studio使ってC#によって翼型datファイルの読み込みと拡大縮小・回転、印刷機能を作ってみました。
これを元に機能を追加することによって汎用的なソフトになるかと思います。
元ネタとしては自分が昔いたサークルで作られていたAirFoilPrinterというソフトです。ソース読んだことないけど、多分これをブラッシュアップしたものな気がします
車輪の再発明すぎるものではありますが、ネット上に挙がってなさそうでしたのでブログにしてみました。
ソースだけ晒しておきますが、もし欲しい人がいればexe形式で公開します。
ポイントとしては
csvなどのカンマ区切りではなく、スペース区切りになっている翼型datファイルを読み込むために数字の部分だけをx座標y座標のコレクションに追加している点。
AeroFoilクラスにDrawメソッドを作って
画面への描写と印刷面への描写に対応しているところです。
ミリメートル単位にしているので解像度など気にしないで大きさを直接指定できます。
あとはGUIを作る部分なので大した意味はないです。
188行目などのグラフィックの単位をミリメートルにした点がポイントです。
e.Graphics.PageUnit = GraphicsUnit.Millimeter;
あと見てもらいたいのは後半のAeroFoilクラスの翼型読み込みのためのメソッドと描写のメソッドの部分です。
AeroFoilクラスのポイント部分だけ抽出すると以下です。全文は一番下にあります。

C#で渦法のプロペラ設計プログラム作ってみた

前置き

ちょっと前まで、オープンソースというか誰かが作ったプロペラ設計プログラムや解析プログラムが中身分かる形でネット上に置いてあるということがなく、不便をしていました。そんな中、これまでAdkins & Liebeckの方法Larrabeeの方法のプロペラ設計法のプログラムを作って公開してきました。

それがXROTORがGPLライセンスになったことでオープンソースということでは便利な形で存在するようにました。これも日本語での説明どころか、英語ですら解説がないものだったので解説記事を書いたりインストール方法を書いたりしていました。

渦法

そうこうしている中、「低レイノルズ数プロペラの設計法」という論文を書かれている原田さんに直接色々と話しを伺うことが出来て、しかも、Matlabで作った論文の重要な部分を抜き出した簡易版のプログラムも頂けました。

渦法は低レイノルズ数領域のプロペラ設計では有用な設計方法ではありますが、(どんな点で有用かは論文や後述の課題見ればわかるはず、あとで説明書くかも)

  1. 論文が難しい
  2. 誰でも使えそうなソフトウェアの形になっているものがない

という点でLarrabeeやAdkins & Liebeckなどの方法が鳥人間でまだまだ使われています。

1.の点に関して、原田さんがMatlabの課題形式で渦法の設計方法のまとめの文章を書いてくださっています。鳥人間の関係者はちゃんとコンタクトを取ってみると課題が受け取れるかと思います。もしくは貰っている人がいるチームの人は担当者から引き継ぐなどできるかと思います。

課題を見ながらなら(しっかり時間をかければ)誰でも渦法の設計は出来るようになるかと思います。難しい論文の内容が自分の手元で組上がってくるのはこれ以上無く、楽しいものです。

2.の点に関して、もし誰でも使えるソフトというものがあれば、設計法の有用性から考えても渦法が広まるかと思っていました。(初期設定がLarrabeeの手法になっている)XROTORや私作成のプログラムの「渦法」版のソフトがあればいいなって話しです。

マニア向けに言うとXFOILの時代からXFLR5が出てきたことによる翼型解析の普及のようなことが起こればいいなと思うわけです。

なので作ってみた

C#ってラーメンタイマーみたいな練習以外で使うのは初めてでしたが、3日かけて(そのうち1日分の作業はブルースクリーンで闇に消えた)C#で作ってみました。パラメータを入れてCalculationボタンを押すと渦法で設計してくれるっていうものです。Matlabでの同じアルゴリズムの計算の3倍以上速く計算完了します。自分の環境だと10秒程度でした。forループが入れ子になって収束計算しているのでMatlabより早くなります。

今は単に原田さん謹製のMatlabプログラムからC#への移植しただけです。

もし反響があれば、以下のようなアップデートを行いつつ、公開していこうかと思います。

  • 複数データが比較可能へ
  • 性能解析モード搭載
  • 翼型性能をXFLR5から読み込み
  • ファイル入出力
  • etc…

 

気圧高度計(MS5611)をArduinoで使ってみた その1

前回までに以下の様な基板を作っていました。

MPU-6050基板作ってみた(前半)基板設計

MPU-6050基板作ってみた(後半)リフローしてみた

これは加速度ジャイロ地磁気、更に気圧計を載せて10軸の情報を得られるものです。ロケット姿勢・位置の計測用に使う予定です。当初の予定では秋月電子にも売っているMPL115A2というセンサを使う予定でしたが、さらに高精度な気圧計があるということでピンコンパチなこともあって、そちらのセンサを試してみました。

measurement社のMS5611-01BA03というものです。系列としては秋月電子などにモジュールとして売られているパララックス社高度計測モジュールに載っている気圧高度計のセンサ(MS5607)の次世代バージョンのものになります。

高分解で10cmまで測れると宣言していて本当かよ、とにわかに信じがたいので実際の動作を見てみました。

マイコンはArduinoでI2Cで接続しています。SPIでも繋げられるのでSPIの方が良さそうですが、今回の基板がそもそもMPU-6050というI2Cで繋がる基板への追加部品ということでI2Cでつなげています。

写真のように、5VのArduinoと繋がるようにセンサボートの前にI2Cのレベル変換モジュールを入れています、そしてSDカードに記録するために抵抗分圧で3.3Vで駆動するSDカードとつなげています。SDカードはArduinoにライブラリがあるのでそのまま使いました。すごくお手軽にうごくので嬉しいです。

MS5611でも秋月に売っているパララックス社の高度計測モジュールでも動くライブラリは存在しますが、センサがどのように動いているか理解のために自分でコード書いてみました。

Arduinoの既存のライブラリとしては

パララックス社が公開しているものFreeIMUを作っている方のブログに公開されているもの

があります。他にもmbedのライブラリとしてセニオネットワークさんが作っておられるライブラリとサンプルがあります。このmbedのライブラリはプロの仕事なので抽象度も高くキレイにまとまっていて勉強になります。

サンプルライブラリ

上のライブラリの方がオススメですが、自分で作ったものも公開しておきます。今回はシリアル通信でPCに表示させるところまでです。このあと、ログをSDカードに保存したり、生データにフィルタをかませたり、加速度ジャイロと合わせたりします。

/*
MS561101BA用スケッチ
MS5607でも使えるはず(たぶん)
I2Cでセンサデータを取得し、シリアル通信でPCに表示
by ina111
2012/07/16
*/

#include <Wire.h>

//気圧計のアドレス
#define MS5611_ADDR 0x76 //CBR =HIGHの時は0x76,LOWの時は0x77
//#define MS5607_ADDR 0x76

//気圧計で使う変数
uint16_t C_[6] ={40127, 36924, 23317, 23282, 33464, 28312}; //初期値
uint32_t D1=0,D2=0;
int32_t dT,TEMP=0;
int64_t OFF,SENS;
int64_t T2,OFF2,SENS2;
int32_t P=0;
int32_t Height;

unsigned long now;                   //現在時間を入れて変換時間を計算
unsigned long lastD1Conv,lastD2Conv; //最後に変換した時間micros()
unsigned long ConvTime = 10000;      //変換時間[マイクロ秒]
boolean SWD1Conv = true, SWD2Conv = true; //D1,D2どっちを処理しているか
boolean SWD1D2 = true;               //trueでD1の処理,falseでD2の処理

void setup(){
  Wire.begin();
  Serial.begin(9600);
  //Serial.begin(115200);
  delay(100);
  readPROM();
}

void loop(){
  getD1();
  getD2();
  getPressTemp(D1,D2);
  getHeight(P);

  Serial.print(millis());Serial.print("\t");
  Serial.print(TEMP);Serial.print("\t");
  Serial.print(P);Serial.print("\t");
  Serial.println(Height);
}

/*
 * func name  : getD1,getD2
 * processing : MS5611の生データD1,D2を取得
 * param      : 
 * summary    : 後述のstartConvとreadADC使用
 *              変換の時間前だと変換、変換後なら読み出しを行う
 *              D1とD2は交互に読み出されるようにスイッチしている
 * return     : 
 */
void getD1(){
  now = micros();
  //SWD1D2==trueだとD1,falseだとD2
  //SWD1Conv==trueだと変換、falseかつ時間経過後は読み出し
  if(SWD1Conv == true && SWD1D2 == true){
    startConv(0x48); //D1のとき0x48,D2のとき0x58,OSRによって変化
    lastD1Conv = micros();
    SWD1Conv = false;
  }else if(now - lastD1Conv >= ConvTime && SWD1D2 == true){
    D1 = readADC();
    SWD1Conv = true;
    SWD1D2 = false;
  }
}

void getD2(){
  now = micros();
  if(SWD2Conv == true && SWD1D2 == false){
    startConv(0x58); //D1のとき0x48,D2のとき0x58,OSRによって変化
    lastD2Conv = micros();
    SWD2Conv = false;
  }else if(now - lastD2Conv >= ConvTime && SWD1D2 == false){
    D2 = readADC();
    SWD2Conv = true;
    SWD1D2 = true;
  }
}

/*
 * func name  : getPressTemp,getHeight
 * processing : MS5611の生データから温度TEMP、気圧P、高度Heightを計算
 * param      : ADC1,ADC2    / getD1,getD2で得られたD1,D2
 *              hPa          / getPressTempで得られたPをhPaにしたもの
 * summary    : MS5611のデータシートによる計算,と高度気圧式から線形化
 * return     : 
 */
void getPressTemp(uint32_t ADC1, uint32_t ADC2){
  dT   = (int32_t)(ADC2 - ((int32_t)C_[4] << 8));
  TEMP = 2000 + ((dT * (int64_t)C_[5]) >> 23);
  OFF  = (((int64_t)C_[1]) << 16) + (((int64_t)C_[3] * dT) >> 7);
  SENS = (((int64_t)C_[0]) << 15) + (((int64_t)C_[2] * dT) >> 8);
  P    = ((ADC1 * SENS) >> 21) - OFF >> 15;

  if (TEMP < 2000) {
    T2    = (dT * dT) >> 31;
    OFF2  = 5 * (TEMP - 2000) * (TEMP - 2000) >> 1;
    SENS2 = 5 * (TEMP - 2000) * (TEMP - 2000) >> 2;
    TEMP = TEMP - T2;
    OFF  = OFF - OFF2;
    SENS = SENS - SENS2;
  }  
}

void getHeight(int32_t hPa){
  //Height = -938.502 * hPa/100.0 + 948697; //t0=30[deg]で1000mまでの線形近似
  float t0 =30.0;float P0 = 1013.25;
  Height = 153.8*(t0+273.2)*(1-pow((hPa/100.0/P0),0.1902));
}

/*
 * func name  : startConv,readADC
 * processing : MS5611の内部ADCの変換開始関数
 *              変換後OSRの時間はアクセス不可
 *              ADCの結果を読み込む関数。returnでuint32_tが出てくる
 * param      : command    / D1,D2,OSRによってアドレスが異なる
 *                         / OSR(Over Sampling Ratio)             
 * summary    : startConv(command) -> delay(10) -> readADC()
 *              変換コマンドを送信
 *              値呼び出しコマンドを送信、受信
 * return     : 
 *              conversion    / ADCした値.D1,D2のこと
 */
void startConv(uint8_t command){
  Wire.beginTransmission(MS5611_ADDR);
  Wire.write(command);
  Wire.endTransmission();
}

uint32_t readADC(){
  uint32_t conversion = 0;
  // start read sequence
  Wire.beginTransmission(MS5611_ADDR);
  Wire.write(0x00);
  Wire.endTransmission();

  Wire.beginTransmission(MS5611_ADDR);
  Wire.requestFrom(MS5611_ADDR, 3); //3byteリクエスト
  if(Wire.available()){
    conversion = Wire.read() * 65536 + Wire.read() * 256 + Wire.read();
  }
  return conversion;
}

/*
 * func name  : readPROM
 * processing : MS5611の工場校正値を読み込み
 * param      : 
 * summary    : C_[i]に0~5読み込む(C1~C6とは一つズレている)
 * return     : 
 */
void readPROM(){ 
  for(int i=0; i<6; i++){
    Wire.beginTransmission(MS5611_ADDR);
    Wire.write(0xA2 + i*2); //PROMの最初 データシートでは0xA0に見えるが0xA2から   
    Wire.endTransmission();

    Wire.beginTransmission(MS5611_ADDR);
    Wire.requestFrom(MS5611_ADDR,2);
    if(Wire.available() >= 2){
      C_[i] = Wire.read() * 256 + Wire.read();
    }
  }
}

 

人力飛行機が飛ぶところを360度カメラで撮ってみた

正確には自分はカメラを貸しただけで『撮ってもらった』ものです。

googleストリートビューみたいな360度方向が見れるものが大好きで胸がワクワクします。ストリートビュー見ていて一日が終わったとかそういう感じで毎日を過ごしています()。そんなこともあってソニーのカメラでBloggieという360度レンズを付けることのできるカメラを持っています。


http://www.sony.jp/bloggie/

ソニー公式サイトの作例もあるように、360度動画の画質はかなり悪いです。ナチュラルモザイクな雰囲気で怪しげな動画が撮れます(笑)。

今回は人力飛行機の方のサークル(東京工業大学”Meister”)の後輩にお願いして、鳥人間コンテストに向けて最終調整している試験飛行の様子を動画に撮ってもらいました。動画はソニー純正の画像ソフトPMBというソフトを使うことで、マウスを動かすとそっちの方向を向くことが出来る360度動画を作ることができます。PMB上でマウスグリグリして遊んでいるところをキャプチャーして動画にしてみました。エンコードやyoutubeにアップロードしての画質の劣化は最小限にしました。撮影の段階でこの程度のモザイク画質です。

まるで自分が(視力が悪い状態で)パイロットになったかのような迫力のある動画ができて楽しいし、飛行の軌跡もわかります。何かトラブルがあった場合も検証用動画として有効じゃないかと思いました。