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

Размещено на http://www.allbest.ru/

Таблица. Исходные данные и полученные результаты:

Калибр d, м	Полный вес G, кг	Вес топлива , кг	Время горения топлива t,с	Длина направляющих Sн, м	Коэффициент i	Угол подъема направляющих , 0	Эффективная скорость истечения Ue, м/с
0.15	60	12	0.6	4	1.6	45	2000

Таблица. Дополнительные параметры

Секундный массовый расход Q, кг/с	Тяга R, кг*м/с2	Коэффициент лобового сопротивления Cx	Площадь миделя S, м2
20	40000	0.26	0.01767

Таблица. Движение по направляющим

Дульная скорость Vd, м/с	Дульное время td, с	Угол наклона d, град
39.514	0.19445	45

Таблица. Активный участок траектории ЛА

Обозначение	Va ,м/с	ta, c	a , рад	Xa , m	Ya, m	M, kg
Метод Эйлера	356,305	2	0,687	274,737	234,633	100
Метод Р-К	357,048	2	0,681	270,014	228,664	100

Таблица. Пассивный участок траектории ЛА.

Обозначение	Vp,м/с	tp,c	p,рад	Xp,m	Yp,m
Метод Эйлера	46,05	285,498	-0,813	10503,017	0
Метод Р-К	45,9	286,014	-0,808	10534,668	0

1.Разработка компьютерной модели расчета траектории неуправляемого летательного аппарата

1.1 Постановка задач. Анализ и представление исходных данных

компьютерный летательный аппарат

В этой главе рассмотрены следующие задачи:

Разработка компьютерной модели расчета траектории неуправляемого ЛА, используя методы численного интегрирования Эйлера и Рунге-Кутта.

Расчет траектории пассивного участка полета неуправляемого ЛА с использованием параболической теории.

Траектория полета неуправляемого ЛА (снаряда РСЗО) состоит из трех характерных участков:

движение по направляющим;

активный участок;

пассивный участок.

Расчет дополнительных параметров:

1. Секундный массовый расход (q) рассчитан по формуле:

2. Тяга двигателя (R) рассчитана по формуле:

R=Ue*q

3. Коэффициент лобового сопротивления рассчитан по формуле:

Cx=Cx(M)*i , где Cx=Cx(M)-коэффициент лобового сопротивления эталонного ЛА. Принимаем Cx =0.3.

4. Площадь миделя ЛА рассчитана по формуле:



5. Закон изменения массовой плотности воздуха от высоты полета принимаем в виде, где ro= 1.23- плотность воздуха на уровне моря; высота y в формуле берется в метрах.

6. Закон изменения массы принимаем в виде m= m0-q*t.

1.2 Модель динамики неуправляемого ЛА

1.2.1 Движение ЛА по направляющим

При расчете движения по направляющем принимается во внимание действие на ЛА следующих сил: R- тяга, G- вес ЛА, F- реакция направляющих, которая состоит из силы трения Fтр и нормальной составляющей N. Аэродинамическими силами при расчете движения по направляющим обычно пренебрегают.

Кроме того, в идеальном случае направляющие принимаются жесткими и недеформированными. Расчетная схема движения представлена на рисунке 1.1.

Рис 1.1.

Уравнение движения ЛА на направляющих в проекции на направление скорости запишется в виде

где f- коэффициент трения скольжения, m=m0-Q*t - принимаемый в работе закон изменения массы ЛА, G=m*g. Записывая уравнение движения в нормальной форме Коши, получаем

Последовательно интегрируя уравнение с учетом соотношения ds/dt=v, находим

Полагая в формулах t=td, получаем выражение для дульной скорости Vd и длины направляющих:

Если задана длина направляющих Sh, то, решая численно уравнение, определим td. Для этого уравнение запишем в виде

Если пренебречь изменением массы при движении ЛА по направляющим, то, заменяя массу ЛА средним значением m= mср,

получаем

По этим формулам можно приближенно найти vd и td.

Примечание: 1. Коэффициент трения скольжения ЛА об направляющие f равен 0.15 (f=0.15).

