Эффективный коэффициент теплопроводности композита при неидеальном контакте шаровых включений и матрицы Текст научной статьи по специальности «Математика»

УДК 541.124

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

ЭФФЕКТИВНЫЙ КОЭФФИЦИЕНТ ТЕПЛОПРОВОДНОСТИ КОМПОЗИТА ПРИ НЕИДЕАЛЬНОМ КОНТАКТЕ ШАРОВЫХ ВКЛЮЧЕНИЙ И МАТРИЦЫ

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

E-mail: zarubin@bmstu.ru, gnk1914@mail.ru, inga_fn2@mail.ru

Ключевые слова: композит с включениями шаровой формы, эффективный коэффициент теплопроводности, наноструктурные элементы

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

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

наружным радиусом R1, окруженный слоем материала матрицы. Примем, что такая составная частица с наружным радиусом R2 является представительным элементом структуры композита и в тепловом отношении взаимодействует с неограниченным массивом однородного материала, коэффициент теплопроводности Л которого подлежит определению как эффективная характеристика композита. Таким образом, модель композита содержит три фазы: включение, слой матрицы и неограниченный массив однородного материала. При этом отношение R3/R3 равно объемной концентрации CV включений в композите. Такая модель формально применима во всем промежутке CV £ [0, 1], но ее использование корректно до таких значений CV < 1, при которых влияние теплового взаимодействия между соседними включениями можно считать мало существенным.

Рассмотрим тепловое взаимодействие отдельно взятой составной частицы и окружающего ее однородного материала, полагая коэффициенты теплопроводности Л1 и Л2 материалов соответственно полого шара и матрицы заданными. Термическое сопротивление между включением и матрицей, характеризующее неидеальный тепловой контакт на разделяющей их поверхности радиусом R1 , является величиной, обратной коэффициенту а контактного теплообмена. При а ^ ж тепловой контакт на этой поверхности становится идеальным. Тепловой контакт на сферической поверхности радиусом R2, отделяющей составную частицу от массива однородного материала, примем идеальным.

Центр полого шара с внутренним радиусом R0 и наружным радиусом R1 < R2 поместим в начале сферической системы координат. Примем, что на большом расстоянии r от начала координат задан вектор градиента температурного поля в однородном материале, направленный по оси сферической системы координат, от которой происходит отсчет угловой координаты 9, т.е. при r ^ ж установившееся распределение температуры в этом материале описывает функция Гте(г, 9) = Gr cos 9, где G — модуль вектора градиента. Эта функция удовлетворяет уравнению Лапласа, которое в сферических координатах имеет вид

1 д ( 2дТ\ 1 д ( . пдТ\ 1 д2Т п ^ Б" r^ sin9— +—-= °. (1)

r2 dr V dr ) + r2 sin 9 д9 V 09 ) ' r2 sin2 9 дф2

В данном случае благодаря параллельности заданного вектора градиента температурного поля и оси отсчета угловой координаты 9 распределение температуры симметрично относительно этой оси и не зависит от угловой координаты ф, т.е. д2T /дф2 = 0.

По мере приближения к составной шаровой частице температурное поле в однородном материале претерпевает возмущение, описываемое

также удовлетворяющим уравнению (1) дополнительным слагаемым AT(r, 0) = (B/r2) cos 0, где B — подлежащий определению постоянный коэффициент. Таким образом, температурное поле в однородном материале, удовлетворяющее заданному условию при r ^ то и уравнению (1), описывает функция

T(r, 0) = Тте(г, 0) + AT(r, 0) = (Gr + B/r2) cos 0. (2)

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

Ti(r,0) = (Air + Bi/r2)cos 0 (3)

и в слое материала матрицы

T2(r, 0) = (A2r + B2/r2) cos 0. (4)

В равенства (2)... (4) входят 5 неизвестных коэффициентов B, A1, B1, A2 и B2, которые необходимо найти из граничных условий на сферических поверхностях с радиусами R0, Ri и R2. При r = R0 из условия отсутствия теплообмена в полости шарового включения с учетом равенства (3) получим

dTi dr или

Ai = 2Bi/R03. (5)

При r = Ri из условия непрерывности плотности теплового потока следует

dT2

= (Ai - 2Bi/R3) cos 0 = 0

r=Ro

A

2"

dT

= ^T2(Ri,0) - Ti(Rb0)) = Ai 1

r=Ri dr

r=Ri

дг

Отсюда с использованием равенств (3) и (4) находим

А2 - 2£2/Д? = - А! + (В - ВО/Я?) = А(А! - 2В^Я?), (6)

где в = аЯ^А2) и А = А1/А2. Наконец, при г = Я2 условия идеального теплового контакта, соответствующие непрерывности не только плотности теплового потока, но и распределения температуры, позволяют записать

' дТ2

ЛдГ

дг

= A

r=R2 2 дг

, Т (Я2,0)= ^2(^2,0).

г=Л 2

После подстановки в эти равенства соотношений (2) и (3) получим

А(Я - 2В/Я?) = А2 - 2В2/Я? и С + В/Я? = А2 + В2/Я3З, (7) где А = А/А2.

Последовательным исключением неизвестных из равенств (5)... (7) находим

2В = 0А(2С1 + С2Оу) - 2(С - С2Оу) Щ3 А(2С1 + С2СУ) + С - С2Су '

где

С1 = А(2+в)(1 - Щ) + в (2 + Щ), С2 = 2А(1 - в)(1 - Щ0) + в (2 + ^2)

и Що = Яо/Я1. в случае включения в виде сплошного шара Щ0 = 0.

Замена составной шаровой частицы равновеликим шаром радиусом Щ2 с искомым коэффициентом теплопроводности А приведет к исчезновению возмущения температурного поля в окружающем ее однородном материале с тем же значением А. Тогда в равенстве (2) следует положить ДТ(г, 9) = 0, что равносильно условию В = 0, которое с учетом формулы (8) позволяет записать

А = 2(С1 - С2Су) (9)

2С1 + С2СУ

В частном случае идеального теплового контакта на сферической поверхности радиусом Щ1, разделяющей включение и матрицу, в ^ ж и равенство (9) при Щ0 = 0 переходит в известную формулу Максвелла [2]

А-2 + Х - 2(1 - Х)Су Л0)

А = 2 + А + (1 - А)Су , (10)

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

А = 2(1 ' = 1 (11)

2 + СУ 2 + СУ' V 7

что соответствует пористому материалу с коэффициентом теплопроводности А2, объемная концентрация пор в котором равна СУ. Равенство (9) переходит в соотношение (11) и в случае абсолютно нетеплопроводных включений (А = 0). Наоборот, при абсолютно теплопроводных включениях (А ^ ж) из равенства (9) получим формулу

- = 2 + в - 2(1 - в)Су

2 + в +(1 - в)Су , ( )

аналогичную формуле (10) Максвелла, но с заменой параметра А на параметр в. Наконец, предельным переходом в правой части формулы (12) при в ^ ж, что соответствует идеальному тепловому контакту между матрицей и идеально теплопроводным включением, приходим

к соотношению

~ 1 + 2Cy _ 3Cv

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

Используем двойственную вариационную формулировку задачи стационарной теплопроводности [3, 4] для получения двусторонних оценок эффективного коэффициента теплопроводности рассматриваемого композита. Область V, содержащую представительный элемент структуры композита в виде половины составной частицы радиусом Д2, выберем в форме прямого цилиндра с достаточно большой площадью S0 параллельных оснований, одно из которых соответствует в сферических координатах значению в = п/2, а точки второй имеют координаты r cos в = H, т.е. высота цилиндра равна H, причем H ^ R2. Боковую поверхность цилиндра примем идеально теплоизолированной, температуру основания при в = п/2 положим равной нулю, а на втором основании зададим температуру GH. Однородный материал в части области вне составной частицы имеет коэффициент теплопроводности Л. Таким образом, в неоднородной цилиндрической области объемом V0 = HS0, ограниченной поверхностью S, распределение температуры T(M) и коэффициент теплопроводности Л(М) являются функциями координат точки M Е V, причем функция Л(М) кусочнопостоянная и принимает значения Л1 при R0 ^ r ^ R, Л2 при R < r < R2 и Л при r > R2.

Для минимизируемого функционала [4]

где V — дифференциальный оператор Гамильтона, а ДТк(0) — разность температур на контактной полусферической поверхности радиусом Дь примем при г ^ Д1 в качестве допустимого линейное по высоте цилиндра распределение температуры с постоянной составляющей градиента С, т.е.

о

T*(г, 0) = Gr cos 0, а в шаровом включении аналогично соотношению (3) T*(r,0) = (Air + Bl/r2) cos 0.

,2

(14)

(15)

Из условия отсутствия теплообмена в полости шарового включения с учетом соотношения (15) получим равенство

А1 = 2В1/Л03, (16)

аналогичное формуле (5). При г = Щ1 из условия непрерывности плотности теплового потока с использованием соотношений (15) и (16) следует

в (С - А1(1 + Щ3/2)) = АА1(1 - Щ3).

Отсюда находим

А = Св (17)

1 А(1 - Щ3) + в(1 + Щ3/2) ( )

и

ДТ (9) СЩ1А(1 - _Щ3)со5 9 (18)

ДТк(9) = А(1 - Щ3) + в(1 + Щ3/2) - (18)

Составляющие градиента температуры в шаровом включении, согласно соотношению (15), равны

дТ * 1 дТ*

д-1 = (А - 2вВ1/г3)со89 и --1 = -(^41 + вВ1/г3)81п9, дг г д9

что позволяет вычислить

(V-*)2 = (А41 + 2А1В1)(1 - 3 соэ2 9)/г3 + В2(1 + 3сой2 9)/г6. (19)

Тогда из формулы (13) с учетом равенств (14), (16), (18) и (19) получим

п2 о-гг г?3 п2 г?3 г?3 п2 г?3 г?3

т Гт, „ С „ „ 2ПЩ „ С Що — Л1 , С Щ1 — Щп

71 [Т ] = А—Я50--г-2 А— + 2п——1А^— + 2п——0 х

2 0 3 2 3 _ 2 2 _ 3

Л Ä2 -3/ч л3 А2в(1 - R)2

х Л1 -2-(1 + л3/2) + 2^А2 (-(1 - -(1^3-i?3/2))

Примем для максимизируемого функционала [4]

I[я] = -2 / ^(м) - IТ(Р)я(Р) • п(Р) ¿5(Р) -

У я

тг/2

- 2пЩ2ДТк(9)вт9^9, Р е 5, (21)

о

где п — единичный вектор внешней нормали к поверхности Б, в качестве допустимого распределения вектора плотности теплового потока я при г ^ Щ1 постоянное значение д = -АС единственной составляющей этого вектора, перпендикулярной основаниям цилиндра.

В шаровом включении скалярный квадрат вектора плотности теплового потока представим с учетом соотношения (19) в виде

q2 = Л?(УТ°)2 =

= (A2 + 2a4iBi)(1 - 3 cos2 e)/r3 + B2(1 + 3cos2 e)/r6, (22)

соответствующим распределению температуры во включении, определяемому формулой

T0(r, в) = (Air + B?i/r2) cos в, (23)

аналогичной формуле (15). Из условия отсутствия теплообмена в полости шарового включения с учетом соотношения (23) находим

Ai = 2Bi/R3, (24)

а при r = R1 из условия непрерывности плотности теплового потока с использованием соотношений (23) и (24) следует

q = -ЛЯ = -Л^А - 2B?i/R3) = -Л1 Aii(l - Rr3)-

Отсюда находим

'i Л G

A = ЛЛ1l-M • (25)

+ + ) + ЛG2HSо• (26)

В данном случае ДТк(в) = (Л/а)С cos в и формула (21) с учетом соотношений (22), (24) и (25) примет вид

(лс)2 /HSo-2nR3/3 R3-R3

