Numerical modelling of the fluid flow above the bottom topography.pdf Для корректного построения математических моделей и численных методов расчета течений в реках важно знать, каким образом описываются основные процессы, происходящие в водоеме. Для этого строятся более простые модели для течения в каналах, учитывающие, однако, основные характеристики течения в природной среде. Уравнения мелкой воды (МВ) - широко известное приближение, на основе которого проводится численное моделирование разнообразных задач экологии, течений в реках и водохранилищах, прибрежных зонах морей и океанов. При выводе данных уравнений предполагается, что среда представляет собой достаточно тонкий слой, глубина которого много меньше его продольного размера, поэтому вертикальной составляющей скорости в слое можно пренебречь и полагать, что продольные скорости постоянны по толщине слоя [1]. Исходя из этого, данное приближение успешно применяется для описания течений, где значительно влияние на поток свободной поверхности и рельефа дна. Цель данной работы - построение и численная реализация математической модели на основе одномерных уравнений мелкой воды, описывающей течение в канале и учитывающей факторы, оказывающие влияние на характер течения в естественном водоеме. Одномерным можно считать течение жидкости в канале с плавно изменяющимся поперечным сечением и малой кривизной оси. При этом принимается, что все параметры течения постоянны в поперечном сечении потока. Полученные в рамках такой простейшей модели решения, естественно, носят приближенный характер, но во многих случаях хорошо совпадают с экспериментальными данными. Физическая постановка задачи Рассматривается течение идеальной жидкости на участке прямолинейного бесконечного в обе стороны канала со свободной поверхностью и дном переменной высоты. Zb(x) Рис. 1. Рассматриваемая область Здесь Hs( x, t) - уровень свободной поверхности, Zb( x) - уровень дна, h( x, t) - глубина потока Математическая постановка задачи Для одномерных течений уравнения МВ имеют вид (a < x < b, t > 0) dh d(uh) Л - + --- = 0; dt dx (1) d(hu) d(hu2) , dh , dZb --- + --- = -gh--ghdt dx dx dx c соответствующими начальными и граничными условиями: h(x,0) = h0; u(x,0) = u0; dh a0 h& t) + a1 - (I, t) = H; dx du 50 u(|, t) + 51 - t) = U; dx где | = a или | = b. Неизвестными величинами в системе (1) являются уровень воды h(x, t), измеряемый от отметки дна, профиль которого задается известной функцией Zb(x), и скорость потока u(x, t). Здесь g - ускорение свободного падения. Численный метод решения Дискретизация уравнений проводится на равномерной сетке rasT ={(xi,tn)|xt = a + is;tn = nx, i = 0"N, n = 0"M} , где s, т - шаги по пространству и времени соответственно. Для численного решения уравнений (1) будем использовать явные по времени разностные схемы первого и второго порядков точности. Первая из них отличается простотой реализации и небольшой трудоемкостью, вторая - более высокой точностью. Схема 1. Схема с разностями против потока. Для построения схемы исходные уравнения аппроксимируются с помощью метода конечных объемов. В центрах ячеек [xi, xi+1] определим дополнительные точки с координатами xct = xt+1/2, которые будут являться гранями конечных объемов. Hs(x,t) h(x,t) Значение скорости потока u будем определять в полуцелых узлах - xi а i-1/2 : глубину h - в xi (см. рис.2); значения в полуцелых (целых) узлах в этом случае вычисляются как среднее арифметическое значений в соседних точках. U h U h U h max(un ,0)h" - max(-un ,0)h"+1 - max(un ,0)hin-1 + max(-un ,0)hf i+1 i+1 i-1 i-1 +-2-2-2-2-= 0; Л Л - max 1 2 V un 3 +un J i +- i +- V 2 2 У ,0 у wn 3 + max i +- 2 1 2 V un 3 +un J i +- i +- V 2 2 У ,0 у 1 г +- 2 x h"+1 - h" x wn+l - wn 1 i +- i +- 2_2_ Л Л - max 1 unj +unj ,0 wn j + max 1 un J +unJ ,0 2 i +- i- i +- 2 i +- i - V V 2 2 У у 2 V V 2 2 У у 1 i - 2 s -2 (Ж1 - hin+1hin+1)-g (( + h+1 )(Zbi+1 - Zbt) = 0, i = 0,N -1, n = 1,2,... s w"! = uj + h"++\ )), i = 0, N -1; '+ 2 '+ 2 V У h0 = h0; u0 = u0, i = 0, N; h0n = H; u0n = U, n = 1M. В результате исследований получено, что схема имеет погрешность аппрок-симмации O(s, т), а исследование по фон Нейману дает в качестве условия устой чивости метода критерий Куранта C = - u < 1 [2]. s Схема 2. Двухшаговый метод Лакса - Вендроффа. Запишем исходную систему (1) в векторном виде: dU dF - + - = S, dt dx где U = fh 1, F = V uh f uh hu2 + gh2 л 1 _ f 0 S = , dZb - 8h~T dx у Начальные и граничные условия для дифференциальной задачи имеют вид U (x,0) = U 0; x 0 U (I, t)+x ^L (I, t) = W, dx где | = a или | = b. Схема основана на разложениях неизвестных функций в ряд Тейлора по времени и пространству до членов второго порядка включительно и представляет собой метод типа «предиктор - корректор»: на первом шаге вычисления производятся по классической схеме Лакса и значения считаются предварительными, а на втором по схеме «чехарда» вычисляются окончательные значения [3]. Разностная аппроксимация системы (1) схемой Лакса - Вендроффа проводилась на равномерной сетке ras т = {(xt, tn)| xt = a + is;tn = пт, i = 0, N, n = 0,M}. Вычислительные формулы метода имеют вид Ur1 = 2[Ui+1 + Uh]-т +TSn; 2 2s Fn+1 - Fn+1 _ ,1 U"+2 = U" - 2т F+1 F-1 + 2tS"+1 , г г 2s i = 1, N -1, n = 0, M -1; U0 = U0, i = 0,N; -n+1 -0 U0 = U , Un = Un-1, n = 0,M -1; (3) -n+ 2 -0 -n+ 2 -n+2 U00 = U , Un = Un-1, n = 0,M - 2. т Описанный метод устойчив при условии C =- u < 1 и имеет второй порядок s аппроксимации по пространству и времени O(s2, т2). Свойством монотонности описанная разностная схема не обладает, что характерно для схем второго порядка в целом. Приведенная разностная схема имеет незначительную схемную вязкость, но высокие дисперсионные свойства и потому вблизи разрывов решение имеет характерные осцилляции. Для уменьшения их амплитуды в схему, как правило (например, в [1], [2]) добавляется искусственная вязкость, которая, однако, накладывает дополнительные ограничения на величину шага по времени. Оптимальным в этом случае будет применение алгоритма сглаживания, позволяющего уменьшить влияние на численное решение свойств схемы. Такой алгоритм был предложен В.П. Колганом [4]. Сглаживание здесь достигается добавлением к значениям искомых переменных в узле с индексом i некоторого значения Д : =n+1 n+1 -n+1 U =U +Д ; =n+2 -n+2 -n+2 U =U +Д ; -n+к Д = min(0.25 - n+ к - n+к Ui+1 - Ui - n+k - n+к - n+ к - n+ k J_ 16 Ui+2 - 3Ui+1 + 3Ui - Ui -1 X sgn(Ui+1n+к - Ur"+к) - min(0.25 - n+к - n+к - n+к - n+к Ui+1 - 3Ui + 3Ui-1 - Ui-2 - n+ к - n+ к 1 16 Ui - Ui-1 - n+к - n+к ) sgn(Ui - Ui-1 ), к = 1,2. Выбор формулы для вычисления Д определяется порядком схемы, что позволяет сохранить порядок точности аппроксимации исходной дифференциальной задачи неизменным. Этот факт, доказанный В. П. Колганом [4], подтверждается расчетами, проведенными в данной работе. Для реализации численных методов на ПЭВМ на языке С++ (MS Visual Studio 2010) написан комплекс программ с возможностью распараллеливания для вычислительной системы с общей памятью с помощью технологии OpenMP. Приведем решения некоторых известных тестовых задач, полученные с помощью вышеописанных методов. Результаты расчетов Задача 1. Задача о распаде разрыва (задача о разрушении плоской плотины). Рассматривается плоское одномерное течение жидкости в канале длины L с плоским дном Zb(x) = const. В начальный момент в центре области задается разрыв уровня воды, разделяющий два однородных состояния с высотой уровня h = h1 слева от разрыва и h = h2 - справа. В начальный момент времени справа и слева от разрыва жидкость неподвижна, u1 = u2 = 0. Далее представлены расчеты, соответствующие данным, взятым из [1]: L = 2000 м, h1 = 10 м, h2 = 0,1 м . Граничные условия для данного случая имеют вид h(0, t) = H; u(0, t) = U; dh _ т ч „ du /T s - (L, t) = 0; - (L, t) = 0. dx dx Результаты приведены для момента времени t = 50 с. Рис. 3 иллюстрирует распространение волны, образовавшейся при распаде разрыва начальных данных. В целом из рисунка видно, что решения, найденные в данной работе, хорошо согласуются между собой. Можно отметить, что профиль высоты жидкости, полученный с помощью схемы второго порядка, более сглаженный из-за применения алгоритма Колгана. При сгущении сетки решения, полученные по рассмотренным схемам, становятся практически идентичными и сходятся к аналитическому. Кроме того, численные решения согласуются также с результатами, полученными в [1] с помощью «явной по времени схемы с аппроксимацией пространственных производных центральными разностями». Задача 2. Течение над негладким дном. Данный пример позволяет продемонстрировать влияние на поток профиля дна, что является одним из основных факторов в течениях, описываемых уравнениями мелкой воды. Т"1 \ h, м 86h, м 6 4 4 Д 2 2 0 0 "Г "Г т т Т "Г \ 0 400 800 1200 1600 0 x, м "Г 400 800 1200 1600 x,м Рис. 3. Сплошные линии - профили высоты жидкости в момент времени t = 50 c (а - про-тивопотоковая схема, б - схема Лакса - Вендроффа, - аналитическое решение,---решения на очень мелкой сетке (шаг по пространству s « 0,001)) h, м 0,4 0,30,20,1h, м 0,4" 0,3 0,2" 0,1 •*-*-*- -* л 1Л

Скачать электронную версию публикации