Численное моделирование процессов тепломассопереноса и кинетики напряжений в термодеструктирующих композитных оболочках Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 17, № 2, 2012

Численное моделирование процессов тепломассопереноса и кинетики напряжений в термодеструктирующих композитных оболочках*

Ю.И. Димитриенко1, В. В. Минин1, Е. К. Сыздыков2 1 Московский государственный технический университет им. Н.Э. Баумана, Россия 2Государственное машиностроительное конструкторское бюро "Радуга" им. А.Я. Березняка, Москва, Россия e-mail: [email protected]

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

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

Оболочки и пластины из полимерных композитных материалов достаточно активно применяют при создании теплонагруженных элементов конструкций летательных аппаратов. Особый интерес имеет проблема исследования термомеханических эффектов взаимодействия излучений с композитными материалами. В отличие от металлических материалов в композитах даже при сравнительно невысоком уровне нагрева (до 1000 °С), когда процессы испарения еще не реализуются, возникает явление термодеструкции полимерной матрицы, сопровождающееся внутренним газообразованием и образованием пиролитической фазы. Следствием внутренних физико-химических превращений матрицы являются специфические эффекты термомеханического поведения композитов: изменение знака термического расширения на отрицательный (химическая усадка), резкое уменьшение значений модулей упругости в зоне нагрева, возникновение локальных пиков порового давления [1-5]. В [1, 2] была разработана теория четырех-фазного термодеструктирующего композита, описывающая основные закономерности и особенности поведения композитов при нагреве до высоких температур. Целью настоящей работы является применение этой теории для моделирования термомеханического

* Работа выполнена при финансовой поддержке РФФИ (грант № 09-08-00323-а).

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

1. Уравнения теории термодеструктирующих композитных оболочек

Система уравнений теории термодеструктирующих композитных оболочек в ортонор-мированных координатах дх, д2, д3, связанных со срединной поверхностью оболочки (ось Од3 направлена по нормали к срединной поверхности), имеет следующий вид [1]:

00 _ __д_ / А дв_\ __д_ / Ах \ А л \

рст = а1А221Ао^) + ААд&V 2Ад^у + дц~Л 3+

, Кг дд dpg К2 дд dpg дд dpg г Q J

+cA A dq; dqg + а d^ dqg + Кз dq~3 dqg1 - AeJ,

dPg _ 1 д ( A2K1 dPg \ + __д- ( AlK2 dPg ) + — (KidP^\ +rJ

dt A1A2 дд;\ A; дд;) A1A2 дд^ A2 дд2/ дд3\ дд3 дрь j j j ( Ea

Pb-

m -J J = J0<fb exp Ed), Pg = Rgд Pg/^g - (1)

уравнения внутреннего тепломассопереноса в оболочке;

д д дА дА дР

-(АрТаа) + — (АаТав) - Твв + ^Тав + АаА^каЯа - Ар= 0,

о V -"-в -*- аа) | о V а -*- ав / о вв 1 о -*- ав \ ^ • ^а^о а г\

дда ддв дда ддв дда

д д дA дA д^Д

(AßМаа) + Tj ^аМав) - -АMßß + —^Мав - АаАвQа - Aß = 0,

-A;A2(k;Tii + k2T22) + + - PeAlA2 - (к; + k2)A;A2'^gPg = 0,

дда ддр дда ддр дд,

дА^Ях + дАЯ ддх дд

а, в = 1, 2; а = в -

уравнения равновесия оболочки;

1 диа 1 дАатт , Т1Г 1 дШ

+ --о-ив + каШ, 2ва3 = -¡-Т.--+ 7а - каиа,

А А а дда

а

Aа дда Ai A2 ддв ' Aа дд,

(U; + ^ U2)

V дд2 дд; /

2 = ± -UI + дЦ2 - /-AIЩ + -AIU2 12 A2 дд2 A; дд; A;A2 V дд2 ; ° 2

1 д^а , 1 дАа

Каа = --+ л л ~Ö-1в, 2Ка3 = -ка^а,

Аа дда А А ддв

2 = —du + 1 (-А; + -Ai \

К;2 = А2 дд2 + А; дд; А;А^ дд2 Ъ + дд; Ъ) кинематические соотношения;

2

Таа (Савевв + Квв) - pgа - Та,

2

Maa = (N*/3 eee + Dae К/З/З) - Mga — Mt

a i

¡3=1

(4)

