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

УДК 532.5:621.694

Золотоносов А.Я. - аспирант Хайруллин М.Р. - инженер

Золотоносов Я.Д. - доктор технических наук, профессор

E-mail: [email protected]

Казанский государственный архитектурно-строительный университет АЛГОРИТМ ЧИСЛЕННОЙ РЕАЛИЗАЦИИ СОПРЯЖЕННОЙ ЗАДАЧИ ТЕПЛООБМЕНА НА ОСНОВЕ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ АННОТАЦИЯ

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

КЛЮЧЕВЫЕ СЛОВА: метод конечных элементов, разбиение области, базисные функции элемента, ансамблирование.

NUMERICAL IMPLEMENTATION ALGORITHM OF CONJUGATE HEAT INTERCHANGE PROBLEM BASED ON FINITE-ELEMENT METHOD

Zolotonosov A.Ya. - post-graduate student

Khairullin M.R. - engineer

Zolotonosov Ya.D. - doctor of technical sciences, professor

Kazan State University of Architecture and Engineering ABSTRACT

Numerical decision algorithm of conjugate heat interchange problem in revolving channels with irregular geometry are developed. It has allowed to determine speed components and temperature in the flowing part of channel depending on spin, Reynolds and Peclet numbers. Obtained results have increased representation about fluid dynamics and heat interchange processes of the viscous fluid flow in channels with irregular shape.

KEYWORDS: finite-element method, partition of region, basic functions of element, combining.

В работах [1-3] была предложена математическая модель задачи сопряженного теплообмена во вращающихся вокруг своей оси криволинейных каналах с оребренной и неоребренной проточной частью типа «конфузор-диффузор», представляющая собой систему дифференциальных уравнений в частных производных второго порядка. В качестве условий однозначности задаются начальные распределения скорости, давления, температуры на входе в канал и граничные условия на стенках канала и ребра.

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

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

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

О = {(р,~2,г )|о < р < р0,0 £ г < 1, 0 < г < 1}

( е)

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

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

а) б)

Рис. 1. Проекции разбиения области в осевом и радиальном сечении: а) осевое сечение; б) радиальное сечение

Внутри каждого элемента е проводится замена неизвестных величин /, О, Н, Р, ґж, ґс, ґр ее пробной аппроксимацией [4]:

/') ( V,г,г ) = Іи(*Ф1( V,г); О'*1 (V,г,г ) = Іп,<*1Ф(*1 (V,г);

Н'■1 (V, г ) Л *1Ф ;*1 (V, г, г), С1 (V, г, г) = £е>Ф ;*1 (,, г );

*1 ( V. г, г ) = І С1Ф ;*1 ( V, г ); /;1 (V, г, г) = I «> ( *1 (V, г г);

где Ф(*1 - базисные функции элемента е , удовлетворяющие, в каждой внутренней точке элемента, условиям:

І Ф( *1 (ф, г, г) = 1, ^ (р., , г )Л°’ г ~]

1=1 1°, ,* ]

Для получения системы алгебраических уравнений (дискретного аналога исходных дифференциальных уравнений} используется метод взвешенных невязок Галеркина [5]:

}ФХ ( X ) й П = °,

где X - искомая величина, Ь (X1 - дифференциальное уравнение, определяющее X, Ф, -базисные функции.

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

Рассмотрим элемент, ограниченный линиями:

р = рт; р = рт + ар; г г = гт + Аг; Г = г(г);Г = г(г) •

Применим линейное отображение [6]:

Х =2 (р — рт )/АР — 1;

Ф :=<Т1 = 2 (г — ^ )/Аг — 1;

С = 2(Г — Г1 (г))/( г2 (г) — Г1 (г) ) — 1.

Тогда отображенный элемент будет единичным кубом (рис. 2):

— 1 < ч < 1, — 1 < С < 1, — 1 < X < 1.

Функции формы в локальной системе координат (%, Ч, С ) представим в виде [4]:

ФЬ = 8 (1 + Х0)(1 + Ч0)(1 + С 0)

Рис. 2. Общий вид одного элемента при линейном отображении Ф

Введение локальной системы координат позволяет использовать квадратурные формулы Гаусса при вычислении коэффициентов матрицы жесткости.

Рассмотрим интегральную матрицу вида

1 [ р уу,

где [ Р ] зависит от Фь (Р,г,Г) или их производных по глобальным координатам

Эф, Эф Эф

дф дг дг

Согласно правилу дифференцирования сложной функции имеем

дФь дФь дV дФь дг дФь дг

—- = —-— + — -------------+ —- —;

д£ др д£ дг д£ дг д£

дФь дФь дV дФь дг дФь дг

—— = — ---------+ — ------+ — -------;

дц дV дц дг дц дг дц

дФр дФр дV дФр дг дФр дг

д£ ~ ду д£ дг д£ дг д£'

Искомые производные по глобальным координатам могут быть получены отсюда по правилу:

1 Э в •со 1 ЭФЬ Эф Эг Эг

Эф ЭХ ЭХ ЭХ ЭХ

дФр т-і дфр , где I = Эф Эг Эг

= I —

Эг Эц Эц Эц Эц

ЭФь ЭФЬ Эф Эг Эг

Эг _ЭС _ Ьс ЭС ЭС 1

при условии, что матрица Якоби преобразования координат I не вырождена.

При этом переход к интегрированию по новым координатам осуществляется по формуле dV = det (I) dX dh dZ , где V - объем произвольного конечного элемента.

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

} I }[>(X,Л,С)]det(г)dg dh dZ ,

-1 -1 -1

которая интегрируется с помощью квадратурной формулы Г аусса [6]:

1 1 }[>(X,л,С)]det(I)dX dh dС = Х X nWkF(,лj,Ск).