2. Изменением массы при движении ЛА по направляющим пренебрегаем.

1.2.2 Активный участок полета ЛА

При расчете активного участка полета ЛА принимают следующие допущения:

1. угол атаки равен нулю.

2. тяга R постоянна.

При сделанных допущениях уравнения для расчета траектории ЛА на активном участке полета запишутся в виде следующей системы дифференциальных уравнений первого порядка:

Здесь X=0.5*Cx*Ro*S*V - сила лобового сопротивления. Остальные элементы системы были указаны выше.

1.2.3 Пассивный участок полета ЛА

При расчете пассивного участка полета ЛА также используется система (1.1). При этом:

тяга R равна нулю (R=0).

поверхность Земли представляется в виде бесконечной плоскости.

масса ЛА постоянна (m= m-mt=const).

1.3 Алгоритмы численного интегрирования на основе методов Эйлера и Рунге-Кутта

Рассмотрим метод Эйлера - линейное приближение, использующее первые два члена ряда Тейлора. Здесь искомая интегральная кривая аппроксимируется ломаной линией. Если шаг h достаточно мал, то интеграл в формуле можно вычислить используя теорему о среднем, т.е. вынося подынтегральную функцию из- под знака интеграла средним значением.

В методе Эйлера подынтегральная функция выносится при нижнем пределе интегрирования: . Это приближение геометрически соответствует движению от точки x к точке х+h по касательной к кривой y(x) в точке х. Запишем расчетные формулы метода Эйлера:

yk+1=yk+f(xk,yk)h, xk=xk-1+h

y(x0)=y0, yk=y(xk)

В усовершенствованном методе Эйлера- Коши в первом приближении полагается:

а во втором

Погрешность метода Эйлера определяется остаточным членом ряда Тейлора

т.е. R~h2 на каждом шаге вычислений. Для обеспечения сходимости шаг h следует выбирать достаточно малым. Для метода Эйлера- Коши погрешность имеет порядок h2.

Рассмотрим метод Рунге и Кутта. В основе получения вычислительных схем этого метода лежит разложение функции y(x) в ряд Тейлора с последующим преобразованием отрезка ряда к виду, не содержащему производных. На шаге h производная dy/dx=f(x,y) аппроксимируется параболой второго порядка. Здесь функция ?(x,h) определяется формулой парабол Симпсона (формула Ньютона - Котеса для трех узлов):

Рассмотрим дифференциальное уравнение при начальном условии (хА,уА). Выполним следующие операции:

По известным начальным условиям (хА,уА) определим значение производной в начальной точке А: .

Из начальной точки А проведем прямую (рис 1.2.)

и отметим значение ее ординаты в середине шага интегрирования h (точка В с координатами

).

Рис 1.2.

2) Найдем значение производной по формуле в точке В: и проведем из точки А прямую . Отметим значение ординаты этой прямой в середине шага интегрирования h (точка С с координатами ).

3) По уравнению найдем значение производной в точке С: и проведем из точки А прямую . Отметим значение ординаты этой прямой в конце шага интегрирования h (точка D с координатами ).

4) По уравнению найдем значение производной в точке D: .

В результате построений найдем значение производных в точках А, В,С и D. Отложим эти значения на графике рис 1.3. Как видно из графика, в точке с абсциссой получены два значения производной вместо одного. Это следствие приближенности метода. Примем в этой точке среднее значение производной: . Отложив на графике (рис 1.3.) ординату , получим точку М.

Рис 1.3.

Будем считать, что кривая, изображающая зависимость должна проходить через точки A,M и D. Проведем через эти три точки параболу, уравнение которой:

.

Значения коэффициентов a,b и с выбираются из условия прохождения параболы через точки А, М, и D. Коэффициент . Из уравнения параболы имеем систему:

Решив эти уравнения, найдем:

Проинтегрируем теперь уравнение параболы в пределах от x=xA до x=xA+h. Значение этого интеграла является приращением искомой функции y при изменении х на величину h. Таким образом Подставив сюда полученные выше выражения для a,b,c, после приведения подобных членов для общего случая () получим:

Как видно, приращение искомой функции на шаге h при помощи описанных построений удалось представить через значения первых производных функции в четырех точках, лежащих в пределах шага интегрирования h.

Запишем расчетные формулы метода Рунге- Кутта:

При условии существования у функции производных четвертого порядка погрешность метода является величиной порядка h5.

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

1.4 Решение контрольной задачи для одного шага интегрирования методами Эйлера и Рунге-Кутта (Проведение расчетов без ПК)

1.4.1 Начальные условия

Для решения тестового примера на 1 шаг интегрирования системы (1.1) методами Эйлера и Рунге- Кутта задаем следующие начальные условия:

1.Шаг интегрирования h принимаем равным 0.05 секунды (h=0.05 сек.);

2. V=Vg, где Vg- скорость ЛА в момент схода с направляющих, рассчитываемая по формуле:

, где

q_g - угол наклона направляющих; q_g =450=0,785 рад;

tg- дульное время, рассчитываемые по формуле:

Vg = и tg = , где

=

Vg =

3. q=q_g - угол наклона направляющих; ??g=450=0,785 рад.

4. х- начальная абсцисса, рассчитываемая по формуле:

5. у- начальная ордината, рассчитываемая по формуле:

1.4.2 Расчет системы уравнений методом Эйлера (1 шаг)

1. Нахождение коэффициентов К1 для каждого уравнения.

По определению коэффициенты К1- это правые части уравнений, вычисленные при начальных условиях.

а. плотность при начальной высоте :

б. лобовое сопротивление при начальной (дульной) скорости:

X=0,5*0.26*0.01767*1.22967*(72.9421)2=15.0288

в. Масса при начальном (дульном) времени:

m=60-20*0.1095=57.81

г. К1(1)= (40000 - 9.81*57.81* 0.785 - 15.0288) / 57.81 = 683.96 м/с2

К1(2)=-9,81*cos450/ 22,708= -0.7761 рад/с

К1(3)=22,708*cos450= 3.51304 м/с

К1(4)=22,708*sin450= 3.4857 м/с

2. Результаты одного шага интегрирования по методу Эйлера:

1.4.3 Расчет системы уравнений методом Рунге-Кутта (1 шаг)

1. Нахождение коэффициентов К1 для каждого уравнения.

По определению коэффициенты К1- это правые части уравнений, вычисленные при начальных условиях.

а. плотность при начальной высоте :

б. лобовое сопротивление при начальной (дульной) скорости:

X=0,5*0.26*0.01767*1.22967*(72.9421)2=15.0288

в. Масса при начальном (дульном) времени:

m=60-20*0.1095=57.81

г. К1(1)= 8,1224487289

К1(2)= -0,0112446838

К1(3)= 1,091383734

К1(4)= 1,0905149817

2. Нахождение коэффициентов К2 для каждого уравнения.

По определению коэффициенты К2- это правые части уравнений, вычисленные при полученных новых значениях параметров.

К2(1)= 8,1419713966

К2(2)= -0,009992507

К2(3)= 1,241945517

К2(4)= 1,2270805529

3. Нахождение коэффициентов К3 для каждого уравнения

По определению К3- это правые части уравнений, вычисленные при полученных новых значениях параметров.

К3(1)= 8,1417519898

К3(2)= -0,0099835329

К3(3)= 1,2415239848

К3(4)= 1,2282011286

4. Нахождение коэффициента К4 для всех уравнений

По определению К4- это правые части уравнений, вычисленные в конце шага

К4(1)= 8,1610849832

К4(2)= -0,0089854278

К4(3)= 1,3930436092

К4(4)= 1,364415478

5. Результаты одного шага по методу Рунге- Кутта.



1.4.4 Сравнение решений тестового примера на 1 шаг, полученных вручную (п.1.4.1 и 1.4.2) и с использованием программного комплекса

Сравнение результатов производится в соответствии с данными, представленными в таблице 1.3.:

Таблица 1.3

