Анализ проекционного подхода при решении уравнений Навье-Стокса в задаче моделирования процесса изготовления СТМ-зондов Текст научной статьи по специальности «Физика»

НАНОЭЛЕКТРОННЫЕ ПРИБОРЫ И УСТРОЙСТВА

УДК 621.385.833

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

ШЕЛКОВНИКОВ ЕЮ., ТЮРИКОВ А.В., ГУЛЯЕВ П.В., ЖУЙКОВ Б.Л., ЛИПАНОВ СИ. Институт механики УрО РАН, 426067, г. Ижевск, ул. Т.Барамзиной, 34

АННОТАЦИЯ. Рассмотрены вопросы применения проекционных методов в задаче моделирования процесса изготовления СТМ-зондов. Приведены результаты расчета производительности при различных плотности узлов сетки и числах Рейнольдса для итерационного алгоритма SIMPLE и безытерационного проекционного метода. Показано, что безытерационный метод выигрывает у итерационного по времени расчета при отличии результатов вычислений для обоих методов не более, чем на 5 %.

КЛЮЧЕВЫЕ СЛОВА: сканирующий туннельный микроскоп, зондирующая игла, химическое травление, уравнения Навье-Стокса.

ВВЕДЕНИЕ

Зондирующая игла (ЗИ) является одним из наиболее важных компонентов сканирующего туннельного микроскопа (СТМ), во многом определяющим его технические характеристики, в частности, атомарное пространственное разрешение, необходимое при изучении физико-химических и геометрических свойств поверхности [1]. Исследования геометрических параметров ЗИ проводятся в ИМ УрО РАН в течение продолжительного времени и касаются как наноскопических особенностей острия [2 - 4], так и его макроскопических характеристик [5 - 9]. Изготовление вольфрамовых ЗИ комбинированным методом электрохимического и химического травления [10] включает разработку и анализ математической модели этого процесса, позволяющей корректно описать происходящие в изучаемой системе явления гидродинамического электромассопереноса. Реализация такой модели невозможна без численного моделирования течения жидкости, сопряженного с решением уравнений Навье-Стокса. Их решение в системе с уравнениями диффузии и теплопроводности [8, 9, 11], особенно в трехмерном случае и при наличии сложных геометрических характеристик исследуемой системы, существенно затруднено. Поэтому необходимы анализ и оптимальный выбор методов, применение которых позволяет значительно уменьшить время расчета (не снижая при этом его точности) и упростить анализ модели процесса изготовления СТМ-зондов. Проблемы численного решения уравнений Навье-Стокса в естественных переменных «скорость-давление» хорошо известны [11 - 13] и связаны, в основном, с их нелинейностью (относительно компонент гидродинамической скорости), а также с трудностями в определении поля давления. Для разрешения указанных трудностей применяются различные методы, в основном подразделяемые на проекционные, основанные на принципе расщепления неизвестных скорость-давление по физическим факторам (использующие, например, алгоритм SIMPLE и его модификации [11]), и прямые методы решения. Данная работа посвящена анализу применения проекционных методов в задаче двумерного моделирования течения несжимаемой жидкости в системе, используемой для моделирования процесса химического и электрохимического травления вольфрамовых заготовок СТМ-зондов.

ЧИСЛЕННАЯ МОДЕЛЬ

Моделирование процесса химического травления вольфрамовых ЗИ возможно производить в рамках приближения Буссинеска [13], целесообразность применения которого обоснована в [9]. При этом модель несжимаемой жидкости заменяется моделью слабосжимаемой жидкости, описываемой системой уравнений Навье-Стокса, теплопроводности и диффузии:

— + (д • У)д = vАV--Vр — ¡ЗТ ' д;

V ^ р0

дТ'

