Трушков А.С. Учебно-методический комплекс по дисциплине Компьютерное моделирование

Подождите немного. Документ загружается.

Вариант 18. Промоделировать траекторию движения малого космического

аппарата, запускаемого с борта космической станции, относительно Земли. Запуск

осуществляется путем толчка в направлении, противоположном движению стан-

ции, по касательной к ее орбите.

Вариант 19. Промоделировать траекторию движения малого космического

аппарата, запускаемого с борта космической станции, относительно Земли. Запуск

осуществляется путем толчка в направлении, совпадающим с направлением дви-

жения станции, по касательной к ее орбите.

Вариант 20. Как будет выглядеть полет искусственного спутника Земли,

если учесть возмущающее действие Луны?

Вариант 21. Разработать и реализовать модель движения искусственного

спутника Земли при учете воздействия на него малой постоянной силы, обуслов-

ленной «солнечным ветром».

Вариант 22. Считая, что движение Луны вокруг Земли происходит практи-

чески по круговой орбите, проанализировать воздействие на эту орбиту со сторо-

ны Солнца для малого участка движения.

Вариант 23. Проанализировать особенности движения искусственного

спутника Земли, движущегося практически по круговой орбите на высоте порядка

300 км, связанные с малым сопротивлением атмосферы.

Вариант 24. Проанализировать изменение круговой орбиты астероида,

движущегося вокруг Солнца, под влиянием вулканического выброса с его по-

верхности.

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

Задание. Рассмотреть движение заряженной частицы в электростатическом

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

частицы при разных соотношениях величин и знаков неподвижных зарядов, а

также начальных траекторных условий. Рассмотреть движение частицы в плоско-

сти неподвижных зарядов.

αα

, y

)

, y

)

, y

)

, y

)

, y

)

Расчетная схема

Будем считать, что частица имеет массу m и заряд q, неподвижные заряды

имеют величины Q

, Q

. Координаты расположения 3-х зарядов равны (x

), (x

, y

), (x

, y

) соответственно. Начальные координаты движущейся частицы

, y

), начальная скорость определяется ее модулем v

и углом наклона

αα

к по-

ложительному направлению оси Ох. Угол

αα

считается положительным, если он

отсчитывается от положительного направления оси Ох против часовой стрелки.

Текущие координаты частицы равны (x, y). Вектора расстояния между движущей-

ся и неподвижными частицами равны

(

((

( )

))

)

111

yy,xxr −

−−

−=

→

→→

→

(

((

( )

))

)

222

yy,xxr −

−−

−=

→

→→

→

(

((

( )

))

)

333

yy,xxr −

−−

−=

→

→→

→

. Расчетная схема приведена на рисунке.

Электростатические силы равны:

3,2,1i,r

=−

−−

−=

→

→→

→

πε

πεπε

πε

Уравнения движения имеют вид:















































→

→→

→

→→

→

→→

→

→→

→

∑

∑∑

∑

Здесь

(

((

( )

))

)

y,xr =

→

→→

→

- радиус-вектор движущейся частицы.

После проектирования на оси координат получим:































































−

−−

−+

+−

−−

−

−−

−

−−

−=

−

−−

−+

+−

−−

−

−−

−

−−

−=

∑

∑∑

∑

∑∑

