Scientific journal

Fundamental research

ISSN 1812-7339

"Перечень" ВАК

ИФ РИНЦ = 1,674

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

Данная модель описывает перенос вещества в порах диаметра от 1 до 1000 нм, длиной от 0,01 до 200 мкм, при давлениях от 0,001 до 2 атм. При малых размерах пор перенос вещества происходит как за счет столкновения молекул друг с другом (молекулярная диффузия), так и за счет соударений молекул со стенкой поры (диффузия Кнудсена), причем второй тип преобладает, так как вероятность соударения молекул друг с другом при низком давлении довольно мала. Расчет количества молекул в поре производится после задания начальных условий: длина, диаметр поры, температура, давление газа и его молекулярная масса. Средняя скорость молекул газа рассчитывается из распределения Максвелла-Больцмана. После вычисления средней скорости молекул газа, каждой молекуле присваивается определенная скорость с учетом вероятности из распределения Максвелла.

Для описания движения и взаимодействия молекул рассматривается два варианта. В первом варианте молекулы движутся по законам классической механики и взаимодействуют друг с другом, а также со стенкой, по принципу абсолютно-упругого удара. Использован механизм Кнудсена-Смолуховского для описания взаимодействия со стенкой поры. Принимается, что частицы могут сталкиваться со стенкой поры двумя разными способами: зеркальным и диффузным, причем для каждого конкретного соударения способ определяется случайным образом, и соотношение количеств соударений по обоим способам является одним из параметров модели.

Во втором варианте взаимодействие молекул описывается с применением межмолекулярного потенциала взаимодействия Леннарда-Джонса:

 (1)

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

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

 (2)

где D - коэффициент диффузии; - среднее значение скорости частиц в системе; - средняя длина свободного пробега частиц. Эти данные легко извлекаются в процессе моделирования, так как известны параметры движения каждой отдельной частицы.

Расчетные значения коэффициента диффузии, определяемые по соотношению (2), сравниваются с теоретической формулой коэффициента диффузии [1]:

 (3)

где D_Т - теоретический коэффициент диффузии; D_К - коэффициент диффузии Кнудсена, учитывающий только столкновения со стенкой и определяемый соотношением [2]:

 (4)

где r_p - радиус поры; R - универсальная газовая постоянная; T - температура; M - молярная масса вещества (в случае многокомпонентной смеси используется средняя молярная масса смеси); D_М - коэффициент молекулярной диффузии бинарной смеси, учитывающий столкновение молекул с друг другом и описываемый соотношением [3]:

 (5)

где σ_AB - средний поперечный размер частицы, который определяется из соотношения σ_AB = (σ_A + σ_B)/2, где σ_A, σ_B - диаметры молекул; k - постоянная Больцмана, m^* - средняя молекулярная масса смеси, которая определяется из соотношения 1/m^* = 1/m_A +1/m_B, где m_A, m_B - масса молекул; p - давление.

Количество молекул в поре рассчитывается в соответствии с законом Менделеева-Клапейрона:

 (6)

где N_A - число Авогадро; V - объем поры.

Средняя начальная тепловая скорость молекул определяется из соотношения:

 (7)

где m - масса молекулы вещества.

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

  (8)

где v_x, v_y, v_z - проекции скоростей на координатные оси x, y, z.

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

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

Организация параллельных вычислений осуществляется с применением технологии CUDA. Преимущество использования технологии параллельного программирования заключается в том, что при моделировании необходимо проводить множество однотипных операций над всеми частицами (проверка на столкновения, численное решение уравнений движения и т.д.). Для такой задачи хорошо подходят решения на основе принципа SIMD (single instruction-stream, multiple data-stream - один поток команд с несколькими потоками данных), когда вычислительное устройство выполняет операции сразу над несколькими элементами данных за один такт [4]. На этом принципе построена и используемая в данной реализации технология CUDA.

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

На рис. 1 представлены зависимости коэффициента диффузии аргона и метана от изменения температуры, определяемые по формуле (3) и рассчитанные по первому варианту взаимодействия частиц, с привлечением соотношения (2).

Расчеты проведены при диаметре поры 40 нм, длине поры 2 мкм, давлении 0,02 атм. Из рис. 1 видно, что при выполнении условия Кнудсена рассчитанные значения коэффициента диффузии по варианту 1 полностью совпадают с данными по теоретической формуле (3), этот результат свидетельствует об адекватности математической модели в условиях Кнудсена. Следует отметить, что для расчета столкновений молекул со стенкой поры наилучшее совпадение для коэффициента диффузии дал механизм диффузного рассеяния молекул при соударениях со стенками пор.

