Численное моделирование гидродинамики весомой жидкости с помощью программно-исследовательского комплекса Текст научной статьи по специальности «Физика»

Уфа :УГАТУ, 2007 уГАТу • ^раёлШМ, и И т. 9, №5 (23). C. 108-113

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ

А. А. РОЗЕНМАН, С. С. ПОРЕЧНЫЙ

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ГИДРОДИНАМИКИ ВЕСОМОЙ

ЖИДКОСТИ С ПОМОЩЬЮ ПРОГРАММНО-ИССЛЕДОВАТЕЛЬСКОГО КОМПЛЕКСА

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

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

Аналогичная задача для невесомой жидкости решена Уайтхедом [1].

Численное исследование задачи обтекания весомой жидкостью диполя, расположенного на дне, проведено в [2].

Задача решается в системе координат, связанной с диполем. Пусть диполь с моментом N расположен в точке А (рис.1, а). Ось диполя ориентирована горизонтально. Скорость на бесконечности направлена в положительном направлении оси . В точках и , координаты которых зависят от параметров течения, скорость равна нулю. При таким течением может моделироваться бесциркуляционное обтекание тела, если его размеры малы по сравнению с расстоянием до дна АВ и до свободной поверхности АС. Расстояние между точками и (которые при N<0 располагается слева и справа от точки А) можно считать размером обтекаемого тела.

точки у уравнения Бернулли: 2 у

Fr"2h

= const, Fr =

Vo

Vgh'

(1)

На других участках границы имеют место следующие краевые условия: Яе г = 0 на АВ и АС, Иег = 0 на ВБ, где г = х + ¡у — комплексные координаты точки, А — точка расположения диполя.

В

Q

©

D

1. ПОСТАНОВКА ЗАДАЧИ

1.1. Граничные условия

Пусть скорость жидкости на бесконечности — , асимптотическая толщина струи — , сила тяжести действует вертикально вниз. Модуль вектора скорости жидкости V на свободной поверхности связан с высотой

с Ч2У

в ! F- А 9

Рис. 1. Формы области течения на различных плоскостях: а — физической плоскости; б — плоскости комплексного потенциала; в — плоскости изменения параметрического переменного

а

б

