Уточненный расчет распределений капель при моделировании атмосферных многофазных сред

При моделировании процессов в атмосфере часто рассматриваются многофазные системы, включающие многокомпонентные несущую и капельную фазы, а также фазу твердых частиц ^[1]. Существенной проблемой является построение адекватного распределения капель (особенно затруднительно моделировать распределения с разрывами, возникающими при межъячеечном переносе капель [задача еще более усложняется, если капли достаточно велики и, как следствие, нельзя пренебречь их собственными скоростями] или резком изменении их концентраций в результате, например, техногенных выбросов). В полной мере данная задача не решена, поэтому исследования в данном направлении актуальны. Известны следующие основные способы решения проблемы моделирования капель со сложными распределениями ^[2]:

1. Дискретные модели ^[3], к которым относятся варианты лагранжева подхода с секционированием, а также метод крупных частиц (упоминается, например, в работе ^[4]). Это достаточно простой подход, но для получения адекватной точности требуется моделирование большого количества капель, что дает очень высокую трудоемкость расчета.

2. Эйлерово-эйлеровы модели (модели сплошных сред). Выделим пять основных подходов:

а) применение фиксированных распределений капель (дает недостаточную точность ^[3]);

б) подход А.Е.Алояна ^[5] с решением кинетического уравнения (имеет высокую трудоемкость);

в) атмосферные модели (весьма приближенные, имеющие много допущений о характерных распределениях капель, см., например, ^[6]);

г) метод Хилла (метод моментов ^{[7, 8]}), в котором рассматриваются лишь интегральные характеристики распределений (сами распределения явно не рассчитываются). При всей привлекательности данного подхода имеется ряд серьезных допущений, снижающих его ценность, в частности, капли не имеют собственной скорости ^[9];

д) композитные методы, использующие элементы метода моментов для расчета локальных распределений в каждой расчетной ячейке. Здесь можно выделить метод аппроксимации распределений суммой δ-функций ^{[8, 10]} и предложенный автором ранее метод на базе поиска кусочно-линейных распределений ^{[1, 11]}. Эти методы отличаются сравнительно невысокой трудоемкостью и позволяют аппроксимировать произвольные, в том числе разрывные, распределения. Метод ^[11] позволяют ввести собственные скорости капель, однако в случаях очень активной конденсации или испарения может иметь недостаточную точность, что объясняется, преимущественно, погрешностями расчета межкомпонентных переходов в капельной фазе.

Таким образом, целесообразно рассмотреть в качестве базового подход, изложенный в ^[11] и поставить целью повышение его точности при сохранении сравнительно невысокой трудоемкости. Для достижения данной цели поставим задачи: а) перехода от кусочно-линейных распределений к кусочно-полиномиальным (кусочно-сплайновым); б) проверки адекватности полученной модели в численном эксперименте; в) сравнения точности подхода ^[11] и его предлагаемой модификации.

Расчет предварительных кусочно-линейных распределений капель

Используем следующий подход ^{[1, 11]}:

1. Секционируем капельную среду по диапазонам диаметров D. Начнем с вычисления предварительных кусочно-линейных распределений вида

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

2. Распределения в каждой ячейке расчетной сетки для каждой секции будем восстанавливать (интерполировать) по интегральным плотностям и концентрациям капель i-й секции на каждой итерации. По двум характеристикам можно определить:

а) параметры линейного распределения (при полном заполнении секции, то есть для капель с диаметрами от d_i до d_i+1);

б) параметры равномерного распределения (при этом a_i = 0) и координату начала x_i или конца y_i незаполненного участка (при частичном заполнении секции). На краях распределений знание положения начала/конца заполненного участка важнее, чем построение точного приближения верхней границы.

Такой предварительный расчет простых приблизительных распределений обладает примечательным свойством – он однозначен. Попытка прямого использования более сложных распределений (по большему количеству физических характеристик) приводит к появлению паразитных «альтернативных» решений.

Возможен расчет распределения с разрывами: тогда некоторые секции будут пустыми:

.

Определим функцию mode(i), определяющую тип процесса (конденсация (К), испарение (И) или стабилизация (С)) для каждой i-й секции:

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

определяет знак градиента концентрации, причем C_П — концентрация пара в среде, — концентрация пара на поверхности капли диаметром d в i-й секции.

Пустой считаем секцию, для которой выполняется условие

,

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

Частично заполненной будем считать i-ю секцию, для которой выполняется условие

все прочие непустые секции считаем полностью заполненными.

Логические функции forw(i) и back(i) выделяют два особых случая (промежуточные секции, содержащие «пограничные» капли):

;

.

Прочие основные соотношения приведены в работе ^[11].

Рассчитаем предварительные кусочно-линейные распределения в полностью заполненной секции [x_i; y_i], соответственно

;

.

Запишем выражения для плотности и концентрации секции:

Выполнив несложные преобразования, получим систему линейных уравнений

решением которой будут искомые значения a_i и b_i.

Приведем формулы для нахождения координаты начала x_i или конца y_i незаполненного участка (при частичном заполнении секции). Для нахождения x_i выпишем выражения для плотности и концентрации:

из которых после преобразований получим нелинейное уравнение

,

которое, как легко показать, имеет единственное решение в вещественных числах. Это уравнение решается, если в i-й секции преобладает конденсация. При этом

.

Аналогично выводится уравнение для нахождения x_i:

,

которое также имеет единственное решение в вещественных числах. Это уравнение решается, если в i-й секции преобладает испарение. В данном случае

.

Параметры распределения для частично заполненной секции:

;

.

Расчет уточненных кусочно-сплайновых распределений

Предлагается следующий подход:

1. Определим для каждой пары (i)-(i+1) соседних секций в каждой ячейке уточненные кусочно-сплайновые распределения:

,

которые справедливы при

c_i ≤ D ≤ c_i+1,

d_i ≤ c_i ≤ d_i+1,

где c_i – некоторая (обычно средняя или начальная/конечная) координата в локальном распределении, определяемая формулой

При этом возможны дополнительные ограничения на величину c_i, учитываемые при расчете межкомпонентных (межсекционных) переходов. В частности, c_i должна быть больше максимального диаметра капель, переходящих из i-й в (i-1)-ю секцию, и меньше минимального диаметра капель, переходящих из (i+1)-й в i-ю секцию.

2. Распределения в каждой ячейке расчетной сетки для каждой пары секций будем восстанавливать (интерполировать) по двум точкам n_i(c_i) и n_i+1(c_i+1), а также по сохранению интегральных плотности (массы) и концентрации (количества) при переходе от кусочно-линейных к кусочно-сплайновым распределениям. Этих четырех характеристик достаточно, чтобы определить кубические сплайны путем решения линейной системы из четырех уравнений относительно неизвестных коэффициентов m_i, p_i, q_i, r_i:

где интегралы могут быть легко вычислены аналитически.

Апробация

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

Были проведены три серии экспериментов, различающихся начальным распределением капель по диаметрам (нормальное Г(d), перевернутое нормальное , линейное Л(d)). В каждой серии варьировались:

а) значение относительной влажности воздуха в диапазоне 0,8÷1,2, что позволило рассмотреть как режимы конденсации, так и испарения (с поправкой Кехлера);

б) количество секций — от двух до десяти. Все распределения капель были заданы в диапазоне диаметров [10^‑6; 10^‑5], применялось равномерное разбиение на секции.

Эксперименты проводились при температуре 283 К, на интервале времени 0÷0,01, с шагом интегрирования по времени t = 0,00001. При расчетах учитывалось наличие загрязнителя SO₂, фоновая концентрация которого равнялась 0,88×10^‑2 моль/м³ (200 млн^‑1). При прямом моделировании рассчитывалось множество из 1000 капель.

В каждом j-ом эксперименте вычислялась средневзвешенная (по концентрациям секций) накопленная относительная погрешность результатов (отдельно по плотностям капельных секций, отдельно по концентрациям поглощенного SO₂ в секциях), полученных при использовании нашего подхода, от результатов с прямым моделированием капель:

,

где eps_j,i – относительная погрешность по i-й секции в j-ом эксперименте. Полученные погрешности eps_j усреднялись по схожим экспериментам (всего проводится K таких экспериментов), в которых варьировалась только влажность.

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

по плотности (отложенной по оси ординат, приведены результаты для исходного подхода ^{[1, 11]} и предложенного в данной работе уточненного подхода) в зависимости от числа секций (отложенного по оси абсцисс). График (а) соответствуют нормальному распределению; (б) — перевернутому нормальному; (в) — линейному. На рис. 2 показаны аналогичные графики, но для погрешности расчета концентраций поглощенного каплями SO₂.

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

Рис. 2. Средние относительные погрешности по концентрациям поглощенного SO₂ (для различных начальных распределений) в зависимости от числа секций. Сплошная линия – предлагаемый подход, пунктирная линия – исходный подход

Из графиков очевидно, что в случае линейных исходных распределений результаты расчетов по исходному и предложенному подходам очень похожи, различие составляет 1÷2% (исходный подход чуть более точен по плотности и несколько менее точен по поглощенному SO₂). Это закономерно, поскольку в таком случае, видимо, кусочно-сплайновые распределения вырождаются в близкие к кусочно-линейным.

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

Любопытно заметить, что

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

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

Дополнительно были проведены эксперименты по дальнейшему увеличению количества секций (до 20 и до 40), показавшие падение погрешности по плотности приблизительно в 2÷2,4 раза (до 3÷4%) и по загрязнителю в 4÷5 раз (до 3÷4%) по сравнению с 10-секционным случаем. Это позволяет утверждать, что с ростом числа секций погрешность предложенного подхода все же имеет тенденцию к падению.

Проведенные замеры времени моделирования показали, что расчет с применением предложенного нами подхода требует в 25 раз меньше времени, чем аналогичный эксперимент с прямым моделированием 1000 капель. При этом данный подход требует времени лишь на 10% больше, чем исходный подход ^{[1, 11]}, что объясняется дополнительными временными затратами на расчет сплайнов.

Выводы

Предложен новый, уточненный подход к расчету распределений капельных фаз, отличающийся: а) более точным (улучшение до 30% в случае сложных нелинейных распределений), по сравнению с исходным подходом, расчетом межкомпонентных переходов за счет улучшенной кусочно-сплайновой аппроксимации распределений (средняя относительная погрешность на 10 секциях имеет величину до 15÷16% по плотности и до 27% по загрязнителю, на 20÷40 секциях -- до 3÷4% как по плотности, так и по загрязнителю), б) низкими вычислительными трудозатратами (в 25 раз меньше по сравнению с лагранжевой моделью на 1000 капель), в) возможностью ввода большого количества секций (с ростом точности при увеличении числа секций), что снижает погрешность моделирования динамики капельных фаз, некоторые компоненты которых имеют собственные скорости.