free fall 2

自由落下の続き

　自由落下の追加の検討です。
自由落下では、空気抵抗が、速度の二乗に比例する場合、速度に比例する場合、空気の抵抗がない場合等の計算を行っていますが、此処では、速度の二乗に比例して、空気抵抗が発生する場合で、特に、鉛直方向で、下方に初速度をもって投射された場合について、検討します。
此処でのプログラムは、放物線の下向き投射(０～－９０°)の計算に利用します、放物線のけいさんで、下向き投射に対する計算をインターネットで見つける事が出来ませんでした。

鉛直　下向き投射の時は、計算を二つに分けて計算します。

　まず、投射速度が、落下時間が無限大の時の落下速度より早いときです。
この時は、減速が働き、落下時間が無限大の時の速度に近づきます。

　　次は、初速度(V0)が、落下時間が無限大の時の速度V(t∞)より遅い時です。
この時は、落下と共に加速され、無限大の時の速度V(t∞)に近づきます。
この計算は、停止状態から自然落下した時の計算を利用すれば、落下の計算ができます。

　

　問題は、落下(投射)の初速度(V0)が、がV(t∞)より速い時です。
解析解を出すのが難しかったので、分割計算をすることにしました。

dm:　落下時間∞時抗力　　　　　　　dm = V(t∞) * V(t∞) * k
dv: 　抗力
k::　空気抵抗係数 [kg/m]
m:　質量　　　　　　[kg]
V:　初速度　　　　[m/s]
mdg:　加速度　　　[m/s²]
fV: 減速後速度　[m/s]
Δt:　分割時間 [Δs]　　　　　　　　Δt = t　/ 分割数　
L:　距離　　　　　　[m]

１．　dv = V * V * k
２．　mdg = (dv - dm) / m
３．　fV = V - mdg × Δt
４．　L = L + (V + fV) / 2 * Δt;
５．　V = fV

Δtの値を十分に小さくとり、１～５の計算を繰り返せば、目的の値を求める事が出来ます。
実際に計算をしてみると、分割数1000位では、結構誤差が大きく、10万分割位必要となり演算に時間がかかります。
少し改良して、分割数が小さくて済むようにします。
加速度の平均値を計算し、加速度の平均値でΔt後の速度とします。

ndv: 予抗力
nmdg: 予加速度

１．　dv = V * V * k
２．　mdg = (dv - dm) / m
３．　fV = V - mdg × Δt
４．　ndv = fV * fV * k
５．　nmdg = (ndv - dm) / m
６．　mdg = (nmdg + mdg) / 2
７．　fV = V - mdg * Δt
８．　L = L + (V + fV) / 2 * Δt
９．　V = Vf

若干、値が値が大きくなる傾向がありますが、1000分割程度で、4桁の精度を得ることが可能です。

　プログラムには、水平方向で、重力加速度ゼロの計算も、分割計算の検算の為に追加しました。

プログラム

unit main;

interface

uses
  Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
  Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.StdCtrls, Vcl.Buttons, Vcl.ExtCtrls,
  VclTee.TeeGDIPlus, VCLTee.TeEngine, VCLTee.Series, VCLTee.TeeProcs,
  VCLTee.Chart;

type
  TForm1 = class(TForm)
    GroupBox1: TGroupBox;
    rhoEdit: TLabeledEdit;
    sumEdit: TLabeledEdit;
    KcalcBtn: TBitBtn;
    kEdit: TLabeledEdit;
    CdEdit: TLabeledEdit;
    GroupBox2: TGroupBox;
    mEdit: TLabeledEdit;
    tEdit: TLabeledEdit;
    gEdit: TLabeledEdit;
    vhBtn: TBitBtn;
    V0Edit: TLabeledEdit;
    Memo1: TMemo;
    Chart1: TChart;
    Series1: TLineSeries;
    Series2: TLineSeries;
    Chart2: TChart;
    Series3: TLineSeries;
    Series4: TLineSeries;
    jEdit: TLabeledEdit;
    CheckBox1: TCheckBox;
    procedure KcalcBtnClick(Sender: TObject);
    procedure vhBtnClick(Sender: TObject);
  private
    { Private 宣言 }
    procedure down(v0, t, k, m, g: extended; j: integer);
    procedure downN(v0, t, k, m, g: extended; j: integer);
    procedure Horizontal(v0, t, k, m: extended; j: integer);
    procedure HorizontalN(v0, t, k, m: extended; j: integer);
  public
    { Public 宣言 }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}
