Пример решения задачи.

Задание. Торпеда, снабженная разгонным двигателем, нацеливается с лежащей на дне подводной лодки на поражение движущегося надводного корабля. Пуск торпеды производится в момент прохождения корабля над лодкой. Исследовать связь между глубиной залегания лодки, временем поражения цели, условиями поражения (энерговооруженностью двигателя и параметрами старта).

Схема задачи приведена на рисунке. Принимается, что торпеда запускается в момент прохождения над ней корабля. Запуск производится с начальной скоростьюпод углом_н к горизонтальной оси. Глубина залегания ракеты составляет L, скорость движения корабля равна . Принимается, что тяга двигателя постоянна, равнаи направлена по касательной к траектории движения торпеды по ходу движения. На торпеду действует сила сопротивления среды, равная по модулю, где - модуль текущей скорости торпеды, которая направлена по касательной к траектории движения торпеды против хода ее движения. Коэффициенты сопротивления линейного и квадратного слагаемых будем вычислять по формулам:

,

где  - динамическая вязкость среды,  - плотность среды, S – площадь поперечного сечения торпеды, r – ее радиус поперечного сечения, с – коэффициент лобового сопротивления. Если поперечное сечение имеет некруглую форму, то в качестве его радиуса берется приведенная величина:

.

На торпеду при движении действует сила тяжести с ускорением свободного падения g, направленная вертикально вниз. Угол  наклона касательной к траектории движения связан с компонентами вектора скорости _х и _у следующими соотношениями:

,

которые используются для проектирования исходных уравнений движения:

на координатные оси. Здесь - вектор перемещения торпеды.

Тогда с учетом сделанных предположений система уравнений принимает вид:

(П1)

Для отладки программы принимаются следующие значения исходных данных: плотность воды  = 1000 кг/м³, динамическая вязкость среды  = 1,002 Нс/м², коэффициент лобового (аэродинамического) сопротивления с = 0,1, ускорение свободного падения g = 9,81 м/с². Принимаются следующие параметры торпеды и параметров старта: начальная скорость торпеды v_н = 1 м/с, угол старта  = 80, радиус поперечного сечения торпеды r = 0,3 м, масса торпеды m = 250 кг, сила тяги двигателя F = 30000 Н, глубина залегания торпеды L = 200 м, скорость движения корабля v_к = 10 м/с. Площадь поперечного сечения, представляющего собой круг, вычисляется по формуле S = r² . Расстояние смещения корабля после старта торпеды вычисляется по формуле L_к = v_к t. Для полного определения задачи Коши задаем начальные условия следующим образом:

v_x (0) = v_н cos _н , v_y (0) = v_н sin _н , x(0) = 0, y(0) = 0. (П2)

Вычисления прекращаются, когда координата y достигнет величины глубина залегания торпеды L.

Для численного интегрирования системы обыкновенных дифференциальных уравнений будем использовать классический метод Рунге-Кутта 4-го порядка. Постановка задачи Коши для системы n дифференциальных уравнений первого порядка имеет вид - требуется найти решение системы:

(П3)

при начальных условиях

y₁ (x₀ ) = y₁₀ , y₂ (x₀ ) = y₂₀ , … , y_n (x₀ ) = y_n0 .

Если правые части дифференциальных уравнений системы – непрерывны вместе со своими частными производными по переменным y_i в некоторой области D (n+1)-мерного пространства, то для любой точки (x₀ , y₀₁ , …, y_n ) D система имеет единственное решение, удовлетворяющее начальным условиям

y₁ (x₀ ) = y₁₀ , y₂ (x₀ ) = y₂₀ , … , y_n (x₀ ) = y_n0 . (П4)

Введем векторные обозначения:

Тогда задача Коши в векторной форме имеет вид:

. (П5)

Численное решение задачи Коши состоит в том, что на отрезке [a, b] требуется получить приближенные значения координат вектора в узлах сеткиx_i , i = 1, 2, … , m. Обозначим вектор, аппроксимирующий решение, через , а его координаты – черезy_ki , k = 1, 2, … , n, , i = 1, 2, … , m так, что y_ki = y_k (x_i ) или

. (П6)

Будем искать решение на равномерной сетке с шагом h = (b – a) /m.

Погрешность численного метода оценивается величиной , гдеd_i – погрешность решения на сетке с шагом h в точке x_i :

. (П7)

Практически погрешность в точке x_i оценивают по формуле Рунге. Пусть: . – значения численного решения в точкеx_i , полученные для шагов h и h/2 соответственно. тогда погрешность d_i в точке x_i для вычислений с шагом h/2 выражается приближенным равенством:

(П8)

где р – порядок точности численного метода. Для метода Рунге-Кутта 4-го порядка р = 4.

Векторная форма алгоритма метода Рунге-Кутта для задачи задачи Коши имеет вид:

,

,

,

, (П9)

,

,

где

Для данной задачи введенные обозначения конкретизируются следующим образом:

, (П10)

. (П11)

Аргумент задачи – время t. Примем в качестве шага по времени величину h = 0,1.

