Обзор
Самый общий подход к численному решению набора связанных PDE (partial differential equation – Дифференциальное уравнение в частных производных) показан на рис. 2.11. В случае гидродинамики это уравнения, введенные ранее. Численное решение этих задач называется вычислительной гидродинамикой (CFD).
Рис. 2.11. Обзор общего подхода в вычислительной гидродинамике
Основным аспектом решения PDE является аппроксимация задействованных производных. Эта аппроксимация выполняется в дискретизированной области и приводит к набору алгебраических уравнений, в конечном итоге к системе линейных уравнений. Существует множество подходов к аппроксимации с различными специализациями и свойствами, где эти три являются наиболее заметными:
Поскольку FDS является решателем на основе FDM, основное внимание в этом разделе уделяется методу конечных разностей. Тем не менее, будет также рассмотрено несколько отдельных аспектов FVM.
Аппроксимация производных
Ряд Тейлора
Ряд Тейлора может аппроксимировать любую C^\infty функцию в точке разложения x_0. Аппроксимация f(x) производится в зависимости от h, являющейся переменной координатой вблизи x_0.
f(x_0 + h) =\sum_{i=0}^{\infty}\frac{1}{i!}\left.\frac{d^if}{dx^i}\right|_{x_0}\cdot h^i =
= f(x_0) + f'(x_0) \cdot h + \frac{1}{2} f''(x_0) \cdot h^2 + \frac{1}{6}f'''(x_0) \cdot h^3 + \cdots \quad (2.11)
На практике расширение прерывается в заданном порядкеr \mathcal{O}(h^n),где e n указывает количество рассмотренных терминов Тейлора, т.е.
f(x_0 + h) =\sum_{i=0}^{n-1}\left(\frac{1}{i!}\left.\frac{d^if}{dx^i}\right|_{x_0}\cdot h^i\right) + \mathcal{O}(h^n)\quad (2.12)
Обозначение Большое "О" \mathcal{O} указывает на асимптотическое поведение функции. При приближении производных, показанных в этом разделе, это указывает на функциональное уменьшение ошибки аппроксимации для h \rightarrow 0.
Расширение до третьего порядка \mathcal{O}(h^3) имеет следующую форму:
f(x_0 + h) = f(x_0) + f'(x_0)\cdot h + \frac{1}{2} f''(x_0)\cdot h^2 + \mathcal{O}(h^3)\quad (2.13)
Пример разложения Тейлора
В качестве примера приведена аппроксимация функции
f(x) = \frac{x}{10} + \frac{1}{1+x^2} \quad (2.14)
В точке оценки x_0 = 2.0.
На следующих рисунках показаны соответствующие аппроксимации Тейлора до значения n=3.
f(x)
n=0
n=1
n=2
n=3
n in[0, … ,3]
Производные
Для аппроксимации первой производной функции f(x) на x=x_0, т.е. f'(x_0), разложение Тейлора вплоть до n=2 могут быть использованы
f(x_0 + h) = f(x_0) + f'(x_0)h + \mathcal{O}(h^2)\quad (2.15)
После перегруппировки это приводит к прямой дифференциации первого порядка
f'(x_0) = \frac{f(x_0 + h) - f(x_0)}{h} + \mathcal{O}(h)\quad (2.16)
Таким образом, с помощью этой формулы становится возможным аппроксимировать значение производной на основе двух значений функции, приведенных здесь в x_0 and x_0 + h.
Графическое представление формулы (2.16) показано ниже. Оно демонстрирует формулу для h=1.
Рис. 2.12. Формула прямой разности первого порядка
Хотя в приведенной выше формуле использовалось значение при x_0 + h и, таким образом, называется формулой прямых различий, значение функции можно использовать при x_0 - h. Это приводит к формуле обратных различий
f'(x_0) = \frac{f(x_0) - f(x_0 - h)}{h} + \mathcal{O}(h)\quad (2.17)
Рис. 2.13. Формула обратной разности первого порядка
Приближений первого порядка в большинстве случаев недостаточно, и необходимы приближения более высокого порядка. Используя большее количество значений функции, т. е. комбинацию нескольких рядов Тейлора, можно вывести формулу более высокого порядка, например, для формулы центральных различий второго порядка.
В качестве отправной точки используются следующие два ряда Тейлора
f(x_0+h) = f(x_0) + f'(x_0)\cdot h + \frac{1}{2}f''(x_0)\cdot h^2 + \mathcal{O}(h^3)
f(x_0-h) = f(x_0) - f'(x_0)\cdot h + \frac{1}{2}f''(x_0)\cdot h^2 + \mathcal{O}(h^3) \quad (2.18)
Вычитание первого уравнения из второго
f(x_0+h) - f(x_0-h) = 2f'(x_0)\cdot h + \mathcal{O}(h^3) \quad (2.19)
результаты в формуле центральной разности для первой производной
f'(x_0) = \frac{f(x_0+h) - f(x_0-h)}{2h} + \mathcal{O}(h^2)\quad (2.20)
Формула центральных различий графически представлена на следующем рисунке.
Рис. 2.14. Формула центральной разности второго порядка
Все три представленных метода показаны на совмещенном графике:
Рис. 2.15. Сравнение всех приведенных выше разностных формул
Также возможно вывести формулы направления (прямого, обратного) второго порядка. Вычисление приводит к следующим формулам
f'(x_0) = \frac{-3f(x_0) + 4f(x_0+h) - f(x_0+2h)}{2h} + \mathcal{O}(h^2)
f'(x_0) = \frac{ 3f(x_0) - 4f(x_0-h) + f(x_0-2h)}{2h} + \mathcal{O}(h^2) \quad (2.21)
Задание. Выведите формулу дискретизации для второй производной
f''(x_0) = \frac{f(x_0-h) - 2f(x_0) + f(x_0+h)}{h^2} + \mathcal{O}(h^2) (2.22)
Дискретизация
В приведенном выше разделе производная функции была аппроксимирована в одной точке, здесь x_0. Для того, чтобы вычислить производную функции в заданном интервале x \in [x_b, x_e], общий подход заключается в дискретизации интервала. Это означает, что он представлен конечным числом точек. Учитывая n_x по мере увеличения количества представляющих точек интервал становится следующим конечным множеством:
x \in [x_0, \dots, x_i, \dots, x_{n_x-1}]
с
x_0 = x_b; \quad x_{n_x-1} = x_e
Если представляющие точки расположены на равном расстоянии друг от друга, то расстояние ce \Delta x является
\Delta x = \frac{x_e - x_b}{n_x-1}\quad
и значение x_i рассчитывается
x_i = x_b + i\cdot \Delta x \quad
В двух- и трехмерных случаях эта дискретизация может быть представлена в виде сетки, см. рис. 1.19. Таким образом \Delta x часто называется интервалом между сетками.
Общим индексированием для пространственно дискретизированных функций является
f(x_i) = f_i\quad
где индекс, указывающий положение на пространственной сетке, помещается в виде нижнего индекса.
Приведенные выше определения и следующие примеры предназначены для одномерной установки. В общем случае моделирование пожара решается в трехмерной области с соответствующими величинами для y- и z-направлений, т.е. \Delta y, \Delta z.
Пример дискретизации: приведен n_x = 20, функция в формуле (2.14) принимает следующее дискретное представление.
Рис. 2.16. Одномерная дискретизация
Дискретные производные
Учитывая формулы аппроксимации и дискретизированную функцию, становится возможным вычислить производную функции только с использованием ее значений. Формулы аппроксимации оцениваются с помощью h = \Delta x.
Для функции, показанной выше, и использующей центральные разности второго порядка (2.20), плюс производные по направлению (2.21) на границах, аппроксимированная производная показана на следующем рисунке.
Рис. 2.17. Аппроксимация первой производной
Задание. Создайте график, аналогичный Рис. 2.17 для функции
f(x) = x^4 + 3x^3 + 2 \quad
- Дискретизируйте функцию с помощью n_x = 20 между -2 и 4 и постройте график непрерывного и дискретного представления.
- Используйте эту дискретизацию для аппроксимации первой производной. Постройте точную и аппроксимированную производную.
- Проделайте то же самое для второй производной.
Ошибка аппроксимации
Если известна точная производная, то можно вычислить точное отклонение численного приближения от точных значений. Общим подходом к вычислению разницы между двумя наборами данных является L2-норма или среднеквадратичная ошибка (RMSE). Учитывая разницу значений функции \Delta f, здесь
\Delta f_i = f'_{numeric,i} - f'_{exact,i}
в тех же местах x_i с i \in [0, \dots, n_x-1], значение RMSE равно
RMSE(\Delta f) = \sqrt{\frac{1}{n_x}\sum_{i=0}^{n_x-1}\left(\Delta f_i\right)^2}\quad
Для примера функции из предыдущего раздела ошибка аппроксимации может быть вычислена для различных значений \Delta x, как показано на следующем рисунке. Как и ожидалось для схемы аппроксимации второго порядка, ошибка ведет себя аналогично квадратичной функции. Однако из-за конечного представления чисел с плавающей запятой сходимость прекращается, когда начинает преобладать ошибка округления.
Рис. 2.18. Погрешность аппроксимации в зависимости от расстояния между ячейками
Задание. Почему ошибка аппроксимации не уменьшается в дальнейшем для \Delta x \lesssim 5\cdot 10^{-6}? Примите во внимание уравнение (2.20) и ошибку округления для 64-разрядных чисел с плавающей запятой.
Интегрирование по времени
Траектории
Поскольку PDE, рассматриваемые в гидродинамике, часто содержат время как одно измерение, интеграция по времени является другим аспектом, помимо пространственных производных.
Аналогично пространственной дискретизации дискретизируется измерение времени. Таким образом, временной интервал [t_b=0, t_e] разделен на n_t временные этапы:
t_n = n\cdot\Delta t\quad with \quad \Delta t = \frac{t_e}{n_t-1} \quad and \quad n \in [0, n_t-1]\quad
Для того, чтобы отличать временной индекс от пространственного, он помещается как надстрочный индекс и не должен быть смешан со степенями:
u(t_n) = u^n\quad
Дана производная функции по времени, т.е. правая часть (rhs),
\partial_t u = \dot{u} = rhs(u, t)
и начальное значение для u на t=0, необходимо вычислить временную эволюцию, также называемую траекторией. Интегрирование по времени направлено на аппроксимацию следующего значения решения u^{n+1}.
Ниже пункт оценки rhs сокращен до
rhs(u(t_n), t_n) = rhs^n \quad
Рис. 2.19. Интегрирование по времени
Явный метод Эйлера
На основе, опять же, ряда Тейора, прямое приближение следующего значения решения основано на
\frac{u^{n+1} - u^n}{\Delta t} = rhs^n + \mathcal{O}(\Delta t)\quad
Это ведет непосредственно к
u^{n+1} = u^n + \Delta t \cdot rhs^n + \mathcal{O}(\Delta t)\quad (2.23)
Уравнение (2.23) является явным / прямым методом Эйлера. Графически это продемонстрировано на (2.20).
Рис. 2.20. Интегрирование по времени – явная схема Эйлера
Неявный метод Эйлера
Оценка производной в (2.23) при t_n то только одна из возможностей. Тем же способом, что и пространственные производные по направлению, оценка может быть выполнена на t_{n+1} и приводит к следующей итерационной формуле:
u^{n+1} = u^n + \Delta t \cdot rhs^{n+1} + \mathcal{O}(\Delta t)\quad (2.24)
Уравнение (2.24) является неявным / обратным методом Эйлера. Графически это продемонстрировано на рис. 2.21. Это неявный подход, в отличие от явной схемы Эйлера, который не позволяет напрямую вычислять следующий шаг итерации, основываясь только на решении текущих значений u^n. Однако, поскольку rhs может быть вычислен для любого значения u, соответствующее значение u^{n+1} можно найти, что они удовлетворяют уравнению (2.24).
Рис. 2.21. Интегрирование по времени – неявная схема Эйлера
Метод Эйлера-Хойна
Метод Эйлера-Хойна иллюстрирует схему интегрирования по времени второго порядка.
Сначала выполняется явное предсказание
u_p^{n+1} = u^n + \Delta t\cdot rhs^n\quad
Это промежуточное решение используется для скорректированной аппроксимации
u^{n+1} = u^n + \frac{1}{2} \Delta t\cdot \left(rhs^n + rhs(u_p^{n+1}) \right) + \mathcal{O}(\Delta t^2)\quad
Рис. 2.22. Интегрирование по времени – схема предсказателя-корректора Эйлера-Хойна
Методы Рунге-Кутты
Существует широкий спектр схем интегрирования по времени. Одним из классов из них являются методы Рунге-Кутты. Классической схемой является следующий метод четвертого порядка:
u^{(1)} = u^n + \frac{1}{2}\Delta t\cdot rhs^n
u^{(2)} = u^n + \frac{1}{2}\Delta t\cdot rhs\left(u^{(1)}\right)
u^{(3)} = u^n + \Delta t\cdot rhs\left(u^{(2)}\right)
u^{n+1} = u^n + \frac{1}{6}\Delta t\cdot \left (rhs^n + 2 rhs \left(u^{(1)}\right) + 2 rhs \left(u^{(2)}\right) + rhs \left(u^{(3)}\right)\right) + \mathcal{O}(\Delta t^4) \quad
Метод конечных разностей
Подход к решению
Рассмотрим простое уравнение диффузии:
\partial_t \phi = \lambda\partial_{xx} \phi
Форвард Эйлер: Эта схема может быть оценена непосредственно – явно.
\frac{\phi^{n+1}_i - \phi^n_i}{\Delta t} = \lambda\frac{\phi^n_{i-1} - 2\phi_i^n + \phi_{i+1}^n}{\Delta x^2}
\phi^{n+1}_i = \phi^n_i + \Delta t \lambda \frac{\phi^n_{i-1} - 2\phi_i^n + \phi_{i+1}^n}{\Delta x^2}
Рис. 2.23. Явная схема Эйлера
Обратный Эйлер: здесь система линейных уравнений должна быть решена для phi^{n+1}.
\frac{\phi^{n+1}_i - \phi^n_i}{\Delta t} = \lambda \frac{\phi^{n+1}_{i-1} - 2\phi_i^{n+1} + \phi_{i+1}^{n+1}}{\Delta x^2}
\phi^{n+1}_i - \frac{\Delta t}{\Delta x^2}\lambda\left( \phi^{n+1}_{i-1} - 2\phi_i^{n+1} + \phi_{i+1}^{n+1} \right) = \phi^n_i
A\phi^{n+1} = \phi^n
Рис. 2.24. Неявная схема Эйлера
Устойчивость
Распространение информации
Важным свойством временного интегратора является его стабильность. Учитывая решаемую PDE и временной интегратор, можно провести простой (линейный) анализ устойчивости, например анализ устойчивости по фон Нейману. Результатом является коэффициент роста, который указывает на скорость роста малых возмущений. Если этот коэффициент больше 1, то схема нестабильна, поскольку колебания будут бесконечно возрастать.
- Явные схемы, как правило, нестабильны или условно устойчивы, т.е. если выполняется условие для временного шага
- Неявные схемы, как правило, безусловно стабильны, однако они имеют тенденцию ослаблять решение
Задача адвекции с постоянной скоростью демонстрирует условие устойчивости.
Низкая скорость распространения
- Модельное уравнение
\partial_t \phi + \partial_x(v_0\phi) = 0
- Небольшое значение \Delta t
- Пройденное расстояние за \Delta t:
\delta x = v_0 \Delta t \lt \Delta x
- Информация перемещается меньше, чем \Delta x с шагом по времени step \Delta t и фиксируется соседней точкой сетки
- Информация о решении не теряется во время интегрирования по времени
Рис. 2.25. Медленное распространение информации
Высокая скорость распространения
- Большое значение \Delta t
- Пройденное расстояние за \Delta t:
\delta x = v_0 \Delta t \gt \Delta x
- Информация перемещается быстрее, чем \Delta x шагом по времени \Delta t и поэтому больше не могут быть зафиксированы
- Информация теряется во время интегрирования по времени
Рис. 2.26. Быстрое распространение информации
Условие Куранта-Фридрихса-Леви (CFL)
В общем случае существуют условия устойчивости для явных схем, которые соотносят максимальную скорость распространения информации v_{max} и скорость сетки v_g = \Delta x / \Delta t:
v_{max} \lt CFL\cdot v_{g} \quad
Учитывая сетчатое разрешение \Delta x вышеуказанное условие ограничивает временной шаг при максимальных скоростях:
\Delta t \lt CFL\cdot \frac{\Delta x}{v_{max}}
Примечания:
- Скорость потока может быть вычислена различными способамиys ( L_\infty, L_1, L_2 нормы) скорость диффузии \propto 1/\Delta x
- Существуют и другие ограничения на \Delta tограничение массовой плотности и объема.
- Значение числа CFL зависит от схемы интегрирования по времени.
Последствия состояния CFL в FDS:
Условие для \Delta t в FDS используется
0.8 \leq CFL = \Delta t \left(\frac{\|\vec{v}\|}{\Delta x} + |\nabla \cdot \vec{v}|\right) \leq 1.0
- Если число CFL растет выше верхнего предела, временной шаг устанавливается равным 90% от допустимого, если оно падает ниже нижнего предела, то его увеличивают на 10%.
- Как правило, временной шаг сокращается в 2 раза, когда расстояние между ячейками уменьшается в 2 раза, т. е. общая вычислительная нагрузка увеличивается в 16 раз.
Дополнительная заметка о том, как выполнялись вычисления в прошлом
Краткий раздел из презентации Ричарда Фейнмана “Лос-Аламос снизу”. Он описывает, как отдельные люди выполняли численные вычисления для каждого узла сетки и типа вычислений.