Сильно неявный попеременно-треугольный метод с использованием B-сплайнов для расчета течений вязкой жидкости Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ Ц А Г И

Т о м XV 1 98 4 №5

УДК 517.9 : 532.5 : 533.7

СИЛЬНО НЕЯВНЫЙ ПОПЕРЕМЕННО-ТРЕУГОЛЬНЫЙ МЕТОД С ИСПОЛЬЗОВАНИЕМ Я-СПЛАЙНОВ ДЛЯ РАСЧЕТА ТЕЧЕНИЙ ВЯЗКОЙ ЖИДКОСТИ

В. А. Казаков

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

Большинство методов решения уравнений Навье—Стокса для несжимаемой жидкости являются итерационными, а уравнения для вихря и функции тока или для скорости и давления решаются отдельно по явным или неявным схемам. При этом возникает необходимость так или иначе задавать физически неопределенные значения вихря или давления на твердой поверхности, что приводит к дополнительному процессу и значительно снижает устойчивость и скорость сходимости численного решения (см., например, [1, 2]). Преодолению этой трудности способствовало введение «искусственной сжимаемости» [3, 4]. Развитие этой идеи описывается в [5]. Получающаяся система решается с помощью расщепления по физическим переменным, что в общем также требует итерационного расчета поля скоростей и поля давления.

В работе |[6] предложен сильно неявный метод, в котором функция тока и вихрь связаны на границе, а уравнения движения, аппроксимируемые с помощью центральных разностей со вторым порядком точности, решаются совместно. Данный алгоритм основывается на сильно неявном попеременно-треугольном методе неполной факторизации, предложенном для решения параболических и эллиптических уравнений [7—11] и развитого в других работах (указания на них имеются в [6]). В |[12] этим методом отдельно решались уравнения для вихря и функции тока в задаче о расчете течения вязкой жидкости около цилиндра. И в этом случае численное решение получалось в 2,5 раза быстрее, чем, например, при использовании метода переменных направлений. Совместное решение уравнений сильно неявным попеременно-треугольным методом еще более эффективно: повышается

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

Дальнейшее улучшение свойств численных схем может быть связано с увеличением точности аппроксимационных формул. В настоящей работе разрабатывается сильно неявный метод, в котором приближенное решение ищется в виде двумерных кубических сплайнов, представленных через В-сплайны [13]. Как показано, например, в [14— 17], применение сплайнов позволяет получать более точное решение и использовать более крупные сетки. Кроме того, существенно облегчается постановка граничных условий, связанных с производными. Число итераций и время счета уменьшаются. Основным недостатком описываемого алгоритма является необходимость дополнительных массивов для обращения матрицы.

1. Рассмотрим сначала сильно неявный метод для уравнения Лапласа иХх + Иуу = О в Прямоугольнике с расчетной сеткой Х0<*1< ... <Хм, Уо<У1< ■ ■ ■ <Ук. Приближенное решение будем искать в виде кубического сплайна, представленного через В-сплайны. Для их построения дополним исходную сетку узлами с симметричными относительно х0 и Ху значениями х:

ДГ_3 < Х-2 <С Х—1 < X, Хм < ЛГдг41 < Хц + 2 <С -Клг + 3>

и аналогично по оси у. Обозначим систему В-сплайнов на сетке {хг} через В1 (х), I — — 1,0,..., Ы, на сетке {у ^ через В] (у), у'= = —1, 0, . . . , К+1. Для функции и (х, у) можно записать

лч-1 к+1

и(х,у£ ^В1{х)В]{у). (1)

<=-1 /=-1

Значения приближенного решения уравнения в точке (хИ уу)

обозначим Иц. Так как каждый В-сплайн В1(х) отличен от нуля

лишь на интервале (л:/_2, -*7+2) [аналогично и для В! (у)], то 1 1

иц=‘ ^ ^ $'1+1.1+т В1+1 (л;г) В,-+т (у^.

/ = — 1 т- — 1

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

1 1

У. 2 ^‘+1, 1+т [^‘+г (ъ) &1+т (у]) -\- &1 + 1 {х1)^]+т (_у;.)] = 0, (2)

/ = -1 ГП=-1

1=1, 2,..., N — 1, У = 1, 2,..., к- 1.

В общем случае вместо (2) можно рассмотреть систему уравнений относительно неизвестных следующего вида:

Ф1-1, /-1 + А2 Ф1, /_ 1 + А3 Ф/+1, /_1 -Ь А4 Ф,_1, / 4- А5 Ф^- +

