ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ УРАВНЕНИЙ МЕЛКОЙ ВОДЫ ПОВЫШЕННОЙ ТОЧНОСТИ НА ОСНОВЕ МОДИФИЦИРОВАННОЙ СХЕМЫ КАБАРЕ

Математические заметки СВФУ Июль—сентябрь, 2023. Том 30, № 3

УДК 519.6

ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ УРАВНЕНИЙ МЕЛКОЙ ВОДЫ ПОВЫШЕННОЙ ТОЧНОСТИ НА ОСНОВЕ МОДИФИЦИРОВАННОЙ СХЕМЫ КАБАРЕ Д. Г. Асфандияров, О. С. Сороковикова

Аннотация. Представлен численный метод на основе балансно-характеристической схемы КАБАРЕ для моделирования в приближении мелкой воды нестационарного течения жидкости на произвольной топографии. Разработанный метод позволяет рассчитывать различные режимы течения, включая транскритические. Для моделирования транскритических переходов используется гибридный подход на основе решения локальной задачи Римана, как это делается в схемах по типу Годунова. Представленный численный метод обладает условием хорошей сбалансированности (well-balance) — выполнение условия гидростатического равновесия или условия покоящейся жидкости на неровном рельефе дна. Помимо этого учитывается возможность сквозного расчета динамических границ, разделяющих жидкость и сухое дно, обусловленных процессами затопления и обмеления, а также ряда физических процессов (трение о дно, осадки). Апробация метода проводится на серии верификационных тестов, допускающих точное решение, и классическом эксперименте, имитирующем разрушении плотины.

DOI: 10.25587/SVFU.2023.33.65.009 Ключевые слова: уравнения мелкой воды, балансно-характеристический подход, схема КАБАРЕ, границы раздела вода/сухое дно, сложная топография дна и береговой линии.

1. Введение

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

© 2023 Асфандияров Д. Г., Сороковикова О. С.

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

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

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

В данной работе для расчета подобного рода задач используется идеология балансно-характеристической схемы КАБАРЕ [6-9]. Ранние работы по моделированию в рамках классического приближения уравнений мелкой воды по схеме КАБАРЕ [8, 9] ограничивались в основном расчетом одномерных тестов без учета практической составляющей, в том числе наличия границы раздела вода/сухое дно на неоднородном рельефе дна, что не позволяло решать большой класс задач (набегание прибрежный волн, наводнения, в том числе паводковые процессы, и т. д.). Данная работа призвана восполнить этот пробел.

В работе приводится обобщение методики КАБАРЕ с учетом требований

применимости к широкому диапазону режимов течения, наличия уреза воды, а также ряда физических процессов (трение о дно, осадки). Под различными режимами течения подразумевается возможность расчета как докритических и сверхкритических течений (аналог дозвукового и сверхзвукового течений в газовой динамике), а также транскритических переходов между этими двумя состояниями. С транскритическими переходами связана так называемая проблема звуковых точек (по аналогии с газовой динамикой). Для схемы КАБАРЕ разработан ряд алгоритмов учета звуковых точек, ограничивающихся частными случаями, например [9,10]. В данной работе используется более общий относительно класса решаемых задач гибридный подход обработки звуковых точек на основе решения локальной задачи Римана, как это делается в схемах по типу Годунова [1,2].

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

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

Нестационарная система уравнений теории мелкой воды, описывающая двумерные в плане движения жидкости, имеет вид

дК ^ дНи ^ дку дЬ дх ду >>>

дНи ^ дки2 ^ дкиу ^ 1 ддК2 ^ ^ЭЬ дЬ дх ду 2 дх дх

д1гу ^ дкуи ^ дНу1 ^ 1 ддК2 ^ ^дЬ дЬ дх ду 2 ду ду

Здесь ¡г = ¡(4, х, у) — глубина жидкости — расстояние от свободной поверхности до дна, м; 4 — время, с; (х, у) — декартовы координаты на горизонтальной плоскости, м; (и, у) — осредненные по глубине компоненты скорости в направлении х и у соответственно, м/с; К(х, у, 4) > 0 — интенсивность осадков, м/с; Ь = Ь(х, у) — отметка рельефа дна относительно произвольного горизонтального уровня, м; (/х, /у) — недифференциальные члены, определяющие действие силы трения в направлении х и у соответственно, м2/с2 Данные уравнения описывают законы сохранения массы жидкости и сохранения (изменения) ее импульса при условии постоянства плотности жидкости с учетом внешних сил. Сила трения задается формулой Маннинга — Штриклера [4]:

/=[/х,Л1Т = -^АУ|У|, У=К«]Т, \ = 2дпЧ-1/3,

где А — коэффициент гидравлического трения, п — коэффициент шероховатости дна по Маннингу,

Для учета свойства хорошей сбалансированности перепишем систему (1), сделав замену Ь = Н — Ь, где Н — уровень свободной поверхности, в последних двух членах второго и третьего уравнения. После преобразования получим

дЬ дЬи дЬу дЬ дх ду

0,

дНи ^ дНи2 ^ дкиу ^ д дН2 ^ дН дЬ дх ду 2 дх дх х'

д1гу ^ дкуи ^ Ну1 ^ д дН2 ^дН дЬ дх ду 2 ду ду у

