Численный метод

Аппроксимация данного дифференциального уравнения в частных производных с заданными граничными и начальными условиями выполняется с использованием неявной конечно-объёмной схемы. Используется регулярная трехмерная сетка из ячеек-параллелепипедов с индексами , где IM, JM, KM - количество ячеек по каждой оси координат. – координаты центров ячеек.

Интегрирование по объёму ячейки приводит к следующей схеме:

где l - индекс стороны ячейки, — единичный вектор внешней нормали. В случае параллелепипедной сетки каждый компонент может быть равен 0, +1 или -1, и ровно два компонента из трёх всегда равны 0 для любой стороны ячейки.

Как отмечалось ранее, уравнение нелинейное. Для его решения используется подход предиктор-корректор. На шаге предиктора мы аппроксимируем потоки в правой части, обозначенные верхним индексом P, используя значения k и dT/dE на текущем временном шаге n, но используем значения энтальпии на следующем (промежуточном) шаге n+1/2:

На шаге корректора мы аппроксимируем потоки в правой части, обозначенные верхним индексом C, используя значения k и dT/dE на промежуточном временном шаге n+1/2, которые могут отличаться от значений на временном шаге n, тогда как значения энтальпии берутся на следующем временном шаге n+1:

Здесь — объём ячейки (i,j,k); — площадь l-й стороны ячейки (i,j,k). Для ячеек на границе расчётной области потоки известны из граничных условий.

Описанная неявная схема, допускающая значительно больший шаг по времени Δt, чем позволяет условие устойчивости для явной схемы, приводит к необходимости решения больших разреженных систем линейных алгебраических уравнений вида Aq = b. Здесь неизвестный вектор q обозначает все неизвестные значения энтальпии E, известный вектор правой части b зависит от известных значений энтальпии на текущем шаге по времени, а A — матрица разреженной системы. Для решения этой системы можно использовать итерационный метод. Одно из его преимуществ заключается в том, что нет необходимости явно формировать матрицу: достаточно уметь эффективно вычислять произведение матрицы на вектор. Мы использовали безматричный итерационный метод сопряженных градиентов для несимметричных матриц. Его вычислительные формулы имеют вид (с дополнительными векторами s и g), где угловые скобки 〈.,.〉 обозначают скалярное произведение: