Численное моделирование потока жидкости над рельефом дна Текст научной статьи по специальности «Математика»

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2014 Математика и механика № 1(27)

УДК 519.715

В.В. Чуруксаева, М.Д. Михайлов ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПОТОКА ЖИДКОСТИ НАД РЕЛЬЕФОМ ДНА1

Исследуется течение невязкой несжимаемой жидкости на прямолинейном участке канала с негладким дном. Строится несколько математических моделей процесса на основе уравнений мелкой воды. Конструируется математическая модель процесса распространения примеси в потоке. Приводятся результаты их численной реализации.

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

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

Уравнения мелкой воды (МВ) - широко известное приближение, на основе которого проводится численное моделирование разнообразных задач экологии, течений в реках и водохранилищах, прибрежных зонах морей и океанов.

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

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

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

Физическая постановка задачи

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

1 Работа выполнена по заданию Министерства образования и науки РФ № 8.4859.201 и при финансовой поддержке ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 (соглашение № 14.B37.21.0667).

Ня(х/)

к(х^)

2Ь(х)

Рис. 1. Рассматриваемая область

Здесь Ш'( х, /) - уровень свободной поверхности, ZЬ( х) - уровень дна, к( х, /)

- глубина потока

Математическая постановка задачи

Для одномерных течений уравнения МВ имеют вид (а < х < Ь, t > 0)

дк д(иН)

— + ——- = 0;

д1 дх

(1)

д(ки) д(ки2) , дк , dZЬ

——- + —-------- = -як-----як-----

д1 дх дх dx

с соответствующими начальными и граничными условиями:

к(х,0) = к0; и(х,0) = и0;

дк