определяющие соотношения. Здесь Taa, Tap, Maa, Map — усилия и моменты в оболочке; Qa — перерезывающие усилия; Pg, Mg — усилие и момент порового давления в оболочке; eaa, еаз, ei2 — деформации срединной поверхности оболочки; Kaa, ка3, к12 — искривления срединной поверхности; Ua, ja, W — перемещения, углы искривления и прогиб срединной поверхности; Aa, ka — параметры первой квадратичной формы и главные кривизны срединной поверхности оболочки [3]; в — температура; pb, pp, pg, pf — объемные концентрации фаз композита (полимерной фазы матрицы, пиролитической фазы, газовой фазы и фазы армирующего наполнителя), причем pg = 1 — pb — pp — Pf, Pf = const (полагаем, что волокна не деструктируют), pp = (pb(0) — pb(t))(pb/рР) (1 — Г); pb, pp, Pf — плотности твердых фаз (полагаем постоянными); pg = р0pg — приведенная плотность газовой фазы; pg — истинная плотность газовой фазы; pg — давление газообразных продуктов в порах композита; cb, cp, cg, Cf — теплоемкости фаз; p = pb pb+pppp + pg + Pf Pf — плотность композита; c = (pbCb'b + ppcppp + Pg Cg + Pf cf pf )/p — теплоемкость композита; Г — коэффициент газификации; J — массовая скорость термодеструкции; Ав° — удельная теплота терморазложения; Ea — энергия активации, Аа и Ka — коэффициенты теплопроводности и газопроницаемости композита; Rg — газовая постоянная.

На внешней поверхности оболочки q3 = h/2 действует локальный нагрев излучением, и соответствующее тепловое граничное условие к системе (4) имеет вид q\ = q0+, где q\ = А3 (de/dq3) — тепловой поток, затрачиваемый на нагрев оболочки; q0+ = qC + qR — qBL — qnw — суммарный тепловой поток к поверхности; qC = (aT/cge)(cgeee — cS) — конвективный тепловой поток; aT — коэффициент теплообмена; qu = qR exp(—r/8) — лучистый тепловой поток, подводимый к оболочке за счет локального нагрева; 8 — эффективный радиус пятна нагрева; r = \Jql + q|, qRW = в4 — лучистый тепловой поток, отводимый от нагретой поверхности композита; asB — постоянная Стефана — Больцмана; — интегральный коэффициент излучения; qBL = YbiPgvg3(cgeSe — cS) — конвективный тепловой поток, отводимый от поверхности вследствие вдува продуктов через поры; Ybl — коэффициент вдува; Vg3 = —RgK^d(PgS)/dq 3 — нормальная скорость фильтрации газообразных продуктов в порах; ве и cge — температура и теплоемкость внешней газовой среды.

Мембранные, смешанные и изгибные жесткости оболочки Cap, Nap, Dap вследствие размягчения полимерной матрицы и ее термодеструкции изменяются при нагреве. Это изменение для ортотропных композитных оболочек тканевой структуры, согласно модели [1], описывается с помощью двух функций а$1 и ав2:

Ca+3 C a+3 а<)2 , a = 1, 2,

h/2

(0)

h/2

где Сав — мембранные жесткости при нормальной температуре в = в0 = 293 К. Функции а$к = &ек(в,рь) зависят от температуры в и концентрации полимера рь, их выражение приведено в [1].

о о

Усилия Та и моменты Ма тепловых напряжений зависят от тепловой деформации

о

£а композита:

3 3

Та = ^^ Сав £в ) , Ма = ^^ С ав £]з ^,

в=1 в=1

Н/2 Н/2

Щ> = ]ап £в ,3 ^ £3" = /а„2 £з ,3 } = 0, ,, а = X (6)

—Н/2 —Н/2

г

°е1 = (а] р]Б1 + аьрь^-у)(в - в0) + (в(Ь) - в(т))рр ¿т - , 7 = 1, 2, 3. (7)

о

Здесь а], аь, ар — коэффициенты теплового расширения волокна, полимера и коксового остатка термодеструкции полимера соответственно; вр — коэффициент усадки; Б1, — коэффициенты, зависящие от расположения волокон в композите [1].

Усилия Рда и моменты Мда межфазного взаимодействия, а также усилия Рд и моменты Мд порового давления в оболочке вычисляются следующим образом:

Н/2 Н/2

Рда = J 'РдIа ¿,3, Мда = ^ 'Рд¡а,3 ¿,3, —Н/2 —Н/2

Н/2 Н/2

Рд = J Рд Рд ¿,3, Мд = J Рд'Рд ,3 ¿,3, (8)

—Н/2 —Н/2

где рд — пористость композита: рд = 1 - рь - рр - р], 1а — коэффициент межфазного взаимодействия [1].

