автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Моделирование внутренних течений вязкой несжимаемой жидкости методом конечных элементов с использованием противопотоковых схем

На правах рукописи

Гобыш Альбина Владимировна

Моделирование внутренних течений вязкой несжимаемой жидкости методом конечных элементов с использованием противопотоковых схем

05 13 18 —математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРЕТ диссертации на соискание ученой степени кандидата физико-математических наук

00316155Т

Новосибирск 2007

003161557

Работа выполнена в Новосибирском государственном техническом университете

Научный руководитель доктор технических наук,

профессор Элла Петровна Шурина

Официальные оппоненты доктор физико-математических наук,

доцент Бердников Владимир Степанович,

доктор физико-математических наук Федорук Михаил Петрович

Ведущая организация Институт вычислительной математики и

математической геофизики СО РАН, г Новосибирск,

проспект Академика М А Лаврентьева, б

Защита состоится 007 г в 9 часов на заседании

диссертационного совета Д003.046 01 при Институте вычислительных технологий СО РАН по адресу 6300090, г Новосибирск, проспект Академика М А Лаврентьева, 6

С диссертацией можно ознакомиться в специализированном читальном зале вычислительной математики и информатики ГПНТБ СО РАН (проспект Академика М А Лаврентьева, 6)

Автореферат разослан" 0_» О&ШЛЯр^ 2007 г

Ученый секретарь диссертационного совета, А

доктор физико-математических наук, профессор/¡//и/^/^ Б Чубаров

Общая характеристика работы

Актуальность работы. Математическое моделирование процессов становится одним из наиболее распространенных методов исследования объектов и явлений различной природы. Широкий круг задач, стоящих перед современной наукой и техникой, связан с решением уравнений движения жидкости Возникает потребность в моделировании внутренних течений несжимаемой жидкости с преобладанием конвективного переноса в каналах сложной геометрии, являющихся элементами различных технических устройств

Движение вязкой несжимаемой жидкости описывается уравнениями Навье - Стокса При изучении движения несжимаемой жидкости необходимо наличие программных реализаций не только сложных моделей, но и более простых приближенных моделей, с помощью которых можно получить предварительное представление о характере течения В случае больших значений кинематической вязкости или при малых значениях вектора скорости решение уравнений Навье - Стокса может быть заменено решением уравнений Стокса Решения уравнений Эйлера, описывающих течение идеальной несжимаемой жидкости, могут служить начальным приближением для итерационных методов решения уравнений Навье - Стокса Поэтому актуальной задачей является разработка и программная реализация вычислительных алгоритмов, позволяющих единообразно решать уравнения данного класса задач

Для решения уравнений Навье - Стокса (Стокса и Эйлера), описывающих течения несжимаемой жидкости, существуют различные численные методы Одним из распространенных методов аппроксимации несжимаемых течений является метод конечных элементов (МКЭ) Важными требованиями к вычислительным схемам на базе метода конечных элементов являются эффективные способы учета (аппроксимации) уравнения неразрывности и линеаризации конвективных членов в уравнении движения Поэтому построение вычислительных алгоритмов с указанными свойствами является на сегодняшний день актуальной задачей

Цель работы состоит в моделировании внутренних течений несжимаемой жидкости с преобладанием конвективного переноса в каналах сложной геометрии, являющихся элементами различных технических устройств

Научная новизна работы.

• Разработана технология реализации противопотоковых схем на неструктурированных сетках с расчетом неизвестных в узлах конечных

элементов и серединах их ребер дня моделирования течений с преобладанием конвективного переноса В контексте метода конечных элементов предложен новый способ аппроксимации потоков через границы ячеек двойственной сетки

• Разработана и программно реализована вычислительная схема на основе смешанных постановок с введением функции штрафа, позволяющая единообразно моделировать стационарные внутренние течения вязкой несжимаемой жидкости (уравнения Навье — Стокса, Стокса) и идеальной несжимаемой жидкости (уравнения Эйлера) Предложен способ учета краевого условия u-n|r =h (n - вектор нормали к границе Г, h - заданная функция) в вариационных постановках уравнений Эйлера Разработана структура данных для хранения элементов матриц системы линейных алгебраических уравнений с расчетом неизвестных компонент вектора скорости и давления в гипервекторе

• Реализованы стабилизированные схемы Галеркина для моделирования трехмерных течений на тетраэдральном разбиении, в которых используется метод стабилизации давления и в качестве стабилизирующего оператора выбирается конвективная часть исходного оператора

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

Обоснованность и достоверность результатов, полученных в диссертации, обеспечивается корректностью постановок рассматриваемых задач и методов их решения, основывается на расчетах широко известных и рекомендуемых тестовых задач и сопоставлении результатов численных, расчетов с результатами, полученными другими авторами

Апробация работы. Созданный комплекс программ внедрен в расчетную практику в Сибирском центре научно-технического обеспечения AHO "Промбезопасность - Сибирь" Результаты решения задач с доминированием конвекции применяются в учебном процессе на факультете летательных аппаратов Новосибирского государственного технического университета в курсе лекций "Экологические проблемы энергетики" и "Численные методы в задачах экологии".

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

Новосибирск, 2002 г ), региональной научной конференции студентов, аспирантов и молодых ученых "Наука. Техника Инновации" (г Новосибирск, 2002 г ), международной конференции "Вычислительные и информационные технологии в науке, технике и образовании" (г Усть-Каменогорск, 2003 г ), IV всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям (г Красноярск, 2003 г ), всероссийской научной конференции студентов, аспирантов и молодых ученых "Наука Технологии Инновации" (г Новосибирск, 2003 г, 2004 г ); международной конференции по вычислительной математике МКВМ-2004 (г. Новосибирск, 2004 г ), V всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям с участием иностранных ученых (г Новосибирск, 2004 г ), семинарах Новосибирского государственного технического университета, Института вычислительных технологий СО РАН и Института вычислительной математики и математической геофизики СО РАН.

Публикации. По теме диссертации опубликовано 6 работ, в том числе (в скобках в числителе указан общий объем этого типа публикаций, в знаменателе - объем, принадлежащий лично автору) 1 статья в издании, рекомендованном ВАК для представления результатов докторских диссертаций (0 62/0 4 печ л ), 2 - в трудах международных конференций (0 68/0 5 печ. л), 3 - в сборниках научных трудов (1 37/1.2 печ. л )

Личный вклад автора. В публикации [1] автор участвовала в постановке задачи, осуществляла программную реализацию вычислительных алгоритмов и сравнение полученных результатов с конечно-разностным решением рассматриваемых задач В работе [2] диссертант осуществляла разработку алгоритмов для решения уравнений Навье — Стокса, основанных на методе конечных элементов В публикации [3] автору принадлежат конструирование и реализация алгоритмов решения задач конвективно-диффузионного переноса с использованием противопотоковых схем В публикации [6] автор участвовала в постановке задач, проведении расчетов и анализе результатов

Структура и объем диссертации. Диссертация состоит из введения, трех глав, заключения, списка использованных источников (131 наименование) и приложения Полный объем диссертации составляет 138 страниц, включая 23 таблицы и 40 иллюстраций

Автор выражает искреннюю признательность и глубокую благодарность научному руководителю д т н., профессору Элле Петровне Шуриной

икф-мн,снс ИВТ СО РАН Нине Юрьевне Шокиной за помощь и поддержку при работе над диссертацией

Основное содержание работы

Во введении обоснована актуальность темы, сформулированы постановка исследуемой проблемы, цели и методы исследований, представлены научная новизна и значимость работы Анализируются современные подходы в моделировании течений несжимаемой жидкости, их преимущества и недостатки Кратко описаны структура и основное содержание диссертации по главам

