Лекция 19.
Уравнения в частных производных

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

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

В качестве примера приведем проект для изучения теплопроводности стены дома. На рис. 19.1 представлена реализация такой системы в среде «Stratum-2000». Схема собрана из элементов (имиджей). Элементы двухмерной системы имитируют кирпичи, из которых сложена стена. На рис. 19.1 видны кирпичи, сделанные из разного материала (что обозначено желтым, голубым, красным цветами) и уложенные в определенном порядке, что, по мнению автора схемы, улучшает теплоизоляцию помещения. Серые и зеленые элементы на краях стены задают краевые условия, имитируя окружающую среду. Рядом показано для примера, как элементы в такой системе связаны друг с другом. По связям соседние элементы в процессе расчета обмениваются информацией между собой о температуре. Заметим, что в данном проекте, который имеет вид конструктора, пользователь схемы может менять параметры элементов системы, условия, конструкцию стены, порядок укладки (топологию), что при решении задачи аналитическими методами крайне затруднительно и требует применения методов моделирования.

Рис. 19.1. Двухмерная система с распределенными параметрами для исследования свойств теплопроводности стены дома

Существуют типовые уравнения, описывающие отдельные свойства систем с распределенными параметрами. Рассмотрим их.

Уравнение диффузии

Уравнение диффузии описывает распространение (растекание) со временем по протяженному телу некоторой субстанции, например, тепла или концентрации. В одномерном случае тело представляется протяженным вдоль оси x.

На рис. 19.2 показан пример распределения вдоль оси x такого параметра как температура T. Из обычного опыта хорошо известно, что в каждый момент времени t температура T на разных участках тела x имеет разные значения, то есть меняется в зависимости от участка и времени. То есть должен существовать закон, по которому изменяется величина этого параметра T как функции от (x, t). Для температуры этот закон чаще всего задается уравнением диффузии.

Рис. 19.2. Пример применения уравнения диффузии к описанию процесса нагрева протяженного одномерного тела

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

и обычно дополняется условиями — значениями переменной y на краях и границах: на левом краю x = 0, на правом краю x = L, на границе — начальные условия (t = 0):

y(x, 0) = f₁(x),
y(0, t) = f₂(t),
y(L, t) = f₃(t),
где f₁(x), f₂(t) и f₃(t) — заданные функции.

На рис. 19.3 представлен схематически вид области, для которой определены граничные и начальные условия. Функции f₁(x), f₂(t), f₃(t) и само уравнение диффузии предопределяют поведение функции y(x, t) внутри этой области, чей полный вид обычно надо определить. Если на схеме дополнительно построить ось y (см. рис. 19.4), то визуально на рисунке можно отобразить и сам вид функций. На рисунке четко видно, что в углах схемы значения задаваемых функций должно совпадать.

Рис. 19.3. Схема задания начальных и краевых условий для систем с распределенными параметрами в осях x, t (пример)

Рис. 19.4. Схема задания начальных и краевых условий для систем с распределенными параметрами в осях x, t, y (пример)

Коэффициент α имеет смысл коэффициента теплопроводности; f(x, t) имеет смысл функции, описывающей работу источников и стоков тепла.

Величина y, описывающая распределение температуры, является функцией двух переменных — протяженности тела x и времени t: y(x, t). Графически функция представляется поверхностью (см. рис. 19.5) или набором изолиний (см. рис. 19.6), вид которых обычно требуется определить.

Рис. 19.5. Расчет распространения тепла в одномерном стержне со временем

Рис. 19.6. Распространение тепла по элементам двухмерной системы (расчет статической задачи, представленной на рис. 19.1). В отдельном окне видна оценка точности расчета в зависимости от количества итераций

Если заменить выражения производных их дискретным аналогом, то в разностном виде уравнение будет выглядеть так:

или, выражая неизвестное через известные величины:

или

В результате получена расчетная формула, реализуемая на цифровой вычислительной машине. Благодаря этой формуле можно рассчитать значение параметра y в любой точке (x, t).

Назовем значение y(x, t) узлом расчета. Тогда схематично расчет выглядит как сетка узлов на поле, составленном из частей тела и отрезков времени (см. рис. 19.7). Сама формула расчета одного узла зависит от состояния трех узлов (левого y(x – Δx, t – Δt), правого y(x + Δx, t – Δt), собственного y(x, t – Δt)) в предыдущий (t – Δt) момент времени и напоминает треугольный шаблон. До начала расчета известны состояния всех узлов для t = 0. Применяя формулу последовательно ко всем узлам для следующего момента времени, можно определить температуру во всех узлах следующего временного слоя (t + Δt). Кроме самого левого и самого правого узлов — их состояние вычислено быть не может, но оно задано краевыми условиями.

Рис. 19.7. Схема расчета одномерной динамической системы с распределенными параметрами методом конечных разностей с использованием явного шаблона

