1.2 Численные методы решения задач вычислительной гидродинамики
1.2.1 Дискретизация определяющих уравнений
Выбор численной схемы дискретизации закона сохранения зависит от формы в которой этот закон представлен [17].
В качестве примера рассмотрим модельное уравнение переноса величины . Интегральная форма записи этого закона сохранения будет иметь вид
(20)
Применяя теорему о дивергенции, и устремляя объем области к нулю, можно перейти к дифференциальной форме записи
(21)
Домножая это выражение на произвольную функцию и интегрируя по частям, можно получить вариационную форму закона сохранения
(22)
Дискретизация дифференциальной формы записи осуществляется при помощи конечно-разностных методов. Интегральная форма записи дискретизируется при помощи метода конечных объемов, а вариационная форма --- при помощи метода конечных элементов.
2.1.1 Метод конечных разностей: FDM
Метод конечных разностей (FDM --- Finite Differences Method) является одним из самых первых и самых простых методов решения дифференциальных уравнений, особенно в случае его использования в задачах с простой геометрией.
Наиболее часто метод применяется на регулярных сетках, когда линии координатной сетки служат локальными координатными линиями. В методе конечных разностей исходное дифференциальное уравнение аппроксимируются системой линейных алгебраических уравнений, где неизвестными являются значения переменных решения в узлах сетки. При этом Каждое слагаемое исходного дифференциального уравнения представляется соответствующим конечно-разностным отношением. В результате получается система алгебраических уравнений относительно неизвестных узловых значений.
В принципе, метод конечных разностей может быть применен к любому типу сетки. На регулярных сетках метод конечных разностей очень прост и эффективен. Особенно просто в этом случае получить схемы более высокого порядка точности. Недостатком конечно-разностных методов является то, что законы сохранения не учитываются без специальной заботы об этом. Кроме того, имеет место ограничение по сложности геометрии расчетной области.
Метод конечных объемов: FVM
Метод конечных объемов (FVM --- Finite Volume Method) представляет собой главный способ решения связанных уравнений переноса импульса и турбулентности [22].
В отличие от метода конечных разностей, данный метод использует формулировку уравнений в интегральной форме. Расчетная область разбивается на определенное количество контрольных объемов (ячеек), каждому из которых сопоставляется неизвестная величина, представляющая собой среднее значение переменной по этому объему. Для того, чтобы получить алгебраическое уравнение, соответствующее интегральному, записанному для некоторой контрольной ячейки, необходимо необходимо осуществить два этапа аппроксимации
• Приближенные значения интегралов, входящих в уравнение, при помощи квадратурных формул, выражаются через значения подынтегральных выражений в точках границы.
• Значения переменных в точках границы ячейки интерполируются по их значениям, заданным в узловых точках.
Интегральное уравнение выполняется как для каждого отдельного контрольного объема в отдельности, так и для расчетной области в целом. Таким образом метод конечных объемов обладает свойством глобального сохранения, что является важным преимуществом этого метода.
Метод конечных объемов может применяться с любым типом сетки, так что он применим для сложных геометрий. Сетка определяет только границы контрольного объема и не нуждается в привязке к системе координат. По сравнению с методом конечных элементов, метод конечных объемов более приемлем для большинства программистов, менее сложен с математической точки зрения и требует меньшей памяти компьютера при том же числе расчетных узлов.
Метод обладает преимуществами несложного программирования, математической простоты и физической адекватности. Вследствие этих достоинств, большинство разработанных коммерческих программ численного решения задач гидродинамики используют метод конечных объемов.
Метод конечных элементов. FEM
Метод конечных элементов (FEM --- Finite Element Method) во многом подобен методу конечных объемов [20]. Область разбивается на раздельные объемы, или конечные элементы, которые в общем случае являются неструктурированными. В двумерном пространстве это обычно треугольники или четырехугольники, а в трехмерном наиболее часто используются тетраэдры или шестигранники.
В методе конечных элементов искомая функция аппроксимируется линейной комбинацией координатных функций. Для получения уравнений метода, исходные уравнения интегрируются по всей расчетной области с некоторым весом, в качестве весовых функций принимаются координатные функции. В самых простых методах используются функции формы, линейные в пределах каждого элемента, что гарантирует непрерывность решения на границах элементов.
Важным преимуществом методов конечных методов является их применимость для задач со сложными пространственными конфигурациями. Данные методы относительно просто анализировать математически.
1.2.2 Разбиение расчетной области на элементы
Для численной реализации того или иного метода необходимо сгенерировать сетку для дискретизации определяющих уравнений. Процесс построения сетки относится к ключевым моментам численного моделирования, так как рациональным выбором сетки можно значительно упростить решение уравнений [20].
Регулярные сетки
Традиционно, при решении задач вычислительной гидродинамики применялись регулярные сетки (четырехсторонние на поверхности и шестигранные в пространстве). Регулярность заключается в том, что сетка представляет собой упорядоченную по определенным правилам структуру данных с явно выраженными сеточными направлениями, которые, в общем случае, образуют криволинейную систему координат. В преобразованном (вычислительном) пространстве ячейки сетки являются топологическими прямоугольниками (двумерные задачи) или параллелепипедами (трехмерные задачи). Для дискретизации уравнений Навье-Стокса на регулярных сетках используют, как правило, метод конечных разностей.
При построении регулярных сеток в геометрически сложных областях применяют преобразование координат общего вида, связанных с поверхностью тела. Основная цель данного преобразования заключается в построении равномерной расчетной сетки в преобразованном пространстве таким образом, чтобы координатные линии совпали с границами физической области. Введение криволинейной системы координат общего вида повышает эффективность расчетов в вычислительной области, так как система координат является прямоугольной. Вместе с тем, при записи уравнений в криволинейных координатах, появляются дополнительные члены --- параметры преобразования, определяющие отображение физической области на пространство обобщенных координат.
Следует особо отметить случаи ортогональных и конформных сеток. В первом случае при дискретизации модельных уравнений обнуляются некоторые параметры преобразования --- компоненты метрического тензора преобразования (матрицы Якоби), находящиеся не на главной диагонали данного тензора. Следствием является уменьшение погрешности, и, следовательно, повышение точности решения. Использование конформных преобразований позволяет сохранить такую же структуру модельных уравнений, записанных в вычислительной системе координат, как и в декартовом пространстве.
Для построения регулярных сеток применяются следующие методы
• Алгебраические методы
• Дифференциальные методы
• Методы теории функций комплексного переменного
Из анализа литературы последних лет следует, что наиболее широко распространенными являются методы построения сеток на основе решения эллиптических уравнений и алгебраические методы --- трансфинитная интерполяция и метод многих поверхностей. Однако, ни один из вышеперечисленных подходов не является универсальным по отношению ко всему спектру задач вычислительной гидродинамики.
Блочные сетки
Для структурированных сеток сравнительно легко реализуются вычислительные алгоритмы на основе современных методов высокого порядка точности. Однако, диапазон геометрических объектов, описываемых структурированными сетками, ограничен. Как правило, невозможно построить единую сетку для всей расчетной области, в связи с чем производится разделение поля течения на подобласти, в каждой из которых генерируется своя сетка регулярной структуры. Блочный подход предоставляет широкие возможности для использования эффективных численных методов внутри отдельных блоков. Основной недостаток блочного подхода состоит в достаточно сложной процедуре сшивки решений, полученных в различных подобластях.
Можно выделить следующие подходы
• Многоблочные структуры. Физическая область разбивается на несколько зон или блоков, причем границы блоков могут не соответствовать границам физической области. Затем для каждого блока отдельно строится сетка в соответствии с граничными условиями для каждой подобласти. Сетки из разных блоков могут как стыковаться точно по поверхностям раздела физической области на зоны, так и пересекаться между собой.
• Иерархические блочные структуры. Данный метод подразумевает иерархическую вложенность блоков сетки друг в друга. Нижестоящие по иерархии сетки ''погружены'' в вышестоящие.
Неструктурированные сетки
Характерной особенностью неструктурированных сеток является произвольное расположение узлов сетки в физической области. Произвольность расположения узлов понимается в том смысле, что отсутствуют выраженные сеточные направления и нет структуры сетки, подобной регулярным сеткам. Число ячеек, содержащих каждый конкретный узел, может изменяться от узла к узлу. Узлы сетки объединяются в многоугольники (двумерный случай) или в многогранники (трехмерный случай). Как правило на плоскости используются треугольные и четырехугольные ячейки, а в пространстве --- тетраэдры и призмы. Основное преимущество неструктурированных сеток перед регулярными состоит в большей гибкости при дискретизации физической области сложной формы, а также в возможности полной автоматизации их построения. Для неструктурированных сеток легче реализуются локальные сгущения и адаптация сетки в зависимости от поведения решения. Современные программы генерации сеток позволяют за приемлемое время строить сетки для сколь угодно сложных геометрических объектов. Для дискретизации уравнений Навье-Стокса на неструктурированных сетках применяются методы конечных элементов и конечных объемов. Метод конечных разностей к таким сеткам неприменим.
Недостатком неструктурированных сеток по сравнению с регулярными сетками является необходимость хранить информацию о структуре сетки, что повышает требования к вычислительным ресурсам.
Тем не менее, данный подход получает все более широкое распространение по следующим причинам
• Процесс генерации сетки сравнительно легко формализуется, и автоматизируется. Генераторы неструктурированных сеток обладают универсальностью по отношению к широкому диапазону прикладных задач.
• В случае очень сложной геометрии расчетной области, время генерации неструктурированной сетки на порядок меньше времени генерации регулярной сетки.
• Измельчение произвольных участков сетки осуществляется естественным образом, что дает преимущество при использовании адаптивных и многосеточных методов.
Как показывает практика, выбор структуры для представления триангуляции оказывает существенное влияние на трудоемкость алгоритмов, использующих данную структуру, а также на скорость конкретной реализации [13]. Рассмотрим структуры данных, применяемые для двумерных треугольных сеток.
• ''Узлы с соседями''. Для каждого узла триангуляции хранятся его координаты на плоскости и список указателей на соседние узлы, с которыми есть общие ребра. Треугольники при этом не представляются вообще, что является обычно существенным препятствием для дальнейшего применения триангуляции.
• ''Узлы, ребра и треугольники''. В явном виде задаются все виды объектов триангуляции узлы, ребра и треугольники. Недостатком данной структуры является большой расход памяти.
• ''Узлы и треугольники''. Для каждого треугольника хранятся три указателя на образующие его узлы и три указателя на смежные треугольники. Данная структура триангуляции наиболее часто применяется на практике в силу своей компактности и относительного удобства в работе.
Гибридные сетки
Гибридные сетки комбинируют регулярные и неструктурированные области сетки. Данный подход позволяет сочетать достоинства и снизить влияние недостатков, присущих каждому типу сеток. Простейший пример гибридной сетки --- течение около системы профилей. Область около профилей покрывается ортогональной регулярной сеткой, а области между профилями и далекое поле --- неструктурированной. Вычислительный алгоритм содержит процедуру переключения вычислительных схем на различных сетках, и, при необходимости, переноса информации с одного типа сетки на другой.
1.2.3 Уравнения Навье-Стокса для несжимаемой жидкости
Для несжимаемой жидкости выполняются свойства
где -- число Маха, -- местная скорость звука. Предположение о несжимаемости является также хорошим приближением, например для воздуха при скорости , или числе Маха .
При рассмотрении задач о несжимаемой жидкости используется частный случай уравнений для сжимаемой жидкости. Эти уравнения, в случае отсутствия массовых сил и подвода тепла извне, записываются следующим образом
(23)
Эти уравнения образуют смешанную эллиптически-параболическую систему относительно неизвестных . Моделирование вихревых течений на основе численного решения уравнений Навье-Стокса для несжимаемой жидкости сопровождается рядом трудностей математического характера.
Уравнение неразрывности для несжимаемой жидкости содержит лишь составляющие скорости, в связи с чем нет прямой связи с давлением, которая для сжимаемых течений осуществляется через плотность. Вследствие этого одной из главных проблем является наложении условия несжимаемости, и возможное возникновение неустойчивости решения при наложении условия несжимаемости.
Другой проблемой при решении системы уравнений Навье-Стокса является нелинейность, связанная с конвективными слагаемыми в уравнениях движения, которая может приводить к появлению осцилляций решения в областях с большими градиентами. В случае превалирования конвекции над диффузией происходит ухудшение решения из-за жесткости, которую вносят в систему уравнений конвективные члены, и несимметричности матрицы системы линейных алгебраических уравнений. Наиболее эффективными методами учета конвективных слагаемых являются противопотоковые схемы, реализованные методом конечных объемов.
Метод завихренности и функции тока
В методе завихренности и функции тока [14], выполняется замена переменных для перехода от компонент скорости к завихренности и функции тока , которые в двумерных декартовых координатах определяются так
(24)
Используя новые независимые переменные, можно скомбинировать уравнения сохранения импульса, исключая из них давление, что дает
(25)
Это параболическое уравнение в частных производных называется уравнением переноса завихренности.
Зависимости (24) позволяют получить уравнение для и
(26)
Это выражение представляет собой уравнение Пуассона, являющееся эллиптическим.
Итак, в результате замены переменных, смешанная эллиптически-параболическая система уравнений разделилась на одно параболическое и одно эллиптическое уравнение. Обычно эти уравнения решают методом установления по времени.
Чтобы определить давление, необходимо решить уравнение Пуассона для давления
(27)
Распространение подхода с использованием завихренности и функции тока в качестве независимых переменных на трехмерные задачи осложнено тем, что для действительного трехмерного течения нельзя ввести функцию тока.
Проекционные методы и метод коррекции давления
Конечно-разностные и конечно-объемные методы решения делятся на методы, использующие процедуру коррекции давления (pressure-based algorithm), и методы, основанные на принципе расщепления неизвестных (pressure-velocity coupling) [26].
Метод расщепления по физическим факторам подразумевает разделение системы уравнений Навье-Стокса на последовательность более простых уравнений, таких как уравнения диффузии, переноса и уравнение Пуассона. Разработка численных методов для этих уравнений оказывается значительно проще, чем непосредственного для уравнений Навье-Стокса.
Фактически, производится расчет поля скоростей в два этапа. На первом этапе вычисляется промежуточное поле скоростей без учета уравнения неразрывности. На втором этапе осуществляется коррекция поля скоростей, чтобы обеспечить выполнение уравнения неразрывности. Тем самым осуществляется ''проекция'' поля скоростей на пространство векторов с нулевой дивергенцией, поэтому методы называются проекционными.
Схема расщепления по физическим факторам (метод проекции), применяемая для дискретизации уравнений Навье-Стокса, записанных в физических переменных [7].
Пусть в момент времени известны поле скорости и поле давления . Тогда для расчета неизвестных функций в момент времени используется схема расщепления [6], состоящая из трех этапов.
На этапе 1 предполагается, что перенос количества движения осуществляется только за счет конвекции и диффузии
(28)
Несмотря на то, что промежуточное поле скорости не удовлетворяет уравнению неразрывности, оно имеет физический смысл, так как сохраняет вихревые характеристики .
На этапе 2 по найденному промежуточному полю скорости , с учетом соленоидальности вектора скорости , рассчитывается поле давления
(29)
Для решения уравнения Пуассона (43) на каждом шаге по времени могут использоваться как итерационные, так и прямые методы.
На этапе 3 предполагается, что перенос количества движения осуществляется только за счет градиента давления (конвекция и диффузия отсутствуют)
(30)
Следует отметить, что уравнение Пуассона (29) получено путем взятия дивергенции от обеих частей равенства (30) с учетом уравнения неразрывности .
Метод искусственной сжимаемости
Рассматривается формулировка метода искусственной сжимаемости, в которой в уравнение неразрывности включено слагаемое, содержащее производную давления по искусственному времени [26, 22]. Полученное уравнение вместе с уравнениями сохранения количества движения образуют гиперболическую систему уравнений, для которой скорость распространения возмущений давления конечна.
Преимущество метода при решении задач стационарного течения, заключается в том, что он не требует получения поля скоростей с нулевой дивергенцией на каждой итерации. Дивергенция поля скоростей автоматически обращается в ноль, когда решение устанавливается во времени. Кроме того, данный метод может быть обобщен на случай нестационарного движения.
За счет введения искусственного уравнения состояния, уравнение неразрывности приводится к виду
(31)
Параметр , представляющий собой искусственную скорость звука должен быть выбран с учетом двух факторов. С одной стороны, искусственная скорость звука должна быть близка к скорости движения жидкости, обеспечивая лучшую обусловленность системы. С другой стороны, искусственная скорость звука должна быть достаточно велика, чтобы возмущения поля давлений успевали распространяться достаточно далеко до того, как произойдет их затухание. Тем самым осуществляется переход к установившемуся состоянию, при котором выполняется условие несжимаемости. Для обезразмеренных уравнений этим условиям удовлетворяет диапазон значений от до , предложенный в работе [29].
С точки зрения линейной алгебры, система линейных уравнений, возникающая при дискретизации уравнений метода искусственной сжимаемости, оказывается хорошо обусловленной (при правильном выборе параметра ) по сравнению с исходными уравнениями. Любые другие методы решения должны использовать специальные приемы для получения корректных распределений давления.
В случае применения метода искусственной сжимаемости для расчета нестационарных течений, применяется метод двойных шагов по времени (dual time step) [33].
В этом случае интегрирование по физическому времени осуществляется только для уравнения сохранения количества движения. На каждом шаге физического времени осуществляется итерационный процесс интегрирования по искусственному времени следующей системы
(32)
Таким образом, на каждом шаге физического времени обеспечивается выполнение условия несжимаемости.
1.2.4 Методы решения СЛАУ
Одним из наиболее емких, с точки зрения затрат машинных ресурсов, этапов вычислительной процедуры является решение системы линейных алгебраических уравнений (СЛАУ) большой размерности. Необходимость решения такой системы возникает в следующих двух случаях
• в результате применения неявных схем для дискретизации определяющих уравнений по времени,
• в результате дискретизации уравнения Пуассона при использовании проекционных методов наложения условия несжимаемости.
Рассмотрим систему линейных алгебраических уравнений, записанную в виде
где - матрица системы, имеющая размерность , вектор - вектор решения, - вектор правых частей.
Численные методы решения систем данного вида [3, 8] принято разделять на два класса: прямые методы и итерационные.
Методы, позволяющие получить решение системы уравнений за конечное число арифметических операций, называются прямыми. К ним относятся метод Крамера, метод исключения Гаусса, метод -разложения и ряд других методов. Основным недостатком прямых методов являются жесткие требования к быстродействию и памяти. Например, метод Гаусса требует выполнения порядка арифметических операций и хранения порядка переменных. Кроме того, если велико, то машинная погрешность вычислений будет оказывать значительное влияние на конечный результат.
Итерационный метод общего вида [2, 8] основан на последовательном улучшении начального приближения решения. Построение последовательности приближений осуществляется посредством единообразного процесса, называемого процессом итераций. Вычислительный процесс заканчивается, когда изменение решения при переходе к следующей итерации становится достаточно малым, или когда невязка уменьшается до заданного значения. Итерационные методы требуют хранения порядка переменных, а время решения зависит от обусловленности матрицы и качества начального приближения. Например для методов Якоби и Гаусса-Зейделя количество арифметических операций составляет порядка , где -- количество проведенных итераций. Особенностью итерационных методов является необходимость исследования сходимости каждого метода.
Общая структура итерационных методов связана с представлением матрицы в виде и видоизмененной формой исходного уравнения . Различные итерационные методы отличаются друг от друга способом выбора матрицы . Примером обычных итерационных методов могут служить метод Якоби (метод простых итераций), метод Зейделя, метод последовательной релаксации. К отдельному классу следует отнести вариационные итерационные методы: метод наискорейшего спуска, минимальных невязок, минимальных поправок, минимальных погрешностей и другие.
Классические итерационные методы
Рассмотрим для системы (53) итерационный процесс общего вида, заданный следующим образом
где -- невырожденная матрица. Перепишем выражение в виде
(33)
где -- оператор -го шага итерационного процесса, -- единичная матрица.
Пусть -- точное решение системы. Введем обозначения
где -- вектор ошибки, -- вектор невязки. В таком случае имеют место соотношения
Итерационный процесс (33) называется сходящимся, если последовательность сходится к при любом начальном приближении . Итерационный процесс называется стационарным, если не зависит от номера итерации . В противном случае процесс называется нестационарным.
Методы вариационного типа
В методах вариационного типа [2, 8] решение линейной алгебраической системы заменяется эквивалентной экстремальной задачей. Пользуясь обозначениями (53), образуем квадратичный функционал следующего вида
(34)
где символ обозначает скалярное произведение двух векторов. Можно показать, что если матрица симметрична и положительно определена, то задача решения системы линейных уравнений и задача минимизации функционала (55) эквивалентны. Для случая произвольной невырожденной матрицы , не являющейся симметричной, можно сформулировать аналогичную задачу минимизации функционала.
Конструирование итерационного процесса осуществляется путем применения к задаче (34) различных методов численной минимизации функционала.
Итерационные методы Крылова
[8] Наиболее эффективными и устойчивыми среди итерационных методов являются проекционные методы, связанные с проектированием на подпространства Крылова (Krylov subspace methods). По сравнению с классическими итерационными методами, они не содержат эмпирически подбираемых параметров и позволяют получить более высокую скорость сходимости, несмотря на увеличение числа операций на каждой итерации.
В семейство итерационных методов Крылова входят, в частности
• обобщенный метод минимальных невязок (Generalized Minimum Residual, GMRES),
• метод сопряженных градиентов (Conjugate Gradients, CG),
• метод квадратичных сопряженных градиентов (Conjugate Gradients Squared, CGS),
• метод бисопряженных градиентов (Bi-Conjugate Gradients, BiCG),
• метод бисопряженных градиентов со стабилизацией (BiCGStab).
В программном коде «SmartFlow» используются обобщенный метод минимальных невязок (GMRES) и метод бисопряженных градиентов (BiCG).
1.2.5 Выводы
Исходя из поставленной задачи расчета течений в сложных областях с подвижными границами, для разделения расчетной области на элементы наиболее перспективным, с точки зрения автора, оказывается неструктурированный тип сетки.
Для пространственной дискретизации определяющих уравнений выбран метод конечных объемов с высоким порядком аппроксимации, который, наряду с другими достоинствами, хорошо зарекомендовал себя при расчетах на неструктурированных сетках.
В качестве модели движения вязкой несжимаемой жидкости выбраны осредненные по Рейнольдсу уравнения Навье-Стокса. В данной работе моделирование турбулентности не рассмотренно, поэтому обезразмеренная турбулентная кинематическая вязкость принята равной 0.
Для наложения условия несжимаемости в уравнениях Навье-Стокса выбран подход с искусственной сжимаемостью. Результирующие уравнения дискретизируются по времени при помощи неявной схемы второго порядка точности. Возможность применения метода искусственной сжимаемости для решения нестационарных задач обеспечивается за счет применения метода двойных шагов по времени.
Учитывая, что система линейных алгебраических уравнений, возникающая в результате применения неявной схемы дискретизации по времени, имеет большую размерность и обладает разреженной несимметричной матрицей, для решения этой системы выбран итерационный метод бисопряженных градиентов и обобщенный метод минимальных невязок.
|