а0 к(5, () + а! — (5, () = Н;

дх

ди

§0 и(5, () + §1 — (5, 0 = и;

дх

где 5 = а или 5 = Ь.

Неизвестными величинами в системе (1) являются уровень воды к(х, t), измеряемый от отметки дна, профиль которого задается известной функцией ZЬ(х), и скорость потока и(х, t). Здесь я - ускорение свободного падения.

Численный метод решения Дискретизация уравнений проводится на равномерной сетке

й51= {(х/, (п )| х/ = а + /5; ^ = пт, / = 0, N, п = 0,М} ,

где 5, т - шаги по пространству и времени соответственно.

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

Схема 1. Схема с разностями против потока.

Для построения схемы исходные уравнения аппроксимируются с помощью метода конечных объемов.

В центрах ячеек [х/, х/+1] определим дополнительные точки с координатами хс/ = х/+1/2, которые будут являться гранями конечных объемов.

Значение скорости потока u будем определять в полуцелых узлах - xi

г-1/2 :

глубину к - в х/ (см. рис.2); значения в полуцелых (целых) узлах в этом случае вычисляются как среднее арифметическое значений в соседних точках.

нинин

и

хг—1 хг—1/2 хі хі+1/2 хг+1

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

Проинтегрируем уравнения системы (1) в соответствующих «конечных объемах» и по временному интервалу |^п, ^+' ^ и, подставляя полученные приближенные значения интегралов в уравнения (1), запишем разностную схему в виде

тах(и п , 0)к” - тах(-и п ,0)кг” 1 - тах(и п , 0)кг” 1 + тах(-и п ,0)к”

/+1 /+1 /-1 /-1 +---------2----------------2----------------2-----------------2------ = 0;

нп+1 - нп

wn+11 - wn 1

і +— і +—

2___________2_

- тах

( і(

г +— г +—

У 2______________2У

( і(

л л

г +— г +—

У 2______________2У

- тах

( і (

w і + тах

і +—

( 1 (

* 1 г —

-§ (Ж - нп+1нп+1)-) (+1 + нп+і1)+і - іь,)

= 0, і = 0, N -1, п = 1,2,...

"”++1 = и”+1 /(1( + кГ+1')), / = 0, N -1;

1+ 2 1+ 2 ^ ^

к/ = к0; и01 = и0, / = 0, N;

к0п = Н; и0п = и, п = 1М.

В результате исследований получено, что схема имеет погрешность аппрок-симмации 0(5, т), а исследование по фон Нейману дает в качестве условия устой-

т

чивости метода критерий Куранта С = — и < 1 [2].

Схема 2. Двухшаговый метод Лакса — Вендроффа. Запишем исходную систему (1) в векторном виде:

дй дР

где й = (І7,1, Р = У ин

( ин

ни2 + н2

- + — = £,

дґ дх

(0 ^

, ё1Ь

- 8н~Г

их ;

£ =

а

т

т

дт

Начальные и граничные условия для дифференциальной задачи имеют вид

и (х, 0) = и 0;

х 0 и (5, t)+\ ^ (5, t) = ¥,

дх

где 5 = а или 5 = Ь.

Схема основана на разложениях неизвестных функций в ряд Тейлора по времени и пространству до членов второго порядка включительно и представляет собой метод типа «предиктор - корректор»: на первом шаге вычисления производятся по классической схеме Лакса и значения считаются предварительными, а на втором по схеме «чехарда» вычисляются окончательные значения [3].

Разностная аппроксимация системы (1) схемой Лакса - Вендроффа проводилась на равномерной сетке т = {(х{, ^)| х{ = а + /5;^ = пт, / = 0, N, п = 0,М}.

Вычислительные формулы метода имеют вид

иг1 = ■2[й'+1 + й1х]-т -+1---1 +т8п;

2 25

гп+1 - рп+! _ +1

игп+2 = и" - 2т—+-----^ + 2тБп ,

г г 25

/= 1, N -1, п = 0, М -1;

и0 = и0, /= 0,N;

-п+1 —0

и 0 = и , UN = и N-1, п = 0, М -1; (3)

-п+ 2 —0 —п+ 2 —п+2

и 0 = и , UN = и N-1, п = 0, М - 2.

т

Описанный метод устойчив при условии С = — и < 1 и имеет второй порядок

аппроксимации по пространству и времени 0(52, т2). Свойством монотонности

описанная разностная схема не обладает, что характерно для схем второго порядка в целом.

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

В.П. Колганом [4].

Сглаживание здесь достигается добавлением к значениям искомых переменных в узле с индексом / некоторого значения Д :

=п+1 —п+1 —п+1

и =и +Д ;

=п+2 —п+2 —п+2

и =и +Д ;

—n+к

Д

— n+ к — n+к

min(0.25 Ui+і - Ui ,

іб

) х

Ut ) - min(0.25 Ut - Ur -і

— n+ к

n+ к

n+к — n+к — n+к — n+к — n+к — n+ к

- 3Ur + 3Ur-1 - Ur-2 ) sgn(Ur - Ur-1 ),

іб

к = 1, 2.

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

Для реализации численных методов на ПЭВМ на языке C++ (MS Visual Studio 2010) написан комплекс программ с возможностью распараллеливания для вычислительной системы с общей памятью с помощью технологии OpenMP. Приведем решения некоторых известных тестовых задач, полученные с помощью вышеописанных методов.

Задача 1. Задача о распаде разрыва (задача о разрушении плоской плотины).

Рассматривается плоское одномерное течение жидкости в канале длины L с плоским дном Zb(x) = const. В начальный момент в центре области задается разрыв уровня воды, разделяющий два однородных состояния с высотой уровня h = h1 слева от разрыва и h = h2 - справа. В начальный момент времени справа и слева от разрыва жидкость неподвижна, u1 = u2 = 0. Далее представлены расчеты, соответствующие данным, взятым из [і]: L = 2000 м, h1 = 10 м, h2 = 0,1 м . Граничные условия для данного случая имеют вид

Результаты приведены для момента времени t = 50 с. Рис. 3 иллюстрирует распространение волны, образовавшейся при распаде разрыва начальных данных.

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

При сгущении сетки решения, полученные по рассмотренным схемам, становятся практически идентичными и сходятся к аналитическому.

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

Задача 2. Течение над негладким дном.

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

Результаты расчетов

h(0, t) = H; u(0, t) = U;

Рис. 3. Сплошные линии - профили высоты жидкости в момент времени ґ = 50 с (а - про-тивопотоковая схема, б - схема Лакса - Вендроффа, — аналитическое решение,------------ре-

шения на очень мелкой сетке (шаг по пространству 5 и 0,001))

Рис. 4. Решение задачи о течении над неровностями дна - профиль высоты жидкости в момент времени t = 70 с (а - противопотоковая схема, б - схема Лакса - Вендроффа, «*» -решение, полученное в [1])

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

В соответствии с расчетом, представленным в [1], рассмотрим течение на участке канала длиной Ь = 25 м, форма дна которого описывается функцией

2Ь(х) = тах(0, 0,2 -0,05(х-10)2). Математическая модель данного процесса

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

Н(х,0) = НЬх- 1Ь(х); НЬх = 0,33 м; и(х,0) = 0 м/с;

Ни(0, ^ = 0,18 м2/с; — (0,0 = 0; И(Ь, 0 = НЬх - 1Ь{х); и(Ь, 0 = 0. дх

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

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

Задача 3. Рассматривается одномерное нестационарное течение жидкости в канале длины Ь = 2000 м с профилем дна, описываемым функцией

Щх) = Р м Х <Ь/2;

0м

40,1

м, х > Ь /2.

Для решения данной задачи начальные и граничные условия для системы (1) задаются следующим образом:

Н(х,0) = 0,2 м - 1Ь(х); и(х,0) = 1 м/с;

дН

— (0, t) = 0; и(0, t) = 1 м/с.

дх

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

Рис. 5 иллюстрирует профили высоты жидкости, рассчитанные для момента времени t = 20 с.

к, м' 0,2 ■ 0,15 0,1 0,05-

0

т

т

400 800

—I-----1----1---1----1

1200 1600 х, м

Рис. 5. Решение задачи о течении над неровностями дна - профиль высоты жидкости (а -противопотковая схема, б - схема Лакса - Вендроффа, «*» - стационарное решение, взятое из [5])

Важное прикладное значение имеет исследование распространения загрязнений, попадающих в водоем.

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

ди др _ -

дг дх

(к >

где и _ ик , р _

1кс V

( ык

кы2 +—к2 2

икс

dZb

, £ _ -—к——

dx

0

(4)

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

п( х - 500)

и (х, 0) _ и0 и (0, г) _ и0 12

Zb(x) _ тах(0,0, 5,0соб

500

Характеристики течения, полученные с помощью модели (4), (5) и используемые в качестве начальных данных для расчета распространения примеси, представлены на рис. 6.

Результаты расчетов с использованием схем первого и второго порядка представлены на графиках (рис. 7, а и б соответственно).

0

(5)

200 400

600

1—1—I

800 х, м

Рис. 6. Характеристики течения при г = 50 с: профили скорости и(х,г) и высоты жидкости к(х,г), рельеф дна Zb(x,t)

8

Рис. 7. Распределение загрязняющего вещества в потоке в указанные моменты времени (а - противопотковая схема, б - схема Лакса - Вендроффа)

Из рисунка видно, что на участке канала с негладким дном процесс распространения примеси ускоряется, что связано с увеличением скорости потока при перетекании через преграду. Также можно заметить, что схема первого порядка сильнее сглаживает решение около скачка, в то время как решение, полученное по схеме второго порядка, ближе к аналитическому.

Фоновая концентрация загрязняющего вещества полагалась равной нулю, что не характерно для естественных условий, однако этот факт не оказывает существенного влияния на характер описываемого процесса. Концентрация, задающая поступление примеси, равна 0,5 мг/л, что соответствует измерениям концентрации нефтепродуктов, проводившимся на р. Томь.

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

Выводы

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

В работе были представлены некоторые модельные задачи, являющиеся приближением задачи о течении в реке, а также распространения примеси в речном потоке. В работе использовались разностные схемы различного порядка точности, одна из которых была построена методом конечных объемов, другая - взята из [3].

Результаты численных расчетов представлены в виде графиков изменения глубины и скорости потока (рис. 3-7). Сравнение полученных приближенных решений с результатами из [1, 5] показало совпадение для обоих численных методов, что говорит об адекватности построенных моделей описываемому физическому процессу, о корректности расчетов и делает возможным обобщение простейших моделей на более сложный случай, когда система уравнений мелкой воды позволяет более детально описать течение и учесть множество факторов, определяющих его в естественных условиях.

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

ЛИТЕРАТУРА

1. Булатов О.В., Елизарова Т.Г. Регуляризованные уравнения мелкой воды и эффективный метод численного моделирования течений в неглубоких водоемах // Журнал вычислительной математики и математической физики. 2011. Т. 51. № 1. С. 170-184.

2. Рихтмайер Р., Мортон К. Разностные методы решения краевых задач. М.: Мир, 1972. 418 с.

3. Роуч П. Вычислительная гидродинамика. М.: Мир, 1976. 612 с.

4. Колган В.П. Применение операторов сглаживания в разностных схемах высокого порядка точности // Вычислительная математика и математическая физика. 1978. Т. 18. № 5. С. 1340-1345.

5. Петросян А.С. Дополнительные главы гидродинамики тяжелой жидкости. М.: ИКИ РАН, 2010. 127 с.

6. Чуруксаева В.В., Михайлов М.Д. Математическое моделирование процессов самоочищения реки // Материалы юбилейной 50-й Международной научной студенческой конференции «Студент и научно-технический прогресс». Новосибирск: Изд-во НГТУ, 2012.

С. 289.

Статья поступила 26.09.2013 г.

Churuksaeva V.V., Mikhailov M.D. NUMERICAL MODELLING OF THE FLUID FLOW ABOVE THE BOTTOM TOPOGRAPHY. This paper presents an investigation of an inviscid incompressible fluid flow in a straight section of a channel with an irregular bottom as a closure of river stream model. Mathematically, the problem is written as a boundary-value problem for shallow water equations.

Three test computational examples for a steady and unsteady flow above regular and irregular bottom have been carried out to study the model and possibilities of its applications.

The computed solutions are obtained using the finite-difference method with the first order UPWIND scheme and two-step Lax-Wendroff scheme, which is second-order accurate in both space and time. To suppress dispersion characteristics which are the feature of second-order schemes, Kolgan’s surfacing algorithm is used. Numerical solutions obtained by the aforesaid schemes well agree with each other and become equivalent upon mesh clustering.

In addition, a model of the contaminant dispersion in a stream over an irregular bottom is constructed. The computed distribution of the contaminant is in a good agreement with the physical flow pattern.

Keywords: mathematical model, shallow water equations, approximation error, solution stability, solution smoothing.

REFERENCES

1. Bulatov O.V., Elizarova T.G. Regulyarizovannye uravneniya melkoy vody i effektivnyy metod chislennogo modelirovaniya techeniy v neglubokikh vodoemakh // Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki. 2011. V. 51. No. 1. P. 170-184 (in Russian).

2. Rikhtmayer R., Morton K. Raznostnye metody resheniya kraevykh zadach. Moscow: Mir, 1972. 418 p. (in Russian).

3. Rouch P. Vychislitel'naya gidrodinamika. Moscow: Mir, 1976. 612 p. (in Russian).

4. Kolgan V.P. Primenenie operatorov sglazhivaniya v raznostnykh skhemakh vysokogo pory-adka tochnosti // Vychislitel'naya matematika i matematicheskaya fizika. 1978. V. 18. No. 5. P. 1340-1345 (in Russian).

5. Petrosyan A.S. Dopolnitel'nye glavy gidrodinamiki tyazheloy zhidkosti. Moscow: IKI RAN, 2010. 127 p. (in Russian).

6. Churuksaeva V.V., Mikhaylov M.D. Matematicheskoe modelirovanie protsessov samoochish-cheniya reki // Materialy yubileynoy 50-y Mezhdunarodnoy nauchnoy studencheskoy konfer-entsii «Student i nauchno-tekhnicheskiy progress». Novosibirsk: Izd-vo NGTU, 2012. P. 289 (in Russian).

CHURUKSAEVA Vladislava Vasilievna (Tomsk State University, Tomsk, Russian Federation) E-mail: [email protected]

MIKHAILOV Mikhail Dmitrievich (Tomsk State University, Tomsk, Russian Federation)

E-mail: [email protected]