Нелокальная математическая модель теплопроводности в твердых телах Текст научной статьи по специальности «Математика»

ПРИКЛАДНАЯ МАТЕМАТИКА И МЕТОДЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ

УДК 539.3:536.2

В. С. Зарубин, Г. Н. Кувыркин, И. Ю. Савельева

НЕЛОКАЛЬНАЯ МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ТЕПЛОПРОВОДНОСТИ В ТВЕРДЫХ ТЕЛАХ

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

E-mail: Inga.Savelyeva@gmail.com

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

Успехи в создании новых материалов с использованием нанотех-нологий [1,2] стимулируют разработку и новых математических моделей, описывающих температурные поля в деталях и элементах конструкций из консолидированных наноструктурных материалов. Эти математические модели должны учитывать две существующие противоположные концепции описания структуры любого твердого тела — концепции непрерывности и дискретности [3]. С одной стороны, большинство технических объектов представляются состоящими из сплошного материала и для определения температурных полей в них используют "макроскопические" функции состояния — абсолютную температуру T (x,t), массовую плотность энтропии h (x,t), массовую плотность внутренней (u (x,t)) и свободной (A (x,t)) энергии, вектор плотности теплового потока q (x,t) с составляющими q (x,t) (i = 1, 2, 3) и, возможно, некоторые другие функции, зависящие от координат, заданных радиус-вектором x, и времени t. С другой стороны, для наноразмерных объектов процесс распространения теплоты ассоциируется с такими понятиями, как вектор плотности распределения фононов при фононном (решеточном) механизме распространения теплоты, и аналогичными представлениями об электронном и некоторых других механизмах переноса энергии. Аккумуляцию теплоты ассоциируют с увеличением кинетической энергии колебаний атомов или молекул наноструктурных элементов — возрастанием частоты и амплитуды их колебаний. Поэтому для построения математи-

ческой модели теплопроводности консолидированных наноструктур-ных материалов на макроуровне необходимо объединение двух этих подходов.

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

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

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

¿Т« + А44к = К, (1)

где ¿Т — время релаксации параметра состояния; () = д( )/д£; А44 — параметр, характеризующий размер наночастиц (0 < А44 ^ 1); при увеличении размера наночастиц А44 ^ 1; К - равновесное значение параметра к. Для полимерных материалов имеет смысл учитывать спектр времен релаксации, в случае кристаллических материалов можно ограничиться рассмотрением одного значения времени релаксации [7].

Закон сохранения энергии в предлагаемой модели имеет вид [4, 6]

рт| = — £ + -, , = 1 2 , (2)

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

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

tq к г + Аг^ Кг = ^ ] = 1 3, (3)

где V — время релаксации параметра К; А^ = А.,-г; ёе! А^ > 0; Кг — равновесные значения составляющих вектора К.

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

