Введение

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

В последние десятилетия интенсивно развиваются технологии, связанные с управлением распределением наночастиц и микроскопических объектов в прозрачных средах. Управляемое разделение частиц в наножидкостях находит применение в самых разных областях – от медицины [1] до космоса [2]. В обзоре [3] приведен сравнительный анализ современных методов, которые обеспечивают неразрушающее манипулирование частицами для аналитических и биотехнологических применений. Среди них – оптический пинцет [4], а также методы магнитной и термической сепарации [5; 6].

Одной из ключевых групп методов манипулирования частицами является использование оптических методов, среди которых выделяется лазерное облучение. Примером успешного внедрения метода являются работы по созданию оптических ловушек для эффективной сортировки и захвата частиц определенных размеров в зависимости от плотности мощности излучения [7; 8].

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

Рассматриваемая математическая начально-краевая задача имеет сингулярность решения, обусловленную малым параметром при старшей производной. Класс задач с сингулярностью исследуется в работах [9–11]. V. A. Rukavishnikov и E. I. Rukavishnikova для класса задач с вырождением на всей границе области получили оценки скорости сходимости, что непосредственно затрагивает рассматриваемую проблему [11].

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

Материалы и методы исследования

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

Рис. 1. Схема кюветы с жидкостью и падающим лазерным излучением. Примечание: составлен авторами на основе источника [12]

Физическая постановка предполагает одномерное распределение концентрации n наночастиц вдоль вертикальной оси кюветы, что позволяет свести процесс к математической модели, описываемой начально-краевой задачей для параболического уравнения непрерывности, являющегося следствием уравнения Эйштейна – Фоккера – Планка, сформулированной в [12; 13]:

с граничными условиями

при или

и начальным условием

где n – концентрация наночастиц; и – пространственная координата вдоль оси цилиндра и время; D и v – коэффициент диффузии и скорость систематического смещения частиц под действием излучения соответственно; l – высота кюветы.

Введем безразмерные величины

Неизвестной функцией является концентрация частиц n(z,t). В новых переменных исходная задача [12] формулируется как начально-краевая задача для уравнения

(1)

с краевыми условиями нулевого потока

, при z = 0 или z = 1

и начальным условием однородного распределения

Точное аналитическое решение задачи (1) получено в [12; 13] и представлено в виде бесконечного сходящегося ряда

(2)

где .

Для удобства дальнейшего анализа введем нормированную безразмерную функцию

.

При малых значениях безразмерного параметра скорости дрейфа частиц (u ≤ 10) ряд (2) сходится быстро, и его численное суммирование не представляет вычислительных сложностей [14]. Однако практический интерес представляют режимы интенсивного осаждения под действием лазерного излучения, соответствующие u > 100. В этом случае в уравнении (1) при производной по времени и старшей производной по координате появляется малый параметр 1 / u, что придает задаче сингулярно возмущенный характер. Краевые задачи для уравнений сингулярного типа исследовались в [9], а начально-краевые задачи конвекции-диффузии – в [10].

Для расчета при больших значениях u (или, эквивалентно, малых 1 / u) в данной работе предложено два взаимодополняющих подхода:

1) аналитический подход: суммирование ряда (2) с использованием большого числа слагаемых (~105) и арифметики произвольной точности, реализуемой в системе компьютерной алгебры Maple. Заметим, что для z = 1/2 суммирование можно выполнить для «бесконечного» (infinity в Maple) числа слагаемых;

2) численный подход: численное решение начально-краевой задачи (1) методом прямых в среде численного анализа MATLAB.

В работах [14; 15] проанализировано качественное поведение решения (2) и выявлены основные вычислительные трудности при u > 40, обусловленные двумя факторами, а именно:

1) потерей значащих цифр при суммировании быстро убывающих или быстро растущих экспоненциальных слагаемых;

2) ограничением на максимальное представимое число в формате с плавающей запятой двойной точности (≈1,8·10308) [16].

