Журнал вычислительной математики и математической физики, 2021, T. 61, № 4, стр. 644-657

1. ВВЕДЕНИЕ

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

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

В последние годы исследователи начали использовать методы вычислительной гидродинамики для эффективного моделирования столкновений пузырьков и частиц во флотационной машине [6], [7]. Несмотря на то что эти модели имеют прикладное значение, многие фундаментальные физические явления, такие как движение пузырьков и частиц, возникающие межфазные силы и процесс дренирования пленки, полностью не изучены. Однако учет межфазных сил является важным фактором при взаимодействии между частицей и пузырьком. В работе [8] теоретически исследована роль гидродинамических и поверхностных сил при взаимодействии пузырька и частицы, и обнаружено существование гидрофобных сил притяжения между поверхностностями, покрытыми поверхностно-активными веществами. В [14] на основе теоретической модели рассчитана гидродинамическая силa, действующая со стороны жидкости на частицу. С изобретением атомно-силового микроскопа многие ученые начали экспериментально изучать межмолекулярные силы взаимодействия, возникающие между пузырьком и частицей [10], [11]. В [12] представлен численно-экспериментальный анализ сил взаимодействия при столкновении всплывающего пузырька с падающей твердой частицей в дистиллированной воде и растворе глицерина. В работах [13], [14] экспериментально изучали процесс прилипания частицы к поверхности большого воздушного пузырька. В [15] проведены аналогичные эксперименты по изучению траектории пузырьков на поверхности большой неподвижной частицы.

Для более полного понимания механизмов взаимодействия пузырьков и частиц при флотационном процессе необходимо рассматривать не только динамику нескольких микрообъектов, но и их кластера. Однако в большинстве теорий, описывающих совместную динамику пузырьков и частиц, трехмерными эффектами пренебрегается, (например, [16]), не учитывается деформация пузырьков, тем самым поток движения жидкости вокруг пузырька и частицы сильно идеализирован. В связи с этим актуальны создание математических моделей и реализация соответствующих программных модулей на базе эффективных методов и алгоритмов в трехмерном случае, которые более адекватно описывают рассматриваемые физические процессы. Основным масштабируемым алгоритмом, используемым в работе, является быстрый метод мультиполей (БММ), который признан одним из 10 лучших алгоритмов 20-го века [17]. Этот метод был разработан в конце 1980-х годов для решения задачи гравитации N-тел и имеет сложность O(N) или O(N logN) [18]. БММ стал одним из самых мощных алгоритмов для решения уравнений Лапласа, Гельмгольца, Стокса, Максвелла и других классических уравнений математической физики. При решении краевых задач он может быть эффективно использован с методом граничных элементов (МГЭ) [19], [20], преимущество которого заключается в возможности отслеживать динамику межфазных границ путем дискретизации поверхности раздела фаз треугольными элементами, что является особенно ценным при моделировании деформируемых пузырьков. БММ допускает возможность получения дополнительного ускорения за счет распараллеливания на графических и центральных процессорах [21]. Наиболее современные работы по применению МГЭ для изучения взаимодействия твердой частицы с пузырьком представлены в [22], [23]. Однако в рассмотренных работах использовалась двумерная постановка и рассматривались большие пузырьки, возникающие при подводном взрыве.

В [24] успешно использовалась комбинация МГЭ и БММ для изучения динамики кластера пузырьков, содержащего более 2 млн расчетных узлов. В отличие от [24] гранично-интегральная формулировка для пузырькового кластера с примесями твердых частиц является более сложной [25]. Однако в ее основе по-прежнему лежат формулы Грина для оператора Лапласа, что позволяет переписать гранично-интегральные уравнения в терминах потенциала простого и двойного слоя. Далее будем называть ядра интегрального представления решений уравнения Лапласа через потенциалы простого и двойного слоя ядрами Лапласа. Таким образом, решение матрично-векторных уравнений для рассматриваемой задачи можно свести к вычислению ядер Лапласа, для которых наиболее эффективно применение БММ.

