Трехмерная модель распространения тумана в воздушном бассейне над водохранилищем

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

Том 6, № 1, 2001

ТРЕХМЕРНАЯ МОДЕЛЬ РАСПРОСТРАНЕНИЯ ТУМАНА В ВОЗДУШНОМ БАССЕЙНЕ НАД ВОДОХРАНИЛИЩЕМ *

Е. В. Квон, Г. С. Ривин Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: rivin@ict.nsc.ru

A three-dimensional model is suggested for the solution of the problem on propagation and structure of vapour fogs. The turbulent exchanged is described by (K-e) — turbulence model. The flows of heat and moisture from the reservoir surface are given by complete formulae. Calculations of microclimatic changes are presented based on meteorological data in the vicinity of Haranor hydropower station.

Введение

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

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

Решению этой проблемы посвящен ряд работ. Например, в [4] с помощью двумерной стационарной модели изучалась структура речных туманов. C применением двумерных [5] и трехмерных моделей [4] исследованы вертикальные и горизонтальные масштабы распространения возмущений от локальных тепловых источников (в частности, водоемов-охладителей). В [6] в квазистационарном случае с помощью трехмерной модели изучен режим течений в пограничном слое атмосферы и распределение температурных и влаж-ностных полей над поверхностью сложной структуры. Для описания турбулентного обмена в работах [6-9] исследовалась K-модель турбулентности.

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

* Работа выполнена при частичной финансовой поддержке Гранта поддержки ведущих научных школ №00-15-98543 и Федеральной целевой программы "Интеграция" (проект №274).

© Е. В. Квон, Г. С. Ривин, 2001.

Кроме того, приведены уравнения модели атмосферы с учетом образования тумана, К-е-модель турбулентности, начальные и граничные условия, описываются численный метод решения и результаты четырех расчетов для трехмерной модели с использованием метеоданных для территории Харанорской ГРЭС.

1. Система уравнений

Математическая модель атмосферы с учетом образования тумана описывается с помощью следующей системы уравнений [10-14]: уравнения движения

дщ + дпап1

дЬ

дха

д дщ . . .. оч д

тх3 Кох:+Л^Ащ +1 (п2 - щ2)+9дХ1

в

—

в,„

хз

Н

дп2 дпап2 д ггдп2 . . , д С в . . , ,

2 + -т^-2 = — К^ + ЛьАп2 - !(щ - Щ) + 9— —(1 + 0.61^)^x3;

дХ2 ] ву

хз

дЬ дха дх3 дх3 уравнение неразрывности

дща дха

0;

уравнения тепло- и влагопереноса

д П дпаП дЬ дха

дБ + дпаБ дЬ дха

д К,Л.АП,

дх3 д

дх3

_ _ д Б . „. ^ К,— + Л, АБ;

дх3 дх3

уравнения для энергии турбулентности, ее диссипации и масштаба турбулентно-

сти

де дпае

--1--—

дЬ дха

дх3

К

де

е— + р - С - СКе—,

де дпае

--1--—

дЬ дха

дх3

Ке

е

д де

ТЫ К ^ + ^ — (Р - (1 - С3^)С) - — •

В этих уравнениях и в дальнейшем использованы следующие обозначения: С, Се, С1е, С3е, ат, а4, а,, ае, ае, Н, кь, к^, к% — эмпирические константы; п1, п2, п3 — проекции скорости ветра на оси декартовых координат х1, х2 и х3, соответственно; Щ, Щ — составляющие скорости геострофического ветра; Л^, Л, — коэффициенты турбулентного обмена; Ь — время; 9 — ускорение свободного падения; f — параметр Кориолиса; в — потенциальная температура воздуха; ву — фоновое значение в; д — удельная влажность воздуха; Т — температура воздуха; р — давление воздуха; р — плотность воздуха; е — энергия турбулентности; е — диссипация энергия турбулентности; — — масштаб тур-

булентности, — = Се

33/2

р

дщ

+

дп2

С

ат VI9 д ву дх:.

в(1 + 0.61д); V

дх3 дх3

молекулярная вязкость; VI — коэффициент турбулентной вязкости; дт — упругость насыщенного водяного пара; 8 - водность; — скрытая теплота испарения воды; К — газовая постоянная; ср — удельная теплоемкость воздуха при постоянном давлении; с€ — удельная

2

2