Обозначения	Vд	д	Xд	Yд	Vt	t	Xt	Yt
Компьютерная модель метода Эйлера	30,85671	0,785398	2.121	2.121	38,979	0,776	3,512	3,487
Ручной счет	30,85671	0,785398	2.121	2.121	38,979	0,776	3,513	3,486
Компьютерная модель метода Рунге-Кутта	30,85671	0,785398	2.121	2.121	38,998	0,775	3,363	3,349
Ручной счет	30,85671	0,785398	2.121	2.121	38,999	0,775	3,364	3,348

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

1.5 Результаты расчета траекторных параметров неуправляемого ЛА

1.5.1 Движение по направляющим

Результатами расчета траекторных параметров данного участка траектории являются значения дульной скорости (vd) и дульного времени (td).

Результаты расчета приведены в таблице 1.4.:

Таблица 1.4

Дульная скорость vg, м/с	Дульное время tg, с	Угол наклона , град
30,85671	0,19445	45

1.5.2 Активный участок траектории ЛА

Расчет траекторных параметров на конце активного участка полета ЛА произведен с помощью программного комплекса с шагом интегрирования равным 0,05

Таблица 1.5

Обозначние	Va ,м/с	ta, c	a , рад	Xa , m	Ya, m	M, kg
Метод Эйлера	356,305	2	0,687	274,737	234,633	100
Метод Р-К	357,048	2	0,681	270,014	228,664	100

1.5.3 Пассивный участок траектории ЛА

Расчет траекторных параметров на конце пассивного участка полета ЛА произведен с помощью программного комплекса с шагом интегрирования равным 0.05

Таблица 1.6.

Обозначение	tp,c	Vp,м/с	p,рад	Xp,m	Yp,m
Метод Эйлера	46,05	285,498	-0,813	10503,017	0
Метод Р-К	45,9	286,014	-0,808	10534,668	0

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

Результаты расчета траекторных параметров неуправляемого ЛА в данном пункте приводятся для характерных точек траектории, показанных на рисунке 1.11.

Рис 1.11.

где A- точка начала движения по направляющим;

B- точка конца движения по направляющим (начало активного участка);

C- точка конца активного участка (начало пассивного участка);

E- наивысшая точка траектории;

D- точка контакта снаряда с поверхностью;

Результаты расчета представлены в таблице 1.7.

Таблица 1.7

Точка	t, c	V, м/с	?, рад	Х, м	Y, м
Метод Эйлера
А	0	0	0.785	0	0
В	0,194	30,857	0.785	2.121	2.121
С	2	356,305	0,687	274,737	234,633
D	22,85	232,458	0	5514,226	2495,135
Е	46,05	285,498	-0,813	10503,017	0
Метод Р-К
A	0	0	0,785	0	0
B	0,194	30,857	0.785	2.121	2.121
C	2	357,048	0,681	270,014	228,664
D	22,75	234,079	0	5520,904	2475,650
E	45,9	286,014	-0,808	10534,668	0

Примечание: Результаты расчета получены при помощи программного комплекса

1.6 Расчет траектории пассивного участка полета ЛА с использованием параболической теории

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

Основные формулы параболической теории:

Здесь t0=0, y0=yc=228,664 м, x0=xc=270,014 м,



-параметры конца активного участка рассчитанные по данным п 1.5.2

Ниже производится расчет пассивного участка траектории ЛА с использованием параболической теории для характерных точек параболической траектории.

Для вершины траектории (точка D) время вычисляется по формуле:

Подставляя tD в формулы 3 и 4, получим X(tD)=6626.495 м Y(tD)=2804.056 м

Для точки Е (конца траектории) имеем Y(tE)=0. Тогда

Результаты расчетов представлены в сравнении с расчетами реальной траектории в таблице 1.8.

Таблица 1.8

Точка	Координата Х, м (числ/ П.Т.)	Координата Y, м (числ/ П.Т.)
C	270,014/270,014	228,664/228,664
D	5520,904/6626.495	2475,650/2804.056
E	10534,668/19615.69	0/0

Приложение 1

Графики

1. V(t) методами численного интегрирования и по параболической теории.

2. и(t) методами численного интегрирования и по параболической теории.

3. y(x) (Метод Эйлера и параболическая теория).

4. y(x) (Метод Рунге-Кутта и параболическая теория).

