斜方投射 空気抵抗が速度に比例する場合
空気抵抗が速度に比例する場合の斜方投射については、"斜方投射"
でインターネットで検索すると、詳しい説明があるのでそちらを参照してください。
計算式についても、インターネットで参照して下さい。
空気抵抗が、速度に比例する条件は投射速度が非常に遅い場合で、発生することが少ないのですが、何故か、インターネットで検索すると、速度に比例する例が沢山出てきます。
水平より下方に対する投射の計算はプログラムにありません。
投射する物体の形状は、球以外には対応していません。
必用であればプログラムを修正してください。
γ:空気抵抗係数計算ボタンをクリックすると、空気粘性係数と、球の半径から、空気抵抗係数を計算します、その時、memo表示欄に、レイノルズ数が１以下になる速度が表示されます。
レイノルズ数が１以上になる初速度を与えると、空気抵抗が、速度に比例しなくなります。
球の場合半径が、０．１mmで速度、０．１７ｍ／ｓ程度です。
球の半径を変更して、質量の変更が必要な場合は、m:質量計算を実行してください。
雨粒の半径は、０．１～８mmで平均１mm程度らしいので、雨粒の落下は、速度の二乗に比例した空気の抵抗が発生する慣性抵抗によるものです。
計算してみると、実際の落下速度に近似します。
空気抵抗が速度に比例する場合は、速度の二乗に比例する場合と違って、粘性による抵抗計算なので、時間無限大の時の到達点があります、粘性により、水平方向への移動が停止する計算となります。
プログラム
unit Main;
interface
uses
Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
Vcl.Controls, Vcl.Forms, Vcl.Dialogs, VclTee.TeeGDIPlus, Vcl.StdCtrls,
Vcl.Buttons, Vcl.ExtCtrls, VCLTee.TeEngine, VCLTee.Series, VCLTee.TeeProcs,
VCLTee.Chart;
type
TForm1 = class(TForm)
Chart1: TChart;
Series1: TLineSeries;
GroupBox1: TGroupBox;
etaEdit: TLabeledEdit;
RadiusEdit: TLabeledEdit;
gammaEdit: TLabeledEdit;
GammaBtn: TBitBtn;
Memo1: TMemo;
GroupBox2: TGroupBox;
mEdit: TLabeledEdit;
VEdit: TLabeledEdit;
angleEdit: TLabeledEdit;
CalcBtn: TBitBtn;
gEdit: TLabeledEdit;
tEdit: TLabeledEdit;
Series2: TLineSeries;
RadioGroup1: TRadioGroup;
GroupBox3: TGroupBox;
RelativeEdit: TLabeledEdit;
mBtn: TBitBtn;
CheckBox1: TCheckBox;
Series3: TLineSeries;
Series4: TLineSeries;
Label1: TLabel;
Label2: TLabel;
Label3: TLabel;
Label4: TLabel;
procedure GammaBtnClick(Sender: TObject);
procedure FormCreate(Sender: TObject);
procedure CalcBtnClick(Sender: TObject);
procedure mBtnClick(Sender: TObject);
private
{ Private 宣言 }
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
uses system.Math;
// γ:空気抵抗係数計算 (gamma)
// 空気抵抗が速度に比例するレイノルズ１の速度の計算
procedure TForm1.GammaBtnClick(Sender: TObject);
var
a : double; // 球半径
eta : double; // 空気粘性係数
gamma : double; // 空気抵抗係数
ch : integer;
rv : double; // レイノルズ数1の速度
begin
memo1.Clear;
val(etaEdit.Text, eta, ch);
if ch <> 0 then begin
memo1.Text := 'η:空気粘性係数に間違いがあります。';
exit;
end;
val(RadiusEdit.Text, a, ch);
if ch <> 0 then begin
memo1.Text := 'a:球の半径に間違いがあります。';
exit;
end;
gamma := 6 * pi * a * eta;
gammaedit.Text := floatTostr(gamma); // 空気抵抗係数
rv := eta / a;
memo1.Lines.Append('レイノルズ数 Re <= 1 の速度' + floatTostr(rv) + ' m');
end;
// 球の質量計算
// 計算は簡単だが手間を省くため
procedure TForm1.mBtnClick(Sender: TObject);
var
a : double; // 球半径
r : double; // 密度
m : double; // 質量
ch : integer;
begin
val(RadiusEdit.Text, a, ch);
if ch <> 0 then begin
memo1.Text := 'a:球の半径に間違いがあります。';
exit;
end;
val(RelativeEdit.Text, r, ch);
if ch <> 0 then begin
memo1.Text := '比重に間違いがあります。';
exit;
end;
m := a * a * a * pi * 4 / 3 * r; // 球の質量計算
medit.Text := floatTostrF(m * 1000, ffGeneral, 10,5);
end;
// 斜方投射計算 下方への投射は無し
// (放物線計算)
procedure TForm1.CalcBtnClick(Sender: TObject);
var
gamma : double; // γ 空気抵抗係数
m : double; // m 質量
v0 : double; // 初速度
deg : double; // 角度 deg
rad : double; // 角度 rad
g : double; // 重力加速度
Lg : double; // Lg = γ/m
v0x : double; // x方向速度
v0y : double; // y方向速度
vxt : double; // t秒後x方向速度
vyt : double; // t秒後y方向速度
xt : double; // t秒後x位置
yt : double; // t秒後y位置
vxtn : double; // 空気抵抗なしt秒後x方向速度
vytn : double; // 空気抵抗なしt秒後y方向速度
xtn : double; // 空気抵抗なしt秒後x位置
ytn : double; // 空気抵抗なしt秒後y位置
th : double; // 最高点到達時間
h : double; // 最高点
hx : double; // 最高点x
h0x : double; // 最高到達点
vyf : double; // 時間∞落下速度
ch : integer;
dt : double; // Δt
j : integer; // 分割数
t : double; // 時間
tm : double; // グラフ最大時間
y : double; // y位置
sc : double; // スケール
i: Integer; //
function calc_vxt(t: double): double; // t秒後x方向速度
begin
result := v0x * exp(-Lg * t);
end;
function calc_vyt(t: double): double; // t秒後y方向速度
begin
result := -g / Lg + 1 / Lg * (g + Lg * v0y) * exp(-Lg * t);
end;
function calc_xt(t: double): double; // t秒後x位置
begin
result := v0x / Lg * (1 - exp(-Lg * t));
end;
function calc_yt(t: double): double; // t秒後y位置
begin
result := -g * t / Lg + 1 / Lg / Lg * (g + Lg * v0y) * (1 - exp(-Lg * t));
end;
function calc_th: double; // 最高点到達時間
begin
result := 1 / Lg * ln((g + Lg * v0y) / g);
end;
function calc_h: double; // 最高点
begin
result := -g / Lg / Lg * ln((g + Lg * v0y) / g) + (v0y / Lg);
end;
function calc_hx: double; // 最高点x
begin
result := v0x * v0y / (g + Lg * v0y);
end;
function calc_h0x: double; // 時間∞到達点
begin
result := v0x / Lg;
end;
function calc_vyf: double; // 時間∞落下速度
begin
result := m * g / gamma;
end;
function xtoy(x: double): double; // X->y 変換
begin
if x = 0 then begin
result := 0;
exit;
end;
if 1 - x * Lg / v0x < 0 then begin // ln演算エラー回避
result := x / v0x * (v0y + 1 / Lg * g) + 1 / (Lg * Lg) * g;
exit;
end;
result := x / v0x * (v0y + 1 / Lg * g) + 1 / (Lg * Lg) * g * ln(1 - x * Lg / v0x);
end;
function no_drag_x(t: double): double; // 空気抵抗無x距離
begin
result := v0x * t;
end;
function no_drag_y(t: double): double; // 空気抵抗無y距離
begin
result := -1 / 2 * g * t * t + v0y * t;
end;
function no_drag_vy(t: double): double; // 空気抵抗無y速度
begin
result := -g * t + v0y;
end;
begin
memo1.Clear;
val(gammaedit.Text, gamma, ch);
if ch <> 0 then begin
memo1.Text := 'γ:空気抵抗係数に間違いがあります。';
exit;
end;
val(medit.Text, m, ch);
if ch <> 0 then begin
memo1.Text := 'm:質量に間違いがあります。';
exit;
end;
val(VEdit.Text, v0, ch);
if ch <> 0 then begin
memo1.Text := 'v:初速度に間違いがあります。';
exit;
end;
val(angleedit.Text, deg, ch);
if ch <> 0 then begin
memo1.Text := 'θ:投射角度に間違いがあります。';
exit;
end;
if deg > 90 then begin
memo1.Text := 'θ:投射角度は90°以内にして下さい。';
exit;
end;
if deg < 0 then begin
memo1.Text := 'θ:投射角度は0(ゼロ)以上にしてください。';
exit;
end;
val(gEdit.Text, g, ch);
if ch <> 0 then begin
memo1.Text := 'g:重力加速度に間違いがあります。';
exit;
end;
val(tEdit.Text, tm, ch);
if ch <> 0 then begin
memo1.Text := 't:グラフ時間に間違いがあります。';
exit;
end;
rad := deg / 180 * pi; // 投射角度 rad単位
Lg := gamma / m; // 空気抵抗係数Γ
v0x := cos(rad) * v0; // x方向初速度
v0y := sin(rad) * v0; // y方向初速度
h := calc_h;
memo1.Lines.Append('非常に速度が遅い場合のみ速度に抵抗が比例します。');
memo1.Lines.Append('計算結果は、理論値です。');
memo1.Lines.Append(' 最高点 ' + floatTostr(h) + ' m');
hx := calc_hx;
memo1.Lines.Append(' 最高点x ' + floatTostr(hx) + ' m');
th := calc_th;
memo1.Lines.Append(' 最高点到達時間 ' + floatTostr(th) + ' sec');
h0x := calc_h0x;
memo1.Lines.Append(' 時間∞到達点x ' + floatTostr(h0x) + ' m');
vyf := calc_vyf;
memo1.Lines.Append(' 時間∞落下速度 ' + floatTostr(vyf) + ' m/s');
j := 1000; // 計算分割数
dt := tm / 1000; // Δt
sc := 1000; // ミリ単位スケール
Series1.Clear;
Series2.Clear;
Series3.Clear;
Series4.Clear;
case RadioGroup1.ItemIndex of
0: begin
Chart1.LeftAxis.Title.Caption := 'mm';
Chart1.BottomAxis.Title.Caption := 'msec';
end;
1: begin
Chart1.LeftAxis.Title.Caption := 'mm';
Chart1.BottomAxis.Title.Caption := 'mm';
end;
2, 3, 4: begin
Chart1.LeftAxis.Title.Caption := 'mm/ms';
Chart1.BottomAxis.Title.Caption := 'msec';
end;
end;
Chart1.LeftAxis.MaximumOffset := 10;
Label1.Visible := false;
Label2.Visible := false;
Label3.Visible := false;
Label4.Visible := false;
case RadioGroup1.ItemIndex of
0: begin
Label1.Visible := true;
Label2.Visible := true;
Label1.Caption := '抵抗有X';
Label2.Caption := '抵抗有Y';
if checkbox1.Checked = true then begin
Label3.Visible := true;
Label4.Visible := true;
Label3.Caption := '抵抗無X';
Label4.Caption := '抵抗無Y';
end;
end;
1: begin
Label1.Visible := true;
Label2.Visible := true;
Label1.Caption := '抵抗有Y';
Label2.Caption := '抵抗有Y';
if checkbox1.Checked = true then begin
Label3.Visible := true;
Label3.Caption := '抵抗無Y';
end;
end;
2: begin
Label1.Visible := true;
Label2.Visible := true;
Label1.Caption := '抵抗有X速度';
Label2.Caption := '抵抗有Y速度';
if checkbox1.Checked = true then begin
Label3.Visible := true;
Label4.Visible := true;
Label3.Caption := '抵抗無X速度';
Label4.Caption := '抵抗無Y速度';
end;
end;
3: begin
Label1.Visible := true;
Label1.Caption := '抵抗有X速度';
if checkbox1.Checked = true then begin
Label3.Visible := true;
Label3.Caption := '抵抗無X速度';
end;
end;
4: begin
Label2.Visible := true;
Label2.Caption := '抵抗有Y速度';
if checkbox1.Checked = true then begin
Label4.Visible := true;
Label4.Caption := '抵抗無Y速度';
end;
end;
end;
vxtn := v0x; // 空気抵抗無 x方向速度は変化なし
for i := 0 to j do begin
t := i * dt;
xt := calc_xt(t); // t秒後x位置
yt := calc_yt(t); // t秒後y位置
vxt := calc_vxt(t); // t秒後x方向速度
vyt := calc_vyt(t); // t秒後y方向速度
y := xtoy(xt); // x位置->y位置
xtn := no_drag_x(t); // 空気抵抗無t秒後x位置
ytn := no_drag_y(t); // 空気抵抗無t秒後y位置
vytn := no_drag_vy(t); // 空気抵抗無t秒後y速度
case RadioGroup1.ItemIndex of
0: begin
Series1.AddXY(t * sc, xt * sc);
Series2.AddXY(t * sc, yt * sc);
if checkbox1.Checked = true then begin
Series3.AddXY(t * sc, xtn * sc);
Series4.AddXY(t * sc, ytn * sc);
end;
end;
1: begin
Series1.AddXY(xt * sc, yt * sc);
Series2.AddXY(xt * sc, y * sc);
if checkbox1.Checked = true then begin
Series3.AddXY(xtn * sc, ytn * sc);
end;
end;
2: begin
Series1.AddXY(t * sc, vxt * sc);
Series2.AddXY(t * sc, vyt * sc);
if checkbox1.Checked = true then begin
Series3.AddXY(t* sc, vxtn * sc);
Series4.AddXY(t* sc, vytn * sc);
end;
end;
3: begin
Series1.AddXY(t * sc, vxt * sc);
if checkbox1.Checked = true then begin
Series3.AddXY(t * sc, vxtn * sc);
end;
end;
4: begin
Series2.AddXY(t * sc, vyt * sc);
if checkbox1.Checked = true then begin
Series3.AddXY(t * sc, vytn * sc);
end;
end;
end;
end;
end;
procedure TForm1.FormCreate(Sender: TObject);
begin
memo1.Clear;
end;
end.
Projectile_motion_One.zip