в частных производных

8.3.1. Построение конечно-разностных схем

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

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

;

, .

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

;

, ;

или

;

, .

Аппроксимация уравнений может осуществляться по явным или по неявным конечно-разностным схемам. При записи конечно-разностных схем полагается, что при определении параметров функции f в момент времени известна информация о значении сеточной функции (соответствует моменту времени ) для всех значений пространственной координаты ().

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

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

Рассмотрим наиболее известные варианты конечно-разностных схем, используемые при решении задач математической физики (например, [Андерсон,Дьяченко,Емельянов, ГодуРябенький, Роуч, Флетчер.]). Для упрощения будем рассматривать уравнения в частных производных, записанных в нестационарной одномерной постановке.

Уравнение переноса

(8.21)

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

1. Явные схемы первого порядка точности по времени и по пространственной координате (шаблоны, представленные на рисунках 8.7а, 8.7б):

; (8.22)

. (8.23)

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

2. Схема «тренога» - явная схема первого порядка точности по времени и второго порядка по пространственной координате (шаблон представлен на рисунке 8.7в):

; (8.24)

3. Схема Лакса - явная схема первого порядка точности по времени и второго порядка по пространственной координате (шаблон представлен на рисунке 8.7г):

; (8.25)

4. Схема «чехарда» - явная схема второго порядка точности по времени и второго порядка по пространственной координате (шаблон представлен на рисунке 8.7д):

; (8.26)

5. Схема Лакса-Вендрофа - явная схема типа «предиктор-корректор» второго порядка точности по времени и второго порядка по пространственной координате:

первый шаг (предиктор) – осуществляется по формуле Лакса (8.25);

второй шаг (корректор) – выполняется по формуле

; (8.27)

6. Схема Мак-Кормака - явная схема типа «предиктор-корректор» второго порядка точности по времени и второго порядка по пространственной координате:

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

;

на втором шаге (корректор) вычисляется окончательное значение сеточной функции по формуле вида (8.23)

. (8.28)

7. Неявная схема первого порядка точности по времени и второго порядка по пространственной координате (шаблон представлен на рисунке 8.7е):

(8.29)

Уравнение диффузии (теплопроводности)

. (8.30)

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

1. Явная схема первого порядка точности по времени и второго порядка по пространственной координате (шаблон представлен на рисунке 8.8а):

; (8.31)

2. Явная схема Дюфорта-Франкела первого порядка точности по времени и второго порядка по пространственной координате (шаблон представлен на рисунке 8.8б):

; (8.32)

3. Неявная конечно-разностная схема запишется в виде (шаблон представлен на рисунке 8.8в):

; (8.33)

4. Неявная конечно-разностная схема Кранка-Николсона (шаблон представлен на рисунке 8.8г):

. (8.34)

В формуле (8.34) - . В частности, можно принять .

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

Рассмотрим решение уравнения переноса, выполняемого по явной схеме (8.22). Уравнение приводится к виду

. (8.35)

Определение сеточной функции по уравнению (8.35) устанавливается следующим алгоритмом:

- примем значение n= 0 и установим шаг интегрирования , обеспечивающий устойчивость последующих вычислений (порядок выбора шага интегрирования уравнений по времени, обеспечивающих устойчивый расчет конечно-разностных уравнений, рассматривается ниже в разделе 8.3.3);

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

- в соответствии с граничными условиями устанавливаем значение сеточной функции ;

- начиная с до в соответствии с уравнением (8.35) вычисляем значения ;

- увеличиваем значение n на единицу и возвращаемся к началу алгоритма.

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

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

. (8.36)

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

;

;

. (8.37)

В записанной системе уравнений граничные условия, соответственно, на левой и правой границах расчетной области (рисунок 8.8?). Количество уравнений совпадает с числом неизвестных значений сеточной функции , поэтому система уравнений имеет единственное решение, если известны значения сеточной функции , соответствующие моменту времени . Система уравнений решается многократно, начиная со значения (соответствует моменту времени ), до значения (соответствует моменту времени ). При решение системы уравнений (8.37) осуществляется по заданным начальным условиям (соответствуют сеточным значениям функций ). Решение системы линейных уравнений (8.37) на любом очередном шаге интегрирования может выполняться методами, изложенными в главе 2 (раздел 2.2.2). Однако для подобных систем существуют более эффективные алгоритмы, и некоторые из этих алгоритмов будут рассмотрены ниже.

Следует дополнительно отметить следующее.

1. Выше рассматривались варианты конечно-разностной аппроксимации уравнений, в которых присутствуют частные производные по времени и по пространственным координатам от одной и той же функции. Однако на практике возникает необходимость в решении уравнений вида (8.7), в которых каждая из частных производных, входящая в уравнение, определяется от своей функции. Эти функции, как правило, связаны между собой. Для подобных уравнений применимы все рассмотренные выше приемы конечно-разностной аппроксимации. Например, при решении уравнения