теплоемкость воздуха при постоянном объеме; Н — верхняя граница расчетной области; ^о — уровень шероховатости подстилающей поверхности; П, Б — переменные Матвеева; П0 — значение П на подстилающей поверхности; Б0 — значение Б на подстилающей поверхности; а = 1, 2, 3; и4 = а4^/ёЬ; С2£ = 2.196[1 — 0.3ехр(—Ю.ПбКе^)], И,е4£ = л/ёЬ/и; К = и + и4; Кв = и + авКе = и + ае^; К£ = и + а^; в = Т(1000/р)х, х = (ср — ср)/ср; А = д2/д х" + 52/Ох"2 — оператор Лапласа.

Переменные Матвеева П и Б, входящие в уравнения (4), (5), связаны с потенциальной температурой и удельной влажностью зависимостями

где q/

д при qm при

6 = 0, 6 > 0;

П = в + — qf, Б = qf + 6,

дт =°^6.4 ехр( 17.55Т-—2^

р \ Т — 31.25

Преимущество введения переменных П и Б в том, что в уравнениях для них отсутствуют члены, учитывающие фазовые переходы воды, в отличие от исходных уравнений (для температуры, влажности и водности).

Распределение давления по вертикали определяется гидростатическим законом

р = р0 ехр

хз

ЕТ

Таким образом, необходимо решить систему из семи уравнений с семью неизвестными функциями п\, и2, п3, П, Б, е, е.

Начальные данные. В качестве начальных условий для системы уравнений

полагаем иа

= па

П

4=0

П0

Б

4=0

4=0

= Б0,

= е°

4=0

4=0

= е0, где в правой

части располагаются известные функции.

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

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

де

дхз

П1

хз=Н

хз=Н

0, е

и1, п2

хз=Н

е3/2

П2, Из

= С —, П

хз=Н Н

хз=Н

П, Б

0,

хз=Н

Б,

хз=Н

где е = е

хз=Н

; П1, П2, П, Б являются известными функциями.

Условия на нижней границе. На нижней границе (на уровне шероховатости г0

р

д

е

е

подстилающей поверхности) задаются условия

К

дщ

дх3

кь1п1щ, г = 1, 2; де

дх3

0, е

Се

3/2

¿о

Над сушей потоки тепла и влаги в переменных Матвеева задаются с помощью упрощенных формул

К,

дП

дх3

4(П - По), к

дБ

дх3

= к0(Б - Бо).

В случае, если параметры распространения тумана измеряются (например, при реконструкции электростанции), то можно подобрать эмпирические константы к° и к°. Но чаще бывает, что таких данных, к сожалению, просто нет. Тогда возникает сложность определения потоков тепла и влаги эмпирическим путем, так как испарение с водной поверхности происходит как в результате вынужденной конвекции (вызванной ветром), так и вследствие свободной конвекции (вызванной эффектами плавучести). При отсутствии дополнительных источников тепла вклад свободной конвекции пренебрежимо мал по сравнению с влиянием вынужденной. Однако с ростом температуры воды, вызванным, например, увеличением тепловой нагрузки электростанции, свободная конвекция становится заметной и должна тоже быть принята во внимание. Поэтому потоки тепла и влаги с водной поверхности удобнее задать по полной формуле через потоки тепла Фе и влаги Фс [15], обусловленные соответственно испарением и теплопроводностью и учитывающие влияние свободной и вынужденной конвекции:

К,

дП

дх3

1000 у

где К,

дд

р

1

К, дТ3+К 4

+

—

с

дд

V

дх

К,

3

дБ

дх3

^ д8 ^ дд Кв— + Кв- 4

дх3

дх3

Фе

К,

дТ

1

Фс

дх3 3.6  103 рср дх3 3.6  103 рЬш Таким образом, предложенные уравнения описывают образование тумана парения с теплой водной поверхности в трехмерном случае, используя К-е-модель турбулентности.

2. Численный метод решения

Численное решение для поставленной задачи строится на С-сетке классификации Аракавы (рис. 1). Скалярные сеточные функции П^т, БНт, еНт, еНт определяются в целых узлах сетки, а компоненты скорости п^т и Щ,т — в полуцелых узлах, при этом области определения сеточных функций пи Щ,т различны.

Уравнение для <р = (П^т, БНт, еНт ,еНт} записывается в виде

- )/т=т!)

к

к

кп

(3) , 7 (3) " ^Уг,з,к+1

кк + к

+

к(3) кк+1

к+1

кк3) + к

—Кп

(3) Уг,],к к+1

кк + к

к+1

к

(3) к—1

кк3) + к

-Кп

Уг,3,к

+

к

(3) к

к—1

кк3) + к

-Кп

Уг,],к-1

ФП+1 1

к-1

кк3) + к

к-1

п,

0

к

Рис. 1. Ячейки сетки по горизонтали (С-сетки Аракавы) и вертикали.

где к = ^

пространственные шаги по г-м, ]-м и

к+1/2

Хз

к_1/2

; через

к-м направлениям вдоль осей х1? х2 и х3 соответственно; кк3) = х3 Л^ обозначены слагаемые, аппроксимирующие все оставшиеся члены уравнений; К^ отвечающий ^ коэффициент вертикального турбулентного обмена.

Горизонтальные компоненты скорости п1 и п2 вычисляются с помощью метода дробных шагов, основанного на расщеплении уравнения по физическим процессам [16].

На первом дробном шаге решаются уравнения движения в целых узлах С-сетки Аракавы с использованием полунеявной схемы:

(Ф

п+1/2

— )/т =

2

к

к

к(3) + к(3) кк + кк+1

К

к

(3) к—1

кк3) + к

-К2

(3) ^ ^г.^к к—1

+

к

кк3) + к

к:

+

к

кк+1 к" кк3) + к(3) ^

ф2+У2 _ , п+

¿77,к

к

к+1

кк3) + к

к+1

(3) ^ ^г^к — к—1

./.2+1/2 _ Ф П+1/2 к

кк + кк—1

+ лф;

п

'¿77,к 5

где ф = (п1т,п2т}; = 0.5(п1г+1/2^,к + п1г_1/2,^,к); п2г,^,к = 0-5(п2г,^+1/2,к + ^-1/2^);

0.25(п1.-1/2 .'к + П1.+1/2 .'к + П1.-1/2'.+1к + П1.+1/2'.+1к); П2и+1/2,к = 0.25(и2г._1/2к + п2г ^.+1/2 к + п2.+1 ^._1/1 к + п2.+1 ^.+1/2 к). Таким образом, на первом дробном шаге решается задача адвективного переноса и турбулентной диффузии. Уравнения (10), (11) вычисляются методом прогонки [16].

На втором дробном шаге решается задача адаптации, где определяются скорости в полуцелых узлах сетки по схеме

;п+1

1г+1/2,^',к

п+1/2

1г + 1/2,^',к

д

Н

в2т

т

;"++ , — ^п+1/2 2г,^' + 1/2,к 21,] + 1/2,к

т

= д

(1 + 0.6^2т )^х3

дх

и в2т

хз

=д

Н

д Г в2т„

дх

хз

(1 +0.61q2r )с£х'3

2 + 1

¿+1/2 ,,,к

п+1

¿,,+1/2,к

где = п2т, г>2 = ИГ.

Чтобы справиться с неустойчивостью, вызванной нелинейными эффектами, используется двухшаговый метод, модифицированный в [17]. Эта модификация осуществляется в три этапа с помощью разностной схемы для зависящих от времени задач аналогично

к

к

методу предиктора-корректора. На первом шаге проводится расчет сеточных функций для нового временного уровня ¿п+1 с шагом по времени т обычным способом с помощью разностной схемы. На втором — осредняются функции по времени на интервале [¿п,£п+1]. Затем на третьем этапе по осредненным величинам рассчитываются функции по той же процедуре, как и на первом этапе, но на этот раз начиная с уровня ^п+1/2, переходя к функциям для п +1 с шагом по времени т/2.

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

По модели (1) - (7) были проведены три варианта расчета для одних и тех же метеоусловий февраля [18] (температура воздуха равна -22.9 °С, влажность воздуха — 0.9 мбар, скорость ветра — 3.6 м/с, облачность — 3.7 балла, значения параметров для суши = 0.00147, кП = 0.00192).

Общее количество узлов сетки равно 21 х 21 х 21. Шаг сетки по х1 и х2 принят равным 1 км. В вертикальном направлении (ось х3) используется неравномерная сетка с шагами ^13) = 2 м, = 3 м, ^33) =5 м, Л<3) = 10 м, = 20 м, = ... = ^ = 30 м, = 15 м, = 5 м. Шаг по времени в расчетах т = 60 с.

На рис. 2 показана схема фрагмента области интегрирования в плане (горизонтальная плоскость х3 = 0).

Результаты, полученные в описанных ниже экспериментах, приведены в графическом виде. Все рисунки соответствуют вертикальному сечению АВ при х1 = 15.5 км.

Рис. 2. Фрагмент территории расчетной области с разностной сеткой.

В начальный момент времени задаются вертикальные распределения и\, и2, П и Б, остальные функции полагаются равными нулю.

В первом эксперименте потоки тепла и влаги в переменных Матвеева задавались по упрощенной формуле (8). В уравнениях модели турбулентности (6) и (7) не учитывался адвективный перенос. При этом к00 = 0.0065, кП = 0.00192 для водохранилища Харанорской ГРЭС с учетом воздействия подогретой до 4.6 °С воды. Это соответствует полной тепловой нагрузке: перепаду температуры на конденсаторах ДТ = 10.8 °С, циркуляционному расходу Q = 40.3 м3/с [18]. Таким образом, в данном расчете были приняты к00 = 0.0065 для воды и к00 = 0.00147 для суши, кП = 0.00192 для воды и суши.

Сопоставляя полученные результаты с результатами предыдущего варианта модели [7] (К-модели турбулентности без учета адвективного переноса), можно отметить, что наиболее заметные различия наблюдаются в значении коэффициента турбулентной вязкости, структуре и области распределения тумана. Протяженность распространения тумана, которая является одной из важных характеристик тумана, в случае применения К-е-модели турбулентности (рис. 3) оказалась больше примерно на 18%, чем при использовании К-модели. Исходя из полученных результатов, заметим, что коэффициент вертикальной турбулентной вязкости, являющийся характеристикой турбулентности, при использовании К-модели определяется как линейная функция. При применении К-е-модели коэффициент вертикальной турбулентной вязкости нелинеен (рис. 4). Остальные характеристики климата практически не различаются.

Далее был проведен второй расчет, в котором в отличие от первого в уравнениях модели турбулентности (6) и (7) учитывался адвективный перенос [19]. На рис. 5, а - 9, а представлены результаты второго варианта расчета изменений микроклимата территории под воздействием водохранилища — охладителя Харанорской ГРЭС — с использованием нелинейной К-е-модели турбулентности.

й 205^ -

СО

Н 137

68-

0 7000 14000 х<2 м 21000

Рис. 3. Изолинии водности (г/кг) для К-е-модели турбулентности без адвективного переноса.

§ 205

Рис. 4. Изолинии коэффициента вертикальной турбулентной вязкости (м2/с) для К-е-модели турбулентности без адвективного переноса.

Видно, что при заданных метеоданных (напомним, что скорость ветра равна 3.6 м/с) туман распространяется на несколько километров по направлению ветра от теплового

Рис. 5. Изолинии водности, г/кг: а — при использовании формулы (8) для задания потоков тепла и влаги; б — расчет с усредненными данными; в — вычисления по формуле для потоков тепла и влаги (9).

источника (рис. 5, а). Отметим, что наличие тумана принято считать при водности, превышающей 0.01 г/кг. Влажность воздуха (рис. 6, а) изменяется по направлению ветра на расстоянии более 15 км. На рис. 7, а показано распределение температуры воздуха, на рис. 8, а — распределение турбулентной вязкости. Поскольку отклонение поля скоростей ветра от фонового значения мало, влияние водохранилища на поле скоростей незаметно. Поэтому, чтобы детально рассмотреть влияние водоема-охладителя на скоростную структуру, на рис. 9, а показано поле отклонений скорости ветра от его фонового значения.

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

В третьем расчете (наиболее физически полная модель) потоки тепла и влаги в пере-

Рис. 6. Изолинии удельной влажности, г/кг; здесь и ниже обозначения а - в соответствуют рис. 5.

менных Матвеева с поверхности водоема задавались через потоки тепла и влаги, обусловленные испарением и теплопроводностью по полной формуле (9) и в уравнениях модели турбулентности (6) и (7) учитывался адвективный перенос [20]. При этом потоки тепла и влаги с поверхности суши задавались, как и раньше, с помощью упрощенной формулы (8).

Рис. 7. Изолинии температуры воздуха, °С.

На рис. 5, в - 9, в показаны результаты третьего варианта расчета изменений микроклимата территории под воздействием водохранилища — охладителя Харанорской ГРЭС. Видно, что при заданных метеоданных туман распространяется на несколько километров по направлению ветра (рис. 5, в). Влажность воздуха (рис. 6, в) изменяется по направле-

Рис. 8. Изолинии коэффициента вертикальной турбулентной вязкости, м2/с.

нию ветра на расстоянии более 15 км. На рис. 7, в показано распределение температуры воздуха, на рис. 8, в — распределение турбулентной вязкости, на рис. 9, в — отклонение от фонового значения поля векторов скорости ветра.

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

СО

H 411

342 274 205 137 68 0

4- afr- * * * * -Ji * --a j

аа" ау .у -V -У -V -V -Y -V -j .J j -K

аУ -V -V ay -a/ -V -y -y j j T

V -V -V аУ J/ -y УУ -y yj л ^ j al -Л

аУ -У J/ -y -y -У yj yl j al i >

Уа -y -y -y -У аУ -У yj al aJ i > > -

-а -У -У -J -У У/ ^ al al al i * >

-а -У -J _J J У/ a/ al al al ^ > У

-а -У J J J J yl yl al a> a> al ai * ^

J J J y> yl ^ •J a> il al ai ai ai J

-а J > i al ai I» ^ a> i J

-ч > i yl ^ a» y> a> ai al Jj

Уа уу J J J J J al al aJ -V

уа < < < « < < t M- ч

г < 4 « f f f f f f f f <- <- <r <r t y

f С- Ï- f- <— <— f- f- e

4--- si а> < f à k k k k k k k л J

7000

14000

ж2, M

21000

^ 480

m

Ы 411

342 274 205 137 68 0

6

^ V ^ a\- Ф Ф * -4/ —X— л-

•v -V ^ -V -V -V -V

v -J ^ -J •V -V w -J -J J Л

^ w ^ J -V -V -J -J J > 1

-V -V -J -J -V Jl J J è -

Л -J -J -J ^ J J л -

V -Л -Л J -J -J -J ^ J J J л J

^ ^ -Л > J J al J J л £ Л J

^ п J J J J J J ^ £ J Л J

-a -v л- r г f L J j J J J\ J J

^ г Г f f С L_ -J -J J ^

^ w -f г г f f r r -J -J -J

•У v- L. Г г г r r г г Г r r- -J -V

4a V- L r г- г- Л- r г r r r Г г r

u V. a- r A- Г r Г r f f

i L L- г г г Г F ^ f- f- f- f- f- <a t

У ai ч à» aà» ,2> л чЭ r A .1 л л л Л Л i .J

7000

14000

Х2, M

21000

^ 480 5? 411 342 274 205 137 68 0

aa 'Л ata * ч> •i/ Ф -4/ aia -J -J

-ta -V -у -У аУ аУ -y ai -J

yv ta И. -V -у аУ аУ -y -if -y -J

^ ^ аУ -У аУ -Jj -y -J -J -J

v л- Г Г аУ аУ J ai -J -ai аУ

-a A- ЛГ s f Î. V v -y -V аУ аУ

к F ^ f t. L Уа ay -V аУ -V

-a V r ra РГ PC г r г •Г V -V w

С r aa /V /V ;а fa г- Г г r •w аУ

-y ^ u /a Л> /V- л- Да Л- Л- аа Ла ya aa aa J аУ

-Г ^ ^ л- Л- л- а- aa aa V

-J -Ja -Г- -г a^ ^ Л- л- д- r Г V

-V J ф ^ ла Па IV- Г ^ F t 1

L- -V > -Л /а- F Р ^ ^ i- 4- <

a i - 4 -à Л .1 а ^ 1 Л ?ï Г а ^ i É i É .Î i a ^

7000

14000

ж2, м

I

21000

Рис. 9. Поле отклонений скоростей ветра от фонового значения, м/с

не всегда удается хорошо подобрать.

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

Видно, что чем численно больше потоки тепла и влаги с поверхности воды, тем сильнее меняется микроклимат территории и увеличивается площадь распространения тумана.

Заключение