Здесь мы не вводим каких-либо дополнительных переменных для моделирования, обозначение для уровня свободной поверхности Н = Ь + Ь представлено для наглядности преобразования. При выборе такой формы записи исходной системы уравнений для моделирования условия гидростатического равновесия выполняются автоматически и не требуют наложения дополнительных условий при аппроксимации рельефа дна [10].

Для построения численного алгоритма на основе балансно-характеристиче-ского подхода используется также характеристическая форма записи линеаризованной исходной системы уравнений (1) [7], которая имеет вид

д1у . д1у

дГ дГ дЬ х дх

С? х

д1 + ду

ССу.

Входящие в левую часть уравнений величины определяют квази-инварианты Римана (определенные вдоль соответствующих направлений) и соответствующие им характеристические скорости:

1у =

и + 2с

и — 2c

у

'у + 2сх у — 2с и

Xх =

ху =

XX

XX

XX.

'х?'

х^

х33.

и + с и— с и

'у + с '

у—с

у

Лх

Лу

хх 0 0

0 хх 0

0 0 хх

х? 0 0

0 х2 0

0 0 х

с = \Z~gh.

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

С? х

£2Х Ау —

Ф =

Ах —

ССу = Пу

и Ь 0 д и 0 0 0 и/

К

дх

Ау =

у0Ь 0 у 0 д0у

—ддЬ/ дх — /х/Ь — иК/Ь —ддЬ/ду — /у/Ь — уК/Ь;

х

Уравнения, входящие в систему (4), можно трактовать как неоднородные локально-одномерные уравнения переноса в направлениях х и у.

3. Численный алгоритм

3.1. Схема КАБАРЕ. Численный алгоритм строится на основе схемы КАБАРЕ. Схема КАБАРЕ оперирует двумя типами переменных: консервативными, относящимися к центрам ячеек, и потоковыми, которые относятся к центрам граней. Для расчета по схеме КАБАРЕ применяется так называемый балансно-характеристический подход, в рамках которого для вычисления консервативных переменных используются уравнения движения в консервативной форме, а для вычисления потоковых величин — в характеристической [6].

Введем в прямоугольнике [XI, X2] х [У1, У"2] неравномерную расчетную сетку с координатами узлов

:

х1 = XI, х1 < х2 <    < х№г + 1 = X2, Дхг+1/2 = хг+1 - хг

у1 = у1 < у2 <    < у^у + 1 = Ду^'+1/2 = у^+1 - уз.

Данный набор узлов порождает Мх х Му расчетных ячеек, Мх х (Жу + 1) продольных и (Жж + 1) х Ыу поперечных граней. Центрам ячеек соответствуют консервативные переменные, обозначаемые двумя полуцелыми индексами (г + 1/2, ] + 1/2), центрам горизонтальных и вертикальных граней соответствуют потоковые переменные, обозначаемые индексами (г + 1/2,и (г,^ + 1/2) соответственно. Схема оперирует тремя временными слоями: начальным п, полуцелым п + 1/2 и новым п + 1. Для шага по времени введем обозначение

тп+1/2 = _ щ

Схема КАБАРЕ для численного решения уравнений движения (3) (без учета силы трения) на конечно-разностной сетке, неоднородной по пространству и по времени, записывается следующим образом:

тт«+1 _ ттп 7=7 7=7 тт тт

*Л+1/2,3+1/2 1/2,3+1/2 + ^»+1,3 + 1/2 3 + 1/2 "»+1/2,3+1 "»+1/2,з _ ^

Тп+1'2 Лжг+1/2 А %+1/2 ~ '

и = [¡, и^, у ¡]Т, С = [¡и, ¡и2 + д(Н2 - 2ЬН)/2, ¡уи]Т, Н = [¡у, ¡иу, ¡у2 + д(Н2 - 2ЬН)/2]Т.

где И — вектор консервативных величин, С и Н — векторы потоков, относящиеся к поперечным и продольным граням соответственно. Верхний индекс обозначает номер временного слоя, а нижний — пространственное положение. Чертой сверху обозначается осреднение по двум временным слоям г] = (77™+1 +1]п)/2. Разностная схема (8) обладает вторым порядком аппроксимации как по времени, так и по пространству на неравномерных расчетных сетках [6, 7,12].

Один временной шаг схемы КАБАРЕ состоит из трех фаз. На первой фазе алгоритма вычисляются значения консервативных переменных на п + 1/2 временном слое. Во второй фазе с использованием значений промежуточных

консервативных величин и характеристической формы записи исходных дифференциальных уравнений вычисляются значения потоковых переменных на п + 1 временном слое. В третьей фазе по найденным на предыдущих шагах значениям консервативных и потоковых величин вычисляются значения консервативных переменных на п + 1 временном слое. При заданных начальных и граничных условиях по известным значениям консервативных и потоковых переменных на временном слое п однозначно определяются значения всех консервативных и потоковых переменных на слое п +1, после чего цикл, состоящий из трех описанных фаз, повторяется. Шаг по времени определяется из условия устойчивости Куранта — Фридрихса — Леви:

т "+1/2 = огь г™

