Поверхность. Рентгеновские, синхротронные и нейтронные исследования, 2020, № 3, стр. 95-102

ВВЕДЕНИЕ

Экспериментальные и теоретические исследования диффузионных свойств пористых веществ и материалов имеют большое практическое значение в различных областях науки и технологий. Диффузионные характеристики, включая коэффициент диффузии и проницаемость, существенно зависят от таких свойств материала, как пористость, средний размер и связность пор, локализация пор по объему образца, распределение по величине сечения пор, степень извилистости пор и других [1]. Отличительной особенностью пористых систем является неупорядоченность их структуры [2]. Точное аналитическое решение газодинамических уравнений, описывающих процессы взаимодействия пара или газа с пористой средой, практически невозможно. Поэтому при теоретическом исследовании заменяют реальную пористую среду упрощенными упорядоченными моделями с эквивалентными гидравлическими свойствами. Точные решения уравнений газовой динамики и гидродинамики могут быть получены для случая течения вязкой жидкости или плотного газа по цилиндрической трубе соответствующих размеров. Данное обстоятельство было использовано при конструировании моделей для описания диффузии вязких жидкостей и плотных газов [3–5]. Другой подход к описанию процессов взаимодействия газообразного вещества с пористыми материалами основан на решении, например, нелинейного уравнения диффузии со ступенчатой аппроксимацией полученного в результате экспериментальных измерений коэффициента диффузии [6]. В связи с интенсивным развитием нанотехнологий [7] широкое распространение получили материалы с равномерно распределенными порами заданного поперечного сечения [8], что привело к значительному повышению актуальности развития моделей, основанных на использовании индивидуальной поры соответствующей геометрической конфигурации для исследования диффузионных характеристик взаимодействия пористого материала с газом или паром. Полученные с использованием упомянутых моделей диффузионные характеристики могут быть использованы в дальнейшем для исследования крупных образцов пористых материалов [9, 10]. Данные исследования особенно актуальны для решения практических задач фильтрации в экологии, защите окружающей среды от воздействия потенциально опасных, с точки зрения химического и радиоактивного загрязнения, объектов промышленного и научно-исследовательского назначения. Поскольку используемые в научных исследованиях, строительстве и промышленном производстве пористые материалы содержат в себе микро- и мезопоры, очень редко макропоры, при рассмотрении физической картины процессов взаимодействия газа с порой или стенками поры должна учитываться молекулярная структура протекающего через пору газа, так как длина свободного пробега молекулы газа не превышает размер индивидуальной поры [11]. Поэтому в настоящей работе моделирование процессов взаимодействия водяного пара с порой проведено методом молекулярной динамики [12]. В качестве модели была использована сквозная пора цилиндрического типа. Было проведено численное исследование зависимости величины проницаемости индивидуальной цилиндрической поры от ее радиуса, длины и характеристик процесса взаимодействия молекул воды со стенками поры.

МОЛЕКУЛЯРНО-ДИНАМИЧЕСКАЯ МОДЕЛЬ ВЗАИМОДЕЙСТВИЯ ВОДЯНОГО ПАРА И СКВОЗНОЙ ПОРЫ

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

Здесь i – номер частицы $\left( {1 \leqslant i{\text{\;}} \leqslant N} \right),$ N – полное число частиц, ${{m}_{i}}$ – масса частицы, ${{{\mathbf{f}}}_{i}}$ – равнодействующая всех сил, действующих на частицу: где U – потенциал взаимодействия между частицами, ${{{\mathbf{f}}}_{i}}^{{{\text{ex}}}}$ – сила, обусловленная внешними полями.

Для моделирования диффузии водяного пара через сквозную пору была использована трехмерная модель. Молекулы воды – материальные точки, движение которых описывается классическими уравнениями Ньютона (1). Модель поры представляет собой цилиндр длиной $L$ и диаметром $d$ (рис. 1). Боковая поверхность цилиндра непрозрачна для молекул воды. Основания цилиндра открыты для движения молекул в обе стороны. Внешнее по отношению к поре пространство моделируется прямоугольными параллелепипедами с размерами $d \times d \times L,$ примыкающими квадратными гранями к обоим основаниям цилиндра.

При моделировании методом молекулярной динамики использовали следующие граничные условия. Боковая поверхность цилиндра Γ_p = $ = \left[ {{{x}^{2}} + {{y}^{2}} = {{r}^{2}};\,\,0 \leqslant z \leqslant L} \right]$, где r – радиус основания цилиндра, абсолютно непрозрачна для молекул воды. Основания цилиндра абсолютно прозрачны для молекул воды. Зафиксированы периодические условия на границах ${\Gamma }_{1}^{{lo}},{\Gamma }_{2}^{{lo}},{\Gamma }_{3}^{{lo}},$ ${\Gamma }_{4}^{{lo}},{\Gamma }_{5}^{{lo}},~{\Gamma }_{6}^{{lo}}$ и ${\Gamma }_{1}^{{ro}},{\Gamma }_{2}^{{ro}},~{\Gamma }_{3}^{{ro}},{\Gamma }_{4}^{{ro}},{\Gamma }_{5}^{{ro}},{\Gamma }_{6}^{{ro}},$ которые являются, соответственно, гранями левого и правого параллелепипедов и моделируют внешнее по отношению к поре пространство.

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

с постоянными $~\varepsilon = 6.74 \times {{10}^{{ - 3}}}$ эВ и $\sigma = 3.17~\,\,{\text{{\AA}}}{\text{.}}$ Здесь r – расстояние между центрами частиц, ε – глубина потенциальной ямы, σ – расстояние, на котором энергия взаимодействия становится равной нулю. Параметры ε и σ являются характеристиками атомов соответствующего вещества. Минимум потенциала находится в точке ${{r}_{{{\text{min}}}}} = 2{{\sigma }^{{{1 \mathord{\left/ {\vphantom {1 6}} \right. \kern-0em} 6}}}}.$ Потенциал Леннарда-Джонса широко используется при молекулярно-динамическом моделировании воды и насыщенного водяного пара [14].

Взаимодействие молекул воды со стенками поры моделируется также с использованием потенциала Леннарда-Джонса, но с параметрами ${{\sigma }_{{wf}}}$ и ${{\varepsilon }_{{wf}}},$ полученными для учета взаимодействия молекул газа с молекулами стенки [15]:

Здесь ${{\sigma }_{{wf}}}$ и ${{\varepsilon }_{{wf}}}$ – параметры длины и энергии, используемые для учета взаимодействия молекул воды с молекулами материала поры без учета молекулярной структуры стенок поры. Переменная ξ представляет собой расстояние от молекулы воды до стенки поры. В случае цилиндрической модели поры ξ – это расстояние от молекулы воды до боковой поверхности цилиндра (рис. 1).

Для интегрирования уравнений движения (1) был использован модифицированный метод Верле [16]. Интегрирование проводится по следующей схеме. В начале каждого шага задают или рассчитывают на предыдущем шаге по времени t следующие величины: ${\mathbf{r}}\left( t \right),~{\mathbf{v}}\left( t \right),~{\mathbf{f}}\left( t \right).$ Затем вычисляют скорости частиц в момент времени $t + \frac{{\Delta t}}{2}$ и координаты нового местонахождения частиц:

После этого пересчитывают силы, действующие на частицу в момент времени $t + $ $\Delta t{\kern 1pt} :$Далее рассчитывают значения скоростей на следующем шаге:

Для контроля температуры при обмене энергии между моделируемой системой и окружающей средой в молекулярной динамике используют специальные алгоритмы, называемые термостатами. В настоящей работе для контроля температуры моделируемой системы был использован термостат Берендсена [17]. Он реализуется в уравнениях движения (1) путем включения в результирующую силу дополнительного переменного нелинейного трения. В результате уравнения движения преобразуются к следующему виду:

Коэффициент пересчета скорости в каждый момент времени $t$:

зависит от так называемого “времени нарастания” ${{\tau }_{{\text{B}}}}$ термостата, значение которого находится в закрытом интервале [0.1, 2] пс. Параметр ${{\tau }_{{\text{B}}}}$ описывает силу связи моделируемой системы с гипотетическим тепловым резервуаром, моделирующим бесконечную внешнюю среду. При увеличении ${{\tau }_{{\text{B}}}}$ сила связи системы с тепловым резервуаром уменьшается, что приводит к увеличению времени достижения заданной температуры ${{T}_{0}}$ в системе, если в данный момент времени $t$ система находится в термодинамическом состоянии при температуре T(t). Алгоритм Берендсена, будучи достаточно простым в численной реализации, очень эффективен для достижения целевой температуры в состояниях системы, далеких от термодинамического равновесия.

Для моделирования учета взаимодействия молекул воды со стенками поры вплоть до ее конденсации на стенках был использован термостат Андерсена [18]. Помимо возможности моделировать процесс прилипания молекул воды к стенкам поры термостат Андерсена (один из самых простых термостатов) позволяет корректно моделировать канонический ансамбль с сохранением числа частиц ($N$), объема системы ($V$) и температуры ($T$). Он действует на молекулы воды в тонком слое толщиной $\Delta r$ около стенки поры (рис. 1). В каждый момент времени скорости (импульсы) всех частиц, попавших в упомянутый выше слой, приравниваются случайной величине с равномерным распределением при заданной температуре. Данный алгоритм имитирует столкновение молекул воды с частицами теплового резервуара при определенной температуре T. Сила сцепления с тепловым резервуаром определяется частотой столкновений ν. Для каждой частицы случайная переменная, пропорциональная числу столкновений, выбирается из интервала от нуля до единицы. Если значение этой переменной меньше, чем $\nu \Delta t$, где $\Delta t$ – временнóй шаг, то скорость частицы ${\mathbf{v}}\left( {t + \Delta t} \right)$ становится равной случайному числу $\zeta \in N\left[ {0;T} \right]$, где $N\left[ {0;T} \right]$ – распределение Гаусса с дисперсией Т. В этом случае скоростная схема Верле с учетом термостата Андерсена может быть представлена в следующем виде: для момента времени $t + ~\,\Delta t$ вычисляются координаты положения частицы

затем вычисляются значения скоростей частиц в момент времени $t + \frac{{\Delta t}}{2}$: после этого вычисляется скорость частицы в момент времени $t + \Delta t$: Скорость частицы выбирают в зависимости от частоты столкновений ν. Если случайно выбранная величина $\xi < \Delta t,$ $\xi \in U\left[ {0;1} \right],$ где $U\left[ {0;1} \right]$ – равномерное распределение вещественной случайной величины, то $\left| {{\mathbf{v}}{\text{\;}}\left( {t + \Delta t} \right)} \right| = \zeta {\text{,}}$ где случайная величина $\zeta \in N\left[ {0;T} \right],$ $N\left[ {0;T} \right]$ – распределение Гаусса с дисперсией Т. Для вычисления давления в поре используется формула, следующая из вириального уравнения [19]: Здесь $V$ – объем поры, $\left\langle {2K} \right\rangle $ – удвоенная кинетическая энергия, усредненная по ансамблю частиц, ${\mathbf{f}}\left( {{{{\mathbf{r}}}_{{ij}}}} \right)$ – сила взаимодействия между частицами $i$ и $j$ на расстоянии ${{{\mathbf{r}}}_{{ij}}}$.

ОПИСАНИЕ МОДЕЛИРУЕМОЙ СИСТЕМЫ

