0 / 0 / 0
Регистрация: 11.12.2013
Сообщений: 16
1

Расчет траектории снаряда с учетом неоднородности плотности воздуха по высоте

18.11.2015, 00:43. Показов 4449. Ответов 12

Author24 — интернет-сервис помощи студентам
Здравствуйте, в универе задали решить задачу: нужно найти угол возвышения ствола, при котором дальность выстрела будет максимальной. Сначала написал программу без учета неоднородной плотности воздуха по высоте:
Matlab M
1
2
3
4
function u=airpoint(~,v,m)
g=10; ro=1; d=7.62*10^-3; s=pi*(d/2)^2; 
k=ro*s/2/m;
u=[-k*sqrt(v(1)^2+v(2)^2)*v(1);-g-k*sqrt(v(1)^2+v(2)^2)*v(2)];
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
function dynairpoin_max()
mas = 0.0079;
n=1; %стрельба начинается с 1грд
d=0;
Vin=806; %начальная скорость
tt=sqrt(Vin);
for angle=n:90
a=Vin*cosd(angle);
b=Vin*sind(angle);
[t,h]=ode113(@airpoint,[0,tt],[a,b], [], mas);
vx=h(:,1); vy=h(:,2);
m=length(t);
x=0; y=0;
x0=0; y0=0;
x(1)=x0; y(1)=y0;
for i=2:m
    x(i)=x(i-1)+vx(i-1)*(t(i)-t(i-1));
    y(i)=y(i-1)+vy(i-1)*(t(i)-t(i-1));
end;
ym=abs(y);
ym(1:ceil(m/2))=mean(ym);
[~,k] = min(ym);
distance = x(k);
disp(distance)
dd(n)=distance;
n=n+1;
end
[d,n] = max(dd);
fprintf('При скорости %.0fм/с, м дальность выстрела(%.1fм) будет при угле %.0f град', Vin, d, n);
disp('.');
a=Vin*cosd(n);
b=Vin*sind(n);
[t,h]=ode113(@airpoint,[0,tt],[a,b], [], mas);
vx=h(:,1); vy=h(:,2);
m=length(t);
x=0; y=0;
x(1)=0; y(1)=0;
for i=2:m
    x(i)=x(i-1)+vx(i-1)*(t(i)-t(i-1));
    y(i)=y(i-1)+vy(i-1)*(t(i)-t(i-1));
end;
figure
plot(x,y,'k.',x,y);
grid on;
axis([0, d, 0, d]);
А вот с плотностью воздуха пока не получается. Мне кажется, что в процессе работы решателя нужно чтобы он находил высоту с каждой итерацией расчета и результат учитывал при расчете следующей итерации. Как это реализовать и есть ли другие алгоритмы для решения этой задачи?
0
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
18.11.2015, 00:43
Ответы с готовыми решениями:

Нахождение траектории движения (с учетом силы сопр. воздуха) снаряда
Здравствуйте, помогите решить задачу! Только по этой задаче надо найти не то, что указано в...

Расчет траектории движении тела с учетом сопротивления воздуха
Подобная задача, без учета сопротивления среды была только, что мною решена. если поможет, могу...

Расчет траектории полета осколков снаряда
снаряд, вылетевший из орудия под углом 60 градусов со скоростью v0. разорвался в верхней точке на...

Интерполяция значения плотности воздуха
Подскажите, как в C# с помощью интерполяции вычислить значение плотности воздуха если известна...

