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

ФИЗИКА

УДК 539.411

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

А. В. Станиславович

ON SOME PROBLEMS OF THERMAL CONDUCTIVITY OF UNIDIRECTIONALL Y REINFORCED MATERIALS

A. V. Stanislavovich

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

Structural model of thermal conductivity of unidirectionally reinforced layer is under consideration. The assumption that temperatures are equal within represented element of reinforced layer is addressed in the present research. In order to answer that question the heat equation with discontinuous coefficients on the phase boundary is solved. On phase boundaries perfect thermal contact conditions are set. Then differential problem is approximated by difference scheme. Effective difference schemes have been used. To get numerical approximation on real-life data fiberglass reinforced thermoplastics are under consideration. Numerical results show that temperature in composition phases equals quite fast (about 10-3 seconds), which is significantly less than the typical time of service of constructions made of composite materials. This allows assume that temperatures are equal within the represented element of reinforced layer.

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

Keywords: thermal conductivity, reinforced materials, partial differential equations, numerical methods, difference schemes.

Введение

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

характера которых необходимо для всестороннего анализа несущей способности конструкции.

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

При выводе определяющих уравнений термоупругого поведения армированного слоя предполагают, что температура в связующем и армирующих элементах в пределах представительного элемента армированного слоя одинакова, например [1; 2; 6]. Однако вопрос о том, как быстро выравнивается температура в пределах представительного элемента армированного слоя при тепловом нагружении композитного материала волокнистой структуры, остается открытым. Исследованию данного вопроса посвящена настоящая работа.

150

А. В. Станиславович

ФИЗИКА

Рис. 1. Представительный элемент армированного слоя

Рассматривая модель однонаправленно-армированного слоя, описанную в [1; 2], в качестве представительного элемента однонаправленно-армированного слоя будем рассматривать параллелепипед ABCDD'C'B'A' (рис. 1). На рис. 1 заштрихованной области соответствует армирующее волокно.

В настоящей работе, также как и в [1; 2], при построении математической модели принимается допущение о прямоугольной форме сечения армирующих волокон.

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

Стеклопластики - композитные материалы, в которых армирующим волокном выступает стекловолокно, а связующим - высокотемпературные пластики (полиимид, полиэфирсульфон, полисульфон и другие) [5; 9; 10].

Существующие в настоящее время значения параметров формирования композитных материалов таковы: характерный диаметр нити стекловолокна порядка 13 мкм [10], массовая доля (Vm) стекловолокна в композите составляет 15 - 50 % [10]. Отсюда следуют характерные размеры представительного элемента (достаточно малого элемента, чтобы нельзя было рассматривать его как квазиоднородную среду) композитного материала.

Физико-механические характеристики (плотность, теплопроводность, теплоемкость) рассматриваемых материалов представлены в таблице 1 [10; 11; 12; 13].

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

Таблица 1

Физико-механические характеристика материалов

кг Плотность, —- м Дж Теплоемкость, кг х С ° Вт Теплопроводность, мх С°

Стекловолокно 2550 800 1,2

Полиимид 1420 1090 0,12

Полиэфирсульфон 1370 1100 0,18

Полисульфон 1240 1100 0,26

1. Одномерная задача

Постановка задачи

Задача в одномерной постановке является модельной и ее можно рассматривать в качестве первого приближения.

Пусть G = |0 < X < lx} - отрезок длиной lx,

Г = | X = 0} U {X = lx } - его граница. Тогда

G = G + Г (рисунок 2). Также введем следующие области:

0 li h lx X

• • • • ►

Рис. 2. Одномерная постановка задачи

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1 151

ФИЗИКА

G1 = {0 < х < /х}и {/2 < x < lx }, G2 = {/1 < Х < l2 },

Г1 ={х = l1}, Г2 ={х = l2}-Тогда G* = G1 + G2,

G = G1 + G2 + Г1 + Г2.

Рассматривается одномерное уравнение теплопроводности на отрезке, составленном из трех отрезков с различными теплофизическими характеристиками [4, с. 17 - 21]:

рС = Ци, х е G*, t е(0, T ],

г 'и

Ци = А------^.

1 rs-,2

дх2

Здесь X - коэффициент теплопроводности, р - плотность, С - удельная теплоемкость. Параметры X, р, С терпят разрывы на границах раздела фаз (области Г1, Г2):

А = Р =

А; хе G1

А.2; Х е G2

р; х е G1

р2; х е G2

С =

Q; х eGx

\G2