Рис.

Рис.



Рис.

Рис.

Рис.



Рис.

Приложение 2

Текст программы, разработанной для интегрирования методами Эйлера и Рунге-Кутта. (Программа выполнена на Visual Basic 6.0 Professional Edition)

Рис.

Option Base 1

Const r0 = 1.23

Dim r As Double

Dim y As Double

Dim x As Double

Dim V As Double

Dim Q As Double

Dim Xx As Double

Dim m As Double

Dim Yg As Double

Dim Xg As Double

Dim Cx As Double

Dim Vg As Double

Dim S As Double

Dim Mn As Double

Dim Rashod As Double

Dim Tg As Double

Dim Rs As Double

Dim T As Double

Dim h As Double

Dim Qg As Double

Dim i As Double

Dim K(4, 4) As Double

Dim str As String

Dim j As Double

Private Sub Command1_Click()

lbl.ForeColor = &HFF&

lbl.Caption = "идет расчет"

DoEvents

Yg = CDbl(Text4.Text) * Sin(CDbl(Text3.Text))

Xg = CDbl(Text4.Text) * Cos(CDbl(Text3.Text))

Cx = CDbl(Text5.Text)

Vg = CDbl(Text1.Text)

S = CDbl(Text6.Text)

Mn = CDbl(Text17.Text)

Rashod = CDbl(Text18.Text)

Tg = CDbl(Text2.Text)

Rs = CDbl(Text19.Text)

T = CDbl(Text16.Text)

h = CDbl(Text7.Text)

Qg = CDbl(Text3.Text)

Text33.Text = ""

Text34.Text = ""

Text36.Text = ""

Text40.Text = ""

Text38.Text = ""

'----------------------------------

'Активный участок траектории полета

'Эйлер

x = Xg

y = Yg

V = Vg

Q = Qg

j = Tg

For i = Tg To T Step h

r = r0 * Exp(-y / 7800)

Xx = 0.5 * Cx * r * S * V ^ 2

m = Mn - Rashod * i

V = V + h * (Rs - 9.81 * m * Sin(Q) - Xx) / m

Q = Q - h * 9.81 * Cos(Q) / V

x = x + h * V * Cos(Q)

y = y + h * V * Sin(Q)

If Option1.Value = True And Round(i, 2) = Round(j, 2) Then

j = j + 1

Text33.Text = Text33.Text & Round(i, 2) & vbCrLf

Text34.Text = Text34.Text & Round(V, 1) & vbCrLf

Text36.Text = Text36.Text & Round(Q, 3) & vbCrLf

Text40.Text = Text40.Text & Round(x, 0) & vbCrLf

Text38.Text = Text38.Text & Round(y, 0) & vbCrLf

End If

If i = Tg Then

str = "Значения коэфициентов" & vbCrLf & "на первом шаге Э" & vbCrLf

str = str & "V = " & Round(V, CLng(Text21.Text)) & vbCrLf

str = str & "Q = " & Round(Q, CLng(Text21.Text)) & vbCrLf

str = str & "X = " & Round(x, CLng(Text21.Text)) & vbCrLf

str = str & "Y = " & Round(y, CLng(Text21.Text)) & vbCrLf & vbCrLf

End If

Next i

Text8.Text = Round(V, CLng(Text21.Text))

Text9.Text = Round(Q, CLng(Text21.Text))

Text10.Text = Round(x, CLng(Text21.Text))

Text11.Text = Round(y, CLng(Text21.Text))

'-----------------------------

'Рунге-Кутт

x = Xg

y = Yg

V = Vg

Q = Qg

j = Tg

For i = Tg To T Step h

r = r0 * Exp(-y / 7800)

Xx = 0.5 * Cx * r * S * V ^ 2

m = Mn - Rashod * i

K(1, 1) = h * (Rs - 9.81 * m * Sin(Q) - Xx) / m

K(1, 2) = -9.81 * h * Cos(Q) / V

K(1, 3) = h * V * Cos(Q)

K(1, 4) = h * V * Sin(Q)

r = r0 * Exp(-(y + 0.5 * K(1, 4)) / 7800)

