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

УДК 532.526

МОДЕЛИРОВАНИЕ ДВУМЕРНОГО ЛАМИНАРНОГО ПОГРАНИЧНОГО СЛОЯ С ПОМОЩЬЮ МЕТОДА «ВИХРЬ В ЯЧЕЙКЕ»

© 2006 В. В. Никонов, В. Г. Шахов Самарский государственный аэрокосмический университет

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

1. Математическая формулировка метода ВЯ

Метод «вихрь в ячейке» (ВЯ) часто используется для расчета отрывных течений в рамках модели вязкой несжимаемой жидко -сти. Данный метод относится к группе методов дробных шагов для решения уравнения Навье-Стокса в безразмерных переменных завихренность-скорость, которое для двумер -ного случая имеет вид

где ю- завихренность О V О и

(1)

ю =

О х О у ’

(2)

і + Ь/2

Г',.= I

X; - Ь/2

у^ Ь/2

І ю(х,у)ёхёу.

уД/2

(3)

Поле скорости в вихревых методах представляется в виде суммы постоянной и ¥ и соленоидальной и у составляющих

и = и ¥ + и у ,

(4)

где и ¥ - скорость набегающего потока, а и у определяется как

и у =УХ (уег) .

(5)

Здесь у - функция тока течения, ег -единичный орт оси г.

Соленоидальная составляющая поля скорости находится из решения уравнения Пуассона для функции тока

Ау = -ю

(6)

1 - время, и - вектор скорости, V - безразмерная кинематическая вязкость п=1/Яе, Яе -

число Рейнольдса, V - оператор Гамильтона, А - оператор Лапласа.

В численной схеме метода ВЯ вместо завихренности используется вихревая интенсивность, по размерности совпадающая с циркуляцией

методом быстрого преобразования Фурье [1, 2] с последующим взятием операции ротора (5). Число необходимых операций составляет при этом порядок N 1о§2(К), N = = пхпу - число ячеек расчетной сетки. Граничные условия (ГУ) для уравнения (6) определяются для функции тока течения, вызванного потенциальными вихрями, находящимися в ячейках

1

у(х) = -X гу1пх-^

2р j=l ;=1

(7)

П, П

для рассматриваемой точки границы х. Здесь X ц - радиус-вектор вихря в ячейке Г^.

В некоторых случаях в методе ВЯ [3, 4] соленоидальная составляющая находится как сумма полей скорости, индуцированных потенциальными вихрями (закон Био-Савара)

uy (x) = _! £ £ Flj£-X (Х 2Xj)

^ 2p£ £ (x-Xij)2

(8)

Число необходимых операций при применении схемы (8) составляет порядок N2 где N = пп.

х У

Для сокращения времени вычисления ГУ (7) или определения скорости по (8) применяется кластеризация и метод мультипольного разложения [2].

При решении уравнения (1) методом ВЯ процессы конвекции, диффузии завихренности в свободном потоке и с поверхности тела на каждом шаге по времени рассматриваются отдельно [1, 5].

Алгоритм метода ВЯ состоит в решении на каждом шаге по времени следующих задач:

1. Определение вихревой интенсивности на поверхности тела.

2. Диффузия завихренности в свободном потоке.

3. Диффузия вихревой интенсивности с поверхности обтекаемого тела.

4. Определение скорости и конвекция жидких частиц.

5. Сохранение результатов расчета и переход к следующему шагу.

1.1. Определение вихревой интенсивности на поверхности тела

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

где и! - величина скорости в точке 1 на поверхности тела, рассчитываемая по (4); Г1 -циркуляция присоединенного дискретного вихря, определяемая из решения системы линейных алгебраических уравнений, которая получается при удовлетворении условиям не-протекания; 11 - размер вихревой панели присоединенного вихря Г1 (расстояние между соседними контрольными точками).

1.2. Диффузия в свободном потоке Для моделирования процесса диффузии в свободном потоке используется схема «донор-акцептор» (Д-А) [7]:

Я еМр

