CC BY
108
Вычислительные технологии
Том Т, № 5, 2002
ОБОСНОВАНИЕ ЧИСЛЕННОГО МЕТОДА ДЛЯ МОДЕЛИРОВАНИЯ НАКАТА ДЛИННЫХ
В данной работе описан подход к численному моделированию взаимодействия длинных волн с берегом. Изучение этого явления в первую очередь представляет прикладной интерес, так как выход длинных волн различного происхождения на сушу часто вызывает катастрофические последствия. Прогресс в развитии вычислительных технологий способствует тому, что при проектировании гидротехнических и гражданских сооружений активно учитываются прогнозы, полученные с помощью численного моделирования. Исследование подобных вопросов интересно также с теоретической точки зрения. Среди разнообразных задач о поведении волн на воде именно задачи о трансформации волн в прибрежной полосе с затоплением береговой зоны или осушением дна являются наиболее сложными ввиду выраженной нелинейности процесса и наличия подвижной границы, на которой вырождаются основные характеристики течения: полная глубина жидкости Н и, как следствие, скорость распространения возмущений с = (дН) 2 обращаются в нуль. Различные постановки задач такого рода и методы их решения рассмотрены в большом количестве публикаций. Укажем лишь посвященную этой проблеме книгу [1], где приведена довольно полная библиография. Обзор работ на эту тему показывает, что наиболее популярной математической моделью, применяемой для изучения динамики длинных волн и позволяющей получить интересные с практической точки зрения результаты, являются нелинейные уравнения мелкой воды.
В настоящее время для моделирования волнового режима в реальных акваториях со свободным выходом волн на побережье в рамках нелинейных уравнений мелкой воды чаще
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №00-01-00899.
© З. И. Федотова, 2002.
З. И. Федотова Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: zf@@ict.nsc.ru
The numerical method of the modelling of a long waves runup on a coast is investigated. The purpose of the work consists in the method validation and estimation of its accuracy using numerical experiments. The model problem of long wave runup on a uniform slope having the analytical solution within the framework of the nonlinear equations of shallow water is considered as test problem.
Введение
других используются методы конечных разностей (см., например, книги [1, 2], содержащие обширный материал о методах расчета взаимодействия волн с берегом).
В данной работе описан конечно-разностный метод сквозного счета, использующий прямоугольную сетку в фиксированной области. Логическая простота алгоритма решения и устойчивость вычислительного процесса делают предложенный метод пригодным для решения широкого класса задач волновой гидродинамики в бассейнах с произвольным рельефом дна и сложными очертаниями береговой линии. Как правило, “универсальность” метода влечет определенную потерю точности. Однако для сходящихся конечноразностных методов существует принципиальная возможность повышения точности путем измельчения шага сетки в наиболее важных областях течения.
Методами сквозного счета, применяемыми для решения задач с нестационарной областью течения (т. е. с внешней подвижной границей), мы будем называть методы, в которых область поиска решения распространяется до прямоугольной области, включающей водную акваторию и прилегающую к берегу часть суши. Этот прямоугольник покрывается сеткой, на которой определена дискретная функция со значениями, равными либо глубине невозмущенной жидкости, если узел сетки принадлежит водной акватории, либо взятому со знаком “минус” возвышению суши над невозмущенным уровнем воды, если узел принадлежит берегу. Для построения численной модели требуется выбрать подходящую интерпретацию среды, соответствующей области, занимаемой сушей. Известно несколько способов решения этого вопроса.
Способ 1. На сушу “натягивается тонкая пленка воды, удерживаемая но откосах силой трения” [3].
Способ 2. В области, занятой сушей, выстраивается искусственное ровное дно, т. е. водная акватория продолжается каналом, с тем чтобы проводить вычисления в области простой геометрии.
Способ 3. В точках суши скорость и полная глубина жидкости полагаются нулевыми.
Что касается способа 1, то в литературе он описан очень кратко (авторы ограничиваются лишь формулировкой идеи на физическом уровне). Исследование вопроса об адекватности используемой математической модели в публикациях отсутствует. Понятно, что в обосновании этого способа главным является вопрос о математическом описании механизма трения, который не допускал бы осушение берегового склона и при этом гарантировал выполнение закона сохранения массы и энергии в пределах допустимой ошибки.
Способ 2 предложен В. А. Новиковым, а также в работе [4], и рассмотрен нами в одномерном случае (см. [5]).
В данной статье речь пойдет о способе 3, так как он имеет простую физическую интерпретацию, а сформулированную на его основе математическую модель удалось обосновать и использовать для построения численного метода, пригодного для решения задач о динамике поверхностных волн, включая движение волны по сухому берегу, в акваториях с произвольным рельефом дна и сложными очертаниями береговой линии. Этот метод был успешно использован при моделировании ряда задач, воспроизводящих лабораторные эксперименты, а также применен для моделирования реального цунами, случившегося в 1993 году вблизи японского острова Окушири [6, 7].
1. Математическая модель
Целью данного параграфа является построение математической модели, на основе которой в следующем параграфе будет построен численный алгоритм для описания процесса
выхода длинных волн на берег. Здесь, как и в подавляющем количестве работ, посвященных решению подобных задач (см., например, библиографию из книги [1]), в качестве основных дифференциальных уравнений взяты классические нелинейные уравнения мелкой воды
дШ дГ 80 ^
+ дХ + 8У = ф (1)
где
W
F
дt
/
і
V
Hu
\
Hu2 + 2 H2 Huv
G
/
Hv
Huv
Hv2 + g-H2 2
Ф = | дИкх — ги(и2 + V2)1/2 кдИку — то(и2 + v2)1/2y
х,у — точки на горизонтальной плоскости ХОУ; Ь — время; и = и(х,у,Ь)^ = v(x,y,t) — компоненты усредненной по вертикали скорости жидкости в направлениях х, у соответственно; И = И(х,у,Ь) — полная глубина жидкости; к = к(х,у) — глубина невозмущенной жидкости; г — коэффициент трения; д — ускорение свободного падения.
Проинтегрируем уравнения (1) по произвольному объему ш(Ь) из декартова пространства (х,у,Ь), в котором определены гидродинамические переменные, входящие в уравнения (1):
Ш (£+%+£Н>=///«*). (2)
ш(Ь) ш(Ь)
Применяя формулу Остроградского к левой части интеграла по объему, получим
/д^№ дF д^\ I I
+ д—+ ~8~)du(t) = 11 [W cos(n,t) + F cos(n,x) + G cos(n, y)] da (и),
(t)
а(ш)
где а(ш) — поверхность, ограничивающая объем ^(і); п — внешняя нормаль к поверхности а(ш). Переходя к элементу площади в прямоугольной декартовой системе координат, полученный интеграл перепишем в виде
JJ [Ш сов(п,і) + Г сов(п, х) + 0 соs(n,y)^dа(u) = JJ Wdxdy + Fdtdy + Gdtdx,
а(со)
где интегрирование производится по внешней стороне поверхности. Наконец, переходя в правой части интегрального уравнения (2) к элементу объема dxdydt в прямоугольных координатах, получаем уравнение
W dxdy + Fdtdy + Gdtdx =
Фdxdydt.
(З)
а(ш)
ш(Ґ)
о
Учитывая специфику переменной t, объем u(t) можно рассматривать как прямое произведение u(t) = и* х [ti,t2], где и* — сечение u(t) гиперплоскостью t = const, t1,t2 £ [0,T]. Обоснование проведенных преобразований можно найти в книге [8] (см. гл. 3, §§ 6, 9).
Отметим, что в работе [9] уравнения в интегральной форме (3) были выведены из общих законов сохранения импульса и массы несжимаемой жидкости при предположении о гидростатичности давления. Применение к (3) формулы Гаусса — Остроградского приводит после ряда известных преобразований к системе дифференциальных уравнений (1). Соотношения (3) называют интегральными законами сохранения для уравнений мелкой воды.
Условия эквивалентности интегральной и дифференциальных формулировок известны: если векторная функция Ш(х,у,Ь) £ С1 удовлетворяет интегральным законам сохранения (3) для любых замкнутых поверхностей о(ш) и ограничиваемых ими объемов ш(Ь), а Г £ С1, О £ С1, Ф £ Со, то функция Ш(х,у,Ь) удовлетворяет системе дифференциальных уравнений (1), и наоборот, всякое решение Ш(х,у,Ь) системы (1) (Ш(х,у,Ь) £ С1) удовлетворяет интегральным законам сохранения (3). В теории дифференциальных уравнений гиперболического типа интегральные соотношения вида (3) применяют для введения обобщенного решения системы (1) [10]. В области численных методов интегральный вид уравнений используется для построения консервативных разностных схем. В данной работе интегральные соотношения (3) использованы для формулировки математической модели наката длинных волн на берег, конечно-разностная аппроксимация которой привела к экономичному численному алгоритму. Построение модели основано на введении некоторой расширенной области и доопределении на ней всех функций, входящих в уравнения мелкой воды.
Рассмотрим прямой цилиндр П* = П х [0, Ь], где П — произвольная фиксированная область в плоскости переменных х, у, физически соответствующая проекции водной акватории с прилегающим участком берега и островными системами на горизонтальную плоскость. Область П можно рассматривать как объединение подобластей, одна из которых (вообще говоря, многосвязная) соответствует “воде”, а остальные — “суше”. Для
математического задания подобластей используем функцию И(х,у,Ь) = к(х, у) + п(х, у, Ь),
описывающую полную глубину жидкости и определенную в области, занятой водой (обозначим эту область символом П(1)). Согласно физическому смыслу функция И(х,у,Ь) неотрицательна, причем И(х,у,Ь) = 0 определяет линию уреза, т. е. линию, вдоль которой граничат водная акватория и берег (правильнее говорить о совокупности линий уреза). Распространим функцию И(х,у,Ь) на всю область П, доопределив ее в точках “суши” П(2) = П \ П(1) нулем:
И/ \ И(х,у,Ь), (х,Ь) £ П(1), . ,
И(х, ул) = | о, <хи £ П<2>. (4)
По построению область П не зависит от времени, но является объединением зависящих от времени областей, граничащих в каждый момент времени вдоль линии уреза, которую можно трактовать как нестационарную внутреннюю границу области П. На рис. 1 показан простейший случай, когда каждая из областей П(г), г = 1, 2, является односвязной.
Таким образом, П = П(1)(Ь) и П(2)(Ь), где П(1)(Ь) = (х, у : И(х,у,Ь) > 0), П(2)(Ь) = (х, у : И(х, у, Ь) = 0). Для функций и, V также построим расширения и, V, определив их в области П(2)(Ь) нулями. По поводу вида и свойств коэффициента трения в научной литературе ведется дискуссия. В ряде работ [4, 11] приводятся соображения, что вблизи линии уреза он ограничен. Поэтому есть основания приписать членам трения при нулевой скорости также нулевые значения.
Наконец, во всех точках области, включающей водную акваторию и прилегающие к ней берега, определим функцию, характеризующую рельеф дна и суши. Пусть П = П(1) У П(2),
Рис. 1. Область моделирования — прямоугольник О = О(1) (£) У О(2)(£): 1(0) — начальное положение линии уреза при £ = 0, 1(£) —положение линии уреза при £ > 0, О(1)(£) — область, занятая
водой, О(2) (£) — суша.
где П(1) — область, занятая находящейся в покое жидкостью. Тогда на П(1) определена неотрицательная функция к = к(х, у), равная глубине покоящейся жидкости (“подводная” точка). Зададим на области П(2) неположительную функцию к = к(х,у) со значениями, равными взятому с противоположным знаком возвышению суши над спокойным уровнем воды. Кривая к = к(х,у) = 0 соответствует береговой линии. Без ограничения общности можно считать, что кривые к = к(х,у) = 0 и к = к(х,у) = 0 полностью совпадают. Таким образом, на всем прямоугольнике П можно рассмотреть функцию
И( ч ( к(х,у), (х,ь) £ п(1), , .
^Нмх.у), Ш) £ а*, (5)
причем И = к(х,у) = 0 определяет линию уреза невозмущенной воды.
Как правило, функцию, задающую рельеф всей области течения, включая зону наката, обозначают той же буквой, что и глубину невозмущенной жидкости: к = к(х,у) для всех (х,у) £ П(1)(Ь) ( что нами уже было использовано при написании уравнений (1) и (3)).
Пусть в некоторый момент времени (для определенности Ь = 0) в окрестности берега вода находится в невозмущенном состоянии. Тогда П(1)(Ь)|*=0 = П(1), П(2)(Ь)|*=0 = П(2).
Итак, в области П определены все функции, стоящие под знаком интеграла (3).
Обозначим
IV = \Нй\ , Р = Г(IV), (0 = О(Ш), Ф = Ф(Ш, к).
\И V
Выпишем следующее интегральное соотношение, аналогичное (3):
IV dxdy + Г dtdy + Gdtdx = Ф dxdydt, (6)
<т(П) О
где 0 = 0 х [іі, і2], 0 = 0(1)(і) У 0(2)(і), 0(1)(і) = шь.
Справедливо следующее утверждение.
Утверждение. Если векторная функция Ш(х,у,і) удовлетворяет интегральным законам сохранения (3), то векторная функция Ш(х,у,і) удовлетворяет (6), где 0(1)(і) = и наборот.
Доказательство. Доказательство следует из вышеизложенного. В самом деле, из определения функций с “крышками” вытекает, что всюду в области П(2) (t) вплоть до границы вектор-функции W, F,G, Ф обращаются в нуль и, следовательно, интегралы вида
(6) по объему H(2)(t) х [t1,t2] и его поверхности обращаются в нуль. Тогда из свойства аддитивности (по областям интегрирования) кратных интегралов вытекает доказательство прямого и обратного утверждений.
В дальнейшем символы, обозначающие функции, будем писать без “крышек”.
Заметим, что функции H = 0,u = 0,v = 0 удовлетворяют уравнениям мелкой воды, записанным в дивергентной форме. Однако эта тройка функций не удовлетворяет недивергентным уравнениям
Ht + (Hu)x + (Hv)y = 0, ut + uux + vuy + gHX ghX^ Vt + uvx + vvy + gHy = ghy,
если только h = const (здесь трение не учитываем, так как в противном случае из-за деления на H нужны дополнительные рассуждения). Другими словами, вводя расширение функций, мы переходим к определению некоторого слабого решения, связанного с интегральной формулировкой законов сохранения. В терминах дифференциальных уравнений решение задачи о накате с применением расширенной области (вода — суша) теперь можно сформулировать как склеивание на внутренней границе областей H(1)(t) и tt(2)(t) решения системы дивергентных уравнений, описывающих распространение волны по воде, с нулевым решением, при этом граница области сама подлежит определению. В конечном результате нас интересует только та часть слабого решения, которая определена в области Q(1)(t), включая ее границу.
Отметим, что можно рассмотреть модель наката, использующую недивергентную форму уравнений мелкой воды. Для этого определим функцию h = h(x,y) следующим образом: пусть в области H(1)(t) функция h = h(x,y) определена так же, как и в предыдущем случае (как функция, описывающая реальный рельеф), а в области H(2)(t) — h = h(x,y) = const = 0. В этом случае функция, описывающая рельеф поверхности (т. е. рельеф дна и затапливаемой суши), становится зависящей от времени. Значения полной глубины и скорость в точках суши (h = h(x, y) = 0) по-прежнему будем полагать равными нулю. В рамках такой интерпретации основных функций можно сформулировать аналогичную описанной выше задачу о склеивании решений недивергентной системы уравнений с нулевым решением.
Основной результат вышеизложенных построений можно кратко сформулировать следующим образом: определив одним из указанных способов функции u,v, H на H(2)(t), мы переходим от задачи с нестационарной областью H(1)(t) к задаче с фиксированной областью П. Для ее решения можно применять однородные конечно-разностные схемы, при этом линия уреза в такой постановке становится внутренней границей раздела двух сред (воды и суши). Для расчета обрушающихся волн более приемлемой является интегральная формулировка (или, что то же, уравнения в дивергентном виде). Наиболее простым вариантом, сохраняющим общность постановки задачи, является выбор области П в виде прямоугольника.
2. Численный алгоритм
Для построения численного алгоритма выпишем аппроксимацию интеграла (3) по параллелепипеду &тк х [пАї, (п + 1)Ді], где &тк — ячейка прямоугольной сетки размером
-^(АХт + АХт+і) х (Аук + Аук+і), АХт = Хт - Хт-1, Аук = ук - Ук-1, Аї — шаг по времени:
Хт+1 Ук+2 t+АtУk+1
(7)
J у 1 - )лхё.у + у у (рт+1Л - рт-1Л^^у+
х 1 У 1 Ь У. 1
'т — 1 к — 2 к — 2
Л X . 1
Ь+АЬ т+1
+ / / 1 - ^к-1 )ЛЫх _ Л! фткАЫЫу.
Ь Хгп — 1 ^тк
Звездочками отмечены усредненные по ячейкам величины. Выбор этих величин позволяет рассмотреть широкий класс разностных схем. В данной работе взяты следующие аппроксимации:
р * _ 1(Рга+1 + Рп ) П * _ 1(/уп+1 + пп )
Р т+1 к 2 т+1к + Р тк ), П тк+1 2 тк+1 + П тк)>
С+1 _ р1), СП+1 _ а(\¥п+1).
Для правой части Ф нужно использовать согласованную аппроксимацию, такую, чтобы тривиальное решение дифференциальных уравнений, описывающее состояние покоя жидкости, являлось также точным решением разностной схемы. Величины, помеченные тильдой, вычислены на шаге предиктор при аппроксимации интегралов по ячейкам Птк размером Ахт х Аук согласно формуле
Хт Ук Ь+АЬ Ук
/ / №П+1 - Я'Пк)йхЧу +1 I (РПк - К-1к)Мс1у+
Хт — 1 Ук — 1 Ь Ук — 1 /о\
(8)
Ь+АЬ Хт
+ I I (Птк — ^-1^ _ !!! фПтк&Ых<1у.
Ь хт—1 ^тк
На равномерной сетке (Ахт _ Ах, Аук _ Ау для всех к,т) из выписанных интегральных соотношений нетрудно получить известную схему Мак-Кормака:
(1-КП.+1 — ЖПк)АхАу + (РПк — РП-1к)А*Ау + (ст,к — Ст-'АА _ Ф"ткАЬАхАу,
Ю1 - "Ск )АХАУ + (С+1 к - К-1 к )АА + (С'тк+2 - С'тк-і )АА = Ф 'тк АААУ.
^ ^
т+ 2 к 2
77 * __ 1 (рга+1 I Т?п \ (-1 * _ 1 (Р^и+1 | /го А
^т+1 к 2 V т+1к + ^тк] , ^ тк+2 2 V тк+1 + ^тк ]
(фШ1 + ФЦ*) , (9)
где
Фптк
\
днпт _
п
т_1/2к
днпк_
п
тк_1/2
Фп+1
тк
в«т+1
т+1/2к
'П+1
тк+1/2
дЙП+1
к тк — кт _1к
Ах
к тк — к тк_ 1
Ау
к т+1к к тк
Ах .
к т к+1 к т к
— ГП 11П (иП
' т 1 /9£>п^ 1 /О£> \ ( и
т_ 1/2к ит_ 1/2к \ ( ит-
{(ит_1/2к)2 + {ут_1/2к^
— ГПтк_1/2<к_1/2 ((итк_1/2)2 + (^к_1/2)2) ^
— ~
п+1
0
((иПт++\/2кУ + (*Пт++\/2к
т+1/2к т+1/2к
Ау
— гтк+1,2Ут\:/2(ттг,^)2 + (С++1/2
~п+1 ~ п+1
тк+1/2 ~ тк+1/2
п+1 )
тк+1/2)
)2 1/2 )2)1/2
/
Здесь индексы т ± 1/2, к ± 1/2 означают, что помеченные ими величины представляют собой полусуммы соответствующих величин из соседних узлов. В расчетах обычно используют более компактный вид схемы Мак-Кормака:
а-С+1 — ^к )АхАу + (гпк — К, _1к ААу + (стк — Ст^ААх = Фттк ААхАу (»;п.+1 — 2 (ЙС+1 + ИС.)) АхАу + 1 ((Гп+1к — гп+ 1)АЬАу + (Сп,++1 — С'Ш^АЬАх
Ф’^АгАхАу.
(10)
Основные свойства этой разностной схемы хорошо известны. Она имеет второй порядок аппроксимации, условием устойчивости в одномерном случае является обычное условие Куранта. Изучению устойчивости схемы в двумерном случае посвящено много работ, что связано в основном с неперестановочностью матриц А = дьГ/дШ и В = д,С/дШ. Подробный обзор по исследованию устойчивости этой схемы дан в [12, 13]. К свойствам разностной схемы (10) мы вернемся в следующем параграфе, а здесь формально опишем алгоритм, применяемый для численного решения задачи о накате.
Подчеркнем, что с помощью разностной схемы (10) определяются значения полной глубины Н'^к и потоков (Ни)птк, (Ну)птк, тогда как компоненты скорости вычисляются по формулам (Ни)тк/Нпк, (Ну)тк/Нпк. Замечено, что при Н ^ 0 профили восстановленных скоростей могут оказаться немонотонными при монотонно изменяющихся Н'^к, (Ни)тк, (Н^)тк. Более того, при Н ^ 0 возможно появление высокочастотных осцилляций на фронте волны, что не только искажает физический портрет решения, но может приводить к развитию неустойчивости расчета (даже если условия устойчивости, полученные на основе спектрального анализа, формально выполнены). Избежать неустойчивости помогает регуляризация величин в окрестности линии уреза, суть которой в том, что при Н'^к < Нт;п (Нт;п — некоторая малая величина, имеющая, например, порядок размера донной шероховатости) выполняется присваивание нулевых значений Н^ = 0,и'тк = 0,у'пк = 0, т. е. точка с номером узла (т, к) объявляется принадлежащей суше. Допускаемую при этом погрешность в тестовых расчетах можно исследовать сравнением результатов для убывающей последовательности Нт;п: Нт;п = 10_2,10_3,10_4,... (в случае обезразмеренной задачи).
Таким образом, алгоритм решения задачи состоит из следующих шагов:
1. Во всех внутренних узлах сетки, покрывающей область П, по разностной схеме находятся НЩ+ \ (Ни)п+1, (НгОпк1.
0
2. Проверяется условие НП+1 > Нтіп; если оно выполнено, то и';+1 = (Ни)"*1 /Щ£1, = (НгГА1 /Н'ПХ1, иначе НЩ1 = О^и^1 = 0,<*+1 = 0.
3. Линия уреза определяется как граница области, в узлах которой Н^Х1 = 0.
4. Если необходимо, то для монотонизации решения применяется процедура сглаживания, не понижающая общего порядка аппроксимации.
Описанный алгоритм ранее изучался на ряде тестовых задач [5-7]. Однако результаты, приведенные в указанных работах, носят в основном качественный характер. Ниже будет рассмотрено аналитическое решение уравнений мелкой воды, сравнение с которым позволило более детально охарактеризовать построенный численный метод.
3. Аналитическое решение одной задачи о накате
В данном параграфе приведено аналитическое решение задачи о накате длинной волны на однородный бесконечный откос, впервые описанное в работе [14]. Это решение нелинейных уравнений мелкой воды часто используется для изучения поведения волн вблизи берега.
Пусть система координат выбрана так, что ее начало совпадает с точкой уреза невозмущенной водной поверхности, ось х направлена в сторону берега, вдоль вертикальной координаты измеряется высота свободной поверхности п — превышение уровня жидкости над ее невозмущенным состоянием. С определенной долей условности схема задачи показана на рис. 2.
В горизонтальной плоскости течение предполагается одномерным. Рассмотрим безразмерные переменные (размерные величины помечены тильдой):
и = и/и0, п = п/аіо, х = Х/10, і = Ї/Т.
Здесь п — амплитуда, и — усредненная по глубине скорость, Т = 10/ад, и0 = (д10а)1/2, Н = (к + п)/даіо, к = -аХ соответствует глубине невозмущенной жидкости при Х < 0, а = tgв, в — угол наклона откоса, 10 — характерная длина. Одномерный вариант уравнений мелкой воды (1) при г = 0 в этих переменных можно переписать в следующем виде:
п + [и(п - х)]х = 0, иг + иих + пх = 0. (11)
т]к
Рис. 2. Схематическое изображение постановки задачи: 1 — фрагмент начального возмущения водной поверхности; 2 — фрагмент заплеска волны на откос. Точка подвижного уреза х8 движется по оси Ох, Я — высота заплеска волны на откос.
Замена переменных
ф
и = —, х = фх/4 — а2/16 — и2/2, Ь = Х/2 — и, п = фх/4 — и2/2 (12)
а
позволяет записать уравнения (11) в виде линейного уравнения второго порядка относительно функции ф:
(афа )а — афхх = 0. (13)
Переход к переменным (12) приводит по крайней мере к двум достижениям: нелинейные уравнения становятся линейными, и задача о накате с подвижной границей становится задачей с фиксированной границей а = 0.
Частное решение уравнения (13) дается формулой ф = А.]0(а) сов(Х), где .]0 — функция Бесселя нулевого порядка, А — константа. Самый важный для применения полученного решения вопрос состоит в том, как дополнить уравнение (13) начальными или граничными условиями, что весьма нетривиально ввиду нелинейности преобразования. В работе [14] найдено общее решение задачи с начальными данными
п(х, 0) = По(х), и(х, 0) = и0 = 0.
В этом случае в плоскости переменных а, Х удается сформулировать следующую задачу с начальными данными:
а(и\\ — иаа) — 3иа = 0, и(а, 0) = 0, и\(а, 0) = 1/2 + 4а-1ха. (14)
Аналитическое решение задачи (14) имеет вид
и = а-13\(га) зт(тХ)1 (т)йт, I(т)= ^^3\(та0)1 (а0)йа0, (15)
оо
где ,11 — функция Бесселя первого порядка. В качестве примера в [14] при Ь = 0 было выбрано следующее однопараметрическое семейство начальных волн:
П = £ а2
х = — 16 + £
2 (а2 + а2)3/2 2 (а2 + а2)5/2
5 а3 3 а5
1 — ^^----------+ -■
2 (а2 + а2)3/2 2 (а2 + а2)5/2
(16)
(17)
где а = 1.5(1 + 0.9£)1/2. Эти волны имеют максимальную высоту, равную £, и высоту 0.9£ при х = —1. Угол подхода волн к линии уреза равен нулю. Начальный профиль волны показан на рис. 3 (горизонтальная линия — невозмущенный уровень воды).
Соответствующая физическая задача такова, что уровень воды к берегу понижается, течение в начальный момент времени удерживается неподвижным, а затем развивается движение: вначале волна взбирается по сухому откосу, а затем начинается откат до полного выравнивания поверхности воды.
Из формул (16), (17) нетрудно вычислить ха и, следовательно, f (а), а затем после подстановки f (а) в формулы (15) и выполнения нужных выкладок находится аналитическое
решение
и =
8е?
а
3
1 — гЛ
х = ———
{(1 — гЛ)2 + а2}3/2 4 {(1 — а)2 + а2}5/2
5/4 — гЛ 3
а2 а2 16
+
1 — 2
(1 — гЛ)2
{(1 — гЛ)2 + а2 }3/2
2 {(1 — гЛ)2 + а
;}б/2
t = -аЛ — и, п = х + — а2а2. 2 16
(18)
(19)
(20)
Здесь г — мнимая единица, ^[ ], ^[ ] обозначают взятие мнимой и действительной частей соответственно от выражений, стоящих в скобках.
На рис. 4 показано изменение профиля свободной поверхности на разные моменты времени, вычисленное в соответствии с вышеприведенными формулами. Для установления зависимости п(х,£) согласно формулам (16)-(20) был применен метод Ньютона вначале для одного нелинейного уравнения (исключение переменной а в формулах (16) - (17) в начальный момент времени £ = 0), а затем для комплекснозначной системы уравнений (18), (20), чтобы исключить а и Л.
Полезными для понимания процесса наката волны являются формулы, описывающие
1
2
Рис. 4. Динамика профиля свободной поверхности на моменты времени £п = £п-1 + 0.5, п = 1,..., 8;
£о = 0.
Рис. 5. Аналитическое решение: скорость и траектория движения линии уреза.
динамику линии уреза а = 0:
и
х
8е 5 Л3 - Л5
"а (1 + л2)4 ’
2 и
1
+ £ -
(1 + Л2 )3
(1 + 3Л2 - 2Л4),
2
аЛ — и.
(21)
(22)
(23)
Поведение скорости и траектории точки уреза ^х/^£ = и3 показаны на рис. 5. Решение ведет себя следующим образом. Вначале скорость движения уреза растет, вода поднимается и “подтапливает” сухой откос, при этом максимальная величина заплеска достигается при Л2 = 5 (т. е. когда скорость движения уреза ^х/^£ = и5, вычисляемая по формуле (21), обращается в нуль). Затем развивается слабое течение в противоположном направлении и начинается выравнивание поверхности жидкости, а скорость уреза асимптотически стремится к нулю.
Отметим, что процесс наката имеет чисто нелинейную природу. Постановка аналогичной задачи в рамках линейных уравнений мелкой воды приводит к тому, что траектория точки уреза задается формулой ^х/^£ = 0, т. е. точка уреза всегда неподвижна.
Приведенное аналитическое решение использовано для тестирования предложенного численного метода.
1
£
4. Тестирование численного алгоритма
Напомним, что описанный выше метод расчета наката длинных волн на берег предназначен для описания течений в сложных областях, где применение более точных методов, использующих априорную информацию о решении или геометрические преобразования области расчета, затруднительно. Очевидно, что при решении одномерных задач с однородным откосом он может уступать по точности методам, использующим криволинейные подвижные координаты.
Для изучения свойств разработанного метода сквозного счета были рассмотрены различные подходы. Прежде всего, для аппроксимации дифференциальных уравнений был применен целый ряд разностных схем с целью выбрать наиболее подходящую. Во-вторых, для оценки эффективности предложенного метода, использующего стационарную сетку, был рассмотрен конкурирующий метод, основанный на отображении нестационарной
области течения в фиксированный прямоугольник (метод подвижных сеток). При этом в обоих методах для аппроксимации были взяты идентичные разностные схемы. Наконец, для получения количественной оценки сходимости метода использовано описанное выше аналитическое решение нелинейных уравнений.
Выбор разностной схемы основан на следующих соображениях. При изучении наката волны на берег одной из главных практических проблем является определение границ затопления побережья. В математической постановке эта цель формулируется как вычисление величины максимального заплеска. Другими словами, требуется как можно точнее определить “крайнее” положение фронта движущейся по берегу волны. В терминах рассматриваемых уравнений это сводится к вычислению границы области, в каждой точке которой хотя бы для некоторого момента времени £ > 0 полная глубина жидкости Н была строго положительной. Известно [10], что при использовании разностных методов скорости движения “разностного” и “истинного” фронтов волны различаются. Скорость “разностного” фронта волны в явной разностной схеме определяется условием Куранта, так как за время Д£ фронт волны “захватывает” соседний узел расчетной сетки, а “истинный” фронт волны отстает (что также следует из условия Куранта). Тем самым линия фронта в разностном счете заменяется некоторой “зоной размытости (размазывания)”, внутри которой содержится истинный фронт. Ввиду необходимости как можно точнее локализовать фронт волны мы не применяли схемы первого порядка точности, а также неявные методы, остановившись на схемах второго порядка аппроксимации [5].
Наиболее подробно был изучен алгоритм, базирующийся на схеме Мак-Кормака. Кроме того, для сравнения мы использовали две схемы на разнесенных сетках. Ввиду немонотонности рассмотренных схем для устранения высокочастотных осцилляций были построены специальные процедуры выборочного сглаживания, теоретически не понижающие порядок аппроксимации. Наилучшие результаты показала схема со сглаживанием, построенная на основе метода Мак-Кормака, как в случае стационарной сетки (описанный выше алгоритм сквозного счета), так и на подвижной сетке. На рис. 4 показано изменение профиля свободной поверхности на разные моменты времени. При выбранном шаге по пространству численные решения, полученные на подвижной и стационарной сетках, и аналитическое решение в масштабе рисунка не различимы. На рис. 6 приведен фрагмент области в увеличенном масштабе, и можно увидеть разницу в поведении профилей, более заметную при расчете на грубой сетке. Здесь уместно проанализировать стратегию выбора шага по пространству Дх. С точки зрения теории сходящихся дискретных методов ясно, что чем мельче шаг, тем точнее расчет. На практике вычислительные ресурсы всегда ограничены, поэтому возникает вопрос об эффективности методов. В волновой гидродинамике для оценки конкурентноспособности численного метода часто используют такую характеристику, как количество расчетных точек, приходящихся на длину волны, необходимое для достижения заданной точности расчета (см., например, [5]). В рассматриваемой задаче ввиду ее специфики определение длины волны затруднительно. Однако видно, что основное изменение профиля волны (см. рис. 3) происходит на отрезке [—1.0, 0.0], при этом глубина проникновения на берег равна 0.2-0.3 единицам при выбранной амплитуде £ = 0.2. Показанный на рис. 4 расчет выполнен при Дх = 0.01, т. е. на расчетную область [—5.0, 1.0] приходится 601 точка. При этом на “длину сухого откоса” укладывается всего 20-25 точек (точнее, 23 точки, так как хтах = Я/а, а =1 (согласно обезразмериванию), Я = 1.157£ (из аналитического решения), £ = 0.2 (см. также рис. 2).
При таком выборе шага Я = 0.231250 в случае расчета предложенным методом, Я = 0.230979 при использовании подвижной сетки и Я = 0.2314 — аналитическое значение
(относительная ошибка составляет меньше 6 х 10-4). Видно, что на крупной сетке (5-6 точек на длину пробега по сухому берегу) численное решение, вычисленное по предложенному методу, вблизи берега довольно сильно искажается, тем не менее величина за-плеска Я = 0.226836 вычисляется достаточно хорошо (относительная ошибка около 2%). Справедливости ради надо заметить, что в практических задачах при расчетах реальных побережий оцифровка, как правило, обеспечивает только грубую сетку [2]. На рис. 7 показаны величины заплеска, полученные на сетках с шагами Ах = 0.04, 0.02, 0.01, 0.005, что соответствует 150, 300, 600, 1200 узлам на всю область.
На рис. 8 для сравнения показаны траектории точки уреза, вычисленные по схеме Мак-Кормака на стационарной и подвижной сетках, в сравнении с аналитической траекторией. На рис. 8, а показаны траектории, полученные описанным выше методом на последовательности вложенных сеток, п — число точек, приходящихся на единицу длины области (Ж = 6п). Рисунок иллюстрирует особенность данного метода, состоящую в том, что вычисленная траектория представляет собой кусочно-ломаную линию, которая монотонно стремится к аналитической траектории при уменьшении шага по пространству (рис. 8, б). Заметим, что “ступенчатость” не связана с выбором Нт;п. В данном случае Нтіп = 2.0 х 10-5 для всех случаев. Очевидно, трактовка скорости точки уреза как про-
Рис. 6. Сравнение решений, вычисленных разными методами, с использованием двух сеток: а — при Ах = 0.04, б — при 0.02.
Рис. 7. Сравнение величин заплеска, полученных по схеме Мак-Кормака на подвижной и стационарной сетках, с аналитическим значением.
--------------- численное
0.31 -------♦------- аналитическое 0.3
Рис. 8. Сравнение траекторий точки уреза, вычисленных разными способами.
изводной от разрывной функции несостоятельна и требует особого подхода. Это плата за чрезвычайную простоту и однородность численного алгоритма. Но это обстоятельство не ограничивает применимость метода для моделирования реальных событий. Как правило, в практических расчетах требуется определить лишь максимальную область затопления. Если все-таки нужно оценить скорость переднего фронта потока воды при его движении по сухому берегу, то это можно сделать экстраполяцией значений скорости изнутри области.
Наконец отметим (рис. 8, в), что применение подвижных сеток даже при больших шагах по пространству хорошо описывает движение линии уреза (здесь преобразование нестационарной области течения отображается на фиксированный отрезок на каждом временном шаге).
Полученные результаты позволяют сделать вывод о том, что точность расчетов при использовании как стационарных, так и нестационарных сеток увеличивается при измельчении сетки. На крупных сетках (в данном случае при Дх = 0.08, 0.04) метод с подвижной сеткой дает в целом более точные результаты. Тем не менее, что касается локальной характеристики — величины заплеска, то разработанный метод при решении данной задачи показывает более высокую точность (см. рис. 7).
Для изучения точности предложенного метода в количественном выражении были выполнены расчеты на последовательности измельчающихся сеток. Длина расчетной области равнялась шести единицам, а число узлов N варьировалось: N = 75• 2т + 1, где т = 1,..., 6. Порядок точности р определялся по формуле Рунге
||^ - И1
Р =
1^2М - ^||
где || • || — норма в пространстве 12; ^, ^2ж — численные решения на сетках с числом узлов N и 2N соответственно; ^ — аналитическое решение. Здесь в качестве функции
^ мы рассмотрели полную энергию Е = — Н(и2 + с2), с2 = дН, а норму вычисляли в
2
пространстве Ь2 согласно формуле
1|Е|| = {/ 1Н(и2 + с2)}'72 = { | 2Н(и2 + с2)}1/2,
П П(1)
что приводит к следующей норме в дискретном пространстве 12:
1 N 1/2 ^ 1 N 1/2
||Е 1 = 1 X! 2Н*(и2 + с2м Ах = 1 ^ 2Н*(и2 + с2м Ах
^ 1 ^ Н (г)>0 ^
В табл. 1 приведена погрешность разностного решения ||Е^ — Е || при измельчении сетки, а в табл. 2 дан порядок точности метода.
Таблица 1
Ошибка по норме 12 разностного решения на последовательных сетках
і 75 150 300 600 1200 2400
0.5 0.1026757 0.0247440 0.0062114 0.0015559 0.0003893 0.0000973
1.0 0.1025273 0.0247435 0.0062112 0.0015560 0.0003893 0.0000974
1.5 0.1021952 0.0247413 0.0062107 0.0015557 0.0003893 0.0000974
2.0 0.1015975 0.0247389 0.0062099 0.0015557 0.0003893 0.0000974
2.5 0.1007275 0.0247367 0.0062089 0.0015554 0.0003893 0.0000974
3.0 0.0996548 0.0247343 0.0062074 0.0015552 0.0003894 0.0000976
3.5 0.0991784 0.0247222 0.0062033 0.0015543 0.0003894 0.0000978
4.0 0.0992288 0.0246753 0.0061938 0.0015537 0.0003900 0.0000985
Таблица 2 Порядок точности р для последовательных пар сеток
і 75 и 150 150 и 300 300 и 600 600 и 1200 1200 и 2400
0.5 2.053 1.994 1.997 1.999 1.999
1.0 2.051 1.994 1.997 1.999 1.999
1.5 2.046 1.994 1.997 1.999 1.999
2.0 2.038 1.994 1.997 1.998 1.999
2.5 2.026 1.994 1.997 1.998 1.998
3.0 2.010 1.994 1.997 1.998 1.996
3.5 2.004 1.995 1.997 1.997 1.993
4.0 2.008 1.994 1.995 1.994 1.985
Поведение погрешности при измельчении сетки проиллюстрировано рис. 9. Приведена абсолютная ошибка, но так как полная энергия данного течения в рассматриваемой области течения плавно изменяется от 0.97638 при £ = 0.5 до 0.96711 при £ = 4.0 (аналитические вычисления), то относительная ошибка ведет себя аналогичным образом. Анализ поведения полной энергии показывает, что выбранная разностная схема аппроксимирует полную энергию со вторым порядком точности и практически обладает свойством полной консервативности, так как с течением времени энергия диссипирует очень медленно.
0.10 І 0.05
а
° 0.00
0.5
І
Рис. 9. Поведение погрешности при измельчении сетки на последовательные моменты времени.
Заключение
В работе исследован численный метод решения задачи о накате длинных волн на берег, впервые предложенный автором в работе [15]. В последние годы метод неоднократно применялся для решения задач, связанных с численным моделированием распространения волн цунами. Сравнение данного метода с другими, описанными в литературе, показал, что построенный алгоритм является оптимальным в том смысле, что, с одной стороны, прост в применении и допускает очевидное обобщение на двумерный случай (для других методов это не всегда возможно либо требует существенных затрат), а с другой стороны, полученные с его помощью значения высоты заплеска (области затопления суши) хорошо коррелируют с результатами, полученными с помощью других более сложных методов.
Рис. 10. Обтекание конусообразного острова уединенной волной с заплеском на откос.
В заключение приведем два рисунка, иллюстрирующих применение метода к решению двумерных задач со сложной геометрией. На рис. 10 показано обтекание с накатом на сухой берег конического острова, а рис. 11 демонстрирует фрагмент заплеска на реальное побережье волны цунами.
Автор выражает благодарность д-ру ф.-м. н. Л. Б. Чубарову за постоянный интерес к работе и плодотворное сотрудничество при постановке и решении задач динамики поверхностных волн.
Рис. 11. Подтопление мыса волной цунами.
Список литературы
[1] ВольцингЕр Н.Е., Клеванный К. А., ПЕлиновский Е. Н. Длинноволновая динамика прибрежной зоны. Л.: Гидрометеоиздат, 1989.
[2] Long-wave runup models / H. Yeh, Ph. Liu, C. Synolakis (Eds). Singapore: World Sci. Publ., 1996.
[3] Железняк М. И., Сидорчук В. Н. К численному расчету трансформации волнения в прибрежной зоне // Волны в сплошных средах. К. 1978. С. 13-19.
[4] Egnesund L. Shallow-water equations in a water-course with sloping shores extended with an artificial water domain. Uppsala Univ., Goteborg, Department of Sci. Comput.
1989. Report No. 124.
[5] Кузьмичева Т. В., Новиков В. А., Федотова З. И. Сравнение некоторых методов сквозного счета при моделировании наката длинных волн на берег // Всесоюз. совещание по численным методам в задачах волновой гидродинамики, 23-25 сент.,
1990, Ростов-на-Дону / Под ред. Ю. И. Шокина. ВЦ СО АН СССР, Красноярск, 1991. С. 9-14.
[6] Завьялов В. К., Кузьмичева Т. В., Лаппо А. Д. и др. Численное и экспериментальное моделирование наката длинных волн на пологий берег с использованием линейно-нелинейной модели // Вычисл. технологии. 1992. Т. 1, №3. С. 205-212.
[7] Федотова З.И., Чубаров Л. Б. Численное моделирование наката цунами // Вычисл. технологии: Тр. Междунар. конф. RDAMM-2001. 2001. Т. 6, Ч. 2 (Специальный выпуск). С. 380-396.
[8] Смирнов В. И. Курс высшей математики. Т. 2. М.: Наука, 1974.
[9] Лятхер И. М., МилитЕЕв А. Н., Школьников С. Я. Расчет наката волн цунами на берега // Изучение цунами в открытом океане / Под ред. Е. Ф. Саваренского, С. Л. Соловьева. М.: Наука, 1978. С. 48-55.
[10] Рождественский Б. Л., Яненко Н. Н. Системы квазилинейных уравнений и их приложения к газовой динамике. М.: Наука, 1978.
[11] ГогодзЕ И. К., Попов Ю.П., Хуцишвили В. В. Непрерывные автомодельные и периодические решения уравнений мелкой воды // Накат цунами на берег / Под ред. Е. Р. Пелиновского. Горький: Изд-во ИПФ СО АН СССР, 1983. С. 64-74.
[12] Шокин Ю.И., ЯнЕнко Н. Н. Метод дифференциального приближения. Применение к газовой динамике. Новосибирск: Наука, 1985.
[13] Федотова З. И. Исследование устойчивости разностных схем с неперестановочными матрицами // Математические модели и методы решения задач механики сплошной среды. Красноярск: ВЦ СО АН СССР, 1986. С. 88-96.
[14] Carrier G. F., Greenspan H. P. Water waves of finite amplitude on a sloping beach // J. Fluid Mech. 1958. No. 4. P. 97-109.
[15] Fedotova Z. I. The simple numerical method for a long wave on a beach // XlX-th Biennial Symp. on on Advanced Problems and Methods in Fluid Mech. Abstracts on papers. Warszawa, 1989. P. 200-201.
Поступила в редакцию 14 февраля 2002 г., в переработанном виде — 14 июня 2002 г.