В работе [25] представлена модификация гранично-интегральной постановки для случая совместной динамики пузырьков и твердых сферических частиц, а также реализовано одно из решений по ускорению кода с применением графических процессоров. Однако распараллеливание кода только на графических процессорах не позволило достичь желаемого результата из-за квадратичной сложности алгоритма. Более того, при реализации кода в [25] не учитывалось представление через ядра Лапласа, а отдельно были построены схемы расчета для левой и правой части системы линейных алгебраических уравнений (СЛАУ). В связи с этим в данной работе для ускорения расчетов и увеличения масштаба задачи совместной динамики пузырьков и частиц в акустическом поле предложен новый эффективный подход, основанный на представлении левой и правой части СЛАУ через потенциалы двойного и простого слоя. Преимущество такого подхода заключается также в возможности применения БММ, который существенно уменьшает сложность вычисления матрично-векторных произведений (МВП) с матрицами специального вида.

2. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим движение пузырька (обозначим индексом “b”) объема V_b, ограниченного поверхностью S_b, и твердой частицы (обозначим индексом “p”) объема V_p, ограниченного поверхностью S_p, в неограниченной идеальной несжимаемой жидкости. Движение жидкости в гравитационном поле ускорения g описывается уравнениями Эйлера

где p, v, ρ_l – давление в жидкости, скорость и плотность соответственно. Далее предполагается, что жидкость покоится на бесконечности, давление на бесконечности определяется согласно действующему акустическому полю P(t) = p₀+ P_asin(ωt), где p₀, P_a, ω – начальное давление в жидкости, амплитуда и циклическая частота акустического поля соответственно.

Обозначим через r_b и r_p радиус-векторы точек на поверхностях S_b и S_p, а через n_b и n_p – нормали к этим точкам (направленные в жидкость) соответственно. Тем самым можно записать кинематические условия на поверхности пузырька и частицы

Будем считать, что газ внутри пузырька ведет себя политропно с показателем политропы κ, тогда давление на границе пузырька вычисляется по формуле

где γ – коэффициент поверхностного натяжения, k – средняя кривизна поверхности, p_g – давление в газе (индексом “0” обозначено начальное состояние).

Динамика твердого тела может быть описана следующим образом. Обозначим координаты и скорость центра масс частицы как R_p и V_p соответственно. Тогда поступательную скорость сферической частицы запишем в виде

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

Кроме того, предполагается, что течение потенциальное, т.е. v = ∇ϕ, где ϕ – потенциал скорости, который удовлетворяет уравнению Лапласа ∇²ϕ = 0. Тогда уравнение (1) для жидкости, покоящейся на бесконечности, можно переписать в форме интеграла Коши–Лагранжа, а потенциал скорости на границе пузырька и частицы определяется согласно динамическому условию

Тогда кинематические граничные условия (2) в случае потенциального течения запишутся в следующем виде:

В отсутствиe частицы задача сводится к интегрированию формулы, которая позволяет определить давление на границе пузырька, если известны текущие позиции точек на поверхности пузырька. Если в некоторый момент времени известен потенциал на поверхности пузырька, то можно решить задачу Дирихле для уравнения Лапласа и определить нормальную производную потенциала скорости согласно условию (6). На основе чего можно найти положение поверхности пузырька на следующем шаге по времени. Потенциал скорости на следующем временном шаге можно вычислить согласно интегралу Коши–Лагранжа (5). Из первого условия (2) можно найти скорость движения поверхности пузырька. Однако разработанная схема не работает для твердого тела, так как давление на поверхности частицы не определяется только ее положением, как в случае с пузырьком. Интеграл Коши–Лагранжа может быть использован, чтобы определить давление на поверхности, но тогда сила, действующая на твердое сферическое тело, должна быть вычислена с помощью формулы (4). Таким образом, сила, действующая на частицу, зависит от нормальной производной потенциала скорости, а это означает, что потенциал на поверхности частицы не может быть найден, как в случае с пузырьком. Поэтому на каждом временном шагe необходимо решать задачу для нахождения потенциала скорости частицы.

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

Подставив в формулу (4) давление, выраженное через потенциал в соответствии с интегралом Коши–Лагранжа (5), получим

Используя (7), уравнение движения частицы можно записать в гранично-интегральной форме

Введя новую переменную

гранично-интегральное уравнение (8) можно переписать в виде

Соотношение (6) с учетом (2) и (9) может быть приведено к виду

Уравнение (11) позволяет исключить производную потенциала скорости по нормали к поверхности частицы и приводит к следующему граничному интегральному уравнению:

где $L[q]$ и $M[\phi ]$ – потенциалы простого и двойного слоя соответственно где $G\left( {{\mathbf{y}},{\mathbf{x}}} \right)$ и $\frac{{\partial G\left( {{\mathbf{y}},{\mathbf{x}}} \right)}}{{\partial n\left( {\mathbf{x}} \right)}}$ – фундаментальное решение уравнения Лапласа и его нормальная производная:

Если потенциал и его нормальная производная известны в момент времени t на поверхности S(t) = S_b(t) ∪ S_p(t), то новые координаты точек и потенциал скорости на поверхности пузырька могут быть вычислены на следующем временном шагe. Неизвестные значения нормальной производной на поверхности пузырька ${{q}_{b}}$ и значения потенциала на поверхности частицы ${{\phi }_{p}}$ могут быть найдены из интегрального уравнения (12), если известны значения потенциала на поверхности пузырька ${{\phi }_{b}}$ в соответствующих точках

Таким образом, задача сводится к решению системы дифференциальных уравнений (СДУ) (5) и (6) для пузырька, (3) и (10) для частицы, которые на первых временных шагах решаются методом Рунге–Кутты 4-го порядка, а затем методом Адамса–Башфорта 4-го порядка. Правая часть СДУ находится из решения гранично-интегрального уравнения (16).

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

3. ДИСКРЕТНЫЙ АНАЛОГ ГРАНИЧНО-ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ

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

где ${\mathbf{y}} = {{{\mathbf{x}}}_{j}}$, ${{N}_{b}},\;{{N}_{p}}$ – количество расчетных узлов на поверхности пузырьков и частиц соответственно, $N = {{N}_{b}} + {{N}_{p}}$ – общее количество расчетных узлов на поверхности объектов, ${{I}_{{ij}}}$ – тензор второго ранга, где верхний индекс соответствует оператору из (16), $j = 1,...,{{N}_{b}}$ относятся к узлам на поверхности пузырьков, $j = {{N}_{b}} + 1,...,N$ – к узлам на поверхности частиц.

Суть нового подхода заключается в представлении уравнения (16) через ядра уравнения Лапласа (14) и использовании ранее реализованных модулей для их расчета, которые показали свою эффективность для моделирования динамики структурированного пузырькового кластера [24].

Поскольку для расчета операторов ${{Q}_{p}}$ и ${{H}_{p}}$ в (13) используется ядро Лапласа $G({\mathbf{y}},{\mathbf{x}})$ (монополь), то представляется возможным переписать операторы в терминах оператора ${{L}_{p}}$

где Тогда тензоры второго ранга в (17), соответствующие операторам ${{Q}_{p}}$ и ${{H}_{p}}$, можно записать через диадик простого слоя где ${{n}_{{ik}}}$ есть k-компонента вектора нормали i-го узла, ${{S}_{i}}$ – площадь, относящаяся к i-му узлу, ${{\widetilde U}_{{ik}}}$ есть k-компонента вектора, отвечающего за скорость движения i-го узла, которая определяется согласно формуле (9). Стоит отметить, что скорость движения каждого узла на поверхности частицы совпадает со скоростью движения центра масс этой частицы.

Диагональные (сингулярные) элементы, соответствующие ядрам $G$ и $\partial G{\text{/}}\partial n$, вычисляются методом вычитания сингулярности на основе интегральных тождеств. Рассмотрим гармоническую и регулярную функцию $\Phi $ в замкнутой области $V$. Для такой функции справедлива формула Грина

В качестве такой функции можно взять $\Phi \left( {\mathbf{y}} \right) \equiv 1$. В этом случае $Q = 0$ и интегральное уравнение (20) перепишется в виде

Таким образом, имеет место следующее тождество для определения сингулярных элементов ядра $\partial G{\text{/}}\partial n$ методом вычитания сингулярности

Используя тождество (22), сингулярную и регулярную часть тензора $I_{{ij}}^{{(\partial G/\partial n)}}$, можно вычислить следующим образом:

Для определения сингулярных элементов ядра $G$ рассматриваются следующие тестовые функции ${{\Phi }_{1}}\left( {\mathbf{y}} \right) \equiv x,$ ${{\Phi }_{2}}\left( {\mathbf{y}} \right) \equiv y,$ ${{\Phi }_{3}}\left( {\mathbf{y}} \right) \equiv z$, которые также являются гармоническими и регулярными в замкнутой области $V$, где $\left( {x,y,z} \right)$ – декартовы координаты точки ${\mathbf{y}}$. Тогда нормальная производная данных функций по поверхности $S$ будет следующей ${{Q}_{1}}\left( {\mathbf{x}} \right) \equiv {{n}_{x}}\left( {\mathbf{x}} \right),$ ${{Q}_{2}}\left( {\mathbf{x}} \right) \equiv {{n}_{y}}\left( {\mathbf{x}} \right),$ ${{Q}_{3}}\left( {\mathbf{x}} \right) \equiv {{n}_{z}}\left( {\mathbf{x}} \right)$ (т.е. компоненты нормали ${\mathbf{n = }}\left( {{{n}_{x}},{{n}_{y}},{{n}_{z}}} \right)$ в точке ${\mathbf{x}}$). Таким образом, используя формулу Грина, получаем тождества

Приведя тождества (24) к дискретному виду, получим переопределенную систему для расчета сингулярных элементов тензора $I_{{ij}}^{{(G)}}$: где ${{F}_{j}}$ и ${{C}_{j}}$ определяются в виде в (26) $\widetilde I_{{ij}}^{{(G)}}$ – элементы тензора $I_{{ij}}^{{(G)}}$ с нулевыми диагональными элементами. Несингулярные элементы тензора $I_{{ij}}^{{(G)}}$ вычисляются как

Для расчета геометрических характеристик поверхности (нормали, площади, кривизна) используются алгоритмы и модули, описанные в [24]–[26].

4. ПРИМЕНЕНИЕ БЫСТРОГО МЕТОДА МУЛЬТИПОЛЕЙ

Тензорная формулировка (17) и (19) сводится к СЛАУ относительно неизвестных значений потенциала скорости ${{\phi }_{p}}$ на поверхности частиц и нормальной производной потенциала скорости ${{q}_{b}}$ на поверхности пузырьков:

где A – расчетная матрица системы, X – вектор неизвестных и b – правая часть системы. Таким образом, при моделировании динамики K_b пузырьков и K_p частиц с общим числом узлов сетки $N = ({{K}_{b}} + {{K}_{p}}) \cdot \Delta N$ требуется решать СЛАУ размером N × N, где $\Delta N$ – количество узлов дискретизации одного дисперсного включения.

Очевидно, что с увеличением количества объектов в кластере прямые методы решения СЛАУ с кубической сложностью алгоритма не справляются. В работе [25] было предложено не хранить расчетные матрицы, а решать полученную СЛАУ с помощью итерационного метода минимальных невязок (GMRES) с расчетом МВП на графических процессорах с использованием технологии CUDA. Однако в данной реализации не использовалась декомпозиция всех МВП через ядра Лапласа, а отдельно рассчитывались левая и правая части (28). Более того, реализованный алгоритм расчета МВП имел квадратичную сложность, что делало его неэффективным для моделирования систем с количеством расчетных узлов более 38 000 [23]. В настоящей работе предложенная декомпозиция (17) и (19) позволяет ускорить расчеты и увеличить масштаб задачи за счет применения многоуровневого БММ, реализованного на гетерогенных вычислительных системах [21].

С учетом представления (17) и (19) для формирования левой части СЛАУ (28) необходимо четыре раза вызвать ядро Лапласа для потенциала простого слоя и один раз – для потенциала двойного слоя. Для расчета правой части СЛАУ (28) требуется три вызова ядра Лапласа для потенциала простого слоя и один вызов – для потенциала двойного слоя.