г,з \ А^-+1/2 /( ^+1/2^+1/2 I + Сг+1/2,^+1/2; / Все расчеты в работе приводятся при числе Куранта СЕЬ = 0.3.

Рассмотрим алгоритм более подробно. Дискретизация уравнений (3) без учета силы трения в первой фазе (шаг предиктор) алгоритма имеет вид

т тп+1/2 _ т тп

^,+ 1/2^+1/2 тг+1/2,^+1/2

п+

1/2/2

рп _ рп нп _ ип

+ ^¿+1,^ + 1/2 ^г,3 + 1/2 + "¿+1/2,^+1 "¡+1/2,3 = 0 ^^

Ажг+1/2 1/2

Из выражения (10) получаются значения консервативных величин на полуцелом временном слое.

Во второй фазе используется характеристическая форма записи (4) исходных уравнений (1). Рассмотрим алгоритм нахождения потоковых переменных на новом временном слое по одному из направлений (вдоль оси ж). Вдоль другого направления алгоритм является полностью симметричным. Локальные инварианты Римана определяются как для потоковых, так и консервативных переменных. Будем использовать гидростатическую реконструкцию для глубины потока [13] при записи инвариантов для обеспечения свойства хорошей сбалансированности. В этом случае инварианты на используемом сеточном шаблоне (рис. 1) определяются так:

(7Г)д =< + 2\Ы/г?)' №)пч=<(7з)^=<'

(1гГт+1/2=<+1/2+№о:+1/2=<+1/2-

(И)Г/2 = <+1/2, ^п = ^ - (Ьо - ьд), д = -2,-1, 0,1, 2, т = -1,1.

Здесь индексы д и т определяют пространственное распределение вычисляемых величин на сеточном шаблоне согласно рис. 1.

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

Рис. 1. Сеточный шаблон для вычисления потоковых переменных.

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

т- <

где максимальные и минимальные значения определяются так:

т- = ш1п{(/Х)-2, (/Х)-1, (/Х)^-} + тп+1/2О-,

М- = шах{(/Х)2, (4Х)-1, (/Х)П-} + тп+1/2О-,

т+ = ш1п{

М+ = шах{(/Х)П, (/?)?, (4Х)П+} + тп+1/2О+, к = 1, 2, 3.

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

( т ж)п+1 / 2 _ (ГХ)И (ГХ)И _ (ГХ)П

— \п+1/2 \Тк )-1 \тк )-1 ( Х)п+1/2 \тк /0- \тк )-2

- тп+1/2/2 Дж_1

,+ \п+1/2 _ (/Х)П+1/2 — (/Х)П (\х)П+1/2 V к/2 V к )о+

к/1 - тп+1/2/2 д^

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

Локальные инварианты Римана на новом временном слое определяются следующим образом:

если [АГ+1/2 + [АХ]?+1/2 > 0, то (^Г = (4

иначе

Потоковые значения глубины и скорости на новом временном слое находятся из полученных инвариантов:

у \ /

<+1 = \т:+1+(фг1). = (^)Г1-

В области транскритического перехода (случай так называемых звуковых точек) выполняется [А^]— [А?]^2 < 0 или [А?]—+1/2  [А?]п+1/2 < 0, т. е. собственные значения справа и слева имеют разные знаки. В данном случае используется гибридный подход — значения потоковых переменных на новом временном слое находятся из решения задачи о распаде разрыва аналогично тому, как это делается в методах типа Годунова [1]. Начальные данные для задачи распада разрыва глубины и скорости потока и = [и, Н] слева и справа от грани (индексы — и + соответственно) находятся из линейной реконструкции глубины и скорости потока в ячейках слева и справа с использованием реконструкции-ограничителя тштс^ [1]:

и™— = и—1/2 + 0.59— Дхг—1/2, иП+ = и+1/2 — 0.5д+Дхг+1/2, 9- = 0.5  (sign(a) + sign(Ь))  тт(|а|, |Ь|), а = (иг+1/2 — иг_1/2)/(0.5  (Джг+1/2 + Джг_ 1/2)), Ь = (иг—1/2 — иг—3/2)/(0.5  (Джг—1/2 + Джг—3/2)), 9+ = 0.5  (sign(c) + sign(d))  тт(|с|, |Л|), С = (иг+3/2 — иг+1/2)/(0.5  (Джг+з/2 + Джг+1/2)), Л = (иг+1/2 — иг—1/2)/(0.5  (Джг+1/2 + Джг—1/2)).

Направление распространения инварианта /ЗХ определяется по знаку продольной компоненты скорости и, полученной в ходе решения задачи Римана на данной границе.

Третья фаза (шаг корректор) аналогична первой, по вычисленным на предыдущих этапах значениям консервативных и потоковых величин определяются значения консервативных переменных на (п + 1)-м временном слое:

тти+1 _ ттП+1/2 ри+1 _ ри+1

¿+1/2,3+1/2 иг+1/2,3+1/2 ^»+1,3+1/2 ^ г, 3+1/2

т«+1/2/2 + Ажг+1/2

НП+1 _ НП+1

+ 1 = 0.

ДУз+1/2

3.2. Учет трения о дно. Методом расщепления по процессам скорректированные значения консервативных переменных с учетом силы трения вычисляются так:

т

п+1 _

™+1 ¡и \"+11Т

п+1 (/ш)"+1 (Н

1 + те 1 + те

з = дп2{Нп+1) 4/3у/(И«+1)2 + (г;«+1)2.

3.3. Граничные условия. Схема допускает задание ряда физически значимых граничных условий таких, как вход, выход, стенка с проскальзыванием, стенка с прилипанием, периодические граничные условия. Реализация всех типов граничных условий в рамках схемы КАБАРЕ объединена общим подходом — дополнением приграничных граней недостающими значениями инвариантов с дальнейшим вычислением потоковых величин по общему алгоритму во второй фазе. Отличие всех граничных условий сводится к способу дополнения недостающих значений инвариантов [7, с. 196-202].

3.4. Алгоритм учета динамических границ вода/сухое дно. Во

многих практических задачах важным является определение уреза воды — линии пересечения водной поверхности бассейна реки или водоема с поверхностью суши. Данная граница может резко меняться во времени, например, при наводнении, в частности при паводковых процессах. Таким образом, значимым дополнением вычислительного метода является учет динамических границ вода/сухое дно. В работе за основу взят робастный алгоритм, подробно описанный в [2]. Критерий минимальной глубины применяется для классификации следующих четырех типов граней:

(1) «сухая» грань — в соседних ячейках выполняется условие Нь < е, Нд < е (рис. 2(а));

(2) «водная» грань — Нь > е, Нд > е (рис. 2(б));

(3) «частично смоченная» грань без потока, например, Нь > е, Нд < е, Нь < Нд (рис. 2(в));

(4) «частично смоченная» грань с потоком, например Нь > е, Нд < е, Нь > Нд (рис. 2(г)).

Исходя из данных типов граней, ячейки делятся на три типа:

(1) «водная» ячейка — грани ячейки принадлежат типу 2 или 4;

(2) «сухая» ячейка — грани ячейки принадлежат типу 1 или 3;

(3) «частично смоченная» ячейка — все остальные случаи.

«Сухие» ячейки исключаются из расчета, «водные» рассчитываются по базовому алгоритму, в «частично смоченных» ячейках учитывается только поток массы. В численных расчетах, представленных в данной работе, параметр е выбирается равным 10-6 м.

3.5. Выводы. В данном разделе был представлен вычислительный алгоритм для моделирования различных задач на основе уравнений мелкой воды. Основой численной методики является балансно-характеристическая схема КАБАРЕ. Алгоритм позволяет рассчитывать как докритические, так и сверхкритические режимы течения, равно как и транскритические переходы. Для решения проблемы звуковых точек (транскритических переходов) используется гибридный подход на основе решения локальной задачи Римана. Существенным дополнением вычислительного алгоритма является возможность учета динамических границ вода/сухое дно на неоднородном рельефе дна, что является определяющим для многих практически значимых задач. Добавлен алгоритм

нв= Ья

(а) сухая грань I я

(б) водная грань

1_ я

нй= Ьй HL

Нк= Ь„

(в) частично смоченная грань без потока (г) частично смоченная грань с потоком Рис. 2. Схематическое изображение границ раздела вода/сухое дно.

н

н

н

учета трения о дно. Представлен вариант учета условия хорошей сбалансированности схемы.

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

4. Результаты моделирования

Для верификации предложенного численного метода рассматривается ряд тестов, имеющих точное решение. Данные тесты включены в набор задач по аналитическим решениям [4] для верификации численных методик, основанных на приближении мелкой воды. Во всех тестах для оценки точности получаемых численных решений для невязки используется сеточный аналог нормы в Ь2: / N X 1/2

4.1. Вода в состоянии покоя на неровном рельефе дна. В данном тесте проверяется свойство хорошей сбалансированности схемы. Хорошо сбалансированные схемы должны удовлетворять условию отсутствия появления нефизичных перемещений воды над неровным дном из-за дискретизации исходных дифференциальных уравнений в частных производных. Рассматривается стационарная задача стоячей воды с неровным дном и границами раздела вода/сухое дно [14]. Расчетная область представляет собой закрытый квадратный бассейн размером 1м х 1м. В центре бассейна расположен симметричный холм, поверхность которого задается формулой Ь(ж, у) = тах[0, 0.25 - 4.8((ж - 0.5)2 + (у - 0.5)2)]. Трение о дно отсутствует, начальный

уровень в канале задается равным Н = 0.2м, и = V = 0м/с. Расчетное время задачи составляет Ь = 60 с. Расчеты проводились на равномерных сетках с количеством ячеек Nх х Му = 25 х 25 и 50 х 50. Топография дна и уровень поверхности воды на начальный момент времени представлены на рис. 3.

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

Полученные значения невязки в норме Ь2 для глубины воды Н и величин Ни, ^ представлены в табл. 1. Видно, что метод хорошо сохраняет стационарное решение.

Таблица 1. Невязки Н, Ни и Ни в норме Ь2. Задача о покоящейся воде на неровном рельефе дна

Nx х Ny V-

25 х 25 6.71 х Ю-13 7.92 х 10"15 7.92 х 10"15

50 х 50 4.81 х Ю-13 4.96 х 10"15 4.96 х 10"15

4.2. Задача Римана на горизонтальном дне. Рассматриваются классические одномерные тесты о распаде произвольного разрыва на горизонтальном дне, предложенные в работе [15], в квазидвумерной постановке. Начальные условия для тестов представлены в табл. 2, где Hl, ul, Hr и ur — начальные значения глубины и скорости слева и справа от разрыва соответственно, X0 — положение разрыва, tout — расчетное время. Длина расчетной области L = 50 м. Расчет проводится на двух сетках с количеством ячеек Nx = 500 и 1000 по оси x.

На рис. 4 представлены результаты расчета на подробной сетке (Nx = 1000). Численный метод дает хорошее совпадение с аналитическими решениями. В первом тесте размытие скачка происходит на 2-3 ячейки, что говорит о том, что скачок разрешается достаточно точно. Полученные значения невязки в норме L2 для глубины воды H для данных тестов на разных сетках представлены в табл. 3. Если говорить о порядке сходимости, то схема обеспечивает

Таблица 2. Данные тестовых задач о распаде произвольного разрыва

N Х0

1 10 1.0 2.5 0.1 0.0 7.0

2 25 1.0 -5.0 1.0 5.0 2.5

3 20 1.0 0.0 0.0 0.0 4.0

Таблица 3. Невязка Л в норме Ь2. Тесты о распаде произвольного разрыва

N Ь2(Н), Мх = 500 Ь2(Н), Мх = 1000

1 1.81 х 10"1 9.61 х 10~2

2 5.70 х 10~2 3.37 х 10~2

3 1.98 х 10~2 1.06 х 10~2

второй порядок только на гладких решениях. В данном случае решения не являются всюду гладкими (равно как и в остальных тестах), и в областях с большими градиентами моделируемых величин по теореме Годунова [1] могут возникать немонотонные осцилляции, поэтому в данных областях срабатывает ограничитель (13), который локально понижает порядок до первого. Исследование на порядок сходимости схемы КАБАРЕ на гладких решениях подробно представлено в статьях [6,16].

4.3. Тесты с обтеканием препятствия. Рассматриваются одномерные тесты в квазидвумерной постановке с обтеканием препятствия при различных режимах течения, представленные в [4]. Препятствие имеет вид параболического холма, поверхность которого задается формулой Ь(х) = max[0, 0.2 — 0.05(х — 10)2]. Длина канала Ь = 25м, трение о дно отсутствует. Как и во всех остальных одномерных тестовых примерах, размеры области по направлению у определяют только задаваемый расход и не влияют на результаты моделирования, равно как и количество ячеек по данному направлению. Начальные и граничные условия для тестов задаются следующим образом.

1. Докритический поток. Начальные условия: Н0 + Ь = 2м и и0 = у0 = 0м/е. В качестве граничных условий вверх по течению задается постоянный расход д = Ньиь = 4.42 м2/с (в 2Б постановке Нь(д) и соответственно иь находятся из аналитического решения), вниз по течению фиксируется глубина потока Ня = 2 м.

2. Транскритический поток без скачка. Начальные условия Н0 + Ь = 0.66 м и и0 = у0 = 0м/е. Граничные условия: вверх по течению д = Ньиь = 1.53 м2/с (Нь и иь определяются как и в предыдущем случае), вниз по течению фиксируется глубина потока Нд = 0.66 м пока поток докритический, затем устанав-

(а) левая критическая волна разряжения и правый скачок

1.2 1

0.8

¿ 0.6

0.4 0.2 0,

t = 2.5 с

, , , ,

10 20 30 40

X

50

б) две волны разряжения вблизи дна

1.2

0.8 Ж 0.6

X

0.4 0.2 0

0

10

20

X

30

t = 4 с - 1

i i i i i i i i . . . . . . . . I I I I

40

50

(в) распад разрыва на сухое дно Рис. 4. Задача о распаде разрыва. Уровень поверхности воды: 1 — точное решение, 2 — численное решение.

1

ливаются граничные условия на свободный выход.

3. Транскритический поток со скачком. Начальные условия Н0 + Ь = 0.33 м и и0 = V0 = 0м/с. Граничные условия: вверх по течению д = Н^и^ = 0.18 м2/с (Нь и и^ находятся аналогично предыдущим двум случаям), вниз по течению Нд = 0.33 м.

Расчеты проводятся на двух сетках с количеством ячеек Мх = 500 и 1000 по оси х. Для построения аналитических решений используется сетка с количеством ячеек Мх = 5000. Результаты расчета на подробной сетке (Мх = 1000) представлены на рис. 5, где для каждого из тестов также показана критическая высота 1гс = (д/у/д)2^3 для определения области докритического и сверхкритического потока. Результаты расчетов показывают очень хорошее совпадение с аналитическими решениями. В случае транскритического перехода со скачком скачок разрешается точно (рис. 5(г) — значения представлены во всех расчетных узлах). Для сравнения, выполнение данного теста в работе [4] приводит к размытию скачка на 4 ячейки. Полученные значения невязки в норме Ь2 для глубины потока Н для данных тестов на разных сетках представлены в табл. 4.

Таблица 4. Невязка Ь в норме Ь2.

Тесты с обтеканием препятствия

N Ь2(1г),Мх = 500 Ь2(Н), Мх = 1000

1 1.92 х 10"4 6.69 х 10~5

2 5.71 х 10"4 1.99 х 10~4

3 8.57 х 10~4 1.53 х 10~4

4.4. Тесты в каналах при наличии трения. Рассматриваются два теста из серии одномерных задач установившегося течения в открытых каналах при наличии трения, широко представленных в [4,5]. Особенность построения решения в этих задачах заключается в нахождении при заданном расходе воды в канале и высоты столба жидкости Н(х) в установившемся течении рельефа дна, которому это течение будет соответствовать. В установившемся течении с учетом параметризации гидравлического трения по формуле Маннинга — Штриклера (2) уклон дна определяется из системы уравнений (1) как

, (Ех + д)2\иЧ , 2(Пх + д)П , (Их + д\V

зд = -ахь=[ 1 -

Здесь д = Ни — задаваемый расход, м2/с. Сам уровень подстилающей поверхности Ь(х) можно получить численно, например, используя равномерную сетку с шагом Ах = Ь/М и соотношение Ь^+х = Ь + АхдхЬ, г = 0,... , N, Ь0 = Ь(Ь), где Ь — длина расчетной области, N — количество ячеек сетки. Затем на полученном рельефе дна воспроизводится модельное значение глубины потока.

2.5 2 1.5 X 1 0.5 0

У / 4 \ \

- 2 ---3 О 4

5 10 15 20 25 X

(а) докритический поток

1.2 1 0.8 — 0.6 X 0.4 0.2 °С

— 2 ---3 О 4

/ / / V

\

V.

5 10 15 20 25 X

(б) транскритический поток без скачка

(в) транскритический поток со скачком (г) область вблизи скачка

Рис. 5. Тесты с обтеканием препятствия. 1 — топография дна, 2 — аналитическое

решение, 3 — критическая высота, 4 — численное решение.

Тест 1. Рассматривается канал длиной Ь = 5000 м. Поток на входе и на выходе докритический. Глубина воды в установившемся состоянии задается функцией Нех(х) = 1.125 + 0.258ш(пх/500). Начальные условия: Н(х) = тах(0, Нех(5000) + 6(5000) — 6(х)), и0 = г>0 = 0 м/с. Граничные условия: вверх по течению д = Н^и^ = 2 м2/с, где Н^ = Нех(0), вниз по течению Нд = Нех(5000). Коэффициент шероховатости дна по Маннингу п = 0.03м-1/3с, интенсивность осадков К = 0 м/с. Расчет проводится на двух сетках с количеством ячеек Nх = 5000 и 10000 вдоль оси х. Для построения аналитического решения используется сетка с количеством ячеек N = 20000.

Тест 2. Рассматривается канал, течение в котором в установившемся состоянии является всюду докритическим. Длина канала Ь = 1000 м. Глубина воды в установившемся состоянии задается функцией

Нех(х) = (4/д)1/3(1 + 0.5ехр(—16(0.001х — 0.5)2)). Начальные условия: Н(х) = 0 м, и0 = у0 = 0 м/с. Граничные условия: вверх по

(а) глубина воды — периодическая (б) тест при наличии осадков

функция в пространстве

Рис. 6. Тесты в каналах при наличии трения: 1 — топография дна, 2 — аналитическое решение, 3 — численное решение.

течению д = Н^и^ = 2м2/с, где Н^ = Нех(0), вниз по течению Нд = /гех(1000). Коэффициент шероховатости дна по Маннингу п = 0.033 м-1/3с, интенсивность осадков К = 0.001м/с. Расчет проводится на двух сетках с количеством ячеек Ых = 1000 и 2000 вдоль оси х. Для построения аналитического решения используется сетка с количеством ячеек Ых = 4000.

Результаты расчета на подробных сетках для двух тестов представлены на рис. 6. Результаты расчетов показывают хорошее совпадение с аналитическими решениями. Метод хорошо справляется с границами раздела вода/сухое дно на неровной поверхности и при наличии трения. Полученные значения невязки в норме Ь2 для глубины потока Н для данных тестов на разных сетках представлены в табл. 5.

Таблица 5. Невязка Л в норме Ь2.

Тесты в каналах при наличии трения

N Ь2(/г), Жх = 1000 Ь2(/г), Жх = 2000 Ь2(/г),Жх = 5000 Ь2(/г),Жх = 10000

1 — — 1.07 х 10° 1.07 х 10°

2 4.61 х Ю-1 4.60 х Ю-1 — —

4.5. Задача о разрушении двумерной плотины. Рассматривается ва-лидационная задача о прорыве плотины [11]. Расчетная область представляет собой прямоугольник размером 3 м х 2 м. Слева находится резервуар размером 1м х 2 м с непротекаемыми стенками. На правой границе резервуара имеется тонкая перегородка с симметричным относительно оси абсцисс отверстием шириной 0,40 м (рис. 7(а)). В начальный момент времени отверстие в перегородке закрыто и вся жидкость располагается слева от нее. Все условия проведенного

Шя

НЮЪ

Ш2 ШЗ

X

(а) схематичное изображение расчетной области и расположение контрольных точек

Рис. 7. Задача о разрушении двумерной плотины.

(б) 3Б-визуализация моделируемого течения на момент времени £ = 0.837 с

численного эксперимента взяты в соответствии с условиями лабораторного и численного экспериментов, представленных в [11]. Глубина воды в резервуаре в начальный момент времени равна 0.6 м, уклон дна составляет 0%. Трение о дно не учитывается. Количество ячеек расчетной сетки Ых х Ыу = 150 х 50. Расчетное время составляет 10 с.

Динамика глубины и скорости потока замеряется в контрольных точках, расположение которых показано на рис. 7(а). На рис. 7(Ь) приведена 3Б-визуа-лизация моделируемого течения (глубины потока) на момент времени Ь = 0.837 с.

Результаты расчета в контрольных точках глубины (точки Р0, Р1, Р2, Р3) и продольной компоненты скорости потока (точки Р0 и Р4) представлены на рис. 8. Также приводится сравнение с результатами лабораторного эксперимента и расчета по авторской численной методике на основе WAF-метода [11]. Экспериментальные данные для продольной компоненты скорости представлены осредненными по глубине.

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

5. Заключение

1

1.0

1.0

1.322

1.722

0.58

1.0

1.16

1.0

1.0

1.2

О

В работе представлен численный метод на основе модифицированной схемы КАБАРЕ для моделирования нестационарных течений в приближении мелкой

0.6 0.5 0.4 - 0.3 0.2 0.1 0,

4 6 8

t -С)

(а) глубина потока в точке Р0

1.6 1.4 1.2 1

¿0.8 0.6 0.4 0.2 0,

- 2

--3

t -С)

(б) продольная компонента скорости в точке Р0

0.6 0.5 0.4 - 0.3 0.2 0.1 0,

---1 - 2

Л ~ - \ - - 3

- Х7\\

0 2 4 6 8 10

t -С)

(в) глубина потока в точке Р1

0.35 0.3 0.25 — 0.2 Х 0.15 0.1 0.05 0,

— 1

11 ^ 1' \ - - 3

^ \ ! \

\ \

\ Ч

0

t -С)