Коэффициенты теплопроводности Аа = Аа(в, рь, р], рд) являются функциями температуры и содержания фаз композита и вычисляются по формулам, приведенным в работе [1]. Коэффициенты теплоемкости фаз Сь, ср, сд, с] полагаются зависящими от температуры [1]. Коэффициенты газопроницаемости Ка вычисляются по формуле

Ка = Ка ехр^рУ3), где Ка и в — константы.

2. Вариационная постановка задачи

Сформулируем вариационную постановку задачи механики (2)-(4), основанную на вариационном принципе типа принципа Хеллинера — Рейсснера [6]. Вариационное уравнение типа уравнения Хеллингера — Рейсснера для задачи (2)-(4) имеет вид 5.]^(и, е) = 0, где 3£(п,е) — следующий функционал:

Ыще) = I {е}т[Е][Ь] {и} ¿Е - ^ {е}т[Е] {е} ¿Е + А£, (9)

£ £

AS — работа суммарных внешних сил:

Ae

J {S} {u} dl - J ({T} + {F}) {u} dE.

dS s

Здесь

{u} = (Ui, U2, W, Yi, Y2)T, {e} = (en, 022, ei2, кп, K22, Ki2, e13, е2з)Т, {T} = {Tii, T22, T12, Mii, M22, M12, Qi, Q2}T,

{S} = {Tii,T22, Mii,M22, Q}T —

:iq)

координатные столбцы перемещений, деформаций, усилий и заданных усилий на контуре дЕ срединной поверхности; ре — перепад давлений на оболочке,

{Р} = {дРд/ддг, дРд/дд2, (кг + к2)Рд, дЫд/ддъ дЫд/дд2}т,

{T} = {T u, T22, 0, M a, M22, 0, 0, 0}

T

Выше обозначены: [E] — обобщенная матрица упругости размером 8 х 8, связывающая столбцы усилий и деформаций {T} = [E] {e}; эта матрица образована из матриц [C], [N], [D], [K]:

'[C] [N] 0\ (Си

[E] = I [N] [D] 0 ! , [С] = I Ci2

0 0 [K]

Nii Ni2 0

[N] = I Ni2 N22 0 0 0 N66j

[D]

Dii Di2 Di2 D22 00

0

0 0

Dm.

Ci2 0 C22 0 0 C66,

[K]

C44 0

0

C55

:i2)

Матрица частных производных [L] имеет вид

[L]

(8x5)

lii l12 ki 0 0

l2i l22 k2 0 0

(I22 - l2i)/2 (lii - li2)/2 0 0 0

0 0 0 lii l12

0 0 0 l2i l11

0 0 0 (l22 - l2i)/2 (l11 - l12

-ki/2 0 lii/2 1/2 0

0 -h/2 l22/2 0 1/2

:i3)

где дифференциальные операторы

1 1

Aa dga

i2

1 dAi AiA2 dg2

2i

1 ЗЛ2 ЛЛ dgi

Полагая вариации ${и} и 5{е} независимыми, из (9) получаем две группы вариационных уравнений

{e}T[D] [L] ¿{u} dE + = 0, I 5{e}T[D] ([L] {u} - {e}) dE = 0.

:i4)

3. Применение метода конечного элемента для моделирования термонапряжений в композитных оболочках

Решение вариационных уравнений (14) будем искать методом конечного элемента. В качестве конечного элемента (КЭ) выберем треугольный шестиузловой КЭ с независимой аппроксимацией перемещений и деформаций

{и} = [Ф] {V}, {е} = [ш] {Ь},

:15)

где — столбец перемещений в узлах размером 30 х 1; {Ь} — столбец деформаций в узлах размером 24 х 1; [Ф], [ш] — матрицы функций формы с размерами 5 х 30 и 8 х 24 соответственно. Подставляя (15) в систему (14), получим уравнения для нахождения {V} и {Ь}:

здесь

[С] [Я]-1[С] {V} + {/} = 0, {Ь} = [Я]-1[С]>}

[С] = /[ш]т[£][В] ¿Е, [Я] = /[ш]т[£][ш] ¿Е,

'16)

[В] = [Ь] [Ф], {/} = - у ({Т} + {^}) ¿Е - ^ {5} ^ (17)

Е ЙЕ

Матрицу функций формы [Ф] представим в блочном матричном виде:

(5 Ф3П) = [Ф(1), Ф(2), . . . , Ф(4

(5x30)

где Ф(г) = фг I , I — единичная матрица размером 5 х 5, фг (г = 1, ..., 6) — квад-

(5x5)

(5x5)

ратичные функции формы, для представления которых удобно воспользоваться естественными безразмерными (барицентрическими) координатами [7]:

ф1 = Ь^Ь - 1), ф2 = ¿2(2^2 - 1), фз = Ьз(2Ьз - 1),

ф4 = 4Ь^2, ф5 = 4Ь2^з, фб = 4ЬзЬь

Нумерация функций формы соответствует нумерации узлов элемента. Барицентрические координаты Ьг определяются через координаты 91, 92 оболочки следующим образом:

Ьг

1

(а(г) + Ь(г)^1 + 0(г)?2) , г = 1, 2, 3

где

а(1) = 91(2)92(з) - ?1(з)?2(2), Ь(1) = 92(2) - 92(з), С(1) = 91(з) - 91(2), 25 = Ь(1)С(2) - Ь(2)С(1),

91(г), 92(г) (г =1, 2, 3) — координаты узлов треугольного КЭ. Остальные коэффициенты ау), 6(/), 0(1) находятся с помощью круговой перестановки индексов, заключенных в круглые скобки.

Матрица функций формы [ш] имеет блочную диагональную структуру

ш

(8x24)

Ь1 0 0 Ь2 0 0 Ьз 0 0 0 ¿1 0 0 ¿2 0 0 Ьз 0 0 0 Ь1 0 0 Ь2 0 0 Ьз

Матрица [В] имеет размер 8 х 30 и с учетом блочной структуры матрицы функций

форм Ф ее можно представить в виде В = [Вт, ... , В(б)], где матричные блоки В(г)

(8хз0)

(г = 1, ... , 6) определяются соотношением В(г) = [Ь]Ф(г). Подставляя выражение (13)

(8x5)

для матрицы [Ь] в это соотношение, получим

В(г) (8x5)

7 . 7жг 0 0 0 0

0 7 7уг 0 0 0

7 7уг 7 7жг 0 0 0

0 0 0 7 7жг 0

0 0 0 0 7 7уг

0 0 0 7 7уг 7 7жг

0 0 7 7жг фг 0

0 0 7 7уг 0 фг

Здесь

:19)

= (Ь(г)/25) (4Ьг - 1), = (0(г)/25) (4Ьг - 1), г = 1, 2, 3,

2

