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

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

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

そこで色々とGoogle先生に聞いてみるも、ほぼ答えは
$$ D[N] = \frac{ 1 }{ 2 } \rho v^2 S C_d [kg/m^3]･[(m/s)^2]･[m^2] $$
(ρは空気の密度、vは速度、Sは前方投影面積、Cdは抗力係数)
この記述以外がなかなか見当たらねぇ…。と悩んでいたそんなところ、
こんなサイトがあるじゃありませんか。

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

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

とのこと。

じゃあ、詰んだの？

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

固有音響インピーダンス

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

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

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

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

とエビデンスなしに述べております。これには疑問符が付きました。
恐らく、長谷川様の仮説のバックボーンとなっている本、
"Understanding Flight"にそう書いてあったのを鵜呑みにしたのでは、と思います。
その点について長谷川様に問い合わせてみても、

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

という回答で、いまいちピンと来ません。
音響インピーダンスは
$$ \frac{ p }{ Sv } ( \frac{ [Pa] }{ [m^2]･[m/s] } ) $$
で求まるのですが、長谷川式揚力の計算式
$$ L[N] = \frac{ v }{ c } ･ \sin \theta ･ \frac{ p }{ 1atm } ･ 1atm ･ S ･ \cos \theta ( \frac{ [m/s] }{ [m/s] }･\frac{ [Pa] }{ [Pa] }･[Pa]･[m^2] ) $$
(v:流速、c:その場の音速、p:雰囲気圧、S:表面積、θ:表面の傾き)
のどこにも音響インピーダンスを考慮した部分が出てきていません。

2022/01/24追記:

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

とのことです。

さて、何故「音響インピーダンス」にたどり着いたかという話をしていませんでしたが、
気圧差を音速で解消しようとするけれども、「音響インピーダンス」に邪魔されて
実際には音速でない速度で解消されていると考えたからです。
(例えば、団扇で扇いでから実際に風が来るまで意外とタイムラグがあると思います)

話を戻して、数式に空気の固有音響インピーダンスを登場させましょう。
空気の固有音響インピーダンス(厳密には疎密波が正弦波の場合)は、
$$ Z = \rho c ([kg/m^3]･[m/s]) = \frac{ p_a }{ v_a } (\frac{ [Pa] }{ [m/s] }) $$
(ρ:その場の密度、c:その場の音速、pa:その場の気圧、va:その場の速度)
で求まるので、その場の速度vaを求めたければ
$$ v_a = \frac{ p_a }{ \rho c } $$
となります。ところで、長谷川様の仮説では物体の前方で上昇した気圧(＝物体の後方で降下した気圧)は
$$ \frac{ v }{ c } \sin \theta ･ p \sin \theta = \frac{ vp \sin^2 \theta }{ c } $$
となるため、式をまとめると
$$ v_a = \frac{vp \sin^2 \theta}{\rho c^2} $$
となります。僕はこのvaが物体によって曲げられた空気の流速だと判断しました。
(長谷川様の仰る「原理2」において、空気の流れが曲げられる際の作用に相当する流速)
つまり、この流速から力を導けば「原理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] $$
(h:物体の高さ、Δx:x軸上の微小変化量、Δy:y軸上の微小変化量)
と立式してみました。
円柱の場合は、「原理2」の影響が最大になるとのことなので、
おそらくほぼ全ての抗力がこの式で導けると思います。

レッツ実践

グラフの描画にはplotly.jsを使用しています。

cylinder.js

// 定数群
const N = 1000;                                 // 分割数[個]
const d0 = 1;                                   // 円柱の直径[m]
const l0 = 1;                                   // 円柱の長さ[m]
const cd = 0.63;                                // 直径:高さが1:1の円柱の抗力係数
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( d0 );

// 作った円を点にプロット
const xyList = 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: []
} );

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

// 従来の式での抗力のリストを作る
const yCDragList = arr.map( ( i ) => xSpeedList[ i ] ** 2 * density * d0 * l0 * cd / 2 );

