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

Вычислительные технологии

Том 7, № 3, 2002

ОБ ОДНОМ МЕТОДЕ РЕШЕНИЯ СТАЦИОНАРНОЙ ЗАДАЧИ ОБТЕКАНИЯ

Ю.Н. Захаров Кемеровский государственный университет, Россия e-mail: [email protected]

Problem of liquid flow in infinite chanel is given. Approximate edge conditions are proposed. These conditions are sequents of Navier — Stoks equations set. Difference schemes are approximating this problem. They are bilinear equations set. Method of minimum residuals with multi-step and component-wise optimization is proposed.

Введение

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

Второй проблемой, возникающей при численном решении задачи обтекания в переменных скоростей и давления, является отсутствие в уравнении неразрывности слагаемых с давлением. Поэтому для численного решения стационарных задач можно использовать или аппроксимацию разностными схемами так называемых нестационарных е-систем [5] или рассматривать разностную схему, аппроксимирующую дифференциальную задачу как систему нелинейных уравнений. Применение е-систем приводит к обычным нестационарным разностным схемам, у которых временной т и пространственный h дискретные шаги связаны условием устойчивости. Даже в случае неявных схем т выбирается из условия аппроксимации производных по времени, а не из условий наиболее быстрой сходимости. Второй подход встречается существенно реже. Это связано с тем, что разностные задачи являются системами большой размерности. А это обстоятельство затрудняет применение неявных итерационных методов типа Ньютона — Канторовича из-за необходимости точного обращения на каждом итерационном шаге матрицы Якоби. И, во-вторых, без существенных ограничений на оператор нелинейной системы трудно получить оптимальный итерационный параметр [6].

© Ю. Н. Захаров, 2002.

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

1. Интегральные краевые условия в задаче обтекания

Для системы уравнений Навье — Стокса

где и = (их,и2, и3), Щ = Пг(х, у, г), г = 1, 2, 3; р = р(х, у, г) — давление; ^ > 0; х = (х, у, г), рассмотрим следующую задачу обтекания (рис. 1):

Здесь Г4 — боковая грань параллелепипеда за Г3; X — значение координаты х, задающее прямоугольник Г3, на котором будет задаваться краевое условие.

Будем считать, что задача (1), (2) имеет решение (см. [7]). Проинтегрируем (1а) по области Предполагаем, что граница Г4 так далеко от обтекаемого контура Г2, что на этих плоскостях поток уже невозмущенный [4]. Тогда можно получить следующие равенства:

иУи + Ур = ^Ди,

и = 0,

их|Г1 = и0, и2|Г1 = 0, и3|Г1 = 0, и|Г2 = 0, щ ^ и0, и2 ^ 0, и3 ^ 0 при |х| ^ то. (2)

Рис. 1.

Z 2 Y2 Z 2 Y2

J J uj(X,y,z)dydz + J Jp(X, y, z)dydz = (Y - Yl)(Z2 — Zl)(u°l)2

Zi Yi Zi Yi

Z2 Y2 Z2 Y2

f du2(x, y, z)

ui(X, y, z)u2(X, y, z)dydz — ^

dx

Zi Yi Zi Yi

Z2 Y2 Z2 Y2

f du3(x,y,z)

u1(X, y, z)u3(X, y, z)dydz — ^

dx

dydz = 0, (3)

x=X

dydz = 0.

x=X

Zi Yi Zi Yi

Таким образом, равенства (3) являются следствием системы (1) и явным образом учитывают краевые условия на бесконечности. Следовательно, эти равенства можно использовать в качестве краевых условий на Г3 при решении задачи (1), (2) методом сеток.

2. Итерационная схема

Введем в области П неравномерную согласованную с границей Г2 сетку и на этой сетке аппроксимируем систему (1) разностной схемой второго порядка. Заменяя интегралы в (3) общепринятыми квадратурами не ниже 2-го порядка точности, на границе Г3 получим краевое условие для построенной разностной схемы. В итоге для нахождения приближенного решения в области П получим билинейную систему уравнений вида

A(u,u) = f, (4)

