Моделирование физических процессов

Физика — наука, в которой математическое моделирование является чрезвычайно важным методом исследования. Наряду с традиционным делением физики на экспериментальную и теоретическую сегодня уверенно выделяется третий фундаментальный раздел — вычислительная физика (computational physics). Причину этого в целом можно сформулировать так: при максимальном проникновении в физику математических методов, порой доходящем до фактического сращивания этих наук, реальные возможности решения возникающих математических задач традиционными методами очень ограниченны. Из многих конкретных причин выделим две наиболее часто встречающиеся: нелинейность многих физических процессов (примеры — ниже в тексте) и необходимость исследования совместного движения многих тел, для которого приходится решать системы большого числа уравнений. Часто численное моделирование в физике называют вычислительным экспериментом, поскольку оно имеет много общего с лабораторным экспериментом.

Аналогии между лабораторным и вычислительным экспериментами

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

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

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

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

а сам второй закон динамики приобретает вид

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

Произведем дискретизацию по времени простейшим возможным способом. Если в некоторый момент времени t_o величина s имеет значение s_o, а величина — значение , то в некоторый последующий момент времени будем иметь

Здесь индекс «0» означает величины в начальный момент времени.

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

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

а) меньше вычислений требуется для того, чтобы пройти весь заданный временной интервал;

б) меньшая точность в передаче значений непрерывных функций и их дискретными представлениями — наборами чисел и .

Вопрос о точности результатов является в описываемом моделировании одним из центральных. Он распадается на два: как оценить эту точность и можно ли, уменьшая , достигать все большей точности?

Остановимся вначале на первом. Формулы (4), (5) представляют собой применение метода Эйлера для приближенного решения системы дифференциальных уравнений (3). Наиболее приемлемой при использовании этого и родственного ему методов (например, Рунге — Кутты) является эмпирическая оценка точности. Для этого отрезок проходится с некоторым шагом , а затем с существенно меньшим (например, в два раза) шагом. Сравнение результатов в точках t1, t2, ..., T позволяет составить представление о реальной точности результатов. Если она недостаточна, то следует повторить процесс с еще меньшим шагом.

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

Более эффективными при моделировании процессов, описываемых дифференциальными уравнениями являются методы Рунге — Кутты более высокого порядка аппроксимации, чем метод Эйлера, неявные методы, методы типа "предиктор — корректор", отличающиеся повышенной устойчивостью, и другие, описанные в специальной литературе.

КОНКРЕТНЫЕ ФИЗИЧЕСКИЕ МОДЕЛИ

Свободное падение тела с учетом сопротивления среды. При реальных физических движениях тел в газовой или жидкостной среде трение накладывает огромный отпечаток на характер движения. Каждый понимает, что предмет, сброшенный с большой высоты (например, парашютист, прыгнувший с самолета), вовсе не движется равноускоренно, так как по мере набора скорости возрастает сила сопротивления среды. Даже эту, относительно несложную, задачу нельзя решить средствами “школьной” физики: таких задач, представляющих практический интерес, очень много. Прежде чем приступать к обсуждению соответствующих моделей, вспомним, что известно о силе сопротивления.

Закономерности, обсуждаемые ниже, носят эмпирический характер и отнюдь не имеют столь строгой и четкой формулировки, как второй закон динамики. О силе сопротивления среды движущемуся телу известно, что она, вообще говоря, растет с ростом скорости (хотя это утверждение не является абсолютным). При относительно малых скоростях величина силы сопротивления пропорциональна скорости и имеет место соотношение, где определяется свойствами среды и формой тела. Например, для шарика — это формула Стокса, где — динамическая вязкость среды, r — радиус шарика. Так, для воздуха при t = 20°С и давлении 1 атм = 0,0182 H^.c^.м^-2 для воды 1,002 H^.c^.м^-2 , для глицерина 1480 H^.c^.м^-2.

Оценим, при какой скорости для падающего вертикально шара сила сопротивления сравняется с силой тяжести (в движение станет равномерным).

Пусть r= 0,1 м, = 0,8 кг/м (дерево). При падении в воздухе м/с, в воде 17 м/с, в глицерине 0,012 м/с.

На самом деле первые два результата совершенно не соответствуют действительности. Дело в том, что уже при гораздо меньших скоростях сила сопротивления становится пропорциональной квадрату скорости: . Разумеется, линейная по скорости часть силы сопротивления формально также сохранится, но если , то вкладом можно пренебречь (это конкретный пример ранжирования факторов). О величине k₂ известно следующее: она пропорциональна площади сечения тела S, поперечного по отношению к потоку, и плотности среды и зависит от формы тела. Обычно представляют k₂ = 0,5сS, где с — коэффициент лобового сопротивления — безразмерен. Некоторые значения с (для не очень больших скоростей) приведены на рис.1.

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

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

		0 Диск	с = 1,11
		Полусфера Полусфера	с = 1,33 с = 0,55
		Шар	с = 0,4
		Каплевидное тело Каплевидное тело	с = 0,045 с = 0,01