(г) глубина потока в точке Р2

(д) глубина потока в точке Р3 (е) продольная компонента

скорости в точке Р4

Рис. 8. Задача о разрушении двумерной плотины. 1 — данные эксперимента [11], 2 — результаты моделирования, 3 — результаты моделирования по ^АЕ-методу [11].

2

4

6

8

2

4

6

8

воды.

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

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

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

ЛИТЕРАТУРА

1. Куликовский А. Г., Погорелов Н. В., Семенов А. Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001.

2. Huang Y., Zhang N., and Pei Y. Well-balanced finite volume scheme for shallow water flooding and drying over arbitrary topography // Eng. Appl. Comput. Fluid Mech. 2013. V. 7, N 1. P. 40-54.

3. Bermudez A., Dervieux J. A., Vázquez M. E. Upwind schemes for the two dimensional shallow water equations with variable depth using unstructured meshes // Comput. Methods Appl. Mech. Eng. 1998. V. 155. P. 49-72.

4. Delestre O., Lucas C., Ksinant P.-A., Darboux F., Laguerre C., et al. SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies // Int. J. Numer. Methods Fluids. Wiley, 2013. V. 72, N 3. P. 269-300.

5. Ngoc Tuoi Vo Thi. One dimensional Saint-Venant system. Anal. PDEs [math.AP]. 2008. https://dumas.ccsd.cnrs.fr/dumas-00597434.

6. Karabasov S., Goloviznin V. Compact accurately boundary-adjusting high-resolution technique for fluid dynamics // J. Comput. Phys. 2009. V. 228. P. 7426-7451.

7. Головизнин В. М., Зайцев М. А., Карабасов С. А., Короткин И. А. Новые алгоритмы вычислительной гидродинамики для многопроцессорных вычислительных комплексов. М.: Изд-во Моск. ун-та, 2013.

8. Goloviznin V. M., Isakov V. A. Balance-characteristic scheme as applied to the shallow water equations over a rough bottom // Comput. Math. Math. Phys. 2017. V. 57. P. 1140-1157.

9. Afanasiev N., Goloviznin V. A locally implicit time-reversible sonic point processing algorithm for one-dimensional shallow-water equations //J. Comput. Phys. 2021. V. 434. Paper ID

110220.

10. Головизнин В. М., Соловьев А. В., Исаков В. А. Аппроксимационной алгоритм обработки звуковых точек в схеме «КАБАРЕ» // Вычисл. методы и программирование. 2016. Т. 17, вып. 2. С. 166-176.

11. Fraccarollo L., Toro E. F. Experimental and numerical assessment of the shallow water model for two-dimensional dam-break type problems //J. Hydraul. Res. 1995. V. 33, N 6. P. 843864.

12. Kostrykin S. V. A variant of the multidimensional generalization of the cabaret scheme // Math. Models Comput. Simul. 2018. V. 2. P. 564-573.

13. Delestre O., Cordier S., Darboux F., Du M., James F., Laguerre C., Lucas C., and Planchon O. FullSWOF: A free software package for the simulation of shallow water flows. 2014. arXiv:1401.4125.

14. Brufau P., Garcia-Navarro P. Unsteady free surface flow simulation over complex topography with a multidimensional upwind technique //J. Comput. Phys. 2003. V. 186, N 2. P. 503-526.

15. Toro E. T. Shock-capturing methods for free-surface shallow flows. Chichester: Wiley, 2001.

16. Chintagunta A., Naghibi S. E., Karabasov S. A. Flux-corrected dispersion improved CABARET schemes for linear and nonlinear wave propagation problems // Comput. Fluids. 2017.

Поступила в редакцию 1 августа 2023 г. После доработки 1 сентября 2023 г. Принята к публикации 4 сентября 2023 г.

Асфандияров Данил Гамилевич, Сороковикова Ольга Спартаковна Институт проблем безопасного развития атомной энергетики РАН, ул. Большая Тульская, 52, Москва 115191 dasfandiyarov@ibrae.ас.ru, olga_sorokov@mail. ru

