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

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

Том 14, № 4, 2009

Метод расщепления в задаче численного моделирования турбулентного следа за буксируемым телом в стратифицированной жидкости*

Н. П. Мошкин

Учреждение Российской академии наук Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: [email protected]

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

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

Введение

Турбулентные следы за телами вращения в стратифицированной жидкости рассматривались во многих публикациях, где содержится достаточно полный обзор ссылок по этой проблеме. Укажем лишь несколько последних, не претендуя на полноту [1-4]. Анализируя цитированную литературу, можно сделать вывод о недостаточной полноте численных моделей динамики турбулентного следа за буксируемым телом в линейно стратифицированной среде (прежде всего для расстояний, соответствующих временам после прохода тела t < 10T, T — период Вяйсяля—Брента). Недостаточно и сопоставление параметров турбулентных следов за буксируемым и самодвижущимся телами в линейно стратифицированной среде. В настоящей работе численный подход, основанный на конечно-разностных методах расщепления, используется для построения численных моделей дальнего турбулентного следа за буксируемым телом в линейно стратифицированной среде.

1. Постановка задачи

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

* Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (грант № 07-01-00363), гранта НШ 9886.2006.9 Президента РФ и СО РАН (интеграционный проект 26).

© ИВТ СО РАН, 2009.

(1)

dV

dV

+ V— + W— = -

dV 1 d (pi) d

d

dx

dy

dz

TT dW dW flW

U0— + V— + W-

P0 dy 1 d (pi)

- eyV2>- dZ (v'w')

dx

dy

dz

Po

dz

d , , a d . (Pi)

- TT(v w) - TT \w /- g—■

dy dz x ' p0

Uo+ V^M + W^ + Wds = - d (v'p') - d Ш

ox dy dz dz dy dz

dV dW dy dz

dUd dx

(2)

(3)

(4)

(5)

В уравнениях (1)-(5) величина Ud = U0 - U — дефект осредненной продольной компоненты скорости; U,V,W — компоненты скорости осредненного движения в направлении осей x,y,z соответственно; (pi) — отклонение давления от гидростатического, обусловленного стратификацией ps(z); U0 — скорость набегающего невозмущенного потока; g — ускорение силы тяжести; (pi) — осредненный дефект плотности: Pi = p — ps, ps = ps(z) — плотность невозмущенной жидкости, которая полагается линейной: ps(z) = p0(1 — az),a = const > 0; штрихом обозначены пульсационные компоненты; символ () — осреднение. Плотность жидкости считается линейной функцией температуры, стратификация предполагается слабой. Система координат связана с перемещающимся телом таким образом, что скорость его движения равна —U0 и ось z направлена вертикально вверх противоположно силе тяжести. В уравнениях (1)-(4) отброшены в предположении малости члены с молекулярной вязкостью и диффузией. Отброшены также производные по x в правых частях.

Для замыкания системы (1)-(5) используется иерархия моделей турбулентности второго порядка. В модели 1 неизвестные величины компонент тензора рейнольдсовых напряжений (u'j2), i = 1, 2, 3, (u'v') = (ulU), (u>w') = (u^n's) и турбулентных потоков (uip'), i = 1, 2, 3, определены алгебраическими соотношениями (см. [4, 5]):

(u'iu'j)

oSij +

1 - С2 Ci

Pij 2 P

- - 3j

£

£

+

1 - Cs Ci

Gij 2 с G

— - о 6ij~

£ 3 £

(6)

I , ' ', dUj . ' '. dUi

pij = -1 (uiuk) + (uj uk) dx

Gij = — ((uip' )gj + (u'j p' )gi), P0

k dxk

g = (0, 0, -g), 2P = Pii, 2G = Gii, Ui = U, U2 = V, Us = W,

e

-(u'p')

CiT £

, ' Л d(p) . . ' '. dU

(uw )+ (1 - C2T)(w p ) —

-(v'p')

(v'2) ed(p)_T^ d(p) 2 eu.JJXd(p)

CiT £ dy

K

ev

dy

(P' ) = - CT' ~£ (w'P'] dz

-(w'p') =

Cit £

{w'2) dh

+ (1 - C2T) g (p'2)

p0

e(w'2)

d (P)

c i i О1 - C2T g e2 d(p)\ dz CiT£l 1 - 2--~2lh)

K

pz

d(p) dz

(7)

(8) (9)

(10)

(11)

С1Т Ст Ро £

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

e

e

рейнольдсова напряжения (и'ад') привлекаются дифференциальные уравнения переноса [5]:

тт де ТДде ттгде д ^ де д ^ де ^

и0— + У^- + Ж— = тг-Кеу— + — Кег— + Р + С - е, (12)

ду дг ду ду дг

дх

'дг

де

.де

где д де д

де

Цо^ + У^- + Ж— = — К£у— + — Ке^ + С£1-(Р + С) - с£2-

дх

дУ

О г\ еу О 1 О с2 о

дг ду ду дг дг

(13)

Цо ' + V д 7 + Ж х 7 = — Кеу ' +

дх ду дг ду ду

д д('у'ад') е

+ —Кег ' + (1 - С2)Р23 + (1 - Сз)С23 - С1-),

дг дг е

(14)

где коэффициенты турбулентной вязкости определены по упрощенным уравнениям (6) следующим образом:

1 - С2 е(у'2)

К

еу

К

К

еу

С1

еу

а

Ке

/1 \ / /2\ (1 - сз)(1 - С2Т) е2 0 , ' / (1 - С2)е(^'2)----(ад р)

С1Т

е ро

Поэтому

, 1 - (1 - Сз) дШ

1 С1С1Т ро е2 дг /

/ ' л К дЦ К дЦ

) = Кеу) = Кех"ду.

Ке

Ке

а

Величины С1, С2, С3, с1т, С2т, Ст, се1, Се2, а — общепринятые эмпирические константы. Их величины равны 2.2, 0.55, 0.55, 3.2, 0.5, 1.25, 1.45, 1.9, 1.3 соответственно. Выбор этой модели турбулентности обусловлен следующими причинами: модель близка к стандартной (е - е)-модели турбулентности и в ее рамках мы можем учитывать анизотропию турбулентных характеристик в следах в стратифицированной жидкости.

е(у'2) е(ад'2)

Модель 2 подобна модели 1, но Кеу = С3-, Кех = С3-, С3 = 0.25.

ее

Модель 3 основана на локально равновестных соотношениях (9)-(11) и дифференциальных уравнениях для определения (и^Ц) г,^ = 1, 2, 3. Величины («2) (г = 1, 2, 3) определяются из следуюших уравнений:

ио ^ + У^ + Ж? <«">

дх

ду

дг

д К д(м'2) + д К д(и'2) + р + С

— Кеу —--+ 7Т" К ех ^--+ Р11 + С11

ду ду дг дг

- 3 е - С е О«'2)" 3 е) - ЧР» - 2 Р) - 4 - 2 (15)

2

е

д(^'2) т,д(^'2) д (V'2) 5 д (V'2) 5 д (V'2)

+ + ^^^ = ТТ Кеу 7 + 77" Кех 7

дх ду дг ду ду дг дг

- 3е - Се ((V2) - 2^ - Р22 - 2р) - С2 (С22 - 2С) , (16)

г д(и'2) + д(и,'2) + д(и'2) д к д(и'2) + д к д(и'2) + + г

ио—^--+ V—^--+ №—-- = — Кеу—---+ —- Кех—---+ Рзз + Ьзз-

дх ду дх ду ду дх дх

2 £ - -С1 е [(и'2) - 2 е) - с^ Рзз - 2 р) - С^ Сзз - | О), (17)

Кеу = С £ (V'2), Кех = С е {и'2), е = {{и'2) + (V'2) + (и'2))/2;

Р11 = 2^Ку(ду) + К^д!) ) =2Р, О11 = 0, Р22 = 0, О22 = 0, Рзз = 0, п о9 / ' '\ оп Р / ' '\ди^ / ' '\дил ЪТ (диЛ2 I ь- (диЛ2

033 = -27о(шр) = 20' Р =(и""^ +(иш^ = КЛ~ду) + К*Ы) '

О = - ^ (и'р') = дКр. Ш.

Ро Ро дх

Модели 1-3 аналогичны моделям, рассмотренным в [4].

Маршевая переменная х в уравнениях (1)-(4), (12)—(17) играет роль времени. На расстоянии х = хо от тела были заданы следующие начальные условия:

иа(хо,у,х) = @1(т), е(хо,у,х) = @2(г), е(хо,у,х) = @з(г),

2

(¿'и') = Р) = V = № = 0, (и'2) = (V'2) = (и'2) = 2 е,

3

г2 = у2 + х2, 0 < г < ж; -ж < х < ж, -ж < у < ж, х = хо.

Здесь в^г), 02(г) и Оз(г) — функции, согласованные с экспериментальными данными Линя и Пао [6, 7] в однородной жидкости. Граничные условия показаны только для моделей 1 и 2 (для примера).

При г2 = у2 + х2 ^ ж ставились условия невозмущенного потока

Щ = V = № = (р1) = е = е = (V Ш) = 0, х > хо. (18)

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

/ ' л д (Р1) „ д№ дил де де (V'и') = = V = = ^ = ^Т = 0, у = 0, х > 0;

ду ду ду ду ду

. . . . ттг дV ди<1 де де

> = Ы = № = & = аТ = & = дх = х = °' у > 0

При численном решении задачи нулевые краевые условия, соответствующие г ^ ж, сносились на границы достаточно большого прямоугольника 0 < у < у*; 0 < х < х*.

Переменные задачи могут быть обезразмерены с применением масштаба длины О — диаметра тела и масштаба ио — скорости невозмущенного потока. В результате в обезразмеренных уравнениях вместо д появится величина 4П2/Ег2, где Е — плотностное число Фруда, определяемое равенством

Ег= ЩТ т = 1 О' л/ад N

_ _ dps

,Po/ dz '

где T, N — период и частота Вяйсяля—Брента. Для удобства интерпретации результатов расчетов ниже используется время t, связанное с расстоянием от тела следующим соотношением:

жX ^ t ^С

t _ Uo' t _ T _ UoDT _ Fr'

2. Алгоритм решения

Для построения конечно-разностного алгоритма решения задачи, по аналогии с [4], введем новые независимые переменные

ж' _ ж, С _ Xi(y)' V _ X2(z) (ж _ У _ z _ ф2(п))' (19)

Это необходимо, чтобы преобразовать неравномерную сетку в физический области (ж, у, z) в равномерную прямоугольную сетку в вычислительной области (ж', С, п). Уравнения (1)-(5) и (12)—(17) преобразуются согласно (19). Функции ф1 и ф2 устанавливают взаимно однозначное соответствие узлов равномерной сетки в вычислительной области и узлов неравномерной сетки в физической области. Функции ф1(С), ф2(п) задаются таблично. Выбор этих функций позволяет сгущать расположение узлов сетки в окрестности турбулентного следа. Метрические коэффициенты вычисляются с использованием конечных разностей. В вычислительной области (ж',£,п) узлы сетки в (С, п)-плоскости распределены равномерно: С _ ¿'АС, Vj _ j' An, i _ 0,... , M1, j _ 0,... , M2,

^1(CMx ) _ У*,^2(ПМ2 ) _ z*'

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

— скалярные функции, такие как (p1), (p1), Ud, e, e и компоненты рейнольдсовых напряжений ), расположены в узлах основной сетки (в центре ячеек) (^^j) _

(i -АС, j -An);

— значения горизонтальной (V) и вертикальной (W) компонент вектора скорости расположены в узлах сеток, сдвинутых на половину шага в горизонтальном и в вертикальном направлениях соответственно (другими словами, значения горизонтальной и вертикальной компонент вектора скорости отнесены к серединам соответствующих сторон ячеек разностной сетки).

Условимся, что верхний индекс n будет соответствовать сечению ж _ xn _ жп-1 + Ажп, а нижние индексы ¿, j — узлу сетки (Сг,П?) в сечении ж _ const (соответственно нижние индексы i + 1/2, j + 1/2 — узлу сдвинутой сетки Сг+1/2 _ (i + 1/2)АС, ^j+1/2 _ (j + 1/2)An). На рис. 1 приведен пример расчетной сетки и показано расположение неизвестных в узлах сдвинутой сетки.

Алгоритм решения задачи основан на неявном методе расщепления по пространственным переменным для уравнений (1), (4), (12)—(17) и на явном методе расщепления по физическим процессам для системы уравнений (2), (3) и (5). Пусть все неизвестные функции уже определены в сечении жп. Решение на следующем маршевом сечении ж _ жп+1 находится по предложенному алгоритму.

1. Из уравнения (1) получаем значения дефекта осевой скорости [/¿П+1. При этом используется неявная схема расщепления по пространственным переменным [8]. Все остальные неизвестные берутся с предыдущего слоя n.

Щ, ]+1/2

1/2, ]

Щ+1, ]+1/2

1.5

Щ, ] -1/2 ▲ Щ+1, ] -1/2 ▲

V г г + 1

0.5

0 1 2 3 4 5

а б

Рис. 1. Пример расчетной сетки (а) и расположение неизвестных в узлах сдвинутой сетки (б)

2. Компоненты скорости осредненного движения Уга+1, Жга+1 и отклонение давления от гидростатического, (р1)п+1, вычисляются с использованием уравнений (2), (3) и (5). Здесь мы используем идею метода расщепления, в котором вычислительный процесс разделен на три стадии [9, 10]. Первая стадия состоит в вычислении предварительных значений компонент вектора скорости У *, Ж * с помощью явной аппроксимации уравнений (2), (3). Затем, во второй стадии, вычисляем (р1)п+1 как решение уравнения Пуассона с граничными условиями Неймана на осях симметрии и граничными условиями Дирихле на удаленных границах. Метод стабилизирующей поправки [8] используется на этой стадии. Наконец, на третьей стадии новые компоненты вектора скорости Уп+1, Жга+1 вычисляются из требования выполнения условий несжимаемости (5).

3. С применением метода расщепления по пространственным переменным [8] из уравнений (4), (12)—(17) последовательно определяются сеточные функции (р1)П'+1 ,

е?1, еЙГ1, («V)? .

Из соображений простоты реализации алгоритма используется известная идея метода последовательных итераций Гаусса—Зейделя. При вычислении функций (р1)п+1, еп+1, ега+1, (г/ад')га+1 мы использовали величины, уже известные в этом сечении х = жга+1, остальные функции взяты с предыдущего сечения х = хп.

Как пример приведем детали временной и пространственной аппроксимаций, используемых для численного решения уравнений переноса (1), (4), (12)—(17) по схеме расщепления. Рассмотрим типичное уравнение переноса, записанное в новых координатах (19):

/ 1 ад + ^¿т = 1,3 + 3 мт (20)

дх 3 \ д£ дп / 3 \д£ \ уг д£ / дп \ дп / /

д(х, у, х)

Здесь х„, уг обозначают производные по нижним индексам; 3 = ———- — якобиан

д(х',4,п)

преобразования (19); / обозначает одну из неивестных функций Ц^, е, е, (р1) и т.д.; Q — правые части соответствующих уравнений. Для упрощения записи введем обозначения для величин в узлах разнесенных сеток, получаемых с помощью интерполяции: Л.7+1/2 = (Л.7+1 + Д.) х 0.5, /¿+1/2'^- = (/¿+1'^- + Д.) х 0.5. Также введем обозначения

производных от преобразования координат (хп.+1/2

х.+1 -кп

, (уг)

уг+1 - у

¿+1/2'.

к

г

2

1

1+1,]

Схема расщепления для уравнения (20) записывается следующим образом:

(п+1/2 еп , лтп £п+112\ („, лтп еп-

Щ - % + 1 УпГ+1/2)г+1/2,3 - Vп/п+1/2=

т>У ( гП+1/2 _ гП+1/2\ _ К>у ( гП+1/2 _ Яп+1/2\

1 Кг+1/2,3 Я%,3 ) Кг-1/2,3 Яг,3 Я г-1,3 )

I Ь2 + , (21)

'г,з ' *

Я?1 - С1/2 . 1 (у*ЖпГ+1)г,]+1/2 - (у*Шпр+\1-1/2 = ь + I ь

Их °%,3 ,1п

1 Кг ( Яп+1 _ Яп+1\ _ ту-г ( яп+1 _ Яп+1 А

1 К %,3+1/2 \Я%,]+1 Яг,3 ) К %,3-1/2 \3г,3

I Ь2

1г,3 ьп

(22)

х ^ у*

Здесь использованы обозначения Ку = —Ку, Кг = —Кг. Уравнения (21) и (22) реша-

у* хп

ются последовательно с использованием метода прогонки для решения трехдиагональ-ной системы алгебраических уравнений (см., например, [8]).

Уравнения (2), (3) и (5) для поля дефекта давления и компонент вектора скорости в поперечном сечении следа подобны двухмерным несжимаемым уравнениям Навье— Стокса, в которых переменная х играет роль времени. В новой системе координат обез-размеренные уравнения (2), (3) и (5) принимают следующую форму

^ +1 д(ХпV2) +1 д(у*ЖУ) дхп Р) +

5Х + I~дТ + ' дп =--&Г~ + Ы{Р1) ^ ^<* Ш (23)

дЖ + 1 + 1 дЩЖ2) = - д-уф) + ^2((Р1) , е, в, )), (24)

дх I д^ I дп дп

I1 У)+1 ^ <* Ж )=дх. (25)

Здесь Г1, Г2 используются для обозначения правых частей уравнений (2), (3). Частные

диа й

производные —— в (25) аппроксимируются односторонней разностью вперед дх

дилп+1 тпг - т

~ __4 = ап+1

дх ~ Нпх+1 '

Рассмотрим три стадии вычислительного процесса.

— Сначала уравнения (23), (24) решаются без учета градиента давления. Выражения Г1, Г2, стоящие в правых частях, вычисляются по значениям неизвестных на п-м слое по переменной х:

Уг+1/2,3 -УП+1/2,3 + 1 (ХпV2)г+1^ - (ХпУ2)?,3 ,

Ьп+1 '+1/2,3 Ь*

1 (у*ЖУ)п+1/2,3+1/2 - (у*ЖУ)п+1/2,3-1/2

'г+1 /2,3 Ьп

(Юг+гЩ, (26)

Wij+1/2 -Wj+1/2 + 1 (^ПVW)П+Ц/2 j+1/2 - VW)?-1/2j+1/2 +

hx

J,

ij+1/2

1 (y«w2)nnj+1 - (y«w2)

2\ n

Ji,

h«

(F2)

2) ij+1/2 •

