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

三行で

既知の二次元図形にかかる空気抵抗を計算したい！！！

だから自分なりに方法を考えてみました！！！

結果はそこそこいい線行ってるんじゃないかとは思いますが実践で使えるかは正直微妙なところです！！！

注釈

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

思想の話

何故にそんなことを考えたの？

当初の目的は、「空気抵抗を最小限にするにはどんな断面がいいの？」という疑問に対する答えでした。
そのためにこの記事を参考にジュコーフスキー翼を作ってみて、その抗力を調べようと思ったのですが、
「抗力係数が分からねぇ…」という根本的な問題にぶち当たりました(当然ですが)。

で、それじゃあどうやって抗力を求めるの？

そこで色々とGoogle先生に聞いてみるも、ほぼ答えは

$$D[N] = \frac{ 1 }{ 2 } \rho{} v^2 S C_D ( [kg/m^3] \cdot{} [(m/s)^2] \cdot{} [m^2] )$$

( $\rho{}$ は空気の密度、 $v$ は速度、 $S$ は前方投影面積、 $C_D$ は抗力係数)
この記述以外がなかなか見当たらねぇ…。と悩んでいたそんなところ、こんなサイトがあるじゃありませんか。

https://h-a-labo.amebaownd.com/posts/2558036

(現在は撤回されていますが、僕がアクセスした当時は抗力の試算方法も載っていました)
こりゃすげぇ!!と早速試しに円柱の空気抵抗を計算してみたら、既存の式に比べて値が異様に高い…。
(従来の式で約48[N]に対して、上のサイトの方式(以降、「長谷川式」)の抗力は約5023[N])
(円柱は半径0.5[m]の高さ1[m]、標高0[m]で15[℃]の状態、風速は40[km/h]を想定しました)
上のサイトの管理者(以降、長谷川様)に問い合わせても、

この２、３日熟慮しましたが、抗力計算結果の間違いの理由が解明できておりません。

とのこと。

じゃあ、詰んだの？

んな訳はあるめぇ、他に方法があるはずだ、と思い色々な値を見てみました。
その紆余曲折は長い上に面白くないので割愛しますが、ある値に行きつきました。

固有音響インピーダンス

簡単に言えば音の伝わりづらさです。交流電圧とインピーダンスの関係と同じです。
長谷川様の仮説では、

物体が空気中を移動するとき、物体の前面には空気が圧縮されて高気圧域が、後面には空気が膨張して低気圧域が発生する。

ということになっています。これについては僕はとても正しいと思います。ですが同時に、

空気中に生じた気圧差は音速で解消されるのですが、物体の移動が続くと物体の表面に気圧差は発生し続けるのです。

とエビデンスなしに述べております。これには疑問符が付きました。
その点について長谷川様に問い合わせてみても、

小生は「気圧差の解消速度はその場の音の伝播速度」と考えています。もちろん流速でも移動速度でもありません。
「その場の音の伝播速度」とはその場の「音響インピーダンス」を考慮してのことで、上空3000mの音速を考慮してセスナの揚力を計算しています。

という回答で、いまいちピンと来ません。
音響インピーダンスは

$$\frac{ p }{ Sv } ( \frac{ [Pa] }{ [m^2] \cdot{} [m/s] } )$$

で求まるのですが、長谷川式揚力の計算式

$$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] )$$

( $v$ :流速、 $c$ :その場の音速、 $p$ :雰囲気圧、 $S$ :表面積、 $\theta{}$ :表面の傾き)
のどこにも音響インピーダンスを考慮した部分が出てきていません。

「その場の音速」と「雰囲気圧」、で計算することでその場の音響インピーダンスを考慮したつもりでした。

とのことですが…。
さて、何故「音響インピーダンス」にたどり着いたかという話をしていませんでしたが、
気圧差を音速で解消しようとするけれども、「音響インピーダンス」に邪魔されて
実際には音速でない速度で解消されていると考えたからです。
(例えば、団扇で扇いでから実際に風が来るまで意外とタイムラグがあると思います)
話を戻して、数式に空気の固有音響インピーダンスを登場させましょう。
空気の固有音響インピーダンス(厳密には疎密波が正弦波の場合)は、