Для решения задачи о распространении и структуре туманов парения предложена трехмерная гидротермодинамическая модель атмосферы, в которой турбулентный обмен описывается с помощью К-е-модели турбулентности. Проведены расчеты изменения микроклимата по метеоданным вблизи водохранилища Харанорской ГРЭС. Показана зависимость микроклимата от мощности потоков тепла и влаги. Полученные результаты позволяют сделать вывод о применимости разработанной модели микроклимата вблизи подогреваемого водоема с нелинейной К-е-моделью турбулентности и с учетом влияния водоема для достаточно достоверного прогноза распространения и структуры туманов парения.

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

Авторы благодарят В. В. Кириллова за оказанную помощь при поиске данных по Харанорской ГРЭС, В. И. Квона за ценные советы и поддержку при проведении работы.

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

[1] БЕлолипЕцкий В.М., Шокин Ю. И. Математическое моделирование в задачах окружающей среды. Новосибирск: ИНФОЛИО-пресс, 1997.

[2] Марчук Г. И. Математическое моделирование в проблеме окружающей среды. М.: Наука, 1982.

[3] Шокин Ю.И., Ривин Г. С., ХАкимдзянов Г. С., Чубаров Л. Б. Вычислительный эксперимент как инструмент для исследования природных явлений // Вычисл. технологии. 1992. Т. 1, №3. С. 12-33.

[4] Берлянд М. Е. Современные проблемы атмосферной диффузии и загрязнения атмосферы. Л.: Гидрометеоиздат, 1975.

[5] Вагнер Б. Г., Надежина е. Д. Пограничный слой атмосферы в условиях горизонтальной неоднородности. Л.: Гидрометеоиздат, 1979.

[6] Иванова Л. А., Надежина е. Д. Моделирование пограничного слоя на побережье нагретого водоема // Метеорология и гидрология. 1991. №8. С. 49-55.

[7] Квон В. И., Квон Е.В., Филатова Т.Н. Численное моделирование пространственно-временных изменений микроклимата территории вблизи водохранилища — охладителя Харанорской ГРЭС // Тез. докл. Второго Сибирского конгресса по прикл. и индустр. математике (ЩРШМ-96). Новосибирск: Изд-во Института математики СО РАН. 1996. С. 283.

[8] Хворостьянов В. И. Моделирование возмущений микроклимата от водоемов-охладителей АЭС, ТЭС, бьефов ГЭС и искусственного рассеяния туманов // Метеорология и гидрология. 1991. №11. С. 27-35.

[9] Черных Г. Г. Введение в численное моделирование свободных турбулентных течений. Новосибирск: Изд-во НГУ, 1996.

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

[11] Лайтхман Д. Л. Физика пограничного слоя атмосферы. Л.: Гидрометеоиздат. 1961.

[12] Матвеев Л. Т. Динамика облаков. Л.: Гидрометеоиздат. 1981.

[13] Монин А. С., Яглом А. М. Статистическая гидромеханика. Ч. 1. М.: Наука. 1965.

[14] Пененко В. В., Алоян А. Е. Модели и методы для задач охраны окружающей среды. Новосибирск: Наука, 1985.

[15] Ryan P. J., Harleman D.R. F., Stolzenbach K. D. Surface heat loss from cooling ponds // Water Resources Research. 1974. Vol. 10, No. 5. P. 930-938.

[16] Марчук Г. И. Методы вычислительной математики. М.: Наука. 1989.

[17] Duwe K. C., Hewer R. R., Backhaus J. O. Results of a semi-implicit two-step method for the simulation of markedly nonlinear flow in coastal seas // Continental Shelf Research. 1983. Vol. 2, No. 4. P. 255-274.

[18] Гидротермические исследования охлаждающей способности водохранилища Харанор-ской ГРЭС без учета и с учетом брызгательной установки для варианта при оборотной системе охлаждения. Отчет о НИР/ВНИИГ им. Б. Е. Веденеева. №240-6288. СПб, 1993.

[19] Квон Е. В. Моделирование турбулентной конвекции и переноса тепла и влаги в приводном слое атмосферы // Тез. докл. Третьего Сибирского конгресса по прикл. и индустр. математике (INPRIM-98). Новосибирск: Изд-во Института математики СО РАН, 1998. С. 63.

[20] Квон Е. В. Численное моделирование распространения тумана парения вблизи подогреваемого водоема // Мат. XXXVI Междунар. науч. студенч. конф. "Студент и научно-техн. прогресс". Новосибирск: Изд-во НГУ, 1998. С. 54.

Поступила в редакцию 28 марта 2000 г. в переработанном виде —16 июня 2000 г.