[С2; х<

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

- ди = - ди

/Li /L«л 1—'

" дх 2 дх при хеГ1,

lim и = lim и

х——/—0 х——/1+0 ,

ди = . ди дх дх

при

х е Г2

lim и = lim и.

х——/2—0 х——/2+0

Начальные и граничные условия задаются в следующем виде:

4=0 = и( х, y,0) = T0, (1.1)

и\Г = и(0, t) = и (/х, t) =

T +

1—-

_1____

(1 +1051 )5

(1.2)

Подобный выбор граничного условия (1.2) позволяет удовлетворить начальному условию (1.1).

Ставится задача определить поле температур в области G и, в частности, рассчитать время T выравнивания температуры в области G.

Метод решения

Зададим Nx - число разбиений по пространственной координате х. В области G построим равномерную сетку G)h с шагами h = /х / Ых.

О ={(G), х = ihx, 0 <i < Ых}.

Дифференциальный оператор L1u аппроксимируем следующим разностным оператором [3, с. 381]:

ЛЛ =

X,—12 ( Ui—1 Ui )

—Xi+1/2 (иi — Ui+1 )

(1.3)

Здесь Xi—1j2 , Xi+1/2 - разностные коэффициенты

теплопроводности. Так как коэффициент теплопроводности X непрерывен в G всюду, кроме точек разрыва, то положим:

X+1/2 Ai +12

А + А,

X—1/2 определяется аналогично. Под узловыми

значениями разрывных величин надо понимать соответствующие односторонние пределы [3, с. 382].

Производную по времени аппроксимируем следующим разностным выражением [7, с. 18]:

ди uk+1 — и

____ i______f_

.

дt т

Пусть т - шаг по времени, тогда t^ = кт - точки разбиения по времени. Будем считать, что

и *и (х, tk).

Для решения разностной задачи используем разностную схему, которая при использовании разностных операторов вида (1.3) называется «наилучшей» схемой [3, с. 380 - 384; 7, с. 185 - 201; 8, с. 385 - 389]. Данная схема позволяет удовлетворить условиям идеального теплового контакта, она абсолютно устойчива и имеет второй порядок аппроксимации по временной и пространственным координатам даже в случае разрывного коэффициента теплопроводности:

O(T2 + h2) [3, с. 382 - 384; 8, с. 388].

Сик+1—ик Л

РгСг~----- = Л1и

т

к+1

1 < i < Nx — 1. (1.4)

Здесь р, с, - значения плотности и теплоемкости в соответствующих узлах сетки. В точках разрыва их значения определяются по формулам:

= р+1 + р—1 С = Ci+1 + Ci—1

2 ’ i 2 '

Сетка строится таким образом, чтобы точки разрыва приходились на узлы сетки (Oh .

Обозначим:

Р

152 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

M = T +

1—

J_____

(1 +105 kr)5

Ti.

Граничные условия получаются из (1.2):

f h+1 , к+1

J uo = M

1 h+1 , к+1 .

I UNx = M

(1.5)

На каждом временном шаге получается система Nx -1 линейных алгебраических уравнений с трех-

диагональной матрицей. Вид этой системы может быть легко получен путем анализа выражений (1.4) -(1.5). Система решается методом прогонки [8, с. 35 -38]. Численное интегрирование ведется пока не будет выполнено условие остановки итерационного процесса.

Численные результаты

Введем линейную интенсивность армирования (рисунок 2):

W =

l_2_____lj_

1

Можно показать, что зависимость интенсивности армирования w от массовой доли армирующего волокна Vm выражается в следующем виде:

W =

VmPc Ра + Vm {Pc

Здесь рс, ра - плотность связующего и армирующего соответственно.

В таблицах 1 - 3 представлено характерное время выравнивания температуры в области G в зависимости от использующегося связующего и массовой доли

стекловолокна в композите. Здесь области G1 соот-

ветствует связующее, а области

G2 - армирующее

стекловолокно.

В качестве критерия остановки процесса принималось условие:

dT =

max(uf- ul г)

1<г <Nx -1

<6 .

Расчёты проводились при следующих начальных

и граничных условиях: T = 100С°, Т0 = 0C°,

6 = 10-7 . Число узлов сетки Nx = 100, шаг по

времени Г = 10 6 с. Таким образом Г ~ h, что хорошо согласуется с литературой [3; 8].

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

ФИЗИКА

Таблица 1

Связующее полиимид___________

Vm 1 10-6 ч м l2 l1 У 10-6 d w = — l T, 10-3 с.

0,15 145 13 0,09 385,8

0,3 67 13 0,19 81,7

0,4 48 13 0,27 41,2

0,5 36 13 0,36 22,2

Таблица 2

Связующее полиэфирсульфон________

Vm 1 10-6 b м l2 l1 У 10-6 d w = — l T, 10-3 с.

0,15 150 13 0,09 268

0,3 69 13 0,19 56,3

0,4 49 13 0,26 27,9

0,5 37 13 0,35 15,5

Таблица 3

Связующее полисульфон__________

Vm l, 10-6 м l2 l1 У 10-6 d w = — l T, 10-3 с.

0,15 164 13 0,08 200,8

0,3 75 13 0,17 41,7

0,4 53 13 0,24 20,6

0,5 40 13 0,33 11,5

2. Двумерная задача

Постановка задачи

Задача в двумерной постановке также является модельной. Ее можно рассматривать как следующий шаг на пути к рассмотрению задачи в трехмерной постановке.

Пусть G = {0 < X < lx ;0 < z < lz } - прямоугольник со сторонами lx, lz.

Г = {x = 0;

0 < Z < к} и Iх = к ;0 < Z < к }U

{0 < x < lx; z = 0}

{0 < x < lx;z = lz} -

его граница. Тогда G = G + Г (рисунок 3).

Прямоугольник G представляет собой проекцию параллелепипеда ABCDD'C'B'A' (рисунок 1) на плоскость Oxz.

Также введем следующие области:

G1 ={0 < x < lx;0 < z < lz} U

{ <x<lx;0<z<lz}U

{С <x < ll;0 < z <l1}U {lx < x < lj;lz < z <lz},

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

153

ФИЗИКА

G2 = {] < х < l] < z < l] },