Математические заметки СВФУ Июль—сентябрь, 2023. Том 30, № 3

UDC 519.6

A HIGH-RESOLUTION NUMERICAL METHOD FOR SOLVING THE SHALLOW WATER EQUATIONS BASED ON THE MODIFIED CABARET SCHEME D. G. Asfandiyarov and O. S. Sorokovikova

Abstract: A numerical method based on the CABARET scheme for modeling unsteady flow over arbitrary topography in the shallow water approximation is developed. The method allows simulating a wide range of flow conditions, including transcritical. To model transcritical transitions, a hybrid approach is used based on solving the local Riemann problem, as is done in Godunov-type schemes. The presented numerical method has a well-balance condition—the fulfillment of the condition of hydrostatic equilibrium or the condition of a fluid at rest on an uneven bottom topography. A robust technique is used to simulate the movement of wet/dry fronts caused by flooding or recession. A number of physical processes are taken into account, such as bed friction and rain. Numerical results are compared with analytical solutions and data from the dam-break experiment.

DOI: 10.25587/SVFU.2023.33.65.009

Keywords: shallow water equations, balance-characteristic approach, wet/dry fronts, CABARET scheme, complex topography.

REFERENCES

1. Kulikovsky A. G., Pogorelov N. V., and Semenov A. Yu., Mathematical Questions of the Numerical Solution of Hyperbolic Systems of Equations [in Russian], FIZMATLIT, Moscow (2001).

2. Huang Y., Zhang N., and Pei Y., "Well-balanced finite volume scheme for shallow water flooding and drying over arbitrary topography," Eng. Appl. Comput. Fluid Mech., 7, No. 1, 40-54 (2013).

3. Bermudez A., Dervieux J. A., and Vázquez M. E., "Upwind schemes for the two dimensional shallow water equations with variable depth using unstructured meshes," Comput. Methods Appl. Mech. Eng., 155, 49-72 (1998).

4. Delestre O., Lucas C., Ksinant P.-A., Darboux F., Laguerre C. et al., "SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies," Int. J. Numer. Methods Fluids, 72, No. 3, 269-300 (2013).

5. Ngoc Tuoi Vo Thi, One dimensional Saint-Venant system, Analysis of PDEs [math.AP]. 2008. https: //dumas.ccsd.cnrs.fr/dumas-00597434

6. Karabasov S. and Goloviznin V., "Compact accurately boundary-adjusting high-resolution technique for fluid dynamics," J. Comput. Phys., 228, 7426-7451 (2009).

7. Goloviznin V. M., Zaitsev M. A., Karabasov S. A., and Korotkin I. A., New Algorithms of Computational Fluid Dynamics for Multiprocessor Computing Systems [in Russian], Izdat. Moskov. Univ., Moscow (2013).

8. Goloviznin V. M. and Isakov V. A., "Balance-characteristic scheme as applied to the shallow water equations over a rough bottom," Comput. Math. Math. Phys., 57, 1140-1157 (2017).

© 2023 D. G. Asfandiyarov, O. S. Sorokovikova

9. Afanasiev N. and Goloviznin V., "A locally implicit time-reversible sonic point processing algorithm for one-dimensional shallow-water equations," J. Comput. Phys., 434, paper ID 110220 (2021).

10. Goloviznin V. M., Solovjev A. V., and Isakov V. A., "An approximation algorithm for the treatment of sound points in the CABARET scheme [in Russian]," Numer. Methods Program., 17, 166-176 (2016).

11. Fraccarollo L. and Toro E. F., "Experimental and numerical assessment of the shallow water model for two-dimensional dam-break type problems," J. Hydraul. Res., 33, No. 6, 843-863 (1995).

12. Kostrykin S. V., "A variant of the multidimensional generalization of the CABARET scheme," Math. Models Comput. Simul., 2, 564-573 (2010).

13. Delestre O., Cordier S., Darboux F., Du M., James F., Laguerre C., Lucas C., and Planchon O., FullSWOF: A free software package for the simulation of shallow water flows, arXiv:1401.4125 (2014).

14. Brufau P. and García-Navarro P., "Unsteady free surface flow simulation over complex topography with a multidimensional upwind technique," J. Comput. Phys., 186, No. 2, 503-526 (2003).

15. Toro E. F., Shock-Capturing Methods for Free-Surface Shallow Flows, Wiley, Chichester (2001).

16. Chintagunta A., Naghibi S. E., and Karabasov S. A., "Flux-corrected dispersion improved CABARET schemes for linear and nonlinear wave propagation problems," Comput. Fluids (2017).

Submitted August 1, 2023 Revised September 1, 2023 Accepted September 4, 2023

Danil G. Asfandiyarov, Olga S. Sorokovikova Nuclear Safety Institute of Russian Academy of Sciences, 52 Bolshaya Tulskaya Street, 115191 Moscow, Russia dasfandiyarov@ibrae.ac.ru, olga_sorokov@mail.ru