УДК 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