В связи с принятыми допущениями об идеальности среды течение соленоидально и потенциально, тогда решение задачи можно получить в виде аналитических функций комплексного переменного г(0, «КО, где £ — параметрическая переменная, областью изменения которой является полукруг (рис. 1,в), т — комплексный потенциал [3].

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

Область течения на плоскости го для данного течения приведена на рис. 1,6, зависимость го(0 определяется формулой [1,2]:

2hV0 1 + С W (С) = - 1п

7Г

1-е

Ni 2тг

1

1

1

гр

С + ip Р2 (С - >/р) Р2 (С + >/р)

(2)

где го (С) = ^ 1п ~~ функция, отображающая полуокружность на полосу, ^(С) — линейная комбинация степенных функций, коэффициенты которой подбираются так, чтобы выполнялось условие уравнения Бернул-ли (1), ^2(0 и 2з(0 — функции, введенные для учета особенностей решения при и в точке излома свободной поверхности.

Функция г1(С) = г(С)Д - «о(0 - МО -

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

^ (о = c,2m+iC

2m+l

m=0

Функции Z2(0 и ^з(0 определяются аналогично [2]

Ni =

n

clz

УсК

(3)

где — дипольный момент, — образ точки расположения диполя на плоскости . Производная (2) равна:

dw

rfc

2hV0

7Г

П1 ' 2

C-p

c2 ,2

2 '2

Cp

(C2 +p2)2 ' {p2c2 + 1)'¿

(4)

где ni = Ni/hVo.

1.2. Метод прямого конформного отображения

Функция z(0 является аналитической в полукруге , , нечетной (в

силу симметрии течения относительно оси у) функцией, имеющей нулевые значения мнимой части на горизонтальном диаметре Im£ = 0, 0<Re£<1 и удовлетворяющей краевому условию (1) на полуокружности ABC (при этом V = у = Imz).

При решении с помощью конформных отображений функцию z(0 будем искать в виде суммы степенного ряда и некоторых функций, учитывающих заданные особенности в точках

Z(0 = MZO(0 + Zl(0 + *2(0 + Z3(0].

(5)

zi(t) = Al

Z2( О = iA-2

1-$

1 + ií

1 + í

(6)

(7)

ü < /3 < 1.

где a — первый корень трансцендентного уравнения Стокса

7Г 7Г 1

a2Ctg02 = F^

ü < а < 1.

(8)

Наличие степенной особенности с показателем 2/3 означает, что свободная граница в точке излома образует угол, равный 2 , что соответствует известной особенности волны Стокса. Это значение получено исходя из предположения, что предельное значение угла при подходе к точке существует и не равно нулю. То есть при нулевом предельном значении можно ожидать существование других режимов с критической точкой .

Здесь также учитывалось, что, как и г(0, функции и должны быть нечетными и иметь действительные значения на действительной оси.

Тем самым определен вид функции (5).

1.3. Обтекание диполя невесомой жидкостью

Функция Жуковского ш(С) = =

( — угол наклона вектора скорости к оси абсцисс, т для данного те-

чения строится аналогично задаче, описанной в [1]. Для аналитической функции имеют место следующие краевые условия:

ГОна В А

0 на

в = 0 на РС, т = 0 на CD, (9) 7г на АК — 7Г на АР\,

где — точка расположения диполя, и — критические точки. Такие условия имеют место при щ > 0. Этим условиям удовлетворяет функция

00 = 'I 1п

(р2С2 + I)2 (С2 + д2) (с2 + С>

(С2

■Р2)2(?2с2

1) (чК2 + 1)

(10)

где — образ точки расположения диполя, и — образы точек и . Параметры и определяются из решения уравнения = 0. В частном случае 0 (диполь расположен на дне) [2]:

00 = 'I 1п

(С2 + д2) С2(д2С2 + 1У

дм

1к

2/гТ'о

7Г

1

5__^и _1

С2 4 ^ с2

(11)

(12)

В формуле (12) момент берется с половинным весом, так как в частном случае =0 на плоскости £ совмещаются два диполя, а их суммарный момент связывается с моментом диполя на физической плоскости.

Для связи дипольных моментов на физической и параметрической плоскостях необходимо найти производную йг/йС, в точке £=0:

с1г

1 , (1(1)

сК и/ ПС Подставим (11) и (12) в (13):

с1г

Ж

2к С2 (^С2 + 1)

7Г

С2

П\

Ь1 ^ с2

(13)

1

Переходя к пределу при , получим:

с1г 2/г 1 п\

< ' =

Используя (14) в (3) найдем уравнение, связывающее и ,

щ 2

п\

2,/1

П1 2

7Г

(15)

которое можно решить численно.

В случае обтекания диполя весомой жидкостью для пересчета дипольных моментов по формуле (3) необходимо вычислить производную / в точке :

^ (¿р) = Ь,. (К \тг{1+р2)

-^((1-р)

-1/3

(1 + р)1'3) + (2т + 1) С2п,+1 (-1)т р2т

т=0

-2А1«\/1 + р2 соя (а — 1) aIctgp) . (16)

2. ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ 2.1. Численный алгоритм

Задача сводится к решению уравнения (1), решать ее будем методом коллокаций [2]. В бесконечной сумме сохраняется конечное число слагаемых, а равенство (1) удовлетворяется в дискретных точках ат = = 7гт/(2Ж), т = Полученные N

нелинейных уравнений образуют систему, которая решается методом Ньютона с регулированием шага относительно параметров С2то+1; т = 0, N — 4, ^-Рг. Для волн типа Стокса является искомым параметром.

Оценка погрешности производится по усовершенствованному правилу Рунге [4] путем сравнения значений параметров (например, числа , координаты точки и др.), полученных при последовательном возрастании N, а также по максимальной невязке Дтах уравнения (1), рассчитанной в промежуточных точках между узлами коллокаций .

2.2. Численные результаты

На рис. 2 представлены зависимости числа, обратного числу Фруда, от высоты гребня волны для разных высот расположения диполя. Кривые начинаются с решения, соответствующего обтеканию невесомой жидкостью (Рг = оо), и заканчиваются волной предельной высоты (типа Стокса) с критической точкой на гребне и изломом свободной поверхности (с внутренним углом 120°).

Для волн с критической точкой на гребне справедливо равенство, следующее из уравнения Бернулли

Рг2

Ус = 1 + — •

ш

1.2 1.4 1.6 1.8 Ус

Рис.2. Зависимости величины 1/Рг от высоты

гребня ус* для разных высот расположения диполя: 1 — ул = 0,05; 2 - ул = 0,5; 3 - ул = 0,75;

4 - уА = 0,9

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

При обтекании диполя такое явление не может иметь места, так как отсутствует как горизонтальная, так и вертикальная составляющая силы взаимодействия диполя с набегающим потоком. Исследование показало, что приближение диполя к свободной поверхности при малых ограничено возникновением критического режима, при котором число Фруда Рг=1, где Т^, Л, — скорость и ордината точки свободной поверхности на бесконечности.

На рис. 3 представлены формы свободной поверхности волн типа Стокса при различных расстояниях диполя от дна у.д. для п = — 1. Верхняя кривая соответствует диполю, расположенному на дне, нижняя кривая - Рг =

.

На рис. 4 даны зависимости от ординаты ус при п=еош^0. При п ^ — 1 кривые начинаются с предельного значения ( ) и заканчиваются при .

У 2

1

Рис. 3. Формы свободной поверхности волн Стокса при различных высотах расположения диполя (п = — 1)

Ул 1.0

0.5

Рис. 4. Зависимости высоты расположения диполя от высоты гребня (для волн типа Стокса): 1 — п = --0,25; 2 — п = —0,5;

3 — п = —1,4 — п= —1,5

Для предельный режим

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

3. ОПИСАНИЕ ПРОГРАММНО-ИССЛЕДОВАТЕЛЬСКОГО КОМПЛЕКСА ДЛЯ РЕШЕНИЯ ЗАДАЧ ВЫЧИСЛИТЕЛЬНОЙ МЕХАНИКИ

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

Была сделана попытка разработать программный комплекс, позволяющий сделать

< А > /

г /

1 .5 1 .75 2.! Ус

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

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

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

Закрытая система

Вычисления

Исходные данные

1 ^Рзультаты^

^^Интерфейс^^

Рис. 5. Структура закрытой системы

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

да и порожденного каждым программным кодом данных (см. рис. 6).

Задача 1.

^^Трансляция^^

^^Результаты^^

Задача N.

К

Вычислитель

С

Данные

Рис. 6. Состав базы данных результатов решения задачи

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

Рис. 7. Структура исследовательского комплекса

Численный эксперимент разбивается на несколько этапов: построение дискретного аналога дифференциальной задачи, отладка и тестирование пользовательских программ численного решения дискретной задачи на серии тестовых задач, проведение расчетов.

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

Проведение расчетов делится на следующие этапы: постановка задачи численного исследования, подготовка начальных данных, а также корректировка приложения-решателя

для проведения данного исследования; трансляция приложения-решателя; запуск приложения-решателя; анализ полученных данных.

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

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

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

Данные в каждом проекте делятся на три типа: файлы с исходными данными; файлы с программным кодом, исходные тексты программы решателя; результаты — данные, полученные после прогона исходных данных через оттранслированную программу-решатель.

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

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

ВЫВОДЫ

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

Наиболее интересным является два типа найденных в работе предельных конфигураций при приближении диполя к свободной поверхности: при предельный режим

, в противном случае — выход критических точек на свободную поверхность.

СПИСОК ЛИТЕРАТУРЫ

1. Whitehead, L. G. Two-dimensional windtunnel interference / L. G. Whitehead // ARC Rep. and Memo. 1950. № 2802.

2. Житников, В. П. Обтекание диполя под поверхностью весомой жидкости с образованием солитона / В. П. Житников, Н. М. Шерых-алина, О. И. Шерыхалин // Динамика сплошной среды. Новосибирск : Ин-т гидродинамики, 1999. Вып. 114. С. 31-34.

3. Гуревич, М. И. Теория струй идеальной жидкости / М. И. Гуревич. М.: Наука. 1979. 536 с.

4. Житников, В. П. Методы верификации математических моделей в условиях неопределенности / В. П. Житников, Н. М. Шерыхалина // Вестник УГАТУ. 2000. № 2. С. 53-60.

ОБ АВТОРАХ

Розенман Алексей Александрович, асп. каф. компьют. матем. Дипл. инж. (УГАТУ, 2000). Готовит дис. в обл. прогр. комплексов.

Поречный Сергей Сергеевич, асп. каф. компьют. матем. Дипл. инж. и магистр по инфор. и прогр. обесп. САПР (УГАТУ, 2006). Готовит дис. в обл. мат. моделирования физ. процессов.