В диссертационной работе исследованы стационарные, внутренние, однофазные, химически однородные течения Уравнения Навье — Стокса описывают движение вязкой несжимаемой жидкости Для полноты исследования явлений в рассмотрение включены уравнения Стокса и Эйлера В случае больших значений кинематической вязкости или при малых значениях вектора скорости решение уравнений Навье - Стокса может быть заменено решением уравнений Стокса, а решения уравнений Эйлера, описывающие течение идеальной жидкости, могут служить начальным приближением для итерационных методов решения уравнений Навье - Стокса

Существуют несколько формулировок уравнений движения жидкости в естественных переменных вектор скорости - давление, в переменных векторный потенциал - вектор вихря, в переменных вектор скорости -вектор вихря Достоинства и недостатки перечисленных формулировок приведены в работах С Патанкара, Р Темама, К Флетчера, F Brezzi, M Fortin, ZUA Warsi и др Отмечено, что при использовании моделей в естественных переменных вектор скорости - давление на уровне вычислительной схемы возникают трудности корректного учета условия, накладываемого на дивергенцию аппроксимируемого поля скоростей, так называемого дивергентного условия

На сегодняшний день ведущие позиции при численном решении дифференциальных уравнений занимают сеточные методы Наиболее популярными для решения уравнений Навье — Стокса, Стокса и Эйлера являются методы конечных разностей, конечных объемов и конечных элементов Численный анализ вычислительных схем, основанных на МКЭ, представлен в работах Р Темама, К Флетчера, F Brezzi, M Fortin, T J R Hughes, D L Coulliette, M Koch, L J P Timmermans, T Tezduyar и др

Несмотря на ряд достоинств МКЭ при решении задач гидродинамики (например, возможность конструирования конечно-элементных аппрокси-

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

В диссертационной работе представлен обзор схем аппроксимации нелинейности, связанной с конвективными членами в уравнении движения При доминировании конвекции используются противопотоковые схемы, стабилизированные и разрывные схемы Галеркина Среди представленных в научной литературе противопотоковых схем выделены схемы с расчетом неизвестных в узлах конечного элемента, серединах его ребер и в центрах конечного элемента

Несмотря на то, что диссертационная работа посвящена стационарным течениям, проведено исследование разработанных вычислительных схем на нестационарных задачах В работе представлен обзор схем аппроксимации по времени, используемых при численном решении дифференциальных уравнений явных, неявных и полунеявных схем Рассмотрены способы диагонализации матрицы массы в схемах аппроксимации по времени Преобразование элементов матрицы с целью ее диагонализации применяется для обращения матрицы массы в вычислительных схемах расчета компонент вектора скорости и давления

В диссертационной работе проведен обзор методов решения системы линейных алгебраических уравнений (СЛАУ) Несамосопряженность оператора порождает ряд сложностей при применении многих численных методов В сеточных методах несамосопряженным операторам, как правило, соответствуют несимметричные матрицы СЛАУ, что затрудняет решение таких систем Для решения СЛАУ, получаемых при дискретизации нелинейных задач, одними из наиболее эффективных методов являются проекционные методы с использованием проектирования на подпространства Крылова1

В первой главе сформулированы математические постановки задач о внутренних стационарных течениях вязкой и идеальной несжимаемой жидкости в естественных переменных вектор скорости - давление На основе этих математических моделей построены вариационные постановки в форме Галеркина, эквивалентные исходным задачам. Глава состоит из трех параграфов

1 Saad Y Iterative methods for sparse linear systems - PWS Publishing Company, 1996

В параграфе 1.1 приведена математическая постановка задачи об установившимся движении вязкой несжимаемой жидкости через ограниченную односвязную область , которая заключается в определении вектора скорости и и давления р, удовлетворяющих в О. уравнениям Навье - Сто-кса

уравнение движения

—уУ2и + (и• \7)и + Ур = Г, хе£2, уравнение неразрывности

У-и = 0, хеП,

и краевым условиям на границе области Г. заданному вектору скорости на входе в область, условиям непротекания на твердых стенках и заданной нормальной компоненте тензора напряжения на выходе. Уравнения Навье - Стокса с введенными на границе области краевыми условиями имеют решение, определяемое с точностью до произвольной постоянной для давления

В случаях больших значений кинематической вязкости или при малых значениях вектора скорости решение уравнений Навье — Стокса может быть заменено решением уравнений Стокса: уравнение движения

-л>У211 + Ур = £, хей, уравнение неразрывности

V и = 0, хе£2, с соответствующими краевыми условиями на границе области

В параграфе 1.2 приведена математическая постановка задачи об установившимся протекании идеальной несжимаемой жидкости через ограниченную односвязную область О, которая заключается в определении вектора скорости и и давления р, удовлетворяющих в О уравнениям Эйлера

уравнение движения

(и + , хеО, уравнение неразрывности

У-и = 0, хеП,

и краевым условиям на границе области заданному вектору скорости на входе в область, условиям непротекания на твердых стенках и заданной нормальной компоненте вектора скорости на выходе

В параграфе 1.3 для поставленных задач сформулированы вариационные постановки в форме Галеркина, эквивалентные исходным задачам

Вторая глава посвящена построению дискретных аналогов вариационных постановок в форме Галеркина, полученных в первой главе Представлена разработанная вычислительная схема расчета неизвестных компонент вектора скорости и давления в гипервекторе, позволяющая единообразно моделировать стационарные внутренние течения вязкой несжимаемой жидкости (уравнения Навье - Стокса, Стокса) и идеальной несжимаемой жидкости (уравнения Эйлера). Выписаны стабилизированные схемы Галеркина для моделирования трехмерных течений на тетраэдральном разбиении Построены противопотоковые схемы с расчетом неизвестных в узлах конечных элементов и серединах их ребер Для аппроксимации потоков через ребра двойственной сетки предложено использовать формулы интегрирования базисных функций по отрезкам медиан конечного элемента Предложен способ учета краевого условия и • п|г — Ъ (п - вектора

нормали к границе Г, к - заданная функция) в конечно-элементных аппроксимациях уравнений Эйлера

В параграфе 2.1 введены двойственные по отношению к первичной сетке разбиения 1) разбиение на ячейки, построенные вокруг вершин симплексов первичной сетки, 2) разбиение на ячейки, построенные вокруг середин сторон симплексов первичной сетки.

В параграфе 2.2 построены дискретные аналоги полученных в первой главе вариационных постановок

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

В параграфе 2.3 отмечены допустимые дискретные пространства базисных функций для смешанных конечно-элементных постановок Пространства выбираются в соответствии с условием Ладыженской - Бабушки

- Бренди2, которое является одним из необходимых условий существования и единственности решения вариационной задачи.

(Уу*,?*)

inf -! sup

Л IIIIW*

•>p>Q,

(1)

где вектор скорости принадлежит пространству Vй, а давление - пространству Р11, /3 - константа, не зависящая от к, ( , ) - скалярное произведение Показано, что пары пространств Тейлора - Худа и Крузея -Равьяра соответствуют условию Ладыженской - Бабушки - Брецци

В параграфе 2.4 представлены вычислительные схемы расчета неизвестных компонент вектора скорости и давления в гипервекторе В блочно-матричном виде система уравнений, основанная на смешанной конечно-

„ п+1

элементной постановке, для неизвестных компонент вектора скорости и и давления р"+1 имеет вид 1

a + -qm;qt

£ Р

Q

гм,

р у

[~un+r Mf

Uw+ij eMvp\

(2)

где р° — заданное начальное приближение давления, е > 0 — штрафной параметр, п > 0 - шаг итерации, в - диффузионная матрица, с - конвективная матрица, матрица а = в + с соответствует аппроксимации уравнений Навье - Стокса, матрица а = 1) аппроксимации уравнений Сто-кса, матрица а = с - аппроксимации уравнений Эйлера, матрица соответствует дискретному уравнению неразрывности, мр - матрица массы

в дискретном уравнении неразрывности, м - матрица массы в дискретном уравнении движения, £ - вектор правой части Для обращения матрицы массы мр выполнена ее диагонализация

Другим способом учета уравнения неразрывности является алгоритм Удзавы Матричный вид модифицированного алгоритма Удзавы для нахождения вектора скорости и"+1 и давления р"+> представлен выражениями

