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

Вычислительные технологии

Том 12, № 5, 2007

НОВЫЙ ПОДХОД К МОДЕЛИРОВАНИЮ ВЯЗКОСТИ В МЕТОДЕ ДИСКРЕТНЫХ ВИХРЕЙ*

, О. А. Шмагунов Институт теоретической и прикладной механики СО РАН

Новосибирск, Россия e-mail: shmag@itam.nsc.ru

A modification of the discrete vortices method is suggested that accounts for viscosity effects by means of corrections of the equations of motion. The correction is determined by the requirement that variations of the motion invariants have to coincide with the variations generated by a given viscosity. A good agreement with experiment is obtained.

Введение

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

1. Основная идея метода

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-01-00080-а).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

Б.Ю. Скобелев

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

1.1. Гипотеза вязкости

В двумерных течениях покоящейся на бесконечности невязкой жидкости энергия Ш инвариантна [4]:

где р — плотность, и = и(х, у), и' = и'(х', у') — распределения завихренности, r = (х, у), интегралы берутся по всей области течения.

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

Распределение завихренности внутри кругового вихря определяется бесконечно дифференцируемой функцией и(р), которая постоянна внутри круга меньшего радиуса (р — расстояние от центра вихря, и(р) = const при р < r0 — е). Можно показать [5], что

(1)

со Л

50=1

Го Р

Рис. І. Дискретизация поля завихренности

в случае вязкого течения Ш не сохраняется и изменяется во времени в соответствии

с законом

(2)

где ^р — коэффициент вязкости, Жр определяется формулой (1), интегралы берутся по области 5р.

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

N

ті и2 <і5 = г!,

ПГР

^0

к=1

где Г& — циркуляция к-го вихря, N — число вихрей в области 50.

Энергия завихренности Ш0 может быть разложена на три компоненты:

Шо = /1 + /2 + /3,

(3)

(4)

где

її = — — 1пгр /[ии'^Ж', /2 = —— [[ии' 1п ^^Ж',

4п ]] 4^^ гр

/3 = — 4“ УJ ии' 1п |г — Г| ^5^5, П = {(г, Г) : г Є £р, — Є 5р, |г — Г| ^ 2гр} .

й0\п

При сравнительно малом гр величина /1+/2 является суммой собственных кинетических энергий N круговых вихрей, где /1 — сингулярная часть суммы, /2 — регулярная часть. Можно показать [6], что /1 =сопэ1 и

32п

(5)

І=1

Интеграл /3 описывает энергию взаимодействия вихрей и при е ^ 0 совпадает с энергией взаимодействия точечных вихрей:

N

/з = -4^5]Г*Г^'1п И - о1>

где гг = (хг,уг) — координаты центра г-го вихря.

С учетом формул (3)-(6) уравнение (2) принимает следующую форму:

(6)

N

= ^р^р 1 г2 е = р

^ пгр 5р г

р г=1

■ N 1 N

32лЁ г2 — г<г;1п |гі — г I

І=1

і=3

Пусть гр ^ 0, и одновременно потребуем, чтобы ^р ^ 0 в соответствии с законом

^р = у-

пгр

5р

где ^ — заданное значение коэффициента вязкости. Тогда, учитывая равенство S0 = 1, получаем

dE N

§ = £ г?- (8)

i=1

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

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

1.2. Неинерциальное преобразование системы координат

Уравнение движения N точечных вихрей в гамильтоновой форме записывается как ^ dxk dH ^ dyk dH

Tk-k = ^-, Tfc-f = - — , k =1, 2,..., N. (9)

dt dyk dt dxk

Гамильтониан H совпадает с (6), если р =1. Известно [4], что координаты X, Y центра завихренности, гамильтониан H и дисперсия завихренности D? являются инвариантами уравнений (9). Вместо D? рассмотрим связанный с ней инвариант:

N

L? = £ Г,(х? + у?) (10)

i=1

(с точностью до численного множителя L? совпадает с моментом импульса). Интегральные аналоги величин X, Y, L? являются инвариантами невязкого течения с непрерывным распределением завихренности. Можно показать [4, 5], что при заданной вязкости ^0 величины X, Y по-прежнему сохраняются, а интегральный аналог L? изменяется по закону

