CC BY
58
Вычислительные технологии
Том 16, № 2, 2011
Об использовании градиентных итерационных методов при решении начально-краевых задач для трехмерной системы уравнений Навье — Стокса
Ю.Н. Захаров, К. С. Иванов Кемеровский государственный университет, Россия e-mail: zyn@kemsu.ru, topspin83@mail.ru
Численно решается нестационарная трехмерная система уравнений Навье — Стокса, описывающая пространственное течение вязкой однородной несжимаемой жидкости, записанная как в естественных переменных, так и в переменных вихрь — векторный потенциал. Для решения уравнений Пуассона на каждом шаге по времени используется метод минимальных невязок с групповой оптимизацией параметров. Приводятся результаты решения некоторых тестовых задач, полученные с применением различных формулировок исходных задач на дифференциальном уровне.
Ключевые слова: численное моделирование, вычислительная гидродинамика, трехмерные уравнения Навье — Стокса, вязкая однородная несжимаемая жидкость, пространственные течения.
Введение
В течение уже многих лет усилия математиков и механиков всего мира направлены на исследование задач о течении вязкой несжимаемой жидкости, описываемых уравнениями Навье — Стокса. Нелинейность этих уравнений, их пространственно эллиптический характер, наличие малого параметра при старших производных создают серьезные трудности как при аналитическом изучении, так и при численном моделировании. В прошлом столетии были достигнуты значительные успехи в изучении движения вязкой несжимаемой жидкости применительно к ламинарным режимам течений в двумерном случае. Доказаны важные теоремы существования и единственности [1], разработаны эффективные численные алгоритмы (см. обзор в [2]), получены точные количественные результаты, согласующиеся с физическим экспериментом.
В последние десятилетия в связи с бурным развитием информационных технологий огромный интерес исследователей вызвали методы численного моделирования пространственных задач. Появление суперЭВМ и параллельных вычислений позволило получить количественные результаты при решении практически важных задач [3-6]. Несмотря на достигнутые успехи, решение трехмерных задач требует значительных усилий. Это объясняется тем, что помимо трудностей в теоретических исследованиях существуют и алгоритмические проблемы. В частности, ни один из обычно применяемых на практике и включенных в современные вычислительные пакеты численных методов решения трехмерных систем уравнений Навье —Стокса не является в достаточной степени универсальным. Выбор подхода во многом зависит от условий конкретной задачи. При этом, как известно, опыт решения двумерных задач не всегда находит
аналогию в трехмерном случае [7]. Поэтому проблема построения эффективных алгоритмов и их применения для решения пространственных уравнений Навье —Стокса остается весьма актуальной,
В настоящее время имеется известно достаточно большое количество методов решения пространственных задач гидродинамики (см, обзоры в [8, 9]), отличающихся как исходной дифференциальной постановкой, так и численной реализацией, В большинстве случаев для моделирования трехмерных течений используются уравнения Навье—Стокса, записанные в физических переменных. Главным преимуществом такой формулировки являются точное задание граничных условий для компонент вектора скорости, известных из физической постановки задачи, и относительная простота численных алгоритмов. Основная трудность при численном решении системы в исходной постановке связана с расчетом поля давления и выполнением уравнения неразрывности. Первый значительный успех в преодолении отмеченных трудностей достигнут благодаря идее искусственной сжимаемости [10]. Были построены эффективные схемы расщепления для решения регуляризованной системы уравнений Навье —Стокса [11-15] (см, также обзор [16]), Подобный подход предлагается использовать главным образом для получения стационарного решения [17],
Другой способ определения давления заключается в применении оператора дивергенции к уравнению сохранения количества движения, в результате чего получается уравнение эллиптического типа для давления. Такой подход вначале был введен при разработке метода маркеров и ячеек [18], а в дальнейшем получил развитие в методах расщепления по физическим факторам [19, 20] и методах семейства SIMPLE [21], Существенные проблемы при этом подходе — задание недостающих краевых условий для давления и обеспечение на каждом шаге по времени соленопдальностп поля скорости. Вторую группу составляют методы решения трехмерных задач в переменных вихрь — векторный потенциал [22], Важнейшим преимуществом такой формулировки является то обстоятельство, что на каждом шаге по времени уравнение неразрывности выполняется автоматически как на дифференциальном, так и на дискретном уровне. Однако у данного подхода существуют серьезные недостатки. Во-первых, количество дифференциальных уравнений в новой системе возрастает до шести. Во-вторых, существенной проблемой является задание граничных условий для вихря на твердых стенках. В-третьих, в пространственном случае в отличие от плоского возникает значительная трудность постановки краевых условий для векторного потенциала.
Из многих известных способов дискретизации дифференциальной модели мы будем использовать метод конечных разностей как универсальный инструмент численного моделирования, Однако при любом выборе способа дискретизации исходной дифференциальной задачи наряду с преодолением указанных выше методологических трудностей неизбежно возникает проблема построения эффективных методов решения систем алгебраических уравнений большой размерности, к которым сводится дискретная модель. Эта проблема становится особенно актуальной в нестационарном случае. При решении нестационарных пространственных задач как в естественных переменных, так и в переменных вихрь — векторный потенциал в большинстве случаев с целью поэтапного определения неизвестных функций используется многоступенчатый алгоритм [19, 22], На первом этапе, решая системы уравнений количества движения или переноса вихря, находят компоненты векторов скорости или вихря. На втором этапе из решения уравнений Пуассона определяют давление или компоненты векторного и скалярного потенциалов. Для решения систем уравнений количества движения или переноса г,их-
ря в основном применяются устойчивые схемы расщепления, и поэтому данный этап обычно не вызывает серьезных проблем. Основные трудности возникают при решении уравнений Пуассона, что связано, во-первых, с постановкой краевых условий и, во-вторых, с решением системы линейных алгебраических уравнений большой размерности с матрицей, во многих случаях не обладающей необходимыми для применения хорошо разработанных итерационных методов решения свойствами, В связи с этим целью настоящей работы является сравнение методов решения трехмерной нестационарной системы уравнений Навье — Стокса, записанной в различных формулировках, и реализация явной градиентной итерационной схемы с групповой оптимизацией параметров для решения разностных уравнений Пуассона,
1. Постановка задачи
Основные уравнения, описывающие пространственное течение однородной вязкой несжимаемой жидкости, имеют следующий вид (все величины обезразмерены):
dV 1
_ + (V.V)V=-VP+—ДУ, (1)
divV = 0, (2)
или в переменных вихрь — векторный потенциал
^ + (V • V) w - (w • V) V = ^Ао;, (3)
Дф = -w, (4)
где V — вектор скор ости, P — давление, Re — число Ре йнольдса, w = го t V — вектор вихря, ф — векторный потенции, определяемый соотношением V = rot ф. Будем счи-
тать, что течение жидкости происходит в некоторой односвязной области С с границей дС па временном промежутке [0, Т],
Непосредственное интегрирование системы уравнений (1), (2) затруднено, поэтому обычно при численном подходе вместо уравнения (2) используют уравнение для давления
АР = ёЬ'((У -У)У), (5)
которое получается путем применения оператора дивергенции к уравнению (1), При этом численный метод решения новой системы уравнений (1), (5) должен обеспечивать соленондальность поля вектора скорости.
Для системы (1), (5) ставятся начальные и краевые условия:
У|4=0 = Уо(ж), х е С, (6)
= VI(х,*), * е [0,Т], (7)
где V0(x) и VI (х,*) — известные вектор-функции своих аргументов. Граничные условия для давления в исходной постановке задачи отсутствуют, поэтому при численном решении их получают из уравнения (1) (см., например, [2]),
Для системы уравнений (3), (4) постановка задачи несколько осложняется формулировкой граничных условий для составляющих векторного потенциала и вихря. Граничные условия для вихря явно не присутствуют в дифференциальной постановке задачи,
в силу чего при численном решении их получают из соотношения ш = rot V с учетом информации о поле скорости па границе и во внутренних точках области G, Задание краевых условий для векторного потенциала также явно не вытекает из краевых условий для скоростей, хотя имеется ряд работ [23, 24], в которых приводится их общая постановка. Мы рассмотрим отдельно постановки краевых условий для течений с протеканием и без протекания через границу, В случае отсутствия протекания граничные условия для векторного потенциала определяются достаточно просто:
фп \эс = Фт2\ас = 0, (8)
где ti , т2 — касательные век торы, n — внешняя нормаль к 0G, фТ1 , фТ2 — касательные составляющие векторного потенциала. Из равенства div ф \ dG = 0 находится условие для нормальной составляющей векторного потенциала фп. Например, в прямоугольной области отсюда следует соотношение
fl^o. (0)
Для задач протекания будем применять разложение векторного поля скорости на потенциальную и вихревую составляющие:
V = rot^ + Vp, (10)
где ф — векторный, р — скалярный потенциалы. Тогда ф определяется из решения задач (3), (4), (8), (9), а для р в силу div V = 0 получается задача
Ар = 0, (11)
= Vn\dG, (12)
где Vn — нормальная составляющая вектора скорости, В случае использования скалярного потенциала возникает необходимость в решении дополнительной задачи Неймана для уравнения Пуассона на каждом физическом временном шаге.
2. Метод решения
Введем в области ^ с границей сетку Gh■ Например, для областей,
состоящих из прямоугольных параллелепипедов, сетка будет иметь следующий вид:
Сн = {(гДж, кАг), г = О^М, ] = 0^1, к = (Щ}. (13)
Для численного решения нестационарной системы уравнений Навье — Стокса, записанной в "естественных" переменных, будем использовать хорошо известный неявный метод расщепления по физическим процессам [25-28], Аппроксимируя задачу (1), (5), (6), (7) па сетке Gh разностной схемой, получим разностные задачи относительно дискретных значений компонент вектора скорости и давления. Для решения уравнения количества движения будем использовать схему стабилизирующей поправки [11], Таким образом, основные трудности при данном подходе связаны с решением уравнения Пуассона для давления на каждом шаге по времени.
Как было отмечено ранее, основным преимуществом формулировки системы уравнений Навье — Стокса в переменных вихрь — векторный потенциал является явное удовлетворение уравнению неразрывности на каждом шаге по времени. Одно из положительных свойств численных моделей при данном подходе — выполнение уравнения неразрывности и на разностном уровне. Однако в этом случае необходимо использовать специальные разнесенные сетки [22], Рассмотрим построение таких сеток на примере, когда область С представляет собой прямоугольный параллелепипед:
С = {(х, у, г) : 0 < х < X, 0 < у < У, 0 < г < Z}. Наряду с основной сеткой СН введем следующие сетки:
(14)
СН
С2
С3
гч4 СН
С 5 Н
С 6 Н
СН
3
1
1
гАж, ( 2 + - 1 Ау, ( к + - ) Ах ) : г = О, Ж, ] = О, Ь - 1, к = О, М - 1
г + - 1 Ах,]Ау, ( к + - ] Ах ] : г = О, N - 1, ] = О, Ь, к = О, М - 1
г + -у Ах, ^ + -у Ау, /сД^ : г = О, N - 1, ] = О, Ь - 1, к = О, М ^г + 0 Ах,]Ау, кА^ : г = 0,ЛГ- 1, ^ = О^Х, Л; = 0,М
гАж, ( ^ + - ) Ау, /сАг ) : г = О, Ж, ] = О, Ь - 1, к = О, М
1Ах,]Ау, ( к + - ) Ах ) : г = О, Ж, ^ = О, Ь, к = О, М - 1
г + -1Аж, Н + -]Ау, (/с + -)Аг) : г = 0,^-1,
0,£ - 1, к = 0,М - 1 !>.
Дивергенцию скорости и скалярный потенциал определим в точках сетки С7Н, компоненты скорости — на сетках СН, СН, С'Н, компоненты вихря и векторного потенциала —
в точках сеток СН, СН, СН и, наконец, дивергенцию вихря — в точках и сходной сетки С7, ННН
ности автоматически выполняется на каждом дискретном временном шаге в центрах ячеек расчетной области. Аппроксимируя задачи (3), (4), (8), (9), (10), (11), (12) на построенных сетках разностной схемой, получим разностные задачи относительно дискретных значений компонент вектора вихря, векторного и скалярного потенциалов. Для решения уравнения переноса вихря, как и для решения уравнения количества движения, мы использовали схему стабилизирующей поправки [11], Это позволяет вновь сосредоточить основные усилия на решении уравнений Пуассона для векторного и скалярного потенциалов.
На каждом дискретном временном шаге разностные задачи для давления, векторного и скалярного потенциалов можно записать как системы линейных алгебраических уравнений вида
Аи = /, (15)
где и, / е Ят, А — линейный оператор: Кт ^ Кт [29, 30], При решении разностных задач для уравнений (5), (4), (11) в непрямоугольных областях и на неравномерных сетках
от оператора вряд ли можно ожидать таких "хороших" свойств как самосопряженность, знакоопределенность и неособенность. Данный факт затрудняет выбор эффективного метода решения системы линейных алгебраических уравнений (15) [29, 30], Именно поэтому в некоторых случаях для решения данной системы рекомендуется использовать прямые методы, что однако накладывает серьезные ограничения на шаг сеток по пространственным координатам. По существу, применение прямых методов ограничивается расчетами, цель которых — выявление лишь качественной картины течения. Понятно, что для получения адекватных количественных показателей система (15) должна иметь достаточно большую размерность, и в этом случае необходимо использовать итерационные методы решения систем линейных алгебраических уравнений, которые, с одной стороны, были бы эффективными с точки зрения сходимости итерационных процессов, а с другой — использовали бы минимальные сведения об операторах этих систем, В настоящей работе для решения разностных задач, аппроксимирующих уравнение Пуассона, предлагается применять итерационный метод, впервые использованный в [31] и затем развитый в [32], который показал высокую эффективность при решении стационарных уравнений Навье — Стокса, Для решения системы (15), следуя [32], построим следующий итерационный процесс:
и
п+1/2
и
- Тп+1 [Аип - /],
и
п+1
и
п+1
ип+1/2,
ип+1 _ ик-1
Е
г€Як
а
п+1
7п к — 1 2 1
гг , к ^ , ^ , ••• , 1 ,
и
п+1
и
п+1
(16)
где и0 — произвольный век тор из Ят, тп+1; аП+1 — итерационные параметры, Бк — непересекающиеся подмножества (группы) множества всех индексов такие, что и Бк — Б,
к
г™ — векторы из Ят с одной ненулевой г-й компонентой. Оптимальные итерационные параметры тп+1 и аП+1 в (16) определяются из условия минимума норм векторов невязок гп+1/2 и ГП+1 соответственно. Как показано в [32], при таком выборе итерационный процесс (16) является сходящимся при любом начальном приближении и произвольном неособенном операторе А Для определения оптимальных значений а^ +1 при каждом к — 1, 2,..., I необходимо решить систему линейных алгебраических уравнений вида
(
(А^ ,Агп) (А^ ,А*;
Агп ,
кРк
Агпх
Агп ,
кРк
к2
А гп Агк 2
Агп Агп 1 Агк1, Агк
Агп Агп Агк2, Агк
Агп ,
< кРк
Агп
\
Рк,
х
( ак+1 \
а
к1 +1 к2 +1
V а/+к1
( Агп , гп+1 \
к1к
Аг
к2 , гк
+1
Агп
Рк
, г +1 , гк
(17)
где кг е Бк, г
1, 2,...,Рк,
^ рк — т. Отметим, что чем большее количество индексов
к=1
кг участвует в группах Бк (т.е. чем болыне рк или меньше /), тем большую эффективность алгоритма (16) следует ожидать. Вместе с тем при увеличении размерности Бк
по степени сложности эквивалентно решению исходной системы (15), Однако оператор системы (15) не является произвольным, а есть некоторая аппроксимация дифференциального оператора Лапласа, и, следовательно, его матрица имеет блочно-ленточную
структуру. Например, в простейшем случае, когда область расчета — параллелепипед, матрица оператора имеет только семь ненулевых диагонадей. Очевидно, что в силу такой структуры матрицы оператора многие из величин (АгП,Агп) , г,3 = 1, 2,...,т, равны нулю. Следовательно, всегда возможно выбрать группы Бк так, что для каждого к матрица системы (17) будет иметь диагональный вид и нахождение агп +1 не составит труда [32], Заметим, что при увеличении размерности исходной системы (15) и неизменности структуры ее матрицы количество групп Бк, обеспечивающих диагональный вид матрицы системы (17), не возрастает (/ = сош^, то увеличиваются значения рк — количество индексов, входящих в каждую из групп. Таким образом, в этом смысле эффективность метода при увеличении размерности исходной системы не понижается. Рассмотрим отдельно ситуацию, когда имеется только одна группа индексов, совпа-
Б
сойдется к точному решению системы (15) за одну итерацию. Для определения оптимальных итерационных параметров аП+1, г = 1, 2,...,т, необходимо решить систему (17), размерность которой здесь совпадает с размерностью исходной системы, В данной постановке задача эквивалентна исходной. Но в качестве гП можно брать не единичные векторы, а такие, что (А^П, Агп) =0, г,3 = 1, 2,..., т, г = 3, при этом матрица системы (17) вновь будет иметь диагональный вид. Иначе говоря, необходимо построить ортогональную систему векторов гп, 3 = 1, 2,..., т, в пространстве со скалярным произведением (А*Аг, г) , г е Ят. Как известно, задача ортогонализации системы векторов в общем случае эквивалентна по степени сложности задаче решения системы линейных уравнений. Однако такой подход в нашем случае обладает рядом значительных преимуществ: во-первых, в силу упомянутой выше блочно-ленточной структуры матрицы исходного оператора А часть векторов гП уже ортогональны в во-вторых, если однажды построить указанную ортогональную систему векторов, то ее можно применять и в последующем для решения системы (15) при условии, что оператор А не изменился; в-третьих, можно решить задачу построения векторов, обеспечивающих не диагональный, а, например, трехдиагональный (не представляющий труда для обращения) вид матрицы системы (17), что значительно экономит вычислительные затраты. Отметим, что разностный оператор систем уравнений для давления, векторного и скалярного потенциалов не изменяется при переходе от слоя к слою по времени, а также сохраняет структуру своей матрицы при увеличении количества ячеек расчетной области, Этот факт позволяет построить достаточно эффективный алгоритм (16) с помощью оптимального выбора групп Бк, включив в каждую из них максимально возможное число индексов. При проведении тестовых расчетов очень экономичным оказывается применение ортогональной системы векторов гп, 3 = 1, 2,...,т, обеспечивающей сходимость процесса (16) за одну итерацию к точному решению системы (15), Построив такую систему ортогональных векторов на первом шаге по времени и затратив на это существенное количество ресурсов, ее можно использовать не только для быстрого решения системы уравнений (15) на последующих временных слоях, но и для проведения новых расчетов, например, при других числах Рейпольдса или краевых условиях.
3. Результаты расчетов
Для проверки эффективности предложенного метода решения диекретизированных систем уравнений (5), (4), (11) были проведены численные расчеты классической мо-
дельной задачи о течении вязкой однородной несжимаемой жидкости в кубической каверне с неравномерно движущейся верхней крышкой и задачи о течении вязкой однородной несжимаемой жидкости в трехмерном канале. Эти задачи довольно хорошо исследованы многими авторами [33-41] как в плоском, так и в пространственном случае, Однако в пространственном случае существенно проявляется трехмерный характер течения, что принципиально отличает его от плоского аналога. Поэтому в сочетании с наличием лабораторных экспериментов данные задачи являются содержательным тестом для апробирования и проверки работоспособности численных алгоритмов. Несмотря на заранее известный симметричный характер течений, для более надежной оценки правильности полученных результатов проводились расчеты полных задач без использования условий симметрии. Первая серия расчетов была проведена для исследования внутренней задачи. Область С в этом случае представляет собой куб (каверну) с длиной ребра Ь = 1. Верхняя крышка каверны движется вдоль оси X в сторону возрастания значений координаты со скоростью 11(1) = 1 — ^ ^. Остальные стенки каверны
неподвижны. Вся граница дС является непроницаемой. Для оценки точности и эффективности алгоритма расчеты проводились на последовательности сеток с числом узлов 21 х 21 х 21 41 х 41 х 41 65 х 65 х 65 при шаге по времени т = 0.01 и умеренных числах Рейнольдса Ее = 100 и 250, При шаге сетки к > 0.05 для решения уравнений Пуассона использовался алгоритм (16) с полной оптимизацией по всем параметрам одновременно с помощью построения ортогональной системы векторов. При более мелкой сетке применялся итерационный процесс (16) с групповой оптимизацией параметров. При этом группы Бк ,к = 1, 2,..., 6, были выбраны так, чтобы матрицы систем (17) имели диагональный вид,
= 100, к =
0.015625 в трех центральных сечениях каверны: у = 0.5 (а), х = 0.5 (б), г = 0.5 (в). Как и предполагалось, сечение у = 0.5 представляет собой плоскость симметрии, что подтверждается расчетами (см, рис, 1,6, в). Имеет место хорошее качественное совпадение полученных результатов с данными, приведенными в [37], На рис, 2 представлены распределения вертикальной компоненты скорости по поперечной (рис, 2, I) и продольной (рис, 2, II) координате при Ее = 100 в начальный момент времени (о), в развивающейся (б) и установившейся (в) стадиях. Линии 1 соответствует решение при к = 0.05, линии 2 — при к = 0.025, линии 3 — при к = 0.015625, Наблюдается сходимость реше-
Рис. 1. Линии тока установившегося
течения
0.03 0.028 0.026 0.024 0.022 0.02 0.018 0 016 0 014 0.012 0.01 0.008 0.006 0.004 0.002
Рис. 2. Распределения вертикальной компоненты скорости по поперечной (I) и продольной (II) координатам
ния разностной задачи, а также хорошее количественное совпадение графиков функций распределения компонент вектора скорости (с точностью до 10-5) с аналогичными результатами, приведенными в [34].
На рис. 3 приведен процесс развития течения на начальной, промежуточной и конечной стадиях при к = 0.025, Ее = 200. В начальной момент времени (рис. 3, о) слои жидкости движутся параллельно друг другу аналогично плоскому случаю. На стадии развития течения (рис. 3, б) слои жидкости постепенно переходят в соседние положения, образуя при этом торовидные вихревые зоны. Наконец, на стадии установления течения (рис. 3, в) наблюдаются существенные перемещения частиц жидкости вдоль третьей координаты. Частицы со стенок каверны закручиваются, перемещаясь в центр, затем раскручиваются и продолжают движение. На рис. 3, г для более наглядного анализа симметрии движения жидкости показан вид спереди полностью установившегося течения.
1
г 0.8
0.6
0.4
0.2
о
х То аз
0.6 « 1
Рис. 3. Процесс развития течения на начальной (а), промежуточной (б) и конечной (е. г) стадиях при Н = 0.025, Ые = 200
Вторая серия расчетов проводилась для исследования задачи протекания. Область О в этом случае представляет собой куб (канал) с длиной стороны Ь = 1, Жидкость в канале движется вдоль оси X в сторону возрастания значений координаты. Входная граница канала имеет верхнее отверстие высотой Н = 0.5, в которое горизонтально
втекает жидкость со скоростью и(Ь)
1 -
1
Все остальные стенки канала, кроме
1 + г5'
выходной границы, неподвижны и непроницаемы. На выходной границе задано условие свободного вытекания жидкости, определяемое постановкой мягких краевых условий па компоненты векторов скорости и вихря. Точность и эффективность алгоритма исследовались при тех же условиях, что и в предыдущей серии расчетов. На рис. 4 приведен процесс развития течения па начальном, промежуточном и конечном этапах при к = 0.025, Не = 100.
На начальном этапе (рис. 4, а, 6) слои жидкости движутся параллельно друг другу аналогично плоскому случаю. На этапе развития течения (рис. 4, в, г) слои жидкости по-
1
0.8 0.6 0.4 0.2
0 0.2 0.4 0.6 0.3 1 %
1
0.3 0.6 0.4 0.2
0 0.2 0.4 0.6 0.8 1 У
X
Рис. 4. Процесс развития течения на начальной (а. б), промежуточной (е. г) и конечной (д, е) стадиях при Н = 0.025 К-е = 100
степенно переходят в соседние положения, образуя в центре входного отверстия канала торовидную вихревую зону. На этапе установления течения (рис, 4, д, е) наблюдаются существенные перемещения частиц жидкости вдоль третьей координаты. Частицы со стенок канала закручиваются, перемещаясь к центру входного отверстия канала, затем раскручиваются и продолжают движение по направлению к выходу из канала.
На основании проведенных численных экспериментов следует вывод, что для внутренних задач наиболее эффективным является применение исходной дифференциальной постановки в переменных вихрь — векторный потенциал, в то время как для задач протекания оптимально использование естественных переменных. Такое утверждение базируется на ряде соображений, В первом случае не возникает особых трудностей с постановкой граничных условий для векторного потенциала, и решение трех уравнений Пуассона со смешанными краевыми условиями с помощью изложенного алгоритма не вызывает проблем. При использовании естественных переменных всегда появляется необходимость в решении задачи Неймана для давления. При этом матрица системы линейных алгебраических уравнений, получающейся в результате дискретизации дифференциальной задачи, всегда оказывается вырожденной. Решение систем с произвольной вырожденной матрицей на каждом шаге по времени представляет собой достаточно трудную задачу, поэтому, по существу, применение без особых трудностей методов с естественной постановкой возможно лишь в самосопряженном случае, в то время как при использовании постановки вихрь — векторный потенциал для решения внутренних задач такого ограничения нет. При численном решении задач протекания с применением дифференциальной постановки вихрь — векторный потенциал появляются серьезные проблемы с заданием граничных значений векторного потенциала, в силу чего в основном прибегают к введению скалярного потенциала, однако здесь возникает дополнительная задача Неймана, Таким образом, уравнение Пуассона с краевыми условиями второго рода на всей границе в задачах протекания необходимо решать при использовании любой из дифференциальных формулировок, В случае же переменных вихрь — векторный потенциал на каждом шаге по времени в дополнение к задаче Неймана нужно решать еще три уравнения Пуассона со смешанными граничными условиями для компонент векторного потенциала. Поэтому решение задач протекания с использованием физических переменных представляется более целесообразным.
Заключение
В работе исследованы возможности применения различных формулировок исходных дифференциальных уравнений Навье — Стокса к моделированию нестационарных пространственных задач разного типа (внутренних, протекания). Реализован итерационный метод неполной аппроксимации с групповой оптимизацией параметров для решения систем линейных алгебраических уравнений, сходящийся независимо от свойств оператора системы. Метод показал достаточную эффективность при решении разностных аналогов уравнений Пуассона для давления, векторного и скалярного потенциалов. Реализован также алгоритм, основанный на оптимизации по всем итерационным параметрам одновременно, с помощью которого можно получить точное решение системы за одну итерацию. Такой подход весьма экономичен при проведении тестовых серийных расчетов нестационарных задач.
Список литературы
[1] Ладыженская O.A. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970. 288 с.
[2] Роуч П. Вычислительная гидродинамика: Пер. с англ. М.: Мир, 1980. 616 с.
[3] Гущин В.А., Матюшин П.В. Математическое моделирование пространственных течений несжимаемой жидкости // Матем. моделирование. 2006. Т. 18, № 5. С. 5-20.
[4] Гущин В.А., Матюшин П.В. Механизмы формирования вихрей в следе за сферой при 200 < Re < 380 // Механика жидкости и газа. 2006. Л*8 5. С. 135-151.
[51 Tychkov S.A., Chernov V.V., Chernykh G.G. Numerical modeling of 3D convection in the Earth mantle // Russ. J. Numer. Anal. Math. Modeling. 2005. Vol. 20, No. 5. P 483-500.
[6] Червов B.B., Черных Г.Г., Червов A.B. Численное моделирование трехмерной конвекции под кратонами Центральной Азии // Вычисл. технологии. 2009. Т. 14, № 5. С. 114-121.
[7] Балаганский М.Ю., Захаров Ю.Н. Численное сравнение дву- и трехмерных течений вязкой несжимаемой жидкости // Там же. 2006. Т. 11. Спецвыпуск, посвященный 85-летию со дня рождения H.H. Яненко. Ч. 1. С. 38-45.
[8] Пасконов В.М., Полежаев В.П., Чудов Л.А. Численное моделирование процессов тепло- и массообмена. М.: Наука, 1984. 288 с.
[9] Prudhomme S. On the Flow Stability of Finite Elements Approximation of the Xavier Stocks Equations. TICAM REPORT, 1998. 52 p.
[10] Владимирова H.H., Кузнецов Б.Г., Яненко H.H. Численный расчет симметричного обтекания пластинки плоским потоком вязкой несжимаемой жидкости // Некоторые вопросы вычислительной и прикладной математики. Новосибирск, 1966. С. 186-192.
[11] Яненко H.H. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука. Сиб. отд-ние, 1967.
[12] Темам R. Sur l'approximation de la solution des equations des Xavier Stocks // Bull. Soc. Matem. des France. 1968. Vol. 96. P. 115-152.
[13] Горовая E.H. О решении пространственных задач для уравнений Навье — Стокса по устойчивым разностным схемам на ЭВМ // Труды IV Всесоюзного семинара по численным методам механики вязкой жидкости. Новосибирск: ВЦ СО РАН, 1973. С. 95-100.
[14] Кузнецов Б.Г., См агулов Ш. Об аппроксимации уравнений Навье^Стокса // Численные методы механики сплошных сред. 1975. Т. 6, № 2. С. 70-79.
[15] Кузнецов Б.Г., См агулов Ш. О сходящихся схемах дробных шагов для трехмерных уравнений Навье^Стокса // Там же. 1984. Т. 15, № 2. С. 69-80.
[16] Ветлуцкий В.Н., Колобов Б.П., Кузнецов Б.Г., Черных Г.Г. Численные методы в динамике вязкой жидкости // Моделирование в механике. 1987. Т. 1(18), № 4. С. 22-45.
[17] Численное моделирование течений в турбомашинах / С.Г. Черный, Д.В. Чирков, В.Н. Лапин, В.А. Скороспелов, C.B. Шаров. Новосибирск: Наука, 2006. 202 с.
[18] Chorin A.J. Numerical solution of Navier — Stokes equations // Math. Comput. 1968. Vol. 22. P. 745-762.
[19] Белоцерковский О.M. Численное моделирование в механике сплошных сред. М.: Наука, 1984. 519 с.
[20] Гущин В.А. Развитие метода расщепления по физическим факторам для расчета течений несжимаемой жидкости // Численное моделирование в аэрогидродинамике. \!.. 1986. С. 90-97.
[21] Патанкар С. Численные методы решения задач теплообмена и динамики жидкости: Пер. с англ. М.: Энергоатомиздат, 1984. 152 с.
[22] Белолипецкий В.М., Костюк В.Ю., Шокин Ю.И. Математическое моделирование течений стратифицированной жидкости. Новосибирск: Наука, 1991. 170 с.
[23] hlrasakl G.J., Hellums J.D. A general formulation of the boundary conditions of the vector potential in three-dimensional hydrodynamics // Quart. Appl. Math. 1968. Vol. 26, No. 3. P. 331-342.
[24] HlRASAKl G.J., Hellums J.D. Boundary conditions on the vector and scalar potential in viscous three-dimensional hydrodynamics // Ibid. 1970. Vol. 28, No. 2. P. 293-296.
[25] Рыков B.B. Численное моделирование нестационарных течений несжимаемой вязкой жидкости // Журн. вычисл. математики и мат. физики. 1985. Т. 25, № 5. С. 789-793.
[26] Ковеня В.Н. Модификация метода расщепления для численного решения уравнений Эйлера и Навье — Стокса // XVI Международная школа-семинар по численным методам механики вязкой жидкости. Новосибирск, 1998.
[27] Волков К.Н. Реализация схемы расщепления на разнесенной сетке для расчета нестационарных течений вязкой несжимаемой жидкости // Вычисл. методы и программирование. 2005. Т. 6, № 2. С. 146-159.
[28] Воеводин А.Ф., Гончарова О. Н. Реализация метода расщепления по физическим процессам для численного решения трехмерных задач конвекции // Вычисл. технологии. 2009. Т. 14, № 1. С. 21-33.
[29] Марчук Г.И. Методы вычислительной математики. М.: Наука, 1977. 535 с.
[30] Самарский A.A., Николаев Е.С Методы решения сеточных уравнений. М.: Наука, 1978. 592 с.
[31] Захаров Ю.Н. Многошаговые схемы с вариационной оптимизацией итерационных параметров. Новосибирск, 1980 (Препр. ИТПМ СО АН СССР; № 47).
[32] Захаров Ю.Н. Градиентные итерационные методы решения задач гидродинамики. Новосибирск: Наука, 2004. 238 с.
[33] Белолипецкий В.М., Костюк В.Ю. Численное исследование рециркуляционных течений в трехмерной каверне // Журн. прикл. механики и техн. физики. 1990. № 1. С. 100-104.
[34] Исаев С.А., Судаков А.Г., Лучко H.H. и др. Численное моделирование ламинарного циркуляционного течения в кубической каверне с подвижной гранью // ИФЖ. 2002. Т. 75, № 1. С. 49-53.
[35] Исаков А.Б. К численному решению задачи о движении вязкой несжимаемой жидкости в кубической каверне при Re = 1000 // Моделирование в механике. 1990. Т. 4, № 2. С. 64-76.
[36] Iwatsu R., Kawamura T., Kuwahara К., Min Hyun J. Numerical simulation of three-dimensional flow structure in a driven cavity // Fluid Dynamics Res. 1989. Vol. 5, No. 3. P. 173-189.
[37] Елизарова Т.Г., Милюкова О.Ю. Численное моделирование течения вязкой несжимаемой жидкости в кубической каверне // Журн. прикл. механики и техн. физики. 2003. Т. 43, № 3. С. 453-466.
[38] Пановко М.Я. Численное моделирование пространственных течений вязкой несжимаемой жидкости в канале с уступом // Теплофизика высоких температур. 1989. Т. 27, № 6. С. 1126-1131.
[39] Белолипецкий В.М., Костюк В.Ю. Численное решение задачи протекания для системы уравнений неоднородной жидкости // Численные методы механики сплошных сред. 1986. Т. 17, № 2. С. 3-9.
[40] Грабарник С.Я., Цепов Д.С. Численный метод расчета вязкого течения в трехмерном канале произвольной формы // Матем. моделирование. 1998. Т. 10, № 10. С. 103-111.
[41] Данилов Ю.М., Ильина И.М., Сидтикова И.П., Бергман А.Н. Решение трехмерных течений вязких несжимаемых жидкостей в каналах прямоугольной формы // Естеств. и техн. науки. 2003. № 3. С. 88-97.
Поступила в редакцию 21 июня 2010 г., с доработки — 22 октября 2010 г.