(8.38)

можно применить схему Лакса в следующей конечно-разностной форме

. (8.39)

В формулах (8.38), (8.39) g – может быть функцией переменных t, x () или t, x, f ().

2. Аналогичные рассмотренным могут быть построены явные конечно-разностные схемы решения многомерных задач, описываемых уравнениями в частных производных. Сеточную функцию в этом случае обозначают, например, в виде . Здесь i, j, k – целочисленные координаты сеточной функции по трем координатным направлениям.

Рассмотрим трехмерный вариант уравнения переноса

и запишем для его решения конечно-разностную явную схему Лакса

.

Решение последнего уравнения остается столь же простым, что и решение уравнения (8.25).

3. Применение неявных схем для многомерных задач существенно усложняет их решение. Рассмотрим, например, многомерный вариант уравнения диффузии

. (8.40)

Для этого уравнения можно записать следующую неявную конечно-разностную схему

(8.41)

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

.

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

Существуют эффективные вычислительные методы решения многомерных задач, называемые методами пространственного расщепления, методами дробных шагов, суть которых состоит в сведении решения многомерной задачи к совокупности одномерных задач. Для уравнения (8.40) применение этого метода требует решения трех систем уравнений – (8.42), (8.43), (8.44).

, (8.42)

, (8.43)

. (8.44)

Каждая из систем уравнений вместе с начальными и граничными условиями может быть приведена к виду (8.33) с трехдиагональной матрицей А.

4. Исходное дифференциальное уравнение в частных производных можно привести к виду системы обыкновенных дифференциальных уравнений. Например, для случая уравнения диффузии (8.30) при явной аппроксимации второй производной по пространственной координате такая система может быть записана в виде

. (8.45)

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

5. В последние годы особый интерес приобретают конечно-разностные схемы высоких порядков точности по времени и по пространственным координатам. Построение таких схем рассматривается, например, в работах [Марчук, МарчукШайдуров, Липанов].

8.3.2. Анализ сходимости конечно-разностных схем

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

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

Рассмотрим примеры, в которых будем использовать равномерную по времени и по пространственной координате сетку (, рисунок 8.9).

В соответствии с формулами (2.30), (2.31) и [Шуп] можно записать следующие соотношения

;

.

Выше записанные соотношения позволяют установить значения производных, записанных в конечно-разностном виде, с учетом остаточных членов, полученных из разложений Тэйлора:

- аппроксимация производной конечной разностью первого порядка точности по времени

; (8.46)

- аппроксимация производной конечной разностью по схеме Лакса

; (8.47)

- аппроксимации производной конечными разностями первого и второго порядков точности по пространственной координате

; (8.48)

; (8.49)

; (8.50)

- аппроксимации производной конечной разностью второго порядка точности по пространственной координате

; (8.51)

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

. (8.52)

Соотношения (8.46)…(8.52) можно применить для оценки точности аппроксимации уравнений по различным конечно-разностным схемам. В уравнениях, записываемых ниже, в силу произвольности времени и пространственной координаты опускаются индексы i и n:

- схема (8.22) для уравнения переноса (8.21)

;

- схема (8.23) для уравнения переноса (8.21)

.

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

- схема Лакса (8.25) для уравнения переноса (8.21)

;

- схема (8.31) для уравнения диффузии (8.30)

;

- схема Дюфорта-Франкела (8.32) для уравнения диффузии (8.30)

. (8.53)

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

Из выше приведенных конечно-разностных аппроксимаций особый интерес представляют аппроксимация по схеме Лакса (она содержит остаточные члены, содержащие отношение ) и аппроксимация по схеме Дюфорта-Франкела (она содержит остаточные члены, содержащие отношения ).

Аппроксимация дифференциальных уравнений с порядком - это одно из условий сходимости конечно-разностного уравнения к дифференциальному. Второе условие – это устойчивость вычислительного алгоритма, заложенного в основу конечно-разностной схемы. Вопросы устойчивости рассматриваются, например, в [СамарскийГулин]. Ниже рассматриваются отдельные методы оценки устойчивости в применении к модельным уравнениям переноса (8.21) и диффузии (8.30).

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

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

Анализ устойчивости с использованием степенной функции

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

С учетом выше записанного условие устойчивости может быть представлено в виде

. (8.54)

Рассмотрим применение изложенного подхода для решения уравнения переноса (8.21), записанного в виде (8.22). Для этого уравнения значение устанавливается формулой

.

Подставим в уравнение значения , и установим значение сеточной функции

.

Условие (8.47) будет выполняться лишь в том случае, если

или

.

Окончательно условие устойчивости запишется в виде

. (8.55)

Из условия (8.55) следует важный вывод – конечно-разностная схема (8.22) устойчива лишь в случае, если значение коэффициента a в уравнении переноса (8.21) положительно (). Можно показать, что в случае отрицательных значений коэффициента a () следует применять конечно-разностную схему (8.23), а условие устойчивости в этом случае запишется в виде

. (8.56)

Рассмотренный метод относительно прост, однако и возможности его ограничены. В частности, применение метода для схем (8.24), (8.25) не позволит получить практически полезного результата.

Спектральный метод (метод Неймана) оценки устойчивости

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

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

Подстановка этого выражения в уравнение переноса (8.21), записанное в виде (8.22), позволяет записать следующее

.

После упрощений последнее уравнение может быть переписано в виде

.

Из этого следует зависимость для q

. (8.57)

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

Рассмотрим другие примеры:

- аппроксимация уравнения переноса (8.21) схемой «тренога» (8.24)

.

После упрощений формула для значения q имеет вид

и одно из решений этого уравнения - означает, что схема «тренога» абсолютно неустойчива;

- аппроксимация уравнения переноса (8.21) неявной схемой (8.29)

.

После упрощений получаем

или

.

Из последнего соотношения следует, что при любых значениях и при любых значениях отношение . Этот факт свидетельствует о том, что неявная конечно-разностная схема (8.29) всегда устойчива к накоплению вычислительных ошибок.

Конечно-разностные схемы, устойчивые при любых значениях , называются абсолютно устойчивыми. Обычно неявные конечно-разностные схемы являются абсолютно устойчивыми. Однако за это их свойство приходится расплачиваться применением (в сравнении с явными конечно-разностными схемами) более сложных вычислительных алгоритмов;

- аппроксимация уравнения диффузии (8.30) явной схемой (8.31)

.

После упрощений выражение для q принимает вид

,

а условие устойчивости вычислительного алгоритма () запишется как

.

Из последнего уравнения следует, что конечно-разностная схема (8.31) условно устойчива, а условие устойчивости записывается формулой

;

- аппроксимация уравнения диффузии (8.30) неявной схемой (8.33)

После упрощений выражение для q принимает вид

- аппроксимация уравнения диффузии (8.30) схемой Дюфорта-Франкела (8.32)

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

.

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

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

Метод анализа устойчивости с использованием собственных значений

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

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

(8.58)

можно представить в виде суммы [Камке]

. (8.59)

В уравнении (8.58) приняты обозначения:

; .

В уравнении (8.59) собственное значение матрицы, которое может быть найдено решением алгебраического уравнения

.

Из (8.59) следует, что решение уравнения (8.58) будет убывающим лишь в том случае, когда все собственные значения по абсолютному значению будут меньше единицы .

В уравнении (8.58) под вектором y можно понимать как решения для функции f, определяемые решением конечно-разностных уравнений, так и ошибки вычислений (здесь - соответственно численное и точное решения уравнений). В обоих случаях матрица A для этих уравнений будет одна и та же.

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

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

Рассмотрим пример из практики [АлиевСуворов].

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

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

(8.60)

Здесь в уравнении (8.60) и ниже прияты следующие обозначения:

r, с, k – плотность, удельная теплоемкость, коэффициент теплопроводности; T – температура; τ – время; r – радиальная координата; z – осевая координата; Q – тепловой эффект реакции; j – доля прореагировавшего вещества; K ₀ – предэкспонент; E _а – эффективная энергия активации процесса разложения; R – газовая постоянная.

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

Для решения задачи теплопроводности вся расчетная область разбивается на совокупность элементарных объемов (ниже нумеруются двумя индексами и ), для каждой из которых выполняется конечно-разностная запись уравнения теплопроводности.

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

. (8.61)

Здесь обозначено:

; ;

; ; .

Конечно-разностная аппроксимация уравнений (8.61) выполняется по явной схеме, что не обеспечивает абсолютной устойчивости вычислительного алгоритма. Из-за неоднородности теплофизических свойств материала анализ устойчивости с использованием, например, спектрального метода представляется трудоемким. Для анализа устойчивости применяемого вычислительного алгоритма наряду с задачей (8.61) будем рассматривать задачу для ошибки вычислений , определяемую разностью его численного решения для температуры и точного значения в произвольный момент времени ( - характерное время в решаемой задаче).

С учетом (8.61) система уравнений для ошибки численного решения () запишется в виде:

(8.62)

В уравнениях (8.62) индекс m связан с индексами i и j по формуле , а значения коэффициентов - по формуле .

