Управление температурным полем и его прогнозирование в нанокомпозиционных материалах Текст научной статьи по специальности «Физика»

УДК 519.63

Управление температурным полем и его прогнозирование в нанокомпозиционных материалах

© Н.И. Сидняев, Ю.С. Ильина, Д. А. Крылов МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

Ключевые слова: нанокомпозит, температурное поле, уравнение, нанотрубки, управление, математическое моделирование, численные методы.

Введение. Исследование свойств вещества, находящегося в нано-размерном состоянии, вызывает повышенный интерес и имеет большое значение как для развития фундаментальной науки, так и для практического применения таких наноматериалов в устройствах нано- и микросистемной техники [1-3]. Уменьшение характерных размеров частиц до величин, соответствующих значению длины волны де Бройля в твердом теле, приводит к квантованию энергетических уровней и сильному изменению поляризуемости частиц. Происходящее при этом возрастание роли релаксации поверхностных атомов сопровождается изменением электронной структуры точечных дефектов. При уменьшении диаметра частиц отношение площадей их поверхностей к внутреннему объему повышается, растет доля поверхностных атомов. Это приводит к смене условий фазовых равновесий, уменьшению значений температуры плавления, изменению пределов растворимости, сдвигу

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

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

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

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

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

На рис. 2 представлены исследования нанокомпозитов методом атомно-силовой микроскопии [3]. Так, например, при протекании прак-

тически необратимых реакций поликонденсации происходили рост фрактальных агрегатов и уменьшение коли -чества возможных перестановок между частицами при за- Рис. 1. Структура нанотрубки креплении их на остове фракталов. Это приводило к уменьшению энтропии смешения, что, в свою очередь, обусловливает увеличение значения энергии смешения (изменение энергии Гиббса) [5, 6].

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

В разработанной конструкции охладителя вода играет роль хладагента, подобно фреону в холодильных установках. Из-за малого (около 50 нм) диаметра капилляров, в роли которых выступают углеродные нанотрубки, вода, проходящая сквозь них, полностью испаряется, сни-

Рис. 2. Общий вид фазовых диаграмм для полимерной системы 1 — граница метастабильной области; 2—граница области неустойчивости; 3 — бинодальная область; 4 — синодальная область

мая таким образом с чипа большее количество тепла (рис. 3, 4). Избыточное давление получившегося пара, который затем конденсируется в холодильнике, и является той движущей силой, которая заставляет циркулировать воду по всей системе охлаждения [7].

Эта система охлаждения является уже не первой системой, созданной учеными университета Пурду для охлаждения горячей электроники. Некоторое время назад уже была разработана система воздушного охлаждения, поверхность радиаторов которого представляла собой «лес» из углеродных нанотрубок. Благодаря этому эффективная площадь радиаторов увеличивалась в несколько сотен раз, позволяя системе сбрасывать в воздух большое количество теплоты и повышая тем самым ее эффективность.

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

Конденсация

Адиабатический перенос тепла

Безпузырьковое кипение жидкости

Рис. 3. Движение воды в углеродной нанотрубке, приводящее к отводу теплоты с поверхности нагреваемого элемента

Рис. 4. Схема системы охлаждения на основе углеродных нанотрубок

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

На границах расчетного объема могут быть заданы граничные условия первого, второго или третьего типов. Кроме того, будем рассматривать наноматериалы с внедренными в них нанотрубками, которые будут позиционироваться как дополнительные источники (стоки) теплоты.

Базовые уравнения и постановка задачи. Основные уравнения.

Базовым уравнением для разработки методов прогноза процессов распределения температурных полей является уравнение баланса теплоты в интегральной форме [8]:

д

—Шу = - Ц цп а а+ш еу. (1)

дт V Е V

Здесь V — контрольный объем; Е — окружающая его поверхность; Т — время; Н = Н(и) — энтальпия (внутренняя энергия) на единицу объема; и — температура; ц — вектор потока тепла; п — внешняя нормаль к рассматриваемому объаму; Е — приток теплоты на единицу объема. В грунте поток теплоты ц определяется по закону Фурье:

Ц = -Х §гаё и, (2)

где X — коэффициент теплопроводности. Фазовые переходы и скачки коэффициентов теплопроводности на границах наноматериалов различных типов или элементов строительных конструкций порождают разрывы в решениях. При отсутствии таких разрывов на основе уравнений (1) и (2) выводится уравнение распространения теплоты в дифференциальной форме:

д д д д д д д Су~и = — — и) + — & — и ) + —& — и) + Е (Х У, г), (3) дТ дх дх дУ ду дг дг

где су = р с — объемный коэффициент теплопроводности; с — удельная теплоемкость; р — плотность; X — коэффициент теплопроводности, и(х, у, г, Т) — температура среды. Здесь считается, что в случае наличия фазового перехода, где энтальпия терпит скачок, объемный коэффициент теплопроводности имеет вид

су = р с + QЬ(u -и*),

где и* — температура фазового перехода; Q — теплота фазового перехода; 5(и - и*) — дельта-функция.

Энтальпия в общем случае — нелинейная, монотонно возрастающая функция от температуры. В простейшем случае можно считать, что

Н = Н0 + с (и - и*), и > и*,

Н = Н0 + Су (и - и*) - Q, и < и*.

Здесь сл — объемная теплоемкость жидкой фазы; су — объемная теплоемкость твердой фазы; Н0 — константа определения энтальпии. Если Н0 = 0, энтальпия равна нулю при и = и*.

В случае, когда коэффициент теплопроводности — постоянная величина, из уравнения (3) выводится уравнение теплопроводности:

Су д д2 д2 д2 Р(х,у,г) (4)

--и =—-и + —-и + —-и + у у . (4)

X дл дх2 ду2 дх2 X

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

Расчеты разрывов. Применяют два разных аналитических подхода и соответственно два разных типа численных методов для решения задач с разрывами. При первом подходе разрывы выделяются, а непрерывные области описываются уравнениями (1) или (2). На разрывах ставятся некоторые граничные условия, выводимые из (1) и (2). Для фазовых переходов граничные условия следующие [8, 9]:

и(Н2 - Н1) = XIГи | - X2 и | , и = 0, и2 = 0. (5)

Здесь и — скорость распространения поверхности разрыва в направлении ее нормали индексы 1 и 2 обозначают значения соответствующих величин по разные стороны разрыва; нормаль ? направлена от области 1 к области 2. Предполагается, что на таком разрыве терпит скачок величина энтальпии при и = 0 в связи с наличием теплоты фазового перехода. Учитывается также скачок теплопроводности для твердого и жидкого вещества, но для существования разрыва наличие такого скачка не обязательно.

Для неподвижных разрывов, связанных со скачком коэффициента теплопроводности на границе раздела двух разных сред, граничные условия следующие [9, 10]:

Хх1— и | -X21— и | = 0, щ = и2. (6)

^д? ) \дэ )2

В одномерном случае, когда все зависит только от одной координаты х, эти соотношения принимают вид

и(Н - И2) = Хх и) -X2^и) , щ = 0, щ = 0.

х (1и )1-х 2 [1и и1=и2'

Граничные условия на границе расчетного объема и их применение для расчетов сооружений. На границе расчетного объема ставят три типа граничных условий.

Граничное условие первого типа:

и = ц(0. (7)

Пример реализации — поверхность с заданной температурой. Граничное условие второго типа:

д

-X —и = у(0. (8)

дп

Примеры применения: плоское терморегулирующее устройство, позволяющее регулировать поток теплоты с поверхности, при у(0 = 0 — теплоизолирующая поверхность, стандартное граничное условие на боковой границе расчетного объема.

Граничное условие третьего типа:

д

X~dnU = ~°(Ы У ■ (9)

Применяется для описания области контакта с двусторонней тепло-проводящей поверхностью, для которой, с одной стороны, температура совпадает с температурой наноматериала, а с другой — поддерживается заданная температура |i(t). Данное условие моделирует тонкий слой теплопроводящего наноматериала. Величина о рассчитывается по формуле о = (km — теплопроводность наноматериала, h — его толщи-h

на). Пример применения — наноконструкция, подверженная температурным перепадам с внешней и внутренней стороны.

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

Пусть температура наноматериала определяется функцией g(t). Другими словами, мы решаем задачу

(УЧП) ut = uxx, 0 < x < 1, 0 < t <ю, \u (0, t) = 1,

(ГУ) \ r , 0 < t <ю. (10)

( У) \ux(1,t) = -[u(1,t)-g(t)], (10)

(НУ) u(x, 0) = 0, 0 < x < 1.

Для ее решения методом конечных разностей [2] построим прямоугольную сетку, узлы которой определяются формулами (см. рис. 5)

Xj = jh, j = 0, 1, 2, ..., n.

у, = ik, i = 0, 1, 2, ..., m.

Отметим, что значения Uj на левой и нижней сторонах сетки на рис. 5 известны из граничных и начальных условий и наша задача состоит в отыскании остальных значений щ. Для этого заменим частные производные в уравнении теплопроводности их конечно-разностными аппроксимациями:

щ =1 [и(х, г + к) - и(х, г)] = 1 [и,} - и,- } ], к к

ихх = -1 [и(х + Л, г) - 2и(х, г) + и(х - И, г)] = -1 [и,,]+ - 2и,,] + и,,}-1 ].

И И

Подставим эти выражения в уравнение иг = ихх и решим получившееся уравнение относительно значений функции на верхнем временном слое. В результате получаем

к

и,+1, ] = и,, ] + — [и,, ]+1 - 2и,, ] + и,, ]-Л. (11)

И

Это и есть искомая формула, поскольку она выражает решение в данный момент времени через решение в предыдущий момент времени (индекс г относится к временной переменной). На рис. 5 выделены те значения, которые входят в данную формулу.

Теперь можно приступить к вычислениям. Однако сначала необходимо аппроксимировать производную в граничном условии на правом конце:

их (1, г) = -[и (1, г) (г) ].

В результате аппроксимации получаем

И [и, ,„ - и, ,„-1 ] = - [и, ,„ - g1], (12)

где значения gi = g(ik) известны. Здесь мы заменили их(1, г) левой разностной производной, поскольку правая разностная производная потребовала бы значений функции за пределами сетки.

Из (12) находим

и, ,„-1 - hgi

. 03)

Формулы (11) и (13) позволяют начать вычисления.

Алгоритм вычислений по явной схеме. Шаг 1. Находим решение на сеточном слое г = Дг, используя явную формулу (рис. 6):

к

и2,] = их,] + — [их,]+1 - 2их,] + их,]-1 ], ] = 2,3,..., п -1.

Шаг 2. Величину и2,п находим по формуле (13)

и.

Чй

и.

2,и

п-1 + 2

1 + И

Рис. 5. Сетка для уравнения теплопроводности

1 = Т

11=1/ (п-1) к=Т/(т-1)

После совершения шагов 1 и 2 получаем решение для Л = ДЛ. Для получения решения при Л = 2ДЛ (вторая строка снизу на рис. 6) повторяем шаги 1 и 2, поднявшись на одну строку вверх, т. е. увеличив г на 1 и используя и^ с предыдущей строки. Аналогично вычисляется решение и в последующие моменты времени Л = 3ДЛ, 4ДЛ, ...

Замечания. 1. У явной схемы есть серьезный недостаток. Если шаг по времени оказывается достаточно большим по сравнению с шагом по х, погрешности округления могут стать настолько большими, что полученное решение потеряет смысл. Отношение шагов по Л и х зависит от уравнения и граничных условий, но в общем случае шаг по времени должен быть много меньше шага по координате. Доказано, что для применимости явной схемы должно выполняться условие к/к2< 0,5.

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

3. Для гиперболической задачи

Рис. 6. Иллюстрация явной разностной схемы:

п — число узлов сетки вдоль оси х; т — число узлов сетки вдоль оси у; И = 1/(и - 1), к = Т/(т - 1)

Уравнение в частных производных (УЧП) иа = ихх, 0 < х < 1, 0 < Л <<х>,

\и (0, Л) = gl(t) Граничные условия (ГУ) 1 ^ ^) = ^) 0 <Л < ^

Начальные условия (НУ) \и (х,0) = х),

к (х, 0) = у(х), 0 < х <1,

также можно построить явную разностную схему. Для этого аппроксимируем иа и Ыхх центральными разностными производными:

иа = -1 [и(х, г + к) - 2и(х, г) + и(х, г - к)], к

ихх = [и(х + к, г) - 2и(х, г) + и(х - к, г)], а начальное условие — по схеме

1 г , ^ . 1

и,

(х, 0) = — [и(х, к) - и(х, 0)] = — [и(х, к) - ф(х)].

В результате для вычисления величины и(х, г + к) получаем следующую явную схему:

и (х,г + к) =

= 2и(х, г) - и(х, г - к) + ^к^ [и(х + к, г) - 2и(х, г) + и(х - к, г)]. (14)

Из (14) видно, что для вычисления на данном временном слое необходимо знать решение на двух предыдущих слоях. Значит, для начала отсчета необходимо воспользоваться начальным условием на скорость:

1 [и( х, к) -ф( х)] = у( х),

из которого получаем и(х, к) = ф(х) + ку(х), т. е. значение при г = Дг. Решение в последующие моменты времени можно находить по явной формуле (14).

Неявная схема для уравнения теплопроводности. Рассмотрим задачу [2, 7]:

(УЧП) щ = их, 0 < х < 1, 0 < г < да,

Ги(0, г) = 0, 0 (ГУ) { 0 < г < да,

К (1, г) = 0,

(НУ) и(х, 0) = 1, 0 < х < 1.

(15)

Воспользуемся следующими конечно-разностными аппроксимациями для частных производных иг и ихх:

щ (х, г) = 1 [и(х, г+к) - и(х, г)], X

ихх(х, г) =—|и(х + И, г + к) - 2и(х, г + к) + и(х - И, г+к)] + И2

(1 — А )

+ [и(х + И, г) — 2и(х, г)+и(х — И, г)],

где X выбирается из отрезка [0, 1]. Отметим, что ихх аппроксимируется взвешенным средним центральных разностных производных в момент г и г + к. При X = 0,5 получается обычное среднее этих двух центральных производных, а при X = 0,75 одна из разностных производных берется с весом 0,75, вторая — с весом 0,25. При X = 0 получается обычная явная схема.

После замены частных производных щ и ихх в задаче (15) получаем разностную задачу.

Разностное уравнение

1 (и1+1,} — и,,}) = к

X, , (1 — X),

= — (иг-+1, ]+1 — 2иг+1,} + иг+1, ;—1) + 7 2 и ]+1 — 2иг-} + и,, ;—1), (16) И И

[и, 1 = 0,

(ГУ) \ ' ' г = 1,2,..., т,

[и,, п = 0,

(НУ) щ,, = 1, ] = 2,..., п — 1.

Перенесем все неизвестные значения и верхнего временного слоя (с индексом г + 1) в левую часть уравнения (16) и получим

—Хги,+1, ]+1 + (1 + 2 гХ) и, +1,} — Хги, +1, ]—1 = = г (1 — X) и,г ]+1 + [1 — 2 г (1 — X) ] ] + г (1 — X) и,, ]■_ 1, (17)

где введено обозначение г = к/И2. Отметим, что если г фиксирован, а] изменяется от 2 до п - 1, соотношения (17) определяют систему п - 2 уравнений с п - 2 неизвестными иг+12, иг+13, иг+14, ..., иг+1; п-1, которые являются решением задачи во внутренних узлах сетки на временном слое г = (г + 1)Дг. Наглядное представление о структуре каждого уравнения системы (17) дает рис. 7.

Перейдем к решению системы (17).

г

г + 1 г

Шаг 1. Выбираем некоторое значение Х(0 < X < 1). Если X = 0, то уравнения (17) переходят в явные формулы из предыдущего раздела.

Шаг 2. Пусть, например, И = Дх = 0,2 и к = Дt = 0,08 (при этом г = = к/И2 = 2). В данном случае сетка содержит шесть узлов вдоль оси х (четыре внутренних узла). Возьмем весовой параметр X = 0,5 (получающаяся при этом схема называется схемой Кранка — Николсон). В соответствии с вычислительным шаблоном на рис. 7, двигаясь слева направо (/ = 2, 3, 4, 5) по первым двум слоям (/ = 1), получаем следующие четыре уравнения:

—М21 + 3Ы22 — Щз = Мц — и\2 + Щз = 1, —М22 + 3М23 — М24 = М12 — М13 + М14 = 1, —М23 + 3М24 — М25 = М14 — М15 + М16 = 1, —М24 + 3М25 — М26 = М14 — М15 + М16 = 1.

Перепишем их в матричной форме:

" 3 —10 0 —13 —10 0 —13 —1 0 0 —13

Матрица этой системы называется трехдиагональной. Для того чтобы решить трехдиагональную систему

¿1 С1 0 0 ••• 0 Х1

¿2 С2 0 ••• 0 Х2 й2

0 а2 ¿3 С3 ... 0 Х3 = й3

0 0 0 ап-\ Ъп хп йп

-г% + 1+2 гк + -г% г(1-Х) + 1 -2г(1 -А,) + г(1 -А.)

_I_I_I_

У-1 У 7+1 Л

Рис. 7. Шаблон для неявной схемы

М22 1

М23 1

М24 1

М25 1

(18)

преобразуем ее в эквивалентную систему вида

1 * С1 0 0 ... 0 Х1

0 1 * С2 0 ... 0 Х2 "2

0 0 1 * сз ... 0 Хз = "з*

0 0 0 0 1 хп И * йп

» С1 » С ] + 1 1 о о

где С1 = —, с = ---, ] = 1,2,..., п - 2,

¿1 Ь}+1 - а]с]

" - = й " * = +1 - 1 = 19 п 1

"1 = —, а]+1 =~-г, •/ = п -1

¿1 ¿;+1 - а;су

Новая система полностью эквивалентна исходной, а ее матрица устроена так, что система решается очень легко. Работая с уравнениями последовательно сверху вниз, получаем

хп = й-, х1 = й) -с-Ху+ь} = п -1, п -2, ..., 2,1.

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

Конечно-разностные численные методы сквозного счета температурных полей в двухфазных средах. Для расчета непрерывных областей применяется любая сходящаяся численная схема, аппроксимирующая уравнения (3) или (4), но схемы, выведенные на основе уравнения (1), являются консервативными и потому более точными. Реализация условий (5) требует либо введения сетки, привязанной к разрыву, так называемой плавающей сетки, либо специальных формул расчета для ячеек вблизи поверхности разрыва. В связи со сложностями программной реализации в двумерном и трехмерном случае такой подход для расчета фазовых переходов обычно не используется.

При втором аналитическом подходе разрывы специально не выделяются и рассматриваются обобщенные решения (1). Численные методы, основанные на этом принципе, называются методами сквозного счета [2, 7]. При построении разностных схем в этих методах всегда применяется интегральное уравнение (1). В свою очередь, численные методы сквозного счета подразделяются на методы с регуляризацией и без нее. Регуляризация означает, что расчет разрывных решений заме-

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

Общие требования к конечно-разностным численным схемам. После разложения величин в расчетных ячейках в ряд Тейлора относительно некоторого расчетного узла сетки и подстановки их в исходное уравнение должна выполняться аппроксимация этого уравнения. Порядок остатка от такой подстановки — порядок аппроксимации схемы. Ошибки расчета не должны нарастать с течением времени, что проверяется с использованием спектрального признака или экспериментально. При измельчении шагов сетки численное решение в пределе, как правило, дает точное решение исходных уравнений. Для линейных уравнений доказано [2, 7], что при выполнении устойчивости и аппроксимации сходимость достигается. В нелинейных задачах, к которым относятся задачи с фазовыми переходами, сходимость обычно проверяется экспериментально путем измельчения сетки и сравнения результатов с известными решениями.

При исследовании задач, содержащих разрывы [7, 11] (в данном случае это фронты перехода фаз и скачки коэффициента теплопроводности), необходимо пользоваться консервативными схемами, т. е. при выводе схем рекомендуется использовать интегральные законы сохранения (в данном случае это уравнение (1)). При расчетах задач с разрывами и большими градиентами значимо свойство монотонности численной схемы. Оно означает, что схема не создает максимумов и минимумов, не имеющихся в точном решении. В немонотонной схеме при наличии разрывов и больших градиентов возникает эффект «пилы», образования множества максимумов и минимумов, приводящего к существенному снижению точности расчета или даже к потере устойчивости из-за переполнения.

В численной схеме, используемой в программном комплексе, учтены факторы, описанные в данном разделе.

Анализ алгоритмов и методов, численная схема. Рассмотрим консервативную численную схему [12, 13]. Расчетная ячейка имеет вид параллелепипеда. Узлы сетки, соответствующие его вершинам, обозначаются целыми числами. Значения температуры, коэффициентов

теплопроводности, теплоемкости относятся к серединам расчетных ячеек, эти узлы обозначаются полуцелыми индексами. В них задаются координаты х,, у,, г,, шаги сетки Ихг+1/2 = х,+1 - х,, Иу+1/2 = у,+1 - у,, Игк+1/2 = гк+1 - гк. После применения уравнения (1) к расчетной ячейке получается соотношение:

Их1+1/2 Иу]+1/2 Игк +1/2 (Нг?+1/2, ]+1/2,к +1/2 — 1/2, ]+1/2,к +1/2) =

= {[Иу]+1/2 Игк+1/2 (^,+1, ] +1/2,к+1/2 — ] +1/2,к+1/2) + Иу+1/2 Ик+1/2 (Ц,+1/2, ] +1,к+1/2 — (19)

+1/2, ] ,к+1/2 ) + Их,+1/2 Иу]+ 1/2 +1/2, ] +1/2,к+1 — Ц+1/2, ]+1/2,к )] + +ИхI +1/2Иу]+1/2Игк+1/2^+1/2, ]+1/2,к + 1/2}^

После деления на И^НуИ^ получаем шаблон для численных схем:

,]+1/2,к+1/2 Н,+1/2,]+1/2,к +1/2) / ^ =

= [(Ц+1 ,]+1/2,к +1/2 — ]+1/2,к+1/2)/ Их+1/2, +

+(<?,+1/2 ,]+1,к+1/2 — $,+1/2,],к+1/2)/ Иу+1/2 +

(20)

+(Ц +1/2,] +1/2,к +1 — Ц + 1/2,] +1/2,к )/ И2к+1/2] +

++1/2, ] +1/2,к+1/2.

Здесь приведена разностная схема для неоднородной сетки, реализованная в программном комплексе. Расчет ведется в два этапа: сначала вычисляются потоки, а затем рассчитывается энтальпия, что удобно при постановке второго и третьего граничных условий. Различные консервативные численные методы отличаются один от другого различными способами вычисления потока ц. Потоки вычисляются путем аппроксимации закона Фурье:

Ц,,]+1/2,к+1/2 = X/,]+1/2,к +1/2 (и,+ 1/2,]+1/2,к+1/2 — и,—1/2,]+1/2,к+1/2 V ^ = (21) = (И,+1/2 + И,—1/2)/2.

Здесь в качестве примера приводится аппроксимация потока только в направлении оси х, остальные потоки аппроксимируются аналогично.

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

и одинаковыми шагами по всем пространственным переменным для явной схемы — это условие т < 0,5(ИХ + Иу + И2) /3.