Особенность БММ состоит в декомпозиции всех матриц, формирующих СЛАУ (28) на разреженную и плотные матрицы:

В (29) разреженная матрица ${{A}^{{{\text{sparse}}}}}$ учитывает ближнее взаимодействие источников ${{{\mathbf{x}}}_{i}}$ и приемников ${{{\mathbf{y}}}_{j}}$, удовлетворяющих условию ${\text{|}}{{{\mathbf{y}}}_{j}} - {{{\mathbf{x}}}_{i}}{\kern 1pt} {\text{|}} < {{r}_{c}}$, где ${{r}_{c}}$ соответствует радиусу, в котором оценивается ближнее взаимодействие. Значение ${{r}_{c}}$ выбирается из необходимой точности и скорости расчета, порядок ${{r}_{c}}$ определяется расстоянием между расчетными узлами. Плотная матрица ${{A}^{{{\text{dense}}}}}$ учитывает дальнее взаимодействие ${\text{|}}{{{\mathbf{y}}}_{j}} - {{{\mathbf{x}}}_{i}}{\kern 1pt} {\text{|}} \geqslant {{r}_{c}}$. Отметим, что в МГЭ источники ${{{\mathbf{x}}}_{i}}$ и приемники ${{{\mathbf{y}}}_{j}}$ совпадают и являются узлами дискретизации. МВП ${{A}^{{{\text{sparse}}}}}{\mathbf{X}}$ вычисляется напрямую по формулам, представленным выше, в то время как МВП ${{A}^{{{\text{dense}}}}}{\mathbf{X}}$ вычисляется приближенно согласно мультипольному разложению в области с центром в точке ${{{\mathbf{x}}}_{*}}$ [28]

где ${{S}_{l}}$ – сингулярные базисные функции, ${{R}_{l}}$ – регулярные базисные функции, которые выбираются из соображений удобства и скорости сходимости рядов [28]. В (30) рассматриваются только первые $P_{{{\text{БММ}}}}^{2}$ членов разложения, где $P_{{{\text{БММ}}}}^{{}}$ – параметр усечения, который определяет точность разложения $\varepsilon = \varepsilon (P_{{{\text{БММ}}}}^{2})$ и самого БММ.

При реализации многоуровневого БММ, кроме построения $S$ и $R$ разложений, используются также операторы трансляций. Трансляция разложений применяется для перехода от представления функции в одной области $\Omega $ с центром ${\mathbf{x}}{\kern 1pt} *$ к представлению функции в другой области $\Omega {\kern 1pt} '$ c центром ${\mathbf{x}}{\kern 1pt} *{\kern 1pt} ' = {\mathbf{x}}{\kern 1pt} * + \;{\mathbf{t}}$, где ${\mathbf{t}}$ – вектор трансляции. Применяя оператор трансляции $T({\mathbf{t}})$ к коэффициентам известного разложения в области $\Omega $ можно получить коэффициенты искомого разложения в области $\Omega {\kern 1pt} '$. В БММ различают три типа операторов трансляции $R\,{\text{|}}\,R$, $S\,{\text{|}}\,S$ и $S\,{\text{|}}\,R$, где первая буква означает тип начального разложения, а вторая – тип конечного разложения.

Для применения многоуровневого БММ необходимо сформировать структуру данных на основе иерархического алгоритма разбиения трехмерного пространства с помощью восьмеричного дерева. Начальный куб уровня 0 разбивается на 8 дочерних кубов, которым присваивается уровень 1. Процесс восьмеричного деления продолжается до некоторого максимального уровня, определяемого из соображений оптимизации. Детальное описание БММ можно найти в публикации [28].

