CC BY
274
СПИСОК ЛИТЕРАТУРЫ
1. Русаков, М.Р. Основные направления совершенствования руднотермической электроплавки медно-никелевого сырья [Текст] / М.Р. Русаков, К.И. Мосиондз, Ю.С. Жуковский // Цветные металлы,- 1998,- № 2,- С." 36-39.
2. Мосиондз, К.И. Теоретические основы повышения производительности электропечей рудной плавки [Текст] / К.И. Мосиондз, М.Р. Русаков, Ю.С. Жуковский [и др.] // Совершенствование технологии, аппаратуры и методов исследования в производстве тяжелых цветных металлов: сб. науч. тр.— СПб: Гипроникель, 1992,- С. 9-15.
3. Русаков, М.Р. О возможности улучшения энерготехнологических показателей работы электропечей при интенсификации плавки [Текст] / М.Р. Русаков, К.И. Мосиондз, Ю.С. Жуковский / / Цветные металлы,- 1995,- № 12,- С. 9-11.
4. Лыков, А.Г. Исследование электрических и тепловых режимов электропечи для выплавки медного штейна [Текст] / А.Г. Лыков, А.Г. Лунин, И.Ф. Овчинников |и др.] // Цветные металлы,— 1993,- №1,- С. 22-26.
5. Leonard, В.Р. Why you should not use 'hybrid', 'power-law' or related exponential schemes for con-vective modeling— there are much better alternatives [Text] / B.P Leonard, J.E. Drummond // International
Journal for Numerical Methods in Fluids.— 1995.— Vol. 20,- P. 421-442.
6. Плетнев, A.A. Математическое моделирование тепловой работы многошлаковой руднотермической электропечи [Текст] / A.A. Плетнев, В.А. Талалов // Advances of Heat Transfer : сб. докл. конф., Санкт-Петербург, 19— 21 сент. 2007 г.— СПб: Изд-во СПбГТУ, 2007,- Т. 2- С. 449-455.
7. Купряков, Ю.П. Шлаки медеплавильного производства и их переработка / Ю.П. Купряков.— М.: Металлургия, 1987,— 200 с.
8. Ванюков, A.B. Шлаки и штейны цветной металлургии / A.B. Ванюков, В.Я Зайцев,— М.: Металлургия, 1969,— 408 с.
9. Плетнев, A.A. Численное моделирование электрического поля и сопротивления ванны многошлаковой руднотермической печи [Текст] / A.A. Плетнев, В.А. Талалов, М. Р. Русаков // Компьютерное моделирование при оптимизации технологических процессов электротермических производств: сб. науч. тр.— СПб: СПбГТИ(ГУ), 2000,- С. 317-323.
10. Нус, Г.С. Теплогенерация и теплообмен в рудно-термических электропечах металлургии тяжелых цветных металлов [Текст] / Г.С. Нус // Электрометаллургия,— 2000,— N° 9,— С. 7—14.
УДК 536.423
А.Ю. Снегирёв, С.С. Сажин, В.А. Талалов
МОДЕЛЬ И АЛГОРИТМ РАСЧЕТА ТЕПЛООБМЕНА И ИСПАРЕНИЯ КАПЕЛЬ ДИСПЕРГИРОВАННОЙ ЖИДКОСТИ
Процессы испарения диспергированных жидкостей играют важную роль при смесеобразовании в двигателях внутреннего сгорания, авиационных двигателях, топках котлов, работающих на жидком топливе, в системах пожаротушения и охлаждения, в технологиях сушки суспензий, порошковой металлургии, при аварийных выбросах сжиженных (в том числе горючих) газов, а также в атмосферных явлениях.
Испарение диспергированных жидкостей протекает, как правило, в условиях турбулентной газокапельной струи. Составной частью математической модели газокапельной струи является модель испарения одиночной капли, в которой при минимальных вычислительных затратах не-
обходимо учесть все основные факторы, определяющие динамику испарения. Минимизация вычислительных затрат, особенно важная с учетом того, что в газокапельной струе приходится одновременно рассматривать 105—106 и более вычислительных «частиц» (групп капель), требует использования упрощенных подходов к описанию процессов переноса внутри капли.
Развитие и современное состояние теории испарения капель подробно изложено в обзоре [1]. Несмотря на многочисленные исследования, направленные на анализ многообразных и сложных физических процессов, сопровождающих испарение капли, практические расчеты часто ограничиваются простейшими версиями моде-
лей, не учитывающими неравномерность поля температуры в капле, циркуляцию жидкости в ней, влияние теплового излучения, эффекты неравновесности и другие факторы. Это в полной мере относится и к современным версиям коммерческого программного обеспечения для гидродинамических расчетов.
В связи с этим возникает необходимость построения и анализа такой модели испарения капли, которая, с одной стороны, учитывала бы все основные физические явления, определяющие динамику испарения, а с другой — оказалась бы достаточно экономичной для использования в качестве компонента вычислительного кода, предназначенного для численного моделирования турбулентных газокапельных струй вихре-разрешающими методами [2].
Цель данной работы — сформулировать математическую модель, отвечающую перечисленным требованиям, выделить ключевые безразмерные критерии и построить устойчивый и экономичный численный алгоритм.
Описание модели
Нагрев и испарение капли. Скорость убыли массы капли определяется диффузией пара от поверхности капли в окружающую среду [1]:
= -п&Ъ*рТ\п{1 + Вм),
(1)
вм =
У (Т)-У 1-У (7*)
(2)
(— температура поверхности капли, У м — массовая доля пара в окружающей среде).
Массовая доля пара у поверхности капли равна
где — парциальное давление пара у поверхности капли, которое вычисляется в предположении межфазного равновесия (давление насыщенного пара); Р() — давление в окружающей среде;
М = (Р^/Ръ )Мшр + {\-Ршр^Р,) МЕ -
молярная масса смеси (М и — молярная масса пара и газа).
Значение модифицированного числа Шервуда 8Ь* учитывает массопередачу в пограничном слое, формирующемся вокруг капли при обтекании потоком окружающего газа. Если окружающая среда неподвижна, то 8Ь* = 2.
Среднеобъемная температура капли Т определяется из уравнения теплового баланса
с1Т _
_УсопV + ^ар + Чгас! >
(3)
где С/ — теплоемкость жидкости, т — масса капли.
В правой части (3) указаны источники тепла, полученного каплей в результате теплообмена с окружающей средой. Конвективный теплообмен моделируется с учетом расходования тепла на нагрев образующегося пара [3, 1]:
1п^ + дт)у (—ТЛ, (4)
Вт
где т, с! — масса и диаметр капли, р — плотность парогазовой смеси, т — коэффициент диффузии пара, 8Ь — число Шервуда, модифицированное (как показано ниже) с учетом испарения, Вм — массовое число Сполдинга;
а — —т
Чсопу
где — коэффициент теплопроводности парогазовой смеси, Те — температура окружающей среды, В'1 — тепловое число Сполдинга для нестационарного режима нагрева (или охлаждения) капли;
Вт= (1+ ВМ )ф-1. (5)
В равенстве (5) использовано обозначение
где = Т — число Льюиса; 8Ь*,
N11* — числа Шервуда и Нуссельта, модифицированные с учетом испарения.
Модификация выполнена в соответствии с рекомендациями работы [3]:
БЬ =2 +
(6)
1Ми = 2 +
ЬО"2 р = (+п \0,71п(1 + Дм). Гм-{ I+ ЯМ) 7. >
гм "м
(1 + ВТГ^11 (7) Вт
1Ми0-2
р
=
'Т "Г
Числа Шервуда и Нуссельта в отсутствие испарения определяются из эмпирических соотношений для сфер (см. например [1]):
811о=2 + 0,5711е1/28с1/3,
К(ио=2 + 0,5711е1/2Рг1/3. (8)
В равенствах (8) число Рейнольдса Яе для капли вычисляется через ее диаметр и разность между скоростями капли и потока.
Расход тепла на испарение жидкости описывается равенством
1 шр=9Кар т, (9)
где АЛ — теплота испарения при температуре
поверхности капли.
Теплообмен излучением моделируется в предположении о непрозрачности капли:
(7?-С), (10)
где а — постоянна Стефана — Больцмана, е — эффективная степень черноты, Тгас, — радиационная температура окружающей среды.
Если окружающий газ участвует в теплообмене излучением, то его радиационная температура определяется равенством
где С — полная интенсивность излучения, приходящего в данную точку пространства со всех направлений (определяется в ходе расчета переноса теплового излучения).
Если окружающий газ прозрачен, то в качестве Тгас1 следует использовать температуру поверхностей на границах расчетной области (например, температуру стен). В данной модели предполагается, что падающий на каплю радиационный тепловой поток поглощается на ее поверхности. Предполагается также, что численное значение эффективной степени черноты е приближенно учитывает как отражение излучения на поверхности капли, так и пропускание излучения через ее объем. Данный подход является приближенным и не учитывает поглощения излучения в объеме капли, вызванного ее частичной прозрачностью (описание моделей, в которых эти допущения не используются, приведено в обзоре [1]). Несмотря на то, что в настоящей работе вкладом радиационного теплообмена мы пренебрегали, соответствующее слагаемое сохранено, так как оно используется в последующих расчетах с помощью данной модели. В частности, поглощение теплового излучения,
испускаемого пламенем, учитывается при моделировании нагрева и испарения капель распыленной огнетушащей жидкости [2].
Для удобства численных расчетов два первых слагаемых в правой части уравнения (3) преобразуются к следующему виду:
Т -Т
4со»у + 4уар=-С1т--~>
^СОПУ
где Т№Ь — равновесная температура поверхности капли, которая достигается в случае равенства конвективного теплоподвода и теплопотерь на испарение и нагрев пара в отсутствие теплового излучения (температура мокрого термометра). Величина хсот представляет характерный масштаб времени релаксации температуры поверхности капли к ее равновесному значению.
В соответствии с выражениями (11) и (4) конвективный масштаб времени %сот определяется равенством
т =_с1т(Т* ь)_
сот Лп(1 + ВТ) , , Л '
"т
Отметим, что в отсутствие испарения т = 0,
и в качестве Т№Ь следует полагать температуру окружающего газа Те.
Температура мокрого термометра определяется из условия теплового равновесия:
ЯсопЛТ) + ^ар{Т) = 0, (13)
которое дает трансцендентное уравнение относительно Т (температура мокрого термометра ТкЬ является решением этого уравнения).
Радиационное слагаемое также преобразуется к линеаризованному виду:
Т -Т
Ч Г/7// -— >
где
Хма=_^_. (14)
С учетом выражений (11) и (14) для средне-объемной температуры капли будем иметь уравнение
Т.-Т^ Т.-Т.
гас!
йТ Л
которое эквивалентно уравнению вида
(15)
йТ Л
т-т
«7
где
Т _
Тууь/^сопу /^гаё ■
-ед
1
1 1 -+ —
С/Р/'
дТ_
1 д 2. дТ -г—г Х!—,
г дг дг
(17)
где р, — плотность и коэффициент теплопроводности жидкости.
Поле температуры в капле удовлетворяет следующим граничным условиям:
1 дг
(16)
Таким образом, рассматриваемая модель испарения капли состоит из уравнения (1) для массы и уравнения (16) для среднеобъемной температуры капли при соответствующих начальных условиях. В уравнении (16) необходимо определить связь между среднеобъемной температурой Т и температурой поверхности капли Т5.
Поле температуры внутри капли. К настоящему времени используется два подхода к моделированию испарения капель. В первом случае полагают, что поле температуры внутри капли близко к равномерному и потому Т = Т5 .Данный подход называют моделью бесконечной теплопроводности, или моделью быстрого перемешивания. Более универсален подход, учитывающий градиент температуры внутри капли. Модели такого типа называют моделями конечной теплопроводности. Наиболее точный метод учета неравномерности поля температуры внутри капли заключается в решении уравнения переноса энергии совместно с уравнением движения (движение жидкости внутри капли развивается благодаря касательным напряжениям на ее поверхности, которые возникают при движении капли относительно окружающего газа). Однако такой подход представляется недопустимо затратным при моделировании очень больших коллективов капель, как это требуется в расчетах турбулентных газокапельных струй. Упрощенный метод учета циркуляции жидкости внутри капли — использование модели эффективной теплопроводности.
Рассмотрим сначала поле температуры внутри капли в отсутствие циркуляции жидкости в капле. Кроме того, будем полагать, что поле температуры симметрично относительно центра капли:
-X —
1 дг
г=с1/ 2
= 0,
г=0
^сот ^ Ууар ^гай
пё1 '
(18)
Интегрирование уравнения (17) по объему капли приводит к уравнению баланса тепла (3), где
Чшр^Кср-^Т ~Т\г=ф
т,
а среднеобъемная температура Т определяется равенством
Т =
1
ф
п^/б
|Г(г)4
Л.
Отметим, что в работах [1,4] построена и апробирована модель на основе аналитического решения уравнения (17). Преимущество такой модели заключается в ее применимости на начальной стадии нагрева или охлаждения капли. Однако из-за значительного объема необходимых вычислений применение указанной модели представляется обоснованным в расчетах испарения лишь небольшого количества капель. В то же время на практике требуется расчет газокапельных струй с числом вычислительных частиц порядка 105—106, что может привести к неприемлемо большим вычислительным затратам в случае использования такой модели. В связи с этим в данной работе (вслед за статьями [5—7] и др.) использован приближенный метод решения уравнения (17), в котором предполагается параболический профиль температуры внутри капли:
Т{г) = с0-с2г .
(19)
Равенство (19) предполагает установление квазистационарного поля температуры в капле, которое достигается по истечении времени, превышающего масштаб времени релаксации температуры внутри капли.
Выполнение граничного условия (18) на поверхности капли означает, что с2 (У/Я3 )>
ГДе Я = йСОт + 4шр + 4гас! ■
Подстановка профиля (19) в выражение для среднеобъемной температуры дает выражение
Т = С0-С2Зс12/20. Далее, для температуры поверхности будем
иметь:
10 П(ТУ
5 - 3
т=-т--т.
'2 2 5
(20)
Из выражений (11) и (14) следует, что
Т—Т Т—Т
? аф е г
д = —С[т—--—¡т-
^сот ^гас!
гас!
= —с1т
—
-— .(21)
^ед
Подставляя выражение (21) в первое равенство (20), получим:
Т =
1+р
= т+-
р
1+р
(( -Т), (22)
Р =
15 х
ее!
Т; =
Ш1.
Х,/с,р! '
(23)
(24)
йТ
сН
т-т„
ее!
-ед
а+р г
Отметим, что при р«1 поле температуры в капле равномерно ( Т^ Тс « Т и А Г~ 0 ), в то время как при р» 1 имеем Т5 « Тес1,
Т =
г
(5Т-ЪТе?)/2; ЛГ = 5(Т-Тщ)/2,
Тс =с0
эти соотношения, получим:
где АТ — максимально возможный перепад температуры внутри капли (с течением времени
он стремится к нулю, поскольку Т ^Т„п ).
Таким образом, использование параболического профиля температуры (19), содной стороны, позволяет приближенно учесть неравномерность температурного поля внутри капли, а с другой — практически не приводит кдополнительным вычислениям. Роль градиента температуры в капле
определяется численным значением безразмерно-р
Т//15 — времени релаксации поля температуры внутри капли и — времени установления равновесной температуры капли в результате внешнего теплообмена. Легко показать, что в отсутствие испарения и теплового излучения
/(л^Чиу,)
и
где р — безразмерный параметр, — масштаб времени релаксации поля температуры внутри капли;
х
Р = В1/5,
Перепад температуры между центром и поверхностью капли Тс - Т5 пропорционален перепаду между среднеобъемной и равновесной температурами:
ЛН^Ь—^-^)- (25)
Наконец, подстановка равенства (22) в правую часть уравнения (16) дает уравнение для среднеобъемной температуры капли
где В1 = —число Био.
р
называть критерием Био, модифицированным с учетом испарения.
Подчеркнем, что равенство (25) и уравнение (26) справедливы по истечении начального промежутка времени, продолжительность которого пропорциональна т1. Уравнение (26) показывает, что при наличии градиента температуры в капле масштабом времени релаксации среднеобъемной температуры капли к ее равновесному значению Теч является величина 1 (1 + р). В соответствии с (25) и (26) различие между Тс< Т!. и Т становится малым за время, пропорциональное 1 (1 + р ),и можно приближенно считать, что перепад температуры в капле снижается до весьма малых значений за время, равное Ъхеч (1 + Р).
р
чина (1+ Р) представляет и масштаб времени выравнивания температуры внутри капли. Напомним, что в отсутствие излучения
leg (! + ß) = ^v+T//15'
а в общем случае teq (l + ß) = тед +т,/15.
Модель (20)—(26) на основе параболического температурного профиля (19) имеет важные ограничения. Во-первых, при сколь угодно малом t значения температуры в центре и на поверхности различны, даже если при t = 0 распределение температуры равномерно и TC = TS. Это связано с тем, что данная модель несправедлива на некотором начальном промежутке времени, пропорциональном Т//15- В самом деле, при t << X/ квазистационарный профиль температуры, который в данном случае аппроксимируется параболой (19), еще не успевает сформироваться. Если длительность указанного временного промежутка намного меньше других масштабов времени, то данное ограничение не представляется существенным. Для случая, когда это не так, в работе [5] предложен метод искусственного сглаживания зависимости Ts(t) при малых значениях /, который позволяет с хорошей точностью воспроизвести зависимость температуры поверхности от времени, полученную точным (численным или аналитическим) решением уравнения теплопроводности (17). С учетом рекомендаций [5] в данной работе вместо равенства (22) использовано следующее выражение:
Т=Т + - ß
1 + ß
(Т«,~Т>
1-ехр
-iooVß
т,/15.
(27)
Для температуры центра Tc[t) использован следующий ограничитель:
^) =
|max(7;(/),7-0), Tg>T0;
min
(«То),
Т„<П.
(28)
характеризующих изменение внешних тепловых потоков. Приведенные рассуждения не учитывают циркуляцию жидкости внутри капли.
Циркуляция жидкости внутри капли. Дальнейшее развитие модели позволяет приближенно учесть наличие циркуляции жидкости внутри капли. При этом используется известное аналитическое решение для поля скорости и температуры при осесимметричном течении жидкости внутри сферы — вихрь Хилла (см. например статьи [3] и [1]). В этих работах показано, что увеличение интенсивности теплопередачи от поверхности капли к внутренним слоям жидкости из-за циркуляции жидкости определяется числом Пекле
Ре:
u-uji/
(29)
у Дс/Р/)'
В работе [3] предложено вместо коэффициента теплопроводности жидкости у использовать величину
V =х(Ре)уь (30)
где множитель х(Ре), который определяется как
' Ре"
p (Ре) = 1,86 + 0,86 th
2,2451g
30
(31)
изменяется от 1,00 до 2,72 при изменении числа Пекле от нуля до бесконечности.
Полезно отметить, что равенство (31) почти идентично равенству
р(Ре) = 1,86 + 0,86(Ре/30) (32)
(Ре/30) +1
Поскольку, в отличие от (31), равенство (32) не требует вычисления трансцендентных функций, оно и используется в данной работе. Введение эффективной теплопроводности означает, что масштаб времени релаксации поля температуры в капле (24) принимает вид:
где температура центра Тс в правой части равенства определяется из второго равенства (20).
Другое ограничение связано с тем, что в данной модели предполагается монотонный профиль температуры, в то время как монотонность может нарушаться, например, если капля быстро перемещается между областями с существенно разной температурой газа. В этом случае применимость рассматриваемой модели должна быть оценена с учетом временных масштабов,
Т; =
№
iPi
(33)
Движение капли в газовом потоке. Скорость движения капель в газовом потоке определяется прежде всего силами тяжести и трения. Силы Бассе, Саффмена, Магнуса и роль присоединенной массы не учитываются, поскольку плотность жидкости намного превосходит плотность
газа (более полная совокупность сил, действующих на движущуюся каплю, анализируется, например, в работе [9]). Компоненты скорости и координаты каждой движущейся частицы описываются следующими уравнениями:
du¡ di
.—±CD\ut-ug^(ut-ugJ)-
dx¡ dt
1--Р-
. Р/,
■ = u¡.
(34)
(35)
Коэффициент сопротивления С0 мы определяли стандартным образом в предположении сферичности капель:
0> =
24 (, U 2/3
1 + —Re 6
Re
0,424,
Ре <1000; Ре >1000,
(36)
^velo
4 d
P¡
3 pCJu-u
(37)
и установившейся скорости осаждения частицы в поле тяжести
позволяет записать уравнение (34) в следующем виде:
du¡ dt
(38)
1 velo
Т =TS,k/f
Ч-е/о ^velo I J '
Sik
где x veIo 2
Ре = р5 |и - и^ | ¿//ц, — число Рейнольдса для
капли, вычисленное через ее диаметр и разность между скоростями частицы и потока, р„, и. —
о Г б
плотность и вязкость газа. Введение времени динамической релаксации частиц
М2 8^ )•
Масштабы времени и безразмерные критерии.
В рассматриваемой модели нагрева и испарения капель можно определить три характерных масштаба времени: X/ (33) — время релаксации поля температуры внутри капли; (1 + р)хщ — время прогрева капли до равновесной температуры, определяемое равенствами (16) и (23) (оно же — время выравнивания температуры внутри капли, если х достаточно мал о) и тешр —время жизни капли до ее полного испарения.
Отметим, что в определении масштаба времени X/ неявно сделано допущение о постоянстве коэффициента р , учитывающего увеличение эффективной теплопроводности за счет внутренней циркуляции жидкости в капле. Это означает, что поле скорости в капле считается установившимся, а это, в свою очередь, эквивалентно допущению о малости характерного времени установления поля скорости в капле. Например, для капель //-декана диаметром 0,1 мм в работе [3] показано установление поля скорости за время порядка 0,001 с. Время установления поля скорости в капле определяется размером капли и вязкостью жидкости, и может оказаться, что оно сопоставимо со временем релаксации поля температуры; детальное описание данной стадии процесса требует решения не только уравнения переноса тепла, но и уравнений движения внутри капли. Такой подход не представляется оправданным для целей данной работы, поскольку на стадии установления поля скорости испаряется пренебрежимо малая часть массы капли.
Время жизни капли можно оценить как
Отметим, что формула для коэффициента сопротивления (36) может быть представлена в виде
Сл = С^/(Ре), где С5/ = 24/Ре — сопротивление по Стоксу,
/(Яе) = 1 + Ре2/3/б.
Тогда время релаксации скорости (37) можно записать как
= U
-2
d^/dt
1 е\'ар ^ "0
Используя это определение и равенство (1), будем иметь:
Р Л
"evap
4Sh рШп (1+ ВМ)
(39)
Важно, что в большинстве случаев начальная температура капли не определяет характерную скорость испарения. Поэтому оценки по формуле (39) следует выполнять для температу-
ры мокрого термометра, при которой испаряется большая часть массы капли.
На основе перечисленных масштабов времени можно получить два независимых безразмер-
р
деляемое равенством (23). Как отмечено выше,
р
ределяет величину градиента температуры внутри капли. Еще один критерий можно определить как отношение времени прогрева до равновесной температуры (1 + Р) к времени жизни капли тетр. Можно показать, что в отсутствие теплового излучения указанное отношение выражается как
^Л1+Р) = 2^ сДт;-т„„)р+р) ^ (40)
Х^ар 3 СР шр ^ ) + 9шрВт
где ДЛЯ В] имеет место равенство (5). Если отношение хсоп„(1 + р)/хешр мало, то мал и начальный период, в течение которого температура капли изменяется, а значит большую часть времени капля испаряется при равновесной температуре. В последнем случае для неподвижной капли (К(и~2) имеет место классический закон й , а для быстро летящей капли (N11 ~ Яе^2 ~ — закон [8].
Уравнение движения капли (34) порождает еще один безразмерный критерий — число Сто-
кса 81к = туе/о/т0, равное отношению времени
у-' Slk
динамическои релаксации по Стоксу хш1о к масштабу времени т0, характеризующему скорость изменения скорости газа. Значение числа Стокса определяет инерционность капель и характер их взаимодействия с турбулентными пульсациями скорости (см. например [10]).
Численный алгоритм
Расчет теплообмена и испарения капли выполняется в следующем порядке.
1. Вычисляются давление и массовая доля пара у поверхности капли (при температуре поверхности).
2. Определяются теплофизические свойства жидкости в капле и парогазовой смеси. Свойства парогазовой смеси определяются для температуры и состава, рассчитанных по так называемому правилу одной третьей:
Tref = Ts+(Tg-Ts )/3;
vap ~^mp,s ^ ( ^vap, g s
3. Определяются числа Re, Sc, Pr и Sc0, Nu0 (по формулам (8)).
4. Вычисляется массовое число Сполдинга Вм по формуле (2).
5. Определяются числа Sh\ Nu" по формулам (6) и (7). Используются значения //у в предыдущий момент времени.
6. Вычисляется тепловое число Сполдинга Вт по формуле (5). Шаги 5 и 6 повторяются до сходимости по Вт.
7. Определяется скорость убыли массы m по формуле (1) и новое значение массы капли. Масса капли на новом шаге по времени вычисляется по формуле
mn+l=mn + m" At.
Если при испарении капли оказывается, что mn <At\ т|я,то шаг по времени корректируется так,
что A = mnjпосле этого масса капли становится равной нулю и расчет прекращается.
8. Определяется новое значение диаметра
1 /3
капли dn+x ={втп+хIяр/) . Дальнейшие расчеты на данном шаге выполняются для d = (dn +dn+{)/2;
т = (тп +тп+х )jl.
9. Рассчитывается среднеобъемная температура T"+l по формулам (41) или (42), а также новые значения температуры центра и поверхности капли (по формулам (22) и (20)).
10. Определяются проекции скорости м"+| и
координаты x"+l капли по формуле (43).
В связи с большим различием численных значений временных масштабов ieq и хетр характерных для процесса прогрева и испарения капель, система дифференциальных уравнений (1) и (26) становится жесткой. Абсолютная устойчивость, монотонность и физически обоснованное асимптотическое поведение численного решения уравнения (26) при больших шагах по времени обеспечивается при использовании не-
явной схемы первого порядка. Если положить, что значения Тец, хщ и р слабо меняются на одном шаге, то для среднеобъемной температуры в следующий момент времени будем иметь
т" + тед 9(1+р)) ^ (41)
1 + 9 ((( + р))
где Т" —температура на предыдущем шаге.
Альтернативный подход заключается в аналитическом интегрировании уравнения (26) на одном шаге по времени в предположении «замороженных» значений Тщ, т и р:
Тп+Х = тея + (тп-Теq |ехр
/
At
(1 + i I
(42)
eg J
Численные расчеты показали, что при достаточно малом шаге по времени 9 формулы (41) и (42) дают практически неразличимые результаты. Отметим, что применение указанных формул требует предварительного расчета асимптотического значения температуры Т . Уравнения движения (34), (35) аппроксимируется аналогично:
«Г =
1/А' + 1hvelo
veto
X,
n+l=x"+A t
uf+ыГ
(43)
Численный алгоритм реализован в составе модуля Spray программы Fire3D [2].
Для адекватного численного разрешения всех временных масштабов ((1 + р) xvelo и тешр) шаг по времени должен быть существенно меньше , чем наименьший из них:
Af<mi п((1 + РК9,хуе/о,хеш/)).
Поскольку, как правило, (1 + Р) xeil*Tvel0
и
а+р)
X << X
'■eg ^^ ^evap1
использования в составе гидродинамического решателя, совместно с которым выполняются расчеты очень большого количества капель.
В самом деле, малый шаг по времени
Д?<тш((1 + р)хе9,хуе/0)
необходим лишь в ходе указанной релаксации, после чего шаг по времени следует адаптировать к времени испарения капли при равновесной температуре. Поэтому в данной работе шаг по времени определяется из условия
+ ^ К (44)
^ evap
выполнение указанного условия потребует выполнения неоправданно большого числа шагов по истечении времени релаксации температуры и скорости к своим равновесным значениям. Это особенно важно, если модель предназначена для
_L-1 _
At~c]rni n( (1 + р ),xveto) где С < 1.
Значения xeq (1 + р), хуе/о и хешр оцениваются
в начале расчета. Параметр 5 = |г- - ^ед\
изменяется от 1 до 0 и указывает на степень близости решения (температуры) к асимптотическому пределу. Поскольку (1+ Р) тед~ xvelo, в случае решения уравнений движения капли полагаем, что этот параметр 5 характеризует также и степень релаксации скорости капли. В начале нагрева (или охлаждения) капли 5 ~ 1 и 9t^Cieq (1 + Р), а при достижении равновесной температуры имеем 5^ 0 и 9 « Cxevap. В результате такой адаптации шага по времени к решению оказывается, что при заданном значении С качество численного решения слабо зависит от соотношения хед (1 + Р) и времени жизни капли xevap.
В качестве примера рассмотрим испарение капли воды с начальным диаметром 0,1 мм и начальной температурой жидкости 20 °С. Параметры окружающего газа будем полагать неизменными. При температуре окружающего воздуха 20 °С имеем
Xg (1 + P)/Xevap =0,0064,
а при 1400 °С указанное отношение равно 0,072.
Сравнение рис. 1, а и рис. 2, а, где приведены результаты расчетов скорости испарения капли при высокой (1400 °С) и низкой (20 °С) температурах окружающего воздуха, показывает качественное различие зависимостей скорости испарения dm/cit от времени. Это связано с тем, что температура мокрого термометра в первом случае оказывается выше, а во втором — ниже на-
чальной температуры капли. В первом случае (рис. 1, а) скорость испарения капли меняется немонотонно. На начальном этапе скорость испарения растет, поскольку температура капли увеличивается. Затем, после разогрева капли до температуры мокрого термометра, скорость испарения снижается из-за уменьшения площади поверхности капли. Во втором случае (рис. 2, а) скорость испарения монотонно снижается как на начальном этапе охлаждения капли до температуры мокрого термометра, так и при последующем квазистационарном испарении.
Результаты расчетов (рис. 1—3) показывают, что удовлетворительное качество решения можно получить при С < 0,1, причем решение практически перестает зависеть от шага по вре-
мени при С <0,05. При температуре окружающего воздуха 20 °С отношение характерных масштабов времени тед (1+ i)/тетр значительно меньше, чем при 1400 °С, что делает задачу более жесткой. Несмотря на это, даже при значительной деградации численного решения (в первую очередь это относится к расчетной скорости испарения), которая наблюдается при О 0,2, возможная ошибка расчетного времени жизни капли не превышает 20—30 %. Это объясняется тем, что на начальной стадии испарения (представлена отдельно на рис. 3), когда все параметры изменяются особенно быстро (характерный масштаб времени их изменения равен Teq 0 + i)) испаряется лишь малая часть массы капли.
dm/dt, кг/с
2.5Е-08
1
ö) dm/dt, кг/с
Рис. 1. Результаты численного решения задачи теплообмена и испарения капли (высокая температура окружающей среды): скорость испарения (а) и диаметр капли (б) в зависимости от шага по времени.
Показаны зависимости, полученные при С= 0,002 (/); 0,005; 0,010; 0,020; 0,05; 0,10; 0,20; 0,50 (<У). Вода, d„ = 0,1 мм, Т„ = 20 Т\ Ttlir = 1400 'С
Рис. 2. Результаты численного решения задачи теплообмена и испарения капли (низкая температура окружающей среды): скорость испарения (а) и диаметр капли (б) в зависимости от шага по времени.
Таи. = 20 °С; значения С, с!п, Т() даны в подписи к рис. 1
Рис. 3. Зависимости расчетной скорости испарения капли на его начальной стадии от шага по времени при высокой (а) и низкой (б) температурах окружающей среды (ТаЬ. = 1400 и 20 "С соответственно);
значения С, cln , Т() даны в подписи к рис. 1
Итак, в данной работе сформулирована математическая модель для описания нестационарного и неизотермического испарения капель диспергированной жидкости с учетом неравномерного поля температуры внутри капли и циркуляции в ней жидкости. На основе трех характерных масштабов времени (время релаксации поля температуры в капле, отношение времени релаксации температуры поверхности капли к значению температуры мокрого термометра и время жизни капли до полного испарения) выделены два безразмерных критерия, определяющих величину перепада температуры в капле и его возможное влияние на динамику испарения. Построен абсолютно устойчивый числен-
СПИСОК J
1. Sazhin, S.S. Advanced models of fuel droplet heating and evaporation [Text] / S.S. Sazhin // Progress in Energy and Combustion Science.— 2006,— Vol. 32.-No 2,- P. 162-214.
2. Snegirev, A. Flame suppression by water sprays: flame-spray interaction regimes and governing criteria |Text| / A. Snegirev, A. Lipjainen, V. Talalov // Proc. of the 12th International conference interflam 2010 (Nottingham, UK, 5—7 July 2010). interscience Comm., London.- 2010,- Vol. 1,- P. 189-199.
3. Abramzon, B. Droplet vaporization model for spray combustion calculations [Text] / B. Abramzon, W.A. Sirignano // International Journal of Heat and Mass Transfer.- 1989,- Vol. 32,- No 9,- P. 1605— 1618.
4. Sazhin, S.S. A simplified model for bi-compo-nent droplet heating and evaporation [Text] / S.S. Saz-
ный алгоритм, позволяющий при достаточно малом числе шагов по времени получать численное решение, слабо зависящее от жесткости системы (она определяется отношением времени прогрева к времени жизни капли). Указанный алгоритм предназначен для использования в составе гидродинамического кода Fire3D [2] при численном моделировании турбулентных газокапельных струй.
Работа выполнена в рамках проекта «Универсальная модель газокапельной струи для инженерных приложений» при поддержке РФФИ (грант 10-08-92602-К0_а), Royal Society (Великобритания, грант JP090548) и ООО «Гефест» (Санкт-Петербург).
hin, A. Elwardany, PA. Krutitskii [et al.( // International Journal of Heat and Mass Transfer.— 2010.— Vol. 53,- No 21-22,- P. 4495-4505.
5. Dombrovsky, L.A. A parabolic temperature profile model for heating of droplets [Text] / L.A. Dombrovsky, S.S. Sazhin // ASME Journal of Heat Transfer.- 2003,- Vol. 125,- P. 535-537.
6. Dombrovsky, L.A. A simplified non-isothermal model for droplet heating and evaporation [Text] / L.A. Dombrovsky, S.S. Sazhin // International Communications in Heat and Mass Transfer.— 2003,— Vol. 30,- No 6,- P. 787-796.
7. Watkins, A.P. Modelling the mean temperatures used for calculating heat and mass transfer in sprays | Text| / A.P. Watkins // International Journal of Heat and Fluid Flow.- 2007,- Vol. 28,- No 3,-P. 388-406.
8. Терехов, В.И. Тепломассообмен при испарении капель бинарных растворов [Текст] / В.И. Терехов, Н.Е. Шишкин // Тр. 5-й Рос. нац. конф. по тепломассообмену РНКТ-5,— М.: Изд-во МЭИ, 2010,- Т. 4,- С. 302-305.
I Text j / C.T. Crowe // Boca Raton, USA: CRC Press, 2006,- 1156 p.
10. Sazhin, S.S. Particle grouping in oscillating flows [Text] / S.S. Sazhin, T. Shakked, V. Sobolev, D. Katoshevski // European Journal of Mechanics
9. Crowe, C.T. (Ed.) Multiphase flow handbook B/Fluids.- 2008,- Vol. 27,- No 2,- P. 131-149.
УДК 519.245, 535.34, 535.36
A.M. Кривцун, А.Ю. Сетейкин
АНАЛИЗ ПРОЦЕССОВ РАСПРОСТРАНЕНИЯ ОПТИЧЕСКОГО ИЗЛУЧЕНИЯ В БИОЛОГИЧЕСКИХ СРЕДАХ С ИСПОЛЬЗОВАНИЕМ ВЫЧИСЛЕНИЙ НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ
Моделирование процесса распространения оптического излучения в биологических средах — актуальная задача биомедицинской оптики. Транспорт фотонов в мутных средах, таких как биологические ткани, достаточно хорошо описывается интегрально-дифференциальным уравнением теории переноса излучения [ 1]. Однако аналитическое его решение без введения ряда приближений и упрощений не является тривиальным [2]. Поэтому в оптике биотканей широкое распространение получили численные методы моделирования распространения света [3].
Одним из таких методов является метод Монте-Карло, который широко используется для численного моделирования сложных систем во многих областях науки. К основным его преимуществам над другими методами относятся высокая точность, сравнительная простота реализации, а также возможность работы со средами произвольной пространственной конфигурации. Тем не менее, поскольку метод Монте-Карло стохастический и основывается на многократном повторении испытаний для большого числа частиц, то наиболее явный его недостаток состоит в низкой вычислительной эффективности.
Текущие достижения в области разработки систем на базе многоядерных процессоров сделали возможным преодоление основного недостатка метода [4], позволив тем самым превратить его в удобный инструмент биомедицинской диагностики [5]. В частности, представляется
целесообразным использование методики массивно параллельных вычислений на графических процессорах (ГП). Поскольку расчеты транспорта фотонов в мутных средах осуществляются независимо (частицы в процессе своего распространения не взаимодействуют друг с другом), данный метод хорошо подходит для реализации с использованием массивно параллельных вычислений на ГП.
Метод Монте-Карло
Этот метод широко применяется для численного решения ряда физических и математических задач, так как предлагает гибкий и в то же время достаточно точный подход к моделированию транспорта фотонов в биологических средах [6]. Суть метода состоит в установлении набора общих правил, описывающих поведение фотонов в тканях и определении математических выражений для функций распределения вероятности двух основных характеристик: длины свободного пробега фотонов между двумя последовательными актами взаимодействия, а также направления, в котором фотоны будут продолжать свое движение после акта рассеяния.
Особенно важный этап реализации метода, определяющий его вычислительную эффективность и точность, — это выбор представления моделируемой среды. Наиболее часто биоткань рассматривают в виде многослойной структуры, при этом оптические параметры в пределах каж-
