0 リーマンサットについて
趣味で宇宙開発を行う団体「リーマンサット・プロジェクト」がお送りする新春アドベントカレンダーです。
リーマンサット・プロジェクトは「普通の人が集まって宇宙開発しよう」を合言葉に活動をしている民間団体です。
他では経験できない「宇宙開発プロジェクト」に誰もが携わることができます。
興味を持たれた方は https://www.rymansat.com/join からお気軽にどうぞ。
1 この記事を書くに至った経緯
リーマンサットって宇宙の団体でしょ？何故航空機？と思った方もいらっしゃるかもしれませんが、火星を飛ぶ飛行機や、高高度を飛ぶドローンの開発なども様々な研究機関で進められており、航空機が宇宙開発の一つのソリューションになってくるのではないかと思い、アドベントカレンダーにこのトピックを投稿しました。
#2 はじめに
※機械系・物理系の学生が習う、非圧縮性流体力学の基礎を学習していることを前提にしてます。その辺が未習の方は、流体力学の教科書等をご参照ください
近年は計算機の発達により、CFD（数値流体力学）での流体解析が一般的になってきました。今回紹介する流体解析法は、ポテンシャル解法と呼ばれる、境界層や剥離の影響を考慮しない手法です。模型飛行機なんかの空力設計に使われたりしてます。注意してほしいのは、流体の粘性・圧縮性・剥離の影響を無視できる範囲内でのみこの解法が適用できます。したがって、粘性の高い低レイノルズ数流体や、Mach数の高い圧縮性流体、高迎角域での翼型の計算等には適用できません。
#3 参考にした文献
Katz,Plotkin Low Speed Aerodynamics
https://www.cambridge.org/core/books/lowspeed-aerodynamics/077FAF851C4582F1B7593809752C44AE
という航空力学の本を参考にしました。この文献は速度ポテンシャルの導出から3次元の流体解析まで網羅的に書かれており、かなり勉強になります。
#4 計算方法
#4.1 ポテンシャル解法の考え方
翼表面をN個の要素に分割し、分割した各要素上に二重湧き出しを配置して、その時の速度場を計算するやり方です。以下の図のように、翼表面上のi番目の微小要素に、二重湧き出しの要素$\mu_i$を配置します。この時の微小要素の速度ポテンシャルは、以下の式の様になります。（グリーンの定理などから導出することができます。詳しくは参考文献のchapter3を読んでください）（説明不足で申し訳ございません泣）
$$ \Phi^{*}_i(x,z) = \frac{1}{2π}\int[σ\ln r - µ \frac{\partial}{\partial n}(\ln r)]dS + \Phi \tag{1} $$
#4.2 ディリクレ条件の適用
4.1でもとめた微小要素での速度ポテンシャルですが、何か境界条件を加えないとこの式から微小速度を算出することができません。そこで、今回は境界条件として、「翼表面での速度ポテンシャルは一定」というディリクレ条件を適用します。
翼表面のポテンシャルは、翼表面に配置した二重湧き出しによる擾乱ポテンシャル$Φ_i$と、一様流のポテンシャル$Φ$を足し合わせたものなので、これが一定になる時、以下の式が成立します。
$$ \Phi_i = \Phi_\infty + \Phi_i = const. \tag{2} $$
ここで、4.1の考え方に則り、微小要素上の二重湧き出しの式を適用すると、
$$ Φ^{*}_i(x,z) = \frac{1}{2π}\int[σ\ln r - µ \frac{\partial}{\partial n}(\ln r)]dS + Φ = const. \tag{3} $$
という式になります。やや複雑な式なので、以下のように係数を定義して簡易化します。
$$ \frac{1}{2π}\int(\ln r)dS|_j \equiv B \tag{4} $$
$$ -\frac{1}{2π}\int \mu \frac{\partial}{\partial n}(\ln r)dS|_j \equiv C \tag{5} $$
これをN個の要素全体について足し合わせると、以下の式が成立します。
$$ \sum_{j=1}^NB_j\sigma_j + \sum_{j=1}^NC_j\mu_j + \Phi_\infty = const. \tag{6} $$
#4.3 Linear Doublet Method による計算
4.2では湧き出しも計算に入れてますが、これをさらに単純化して、二重湧き出しのみの分布を仮定します。Low Speed Aerodynamcisでは、linear double method と呼ばれており、以下の式が成立します。
$$ \sum_{j=1}^NC_j\mu_j + \Phi_\infty = 0 \tag{7}$$
ここで上式中の主流の速度ポテンシャル$\Phi$について、一様流のポテンシャルの式から、以下のように定義できます。
$$ \Phi_\infty = U_\infty x +W_\infty z \tag{8} $$
次に、係数$C_j$の算出方法についてですが、以下のようにして求めます。
4.1の図の翼表面上の座標系において、極座標系$(r,\theta)$における任意の点Pでの、二重湧き出しの影響によるポテンシャルは、
$$ \Phi = \frac{\mu_j}{2\pi}(\theta_{j+1} - \theta_{j}) + \frac{1}{4\pi}\frac{\mu_{j+1}-\mu_j}{x_{j+1}-x_j}[2x(\theta_{j+1}-\theta_j) + z\ln \frac{r_{j+1}}{r_j}] $$
となります。（これも導出が難しいので割愛します泣）
この式を、$\mu_j$と&\mu_{j+1}&についての項に分けます。すなわち、
$$ \Phi = \Phi_a +\Phi_b \tag{9}$$
$$ \Phi_a = \frac{\mu_j}{2\pi}(\theta_{j+1} - \theta_{j} + \frac{1}{x_{j+1}-x_j}[x(\theta_{j+1}-\theta_j) + \frac{z}{2}\ln \frac{r_{j+1}}{r_j}]) \tag{10}$$
$$ \Phi_b = \frac{1}{2\pi}\frac{-\mu_{j+1}}{x_{j+1}-x_j}[x(\theta_{j+1}-\theta_j) + \frac{z}{2}\ln \frac{r_{j+1}}{r_j}] \tag{11}$$
係数の計算についてですが、行列$C$の各要素$c_{ij}$について考えると、i番目の微小要素上の二重湧き出しによってj番目の要素に発生したポテンシャルを意味しているので、これを1~N番までの各要素で計算するので、$c_{ij}$は、
$$ c_{ij} = \Phi_{b(i,j-1)} + \Phi_{a(i,j)} \tag{12}$$
となります。これにより、となりあうj-1番目とj番目の二重湧き出しが及ぼす影響が計算できます。
このようにして各要素での$C_j$を求まったら、上の式から、二重湧き出しの強さ$\mu_j$を求めることができます。具体的には、$\sum_{j=1}^NC_j$は、i番目の要素の二重湧き出しがj番目の要素に及ぼす影響を表しているので、N+1×N+1の正方行列となります。
\sum_{j=1}^{N+1}C_j\mu_j = \left(\begin{array}{ccc}
c_{11} & \cdots & c_{1i} & \cdots & c_{1,N+1} \\
\vdots & \ddots & & & \vdots \\
c_{i1} & & c_{ii} & & c_{i,N+1} \\
\vdots & & & \ddots & \vdots \\
c_{N+1,1} & \cdots & c_{N+1,i} & \cdots & c_{N+1,N+1} \\
\end{array} \right)
\left( \begin{array}{ccc}\mu_1 \\ \mu_2 \\ \vdots \\ \mu_{N+1} \end{array} \right)
=\left( \begin{array}{ccc}\Phi_1 \\ \Phi_2 \\ \vdots \\ \Phi_{N+1} \end{array} \right) \tag{13}
Cの行列の各要素を計算した後、逆行列を解いて、二重湧き出しの強さ$\mu_j$を求めます。数値計算する際に要素数が大きくなる時は、SOR法などの反復解法を使うと良いです。
#4.4 圧力分布の算出
二重湧き出しの強さ$\mu_j$が求まったら、各要素上での、翼表面上での接線方向の速度が以下の式から求まります。
$$ Q_{tj}=\frac{\mu_j-\mu_{j+1}}{x_{j+1}-x_j} \tag{14}$$
ここで求めた速度から、各要素上の圧力係数が求まります。
$$ Cp_j = 1 - \frac{Q^{2}_{tj}}{U^{2}} \tag{15}$$
これを翼の各要素の座標との対応をとることで、圧力分布を解析することができます。
この圧力分布から、翼の性能を示す指標である、揚力係数を求めることができます。
#5 Matlabによる実装
MATLAB R2017aで作成してみました。ToolBoxとかは特に使ってません。
NACA0012翼型の座標データを持ってきてdatファイルなどで読み込んでおけば使えます。
%%%%%2次元パネル法による翼の空力計算
%%%%%二重湧き出しを配置
%%%%%Low Speed Aerodynaimcsの11.5.2を参照
%N1=102;
%%%%%%変数の初期化
A=zeros(N1,N1);CP=zeros(N1-1,1);RHS=zeros(N1,1);dCL=zeros(N1-1,1);dCD=zeros(N1-1,1);dCM=zeros(N1-1,1);
EP=zeros(N1-1,2);PT1=zeros(N1-1,2);PT2=zeros(N1-1,2);CO=zeros(N1,2);
Theta=zeros(N1,1);DL=zeros(N1,1);
PI=3.14159265;
%%%%解析条件
aoa=aoa/57.2958;%rad(迎え角)
%c=1.0;%m 翼弦長
%U=1.0;%(m/s) 主流速度
M=N1-2;%%%%%%パネルの要素数。datファイルの要素数-1に合わせること！！
N=M+1;
%%%%%%%%翼型datファイルの読み込み
Airfoil=dlmread('NACA0012.dat');
%%%%%%%プロットを時計まわりに変換。翼弦長倍する
EP1=Airfoil(:,1)*c;
EP2=Airfoil(:,2)*c;
for i=1:N
EP(i,1)=EP1(N-i+1);
EP(i,2)=EP2(N-i+1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%パネル終端点の設定
for i=1:M
PT1(i,1)=EP(i,1);
PT2(i,1)=EP(i+1,1);
PT1(i,2)=EP(i,2);
PT2(i,2)=EP(i+1,2);
end
%%%%%%%%%%
%%%%%%パネル角度
for i=1:M
DZ=PT2(i,2)-PT1(i,2);
DX=PT2(i,1)-PT1(i,1);
Theta(i)=atan2(DZ,DX);
%Theta2(i)=Theta(i)*180/PI;
end
%%%%%%%%
%%%%%%パネル参照点の設定。（各パネルの中央点）
%%%%%%%%%%%%%
for i=1:M
CO(i,1)=(PT2(i,1)-PT1(i,1))/2 + PT1(i,1);
CO(i,2)=(PT2(i,2)-PT1(i,2))/2 + PT1(i,2);
end
%%%%%%%影響係数の算出
for i=1:M
for j=1:M
%%%%%%%%%%%%%パネル上の局所座標系への変換
XT=CO(i,1)-PT1(j,1);
ZT=CO(i,2)-PT1(j,2);
X2T=PT2(j,1)-PT1(j,1);
Z2T=PT2(j,2)-PT1(j,2);
%%%%%%%%%%%%%
X=XT*cos(Theta(j))+ZT*sin(Theta(j));
Z=-XT*sin(Theta(j))+ZT*cos(Theta(j));
X2=X2T*cos(Theta(j))+Z2T*sin(Theta(j)); %%%%%%パネルの長さ
Z2=0;
%%%%%パネル長さの保存
if(i == 1)
DL(j) = X2;
end
%%%%%%
%%%%%%影響係数算出のための距離や角度の算出
R1=sqrt(X^2 + Z^2);
R2=sqrt((X-X2)^2 + Z^2);
TH1=atan2(Z,X);
TH2=atan2(Z,X-X2);
if(i == j)
PH1=-0.5*(X/X2 - 1);
PH2=0.5*(X/X2);
else
PH1=0.15916*(X/X2*(TH2-TH1)+Z/X2*log(R2/R1)-(TH2-TH1));
PH2=-0.15916*(X/X2*(TH2-TH1)+Z/X2*log(R2/R1));
end
if(j == 1)
A(i,1) = PH1;
HOLDA = PH2;
elseif(j == M)
A(i,M)=HOLDA+PH1;
A(i,M+1) = PH2;
else
A(i,j)=HOLDA+PH1;
HOLDA=PH2;
end
end
%%%%%wakeの影響A(M+2)
XW=CO(i,1)-PT2(M,1);
ZW=CO(i,2)-PT2(M,2);
DTHW=-atan(ZW/XW);
A(i,M+2)=-0.15916*DTHW;
RHS(i)=(CO(i,1)*cos(aoa) + CO(i,2)*sin(aoa));
%RHS(i)=U*(CO(i,1)*cos(aoa) + CO(i,2)*sin(aoa));
end
%%%%%%二重湧き出しの勾配条件（左側行列の一番下の行）
A(M+1,1)=-1;
A(M+1,2)=1;
A(M+1,M)=1;
A(M+1,M+1)=-1;
%%%%%%Kutta条件
A(M+2,1)=-1;
A(M+2,M+1)=1;
A(M+2,M+2)=-1;
N=M+3;
%%%%%%%行列の計算
G=A\RHS;
%%%%%二重湧き出しの強さを翼表面に沿った速度に変換
%%CP分布も算出
CL=0.0;CD=0.0;CM=0.0;
for i=1:M-1
R=(DL(i)+DL(i+1))/2;
T1=(G(i)+G(i+1))/2;
T2=(G(i+1)+G(i+2))/2;
VEL=(T2-T1)/R;
CP(i)=1-((U*cos(aoa+Theta(i))+VEL)^2/U^2);
%%%%%揚力係数の計算。各要素の揚力分布から足し集めて計算する
dCL(i)=-CP(i)*DL(i)*cos(Theta(i))/c;
dCM(i)=-dCL(i)*(CO(i,1)-0.25*c)*cos(aoa);
CL=CL+dCL(i);
CM=CM+dCM(i);
end
#6 あとがき
だいぶ駆け足になってしまい、説明が至らぬ点があるかもしれませんが、航空宇宙に興味がある方は、是非リーマンサットにお越しください！！
次回は次回は@prestissimoさんの「趣味の宇宙開発の広報のおしごと」です!