МВП для ближнего взаимодействия вычисляется напрямую на GPU c использованием технологии CUDA, а МВП для дальнего взаимодействия вычисляется через оценку разложений и операторов трансляций на CPU с использованием технологии OpenMP. Точность и время расчета плотной матрицы зависят от количества элементов в разложении P_БММ, чем больше значение этого параметра, тем больше времени требуется на вычисление этой части, но выше точность расчета дальнего взаимодействия [24], [26]. Реализация МВП на GPU позволяет снижать глубину иерархического дерева структуры данных, что обеспечивает дополнительное ускорение алгоритма [21].

5. РЕЗУЛЬТАТЫ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ

Все программные модули, кроме БММ, разработаны в среде Matlab. Расчеты проводились на персональном компьютере, оснащенном CPU Intel Xeon 5660, 2.8GHz, 12 GB RAM с 12 физическими и 12 виртуальными ядрами и GPU NVIDIA Tesla K20 (Kepler), 5GB глобальной памяти.

Расчеты выполнены для безразмерных параметров с характерными величинами: $L{\kern 1pt} * = {{r}_{{b0}}}$, $P{\kern 1pt} * = {{p}_{0}}$, $\rho {\kern 1pt} * = {{\rho }_{l}}$, $U{\kern 1pt} * = \sqrt {{{p}_{0}}{\text{/}}{{\rho }_{l}}} $, $T{\kern 1pt} * = L{\kern 1pt} *{\text{/}}U{\kern 1pt} * = {{r}_{{b0}}}{\text{/}}\sqrt {{{p}_{0}}{\text{/}}{{\rho }_{l}}} $. Согласно ограничениям математической модели, описывающей физические процессы для больших чисел Рейнольдса и малых капиллярных чисел, были выбраны следующие безразмерные параметры: $r_{{b0}}^{'} = 1$, $r_{p}^{'} = 1$, $\rho _{l}^{'} = 1$, $\rho _{p}^{'} = 2$, $p_{0}^{'} = 1$, $P_{a}^{'} = 1$, $\omega {\kern 1pt} '{\text{/}}(2\pi ) = 0.2$, $\gamma {\kern 1pt} ' = 0.073$, $\kappa {\kern 1pt} ' = 1.4$, $g{\kern 1pt} ' = 0$, обозначения которых соответствуют размерным параметрам, описанным выше.

Для тестирования производительности реализованного подхода рассматривался структурированный кластер пузырьков и частиц, сформированный из n × n × n кубов, каждый из которых содержит четыре пузырька и четыре частицы. Центры дисперсных включений расположены на расстоянии четырех радиусов. Пример такого куба, который соответствует кубическому кластеру с параметром n = 1, представлен на фиг. 1а, где серым цветом раскрашены пузырьки, а черным – частицы. На фиг. 1б представлена форма кластера, состоящая из 27 кубов (n = 3). Количество пузырьков и частиц в кластере размера n можно рассчитать по формуле ${{K}_{b}} = {{K}_{p}} = 4{{n}^{3}}$.

Для оценки производительности кода в зависимости от точности расчетов рассматривалось 6 случаев.

1. БММ/CPU: расчет всех МВП БММ с вычислением ближнего взаимодействия на CPU, дальнего взаимодействия с точностью $\;{{P}_{{{\text{БММ}}}}} = 12$, решение СЛАУ с точностью$\;{{\varepsilon }_{{{\text{GMRES}}}}} = {{10}^{{ - 5}}}$.

2. БММ/CPU: расчет всех МВП БММ с вычислением ближнего взаимодействия на CPU, дальнего взаимодействия с точностью $\;{{P}_{{{\text{БММ}}}}} = 8$, решение СЛАУ с точностью$\;{{\varepsilon }_{{{\text{GMRES}}}}} = {{10}^{{ - 4}}}$.

3. БММ/GPU(double precision): расчет всех МВП БММ с вычислением ближнего взаимодействия на GPU c двойной точностью, дальнего взаимодействия с точностью $\;{{P}_{{{\text{БММ}}}}} = 12$, решение СЛАУ с точностью$\;{{\varepsilon }_{{{\text{GMRES}}}}} = {{10}^{{ - 5}}}$.