Ji [q] =-(—л + "Злт1+

(l + R3/2) + 2 _R|_ ■3Л1(1-^23) + зЛ2в

Принятые в качестве допустимых распределения температуры и плотности теплового потока для неоднородной области отличаются от действительных и поэтому значения J1 [T] и J1 [q] не будут совпадать, причем J1 [T] > I1[q]. В промежутке между этими значениями должно быть расположено и значение J0 = ^/2)G2HS0 минимизируемого функционала (13) для однородной области с коэффициентом теплопроводности Л. Тогда при (R1/R2)3 = CV с учетом формулы (20) из условия J1[T] ^ J0 получим

л < (1 г ) + la Л + /ß(1 + R/2)(T - R) - л

А < (1 - Г) + ^ (Л(1 - R0) + ß(1 + R8/2))' - Л

а при использовании формулы (26) из условия Ii [q] < J0 найдем

А ^ 1 - Cy + (Cy/Л)(1 + R0/2)/(1 - R3) + Cy/ß = Л-.

Рис. 1

На рис. 1 для случая Щ0 = 0 и в = 1 при различных значениях А приведены графики зависимостей от СУ верхней А+ (штрихпунк-тирные линии) и нижней А- (штриховые линии) оценок отношения А = А/А2. Сплошными линиями представлены графики зависимостей А, построенные по формуле (9), которая в случае Щ0 = 0 принимает вид

2А + в(2 + А) - 2(в(1 - А) + А)Су

А =

2Л + в(2 + А) + (в(1 - А) + A)CV

(27)

Из этого рисунка следует, что разность А+ - А- уменьшается по мере увеличения значения А. При А = 10 и А = 100 кривые для А-, отмеченные соответственно светлыми кружками и квадратами, практически совпадают со сплошными кривыми для А и штрихпунктирными кривыми для А+, поскольку различие в значениях А-, А и А+ очень мало. Например, для СУ = 0,5 при А = 10 А- = 0,9524, А = 0,9538 и А+ = 0,9545, а при А = 100 А- = 0,995025, А = 0,995041 и А+ = 0,995050.

С увеличением параметра в тепловой контакт между включениями и матрицей улучшается. Поэтому при прочих равных условиях значения А, А- и А+ возрастают. На рис. 2 представлены зависимости, аналогичные графикам на рис. 1, но при в = 10. Следует отметить, что в случае сплошных шаровых включений (Щ0 = 0) оценки А- и А+, а также правая часть формулы (23) не изменяют своих значений при перестановке значений в и А. Например, графики на рис. 2 для А = 1 идентичны графикам на рис. 1 для А = 10.

Рис.2

Из рис. 1 видно, что при в = 1 ширина полосы, ограниченной кривыми для верхней А+ и нижней А+ оценок, увеличивается по мере уменьшения значения Л несмотря на совпадение значений Л, А+ и А+ при CV = 0 и CV = 1. При в = 10 ширина такой полосы растет по мере отклонения значений Л от единицы. Возможная причина состоит в том, что для таких сочетаний в и Л использованные выше допустимые для функционалов распределения плотности теплового потока и температуры становятся достаточно грубыми. Уточним эти распределения с учетом переменных плотности теплового потока и градиента температуры не только в шаровом включении, но и в шаровом слое матрицы, входящем в представительный элемент структуры композита.

В дополнение к формуле (15), описывающей распределение температуры в шаровом включении, допустимое для минимизируемого функционала (13) распределение температуры в шаровом слое матрицы примем аналогично соотношению (4) в виде

T2* (r, в) = (A2r + B2/r2) cos в. (28)

Для нахождения коэффициентов A2 и B2 и уточненного соотношения для коэффициента Ai в формуле (15) используем с учетом равенства (16) условие непрерывности плотности теплового потока при r = Rl в виде

A2 - 2B2/R3 = в(A2 + B2/R3 - Ai(1 + R3/2)) = АAi(1 - R3) (29)

и граничное условие для температуры при r = R2 в форме

a42 R + B2/R2 = GR2.