+ •Лб Ф<+1, / + А-! Ф,_х, /+1 + А8 Фл /_]_! + Ад Ф,+1, /+1 = (3)

Здесь у коэффициентов АI — А§ индексы г, у опущены, для уравнения Лапласа <3,7 = 0.

Определим вектор Ф = (Ф00, Ф

10» • • ■ > Флго, Ф0|» ■ ■ ■ > Фдтс)г И аналогично вектор Тогда данную систему можно записать в матричной форме

А Ф - <?, (4)

где матрица А имеет вид, показанный на рис. 1.

Прибавив к обеим частям уравнения (4) вектор РФ, получим общую итерационную схему

(А + Р)Фп+' = 0_ + РФ". (5)

Матрица Р выбирается таким образом, чтобы матрица А + Р факторизовалась — представлялась в виде произведения нижней Ь и верхней и треугольных разреженных матриц: А + Р = Ш. В результате получается система

ьи <Ьп+1 = (%=•(} + РФ*, (6)

которая решается за два шага:

££=(?', иФя+х = 2.

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

Ао А,

I и

Рис. 1

матрицу Ь так, чтобы ее ненулевые элементы находились на тех же диагоналях, что и в нижней части исходной матрицы А (рис. 1). Аналогично составим верхнюю матрицу II. Таким образом, получается следующий алгоритм решения системы (6):

^у=:Я5 (Сіі; Д1 /-1 + Ь/—1 вц/їі — 1, у), І

Фп = -£,7 + 51 п Ф/+і, і +■ 52г■ Ф/_і, /-|-і + 53, - Ф/, /+1 + 54,7 Ф/-і, /4і- І

Рекуррентные соотношения для коэффициентов ах—аъ, 51^ — 54^ и Zlj приведены в приложении.

Произведение Ш содержит ненулевые элементы также на дополнительных диагоналях, которых нет в А (см. рис. 1.). Обозначим матрицу с этими ненулевыми диагоналями 5'. Ее элементы являются коэффициентами в выражении Ц]Ф при неизвестных Ф/+2,/-ъ Ф/-2. /, Ф<+2, /, Ф/-2, /+ь ЧТО легко увидеть, выразив Zi} из второго уравнения (7) и подставив соответствующие значения в первое уравнение. Сравнение результата данной подстановки с уравнением (3) показывает, что исходный девятиточечный шаблон расширился. Разложив получившиеся неизвестные в ряд Тейлора с точностью О (Дл: + Ду) относительно ближайших неизвестных Ф с ненулевыми коэффициентами матрицы А, которые входят в уравнение (3), получим У Ф = 5 Ф + О (Дх + Ду).

Положив Р = 3' — а5, где а — итерационный параметр, 0<а<;1, будем иметь следующую систему для нахождения элементов матриц X и и:

Вместо системы (6) целесообразнее решать уравнения относительно приращения ДФ —Фп+1— Ф”:

2. Уравнения Навье—Стокса для вязкой несжимаемой жидкости в переменных вихрь о» — функция тока Ч7 имеют вид

где М = фу) Ъ = — 'рх.

Численное решение этих уравнений приближенно представим с помощью сплайнов

Нелинейные конвективные члены квазилинеаризуем, сохранив тем самым члены с производными от функции 1|) на (л+1)-м временном слое. Удовлетворяя эти уравнения в узлах сетки и используя выражения (11) для со и о|), получим систему линейных алгебраических уравнений относительно коэффициентов ац, р^- следующего вида:

афп.

(8)

(9)

(10)

ш (X, у) ~ 2 X В‘в) ^

;=-1 /=-1 Л/ + 1 к+1

(10

/=-1 /=-1 ) Воспользуемся неявной схемой

I 1 т =— 1

(12)

3—«Ученые записки» № 5

33