22 7ж4 = 5 (Ь(1)Ь2 + Ь(2) Ь1) , 7ж5 = 5 (Ь(2)Ьз + Ь(з) Ь2) , 7жб = 5 (Ь(з)Ь1 + Ь(1)Ьз^

2

2

7у4 = 5 (о(1)Ь2 + 0(2)Ь1) , 7у5 = 5 (о(2)Ьз + 0(з)Ь2) , Матрица Н для отдельного выбранного КЭ имеет вид

2

5 (о(з)Ь1 + 0(1) Ьз).

Н

(24x24)

Ь1Ь1В Ь1Ь2В Ь1ЬзВ Ь^В ¿2^В Ь2ЬзВ Ь1ЬзВ Ь2ЬзВ ЬзЬзВ

¿Е.

(20)

Для решения систем линейных уравнений (16) был применен метод сопряженного градиента.

4. Метод решения уравнений тепломассопереноса в оболочке

Непосредственное решение системы уравнений тепломассопереноса (1) в трехмерной постановке для тонкостенных оболочек связано с определенными трудностями при использовании как конечно-элементных, так и конечно-разностных методов. Для того чтобы получить достаточно высокую точность решения задачи (обычно важно иметь достаточно точное распределение полей температуры и порового давления по толщине), необходимо применять очень мелкие сетки, с числом N элементов по толщине, как правило, не менее 20-30. Это приводит к необходимости генерирования сеток с большим числом трехмерных КЭ для всей оболочки — О(^(/0/Л,)2) ~ 3 • 108, где /0 — характерный размер оболочки по координатам 91, 92. Учитывая, что система (1) является нелинейной и для ее решения приходится применять итерационные методы, использование больших сеток приводит к значительному времени вычислений. Аналогичные проблемы возникают и в случае применения конечно-разностных методов.

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

разделена на две части: Ее1 (в которой задана "основная часть" плотности лучистого теплового потока дп) и Ее2 (оставшаяся часть внешней поверхности, где плотность теплового потока дп существенно меньше). Тогда решение системы (1) можно также разделить на два этапа: решение в области У\, граница которой содержит поверхность Ее1, и решение в оставшейся области У2 задачи (VIиУ2 = V). В области VI будем решать общую систему уравнений (1) в трехмерной постановке, а в области V2 — упрощенную одномерную (по координате д3) систему уравнений, решение которой зависит от координат д2 только параметрически. В строгой постановке решение этих задач связано за счет граничных условий на поверхности раздела Е12 областей V1 и V2. Однако если часть поверхности Ее1 выбрана таким образом, что на ее границе дЕе1 плотность теплового потока дп существенно уменьшается по сравнению с максимальным его значением, то можно использовать приближенную постановку, в которой граница раздела Е12 считается теплоизолированной. В результате задачи для областей V1 и V2 становятся не связанными. Размеры части Ее1, обеспечивающие требуемые условия уменьшения плотности теплового потока на границе Е12, подбираются экспериментально путем численного решения задачи для области V1. Для рассматриваемых относительно малых времен прогрева фронт прогрева от локального источника не успевает дойти до границы локальной области Е12, поэтому условие тепловой изоляции для этих моментов времени выполняется с высокой степенью точности.

Если на поверхности оболочки действует только плотность лучистого теплового потока дп, то можно ограничиться решением задачи для области V!, что далее и было осуществлено. Если кроме локального теплового потока излучением на внешней поверхности действует "медленно изменяющийся" конвективный тепловой поток дс, то необходимо решение обеих задач: для областей V1 и V2.

Будем рассматривать систему (1) для тонкой оболочки и запишем ее в безразмерном виде. Для этого введем в рассмотрение характерные значения: времени ¿о, продольных координат д1 и д2 — 1о, поперечной координаты Н (толщина оболочки), коэффициентов квадратичной формы срединной поверхности Л^, плотности ро, удельной теплоемкости со, температуры во, теплопроводности ко, газопроницаемости Ко, скорости тепломассо-переноса 1о, а также соответствующие им безразмерные величины

г = г/и, % = д1 л°/1о, сз = дз/н, Л„ = Ла/Л°а,

Рд = Рд <Рд/Ро, Рь = Рь/Ро, Рр = Рр/Ро, р/ = Р/ /ро, р = р/ро, в = в/во, % = Сг/со (г = ¡,Ъ, р, д), с = с/со,

Ха = Ха/Хо, Ка = Ка/Ко, 1 = 1/1, Ре± = ^ п , % ± = . (21)

РоЩ во Хово

Тогда систему уравнений тепломассопереноса (1) в безразмерном виде можно представить следующим образом:

РЬ ^ = —Рд 3, Рд = Рдв/>£д , (22)

дрд = (ЛК дрЛ + А (ЛК. дрЛ\ + ¥^(Кд_р

дс Л 1Л2\дС1\ Л 1 дС1) + д%2 \ Л2 дд2)) + гдсз V 3дс

__дб = Foß2 //Л1А2 Зб-

= A?A V дл V ~АГ дЛ

+

д /%A? дб

дq2 V A2 дq2

^ д /л дб . + Fo— Л3—- ] +

длз V ддз,

+ß 2 Fr Л

к? дб др K2 дб дрл

+

A? дql дql A2 дq2 дq2

+ Fr Cg K3 —

дб др

где

Fo

Лoto F Pocoh2' r

KoRg 6oto

Fg

Joto

Ft

дqз дqз

Aeo

Л2 ' ^ ро ' ^ оо^о; безразмерные параметры (критерии). Граничные и начальные условия, налагаемые на эту систему, в безразмерном виде записываются как

- FtFg J

ß=h— lo

(24)

(25)

_ _ l e дб _

q3 = 7: : k3^~ = Фж

2 ддз

l e дб - дрб = Лз = -2: -k3д- =Co-, Кз дЛ

O,

t = 0: Pg = PgO, б = 00, ^p = ^pü. (26)

Задачу (22)-(24), (26) для области Vi будем называть локальной задачей тепломассопереноса. Для ее численного решения применим конечно-разностный метод в сочетании с пошаговым методом линеаризации. Введем конечно-разностную сетку для области V1, узлы которой обозначим g^, g2j, g3k (i = 1, ... , j = 1, ... , N2, k = 1, ... , N3). Обозначим значения функций в узлах разностной сетки на m-м временном шаге 6™fc, P^^jfe и введем разностные операторы

л^1

Fr ß 2 / /Â2K1

(A1A2 )ijk Aq2U Al

i+1,jk

(ífli+ljk iflij+k )

A2K1

A?

ijk

(V^1 — ñm+1 Л

V^gijfc q i 1 jfc/

Лр2РТ+1

Fr ß2

(A iA2)ijk Aq;

AiK

'A 1K2

A2

ij+1,k

Aem+1 _ pm+n _ V^gij+^k pgijk J

- i V^gijfc ^gij A2 /ij—1,k

^ - jk) ь

