Модернизация приближения Обербека-Буссинеска для решения нестационарных задач конвекции в замкнутой полости Текст научной статьи по специальности «Математика»

Модернизация приближения Обербека-Буссинеска для решения нестационарных задач конвекции в замкнутой полости

H.A. Винниченко0, A.B. Уваров6

Московский государственный университет имени М. В. Ломоносова, физический факультет, кафедра молекулярной физики. Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2.

E-mail: а nickvinn@yandex.ru, ь uuarou@phys.msu.ru Статья поступила 17.02.2010, подписана в печать 25.05.2010

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

Ключевые слова: нестационарная конвекция, приближение Обербека-Буссинеска.

УДК: 533.6.011.6:536.24. PACS: 47.55.pb, 44.25,+f.

Введение

Приближение Обербека-Буссинеска (ОБ) [1] оказалось очень продуктивным для анализа конвективной устойчивости. Несмотря на ограничения, накладываемые характерными размерами системы и градиентами параметров [2, 3], это приближение позволяет существенно упростить систему уравнений гидродинамики и провести в ряде случаев даже аналитический анализ. Развитие вычислительных методов только укрепило роль данного приближения, поскольку из-за большой разницы времен распространения звуковых волн и характерных времен, связанных с конвективным движением среды, решение полной системы уравнений сопряжено со значительными затратами времени. Действительно, при моделировании конвективного течения со скоростью 1/0 ~ м/с требуется вести расчет, по крайней мере в течение времени ¿/^о, где Ь — характерный размер системы. При ограничении, накладываемом на шаг по времени для явной схемы (И ^ йг ■ Ь/сц, где Сц — скорость звука, а ¿г -безразмерный шаг сетки, это составляет э • йг) шагов по времени. Таким образом, для типичной сетки, содержащей около 200 узлов по каждой координате, подобная оценка дает 7 • 106 шагов по времени. Применение приближения ОБ и других упрощенных моделей учета сжимаемости среды позволяет использовать более эффективные вычислительные методы и увеличить шаг по времени примерно в 500 раз, соответственно сократив время расчета. Конечно, при значительном увеличении чисел Рэлея даже в стандартных задачах прохождения теплового потока через объем приближение ОБ становится неприменимо (см., например, [4], где сравнение с точным двумерным решением в задаче с нагревом и отводом тепла через боковые стенки показало, что для чисел Рэлея Иа > 3 • 108 отклонения от приближения ОБ становятся существенными), но для не очень больших чисел Рэлея этот метод продолжает

оставаться очень важным и практически очень широко используемым приближением.

В то же время существует класс задач, где приближение ОБ оказывается формально неприменимым даже при Иа<сЗ-108. Такая ситуация возникает при расчете нестационарной конвекции в замкнутых системах, например в детекторе утечек в трубах [5]. Определение утечки основано на измерении разности давлений в двух одинаковых камерах — измерительной и контрольной, в которые газ поступает из основной трубы. Поскольку в ходе процесса набора газа в камеру газ заметно нагревается, то непосредственно после закрытия впускного клапана происходит его охлаждение за счет теплообмена со стенками камеры и возникает конвективное течение. Пока температура газа не выровнялась, его давление также продолжает меняться. Таким образом определяется минимальное время, которое должно пройти с момента набора газа в камеры до момента измерения разности давлений, чтобы измерение было достаточно точным и отражало разность количества газа в камерах, а не разность температур. Для нахождения этого времени в работе [5] используется упрощенный подход с введением эффективного коэффициента теплопередачи, который определяется экспериментально. Современные вычислительные методы позволяют решить эту задачу точно, с учетом уравнений газовой динамики. Далее, на примере этой задачи, которая является, с одной стороны, важной с точки зрения технических приложений, а с другой стороны — классической задачей о конвекции в замкнутой полости, будет показано, что приближение ОБ не описывает зависимость давления, плотности и температуры от времени в случае, если изменением среднего (по пространству) давления нельзя пренебречь. Аналогичная ситуация возникает в задачах, в которых при объемном нагреве замкнутого объема давление внутри системы растет. В настоящей работе предложен простой способ модернизации систе-

ТЕОРЕТИЧЕСКАЯ И МАТЕМАТИЧЕСКАЯ ФИЗИКА

29

мы ОБ, который позволяет сохранить все преимущества метода и в то же время существенно повысить точность расчета.

Приближение ОБ и упрощенная модель сжимаемой жидкости

Рассмотрим полную систему уравнений газовой динамики с учетом вязкости и теплопроводности

^ + (Цу(рс) = О,

йч) 91

Р~сй = ~ Р + ^ + 3 ® + рё'

йТ йр . . _ т] (дт)( дук 2. <9оД2 СрР~сй + \ \дхк + 3 1кдх~1)

ЯТ'

(1)

где р, V, р, Т — плотность, скорость, давление и температура газа, г) и Л — коэффициенты вязкости и теплопроводности, р и ср — молярная масса газа и теплоемкость при постоянном давлении, Я — газовая постоянная, g — ускорение свободного падения, — символ Кронекера, й/сН = д/д1 + (»V).

В приближении ОБ считается, что изменением давления и связанными с ним изменениями плотности можно пренебречь во всех уравнениях системы (1) за исключением уравнения движения. Кроме того, поскольку скорость конвективных течений мала по сравнению со скоростью звука, пренебрегают третьим членом в правой части уравнения движения и вторым — в правой части уравнения энергии. Работой силы тяжести р£а можно пренебречь, если gL/{cpв) <С 1, где Ь — характерный размер системы, а в — характерная разность температур. Таким образом, получаем уравнения газовой динамики в приближении ОБ:

(Ну® = О,

йV

аТ \ ЛТ

(2)

(И

Ш р КТ,

где ро — значение давления в невозмущенной среде.

Именно это приближение обычно используется для расчетов конвективных течений, так как оно не учитывает процесс распространения звуковых волн, но достаточно точно описывает его результат — унос массы из нагреваемых областей — при условии, что давление остается примерно постоянным и одинаковым (с точностью до гидродинамического перепада давления ро2<Сро) во всей области. Величина шага по времени при решении уравнений (2) определяется не скоростью распространения звуковых волн с5, как в случае системы (1), а скоростью конвективного течения о0 -С Сц, что позволяет значительно сократить время расчета.

В случае нестационарной конвекции в замкнутой полости пространственные градиенты давления малы, однако изменением давления с течением времени пренебречь нельзя, поэтому прямое применение приближения ОБ приводит к ошибочным результатам. Дей-

ствительно, использование уравнения состояния в виде р = рор/(ЯТ), например в рассматриваемой задаче об остывании газа в полости с холодными стенками, приводит к нефизическому увеличению массы газа, содержащегося в полости. Приближение ОБ предполагает «в среднем» изобарный процесс, в то время как охлаждение или нагрев газа в замкнутой полости — «в среднем» изохорный процесс. Легко оценить величину ошибочного изменения массы 5т при описании с помощью приближения ОБ поведения газа, которому было передано количество тепла С}:

РоУ=-ЙТи р0У = Я = с0т(Т2-Т1),

р р

где V — объем полости, пг — масса газа, Т\ и Т2 — начальная и конечная температура газа, с0 — теплоемкость при постоянном объеме. Из записанных уравнений состояния и теплового баланса следует, что 6т = т(Т2 — Т\)/Т2 = (¡/(с0Т2). Заметим, что относительная погрешность 8т/т не зависит от размеров полости. Таким образом, даже для большой, но замкнутой полости применение приближения ОБ приводит к ошибочным результатам.

Нарушение сохранения массы при расчете нестационарной конвекции в замкнутой полости характерно и для более точных моделей, в которых пренебрегают полной производной давления по времени в уравнении энергии, например для упрощенной модели сжимаемой жидкости (УМСЖ):

др

+ (НУ(/Э®) = О,

(3)

т

йч) 91

= ^ р + ^ + з " + ре'

йТ . . _ 91 (дюк 2 . дюЛ

СрР~Ш = + 2 \dx~k + ~дх~1 ~ 3 к~дх~1) + р=РР

Р кТ

При решении системы уравнений (3) на каждом шаге по времени сначала из уравнения энергии определяется температура. Затем из уравнения состояния находится оценка плотности на новом шаге, которая используется для аппроксимации производной по времени в уравнении непрерывности. Такой подход позволяет не описывать распространение звуковых волн и пользоваться всеми преимуществами приближения ОБ, одновременно учитывая часть эффектов, связанных со сжимаемостью газа. УМСЖ дает результаты, очень близкие к полученным с использованием полной системы уравнений, при описании течений, связанных с выделением тепла, в неограниченных областях, например при расчете влияния релаксации на структуру вихревой дорожки Кармана в колебательно-неравновесном газе [6]. Для замкнутой области, давление в которой заметно меняется при нагреве или охлаждении, пренебрежение членом др!д1 в уравнении энергии (членом (г>У)р можно пренебречь) приводит к таким же ошибочным результатам, как и в случае приближения ОБ.

Следует также отметить, что неверными оказываются не только распределения плотности и давления, полученные при расчете конвекции в замкнутой по-

15 ВМУ. Физика. Астрономия. М 6

лости в рамках упрощенных моделей, но и распределение температуры, а следовательно, и поле скорости конвективного течения. Действительно, «в среднем» изохорный процесс изменения температуры газа в замкнутой полости не описывается уравнением энергии СррйТ/еИ = АДГ, принятым в приближении ОБ.

Модификация упрощенных моделей

Выше было показано, что приближения ОБ и УМСЖ при расчете нестационарной конвекции в замкнутой полости приводят к ошибочным результатам, а именно:

1) не сохраняется масса газа, заключенного в полости;

2) уравнение энергии вида сррс1Т/сИ = АДГ неправильно описывает процесс изменения температуры при постоянном объеме. Соответственно необходимо модифицировать эти приближения так, чтобы при расчете течений с небольшим числом Маха получались достаточно точные значения. При этом необходимо сохранить общую структуру уравнений и вычислительного алгоритма, так как она определяет допустимую величину шага по времени и эффективность расчета.

