Явная схема Эйлера

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

Построение разностной схемы

Используем для решения уравнения теплопроводности шаблон, изображенный на рис. 13.6. Для дискретизации второй производной по пространственной координате необходимо использовать три последовательных узла, в то время как для разностной записи первой производной по времени достаточно двух узлов. Записывая на основании данного шаблона дискретное представление для (i,k)-ro узла, получим разностное уравнение:

Рис. 13.6. Шаблон аппроксимации явной схемы для уравнения теплопроводности

Приведем в разностной схеме (8) подобные слагаемые, перенеся в правую часть значения сеточной функции с индексом k (как часто говорят, с предыдущего слоя по времени), а в левую - с индексом k-t-i (т. е. со следующего временного слоя). Кроме этого, введем коэффициент с , который будет характеризовать отношение шагов разностной схемы по времени и пространству Несколько забегая вперед, заметим, что значение параметра с, называемого коэффициентом Куранта, имеет большое значение для анализа устойчивости разностной схемы. С учетом этих замечаний, разностная схема (8) запишется в виде:

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

Несложно убедиться в том, что для получения замкнутой системы разностных алгебраических уравнений систему (9) необходимо дополнить дискретным представлением начального и граничных условий (6) и (7). Тогда число неизвестных будет в точности равно числу уравнений, и процесс формирования разностной схемы будет окончательно завершен.

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

Если присмотреться к разностным уравнениям (9) повнимательнее, то можно сразу предложить несложный алгоритм реализации этой разностной схемы. Действительно, каждое неизвестное значение сеточной функции со следующего временного слоя, т. е. левая часть соотношения (9) явно выражается через три ее значения с предыдущего слоя (правая часть), которые уже известны. Таким образом, в случае уравнения теплопроводности нам очень повезло — для расчета 1-го слоя по времени следует попросту подставить в (9) начальное условие (известные значения и с нулевого слоя в узлах сетки), для расчета 2-го слоя достаточно использовать вычисленный таким образом набор и с 1-го слоя и т. д. Из-за того, что разностная схема сводится к такой явной подстановке, ее и называют явной, а благодаря пересчету значений с текущего слоя через ранее вычисленные слои — схемой бегущего счета.

Линейное уравнение

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

Решение системы разностных уравнений (9) для модели без источников тепла, т. е. ф(х, T, t)=o и постоянного коэффициента диффузии D=const приведена в листинге 13.1. В его первых трех строках заданы шаги по временной и пространственной переменным т и Д, а также коэффициент диффузии D, равный единице. В следующих двух строках заданы начальные (нагретый центр области) и граничные (постоянная температура на краях) условия, соответственно. Затем приводится возможное программное решение разностной схемы, причем функция пользователя v(t) задает вектор распределения искомой температуры в каждый момент времени (иными словами, на каждом слое), задаваемый целым числом t.

Листинг 13.1. Явная схема для линейного уравнения теплопроводности

Рис. 13.7. Решение линейного уравнения теплопроводности (листинг 13.1)

Начальное распределение температуры вдоль расчетной области и решение для двух моментов времени показано на рис. 13.7 сплошной, пунктирной и штриховой линиями, соответственно. Физически такое поведение вполне естественно — с течением времени тепло из более нагретой области перетекает в менее нагретую, а зона изначально высокой температуры остывает и размывается.

Нелинейное уравнение

Намного более интересные решения можно получить для нелинейного уравнения теплопроводности, например с нелинейным источником тепла Ф(u)=10³(u-u³). Заметим, что в листинге 13.1 мы предусмотрительно определили коэффициент диффузии и источник тепла в виде пользовательских функций, зависящих от аргумента и, т. е. от температуры (если бы собирались моделировать явную зависимость их от координат, то следует ввести в пользовательскую функцию в качестве аргумента переменную х, как это сделано для источника тепла ф). Поэтому нет ничего проще замены определения этих функций с констант D(u)=1 и ф(х,u)=0 на новые функции, которые станут описывать другие модели диффузии тепла. Начнем с того, что поменяем четвертую строку листинга 13.1 на ф(х,u)=10³(u-u³), не изменяя пока постоянного значения коэффициента диффузии.

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

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

Еще более неожиданные решения возможны при нелинейности также и коэффициента диффузии. Например, если взять квадратичный коэффициент диффузии D(x,u)=u² (что с учетом его умножения на неизвестные функции создаст кубическую нелинейность уравнения), а также ф(х,u)=10³u^{3 5}, то Вы сможете наблюдать совсем иной режим горения среды. В отличие от рассмотренного эффекта распространения тепловых фронтов, горение оказывается локализованным в области первичного нагрева среды, причем, температура в центре нагрева со временем возрастает до бесконечной величины (рис. 13.9). Такое решение описывает так называемый режим горения «с обострением».

Рис. 13.8. Решение уравнения теплопроводности с нелинейным источником (тепловой фронт)

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

Рис. 13.9. Решение уравнения теплопроводности с нелинейным источником и коэффициентом диффузии (режим локализации горения)

Устойчивость

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

Слегка изменим соотношение шагов по времени и пространственной координате, произведя расчеты сначала с т=0.0005 (эти результаты показаны на рис. 137 выше) и также с т=0.0010 и т=0.0015. Результат применения разностной схемы Эйлера для т=0.0010 примерно тот же, что и для меньшего значения т, приведенного на рис. 13.7. А вот следующее (казалось бы, незначительное) увеличение шага по времени приводит к катастрофе (рис. 13.10). Вместо ожидаемого решения получается совершенно неожиданные профили температуры, которые быстро осциллируют вдоль пространственной координаты, причем амплитуда и число пиков этих осцилляции быстро увеличиваются от шага к шагу. Совершенно ясно, что полученное решение не имеет ничего общего с физикой моделируемого явления, а является следствием внутренних свойств самой разностной схемы, которые до этого были для нас скрыты.

Рис. 13.10. Численное решение уравнения теплопроводности при помощи явной схемы Эйлера (см. листинг 13.1 ниже с временным шагом т=0.0015)

Характерная «разболтка» решения как раз и является проявлением неустойчивости явной схемы Эйлера для выбранного соотношения шагов по времени и пространству. В теории численных методов показывается, что явная схема Эйлера для уравнения теплопроводности устойчива при значениях коэффициента Куранта, меньших 1, и неустойчива в противоположном случае. Иными словами, существует ограничение для выбора соотношения шагов, заключающееся в том, что для расчета на более частых пространственных сетках необходимо использовать также и малые шаги по времени

Как несложно убедиться, для т=0.0005 коэффициент Куранта C=0.4, для т=0.0010 он все еще меньше единицы C=0.8, а для т=0.0015 решение уже больше единицы: C=1.2, в связи с чем схема становится неустойчивой (рис 13.10).