Расчет течения вязкой жидкости с помощью В-сплайнов Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ ЦАГИ Том XIV 1983

№ 4

УДК 517.9:532.5:533.7

РАСЧЕТ ТЕЧЕНИЯ ВЯЗКОЙ ЖИДКОСТИ С ПОМОЩЬЮ

5-СПЛАЙНОВ

В. А. Казаков

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

В работе [1] для численного решения уравнений Навье—Стокса, описывающих течение вязкой несжимаемой жидкости, было предложено использовать сплайны [2—4]. В дальнейшем метод сплайн-коллокации был развит в работах [5—8] при построении схем повышенного порядка точности. Благодаря хорошим аппрок-симационным свойствам сплайнов для расчета течений требуется меньшее число узлов сетки и, следовательно, меньшее время счета. В настоящей работе для расчета течений вязкой жидкости предлагается использовать кубические 5-сплайны, которые позволяют еще более повысить эффективность решения уравнений Навье — Стокса: упрощаются алгебраические уравнения для численной схемы, уменьшается требуемая машинная память и время счета. В качестве примера рассчитывается течение вязкой жидкости в каверне.

1. Уравнения Навье —Стокса, описынающие плоское течение несжимаемой жидкости, можно записать следующим образом:

и 'а>х + V , (1)

Ьх + Ьу = ш> (2)

здесь и = фу, V=—'Ьх — составляющие вектора скорости в направлениях х и у, Ке — число Рейнольдса.

Пусть имеется некоторая расчетная область с сеткой х0< х1 < < . , . < х.\г, У0<^1< • ■ • <СУк- Для построения Б-сплайнов дополним ее произвольными узлами по оси х\ х^з < х_2 < Х—\ < х0, Хк<СХд>+1 <Сх;\."+2 <СХу+з и так же —по оси у. 5-сплайны являются финитными функциями, отличными от нуля только на некотором

интервале. Будем обозначать кубический 5-сплайн 5, (х)—по среднему узлу его интервала-носителя (х;_2, хг+2), I — — 1, 0, . . . , N + 1. Аналогично определим (£), 1, 0, ... , К+ 1.

Кубический сплайн 5 (х) можно представить [3] в виде

Л7+1

$(*)=£& В1(х),

1=-1

где — постоянные коэффициенты.

Вследствие финитности 5-сплайнов в каждой точке сетки отличны от нуля только три 5-сплайна, поэтому для 5 (х) и его производных порядка &

5(А) (*,) = £ Р;+т М*1» и,), & = 0, 1, 2.

т =—1

При использовании 5-сплайнов для решения дифференциальных уравнений методом коллокации требуется знать значения 5-сплайнов и их производных только в точках коллокации, а поскольку в процессе счета эти величины используются многократно, целесообразно вычислить их один раз и образовать специальные массивы. Если в качестве точек коллокации берутся узлы сетки, для каждой функции В1 (х), а также их производных, вычисляются значения только в точках Х;_1, хь Х1+\. В результате составляются массивы вида Ь1т = В1{Х1+т), г = — 1, О,..., N-\- 1, т, = — 1, 0, 1, размерности (ЛГ + 3)><3. В случае равномерной сетки вместо массивов получаются постоянные величины.

