Расчет обтекания острова с использованием адаптивных сеток Текст научной статьи по специальности «Физика»

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

Том 8, № 2, 2003

РАСЧЕТ ОБТЕКАНИЯ ОСТРОВА С ИСПОЛЬЗОВАНИЕМ АДАПТИВНЫХ СЕТОК*

Г. С. Хлкимзянов, Н. Ю. ШокинА Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: khak@ict.nsc.ru, nina.shokina@ict.nsc.ru

Iterative algorithm for numerical modelling of stationary fluid flows with surface gravitational waves in river-beds with islands using adaptive grids is presented. Plane shallow water model is used taking into account bottom irregularity, frictional force and Coriolis force. The problem is formulated in new dependent variables, i.e. stream function and vorticity function. Modification of equidistribution method for curvilinear grid generation in multiply connected domains is given. Results of numerical modelling of flow around island in the river-bed with complicated shore line are presented.

Введение

Проблема моделирования установившихся и нестационарных течений в руслах рек и их устьях является актуальной [1-3] ввиду возможности использования результатов расчетов для прогнозирования размыва берегов и переноса загрязняющих веществ.

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

В настоящей статье методика работы [4] обобщена на случай многосвязных областей, что позволило применить ее для исследования обтекания островов со сложными очертаниями береговой линии. Сделаны акценты на отличия алгоритма для многосвязных областей от алгоритма, описанного в [4] для односвязных областей. В частности, приведена модификация метода эквираспределения [6] для построения криволинейной сетки в многосвязной области и разработан алгоритм для нахождения значения функции тока на контуре острова. Отметим, что здесь возникает та же проблема, что и при моделировании течений вязкой жидкости в многосвязной области, заключающаяся в необходимости определения на каждой итерации значений функции тока на внутренних контурах (см., например, [7, 8]).

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №00-01-00899, и Программы интеграционных фундаментальных исследований СО РАН, проект №1.

© Г. С. Хакимзянов, Н.Ю. Шокина, 2003.

1. Математическая модель

Для описания течения используется модель мелкой воды первого приближения для идеальной несжимаемой жидкости со свободной границей. Для простоты изложения будем далее предполагать, что имеется только один остров, так что область течения П — дву-связная область в плоскости декартовых координат x1Ox2. Внешнюю границу области П будем обозначать через Г = дП = Г! U Г0 U Г2, а для контура острова, лежащего внутри Г, введем обозначение Г^стр. Здесь через Г! обозначен вход в П, через Г2 — выход из П, через Г0 обозначена непроницаемая часть внешней границы Г. Будем предполагать, что граница острова — вертикальная непроницаемая стенка.

Математическая постановка задачи об установившемся протекании жидкости с поверхностными гравитационными волнами через область П заключается в определении вектора скорости u и функции п(х) — возвышения поверхности жидкости над невозмущенным уровнем, удовлетворяющих в П уравнениям неразрывности и движения

V(Hu) = 0, х = (x!,x2) G П, (1)

H 2 g

(Hu ■ V)u + gV — = gHVh — kH[ез x u] - -g-u|u|, (2)

2 C

краевым условиям на внешней границе

Hu

xerx

vi(x), u ■ n =0, Hu ■ n = v2(x) (3)

xero

хеГ2

и условию непротекания через границу острова

0. (4)

un

с€Г0стр

В этих формулах u = (u1,u2) — вектор скорости, ua (а = 1, 2) — его компоненты в направлении осей Oxa; V = (d/dx1, d/dx2); n — внешняя нормаль к границе Г U Г^стр, v 1 ■ n < 0 на входе, v2 > 0 на выходе; H = n + h — полная глубина; z = —h(x) — функция, описывающая дно; g — ускорение свободного падения; е3 = (0, 0, 1) — базисный вектор, направленный вдоль вертикальной оси Oz, перпендикулярной к плоскости x!Ox2; C — коэффициент Шези; k — параметр Кориолиса; плотность жидкости полагается равной единице.

Система уравнений (1), (2) записывается сначала в безразмерных переменных, а затем в новых зависимых переменных — функции тока ф и функции вихря ш [4].

Условие (4) приводит к следующему граничному условию для функции тока на контуре острова:

ф = фостр = const, (5)

хег0стр

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

дф дТ

= Hu n

с€Г0стр

(6)

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

0

ф на Г0 также определяются однозначно, а на участке Г0стр значение фостр не задано и не может быть определено заранее, до решения задачи. Эта константа будет находиться в ходе решения задачи вместе с решением внутри области.

С помощью взаимно-однозначного невырожденного отображения единичного квадрата Q, лежащего в плоскости координат д1 Од2, на область ^

= (д1 ,д2) (7)

система уравнений для ф и и записывается в новых независимых переменных да, а = 1, 2. При вычислении полной глубины используются уравнения движения в форме Громеки — Лэмба, записанные в криволинейных координатах [4].

