CC BY
40
Вычислительные технологии
Том 1, № 2, 1996
НЕКОТОРЫЕ ПРОБЛЕМЫ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ВОЛНОВЫХ РЕЖИМОВ В ОГРАЖДЕННЫХ АКВАТОРИЯХ^
В. Б. Барахнин, Г. С. Хлкимзянов, Л. Б. Чубаров, Д. А. Шкуропацкий Институт вычислительных технологий СО РАН Новосибирск, Россия
Рассматриваются некоторые аспекты построения компьютерной модели волнового режима огражденной акватории на примере верхнего водоема ГАЭС (гидроаккуму-лирующей электростанции). Вычисление общего плана течения — уровня свободной поверхности и поля скоростей — осуществляется в ходе сценарных расчетов, воспроизводящих различные режимы работы ГАЭС с учетом варьирования основных параметров энергетических установок, изменения геометрии водоема, в том числе предусматривающего неодновременный ввод объекта в эксплуатацию, и связанного с этим сооружения временной дамбы, а также ее размыв, устранение. Изложенные в статье результаты могут быть использованы также и при исследовании гидрофизических процессов, связанных с проникновением длинных океанических волн цунами в бухты, гавани и т.п.
1. Введение
Верхние водоемы гидроаккумулирующих электростанций имеют два характерных режима работы: насосный (наполнение водоема в период спада энергопотребления) и турбинный (опорожнение, срабатывание водоема в период пикового энергопотребления). Таким образом, специфической особенностью водохранилищ ГАЭС (по сравнению с водохранилищами ГЭС) является более напряженный режим их использования вследствие регулярного периодического изменения основных гидрофизических полей (уровень поверхности и скорость течения) при переключении режимов. Поэтому как на этапе проектирования объекта, так и на этапах его строительства и эксплуатации крайне важным является понимание динамики процессов, происходящих в водоеме, и, в частности, их критических режимов с возможными катастрофическими последствиями: переполнение водоема и разрушение ограждающих дамб, размыв дна, осушение зоны водозаборных сооружений и т.п.
*© В. Б. Барахнин, Г. С. Хакимзянов, Л. Б. Чубаров, Д. А. Шкуропацкий, 1996.
^ Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №94-05-16281.
Основной задачей рассматриваемой компьютерной модели является определение общего плана течения — уровня свободной поверхности и поля скоростей. Вычисление указанных полей осуществляется в ходе сценарных расчетов, воспроизводящих различные режимы работы ГАЭС с учетом варьирования основных параметров энергетических установок, а также изменения геометрии водоема, связанного с выбором в ходе проектирования оптимальных форм ограждающей дамбы и рельефа дна.
Прообразом изучаемого объекта является водоем Днестровской ГАЭС, сооружение которого выполняется из экономических соображений в два этапа и позволит сразу после возведения первой очереди водоема и временной дамбы ввести в эксплуатацию энергетические агрегаты в режиме, близком к проектному. Поэтому перед проектировщиками встает также проблема устранения такой разделяющей дамбы после полного завершения строительства водоема. При этом рассматриваются различные варианты — принудительный размыв дамбы без остановки агрегатов, взрывная или механическая разборка с временным прекращением работы станции и осушением водохранилища. Главным критерием выбора варианта являются соображения безопасности процесса.
С точки зрения математического моделирования, первый вариант может быть рассмотрен в рамках классической задачи о разрушении плотины. Он представляет несомненный интерес для проектировщиков и должен быть обеспечен соответствующими средствами исследования с помощью создаваемой компьютерной модели.
Задача построения такой компьютерной системы в рамках концепции вычислительного эксперимента требует обоснованного выбора математической модели, конструирования адекватных вычислительных алгоритмов, обоснования и реализации структуры программного обеспечения, включающего сервисные и вычислительные модули. Обсуждение указанных аспектов вычислительного эксперимента в связи с конкретным кругом задач и составляет содержание настоящей статьи.
2. Математические модели
Выбор математических моделей для описания конкретного физического явления требует предварительного анализа и учета характерных параметров этого явления и сравнительной оценки их значимости. Так, в частности, геометрия области течения верхнего водоема ГАЭС (рис. 1) характеризуется значительным преобладанием горизонтальных размеров по сравнению с вертикальным (Нтах & 20 м, Ьх & 3000 м, Ьу & 1000 м)1, вытянутостью водоема вдоль одного из горизонтальных направлений, малыми изменениями поля глубин. Другим определяющим параметром является характерное время процесса, которое задается продолжительностью режимов наполнения и опорожнения — около 5 часов (рис. 2). Отметим также вопрос о важности учета геометрии ограждающей дамбы (профиля ее сечения) на динамику уреза свободной поверхности.
Влияние характеристик подстилающей поверхности (материалов, ее составляющих), которые определяют эффекты донного трения, и образующегося в зимний период ледового покрова не столь очевидно.
Упомянутое выше соотношение характерных размеров водоема делает возможным изучение происходящих в нем гидрофизических процессов в терминах теории мелкой воды (возможно, и в простейшей линейной постановке — параметр нелинейности а/Н не пре-
1 Здесь и далее используются характерные размеры Днестровской ГАЭС [1].
Рис. 1. Геометрия области течения верхнего водоема ГАЭС.
вышает значения 0.05), уравнения которой имеют вид
ut + + VUy + gnx = fi, Vt + uVx + vvy + д'Пу = f2,
nt + (uh)x + (vh)y = fз,
(1)
где и, V — компоненты "горизонтальной" скорости, к = Н+п — полная глубина, Н — глубина невозмущенного слоя жидкости, п — смещение свободной поверхности, д — ускорение силы тяжести, правые части / в уравнениях движения содержат члены, описывающие воздействие внешних факторов (сил Кориолиса, донного и ветрового трения):
fi = lv - д
uVu2 + V2
C 2h
+ V (uxx + uyy) + f(x), f2 = lu -
vVu2 + v2
д
C2 h
+ V (vxx + vyy) + f(y),
l = sin ^ — угловая скорость вращения Земли, ^ — географическая широта, C — ко-
R1/6
эффициент Шези, определяемый из соотношения C =- (формула Маннинга), R —
n
гидравлический радиус, обычно принимаемый для широкого открытого русла равным глубине H, n — коэффициент групповой шероховатости, v — коэффициент динамической вязкости, f(x), f(y) — слагаемые, учитывающие эффекты ветрового трения на свободной поверхности, f3 — правая часть в уравнении неразрывности, отвечающая за внешние источники массы.
Укажем также на другую форму записи этих уравнений, которая оказывается предпочтительнее при решении задач, связанных с расчетом разрывных полей ("прорыв плотины"):
Wt + Fx + Gy = E, (2)
w frrT,n t, í U2 h2 UV TT\ „ (UV V2 h2 где W = {U,V,h}, F = |T + gy,—,Uj, G = {~,T +
í U (U2 + V2) V (U2 + V2) i
E = <¡ ghHx + IV - g V ; + Ai, ghHy - lU - g V ; + A2, 0
C2 h2
Ah
C 2h2
U
A = {Ai,A2> = vAU - vU— - v V— -Vh, U = {U, V}, U = hu, V = hv.
h
h
суббота воскресенье
Рис. 2. Характерный график режимов наполнения и опорожнения водоема.
Боковая граница Г(^) области течения Q(t) состоит из двух фрагментов, один из которых Г1 (^ соответствует "водозабору", и на этом фрагменте задается расход Q(t), зависящий от времени (от числа работающих агрегатов и режима их работы), график изменения которого приведен на рис. 3.
Рис. 3. Типичный график изменения расхода.
Другой фрагмент границы Г2^) соответствует ограждающей дамбе и, в зависимости от детальности моделирования, описывается либо вертикальной, либо наклонной (с постоянным углом наклона а) поверхностью, тогда и можно сказать, что задача решается в области с подвижной границей, перемещающейся со скоростью, равной скорости частиц жидкости вблизи нее, и на этой границе формулируется условие равенства нулю полной глубины. Известны различные способы реализации таких граничных условий (см., например, [1]—[3]). Некоторые из них будут рассмотрены ниже при обсуждении численных алгоритмов.
Заметим также, что упомянутая выше вытянутость водоема позволяет на предварительном этапе исследований проводить расчеты по одномерным моделям теории мелкой воды [4], а также с использованием хорошо известной в гидравлике модели Сен-Венана. Уравнения Сен-Венана занимают более высокую ступень в иерархии приближенных гид-
родинамических (гидравлических) моделей и позволяют в одномерной постановке учесть эффекты, обусловленные изменением ширины русла. По мнению авторов, эти простые, в том числе и по своей реализации, модели также должны быть включены в компьютерную модель.
Подробный вывод уравнений модели Сен-Венана и дискуссия о связанных с ней допущениях и упрощениях приведены в монографиях М. Б. Эббота [5] и М. С. Грушевского [6]. Отметим здесь главные из допущений:
1) поперечные составляющие скорости малы по сравнению с продольной, центробежный эффект, обусловленный кривизной русла, не учитывается;
2) уклон дна мал;
3) силы сопротивления вводятся в уравнения в том же виде, что и для равномерного движения; предполагается, что суммарное влияние сил трения и турбулентности можно учесть включением в уравнение некоторой силы сопротивления.
При таких предположениях уравнения Сен-Венана имеют вид
^ + д ( Q2N + дп + Q|Q| 0
Иъ + дх{а-л) + дЛдХ + 9с2лп — щ е°8 Р = 0'
дл + дЯ =
дъ дх
здесь Л = к ¿у — площадь полного сечения, Q = к й(1у = иЛ — расход, к — локальное значение полной глубины, м — осредненная по глубине скорость, у — поперечная горизонтальная координата, В — полная ширина живого сечения русла, д — расход боковой приточности, направленной под углом р к линиям тока, V — ее средняя скорость, Л2
а = —- и?в,Л — коэффициент распределения скорости для всего сечения (коэффици-^ J А
Л
ент Буссинеска, обычно принимается равным единице), ' = т^ — Н, уточненная формула
В
для вычисления гидравлического радиуса К имеет вид К / /(к )3/2ё,зу, где / / —
л Л)
отношение коэффициентов Шези в рассматриваемом сечении и в главном русле. Как уже было отмечено выше, обычно принимают К = Н.
В более полной постановке задача может рассматриваться в терминах трехмерного потенциального течения идеальной жидкости со свободной границей, которая в произвольной системе координат сводится к определению потенциала скорости р и функции п(д1,д1'Ъ), описывающей свободную поверхность, удовлетворяющих нелинейному уравнению эллиптического типа в единичном кубе Q
£ ^ (*-=о. д = (д1.д2-д3) € а, (4)
нелинейным краевым условиям на границе Г = Г0 и Г1 (Г0 — свободная поверхность, Г1 — твердая стенка или участки границы, через которые происходит наполнение или опорожнение водоема)
V 2
Пг + v1пql + v2пq2 = из- рг = Vо — 1 — — — дп- д € Го-
2^кав ^в в=1 4
= ? е Г1,
qa=const
начальным условиям при £ = 0
0) = ^(?), 5 е <3, п(?1 ,52, 0) = п(?\?2), ? = (?\?2,1) е Го,
где иа и уа — декартовы и контравариантные компоненты вектора скорости V, V = |У|, < — обезразмеренный расход жидкости на участке водозабора,
иа = (д7в^в), = ^ + дав^, кав = ,
7 — якобиан преобразования £ = ха = жа(дв,£), а, в = 1, 2, 3, дав — элементы матрицы, обратной к матрице {дав} ковариантных компонент метрического тензора дав = 5ж7 5ж7
/ т;--тп;, V« — ковариантные компоненты вектора скорости,
^=1 5дв
^в = ф,«о = 1 + -^и«, -д- = (-1)а ае^д^^^ / ^ = 0,1, 2, 3, ^ = а.
Здесь всюду предполагается, что по повторяющимся индексам а, в, 7 = 1, 2, 3 производится суммирование.
Использование трехмерной модели позволяет учесть изменение параметров течения по глубине и определить распределение скорости на любом заданном горизонте. Включение этой компоненты в компьютерную модель значительно повышает уровень достоверности результатов и оказывается необходимым для верификации приближенных гидродинамических моделей. Такое решение должно быть принято, несмотря на очевидный рост объема вычислительных затрат.
3. Вычислительные алгоритмы
Одним из принципиальных вопросов является обеспечение быстродействия алгоритмов наряду с их высокой точностью, позволяющее адекватно воспроизводить медленные процессы в течение достаточно продолжительного времени. В связи с этим обстоятельством для аппроксимации уравнений (1) и их модификаций предлагается использовать неявные (полунеявные) конечно-разностные схемы на подвижных адаптивных сетках. В то же время для моделирования быстропротекающих процессов (упомянутая выше задача об удалении временной дамбы) будут применяться явные схемы на различных типах сеток (равномерных и неравномерных, подвижных и неподвижных криволинейных сетках).
Численные алгоритмы, предназначенные для реализации описанных выше моделей, в связи с обсуждением соответствующих модельных расчетов будут рассмотрены в основном в одномерном варианте. Некоторые аспекты двумерной реализации этих алгоритмов затронуты вкратце, а описание соответствующих методов расчета трехмерных течений приведено в работе [7].
3.1. Конечно-разностные схемы
Начнем рассмотрение с одномерных уравнений теории мелкой воды в их простейшем варианте, т.е. без учета дополнительных факторов:
и
к2
Wг + Ех = Е, W = {и-к}' Е = < — + д—^}- Е = {дкНх-0}- И = ки
к
2
(5)
Одной из наиболее эффективных конечно-разностных схем для этой модели является явный аналог схемы предиктор-корректор Годунова [8] с автоматически настраиваемой аппроксимационной вязкостью. Заметим, что на предварительном шаге (предиктор) аппроксимируются уравнения (5), записанные в недивергентной форме:
Рг + БРх = Т- Р
и к
Б
ид ки
Т
дНх 0
На неподвижной равномерной сетке эта схема будет иметь следующий вид. Предиктор:
+
2
Корректор:
Р,
"рга I "рга Р,+1/2 + Р 3 — 1/2
3™
Г»™
+ Б—1/2
рп
Р,+1/2
рп
Р,-1/2
Ах
+
т™
т з-
М- 1.
(6)
Wn+1, - Wn , Е Е ' ,+1/2 ' ,+1/2 + Е 3+1 — Е 3
тп
Ах
= Е3+1/2- 3
1-
М- 1.
Помимо этой схемы были рассмотрены также неявная схема предиктор-корректор Годунова [8], требующая для своей реализации применения матричной прогонки, схема Мак-Кормака в явном и неявном вариантах [9]. Причем неявная схема Мак-Кормака оказывается достаточно простой в реализации и основана на явной схеме, однако при локальном нарушении условия устойчивости в нее вводится неявная стабилизирующая добавка, рассчитываемая методом "бегущего счета".
Что касается алгоритмов для модели Сен-Венана, то для их конструирования используются те же принципы, что и для уравнений теории мелкой воды. Также, как и в предыдущем случае, конечно-разностная аппроксимация будет основана на упрощенном варианте исходной дифференциальной системы (3) без учета внешних сил и источников.
Предиктор:
Хг + У^Хх — ^ - X
Q
Л
у
2Q 1
Q л 0
Z
—длпх 0
здесь
X, — X™
3 т
+ У™ ~ з
■у"п _ "V"™
п 3+1/2 З—1/2
Ах
/,
2 (/3+1/2 + /3—1/2)-
Ах
3 = 2-
М- 1.
(7)
2
1
2
з
2
Корректор:
Х + Ь = ^ Ь =
<2 А
уп+1 _ ХП т т
Х.7+1/2 Х,7+1/2 -Ц+1 - Ь
+
Аж
= Z,■
5+1/2,
здесь
Z
-дА^+1/2 ^
'¿+1/2
0
^¿+1/2 = 2(А+1 + А), ^ = 1, ..., N - 1.
Двумерный вариант схемы (6) подробно рассмотрен в работе [10]. Так же, как и в одномерном случае, величины и, V, Н вычисляются на шаге корректор в центрах ячеек, а на шаге предиктор — на боковых сторонах этих ячеек. Управляющий "сдвигом по времени" параметр и определяется градиентом свободной поверхности Уп и обеспечивает увеличение схемной вязкости в областях резкого изменения уровня.
Применение двумерных неявных схем потребует, естественно, использования подходов, связанных с тем или иным методом расщепления на последовательность "одномерных" задач с сохранением преимуществ простоты реализации.
3.2. Расчетные сетки
Остановимся на способах построения расчетных сеток. На рис. 4 показано дно водохранилища в плане, сетка, лежащая на дне, а также контур ограждающей дамбы.
Рис. 4. Криволинейная сетка и контур ограждающей дамбы.
Координаты сетки определяются в результате итерационного решения двух нелинейных уравнений
д ( джа\ д ( джа\
а? Гд22 ^ + д? ~эФ) = «• а =12 (8)
где ю — управляющая функция, которая может зависеть от глубины Н, возвышения свободной поверхности п и его градиента Уп,
дж дж ду ду даа = "ЁПЮТГЮ +
дда дда дда дд0
Построенная сетка является криволинейной, адаптирующейся к границам бассейна, квазиортогональной. Выбирая управляющую функцию ю(ж,у) в виде ю = 1 + |Уп|, можно обеспечить движущееся сгущение сетки в окрестности быстрого изменения свободной поверхности и разрежение сетки в областях малых возмущений. Такое избирательное измельчение сетки повышает точность расчетов и позволяет получать адекватные результаты на сравнительно небольшом количестве узлов. Например, в тестовых одномерных расчетах для получения результатов одинаковой точности на неравномерной сетке требуется на порядок меньшее количество узлов, чем на равномерной.
При решении задач в двумерной постановке с учетом наклона ограждающей дамбы расчетная сетка на каждом временном шаге строится для области, ограниченной линией уреза. Такая сетка становится подвижной, т.е. перемещаемой в соответствии с движением границы области. Для адекватного описания граничных условий, связанных с распределенным расходом (в зависимости от числа работающих агрегатов N), количество узлов сетки на такой границе выбирается кратным N.
В одномерном случае сетки строятся на основе уравнения, аналогичного (8). В трехмерной постановке сетка на дне строится так же, как и в плановой (двумерной) задаче, а далее из каждого узла, лежащего на дне, проводится наклонная прямая до пересечения со свободной поверхностью. Наклоны прямых подстраиваются под наклоны берегов, и узлы сетки движутся вдоль этих прямых.
3.3. Граничные условия
При расчете одномерных течений для случая вертикальной стенки граничные значения
0 и следствия дифференци-
определяются на основе заданного краевого условия и
альных уравнений (5) пх
0.
В случае наклонного берега точка уреза ж = ¿(¿) движется вдоль стенки по упомя-^ж
нутому выше закону — = и .В подвижной системе координат движущейся точке
ж = ¿(¿) соответствует точка д =1 неподвижной системы координат, в которой и вычисляются граничные значения. Для их определения здесь привлекается заданное граничное условие дифференциальной задачи Н = 0 (Н — полная глубина) и следствие исходных уравнений
, Пд
и + 7
0,
(9)
д=1
из которого вычисляется значение и |д1, используемое для нахождения нового положения точки уреза с помощью второго краевого условия дифференциальной задачи
^ж
и
и
(10)
д=1
На участке "водозабора" граничное условие определяется заданием расхода жидкости Заметим, что так как течение имеет докритический характер, то в дифференциальной задаче на такой границе ставится только одно условие, а необходимые для замыкания вычислительного алгоритма дополнительные значения определяются численным интегрированием по характеристике, приходящей на границу изнутри области.
При решении двумерных (плановых) задач граничные значения на вертикальной стенке определяются также с использованием условия непротекания и с привлечением соответствующих исходных дифференциальных уравнений для вычисления производной по нормали к границе от уровня свободной поверхности ' [10]. Для реализации граничных условий на наклонной стенке привлекаются двумерные аналоги соотношений (9), (10). При этом линия уреза определяется системой из двух уравнений типа (10), решение которой параметрически описывает искомую кривую.
Основные проблемы расчета граничных значений в задачах о потенциальных течениях жидкости в двух- и трехмерных постановках детально описаны в работах [2, 7].
4. Тестовые и модельные расчеты
Описываемые здесь расчеты имели целью, во-первых, выбор и апробацию вычислительных алгоритмов на данном классе задач и, во-вторых, самое предварительное определение характера течений в изучаемом водоеме. На этой стадии использовались наряду с новыми, реализованными специально для исследуемой задачи, также и алгоритмы, успешно применявшиеся ранее для решения других задач волновой гидродинамики, в том числе и моделирования трансформации волн цунами в открытом океане и в прибрежных акваториях.
4.1. Тестовые расчеты
На первом этапе рассматривалась линейная модель теории мелкой воды и в качестве вычислительного теста использовалась задача о распаде начального косинусообразного возвышения на фоне покоящейся жидкости:
Сравнивались результаты, полученные с помощью трех неявных схем (Мак-Кормака, предиктор-корректор Годунова, Кранка—Николсона [11]) и их явных аналогов (Мак-Кор-мака, Годунова).
Сравнение уровней свободной поверхности и скоростей (рис. 5), рассчитанных по неявным схемам, с соответствующим точным решением показывает предпочтительность схем Годунова и Кранка—Николсона, обеспечивающих приемлемую точность даже при достаточно больших числах Куранта.
В этих же расчетах исследовалась зависимость свойств схем от параметра сглаживания а и параметра ш, отвечающего за величину аппроксимационной вязкости. Определенные здесь оптимальные значения этих параметров аопт = 0.5 и шопт=0.6—0.7 оказались пригодными и для решения рассматриваемых ниже нелинейных задач.
Классическая задача о распаде плотины, обладающая в случае ровного дна известным аналитическим решением [12]
и(х- 0) = 0.
и(х- 0) = 0-
Рис. 5. Сравнение неявных схем на линейной задаче о распаде косинусоидального начального возвышения. Распределения на момент времени Ь = 4 полной глубины Н (верхние кривые) и потока Ни (нижние кривые). Результаты получены по схемам Мак-Кормака (а), Годунова (б), Кранка—Николсона (в). Кривые на графике (г) показывают изменение со временем ошибки по глубине в максимальной норме.
Расчеты выполнялись при следующих значениях параметров вычислительных алгоритмов: Ах = 0.2, ш = 0.6, коэффициент сглаживания после шага предиктор а\ = 0.5 и после шага корректор а2 = 0.2, ж = АЬ/Ах = 4.0 (----Годунов,----Мак-Кормак,......Кранк-Николсон).
решалась в рамках нелинейной модели мелкой воды.
В численной реализации движение по сухому руслу заменялось движением по тонкой пленке (толщина пленки £ = 5к0, где ко — начальная глубина слоя воды за плотиной, а 5 ~ 0.001). Расчеты проводились с помощью явных и неявных вариантов схем предиктор-корректор Годунова и Мак-Кормака.
Рисунок 6 демонстрирует сравнительные характеристики неявного и явного вариантов схемы предиктор-корректор Годунова, которая оказывается несколько предпочтительнее схемы Мак-Кормака.
Надо заметить, что численное моделирование в терминах "скорость — полная глубина" приводит к появлению на переднем фронте нефизического разрыва типа ударной волны при любом выборе 5. Переход к переменным "поток — полная глубина" позволяет избавиться от этого эффекта при условии правильного выбора параметра ш. Однако при счете с большими числами Куранта упомянутый здесь эффект появляется вновь и предшествует проявлению нелинейной неустойчивости схем.
Подводя итоги обсуждения тестовых расчетов, можно заметить, что при воспроизведении такого относительно быстрого процесса, как распад плотины, вполне допустимым является использование явных конечно-разностных схем — безусловно простых в реализации и обеспечивающих необходимую точность за приемлемое время.
Что касается рассмотренных неявных схем, то при их использовании не удается зна-
a
б
-0,5
-0,5
е
г
-10
к
-0,5
-0,5
Рис. 6. Сравнение явного и неявного вариантов схемы Годунова для задачи о распаде плотины. Аналитическое решение (сплошная линия) задается соотношением (11). Распределения на момент времени t = 4 полной глубины h (верхние кривые) и потока hu (нижние кривые). Расчеты выполнялись при следующих значениях параметров вычислительных алгоритмов: Ax = 0.2, ш = 0.6, коэффициент сглаживания после шага предиктор а\ = 0.5 и после шага корректор а2 = 0.2; для явной схемы ж = At/Ax = 0.4 (а), для неявной — ж = 0.4 (б), 1.0 (в), 2.0 (г).
чительно увеличивать шаг по времени с одновременным сохранением требуемой точности описания переходной зоны. Решение этой проблемы требует привлечения более сложных в реализации монотонных неявных схем повышенного порядка аппроксимации. Как показывают результаты тестовых расчетов, простые неявные схемы могут тем не менее использоваться в реальных расчетах для получения предварительных оценок.
4.2. Модельные расчеты
Одномерная задача о распаде плотины в предварительно заполненную часть водоема (5 = 0.5) для близкого к реальному неоднородного распределения глубин вдоль основной оси проектируемого водоема (рис. 7) решалась по явной схеме предиктор-корректор с использованием неравномерной адаптированной к распределению глубин фиксированной сетки.
Этот расчет (рис. 8) позволил получить предварительную картину возможного развития процесса при разрушении разделяющей дамбы. Образовавшаяся в результате распада плотины волна распространяется в сторону ограждающей дамбы и, отразившись от нее, движется в сторону "водозабора".
Сознавая приблизительность полученных одномерных результатов, авторы провели двумерное "плановое" моделирование рассматриваемого явления с помощью явной схемы предиктор-корректор, аппроксимирующей на криволинейной неравномерной сетке нелинейные уравнения мелкой воды (2) без учета сил Кориолиса и диссипативных эффек-
-20
-10
20
10
о
Рис. 7. Изолинии глубин проектируемого водоема и его основная ось (штриховая линия). Жирной линией отмечено положение плотины в момент времени Ь = 0.
тов. Модельность ситуации определялась предположением о вертикальности ограждающей дамбы и заменой "водозабора" участком такой же вертикальной стенки.
Как показывают результаты моделирования (рис. 9), основные процессы, происходящие в водоеме, носят ярко выраженный двумерный характер, определяющийся неравномерностью распределения глубин и сложной геометрией ограждающей дамбы — границы расчетной области. Так, при движении волны от линии прорыва плотины к вершине дамбы образуются две группы волн, распространяющиеся вдоль береговой границы. Затем эти группы волн взаимодействуют у вершины и создают волновой объект значительной амплитуды, который инициирует движение в обратную сторону, к "водозабору" — горловине водоема. С очевидностью здесь проявляется преимущественное распространение волновой энергии от правого берега к левому ("разворот" потока против часовой стрелки), приводящее к всплеску у левой границы зоны сужения.
Попытка проверки этих, также предварительных, оценок была предпринята с помощью трехмерной модели потенциального течения жидкости (4) при тех же начальных условиях и с той же геометрией границ. Результаты этого расчета в целом подтверждают изложенную выше картину течения, дополняя ее немаловажными деталями. Так, анализ векторных полей скорости на поверхности жидкости (рис. 10) показывает, что при дальнейшем распространении волны по горловине бухты происходит концентрация кинетической энергии, приводящая к значительному всплеску непосредственно в зоне "водозабора".
Достаточно нетривиальная картина течения позволяет уже на стадии предварительно-
Рис. 8. Динамика свободной поверхности в одномерной задаче о распаде плотины с учетом реального рельефа дна (теория мелкой воды).
Рис. 9. Динамика свободной поверхности в двумерной задаче о течении, инициированном разрушением временной дамбы. Изолинии свободной поверхности и ее пространственные изображения на различные моменты времени.
г=151
Г=227
Рис. 10. Трехмерное моделирование течения, инициированного разрушением временной дамбы. Векторные поля скорости на свободной поверхности и ее пространственные изображения на различные моменты времени.
Рис. 11. Профили свободной поверхности (I а, II а) и скоростей (I б, II б), рассчитанные на различные моменты времени при насосном режиме работы агрегатов ГАЭС для модели Сен-Венана (I) и одномерной модели мелкой воды (II).
----входящая волна,----заплеск на стену, ......волна у водозабора.
го оценочного моделирования предположить наличие по крайней мере трех критических по заплеску зон — вершины ограждающей дамбы, левой границы зоны сужения и зоны "водозабора".
Следующая модельная задача, рассматриваемая авторами, связана с воспроизведением течений в водохранилище, обусловленных насосным и турбинным режимами работы агрегатов ГАЭС. Одномерные расчеты выполнялись в рамках моделей нелинейной мелкой воды и уравнений Сен-Венана с использованием тех же расчетных сеток, вычислительных алгоритмов (для уравнений Сен-Венана использовалась конечно-разностная схема (7)) и "реального" рельефа дна, что и в задаче о разрушении временной дамбы. Начальный уровень свободной поверхности к0 в насосном режиме соответствовал уровню мертвой отметки, в турбинном — нормальному подпорному уровню, условие на границе расчетной области, соответствующей "водозабору", определялось фактическим расходом (см. рис. 3), отнесенным к единице площади входного сечения.
Результаты расчетов, приведенные на рис. 11,12, получены с использованием равномерных сеток (Ы=400). Здесь сплошной линией показаны графики, соответствующие начальному этапу процесса (волна у зоны расширения водоема), штриховой линией — моменту
заплеска на ограждающую дамбу, пунктиром — моменту подхода волны к водозабору. На врезке схематически изображены распределения глубин вдоль оси водоема и его контур (в метрах).
Как показывают эти результаты, обе используемые модели приводят к близкому качественному описанию процесса, воспроизводя движение скачка уровня по акватории, его взаимодействие с вертикальной стенкой и распространение отраженной волны. Однако квазидвумерность модели Сен-Венана позволяет учесть неоднородность ширины водоема, приводящую к уточнению качественных и количественных характеристик течения.
Из рис. 11,12 видно, что в зонах сужения водоема как уровни свободной поверхности, так и значения скоростей, рассчитанные по модели Сен-Венана, превышают соответствующие значения, определенные по одномерной модели мелкой воды (воспроизводящей движение по акватории постоянной ширины). Заметим, что незначительные осцилляции вычислительного характера за фронтом волны связаны с неоднородностью поля глубин и контура водоема, а также с разрывным типом воспроизводимого течения.
Рис. 12. Профили свободной поверхности (I а, II а) и скоростей (I б, II б), рассчитанные на различные моменты времени при турбинном режиме работы агрегатов ГАЭС для модели Сен-Венана (I) и одномерной модели мелкой воды (II).
----входящая волна,----заплеск на стену, ......волна у водозабора.
Приближая постановки задач к условиям реального моделируемого объекта, авторами были выполнены расчеты по воспроизведению различных режимов эксплуатации водое-
ма ГАЭС в рамках плановой модели мелкой воды. Применяемая модель позволяет учесть пространственное распределение насосных (турбинных) агрегатов на входной границе, что в сочетании с неоднородностью поля глубин в этой части акватории, по существу, определяет гидродинамические характеристики начальных этапов развития процесса.
Вычислительные эксперименты проводились с использованием явной схемы предиктор-корректор на неравномерной криволинейной сетке, адаптирующейся к геометрии водоема.
На рис. 13 показаны рассчитанные поля вектора скорости при насосном режиме работы агрегатов. В отличие от расчетов по модели Сен-Венана здесь реально учитывается геометрия акватории, в том числе искривленность и несимметричность ограждающей дамбы, что проявляется в существенно выраженном двумерном характере воспроизводимого течения, где явно видны зоны зарождения водоворотов, струйных течений и застойные зоны. В то же время качественная картина потока близка к картине, предсказываемой по результатам одномерных расчетов.
Рис. 13. Поля скоростей в задаче о наполнении водоема. Расчет по двумерной модели мелкой воды.
Сравнение результатов, полученных с использованием модели Сен-Венана, и осреднен-ных по площади сечения уровней свободной поверхности, рассчитанных в ходе двумерного моделирования по уравнениям теории мелкой воды, позволило определить степень достоверности квазидвумерной гидравлической модели. Как видно из рис. 14, модель Сен-Венана достаточно адекватно воспроизводит динамику свободной поверхности, учитывая при этом изменение геометрии водоема.
Необходимо заметить, что в целом процесс характеризуется наличием существенного скачка полной глубины, распространяющегося в начальный период по зоне покоя, а после отражения — навстречу потоку.
Обсуждавшиеся выше результаты двумерного и трехмерного моделирования демонстрируют наличие довольно сложных по своей природе гидродинамических процессов, развивающихся в области водоема, непосредственно примыкающей к вершине ограждающей дамбы. Именно с этими процессами связывается нарастающее различие описываемых
здесь моделей, которое приводит, в частности, к значительному фазовому сдвигу, наблюдаемому в поведении отраженных от дамбы волн, распространяющихся в направлении водозабора.
Проведенные модельные расчеты позволяют надеяться на успех применения двумерной плановой модели мелкой воды при решении рассматриваемой задачи, в то же время очевидна необходимость включения в компьютерную модель квазидвумерной модели Сен-Венана и одномерной модели мелкой воды для получения предварительных оценок, а также трехмерных полных гидродинамических моделей для проведения контрольных уточняющих расчетов.
Рис. 14. Уровни свободной поверхности в задаче о наполнении водоема. Расчет по модели Сен-Венана (сплошная линия) и осредненные по площади сечения результаты, полученные в ходе двумерного моделирования (теория мелкой воды). Моменты времени соответствуют начальной стадии процесса — волна у начала зоны расширения водоема (1), приближению к ограждающей дамбе (2) и подходу отраженной волны к зоне водозабора (3).
Обсуждая дальнейшую программу модельных испытаний, следует указать на целесообразность анализа применимости двумерной линейной модели, некоторых двумерных и квазидвумерных нелинейно-дисперсионных моделей мелкой воды и на необходимость исследования вертикальной структуры процесса с помощью трехмерной модели для определения вертикальной изменчивости рассчитываемых полей.
5. Структура программной реализации компьютерной модели
Программная реализация компьютерной модели требует выбора адекватного инструментария и средств создания ее отдельных компонент.
Структура модели видится состоящей прежде всего из управляющей программы, блока пользовательского интерфейса, блока ввода и редактирования начальных данных, блока описания и модификации геометрии расчетной области, блока расчета сетки, собственно набора вычислительных модулей, реализующих комплекс гидродинамических моделей, и блока визуализации промежуточных и окончательных результатов.
Управляющая программа должна предоставить пользователю возможность описания исходной ситуации и сценария процесса, в том числе одного из стандартных, хранящихся в системе, описания геометрии расчетной области, выбора математической модели и алгоритма ее реализации, причем альтернативой любому выбору должен быть один из стандартных, рекомендуемых вариантов.
5.1. Обоснование выбора среды разработчика для создания компьютерной модели
Выбор программной среды разработчика определяют, в частности, следующие факторы: возможность адекватной реализации вычислительных алгоритмов и пользовательского интерфейса,
простота написания программ (технологичность создания текста), возможность их модификации и расширения,
оптимальность генерируемого кода,
универсальность среды, портабельность, (переносимость, адаптируемость на различные типы компьютеров, многоплатформность),
простота в освоении и эксплуатации (с учетом квалификации пользователей). Удовлетворить этим требованиям одновременно практически невозможно в силу существования объективных противоречий между отдельными положениями. Так, естественное желание разработчика на близком ему языке легко определять объекты и отношения между ними может привести к существенным затратам системных ресурсов для реализации алгоритмов в командах и типах данных процессора и внешних устройств. Универсальность среды может препятствовать ее быстрому освоению и эффективному применению, а модификация компьютерной модели может потребовать от разработчика использования новых программных и вычислительных технологий.
Анализ имевшегося в распоряжении авторов программного инструментария привел к выводу о невозможности в разумный срок реализовать проект с мощными вычислительными возможностями и развитым интерфейсом в рамках какой-либо одной из рассмотренных технологий:
Библиотека Быстродействие Простота
интерфейса кода освоения
1. Borland C++ v.4.0 + + +
2. Borland Pascal v.7.0 + + +
3. MS Visual C++ v.2.0 ! + +
4. MS Visual Basic v.3.0 Pro ! - !
5. MS Fortran v.5.1 - ! +
6. Watcom C/C++ v.10.0 ! + +
7. Watcom Fortran 77 v.9.5 - ! + здесь "!" — отлично, "+" — хорошо, "-" — плохо.
С точки зрения авторов, комбинация MS Visual Basic v.3.0 Pro и MS Fortran v.5.1 оказывается хорошим выбором, позволяющим эффективно использовать достоинства каждой из этих программных сред и компенсировать присущие им недостатки. Таким образом, модули интерфейса (пользовательского и программного), от которых не требуется высокого быстродействия, разрабатывались в среде Visual Basic, а вычислительные алгоритмы — в среде MS Fortran, являющейся, по существу, стандартом в этой проблемной области.
Концепция объектно-ориентированного программирования применяется в Visual Basic
исключительно для создания пользовательского интерфейса, а в остальном этот язык остается, так же как и Fortran, процедурно-ориентированным. В качестве операционной среды для функционирования компьютерной модели водохранилища выбрана графическая оболочка MS Windows.
5.2. Принципы программной реализации компьютерной модели
Создание сложной компьютерной модели требует адекватного решения вопросов согласования и интеграции различных модулей в единую программную систему. При этом должны быть выполнены следующие условия:
совместимость сред разработчика на уровне выходных модулей, единый стиль программирования для написания сопрягаемых модулей, обеспечение расширения возможностей программного комплекса.
Что касается совместимости сред, то здесь гарантией стали общие принципы создания используемого программного инструментария и тщательное всестороннее его тестирование. Для согласования программных модулей, написанных на разных языках, использовалась технология динамически подключаемых библиотек (DLL).
Можно сказать, что модули рассматриваемой здесь программной системы представляют собой наборы подпрограмм (функций и процедур), объединенных по типам выполняемых операций.
Взаимодействие блока интерфейса с вычислительными осуществляется через модуль согласования, обеспечивающий связь с выбранным алгоритмом, и через соответствующий модуль адаптации, выполняющий передачу, сохранение и графический вывод данных. Каждый вычислительный алгоритм реализован в виде отдельной библиотеки DLL (например, MODEL1D.DLL — библиотеки расчета по одномерной модели мелкой воды), поэтому модификация какого-либо алгоритма (или добавление новых) не влияет на работу остальных — достаточно изменить (или создать новый) модуль адаптации.
Вычислительный блок выполняет три основные функции: инициализацию расчета, расчет параметров вычислительного алгоритма, расчет одного шага по времени.
Процесс инициализации предполагает расчет сетки и начальных значений искомых функций (полная глубина, поле скоростей). Для каждого из моделируемых режимов (наполнение — опорожнение) предусмотрен соответствующий набор данных, описывающих начальные и граничные условия, которые определяют начальный уровень невозмущенной поверхности, начальное поле скоростей, характеристики потока на входной границе.
Параметры алгоритма (смещение по времени, пересчет сетки и т.д.) могут вычисляться не только при инициализации расчета, но и на каждом временном шаге или через несколько шагов в зависимости от характеристик моделируемых явлений и эффективности текущего расчета.
Для обеспечения гибкости расчета вычислительный блок возвращает управление головной программе после выполнения каждого временного шага. Это обусловлено, в частности, псевдо-многозадачностью MS Windows, которая требует от разработчика предусмотреть явную передачу управления операционной системе. Таким образом, интерфейсная и вычислительная части программной системы существуют не как две взаимодействующие программы, а как модули одной программы.
Интерфейсная часть обсуждаемой программной системы обеспечивает описание вычислительной области, начальных и граничных условий; выбор алгоритма расчета; вывод результатов; управление расчетом.
Для конструирования геометрии водоема в состав компьютерной модели включен так называемый редактор батиметрии, который создает образ акватории.
В целом геометрия водоема описывается набором опорных точек, опорных сечений, и опорных констант (см. рис. 1). Часть опорных точек вместе с заданными в них углами определяет набор поперечных опорных сечений водоема, на которых, в свою очередь, устанавливаются опорные точки с заданными опорными константами — глубинами и опорные точки, определяющие границу шва (дноограждающая дамба), в каждой из таких граничных опорных точек указываются углы наклона нормалей к линии шва по отношению к горизонтальной плоскости. Откладывая на каждой такой нормали заданную опорную константу — высоту ограждающей дамбы, можно полностью завершить описание геометрии водоема ГАЭС. Предложенная структура данных, идентифицирующая структуру водоема, позволяет эффективно организовать процессы ввода, хранения, визуализации и модификации геометрии изучаемого объекта.
Геометрия водоема может создаваться заново или модифицироваться в терминах интерактивного редактирования либо векторного изображения, либо набора числовых значений. В основной программе задание акватории осуществляется указанием файла, в котором (в специальном формате) хранится описывающая ее информация (созданная редактором батиметрии).
Управляя расчетом, пользователь получает возможность не только выбора модели, алгоритма, точности, параметров сетки, но и установки сечений, точек и времен фиксации состояния системы с последующей записью в файл. Кроме того, результаты компьютерного моделирования выводятся на дисплей и принтер непосредственно в процессе расчета. При этом можно задавать режим визуализации: изолинии и пространственные изображения свободной поверхности, векторные поля скоростей (или потоков), графики изменения рассчитываемых величин в точках, эпюры скоростей и т. д. Расчет может быть приостановлен, продолжен и окончательно завершен.
Внешний графический модуль предназначен для высококачественной визуализации результатов расчетов в виде векторных полей, полей изолиний, изображения пространственных картинок и графиков изменения величин (уровня) в заданных точках с возможностью вывода на экран.
6. Заключение
Подводя итоги сказанному, отметим, что процесс создания компьютерной модели находится в начальной стадии, в частности, не рассматривались вопросы моделирования режимов работы гидроузла в зимнем режиме с учетом ледяного покрова. В дальнейшем авторы предполагают расширить круг используемых математических моделей, включив нелинейно-дисперсионные приближения мелкой воды, модели плоских вихревых течений идеальной жидкости, а также предусмотреть возможность применения неявных алгоритмов для исследования характерных волновых режимов водоема. Предстоит определить и оптимальные стратегии построения неравномерных сеток.
Авторы выражают признательность В. Х. Давлетшину за постановку этой задачи и полезные обсуждения отдельных ее аспектов.
Список литературы
[1] ДАВЛЕТШИН В. Х. Волны в водоемах ГАЭС и защита их берегов. Вычислительные технологии, ИВТ СО РАН, Новосибирск, 4, №11, 1995, 104-107.
[2] ВольцингЕр Н.Е., Клеванный К. А., ПЕлиновский Е. Н. Длинноволновая динамика прибрежной зоны. Гидрометеоиздат, Л., 1989.
[3] Новиков В. А., Федотова З. И., Кузьмичева Т. В. О некоторых проблемах моделирования наката длинных волн на берега сложных очертаний. Вычислительные технологии, ИВТ СО РАН, Новосибирск, 2, №4, 1993, 196-209.
[4] ХАкимзянов Г. С. О краевых условиях для конечно-разностной схемы расчета потенциальных течений жидкости со свободной границей. В "Актуальные проблемы современной математики", НИИ МИОО НГУ, Новосибирск, 1, 1995, 164-174.
[5] Барахнин В. Б., ХАкимзянов Г. С. Численная реализация краевых условий в одномерных задачах теории мелкой воды. Там же, 18-30.
[6] Марчук Ан. Г., Чубаров Л. Б., Шокин Ю. И. Численное моделирование волн цунами. Наука, Новосибирск, 1983.
[7] Эббот М. Б. Численная гидравлика. Гидравлика открытого потока. Энергоатомиз-дат, М., 1983.
[8] Грушевский М. С. Неустановившееся движение воды в реках и каналах. Гидроме-теоиздат, Л., 1982.
[9] Шокин Ю. И., ХАкимзянов Г. С. Конечно-разностный метод расчета трехмерных потенциальных течений жидкости со свободной поверхностью. Вычислительные технологии, ИВТ СО РАН, Новосибирск, 1, №1, 1992, 154-176.
[10] Годунов С. К., Забродин А. В., Иванов М. Я., Крайко А. Н., Прокопов Г. П.
Численное решение многомерных задач газовой динамики. Наука, Гл. редакция физ.-мат. лит., М., 1976.
[11] Маккормак Р. В. Численный метод решения уравнений вязких течений. Аэрокосмическая техника, 1, №4, 1983, 114-123.
[12] Barakhnin V. B., Khakimzyanov G. S. Adaptive-grid numerical solution of one-dimensional and two-dimensional problems for the shallow-water equations. In "Advanced Mathematics: Computations and Applications": Proc. of AMCA-95, A.S.Alexeev, N.S. Bakhvalov, eds., Novosibirsk, 1995, 144-153.
[13] Барахнин В. Б., ХАкимзянов Г. С. Численная реализация краевых условий в плановых задачах для нелинейных уравнений мелкой воды. В "Актуальные проблемы современной математики", НИИ МИОО НГУ, Новосибирск, 2, 1996, 3-12.
[14] Флетчер К. Вычислительные методы в динамике жидкостей: В 2-х томах. Мир, М., 1991.
[15] Стокер Дж. Дж. Волны на воде. Математическая теория и приложения. Иностр. лит., М., 1959.
Поступила в редакцию 30 мая 1996 г.
