Simulation of the absorption of gaseous SO2 by fog droplets using a refined interpolation-sectional droplet model

Pekunov Vladimir Viktorovich Doctor of Technical Science Software Engineer, JSC "Informatika" 153000, Russia, Ivanovskaya oblast', g. Ivanovo, ul. Tashkentskaya, 90 pekunov@mail.ru Other publications by this author

Введение

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

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

Существуют «тяжелые» пакеты вычислительной гидродинамики (FlowVision ^{[1, 2]}, FLUENT ^[3], GAS DYNAMICS TOOL, PHOENICS, Star-CD (см., например, ^[4]) и другие ^[5]), которые способны в полной мере решить сформулированную задачу. Однако универсальность данных пакетов означает избыточность расчетной модели и, как следствие, не всегда наилучшую скорость счета. В связи с этим представляется интересным применить специализированную систему моделирования экологических процессов AirEcology-P, расширив используемую в ней компромиссную математическую модель за счет учета дополнительного фактора слияния капель (данный момент нельзя не учитывать при моделировании тумана, в котором отмечаются высокие концентрации мелких и мельчайших капель). Данная система реализует параллельный счет, при этом решающий программный блок генерируется системой порождения программ таким образом, чтобы включить только необходимые для текущего расчета соотношения, что может в достаточной степени снизить временные затраты на моделирование.

Итак, целью данной работы является повышение точности базовой модели системы AirEcology-P при моделировании взаимодействия тумана с выбрасываемым в атмосферу газообразным диоксидом серы. Сформулируем следующие задачи: а) описать базовый компонент многофакторной двухфазной математической модели, следуя, в основном, работе ^[6], б) расширить математическую модель, включив в нее фактор слияния капель в результате соударений, в) проверить адекватность полученной подмодели капельной фазы, г) провести численные эксперименты по моделированию взаимодействия капель тумана с диоксидом серы.

1. Базовая математическая модель

Следуя работе ^[7], в качестве базовой выберем трехмерную двухфазную многокомпонентную модель распространения загрязнений, учитывающую факторы турбулентности, переноса тепла, прямого и диффузного излучения, динамику водяного пара и капель, конденсацию и испарение капель, поглощение газов каплями. Воспользуемся правой системой координат (x₁, x₂, x₃) с вертикальной осью Ox₃.

где U — вектор скорости основной фазы, P — давление, — молекулярная и турбулентная вязкости, T — температура, — теплопроводность, удельная теплоемкость, плотность и термический коэффициент расширения воздуха, . Коэффициенты отражают влияние турбулентности на физические процессы. Запишем модель турбулентности K-E (RNG):

Выбор модели K-E (RNG) был обусловлен ее весьма хорошими показателями ^[6] при сравнении результатов численного эксперимента с опытными данными для стандартной задачи об обтекании препятствия (CEDVAL A-1).

где C_SO2 — концентрация газообразного диоксида серы, W_SO2 и D_SO2 — скорость витания и коэффициент диффузии диоксида серы.

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

где — коэффициент диффузии пара, M_k — молярная масса воды.

Для тумана также характерна достаточно высокая степень поглощения солнечного излучения, также оказывающая влияние на тепловые эффекты (в частности, члены из приведенного выше уравнения распространения тепла отражают тепловой поток в среду при поглощении диффузного и прямого излучения). Поэтому, целесообразно включить в модель уравнения переноса диффузного излучения в N_rad диапазонах (в последнем диапазоне — тепловое инфракрасное излучение), согласно модели K.N.Liou:

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

Для нахождения используется обратная трассировка луча: из каждой ячейки трассируются лучи (в направлении, обратном вектору падения прямого солнечного излучения) до ближайшей ячейки с известной величиной , причем для каждого луча рассчитывается ослабление по закону Бугера-Ламберта, а результаты усредняются (это повышает точность расчета частично затененных ячеек).

Для моделирования фазы капель из Z компонентов воспользуемся базовой моделью ^[8], в которую будут добавлены новые члены , отражающие изменение плотностей и концентраций компонентов в результате соударения капель, которое нельзя не учесть в случае высоких концентраций капель, характерных для тумана. Уравнения капельной фазы имеют вид:

где — плотность i-го компонента, — концентрация i-го компонента, — вектор скорости i-го компонента, — коэффициенты, — сила сопротивления капель потоку (в случае тумана капли достаточно малы и деформацию капель можно не учитывать), — шаг интегрирования по времени, — потоки пара и SO2 между i-м компонентом и средой, — плотность вещества i-го компонента, — концентрация SO2 в каплях i-го компонента, — концентрация пара на поверхности капли (с поправками Кельвина и Кехлера).

Алгоритм расчета параметров кусочно-сплайновых распределений n_i(D) путем интерполяции по значениям был ранее предложен автором в работе ^[9] и успешно работает даже с разрывными распределениями, которые могут недостаточно корректно обрабатываться при использовании традиционных подходов ^{[3, 10, 11]}. Данный алгоритм используется без изменений и здесь не приводится всвязи с большим объемом.

Граничные условия, в целом, ставятся, следуя работе ^[6].

Для численного интегрирования уравнений переноса используется локально-одномерное расщепление (неявная разностная схема), для уравнений Пуассона и Гельмгольца — верхняя релаксация (с «шахматным» порядком обхода узлов и чебышевским ускорением). Отметим, что возможность применения метода верхней релаксации для уравнения Гельмгольца (для диффузного излучения) подтверждена на практике в расчетах автора данной работы и авторов модели переноса излучения ^[12].

2. Новые компоненты модели

В базовую математическую модель были введены члены , описывающие изменение плотности и концентрации i-го компонента вследствие возникновения новых капель при столкновениях. Известно ^[13], что число соударений f на единицу объема в секунду для двух множеств сферических частиц с характерными диаметрами d₁ и d₂ и концентрациями n₁ и n₂ равняется

где u₁ и u₂ – характерные скорости множеств. Переходя к принятым в данной работе параметризованным распределениям капель n_i(D), можем записать выражение

характеризующее число образующихся новых капель на единицу объема в секунду, при соударении капель i-го и j-го компонентов (двойной интеграл может быть легко вычислен методом Симпсона). Это число также равняется числу «уходящих» в результате соударения капель (на единицу объема в секунду) как из i-го, так и из j-го компонентов. Соответственно, можно записать, что

где x_i и y_i – начальный и конечный диаметры капель заполненного участка i-го компонента (см. ^{[6, 8]}).

Была проведена проверка адекватности предложенных модифицированных уравнений, описывающих кинетику капельной фазы, в простом эксперименте по моделированию популяции капель, находящихся в одной ячейке. Результаты экспериментов с использованием предложенных уравнений сравнивались с результатами, полученными с применением лагранжевого подхода. Обнаружилось, что изложенный в данной работе метод моделирования капель при низких вычислительных трудозатратах имеет относительную погрешность около 1% с тенденцией к падению погрешности с увеличением числа секций (см. рис. 1) на гауссовском и линейном начальных распределениях капель по диаметрам. Приблизительно такую же величину имела погрешность на «перевернутом» гауссовском распределении, однако здесь увеличение числа секций не дало повышения точности (при таком распределении капли сконцентрированы преимущественно в первом и последнем компонентах, в которых обычно погрешность моделирования наиболее велика, поскольку для них кусочно-сплайновое распределение вырождается в кусочно-линейное ^[9]).

Рис. 1. Зависимость относительной погрешности от числа секций при гауссовском (а), «перевернутом» гауссовском (б) и линейном (в) начальных распределениях капель

3. Численные эксперименты

Моделировался трехмерный участок (см. рис. 2), включающий улицу (Ганновер, Германия) ^[14]. Сплошной линией на рисунке показано положение фронтального сечения, для которого далее будут приведены распределения значений различных переменных. Задавалось начальное распределение температуры от +27 градусов Цельсия сверху до +23 градусов Цельсия снизу. На верхней границе (высота 44 метра) был задан постоянный вектор скорости ветра и постоянные значения интегральной освещенности и направления прямого солнечного излучения. На боковых границах ставились циклические граничные условия для всех переменных.

Рис. 2. Схема расчетной области (толстой сплошной чертой обозначено положение фронтального сечения)

На нижней границе находились источники выделения SO2 (транспорт), рассчитанные по модели MOBILE 6.2 (интенсивность выделения 1,114×10^-14 моль/см³/с). Фоновая концентрация SO2 составляла 200 млн^-1 (8,8×10^-9 моль/см³). Для давления на верхней и нижней границах ставились граничные условия второго рода. Учитывались тепловые эффекты на границах, связанные с поглощением излучения, лучистым теплообменом со стенами зданий – фактически, в поверхностных слоях зданий решалось двумерное уравнение теплопроводности.

