УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
Том 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