Как показано в [14; 15], применение специализированных алгоритмов позволяет преодолевать эти ограничения и проводить расчеты, например, для u = 2000, а при z = 1/2 для u = 70300, однако получаемые значения функции F(z,t) остаются доступными лишь для фиксированных координат z, что затрудняет построение полных пространственно-временных зависимостей.

В настоящей работе разработан адаптивный алгоритм для MATLAB, который автоматически выбирает стратегию вычисления в зависимости от величины параметра u. Данный подход позволяет строить трехмерные графики функции F(z,t) в широком диапазоне параметров. При проведении исследований использована версия MATLAB R2021a.

Для численного решения начально-краевой задачи (1) реализована вычислительная процедура, в которой используется встроенная функция pdepe(). Формат вызова функции pdepe() следующий:

pdepe(0, @pdefun, @icfun, @bcfun, zspan, tspan),

где zspan и tspan задают пространственно-временную сетку, а указатели @pdefun, @icfun, @bcfun ссылаются на функции, описывающие соответственно дифференциальное уравнение (1), начальное и краевые условия. После получения решения от pdepe() проводится его постобработка и визуализация результатов вычислений.

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

– пространственная сетка zspan по безразмерному параметру z ∈ [0;1]:

,

где Nz = 1800;

– сетка по времени tspan строится в интервале от 0 до tend:

,

где Nt = 18000.

Заметим, что решатель ode15s является адаптивным, то есть может изменять шаг по времени для обеспечения устойчивости в заданных пределах.

Вычисления проведены с заданными по умолчанию параметрами options для решателя pdepe, в том числе использованы следующие значения для допусков: RelTol = 10-3, AbsTol = 10-6.

Для проверки численного метода проведено сравнение с известными значениями аналитического решения при u = 17,2. Также в работе представлены результаты расчетов для практически значимых значений параметра u = 8770 и u = 70300, соответствующих условиям реального эксперимента. Заметим, что значение безразмерного параметра u = 17,2 относится к области малых скоростей дрейфа наночастиц, а u = 8770 и u = 70300 являются высокими и даже экстремальными величинами. Для таких скоростей дрейфа характерно быстрое затухание решения, то есть при времени t → tend функция F(z,t) экспоненциально убывает к нулю.

Время расчета tend определяется моментом выхода решения на стационарный режим:

Далее обозначено: Fmat – численное решение; Fser – аналитическое решение (2), методы вычисления которого изложены в [14; 15].

Сравнение решений Fmat и Fser проводится на отрезке z ∈ [0;8/9]. Численное решение Fmat получено для всей области определения z ∈ [0;1]; однако верификация численного метода выполнялась в указанной подобласти, где эталонный ряд (2) вычисляется прямым суммированием при N = 105. Верхняя граница по координате выбрана как z = 8/9, так как при экстремальных значениях параметра дрейфа точное аналитическое решение вблизи границы z = 1 становится сингулярным для фиксированного числа слагаемых ряда из-за катастрофической потери значащих цифр при суммировании экспонент. Ограничение подобластью позволяет вычислить значение ряда (2) прямым суммированием без привлечения процедур экстраполяции.

Для оценки точности численного решения использовались три величины:

1) разность решений Fmat – Fser;

2) относительная погрешность, вычисляемая в каждом узле сетки традиционным способом:

(3)

3) гибридная погрешность, которая отличается от (3) наличием регуляризатора в знаменателе.

Выбор метода вычисления погрешности зависел от поведения решения на заданных интервалах. При u = 17,2 функция F(z,t) не стремится экспоненциально к нулю, поэтому может быть использована формула (3).

Для режимов с высокими значениями параметра дрейфа формула (3) неприменима из-за деления на близкие к нулю значения функции F(z,t). В этом случае используется формула вычисления гибридной погрешности, аналогичная применяемой для контроля шага в схемах Рунге – Кутты [17]:

