Реализация метода обмена интенсивностями вортонов-отрезков для учета вязкости в методе вихревых элементов Текст научной статьи по специальности «Физика»

УДК 532.5.013 DOI: 10.18698/1812-3368-2018-3-48-67

РЕАЛИЗАЦИЯ МЕТОДА ОБМЕНА ИНТЕНСИВНОСТЯМИ ВОРТОНОВ-ОТРЕЗКОВ ДЛЯ УЧЕТА ВЯЗКОСТИ В МЕТОДЕ ВИХРЕВЫХ ЭЛЕМЕНТОВ

О.С. Коцур oskotsur@gmail.com

Г.А. Щеглов shcheglov_ga@bmstu.ru

МГТУ им. Н.Э. Баумана, Москва, Российская Федерация

Аннотация

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

Ключевые слова

Метод вихревых элементов, метод обмена интенсивностями, вортоны-отрезки, вязкость, завихренность, уравнение эволюции завихренности

Поступила в редакцию 14.04.2017 © МГТУ им. Н.Э. Баумана, 2018

Работа выполнена при поддержке гранта РФФИ (проект № 17-08-01468 А)

Введение. Широко известны сеточные методы вычислительной гидродинамики, например методы контрольных объемов [1], конечных элементов и конечных разностей, которые являются стандартом де-факто в сфере инженерных расчетов. Метод вихревых элементов (МВЭ), не требуя пространственной сетки для расчета, представляет собой альтернативу сеточным методам в большом классе задач в несжимаемой постановке [2, 3], где использование сетки нежелательно либо неэффективно, например, для задач взаимодействия упругого тела и жидкости (Fluid-Structure Interaction, FSI) в безграничной области. Использование деформируемой сетки в задачах FSI значительно замедляет расчет, по-

скольку кроме уравнений гидродинамики также решаются уравнения деформации сетки [4].

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

V-V = 0; (1)

—+ (V-V)© = (ffl-V)V + vA© (2)

c начальными и граничными условиями ©(x,t0) = ©0(x), lim V(x,t) = V,

lim ©(x, t) = 0, где V=V (x, t), © =VxV — поля скорости и завихренности,

заданные в R3, t е[0,Т]; ©0 (x) — распределение завихренности в начальный момент времени t0 = 0; V — скорость на бесконечности; v — кинематическая вязкость.

Классический МВЭ можно рассматривать как обобщение метода частиц [5] на уравнения движения невязкой несжимаемой жидкости, записанные в виде (1), (2) при v = 0. Он подходит для решения инженерных задач с большим числом Рейнольдса. В таких задачах начальная завихренность движется по траекториям жидких частиц. В случае рассмотрения задачи расчета обтекания тел генерация завихренности вблизи стенки учитывается в соответствии с подходом Прандтля о разделении течения на пристеночную область, где влияние вязкости существенно, и область внешнего течения, где вязкость не оказывает значительного влияния на параметры потока.

Попытки учесть вязкость в области течения в МВЭ привели к появлению множества подходов [2]. Одним из них является «расщепление» уравнения (2) по физическим процессам (splitting algorithm) [6]. При этом за полшага по времени происходит только перенос завихренности (v = 0), а за следующие полшага — поправка поля завихренности за счет диффузии. Диффузия моделируется приданием частицам (вихревым элементам, ВЭ) «случайных блужданий» по типу броуновского движения (random-walk) [6].

Известны также подходы, сохраняющие детерминированный характер метода частиц. В методе диффузионных скоростей (Diffusion Velocity Method, DVM) для вязких течений вихревые частицы движутся не по траекториям жидких частиц, а по траекториям, определяемым так называемой диффузионной скоростью (иногда называемой скоростью Фридмана [7]), сохраняя свою интенсивность [8, 9]. Метод DVM основан на представлении уравнения эволюции завихренности (2) в форме, аналогичной форме уравнения конвекции, которое эффективно решается методом частиц. В явном виде такое преобразование можно сделать лишь в двумерной [9], либо осесимметричной [7], постановке. В трехмерной постановке в явном виде выделить диффузионную скорость не удается.

Более подходящим для трехмерного случая является метод обмена интен-сивностями (Particle Strength Exchange, PSE) [2, 10]. Здесь в отличие от метода DVM вихревые частицы движутся вместе с жидкими частицами по полю скоростей, как и в невязком случае. При этом за счет действия диффузии происходит перераспределение интенсивностей близкорасположенных частиц между собой.

Существуют также попытки объединения методов DVM и PSE [11]. При этом из диффузионного члена уДю выделяется в некотором смысле «максимальная» часть, которая может быть преобразована к слагаемому, содержащему аналог диффузионной скорости. Вклад оставшейся части диффузионного члена учитывается в рамках PSE.

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