r, = l < x < ii\z=i;},

Г] =Ц; < X < IX; z = l },

г = {x=I, i;< z < i},

Г4={x=ix i;< z < c}.

Тогда G = G, + G] ,

Г* =r, + r2 + r3 + Г4, G = G* + r*.

Рассматривается двумерное уравнение теплопроводности в плоской области, составленной из двух областей (G,, G2) с различными теплофизическими характеристиками [4, с. 17 - 21]:

pC — = L,u + L2u , х,z е G*,t e(0,T]; dt

L,u = Л

д 2u ~dx2

; L2u = Л

д 2u dz2

C =

C,; (x, z) е G, C2; (x,z) е G]

На границах раздела фаз ставятся условия идеального теплового контакта:

- du du *

Л, — = Л2— при х, z <=! ,

дп дп

lim u = lim u.

x,z^r -0 x,z^r +0

Начальные и граничные условия задаются в виде:

u

It=0

= u( x, z, 0) = To:

u| r = T0 +

,--

,

(! + Ш51 )5

T,.

(2.1)

(2.2)

Подобный выбор граничного условия (2.2) позволяет удовлетворить начальному условию (2Л).

Ставится задача определить поле температур в области G и, в частности, рассчитать время T выравнивания температуры в области G.

и

▲

x ►

Метод решения

Зададим Nx, Nz - число разбиений по пространственным координатам x, y. В области G построим равномерную сетку Оо h с шагами

^ = I , h. = ^

О =

N/ ‘ Nz

f(x>z j)> x = ihx, у, = A,