(4)

добавка C в знаменатель реализует стандартный для решателей ОДУ [18] смешанный допуск (аналог RelTol + AbsTol), предотвращающий нефизический рост погрешности вблизи нулей функции. Теоретическое обоснование подобных оценок для разностных схем приведено в [19, с. 367–379, 392–400].

Чтобы обеспечить воспроизводимость расчетов и математическую строгость оценок, параметры регуляризации определяются априори на основе свойств эталонного аналитического ряда (2). Величина параметра C вычисляется по формуле

Характерный масштаб решения Fmax определяется как супремум модуля нормированной концентрации по области, выбранной для сравнения Fmat и Fser:

при фиксированном u. В начальный момент времени t = 0 согласно постановке задачи концентрация однородна: F(z,0) = 1 для всех z ∈ [0;1]. Поскольку физический процесс описывает направленное осаждение частиц под действием дрейфа, функция F(z,t) монотонно убывает во времени в любой точке внутри рассматриваемой области. Следовательно, глобальный максимум достигается при t = 0, что дает строгое равенство Fmax = 1,0. Следовательно,

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

Сформулируем правило выбора ε. Для заданных параметров u и сетки ωz × ωt проводится тестовая серия вычислений ряда (2). Выбирается базовое число слагаемых N. Для каждого узла сетки рассчитывается разность между частичной суммой при N и 2N слагаемых. Максимальное относительное расхождение по всей области принимается за порог неопределенности эталона:

где F(N) – частичная сумма ряда (2) из N слагаемых.

Данная процедура выполняется однократно до начала сравнения с численным методом MATLAB. Для N = 105 получено ε = 0,051 при u = 8770 и ε = 0,31 при u = 70300.

Добавка C предотвращает деление на близкие к нулю значения функции F(z,t) на поздних этапах процесса, когда концентрация становится экспоненциально малой. Величина C интерпретируется как «уровень шума» или фундаментальная точность эталона. Если фактическая разница между численным и аналитическим решением оказывается меньше ε, это означает, что численный метод достиг предела различимости относительно выбранного аналитического решения, а не просто приблизился к нему. Выбор конкретного N = 105 обусловлен балансом между машинным временем расчета суммы ряда и насыщением кривой ε(N) – дальнейшее увеличение числа слагаемых сверх этого порога не приводит к существенному снижению ε в формате двойной точности.

Агрегация погрешности по всему расчетному интервалу выполняется путем поиска максимума погрешности

,

что отражает наихудший сценарий расхождения алгоритмов.

Результаты исследования и их обсуждение

Для верификации предложенного численного метода выполнены расчеты функции F(z,t) при значении параметра дрейфа u = 17,2. Укажем физические параметры безразмерности. Значение параметра дрейфа u = 17,2 найдено из параметров, приведенных в [20]: l = 0,03 м, D = 1,27·10–12 м2/c, v = 7,28·10–10 м/c. Конечное время tend = 0,08 соответствует физическому времени tф ≈ 5,67·107 с.

На рис. 2 приведено сравнение результатов, полученных двумя независимыми способами: суммированием аналитического ряда (2) в Maple и прямым численным решением задачи (1) в MATLAB. Сравнение выполнено для трех характерных координат: z = 0, z = 1/2 и z = 8/9.

При u = 17,2 на интервале безразмерного времени 0 ≤ t ≤ 0,08 максимальная относительная погрешность не превысила 6 %, что подтверждает адекватность численной схемы.

Рис. 2. Сравнение численного и аналитического решений при u = 17,2 Примечание: составлен авторами по результатам данного исследования

Рис. 3. Трехмерный график численного решения при u = 100 Примечание: составлен авторами по результатам данного исследования

