自作プロペラ設計プログラムを公開してみた

こんにちは。yuukivelです。

皆さん、体調には気を付けましょう。
数日前、A型インフルエンザと診断されました

インフルエンザの薬もだいぶ進歩したようですね。3年前、修学旅行の前日にインフルエンザ診断を受けたときは3日ぐらいきつかったのに、今回は薬を服用して1日で熱は下がってしまいました。医学の進歩に感謝です。

さて、現在プロペラ設計の話題がわき上がっているようで、今回はそれに対応して昔作ったプロペラプログラムの一部を公開して見ようかと思う。

人力飛行機開発チームにおいてプロペラを開発する際に、最も重要なのは何か。作れることである。
プロペラを支える桁がプロペラ翼面の中に収まること
プロペラを設計してその結果を図面に示すこと
が重要だと私は考える。

（まぁ、桁径のインプットミスって盛大に入らなかったことが1度ありましたがね）

今回公開するプログラムは
「できる限り回転中心近くまでプロペラ翼面を持ってくるように翼型混合を行う」
「翼型混合の際には、その設計法における最適化を崩さないようにする」
ということを意識して作成した。設計法としてはLarrabee法にもとづいた方法で計算している。

言語はフリーの数値計算ソフト"octave"

下記のリンクにプログラムをまとめたものを公開した。octaveの作業ディレクトリをこのフォルダにして、「empeller」とoctaveに入力すればプログラムが起動する。
https://www.dropbox.com/s/e5yaaf0nkh486h4/Empeller_yuukivel.zip

プログラムの改変等は自己責任で自由に行って下さい。ただし再配布の際は一声下さい。
バグ等あったら一報もらえると非常に助かります。

まだまだプログラムは上手くないことを実感。今だったらもう少し上手くかけるかな。

まずはinput

%----ペラ設定値
B=2;						%ブレード数B
dr=0.01;					%ペラ分割幅[m]
R=input('Radious of Propeller: ');	%ペラ半径R[m]
rpm=input('rpm: ');			%ペラ回転数[rpm]
rpm0=rpm;
omega=rpm/60*2*pi;			%ペラ回転数[rad/s]
V=input('Velocity: ');			%機体速度[m/s]
rho=1.184;				%空気密度ρ[kg/m^3]
nu=1.540*10^(-5);			%動粘性係数[m^2/s]
lambda=V/omega/R;			%進行率
T=input('Thrust: ');			%推力[N]
cl1=input('Cl (at r=dr): ');		%翼根での揚力係数
cl2=input('Cl (at r=20dr): ');		%翼素番号20での揚力係数
cl25=input('Cl (at r=30dr): ');	　　　　%翼素番号30での揚力係数
cl3=input('Cl (at r=R): ');		%翼端での揚力係数