(27)

¿3+1/2 '"4

— На втором этапе поле дефекта давления (р1)га+1 вычисляется из разностного аналога уравнения Пуассона

д_ Ж

Zn d К+1) y«

+ —

y« д (p?+ х)

ön

zn

_L + ^ - jq«+.\

hx I dn I

ön

(28)

Как упомянуто выше, уравнение (28) решается с применением итерационной схемы стабилизирующей поправки. В качестве краевых условий для дефекта давления на уда-

ленной границе z = z*, y = y* задаются соотношения (p1)i

0, (P1)

1)N

y j

0. На

осях симметрии у = 0, г = 0 ставились конечно-разностные аналоги условий Неймана, являющиеся следствием дифференциальных уравнений (2), (3) и аппроксимированные следующим образом:

(Р1)2,. -(Р1)1,. п (Р1 )г,2 -(Р1)г,1 „

— 0, —н-- = 0.

h

«

hn

— На последней (третьей) стадии значения компонент вектора скорости Жга+1, Уга+1 определены из требования удовлетворения разностному аналогу уравнения несжимаемости (5)

1

(Zn (P1))n+1j- - (Zn (Р1))

V+1/2=Vi+1/2,j -h"+1 J I h

+ ' ' Ji+1/2j \ h«

n+1

W j1/2 = Wij+1/2 -hn+1

J,

ij+1/2

(y« Ы) j+1 - (y« (p 1))

hn

n+1 ij

+

1

3. Результаты вычислений

С целью проверки эффективности математических моделей и численного алгоритма выполнен ряд численных экспериментов. Вычисления проводились на последовательности сеток и сравнивались с экспериментальными данными Дж. Линя и Ю. Пао [6, 7] по вырождению безымпульсного следа и следа за буксируемым телом в линейно стратифицированной среде. Основные вычисления были выполнены на сетке размерностью 71 х 36 узлов в yz-плоскости. Измельчение размеров ячеек сетки в два раза в окрестности следа приводило к отклонениям величин ^/ёо, Ud0, не превышающим 1-3%. Результаты вычислений для плотностного числа Фруда Fr = 31 представлены на рис. 2 и 3. Изменение во времени обезразмеренных осевых значений энергии турбулентности e0(x) = e(x, 0, 0) и дефекта продольной компоненты скорости UD(x) = Ud(x, 0,0) в безымпульсном турбулентном следе и следе за буксируемым телом также показаны на этих рисунках.

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

1 1 1 1 1 Hassid

- b О g = 0, Lin & Рао "

т i • Fr = 31,Lin&Pao

- g = 0, модель 1 Fr = 31, модель 1

- flL к g = 0, модель 2 -

+ Fr = 31, модель 2

□ g = 0, модель 3

Л Ч >4 ч _ ■ Fr = 31, модель 3

20

40

60

80

x/D

100

ил/и0 1 1 1 1 1

♦ i _ i - - Hassid

0.14 0 g = 0, Lin & Pao .

i • Fr = 31,Lin&Pao g = 0, модель 1 Fr = 31, модель 1

0.1 - Oj « g = 0, модель 2

• V \ + Fr = 31, модель 2

№ 4 W. 4 D g = 0, модель 3

0.06 - ■ Fr = 31, модель 3 _

• - -¥ _ ___

+

0.02 — ч «¿4 -■■■---------- ■J? - n-a 1 " "u 1

20

40

60

80 X/D 100

а

Рис. 2. Сопоставление рассчитанных с применением моделей 1-3 осевых значений энергии турбулентности ео(х) = е(х, 0, 0) (а) и дефекта продольной компоненты скорости По(х) = 0, 0) (б) в турбулентном следе за самодвижущимся телом с экспериментальными данными Дж. Линя и Ю. Пао и результатами расчетов С. Хессида

0.07

4'2/и0

0.05

0.03

0.01

-1-1- -i-Г Hassid ™

о g = 0, Lin & Pao

I • Fr = 31,Lin&Pao

о - g = 0, модель 1

*t\4 — Fr = 31, модель 1 ™

Jb\ Я \ « g = 0, модель 2

+ Fr = 31, модель 2

- D ■ g = 0, модель 3 Fr = 31, модель 3 "

+ t* ~

Цю/и0 0.25

0.15

0.05

Hassid

О g = 0, Lin & Рао1 • Fr = 31,Lin&Pao

0 20 40 60 80 x/D 100 0 50 100 %/D 150

а б

Рис. 3. Сопоставление рассчитанных с применением моделей 1-3 осевых значений энергии турбулентности eo(x) = e(x, 0, 0) (а) и дефекта продольной компоненты скорости Ud(x) = Ud(x, 0, 0) (б) в буксируемом турбулентном следе с экспериментальными данными Дж. Линя и Ю. Пао и результатами расчетов С. Хессида

На рис. 4 представлены линии равной энергии e/e0 = const. Область, занимаемая турбулентным следом за буксируемым телом, значительно больше, чем область турбулентного следа за самодвижущимся телом. Это явление обусловлено тем, что в турбулентном следе за буксируемым телом имеется порождение за счет градиентов осредненной продольной компоненты скорости. В безымпульсном следе вклад порождения за счет градиентов осредненной продольной компоненты скорости несуществен. В следе за самодвижущимся телом вырождение турбулентности происходит быстрее по сравнению со следом за буксируемым телом. В результате, в следе за буксируемым телом

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

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

Рис. 4. Изолинии e/eo = const, eo = e(t, 0, 0) для турбулентных следов за самодвижущимися и буксируемыми телами, Fr = 280, t/T =1, 6

■ • -

< V2 >„ g = о, модель 1 <v'2>0Fr = 280, модель 1 + (w'2)0g = 0, модель 2 * (w'2'\Fr = 280, модель 2 А < V'2 >0 Fr = 280, модель 3 О (V >о модель 3 (avt) A {w'2\Fr = 280,модель 3 а (~w'2\ модель 3 (avt) // « ■ // /0, • . А* л /О X« А+ :

10

<w2>;

ю-4

(и'2)'

10'

10

I lililí

I I I I I 11.

— <»■>,8 = 0' 1 <VFr=280, модель 1

+

x ▲

▲

О

(H'l)0g = 0) модель 2 (i»'2), Fi=280, модель 2 (v l2>. Fr=280, модель 3

<и>' )„Fr=280, модель3 2

{u >0 g = 0, модель 3 - <"'2>o Fi=280,модель3

10'

10

t/T 10

10'

* «i "a

A

к \ % ■

10

t/T 10'

a

Рис. 5. Поведение осевых значений дисперсий турбулентных флуктуаций горизонтальной и вертикальной компонент скорости для безымпульсного следа (а) и следа за буксируемым телом (б); модель 3 (ау^ соответствует результатам расчетов по модели 3 с автомодельными начальными данными, согласованными с экспериментальными данными; модель 2 дает результаты, близкие к модели 1; в безымпульсном следе (и'2) ~ (V2)

х106 х106

а б

Рис. 6. Изменение во времени (t/T = 1.0; 3.0; 5.0) линий (pi)*(y*, z0, t), z0 =2 в безымпульсном следе (а) и в следе за буксируемым телом (б); расчеты проведены по модели 1; модели 2 и 3 дают близкие значения

Динамика турбулентных следов в линейно стратифицированной жидкости характеризуется рис. 6, на котором изображено изменение во времени (t/T = 1.0; 3.0; 5.0) линий (pi)*(y*,Zo,t), = 2 в безымпульсном следе и следе за буксируемым телом. Изолинии избыточного давления носят "волновой" характер. Обращает на себя внимание тот факт, что буксируемому телу соответствуют большие значения (pl)*(y*, z0, t) и это обусловлено большими возмущениями, вносимыми турбулентным следом за буксируемым телом в стратифицированную жидкость. Последнее, как отмечалось выше, связано с наличием в следе за буксируемым телом порождения энергии турбулентности за счет градиентов продольной компоненты скорости.

Заключение

Основные результаты работы могут быть сформулированы следующим образом.

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

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

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

— Рассмотрен вопрос об анизотропном вырождении турбулентности в следе за буксируемым телом в линейно стратифицированной среде.

Автор признателен Г.Г. Черных и А.В. Фоминой за помощь при выполнении работы.

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

[1] Meunier P., Spedding G.R. Stratified propelled wakes //J. Fluid Mech. 2006. Vol. 552. P. 229-256.

[2] Dommermuth D.G., Rottman J.W., Innis G.E., Noyikoy E.A. Numerical simulation of the wake of a towed sphere in a weakly stratified fluid //J. Fluid Mech. 2002. Vol. 473. P. 83-101.

[3] Gourlay M.J., Arendt S.C., Fritts D.C., Werne J. Numerical modelling of initially turbulent wakes with net momentum // Physics of Fluids. 2001. Vol. 13, N. 12. P. 3783-3802.

[4] Chernykh G.G., Voropayeva O.F. Numerical modelling of momentumless turbulent wake dynamics in a linear stratified medium // Comput. and Fluids. 1999. Vol. 28. P. 281-306.

[5] Rodi W. Examples of calculation method for flow and mixing in stratified fluids //J. Geophys. Res. 1987. Vol. 92, N. 5. P. 5305-5328.

[6] Lin J.T., Pao Y.H. Wakes in stratified fluids // Annu. Rev. Fluid Mech. 1979. Vol. 11. P. 317-336.

[7] Hassid S. Collapse of turbulent wakes in stratified media //J. Hydronautics. 1980. Vol. 14. P. 25-32.

[8] ЯнЕнко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.

[9] Белоцерковский О.М. Численное моделирование в механике сплошных сред. М.: Наука, 1977.

[10] Даниленко А.Ю., Костин В.И., Толстых А.И. О неявном алгоритме расчета течений однородной и неоднородной жидкости. М., 1985 (Препринт/ВЦ АН СССР).

[11] Воропаева О.Ф., Мошкин Н.П., Черных Г.Г. Внутренние волны, генерируемые турбулентными следами за буксируемым и самодвижущимся телами в линейно стратифицированной среде // Математ. моделирование. 2000. Т. 12, № 1. C. 77-94.

Поступила в редакцию 2 ноября 2008 г.