CC BY
89
2014 Математика и механика № 5(31)
УДК 519.6:532.516
Б.О. Цыденов, А.В. Старченко
ПРИМЕНЕНИЕ ДВУХПАРАМЕТРИЧЕСКОЙ А-ю-МОДЕЛИ ТУРБУЛЕНТНОСТИ ДЛЯ ИССЛЕДОВАНИЯ ЯВЛЕНИЯ ТЕРМОБАРА1
При математическом моделировании явления термобара впервые применена двухпараметрическая дифференциальная ¿—»-модель для расчёта значений коэффициентов турбулентной вязкости и диффузии. Проведен сравнительный анализ алгебраической модели Холланда П.Р. и др. и двухпарамет-рической дифференциальной ¿-»-модели Уилкокса на примере воспроизведения термобара в канадском озере Камлупс.
Ключевые слова: термобар, модель турбулентности, математическое моделирование, численный эксперимент, озеро Камлупс.
1. Введение
С 1883 г. (классические работы О. Рейнольдса) многие математики и механики занимаются поиском наиболее приемлемых для практики моделей турбулентности. При численном воспроизведении термобара (термобар - это узкая вертикальная зона в озёрах умеренных широт, разделяющая водоём на области с разными значениями теплофизических характеристик воды) часто используются простейшие полуэмпирические модели турбулентности. В частности, многие исследователи явления термобара (например, авторы работ [1, 2]) коэффициенты турбулентной диффузии задают в виде константы.
Из семейства двухпараметрических моделей турбулентности в океанологической практике наиболее известны:
• «классическая» модель Меллора - Ямада (1982) [3], которая состоит из уравнений для кинетической энергии турбулентности (¿) и произведения ¿1 (I - турбулентный масштаб длины);
• ¿-е-модель Лаундера - Сполдинга с поправками на плавучесть Роди (1987) [4], где е представляет собой скорость диссипации ¿;
• ¿-»-модель Уилкокса (1988) [5], здесь » - частота турбулентных пульсаций.
В работе [6] проведён сравнительный анализ вышеперечисленных моделей
турбулентности для океанических течений. Авторы [6] сделали вывод о том, что ¿-»-модель наиболее реалистично описывает процессы перемешивания в океане. Поэтому можно предположить, что применение двухпараметрической ¿-»-модели турбулентности для численного исследования течений, связанных с термобаром, скорее всего, приведёт к более качественным результатам.
Целью данной работы является сравнение результатов моделирования термобара, полученных с помощью двух различных моделей турбулентности: алгебраической модели Холланда П.Р. и др. [7] и двухпараметрической ¿-»-модели Уилкокса [5].
1 Работа выполнена при финансовой поддержке стипендии Президента РФ для молодых учёных и аспирантов, осуществляющих перспективные научные исследования и разработки по приоритетным направлениям модернизации российской экономики (СП-71.2012.5), и Минобрнауки России (государственное задание № 5.628.2014/К).
2. Математическая модель 2.1. Основные уравнения модели
Негидростатическая математическая модель для воспроизведения гидродинамических процессов в глубоком озере, учитывающая влияние силы Кориолиса, связанной с вращением Земли, и записанная в приближении Буссинеска, включает в себя следующие уравнения:
а) уравнения количества движения
ди ди2 дuw 1 др д ( „ ди Л д (ди Л
— +-+-=---- +—I Кх— I Кг— 1 + 2ПгУ - 20 м>; (1)
д( дх д2 р0 дх дх \ дх) д2 \ д2)
dv duv dwv di dv Л di dv Л
— +-+-= — I Kx— I Kz — 1 + 2Q xw - 2Qzu; (2)
dt dx dz dx v dx J dz v dz J
dw duw dw2 1 dp д i „ dw Л д i „ dw Л gp „_
— +-+-=---- +—\ Kx— I Kz— + 2Q yu - 2Qxv; (3)
dt dx dz p0 dz dx v dx J dz v dz Jp0
б) уравнение неразрывности
dx dz
в) уравнение энергии
^ +dw = 0; (4)
dT duT , dwT d f ^ dT V d (^ dT\ 1 дЯsol
+-+-= —\ Dx— \+—\ Dz— \+--(5)
dt dx dz dx v dx J dz v dz J p0cp dz
г) уравнения баланса солёности в озере
dS duS dwS d í dSЛ di dSЛ
— +-+-= —\ Dx— \+—\ Dz— \, (6)
dt dx dz dx v dx J dz v dz J
где u, v - горизонтальные компоненты скорости; w - вертикальная компонента скорости; Qx, Q.y и Qz - компоненты вектора угловой скорости вращения Земли; g - ускорение свободного падения; cp - удельная теплоёмкость; T - температура; S - солёность; p - давление; р0 - плотность воды при стандартном атмосферном давлении, температуре TL и солёности SL (TL и SL - характерная температура и солёность озера соответственно). Коротковолновая солнечная радиация, проникающая в воду, рассчитывается по закону
Hsol = Q (1 - г,)exp (-8absz),
где Q - поток тепла через свободную поверхность, r, - коэффициент отражения воды, eabs - коэффициент поглощения.
В качестве уравнения состояния p = p(T,S,p) выбрано уравнение Чена -Миллеро [8], принятое UNESCO. Данное уравнение состояния связывает плотность воды с температурой, солёностью, давлением и справедливо в диапазоне 0 < T < 30 °C, 0 < S < 0.6 г/кг, 0 <p < 180 бар.
2.2. Модель турбулентности
Для замыкания системы уравнений (1) - (6) используются двухпараметриче-ская Л-ю-модель турбулентности Уилкокса [5], состоящая из уравнений для кине-
тической энергии и частоты турбулентных пульсаций:
дк дик дм>к д ( V т дк Л д ( V т дк Л п п
— +-+-= —I —т--1+—I т--1 + Р + В-е; (7)
дt дх дг дх ^ стк дх) дг ^ стк дг )
дю дию д^ю д (Vт дю Л д (Vт дю Л ю , л л ч
17= дХ| 17 1+7"! "Г 1+к(Р + сюзВ-сю2е) (8)
дt дх дг дх I стт дх ) дг I стт дг ) к
где к - кинетическая энергия турбулентности; ю - частота турбулентных пульсаций; Vт = Т - турбулентная кинематическая вязкость. Числовые значения пара-
тю
метров уравнений (7), (8) представлены в табл. 1.
Таблица 1
Параметры Л-ю-модели турбулентности
Модельная константа Сю1 Сю2 Сю3 стк (с,0)2
Значение 0,555 0,833 0,755 (В>0) -1 (В<0) 2 2 0,307
Член е = (с") кю в (7), (8) отвечает за вязкую диссипацию энергии турбулентности, а члены Р и В выражают генерацию турбулентности за счёт сдвиговых напряжений и действия силы плавучести:
Р = vт
+ 2
( ди д^ Л2 (ду
+ | — + — I +| —
\дг дх) \дх
2 + 2
V 2
В =--^ N2,
Ргт
где N2 - частота Брента - Вяйсяля (частота плавучести), определяемая соотношением
N2 = -Х. Ро дг
Коэффициенты вертикальной диффузии количества движения и скалярной характеристики рассчитываются как
К = vт +V;
V V
=+—,
г Ргт Рг
где V - молекулярная кинематическая вязкость воды; Рг, Ргт - молекулярное и турбулентное числа Прандтля, Ргт = 1.
Коэффициенты горизонтальной диффузии задаются постоянными [7]
Кх = = 2,5 м2/с.
2.3. Начальные и граничные условия
В начальный момент времени система находится в состоянии покоя и удовлетворяет заданным полям температуры, солёности и турбулентных характеристик:
t = 0: u = 0; у = 0; w = 0; T = ^; £ = SL; к = 0; ю = юр,.
Граничные условия имеют следующий вид:
а) на поверхности
ди. = 0; ^ = 0; w = 0; = 0; В,дТ = ; = 0; дю = ^
дг дг дг дг р0 • ср дг дг с°к(г - г + г0)
где к = 0,41 - константа фон Кармана, Lz - глубина озера, г0 = 0,5 м - высота шероховатости.
б) на твёрдых границах (на дне)
дТ „дБ „ дк „ дю 4к
u = 0; v = 0; w = 0; — = 0; — = 0, — = 0;
дп дп дп дп с0к( + г0)
где п - направление внешней нормали к области, г0 = 0,05 м;
в) на границе входа реки (на левой границе) (см. рис. 1, б)
и = ик ; у = 0; w = 0; Т = Тк ; Б = SR , к = кК , ю = юй ,
где ия - скорость речного притока, ТК, БК, кК, юд - температура, солёность, кинетическая энергия турбулентности и частота турбулентных пульсаций в реке соответственно. Начальные значения для кинетической энергии турбулентности и частоты турбулентных пульсаций в реке задаются исходя из предположения изотропии и интенсивности турбулентности 30 %.
г) на правой границе
^ = 0; — = 0; ^ = 0; дТ = 0; = 0; = 0; ^ = 0.
дх дх дх дх дх дх дх
Начальное поле давления определяется из решения уравнений состояния и гидростатики с граничным условием на поверхности р = ра (ра - атмосферное давление) методом Рунге-Кутты четвёртого порядка точности. Уравнение гидростатики выводится из уравнения движения для вертикальной компоненты скорости (при условии w = 0 и отсутствии силы Кориолиса) и имеет вид
др
¥ = "РЯ .
3. Численный метод решения
Решение задачи основано на методе конечного объёма, согласно которому скалярные величины (температура, солёность и т.д.) определяются в центре сеточной ячейки, в то время как компоненты вектора скорости - в средних точках на границах ячеек. В целях приближения расчётной области к прибрежному профилю озера применяется метод блокировки фиктивных областей [9]: приравниваются нулю компоненты скорости в выключенной зоне за счёт использования больших значений коэффициентов вязкости в этой зоне.
Численный алгоритм нахождения поля течения и температуры опирается на разностную схему Кранка - Николсон. Конвективные слагаемые в уравнениях аппроксимируются по противопотоковой схеме QUICK [10] второго порядка.
Для согласования рассчитываемых полей скорости и давления разработана процедура SIMPLED для течений с плавучестью, представляющая собой модификацию метода SIMPLE Патанкара [9]. Алгоритм SIMPLED, основанный на циклической последовательности операций «предположение-коррекция», приведен в ранее опубликованной работе [11].
Системы сеточных уравнений на каждом шаге по времени решаются методом нижней релаксации или явным методом Н.И. Булеева [12].
Тестирование численного алгоритма проводилось для случая квадратной каверны при изотермических боковых границах [13]. Кроме того, серия вычислительных экспериментов «зима», «ранняя весна», «середина весны» и «поздняя весна» [11], связанных с периодом существования термобара, демонстрирует полное соответствие гидродинамических сценариев взаимодействия систем «река -озеро», описанных Кармаком Э.К. и др. на основе натурных наблюдений в озере Камлупс [14].
Озеро Камлупс находится на юго-западе Канады (провинция Британская Колумбия) в 340 км северо-восточнее Ванкувера и расположено между 50°26' -50°45' с. ш. и 120°03' - 120°32' з. д. по течению реки Томпсон, имеет вытянутую форму (см. рис. 1, а).
В качестве исследуемой области выбрано вертикальное сечение озера Камлупс, соответствующее направлению впадения р. Томпсон, начало системы координат
4. Результаты моделирования
4.1 Морфометрия исследуемой области и описание численных расчётов
rf -50-
б
и
x
-150-1
0 ^1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
Длина, м
Рис. 1. Морфометрия оз. Камлупс: а - батиметрия оз. Камлупс, б - расчётная область (продольное сечение)
совпадает с устьем реки (см. рис. 1, а). Расчётная область имеет протяженность 10 км и глубину 150 м (см. рис. 1, б). Участок, через который втекает в озеро речной поток, составляет 15 м от поверхности.
В начальный момент времени система находится в состоянии покоя и удовлетворяет заданным полям температуры и солёности. Начальное распределение температуры в озере Камлупс имеет постоянное значение равное 2,4 °C, в то время как температура воды в реке соответствует 3,6 °C и нагревается на 0,2 °C в день. Река Томпсон впадает в озеро со скоростью 0,01 м/с, минерализация воды в озере и в реке составляет 0,1 г/кг. Поток тепла, поступающий на водную гладь, равен 170 Вт/м2, коэффициент поглощения проникающей в воду солнечной радиации - 0,3 м-1. Ветровое трение отсутствует.
Расчётная область (см. рис. 1, б) покрывается равномерной ортогональной сеткой с шагами hx = 25 м и hz = 3 м . Шаг по времени A = 60 с. Вычисления проводились на суперкомпьютере НИ ТГУ «СКИФ Cyberia».
4.2 Сравнение моделей турбулентности
Результаты расчётов показали (см. рис. 2), что при использовании ю-модели начальная стадия развития термобара существенно не отличается от картины течений, полученной с помощью алгебраической модели Холланда П.Р. и др. [7]:
, 0 0004 + 6.... . N . , м2/с если N2 ^ Лг2
Kx = 2,5 м2/с, Kz =•
,2/„ ^ _) 0,0004 + 6-10-7(Ж2) ' м2/с, если N2 >N^n, [0,02 м2/с, если N2 < N2,
Dx = Kx, Dz = Kz,
где N^ = 9,371-10-10 с-2 (пороговое значение для условия устойчивости вводится для того, чтобы избежать больших значений Kz при N2 ^ 0 ),
N2 = gai-^-Г |, Г = -ga(TC + 273,15) - адиабатический градиент температу-
V dz ) cp
ры, a - коэффициент термического расширения.
Расстояние, м Расстояние, м
Рис. 2. Линии тока на 8-е и 12-е расчётные сутки (сверху вниз): слева - ¿—»-модель, справа - модель Холланда П.Р. и др. [7]
Однако с увеличением времени Л-ю-модель в отличие от модели Холлан-да П.Р. и др. описывает более сложную структуру вертикального перемешивания (рис. 3).
-150Н-1-1-1-1-1-1-1-1-г Н-1-1-1-1-1-1-1-1--
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 0 1000 2000 3000 4000 5000 6000 7000 8000 9000
Расстояние, м Расстояние, м
Рис. 3. Линии тока на 32-е и 40-е расчётные сутки (сверху вниз): слева - Л-ю-модель, справа - модель Холланда П.Р. и др. [7]
Из рис. 3 видно, что при вычислении коэффициентов турбулентной диффузии с помощью Л-ю-модели слева от фронта термобара, в теплоактивной области, формируется приповерхностная локальная циркуляция, которая развивается с течением времени преимущественно в горизонтальном масштабе.
По определению термобар возникает в том месте, где температура воды близка температуре максимальной плотности. Рассмотрим функцию Т-Тмв(р,8), где Тмо(р$) - температура максимальной плотности. Тогда нулевая изолиния функции Т-Тмв(р,8) будет соответствовать температуре максимальной плотности (рис. 4). На основании местоположения температуры максимальной плотности можно сделать вывод о том, что возникающая локальная циркуляция влияет на профиль термобара (рис. 4) и на скорость его горизонтального перемещения (рис. 5). Согласно рис. 5, образование термобара происходит приблизительно на
Рис. 4. Изолинии функции Т-Тмв(р,£) на 32-е и 40-е расчётные сутки (сверху вниз): слева - Л-ю-модель, справа - модель Холланда П.Р. и др. [7]
5-е расчётные сутки. Исходя из графиков зависимости температуры максимальной плотности от времени, представленной на рис. 5, можно заключить, что после 20 суток предсказанные координаты положения термобара начинают различаться: при использовании ю-модели фронт термобара продвигается быстрее по сравнению с результатами, полученными по модели Холланда П.Р. и др. На 40-е расчётные сутки расхождение увеличивается почти на 100 м.
8000 |-
ts
6000 -
4000 -
£ 2000 -
Рч
k-ю-модель модель Холланда и др.
15 20 25 Время, сут
Рис. 5. Динамика эволюции термобара: положение температуры максимальной плотности на глубине 4,5 м от поверхности озера
Заключение
С помощью разработанной негидростатической модели, учитывающей нестационарную и пространственную динамику турбулентной структуры, воспроизведены термогидродинамические процессы, связанные с явлением термобара в озере Камлупс. Численные эксперименты показали, что использование k-œ-модели турбулентности даёт более сложную картину эволюции термобара: в теплоактивной области формируется приповерхностная локальная циркуляция, влияющая на скорость горизонтального перемещения фронта термобара от устья притока.
ЛИТЕРАТУРА
1. Овчинникова Т.Э., Бочаров О.Б. Сезонное влияние вод притока на водообмен в глубоком озере в условиях больших уклонов дна // Вычисл. технологии. 2007. Т. 12. № 6. C. 59-72.
2. Killworth P.D., Carmack E.C., Weiss R.F., Matear R. Modeling deep-water renewal in Lake Baikal // Limnol. Oceanogr. 1996. V. 41. No. 7. P. 1521-1538.
3. Mellor G.L., Yamada Т. Development of a turbulence closure model for geophysical fluid problems // Rev. Geophys. Space Phys. 1982. V. 20. No. 4. P. 851-875.
4. Rodi W. Examples of calculation methods for flow and mixing in stratified fluids // J. Geophys. Res. 1987. V. 92. No. C5. P. 5305-5328.
5. Wilcox D.C. Reassessment of the scale-determining equation for advanced turbulence models // AIAA Journal. 1988. V. 26. No. 11. P. 1299-1310.
6. Umlauf L., Burchard H., Hutter K. Extending the k-œ turbulence model towards oceanic applications // Ocean Modelling. 2003. V. 5. P. 195-218.
7. Holland P.R., Kay A, Botte V. Numerical modelling of the thermal bar and its ecological consequences in a river-dominated lake // J. Mar. Sys. 2003. V. 43. No. 1-2. P. 61-81.
8. Chen C.T., Millero F.G. Precise thermodynamic properties for natural waters covering only limnologies range // Limnol. Oceanogr. 1986. V. 31. No. 3. P. 657-662.
9. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости: пер. с англ. / под ред. В.Д. Виоленского. М.: Энергоатомиздат, 1984. 152 с.
10. Leonard B. A stable and accurate convective modeling procedure based on quadratic upstream interpolation // Comput. Meth. Appl. Mech. Eng. 1979. V. 19. No. 1. P. 59-98.
11. Цыденов Б.О., Старченко А.В. Численная модель взаимодействия систем «река - озеро» на примере весеннего термобара в озере Камлупс // Вестник Томского государственного университета. Математика и механика. 2013. № 5(25). С. 102-115.
12. Булеев Н.И. Метод неполной факторизации для решения двумерных и трехмерных разностных уравнений типа диффузии // Журн. вычисл. матем. и матем. физ. 1970. Т. 10. № 4. С. 1042-1044.
13. Цыденов Б.О., Старченко А.В. Численное моделирование эффекта термобара в озере Байкал в период весенне-летнего прогревания // Вестник Томского государственного университета. Математика и механика. 2011. № 1(13). С. 120-130.
14. Carmack E.C., Gray C.B.J., Pharo C.H., Daley R.J. Importance of lake-river interaction on seasonal patterns in the general circulation of Kamloops Lake, British Columbia // Limnol. Oceanogr. 1979. V. 24. No. 4. P. 634-644.
Статья поступила 06.07.2014 г.
Tsydenov B.O., Starchenko A.V. APPLICATION OF THE TWO-PARAMETRIC k-ю TURBULENCE MODEL FOR STUDYING THE THERMAL BAR PHENOMENON
In this paper, the phenomenon of the thermal bar in Kamloops Lake (Canada) is studied with a nonhydrostatic mathematical model. A thermal bar is a narrow zone in a lake in temperate latitudes where maximum-density waters sink from the surface to the bottom. Two different turbulence models are compared: the algebraic model of Holland P. R. et al. [1] and the two-equation k-ю model of Wilcox D.C. [2]. The two-parameter model of turbulence developed by D.C. Wilcox consists of equations for turbulence kinetic energy (k) and specific dissipation rate (ro).The mathematical model which includes the Coriolis force due to the Earth's rotation, is written in the Boussinesq approximation with the continuity, momentum, energy, and salinity equations. The Chen-Millero equation [8], adopted by UNESCO, was taken as the equation of state. The formulated problem is solved by the finite volume method. The numerical algorithm for finding the flow and temperature fields is based on the Crank-Nicholson difference scheme. The convective terms in the equations are approximated by a second-order upstream QUICK scheme [10]. To calculate the velocity and pressure fields, the SIMPLED procedure for buoyant flows [11], which is a modification of the well-known Patankar's SIMPLE method [9], has been developed. The systems of grid equations at each time step are solved by the under-relaxation method or N.I. Buleev's explicit method [12]. The turbulence models were applied to predict the evolution of the spring thermal bar in Kamloops Lake. The numerical experiments have shown that the application of the k-ю turbulence model leads to new effects in the thermal bar evolution.
Keywords: thermal bar, temperature of maximum density, Boussinesq approximation, numerical experiment, Kamloops Lake.
TSYDENOVBair Olegovich (M. Sc., Tomsk State University, Tomsk, Russian Federation) E-mail: btsydenov@gmail.com
STARCHENKO Alexander Vasilievich (Doctor of Physics and Mathematics, Prof., Tomsk State University, Tomsk, Russian Federation) E-mail: starch@math.tsu.ru
REFERENCES
1. Ovchinnikova T.E., Bocharov O.B. Sezonnoe vliyanie vod pritoka na vodoobmen v glubokom ozere v usloviyakh bol'shikh uklonov dna. Vychisl. tekhnologii, 2007, vol. 12, no. 6. C. 59-72. (in Russian)
2. Killworth P.D., Carmack E.C., Weiss R.F., Matear R. Modeling deep-water renewal in Lake Baikal. Limnol. Oceanogr., 1996, vol. 41, no. 7, pp. 1521-1538.
3. Mellor G.L., Yamada T. Development of a turbulence closure model for geophysical fluid problems. Rev. Geophys. Space Phys., 1982, vol. 20, no. 4, pp. 851-875.
4. Rodi W. Examples of calculation methods for flow and mixing in stratified fluids. J. Geophys. Res, 1987, vol. 92, no. C5, pp. 5305-5328.
5. Wilcox D.C. Reassessment of the scale-determining equation for advanced turbulence models. AIAA Journal, 1988, vol. 26, no. 11, pp. 1299-1310.
6. Umlauf L., Burchard H., Hutter K. Extending the k-ro turbulence model towards oceanic applications. Ocean Modelling, 2003, vol. 5, pp. 195-218.
7. Holland P.R., Kay A., Botte V. Numerical modelling of the thermal bar and its ecological consequences in a river-dominated lake. J. Mar. Sys., 2003, vol. 43, no. 1-2, pp. 61-81.
8. Chen C.T., Millero F.G. Precise thermodynamic properties for natural waters covering only limnologies range. Limnol. Oceanogr., 1986, vol. 31, no. 3, pp. 657-662.
9. Patankar S. Chislennye metody resheniya zadach teploobmena i dinamiki zhidkosti. Moskow, Energoatomizdat Publ., 1984, 152 p. (in Russian)
10. Leonard B. A stable and accurate convective modeling procedure based on quadratic upstream interpolation. Comput. Meth. Appl. Mech. Eng., 1979, vol. 19, no. 1, pp. 59-98.
11. Tsydenov B.O., Starchenko A.V. Chislennaya model' vzaimodeystviya sistem "reka - ozero" na primere vesennego termobara v ozere Kamlups. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika, 2013, no. 5(25), pp. 102-115. (in Russian)
12. Buleev N.I. Metod nepolnoy faktorizatsii dlya resheniya dvumernykh i trekhmernykh raznostnykh uravneniy tipa diffuzii. Zhurn. vychisl. matem. i matem. fiz., 1970, vol. 10, no. 4, pp. 1042-1044. (in Russian)
13. Tsydenov B.O., Starchenko A.V. Chislennoe modelirovanie effekta termobara v ozere Baykal v period vesenne-letnego progrevaniya. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika, 2011, no. 1(13), pp. 120-130. (in Russian)
14. Carmack E.C., Gray C.B.J., Pharo C.H., Daley R.J. Importance of lake-river interaction on seasonal patterns in the general circulation of Kamloops Lake, British Columbia. Limnol. Oceanogr, 1979, vol. 24, no. 4, pp. 634-644.
