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

УДК 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 — удельная объемная теплоемкость; с — удельная теплоемкость материала; р — плотность; Кхх, Куу —

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

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

Куу + q — 0.

ду

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

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

дТ_

дх

= 0,

AD

дТ

дх

= 0.

BC

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

Але = Tgrn .

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

К

дТ

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

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

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

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

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

Кх

дх

+ К

УУ

дТ

ду

+ 2Х—Т Bt

dV +

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

Si S 22

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

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

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

' Т «) ]

Т

Тк «)

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

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

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

Здесь

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

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

Vе s2

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

^ S2

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

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

[се ] = ХАа^

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

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

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

"2 1 0"

J h[ N][ N]T dS =

1 2 0 0 0 0

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

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

[14].

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

ш

J q[ N f dS =

qLijat

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

J hTGRN [ N]т dS =

hTGRNLijat

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

Лт}

[С]-

dt

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

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

ее е

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

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

At I At

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

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

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

^ + ^2Р2

Pk =-

VI + V

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

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

ст1 + С2т2

Ck =-

т1 + т2

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

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

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

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

¿=1 ]=1

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

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

+

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

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

Q = -л-1-------------------------------------,

АВ

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

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

&АВ

Kyy StAT

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

КХУХ + К2У1

Kyy =-

V + V

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

К1т1 + К 2 т2

K

yy

mx + m2

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

n=1

то

Ё Am cos

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