НАУКА и ОБРАЗОВАНИЕ

УДК 539.3

Россия, МГТУ им. Н.Э. Баумана

dimit.bmstu@gmail.com

mkoryakov@bmstu.ru

polardusk@mail.ru

Полимерные композитные материалы до сих пор являются одним из основных типов материалов, применяемых в теплозащитных системах высокоскоростных летательных аппаратов, движущихся в плотных слоях атмосферы Земли, а также других планет (Марса, Венеры). В работах [1-13] была разработана математическая модель процессов внутреннего тепломассопереноса в полимерных теплозащитных материалах при нестационарном нагреве, учитывающая процессы термодеструкции, газообразования и фильтрации газов в порах материала. Модель была применена для расчета одномерных процессов теплопереноса в элементах конструкций с относительно простой геометрической формой: пластинах, цилиндрических оболочках [11-13]. Вопросы моделирования внутреннего тепломассопереноса в композиционных материалах при высоких температурах рассматривались также в современных работах зарубежных авторов [14-26]. В настоящей работе предложен метод численного конечно-элементного расчета процессов внутреннего теплопереноса в трехмерных элементах теплозащитных конструкций из композитов.

Система трехмерных уравнений внутреннего тепломассопереноса в теплозащитных конструкциях.   Под действием высоких температур, обычно свыше 300^оС, в полимерных композиционных материалах происходят внутренние физико-химические превращения – процессы термодеструкции. В результате композит при нагреве представляет собой многофазную систему, в работах [1, 2] была предложена модель 4-х фазной системы, в которой 1-я фаза — наполнитель (волокна, армирующие частицы), 2-я фаза — полимерная матрица; 3-я фаза – твердая пиролитическая фаза, образующаяся из полимерной при высоких температурах, 4-я фаза – газ, также образующийся вследствие термодеструкции исходной фазы 2 и находящийся в порах. Обозначим объемные концентрации фаз как . Эти концентрации фаз  постоянно меняются при нагреве, что приводит к изменению характеристик всего материала: коэффициентов теплопроводности, плотности, теплоемкости. Система уравнений изменения концентраций фаз, а также уравнение энергии всей 4-х фазной системы, образует систему уравнений внутреннего теплопереноса в композите, и записывается в тензорной форме следующим образом [1]:

                                           (1)

где  ,  ,    ,   -  плотности , объёмные доли  и удельные теплоёмкости фаз,  - температура, общая для всех фаз, t– время,- радиус-вектор, - область, занятая композитной конструкцией,  - набла -оператор [27],  - скорости массопереноса (термодеструкции), - теплота термодеструкции, для которых имеют место соотношения:

,  ,  ,  ,

,    .                                                                   (2)

Здесь  - плотность всего композита,  - поровое давление газообразных продуктов термодеструкции,   - удельная теплоёмкость ,   - коэффициент вязкости газа,   - энергия активации терморазложения,   - универсальная газовая постоянная,   - коэффициент газификации.           Плотности всех твердых фаз полагаются постоянными: . Тензор коэффициентов газопроницаемости  и тензор теплопроводности  композита являются функциями от объемных долей фаз и температуры: 

,     ,                                                              (3)

где   , - константы, а  - метрический тензор. Конкретные выражения для  функций  приведены в [1].

На нагреваемой части внешней границы области Vзаданы тепловой поток и давление, на части поверхности  заданы температура и давление, а на оставшейся части  заданы условия герметичности и теплоизоляции:

:       ,                                                   (4)

 :            ,                                                                                (5)

     ,:      ;       ,                                                           (6)

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

Начальные условия к системе (1)-(3) имеют вид

:   ,  ,  .                                                                      (7)

         В итоге получаем постановку задачи внутреннего тепломассопереноса (1)-(7) относительно 3-х неизвестных скалярных функций: , зависящих от .

Вариационная постановка задачи внутреннего тепломассопереноса.Дадимвариационную постановку задачи внутреннего тепломассопереноса (1)-(7). Введем вариации 3-х неизвестных функций: ,  и , которые удовлетворяют граничным условиям 1-го рода на внешних поверхностях:

:        ,                                                                                                (8)

 :            ,                                                                                (9)

      Домножим первое уравнение системы (1) на , второе на , а третье – , и проинтегрируем получившиеся уравнения по области V, в результате получим:

,                                                                                  (10)

                                                    (11)

    (12)

Преобразуем в этих уравнениях интегралы, содержащие ковариантные производные, следующим образом:

 ,                    (13)

,

