Методичка по МС

Компьютерное моделирование физических процессов

4 Интегрирование простейших дифференциальных уравнений

4.1 Движение тела в поле тяжести Земли без учета сил вязкого трения

Рассмотрим тело, брошенное под углом α к поверхности Земли, с начальной скоростью v0 . На тело

действует одна сила — сила тяжести (рисунок 4.1). Под действием этой силы тело движется с ускорением свободного падения g . Начальными условиями являются

x0 = 0 , y0 = 0, v0 x = v0 cosα , v0 y = v0 sinα .

Согласно основному закону динамики:

ma = mg .

В проекции на оси координат

Рисунок 4.1

(4.1)

m

d 2 x

= 0 , m

d 2 y

= −mg .

(4.2)

dt2

dt

2

После сокращения массы и понижения порядка дифференциального уравнения

d 2 x

=

dv

x

= 0 ,

d 2 y

=

dvy

= −g .

(4.3)

dt2

dt

dt2

dt

Согласно методу Эйлера, выберем малое значение ∆t

в качестве кванта вре-

мени, тогда, вместо дифференциалов скоростей запишем малые конечные приращения ∆vx , ∆vy . По прошествии первого интервала времени ∆t , проекции скоростей

изменятся на величины ∆v0 x , ∆v0 y и будут равны

v1x = v0 x +∆v0 x , v1y = v0 y +∆v0 y .	(4.4)
Приращения скоростей находим из (4.3)
∆vx = 0 ∆t , ∆vy = −g ∆t .	(4.5)

Эти приращения в данном случае не зависят от скорости, поэтому числовые индексы можно не указывать. Из (4.5) видно, что приращение горизонтальной скорости равно нулю, а вертикальная скорость изменяется равномерно (равноускоренное движение). Скорость в любой момент времени может быть вычислена при помощи рекурсивной функции

	vi+1x = vix , vi+1y = viy − g∆t .		(4.6)
Для отыскания координат тела необходимо
на каждом шаге определения скорости вычислять
и координату, считая, вследствие малости отрезка
времени	∆t ,	движение	равномерным.
Аналогичным способом, как и для скорости, полу-
чаем, что	xi+1 = xi +∆xi = xi +vi ∆t , yi+1 = yi +∆yi					Рисунок 4.2
				= yi +vi ∆t .		(4.7)

11

Компьютерное моделирование физических процессов

Таким образом, формулы (4.6) и (4.7) — рекурсивные выражения для скорости точки и ее координат.

В качестве алгоритма решения этой задачи можно использовать блок-схему — рисунок 4.2.

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

v		= v	cosα , x = x	+v cosαt , v		= v sinα − gt ,	y = y	+v sinαt −	gt	2	(4.8)
	x				y					.
		0	0	0		0	0	0	2

4.2 Движение тела в поле тяжести Земли с учетом сил вязкого трения

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

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

сопротивления, направленная в сторону противоположную

движению (сила сопротивления противоположна вектору

скорости) — рисунок 4.3.

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

пишется в виде

Рисунок 4.3

ma = mg + R .

(4.9)

Будем считать, что сила сопротивления пропорциональна скорости движения

(4.10)

R = −kv .

В проекции на оси координат уравнение (4.9) с учетом (4.10) примет вид

m

d 2 x

= m

dv

x

= −kvx , m

d 2 y

= m

dvy

= −mg −kvy

(4.11)

dt2

dt2

dt

dt

Вновь, вместо бесконечно малого приращения времени dt будем использовать

малое конечное приращение времени ∆t ,

тогда приращения проекций скорости на

оси координат будут равны

k

k

∆vx = −

vx ∆t , ∆vy = − g +

vy ∆t .

(4.12)

m

m

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

		k		k
	vi+1x = vix +∆vix = vix −		vix ∆t , vi+1 y = viy +∆viy = viy − g +		viy	∆t .	(4.13)
		m		m
	Координаты в этом случае будут равны
	xi+1 = xi +∆xi = xi +vi ∆t , yi+1 = yi +∆yi = yi +vy ∆t .						(4.14)

При расчетах координат вновь предполагаем, что вследствие малости интервала времени ∆t движение в течение каждого временного интервала является равномерным. Такое приближение будет тем более точным, чем меньший промежуток времени ∆t выбран.

12

Компьютерное моделирование физических процессов

4.3 Колебания математического маятника

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

В идеальном случае, тело совершает колебания под действием двух сил — си-

лы тяжести G и силы натяжения нити T					(рисунок 4.4). Запишем
для него основной закон динамики вращательного движения
относительно точки подвеса маятника в проекции на ось z
Iε = −Gl sinϕ					(4.15)
Преобразуем полученный закон					(4.16)
2
ml2 d ϕ + mgl sinϕ			= 0.
dt2
И введем новое обозначение ω2	=	g		и получим
0		l
					(4.17)