г„(1+а 1)=г„(о + х (г,(1)о;,(х)о;,(у) -

-rp(t)G;(x)G;(y)),

(10)

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

4+1;

erf

/zp + h/2 - Zq Л

V4vAt

erf

/zp - h/2 - Zq4

V4VAt

(11)

где z-xp или yp - координаты центра ячеики p-го вихря, erf - интеграл вероятности (функция ошибок) [8]. В работе [9] показано, что для достижения заданной точности шаг по времени для метода Д-А определяется следующим соотношением:

xt = kdh / v ,

(12)

gi± = Ш • S ±Г; /(2li),

(9)

где кё - константа, зависящая только от радиуса пё «диффузионной молекулы» и для пё = 1 удовлетворяющая условию: 0,2< кё < 0,21.

q

Здесь принимается ка = 0,21, так как в данном случае [5] ошибки методов Д-А и моделирования процесса конвекции при использовании подхода Лагранжа с последующим перераспределением вихрей в ячейки расчетной сетки будут иметь разные знаки и компенсировать друг друга.

1.3. Диффузия с поверхности обтекаемого тела Диффузия завихренности с поверхности тела в связанной с вихревой панелью системе координат определяется [1, 10] следующим образом:

Г =-т± <л)ь

егі"

л/4^

-егГ

(13)

Гij = Г(х^ уі )Л(хі - хк )Л^ - уі; ,

(14)

Л1(ъ)=

(15)

где ъ* = |ъ|/Ь.

1.4. Конвекция вихревых частиц После расчета поля скоростей (4) новые координаты «вихрей в ячейках» получаются

(аналогично методу дискретных вихрей) численным интегрированием системы обыкновенных дифференциальных уравнений движения методом Эйлера

-1+а 1 -і -

х; = х; + и А1.

(16)

Новое местоположение вихрей не обязательно совпадет с координатами расчетной сетки, поэтому используется процедура (14) перераспределения их интенсивностей в ячейки сетки. При этом количество вихревых частиц не растет с течением времени, как в бессеточных вихревых методах. В качестве интерполяционных функций в данной работе применяется формула М4’ [11]

ЛзМ4'(ъ) =

^ 5 *2 3 *3

2 2 0 < ъ < 1

(1 -ъ*)(2 - ъ*)2 / 2, 1 < ъ* < 2 (17) 0, ъ* > 2

где у± ^) находится из выражения (9), И -продольный размер ячейки, связанной с панелью сетки. Вертикальный размер ячейки сетки, связанной с панелью, равен шагу глобальной сетки И.

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

где Л - интерполяционная функция; Г(хк,у1) -интенсивность перераспределяемого вихря, находящегося в произвольной точке с координатами (хк,у1); Г - циркуляция, получаемая вихрем в ячейке (і, _]) от перераспределяемого вихря. В качестве интерполяционной функции используется «облако в ячейке» [11, 12]:

для перераспределения вихрей, которые находятся на расстоянии двух размеров ячейки сетки от поверхности тела. Если перераспределяемый вихрь находится вблизи тела, то используется формула «облако в ячейке» (15). Формула (17) имеет меньшую численную диффузию, чем (15). Однако последняя использует минимальное количество ячеек для интерполяции и не вносит пульсаций в поле перераспределяемой величины.

2. Прямое численное моделирование ламинарного пограничного слоя на плоской пластине

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

и величины вертикальной компоненты ско -рости

v u :x

Vf=u:Vv

профили компонент скорости u и vf в разных сечениях пластины x = const будут совпадать, что является удобным для сравнения получаемых данных.

При рассмотрении продольного обтекания плоской пластины в выражениях (13) и (24) в данной работе для удобства принималось hs = h = l;.

Сравнение профилей скорости с решением Блазиуса производилось для четырех сечений на расстояниях: xx » 0,25, x2 » 0,5, x3 » 0,75 и x4 » 0,9 от переднего края пластины .

При численном моделировании расчетная область имела следующие размеры: длину Ь = 4, высоту Н = 2. Расстояние от левой границы до центра пластины составляло х = -1,0, от нижней - у = 1,0. Сетка имела