Введя векторы Фи — (аи, ру)г, С}и = (С1ц, ($;)т и матрицы коэффициентов

систему (12) можно записать в матричном виде аналогично (3)

/=—1т-—1

и решать ее сильно неявным попеременно-треугольным методом, как было описано выше. Двухшаговый алгоритм дает решение

Рекуррентные соотношения для вычисления коэффициентов приведены в приложении.

В данном алгоритме требуется 18 вспомогательных двумерных массивов коэффициентов — в два раза больше, чем для сильно неявного метода с аппроксимацией с помощью центральных разностей [6]. Однако благодаря высокой точности приближения с помощью сплайнов общее количество расчетных точек можно уменьшить в четыре раза. В результате общий объем памяти сокращается, стационарное решение получается за небольшое число итераций, а время счета невелико.

3. В качестве примера решалось уравнение Лапласа в квадрате 0<х, у<1 с заданными граничными значениями функции. Для полного определения сплайна необходимо также второе краевое условие, например задание нормальной второй производной, которую можно выразить из уравнения. Так, при у = О

Из этих уравнений исключаются коэффициенты Рі.—і и таким образом получаются замыкающие уравнения для системы (2). Начальное приближение полагалось нулевым.

Скорость сходимости зависит от значений итерационного параметра а. Наименьшая скорость получается при постоянном а. На рис. 2 показана сходимость численного решения и к точному решению ы0, шах «о=1. При а=0 на равномерной сетке с числом 15X15 численное решение через 18—25 итераций отличается от точного на 10“3—10~4; на сетке 30X30 — через 70—100 итераций. При циклическом равномерном уменьшении с итерациями величины а от ад=1 до аі = 0 сходимость в несколько раз быстрее: для сетки 15X15 точность 10_3—10~4 достигается за 5—9 итераций, на сетке 30X30 — за 20—30 итераций.

Благодаря лучшей аппроксимации с помощью сплайнов скорость сходимости численного решения к точному оказывается значительно выше, чем с использованием той же схемы с центральными разностями (рис. 2). При а = 0 на сетке 30X30 указанная выше точность конеч-но-разностного решения достигается за 160—220 итераций.

А}1 А?

, I— 1, 2, . . . , 9, 6, р=1, 2,