— [ ш(х? + y?)dS = [ wdS. (11)

dt J р J

Как и в предыдущем разделе, аппроксимируем завихренность ш дискретными круговыми вихрями и используем гипотезу вязкости: ^0 = ^пг?. Далее, можно показать [6], что величина L? (10) должна оставаться постоянной в процессе моделирования вязкого течения точечными вихрями:

f =« „2,

Поскольку циркуляция точечных вихрей и гамильтониан H сохраняются в процессе движения, то E = const на точных решениях уравнений (9) и, следовательно, условие (8) не выполняется. С другой стороны, дискретизация уравнений (9) ведет к появлению погрешности, а E и L? в численных расчетах изменяются случайным образом. Условия (8) и (12) также не выполняются. Чтобы выполнить условия (8) и (12), примем следующую гипотезу.

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

Рассмотрим конечно-разностные формы уравнений (8) и (12):

N

- \ Л т^2

ДЕ = -^Д*2_, Г2, ДЬ2 = 0, (13)

г=1

где V — кинематическая вязкость, а Е обозначает соотношение (7) при р =1. Обозначим через £г)П+1, Уг,га+1 координаты г-го вихря, полученные на (п + 1)-м шаге интегрирования уравнений (9), а через Г^ — соответствующее значение циркуляции (оно не меняется при переходе с п-го на (п+1)-й слой по времени). Преобразование в физическую систему координат будем искать в следующем виде:

хг,га+1 (1 + ДЬга+1) хг,го+1) (14)

уг,га+1 (1 + ДЬга+1) уг,п+1 ?

Гг,га+1 = (1 + ДГп+1) Гг,га, г = 1, 2, . . . , (15)

где ДЬп+1, ДГга+1 — некоторые параметры. Из (9) видно, что уравнения движения в

новых переменных сохраняют прежний вид, если шаг по времени Д£, используемый в интегрировании, преобразуется по формуле

Л+ =(1+ДЬп+1)2 (16)

м'“+1 = (1 + ДГ„+1) (16)

Значения параметров ДЬп+1, ДГп+1 определяются из условий (13) на приращения ДЬП+1 = ЬП+1 - и ДЕп+1 = Еп+1 — Еп, в которые подставляются выражения (14) и (15):

(1 + ДГ)(1 + ДЬ)2/П+1 — = 0, (17)

1 N N

(1 + ДГ)2е„+1 — — (1 + ДГ)2 1п(1 + ДЬ) Гг,гаГ^',га — Еп = — VД^У^ (18)

1= г=1

где 1П+1 и еп+1 — дисперсия и энергия, рассчитанные для £г)П+1, Уг,п+1. Система уравнений (17)—(18) решается численно, например, из первого уравнения выражается ДГ и подставляется во второе, которое решается методом итераций. После преобразований (14) и (15) мы получаем величины координат и циркуляций идеальных точечных вихрей, моделирующих течение с заданной кинематической вязкостью V.

Преобразования (14) и (15) удобно выполнять в два этапа. На первом этапе делается только преобразование координат. Параметр преобразования получается из (17) и имеет следующий вид:

Г 2 \

