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

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

Том 151, кн. 3

Физико-математические пауки

2009

УДК 532.516.5

ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧ ОБТЕКАНИЯ СИСТЕМЫ ТЕЛ В ПЕРЕМЕННЫХ ФУНКЦИЯ ТОКА - ЗАВИХРЕННОСТЬ

ЕЛ. Калинин, А.Б. Мазо

Аннотация

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

Ключевые слова: численное моделирование, обтекание системы тел. уравнения Навье Стокса в переменных функция тока завихренность, граничные условия, интегральные условия Пирсона.

Введение

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

Значения ф на границах можно задать из геометрических соображений лишь для простейших симметричных стационарных течений. Однако при моделировании сложных нестационарных течений требуются дополнительные соотношения, определяющие значения функции тока на каждом контуре. В работах [1 3] и ряде других применяются вычислительные схемы, использующие нелокальные условия. полученные интегрированием условий Пирсона [4]: теоретические аспекты рассмотрены в монографии [5].

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

1. Постановка задачи

Рассмотрим течение жидкости в области В = [О, Ь] х [О, Н]. Границы области течения обозначим следующим образом: входное и выходное сечения 7;п и 7oi.it < верхняя стенка канала 70, нижняя стенка 7дг+1. границы обтекаемых

Ю

УN+1

Рис. 1. Схема течения

тол 7.;,, г = 1,..., N (рис. 1). В декартовой системе координат система уравнений Навье Стокса в переменных ф — си имеет вид

+ ^ (1)

д , дг>у д1'х дф дф

—Аф = со, и = — -— Г, - = (2)

где г>х, г'у горизонтальная и вертикальная компоненты скорости жидкости V. 11е число Ройнольдса. слагаемое ^ в правой части уравнения (1) определяется через вектор массовых сил g = (дх,ду) как

дх ду

На твердых стенках зададим условие прилипания

х,у£ = Ф = С\(1), г = 0,...,ЛГ+1, (4)

где Сг набор неизвестных функций, а Ц, заданная касательная скорость на контуре 7;. Для определения Сг воспользуемся соотношением Пирсона

др 1 дю

х,у&л: тг = "5_т--г = 0,..., N + 1, о)

ов Ке оп

в котором п и в внешняя нормаль и касательная к границе. Следуя [5]. проинтегрируем (5) по границе 7¿. Будем иметь

дои \

£~КедЛ <1з = Ве\р], г = 0,..., ЛГ + 1, (6)

где через [р] обозначен скачок давления вдоль контура 7;,. На замкнутых контурах, очевидно, [р] = 0. а на стенках канала величина [р] равна перепаду давления на входном и выходном сечениях. Систему соотношений (6) будем использовать для определения неизвестных С г, г = 0,..., N 1. Не нарушая общности, можно положить значение функции тока на нижней стенке равным нулю: Сдг+1 = 0. Тогда величина Со. очевидно, будет равна расходу жидкости в канале.

В настоящей статье рассматриваются течения в прямоугольных областях, где поперечная составляющая скорости набегающего потока в окрестности входного сечения равна нулю: г>у = 0. дг>у/дх « 0. На границе 7;п скорость г>х, функция тока ф. завихренность ю и расход Со будут связаны соотношениями

н

п [ , дф д2ф дпх

( .Г''11»- —„• (,)

о

Представим эпюру входной скорости в виде г'х = г>хСо. где функция г>х(у) нормированный профиль скорости, интеграл которого по границе 7;п равен единице. В дальнейшем изложении на входном сечении будут рассматриваться однородный поток жидкости, для которого

МУ) = (8)

а также течения с параболическим профилем

6 9 6

г'ЛУ) = "язУ" + -±рУ- (9)

Используя соотношения (7). можно записать граничные условия на границе 7;п в виде

х, У € 7ш : Ф = Рф (у) С0, ш = Ри (у) Со, (10)

у

/Зл)

<'< '/.</• Рш(у) = —

о

На выходном сечении поставим «мягкие» граничные условия

дф дю

х, У € 7ой, : = 0, — = 0. (И)

Требуется найти решение системы уравнений (1). (2) с граничными условиями (4). (10). (11) и дополнительными нелокальными соотношениями (6).

2. Метод численного решения