Неявные схемы более устойчивы, а чисто неявная двухслойная схема абсолютно устойчива. Неявные схемы требуют специальных алгоритмов для решения неявных уравнений [2, 10]. При их программной реализации решаются системы неявных уравнений с помощью итерационных алгоритмов, которые нередко сводятся к решению уравнения переноса тепла с использованием явных схем. При применении равномерных сеток это обычно не оправдано из-за большого расхода машинного времени, сопоставимого с расчетами с применением явных схем с мелким шагом. Расчеты по явным схемам за счет мелкого шага оказываются более точными. Введение неявности затрудняет модификацию программ и включение в них различных дополнительных блоков. Но при наличии в расчете областей с существенно различающимися значениями теплопроводности и теплоемкости или в случае неравномерных сеток при наличии областей с существенно отличными пространственными шагами применение неявных схем становится целесообразным.

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

Аналогичные схемы [7, 14] существуют и для уравнения (3). Однако при расчетах фазовых переходов с помощью уравнения (3) регуляризация необходима всегда, что заставляет для достижения точности делать мелкий временной шаг. В силу квазилинейного характера уравнения (3) при использовании неявных схем здесь возникают системы линейных уравнений. Это позволяет построить неявную схему, называемую методом переменных направлений с последующей прогонкой. Расход машинного времени при выполнении одного временного шага здесь имеет тот же порядок, что и у явных схем. Однако включение в такую схему дополнительных рассчитываемых объектов, на которых надо выставлять граничные условия, сильно усложняет алгоритм. По этой причине в программном комплексе реализованы явные схемы.

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

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

В программном комплексе реализован способ, являющийся оптимальным способом вычисления коэффициентов теплопроводности [7]:

Я(,]+1/2,к+1/2 -

- 2(мг+1/2,]+1/2,к+1/2 _ и-1/2,]+1/2,к+1/2)/

^ Ихг+1/2, ]+1/2,к+1/2 + ИХ(-1/2, ]+1/2,к + 1/2^

^ X/ +1/2,]+1/2,к+1/2 Хг_1/2,]+1/2,к+1/2 )

.(22)

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

Разностные методы, в которых пространственные шаги зависят от положения расчетного узла в пространстве, называются неравномерными. В простейшем случае, который рассматривался выше,

Их - Их(х), Иу - Иу (у), к - к (г). (23)

Более сложный вариант

Их - Их (х, у, 2), Иу - Иу (х, у, 2), Иг - И г (X, у, 2). (24)

Неравномерные сетки позволяют осуществить привязку сетки к расчетным объектам и расчет в отдельных областях с мелким шагом. Однако это не является универсальным методом [7]: равномерные сетки обычно имеют более высокий порядок точности аппроксимации по пространству, при достаточно мелкой сетке привязка к объектам становится неактуальной [2]. Первый вариант (23) наиболее эффективен в одномерных расчетах, в двумерном и трехмерном случаях такой подход приводит к тому, что шаг оказывается мелким там, где это не требуется. Второй вариант (24) разрабатывается в последние годы [9-12], он требует существенно более сложных алгоритмов при его реализации. Один из способов состоит в том, что расчетный шаг от ячейки к ячейке последовательно уменьшается в два раза, это может осуществляться автоматически при увеличении градиента температуры. Применяется

также комбинированный подход: равномерная сетка в целом сращивается с мелкой сеткой в определенной области.

В случае равномерных сеток явная схема приобретает более простой вид. Для однородной среды

ТТП+1 _ ТТ" п" + п" _ 2п" и" + и" _ 2и"

п /,у,к п /,у,к _ иг+\,у,к /_1,у,к у,к иЛ,у+1,к ^ ¡-у_1,к у,к +

ХАг Ах2 Ау

+ Ри у к. (25)

и" + и" _ 2и"

/, у ,к+1 + /1, у ,к_1 /, у ,к

Ах2

Для неоднородной среды

" " " "

Л и/_1,у ,к _ Иг,у ,к - Иг+1,у ,к _ Иг,у ,к

И"+1 _ И" _1/2,у,к . + Х/+1/2,у,к .

п,, у ,к у ,к = "_Ах___Ах +

Аг Ах

Х и"!у_1,к _ и"у,к +Х и"у+1,к _ и"у,к Х/, у _1/2,к-:--+ Х, у +1/2,к-:-

+_Ау_Ау +

Ау

Х К у ,к _1 _ И" у ,к и1 у,к +1 _ К у,к

Х/, у ,к _1/2----+ Х/, у,к+1/2---

+-А-:-А-+ , у . (26)

Аг

Здесь для простоты записи численной схемы температуры отнесены к целым узлам, а коэффициенты теплопроводности — к полуцелым.

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

+1/2, у+1/2. к +1/2 = М- (27)

Граничное условие поставлено в полуцелом узле с индексом / + 1/2. Это несущественно, если сетка мелкая или это условие выставлено на дне расчетного объема. При необходимости поставить это условие в целом узле с индексом / применяется другой способ:

_ п., ч Х/_1/2,у+1/2,к+1/2

,у+1/2,к+1/2 = 2(М_ И/_1/2,у+1/2,к+1/2)"

И

1/2

_ п.( Ч Х/+1/2,у+1/2,к+1/2

,у+1/2,к+1/2 = _2(М _ и/+1/2,у + 1/2,к+1/2 ) "

(28)

И

/+1/2

В первом случае расчетный объем слева от поверхности, а во втором — справа.

Граничное условие второго рода в узле с индексом / ставится так:

Ц/,]+1/2,к+1/2 - V, Ц/,]+1/2,к+1/2 - -V . (29)

Однородное граничное условие второго рода в узлах с целым индексом / ставится и через температуру, посредством добавления дополнительных краевых узлов:

+1/2,] +1/2,к+1/2 - Ы/-1/2,] +1/2,к +1/2. (30)

Это условие одновременно является условием симметризации: оно применяется для расчетов, где используются свойства симметрии постановки задачи.

Граничное условие третьего рода в узле с индексом / ставится так:

Ц,]+1/2,к +1/2 -

ч 2

- 2(Д" Ы/-1/2,] +1/2,к+1/2)~-—-Г"—---

И/-1/2 / А/-1/2,]+1/2,к+1/2 + 2(1 / апов + -пов ) _ (31)

Ц,]+1/2,к +1/2 -

2

- -2(д - Ы/+1/2,]■ +1/2,к+1/2)

И/+1/2 / А/+1/2,] +1/2,к +1/2 + 2(1 / а пов + —пов )

Здесь -пов — коэффициент термического сопротивления поверхности; апов — коэффициент конвективного теплообмена.

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

В качестве расчетного объема рассматривается двухфазный нано-материал. Характерные размеры и температура изначально обезраз-мериваются, что позволяет привести задачу к универсальному виду. Обезразмеривание проводится таким образом, что критической температурой (температурой фазового перехода) является 0.

Первый этап расчетов — наблюдение за тем, как изменяется распределение температурных полей в на-номатериале при наличии внешнего поверхностного источника тепла. Пространственные шаги расчета одинаковы по всем осям. Считаем, что в пределах этой области на верхней границе расчетного объема будет иметь место первое граничное условие: постоянно заданная температура. Рассмотрим, как меняется температура в структуре наномате-риала с течением времени. На рис. 8 и 9 представлены графики изменения температуры в точках наблюдения. Точки 1, 2, 3 находятся на оси симметрии расчетного объема, на углу области, подвергающейся воздействию поверхностной температуры, и в области без влияния внешнего тепла соответственно. Точки 4, 5 и 6 выбраны в тех же областях, но, в отличие от точек 1, 2 и 3, они располагаются глубже относительно поверхности наноструктуры.

На рис. 10 отражена цикличность изменения температурного распределения в структуре наномате-риала в случае, когда поверхностная температура меняется во времени.

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

Рис. 8. Изменение температуры в точках наблюдения 1, 2, 3 за определенный интервал времени

Рис. 9. Изменение температуры в точках наблюдения 4, 5, 6 в течение определенного отрезка времени

Рис. 10. Цикличность изменения температуры в точках наблюдения 1, 2, 3 с течением времени

Рис. 11. Изменение температуры в точках наблюдения 1, 2, 3 в течение времени

бесконечно тонкий источник (сток) тепла с постоянной мощностью.

Рассмотрим температурные профили в точках наблюдения. Точки наблюдения взяты аналогично первому этапу исследования. На рис. 11 представлены результаты численного моделирования для температурного распределения в точках 1, 2 и 3 — это точки вблизи поверхности, пространственное распределение которых было указано ранее. Влияние стоков тепла на температуру наноструктуры под областью с постоянным подогревом ярко выражено: если раньше объем прогревался до температуры, подразумевающей фазовый переход, то теперь на оси симметрии температура не превышает критического значения.

Наконец, сравним показатели температуры в расчетном объеме без применения и с применением нанотрубок. На рис. 12 приведены Рис. 12. Показатели температур в температурные профили точек, расточках наблюдения: полагающихся на оси симметрии

1 2 — т°чки на оси симметрии и на угау (точка 1) и на углу расчетного объ-расчетного объема (без учета стоков . тт

тепла); Г, 2' — с учетом стоков тепла ема (точка 2). Черным цветом выделены графики температур в том

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

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

тепла в наноструктуре). На рис. 13 представлено сравнение динамики температур в точках наблюдения 1 и 2 при периодическом цикле работы нанотрубок.

