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

УДК 536.2

Определение эффективных теплофизических характеристик композиционного материала

П.А. Люкшин, Б.А. Люкшин1, Н.Ю. Матолыгина, С.В. Панин

Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия 1 Томский университет систем управления и радиоэлектроники, Томск, 634050, Россия

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

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

Determination of effective thermophysical characteristics of a composite material

P.A. Ljukshin, B.A. Ljukshin1, N.Yu. Matolygina, and S.V. Panin

Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia 1 Tomsk State University of Control Systems and Radioelectronics, Tomsk, 634050, Russia

In the paper we propose an approach to determining effective thermophysical quantities for a particle reinforced composite. This approach is similar in many respects to a method of determining strain and strength characteristics for such materials. In both cases the task on the detailed distribution of state parameters in the calculation area (over a heterogeneous medium simulating the composite) is solved. For thermophysical characteristics this is temperature distribution and for strain and strength ones the distribution of parameters of the stress-strain state, namely, displacements, strains, and stresses. Different averaging procedures are further performed and the parameters averaged over a heterogeneous calculation area are compared with the corresponding data for a conditionally homogeneous medium. This allows effective characteristics to be estimated. Note that the calculation of the heat conductivity coefficient according to the mixture theory leads to large errors.

Keywords: nonstationary problem of thermal conductivity, finite-element method, Dirichlet boundary conditions, implicit difference scheme, isotherms, thermal conductivity coefficient of a composite material

В работах авторов [1-7] и других исследователей, например [8-12], обсуждался процесс получения эффективных деформационно-прочностных характеристик, когда на основе анализа параметров напряженно-деформированного состояния (перемещений, деформаций и напряжений) в расчетной области отыскивается смещение границы расчетной области (при заданных нагрузках) или распределение напряжений (при задан-

1. Введение

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

Аналогичный подход предлагается использовать для определения таких теплофизических характеристик, как

© Люкшин П.А., Люкшин Б.А., Матолыгина Н.Ю., Панин С.В., 2008

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

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

2. Расчет температурного поля

Задача о распределении температуры в различные моменты времени решается в плоской постановке. Для нахождения температурного поля в прямоугольной расчетной области (пластинке) ABCD (рис. 1) решается двумерная нестационарная задача теплопроводности на основе уравнения [13]

.ЪТ — к д 2Т + к д 2Т

X---— Кхх — + К ------—,

Эt ” дх2 ^ ду2

X = ср,

где X — удельная объемная теплоемкость; с — удельная теплоемкость материала; р — плотность; Кхх, Куу —

(1)

коэффициенты теплопроводности в соответствующих направлениях.

На линии DC задан поток тепла #, а соответствующие граничные условия записываются в виде:

Куу + q — 0. (2)

ду

В некоторых примерах на линии DC задавалась постоянная температура (условие Дирихле).

На линиях AD и ВС задаются условия симметрии (так называемые естественные граничные условия [15]):

дТ_

дх

= 0,

AD

дТ

дх

= 0.

(3)

BC

На линии АВ задается постоянная температура (условие Дирихле), равная температуре TCRN окружающей среды:

Але = Tgrn . (4)

В некоторых примерах на линии АВ задавались условия, определяющие конвективный обмен тепла с окружающей средой:

К

дТ

уу ду + h(T TGRN) = 0,

где И — коэффициент теплообмена.

В качестве начального условия задается поле температуры во всей расчетной области в момент времени t = 0:

Т (х, у, 0) — Т 0( х, у).

Предлагается решение задачи теплопроводности методом конечных элементов. В этом случае решение системы, состоящей из уравнения (1) и граничных условий (2)-(4), сводится к минимизации функционала [14]:

Кх

дх

+ К

УУ

дТ

ду

+ 2Х—Т Bt

dV +

+ J qTdS + J -(T - TGRN )2ds,

Si S 22

(5)

где ^ — площадь поверхности, на которой задан поток тепла; — площадь поверхности, где происходит конвективный обмен тепла.

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

Т(е ) = N (х, у )Т ^) или в матричном виде:

' Т «) ]

Т(е) = [ N (х, у) Nj (х, у) N. (х, у )]| Т} (І)

Тк «)

(6)

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

