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

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

Том 1, № 1, 1996

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ СЕЗОННЫХ ИЗМЕНЕНИЙ ТЕМПЕРАТУРЫ ВОДЫ В ТЕЛЕЦКОМ ОЗЕРЕ*Т

Д. В. Квон

Институт водных и экологических проблем СО РАН Барнаул, Россия

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

1. Введение

Телецкое озеро — проточный водоем вытянутой формы. Длина озера 77.8 км, максимальные ширина 5.2 км, глубина 325 м. В озеро впадают около 70 рек. Основная часть стока в озеро поступает через р.Чулышман (70-75 %) в южном его конце и вытекает в противоположном через р. Бию. Из-за значительной протяженности в продольном направлении и относительно малой ширины основные изменения в сезонной динамике термической структуры озера происходят в продольно-вертикальном направлении. Поэтому при моделировании этого озера здесь используется предложенный в [1] подход с осреднением по ширине озера уравнений импульсов, сохранения массы, переноса тепла и уравнений для энергии турбулентности и скорости ее диссипации. Отметим, что идея осреднения гидродинамических уравнений по ширине для вытянутых водоемов получила затем развитие в ряде работ (см., например, библ. работы [2]). В них для определения коэффициентов турбулентного обмена в реальных расчетах используются эмпирические формулы.

В озере дважды в год в весенне-летний и осенне-зимний периоды наблюдаются термобары, продвигающиеся к средней части озера с обоих его концов в продольном направлении, продолжительностью 1-1.5 месяца.

В настоящей работе описываются математическая модель и результаты расчетов по ней динамики термической стуктуры озера с сопоставлением с данными измерений полей температур в нем. Модель непосредственно учитывает влияние следующих основных боковых притоков озера: рек Чеченек, Кокши, Чири, Большие Чили, Колдор. Влияние остальных,

*© Д. В. Квон, 1996.

^ Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №96-01-01940.

более мелких боковых притоков в модели учитывается добавлением их расходов к расходам основных притоков в соответствии со схемой районирования водосборного бассейна [3]. Модель численно реализована с использованием полунеявной конечно-разностной схемы. Двумерная вертикальная численная модель для Телецкого озера разработана здесь впервые.

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

Осредненные по ширине уравнения импульсов, неразрывности и переноса тепла [1, 2, 4] с учетом притоков имеют вид

дЪщ дЪп1п1 дЪп1п2 ,

+—-------1—-— = -дЪ

т

дхл

дхо

д_

дхл

г +-------I рйх

Ро 3

Х2

+

д дп1 ,

ЪКН------+

дх

дхл дх

д иы дпл , с

ЪКу —----------т + Ьп\,

дхо

дЪщ дЪп2

--------1-------

дхл дх2

д; р = ро[1 - 6.8 ■ 10 ■ (Т - 4)2]

дЪТ дЪщ Т дЪщТ

+ ^т^ + 2

дг

дхл

дх2

д дТ д дТ ЪКнт + т;—ЪКут — + БТ,

дхл

дхл дх2

дх2

(1)

(2)

(3)

где

Бу = дшфш - допьУ, У = (щ,Т), д^а = 0.5(д + |д|), допь = 0.5(д — |д|),

Ъ = Ъл + Ъ2, т = кь1щ1щ

г=1

1 +

дЪг

дхл

+

дЪг

дх2

1/2

г — время; хл и х2 — оси декартовой системы координат, причем ось х2 направлена вертикально вверх; щ и п2 — компоненты скорости по хл и х2 соответственно; д — боковой приток; р и р0 — плотность воды; Т — температура воды; Ъ — “ширина"озера; д — ускорение силы тяжести; г — отметка уровня водной поверхности; Кн и Ку (Кнт и Кут) — коэффициенты суммарной (турбулентной и молекулярной) вязкости (температуропроводности) соответственно в горизонтальном (Н) и вертикальном (V) направлениях; т — сопротивление трения боковой поверхности озера.

Коэффициенты вертикального турбулентного обмена определяются с привлечением уравнений для энергии турбулентности е и скорости ее диссипации е [5]:

дЪе дЪще дЪп2е

дг + дхл + дх2

д ЪК де ЪКне

дЪе дЪще дЪп2е д

дг + дхл

дх2

дх

ЪК

дхл

де

дхл

Не

д

дхл дх

д де

+ —ЪКуе^~ + Ъ(Р - С) - Ъе + Бе, (4)

дхо дхо

де е е2

ЪКует:----+ еле—Ъ(Р- (1 - сзе)С) - С2еЪ—+ Бе. (5)

дх2 е е

Здесь Бу = дгаУга - допЬУопЬ, У = (е,е); К Не (К Не) и Куе (Куе) — коэффициенты суммарной диффузии энергии турбулентности (скорости ее диссипации) соответственно в горизонтальном (Н) и вертикальном (V) направлениях,

Р=К

дщ

дх2

+

дщ

дх2

п ав „ дР С = -д—Ктт, Ро дг

z

1

2

2

2

1

2

2

2

и3 — поперечная составляющая скорости, вычисляемая по одномерной вертикальной модели [6].

Коэффициент турбулентной вязкости К определяется по формуле

е2

К = с(1

Остальные коэффициенты обмена вычисляются по формулам [5, 7]

Ку = v + К, Kvt = vt + ат К, Куе = v + аеК, Ky£ = v + а£К,

где v и vt молекулярные вязкость и коэффициент температуропроводности воды, с^ = 0.09, ат = 0.8, ае = 1.0, а£ = 0.77. Значения остальных констант в уравнениях (4) и

(5) приняты стандартными: с\£ = 1.44, с2£ = 2.0[1, 0 — 0, 3exp (—ReT)], ReT = e2/(ve),

c3£ = 0.8. Суммарные коэффициенты горизонтального обмена Кн, Kht, Кне, Kh£ принимаются постоянными. Их численные значения приводятся ниже при описании численных результатов.

Для системы уравнений (1)-(5) принимаются следующие краевые условия: на входном сечении озера (р.Чулышман) задаются полные потоки импульса, тепла, энергии турбулентности, ее диссипации и расход реки Q(t)

bui — ЬКнUi = fub(t), (bui — ЬКн p = fv(t), Q = QL(t), (6)

где p = (T, e, e), QL(t) — расход реки Чулышман; на выходном сечении озера (р. Бия)

дд bui - bKu) ui = /uR(t), -jxV = 0, Q = QR(t),

где /иь, /ия, 1<р — заданные функций от х2 и 1, QR(t) — расход реки Бии; на водной поверхности при х2 = £

(Т)

дz дz дt + Ul дxl

y

дщ = Tw

дx2 ро

, cpPK

дТ

дx2

Ф,

K

ye

де

дx2

Tw

р0

3/2

є = Се

е3^ У0 '

(S)

где Ф — поток тепла через водную поверхность, — напряжение ветра на водной поверхности, ср — удельная теплоемкость воды; на дне при х2 = 0

д^ дТ

u2 = 0, Ky—— = kblullul,

дx2

дx2.

0, u2 = 0,

де

дx2.

0, є = Се

Уо '

(9)

Здесь у0 и у0 — шероховатости водной поверхности и дна соответственно, кь = 0.14, с£ = 0.314 (см. [8, 9]).

Помимо краевых условий (6) - (9) для системы (1) - (5) необходимо задать начальные условия, соответствующие состоянию покоя и начальному распределению температуры.

є

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

Численное решение поставленной задачи строится на разнесенной C-сетке Аракавы [10]. При этом скалярные величины T, е, е определяются в центральных точках разностной сетки, а компоненты скорости — на гранях элементов этой сетки. Уравнения для р = (T, е,е) записываются в виде

^ = £-Кг^ + Лр”, (10)

Т ОХ2 дх2

где через Лр обозначены слагаемые, аппроксимирующие все оставшиеся члены уравнений, индексы i, к соответствуют точкам конечно-разностной сетки по xi, х2 и n — шагу по времени. Способ аппроксимации источниковых слагаемых в уравнениях (4) и (5) с его тестированием на простом примере течения приведен в [8].

Горизонтальная компонента скорости вычисляется с использованием метода расщепления [11]. На первом дробном шаге в центральных точках производится расчет изменений горизонтального импульса за счет адвекции и диффузии без учета сил давления. При этом используется полунеявная схема (10), но без члена давления:

,,n+1/2 о о n+i/2

- Ui* = AKv + Ли” (11)

Т дх2 ох2

На втором дробном шаге производится пересчет по схеме

n+1/2 z

UU+1/2■k - U“+1/2* = -g(^z”+1 + ± / pn+1 dx)W/2k, (12)

Т ОХ1 Po J

X2

где uci+1/2■ k = 0.5(uci * + uci+1,*). Этот шаг соответствует адаптации полей скорости к распределению давления на гранях элементов разностной сетки. Напомним, что распределе-

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

z

dz д Г

_ +—ju1dx = °, (13)

Z0

в которое подставляется u”+l из уравнения (12). Тогда получим уравнение для определения z на n +1 слое вида

zn+1 - zn д dzn+1

---------= Ъ— Khz^--------+ F, (14)

т дх1 дх1

в котором выражения для Khz и F определяются при указанной выше подстановке. Найденные таким образом значения zn+1 затем подставляются в (12) для вычисления скорости на n +1 слое.

Далее, коэффициенты вертикального турбулентного обмена находятся явно (с известного слоя по времени). Уравнения переноса (10), (11) и уравнение для водной поверхности (14) решаются методом прогонки, при этом полагается, что

U”i, k = °-5(uli+1/2 ,k + Uli-1/2 ,k )• (15)

Значения суммарных коэффициентов горизонтального обмена Kh , Kht , Кне и Кн£ полагались равными между собой и вычислялись по формуле Ричардсона Kh = cl4/3, где параметр l был принят равным шагу сетки в горизонтальном направлении (параметризация процессов подсеточных масштабов). Числовой коэффициент в этой формуле с = 0,01 см2/3/с [12].

Расчеты термической структуры Телецкого озера проведены для промежутка времени от 20 мая до 30 августа 1968 года, данные измерений по которому наиболее полно представлены в [3]. Расчеты были начаты с состояния покоя и с равномерным распределением начальной температуры по глубине со значением, равным 2.3 градуса. Влияние притоков учитывалось заданием только их расходов и температур воды. Динамика расходов задавалась в соответствии с ежемесячным водным балансом озера. Кроме того, в расчетах задавались среднемесячные метеоданные 1968 года, причем влияние ветра учитывалось только при определении турбулентной вязкости, влияние его на продольную составляющую скорости не учитывалось. Результаты расчетов качественно верно отражают наблюдаемую динамику развития вертикально-продольной термической структура Телецкого озера. Как и в наблюдениях, результаты расчетов обнаруживают формирование термобар на обоих концах озера и их продвижение к центральной части озера с последующим слиянием в средней части озера. На рис. 1 представлены результаты расчетов температур в виде изотерм и полей скорости на 29 июня. Справа втекает р. Чулышман, порождая термобарическую циркуляцию своими более теплыми водами. Слева вытекает из озера р. Бия, создавая стоковое течение в относительно мелководной части озера. Несмотря на стоковое течение, выносящее воды и тепло из озера, фронт тепла перемещается в обратную сторону, к средней части озера. Перед фронтом термобара имеет место вертикальная циркуляция, охватывающая почти всю толщу воды.

Результаты расчетов, представленные на рис. 2, соответствуют 29 августа. Распределение температуры почти однородно по длине озера. Изотерма с температурой 4 градуса достигла глубины 100 - 120 м, вертикальные циркуляции исчезли, имеет место только стоковое течение. Изотермы, представленные на рис. 1 и рис. 2, качественно правильно отражают соответствующие данные измерений, приведенные в работе [3]. На рис. 3 даются как вычисленные, так и измеренные [3] вертикальные распределения температур на рейдовой вертикали №26 (р. Б.Корбу), соответствующие 20 июня, 20 июля и 20 августа.

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

[1] Vasiliev O. F., KvON V. I., Chernyshova R. T. Mathematical modelling of the thermal pollution of a water body. In “Proc. XV IAHR Congress”, Istanbul, 2, 1973.

[2] Архипов Б. А., СОЛБАКОВ В. В. Расчет термогидродинамического режима водоема по двумерной модели. Известия РАН, Физика атмосферы и океана, 30, №5, 1994, 671-685.

[3] Селегей В. В., Селегей Т. С. Телецкое озеро. Гидрометеоиздат, Л., 1978.

[4] ВАСИЛЬЕВ О. Ф., Квон В. И. О теоретическом описании гидротермических явлений в водоемах-охладителях. В “Проблемы теплофизики и физической гидродинамики”, Наука, Новосибирск, 1974, 100-111.

[5] РОДИ В. Модели турбулентности окружающей среды. В “Методы расчета турбулентных течений”, Мир, М., 1984, 227-322.

[6] Квон В. И., Игнатова Г. Ш. Одномерная модель сезонного термоклина в озерах. Водные ресурсы, 6, 1979.

[7] JONES W. P., Launder B. E. The calculation of low Reynolds-number phenomena with a two-equation model of turbulence. Int. J. Heat Mass. Transfer, 16, №6, 1973, 1119-1130.

[8] Игнатова Г. Ш., Квон В. И. О гидродинамической схеме скольжения при турбулентном течении. Метеорология и гидрология, №7, 1987, 50-54.

[9] Квон В. И. Об условиях скольжения на дне водотока. Численные методы механики сплошной среды, 6, №3, 1975, 38-44.

[10] BACKHAUS J. O. A semi-implicit scheme for the shallow water equations for application to shelf sea modelling. Continental Shelf Research, 2, №4, 1983, 243-254.

[11] МАРЧУК Г. И. Численное решение задач динамики атмосферы и океана. Гидроме-теоиздат, Л., 1974.

[12] Озмидов Р. В. Горизонтальная турбулентность и горизонтальный обмен в океане. Наука, М., 1968.

Поступила в редакцию 8 июля 1996 г.

Рис. 1. Вычисленные изотермы и поле скорости в Телецком озере (продольное сечение) на 29 июня 1968 г.

Рис. 2. Вычисленные изотермы и поле скорости в Телецком озере (продольное сечение) на 29 августа 1968 г. Над полем скорости помечено местоположение рейдовой вертикали № 26.

Рис. 3. Измеренные (□) и вычисленные (—) вертикальные распределения температур на рейдовой вертикали № 26 (см. поле скорости рис. 2).