ДГо.п+1 = ( - 1.

^2 \ 1/2 ^°,п+і = ^ у2“ I

ч ^га+1/

После этого значение энергии е0;П+1 и величина /2„+1 рассчитываются в преобразованных координатах.

На втором этапе выполняются полные преобразования (14) и (15). По определению, 12п+і = ГП, следовательно, получаем из (17):

ДГі = (-±-Г -1

1 V і + ДГі/

Подставляя это выражение для ДГ1 в (18), получаем нелинейное уравнение на ДГ1, которое решается методом итераций.

2. Моделирование обтекания плоской пластины

Рассмотрим традиционную для методов дискретных вихрей задачу обтекания плоской пластины потоком жидкости с постоянной скоростью и0 (рис. 2).

В методе дискретных вихрей пластина замещается системой неподвижных точечных вихрей, расположенных на отрезке прямой. Точки отрыва свободных вихрей располагаются на фиксированном расстоянии от краев пластины. Два свободных вихря порождаются перед каждым шагом интегрирования уравнений движения. Циркуляции неподвижных и свободных вихрей определяются из двух условий: условия непротека-ния в контрольных точках, расположенных между неподвижными вихрями и на краях пластины, и равенства нулю полной циркуляции. Таким образом, есть два типа вихрей в задаче перед переходом к (п + 1)-му шагу по времени: N неподвижных вихрей и 2п свободных вихрей. Свободные вихри движутся в соответствии с уравнениями (9), дополненными вкладом от неподвижных вихрей и внешнего течения.

Преобразования координат и циркуляций, описанные в предыдущем разделе, выполняются на каждом шаге по времени, начиная с некоторого временного слоя п0, перед порождением двух новых вихрей. При этом принимается во внимание, что N свободных вихрей, находящихся на (п + 1)-м временном слое ^ = 2п + 2), движутся во внешнем поле скорости, порождаемом набегающим потоком и неподвижными вихрями. Влияние внешней скорости и0 исключается переходом в систему координат, движущуюся с внешним потоком.

Скорость и = (их,иу), порожденная неподвижными вихрями, приводит к дополнительным изменениям энергии Е и величины Ь2 [5]. На интервале времени Д£ эти изменения имеют вид

N N

ДЕ = Д£ ^Г^ПхЛиуг - Пуг ^г), ДЬ2 = Д? ^ ^(2^x1 + Уг^), (19)

г=1 г=1

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

В проведенных численных экспериментах пластина замещалась N0 = 20 неподвижными вихрями. Точки отрыва свободных вихрей располагались на расстояниях 1/(2^) от краев пластины, ширина которой была взята за единицу. Неодновременность начала движения двух новых свободных вихрей оказалась существенной для формирования вихревой дорожки Кармана. Экспериментально было установлено, что тот из вихрей,

у 4 1.1

и0 ;5!

» >0'

!х!

■о;

■ х>

!о!

гг

Рис. 2. Расчетная область

циркуляция которого была меньше, должен был отрываться первым. Таким образом, уравнения движения (9) интегрировались явно-неявным методом Эйлера следующего вида:

Гаусса—Зайделя для решения системы линейных алгебраических уравнений. Вихрю с меньшей циркуляцией перед интегрированием присваивался четный номер. Численное интегрирование производилось с шагом по времени Ді = 0.1. Дальнейшее уменьшение шага существенно не отразилось на картине течения, что свидетельствует в пользу того, что метод успешно компенсирует погрешность дискретизации поля завихренности. Эффективность компенсации погрешности дискретизации уравнений движения оценивалась при помощи тестовых расчетов с использованием метода Рунге—Кутты второго порядка. Как и в случае с варьированием шага по времени, повышение порядка метода интегрирования не отразилось существенно на результате. Для регуляризации численного моделирования точечные вихри замещались круговыми вихрями с радиусом дискретизации г0 = 1/40. Внутри круга индуцированная скорость линейно уменьшалась до нуля в центре. Дальнейшее уменьшение радиуса дискретизации, в частности до г0 = 1/80, не повлияло на результат.

На расстоянии четырех длин пластины вниз по потоку были расположены два счетчика завихренности Сі, С2 (прямоугольники размером 1 х 0.3, суммирующие циркуляции всех попадающих в них вихрей), а также измеритель V поперечной скорости V в точке х = 5, у =1 (колебания скорости усреднялись по пяти точкам).

Результаты показали, что предлагаемый метод учета вязкости хорошо моделирует вихревую дорожку Кармана и зависимость числа Струхаля БЬ от числа Рейнольдса И,е. На рис. 3 приведены примеры вихревой дорожки Кармана в момент і = 70 для Ие = 100 и 1000.

Графики циркуляций, насчитанных счетчиками Сі, С2, и график поперечной скорости при И,е = 50 даны на рис. 4.

Сравнение значений числа Струхаля БЬ с экспериментальными данными [7] в диапазоне 50 ^ И,е ^ 1000 даны на рис. 5. Сплошные линии показывают экспериментальную зависимость БЬ от И,е, кружки обозначают результаты численного моделирования.

(20)