а 5 5 а 5

400x200 ячеек (И = 0,01). Рассматривалась пластина единичной длины, и на ней располагалось 100 вихревых особенностей. Радиус ядра присоединенных вихрей равнялся

о = 0,5 1. .

5 1

Результаты расчета, полученные для числа Рейнольдса Яе = 10 3 (таблица 1, рис. 1, 2), показывают, что наблюдается неплохое соответствие с решением Блазиуса в окрестности поверхности пластины для продольной компоненты скорости (рис. 1). В то же время в верхней части пограничного слоя для полученного профиля скорости наблюдается больший «разгон» потока. Такой же «разгон» наблюдается и для данных, приведенных в [14]. Профиль скорости, полученный в [15], имеет хорошее согласование с решением Блазиуса в верхней части погра-

Таблица 1

Максимальная погрешность численного решения метода ВЯ в каждом из рассматриваемых сечений в момент времени 1 = 16,8 (Яе = 10 3, И = 0,01)

Номер сечения Расстояние от начала пластины до рассматриваемого сечения 8u 8v

1 0,245 0,02 0,11

2 0,495 0,04 0,18

3 0,745 0,07 0,40

4 0,895 0,10 1,14

0 1 2 3 4 5 6

h

Рис. 1. Профиль продольной скорости на плоской пластине в сравнении с результатами других авторов, ВЯ, Д-А, Re = 103, Atc = Atd = 0,021, t = 16,8, h = 0,01, LxH = 4x2;

----- - профиль Блазиуса [13],---- - Wu [14] (xn = 0,5), ..-Ota [15] (x = 0,9),

настоящая работа: ^ - хп = 0,25, + - хп = 0,5, Q - хп = 0,75, X - хп = 0,9

Vf

0.9

0.7

0.5

0.3

0.1

-0.1

4- 4-

y' + '4 4 T

/Ї cr □ □ □ □ □ с □ □

уf Г n ] v X X X X XXX

x X X X X X X X x

0

3

4

5

6

h

Рис. 2. Профиль вертикальной скорости на плоской пластине в сравнении с результатами других авторов, ВЯ, Д-А, Re = 103, А! = А! = 0,021, t = 16,8, к = 0,01, ЬхИ = 4x2;

■ профиль Блазиуса [13],

настоящая работа: ^ - хп = 0,25, + - хп = 0,5, Q - хп = 0,75, X - х„ = 0,9

ничного слоя, но имеет сильное отклонение в районе поверхности пластины. Профиль вертикальной скорости (рис. 2) хорошо согласуется с результатами решения Блазиуса только для первого сечения, расположенного при хх » 0,25. Для последнего сечения х4 » 0,9 наблюдается область с отрицательной скоростью, что качественно отличается от решения Блазиуса. Отметим, что в работах [14, 15] профиль вертикальной компоненты скорости не приводится.

В таблице 1 приводятся относительные погрешности профилей продольной 8ц и вертикальной 8у скорости. Результаты, показанные на рис. 1, 2 и в таблице 1, получены с помощью схемы Д-А моделирования процесса диффузии (Д^ = Д^ = 0,021). Точно такие же результаты (для Яе = 10 3) были получены с помощью схемы Д-А при интегрировании с разными шагами по времени (ИРШ) для Д^ = 0,001 и Д^ = 0,021, и поэтому они здесь не приводятся.

Для числа Рейнольдса Яе = 100 (V = 10 -2) согласование с решением Блазиуса сильно ухудшается. Измельчение расчетной сетки к положительному результату не приводит. При Яе = 10 4 (V = 10 -4) результаты численного моделирования согласуются с решением Блазиуса еще хуже. По этой причине результаты для Яе = 100 и 10 4 здесь не показаны. Заме-

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

Полученные результаты позволяют сделать следующие выводы:

- при применении стандартной схемы вычисления скорости метода «вихрь в ячейке» удается получить, как и в [3, 14, 15], удовлетворительные результаты только для продольной компоненты скорости u для числа Рейнольдса Re = 10 3 и в его некоторой окрестности;

