🚲
【マサカリ求む】二次元の物体にかかる空気抵抗が完全に理解できたかもしれない話

2024/11/25に公開
 三行で既知の二次元図形にかかる空気抵抗を計算したい！！！
だから自分なりに方法を考えてみました！！！
結果はそこそこいい線行ってるんじゃないかとは思いますが実践で使えるかは正直微妙なところです！！！

 注釈この記事は三年ほど前に投降した記事のブラッシュアップ版となります。

 思想の話
 何故にそんなことを考えたの？当初の目的は、「空気抵抗を最小限にするにはどんな断面がいいの？」という疑問に対する答えでした。

そのためにこの記事を参考にジュコーフスキー翼を作ってみて、その抗力を調べようと思ったのですが、

「抗力係数が分からねぇ…」という根本的な問題にぶち当たりました(当然ですが)。

 で、それじゃあどうやって抗力を求めるの？そこで色々とGoogle先生に聞いてみるも、ほぼ答えは
D[N]=12ρv2SCD([kg/m3]⋅[(m/s)2]⋅[m2])
D[N] = \frac{ 1 }{ 2 } \rho{} v^2 S C_D ( [kg/m^3] \cdot{} [(m/s)^2] \cdot{} [m^2] )
D[N]=21​ρv2SCD​([kg/m3]⋅[(m/s)2]⋅[m2])( ρ\rho{}ρ は空気の密度、 vvv は速度、 SSS は前方投影面積、 CDC_DCD​ は抗力係数)

この記述以外がなかなか見当たらねぇ…。と悩んでいたそんなところ、こんなサイトがあるじゃありませんか。
https://h-a-labo.amebaownd.com/posts/2558036
(現在は撤回されていますが、僕がアクセスした当時は抗力の試算方法も載っていました)

こりゃすげぇ!!と早速試しに円柱の空気抵抗を計算してみたら、既存の式に比べて値が異様に高い…。

(従来の式で約48[N]に対して、上のサイトの方式(以降、「長谷川式」)の抗力は約5023[N])

(円柱は半径0.5[m]の高さ1[m]、標高0[m]で15[℃]の状態、風速は40[km/h]を想定しました)

上のサイトの管理者(以降、長谷川様)に問い合わせても、
この２、３日熟慮しましたが、抗力計算結果の間違いの理由が解明できておりません。
とのこと。

 じゃあ、詰んだの？んな訳はあるめぇ、他に方法があるはずだ、と思い色々な値を見てみました。

その紆余曲折は長い上に面白くないので割愛しますが、ある値に行きつきました。

 固有音響インピーダンス簡単に言えば音の伝わりづらさです。交流電圧とインピーダンスの関係と同じです。

長谷川様の仮説では、
物体が空気中を移動するとき、物体の前面には空気が圧縮されて高気圧域が、後面には空気が膨張して低気圧域が発生する。
ということになっています。これについては僕はとても正しいと思います。ですが同時に、
空気中に生じた気圧差は音速で解消されるのですが、物体の移動が続くと物体の表面に気圧差は発生し続けるのです。
とエビデンスなしに述べております。これには疑問符が付きました。

その点について長谷川様に問い合わせてみても、
小生は「気圧差の解消速度はその場の音の伝播速度」と考えています。もちろん流速でも移動速度でもありません。

「その場の音の伝播速度」とはその場の「音響インピーダンス」を考慮してのことで、上空3000mの音速を考慮してセスナの揚力を計算しています。
という回答で、いまいちピンと来ません。

音響インピーダンスは
pSv([Pa][m2]⋅[m/s])
\frac{ p }{ Sv } ( \frac{ [Pa] }{ [m^2] \cdot{} [m/s] } )
Svp​([m2]⋅[m/s][Pa]​)で求まるのですが、長谷川式揚力の計算式
L[N]=vc⋅sin⁡θ⋅p1atm⋅1atm⋅S⋅cos⁡θ([m/s][m/s]⋅[Pa][Pa]⋅[Pa]⋅[m2])
L[N] = \frac{ v }{ c } \cdot{} \sin{} \theta{} \cdot{} \frac{ p }{ 1atm } \cdot{} 1atm \cdot{} S \cdot{} \cos{} \theta{} ( \frac{ [m/s] }{ [m/s] } \cdot{} \frac{ [Pa] }{ [Pa] } \cdot{} [Pa] \cdot{} [m^2] )
L[N]=cv​⋅sinθ⋅1atmp​⋅1atm⋅S⋅cosθ([m/s][m/s]​⋅[Pa][Pa]​⋅[Pa]⋅[m2])( vvv :流速、 ccc :その場の音速、 ppp :雰囲気圧、 SSS :表面積、 θ\theta{}θ :表面の傾き)

のどこにも音響インピーダンスを考慮した部分が出てきていません。
「その場の音速」と「雰囲気圧」、で計算することでその場の音響インピーダンスを考慮したつもりでした。
とのことですが…。

さて、何故「音響インピーダンス」にたどり着いたかという話をしていませんでしたが、