$$Z = \rho{} c ([kg/m^3] \cdot{} [m/s]) = \frac{ p_a }{ v_a } ( \frac{ [Pa] }{ [ m/s ] } )$$

( $\rho{}$ :その場の密度、 $c$ :その場の音速、 $p_a$ :その場の気圧、 $v_a$ :その場の速度)
で求まるので、その場の速度 $v_a$ を求めたければ

$$v_a = \frac{ p_a }{ \rho{} c }$$

となります。ところで、長谷川様の仮説では物体の前方で上昇した気圧(＝物体の後方で降下した気圧)は

$$\frac{ v_0 }{ c } \sin{} \theta{} \cdot{} p_0 \sin{} \theta{} = \frac{ v_0 p_0 \sin{}^2 \theta{} }{ c }$$

となるため、式をまとめると

$$v_a = \frac{ v_0 p_0 \sin{}^2 \theta{} }{ \rho{} c^2 }$$

となります。僕はこの $v_a$ が物体によって圧縮された空気の基本的な流速だと判断しました。

理論構築の話

圧力に関する立式

上で求めた流速からベルヌーイの定理を用いて圧力を導けば圧力に関する力が導けるはずです。
二次元における非圧縮のベルヌーイの定理では

$$\frac{ p }{ \rho{} } + \frac{ v^2 }{ 2 } = const.$$

なので、仮に流れの外側と内側が同じ流体で構成されていて、内側には一切の流速がなく、かつその境界が無限に薄い何かで覆われていると考えれば下の式が成り立ちます。

$$\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}$$

( $p_i$ : 内側の流体圧、$\varDelta{ D_p }$ : 微小部位の圧力抵抗、 $\varDelta{ x }$ : $x$ 軸上の微小変化量、 $\varDelta{ y }$ : $y$ 軸上の微小変化量)
これで圧力に関する抗力は計算できるはずです。

摩擦力に関する立式

空気抵抗には圧力抵抗だけではなくもう一つ重要な要素があります。摩擦抵抗です。
ですが、摩擦抵抗を求めようにも流体は固体と違って塊ではやってこないので少し考え方を工夫する必要があります。
まず流体を塊として捉えるために、適当に流体を切り取ったものが境界表面を滑っていくことを考えます。
そして先の圧力抵抗の式を変形させて、まず垂直抗力を求めます。

$$\begin{aligned} \varDelta{ N } &= [ \frac{ kg }{ m^3 } ] \cdot [m^3] \cdot [ ( \frac{ m }{ s } )^2 ] \cdot [ m^{-2} ] より\\ \varDelta{ N } &= \frac{ 1 }{ 2 } \rho{} | \varDelta{ x } \varDelta{ y } | ( \frac{ v_0 p_0 \sin{} \theta{} }{ \rho{} c^2 } )^2 \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}$$

( $\varDelta{ N }$ : 微小部位の垂直抗力 )
そしてここに「仮想摩擦係数」とでも呼ぶべきものを導入します。
(文字では $\mu{}$ とおくことにします)
そうすれば、自ずと摩擦抵抗の式が出来上がります。
ただし、摩擦抵抗は摩擦する向きと抗力にカウントされる力の向きは必ずしも一致しない点について注意してください。
そのため、摩擦力に $\cos{} \theta{}$ を噛ませて対応します。
よって摩擦力による微小部位の空気抵抗 $\varDelta{ D_f }$ は

$$\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 }$$

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

全抗力に関する立式

この $\textcircled{ 1 }$ (微小部位の圧力抵抗)と $\textcircled{ 3 }$ (微小部位の摩擦抵抗)の和が、微小部位の全抗力 $\varDelta{ D }$ を構成します。
そしてこの $\varDelta{ D }$ を物体の全周にわたって積算した結果が、その物体が受ける空気抵抗の総和 $D$ となります。
そして $D$ から流体密度や速度などの余計なものを飛ばした値が、最終的な抗力係数面積 $C_D A$ となります。
式に表すと以下のようになります。