4. БММ/GPU(single precision): расчет всех МВП БММ с вычислением ближнего взаимодействия на GPU c одинарной точностью, дальнего взаимодействия с точностью$\;{{P}_{{{\text{БММ}}}}} = 8$, решение СЛАУ с точностью$\;{{\varepsilon }_{{{\text{GMRES}}}}} = {{10}^{{ - 4}}}$.

5. GPU(double): расчет всех МВП на GPU c двойной точностью, решение СЛАУ с точностью$\;{{\varepsilon }_{{{\text{GMRES}}}}} = {{10}^{{ - 5}}}$.

6. GPU(single): расчет всех МВП на GPU c одинарной точностью, решение СЛАУ с точностью$\;{{\varepsilon }_{{{\text{GMRES}}}}} = {{10}^{{ - 4}}}$.

В табл. 1 представлено среднее время расчета левой части, правой части и решения СЛАУ для пузырькового кластера с примесями для структур n = 4, 5, 6, 7. Из таблицы видно, что наибольшее ускорение получено в четвертом случае для БММ с расчетом ближнего взаимодействия на GPU с одинарной точностью. Однако стоит отметить, что ускорение в этом случае по сравнению с третьим достигается за счет потери точности вычислений из-за операций на GPU с одинарной точностью, меньшего количества элементов в разложении по сингулярному базису и точности GMRES. Так, например, для структуры из 1372 пузырьков и 1372 частиц, что соответствует $N = 1\,761\,648$ расчетным узлам, расчет правой части занимает 4 с, расчет левой части – 5 с, GMRES сходится за 9 итераций, на решение СЛАУ требуется 51 с, что в 60 раз быстрее, чем требуется на аналогичный расчет на GPU c одинарной точностью. На фиг. 2 представлено время вычисления левой и правой части СЛАУ в логарифмической системе координат, откуда видно, что сложность модулей, реализованных с использованием БММ, (случай 1–4) линейная, а прямой расчет на GPU (cлучай 5, 6) имеет квадратичную сложность. Заметим, что для $N \leqslant 40\,000$ эффективнее использовать модули 5 и 6.

На фиг. 3 представлены формы пузырькового кластера с примесями твердых сферических частиц, состоящего из четырех пузырьков и четырех частиц (фиг. 1а) в фазе расширения (t = 0.6T, фиг. 3а) и сжатия (t = 0.85T, фиг. 3б). Из фигуры видно, что при сжатии в пузырьках образуются струи, направленные в центр кластера. Проведен анализ деформации углового пузырька для различных размеров кластера (фиг. 4). Анализ показал, что в случаях n = 1, 2 и 3 в крайних пузырьках формируются струи, направленные в центр кластера. Однако с увеличением размера кластера деформация крайних пузырьков уменьшается. Так, например, в случае структуры кластера, содержащего 854 пузырька и 854 частицы, что соответствует n = 6, форма углового пузырька остается сферической. С увеличением размера кластера также уменьшается мобильность пузырьков (фиг. 5а) и частиц (фиг. 5б), что обусловлено незначительными изменениями объема отдельных пузырьков (фиг. 6а) и кластера в целом (фиг. 6б), следовательно, и слабыми гидродинамическими потоками, создаваемыми осциллирующими пузырьками.

6. ЗАКЛЮЧЕНИЕ

Для проведения расчетов высокой вычислительной сложности требуются высокопроизводительные вычисления, которые включают в себя две составляющие: усовершенствованные масштабируемые алгоритмы и вычислительные системы с большим количеством ядер. В данной работе основу таких алгоритмов составил быстрый метод мультиполей (БММ), который позволяет решать СЛАУ с миллионами неизвестных на относительно недорогих персональных суперкомпьютерах, оборудованных графическими процессорами (GPU) и уже показал свою эффективность для прямого моделирования динамики пузырькового кластера [24]. Указанный подход был использован для ускорения метода граничных элементов, позволяющего эффективно рассчитывать динамику больших пузырьковых кластеров с примесями твердых сферических частиц.

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

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

Авторы выражают благодарность Н.А. Гумерову (Институт компьютерных исследований штата Мэриленд, США) за постановку задачи, компании Fantalgo, LLC (Мэриленд, США) за предоставление библиотеки БММ.