Xx = 0.5 * Cx * r * S * (V + 0.5 * K(1, 1)) ^ 2

m = Mn - Rashod * (i + h / 2)

K(2, 1) = h * (Rs - 9.81 * m * Sin(Q + K(1, 2) / 2) - Xx) / m

K(2, 2) = -9.81 * h * Cos(Q + K(1, 2) / 2) / (V + K(1, 1) / 2)

K(2, 3) = h * (V + K(1, 1) / 2) * Cos(Q + K(1, 2) / 2)

K(2, 4) = h * (V + K(1, 1) / 2) * Sin(Q + K(1, 2) / 2)

r = r0 * Exp(-(y + 0.5 * K(2, 4)) / 7800)

Xx = 0.5 * Cx * r * S * (V + 0.5 * K(2, 1)) ^ 2

m = Mn - Rashod * (i + h / 2)

K(3, 1) = h * (Rs - 9.81 * m * Sin(Q + K(2, 2) / 2) - Xx) / m

K(3, 2) = -9.81 * h * Cos(Q + K(2, 2) / 2) / (V + K(2, 1) / 2)

K(3, 3) = h * (V + K(2, 1) / 2) * Cos(Q + K(2, 2) / 2)

K(3, 4) = h * (V + K(2, 1) / 2) * Sin(Q + K(2, 2) / 2)

r = r0 * Exp(-(y + K(3, 4)) / 7800)

Xx = 0.5 * Cx * r * S * (V + K(3, 1)) ^ 2

m = Mn - Rashod * (i + h)

K(4, 1) = h * (Rs - 9.81 * m * Sin(Q + K(3, 2)) - Xx) / m

K(4, 2) = -9.81 * h * Cos(Q + K(3, 2)) / (V + K(3, 1))

K(4, 3) = h * (V + K(3, 1)) * Cos(Q + K(3, 2))

K(4, 4) = h * (V + K(3, 1)) * Sin(Q + K(3, 2))

V = V + (K(1, 1) + K(4, 1) + 2 * (K(2, 1) + K(3, 1))) / 6

Q = Q + (K(1, 2) + K(4, 2) + 2 * (K(2, 2) + K(3, 2))) / 6

x = x + (K(1, 3) + K(4, 3) + 2 * (K(2, 3) + K(3, 3))) / 6

y = y + (K(1, 4) + K(4, 4) + 2 * (K(2, 4) + K(3, 4))) / 6

If Option2.Value = True And Round(i, 2) = Round(j, 2) Then

j = j + 1

Text33.Text = Text33.Text & Round(i, 2) & vbCrLf

Text34.Text = Text34.Text & Round(V, 1) & vbCrLf

Text36.Text = Text36.Text & Round(Q, 3) & vbCrLf

Text40.Text = Text40.Text & Round(x, 0) & vbCrLf

Text38.Text = Text38.Text & Round(y, 0) & vbCrLf

End If

If i = Tg Then

str = str & "Значения коэффициентов" & vbCrLf & "на первом шаге Р-К" & vbCrLf

str = str & "K(1,1) = " & Round(K(1, 1), CLng(Text21.Text)) & vbCrLf

str = str & "K(1,2) = " & Round(K(1, 2), CLng(Text21.Text)) & vbCrLf

str = str & "K(1,3) = " & Round(K(1, 3), CLng(Text21.Text)) & vbCrLf

str = str & "K(1,4) = " & Round(K(1, 4), CLng(Text21.Text)) & vbCrLf

str = str & "K(2,1) = " & Round(K(2, 1), CLng(Text21.Text)) & vbCrLf

str = str & "K(2,2) = " & Round(K(2, 2), CLng(Text21.Text)) & vbCrLf

str = str & "K(2,3) = " & Round(K(2, 3), CLng(Text21.Text)) & vbCrLf

str = str & "K(2,4) = " & Round(K(2, 4), CLng(Text21.Text)) & vbCrLf

str = str & "K(3,1) = " & Round(K(3, 1), CLng(Text21.Text)) & vbCrLf