Выводы. На основе теории численных методов и математического моделирования решена задача расчета и прогноза распределения температурного поля в двухфазной нанокомпозиционной среде. Сформулирована математическая постановка задачи в виде интегрального уравнения теплового баланса с учетом теплового потока, изменяющегося по закону Фурье, учтены скачки энтальпии и коэффициента теплопроводности. Рассмотрены разнообразные численные схемы и методы и выбран оптимальный — метод кон -трольного объема. С помощью программного комплекса проведен расчет динамики температурного поля в наноструктуре.

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

ЛИТЕРАТУРА

[1] Белов В.В, Доброхотов С.Ю., Тудоровский Т.Я. Асимптотические решения нерелятивистских уравнений квантовой механики в искривленных нано-трубках. Теор. мат. физ, 2004, т. 141, № 2, с. 267-303.

[2] Калиткин Н.Н. Численные методы. Москва, Наука, 1978, 512 с.

[3] Климов Д.М., Васильев А.А., Лучинин В.В., Мальцев П.П. Перспективы развития микросистемной техники в XXI веке. Нано- и микросистемная техника, 1999, № 1, с. 3-6.

[4] Кособудский И. Д., Ушаков Н.М., Юрков Г.Ю. Введение в химию и физику наноразмерных объектов. Саратов, СГТУ, 2007, 182 с.

[5] Маслов В.П. Комплексный метод ВКБ в нелинейных уравнениях. Москва, Наука, 1994, 384 с.

[6] Маслов В.П., Федорюк М.В. Квазиклассическое приближение для уравнений квантовой механики. Москва, Наука, 1981, 296 с.

Рис. 13. Расчет температуры с горизонтальной нанотрубкой со сдвигом цикла функционирования:

1, 2 — температура в первой и второй точке наблюдения без применения на-нотрубки; 3, 4 — температура в первой и второй точке при применении нано-трубки

[7] Dexter J. Carbon Nanotubes Enable Pumpless Liquid Cooling System for Computers. IEEE Spectrum's nanotecmnology blog. URL: mttp://spectrum.ieee. org/nanoclast/semiconductors/nanotecmnology/carbon-nanotubes-enable-pumpless-liquid-cooling-system-for-computers.

[8] Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. Москва, Эдиториал УРСС, 2009, 784 с.

[9] Сидняев Н.И. Теория фазовых переходов и статистические явления механики наноструктурированных веществ. ВестникМГТУим. Н.Э. Баумана. Спец. выпуск «Наноинжененрия». Приборостроение, 2010, с. 9-22.

[10] Сидняев Н.И., Мельникова Ю.С., Храпов П.В., Гласко А.В. Влияние температурного режима вечномерзлых грунтов на надежность оснований. Мат. Междунар. науч.-практ. конф. по инженерному мерзлотоведению, посвященной ХХ-летию создания ООО НПО «Фундаментстройаркос». Тюмень, Издательство «Сити-пресс», 2011, с. 247-256.

[11] Сидняев Н.И., Храпов П.В., Мельникова Ю.С. Основы математического моделирования распределения температурных полей в многофазных средах. Сб. докладов IV Всероссийской молодежной научно-инновационной школы «Математика и математическое моделирование» 19-22 апреля 2010 г. Саров, ФГБОУ ВПО СарФТИ НИЯУ МИФИ, 2010, с. 85-93.

[12] Сидняев Н.И. Федотов А. А., Мельникова Ю.С. Управление распределением температурных полей в криолитозоне. Academia. Архитектура и строительство. Москва, НИИСФ РААСН, 2010, № 3, с. 372-374.

[13] Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. Москва, Энергоатомиздат, 1984, 152 с.

[14] Сидняев Н.И., Федотов А.А., Крылов Д.А. Интегральный метод в задачах математического моделирования распределения температурных полей. Сб. докладов IVВсероссийской молодежной научно-инновационной школы «Математика и математическое моделирование» 19-22 апреля 2010 г. Саров, ФГБОУ ВПО СарФТИ НИЯУ МИФИ 2010, с. 72-76.

Статья поступила в редакцию 16.07.2013

Ссылку на эту статью просим оформлять следующим образом:

Сидняев Н.И., Ильина Ю.С., Крылов Д.А. Управление температурным полем и его прогнозорование в нанокомпозиционных материалах. Инженерный журнал: наука и инновации, 2013, вып. 6. URL: .http://engjournal.ru/catalog/nano/hidden/812.html

Сидняев Николай Иванович — д-р техн. наук, профессор, заведующий кафедрой «Высшая математика» МГТУ имени Н.Э. Баумана. e-mail: Sidn_ni@mail.ru

Ильина Юлия Сергеевна — аспирантка, ассистент кафедры «Высшая математика» МГТУ им. Н.Э. Баумана. e-mail: Jm.bmstu@yandex.ru

Крылов Дмитрий Алексеевич — аспирант, ассистент кафедры «Высшая математика» МГТУ им. Н.Э. Баумана. e-mail: dmitrykrylov@rambler.ru