-1-1 -^ ^ 1=0 j=0 k=0

При решении уравнений движения для исключения безразмерного параметра давления P из числа неизвестных используем метод штрафа. Для этого правую нулевую часть уравнения неразрывности заменяем произведением давления на малый штрафной параметр е (е ®0, DivV ®0) [7]:

/ Лгг Лгг \

N ЭG f

+----------+ з = - Ре

Я Эф г

д/ дИ _дИ

-^-- Я'(г)— - Я—

дг ^ дг дг у

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

д1 г ю Л*1 ^

д2Х

с1У = |[ф| — 1Г&

дх

I

-(IV,

V дх я дх у дх дх

где 5” - граница области, I - внешняя нормаль к границе. А в силу того, что базисные функции

Фь обращаются в нуль на границе элемента, то интеграл по границе области 5 равен нулю.

Получившуюся систему уравнений можно представить в виде:

К ■

и В К11 К12 К13 К14 _

V В2 % где К = К 21 К 22 3 2 К 4 2 К

< > = <

w В3 К31 2 3 К 3 3 К 4 3 К

Ґ ж В4 К41 К42 К43 1 4 ьГ

где К - локальная матрица элемента.

Для получения глобальной матрицы жесткости, построенной для всей расчетной области

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

Для этого строкам и столбцам локальной матрицы элемента приписываются номера соответствующих узлов в глобальной матрице жесткости. Порядок расположения узлов соответствует обходу элемента по часовой стрелке, начиная с /-го узла [5]. Т.к. мы имеем систему из 4 уравнений с 4 неизвестными, то глобальная матрица будет размером 4К X 4К, где N - число узлов в расчетной области. Следовательно, каждый узел имеет 4 степени свободы, в

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

К =

К„ К12 К13 К14

К 2! К 22 К 23 К 24

К31 К32 К33 К34

К41 К42 К43 К44

где кажый К. представляет собой матрицу размера NXN.

1

1 Кг со

II ^ . кг1 ь- 22 ь 28

_ ^81 82 со

^2 —

% = е

N

К1]

ц

ъгу

N1

е1 К ' 2 у К ЙЙ1 £ і к 1е2 и 1е3

..К н ещ уки...К1} а И е^е2 -к,...К13 А 12 ' де8 ■ ^18 ■ .

. _К ‘ н «24 \-кп...К н 21 е2^2 22 ^8 ■ к28..

-К- % н£,...Г ч 51 ^8^2 -кю...К ^ ^88 ■ ■

/ $ < К ^ А "у

N

К:;

^ ХУ

Полученная глобальная система уравнений движения и энергии является нелинейной.

Кіі К12 К13 К14 и В1

К 21 К31 К 22 К32 К 23 К33 К24 К34 V w > - < м В3

К41 К42 К43 К44 Ґ В4

> = 0,

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

