Расщепление по физическим процессам для расчета задач конвекции в двумерных областях с криволинейной границей Текст научной статьи по специальности «Физика»

УДК 532.516

О.Н. Гончарова

Расщепление по физическим процессам для расчета задач конвекции в двумерных областях с криволинейной границей

Для исследования конвективных движений жидкости в двумерных областях с твердой непроницаемой границей, часть которой является криволинейной вставкой, предлагается численный метод, основанный на идее расщепления по физическим процессам. Расщепление по физическим процессам в задачах гидродинамики базируется на методе слабой аппроксимации и обосновании аддитивности этих процессов при достаточно малых шагах по времени [1-3].

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

Тестирование метода проводится на известных задачах о свободной конвекции в полости при подогреве сбоку.

Введение. Для моделирования тепловой гравитационной конвекции в замкнутых областях О с твердой непроницаемой границей £ применяются уравнения Обербека-Буссинеска, которые в безразмерной форме имеют следующий вид:

щ + V • Ур = —Ур' + —ДV + р, (1)

Кв

ШуУ=0, (2)

Тг + Р • УТ = —^—АТ. (3)

КвРг

Задача состоит в определении скорости жидкости V, давления р и температуры Т, удовлетворяющим при Ь = О начальным условиям

V = 0, Т = Т0(х), х € О,

на границе области £ — условиям прилипания для скорости

V = 0, х € £, Ь € [0,Ь*\

и условиям первого, второго или третьего рода для температуры:

дТ

т = тв( х,г), = Тв( х,г),

дТ

— + а(Т — Тв) = 0; х € £, Ь ^0,^].

дп

р'

клонение от гидростатического давления); р =

— {Ог/Кв)рТ - сила плавучести; р - единичный вектор ускорения силы тяжести; а - безразмерный коэффициент, характеризующий теплообмен с внешней средой.

Конвективный и диффузионный переносы выделяются в уравнениях движения (1) и в уравнении переноса тепла (3). Идея расщепления по физическим процессам, как основа для организации расчета свободной конвекции, реализуется на каждом временном слое. Метод расщепления по физическим процессам базируется на методе слабой аппроксимации [1], обосновании аддитивности процессов для достаточно малых шагов по времени [2] и доказательстве суммарной аппроксимации исходного уравнения вследствие расщепления. Общая теория расщепления наиболее полно изложена в [3], а особенности расщепления для задач конвекции в прямоугольных областях и параллелепипедах - в [46].

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

уравнений. На этапе диффузии осуществляется переход к новым искомым функциям: вихрю и функции тока. Данный переход позволяет исключить расчет градиента давления, обеспечить автоматическое выполнение условия со-леноидальности вектора скорости, так как диффузионная скорость принимается в итоге в качестве искомой. Для функции тока считаются выполненными следствия условий прилипания для физической скорости. В дальнейшем будем опускать кавычки в терминах конвективная или диффузионная скорость. Отметим, что принципиальным моментом в организации расчета является введение смещенных сеток [7].

Тестирование метода проводится на известной задаче о свободной конвекции в квадратной полости при подогреве сбоку [8].

Постановка разностной задачи. Пусть искомые функции У, У, Т удовлетворяют в области О = {(х,у) : 0 < х < 1, 0 < у < Цх)} системе уравнений (1)-(3). Перепишем данные уравнения в виде

dv -»

——|- Kv = —Vp' + Dv + F, dt

div v = 0,

dT

— + KT = DT. dt

(4)

(5)

(6)

—Kv,

(7)

где V - конвективная скорость. При построении разностной аппроксимации оператор конвективного переноса будем рассматривать в виде суммы двух одномерных кососимметричных операторов: КУ = КУ + К2У [3].

Этап диффузии определяется уравнением

V — V

-----= —Ур' + БУ + Е,

которое в результате применения операции rotor преобразуется к виду:

Du + (rot F, k).

(8)

Здесь под V понимается диффузионная скорость процесса и используются обозначения и = (пй V, к),и = (пЛУ,к), V = (щ,й2), й = дф/ду, й2 = — дф/дх. При этом функция тока ф удовлетворяет в области О уравнению Пуассона:

Дф = —и,

(9)

которое решается в дальнейшем методом установления [9]. На границе области имеют место условия

дф

ф = 0, ^ = 0,

дп

(10)

