Программный комплекс "Тепловой взрыв" (стр. 84 )

Аналогичная производная, действующая по y, аппроксимируется точно так же. Для поиска на границе области используется односторонняя разность 2-го порядка, аналогичная (12.54).

12.3.4. Оценка временного шага

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

Характерные масштабы времени, связанные с теплопроводностью и диффузией имеют вид:, где Δr –минимальный шаг по пространственной координате. Если нет никаких других процессов, кроме теплопроводности и диффузии, то приведенный выше алгоритм устойчив при любом временном шаге. Величина шага влияет только на точность решения. Опыт расчетов показывает, что для подобных задач, удовлетворительные результаты получаются при временных шагах порядка , где коэффициенты К изменяются от 1 до 10, в зависимости от требуемой точности. Если коэффициенты переноса зависят от температуры, то для оценки временного шага выбирается максимальная величина.

Характерное время, связанное с конвекцией, определяется как , где V – максимальное значение проекций скорости. Поскольку в начальный момент времени жидкость неподвижна, а условия на границе области и механизм химического превращения могут быть самыми разнообразными, то для первого шага используется самая жесткая оценка:, где Н – максимальная протяженность расчетной области. Опыт проведенных расчетов задач с конвекцией показывает, что описанный выше алгоритм, устойчив при достаточно больших временных шагах порядка , где KV может достигать 50 –80. Однако, точность расчета нестационарных процессов при таких больших шагах по времени очень низкая. Удовлетворительные результаты получаются для KV не более 5.

Теперь рассмотрим, какие ограничения на шаг по времени накладывает наличие химической реакции. Учитывая, что алгоритм включает расщепление на "химический" и "тепловой" шаги и содержит всего одну итерацию по нелинейностям, потребуем, чтобы на "химическом" шаге относительные приращения концентраций и температуры не превышали некоторой заданной величины Kchem и чтобы при дроблении шага невязка, построенная по температуре и концентрации так же не превышала заданных значений KT, KC. Чтобы сократить объем вычислений, все перечисленные условия можно проверять только для точки с максимальным тепловыделением.

Таким образом, стратегия автоматического выбора шага следующая: сначала оцениваются τλ, τD и τV и из них выбирается минимальный шаг. Затем, определяется ячейка сетки с максимальным тепловыделением и на "химическом" шаге для этой точки проверяются следующие условия:

               (12.64)

Если неравенство (12.64) не выполняются, то временной шаг дробится до тех пор, пока это неравенство не станет справедливым. Следует отметить, что если в результате расчетов требуется увеличить шаг по времени, то нельзя это делать слишком резко, например, шаг не должен возрастать более, чем в 1.5 –2 раза.

12.4. Верификация алгоритма моделирования

С целью оценки адекватности и точности моделирования теплопереноса при естественной конвекции химически инертных жидкостей были рассмотрены 5 задач. Результаты решения двух из них (о конвекции в замкнутой прямоугольной полости при Ra=2.2·104, H/L=1, Pr=0.71 и в вертикальном цилиндре, подогреваемом сбоку, при Ra=1·109, H/Rin=2, Pr=1) сопоставлялись с опубликованными данными численных расчетов [65, 36]. Решения трех других задач (о конвекции в прямоугольной вертикальной щели при Ra=4·105, H/L=20, Pr = 1000 и о конвекции в вертикальном цилиндре при Ra=1·1013, H/Rin=6, Pr ≈ 5, Ra = 1·1011, H/Rin=2, Pr=319) сравнивались с результатами экспериментов [66-68].

На примере решения задачи об определении критических условий теплового взрыва в условиях свободной конвекции в вертикальном цилиндре проведена комплексная проверка разработанных моделей и алгоритмов. Данная задача решена с использованием граничных условий 1-го рода для реакции 1-го порядка при Ra=104 – 107, H/Rin=2, Pr=300. Полученные при этом результаты сопоставлены с интерполяционной формулой, обобщающей данные экспериментов.

12.4.1. Тест: нестационарная тепловая конвекция в прямоугольной полости

Подавляющее большинство опубликованных как экспериментальных, так и численных исследований течений со свободной конвекцией выполнено для прямоугольных областей. В связи с этим, в качестве эталонных данных для тестирования разработанного алгоритма были выбраны результаты работы [64], в которой в рамках полных уравнений Навье-Стокса для сжимаемого газа численно рассмотрено течение в замкнутой квадратной полости в поле силы тяжести, направленной вдоль вертикальной оси, при малых числах Маха. На боковых стенках полости заданы постоянные различные температуры T2>T1, а торцы полости считаются теплоизолированными. В начальный момент времени газ покоится, а распределение температуры предполагается линейным, (соответствует стационарному распределению в отсутствии силы тяжести).

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

               

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

Задача решалась на равномерной пространственной сетке 40 Ч 40 для следующих значений перечисленных параметров: , , .

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

На рис. 12.3. приведены результаты расчета числа Нуссельта от времени, полученные в настоящей работе при различных значениях числа Куранта KC и аналогичная зависимость из работы [64].

Рис. 12.3. Зависимость среднего числа Нуссельта от времени: 1 – результаты расчетов для сжимаемого газа из [64], 2 – KC = 5, 3 – KC = 50.

Видно, что при больших значениях числа Куранта нестационарные участки кривых сильно различаются, но стационарные значения числа Нуссельта, полученные в настоящей работе, практически совпадают результатами [64]. Что касается расчета переходного процесса, то для его достаточно точного описания с помощью разработанного алгоритма значение числа Куранта не должно превышать 1 – 1.5.

12.4.2. Тест: тепловая конвекция в вертикальной прямоугольной щели

Данное течение исследовалось экспериментально в работе [64]. В качестве рабочей среды в экспериментах использовались медицинский парафин и силиконовое масло, которые помещались в прямоугольную кювету (см. рис. 12.4). Боковые стенки кюветы поддерживались при постоянных, но различных температурах.

Рис. 12.4. Схема кюветы для исследования конвекции в вертикальной щели в экспериментах [64]

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

Высота уровня жидкости H в кювете была намного больше расстояния между боковыми стенками L и глубины кюветы B, а отношение L/B было меньше 1. Как отмечается в [64] при таком соотношении пропорций кюветы, при ламинарных режимах конвекции трехмерными эффектами (изменением параметров по оси Z) можно пренебречь. С учетом этого обстоятельства расчет описанного течения в настоящей работе проводился в рамках двумерной постановки задачи. При этом использовалась равномерная сетка по x (число узлов по равнялось 40) и неравномерная (со сгущением у стенок) сетка по y (ее минимальный шаг был равен шагу по x, параметр сгущения составлял 1.1, а полное число узлов – 96).