🚲

【マサカリ求む】高さが1mの円柱にかかる空気抵抗を完全に理解したかった

2022/06/18に公開

前回のおさらい

https://zenn.dev/onlinsanity/articles/2c62f247c6525e
ここを参照していただければ大体分かるとは思うのですが、

  • 直径:高さが1:1な円柱の空気抵抗について調べてみた
  • 音の性質を勘案して立式してみたところそれっぽいグラフが出た
  • 直径:高さを変えたらどうなるか分からない

ってなところです。
立式してみたものが以下です。

v_a[m/s] = \frac{vp \sin^2 \theta}{\rho c^2} ([m/s]・[Pa])/([kg/m^3]・[(m/s)^2])
D[N] = \sum \frac{ 1 }{ 2 } \rho v_a^2 h \sqrt{\varDelta x^2 + \varDelta y^2} [kg/m^3]・[(m/s)^2]・[m^2]

遅すぎた気付き

直径:高さを変えたらどうなるか分からないとは書きましたが、
その大きな理由が 高さ の方を変更していたからです。
故に計算しても結果がまともに現れず、
「あれれ~、おかしいぞ~(某名探偵風に)」 となっていました。
そしてつい先日、とうとう気が付いたのです。
「これって比率さえ変わりゃ良いんだから直径を縮めればいいんじゃね…?」
何故こんな簡単なことに気が付かなかったのか… orz

いざ実践

気付いてしまえばなんてことはありません。
前回と同じ手法で手を動かしてしまえば良いのです。
ソースコードは以下に載せますが、わざと長ったらしく書いたので縮めておきます。

