Приводятся результаты имитационного численного моделирования течения и сопряженного тепломассопереноса в трехмерной модели испарителя контурной тепловой трубы спутника TacSat-4. Установлено, в частности, что отводящие пар аксиальные канавки испарителя работают в существенно разных условиях - расходы пара отличаются в несколько раз. Математическая модель основана на общих уравнениях баланса массы, импульса и энергии в потоке рабочей жидкости и пара, в пористом фитиле и в твердых элементах конструкции. Для решения уравнения Пуассона использован эффективный многосеточный метод на основе оригинального алгоритма агломерации ячеек. Ключевые слова: численное моделирование, сопряженный тепломассоперенос, испаритель контурной тепловой трубы, многосеточный метод, неструктурированные сетки.
Численное моделирование сопряженного тепломассопереноса в испарителе контурной тепловой трубы Внутри расчетной области на границе между твердой стенкой и паром, а также на границе между твердой стенкой и пористым телом ставится условие сопряженного теплообмена – равенство тепловых потоков и температур по обе стороны от границы, а для уравнения движения – условие прилипания. На границе фитиль-пар температура обеих фаз полагается одинаковой и равной температуре насыщенного пара при давлении, равном давлению в паровой фазе. Связь между давлением и температурой на линии насыщения определяется по известной формуле Антуана. Нормальная компонента скорости на границе фитиль-пар вычисляется исходя из скорости испарения, определяемой по величине подводимого к границе теплового потока. Для пара касательная составляющая скорости полагается равной нулю. Для жидкости в пористом теле условия на касательную составляющую скорости не ставятся, а условия для нормальной составляющей выражаются через градиент давления с помощью закона фильтрации Дарси. 3. Программная реализация. Алгоритм агломерации
Для проведения расчетов используется программный код Flag-S [
Как показывает практика гидродинамических расчетов, наиболее ресурсоемкой частью
алгоритма SIMPLEC является решение уравнения Пуассона для поправки давления. Уравнение
Пуассона описывает также теплопроводность в твердых элементах конструкции и
фильтрационное движение рабочей жидкости в пористом фитиле. Таким образом, использование в
линейном солвере специализированного алгоритма, обеспечивающего эффективное решение
уравнения Пуассона, способно в значительной степени снизить затрачиваемое на решение задачи
время и, соответственно, заметно повысить эффективность использования имеющихся
вычислительных ресурсов. При решении линейных систем большой размерности наибольшую
эффективность демонстрируют многосеточные методы, т.к. их использование позволяет добиться
линейного роста времени решения задачи, при увеличении ее размерности. Основная трудность
применения многосеточных методов на неструктурированных сетках заключается в построении
сеток для грубых уровней. Одним из наиболее перспективных считается агломерационный
метод [10], в рамках которого ячейки грубых расчетных сеток получаются путем агломерации
(объединения) ячеек исходной сетки. Вопрос о наиболее эффективном способе агломерации
ячеек до сих пор остается открытым и является предметом ряда современных исследований,
например [
В рамках настоящей работы в код Flag-S был внедрен оригинальный агломерационный
многосеточный метод. Способ агломерации базируется на предложенном в работе [
1) Построение базиса. Ячейка из очереди объединяется с двумя (для двумерного случая) или тремя (для трехмерного случая) соседними ячейками, у которых имеется общий узел, и расстояние до которых не превышает расстояние до ближайшей ячейки более чем в 4 раза. Если имеется несколько вариантов объединения, то выбирается тот, который содержит большее количество уже находящихся в очереди ячеек (если вариантов, по-прежнему, несколько, то берется первый из них). Если подходящего набора ячеек не находится и расстояние хотя бы до одной из соседних ячеек превышает расстояние до ближайшей ячейки более чем в 4 раза, то первый шаг алгоритма повторяется, с попыткой построить базис, содержащий на одну ячейку меньше. Если в итоге удается сформировать базис, то осуществляется переход к шагу 2; иначе – к рассмотрению следующей ячейки из очереди. 2) Добор ячеек. Просматриваются ячейки, соседствующие с полученным агломератом. Если среди них есть ячейка, соседствующая не менее чем с двумя ячейками из агломерата, и добавление этой ячейки уменьшает значение параметра Ar (aspect ratio, определяемого как отношение линейных размеров, характеризующих площадь поверхности и объем агломерата), то эта ячейка добавляется в агломерат; иначе осуществляется переход к следующему шагу. Если подходящих ячеек несколько, то выбирается та из них, добавление которой приводит к меньшему значению параметра Ar. Если число ячеек, вошедших в агломерат, меньше 4 (двумерный случай) или 8 (трехмерный случай), то шаг повторяется.
3) Пополнение очереди. Ячейки, которые соседствуют с полученным агломератом и еще не были объединены с другими ячейками, добавляются в очередь (если они еще не были в нее добавлены). Затем осуществляется переход к рассмотрению следующей ячейки из очереди (шаг 1).
4) «Залатка дыр». Все необъединенные ячейки маркируются. Выбранная необъединенная ячейка объединяется с одним из соседних агломератов. При этом ставится запрет на добавление необъединенной ячейки к агломерату в случае, когда все соседи этой ячейки, входящие в рассматриваемый агломерат, были маркированы (т.е. тоже были необъединенными ячейками). Среди разрешенных для объединения агломератов выбирается тот, для которого добавление этой ячейки сильнее всего уменьшит (или слабее всего увеличит) значение параметра Ar. Затем осуществляется переход к следующей необъединенной ячейке.
Число сеточных уровней подбирается автоматически – построение новых сеток
продолжается пока в каждом расчетном блоке не останется по 1 ячейке. Разностный оператор на грубых
сетках формируется путем суммирования уравнений с масштабированием коэффициентов [
Эффективность разработанного многосеточного метода была проверена на ряде тестовых задач. Первая из них – это задача теплопроводности в кубической области с постоянной теплопроводностью. Рассматривался случай, когда почти все границы области – адиабатические, за исключением двух небольших, прилегающих к вершинам участков, на одном из которых фиксировалось значение температуры, а на втором – теплового потока. Безразмерная постановка задачи иллюстрируется на рисунке 1а. Сопоставительные расчеты проводились с использованием кода Flag-S и кода ANSYS Fluent 14.0 на равномерных (изотропных) гексаэдральных и тетраэдральных расчетных сетках разного пространственного разрешения, а также на полиэдральных сетках, полученных из тетраэдральных с помощью конвертации средствами пакета Fluent. При расчетах по коду Flag-S в качестве линейного солвера применялись: разработанный многосеточный метод, метод GMRES (с предобуславливателем SGS) и алгебраический многосеточный метод, реализованный в соответствующей процедуре открытой библиотеки HYPRE. Из опций, доступных в коде Fluent, выбран наиболее эффективный для данной задачи (согласно результатам предварительных тестов) линейный солвер – метод BCGSTAB с многосеточным предобуславливателем, использующим W-цикл. В ходе тестирования оценивалось время, необходимое для решения задачи (задача считалась решенной, когда значение относительного интегрального небаланса теплового потока становится меньше 10-4).
На рисунке 1б представлена зависимость времени решения задачи от числа ячеек расчетной сетки при использовании различных кодов и численных методов. Видно, что при использовании метода GMRES время решения задачи растет пропорционально количеству ячеек в степени 1.5. Вместе с тем, применение любого из указанных выше многосеточных алгоритмов обеспечивает, как и ожидалось, линейный рост времени решения задачи при увеличении размерности расчетной сетки. Разработанный многосеточный метод в несколько раз (до одного порядка) эффективнее, чем библиотечный алгебраический многосеточный метод. На полиэдральной сетке разработанный и реализованный в коде Flag-S многосеточный метод позволил решить данную задачу примерно в 2 раза быстрее, чем при использовании кода Fluent; в случае гексаэдральных и тетраэдральных сеток на решение задачи по двум кодам требовалось примерно одинаковое время.
T=0 а) q=1
Для проверки работоспособности многосеточного метода в более сложных условиях была рассмотрена трехмерная стационарная задача теплопроводности с трансляционной периодичностью, иллюстрируемая на рисунке 2. Расчетная область имеет ступенчатые вырезы, что позволило увеличить число (усложняющих задачу) разрывов в граничных условиях первого рода – на разных изотермических участках границы области задавались различные температуры (использовался набор из трех значений). Другая отличительная особенность данного вычислительного теста заключается в использовании сеток, содержащих сильно вытянутые ячейки – максимальное удлинение (отношение наибольшего размера ячейки к наименьшему) составляло 5000. Использовались две расчетные сетки: гибридная – с тетраэдральными ячейками внутри области и сеточными пограничными слоями из призматических ячеек у изотермических границ, и полученная из нее полиэдральная сетка. Гибридная сетка содержала 1 миллион ячеек, полиэдральная – 300 тысяч ячеек. Фрагмент полиэдральной сетки показан на рисунке 2.
Рис. 2. Схема к постановке второй тестовой задачи и фрагмент полиэдральной расчетной сетки Время решения задачи при использовании разработанного многосеточного метода составило 54 и 30 секунд на гибридной и полиэдральной сетках соответственно. При использовании библиотечного алгебраического многосеточного метода затрачиваемое на решение задачи время превышало указанное, соответственно, в 4.4 и 3.3 раза.
Результаты представленных и других тестов позволяют заключить, что разработанный метод показывает высокую эффективность и обеспечивает устойчивость счета при использовании различных расчетных сеток, включая случаи сильно вытянутых ячеек и быстрых изменений шага сетки по пространству. 5. Задача тепломассопереноса
Рассматривается модель испарителя, близкого по форме к установленному на американском микроспутнике TacSat-4 [8], при этом часть геометрических параметров была оценена по приведенным в [8] фотографиям экспериментального стенда. С учетом симметрии задачи рассматривается половина испарителя. Расчетная область, представленная на рисунке 3, включает в себя корпус испарителя, пористый фитиль и занятую паром область, включающую паровые канавки, сборный пленум и выходной патрубок. Модель испарителя имеет следующие геометрические характеристики: высота корпуса испарителя составляет 3 см, ширина нижнего и верхнего оснований корпуса – 6.7 см и 2.2 см соответственно, внешний и внутренний диаметры фитиля – 2.5 см и 1.6 см, высота и ширина паровых канавок – 0.7 мм и 1.5 мм соответственно, общая длина фитиля – 30.5 см, длина торца фитиля – 1.25 см, длина сборного пленума – 4 мм, длина выходного патрубка – 1.4 см. В представленном ниже расчете на всей внутренней поверхности фитиля задавалось условие постоянной температуры в 255.5 К, соответствующей температуре в компенсационной камере по данным [8], а для течения жидкости данная граница считается «свободной» – на ней задавалось условие постоянного давления. На нижней поверхности полки корпуса испарителя задавалось линейное распределение температуры вдоль испарителя – от 270 К на стороне выхода пара до 272 K на противоположном конце корпуса. Все остальные границы полагались адиабатическими. Уровень давления на выходе из испарителя задавался равным 2.22∙105 Па, в соответствии с полученной в экспериментах [8] температурой пара на выходе из испарителя, составляющей 256.7 К (в предположении, что в этой области пар находится в состоянии насыщения).
При проведении расчетов использовались следующие физические параметры сред: корпус испарителя (алюминий): λs = 200 Вт/(м∙К); фитиль (спеченный никелевый порошок): пористость ε = 0.75, коэффициент эффективной теплопроводности λeff = 5.5 Вт/(м∙К), проницаемость K = 2∙10−14 м2, жидкость в фитиле (жидкий аммиак): ρl = 660 кг/м3, Cpl = 4510 Дж/(кг∙К), μl = 2∙10−4 Па∙с; пары аммиака: ρv = 1.89 кг/м3, Cpv = 2400 Дж/(кг∙К), μv =8.65∙10−6 Па∙с, λv = 2.1∙10-2 Вт/(м∙К); удельная теплота парообразования 1.31∙106 Дж/кг.
Расчетная сетка содержала 1.3 миллиона ячеек: на паровую зону, фитиль и корпус испарителя приходится 200, 600 и 500 тысяч ячеек соответственно. Как видно из рисунка 3, узлы сетки сильно сгущены к границе между фитилем и паровой зоной – для качественного разрешения температурных градиентов и подводимого к межфазной границе теплового потока.
По результатам расчетов величина подводимого к испарителю теплового потока составила 780 Вт, что примерно соответствует максимальной подводимой тепловой мощности, указанной в [8] для реального прототипа. На рисунке 4а представлено типичное распределение температуры в поперечном сечении испарителя. Наибольшие градиенты температуры наблюдаются в фитиле, вблизи нижней паровой канавки. Вблизи других канавок градиенты температуры заметно ниже и уменьшаются с увеличением расстояния от нагреваемой поверхности испарителя до канавки. На рисунке 4б показано распределение векторов скорости на входной границе в том же сечении. Локальные максимумы скорости наблюдаются напротив паровых канавок, при этом основная часть потока приходится на нижнюю часть устройства. По результатам расчетов расход пара в канавках, в долях от общего расхода, составляет: 7%, 211%, 222%, 27% (от верхней канавки к нижней). а) б) 16 .12 a P ,k 8 P 4 0
в) pфверх pфниз pгверх pгниз 0 0.1 0.2
0.3 Z, m.
Рис. 4. Результаты численного моделирования работы испарителя: а) распределение температуры, б) векторы скорости на входной границе, в) падение давления в газе и в фитиле в зависимости от осевой координаты На рисунке 4в показано падение давления в фитиле (при движении жидкости от входной границы к канавкам) и в паровой области (при движении пара по канавкам до пленума) в зависимости от осевой координаты; данные приведены для верхней и нижней канавок. Падение давления вдоль нижней канавки достигает 14кПа. Если к этому добавить примерно 3 кПа падения давления по фитилю и 1 кПа падения давления по контуру конденсатора (последнее – по данным [8]), то общий перепад давления составит весьма большую величину – 18 кПа. В [8] указывается, что размер пор фитиля составляет 1.4 мкм, оценка по формуле Лапласа предельного капиллярного напора по физическим свойствам аммиака для данных тепловых условий дает величину в 25 кПа, что, как и ожидалось, превышает оцененный по расчетам перепад давления по контуру. При повышении подаваемой на конденсатор тепловой нагрузки расходы рабочего тела и, соответственно, перепады давления увеличатся, и лимитирующим максимальную тепловую нагрузку фактором окажется уровень давления в нижней паровой канавке. Таким образом, работа паровых канавок в существенно разных условиях может приводить к заметному ограничению максимальной тепловой мощности, отводимой испарителем. 6. Заключение
Проведено имитационное численное моделирования течения и сопряженного тепломассопереноса в трехмерной модели испарителя контурной тепловой трубы. Математическая модель базируется на применении трехмерных уравнений Навье-Стокса для паровой фазы и уравнении Пуассона для описания фильтрационного движения рабочей жидкости в пористом фитиле и теплопроводности в твердых элементах. Показано, что при подводе тепла к испарителю только со стороны полки паровые канавки работают в существенно разных условиях; это может значительно ограничивать максимальную тепловую мощность, отводимую испарителем.
Для повышения эффективности численного алгоритма разработан оригинальный агломерационный многосеточный алгоритм. Продемонстрирована высокая эффективность разработанного метода.
В ходе дальнейших исследований планируется выполнить моделирование усложненной модели испарителя, включающей вторичный фитиль и подводящую байонетную трубку, с учетом прогрева рабочей жидкости при ее движении из компенсационной камеры к внутренней поверхности фитиля. Литература 10. Mavriplis D.J. Multigrid techniques for unstructured meshes // ICASE Report No 95-27. – 1995. – 61 p. 14. Patel A. Development of an adaptive RANS solver for unstructured hexahedral meshes // PhD
Thesis. – ULB, Brussels, Belgium. – 2003. – 266 p. 15. Пожилов А.А. Численное решение трехмерных задач теплопроводности на полиэдральных сетках с применением многосеточных алгоритмов: Диссертация на соискание ученой степени магистра: 010900 – Прикладные математика и физика / Санкт-Петербургский государственный политехнический университет, 2014. 48 с. Numerical simulation of conjugate heat and mass transfer in a loop heat pipe evaporator
The results of 3D numerical simulation of conjugate heat and mass transfer in a model of the TacSat-4 satellite loop heat pipe evaporator are presented. In particular, it is found that the vapor axial grooves operate under essentially different conditions – the flow rates differ several times. The mathematical model is based on the mass, momentum and energy conservation laws applied to the fluid and vapor flow, the porous wick, and solids. An efficient multigrid method based on an original agglomeration strategy is used to solve the Poisson equation.
Keywords: numerical simulation, conjugate heat and mass transfer, loop heat pipe evaporator, multigrid method, unstructured grids. 10. Mavriplis D.J. Multigrid techniques for unstructured meshes // ICASE Report No 95-27. – 1995. – 61 p.