где u = {uUjk,u2ijk,u3ijk,Pijk} = (ui,...,um)T, i £ I, j £ J, к £ K — вектор размерности m, равной учетверенному количеству расчетных точек, принадлежащих области П; f = (fi,... , fm)T — вектор правой части, зависящий от краевых условий исходной дифференциальной задачи; I, J, K — множества индексов, зависящие от величины и характера построенной сетки.

Можно показать, что

A(u, u) = A1(u, u) + A2u,

где A]^(u, u) — билинейный оператор, аппроксимирующий uvu и нелинейную часть (3); A2 — линейный оператор, аппроксимирующий линейные части системы (1) и (3).

Для решения системы (4) мы использовали итерационную схему неполной аппроксимации (см. [8])

un+1/2 = un — Тп+1В[A(un, un) — f ], (5a)

un+1 = un+1/2 + an+1zn, n = 0,1, 2,..., (5b)

где В — заданная матрица с элементами, зависящими от un; zn £ Em — некоторый вектор; u0 — произвольное начальное приближение из области определения оператора A; тп+1 = const, an+1 — матрица итерационных параметров; Em — m-мерное евклидово пространство.

Первый шаг схемы (5) — это какая-либо известная схема с некоторым итерационным параметром тга+1. Элементы матрицы ага+1 мы выбирали с помощью следующего алгоритма.

Пусть ага+1 — диагональная матрица с элементами а(;+1), г = 1,... , т. Перепишем (5Ь) в виде

u

п+1 _ yn+l/2 + a(n+l)e(i)

У(Р-1) + L, ar e( , P (6)

(p

г=р

где у;;-!)2 = и;+1/2 + а(п+1)е(1) + ■ ■ ■ + у^2 = и;+1/2, в« — вектор с одной

ненулевой г-й компонентой — е(г) = (0,... , 0,1, 0,... , 0).

Параметры а(;+1) будем выбирать последовательно из условия

min |r;+1/21| _ min I|A(y£1/2 + a^e« y£1/2 + a^e«) - f II, i _ 1,... ,m. (7)

a - -

гг

Очевидно, что гП+)1/2 _ rn+1 _ A(un+1 ,un+1) — f — невязка схемы (5) (более подробно см. [9]).

m

3. Численные расчеты

С целью проверки скорости сходимости схемы (5) мы ее использовали для нахождения картины течения в плоской квадратной каверне с крышкой, движущейся с постоянной (и1 = 1, и2 = 0) скоростью. Система (1) аппроксимируется разностной схемой 2-го порядка аппроксимации на квадратной сетке с шагом к = 1/Ж. В таблице приведены некоторые результаты расчетов для N =10 различными итерационными схемами. Итерации продолжались до выполнения условия ||гга||/||г°|| < £°.

м Схема ео Число итераций

10 СП 10-6 836

ЯС, Тп 9 х 10-6 81

ЯС, тп, a 0.8 х 10-6 15

10-1 СП 10-3 > 10000

ЯС, Тп 10-3 63

ЯС, тп, a 10-3 146

10-2 СП 10-1 > 10000

ЯС, тп, a 10-3 49

10-3 СП 10-1 530

ЯС, тп, a 10-1 18

В таблице приняты следующие обозначения: СП — схема стабилизирующей поправки для е-системы с т _ const [5], выбираемым экспериментально из условия максимальной скорости сходимости; ЯС, тп — явная схема un+1 _ un — Tn+1rn, rn _ A(un, un) — f — невязка; ЯС, тп, a — явная схема

где тп+1 выбирается из условия min ||rn+1/2||, а ап+1 — диагональная матрица, элементы которой удовлетворяют условиям (7).

Как видно, явная схема (8) с матрицей итерационных параметров сходится существенно быстрее схемы СП.

Качество использования краевого условия (3) мы проверяли на примере двумерной задачи обтекания прямоугольного препятствия плоским равномерным потоком вязкой несжимаемой жидкости.

На рис. ??, ?? приведена картина течения, полученная алгоритмом из п. 2 в области П различной длины. Во всех нижеприведенных расчетах ß = 0.002, размеры прямоугольника 0.2 х 0.3. Он расположен вертикально по центру. Как следует из проведенных расчетов, условие

