🚲

【マサカリ求む】任意の円柱にかかる空気抵抗を完全に理解した。ただし理屈は不明。

2022/11/23に公開約7,600字

前回のおさらい

https://zenn.dev/onlinsanity/articles/19b7e0d17f8122
ここを参照していただければ大体わかるのですが、

  • 直径:高さが1:1でない場合の空気抵抗値を計算出来るかやってみた
  • 直径:高さのアスペクト比が大きくなるにつれて誤差も増えた
  • 何でかはよく分からなかった

ってなところです。

式に潜んでいた致命的間違い

アスペクト比が開くと誤差が開く理由について色々と頭を悩ませていたのですが、
ある時、閃いて計算した結果から抗力係数を逆算してみました
するとどうでしょう。
アスペクト比が変わっても抗力係数が全く変化していないじゃん。
何という大ポカ。
そもそも、数式を注意深く見れば発見出来た問題だったのです。

\begin{aligned} v_a&=\frac{vp\sin{}^2\theta{}}{\rho{}c^2} \\ &=\frac{vp}{\rho{}c^2}\cdot{}\frac{\varDelta{y^2}}{\varDelta{x^2}+\varDelta{y^2}} \end{aligned}

わざと現実的な計算に使いやすい\varDelta{x}\varDelta{y}を使った表現に書き換えましたが、
これって図形が等比変形したらv_aが結果的に同じになっちゃいますよね
道理でアスペクト比が変わっても抗力係数が一緒な訳だ…orz
一応確認はしたのですが、円柱の直径でなく高さを変えてアスペクト比を変化させた場合においても似たような結果が得られました。残念。

既知の結果から逆算してみよう

ならばと思い、既知の抗力係数とアスペクト比をプロットして、近似した関数を探したところ、厳密には違いましたが大体以下のような関数が得られました。

0.1\cdot{}\ln{(\frac{h}{d})}+0.63

ここでhは円柱の高さ、dは円柱の直径です。
0.63はアスペクト比が1:1の場合の抗力係数です。
ならば、これを式に組み込めば行けるんと違いますか?
とりあえず手を動かしてみましょう。

いざ実践

ソースコードは以下に載せます。
前回よりは省略できる所を省略したので賢そうなソースになったと思います。

ソースコード
cylinder3.js
// 定数群
const N = 1000;                                                         // 分割数[個]
const ds = [ 1 / 40, 1 / 10, 1 / 5, 1 / 2, 1, 1, 1, 1, 1 ];             // 円柱の直径[m]
const hs = [ 1, 1, 1, 1, 1, 2, 5, 10, 40 ];                             // 円柱の長さ[m]
const cds = [ 0.98, 0.82, 0.74, 0.68, 0.63, 0.68, 0.74, 0.82, 0.98 ];   // 抗力係数(出典:http://skomo.o.oo7.jp/f28/hp28_63.htm)
const sonic = 340.6520138;                                              // 15℃時の音速[m/s]
const density = 1.225;                                                  // 空気の密度[kg/m^3]
const pressure = 101325;                                                // 1気圧[Pa]
const graphWidth = 720;                                                 // 表示するグラフの横幅[px]
const graphHeight = 480;                                                // 表示するグラフの高さ[px]

// 共通で使う配列
const nar = [ ...Array( N ) ].map( ( k, i ) => i );
const lar = hs.map( ( k, i ) => i );

// 円を描く関数
const cylinder = d => {
    const r = d / 2;
    return ( t ) => {
        return {
            x: r * Math.cos( t ),
            y: r * Math.sin( t )
        };
    };
};

// 円を作る
const foils = lar.map( n => cylinder( ds[ n ] ) );

// 作った円を点にプロット
const xyLists = lar.map( n => nar.reduce( ( p, i ) => {
    const tmp = foils[n]( i * 2 * Math.PI / N );
    p.x.push( tmp.x );
    p.y.push( tmp.y );
    return( p );
}, {
    x: [],
    y: []
} ) );

// 速度のx軸を作る(0~音速まで)
const xSpeedList = nar.map( i => i * sonic / N );