ϕ +ω02 sinϕ = 0 .							Рисунок 4.4

Полученное уравнение не имеет решения в виде конечных сумм, и, поэтому, сделаем упрощение sinϕ ≈ϕ , при условии малости угла ϕ .

ϕ +ω2ϕ = 0 .	(4.18)
0

Решением этого дифференциального уравнения является гармоническая функция

ϕ =ϕ0 sin (ω0t + β ),

(4.19)

где ϕ0 и β определяются из начальных условий

ϕ

=

ϕ2

+

ϕ2

, sin β =

ϕ

н , cos β =

ϕ

н

.

(4.20)

н

ϕ

ϕ ω

0

н

ω2

0

0

0

0

Сила натяжения, возникающая в точке подвеса нити, определяется из второго закона Ньютона в проекции на нормальную ось (рисунок 4.4)

	m	v2	= −mg cosϕ + N .	(4.21)
		l

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

d 2ϕ = dω dϕ =ω dω , dt2 dt dϕ dϕ

иподставляем ее в дифференциальное уравнение угла поворота

ωddωϕ = − gl sinϕ .

(4.22)

(4.23)

Решаем полученное дифференциальное уравнение методом разделения переменных

g

ω

g

ϕ

g

(cosϕ −cosϕн ).

ωdω = −

sinϕdϕ ,

∫

ωdω = −

∫sinϕdϕ,

ω2 =ωн2 +2

(4.24)

l

ω

н

l ϕ

н

l

Умножаем правую и левую части уравнения на квадрат длины нити

13

Компьютерное моделирование физических процессов
l2ω2 = l2ωн2 +2gl (cosϕ −cosϕн ), v2		= vн2 +2gl (cosϕ −cosϕн ),	(4.25)
и подставляем в уравнение (4.21)
N = m v2	+mg cosϕ = m (vн2 +2gl (cosϕ −cosϕн ))+mg cosϕ =
l	l		(4.26)
	2
	= mg vн +3cosϕ −2cosϕн
	gl

Уравнение (4.26) получено без учета упрощения sinϕ ≈ϕ и, поэтому, справед-

ливо для любых углов отклонения.

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

	Iϕ = −mgl sinϕ −αϕ , или ϕ +		α	ϕ +	g	sinϕ = 0 ,			(4.27)
			ml2		l
							g
	с учетом введения новых обозначений 2δ =		α	и ω2		=
			ml2	0			l
	ϕ +2δϕ +ω2 sinϕ = 0 .								(4.28)
	0

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

чальное угловое ускорение:

ϕн = −(2δϕн +ω02 sinϕн ).

(4.29)

Решение уравнения (4.28) в виде конечных сумм не существует, поэтому будем решать его численными методами. Для этого понизим степень дифференциального уравнения, переходя к угловой скорости

dϕ +2δϕ +ω02 sinϕ = 0 .	(4.30)
dt
Теперь, используя метод Эйлера, определяем приращение угловой скорости
∆ϕ1 = −(2δϕн +ω02 sinϕн )∆t .	(4.31)

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

	ϕ1 =	∆ϕ1	, ϕ1 =ϕн +∆ϕ1 , ϕ1 =ϕн +ϕ1∆t .	(4.32)
		∆t

Для угла поворота, угловых скорости и ускорения в любой момент времени, обобщим проделанные операции на любое количество итераций

∆ϕi+1 = −(2δϕi +ω02 sinϕi )∆t , ϕi+1 =	∆ϕi+1	, ϕi+1 =ϕi +∆ϕi+1 , ϕi+1 =ϕi +ϕi+1∆t .	(4.33)
	∆t

Построенная модель позволяет решить практически любую задачу колебаний маятника. При условии, что α = 0 (δ = 0 ) получим колебания идеального маятника.

При малых углах отклонения (менее 100 ), колебания соответствуют классической теории математической маятника. При большом значении коэффициента затухания

14

Компьютерное моделирование физических процессов

δ колебания как таковые не наблюдаются, а движение становиться апериодическим.

Решение дифференциальных уравнений методами математического анализа, во-первых, не позволяет получить точное решение уравнений (4.17) и (4.28), вовторых, при переходе к апериодическому движению изменяется способ решения уравнений и в итоге получается три уравнения, описывающих колебания мятника при различных соотношениях между коэффициентом затухания и собственной частотой колебаний системы.

Сила натяжения нити не зависит явно от силы сопротивления и в любом случае определяется выражением (4.26) как функция угла.

Рассмотрим колебания математического маятника в переменном электрическом поле (рисунок 4.5). Переменное электрическое поле создается двумя параллельными пластинами — конденсатором. Напряжение на обкладках конденсатора меняется по гармоническому закону и тогда напряженность поля внутри конденса-