В настоящей работе исследование процессов диффузии водяного пара через пористую поверхность проводили путем моделирования потока молекул воды через сквозную цилиндрическую пору объемом $V = {{L\pi {{d}^{2}}} \mathord{\left/ {\vphantom {{L\pi {{d}^{2}}} {4,}}} \right. \kern-0em} {4,}}$ где $L = 400$ нм – высота цилиндра, $d = 40$ нм – диаметр его основания. Моделирование осуществляли при температуре окружающей среды Т = 35°С слева и справа от поры и определенных давлении и концентрации на открытых концах поры. Слева от поры (рис. 1) находится 100% насыщенный пар при атмосферном давлении ${{p}_{L}} = 5.622~\,\,{\text{кПа,}}$ концентрация ρ_L = $ = 0.03962~\,\,{\text{кг/}}{{{\text{м}}}^{3}}.$ Справа от поры – 20% насыщенный водяной пар при атмосферном давлении ${{p}_{R}} = 1.124~\,\,{\text{кПа}},$ ${{\rho }_{R}} = 7.9~\, \times {{10}^{{ - 3}}}{\text{\;кг/}}{{{\text{м}}}^{3}}.$ В приведенных единицах имеем следующие значения: $L* = 1264,$ $d* = 126.4,$ $T* = 3.94,~$ $p_{L}^{*}$ = = $1.65183\,\, \times \,\,{{10}^{{ - 4}}},$ $p_{R}^{*} = 3.3\,\, \times \,\,{{10}^{{ - 5}}}~.$

В настоящей работе при исследовании взаимодействия поры и водяного пара выбирали начальные и граничные условия на открытых концах поры. В начальный момент времени давление и концентрация на левом и правом концах поры равны ${{p}_{L}},$ ${{\rho }_{L}}$ и ${{p}_{R}},$ ${{\rho }_{R}}$ соответственно. Внутри поры водяной пар находится в состоянии термодинамического равновесия. То есть оба конца поры не изолированы от внешней среды. В этом случае давление вдоль оси z зависит линейно от z:

Разность значений давления насыщенного пара на левом и правом концах поры приводит к возникновению внешней силы, действующей на каждую молекулу воды внутри поры. Эта сила также является линейной функцией координаты $z.$ В общем случае модуль силы определяется как $F = pS,$ где $p$ – давление, действующее на плоскую поверхность с площадью $S.$ В случае частиц внутри поры площадь, к которой приложено давление, равна поперечному сечению молекул воды $S = {{\pi {{d}^{2}}} \mathord{\left/ {\vphantom {{\pi {{d}^{2}}} 4}} \right. \kern-0em} 4},$ или в приведенных единицах $S* = {\pi \mathord{\left/ {\vphantom {\pi {4,}}} \right. \kern-0em} {4,}}$ и внешняя сила может быть представлена в следующем виде:

Таким образом, равнодействующая всех сил, действующих на частицу внутри поры, будет включать в себя не только векторную сумму всех сил, порождаемых потенциалом Леннарда-Джонса межмолекулярного взаимодействия (3) и потенциалом взаимодействия со стенками поры (4), но и внешнюю силу, обусловленную разницей давлений и градиентом концентраций между открытыми левым и правым концами поры. В работе при численном моделировании в качестве материала стенки поры использовали диоксид кремния с параметрами $\sigma = 0.42{\text{\;нм}}$ и $\varepsilon = 32.59 \times {{10}^{{ - 23}}}\,\,{\text{Дж,}}$ или в приведенных единицах $\sigma * = 1.33$ и $\varepsilon * = 0.3.$ Применяя процедуру Лоренца–Бертло [20–22], можно получить следующие значения для параметров Леннарда-Джонса (4) при взаимодействии молекул воды с молекулами материала стенки поры: $\sigma _{{wf}}^{*} = 1.16$ и $\varepsilon _{{wf}}^{*} = 0.55.$ Вычисления проводили для $2 \times {{10}^{6}}$ шагов по времени с шагом $\Delta t* = 0.005,$ что в абсолютных единицах соответствует 16 663 пс. Параметр $\tau _{{\text{B}}}^{*}$ термостата Берендсена принимали равным $\tau _{{\text{B}}}^{*} = 0.05,$ параметр ν термостата Андерсена – единице. Взаимодействие со стенкой поры учитывали для частиц, попавших при движении внутри поры в слой толщиной $\Delta r = 0.025d$ вблизи стенки поры.

РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ И ИХ ОБСУЖДЕНИЕ

Моделирование проводилось для нескольких вариантов цилиндрической поры длиной L = 200, 400, 800 нм и диаметром d = 20, 30, 40 нм. Как видно из рис. 2, коэффициенты диффузии $D$ и проницаемости $K$ резко уменьшаются во временнóм интервале [0; 833] пс. Затем они очень медленно эволюционируют к некоторым значениям, обусловленным начальными и граничными условиями задачи. Такой характер эволюции $D$ и $K$ не зависит от размеров поры. Качественно поведение усредненных по объему поры коэффициентов диффузии и проницаемости можно объяснить, если проанализировать процесс конденсации молекул воды на стенках поры (рис. 3). К моменту времени $t = 833~\,\,{\text{пс}}$ бóльшая часть частиц ($ \approx {\kern 1pt} 96\% $), находящихся внутри поры, конденсируется на ее стенках. Поэтому лишь малое количество частиц, находящихся в данный момент внутри поры, дает вклад в процесс диффузии. Именно поэтому (рис. 2) резко падают усредненные коэффициенты диффузии и проницаемости в интервале [0; 833] пс.

Результаты визуализации эволюции во временнóм интервале [833; 16 663] пс усредненного по объему коэффициента диффузии $D\left( t \right)$ при течении водяного пара через пору представлены на рис. 4 и 5. Как видно из рисунков, с ростом диаметра поры $d$ увеличивается зависимость коэффициента диффузии $D\left( t \right)$ от длины поры $L$ (рис. 4). Аналогичная картина наблюдается для зависимости коэффициента диффузии $D\left( t \right)$ от диаметра поры $d$ при увеличении длины поры $L$ (рис. 5).

Результаты визуализации эволюции во временнóм интервале [833; 16 663] пс усредненного по объему коэффициента проницаемости $K\left( t \right)$ при течении водяного пара через пору представлены на рис. 6 и 7. Рисунки демонстрируют зависимость проницаемости $K\left( t \right)$ от времени, качественно отличающуюся от зависимости $D\left( t \right)$. C ростом диаметра поры $d$ зависимость коэффициента проницаемости $K\left( t \right)$ от длины поры $L$ уменьшается (рис. 6). Аналогичная картина наблюдается для зависимости коэффициента проницаемости $K\left( t \right)$ от диаметра поры $d$ при увеличении длины поры $L$ (рис. 7).

ЗАКЛЮЧЕНИЕ

В работе проведено исследование методом классической молекулярной динамики эволюции во времени коэффициента диффузии $D\left( t \right)$ и проницаемости $K\left( t \right)$ в процессе взаимодействия водяного пара и сквозной цилиндрической поры в случае ненулевой разности давлений на открытых концах поры. Построена молекулярно-динамическая модель, учитывающая процесс конденсации молекул воды на стенках поры, что позволяет моделировать степень гидрофильности/гидрофобности внутренней поверхности поры в зависимости от химического состава ее материала. Методами численного моделирования установлена зависимость коэффициента диффузии $D\left( t \right)$ и проницаемости $K\left( t \right)$ от длины и диаметра поры. Показано, что скорость резкого уменьшения коэффициента диффузии $D\left( t \right)$ и проницаемости $K\left( t \right)$ на начальных этапах эволюции $t \in \left[ {0;833} \right]~\,\,{\text{пс}}$ обусловлена тем, что большая часть частиц ($ \approx {\kern 1pt} 96\% $), находящихся в данный момент внутри поры, конденсируется на ее стенках. Данный эффект полностью определяется характером взаимодействия молекул воды и атомов материала внутренней поверхности поры. Результаты данного исследования могут быть использованы для моделирования тепломассопереноса в однородных пористых материалах в различных областях науки и технологий, связанных с производством строительных материалов и защитой окружающей среды.