здесь использована формула Гаусса-Остроградского, а также граничные условия (4)-(6), (8), (9). Подставляя формулы (13) в (11) и (12), получаем искомую вариационную остановку задачи внутреннего тепломассопереноса

,                                                                                  (14)

                                                    (15)

    (16)

Производная по времени в (15) преобразуется следующим образом:

                                                                                (17)

Конечно-элементная аппроксимация задачи тепломассопереноса. Согласно общей концепции метода конечных элементов [28-30], область V решения вариационных уравнений (14)-(16) разобьем на систему большого числа подобластей (конечных элементов) V_e , в каждом из которых решение системы (14)-(16) ищем в виде аппроксимаций по узловым функциям температуры , порового давления  и концентрации полимерной фазы  (значения функций в узлах КЭ):

            ,      ,     ,   (18)

где ,, - координатные столбцы функций форм КЭ [29], все имеющие одинаковую структуру и выражающиеся через барицентрические координаты КЭ . Тогда вариации неизвестных функций имеют вид:

,   , .                                       (19)

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

,                                                                       (20)

где  - координатный столбец дифференциальных операторов, m– размерность пространства, а - матрица, полученная внешним умножением координатных столбцов. Аналогично вычисляем

,                                                                     (21)

где .

Подставляя  выражения (18)-(21) в систему (14)-(16), имеем:

,                                          (22)

        (23)

                                                                (24)

где обозначена функция

             (25)

Вынося вариации неизвестных функций за скобки, преобразуем эту систему к виду

                                                (26)

где обозначены матрицы и координатные столбцы:

  ,  ,               (27)

,  ,

  , ,

 ,     

.

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

   (28)

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

Линеаризованная система (28) может быть записана в обобщенном матричном виде

,                                                    (29)

где введена обобщенная матрица жесткости

,                      (30)

обобщенный столбец неизвестных и - столбец правых частей.

Для всех функций будем рассматривать линейную аппроксимацию с использованием 4-х узлового КЭ, для него барицентрические координаты  имеют вид [29]:

,                                                               (31)

            ,,,

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

.                                    (32)

Тогда матрицы  и  могут быть вычислены без расчета интегралов по КЭ:     , . Для вычисления матриц , ,,,  воспользуемся формулами [29]

                                                            (33)

,      ,      .

Тогда, вычисляя переменные  ,  и аналогичные другие, участвующие в выражениях (27) в матрицах, как средние по 4-м узловым значениям, получаем:

  ,  ,     ,       ,   ,  (34) где обозначены следующие величины

, , ,  ,    ,

,              .                                           (35)

Алгоритм конечно-элементного метода решения задачи внутреннего тепломассопереноса.  Общая схема алгоритма конечно-элементного метода решения задачи внутреннего тепломассопереноса состоит из следующих этапов:

1.Разбивка области на КЭ.

2.Нумерация:  границ области; узлов на границе; КЭ; узлов КЭ; узловых степеней свободы.

3.Задание свойств материалов: констант  и других.

4.Задание граничных условий.

5. Задание начальных условий.

6. Организация итерационного процесса решения задачи, на каждой итерации:

6.1. вычисляются локальные матрицы «жесткости» и матрицы «масс» конечных элементов, для этого используются значения температуры и давления для конечного элемента из предыдущей итерации (на 1 – начальные значения);

6.2 формируются на их основе глобальная матрица «жёсткости» и векторы «узловых нагрузок»;

6.4. формируется разрешающая СЛАУ с учетом граничных условий;

6.5. осуществляется решение СЛАУ методом локального оптимума.

7. Проверяется условие окончания процесса нагрева конструкции.

8. Осуществляется вывод и обработка результатов расчётов.

Хотя все матрицы (27) являются симметричными, но несимметричной является обобщенная матрица линеаризованной системы (30), поэтому для решения системы (29) использовались методы бисопряженного градиента.

Генерирование трехмерной КЭ сетки осуществлялось с помощью функций свободно распространяемой библиотеки tetgen. Программная реализация разработанного конечно-элементного метода решения задачи осуществлялась в программной среде  VS-2010. Был написан специализированный программный модуль FemHMTCom, реализующий конечно-элементный метод решения задачи.  

 Тестовые расчеты.  Для тестирования разработанного программного обеспечения использовалось решение одномерной задачи нестационарной теплопроводности об одностороннем нагреве цилиндра

                                            (36)