Все полученные интегралы вычисляются численно на равномерной сетке 2 X 2 X 2 по квадратурным формулам Гаусса.

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

Представленный алгоритм, реализованный в виде программы в системе Ма1:ЬаЬ, позволил получить профили и значения безразмерных осевых, радиальных и окружных компонентов скоростей в конфузорных и диффузорных элементах, поля температур в проточной части канала, представленные на рис. 3-11.

На рис. 3-8 представлены графики формирования профиля осевой, радиальной и окружной компоненты скорости по длине вращающегося диффузора и конфузора с прямолинейной, криволинейной и криволинейной оребренной проточной частью при изменении температуры в радиальном сечении (ю = 50с-1). Как следует из расчетных данных, во всем диапазоне исследуемых режимов течения наблюдается параболический закон распределения скоростей по сечению трубы (у стенки скорость убывает, а в центре сечения возрастает), что, в целом, справедливо для ламинарного режима движения жидкости в трубах, в том числе и конфузорно-диффузорных. В случае каналов «конфузор-диффузор» с криволинейной теплообменной поверхностью рост скорости на оси (в сравнении с трубой, выполненной в виде чередующихся усеченных прямых конусов) обусловлен общим снижением гидравлического сопротивления за

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

О 0.2 0.4 0.6 0.8 1

гЩг)

0.2 0.4 0.6 0.8

Г/Я(2)

а)

б)

в)

Рис. 3. Формирование профиля: а) осевой; б) радиальной; в) окружной компоненты скорости при изменении температуры в радиальном сечении

по длине криволинейного диффузора (ю = 50с_1)

а) б) в)

Рис. 4. Формирование профиля: а) осевой; б) радиальной; в) окружной компоненты скорости при изменении температуры в радиальном сечении по длине прямолинейного диффузора (ю = 50с_1)

а)

б)

в)

Рис. 5. Формирование профиля: а) осевой; б) радиальной; в) окружной компоненты скорости при изменении температуры в радиальном сечении

по длине криволинейного оребренного диффузора (ю = 50с_1)

а)

б)

в)

Рис. 6. Формирование профиля: а) осевой; б) радиальной; в) окружной компоненты скорости при изменении температуры в радиальном сечении по длине криволинейного конфузора (ю = 50с_1)

гЩг) т/Я(г) г/И(2)

а) б) в)

Рис. 7. Формирование профиля: а) осевой; б) радиальной; в) окружной компоненты скорости при изменении температуры в радиальном сечении по длине прямолинейного конфузора (ю = 50с_1)

а) б) в)

Рис. 8. Формирование профиля: а) осевой; б) радиальной; в) окружной компоненты скорости при изменении температуры в радиальном сечении

по длине криволинейного оребренного конфузора (ю = 50с-)

