Научная статья на тему 'Численное решение задач обтекания системы тел в переменных функция тока-завихренность'

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

CC BY
334
73
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ / ОБТЕКАНИЕ СИСТЕМЫ ТЕЛ / УРАВНЕНИЯ НАВЬЕ-СТОКСА В ПЕРЕМЕННЫХ ФУНКЦИЯ ТОКА - ЗАВИХРЕННОСТЬ / ГРАНИЧНЫЕ УСЛОВИЯ / ИНТЕГРАЛЬНЫЕ УСЛОВИЯ ПИРСОНА / NUMERICAL SIMULATION / FLOW AROUND A SYSTEM OF BODIES / NAVIER-STOKES EQUATIONS IN STREAM FUNCTION / VORTICITY VARIABLES / BOUNDARY CONDITIONS / INTEGRAL PEARSON CONDITIONS

Аннотация научной статьи по физике, автор научной работы — Калинин Евгений Игоревич, Мазо Александр Бенцианович

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

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

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.

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

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

Том 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.

дп дп

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Функции из*, из'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.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

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

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

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

i Надоели баннеры? Вы всегда можете отключить рекламу.