тора будет равна

E = E0 sin(ωt + βE ) .

(4.34)

Если колеблющееся тело математического маятника

несет на себе заряд q , тогда кроме силы тяжести и

натяжения нити на тело будет действовать еще и пе-

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

F = E0qsin (ωt + βE ),

(4.35)

и основной закон динамики вращательного движения для

тела примет вид

Iϕ = −mgl sinϕ −αϕ + E0 qsin(ωt + βE )cosϕ , или

(4.36)

ϕ +

α

ϕ +

g

sinϕ −

E0q sin(ωt + βE )cosϕ = 0

ml2

Рисунок 4.5.

l

ml2

g

с учетом введения новых обозначений 2δ =

α

, ω2

=

и

f = E0 q

ml2

0

l

ml2

ϕ +2δϕ +ω02 sinϕ − f sin (ωt + βE )cosϕ .

(4.37)

Как и уравнения (4.17) и (4.28) уравнение (4.37) не имеет решения в конечном виде, поэтому его решение следует искать численными методами. Решение получается аналогично решению уравнения (4.28).

По начальным условиям ϕн и ϕн вычисляем начальное угловое ускорение
ϕн = −(2δϕн +ω02 sinϕн )+ f sin βE cosϕн .	(4.38)
Затем, вычисляем изменение скорости через время ∆t
∆ϕ1 = (−2δϕн −ω02 sinϕн + f sin βE cosϕн )∆t .	(4.39)

По полученному изменению угловой скорости находим значение угла поворота и углового ускорения через интервал времени ∆t

∆ϕi+1 = (−2δϕi −ω02 sinϕi

+ f sin (ωti + βE )cosϕi )∆t ,

(4.40)

ϕ

=

∆ϕi+1

, ϕ

=ϕ

+∆ϕ

, ϕ

=ϕ

+ϕ

∆t , t

= t +∆t .

∆t

i+1

i+1

i

i+1

i+1

i

i+1

i+1

i

15

Компьютерное моделирование физических процессов

При этом необходимо фиксировать не только интервал времени ∆t , но и само время t , потому что при расчете электрической силы необходимо знать не приращение времени, а его значение.

Рассмотрим колебания математического маятника с учетом собственного движения точки подвеса (маятник Капицы) рисунок 4.6.

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

Дифференциальное уравнение движения колебания точки имеет вид

ϕ +δϕ +

ω

2

+

Bω2

cosωt +δ

Bω

sinωt

sinϕ = 0,

(4.41)

0

l

l

g

где δ — параметр, характеризующий затухание в системе, ω2

=

— Рисунок 4.6.

0

l

квадрат собственной частоты колебаний системы, B,ω— амплитуда и частота вер-

тикальных колебаний точки подвеса.

Записанное уравнение описывает движение маятника при наличии силы трения, пропорциональной скорости тела. Колебания подобных систем называются параметрическими. В некоторых диапазонах частот вынуждения амплитуда колебаний возрастает и наступает параметрический резонанс.

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

	2gl	(4.42)
ω >	B .

16

Компьютерное моделирование физических процессов

5 Модели нескольких тел

5.1 Движение тел Солнечной системы

Рассмотрим движение тел солнечной системы под действием сил гравитационного притяжения.

Основной целью построения такой модели является:

•предсказание природных явлений,

•проверка табличных данных,

•демонстрация отклонений траекторий тел от круговых,

•колебания Солнца,

•доказательство законов Кеплера.

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

Начальными условиями являются местоположение и скорость тел солнечной системы и их спутников. Из спутников будем рассматривать только Луну. В итоге получим взаимное движение 9 планет, Солнца и Луны, всего 11 тел.

Тела движутся по законам Ньютона, на каждое тело будет действовать 10 сил гравитационной природы:

11
mj aj = ∑Pij ,	(5.1)
i=1
где m j ,aj — масса и ускорение j -ого тела, Pij — сила действующая на	j -е тело со

стороны i -ого тела (при составлении уравнения необходимо пропускать одно сла-

гаемое — силу Pij	при i = j ). В проекциях на оси координат для каждого тела будет
получено по два уравнения (считаем, что Солнечная система
является плоской)	11	11	(5.2)
mj xj = ∑Xij , mj yj = ∑Yij .
	i=1	i=1
Эти дифференциальные уравнения для каждого тела с учетом
начальных условий однозначно будут определять дальнейшее
движение тела.		j -ого тела (рисунок 5.1).
Рассмотрим	движение
Систему координат выберем так, чтобы Солнце
располагалось в начальный момент времени в начале коор-
динат. Координаты тел будут равны xi , yi , их скорости xi , yi .
					Рисунок 5.1