Для решения уравнений (1), (2) воспользуемся схемой переменных направлений [1]. Функции миф будем аппроксимировать кубическими сплайнами двух переменных, которые можно представить в виде „тензорного произведения11 одномерных сплайнов [3]. Для функции (о

ЛГ+1 К+1

«>(*, ,У)~Е 2 ГЬВ1{х)В1(у). (3)

I = -!/=-!

Значения сплайна и его производных в точках (хь, у}) обозначим со(-у, (Шд.)г-у, (ш*х)г./> (шуу);7-

На первом шаге для уравнения (1)

на втором

фЯ + 1 ___ ф/г+1/2 I

и ч' 1

Обозначим в (3)

ЛГ+1

(*) = 2 в‘ач = (х‘)>

;*-1

тогда

К+1

>(х,

,=-1

К+\

і

И В точке {хь _Уу)

«(*„ ^)«0>у = 2 %1кВ к (Уу) = X “//+«£/+» (-У/)'

&=—1 т = — 1

Таким образом, для временного слоя я+1/2 можно записать 1 1

«?/и- 2 *гй?а№.<ул <»л+,,-2 «?дтв;+.ад,

т =—1

ГП-—1

1

(со )Л + 1/2: ' УУПІ

V п+1 2

Т . .

т = — 1

я<у +ш В/+т (У у).

(6)

При вычислении Ш7+1 на втором шаге будем использовать представления через 5-сплайны по л:

ил + 1

і І

1

V

ач-ш, у5(+т(л:г),

I

(“лг)”}'"1 ^ ^ а/+т, І Ві + т (X;);

т = — 1 1

= 2 яг+т, / 5, + т (х;).

т = — 1 1

(7)

т=—1

В вычислительной программе достаточно одного массива коэффициентов {«/у}. Индексами п, п ■+- 1/2, /г -)- 1 будем обозначать их значения на соответствующих временных слоях.

На каждом шаге вместо двумерного массива {а<7} целесообразно использовать одномерный. На первом шаге для каждого фиксированного Х1 будем полагать в (6) Ру — а?^1/2, у = — 1, О, . . . , А’ + 1. Подставив эти выражения в уравнение (4), получим для каждого столбца х1 трехдиагональную систему уравнений

аі Р/—і Ь] Ру + сі Р;+і — 1 — 1, 2, . . . , /С — 1,

(8)

где

аі = 5/-і (Уу) +

Дг1

С;

Д£

;- = 5ж(з/у) 4 2

'' Не

Я/-1 (Уу)

0>у)

■Ві+ііу^

Ие

*у

= со'.1. + — Г—1—(со„)?. — И" (а»,)?.

Производные (о)Л)« и (шхх)?}, входящие в правые части уравнений, определяются по известным коэффициентам {а"} на п-м временном слое по формулам (7). После вычисления величин на столбце х1 можно полагать а._1 у = Ру, *где Ру— решение системы

(8) на предыдущем столбце л/_ь В результате такого обмена достигается в общей сложности большая экономия машинной памяти.

Решив систему (8) для каждого столбца хіг получим массив коэффициентов (а”+1/2), которые определяют решение на слое

Таким же образом на втором шаге находятся коэффициенты {а"+1}, определяющие по формулам (7) решение на новом временном слое.

На равномерной сетке данная схема в линейном приближении абсолютно устойчива.

После вычисления значений со?*1 на (га+ 1)-м слое решается уравнение (2) для функции Ф. Численное решение ищется методом установления по фиктивному времени т. Используется двухшаговая схема переменных направлений:

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

ДЛЯ ф.

Получив установившееся решение {фу} непосредственно по соответствующим коэффициентам сплайна, вычисляются значения компонента скорости — (ФД,-. Для определения величин

и1} = (фу),7 по значениям ф/;- строятся одномерные сплайны по столбцам для каждого х{.

Приведенные схемы аппроксимируют исходные уравнения со вторым порядком на неравномерной сетке.

2. Рассмотрим задачу о расчете течения вязкой жидкости в каверне (рис. 1). Жидкость приводится в движение верхней стенкой, которая перемещается в направлении оси д; со скоростью

л + 1/2.

(10)

(9)

0 0,2 0,4 0,6

Рис. 1

0,8

и Определим л;0 = у0 = 0, хм=\\ для квадратной каверны ук=\. Число Рейнольдса, таким образом, определяется вязкостью: Ие = 1^.

На границах каверны положим ф = 0, на неподвижных стенках « = ^ = 0, на движущейся стенке и — '1>у= 1. В качестве начальных данных зададим состояние покоя (ф = ш = 0).

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

Остановимся подробнее на постановке граничных условий для сплайнов и выводе замыкающих уравнений для систем уравнений относительно коэффициентов 5-сплайнов. Для функции ш на нижней границе из уравнения (2) следует (фуу)/0 = “/о- Отсюда, используя известное соотношение для кубического сплайна

(V.. =6 ~ <«<»- т- «Л >•

п о п0 п0

где Ь0 = у1 —_у0, и учитывая, что фг 0 — (фу)г 0 = 0 и (Ф,)м = «/ и получим

(П)

Это условие имеет второй порядок точности, так как вторая производная 'Ьуу аппроксимируется кубическим сплайном со вторым порядком. Непосредственное использование разложения функции ф в ряд Тейлора дает граничное условие для а>го лишь первого порядка:

Фг 1 + О (Н0).

По формуле (11) вычисляются значения (««+1. величина находится интерполяцией

“Го1/2 = 4“ (“"о1 + Ьш?о — :) + 0 (д^2)- (12)

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

( Шуу)« + >/2 = 2Ие К+1/2 _ Шп0)/М - («^)»0.

Таким образом, для каждого х1 при / = 0 можно записать 1 1

2 ^пгВт(Уо) = ^Г12, 2 &Х(Уо)=(. “уу)?о1/2- (13)

т — — \ т— — 1

Исключив из этих уравнений величину получим уравнение

вида

Ро + Рх = ^о- (14)

На верхней движущейся границе таким же образом выводится второе замыкающее уравнение

ак $к-\ + Ьк $к = йк. (15)

В результате для каждого хг(г = 0, 1, . . . , И) получается трех диагональная система алгебраических уравнений (8), (14), (15) относительно величин Ъ:, ./ = 0, 1, К. Решив систему, из (13) определим коэффициент из аналогичных соотношений при

7 = К находится ^*-+1.

Так же получаются замыкающие уравнения на втором шаге. При выводе соответствующих уравнений для функции ф используются условия: на неподвижных стенках 4 = 0, и = V = 0, на движущейся— = 0, фуу = ш. В верхних угловых точках функция со имеет два значения, замыкающие уравнения выводятся из условия параболической экстраполяции.

Вследствие того что в начальный момент около движущейся стенки возникают значительные градиенты, устойчивость счета может ухудшиться, поэтому первые несколько шагов делались с условием проскальзывания: значения и!К за 20 шагов по времени увеличивались до 1.

Установление решения контролировалось с помощью критерия тах I <о"+1 — (о?. I <■ г, величина е = 10~3 -г-10—5.

I/

3. Задача о расчете течения вязкой жидкости в каверне решалась многими исследователями. В работе [1] для чисел Ие = 10 и 100 приводится сравнение результатов, полученных конечноразностными методами (для дивергентных и недивергентных уравнений) и методом сплайн-коллокации с использованием кубических сплайнов в кусочно-многочленном представлении. Показано, что применение сплайнов позволяет сократить число точек расчетной сетки в четыре и более раз (в два раза по каждой координате). В настоящей работе получены такие же результаты.

На рис. 2 приводятся профили составляющей скорости и при .* = 0,5 и в сечении, проходящем через центр вихря с координа-

Рис. 2. Данные работы [1] на сетках 15x15, 65x65 и [14] получены конечно-разностными методами

тами хв, ув (л:в = 0,63, ув = 0,75 при Ие=100). Расчеты выполнены на равномерных сетках 15X15 (Н = 0,071) и 29X29 (Л = 0,0357). Максимальное значение функции тока | ф |Шах = Р, 105, величина ш в средней точке движущейся стенки сото = 6,35 (в работе [1] | ф |шах=0,1043, = 6,687). Эти данные совпадают в пределах 1—5% с результатами, полученными различными авторами на более мелких сетках (в работах [9, 10] имеются подробные сравнения).

Для чисел Ие = 103 устойчивого счета на сетке 15X15 достичь не удается, на неравномерной сетке 29X29 (Нт1п = 0,01 -ь 0,02, Лтах ~ 0,06 -г- 0,08) при достижении установления решения в распределении св на движущейся стенке имеются некоторые колебания, которые исчезают для более мелкой сетки 35X35 (Нтт = 0,015, Лшах~0,04). В работах [6,9, 10] отмечается, что для Ие^400,устойчивый счет на относительно грубых сетках возможен при использовании уравнений в дивергентном виде, В таблице приводится

Работы Сетка '■Рв (Ов (хв, Ув)

[П] 41X41 0,1071 1,891 (0,55; 0,55)

[П] 51X51 0,1090 1,986

[П] 81X81 0,1132 2,080 (0,52; 0,58)

[6] 65X65 0,114 1,985 (0,52; 0,56)

[6] 17X17 0,115 1,828 (0,56; 0,56)

Настоящая работа 35X35 0,117 2,00 (0,525; 0,56)

сравнение некоторых результатов вычислений для Ке=103 с данными, полученными конечно-разностными методами [11,6] на мелких сетках и неявным методом сплайн-коллокации четвертого порядка для уравнений в дивергентном виде на сетке 17X17 [6].

На рис. 1 показана картина линий тока течения при Ке=103. В нижних углах каверны наблюдаются вторичные вихри. В центральной части области имеется отчетливо выраженное ядро течения с постоянной завихренностью (рис. 3), что соответствует модели Прандтля — Бэтчелора [12, 13]. Величина скорости | V | движения жидкости в этой области линейно зависит от расстояния до центра вращения (рис. 4). Из этого рисунка видно также, что с увеличением числа Ие размер ядра растет, а центр вихря смещается к середине каверны. Это отмечается также в работах [11, 14]. Линейная протяженность области, где величина <о отличается от значения в центре вихря % не более чем на 2%, при Ре= 108 составляет 0,4, на 5% —0,6.

На вычислительной машине с быстродействием около 1 млн. операций в секунду время расчета для Не=100 на сетке 15X15 (А = 0,071, Д* = 0,1, Дх = 0,01) составляет около 0,5 мин, для Не = 400 на сетке 29X29 (Лтт = 0,03, Д£ = 0,1, Д-с = 0,05) при в= 10—3+-10—4 около 5—7 мин, для Не= 103 на сетке 29X29 (/гш1п = = 0,02, Д^ = 0,1, Дт = 0,05) — 12—18 мин. В работе [6] решение с помощью кусочно-многочленных сплайнов для Ие = 400 на сетке 17X17 получалось за 4 мин на вычислительной машине с быстродействием в 3—8 раз большим.

Рис. 3

Рис. 4

Автор благодарит Вик. В. Сычева за полезные замечания при обсуждении результатов работы.

ЛИТЕРАТУРА

1. Rubin S. G., Graves R. A. Viscous flow solutions with

a cubic spline approximation. Comp, and Fluids, 3, N 1, 1975.

2. Алберг Д ж., Нильсон Э., Уолш Дж. Теория сплайнов и ее приложения. М., „Мир*, 1972.

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

4. Стечкин С. Б., Субботин Ю. Н. Сплайны в вычислительной математике. М., „Наука", 1476.

5. Rubin S. G., Khosla Р. К. Higher-order numerical solutions using cubic splines. „А1АА J.“, 14, N 7, 1976.

6. Rubin S. G., Khosla P. K. Polynomial interpolation methods for viscous flow calculations. BJ. Comp. Phys.“, 24, p. 1, 1977.

7. Rubin S. G.. Khosla P. K. A simplified spline solutions

procedure. Proc. 6-th Int. Conf. Num. Meth. Fluid Dyn. Lecture Notes in

Ph>sics. Berlin — Heidelberg — New York. Springer, N 90, 1979.

5— «Ученые записки ЦАГИ» № 4

65

8. Rubin S. G., K h o s 1 a P. K. Navier — Stokes calculations with

a coupled strongly implicit method. Part II: spline differed — corrector

solutions. Lecture Notes in Math.'Berlin — Heidelberg — New York, Springer, N 771, 1979.

9. Bontoux P., Forestier B., Roux B. Analysis of higher-

order methods for the numerical simulation of confined flows. Proc. 6-th. Int. Conf. Num. Meth. Fluid Dyn. Lecture Notes in Physics. Berlin — Heidelberg — New York. Springer, N 90. 1979.

10. Roux B., Bontoux P., Daube O., Phuoc Loc T. Optimization of Hermitian methods for Navier — Stokes equations in the vor-

ticity and stream-function formulation, Lecture Notes in Math. Berlin — Heidelberg — New York, Springer, N 771, 1979.

11. Shay \V. A. Development of a second order approximation for

the Navier — Stokes equations. „Comp, and Fluids", vol. 9, 1981.

12. Prandtle L. Ober Fliissigkeitsbewegung bei sehr kleiner Rei-bung Verhandlung, d. Ill Int. Math. Kongr. Heidelberg. 1904. Teubner, Leipzig, 1905.

13. Batchelor G. K. On steady laminar flow with closed streamlines at large Reynolds number. „J. Eluid Mech.“. 1, p. 2, 1956.

14. Burgraff O. R. Analytical and numerical studies of the structure of steady separated flows. „J. Fluid Mech.“, 24, p. 1, 1966.

Рукопись поступила 2611 1982