Проведем дискретизацию с шагом г по времени уравнения (1) таким образом, чтобы избавиться от нелинейности уравнения переноса завихренности:

ш - -^-аи; = -туг ■ уси + т^1 + си. (12)

Ке

—Аф = Со. (13)

Здесь использовано стандартное обозначение й(1) = и(1 — т). На каждом временном слое уравнения (12). (13) решаются с граничными условиями (4) на твердых стенках 7;,. причем условие Дирихле используется для уравнения (13). а с помощью условий Неймана для ф формулируются граничные условия первого рода для ю

х,у€ц: си = и>, г= 0,...,М+1 (14)

Для отыскания функции ¿о решается линейная задача [6]

¿> = -Аф, х, у € ъ : ^ = г = 0,..., N + 1. (15)

Заметим, что функция тока ф на текущем временном слое определяется набором констант Сг. поэтому и завихренность ю будет зависеть от этих констант. Задачи (12). (13). (15) линейны, следовательно, можно записать линейное представление искомых функций через С г:

n n

ф = Ф* + ^ ш =+12ш<с<> (16)

¿=0 ¿=0

В котором функции ф* , фг , и!* , Шг подложат определению. Множители С г будут найдены из системы уравнений

n

X1и/; • /' • ¿ = о,...,лг, (17)

3=0

полученной с помощью подстановки представления (16) в интегральные соотношения (6). Коэффициенты алгебраической системы (17) равны

Для определения функций ф* , фг сформулируем вспомогательные задачи. Подставим представление функции тока из (16) в линейное уравнение Пуассона (13). Будем иметь:

n