// 従来の式での抗力のリストを作る
const yCDragLists = lar.map( n => nar.map( i => xSpeedList[ i ] ** 2 * density * ds[ n ] * hs[ n ] * cds[ n ] / 2 ) );

// 形状から抗力係数を計算する
const kCd = ( n ) => nar.map( i => {
    const tmpX = xyLists[ n ].x[ i ] - xyLists[ n ].x[ i < N - 1 ? i + 1 : 0 ];
    const tmpY = xyLists[ n ].y[ i ] - xyLists[ n ].y[ i < N - 1 ? i + 1 : 0 ];
    const delta = Math.sqrt( tmpY ** 2 + tmpX ** 2 );
    return ( pressure / ( density * sonic ** 2 ) ) ** 2 * ( ( tmpY ** 4 ) / ( delta ** 3 * ds[ n ] ) );
} ).reduce( ( p, c ) => p + c, 0.1 * Math.log( hs[ n ] / ds[ n ] ) );
// 速度対樺沢式抗力のリストを作る
const kDrag = ( n, i ) => xSpeedList[ i ] ** 2 * density * ds[ n ] * hs[ n ] * kCd( n ) / 2;
const yKDragLists = lar.map( n => nar.map( i => kDrag( n, i ) ) );

// プロットした速度対抗力を線グラフとして繋ぐ
const sCDrags = lar.map( n => {
    return( {
        x: xSpeedList,
        y: yCDragLists[ n ],
        type: 'scatter',
        mode: 'lines',
        name: 'Conventional(' + hs[ n ] + ':' + ds[ n ] + ')'
    } );
} );
const sKDrags = lar.map( n => {
    return( {
        x: xSpeedList,
        y: yKDragLists[ n ],
        type: 'scatter',
        mode: 'lines',
        name: 'Kabasawa(' + hs[ n ] + ':' + ds[ n ] + ')'
    } );
} );

// 表示するグラフの配列を入れる
const veloDatas = lar.map( n => [ sCDrags[ n ], sKDrags[ n ] ] );

// レイアウト設定
const veloLayouts = lar.map( n => {
    return( {
        height: graphHeight,
        width: graphWidth,
        xaxis: {
            range: [ Math.min( ...xSpeedList ), Math.max( ...xSpeedList ) ],
            title: 'velocity[m/s]'
        },
        yaxis: {
            range: [ Math.min( ...yCDragLists[ n ] ), Math.max( ...yCDragLists[ n ] ) ],
            title: 'aerodynamic drag[N]'
            },
        title: 'Aerodynamic Drag on the Cylinder(' + hs[ n ] + ':' + ds[ n ] + ') vs. Velocity'
    } );
} );

// 実際の描画(velo[n]の名が付いたdivタグの内容として表示)
lar.map( n => Plotly.newPlot( 'velo' + n, veloDatas[ n ], veloLayouts[ n ], {
    editable: true,
    scrollZoom: true,
    showLink: false,
    displaylogo: false,
    modeBarButtonsToRemoie: ['sendDataToCloud']
} ) );
cylinder3.html
<!doctype html>
<html lang="ja">
    <head>
        <script src="https://cdn.plot.ly/plotly-latest.min.js"></script>
    </head>
    <body>
        <div id="velo0"><!-- ここに速度対抗力のグラフを表示する --></div>
        <div id="velo1"><!-- ここに速度対抗力のグラフを表示する --></div>
        <div id="velo2"><!-- ここに速度対抗力のグラフを表示する --></div>
        <div id="velo3"><!-- ここに速度対抗力のグラフを表示する --></div>
        <div id="velo4"><!-- ここに速度対抗力のグラフを表示する --></div>
        <div id="velo5"><!-- ここに速度対抗力のグラフを表示する --></div>
        <div id="velo6"><!-- ここに速度対抗力のグラフを表示する --></div>
        <div id="velo7"><!-- ここに速度対抗力のグラフを表示する --></div>
        <div id="velo8"><!-- ここに速度対抗力のグラフを表示する --></div>
        <script src="cylinder3.js"></script>
    </body>
</html>

\結果発表~!!/