где /, — правые части уравнений (9). Метод (20) подобен итерационному методу

а

б

Рис. 3. Вихревая дорожка Кармана: і = 70, Ие = 100 (а) и 1000 (б)

Рис. 4. Графики циркуляций и поперечной Рис. 5. Зависимость числа Струхаля ЯЬ от скорости числа Рейнольдса Ие

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

Тестовые расчеты показали, что начальная стадия движения пластины 0 < і < 10 может быть достаточно точно смоделирована с временным шагом Ді = 0.04. Для исследования аэродинамических характеристик при і > 10 уравнения движения можно моделировать с шагом по времени Ді = 0, 1.

Расчеты, проведенные при числах Рейнольдса И,е = 50 и И,е = 2000, показали, что при больших углах атаки а = 60° и а = 90° (поперечное обтекание) вид функций сп(і) и ха(і) при малых временах 0 < і < 10 практически не зависит от числа Рейнольдса. При меньших углах атаки а = 30° вязкость играет более существенную роль (рис. 6).

Рис. 6. Зависимость функции сп(і) и Ха(і) от числа Рейнольдса: а = 30°, Ие = 50 (а) и 2000 (б)

Рис. 7. Зависимость переходного режима сп от различных чисел Рейнольдса: а = 90°, Ие = 50 (а) и 2000 (б)

При И,е = 50 и і > 5 изменение сп монотонно, но при И,е = 2000 сп осциллирует. Усредненные на интервале 7.5 < і < 10 значения сп также существенно различаются: сп = 1.21 при И,е = 50 и сп = 1.03 при И,е = 2000.

При больших углах атаки, в частности при а = 90°, влияние вязкости становится существенным на больших временах.

Рисунок 7 показывает, что переходный режим для сп (17 < і < 52) качественно различается при различных числах Рейнольдса. Периоды установившихся при і > 52 осцилляционных режимов также зависят от числа Рейнольдса, хотя усредненные на интервале 52 < і < 70 значения сп мало различаются: сп = 3.27 при И,е = 50 и сп = 3.2 при И,е = 2000.

Заключение

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

Аналогичный подход может быть развит для трехмерных течений. Трехмерные идеальные течения имеют следующие инварианты: полная завихренность, импульс, момент импульса, кинетическая энергия и спиральность. Для вязких течений первые три величины по-прежнему сохраняются, а энергия и спиральность изменяются по известным законам [4, 8]. Непрерывное поле завихренности может быть аппроксимировано набором коротких сегментов вихревых нитей. Тогда появляется возможность разработать численный алгоритм для моделирования вязкости, подобный алгоритму для двумерных течений.

Список литературы

[1] Belotserkoysky S.M., Lifanoy I.K. Method of Discrete Vortices. USA: CRC Press, 1997

[2] Chorin A.J. Numerical study of slightly viscous flow // J. Fluid Mech. 1973. Vol. 57. P. 785-796.

[3] Winokelmans G.S., Leonard A. Contributions to Vortex Particle Methods for the Computation of Three-Dimensional Incompressible Unsteady Flows // J. Comp. Phys. 1993. Vol. 109. P. 247-273.

[4] Бэтчелор Дж. Введение в динамику жидкости. М.: Мир, 1973.

[5] Шавалиев М.Ш. Законы изменения моментов распределения завихренности под влиянием вязкости, внешнего поля скорости и наличия твердых границ // Механика неоднородных и турбулентных потоков: Сб. науч. тр. М.: Наука, 1989. С. 63-69.

[6] Belotserkoysky S.M., Soobeley B.Yu., Shmagunoy O.A. Viscosity simulation in the method of discrete vortices // Computational Fluid Dynamics’96. Proceedings of the Third ECCOMAS Computational Fluid Dynamics Conference, 9-13 September 1996. P., France. P. 791-796.

[7] Roshko A. On the Development of Turbulent Wakes from Vortex Streets // NACA TN. 1953. N 2913.

[8] Moffat H.K. The degree of knotedness of tangled vortex lines // J. Fluid Mech. 1969. Vol. 35. P. 117-129.

Поступила в редакцию 24 июля 2006 г., в переработанном виде — 6 июня 2007 г.