Благодаря усовершенствованному алгоритму, адаптированному к сингулярно возмущенной структуре уравнения (1), удалось получить полные трехмерные графики функции F(z,t) в широком диапазоне параметров. Ранее построение таких графиков по аналитическому ряду (2) было возможно лишь для фиксированных значений координаты z из-за вычислительных ограничений [14]. На рис. 3 представлен пример трехмерного распределения концентрации при u = 100, полученный численным методом.

Параллельно разработан модифицированный алгоритм суммирования ряда (2), позволяющий строить 3D-графики вплоть до u = 500. В ходе расчетов установлено, что при u ≥ 83 в окрестности точки (z = 1, t = 0) функция F(z,t) проявляет вычислительную особенность, связанную с потерей точности при суммировании быстро меняющихся экспоненциальных слагаемых для малых времен. Для устранения этого эффекта алгоритм дополнен процедурой сглаживания, основанной на асимптотическом разложении функции F при t → 0 [13; 15].

Сопоставление численного решения с аналитическим рядом (2) для режимов u = 8770 и u = 70300 выполнено при параметрах регуляризации C = 0,051 и C = 0,31, определенных в разделе «Материалы и методы исследования». При этих фиксированных C вычислена гибридная погрешность δ(t) по формуле (4) на всем интервале t ∈ [0;tend]. Ее максимум

достигается вблизи tend и составляет:

u = 8770: δmax = 0,0503

на интервале t ∈ [5,5·10–5; 6,5·10–5];

u = 70300: δmax = 0,3022

на интервале t ∈ [6,0·10–6; 7,5·10–6].

В обоих режимах δmax < C, то есть фактическая погрешность не превышает порога регуляризации. Близость δmax к C означает, что расхождение численного и аналитического решений находится на уровне собственной погрешности эталонного ряда (2) и не может быть измерена точнее в рамках выбранного N. Таким образом, численный метод согласован с эталоном в пределах его точности. При u = 70300 и z = 1/2, начиная с t = t0 = 7,213·10–6, функция принимает малые значения |F| < C, регуляризатор C ограничивает знаменатель формулы (4) и предотвращает нефизический рост δ(t).

Приведем параметры безразмерности. Значение параметра дрейфа u = 8770 найдено из параметров, указанных в [20]: l = 0,03 м, D = 4,5·10–13 м2/c, v = 1,31·10–7 м/c. Конечное время tend = 6,5·10–5 соответствует физическому времени tф ≈ 1,3·105 с.

Параметр u = 70300 выбран авторами для проверки численного метода в экстремальном режиме.

Рис. 4. Сравнение численного и аналитического решений уравнения (1) для режима сильного дрейфа (u = 8770) в точке наблюдения z = 1/2. Параметры расчета: пространственная сетка Nz = 1800, временная сетка Nt = 18000. Допуски встроенного решателя MATLAB pdepe: RelTol=10-3, AbsTol=10-6. Безразмерные переменные (z,t,F). Порог гибридной погрешности C = 0,051 обозначен горизонтальной линией. Панели: (а) графики функции F(z,t); (б) разность решений Fmat – Fser; (в) гибридная погрешность δ(t) Примечание: составлен авторами по результатам данного исследования

Рис. 5. Сравнение численного и аналитического решений уравнения (1) для экстремального значения (u = 70300) в точке наблюдения z = 1/2. Параметры расчета: пространственная сетка Nz = 1800, временная сетка Nt = 18000. Допуски встроенного решателя MATLAB pdepe: RelTol=10-3, AbsTol=10-6. Безразмерные переменные (z,t,F). Порог гибридной погрешности C = 0,31 обозначен горизонтальной линией. Панели: (а) графики функции F(z,t); (б) разность решений Fmat – Fser; (в) гибридная погрешность δ(t) Примечание: составлен авторами по результатам данного исследования

Результаты вычислений приведены на рис. 4 и 5 для z = 1/2. Представлены результаты, полученные для малых значений времени t: рис. 4 для t ∈ [5,5·10–5; 6,5·10–5]; рис. 5 для t ∈ [6,0·10–6; 7,5·10–6].