気圧差を音速で解消しようとするけれども、「音響インピーダンス」に邪魔されて

実際には音速でない速度で解消されていると考えたからです。

(例えば、団扇で扇いでから実際に風が来るまで意外とタイムラグがあると思います)

話を戻して、数式に空気の固有音響インピーダンスを登場させましょう。

空気の固有音響インピーダンス(厳密には疎密波が正弦波の場合)は、
Z=ρc([kg/m3]⋅[m/s])=pava([Pa][m/s])
Z = \rho{} c ([kg/m^3] \cdot{} [m/s]) = \frac{ p_a }{ v_a } ( \frac{ [Pa] }{ [ m/s ] } )
Z=ρc([kg/m3]⋅[m/s])=va​pa​​([m/s][Pa]​)( ρ\rho{}ρ :その場の密度、 ccc :その場の音速、 pap_apa​ :その場の気圧、 vav_ava​ :その場の速度)

で求まるので、その場の速度 vav_ava​ を求めたければ
va=paρc
v_a = \frac{ p_a }{ \rho{} c }
va​=ρcpa​​となります。ところで、長谷川様の仮説では物体の前方で上昇した気圧(＝物体の後方で降下した気圧)は
v0csin⁡θ⋅p0sin⁡θ=v0p0sin⁡2θc
\frac{ v_0 }{ c } \sin{} \theta{} \cdot{} p_0 \sin{} \theta{} = \frac{ v_0 p_0 \sin{}^2 \theta{} }{ c }
cv0​​sinθ⋅p0​sinθ=cv0​p0​sin2θ​となるため、式をまとめると
va=v0p0sin⁡2θρc2
v_a = \frac{ v_0 p_0 \sin{}^2 \theta{} }{ \rho{} c^2 }
va​=ρc2v0​p0​sin2θ​となります。僕はこの vav_ava​ が物体によって圧縮された空気の基本的な流速だと判断しました。

 理論構築の話
 圧力に関する立式上で求めた流速からベルヌーイの定理を用いて圧力を導けば圧力に関する力が導けるはずです。

二次元における非圧縮のベルヌーイの定理では
pρ+v22=const.
\frac{ p }{ \rho{} } + \frac{ v^2 }{ 2 } = const.
ρp​+2v2​=const.なので、仮に流れの外側と内側が同じ流体で構成されていて、内側には一切の流速がなく、かつその境界が無限に薄い何かで覆われていると考えれば下の式が成り立ちます。
piρ=paρ+va22pi−paρ=va22ΔDp=ρ12(v0p0sin⁡2θρc2)2∣Δy∣=12ρv02(p0ρc2)2∣Δy∣5Δx2+Δy24⋯1◯
\begin{aligned}
\frac{ p_i }{ \rho{} } &= \frac{ p_a }{ \rho{} } + \frac{ v_a^2 }{ 2 } \\
\frac{ p_i - p_a }{ \rho{} } &= \frac{ v_a^2 }{ 2 } \\
\varDelta{ D_p } &= \rho{} \frac{ 1 }{ 2 }( \frac{ v_0 p_0 \sin{}^2 \theta{} }{ \rho{} c^2 } )^2 |\varDelta { y }| \\
&= \frac{ 1 }{ 2 } \rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2 \frac{ |\varDelta{ y|^5 } } {\sqrt{ \varDelta{ x^2 } + \varDelta{ y^2 } } ^4 } \cdots \textcircled{ 1 }
\end{aligned}
ρpi​​ρpi​−pa​​ΔDp​​=ρpa​​+2va2​​=2va2​​=ρ21​(ρc2v0​p0​sin2θ​)2∣Δy∣=21​ρv02​(ρc2p0​​)2Δx2+Δy2​4∣Δy∣5​⋯1◯​( pip_ipi​ : 内側の流体圧、ΔDp\varDelta{ D_p }ΔDp​ : 微小部位の圧力抵抗、 Δx\varDelta{ x }Δx : xxx 軸上の微小変化量、 Δy\varDelta{ y }Δy : yyy 軸上の微小変化量)

これで圧力に関する抗力は計算できるはずです。

 摩擦力に関する立式空気抵抗には圧力抵抗だけではなくもう一つ重要な要素があります。摩擦抵抗です。

ですが、摩擦抵抗を求めようにも流体は固体と違って塊ではやってこないので少し考え方を工夫する必要があります。

まず流体を塊として捉えるために、適当に流体を切り取ったものが境界表面を滑っていくことを考えます。