(30)

Отсюда получим

A = 3Gb/Z, A = g(2в(1 + R3/2) + Л(2 + в)(1 - R0)) /Z, -2/ЯЗ = с(в(1 + R3/2) + Л(1 - в)(1 - R3))/Z,

где Z = в(1 + R0/2)(2 + Cv) + Л(1 - R0) (2 + в + (1 - в) Су). В рассматриваемом случае с учетом формул (15) и (28)

ДТк(0) = T2* (Ri, 0) - T1*(Ri, 0) =

(D ) 1 _ R>3

A + -0 - Ai(1 + R0/2)) Ri cos 0 = 3СЛ——0R cos 0. (31) Rl^ / Z

Тогда из формулы (13) получим

■ [Т] = Л^г(Я5о - 2пД3/3) + Л2(Д2 + д^) +

+ 2пДЬД Л, Д22 (1 + до/2) + 2пД Л2 в (

что с учетом условия (Т) > </0(Т) = (Л/2)С2Я50 приводит к неравенству

А « ^(а*2+2С^ДАДЦД^+Скв(злЦ^)2 = Л+.

При построении допустимого для максимизируемого функционала (21) распределения плотности теплового потока используем прежние формулы (15) и (28), обозначив в них коэффициенты через А2, В2 и А2, В22 соответственно. Для нахождения этих коэффициентов остаются справедливыми равенства вида (16) и (29), а граничное условие (30) следует заменить на условие q = —ЛС = — Л2(Д — 2В2/Л3) непрерывности нормальной к сферической поверхности радиусом Д2 составляющей вектора плотности теплового потока. В итоге получим

Ах = А = 3в/^ А = Д = 2в(1 + Д3/2) + Л(2 + в )(1 — Д3)

1 СЛ 1' 2 СЛ '

^/Д = В2*/(Л3СЛ) = (в (1 + Д 3/2) + Л(1 — в )(1 — Д3)) /^,

где = 2(1 — ) (в(1 + ДД0/2) + Л(1 — Д*)) + (1 + 2СУ)лв(1 — ДД3) и по аналогии с формулой (16) В2/Д3 = А/2.

Теперь в правой части равенства (31) для разности температур на контактной поверхности при г = Д1 необходимо заменить знаменатель Z на и добавить множитель Л. Тогда из формулы (21) получим

(АС)2 ЯБс-2пД3/3 2 Д3-Д? ( А" В2 ^ 11 и = --А--2п"зАТ^- ЩЩ)(АС) -

- | Щ? ((1 - Я 3)А1 | (1 + № вА"( )2)х

х (АС)2 + АС2ИБо.

Отсюда с учетом условия /1(д) ^ 30(Т) = (А/2)С"ИБ0 находим

Л >_1_= ДО

^ (1 - Су )М + Су Б

где

М=А" + 2СуВ1/Щ\ и Б=Ла4"(1 - Я3)(1 + Щ/2) + в (3А(1 - Щ)/^)".

Непосредственная проверка показывает, что при Су € [0, 1], Щ0 € [0, 1) и любых положительных значениях параметров в и Л значения д+ и Л° совпадают между собой и со значением А, вычисляемым по формуле (9). Это означает, что в рамках использованной математической модели с выбранным представительным элементом структуры рассматриваемого композита формула (9) обеспечивает получение достоверных результатов при таких значениях Су < 1, когда можно пренебречь тепловым взаимодействием между соседними включениями.

Работа выполнена по гранту НШ-255.2012.8 программы государственной поддержки ведущих научных школ.

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

1. К а ц Е. А. Фуллерены, углеродные нанотрубки и нанокластеры. Родословная форм и идей. М.: Изд-во ЛКИ, 2008. - 296 с.

2. К а р с л о у Г., Е г е р Д. Теплопроводность твердых тел: Пер. с англ. - М.: Наука, 1964.-488 с.

3. Зарубин В. С. Инженерные методы решения задач теплопроводности. - М.: Энергоатомиздат, 1983. - 328 с.

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

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