CC BY
2
Math-Net.Ru
К. С. Колегов, Подавление пилообразных осцилля-ций при использовании разностной схемы для моделирования массопереноса в высыхающей на подложке капле в приближении тонкого слоя, Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2023, номер 2, 309-335
001: 10.14498/уБ^и1994
Использование Общероссийского математического портала Math-Net.Ru подразумевает, что вы прочитали и согласны с пользовательским соглашением
http://www.mathnet.ru/rus/agreement
Параметры загрузки:
IP: 109.252.33.182
29 сентября 2024 г., 12:09:41
Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. 2023. Т. 27, № 2. С. 309-335 ISSN: 2310-7081 (online), 1991-8615 (print) d https://doi.org/10.14498/vsgtu1994
EDN: HZRAJS
УДК 519.635.4
Подавление пилообразных осцилляций при использовании разностной схемы для моделирования массопереноса в высыхающей на подложке капле в приближении тонкого слоя
К. С. Колегов1'2
1 Астраханский государственный университет имени В. Н. Татищева, Россия, 414056, Астрахань, ул. Татищева, 20 а.
2 Тюменский государственный университет, Россия, 625003, Тюмень, ул. Ленина, 25.
Аннотация
Испаряющиеся капли и пленки используются в приложениях из разных областей. Особый интерес представляют различные методы испарительной самосборки. В работе описана математическая модель массо-переноса в высыхающей на подложке капле на базе приближения тонкого слоя. Модель учитывает перенос растворенного или взвешенного вещества капиллярным потоком, диффузию этого вещества, испарение жидкости, формирование твердого осадка, зависимость вязкости и плотности потока пара от концентрации примеси.
Рассматривается случай, когда трехфазная граница «жидкость-подложка-воздух» закреплена. Для уравнений модели разработаны явные и неявные разностные схемы. Предложена модификация численного метода, в которой комбинируется расщепление по физическим процессам, итерационный метод явной релаксации и метод прогонки. Описан практический рецепт подавления пилообразных осцилляций на примере конкретной задачи.
Разработан программный модуль на языке С++, который в дальнейшем можно использовать для задач испарительной литографии. С помощью этого модуля проведены численные расчеты, результаты которых сравнивались с результатами, полученными в пакете Maple.
Численное моделирование предсказало случай, когда направление капиллярного потока с течением времени меняется на противоположное из-за изменения знака градиента плотности потока пара. Это может
Математическое моделирование, численные методы и комплексы программ Научная статья
© Коллектив авторов, 2023 © СамГТУ, 2023 (составление, дизайн, макет)
3 ©® Контент публикуется на условиях лицензии Creative Commons Attribution 4.0 International (https://creativecommons.org/licenses/by/4.0/deed.ru) Образец для цитирования
КолеговК. С. Подавление пилообразных осцилляций при использовании разностной схемы для моделирования массопереноса в высыхающей на подложке капле в приближении тонкого слоя // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2023. Т. 27, № 2. С. 309-335. EDN: HZRAJS. DOI: 10.14498/vsgtu1994. Сведения об авторе
Константин Сергеевич Колегов © https://orcid.org/0000-0002-9742-1308 кандидат физико-математических наук; старший научный сотрудник; лаб. математического моделирования и информационных технологий в науке и образовании1; центр природо-вдохновленного инжиниринга2; e-mail:[email protected]
приводить к замедлению выноса вещества на периферию, что в результате будет способствовать формированию более или менее равномерного осадка по всей площади контакта капли с подложкой. Данное наблюдение полезно для совершенствования методов подавления кольцевых осадков, связанных с эффектом кофейных колец и нежелательных для некоторых приложений, как, например, струйная печать или нанесение покрытий.
Ключевые слова: испаряющаяся капля, массоперенос, капиллярный поток, разностная схема, пилообразная осцилляция, подавление эффекта кофейных колец.
Получение: 18 января 2023 г. / Исправление: 29 апреля 2023 г. / Принятие: 3 мая 2023 г. / Публикация онлайн: 27 июня 2023 г.
Введение. Испаряющиеся капли и пленки используются в приложениях из разных областей, например, охлаждение нагретых поверхностей электронных приборов, диагностика в медицине, формирование прозрачных электропроводных покрытий на гибкой подложке, структурирование поверхности [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]. Но эти модели либо не учитывают гидродинамику, либо используют простые аналитические решения для частного случая. Для улучшения этих моделей, алгоритмов и программ необходимо разрабатывать дополнительные модули, позволяющие учесть различные эффекты, влияющие на формирование осадков, в том числе и гидродинамику. В коммерческом пакете ^швс! Multiphysics есть модули, позволяющие рассчитывать динамику частиц и гидродинамику [46], но цена на этот пакет относительно высокая. Это касается и дополнительных тулбоксов к Matlab [47]. Недавно для нашей страны был ограничен доступ к официальному сайту свободного пакета для молекулярной динамики LAMMPS (https://www.lammps.org/). Зарубежные коммерческие пакеты в любой момент могут оказаться недоступными из-за политической ситуации. К примеру, на момент написания данной статьи приостановлены продажи пакета Ansys в РФ. В условиях импортозамещения очень важной является разработка отечественного программного обеспечения для моделирования, написание своих кодов и библиотек. Для проведения исследований в области испарительной литографии требуется разработка программного комплекса, включающего ряд модулей (рис. 1). Пустой прямоугольник на схеме обозначает дополнительные модули, которые необходимо создать. К примеру, в случае необходимости учета таких важных факторов, как взаимодействие частиц с подложкой (адсорбция или адгезия) или со свободной поверхностью капли (капиллярные силы).
Разработка модулей для моделирования диффузии и капиллярного притяжения сферических частиц обсуждалась ранее в [42, 44, 45]. Также был рассмотрен случай для бинарной смеси частиц разного размера [13]. Цель текущей работы — разработать численный алгоритм и программный модуль для расчета гидродинамики в высыхающей капле для дальнейшего использования в программном комплексе, ориентированном на решение задач в области испарительной литографии. Но при разработке необходимо учитывать одно требование. Алгоритм должен быть рассчитан на большое количество предельно малых временных шагов (минимальные вычислительные затраты на каждом временном шаге), что позволит корректно работать уже имеющемуся модулю «Диффузия частиц», основанному на методе Монте—Карло, для явного предсказания динамики частиц, так как в этом случае выполняется соотношение Эйнштейна—Смоулховского [45]. Величина временного шага
косвенно зависит от размера частиц, например, для частиц с радиусом
Рис. 1. Схема планируемого программного комплекса для решения задач в области испарительной литографии [Figure 1. Scheme of the planned software complex for solving problems
in the field of evaporative lithography]
0.35 мкм в [45] использовался шаг = 10-4 с.
Также здесь стоит отметить, по какой причине не подходят существующие готовые решения. Описанное в [48] аналитическое решение базируется на смеси кинематического подхода и приближения смазки, которое еще известно как приближение тонкого слоя. Под кинематическим подходом здесь понимается нахождение усредненной радиальной скорости потока из закона сохранения массы. Такой подход не объясняет природу возникновения потока, как например, капиллярные силы, когда в задаче явным образом учитывают градиент давления Лапласа и кривизну свободной поверхности. Как видно из полученного решения [48], оно явным образом не зависит от ряда физических параметров, как, например, вязкость жидкости. Кроме того, это решение опирается на частный случай с определенным видом функции плотности потока пара ■](г,Ь), что не позволяет решать задачи, связанные с испарительной литографией. Здесь г — радиальная координата и £ — время. Различные маски, излучатели и прочее внешнее воздействие оказывают влияние на концентрацию пара вблизи свободной поверхности, что в итоге влияет на капиллярный поток жидкости [2]. Таким образом, для определения .](г, ¿) необходимо численно моделировать перенос пара в воздухе. Предложенная в [49] неявная разностная схема является неустойчивой. Вблизи границ появляются осцилляции, которые со временем разрастаются на всю область. Авторы [49] подробно не разъясняют способ численной реализации граничных условий, ссылаясь на приложенный код в дополнении. Но восстановить некоторые граничные условия из кода крайне сложно. В программе [49] на каждом временном шаге выполняется искусственное сглаживание осциллирующих точек с помощью специального фильтра. Кроме того, на каждом шаге необходимо выполнять такие вычислительно затратные операции, как нахождение обратной матрицы и умножение матрицы на вектор. По этим причинам данный программный модуль [49] нам не подходит. В отличие от [27, 28, 49] здесь дополнительно будет описана численная реализация расчета поля скорости потока жидкости, а не только усредненной по толщине жидкого слоя радиальной компоненты скорости.
1. Физическая постановка задачи. Капля жидкости испаряется на твердом горизонтальном непромокаемом основании при нормальных комнатных условиях. Считаем, что жидкость несжимаема. Трехфазная граница закреплена, поэтому радиус контакта капли и подложки К — постоянная величина. Рассматриваем случай относительно малого объема жидкости, когда сила поверхностного натяжения преобладает над силой тяжести, число Бонда Во = рЪ0д/а ~ 1.4 ■ 10-3 ^ 1, где д — ускорение свободного падения, р — плотность жидкости, а — коэффициент поверхностного натяжения и Но = Н(0, 0) —начальная высота капли. Здесь функция Н = Н(г,Ь) задает профиль свободной поверхности жидкости (граница «жидкость - воздух»). В этом случае форма капли близка форме сферического сегмента. По причине осевой симметрии удобно использовать цилиндрическую систему координат (г, г). Движение жидкости вызывает локальные изменения кривизны поверхности в процессе испарения жидкости, приводящие к возникновению градиента давления Лапласа (капиллярный поток). В случае малых значений краевого угла смачивания в капиллярный поток преобладает над потоком Марангони (см. оценку в [44]). Здесь речь идет о тонких каплях, для которых аспектное отношение е = Но/К ^ 1. Значения параметров задачи приведены в табл. 1.
Распределение растворенного или взвешенного в жидкости вещества будем описывать функцией (С)(г,Ь). Предполагаем, что массовая доля вещества (С) не зависит от координаты г. Если время диффузионной релаксации вещества много меньше времени полного испарения, т.е. ^ tf, то устанавливается равномерная концентрация вещества вдоль вертикального направления. Здесь время = Но/О = 100 с и время tf = 450 с. Таким образом, величина (С) является усредненной по высоте жидкого слоя.
2. Математическая модель. Запишем уравнение неразрывности и стационарные уравнения Навье—Стокса без учета инерционных слагаемых в цилиндрической системе координат:
1 д (ги) дад
г дг дг '
д / ^д(ги) \ д / ди \
дг\г дг ) дг \ дг)
1 д / дад \ д / дад \
г дг
дг ) + дг
дР
дг ' дР
дг '
(1) (2) (3)
где ц = ц((С)) —динамическая вязкость жидкости. Давление Р, горизонтальная и вертикальная компоненты вектора скорости V = (и, ад) являются функциями, зависящими от ¿, г и г. Система уравнений (1)-(3) справедлива для несжимаемой жидкости и малых значений числа Рейнольдса Ре = рисНо/^о. Как правило, характерная скорость потока жидкости ис ~ 1 мкм/с в испаряющейся при комнатных условиях капле воды. В таком случае число Рейнольдса Ре « 10-4 ^ 1.
Для перехода к безразмерным записям рассмотрим масштабные соотношения: ц = щг}, 3 = 33с, и = иси, ад = еисад, £ = ¿Д Р = РСР, г = Кг и г = Ног. Здесь знаком «тильда» помечаем безразмерные величины. Заметим, что горизонтальные и вертикальные составляющие (размер, скорость)
Таблица 1
Физические и геометрические параметры задачи [Physical and geometric parameters of the problem]
Обозначение Расшифровка Величина Единицы измерения
[Designation] [Description] [Value] [Units of Measurement]
h0 высота капли [drop height] 10-4 m
R радиус основания капли [drop base radius] 1(Г3 m
hf толщина жидкого слоя в районе контактной линии [thickness O.Olho m
of the liquid layer in the area of the contact line]
e = ho/R аспектное отношение [aspect ratio] 0.1 -
в « 2e краевой угол смачивания [contact angle] 0.2 radians
C0 начальная массовая доля [initial mass fraction] 0.05 -
cg концентрация фазового перехода, например, гелеобразования 0.7 -
[phase transition concentration, e.g. gelation] ю-10
D коэффициент диффузии растворенного или взвешенного ве- m2/s
щества [diffusion coefficient of a dissolved or suspended matter]
m вязкость чистой жидкости [viscosity of a pure liquid] 10~3 Pa • s
p плотность жидкости [liquid density] 103 kg/m3
a коэффициент поверхностного натяжения [surface tension 0.072 N/m
coefficient]
Ь время полного испарения [total evaporation time] 450 s
Dv коэффициент диффузии пара [vapor diffusion coefficient] 2.4 • 10-5 m2/s
Cv плотность насыщенного пара [saturated steam density] 2.32 • 10-2 kg/m3
H относительная влажность [relative humidity] 0.4 -
масштабируются по-разному. Такой подход лежит в основе приближения смазки.
Зададим характерные величины: скорость ис = щ/(рЬо) ~ 0.01 м/с, время ¿с = К/ис & 0.1 с, плотность потока пара ,]с = еиср & 1 кг/(м2 ■ с) и давление Рс = Кисщ/к2 & 1 Па. Запишем уравнение неразрывности (1) в безразмерном виде:
+ Оы =о. (4)
г дг дг
Заметим, что уравнение (4) не отличается по форме записи от (1). Уравнения (2) и (3) записываются так:
д /fjd (ru) \ д /дй \ дР
дг\г дг ) dz \ dz) дг'
31 д дй\ 2 д /дй\ дР „
£'i-^z[Vii'' -трг + t^ziv тгг = ^гг. (6)
£ дг
г дг\ дг J дг\дг/ дг Пренебрегая слагаемыми, при которых стоят е2 и е3, в (5) и (6), получаем
дР д ( дй \
& - Ш\Рд2 )=°' (7)
д Р
дР = 0 (8)
В результате получена упрощенная система уравнений (4), (7) и (8).
Из (8) следует, что давление Р не зависит от z, тогда вместо Р будем рассматривать усредненную по высоте жидкого слоя величину (Р) = (Р)(г, t). Интегрируя левую и правую части (7) и учитывая, что вязкость "q здесь не зависит явно или косвенно от ,
^J -fl J^d~z = j0 d~z'
получаем
(* + Ci) -<I + C) = Сз.
Как было отмечено ранее, термокапиллярные потоки здесь не рассматриваются, поэтому на свободной границе z = h выполняется баланс касательных напряжений дй/дг = 0, тогда
д(Р) г д(РЪ ^ _ h + AtzLC i =Сз + г]С2. д д
Для г = 0 выполняется условие осевой симметрии д(Р)/дг = 0, тогда С3 = — fj(0' t)C2. В действительности вязкость fj зависит от г и t неявным образом, через массовую долю (С), но здесь для краткости пишем fj(г, t). На границе г = R горизонтальная компонента скорости является постоянной величиной, не зависящей от z (в силу прилипания к подложке й = 0).
Тогда в этой точке д2й/дг2 = 0. Учитывая это, из (7) получаем граничное условие д(Р)/дг = 0 для г = К. С таким условием приходим к соотношению Сз = —г](К, ¿)С*2. На разных границах получили одно и то же соотношение, хотя в общем случае г/(0, Ь) = г/(К, ¿). Приходим к выводу, что С2 = Сз = 0. В таком случае константа интегрирования С\ = —Н и последнее выражение принимает вид
д(Р) . ~. ди
у. и — К) — „_ = 0.
Интегрируем полученное еще раз:
[ 1й1 — Н [ йг — г! [^йг = 10йх, от ) от ) ) дг )
получаем
Щ) (I2 + ^4 — } + С 55) — П(й + Сб) = Ст.
Учитывая, что для г = 0 выполняется условие прилипания и = 0 и для г = 0 выполняется условие д(Р)/дг = 0, получаем Ст = — г/(0, ¿)С6. Учитывая аналогичные граничные условия для г = К, получаем Ст = — ^¡(К, £)Сб. Приходим к выводу, что С б = Ст = 0. Затем, учитывая лишь условие и = 0 для г = 0, получаем С4 = — С5. В итоге выражаем горизонтальную компоненту вектора скорости:
и = НМ Ф- ( " — НЛ (9)
Г} дг 4 2/ У 7
В некоторых уравнениях модели, в том числе и в (9), дополнительно используется аналитическое приближение функции Хевисайда На, чтобы исключить некоторые виды массопереноса при возникновении фазового перехода, как, например, «золь-гель», «жидкость-стекло» и т.п.:
На(х) = 1
1 + ехр(—2ких)'
где х = Са — (С) — 5х, кь = 10/5х, 5х — ширина области перехода. Теперь, подставляя (9) в (4) и интегрируя полученное:
Г 1 д г^Шт —Нг))йг + I ^йг = [ 0йг,
} г дг V г] дг \ 2 П } дг }
получаем
1 д_ /ша(х) д(Р) /+С8\\ + Сд = 0
г дг\ V дг I 6 2 //
г] дг V 6 2 Запишем полученное выражение в развернутом виде:
На(х) д(Р) /_ Н^2ч + На(х) /Н2ч д2(Р) На(х)г2 д(Р) дк + Щ дг I 6 2 ) г/ \ 6 2 ) дг2 2г! дг дг
+ с8на(х) д(Р) + с8на(х) д2(Р) +
г} дг } дг2
+ ЩМ1 6 - О* +сл + г + с = 0.
дг г дг \ 6 2 /
} дг V 6 2
Учитывая условие непротекания го = 0 для Г = 0 и условие д(Р)/дг = 0 для Г = 0 и г = П, получаем два соотношения:
С8На(х) д2(Р )(0, {) + с 0 С8На (х) д2(Р )(КЛ) 0
гЛ----+ С =0 и —^--—2--+С = 0.
}(0, г) дг2 Г}(К, ¿) дг2
Из чего следует, что С8 = Сд = 0. В таком случае получаем выражение для вертикальной компоненты вектора скорости:
г = -1 ^Г-ННМ _ О*)). (10)
г дг\ Г дг \ 6 2//
Движение свободной поверхности капли в процессе испарения описывается законом сохранения (см. ссылки в
дО 1 д (гО(и)) _ ¿Г / дО \ 2 дЪ г дг р V \дг / '
где (и) — усредненная по высоте жидкого слоя скорость радиального потока жидкости. В безразмерном виде (11) записывается так:
д° + 1 д (г°(и)) = сЯ г дг '
где пренебрегаем ^ 1 + е2(дО/дг)2 в правой части как величиной второго порядка малости. С учетом (9) для усредненной радиальной скорости (и) имеем
¡■к
(и) = 1 Г и <12 = _На(х)02 ^. (13)
О]о 3 } дг
Капиллярное давление зависит от локальной кривизны поверхности:
/Г\ 1 1 д (_дО\
(Р) = _ Са1 дТгы, (14)
что следует из формулы для давления Лапласа [28] с учетом приближения смазки. Здесь капиллярное число Са = г}0ис/(ае3) = г}2/(е3рк0а) & 0.14. В [49] выражения (13) и (14) подставляются в уравнение (12) и затем строится разностная схема для уравнения с производной четвертого порядка по относительно О. Чтобы избежать таких громоздких разностных схем, здесь (13) и (14) рассматриваются как отдельные уравнения.
Для описания переноса растворенного или взвешенного вещества используем уравнение конвекции-диффузии
сС + {U дг rh дг ОГh дг ) + ph + ( дг ) ' ( 5)
где D — коэффициент диффузии растворенного или взвешенного вещества. Подробный вывод уравнения (15) описан в [50]. В безразмерном виде уравнение (15) с учетом приближения смазки принимает вид
ддС)+(й)дд^1 = eML (rhm \+т. (16)
dt + W дГ = Pe fhdA дГ )+ h ' (16)
где Pe — число Пекле, Pe ~ Ruc/D ~ 105. В итоге система уравнений включает в себя уравнения (9), (10), (12), (13), (14) и (16).
3. Вязкость и плотность потока пара. Теперь необходимо ввести замыкающие соотношения для некоторых величин. В предложенной здесь модели вязкость зависит от массовой доли вещества. Для описания этой зависимости используем формулу Муни
_ ( S {С) \
где значения эмпирических параметров S = 1.692 и К = 1.236 возьмем из [51]. Также считаем, что плотность потока пара зависит от концентрации вещества и толщины жидкого слоя:
J = Jo
1 - {С)2/С2
к + h
где Сд — критическая массовая доля, при которой происходит фазовый переход [27] (золь-гель, стеклообразование, кристаллизация и т. п.), полагаем Сд & 0.7 [51]. Для расчетов будем использовать значение параметра к = 1. Это некоторая аппроксимация для учета испарения, поэтому в дальнейшем планируется разработать модуль численного расчета переноса пара в воздухе (рис. 1). Величина Jо рассчитывается так:
,!о = °ьСьт(1~ Н) (0.2762 + 1.3)(0.6381 — 0.2239(6 — ж/4)2), исн
где Н — относительная влажность, Иу — коэффициент диффузии пара иС„ — плотность насыщенного пара [48].
4. Начальные и граничные условия. Запишем начальные условия задачи. Для тонких капель, размер которых меньше капиллярной длины, допустимо использовать параболическое приближение начальной формы свободной поверхности г, 0) = 1 — г2 [48]. Начальную массовую долю вещества опишем с помощью выражения [27]
(С) (?0) С2 — °о + 2(°о — 1) (°)(Г, 0) = °д 1+ехр(Мг — 1)),
где параметр регулирует начальную ширину кольцевого осадка (здесь в расчетах используем = 30). Такое условие предполагает, что вещество в начале процесса распределено практически равномерно вдоль радиальной координаты, за исключением района периферии (г ~ К = 1), где концентрация устремляется к пороговому значению Сд, что обеспечивает закрепление трехфазной границы «жидкость-подложка-воздух».
Для = 0 по причине осевой симметрии задаются следующие граничные условия:
д(Р) _дк _ д(С)
г=о дг
д
г=о дг
дад
г=о дг
= й\~ „ = (й)\. п = 0. г=о 1г=о х '1г=о
На границе г = К выполняются условия д( Р)
д
= Щ?=К = й\г=К = (П)\г=к = 0, Ъ\г=К = ^, (С) \ г=К = С9.
Условие для (Р) уже обсуждалось в п. 2. Вертикальная и горизонтальная компоненты вектора скорости потока жидкости записываются из соображений непротекания и прилипания. Здесь Ъ^ — толщина жидкого слоя в районе контактной линии.
Граничные условия вдоль подложки следующие:
„ = й\~ с = 0,
=о =о
по причине непротекания и прилипания. Вдоль свободной границы жидкости условия принимают вид
дп
= ьз\- г =0
д п д
х=Ь ,Х=Н
в предположении, что на этой границе нет вязкого трения и отсутствует поток жидкости через границу «жидкость-воздух» (из-за медленного испарения влияние подвижной границы и отдачи пара не учитываем). В приближении смазки на этой границе ад и П соответствуют нормальной и касательной компонентам вектора скорости. Последние условия для скорости потока на свободной поверхности капли являются излишними, так как значения П и ад на этой границе несложно рассчитать с использованием выражений (9) и (10).
5. Численный метод решения задачи. В этом пункте рассматриваются только безразмерные величины, поэтому знак «тильда» опускаем для краткости. Для дискретной аппроксимации уравнений модели с помощью разностных схем введем пространственно-временную сетку с координатами гп = пАг, гт = тАг и ^ = кА1, где п = 0,1,... ; т = 0,1,..., М; к = 0,1,..., К. Параметры А1, Аг и Аг обозначают шаги по времени и пространству соответственно. Кроме того, сетка характеризуется такими величинами, как количество отрезков ( N, М и К) и узлов (Ы + 1, М + 1 и К + 1). Параметр К = [£f /АЬ], здесь квадратные скобки обозначают целочисленное деление. Шаги по пространству определяются как Аг = К/Ы и Аг = Ъо/М.
Перепишем уравнения (9), (10), (12), (13), (14) и (16) в дискретном виде:
к = На(хп) (Р)П+1 — (Р)П-1 >пт = ^ 2 Аг \2 П^т
— ЪПит), (17)
к = _ На{хп+1/2) Гп+1/2 (Р)кп+1 - (Р)к / Ъп+1/2\ +
]щт Гп А Г Г]^ Аг \б 2 ) +
+ На(хп-1/2) Гп-1/2 (Р)к - (Р)кп-х / Ъп-1/2\ (18)
гп Аг г]кп-1/2 Аг V 6 2 )' ()
Ък+1 - К , 1 гп+1ккпХ\ (и)кпХ\ - гп-1ккп-\(и)к-\ __
А* + гп 2А г = ' ( )
и.\к _ и л(Ъп)2 (Р)п+1 -(Р) п—1 /оп\
(и)п = -НаЫ^Г-2Д-' (20)
( С)к+1 - (С)к + \(и)к+1| + (и)к+1 (С)кп - (С)к-1
Дt 2 Аг
\(и)кп+1\ - (и)к+1 (С)кп+1 - (С)к
2А
3п (С)п +1 На(хп) ( ък+1 (С)п+1 - (С
нк+1 + Ре Гплк+1А г V Гп+1/Фп+1/2 А г
- Гп-1^п-\/2 (С)п А (С)п-1 ), (21)
(Р )п = - са а^ I ^^—аг— ^р—дт-) > (22)
где хп = Са - (С)к - 5х и к^ = 10/6х = 2000. Промежуточные значения между узлами вида ¡п±1/2 рассчитываются как (¡п±1 + ¡п)/2. Стоит заметить, что в (21) для дискретизации конвективного слагаемого используется схема с разностями против потока. Для приближения пространственной производной первого порядка в (17), (19) и (20) используется центральная разностная схема. Это дает второй порядок аппроксимации по А г, как и аппрокисма-ция производной второго порядка в (18), (21) и (22). Здесь используется как явная разностная схема (21), так и неявная — (19).
Теперь запишем начальные и граничные условия из п. 4 в дискретном виде. Начальная форма свободной поверхности в дискретном виде задается как = 1 - г2. Формула начального распределения массовой доли частиц принимает вид
.о „ 2 -Со + 2(Со - 1)
(С )п = С9
1 + ехр(адс(гп - 1))'
С граничными условиями первого рода все просто, /кт = с, где £ — искомая функция и с — некоторая константа (здесь п = 0 или п = N). Для граничных условий второго рода вида д//дг = 0 используются аппроксимации второго порядка точности: /0,т = (4/1,т - ¡2,т)/3, ¡м,т = 4¡м-1,т - 3¡м-2,т. Для величин, не зависящих от х, индекс т в этих формулах отсутствует.
Дискретное уравнение (19) можно расщепить по физическим процессам на два уравнения. Одно из них описывает изменение Ъ в результате испарения
ъЬтр Ъ к
Ъп_'Ъп. = _ тк
А1 = °п,
в результате конвективного переноса массы: ЪП+1 — ЪОТ , 1 Гп+1ЪП%\(п)П%\ — Тп-1 Ъкп-\(и)кп-\
к
другое
А
+
2А
0.
(23)
(24)
Здесь (23) и (24)—явная и неявная разностные схемы соответственно. Временное значение высоты капли Ъ"9 определяется из (23). Затем, после пересчета давления (22) и усредненной скорости (20), решается (24) и находится итоговое значение 1 на новом временном слое. Уравнение (24) можно переписать в виде
апЪкп-\ + Ъп Ък+1 + СпЪкп\ \ = б,п, (25)
где ап = —А1гп-1(и)п-\/(2ГпАг), Ъп = 1, Сп = \(и)пХ\/(2гпАг) и бп =
= Ъ"9. В совокупности для всех узлов п из (25) получаем СЛАУ, которую можно представить в матричной форме:
о о 0 0
а1 Ь С1 0
0 а2 Ь2 С2
0 0 ам-
0 0 0
\
ам
Ъм
/Ък+1\ Ък+\
Ъм-1
( бо \
бм -1 V бм /
Запишем последнее выражением в более компактном виде:
а ■ ь = а,
где А — матрица коэффициентов, а — вектор правых частей и Ь — искомый вектор. Обратим внимание, что коэффициенты ао и см в вычислениях не задействованы, квадратная матрица А характеризуется ленточной трехдиаго-нальной структурой ненулевых элементов и такая СЛАУ может быть эффективно решена методом прогонки. На основании граничных условий запишем значения коэффициентов Ъо, со, бо, ам, Ъм и бм. Чтобы не нарушить трехдиа-гональную структуру матрицы А, для граничного условия на оси симметрии дЪ/дг = 0 используем аппроксимацию первого порядка (Ъ,{+1 — Ъ/к+1 )/Аг & 0. В таком случае получаем коэффициенты Ьо = 1, со = — 1 и бо = 0. Условие на трехфазной границе Ъ = Ъ/ приводит к коэффициентам ам = 0, Ъм = 1 и бм = Ъ^.
Значения Ъп+1, (и)^1 и (Р)к+1 необходимо итеративно уточнять. На первой итерации (г = 1) задаем отправное значение Ъп^1'1 = Ъп"'9. Относительно резкое изменение давления за дискретный временной шаг в результате испарения согласно уравнению (23) приводит к неустойчивости вычислений. По этой причине здесь используется метод явной релаксации
( Р )п+1'г = (Р)
к+1,%-1
+ а((Р )п+1'г — (Р)
\к+1,г-1
),
п
п
п
где г —номер итерации, а а — фактор релаксации. Предварительно (Р)^1,г рассчитывается по формуле (22). Значение фактора релаксации подобрано эмпирическим путем. Здесь в расчетах использовалось значение а = 10-3. Малое значение этого параметра улучшает устойчивость численного счета, но при этом увеличивает количество итераций, необходимых для сходимости решения. Вблизи границ при расчете усредненной скорости (и) по формуле (20) возникают осцилляции, которые за несколько временных шагов приводят к «разрушению» вычислений. С целью предотвращения таких ос-цилляций для приграничных узлов использовалась линейная интерполяция. Учитывая, что (и) ^ 0 при г ^ 0 или при г ^ К, получаем (и)1 = (и)2/2 и (и)м-1 = (и)м-2/2.
Механизм остановки итераций включает два возможных события: достигнуто максимальное число итераций 1тяж = 100 или максимальная невязка стала меньше заданной погрешности гтяж < ент, где е нт = 10-17 и
Утяж —
шах( ЪкХ1 - Ът +
п V
А гп+фкпХ\ (и)кХ\ - Гг-фЫ (и)кХ\ -
2А г
согласно уравнению (24).
После того как итерации завершились, полученные значения (и)кХ1
и (Р)п+1 используются для нахождения (С)к+1, и!+ и м^Х^ через явные зависимости (21), (17) и (18). Но предварительно необходимо подавить пилообразные осцилляции, возникающие при расчете давления. Этот вид осцил-ляций не приводит к «разрушению» численного счета со временем, но поле скорости потока жидкости получается некорректным. Чтобы разрешить эту проблему, сделаем следующее. Выражение (13) запишем в виде
д(Р) = -3 МЧ
дг Ъ2
и проинтегрируем его по :
(и)Г]
гг
(Р)\г - (Р)\о = -3
о
Ъ2
(г.
Слагаемое с интегралом аппроксимируем методом прямоугольников и получим формулу для восстановления давления:
(Р^ = -3А г £ Мт + ( РуоРЯк. (26)
п=0
Учитывая пилообразную форму осцилляций, (р у0ера1г можно аппроксимировать следующим образом:
(Р« ((Р)1 + (Р)2)/2.
Со временем такие пилообразные осцилляция могут передаваться усредненной скорости (и) и профилю капли Ъ. Здесь эта проблема решена за счет использования сглаживающего фильтра:
/Г = Р к+1 +7 к + 5 ^-1. (27)
п
Параметры подобраны эмпирическим путем: р = 5 = 0.25 и 7 = 0.5. В отличие от [49] используемый здесь фильтр прост в реализации, занимает в разы меньше памяти и вычислительных ресурсов. Кроме того, количество итераций фильтра из [49] заранее неизвестно. В нашем случае достаточно одного прохода по узлам п = 1, 2,..., N — 1. Узлы п = 0 и п = N рассчитываются из граничных условий.
Блок-схема с алгоритмом1 расчета представлена на рис. 2. В ней используются следующие дополнительные обозначения: рс — счетчик временных периодов, рп — количество временных периодов, р3 — размер временного периода во временных шагах и Нп = На(хп). Записи вида = следует трактовать как векторные операции.
Begin J^)
Pc = Pc + 1
Writing data to a file
Calculation ui+i and «¿+1
Problem parameters and initial conditions
Pc = 0, p„ = 45, ps = tf/(Atpn)
<
: = 1> P.i 1
>
Calculation himp, Jt and Hn
Filtering (u) and h
Recovery {.P) from (u)
С
Writing calculation parameters to a file
End
(C)kn = (С)к„+\
ftj = hk„+1
Calculation
(C) *+l n
i = ^max 1, Elim + 1
Calculation rmax, i = i + 1
SLAE solution A • h = d
Calculation Calculation / \ft+i
Рис. 2. Алгоритм численного расчета [Figure 2. Numerical calculation algorithm]
Для сравнения численный расчет также выполнялся в коммерческом пакете Maple 18 с помощью модуля pdsolve с параметром method = Theta[l] (неявная разностная схема). Для того чтобы pdsolve справился с задачей, систему уравнений пришлось переписать в другом виде. Вместо (13) и (14) использовались уравнения
ib- 1 dh Ca dr '
(28)
хКод программы на языке С++ доступен по ссылке https://github.com/kolegovk/ Suppression-of-sawtooth-oscillations.git
где ф — вспомогательная функция. Причем модуль pdsolve не требует граничных условий для (28). Скорее всего, модуль определяет их автоматически из заданных условий для h. Для уравнения (29) на границах задается равенство скорости нулю (см. п. 4). Таким образом, в Maple решалась система уравнений (12), (16), (28) и (29).
6. Результаты и обсуждение. Результаты расчетов средствами двух программ сравниваются на рис. 3. На большинстве графиков расхождение практически незаметно. Максимальное отличие результата, полученного с помощью программы, написанной на С++, от результата, полученного с помощью модуля pdsolve в Maple 18, наблюдается для (и) на времени t = 220 с (рис. 3, c). В абсолютных величинах наибольшее расхождение соответствует точке г ~ 0.75 мм и составляет примерно 0.046 мкм/с, в относительных величинах это около 10 %.
Расчет выполнялся для N = 75. Значения временных шагов для разработанной программы и для модуля pdsolve брались разные, так как в них используются разные алгоритмы расчета, соответственно и требования к временным шагам накладываются тоже разные. В Maple 18 для численного решения СЛАУ используется итерационный метод Ньютона, но технические
0.70.60.50.40.3 0.2 0.10.0-
0.4 0.6
г, mm
b
0.4 0.6
С
Рис. 3. Сравнение результатов расчетов в двух программах для нескольких последовательных моментов времени: (а) высота капли, (b) массовая доля растворенного или взвешенного вещества и (с) радиальная скорость потока, усредненная по толщине жидкого слоя
[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]
a
детали реализации закрыты для пользователя, есть лишь краткая справка по разностным схемам. Скорость вычислений с относительно малым временным шагом, конечно же, ниже, но это частично компенсируется тем, что C++ является компилируемым языком, а язык Maple — интерпретируемый. Для расчетов в Maple использовалось значение At = 0.1 с, в разработанной программе— At = 10-6 с. Кроме того, для программы на C++ выполнены тестовые расчеты для N = 25 и N = 50 с целью проверки сходимости по сетке (во всех расчетах M = N). Для N = 25 приближенно подобрано максимальное значение временного шага At = 10-4 с, при котором расчет не «разрушается». Значению N = 50 соответствует At = 5 ■ 10-6 с. Таким образом, небольшое уменьшение Ar = R/N приводит к тому, что необходимо значительно уменьшать временной шаг. Эти наблюдения в вычислительных экспериментах косвенно свидетельствуют о том, что предложенная разностная схема является условно устойчивой. Длительность расчета в Maple составила около 10 мин для N = 75 (объем занимаемой памяти ~ 217 Мб). Время расчета в программе на С++ примерно составило 10 мин для N = 25, 5 ч для N = 50 и 29 ч для N =75 (объем занимаемой памяти 3 Мб). Расчеты проводились на компьютере с ЦП Intel Xeon CPU E3-1230 v3 c максимальной частотой 3.6 ГГц под управлением ОС Windows 10. Для получения исполняемого файла программы на С++ использовался компилятор, встроенный в Visual Studio 2022 Community Edition.
Проверка закона сохранения массы для растворенного или взвешенного вещества для N = 75 показала погрешность около 0.27 %, что экспериментально подтверждает свойство консервативности в предложенной разностной схеме. Для этого были рассчитаны и сопоставлены значения
N
'У ^ К(С)п
п=0
на моментах времени t = 0 и t = tf. Для расчетов в Maple также проводилась проверка закона сохранения массы. Выполнены тестовые расчеты и определены погрешности для разных временных шагов в диапазоне, для которого численное решение итерационно сходится: At = 0.01 с —0.43%, At = 0.1 с — 0.1 %, At = 1 с — 3.11 % и для At = 10 с — 32.54 %. Минимальная погрешность в 0.1 % соответствует временному шагу At = 0.1 с, который и был выбран для расчетов в Maple.
Толщина жидкого слоя с течением времени уменьшается (рис. 3, a). Из-за пространственно неравномерного испарения происходит отклонение формы капли от равновесного состояния, возникает градиент капиллярного давления, что порождает капиллярный поток (рис. 3, с), переносящий вещество к периферии капли. В связи с этим массовая доля вещества, содержащегося в жидкости, увеличивается вблизи трехфазной границы со временем (рис. 3, b). Под конец процесса массовая доля {С) распределена равномерно, так как достигнуто критическое значение Са, при котором произошел фазовый переход. Дальнейшее испарение жидкости из твердой фазы здесь не рассматривается.
Финальная расчетная толщина слоя h является почти равномерной в пространстве. Вблизи г ~ 0.9 мм наблюдается незначительное возвышение, которое связано с эффектом кофейных колец [3]. Возникающий под конец про-
цесса обратный поток жидкости, направленный в центральную область, смазывает этот эффект. Это можно использовать для улучшения существующих методов подавления эффекта кофейных колец (см. ссылки в [1]).
Также стоит заметить, что поначалу усредненная радиальная скорость потока с течением времени увеличивается в связи с увеличением градиента капиллярного давления в процессе неравномерого испарения, преобладающего вблизи периферии. Затем в области роста массовой доли вещества, приводящего к увеличению вязкости, поток начинает замедляться. Вблизи трехфазной границы появляется противоток (рис. 3, с для £ = 220 с), что ранее наблюдалось в эксперименте [52]. Далее по времени направление потока постепенно меняется во всей области (рис. 4, а). Это связано с тем, что в начале процесса плотность потока пара преобладает в районе периферии, но со временем .] начинает увеличиваться в центральной области из-за понижения Ъ и уменьшаться вблизи контактной линии из-за роста (С) (рис. 4, Ь). Когда во всей области достигается значение (С) ~ Сд, модель предсказывает отсутствие потока и испарения. Поле скорости потока для различных моментов времени представлено на рис. 5. В отличие от результатов расчетов [53] здесь на поздних временных этапах наблюдается изменение направления потока. Это связано с тем, что в настоящей работе учитывается зависимость .] и г] от массовой доли (С). Кроме этого, здесь при заданных параметрах градиент плотности потока пара с течением времени меняет знак. Технические детали численного решения в [53] не раскрыты. В [27,28,49] поле скорости потока не рассчитывалось.
Здесь поле скорости потока удалось рассчитать по формулам (17) и (18) лишь после того, как была решена проблема с пилообразными осцилляциями, возникающими при расчете капиллярного давления и постепенно с временными шагами накапливающимися в (и) и Ъ. Для этого использовались формулы (26) и (27). Восстановленное из скорости потока (и) по формуле (26)
Рис. 4. Результаты расчета в Maple: пространственно-временная зависимость (а) усредненной по толщине жидкого слоя скорости радиального потока и (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]
a
b
", рш/в
V, pm/s
0.4 0.6
a
v, pm/s
0.4 0.6
b
v, pm/s
0.4 0.6 0.8
Рис. 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]
27.89527.89027.88527.880-• 27.87511 27.87027.86527.86027.855-
0.2 0.4 0.6 0.8
a 20-fu
a,"
0.0 0.2 0.4 0.6
T, ШШ b
Рис. 6. Результаты расчета в программе, написанной на С++: (a) капиллярное и (b) возможные пилообразные осцилляции
[Figure 6. The results of the calculation in program written with C++: (a) capillary and (b) possible sawtooth oscillations]
давление pressure
d
c
a
капиллярное давление представлено на рис. 6, а. Если такое восстановление не проводить, то получаем осцилляции, представленные на рис. 6, Ь. Особенность этих осцилляций заключается в стабильности на временных шагах, то есть они не приводят к полному «разрушению» расчета. Пилообразные осцилляции теоретически изучаются в [54], но здесь предложено практическое решение проблемы на конкретном примере.
Заключение. Для решения практических (прикладных) задач в области испарительной самосборки и испарительной литографии требуется разработка математических моделей, численных методов и комплекса программ. При разработке необходимо учитывать особенности таких задач, ведь сопутствующие процессы весьма сложны. Вязкость раствора, градиент плотности потока пара и локальная кривизна свободной поверхности влияют на структуру потока в сидячей капле. В работе показан случай, когда направление потока жидкости может меняться со временем в процессе испарения даже в отсутствие эффекта Марангони и какого-либо внешнего воздействия на систему. Это может приводить к замедлению выноса вещества на периферию, что в результате будет способствовать формированию более или менее равномерного осадка по всей площади контакта капли с подложкой. Данное наблюдение полезно для совершенствования методов подавления кольцевых осадков, связанных с эффектом кофейных колец и нежелательных для некоторых приложений, как, например, струйная печать или нанесение покрытий. Моделирование гидродинамических потоков важно, так как это один из основных механизмов, влияющих на массоперенос растворенного или взвешенного вещества. При численной реализации подобных задач может появиться проблема, связанная с возникновением пилообразных осцилляций, устойчивых на временных шагах, но создающих проблему в расчете поля скорости потока жидкости. Здесь предложены практические рецепты, позволяющие эффективно бороться с такими осцилляциями при разработке программного обеспечения для вычислительной гидродинамики.
Конкурирующие интересы. Заявляю, что в отношении авторства и публикации этой статьи конфликта интересов не имею.
Авторская ответственность. Я несу полную ответственность за предоставление окончательной версии рукописи в печать. Окончательная версия рукописи мною одобрена.
Финансирование. Часть работы, относящаяся к разработке математической модели, выполнена при поддержке программы развития университета «Приоритет 2030» проект № 122112500011-0 «Природовдохновленные оптические технологии» (Центр природовдохновленного инжиниринга ТюмГУ). Исследование, связанное с разработкой численного метода и программы для ЭВМ, выполнено за счет гранта Российского научного фонда № 22-79-10216 в АГУ им. В. Н. Татищева, Ь^р8://гБс^ ги/рго^ес-Ъ/22-79-10216/.
Благодарность. Автор выражает благодарность профессору, доктору физико-математических наук Ю. Ю. Тарасевичу за полезные рекомендации по оформлению графического материала.
Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki
[J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2023, vol. 27, no. 2, pp. 309-335
d https://doi.org/10.14498/vsgtu1994
ISSN: 2310-7081 (online), 1991-8615 (print)
MSC: 65M06, 80A20
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
K. S. Kolegov1'2
1 Astrakhan Tatishchev State University,
20 a, Tatischev st., Astrakhan, 414056, Russian Federation.
2 Tyumen State University,
25, Lenin st., Tyumen, 625003, Russian Federation.
Abstract
Evaporating droplets and films are used in applications from different fields. Various methods of evaporative self-assembly are of particular interest. The paper describes a mathematical model of mass transfer in a droplet drying on a substrate based on the lubrication approximation. The model takes into account the transfer of a dissolved or suspended substance by a capillary flow, the diffusion of this substance, the evaporation of liquid, the formation of solid deposit, the dependence of the viscosity and the vapor flux density on the admixture concentration.
The case with pinning of the three-phase boundary ("liquid-substrate-air") is considered here. Explicit and implicit finite-difference schemes have been developed for the model equations. A modification of the numerical method is proposed, in which splitting by physical processes, the iterative method of explicit relaxation and Thomas algorithm are combined. A practical recipe for suppressing sawtooth oscillations is described using the example of a specific problem.
A software module in C++ has been developed, which can be used for evaporative lithography problems in the future. With the help of this module, numerical calculations were carried out, the results of which were compared with the results obtained in the Maple package.
Numerical simulation predicted the case in which the direction of the capillary flow changes to the opposite over time due to a change in the sign
Mathematical Modeling, Numerical Methods and Software Complexes Research Article
© Authors, 2023
© Samara State Technical University, 2023 (Compilation, Design, and Layout)
3 ©® The content is published under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/) Please cite this article in press as:
Kolegov K. S. 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, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2023, vol. 27, no. 2, pp. 309-335. EDN: HZRAJS. DOI: 10.14498/vsgtu1994 (In Russian).
Author's Details:
Konstantin S. Kolegov © https://orcid.org/0000-0002-9742-1308
Cand. Phys. & Math. Sci.; Senior Researcher; Lab. of Mathematical Modeling and Information Technologies in Science and Education1 ; Center for Nature-inspired Engineering2; e-mail: [email protected]
of the gradient of the vapor flux density. This can lead to a slowdown in the transfer of the substance to the periphery, which as a result will contribute to the formation of a more or less uniform precipitation over the entire contact area of the droplet with the substrate. This observation is useful for improving methods of annular deposit suppression associated with the coffee-ring effect and undesirable for some applications, such as inkjet printing or coating.
Keywords: evaporating droplet, mass transfer, capillary flow, finite-difference scheme, sawtooth oscillation, suppression of the coffee-ring effect.
Received: 18th January, 2023 / Revised: 29th April, 2023 / Accepted: 3rd May, 2023 / First online: 27th June, 2023
Competing interests. I declare that I have no conflicts of interest regarding the authorship and publication of this article.
Author's Responsibilities. I take full responsibility for submitting the final version of the manuscript for printing. The final version of the manuscript has been approved by me.
Funding. Part of the work related to the development of the mathematical model was carried out with the support of the University Development Program "Priority 2030", project no. 122112500011-0 "Nature-inspired optical technologies" (Center for Nature-inspired Engineering, Tyumen State University). The research related to the development of the numerical method and computer program was supported by the grant from the Russian Science Foundation, project no. 22-79-10216 at Astrakhan Tatishchev State University, https://rscf.ru/en/project/22-79-10216/.
Acknowledgments. The author is grateful to Prof. Yu. Yu. Tarasevich, Dr. Phys. & Math. Sci., for valuable recommendations regarding the presentation of graphic material.
References
1. Zang D., Tarafdar S., Tarasevich Yu. Yu., et al. Evaporation of a droplet: From physics to applications, Phys. Rep., 2019, vol.804, pp. 1-56. EDN: HMEGXP. DOI: https://doi.org/ 10.1016/j.physrep.2019.01.008.
2. Kolegov K. S., Barash L. Yu. Applying droplets and films in evaporative lithography, Adv. Colloid Interf. Sci., 2020, vol. 285, 102271. EDN: BEBGWV DOI: https://doi.org/10.1016/j-cis.2020.102271.
3. Deegan R. D., Bakajin O., Dupont T F., et al. Capillary flow as the cause of ring stains from dried liquid drops, Nature, 1997, vol.389, no. 6653, pp. 827-829. DOI: https://doi. org/10.1038/39827.
4. Kim H., Yang J., Kim S., et al. Numerical simulation of the coffee-ring effect inside containers with time-dependent evaporation rate, Theor. Comput. Fluid Dyn., 2022, vol.36, no. 3, pp. 423-433. DOI: https://doi.org/10.1007/s00162-021-00602-x.
5. Baba H., Yoshioka R., Takatori S., et al. Transitions among cracking, peeling and homoge-nization on drying of an aqueous solution containing glucose and starch, Chem. Let., 2021, vol.50, no. 5, pp. 1011-1014. DOI: https://doi.org/10.1246/cl.210009.
6. Wang W., Wang Q., Zhang K., et al. On-demand contact line pinning during droplet evaporation, Sens. Act. B: Chem., 2020, vol.312, 127983. DOI: https://doi.org/10.1016/j. snb.2020.127983.
7. Saroj S. K., Panigrahi P. K. Magnetophoretic control of diamagnetic particles inside an evaporating droplet, Langmuir, 2021, vol.37, no. 51, pp. 14950-14967. DOI: https://doi. org/10.1021/acs.langmuir.1c02968.
8. Al-Muzaiqer M. A., Kolegov K. S., Ivanova N. A., Fliagin V. M. Nonuniform heating of a substrate in evaporative lithography, Phys. Fluids, 2021, vol. 33, no. 9, 092101. DOI: https:// doi.org/10.1063/5.0061713.
9. Du F., Zhang L., Shen W. Controllable dried patterns of colloidal drops, J. Col. Int. Sci., 2022, vol.606, pp. 758-767. DOI: https://doi.Org/10.1016/j.jcis.2021.08.089.
10. Cedeno R., Grossier R., Lagaize M., et al. Nucleation in sessile saline microdroplets: Induction time measurement via deliquescence-recrystallization cycling, Faraday Discuss., 2022, vol.235, pp. 183-197. DOI: https://doi.org/10.1039/d1fd00090j.
11. Perkins-Howard B., Walker A. R., Do Q., et al. Surface wettability drives the crystalline surface assembly of monodisperse spheres in evaporative colloidal lithography, J. Phys. Chem. C, 2022, vol.126, no. 1, pp. 505-516. DOI: https://doi.org/10.1021/acs.jpcc. 1c07098.
12. Jose M., Mayarani M., Basavaraj M. G., Satapathy D. K. Evaporative self-assembly of the binary mixture of soft colloids, Phys. Chem. Chem. Phys., 2021, vol. 23, no. 12, pp. 71157124. DOI: https://doi.org/10.1039/d1cp00440a.
13. Zolotarev P. A., Kolegov K. S. Monte Carlo simulation of particle size separation in evaporating bi-dispersed colloidal droplets on hydrophilic substrates, Phys. Fluids, 2022, vol. 34, no. 1, 017107. DOI: https://doi.org/10.1063/5.0072083.
14. Kirner F., Sturm E. V. Advances of nonclassical crystallization toward self-purification of precious metal nanoparticle mixtures, Cryst. Growth Des., 2021, vol.21, no. 9, pp. 51925197. DOI: https://doi.org/10.1021/acs.cgd.1c00544.
15. Hossain M. T., Gates I. D., Natale G. Dynamics of Brownian Janus rods at a liquidliquid interface, Phys. Fluids, 2022, vol. 34, no. 1, 012117. DOI: https://doi.org/10.1063/ 5.0076148.
16. Mustakim M., Anil Kumar A. V. Depletion induced demixing and crystallization in binary colloids subjected to an external potential barrier, J. Phys. Chem. B, 2021, vol. 126, no. 1, pp. 327-335. DOI: https://doi.org/10.1021/acs.jpcb.1c08591.
17. Nozawa J., Uda S., Toyotama A., et al. Heteroepitaxial fabrication of binary colloidal crystals by a balance of interparticle interaction and lattice spacing, J. Col. Int. Sci., 2022, vol.608, pp. 873-881. DOI: https://doi.org/10.1016Zj.jcis.2021.10.041.
18. Galy P. E., Guitton-Spassky T., Sella C., et al. Redox control of particle deposition from drying drops, ACS Appl. Mater. Interfaces, 2022, vol. 14, no. 2, pp. 3374-3384. DOI: https:// doi.org/10.1021/acsami.1c18933.
19. Inoue K., Inasawa S. Drying-induced back flow of colloidal suspensions confined in thin unidirectional drying cells, RSC Adv. ,2020, vol. 10, no. 27, pp. 15763-15768. DOI: https:// doi.org/10.1039/d0ra02837a.
20. Homede E., Manor O. Deposition of nanoparticles from a volatile carrier liquid, J. Col. Int. Sci., 2020, vol. 562, pp. 102-111. DOI: https://doi.org/10.1016/jjcis.2019.11.062.
21. Li D., Chen R., Zhu X., et al. Light-fueled beating coffee-ring deposition for droplet evaporative crystallization, Anal. Chem., 2021, vol. 93, no. 25, pp. 8817-8825. DOI: https://doi. org/10.1021/acs.analchem.1c00605.
22. Shao X., Hou Y., Zhong X. Intense jet flow with symmetric vortices induced by saline concentration gradient at free surface of a drying saline droplet, Int. Comm. Heat Mass Transf., 2021, vol.128, 105573. DOI: https://doi.org/10.1016/j.icheatmasstransfer. 2021.105573.
23. Hegde O., Chatterjee R., Rasheed A., et al. Multiscale vapor-mediated dendritic pattern formation and bacterial aggregation in complex respiratory biofluid droplets, J. Col. Int. Sci., 2022, vol. 606, pp. 2011-2023. DOI: https://doi.org/10.1016/j.jcis.2021.09.158.
24. Corletto A., Shapter J. G. Thickness/morphology of functional material patterned by topographical discontinuous dewetting, Nano Select, 2021, no. 9, pp. 1723-1740. DOI: https:// doi.org/10.1002/nano.202000301.
25. Bayat F., Tajalli H. Nanosphere lithography: the effect of chemical etching and annealing sequence on the shape and spectrum of nano-metal arrays, Heliyon, 2020, vol. 6, no. 2, e03382. DOI: https://doi.org/10.1016/j.heliyon.2020.e03382.
26. Bhardwaj R., Fang X., Somasundaran P., Attinger D. Self-assembly of colloidal particles from evaporating droplets: role of DLVO interactions and proposition of a phase diagram, Langmuir, 2010, vol. 26, no. 11, pp. 7833-7842. DOI: https://doi.org/10.1021/la9047227.
27. Tarasevich Yu. Yu., Vodolazskaya I. V., Isakova O. P. Desiccating colloidal sessile drop: Dynamics of shape and concentration, Coll. Pol. Sci., 2011, vol.289, no. 9, pp. 1015-1023. EDN: OHRMJP. DOI: https://doi.org/10.1007/s00396-011-2418-8.
28. Kolegov K. S., Lobanov A. I. Numerical study of mass transfer in drop and film systems using a regularized finite difference scheme in evaporative lithography, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2018, vol.22, no. 2, pp. 344-363 (In Russian). EDN: XWXSMX. DOI: https://doi.org/ 10.14498/vsgtu1601.
29. Petsi A. J., Kalarakis A. N., Burganos V. N. Deposition of Brownian particles during evaporation of two-dimensional sessile droplets, Chem. Eng. Sci., 2010, vol. 65, no. 10, pp. 29782989. DOI: https://doi.org/10.1016Zj.ces.2010.01.028.
30. Hu S., Wang Y., Man X., Doi M. Deposition patterns of two neighboring droplets: Onsager variational principle studies, Langmuir, 2017, vol.33, no. 23, pp. 5965-5972. DOI: https:// doi.org/10.1021/acs.langmuir.7b01354.
31. Yang J., Kim H., Lee C., et al. Phase-field modeling and computer simulation of the coffee-ring effect, Theor. Comput. Fluid Dyn., 2020, vol.34, no. 5-6, pp. 679-692. DOI: https:// doi.org/10.1007/s00162-020-00544-w.
32. Seo H. W., Jung N., Yoo C. S. Oscillation dynamics of colloidal particles caused by surfactant in an evaporating droplet, J. Mech. Sci. Technol., 2020, vol.34, no. 2, pp. 801-808. DOI: https://doi.org/10.1007/s12206-020-0128-1.
33. Kim H.-S., Park S. S., Hagelberg F. Computational approach to drying a nanoparticle-suspended liquid droplet, J. Nanopart. Res., 2010, vol. 13, no. 1, pp. 59-68. DOI: https:// doi.org/10.1007/s11051-010-0062-8.
34. Crivoi A., Duan F. Three-dimensional Monte Carlo model of the coffee-ring effect in evaporating colloidal droplets, Sci. Rep., 2014, vol.4, no. 1, 4310. DOI: https://doi.org/ 10.1038/srep04310.
35. Zhang H., Shan Y. G., Li L., et al. Modeling the self-assembly of nanoparticles into branched aggregates from a sessile nanofluid droplet, Appl. Therm. Eng., 2016, vol. 94, pp. 650-656. DOI: https://doi.org/10.1016/j.applthermaleng.2015.10.160.
36. Ren J., Crivoi A., Duan Fe. Disk-ring deposition in drying a sessile nanofluid droplet with enhanced Marangoni effect and particle surface adsorption, Langmuir, 2020, vol. 36, no. 49, pp. 15064-15074. DOI: https://doi.org/10.1021/acs.langmuir.0c02607.
37. Ren J., Crivoi A., Duan F. Dendritic nanoparticle self-assembly from drying a sessile nanofluid droplet, Phys. Chem. Chem. Phys., 2021, vol.23, no. 29, pp. 15774-15783. DOI: https://doi.org/10.1039/d1cp01181b.
38. Lebedev-Stepanov P., Vlasov K. Simulation of self-assembly in an evaporating droplet of colloidal solution by dissipative particle dynamics, Colloids Surf. A Physicochem. Eng. Asp., 2013, vol.432, pp. 132-138. EDN: RFHBFJ. DOI: https://doi.org/10.1016/j.colsurfa. 2013.05.012.
39. Breinlinger T., Kraft T. A simple method for simulating the coffee stain effect, Powder Technol., 2014, vol.256, pp. 279-284. DOI: https://doi.org/10.1016/j.powtec.2014.02. 024.
40. Liu W., Midya J., Kappl M., et al. Segregation in drying binary colloidal droplets, ACS Nano, 2019, vol.13, no. 5, pp. 4972-4979. DOI: https://doi.org/10.1021/acsnano. 9b00459.
41. Andac T., Weigmann P., Velu S. K. P., et al. Active matter alters the growth dynamics of coffee rings, Soft Matter, 2019, vol.15, no. 7, pp. 1488-1496. DOI: https://doi.org/ 10.1039/c8sm01350k.
42. Kolegov K. S. Monte Carlo simulation of colloidal particles dynamics in a drying drop, J. Phys.: Conf. Ser., 2019, vol. 1163, 012043. EDN: MSBXPO. DOI: https://doi.org/10.1088/ 1742-6596/1163/1/012043.
43. Lebovka N. I., Tarasevich Yu. Yu., Bulavin L. A., et al. Sedimentation of a suspension of rods: Monte Carlo simulation of a continuous two-dimensional problem, Phys. Rev. E, 2019, vol.99, no. 5, 052135. EDN: XNZLXM. DOI: https://doi.org/10.1103/physreve.99.052135.
44. Kolegov K. S., Barash L. Yu. Joint effect of advection, diffusion, and capillary attraction on the spatial structure of particle depositions from evaporating droplets, Phys. Rev. E, 2019, vol. 100, no. 3, 033304. EDN: RUWORY. DOI: https://doi.org/10.1103/physreve.100. 033304.
45. Zolotarev P. A., Kolegov K. S. Average cluster size inside sediment left after droplet desiccation, J. Phys.: Conf. Ser., 2021, vol.1740, 012029. EDN: OPMGGJ. DOI: https://doi.org/ 10.1088/1742-6596/1740/1/012029.
46. Song R., Lee M., Moon H., et al. Particle dynamics in drying colloidal solution using discrete particle method, Flex. Print. Electron., 2021, vol.6, no.4, 044007. DOI: https://doi.org/ 10.1088/2058-8585/ac428e.
47. Marinaro G., Riekel C., Gentile F. Size-exclusion particle separation driven by micro-flows in a quasi-spherical droplet: Modelling and experimental results, Micromachines, 2021, vol. 12, no. 2, 185. DOI: https://doi.org/10.3390/mi12020185.
48. Hu H., Larson R. G. Analysis of the microfluid flow in an evaporating sessile droplet, Langmuir, 2005, vol.21, no. 9, pp. 3963-3971. DOI: https://doi.org/10.1021/la047528s.
49. Park Y., Park Y, Lee J., Lee C. Simulation for forming uniform inkjet-printed quantum dot layer, J. Appl. Phys., 2019, vol.125, no. 6, 065304. DOI: https://doi.org/10.1063/1. 5079863.
50. Kolegov K. S. Simulation of mass transfer in droplet-film systems using a regularized difference scheme in evaporative lithography, Thesis of Dissertation (Cand. Phys. & Math. Sci.). Astrakhan, Astrakhan State University, 2018, 162 pp. (In Russian). DOI: https://doi.org/ 10.13140/RG.2.2.29663.97446
51. Kolegov K. S. Simulation of patterned glass film formation in the evaporating colloidal liquid under IR heating, Microgravity Sci. Technol., 2018, vol.30, no. 1-2, pp. 113-120. EDN: UXPYTS. DOI: https://doi.org/10.1007/s12217-017-9587-0.
52. Bodiguel H., Leng J. Imaging the drying of a colloidal suspension, Soft Matter, 2010, vol. 6, pp. 5451-5460. DOI: https://doi.org/10.1039/C0SM00323A.
53. Fischer B. J. Particle convection in an evaporating colloidal droplet, Langmuir, 2002, vol. 18, no. 1, pp. 60-67. DOI: https://doi.org/10.1021/la015518a.
54. Dorodnitsyn L. V. Grid oscillations in finite-difference scheme and a method for their approximate analysis, Comput. Math. Model., 2016, vol.27, no.4, pp. 472-488. EDN: YUUJUV. DOI: https://doi.org/10.1007/s10598-016-9337-y.