そして先の圧力抵抗の式を変形させて、まず垂直抗力を求めます。
ΔN=[kg/m3]⋅[m3]⋅[(m/s)2]⋅[m−1]　よりΔN=ρ⋅12∣ΔxΔy∣⋅(v0p0sin⁡θρc2)2⋅1Δx2+Δy2=12ρv02(p0ρc2)2∣ΔxΔy3∣Δx2+Δy23⋯2◯
\begin{aligned}
\varDelta{ N } &= [ kg / m^3 ] \cdot [m^3] \cdot [ ( m/s )^2 ] \cdot [ m^{-1} ]　より\\
\varDelta{ N } &= \rho{} \cdot{} \frac{ 1 }{ 2 } | \varDelta{ x } \varDelta{ y } | \cdot{} ( \frac{ v_0 p_0 \sin{} \theta{} }{ \rho{} c^2 } )^2 \cdot{} \frac{ 1 }{ \sqrt{ \varDelta{ x^2 } + \varDelta{ y^2 } } } \\
&= \frac{ 1 }{ 2 } \rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2 \frac{ | \varDelta{ x } \varDelta{ y^3 } | } { \sqrt{ \varDelta{ x^2 } + \varDelta{ y^2 } } ^3 } \cdots{} \textcircled{ 2 }
\end{aligned}
ΔNΔN​=[kg/m3]⋅[m3]⋅[(m/s)2]⋅[m−1]　より=ρ⋅21​∣ΔxΔy∣⋅(ρc2v0​p0​sinθ​)2⋅Δx2+Δy2​1​=21​ρv02​(ρc2p0​​)2Δx2+Δy2​3∣ΔxΔy3∣​⋯2◯​( ΔN\varDelta{ N }ΔN : 微小部位の垂直抗力 )

そしてここに「仮想摩擦係数」とでも呼ぶべきものを導入します。

(文字では μ\mu{}μ とおくことにします)

そうすれば、自ずと摩擦抵抗の式が出来上がります。

ただし、摩擦抵抗は摩擦する向きと抗力にカウントされる力の向きは必ずしも一致しない点について注意してください。

そのため、摩擦力に cos⁡θ\cos{} \theta{}cosθ を噛ませて対応します。

よって摩擦力による微小部位の空気抵抗 ΔDf\varDelta{ D_f }ΔDf​ は
ΔDf=12μρv02(p0ρc2)2∣Δx2Δy3∣Δx2+Δy24⋯3◯
\varDelta{ D_f } = \frac{ 1 }{ 2 } \mu{} \rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2 \frac{ | \varDelta{ x^2 } \varDelta{ y^3 } | } { \sqrt{ \varDelta{ x^2 } + \varDelta{ y^2 } } ^4 } \cdots{} \textcircled{ 3 }
ΔDf​=21​μρv02​(ρc2p0​​)2Δx2+Δy2​4∣Δx2Δy3∣​⋯3◯となります。

μ\mu{}μ が具体的にどのような値になるのかについては「 μ\mu{}μ は具体的にいくつ？」の章で計算します。

 全抗力に関する立式この 1◯\textcircled{ 1 }1◯ (微小部位の圧力抵抗)と 3◯\textcircled{ 3 }3◯ (微小部位の摩擦抵抗)の和が、微小部位の全抗力 ΔD\varDelta{ D }ΔD を構成します。

そしてこの ΔD\varDelta{ D }ΔD を物体の全周にわたって積算した結果が、その物体が受ける空気抵抗の総和 DDD となります。

そして DDD から流体密度や速度などの余計なものを飛ばした値が、最終的な抗力係数面積 CDAC_D ACD​A となります。

式に表すと以下のようになります。
ΔD=12ρv02(p0ρc2)2∣Δy5∣+μ∣(Δx2Δy3)∣Δx2+Δy24D=12ρv02(p0ρc2)2∑∣Δy5∣+μ∣(Δx2Δy3)∣Δx2+Δy24CDA=(p0ρc2)2∑∣Δy5∣+μ∣(Δx2Δy3)∣Δx2+Δy24
\begin{aligned}
\varDelta{ D } &= \frac{ 1 }{ 2 } \rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2 \frac{ | \varDelta{ y^5 } | + \mu{} | ( \varDelta{ x^2 } \varDelta{ y^3 } ) | }{ \sqrt{ \varDelta{ x }^2 + \varDelta{ y }^2 }^4 } \\
D &= \frac{ 1 }{ 2 } \rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2 \sum \frac{ | \varDelta{ y^5 } | + \mu{} | ( \varDelta{ x^2 } \varDelta{ y^3 } ) | }{ \sqrt{ \varDelta{ x }^2 + \varDelta{ y }^2 }^4 } \\
C_D A &= ( \frac{ p_0 }{ \rho{} c^2 } )^2 \sum \frac{ | \varDelta{ y^5 } | + \mu{} | ( \varDelta{ x^2 } \varDelta{ y^3 } ) | }{ \sqrt{ \varDelta{ x }^2 + \varDelta{ y }^2 }^4 } \\
\end{aligned}
ΔDDCD​A​=21​ρv02​(ρc2p0​​)2Δx2+Δy2​4∣Δy5∣+μ∣(Δx2Δy3)∣​=21​ρv02​(ρc2p0​​)2∑Δx2+Δy2​4∣Δy5∣+μ∣(Δx2Δy3)∣​=(ρc2p0​​)2∑Δx2+Δy2​4∣Δy5∣+μ∣(Δx2Δy3)∣​​
 手を動かしてみる
 μ\mu{}μ は具体的にいくつ？では μ\mu{}μ の値を計算していきましょう。

