Применение сплайнов к численному решению уравнений Навье-Стокса при больших числах Рейнольдса Текст научной статьи по специальности «Физика»

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

Том XIII 198 2 №5

УДК 532.5.032

ПРИМЕНЕНИЕ СПЛАЙНОВ К ЧИСЛЕННОМУ РЕШЕНИЮ УРАВНЕНИЙ НАВЬЕ—СТОКСА ПРИ БОЛЬШИХ ЧИСЛАХ

РЕЙНОЛЬДСА

В. А. Казаков

На примере численного решения уравнений Навье—Стокса в задаче об обтекании плоской полубесконечной пластины вязкой несжимаемой жидкостью при больших числах Рейнольдса (1(Р—105) исследуется возможность применения сплайнов для аппроксимации функций с большими градиентами при численном решении нелинейных задач.

Во многих задачах, описывающихся нелинейными уравнениями, встречаются области с большими градиентами функций. Это обстоятельство приводит к значительным трудностям при решении задач-Узкие области с большими градиентами функций образуются, например, в поле течения вязкой жидкости при больших числах Рейнольдса. Точность численного решения при этом зависит от мелкости расчетной сетки и, следовательно, числа узлов в зоне больших градиентов. Для того чтобы численное решение правильно описывало течение, величина шага сетки в этих зонах должна быть достаточно мала. Вместе с тем общее число узлов сетки, исходя из возможностей вычислительных машин, не должно быть слишком велико. Чтобы удовлетворить этим требованиям, в ряде случаев используются преобразования координат, растягивающие узкие зоны. Эти преобразования можно производить либо заранее, если расположение области с большими градиентами приблизительно известно (например, [1, 2]), либо в процессе численного решения, когда сгущение сетки в исходных физических координатах связано с появлением в этой области повышенных градиентов функций, [3, 4].

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

Представляет интерес рассмотреть возможность их применения и при численном решении более сложных уравнений. В настоящей работе численно решаются уравнения Навье—Стокса в задаче об обтекании плоской полубесконечной пластины вязкой несжимаемой жидкостью при больших числах Рейнольдса (от 103 до Ю5). '

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

Re

(1)

(2>

скорости в направлении х и у, Re — число

Рейнольдса, соответствующее единичной длине.

Численные решения будем искать в ограниченной расчетной области л0-<х<;1, 0 = у0 <; у у* (рис. 1). Пластину расположим вдоль положительной оси х с нулевой передней кромкой.Введем сетку (Л^+ 1) X (ЛГ+ 1)( *0<*1 • • .<Хн = 1, 0 = _У0<_У! . . . <ук так, чтобы один из узлов хк совпадал с началом координат. Величины х0, Ук выбираются из условия малости их влияния на численное решение.

M,j

х а

X

•■к

Рис. I

При больших числах Re около пластины образуются области с большими градиентами параметров течения, поэтому используется неравномерная сетка со сгущением узлов в этих областях. Около носика пластины задается минимальный шаг сетки остальные узлы Xj, у} определяются по геометрическим прогрессиям. Их показатели в общем случае различны по разным направлениям от носика пластины и вычисляются по значениям Amin, х0, XNt ук,

Стационарное решение задачи будем получать как предельное установившееся решение уравнений (1), (2) методом переменных направлений, описанным в работе [9], с применением кубических сплайнов для приближения функций и их производных. Обозначим сплайн-аппроксимацию функции / в точке у}) fj}, ее первых