// 空気抗力係数から空気抵抗計算kg/m計算
procedure TForm1.KcalcBtnClick(Sender: TObject);
var
  k   : extended;   // 空気抵抗係数  kg/m
  Cd  : extended;   // 空気抗力係数
  r   : extended;   // 空気密度   kg/m^3
  S   : extended;   // 代表面積   m^2
  ch  : integer;
begin
  val(CdEdit.Text, Cd, ch);
  if ch <> 0 then begin
    Memo1.Text := 'Cd:空気抗力係数に間違いがあります。';
    exit;
  end;
  val(rhoEdit.Text, r, ch);
  if ch <> 0 then begin
    Memo1.Text := 'ρ:空気密度に間違いがあります。';
    exit;
  end;
  val(sumEdit.Text, S, ch);
  if ch <> 0 then begin
    Memo1.Text := 'S:代表面積に間違いがあります。';
    exit;
  end;
  k := Cd * r * s / 2;
  kEdit.Text := floatTostr(k);
end;

// 水平投射 補正有
procedure TForm1.HorizontalN(v0, t, k, m: extended; j: integer);
var
  dv, dt, mdg, fv : extended;
  ndv, nmdg : extended;
  ti : extended;
  l,  v : extended;
  i: Integer;
begin
  series3.Clear;
  series4.Clear;
  l := 0;
  fv := 0;
  v := abs(v0);
  dt := t / j;
  series3.AddXY(0 , v);
  series4.AddXY(0 , l);
  for i := 1 to j do begin
    ti := i * dt;                  // 時間
    dv := v * v * k;               // 抗力
    mdg := dv / m;                 // 加速度
    fv := v - mdg * dt;            // 予速度
    ndv := fv * fv * k;            // 予抗力
    nmdg := ndv / m;               // 予加速度
    mdg := (mdg + nmdg) / 2;       // 平均加速度
    fv := v - mdg * dt;            // 新速度
    if fv < 0 then break;
    l := l + (v + fv) / 2 * dt;    // 距離
    v := fv;
    series3.AddXY(ti , v);
    series4.AddXY(ti , l);
  end;
  if fv < 0 then begin
    memo1.Lines.Append('新減速計算 分割数不足、時間が長い為計算出来ません。');
    exit;
  end;
  memo1.Lines.Append('新速度 ' + floatTostr(v) + ' m/s');
  memo1.Lines.Append('新距離 ' + floatTostr(l) + ' m');
end;

// 水平投射 補正無
procedure TForm1.Horizontal(v0, t, k, m: extended; j: integer);
var
  dv, dt, mdg, fv : extended;
  ti : extended;
  l,  v : extended;
  i: Integer;
begin
  series1.Clear;
  series2.Clear;
  l := 0;
  fv := 0;
  v := abs(v0);
  dt := t / j;
  series1.AddXY(0 , v);
  series2.AddXY(0 , l);
  for i := 1 to j do begin
    ti := i * dt;                  // 時間
    dv := v * v * k;               // 抗力
    mdg := dv / m;                 // 加速度
    fv := v - mdg * dt;            // 速度
    if fv < 0 then break;
    l := l + (v + fv) / 2 * dt;    // 距離
    v := fv;                       // 新速度
    series1.AddXY(ti , v);
    series2.AddXY(ti , l);
  end;
  if fv < 0 then begin
    memo1.Lines.Append('新減速計算 分割数不足、時間が長い為計算出来ません。');
    exit;
  end;
  memo1.Lines.Append('旧速度 ' + floatTostr(v) + ' m/s');
  memo1.Lines.Append('旧距離 ' + floatTostr(l) + ' m');
end;

// 鉛直降下  補正無
procedure TForm1.down(v0, t, k, m, g: extended; j: integer);
var
  vmax, ti : extended;
  h, v: extended;
  dv, dm, dt, mdg, fv : extended;
  i :  integer;
begin
  series1.Clear;
  series2.Clear;
  h := 0;
  dt := t / j;                      // Δt
  v := -v0;
  fv := 0;
  vmax := sqrt(m * g / k);          // t∞速度
  dm := vmax * vmax * k;            // t∞抗力
  series1.AddXY(0 , -v);
  series2.AddXY(0 , -h);
  for i := 1 to j do begin
    ti := i * dt;                   // 時間
    dv := v * v * k;                // 抗力
    mdg := (dv - dm) / m;           // 加速度
    fv := v - mdg * dt;             // 減速
    if fv < 0 then break;
    h := h + (v + fv) / 2 * dt;     // 高さ + 平均速度*時間
    v := fv;                        // 新速度
    series1.AddXY(ti , -v);
    series2.AddXY(ti , -h);
  end;
  if fv < 0 then begin
    memo1.Lines.Append('旧減速計算 分割数不足、時間が長い為計算出来ません。');
    exit;
  end;
  memo1.Lines.Append('旧速度 ' + floatTostr(-v) + ' m/s');
  memo1.Lines.Append('旧高さ ' + floatTostr(-h) + ' m');
