Моделирование движения  спускаемого аппарата в атмосфере

Интегрированные математические пакеты

Самарский университет. Кафедра теоретической механики

Юдинцев В. В.

Модель

Рассматривается движение спускаемого аппарата в атмосфере Земли.

Условия

Вращение Земли не учитывается

Спускаемый аппарат рассматривается как материальная точка постоянной массы

Рассматривается движение спускаемого в одной плоскости

Движение спускаемого аппарата происходит под действием силы тяжести и сил аэродинамического сопротивления

Аэродинамические силы

На спускаемый аппарат при движении в атмосфере действуют аэродинамические силы и моменты. 

Аэродинамические силы и моменты зависят от следующих параметров:
- угла атаки спускаемого аппарата;
- плотности воздуха на высоте полета;
- размеров спускаемого аппарата;
- формы спускаемого аппарата.

Аэродинамические силы

Представление аэродинамических сил

(Subscript[F, D])=(Subscript[C, D])(M,α)×S×((ρ (V^2))/2)
(Subscript[F, L])=(Subscript[C, L])(M,α)×S×((ρ (V^2))/2)

S 	- площадь Миделя (характерная площадь) спускаемого аппарата
V 	- скорость набегающего потока воздуха
ρ 	- плотность воздуха (зависит от высоты полёта)
M 	- число Маха: отношение скорости набегающего потока воздуха к местной скорости  звука
(Subscript[C, D]) 	-  коэффициент лобового сопротивления, зависящий от M, α, условий обтекания 
(Subscript[C, L]) 	-  коэффициент подъемной силы сопротивления, зависящий от M, α, условий обтекания

Скоростной напор (имеет размерность давления)

q = ((ρ (V^2))/2) [Па]

Аэродинамические характеристики

Спускаемый аппарат космического корабля “Союз”

Масса: 3000 кг

Claus Weiland 
Aerodynamic Data of Space Vehicles
Springer Heidelberg New York Dordrecht London, 2014.

Аэродинамические характеристики

Плотность воздуха

Зависимость плотности воздуха от высоты может быть определена при помощи термодинамической модели атмосферы, экспериментальных данных.

StandardAtmosphereData[Quantity[10000,"Feet"],"Density",Method->"USStandardAtmosphere"]
%//QuantityMagnitude

ρ=.;
ρ[h_]:=QuantityMagnitude[StandardAtmosphereData[Quantity[h,"Meters"],"Density",Method->"USStandardAtmosphere"]];

Plot[ρ[h],{h,0,100000},PlotRange->All,PlotTheme->{"Presentation","FrameGrid"},FrameLabel->{ "h, м","ρ, кг/\!\(\*SuperscriptBox[\(м\), \(3\)]\)"} ,ImageSize->Large]

Quantity[a, "Meters"]*Quantity[b, "Newtons"]

Загрузка таблицы из файла

SetDirectory[NotebookDirectory[]];

data=Import["atm.txt"];

data=Import["atm.txt","table"];

(data[[3;;,{1,7}]])

ρ=Interpolation[%]

Plot[ρ[h],{h,0,86},PlotRange->All]

Исправить (самостоятельно)
- Таблица слишком “короткая” (до 86 км)
- Высота в километрах. В уравнениях  движения высота в метрах.

SetOptions[Plot,PlotTheme->{"Presentation", "FrameGrid"},ImageSize->Large];
SetOptions[ParametricPlot,PlotTheme->{"Presentation", "FrameGrid"},ImageSize->Large];
SetOptions[ListPlot,PlotTheme->{"Presentation", "FrameGrid"},ImageSize->Large];

Уравнения движения

Уравнения движения

q={V[t],r[t],θ[t],s[t]};

eq={
V'[t]==-g Sin[θ[t]]-(Subscript[c, d]) S(( ρ[r[t]-Rz] ((V[t])^2))/(2 m)),
V[t] θ'[t]==g((((V[t])^2)/(r[t] g))-1)Cos[θ[t]]+(Subscript[c, L])S(( ρ[r[t]-Rz] ((V[t])^2))/(2 m)),
s'[t]==((V[t])/(r[t]))Cos[θ[t]] Rz, 
r'[t]==V[t] Sin[θ[t]]
};

