人力飛行機を割と本気で設計してみた。第一報
明けましておめでとうございます。
今年もよろしくお願いします。
去年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の急激な低下を起こさないようなワイヤー張力を探したり、とにかく自由に設計できてしまう。
恐ろしい。そしてこの手法は本当に楽しい…
では次の第二報の記事にて、翼型の設計とかの話をしましょ。