ди1 дх

на правой границе не выполняется:

0 (9)

|unNj иШ-1 j 1 п г

max-—-— = 0.5,

h

N

где N — количество узлов по х. Поэтому численное решение задачи (1), (2) с условием (9) на правой границе возможно только в том случае, если она достаточно далека от обтекаемого тела. Использование же равенств (3) позволяет уменьшить область решения без существенного изменения самого приближенного решения. В наших расчетах уменьшение области решения на 25 % привело к изменению решения на 2-3 %.

Наряду с уравнениями (3) в некоторых задачах удобно использовать и другое условие замыкания разностной схемы на границе конечной области Г3 (см. также [10]).

Аппроксимируем в П систему (1) какой-либо разностной схемой. Затем, так как за обтекаемым телом также выполняется (1), на тех границах, на которых мы не можем

Рис. 2.

Рис. 3.

Рис. 4.

разумно поставить краевое условие, аппроксимируем (1) внутрь П с тем же порядком, что и исходная разностная схема. Полученная система нелинейных уравнений является системой вида (4). Эта система, как и в предыдущем случае, решалась итерационным методом (8).

Для проверки работоспособности данного подхода была рассмотрена задача обтекания потоком жидкости, набегающим на квадратное препятствие неперпендикулярно к стороне квадрата, т.е. для (1) мы имеем следующие краевые условия:

их |Г1 = 1, и2|Г1 = -0.1, и|Г2 = 0, их ^ 1, и2 ^ -0.1 при г = \/х2 + у2 ^ то. (10)

На рис. 4 приведена картина течения решения задачи (1), (10) расчетной схемы 2-го порядка аппроксимации, когда за краевое условие на прямых Г3 была взята аппроксимация (1) внутрь области.

Как видно из приведенного расчета, под обтекаемым телом поток из наклонного переходит в плоскопараллельный, а на правую границу линии тока подходят с различными углами. Увеличение области решения, как и в предыдущем случае, не приводит к существенному изменению решения в меньшей области.

В заключение автор благодарит С. А. Кривушина, выполнившего часть расчетов.

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

[1] Роуч П. Вычислительная гидродинамика. М.: Мир, 1990. 616 с.

[2] ГОРЕЛОВ Д. Н., ЩЕПАНОВСКАЯ Г. И. Об одном методе решения задачи обтекания тел вязкой жидкостью при больших числах Рейнольдса // Изв. АН СССР. Механика жидкости и газа, 1973. №4.

[3] Orlanski I. A simple boundary condition for unbounded hyperbolice flow // J. Comput. Phys. 1976. Vol. 21, No. 3. P. 251-269.

[4] Пухначев В. В. К асимптотике осесимметрического решения задачи обтекания для уравнений Навье — Стокса // Тр. Всес. сем. по численным методам механики вязкой жидкости. Новосибирск: Наука, 1969. С. 208-211.

[5] ЯнЕнко Н. Н. Метод дробных шагов решения многомерных задач математической физики // Новосибирск: Наука, 1967. 195 с.

[6] Самарский А. А., Николаев е. С. Методы решения сеточных уравнений. М.: Наука, 1978.

[7] Ладыженская О. А. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970.

[8] Захаров Ю. Н., Нагорнова О. Н. Итерационная схема минимальных невязок решения стационарной системы уравнений Навье — Стокса // Проблемы динамики вязкой жидкости. Новосибирск: ИТПМ СО АН СССР, 1985. С. 156-159.

[9] Захаров Ю.Н., Егорова Е. Ф., Толстых М.А., Шокин Ю. И. Метод минимальных невязок решения одного класса нелинейных уравнений // Красноярск, 1991. (Препр. / СО АН СССР. ВЦ; №9).

[10] Захаров Ю. Н. Об одном методе решения уравнений с краевыми условиями на бесконечности // Вычислит. технологии. 1993. Т. 2, №4. C. 56-68.

Поступила в редакцию 16 мая 2002 г.