Программа по моделированию полета тела, брошенного под углом к горизонту (артиллерийского снаряда или ядра).
Условия задачи:
Выстрел произведен с начальной скоростью V0, под углом к горизонту α.
Считать, что снаряд имеет форму шара с радиусом r, изготовлен из материала, имеющего определенную плотность ρ.
Построить траекторию полета снаряда Y(x) ,
указать максимальную высоту полета Hk ,
дальность падения снаряда Xk и время полета tk ,
построить график скорости V(t) на отрезке [0,tk].
От себя: Целесообразно было бы добавить в исходные данные начальную высоту расположения орудия H0 , т.к. в реальности огонь может вестись с холма или из капонира (т.е. ниже уровня земли).
Таким образом, исходные данные, которые пользователь может задать на форме:
Начальная скорость V0, м/с
Начальная высота H0, м
Угол выстрела α, ° (град)
Плотность материала ρ, кг/м3
Радиус r, м
При построении математической модели условимся, что ось Оx системы координат направлена горизонтально в направлении выстрела, а ось Oy - вертикально вверх. Вектор скорости снаряда V(t) за время полета будет изменяться как по величине, так и по направлению, поэтому в модели рассматриваем его проекции на координатные оси. Горизонтальную составляющую скорости в момент времени t обозначим Vx(t), а вертикальную – Vy(t).
Пусть поверхность Земли плоская. Согласно законам механики, при сделанных предположениях движения тела в горизонтальном направлении является почти равномерным, а в вертикальном – равнозамедленным или равноускоренным с ускорением свободного падения g.
Если с силой тяжести FT все достаточно просто (она свой вектор не меняет ни по величине, ни по направлению),
то сила лобового сопротивления FC , действующая на снаряд, пропорциональна квадрату скорости движения тела.
Обозначим через FX и FY горизонтальную и вертикальную проекции вектора силы
лобового сопротивления,
причем FX/F= VX/V,
FY/F= VY/V.
Значение силы лобового сопротивления F= -b·V² (пропорционально квадрату скорости тела). Коэффициент b=0.5·C·S·ρ, где C – коэффициент лобового сопротивления (для многих задач баллистики C≈0.15), S – площадь поперечного сечения (S=πr²), ρ - плотность воздуха (ρ=1,29 кг/м3).
Решение поставленной задачи можно свести к решению системы дифференциальных уравнений
В результате численного решения задачи Коши для системы обыкновенных дифференциальных уравнений будут получены значения функций Vx(ti), Vy(ti) в точках ti=i·h, i=1,2,...; h – шаг метода.
Для решения задачи Коши буду использовать простой в реализации метод Рунге-Кутта 4-порядка, обеспечивающий достаточно высокую точность вычислений.
Для приближенного интегрирования (нахождения дальности),
после получения нового значения скорости Vx(ti)
рассчитывается X(ti)=X(ti-1)+h·Vx(ti), где
h= ti-ti-1 = const.
Кроме того, параллельно рассчитывается значение высоты
Y(ti)=Y(ti-1)+h·Vy(ti),
где h= ti-ti-1 = const по найденным значениям Vy(ti).
Когда будет получено значение Y(ti)<0, вычисления заканчиваются.
Код функции, выполняющей интегрирование по методу Рунге-Кутта, будет выглядеть так:
function ByRungeKutta:integer;
//расчет через интегрирование по методу Рунге-Кутты, формирует массивы X[i] Y[i]
var vzbool:boolean;//взлет
i:integer;
k1,k2,k3,k4,l1,l2,l3,l4,_h:Extended;
begin
Y[0]:=HH; X[0]:=0; Ymax:=Y[0];
vzbool:=true;//взлет
i:=1;
//расчет по модели и заполнение массивов
while (Y[i-1]>-0.00001) or vzbool do
begin
if i>1
then _h:=h
else _h:=h/2; //первый шаг в два раза меньше,
// чтобы уменьшить ошибку метода средних прямоугольников
t[i] :=t[i-1]+_h;
k1 := _h * fvx( Vx[i-1], Vy[i-1]);
l1 := _h * fvy( Vx[i-1], Vy[i-1]);
k2 := _h * fvx( Vx[i-1]+k1/2, Vy[i-1]+l1/2 );
l2 := _h * fvy( Vx[i-1]+k1/2, Vy[i-1]+l1/2 );
k3 := _h * fvx( Vx[i-1]+k2/2, Vy[i-1]+l2/2 );
l3 := _h * fvy( Vx[i-1]+k2/2, Vy[i-1]+l2/2 );
k4 := _h * fvx( Vx[i-1]+k3, Vx[i-1]+l3 );
l4 := _h * fvy( Vx[i-1]+k3, Vy[i-1]+l3 );
Vx[i]:=Vx[i-1] + (k1+2*k2+2*k3+k4)/6;
Vy[i]:=Vy[i-1] + (l1+2*l2+2*l3+l4)/6;
V[i] := sqrt( Vx[i]*Vx[i] + Vy[i]*Vy[i] ); //величина вектора скорости из проекций
X[i]:=X[i-1]+Vx[i]*h; //изменение координат ядра стр.3(методичка2)
Y[i]:=Y[i-1]+Vy[i]*h;
if Y[i]>Ymax then begin Ymax:=Y[i]; imaxY:=i; end //сохранение максимальной высотой
else vzbool:=false;//падение
i:=i+1;
end;
np:=i-1; //количество реальных шагов
Result:=np;
end;
где:
imaxY; //узел с макс.высотой полета
HH; //уровень расположения орудия в момент выстрела
V[],Vx[],Vy[]; //скорость и ее проекции на оси (массивы значений в узлах)
На следующем этапе интерполируем функции y(t), у(х) для нахождения точных значений дальности и времени полета, а также максимальной высоты траектории. Полученные данные используем для построения траектории полета снаряда с более мелким шагом, например h=0,02 (сек).
Для сравнения проведем расчет по упрощенной модели, известной как модель Галилея, пренебрегая сопротивлением воздуха.
скачать exe-файл для тестирования
Другие примеры на тему «РЕШЕНИЕ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ»
Другие примеры на языке «Delphi»
Поделиться в соц сетях: