CC BY
28
УДК 519.642.2:6
Р. Р. Загидуллин1, А. П. Смирнов2, С. А. Матвеев3, Е. Е. Тыртышников4
ЭФФЕКТИВНЫЙ МЕТОД ЧИСЛЕННОГО РЕШЕНИЯ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ПЕРЕНОСА КОАГУЛИРУЮЩИХ ЧАСТИЦ
В работе рассмотрена математическая модель переноса коагулирующих частиц, для которой предложен новый вычислительный алгоритм, основанный на использовании быстрого метода вычисления интегральных операторов, описывающих процесс коагуляции. Алгоритм позволяет существенно снизить вычислительную сложность реализации шагов явной по времени разностной схемы. Вместо 0(ЫМ2) арифметических операций сложность выполнения каждого шага по времени снижена до О(ЫМКЫМ) операций, где N — число узлов вдоль физических координат частиц, М — число узлов в объемной сетке, соответствующей размерам коагулирующих частиц, ай — ранг матрицы, соответствующей значениям функции ядра коагуляции в узлах сетки. Построенный метод позволяет значительно ускорить расчеты в случае ядер с небольшими рангами Я.
Ключевые слова: малоранговые матричные разложения, коагуляция, уравнение переноса, численные методы.
1. Введение. Процессы коагуляции (слияния) и дробления хаотически движущихся частиц лежат в основе разнообразных физических явлений: от роста мелких полимерных цепочек [1, 2] до образования звезд [3]. В пространственно-однородном случае данные процессы могут быть описаны при помощи феноменологических уравнений Смолуховского, для которых известны различные эффективные численные алгоритмы получения приближенных решений [4-9].
В работе рассматривается математическая модель переноса коагулирующих частиц вдоль отрезка, полученная в результате объединения уравнений коагуляции и переноса:
V
^ ' '—- + c(v) ^ '—-—- = - [ К (и, v — и) f(t, ж, и) f(t, x,v — и) du — at ox 2 J
о
ос
—f(t,x,v) J K(u,v)f(t,x,u) du,
0
где значение функции f(t, x, v) соответствует значению концентрации на единицу объема среды частиц размера v в точке с координатой х в момент времени t. Вдоль оси координаты х ^ 0 осуществляется перенос коагулирующих частиц со скоростью c(v). Первый интегральный оператор из правой части соответствует росту концентрации частиц размера v в результате слияний частиц размеров и и v — и, второй соответствует уменьшению концентрации вследствие слияний частиц размера v с другими частицами среды. Скорости процессов слияния частиц определяются значениями функции К (и, v) — ядра коагуляции. При известном граничном условии f(t, 0, v) и заданных начальных условиях /(Q,x,v) имеем задачу Коши для модели переноса коагулирующих частиц.
Дополнительно укажем максимальный допустимый размер частиц VTпх - В терминах модели это означает, что при достижении размера V^ax частицы перестают участвовать в процессе коагуляции, например вследствие гравитации [10, 11]. Таким образом, получаем задачу Коши в области
1 Факультет ВМК МГУ, студ., e-mail: zagidullinrishatQgmail.com
2 Факультет ВМК МГУ, доц.; ИВМ РАН, науч. сотр., к.ф.-м.н., e-mail: sapQcs.msu.su
3 Сколковский институт науки и технологий, мл. науч. сотр.; факультет ВМК МГУ, асп.; ИВМ РАН, мл. науч. сотр., e-mail: s.matveevQskoltech.ru
4 ИВМ РАН, директор; Московский физико-технический институт, проф.; факультет ВМК МГУ, проф., e-mail: eugene.tyrtyshnikovQgmail.com
* Работа выполнена при поддержке Российского научного фонда (проект № 14-11-00806).
(1,х,у) € [0; сю] X [0; сю] х [0;Утах]:
<Э/(£, ж, г>) <Э/(£, ж, г>) 1
^гг^--И с(у) ^ '—-— = - [ К (и, V — и)/(1 ж, и)/(1 ж, V — и) йи —
от ож 2 J
^тах
—/(£, ж, -и) J К (и, ж, и) (1и.
(1)
о
Использование этой математической модели позволяет отказаться от предположения о пространственной однородности процесса коагуляции и расширить класс потенциально решаемых задач.
В п. 2 описывается схема численного решения уравнения переноса коагулирующих частиц. Для численного решения уравнения модели предлагается метод, основанный на явной схеме динамического интегрирования уравнений по времени. Для аппроксимации оператора переноса используется схема, полученная методом конечных объемов с функцией-ограничителем. Для вычисления оператора коагуляции применяются быстрые методы вычисления интегральных операторов [5], входящих в оператор коагуляции Смолуховского. В п. 3 обсуждаются результаты численных экспериментов. В заключении приводятся выводы и предлагаются направления для дальнейшего изучения математических моделей, записанных при помощи уравнений рассматриваемого типа.
2. Численный метод. Для численного решения рассматриваемого уравнения мы приведем схемы вычисления интегральных операторов [4, 5] коагуляции и дифференциальных операторов переноса.
2.1. Метод вычисления интегральных операторов Смолуховского. Опишем схему моделирования процесса коагуляции частиц в рамках уравнения (1). Сначала рассмотрим уравнение коагуляции с предельным размером частиц Ута1Х:
1, V____
<Э/(£, ж, у) 1
т
= ~У К (и, V — «)/(£, ж, «)/(£, ж, V — и) йи — /(¿, ж, у) ^ К(и,у)/(1,х,и)йи.
Для решения этого уравнения можно применить метод трапеций вычисления интегралов и явную схему 1'унго Кутты второго порядка интегрирования по времени [5].
Для ускорения вычисления интегралов построим скелетное разложение ядра коагуляции при помощи крестового алгоритма [12, 13]:
я
К (и, у) = аа(и)Ьа(у).
а=1
Таким образом, для первого интеграла имеем
я
/ К (и, V — «)/(£, ж, «)/(£, ж ,у — и) йи = / аа(и)Ьа(у — «)/(£, ж, «)/(£, ж ,у — и) йи = О О а=1
л "
= / аа(и)/^,Х,и)Ьа(у — и)/(1,Х,У — и) <1и.
а=1 '
Аналогично для второго интеграла
к к
^'тах
К (и, и)ж, и) йи ~ J аа(и)Ьа(у)/(1, ж, и) йи = ^^ Ъа(у) J аа («)/(£, ж, и) йи.
О О а=1 а=1
Заметим, что первый интеграл является суммой из К нижнетреугольных сверток, каждую из которых можно вычислить с использованием алгоритма быстрого преобразования Фурье
за 0(М1пМ) операций. Скелетное разложение ядра коагуляции позволяет вычислить второй интеграл за О(МВ) операций. Таким образом, используя квадратурную формулу трапеций, имеем схему вычисления интегралов Смолуховского, обладающую вторым порядком точности с алгоритмической сложностью О(МКЫМ) операций. Далее в этой работе будем пользоваться быстрой схемой вычисления интегральных операторов из правой части уравнения Смолуховского. Подробности вывода схемы и обоснование ее эффективности приведены в работах [4, 5].
2.2. Аппроксимации уравнения переноса. Рассмотрим уравнение переноса
df(t,x,v) = df(t,x,v) dt дх
и построим аппроксимирующую его схему.
Как известно, при численном решении уравнений переноса с использованием схемы первого порядка точности по времени возникает искусственная диффузия, существенно искажающая расчеты. В то же время в линейных схемах порядка точности выше первого не выполняется сохранение свойства монотонности схем, что приводит к возникновению неестественных колебаний аппроксимируемой функции. В качестве решения проблемы сохранения монотонности и повышения порядка точности был разработан класс схем, представляющих собой комбинацию схем первого и второго порядка точности под названием TVD-схемы (Total Variation Diminishing).
Указанные схемы обладают следующим свойством:
TV(/n+1) < TV (Л, TV(/) = /,-+1 - fj I,
которое обеспечивает монотонность разностной схемы (/" ^ /™+1 —> ^ работе [14]
было показано, что схема для уравнения переноса является Т\'1)-схомой. если она имеет следующий вид:
( .. Г / гп . — Г^
Í fn _ f
Q, I Jj+i -h
3 I jn _ jn
3-1
dt
где /х = с---число Куранта. При этом должно выполняться соотношение
dx
Sd
M)
Ci
f.
n
3+1
C'j-i
^ 1.
fn _ fn /3 ->3-1
Если потребовать выполнения условия Куранта, то получим условие для С-у.
^ о 0 < Сз < 2,
en
'JL
fn
13 + 1
/?+1 - f
fn _ fr,
jj+1 Ji
p
< о -»• Cj = o.
Сеточная функция Cj называется "flux limiter function". Чтобы выяснить практический смысл данной функции, запишем полученную схему следующим образом:
еП+1 ¿j_
dt
с 2 dx
1 + с2
с2 di
сз - Cj-i)fJ - (2 - Cj-i)/"_i]
2 dx2
[Cj,fi+l - (Cj + Cj-i)fj' + ('j ,./'" ,]
Эта схема по структуре напоминает схему Лакса-Вендроффа с коэффициентами, зависящими от сеточной функции Cj . Напомним, что значения Cj определяются необходимостью сохранения свойства монотонности TVD-схемы. Таким образом, полученную схему можно определить как схему Лакса-Вендроффа с ограничителями, сохраняющими ее монотонность.
Вопрос выбора конкретной функции Cj был исследован во многих работах (см., например, [1517]). В результате создано множество функций, удовлетворяющих указанным ограничениям. В данной работе применяется функция, носящая в литературе назавание "Monotonized Central Flux
Limiter" (функция МС) [16]:
Cj = max [О, miliar,, 0.5(1 + гД 2)].
U ~ h-1 fj+1 ~~ fj
В области экстремума функции rj < 0 имеем разностную схему первого порядка. В остальных случаях схема аппроксимирует решение со вторым порядком точности по времени [14].
В качестве условия на правой границе определяется поглощающий слой в виде функции а(х), которая принимает нулевые значения везде кроме небольшого слоя, где переносимый сигнал уничтожается. Подобные слои распространены в работах по исследованию задач электродинамики. Мы реализуем поглощение при помощи PML-слоя (Perfectly Matched Layer). Уравнение переноса заменяется уравнением
с функцией поглощения
= f
at сдх caJ
1)W In 10 fx-I
a = a(x) = < d \ d
x € [I; I + d], иначе,
где I — ширина области расчета без учета поглощающего слоя, д, — ширина поглощающего слоя, Ш ж т — параметры поглощающего слоя.
Выбор функции а(х) определяется требованием сохранения непрерывности искомой функции / и ее производной по пространству. В работе [18] было предложено использовать слои вида
(х^1\т _(т + 1)Шп10
(Т{Х) — Стах I I ; О'тах —
} " ШЙА 7
/ ®
При оптимальных значениях 2 < тп < 3 данный поглощающий слой уменьшает значение переносимой функции на внешнем слое в 10^' раз.
2.3. Полная аппроксимация задачи. Аппроксимируя интегральные операторы в правой части уравнения Смолуховского и уравнения переноса, получаем явную по времени схему для общей модели переноса коагулирующих частиц:
где оператор А(/) соответствует процессу переноса, <5(/) — результат аппроксимации интегралов Смолуховского. Пусть М — число узлов сетки, соответствующей распределению частиц по размерам, а N — число узлов вдоль оси х.
Для каждой фиксированной координаты а.ц воспользуемся быстрым методом вычисления интегральных операторов Смолуховского за О(МКЫМ) арифметических операций, где К — ранг ядра коагуляции. Так как оператор коагуляции и разностный шаблон для оператора переноса необходимо вычислить в каждом узле сетки, соответствующем физическим координатам частицы, то общая сложность одного шага схемы будет составлять 0(Ж + iVMi?lnM) = 0(iVMi?lnM) арифметических операций. Отметим, что без использования быстрого метода вычисления оператора Смолуховского сложность схемы составит О (ММ2) арифметических операций, что сделало бы ее бесполезной для проведения прикладных расчетов. Отметим, что сложность шагов предложенной разностной схемы существенно зависит от значения ранга Д ядра коагуляции. Так, в работах [5, 19] для широкого класса ядер, используемых в приложениях, доказаны оптимистичные оценки на Д в зависимости от точности их приближения. Далее приводятся результаты численных экспериментов для задач Коши с постоянным и баллистическим ядрами коагуляции, для которых гарантированно существуют малоранговые представления с небольшими значениями Д [5, 19].
3. Результаты численных экспериментов. Результаты расчетов представлены на рис. I
4, а также в таблице. Во всех проведенных численных экспериментах использовались нулевые
2
начальные условия и следующее граничное условие: /(£, 0,г>) = . На рис. 1, 2 приведены результаты экспериментов с постоянным ядром коагуляции К (и. г>) = 1 для двух типов скорости
1
переноса частиц: с(г>) = 1 (рис. 1) и с(г>) = —(рис. 2), а на рис. 3. 4 результаты экспериментов с баллистическим ядром
Из полученных результатов следует, что существует такой выбор параметров модели, при котором образуются пространственно неоднородные стационарные распределения частиц по размерам. Сравнивая результаты, представленные на рис. 1,<? и 2,д, а также на рис. 3,<? и 4Д можно утверждать, что выбор скорости переноса существенно влияет на вид решения: получены различные структуры линий уровня при одинаковых начальных условиях. Приведенные на рис. 1,е и 3,в результаты демонстрируют существенное влияние ядра коагуляции на вид решения: наблюдаются различные пространственные распределения концентраций частиц совпадающих размеров при одинаковых начальных условиях.
Рис. 1. Графики и тепловые карты решения: а перенос вдоль оси х частиц размеров V в момент времени £ = 2; б в момент времени £ = 4; в в момент времени £ = 6; я распределение частиц по размерам в момент времени £ = 2; с) в момент времени £ = 4; е в момент времени £ = 6
Юс,У,г)
I
Рис. 2. Графики и тепловые карты решения: а перенос вдоль оси х частиц размеров -о в момент времени £ = 2; б в момент времени £ = 4; в в момент времени £ = 10; я распределение частиц по размерам в момент времени £ = 2; с) в момент времени £ = 4; е в момент времени £ = 10
Время расчетов с использованием предложенной методологии для области у.х € [0,5], I € € [0,10], где Ах = Дг> = А/г = 0.05, = 0.01, приведено в таблице. В случае экспериментов с баллистическим ядром коагуляции ранг полученных скелетных разложений Л = 7. Результа-
4=0.2 -4=2.5 -4=35
4=0.2 4=25 4=3.5
4=0.2 -4=2.5 -4=3.5 -
/ГвдО
Рис. 3. Графики и тепловые карты решения: а перенос вдоль оси х частиц размеров х> в момент времени £ = 2; б в момент времени £ = 4; в в момент времени £ = 6; г распределение частиц по размерам в момент времени I = 2; с) в момент времени I = 4; е в момент времени 1 = 6
4=0.2 -4=2.5 -4=3.5 -
К",
4=0.2 4=25 4=3.5
/Гад
у-0.2 у=2.5 у=3.5
I
Юс,4,г)
Рис. 4. Графики и тепловые карты решения: а перенос вдоль оси х частиц размеров х> в момент времени £ = 2; б в момент времени £ = 4; в в момент времени £ = 10; г распределение частиц по размерам в момент времени I = 2; с) в момент времени I = 4; е в момент времени I = 10
ты численных экспериментов согласуются с приведенной арифметической сложностью алгоритма О(ММТШпМ), где N число узлов в сетке вдоль оси х, М число узлов в сетке размеров г>, Т полное количество шагов по времени, а Л ранг скелетного разложения ядра коагуляции.
Время расчетов экспериментов, с
Ядро коагуляции К(и, V) Скорость переноса Ф>) Параметры сетки
2АН, 2Д£ АН, Д* Д/( лt 2 ' 2 Д/( Лt 4 ' 4
Константное Константное Баллистическое Баллистическое Постоянная Убывающая Постоянная Убывающая 0.737 0.707 2.643 2.522 4.190 4.346 21.176 21.332 35.213 34.718 301.308 270.761 228.024 248.944 1779.738 2087.567
4. Заключение. В данной работе был предложен эффективный численный метод решения уравнения, описывающего перенос коагулирующих частиц. Метод основан на использовании монотонной схемы для аппроксимации оператора переноса и быстрого метода вычисления интегрального оператора Смолуховского [5]. Алгоритмическая сложность одного шага численной схемы составляет 0(ЖЛМ1пМ) операций вместо исходных 0(ЫМ2). Снижение сложности получено путем применения малоранговых матричных аппроксимаций ядра коагуляции и быстрых алгоритмов линейной алгебры.
В результате численных экспериментов продемонстрировано существенное влияние выбора функции скорости переноса коагулирующих частиц c(v) и функции ядра коагуляции К (и, v) на вид численного решения рассматриваемой модели.
В дальнейшем мы планируем провести анализ рангов численных решений рассматриваемой модели и исследовать возможности расширения предложенного быстрого алгоритма для представления решения модели в виде малоранговых разложений. Предложенная численная схема легко обобщается и для моделирования процессов переноса многокомпонентных коагулирующих частиц [7, 19-21].
СПИСОК ЛИТЕРАТУРЫ
1. Poeschel Т., Brilliantov N. V., Frommel С. Kinetics of prion growth // Biophysical J. 2003. 85. P. 3460-3474.
2. Eigen M. Prionics or the kinetic basis of prion diseases // Biophysical Chemistry. 1996. 63. N 1. P. A1-A18.
3. Brilliantov N., Krapi vsky P. L., Bodrova A., et al. Size distribution of particles in Saturn's rings from aggregation and fragmentation // Proceedings of the National Academy of Sciences. 2015. 112. N 31. P. 9536-9541.
4. Матвеев С. А., Тыртышников E.E., Смирнов А.П. и др. Быстрый метод решения уравнений агрегационно-фрагментационной кинетики типа уравнений Смолуховского // Вычислительные методы и программирование. 2014. 15. С. 1-8.
5. Matveev S.A., Smirnov А.P., Tyrtyshnikov Е. Е. A fast numerical method for the Cauchy problem for the Smoluchowski equation //J. Comput. Phys. 2015. 282. P. 23-32.
6. Sorokin A. A., Strizhov V.F., Demin M.N., Smirnov A. P. Monte Carlo modeling of aerosol kinetics // Atomic Energy. 2015. 117. N 4. P. 289.
7. Chaudhury A., Oseledets I., Ramachandran R. A computationally efficient technique for the solution of multi-dimensional PBMs of granulation via tensor decomposition // Computers & Chemical Engineering. 2003. 61. P. 234-244.
8. Debry E., Sportisse В., Jourdain B. A stochastic approach for the numerical simulation of the general dynamics equation for aerosols //J. Comput. Phys. 2003. 184. N 2. P. 649-669.
9. Матвеев С. А. Параллельная реализация быстрого метода решения уравнений агрегационно-фрагментационной кинетики типа уравнений Смолуховского // Вычислительные методы и программирование. 2015. 16. С. 360-368.
10. Fuks N. A. The Mechanics of Aerosols. New York: Dover Publications, 1964.
11. Алоян A.E. Динамика и кинетика газовых примесей и аэрозолей в атмосфере. М.: ИВМ РАН, 2002.
12. Tyrtyshnikov Е. Е. Incomplete cross approximation in the mosaic-skeleton method // Computing. 2000. 64. N 4. P. 367-380.
13. Ж e л т к о в Д. А., Т ы р т ы ш н и к о в Е. Е. Параллельная реализация матричного крестового метода // Вычислительные методы и программирование. 2015. 16. Р. 369-375.
14. Sweby Р. К. High resolution schemes using flux limiters for hyperbolic conservation laws // SIAM J. Numer. Analysis. 1994. 21. N 5. P. 995-1011.
15. К or en B. A robust upwind discretization method for advection, diffusion and source terms // Notes on Numerical Fluid Mechanics. 1993. 45. P. 117-138.
16. Leer B. Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection //J. Comput. Phys. 1977. 23. N 3. P. 276-299.
17. Lyra P.R.M., Morgan K., Peraier J., et al. TVD algorithms for the solution of the compressible Euler equations on unstructured meshes // Intern. J. Numer. Meth. Fluids. 1994. 19. P. 827-847.
18. Berenger J.-P. A perfectly matched layer for the absorption of eletromagnetic waves //J. Comput. Phys. 1994. 114. P. 185-200.
19. Matveev S.A., Zheltkov D. A., Tyrtyshnikov E. E., et al. Tensor train versus Monte Carlo for the mult ¡component Smoluchowski coagulation equation //J. Comput. Phys. 2016. 316. P. 164-179.
20. Matveev S. A., Zheltkov D. A., Tyrtyshnikov E.E., et al. Fast and accurate finite-difference method solving multicomponent Smoluchowski coagulation equation with source and sink terms // Procedia Computer Science. 2016. 80. P. 2141-2146.
21. Palaniswaamy G., Loyalka S. K. Direct simulation Monte Carlo aerosol dynamics: coagulation and collisional sampling // Nuclear Technology. 2006. 156. N 1. P. 29-38.
Поступила в редакцию 07.06.17