ソースコード
cylinder2.js
// 定数群
const N = 1000;                                                 // 分割数[個]
const d = [ 1, ( 1 / 2 ), ( 1 / 5 ), ( 1 / 10 ), ( 1 / 40 ) ];  // 円柱の直径[m]
const l = 1;                                                    // 円柱の長さ[m]
const cd = [ 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 arr = [ ...Array( N ) ].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 foil0 = cylinder( d[ 0 ] );
const foil1 = cylinder( d[ 1 ] );
const foil2 = cylinder( d[ 2 ] );
const foil3 = cylinder( d[ 3 ] );
const foil4 = cylinder( d[ 4 ] );

// 作った円を点にプロット
const xyList0 = arr.reduce( ( p, i ) => {
    const tmp = foil0( i * 2 * Math.PI / N );
    p.x.push( tmp.x );
    p.y.push( tmp.y );
    return( p );
}, {
    x: [],
    y: []
} );
const xyList1 = arr.reduce( ( p, i ) => {
    const tmp = foil1( i * 2 * Math.PI / N );
    p.x.push( tmp.x );
    p.y.push( tmp.y );
    return( p );
}, {
    x: [],
    y: []
} );
const xyList2 = arr.reduce( ( p, i ) => {
    const tmp = foil2( i * 2 * Math.PI / N );
    p.x.push( tmp.x );
    p.y.push( tmp.y );
    return( p );
}, {
    x: [],
    y: []
} );
const xyList3 = arr.reduce( ( p, i ) => {
    const tmp = foil3( i * 2 * Math.PI / N );
    p.x.push( tmp.x );
    p.y.push( tmp.y );
    return( p );
}, {
    x: [],
    y: []
} );
const xyList4 = arr.reduce( ( p, i ) => {
    const tmp = foil4( i * 2 * Math.PI / N );
    p.x.push( tmp.x );
    p.y.push( tmp.y );
    return( p );
}, {
    x: [],
    y: []
} );

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

// 従来の式での抗力のリストを作る
const yCDragList0 = arr.map( ( i ) => xSpeedList[ i ] ** 2 * density * d[ 0 ] * l * cd[ 0 ] / 2 );
const yCDragList1 = arr.map( ( i ) => xSpeedList[ i ] ** 2 * density * d[ 1 ] * l * cd[ 1 ] / 2 );
const yCDragList2 = arr.map( ( i ) => xSpeedList[ i ] ** 2 * density * d[ 2 ] * l * cd[ 2 ] / 2 );
const yCDragList3 = arr.map( ( i ) => xSpeedList[ i ] ** 2 * density * d[ 3 ] * l * cd[ 3 ] / 2 );
const yCDragList4 = arr.map( ( i ) => xSpeedList[ i ] ** 2 * density * d[ 4 ] * l * cd[ 4 ] / 2 );

// 速度対樺沢式抗力のリストを作る
const kDrag = ( v, list ) => arr.map( ( i ) => {
    const tmpX = list.x[ i ] - list.x[ i < N - 1 ? i + 1 : 0 ];
    const tmpY = list.y[ i ] - list.y[ i < N - 1 ? i + 1 : 0 ];
    const delta = Math.sqrt( tmpY ** 2 + tmpX ** 2 );
    const accel = ( ( tmpY / delta ) ** 2 * pressure * v / ( density * sonic ** 2 ) );
    return( accel ** 2 * density * delta * l / 2 );
} ).reduce( ( p, c ) => p + c, 0 );
const yKDragList0 = arr.map( ( i ) => kDrag( xSpeedList[ i ], xyList0 ) );
const yKDragList1 = arr.map( ( i ) => kDrag( xSpeedList[ i ], xyList1 ) );
const yKDragList2 = arr.map( ( i ) => kDrag( xSpeedList[ i ], xyList2 ) );
const yKDragList3 = arr.map( ( i ) => kDrag( xSpeedList[ i ], xyList3 ) );
const yKDragList4 = arr.map( ( i ) => kDrag( xSpeedList[ i ], xyList4 ) );

// プロットした速度対抗力を線グラフとして繋ぐ
const sCDrag0 = {
    x: xSpeedList,
    y: yCDragList0,
    type: 'scatter',
    mode: 'lines',
    name: 'Conventional(1:1)'
};
const sKDrag0 = {
    x: xSpeedList,
    y: yKDragList0,
    type: 'scatter',
    mode: 'lines',
    name: 'Kabasawa(1:1)'
};
const sCDrag1 = {
    x: xSpeedList,
    y: yCDragList1,
    type: 'scatter',
    mode: 'lines',
    name: 'Conventional(2:1)'
};
const sKDrag1 = {
    x: xSpeedList,
    y: yKDragList1,
    type: 'scatter',
    mode: 'lines',
    name: 'Kabasawa(2:1)'
};
const sCDrag2 = {
    x: xSpeedList,
    y: yCDragList2,
    type: 'scatter',
    mode: 'lines',
    name: 'Conventional(5:1)'
};
const sKDrag2 = {
    x: xSpeedList,
    y: yKDragList2,
    type: 'scatter',
    mode: 'lines',
    name: 'Kabasawa(5:1)'
};
const sCDrag3 = {
    x: xSpeedList,
    y: yCDragList3,
    type: 'scatter',
    mode: 'lines',
    name: 'Conventional(10:1)'
};
const sKDrag3 = {
    x: xSpeedList,
    y: yKDragList3,
    type: 'scatter',
    mode: 'lines',
    name: 'Kabasawa(10:1)'
};
const sCDrag4 = {
    x: xSpeedList,
    y: yCDragList4,
    type: 'scatter',
    mode: 'lines',
    name: 'Conventional(40:1)'
};
const sKDrag4 = {
    x: xSpeedList,
    y: yKDragList4,
    type: 'scatter',
    mode: 'lines',
    name: 'Kabasawa(40:1)'
};

// 表示するグラフの配列を入れる
const veloData0 = [ sCDrag0, sKDrag0 ];

// レイアウト設定
const veloLayout0 = {
    height: graphHeight,
    width: graphWidth,
    xaxis: {
        range: [ Math.min( ...xSpeedList ), Math.max( ...xSpeedList ) ],
        title: 'velocity[m/s]'
    },
    yaxis: {
        range: [ Math.min( ...yCDragList0 ), Math.max( ...yCDragList0 ) ],
        title: 'aerodynamic drag[N]'
    },
    title: 'Aerodynamic Drag on the Cylinder(1:1) vs. Velocity'
};

// 実際の描画(velo0の名が付いたdivタグの内容として表示)
Plotly.newPlot( 'velo0', veloData0, veloLayout0, {
    editable: true,
    scrollZoom: true,
    showLink: false,
    displaylogo: false,
    modeBarButtonsToRemoie: ['sendDataToCloud']
} );

// 表示するグラフの配列を入れる
const veloData1 = [ sCDrag1, sKDrag1 ];

// レイアウト設定
const veloLayout1 = {
    height: graphHeight,
    width: graphWidth,
    xaxis: {
        range: [ Math.min( ...xSpeedList ), Math.max( ...xSpeedList ) ],
        title: 'velocity[m/s]'
    },
    yaxis: {
        range: [ Math.min( ...yCDragList1 ), Math.max( ...yCDragList1 ) ],
        title: 'aerodynamic drag[N]'
    },
    title: 'Aerodynamic Drag on the Cylinder(2:1) vs. Velocity'
};

// 実際の描画(velo0の名が付いたdivタグの内容として表示)
Plotly.newPlot( 'velo1', veloData1, veloLayout1, {
    editable: true,
    scrollZoom: true,
    showLink: false,
    displaylogo: false,
    modeBarButtonsToRemoie: ['sendDataToCloud']
} );

// 表示するグラフの配列を入れる
const veloData2 = [ sCDrag2, sKDrag2 ];

// レイアウト設定
const veloLayout2 = {
    height: graphHeight,
    width: graphWidth,
    xaxis: {
        range: [ Math.min( ...xSpeedList ), Math.max( ...xSpeedList ) ],
        title: 'velocity[m/s]'
    },
    yaxis: {
        range: [ Math.min( ...yCDragList2 ), Math.max( ...yCDragList2 ) ],
        title: 'aerodynamic drag[N]'
    },
    title: 'Aerodynamic Drag on the Cylinder(5:1) vs. Velocity'
};

// 実際の描画(velo0の名が付いたdivタグの内容として表示)
Plotly.newPlot( 'velo2', veloData2, veloLayout2, {
    editable: true,
    scrollZoom: true,
    showLink: false,
    displaylogo: false,
    modeBarButtonsToRemoie: ['sendDataToCloud']
} );

// 表示するグラフの配列を入れる
const veloData3 = [ sCDrag3, sKDrag3 ];

// レイアウト設定
const veloLayout3 = {
    height: graphHeight,
    width: graphWidth,
    xaxis: {
        range: [ Math.min( ...xSpeedList ), Math.max( ...xSpeedList ) ],
        title: 'velocity[m/s]'
    },
    yaxis: {
        range: [ Math.min( ...yCDragList3 ), Math.max( ...yCDragList3 ) ],
        title: 'aerodynamic drag[N]'
    },
    title: 'Aerodynamic Drag on the Cylinder(10:1) vs. Velocity'
};

// 実際の描画(velo0の名が付いたdivタグの内容として表示)
Plotly.newPlot( 'velo3', veloData3, veloLayout3, {
    editable: true,
    scrollZoom: true,
    showLink: false,
    displaylogo: false,
    modeBarButtonsToRemoie: ['sendDataToCloud']
} );

// 表示するグラフの配列を入れる
const veloData4 = [ sCDrag4, sKDrag4 ];

// レイアウト設定
const veloLayout4 = {
    height: graphHeight,
    width: graphWidth,
    xaxis: {
        range: [ Math.min( ...xSpeedList ), Math.max( ...xSpeedList ) ],
        title: 'velocity[m/s]'
    },
    yaxis: {
        range: [ Math.min( ...yCDragList4 ), Math.max( ...yCDragList4 ) ],
        title: 'aerodynamic drag[N]'
    },
    title: 'Aerodynamic Drag on the Cylinder(40:1) vs. Velocity'
};

// 実際の描画(velo0の名が付いたdivタグの内容として表示)
Plotly.newPlot( 'velo4', veloData4, veloLayout4, {
    editable: true,
    scrollZoom: true,
    showLink: false,
    displaylogo: false,
    modeBarButtonsToRemoie: ['sendDataToCloud']
} );

// 表示するグラフの配列を入れる
const totalData = [ sCDrag0, sKDrag0, sCDrag1, sKDrag1, sCDrag2, sKDrag2, sCDrag3, sKDrag3, sCDrag4, sKDrag4 ];

// レイアウト設定
const totalLayout = {
    height: graphHeight,
    width: graphWidth,
    xaxis: {
        range: [ Math.min( ...xSpeedList ), Math.max( ...xSpeedList ) ],
        title: 'velocity[m/s]'
    },
    yaxis: {
        range: [ Math.min( ...yCDragList0 ), Math.max( ...yCDragList0 ) ],
        title: 'aerodynamic drag[N]'
    },
    title: 'Aerodynamic Drag on the Cylinders vs. Velocity'
};

// 実際の描画(velo0の名が付いたdivタグの内容として表示)
Plotly.newPlot( 'total', totalData, totalLayout, {
    editable: true,
    scrollZoom: true,
    showLink: false,
    displaylogo: false,
    modeBarButtonsToRemoie: ['sendDataToCloud']
} );

\結果発表~!!/

all.png
全体的な傾向としては割といい感じに見えますよね。
音速時における最大誤差も3kN程度ですし。
ですが、個々を見ていくと次第に「あれ?」と思ってしまう事態になってゆきます。
1_1.png
これは前回の記事と一緒のグラフです。
音速時における最大誤差は3kN弱です。
それでは次のグラフへ。
2_1.png
直径を50cmにしてみました。
音速時における最大誤差は上より若干増えましたがまだ3kN弱です。
あれ、空気抵抗の絶対値は減ったのに誤差があまり変わらない…?
5_1.png
直径が20cmになりました。
音速時における最大誤差は2kN強です。
またしても空気抵抗の絶対値が減った割に誤差が変化しない…。
それに、グラフの開きが酷くなってきたような…。
10_1.png
直径が10cmになりました。
音速時における最大誤差は2kN弱。
この段階においても空気抵抗の絶対値に対して誤差が大きく変わらない…?
そして、どんどん開きが酷くなってまいりました。
40_1.png
直径は2.5cmです。
音速時における最大誤差は約700N。
逆に何故この条件だと誤差がそれっぽく縮まったのでしょうか?
しかしこの開きはえげつないですね。一瞬「半分ぐらい?」と思いかねません。

考察らしきもの

これは勝手な僕の想像なのですが、風洞実験によって確認した抗力係数は、
風洞の床面に棒を立てて計測したのではないか と思います。
故に地面側の端に当たる空気は床にさえぎられるために棒の端より外に逃げられず、
棒に当たることを強制された結果空気抵抗が上がるのでは、と思うのです。
一方で、僕の計算式で出される値に関しては、
風洞の中に棒が浮遊している理想の状態で計測した結果 に近い数値なのかな、と思っています。
もっとも、風洞実験装置なんて持っていないので確認のしようがありませんが…。
さらに言えば、どうやって浮遊させろと。
よしんば浮遊させることが出来たとして、どうやって抗力を求めろと。
計測装置自体が空気抵抗を生むようでは話が変わってきてしまいますからね。

まとめ

以前に提唱した式を直径:高さの比率を変えた上で適用してみたのですが、

  • 何故かほとんどの条件で実測値と計算値の最大誤差が2~3kN前後
  • 直径:高さの比率が上がれば上がるほどグラフの開きが酷い

原因は追って調査します。

【緩募】マサカリ

前回も書きましたが、僕は航空力学についてはずぶのド素人です。
故にいい感じのマサカリは超大歓迎です。
航空力学に詳しい諸兄のご意見、ご感想をお待ちしております。

Discussion