ここでは、このサイトから流れに直交する平板の値を取って使います。
まずは圧力抵抗について求めましょう。1◯\textcircled{ 1 }1◯の式を使うと、まず
12ρv02(p0ρc2)2
\frac{ 1 }{ 2 } \rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2
21​ρv02​(ρc2p0​​)2が板の表と裏で2回必要なので
ρv02(p0ρc2)2
\rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2
ρv02​(ρc2p0​​)2が求まります。

続いて摩擦抗力については 2◯\textcircled{ 2 }2◯ の式から垂直抗力を求めるのですが、流れに直交する平板なので Δx\varDelta{ x }Δx が 000 になってしまいます。一方で Δx\varDelta{ x }Δx の代わりに何の値を代入してもOKであるので Δy\varDelta{ y }Δyを代入して
12ρv02(p0ρc2)2123
\frac{ 1 }{ 2 } \rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2 \frac{ 1 }{ \sqrt{ 2 }^3 }
21​ρv02​(ρc2p0​​)22​31​となるので、摩擦抵抗はそのまま μ\mu{}μ を追加して
12ρv02(p0ρc2)2μ23
\frac{ 1 }{ 2 } \rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2 \frac{ \mu{} }{ \sqrt{ 2 }^3 }
21​ρv02​(ρc2p0​​)22​3μ​となります。

さて、ここで流れに直交する平板の抗力係数は 1.121.121.12 であることから、以下の式を μ\mu{}μ について解くことで μ\mu{}μ が求められます。
1.12=(p0ρc2)2(2+μ23)
1.12 = ( \frac{ p_0 }{ \rho{} c^2 } )^2( 2 + \frac{ \mu{} }{ \sqrt{ 2 }^3 } )
1.12=(ρc2p0​​)2(2+2​3μ​)ここで、
p0=101325[Pa]ρ=1.225[kg/m3]c=340.6520138[m/s]
\begin{aligned}
p_0 &= 101325[ Pa ] \\
\rho{} &= 1.225[ kg / m^3 ] \\
c &= 340.6520138[ m/s ]
\end{aligned}
p0​ρc​=101325[Pa]=1.225[kg/m3]=340.6520138[m/s]​を代入することで、
μ=0.578280037531404⋯
\mu{} = 0.578280037531404 \cdots{}
μ=0.578280037531404⋯が求まりました。
さて、この値は果たして正しいのでしょうか。円柱にかかる抗力係数の計算をすることで実証実験をしましょう。

 円柱にかかる抗力本題です。真打です。以下のコードで実験いたします。
ソースコードcylinder5.js
/*
    パラメータ群
*/
/* 円柱の半径 */
const D0    = 1 / 2;                // 円柱の半径[m]
/* 算出精度のパラメータ */
const N0    = 360;                  // 分割数[個]
/* 抗力係数用比例定数算出に用いるパラメータ */
const P0    = 101325;               // 1気圧[Pa]
const C0    = 340.6520138;          // 15℃時の音速[m/s]
const RHO0  = 1.225;                // 空気の密度[kg/m^3]
const CD0   = 0.63;                 // 実験的に求めた抗力係数
/* 摩擦抵抗算出用の仮想摩擦係数 */
const MU0   = 0.578280037531404;    // 仮想摩擦係数
/* グラフ表示のパラメータ */
const W0    = 720;                  // 表示するグラフの横幅[px]
const H0    = 480;                  // 表示するグラフの高さ[px]

/*
    空配列：描画および計算のメッシュ量に応じたもの
*/
const arr = [ ...Array( N0 ) ].map( ( k, i ) => i );

/*
    形状から圧力の抗力係数を計算する関数
*/
// in  : 円(複素数リスト)
// out : 抗力係数(複素数)(実部：圧力の抗力係数, 虚部：摩擦の抗力係数)
const findKCd = l => math.multiply(
    ( ( P0 / ( RHO0 * C0 ** 2 ) ) ** 2 ),                               // 比例定数
    [ ...Array( l.length ) ].map( ( k, i ) => {
        const dX = l[ i ].re - l[ i < l.length - 2 ? i + 1 : 0 ].re;    // X軸上の変化量
        const dY = l[ i ].im - l[ i < l.length - 2 ? i + 1 : 0 ].im;    // Y軸上の変化量
        const dS = math.sqrt( dX ** 2 + dY ** 2 );                      // 変化前と後の距離
        const dP = dS == 0 ? 0 : dY ** 5 / dS ** 4;                     // 圧力抵抗
        const dF = dS == 0 ? 0 : dX ** 2 * dY ** 3 / dS ** 4 * MU0;     // 摩擦抵抗
        return( math.complex( math.abs( dP ), math.abs( dF ) ) );
    } ).reduce( ( p, c ) => math.add( p, c ), 0 ),
    ( 1 / ( D0 * 2 ) )                                                  // 前方投影長さ
);