Если процедуру повторять, переходя от одной точки тела x к другой, и далее от одного временного слоя к другому, то по данной формуле можно вычислить значение температуры в любой части тела в любой момент времени. Таким образом, расчетом покрывается все поле (L x T) (см. рис. 19.7). Последовательное определение неизвестных значений в данном случае возможно, потому что шаблон имеет вид явного выражения — единственное неизвестное в формуле выражено через ранее вычисленные значения.

Заметим, что при больших значениях производных и больших значениях шагов расчет может дать неверные решения. Решения могут оказаться неточными или даже неустойчивыми (качественно неверными) (см. лекцию 10. «Численные методы интегрирования дифференциальных уравнений. Метод Эйлера»).

Условие устойчивости для треугольного шаблона при решении уравнения диффузии: Δx/Δt > α (см. подробнее рис. 19.12).

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

Рис. 19.8. Схемы типовых шаблонов для решения уравнений с распределенными параметрами

На практике устойчивости, а далее — точности добиваются получением решений с использованием разных шаблонов и разных значений шага. Если значения искомой переменной, вычисленные с шагом h и с шагом h/2, отличаются в узлах с одинаковыми индексами не более чем на 1—5%, то вычисленное значение принимают за приближенное решение задачи. Иначе уменьшают шаг еще в два раза, и процедуру оценки повторяют. (Дополнительно см. лекцию «Умеем ли мы вычислять на компьютере?».)

Свойства уравнения диффузии отражены на рис. 19.9 и заключаются в том, что при возникновении неоднородности в какой-то из частей тела со временем тепло за счет процессов теплообмена перетекает в соседние области. Температуры соседних областей выравниваются, усредняются. Темп процесса зависит от величины коэффициента теплопроводности.

Рис. 19.9. Характерный вид решения уравнения диффузии. На рисунке отражено изменение распределения параметра y вдоль оси x во времени. а) — в двух осях; б) — в трех осях

Если принять условие, что задача стационарная, то есть процессы протекают так долго, что все переходные процессы успели закончиться (производная по времени равна 0), то уравнение диффузии приобретает следующий вид (для случая двухмерного пространства — оси x и z) без источников и стоков:

∂²y/∂x² + ∂²y/∂z² = 0.

В разностном виде уравнение имеет вид:

(Y_i + 1, j – 2 · Y_i, j + Y_{i – 1, j})/Δx² + (Y_{i, j – 1} – 2 · Y_i, j + Y_i, j + 1)/Δz² = 0.

Если принять Δx = Δz, то уравнение примет вид:

4 · Y_i, j – Y_i + 1, j – Y_{i – 1, j} – Y_{i, j – 1} – Y_i, j + 1 = 0.

Легко понять, что шаблон расчета уравнения неявный и имеет вид креста (чтобы рассчитать значение температуры в узле сетки, надо знать температуры его соседей слева, справа, сверху и снизу). Если стена дома имеет размеры 2 метра на 2 метра, а шаг Δx = Δz = 20 мм, то всего для расчета температурного режима стены придется решать систему из 10 000 линейных уравнений c 10 000 неизвестных Y_i, j:

4 · Y_i, j – Y_i + 1, j – Y_{i – 1, j} – Y_{i, j – 1} – Y_i, j + 1 = 0, для i = 1÷100 и j = 1÷100,

к которым следует присоединить 400 штук краевых условий:
Y_0, j = f₁(j);
Y_101, j = f₂(j);
Y_i, 0 = f₃(i);
Y_i, 101 = f₄(i).

Вид решения уравнения показан на рис. 19.6.

Уравнение тепломассопереноса

В ряде физических процессов масса вещества или энергия движутся внутри воображаемой системы, перемещаясь из одних ее частей в другие. Такие процессы описываются уравнением тепломассопереноса. Для одномерного случая это уравнение имеет следующий вид: ∂ρ/∂t + ∂q/∂x = B, где ρ(x, t) — плотность вещества или энергии; q(x, t) — поток вещества или энергии; B(x, t) — функция источников и стоков. ρ(x, t), q(x, t), B(x, t) — функции распределения соответствующих переменных в пространстве x и времени t.

Разобьем поток вдоль оси x на элементарные ячейки, взаимодействующие друг с другом. Каждая ячейка характеризуется плотностью ρ находящегося в ней вещества и потоком вещества q через ее границы (левую и правую). Ячейки пронумеруем индексом i. Ячейку принято называть элементарным объемом.

Рис. 19.10. Схема уравнения теплопроводности (пример)

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

Уравнение и соответствующая ему схема на рис. 19.10 показывает, сколько вещества добавится в элементарном объеме ∂ρ/∂t за время Δt, если через левую границу этого участка x за это время перетекает q(x) · Δt частиц вещества, и если через правую границу этого участка (x + Δx) за это же время перетекает q · (x + Δx) · Δt частиц вещества, и на участке есть источник или сток вещества величиной b.