производных fx и fy — m{j и вторых производных fXI, /уу — M{j, Ly- Численную схему для функции ® запишем следующим образом. Для первого шага

.f/w _.«,?/ +‘^[- (Й»1л;)"+ (mW(0"+'»+i(£r/)"+l;5 +

для второго—

+-±{жг,)л+1

Ке

Здесь верхние индексы п обозначают значения функции в момент времени г = п.М, где М — приращение по времени на каждом шаге п.

Для функции тока ф используется итерационная схема с установлением по фиктивному времени х = ^Дт,

Для первого шага

Как показано в работе [9], приведенные схемы на равномерной сетке абсолютно устойчивы.

Определим граничные условия в задаче. На оси х ф = 0, v = 0; на пластине a— и£0 = const (ut0 = 0 в случае отсутствия проскальзывания). Вне пластины из условия симметрии <о(0 = 0, при этом из уравнения движения следует, что и = L>?o& 0. На пластине

граничное условие для ш получается из уравнения (2) для функции тока и условия 4 = 0 вдоль оси х\

где = —у0. . ;

Выразим отсюда L% и подставим в (7), учитывая, что /f0=(^y)i0= = Mi0 и ’т'го — О; получим

Это граничное условие имеет второй порядок точности относительно шага сетки, поскольку вторая производная ф аппрокси1 мируется кубическим сплайном на неравномерной сетке со вторым порядком, производная фу—с третьим, а сама функция—с четвертым.

Величина ш/0, вычисленная по формуле (8) с условием прилипания на пластине, принимается в качестве а>?0+1, значение <вй)+1''2 определяется интерполяцией:

(5)

для второго —

(G)

0);о — t-io •

(7)

Для кубического сплайна имеет место соотношение

(8)

Л-Н/2 1 /о л-4-1 , /? п П—1\

<»;о ' = — ( бщо' + 6o>i0 — ©го ) . о

Время счета при этом сокращается.

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

33

Для полного определения сплайна требуется задать второе граничное условие. Вне пластины = 0, второе условие на пластине получается с помощью уравнения (3), из которого выражается вторая производная (ЛГо)л+! 2. Для функции 41 на оси х ^ = вне пластины Ь% — 0.

На верхней границе области зададим условие безвихревого течения о).д ==0. В качестве второго условия для сплайна положим /."*• = 0, что, как показывают расчеты, выполняется с большой точностью. Для функции <р можно задавать жесткие условия невозмущенного потока = Ук, $к = Щк—1 или мягкие условия протекания, когда значения ^ могут меняться. В случае задания жестких условий может проявиться эффект запирания: поджатие пограничного слоя, увеличение скорости и сопротивления трения пластины. Поэтому для его устранения необходимо удалять верхнюю границу, что существенно увеличивает время счета и, кроме того, может потребовать введения дополнительных узлов сетки.

При задании мягких условий используется уравнение (5), в качестве второго условия можно взять £/*■==/,?*-_ 1. Условие £^=0 действует для данной схемы, как жесткое.

На левой границе полагаются условия набегающего потока: <Ь=У, й=1, /яФ=—= 0, ш = 0. Второе условие для функции <в получается из уравнения движения (1):

Не

на оси х Жоо = 0.

На правой границе расчетной области используются уравнения (4), (6) и условия М£/ = М%и = 0. Применение в качестве граничных условий при Хк = 1 асимптотических выражений — и>Ц2х) и

Лх $/'(2х) [И, 12] при начальных данных, значительно отличающихся от решения для пограничного слоя на пластине, ведет к неустойчивости. В то же время установившееся решение, полученное С условиями Мяу = Жа'/= 0, удовлетворяет с той же точностью и асимптотическому условию (не требуется делать ни одного дополнительного шага).

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

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

2. В уравнение (3) на первом шаге для функции <и входят неизвестные величины , (/"/)"+12 и (£“/)в+1'2 на временном слое

п + 1/2. Производные удовлетворяют следующим уравнениям [5]:

к} + 2 (к} + £,ч,)Ч- + */+1 ^У+1 = 6 (■(9)

где к]^=у]—у!-1.

Кроме того, для кубического сплайна имеют место соотношения

/+1

АУ+1 >

4- £*7+1) +

<7+1

д N 1 ^Ч-Г^/ I Ь

6 я/+1

(10)

(П)

Исключив из выражений (3), (9)—(11) величины <»,7+1/2 и (///)п+1,2) получим для каждого г, г = 1, . . . , Л;, трехдиагональную систему уравнений относительно (£рл+12

(!?,■+,)л+,;24 Ь} (Ц)п+и2 + с, (1“ж )л+1/2 = <*у, ] == 1, . ■ - , К- 1 (12)

с коэффициентами

—"5/12*'+ ^ (га> + “</-■)]

°»7+1

к 1+1

зд/

2*,+,+ ^ 2 г ф

+ тч+1

Ке6у+1

7+1

где

и правой частью 6

й,

В<7+1 *У+

г И

. п , Д £

/+1 — ^7 + —

я/+1 «+1)я

Щ)*+ &№,)*-