/*
    変数・リスト群
*/
/* 円(複素数リスト) */
const lstCircle = arr.map( n => math.complex( { r: D0, phi: n / ( arr.length / 2 ) * math.PI } ) );
/* 45°傾いた正方形(複素数リスト) */
const lstSquare = [ math.complex( D0, 0 ), math.complex( 0, D0 ), math.complex( -1 * D0, 0 ), math.complex( 0, -1 * D0 ) ];
/* 円の抗力係数(複素数)(実部：圧力の抗力係数, 虚部：摩擦の抗力係数) */
const cmpKCd    = findKCd( lstCircle );
/* 正方形の抗力係数(複素数)(実部：圧力の抗力係数, 虚部：摩擦の抗力係数) */
const cmpHCd    = findKCd( lstSquare );
/* 抗力係数(合計)(実数) */
const fltKCd    = math.round( ( cmpKCd.re + cmpKCd.im ) * 1000 ) / 1000;
const fltHCd    = math.round( ( cmpHCd.re + cmpHCd.im ) * 1000 ) / 1000;
/* 圧力の抗力係数(実数リスト) */
const lstCdp    = [ CD0, cmpKCd.re, cmpHCd.re ]; 
/* 摩擦の抗力係数(実数リスト) */
const lstCdf    = [ 0,   cmpKCd.im, cmpHCd.im ]; 

/*
    実際の描画
*/
/* cdの名が付いたdivタグの内容として表示 */
Plotly.newPlot( 'cd',
    [       // 表示するグラフ
        {
            name:           'D<sub>P</sub>',
            x:              [ 'Conventional', 'Kabasawa Model', 'Hasegawa Model' ], // X軸:抗力係数を求める方式
            y:              lstCdp,                                                 // Y軸:抗力係数リスト
            type:           'bar',                                                  // 棒グラフ指定
            text:           lstCdp.map( String ),                                   // 名前分け
            marker: {
                color:      [                                                       // 色分け
                    'rgba( 219, 64, 82, 0.7 )',                                     // 実験的に求めた抗力係数
                    'rgba( 55, 128, 191, 0.7 )',                                    // 樺沢式圧力抗力係数
                    'rgba( 142, 124, 195, 0.7 )'                                    // 長谷川式圧力抗力係数
                ],
            },
            textposition:   'auto',
            hovertext:      'Conventional: ' + CD0 + '<br>Kabasawa Model: ' + fltKCd + '<br>Hasegawa Model:' + fltHCd
        }, {
            name:           'D<sub>F</sub>',
            x:              [ 'Conventional', 'Kabasawa Model', 'Hasegawa Model' ], // X軸:抗力係数を求める方式
            y:              lstCdf,                                                 // Y軸:抗力係数リスト
            type:           'bar',                                                  // 棒グラフ指定
            text:           lstCdf.map( String ),                                   // 名前分け
            marker: {
                color:      [                                                       // 色分け
                    'rgba( 219, 64, 82, 0.7 )',                                     // 実験的に求めた抗力係数
                    'rgba( 50, 171, 96, 0.7 )',                                     // 樺沢式摩擦抗力係数
                    'rgba( 55, 83, 109, 0.7 )'                                      // 長谷川式摩擦抗力係数
                ],
            },
            textposition:   'auto',
            hovertext:      'Conventional: ' + CD0 + '<br>Kabasawa Model: ' + fltKCd + '<br>Hasegawa Model:' + fltHCd

        }
    ], {    // レイアウト設定
        height:     H0,         // グラフの高さ
        width:      W0,         // グラフの幅
        barmode:    'stack',    // 積み重ね棒グラフ
        showlegend: false,      // 凡例は隠す
        xaxis:  {               // X軸:抗力係数を求める方式
            title:  'finding method'
        },
        yaxis:  {               // Y軸:抗力係数
            range:  [ 0, 1.2 ],
            title:  'C<sub>D</sub>'
        },
        title:      'Drag Coefficiency on a 2D Cylinder'
    }, {    // グラフ描画のプロパティ
        editable:               true,
        scrollZoom:             true,
        showLink:               false,
        displaylogo:            false,
        modeBarButtonsToRemove: ['sendDataToCloud']
    }
);
cylinder5.html
<!doctype html>
<html lang="ja">
    <head>
        <script src="https://unpkg.com/mathjs/lib/browser/math.js"></script>
        <script src="https://cdn.plot.ly/plotly-latest.min.js"></script>
    </head>
    <body>
        <div id="cd"><!-- ここに抗力係数のグラフを表示する --></div>
        <script src="cylinder5.js"></script>
    </body>
</html>

 ＼結果発表～！！！／
赤い棒が実測値、青い棒が算出した圧力抵抗、その上の緑色の棒が算出した摩擦抵抗です。

0.630.630.63 (実測値)vs 約 0.6250.6250.625 (今回の算出値)。 約 0.0050.0050.005 の差です。この差はデカいのでしょうか。

詳しい方は優しめにふんわりとマサカリを投げてください。

