CC BY
89
УДК 532.137.2; 532.5.013.12
МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ В НАНОКАНАЛАХ МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ
В. Я. Рудяк1, А. А. Белкин1, В. В. Егоров1, Д. А. Иванов1
Новосибирский Государственный Архитектурно-Строительный Университет (Сибстрин),
Новосибирск, Россия
valery.rudyak@mail.ru
PACS 47.11.Mn 47.67.-k 62.25.-g
На основе метода молекулярной динамики предложен алгоритм, впервые позволяющий моделировать плоское течение флюида в наноканалах с перепадом давления. Взаимодействие молекул флюида моделировалось потенциалом твердых сфер или потенциалом Леннард-Джонса. Изучены свойства нанотечений. Показано, что структура флюида в наноканале существенно отличается от его структуры в объеме. Представлены данные о зависимости коэффициента трения флюида на стенке от чисел Кнудсена и Рейнольдса. Установлено, что падение давления существенно зависит от коэффициентов аккомодации (для флюида твердых сфер) или от параметров взаимодействия молекул стенки с молекулами флюида (для флюида Леннард-Джонса).
Ключевые слова: наноканалы, падение давления, сопротивление трения, метод молекулярной динамики, структура жидкости.
1. Введение
Активное изучение микротечений жидкости и газа, наблюдающееся в последнее двадцатилетие, мотивировано появлением большого числа микроэлетромеханических систем (МЭМС), а позднее и нанотехнологий. Данная тематика столь актуальна, что уже существуют специализированные монографии. И если в конце 80-х годов, когда интерес к этой тематике стал очевиден, основные ее приложения сводились к созданию различных МЭМС и биохимических систем типа lab-on-chip, то, уже начиная с 90-х, происходит ее диверсификация. Шире становится область применения, в том числе в медицине, фармакологии, биологии, теплоэнергетике, приборостроении, катализе и т.д. Уже в этом столетии проявился активный интерес и к нанотечениям. Связано это не только с созданием нано-технологий различного назначения, но и с исследованиями в достаточно традиционных областях: биологии, геофизике, теплоэнергетике и т.д. До сих пор одним из самых загадочных процессов является перенос питательных веществ в растениях и живых организмах. Активно проводятся исследования по созданию различного рода микропористых покрытий и течений в них. Актуальность для России нефте- и газодобычи сегодня понятна всем, а типичные размеры пор в несущих формациях меняются от десятков нанометров до десятков или даже сотен микрон [1]. Таким образом, и здесь имеют место микро- и нанотечения.
Как известно, в обычных условиях течения жидкостей и не слишком разреженных газов вполне можно описывать методами механики сплошной среды. Однако в микроканалах ситуация существенно меняется. Микроканалами обычно называют каналы, один из характерных размеров которых h (например, высота плоского канала или диаметр цилиндрического) оказывается порядка от 1 до примерно 300 мкм. В этих условиях течения жидкости и газа, как правило, следует описывать по-разному. Действительно, если газ не слишком плотный (до давлений примерно 10-20 атм.) соответствующее число Кнудсена
Кп таких микротечений изменяется в пределах: 10_2 ^ Кп ^ 102. В этом диапазоне чисел Кнудсена течение уже не описывается уравнениями гидродинамики. Точнее, на нижнем пределе все еще можно использовать уравнения Навье-Стокса, но с граничными условиями скольжения. Затем примерно до чисел Кнудсена Кп ~ 10_1 необходимо применять уравнения Барнетта, а потом - кинетическое уравнение Больцмана. При этом стоит иметь в виду, что к использованию уравнений Барнетта необходимо относиться с осторожностью. Строго говоря, они не полные и не учитывают эффектов памяти [2].
Концепция сплошной среды для жидкостей работает, если можно выделить гидродинамический физически бесконечно малый масштаб г/,, чтобы флуктуациями внутри соответствующего объема можно было пренебречь. Для жидкости г/, ~ \fdfi [2], где (1 -характерный размер молекулы жидкости. Если микроканал имеет высоту Н ~ 1 мкм, то гь ~ 2 • 10"6 см, что уже вполне сопоставимо с высотой канала. Поэтому при наличии в таком течении градиентов макроскопических переменных гидродинамическое описание будет давать сбои. Трудности вызывает даже просто введение макроскопических переменных, которые по определению есть средние по физически бесконечно малому объему соответствующих динамических переменных молекулярной системы. Таким образом, начиная примерно с 50 мкм, для моделирования микротечений нельзя применять обычные гидродинамические методы. Какова альтернатива?
При описании газов следует, прежде всего, иметь в виду, что в микроканалах необходимо анализировать несколько параметров подобия. Типичной является ситуация, когда число Кнудсена по ширине равно Кп ^ 10_1, а по высоте — Кп ^ 1, в наноканалах — Кп ^ 10. Казалось бы, для моделирования можно использовать метод прямого статистического моделирования Монте-Карло (ПСМ) [3]. Однако скорости течений в микроканалах обычно не велики, а в этих условиях метод ПСМ работает не удовлетворительно. Реального продвижения можно ожидать, решая для данных задач полное уравнение Больцмана или применяя метод молекулярной динамики (МД). В последнее время достаточно активно развивается метод, в котором для описания течения газа в микроканале используется то или иное модельное кинетическое уравнение.
При моделировании течений плотных газов и жидкостей, как уже указывалось, для достаточно малых микроканалов нельзя использовать обычный гидродинамический подход. Фактически единственным методом моделирования, не вызывающим концептуальных возражений, является метод МД. Активное изучение свойств микротечений этим методом началось с середины восьмидесятых годов прошлого века (см., например, [4] и цитируемую там литературу), за это время было получено достаточно много интересных результатов. В основном объектом моделирования являлись микроскопические аналоги хорошо известных течений Куэтта и Пуазейля.
Течение Куэтта, как правило, моделируется с помощью ячейки, две границы которой представляют собой движущиеся в противоположные стороны пластины. Между ними находятся молекулы жидкости. Для такой системы были получены данные о профилях скорости и других характеристиках течения (см., например, [5]). Говоря о соответствии модели течения Куэтта и реальных течений такого рода, отметим, что скорости сдвига в МД моделировании в сотни раз превосходят реально достижимые, что связано с ограниченным быстродействием современных компьютеров. В остальном МД технология моделирования данного течения соответствует реальному течению Куэтта.
В то же время при МД моделировании течения Пуазейля для его генерации вводится некоторая фиктивная сила (см., например, [4-7] и цитированную там литературу). Часто эту силу называют гравитационной, но она на порядки превышает последнюю. С другой стороны, наличие постоянно действующей силы порождает ускорение молекул в
природе отсутствующее, поэтому приходится использовать различные методы коррекции их скоростей (так называемый термостат и т.п., см., например, [4,7]). В результате использования столь искусственных процедур так и не удается смоделировать реальное течение, возникающее в канале под действием перепада давления.
В данной работе предлагается новый МД алгоритм, позволяющий моделировать реальное плоское течение типа Пуазейля, характеризуемое определенным градиентом давления, и обсуждаются особенности плоских течений в наноканалах.
2. Методика моделирования течения
В данной работе моделируется плоское течение (вдоль оси х, см. рис. 1) молекулярного флюида в канале длиной Ь между двумя параллельными пластинами, расстояние между пластинами равно к. Ячейка моделирования представляет собой прямоугольный параллелепипед, нижняя и верхняя грани которого (перпендикулярные оси г) являются стенками канала. Размер ячейки вдоль оси у равен Ь ^ к ив этом направлении используются обычные периодические граничные условия.
Рис. 1. Схема моделирования течения в плоском канале с псевдопериодическими граничными условиями
Чтобы организовать течение флюида, был разработан алгоритм, использующий специальным образом модиф ицированные периодические граничные условия на левой и правой гранях ячейки (перпендикулярных оси х). Молекулы, находящиеся внутри канала, не могут пересекать левую грань ячейки (например, молекула 1 на рис. 1), взаимодействие с ней задаётся зеркальным или диффузным граничными условиями (см. ниже). Правую грань молекулы могут пересекать, при этом для молекулы, проходящей правую границу (например, молекула 2a на рис. 1), создается копия на левой (2Ь на рис. 1). Координаты х центров молекулы и копии отличаются на Ь, остальные их координаты и скорости одинаковы. При решении уравнений движения такой молекулы учитываются ее взаимодействия с молекулами, находящимися как у левой границы, так и у правой, что соответствует периодическим граничным условиям. Когда молекула, пересекает правую грань ячейки и полностью выходит из канала, ее копия у левой грани остается. Её скорость разыгрывается случайным образом в соответствии с распределением Максвелла
2 *
I
у
где Т — температура левой стенки, к — постоянная Больцмана, т — масса молекулы, V— компонента скорости, V — модуль скорости. При этом проекция скорости молекулы на ось х всегда задается положительной. Таким образом, левую границу канала можно рассматривать как источник молекул, скорости которых распределены по закону (1). Полное число молекул в ячейке в данном алгоритме не меняется, что удобно для его программной реализации.
Для описания взаимодействия молекул флюида использовались потенциал Леннард-Джонса (ЛД)
(2)
или потенциал твердых сфер (ТС)
а- (3)
где а, £ — параметры потенциала ЛД, а — диаметр твердой сферы, моделирующей молекулу.
Взаимодействие молекул флюида со стенкой задавалось двумя различными способами, в зависимости от вида используемого межмолекулярного потенциала. В первом случае, для течения флюида ЛД молекул, каждая стенка моделировалась двумя рядами неподвижных молекул, расположенных в узлах кубической гранецентрированной решетки (принципиально возможно реализовать любую укладку). Взаимодействие молекул флюида с молекулами стенок задается потенциалом (2), параметры которого определялись с помощью комбинационных соотношений: а\2 = ^о\\о22, £12 = у/£ц£22, где <?п, £ц ~ параметры взаимодействия молекул флюида, а а22, 22 —молекул стенок. В расчетах, приведенных в
о
этой работе, изучалось течение аргона: а11 = 3.405 А, £11/к = 119.8 К [8]. Параметры
о
молекул стенок соответствовали углероду: а22 = 3.4 А, е22/к = 28 К [9].
Для течения флюида твердых сфер верхняя и нижняя стенки канала моделировались двумя твердыми параллельными поверхностями. Взаимодействие молекул флюида со стенками происходило мгновенно, их скорость после соударения со стенкой канала определялась либо зеркальным законом, либо диффузным, либо зеркально-диффузным. При зеркальном отражении молекула изменяет знак проекции скорости на нормаль к стенке , проекции ух и ьу не меняются. При диффузном отражении компоненты скорости молекулы разыгрываются по закону (1) с соответствующей температурой. Наконец, при зеркально-диффузном отражении доля молекул в взаимодействовала диффузно, а (1 — в) — зеркально. Отметим, что зеркальные граничные условия соответствуют нулевому значению коэффициента аккомодации импульса стенки, коэффициент аккомодации диффузных стенок равен единице.
В начальный момент времени молекулы флюида размещались в ячейке моделирования равномерно. Их плотность определялась параметром Ван-дер-Ваальса еу = пс13 для ТС флюида и еу = па3для ЛД флюида, где п — числовая концентрация молекул. Скорости молекул в начальный момент задавались согласно распределению Максвелла (1). Затем рассчитывалась эволюция системы. Для системы ЛД молекул использовалась схема Шофилда, для ТС молекул — стандартный алгоритм МД моделирования [10]. Шаг интегрирования по времени равнялся = 10_12с. Радиус действия ЛД потенциала во всех расчетах равен 2.35 а. После завершения релаксационных процессов, время которых зависит от размеров канала, плотности жидкости и условий взаимодействия со стенками, в канале устанавливается течение жидкости. Далее начинается расчёт нужных характеристик (полей концентрации, скорости, давления и т.д.).
Для определения давления ЛД молекул использовалась формула, являющаяся следствием теоремы вириала [11]
Для ТС флюида потенциал взаимодействия и сила Fij сингулярные. В этом случае вириал силы (г^ • Fij) задается импульсом, переданным в столкновениях, и давление определяется формулой
Здесь г*^ — расстояние между молекулами, в момент столкновения оно равно их диаметру, Дрц — изменение модуля импульса г-ой молекулы за время ДЬ за счет ее столкновения с 2-ой.
3. Характерные черты нанотечений в плоском канале
Описанный выше МД алгоритм впервые позволил смоделировать реальное плоское течение, в котором имеет место перепад давления вдоль канала. В процессе эволюции молекулярного флюида в канале устанавливается стационарное течение. Его характеристики, их зависимость от геометрии канала, характера взаимодействия молекул со стенками, структура флюида в канале будут описаны в данном разделе. Во всех случаях изучались течения в наноканалах. Высота канала варьировалась от 6 до 50 диаметров молекул, длина — от 60 до 250 диаметров, а ширина —от 6 до 20 диаметров (для ЛД молекул в качестве их диаметра выбирался параметр а). Таким образом, самый длинный канал имел длину всего лишь чуть больше 70 нм, а его высота варьировалась от 2 до 15 нм. Плотность флюида изменялась в достаточно широких пределах: £у = 0.0014 ^ 0.88.
3.1. Профиль скорости течения
Алгоритм моделирования, описанный в разделе 2, предполагает наличие входного и выходного участков течения, как это имеет место и в реальных течениях. Стационарное течение устанавливается на некотором расстоянии от входа (левой границы канала, представленного на рис. 1). Формируемый в течении профиль скорости обусловлен взаимодействием молекул флюида со стенкой. Для ЛД флюида во всех случаях формируется параболический профиль скорости вида
где А — некоторая константа, а 8 — так называемая длина скольжения, определяемая соотношением и(г = 0, к) = 8 (ди/дг)\г=0 н.
Эффект скольжения наблюдался и ранее в известных работах по МД моделированию нанотечений ЛД флюида. Течение ТС флюида моделируется впервые. Здесь параболический профиль скорости вида (4) фиксируется, только если взаимодействие молекул со стенкой не является зеркальным. При зеркальном взаимодействии по всей длине канала наблюдается ударный профиль скорости.
В реальных макроскопических течениях жидкости обычно считается, что имеет место прилипание флюида на стенке, т.е. скорость течения на стенке равна скорости самой стенки. Скольжение фиксируется лишь для течений разреженного газа. В этом случае 8 ~ I ~ Кп, где I — длина свободного пробега молекулы [12]. Таким образом, обычно считается, что учет скольжения необходим, начиная с чисел Кп ~ 5 • 10_3. Поскольку длина
и = А(г2 — гк — 8к),
(4)
скольжения пропорциональна , то в макроскопических течениях не слишком разреженного газа им, естественно, можно пренебречь. Совершаемая ошибка при малых числах Кнудсена будет мала.
Природа появления скольжения в течениях жидкостей сложнее, чем в течениях газов. И если для газов существует систематическая теория, определяющая длину скольжения, то для жидкостей соответствующего аналога нет. Тем не менее, экспериментально длина скольжения в микротечениях последние годы очень активно изучается. Регистрируются длины скольжения от нескольких нанометров до примерно двадцати микрометров [13]. При использовании гидрофобных и в особенности так называемых ультрагидрофобных покрытий длина скольжения может быть еще больше [14,15]. Ясно, что наличие столь больших значений (впрочем, как и малых, порядка нескольких нанометров) длины скольжения нельзя объяснить с точки зрения кинетической теории.
При моделировании ЛД флюида с помощью описанного в разделе 2 алгоритма профиль скорости установившегося течения определяется лишь параметрами потенциала взаимодействия флюида с молекулами стенки. В проведенных расчетах длина скольжения уменьшалась с ростом плотности флюида. Поверхности моделировались системой молекул, расположенных в узлах кубической гранецентрированной решетки. Однако в общем случае характер взаимодействия молекул флюида со стенкой, а значит, и величина длины скольжения будет зависеть от типа кристаллической решетки твердых поверхностей, образующих канал. В частности, в работе [16] было установлено, что на длину скольжения оказывает влияние тип упаковки молекул стенок канала и угол ориентации потока жидкости относительно базисов образующей стенку кристаллической решетки. Таким образом, изменяя топологию стенок, можно регулировать сопротивление движению жидкости в наноканале.
Длина скольжения для течения ТС флюида зависела в первую очередь от коэффициента аккомодации вис его ростом уменьшалась. Так, если при 9 = 0.5 длина скольжения равнялась 1.7 д, (плотность флюида еу = 0.88), то при увеличении 9 до единицы длина скольжения уменьшалась более чем в три раза и равнялась 0.5 . Естественно, длина скольжения в общем случае увеличивается и с уменьшением плотности флюида. Так, при плотности флюида £у = 0.0014 она составляла около 10с1.
3.2. Падение давления
Обычное стационарное течение Пуазейля —это течение, формирующееся при заданном градиенте давления и, как следствие, с заданными значениями давления на входе и выходе. Падение давления вдоль канала обусловлено трением на стенках. До сих пор методом МД не удавалось смоделировать падение давления вдоль канала. Одним из достоинств предлагаемого алгоритма является реализация течения с линейным градиентом давления вдоль канала. Численные значения градиента давления зависят не только от геометрии канала (высоты и длины), но и в первую очередь от параметров взаимодействия флюида со стенками канала. На рис. 2 приведены типичные данные для падения давления при различных типах взаимодействия молекул флюида между собой и со стенками канала, давление нормировано на соответствующую величину на входе канала. Падение давления вдоль канала связано с сопротивлением, обусловленным взаимодействием молекул флюида со стенками. Поэтому при зеркальном отражении молекул ТС флюида давление вдоль канала остается неизменным (линия 1 на рис. 2). Величина градиента давления растет с увеличением коэффициента аккомодации (сравни линии 1 и 2).
Расход жидкости задается соотношением
н
Q = — J р г/:(л) ЬсЬ,
о
где Б, Н и Ь — соответственно площадь поперечного сечения, высота и ширина канала, р, и — массовая плотность и скорость флюида. Средняя скорость течения, определяемая расходом, равна: и = Ц/р. Она, как и расход, зависит от геометрии канала. Влияние размеров канала обусловлено особенностями используемого алгоритма. Течение в канале создается за счет условий, формулируемых на левой и правой границах канала. С увеличением высоты канала их площадь растет, естественно растет и доля молекул, скорость которых увеличивается при взаимодействии с этими границами. Наоборот, с увеличением длины канала относительное количество таких молекул снижается, что приводит к уменьшению средней скорости течения.
Рис. 2. Зависимость давления от продольной координаты канала. Линия 1 — ТС флюид, 9= 0, линия 2 — ТС флюид, в = 0.5, линия 3 — ЛД флюид. Ь = 60а, Н = 6а, £у = 0.79
3.3. Коэффициент сопротивления
Выявление отличий течения в наноканале от обычного гидродинамического чрезвычайно важно с практической точки зрения. Одной из важнейших характеристик течения является коэффициент гидравлического сопротивления, который определяется по формуле Дарси-Вейсбаха
А =
¿х рй2
В течении, реализованном посредством предложенного здесь МД алгоритма, не известна связь между градиентом давления и расходом. Поэтому сравнение коэффициента трения, полученного в МД моделировании, с гидродинамическим выполнялось при одинаковых перепадах давления или расходах. На рис. 3 представлены результаты сопоставления коэффициента трения, полученного в МД моделировании, с гидродинамическим коэффициентом при одинаковых перепадах давления. Здесь непрерывная кривая соответствует гидродинамическому коэффициенту сопротивления \н = 24/Яе, а МД данным — крестики. Пунктирная кривая получена по МД данным методом наименьших квадратов. Обе кривые качественно подобны, но гидродинамические значения несколько выше. Следует заметить,
Л
600
400
200
0
I-..........1-х
О 0.2 0.4 0.6 0.8 1 Ке
Рис. 3. Зависимость коэффициента сопротивления в наноканале от числа Рейнольдса. Непрерывная кривая —гидродинамическое значение, пунктирная кривая — МД данные
что эти значения получены только для малых чисел Рейнольдса. С другой стороны, именно эти числа Рейнольдса обычно и типичны для нанотечений.
Следующий график (рис. 4) показывает зависимость коэффициента сопротивления от числа Кнудсена. Здесь квадратики соответствуют коэффициенту аккомодации в = 1, а крестики — 0 = 0.5. Как и следовало ожидать, с ростом коэффициента аккомодации увеличивается и коэффициент сопротивления. В этой связи стоит отметить, что число Рейнольдса является более грубым параметром подобия и в некотором смысле более универсальным. Различным коэффициентам аккомодации просто будут соответствовать разные числа Рей-нольдса.
4. О структуре жидкости в плоском канале
Несмотря на отмеченные во введении недостатки предложенных ранее алгоритмов моделирования течений, с их помощью были получены важные результаты по изучению структуры жидкости в наноканале. Было установлено, что вблизи стенок канала наблюдаются максимумы профиля плотности жидкости. Ранее, однако, не моделировались течения ТС флюида. На рис. 5 проведено сопоставление профилей плотности в одинаковом канале для ТС и ЛД флюидов.
В обоих случаях плотность квазипериодически меняется поперек канала. В ТС флюиде максимумы и минимумы профиля плотности проявляется значительно сильнее, существенно выше оказываются и первые максимумы плотности у стенок. Это связано в первую очередь с более явно выраженными зонами экранирования для ТС молекул у поверхности стенок. Наличие этих зон приводит к тому, что эффективный объем, занимаемый ЛД молекулами, больше, чем занимаемый ТС молекулами. По этой причине эффективная плотность ТС флюида несколько выше, чем ЛД флюида. С ростом же плотности эффекты структурирования флюида увеличиваются.
Рис. 4. Зависимость коэффициента сопротивления в наноканале от числа Кнудсена
Р
Р
4
-1 1 1 1 1 1 Г -
1 I / 7
4 П Лх^/ ■1 ^ т X I ^ГГ^Т- ц ¿-
* | 1 1 1 1 V
Рис. 5. Профиль плотности в поперечном сечении канала ТС молекулы (+), ЛД молекулы (х). Ь = 60а, Ь= 6о, £у = 0.79
Рис. 6. Поле плотности в наноканале. Ь = 60й, к= 6й, £у = 0.79
Структурирование флюида в канале является принципиальным фактом. Однако вдоль канала характер этого структурирования меняется. Падение давления в изотермических условиях должно вызывать соответствующее уменьшение плотности флюида. На рис. 6 представлено поле плотности вдоль канала. Хорошо видно, что структура поля плотности заметно меняется. Если вначале канала жидкость хорошо структурирована, то степень этого структурирования вдоль канала снижается. Связано это с тем, что плотность вдоль канала уменьшается. Флюид оказывается существенно сжимаемым. Средний расход флюида остается постоянным в любом поперечном сечении канала, поэтому снижение плотности приводит к росту средней скорости по мере удаления от начала канала. В этом смысле реализуется очень любопытная ситуация. Грубо говоря, параметры течения меняются от сечения к сечению при том, что течение является стационарным, оно повторяемо и реализуется с периодом Ь/и. Таким образом, мы имеем специфическое стационарное течение. Изменение плотности вдоль канала связано, как уже отмечалось, с падением давления. Средняя плотность жидкости остается постоянной по всей длине канала лишь для случая зеркального взаимодействия со стенкой ТС-молекул.
Рис. 7. Профили плотности в наноканалах различной высоты £у
1 - к = 6а, 2 - к = 12а, 3 - к = 24а, 4 - к = 48а
0.88,
Важно подчеркнуть, что упорядоченность жидкости вблизи стенок является характерной чертой течений в наноканалах, она не исчезает при увеличении расстояния между стенками канала. Это иллюстрирует рис. 7, где представлены профили плотности поперек
каналов с разной высотой. Здесь кривая 1 соответствует высоте к = 6а, кривая 2 — к = 12а, кривая 3 — к = 24а и кривая 4 — к = 48а. Структурирование жидкости в канале практически перестает зависеть от его высоты при к > 10а и наблюдается на расстояниях порядка 5-6а.
Профиль плотности дает локальную информацию о структуре жидкости. Но анализ рис. 5-7 показывает, что в канале у жидкости меняется характер ближнего порядка. Это должно наблюдаться, по крайней мере, вблизи стенок канала. Характер структуры жидкости детектируется парной радиальной функцией распределения молекул
д2(г) = Ш (5)
где ¿М — число молекул в радиальном слое толщиной ¿г на расстоянии г от выделенной молекулы. Данная функция показывает, как распределена плотность молекул вокруг произвольной выделенной. Распределение плотности в наноканале неоднородное, поэтому, чтобы получить достаточно объективную информацию о структуре жидкости, радиальные функции распределения следует строить в слоях толщиной dh, параллельных стенкам канала и находящихся на разном расстоянии от них
^ ^ 2ттпг с1гс1к ^ ^
Рис. 8. Радиальная функция распределения д2 в объеме (сплошная линия) и в слое первого максимума плотности жидкости в наноканале (пунктирная линия), к = 6а, £у = 0.88
Характер получающейся радиальной функции распределения, построенной для первого максимума плотности (см. рис. 7), представлен на рис. 8 (пунктирная кривая). Здесь же для сопоставления приведена радиальная функция распределения для флюида той же плотности в объеме (непрерывная кривая). Расстояние на этом рисунке измеряется в единицах а. Сопоставление этих функций показывает, что вблизи поверхности характер структуры жидкости существенно меняется. Если в объеме ближний порядок проявляется на расстояниях
порядка 1 нм (3а-4а), то вблизи поверхности фактически имеет место квазидальний порядок. Затухание радиальной функции распределения происходит на расстояниях, по крайней мере, на порядок больших, чем в объеме. Парная функция распределения (6) определяется для молекул, находящихся на одинаковом расстоянии от границы канала, то есть в зоне практически одинаковой концентрации. Если говорить о структуре жидкости во всем канале, то здесь ситуация существенно более сложная. Упорядоченность молекул обусловлена двумя факторами: наличием ближней структуры вокруг каждой молекулы и взаимодействием жидкости со стенками.
5. Заключение
Алгоритм МД моделирования, предложенный в настоящей работе, впервые позволил смоделировать течение с градиентом давления вдоль канала, аналогом которого является стационарное течение Пуазейля. Вместе с тем реализованное течение является необычным стационарным течением, его параметры меняются вдоль канала. В частности, изменяется плотность флюида, он ведет себя подобно сжимаемому газу даже при плотностях соответствующих жидкости. Его сжимаемость, однако, проявляется только в том случае, когда имеет место падение давления. Оно в свою очередь наблюдается лишь при диффузном или зеркально-диффузном законе отражения молекул от стенки. При зеркальном отражении молекул от стенки профиль скорости вдоль канала не меняется. И если на входе он был ударным, то таким и остается. Падение давления в этом случае также отсутствует. Таким образом, падение давления обусловлено исключительно законом взаимодействия молекул флюида со стенкой. Потери давления в свою очередь индуцируют снижение плотности (температура флюида остается практически постоянной).
Характер течения ЛД и ТС флюидов качественно одинаков. В обоих случаях имеет место перепад давления в канале. Градиент давления в канале с ЛД флюидом определяется параметрами потенциала взаимодействия молекул флюида с молекулами стенки. Существуют параметры, когда характеристики течения обоих флюидов оказываются идентичными. Это дает реальный инструмент для определения коэффициентов аккомодации [17].
Взаимодействие со стенкой существенно меняет и структуру флюида в канале. Вблизи стенок, на расстоянии в несколько размеров молекул плотность флюида значительно выше средней. Кроме того, жидкость на этих расстояниях оказывается структурированной, а вблизи стенки имеет место даже квазидальний порядок. Варьирование высоты канала показывает, что характер этого структурирования от нее практически не зависит. Таким образом, вблизи поверхности существует наноразмерный слой флюида толщиной один —два нанометра, обладающий особыми свойствами. По-видимому, именно наличие этого слоя и определяет в идеальных условиях длину скольжения. Этот слой играет роль слоя Кнуд-сена в разреженном газе. Его наличие, однако, никак не объясняет реально наблюдаемые в эксперименте длины скольжения в десятки и даже сотни нанометров. Из-за чего возникают столь большие длины скольжения? Если отвлечься от возможных «неидеальностей» реального флюида, например, его газонасыщенности, вследствие чего в некоторых случаях возникает «газовая смазка», то наиболее реальной причиной может явиться наличие неоднородностей на стенке. Фиксируемые шероховатости на твердых поверхностях имеют размер от нескольких нанометров. Безусловно, это будет менять характер взаимодействия со стенкой. И изучение этого вопроса является одной из наиболее актуальных задач физики микротечений.
Подводя итог, следует отметить, что разработанный для моделирования микротечений МД алгоритм не лишен недостатков. Главный из них состоит в том, что алгоритм не позволяет задавать непосредственно градиент давления, хотя его удается предсказать с
неплохой точностью по известным размерам канала, свойствам жидкости и стенок. Еще одним недостатком является нефизичное дальнодействующее взаимодействие молекул, находящихся у выхода из канала, с молекулами на его входе. Это влияние не столь драматично, поскольку имеет место лишь в узких релаксационных зонах, на выходе из которых собственно и реализуется течение с градиентом давления.
Работа выполнена при частичной поддержке РФФИ (грант № 10-01-00074) и ФЦП «Научные и научно-педагогические кадры инновационной России на 2009 — 2013 годы» (Гос. контракты П230, 14.740.11.0579).
Литература
[1] Nelson P.H. Pore-throat sizes in sandstones, tight sandstones, and shales // AAPG Bulletin. — 2009. — V. 93(3). — P. 329-340.
[2] Рудяк В.Я. Нелокальное решение уравнения Больцмана // ЖТФ. — 1995. — Т. 65(11). — C. 29-40.
[3] Bird G.A. Molecular gas dynamics. — Oxford: Clarendon Press, 1976. — 415 р.
[4] Karnidakis G., Beskok A., Aluru N. Microflows and nanoflows. — Interdisciplinary Applied Math. 29. Springer Science+Business Media, Inc., 2005. — 817 p.
[5] Thompson P. A., Robbins M.O. Shear flow near solids: Epitaxial order and flow boundary conditions // Phys. Rev. A. — 1990. — V. 41. — P. 6830-6837.
[6] Koplik J., Banavar J.R., Willemsen J.F. Molecular dynamics of fluid flow at solid surfaces // Phys. Fluids A. — 1989.—V. 1. —P. 781-794.
[7] Heinbuch U., Fischer J. Liquid flow in pores: Slip, no-slip, or multilayer sticking // Phys. Rev. A. — 1989. — V. 40. — P. 1144-1146.
[8] Гиршфельдер Дж., Кертисс Ч., Берд Р. Молекулярная теория газов и жидкостей. — М.: Иностранная литература, 1961. — 930 с.
[9] Xu L., Sedigh M.G., Sahimi M., Tsotsis T.T. Nonequilibrium molecular dynamics simulation of transport of gas mixtures in nanopores // Phys. Rev. Lett. — 1998. — V. 80. — P. 3511-3514.
[10] Рудяк В.Я., Белкин А.А., Иванов Д.А., Егоров В.В. Моделирование процессов переноса на основе метода молекулярной динамики. Коэффициент самодиффузии // ТВТ. — 2008. — Т. 46(1). — C. 35-44.
[11] Зубарев Д.Н. Неравновесная статистическая механика. — М.: Наука, 1971. — 415 с.
[12] Рудяк В.Я. Статистическая теория диссипативных процессов в газах и жидкостях. — Новосибирск: Наука, 1987. — 271 с.
[13] Lauga E., Brenner M.P., Stone H.A. Microfluidics: the no-slip boundary condition//Handbook of Exper. Fluid Dynamics. Ch. 19., N.Y.: Springer, 2007. — P. 1219-1240.
[14] Ou J., Perot B., Rothstein J.P. Laminar drag reduction in microchannels using ultrahydrophobic surfaces // Phys. Fluids. — 2004. — V. 16(12). — P. 4635-4643.
[15] Watts E.T., Krim J., Widom A. Experimental observation of interfacial slippage at the boundary of molecularly thin films with good substraite // Phys. Rev. B. — 1990. — V. 41(6). — P. 3466-3472.
[16] Soong C.Y., Yen T.H., Tzeng P.Y. Molecular dynamics simulation of nanochannels flows with effects of wall lattice-fluid interactions // Phys. Rev. E. — 2007. — V. 76. — 036303.
[17] Rudyak V., Belkin A., Egorov V., IvanovD. Modeling of plane flow with a pressure gradient in a nanochannel// Proc. of 1st European Conference on Microfluidics. Bologna, CHF, 2008, ^FLU08-38.
