Применение квазигидродинамических уравнений для математического моделирования теченийвязкой несжимаемой жидкости

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

им. М.В. ЛОМОНОСОВА

ФАКУЛЬТЕТ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И КИБЕРНЕТИКИ

Кафедра вычислительных методов

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

Жериков Андрей Валерьевич

ПРИМЕНЕНИЕ КВАЗИГИДРОДИНАМИЧЕСКИХ УРАВНЕНИЙ ДЛЯ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ТЕЧЕНИЙ
ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ

05.13.18 – "Математические моделирование, численные методы и комплексы программ"

АВТОРЕФЕРАТ

диссертации на соискание ученой степени

кандидата физико-математических наук.

Москва, 2009

Работа выполнена на кафедре вычислительных методов факультета вычислительной математики и кибернетики Московского государственного университета имени
М.В. Ломоносова.

Научный руководитель:	старший научный сотрудник кандидат физико-математических наук Калачинская Ирина Станиславовна
Официальные оппоненты:	доктор физико-математических наук доцент Попов Виктор Юрьевич старший научный сотрудник кандидат физико-математических наук Чурбанов Александр Георгиевич
Ведущая организация:	Институт Математического Моделирования РАН

Защита состоится « » 2009 г. в часов минут на заседании диссертационного совета Д 501.001.43 при Московском государственном университете имени М.В. Ломоносова по адресу: 119991, Москва, Ленинские горы, МГУ, второй учебный корпус, факультет ВМК, ауд. 685.

С диссертацией можно ознакомиться в научной библиотеке факультета ВМК Московского государственного университета имени М.В. Ломоносова.

Автореферат разослан « » _____________2009 г.

Ученый секретарь диссертационного совета Д 501.001.43, доктор физико-математических наук, профессор

Захаров Е.В.

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

Актуальность

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

Нетрадиционным подходом к построению алгоритмов расчета вязких течений является использование квазигидродинамических (КГД) уравнений [1–4], отличающихся от уравнений Навье-Стокса дополнительными слагаемыми дивергентного вида с малым параметром в качестве коэффициента. При построении разностных схем, аппроксимирующих данную систему уравнений, эти дополнительные слагаемые проявляют себя как эффективные регуляризаторы, обеспечивая устойчивость численных алгоритмов. Тем самым КГД система может рассматриваться как система уравнений Навье-Стокса с регуляризатором специального вида. Ранее на ряде примеров было показано, что применение КГД уравнений открывает широкие возможности для построения новых эффективных численных методов для расчета нестационарных течений жидкости в переменных «скорость – давление».

Цель работы

Основные цели работы состоят в построении численных алгоритмов для расчета течений вязкой несжимаемой жидкости, основанных на использовании КГД уравнений, а именно:

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

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

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

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

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

Научная и практическая значимость работы

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

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

Апробация работы

Основные результаты работы докладывались и обсуждались:

• на V Международной Конференции по Математическому Моделированию. 30 сентября – 6 октября 2002 года, г. Дубна.
• на XIV Международной Конференции «Математика. Компьютер. Образование». 22 ­– 27 января 2007 года, г. Пущино.
• на XV Международной Конференции «Математика. Компьютер. Образование». 28 января ­– 2 февраля 2008 года, г. Дубна.
• на научном семинаре в лаборатории математического моделирования в физике кафедры вычислительных методов, факультет вычислительной математики и кибернетики, Московский государственный университет имени М.В. Ломоносова. 25 февраля 2009 года.
• на заседании кафедры вычислительных методов, факультет вычислительной математики и кибернетики, Московский государственный университет имени М.В. Ломоносова. 13 мая
  2009 года.
• на семинаре в Институте математического моделирования РАН.
  1 октября 2009 года.

Публикации

Основные результаты диссертации опубликованы в 9 работах. Список публикаций приведен в конце автореферата.

Объем и структура диссертации

Диссертация состоит из введения, 5 глав, заключения и списка литературы. Текст изложен на 87 страницах, диссертация содержит 9 таблиц, 66 рисунков. Список литературы составляет 99 наименований.

Содержание диссертации

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

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

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

КГД система в приближении Буссинеска может быть представлена в виде