また、紫色(圧力抵抗)と灰色(摩擦抵抗)で長谷川様の形状単純化モデル

(最前端、最後端、最上端、最下端をそれぞれ計測し、それらで作る四角形について計算する)

を取り入れた計算の結果を示しました。値は 0.3010.3010.301 と正直なところあまり正しくない結果が返ってきました。

別の図形(例えば、翼型など)ならばより良い結果が得られるのかの検証は必要ですが、少なくとも単純な円柱においては正しくない結果のようです。

 考察もどき今回は流れに直交する平板から「仮想摩擦係数」を約 0.580.580.58 と算出して、それを円柱にかかる抗力の計算に混ぜ込むというアプローチで攻めてみたわけですが、さすがに値はピタリと一致とはいきませんでしたが割といい線を行っているのではないでしょうか？

(ちなみに手で値をいじって探索した結果、約 0.60.60.6 がピタリと一致するための「仮想摩擦係数」のようです。)

ただ、任意の形状についてはどのように「仮想摩擦係数」を設定すればよいのか、あるいは 0.580.580.58 のままで十分近似値を求めることができるのか、知る方法が現状ありません。

 課題とかまずは「仮想摩擦係数」が 0.580.580.58 のままでよいのかという検証が必要です。

また、この抗力を求める方式を応用して「循環」に依存しない揚力計算も可能なのではないかと考えられるので、その方式を開発することが一つの大きなポイントです。

さらに、上は全て二次元でのお話なので、三次元に拡張した場合どのような変化があるのかについて考える必要があります。

 結論らしきもの以下の考え方を導入して、新たな空気抵抗の計算方法を作成してみました。
長谷川様からインスパイアされた、「音の性質を勘案する」という発想
ベルヌーイの定理から圧力抵抗の値を計算する式の立式
「仮想摩擦係数」の概念を取り入れた摩擦抵抗の値を計算する式の立式
結果としては実用に耐えうるとは断言できない微妙な結果となりましたが、あくまで「概算」程度には使えるのではないかと思います。

 最後に、マサカリ募集のお知らせここまで偉そうにペラペラと理屈を並べましたが、僕は航空力学についてはずぶのド素人です。

故にいい感じのマサカリは超大歓迎です。

航空力学に詳しい諸兄のご意見、ご感想をお待ちしております。

あ、マサカリを投げるときは優しくふんわりめでお願いしますね。痛いので。

 2024/11/27 : 追記初歩的な計算間違いがあったので少し内容を書き換えました。

結果的に修正前と比べて算出した数値の信頼性が上がったという嬉しい誤算はありましたが、それはそれとして間違った内容を掲出してしまい申し訳ございません。

今後もアップデートする内容があれば随時追記いたします。

 2024/12/2 : 追記長谷川様の形状単純化モデルが通用するのかどうかを確かめるべく、グラフに結果を追記しました。

今回の結果としては良くないものでしたが、今後形状の一般化を模索する際においてはそれなりに考慮すべき内容かと思います。

 2024/12/11 : 追記ジュコーフスキー翼の場合どのようになるのかについて、新たにコードを書いたので追記します。

一般的に、翼形状の場合は抗力算出の公式が少し変化し、
D[N]=12ρv2cCD([kg/m3]⋅[(m/s)2]⋅[m2])
D[N] = \frac{ 1 }{ 2 } \rho{} v^2 c C_D ( [kg/m^3] \cdot{} [(m/s)^2] \cdot{} [m^2] )
D[N]=21​ρv2cCD​([kg/m3]⋅[(m/s)2]⋅[m2])( ccc は翼のコード長)

となるようです。何故なのでしょうか。詳しい方は教えてくださると非常にありがたいです。
さて、上記した僕の算出方法では、一つ大きな問題が存在しています。それは、同じ前方投影長さの場合、周囲長が長ければ長いほど算出抗力が青天井状態で増え続けてしまうという点です。

今回のジュコーフスキー翼の場合では、それに対応すべく、算出された CDAC_D ACD​A を P/cP/cP/c ( PPP は周囲長 )で割るという措置を施しました。

正直なところ、先に結果から言ってしまうと正解っぽいのですが、何故正解っぽいのかが上記のコード長を使用する理由と同様全く理解できていません。こちらについても詳しい方はゆるふわなマサカリを投げていただけると幸いです。

 手を動かした部分ソースコードは以下です。