со следующими значениями  параметров: ,  ,  , , , . Число КЭ в сетке, построенной для тестовых расчетов составляло 5095, число узлов в сетке – 1137. Тестирование проводилось путем сравнения расчетов на основе метода МКЭ и расчетов, полученных с помощью неявной конечно-разностной схемы, для решения которой использовался метод прогонки. На рисунке 1 показаны распределения температуры по толщине цилиндра, полученные с помощью обоих методов.  Как видно, расчетные кривые, полученные по обоим методам достаточно хорошо совпадают, волнистый характер графиков температуры, полученных с помощью МКЭ, обусловлен первым порядком аппроксимации температуры в конечном элементе и достаточно грубой КЭ-сеткой.  При использовании более мелких сеток характер кривых температуры становится более гладким, однако при этом существенно возрастает время расчета задачи, что связано с   трехмерностью задачи.

Рис. 1 Результаты тестовых расчетов температуры в цилиндре (гладкие кривые – расчет по конечно-разностному методу, волнистые – МКЭ-расчет), цифры – у кривых – время нагрева t ( с).

Численное конечно-элементное моделирование  тепломассопереноса в теплозащитной конструкции.  МКЭ расчеты проводились для трехмерной цилиндрической конструкции из полимерного композиционного материала при внешнем нагреве. Геометрические параметры конструкции были следующими: длина цилиндра  L=0,15 м, внутренний радиус R1=0,12 м, толщина цилиндра h=0,02 м. КЭ сетка состояла из 49160 элементов и 11848 узлов, элемент  использованной КЭ сетки для конструкции показан рис. 2. Характеристики полимерного композиционного материала (был выбран стеклопластик на эпоксифенольном связующем с углеродными волокнами) взяты из работы [1], температура внешнего потока , давление внешнего потока   .

Результаты расчетов показаны на рисунках 3…5. На рисунке 3 приведены распределения температурного поля  по толщине цилиндра (r– радиус) для различных моментов времени tпрогрева.  С увеличением времени нагрева конструкции происходит возрастание температуры как на внешней нагреваемой поверхности, так и внутри самого материала. На рисунке 4 показаны распределения концентрации полимерной фазы композита  по толщине конструкции. Поскольку на внешней поверхности и в подповерхностном слое конструкции температура композита превышает температуру начала интенсивной термодеструкции полимерной матрицы (примерно ), то в зоне максимального нагрева реализуется процесс термодеструкции матрицы и содержание полимерной фазы снижается. С увеличением времени нагрева  зона уменьшенного содержания полимерной фазы продвигается внутрь конструкции. 

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

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

Рис. 2. Элемент КЭ сетки, использованной для моделирования тепломассопереноса в цилиндре

Рис. 3. Распределение температуры  по толщине цилиндра.

Рис. 4. Распределение концентрации полимерной фазы  по толщине цилиндра.

Рис. 5. Распределение порового давления pпо толщине цилиндра.

Исследование выполнено при поддержке Министерства образования и науки Российской Федерации (Соглашения  №14.B37.21.0448, 14.B37.21.1869, НИР № 1.5433.2011), грантов Президента РФ (МК-6421-2012-9, МК-765.2012.8, МК-6573.2013.3) и гранта РФФИ №12-08-00998.

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

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

2.               Dimitrienko Yu.I. Thermal stresses and heat mass-transfer in ablating composite materials // Int. Journ. of Heat Mass Transfer. 1995. Vol. 38, no. 1. P. 139-146.

3.               Dimitrienko Yu.I. Thermal Stresses in Ablative Composite Thin-Walled Structures under Intensive Heat Flows // Int. Journal of Engineering Science. 1997. Vol. 35, no. 1. P. 15-31.

4.               Dimitrienko Yu.I. Heat- Mass-Transport and Thermal Stresses in Porous Charring Materials // Journal of Transport in Porous Media. 1997. Vol. 27, no. 2. P. 143-170.

5.               Dimitrienko Yu.I. Modelling of Mechanical Properties of Composite Materials under High Temperatures. Part 1. Matrix and Fibres // Int. Journal of Applied Composite Materials. 1997. Vol. 4, no. 4. P. 219-237.

6.               Dimitrienko Yu.I. Modelling of Mechanical Properties of Composite Materials under High Temperatures. Part 2. Unidirectional Composites // Int. Journal of Applied Composite Materials. 1997. Vol. 4, no. 4. P. 239-261.

7.               Dimitrienko Yu.I. Thermomechanical Behaviour of Composite Materials and Structures under High Temperatures. Part 1. Materials Composites. Part A // Applied Science and Manufacturing.  1997. Vol. 28A. P. 453-461.

8.               Dimitrienko Yu.I. Thermomechanical Behaviour of Composite Materials and Structures under High Temperatures. Part 2. Structures. Composites. Part A // Applied Science and Manufacturing.  1997. Vol. 28A. P. 463-471.

9.               Dimitrienko Yu.I. Effect of Finite Deformations on Heat-Mass Transfer in Elastomer Ablating Materials // Int. Journ. of Heat Mass Transfer. 1997. Vol. 40, no. 3. P. 699-709.

10.            Dimitrienko Yu.I. Internal Heat-Mass-Transfer and Stresses in Thin-Walled Structures of Ablating Materials // Int. Journ. of Heat Mass Transfer. 1997. Vol. 40, no. 7. P. 1701-1711.

11.            Димитриенко Ю.И., Минин В.В., Сыздыков Е.К. Численное моделирование процессов тепломассопереноса и кинетики напряжений в термодеструктирующих композитных оболочках // Вычислительные технологии. 2012. Т. 17, № 2. С. 44-60.

12.            Димитриенко Ю.И., Минин В.В., Сыздыков Е.К. Моделирование внутреннего тепломассопереноса и термонапряжений в композитных оболочках при локальном нагреве // Математическое моделирование. 2011. Т. 23, № 9. С. 14-32. 

13.            Димитриенко Ю.И., Минин В.В., Сыздыков Е.К. Моделирование термомеханических процессов в композитных оболочках при локальном нагреве излучением // Механика композиционных материалов и конструкций. 2011. Т. 17, № 1. С. 71-91.

14.            McGurn Matthew, DesJardin P., Dodd A. Numerical simulation of expansion and charring of carbon-epoxy laminates in fire environments // International Journal of Heat and Mass Transfer. 2012. V. 55. P. 272-281.

15.            Mouritz A.P., Mathys Z. Post-fire mechanical properties of marine polymer composites // Composite Structures. 1999. V. 47. P. 643-653.

16.            Mouritz A.P., Mathys Z. Post-fire mechanical properties of glass-reinforced polyester composites // Composites Science and Technology. 2001. V. 61. P. 475-490.

17.            Gibson A.G., Wright P.N.H., Wu Y.S., Mouritz A.P., Mathys Z., Gardiner C.P. The integrity of polymer composites during and after fire // Journal of Composite Materials. 2004. V. 38. P. 1283-1307.

18.            Kardomateas G.A., Liu L., Birman V., Holmes J.W., Simitses G.J. Thermal buckling of a heat-exposed, axially restrained composite column // Composites Part A (Applied Science and Manufacturing). 2006. V. 37. P. 972-980.

19.            Henderson J.B., Wiebelt J.A., Tant M.R. A model for the thermal response of polymer composite materials with experimental verification // Journal of Composite Materials. 1985. V. 19. P. 579-595.

20.            Henderson J.B., Wiecek T.E. A mathematical model to predict the thermal response of decomposing expanding polymer composites // Journal of Composite Materials. 1987. V. 21. P. 373-393.

21.            Sullivan R.M., Salamon N.J. A finite element method for the thermochemical decomposition of polymeric materials–1. theory // International Journal of Engineering Science. 1992. V. 30. P. 431-441.

22.            Sullivan R.M., Salamon N.J. A finite element method for the thermochemical decomposition of polymeric materials–2. carbon phenolic composites // International Journal of Engineering Science. 1992. V. 30. P. 939-951.

23.            Sullivan R.M. A coupled solution method for predicting the thermostructural response of decomposing, expanding polymeric composites // Journal of Composite Materials. 1993. V. 27. P. 408-434.

24.            Pering G.A., Farrell P.V., Springer G.S. Degradation of tensile and shear properties of composites exposed to fire or high temperature // Journal of Composite Materials. 1980. V. 14. P. 54-68.

25.            McManus H.L.N., Springer G.S. High temperature thermomechanical behavior of carbon-phenolic and carbon-carbon composites. 1. Analysis // Journal of Composite Materials. 1992. V. 26. P. 206-229.

26.            McManus H.L.N., Springer G.S. High temperature thermomechanical behavior of carbon-phenolic and carbon-carbon composites. 2. Results // Journal of Composite Materials. 1992. V. 26. P. 230-255.

27.            Димитриенко Ю.И. Тензорное исчисление. М.: Высшая школа, 2001. 575 с.

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

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

30.            Алфутов Н.А., Зиновьев П.А., Попов Б.Г. Расчет многослойных пластин и оболочек из композиционных материалов. М.: Машиностроение, 1984. 264 с.