一つずつ結果を見ていきましょう。
まずは高さ:直径が40:1のもの。上が直径、下が高さを変えてあります。
1vs0.025.png
40vs1.png
お?
次は高さ:直径が10:1もの。これも上が直径、下が高さを変えてあります。
1vs0.1.png
10vs1.png
これは?
次は高さ:直径が5:1のもの。これも(ry
1vs0.2.png
5vs1.png
まさか?
次は高さ:直径が2:1のもの。こr(ry
1vs0.5.png
2vs1.png
来たか?
次は高さ:直径が1:1のもの。何度も見せたやつです。
1vs1.png
これが一番誤差が大きい気がする…。
結論:無限大が出てこなきゃ大体OK!

考察らしきもの

既知の抗力係数から逆算したとは言え、それっぽく値が出たのは喜ばしいことです。
また、ソースコードの数式をよく読んでいただければ分かると思うのですが、
抗力係数(C_d)の算出が、アスペクト比による部分と、流体の特性による部分と、断面形状による部分とに分かれています。

C_d=0.1\cdot{}\ln{(\frac{h}{d})}+(\frac{p}{\rho{}c^2})^2\cdot{}\sum{}\frac{\varDelta{y^4}}{d\cdot{}\sqrt{\varDelta{x^2}+\varDelta{y^2}}^3}

ここでhは円柱の高さ、dは円柱の直径、pは流体の圧力、\rho{}は流体の密度、cは流体内での音速、\varDelta{x}は図形のx軸上の微小変化量、\varDelta{y}は図形のy軸上の微小変化量を指します。
アスペクト比による部分が何故自然対数関数なのか、何故0.1という定数が引っ付いているのかは全く分かりません。
どなたか詳しい方はマサカリをぶん投げてください。お願いいたします。
また、抗力係数が計算出来たということは、実は揚力係数も計算出来ることになります。
何故なら、

  • 抗力と揚力は直交するので、抗力(D)が求まれば揚力(L)はL=D\cdot{}\frac{\cos{\theta{}}}{\sin{\theta{}}}で求まるはず(\theta{}は仰角)
  • 抗力と揚力の算出式は、\frac{1}{2}v^2\rho{}Sまでは共通(最後に抗力係数C_dか揚力係数C_lかをかける)

よって、揚力係数(C_l)は導出した抗力係数に\frac{\cos{\theta{}}}{\sin{\theta{}}}をかけるだけでよく、

\frac{\cos{\theta{}}}{\sin{\theta{}}}=\frac{1}{\tan{\theta{}}}=\frac{\varDelta{x}}{\varDelta{y}}より \\ C_l=0.1\cdot{}\ln{(\frac{h}{d})}\cdot{}\sum{\frac{\varDelta{x}}{\varDelta{y}}}+(\frac{p}{\rho{}c^2})^2\cdot{}\sum{}\frac{\varDelta{x}\varDelta{y^3}}{d\cdot{}\sqrt{\varDelta{x^2}+\varDelta{y^2}}^3}

となります。まだ翼型断面でこの計算式が通用するか分からないので「知らんけど」が付きますが。

課題

円柱が無限大な 夢の後の~♪ 高さを持っている場合、抗力係数は1.2になるはずなのですが、今回提唱したやり方だと成り立ちません。
また、前述したアスペクト比による部分の説明が出来ない点については調査の続行が必要とされます。
あとは、円柱以外にもこの方法論が通用するかどうかについて、検証が必要です。四角柱とかなら計算出来そうですが、翼型断面の場合は?

まとめ

  • 既知の抗力係数から逆算したところ、アスペクト比の変化を計算上の抗力係数に反映させる方法を発見した。
  • 今回提唱した抗力係数の計算式は、アスペクト比による部分と、流体の特性による部分と、断面形状による部分とに分かれている。
  • 多分、式をちょっと操作すれば揚力も求めることが出来るのではないかと考えられる。
  • 誰か僕に風洞実験施設を買ってください。

【緩募】マサカリ

恒例となっている〆ですが、僕は航空力学についてはずぶのド素人です。
故にいい感じのマサカリは超大歓迎です。
航空力学に詳しい諸兄のご意見、ご感想をお待ちしております。

Discussion

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