Рис. 9. Формирование профиля температуры жидкости по длине криволинейного канала (ю

Рис. 10. Формирование профиля температуры жидкости по длине прямолинейного канала (ю

г

Рис. 11. Формирование профиля температуры жидкости по длине криволинейного оребренного канала (ю = 50с-)

Более сложный характер распределения поля скоростей в конфузорно-диффузорных каналах наблюдается в радиальном направлении. Так, в диффузоре под влиянием центробежных сил при движении жидкости к выходному сечению имеет место рост численных значений радиальной скорости в канале. В конфузоре по мере возрастания осевой скорости в направлении к выходному сечению и стесненности потока наблюдается увеличение градиента скорости движения жидкости вдоль конически-сходящегося элемента, что вызывает рост сил внутреннего трения. При этом значения компонентов радиальной скорости (по отношению к диффузору) уменьшаются, а вектор скорости меняет направление на противоположный, что согласуется с данными работы [8]. Эта закономерность, наряду с положительным влиянием на теплообмен центробежного поля, дополнительно интенсифицирует конвективный перенос тепла, что определяется оригинальным последовательно повторяющимся расположением системы «диффузор-конфузор» [1, 2], определяющих общую геометрию канала. Следует отметить, что наибольшие значения величины радиальной скорости наблюдаются в конфузорно-диффузорном канале с оребренной проточной частью, так как в этом случае оребрение выполняет роль радиально расположенных лопаток, снижающих эффект проскальзывания среды относительно стенок трубы, ранее подобное явление было описано в работе [8].

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

На рис. 9-11 представлены результаты численного решения задачи сопряженного теплообмена, позволяющие установить закономерности изменения температуры жидкости по длине исследуемых каналов. Анализ показывает, что за счет профилирования стенок конфузорно-диффузорных элементов растет интенсивность нагрева среды в проточной части канала, причем по отношению к каналу «конфузор-диффузор» с прямыми стенками в среднем в 1,25 раза, а за счет оребрения проточной части практически в 2 раза. Как показали последующие исследования, длина проточной части таких каналов может быть сокращена почти в два раза.

СПИСОК ЛИТЕРАТУРЫ

1. Золотоносов А.Я., Золотоносов Я.Д. Теплообмен в аппарате типа «труба в трубе» с вращающейся теплообменной поверхностью «конфузор-диффузор» и оребренной проточной частью // Известия КазГАСУ, 2010, № 1 (13). - С. 200-211.

2. Золотоносов А.Я. Построение профиля стенок криволинейных теплообменных элементов в трубах «конфузор-диффузор» // Известия КазГАСУ, 2010, № 2 (14). - С. 168-175.

3. Золотоносов Я.Д., Золотоносов А.Я., Белавина Т.В., Хайруллин М.Р. Математические модели ламинарного течения вязкой жидкости во вращающихся каналах типа «конфузор-диффузор» // Сборник научных трудов КазГАСУ. - Казань, 2010. - С. 221-228.

4. Сегерлинд Д. Применение метода конечных элементов. Пер. с англ. - М.: Изд-во Мир, 1979. - 392 с.

5. Формалев В.Ф., Ревизников Д.Л. Численные методы. - М.: Изд-во ФИЗМАТЛИТ, 2004.- 400 с.

6. Багоутдинова А.Г., Золотоносов Я.Д. Численное решение уравнений движения методом конечных элементов в условиях ламинарного течения вязкой жидкости в конвергентном канале // Известия вузов. Проблемы энергетики, 2007, № 3-4. - Казань: Изд-во КГЭУ. - С. 23-27.

7. Темам Р. Уравнения Навье-Стокса. Теория и численный анализ. - М.: Изд-во Мир, 1981. - 408 с.

8. Пантелеева Л.Р. Теплообмен при ламинарном течении вязкой жидкости в теплообменных устройствах типа «труба в трубе» с вращающейся поверхностью «конфузор-диффузор». Дис... канд. техн. наук. - Казань, 2005. - 116 с.

9. Гук Г.А. Процессы и аппараты молочной промышленности. - Промтехиздат, 1995. - 210 с.

REFERENCES

1. Zolotonosov A.Ya., Zolotonosov Ya.D. Heat Interchange in the type device «pipe in pipe» with rotating heat-exchange surface «confusor-diffusor» and ribbed flowing part // Izvestiya KazGASU, 2010, № 1 (13). - P. 200-211.

2. Zolotonosov A.Ya. A construction of profile of walls of curvilinear heat exchange elements of «confusor-diffusor» type pipes // Izvestiya KazGASU, 2010, № 2 (14). - P. 168-175.

3. Zolotonosov A.Ya., Zolotonosov Ya.D., Belavina T.V., Khairullin M.R. Mathematical models of the laminar viscous fluid fllow in the «pipe in pipe» type revolving channels // Sbornik nauchnych trudov KazGASU. - Kazan, 2010. - P. 221-228.

4. Segerlind D. Application of the finite-element method. Eng. tran. - М.: Izd-vo Mir, 1979. - 392 p.

5. Formalev V.F., Reviznikov D.L. Numerical methods. - М.: Izd-vo FIZMATLIT, 2004. - 400 p.

6. Bagoutdinova AG, Zolotonosov JD The numerical solution of equations of motion of the finite element method in a laminar flow of viscous fluid in the convergent channel / / Proceedings of the universities. Energy problems, 2007, № 3-4. - Kazan: Izd. KGEU, - Р. 23-27.

7. Temam R. Navier-Stokes equations. Theory and numerical analysis. - М.: Publishers Мир, 1981. - 408 p.

8. Panteleyev L.R. Heat transfer in laminar flow of viscous fluid in the heat transfer devices such as the «tube in tube» with a rotating surface, «confuser-cone». Dis ...candidate. tech. Science. -Kazan, 2005. - 116 p.

9. Hooke G.A. Processes and devices dairy industry. Promtehizdat, 1995. - 210 p.