$$\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}$$

手を動かしてみる

$\mu{}$ は具体的にいくつ？

では $\mu{}$ の値を計算していきましょう。
ここでは、このサイトから流れに直交する平板の値を取って使います。

まずは圧力抵抗について求めましょう。$\textcircled{ 1 }$の式を使うと、まず

$$\frac{ 1 }{ 2 } \rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2$$

が板の表と裏で2回必要なので

$$\rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2$$

が求まります。続いて摩擦抗力については $\textcircled{ 2 }$ の式から、 $\varDelta{ x }$ の値が0になってしまう一方で、 $\varDelta{ x }$ の代わりに何の値を代入してもOKであるので $\varDelta{ y }$を代入して

$$\frac{ 1 }{ 2 } \rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2 \cdot{} \frac{ 1 }{ \sqrt{2}^3 }$$

の垂直抗力が発生するため、かかる摩擦抗力は

$$\frac{ 1 }{ 2 } \rho{} v_0^2 ( \frac{ p_0 }{ \rho{} c^2 } )^2 \frac{ \mu{} }{ \sqrt{ 2 }^3 }$$

となります。さて、ここで流れに直交する平板の抗力係数は $1.12$ であることから、以下の式を $\mu{}$ について解くことで $\mu{}$ が求められます。

$$1.12 = ( \frac{ p_0 }{ \rho{} c^2 } )^2( 2 + \frac{ \mu{} }{ \sqrt{ 2 }^3 } )$$

ここで、

$$\begin{aligned} p_0 &= 101325[ Pa ] \\ \rho{} &= 1.225[ kg / m^3 ] \\ c &= 340.6520138[ m/s ] \end{aligned}$$

を代入することで、

$$\mu{} = 0.578280037531404 \cdots{}$$

が求まりました。

さて、この値は果たして正しいのでしょうか。円柱にかかる抗力係数の計算をすることで実証実験をしましょう。

円柱にかかる抗力

本題です。真打です。以下のコードで実験いたします。

ソースコード

cylinder5.js

/*
    パラメータ群
*/
/* 算出精度のパラメータ */
const N0    = 360;                  // 分割数[個]
/* 抗力係数用比例定数算出に用いるパラメータ */
const P0    = 101325;               // 1気圧[Pa]
const C0    = 340.6520138;          // 15℃時の音速[m/s]
const MU0   = 0.578280037531404;    // 仮想摩擦係数
const RHO0  = 1.225;                // 空気の密度[kg/m^3]
const CD0   = 0.63;                 // 実験的に求めた抗力係数
/* グラフ表示のパラメータ */
const W0    = 720;                  // 表示するグラフの横幅[px]
const H0    = 480;                  // 表示するグラフの高さ[px]
/* パラメータのリスト */

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

/*
    関数群
*/
/* 円を描く関数 */
// in  : void
// out : 円(複素数リスト)
const makeCircle = () => arr.map( n => math.complex( { r: 1 / 2, phi: n / ( N0 / 2 ) * math.PI } ) );
/* 形状から圧力の抗力係数を計算する */
// in  : 円(複素数リスト)
// out : 複素数(実部：圧力の抗力係数, 虚部：摩擦の抗力係数)
const findKCd = l => math.multiply(
    ( P0 / ( RHO0 * C0 ** 2 ) ) ** 2,   // 比例定数みたいなもの
    arr.map( i => {
        const dX = l[ i ].re - l[ i < N0 - 2 ? i + 1 : 0 ].re;                  // X軸上の変化量
        const dY = l[ i ].im - l[ i < N0 - 2 ? i + 1 : 0 ].im;                  // Y軸上の変化量
        const dS = math.sqrt( dX ** 2 + dY ** 2 );                              // 変化前と後の距離
        const Fp = math.abs( dS == 0 ? 0 : dY ** 5 / dS ** 4 );                 // 圧力抵抗
        const Ff = math.abs( dS == 0 ? 0 : dX ** 2 * dY ** 3 / dS ** 4 * MU0 ); // 摩擦抵抗
        return( math.complex( Fp, Ff ) );
    } ).reduce( ( p, c ) => math.add( p, c ), 0 )
);