Рис. 1. Схема расчетной области

[с" ] ^+[ Ке ]{Т} + {fe} = 0. дt

Здесь

[ce ] = JX[ N][N]r dV,

(8)

[ Ке ] = | [Ве ][ De ][Ве ]Т йГ + | Н[ Ж][ Щт (9)

Vе s2

{П = I ч1 ЩТ- IНТС^ЕМ[М]ТdS, (10)

^ S2

где Vе—объем конечного элемента; [Л/] — матрица, которая содержит функции формы; [ Ве ] — матрица, которая содержит производные от функции формы; [ De ] — матрица свойств материала, содержащая коэффициенты теплопроводности.

Интеграл (8) при решении двумерной задачи теплопроводности (матрица демпфирования) равен "2 1 1"

(11)

[се ] = ХАа^

12

1 2 1 1 1 2

где А — площадь конечного элемента; а( — толщина элемента.

Матрица «теплопроводности» в двумерном случае (при отсутствии конвекции) имеет вид:

[ Ке ] =

kxxat

4 А

kyyat

4 А

bibi ib j k bi bi

j b b b j b

bkbi b k b bkbk _

cici cicj cick

cjci cjcj c jck

ckci ckcj ckck

(12)

где Ьг- = Уj - ук; с{ = Xj - хк; остальные величины получаются круговой перестановкой индексов i, j, k.

Второй интеграл в (9), если конвекции подвержена сторона i, j конечного элемента, равен

"2 1 0"

J h[ N][ N]T dS =

1 2 0 0 0 0

(13)

где Lij — длина стороны элемента между узлами і, j.

Если тепло теряется путем конвекции между сторонами с узлами j, k или ^ і, то матрица (13) изменяется

[14].

Если источники тепла внутри пластины отсутствуют, а приток тепла осуществляется в виде теплового потока q на стороне элемента с узлами і, j, то «вектор нагрузки» элемента (интеграл (10)) равен

ш

J q[ N f dS =

qLijat

(14)

Если на стороне с узлами і, j происходит конвективный обмен тепла, то «вектор нагрузки» при решении задачи теплопроводности равен

J hTGRN [ N]т dS =

hTGRNLijat

Для сетки конечных элементов записывается система обыкновенных дифференциальных уравнений:

Лт}

[С]-

dt

- + [K ]{T} + {F} = 0,

(15)

где [С] = £[ее]; [К] = £[Ке]; [F] = £[Г].

ее е

Заменяя производную по времени в уравнении (15) ее конечно-разностным аналогом, получим неявную разностную схему для решения уравнения теплопроводности методом конечных элементов [15, 16]:

— + [K] ){T}”+1 = ^-{T}” - {F}”+1.

At I At

(16)

Таким образом, если известен вектор температуры {Т}и в момент времени ьп, то температура пластины в момент времени Ьп+1 = ьп + Дь получается в результате решения системы линейных алгебраических уравнений (16) методом Г аусса. В данном примере сетка содержала 7 260 конечных элементов и 3751 узел. Число уравнений в системе при решении задачи теплопроводности равно 3751.

3. Определение плотности и удельной теплоемкости композиционного материала

Плотность композиционного материала рассчитывается по формуле теории смесей

^ + ^2Р2

Pk =-

(17)

VI + V

где V1 иV2 — объем первой и второй фазы; р1 и р2 — плотность материалов первой и второй фазы.

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

ст1 + С2т2

(18)

Ck =-

т1 + т2

где т1 и т2 — масса первой и второй фаз соответственно; с1 и с2 — удельная теплоемкость первой и второй фаз.

4. Определение коэффициента теплопроводности

Коэффициент теплопроводности рассчитывается следующим образом. Количество теплоты, передаваемое пластинке АBCD через сторону DC, равно

Q = £ С^Р^Т + £ С2^-Р2ДТ;, (19)

¿=1 ]=1

где П1 и П2 — число конечных элементов первой и второй фаз (п1 + п2 = 7260); ДТ-, ДTj — изменение температуры в i или j конечном элементе.

Изменение температуры в каком-либо элементе пластинки равно разности между температурой Т(х, у, (),

+

полученной в результате решения задачи теплопроводности, и начальной температурой Т(х, у, 0).

То же самое количество теплоты можно получить исходя из того, что композиционный материал представляет собой однородный материал с некоторым осред-ненным коэффициентом теплопроводности Куу [17]: К БьДТ

Q = -л-1-------------------------------------, (20)

АВ

где £ — площадь стороны пластинки DC; t — время, в течение которого в пластинку передается количество теплоты Q; АТ — разность между температурой на стороне DC и начальной температурой пластинки; 1АВ — расстояние, на которое распространилось тепло за время t.

Если в (20) задаться конкретными значениями всех величин, кроме коэффициента теплопроводности, то легко видеть, что коэффициент теплопроводности определяется соотношением

&АВ

Kyy StAT

(21)

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

КХУХ + К2У1

Kyy =-

V + V

(22)

где К1 и К2 — коэффициенты теплопроводности пер вой и второй фаз, либо формулой, где фигурируют мас сы фаз:

К1т1 + К 2 т2

K

yy

mx + m2

(23)

5. Тестовая задача

На рис. 2 показана пластинка из неоднородного материала вместе с нанесенной на ней сеткой конечных элементов. Включения в пластинке показаны жирными линиями.

На рис. 3 приведено сравнение численного и аналитического решений нестационарной задачи теплопро-

Рис. 2. Сетка конечных элементов, нанесенная на неоднородную пластинку, состоящую из двух различных материалов

Рис. 3. Сравнение численного (а) и аналитического (б) решений нестационарного уравнения теплопроводности. На кромках АВ и АБ задаются условия Дирихле (Т = 0), на СВ и БС—условия симметрии дТ/дх|св = 0, дТ/ду|^с = 0- Начальная температура пластинки Т(х,у, 0) = 1

водности для однородной пластинки. Начальная температура в расчетной области (прямоугольной пластинке) равна Т(х, у, 0) = 1. На кромках АВ и АБ температура равна 0. На кромках ВС и СБ заданы условия симметрии.

Аналитическое решение нестационарного двумерного уравнения теплопроводности при вышеприведенных граничных условиях задается следующим выражением [13]:

Т(х, у, t) =

4(-1) «+1 4(-Г П

гда * = (2ГГП; * = ¿ГГ*; =(2n-»У

п K

Vm = (2m -1)—; a = HL; n = 1, 2, 3, ...; m = 1, 2, 3, ...; 2 cp

R — длина пластинки в направлении x; L — длина пластинки в направлении у; a — коэффициент температуропроводности.

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

Удельная теплоемкость материала пластинки С = = 460 Дж/(кг • K), удельная плотность p = 7 800 кг/м3, ко-(24) эффициент теплопроводности Kxx = Kyy = 60 Вт/(м • K).

\m+1

Ё An cos(^nxR) exp(- цlat/R2 )

n=1

то

Ё Am cos(Vmy/L)exp(-v;maV l2)

m=1

t = 40 С

t = 60 С

t = 80 с

X

Рис. 4. Изотермы в неоднородной пластинке (слева) и изотермы в пластинке из эквивалентного однородного материала (справа), соответствующие различным моментам времени t = 40, 60, 80 с

Время, в течение которого идет процесс теплопроводности, равно 800 с. Для этого момента времени на рис. 3 построены поверхности температуры и изотермы.

6. Расчет коэффициента теплопроводности композиционного материала

На рис. 4 и 5 приводятся изотермы в пластинке из двухфазного материала и изотермы в пластинке из однородного материала для различных моментов времени с теплофизическими характеристиками, вычисленными по формулам (17), (18), (21). При нахождении коэффициента теплопроводности по формуле (21) количество

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

0.000 0.005 0.010 0.000 0.005 0.010

0.000 0.005 0.010 0.000 0.005 0.010

D CD С

0.000 0.005 0.010 0.000 0.005 0.010

А В А В

Рис. 5. Изотермы в неоднородной двухфазной пластинке (слева) и изотермы в пластинке из однородного материала (справа), соответствующие различным моментам времени t = 5, 10, 15 с

0.010

0.005 '

0.000

0.000

А

0.005

C

0.010

в

0.000

0.000

A

Рис. 6. Изолинии и поверхности температуры в неоднородной двухфазной пластинке (а), в пластинке из однородного материала с Куу = = 0.38 (б) и 22 Вт/(м • К) (в). Количество теплоты в неоднородной двухфазной пластинке Q = 0.0057 Дж (а), в пластинке из однородного материала Q = 0.0068 (б) и 0.0198 Дж (в), t = 15 с

Характеристики эквивалентного однородного материала рассчитаны по формулам (17), (18), (21). В результате, пластинка из этого материала имеет следующие характеристики: плотность р к = 3690 кг/м3, удельная теплоемкость Ск = 643 Дж/(кг • К), коэффициент теплопроводности Куу = 0.38 Вт/(м • К). Изотермы в неоднородной пластинке и в пластинке из однородного материала получены из решения задачи теплопроводности (решения дифференциального уравнения в частных производных с граничными условиями Дирихле и Неймана). На кромках АВ и СБ ставятся условия Дирихле (Т = 0; Т = 1), на кромках АБ и ВС ставятся условия симметрии ЭТ/ Эх = 0.

В расчетах принималось, что неоднородная пластинка состоит из двух материалов: железа и полиэтилена. Для железа принимались следующие значения физических характеристик: удельная теплоемкость С = = 460 Дж/(кг • К), плотность р = 7800 кг/м3, коэффициент теплопроводности Куу = 60 Вт/(м • К). Характеристики полиэтилена принимались следующими: удельная теплоемкость С = 1257 Дж/(кг • К), плотность р = 1333 кг/м3, коэффициент теплопроводности Куу = 0.14 Вт/(м • К).

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

теристиками [13]. Изотермы в неоднородной пластинке и в пластинке из однородного материала получены из решения нестационарной задачи теплопроводности (решения дифференциального уравнения в частных производных с граничными условиями Дирихле и Неймана).

На рис. 6 приведены изолинии и поверхности температуры в пластинке из двухфазного материала (рис. 6, а), из однородного материала с характеристиками, вычисленными по формулам (17), (18), (21) (рис. 6, б), из композиционного материала с характеристиками, вычисленными по формулам (17), (18), (22) (рис. 6, в). На кромках АВ и БС задавались условия Дирихле, на кромках АБ и ВС задавались условия симметрии. Нетрудно заметить, что поверхности, характеризующие распределения температуры, и изотермы на рис. 6, а и б весьма близки, а отличия поверхности распределения температуры и изотерм на рис. 6, в от приведенных на рис. 6, а и б имеют принципиальный, качественный характер.

7. Выводы

Из полученных результатов расчетов следует, что скорости распространения тепла в неоднородной пластинке из двухфазного материала и в пластинке из однородного материала с характеристиками, вычисленными по формулам (17), (18), (21), примерно равны. Скорость распространения тепла в пластинке с коэффициентом

теплопроводности, вычисленным по теории смесей по формуле (22), превышает скорость распространения тепла в неоднородной пластинке практически втрое.

Количество теплоты в неоднородной пластинке и в пластинке из однородного материала, характеристики которого рассчитаны по формулам (17), (18), (21), отличается на 20 %. Количество теплоты в неоднородной пластинке и в пластинке из однородного материала с характеристиками, рассчитанными по формулам (17), (18), (22), отличается в 3.5 раза. Следовательно, вычисление коэффициента теплопроводности композиционного материала по теории смеси по формуле (22) или (23) приводит к большим погрешностям.

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

Приведенный пример свидетельствует о том, что при определении эффективных теплофизических характеристик, как, впрочем, и деформационно-прочностных, теория смесей совершенно непригодна. Так, по теории смесей расчет по объемному содержанию включений (около 50 %) дает значение эффективного коэффициента теплопроводности примерно равным Куу = 30 Вт/(м • К); если же взять за основу массовое содержание, оценка будет еще выше. Ни в том, ни в другом случае получаемые значения не соответствуют реальности. Предлагаемый подход к определению эффективного коэффициента теплопроводности в этом отношении представляется вполне оправданным.

Литература

1. Люкшин Б.А., Люкшин П.А. Влияние свойств межфазного слоя на

напряженно-деформированное состояние полимерной композиции в окрестности включения // Механика композиционных материалов и конструкций. - 1998. - Т. 4. - № 2. - С. 56-68.

2. Анисимов И.И., Десятых В.И., Люкшин Б.А., Люкшин П.А. Форми-

рование прочностных характеристик наполненных полимерных

систем на мезоуровне // Механика композиционных материалов и конструкций. - 1998. - Т. 4. - № 4. - С. 74-92.

3. Люкшин Б.А., Люкшин П.А. Прочностной анализ дисперсно-напол-

ненных полимерных систем на мезоуровне // Физ. мезомех. -

1999. - Т. 2. - № 1-2. - С. 57-67.

4. Алексеев Л.А., Гузеев В.В., Липовка М.В., Люкшин Б.А., Люкшин П.А., Матолыггина Н.Ю. Опыт прочностного конструирования наполненной полимерной композиции // Физ. мезомех. -

2000. - Т. 3. - № 1. - С. 59-66.

5. Люкшин Б.А., Люкшин П.А., Матолыггина Н.Ю. Двухэтапный про-

цесс компьютерного конструирования наполненной полимерной композиции // Физ. мезомех. - 2000. - Т. 3. - № 4. - С. 71-77.

6. Дашук ИВ., Люкшин Б.А., Люкшин П.А., Матолыггина Н.Ю. Влия-

ние деформационно-прочностных свойств структурных элементов на характеристики дисперсно-наполненных композиций // Механика композиционных материалов и конструкций. - 2004. - Т. 10.-№ 3. - С. 366-384.

7. Анисимов И.И., Бочкарева С.А., Десятыгх В.И., Люкшин Б.А., Люк-

шин П.А., Матолыггина Н.Ю., Смолянинова Н.В. Эффективные деформационно-прочностные характеристики полимерной композиции с дисперсными включениями разных размеров // Физ. мезо-мех. - 2006. - Т. 9. - № 2. - С. 11-15.

8. Килина О.В., Килин П.С., Кулыков С.Н. Моделирование деформационного поведения пористой керамики // Физ. мезомех. - 2002. -Т. 5. - № 4. - С. 47-54.

9. Сидоренко Ю.Н., Шевченко Н.А. Прогнозирование механических свойств биометаллического материала на основе многоуровневой математической модели // Физ. мезомех. - 1999. - Т. 2. - № 1-2. -С. 37-42.

10. Макаров П.В. Подход физической мезомеханики к моделированию процессов деформации и разрушения // Физ. мезомех. -1998.- Т. 1. - № 1. - С. 61-81.

11. МакаровП.В., БекетовК.А., Атаманов О.А., Кулыков С.Н. Вязкая конструкционная керамика: моделирование эволюции мезострук-туры под нагрузкой // Физическая мезомеханика и компьютерное конструирование материалов. - Новосибирск: Наука, 1995. -Т.2.- С. 153-171.

12. Psakhie S.G., Korostelev S.Yu., Negreskul S.I., Zolnikov K.P., Wang Z., Li S. Vortex mechanism of plastic deformation of grain boundaries. Computer simulation // Phys. Stat. Sol. B. - 1993. - V. 176. - P. K41-K44.

13. Плят Ш.Н. Расчеты температурных полей бетонных сооружений. - М.: Энергия, 1974. - 408 с.

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

15. Зенкевич О., Морган К. Конечные элементы и аппроксимация. -М.: Мир, 1986. - 316 с.

16. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975. - 541 с.

17. КухлингX. Справочник по физике. - М.: Мир, 1982. - 520 с.

Поступила в редакцию 11.03.2008 г., после переработки 21.04.2008 г.

Сведения об авторах

Люкшин Петр Александрович, к.ф.-м.н., старший научный сотрудник ИФПМ СО РАН, petrljuk@ispms.tsc.ru

Люкшин Борис Александрович, д.т.н., профессор, зав. кафедрой механики, графики и управления качеством ТУСУР, borisljuk@mail.ru Матолыгина Наталья Юрьевна, к.ф.-м.н., научный сотрудник ИФПМ СО РАН, ksa@ispms.tsc.ru

Панин Сергей Викторович, д.т.н., зав. лабораторией полимерных композиционных материалов ИФПМ СО РАН, svp@ispms.tsc.ru