-Аф* - Аф'''С.1, = ш.

Поскольку правая часть полученного равенства не зависит от констант . можно разбить его на N + 2 уравнения следующим образом:

-Аф* = -Аф1 = 0, г = 0,..., АГ. (18)

Граничные условия Дирихле для уравнений (18) получаются подстановкой разложения (16) в (4) и (10):

n

х, у € ъ : ф* + Е ¥С-з = Си

3 = 0 n

х, у е 7ш : Ф* + Е ¥С-з = Р-Ф (у) Со-

3 = 0

Приравнивая члены при С г. получим

х, у €= 7< : ф* = 0, фз = 6*,

х,у£ут: ф*= 0, ф? = 5{]Р.ф(у),

г^=0,...,М. (19)

Здесь д? символ Кроиекера. Аналогично получаем граничные условия на выходном сечении:

дф* дф1 х, У € 7ой, : = о, = 0. (20)

Вспомогательные задачи (18) (20) не содержат констант : кроме того, уравнения для фг не зависят от времени, поэтому могут быть решены один раз до начала счета.

Аналогично, подставляя (16) в (12). получим уравнения для вспомогательных функций и>* , и>г:

Аш* =Ъ\ Аш1 = 0, 1 = 0,..., И, (21)

где

А = Е--^—А, Ъ = Со - тч ■ У Со + тЛ Ке

Граничные условия для уравнений (21) получаются в результате подстановки (16) в исходные граничные условия (4). (10). (11) для из и приравнивания коэффициентов при С;,:

и;*=йз*, изг = йз\

х, у £ 7т : о;*=0, шг = б?Рш{у), ^

диз* диз'1

х, У € : = 0, — = 0.

дп дп

Функции из*, из'1 определяются как решения задач

д I * ^ дф*

: -Аф , X, у € 7д- . = -Уз

л • дф1

-Аф\ х, у &Ъ: 0,

(23)

следующих из (15).

В задачах (21) (23) для у', i = 0,..., N. правые части равны нулю, а граничные условия зависят только от фг. которые, как отмечалось ранее, не зависят от времени. Следовательно, и функции из'1 могут быть найдены до начала вычислений: непосредственно на каждом временном слое нужно решить задачи лишь для функций ф* , из*. после чего найти константы С г из системы (17). Тогда значения искомых функций ф, из на слое будут получены по формулам (16).

Для определения коэффициентов Ц , I* системы (17) необходимо найти значения нормальных производных диз'1/дп, диз*/дп на твердых стенках. Ниже представлен способ вычисления этих производных в рамках метода конечных элемен-

Выделнм узел конечноэлементной сетки с координатами (ж.;,у.;), лежащий на границе 7^ . Вычислим значения диз*/дп в этом узле. Для этого домножим первое уравнение (21) па пробную функцию ф^, равную единице в выделенном узле и нулю во всех остальных. После интегрирования по области придем к соотношению

^ = V У {ш*ф*+ '~ Ьф<) сЮ = ' (24)

7 к С

где след базисной функции ф^ на границе 7^. Чтобы найти производную

диз*/дп в выделенном узле г, интеграл в левой части этого уравнения вычислим с помощью квадратурной формулы

J f{s)ds {в^Гц, h^ = Jф'lds,

в которой в; узлы интегрирования, совпадающие с узлами сетки, а Й; весовые коэффициенты. В результате получим

[ди*.^¿Ц* ^Зш*.,. диз* ^

у *>«Е^ Ъ = Е з^* = —>':■ (2°) г г

Заменяя левую часть (24) на (25). найдем:

д'П К;

Аналогично получим выражения для определения узловых значений dujf /дп нормальных производных на твердых стенках:

д-£ = Ъ M = + (27)

d

Итак, алгоритм численного решения нестационарной задачи (1). (2) состоит из шага инициализации и цикла по времени.

На шаге инициализации сначала решается 7V+1 задача (18) (20) для вспомогательных функций фг. Эти функции подставляются в соотношения (15). из которых находятся функции шг, определяющие граничные значения шг на 7j (22). Сами функции w', i = 0,..., N находятся как решения задач (21). (22). По формулам (27) подсчитываются производные дш1 /дп, после чего вычисляются коэффициенты If системы уравнений (17).

На каждом шаге временного цикла сначала решается одна задача (18) (20) для вспомогательной функции ф* . Она подставляется в соотношение (15). из которого находится функция ¿>*. определяющая граничные значения ю* на 7j (22). Сама функция iv* находится как решение задачи (21). (22). По формуле (26) подсчи-тываотся производная дю*/дп, после чего вычисляются коэффициенты I* правой части системы уравнений (17). Решение этой алгебраической системы дает набор констант Cj. которые вместе с найденными вспомогательными функциями фг, и>г, ф*, и>* позволяют найти решение ф, и> задачи на текущем временном слое по формулам (16).

3. Тестирование метода

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

3.1. Течение с периодическим изменением функции тока на теле [3].

Это течение в области D = [—Зл", Зтт] х [—Зл", Зтт] с помещенным в ее центр телом 7 в виде квадрата со стороной 2п определено функциями

фе (X, y,t) = (COS X + COS у + COS X COS у) COS t,

(28)

LVe (x, y, t) = —Афе = (cos X + cos у + 2 COS X cos y) COS t.

При соответствующем выборе слагаемого F эти функции удовлетворяют уравнениям (1). (2). периодическим граничным условиям прилипания

7 : фе = Се(1) = - cost, ^=0, (29)

а также нелокальному соотношению

/§^Ь = 0. ,30,

7

В качестве начальных условий и условий Дирихле на внешних границах области для функций фе и <ле будем использовать точные значения (28).

Табл. 1

Погрешность числеппого решепня при t = '2тг

N IIV' - Фе\\с \\Ф-ФАЬ W - из£\\с W - U3e 2 \С - Се\

880 2.166 • 1СГ1 7.340 • 10"- 2.345 • 10"1 5.370 • 10"- 1.429 • 10"-

8448 2.085 • 10"- 7.074 • 10"3 5.978 • 10"- 7.660 • 10"3 1.394 • 10"3

13120 1.321 • 10"- 4.501 • 10"3 3.731 • 10"- 4.848 • 10"3 8.986 • 10"4

18816 9.118 • 10"3 3.107 • 10"3 2.666 • 10"- 3.376 • 10"3 6.271 • 10"4

Табл. 2

Погрешность числеппого решения при t = 47т

N \\Ф-Фе\\с \\Ф-фа\2 W - изе\\с w - U3e 2 \С -С е\

880 2.255 • 10"1 7.798 • 10"- 2.340 • 10"1 6.420 • 10"- 1.403 • 10"-

8448 2.243 • 10"- 7.656 • 10"3 5.863 • 10"- 9.421 • 10"3 1.385 • 10"3

13120 1.421 • 10"- 4.868 • 10"3 3.819 • 10"- 6.026 • 10"3 8.894 • 10"4

18816 9.804 • 10"3 3.355 • 10"3 2.788 • 10"- 4.211 • 10"3 6.192 • 10"4

Тестовый расчет производился с числом Рейнольдса Re = 102 на регулярной сетке билинейных конечных элементов с шагом по времени г = 0.02. На каждом временном слое значение С подсчитывал ось с помощью условия (30) по описанной ранее процедуре и сравнивалось с точным значением Се = — cos t.

В табл. 1. 2 представлены максимальные || • ||с и среднеквадратичные || • Ц2 отклонения сеточных решений ф, из от точных фе, изе. а также ошибки определения па моменты t = 2л", 4л". Видно, что погрешность в определении константы С не возрастает со временем и составляет от 1.5% для сетки с числом узлов N = 880 до 0.06% для сетки с числом узлов N = 18816.

3.2. Обтекание системы вращающихся круговых цилиндров [9]. Это

вынужденное течение в области D = [0, 80] х [0,40] с помещенными в нее двумя цилиндрами единичного диаметра, с центрами в точках (40,38.75). (40,41.25). На входе задан однородный поток жидкости (8). Скорость вращения цилиндров задана таким образом, что V\ = — Vn, при этом вращение нижнего цилиндра направлено против часовой стрелки.

Был проведен ряд тестовых расчетов при Re = 102. 0 < \Vj\ < 3. При этом течение имеет автоколебательный характер. Для сопоставления результатов расчета с полученными в [9] определялись осредиениые по периоду колебаний коэффициенты сопротивления Cci и подъемной силы С;, а также число Струхаля Sh. В работе [9] эти параметры течения рассчитывались с помощью метода конечных объемов для решения уравнений Навье Стокса в естественных переменных.

Задача решалась с шагом по времени г = 0.005 на неструктурированной сетке из 31715 билинейных элементов (N = 31452) с локальным сгущением узлов около обтекаемых тел. В качестве начального условия выбиралось состояние покоя. При заданной скорости на входном сечении расход жидкости в канале Со явно вычисляется по формуле (7). поэтому в задачах вынужденного обтекания применять условие (6) для границы 70 не нужно, достаточно заменить соответствующее уравнение системы (17).

Картина течения при \Vi \ = 1.3 представлена на рис. 2. Изменение величин С; со временем для различных скоростей вращения показано на рис. 3.

Отклонение полученных значений коэффициента сопротивления Cci от результатов. представленных в [9]. составило не более 4%. подъемной силы С; не более 0.86%. числа Струхаля Sh не более 1.26%.

■ЯЁШЯВИИВ

■НЗйЭвИв

Рис. 2. Линии тока и поле завихренности течения вокруг пары вращающихся круговых цилиндров при Re = 100, \ V¿\ = 1.3

Рис. 3. Изменение величин Сг со временем при 11е = 100: а) \Уг\ = 0; б) \Уг\ = 1.3; в) | Уг | = 2. Сплошная линия соответствует верхнему цилиндру, пунктирная - нижнему

3.3. Естественная конвекция в канале около нагретого цилиндра

[10]. Для иллюстрации работы алгоритма в случае, когда расход жидкости Со неизвестен заранее и меняется со временем, моделировалась естественная конвекция в вертикальном канале Б = [0,20] х [0,2.5] с помещенным в него квадратным нагревателем. На поверхности нагревателя температура Т = 1, стенки канала теплоизолированы, а на входном сечении Т = 0.

На стенках канала и поверхности нагревателя ставятся условия прилипания (4), а на входном сечении задается параболический нормированный профиль скорости (9). Скачок давления вдоль стенки канала в интегральном условии (6) принимается равным нулю. Вектор массовых сил, входящий в систему уравнений На-вье-Стокса (1)-(3) и граничные условия (5), (6), в приближении Буссинеска равен g = (Т, 0) [11]. Для определения температурного поля дополнительно решается уравнение конвективной теплопроводности

дТ ~dt

+ v•VT =

1

Pr Re

-А Т.

Здесь Рг — число Прандтля.

Расчеты проводились при Рг = 0.72, 10 < Re < 500 на сетке, содержащей 16690 узлов. Для аппроксимации конвективного слагаемого в уравнении переноса завихренности использовался ТУБ-иодход с ограничителем МС, описанный В [7, 8].

Наблюдается хорошее согласование картины установившегося течения со стационарным решением, полученным в работе [10] (рис. 4, в, г). Разница в значениях

Рис. 4. Линии тока и поле температур конвективного течения при 11е = 100, Рг = 0.72 на моменты: а) £ = 10; б) £ = 20; в) £ = 50 в сравнении со стационарным решением [10] (г)

расхода Со составляет не более 2%, а отклонения чисел Нуссельта

Г дТ 7 N11 = / т—сЬ

} дп

7

не превышает 2.5%.

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

Работа выполнена при финансовой поддержке РФФИ (проекты № 08-01-00548-а, 08-01-00163-а, 07-01-00499-а).

Summary

E.I. Kalinin, A.B. Mazo. Numerical Simulation of Flow around a System of Bodies in Stream Function - Vorticity Variables.

A numerical algorithm for simulating non-stationary viscous 2D flows based on solving Navier - Stokes equations in stream function - vorticity variables is developed. The algorithm allows to model the forced flows as well as the complex convective flows around a system of bodies in channel. Problem definition is amplified with nonlocal integral equalities based on Pearson conditions in order to obtain boundary values of stream function.

Key words: numerical simulation, flow around a system of bodies, Navier - Stokes equations in stream function - vorticity variables, boundary conditions, integral Pearson conditions.

Литература

1. Mizukami A. A stream fuuctiou-vorticit.y finite element formulation for Navier Stokes equations in mult.iconnect.ed domain // Int. J. Numer. Methods. Eng. 1983. No 19. P. 1403 1409.

2. Tezduyar T. Finite element formulation for the vorticit.y-st.ream function form of the incompressible Euler equations on multi-connected domains // Comput. Methods. Appl. Mech. Eng. 1989. V. 73, No 3. P. 331 340.

3. Liu J.-G., Wang C. High order finite difference methods for unsteady incompressible flows in multi-connected domains // Computers & Fluids. 2004. No 33. P. 223 255.

4. Флетчер К. Вычислительная гидродинамика. Ч. 2. М.: Мир. 1991. 552 с.

5. Gluvinski R. Finite element methods for incompressible viscous flow // Handbook of numerical analysis. V. 9. Numerical Methods for Fluids (Part 3). Elsevier Science B.V., 2003. 1176 p.

6. Мазо А.Б., Даутоо Р.З. О граничных условиях для уравнений Навье Стокса в переменных функция тока завихренность при моделировании обтекания системы тел // Ипжеперпо-физ. жури. 2005. Т. 78, Л' 2. С. 75 79.

7. Kuzmin D., Turek S. High-resolution FEM-TVD schemes based on a fully multidimensional flux limiter // J. Сотр. Pliys. 2004. No 198. P. 131 158.

8. Мазо А.В., Калинин Е.И. Решение задач обтекания в переменных «функция тока вихрь» методом конечных элементов с применением TVD-подхода // Модели и методы аэродинамики: Материалы 8-й междупар. школы-семипара. М.: МЦНМО, 2008. С. 100 101.

9. Hyun Sik Yuan, Но Hwan Chun, Jeung Hu Kim, I.L. Ryung Park Flow characteristics of two rotating side-by-side circular cylinder // Computers & Fluids. 2009. No 38. P. 466 474.

10. Khodary K., Bhattacharyya Т.К. Optimum natural convection from square cylinder in vertical channel // Int. J. Heat Fluid Flow. 2006. No 27. P. 167 180.

11. Мазо А.Б. Численное моделирование свободной конвекции вязкой жидкости в канале с нагретым цилиндром // Учен. зап. Казан, ун-та. Сер. Физ.-матем. пауки. 2005. Т. 147, кп. 3. С. 141 147.

Поступила в редакцию 27.04.09

Калинин Евгений Игоревич аспирант кафедры аэрогидромехапики Казанского государственного университета. Е-шаП: kalinineieyantlex.ru

Мазо Александр Венцианович доктор физико-математических паук, профессор кафедры аэрогидромехапики Казанского государственного университета. Е-шаП: amazoQksu.ru