дt +(д • V Т ') = хАТ '; (1)

^ + ( )с. = Б У с. + 5

дЬ \ ! 3 33 3

дополненных уравнением неразрывности, записанным в предположении о слабой сжимаемости жидкости:

й-IV (V ) = 0. (2)

При этом зависимость плотности от температуры слабая и полагается:

р (Т) = ро (1 - вТ'). (3)

В системе уравнений (1), (2) введены следующие обозначения: р0 - плотность при

некоторой равновесной температуре Т0 (за Т0 может быть принята начальная температура

травящего раствора); V - гидродинамическая скорость вязкой жидкости; V - кинематическая вязкость; Т' - отклонения температуры от равновесного значения; д - ускорение

свободного падения; х - коэффициент температуропроводности; Б3 - коэффициент диффузии у-го компонента раствора; Б - источниковый член (позволяющий описать,

например, приток «свежего» раствора, его сток, а также «уничтожение» компонентов раствора на границе протекания химической реакции (на границе заготовки ЗИ)); в - коэффициент температурного расширения жидкости.

Учитывая масштабы искомых величин, система (1), (2) может быть приведена к безразмерному виду [13]:

дV ¡д д\д 1 .д ^ д От „

--+ (V • V ) V = —А V — Ур — п-В;

8Ь \ > Ие У Ие2

—+ (V •У 0) = -^- Ав; (4)

дЬ V ' ИеРг

8С. / д\ 1

3 ■ 1 >

+ (дУ)с. =-Ус. + е. ,

V / 3 5С 3 3

дЬ V ' 3 Ие 5с

где 0, Е3 - соответственно, безразмерная температура и источник; Ке,От, Рг, 5с - числа

Рейнольдса, Грасгофа, Прандтля и Шмидта, соответственно; п - направление силы тяжести. Система (4), записанная в декартовых координатах, в двумерном случае может быть представлена в следующем виде (индексу в с опущен для краткости):

du i Arr TT3U au dp n

---AU + U--|-V--+— = 0;

dt Re dx dy dx

dV 1 dV dV dp Gr

—---- AV + U — + V — + dp = — в;

dt Re dx dy dy Re2

дв ггдв тгдв 1

-= U--|-V-=-

dt dx dy Re Pr

Ав;

dc rr dc dc 1

— + U— + V— =-

dt dx dy Re Sc

Ac;

dU

dV л +-= 0,

dx dy где U, V - компоненты скорости.

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

Исследование и анализ эффективности итерационного и безытерационного подходов к задаче решения уравнений Навье-Стокса производились для случая двумерного течения жидкости в сосуде, одной из границ которого являлась поверхность заготовки ЗИ (рис. 1). На рис. 1 также представлены типы границ, для которых должны быть определены граничные условия: I - граница сосуда и заготовки ЗИ, на которой присутствуют условия прилипания; II, III - части границы, на которых организован источник и сток травящей жидкости. На рис.2 приведен пример сетки, используемой для реализации МКЭ рассмотренной системы.

I - границы с условиями прилипания; II, III - границы с «втеканием» и «вытеканием» травящей жидкости, соответственно

Рис. 1. Схематичное изображение исследуемой системы

Рис. 2. Пример сетки, используемой для реализации МКЭ

Начальные и граничные условия для границ, показанных на рис. 1, приведены ниже:

и = V = 0; р = 0; в = с = (6)

UI = VI = 0;

dc dn

= 0;

dU дП

дс

дП

= 0;

II

dV дП

= a,

V>

II

c=c

др дП

= 0;

i,ii,iii

(7)

= y;

Qii = в •

ii

Граничные условия (7) для компонент скорости на границе I описывают условия прилипания, на границе II - присутствие «втекания», на границе III - граничные условия считаются «свободными». Аналогично заданы граничные условия для концентрации примеси (травящего компонента жидкости). Втекание жидкости определенной температуры обусловлено наличием соответствующего граничного условия для в.

Для реализации итерационного метода при получении системы линейных алгебраических уравнений (являющихся дискретными аналогами уравнений системы (5)) применяется метод Галеркина с использованием квадратичных функций формы для компонент скорости, и линейных - для остальных переменных [13]. Ортогонализация невязки к функциям формы приводит к линейным уравнениям, записанным для к-го узла сетки (по парным индексам подразумевается суммирование):

m,

U-и0

m

т

V-V0

т

+ ClkUl + dmkPm = 0

+ ClkVl + fmkPm Smk вn

(8)

(9)

enU + grnVr = o. (10)

Уравнения (8) - (10) являются дискретными аналогами уравнений для компонент скорости и уравнения неразрывности, соответственно. Уравнения для определения температуры и концентрации записываются аналогичным способом. Коэффициенты mlk, clk, dmk, fmk, smk, elm, glm определяются геометрическими характеристиками расчетной сетки, видом функций формы и значений U и V на предыдущей итерации (или предыдущем временном шаге); U0,V0- переменные, определенные на предыдущем временном шаге; т - шаг по времени. Итерационная процедура (алгоритм SIMPE и его модификации[11]) решения трех уравнений на каждом временном шаге (8) - (10) заключается в следующем:

1. Определение полей приближенных компонент U* и V* из (8) с учетом начального приближения р * для давления;

2. Нахождение поправки р' к давлению с учетом удовлетворения скорректированных

* *

скоростей U + U' и V +V' уравнению неразрывности (10), вычисление поля давления

*

р = р + р';

3. Переход к пункту 1 и повторение процедуры определения U *, V * до достижения сходимости;

*

4. Переход следующему временному шагу и учет в качестве р значения давления, определенного на предыдущем шаге по времени;

5. Определение полей температуры в и концентрации c путем решения соответствующих уравнений с учетом найденных для данного временного шага полей U, V, p.

Процедура алгоритма SIMPLE может быть существенно упрощена, если для каждого

временного шага производить следующие действия:

, *

1. Вычисление промежуточной скорости v в предположении, что перенос количества движения в (4) происходит лишь вследствие диффузии и конвекции и определяется уравнением:

dv ~dt

= -(v

(v * • V)

\ v* . v*

v)v + Av .

(11)

При этом схема решения (11) может быть и полностью явной [12], описываемой при шаге по времени т как:

v * v 0

v — v

т

= —( v

(v0 V) v0 + Av0,

(12)

где V0 - поле скорости, определенное на предыдущем временном шаге. При этом поправка к

далее будет определена уже с учетом градиента давления:

v — v

т

= —Vp.

2. Определение поля давления согласно уравнения:

vv*

vv*

= —Ap,

(13)

(14)

т

которое формируется из (13) взятием от обеих частей дивергенции и учетом того, что поле v на данном временном шаге удовлетворяет уравнению неразрывности.

3. Вычисление v из (13) с учетом найденного на предыдущем шаге поля давления p, а также поля концентрации и температуры.

4. Переход к следующему по времени шагу.

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

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

При решении задачи, описываемой системой уравнений (5), с учетом начальных и граничных условий (6), (7) методом МКЭ для итерационной и безытерационной схем были получены следующие результаты. На рис. 3 приведены распределение концентрации травящей жидкости (представлено оттенками серого цвета), определяющееся ее течением в системе, показанной на рис. 1, а также поле скоростей (представлено векторами).

а) б)

Рис. 3. Распределение концентрации травящего компонента и поле гидродинамической скорости (векторы), определяемые решением (5) при установившемся течении: а — путем реализации алгоритма SIMPLE; б — безытерацонным методом

*

*

Значения концентрации вдоль профиля L, приведенного на рис. 3, а, отражены на графике (рис. 4). В таблице представлены результаты расчета производительности (времени вычислений до установления течения) при различной плотности узлов сетки и числах Рейнольдса для итерационного алгоритма SIMPLE и безытерационной проекционной схемы.

0 -

Рис. 4. Нормализованные значения концентрации вдоль профиля L,

полученные при установившемся течении (■ — итерационный алгоритм SIMPLE; ▲ - безытерационная схема)

Таблица 1

Результаты расчета производительности при различной плотности узлов сетки и числах Рейнольдса

№ Узлы Треугольники Re Безытерационная схема, время, с SIMPLE, время, с

1 4128 8001 20 8,1 8,6

2 16448 32385 20 39,6 44,6

3 65664 130305 20 297,9 349,6

4 262400 522753 20 864,8 1745,5

5 1049088 2094081 20 4255,3 9731,2

6 65664 130305 50 305,4 361,3

7 65664 130305 100 321,9 397,1

Из анализа расчетных данных следует, что безытерационная схема выигрывает у итерационной по времени расчета. При этом можно ожидать, что при решении трехмерной задачи этот выигрыш станет гораздо более существенным. Точность же хорошо зарекомендовавшего себя алгоритма SIMPLE (реализованного в рамках МКЭ) согласно рис. 4 сравнима с результатами, формируемыми безытерационной схемой, отличаясь от них не более, чем на 5%. Несмотря на то, что введение уточняющих итераций на каждом шаге по времени позволяет осуществлять моделирование в условиях поставленной задачи с большей точностью, использование в рамках проведения численного эксперимента безытерационного метода при решении задачи травления (особенно, трехмерной) вольфрамовой заготовки ЗИ в рамках приближения Буссинеска представляется более перспективным.

СПИСОК ЛИТЕРАТУРЫ

1. Binnig G., Rohrer H. Scanning tunneling microscopy // Helvetica Physica Acta. 1982. V. 55. P. 726-735.

2. Липанов А.М., Тюриков А.В., Шелковников Е.Ю., Гудцов Е.В. Исследование разрыва «шейки» заготовки зондирующей иглы СТМ при ее изготовлении методом химического травления // Химическая физика и мезоскопия. 2005. Т. 7, № 2. С. 162-168.

3. Липанов А.М, Тюриков А.В., Шелковников Е.Ю., Гудцов Д.В., Гуляев П.В. Численные исследования микротопологии острия зондирующей иглы СТМ при его формировании электрохимическим методом // Ползуновский Альманах. 2006. № 4. С. 45-46.

4. Липанов А.М., Тюриков А.В., Шелковников Е.Ю., Гуляев П.В., Осипов Н.И. Моделирование влияния остроты зондирующего эмиттера на диффузионную составляющую тока Фарадея при электрохимических СТМ-исследованиях // Химическая физика и мезоскопия. 2012. Т. 14, № 2. С. 292-295.

5. Липанов А.М., Тюриков А.В., Шелковников Е.Ю., Кизнерцев С.Р., Мышкин О.И. Методика моделирования процесса электрохимического травления зондирующих игл сканирующего туннельного микроскопа // Химическая физика и мезоскопия. 2010. Т. 12, № 2. С. 281-285.

6. Липанов А.М., Тюриков А.В., Гуляев П.В., Шелковников Е.Ю., Осипов Н.И. Исследования влияния геометрических характеристик зондирующего острия электрохимического СТМ на величину регистрируемого им фарадеевского тока // Химическая физика и мезоскопия. 2011. Т. 13, № 2. С. 215-219.

7. Липанов А.М., Тюриков А.В., Суворов А.С., Шелковников Е.Ю., Гуляев П.В., Кизнерцев С.Р., Жуйков Б.Л. Метод исследования химического травления заготовок измерительных игл туннельного микроскопа // Химическая физика и мезоскопия. 2007. Т. 9, № 2. С. 172-182.

8. Шелковников Е.Ю., Тюриков А.В., Гуляев П.В, Жуйков Б.Л., Липанов С.И. Исследование трехмерной диффузионно-конвективной модели химического травления зондирующих острий СТМ // Химическая физика и мезоскопия. 2013. Т. 15, № 2. С. 204-209.

9. Шелковников Е.Ю., Тюриков А.В., Гуляев П.В., Жуйков Б.Л., Липанов С.И. Схема численного исследования влияния тепловой гравитационной конвекции на процесс травления зондов СТМ // Химическая физика и мезоскопия. 2013. Т. 15. № 4. С. 645-649.

10. Шелковников Е.Ю., Тюриков А.В., Кизнерцев С.Р., Липанов С.И. Методика моделирования процесса комбинированного травления нанозондов // Химическая физика и мезоскопия. 2012. Т. 14, № 3. С. 143-146.

11. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М. : Энергоиздат, 1984. 151 с.

12. Белоцерковский О.М. Численное моделирование в механике сплошных сред. М. : Физматлит, 1994. 448 с.

13. Полежаев В.И., Бунэ А.В., Верезуби Н.А. и др. Математическое моделирование конвективного тепломассообмена на основе уравнений Навье-Стокса. М. : Наука, 1987. 272 с.

14. Сегерлинд Л. Применение метода конечных элементов. М. : Мир, 1979. 392 с.

ANALYSIS OF A PROJECTIVE APPROACH FOR SOLVING THE NAVIER-STOKES EQUATION FOR PROCESS OF MODELING THE MANUFACTURING STM PROBES

Shelkovnikov E.Yu., Tyurikov A.V., Gulyaev P.V., Zhuikov B.L., Lipanov S.I.

Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russian

SUMMARY. The aspects of projection methods in the problem of modeling the manufacturing process STM probes are observed. Results of calculation of performance under varying density of grid points and Reynolds numbers for the iterative algorithm SIMPLE and noniterative projection method are presented. It is shown that noniterative method outperforms the iterative calculation of time with difference in the results of calculations within 5%.

KEYWORDS: scanning tunnel microscope, the probe, chemical etching, Navier-Stokes equations.

Шелковников Евгений Юрьевич, доктор технических наук, профессор, зав. лабораторией ИМ УрО РАН, профессор кафедры «Вычислительная техника» ИжГТУ имени М.Т. Калашникова, e-mail: iit@udman.ru

Тюриков Александр Валерьевич, кандидат физико-математических наук, старший научный сотрудник ИМ УрО РАН, e-mail: alex.tyurikov@mail.ru

Гуляев Павел Валентинович, кандидат технических наук, старший научный сотрудник ИМ УрО РАН, e-mail:lucac@e-izhevsk. ru

Жуйков Богдан Леонидович, аспирант ИМ УрО РАН, e-mail: fastblood@mail.ru Липанов Святослав Иванович, аспирант ИМ УрО РАН, e-mail: HPotterS@yandex.ru