Будем решать задачу в табличном процессоре MS Excel. Оформим рабочий лист Расчет для ввода исходных данных и оформим строку заголовков для таблицы результатов, как показано на рис. 1. Вычисления организуем в макросе Visual Basic for Applications (VBA), текст которого с комментариями приведен ниже.

' Инструкция об обязательном объявлении переменных

Option Explicit

' Объявление именованной константы

Const pi = 3.1415926

' Объявление объектовых переменных

' ID – массив исходных данных

' Rez – массив результатов

Dim ID As Object

Dim Rez As Object

' Объявление переменных

' Vx, Vy, x, y – проекции скорости, траекторные координаты

' Vn, alfa, r, Ro – начальная скорость и угол старта, радиус торпеды, плотность воды

' mu, c – динамическая вязкость, коэффициент лобового сопротивления

' m, F – масса торпеды, тяга двигателя

' g, L – ускорение свободного падения, глубина залегания торпеды

' Vk, h – скорость корабля, шаг интегрирования

' k1, k2 – линейный и квадратный коэффициент сопротивления

' t, Lk – время, расстояние смещения корабля

Dim Vx As Double, Vy As Double, x As Double, y As Double

Dim Vn As Double, alfa As Double, r As Double, Ro As Double

Dim mu As Double, c As Double, m As Double, F As Double

Dim g As Double, L As Double, Vk As Double, h As Double

Dim k1 As Double, k2 As Double, t As Double

Dim Lk As Double

' Функции для вычисления компонент вектора правых частей (П11)

Function Vxp(Vx As Double, Vy As Double) As Double

Vxp = (F / Sqr(Vx ^ 2 + Vy ^ 2) - k1 - k2 * _

Sqr(Vx ^ 2 + Vy ^ 2)) * Vx / m

End Function

Function Vyp(Vx As Double, Vy As Double) As Double

Vyp = (F / Sqr(Vx ^ 2 + Vy ^ 2) - k1 - _

k2 * Sqr(Vx ^ 2 + Vy ^ 2)) * Vy / m - g

End Function

Function xp(Vx As Double) As Double

xp = Vx

End Function

Function yp(Vy As Double) As Double

yp = Vy

End Function

' Процедура для реализации метода Рунге-Кутта

Sub Расчет()

' Объявление переменных

' i – счетчик шагов цикла

' kvx1, kvy1, kx1, ky1 – компоненты вектора k₁ в (П9)

' kvx2, kvy2, kx2, ky2 – компоненты вектора k₂ в (П9)

' kvx3, kvy3, kx3, ky3 – компоненты вектора k₃ в (П9)

' kvx4, kvy4, kx4, ky4 – компоненты вектора k₄ в (П9)

Dim i As Integer

Dim kvx1 As Double, kvx2 As Double, kvx3 As Double

Dim kvx4 As Double, kvy1 As Double, kvy2 As Double

Dim kvy3 As Double, kvy4 As Double, kx1 As Double

Dim kx2 As Double, kx3 As Double, kx4 As Double

Dim ky1 As Double, ky2 As Double, ky3 As Double

Dim ky4 As Double

' Связывание объектовой переменной ID с диапазоном исходных данных

' листа Расчет

Set ID = Worksheets("Расчет ").Range("A2:A13")

' Связывание объектовой переменной Rez с диапазоном результатов расчета

' листа Расчет

Set Rez = Worksheets("Расчет").Range("A17:F10000")

' Обнуление содержания объектовой переменной Rez

Rez.ClearContents

' Присвоение значений исходным данным из объектовой переменной ID

Vn = ID(1)

alfa = ID(2)

alfa = alfa * pi / 180

r = ID(3)

Ro = ID(4)

mu = ID(5)

c = ID(6)

m = ID(7)

F = ID(8)

g = ID(9)

L = ID(10)

Vk = ID(11)

h = ID(12)

' Вычисление линейного и квадратного коэффициентов сопротивления (П1)

k1 = 6 * pi * mu * r

k2 = c * pi * r ^ 2 * Ro / 2

' Вычисление начальных условий (П2)

Vx = Vn * Cos(alfa)

Vy = Vn * Sin(alfa)

x = 0

y = 0

t = 0

Lk = 0

' Вывод результатов в начальный момент времени

i = 1

Rez(1, 1) = t

Rez(1, 2) = x

Rez(1, 3) = y

Rez(1, 4) = Vx

Rez(1, 5) = Vy

Rez(1, 6) = Lk

' Начало цикло вычислений, который заканчивается, когда координата y превысит

' глубину залегания

Do Until y > L

' Вычисление компонент векторов k1 – k4 по (П9) и (П11)

kvx1 = Vxp(Vx, Vy)

kvy1 = Vyp(Vx, Vy)

kx1 = xp(Vx)

ky1 = yp(Vy)

kvx2 = Vxp(Vx + kvx1 * h / 2, Vy + kvy1 * h / 2)

kvy2 = Vyp(Vx + kvx1 * h / 2, Vy + kvy1 * h / 2)

kx2 = xp(Vx + kvx1 * h / 2)