[0 < i < Nx, 0 < j < Nz

Дифференциальные операторы L,u, L2u аппроксимируем следующими разностными операторами [3, с. 394]:

Г.. (. . \ \

0 Т Т l

, lx lx lx

Рис. 3. Двумерная постановка задачи

Здесь X - коэффициент теплопроводности, р -плотность, C - удельная теплоемкость. Параметры X, р, С терпят разрывы на границах раздела фаз (области

r, r2 , r3 , r4):

A'"i,j = и

Л 2u-j = ¥

XiД/2, j (ui-,, j ui, j )

Xi+l/2, j (ui, j - "i+1,^)_

rx,j-М2 (-, - ui.j К

-%U j+V2 (ui, j - ui, j+, ^

Здесь Xi j - разностный коэффициент теплопроводности. Так как коэффициент теплопроводности Л непрерывен в G всюду, кроме линий разрыва Г , то

Л=

р=

Л; (x,z)eg,

Л(x,z)eG2

p; (x,z)eg,

P2; (x,z) e G2

положим:

Xi+1|2, j Л+,2,

Л, j + Л+и

2

XM/2, j , Xi j-1|2, Xi, j+,2 определяются аналогично. Под узловыми значениями разрывных величин

l

z

2

l

z

,54 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

ФИЗИКА

надо понимать соответствующие односторонние пределы [3, с. 382].

Производную по времени аппроксимируем следующим разностным выражением [7, с. 18]:

dt т

Пусть т - шаг по времени, тогда t^ — кт - точки разбиения по времени. Будем считать, что:

uki,j ~ u(X, zj, tk).

Для решения разностной задачи воспользуемся схемой переменных направлений, состоящей из двух шагов [3, с. 391 - 394]. Данная схема имеет второй порядок аппроксимации по временной и пространственным координатам 0(т2 + к2).

PjC

к+ 0.5 к

u, j - Ч j

i, J i, J

0.5т

— Л,м

к j 0.5

1“ i, j

Л 2u

2“ i, j’

1 < i < Nx -1;

(2.3)

ukj1 - uk j 0.5

p,. c,. , ———ли j1 j л 2ukj0.5

*, j ,, j

0.5т

1 < j < Ny -1.

1 *, J

2 *, J

(2.4)

Здесь p., Cij - значения плотности и теплоемкости в соответствующих узлах сетки соответственно. На линиях разрыва их значения определяются по формулам:

PJ ='

Р

i j1, J

р,

i-1, J

или pt . —■

p,

i, J j1

Р,

*, J-1

2

2

в зависимости от того по какой переменной выполняется прогонка. С. на линиях разрыва определяется аналогично.

Первая прогонка явная по z и неявная по x. Вторая прогонка явная по x и неявная по z.

Сетка строится таким образом, чтобы точки раз-

рыва приходились на узлы сетки 6^ к . Обозначим:

1

(

Ч — T j

Л

1—

V

(1 j 105 кт)5

Граничные условия получаются из (2.2):

и,

,к j0.5 0, J

2 (w )-

-т- Л 2 ( +,/), 1 < J < Nz -1

kj0.5 1 / kj1 . к\

UNX,j — -( jv )-

-4 Л 2 (j1 j^k), 1 < j < Nz -1

{ ukj —Ukjl, 0 < i < Nx lUN —Pkj\ 0 < i < Nx.

(2.5)

(2.6)

Граничное условие (2.5) записывается в таком виде, чтобы компенсировать погрешность на промежуточном временном слое и обеспечить второй порядок

точности по временной координате 0(т2) [3,

с. 393].

Для перехода на полуцелый временной слой к kj0.5

(U —— U ) решается последовательность Nz - 1

систем, состоящих из Nx - 1 линейных алгебраических уравнений (2.3), (2.5) с трехдиагональными матрицами. Для перехода на новый временной слой ( kj0.5 kj1

u u ) решается последовательность Nx - 1

систем, состоящих из Nz - 1 линейных алгебраических уравнений (2.4), (2.6) с трехдиагональными матрицами. СЛАУ решается методом прогонки [8, с. 35 - 38]. Численное интегрирование ведется, пока не будет выполнено условие остановки итерационного процесса.

Численные результаты

Введем линейные интенсивности армирования в направлении осей Ox, Oz (рисунок 3):

wx

к -12

xx

w„

i

X

к -1

2

z

l

z

Примем допущение, что wx = wz.

Можно показать, что с учетом принятого выше допущения зависимость линейных интенсивностей армирования wx, wz от массовой доли армирующего волокна Vm выражается в следующем виде:

Wx — Wz = ------VV(P--------Т- (2.7)

\Ра j Vm (c Ра )

Здесь рс, ра - плотность связующего и армирующего соответственно.

Ниже в таблицах 6 - 8 представлены численные данные, позволяющие судить о характерном времени выравнивания температуры в области G в зависимо -сти от использующегося связующего и массовой доли стекловолокна в композите. Здесь области G1 соответствует связующее, а области G2 - армирующее стекловолокно.

Также введем следующие обозначения:

d — ll -12 d — l1 - i2

x x x ’ ы z lz lz ■

Расчёты проводились при следующих начальных и граничных условиях: Tx — 100С°, T0 — 0С °,

S —10 7.

В качестве критерия остановки процесса принималось условие:

dT —

max(ulj - u |г)

1 < i < N x -1 1 < j < N z -1

< S

Число разбиений Nx — Nz — 100, шаг по време-

ни т — 10 6 с. Таким образом т ~ h , что согласуется с [3; 8].

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1 155

ФИЗИКА

Таблица 6

Связующее полиимид

Vm lx, 10-6 м. lz, 10-6 м. dx, 10-6 м. dz, 10-6 м. wx wz T, 10-3 с.

0,15 43 43 13 13 0,3 0,3 27,237

0,3 30 30 13 13 0,44 0,44 13,045

0,4 25 25 13 13 0,52 0,52 8,670

0,5 22 22 13 13 0,6 0,6 6,219

Таблица 7

Связующее полиэфирсульфон

Vm lx, 10-6 м. lz, 10-6 м. dx, 10-6 м. dz, 10-6 м. wx wz T, 10-3 с.

0,15 44 44 13 13 0,28 0,28 18,674

0,3 30 30 13 13 0,43 0,43 8,713

0,4 25 25 13 13 0,51 0,51 5,875

0,5 22 22 13 13 0,59 0,59 4,289

Таблица 8

Связующее полисульфон

Vm lx, 10-6 м. lz, 10-6 м. dx, 10-6 м. dz, 10-6 м. wx wz T, 10-3 с.

0,15 46 46 13 13 0,28 0,28 13,294

0,3 31 31 13 13 0,42 0,42 6,234

0,4 26 26 13 13 0,49 0,49 4,353

0,5 23 23 13 13 0,57 0,57 3,292

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

3. Трехмерная задача

Постановка задачи

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

Пусть G = |0 < х < lx;0 < y < ly;0 < z < lz }

параллелепипед со сторонами lx, ly, lz.

Г = {х = 0;0 < y <ly ;0 <z <lz }U

{х={ ;0 < y < ly ;0 <z < lz }U

{0 < x < lx; y = 0;0 < z < lz }U {0 < x <lx;y = ly;0 < z <lz}U {0 < x < lx ;0 < y < ly; z = 0}

{0 < x < lx ;0 < y < ly;z = lz}

- его граница.

Тогда G = G + Г (рисунок 4).

Также введем следующие области:

G1 ={0 < x < /';0 < y < ly;0 < z < lz} U

{ <x<{;0<y <ly;0<z<h}U { <x<{;0<y <U;0<z<lz}U { <x <{;{ <y <ly;0< z <lz}

G2 ={/J <x<Z2;% <y </2;0<z<lz},

Г, = {1 < x < l2; y = lj ;0 < z < lz}, Г2 = {lj < x < l2; y = 12;0 < z < l:},

Гз ={x=l\j'> < y < ll;0 <z < lz},

Г4 = = lx;l'y < y < £° < z < lz }.

156 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

ФИЗИКА

Рис. 4. Трехмерная постановка задачи

Тогда G — G + G2 ,

Г* —Г1 +Г2 + Г3 + Г4, G — G* + Г*.

и\Г — Т0 +

1 --

1

5 >\5

(1 +1051)

Ti-

(3.2)

Подобный выбор граничного условия (3.2) позво-Рассматривается трехмерное уршияле теплопро- ляет удовлетворить начальному условию (3.1).

водности в области, составленной из двух областей

Ставится задача определить поле температур в

GG с различными теплоф™ескими характери- области G и, в частности, рассчитать время T выравнивания температуры в области G.

стиками.

X.

дг^ — + ju + ти,

51 2

, у, z е G*, 1 е( 0, T ];

.д 2и Т — Лд 2и д 2и

L^U — Л ; Т^и — у. _ 0 ; T3U — Л

дх

22

ду2

5z

2 •

У —

Р —

C —

У; (х,у,z) е G1 .2;(хy,z) е G2

а; (х,y,z)еG1;

р2; (х,у,z)еG2’ f Q; (х, y, z) е 61 IQ; (х, У, z) е G2

На границах раздела фаз ставятся условия идеального теплового контакта:

ди ди *

Л — — Л2~ при х У, z еГ ,

дп дп

lim* и — lim* и;

х,у,2^Г -0 х,у, 2^Г +0

Граничные и начальные условия задаются в виде:

и

\1—0

— и( х, у,0) — T0,

(3.1)

Метод решения

Зададим Nx. Ny, Nz - число разбиений по пространственным координатам х, у, z. В области G построим равномерную сетку COh h h с шагами

Здесь X - коэффициент теплопроводности, р -плотность, C - удельная теплоемкость. Параметры X, р, C терпят разрывы на границах раздела фаз (области Г1, Г2, Г3, Г4):

h —

L

h

l

h_ —

N^ N^ Nz

(х,у3, z.), х — k,

®ККуК ^

>.

y} — h, z.— shz,

0 < i < Nx, 0 < j < Ny

0 < s < Nz

Дифференциальные операторы L1u, L2u, L3u аппроксимируем следующими разностными операторами [3, с. 397 - 398]:

f„L \ \

1

Лл- •— к-

Л Aj,, — K пу

Л —h

Xi-V2, j,. (('-1,(,. и', j,.)

vXi+1/2, j,. (i,(,. иi+1, j,. ) )

Xi, j-1/2,. (иi, j-1,. иi, j,. )

0

Xi, j+1/2,. ((i, (,. иi, j+1,. J J

Xi, j,.-1/2 (ui, j,.-1 - иi, ^ ^

Xi, j,.+1/2 (i,(,,

- и.

у м,j,.+1/2 ^Ч, ],. ■ i, j,.+1 ) у

Здесь /у,. - разностный коэффициент теплопроводности. Так как коэффициент теплопроводности X

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

l

157

ФИЗИКА

непрерывен в G всюду, кроме поверхностей разрыва фаз композита Г , то положим:

, _ л _ Л, j ,s + Л

Xi+l/2, j,s А+1/2, j ,s 2

.'+1, j ,s

Xi-1/2,j,s , Xi, j-1/2,s, Xi, j+1/2,s , Xi, j,s-1/2,

X. j s+1/2 определяются аналогично. Под узловыми

значениями разрывных величин надо понимать соответствующие односторонние пределы [3, с. 382].

Производную по времени аппроксимируем следующим разностным оператором [7, с. 18]:

ды ^2, — ыкj,s

dt т

Пусть т - шаг по времени, тогда t^ _ кт - точки разбиения по времени. Будем считать, что

ы.

!ы( , у, , ^s , tk ).

Для решения разностной задачи воспользуемся локально-одномерной схемой [3, с. 394 - 398]. Данная схема состоит из трех шагов.

ы к+1/3 — ы к

P C -j_______-j _ Л ык+1/3

r'i,j,y-'i,j,s yv1Mi,j,s ’

1 < i < Nx —1;

(3.3)

0i, j,sCi

ык+2/3 - ык+V3

i, J, s

i,J,s i,j,s

i, j ,s _ л ык+2/3

i'-2Mi,j,s ’

1 < j < Ny —1;

(3.4)

к+1 к+2/3

u — ы 1 , ,

0 с -2_________ Л ык+1

r'i, j, y-'i, j, s yV3M i, j ,s’

т

1 < s < N„ —1.

Здесь о , C -

л^1,j.s ’ i,j,s

(3.5)

значения плотности и теп-

лоемкости в соответствующих узлах сетки соответственно. На поверхностях разрыва их значения определяются по формулам:

°+1, J,s + Pi—1, J ,s

о _ —

^i ,J ,s

Pi,j ,s

Pi

i, J+1,s

2

+ pi,, —

i,j—1,s

Pi ,j ,s

Pi, J ,s+1 + Pi ,j ,s+1

2

в зависимости от того по какой переменной выполнятся прогонка. Ci ,■ s на поверхностях разрыва

i, j ,s

определяется аналогично.

Сетка строится таким образом, чтобы точки разрыва приходились на узлы сетки (Ohh h .

Обозначим:

/ик _т0 + 11--------5—- It;.

0 ^ (1 +105кт)5 J 1

Граничные условия получаются из (3.2):

-tt _Мк+;,0 < i < Nx,0 < s < Nz -kh+2,s _S+',0 < i < Nx,0 < s < Nz Ыok,+;,s _s+;,1 < j < Ny —1,0 < s < Nz ' -N+J _ик+\1 < j < Ny —1,0 < s < Nz -к+1 _цк+\1 < i < Nx —1,1 < j < Ny — 1 -k, +,N _Мк+',1 < i < Nx —1,1 < j < Ny — 1.

Для перехода на промежуточный временной слой

к к+3

ы ^ ы 3 решается последовательность (Ny -

1)x(Nz - 1) систем, состоящих из Nx - 1 линейных алгебраических уравнений (3.3), (3.6) с трехдиагональными матрицами.

Для перехода на промежуточный временной слой

к+-

ы 3 ^ ы

к + 2 3

решается последовательность

(Nx - 1)x(Nz - 1) систем, состоящих из Ny - 1 линейных алгебраических уравнений (3.4), (3.6) с трехдиагональными матрицами.

Для перехода на новый временной слой

к+— ,

ы 3 ^ ык+1

решается последовательность (Nx - 1)x(Ny - 1) систем, состоящих из Nz - 1 линейных алгебраических уравнений (3.5), (3.6) с трехдиагональными матрицами. СЛАУ решается методом прогонки [8, с. 35 - 38].

Численные результаты

Введем линейные интенсивности армирования в направлении осей Ox, Oz (рисунок 4):

11 — 12

wx _

И — z

i

wz _

l

Примем допущение, что wx = wz.

Можно показать, что зависимость интенсивности армирования w от массовой доли армирующего волокна Vm выражается в формулой (2.7).

Ниже в таблицах 9 - 11 представлены численные данные, позволяющие судить о характерном времени выравнивания температуры в области G, в зависимости от использующегося связующего и массовой доли стекловолокна в композите. Здесь области G; соответствует связующее, а области G2 - армирующее стекловолокно.

В качестве критерия остановки процесса принималось условие:

dT

maxO^s

1<i< Nx —1

1< J< Ny —1 1 <s< Nz —1

Щ г )

< s.

158 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

ФИЗИКА

Расчёты проводились при следующих начальных и граничных условиях: Т =

£ = 10"7.

Число разбиений Nx = Ny = Nz = 50, шаг по времени т = 10-6 с.

Таким образом т ~ h, что согласуется с [3; 8].

Vm - объемная доля стекловолокна в композите.

100С°, Т0 = 0C °,

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

Таблица 9

Связующее полиимид

Vm lx, 10-6 м. ly, 10-6 м. lz, 10-6 м. dx = £ -£ У 10-6 м. dy=ч - Ч У 10-6 м. d w = — l dy w=ту T, 10-3 с.

0,15 43 43 43 13 13 0,3 0,3 11,05

0,3 30 30 30 13 13 0,44 0,44 4,02

0,4 25 25 25 13 13 0,52 0,52 2,451

0,5 22 22 22 13 13 0,6 0,6 1,69

Таблица 10

Связующее полиэфирсульфон

Vm lx, 10-6 м. ly, 10-6 м. lz, 10-6 м. dx = l1 - lX У 10-6 м. dy = l'y - 1У У 10-6 м. d w = — l dy у = ~y T, 10-3 с.

0,15 44 44 44 13 13 0,28 0,28 8,598

0,3 30 30 30 13 13 0,43 0,43 3,23

0,4 25 25 25 13 13 0,51 0,51 2,04

0,5 22 22 22 13 13 0,59 0,59 1,446

Таблица 11

Связующее полисульфон

Vm lx, 10-6 м. ly, 10-6 м. lz, 10-6 м. dx = 11 - lX, 10-6 м. dy = l'y - ly У 10-6 м. d w = — l dy у=^ T, 10-3 с.

0,15 46 46 46 13 13 0,28 0,28 6,782

0,3 31 31 31 13 13 0,42 0,42 2,804

0,4 26 26 26 13 13 0,49 0,49 1,857

0,5 23 23 23 13 13 0,57 0,57 1,369

Достоверность результатов В этом разделе предпринимается попытка оценить достоверность полученных результатов. Последовательность сеток Для одномерной задачи была построена последовательность сеток. Было исследовано поведение решения задачи в зависимости от числа разбиений Nx.

При этом остальные параметры задачи были зафиксированы.

Расчёты проводились при следующих начальных

и граничных условиях: Т = 100C°, Т0 = 0C °,

£ = 10“7. В качестве критерия остановки процесса принималось условие:

dT =

max(uf+1 - ui)

ь

1<i<NX-1

< £.

Шаг по времени Г = 10 6 с.

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1 159

ФИЗИКА

Рис. 5. Профиль графика температуры в случае одномерной задачи

Интенсивность армирования была принята w = 0,4 .lx = 30 х 10'6 . Это соответствует мас-

совой доли армирующего волокна Vm = 0,3 для композиции, составленной из стекловолокна и полиимида.

На рисунке 5 представлен профиль графика температуры на момент остановки расчета. Линии 1 соответствует решение задачи при N = 10, линии 2 при N = 20, линии 3 при N = 40, линии 4 при N = 80. При увеличении числа узлов сетки наблюдается сходимость решения разностной задачи.

Аналогичная последовательность сеток была построена для трехмерной задачи. На рисунке 6 представлен профиль графика температуры в момент остановки процесса вдоль трех прямых, проходящих через центр параллелепипеда:

У

x

К z = К 2 , z 2

hx

h.

(а);

(б);

x =

h

x

2

к

2

(в).

Линиям 1 соответствует решение задачи при N = 10, линиям 2 при N = 20, линиям 3 при N = 40, линиям 4 при N = 60. При увеличении числа узлов сетки наблюдается сходимость решения разностной задачи.

Также для небольших размерностей сеток (N = 40) результат, полученный с применением экономичной схемы (2.3) - (2.6), сравнивался с результатом, полученным с применением чисто неявной схемы. Результаты представлены в таблице 12.

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

160 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

ФИЗИКА

б)

в)

1 1 1

\

_ . .. . .

// .

\,4.

V s / fy

\v\ у // .

\N\ : / //

■ / /Y

/ /V -

\\ !\ ! f //

\\ V ■ %'s ч% X у .

У // X /'/ -

■ \ i ; X /У . SJ f &

_ ...\V : e / Ar s' у Ж • • -

■ " ''у ■

УЖ

Xs ''Ч., 2 -

4 1 1 1 1 1

15

У, 106 I

Рис. 6. Распределение температуры вдоль координатных осей

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

161

ФИЗИКА

Сравнение чисто неявной и экономичной разностных схем

Таблица 12

Vm lx, 10-6 м. ly, 10-6 м. dx У 10-6 м. dy y, 106 м . w wy ТЭКОНОМ, 10-3 с. ТНЕЯВН, 10-3 с. в, %

0,15 43 43 13 13 0,3 0,3 18,454 17,756 3,93

0,3 30 30 13 13 0,44 0,44 8,778 8,473 3,6

0,5 22 22 13 13 0,6 0,6 4,189 4,05 3,43

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

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

Зависимость времени установления стационарного режима задачи от начальных и граничных условий

Было проведено исследование зависимости времени установления стационарного режима задачи в зависимости от начальных данных. В частности, из-

менялось значение AT = T — T0.

Расчёты проводились с использованием разностной схемы (3.3) - (3.6) при следующих параметрах:

N = 20, е = 10-7, т = 10-6 с.

Связующее - полиимид. Геометрические параметры пластинки следующие:

1. 4 = ly = 4 = 43 х 10-6 м., wx = wz = 0,3.

2. lx = ly = 4 = 30 х 10-6 м., wx = wz = 0,44.

3. lx = ly = lz = 22 х 10-6 м., wx = wz = 0,6.

Таблица 13

Зависимость времени установления стационарного режима задачи от AT, 10-3 с

ДТ 0,1 1 10 50 100

1 6,853 7,956 9,06 9,831 10,163

2 2,779 3,228 3,676 3,99 4,125

3 1,097 1,275 1,453 1,577 1,631

4

2

0

0,1

1

10

50

100

Рис. 7. Время установления стационарного режима в зависимости от AT, 10-3 с

Из таблицы 13 и рисунка 7 видно, что даже значительное изменение ДТ не оказывает существенного влияния на время установления стационарного режима задачи.

Выводы

Задача рассматривалась в трех постановках: одномерной, двумерной и трехмерной. Полученные результаты во всех трех вариантах постановки задачи имеют одинаковый порядок (порядка 10-3 с) и много

162 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

ФИЗИКА |

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

лов. Поэтому можно сделать вывод о правомерности

Литература

1. Андреев А. Н. Дифференциальные уравнения связанной задачи термоупругого деформирования слоистой композитной оболочки // Известия Алтайского государственного университета. (Серия: Математика и механика. Управление, вычислительная техника и информатика). Физика. Барнаул: Издательство Алтайского государственного университета. 2012. № 1/1(73). С. 11 - 13.

2. Андреев А. Н. Упругость и термоупругость слоистых композитных оболочек. Palmarium Academic Publishing, 2013. 93 с.

3. Калиткин Н. Н. Численные методы. М.: Наука, 1978. 512 с.

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

5. Композиционные материалы: Справочник / В. В. Васильев, В. Д. Протасов, В. В. Болотин, Н. А. Алфу-тов, А. И. Бейль, В. А. Бунаков, И. А. Дымков, А. Ф. Ермоленко, И. Г. Жигун, П. А. Зиновьев, Т. Я. Кинцис,

В. В. Клейменов, А. А. Круклиньш, А. А. Кульков, Ф. В. Мануйлов, Б. Г. Попов, Г. Г. Портнов, О. С. Сироткин, А. М. Скудра, И. А. Соловьев, Ю. М. Тарнопольский, Ю. С. Царахов; под общ. ред. В. В. Васильева, Ю. М. Тарнопольского. М.: Машиностроение, 1990. 512 с.

6. Немировский Ю. В., Янковский А. П. Рациональное проектирование армированных конструкций. Новосибирск: Наука, 2002. 488 с.

7. Самарский А. А. Введение в теорию разностных схем. М.: Наука, 1971. 552 с.

8. Самарский А. А. Теория разностных схем. М.: Наука, 1989. 616 с.

9. Справочник по композиционным материалам: в 2 кн. Кн. 1 / под ред. Дж. Любина. М.: Машиностроение, 1988. 448 с.

10. Справочник по композиционным материалам: в 2 кн. Кн. 2 / под ред. Дж. Любина. М.: Машиностроение, 1988. 448 с.

11. Material: Polyimide. Режим доступа: http://www.mit.edu/~6.777/matprops/polyimide.htm (дата обращения:

16.05.2014).

12. Polysulfon (PSU). Режим доступа: http://www.makeitfrom.com/material-data/?for=Polysulfone-PSU (дата обращения: 16.05.2014).

13. Polyethersulfone (PES, Radel A). Режим доступа: http://www.makeitfrom.com/material-data/7foFPolyether-sulfone-PES -Radel-A (дата обращения: 16.05.2014).

Информация об авторе:

Станиславович Алексей Вячеславович - программист «ПКП РИФ», соискатель кафедры алгебры и геометрии КемГУ, [email protected].

Alexey V. Stanislavovich - post-graduate student at the Department of Algebra and Geometry, Kemerovo State University.

Статья поступила в редколлегию 18.09.2014 г.

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1 163