2 Brezzi F , Fortin M Mixed and hybrid finite element methods Springer - Verlag, New York, 1991

(A + 1qMp1Qt)u"+1 sMf + Qp",

* (3)

Pn+l = p" --M;1QTU"+1,

s

где p° - заданное начальное приближение давления, s > 0 - штрафной параметр, п> 0 - шаг итерации, матрицы А, D, С, Q, Мр, М и вектор f определены в алгоритме (2)

В рассмотренных матричных уравнениях несжимаемое слагаемое штрафуется пропорционально давлению, что дает возможность отдельно получить уравнения для расчета компонент вектора скорости и уравнение для вычисления давления

Следующие алгоритмы построены на основе стабилизированной схемы Галеркина, в которой введены два члена Первый член соответствует методу Streamline Upwind Petrov - Galerkin (SUPG), второй - методу стабизи-зации давления Pressure Stabilization Petrov - Galerkin (PSPG) Стабилизирующие члены входят в постановку с соответствующими параметрами стабилизации, позволяющими регулировать величину искусственной вязкости, добавляемой в уравнения для подавления больших нефизических осцилляции В диссертационной работе предлагается определять параметр стабилизации в каждой ячейке ее размерами и заданными в ней коэффициентами уравнений

В блочно-матричном виде метод SUPG/PSPG для решения уравнений Навье - Стокса и Эйлера имеет вид

(4)

А + ~(0 + К^/га )

(<2 +К/это) ^ря-в ^ С

где и - вектор скорости, р - давление, обозначения матриц Б, С, С>, М введены ранее в (2), А = О+С соответствует аппроксимации уравнений Навье - Стокса, А = С - аппроксимации уравнений Эйлера, f - вектор правой части Матрицы стабилизации &зиР0, К, Сзига включают параметр стабилизации тзара, матрицы стабилизации НР8Рв, КР№С содержат параметр стабилизации тр5ра Для решения уравнений Стокса используется метод стабилизации давления РБРО, поскольку в уравнения не входят конвективные члены

Для реализации краевого условия и п| = Й в вариационных постановках уравнений Эйлера предлагается ввести матрицу преобразования

декартовых компонент вектора скорости и в тангенциальную и нормальную компоненты

где т'

^соэ в' -?,твЛ

со5,в'

п^Тй?^*),

=1

матрица преобразования, в' - угол между

осью х и касательным вектором к ребру в г -ом узле, й'1 = (инп1, м^ -нормальная и тангенциальная компоненты вектора скорости, т - число расчетных узлов на конечном элементе, (р1 (х) - базисная функция пространства Vй

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

Таким образом, в диссертационной работе для корректного учета дивергентного условия используются смешанные и стабилизированные постановки и соответствующие им вариационные задачи Учет конвективных членов выполнен противопотоковыми схемами с расчетом неизвестных в узлах и серединах ребер симплексов сетки и стабилизированными схемами Галеркина

Третья глава посвящена описанию разработанного комплекса программ, вычислительным экспериментам и анализу их результатов. Глава состоит из трех параграфов

В параграфе 3.1 описываются структура комплекса программ, структура данных, его основные модули и функциональные особенности Разработана структура данных для хранения элементов матриц СЛАУ, в которых неизвестным вектором является гипервектор, включающий компоненты вектора скорости и давления Исследованы структуры глобальных мат-

риц и их размерности для разных вариационных постановок Процесс решения задачи разделен на этапы, что позволяет в комплекс программ без особых затрат добавлять новые вычислительные алгоритмы и расширять круг задач

Параграф 3.2 посвящен тестированию реализованных методов на ряде модельных задач, проведению вычислительных экспериментов на задачах, имеющих практическую ценность

На двумерной задаче конвекции-диффузии с доминированием конвекции выполнено сравнение следующих методов метода Петрова - Галерки-на стабилизации давления со стабилизацией оператора конвекции в направлении линии тока, метода конечных элементов/конечных объемов, МКЭ и конечно-элементных противопотоковых схем с расчетом неизвестных в узлах и серединах ребер треугольников

Выполнено сравнение точности решения уравнений конвективно-диффузионного типа для трехмерной задачи, полученного методом конечных элементов и стабилизированным методом конечных элементов Для этой же тестовой задачи решено нестационарное уравнение конвективно-диффузионного типа Аппроксимация по времени выполнена схемой Кранка - Николсон

Программная реализация алгоритмов расчета полей вектора скорости и давления (смешанного метода конечных элементов, алгоритма Удзавы и стабилизированных схем Галеркина) верифицирована на модельных двумерных и трехмерных задачах Навье - Стокса (в стационарном и нестационарном случаях), Стокса и Эйлера Для алгоритмов, использующих функцию штрафа, определены оптимальные значения штрафного параметра

В качестве тестовой задачи численно решена задача о каверне Наблюдается хорошее совпадение полученных результатов с опубликованными данными3. В диссертационной работе приведены результаты расчетов каверны при Re = 400, 1000, 2500.

В параграфе 3.3 приведены результаты моделирования течений в двумерных и трехмерных криволинейных каналах с разворотом потока на 180 и 270 градусов Для криволинейного канала с разворотом потока на 180 градусов расчеты выполнены при различном задании сдвига скорости на входе в область, показано возникновение вторичного течения и зависи-

3 Elias R N, Coutinho ALGA., Martins MAD Inexact Newton - type methods for non-linear problems arising from the SUPG/PSPG solution of steady incompressible Navier-Stokes equations//J Braz Soc Mech Sci &Eng -2004 - Vol XXVI, No 3 -P 330-339

мость интенсивности вторичного течения от величины сдвига скорости на входе (рис. 1).

Сравнение результатов расчетов, полученных в переменных вектор скорости - давление методом конечных элементов и в переменных векторный потенциал - вектор вихря методом конечных разностей4, показало хорошее согласование.

а) б)

Рис. 1. Проекция вектора скорости на поперечное сечение канала при у=40, на входе задан вектор скорости со сдвигом: а) иу= 0,75+0,5г; б) иу= 0,9+0,2г.

Трубопроводы, используемые в промышленности, представляют собой систему, состоящую из трубы ("русла") и различных конструктивных включений (сужение, расширение и другое). В диссертационной работе приведены результаты моделирования течения в каналах переменного сечения (с расширением и сужением русла). Определено влияние соотношения узкого и широкого сечений на погрешность вычислений при измельчении сетки (рис. 2).

В приложениях приводятся компоненты локальных матриц и векторов правых частей, вычисленные на треугольных конечных элементах для кусочно-линейных и кусочно-квадратичных базисных функций, а также на тетраэдральных конечных элементах для кусочно-линейных базисных

4 Хакимзянов Г.С., Шокин Ю.И., Барахнин В.Б., Шокина Н.Ю. Численное моделирование течений жидкости с поверхностными волнами. - Новосибирск: Изд-во СО РАН, 2001.

функций. Для полученных в работе вариационных постановок представлены портреты генерируемых матриц СЛАУ.

>

X

а)

б)

Линии тока для каналов переменного сечения: а) расширение русла; б) сужение русла.

Рис. 2. Линии

В заключении сформулированы основные результаты диссертационной работы, выносимые на защиту:

• Разработана технология реализации противопотоковых схем на неструктурированных сетках с расчетом неизвестных в узлах конечных элементов и серединах их ребер для моделирования течений с преобладанием конвективного переноса. В контексте метода конечных элементов предложен новый способ аппроксимации потоков через границы ячеек двойственной сетки.

• Разработана и программно реализована вычислительная схема на основе смешанных постановок с введением функции штрафа, позволяющая единообразно моделировать стационарные внутренние течения вязкой несжимаемой жидкости (уравнения Навье - Стокса, Стокса) и идеальной несжимаемой жидкости (уравнения Эйлера). Предложен способ учета краевого условия и-п| =й (п - вектор нормали к границе Г, к - заданная функция) в вариационных постановках уравнений Эйлера. Разработана структура данных для хранения элементов матриц СЛАУ с расчетом неизвестных компонент вектора скорости и давления в гипервекторе.

