Содержание к диссертации
Введение
ГЛАВА 1. Основные понятия и эмпирические методы описания атмосферы 23
1. Элементы физики атмосферы и понятие турбулентной диффузии 23
2. Основные эмпирические формулы и параметры 28
ГЛАВА 2. Построение ветрового поля 35
1. Уравнения Навье-Стокса 35
2. Один из эмпирических методов построения ветрового поля 37
3. Зануление дивергенции векторного поля с помощью проецирования на пространство соленоидальных векторов 48
ГЛАВА 3. Моделирование распространения загрязнений в атмосфере 59
1. Транспортно-диффузионное уравнение 5
2. Метод расщепления на процессы 6
3. Сеточно-характеристический метод 72
4. Специальный метод точечных и распределенных частиц 77
ГЛАВА 4. Результаты численного моделирования 88
Заключение 105
Библиография 109
Список публикаций автора по теме диссертации 121
- Основные эмпирические формулы и параметры
- Один из эмпирических методов построения ветрового поля
- Зануление дивергенции векторного поля с помощью проецирования на пространство соленоидальных векторов
- Метод расщепления на процессы
Введение к работе
Задачи, связанные с экологией, выходят на первый план во всех сферах человеческой деятельности, находят особенно широкое применение в народном хозяйстве в связи с усилившейся в последние годы ролью химии в промышленном производстве. Интенсивное социально-экономическое, агротехническое и промышленное развитие оказывают на окружающую среду глобальное воздействие. Проблемы выживания человека требуют конкретных ответов на вопросы о происходящих изменениях в окружающей среде. С увеличением количества автотранспорта постоянно растет суммарный объем выбросов в атмосферу, экологическая обстановка в городах ухудшается. Происходят аварии в химической и нефтехимической промышленности, сопровождающиеся выбросом и распространением облаков горючих и токсичных газов. Для выработки правильных решений по предотвращению или ликвидации чрезвычайных ситуаций необходимо верно представлять динамику их развития.
Решение экологических задач проводится на различных уровнях, в том числе и с помощью компьютерного моделирования. Математическое моделирование является наиболее перспективным направлением решения задач экологии по своим возможностям прогнозирования, а также экономичности материальных затрат и безопасности для человека проводимых прогностических экспериментов. По своей природе задачи экологии и оценки состояния окружающей среды не допускают проведения полномасштабных натурных экспериментов, и математическое моделирование является, по существу, единственным методом для оценки ситуационных рисков, изучения динамики природных и техногенных катастроф и прогнозирования их последствий, получения общей картины экологической ситуации.
Одной из важных проблем, связанных с экологией, является прогнозирование распространения загрязнений в воздушной среде. К настоящему времени в области математического моделирования распространения загрязнений в атмосфере и разработки численных методов для него сложилась ситуация, при которой проводимые в мире работы рассматривают, как правило, отдельные явления, но не охватывают их комплекса. Обширный экспериментальный материал, накопленный в мире по проблемам экологического мониторинга окружающей среды, позволяет строить физические модели, адекватные реальным процессам на качественном уровне, но только с развитием современных вычислительных методов и фундаментальных исследований в этой области стало возможным создание визуально-прогностических моделей, обеспечивающих количественную оценку результатов возможных аварий и степени опасности их для людей. Эти модели базируются на фундаментальных разработках специальных вычислительных алгоритмов для решения определенного класса газодинамических задач. В настоящее время подобные исследования проводятся в ряде научных центров мира (Калифорнийский университет, Международный институт системного анализа в Австрии, Германский национальный исследовательский центр информационных технологий). Однако проблемы, отвечающие в полной мере поставленной задаче, требуют разработки новых математических моделей, базирущихся на законах сохранения вещества и уравнениях газовой динамики.
Для адекватного математического описания процессов, происходящих в атмосфере, требуется решить проблему построения ее физической модели, поскольку она существенным образом влияет на построение поля ветра и на описание переноса, происходящего в воздушной среде. Необходимая справочная информация по этому вопросу содержится в ряде научных работ. Так, в работе [1] исследовано поведение ветра с высотой, составлены эмпирические формулы для нахождения коэффициентов турбулентной
5 диффузии, рассмотрено влияние температурной стратификации на ветер и на распространение примесей в атмосфере, проанализировано влияние рельефа на скорость ветра. В работе [2] даны основные понятия о термодинамике атмосферы, рассмотрено явление турбулентной диффузии, изучено поведение давления и температуры с высотой, составлены уравнения движения воздушных масс, и на основе их проанализировано поведение ветра при различных физических условиях, приведен ряд эмпирических формул для вычисления коэффициентов диффузии. В работе [3] дана общая характеристика атмосферного пограничного слоя, рассмотрен ряд методов его аналитического описания, изучено несколько динамических моделей его поведения. В работе [4] экспериментально исследовано влияние подстилающей поверхности на турбулентность в атмосфере. В работе [6] сделаны некоторые замечания о турбулентной диффузии в атмосфере и приведены аналитические решения простейших диффузионных уравнений, описана методика расчетов выбросов из дымовых труб (эффективная высота подъема и угол наклона дымового факела, максимальное значение приземной концентрации вредных выбросов и т.д.), дан обзор основных химических реакций, оказывающих вредное влияние на окружающую среду и здоровье человека, приводятся таблицы предельно допустимых коэффициентов (ПДК) вредных веществ. В работе [5] предложены эмпирические формулы для расчета коэффициентов турбулентной диффузии, где особую ценность представляет формула расчета коэффициента горизонтальной турбулентной диффузии, более нигде в научной литературе не встречающаяся, а также описан один из способов введения в транспортно-диффузионное уравнение поправки, описывающей процесс влажного осаждения. В работе [10] приводятся основные понятия, используемые при описании воздушной среды, в частности, вводятся определения градиентного, геострофического, антитриптического и эйлерианского ветра, а также объясняется связь числа Ричардсона с атмосферной стратификацией. В [11] кратко рассмотрена
структура ветра, причины образования вихрей, шквалов и порывов ветра в атмосфере, картина огибания препятствий и переваливания через препятствия масс воздуха, природа сил трения в воздухе, а также движение воздуха при криволинейных изобарах. В работах [16], [21] приведено множество таблиц, отражающих связь физических параметров в турбулентной атмосфере: класса стратификации, высоты пограничного и приземного слоев, диапазона скоростей ветра, величины флуктуации направления ветра и т.д. В работе [15] дается строгое математическое определение процессов турбулентной диффузии в воздушной среде с применением интегрального исчисления, тензорной алгебры и теории рядов, и предлагается описание теории турбулентных процессов на основе статистической концепции, а также с точки зрения спектральной теории турбулентности; в работе перечислены фундаментальные понятия, модели и экспериментальные методы, применяемые для изучения теории турбулентности. Здесь же для моделирования турбулентных течений предлагается прямое численное решение уравнений Навье-Стокса. В монографии [17] даны теоретические понятия и формулы на базе статистических методов и интегрального исчисления, связанные с описанием турбулентных процессов, приведены основы теории турбулентности, предложены различные эмпирические расчетные методы для моделирования диффузионных процессов в атмосфере, изучены процессы рассеяния примеси в струе при различных метеоусловиях, изложены результаты натурных опытов. В монографии приводятся положения и инженерные формулы, используемые в нормативных документах. В работе [36] проводится анализ химических превращений в атмосфере с использованием эмпирических формул и таблиц: перечисляются важнейшие химические реакции, указываются скорости процессов, формулы для вычисления изменения концентраций различных веществ в атмосфере, даются примеры мониторинга концентрации загрязнений вредными веществами в различных
7 географических пунктах. В работе [37] рассмотрены процессы трансформации веществ при их переносе в атмосфере на большие и средние расстояния, изложены методы и результаты измерений доли загрязнителей от различных источников, участвующих в дальнем переносе, описываются траекторные и эволюционные модели переноса веществ в атмосфере и дается сравнение результатов расчетов с натурными измерениями. В монографии [50] рассматривается строение пограничного слоя атмосферы при некоторых упрощенных условиях, приведены уравнения, описывающие поведение сжимаемого турбулентного потока и использующие понятие пульсаций различных физических параметров, обсуждаются вопросы, связанные с суточными колебаниями метеопараметров.
Применению физических моделей, описывающих состояние воздушной среды и перенос вещества в ней, к решению конкретных задач, а также построению для этой цели математических методов также уделено внимание во многих научных публикациях. Так, в работах [2], [14], [15] движение воздушных масс описывается с помощью системы дифференциальных уравнений Навье-Стокса. [2] предлагает некоторое упрощение системы уравнений Навье-Стокса, сводящее ее к уравнениям Экмана, описывающим вертикальный профиль ветра. Возможно также решение системы уравнений Навье-Стокса напрямую с помощью различных разностных схем, которое использует на сегодняшний день ряд научных коллективов. Например, в работе [25] предлагается решение системы уравнений Навье-Стокса на крупной сетке для нахождения распределения давления в области, а затем переход к более мелкой сетке для решения исходной системы. Указанные методы не могут удовлетворять основным требованиям для программного продукта, используемого в системах мониторинга: методы описания состояния атмосферы, основанные на непосредственном решении уравнений Навье-Стокса, требуют колоссальных затрат вычислительного времени, делая данные модели недееспособными в
8 чрезвычайных условиях, предлагаемые же обычно упрощения не позволяют корректно описывать конкретные физические условия (наличие сложного рельефа местности, изменчивость метеоусловий, поле ветра над возвышениями и в условиях городской застройки), для которых решается поставленная задача.
Недостатки существующих методов побудили к разработке быстрого и эффективного способа построения ветра над местностью со сложным рельефом, описанного в Гл.2. При этом был использован накопленный в мире богатый опыт по построению эмпирических методов для моделирования поля ветра. В частности, за основу была взята идея многоступенчатой процедуры, состоящей из построения начального приближения и последующих его корректировок, изложенная, например, в [47], [57], которая была развита в процессе написания диссертации с учетом особенностей решаемых задач.
Одним из основных требований к построенному ветровому полю является удовлетворение этого поля уравнению неразрывности, для чего был разработан метод зануления дивергенции векторного поля на основе начального приближения. В мире неоднократно предпринимались попытки решения задачи минимизации дивергенции ветрового поля. Так, в [44] для этой цели предлагался итерационный метод. Затем в [48] этот метод был адаптирован к двумерным мезомасштабным ветровым полям- поле тока внутри пограничного слоя интегрировалось по вертикали, а дивергенция согласовывалась от точки к точке с учетом необходимости поддержания значений ветра на метеостанциях фиксированными. В [45] уменьшение трехмерной дивергенции ветра базируется на учете ошибок данных измерения, особенно тех, которые возрастают с высотой. В работах [42], [49], [55] описана процедура конструирования трехмерных согласованных по массе полей, основанная на решении уравнения множителей Лагранжа, с использованием вариационного подхода [53-54]. Влияние на поле ветра
топографии, шероховатости подстилающей поверхности и температурного профиля учтено в работе [56], где для учета вклада различных процессов в дивергенцию поля применяются эмпирические коэффициенты. Основным недостатком перечисленных методов является сильная зависимость ветрового поля от эмпирических констант. В [57] предлагается итерационный метод минимизации дивергенции с использованием специальных подгоночных скоростей, однако он слабо обоснован математически и не обладает универсальной и быстрой сходимостью. В статье [5] приводится экстраполяционный метод построения двумерного бездивергентного ветрового поля по известным значениям ветра в нескольких точках (где расположены метеостанции), основанный на выражении ветра через градиент скалярного потенциала, удовлетворяющего двумерному уранению Лапласа; этот метод пригоден только при наличии плоской подстилающей поверхности и часто дает решение, не согласующееся с требованиями логики— например, если ветер известен в одной единственной точке, то наилучшим решением поставленной задачи является однородное ветровое поле, тогда как упомянутый метод дает и в этом случае достаточно сложную картину распределения ветровых потоков. Методика решения двумерного уравнения неразрывности, предложенная в диссертации, обеспечивающая строгое выполнение этого уравнения при минимальном отклонении от начального приближения, является уникальной и в литературе не встречается.
Существует также множество работ, в которых отражены разные подходы для описания физических процессов, связанных с распространением загрязнений. Так называемые модели рассеяния описывают шлейф от облака, движущегося в направлении «среднего ветра» и расширяющегося под действием турбулентных вихрей в пограничном слое. Наиболее сильное влияние на шлейф оказывают турбулентные вихри сходного со шлейфом размера. Большинство моделей рассеяния написаны для близких и средних
10 (мезомасштабных) расстояний- от 2 до 2000 км [81]. На таких расстояниях моделирование конвекции с учетом влияния особенностей подстилающей поверхности имеет особое значение. При моделировании на дальних расстояниях особенности подстилающей поверхности не рассматриваются, для таких случаев используются так называемые траекторные модели, основным входным параметром которых служит поле ветра. В таких моделях примесь считается равномерно перемешанной по всей высоте пограничного слоя и движущейся по направлению ветра. Для ближних расстояний необходимо учитывать опускание шлейфа от приподнятого источника к земле за счет конвекции.
Среди возможных подходов к моделированию распространения загрязнений- подход с применением статистических моделей, основанных на функции распределения Гаусса [17], [27], [78]. Этот подход является полуэмпирическим и дает удовлетворительные результаты для ровной подстилающей поверхности в случае однородной турбулентности и однонаправленного потока воздуха. Гауссов подход применим на коротких расстояниях и непригоден в условиях мезомасштаба, описанных выше.
Одно из направлений в моделировании распространения примеси над местностью, имеющей сложный ландшафт, и в условиях промышленной застройки также заключается в использовании моделей распространения субстанций, предназначенных для ровной подстилающей поверхности (Гауссовых моделей), которые модифицируются путем введения эмпирических коэффициентов, учитывающих возможное повышение концентрации в застойных зонах вблизи зданий и сооружений. Такой подход использован, например, в документе ОНД-86. Этот метод рекомендуется для установления нормативов ПДК (предельно допустимых концентраций) в Российской Федерации. В упомянутом документе вводится поправочный коэффициент, зависящий от взаимного расположения источника загрязнения атмосферы и близлежащих зданий. Подход практически эквивалентен
введению понятия эффективной геометрии источника, поскольку застройка, расположенная на удалении от источника, не учитывается. Метод корректировки значений горизонтальной дисперсии при использовании Гауссовых моделей так же, как и в ОНД-86, дает возможность оценить вероятные повышения концентраций вблизи зданий.
Распределение концентрации с(х, у, z, t) загрязнителей, выбрасываемых в атмосферу единичным источником, с использованием подхода, основанного на распределении Гаусса, для нестационарного случая выражается формулой
c(x,y,z,t) =
2а.
2а.
Q exp[-((*-4~H<)2]exp[-(^l]
(2л-) сххауа2
{ехр[-
(г-Я)'
2а.2
] + ехр[-
(z + НУ
2а.
]},
а для стационарного случая
Ґ (* + н/л
(z-H)2
c(x,y,z) =
ехр -
+ ехр
2тшо о
У г
2а.
2а.
(У-УоУ 2а.2
J)
где х, у, z - линейные координаты; t - время; (хо,уо) - координаты основания источника; Q - мощность точечного источника; и - скорость ветра на высоте Н вдоль оси X; ах, ау - горизонтальные дисперсии по различным направлениям; az - вертикальная дисперсия; Н - эффективная высота источника (примеры вычисления, например, приведены в [82] и [98]); и -скорость ветра на высоте 10 м. Различные аналитические формулы для вычисления значений дисперсий при разной атмосферной стабильности приводятся, например, в [16]. В работе [82] приводятся формулы для вычисления дисперсий по Бриггсу для сельской и урбанизированной местности, справедливые на расстояниях от 100 м до 10 км.
Гауссовы модели обладают рядом существенных недостатков: они не могут учитывать локальные особенности рельефа и непостоянство в
12 пространстве и во времени метеопараметров; не описывают источники, работающие ограниченное время; в них используются дисперсионные характеристики, полученные для наземных, а не приподнятых источников; не учитывают вертикальную структуру пограничного слоя. Численные и натурные эксперименты показали [114], что Гауссовы модели могут адекватно описывать концентрации загрязнений только в горизонтальном направлении, а для расчета вертикального профиля они применимы только на очень коротких расстояниях.
При моделировании течений в уличных «каньонах» в [27] учитываются только здания, расположенные вблизи источника. Такие же предпосылки вводятся при решении уравнений термической гидродинамики и так называемых транспортно-диффузионных уравнений [28-34]. Моделирование течений в каньонах на основе решения уравнений термической гидродинамики сопряжено с известными математическими трудностями, а также с принципиальными трудностями для всех моделей - заданием входных параметров: условий на границах (нижняя граница- с потоком транспорта, здания со своим обменом с уличным воздухом; параметры верхней границы зависят от многих метеорологических факторов) и начальных значений, которые, как правило, должны зависеть от времени и, в частности, от метеоусловий. Кроме того, метеорологические модели в условиях больших городов могут иметь свои специфические особенности, например, они могут описывать образование острова тепла над промышленными и жилыми районами. Проблема при решении уравнений заключается и в том, что необходимо задавать коэффициент переноса, зависящий от энергии турбулентных движений, являющейся функцией многих величин. Наиболее простой способ определения этой функции следует из уравнения баланса турбулентной энергии. Адекватность приводимых моделей реальным условиям во многом определяется выбором значений эмпирических констант. Для описания формирования полей
13 концентраций примеси часто используется полуэмпирическое уравнение переноса и диффузии. Так, в работе [30] сделана попытка на основе полуэмпирического уравнения переноса и диффузии примеси получить распределение примеси в отдельных уличных каньонах.
Физическое моделирование в аэродинамических трубах, заключающееся в проведении в них физических экспериментов [35], служит проверкой правильности выбора математических моделей. Эксперименты дают возможность оценить некоторые особенности распределения примеси в условиях застройки для таких метеорологических условий, которые можно с той или иной точностью воспроизвести в аэродинамической трубе. Следует отметить, что в трубах невозможно соблюсти подобие течения по достаточному набору критериев, например, задать число Рейнольдса одновременно с числом Росби. В то же время метод физического моделирования в аэродинамических трубах часто является единственным для определения некоторых необходимых для моделирования параметров и дает возможность сравнения модели с измерениями, например, распределения потоков воздуха по улицам при различных направлениях ветра. Моделирование потоков в аэродинамических трубах использовалось в работах Института гигиены и патологии с участием Института глобального климата и экологии РАН для оценки санитарного состояния некоторых городов, например, Кировочепецка. Построение эмпирических моделей позволяет анализировать результаты натурных экспериментов. Результаты численного моделирования и физического моделирования сопряжены с построением параметрических моделей распределения примеси в уличных каньонах в зависимости от метеоусловий: скорости и направления ветра, температурной стратификации атмосферы, влажности и т.п. В параметрических моделях концентрацию загрязняющего вещества в уличном каньоне представляет как сумму концентраций: Cj, поступающих непосредственно от источников самого каньона (в основном, автотранспорт);
14 Cr от сторонних источников (например, примесь от промышленных предприятий, переносимая над данной местностью); Сг, обусловленных явлением рециркуляции внутри каньона. Таким образом, суммарная концентрация С может быть записана в виде C=C
15 В работе [31] распределение примеси над застройкой моделируется транспортно-диффузионным уравнением:
f дС^
\ fy J
\
\
дС,
дС,
+ —
+-\к,
dz\ dz, )
+
дС, _ д{иС) d(vC() d(wC) д
дх j
dt дх ду dz дх,
+ R.+ S.+Ln
где Cj - концентрация і-й компоненты примеси, Rj - скорость генерации і-й
компоненты примеси за счет протекания химической реакции, Sj - мощность
источника і-й компоненты, Lj - скорость генерации і-й компоненты за счет
взаимодействия с поверхностью, u, v и w - компоненты скорости ветра, К и
Kz- коэффициенты диффузии в горизонтальном и вертикальном
направлениях.
Решение транспортно-диффузионного уравнения также требует
быстроты и эффективности. Существующие методы, предполагающие запись
решения транспортно-диффузионного уравнения в виде аналитической
формулы, неприменимы для решения поставленной задачи, поскольку не
отражают всей сложности реальных условий. Например, в [6], [23]
дс приводится аналитическое решение уравнения и — = KAc + QS(r),
дх описывающего картину установившегося распределения концентрации загрязнителя от постоянно действующего точечного источника мощности Q в однородном постоянном горизонтальном ветровом поле со скоростью ветра
Q -—(г-х)
и. Это решение выглядит как с = е 2К , где К-коэффициент
АпКг
турбулентной диффузии, одинаковый по всем направлениям; х - координата
по оси, направление которой совпадает с направлением ветра (начало отсчета
совпадает с источником); г - расстояние от источника. Данная аналитическая
формула является точным решением уравнения, однако в записанном виде
это уравнение не отражает реальной физической картины.
Вообще говоря, моделирование турбулентного переноса аналогично молекулярному, с использованием диффузионных коэффициентов или коэффициентов турбулентной вязкости было предложено в [62] Буссинеском. Им было выдвинуто предположение, что турбулентные потоки связаны со средними градиентами физических величин через коэффициенты, зависящие от свойств потоков. Модели, в которых полный турбулентный поток в атмосфере представляется через средний поток, а локальный перенос физических величин соотносится с их градиентами, описаны также, например, в [68] и [102]. Их называют К-моделями или моделями замыкания 1-го порядка.
О применении для моделирования переноса в атмосфере уравнений Навье-Стокса см. Гл.2 п.1.
При моделировании практически важных турбулентных потоков во избежание трудностей, связанных с большим количеством узлов сетки при численных экспериментах [107], может применяться так называемый метод моделирования крупных вихрей (Large Eddy Simulation, LES), состоящий из явного численного представления крупных и параметризации малых вихрей. Внутри пограничного слоя имеются вихри различных масштабов [95], причем крупные вихри (от 100 м до более 1км) образуются из-за неустойчивости среднего потока, а мелкие (от нескольких см до 100 м) - из-за распада крупных вихрей. При достаточно малых размерах вихри не могут служить переносчиками каких-либо физических характеристик, а лишь диссипируют энергию. Первое применение LES-модели описывается в [72]. LES-модели являются промежуточными между прямым численным моделированием турбулентных потоков и статистической теорией турбулентности, которая использует осреднение искомых физических величин. LES превращается в прямое моделирование при достаточно высоком разрешении. Примеры LES-моделей содержатся в работах [72], [73], [95], [96], [ПО], [122]. Способы генерирования величин сеточного масштаба
17 для LES-моделей описаны в работах [87] и [108]. В [115] моделирование крупных вихрей используется для количественного определения условий образования валиков завихренности на основе исследования конвекции между плоскими пластинами с применением параметризации поверхностного слоя; исследован случай движения пластин. Расчеты показали, что важным параметром является соотношение скорости трения на поверхности к масштабу скорости плавучей конвекции: при нахождении этого соотношения в определенном диапазоне конвекция приобретает вид двумерных валиков. В [111] при широком диапазоне размеров вихрей из-за большого количества узлов сетки в LES осуществлялся расчет среднего потока без детальной информации о мелких вихрях, который показал, что можно рассматривать турбулентность в пограничном слое атмосферы как движение вверх небольшого количества островков тепла (термиков), которые, ударяясь о верхнюю границу пограничного слоя, могут захватывать теплый воздух сверху и вовлекать его в пограничный слой. Вокруг термиков воздух, в основном, медленно опускается.
Существуют также так называемые схемы расчета турбулентности с замыканием 2-го и 3-го порядка. Наиболее важная схема описана в [73], [74], где автор предложил производить явный расчет основной порции турбулентности, а мелкомасштабную турбулентность описывать с помощью аппроксимации замыкания второго порядка. Ввиду того, что схеме требовался большой объем вычислительных ресурсов, были предложены схемы с осреднением турбулентности по ансамблю [58], [93], [121], [123], [124]. Схемы с замыканием 2-го порядка можно найти в работах [75], [77], [92], [112], [117], [120], [124], а схемы с замыканием 3-го порядка- в [58], [63], [67]. В статье [59] используется одномерная схема с замыканием 2-го порядка, однако она дает достаточно реалистичную картину турбулентности за счет особого внимания к членам, относящимся к перераспределению давления. Использование моделей высоких порядков замыкания не требует
18 знания коэффициентов турбулентной диффузии [119], поскольку для описания турбулентных потоков в этих моделях применяются прогностические уравнения. Вывод этих уравнений таков, что они содержат неизвестные корреляции между флуктуационными частями давления и скорости, диссипации n-х моментов и (п+1)-е моменты. Например, в случае использования уравнений Навье-Стокса уравнения, описывающие среднее состояние, вычитаются из уравнений для действительных состояний, а затем умножаются на флуктуационные части физических величин. Нелинейность уравнений приводит при осреднении полученных уравнений к появлению моментов более высокого порядка. Чтобы избежать возникновения моментов высоких порядков прибегают, к параметризации неизвестных выражений на определенном этапе расчетов.
Еще одним видом моделей турбулентности являются траекторные модели. Траекторию можно определить как путь пассивных частиц, переносимых воздухом [89]. Несмотря на сложность траекторий отдельных частиц, в целом вещество в атмосфере движется в направлении среднего ветра— ветра, осредненного за период много больший, чем временные масштабы отдельных вихрей. В работе [89] предлагается рассчитывать траектории не отдельных частиц, а целых их пакетов. Мелкомасштабная турбулентность учитывается через изменение размеров этих пакетов. При этом значения компонент поля ветра хранятся в узлах трехмерной сетки, в следствие чего для расчета ветра в любой точке изучаемой области требуется процедура интерполяции [88]. Пакетная модель тесна связана с так называемыми puff-моделями, где клубы от непрерывного источника движутся в меняющемся поле ветра. При этом ветровое поле может строиться разными способами [83], [85], а дисперсии для клубов можно определять либо через экстраполяцию кривых Пасквилла-Гиффорда из Гауссовых моделей на большие расстояния, либо по эмпирическим формулам, как это сделано в работах [83], [85]. В [71] задача о вертикальном
19 расплывании клубов решается на основе уравнений диффузии.
Для решения уравнений гидротермодинамики и уравнений баланса концентрации примеси, возникающих при построении моделей распространения загрязнений, использующих замыкание различных порядков и LES-модели, используются методы конечного дифференцирования, спектральные и псевдоспектральные схемы, методы конечных элементов и интерполяционные схемы [60], [101]. Большинство мезомасштабных моделей используют метод конечных разностей, однако авторы работы [65] разработали модель конечных элементов, которая была опробована в мезомасштабном моделировании над местностью со сложным рельефом. Спектральная модель с применением ортогональных криволинейных координат описана в [106]. О преимуществах спектрального подхода по сравнению с конечно-разностным дифференцированием см. [79], [80]. Об использовании спектральной модели при расчетах бризов см. также [94].
Адаптация перечисленных моделей к топографическим неровностям может проводиться по-разному: в [109] предлагается использовать давление в качестве 3-й координаты при отсутствии вертикальных ускорений, в [104] -представить рельеф ступеньками сетки по координатным осям. Возможно также преобразование системы координат так, чтобы подстилающая поверхность стала координатной поверхностью (например, [101]). Модель в работе [106] базируется на конформном преобразовании осей координат, а в [109] используется специальная схема генерации ортогональной сетки для моделирования метеоявлений.
На основе сказанного выше можно сделать вывод, что существующие методы непригодны для моделирования транспортно-диффузионных процессов либо в силу чрезмерного упрощения реальной картины, либо в силу больших временных и вычислительных затрат. Для быстрого и, в то же время, качественного решения транспортно-диффузионного уравнения в
20 диссертации предлагается предварительное расщепление исходного уравнения на процессы: адвекцию, диффузию и физико-химические процессы.
Для решения уравнения адвекции в мировой практике разработаны различные методы. Наиболее простыми являются методы с использованием явных и неявных разностных схем [18]. В этой области также хорошо известен так называемый метод характеристик [40]. Однако этот метод обладает существенным недостатком, не являясь консервативным. Другим способом решения уравнений адвекции могут являться явные схемы с использованием компенсационных поправок. Среди них широко известен FCT (flux-corrected п-апзроі1)-метод, описанный в [7], [8], [13]. Однако он также не обладает консервативностью.
Вместо метода характеристик в диссертации используется сеточно-характеристический метод. Этот метод был предложен в свое время известным ученым А.С. Холодовым, однако приобрел свою окончательную форму и впервые нашел конкретное применение лишь в процессе написания представляемой работы. Сеточно-характеристический метод обладает несомненным преимуществом перед более известным методом характеристик в силу своей консервативности.
Для решения уравнений адвекции в диссертации был разработан также специальный метод частиц, обладающий 2 преимуществами перед сеточно-характеристическим методом: отсутствием численной диффузии и отсутствием необходимости разбиения процесса двумерной адвекции на 2 одномерных процесса вдоль каждой из координатных осей.
Отправной моделью для создания специального метода частиц послужил классический метод частиц в ячейке. Однако, хотя в мире известен целый ряд вычислительных методов, связанных с введением в рассмотрение частиц при моделировании процессов переноса, предложенный специальный метод кардинально отличается от всех существовавших ранее. Например,
21 описанный в работе [12] вариант так называемого метода частиц в ячейках вводит в рассмотрение поле давления; предполагает учет удельной внутренней энергии частиц; частицы в этом методе могут изменять свой размер; интерполяция ветрового поля выполняется иначе, чем предложено в специальном методе частиц; в работе [12] не рассматривается возможное наличие неадвективных процессов. Специальный метод частиц не требует знания поля давления, не учитывает удельную энергию частиц и предполагает у частиц наличие постоянных, нулевых (точечная частица) либо ненулевых (распределенная частица) размеров. В работе [19] рассматривается решение дифференциальных уравнений первого порядка конкретного вида, тогда как специальный метод частиц может применяться для решения транспортно-диффузионного уравнения, являющегося дифференциальным уравнением второго порядка. Метод, описанный в работе [19], использует произвольное фиксированное число частиц, причем сами частицы фигурируют в виде так называемых функций ядра; при пересчетах физических параметров с частиц на разностную сетку и обратно используются интерполяционные функции; функции ядра и интерполяционные функции представлены в довольно общем виде. Не оговаривается также способ интерполяции поля скоростей адвекции для моделирования движения частиц. В специальном методе частиц частицы рассматриваются как конкретные физические объекты, их число может изменяется на каждом временном шаге, зависит от параметров сетки и от распределения в рассматриваемой области искомой физической скалярной величины; в этом методе указан конкретный способ интерполяции поля скоростей адвекции для любой точки рассматриваемой области; перенос рассматриваемой физической величины с частиц на разностную сетку и обратно производится не по интерполяционным формулам, а на основе наглядных соображений, следующих из представления частиц как физических объектов, а также на основе принципа сохранения пропорций
22 между вкладами частиц, находящихся в пределах одной ячейки разностной сетки, в соответствующее этой ячейке значение искомой величины до и после моделирования неадвективных процессов. Описанный в работе [20] метод крупных частиц вообще не предполагает разбиения движущейся субстанции на частицы. В силу вышесказанного, специальный метод частиц обладает рядом преимуществ по сравнению с существовашими ранее методами и не имеет аналогов в мировых разработках.
Решение диффузионной части уравнения переноса вещества выполняется с помощью широко известных неявных методов: метода сопряженных градиентов и прогонки,- однако наличие сложного рельефа потребовало создания особого способа заполнения используемых матриц.
Автор выражает глубокую признательность за помощь в написании диссертации своим научным руководителям, сотрудникам ИММ РАН, доктору физико-математических наук, проф. Тишкину В.Ф и кандидату физико-математических наук Клочковой Л.В., а также сотруднику ИГКЭ РАН, кандидату физико-математических наук Беспалову М.С. за ценные консультации.
Основные эмпирические формулы и параметры
При Ri 0 имеет место неустойчивая стратификация, при Ri=0-безразличная, при Ri 0 - устойчивая. При Ri Rio турбулентность затухает, при Ri Ri0- усиливается (Rio- некоторое критическое значениек числа Ричардсона). Критерий Ричардсона учитывает не только температурную, но и так называемую ветровую стратификацию. В работах [128], [130], [125], [103] дается также связь стратификации с масштабом длины Монина-Обухова, широко применяемым в теории подобия. В [126] приводится соответствие между классами и параметрами устойчивости (внешним параметром устойчивости, параметром Казанского-Монина).
Слой атмосферы, в пределах которого происходит наиболее интенсивный турбулентный перенос, называется пограничным слоем атмосферы или слоем перемешивания [2], [90]. Пограничный слой лежит между подстилающей поверхностью и высотой, величина которой зависит от метеоусловий и составляет от 250 м до 2 км. Над верхней границей пограничного слоя турбулентная диффузия малоинтенсивна.
Соответствие между классами стабильности и высотой пограничного слоя дается в Таблице 1.3. п.2 данной главы согласно [16]. В [125] для устойчивых и нейтральных состояний приводятся несколько отличные данные. В работе [128] для расчета высоты пограничного слоя приводится эмпирическая формула при высокой устойчивости атмосферы, в [17] -эмпирическая формула при нейтральной атмосфере в холодный период года и при неустойчивой атмосфере в теплое время года. В [103] для формулы расчета, предложенной в [17], подобран эмпирический коэффициент, позволяющий применять ее при более широком диапазоне метеоусловий. О высоте пограничного слоя см. также [127].
Слой атмосферы, в пределах которого наблюдается наибольший рост интенсивности турбулентной диффузии с высотой, называется приземным слоем. Его нижняя граница также совпадает с подстилающей поверхностью, а высота верхней границы лежит в диапазоне от 100 м до 250 м.
Анемометрическая высота- высота мачты, как правило, равная 10 м, на которой с помощью анемометра- прибора для определения ветра -измеряется ветер на метеорологических станциях.
Шероховатостью z0 называется величина, характеризующая среднюю неровность подстилающей поверхности и измеряемая в единицах высоты. Шероховатость оценивается по Таблице 1.2. или определяется соотношением, согласно которому она составляет 1 /7-И/10 средней высоты неровностей на подстилающей поверхности [1].
В предложенных формулах использованы следующие обозначения: v -коэффициент молекулярной вязкости, Ki - значение коэффициента турбулентной диффузии на высоте z\ (обычно при неустойчивой атмосфере Кі«0.1-Ю.2 м /с при Zi=l м), h - высота приземного слоя, к«0.38 - константа Кармана (см. [1], [2]), zo- шероховатость подстилающей поверхности, V-масштаб Монина-Обухова, 5Т - разность температур на высотах z3 и z2, g -ускорение свободного падения, Та- средняя температура воздуха, сэ2-вертикальная составляющая угловой скорости вращения Земли, є - параметр, зависящий от стратификации (є=0 при безразличной, 0 є 0.5 при устойчивой и -0.5 є 0 при неустойчивой атмосфере), и - динамическая скорость, Н -высота верхней границы пограничного слоя, s - целочисленный параметр, характеризующий устойчивость атмосферы, а=4, b=1.3, с=0.85, d=0.217 эмпирические константы. Наиболее простой и удобной является формула (1.2), которая и используется в дальнейшем в вычислениях. Значения h и Kz(h) (см. [16], [70], [82]) в ней зависят от стратификации и определяются по Таблице 1.3. Различные формулы для вычисления коэффициента вертикальной турбулентной диффузии содержатся также, например, в [16], [21], [64], [75], [84], [86], [105] (взята из работы [129]), [113]. Среди них присутствует формула, в которой искомый коэффициент задается постоянным, однако хотя она и обеспечивает легкое аналитическое решение поставленной задачи, но не отражает реальной физической картины турбулентных процессов [76].
Значения коэффициента горизонтальной турбулентной диффузии вычисляются по эмпирической формуле К = a-crQ -max(0.5,w)-#, где а весовой коэффициент, а - угол горизонтальной флуктуации направления ветра в радианах, зависящий от стратификации и определяемый из Таблицы 1.3, w - скорость ветра в м/с, Н- высота пограничного слоя в м. Предложенная формула является модификацией выражения для расчета коэффициента горизонтальной турбулентной диффузии, приведенного в работе [5].
Один из эмпирических методов построения ветрового поля
Существуют модели [21], в которых, помимо описанных выше уравнений, вводятся в рассмотрение уравнение баланса кинетической энергии, уравнение состояния, уравнение сохранения энтальпии, уравнение, описывающее средние движения. Воздух при этом может считаться сжимаемым, либо несжимаемым, но его плотность зависит от температуры. При решении упомянутых уравнений используются такие понятия из теории подобия, как масштаб температуры, масштаб Монина-Обухова, масштаб толщины пограничного слоя, число Россби [21], [100]. Определение в приземном слое скоростей ветра, а также температуры и влажности на основе теории подобия опирается на уравнения для градиентов температуры, скорости и влажности [100]. При их решении возникает ряд дополнительных задач: определение скорости трения, потенциальной температуры, масштаба Монина-Обухова без учета и с учетом влажности [61], параметров устойчивости. Интегрирование уравнений для градиентов скорости, влажности и температуры дает выражения для их вертикального профиля, использующие так называемые универсальные функции устойчивости, определяемые из экспериментов и отличающиеся у разных авторов [24], [41], [46], [69], [118], [130], [131]. Однако в представляемой работе подобные модели не рассматриваются с целью упрощения решения исходной задачи.
Процедура решения системы уравнений (2.1), (2.2) в общем виде с учетом всех факторов также достаточно сложна, поэтому в диссертации предложен полуэмпирический метод построения ветрового поля, состоящий из нескольких последовательных шагов. При этом в виду малости вертикальных движений в атмосфере над горизонтальной однородной подстилающей поверхностью [1] вертикальная компонента скорости ветра полагается w=0. Будем полагать, что при наличии небольших возвышений рельефа, а также в случае городской застройки данное условие также выполняется.
Поскольку показания с метеостанций снимаются через определенные промежутки времени (метеоэпизоды), в рассматриваемой модели скорость и направление ветра полагаются постоянными в течение каждого метеоэпизода.
При математическом моделировании распространения загрязнений в атмосфере ее удобно разбивать на несколько слоев, между которыми происходит обмен по вертикали благодаря турбулентной диффузии. Это разбиение может производиться двумя различными способами. Если рельеф достаточно гладкий, без резких перепадов высот, то слои повторяют рельеф, приподнимаясь над возвышениями, а сами возвышения не принадлежат трехмерной рассматриваемой области моделирования (Рис.2.1а). Если же перепады высот резки, например, в условиях городской застройки, то верхние и нижние границы слоев располагаются строго горизонтально, а сами слои как бы разрываются возвышениями (Рис.2.16, в, г), при этом возвышения включаются в трехмерную область моделирования, а горизонтальная поверхность, относительно которой ведется отсчет высоты при моделировании, выбирается так, чтобы она проходила через наиболее низко расположенную точку земной поверхности в рассматриваемой области. Будем в дальнейшем для краткости говорить, что слои могут или огибать рельеф (1-й вариант), или разрываться рельефом (2-й вариант). Толщины слоев могут как сохранять свои значения в пространстве (Рис.2.1а, б), так и изменять их (Рис.2.1 в, г). Кроме того, верхняя граница пограничного слоя может выполнять различную роль. На Рис.2.1 в изображен случай, когда верхняя граница пограничного слоя (обозначена пунктиром) не связана с верхней границей расчетной области: она лежит внутри области, а высота ее расположения Н оказывает влияние только на характер изменения коэффициента вертикальной турбулентной диффузии. Для случая, изображенного на Рис.2.1 г, верхняя граница пограничного слоя играет также роль верхней границы расчетной области. Поскольку при этом высота ее расположения Н влияет на толщины атмосферных слоев, изменение величины
Зануление дивергенции векторного поля с помощью проецирования на пространство соленоидальных векторов
Необходимость выполнения уравнения неразрывности для поля ветра потребовала разработки эффективного метода зануления дивергенции векторного поля на основе начального приближения. В работе предложен следующий алгоритм. Пусть задано начальное приближение в прямоугольной области Q=LxxLy для двумерного векторного поля: у/(х,у) = (u(x,y),v(x,y)). (2.6) На разностной сетке оно имеет вид: w;- = fiW W где i=0...Nb j=0...N2. (2.7) ті І J+l ш— І+1 v" Vi+l/2j+l Uij+1/2 A vi+l/2j Ui+lj+l/2 ij i+lj
Задание компонент скорости ветра на разностной сетке. Построим на разностной сетке двумерное векторное поле (Рис.2.3) 9 = Кч,/2 Ч1/2,Л i=0...Nl5 j=0...N2, (2.8) обладающее следующими свойствами: его разностная дивергенция Ui+lj+l/2 Uij+I/2 , Vi+l/2j+l Vi+l/2j divw\. = равна нулю, а отклонение от Ac Ay начального приближения w. минимально для всех i=0...Nb j=0...N2. С этой целью введем следующее семейство сеточных векторных функций, значения которых будут определяться в узлах сетки: (2.9) (0,0,0),i p[Jj k, pk (0,0,А), і = Р П j = к; і, р = 0..JV,, J,к = O..JV2, где A - некоторая произвольная константа. Запишем разностные роторы для данных функций: Up,k-l/2 Up,k+l/2 F =rotAM = \Ур+1/2,к Vp-l/2,kJ (2.10) где uptk_l/2, upMxm vp+1/2jt, vp_U2Jc - скалярные сеточные функции, значения которых задаются на сторонах ячеек сетки, такие что (2.11) 0,i p\Jj k, Up,k+l 2 = « -A/Ay,i = pC\j = к; і, р = 0.. JV,, j, к = O...N2, и р,к-\/2 Vp+l/2,k р-\П,к 0, / Ф р U j к, A/Ay,i = Pr[j = k;i,p = 0...Nl,j,k = 0...N2, 0, / р U j кf AI Ac, j = pf]j = k;i,p = 0...7 , ,j,k = O...N2, 0,/ pl)j k, -A/Ax,i = pC)j = k;i,p = 0...Nx,j,k = O...N Получившийся элементарный вихрь изображен на Рис. 2.4.
Компоненты элементарного ротора на разностной сетке Линейная оболочка всех элементарных роторов на разностной сетке образует некоторое линейное пространство % в котором данные роторы играют роль базиса. Очевидно, любая их линейная комбинация также является ротором некоторой сеточной функции. Если функция rotA(x,y) принадлежит указанному линейному пространству, то ее разложение по базису к имеет вид: (rot А), = (ЕСДД, где i=0...Nlf 0=1...N2. (2.12) Запишем начальное приближение для искомого двумерного векторного поля в виде: w (х,у) = rot А(х,у) + а(х,у), (2.13) где х(х,у) - некоторое векторное поле. Если спроецировать начальное приближение на описанное выше линейное пространство, то для искомого векторного поля будет справедливо соотношение: у/(х,у) = хоХ\(х,у). (2.14)
Начальное приближение для двумерного векторного поля на разностной сетке можно записать как w . = (rot A\ + a.Jt i=0...Nb j=0...N2 тогда искомое векторное поле запишется в виде (2.15) (2.16)
Скалярно помножим обе части уравнения (2.15) на величину F„ и примем во внимание (2.12): ; = F,(SCAl + Л i=0-Ni, J=0...N2. (2.17)
Для того, чтобы норма функции а, являющейся отклонением w,y от начального приближения w., была минимальной, проецирование начального приближения на линейное пространство элементарных роторов должно быть ортогональным (6ШК), то есть, должно выполняться условие ортогональности F.. второму слагаемому в правой части формулы (2.15) для w -w всех IHJ:
Метод расщепления на процессы
Математическое описание процессов распространения примесей в ветровом поле над местностью со сложным рельефом и в условиях городской застройки выполняется на основе модели, выражающей закон сохранения вещества. Изменение концентрации вещества за счет адвективного и турбулентного переноса и физико-химических процессов в такой модели в случае, если предполагается огибание рельефа атмосферными слоями (см. п.2), описывается для трехмерной расчетной области как дс — + div cw - div(AT grad c) + rc = Q, (3.1) a где с- концентрация субстанции, К- матрица некоторых когнитивных коэффициентов (в большинстве случаев можно пользоваться упрощением дс дс дс К grad с = (Кх — ,К —, Kz —), где Кх=Ку=Кс, Ко - коэффициент дх ду dz горизонтальной турбулентной диффузии, в дальнейшем обозначаемый как К, a Kz- коэффициент вертикальной турбулентной диффузии), w - поле скоростей адвекции, Q - поле эмиссии, г - суммарный коэффициент, характеризующий интенсивность физико-химических процессов: сухого поглощения, влажного осаждения, химических превращений, радиоактивного распада и т.п.
Модель процессов распространения загрязнителей при разрыве атмосферных слоев возвышениями может быть выражена следующей системой уравнений где ht(x,y) - функция, которая описывает рельеф области, причем ht(x, у) = hb, если точка (х, у, 0) принадлежит основанию возвышения высотой hb, и ht(x, у) = 0, если в этой точке возвышения нет. Надо также учитывать отсутствие диффузии внутри возвышений и сквозь их границы, например, с помощью наложения на коэффициенты турбулентной диффузии дополнительных условий: К = fl(x,y,z,t),z hl(x,y)r\z hl(x+ lim є,у)Г)z h,(x- lim є,у) 0, z ht ( , у) U z ht (x + lim є, у) U z ht {x - lim є, у) -++0 -»+0 fl(x,y,z,t),z hl(x,y)C\z hl(x,y+ lim є) П z ht( ,у - lim є) ,= (3.2a) 0,z ht (x,y)Uz ht (x,у + lim e)\Jz ht(x,y- lim є) /2 (x, yy z, 0, z ht (x, y) + lim є 0,z hl (x, y) + lim є К. = e- + -+0 где fj, f2, - функции, которыми описываются коэффициенты турбулентной, диффузии вне возвышений (см. п.2 Главы 2).
При моделировании трехмерная расчетная область разбивается на слои, как описано в п.2. Главы 2. С учетом того, что толщины слоев могут изменяться в пространстве и во времени, транспортно-дифузионное уравнение приобретает вид: — + divA,c,w, -div ./z.gradcJ-D, +r.c, = htQ. (3.3) dt где і- номер слоя, hj- толщина і-го слоя, w, =(ul vi) - двумерное поле скоростей ветра, div и grad - двумерные операторы дивергенции и градиента, соответственно, - члены, описывающие вертикальную диффузию, N3 - количество атмосферных слоев. Все разработанные и использованные в данной работе вычислительные методы будут описаны ниже исключительно в применении к уравнению (3.3), выражающему слоистую модель распространения загрязнений в атмосфере.
Следует отметить, что построение ветрового поля при использовании слоистой модели требует для удовлетворения закону сохранения вещества в каждом слое выполнения уравнения неразрывности при несжимаемом воздухе в виде divwA =0 (вместо divw, =0), где і- номер атмосферного слоя. Зануление двумерной дивергенции величины w,ht осуществляется способом, описанным в п.З Главы 2. Начальные условия имеют вид: c(x,y,z,0) = c(x,y,z).
При этом очевидно, что для модели, в которой предполагается разрыв атмосферных слоев возвышениями, внутри возвышений с=0.
Общий вид записи для граничных условий на боковой границе области: дс ааа3с + Р8а— = Гда3 (3-4) где аш , Рт и ут - некоторые коэффициенты, которые, вообще говоря, могут иметь различные значения в разных точках границы, а также в разные моменты времени, й - нормаль к боковой границе трехмерной области Пз Граничные условия на верхней границе области предполагают беспрепятственное проникновение загрязнителей через нее .dz dz \ ,а f--ol \dz ) , причем на нижней границе - их полное отражение от границы под нижней границей области rdown в данном случае понимается: для модели, предполагающей огибание атмосферными слоями рельефа- земная поверхность; для модели, предполагающей разрыв атмосферных слоев рельефом - горизонтальная плоскость на нулевой высоте.
Источники
Рассматриваемая модель предполагает наличие площадных (расположенных на некотором участке земной поверхности) и точечных (например, заводские трубы) источников. Площадной источник характеризуется интенсивностью эмиссии Qs, измеряемой вг/(м2-с) и численно равной массе вещества, выбрасываемой с единицы площади за единицу времени. Точечный источник характеризуется интенсивностью эмиссии Qe, измеряемой в г/(м -с) и численно равной массе вещества, выбрасываемой в единицу объема за единицу времени трубой высоты he и диаметром de. Высота и диаметр трубы, а также скорость истечения газов vc и их температура Те, влияют на эффективную высоту подъема газового облака.
Расчет эффективной высоты и эффективных координат точечного источника выполняется следующим способом [5]: 1. рассчитывается эффективная скорость истечения газов: v =0.5 he 6 м/с; 2. рассчитывается тепловая эмиссия:
