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

УДК 536.8

Численное решение сопряженной задачи газодинамики и теплообмена для воздухозаборной решетки с противообледенительной системой

© Ю.И. Димитриенко, М.Н. Коряков, В.Ю. Чибисов МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

Работа посвящена моделированию сопряженных процессов обтекания потоком холодного воздуха и теплообмена между потоком и корпусом воздухозаборной решетки с противообледенительной системой нагрева. Моделирование осуществляется с помощью программного комплекса Sigma, разработанного на кафедре «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Проведено параметрическое численное исследование режимов обтекания и теплообмена при различных скоростях и температурах потока, в результате которых определены мощности системы обогрева, обеспечивающие поддержание температуры корпуса решетки в заданном режиме. Разработанная методика и результаты численного моделирования могут быть использованы для проектирования противообледенительных систем, в том числе для морских судов, эксплуатирующихся в арктических условиях.

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

Введение. Для морских судов, осуществляющих плавание в высоких арктических широтах, существует проблема борьбы с нарастанием льда на различных судовых конструкциях [1-4]. Особую проблему составляет обледенение воздухозаборных решеток (ВЗР) судовой системы вентиляции. Обледенение решеток сужает поперечное сечение воздухозаборных трактов, уменьшает количество воздуха, поступающего на различные судовые системы, и в отдельных случаях может вывести их из строя. Для борьбы с обледенением ВЗР разрабатывают специализированные противообледенительные система обогрева (ПСО) [3-5]. Вопросы моделирования процессов обледенения рассматривались в [6], однако задача совместного моделирования процессов обтекания ВЗР с ПСО холодным потоком воздуха и процессов теплообмена в этих системах не рассматривались. Настоящая работа посвящена моделированию сопряженных процессов газодинамики и теплообмена в конструкциях ВЗР с ПСО.

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

Рис. 1. Фрагмент ВЗР с ПСО

Для расчета теплообмена в конструкции жалюзи с внутренним обогревом от электронагревателей необходим расчет конвективного теплового потока, отводимого от поверхности жалюзи при воздействии холодного воздушного потока, а также расчет коэффициента конвективного теплообмена между поверхностью корпуса ВЗР и холодным воздушным потоком. Для моделирования процессов теплообмена на ВЗР рассматривается постановка сопряженной задачи газодинамики воздушного потока, обтекающего конструкцию, и теплопроводности стенки конструкции с условиями внутреннего электроподогрева. Эта задача в общей постановке состоит из системы уравнений динамики линейно-вязкого теплопроводного газа (образованной из уравнения неразрывности, уравнений движения, уравнения энергии и определяющих соотношений), описывающей движение холодных воздушных масс на поверхности жалюзи ВЗР [7-9]:

[^У-ру = 0,

ы

[ ^ + У-(ру ®у + рЕ -Т ) = 0, Ы К }

дрЕ + У-((рЕ + р )у-Т-у + я) = 0, X е¥,

I I?

р = Яре, е = суе, Е = е + |у| /2, Ту = д1(У-у)Е + д2(У®у + У®уг ),

Я = -АУе, X е¥1,

а также уравнения теплопроводности корпуса ВЗР 2

(1) (2)

д0

p.c. f = V-(A, . V0), x eF2, (3)

где p — плотность газа; t — время; E — плотность полной энергии

I |2 V

газа: E = су0 + —^~; cv — теплоемкость газа при постоянном объеме;

Л I 12 i

0 — температура газа; v = vvi — квадрат модуля скорости; p — давление; R — газовая постоянная (R = R/ц, ц — молекулярная масса газа; R — универсальная газовая постоянная); E — метрический тензор; Tv — тензор вязких напряжений в газе; q — вектор потока тепла;

2

ц1з ц — коэффициенты вязкости газа (полагаем далее ^ = -—); X —

коэффициент теплопроводности газа; p., c., X. — плотность, удельная теплоемкость и коэффициент теплопроводности материала корпуса ВЗР; V — набла-оператор Гамильтона [10], V1 — условная область воздушного потока, обтекающего ВЗР; V2 — область конструкции ВЗР. Рассмотрим 4 случая граничных условий для системы уравнений

(1)-(3).

1. На границе Е1 раздела областей V1 и V2, представляющей собой твердую непроницаемую стенку (поверхность корпуса ВЗР), ставятся условия прилипания, баланса теплового потока и равенства температур воздушного потока и твердой стенки:

v = 0, XV0|E • n = X„V0|E • n + 8.G0:, 0|E = 0|E . (4)

2. На дозвуковой границе E2 входа воздушного потока в область V1 задаются следующие условия:

p = Р* v = v^ e = 0e. (5)

3. На дозвуковой границе E3 выхода воздушного потока из области V1 задаются условия

dv д0

P = Pe, = 0, Т" = 0, (6)

дп дп

dv

где — = n • V ® v — нормальная производная вектора скорости.

дп

4. На внутренней (нагреваемой) границе Е4 твердой стенки задаются условия теплового притока

-X,V0|S • n = qw, (7)

1 ¿-'А

где qw — заданный тепловой поток, подводимыи за счет электронагрева; е^ — интегральный коэффициент теплового излучения твердой поверхности; о — коэффициент Стефана—Больцмана. Начальные условия к системе (1)-(3) имеют вид

t = 0: р(0,х') = р0, V(0,х') = 0, Е(0,х') = Су0О,

(8)

где р0, 00 — заданные значения.

Метод численного решения сопряженной задачи газодинамики и теплообмена. Численное решение сопряженной задачи (1)-(8) осуществлялось в модельной двумерной постановке — рассмотрено нормальное сечение одной опорной стойки ВЗР ПСО (рис. 2), в которой области решения У1 и У2 ограничены двумя концентрическими эллипсам.

Рис. 2. Модель опорной стойки ВЗР ПОС, обтекаемой воздушным потоком, примененная для решения сопряженной задачи в двумерной постановке

Для численного решения сопряженной задачи был применен следующий метод: вводился итерационной цикл по «медленному» времени t = I/10, соответствующему процессу распространения тепла в стенке конструкции ВЗР, где t0 — характерное время нагрева конструкции. Внутри этого цикла введено «быстрое» время т = t / где tу — характерное время установления течения газового потока. На каждом п-м шаге итерации по медленному времени ^ п) расчет осуществлялся в четыре этапа: 1) выбиралась температура поверхности конструкции на п-м шаге, с этой температурой методом установления решались системы уравнений газодинамики (1), (2) с граничными условиями (5), (6), а вместо (4) задавались условия

V = 0, 0,

0

w( п)'

(9)

2) в результате вычислялась температура в области У1 газа, в том числе температура 0 на внешней условной границе пристеночного вязкого

1-

газового слоя, а также тепловой поток от воздушной среды к твердой поверхности конструкции q- = -А,У0„ • п и коэффициент теплообмена а =-^-; 3) составлялось уравнение теплового баланса на твер-

0е -0^(п)

дой поверхности конструкции ВЗР (второе уравнение системы (4)), которое рассматривалось как нелинейное алгебраическое уравнение для вычисления температуры 0^ поверхности на следующем шаге итерации

а(0 -0(+п) = а (0 ( +1)-00) + в а04( +1). (10)

V е п+1)/ 5 V п+1) 0/ 5 п+1) V у

Здесь использованы уравнения Ньютона для тепловых потоков

q- =а (0e - 0п+1) ) и q+=аs ((п+1) - 00 ) , где — коэффициент теплообмена в твердой стенке, способ вычисления которого изложен ниже; 4) решалось уравнение теплопроводности (3) с граничным условием заданной температуры 0|Е = 0^ п+1). Далее осуществлялся переход на

следующий (п + 1) шаг итерационного цикла.

Уравнение теплопроводности (3) с граничными условиями (7) и 0|Е = 0№( п+1) решалось конечно-разностным методом в криволинейных координатах, одна из которых — X1 — связана с обтекаемой поверхностью корпуса ВЗР, а вторая — X2 — ориентирована по нормали к этой поверхности. Вводятся 10 — характерное значение координаты X1, и h — координаты X2, в качестве которого выбирается толщина конструкции ВЗР, а также A0 — характерные значения коэффициентов квадратичной формы срединной поверхности оболочечной конструкции жалюзи, температуры 00, а также соответствующие им безразмерные величины

X = XlA0/lo, X2 = X2/h, 4 = А, ё = 0-, дк (11)

A1 00

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

д0 Кв

2 ( д ( 1 д0 ^

дt А

1 _д0 А дХ

удХ1 к А1

+К А0г, (12)

0 дХ22 4 7

17 ^А h

где К0 =-2 — критерий Фурье, р = —.

Р^ К

Граничные условия (4), (7) и начальные условия (8) в безразмерном виде следующие:

" 1 ' ' " 1: í-»

2 dX„

X2 -- : 0 _ 0w(и+1)' X2 : ,

(13) t - 0: 0-0о.

Для численного решения (12), (13) применялся конечно-разностный метод в сочетании с пошаговым методом линеаризации. Для решения разностных систем уравнений был использован метод скалярной прогонки по координатным направлениям.

Коэффициент теплообмена а, на твердой стенке рассчитывался с помощью специального метода, предложенного в [11], согласно которому численно-аналитическое решение уравнения теплопроводности (12) находится только для главных членов теплового потока в направлении по нормали к нагреваемой поверхности с граничным условием в виде заданной температуры поверхности:

бё ,б20 -— -F0^=-, 0<X2 < 1; d t dX22 2

2 - (14)

1 — — — 1 d0 —

X2--: 0-0,, X2---: -^L- - 0, t - 0: 0-1.

2 2 s 2 2 dX2

В силу линейности задачи (14) ее решение — безразмерная темпе— — 50 —

ратура 0(X2, t) и тепловой поток q (X2, t) являются линейны-

dX2

ми функциями от входных данных задачи — от температуры внеш—

ней поверхности 0 , тогда значение теплового потока qs - (0,1)

dX2

на поверхности контакта газа и конструкции среды в момент времени t -1 можно представить в виде qs - g(Fo)(0, -1), где g(Fo) — некоторая функция от параметра Фурье, которую находим из формулы g(Fo^)-^^(0,1)/(0, -1). Возвращаясь к размерным величи-

dX2

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

а,(15) h

Для решения задачи газодинамики на одном шаге итерации по «медленному» времени применялся метод установления по «быстро-

му» времени т с использованием модели многомерного пограничного слоя [12], благодаря которой уравнения течения идеального газа и вязкого газа также разделяются. Решение уравнений идеального газа ищется во всей области V1 течения газового потока с граничными условиями непроницаемости на твердой стенке, затем полученное решение идеального потока на твердой стенке для плотности, касательных компонент скорости и температуры: pb, vbT = vb • т, 0b переносится на внешнюю поверхность условного пограничного слоя, на которой формулируются следующие условия для системы уравнений (1)-(3):

^ : Р = Р^ v • п = 0, V • т = vbт, 0 = 0^ (16)

где т — единичный касательный вектор в поверхности пограничного слоя. Численный метод решения задачи газодинамики (1)-(3) для идеального газа на основе схем TVD описан в работах [7, 13, 15]. Решение системы уравнений (1), (2) c условиями (5)-(8), (9) в области V1 проводится до установления.

Численное моделирование по разработанному алгоритму осуществлялось с помощью программного комплекса Sigma, разработанного на кафедре ФН-11 МГТУ им. Н.Э. Баумана [15-17].

Результаты численного моделирования. Были проведены несколько серий вариантного численного моделирования прогрева конструкций створок жалюзи ВЗР с ПОС. Изменялись следующие параметры:

температура холодного воздуха в диапазоне 0е = -10...-50 °С.

скорость обтекающего потока холодного воздуха в диапазоне 5...20 м/с;

мощность электронагревателя W, приходящегося на один элемент жалюзи ВЗР, в диапазоне 100.1000 Вт.

В расчетах варьировались значения плотности теплового потока qw = yW/S, где у — тепловой эквивалент передачи электрической энергии в тепловую; S — площадь нагреваемой поверхности створок жалюзи. В качестве материалов опорной стойки ВЗР была выбрана сталь со следующими характеристиками: ps = 7,8 г/см3, cs = 0,8 кДж/(кг • К), Xs = 5 Вт/(м • К). В результате численных расчетов было получено значение коэффициента теплообмена на поверхности ВЗР: а = 0,4 Вт/(м • К).

На рис. 3-6 показаны некоторые из полученных результатов численного моделирования параметров воздушного потока, обтекающего конструкцию жалюзи ВЗР с ПОС при температуре холодного воздуха 0е = 50° C, скорости ve = 5 м/с и плотности теплового потока обогревателя qW = 210 Вт/м2.

Рис. 3. Расчетное поле плотности воздушного потока, обтекающего конструкцию опорной стойки ВЗР с ПСО, при температуре холодного воздуха 0 = -50 °С, скорости 5 м/с и плотности теплового потока обогревателя

q¡

ш

= 210 Вт/м2

Рис. 4. Расчетное поле числа Маха воздушного потока, обтекающего конструкцию опорной стойки ВЗР с ПСО, при температуре холодного воздуха 0 = -50 °С, скорости 5 м/с, плотности теплового потока обогревателя

. = 210 Вт/м2

1ш

ргеээиге

1.039+05 1.039+05 1.039+05 1.03©+05 1.039+05 1.049+05 1.049+05 1.049+05 1.049+05 1.049+05 1.049+05

Рис. 5. Расчетное поле давления воздушного потока, обтекающего конструкцию опорной стойки ВЗР с ПСО, при температуре холодного воздуха 0 = -50 °С, скорости 5 м/с, плотности теплового потока обогревателя

qW = 210 Вт/м2

Рис. 6. Расчетное температурное поле воздушного потока, обтекающего конструкцию опорной стойки ВЗР с ПСО, при температуре холодного воздуха 0 = -50 °С, скорости 20 м/с, плотности теплового потока обогревателя

дш = 210 Вт/м2

На рис. 7 показан один из полученных результатов численного вариантного моделирования температурного поля в стенках конструкции опорной стойки ВЗР с ПСО. Результаты расчетов показывают, что температурное поле распределяется очень неравномерно: минимальное значение температуры достигается в передней (наветренной) части стойки, выступающей в сторону набегающего потока, а максимальное значение локализуется в задней части стойки.

Рис. 7. Расчетное температурное поле в конструкции створки опорной стойки ВЗР с ПСО при температуре холодного воздуха 0е = -10 °С, скорости 5 м/с, плотности теплового потока обогревателя = 210 Вт/м2

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

При температуре воздуха 0е = -50 °С и скорости потока \е = 5 м/с для обеспечения положительных значений температуры (0 . = 0,04 °С, 0» тах = 97,1 °С) на корпусе опорной стойки ВЗР необходимы значения плотности теплового потока электронагрева ^ = 670 Вт/м2. Таким образом, установлено, что с понижением температуры окружающего воздуха, возрастает не только мощность электронагрева, необходимая для

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

Таблица

Значения максимальной и минимальной температур в корпусе опорной стойки ВЗР с ПСО при различных отношениях температуры и скорости потока и плотности теплового потока электронагрева

№ Температура, °С, и Плотность теплового потока электронагрева дж, Вт/м2

скорость потока, м/с 210 250 300 350

1 -10 53,4 59,2 62,4 62,6

5 0,9 1,5 4,6 7,8

2 -10 56,4 61,2 68,4 72,6

20 1,9 3,5 5,6 8,8

3 -30 83,2 87,6 94,2 112,2

5 0,9 1,8 2,6 4,6

4 -30 89,5 86,5 93,1 107,5

20 0,9 1,8 2,4 4,2

5 -50 89,5 86,5 93,1 107,5

5 0,9 1,8 2,4 4,2

Выводы. Разработана математическая модель тепловых процессов, проходящих в конструкциях ВЗР с ПСО в условиях воздействия холодных воздушных потоков.

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

Проведены серии численных экспериментов по моделированию процессов обтекания и теплового режима ВЗР с ПСО при различных значениях параметров холодного потока и мощности электронагрева, в результате которых установлено, что с уменьшением температуры холодного воздуха от -10 до -50 °С потребная мощность электронагревателя для обогрева жалюзи и плотность теплового потока возрастает от 210 до 670 Вт/м2, при этом изменение скорости ветра от 5 до 20 м/с сказывается менее существенно на потребной мощности электронагрева.

Исследование выполнено при финансовой поддержке Министерства образования и науки Российской Федерации (номер НИР 1.5433.2011), РФФИ (грант № 12-08-00998-а) и Министерства промышленности и торговли РФ (государственный контракт №12411.1007499.09.062).

ЛИTЕPАTУPА

[1] Захаров Ю^. Судовые установки кондиционирования воздуха и холодильные машины. Ленинград, Cудостроение, 1972, 566 с.

[2] Языков B.H. Теоретически основы систем кондиционирования воздуха. Ленинград, Cудостроение, 1967, 234 с.

[3] Feher L., Thumm M. Design of Avionic Microwave De-Anti-Icing Systems. in M.Willert-Porada(ed), Microwave devices. Springer 792 p., ISBN: 354G432523, pp. 695-7G2. 2GG6.

[4] Petrenko V.F., Sullivan C. Methods and Systems for Removing Ice from Surfaces. Patent 6,653,598 B2, US. 2GG3.

[5] Bhamidipat: M. Smart Anti-Ice Coating. URL: http://www.virtualacquistitionshow case.com/document/44G/briefing

[6] Lozowski E., Szilder K., Makkonen L. Computer simulation of marine ice accretion. Philosophical Trans of the Royal Society. Mathematical, Physical and Engineering Sciences, 2GGG, vol. 358, pp. 2811-2845.

[7] Димитриенко Ю.И., Котенев B.n., Захаров А. А. Метод ленточных адаптивных сеток для численного моделирования в газовой динамике. Москва, ФИЗ-МА1ЛИ1, 2G11, 28G с.

[8] Димитриенко Ю.И. Механика сплошной среды. T. 2: Универсальные законы механики и электродинамики сплошной среды. Москва, Изд-во МГХУ им. Н.Э. Баумана, 2G11, 56G с.

[9] Димитриенко Ю.И. Нелинейная механика сплошной среды. Москва, ФИЗ-МА1ЛИ1, 2GG9, 61G с.

[1G] Димитриенко Ю.И. Механика сплошной среды. Т. 1: Тензорный анализ. Москва, Изд-во МГГУ им. Н.Э. Баумана, 2G11, 463 с.

[11] Dimitrienko Yu.I., Efremov G.A., Chernyavsky S.A. Optimal Designing of Erosion-Stable Heat-Shield Composite Materials. International Journal of Applied Composite Materials, 1997, vol. 4 (1), p. 35-52.

[12] Димитриенко Ю.И., Захаров А.А., Коряков М.Н. Модель трехмерного пограничного слоя и ее численный анализ // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. Спец. вып. 2G11, с. 136-15G.

[13] Димитриенко Ю.И., Коряков М.Н., Захаров А.А., Cыздыков Е.К. Развитие метода ленточно-адаптивных сеток на основе схем TVD для решения задач газовой динамики. Вестник МГТУ им. Н. Э. Баумана. Сер. Естественные науки, 2G11, № 2, с. 87-97.

[14] Димитриенко Ю.И., Захаров А.А., Коряков М.Н. Моделирование газодинамических потоков в каналах сверхзвуковых воздухозаборников на основе модели трехмерного пограничного слоя. Инженерный журнал: наука и инновации, 2G12, вып. 2. URL: http://engjournal.ru/catalog/mathmodel/aero/39.html

[15] Димитриенко Ю.И., Захаров А.А. Автоматизированная система для моделирования газовых потоков методом ленточных адаптивных сеток. Информационные технологии, 2GG9, № 6, с. 12-16.

[16] Димитриенко Ю.И., Минин B.B., Cыздыков Е.К. Моделирование внутреннего тепломассопереноса и термонапряжений в композитных оболочках при локальном нагреве. Математическое моделирование, 2G11, т. 23, № 9, с. 14-32.

[17] Димитриенко Ю.И., Минин B.B., Cыздыков Е.К. Численное моделирование процессов тепломассопереноса и кинетики напряжений в термодеструкти-

рующих композитных оболочках. Вычислительные технологии, 2012, т. 17, № 2, с. 44-60.

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

Статья поступила в редакцию 28.06.2013

Ссылку на эту статью просим оформлять следующим образом:

Димитриенко Ю.И., Коряков М.Н., Чибисов В.Ю. Численное решение сопряженной задачи газодинамики и теплообмена для воздухозаборной решетки с противообледенительной системой. Инженерный журнал: наука и инновации, 2013, вып. 9. URL: http://engjournal.ru/catalog/mathmodel/technic/1116. html

Димитриенко Юрий Иванович родился в 1962 г. Окончил МГУ им. М.В. Ломоносова в 1984 г. Д-р физ.-мат. наук, профессор, заведующий кафедрой «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана, директор Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов» МГТУ им. Н.Э. Баумана (НОЦ «СИМПЛЕКС»), действительный член академии инженерных наук. Автор более 250 научных работ в области вычислительной механики, газодинамики, термомеханики композитов, математического моделирования в науке о материалах. e-mail: [email protected]

Коряков Михаил Николаевич родился в 1987 г., окончил МГТУ им. Н.Э. Баумана в 2010 г. Аспирант кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Автор 10 работ в области вычислительной газодинамики.

Чибисов Виктор Юрьевич родился в 1989 г., окончил МГУПБ в 2011 г. Ведущий инженер Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов» МГТУ им. Н.Э. Баумана. Специалист в области теплотехнических технологий. e-mail: [email protected]