представляющие собой следствия условий прилипания для физической скорости.

Расщепление на конвективный и диффузионный этапы в уравнении (3) приводит к последовательному решению следующих разностных уравнений:

TT

TT

KT,

DT.

(Н)

(12)

КК V •У Б - оператор диффузионного переноса Б = РД, V = 1/Де или V = 1/(ДеРт).

В уравнении (4) выделим два различных процесса переноса количества движения: конвективный и диффузионный. Считая эти процессы аддитивными на малых временных интервалах Ьк < Ь < гк+1 (т = ¿й+1 — Ьк), сформулируем задачу конвективного движения жидкости, исходя из нижеследующих двух этапов (опустим стрелочку над векторами). Назовем этапом конвекции движение, описываемое уравнениями

Для решения задачи в области О с криволинейной границей осуществим, как и в [10], ото-

{<

£ < 1, 0 < п < 1}; применив преобразование £ = х, п = у/Цх)- Перепишем уравнения (8), (9), (12) в переменных £,п в следующем общем виде:

ф,= B

* f

Щ + Vv

+ R,

(13)

где Ф = и, Ф = Ф (крышка опускается) или Т

и = Вц Ф5 + В12 Фп , V = В12 Ф5 + В22 Фп,

а коэффициенты вычисляются следующим образом:

B

її

f, B12 — —Vf'