str = str & "K(3,2) = " & Round(K(3, 2), CLng(Text21.Text)) & vbCrLf

str = str & "K(3,3) = " & Round(K(3, 3), CLng(Text21.Text)) & vbCrLf

str = str & "K(3,4) = " & Round(K(3, 4), CLng(Text21.Text)) & vbCrLf

str = str & "K(4,1) = " & Round(K(4, 1), CLng(Text21.Text)) & vbCrLf

str = str & "K(4,2) = " & Round(K(4, 2), CLng(Text21.Text)) & vbCrLf

str = str & "K(4,3) = " & Round(K(4, 3), CLng(Text21.Text)) & vbCrLf

str = str & "K(4,4) = " & Round(K(4, 4), CLng(Text21.Text)) & vbCrLf

str = str & "V = " & Round(V, CLng(Text21.Text)) & vbCrLf

str = str & "Q = " & Round(Q, CLng(Text21.Text)) & vbCrLf

str = str & "X = " & Round(x, CLng(Text21.Text)) & vbCrLf

str = str & "Y = " & Round(y, CLng(Text21.Text))

Text20.Text = str

End If

Next i

Text12.Text = Round(V, CLng(Text21.Text))

Text13.Text = Round(Q, CLng(Text21.Text))

Text14.Text = Round(x, CLng(Text21.Text))

Text15.Text = Round(y, CLng(Text21.Text))

'----------------------------------

'Пассивный участок траектории полета

'Эйлер

V = CDbl(Text8.Text)

Q = CDbl(Text9.Text)

x = CDbl(Text10.Text)

y = CDbl(Text11.Text)

i = T

m = Mn - Rashod * T

j = T

Do While y > 0

r = r0 * Exp(-y / 7800)

Xx = 0.5 * Cx * r * S * V ^ 2

V = V - h * (9.81 * m * Sin(Q) + Xx) / m

Q = Q - h * 9.81 * Cos(Q) / V

x = x + h * V * Cos(Q)

y = y + h * V * Sin(Q)

If Option1.Value = True And Round(i, 2) = Round(j, 2) Then

j = j + 1

Text33.Text = Text33.Text & Round(i, 2) & vbCrLf

Text34.Text = Text34.Text & Round(V, 1) & vbCrLf

Text36.Text = Text36.Text & Round(Q, 3) & vbCrLf

Text40.Text = Text40.Text & Round(x, 0) & vbCrLf

Text38.Text = Text38.Text & Round(y, 0) & vbCrLf

End If

If Abs(Q) < 0.001 Then

str = "Значения параметров" & vbCrLf & "в наивысшей точке траектории" & vbCrLf & "методом Эйлера" & vbCrLf

str = str & "t = " & Round(i, 2) & vbCrLf

str = str & "V = " & Round(V, CLng(Text21.Text)) & vbCrLf

str = str & "Q = 0" & vbCrLf

str = str & "X = " & Round(x, CLng(Text21.Text)) & vbCrLf

str = str & "Y = " & Round(y, CLng(Text21.Text)) & vbCrLf & vbCrLf

End If

i = i + h

Loop

Text25.Text = Round(V, CLng(Text21.Text))

Text24.Text = Round(Q, CLng(Text21.Text))

Text23.Text = Round(x, CLng(Text21.Text))

Text22.Text = Round(y, CLng(Text21.Text))

Text31.Text = Round(i, 2)

'-----------------------------

'Рунге-Кутт

V = CDbl(Text12.Text)

Q = CDbl(Text13.Text)

x = CDbl(Text14.Text)

y = CDbl(Text15.Text)

i = T

j = T

m = Mn - Rashod * T

Do While y > 0

r = r0 * Exp(-y / 7800)

Xx = 0.5 * Cx * r * S * V ^ 2

K(1, 1) = h * (-9.81 * m * Sin(Q) - Xx) / m

K(1, 2) = -9.81 * h * Cos(Q) / V

K(1, 3) = h * V * Cos(Q)

K(1, 4) = h * V * Sin(Q)

r = r0 * Exp(-(y + 0.5 * K(1, 4)) / 7800)