2. Метод эквираспределения для многосвязных областей

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

Рис. 1. Река с островом.

Рис. 2. Область течения и изобаты поверхности дна.

Рис. 3. Вычислительная область.

Рис. 4. Расчетная сетка (81 х 51), а = 10.

строится на основе варианта двумерного метода эквираспределения для многосвязных областей.

Остановимся на модификациях ЕБ2-метода, разработанного ранее для односвязных областей [6]. Предполагая, что граница острова Г0°тр — криволинейный четырехугольник, будем считать, что она является образом при отображении (7) контура 70°тр прямоугольника, лежащего внутри единичного квадрата плоскости д1 Од2, при этом стороны этого прямоугольника, являющиеся прообразами дуг АВ, БС и АБ, ВС, параллельны осям Од1 и Од2 соответственно (рис. 3).

Таким образом, в рассматриваемом случае вычислительная область Q представляет собой единичный квадрат с вырезанным из него прямоугольником с контуром 70°тр. Область Q покрывается равномерной прямоугольной сеткой (^. Стороны прямоугольника 70°тр проходят по некоторым линиям этой сетки. Положение прямоугольника на сетке и длины его сторон определяются задаваемыми изначально числом узлов №°тр на сторонах АВ и БС, числом ^2°тр узлов на сторонах АБ и ВС, а также указанием узла цА, являющегося прообразом угловой точки А.

В отличие от односвязных областей теперь необходимо расставлять узлы не только на внешней границе Г, но и на контуре острова Г°стр. Для этого используется метод эквираспределения для плоских кривых [6] на каждой из дуг АВ, БС, АБ, ВС. Далее на основе

алгебраического метода строится начальное приближение и методом последовательной верхней релаксации решаются уравнения ЕБ2-метода для определения координат внутренних узлов сетки. Сетка, показанная на рис. 4, построена при управляющей функции

Цж1 , ж2) = 1 + , ^ 0. , (8)

р(х1, ж2)

где р — расстояние от точки (ж1, ж2) до центра острова. С ростом параметра а сетка в окрестности острова сгущается.

3. Итерационный метод

Итерационный процесс, как и в односвязной области [4], состоит из трех этапов: расчета потенциального течения, расчета вихревого течения "под крышкой" и расчета вихревого течения со свободной границей.

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

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

Зависимость среднего значения функции тока от глубины отмели

Ло ^ср Ло ^ср

0.5 1.444 0.0625 1.419

0.25 1.433 0.0313 1.417

0.125 1.423 0.0157 1.415

Рис. 5. Река с мелью на глубине 0.0313 м.

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

) = Р0 + J / + HJv2(u + /cor)) dq1 + / - HJv\u + /cor)) dq2 7(qo>qj)

(9)

где р0 = Н/2, Н0 — полная глубина, заданная в некоторой точке х0 £ Гт, прообразом

которой в вычислительной области является точка q0; 3 — якобиан преобразования (7); 1 2

Vт, V2 — контравариантные компоненты скорости:

/а = H

dq0

1 дф HJdq2

|u|2

2 1 дф

v = —■

HJ dq1'

h-

2

v0|u|,

J Ch

а

1, 2,

va — ковариантные компоненты скорости; через /cor и /ch обозначены безразмерные значения параметра Кориолиса и коэффициента Шези: /cor = k\Jh0/g, /ch = C/^/g.

В работе [4] показано, что при специальной аппроксимации соотношения (9), согласованной с аппроксимацией уравнения для вихря, величина p = H 2/2 не будет зависеть от контура интегрирования Y(q0, qj). Для двусвязной области также потребуем выполнения условия независимости p от контура интегрирования, и значение фостр будем определять из этого условия.

Пусть y1 и y2 — кривые из формулы (9), связывающие точку qj, в которой требуется определить p, с точкой q0, в которой задано p0. Предположим вначале, что область D, ограниченная этими кривыми, не содержит в себе YcT"^ (т.е. является односвязной). Границу области D обозначим через L. Тогда условие независимости p от кривых y 1 и y2 сведется к равенству

У (/1 + HJv2(u + /cor)) dq1 + (/2 - HJv\u + /cor)) dq2 = 0.

Используя формулу Грина, это условие можно переписать так:

д д — (HJv1(w + /cor)) + ^ (HJv2(^ + /cor))

D

dq1 dq2

D

/ + / dq2 dq1

dq1dq2.

С учетом уравнения неразрывности

(10)

|т (HJvl) + А (НЛ2) = 0

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

1

v

L

остр 1 2 тт

контура острова y0 , то величина p не зависит от кривых y1 и y2. На разностном уровне данное утверждение доказывается, как ив [4].

Теперь рассмотрим случай, когда кривые y 1 и y2 образуют замкнутую кривую L, охватывающую прообраз границы острова YoCTP, т.е. область D является двусвязной. Тогда

остр