На рис. 2 представлены зависимости изменения коэффициента диффузии аргона от размера поры, определяемые по соотношению (3) и в результате расчетов вариантов 1 и 2 с применением соотношения (2). Расчеты сделаны при температуре 298 К, давлении 0,02 атм. и длине поры 2 мкм. Из рис. 2 видно, что учет сил межмолекулярного взаимодействия не влияет на определение коэффициента диффузии. Следует отметить, что при малых диаметрах пор данные моделирования соответствуют формуле (3). При больших диаметрах наблюдаются отклонения от формулы, которые вызваны тем, что в системе возрастает доля соударений молекул друг с другом по отношению к общему числу соударений (друг с другом и со стенкой).

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

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

Математическое моделирование на основе модели твердых сфер с организацией параллельных вычислений позволяет определить все параметры модели по длине поры (каналу мембраны), в отличие от определения коэффициента диффузии по соотношению (3). Расчеты проводились для определения коэффициента диффузии водорода для диаметра поры 40 нм, давления 0,02 атм, температуры 298 К, длина поры менялась от 1 до 20 мкм. На малых длинах пор ~2 мкм значения коэффициентов диффузии водорода, рассчитанные на основе модели твердых сфер и по соотношению (3), совпали. С увеличением длины поры наблюдается значительное снижение значения коэффициента диффузии водорода с 2,4∙10^-5 до 7,3∙10^-6 м²/с. Это связано с тем, что при увеличении длины в поре накапливаются частицы с низкой скоростью теплового движения, в результате чего средняя скорость теплового движения внутри поры несколько снижается, и в соответствии с расчетной формулой (2) коэффициент диффузии также падает.

Проводились расчеты по массопереносу бинарных смесей: Ar (95 об. %) - Н₂ (5 об. %) и СО₂ - СН₄ (составов 10:90, 20:80, 40:60, 60:40, 80:20, 90:10 соответственно СО₂, СН₄, об. %), параметры математической модели определялись по длине поры: число столкновений молекул между собой и стенкой, средние скорости молекул, распределение числа молекул, концентрация компонентов смеси, фактор разделения. В ходе расчета проводились итерационные процедуры ~160 млн итераций, с шагом итерации Δτ = 5∙10^-13 секунд до тех пор, пока процесс массопереноса не выходил на стационарный режим. Для сравнения с экспериментальными данными рассматривались экспериментальные данные по массопереносу этих смесей, представленные в работе [5] для пор с характеристиками: радиус поры 37 нм, длина поры 70 мкм, при условиях: давление 0,02 атм, температура 298 К. Фактор разделения для смеси Ar - H₂, рассчитанный на основе математического моделирования, составил 1,98, а его экспериментальное значение, приведенное в работе [5], равнялось 1,9 (относительная ошибка расхождения расчетных и экспериментальных данных составляет менее 5 %). Относительная ошибка между фактором разделения для смесей CO₂ - CH₄ (представленных выше составов), рассчитанных на основе математической модели и найденных в результате экспериментальных исследований [5], не превосходит 10 %. Результат сопоставления для различных смесей свидетельствует об адекватности представленной математической модели разделения смеси.

На рис. 4 представлены кривые изменения числа столкновений молекул со стенкой поры и друг с другом. С увеличением расстояния от начала поры в ней накапливаются молекулы с более низкой скоростью теплового движения, и из-за увеличения численности молекул в центральной части длины канала происходит увеличение столкновения молекул друг с другом, что отражено на рис. 5. Падение числа столкновений молекул друг с другом после прохождения центральных участников канала происходит за счет уменьшения общей численности молекул к концу канала из-за падения давления на этих участках. Средняя скорость теплового движения молекул имеет наибольшее значение в начальных зонах канала, что приводит к более усиленному числу столкновений на начальных участках канала, что отражено на рис. 4. В результате моделирования разделения смесей СО₂ - СН₄ показано, что уменьшение радиуса поры и увеличение длины поры приводит к существенному увеличению фактора разделения смесей.

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

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

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

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

• Чулок А.И., д.т.н., профессор ГОУ ВПО МО «Академия социального управления», г. Москва;
• Куркина Е.С., д.ф.-м.н., профессор, в.н.с. факультета ВМК Московского государственного университета им. М.В. Ломоносова, г. Москва.

Работа поступила в редакцию 16.12.2011.