• Реализованы стабилизированные схемы Галеркина для моделирования трехмерных течений на тетраэдральном разбиении, в которых используется метод стабилизации давления и в качестве стабилизирующего оператора выбирается конвективная часть исходного оператора.

• Показана зависимость интенсивности вторичного течения от величины сдвига скорости на входе. Определено влияние соотношения узкого и

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

Список основных работ по теме диссертации

публикация в издании, рекомендованном ВАК

1 Гобьпп А В , Шокина Н Ю Анализ вычислительных схем методов конечных элементов и конечных разностей для моделирования течений несжимаемой жидкости // Вычислительные технологии - 2006 - Том 11.-№6 С 22-31

публикации в трудах международных конференций

2 Гобыш А В , Шокина Н Ю , Шурина Э П Анализ вычислительных алгоритмов для уравнений Навье - Стокса, основанных на методе конечных элементов // Труды международной конференции по вычислительной математике МКВМ-2004 Часть 1 - Новосибирск Изд-во ИВМиМГ СО РАН, 2004 - С 467-472

3 Шурина Э П, Гобыш А.В Исследование вычислительных схем расчета несжимаемых течений на базе ММКО/МКЭ // Вычислительные технологии (2003, Т 7), Региональный вестник Востока (2003, Т 3) - Совместный выпуск по материалам Международной конференции "Вычислительные и информационные технологии в науке, технике и образовании" Ч. IV - Новосибирск - Алматы - Усть-Каменогорск - 2003 -С 13-17

публикации в сборниках научных трудов

4. Гобыш А В Конечно-элементные аппроксимации уравнений Навье -Стокса на базисных функциях различных порядков // Сборник научных трудов НГТУ - Новосибирск- Изд-во НГТУ, 2005 - Вып 1(39), С 101-106

5 Гобыш А В Трехмерные конечно-элементные аппроксимации уравнений Навье - Стокса, Стокса, Эйлера // Сборник научных трудов НГТУ.-Новосибирск Изд-во НГТУ, 2006-Вып 1(43) -С 55-60

6 Shurina Е Р , Gobysh А V Comparative analysis of finite element approximation for the Navier - Stokes equations with basis functions of different orders // Bulletin of Novosibirsk Computing Center Series Numerical Analysis-Novosibirsk NCC Publisher, 2005 -Iss 13 -p 77-86

Отпечатано в типографии Новосибирского государственного технического университета 630092, г Новосибирск, пр К Маркса, 20, тел /факс (383) 346-08-57 формат 84x60x1/16, объем , 0. Ьп.л., тираж 100 экз , заказ № ¿подписано в печать 03 10 2007 г.

ВВЕДЕНИЕ

Глава 1. ПОСТАНОВКА ЗАДАЧИ

1.1. Математическая модель движения вязкой несжимаемой жидкости: уравнения Навье - Стокса, уравнения Стокса.

1.2. Математическая модель движения идеальной жидкости: уравнения Эйлера.

1.3. Вариационные постановки.

Глава 2. ПОСТРОЕНИЕ ДИСКРЕТНЫХ АНАЛОГОВ ВАРИАЦИОННЫХ ПОСТАНОВОК

2.1. Триангуляция расчетной области, построение двойственных сеток

2.2. Дискретные аналоги вариационных постановок

2.3. Пространства интерполяционных функций в смешанных постановках

2.4. Алгоритмы расчета поля скоростей и давления.

2.4.1. Смешанный метод конечных элементов.

2.4.2. Алгоритм Удзавы.

2.4.3. Стабилизированные методы конечных элементов.

2.4.4. Модификация конвективных членов в противопотоко-вых схемах.

2.4.5. Противопотоковая схема с вычислением неизвестных в узлах конечных элементов.

2.4.6. Противопотоковая схема с вычислением неизвестных в серединах ребер конечных элементов.

Глава 3. РЕАЛИЗАЦИЯ АЛГОРИТМОВ МОДЕЛИРОВАНИЯ, РЕЗУЛЬТАТЫ ВЫЧИСЛИТЕЛЬНЫХ ЭКСПЕРИ

МЕНТОВ

3.1. Структура комплекса программ

3.1.1. Структура данных комплекса программ.

3.1.2. Основные модули комплекса программ

3.2. Тестирование методов.

3.2.1. Численное решение уравнений конвективно-диффузионного типа в двумерной области. Проти-вопотоковые схемы

3.2.2. Численное решение уравнений конвективно-диффузионного типа в трехмерной области.

3.2.3. Численное решение нестационарных уравнений конвективно-диффузионного типа в трехмерной области

3.2.4. Численное решение уравнений Навье - Стокса и Стокса в двумерной области.

3.2.5. Численное решение уравнений Навье - Стокса, Стокса и Эйлера в трехмерной области.

3.2.6. Численное решение задачи о каверне

3.3. Моделирование течений в криволинейных каналах и каналах переменного сечения.

3.3.1. Моделирование течения в двумерном криволинейном канале с разворотом потока на 270°.

3.3.2. Моделирование течения в криволинейном канале с разворотом потока на 180°

3.3.3. Моделирование течений вязкой несжимаемой жидкости в каналах переменного сечения.

В современной науке для изучения самых разнообразных природных явлений, технологических процессов и экологических проблем широко используется математическое моделирование [36]. Исходный объект заменяется математической моделью, которая исследуется средствами вычислительной математики. Математическое моделирование сочетает в себе многие достоинства как теории, так и эксперимента [36]. Работа с моделью объекта дает возможность относительно быстро и без существенных затрат исследовать его свойства и поведение в различных ситуациях (преимущества теории). В то же время вычислительные эксперименты с моделями объектов позволяют подробно и глубоко изучать объекты в достаточной полноте, недоступной теоретическим подходам (преимущества эксперимента). Совместное использование физического и вычислительного экспериментов в исследовании объекта позволяет, с одной стороны, уменьшить количество натурных дорогостоящих измерений, а с другой стороны - провести верификацию и усовершенствование математических моделей.

Широкий круг проблем, стоящих перед современной наукой и техникой, связан с изучением движения жидкости. Для описания явлений используются модели стационарных и нестационарных течений, вязкой и невязкой жидкости, пограничных слоев и другие [26, 42].

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

При подборе методов решения необходимо учитывать особенности исследуемых течений и соответствующих математических моделей. В диссертационной работе рассматриваются стационарные, внутренние, однофазные, химически однородные течения. Движение вязкой несжимаемой жидкости описывается уравнениями Навье - Стокса. При создании комплекса программ необходимо наличие не только программных реализаций сложных моделей, но и реализации более простых приближенных моделей, с помощью которых можно получать предварительное представление о характере течения. В случае больших значений кинематической вязкости или при малых значениях вектора скорости решение уравнений Навье - Стокса может быть заменено решением уравнений Стокса. Решения уравнений Эйлера, описывающие течение идеальной несжимаемой жидкости, могут служить начальным приближением для итерационных методов решения уравнений Навье - Стокса. Поэтому в диссертационной работе рассмотрен и включен в комплекс программ набор моделей, что помогает более полно исследовать явление.

Существуют следующие математические формулировки перечисленных уравнений, описывающих движение несжимаемой жидкости [40, 125]:

• в естественных переменных вектор скорости - давление (и - р формулировка);

• в переменных векторный потенциал - вектор вихря (ф - и формулировка);

• в переменных вектор скорости - вектор вихря (и - ш формулировка).

