Suppression of sawtooth oscillations when using a finite-difference scheme for mass transport simulation in a drying droplet on a substrate in the thin layer approximation

Full Text

Введение

Испаряющиеся капли и пленки используются в приложениях из разных областей, например, охлаждение нагретых поверхностей электронных приборов, диагностика в медицине, формирование прозрачных электропроводных покрытий на гибкой подложке, структурирование поверхности [1, 2]. Метод испарительной литографии появился после выяснения связи возникающего при испарении капель коллоидных растворов эффекта кофейных колец [3] с естественным образом формирующимися неоднородными потоками пара с поверхности капли (см. обзор [2]). В методе испарительной литографии контролируемое создание пространственных структур в осадках, остающихся на подложке после высыхания жидкости, достигается при помощи внешних условий, индуцирующих неравномерное испарение с поверхности коллоидной жидкости. Осадки могут оставаться не только на подложке, но и на стенке ячейки [4]. Испарительная литография является частью более широкого направления. Речь идет о самосборке, вызванной испарением (evaporative-induced self-assembly (EISA)). К этому обширному направлению относятся методы на основе процессов, связанных с контактной линией (граница трех фаз «жидкость–подложка–воздух»), методы на основе сил межчастичного взаимодействия и испарительная литография. Как правило, испарительная литография является гибким и одноступенчатым процессом, преимущества которого связаны с простотой, дешевизной и применимостью практически к любой подложке без предварительной обработки. В такой литографии отсутствует механическое воздействие на шаблон, поэтому его целостность в процессе работы не нарушается. Также этот метод полезен для создания материалов с локализованными функциями, такими как скользкость и самовосстановление. По этим причинам испарительная литография привлекает все большее внимание и к настоящему времени имеет ряд достижений. В [2] также проанализированы имеющиеся ограничения рассматриваемого метода и пути его дальнейшего развития. 

Методы испарительной литографии можно разделить на активные и пассивные. Их отличие в том, что первая подгруппа характеризуется наличием ключевых параметров, которые регулируются в режиме реального времени, а вторая подгруппа подразумевает наличие статических ключевых параметров, которые настраиваются до начала процесса [2]. К отдельной подгруппе относятся гибридные методы, которые также могут быть как пассивными, так и активными. Эти методы сочетают в себе испарительную литографию с другими методами, относящимися и не относящимися к испарительной самосборке [2]. Испарительная литография предоставляет больше возможностей для формирования структур различной геометрической формы на микро- и наноуровне по сравнению с испарительной самосборкой. Но, с другой стороны, испарительная самосборка позволяет получать структуры, меньшие по размеру. Разработка новых гибридных методов в испарительной литографии позволит получать структуры осадков или рельефных твердых пленок требуемой формы и морфологии еще меньших размеров для относительно больших площадей. Такие покрытия будут устойчивы к внешним механическим воздействиям, в них будут отсутствовать трещины [5]. Кроме того, эти покрытия можно наделить некоторыми требуемыми функциональными свойствами.

К примеру, для разработки гибридных методов можно использовать дополнительное воздействие через различные факторы: пропускание электрического тока через подложку с испаряющейся каплей [6], воздействие магнитным полем на частицы [7], локальный нагрев подложки [8], направление вектора силы тяжести относительно расположения капли и влажность окружающего воздуха [9, 10]. Напряжение электрического тока в подложке влияет на ее смачиваемость жидкостью и на режим контактной линии: пиннинг (закрепление границы) или режим постоянного краевого угла (скольжение границы) [6]. Эти факторы влияют на геометрию капли, пространственную неоднородность испарения вдоль ее свободной поверхности и поле скорости потока жидкости. Такой способ управления можно использовать в режиме реального времени, получая при этом требуемую форму осадка. Еще один дополнительный способ контроля заключается в создании магнитного поля 
в области капли [7]. В эксперименте [7] в отсутствие магнитного поля формировался кольцевой осадок частиц. Увеличение силы магнитного поля приводило к увеличению количества частиц, осаждающихся в центральной области в виде пятна. Сидячие и висячие капли на подложке изучались в эксперименте при разной влажности воздуха [9]. Направление вектора силы тяжести относительно расположения капли влияет на то, будут ли поток Марангони и объемная тепловая конвекция сонаправлены или нет. В сидячей капле эти потоки противодействуют друг другу, а в висячей капле, наоборот, усиливают друг друга. Также не стоит забывать о капиллярном потоке. Влажность воздуха влияет на скорость испарения. Комбинации этих параметров приводили к возникновению разных структур осадков: кольцо, равномерное пятно, диск или центральное пятно [9]. Коллоидную литографию возможно комбинировать с испарительной литографией [11]. Для получения масок из микрочастиц на подложке можно использовать испарительную литографию, а затем на их основе с помощью коллоидной литографии получать упорядоченные осадки наночастиц. Эксперимент с сохнущими каплями [11] показал, что равномерная упорядоченная морфология осадка лучше получается на гидрофильных подложках, чем на гидрофобных. Это связано с отличием поведения контактной линии, что зависит от смачиваемости поверхности.

Еще одно важное направление заключается в использовании смесей частиц. Эти смеси могут состоять из частиц разного размера [12, 13], формы [14], материала и т.д. Частицы Януса состоят из двух частей, материалы которых отличаются по физико-химическим свойствам [15]. Разделение частиц по размеру вблизи трехфазной границы высыхающей капли моделировалось в [13]. В бинарных смесях частиц разного размера могут возникать силы исключенного объема (энтропийные силы), которые способствуют притяжению крупных частиц друг к другу [16, 17]. Химическое воздействие на частицы с помощью поверхностно-активных веществ также позволяет контролировать форму осадка и его морфологию [18].

Испарение коллоидной жидкости из ячейки Хелли–Шоу также можно отнести к испарительной литографии [19, 20]. Такая ячейка состоит из двух параллельно расположенных пластин, между которыми есть узкая щель, в которой заключен коллоидный раствор. Все боковые стороны ячейки, кроме одной, закрыты. Жидкость испаряется через открытое боковое отверстие, поэтому этот процесс называют направленным испарением. Частицы переносятся капиллярным потоком в сторону направленного испарения. Возможно образование как сплошного осадка [19], так и периодических полос [20]. В [20] изучались такие параметры, как степень адсорбции частиц к подложке и скорость испарения, влияющие на морфологию осадка. Сложная структура возникающих потоков жидкости была показана в эксперименте [19].

На примере капли солевого раствора показана возможность управлять формированием кристаллического осадка с помощью точечного лазерного нагрева локального участка свободной поверхности жидкости [21]. В эксперименте изучалось влияние таких параметров, как мощность лазера и смачиваемость/несмачиваемость подложки. Комбинации значений этих параметров приводили к таким кристаллическим паттернам, как кольцо, спираль, пятно и прочее. Визуализация структуры потока показала, что течение направлено в сторону центральной зоны нагрева от подложки к границе жидкости и воздуха. Вдоль свободной поверхности поток направлен от центра 
к периферии капли. На перенос растворенного вещества в этом эксперименте в большей степени влиял тепловой поток Марангони и капиллярный поток [21]. Во-первых, неравномерный нагрев поверхности приводит к неравномерной локальной плотности потока пара. Во-вторых, перепад температуры влияет на возникновение градиента поверхностного натяжения вдоль свободной поверхности жидкости. Эти два фактора объясняют наблюдаемую структуру потока жидкости. Концентрация соли растет в области нагрева поверхности не только за счет переноса потоком, но и благодаря интенсивному испарению, происходящему в зоне воздействия лазера [21]. Для сравнения в эксперименте [22] при испарении капли солевого раствора без какого-либо внешнего воздействия преобладал концентрационный поток Марангони. С помощью метода PIV (Particle Image Velocimetry) была изучена динамика структуры потока. Исследование показало, что зарождение кристаллов и их рост в процессе испарения жидкости приводит к нарушению осесимметричного потока. Вокруг кристаллов возникают симметричные вихри [22].

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

Моделирование массопереноса в высыхающих каплях важно, так как численные результаты позволяют подобрать необходимые параметры для проведения экспериментальных исследований. Это дает возможность совершенствовать существующие методы и разрабатывать новые приложения. Континуальные модели позволяют описать форму осадка, но не его морфологию [26–28]. Полудискретные модели, в которых частицы описываются точками, тоже не в состоянии предсказывать морфологию осадка [29–32]. Решеточные модели описывают лишь частицы в форме кубоидов [33–37]. Методы молекулярной динамики и диссипативной динамики частиц позволяют делать прогнозы лишь для относительно малого числа частиц [38–40]. Безрешеточные модели на основе метода Монте– Карло, в которых явно отслеживается динамика каждой частицы, лишены упомянутых недостатков [13, 41–45]. Но эти модели либо не учитывают гидродинамику, либо используют простые аналитические решения для частного случая. Для улучшения этих моделей, алгоритмов и программ необходимо разрабатывать дополнительные модули, позволяющие учесть различные эффекты, влияющие на формирование осадков, в том числе и гидродинамику. В коммерческом пакете Comsol Multiphysics есть модули, позволяющие рассчитывать динамику частиц и гидродинамику [46], но цена на этот пакет относительно высокая. Это касается и дополнительных тулбоксов к Matlab [47]. Недавно для нашей страны был ограничен доступ к официальному сайту свободного пакета для молекулярной динамики LAMMPS (https://www.lammps.org/). Зарубежные коммерческие пакеты в любой момент могут оказаться недоступными из-за политической ситуации. К примеру, на момент написания данной статьи приостановлены продажи пакета Ansys в РФ. В условиях импортозамещения очень важной является разработка отечественного программного обеспечения для моделирования, написание своих кодов и библиотек. Для проведения исследований в области испарительной литографии требуется разработка программного комплекса, включающего ряд модулей (рис. 1). Пустой прямоугольник на схеме обозначает дополнительные модули, которые необходимо создать. К примеру, в случае необходимости учета таких важных факторов, как взаимодействие частиц с подложкой (адсорбция или адгезия) или со свободной поверхностью капли (капиллярные силы).

Разработка модулей для моделирования диффузии и капиллярного притяжения сферических частиц обсуждалась ранее в [42, 44, 45]. Также был рассмотрен случай для бинарной смеси частиц разного размера [13]. Цель текущей работы — разработать численный алгоритм и программный модуль для расчета гидродинамики в высыхающей капле для дальнейшего использования в программном комплексе, ориентированном на решение задач в области испарительной литографии. Но при разработке необходимо учитывать одно требование. Алгоритм должен быть рассчитан на большое количество предельно малых временных шагов (минимальные вычислительные затраты на каждом временном шаге), что позволит корректно работать уже имеющемуся модулю «Диффузия частиц», основанному на методе Монте– Карло, для явного предсказания динамики частиц, так как в этом случае выполняется соотношение Эйнштейна–Смоулховского [45]. Величина временного шага $\Delta t$ косвенно зависит от размера частиц, например, для частиц с радиусом 0.35 мкм в [45] использовался шаг $\Delta t= 10^{-4}$ с.

Рис. 1. Схема планируемого программного комплекса для решения задач в области испарительной литографии
[Figure 1. Scheme of the planned software complex for solving problems in the field of evaporative lithography]

Также здесь стоит отметить, по какой причине не подходят существующие готовые решения. Описанное в [48] аналитическое решение базируется на смеси кинематического подхода и приближения смазки, которое еще известно как приближение тонкого слоя. Под кинематическим подходом здесь понимается нахождение усредненной радиальной скорости потока из закона сохранения массы. Такой подход не объясняет природу возникновения потока, как например, капиллярные силы, когда в задаче явным образом учитывают градиент давления Лапласа и кривизну свободной поверхности. Как видно из полученного решения [48], оно явным образом не зависит от ряда физических параметров, как, например, вязкость жидкости. Кроме того, это решение опирается на частный случай с определенным видом функции плотности потока пара $J(r, t)$, что не позволяет решать задачи, связанные с испарительной литографией. Здесь $r$ — радиальная координата и $t$ — время. Различные маски, излучатели и прочее внешнее воздействие оказывают влияние на концентрацию пара вблизи свободной поверхности, что в итоге влияет на капиллярный поток жидкости [2]. Таким образом, для определения $J(r, t)$ необходимо численно моделировать перенос пара в воздухе. Предложенная в [49] неявная разностная схема является неустойчивой. Вблизи границ появляются осцилляции, которые со временем разрастаются на всю область. Авторы [49] подробно не разъясняют способ численной реализации граничных условий, ссылаясь на приложенный код в дополнении. Но восстановить некоторые граничные условия из кода крайне сложно. В программе [49] на каждом временном шаге выполняется искусственное сглаживание осциллирующих точек с помощью специального фильтра. Кроме того, на каждом шаге необходимо выполнять такие вычислительно затратные операции, как нахождение обратной матрицы и умножение матрицы на вектор. По этим причинам данный программный модуль [49] нам не подходит. В отличие от [27, 28, 49] здесь дополнительно будет описана численная реализация расчета поля скорости потока жидкости, а не только усредненной по толщине жидкого слоя радиальной компоненты скорости.

1. Физическая постановка задачи

Капля жидкости испаряется на твердом горизонтальном непромокаемом основании при нормальных комнатных условиях. Считаем, что жидкость несжимаема. Трехфазная граница закреплена, поэтому радиус контакта капли и подложки $R$ — постоянная величина. Рассматриваем случай относительно малого объема жидкости, когда сила поверхностного натяжения преобладает над силой тяжести, число Бонда $\mathsf{Bo} = \rho h_0^2 g/ \sigma \approx 1.4\cdot 10^{-3} \ll 1$, где $g$ — ускорение свободного падения, $\rho$ — плотность жидкости, $\sigma$ — коэффициент поверхностного натяжения и $h_0=h(0, 0)$ — начальная высота капли. Здесь функция $h =h(r,t)$ задает профиль свободной поверхности жидкости (граница «жидкость – воздух»). В этом случае форма капли близка форме сферического сегмента. По причине осевой симметрии удобно использовать цилиндрическую систему координат $(r, z)$. Движение жидкости вызывает локальные изменения кривизны поверхности в процессе испарения жидкости, приводящие к возникновению градиента давления Лапласа (капиллярный поток). В случае малых значений краевого угла смачивания $\theta$ капиллярный поток преобладает над потоком Марангони (см. оценку в [44]). Здесь речь идет о тонких каплях, для которых аспектное отношение $\varepsilon = h_0 / R \ll 1$. Значения параметров задачи приведены в табл. 1.

Таблица 1. Физические и геометрические параметры задачи [Physical and geometric parameters of the problem]

Обозначение [Designation]	Расшифровка [Description]	Величина [Value]	Единицы измерения [Units of Measurement]
$h_0$	высота капли [drop height]	$10^{-4}$	m
$R$	радиус основания капли [drop base radius]	$10^{-3}$	m
$h_f$	толщина жидкого слоя в районе контактной линии [thickness of the liquid layer in the area of the contact line]	$0.01 h_0$	m
$\varepsilon= h_0 / R$	аспектное отношение [aspect ratio]	0.1	–
$\theta \approx 2 \varepsilon$	краевой угол смачивания [contact angle]	0.2	radians
$C_0$	начальная массовая доля [initial mass fraction]	0.05	–
$C_g$	концентрация фазового перехода, например, гелеобразования [phase transition concentration, e.g. gelation]	0.7	–
$D$	коэффициент диффузии растворенного или взвешенного вещества [diffusion coefficient of a dissolved or suspended matter]	$10^{-10}$	m$^2$/s
$\eta_0$	вязкость чистой жидкости [viscosity of a pure liquid]	$10^{-3}$	Pa${}\cdot{}$s
$\rho$	плотность жидкости [liquid density]	$10^3$	kg/m$^3$
$\sigma$	коэффициент поверхностного натяжения [surface tension coefficient]	0.072	N/m
$t_f$	время полного испарения [total evaporation time]	450	s
$D_v$	коэффициент диффузии пара [vapor diffusion coefficient]	$2.4\cdot 10^{-5}$	m$^2$/s
$C_v$	плотность насыщенного пара [saturated steam density]	$2.32\cdot 10^{-2}$	kg/m$^3$
$H$	относительная влажность [relative humidity]	0.4	–

Распределение растворенного или взвешенного в жидкости вещества будем описывать функцией $\langle C \rangle(r,t)$. Предполагаем, что массовая доля вещества $\langle C \rangle$ не зависит от координаты $z$. Если время диффузионной релаксации вещества много меньше времени полного испарения, т.е. $t_d \ll t_f$, то устанавливается равномерная концентрация вещества вдоль вертикального направления. Здесь время $t_d= h_0^2 / D =$ 100 с и время $t_f =$ 450 с. Таким образом, величина $\langle C \rangle$ является усредненной по высоте жидкого слоя.

2. Математическая модель

Запишем уравнение неразрывности и стационарные уравнения Навье– Стокса без учета инерционных слагаемых в цилиндрической системе координат:
\[ \begin{equation}
 \frac{1}{r} \frac{\partial (r u)}{\partial r} + \frac{\partial w}{\partial z} = 0,
\end{equation} \tag{1} \]
\[ \begin{equation}
 \frac{\partial}{\partial r}\Bigl(\frac{\eta}{r}\frac{\partial (r u)}{\partial r} \Bigr) + 
 \frac{\partial}{\partial z}\Bigl( \eta \frac{\partial u}{\partial z} \Bigr)= 
 \frac{\partial P}{\partial r},
\end{equation} \tag{2} \]
\[ \begin{equation}
 \frac{1}{r} \frac{\partial}{\partial r}\Bigl( \eta r \frac{\partial w}{\partial r} \Bigr) + 
 \frac{\partial}{\partial z}\Bigl( \eta \frac{\partial w}{\partial z}\Bigr) = 
 \frac{\partial P}{\partial z},
\end{equation} \tag{3} \]
где $\eta = \eta(\langle C \rangle)$ — динамическая вязкость жидкости. Давление $P$, горизонтальная и вертикальная компоненты вектора скорости $\mathbf{v}=(u, w)$ являются функциями, зависящими от $t$, $r$ и $z$. Система уравнений (1)–(3) справедлива для несжимаемой жидкости и малых значений числа Рейнольдса $\mathsf{Re}=\rho u_c h_0/ \eta_0$. Как правило, характерная скорость потока жидкости $u_c \approx$ 1 мкм/с в испаряющейся при комнатных условиях капле воды. В таком случае число Рейнольдса $\mathsf{Re} \approx 10^{-4} \ll 1$.

Для перехода к безразмерным записям рассмотрим масштабные соотношения: $\eta = \eta_0 \tilde \eta$, $J = \tilde J J_c$, $u= u_c \tilde u$, $w= \varepsilon u_c \tilde w$, $t= t_c \tilde t$, $P= P_c \tilde P$, $r= R \tilde r$ и ${z= h_0 \tilde z}$. Здесь знаком «тильда» помечаем безразмерные величины. Заметим, что горизонтальные и вертикальные составляющие (размер, скорость) масштабируются по-разному. Такой подход лежит в основе приближения смазки. 

Зададим характерные величины: скорость $u_c= \eta_0/ (\rho h_0)\approx 0.01$ м/с, время $t_c= R/ u_c\approx 0.1$ с, плотность потока пара $J_c = \varepsilon u_c \rho\approx 1$ кг/(м$^2{}\cdot{}$c) и давление $P_c = R u_c \eta_0/ h_0^2 \approx 1$ Па. Запишем уравнение неразрывности (1) в безразмерном виде:
\[ \begin{equation}
 \frac{1}{\tilde r} \frac{\partial (\tilde r \tilde u)}{\partial \tilde r} + \frac{\partial \tilde w}{\partial \tilde z} = 0.
\end{equation} \tag{4} \]
Заметим, что уравнение (4) не отличается по форме записи от (1). Уравнения (2) и (3) записываются так:
\[ \begin{equation}
\varepsilon^2 \frac{\partial}{\partial \tilde r}\Bigl( \frac{\tilde \eta}{\tilde r}\frac{\partial (\tilde r \tilde u)}{\partial \tilde r} \Bigr) + 
\frac{\partial}{\partial \tilde z}\Bigl( \tilde \eta \frac{\partial \tilde u}{\partial \tilde z}\Bigr) =
\frac{\partial \tilde P}{\partial \tilde r},
\end{equation} \tag{5} \]
\[ \begin{equation}
\varepsilon^3 \frac{1}{\tilde r} \frac{\partial}{\partial \tilde r}\Bigl( \tilde \eta \tilde r \frac{\partial \tilde w}{\partial \tilde r} \Bigr) + \varepsilon^2 \frac{\partial}{\partial \tilde z} \Bigl( \tilde \eta \frac{\partial \tilde w}{\partial \tilde z}\Bigr) = \frac{\partial \tilde P}{\partial \tilde z}.
\end{equation} \tag{6} \]

Пренебрегая слагаемыми, при которых стоят $\varepsilon^2$ и $\varepsilon^3$, в (5) и (6), получаем
\[ \begin{equation}
 \frac{\partial \tilde P}{\partial \tilde r} - \frac{\partial}{\partial \tilde z}\Bigl( \tilde \eta \frac{\partial \tilde u}{\partial \tilde z}\Bigr) = 0,
\end{equation} \tag{7} \]
\[ \begin{equation}
 \frac{\partial \tilde P}{\partial \tilde z} = 0.
\end{equation} \tag{8} \]
В результате получена упрощенная система уравнений (4), (7) и (8).

Из (8) следует, что давление $\tilde P$ не зависит от $\tilde z$, тогда вместо $\tilde P$ будем рассматривать усредненную по высоте жидкого слоя величину $\langle \tilde P \rangle = \langle \tilde P \rangle (\tilde r, \tilde t)$. Интегрируя левую и правую части (7) и учитывая, что вязкость $\tilde \eta$ здесь не зависит явно или косвенно от $\tilde z$,
\[ \begin{equation*}
\frac{\partial \langle \tilde P \rangle}{\partial \tilde r}\int d\tilde z - \tilde \eta \int \frac{\partial^2 \tilde u}{\partial \tilde z^2}\,d\tilde z = \int 0\,d\tilde z, 
\end{equation*} \]
получаем
\[ \begin{equation*}
\frac{\partial \langle \tilde P \rangle}{\partial \tilde r} (\tilde z + C_1) - 
\tilde \eta \Bigl( \frac{\partial \tilde u}{\partial \tilde z} + C_2 \Bigr) = C_3.
\end{equation*} \]

Как было отмечено ранее, термокапиллярные потоки здесь не рассматриваются, поэтому на свободной границе $\tilde z = \tilde h$ выполняется баланс касательных напряжений $\partial \tilde u/ \partial \tilde z =0$, тогда
\[ \begin{equation*}
\frac{\partial \langle \tilde P \rangle}{\partial \tilde r} \tilde h + \frac{\partial \langle \tilde P \rangle}{\partial \tilde r} C_1 = C_3 + \tilde \eta C_2.
\end{equation*} \]

Для $\tilde r = 0$ выполняется условие осевой симметрии $\partial \langle \tilde P \rangle/ \partial \tilde r = 0$, тогда $C_3  = - \tilde \eta(0, \tilde t) C_2$. В действительности вязкость $\tilde \eta$ зависит от $\tilde r$ и $\tilde t$ неявным образом, через массовую долю $\langle C \rangle$, но здесь для краткости пишем $\tilde \eta(\tilde r, \tilde t)$. На границе $\tilde r = \tilde R$ горизонтальная компонента скорости является постоянной величиной, не зависящей от $\tilde z$ (в силу прилипания к подложке $\tilde u = 0$). Тогда в этой точке $\partial^2 \tilde u/ \partial \tilde z^2 = 0$. Учитывая это, из (7) получаем граничное условие $\partial \langle \tilde P \rangle / \partial \tilde r = 0$ для $\tilde r = \tilde R$. С таким условием приходим к соотношению $C_3  = - \tilde \eta(\tilde R, \tilde t) C_2$. На разных границах получили одно и то же соотношение, хотя в общем случае $\tilde \eta(0, \tilde t) \neq \tilde \eta(\tilde R, \tilde t)$. Приходим к выводу, что $C_2 = C_3 = 0$. В таком случае константа интегрирования $C_1 = - \tilde h$ и последнее выражение принимает вид
\[ \begin{equation*}
\frac{\partial \langle \tilde P \rangle}{\partial \tilde r} (\tilde z - \tilde h) - \tilde \eta \frac{\partial \tilde u}{\partial \tilde z} = 0.
\end{equation*} \]

Интегрируем полученное еще раз:
\[ \begin{equation*}
\frac{\partial \langle \tilde P \rangle}{\partial \tilde r} \int \tilde z\,d\tilde z - \frac{\partial \langle \tilde P \rangle}{\partial \tilde r} \tilde h \int d\tilde z - \tilde \eta \int \frac{\partial \tilde u}{\partial \tilde z}\,d\tilde z = \int 0\,d\tilde z,
\end{equation*} \]
получаем
\[ \begin{equation*}
\frac{\partial \langle \tilde P \rangle}{\partial \tilde r}\Bigl( \frac{\tilde z^2}{2} + C_4 - \tilde h \tilde z +C_5 \Bigr) - \tilde \eta (\tilde u + C_6) = C_7.
\end{equation*} \]

Учитывая, что для $\tilde z = 0$ выполняется условие прилипания $\tilde u = 0$ и для $\tilde r = 0$ выполняется условие $\partial \langle \tilde P \rangle/ \partial \tilde r = 0$, получаем $C_7 = - \tilde \eta(0, \tilde t) C_6$. Учитывая аналогичные граничные условия для $\tilde r = \tilde R$, получаем $C_7 = - \tilde \eta(\tilde R, \tilde t) C_6$. Приходим к выводу, что $C_6 = C_7 = 0$. Затем, учитывая лишь условие $\tilde u = 0$ для $\tilde z = 0$, получаем $C_4 = - C_5$. В итоге выражаем горизонтальную компоненту вектора скорости:
\[ \begin{equation}
\tilde u = \frac{H_a(x)}{\tilde \eta} \frac{\partial \langle \tilde P \rangle}{\partial \tilde r}\Bigl( \frac{\tilde z^2}{2} - \tilde h \tilde z \Bigr).
\end{equation} \tag{9} \]

В некоторых уравнениях модели, в том числе и в (9), дополнительно используется аналитическое приближение функции Хевисайда $H_a$, чтобы исключить некоторые виды массопереноса при возникновении фазового перехода, как, например, «золь–гель», «жидкость–стекло» и т.п.: 
\[ \begin{equation*}
H_a(x) = \frac{1}{1+\exp(-2 k_h x)},
\end{equation*} \]
где $x = C_g - \langle C \rangle - \delta x$, $k_h = 10/\delta x$, $\delta x$ — ширина области перехода. Теперь, подставляя (9) в (4) и интегрируя полученное:
\[ \begin{equation*}
\int \frac{1}{\tilde r} \frac{\partial}{\partial \tilde r} \Bigl( \frac{\tilde r H_a(x)}{\tilde \eta} \frac{\partial \langle \tilde P \rangle}{\partial \tilde r} \Bigl( \frac{\tilde z^2}{2} - \tilde h \tilde z \Bigr) \Bigr)\,d \tilde z + \int \frac{\partial \tilde w}{\partial \tilde z}\,d \tilde z = \int 0\,d \tilde z,
\end{equation*} \]
получаем
\[ \begin{equation*}
\frac{1}{\tilde r} \frac{\partial}{\partial \tilde r} \Bigl( \frac{\tilde r H_a(x)}{\tilde \eta} \frac{\partial \langle \tilde P \rangle}{\partial \tilde r} \Bigl( \frac{\tilde z^3}{6} - \frac{\tilde h \tilde z^2}{2} + C_8 \Bigr) \Bigr) + \tilde w + C_9 = 0.
\end{equation*} \]

Запишем полученное выражение в развернутом виде:
\[ \begin{multline*}
\frac{H_a(x)}{\tilde r \tilde \eta} \frac{\partial \langle \tilde P \rangle}{\partial \tilde r} 
\Bigl( \frac{\tilde z^3}{6} - \frac{\tilde h \tilde z^2}{2} \Bigr) + 
\frac{H_a(x)}{\tilde \eta} \Bigl( \frac{\tilde z^3}{6} - \frac{\tilde h \tilde z^2}{2} \Bigr) 
\frac{\partial^2 \langle \tilde P \rangle}{\partial \tilde r^2} - \frac{H_a(x) \tilde z^2}{2 \tilde \eta} \frac{\partial \langle \tilde P \rangle}{\partial \tilde r} \frac{\partial \tilde h}{\partial \tilde r} + 
{}
\\
{}+ \frac{C_8 H_a(x)}{\tilde r \tilde \eta} \frac{\partial \langle \tilde P \rangle}{\partial \tilde r} + \frac{C_8 H_a(x)}{\tilde \eta} \frac{\partial^2 \langle \tilde P \rangle}{\partial \tilde r^2} + \frac{\partial H_a(x)}{\partial \tilde r} \frac{1}{\tilde \eta} \frac{\partial \langle \tilde P \rangle}{\partial \tilde r} \Bigl( \frac{\tilde z^3}{6} - \frac{\tilde h \tilde z^2}{2} + C_8 \Bigr) + \tilde w + C_9 = 0.
\end{multline*} \]
Учитывая условие непротекания $\tilde w = 0$ для $\tilde z = 0$ и условие $\partial \langle \tilde P \rangle/ \partial \tilde r = 0$ для $\tilde r = 0$ и $\tilde r = \tilde R$, получаем два соотношения:
\[ \begin{equation*}
\frac{C_8 H_a(x)}{\tilde \eta(0, \tilde t)} \frac{\partial^2 \langle \tilde P \rangle(0, \tilde t)}{\partial \tilde r^2} + C_9 = 0 \text{ и } 
\frac{C_8 H_a(x)}{\tilde \eta(\tilde R, \tilde t)} \frac{\partial^2 \langle \tilde P \rangle(\tilde R, \tilde t)}{\partial \tilde r^2} + C_9 = 0.
\end{equation*} \]
Из чего следует, что $C_8 = C_9 =0$. В таком случае получаем выражение для вертикальной компоненты вектора скорости:
\[ \begin{equation}
\tilde w = -\frac{1}{\tilde r} \frac{\partial}{\partial \tilde r}\Bigl( \frac{\tilde r H_a(x)}{\tilde \eta} \frac{\partial \langle \tilde P \rangle}{\partial \tilde r} \Bigl( \frac{\tilde z^3}{6} - \frac{\tilde h \tilde z^2}{2} \Bigr) \Bigr).
\end{equation} \tag{10} \]

Движение свободной поверхности капли в процессе испарения описывается законом сохранения (см. ссылки в [28])
\[ \begin{equation}
\frac{\partial h}{\partial t} + \frac{1}{r} \frac{\partial (rh\langle u\rangle)}{\partial r} = -\frac{J}{\rho} \sqrt{1+\Bigl(\frac{\partial h}{\partial r}\Bigr)^2},
\end{equation} \tag{11} \]
где $\langle u\rangle$ — усредненная по высоте жидкого слоя скорость радиального потока жидкости. В безразмерном виде (11) записывается так:
\[ \begin{equation}
\frac{\partial \tilde h}{\partial \tilde t} + \frac{1}{\tilde r} \frac{\partial (\tilde r \tilde h\langle \tilde u\rangle)}{\partial \tilde r} = - \tilde J,
\end{equation} \tag{12} \]
где пренебрегаем $\sqrt{1+\varepsilon^2 (\partial \tilde h/ \partial \tilde r)^2}$ в правой части как величиной второго порядка малости. С учетом (9) для усредненной радиальной скорости $\langle \tilde u\rangle$ имеем
\[ \begin{equation}
\langle \tilde u\rangle = \frac{1}{\tilde h} \int _{0}^{\tilde h} \tilde u\, d\tilde z = -H_a(x)\frac{\tilde h^2}{3 \tilde \eta} \frac{\partial \langle \tilde P \rangle}{\partial \tilde r}.
\end{equation} \tag{13} \]
Капиллярное давление зависит от локальной кривизны поверхности:
\[ \begin{equation}
 \langle \tilde P \rangle = - \frac{1}{\mathsf{Ca}} \frac{1}{\tilde r} \frac{\partial}{\partial \tilde r} 
 \Bigl( \tilde r \frac{\partial \tilde h}{\partial \tilde r} \Bigr),
\end{equation} \tag{14} \]
что следует из формулы для давления Лапласа [28] с учетом приближения смазки. Здесь капиллярное число $\mathsf{Ca}= \eta_0 u_c/ (\sigma \varepsilon^3) = \eta^2 /(\varepsilon^3 \rho h_0 \sigma) \approx 0.14$. В [49] выражения (13) и (14) подставляются в уравнение (12) и затем строится разностная схема для уравнения с производной четвертого порядка по $r$ относительно $h$. Чтобы избежать таких громоздких разностных схем, здесь (13) и (14) рассматриваются как отдельные уравнения.

Для описания переноса растворенного или взвешенного вещества используем уравнение конвекции–диффузии [28]
\[ \begin{equation}
\frac{\partial \langle C \rangle}{\partial t} + \langle u\rangle \frac{\partial \langle C \rangle}{\partial r} =
\frac{D}{r h} \frac{\partial}{\partial r} 
\Bigl( r h \frac{\partial \langle C \rangle}{\partial r} \Bigr) + 
\frac{J \langle C \rangle}{\rho h} \sqrt{1+\Bigl(\frac{\partial h}{\partial r}\Bigr)^2},
\end{equation} \tag{15} \]
где $D$ — коэффициент диффузии растворенного или взвешенного вещества. Подробный вывод уравнения (15) описан в [50]. В безразмерном виде уравнение (15) с учетом приближения смазки принимает вид
\[ \begin{equation}
\frac{\partial \langle C \rangle}{\partial \tilde t} + \langle \tilde u\rangle 
\frac{\partial \langle C \rangle}{\partial \tilde r} = \frac{H_a(x)}{\mathsf{Pe}} \frac{1}{\tilde r \tilde h} \frac{\partial}{\partial \tilde r} 
\Bigl( \tilde r \tilde h \frac{\partial \langle C \rangle}{\partial \tilde r} \Bigr) + 
\frac{\tilde J \langle C \rangle}{\tilde h},
\end{equation} \tag{16} \]
где $\mathsf{Pe}$ — число Пекле, $\mathsf{Pe}\approx R u_c/ D\approx 10^5$. В итоге система уравнений включает в себя уравнения (9), (10), (12), (13), (14) и (16).

3. Вязкость и плотность потока пара

Теперь необходимо ввести замыкающие соотношения для некоторых величин. В предложенной здесь модели вязкость зависит от массовой доли вещества. Для описания этой зависимости используем формулу Муни
\[ \begin{equation*}
\tilde \eta = \exp \Bigl(\frac{S \langle C \rangle}{1 - K \langle C \rangle}\Bigr),
\end{equation*} \]
где значения эмпирических параметров $S = 1.692$ и $K = 1.236$ возьмем из [51]. Также считаем, что плотность потока пара зависит от концентрации вещества и толщины жидкого слоя:
\[ \begin{equation*}
\tilde J = \tilde J_0 \frac{ 1 - \langle C \rangle^2 / C_g^2}{\kappa + \tilde h},
\end{equation*} \]
где $C_g$ — критическая массовая доля, при которой происходит фазовый переход [27] (золь–гель, стеклообразование, кристаллизация и т. п.), полагаем $C_g \approx 0.7$ [51]. Для расчетов будем использовать значение параметра $\kappa = 1$. Это некоторая аппроксимация для учета испарения, поэтому в дальнейшем планируется разработать модуль численного расчета переноса пара в воздухе (рис. 1). 
Величина $\tilde J_0$ рассчитывается так:
\[ \begin{equation*}
\tilde J_0 = \frac{D_v C_v (1 - H)}{J_c R} (0.27 \theta^2 + 1.3)(0.6381 - 0.2239(\theta - \pi / 4)^2),
\end{equation*} \]
где $H$ — относительная влажность, $D_v$ — коэффициент диффузии пара и $C_v$ — плотность насыщенного пара [48].

4. Начальные и граничные условия

Запишем начальные условия задачи. Для тонких капель, размер которых меньше капиллярной длины, допустимо использовать параболическое приближение начальной формы свободной поверхности $\tilde h(\tilde r, 0)= 1 - \tilde r^2$ [48]. Начальную массовую долю вещества опишем с помощью выражения [27]
\[ \begin{equation*}
\langle C \rangle (\tilde r, 0) = C_g \frac{2 - C_0 + 2(C_0 -1)}{1 + \exp(w_c (\tilde r - 1))},
\end{equation*} \]
где параметр $w_c$ регулирует начальную ширину кольцевого осадка (здесь в расчетах используем $w_c = 30$). Такое условие предполагает, что вещество в начале процесса распределено практически равномерно вдоль радиальной координаты, за исключением района периферии ($\tilde r \approx \tilde R = 1$), где концентрация устремляется к пороговому значению $C_g$, что обеспечивает закрепление трехфазной границы «жидкость–подложка–воздух».

Для $\tilde r=0$ по причине осевой симметрии задаются следующие граничные условия:
\[ \begin{equation*}
\frac{\partial \langle \tilde P \rangle}{\partial \tilde r}
\Big|_{\tilde r=0}
= \frac{\partial \tilde h}{\partial \tilde r}\Big|_{\tilde r=0} 
= \frac{\partial \langle C \rangle}{\partial \tilde r}\Big|_{\tilde r=0}
= \frac{\partial \tilde w}{\partial \tilde r}\Big|_{\tilde r=0} = \tilde u\big|_{\tilde r=0} = \langle \tilde u \rangle\big|_{\tilde r=0} = 0.
\end{equation*} \]
На границе $\tilde r= \tilde R$ выполняются условия
\[ \begin{equation*}
\frac{\partial \langle \tilde P \rangle}{\partial \tilde r}\Big|_{\tilde r= \tilde R} = 
\tilde w\big|_{\tilde r= \tilde R} = \tilde u\big|_{\tilde r= \tilde R} = \langle \tilde u \rangle\big|_{\tilde r= \tilde R} = 0, \quad 
\tilde h \big|_{\tilde r= \tilde R} = \tilde h_f,\quad 
 \langle C \rangle\big|_{\tilde r= \tilde R} = C_g.
\end{equation*} \]
Условие для $\langle \tilde P \rangle$ уже обсуждалось в п. 2. Вертикальная и горизонтальная компоненты вектора скорости потока жидкости записываются из соображений непротекания и прилипания. Здесь $\tilde h_f$ — толщина жидкого слоя в районе контактной линии.

Граничные условия вдоль подложки следующие: 
\[ \begin{equation*}
\tilde w\big|_{\tilde z = 0} = \tilde u \big|_{\tilde z = 0} = 0,
\end{equation*} \]
по причине непротекания и прилипания. Вдоль свободной границы жидкости условия принимают вид 
\[ \begin{equation*}
\frac{\partial \tilde u}{\partial \tilde z}\Big|_{\tilde z = \tilde h} = \tilde w\big|_{\tilde z = \tilde h} = 0
\end{equation*} \]
в предположении, что на этой границе нет вязкого трения и отсутствует поток жидкости через границу «жидкость–воздух» (из-за медленного испарения влияние подвижной границы и отдачи пара не учитываем). В приближении смазки на этой границе $\tilde w$ и $\tilde u$ соответствуют нормальной и касательной компонентам вектора скорости. Последние условия для скорости потока на свободной поверхности капли являются излишними, так как значения $\tilde u$ и $\tilde w$ на этой границе несложно рассчитать с использованием выражений (9) и (10).

5. Численный метод решения задачи.

В этом пункте рассматриваются только безразмерные величины, поэтому знак «тильда» опускаем для краткости. Для дискретной аппроксимации уравнений модели с помощью разностных схем введем пространственно-временную сетку с координатами $r_n = n \Delta r$, $z_m = m \Delta z$ и $t_k = k \Delta t$, где $n= 0, 1, \dots , N$; $m= 0, 1, \dots , M$; $k= 0, 1, \dots , K$. Параметры $\Delta t$, $\Delta r$ и $\Delta z$ обозначают шаги по времени и пространству соответственно. Кроме того, сетка характеризуется такими величинами, как количество отрезков ($N$, $M$ и $K$) и узлов ($N+1$, $M+1$ и $K+1$). Параметр $K = [t_f / \Delta t]$, здесь квадратные скобки обозначают целочисленное деление. Шаги по пространству определяются как $\Delta r= R/N$ и $\Delta z= h_0/M$.

Перепишем уравнения (9), (10), (12), (13), (14) и (16) в дискретном виде:
\[ \begin{equation}
u_{n,m}^k = \frac{H_a(x_n)}{\eta_n^k} \frac{\langle P \rangle_{n+1}^k - \langle P \rangle_{n-1}^k}{2\, \Delta r}\Bigl( \frac{z_m^2}{2} - h_n^k z_m \Bigr),
\end{equation} \tag{17} \]
\[ \begin{equation}
w_{n,m}^k = -\frac{H_a(x_{n+1/2})}{r_n\, \Delta r} 
\frac{r_{n+1/2}}{\eta_{n+1/2}^k} \frac{ \langle P \rangle_{n+1}^k - \langle P \rangle_n^k}{\Delta r} 
\Bigl( \frac{z_m^3}{6} - \frac{h_{n+1/2}^k z_m^2}{2} \Bigr)+\frac{H_a(x_{n-1/2})}{r_n\, \Delta r} \frac{r_{n-1/2}}{\eta_{n-1/2}^k} \frac{ \langle P \rangle_n^k - \langle P \rangle_{n-1}^k}{\Delta r} \Bigl( \frac{z_m^3}{6} - \frac{h_{n-1/2}^k z_m^2}{2} \Bigr),
\end{equation} \tag{18} \]
\[ \begin{equation}
\frac{h_n^{k+1} - h_n^k}{\Delta t} + \frac{1}{r_n} \frac{r_{n+1} h_{n+1}^{k+1} \langle u\rangle_{n+1}^{k+1} - r_{n-1} h_{n-1}^{k+1} \langle u\rangle_{n-1}^{k+1}}{2\Delta r} = - J_n^k,
\end{equation} \tag{19} \]
\[ \begin{equation}
\langle u\rangle_n^k = -H_a(x_n)\frac{(h_n^k)^2}{3 \eta_n^k} \frac{\langle P \rangle_{n+1}^k - \langle P \rangle_{n-1}^k}{2 \Delta r},
\end{equation} \tag{20} \]
\[ \begin{multline}
\frac{\langle C \rangle_n^{k+1} - \langle C \rangle_n^k}{\Delta t} + 
\frac{|\langle u\rangle_n^{k+1}|+\langle u\rangle_n^{k+1}}{2} \frac{ \langle C \rangle_n^k - \langle C \rangle_{n-1}^k}{\Delta r} - \frac{|\langle u\rangle_n^{k+1}|-\langle u\rangle_n^{k+1}}{2} \frac{ \langle C \rangle_{n+1}^k - \langle C \rangle_n^k}{\Delta r} = {} \\
{} = \frac{J_n^k \langle C \rangle_n^k}{h_n^{k+1}} + \frac{1}{\mathsf{Pe}} \frac{H_a(x_n)}{r_n h_n^{k+1} \Delta r} 
\Bigl( r_{n+1/2} h_{n+1/2}^{k+1} \frac{ \langle C \rangle_{n+1}^k - \langle C \rangle_n^k }{\Delta r} - r_{n-1/2} h_{n-1/2}^{k+1} \frac{ \langle C \rangle_n^k - \langle C \rangle_{n-1}^k }{\Delta r}\Bigr),
\end{multline} \tag{21} \]
\[ \begin{equation}
 \langle P \rangle_n^k = -\frac{1}{\mathsf{Ca}} \frac{1}{\Delta r r_n} \Bigl( r_{n+1/2} \frac{h_{n+1}^k - h_n^k}{\Delta r} -r_{n-1/2} \frac{h_n^k - h_{n-1}^k}{\Delta r} \Bigr),
\end{equation} \tag{22} \]
где $x_n = C_g - \langle C \rangle_n^k - \delta x$ и $k_h = 10/\delta x = 2000$. Промежуточные значения между узлами вида $f_{n\pm 1/2}$ рассчитываются как $(f_{n\pm 1} + f_n)/2$. Стоит заметить, что в (21) для дискретизации конвективного слагаемого используется схема с разностями против потока. Для приближения пространственной производной первого порядка в (17), (19) и (20) используется центральная разностная схема. Это дает второй порядок аппроксимации по $\Delta r$, как и аппрокисмация производной второго порядка в (18), (21) и (22). Здесь используется как явная разностная схема (21), так и неявная — (19).

Теперь запишем начальные и граничные условия из п. 4 в дискретном виде. Начальная форма свободной поверхности в дискретном виде задается как $h_n^0 = 1 - r_n^2$. Формула начального распределения массовой доли частиц принимает вид
\[ \begin{equation*}
\langle C \rangle_n^0 = C_g \frac{2 - C_0 + 2(C_0 -1)}{1 + \exp(w_c (r_n - 1))}.
\end{equation*} \]

С граничными условиями первого рода все просто, $f_{n,m}^k = c$, где $f$ — искомая функция и $c$ — некоторая константа (здесь $n=0$ или $n=N$). Для граничных условий второго рода вида $\partial f/ \partial r = 0$ используются аппроксимации второго порядка точности: $f_{0,m} = (4 f_{1,m} - f_{2,m}) / 3$, $f_{N,m} = 4 f_{N-1,m} - 3 f_{N-2,m}$. Для величин, не зависящих от $z$, индекс $m$ в этих формулах отсутствует.

Дискретное уравнение (19) можно расщепить по физическим процессам на два уравнения. Одно из них описывает изменение $h$ в результате испарения
\[ \begin{equation}
 \frac{h_n^\mathrm{tmp} - h_n^k}{\Delta t} = - J_n^k,
\end{equation} \tag{23} \]
другое — в результате конвективного переноса массы:
\[ \begin{equation}
 \frac{h_n^{k+1} - h_n^\mathrm{tmp}}{\Delta t} + \frac{1}{r_n} \frac{r_{n+1} h_{n+1}^{k+1} \langle u\rangle_{n+1}^{k+1} - r_{n-1} h_{n-1}^{k+1} \langle u\rangle_{n-1}^{k+1}}{2\Delta r} = 0.
\end{equation} \tag{24} \]
Здесь (23) и (24) — явная и неявная разностные схемы соответственно. Временное значение высоты капли $h_n^\mathrm{tmp}$ определяется из (23). Затем, после пересчета давления (22) и усредненной скорости (20), решается (24) и находится итоговое значение $h_n^{k+1}$ на новом временном слое. Уравнение (24) можно переписать в виде
\[ \begin{equation}
a_n h_{n-1}^{k+1} + b_n h_n^{k+1} + c_n h_{n+1}^{k+1} = d_n,
\end{equation} \tag{25} \]
где $a_n = - \Delta t r_{n-1} \langle u\rangle_{n-1}^{k+1}/(2 r_n \Delta r)$, $b_n = 1$, $c_n = \Delta t r_{n+1} \langle u\rangle_{n+1}^{k+1}/(2 r_n \Delta r)$ и $d_n  = h_n^\mathrm{tmp}$. В совокупности для всех узлов $n$ из (25) получаем СЛАУ, которую можно представить в матричной форме:
\[ \begin{equation*}
\begin{pmatrix}
 b_0& c_0 & 0 & 0 &\ldots & 0\\
 a_1& b_1 & c_1 & 0 &\ldots & 0\\
 0 & a_2 & b_2 & c_2 &\ldots & 0\\
 \vdots& \vdots& \ddots &\ddots &\ddots & \vdots\\
 0 & \ldots & 0 & a_{N-1} & b_{N-1} & c_{N-1}\\
 0 & 0 & 0 & \ldots & a_N & b_N
\end{pmatrix}
\cdot
\begin{pmatrix}
 h_0^{k+1}\\
 h_1^{k+1}\\
 \vdots\\
 h_{N-1}^{k+1}\\
 h_N^{k+1}
\end{pmatrix}
=
\begin{pmatrix}
 d_0\\
 d_1\\
 \vdots\\
 d_{N-1}\\
 d_N
\end{pmatrix} .
\end{equation*} \]

Запишем последнее выражением в более компактном виде: 
\[ \begin{equation*}
\mathbf{A}\cdot \mathbf{h} = \mathbf{d},
\end{equation*} \]
где $\mathbf{A}$ — матрица коэффициентов, $\mathbf{d}$ — вектор правых частей и $\mathbf{h}$ — искомый вектор. Обратим внимание, что коэффициенты $a_0$ и $c_N$ в вычислениях не задействованы, квадратная матрица $\mathbf{A}$ характеризуется ленточной трехдиагональной структурой ненулевых элементов и такая СЛАУ может быть эффективно решена методом прогонки. На основании граничных условий запишем значения коэффициентов $b_0$, $c_0$, $d_0$, $a_N$, $b_N$ и $d_N$. Чтобы не нарушить трехдиагональную структуру матрицы $\mathbf{A}$, для граничного условия на оси симметрии $\partial h/ \partial r = 0$ используем аппроксимацию первого порядка $(h_1^{k+1} - h_0^{k+1})/ \Delta r \approx 0$. В таком случае получаем коэффициенты $b_0 = 1$, $c_0 = -1$ и $d_0 = 0$. Условие на трехфазной границе $h = h_f$ приводит к коэффициентам $a_N = 0$, $b_N = 1$ и $d_N = h_f$.

Значения $h_n^{k+1}$, $\langle u\rangle_{n}^{k+1}$ и $\langle P\rangle_{n}^{k+1}$ необходимо итеративно уточнять. На первой итерации ($i=1$) задаем отправное значение $h_n^{k+1, 1} = h_n^\mathrm{tmp}$. Относительно резкое изменение давления за дискретный временной шаг в результате испарения согласно 
уравнению (23) приводит к неустойчивости вычислений. По этой причине здесь используется метод явной релаксации 
\[ \begin{equation*}
\langle P\rangle_{n}^{k+1,i} = \langle P\rangle_{n}^{k+1,i-1} + \alpha (\langle P\rangle_{n}^{k+1,i} - \langle P\rangle_{n}^{k+1,i-1} ),
\end{equation*} \]
где $i$ — номер итерации, а $\alpha$ — фактор релаксации. Предварительно $\langle P\rangle_{n}^{k+1,i}$ рассчитывается по формуле (22). Значение фактора релаксации подобрано эмпирическим путем. Здесь в расчетах использовалось значение $\alpha= 10^{-3}$. Малое значение этого параметра улучшает устойчивость численного счета, но при этом увеличивает количество итераций, необходимых для сходимости решения. Вблизи границ при расчете усредненной скорости $\langle u\rangle$ по 
формуле (20) возникают осцилляции, которые за несколько временных шагов приводят к «разрушению» вычислений. С целью предотвращения таких осцилляций для приграничных узлов использовалась линейная интерполяция. Учитывая, что $\langle u\rangle \to 0$ при $r \to 0$ или при $r \to R$, получаем $\langle u\rangle_1 = \langle u\rangle_2 / 2$ и $\langle u\rangle_{N-1} = \langle u\rangle_{N-2} / 2$.

Механизм остановки итераций включает два возможных события: достигнуто максимальное число итераций $I_\mathrm{max}=100$ или максимальная невязка стала меньше заданной погрешности $r_\mathrm{max} < \epsilon_\mathrm{lim}$, где $\epsilon_\mathrm{lim} = 10^{-17}$ и
\[ \begin{equation*}
r_\mathrm{max} = \left| \underset{n}{\max} \Bigl( h_n^{k+1} - h_n^\mathrm{tmp} + \frac{\Delta t}{r_n} \frac{r_{n+1} h_{n+1}^{k+1} \langle u\rangle_{n+1}^{k+1} - r_{n-1} h_{n-1}^{k+1} \langle u\rangle_{n-1}^{k+1}}{2\Delta r} \Bigr)\right|
\end{equation*} \] 
согласно уравнению (24).

После того как итерации завершились, полученные значения $h_n^{k+1}$, $\langle u\rangle_{n}^{k+1}$ и $\langle P\rangle_{n}^{k+1}$ используются для нахождения $\langle C \rangle_n^{k+1}$, $u_{n,m}^{k+1}$ и $w_{n,m}^{k+1}$ через явные зависимости (21), (17) и (18). Но предварительно необходимо подавить пилообразные осцилляции, возникающие при расчете давления. Этот вид осцилляций не приводит к «разрушению» численного счета со временем, но поле скорости потока жидкости получается некорректным. Чтобы разрешить эту проблему, сделаем следующее. Выражение (13) запишем в виде
\[ \begin{equation*}
\frac{\partial \langle P\rangle}{\partial r} = - 3\frac{\langle u\rangle \eta}{h^2}
\end{equation*} \]
и проинтегрируем его по $r$:
\[ \begin{equation*}
\langle P\rangle \big|_r - \langle P\rangle \big|_0 = -3 \int _{0}^{r} \frac{\langle u\rangle \eta}{h^2}\,dr.
\end{equation*} \]
Слагаемое с интегралом аппроксимируем методом прямоугольников и получим формулу для восстановления давления:
\[ \begin{equation}
\langle P\rangle_n^\mathrm{repair} = -3 \Delta r \sum\limits_{n=0}^{N-1} \frac{\langle u\rangle \eta}{h^2} + \langle P\rangle_0^\mathrm{repair}.
\end{equation} \tag{26} \]

Учитывая пилообразную форму осцилляций, $\langle P\rangle_0^\mathrm{repair}$ можно аппроксимировать следующим образом: 
\[ \begin{equation*}
\langle P\rangle_0^\mathrm{repair} \approx (\langle P\rangle_1 + \langle P\rangle_2 )/ 2.
\end{equation*} \]
Со временем такие пилообразные осцилляция могут передаваться усредненной скорости $\langle u\rangle$ и профилю капли $h$. Здесь эта проблема решена за счет использования сглаживающего фильтра:
\[ \begin{equation}
f_n^\mathrm{filt} = \beta f_{n+1} + \gamma f_{n} + \delta f_{n-1}.
\end{equation} \tag{27} \]
Параметры подобраны эмпирическим путем: $\beta = \delta = 0.25$ и $\gamma = 0.5$. В отличие от [49] используемый здесь фильтр прост в реализации, занимает в разы меньше памяти и вычислительных ресурсов. Кроме того, количество итераций фильтра из [49] заранее неизвестно. В нашем случае достаточно одного прохода по узлам $n=1, 2, \dots, N-1$. Узлы $n=0$ и $n=N$ рассчитываются из граничных условий.

Блок-схема с алгоритмом¹ расчета представлена на рис. 2. В ней используются следующие дополнительные обозначения: $p_c$ — счетчик временных периодов, $p_n$ — количество временных периодов, $p_s$ — размер временного периода во временных шагах и $H_n = H_a (x_n)$. Записи вида $f_n^k = f_n^{k+1}$ следует трактовать как векторные операции.

Рис. 2. Алгоритм численного расчета [Figure 2. Numerical calculation algorithm]

Для сравнения численный расчет также выполнялся в коммерческом пакете Maple 18 с помощью модуля pdsolve с параметром $method= Theta[1]$ (неявная разностная схема). Для того чтобы pdsolve справился с задачей, систему уравнений пришлось переписать в другом виде. Вместо (13) и (14) использовались уравнения
\[ \begin{equation}
\psi = - \frac{1}{\mathsf{Ca}} \frac{\partial h}{\partial r},
\end{equation} \tag{28} \]
\[ \begin{equation}
\langle u\rangle = -H_a(x)\frac{h^2}{3 \eta} \frac{\partial}{\partial r}\Bigl( \frac{1}{r} \frac{\partial (r \psi)}{\partial r} \Bigr),
\end{equation} \tag{29} \]
где $\psi$ — вспомогательная функция. Причем модуль pdsolve не требует граничных условий для (28). Скорее всего, модуль определяет их автоматически из заданных условий для $h$. Для уравнения (29) на границах задается равенство скорости нулю (см. п. 4). Таким образом, в Maple решалась система уравнений (12), (16), (28) и (29).

6. Результаты и обсуждение

Результаты расчетов средствами двух программ сравниваются на рис. 3. На большинстве графиков расхождение практически незаметно. Максимальное отличие результата, полученного с помощью программы, написанной на С++, от результата, полученного с помощью модуля pdsolve в Maple 18, наблюдается для $\langle u\rangle$ на времени $t = 220$ с (рис. 3, c). В абсолютных величинах наибольшее расхождение соответствует точке $r\approx 0.75$ мм и составляет примерно 0.046 мкм/с, в относительных величинах это около 10 %.

Рис. 3. Сравнение результатов расчетов в двух программах для нескольких последовательных моментов времени: (a) высота капли, (b) массовая доля растворенного или взвешенного вещества и (c) радиальная скорость потока, усредненная по толщине жидкого слоя
[Figure 3. Comparison of the calculation results in two programs for several consecutive time points: (a) drop height, (b) mass fraction of dissolved or suspended matter, and (c) radial flow velocity averaged over the thickness of the liquid layer]

Расчет выполнялся для $N=75$. Значения временных шагов для разработанной программы и для модуля pdsolve брались разные, так как в них используются разные алгоритмы расчета, соответственно и требования к временным шагам накладываются тоже разные. В Maple 18 для численного решения СЛАУ используется итерационный метод Ньютона, но технические детали реализации закрыты для пользователя, есть лишь краткая справка по разностным схемам. Скорость вычислений с относительно малым временным шагом, конечно же, ниже, но это частично компенсируется тем, что C++ является компилируемым языком, а язык Maple — интерпретируемый. Для расчетов в Maple использовалось значение $\Delta t = 0.1$ с, в разработанной программе — $\Delta t =10^{-6}$ с. Кроме того, для программы на C++ выполнены тестовые расчеты для $N = 25$ и $N = 50$ с целью проверки сходимости по сетке (во всех расчетах $M=N$). Для $N = 25$ приближенно подобрано максимальное значение временного шага $\Delta t =10^{-4}$ с, при котором расчет не «разрушается». Значению $N = 50$ соответствует $\Delta t =5\cdot 10^{-6}$ с. 
Таким образом, небольшое уменьшение $\Delta r= R/N$ приводит к тому, что необходимо значительно уменьшать временной шаг. Эти наблюдения в вычислительных экспериментах косвенно свидетельствуют о том, что предложенная разностная схема является условно устойчивой. Длительность расчета в Maple составила около 10 мин для $N=75$ (объем занимаемой памяти $\approx 217$ Мб). Время расчета в программе на С++ примерно составило 10 мин для $N = 25$, 5 ч для $N = 50$ и 29 ч для $N=75$ (объем занимаемой памяти $\approx 3$ Мб). Расчеты проводились на компьютере с ЦП Intel Xeon CPU E3-1230 v3 c максимальной частотой 3.6 ГГц под управлением ОС Windows 10. Для получения исполняемого файла программы на С++ использовался компилятор, встроенный в Visual Studio 2022 Community Edition.

Проверка закона сохранения массы для растворенного или взвешенного вещества для $N=75$ показала погрешность около 0.27 %, что экспериментально подтверждает свойство консервативности в предложенной разностной схеме. Для этого были рассчитаны и сопоставлены значения 
\[ \begin{equation*}
\sum\limits_{n=0}^{N} h_n \langle C \rangle_n
\end{equation*} \]
на моментах времени $t=0$ и $t=t_f$. Для расчетов в Maple также проводилась проверка закона сохранения массы. Выполнены тестовые расчеты и определены погрешности для разных временных шагов в диапазоне, для которого численное решение итерационно сходится: $\Delta t = 0.01$ с — 0.43 %, $\Delta t = 0.1$ с — 0.1 %, $\Delta t = 1$ с — 3.11 % и для $\Delta t = 10$ с — 32.54 %. Минимальная погрешность в 0.1 % соответствует временному шагу $\Delta t = 0.1$ с, который и был выбран для расчетов в Maple.

Рис. 4. Результаты расчета в Maple: пространственно-временная зависимость (a) усредненной по толщине жидкого слоя скорости радиального потока и (b) плотности потока пара
[Figure 4. The results of the calculation in Maple: the spatial-temporal dependence of (a) the radial flow velocity averaged over the thickness of the liquid layer and (b) the vapor flux density]

Толщина жидкого слоя с течением времени уменьшается (рис. 3, a). Из-за пространственно неравномерного испарения происходит отклонение формы капли от равновесного состояния, возникает градиент капиллярного давления, что порождает капиллярный поток (рис. 3, c), переносящий вещество к периферии капли. В связи с этим массовая доля вещества, содержащегося в жидкости, увеличивается вблизи трехфазной границы со временем (рис. 3, b). Под конец процесса массовая доля $\langle C \rangle$ распределена равномерно, так как достигнуто критическое значение $C_g$, при котором произошел фазовый переход. Дальнейшее испарение жидкости из твердой фазы здесь не рассматривается.

Финальная расчетная толщина слоя $h$ является почти равномерной в пространстве. Вблизи $r\approx 0.9$ мм наблюдается незначительное возвышение, которое связано с эффектом кофейных колец [3]. Возникающий под конец процесса обратный поток жидкости, направленный в центральную область, смазывает этот эффект. Это можно использовать для улучшения существующих методов подавления эффекта кофейных колец (см. ссылки в [1]).

Рис. 5. Векторное поле скорости потока жидкости по результатам, полученным с использованием программы, написанной на C++, для следующих времен: (a) \(t=10\) с, (b) \(t=90\) с, (c) \(t=220\) с и (d) \(t=300\) с
[Figure 5. Velocity vector field of fluid flow according results obtained with using the program written in C++ for the following times: (a) $t= 10$ s, (b) $t= 90$ s, (c) $t= 220$ s, and 
(d) $t= 300$ s]

Рис. 6. Результаты расчета в программе, написанной на С++: (a) капиллярное давление и (b) возможные пилообразные осцилляции
 [Figure 6. The results of the calculation in program written with C++: (a) capillary pressure and (b) possible sawtooth oscillations]

Также стоит заметить, что поначалу усредненная радиальная скорость потока с течением времени увеличивается в связи с увеличением градиента капиллярного давления в процессе неравномерного испарения, преобладающего вблизи периферии. Затем в области роста массовой доли вещества, приводящего к увеличению вязкости, поток начинает замедляться. Вблизи трехфазной границы появляется противоток (рис. 3, c для $t = 220$ с), что ранее наблюдалось в эксперименте [52]. Далее по времени направление потока постепенно меняется во всей области (рис. 4, a). Это связано с тем, что в начале процесса плотность потока пара преобладает в районе периферии, но со временем $J$ начинает увеличиваться в центральной области из-за понижения $h$ и уменьшаться вблизи контактной линии из-за роста $\langle C \rangle$ (рис. 4, b). Когда во всей области достигается значение $\langle C \rangle \approx C_g$, модель предсказывает отсутствие потока и испарения. Поле скорости потока для различных моментов времени представлено на рис. 5. В отличие от результатов расчетов [53] здесь на поздних временных этапах наблюдается изменение направления потока. Это связано с тем, что в настоящей работе учитывается зависимость $J$ и $\eta$ от массовой доли $\langle C \rangle$. Кроме этого, здесь при заданных параметрах градиент плотности потока пара с течением времени меняет знак. Технические детали численного решения в [53] не раскрыты. В [27, 28, 49] поле скорости потока не рассчитывалось.

Здесь поле скорости потока удалось рассчитать по формулам (17) и (18) лишь после того, как была решена проблема с пилообразными осцилляциями, возникающими при расчете капиллярного давления и постепенно с временными шагами накапливающимися в $\langle u\rangle$ и $h$. Для этого использовались формулы (26) и (27). Восстановленное из скорости потока $\langle u\rangle$ по формуле (26) капиллярное давление представлено на рис. 6, a. Если такое восстановление не проводить, то получаем осцилляции, представленные на рис. 6, b. Особенность этих осцилляций заключается в стабильности на временных шагах, то есть они не приводят к полному «разрушению» расчета. Пилообразные осцилляции теоретически изучаются в [54], но здесь предложено практическое решение проблемы на конкретном примере.

Заключение

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

Конкурирующие интересы. Заявляю, что в отношении авторства и публикации этой статьи конфликта интересов не имею.
Авторская ответственность. Я несу полную ответственность за предоставление окончательной версии рукописи в печать. Окончательная версия рукописи мною одобрена.
Финансирование. Часть работы, относящаяся к разработке математической модели, выполнена при поддержке программы развития университета «Приоритет 2030» проект № 122112500011–0 «Природовдохновленные оптические технологии» (Центр природовдохновленного инжиниринга ТюмГУ). Исследование, связанное с разработкой численного метода и программы для ЭВМ, выполнено за счет гранта Российского научного фонда № 22–79–10216 в АГУ им. В. Н. Татищева, https://rscf.ru/project/22-79-10216/.
Благодарность. Автор выражает благодарность профессору, доктору физико-математических наук Ю. Ю. Тарасевичу за полезные рекомендации по оформлению графического материала.

---

¹Код программы на языке С++ доступен по ссылке https://github.com/kolegovk/Suppression-of-sawtooth-oscillations.git

About the authors

Konstantin S. Kolegov

Astrakhan Tatishchev State University; Tyumen State University

Author for correspondence.
Email: konstantin.kolegov@asu.edu.ru
ORCID iD: 0000-0002-9742-1308
SPIN-code: 1872-4975
Scopus Author ID: 57191072458
ResearcherId: T-6123-2017
https://science.asu.edu.ru/index.php/user/9/

Cand. Phys. & Math. Sci., Senior Researcher, Lab. of Mathematical Modeling and Information Technologies in Science and Education; Center for Nature-inspired Engineering

Russian Federation, 414056, Astrakhan, Tatischev st., 20 a; 625003, Tyumen, Lenin st., 25

References

Supplementary files