∑

))yy()xx((

)yy(Q

))yy()xx((

)xx(Q

πε

πεπε

πε

πεπε

πε

В качестве масштабных величин выберем характеристики круговой орбиты

радиуса

ρρ

для заряда q массой m в электростатическом поле неподвижного заряда

Q; скорость равномерного вращения по этой орбите:

ρρ

ρπε

πεπε

πε

период обращения по круговой орбите:

42T

ρρ

πε

πεπε

πεπρ

πρπρ

πρ

Вводим безразмерные переменные:

ττ

ρρ

После подстановки безразмерных переменных в систему и преобразова-

ний получим:







































































−

−−

−+

+−

−−

−

−−

−

−−

−=

−

−−

−+

+−

−−

−

−−

−

−−

−=

∑

∑∑

∑

∑∑

∑

iix

))YY()XX((

)YY(a

))YY()XX((

)XX(a

ππ

ττ

ππ

ττ

ππ

ττ

ππ

ττ

где

ππ

ρρ

πε

πεπε

πε

ππ

ρρ

TQq

−

−−

−=

=−

−−

−=

Будем регулировать величину шага интегрирования в зависимости от траек-

торных параметров с помощью следующего алгоритма. Вычисляются расстояния

от движущейся частицы до всех трех неподвижных зарядов:

3,2,1i,)yy()xx(r

=−

−−

−+

+−

−−

−=

Определяется минимальное из этих расстояний r

min

. В качестве шага по

времени выбирается интервал времени движения по заданной доле круговой тра-

ектории радиуса r

min

с текущей скоростью частицы

VVV +

. В данном

примере эта доля равна n-ой части длины окружности радиуса r

min

. Таким обра-

зом, после каждой итерации выбирается новый шаг по времени, равный:

min

VVn

ππ

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

αα

sinV)0(V,cosV)0(V,Y)0(Y,X)0(X

0x0x00

Будем интегрировать систему обыкновенных дифференциальных уравне-

ний классическим методом Рунге-Кутта 4-го порядка. Вид рабочего лист Рас-

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

на рис. 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, t – проекции скорости, траекторные координаты, время

' x1, y1, x2, y2, x3, y3 – координаты неподвижных зарядов

' a1, a2, a3 – величины неподвижных зарядов

' x0, y0 – начальные траекторные координаты

' V0, alfa – начальная скорость и угол наклона траектории

' Tmax – предельное время счета

' h – шаг по времени

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

Dim t As Double

Dim x1 As Double, y1 As Double, x2 As Double, y2 As Double

Dim x3 As Double, y3 As Double

Dim a1 As Double, a2 As Double, a3 As Double

Dim x0 As Double, y0 As Double

Dim V0 As Double, alfa As Double

Dim Tmax As Double, h As Double

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

Function Vxp(x As Double, y As Double) As Double

Vxp = -2 * pi * (a1 * (x - x1) / Sqr(((x - x1) ^ 2 + _

(y - y1) ^ 2) ^ 3) + a2 * (x - x2) / Sqr(((x - x2) ^ 2

+ (y - y2) ^ 2) ^ 3) + a3 * (x - x3) / Sqr(((x - x3) ^

2 _

+ (y - y3) ^ 2) ^ 3))

End Function

Function Vyp(x As Double, y As Double) As Double

Vyp = -2 * pi * (a1 * (y - y1) / Sqr(((x - x1) ^ 2 + _

(y - y1) ^ 2) ^ 3) + a2 * (y - y2) / Sqr(((x - x2) ^ 2

+ (y - y2) ^ 2) ^ 3) + a3 * (y - y3) / Sqr(((x - x3) ^

2 _

+ (y - y3) ^ 2) ^ 3))

End Function

Function xp(Vx As Double) As Double

xp = 2 * pi * Vx

End Function

Function yp(Vy As Double) As Double

yp = 2 * pi * Vy

End Function

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

Sub Расчет()

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

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

' kvx1, kvy1, kx1, ky1 – компоненты вектора k

в (П9) в лаб/раб № 1

' kvx2, kvy2, kx2, ky2 – компоненты вектора k

в (П9) в лаб/раб № 1

' kvx3, kvy3, kx3, ky3 – компоненты вектора k

в (П9) в лаб/раб № 1

' kvx4, kvy4, kx4, ky4 – компоненты вектора k

в (П9) в лаб/раб № 1

' r1, r2, r3 – расстояния от движущейся частицы до неподвижных зарядов

' r – минимальное из этих расстояний

' tt – имя диаграммы

' n – частота дробления траектории

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, ky4 As

Double

Dim r1 As Double, r2 As Double, r3 As Double, r As Double

Dim tt, n As Integer

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

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

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

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

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

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

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

Rez.ClearContents

' Запоминание имени диаграммы и переопределение диапазона исходных данных

tt = ActiveSheet.ChartObjects(1).Name

ActiveSheet.ChartObjects(tt).Activate

ActiveChart.ChartArea.Select

ActiveChart.SetSourceData Source:=Sheets("Расчет"). _

Range("B22:C20000"), PlotBy:=xlColumns

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

x1 = ID(1)

y1 = ID(2)

x2 = ID(3)

y2 = ID(4)

x3 = ID(5)

y3 = ID(6)

a1 = ID(7)

a2 = ID(8)

a3 = ID(9)

x0 = ID(10)

y0 = ID(11)

V0 = ID(12)

alfa = ID(13)

Tmax = ID(14)

n = ID(17)

' Перевод угловой меры в радианную

alfa = alfa * pi / 180

' Вычисление начальных условий 17

Vx = V0 * Cos(alfa)

Vy = V0 * Sin(alfa)

h = 0.001

x = x0

y = y0

t = 0

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

i = 1

Rez(1, 1) = t

Rez(1, 2) = x

Rez(1, 3) = y

Rez(1, 4) = Vx

Rez(1, 5) = Vy

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

' предельное значение Tmax

Do Until t > Tmax

' Вычисление компонент векторов k1 – k4 по (П9) из лаб/раб № 1

kvx1 = Vxp(x, y)

kvy1 = Vyp(x, y)

kx1 = xp(Vx)

ky1 = yp(Vy)

kvx2 = Vxp(x + kx1 * h / 2, y + ky1 * h / 2)

kvy2 = Vyp(x + kx1 * h / 2, y + ky1 * h / 2)

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

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

kvx3 = Vxp(x + kx2 * h / 2, y + ky2 * h / 2)

kvy3 = Vyp(x + kx2 * h / 2, y + ky2 * h / 2)

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

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

kvx4 = Vxp(x + kx3 * h, y + ky3 * h)

kvy4 = Vyp(x + kx3 * h, y + ky3 * h)

kx4 = xp(Vx + kvx3 * h)

ky4 = yp(Vy + kvy3 * h)

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

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

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

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

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

' Вычисление расстояния от движущейся частицы до неподвижных зарядов

(15)

r1 = Sqr((x - x1) ^ 2 + (y - y1) ^ 2)

r2 = Sqr((x - x2) ^ 2 + (y - y2) ^ 2)

r3 = Sqr((x - x3) ^ 2 + (y - y3) ^ 2)

' Вычисление минимального расстояния до зарядов

If r1 < r2 Then r = r1 Else r = r2

If r3 < r Then r = r3

' Вычисление шага по времени (16)

h = 2 * pi * r / 500 / Sqr(Vx ^ 2 + Vy ^ 2)

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

t = t + h

i = i + 1

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

Rez(i, 1) = t

Rez(i, 2) = x

Rez(i, 3) = y

Rez(i, 4) = Vx

Rez(i, 5) = Vy

' Программный код для графического вывода в процессе счета

' (снять комментарии с 4-х строк)

' ActiveSheet.ChartObjects(tt).Activate

' ActiveChart.ChartArea.Select

' ActiveChart.SetSourceData Source:=Sheets("Расчет").

' Range("B22:C20000"), PlotBy:=xlColumns

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

Loop

End Sub

Для запуска макроса разместим на рабочем листе элемент управления

Кнопка с надписью Пуск и сформируем процедуру обработки события Click

следующего содержания:

Private Sub CommandButton1_Click()

Расчет

End Sub

Для анимации движения частица используется следующий макрос:

Private Sub CommandButton2_Click()

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

' t_name – имя диаграммы

' t_range – текстовая переменная для диапазона исходных данных диаграммы

' ii – текущий номер строки в массиве результатов

' ii0 – номер максимальной строки в предыдущем кадре

' dtv – шаг вывода фрагментов траектории

' dtimer – шаг задержки кадров, с

Dim t_name, t_range, Start

Dim ii As Integer, ii0 As Integer

Dim dtv As Double, dtimer As Double

' связывание объектовых переменных с диапазонами ячеек рабочего листа

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

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

' ввод необходимых исходных данных

Tmax = ID(14)

dtv = ID(15)

dtimer = ID(16)

' Запоминание имени диаграммы

t_name = ActiveSheet.ChartObjects(1).Name

' Определение начальных значений номеров строк

ii = 1

ii0 = 1

' Цикл вывода траектории до момента превышения предельного времени

Do Until Rez(ii, 1) > Tmax - dtv

' Цикл определения максимальной строки в текущем кадре

Do Until Rez(ii, 1) > Rez(ii0, 1) + dtv

ii = ii + 1

Loop

' Запоминание максимальной строки в текущем кадре

ii0 = ii

' Формирование текстовой переменной с диапазоном исходных данных

' диаграммы для текущего кадра

t_range = "B22:C" & CStr(21 + ii)

' Построение текущего кадра диаграммы

ActiveSheet.ChartObjects(t_name).Activate

ActiveChart.ChartArea.Select

ActiveChart.SetSourceData Source:=Sheets("Ðàñ÷åò").

Range(t_range), PlotBy:=xlColumns

' Запоминание текущего времени

Start = Timer

' Цикл задержки текущего кадра

Do While Timer < Start + dtimer

Loop

End Sub

Вид диаграммы для приведенных выше исходных данных приведен на рис.

Используя разработанную программу можно проводить разнообразные ис-

следования.

1) Рассмотрим, как изменятся результаты расчета для тех же исходных дан-

ных, если изменить шаг интегрирования. На рис. 3 приведены результаты расчета

при частоте дробления траектории n = 100 и заметно существенное отличие тра-

екторных координат. На рис. 4 приведены результаты расчета при частоте дроб-

ления траектории n = 1000 и результаты практически не отличаются по сравнению

со случаем n = 500. Таким образом, для данного расчетного случая частота дроб-

ления траектории n = 500 обеспечивает удовлетворительную точность результа-

тов.

2) Рассмотрим, как изменятся результаты для тех же исходных данных, если

изменится только начальная скорость. На рис. 5 представлена траектория движе-

ния частицы при начальной скорости V

= 1,9 и максимальном времени счета T

max

= 1000. Анализ траекторий показывает, что при увеличении времени процесса в

100 раз по сравнению с опорным случаем (рис. 2) число оборотов вокруг системы

зарядов увеличилось с 3 до 5, но с существенно большими размерами и меньшими

скоростями. При увеличении начальной скорости до

= 2,0 заряженная частица не возвращается к данной системе зарядов. На рис. 6

приведены диаграммы скоростей для этих 3-х траекторий, которые подтверждают

сделанные выводы.

Рис. 1.