6 (л л . Д* Г 1 ^ ЛЛШ ЛЛ® \Я Г

I Ш;у ~ °>г/_1 ~2 ’ [ Кё- (' <7' ^>7-0 Т

Для замыкания этой системы используются условия на нижней и верхней границах. Получившаяся система решается методом прогонки.

На равномерной сетке с шагом Л при сеточном числе Кеь = —Revfl<C 2 условие диагонального преобладания матрицы выполняется для любых М, При 1?ел > 2 для этого необходимо выполнение условия

'(■-■У-

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

На втором шаге для ш, а также для функции ф получаются системы, аналогичные (12).

Установление стационарного решения проверялось с помощью

критерия тах [а>Г/и — <о£ | < е, в большинстве расчетов величина

и

э = 10—4. При меньших г коэффициент трения на пластине менялся лишь в четвертом—пятом знаках, с уменьшением е с 10~3 до 10~4 изменение в распределении коэффициента трения составляло десятые доли процента, а время расчета увеличивалось в 1,5 раза.

3. На рис. 2 изображена полученная по результатам расчетов при Ие = 103 картина линий тока течения в окрестности носика пластины. Хорошо видно выталкивающее действие заторможенной жидкости около поверхности пластины, в результате чего линии тока отклоняются вверх. На рис. 3 приводится распределение коэффициента трения на пластине:

т (-*) 2 , т С, = —— —------ ш(.Г, 0),

1 р ицг Ие к ’

полученного из расчетов для Не=Ю3 на сетках 21X14 = 9,

лг0 = — 0,3, у/г = 1, Лщщ = 0,02) и 61X31 (г0 = 20, я:0—— 1, у* = 2,

Лтщ = 0,015) с коэффициентами неравномерности (отношениями соседних шагов), изменявшимися в пределах от 1,02 до 1,05 для мелкой сетки и от 1,12 до1,28 — для грубой (при этом наибольший шаг сетки на пластине А,= 0,039 ч- 0,23). Как видно из графика, значения местных коэффициентов трения на большей части пластины отличаются менее чем на 2 — 3%. При х > 0,4 результаты расчетов близки к решению Блазиуса УИе су = 0,664/1/*. В точке лдг из-за недостаточно точных граничных условий значения коэффициента трения на 4—7% меньше. Коэффициенты сопротивления трения участка пластины единичной длины, вычисленные на различных сетках, отличаются не более чем на 2% и на 4—5% больше соответствующего значения Блазиуса.

Для сравнения на рис. 3 приводятся коэффициенты трения из работ [1, 2], в которых численное решение получалось итерационным разностным методом на сетках 33 X 129 и 120 x 40 в параболической системе координат, оптимальной для данной задачи. В качестве начального приближения н условий на выходной границе расчетной области принималось решение Блазиуса. Особенность для функции вихря в носнке пластины учитывалась путем введения новой переменной вида (л:2-)-_уг) ш, в результате чего можно было уменьшить неточности, возникающие при замене производных конечными разностями. При х >0,1 коэффициенты трения близки к значениям Блазиуса. В работе [14] вводятся безразмерные переменные, соответствующие пограничному слою, где величины имеют порядок: у~1/]/гЯе, г'~1/уг1?е, ф~1/|/Ке. Таким образом,

получаются уравнения с коэффициентом 1/Ке лишь при второй производной по х. Результаты расчетов для Ие=56 000 на сетке 77X75 близки к решению Блазиуса и на рис, 3 не показаны. Как и в настоящей работе, при использовании однородной аппроксимации величина коэффициента треиия в носике пластины (функция ш) конечна и зависит от местного шага сетки ктш- чем меньше шаг, тем больше С; и со (см. рис. 3). При одинаковых Атш значения с} в носике пластины, вычисленные на различных сетках, отличаются не более чем на 1%. В работе [15] итерационным методом с применением аппроксимации на 13 точках решалось одно бигармоничес-кое уравнение для функции тока на сетке 31 X 15 для Ие = 10250. Сопротивление пластины получилось на 10% меньше значения Блазиуса.

Коэффициент сопротивления, вычисленный на сетке 61X31 для ряда значений х0, ук{—2<л:0<—0,3; 0,5<_у^<2), изменялся на 1,4% и составлял 0,0434—0,0440. При жестких условиях непротекания на верхней границе расчетной области проявляется эффект запирания: при Улг = 1 сопротивление пластины на 6—7% больше, чем при мягких граничных условиях, и на 10—12% больше значения для пограничного слоя Блазиуса; величина продольной составляющей скорости и вне пограничного слоя на 6% превышает скорость набегающего потока (в работе [15] для Ке = 10250 — на 16%). При удалении верхней границы эффект запирания исчезает, максимальное значение и равняется 1,018 на грубой сетке и 1,007 — на мягкой. При мягких верхних условиях функция тока на линии у=^ук изменяется на такую же величину, что и при жестких условиях на этом же уровне, но при достаточном удалении верхней границы. Таким образом, оказалось возможным рассчитывать течение практически с той же точностью для различных^ вплоть до _уд- = 0,25,

что лишь в 1,5 раза превышает толщину пограничного слоя в конце пластины, вычисленную по Блазиусу [13]. Время расчета при этом уменьшается в несколько раз.

На рис. 4 приводятся профили продольной составляющей скорости и при некоторых х, полученные на сетке 21X14. Значения и, как функции автомодельной переменной пограничного слоя, в различных сечениях пластины на большей ее части близки к профилю Блазиуса (рис. 5, значения скорости и отнесены к ее наибольшей в данном сечении величине ие).

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

работе [12], целесообразно выделить область с большими градиентами и рассчитывать ее отдельно, сращивая получаемое решение с решениями в других областях. В качестве примера в настоящей работе отдельно (без решения в других областях) рассчитывалось течение при Яе = 105 около носика пластины в области —0,001, Хл/=0,004, _удт = 0,01 с АШ1П—10~4. Вследствие подобия решений полученные результаты аналогичны приведенным для 1?е=103.

На вычислительной машине с, быстродействием около 1 млн. операций в секунду время расчета течения на сетке 21 X 14 при г = Ю~4 и ук = 1 составляло 4,2 минуты, а на сетке 61X31 —на порядок больше.

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

ЛИТЕРАТУРА

1. Van de Vooren A. I., Dijkstг a D. The Navier—Stokes solutions for laminar flow pas! semi—infinite flat plate. J. .Engrg. Math.", vol. 4, N 1, 1970,

2. Yoshizawa A. Laminar viscous flow past a semi—infinite flat plate. J. .Physical Soc‘. Japan, vol. 28, N 3, 1970.

3. Толстык А, И, Об исследовании течений вязкого сжимаемого газа при помощи полных уравнений Навье—Стокса. В сб. „Чнсл. методы механики сплошной среды". Новосибирск, т. 6, № 4, 1975.

4. Ковен я В. М., Яненко Н. Н. Разностная схема на подвижных сетках для решения уравнений вязкого газа. „Ж. вычисл. матем. и матем. физ.', т. 19, № 1, 1979.

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

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

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

8. Гребенников А. И. Методы сплайнов в численном анализе. МГУ, 1979.

9. RubinS. G., Graves R. A. Viscous flow solutions with a cubic spline approximation- .Int. J. Comp, and Fluids', vol. 3, N 1, 1975.

10. Казаков В. А. Применение сплайнов к численному решению нелинейных задач с большими градиентами функций. „Ученые записки ЦАГИ“, т. IX, № 5, 1978.

11. Дородницын А. А., Меллер Н. А. О некоторых подходах к решению стационарных уравнений Навье— Стокса. „Ж. вы-чнсл. матем. и матем. фнз.', т. 8, № 2, 1968.

12. D о го d п 1 с у п A. A. Review of methods for solving the Navler— Stokes eqations. Proc. 3-d. Int. Conf. Num. Meth, Fluid Mech. Lecture Notes in Physics. Berlin—Heidelberg—New York, Springer—Verlag, vol. 18, p. 1, 1969.

13. ШлттингГ. Теория пограничного слон. М., .Наука*, 1969.

14. L о е г S. Von. Eine numerische Methode zur LOsung der Navier— Stokesschen Gleichungen fur die zweidimensionale, incompressible, staiio-nare Stromung ISngs einer diinnen Platte. Ingeneur— Archiv, vol. 41, N 1, 1971.

15. Walter К. Т., Larsen P. S. The FON method for the steady two-dimensional Navler—Stokes equations. .Int. J. Comp, and Fluids", vol. 9, N 3, 1981.

Рукопись поступила 121VII 1979 Переработанный вариант поступил ЩХ/f 1981