Fr

Лlбm+1

ЛР3Р!Г+1 = Aq2 (Kmijfc+i G^jfc+i ^j?) K3mj G^jfc1 ^j—1)) ,

_F_ i / ZCm+l _ Cm+1) _ / Л1А2\ /ñm+^ em+1 Л

(AiA2)ijkAq2 U Al Ji+ij k ^k 6mj k ^ V Al Jijk ^ + 6-1 j kJ

Fo

Л 2A

^ = (nijqf (w k - 6mm+1) - ) (л^- ji,,

ijk

Д flm+1 _ o ( um /nm+l e m+1A / m /nm+l em+l ЛЛ

б = Aq2 Vk3,ijk+H6ijk+1 - 6ijk )- k3,ijk i6ijk - 6ijk-lJJ

fm

ß2Fr /-gК?

4Aq2 V a? y ijk

ß2 Fr / -g K2 4 m

4Aq22 A 22

m m m m

^i+ljk 6í—1,jk )(-Fgi+1,jk pi—1,jk)ч

m _ Cm )(,em _ Л1 )I

ij+1,k "ij—1,k )(-pgij+1,k -Fgij—1,k)4

2 / ijk

m

m

m

m

(сд К3 )™к (вСтк+1 вРтк-1)(ррт,]к+1 рдчк-1) Л^к,

4Ад чк\"цк+1 "г]к-1)\Гдг]к+1 гдг]к-1) ^ г*- д

¡д = Рдк, рдг+к = рдг+квтк/ * дг]к, Рдг]к = рдг]кк/ * дг]к, (27)

где Ад1, Ад2, Ад3, АЬ — шаги разностной сетки. Тогда для системы уравнений (1) можно записать следующую трехслойную разностную схему:

РьШ1 — <р&к) = —АрРд 1(ч>д£1 Щк),

{рдг+к 1 — р дг]к) = А^ЛР1р9 + ^ , (Рдг+к ^ — р дг+к ^ ) = А^ЛР2р9 + ^ ,

(-ш+3/4 _ -т+1/2) = АЬ . -т+3/4 \рдг]к Рдг]к ) = 3 ЛР3рд ,

Шк1 — Р7г3к) = А1{Лр1 + Лр2 + Ар3)рПП+3/4 + ¡7 аь, (рЩк{Щ+1/4 — С) = Ал^114, (рс)Ткт1/2 — Щ+1/4) = АЛ2%п+1//,

(рс^к т3/4 — Щ+1/2) = А Л3ед+3/4)

аь

У

(рЩк ш1 — Щк) = А1(лх + Л2 + а3)0п+3/4 + ! тАь.

Разностная аппроксимация граничных условий (26) имеет вид