В численном эксперименте достигнуто хорошее согласование результатов для функции F(z,t), полученных как по аналитическому ряду (2), так и при прямом численном решении начально-краевой задачи (1). Анализируя полученные результаты, можно сделать вывод, что предложенный адаптивный алгоритм позволяет находить достоверные результаты для больших значений параметра дрейфа u без применения методов «длинной арифметики» [21], что сокращает вычислительные затраты.

Для больших значений параметра дрейфа (u = 8770, u = 70300) качественный вид трехмерных графиков функции F(z,t) соответствует асимптотике точного решения и физическому смыслу задачи. Количественный анализ показал, что при z = 0 численное решение демонстрирует нефизические колебания в окрестности t = 0, что связано с ограниченной разрешающей способностью сетки в пограничном слое. При z = 1/2 и z = 8/9 решение ведет себя корректно в соответствии с физической картиной процесса.

При u = 70300 численное решение остается достоверным до момента окончания переходного процесса: t ≈ 7,2·10–6 при z = 1/2 и t ≈ 1,26·10–5 при z = 8/9. Отметим, что для координаты z = 8/9 расчетный интервал времени был умышленно расширен за пределы базового t ≈ 7,5·10–6, так как стабилизация процесса в этой зоне наступает позже. При этом, в силу накопления машинных погрешностей округления при суммировании аналитического ряда в данной области, верификация решения при t > 7,5·10–6 осуществлялась на основе встроенных механизмов контроля устойчивости, адаптивного выбора шага интегрирования и оценки сходимости среды MATLAB (компоненты пакета pdepe), которые подтвердили корректность и вычислительную стабильность получаемых численных результатов вплоть до достижения стационарного режима.

За пределами указанных временных интервалов начинают проявляться численные колебания, не имеющие физического обоснования. Следовательно, при экстремально больших u численное решение целесообразно использовать преимущественно для исследования начальной стадии перехода от начального распределения к стационарному профилю концентрации. Границы применимости данного подхода детально обсуждаются в [14].

Выводы

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

1. Разработан адаптивный численный алгоритм и реализована программа в среде MATLAB для решения одномерной начально-краевой задачи уравнения диффузии с дрейфовым членом и краевыми условиями третьего рода. Алгоритм автоматически выбирает стратегию вычисления в зависимости от величины безразмерного параметра скорости дрейфа u, что обеспечивает устойчивость расчетов в широком диапазоне параметров.

2. Проведена верификация численного метода: при u = 17,2 получено совпадение с аналитическим решением, представленным в виде бесконечного ряда (2), с относительной погрешностью не более 6 %. Это подтверждает корректность реализации численной схемы и адекватность выбранной пространственно-временной дискретизации.

3. Выполнен анализ погрешностей и разработан алгоритм сглаживания вычислительной особенности функции F(z,t) в окрестности точки (z=1, t=0) при суммировании конечного числа слагаемых ряда (2).

4. Проведены вычислительные эксперименты для больших значений параметра дрейфа (u = 8770 и u = 70300), соответствующих условиям реального эксперимента. Показано, что прямое численное решение в MATLAB корректно описывает переходный процесс формирования стационарного профиля концентрации. Однако при z = 0 и малых t наблюдаются нефизические колебания решения, обусловленные ограниченной разрешающей способностью сетки в пограничном слое. Это определяет границы применимости метода: решение достоверно в интервале 0 < t ≲ ttrans(z,u), где ttrans – характерное время перехода.

Использование адаптивной стратегии позволило отказаться от ресурсоемких библиотек произвольной точности. Работа вносит вклад в развитие вычислительных методов моделирования оптически управляемых коллоидных систем и создает методическую базу для дальнейшего исследования процессов переноса в дисперсных средах под действием внешних полей. Полученные результаты могут быть использованы в учебном процессе при изучении численных методов решения сингулярно возмущенных задач диффузии-дрейфа.