OpenVSPのVSPAEROにおける有害抗力CDoの計算方法について説明する
はじめに
OpenVSPにはそもそも有害抗力（Parasite Drag）を推算するParasite Drag Toolというものが存在し、モデルの濡れ面積などの形状から、様々な推算式をもとに有害抗力CDoの値を計算することができる
一方、OpenVSPに組み込まれたVLMの解析ツールであるVSPAEROの結果を見てみると、何やら有害抗力らしきもの（CDo）が計算されている
（本来VLMはポテンシャル流れの解析ツールなので粘性による抗力は計算できない）
Beta Mach AoA Re/1e6 CL CDo CDi CDtot CDt CDtot_t CS L/D E CFx CFy CFz CMx CMy CMz CMl CMm CMn FOpt
0.000000000000 0.000000000000 -4.000000000000 1.900000000000 1.311552353449 0.114656540231 -0.524407559934 -0.409751019703 0.004250846583 0.118907386814 -0.848423225527 -3.200851957366 -0.077700800475 -0.317263619506 -0.848423225527 1.336940264082 0.030034010109 -6.043119578288 0.087434721292 -0.030034010109 -6.043119578288 -0.087434721292 0.000000000000
0.000000000000 0.000000000000 -2.000000000000 1.900000000000 1.485124859876 0.112280106431 -0.230038998600 -0.117758892169 0.004031172310 0.116311278741 -0.276023624493 -12.611572956579 -0.346662207559 -0.065857046484 -0.276023624493 1.488329888007 0.051486116571 -6.219237031168 0.165293160878 -0.051486116571 -6.219237031168 -0.165293160878 0.000000000000
0.000000000000 0.000000000000 0.000000000000 1.900000000000 1.360375720090 0.103214489196 -0.361169432919 -0.257954943723 0.006611155320 0.109825644516 -0.233660385215 -5.273695089760 -0.132784738418 -0.257954943723 -0.233660385215 1.360375720090 0.040341288894 -4.974911348263 0.167144015762 -0.040341288894 -4.974911348263 -0.167144015762 0.000000000000
0.000000000000 0.000000000000 2.000000000000 1.900000000000 1.810798641752 0.125885128726 -0.127459658082 -0.001574529356 0.012237021468 0.138122150193 -0.287071272037 -1150.057085222452 -38.544625620424 -0.064769531422 -0.287071272037 1.809640601864 0.059235651774 -6.463905493673 0.208365094382 -0.059235651774 -6.463905493673 -0.208365094382 0.000000000000
0.000000000000 0.000000000000 4.000000000000 1.900000000000 1.749575111588 0.243476703343 -0.141744347245 0.101732356097 0.026751563862 0.270228267205 -2.438167441169 17.197823570614 0.556904022954 -0.020559649144 -2.438167441169 1.752409724976 0.053699456754 -4.694084430460 0.227586895443 -0.053699456754 -4.694084430460 -0.227586895443 0.000000000000
0.000000000000 0.000000000000 6.000000000000 1.900000000000 1.651913327009 0.085525892024 -0.077812488415 0.007713403609 0.032862730558 0.118388622582 -1.135235889949 214.161401468523 6.547910444034 -0.165000812747 -1.135235889949 1.643670243187 0.061782976793 -4.028671753943 0.209555362115 -0.061782976793 -4.028671753943 -0.209555362115 0.000000000000
0.000000000000 0.000000000000 8.000000000000 1.900000000000 1.594771960539 0.067143407048 0.009178557182 0.076321964229 0.042904598802 0.110048005849 -0.330970870809 20.895321243925 0.616768119917 -0.146370154952 -0.330970870809 1.589873713880 0.071814936164 -3.499405534240 0.226868480159 -0.071814936164 -3.499405534240 -0.226868480159 0.000000000000
0.000000000000 0.000000000000 10.000000000000 1.900000000000 1.822344295509 0.058323289039 -0.138135923714 -0.079812634675 0.058296897111 0.116620186150 -0.812639491518 -22.832779583411 -0.770129021561 -0.395046867413 -0.812639491518 1.780799472309 0.093955405426 -3.234517064362 0.260759726389 -0.093955405426 -3.234517064362 -0.260759726389 0.000000000000
0.000000000000 0.000000000000 12.000000000000 1.900000000000 2.052305370069 0.096457861911 -0.190078573594 -0.093620711683 0.074232440203 0.170690302114 -0.950820885410 -21.921488666043 -0.832695655229 -0.518273154077 -0.950820885410 1.987992733244 0.054538209829 -3.096267604677 0.210323263635 -0.054538209829 -3.096267604677 -0.210323263635 0.000000000000今回はこれがどのように計算されているのかについて調べてみる
VSPAEROにおける有害抗力の計算方法
OpenVSPのVSPAEROの有害抗力は以下の4ステップで計算されている
- 平板と仮定して局所レイノルズ数から摩擦抗力を計算する
- M>0.6のとき、圧力抗力に対する圧縮性の補正係数を計算する
- NACA0012とNACA2412のデータから算出した圧力抗力（断面揚力による抗力）を計算する
- 1.~3.を足し合わせて有害抗力とする
つまり、数式にすると以下のようになる
\begin{align}
C_{D_{o}}&=C_{D_{f}}+f_{mach}C_{D_{p}} \\\\
C_{D_{f}}&=\frac{1.037}{\left(\log_{10}{Re}\right)^{2.58}} \\\\
f_{mach} &= \left\{\begin{array}{l}
1 & (M \le 0.6)\\
1+\frac{\left(\min\left(M, 1\right)-0.6\right)^2}{2} & (M>0.6)\\
\end{array}\right. \\\\
C_{D_{p}}&=\left(0.00625+0.01 |C_{l_{o}}|\right) \times \left(C_{l}-C_{l_{o}}\right)^2
\end{align}
1つずつ説明していく
↓参考
↓実装部分（17100行目あたりから）
平板の摩擦抗力
LocalVel = SpanLoadData(i).Span_Local_Velocity(k)[3];
// Chord
Length = SpanLoadData(i).Span_Chord(k);
// Local Re number, note that Local_Vel is scaled by Vref_
Re = MAX(1.e3,ReCref_ * LocalVel * Length / Cref_);
Cf = 1.037 / pow(log10(Re),2.58);
// Calculate flat plate viscous force... note density is booked kept as '1'
ViscousForce = 0.5 * Cf * pow(LocalVel * Vref_, 2.) * SpanLoadData(i).Span_Area(k);まず、平板（Flat Plate）を仮定して、二次元翼の上下面あわせた摩擦抗力係数を以下の式で計算する
\begin{align}
C_{D_{f}}&=\frac{1.037}{\left(\log_{10}{Re}\right)^{2.58}}
\end{align}
これは、完全に発達した乱流（fully developed turbulence）に対する平板の片面の摩擦抗力係数についてのPrandtl-Schlichting (1932) の式の約2倍となっている
\begin{align}
C_{D_{f}}&=\frac{1.037}{\left(\log_{10}{Re}\right)^{2.58}}
\simeq \frac{0.455}{\left(\log_{10}{Re}\right)^{2.58}} \times 2 \times 1.14
\end{align}
平板の片面の摩擦抗力係数であるPrandtl-Schlichting (1932) の式を上下面のために2倍するのは分かるが、謎係数1.14がどこからきているのかはよく分からない
まあ、干渉抗力とかなのだろう
圧力抗力に対する圧縮性の補正係数
LocalMach = 0.;
if ( Mach_ < 1. ) {
if ( Machref_ > 0. ) {
LocalMach = Machref_*ABS(LocalVel);
}
else {
LocalMach = Mach_*ABS(LocalVel);
}
LocalMach = MIN(LocalMach,0.95);
}
Fact = 1.;
if ( LocalMach >= 0.6 ) Fact = 1. + pow(MIN(LocalMach,1.) - 0.6,2.)/2.;次に、M＞0.6のとき、圧縮性についての補正係数を以下の式で計算する
\begin{align}
f_{mach} &= \left\{\begin{array}{l}
1 & (M \le 0.6)\\
1+\frac{\left(\min\left(M, 1\right)-0.6\right)^2}{2} & (M>0.6)\\
\end{array}\right.
\end{align}
この式がどこからきているのかは分からないが、趣味の範囲で遷音速以上の機体を設計することは当面予定していないので、これ以上は立ち入らないことにする
圧力抗力（断面揚力による抗力）
Clo_2d_ = 0.;
// Calculate 2D drag due to lift, simple fit to NACA 0012 data
CLv = 0.;
// Curve fit, NACA 0012 and 2412 using numbers split between smooth and standard roughness data
CvCl = 0.00625 + 0.01*ABS(Clo_2d_);
ViscousForce += 0.5*Fact*CvCl*pow(CLv-Clo_2d_, 2.) * pow(LocalVel * Vref_, 2.) * SpanLoadData(i).Span_Area(k);最後に、圧力抗力（断面揚力による抗力）を以下の式で計算する
\begin{align}
C_{D_{p}}&=\left(0.00625+0.01 |C_{l_{o}}|\right) \times \left(C_{l}-C_{l_{o}}\right)^2
\end{align}
Clo_2d_はプログラムで実装されているものの、ゼロが代入されたっきり変化がないので、使われているのかどうかはよくわからない
また、上記の式は、NACA0012（とNACA2412）のデータから引っ張ってきているそうなので、実際にXFLR5で計算した結果と比較してみる
なんとなくRe=1,000,000~10,000,000の間に入りそうな雰囲気である
NACA0012の有害抗力との比較
以上を踏まえて、摩擦抗力と圧力抗力を合計したVSPAEROのCDoの計算結果と、NACA0012のXFLR5のCDの計算結果（M=0，前縁で強制遷移）を比較してみる
↓VSPAEROのCDoの計算式
\begin{align}
C_{D_{o}}&=\frac{1.037}{\left(\log_{10}{Re}\right)^{2.58}}+1.0 \times 0.00625 \times C_{l}^2
\end{align}
まあ、当たらずとも遠からずという感じなので、最初のざっくり推算では十分だと思う（干渉抗力的な謎係数1.14がかけられているのでXFLR5の結果よりも抗力が大きく出ている）
おわりに
OpenVSPのVSPAEROにおける有害抗力CDoの計算方法について説明した
CDoの計算は、平板の摩擦抗力係数とNACA0012の圧力抗力をベースとして計算されており、ざっくり推算ではそれなりにそれっぽい値が出てくるようになっている
設計が進んで、詳細な翼型データが出るようになれば、多少計算結果に補正をかけてもいいかもしれない
↓関連記事
コメント