ソースコードzhukovsky.js
/*
    パラメータ群
*/
/* 算出精度のパラメータ */
const N0    = 361;                  // メッシュ分割数[個]
/* 写像前の円のパラメータ */
const A0    = 1;                    // 写像の定数1(任意)
const R0    = 1.13105;              // 写像前の円の半径(翼厚)(15%になる値)
const X0    = A0 * ( -1 * R0 + 1 ); // 写像前の円の中心(x軸)(コード長)
const Y0    = 0;                    // 写像前の円の中心(y軸)(キャンバー(反り))
/* 抗力係数用比例定数算出用のパラメータ */
const P0    = 101325;               // 1気圧[Pa]
const C0    = 340.6520138;          // 摂氏15度時の音速[m/s]
const RHO0  = 1.225;                // 空気の密度[kg/m^3]
/* 摩擦抵抗算出用の仮想摩擦係数 */
const MU0   = 0.578280037531404;    // 仮想摩擦係数
/* グラフ表示のパラメータ */
const W0    = 720;                  // 表示するグラフの横幅[px]
const H0    = 480;                  // 表示するグラフの高さ[px]
/* グラフ描画のプロパティ */
const PROP0 = {
    editable:               true,
    scrollZoom:             true,
    showLink:               false,
    displaylogo:            false,
    modeBarButtonsToRemove: [ 'sendDataToCloud' ]
};

/*
    共通配列(キーを角度に、インデックスをイテレータに使う)
*/
const arr = [ ...Array( N0 ) ].map( ( n, i ) => i - ( N0 - 1 ) / 2 );

/*
    共通関数群
*/
/* 角度(度数)をラジアンに変換する */
// in  : 角度[deg](実数)
// out : 角度[rad](実数)
const convertDeg2Rad = d => d / ( ( N0 - 1 ) / 2 ) * math.PI;
/* 複素数のリストの実数軸と虚数軸をそれぞれ実数リスト化する */
// in  : 複素数リスト
// out : 複素数リストの実部(実数リスト)
const listRe = l => arr.map( ( n, i ) => l[ i ].re );
// in  : 複素数リスト
// out : 複素数リストの虚部(実数リスト)
const listIm = l => arr.map( ( n, i ) => l[ i ].im );
/* 最大値と最小値の差を返す */
// in  : 実数リスト
// out : 最大値と最小値の差(実数)
const findBiggestDiff = l => math.abs( math.max( ...l ) - math.min( ...l ) );

/*
    翼型の描画
*/
/* 関数群 */
// 写像前の円を描く
// in  : 半径(実数), 円の中心のX座標(実数), 円の中心のY座標(実数)
// out : 円(複素数リスト)
const makeCircle = ( r, x, y ) => arr.map( ( n, i ) => {
    const re = r * math.cos( convertDeg2Rad( n ) ) + x;
    const im = r * math.sin( convertDeg2Rad( n ) ) + y;
    return math.complex( re, im );
} );
// 写像後の翼型を描く
// in : 写像前の円(複素数リスト)
// out: 写像後の翼型(複素数リスト)
const makeFoil = l => arr.map( ( n, i ) => {
    const tmp  = math.pow( math.divide( math.subtract( l[ i ], A0 ), math.add( l[ i ], A0 ) ), 2 );
    return math.chain( 2 * A0 ).multiply( math.divide( math.add( 1, tmp ), math.subtract( 1, tmp ) ) ).done();
} );
/* 変数・リスト群 */
// 写像前の円(複素数リスト)
const lstCircle = makeCircle( R0, X0, math.abs( Y0 ) );
// 写像後の翼型(複素数リスト)
const lstFoil = makeFoil( lstCircle );
// 本来の翼型の座標(実数リスト)
const lstFoilX = listRe( lstFoil );  // X座標
const lstFoilY = listIm( lstFoil );  // Y座標
// グラフ用の変数、条件式
const chord  = findBiggestDiff( lstFoilX );
const thick  = findBiggestDiff( lstFoilY );
const ratio  = chord / thick;
/* 実際の描画 */
// 表示するグラフの配列を入れる
const datFoil = [
    {
        // 写像前の円
        x:      listRe( lstCircle ),
        y:      listIm( lstCircle ),
        type:   'scatter',
        model:  'lines',
        name:   'original circle'
    }, {
        // 写像後の翼型
        x:      lstFoilX,
        y:      lstFoilY,
        type:   'scatter',
        model:  'lines',
        name:   'zhukovsky airfoil'
    }
];
// レイアウト設定
const layFoil = {
    height: H0,
    width:  W0,
    xaxis:  {
        // X軸の範囲：本来の翼型の先端から後端まで
        range: [ math.min( ...lstFoilX ), math.max( ...lstFoilX ) ],
        title: 'chord length: ' + chord * 1000 + '[mm]'
    },
    yaxis:  {
        // X軸とY軸の縮尺率は1:1になるように設定
        scaleanchor:    'x',
        scaleratio:     1,
        title:          'Max thickness: ' + thick * 1000 + '[mm]'
    },
    title:  'Zhukovsky Airfoil ( Ratio: ' + ratio + ' )'
};
// foilの名が付いたdivタグの内容として表示
Plotly.newPlot( 'foil', datFoil, layFoil, PROP0 );

