2014-01-01

人力飛行機を割と本気で設計してみた。第一報

明けましておめでとうございます。
今年もよろしくお願いします。

去年2013は忘れられない年になりそうですね。私の設計した飛行機が琵琶湖上空を飛行し、旋回し、帰還し、記録に残らなかったというのは後々生きていく中で時々に語る事になるでしょう。それ以外にも色々あったなぁ2013年。

では、早速タイトルの通り、内容に移りましょう。

今回作ってみた機体は私の今持つ技術を使ったらどういう機体になるのか試してみたかったということ、そして来たる機体のたたき台として存在する。
また、この機体の設計を通して技術やノウハウを公開していくことによって、人力飛行機界全体の底上げができたらと考えている。

今回は主翼の空力設計及び構造設計について。

3DCADでも主翼だけ作ってみた。フレームは昔作ったTT機を大きさ比較のために置いている。スパン34m

主翼設計のコンセプトについて説明しよう。
出発点はこれ。
構造最優先設計
構造には考慮しなければならない事が非常に多い。曲げ荷重だけでなく捻り荷重、局所座屈（参考 ina111's blog)、またワイヤー
を張るなら圧縮荷重も考慮しなければならず、経験で積層構成を決める場所もおおい。こうしてある程度経験が活かされた積層構成に対して、空力が適合するような設計をすることでより軽く出力の小さい機体を目指そうというものである。
採用した主桁の特徴で目立つものは二つ。

ワイヤー単桁構造　　　　　　　ワイヤーの抵抗は言われるよりも小さいのではないかというところが出発点。
翼根が最も太い桁ではない　　　ワイヤー機において最も曲げモーメントがきついのは翼根ではなくワイヤー結節点

この二つの特徴を持つ主桁に対して、空力の側が採用した技術は以下の通り

最大たわみ量を構造制約とした多角形最小循環分布の探索　（後述）
XGAGを用いた翼型最適化　（後述）
捻り下げ　中央翼-無　内翼-有　外翼①-無　外翼②-有　最外翼-無　翼厚の不連続を出さずにできるだけ翼厚を薄くするため（及び後述のスパン変更のため）

具体的なところは、XFLR5のファイルを公開するのでこちらでを見て欲しい。尾翼周りは割と適当。
benishake_koe_97kg.wpa

また、積層構成や具体的な特徴をまとめたExcelファイルを公開する。桁設計については正直全然詰まってない。
紅鮭・肥　概要

最大たわみ量を構造制約とした多角形最小循環分布の探索

これが空力の側の切り札その一
実際に用いた.m（octaveファイル)を公開する。解答してgamma_designをoctaveで(たぶんMATLABでも可)を実行すればOK
gamma_design
また、このお試し版として、octaveの実行環境がなくても実行できる、pyQtを用いたGUIアプリケーションを開発している。こちらも今月中には公開できると思う。

中身の説明に映っていこう。
何度も紹介している元論文TR797はこちらの記事参照
人力飛行機の主翼設計を銀本より良くするたった一つの方法

インプットやら桁の剛性やら準備やらのところ、読み飛ばしてOKです。

clear all

%---------------二次構造重量を考慮すること

rho=1.184;
U=7.5;
Mass=97;

%-----ワイヤー取付位置(m)
y_wire=6.55;
%-----ワイヤー引っ張り力(N)
N_wire=485;