В первом эксперименте рассчитывался случай отсутствия тумана, соответственно, начальные концентрации капель во всех компонентах были нулевыми. Полученные за интервал времени моделирования (T_мод = 25 секунд) средние значения концентраций диоксида серы считались базовыми. На рис. 3 показано полученное распределение диоксида серы во фронтальном сечении.

Рис. 3. Распределение концентрации SO2 во фронтальном сечении при отсутствии тумана, t = 25 с. Более темным цветом показаны меньшие значения

В прочих четырех экспериментах задавались начальная повышенная влажность (> 100%) и варьируемая начальная концентрация N₀ (2×10⁶, 5×10⁶, 8×10⁶ и 1,14×10⁷ ед/м³) мельчайших капель тумана (с диаметрами в диапазоне от 10^-6 до 10^-5 м, всего же капельная фаза включала семь компонентов с общей вариацией диаметров от 10^-6 до 5×10^-3 м), начальные концентрации капель прочих диаметров считались нулевыми. На рис. 4 показаны результаты эксперимента (во фронтальном сечении, отмеченном на рис. 2) при N₀ = 1,14×10⁷ ед/м³: распределение диоксида серы (a) в том же фронтальном сечении, распределение температуры (б), распределение интенсивности диффузного теплового инфракрасного излучения (в) и концентрация капельного компонента (г) с диаметрами капель от 6×10^-5 до 10^-4 м. Достаточно хорошо заметен эффект прогрева стен зданий солнечным излучением (вполне согласующийся с картиной распределения диффузного теплового излучения), а также эффект снижения средней концентрации SO2 (в нижней части области) с увеличением высоты. Интересно отметить повышенную концентрацию не оседающих/практически не оседающих капель (с диаметрами от 6×10^-5 до 10^-4 м) в верхней части расчетной области, что, видимо, объясняется более благоприятными условиями для конденсации таких капель в данной части по причине меньшей температуры.

Рис. 4. Распределения различных переменных (концентрации диоксида серы (а), температуры (б), интенсивности диффузного теплового инфракрасного излучения (в), концентрации одного из капельных компонентов (г)) во фронтальном сечении расчетной области. Более темным цветам соответствуют меньшие значения

В экспериментах с туманом средняя (осреднение проводилось по каждому горизонтальному слою при некоторой высоте z) концентрация диоксида серы C_SO2(t = T_мод, N₀, z) показала тенденцию к снижению. Величина убыли

достаточно существенна (особенно в нижней части, в которой находятся источники загрязнителя, а также, при некоторых N₀, в верхней части (см. рис. 5), где возможно формирование повышенной концентрации капель, см. рис. 4, г), что говорит о некоторой самоочистке атмосферы от диоксида серы путем его поглощения каплями тумана, испытывающими конденсационный рост и, при достижении определенных размеров, оседающими на поверхность. Это положительный момент, однако нельзя не сказать и о его обратной стороне: часть поглощенного каплями диоксида серы вступает в химическую реакцию с веществом капель (водой) с образованием сернистой кислоты, достаточно слабой, но все же оказывающей некоторый коррозионный эффект на металлические уличные конструкции (эффект, характерный для влажного лондонского смога ^[6] – в данной работе химические реакции не рассматривались и данный эффект не моделировался).

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

Рис. 5. Профили убыли средней концентрации SO2 (горизонтальная ось) с высотой (вертикальная ось) при N₀ = 2×10⁶ (а), N₀ = 5×10⁶ (б), N₀ = 8×10⁶ (в), N₀ = 1,14×10⁷ (г)

Заключение

Получены следующие основные результаты:

1. Предложена трехмерная двухфазная многокомпонентная модель распространения газообразного диоксида серы в условиях тумана. Модель учитывает наиболее значимые факторы и включает нестандартную подсистему уравнений, описывающих динамику и кинетику капельной фазы, в которую включены новые члены, описывающие слияние капель при соударении. Вычисления проведены в программе параллельного численного моделирования экологических процессов AirEcology-P (зарегистрирована в РОСПАТЕНТ, свидетельство №2006611068, 21.03.2006).

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

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