4.1 Метод конечных разностей

Введенное в предыдущем параграфе определение сеточной функции естественным образом распространяется на случай двух и более аргументов. Простейшую двумерную сетку {x_i,j, y_i,j} можно построить как набор из J одномерных сеток, каждая из которых содержит I узлов, при этом должно выполняться условие:

Очевидно, что суммарное число сеточных узлов в этом случае будет IJ. Расчетные сетки, удовлетворяющие условию (19), называются регулярными (упорядоченными). Для упрощения записи в дальнейшем мы будем использовать обозначение Ф_i,j= (х_i,j, у_i,j). Пример сеточной функции двух переменных, построенной на регулярной сетке, приведен на рис. 5.

Рис. 5. Сеточная функция двух переменных, заданная на упорядоченной сетке

Частные производные функции могут быть аппроксимированы аналогично (9). Например, в случае Δx=const, Δy=const¹⁹ при использовании центрального конечно-разностного отношения, будем иметь:

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

где Т – температура (подлежащая определению в результате решения), t – время, α – коэффициент теплопроводности, х – координата.

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

С математической точки зрения, решением уравнения (21) является функция, зависящая от двух переменных (координата х и время t).

Рис. 6. Стержень с теплоизолированной боковой поверхностью

Для решения задачи должно быть заданы:

• граничное условие на левом торце стержня – либо зависимость температуры от времени T₀=T(x₀,t) (условие Дирихле) либо пространственная производная (условие Неймана);

• граничное условие на правом конце стержня – T_N=T(x_N,t) либо ;

• начальное условие – закон распределения температуры по всей длине стержня в начальный момент времени Т(х, 0), где x₀<x< x_N.

Заменим²⁰ производную по времени в уравнении теплопроводности (21) разностным отношением вперед (см. 9b), а вторую производную по пространственной координате – центральным конечно-разностным отношением (см. 20)²¹:

где i – номер рассматриваемого сеточного узла, j – номер рассматри­ваемого шага по времени, Δx – шаг пространственной сетки, Δt – шаг по времени, T_i,j – температура в сеточном узле с координатой x_i=x₀+i∙Δx в момент времени t_j= j∙Δt.

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

Рис. 7. Расчетная сетка, используемая для решения уравнения теплопроводности (21)

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

Алгебраическое уравнение (21) позволяет методом "последовательного обхода" узлов расчетной сетки найти все неизвестные значения сеточной функции T_{i, j}:

Как следует из рассмотрения рис. 8а, в правой части уравнения (22a) упоминаются лишь те узлы расчетной сетки, значения температуры в которых уже известны. Расчетные схемы, обладающие указанным свойством, называются явными. Данное обстоятельство, с одной стороны, существенно упрощает решение задачи об отыскании значений сеточной функции Т_{i, j}, с другой стороны, анализ формулы (22a), показывает, что температура в узле x_i в момент времени t_j+1, зависит только от температуры в узлах x_i-1, x_i, x_i+1 в момент времени t_j, и не зависит от распределения температур внутри стержня в момент времени t_j+1, что не вполне соответствует физическому смыслу задачи.

Для устранения этого противоречия, несколько видоизменим аппроксимацию производной в уравнении (21):

или:

Правая часть уравнения (23a) содержит неизвестные величины T_i-1,j+1 и T_i+1,j+1, поэтому оно не может быть непосредственно (без увеличения объема вычислительной работы) использовано для нахождения температуры T_i,j+1. Расчетный шаблон для неявной расчетной схемы (23a) приведен на рис. 8б. Как показывают оценки [Флетчер], объем вычислений, при применении неявной схемы вместо явной схемы (22a), возрастает, ориентировочно, в два раза. Однако лучшая физическая обоснованность неявной схемы позволяет нам рассчитывать на более высокое качество результатов расчета.

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

где верхний индекс означает номер итерации.

Тогда, с учетом (23а), в первом приближении температура в момент времени t_j+1 может быть определена как

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

Следует также предусмотреть прекращение расчета в том случае, если решения не удастся достигнуть (решение не сойдется) после некоторого "разумного" числа итераций К. Причиной отсутствия сходимости, в частности, могут стать завышенные требования к точности решения ε, чрезмерно большие шаги расчетной сетки Δx и Δt, неудачный выбор начального приближения, а также некоторые физические эффекты²³.

В некоторых случаях, добиться сходимости решения удается за счет использования метода релаксации. Преобразуем последнюю из формул (25) к виду:

где r – коэффициент релаксации (r>0).

Очевидно, что, при r=1, формула (25а) совпадает с последней из формул (25). При значениях 0<r<1 (нижняя релаксация) число итераций, необходимое для достижения сходимости возрастает²⁴, однако при этом снижается вероятность получения расходящегося решения. И наоборот, использование коэффициента релаксации r>1 (верхняя релаксация), может привести к сокращению продолжительности расчета, а может вызвать полный "развал" численной схемы.

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