Интегральная аппроксимация диффузионного члена. Метод PSE сложился как методика интегральной аппроксимации диффузионного члена D (ф) при решении уравнений типа конвекции-диффузии в форме, удобной для применения метода частиц:

-J + (a-V)cp+аоф = D(q>), (3)

где ф = 9(x,t) — искомое поле; D(ф) = V-(Г Уф) — диффузионный оператор; a = a (x, t), a0 = a0(x, t) — заданные векторное и скалярное поля; Г — коэффициент диффузии, x eR3, t е[0,T]. Дифференциальный оператор D (ф) из (3) при постоянном коэффициенте диффузии Г = v = const примет вид v Дф.

Назовем ядром порядка k скалярную функцию r) скалярного аргумента, удовлетворяющую следующим моментным уравнениям для k > 1, k е N:

да 3

J r 4ц(г )dr = —;

0 2%

J гр+2ц(г )dr = 0; (4)

0

да

J rk+4ц(г)dr < да,

о

где 2< р < к +1, р — четное.

В качестве примера ядра порядка к = 2 можно привести гауссову функцию вида

x2

1 - —

^(х) = —¡=-e 2

Обзор других вариантов ядер различных порядков, удовлетворяющих условиям (4), можно найти в работах [10, 13]. Обозначим через ^(х), (х) сферически симметричные функции, заданные на К3 так, что г|(х) = х|);

^8(х) = —3, где 8 — положительное вещественное число, характеризую-

3

в3 \г

щее радиус шара, в котором сосредоточены значения (х). Величину в называют характерным размером ядра.

В работе [10] показано, что в некоторой норме интегральный оператор

Qs(cp(x,t)) = 4 J Пв(х-y)(ф(у,t)-ф(х,t))dy (5)

,2

R

S R3

аппроксимирует диффузионный оператор D (ф) с точностью О (е*) при выборе ядра порядка к.

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

Метод вихревых элементов. Задача о движении невязкой несжимаемой жидкости описывается системой (1), (2), где вязкость V принимается равной нулю. С помощью МВЭ поле завихренности представляется в виде суперпозиции полей отдельных вихревых элементов:

N

®(х, г )=X©* (х - х * (г)),

к=1

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

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

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

Классический МВЭ основан на представлении функций ВЭ в виде дельта-функций

(х, *)=У к (*) 5( х - хк (г)), (6)

характеризующих «концентрацию» завихренности в точке x к (£), в которой находится ВЭ в момент времени £. Весовой векторный коэффициент у к (£) определяет ориентацию и «интенсивность» ВЭ. Такой ВЭ называется сингулярным точечным вортоном. Подобное представление поля завихренности в виде суперпозиции дельта-функций и составляет суть классического метода частиц [5]. Подробный теоретический анализ классического МВЭ, условия сходимости решений, полученных МВЭ, к решению уравнений Эйлера для невязкой жидкости приведены в работах [2, 15].

Использование дельта-функций для аппроксимации поля завихренности некорректно вследствие нарушения свойства соленоидальности базисного поля [14]. Точечный вортон с учетом добавочной завихренности известен как вортон Новикова [16].

Однако для любого дифференцируемого убывающего на бесконечности векторного поля существует возможность дополнения его до соленоидального с помощью теоремы Гельмгольца [17]:

О, = Па +П'.

Здесь 0,а — исходное векторное поле; О, — соленоидальное векторное поле, построенное на базе 0,а; О' — добавочное поле, представляющее собой градиент некоторой скалярной функции и дополняющее исходное поле Оа до соленоидального. Если в качестве поля Оа использовать поле, формируемое функцией ак то добавочное поле О' называется добавочной завихренностью ВЭ.

В настоящей работе рассмотрена формулировка МВЭ, основанная на представлении поля завихренности в виде сингулярных вортонов-отрезков (ВО) [18]. Функция ВЭ к-го ВО получается интегрированием сингулярной функции (6) вдоль заданного в пространстве отрезка 2^ (рис. 1), который в процессе движения может изменять длину и направление:

1

Шк ^, £ ) = у к (£) /б^к (£) + 5 h к (£ )))&,

-1

где xк (£) — радиус-вектор середины отрезка 2^(£), который является маркером ВЭ.

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

Поле завихренности, аппроксимируемое с помощью ВО, имеет вид

N 1

^ ) = 2 У к (О /З^к ^ ) + 5 hk ^ )))&. (7)

к=1 -1

Вектор у к удобно представить в виде у к ^ ) = у к ^ ) h к ^ ), где — коэффициент пропорциональности между векторами, который назовем интенсивностью ВЭ [18].

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

Введение ВО позволяет в лагранжевой постановке осуществить переход от уравнений в частных производных (1), (2) к системе обыкновенных дифференциальных уравнений относительно координат маркеров x к, векторов hk и ин-тенсивностей у к [18]:

dx k dt dhk dt

d jt dt

= V (x k, t);

= [B(xt, t)]• ht; (8)

= 0.

Из первого уравнения (8) видно, что маркеры ВЭ перемещаются по траекториям жидких частиц со скоростью V ( xк, t). Поле скоростей V (х, t) восстанавливается из поля завихренности ю( x, t) с помощью интеграла Био — Сава-ра. Поле скорости, индуцируемое одним вортоном, имеет вид [18]:

Vk (x,t) = Уk (hk Xx

4% h k x x k

( s s ^

S2___si

LV ls1|V

• hk

где s1, S2 — векторы, соединяющие концы вортона с точкой x (см. рис. 1). Общее поле скорости представляет собой суперпозицию полей всех вортонов:

N

V(х, t) =X Ук (х -xк ^)).

к=1

Конкретный вид матрицы Якоби [В (хк, I)], соответствующей в (8) осред-ненному по ВЭ градиенту поля скоростей жидкости V (х, t) в точке хк зависит от типа используемого ВЭ [18, 19].

Интенсивности ук в (8) постоянны, поскольку интенсивность вихревых трубок (поток вектора завихренности через сечение трубки) в невязких течениях не меняется.

Определение начальных условий для системы (8), векторных величин Уко = Ук (to) (интенсивности ук0 и вектора-отрезка Ьк0), соответствующих полю ю0 в начальный момент времени, может быть осуществлено следующим образом. Пусть завихренность ©0 отлична от нуля в некоторой ограниченной области О.. Разобьем область О. на N достаточно малых непересекающихся

N

подобластей О,, так, что О. = и О;. Поместим в центры полученных подобла-

г=1

стей х к 0 по одному ВО. Интегрируя (7) для ю0 по каждой подобласти 0,к, получаем

N 1

I ©0 (хх = | X I Yq05(х-хq0 -=

Ок Qkq=1 -1

1 1 = 11 У к 05 (х - х к 0 - Ьк 0 5 = {у к 0^5 = 2ук0.

-1 ак -1

Таким образом, начальные коэффициенты ук0 должны назначаться по правилу

Начальную длину ВО | 1 можно выбрать равной характерному размеру области . Например, при дискретизации области кубической сеткой с ребром к можно принять |Ь^0| = к /2. Тогда направление Ью и величина ук0 однозначно определяются из (9).

Учет вязкости с помощью метода обмена интенсивностями. Рассмотрим модель вязкой несжимаемой жидкости. Запишем уравнение эволюции завихренности (2) через субстанциональную производную В / В1 и аппроксимируем диффузионный член V Да с помощью интегрального оператора 0е (ф) (5), записанного для ю:

Дискретизация (10), (11) с помощью точечных вортонов приведена в работах [1, 3]. Получим дискретизацию уравнения (10) с использованием ВО. В качестве допущения не будем учитывать добавочную завихренность, полагая ее влияние малым. Такое допущение в точности выполняется, например, при моделировании замкнутых вихревых структур (вихревых колец, рамок), аппроксимируемых ВО одинаковой интенсивности, поскольку в этом случае добавочной завихренности не возникает. Это допущение применимо и в случае, когда в области течения находится множество ВЭ близкой по модулю интенсивности. Однако дать точную оценку вносимой погрешности, обусловленной влиянием добавочной завихренности, представляется затруднительным.

Далее для краткости в выкладках опустим обозначение зависимости различных полей от времени.

Подставляя в (11) вместо ©(у) аппроксимацию поля завихренности (7), получаем

уk0 =уkohk0 = - J Юо (x)dx.

2

(9)

-= ©VV + QS(©

Dt v

(10)

Qs(©)= ^ J Пв(х"У)(®(y,t)-®(x,t))dy.

S R3

(11)

- 1 -| Л

= - У q ¡Це(х - X q - hqS )ds -©(x ) J (x - y) dy .

Л

R3

/

Интеграл J ^s(x - y)dy возьмем приближенно, используя квадратурную

R3

формулу

N

J " y )dy » - X q )oq,

r3 q=1

(12)

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

Подставим теперь (7) вместо ю( х) и проинтегрируем О8 (га) по области Ок, учитывая %(х - хч ) = %(хч -х):

J Qs(ro)dx=

"k

( N

z

v q=1

Yq j ^s(x k - X q ~ hqS )ds

-1

N

-z

q=i

Yk<5q j ^s(xq ~ Xk ~ hkS) ds

-1

(13)

Здесь | % (х—хц — Ьд5) dx — интеграл, также взятый приближенно, учитывая

малость области интегрирования 0,к аналогично (12):

| цЕ (х -хц - Ьц5) dx « Ск% (хк - хц -Ьц5).

"к

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

1 1

Gij=2 1 (xi " x j - hjS )ds.

2 -1

(14)

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

С учетом (14) формула обмена интенсивностями (13) принимает простой симметричный вид

2\ м

| Оs((й)dх = — ^(УцСкОкц ~УкЪсряк). Ок £ 9=1

Имея интегральную аппроксимацию диффузионного оператора О8 (га), получаем модификацию МВЭ с вортонами-отрезками на основе метода РББ. Подставим (7) в (10):

( N 1 ^ N

D_ Dt

X \yk5(x - xk - hkS)ds

\ q=1 -1 у

=Z

q=1

|5(x - xk - hkS) Yk • VVds

.-1

-Qs(ö>).(15)

Субстанциональную производную в левой части можно внести под знаки суммы и интеграла. Для к-го слагаемого в левой части (15) это позволит получить

D N 1 1 dvk

Z J Yk ô(x - xk - hks )ds = J —d— 8(x - x k - hkS )ds.

q=1 -i J -i dt

Далее проинтегрируем левую и правую части (15) по объему Qk, содержащему вортон k, и поменяем порядок интегрирования:

1 N dy k

J Z J ~T~ s(x - xq - hqS )dx ds =

-1 q=1 Qk dt

1 N

= J s

-1 q=1

J ô(x - xq - hqS) Yq • VVdx

Qk

2v N

dS +— Z(ïqakGkq "Y k VqGqk ). (16)

s q=1

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

(16) исчезают. Учитывая, что | ——ds = 2-^-, получаем

_1 dt dt

—у >. 1 1 V м

~т=7 У к I УУ (х - X к - hkS )ds+— £(уч ЯкОщ-У к ). (17)

dt 2 _1 8 q=l

Аналитическое выражение для тензора градиента скорости УУ (х) для совокупности ВО получено в работе [18]. Однако интеграл УУ (х - хк - Ьks) по s предпочтительно брать приближенно.

Составим матрицу Якоби [Вк ] из компонент тензора градиента скорости УУ (х — хк — hkS). Переходя от тензорной нотации к матричной, отметим (у к-V) У = [Вк ]-у к.

Обозначим через [ В к ] осредненную по длине вортона матрицу [Вк ]:

_ 1 1

[Вк ] =-|[УУ (х -хк -hkS(18)

2 -1

Тогда первое слагаемое в правой части (17) можно записать как матричное произведение [Вк ]-у к.

Запишем получившиеся дифференциальные уравнения относительно параметров ВО:

dx k dt

=V (x k, t );

—у, _ V N

-р-=[Вк )ук +— £ (уч-уко^к). (19)

dt в ч=1

Получим теперь уравнение для изменения интенсивностей вихревых элементов у к. Левую часть (19) можно представить в виде

йу к йу к , йЬк

= -*— Ь +Ук——. (20)

йг йг йг

Для нахождения производной й^ требуется определить, как изменяется

йг

со временем длина Ьк, от которой также зависят коэффициенты Оц. Известно (например, [20]), что бесконечно малый материальный отрезок 51 изменяется по закону

—51 = (51 -У)У. йг v '

Полагая длины ВО Ьк достаточно малыми, можно принять, что

^ = (Ьк-У)У (хк), йг К 7

где черта означает осреднение по отрезку Ьк, либо в матричном виде с использованием обозначения (18)

= [Вк ]■ Ьк. (21)

йг -1

Подставляя (21) в (20), затем в (19) и отметив, что члены [Вк]-ук справа и слева взаимно уничтожаются, получаем

йу- Ь к = -у Е {lq Ь ^kОkq -у к Ь к GqОqk ). йг в ^=1

Полагая |Ьк| ^ 0, приведем последнее соотношение в скалярный вид. Для этого скалярно умножим правую и левую части на вектор Ьк/1 Ьк |2 и введем для удобства записи обозначение

Ь к ■ Ь я

= и и

Ьк ■ Ьк

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

=V (x k, t);

dx k dt

=[Bk ]• hk; (22)

dhk r-

йг

йук V Я; п п п \

—¡— = — 'L\Уq akCkqОkq ~У к ^qОqk ). йг В q=l

Система (22) отличается от системы (8), описывающей невязкое течение, только уравнением эволюции интенсивностей. Изменение интенсивностей ВО обусловливается диффузией завихренности, возникающей вследствие наличия члена уД© в уравнении (2). Необходимо отметить, что вязкость в (22) не влияет на перемещение ВО и деформацию его вектора. Это является следствием сделанных допущений о малости дополнительной завихренности.

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

( hm Л

рядок OI v m+i I, где h — размер ячейки сетки, на которой изначально расположены частицы; m — параметр, определяемый принадлежностью искомого поля f пространству Соболева Wm,w (R3). Таким образом, для снижения погрешности необходимо обеспечивать условие h < в, которое в практическом смысле можно интерпретировать как требование сохранения расстояния между частицами (вортонами) меньшим, чем выбранный параметр в для модели PSE. Есть основания полагать, что сделанный вывод остается справедливым и для рассмотренной модели МВЭ с вортонами-отрезками. Выбор длин ВО hk и объемов , приходящихся на один ВО, зависит от геометрии исходного поля завихренности, в частности от кривизны вихревых линий.

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

Исследуем модельную задачу о диффузии бесконечной прямолинейной вихревой трубки. Выделим в бесконечной трубке фрагмент, состоящий из нескольких слоев, составленных из ВО. Каждый слой представляет собой круг радиусом R, равномерно заполненный ВО длиной 2h, векторы которых направлены перпендикулярно плоскости круга (рис. 3). В подобной постановке, когда все ВО параллельны между собой, правая часть второго уравнения (22) [Bk ]• hk обнуляется, и векторы ВО hk не изменяют ни направление, ни длину (h k = h = const) .На каждый ВО в круге приходится одинаковый объем, так что Ok = о = const. Параллельность и равенство длин всех ВО обеспечивает равенства Ckq = 1, Gkq = Gqk .

Рис. 3. Конфигурация набора вортонов в слое (всего 1521 ВЭ в слое)

Движение такой конфигурации ВО описывается системой (22), которая принимает следующий вид:

dx k dt dhk dt

^ = "Г q-У k )Gqk ].

dt 8 q = 1

=V(xk,t); = 0;

(23)

N

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

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

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

V =-. При этом для моделирования трубки достаточно рассмотреть один

2%т

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

Уравнение эволюции завихренности для прямолинейной бесконечной вихревой трубки есть классическое уравнение диффузии [21]:

ÖL>3

ö>2

щ

cto ~dt

-=vAra.

(24)

COl

CO4

Рис. 4. Шаблон типа «крест»

Последнее уравнение в (23) можно рассматривать как разностную аппроксимацию уравнения (24) на сетке из ВЭ. Действительно, проведем разностную аппроксимацию второго порядка для оператора Лапласа Дю в точке х ^ на регулярном шаблоне типа «крест» (рис. 4) с одинаковым расстоянием е между точками [22]:

Щ -2^ +Ю3 <В2 -2(йк +Ю4

Аею k =-:-+-:-. (25)

в2 в2

Подставляя (25) в (24), имеем следующее разностное представление уравнения эволюции завихренности вихревой трубки:

——^-7 -ак), (26)

Ш В

где N = 4 — число соседних точек в шаблоне.

Третье уравнение в (23) и (26) имеют схожую структуру с точностью до коэффициентов. Коэффициенты в (23) определяются выбором ядра , параметра в и взаимным расположением ВО. Кроме того, в (23) суммируются все без исключения ВО, однако наибольшее влияние оказывают только ВО, близко расположенные к вортону к, поскольку ядро , а следовательно, и коэффициенты Gqk быстро убывают с увеличением расстояния.

Рассмотрим тонкую вихревую нить с циркуляцией Г в начальный момент. Завихренность в момент времени к в плоскости, перпендикулярной нити, имеет вид [21]

, Л Г ( r2 >

ю(г, t)=-exp--

4^vt 4vt

(27)

Циркуляция скорости вдоль окружности радиусом Я в момент времени к

Я

( f г>2 ^

Гв (t) = Г

(28)

1 - ехр

V

Предположим, что начальное распределение завихренности вихревой трубки соответствует распределению (27) в момент времени £0. Смоделируем задачу в описанной постановке с помощью ВО с учетом вязкости по методу РББ и сравним результаты с аналитическим расчетом по формулам (27) и (28).

Для моделирования выберем следующие численные параметры: радиус трубки Я = 0,8, разбиваемой на ВО, 2к = 0,1 (всего 1521 ВЭ в одном слое, см. рис. 3); расстояние между вортонами А «0,0021; параметр в модели РББ в = 0,003; вязкость V = 0,01; £0 = 1; циркуляция трубки Г = 1; характерный объем а = 1,32 -10~6 рассчитывался как произведение длины ВО 2к и части площади круга радиусом Я, приходящейся на один ВО. Интеграл в Gqk (14) в этой задаче вычислен аналитически с помощью пакета компьютерной алгебры. В силу громоздкости выражения не приведены. Уравнения (23) интегрировались явным методом Рунге — Кутты четвертого порядка с шагом А£ = 0,01 до достижения конечного времени 1.

Эволюция завихренности вортонов по радиусу трубки в интервале времени к е [0;1] с шагом Д^ = 0,2 приведена на рис. 5. Здесь под завихренностью понимается формальная величина, полученная делением интенсивности вортона на эквивалентную площадь, занимаемую им в круге. Точки на графике (20 точек)

X 1

,2

_ 3 С

А В

0,2

0,4

0,6

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

Изменение со временем приведенной завихренности трех характерных ВО А, В и С, находящихся на расстояниях гА =

Рис. 5. Эволюция профиля завихрен- = 0 (центральный), гв = °Д64, Гс = °,396 ности по радиусу круга при f = 0 (1), (сплошные линии), в сравнении с аналити-0,05 (2) и 1 (3) ческим расчетом по (27) (штриховые ли-

нии) приведено на рис. 6, а. Погрешность нарастает со временем и на момент времени I = 1 составляет 9,6 % для вортона А, 5,2 и 7, 7 % для вортонов В и С. Несмотря на сходство результатов, диффузия, описываемая численной моделью, происходит менее интенсивно, чем необходимо. Расчеты показывают, что погрешность зависит от выбранной длины ВО 2к. Степень влияния длины ВО на точность моделирования диффузии завихренности является объектом дальнейших исследований.

— ____К_ R — U,о = 0,4

' -—\ ,

R = 0,2

| 1 = 0,1

0 0,2 0,4 0,6 0,8 г 0 0,2 0,4 0,6 0,8 I а б

Рис. 6. Сравнение изменений завихренности трех характерных ВО гА = 0, гв = 0,164, гс = 0,396 (а) и изменения циркуляций для четырех кругов (б): сплошные линии — численный расчет; штриховые — аналитический расчет

Изменение циркуляций на четырех кругах радиусами 0,1, 0,2, 0,4, 0,8 (сплошные линии) в сравнении с расчетом по соотношению для циркуляций (28) (штриховые линии) показано на рис. 6, б. Наибольшая погрешность расчета циркуляции возникает для круга радиусом 0,2 и составляет 8,7 % в начальный момент времени. Эту неточность можно объяснить тем, что ВО заполняют круг дискретным образом (в отличие от реальной вихревой трубки с непрерывным распределением завихренности), что влияет на расчет циркуляции по контуру, проходящему вблизи ВО. Некоторым малым смещением контура круга внутрь или наружу начальную погрешность можно полностью устранить.

Анализ результатов, приведенных на рис. 6, показывает сходство исследуемой численной и аналитической моделей для задачи о диффузии вихревой трубки.

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

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

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

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

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

ЛИТЕРАТУРА

1. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 152 с.

2. Winckelmans G.S., Leonard A. Contributions to vortex particle methods for the computation of three-dimensional incompressible unsteady flows // Journal of Computational Physics. 1993. Vol. 109. Iss. 2. P. 247-273. DOI: 10.1006/jcph.1993.1216

3. Cottet G.H., Koumoutsakos P.D. Vortex methods: theory and practice. Cambridge University Press, 2000. 328 p.

4. Kotsur O., Scheglov G., Leyland P. Verification of modelling of fluid-structure interaction (FSI) problems based on experimental research of bluff body oscillations in fluids // 29th Congress of the International Council of the Aeronautical Sciences. 2014. P. 987-996.

5. Raviart P.-A. An analysis of particle methods // Numerical methods in fluid dynamics. Springer, 1985. P. 243-324.

6. Chorin A.J. Numerical study of slightly viscous flow // Journal of Fluid Mechanics. 1973. Vol. 57. Iss. 4. P. 785-796. DOI: 10.1017/S0022112073002016

7. Сизых Г.Б. Эволюция завихренности в закрученных осесимметричных течениях вязкой несжимаемой жидкости // Ученые записки ЦАГИ. 2015. № 3. С. 14-20.

8. Lacombe G., Mas-Gallic S. Presentation and analysis of a diffusion-velocity method // ESAIM: Proc. 1999. Vol. 7. P. 225-233. DOI: 10.1051/proc:1999021

9. Дынникова Г.Я. Движение вихрей в двухмерных течениях вязкой жидкости // Известия РАН. Механика жидкости и газа. 2003. № 5. С. 11-19.

10. Degond P., Mas-Gallic S. The weighted particle method for convection-diffusion equations. Part 1: the case of an isotropic viscosity // Mathematics of Computation. 1989. Vol. 53. No. 188. P. 485-507. DOI: 10.2307/2008716

11. Mycek P., Pinon G., Germain G., Rivoalen E. Formulation and analysis of a diffusionvelocity particle model for transport-dispersion equations // Computational and Applied Mathematics. 2016. Vol. 35. Iss. 2. P. 447-473. DOI: 10.1007/s40314-014-0200-5

12. Богомолов Д.В., Марчевский И.К., Сетуха А.В., Щеглов Г.А. Численное моделирование движения пары вихревых колец в идеальной жидкости методами дискретных вихревых элементов // Инженерная физика. 2008. № 4. С. 8-14.

13. Eldredge J.D., Leonard A., Colonius T. A general deterministic treatment of derivatives in particle methods // Journal of Computational Physics. 2002. Vol. 180. Iss. 2. P. 686-709. DOI: 10.1006/jcph.2002.7112

14. Брутян М.А. Точечная вихревая особенность в N-мерном пространстве // Фундаментальные и прикладные исследования: проблемы и результаты. 2014. № 10. С. 192197.

15. Beale J.T., Majda A. Vortex methods. I. Convergence in three dimensions // Mathematics of Computation. 1982. Vol. 39. No. 159. P. 1-27. DOI: 10.1090/S0025-5718-1982-0658212-5

16. Новиков Е.А. Обобщенная динамика трехмерных вихревых особенностей // ЖЭТФ. 1983. Т. 57. С. 556-562.

17. Басин М.А., Корнев Н.В. Аппроксимация поля завихренности в безграничном объеме // Журнал технической физики. 1994. Т. 64. № 11. С. 179-185.

18. Marchevsky I.K., Shcheglov G.A. 3D vortex structures dynamics simulation using vortex fragmentons // Proc. of ECCOMAS. 2012. P. 5716-5735.

19. Alkemade A.J.Q. On vortex atoms and vortons. PhD Thesis. Delft, Netherlands, 1994. 209 p.

20. Сэффмэн Ф.Дж. Динамика вихрей. М.: Научный мир, 2000. 375 с.

21. Лойцянский Л.Г. Механика жидкости и газа. М.: Дрофа, 2003. 840 с.

22. Самарский А.А. Введение в теорию разностных схем. М.: Наука, 1989. 616 с.

Коцур Олег Сергеевич — ассистент кафедры «Аэрокосмические системы» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, стр. 1).

Щеглов Георгий Александрович — д-р техн. наук, доцент, профессор кафедры «Аэрокосмические системы» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, стр. 1).

Просьба ссылаться на эту статью следующим образом:

Коцур О.С., Щеглов Г.А. Реализация метода обмена интенсивностями вортонов-отрезков для учета вязкости в методе вихревых элементов // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2018. № 3. С. 48-67. БО!: 10.18698/1812-3368-2018-3-48-67

IMPLEMENTATION OF THE PARTICLE STRENGTH EXCHANGE METHOD FOR FRAGMENTONS TO ACCOUNT FOR VISCOSITY IN VORTEX ELEMENT METHOD

O.S. Kotsur G.A. Shcheglov

oskotsur@gmail.com shcheglov_ga@bmstu.ru

Bauman Moscow State Technical University, Moscow, Russian Federation

Abstract

This paper focuses on the modification of the vortex element method for viscous incompressible flow simulations. Viscosity is accounted for using the Particle Strength Exchange (PSE) method by means of approximation of the viscous term in the vorticity evolution equation with an integral operator. Resolution equations for the parameters of the vortex elements (fragmentons) have been deduced: their marker position, direction and intensity. All the derivations are made by the assumption of noneffect of the vortex element additional vorticity. It has been shown that in the scope of this assumption, viscosity affects only the change of intensity of vortex elements. In such a case, in PSE method the elements follow trajectories of fluid particles, unlike in the Diffusive Velocity Method (DVM). Special focus is given to the discretization of initial vorticity onto the distributed system of vortex elements. Finally, the results of the diffusion problem of infinite vortex tube are described. Comparison of the numerical and analytical results confirms the correctness of the viscosity model for the considered problem

Keywords

Vortex element method, particle strength exchange method, fragmentons, viscosity, vorticity, vorticity evolution equation

Received 14.04.2017 © BMSTU, 2018

REFERENCES

[1] Patankar S. Numerical heat transfer and fluid flow. Hemisphere Publishing Corporation, 1980. 197 p.

[2] Winckelmans G.S., Leonard A. Contributions to vortex particle methods for the computation of three-dimensional incompressible unsteady flows. Journal of Computational Physics, 1993, vol. 109, iss. 2, pp. 247-273. DOI: 10.1006/jcph.1993.1216

[3] Cottet G.H., Koumoutsakos P.D. Vortex methods: theory and practice. Cambridge University Press, 2000. 328 p.

[4] Kotsur O., Scheglov G., Leyland P. Verification of modelling of fluid-structure interaction (FSI) problems based on experimental research of bluff body oscillations in fluids. 29th Congress of the International Council of the Aeronautical Sciences, 2014, pp. 987-996.

[5] Raviart P.-A. An analysis of particle methods. Numerical methods in fluid dynamics. Springer, 1985. Pp. 243-324.

[6] Chorin A.J. Numerical study of slightly viscous flow. Journal of Fluid Mechanics, 1973, vol. 57, iss. 4, pp. 785-796. DOI: 10.1017/S0022112073002016

[7] Sizykh G.B. Vorticity evolution in swirling axisimmetrical flows of viscous incompressible fluid. Uchenye zapiski TsAGI, 2015, no. 3, pp. 14-20 (in Russ.).

[8] Lacombe G., Mas-Gallic S. Presentation and analysis of a diffusion-velocity method. ESAIM: Proc., 1999, vol. 7, pp. 225-233. DOI: 10.1051/proc:1999021

[9] Dynnikova G.Ya. Vortex motion in two-dimensional viscous fluid flows. Fluid Dynamics, 2003, vol. 38, iss. 5, pp. 670-678. DOI: 10.1023/B:FLUI.0000007829.78673.01

[10] Degond P., Mas-Gallic S. The weighted particle method for convection-diffusion equations. Part 1: the case of an isotropic viscosity. Mathematics of Computation, 1989, vol. 53, no. 188, pp. 485-507. DOI: 10.2307/2008716

[11] Mycek P., Pinon G., Germain G., Rivoalen E. Formulation and analysis of a diffusionvelocity particle model for transport-dispersion equations. Computational and Applied Mathematics, 2016, vol. 35, iss. 2, pp. 447-473. DOI: 10.1007/s40314-014-0200-5

[12] Bogomolov D.V., Marchevskiy I.K., Setukha A.V., Shcheglov G.A. Numerical modelling of movement of pair vortical rings in the ideal liquid methods of discrete vortical elements. Inzhenernayafizika [Engineering Physics], 2008, no. 4, pp. 8-14 (in Russ.).

[13] Eldredge J.D., Leonard A., Colonius T. A general deterministic treatment of derivatives in particle methods. Journal of Computational Physics, 2002, vol. 180, iss. 2, pp. 686-709. DOI: 10.1006/jcph.2002.7112

[14] Brutyan M.A. Point vortex singularity in N-dimensional space. Fundamental'nye i pri-kladnye issledovaniya: problemy i rezul'taty, 2014, no. 10, pp. 192-197 (in Russ.).

[15] Beale J.T., Majda A. Vortex methods. I. Convergence in three dimensions. Mathematics of Computation, 1982, vol. 39, no. 159, pp. 1-27. DOI: 10.1090/S0025-5718-1982-0658212-5

[16] Novikov E.A. Generalized dynamics of three-dimensional vertical singularities (vortons). JETP, 1983, vol. 57, no. 2, pp. 566-569.

[17] Basin M.A., Kornev N.V. Approximation of vorticity field in unbounded volume. Zhurnal tekhnicheskoy fiziki, 1994, vol. 64, no. 11, pp. 179-185 (in Russ.).

[18] Marchevsky I.K., Shcheglov G.A. 3D vortex structures dynamics simulation using vortex fragmentons. Proc. of ECCOMAS, 2012, pp. 5716-5735.

[19] Alkemade A.J.Q. On vortex atoms and vortons. PhD Thesis. Delft, Netherlands, 1994. 209 p.

[20] Saffman P.G. Vortex dynamics. Cambridge University Press, 1992, 311 p.

[21] Loytsyanskiy L.G. Mekhanika zhidkosti i gaza [Fluid mechanics]. Moscow, Drofa Publ., 2003. 840 p.

[22] Samarskiy A.A. Vvedenie v teoriyu raznostnykh skhem [Introduction into the theory of difference schemes]. Moscow, Nauka Publ., 1989. 616 p.

Kotsur O.S. — Assistant, Department of Aerospace Systems, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation).

Shcheglov G.A. — Dr. Sc. (Eng.), Assoc. Professor, Professor, Department of Aerospace Systems, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation).

Please cite this article in English as:

Kotsur O.S., Shcheglov G.A. Implementation of the Particle Strength Exchange Method for Fragmentons to Account for Viscosity in Vortex Element Method. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Estestv. Nauki [Herald of the Bauman Moscow State Tech. Univ., Nat. Sci.], 2018, no. 3, pp. 48-67 (in Russ.). DOI: 10.18698/1812-3368-2018-3-48-67

В Издательстве МГТУ им. Н.Э. Баумана вышло в свет учебное пособие автора

Н.И. Сидняева

«Статистический анализ и теория планирования эксперимента»

Изложены краткие теоретические сведения по курсу «Теория планирования эксперимента». Представлено введение в статистический анализ и теорию планирования эксперимента. Основные понятия проиллюстрированы примерами практического содержания, рассмотренными с позиций регрессионного анализа. Издание носит справочный характер и поможет студентам старших курсов овладеть методами теории планирования эксперимента, которые широко используются при решении прикладных задач. Для студентов четвертого-шестого курсов инженерных специальностей технических университетов.

По вопросам приобретения обращайтесь:

105005, Москва, 2-я Бауманская ул., д. 5, стр. 1 +7 (499) 263-60-45 press@bmstu.ru www.baiunanpress.ru