Тогда на j -е тело будут действовать силы

Xij =G	mi m	j	cosαij , Yij =G	mi m	j	sinαij	(5.3)
	2			2
	r			r
	ij			ij

17

Компьютерное моделирование физических процессов

где G = 6,67 10−11 Нмкг22 — постоянная всемирного тяготения, mi , mj — массы взаимо-

действующих тел, rij — расстояние между телами, αij — угол между направлением

действующих сил и горизонтальной осью.

Расстояние между телами можно отыскать по теореме Пифагора

rij2 = (xi − xj )2 +(yi − yj )2 ,					(5.4)
а синус и косинус углов также определяется геометрически через координат тел
cosαij =	xi − xj	, sinαij =	yi − yj	.	(5.5)
	rij
			rij

Определение не угла, а сразу соответствующей тригонометрической функции позволяет правильно определить ее знак.

Определив для каждого тела проекции всех сил, следует просуммировать их и определить проекции равнодействующей согласно (5.2). Зная равнодействующую можно определить изменение скорости j -ого тела

∆xj

11

∆t

11

mj

= ∑Xij (i ≠ j) ,

∆xj =

∑Xij (i ≠ j) ,

∆t

m

=

j

=

i 1

i 1

(5.6)

∆yj

11

∆t

11

mj

= ∑Yij (i ≠ j) ,

∆yj =

∑Yij (i ≠ j) .

∆t

m

=

j

=

i 1

i 1

Тогда скорость этого тела и координаты будут равны xj = xj +∆xj , xj = xj + xj ∆t ,

yj = yj +∆yj , yj = yj + yj ∆t .

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

(5.7)

которые были оговорены

5.2 Движение молекул идеального газа

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

стенками сосуда.

При ударе молекул изменяется модуль и направление их скорости. Удар упругий, значит, при взаимодействии будет

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

До взаимодействия частиц они имели скорости v1, v2 , после взаимодействия их скорости стали равны u1, u2 . Для простоты решения будем считать, что массы всех

молекул одинаковы. Законы в векторном виде:
mv +mv = mu +mu ,				mv2	mv2	mu2	mu2	(5.8)
				1 +	2 =	1 +	2 .
1	2	1	2	2	2	2	2
18

Компьютерное моделирование физических процессов

Проводим преобразования законов
							v1 −u1 = u2 −v2 ,											(5.9)
v2	−u2	= v2	−u2	(v −u )(v +u				)= (v	−u	)(v	+u	)	(v +u )= (v			+u	).
1	1	2	2	1	1	1	1	2	2	2	2		1	1	2	2

После проецирования на оси координат x, y, z получаем систему из шести уравне-

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

5.3 Движение молекул реального газа

Рассмотрим модель, более приближенную к реальному газу — молекулы имею различные массы, между молекулами существуют силы притяжения и отталкивания, при этом прямого механического контакта (удара) между молекулами не происходит. Взаимодействие между молекулами определяется наиболее часто по формулам Леннард– Джонса (5.10) или Букингема (5.11) (рисунок 5.3)

U (r) = −ar−6	+br−12	(5.10)
U (r) = −ar−6	+be−cr	(5.11)
коэффициенты в формулах связаны с глубиной и			Рисунок 5.3
положением потенциальной ямы		и определяются из

экспериментальных данных. Для определения силы взаимодействия необходимо воспользоваться соотношением

f (r)= −	∂U (r)	.	(5.12)
	∂r
При движении молекул выполняются законы Ньютона, поэтому можно запи-
сать основное уравнение динамики
N
mj aj = ∑ f (rij ) +mj g .			(5.13)

i=1, i≠ j

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

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

19

Рисунок 5.4

Компьютерное моделирование физических процессов

				N	еслиi ≠ j и r		< ρ, то f (r )
m	a		=	∑i=1		ij	ij	+ m	g .	(5.14)
j		j			еслиi = j или rij ≥ ρ, то 0			j

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

Взаимодействие молекул со стенкой можно рассчитывать по той же формуле, пологая, в простейшем случае, что и коэффициенты потенциальной функции те же. В этом случае следует считать, что сила со стороны стенки действует по нормали к ней, а rj — нормальный вектор к поверхности стенки сосуда, направленный на дви-

жущуюся частицу. Если молекула движется вблизи угла, то воздействий может быть два или три (в зависимости от положения молекулы, рисунок 5.4).

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

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

распределение молекул по высоте позволит получить барометрическую формулу и так далее.

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

5.4 Модель электронного газа

Модель электронного газа абсолютно аналогична модели реального газа, с тем отличием, что кроме взаимодействия с соседними электронами газа происходит взаимодействие еще и с ионами в узлах кристаллической решетки.

20