криволинейный интеграл по кривой L равен криволинейному интегралу по контуру y0 [9]. Следовательно, условие (10) независимости p от кривой L примет вид

У (/ + HJv2(w + /cor)) dq1 + (/2 - HJvV + /Cor)) dq2 = 0. (11)

остр

Yo

С использованием условий непротекания (4), которые в криволинейных координатах имеют вид (см. рис. 2)

' ' ' ' (12)

AD

0, v2

BC

AB

0,

DC

условие (11) можно записать так:

остр

Yo

д_

öq1

|u|

Я-^т ( h - -гV1|u|

f Ch

2

1

dq1 +

^ /Ch V2H

'Ch

dq2 = 0. (13)

Предположим, что глубина дна Н около берегов острова — величина постоянная. Тогда с учетом равенств

V1

|u|2

AB

AB

|u|2

V1

DC

DC

= g11v

V2

AD

= V2

BC

= g22V

g11 (v^2 , |u|2

AD

|u|2

BC

следующих из (12), условие (13) примет вид

/ г2(g11 (v1)2)- /Cch^v^|v11

остр

Yo

g22 v

dq1+

2) 2

+

2^ (g22 И2) - /Chg22v2^22|v2

dq2

0.

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

Сделаем несколько замечаний об особенностях итерационного процесса расчета течения в речном русле с островом по сравнению с алгоритмом, изложенным в [4]. Функция тока ф по-прежнему вычисляется с помощью метода последовательной верхней релаксации, но теперь только во внутренних узлах сетки, покрывающей двусвязную область

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

2

4. Пример расчета

Приведем некоторые результаты расчетов для водоема с неровным дном при тех же значениях параметров, что ив [4], но при наличии острова:

к = 1.19 • 10-4, С = 45 м1/2/с, ивх = 0.5 м/с,

где ивх — скорость втекания жидкости через участок Гт. Вектор потока массы Vт и расход из условий (3) задаются следующим образом:

V т (х)

хеГ1

(к(х)ивх, 0) , ^2 (х)

хеГ2

^(х)ивых,

где постоянная ивых определяется из условия равенства прихода жидкости через Гт и ее расхода через Г2.

Рис. 6. Линии тока установившегося вихревого течения.

Рис. 7. Поле вектора скорости установившегося вихревого течения.

Рис. 8. Изолинии полной глубины.

Рис. 9. Поле вектора скорости установившегося вихревого течения вокруг острова.

Расчеты проводились на криволинейной адаптивной сетке (см. рис. 4). На рис. 6 показаны линии тока установившегося вихревого течения жидкости со свободной границей. На рис. 7 приведено соответствующее поле вектора скорости. На рис. 8 показаны изолинии полной глубины. На рис. 9 приведено поле вектора скорости около острова в увеличенном виде.

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

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

[1] BORTHWICK A.G.L., Barber R.W. River and reservoir flow modelling using the transformed shallow water equations // Intern. J. Numer. Methods Fluids. 1992. Vol. 14. P. 1193-1217.

[2] Kashiyama K., ITO H., Behr M., Tezduyar T. Three-step explicit finite element computation of shallow water flows on a massively parallel computer // Intern. J. Numer. Methods Fluids. 1995. Vol. 21. P. 885-900.

[3] Lin B., Chandler-Wilde S.N. A depth-integrated 2D coastal and estuarine model with conformal boundary-fitted mesh generation // Intern. J. Numer. Methods Fluids. 1996. Vol. 29. P. 819-846.

[4] ХАКИМЗЯНОВ Г.С., ШОКИНА Н.Ю. Численное моделирование установившихся течений жидкости в рамках модели мелкой воды // Вычисл. технологии. 1996. Т. 1, №3. С. 93-105.

[5] КОЧЕРГИН В.П. Теория и методы расчета океанических течений. М.: Наука, 1978. 127. с.

[6] ХАКИМЗЯНОВ Г.С., ШОКИНА Н.Ю. О методе эквираспределения для построения двумерных адаптивных сеток // Вычислительные технологии: Сб. науч. трудов / РАН. Сиб. отд-ние. ИВТ. 1995. Т. 4, №13. С. 271-282.

[7] КУЗНЕЦОВ Б.Г., СИРОЧЕНКО В.П. О постановке задач гидродинамики в многосвязных областях // Вычислительные технологии: Сб. науч. трудов / РАН. Сиб. отд-ние. ИВТ. 1995. Т. 4, №12. С. 209-218.

[8] Суд Э. мл. Численное решение уравнений Навье — Стокса в двусвязных областях для течения несжимаемой жидкости // Ракетная техника и космонавтика. 1974. Т. 12, №5. 1974. С. 76-82.

[9] ФИХТЕНГОЛЬЦ Г.М. Курс дифференциального и интегрального исчисления. Т.3. М.: Наука, 1963.

Поступила в редакцию 24 сентября 2002 г., в переработанном виде — 4 января 2003 г.