パラシュートによる空気抵抗が係数ｋ[kg/m]に対してｋ*P
になる場合と、ｋ*Ｐ２になる場合として計算します。自由落下２の続き
自由落下２の追加の検討です。
ここでは、落下の途中で、空気抵抗係数が大きくなる場合についてです。
落下の途中で、空気抵抗係数が著しく大きくなるのは、スカイダイビングによる落下です。
スカイダイビングで、降下しているときの速度は、５０m/s～６０m/sの速度になるようですが、地表に近づいて、パラシュートを開いたとき、どれ位で減速するのかまず考えてみます。
パラシュートを開く前と、開いた後では空気抵抗係数が、１５０倍程度変化します。
パラシュートの開き方についてはスカイダイブ藤岡のパラシュートオープン解説を読んでください。
急激にオープンすると、パラシュートと、スカイダイバーにダメージを与えるので、パラシュートは徐々に開くように工夫されています。
現代のスカイダイビングで利用されているパラシュートは、小さなパラグライダーであり、降下時、３０km/h～５０km/hで前進するので揚力が働き、降下率は３m/ｓ～４m/sで、前進しない時は、５m/s～６m/s位になるようです。
自衛隊でのパラシュートは、円形であり、滑空の降下はあまりなく、飛行機が低空でとび、飛行機から飛び降りてから、少し間をおいて開いているようですが、スカイダイビングの様に長い時間パラシュートを開かずに落下することは無いようです。
スカイダイビングのパラシュートを開くタイミングは、地上から１０００m位なので、パラシュートが無かったら地上まで２０秒弱の時間となります。
メインのパラシュートが開かなかったら、リザーブのパラシュートを開くので５００mで開いて減速をする必要があります。
落下速度６０m/sから５m/s迄、５００mの高さで等減速度で減速するのが、一番負荷が小さくて済みます。
しかし、空気抵抗が、速度の二乗に比例する為に、等加速度で減速する様に、パ゛シュートを開く細工をすることは、非常に困難かと思われます。
そこで、パラシュートが開くときの空気抵抗係数が、ｓｉｎカーブ変化するものとして、計算してみる事にしました。
実際には、スカイダイビングでの、パラシュートが開くときの速度変化を計測すればよいのでしょうが、データーが無いので、代用計算値を使用します。
(パラシュートを製作しているメーカーには、実測したデーターがあるのですが、公開されていないようなので)
此処での計算は、あくまでも、パラシュートがｓｉｎカーブで開くとしたもので、実際の値とは違いますので注意してください。
パラシュートによる空気抵抗が係数ｋ[kg/m]に対してｋ*P
になる場合と、ｋ*Ｐ２になる場合として計算します。
パラシュートにより減速が働き、パラシュートでの落下時間が無限大の時の速度に近づきます。
k : ダイバー空気抵抗係数
kb: パラシュート空気抵抗係数
m: 質量
t: 時間
g: 重力加速度
j := 1000;
dt := t / j;
Qs2 := -pi / 2;
dq := pi / j;
for i := 0 to j do begin
ti := i * dt;
sq := (1 + sin(Qs2 + dq * i)) / 2; //
Q パラシュートの開き係数sin(-pi/2～pi/2)
if
checkbox1.Checked = true then
mkb := k +
kB * sq * sq // 空気抵抗係数計算 sin(Q)^2
else
mkb := k + kB * sq; //
空気抵抗係数計算 sin(Q)
vmax := sqrt(m * g / mkb); //
t∞速度 m/s
dm := vmax * vmax * mkb; //
t∞抗力 m.kg/s^2
dv := v * v * mkb; //
抗力 m.kg/s^2
mdg := (dv - dm) / m; //
加速度 m/s^2
fv := v - mdg * dt; //
予減速 m/s
ndv := fv * fv * mkb; //
予抗力 m.kg/s^2
nmdg := (ndv - dm) / m; //
予加速度 m/s^2
mdg := (nmdg + mdg) / 2; //
平均加速度 m/s^2
fv := v - mdg * dt; //
減速 m/s
if fv < 0 then break;
h := h + (fv + v) /
2 * dt; // 高さ + 平均速度*時間 m
v := fv; // 新速度 m/s
end;
若干、値が値が大きくなる傾向がありますが、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;
gEdit: TLabeledEdit;
vhBtn: TBitBtn;
Memo1: TMemo;
Chart1: TChart;
Series1: TLineSeries;
jEdit: TLabeledEdit;
Series2: TLineSeries;
GroupBox3: TGroupBox;
sum2Edit: TLabeledEdit;
k2Edit: TLabeledEdit;
cd2Edit: TLabeledEdit;
kcalc2Btn: TBitBtn;
tBEdit: TLabeledEdit;
CheckBox1: TCheckBox;
Series3: TPointSeries;
RadioGroup1: TRadioGroup;
Label1: TLabel;
Series4: TPointSeries;
procedure KcalcBtnClick(Sender: TObject);
procedure vhBtnClick(Sender: TObject);
procedure kcalc2BtnClick(Sender: TObject);
procedure FormCreate(Sender: TObject);
private
{ Private 宣言 }
procedure downN(tB, k, kb, m, g: extended; j: integer);
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
// パラシュート空気抗力係数から空気抵抗計算kg/m計算
// パラシュートが開いた時の代表面積で、開く途中の
// 面積は、降下計算の中でsinカーブで計算します。
procedure TForm1.kcalc2BtnClick(Sender: TObject);
var
kp : extended; // 空気抵抗係数 kg/m
Cdp: extended; // 空気抗力係数
r : extended; // 空気密度 kg/m^3
Sp : extended; // 代表面積 m^2
ch : integer;
begin
val(Cd2Edit.Text, Cdp, ch);
if ch <> 0 then begin
Memo1.Text := 'Cdp:空気抗力係数に間違いがあります。';
exit;
end;
val(rhoEdit.Text, r, ch);
if ch <> 0 then begin
Memo1.Text := 'ρ:空気密度に間違いがあります。';
exit;
end;
val(sum2Edit.Text, Sp, ch);
if ch <> 0 then begin
Memo1.Text := 'Sp:代表面積に間違いがあります。';
exit;
end;
kp := Cdp * r * sp / 2;
k2Edit.Text := floatTostr(kp);
end;
// ダイバー空気抗力係数から空気抵抗計算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;
// 鉛直降下 滑空の計算はなし
// パラシュートがsinカーブで開くものとして計算
// ダイバーの姿勢変化は考慮無し
// tB: 時間 sec
// k: 身体部 空気抵抗係数 kg/m
// kb: パラシュート 空気抵抗係数 kg/m
// m: 身体部荷重(重さ) kg
// g: 重力加速度 m/s^2
procedure TForm1.downN(tB, k, KB, m, g: extended; j: integer);
var
vmax, ti : extended;
h, v: extended;
hb, vb: extended;
dv, dm, dt, mdg, fv : extended;
ndv, nmdg : extended;
i, n : integer;
mkB : extended;
maxmdg : extended;
dq, sq : extended;
Qs2 : extended;
kf : extended;
begin
series1.Clear;
series2.Clear;
series3.Clear;
series4.Clear;
Chart1.RightAxis.AutomaticMaximum := True;
case radiogroup1.ItemIndex of
0: begin
Chart1.LeftAxis.Title.Text := 'G';
Chart1.RightAxis.Title.Text := 'm';
end;
1: begin
Chart1.LeftAxis.Title.Text := 'G';
Chart1.RightAxis.Title.Text := 'm/s';
Chart1.RightAxis.AutomaticMaximum := False;
Chart1.RightAxis.Maximum := 0;
end;
2: begin
Chart1.LeftAxis.Title.Text := 'kg/m';
Chart1.RightAxis.Title.Text := 'm';
end;
3: begin
Chart1.LeftAxis.Title.Text := 'G';
Chart1.RightAxis.Title.Text := 'm';
end;
4: begin
Chart1.LeftAxis.Title.Text := 'G';
Chart1.RightAxis.Title.Text := '％';
end;
end;
n := j div 10;
h := 0; // 高さ or 距離 m
hb := h;
dt := tb / j; // Δt
fv := 0; // 速度 m/s
vmax := sqrt(m * g / k); // t∞速度 m/s
v := vmax;
vb := v; // 初速度
maxmdg := 0; // 最大加速度
Qs2 := - pi / 2; // sin 開始角度 rad
dq := pi / j; // Δ角度 Δrad
sq := 0;
for i := 0 to j + n do begin
ti := i * dt;
if i <= J then begin
sq := (1 + sin(Qs2 + dq * i)) / 2; // Q パラシュートの開き係数sin(-pi/2～pi/2)
if checkbox1.Checked = true then
mkb := k + kB * sq * sq // 空気抵抗係数 sin(Q)^2修正
else
mkb := k + kB * sq; // 空気抵抗係数 sin(Q)修正
end
else
mkb := k + kB;
vmax := sqrt(m * g / mkb); // t∞速度 m/s
dm := vmax * vmax * mkb; // t∞抗力 m.kg/s^2
dv := v * v * mkb; // 抗力 m.kg/s^2
mdg := (dv - dm) / m; // 加速度 m/s^2
fv := v - mdg * dt; // 予減速 m/s
ndv := fv * fv * mkb; // 予抗力 m.kg/s^2
nmdg := (ndv - dm) / m; // 予加速度 m/s^2
mdg := (nmdg + mdg) / 2; // 平均加速度 m/s^2
fv := v - mdg * dt; // 減速 m/s
if fv < 0 then break;
h := h + (fv + v) / 2 * dt; // 高さ + 平均速度*時間 m
v := fv; // 新速度 m/s
kf := (m - v * v * k / g) / m; // 身体部 G
if i = j then begin
hb := h;
vb := v;
end;
if radiogroup1.ItemIndex = 3 then begin
if maxmdg < mdg + kf * g then maxmdg := mdg + kf * g;
end
else
if maxmdg < mdg then maxmdg := mdg;
case radiogroup1.ItemIndex of
0: begin
series1.AddXY(ti , -mdg / g); // 加速度 G
series2.AddXY(ti , -h); // 距離 m
end;
1: begin
series1.AddXY(ti , -mdg / g); // 加速度 G
series2.AddXY(ti , -v); // 速度 m/s
end;
2: begin
series1.AddXY(ti , mkb); // 空気抵抗係数 kg/m
series2.AddXY(ti , -h); // 距離 m
end;
3: begin
series1.AddXY(ti , -mdg /g - kf); // 加速度 G
series2.AddXY(ti , -h); // 距離 m
end;
4: begin
series1.AddXY(ti , -mdg / g); // 加速度 G
if checkbox1.Checked = false then
series2.AddXY(ti , sq * 100) // パラシュート開度％
else
series2.AddXY(ti , sq * sq * 100); // パラシュート開度％
end;
end;
end;
if fv < 0 then begin
memo1.Lines.Append('新減速計算 分割数不足、時間が長い為計算出来ません。');
exit;
end;
if checkbox1.Checked = true then
memo1.Lines.Append('パラシュートの開く速度が時間のsin二乗に比例')
else
memo1.Lines.Append('パラシュートの開く速度が時間のsinに比例');
memo1.Lines.Append('t∞速度 ' + floatTostr(-vmax) + ' m/s');
memo1.Lines.Append(floattostr(tB) + 'sec後 速度 ' + floatTostr(-vb) + ' m/s');
memo1.Lines.Append(floattostr(tB) + 'sec後 高さ ' + floatTostr(-hb) + ' m');
if radiogroup1.ItemIndex = 3 then
memo1.Text := memo1.Text + 'ライザー最大負荷 ' + floatTostr(-maxmdg / g) + ' G'
else
memo1.Text := memo1.Text + '最大加速度 ' + floatTostr(-maxmdg / g) + ' G';
end;
// vh計算実行
procedure TForm1.vhBtnClick(Sender: TObject);
var
k : extended; // 空気抵抗係数 kg/m
kB : extended; // 空気抵抗係数 kg/m
m : extended; // 質量 kg
tB : extended; // 減速時間 sec
g : extended; // 重力加速度 m/sec^2
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(k2Edit.Text, kB, ch);
if ch <> 0 then begin
Memo1.Text := 'kB:空気抵抗係数に間違いがあります。';
exit;
end;
if kB <= 0 then begin
Memo1.Text := 'kB:空気抵抗係数の値は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(tbEdit.Text, tB, ch);
if ch <> 0 then begin
Memo1.Text := 'tB:Reduce speed timeに間違いがあります。';
exit;
end;
if tB < 0 then begin
Memo1.Text := 'tB:Reduce speed timeは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(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;
downN(tB, k, kB, m, g, j);
end;
procedure TForm1.FormCreate(Sender: TObject);
begin
label1.Caption := 'ρ:空気密度[kg/m^3]は' + #13#10 + ' 左の値です。';
end;
end.
free_fall_test2.zip