В зависимости от решаемой задачи уравнение может быть дополнено. Функция b(x, t) обычно задана. Если требуется определить закон движения вещества q(x, t), чтобы обеспечить заданное пользователем распределение вещества по оси x со временем ρ(x, t), то для такого расчета достаточно одного уравнения.

Если задан закон движения вещества q(x, t) и требуется определить, где в результате этого движения, когда и сколько будет вещества ρ(x, t), то также достаточно одного уравнения.

Но часто требуется определить одновременно и ρ(x, t) и q(x, t) при заданных начальных условиях, то есть вопрос состоит в том, как быстро q(x, t) и сколько ρ(x, t) вещества окажется в определенных точках x. Тогда одного уравнения для расчета двух неизвестных недостаточно. И основное уравнение должно быть дополнено вспомогательными выражениями, указывающими, как связаны между собой ρ(x, t), q(x, t). Такое выражение должно указывать, как быстро движется вещество, если задано его количество. В производстве и технике мы называем такое выражение регулятором: как связано количество деталей ρ, хранящихся на участке x, с желаемой скоростью обработки таких деталей q.

В зависимости от того, с каким веществом мы имеем дело, могут существовать различные выражения дополнительной связи (дополнительное уравнение) между потоком q и плотностью ρ. Например, для жидкости, где все частицы движутся одновременно: q = ρ · v, где v — скорость перемещения частиц в потоке. Или более сложное уравнение движения, связывающее скорость движения потока v(x, t) с причиной этого движения, силой F:

∂v/∂t + v · ∂v/∂x = F/ρ.

В свою очередь, сила F может быть вызвана к действию перепадом давления P: F = ∂P/∂x.

Легко увидеть в данном уравнении в качестве его основы второй закон Ньютона, связывающий ускорение ∂v/∂t и силу F через массу ρ. Дополнительное слагаемое v · ∂v/∂x появляется постольку, поскольку при движении появляется приток вещества через границы x и (x + Δx) элементарного участка, который несет за собой импульс v · ∂v/∂x · ρ.

В разностном виде при q = ρ · v уравнение переноса будет иметь вид:

Расчет уравнения производится аналогично уравнению диффузии.

Свойства уравнения переноса заключаются в перемещении неоднородностей в потоке со временем вдоль пространственной оси. Скорость перемещения связана с величинами v или q. На рис. 19.11 показан вид одного из возможных вариантов решения уравнения тепломассопереноса. Видно, что распределение плотности (неоднородность) может перемещаться вдоль оси со временем.

Рис. 19.11. Свойства уравнения тепломассопереноса. а) — в двух осях; б) — в трех осях

Особое внимание при расчете уравнения следует обращать на его устойчивость. Устойчивость связана с видом выбранного разностного шаблона и его параметрами: шагом по t, шагом по x, а также скоростью (потоком) перемещения массы (см. рис. 19.12). Чем больше шаг, тем больше риск получения неустойчивых (качественно неверных) решений. Но чем меньше шаги, тем медленнее происходит компьютерный расчет всего поля, так как объем вычислений возрастает значительно.

Рис. 19.12. Соотношение параметров расчета для получения устойчивого решения

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

Рис. 19.13. Вид возможных решений уравнения тепломассопереноса с элементами диффузии. На рисунке видно проявление нескольких свойств одновременно

К распределенным системам относятся производственные системы (см. лекцию 31. «Моделирование производственных процессов и систем»).

Одной из самых сложных задач, рассматриваемых в классе систем с распределенными параметрами, является прогноз погоды. При этом учитываются температура, давление, скорость ветра, влажность, плотность воздуха, рельеф местности. Все эти величины являются функцией трех пространственных координат и времени. Ежедневно анализу в центрах прогнозирования погоды подвергается 10⁶ исходных данных, не считая фотографий. Расчет производится два раза в сутки: в 00 часов и в 12 часов по Гринвичу. Обмен измерительной информацией между метеостанциями разных стран идет в рамках Всемирного метеорологического общества, куда входит и Российская Федерация. В среднем для расчета одного прогноза учитывается информация от 4 000 приземных источников (включая 800 морских) и 650 аэроисточников. Сеть измеряющих приборов покрывает земной шар неравномерно и это осложняет расчеты.

Нетрудно оценить в среднем шаг расчетной сетки. При радиусе Земли R = 6 400 км площадь Северного полушария равна: S = 260 000 000 км². Таким образом, на одну аэростанцию падает S/650 = 400 000 км² или квадрат со стороной, равной расстоянию от Москвы до Петербурга. На одну наземную станцию падает квадрат со стороной 300 км. Такой шаг сетки расчета ведет к огромным трудностям при борьбе за устойчивость расчета.

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

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