%-----桁情報の読込
s_length=1.450					%桁長さ
s_t=0.0094333;					%桁テーパー比
s_mind=4.7						%長手方向最小直径[mm]
s_margin=1.4					%桁余裕[mm]
s_slength=round((R-s_length)/dr)	%先端ステンレス棒翼素数(翼厚計算を無視する翼素） 
s_d=2						%ステンレス棒直径


%-----半径位置r(ベクトル)を求める
i=0;
do
	i=i+1;
	r(i,1)=dr*i;
until(dr*i>R)
r(i,1)=R-dr*(i-1);
r(rows(r),:)=[];

続いて各翼素の揚力係数を求める。翼素番号20〜30の間は三次元スプライン補間を行う。これにより翼弦長変化がなめらかになる。

Cl(1:20,1)=linspace(cl1,cl2,20)';
Cl(51:rows(r),1)=linspace(cl25,cl3,rows(r)-50)';
Cl(21:50,1)=spline([r(1:20,1);r(51:rows(r))],[Cl(1:20,1);Cl(51:rows(r))],r(21:50));

さて、次からが設計ルーチン。空力解析結果の読込などのアルゴリズムは複雑な上に紹介してもしょうがないので略。

%-----理論はlarrabeeの計算法を使用．計算法はadkins&liebeck法に近い

%----ζを仮定し，ζ探索ルーチン　翼型の形状抗力を無視し、誘導抗力のみの影響をみて推力を求める。
zeta=0.1;
cov_zeta=0;
i=1;
do
	vd=V*zeta;
	xi=r./R;
	x=omega.*r./V;
	sinphi=1./sqrt(1+x.^2);
	cosphi=x./sqrt(1+x.^2);
	ad=1./2.*zeta./(x.^2+1);
	f=B./2.*(sqrt(lambda.^2+1)./lambda).*(1-r./R);
	F=2./pi.*acos(exp(-f));
	ad=1./2.*(vd./V)./(x.^2+1);
	G=F.*x.^2./(1+x.^2);
	gamma=2.*pi.*V.*vd./B.*G./omega;
	dL=gamma.*rho.*sqrt(V^2+*1;

ここで何をやっているかというと「翼型の形状抗力を全く無視して、誘導抗力のみの影響を考えて推力を算出し、その値が大きければ推力を小さく、大きければ小さくなるように値を調整」している。
これによって何がいいかというと、どの翼型を用いようとも、無次元移流速度ζの値を決めることが出来る、ということ。これによって翼根の翼型混合が、誘導エネルギ損失最小となる循環分布を崩さないままに可能となる。

%-----翼厚判定
foila_thn=f_thn(foila,fsa,foilb,fsb,1,0.25);	
foilb_thn=f_thn(foila,fsa,foilb,fsb,0,0.25);
fillet_thn=f_thn(foilb,fsb,fillet,fsf,0,0.25);
	
for i=rows(r):-1:1
	foil_thn(i,1)=f_thn(foila,fsa,foilb,fsb,1,0.25)*chord(i,1)*1000;
	if(i>=rows(r)-s_slength)
		%-----s_mind[mm]ステンレス棒のところ
		spar_thn(i,1)=s_margin*2+s_d;
	else
		spar_thn(i,1)=s_margin*2+s_mind+s_t*(R*1000-dr*(s_slength)*1000-r(i,1)*1000);
	endif
	spar_thn2(i,1)=spar_thn(i,1);
	
	%-----桁が入っているかどうかの判定
	if(spar_thn(i,1)>foil_thn(i,1)&&i*2.^2).*chord/nu;

この部分がこのプログラムの切り札、翼厚判定。まず全てをfoilaとして翼厚を計算する。桁の径を計算しておいて、桁が入らなかったところをチェック。翼弦長はもう決まっているから、foilbとfoilaを桁が入る混合率になるように混ぜる。（これは代数演算で出来る) それでもダメだった場合はfilletとfoilbを混ぜる。ただしこの部分の空力的な影響はないとして計算を進める。

残りは、Larrabeeにおける推力とトルクの演算なので、ソースコードを参照されたい。
設計結果はpelaという関数に、長手方向位置(m)、翼弦長(m)、翼素取付角、翼型混合率、翼厚、の順で格納される。

今回は公開版と言うこともあって、プロペラの詳細な空力計算と、切り札その2であるソリッドワークス出力支援と、PostScript図面作成機能は削除させてもらった。
これらの公開は、気長に待って欲しいと思っている。

ぜひ、触って、遊んでみて、プロペラ設計のおもしろさが伝わってくれればなぁと思う。

最後にこのプログラムを使って設計したプロペラの画像を一枚添付して、この記事の終わりとしよう。

*1:1-ad).*r.*omega).^2).*cosphi*B; L=trapz(r,dL); %-----判定係数k k=1.02; if(L>T*k+0.2) if(L>T*k+3) zeta=zeta-0.01; else zeta=zeta-0.001; endif elseif(L3000) printf('zeta is not converted. zeta:%f\n',zeta); i cov_zeta=1; endif i++; until(cov_zeta==1) zeta0=zeta; %-----chordの決定 chord=4.*pi./B.*lambda.*G./sqrt(1+x.^2).*zeta.*R./Cl; phi=atan(lambda./(r./R).*(1+zeta/2

*2:spar_thn(i,1)/chord(i,1)/1000)-foilb_thn)/(foila_thn-foilb_thn); if(foil_mix(i,1)<0) t_judge(i,1)=-1; foil_mix(i,1)=-1+((spar_thn(i,1)/chord(i,1)/1000)-fillet_thn)/(foilb_thn-fillet_thn); %-----全く桁が入らないところはなかったことにする if(foil_mix(i,1)<-1) chord(i,:)=; foil_mix(i,:)=; foil_thn(i,:)=; x(i,:)=; sinphi(i,:)=; cosphi(i,:)=; ad(i,:)=; f(i,:)=; F(i,:)=; G(i,:)=; gamma(i,:)=; dL(i,:)=; phi(i,:)=; Cl(i,:)=; r(i,:)=; spar_thn(i,:)=; xi(i,:)=[]; t_judge(i,1)=-2; end end else t_judge(i,1)=0; foil_mix(i,1)=1; endif endfor for i=1:rows(r) if(foil_mix(i,1)>0) foil_thn(i,1)=f_thn(foila,fsa,foilb,fsb,foil_mix(i,1),0.25)*chord(i,1)*1000; else foil_thn(i,1)=f_thn(foilb,fsb,fillet,fsf,1+foil_mix(i,1),0.25)*chord(i,1)*1000; end endfor Re = sqrt(V^2+(omega.*r.*(1-ad