// 速度対長谷川式抗力のリストを作る
const hDrag = ( v ) => arr.map( ( i ) => {
	const tmpX = xyList.x[ i ] - xyList.x[ i < N - 1 ? i + 1 : 0 ];
	const tmpY = xyList.y[ i ] - xyList.y[ i < N - 1 ? i + 1 : 0 ];
	const delta = Math.sqrt( tmpY ** 2 + tmpX ** 2 );
	return( ( tmpY / delta ) ** 2 * delta * l0 * pressure * v / sonic );
} ).reduce( ( p, c ) => p + c, 0 );
const yHDragList = arr.map( ( i ) => hDrag( xSpeedList[ i ] ) );

// 速度対長谷川式抗力(前方投影面積版)のリストを作る
const yHDragList2 = arr.map( ( i ) => 0.35 * d0 * l0 * pressure * xSpeedList[ i ] / sonic );

// 速度対樺沢式抗力のリストを作る
const kDrag = ( v ) => arr.map( ( i ) => {
	const tmpX = xyList.x[ i ] - xyList.x[ i < N - 1 ? i + 1 : 0 ];
	const tmpY = xyList.y[ i ] - xyList.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 * l0 / 2 );
} ).reduce( ( p, c ) => p + c, 0 );
const yKDragList = arr.map( ( i ) => kDrag( xSpeedList[ i ] ) );

// プロットした速度対抗力を線グラフとして繋ぐ
const sCDrag = {
	x: xSpeedList,
	y: yCDragList,
	type: 'scatter',
	mode: 'lines',
	name: 'Conventional'
};
const sHDrag = {
	x: xSpeedList,
	y: yHDragList,
	type: 'scatter',
	mode: 'lines',
	name: 'Hasegawa-1'
};
const sHDrag2 = {
	x: xSpeedList,
	y: yHDragList2,
	type: 'scatter',
	mode: 'lines',
	name: 'Hasegawa-2'
};
const sKDrag = {
	x: xSpeedList,
	y: yKDragList,
	type: 'scatter',
	mode: 'lines',
	name: 'Kabasawa-1'
};

// 表示するグラフの配列を入れる
const veloData = [ sCDrag, sHDrag, sHDrag2, sKDrag ];

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

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

cylinder.html

<!doctype html>
<html lang="ja">
	<head>
		<script src="https://cdn.plot.ly/plotly-latest.min.js"></script>
	</head>
	<body>
		<div id="velo"><!-- ここに速度対抗力のグラフを表示する --></div>
		<script src="cylinder.js"></script>
	</body>
</html>

＼結果発表～!!／

青線が従来の式(抗力係数0.63)(Conventional)、
橙線が「原理1」による長谷川式抗力(Hasegawa-1)、
緑線が前方投影面積から求めた「原理1」による長谷川式抗力に0.35をかけたもの(Hasegawa-2)、
赤線が「原理2」のみを考慮した僕の式(Kabasawa-1)です。
青線と赤線がいい感じに近似していますね。
青線と赤線では音速付近で3kN程の差がありますが、それは「原理1」を考慮していないからでしょう。
しかしこれはあくまで「概算」なので、そう考えれば十分許容可能な誤差でしょう。
また、「原理1」のみで計算した緑線が当たらずと言えども遠からずなのも面白いところです。
(長谷川様曰く、緑線の抗力は実測値に基づいて経験的に0.35という係数をかけたそうです)

課題

高さ対直径が変わると抗力係数が変化しますが、それについていけないところです。
また、円柱以外の立体図形に対して有効なのか分かりません。
さらに言えば、長谷川様の仰る「原理1」で求まる揚力/抗力との兼ね合いがどうなるのか、
例えば何か適当な係数を噛ませる等の措置(例えば、真円度など)をすればよいのか、
そこの方法論がサッパリ思いつきません。

追記(2022/01/24)

長谷川様からの返信をいただいた後、自分なりに再考慮して文面を若干修正しました。

【緩募】マサカリ

僕は流体力学や航空力学についてはずぶのド素人なので、いい感じのマサカリが飛んでくると喜びます。
皆様と一緒にこの方式をより良いものとするためにも、是非ご協力をお願いいたします。