/*
    抗力係数と揚力係数のグラフ化
*/
/* 関数群 */
// 翼型の周囲長を求める
// in : 実部(実数リスト), 虚部(実数リスト)
// out: 翼型の周囲長(実数)
const findFoilPerimeter = ( lx, ly ) => arr.map( ( n, i ) => {
    return( math.sqrt( findListDiff( lx, i ) ** 2 + findListDiff( ly, i ) ** 2 ) );
} ).reduce( ( p, c ) => p + c, 0 );
// 複素数平面上の回転を行う
// in : 実部(実数リスト), 虚部(実数リスト), 回転角度(radではなくdeg単位)(実数)
// out: 回転された翼型(複素数リスト)
const rotateFoil = ( lx, ly, d ) => arr.map( ( n, i ) => {
    return( math.complex( {
        r:      math.sqrt( lx[ i ] ** 2 + ly[ i ] ** 2 ),
        phi:    math.atan2( ly[ i ], lx[ i ] ) + convertDeg2Rad( d )
    } ) );
} );
/* リストの現在の要素と次の要素との差を求める */
// in  : 実数リスト, インデックス(実数), 配列の長さ(実数)
// out : リストの現在の要素と次の要素との差(実数)
const findListDiff = ( l, i ) => l[ i ] - l[ i < l.length - 2 ? i + 1 : 0 ];
// 樺沢式抗力係数面積の計算(実数)
// in  : 翼型(複素数リスト)
// out : 翼型にかかる抗力係数面積(実数)
const findKCdA = l => ( P0 / ( RHO0 * C0 ** 2 ) ) ** 2 * arr.map( ( n, i ) => {                    
    const dX  = findListDiff( listRe( l ), i );                     // X軸上の変化量(実数)
    const dY  = findListDiff( listIm( l ), i );                     // Y軸上の変化量(実数)
    const dS  = math.sqrt( dX ** 2 + dY ** 2 );                     // 変化前と後の距離(実数)
    const dDp = dS == 0 ? 0 : dX ** 0 * dY ** 5 / dS ** 4;          // 圧力抗力
    const dDf = dS == 0 ? 0 : dX ** 2 * dY ** 3 / dS ** 4 * MU0;    // 摩擦抗力
    return( math.abs( dDp + dDf ) );
} ).reduce( ( p, c ) => p + c, 0 );
/* 変数・リスト群 */
// 周囲長(実数)
const fltPerimeter = findFoilPerimeter( lstFoilX, lstFoilY );
// 翼型(複素数リストのリスト)
const lstFoilList = arr.map( ( n, i ) => rotateFoil( lstFoilX, lstFoilY, n ) );
// 樺沢式抗力係数(実数リスト)
const lstCd = arr.map( ( n, i ) => findKCdA( lstFoilList[ i ] ) / ( fltPerimeter / chord ) );
/* 実際の描画 */
// 表示するグラフの配列を入れる
const datDrag = [
    {
        // ジュコーフスキー翼の抗力係数(提唱する求め方)
        x:      arr,
        y:      lstCd,
        type:   'scatter',
        model:  'lines',
        name:   'C<sub>D</sub>(original)'
    }
];
// レイアウト設定
const layDrag = {
    height: H0,
    width:  W0,
    xaxis:  {
        // X軸の範囲：arrの最小値から最大値まで
        range:  [ math.min( ...arr ), math.max( ...arr ) ],
        title:  'angles of yaw[deg]'
    },
    yaxis:  {
        // Y軸の範囲：抗力係数の最低値から最高値まで
        range:  [ math.floor( math.min( ...lstCd ) ), math.ceil( math.max( ...lstCd ) ) ],
        title:  'C<sub>D</sub> of Zhukovsky Airfoil'
    },
    title:  'C<sub>D</sub> of Zhukovsky Airfoil'
};
// dragの名が付いたdivタグの内容として表示
Plotly.newPlot( 'drag', datDrag, layDrag, PROP0 );
cylinder5.html
<!doctype html>
<html lang='ja'>
    <head>
        <script src="https://unpkg.com/mathjs/lib/browser/math.js"></script>
        <script src="https://cdn.plot.ly/plotly-latest.min.js"></script>
    </head>
    <body>
        <div id="foil"><!-- ここに翼型のグラフを表示する --></div>
        <div id="drag"><!-- ここに抗力のグラフを表示する --></div>
        <script src="zhukovsky.js"></script>
    </body>
</html>

 結果
上の画像でオレンジ色の線で描かれた図形が今回抗力計算をしたジュコーフスキー翼です。コード長に対して翼厚が15%になるよう調整してあります。
青線のグラフが実際に算出された抗力です。このサイトの下のグラフに載っているNACA0015翼型(厳密には違いますが、コード長に対して翼厚が15%になっているソックリなものです)の抗力係数グラフとほぼ一致していることから、今回の算出は概ね正しそうに見えます。

 考察っぽい部分今回の結果は僕の算出方法を翼型という実用で用いられる図形について一つの正しそうな使い方を提唱できたという点である意味大きいものです。

ただし、何がどう正しいのかの裏付けがないだけでなく、翼型を含む全ての図形についての一般論である「同じ前方投影長さで周囲長が違う場合、どのようにして一般化すればよいのかが不明である」という新たな問題を呼び起こしてしまっていると言えます。

さらなる議論と考察が必要とされるため、しばらく考えてみます。
Discussion

ログインするとコメントできます