end;

// 鉛直降下 補正有
procedure TForm1.downN(v0, t, k, m, g: extended; j: integer);
var
  vmax, ti : extended;
  h, v: extended;
  dv, dm, dt, mdg, fv : extended;
  ndv, nmdg : extended;
  im, itwo : extended;
  i : integer;
begin
  series3.Clear;
  series4.Clear;
  im := 1 / m;
  itwo := 1 / 2;
  h := 0;
  dt := t / j;                      // Δt
  v := -v0;
  fv := 0;
  vmax := sqrt(m * g / k);          // t∞速度
  dm := vmax * vmax * k;            // t∞抗力
  series3.AddXY(0 , -v);
  series4.AddXY(0 , -h);
  for i := 1 to j do begin
    ti := i * dt;
    dv := v * v * k;                // 抗力
    mdg := (dv - dm) * im;          // 加速度
    fv := v - mdg * dt;             // 予減速
    ndv := fv * fv * k;             // 予抗力
    nmdg := (ndv - dm) * im;        // 予加速度
    mdg := (nmdg + mdg) * itwo;     // 平均加速度
    fv := v - mdg * dt;             // 減速
    if fv < 0 then break;
    h := h + (fv + v) * itwo * dt;  // 高さ + 平均速度*時間
    v := fv;                        // 新速度
    series3.AddXY(ti , -v);
    series4.AddXY(ti , -h);
  end;
  if fv < 0 then begin
    memo1.Lines.Append('新減速計算 分割数不足、時間が長い為計算出来ません。');
    exit;
  end;
  memo1.Lines.Append('新速度 ' + floatTostr(-v) + ' m/s');
  memo1.Lines.Append('新高さ ' + floatTostr(-h) + ' m');
end;

// vh計算実行
procedure TForm1.vhBtnClick(Sender: TObject);
var
  k   : extended;   // 空気抵抗係数  kg/m
  m   : extended;   // 質量          kg
  t   : extended;   // 落下時間      sec
  g   : extended;   // 重力加速度    m/sec^2
  v0  : extended;   // 初速度        m/s
  ch  : integer;
  j   : integer;
begin
  Memo1.Clear;
  val(kEdit.Text, k, ch);
  if ch <> 0 then begin
    Memo1.Text := 'k:空気抵抗係数に間違いがあります。';
    exit;
  end;
  if k <= 0 then begin
    Memo1.Text := 'k:空気抵抗係数の値は0(ゼロ)より大きくしてください。';
    exit;
  end;
  val(mEdit.Text, m, ch);
  if ch <> 0 then begin
    Memo1.Text := 'm:質量に間違いがあります。';
    exit;
  end;
  if m <= 0 then begin
    Memo1.Text := 'm:質量の値は0(ゼロ)より大きくしてください。';
    exit;
  end;
  val(tEdit.Text, t, ch);
  if ch <> 0 then begin
    Memo1.Text := 't:時間に間違いがあります。';
    exit;
  end;
  if t <= 0 then begin
    Memo1.Text := 't:時間の値は0(ゼロ)より大きくしてください。';
    exit;
  end;
  val(gEdit.Text, g, ch);
  if ch <> 0 then begin
    Memo1.Text := 'g:加速度に間違いがあります。';
    exit;
  end;
  if g <= 0 then begin
    Memo1.Text := 'g:加速度の値は0(ゼロ)より大きくしてください。';
    exit;
  end;
  val(v0Edit.Text, v0, ch);
  if ch <> 0 then begin
    Memo1.Text := 'V0:初速度に間違いがあります。';
    exit;
  end;
  val(jEdit.Text, j, ch);
  if ch <> 0 then begin
    Memo1.Text := 'j:分割数に間違いがあります。';
    exit;
  end;
  if j <= 2 then begin
    Memo1.Text := 'j:分割数は2より大きくしてください。';
    exit;
  end;
  if j > 10000000 then begin
    Memo1.Text := 'j:分割数が大きすぎます。';
    exit;
  end;
  memo1.Clear;
  if checkbox1.Checked = false then begin
    down(v0, t, k, m, g, j);
    downN(v0, t, k, m, g, j);
  end
  else begin
    Horizontal(v0, t, k, m, j);
    HorizontalN(v0, t, k, m, j);
  end;
end;

end.

　　　free_fall_test.zip

各種プログラム計算例に戻る