Система уравнений (8.62) – это система обыкновенных линейных дифференциальных уравнений. Можно утверждать, что выбранная конечно-разностная аппроксимация уравнения теплопроводности будет обеспечивать устойчивость его решения на интервале времени (или ) в случае, когда модули всех собственных значений матрицы А, составленной из коэффициентов , будут меньше единицы. Такой вывод позволяет свести задачу об анализе устойчивости численного решения уравнения теплопроводности к задаче о нахождении такого значения времени , при котором максимальные по модулю собственные значения матрицы А будут меньше единицы. В рассматриваемой задаче в качестве характерного времени целесообразно использовать значение шага интегрирования уравнений по времени процесса . Для вычисления собственных значений матрицы А могут быть применены любые известные вычислительные методы, в частности, ортогональные методы равносильных преобразований [Бахвалов,Вежбицкий].

Проверка разработанных вычислительных алгоритмов выполнена решением нескольких задач. В частности, выполнены расчеты нагрева однородного материала и многослойных пластин (рисунок 8.11). Нагрев материалов осуществлялся внешним тепловым потоком при отсутствии внутренних источников тепла. Решение осуществлялось на нерегулярной сетке с элементарными объемами, вид которых представлен на рисунке 8.12. В первом варианте (рисунок 8.12а) использовались элементарные объемы квадратного сечения, во втором (рисунок 8.12б) и в третьем (рисунке 8.1?2в) вариантах – треугольного сечения. Центры элементарных объемов (отмечены точками) располагались в центрах тяжести сечений. В расчетах использовалось от 60 до 20000 конечных объемов.

В расчетах неоднородного материала использовалась двухслойная композиция. Первый материал в композиции – алюминий (теплофизические свойства - k = 237 (Дж/м×с×К), c = 900 (Дж/кг×К), r = 2699 (кг/м³)). Второй материал – порох «Н» (теплофизические свойства - k = 0.234 (Дж/м×с×К), c = 1465 (Дж/кг×К), r = 1600 (кг/м³), температура зажигания – 618 К). Решение выполнялось в предположении отсутствия твердофазных химических реакций и фазовых переходов. В первоначальный момент времени температура сложного материала принималась одинаковой во всей расчетной области - К. Нагрев материала осуществлялся на внешней границе материала конвективным тепловым потоком в результате воздействия высокотемпературного потока газа (температура газа - 2500 К), значение коэффициента теплоотдачи α принималось постоянным Вт/м²К.

Результаты выполненных тестовых расчетов представлены на рисунках 8.13…8.15.

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

На рисунке 8.14 представлены результаты расчета прогрева многослойной пластины при использовании элементарных объемов разных типов. На рисунках кривая 1 соответствует расчетам, в которых элементарные объемы имеют прямоугольное сечение (рисунке 8.10а). Кривая 2 соответствует расчетам, в которых элементарные объемы имеют треугольное сечение вида, представленного на рисунке 8.10б. Кривая 3 соответствует расчетам с треугольными элементарными объемами вида, представленного на рисунке 8.10в. Кривая 4 соответствует аналитическому расчету многослойной пластины на стационарном режиме [ИсаченкоТеплопередача]. На рисунке 8.14а представлены результаты, соответствующие моменту времени с. На рисунке 8.14б – моменту времени с.. Анализ показывает, что погрешность расчетов во всех рассмотренных случаях обеспечивается на допустимом для практических расчетов уровне (погрешность не превосходит 2%).

Условие Куранта-Фридрихса-Леви

Оценка устойчивости конечно-разностных схем с использованием условия Куранта-Фридрихса-Леви [Флетчер, Роуч] применяется на практике часто, несмотря на то, что это только необходимое условие устойчивости. Смысл условия прост – реальные и моделируемые физические процессы должны быть согласованы. В частности, если в какой-то точке расчетной области произошло возмущающее воздействие, то скорость его распространения в конечно-разностной схеме не может превосходить физической скорости. Пусть, например, скорость среды в какой-либо ячейке составляет u, длина ячейки . В этом случае время распространения возмущения от одной границы до другой будет составлять . Условие Куранта-Фридрихса-Леви требует, чтобы реально выбираемый в расчетах конечно-разностных задач шаг интегрирования по времени не превосходил значения . Это условие записывается в виде

. (8.63)

В формуле (8.63) Ku – число Куранта. Обычно в расчетах принимается .

В задачах газовой динамики возмущение распространяется со скоростью (). Здесь c – скорость звука в газе. Для задач газовой динамики, как правило, условие устойчивости Куранта-Фридрихса-Леви записывается в виде

. (8.64)

В общем случае для системы дифференциальных уравнений в частных производных вида

условие Куранта-Фридрихса-Леви записывается

. (8.65)

В формуле (8.65) - максимальное по абсолютному значению собственное значение матрицы А.

Применение искусственной вязкости

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

Как пример можно привести вариант простейшей конечно - разностной схемы первого порядка точности В.В. Русанова [Роуч], применявшейся на ранних этапах развития вычислительной гидромеханики. В соответствии с методом В.В. Русанова вместо уравнения (8.38) решается уравнение вида

.

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

; ; .

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