В22 = (1 + п1 '*)/I, В = 1/Ее,

Д = (Ст/Ее2) (дТ/д£ — п(1'/1)дТ/дп),

если Ф = и. Если же Ф = ф, то В = А и Д = Аи (А — итерационный параметр). Если Ф = Т, то В = 1/(ДеРт), а Д = 0.

КК

ных £, п имеют вид:

Kv = —

dv d{y\v)

v'Tí + —

vv

I'

К2'0 = ^

/ и \дv

(щ — П1 ^)дП] +

+ 3^2 — п1 'Щ^

В области О вводится основная разностная сетка (£п,Пт) для расчета температуры и реализации этапа диффузии (9), (12): £п = (и — 1 )Н^, п = 1, 2, ...,Щ Цш = (т — 1 )^^ т = 1,2,М;

= ^, Ж = Ж+1;Ьп = ^! М = М + 1. В расчетной области О вводятся также две вспомогательных разностных сетки: (£,п,пш')> (т' = т + 1/2, и = 1,..., Ж, т' = 1,..., М) - для расчета первой компоненты конвективной скорости VI, и (£п,Пш), и' = п+ 1/2, и' = 1,...,М, т =

, ... , М

вективной скорости Щ . ЗдеСЬ Щ , V2 - компоненты конвективной скорости V.

. 1' ~ к+3/4 п

+ = О,

к+1 _—к+З,/4

-———' + КЧ

/

к+1 I 4

-г3

.

(14)

Здесь КЧ, КЧ представляют собой конечноразностную аппроксимацию конвективных опе-

и, т

(п', т) или (и, т') рассматривается аппроксимация, использующая центральные разности вида

1

т^Ъ ~ -1 ( к -

К1 = V11+1 / 2 3 -г+1,3

—°1 г-1 /2,3,3

КЪ

1*2 Ьу

°2 ¿,3+1/2 тг,3+1

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

1 фп,Ш1+1 фп

1п Ьп

° г,з-1 /г -г,3—

/¿+1 — I— П3 П 3+ Пз+1 П 1з+1

" 2Не 2

йг,3+г\-

/¿+1 — I— пз + Пз-1 П г— ,

+------^----------------------~------------™г,з-\

2Н,

"^2 п+1 /2,ш

ф

п ,ш

- фп

пш

к

I п фп,ш+1 фп

и

Ьп

Для решения уравнений (7), (11) используется двуциклический способ расщепления по направлениям на основе элементарных схем Кранка-Николсона [3]. Представим данную схему для искомой функции й на произвольной сетке (г^), которая может быть основной (и,т) для температуры либо одной из вспомогательных (и',т), (и,т') для компонент конвективной

скорости:

к+1У4 _ --к ¿3 Шг3

т/

КЪ

к+1У4 йк

¿3 т *3

,

к+1/2 _ -к+1/4: г3 —3

т/

к+1^2 ~к+\/А

+ КЪ1~гЗ ^¿3

>+

. I — к^-1 / 4 + й°-г3

к+3/4_ -к+1/2 г3 —3

т/

+ къ1~13

,

к/

-

к+1/2

Значения компонент скорости с предыдущего к

ждой расчетной сетки в зависимости от искомой функции. Заметим, что диффузионная скорость (щ, й2) уже заранее рассчитана на основи, т

ется аппроксимация направленными разностями аналогично тому, как это было сделано в [5]. Заметим также, что при расчете компонент конвективной скорости на смещенных сетках учитываются в коэффициентах условия прилипания при выходе за границу вспомогательной сетки и при попадании на границу основной сетки. Используемая аппроксимация конвективных операторов обеспечивает кососимметричность разностных операторов конвекции. Полученные системы разностных уравнений в (14) решаются с помощью метода прогонки. Хотя для матрицы системы и не выполнены условия диагонального преобладания, метод прогонки устойчив [11]. Вышеизложенные аппроксимации опеК , К

пользуемая схема обладает свойством энергетической нейтральности поля скоростей, т.е. сохранением нормы скорости в 1,2 (для доказательства см.: [4]).

¿3

Реализация этапа диффузии. При переходе на этап диффузии насчитывается вихрь конвективной скорости и на основной сетке

^2 п+1/2,т ^2 п—1 /2,т

ч

/п + 1 1п — 1 1

Пт

/Ь и

^2 п,т+1 /2 _ '^'2 п,т—1 /2

1 ^1 п,т+1/2 ^1 п,т—1 /2

^ Ьп

и именно он берется в качестве начального усло-

ик и

Этап диффузии реализуется с использованием метода стабилизирующей поправки. Представим разностную схему второго порядка аппроксимации для решения уравнения (13) в следующем виде [1, 10]:

фй+І^ _ фк в

Т /

и ке + V к+1/2п

+

фк+і _ фк+іЛ в

т /

и к+\ _ и к£

(15)

Здесь используются представления

ик+1 = Ви Ф^1 + В12 ,

и к = вп фк + в12ф к,

V к+1А2 = вш Фк + в22 ъкг?1/2

и аппроксимация входящих в (15) производных следующего вида:

дик

((в)п„

$к+1п+1,т _фк+\

п+1,т

к

е

_( В)

Фк+1пш- Ък+1п-

п—1 ,т

11)

+

2Не

т

12,

ч

п+1,т

+

к/

п+1,т+1

к/

п+1,т—1

Ч В)

12 .

п — 1 ,т

к/

п—1 ,т+1

к/

п—1 ,т—1

2НГ,

ф

веденная схема для организации итерационного процесса.

Воспользовавшись условиями (10) и уравнением (9), получим аналог условий Тома для вихря (см. также: [10]). В частности, на границе пт щий вид:

фп,1 = 0, ип,1 = — Ь Vn,2 ,

ЬП ¿п

Vn,2 — —Щ{J 0 п( ф^п,2 +

+

1 + п22(/О

/п

/\2

■-(ф

где для первых производных функции тока используются аппроксимации второго порядка. Заметим, что условия превращаются в классические условия Тома в случае прямолинейной границы Кх) = 1.

Результаты численного анализа. Общая схема решения задачи состоит в осуществлении следующих этапов.

1. Переход на новый временной слой Ьк+1

Тк

на основной сетке по схеме (14) для этапа конвекции и по схеме (15) с соответствующими температурными граничными условиями на этапе диффузии.

Тк

расчет конвективной скорости vk+1 на

вспомогательных сетках и вихря скорости ик

ки по схеме (15) вводится внутренний итера-

ф

чания итераций при в = S считается, что с заданной точностью е^ определены зна-фк фк+г = фв .

4. При переходе на следующий временной слой (к этапу 1) искомая скорость полагается равной V (она известна на основной сетке и может быть насчитана на вспомогательных) .

ф

шимся, если выполнен критерий сходимости [5,

8]

тах фп,т _ фп,тI < • тах \фп,т\-

п. т. ’ ’ п. т. ’

к

Рис. 1: Линии тока. Изотермы (справа)

Рис. 2: Линии тока. Изотермы (справа)

В качестве основного теста рассматривается задача о свободной конвекции в полости 0 < х <1; 0 < у < Цх) при подогреве сбоку. При этом Яв = Рг = 1, От = 104, а вектор силы тяжести направлен против оси Оу. В начальный момент времени жидкость покоится и нагрета по закону Т(х,у,§) = х. Граничные условия соответствуют твердым непроницаемым границам с заданной температурой Т(х,у,Ь) = х . Тест заключается в получении стационарного решения до удовлетворения неравенства [4, 5]

\Мйк+1 — Жйк \ < е,

где число Нуссельта Жй определяет интегральный поток тепла на правой (горячей) границе х

е = 10-, еф = 10-, т = то = Н2/4, Н = .

Результаты тестирования алгоритма, а именно: характеристика интенсивности течения фшах и интегральный тепловой поток на правой стенке х = 1 в случае прямолинейной границы Кх) = 1 приводятся в таблице, в последней строке которой представлены предельные значения характеристик из [8].

Сетка фтах Мп

9x9 7.3558 1.9405

11 х 11 7.0415 1.8886

17 х 17 6.6450 1.8057

21 х 21 6.5460 1.7831

33 х 33 6.4360 1.7607

41 х 41 6.4115 1.7539

[8] 6.37 1.752

Устойчивость алгоритма и экспериментальный порядок сходимости проверяются путем вы-

числительного эксперимента на последовательностях сеток 1 lx 11, 21 х21, 41x41, при этом ведутся наблюдения за величиной т\, характеризующей для каждой сетки интенсивность движения max \фп m\. Экспериментальный порядок

п, m ’

сходимости r и приближенно определенная относительная погрешность вычисляются, согласно правилу Рунге (см., например: [12]). Таким способом определяется r « 1.9 и относительная погрешность ё « 2.8%.

Для криволинейных границ, задаваемых как

!(х) = 1 + ОЛвш^х и Кх) = —0.1х + 1.1, результаты расчетов приведены на рисунках 1 и 2, соответственно. На рисунках представлены семейства линий тока и изотерм.

Таким образом, проведено построение и тестирование численного алгоритма для исследования задач конвекции в двумерных областях с криволинейной границей. Отметим также, что реализация этапа диффузии может быть осуществлена на основе метода прогонки с параметрами [5].

Библиографический список

1. Яненко, H.H. Метод дробных шагов решения многомерных задач математической физики / H.H. Яненко. - Новосибирск, 1967.

2. Самарский, A.A. О принципе аддитивности для построения экономичных разностных схем / A.A. Самарский // ДАН СССР. 1965. -Т. 165, №6.

3. Марчук, Г.И. Методы расщепления / Г.И. Марчук. - М., 1988.

4. Воеводин, А.Ф. Метод расчета вязких течений в замкнутых областях / А.Ф. Воеводин, Т.В. Протопопова // Сиб. журнал индустриальной математики. - 2001. - Т. 4, .V' I.

5. Воеводин, А.Ф. Метод расчета двумерных задач конвекции на основе расщепления по физическим процессам / А.Ф. Воеводин, О.Н. Гончарова // Вычислительные технологии. - 2002.

- Т. 7, Ж.

6. Гончарова, О.Н. Метод расщепления по физическим процессам для расчета трехмерных задач конвекции / О.Н. Гончарова // Известия АлтГУ. 2007. - Вып. 1 (53).

7. Белолипецкий, В.М. Математическое мо-

делирование течений стратифицированной жидкости / В.М. Белолипецкий, В.Ю. Костюк, Ю.И. Шокин. - Новосибирск, 1991. 176 с.

8. Тарунин, Е.Л. Вычислительный эксперимент в задачах свободной конвекции / Е.Л. Тарунин. - Иркутск, 1990.

9. Самарский, A.A. Теория разностных схем / A.A. Самарский. - М., 1977.

10. Овчарова, A.C. Метод решения термоконвективной задачи в многослойной среде с криволинейными границами раздела / A.C. Овчарова // Динамика сплошной среды : сб. науч. тр. АН СССР / Сиб. отд-ние. Ин-т гидродинамики. - 1994. - Вып. 106.

11. Воеводин, А.Ф. О корректности метода прогонки для разностных уравнений / А.Ф. Воеводин // Численные методы механики сплошных сред: Сб. науч. тр. АН СССР / Сиб. отд-ние. Ин-т теор. и прикл. механики. - 1972. -Т. 3, №5.

12. Калиткин, H.H. Численные методы / H.H. Калиткин. - М., 1978.