Каждая постановка имеет свои преимущества и недостатки. Исторически описание в переменных функция вихря - функция тока было весьма популярно при расчете двумерных течений вязкой несжимаемой жидкости [28, 30, 39, 40, 94, 125]. Трудности обобщения этого подхода на трехмерный случай связаны с установлением граничных условий для трехмерного аналога функции тока - векторного потенциала [30, 40, 125]. Это является одной из причин, по которой ф - ш формулировки использовались относительно редко при численном моделировании трехмерных течений. Однако, переход к зависимым переменным фиш имеет ряд преимуществ по сравнению с и - р подходом. Так, для ф - ш формулировок уравнение неразрывности на дифференциальном уровне выполняется автоматически, а на разностном -при подходящей аппроксимации компонент вектора скорости и использовании разнесенных сеток. С введением векторного потенциала (функции тока в двумерном случае) давление не входит в число зависимых переменных, поэтому на промежуточных итерациях давление исключается из вычислительного процесса и при необходимости может быть восстановлено после сходимости —* итераций для ф и и. Примеры вычислительных алгоритмов, основанных на ф - и) формулировках, представлены в [40, 121, 127].

В работах [74, 104, 105] используется u — w формулировка, в которой поле скоростей вычисляется путем решения уравнений Пуассона для каждой компоненты вектора скорости. В статье [77] на каждом итерационном шаге компоненты вектора скорости находятся из решения уравнения неразрывности и соотношений, задающих связь между вектором скорости и вектором вихря. В трехмерном случае u — w подход менее эффективен, чем описание в естественных переменных, если только движение вихря, в особенности нестационарного, не представляет особого интереса. Более подробный обзор данных подходов можно найти в книге [40].

В диссертационной работе рассмотрены математические модели в естественных переменных вектор скорости - давление.

Система уравнений Навье - Стокса состоит из уравнений неразрывности и движения. Преобладание конвективного или диффузионного переноса в течении жидкости определяется числом Рейнольдса. Уравнения движения и неразрывности образуют систему связанных уравнений с неизвестными величинами - компонентами вектора скорости и давления. Для несжимаемых течений важно влияние перепадов давления на поле скоростей, способ же расчета давления не очевиден, поскольку давление не является термодинамической переменной и для его определения не существует отдельного уравнения. Отметим, что при использовании моделей в естественных переменных на уровне вычислительной схемы возникают трудности корректного учета условия, накладываемого на дивергенцию аппроксимируемого поля скоростей, так называемого дивергентного условия.

Численное решение уравнений Навье - Стокса, Стокса и Эйлера в естественных переменных предполагает решение следующих задач, определяющих структуру комплекса программ [14, 17, 28]: пространственной дискретизации уравнений; выбора базисных функций для вектора скорости и и давления р, гарантирующих существование решения (и,р)] учета нелинейности, связанной с конвективными членами в уравнении движения; выбора алгоритма расчета неизвестных вектора скорости и давления; решения системы линейных алгебраических уравнений (СЛАУ).

На сегодняшний день ведущие позиции при численном решении дифференциальных уравнений занимают сеточные методы. В зависимости от способа построения приближенного решения можно выделить следующие классы сеточных методов, часто используемые в вычислительной гидродинамике [18, 24, 28, 39, 40, 41, 46]: дифференциальные и вариационные.

Отличительной характеристикой дифференциальных методов является непосредственная аппроксимация исходного дифференциального уравнения, либо соответствующих интегробалансных соотношений, с помощью разностных аналогов производных. В этот класс входят метод конечных разностей (МКР) и метод конечных объемов (МКО).