%-----最大たわみ(構造制限測定点でのたわみpl
max_tawami=2.5

%-----積層構成を入力
%-----片翼桁分割数
num_section=6;

for np=1:num_section
	ply_data{np}=dlmread(strcat('ply_data_',mat2str(np),'.csv'))(:,1:6);
	mandrel_data(np,3)=min(ply_data{np}(:,1));
	mandrel_data(np,2)=mandrel_data(np,3)*3;
	mandrel_data(np,1)=max(ply_data{np}(:,6))/1000;
end	


for i=1:num_section
	%-----桁長さ(m)
	spar_length(i)=mandrel_data(i,1);
	%-----桁外径
	spar_dia(i)=max(ply_data{i}(:,1));
end

%-----スパン
b=sum(spar_length)*2
%-----dS(m)
dy=0.05
for i=1:num_section
	y_section(i)=sum(spar_length(1:i));
end

i=0;
j=1;
do 
	i++;
	y_div(i)=dy*i;
	if y_div(i)>=y_section(j)
		Ndiv_sec(j)=i;
		j++;
	end
	if y_div(i) >=y_wire - dy/2  && y_div(i) <= y_wire + dy/2
		Ndiv_wire=i;
	end
until(max(y_div)>=b/2);



%----------------------------------------------------------------------------------構造関連

%-----カーボン物性	
%-----24tor40t
ply_data{np}(1:4,2)=24;
%-----カーボン物性(2列目24t,3列目40t)kgfmm
E=[0,24,40;0,13000,22000;45,1900,1900;90,900,800];
%-----プリプレグ密度(g/mm^3)
C_D=[24,0.001496;40,0.001559];
%-----形成厚(mm)
C_T=[24,0.125;40,0.111];
%-----二次部材密度(kg/m2)
Snd_rho=0.4;
%-----ポアソン比
Pois=0.3;
		
%-----剛性分布
for n=1:num_section
	if n==1
	%-----座標、積層データ、マンドレルデータ、翼番号、物性を送って各座標点での断面二次モーメントを計算する関数(kgf.m^2)
	I(1:Ndiv_sec(n))=akanenosora_Calc_I(y_div(1:Ndiv_sec(n))'.*1000,ply_data{n},mandrel_data,n,E,C_T)(:,2)';
	%-----座標、積層データ、マンドレルデータ、翼番号、物性を送って各座標点での曲げ剛性を計算する関数(kgf.m^2)
	EI(1:Ndiv_sec(n))=akanenosora_Calc_Hardness(y_div(1:Ndiv_sec(n))'.*1000,ply_data{n},mandrel_data,n,E,C_T)(:,2)'.*9.8;
	%-----各座標での線密度を計算する関数(g/m)
	sigma_wire(1:Ndiv_sec(n))=akanenosora_Calc_Linear_Density(y_div(1:Ndiv_sec(n))'.*1000,ply_data{n},mandrel_data,n,C_D,C_T)';
	else
	%-----座標、積層データ、マンドレルデータ、翼番号、物性を送って各座標点での断面二次モーメントを計算する関数(kgf.m^2)
	I(Ndiv_sec(n-1)+1:Ndiv_sec(n))=akanenosora_Calc_I(y_div(Ndiv_sec(n-1)+1:Ndiv_sec(n))'.*1000,ply_data{n},mandrel_data,n,E,C_T)(:,2)';
	%-----座標、積層データ、マンドレルデータ、翼番号、物性を送って各座標点での曲げ剛性を計算する関数(kgf.m^2)
	EI(Ndiv_sec(n-1)+1:Ndiv_sec(n))=akanenosora_Calc_Hardness(y_div(Ndiv_sec(n-1)+1:Ndiv_sec(n))'.*1000,ply_data{n},mandrel_data,n,E,C_T)(:,2)'.*9.8;
	%-----各座標での線密度を計算する関数(g/m)
	sigma_wire(Ndiv_sec(n-1)+1:Ndiv_sec(n))=akanenosora_Calc_Linear_Density(y_div(Ndiv_sec(n-1)+1:Ndiv_sec(n))'.*1000,ply_data{n},mandrel_data,n,C_D,C_T)';
	end
end

%-----断面係数
for n=1:num_section
if n==1
	Zz(1:Ndiv_sec(n))=I(1:Ndiv_sec(n))./(spar_dia(n)/2);
else
	Zz(Ndiv_sec(n-1)+1:Ndiv_sec(n))=I(Ndiv_sec(n-1)+1:Ndiv_sec(n))./(spar_dia(n)/2);
end
end




coe_tawami=max_tawami/(b/2)^2;
z_div=coe_tawami.*y_div.^2;

for i=1:length(y_div)
	if i!=1
		dS(i)=sqrt((y_div(i)-y_div(i-1))^2+(z_div(i)-z_div(i-1))^2)/2;
	else
		dS(i)=sqrt(y_div(i)^2+z_div(i)^2)/2;
	end
end 





%-----Control Point
for i=1:length(y_div)
	if i!=1
		y(i)=(y_div(i)+y_div(i-1))/2;
		z(i)=(z_div(i)+z_div(i-1))/2;
		phi(i)=atan((z_div(i)-z_div(i-1))/(y_div(i)-y_div(i-1)));
	else
		y(i)=y_div(i)/2;
		z(i)=z_div(i)/2;
		phi(i)=atan(z_div(i)/y_div(i));
	end
end

%-----桁重量
spar_weight=trapz(y',sigma_wire')*2

%-----線密度&ワイヤー引っ張り
sigma_wire(Ndiv_wire)=sigma_wire(Ndiv_wire)+N_wire/(dS(Ndiv_wire)*2);

本質はここから

以前記事にした多角形化行列を作成している

%-----循環分布多角形化行列
for i=1:Ndiv_sec(1)
	pol_g(i,1)=1;
	pol_g(i,num_section)=0;
end

for j=2:num_section
	for i=Ndiv_sec(j-1)+1:Ndiv_sec(j)
		pol_g(i,j-1)=(-(y(i)-y_section(j))/(y_section(j)-y_section(j-1)));
		pol_g(i,j)=((y(i)-y_section(j-1))/(y_section(j)-y_section(j-1)));
	end
end

循環分布に書けることで誘導速度を求める行列Qを作成している。やってることはViot-Savartの積分なだけなんだけれども、コードに落とすと結構複雑なので別ファイル化。またそのあとその行列を多角形化

%-----誘導速度行列
Q_ij=Calc_Q(y',z',phi',dS');
%-----多角形化
Q=transpose(pol_g)*Q_ij*pol_g;

ここから、構造制約の部分。texで書くの面倒くさいから言葉で説明するけれども、数値的な積分が行列で実行できるというのは忘れられやすいところ。この積分はオイラー積分だから精度は良くないけれども、循環分布を求める上では十分な精度だろう。とにかく、循環分布にかけることで、たわみ曲線を出力する行列を作っている。これ、超肝。

%-----せん断力qを求める行列
for j=length(y):-1:1
	for i=j:-1:1
		if j==i
			q(i,j)=dS(j);
		else
			q(i,j)=dS(j)*2;
		end
	end
end

%-----せん断力からモーメントを求める行列
M=q;
q=q.*rho.*U;
%-----たわみ角を求める行列
for i=length(y):-1:1
	for j=i:-1:1
		if j==i
			vd(i,j)=dS(j)/EI(j)*10^(6);
		else
			vd(i,j)=dS(j)*2/EI(j)*10^(6);
		end
	end
end

%-----たわみを求める行列
for i=length(y):-1:1
	for j=i:-1:1
		if j==i
			v(i,j)=dS(j);
		else
			v(i,j)=dS(j)*2;
		end
	end
end

%-----欲しいところの値
B_want(1,columns(y))=0;
B_want(1,Ndiv_sec(5))=1;

%-----構造制約条件行列
B=B_want*(v*vd*M*q*pol_g);
Bel=B_want*(v*vd*M*q);
%-----構造制約条件の値
B_val=max_tawami+B_want*(v*vd*M*q)*sigma_wire'./rho./U;

あとは元の論文の通り、最適化。ちょっと注意すべきは最適化行列を作るときに多角形化行列を右から、その転置を左から書けていること。また比較のために、多角形化、構造制約無し、構造制約無しの循環分布も同時に同時に作っている。あとはグラフィック関連

%-----揚力制約条件行列
C=4.*rho.*U.*(dS*pol_g);
Cel=4.*rho.*U.*(dS);
%-----揚力制約条件の値
C_val=Mass*9.8;
		
%------------------------構造無視
Ael_pol=Q_ij;
for i=1:columns(Q_ij)
	Ael_pol(:,i)=Ael_pol(:,i).*(dS'.*2);
end
Mat_indDel_pol=(Ael_pol.*rho);
Ael_pol=(Ael_pol+transpose(Ael_pol));
Ael_pol=rho.*(transpose(pol_g)*Ael_pol*pol_g);
Ael_pol(rows(Ael_pol)+1,:)=-C;
Ael_pol(1:rows(Ael_pol)-1,size(Ael_pol)(2)+1)=-C';

A_valel_pol(rows(Ael_pol),1)=-C_val;

G_el_pol=(Ael_pol\A_valel_pol)(1:rows(A_valel_pol)-1);
D_indel_pol=trapz(y',rho.*(Q_ij*pol_g*G_el_pol)./2.*(pol_g*G_el_pol))*2;
V_indel_pol=(Q_ij*pol_g*G_el_pol./2);

clear Ael Mat_indDel A_valel
%-----構造無視その2
Ael=Q_ij;
for i=1:columns(Q_ij)
	Ael(:,i)=Ael(:,i).*(dS'.*2);
end
Mat_indDel=(Ael.*rho);
Ael=(Ael+transpose(Ael));
Ael=rho.*(Ael);
Ael(rows(Ael)+1,:)=-Cel;
Ael(1:rows(Ael)-1,size(Ael)(2)+1)=-Cel';

A_valel(rows(Ael),1)=-C_val;

G_el=(Ael\A_valel)(1:rows(A_valel)-1);
D_indel=trapz(y',rho.*(Q_ij*G_el)./2.*(G_el))*2;

clear Ael Mat_indDel A_valel
%------------------------構造制約
%-----最適化行列
A=Q_ij;
for i=1:columns(Q_ij)
	A(:,i)=A(:,i).*(dS'.*2);
end
Mat_indD=(A.*rho);
A=(A+transpose(A));
A=rho.*(transpose(pol_g)*A*pol_g);
A(rows(A)+1,:)=-B;
A(rows(A)+1,:)=-C;
A(1:rows(A)-2,size(A)(2)+1)=-B';
A(1:rows(A)-2,size(A)(2)+1)=-C';

A_val(rows(A)-1,1)=-B_val;
A_val(rows(A),1)=-C_val;

G_sec=(A\A_val)(1:rows(A_val)-2)


D_ind=trapz(y',rho.*(Q_ij*pol_g*G_sec)./2.*(pol_g*G_sec))*2
Lift=C*G_sec
V_ind=(Q_ij*pol_g*G_sec./2);

%-----誘導速度を多角形近似
function res=pol_induced(pol_g,V_ind,Vind_sec)
	res=sum((V_ind-pol_g*Vind_sec).^2);
endfunction

Vind_sec0(1,length(Ndiv_sec))=0;
V_ind_sec=fminunc(@(Vind_sec)pol_induced(pol_g,V_ind,Vind_sec),Vind_sec0')

Vindel_sec0(1,length(Ndiv_sec))=0;
V_indel_sec=fminunc(@(Vindel_sec)pol_induced(pol_g,V_indel_pol,Vindel_sec),Vindel_sec0')

bending_angle = (vd*M*q*(pol_g*G_sec-sigma_wire'./U./rho)).*180./pi;
bending = v*vd*M*q*(pol_g*G_sec-sigma_wire'./U./rho);

figure(1)
subplot(2,1,1)
	plot(y',pol_g*G_sec)
	hold on
	plot(y',G_el,'r--')
	hold off
	grid on
	axis 'equal'
subplot(2,1,2)
	plot(y',atan((Q_ij*pol_g*G_sec./2)./U).*180./pi)
	grid on
	hold on
	plot(y',atan((pol_g*V_ind_sec)./U).*180./pi,'r--')
	grid on
	hold off
figure(2)
subplot(3,1,1)
	plot(y',(M*q*(pol_g*G_sec-sigma_wire'./U./rho))./Zz'.*1000);
	hold on
	plot(y(Ndiv_wire),((M*q*(pol_g*G_sec-sigma_wire'./U./rho))./Zz'.*1000)(Ndiv_wire),'ro')
	hold off
	grid on
subplot(3,1,2)
	plot(y',(vd*M*q*(pol_g*G_sec-sigma_wire'./U./rho)).*180./pi);
	hold on
	plot(y(Ndiv_wire),((vd*M*q*(pol_g*G_sec-sigma_wire'./U./rho)).*180./pi)(Ndiv_wire),'ro')
	hold off
	grid on
subplot(3,1,3)
	plot(y',(v*vd*M*q*(pol_g*G_sec-sigma_wire'./U./rho)));
	hold on
	plot(y(Ndiv_wire),(v*vd*M*q*(pol_g*G_sec-sigma_wire'./U./rho))(Ndiv_wire),'ro')
	plot(y',z','r--')
	hold off
	grid on
	axis 'equal'

一つの記事でぱっとまとめるつもりだったのに、長くなってしまったから、翼型の話とかは次の記事で。

個人的には、この最適化方法は圧倒的な影響力を持っていて、人力飛行機界において解析的に楕円循環分布を作り出す方法（楕円循環分布決めうち）を過去のものにしてしまうほどの力を持っていると思っている。この方法は言い換えると、「その桁で、理想の最大たわみを実現する循環分布を求める方法」であり、上反角が旋回にとって決定的な影響をもたらすタイムトライアル部門だけでなく、揚力損失とサイドスリップの板挟みになる距離部門にも有効だ。また片持ちの桁を持つ機体だけでなく、ワイヤーを張る機体も、どれくらいの張力がいいのか、そのときのたわみはどうなるのかぱっと分かるし、すごくいい。

実際、紅鮭・肥（機体名）の設計を通して、例えばワイヤー張力が足りなくて、循環の値が負になってしまったり、逆にワイヤー内が下に撓んでしまうような状態をぱっと見て判断し、フィードバックして設計する事ができた。また翼効率eの急激な低下を起こさないようなワイヤー張力を探したり、とにかく自由に設計できてしまう。

恐ろしい。そしてこの手法は本当に楽しい…

では次の第二報の記事にて、翼型の設計とかの話をしましょ。

2013-12-05

XGAGアップデート　&　SSS前刷原稿公開

こんにちは。yuukivelです。

XGAGをアップデートしました。元記事にもリンクを加えておきます。

https://dl.dropboxusercontent.com/u/36164244/XGAG_1.01_Win.zip

アップデート内容は

システムの安定化
表示の修正

となります。前者は問題児であった「a0_pwrt.dat」をこまめに消すようにしました。このファイルはXFOILからの結果を受け取るための一時フォルダなのですが、消去が間に合わずに次の解析結果受け取りに悪影響を及ぼすことがあるようです。起動時と解析開始時に消すコマンドを追加しました。

また、スカイスポーツシンポジウムにおける私の講演の前刷り原稿を公開します。こちらは用いた遺伝的アルゴリズムについて詳しく書いてありますので、自分でプログラム組もうと思ったときの参考になれば、と思います。

最適化アルゴリズムによる低レイノルズ数向け翼型設計の紹介

これからドキュメンテーションに積極的に取り組んでいきます。今回GUIの構築にpythonのGUIモジュールであるpyQtを用いましたが、こちらの技術紹介とか…
pyQt、日本語の参考サイトがメチャメチャ少ないんですよ。せっかく得たノウハウなので、あとで記事にまとめたいと思います。

あと、昔某SNSにて、勢いでTR797の改良版も公開するよ、と言いました。こちらもやります。

そうだ、せっかくだから、「構造制約を考慮した多角形最適循環分布」を求めるGUIアプリケーションの作成を通してpyQtの技術解説をしようか。

2013-11-30

XGAG　遺伝的アルゴリズムによる翼型設計GUIアプリケーション

翼型について

ご無沙汰しておりました。yuukivelです。

今回の記事は、遺伝的アルゴリズムによる翼型設計GUIアプリケーション
XGAG　XFOIL Genetic Algolithm Graphical user interface airfoil design tool
の公開、および使い方の解説です。

XGAGは以下のレポジトリよりダウンロードできます。Window向けに.exeを用意しています。Mac用.appは少々お待ち下さい　pythonを入れている人はソースコードも使用できます。（モジュールの要求あり）

Windows
https://dl.dropboxusercontent.com/u/36164244/XGAG_1.01_Win.zip

(前バージョン : https://dl.dropboxusercontent.com/u/36164244/XGAG_1.00_Win.zip)

ソースコードはこちらより閲覧できます。今後より読みやすくする予定です。
https://github.com/NaotoMorita/NFoilDesign/

また、11/30に行われたスカイスポーツシンポジウムで用いたスライドを公開します。
https://dl.dropboxusercontent.com/u/36164244/NaotoMORITA19thSSS%E7%99%BA%E8%A1%A8%E3%83%91%E3%83%AF%E3%83%9D.pptx

XGAGは、自分の機体にあった翼型を設計したい飛行機設計者だけでなく、翼型のことを理解したい鳥人間、および模型飛行機製作者など、わかりやすく、幅広く使えるように配慮しました。
特に人力飛行機開発チームにおいては、設計者が指定したパラメータでXGAGを用いて、各自できた翼型を持ち寄るといった使い方もできるでしょう。

もちろんこのソフト単体でも使うことができますが、既存の翼型解析アプリケーションであるXFLR5と組み合わせると非常に有用です。個人的には、XGAGがとっかかりとなって翼型や翼の事についての理解が深まる人が増えてくれればいいなと思っています。

では早速使い方にうつっていきましょう。

XGAGの使い方-設定編

まず翼型を集めましょう。以下のサイトをオススメします。
http://www.airfoildb.com/
・欲しい翼型を検索して、そのページを表示。Other formatのDATの中身を保存しましょう。DATの中身をメモ帳とかに貼って.datの拡張子で（分からなかったら.txtでもよい）保存するのでも結構です。翼型ファイルの本質は、縦にダーっとならんだXY座標ですので、それがあれば大丈夫です。

上のサイトは時々アクセスできなくなったりするので、有名どころも載せておきます。
http://aerospace.illinois.edu/m-selig/ads/coord_database.html

余談ですがUIUCのM-Selig先生は翼型界の超大御所です。プロペラの人とか知っている人も多いでしょうSD7037とかはSelig先生とDonovan先生の開発です。

DAE-11 //// DAEシリーズのなかで最も太い翼型です。31、21と比べると空力的性能はちょっと微妙。。。
DAE-21 //// DAEシリーズのなかでよく使われる翼型です。太さと空力性能がちょうどいいです。
DAE-31 //// 言わずと知れた名作。細いですが空力的な性能が良いです。
DAE-41 //// 影に隠れてあまりぱっとしないですが、低揚力係数での抗力係数が小さく、翼端向けの翼型としてよい翼型です。
DAE-51 //// プロペラ向けの翼型です。低レイノルズ数で使いやすいです。細いですが。
FX76MP120 //// 高いレイノルズ数、高い揚力係数では、圧倒的な揚抗比を誇ります。1度解析して取り憑かれる人が多いですが、高すぎる揚力係数、凄まじい捻りモーメント係数など欠点も多いです。
FX76MP140　//// FX76MP120のとんでもない性能には劣りますが、その代わり翼型の中ではかなり太い部類に入ります。それでも十分に良い性能です。
FX76MP160 //// 太さがとんでもないです。太すぎて二段失速の傾向が見えていますが、性能はよいです。
SD7037 //// うちのチームでプロペラに使っていました。下面が直線で作りやすいです。性能は普通です。探せば風洞試験データが手に入るようです。
epper66 //// 低レイノルズ数での揚抗比がとんでもないです。ちょっと細いですが使いやすいです。

保存した翼型はXGAG内の「FOILS」に保存しておくと、後々選択が楽です。このデフォルトの翼型フォルダは、プログラム内で変更することができます。

ではソフトを起動しましょう。

起動すると以下のような画面になります。これがターミナルです。（画像は開発中のものです。レイアウト等は今後変わることがあります）

初回起動したらOptionタブからDefault Directory & Foilsを選択して翼型ディレクトリとデフォルト翼型を設定して下さい。次回起動時からその翼型が選択された状態で起動します。
（初回起動時でもFileタブでNewを選ぶと、デフォルトがすべてBase Airfoilに設定されます）

初めデフォルトに設定されているXGAG内のFOILSディレクトリには、サンプルとしていくつか翼型を入れてあります。このまま翼型ディレクトリとして使うのもいいと思います。

デフォルト設定が終わったら早速解析の設定を行います。翼型はデフォルトから変えることも簡単です。.datもしくは.txtの翼型座標ファイルを選択するだけです。

まずは目標値を設定します。

解析する迎角α(deg) : XFOILで解析する迎角を設定します。

Reynolds数 : 解析するレイノルズ数です。初めのうちは高レイノルズ数=500000、低レイノルズ数=200000としておきましょう。定義は $Re = Chord * Velocity / \nu$ $\nu$ : 動粘係数（空気だと1.512×10^(-5)くらい）

揚力係数 : 目標となる揚力係数を設定します。揚力係数とは、単位面積の板が発生する揚力が、その板を風に垂直に立てた時の何倍になるかを表す係数です。

翼厚 : 翼の厚さが翼弦長に対して何%になるかを表しています。13〜15%くらいは太い翼型、8〜9%くらいが細身の翼型と呼ばれます。

翼厚計算位置 : 翼の厚さをどこで計算するか指定します。翼にかかる捻りモーメントができるだけ0になるように設定します。だいたい35〜40%が多いです。（別の考え方もあります）

最小抗力係数 : 計算で現れる異常値を弾くために設定します。抗力係数がこれ以下になると計算しません。とりあえずデフォルト値で大丈夫です。

続いて、評価関数について。評価関数は解析に書けた翼型がどれだけ目標に近づいているかの”点数”を表しています。高ければ高いほどその評価関数の中でよい翼型です。
抗力を小さくする場合は1/Cdの係数を大きく、捻りモーメントを小さくするには1/Cmの係数を大きくして下さい。となりの指数関数項は目標値から離れれば離れるほど小さくなります。いわば”罰則”です。
きつい制限を加えたければ罰則の項を大きくして下さい。大きくしすぎると抗力やモーメントの低減効果も小さくなってしまいます。

これで設定が済んだので、早速実行してみましょう！

XGAGの使い方-実行編

実行は　startボタンです。途中で止めるのはstopでできます。restartボタンはその設定のまま初めから計算し直します。

ある程度計算をすすめると、下のような画面になります。解析をすすめて、これはとっておきたい！という翼型ができたら、エクスポート（後述）するか、計算のセーブをしましょう。Ctrl+SもしくはファイルタブのSave,Save as で計算状況が.gagで保存されます（.gagの中身は基本的にCSVですが、エクセルで開いてしまうと壊れるので独自拡張子にしました。メモ帳で中身が見られます)

さて、計算をストップすると、以下のようなことができるようになります。

新規プロジェクトの開始 (Fileタブ→New)
プロジェクトのオープン (Fileタブ→Open)
翼型のエクスポート (翼型表示の下の世代選択ボタンとexprot foilボタン
翼型のロールバック

XGAGの使い方 -翼型のエクスポートについて

翼型のエクスポートは、解析を実行した全ての各世代で、最も良かった翼型を出力できます。出力したい翼型の世代を世代選択ボタンで選択して、export foilボタンを押して下さい。そうするとこんなウィンドウが出てきます。

このうち、速度分布平滑化とは

こんな圧力分布（速度分布）を

こんな風にすることです。

翼型の速度分布の振動は、翼型座標の２階微分値の不連続によって現れます。様々な２階微分値をもつ翼型を混ぜているのですから、速度分布の振動がでるのはある意味当然です。（元となる翼型の座標点ファイルの点数を増やしたり、翼弦長の正規化等をすれば多少は改善します。）速度分布振動をXFOILの機能であるMDES（複素速度マッピング)のFILTER機能を用いてある程度取り除くかどうか聞いているのが先のウインドウです。

上の段落よんで、訳分からん、という人はとりあえず平滑化かけておけばいいでしょう。ただし、その場合指定した揚力係数や翼厚からは少々ずれます。厳密に指定値を使いたい設計者さんとかは、平滑化せずにXFLR5等を用いて自分で速度分布いじってください。ついでに、その際のヒントは
「揚力係数の値は圧力係数Cpを用いて
$CL = \int_0^1Cpdx * cos(\alpha)$
（上面下面両方積分して下さい。つまりは翼型座標のxを横軸にCpを縦軸にとったグラフの囲む面積に迎え角のコサインをかけたものです）
このことは基本中の基本なのに意外と知られていないです。僕も翼型始めるまで知りませんでした。

要は面積一定になるように速度分布いじれって事です。

XGAGの使い方 -翼型のリバートについて

遺伝的アルゴリズムはもともと非効率な最適化手法です。少しでも効率を良くするためにXGAGは以下のような手法を取り入れています。

ルーレットルール　(評価の低い翼型も生き残る可能性を残すことによって系の多様性を維持する)
突然変異　(系の多様性を維持)
エリート戦略

このうちエリート戦略とは今までで最も良い個体の遺伝子を毎世代必ず1つ投入することによって、系の収束を早める手法の事をいいます。すごく重要なのですがこれがちょっと困ったちゃんなのです。
つまり、解析に異常値がでて最良個体の更新が行われてしまった場合あまり良くない遺伝子が毎世代投入されることになり、収束に悪影響を与える、ということです。XFOILは時々とんでもないCDを吐いたりするので、これを防ぐ必要があります。

対策として、最小抗力の設定を行っていますが、これだけでは不十分です。最適化をかけて放っておいて、異常値がトップに来てるのは哀しいものがあります。

そこで、リバートを使います。この機能は、セーブされている最良個体の遺伝子を、何世代か前のものに戻すことができます。
これはおかしい、と感じたら、履歴を見ながら世代選択ボタンで世代を選択し、その世代まで最良個体を戻すのが良いでしょう。変わるのは最良個体一つだけです。

終わりに

操作説明はこんなところでお終いです。何か不具合やバグ等見つけたら、コメント等で知らせていただけると幸いです。

アプリケーションを開発した身としては普及して、もっともっと自由に、飛行機創りを楽しんでくれるといいなと思います。

かつて、DAE Seriesというユートピアがありました。それの使用に理由は必要なく、賛美の歌を歌っていればよかった。
このXGAGが新しいユートピアを作り出してくれることを望みます。
そして、もしユートピアが完成したら、そのディストピアで、自由を求める主人公たちが現れることを、切に望んでいます。

つまらないギャグで笑うなんて、そんなのジョークだろ、と笑い飛ばしてもらえることを。

2013-09-13

多角形化行列とTR-797改

[人力飛行機設計]ご無沙汰しております。yuukivelです。

SSSで遺伝的アルゴリズムを用いた翼型設計について講演をすることになりました。
これからこのブログでも書きますが、私の翼型設計に関する知識の今のところ集大成になります。

閑話休題。本題に行きましょうか。

さて、表題の通り、今回は多角形化行列（勝手に命名）について解説する。この行列を作成することによって、多くのチームで行われている、循環分布の線形化を、行列演算だけで実現できるだけでなく、たとえば以前私が紹介した論文(TR-797,"非平面翼の最適設計-揚力と翼根曲げモーメントを与えた時の最小誘導抵抗-"）の誘導抗力最小化の演算を、多角形の循環分布で行うことができる。

つまり、こういうなめらかな循環分布ではなく…

下の図の上のようなカクカクした循環分布で、誘導抗力最小となる循環分布を求めてやろうということ。青が最大たわみ2.5mの構造制約（別の記事で詳しく述べます。 $\beta$ は使いません）　赤が構造制約無し。（ついでに下のグラフは誘導角度（青)と誘導速度の最小二乗法による多角形フィッティング(赤))

循環分布を多角形化すると何がいいか、っていうと、いろいろなものが線形化できること。
1つは捻り下げや翼型変化といった主翼空力の根幹を担う値を線形化することで、解析の時間を短縮したり、解析の精度を上げたりできる。
加えて、実際に飛行機を作ったときに、リブとリブの間のプランクやストリンガーは基本的には線形補完となる。つまり、製作の際にも設計再現性の向上を図ることができる。

いままで、私の知る限りでは、「楕円循環を作る」→「テーパー変化点を探す」→「点と点を直線で結ぶ」という、何というか、スマートじゃない方法で多角形化が行われてきたように思う。

そこで今回は、もっとスマートにやりましょう。

-多角形化行列の作り方
多角形化行列POLは、区切り位置の値ベクトル $\Gamma_{ysec}$ に左からかけることで、各翼素の値に”展開”される。すなわち
$\Gamma_y=POL*\Gamma_{ysec}$
ただし、 $\Gamma_y$ は各翼素の循環の値である。

基本となる式は、線形補完の式
$\Gamma_y=(\Gamma_{ysec2}-\Gamma_{ysec1})/(ysec2-ysec1)*(y-ysec1)+\Gamma_{y1}$ 　
による。ysec1は多角形化の区切りの位置を意味する。

この式を $\Gamma_{ysec1}$ と $\Gamma_{ysec2}$ に関する項に分ける。
[tex:\Gamma_y=(-(y-ysec2)/(ysec2-ysec1))*\Gamma_{ysec1}+*1*\Gamma_{ysec2}]

こんな感じで行列表記すると
$\begin{bmatrix} \Gamma_y1 \\ \Gamma_y2 \\ \Gamma_y3 \\ \Gamma_y4 \\ \Gamma_y5 \\ \Gamma_y6 \end{bmatrix}=\begin{bmatrix} -(y1-ysec2)/(ysec2-ysec1)&(y1-ysec1)/(ysec2-ysec1)&0 \\ -(y2-ysec2)/(ysec2-ysec1)&(y2-ysec1)/(ysec2-ysec1)&0 \\ -(y3-ysec2)/(ysec2-ysec1)&(y3-ysec1)/(ysec2-ysec1)&0 \\ -(y4-ysec2)/(ysec2-ysec1)&(y4-ysec1)/(ysec2-ysec1)&0 \\ 0&-(y5-ysec3)/(ysec3-ysec2)&(y5-ysec2)/(ysec3-ysec2) \\ 0&-(y6-ysec3)/(ysec3-ysec2)&(y6-ysec2)/(ysec3-ysec2) \end{bmatrix}*\begin{bmatrix} \Gamma_ysec1 \\ \Gamma_ysec2 \\ \Gamma_ysec3 \end{bmatrix}$

おわかり戴けただろうか。上の例ではy4とy5の間にセクションの切れ目がある。つまり、ここでカクッと曲がる。
これを私は勝手に多角形化行列と名付けた。

octaveで書いた多角形化行列の作り方を載せておく。num_sectionはセクションの総数を表し、yはそのままy座標、Ndiv_secはセクションの切れ目のy座標が変数yの何行目とするかを示している。

%-----循環分布多角形化行列
for i=1:Ndiv_sec(1)
	pol_g(i,1)=1;
	pol_g(i,num_section)=0;
end

for j=2:num_section
	for i=Ndiv_sec(j-1)+1:Ndiv_sec(j)
		pol_g(i,j-1)=(-(y(i)-y_section(j))/(y_section(j)-y_section(j-1)));
		pol_g(i,j)=((y(i)-y_section(j-1))/(y_section(j)-y_section(j-1)));
	end
end

-多角形化行列の応用
例として、先に書いたとおり多角形化行列を用いてTR797を改良する。
TR797の中には、ラグランジュの未定乗数法を用いて最小誘導抗力を求めるのに、
評価関数
$H({g};\mu_1 \mu_2)={g}^{T}*{A}*{g}-\mu_1*({c}^{T}*{g}-1)-\mu_2*({b}^{T}*{g}-\beta)$
の極値を求める。このとき、解を求めるのに逆行列問題
$A*{g}={a}$
に帰着するが、この $A$ の誘導抗力部分（ $Q_{ij}$ だった部分)に多角形化行列の転置行列を左から、多角形化行列を右からかけたものに変更すると、なんと、多角形の状態で最適化計算ができるのだ！

つまりはこういうこと
$\begin{bmatrix} \begin{array}\\ &Q+Q^{T}& \\ \end{array} & \begin{array} -c_1&-b_1\\ \vdots&\vdots \\-c_n&-b_n \end{array} \\ \begin{array} -c_1&\cdots&-c_n \\ -b_1&\cdots&-b_n \end{array} & \begin{array} 0&0\\0&0 \end{array} \end{bmatrix}*\begin{bmatrix} g_1\\ \\\vdots\\ \\g_n \\\mu_1 \\\mu_2 \end{bmatrix}=\begin{bmatrix} 0\\ \\\vdots\\ \\0 \\1 \\\beta \end{bmatrix}$

となる逆行列問題の各行列に以下のように多角形化行列をかけてやる。mはセクションの総数であり、行列の大きさはグッと減っている（セクションの総数＋2の行数になる)
$\begin{bmatrix} \begin{array}\\ &POL^{T}*(Q+Q^{T})*POL& \\ \end{array} & \begin{array} -C_1&-B_1\\ \vdots&\vdots \\-C_m&-B_m \end{array} \\ \begin{array} -C_1&\cdots&-C_m \\ -B_1&\cdots&-B_m \end{array} & \begin{array} 0&0\\0&0 \end{array} \end{bmatrix}*\begin{bmatrix} g_{sec1}\\ \\\vdots\\ \\g_{secm} \\\mu_1 \\\mu_2 \end{bmatrix}=\begin{bmatrix} 0\\ \\\vdots\\ \\0 \\1 \\\beta \end{bmatrix}$

ここで
$C=(\begin{array}c_1&\cdots&c_n \end{array})*POL$ (揚力制約の行列に右から多角形化行列をかける）
$B=(\begin{array}b_1&\cdots&b_n \end{array})*POL$ (構造制約の行列に右から多角形化行列をかける）
である。

このようにしてラグランジュの未定乗数法で解くことで、簡単に多角形で誘導抵抗最小な循環分布が探索可能なのだ。

さて、構造制約を外して、多角形化した場合としなかった場合の誘導抗力の大きさの違いを調べてみた。
もちろん、多角形化しなければ誰もが知っている「楕円循環分布」となる。その結果がこちら（スパン34.6m 機速7.5m/s 揚力90kgf)

見て分かるとおり、上が多角形化した楕円循環、真ん中が多角形化しない楕円循環、一番下が重ねてみたグラフ、である。

このときの誘導抗力は
多角形化したものが、5.9667N
多角形化しないものが、5.9599N
その差、0.0068695N

という結果になった。多角形化して最適化すれば、多角形化しない楕円循環分布とほとんど変わらない値となる。これで楕円循環分布より作りやすいとくれば、ぜんぜん損はない。

実はTR797の改良の真骨頂は、工夫した構造制約のところにあるのだけれども、今回は多角形化の記事にした。
この多角形化行列、もっともっと応用すれば、設計作業（最適化作業)はもっと楽になる。なんせ、変数がグッと減るのだから。

*1:y-ysec1)/(ysec2-ysec1

2013-06-02

粘性の世界に足を踏み入れてみた

翼型について

皆さんこんにちは。yuukivelです。

部屋の模様替えをして、買ってしまいました。人をダメにするソファ。

ビーズクッション CUBE L ブラウン

出版社/メーカー: フレックス販売
メディア: ホーム&キッチン
購入: 3人クリック: 14回
この商品を含むブログを見る

正確には無印良品のアレではないんですが、読書するのにちょうど良くて重宝してます。
カーペット敷いて一人暮らしの部屋の床に座れるようになると、なかなか捗る気がします。ベッドに寝転がりながらもできないし、机の前に座ってやるようでもないものをやるのにちょうどいいかなと。

今後のプログラム作成とかももっと捗ればいいなぁ。

さて、今回は前回の続きとなる。前回のポテンシャル流れ解析に粘性の影響を入れてみた。Navier-Stokesを離散化して、乱流モデルをつっこむような粘性考慮ではなく、いわゆるパネル法で求めたポテンシャル流れでの翼型周りの速度を境界層外縁速度として、境界層方程式を解いて翼型抗力を求める方法、である。用いた表面摩擦係数にかんする関係式は層流が「the Falkner-Skan one-parameter profile」、乱流が「Ludwieg-Tillmannの式」を用いている。乱流遷移位置判定法は通称「e^N法」を用いた。

NaGAf32_412_12.dat
上の翼型は前回記事にした遺伝的アルゴリズムで新たに作ってみた翼型。Clを1.2(理由は後述するが実は1.2を0.88で除した値)、翼厚を12%として外れたら罰則を加え、迎え角4° レイノルズ数が55万という条件の下Cdが最小になるように最適化をかけている。条件設定は、スパン20-22m位の速度競技型人力飛行機の翼根に用いる翼型を想定している。こちらも実はXFLR5のInverse Designで速度分布の平滑化を行っている。

XFOILはパネル法で翼型周りのポテンシャル流れの速度を求め、それを用いて排除厚を求めた後、その排除厚の分翼厚を増してまたパネル法で翼型周りのポテンシャル流れを解く、というプロセスを繰り返して、最終的に都合の良い粘性流れを求めているらしい。だからポテンシャル流れと粘性流れの速度分布に違いが生まれている。

今回作成したプログラムは収束計算を入れずに、前方よどみ点から粘性の影響を積分していって求めている。だからClはポテンシャル流れのものと同じになる。西山哲男先生の「翼型学」によると、粘性流れの揚力係数は、非粘性流れの揚力係数の約0.88倍になるという。結構大きく出る。このように収束計算を入れないで粘性の影響を求めるプログラムで有名なものとしてR eppler先生の「epplerコード」が挙げられる。

翼型学

作者: 西山哲男
出版社/メーカー: 日刊工業新聞
発売日: 1992/01
メディア: 単行本
この商品を含むブログを見る

抗力係数に関してはそこそこXFOILと一致している。といってもXFOILでのDAE31 迎え角5° レイノルズ数500000での解析の値が0.00921 私の書いたプログラムで0.00902という結果になった。この値を一致しているというかずれているというかは微妙なところだ。翼型によってもどれほど一致するかは変わってくるし、そもそもともに数値計算の出力した値であって実験値ではないので比較はできない。いい感じの値がでた、としか言えないのが残念なところ。

そして私が最も重要に感じたのは「遷移点判定」だ。今回はXFOILと同様のe^N法を用いた。e^N法とは、跡部隆先生の「e^N 法に基づく境界層の遷移予測とその検証」 http://airex.tksc.jaxa.jp/pl/dr/AA0001981000　によると、層流内における流れの擾乱の振幅が、初期振幅のe^N倍になったときに乱流に変化するというもの。このNの値には9とか10が一般に用いられているが、理論より導き出されたものではなく実験値らしい。つまり表面の状態が異なればこのNの値も変わってくるものと思われる。
このNという値だがXFLR5に触っている人はもう目にしたことがあるはずだ。これはまさしくdifine analysisのところに出てくるNCritのことだ。表面の状態によってNの値を変えることでより正確な遷移点予測ができるということである。

…とはいっても、Nの値を調整できる人は、特に人力飛行機界にかぎって言えばまずいないだろう。これに直面して感じるのは、「やっぱり人力飛行機開発には実験が全然足りないなあ」ということ。足りないというよりもこの部分には諦めが必要なのかな、と。数値計算は数値計算として割り切ってしまわなければやっぱりなにも変化していかなくなってしまうし、なんというか、「人力飛行機の」空力と付き合うのであれば、設計と現実の整合性を考えてもナンセンスな気さえしてきた。

ああ、うまくまとめられない。

で、今回はお約束（？）の通り、コードを公開してみようと思う。いつもの通りoctaveのプログラムで、GUIも結果のグラフ表示もない簡素なものだけれども、何かの役に立ってくれれば。
また、間違っているところとかあれば教えていただけると嬉しいです。
　
NaFoil.m

作業ディレクトリに翼型の座標ファイル(.datでも.txtでも可。ただしファイルの1行目は翼型名だとして読み込まないので注意)を入れてNaFoilを起動。翼型名「dae31.dat」みたいな感じで入力して条件を打ち込めばOK。

主に参考にした文献は
西山哲男「翼型学」日刊工業新聞社 1992年
Mark Drela and Michael B. Gilest,Viscous_Inviscid_Analysis_of_Transonic_and Low_Reynolds_Number_Airfoils
跡部隆,e(exp N)法に基づく境界層の遷移予測とその検証
河崎俊夫 ; 石田洋治,低マッハ数における翼型の翼型抗力の計算
石田洋治,圧縮性流れにおける翼型抗力の計算

他にも参考にした文献はたくさんあるが、主なものを抜粋して紹介した。
そもそも付け焼き刃な知識も多いので、いたらないところも多かった。ただなんというか、粘性はなんだかつかみ所がないようなところもあって、こりゃ学部低学年で踏み込むのは無理があるなぁとも感じた。

今回やってみて思ったことは、XFOILにしろなんにせよ、数値計算は数値計算だということ。遷移点の予測が狂えば揚抗比なんて乱高下するし、そもそもe^9倍になったとことで乱流に遷移するような表面が作れるとは限らない。外皮の薄さが少し変化しただけでも圧力分布はガタガタになるし、うーん、って感じ。設計だけするのであれば、粘性の数値計算に踏み込まない方が幸せかもしれない。ただし、これをやったおかげで翼型の圧力分布はだいぶ見ることができるようになったのと、翼型作れるようになったのは良かったかな。

これ以上粘性の世界に踏み込むのは、一生を粘性、もとい乱流に捧げてしまって泥沼にはまりそうな気がするし、とりあえずXFOILは素晴らしい、ってことが分かったし、やめておこうかなと思う。次は機体運動学的で最適化的なことに踏み込みたいなあ。

でもそろそろ自分で翼型も作れるようになったことだし、私が2年半の間勉強したことの総大成として、また人力距離競技機を描いてみたいなぁ

2013-05-08

2次元翼型まわりの流れ解析・翼型最適化　実践編

翼型について

引き続いて　皆さんこんにちは。yuukivelです。

前回の記事に引き続いて、2次元翼型まわりの流れ解析・翼型最適化　実践編　ということで、早速本題に入っていきましょう。

さて、今回は先にも述べてきたとおり、２次元翼型まわりのポテンシャル流れ（非粘性非圧縮)のCFDプログラムを作成した。
この手のソフトで有名なものとして、このブログで何度も出てきているXFLR5（中身のプログラムはXFOIL)が挙げられる。
XFLR5の有用性は揺るぎないものとなっているが、自分で書いた理由として

XFLR5では最適化問題における流れ解析に使いづらい。逆問題も限られた状態でしか解けない。

を挙げておこう。
つまりは

「自分のプログラムに組み込める流れ解析ソフトが欲しかった」

この理由につきる。

では、このプログラムの有用性を見てみよう。以下のグラフはDAE31（迎え角5°)のポテンシャル流れの圧力分布をXFLR5で解析したものと、自作プログラム「NaFoil」で解析した結果の比較だ。
青い菱形が自作ソフト。

非常に良く一致しているのが分かる。理論的には違ったものとなっているが（この記事の付録「理論編」で解説）感動するくらい一致している。
これが出てきたとき嬉しくて飛び上がりそうになった。

ここまでできると、次にやれそうなことが広がってくる。より実用的にするためには粘性の影響を含めなければならない。Navier Stokes方程式は壁が厚いので境界層方程式を解いて粘性の影響を含められないかと考えている。（NS方程式もそのうちやろうと考えている）
さらに、翼型の逆設計もこのプログラムによって可能になる。事実、重見先生の提唱する逆設計法はすぐにでも実装できると思われる。収束も早そうで、圧力分布を自在に操れるようになれば、かなり有用なツールとなろう。

で、このプログラムをポテンシャル流れの解析にしか使わないのはもったいないので、
「遺伝的アルゴリズムによる翼型設計」
やってみた。

遺伝的アルゴリズムを用いた空力設計の分野では東北大学の大林茂先生が熱心に取り組まれていて、今回の設計も大林先生のホームページにアップロードされている、
多目的遺伝的アルゴリズムによる空力最適設計
を参考にアルゴリズムを組んでみた。

その結果がこちら。

NaGAf54-512.dat
この翼型は迎え角5°のもと、DAE31,FX76MP120,FX76MP140,eppler66を基として遺伝的アルゴリズムにより翼厚12%縛りで揚力係数を最大化したものである。
XFLR5で解析した結果を示す。参考用にFX76MP120も合わせて解析したものも示す。レイノルズ数はともに50万である。太い緑が作った翼型「NaGAf54-512」

FX76MP120に比べて揚力係数が低く、抵抗が小さい。これのため最大揚抗比ではFX76MP120に匹敵し、さらに迎え角が小さい領域ではFX76MP120を凌駕している。
近年の日本の人力飛行機は使用する揚力係数の範囲が1前後と小さく、NaGAf54-512の方が使いやすいと思われる。さらにモーメント係数がFX76MP120より穏やかなため、正直なはなし、私なら自分の設計する人力飛行機に使う。DAE31に似た性能を持っているが、最大揚抗比ではこちらも凌駕している。

この翼型と似た性能を示し、揚力係数が大きい、NaGAf40-412という翼型も作りましてよ。

実は以上２つの翼型はXFLR5のInverse Design で速度分布を平滑化している。平滑化する前とほとんど性能は変わっていないのでご容赦いただきたい。

さて、遺伝的アルゴリズムをやってみて考えたことなど。。。
遺伝的アルゴリズムはその名前のかっこよさから、大学2年生がやってみたいと思ってしまう大二病の症状の１つと言われている（適当）
しかし、遺伝子の作り方や交配のさせ方によっては効率の悪い最適化になるどころか解にすらたどりつけない事もありうる。
計算回数が多くなりがちで私のVAIOが唸るなか待ち続けなければならない。

とは言っても、プログラムが簡単で、直感的でわかりやすいことを鑑みれば、最適化を本格的にやり始めるきっかけとしてはとてもいいのかもしれない。
交配のさせ方とか、考えるのが楽しい方法だと感じた。

- - - -

余談ではあるのですが、やってみて感じた事として、遺伝的アルゴリズムに関わるプログラミングは、「私情」が入ります（笑）
具体的に言えば、交配のさせ方は、収束のためと言い聞かせながらも、この交配のさせ方はなんかなーと思ったりします。同一世代で成績優秀個体がたくさん遺伝子を残せるとかね（笑）
他にも似たものどうし掛け合わせていくと全体として成績が下がっていったり、どっかの動物の事を連想してしまってなかなか楽しい（）プログラミングでした。

- - - -

ということで、実践編では作ったプログラムから作ったものを紹介してみた。機会があれば今回作ったプログラムは公開してみようかと思っている。特に日本の人力飛行機界隈では３次元翼の最適化は熱心に行われていても、２次元翼に関してはその取っつきにくさも相まってかほとんど手が付けられていないように思われる。今回の記事を通して、特にこれから人力飛行機に関わっていく人たちが、翼型設計とかに興味を持ってもらえれば嬉しい。

- - - - -

付録　2次元翼型まわりの流れ解析・翼型最適化　ざっと理論編

はてなダイアリーでは数式が書きづらいので、あんまり式を使わない感じで今回作成したソフト「NaFoil」「GA_Foil」について解説していきます。

NaFoil

ほとんど「流れの数値解析入門」内のアルゴリズムに従っています。
流れは非粘性非圧縮を仮定したポテンシャル流れですので、支配方程式は最も簡単な「ラプラス方程式」です。
これを数値的に解くため、翼型を多角形で近似し、各辺に分布渦を配置していきます。そして一様流と分布渦によって誘導される流れ（原理的にはビオ・サバールの積分です）を足し合わせて翼型線素に対する法線方向の速度がゼロとなる条件と翼型後縁から流れがなめらかに流れ去るクッタの条件を連立して分布渦強度に関する連立方程式を作り、逆行列問題を解いて分布渦強度を決定します。分布渦強度が決定されれば線素接線方向の速度はすぐに分かります。（この方法では分布渦強度にマイナス付けたものが接線方法速度になります）たったこれだけです。プログラムの行数は60行弱ととっても簡単です。

XFOILは翼線素に渦と湧き出しを分布させて法線方向速度がゼロになる流れを実現しています。これに比べて今回の方法は少し簡単です。しかし、粘性の影響を考慮する上ではXFOILの取る方法の方が有利らしいです。

このプログラムでは、XFLR5のCl計算値に比べて、揚力係数が高く出ます。どちらが正しいのか正直分かりませんが、もしかしたらXFLR5の揚力係数の値が実際より小さくなっている可能性も否定できないように感じます。検証してみる価値があるのかもしれません。この結果は大阪府立大学の方が言っていた、wortmann系列は設計速度より低速で飛びにくくて、DAEの方が低速で飛びやすいという経験に合致します。DAEの方がNaFoilとXFLR5の揚力係数計算値の差が小さいのです。

話を戻して遺伝的アルゴリズムによる翼型設計プログラム「GA_Foil」に行きましょう。名前適当です（笑）

GA_Foil

まず翼型を４つ準備します。探索する翼型のY座標 $Y_{con}$ を、準備した翼型のY座標 $Y_i$ を用いて

$Y_{con}=a_1Y_1+a_2Y_2+a_3Y_3+a_4Y_4$ $-0.6<=a_i<=0.6$

とします。そして交配させる個体の遺伝子を $a_i$ の65%と35%に分割したものとし、計8個の遺伝子とします。

評価関数は

$F_{con}=C_L*exp(-100*|(thickness)-(wanted.thickness)|)$

です。expの部分は翼厚と目標翼厚の差が大きければ大きいほど小さくなり、罰則となります。この評価関数を最大化するように解を探索します。

交配方法に移ります。最も良い値をたたき出した個体を保存し、これを越える個体が生じるまで保存します。さらにこの個体のクローンを毎世代投入しています。

そして上位より2個体ずつ一定の確率で遺伝子を交差させます。また、上位2個体は下位個体に対して一定確率で遺伝子を残します。
0.5%の確率で各個体に突然変異が起こります。

以上繰り返して最適解を探索していきます。

試行回数が一定世代となったところ（私が待ちきれなくなったところ）でプログラムを終了させ、最も評価関数の高かった個体を解として出力します。

ざっと説明するとこんなところです。長くなってきてしまったのでそろそろこの記事を締めくくりたいと思います。

最後に決めぜりふをば。

なぜ最適化をするのか。そこに最適解という頂があるからだ

(ﾄﾞﾔ

2013-05-08

2次元翼型まわりの流れ解析・翼型最適化　参考文献編

翼型について

こんにちは。yuukivelです。

最近、色々あった理由からとうとう麻雀に手を出し、どっぷり浸かっています。今日の記事の内容、2次元翼型まわりの流れ解析・逆設計・翼型最適化と相まって僕の時間がどんどん食われている状況です（笑)
遺伝的アルゴリズム（GA:Genetic Algorithm）の最適化を回しながらスマートフォンで牌を切る日々が続いております。（といってもここ数日ですが）

閑話休題。さて、本題に行きましょうか。

今回の内容はタイトル、そして前置きの中にあったとおり、2次元翼型まわりの流れ解析・逆設計・翼型最適化である。
イメージでイメージ（画像）を乗せればこんなやつ。この図は翼型の接線方向の速度をグラフに示したもの。ポテンシャル流れ（後述）を仮定している。

そもそも2次元翼型に手を出したのは、純粋な興味8割、人力飛行機への応用の必要性2割だった。このうち後者を具体的に述べれば、人力距離競技機における低迎角で揚抗比のよい翼型の必要を感じたからだ。

というわけで、確か3月中頃くらいからだっただろうか、2次元翼の設計や解析について論文や参考文献をあさっていた。最初のうちは英語論文ばかりで、とっても読むのが大変だった。しかも某国の航空宇宙学会の論文は有料なものが多く、手に入れるのに苦労した論文もあった。このとき見つけた英語論文で参考になるものをいくつか挙げていこう。

まずはUIUCのMichael S. Selig先生の論文。（英語が読めれば）かなり参考になる論文がUIUCの応用航空力学グループのホームページにアップロードされている。(http://www.ae.illinois.edu/m-selig/publications.html)
その中でも、翼型の逆設計（圧力分布を決めて翼型の形を決める方法）の論文として
"Multipoint Inverse Airfoil Design Method Based on Conformal Mapping" PDF

"Generalized Multipoint Inverse Airfoil Design" PDF

の2つがわかりやすい論文。この論文の中身は"Conformal Mapping"といって円柱周りの速度分布を翼型周りの流れに写像して任意の速度分布の翼型を得る方法。Selig先生が提唱する方法では、ある部分はこの迎え角でこの速度分布、またある部分は別の部分でまた違った速度分布を、といった形で速度分布を指定できるため、自在な翼設計ができる。
…といっても圧力分布を自在に扱えなければ意味がないのだが…その実力は僕はまだ持っていません（泣）

で、圧力分布に関する有名な論文を１つ。プロペラ設計法の１つ、Adkins&Liebeck法で有名なRobert H. Liebeck先生の
"A Class of Airfoils Designed for High Lift in Incompressible Flow" (有料論文)
が挙げられよう。

こんなところを（苦労して）読みながら、なんとか日本語論文はないかと探っていたところ、あった。今回作成した２次元翼型解析ソフトは以下２つの論文を参考とした。
"多翼素翼型の逆問題の解法" 重見　仁　http://repository.tksc.jaxa.jp/pl/dr/NALTR0571000
"低速で高揚力を得るための多翼素翼型の設計" 重見　仁 http://repository.tksc.jaxa.jp/pl/dr/NALTR0631000

そしてこの本があって２次元ポテンシャル流れのCFDは完成した。
「流れの数値解析入門」水野明哲　朝倉書店

記事が長くなりそうなので、二分割して書こうかと思う。続きは、2次元翼型まわりの流れ解析・翼型最適化　実践編　で！