Вывод уравнений приведен в книге:
В. А. Ярошевский Вход в атмосферу космических летательных аппаратов
М.: Наука, 1988, стр. 24.

Параметры системы

params={(Subscript[c, d])->1.5,(Subscript[c, L])->0,S->((π (2.2^2))/4),m->3000.0,g->9.807,Rz->6371000.0,h0->100000.0}

nu={V[0]==7000,θ[0]==0,r[0]==Rz+h0,s[0]==0};

Круговая скорость на высоте h0

(Sqrt[(μ/((Rz+h0)))])/.μ->398600.4415 (10^9)//.params

params={(Subscript[c, d])->1.3,(Subscript[c, L])->0.3,S->((π (2.2^2))/4),m->3000.0,g->9.807,Rz->6371000.0,h0->100000.0,μ->398600.4415 (10^9)}
nu={V[0]==(Sqrt[μ/(Rz+h0)]),θ[0]==0,r[0]==Rz+h0,s[0]==0};

Ускорение свободного падения зависит от высоты полета. 
Например, на высоте h0:

((μ )/(((Rz+h0))^2))//.params

ρ[100000]

Результаты интегрирования

Проинтегрируем уравнения движения до 450 секунды полета

tk=1040;
sol=NDSolve[{eq,nu}//.params,q,{t,0,tk}]

Plot[(r[t]-Rz)*0.001/.params/.sol,{t,0,tk},FrameLabel->{"t, c","h, км"}]

Конечная высота

(r[t]-Rz)*0.001/.params/.sol/.t->tk

Результаты интегрирования

Траектория от дистанции полета (по поверхности Земли)

ParametricPlot[{s[t],(r[t]-Rz)}*0.001/.params/.sol,{t,0,tk},FrameLabel->{"s, км","h, км"},AspectRatio->1/GoldenRatio,PlotRange->All]

Результаты интегрирования

Траектория от дистанции полета (по поверхности Земли)

Animate[ParametricPlot[{s[t],(r[t]-Rz)}*0.001/.params/.sol,{t,0,tk},FrameLabel->{"s, км","h, км"},AspectRatio->1/GoldenRatio,Epilog->{Red,PointSize[0.02],Point[{s[t],(r[t]-Rz)}*0.001/.params/.sol/.t->tt]}],{tt,0,tk}]

Результаты интегрирования

Скорость от высоты

ParametricPlot[{(r[t]-Rz), V[t]}*0.001/.params/.sol,{t,0,tk},FrameLabel->{"h, км","V, м/c"},AspectRatio->1/GoldenRatio]

Управление процессом интегрирования

Остановка процесса численного интегрирования при достижении нулевой высоты.
Функция WhenEvent

tk=1100;
sol=NDSolve[
{
eq,nu,WhenEvent[r[t]<Rz//.params,"StopIntegration"]
}//.params,q,{t,0,tk}
]

tk=sol[[1,1,2,0,1,1,2]]

Конечная высота

r[t]-Rz/.sol/.params/.t->tk

Управление процессом интегрирования

Остановка процесса численного интегрирования при достижении нулевой высоты.
Функция WhenEvent

conditions=r[t]<Rz+#&/@Range[0, 50000, 10000]

tk=1100;
{sol,velocities}=Reap[NDSolve[
{
eq,nu,WhenEvent[Evaluate[conditions//.params],{Sow[V[t]],If[Abs[r[t]-Rz]<10//.params,"StopIntegration","CrossDiscontinuity"]}]
}//.params,q,{t,0,tk}
]
]

velocities

tk=sol[[1,1,2,0,1,1,2]]

Reap[Sum[Sow[i^2]+1,{i,10}]]

Траектория спуска

ParametricPlot[{s[t],(r[t]-Rz)}*0.001/.params/.sol,{t,0,tk},FrameLabel->{"s, км","h, км"},AspectRatio->1/GoldenRatio]

Аэродинамическая сила

Аэродинамическая сила лобового сопротивления

ParametricPlot[{(r[t]-Rz)*0.001,(Subscript[c, d]) S(( ρ[r[t]-Rz] ((V[t])^2))/2)}//.params/.sol,{t,0,tk},AspectRatio->1/GoldenRatio,FrameLabel->{"h, км","(Subscript[F, D]), Н"},PlotRange->All]

Аэродинамическая сила

Максимум аэродинамической силы на траектории спускаемого аппарата

FindMaximum[(Subscript[c, d]) S(( ρ[r[t]-Rz] ((V[t])^2))/2)/.params/.sol/.t->tt,{tt,100.0}]
Famax=%[[1]];

ParametricPlot[Flatten[{s[t],(r[t]-Rz)}]*0.001/.params/.sol,{t,0,tk},FrameLabel->{"s, км","h, км"},AspectRatio->1/GoldenRatio,ColorFunction->(ColorData["TemperatureMap"][(Subscript[c, d]) S(( ρ[r[t]-Rz] ((V[t])^2))/(2 Famax))//.params/.sol[[1]]/.t->#3]&),ColorFunctionScaling->False,PlotStyle->Thickness[0.01]]

ColorData["AvocadoColors"][0.5]

ColorData["TemperatureMap"][1]

Перегрузка

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

В рассматриваемой задаче негравитационные сила только одна -- это аэродинамическая сила

n=((Sqrt[((Subscript[F, D])^2)+((Subscript[F, L])^2)])/(m g))

((Norm[{(Subscript[c, d]),(Subscript[c, L])}S(( ρ[r[t]-Rz] ((V[t])^2))/2)])/(g m));
ParametricPlot[{(r[t]-Rz)*0.001,%}//.params/.sol,{t,0,tk},AspectRatio->1/GoldenRatio,FrameLabel->{"h, км","n, ед"},PlotRange->All]

Задание

1. Модифицируйте программу интегрирования уравнений движения спускаемого аппарата в атмосфере, добавив зависимость плотности воздуха от высоты, загружаемой из файла.

2. Модифицируйте программу интегрирования уравнений движения спускаемого аппарата в атмосфере, добавив зависимость ускорения свободного падения g от высоты. На сколько отличается дальность полета спускаемого аппарата при учете и без учета изменения ускорения свободного падения от высоты полета? 

3. Постройте зависимость максимальной перегрузки, которую испытывает экипаж спускаемого аппарата в зависимости от начального угла входа в атмосферу ((Subscript[θ, 0])) в диапазоне от минус 15 до 0 градусов.

4. Постройте зависимость дальности точки приземления от начального угла входа в атмосферу ((Subscript[θ, 0])) в диапазоне от минус 15 до 0 градусов. Нарисуйте на одном графике траектории движения спускаемого аппарата при (Subscript[θ, 0])=-(15^( o)) и при (Subscript[θ, 0])=-(2^( o)).

Пример

Функция расстояния до точки приземления в зависимости от начального угла наклона траектории.

dist[θ0_?NumberQ]:=Module[{params,nu,neq,tk,sol},
params={(Subscript[c, d])->1.3,(Subscript[c, L])->0.0,S->((π (2.2^2))/4),m->3000.0,g->9.807,Rz->6371000.0,h0->100000.0,μ->398600.4415 (10^9)};
nu={V[0]==(Sqrt[μ/(Rz+h0)]),θ[0]==θ0,r[0]==Rz+h0,s[0]==0};
neq=eq/.params;
sol=NDSolve[{eq,nu,
WhenEvent[r[t]<=Rz//.params,"StopIntegration",DetectionMethod->"Sign"]}//.params,
q,{t,0,5000}];
tk=sol[[1,1,2,0,1,1,2]];
s[t]/.sol/.t->tk
];

dist[0.0°]*0.001

Пример

Plot[dist[θθ °]*0.001,{θθ,-15,0},PlotRange->All,PlotPoints->2]

Flatten[{#*(180/π),dist[# ]*0.001}]&/@Range[-15.0°,0,1°]
ListPlot[%,FrameLabel->{"(Subscript[θ, 0]), ((...)^o)","s, км"},Joined->True]