Будем считать, что давление в полости примерно одинаково во всех точках (ри2 <сро), но зависит от времени. То есть пространственными градиентами давления и связанными с ними изменениями плотности, как и в приближении ОБ, можно пренебречь, но частной производной давления по времени пренебрегать нельзя. Если также учесть, что при изохорном процессе среднее значение давления р~Г и др/д1« Шр/р) ■ дТ!д1, то левая часть уравнения энергии из системы (1) примет вид

dT dp дТ dp

дТ

» cvp— + CpP{vV)T. (4)

Таким образом, уравнение энергии, как и в приближениях ОБ и УМСЖ, можно разрешить относительно температуры, но оно учитывает изменение среднего давления с течением времени.

Чтобы найти поправки к плотности и давлению, обеспечивающие сохранение массы, представим плотность и давление в виде p(r,t) = pi(r,t) + 8p(t), p(r, t) = pi (r, 0 + Sp(t), где pi (r, t) и pi (r, t) подчиняются уравнениям приближения ОБ (или, если более полно учитывать влияние сжимаемости газа, УМСЖ), 8р(0) = 8р(0) = 0. Чтобы с учетом поправок 8p(t) и 8p(t) выполнялись закон сохранения массы и уравнение состояния, необходимо выполнение равенств

(pi(r, t) + 8p(t))dV =

pinit(r)dV,

8p(t) dV =

' fJ'Spjt) RT(r, t)

d V,

где рш\{г) — распределение плотности в начальный момент времени. Получаем

8p(t) =

SMr)-Pl(r,t))dV JdV '

_ R S(p-m«(r)^pi(r,t))dV

,, r dV ■

M J W7i

(5)

Таким образом, предлагаемая модификация приближений ОБ и УМСЖ при расчете нестационарной конвекции в замкнутой полости заключается 1) в использовании левой части уравнения энергии вида (4), 2) в добавке к найденным распределениям плотности и давления в конце каждого шага по времени поправок, определяемых по формулам (5).

Постановка задачи и численные методы

Для проверки точности упрощенных моделей было проведено двумерное численное моделирование охлаждения газа в квадратной полости с холодными стенками, соответствующее упрощенному описанию процессов, происходящих в камерах детектора утечек [5]. Для решения полной системы уравнений (1) был применен метод Годунова второго порядка точности на основе схемы Родионова [7] с точным решением задачи распада произвольного разрыва и Т\Ю-ограничителем Ван Лира. Решение уравнений приближений ОБ (2) и УМСЖ (3) было выполнено с помощью полунеявной схемы [6, 8]. Интегрирование по времени осуществлялось с помощью метода Рунге-Кутта третьего порядка, а также схемы Кранка-Николсона для диссипативных членов. При аппроксимации пространственных производных использовались центральные разности. Для одновременного решения уравнений непрерывности и движения в случае упрощенных моделей был применен проекционный метод определения давления.

Были заданы следующие параметры задачи: постоянная температура стенок = 300 К, максимальная начальная температура газа Гтах = 360 К, сторона полости Ь = 5.2 мм, начальное давление газа рШ\ = 5 атм, все параметры среды соответствовали воздуху. Длина стороны полости была специально уменьшена в 10 раз по сравнению с параметрами детектора утечек [5], чтобы понизить число Рэлея и от турбулентной конвекции, имеющей место в эксперименте, перейти к ламинарной. Таким образом, вопрос о точности используемых уравнений отделяется от проблемы моделирования турбулентности в реальном течении. На стенках полости были заданы следующие граничные условия: v = 0, Т = Гтап, в случае полной модели (1) требовалось дополнительное условие для давления др/дп = 0. Поскольку, в отличие от метода Годунова, полунеявная схема [6, 8] теряет устойчивость в областях с большими градиентами параметров, начальные условия во всех случаях были заданы в виде

V = 0, Т(х, у) = Ттш - (Гтах - rwaü) ,

Р = Рш, р{х,у)=р

RT(x,y)'

(6)

где 8(х,у) — расстояние до ближайшей стенки, d = 0.2L — толщина переходного слоя. Моделирование было выполнено для сеток 100x 100, 200x200 и 300 х 300 узлов для упрощенных моделей и для сеток 50 х 50, 100 х 100 и 170 х 170 узлов для полной системы уравнений (1). Шаг сетки во всех случаях уменьшался при приближении к стенкам полости для более точного описания областей с большим градиентом температуры. Полунеявная схема [6, 8] характери-

ТЕОРЕТИЧЕСКАЯ И МАТЕМАТИЧЕСКАЯ ФИЗИКА

31

зуется значительно меньшей, по сравнению с методом Годунова, численной диссипацией, поэтому для нее результаты, полученные на разных сетках, имеют лишь незначительные отличия. При расчете на разных сетках по методу Годунова, напротив, наблюдаются значительные отличия: величина скорости конвективного течения при расчете на грубой сетке оказывается значительно меньше из-за влияния численной диссипации. Даже при расчете на максимально подробной из использованных сеток 170 х 170 узлов с минимальным безразмерным шагом 2.5 • 10~3 это влияние еще заметно, однако дальнейшее увеличение числа узлов невозможно из-за слишком большого времени расчета (расчет полной системы уравнений (1) занимает более трех недель на персональном компьютере).

Результаты численного моделирования

На рис. 1 показано изменение давления в центре полости с течением времени для разных моделей и сеток. Приближения ОБ и УМСЖ без модификаций (кривая 1) не описывают изменение давления, в них среднее давление в полости остается постоянным за счет нефизического увеличения массы газа. Кривые 2 и 3, полученные на разных сетках для полной системы уравнений (1), практически совпадают, так как величина численной диссипации оказывает слабое влияние на эту зависимость. Максимальное отклонение от полной модели кривых, полученных для приближений ОБ и УМСЖ с учетом предлагаемых модификаций (кривые 6 и 7), составляет 3% от наблюдаемых

0.02 0.04

0.14 t, с

Рис. 1. Зависимость давления в центре полости от времени: кривая 1 — УМСЖ без модификаций; 2 — полная система (1), сетка 50 х 50; 3 — полная система (1), сетка 170 х 170; 4 — приближение ОБ с учетом (5), но без учета (4); 5 — УМСЖ с учетом (5), но без учета (4); 6 — приближение ОБ с учетом (4) и (5);

7 - УМСЖ с учетом (4) и (5)

изменений давления и может объясняться как приближенным характером уравнений, так и различиями использованных численных методов и сеток. Оценить влияние использования левой части уравнения энергии в виде (4) можно по отклонению кривых для моделей ОБ и УМСЖ с коррекцией плотности и давления на каждом шаге по времени, но без модификации уравнения энергии (кривые 4 и 5). Это отклонение составля-

0 0.2

0 0.2

Рис. 2. Распределение вертикальной компоненты скорости (м/с) в момент / = 0.14 с: а - УМСЖ с учетом

модификаций; б — полная система (1), сетка 170 х 170

Рис. 3. Распределение температуры (К) в момент / = 0.14 с: а - УМСЖ с учетом модификаций; б — полная

система (1), сетка 170 х 170

16 ВМУ. Физика. Астрономия. JS» 6

ет 9%, что подтверждает необходимость использования данной модификации для точного описания скорости теплообмена. Поскольку число Маха образующегося конвективного течения очень мало (Мага 6- 10^5), то результаты, полученные при использовании приближений ОБ и УМСЖ, практически не отличаются.

На рис. 2 и 3 приведены рассчитанные поля вертикальной компоненты скорости и температуры для момента ¿ = 0.14 с, когда скорость конвективного течения достигает максимального значения. Результаты, полученные с помощью модифицированной упрощенной модели сжимаемой жидкости и при решении полной системы уравнений (1), очень близки. Разность максимальных значений температуры (333.3 К для полной модели и 330.7 К для УМСЖ) составляет 4.3% начальной разности температур между газом и стенками полости. Несколько меньшее значение максимальной вертикальной скорости для полной модели (0.0178 м/с) по сравнению с моделью УМСЖ (0.0206 м/с) объясняется, вероятно, влиянием численной диссипации при решении полной системы (1).

Выводы

Использование упрощенных моделей сжимаемости среды в случае медленных конвективных течений позволяет избежать необходимости описывать быстрый процесс распространения звуковых волн, что значительно (в 500 раз) снижает время расчета. Показано, что применение упрощенных моделей, в которых прене-брегается производной давления по времени в уравнении энергии, приводит к ошибочным результатам в случае расчета нестационарной конвекции в замкнутой полости, когда нельзя пренебречь изменением среднего

давления. В частности, показано, что в приближениях ОБ и УМСЖ не сохраняется масса газа в полости и недостаточно точно описывается скорость теплообмена. Предложена модификация этих приближений, устраняющая указанные недостатки. Проведено численное моделирование охлаждения газа в квадратной полости с холодными стенками в рамках различных моделей. Отличия результатов, полученных для упрощенных моделей с учетом предложенных модификаций, от результатов для полной системы уравнений составляет 3-4%, что находится в пределах погрешности использованных численных методов.

Работа выполнена при частичной финансовой поддержке РФФИ (грант 09-08-00961-а).

Список литературы

1. Гершуни Г.З., Жуховицкий Е.М. Конвективная устойчивость несжимаемой жидкости. М., 1972.

2. Уваров A.B., Осипов А.И. // Вестн. Моск. ун-та. Физ. Астрой. 1998. № 5. С. 42.

3. Zeytounian R.Kh. // Comptes Rendus Mecanique. 2003. 331, N 8. P. 575.

4. Mlaouah H., Tsuji T., Nagano Y. // Int. J. Heat and Fluid Flow. 1997. 18, N 1. P. 100.

5. Harus L.G., Cai M., Kawashima K., Kagawa T. // Meas. Sei. Technol. 2006. 17, N 2. P. 411.

6. Osipov A.I., Uvarov A.V., Vinnichenko N.A. // Phys. Fluids. 2006. 18, N 10. P. 105106.

7. Куликовский A.Г., Погорелое H.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М., 2001.

8. Kim J., Kim D., Choi H. 11 J. Comput. Phys. 2001. 171, N 1. P. 132.

Improvement of Oberbeck-Boussinesq approximation for solving unstationary convection problems in a closed cavity

N.A. Vinnichenko0, A.V. Uvarov'

Department of Molecular Physics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia.

E-mail: a nickvinn@yandex.ru, b uvarov@phys.msu.ru.

Validity of using simplified models of compressible fluid for calculations of unstationary convection flows inside closed cavities is tested for the case study problem of hot gas cooling in a rectangular cavity with cold walls. Oberbeck-Boussinesq approximation and other simplified models are shown to yield incorrect values of pressure, density and temperature if average pressure change is not negligible. Modification of numerical method is proposed, which allows to describe temporal dependence of pressure correctly without loss of computational efficiency. Results obtained by using modified Oberbeck-Boussinesq approximation and the complete equations for fully compressible fluid are compared to each other.

Keywords: unstationary convection, Oberbeck-Boussinesq approximation. PACS: 47.55.pb, 44.25.+f. Received 17 February 2010.

English version: Moscow University Physics Bulletin 6(2010).

Сведения об авторах

1. Винниченко Николай Аркадьевич — канд. физ.-мат. наук, науч. сотр.; тел.: (495) 939-27-41, e-mail: nickvinn@yandex.ru.

2. Уваров Александр Викторович — докт. физ.-мат. наук, профессор, профессор; тел.: (495) 939-26-94, e-mail: uvarov@phys.msu.ru.