Шм3 = 1 : —(Х^Л Щ+Т — (сППЖ-1) = (Ас3)(со+)?, С = {*д Рел/О X

2" г31\з\" %]М'3 г],М'з — 1) Гдг3,Мз угд^е1/») цм^

Ш1 = — 1: (к3)Тп Ш3/4 — с3/4) = т^д,

(РПт+3/4от/*ПП)г1,2 — (рт+3/4°т /*тС = о. (29)

Поверхность раздела Е12 областей V1 и V2 выберем в виде совокупности четырех координатных поверхностей р1 = р1± и р2 = р2±, которые обозначим как Е12д± и Е12,2±. На этих поверхностях заданы условия теплоизоляции и герметичности:

% = с 1-: (ртт+^кс — (рт+1/4вт/*пт)1]к = о, с1/4 — с1/4 = %

с1 = р1+: (рт+1/4°т/*т)с — (ррт+1/4°т/*7)м1-с = о, с — вд-4 к =

[2

д 1 ^д ) N-1] к V д 1 г д ) N-1-1,3 к

Р2-: (рр?+1/2ет/*т)2к — (Рт+1/2)ак = Рт+1/2 — 1/2 = о,

р2 = р2+: (рт+1/2р>т/€)^ к — (рт+1/2етЮмк = о, & — = о. т

Для решения разностных систем уравнений (28)-(30) применялся метод матричной прогонки [8] по координатным направлениям г, ], к.

После определения полей температуры в(да, Ь), содержания фаз *р(да, Ь), *д(да, Ь) и плотности газовой фазы рд (да, Ь) рассчитывались поровое давление рд = Крвд и тепловые деформации композита по формулам (6), (7), а затем усилия и моменты (8) и изменение жесткостей оболочки при нагреве по формулам (5), затем определялись перемещения и деформации {е} в КЭ.

5. Расчет усилий, моментов и напряжений в оболочке

о

Усилия и моменты в КЭ оболочки определялись по формуле {Т} = [Е] {е} - {Т}. После этого напряжения вычислялись обычным образом: = Тав/к. Для расчета межслойных ааз и поперечных азз напряжений в оболочке применялась аппроксимация

^«з = ^«(9ъ 92) П(9з), ^зз = -р- - (р+ - Р-)£ + ^з(9ъ 92) П(9з^

) = 1 - Ч(9з) = 1 - (|)2. а = 1. 2. (.И)

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

= 12 Са+з,а+зе«з, ^з = 6(С^1з + С2зе2з) + 3(р+ + р-). (32)

6. Результаты численного моделирования

Численные расчеты были проведены для осесимметричной оболочки, имеющей форму усеченного параболоида, у которой индекс а = 1 соответствует направлению вдоль оси симметрии, а а = 2 — окружному направлению. Рассмотрен случай локального нагрева оболочки излучением, одновременно воздействовало внутреннее давление ре- = 0.3 МПа. Торцы оболочки свободны от нагрузки. Характеристики композита были выбраны следующими [1]:

С11 = 22 = 20 ГПа, С 1з = ¿гз = 2.2 ГПа, С44 = 4.2 ГПа,

а/ = 2 • 10-6 К-1, аь = 20 • 10-6 К-1, ар = 2 • 10-6 К-1, вР = 0.1, = 0.5, рь = 1200 кг/мз, рр = 2000 кг/мз, ^ = 0.15,

Лз = 0.3 Вт/(м • К), Л1 = 0.8 Вт/(м • К) (при В = 293 К), К1 = 10-23 с,

5 = 50, оь =1.3 кДж/(кг • К), 0р = 0.6 кДж/(кг • К), = 3.6 • 105 кг/(мз • с), К/Я = 8300 К, Г = 0.65, Де° = 1 МДж/кг.

Длина оболочки Ь = 3 м, максимальный радиус Я =1 м, толщина к = 5 • 10^ м, радиус пятна нагрева 8 = 0.02 м. Максимальная температура нагрева составляла 580 оС.

На рис. 1 показаны результаты расчета температурного поля на внешней поверхности оболочки в виде двумерных графиков В(^1, д2, к/2, £) распределения температуры по координатам ^ = 91/Ь, д2 = 92/Ь в моменты времени £ = 18 и 26 с. В рассмотренном случае плотность теплового потока 9д вне зоны эффективного нагрева быстро убывала, а конвективный нагрев отсутствовал, поэтому расчет тепломассопереноса проводился только с помощью решения локальной задачи для области У1.

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

Рис. 1. Распределение температуры в плоскости координат 51, 52 в моменты времени г = 18 с (а) и г = 26 с (б)

к тыльной поверхности оболочки. При достижении условной температуры начала термодеструкции в том слое, до которого дошел температурный фронт, резко интенсифицируется термодеструкция полимерной матрицы. В результате в композите возникает зона интенсивной термодеструкции, которая распространяется в глубь толщины оболочки. При £ =16 с уровень температуры примерно в 3/4 толщины оболочки превышает условную температуру начала термодеструкции материала.

В области термодеструкции содержание полимерной фазы композита резко снижается, а доля пиролитической фазы и пористости, наоборот, резко возрастает (рис. 3). Поровое давление газов, появление которых также обусловлено термодеструкцией, тоже в указанной зоне резко возрастает, достигая максимума в 65 атм (см. рис. 2) на некоторой поверхности q3 = ^ оболочки, причем это значение ^з с течением времени смещается по направлению к ненагреваемой поверхности оболочки.

На рис. 4 показано распределение температурной деформации е2 по толщине оболочки в различные моменты времени в точке Т0. Эти распределения в области умеренного нагрева имеют положительные значения, поскольку температурные деформации в данном случае вызваны только тепловым расширением композита, а в зоне термодеструк-

. ООО

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

Рис. 2. Распределение температуры (а) и порового давления (б) по толщине оболочки для центра пятна нагрева в различные моменты времени с (цифры у кривых)

</>ь 0.50

0.48

0.46

0.44

^16"

\20

28^ 24' \

-2.5

0

а

Ч>9

0.3

0.2

0.1

28 ^ / /

24 /

/ 20 16 у

2,5 <73, мм

-2.5

0 б

2.5 мм

Рис. 3. Распределение концентраций полимерной фазы (а) и пористости (б) по толщине оболочки в различные моменты времени

£2, 0.2

0

-0.2 -0.4 -0.6 -0.8 -1.0

/8 .

-2.5 124 1б\

20\

24\

28 \

<Й, мм

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

Картина распределения прогиба Ш вдоль координатной линии оболочки для моментов времени £ = 12 и 28 с приведена на рис. 5. Зависимость W(51) имеет волнообразный вид, характерный для краевого эффекта в оболочках. Прогиб достигает максимума на некотором расстоянии от края оболочки. При £ = 28 с, когда интенсифицируется процесс термодеструкции, появляется новый локальный максимум прогиба в центре пятна нагрева оболочки, направленный в сторону от источника нагрева, который по абсолютной величине превышает значение локального максимума, обусловленного краевым эффектом.

Картина распределения тангенциального напряжения а2 в срединной поверхности осесимметричной оболочки обусловлена двумя факторами: внутренним давлением и локальным нагревом, и существенно меняется с течением времени нагрева (рис. 6). При £ < 2 с на фоне не зависящего от поля напряжения а2 возникает двойной максимум растягивающих напряжений в зоне пятна нагрева. При дальнейшем нагреве эти максимумы "разделяются", но остаются в окрестности зоны нагрева. Интенсификация процесса термодеструкции при £ > 20 с приводит к тому, что локальные максимумы значительно увеличиваются (до 284 МПа) и превышают величину предела высокотемпературной прочности композита при растяжении [1], что вызывает локальное растрескивание композита с ориентацией микротрещин в зоне нагрева, перпендикулярной окружному направлению.

Картина распределения касательных напряжения а13 содержит многие из приведенных выше эффектов. На рис. 7 показан график зависимости напряжения от времени в семи контрольных точках: точка 1 (То) находится в центре пятна нагрева, точки 2 — 4 расположены вдоль луча, выходящего из точки центра пятна в направлении оси Од1, а точки 5 — 7 — вдоль линии Од2. С течением времени во всех точках вначале реализуются положительные напряжения а13, а затем функция а13(£) меняет знак и ее абсолютные значения увеличиваются. В центре пятна нагрева при £ > 26 с касательные напряжения еще раз меняют знак и становятся положительными. В целом наблюдается эффект резкого увеличения касательных напряжений в момент времени, когда в рассматриваемую точку приходит фронт интенсивного процесса термодеструкции. Результаты расчетов показывают, что при сравнительно невысоком уровне лазерного нагрева наиболее вы-

IV, мм УУ, мм

Рис. 5. Распределение прогиба Ш вдоль координатных линий в моменты времени Ь = 12 с (а) и 28 с (б) для д2/Ь = 0 (1), 0.25 (2), 0.27 (3), 0.53 (4), 0.79 (5)

2.5- Ю-1 2.8- 10"1

Рис. 6. Распределение тангенциального напряжения (ГПа) на срединной поверхности осе-симметричной оболочки в моменты времени Ь = 20 с (а) и 28 с (б)

ffis, ГПа

7/

5 9 5

<3

\

4-Л x

\

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

сокими являются межслойные напряжения а\3, максимальные значения которых в зоне пятна нагрева составляют 0.48 МПа. Учитывая, что прочность композита при межслой-ном сдвиге для рассматриваемой температуры равна примерно 0.5 МПа [1], указанные напряжения будут приводить к расслоению композитной оболочки. Этот вывод согласуется с результатами расчетов [2], проведенных для плоской пластины.

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

Список литературы

[1] Димитриенко Ю.И. Механика композиционных материалов при высоких температурах. М.: Машиностроение, 1997. 362 с.

[2] Dimitrienko Yu.I. Termomechanics of Composites under High Temperatures. Dordrecht; Boston; London: Kluwer Acad. Publ., 1999. 347 p.

[3] Dimitrienko Yu.I. A structural thermomechanical model of textile composite materials at high temperatures // Composite science and technologies. 1999. Vol. 59. P. 1041-1053.

[4] Dimitrienko Yu.I. Modelling of mechanical properties of composite materials under high temperatures. Part 3. Textile Composites // Intern. J. of Appl. Composite Materials. 1998. Vol. 5, No. 4. P. 257-272.

[5] Dimitrienko Yu.I. Thermomechanical behaviour of composites under local intense heating by irradiation // Composites. Part A. 2000. Vol. 31. P. 591-598.

[6] Попов Б.Г. Расчет многослойных конструкций вариационно-матричными методами. М.: Изд-во МГТУ им. Н.Э. Баумана, 1993. 297 с.

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

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

Поступила в 'редакцию 20 января 2011 г., с доработки —11 ноября 2011 г.