ky2 = yp(Vy + kvy1 * h / 2)

kvx3 = Vxp(Vx + kvx2 * h / 2, Vy + kvy2 * h / 2)

kvy3 = Vyp(Vx + kvx2 * h / 2, Vy + kvy2 * h / 2)

kx3 = xp(Vx + kvx2 * h / 2)

ky3 = yp(Vy + kvy2 * h / 2)

kvx4 = Vxp(Vx + kvx3 * h, Vy + kvy3 * h)

kvy4 = Vyp(Vx + kvx3 * h, Vy + kvy3 * h)

kx4 = xp(Vx + kvx3 * h)

ky4 = yp(Vy + kvy3 * h)

' Вычисление значений переменных на новом слое по (П9)

Vx = Vx + h * (kvx1 + 2 * kvx2 + 2 * kvx3 + kvx4) / 6

Vy = Vy + h * (kvy1 + 2 * kvy2 + 2 * kvy3 + kvy4) / 6

x = x + h * (kx1 + 2 * kx2 + 2 * kx3 + kx4) / 6

y = y + h * (ky1 + 2 * ky2 + 2 * ky3 + ky4) / 6

' Вычисление новых значений времени, смещения корабля и номера строки

' вывода результатов

t = t + h

Lk = t * Vk

i = i + 1

' Вывод результатов в текущий момент времени

Rez(i, 1) = t

Rez(i, 2) = x

Rez(i, 3) = y

Rez(i, 4) = Vx

Rez(i, 5) = Vy

Rez(i, 6) = Lk

' Оператор конца цикла

Loop

End Sub

Для запуска макроса разместим на рабочем листе элемент управления Кнопка и сформируем процедуру обработки события Click следующего содержания:

Private Sub CommandButton2_Click()

Расчет

End Sub

Результаты расчета показывают, что для приведенных исходных данных время движения торпеды до поверхности воды составило t = 5,4 с, за которое корабль сместится на L_к = 54 м. При этом координата торпеды в момент выхода на поверхность составляет х = 94,691 м, то есть цель не поражена.

Используя разработанную программу, проведем серию расчетов по моделированию изучаемого процесса.

1. Построим траекторию движения торпеды для разных углах старта. Данные расчетов представлены на рис. 2. Из приведенных данных видно, что начиная с некоторого угла старта, торпеда не достигает поверхности воды.

2. Проведем серию расчетов для углов старта  от 90 до 80 через 1 и трех значений тяги двигателя F = ( 30 кН, 50 кН, 100 кН ) для глубины залегания L = 500 м. В расчетах определялись смещение x_вых торпеды по оси Х на момент выхода на поверхность воды и время движения t_вых до поверхности. Результаты расчетов приведены в таблице.

, 	F = 30 кН		F = 50 кН		F = 100 кН
	xвых	tвых	xвых	tвых	xвых	tвых
90	0	11,7	0	8,9	0	6,2
89	57,6	11,8	24,4	8,9	14,0	6,2
88	123,1	12,2	49,3	9,0	27,9	6,2
87	219,8	13,3	73,6	9,0	42,8	6,3
86	-	-	99,7	9,1	56,9	6,3
85	-	-	126,3	9,2	71,1	6,3
84	-	-	156,3	9,4	85,1	6,3
83	-	-	187,7	9,6	101,3	6,4
82	-	-	224,1	9,9	115,4	6,4
81	-	-	262,0	10,2	129,4	6,4
80	-	-	316,5	10,8	146,5	6,5

Из приведенной таблицы видно, что при малой энерговооруженности двигателя и малых углах старта торпеда не достигает поверхности воды. По данным проведенных расчетов построена номограмма определения угла старта  торпеды при заданных тяге F двигателя при заданной глубине залегания L = 500 м для трех значений скорости корабля v = ( 5 м/с, 10 м/с, 15 м/с ). Величина  определяется как абсцисса пересечения графиков функций x_вых ( ) и L_к( )(рис 3). Результаты обработки номограммы приведены в таблице.

	v = 5 м/с	v = 10 м/с	v = 15 м/с
F = 30 кН	89,1	88,0	87,2
F = 50 кН	88,2	86,2	84,3
F = 100 кН	87,9	85,8	83,7

3. Оценим точность полученных результатов используя правило Рунге. Для этого воспользуемся результатами расчетов предыдущего пункта при F = 100 кН и  = 80. Рассмотрим результаты расчета параметров задачи на момент времени t = 6 c, полученные при расчете с шагом h = 0,1 с и h = 0,05 с и оценим абсолютную точность полученных результатов по формуле Рунге, учитывая 4-ый порядок точности используемого метода:

Относительную погрешность определим по формуле:

.

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

	x, м	y, м	vx , м/с	vy , м/с
h = 0,1	130,725	466,604	30,685	77,034
h = 0,1	134,764	465,432	31,602	76,667
d(h/2)	0,2613	0,0761	0,0611	0,0232
 , %	0,199	0,018	0,195	0,032

Таким образом, точность выполненных расчетов не менее 0,2 %.