Идея МКР состоит в замене расчетной области множеством дискретных точек (сеткой) и построении дискретизации дифференциальных операторов в узлах этой сетки. Существует более десятка конечно-разностных схем для решения уравнений Навье - Стокса, отличающихся различными признаками [18, 24, 28, 35, 39, 40, 41, 46, 92]. Особенностью этих схем является наличие различных сеточных параметров, рациональный выбор которых трудоемок [28]. К недостаткам метода относятся также сложность построения конечно-разностных операторов, аппроксимирующих уравнения Навье - Стокса на произвольном сеточном шаблоне, учет краевых условий и адаптация регулярных сеток к областям с произвольной криволинейной границей [28]. Среди преимуществ МКР - низкие вычислительные затраты и возможность применения к широкому спектру решаемых задач. Применение в контексте метода конечных разностей разнесенных („шахматных") сеток для вычисления компонент вектора скорости и давления, впервые предложенное в 1965 г. Харлоу и Велч [81], позволяет связать значения вычисляемых неизвестных в соседних точках. Осцилляции, как правило, возрастают при увеличении числа Рейнольдса, поскольку диссипативные члены, посредством которых осуществляется связь значений компонент вектора скорости в соседних точках, в этом случае малы. Случай разнесенных расстановок переменных эквивалентен применению разных порядков базисных функций для обеспечения условия существования и единственности решения (и,р) [30, 91].

В отличие от МКР, метод конечных объемов [34] связан с построением дискретизации интегральных формулировок на заданной сетке. Преимуществом конечно-объемных аппроксимаций являются простота реализации и возможность работы с неструктурированными сетками. Для построения устойчивых конечно-объемных аппроксимаций в ряде работ используются пространственные „шахматные" сетки. В работах [110, 116] была предложена разнесенная дискретизация для неструктурированных сеток. В работе [110] для компонент вектора скорости используются одни и те же контрольные объемы, а вклад градиента давления в уравнение движения определяется с помощью специальной реконструкции давления на разнесенной сетке.

Необходимо отметить, что использование „шахматных" сеток в МКР и МКО позволяет повысить устойчивость схем и на уровне вычислительной схемы корректно учесть дивергентное условие.

Другой класс сеточных методов - методы, основанные на использовании вариационных постановок, эквивалентных исходной дифференциальной краевой задаче. В этот класс входят различные варианты метода конечных элементов (МКЭ) [28, 29, 37, 38, 45, 51, 57, 58, 59, 60]. Одними из важных достоинств МКЭ для задач гидродинамики является эффективность конструирования конечно-элементных аппроксимаций на неструктурированных сетках и возможность практически полной автоматизации всех этапов решения задачи - от построения сетки до сборки глобальной матрицы результирующей СЛАУ [45, 51]. Высокую вычислительную эффективность МКЭ обеспечивает выбор функций кусочно-полиномиального вида [29, 38].

В последние годы создано множество методов моделирования течений несжимаемой жидкости с использованием МКЭ. Важным критерием эффективности таких схем является точность учета дивергентного условия. Наиболее распространенным подходом, позволяющим учесть дивергентное условие, является введение смешанных вариационных постановок [15, 128, 129]. Под термином „смешанные постановки" понимаются такие постановки, в которых аппроксимация происходит на паре специально подобранных пространств. Пары пространств, используемых в смешанных вариационных постановках, необходимо подбирать в соответствие с условием Ладыженской - Бабушки - Брецци (ЛББ), которое является одними из необходимых условий существования и единственности решения [15, 27, 39]. Интерполяционные функции в смешанных постановках, удовлетворяющие условию ЛББ, называются div-устойчивыми (равномерно устойчивыми) [15].

Имеются и другие подходы [40] к удовлетворению дивергентного условия. и

Например, появление методов штрафа дало начало новому направлению в области моделирования течений вязкой несжимаемой жидкости методом конечных элементов - численным алгоритмам, допускающим произвольные сочетания интерполяционных полиномов решения. Соответствующие методы называются методами, „обманывающими" условие ЛББ. Применение этих методов [39, 56, 57, 64] позволяет исключить явное наличие давления в уравнении движения.

К классу методов с произвольным сочетанием интерполяционных полиномов решения также относятся стабилизированные схемы Галеркина [61, 83, 84, 97, 98, 120]. Назначение параметров стабилизации в этих методах - регулировать величину искусственной вязкости, добавляемой в уравнения для подавления больших нефизических осцилляций. Способы задания параметров стабилизации, зависящих от пространственных и временных характеристик, приводятся в [120]. Получение оптимального значения численной диффузии в направлении линий тока также исключает возможность возникновения осциллирующего поля давления.

Отметим также, что хорошие результаты дает алгоритм Удзавы [51, 52, 64], традиционно используемый при решении задачи о седловой точке. Для улучшения свойств сходимости алгоритма Удзавы в уравнение неразрывности вводится слагаемое, содержащее давление и штрафной параметр [52, 64, 117, 118, 120].

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

Для течений, в которых доминирует конвекция, использование центральных разностей приводит к появлению в решении сильных нефизических осцилляций [30]. При попытке стабилизировать решение путем дискретизации конвективных членов двухточечными разностями против потока или суммами с весами центрально-разностных выражений и разностей против потока получается неточное решение, особенно в случаях совпадения локального направления вектора скорости с ориентацией (левой, правой) сетки и больших локальных градиентов вектора скорости.

Более точное решение получается, если для аппроксимации конвективных членов используются разности высокого порядка точности против потока, подобные четырехточечным формулам. В комбинированной схеме [30] в зависимости от абсолютной величины сеточного числа Пекле сочетаются центрально-разностная и противопоточная схемы аппроксимации конвективных членов уравнения. Схема квадратичной интерполяции против потока [96] (схема Леонарда) имеет второй порядок точности, но менее устойчива при расчетах, чем комбинированная схема. В работе Дэвиса и Мура [65] для решения нестационарных уравнений Навье - Стокса используются многомерные разности третьего порядка точности против потока, которые являются обобщением одномерной квадратичной интерполяции против потока Леонарда [96]. Схема с косыми разностями против потока Рейсби позволяет учитывать направление вектора скорости потока по отношению к грани контрольного объема [110].

Далее перечислим противопотоковые схемы на базе методов конечных элементов и конечных объемов. Наиболее исследованными и часто применяемыми в случае неструктурированных сеток являются противопотоковые схемы, основанные на методе конечных объемов с расчетом неизвестных в центрах ячеек [53, 69]. Несмотря на то, что грани элементов разбиения не параллельны координатным осям, данные схемы в большинстве случаев являются одномерными [111]. Расчеты по подобным схемам не воспроизводят многомерную структуру потока и обладают чрезмерной численной диффузией [69]. Для построения противопотоковых схем второго порядка необходимо расширение шаблона, что для неструктурированных разбиений приводит к усложнению структур данных [53, 69].

Противопотоковые схемы с расчетом неизвестных в узлах сетки и серединах ее ребер немногочисленны [99, 108, 111]. В ряде случаев противопото-ковый принцип аппроксимации сводится к использованию одного значения скалярной субстанции (в узле симплекса, лежащего против потока) [93], либо двух взвешенных значений (на концах ребра симплекса, лежащих против потока) [99, 111].

Схема, учитывающая направления потока (Flow Oriented Upwind Scheme), использует преимущество расчета неизвестных в узлах - возможность построения асимметричных функций профиля. Но расчеты по данной схеме признаны неудовлетворительными и итерационные процессы часто расходятся [99], поскольку схема не обладает свойством положительности.

Конечно-элементный аналог схем против потока, предложенный Кристи [62], состоит в использовании функций формы, модифицированных с тем, чтобы вклады каждого из элементов распределялись в зависимости от доминантного направления течения на элементе, а не равномерно по его узлам. Однако, предложенные модификации вносят чрезмерную численную диффузию в направлении перпендикулярном вектору средней скорости на элементе.

В работах Келли и соавторов [89] в методе Петрова - Галеркина со стабилизацией оператора конвекции в направлении линии тока [83, 84, 122] предложено использовать функции формы, искривленные в направлении средней скорости на элементе. Другим подходом является метод стабилизирующей тензорной вязкости [78], позволяющий работать с конечно-элементными пространствами младших порядков. Выражения для стабилизирующей тензорной вязкости учитывают направления средней скорости на элементе. Анализ этого метода показывает [78], что в нем не отражается многомерная структура потока, поэтому привносимая численная диффузия велика.

Другим эффективным подходом к решению задач с преобладанием конвекции является смешанный метод конечных элементов/объемов, в котором используются метод Галеркина с симметричными тестовыми функциями для самосопряженной части дифференциальных операторов и противопотоковые схемы МКО - для несимметричной их части [67, 69].

Метод конечных суперэлементов (МКСЭ) [19, 20, 21] был предложен как метод дискретизации дифференциальных уравнений, учитывающих мелкомасштабные (относительно шага сетки) неоднородности задачи (стержни в ядерном реакторе, инородные включения в композитных материалах и т. п.) или, если явных неоднородностей нет, то присутствие различных пространственных масштабов решения. Основные вычислительные затраты МКСЭ связаны с построением базиса на элементе разбиения. С другой стороны, при использовании схем высокого порядка МКСЭ включает в себя процедуру локального исключения неизвестных, ассоциированных с внутренними узлами элемента. Такое исключение приводит к уменьшению размерности результирующей линейной системы уравнений, улучшению ее обусловленности и уменьшению вычислительных затрат на ее решение.

В разрывном методе Галеркина [50, 63, 66, 70] решение аппроксимируется на конечно-элементной сетке кусочно-полиномиальными функциями, разрывными на межэлементной границе. Главной особенностью этих методов является то, что решение удовлетворяет закону сохранения локально на каждом элементе (поэлементно), степени свободы решения ассоциируются с элементами пространственной дискретизации, а не с вершинами элементов как в классическом методе Галеркина. Несмотря на большой интерес, проявляемый к методу и его модификациям, для него еще не создано единой вычислительной технологии [70].

Основная особенность класса противопотоковых конечно-элементных схем, использующих аналог конечно-объемной двойственной сетки, состоит в построении ячеек двойственной сетки, на которых используется кусочно-постоянный оператор модификации конвективных членов. В данном классе выделяются противопотоковые схемы с расчетом неизвестных в вершинах симплексов [80] и в серединах ребер симплексов [103] первичного разбиения расчетной области.

Методы аппроксимации по времени. Несмотря на то, что диссертационная работа посвящена стационарным течениям, проведено исследование разработанных вычислительных схем на нестационарных задачах. По этой причине рассмотрим схемы аппроксимации по времени, используемые при численном решении дифференциальных уравнений. В литературе (см., например, [16, 35]) схемы аппроксимации по времени классифицируются на явные, неявные и полунеявные.

К классу явных схем относятся, например, схемы Адамса - Башфорта и Рунге - Кутта [122]. В явных схемах эллиптическая часть исходного уравнения вычисляется на предыдущих временных слоях, и матрица массы имеет ненулевые элементы только на главной диагонали. В данном случае нет необходимости решать СЛАУ на каждом временном слое, так как компоненты вектора решения явным образом выражаются через компоненты вектора решения на предыдущем временном слое. Диагональную матрицу массы можно получить при использовании подходящей дискретизации. В работе [127] предлагается набор базисных функций, которые позволяют получить диагональную матрицу массы. Показано, что при использовании таких базисных функций скорость расчета увеличивается в три раза, однако количество базисных функций и размерность СЛАУ при этом также увеличивается в три раза.

Если матрица массы, получаемая при дискретизации, не является диагональной, то выполняется ее диагонализация [113], состоящая в поэлементном преобразовании матрицы массы, при котором диагональным элементам присваиваются значения суммы всех элементов в строке, а внедиагональные элементы полагаются нулевыми. В диссертационной работе поэлементное преобразование матрицы массы с целью ее диагонализации применяется для обращения матрицы массы в вычислительных схемах расчета вектора скорости и давления.

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

Неявные схемы характеризуются вычислением матрицы жесткости на текущем временном слое, следовательно, на каждом временном слое требуется решение СЛАУ. Однако неявные схемы более устойчивы, чем явные, и порой оказываются более экономичными по времени счета, поскольку для устойчивости явных схем требуется большее число временных слоев. К классу неявных схем относятся, например, схемы с обратным дифференцированием и схемы типа Адамса - Моултона, среди которых можно выделить неявную схему Эйлера и схему Кранка - Николсон [25, 122].

Полунеявные схемы сочетают в себе достоинства явных и неявных схем. К этим схемам можно отнести методы типа „предиктор-корректор", в которых в качестве начального приближения на каждом временном шаге неявной схемы служит решение на этом же шаге, полученное с помощью явной схемы. В работе [55] предложен новый метод подобного класса, использующий для решения СЛАУ только некоторое фиксированное число итераций метода обобщенных минимальных невязок. Число итераций обычно выбирается небольшим, что значительно сокращает вычислительные затраты на решение СЛАУ. Полученная схема может рассматриваться как явная, поэтому приходится накладывать ограничения на шаг по времени.

Одним из способов повышения эффективности схем аппроксимации по времени является метод Тейлора - Галеркина [48, 71]. Эта схема основана на разложении искомой функции в ряд Тейлора по времени, последующей замене временных производных пространственными, и применении к полученной системе метода Галеркина. На основе метода Тейлора - Галеркина предложен ряд явных и неявных схем [71]. Изначально метод предлагался для задач конвективного переноса, но в последствии был применен и в других приложениях [71].

Несмотря на разнообразие подходов к аппроксимации дифференциальных задач, четких рекомендаций по выбору оптимальной схемы для конкретной задачи не существует. Этот выбор должен делаться на основании многих факторов, таких как физические и геометрические особенности решаемой задачи, метод пространственной дискретизации, вычислительные ресурсы и другие.

Решение систем линейных алгебраических уравнений. Решение систем линейных алгебраических уравнений, получаемых в результате дискретизации нелинейных задач, является важным этапом в математическом моделировании. Несамосопряженность оператора исходной дифференциальной постановки порождает ряд сложностей при применении многих численных методов. В сеточных методах несамосопряженным операторам, как правило, соответствуют несимметричные матрицы СЛАУ, что затрудняет решение таких систем.

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

Более эффективным и устойчивым среди итерационных методов решения таких систем уравнений является класс проекционных методов с использование проектирования на подпространства Крылова [23, 100, 113, 114, 123]. Общий подход к построению проекционных методов заключается в применении условия Петрова - Галеркина, которое означает ортогональность вектора невязки некоторому n-мерному подпространству. Наиболее известными являются следующие методы [23, 100, ИЗ, 114, 123]: метод сопряженных градиентов и метод полной ортогонализации, применяемые для решения СЛАУ с симметричной матрицей, а также метод обобщенных минимальных невязок, метод бисопряженных градиентов и метод бисопряженных градиентов со стабилизацией - для случая несимметричных матриц.

О предмете и содержании диссертации. Диссертация состоит из введения, трех глав, заключения, списка использованных источников (131 наименование), приложения. Работа изложена на 138 страницах, включая 23 таблицы и 40 иллюстраций.

ЗАКЛЮЧЕНИЕ

Математическое моделирование процессов становится одним из наиболее распространенных методов исследования объектов и явлений различной природы. Широкий круг задач, стоящих перед современной наукой и техникой, связан с решением уравнений движения жидкости. В диссертационной работе рассматриваются стационарные, внутренние, однофазные, химически однородные течения. Движение вязкой несжимаемой жидкости описывается уравнениями Навье - Стокса. Для полноты исследования явлений в рассмотрения включены уравнения Стокса и Эйлера. Для численного решения уравнений Навье - Стокса, Стокса и Эйлера разработаны вычислительные алгоритмы и получены результаты, соответствующие поставленным целям.

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

Моделирование течения идеальной жидкости требует задание нормальной компоненты вектора скорости на выходе. В дискретных вариационных постановках уравнений Эйлера при использовании базисных функций пространств Тейлора - Худа в явном виде не задается краевое условие для нормальной компоненты вектора скорости на границе области. В диссертационной работе предложен способ учета краевого условия (нормальной компоненты вектора скорости на выходе), заключающийся в построение матрицы преобразования декартовых компонент вектора скорости в тангенциальную и нормальную компоненты скорости.

При численной аппроксимации уравнений Навье - Стокса возникает трудности с учетом конвективных членов. Разработаны конечно-элементные про-тивопотоковые схемы линеаризации конвективных членов на треугольных сетках с расчетом неизвестных в узлах конечных элементов и серединах их ребер. Для аппроксимации потоков через границы ячеек двойственной сетки предложено использовать формулы интегрирования базисных функций по отрезкам медиан конечного элемента.

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

Проведено моделирование поля течения в каверне при Re = 400,1000,2500. Полученных результаты хорошо согласуются с данными, опубликованными другими авторами.

Выполнены расчеты поля скоростей и давления в криволинейных каналах с разворотом потока на 180° и 270°. Расчеты в криволинейном канале с разворотом потока на 180° выполнены при задании сдвига скорости на входе. Показана зависимость интенсивности вторичного течения от величины сдвига скорости на входе, при увеличении сдвига скорости на входе увеличивается интенсивность вторичного течения, характер течения существенно не меняется. Выполнено сравнение полученных результатов с данными, опубликованными другими авторами.

Трубопроводы, используемые в промышленности, представляют собой систему, состоящую из трубы („русла") и различных конструктивных включений (заслонки, резкие сужения, расширения и другое). Проведена серия расчетов для каналов с расширением и сужением русла. Представлены результаты расчетов для каналов переменного сечения, имеющих различное соотношение между узким и широким сечениями. Определено влияние соотношения узкого и широкого сечений на погрешность вычислений при измельчении сетки в трубопроводах, состоящих из различных конструктивных включений (элементов сужения и расширения).

1. ГОБЫШ А. В. Расчет тепловых процессов смешанным МКЭ и ММКО // Наука. Техника. Инновации. Региональная научная конференция студентов, аспирантов и молодых ученых: Тез. докладов в 5 частях. Часть 1. - Новосибирск: Изд-во НГТУ, 2002. - С. 147.

2. ГОБЫШ А. В. Анализ алгоритмов Удзавы для решения системы уравнений Навье Стокса //IV Всероссийская конференция молодых ученых по математическому моделированию и информационным технологиям: программа и тезисы докладов. - Красноярск. - 2003. - С. 19.

3. ГОБЫШ А. В. Учет давления при решении уравнений Навье Стокса с использованием МКЭ // Наука. Технологии. Инновации. - Всероссийская научная конференция молодых ученых: Тез. докладов в б частях. Часть 1. - Новосибирск: Изд-во НГТУ, 2003. - С. 103-104.

4. ГОБЫШ А. В., ШОКИНА Н. Ю. Анализ вычислительных схем методов конечных элементов и конечных разностей для моделирования течений несжимаемой жидкости // Вычислительные технологии. 2006. -Том 11. - Ш. С. 22-31.

5. Багаев Б. М., Карепова Е. Д., Шайдуров В. В. Сеточные методы решения задач с пограничным слоем. Новосибирск: Наука, 2001. -240 с.

6. Госмен А. Д., Пан В. М., Ранчел А. К., Сполдинг Д. Б., Вольфштейн М. Численные методы исследования течений вязкой жидкости.- М.: Мир, 1972. 323 с.

7. Жуков В. Т., Новикова Н. Д., Феодоритова О. Б. Итерационный метод для конечно элементных схем высокого порядка. Часть 2. М.: Препринт ИПМ им. М.В. Келдыша, № 43, 2003. - 32 с.

8. ИГНАТЬЕВ А. А. Разностная схема для вязкости Навье Стокса // Математическое моделирование. - 2001. - Т. 13, № 8. - С. 107-116.

9. Ильин В. П. Методы неполной факторизации для решения алгебраических систем. М.: Физматлит, 1995. - 288 с.24. ковеня в. М., яненко Н. Н. Метод расщепления в задачах газовой динамики. Новосибирск.: Наука. - 1981. - 304 с.

10. Математическое моделирование конвективного тепломассообмена на основе уравнений Навье-Стокса / В.п. полежаев, А.В. бунэ, Н.А. верезуб и др. М.: Наука, 1987.

11. РОУЧ П. Вычислительная гидродинамика. М.: Наука, 1977. - 656 с.

12. Самарский А. А., Вабищевич П. H. Математическое моделирование и вычислительный эксперимент. М.: ИММ РАН, 2000. 409 с.

13. СЕГЕРЛИНД Д. Применение метода конечных элементов. М.: Мир, 1979. - 392 с.38. сьярле Ф. Метод Конечных элементов для эллиптических задач. -М.: Мир, 1980. 512 с.

14. Темам Р. Уравнения Навье Стокса. Теория и численный анализ. -М.: Мир, 1981.

15. ФЛЕТЧЕР К. Вычислительные методы в динамике жидкостей: В 2 х т.: Пер. с англ. - М.: Мир, 1991.

16. Шурина Э. П., Солоненко О. П., Войтович Т. В. Новая технология метода конечных объемов на симплициальных сетках для задач конвективно-диффузионного типа. Новосибирск, 1999. (Препринт ИТ-ПМ СО РАН; Щ.

17. Шурина Э. П., Соловейчик Ю. Г. Решение краевых задач в составных областях: Учеб. Пособие/ Новосиб. Электротех. Ин-т. Новосибирск, 1986. - 73 с.

18. Численное решение многомерных задач газовой динамики / С. К. годунов, а. в. Забродин и др. м.: Наука, 1976. - 400 с.

19. Хакимзянов Г. С., Шокин Ю. И., Барахнин В. В., Шокина Н. Ю. Численное моделирование течений жидкости с поверхностными волнами. Новосибирск: Изд-во СО РАН, 2001.

20. Ambrosi D, Quartapelle L.A. Taylor-Galerkin method for simulating nonlinear dispersive water wave //J. Comput. Phys. 1998. - Vol. 146. -P. 546-569.

21. Automatic mesh domain partitioning for adaptively refined grids/ Leender Willem Vijfvinkel, Ph. D. Thesis Katholieke universiteit Nijmegen, 2001.

22. BARTON I.E. The entrance effect of laminar flow over a backward facing step geometry // Int. J. Numer. Methods Fluids. 1995. Vol. 25. P. 633-644.

23. Brown D. L., Cortez R., Minion M. L. Accurate projection method for the incompressible Navier-Stokes equations // J. Comput. Phys. 2001.- Vol. 168. P. 464-499.

24. Canuto C., Hussaini M. Y., Quarteroni A., Zang T. A. Spectral methods in fluid dynamics. Springer Verlag, New York - Berlin, 1987.

25. Carey G. F., Pehlivanov A. I., Shen Y., Bose A., Wang К. C. Least-squares finite elements for fluid flow and transport // Int. J. Numer. Methods Fluids. 1998. - Vol. 27. - P. 97-107.

26. Dawis R. W., Moore E. F. // J. Fluid Mech. 1982. - Vol. 116. -p. 475-506.66. dawson C., Proft J. Coupling of continuous and discontinuous Galerkin methods for transport problems // Comput. Methods Appl. Mech. Engrg.- 2002. Vol. 191. - P. 3213-3231.

27. Debiez C., Dervieux A., Mer K. and Nkonga B. Computation of unsteady flows with mixed finite volume/finite element upwind methods // Int. J. Numer. Methods Fluids. 1998. - Vol. 27. - P. 193-206.

28. Deconinck H., Degrez G. Multidimensional upwind residual distribution schemes and applications // Proc of Second Intern. Symp. On Finite Volumes for Complex Applications, 1999, Duisburg, Germany.- HERMES Science Publications, Paris, P. 27-40.

29. DELANAYE M., ESSERS J. A. Quadratic-reconstruction Finite volume scheme for Compressible flows on unstructured adaptive grids // AIAA Journal, 1997. Vol. 35. P. 631-639.

30. Dolejsi V., Feistauer M. and Schwab C. On some aspects of the discontinuous Galerkin finite element method for conservation laws. Preprint No. MATH-KNM-2001/5.

31. Franca L. P., Russo A. Unlocking with residual free bubbles // submitted to: Comput. Methods in Appl Mech. and Engrg.

32. Gatcki Т. В., Grosch С. Е., Rose М. Е. The numerical solution of the Navier Stokes equations for 3 - dimensional, unsteady, incompressible flow by compact schemes //J. Comput. Phys. - 1989. - Vol. 82. - P. 298-329.

33. Gresho P. M., Chan S. Т., Lee R. L., Upson G. D. A modified finite element method for solving the time-dependent incompressible Navier-Stokes equations. Part 2. Applications // Int. J. Numer. Methods Fluids. 1984. Vol. 4. P. 619-640.

34. Girault v., Raviart P. A. Finite element methods for Navier Stokes equations. Springer - Verlag, New York - Berlin, 1986.

35. Huo-Yuan Duan, GuO-Ping Liang Nonconforming elements in least-squares mixed finite element methods // Mathematic of computation. Vol. 73, Num. 245, pp. 1- 18, article electronically published on March 27, 2003.

36. J. van Kan A second order accurate pressure - correction scheme for viscous incompressible flow // SIAM J. Sci. Stat. Сотр. - 1986. - Vol. 7(3). - pp. 870-891.

37. Kelly D. W., Nakazawa S., Zienkiewicz О. C. and Heinrich J. C. A note on upwinding and anisotropic balancing dissipation in finite element approximations to convective diffusion problems //Int. J. Numer. Methods Eng. 1980, Vol. 15. 1705-1711.

38. Kim J., moin P. Application of a fractional-step method to incompressible Navier Stokes equations // J. Comput. Phys. 1985. - Vol. 59. P. 308-323.

39. KHAN L. A., Liu P. L. F. An operator splitting algorithm for the three-dimensional advection-diffusion equation // Int. J. Numer. Methods Fluids. - 1998. - Vol. 28. - P. 461-476.

40. KUPFERMAN R. A numerical study of the axisymmetric Couette Taylor problem using a fast high-resolution second-order central scheme // SIAM J. Sci. Comput. - 1998. - Vol. 20 (3). - pp. 858-877.

41. LEONARD B. P. A. Stable and accurate convective modeling procedure based on quadratic upstream interpolation // Comput. Methods Appl. Mech. Engrg. 1979. - Vol.19, N.l. - P. 59-98.

42. Liu X. D., Osher S., Chan T. Weighted essentially non-oscillatory schemes // J. Comput. Phys. 1994. - 115, No 2 - pp. 200-212.

43. Rida S., McKenty F., Meng F. L., Reggio M. A. A staggered control volume scheme for unstructured triangular grids // Int. J. Numer. Methods Fluids. 1995. - Vol. 25. - P. 697-717.

44. ROE P. L. Approximate Riemann Solvers, parameter vectors and difference schemes // J. Сотр. Phys. 1981. Vol.43. P. 357-372.

45. Shih Т. M., Tan С. H. and Hwang В. C. Effects of grid staggering on numerical schemes. // Int. J. Numer. Methods Fluids. 1989- Vol. 9. -P. 193-212.

46. VUIK C. Fast iterative solvers for the discretized incompressible Navier -Stokes equations // Int. J. Numer. Methods Fluids. 1996. - Vol. 22. -P. 195-210.

47. Warburton Т., Pavarino L. F., Hesthaven J. S. A pseudo-spectral scheme for the incompressible Navier-Stokes equations using unstructured nodal elements // J. Comput. Phys. 2000. - Vol. 164. - P. 1-21.

48. WARSI Z.U.A. Fluid dynamics. Theoretical and computational approaches. 3rd ed. Boca Raton, FL: Taylor к Francis, CRC Press, 2005.

49. Younes A., Mose R., Ackerer P., and Chavent G. A new formulation of the mixed finite element method for solving elliptic and parabolic PDE with triangular elements //J. Comput. Phys. 1999. -Vol. 149. - P. 148-167.

50. ZHANG S. A new family of stable mixed finite elements for the 3D Stokes equations // Mathematic of computation. article electronically published on August 31, 2004.130. http://www.itc.nsc.ru/linpar131. http://www.hpfem.jku.at/netgen