Двухфазные парокапельные течения часто встречаются в теплотехнике и энергетике, где водяной пар используется в качестве рабочего тела. Наиболее широко такие течения представлены в элементах технологического оборудования АЭС, так как отличительной чертой АЭС из-за конструктивных особенностей атомных реакторов является использование слабоперегретого, насыщенного и влажного пара, что приводит к интенсивному влагообразованию (см., например, [26]). Процессы, происходящие во влажном паре, играют чрезвычайно важную роль в работе паровых турбин и другого технологического оборудования, и поэтому без адекватного моделирования парокапельных течений невозможно четкое и детальное представление о процессах, происходящих в этих устройствах. С другой стороны, эффекты, связанные с наличием в паре дисперсной влаги, довольно многочисленны и порой чрезвычайно сложны: фазовые переходы, нук-леация (образование в переохлажденном паре ядер конденсации), скольжение (разность скоростей пара и капель), дробление и коагуляция (слияние) капель, осаждение капель на твердые поверхности с образованием пленки, последующий отрыв этой пленки от задних кромок турбинных лопаток в виде крупных капель и т. д. На многие из перечисленных эффектов существенно влияет также турбулентность несущего потока.
Математические модели, аналогичные используемым для исследования парокапельных потоков, находят свое применение также и в других областях: моделирование движения газа с твердыми частицами или жидкими каплями во всевозможных влагои пылеуловителях, соплах реактивных двигателей, топочных устройствах и промышленных химических установках.
Цель данной работы состояла в исследовании некоторых моделей двухфазных сред, а также в создании универсальной численной методики, пригодной для расчета как можно более широкого круга моделей и задач, связанных с двухфазными течениями в геометрически сложных областях.
Применяемые математические модели парокапельных течений весьма разнообразны и различаются учетом или неучетом тех или иных физических эффектов из числа приведенных выше, так как в зависимости от конкретных приложений и конкретных параметров потока те или иные эффекты становятся определяющими или, наоборот, мало влияющими. Кроме того, модели различаются по способу описания дисперсной фазы, по введению тех или иных упрощающих предположений.
Наиболее общее описание состояния дисперсной фазы получается, если ввести функцию распределения капель по радиусам /(х, Г, I) (или по массам Рт (х, V, т, Т, ?)), где х, V, г, т и Т есть радиус-вектор, скорость, радиус, масса и температура капли,? — время. (У функции / или Рт могут быть также и другие аргументы, например, скорость вращения капли, параметры, характеризующие ее форму, и т. п.) Модель многокомпонентой среды в этом случае представляет собой комбинацию уравнений движения несущей фазы (пара) и кинетического уравнения для функции распределения капель: д р + сИу (ръ) = Др, дръ? дГ сНу (/т <8> = — grad Р + Я р W) д т ре + w сНУ.
В.1).
Здесь в правых частях стоят выражения, описывающие взаимодействия между различными каплями и между каплями и несущей фазой. Через р, и? в (В.1) обозначены плотность, скорость и удельная внутренняя энергия несущей фазы, Р — давление. Заметим, что (В.1) уже несколько ограничивает класс возможных моделей, так как в качестве уравнений для несущей фазы здесь фигурируют уравнения Эйлера. При необходимости учета более сложных эффектов в несущей фазе (вязкость, турбулентность) эти уравнения должны быть заменены на уравнения Навье — Стокса или на уравнения той или иной модели турбулентности.
Однако модель (В.1) непригодна для практических вычислений из-за слишком большого числа независимых переменных и может рассматриваться лишь как базовая для вывода более простых моделей. Один из классов таких моделей получается посредством аппроксимации функции распределения Рто (х, V, т, Т, I) на лагранжевой сетке. Хотя, на первый взгляд, это не уменьшает трудностей при реализации модели, есть надежда, что при этом удастся отследить основные, наиболее характерные группы капель, имеющиеся в потоке, не прибегая к чрезмерно большому количеству лагранжевых ячеек.
Наиболее последовательное применение этого метода, когда лагранжева сетка вводится во всем пространстве переменных (х, у, т, Т), приводит к моделям типа «больших частиц», аналоги которых широко используются для моделирования движения заряженных частиц в электромагнитном поле. Суть метода состоит в том, что на неподвижной разностной сетке, аппроксимирующей пространственную расчетную область, решаются лишь уравнения движения несущей фазы, а дисперсная фаза представляется ограниченным набором групп капель, выбранных так, что капли из одной группы имеют относительно близкие координаты и другие параметры. Это позволяет отследить движение всей группы капель («большой частицы»), проинтегрировав уравнения движения для одной капли, имеющей характерные параметры для этой группы. При этом, разумеется, необходимо, с одной стороны, интерполировать параметры несущего потока с сеточных узлов в точку, где находится «пробная» капля, а с другой стороны, в каждом узле сетки вычислять усредненные параметры дисперсной фазы, чтобы учесть обратное влияние капель на пар. Подобную методику применяли авторы [59].
Другая модель в рамках этого подхода получится, если ввести лагранжеву сетку только по переменным V, т и Т. Получающаяся при этом модель оперирует с несколькими группами капель, каждая из которых занимает, вообще говоря, всю расчетную область. Таким образом, при реализации такой модели нет необходимости переинтерполировать величины из эйлеровых узлов в лагранжевы и обратно, но для каждой группы капель необходимо решать уравнения баланса, похожие на уравнения несущей фазы. Модели этой группы использовались, например, в [33]. Результаты этих расчетов приведены также в [18, 29]. Пример такой модели и расчета по ней (изложенный также в [8]) приведены в настоящей работе, см. 3.1.
Широкий класс моделей получается из (В.1) в предположении, что капли малы и, следовательно, их скорость равна скорости пара (о дно скоростные модели). Часто при этом также предполагают, что температура капель равна температуре пара или некоторому другому, но также точно определенному, значению, оправдывая это предположение также малостью капель. Тогда среди аргументов функции распределения, кроме пространственных координат х и времени остается лишь один «лишний» параметр — масса капли т (или радиус ?¦). Область применимости таких моделей — мелкодисперсные парокапельные среды, как правило, возникающие при спонтанной конденсации в переохлажденном паре, см. [29, 30, 27, 35, 36, 18].
Односкоростные (и однотемпературные) модели могут быть получены из (В.1) с помощью подстановки туда функции распределения Рт специального вида:
Рт (х, у, т, Т,*) = Р^(х, т, Т, г) • 6(у или.
Рт (х, V, т, Г, *) = р-(х, т, I) ¦ 6{V — • 6(Т — Г0), где Го — заданная температура капель. Штрихи у Р’т после этого опускаются.
Односкоростное приближение позволяет, прежде всего, упростить систему уравнений за счет того, что вместо уравнений несущей фазы можно записать уравнения баланса массы, импульса и энергии для всей многокомпонентной среды. При этом в правых частях этих уравнений взаимно сокращаются члены, описывающие межфазовый обмен массой, импульсом и энергией.
Среди односкоростных и однотемпературных моделей можно выделить метод моментов (изложенный и широко используемый, например, в [29]), основанный на замене функции распределения /(х, г, ?) ее моментами [35, 36, 34] оо.
Пп (М) = I Дх, М) гп^г, п = 0,1,2,. о.
При этом величина оказывается связанной простой зависимостью с содержанием жидкой фазы (средний куб радиуса пропорционален среднему объему, а следовательно, и средней массе капли). Поэтому при некоторых дополнительных предположениях (что скорость роста капли г не зависит от ее радиуса) система, содержащая уравнения баланса для величин Г20, и 03, становится замкнутой (естественно, вместе с уравнениями баланса массы, импульса и энергии всей многокомпонентной смеси).
Кроме того, к функции распределения Рт (х,(или /(х, г,?)) в полной мере применимы и те подходы, которые излагались выше для функции распределенияРт (х, у, т, Т, ?) (соответственно /(х, v, г, Г, ?)). При этом получаются либо методы типа «больших частиц», но только односкоростные (см. [59]), либо модели, основанные на делении всего спектра размеров капель на фракции (методы лагранжевых фракций, в качестве примера использования можно привести работу [44]).
И, наконец, в односкоростном приближении возможна аппроксимация функции распределения на эйлеровой сетке (поскольку дополнительная независимая переменная всего одна, получающаяся задача вполне по силам существующей вычислительной технике). При этом спектр радиусов капель также делится на фракции, но теперь границы этих фракций остаются неподвижными, а капли могут переходить из фракции во фракцию по мере их испарения или конденсации (метод эйлеровых фракций [13, 7, 3]).
Выбор того или иного численного метода для решения системы дифференциальных уравнений какой-либо из моделей в значительной мере произволен, так как на этом этапе моделирования уже не делается никаких новых предположений о преобладании тех или иных физических процессов, и, следовательно, в какой-то мере результаты моделирования должны быть независимы от используемого численного метода. Однако численные методы могут различаться по таким параметрам, как точность, надежность, эффективность на данном классе задач, быстродействие, удобство в реализации, универсальность (возможность легкой модификации при изменении модели). Вместе с тем спектр используемых численных методик отличается большим многообразием. Так, например, авторы [59] для решения уравнений несущей фазы использовали четырехшаго-вую схему интегрирования по времени на основе метода Рунге — Кутта, а для интегрирования уравнений движения капли ими был разработан оригинальный подход, включающий в себя отслеживание траектории капли назад по времени.
Для двухскоростных моделей, основанных на лагранжевой аппроксимации по V, т и Г, в [10] предложен смешанный эйлерово-лагранжев (СЭЛ) подход (см. [25]), предполагающий использование различных сеток для интегрирования уравнений несущей фазы и частиц (т. е. с выделением особенности — сепаратрисы, разделяющей области чистого пара и двухфазной среды). Аналогичные уравнения решались в [33] (результаты см. также в [18, 29]) на неподвижной сетке (без выделения особенностей) методом Годунова.
Следует заметить, что метод Годунова [15, 16] широко используется в работах [33, 29] как для многоскоростных, так и для односкоростных моделей на основе метода моментов. При этом проводилось исследование влияния частиц на распад газодинамического разрыва, и формулы расчета распада разрыва соответствующим образом модифицировались. Учет источниковых членов, не содержащих пространственных производных (т. е. межфазного взаимодействия), проводился явным либо неявным образом.
Для расчета стационарных режимов течения в [29] применялись также метод характеристик и (для одномерного случая) метод Рунге — Кутта. Однако эти методы имеют ограниченную применимость, так как позволяют рассчитывать только сверхзвуковые гладкие течения, в то время как в основном интерес представляют трансзвуковые течения, в которых при некоторых режимах конденсации возможны разрывы и нестационарности.
Численная методика, построенная в настоящей работе, основана на раздельном учете гидродинамических процессов и межфазового взаимодействия (метод расщепления, или метод суммарной аппроксимации, см. об этом [31, 32]). Такой подход определяется принципиально разной природой этих процессов, различными характерными временами и различными подходами к численному моделированию этих процессов.
На этапе учета гидродинамических процессов определяющим фактором является необходимость решения гиперболической системы уравнений в области, зачастую имеющей сложную геометрическую форму. Кроме того, некоторые особенности решаемых задач делают крайне желательной возможность адаптации сетки, сгущения ее в том или ином месте области. Это делает привлекательным использование неструктурированных разностных сеток. Это направление в настоящее время бурно развивается, поскольку неструктурированные сетки обладают рядом преимуществ перед структурированными, сводящихся в основном к упрощению алгоритмов автоматизированного построения и адаптации таких сеток, а также к большим возможностям дискретизации областей сложной формы. Кроме того, в целях повышения точности расчета и лучшего разрешения особенностей течения целесообразно применять разностные схемы повышенного порядка аппроксимации.
Как известно, при численном решении гиперболических систем уравнений разностные схемы первого порядка имеют большую схемную вязкость, что выражается в сильном (зачастую неприемлемом) «размазывании» особенностей решения, в то время как разностные схемы второго порядка, полученные на основе аппроксимации центральными разностями, немонотонны (что, как правило, тоже неприемлемо). Проблема уменьшения схемной диффузии при сохранении монотонности схемы по этой причине издавна стояла перед вычислителями, и были предложены разные подходы к ее решению (см. [37, 17, 42, 22, 45, 57, 52, 14, 58, 43, 41] и многие другие работы). В настоящее время чаще всего используются разностные схемы класса TVD (total variation diminishing — схемы с невозрастанием полной вариации, см., например, [46]). Однако исследуются также и расширенные классы квазимонотонных схем: TVB (с ограниченной полной вариацией, см. [54]) и ENO (без существенных осцилляций, см. [55, 56]). Применение разностных схем класса TVD выглядит более предпочтительно, так как остальные схемы не устраняют совсем нефизические осцилляции решений, а лишь делают их сколь угодно малыми (не превышающими заранее заданного порога).
Наиболее распространенным способом построения TVD схем (многие другие способы ему эквивалентны в том смысле, что в итоге дают те же самые схемы) является введение в схему первого порядка специально ограниченных антидиффузионных потоков [57, 14], и именно это направление было выбрано для построения рассматриваемой здесь численной методики. Вопросам обобщения одномерной схемы с ограничителями антидиффузии на случай неструктурированной сетки в многомерной области посвящены работы [12, 9, 11, 53] и др. Построению методики повышенного порядка на неструктурированных сетках посвящены также работы [49, 48, 60], но в них в основном развиваются конечноэлементные подходы.
На этапе учета межфазного взаимодействия на первый план выходит проблема возможной жесткости решаемых систем уравнений, то есть слишком малых времен релаксации процессов взаимодействия фаз по сравнению с используемыми временными шагами. Хотя, например, в [29] отмечалось, что, как правило, межфазное взаимодействие не приводило к необходимости уменьшения шага по времени по сравнению с выбранным исходя из гидродинамического условия Куранта, тем не менее были зафиксированы случаи, когда жесткость системы уравнений кинетики фазового перехода все же проявлялась. Известные способы решения жестких систем сводятся либо к дроблению шага по времени для учета межфазного взаимодействия по сравнению с гидродинамическим, либо к применению асимптотически устойчивых методов решения обыкновенных дифференциальных уравнений (все такие методы являются неявными). Однако первый из этих способов непригоден при очень малых временах релаксации, когда для устойчивого счета пришлось бы раздробить шаг по времени на слишком большое число подшагов, а второй сложен в реализации из-за сильно нелинейной и сложно вычисляемой правой части системы уравнений. В настоящей работе показано, что для определенного класса жестких систем (а именно для таких, которые описывают одиночный релаксационный процесс) применение таких подходов, хотя в принципе и возможно, но излишне: в вычислительном плане оказывается проще при обнаружении жесткости системы найти равновесное состояние.
Диссертация состоит из трех глав. Глава 1 посвящена последовательному и строгому изложению применяемых численных алгоритмов как для учета конвективных процессов, так и межфазного взаимодействия. Для учета конвективных процессов используется изложенная в 1.1 квазимонотонная разностная схема повышенного порядка точности на неструктурированных сетках. Эта. разностная схема является обобщением подхода [51, 12] по способу аппроксимации разностных выражений и метода крупных частиц [4] по способу вычисления потоков. В 1.1.3 проведено теоретическое исследование монотонности и устойчивости данной методики применительно к линейному уравнению переноса. В 1.1.6 приведен пример решения реальной задачи (не двухфазной) по этой методике с использованием неструктурированных сеток.
Глава 2 посвящена некоторым вопросам алгоритмической реализации приведенной здесь и других подобных методик на неструктурированных сетках.
Глава 3 содержит подробные описания некоторых моделей двухфазной среды и результаты применения этих моделей к конкретным задачам. В 3.1 рассматривается многоскоростная модель для перегретого пара с крупнодисперсной влагой, а в 3.2 — односкоростные методы моментов и эйлеровых фракций. В 3.2.5 приведены результаты расчетов стационарных и нестационарных режимов течения в соплах со спонтанной конденсацией.
Некоторые результаты, вошедшие в диссертацию, докладывались на:
• Второй международной научно-технической конференции «Актуальные проблемы фундаментальных наук», Москва, 1994 г.
• III Минском международном форуме по теплои массообмену, Минск,.
1996 г.
• Третьем международном симпозиуме по экспериментальной и вычислительной аэротермодинамике внутренних течений, Пекин, 1996 г.
• Международном симпозиуме по физике теплопередачи при кипении и конденсации, Москва, 1997 г.
• Семинаре ИММ РАН под руководством академика А. А. Самарского,.
1997 г.
• Совместном семинаре отделов 5, 8 ИММ РАН и кафедры математического моделирования МФТИ под руководством проф. Е. И. Леванова, 1999 г.
Основные результаты.
• Для квазилинейных гиперболических уравнений типа конвективного переноса разработана и изучена квазимонотонная разностная схема повышенного порядка точности на неструктурированных сетках, использующая т. н. комплексное описание сетки. Предложена и используется система обозначений, существенно облегчающая работу с конечноразностными выражениями на таких сетках.
• Предложен универсальный подход к построению численных алгоритмов, реализующих модели двухфазных течений, основанный на раздельном учете межфазных взаимодействий (в т. ч. кинетики фазовых переходов) и газодинамических (конвективных) процессов. Для учета межфазных взаимодействий, носящих, как правило, характер релаксационных процессов, ведущих к установлению равновесного состояния, предложен простой алгоритм, позволяющий без использования неявных схем учитывать как-быстрые, так и медленные процессы.
• Создана континуальная математическая модель парокапельного течения в быстродействующей редукционно-охладительной установке (БРОУ), основанная на разбиении капель на группы по их начальному радиусу. Построен алгоритм ее численной реализации и проведен вычислительный эксперимент, показавший практическую применимость этой модели.
• На основе представления функции распределения частиц в парокапель-ном потоке эйлеровыми фракциями разработана методика моделирования трансзвуковых течений влажного спонтанно конденсирующегося пара, позволяющая учесть гомогенные и гетерогенные фазовые переходы и зависимость температуры капли от радиуса. Выполнены расчеты течений в экспериментальных сопловых аппаратах. Показано, что на рассмотренном классе течений методика дает результаты, сопоставимые с традиционно применяющейся методикой, основанной на решении уравнений для моментов.