У У, Аі-^Зт+5 Фі+1, І + Ш = (2ц

(13)

2 2 ^+і,'пВі+і(хі)Вт(0) = иі о,

1~— 1 т=—1 1 1

(14)

^ ^ Рі+', т &І+1 (Хі) Вт (0) — (Чуу)і 0.

I = — 1 т=*~ 1

----------В-сплайны,-------конечные разности;

1 — 4 — сетка 30X30; 5, 6-15x15; 1, 3, 5 — а = 0;

2, 4,6 — равномерное распределение , ав, . . . , а,,

0 < < 1

Рис. 2

Используя сплайны, можно получать приближенное решение той же точности, что и с конечными разностями, но на сетке с числом узлов, по каждому направлению в два раза меньшим [14]. Поэтому, хотя время счета на одной итерации для схемысВ-сплайнамивдва-три

раза больше, чем с конечными разностями, производя вычисления на

более крупной сетке, можно сократить общее время счета в 13 раз. На вычислительной машине с быстродействием 1 млн. операций в секунду на сетке 30X30 на одну итерацию затрачивается около 1 с.

4. Уравнения Навье—Стокса решались в задаче о расчете течения в квадратной каверне. Жидкость приводится в движение за счет вязкости верхней стенкой, которая перемещается со скоростью и= 1. На остальных стенках u = v = 0. Таким образом, для функции -ф можно записать граничные условия вида (14). Значения вихря и функции тока на стенках в сильно неявном методе связываются точно с помощью уравнения (10). Например, на нижней стенке

— (^уу)г о»

откуда

2 X [ <#*! » В1+1 (*,) Вт (У0) - р?#. т В,+, (*г) В"т (у0) ] = 0.

г=-1т=—1

Второе условие для исключения величин а,-1 получается из уравнения для со и условия прилипания.

В начальный момент принималось £0=113 = 0.

Расчеты проводились на сетках 14X14 и 28X28 для чисел Не=10, 100, 400 при ,М=1, 10, 106. При Ие —400 первые 10 шагов делались с условием проскальзывания на движущейся стенке. Полученные решения близки к известным (см., например, [6, 14, 17]). При больших решение получалось менее точным. Установление решения (| Асо/со | < <10~3) для Ие=100 на равномерной сетке 14X14 достигалось за 19 итераций, для Ие = 400 на неравномерной сетке 28x28 с минимальным шагом около стенки 0,02 и М=1—за 49, — значительно быстрее, чем

другими методами [1, 6, 14, 15, 16, 17], в которых требуются сотни итераций. Время счета составляло 0,5—5 мин. Сокращение числа итераций на порядок происходит благодаря сильной неявности метода, связи со и г|з на границах и высокой точности аппроксимации функций сплайнами.

ПРИЛОЖЕНИЕ

1. Рекуррентные формулы для решения систем (3), (8):

а1~Аи ai = Ai-j-alSli-h j-i,

«з = (1 — aSl/+i, j-1)-1 (Л3 + a2 51 г, /— i), а4 = Л4-)-а1 (S3i_i, j— i + a S2,_i. j—i) + a252/, y_i, аь = Л5 + 54/_i, j—\ 4- a2 53/, /_i -j- o3 52,-_|_i, y_i + #4 51/_i, /,

51 ij = as1 [Л 6 + «2 54/, /—i + a3 (53;+i, ,-_i + a 54/+i, /_i) ],

S2tj = «5 [Л, a4 (53;_i,a 52,_i, /) ],

53i;- = as1 (Л8 -f 54/—i, j), S4ij = a51 Л9,

Zij~~ &Ъ (Gij <X\ Z/—1, j—j Z/, j—1 Zi + i' j—\ <z4 Z<_i, у).

2. Решение систем (12), (13) для уравнений производится с помощью приведенных выше формул, в которых коэффициенты теперь являются матрицами:

Л^Л?-'], 1=1, 2, ..., 9, a, = la?*'], / = 1, 2, . . . , 5,

51/;- = [51?П , . . . , S4i; = [S4?/], k, р= 1, 2,

Zij=(Zh> Z2ij)T> {G, = G}h Gl)T.

ЛИТЕРАТУРА

1. Burgraff О. R. Analytical and numerical studies of the structure of steady separated flows. — J. Fluid Mech., 1966, vol. 24, p. 1.

2. Дайковский А. Г., Полежаев В. И., Федосеев А. И.

О расчете граничных условий для нестационарных уравнений Навье—

Стокса в переменных «вихря, функции тока». — В сб.: Численные методы механики сплошной среды, 1979, т. 10, № 2.

3. Владимирова Н. Н., Кузнецов Б. Г., Яненко Н. Н. Численные расчеты симметричного обтекания пластинки плоским потоком вязкой несжимаемой жидкости. — В сб.: Некоторые вопросы вычисл. и прикл. матем. Новосибирск, СО АН СССР, 1966.

4. С h о г i n A. J. A numerical method for solving incompressible viscous flow problems. — J. Comput. Phys., 1967, vol. 2, N 1.

5. Белоцерковский О. М., Гущин В. А., Щенников В. В.

Метод расщепления в применении к решению задач динамики вязкой несжимаемой жидкости. Ж- вычисл. матем. и матем. физ., 1975, т. 15, № 1.

6. R u b i n S. G., Khosla Р. К. iNavier—Stokes calculations with a coupled strongly implicit method. Part 1: Finite-difference solutions.—

AIAA Paper, N 0011, 1979.

7. Б у л e e в H. И. Численный метод решения двумерных и трехмерных уравнений диффузии. — Матем. сб., 1960, т.! 51, № 2.

8. Б у л е е в Н. И. Метод неполной факторизации для решения двумерных и трехмерных уравнений типа диффузии. — Ж. вычисл. матем. и математ. физ., 1970, т. 10, № 4.

9. Stone H. L. Iterative solution of implicit approximations of multidimensional partial equations. — SIAM J. Numer. Anal., 1968, vol. 5.

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

11. Марчук Г. И. Методы вычислительной математики. — М.: Наука, 1980.

12. Lin С. L., Pepper D. W., Lee S. С. Numerical methods for

separated flow solutions around a circular cylinder. — AIAA J., 1976,

vol. 14, N 7.

13. Завьялов Ю. С., Квасов Б. И., Мирошниченко В. Л. Методы сплайн-функций. — М.: Наука, 1980.

14. Rubin S. G., Graves R. A. Viscous flow solutions with a cubic spline approximation. — Int. J. Comp, and Fluids, 1975, vol. 3, N 1.

15. Rubin S. G., Khosla P. K. Navier—Stokes calculations with a coupled strongly implicit method. Part II: spline differred-corrector solutions. Lecture Notes in Math., N 771, Berlin efc., Springer—Verbog, 1980.

16. Казаков В. А. Применение сплайнов к численному решению уравнений Навье—Стокса при больших числах Рейнольдса. — Ученые записки ЦАГИ, 1982, т. XIII, № 5.

17. Казаков В. А. Расчет течения вязкой жидкости с помощью В-сплайнов. — Ученые записки ЦАГИ, 1983, т. XIV, № 4.

Рукопись поступила 22/111 1983