Содержание к диссертации
Введение
1 Точность методов с выделением разрывов в случае скалярного закона сохранения 14
1.1 Алгоритм расчета изолированной ударной волны 14
1.2 Задача Коши с разрывными начальными данными 15
1.3 Задача Коши с гладкими начальными данными (висячий скачок) 17
1.4 Асимптотика функции D(x,t) в окрестности точки образования ударной волны 22
1.5 Задача Коши с гладкой начальной функцией, имеющей при минимуме производной точку перегиба более высокого порядка 27
1.6 Асимптотика функции D(x,t) в случае начальных данных, имеющих при минимуме производной точку перегиба п-ого порядка 30
2 Точность методов с выделением разрывов в случае системы законов сохранения 35
2.1 Система уравнений изоэнтропической газовой динамики с показателем адиабаты 7 = 3 35
2.2 Задача Коши с разрывными начальными данными 39
2.3 Задача Коши с гладкими начальными данными 42
3 Моделирование волновых процессов в задаче о разрушении плотины с сухим руслом в нижнем бьефе 50
3.1 Лабораторное моделирование задачи о разрушении плотины с сухим руслом в нижнем бьефе 50
3.2 Базисные законы сохранения теории мелкой воды 52
3.3 Прерывные волны, распространяющиеся по сухому руслу, и их устойчивость 54
3.4 Модифицированный закон сохранения импульса, допускающий распространение прерывных волн по сухому руслу 57
3.5 Задача о разрушении плотины 59
4 Численное моделирование процесса распространения прерывных волн по сухому руслу 65
4.1 Разностная схема 65
4.2 Алгоритм распространения воды по сухому руслу 68
4.3 Результаты одномерных тестовых расчетов 69
4.4 Плановые уравнения теории мелкой воды 75
4.5 Разностная схема в двумерном случае 76
4.6 Численные расчеты в двумерном случае 78
Заключение 80
Список литературы
- Асимптотика функции D(x,t) в окрестности точки образования ударной волны
- Задача Коши с разрывными начальными данными
- Модифицированный закон сохранения импульса, допускающий распространение прерывных волн по сухому руслу
- Результаты одномерных тестовых расчетов
Введение к работе
Настоящая работа посвящена анализу реальной точности расчета нестационарных ударных (прерывных) волн в методах с выделением разрывов, а также разработке численного алгоритма, позволяющего в рамках модифицированных уравнений первого приближения теории мелкой воды моделировать процесс распространения прерывных волн по сухому руслу.
Актуальность. В настоящее время широкое распространение получили разностные схемы повышенной точности для сквозного расчета разрывных решений гиперболических систем законов сохранения [18, 33, 80, 82] (в частности законов сохранения газовой динамики [54] и гидравлики [52, 56]). Однако в большинстве работ, посвященных построению таких схем (см., например, [19, 24, 31, 53, 55, 58, 62, 67, 73, 74, 76, 81, 90]), под точностью схемы понимается порядок ее тейлоровского разложения на гладких решениях, что, как показано в [41, 43] не гарантирует аналогичного повышения порядка слабой аппроксимации на разрывных решениях и, следовательно, не обеспечивает повышенной точности при передаче условий Гюгонио через "размазанные фронты" нестационарных ударных волн [45]. Несмотря на это долгое время было распространено ошибочное мнение о том, что указанные схемы сохраняют повышенный порядок сходимости во всех гладких частях рассчитываемых обобщенных решений. Однако, в работах [20, 44, 65, 71,47] было показано что все эти схемы имеют не более чем первый порядок сходимости в области влияния нестационарной ударной волны (т.е. ударной волны, распространяющейся с переменной скоростью) и тем самым по существу схемами повышенной точности не являются.
Такое снижение реальной точности разностных схем сквозного счета (shock-capturing schemes) привело к тому, что в последние годы заметно повысился интерес к численным методам с выделением разрывов (shock-fitting methods), в основе которых лежит классический метод характеристик [26, 54], часто используемый одновременно с некоторой разностной схемой повышенной точности, применяемой только в областях гладкости рассчитываемого решения [68, 77]. Такие методики, позволяющие отслеживать взаимодействие движущихся разрывов, а также их возникновение и исчезновение с течением времени, иногда называют выделением плавающих разрывов [6, 33, 84]. При этом широко распространен подход, особенно в стационарных зада-
чах обтекания решаемых методом установления, при котором выделяются только основные разрывы, существование которых заранее известно [34, 36, 59]; остальные разрывы рассчитываются по схемам сквозного счета. Эффективным способом выделения разрывов является применение подвижных сеток [23] в сочетании с методом Годунова [22], основанном паточном решении задачи Римана о распаде произвольного разрыва, а также самоподстраивающихся сеток в схемах сквозного счета [30, 75], использующих метод Роу [86] приближенного решения задачи Римана. Несмотря на достаточно широкое распространение численных методов с выделением разрывов до последнего времени отсутствовал детальный анализ их реальной точности при расчете нестационарных ударных волн и прежде всего висячих скачков возникающих внутри расчетной области в результате градиентной катастрофы.
Уравнения первого приближения теории мелкой воды [18, 35, 40], которые аналогичны уравнениям изоэнтропической газовой динамики [54] с показателем адиабаты 7 = 2, представляют собой простейший пример сильно нелинейной строго гиперболической системы законов сохранения [52, 33, 80], имеющей четко выраженный физический смыл [56]. Эти уравнения широко применяются при численном моделировании процесса распространения прерывных волн [4, 16, 17, 42] (гидравлических боров [67, 73, 91]), возникающих при полном или частичном разрушении плотины гидросооружения, или при выходе крупных морских волн типа цунами [61] на мелководье. Известно, что классическая система базисных законов сохранения теории мелкой воды, состоящая из законов сохранения массы и полного импульса [18,33, 52], правильно передавая параметры прерывных волн, распространяющихся по жидкости конечной глубины [56], не допускает распространения прерывных волн по сухому руслу. Точные решения, описывающие в рамках этой системы процесс течения воды по сухому руслу, являются непрерывными волнами понижения (простейший пример такой волны, возникающей при разрушении плотины над горизонтальным дном, впервые был построен в работе [85]). Учет донного трения, которое представляет собой распределенный источниковый член, не дающий вклада в условия Гюгонио на фронте прерывной волны, никак не может изменить эту ситуацию [57, 70, 93]. В то же время лабораторные эксперименты [8, 12, 13, 69, 83, 89] и результаты натурных наблюдений [61] показали, что эти непрерывные решения существенно завышают
скорость распространения передней кромки волны и заметно искажают профиль ее поверхности. В опытах фронт волны, распространяющейся по сухому руслу, является существенно более крутым и на нем происходят обрушения, характерные для прерывных волн.
Необходимо также отметить, что численное моделирование волновых течений по мелководью и сухому руслу в рамках классических законов сохранения массы и полного импульса теории мелкой воды сопряжено с трудностью обеспечения в процессе расчетов неотрицательности глубины жидкости Д. Связано это с тем, что после нахождения расхода q из уравнения для полного импульса скорость потока определяется по формуле v = qjh, которая имеет особенность при h —> 0. В результате практически ни один из известных в настоящее время численных алгоритмов не обеспечивает положительности глубины жидкости h без какого-либо специального сглаживания, искусственной коррекции или вынужденного измельчения шагов по времени или пространству. Так, например, с целью обеспечения неотрицательности глубины применяют специальные итерационные алгоритмы [38]; используют бессеточный метод, основанный на радиальных базисных функциях, с помощью которого в [94] решается одномерная начально-краевая задача с одной фиксированной и одной свободной границей; применяют разбивку области течения на произвольные четырехугольные элементы [92], что позволяет рассматривать каналы не только с сухим дном, но и с нерегулярными границами; используют фиксированную систему конечных элементов, различая сухие, частично смоченные и полностью наполненные жидкостью элементы [64], вводя в уравнение неразрывности функцию, характеризующую наполнение конечного элемента жидкостью. Применяют также нестандартные методы, в которых движение ограниченного объема жидкости описывается не моделью сплошной среды, а механической моделью с конечным числом степеней свободы [1, 79]. При этом при больших значениях глубины практически все известные ранее методы дают вполне удовлетворительные результаты.
Математическое моделирование процесса распространения волн по сухому руслу представляет собой достаточно сложную задачу не только в рамках осредненных по вертикали уравнений мелкой воды, но и в рамках полных уравнений гидродинамики. Как показано в [66], ее решения, получаемые при помощи численных методов,
основанных на уравнениях Эйлера и Навье-Стокса, существенно завышают скорость распространения передней кромки волны. В численных расчетах задачи о разрушении плотины, приведенных в [66], эта скорость асимптотически выходит на скорость распространения соответствующей волны понижения, получаемой из классических уравнений мелкой воды. При этом указанные численные решения могут заметно искажать профиль головной части волны, распространяющейся по сухому руслу. Кроме того на основе полных уравнений гидродинамики, позволяющих описывать волновые течения только в регулярных областях, в принципе нельзя моделировать процесс обрушения волн, приводящий к нарушению односвязности расчетной области [5, 25, 32]. В уравнениях мелкой воды такой процесс обрушения приближенно описывается разрывом параметров течения, т.е. прерывной волной, на фронте которой происходит потеря полной энергии набегающего потока [52] (в рамках теории мелкой воды потеря полной энергии на прерывных волнах означает, что в реальном течении некоторая ее часть переходит в энергию движения вертикальных вихрей). В результате уравнения первого приближения теории мелкой воды представляют собой одну из наиболее распространенных и широко применяемых на практике моделей для описания волновых процессов [4, 60, 67, 73, 88] (в том числе и катастрофических [16, 17]) в системах открытых русел и в прибрежных зонах морей и океанов. Не смотря на это до последнего времени отсутствовали эффективные численные алгоритмы, которые бы позволяли в рамках этой модели описывать процесс распространения прерывных волн по мелководью и сухому руслу.
Целью диссертационной работы является анализ реальной точности методов с выделением разрывов и разработка численного алгоритма, позволяющего в рамках теории мелкой воды моделировать процесс распространения прерывных волн по сухому руслу.
Научная новизна:
Впервые теоретически и численно проведен анализ максимально возможной точности, которую могут допускать численные методы с выделением разрывов при расчете по ним нестационарных ударных волн и прежде всего висячих скачков, возникающих внутри расчетной области в результате градиентной катастрофы.
Впервые предложена и верифицирована путем согласования с экспериментом
модифицированная система базисных законов сохранения теории мелкой воды, допускающая распространение по сухому руслу прерывных волн конечной амплитуды. З.На основе метода расщепления по физическим процессам и метода фиктивных ячеек разработан эффективный численный алгоритм сквозного счета, позволяющий в рамках модифицированных уравнений теории мелкой воды моделировать процесс распространения прерывных волн по мелководью, их выход на берег и перетекание через различные береговые препятствия.
Краткое содержание работы. Работа состоит из введения, четырех глав и заключения.
В первой и второй главах изучается максимальная точность, которую могут обеспечить численные методы с выделения разрывов при расчете нестационарных ударных волн и, прежде всего, висячих скачков, возникающих внутри расчетной области при гладких начальных данных. В первой главе точность этих методов изучается в случае скалярного закона сохранения, а во второй главе — в случае системы законов сохранения. В пункте 1.1 описан алгоритм расчета пространственно одномерной изолированной ударной волны с помощью классических методов выделения разрывов и процедура проверки реальной точности этих методов путем проведения трех численных расчетов на сетках с временными шагами т, т/2 и т/4. Проверка точности методов выделения разрывов будет проводиться на этапе численного интегрирования методами Рунге-Кутта 1-го, 2-го и 4-го порядков обыкновенного дифференциального уравнения dx/dt = D(x,t), где D(x,t) - скорость распространения прерывной волны, в предположении что на всех предшествующих этапах вычисления проведены с высокой точностью. В пункте 1.2 методом выделения разрыва решается задача Коши для нелинейного уравнения переноса
щ + (и2/2)х = 0, u(0,x) = v(x)
с кусочно-линейными начальными данными. Вычисляются реальные порядки сходимости и ошибки разностного решения, из которых видно, что в данном случае реальная точность методов Рунге-Кутта совпадает с их формальной точностью. В пункте 1.3 рассматривается та же задача Коши для нелинейного уравнения перено-
са, но уже с гладкой начальной функцией v(x,t) со следующими свойствами:
v{0) = c>0, v'(x)>v'(0) = -b<0 Ух ^0, i/'(0) = 0, i/"(0) = d>0,
из которых следует, что производная функции и (г) достигает своего абсолютного отрицательного минимума в точке х = 0, являющейся для нее точкой перегиба первого порядка. Для функции с такими свойствами описана процедура точного вычисления координат точки, в которой внутри расчетной области происходит градиентная катастрофа, и методика построения начальной функции такого вида. Вычисляются реальные порядки сходимости и ошибки разностного решения для данной задачи, из которых видно, что в данном случае методы Рунге-Кутта обеспечивают скорость сходимости к линии фронта ниже их формальной точности. Для определения причины снижения точности в пункте 1.4 строится асимптотика функции D(x, t) в окрестности точки образования ударной волны в случае начальной функции, имеющей в нуле точку перегиба первого порядка, и показано, что функция D(x, t) имеет в окрестности это точки корневую особенность, что и объясняет снижение точности методов Рунге-Кутта при вычисление висячих скачков, так как эти методы сохраняют повышенный порядок сходимости только при условии достаточной гладкости интегрируемой функции. В пункте 1.5 рассматривается задача Коши для нелинейного уравнения переноса с гладкой начальной функцией, имеющей в нуле точку перегиба п-го порядка:
v(0) = c>0, v'(x)>v'(0) = -b<0 Vx^O,
v{k)(0) = 0, к = 2^п, и(2п+1)(0) = d > 0.
Построены начальные функции, имеющие в нуле точку перегиба 2-го и 3-го порядков. Для данных начальных функций вычисляются реальные порядки сходимости и ошибки разностного решения при г = 0.06, из которых видно заметное снижение, по сравнению со случаем п = 1, реальной скорости сходимости метода Рунге-Кутта четвертого порядка: при п = 2,3 реальный порядок точности Д* этой схемы не только становится заметно ниже двух, но оказывается даже ниже реальной скорости сходимости R2 метода Рунге-Кутта второго порядка. Для объяснения кажущегося несоответствия проведены дополнительные расчеты этих задач на более мелкой сетке с шагом г = 0.006. В пункте 1.6 показано, что дальнейшее снижение при
n = 2,3 (по сравнению с n = 1) реальной точности методов Рунге-Кутта связано с увеличением степени корневой особенности у функции D(x, t) в окрестности точки градиентной катастрофы, из которой образуется висячий скачок, что обнаруживается путем построения ее асимптотики в окрестности этой точки.
В пункте 2.1 рассматривается система уравнений изоэнтропического газа [54] с показателем адиабаты 7 = 3
Pt + qx = 0, qt + (qu + p3/3)x = 0,
которая в области гладкости своего решения распадается на два линейно-независимых уравнения переноса для инвариантов системы г — и — р и s = и + р :
П + (г2/2)х = 0, st + (s2/2)x = О,
что удобно для проведения численных тестовых расчетов. Проводится анализ этой системы, описывается алгоритм решения задачи Коши для этой системы методом выделения разрывов. В пункте 2.2 для системы уравнений изоэнтропического газа с показателем адиабаты 7 = 3 решается задача Коши с разрывными начальными данными. Приведены порядки сходимости и ошибки вычисления линии фронта ударной волны и "проходящего" инварианта на отрезке, лежащем в области влияния ударной волны, из которых видно, что, также как и в случае скалярного закона сохранения, при разрывных начальных данных метод выделения разрыва обеспечивает реальный порядок сходимости, совпадающий с формальной точностью соответствующей ему схемы Рунге-Кутта, как при расчете линии фронта ударной волны, так и при вычислении инварианта в области ее влияния. В пункте 2.3 для системы уравнений изоэнтропического газа с показателем адиабаты 7 = 3 решается задача Коши с гладкими начальными данными
г(0, х) - -1, s(0, х) = vn(x).
По аналогии со скалярным случаем рассматривается три варианта, а именно когда начальная функция vn(x) имеет в нуле точку перегиба первого, второго и третьего порядка. Из проведенных расчетов видно, что, также как и в случае скалярного закона сохранения, при решении задачи Коши с гладкими начальными данными методы Рунге-Кутта второго и четвертого порядков обеспечивают скорость сходимости
как к линии фронта ударной волны, так и к инварианту г в области ее влияния, заметно ниже их формальной точности. Проведен детальный сравнительный анализ всех численных расчетов, который показывает, что при переходе от скалярного уравнения к системе при всех п = 1,2,3 реальная точность методов Рунге-Кутта второго и четвертого порядков становится несколько ниже, а метода первого порядка - несколько выше.
В третьей главе предложен метод, позволяющий в рамках первого приближения теории мелкой воды моделировать процесс распространения прерывных волн по сухому руслу. В пункте 3.1 демонстрируется несоответствие результатов лабораторных экспериментов и натурных наблюдений, в которых прерывные волны распространяются по сухому руслу, с результатами численного моделирование аналогичных процессов в рамках классических уравнений теории мелкой воды. В пункте 3.2 рассматривается система базисных законов сохранения теории мелкой воды (система уравнений Сен-Венана) в случае прямоугольного русла постоянной ширины и переменной глубины, состоящая из закона сохранения массы и полного импульса. Показано, что эта система не допускает распространения прерывных волн конечной амплитуды по сухому руслу. В пункте 3.3 обнаружено, что скорость прерывных волн, распространяющихся по сухому руслу, совпадает со скоростью жидкости за их фронтом (D = щ), что на таких волнах происходит потеря полного импульса (5q < 0) и полной энергии (5е < 0) набегающего потока, что такие волны удовлетворяют энтропийному и эволюционному критерию устойчивости. В пункте 3.4 получена система законов сохранения массы и модифицированного импульса, представляющего собой линейную комбинацию законов сохранения полного и локального импульсов, в который входит эвристический параметр, подбираемый путем согласования с экспериментом. Показано, что получаемое из нее условие Гюгонио с одной стороны допускает распространение прерывных волн по сухому руслу, а с другой стороны с ростом глубины ho перед фронтом прерывной волны достаточно быстро асимптотически переходит в классическое соотношение. В пункте 3.5 проводится сравнение решения задачи о разрушении плотины по классической и модифицированной системе уравнений мелкой воды как в случае наличия жидкости, так и в случае ее отсутствия в нижнем бьефе. Показано, что заметные различия между классическим
и модифицированным решением имеют место только при do < 0.14 и особенно при do —> 0, где d0 - отношение начальных глубин жидкости в нижнем и верхнем бьефе плотины.
В четвертой главе предложен эффективный численный алгоритм сквозного счета, позволяющий на основе модифицированных базисных законов сохранения первого приближения теории мелкой воды численно моделировать процесс распространения прерывных волн по сухому руслу. В пункте 4.1 для модифицированной системы уравнений мелкой воды методом расщепления по физическим процессам на разнесенной по пространству сетке строится условно устойчивая разностная схема. Данная схема с учетом искусственных вязкостей, коэффициенты которых выбираются из тестовых расчетов, имеет первый порядок аппроксимации как по времени, так и по пространству и обладает заметно более высокой разрешимостью по сравнению с монотонными схемами формально повышенной точности (типа TVD или ENO) при расчете сложных волновых течений. В пункте 4.2 описан алгоритм, позволяющий моделировать движение свободной границы жидкости внутри расчетной области на основе метода "фиктивных ячеек". В пункте 4.3 приводятся результаты численного моделирования процесса формирования и распространения по сухому руслу прерывной волны, возникающей в результате разрушения плотины и резкого локального подъема дна (волны типа цунами). Рассматриваются случаи горизонтального и наклонного дна, а также дна с локальным препятствием в нижнем бьефе. Значения коэффициентов искусственных вязкостей выбираются как минимальные среди тех значений, которые обеспечивают отсутствие осцилляции за фронтами прерывных волн. В пункте 4.4 по аналогии с одномерным случаем получены пространственно двумерные уравнения модифицированного импульса, которые будут использоваться для численного моделирования задач о распространении прерывных волн по сухому руслу. В пункте 4.5 строится разностная схема для модифицированной системы уравнений мелкой воды. В пункте 4.6 приведены результаты численного моделирования задачи о распространении прерывной волны, возникающей в результате частичного разрушения плотины в русле с горизонтальным и трапециевидным дном, горизонтальным дном, имеющим в нижнем бьефе препятствие в виде полуэллипса, задачи о набегании на наклонный берег прерывной волны, возникающей в результате частич-
ного разрушения плотины и задачи о набегании на наклонный берег, на котором расположено локальное препятсвие в виде полуэлипсоида, прерывной волны типа цунами, возникающей в результате резкого подъема локального участка дна. В заключении приведены основные результаты работы.
Апробация работы. Основные положения и результаты диссертационной работы докладывались и обсуждались на следующих конференциях и семинарах: XXXIX, XLII Международных научных конференциях "Студент и научно-технический прогресс" (Новосибирск, 2001; 2004); Четвертой Сибирской школе-семинаре "Математические проблемы механики сплошных сред" (Новосибирск, 2001); Международной конференции "Вычислительные и информационные технологии в науке, технике и образовании" (Казахстан, Усть-Каменогорск, 2003); Всероссийской конференции "Новые математические модели в механике сплошных сред: построение и изучение" (Новосибирск, 2004); Международной конференции по вычислительной математике МКВМ-2004 (Новосибирск, 2004); на семинарах института математического моделирования РАН, института гидродинамики им. М.А. Лаврентьева СО РАН; института вычислительных технологий СО РАН, института вычислительной математики и математической геофизики СО РАН, кафедры дифференциальных уравнений НГУ.
Асимптотика функции D(x,t) в окрестности точки образования ударной волны
Для введения несимметрии, не нарушающей свойств (1.13), (1.21), добавим к функции (1.19) слагаемое г)(х) = 5х4еахх{—х) четвертого порядка малости по х, в котором 5 = 0.1 и (1.22) Х(х) = 0, х 0 1, х - функция Хевисайда. В результате получим следующую несимметричную функцию / ч / \ , / \ і -х(1 -ах-5х3)еах, х 0 v(x) = w(x) + ф) = { _х{1 + ах) е-ах х 0 (1.23) аналитическую при х Ф 0, имеющую непрерывную третью производную при х = 0, бесконечно малую при х —» со и удовлетворяющую условиям (1.13), (1.21).
Положение висячего скачка, возникающего при решении задачи Коши (1.6), (1.23), определялось в результате численного интегрирования уравнения (1.1), правая часть которого D(x,t) находилась из условия Гюгонио (1.9), в котором и0{х, t) = v{y+(x, t)), щ(х, t) = v{y.{x, t)) (1.24) - значения искомой функции, переносимые вдоль характеристик (1.14). При этом значения у+ и у_ на каждом временном шаге находились в результате числен ного решения с высокой точностью алгебраических уравнений (1.14). Поскольку в данном случае положение линии фронта заранее неизвестно, то порядок сходимости R и ошибка разностного решения 8ХТ определялись по формулам (1.5). Резуль таты численных расчетов для сетки с шагом т = 0.06 приведены в таблице 2, из которой видно, что в данном случае методы Рунге-Кутта 2-го и 4-го порядков обеспечивают скорость сходимости к линии фронта ниже их формальной точности. Как будет показано в следующем пункте причина этого заключается в наличии корневой особенности у функции D(x, t) в окрестности точки градиентной катастрофы.
Для того, чтобы объяснить снижение скорости сходимости методов Рунге-Кутта второго (1.3) и четвертого (1.4) порядков при решении задачи (1.6) с гладкими несимметричными начальными данными (1.23), удовлетворяющими условиям (1.13), построим асимптотику функции D(x,t) в окрестности точки образования ударной волны. Для этого, следуя [29], удобно сделать следующие преобразования переменных: а) совместное растяжение независимых переменных x,t я искомой функции и(х, t) х = рх, t = bt, й = и, P=\hr, (1-25) о V оо в результате которого для новой искомой функции й(х, Г) = й((3х, U) = й(х, t) = -и(х, t) условия (1.13) примут вид й$(0,0) = -их{0,0) = -1, uS2s{0,0) = Гд2«м (0,0) = б, откуда, с учетом (1.16), следует, что ударная волна в новых переменных образуется в момент времени t = 1; б) переход в подвижную систему координат й(х,ї) = й(х,ї)-й(0,0), (1.26) в которой начальная скорость ударной волны с = й(0,0) = 0 и образуется эта. ударная волна, с учетом (1.17), в точке х = 0; в) сдвиг по оси t t = t = t-l, (1.27) переводящий точку градиентной катастрофы (0,1) в начало координат, а линию, на которой задаются начальные данные, в прямую t = —1. В результате в новых переменных х, Ї, и, за которыми сохраним исходные обозначения x,t,u, задача Коши (1.6), (1.13) примет вид щ + (и2/2)х = 0, u(x,-l) = v(x), (1.28) где i»(0) = 0, v (x) i/(0) = -1 Vx O, и"(0) = 0, г/"(0) = 6. (1.29)
С учетом свойств конкретной несимметричной функции (1.23), будем дополнительно предполагать, что функция v(x) аналитическая при всех х/0 ив точке х — О имеет непрерывную третью производную.
При решении задачи (1.28), (1.29) ударная волна зарождается в начале координат х = t = 0 ив момент образования имеет нулевую скорость распространения, т.е. D(Q) = 0. Из (1.28), (1.29) следует, что при t 0 функция D(x,t) получается из условия Гюгонио (1.9) , в котором значения (1.24) переносятся в точку (x,t) вдоль характеристик x = y+(x,t) + (l + t)v(y+(x,t)), x = y-(x,t) + (l + t)v(y4x,t)), (1.30) выходящих из точек у+ 0 и j/_ 0, лежащих на линии начальных данных t = -1 (см. рис. 3). Подставляя (1.30) в (1.9), получим г { .+\ 2x-y-{x,t)-y+(x,t) , . D{x t) = —щтт) (U1) где у+ и у- при малых t однозначно определяются из уравнений (1.30). Построим асимптотику функции (1.31) в окрестности начала координат. Из условий (1.29) следует, что разложение в нуле начальной функции v(y) имеет вид v{y) = -y + g , (1-32) где у Є (0, у). Поскольку точка у в общем случае определена неоднозначно, то через у (у) обозначим такую из них, для которой величина \у \ принимает минимальное значение
Задача Коши с разрывными начальными данными
В данной главе предложен метод, позволяющий в рамках первого приближения теории мелкой воды моделировать процесс распространения прерывных волн по сухому руслу. В основе этого метода лежит модифицированный закон сохранения полного импульса, в котором учитываются сосредоточенные потери импульса, связанные с образованием локальных турбулентно-вихревых структур в поверхностном слое жидкости на фронте прерывной волны. Эвристический параметр, входящий в этот модифицированный закон сохранения, подбирается путем согласования с результатами лабораторных экспериментов. Показано, что в рамках уравнений мелкой воды на фронтах прерывных волн, распространяющихся по сухому руслу, из закона сохранения массы следуют согласованные потери полного импульса и полной энергии набегающего потока. Обсуждается физическая природа этих потерь. Исследована устойчивость прерывных волн, распространяющихся по сухому руслу. В качестве примера построено решение задачи о разрушении плотины с сухим руслом в нижнем бьефе и проведено сравнение этого решения с результатами лабораторных экспериментов.
Уравнения первого приближения теории мелкой воды [18, 35, 40] широко применяются при численном моделировании процесса распространения прерывных волн [4, 16, 17, 42] (гидравлических боров [67, 73, 91]), возникающих при полном или частичном разрушении плотины гидросооружения, или при выходе крупных морских волн типа цунами [61] на мелководье. Известно, что классическая система базисных законов сохранения теории мелкой воды, состоящая из законов сохранения массы и полного импульса [18, 33, 52], правильно передавая параметры прерывных волн, распространяющихся по жидкости конечной глубины [56], не допускает распространения прерывных волн по сухому руслу. Точные решения, описывающие в рамках этой системы процесс течения воды по сухому руслу, являются непрерывными волнами Ь) Рис. 9: Фотография волны в нижнем бьефе при лабораторном моделировании задачи о разрушении плотины, (а): 1 - положение плоского щита, задающего начальный перепад уровней, 2 - головная часть волны, (б): фотография головной части волны, распространяющейся по сухому руслу. АВ - тонкий слой, ВС - гидравлический прыжок
понижения (простейший пример такой волны, возникающей при разрушении плотины над горизонтальным дном, впервые был построен в работе [85). Учет донного трения, которое представляет собой распределенный источниковый член, не дающий вклада в условия Гюгонио на фронте прерывной волны, никак не может изменить эту ситуацию [57, 70, 93). В то же время лабораторные эксперименты (8, 12,13, 69, 83, 89] и результаты натурных наблюдений [61] показывают, что эти непрерывные решения существенно завышают скорость распространения передней кромки волны и заметно искажают профиль ее поверхности. В опытах фронт волны, распространяющейся по сухому руслу, является существенно более крутым и на нем происходят обрушения, характерные для прерывных волн.
На рис. 9 приведена фотография волны, получающейся при лабораторном моделировании задачи о разрушении плотины с сухим руслом в нижнем бьефе. Это волновое течение возникло в результате резкого удаления плоского щита, ограничивающего неподвижную жидкость в верхнем бьефе (подробное описание результатов этого эксперимента дано в [8, 13]). Опыты показали, что профиль головной части волны в процессе ее движения является пульсирующим и квазистационарным. Гидравлический прыжок ВС, набегая на распространяющийся впереди него тонкий слой воды АВ, постепенно поглощает некоторую его часть. При этом крутизна переднего фронта прыжка ВС уменьшается, в результате чего прыжок почти полностью вырождается и профиль свободной поверхности в головной части волны становится достаточно пологим. После этого процесс начинает развиваться в обратном направлении: профиль свободной поверхности в головной части волны на некотором расстоянии от ее кромки постепенно становится все более крутым, что приводит к формированию нового гидравлического прыжка вида ВС. В то же время полоса воды АВ перед его фронтом удлиняется и становится более тонкой. В результате профиль передней части волны вновь принимает вид, показанный на рис. 9, б. Далее этот процесс циклически повторяется.
Далее описанный выше процесс течения по сухому руслу головной части волны моделируется в рамках первого приближения теории мелкой воды как прерывная волна, распространяющаяся по сухому руслу.
Модифицированный закон сохранения импульса, допускающий распространение прерывных волн по сухому руслу
При рассмотрении уравнений мелкой воды с точки зрения общей теории гиперболических систем законов сохранения с выпуклым расширением [72] закон сохранения полной энергии представляет собой замыкающий выпуклый закон сохранения [21], в силу чего потеря полной энергии на разрывах (5е 0) представляет собой энтропийный [80] (энергетический [18], [52]) критерий устойчивости прерывных волн. В рамках теории мелкой воды потеря полной энергии на прерывных волнах означает, что в реальном течении некоторая ее часть переходит в энергию мелкомасштабных вихрей, которые в классическом случае сносятся вверх по потоку от фронта волны, где постепенно затухают под действием вязкости. В течениях, которые описываются прерывными волнами, распространяющимися по сухому руслу, это вихревое движение, в силу условия (3.12), локализуется на фронте волны и переходящая в него энергия 8е с учетом равенства бе = Ui6q, которое следует из (3.15) и (3.19), целиком расходуется на торможение набегающего потока, связанное с потерей его полного импульса.
Прерывные волны (3.11), (3.12) удовлетворяют также эволюционному [33] (характеристическому [52], [80]) условию устойчивости \] D \s = \r, D \\, (3.20) в котором А = Aj? = щ - совпадающие скорости г- и s-характеристик перед фронтом волны в зоне сухого русла ho = щ = О, А = щ — С\ Х\ = щ + d -скорости г— и 5-характеристик за фронтом прерывной волны. Вдоль г-характеристик сохраняется r-инвариант г — и — 2с, а вдоль 5-характеристик - 5-инвариант 5 = и + 2с, где с = л/gh - скорость распространения малых возмущений (скорость звука) в теории мелкой воды [33], [52]. Неравенства (3.20) означают, что на фронт прерывной волны (3.11), (3.12) приходят три характеристики А, Л и X], а уходит с него только одна А -характеристика.
3.4 Модифицированный закон сохранения импульса, допускающий распространение прерывных волн по сухому руслу
Преобразуем закон сохранения полного импульса (3.2) таким образом, чтобы получаемое из него условие Гюгонио с одной стороны допускало распространение прерывных волн (3.11), (3.12) по сухому руслу, а с другой стороны с ростом глубины h0 перед фронтом прерывной волны достаточно быстро асимптотически переходило в классическое соотношение (3.6). Для этого рассмотрим закон сохранения локального импульса (3.17) с учетом (3.4), для которого условие Гюгонио имеет вид
Bi B = hohi (h) Ah0hx(hY то условие Гюгонио (3.6) и (3.21) для законов сохранения полного (3.2) и локального (3.17) импульсов с учетом условия Гюгонио (3.5) для закона сохранения массы (3.1) несовместны на любых разрывах конечной амплитуды [h] ф 0. Поскольку соотношение (3.23), в отличие от соотношения (3.10), при ho - 0 дает конечную зависимость щ = y/2ghi, то система законов сохранения массы (3.1) и локального импульса (3.17), в отличие от классической системы (3.1), (3.2), допускает распространение по сухому руслу прерывных волн конечной амплитуды. С учетом этого рассмотрим следующую линейную комбинацию законов сохранения полного (3.2) и локального (3.17) импульсов qt+(qu + -j +4\ut+\ + ghj)=R-ghbx + 1(--gbxj, (3.24) 7 = SH = const, где Я - некоторая характерная для данного течения постоянная глубина, 6 - малый безразмерный параметр, выбираемый путем согласования с результатами лабораторных экспериментов. Вводя обозначения Q = q + fu = hu, /г = /і + 7) h = h + -, /1 = /1 + 27, (3.25) R=(l + l)R= g Q\u\ /i4/3 (3.26) перепишем уравнение (3.24) в виде закона сохранения модифицированного импульса Q: LL (3.27) Qt + I hu2 + д— I = R - ghh Рассмотрим систему законов сохранения массы (3.1) и модифицированного импульса (3.27) и выпишем определитель матрицы Якоби преобразования v = v(u), где u = (h,q), v - (h,Q), переводящий систему (3.1), (3.2) в систему (3.1), (3.27): - 4 1 0 (3.28) Vn = и л Ч 1 + 1 Из формулы (3.28) следует, что на гладких решениях системы (3.1), (3.2) и (3.1), (3.27) эквивалентны при всех h ф 0, т.е. во всей области гиперболичности h 0 модели мелкой воды.
Результаты одномерных тестовых расчетов
Если в результате расчета в некоторой ячейке 1а С 1п оказалось, что /г+1 є, то в ней полагается /І+1 = є. Если при этом во всех смежных с ней ячейках (т.е. имеющих с ней общий целый узел ХІ или ХІ+І, где і = а — 1/2) глубина hn+1 = є или уровень zn+1 z%+1, то ячейка 1а на (n + 1)-м временном слое становится "фиктивной", т.е. Ia . In+1. Наоборот, если в некоторой ячейке, смежной с ячейкой 1а 1п, оказалось, что hn+1 є и zn+1 z%+l, то ячейка Ia на (n + 1)-м временном слое переходит в разряд "рабочих", т.е. 1а С In+1. Таким образом происходит изменение со временем внутренней расчетной области, что позволяет при помощи сквозного счета моделировать набегание волн на берег, стекание воды с берега, а также перетекание воды через различные береговые препятствия.
Замечание. Предложенный алгоритм распространения воды по сухому руслу является корректным только при непрерывном изменении отметки дна Ь{х) в русле. Как показано в работах [14, 15, 50, 51] для правильного моделирования волновых течений над скачком отметки дна необходима дополнительная модификация уравнений мелкой воды, связанная с заданием специальных условий на этом скачке. 4.3 Результаты одномерных тестовых расчетов
В этом пункте приводятся результаты численного моделирования процесса формирования, распространения и трансформации прерывной волны, возникающей в результате разрушения плотины и резкого локального подъема дна. Основное внимание уделяется расчетам, в которых прерывная волна распространяется по сухому руслу. Рассматриваются случаи горизонтального и наклонного дна, а также дна с локальным препятствием в нижнем бьефе. Характерная глубина Я выбирается равной начальной глубине верхнего бьефа в створе плотины. Численные расчеты проводятся по разностной схеме, описанной в пункте 4.1, в которой шаг по пространству А = 10 м, коэффициент запаса в условии устойчивости (4.43) Л = 0.5, коэффициент шероховатости русла г = 0.02 и параметр 5 = 0.05.
Поскольку введение искусственных вязкостей снижает точность разностной схемы на гладких решениях, то оптимальные значения их коэффициентов СІ выбираются как минимальные среди тех значений, которые обеспечивают отсутствие осцилляции за фронтами прерывных волн. При расчете начальной прерывной волны, возникающей после разрушения плотины и распространяющейся по горизонтальному или наклонному дну, такие оптимальные значения СІ получаются при условиях Сі = 0, С\ 0. Использование коэффициента Сч 0 необходимо при расчете водосливов, а также стоячих скачков или близких к ним прерывных волн, распространяющихся с малой скоростью. На таких разрывах модифицированный расход Q — q + SHu меняется незначительно, в силу чего вязкость (4.42) становится малой и слабо влияет на процесс подавления осцилляции.
На рис. 10, И кружками показаны результаты численных расчетов задачи о разрушении плотины над горизонтальным дном, выполненные по схеме, в которой С2 = 0, Сі 0. На рис. 10 и 11 расчеты приведены над горизонтальным дном в сравнении с соответствующими точными автомодельными решениями. Для разностного решения, показанного на рис. 10, коэффициент С\ = 0.85, а для решения, показанного на рис. И - С\ = 0.27. Поскольку во втором случае оптимальный коэффициент вязкости почти в три раза меньше, чем в первом случае, то на рис. 11 центрированная волна понижения воспроизводится разностным решением с гораздо более высокой точностью. Для сравнения на рис. 116 треугольниками изображено разностное решение, получаемое по схеме, аппроксимирующей классическую систему (3.1), (3.2) с коэффициентами искусственной вязкости Сі = 0.27, Сг = 0. Поскольку эта система не допускает распространения прерывных волн по сухому руслу, то в данном решении повышение уровня жидкости в головной части волны происходит гораздо более плавно и связано это повышение с действием донного трения (3.3), имеющего особенность при h — 0. Отметим, что, в отличие от этого решения, учет донного трения на столь малых временах визуально практически не влияет на профиль разностных решений, полученных путем аппроксимации модифицированной системы (3.1), (3.27) и показанных кружками на рис. 10 и 11.
На рис. 13 на четыре последовательных момента времени приведены результаты расчета задачи о трансформации прерывной волны на локальном береговом препятствии типа дамбы, расположенном в нижнем бьефе плотины. Прерывная волна возникает в результате решения задачи о разрушении плотины (3.34), в которой #о = 15 м, ho = 0. Сплошной линией показан профиль волны при Т = 50 с, когда она ещё не дошла до препятствия, штрих-пунктирной линией — при Т = 100 с, когда волна начинает натекать на него, штриховой линий — на момент времени Т = 200 с, когда волна уже перетекла через препятствие, и точками — при Т = 400 с, когда волна начинает распространятся по наклонному берегу за препятствием.
На рис. 13а приведены расчеты, в которых Сі = 0.3, Сг = 0. Такая вязкость эффективно подавляет осцилляции на фронте головной волны, однако не может их сгладить на фронте отраженной волны, в области водослива и на фронте вторичной волны за препятствием. На рис. 136 приведены расчеты, в которых Сі = 0.3, Сг = 0.5. Введение положительного коэффициента вязкости Сг позволяет подавить эти осцилляции.
Рассмотрим далее задачу о набегании прерывной волны на наклонный берег. Для решения данной задачи по модифицированной системе возьмем 6 — 0.3, что обеспечивает хорошее согласование с экспериментальными профилями прерывных волн, набегающих на наклонный берег [78].