12
Эксперт по математике/физике
3390 / 1913 / 571
Регистрация: 09.04.2015
Сообщений: 5,365
18.11.2015, 07:31 2
Составь 2 дифф уравнения второй степени по расстоянию и по высоте и решай их солвером.
А угол возвышения ствола в 90 градусов, это здорово
1
373 / 343 / 42
Регистрация: 14.07.2015
Сообщений: 2,890
18.11.2015, 10:51 3
Хорошая задача, при определенной начальной скорости ядро и на орбиту сможет выйти (дальность полета бесконечна). Интересно будет рассчитать такой случай
1
0 / 0 / 0
Регистрация: 11.12.2013
Сообщений: 16
18.11.2015, 12:12  [ТС] 4
SSC, можете помочь их составить? Я диффуры прохо шарю(
90 градусов - это отлично
bobah16, да, интересный случай.
0
373 / 343 / 42
Регистрация: 14.07.2015
Сообщений: 2,890
18.11.2015, 13:33 5
super_bum, тут не диффуры надо знать, а физику. В векторном виде ваше уравнение выглядит как: https://www.cyberforum.ru/cgi-bin/latex.cgi?m\vec{a}=\sum \vec{{F}_{i}}. Ускорение это вторая производная координаты по времени, в правой части силы, действующие на ядро, т.е. сила тяжести и сила трения об воздух. Расписав проекции на вертик. и гор. ось, получим 2 дифф. уравнения. Ну и само собой неоднородность плотности воздуха входит в силу трения.
Про силу сопротивления воздуха можно почитать тут -http://studopedia.org/3-139617.html
0
Эксперт по математике/физике
3390 / 1913 / 571
Регистрация: 09.04.2015
Сообщений: 5,365
18.11.2015, 15:59 6
Если появится премя, то заготовку для MATLAB сделаю.

По поводу выхода на орбиту, то расчеты надо вести в полярных координатах.
Если считать в прямоугольных, то хоть с третьей космической запускай - все равно назад Земля притянет
0
0 / 0 / 0
Регистрация: 11.12.2013
Сообщений: 16
18.11.2015, 16:39  [ТС] 7
Написал что-то такое. Правильно?
Matlab M
1
2
3
4
5
6
function u=airpoint_2por(~,x)
g=10; ro=1; d=7.62*10^-3; s=pi*(d/2)^2; m=0.0079; k=s/2/m;
u=[x(3);x(4);...
    -k*ro*exp(-29*10^-3*9.8*(x(2)/(8.31*300)))*sqrt(x(3)^2+x(4)^2)*x(3);...
    -g-k*sqrt(x(3)^2+x(4)^2)*x(4)];
%ro*exp(-29*10^-3*9.8*(x(2)/(8.31*300))) - зависимость плотности от высоты
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
function dynairpoint_2por()
Vin = 806;
angle = 45;
a=Vin*cosd(angle);
b=Vin*sind(angle);
[t,h]=ode45(@airpoint_2por,[0,30],[0,0,a,b],...
    odeset('Reltol', 10e-9 ));
m=length(t);
x=h(:,1); y=h(:,2);
yabs = abs(y);
yabs(1:ceil(m/2))=mean(yabs);
[~,k] = min(yabs);
dist = x(k);
disp(dist);
figure
plot(x,y,'k.',x,y);
grid on
0
373 / 343 / 42
Регистрация: 14.07.2015
Сообщений: 2,890
18.11.2015, 22:34 8
Цитата Сообщение от SSC Посмотреть сообщение
По поводу выхода на орбиту, то расчеты надо вести в полярных координатах.
Если считать в прямоугольных, то хоть с третьей космической запускай - все равно назад Земля притянет
Не совсем так. Если в прямоугольных координатах учесть кривизну Земли, то все будет ок. Решение не зависит от выбора с.к. Если это так, где-то ошибка.
Да и при далеких полетах также надо учитывать силу Кориолиса и вращение Земли вокруг своей оси.

Добавлено через 5 часов 26 минут
Цитата Сообщение от super_bum Посмотреть сообщение
Правильно?
В последнем уравнении также необходимо учесть, что плотность изменяется. Неплохо бы еще учесть и то, что плотность зависит не только от высоты, но и от температуры. Формулу такой зависимости несложно найти. Хотя, если высоко не пулять, то и так сойдет.
P.S. У меня примерно 23 градуса получилось
0
0 / 0 / 0
Регистрация: 11.12.2013
Сообщений: 16
19.11.2015, 02:32  [ТС] 9
Цитата Сообщение от bobah16 Посмотреть сообщение
В последнем уравнении также необходимо учесть, что плотность изменяется.
А слона то я и не заметил

Добавлено через 46 минут
В итоге, код такой:
Попробовал сначала в airpoint_2por(~,x) подставлять в систему не всю формулу расчета плотности, а только идентификатор с готовой формулой, Но почему-то матлаб возвращал нули вместо ответа.
Matlab M
1
2
3
4
5
function u=airpoint_2por(~,x)
g=9.8; d=7.62*10^-3; s=pi*(d/2)^2; m=0.0079; k=s/2/m;
u=[x(3);x(4);...
    -k*((10^5*(1+(-0.0065*x(2))/288)^((-g*29*10^-3)/(8.31*(-0.0065))))*29*10^-3)/(8.31*(288+0.0065*x(2)))*sqrt(x(3)^2+x(4)^2)*x(3);...
    -g-k*((10^5*(1+(-0.0065*x(2))/288)^((-g*29*10^-3)/(8.31*(-0.0065))))*29*10^-3)/(8.31*(288+0.0065*x(2)))*sqrt(x(3)^2+x(4)^2)*x(4)];
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
function dynairpoint_2por()
tic
Vin = 806;
tol = 10^-4;
for angle=1:90;
a=Vin*cosd(angle);
b=Vin*sind(angle);
[t,h]=ode45(@airpoint_2por,[0,sqrt(Vin)],[0,0,a,b],...
    odeset('Reltol', tol));
m=length(t);
x=0; y=0;
x=h(:,1); y=h(:,2);
yabs = abs(y);
yabs(1:ceil(m/2))=mean(yabs);
[~,k] = min(yabs);
dist(angle) = x(k);
end;
[d,angle]=max(dist);
disp(d);
disp(angle);
a=Vin*cosd(angle);
b=Vin*sind(angle);
[~,h]=ode45(@airpoint_2por,[0,70],[0,0,a,b],...
    odeset('Reltol', tol));
x=0; y=0;
x=h(:,1); y=h(:,2);
figure
plot(x,y,'k.',x,y);
grid on
axis([0, d, 0, d]);
toc
Добавлено через 49 секунд
А ведь еще и ускорение свободного падения меняется...
0
Эксперт по математике/физике
3390 / 1913 / 571
Регистрация: 09.04.2015
Сообщений: 5,365
19.11.2015, 07:28 10
По поводу прямоугольной системы координат я не совсем корректно написал.
Более точно - начало координат прямоугольной системы надо помещать в центр Земли,
и главное направление силы тяжести всегда направлять к началу координат.

Писать всю формулу в одну строку крайне не рационально (стр.4,5).
1. Очень трудно найти ошибку
2. Ряд вычислений повторяется, а это прямая потеря производительности

Предлагаю после стр.2 провести предварительные вычисления текущих вспомогательных величин, а дальше уже более упрощенные формулы с лучьшим сохранением физического смысла.

Мне совершенно не понятен физический смысл следующих фрагментов формул
....*sqrt(x(3)^2+x(4)^2)*x(3)
....*sqrt(x(3)^2+x(4)^2)*x(4)
Единицы изменения этих фрагментов м^2/c^2
Нам надо определить силу трения, действующую в направлении оси координат
Я так думаю что надо сначала определить полную силу трения, которая пропорциональна полной скорости ( sqrt(...) ), или какой-нибудь степени этой скорости.
А потом, умножив на sin или cos угла, найти соответствующую составляющую силы трения в каждом из направлений оси координат
(Умножение модуля скорости на координатную составляющую дает, по моему, несколько другой результат)
1
373 / 343 / 42
Регистрация: 14.07.2015
Сообщений: 2,890
19.11.2015, 09:12 11
Цитата Сообщение от SSC Посмотреть сообщение
Более точно - начало координат прямоугольной системы надо помещать в центр Земли,
и главное направление силы тяжести всегда направлять к началу координат.
Да, так будет удобно, но можно н.к. поместить куда угодно, результат опять же от этого не зависит. Что, если с Луны смотреть на пулю, она по-другому двигаться будет? Естественно формулы будут посложнее. Это ведь физика Не важно какая с.к. выбрана и где ее начало, просто выбор естественной с.к. для конкретной задачи сильно облегчает решение.
Цитата Сообщение от SSC Посмотреть сообщение
Мне совершенно не понятен физический смысл следующих фрагментов формул
....*sqrt(x(3)^2+x(4)^2)*x(3)
....*sqrt(x(3)^2+x(4)^2)*x(4)
Единицы изменения этих фрагментов м^2/c^2
Это https://www.cyberforum.ru/cgi-bin/latex.cgi?{V}^{2}*cos(angle) и https://www.cyberforum.ru/cgi-bin/latex.cgi?{V}^{2}*sin(angle). Единицы измерения соответственно верные. Квадрат скорости взялся из силы трения.
https://www.cyberforum.ru/cgi-bin/latex.cgi?{V}^{2}={{V}_{x}}^{2}+{{V}_{y}}^{2}
https://www.cyberforum.ru/cgi-bin/latex.cgi?cos(a)={V}_{x}/V={V}_{x}/(sqrt({V}^{2}))={V}_{x}/(sqrt({{V}_{x}}^{2}+{{V}_{y}}^{2}))
https://www.cyberforum.ru/cgi-bin/latex.cgi?{V}^{2}*cos(a)=({{V}^{2}}_{x}+{{V}^{2}}_{y})*{V}_{x}/sqrt({{V}^{2}}_{x}+{{V}^{2}}_{y})=sqrt({{V}^{2}}_{x}+{{V}^{2}}_{y})*{V}_{x}

Ну и x(1)=x, x(2)=y, x(3)=Vx, x(4)=Vy

Цитата Сообщение от SSC Посмотреть сообщение
Нам надо определить силу трения, действующую в направлении оси координат
Я так думаю что надо сначала определить полную силу трения, которая пропорциональна полной скорости ( sqrt(...) ), или какой-нибудь степени этой скорости.
А потом, умножив на sin или cos угла, найти соответствующую составляющую силы трения в каждом из направлений оси координат
Да, все так, вот только угол выражен через скорость. Хотя в функцию можно и угол передать и не выражать его никак. Тут кому как больше нравится.

Цитата Сообщение от SSC Посмотреть сообщение
Писать всю формулу в одну строку крайне не рационально (стр.4,5).
1. Очень трудно найти ошибку
2. Ряд вычислений повторяется, а это прямая потеря производительности
Полностью согласен.

Добавлено через 2 минуты
Цитата Сообщение от super_bum Посмотреть сообщение
А ведь еще и ускорение свободного падения меняется...
Ну так не проблема и это учесть Сложнее будет, когда g будет зависеть и от высоты и от широты

На самом деле, задача состоит не в расчете траектории пули, а в нахождении угла при максимальной дальности. Это уже вариационная задача, которую также можно численно решить и найти угол за 1 раз, а не строить кучу траекторий и примерно определить угол
2
Эксперт по математике/физике
3390 / 1913 / 571
Регистрация: 09.04.2015
Сообщений: 5,365
19.11.2015, 09:19 12
Цитата Сообщение от bobah16 Посмотреть сообщение
Это *cos(angle) и *sin(angle).
Теперь понятно, это частное решение для квадратичной зависимости.
К тому же есть доп. плюс, отслеживаются знаки при движении вниз
0
0 / 0 / 0
Регистрация: 11.12.2013
Сообщений: 16
19.11.2015, 11:36  [ТС] 13
Цитата Сообщение от SSC Посмотреть сообщение
Писать всю формулу в одну строку крайне не рационально (стр.4,5).
1. Очень трудно найти ошибку
2. Ряд вычислений повторяется, а это прямая потеря производительности
Да, я понимаю это, но попробовал сначала написать вот так:
Matlab M
1
2
3
4
5
6
7
8
9
function u=airpoint_2por(~,x)
g=9.8; d=7.62*10^-3; s=pi*(d/2)^2; m=0.0079; k=s/2/m;
p0=10^5; T0=288; L=-0.0065; R=8.31447; M=28.9644*10^-3;
p=p0*(1+(L*x(2))/T0)^((-g*M)/(R*L));
T=T0*L*x(2);
ro=(p*M)/(R*T);
u=[x(3);x(4);...
    -k*ro*sqrt(x(3)^2+x(4)^2)*x(3);...
    -g-k*ro*sqrt(x(3)^2+x(4)^2)*x(4)];
Решатель возвращал нули вместо ответа.

Добавлено через 6 минут
Цитата Сообщение от SSC Посмотреть сообщение
Мне совершенно не понятен физический смысл следующих фрагментов формул
....*sqrt(x(3)^2+x(4)^2)*x(3)
....*sqrt(x(3)^2+x(4)^2)*x(4)
Это https://www.cyberforum.ru/cgi-bin/latex.cgi?\sqrt{{Vx}^{2}+{Vy}^{2}}
Где Vx - горизонтальная скорость, Vy - вертикальная скорость. Vx - 3 столбец массива, Vy - 4 столбец.
0
19.11.2015, 11:36
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
19.11.2015, 11:36
Помогаю со студенческими работами здесь

Напечатать таблицу зависимости плотности воздуха от высоты
Плотность воздуха убывает с высотой по закону р = р0 е-bz, где р — плотность на высоте b...

Определить плотность воздуха на высоте х
Плотность воздуха убывает с высотой по закону: p=p0*e^(-hz). Считая, что p0=1,29 кг/м^(3),...

Определить на какой высоте плотность воздуха будет меньше
Написать программу с оператором do..while Я конечно понимаю что у меня много ошибок, пытался...

На какой высоте давление воздуха будет в n раз меньше?
Давление воздуха у поверхности Земли равно р0. На какой высоте давление воздуха будет в n раз...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
13
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin
Copyright ©2000 - 2024, CyberForum.ru