- распределение вертикальной компоненты скорости v качественно отличается от аналитического решения наличием течения к поверхности пластины в ее задней части.

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

1. Taranov A., Kornev N., Leder A. Development of the Computational Vortex Method for Calculation of Two-Dimensional Ship Sections with Flow Separation, Schiffbauforschung, 39, 2000, 2, pp. 95-105.

2. Kornev N., Leder A., Mazaev K. Comparison of two fast algorithms for the calculation of flow velocities induced by a threedimensional vortex field, Schiffbauforschung,

40, 2001, 1, pp. 47-55.

3. Nakamura H., Kamemoto K., Igarashi T Analysis of unsteady heat transfer in the wake behind a circular cylinder in a uniform flow by a

1

2

vortex and heat element method, Proceedings of the Second International Conference on Vortex Methods, Sept. 26-28, Turkey, 2001, pp. 235242.

4. Iida A., Kamemoto K., Ojima A. Prediction of aerodynamic sound spectra by using an advanced vortex method, Proceedings of the Second International Conference on Vortex Methods, Sept. 26-28, Turkey, 2001, pp. 235242.

5. Никонов В. В., Шахов В. Г. Модификация схемы «донор-акцептор» для расчета диффузии завихренности и ее применение в методе «вихрь в ячейке» // Вестник СГАУ, № 1 (3), Самара, 2003. - С. 38-46.

6. Белоцерковский С. М., Ништ М. И. Отрывное и безотрывное обтекание тонких крыльев идеальной жидкостью. - М.: Наука, 1978.

7. Basin M., Kornev N. Beruecksichtigung der Reibung in der Wirbelmethode, ZAMM, 78, 1998, 5, pp. 335-344.

8. Справочник по прикладной статистике в 2-х т.: Под ред. Э. Ллойда, У. Ледерма-на. - М.: Финансы и статистика, 1989.

9. Nikonov V., Kornev N., Leder A. The

Ratio between Spatial and Time Resolutions for the Diffusion Substep in 2D Computational Vortex Methods, Schiffbauforschung, 2002, vol.

41, N 3/4. pp. 5-12.

10. Koumoutsakos P., Leonard A. Highresolution simulations of the flow around an impulsively started cylinder using vortex methods, J. Fluid Mech., 296, 1995, pp. 1-38.

11. Cottet G.-H., Koumoutsakos P Vortex methods: theory and practice, Cambridge University Press, 2000.

12. Григорьев Ю. Н., Вшивков В. А. Численные методы «частицы-в-ячейках». -Новосибирск: Наука, Сибирская издательская фирма РАН, 2000.

13. Шлихтинг Г. Теория пограничного слоя: Пер. с нем. Г. А. Вольперта, под. ред. Лойцянского Л. Г. - М.: Наука, 1974.

14. Wu J.C. Numerical boundary conditions for viscous flow problems, AIAA Journal, 14, 1976, pp. 104-1049.

15. Ota S., Kamemoto K. Study on higher resolution of vorticity layer over a solid boundary for vortex methods, Proc. of The Second Intern. Conf. on Vortex Methods, Istanbul, Turkey, September 26-28, 2001, pp. 33-40.

SIMULATION OF A TWO-DIMENSIONAL LAMINAR BOUNDARY LAYER USING THE METHOD «WHIRL IN A CELL»

© 2006 V. V. Nikonov, V. G. Shakhov Samara State Aerospace University

The paper deals with the use of the «whirl in a cell» (WC) method for direct mathematical simulation of a laminar boundary layer. According to the method the processes of convection and those of diffusion both in the free flow and from the surface of the body are considered separately at each step in time. Because of the difference in velocities of diffusion and convection processes integration with different steps in time is used. The velocity profiles obtained are compared with the Blazius solution and the results of other authors. It is shown that using the «WC» method good agreement with Blazius’ results is obtained only in a certain narrow area of Reynolds’ numbers.