Рис1. Значения коэффициента лобового сопротивления для некоторых тел, поперечное

сечение которых имеет указанную на рисунке форму

Примем r = 0,1 м, =0,8^.10³ кг/м³ (дерево). Тогда для движения в воздухе (= 1,29 кг/м³ ) получаем 18 м/с, в воде(= 1^.10³ кг/м³ ) 0,65 м/с, в глицерине (= 1,26^.10³ кг/м³ ) 0,58 м/с.

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

Движение является одномерным; проецируя векторное уравнение на ось, направленную вертикально вниз, получаем

Вопрос, который мы будем обсуждать на первом этапе, таков: каков характер изменения скорости со временем, если все параметры, входящие в уравнение (7) заданы? При такой постановке модель носит сугубо дескриптивный характер. Из соображений здравого смысла ясно, что при наличии сопротивления, растущего со скоростью, в какой-то момент сила сопротивления сравняется с силой тяжести, после чего скорость больше возрастать не будет. Начиная с этого момента, , и соответствующую установившуюся скорость можно найти из условия =0, решая не дифференциальное, а квадратное уравнение. Имеем

(второй — отрицательный — корень, естественно, отбрасываем). Итак, характер движения качественно таков: скорость при падении возрастает от до . Как и по какому закону – это можно узнать, лишь решив дифференциальное уравнение (7).

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

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

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

Таблица 2

Зависимость перемещения и скорости падения “безпарашютиста” от времени (от 0 до 15 с)

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

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

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

Приведем конкретный пример решения задачи о свободно падающем теле. Герой знаменитого фильма “Небесный тихоход” майор Булочкин, упав с высоты 6000 м в реку без парашюта, не только остался жив, но даже смог снова летать. Попробуем понять, возможно, ли такое на самом деле или же подобное случается только в кино. Учитывая сказанное выше о математическом характере задачи, выберем путь численного моделирования. Итак, математическая модель выражается системой дифференциальных уравнений.

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

Так как постановка задачи должна быть конкретной, мы примем соглашение, каким образом падает человек. Он опытный летчик и наверняка совершал раньше прыжки с парашютом, поэтому, стремясь уменьшить скорость, он падает не “солдатиком”, а лицом вниз, “лежа”, раскинув руки в стороны. Рост человека возьмем средний — 1,7 м, а полуобхват грудной клетки выберем в качестве характерного расстояния — это приблизительно 0,4 м. для оценки порядка величины линейной составляющей силы сопротивления воспользуемся формулой Стокса. Для оценки квадратичной составляющей силы сопротивления мы должны определиться со значениями коэффициента лобового сопротивления и площадью тела. Выберем в качестве коэффициента число с=1,2 как среднее между коэффициентами для диска и для полусферы (выбор дня качественной оценки правдоподобен). Оценим площадь: S = 1,7 ∙ 0,4 = 0,7(м²).

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

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

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

Вычисления производятся до тех пор, пока “безпарашютист” не опустится на воду. Примерно через 15 с после начала полета скорость становится постоянной и остается такой до приземления (рис7). Отметим, что в рассматриваемой ситуации сопротивление воздуха радикально меняет характер движения. При отказе от его учета график скорости, изображенный на рисунке 1, заменился бы касательной к нему в начале координат.

Рис. 1. График зависимости скорости падения “безпарашютиста” от времени

А теперь ответим на вопрос, поставленный в задаче, Известен такой факт: один из американских каскадеров совершил прыжок в воду с высоты 75 м (Бруклинский мост), и скорость приземления была 33 м/с. Сравнение этой величины с получившейся у нас конечной скоростью 37,76 м/с позволяет считать описанный в кинофильме эпизод вполне возможным. Обсуждаемой модели можно придать черты оптимизационной, поставив задачу так: парашютист прыгает с некоторой высоты и летит, не открывая парашюта; на какой высоте (или через какое время) ему следует открыть парашют, чтобы иметь к моменту приземления безопасную скорость? Или по-другому: как связана высота прыжка с площадью поперечного сечения парашюта (входящей в k₂ чтобы скорость приземления была безопасной? Выполнение таких исследований многократно более трудоемко, нежели просто изучение одного прыжка при заказанных условиях.

Взлет ракеты

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

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

Однако при самом простом моделировании данного процесса необходимо принять во внимание, что плотность воздуха , входящая в коэффициент k₂, убывает по мере подъема ракеты по закону , где h — высота, — иначе модель будет совершенно неадекватна реальности. Таким образом, модель будет описываться системой двух дифференциальных уравнений для функций v(t) и h(t):

• — начальная масса ракеты, заправленной топливом;

• — остаточная масса после полного выгорания топлива;

• величины, определяющие , — коэффициент сопротивления воздуха (линейной составляющей силы сопротивления можно заведомо пренебречь);

• — сила тяги двигателя (принять постоянной).

Движение тела, брошенного под углом к горизонту

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

• a — угол начального наклона вектора скорости к горизонту;

• величины, определяющие коэффициенты сопротивления среды k₁ и k₂.