Xx = 0.5 * Cx * r * S * (V + 0.5 * K(1, 1)) ^ 2

K(2, 1) = h * (-9.81 * m * Sin(Q + K(1, 2) / 2) - Xx) / m

K(2, 2) = -9.81 * h * Cos(Q + K(1, 2) / 2) / (V + K(1, 1) / 2)

K(2, 3) = h * (V + K(1, 1) / 2) * Cos(Q + K(1, 2) / 2)

K(2, 4) = h * (V + K(1, 1) / 2) * Sin(Q + K(1, 2) / 2)

r = r0 * Exp(-(y + 0.5 * K(2, 4)) / 7800)

Xx = 0.5 * Cx * r * S * (V + 0.5 * K(2, 1)) ^ 2

K(3, 1) = h * (-9.81 * m * Sin(Q + K(2, 2) / 2) - Xx) / m

K(3, 2) = -9.81 * h * Cos(Q + K(2, 2) / 2) / (V + K(2, 1) / 2)

K(3, 3) = h * (V + K(2, 1) / 2) * Cos(Q + K(2, 2) / 2)

K(3, 4) = h * (V + K(2, 1) / 2) * Sin(Q + K(2, 2) / 2)

r = r0 * Exp(-(y + K(3, 4)) / 7800)

Xx = 0.5 * Cx * r * S * (V + K(3, 1)) ^ 2

K(4, 1) = h * (-9.81 * m * Sin(Q + K(3, 2)) - Xx) / m

K(4, 2) = -9.81 * h * Cos(Q + K(3, 2)) / (V + K(3, 1))

K(4, 3) = h * (V + K(3, 1)) * Cos(Q + K(3, 2))

K(4, 4) = h * (V + K(3, 1)) * Sin(Q + K(3, 2))

V = V + (K(1, 1) + K(4, 1) + 2 * (K(2, 1) + K(3, 1))) / 6

Q = Q + (K(1, 2) + K(4, 2) + 2 * (K(2, 2) + K(3, 2))) / 6

x = x + (K(1, 3) + K(4, 3) + 2 * (K(2, 3) + K(3, 3))) / 6

y = y + (K(1, 4) + K(4, 4) + 2 * (K(2, 4) + K(3, 4))) / 6

If Option2.Value = True And Round(i, 2) = Round(j, 2) Then

j = j + 1

Text33.Text = Text33.Text & Round(i, 2) & vbCrLf

Text34.Text = Text34.Text & Round(V, 1) & vbCrLf

Text36.Text = Text36.Text & Round(Q, 3) & vbCrLf

Text40.Text = Text40.Text & Round(x, 0) & vbCrLf

Text38.Text = Text38.Text & Round(y, 0) & vbCrLf

End If

If Abs(Q) < 0.001 Then

str = str & "Значения параметров" & vbCrLf & "в наивысшей точке траектории" & vbCrLf & "методом Рунге-Кутта" & vbCrLf

str = str & "t = " & Round(i, 2) & vbCrLf

str = str & "V = " & Round(V, CLng(Text21.Text)) & vbCrLf

str = str & "Q = 0" & vbCrLf

str = str & "X = " & Round(x, CLng(Text21.Text)) & vbCrLf

str = str & "Y = " & Round(y, CLng(Text21.Text))

Text30.Text = str

End If

i = i + h

Loop

Text29.Text = Round(V, CLng(Text21.Text))

Text28.Text = Round(Q, CLng(Text21.Text))

Text27.Text = Round(x, CLng(Text21.Text))

Text26.Text = Round(y, CLng(Text21.Text))

Text32.Text = Round(i, 2)

lbl.ForeColor = &HC000&

lbl.Caption = "Расчет выполнен"

End Sub

Private Sub Command2_Click()

End

End Sub

Список литературы

1. Бахвалов Н.С. Численные методы. т.1., М.: Наука, 1975.

2. Волков Е.А. Численные методы. М.: Наука, 1982

3. Шалыгин А.С. и др. Автоматизация расчета траекторий ЛА: Учебн. пособие. - Л., 1990.

4. http://a5.bstu.spb.su

Размещено на Allbest.ru