/*
    変数群
*/
const cmpKCd = findKCd( makeCircle() );
const kCd    = math.round( ( cmpKCd.re + cmpKCd.im ) * 100 ) / 100;

/*
    リスト群
*/
const lstCdp = [ CD0, cmpKCd.re ];  // 圧力の抗力係数(実数リスト)
const lstCdf = [ 0,   cmpKCd.im ];  // 摩擦の抗力係数(実数リスト)

/*
    実際の描画
*/
/* cdの名が付いたdivタグの内容として表示 */
Plotly.newPlot( 'cd',
    [       // 表示するグラフ
        {
            name:           'D<sub>P</sub>',
            x:              [ 'Conventional', 'Kabasawa 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)'                          // 樺沢式圧力抗力係数
                ],
            },
            textposition:   'auto',
            hovertext:      'Conventional: ' + CD0 + '<br>Kabasawa Model: ' + kCd
        }, {
            name:           'D<sub>F</sub>',
            x:              [ 'Conventional', 'Kabasawa 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)'                         // 樺沢式摩擦抗力係数
                ],
            },
            textposition:   'auto',
            hovertext:      'Conventional: ' + CD0 + '<br>Kabasawa Model: ' + kCd
        }
    ], {    // レイアウト設定
        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 1:1 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.63$ (実測値)vs $0.62$ (今回の算出値)。 $0.01$ の差です。この差はデカいのでしょうか。
詳しい方は優しめにふんわりとマサカリを投げてください。

考察もどき

今回は流れに直交する平板から「仮想摩擦係数」を約 $0.58$ と算出して、それを円柱にかかる抗力の計算に混ぜ込むというアプローチで攻めてみたわけですが、さすがに値はピタリと一致とはいきませんでしたが割といい線を行っているのではないでしょうか？
(ちなみに手で値をいじって探索した結果、約 $0.6$ がピタリと一致するための「仮想摩擦係数」のようです。)
ただ、任意の形状についてはどのように「仮想摩擦係数」を設定すればよいのか、あるいは $0.58$ のままで十分近似値を求めることができるのか、知る方法が現状ありません。

課題とか

まずは「仮想摩擦係数」が $0.58$ のままでよいのかという検証が必要です。
また、この抗力を求める方式を応用して「循環」に依存しない揚力計算も可能なのではないかと考えられるので、その方式を開発することが一つの大きなポイントです。
さらに、上は全て二次元でのお話なので、三次元に拡張した場合どのような変化があるのかについて考える必要があります。

結論らしきもの

以下の考え方を導入して、新たな空気抵抗の計算方法を作成してみました。

• 長谷川隆様からインスパイアされた、「音の性質を勘案する」という発想
• ベルヌーイの定理から圧力抵抗の値を計算する式の立式
• 「仮想摩擦係数」の概念を取り入れた摩擦抵抗の値を計算する式の立式

結果としては実用に耐えうるとは断言できない微妙な結果となりましたが、あくまで「概算」程度には使えるのではないかと思います。

最後に、マサカリ募集のお知らせ

ここまで偉そうにペラペラと理屈を並べましたが、僕は航空力学についてはずぶのド素人です。
故にいい感じのマサカリは超大歓迎です。
航空力学に詳しい諸兄のご意見、ご感想をお待ちしております。
あ、マサカリを投げるときは優しくふんわりめでお願いしますね。痛いので。

2024/11/27 : 追記

初歩的な計算間違いがあったので少し内容を書き換えました。
結果的に修正前と比べて算出した数値の信頼性が上がったという嬉しい誤算はありましたが、それはそれとして間違った内容を掲出してしまい申し訳ございません。
今後もアップデートする内容があれば随時追記いたします。