(1)

(2)

(3)

При этом величины и вычисляются по формулам

(4)

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

Параметр может быть вычислен по формуле , где cs – скорость звука в жидкости при температуре T0. величина интерпретируется как пространственно-временной средний импульс единицы объема.

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

В записи системы (1) (4) использованы стандартные обозначения тензорного анализа. При вычислении дивергенции от несимметричного тензора П свертка осуществляется по его первому индексу.

Во второй главе проводится моделирование термокапиллярной конвекции полупроводникового расплава. Такая задача возникает при получении кристаллов методом безтигельной зонной плавки в условиях невесомости, когда конвективное движение расплава определяется процессами термокапиллярной конвекции, или конвекцией Марангони.

Ламинарные режимы конвекции Марангони хорошо исследованы экспериментально и теоретически.

В первом параграфе рассматривается конвекция в прямоугольной каверне. Исследуется влияние электромагнитного поля различной интенсивности и направления. В безразмерном виде интенсивность магнитного поля представляется числом Гартмана, которое вычисляется по формуле , где L – размер области, H0 – напряженность магнитного поля, с – скорость света, – электропроводность жидкости и – динамическая вязкость.

Также в первом параграфе рассматривается влияние шага расчетной сетки на результаты моделирования.

Рис. 1. Линии тока для течения под действием
магнитного поля в квадратной каверне
для чисел Гартмана Ha=0, 50 и 100 соответственно

На рис. 1 изображены линии тока установившегося течения полупроводникового расплава при отсутствии магнитного поля, т.е. Ha=0, (слева), На=50 (в центре) и при Ha=100 (справа). Расчеты проводились для числа Марангони Ma=1000 и числа Прандтля Pr=0.018.

Во втором параграфе проводится исследование термокапиллярной конвекции в цилиндрическом образце. Моделирование проводилось в два этапа: сначала моделировалось нестационарное течение при отсутствии магнитного поля (то есть при числе Гартмана Ha=0) и затем после установления колебательного режима включалось магнитное поле.

Рис. 2. Линии тока для течения в цилиндрической области
для чисел Гартмана Ha=0 и 50 соответственно

На рис. 2 изображены линии тока в установившемся колебательном режиме течения полупроводникового расплава при отсутствии магнитного поля, т.е. Ha=0, (слева) и при Ha=50 (справа). Расчеты проводились для числа Марангони Ma=10000 и числа Прандтля Pr=0.018.

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

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

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

Рис. 3. Треугольная сетка

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

Производные по времени аппроксимируются с помощью метода Эйлера. Частные пространственные производные определяются на основе производных по направлению или с использованием формулы Грина (параграф 3.3).

В параграфе 3.4 приведен алгоритм построения фиктивных узлов и построены аппроксимации граничных условий. В параграфе 3.5 описана аппроксимация уравнения Пуассона для давления, при которой уравнение Пуассона интегрировалось по контрольному объему, а частные производные определялись с помощью формулы Грина. В результате чего была получена система линейных алгебраических уравнений для давления с симметричной матрицей. Эта система уравнений решалась как итерационно методом наискорейшего спуска, так и методом Холецкого.

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

На рис. 4 приведено распределение скорости в задаче обтекания кругового цилиндра для числа Рейнольдса Re=150 (число Рейнольдса рассчитывается по диаметру цилиндра). В данной задаче жидкость втекает слева с постоянной скоростью и свободно вытекает справа.

Рис. 4. Дорожка Кармана в задаче обтекания цилиндра

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

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

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

Рис. 5. Линии тока для нестационарного течения в канале

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

Исследовалась структура течения при числах Рейнольдса Re=10, 100 и 500. Результаты моделирования показали, что при числах Рейнольдса Re=10 и 100 течение является стационарным, а при Re=500 – нестационарным. На рис. 5 приведены линии тока для нестационарного течения при Re=500. Проведено сравнение распределения осредненной скорости в вертикальном сечении в двух метрах после уступа с результатами натурного эксперимента. Сравнение показало совпадение профилей скорости, полученных в результате моделирования, с результатами натурного эксперимента с точностью 10%.

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

Основные результаты