Яг = 1 Фг3(|ж' - ж|)к(ж',*)^(ж'), (4)

V

где ф^(|ж' — ж|) =0 при (|ж' — ж|) € V' и ф^(|ж' — ж|) = 0 при (|ж' — ж|) € V', ж' € V', V' — объем наночастицы, ж' = ж; V — характерный объем среды.

Далее равновесные значения переменных микроуровня К и Кг положим зависящими от Т (ж, ¿) и дТ (ж, £)/дх следующим образом:

к(ж,*)= / фт(|ж' — ж|)Т(ж', (ж'),

V

dT(x',t) (5)

К (ж, ¿) = —у % (|ж' — ж|) дж/ ¿V (ж'),

V ^

где фт(|ж' — ж|) = 0 и ¿у(|ж' — ж|) = 0 только при (|ж' — ж|) € V' и ж' € V'.

Соотношения, аналогичные (4) и (5), используют в механике деформируемого твердого тела при построении нелокальных зависимостей компонент тензора напряжений от компонент тензора деформации [8, 9]. В общем случае эти соотношения определяют зависимость активных переменных от интегральных (осредненных) значений реактивных переменных, в данном случае Т, к, кг и дТ/дхг. Отметим, что если принять

фгз(|ж' — ж|) = ф° 5(|ж' — ж|), фт(|ж' — ж|) = фТ 5(|ж' — ж|), % (|ж' — ж|) = ¿0.5(|ж' — ж|),

где £(|ж' — ж|) — дельта-функция Дирака, то (4) и (5) переходят в известные соотношения [10].

Если решить систему уравнений (3) при заданном виде Кг из (5), то получим функциональную зависимость составляющих вектора К от

градиента абсолютной температуры в виде

dT (Ж

' £

Kj = Kj | - Z3k(|x' - x|)d ,b) dV(x') I , k = 1, 2, 3.

V

В состоянии, близком к равновесию (t* Kj ~ OJ, из (3) и (5) имеем

dT (X,t) ' i"

Kj (x,t) = - Yjk (|x' - x|) -Ti-JLi dV (ж'), (6)

где У,к(|ж' - ж|) = BJгZгk(|ж' - ж|), Вц = (Лц)-1.

Из (1) при ¿ТК ~ 0 с учетом первого равенства (5) получим

к(ж,Ь) = [ т(|ж' - ж|)Т(ж',Ь)йУ (ж'). (7)

Л44 ]

V

В дальнейшем будем рассматривать только равновесные процессы и сравнительно небольшие отклонения значений абсолютной температуры от температуры Т0 естественного состояния. Тогда, задав Л(Т, к) = Т • Л0(к), получим

дк ЯЛ0 дк ^Л0 дк дк

рТт = -рТ^К Ж * -рТ01К дЬ = ^ (8)

дЛ0

где с = —Т0—--удельная массовая теплоемкость, обусловленная из-

йК

менением вида и границ фононного спектра, т.е. изменением функции распределения частот колебаний атомов. Подставив в (8) равенство (7), получим

РЛ44 / - (|ж' - ж|)^ йУ (ж') = - £ + „.

V

Комбинируя (4) и (6), имеем ?г(ж, Ь) =

= 4>ц (|ж' - ж| )/ Уцк (IX"7' - ж'|) йУ (Х')йУ (ж'). (9)

V V

С учетом вышеизложенного уравнение теплопроводности примет

вид

фТ(|ж- (ж) = = дХ/фу(|х' - *1)/У(X - ж'|)¿V(X'(ж') + ду.

V V

(10)

Уравнение (10) принципиально отличается от известного, оно позволяет рассматривать процесс теплопроводности на макроуровне с учетом процессов микроуровня. При этом интегрально учтено взаимное влияние процессов макро- и микроуровня. Будем полагать, что функции влияния можно всегда представить в виде фт(|ж' - х|) = фтфо(|ж' - ж|) ф»,(|ж' - ж|) = ф°ф1(|ж' - ж|)

%(|ж' - ж|) = Ф2(|х' - ж|) и У,(|ж' - ж|) = ^фз(|ж' - ж|), где фТ, фу, , — константы; фо(|ж' - ж|), ф1(|ж' - ж|), ф2(|ж' - ж|), ф3(|ж' - ж|) — единичные функции влияния такие, что

У фо(|ж' - ж|)^' =1^ ф1(|ж' - ж|)^' =1^ ф2(|ж' - ж|)^' = 1,

V V V '

У фз(|ж' - ж|)^' = 1, (|ж' - ж|) е V', ж' = ж.

V

Отметим, что осреднив подынтегральные выражения в (10) по объему сферы единичного радиуса, получим при JфТ(|ж'-ж|)^(ж') = Н(ж),

V

I фу (|ж' - ж|)^(ж') = фО,Н(ж) и I ,(|ж' - ж|)^(ж') = ,Н(ж),

V V

где Н(ж) — функция Хевисайда, ж е V', результаты, совпадающие с известными [10].

Поскольку тах |ж' - ж| ^ /, где I — характерный размер тела, то Т(ж', ¿) и дТ(ж', £)/дж'к можно разложить в ряд Тейлора в окрестности точки ж:

т, / N . , ! . дТ(ж,*) 1 . . .. . д2Т(ж,*)

Т (ж', *) = Т (ж, *) + (х< - х) + - (х' - хг)(х, - Ху ) +

!( ' _ )( ' _ )( ' _ ) д3Т(ж,£)

+ з! (хг )(Хк Хк) дх^дх,дхк + ... ;

dT (ж' ,t)_ dT (x, t)

dx'k dxk

( ^ )дТ2(ж,^) 1 '_ )( '_ ) д3Т(ж,Ь)

+ (Хг Хг) дХгдХк +2!(Хг Хг)(Хц Х) дХгдХц дХк +

\ , ! \(! \( ! \ д Т(ж, Ь) + 3!(Х г-Хг)(Хц-Хц )(Х к - Хк ) дх гдх, дх к дх 1 + ... .

Уравнение теплопроводности (10) в этом случае принимает следующий вид:

pA-/ (Iх' - xi) |(T <x,t) + (x; - *,)

V

+ 2j(xi - xi)(xj - Xj) дxidx + -J dV(x') =

2| i j

2| д x idxj

д L,ч ^ „i /dT(x,t)i0, ' _ ,dT2(x,t) +

= < |x'-x 1 >/j( ^-x' 1 КПЙГ^-*> dxpdxk

VV

+2(xq - xq)(xp - xp)^ + ...) dV(x'')dV(x') + qv. (11) В нулевом приближении из (11) имеем

PA44 / ^( 1 x'- x 1 )dV <x'> f =

V

^ | A(P/ Ы|x' - x I)dV(x')/ ^a(x'' - x'|)dV(x'')

ik

VV

где л£) = у0 — компоненты тензора теплопроводности, и, наконец, если положить У (| ж' - ж| )йУ (ж') = 1, то уравнение теплопровод-

V

ности принимает известную форму

_дТ д / (т) дТ \

= ^ дХ^ + ^,

где с = с/Л44 эффективная удельная массовая теплоемкость консолидированного наноструктурного материала.

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

^(Iх' - х|) = ^(|х' - х|) =

= - |х' - , |х' - х| < а, р = 0, 1, 2, 3, (12)

где а — характерный размер наноструктурного элемента. Тогда уравнение теплопроводности (10) примет вид:

/ ф(|х' - х|) ^^^^(х') = = Л(Т)ф(|х' - х|)/ ф(|х'' - х'|)^Т^¿V(х'')^(х'), (13)

V V

где Л(т) — теплопроводность тела.

Краевые условия для задачи поверхностного нагрева с учетом (9) запишем следующим образом:

г = 0, Т(х, 0) = То;

х = 0, - Л(т) У ф(|х' - х|)/ф(|х'' - х'|) ¿V(х'')^(х') = дя;

V V

a

х = Ь, - Л(т) / ф(|х' - х|)^ ф(|х'' - х'|) дТ(х''';г) ¿V(х'')^(х') = 0,

V V

(14)

где = ВМгт ехр (-т£Д0), т > 1, т € N.

Применив к уравнению (13) и граничным условиям из (14) конечно-элементную процедуру в форме метода Бубнова-Галеркина, получим систему обыкновенных дифференциальных уравнений первого порядка:

[С] {т} + [К] {Т} = (Т} ,

где [С], [К] — матрицы, характеризующие теплоемкость и теплопроводность исследуемого тела; (Т}, {^Т71 — векторы неизвестных узловых значений температуры и скорости ее изменения; {Т} — вектор тепловой нагрузки. Компоненты вектора {Т} и , матриц

[C], [K] определяются соотношениями

Fp = qs (t)5pi,

Cpq = E / ( / ^(Ix' - x|)Nq(x')dx' | Np(x)dx;

6=1 V (e) \V

V (15)

K =

£ IAIT)( fi>(|*'-z|) f-x'|)iNxfLdx''dx1dNp(xldx,

6-1 V (е) XV V

где Е — число конечных элементов; V(е) — объем конечного элемента; ^р, — одномерные квадратичные функции формы, зависящие от координаты х; р, д — номера узлов сетки, р, д = 1, п;

,г , , (х - х2)(х - х3)

^(х) =-—-, XI < х < Хз;

ЛГ / Ч (х - х2к-1)(х - х2к+1) ^2к (х) =--¿2 ,

x2k-1 < x < x2k+i, k = 1, (n - 1)/2; (x - x2k-l)(x - x2k)

2d2

N2k+1(x) = ^ (x - x2k+2)(x - x2k+3)

2d2

, x2k-1 ^ x ^ x2k+1, , x2k+1 < x < x2k+3,

k = 1, (n - 3)/2;

дт- / \ _ (х хп—2)(х хп—1) , ,

^'п(х) 2^2 , хп—2 ^ х ^ хп;

^ = хр+1 - хр — шаг сетки по пространству, длина конечного элемента в случае использования одномерных квадратичных функций равна 2^.

Отметим, что использование линейных функций формы приведет к исчезновению зависимости элементов Крд от функций влияния <^(|х' - х|) и, соответственно, от характерного размера наноструктур-ного элемента а, что не соответствует известным экспериментальным данным [1, 2].

Значения температуры в узлах расчетной сетки конечных элементов в (к + 1)-й момент времени могут быть получены из решения

системы алгебраических уравнений при аппроксимации {Т^ двухслойной разностной схемой [11]:

([Ск+1] + Д 1(к) [Кк+1]) {Т(к+1)} =д 1(к){^(к+1)} + [Ск+1] {т(к)},

где Д ¿(к) = ¿к+1 — ¿к — шаг по времени.

Для вычисления элементов Си, Кп в (15) рассмотрим два варианта расположения функции влияния относительно конечного элемента (рис. 1). В первом варианте (рис. 1, а) максимальное значение функции влияния достигается в центре наноструктурного элемента, во втором (рис. 1, б) — на границах наноструктурного элемента.

Чтобы оценить влияние характерного размера а на распределение температуры, рассмотрим случаи а1 = й, а2 = й/2, а3 = й/4.

При построении графиков использовали следующие безразмерные параметры и переменные [10]:

z = x/Vato, t = t/to, в = (T - To)/T*,

t * = bto

а = Л(т)/(ре), яа(^) = МГехр(шГ), М = шгп/{ш — 1)!.

На рис. 2 изображены распределения температуры по глубине тела для различных моментов времени (а-в) и расчет температуры поверхности нагреваемого тела (г) при ш = 2 в сравнении с температурой, вычисленной в случае параболического уравнения теплопроводности.

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

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

Рис. 1. Расположение функции влияния относительно конечного элемента:

1, 2, 3 — соответственно а\,а2, а3

Рис. 2. Распределения температуры по глубине тела (а-в) и во времени для нагреваемой поверхности тела (г):

1 — параболическое уравнение теплопроводности; 2 — первый вариант расположения функций влияния относительно функций формы (см. рис. 1,а); 3 — второй вариант (см. рис. 1, б); а — ? = 1, б — ? = 3, в — ? = 5; г — г = 0

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

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

1.Андриевский Р. А., Р а г у л я А. В. Наноструктурные материалы. - М.: Изд. центр. "Академия", 2005. - 192 с.

2. Г у с е в А. И. Наноматериалы, наноструктуры, нанотехнологии. - М.: ФИЗ-МАТЛИТ, 2005.-416 с.

3. Кривцов А. М. Деформирование и разрушение твердых тел с микроструктурой. - М.: ФИЗМАТЛИТ, 2007. - 304 с.

4. Зарубин В. С., Кувыркин Г. Н. Математические модели механики и электродинамики сплошной среды. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2008. -512 с.

5. З а р у б и н В. С., К у в ы р к и н Г. Н. Математическое моделирование термомеханических процессов при интенсивном тепловом воздействии // Теплофизика высоких температур. - 2003. - Т. 4, № 2. - С. 300-309.

6. З а р у б и н В. С., К у в ы р к и н Г. Н. Математические модели термомеханики. - М.: ФИЗМАТЛИТ, 2002. - 168 с.

7. Новик А., Б е р р и Б. Релаксационные явления в кристаллах: Пер. с англ. -М.: Атомиздат, 1975. - 472 с.

8. Введение в микромеханику / Онами М. и др.: Пер. с япон. - М.: Металлургия, 1987.-280 с.

9. К у н и н И. А. Теория упругих сред с микроструктурой. Нелокальная теория упругости. - М.: Наука, 1975. - 416 с.

10. К у в ы р к и н Г. Н., Савельева И. Ю. Математическая модель теплопроводности новых конструкционных материалов // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. - 2010. - № 3. - С. 72-85.

11. Кувыркин Г. Н. Термомеханика деформируемого твердого тела при высокоинтенсивном нагружении. - М.: Изд-во МГТУ, - 1993. - 142 с.

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

Владимир Степанович Зарубин родился в 1933 г., окончил МВТУ им. Н.Э. Баумана в 1957 г. Д-р техн. наук, профессор кафедры "Прикладная математика" МГТУ им. Н.Э. Баумана. Заслуженный деятель науки и техники РФ, лауреат премии Правительства РФ в области науки и техники. Автор более 250 научных работ в области математического моделирования термомеханических процессов в материалах и элементах конструкций.

V.S. Zarubin (b. 1933) graduated from the Bauman Moscow Higher Technical School in 1957. D. Sc. (Eng.), professor of "Applied Mathematics" department of the Bauman Moscow State Technical University. Honored Science and Technology Worker of the Russian Federation, Laureate of RF Government Prize in Science and Technology. Author of more than 250 publications in the field of mathematical simulation of thermomechanical processes in materials and construction members.

Георгий Николаевич Кувыркин родился в 1946 г. Окончил МВТУ им. Н.Э. Баумана в 1970 г. Д-р техн. наук, профессор, зав. кафедрой "Прикладная математика" МГТУ им. Н.Э. Баумана. Лауреат премии Правительства РФ в области науки и техники. Автор более 160 научных работ в области прикладной математики и математического моделирования термомеханических процессов в материалах и элементах конструкций.

G.N. Kuvyrkin (b. 1946) graduated from Bauman Moscow Higher Technical School in 1970. D. Sc. (Eng.), professor, head of "Applied Mathematics" department of the Bauman Moscow State Technical University. Laureate of RF Government Prize in Science and Technology. Author of more than 160 publications in the field of applied mathematics and mathematical simulation of thermomechanical processes in materials and construction members.

Инга Юрьевна Савельева родилась в 1985 г., окончила в 2008 г. МГТУ им. Н.Э. Баумана. Аспиранка кафедры "Прикладная математика" МГТУ им. Н.Э. Баумана. Автор 6 научных работ в области моделирования нестационарной теплопроводности.

I.Yu. Savelyeva (b. 1985) graduated from Bauman Moscow State Technical University in 2008. Post-graduate of "Applied mathematics" department of Bauman Moscow State Technical University. Author of 6 publications in the field of mathematical simulation of non-stationary heat conduction.