Известия РАН. Физика атмосферы и океана, 2020, T. 56, № 3, стр. 360-372

Трехмерное численное моделирование морских волн

Д. В. Чаликов ab*

a Институт океанологии им. П.П. Ширшова РАН
117997 Москва, Нахимовский пр., 36, Россия

b University of Melbourne
3010 Victoria, Australia

* E-mail: dmitry-chalikov@yandex.ru

Поступила в редакцию 24.10.2019
После доработки 25.12.2019
Принята к публикации 05.02.2020

Полный текст (PDF)

Аннотация

Даются характеристики численных моделей поверхностных волн, основанных на полных уравнениях движения жидкости со свободной поверхностью в потенциальном приближении. Рассматриваются их возможности и области применения. Более подробно описана модель, использующая связанную с поверхностью систему координат. Уравнение Лапласа для потенциала в такой системе координат превращается в полное эллиптическое уравнение, которое решается как уравнения Пуассона с итерациями по правой части. Предложен метод решения основанный на разделении решения на аналитическое и отклонение от него. Этот подход значительно ускоряет вычисления. Задача в целом решается для периодической по обоим горизонтальным направлениям области Фурье методом с расчетом нелинейных членов на густой сетке. Проведено моделирование развития волнового поля под действием ветра и диссипации. Показано, что трансформация основных интегральных характеристик решения удовлетворительно согласуется с известными данными. Перечисляются возможные приложения развиваемого подхода.

Ключевые слова: моделирование, ветровые волны, приток энергии к волнам, диссипация энергии, развитие волн под действием ветра

1. ВВЕДЕНИЕ

Термины, употребляемые в данной статье, заслуживают пояснения. Волны на поверхности жидкости могут быть одномерными (однонаправленными) или двумерными, т.е. распространяющиеся в разных направлениях. Одномерные волны инвариантны по отношению к поперечному сдвигу и поэтому они описываются двумерными уравнениями, решаемыми в вертикальной плоскости $\left( {x - z} \right).$ В конформных координатах [1] вертикальная структура решения известна, поэтому допустимо называть конформную модель одномерных волн также одномерной моделью. Двухмерное волновое поле описывается трехмерными уравнениями, и пока никому не удалось свести эту задачу к двухмерной. Метод поверхностного интеграла [2] оперирует только с поверхностными переменными, но поскольку метод основан на использовании функции Грина, его трудно отнести к двухмерному.

Первые модели по воспроизведению движения жидкости со свободной поверхностью на основе некоторой разновидности Лагранжева подхода были сформулированы в работах [37]. Метод был основан на прослеживании положения поверхности на фиксированной сетке с различной степенью точности. Этот метод оказался все же недостаточно точным, поскольку его можно было применять лишь к коротким отрезкам времени. Тем не менее, метод ценен тем, что он может быть использован для исследования нестационарного потенциального движения со свободной (не обязательно однозначной) поверхностью на основе уравнений для двухфазной жидкости (см. [8]).

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

Сложность моделирования волн зависит от принятых предположений. Условие периодичности существенно упрощает формулировку модели и численную схему благодаря тому, что в этом случае можно использовать представление Фурье, обеспечивающее точное дифференцирование. Мы не склонны приписывать модам Фурье высокий физический смысл, поэтому рассматриваем их как один из способов интерпретации результатов. Преобразование Фурье используется, главным образом, как способ численного решения уравнений. Периодичность переменных можно интерпретировать как свойство объекта малого по сравнению с большим водоемом, имеющим медленно меняющиеся по пространству характеристики. Периодичность позволяет использовать точный сеточный-Фурье (Fourier-transform) метод, что дает возможность вести интегрирование на длительные сроки и исследовать трансформацию волнового поля за счет нелинейных взаимодействий, притоков и стоков энергии. Для непериодических объектов, таких, как небольшие замкнутые бассейны или малые части больших бассейнов с известными граничными условиями, глобальное Фурье представление неприменимо, и численная схема должна использовать конечные разности. В этом случае возможен учет специфических граничных условий (поглощение, отражение, частичное отражение), реальной батиметрии и донных процессов. В целом, такие модели не могут быть очень точными. Тем не менее, этот подход дает гораздо более подробные результаты, чем спектральные модели.

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

В настоящее время наиболее популярной моделью является HOS (High Order Scheme, [19, 20]) модель, основанная на работах [21, 22] и формулировке уравнений, предложенной Захаровым [23]. В модели используются две координатные системы: следующая поверхности криволинейная координата для кинематического и динамического граничных условий и декартова система координат, в которой известно аналитическое решение уравнения Лапласа для потенциала скорости. Переход от одной системы координат к другой происходит с использованием ряда Тейлора в декартовой системе. Точность этого метода и его эффективность зависят от числа используемых членов ряда Тейлора. Для волн малой амплитуды или узкого волнового спектра точность метода может быть очень высока. Моделирование волн с широким спектром требует много членов ряда Тейлора, что, разумеется, может быть причиной вычислительной неустойчивости. Трудность возникают из-за того, что возмущения потенциала скорости для малых волн, распространяющихся по поверхности больших волн, затухают с глубиной очень быстро на глубине, сопоставимой с высотой большой волны. Моделирование волнового поля с помощью модели HOS с коротким отрезком ряда Тейлора осуществляется благополучно, поскольку происходит фильтрация коротких волн, что повышает устойчивость схемы. Модель HOS широко используется в исследованиях (см. например [2426]). Недавно (Ecole Centrale Nantes, LHEEA Laboratory CNRС) было объявлено, что HOS модель стала доступной для общего использования [27].

Другая группа моделей, в противоположность HOS модели, основана на прямом решении трехмерного уравнения для потенциала скорости, записанного в следующей поверхности криволинейной системе координат. Наиболее универсальный подход успешно развивается в Датском Техническом Университете (Technical University of Denmark, TUD) (см. [28]). ModelWave3D, созданная в TUD, предназначена для исследования широкого круга проблем, включая проблему взаимодействия волн с плавающими и погруженными объектами. Благодаря своей универсальности, модель может использоваться для океанологических исследований, например, для воспроизведения волнового режима в небольших водоемах с реальной формой и батиметрией. Влияние рельефа дна введено использованием так называемой сигма-координаты, выпрямляющей дно и свободную поверхность. На береговой линии используются граничные условия поглощения или отражения или смешанные условия. Сравнение TUD модели с HOS моделью для сравнительно простых волновых процессов показало хорошее согласование.

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

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

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

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

Важную роль в исследовании волновых процессов играет сравнение лабораторных измерений в волновых или ветроволновых каналах с результатами моделирования. Этот метод, например, использовался в [24] для исследования неустойчивости и разрушения одиночных волн. В работе [21] решение на основе интеграла Коши сравнивалось с точными измерениями в канале, и получено очень хорошее совпадение.

2. ОСНОВНЫЕ РЕЗУЛЬТАТЫ ПРЯМОГО МОДЕЛИРОВАНИЯ

2.1. Одномерное моделирование

Нестационарное конформное преобразование для конечной глубины позволяет переписать уравнения потенциальных волн в следующей поверхности координатной системе. Благодаря конформности уравнение Лапласа для потенциала скорости приобретает стандартное представление, основанное на разложении Фурье на поверхности. В результате, система уравнений без упрощений превращается в систему двух эволюционных уравнений, которые можно решать традиционными методами, в частности, методом Фурье. В этом случае предположение о потенциальности упрощает проблему столь кардинально, что уравнения можно решать без применения конечных разностей, когда пространственные производные вычисляются аналитически, а нелинейные члены – на густой сетке с оцениваемой точностью. Для ограниченного порядка нелинейности этот метод также является точным. Модель дает пример редкого в гидродинамике случая, когда уравнения описывают реальный процесс с высокой достоверностью. Это утверждение точно, если крутизна поверхности в декартовых координатах не обращается в бесконечность. Теоретически в случае односвязной области конформное преобразование существует, и поверхность остается однозначной, хотя для поддержания точности требуется быстро возрастающее число мод. Практически оказывается, что поверхность может наклониться на угол, превышающий 90°. После этого начинается физическая неустойчивость, связанная с падением жидкости, не поддерживаемой давлением. Если не принимаются особые искусственные меры (см., например, [29]), возникает быстрая неустойчивость гребня [30], проявляющаяся в конечном итоге разделением на две фазы. Легко понять, что эти процессы не потенциальны. Следовательно, если вычисления желательно продолжить, то, как и во многих разделах геофизической гидродинамики, необходимо принять меры для предотвращения вычислительной неустойчивости и параметризовать процессы, связанные с явлением обрушения. После этого модель теряет прозрачность и математическую строгость, но оказывается приближенной к реальным процессам.

Первым результатом, полученным с конформной моделью, можно считать разработку методов расчета стационарных решений уравнений потенциального движения со свободной поверхностью: гравитационных волн на глубокой воде (волн Стокса, [1]), гравитационных волн на мелкой воде, гравитационно-капиллярных волн и капиллярных волн (волн Краппера, [31]). Эти результаты изложены в работах [1, 32]. Разработанные ранее многочисленные методы расчета стационарных решений основывались на плохо сходящихся разложениях и были очень неэкономичны. Схемы, основанные на конформном преобразовании, позволяют рассчитывать стационарные решения с неограниченной точностью и скоростью на два порядка выше, чем традиционные схемы. Существование трех видов стационарных решений для полных уравнений потенциального движения со свободной поверхностью представляет уникальную возможность проверки нестационарных моделей. Для расчетов в качестве начальных условий задается одно из стационарных решений, и далее проводится интегрирование по времени с полной стационарной моделью. Если модель (включая численную схему и программирование) верна и стационарное решение устойчиво к возмущениям, заданная волна должна двигаться с верной фазовой скоростью, не изменяя своей формы.

В последнее время возникло мнение, что волны Стокса являются не просто интересным математическим объектом и средством проверки численных моделей. Они начинают играть важную роль в исследовании физики поверхностных волн. Традиционный подход к исследованию нелинейных свойств многомодового волнового процесса основан на предположении, что реальное волновое поле можно представить как суперпозицию линейных волн. Такой подход, используемый в знаменитой работе Хассельманна [33], лежит в основе моделирования нелинейной трансформации волнового спектра. Этот существенно спектральный подход предполагает, что фазы всех мод отвечают линейному дисперсионному соотношению и распределены равномерно и случайно. В действительности фазовые скорости некоторых мод равны фазовой скорости более крупных волн, поскольку волновое поле является смесью свободных волн и присоединенных мод. Присоединенные моды не являются волнами, а всего лишь “кирпичиками”, используемыми для построения нелинейной формы (поэтому термин “окаймляющие” волны (bound waves) затемняет суть дела). Даже визуальные наблюдения морской поверхности подтверждают, что реальные волны заострены, а их подошвы сглажены.

Принято считать, что быстрый (по сравнению с хассельмановским) механизм трансформации волнового спектра предложен теорией модуляционной неустойчивости, известной в оригинале как теория неустойчивости Бенджамина-Фейера [34]. Суть этой теории весьма проста: в одномерном цуге волн развиваются новые моды, что в конечном итоге приводит к развитию плотного волнового спектра. Согласно результатам авторов теории и численных расчетов [35] развитие новых мод происходит сравнительно медленно. Общепринятый сценарий развития экстремальной волны предполагает, что одна из мод спектра по неизвестным причинам начинает отбирать энергию из спектрального окружения, медленно растет и, в конце концов, достигает очень больших значений. Никаких данных о том, что такой процесс реализуется в природе или лаборатории, не существует. Численные эксперименты с точной трехмерной моделью, наоборот, показывают, что формирование экстремальной волны происходит очень быстро – за время порядка одного периода.

Численное исследование процесса опрокидывания волн [37] показало, что ни одна из динамических, кинематических или геометрических характеристик не может быть надежным предиктором опрокидывания волны, так что эволюция этих характеристик характеризует само опрокидывание, но не объясняет его механизм. Между тем, даже визуальные наблюдения показывают, что опрокидыванию всегда предшествует заострение волновых гребней. Этот эффект нельзя объяснить гребневой неустойчивостью [30], которая относится лишь к волнам предельной крутизны, тогда как заострение гребней, сопровождаемое фокусировкой энергии, происходит даже при небольшой крутизне.

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

2.2. Трехмерное моделирование

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

С точки зрения физики проблему прямого (фазоразрешающего) моделирования можно разделить на две группы: адиабатическое и не адиабатическое моделирование. Простая адиабатическая модель предполагает отсутствие притока и стока энергии. Этот подход, не будучи полностью свободным от ограничений, позволяет исследовать волновые процессы с помощью точных исходных уравнений. Учет эффектов притока энергии и ее диссипации чаще всего связан с гипотезами, противоречащими предположению о потенциальности, т.е. к уравнениям добавляются члены, которые можно отнести к феноменологическим. По этой причине формулировки не адиабатического подхода часто оказываются различными.

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

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

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

Основные типы трехмерных моделей перечислены в предисловии. Здесь несколько более подробно будет рассмотрена модель, наиболее удачно на наш взгляд сочетающая сравнительную простоту с приемлемой эффективностью. Модель [32, 41] предназначена для интегрирования трехмерных полных уравнений движения жидкости в потенциальном приближении со свободной поверхностью с периодическими граничными условиями в горизонтальных направлениях. Модель многократно обсуждалась в публикациях, поэтому здесь дается лишь краткое ее описание. Используется связанная с поверхностью $\eta \left( {x,y,t} \right)$ неортогональная система координат:

(1)
$\xi = x,\,\,\,\,\vartheta = y,\,\,\,\,\zeta = z - \eta (\xi ,\vartheta ,\tau ),\,\,\,\,\tau = t,$
где $\eta (x,y,t)$ = $\eta (\xi ,\vartheta ,\tau )$ – движущаяся, однозначная периодическая поверхность. Уравнения движения имеют вид:
(2)
${{\eta }_{{\tau }}} = - {{\eta }_{{\xi }}}{{\varphi }_{{\xi }}} - {{\eta }_{\vartheta }}{{\varphi }_{\vartheta }} + \left( {1 + \eta _{{\xi }}^{2} + \eta _{\vartheta }^{2}} \right){{\Phi }_{{\varsigma }}},$
(3)
${{\varphi }_{{\tau }}} = - \frac{1}{2}\left( {\varphi _{{\xi }}^{2} + \varphi _{\vartheta }^{2} - \left( {1 + \eta _{{\xi }}^{2} + \eta _{\vartheta }^{2}} \right)\Phi _{{\zeta }}^{2}} \right) - \eta - p,$
(4)
${{\Phi }_{{{\xi \xi }}}} + {{\Phi }_{{\vartheta \vartheta }}} + {{\Phi }_{{{\zeta \zeta }}}} = \Upsilon \left( \Phi \right),$
где $\Upsilon $ – оператор:

(5)
$\begin{gathered} \Upsilon \left( {} \right) = 2{{\eta }_{{\xi }}}{{\left( {} \right)}_{{{\xi \zeta }}}} + 2{{\eta }_{\vartheta }}{{\left( {} \right)}_{{\vartheta {\zeta }}}} + \\ + \,\,\left( {{{\eta }_{{{\xi \xi }}}} + {{\eta }_{{\vartheta \vartheta }}}} \right){{\left( {} \right)}_{{\zeta }}} - \left( {\eta _{{\xi }}^{2} + \eta _{\vartheta }^{2}} \right){{\left( {} \right)}_{{{\zeta \zeta }}}}, \\ \end{gathered} $

заглавная буква $\Phi $ используется для области $\zeta < 0,$ а строчная $\varphi $ относится к поверхности $\zeta = 0.$ Член $p$ описывает давление на поверхности $\zeta = 0.$ Уравнения записываются в безразмерном виде, формально следующим из предположения, что ускорение свободного падения. В работе [41] предлагается представить потенциал скорости $\Phi \;\left( \varphi \right)$ как суммы аналитической $\bar {\Phi },$ $\left( {\bar {\varphi } = \bar {\Phi }\left( {\xi ,\vartheta ,0} \right)} \right)$ и произвольной нелинейной $\tilde {\Phi },$ $\left( {\tilde {\varphi } = \tilde {\Phi }\left( {\xi ,\vartheta ,0} \right)} \right)$ компонент.

(6)
$\varphi = \bar {\varphi } + \tilde {\varphi },\,\,\,\,\Phi = \bar {\Phi } + \tilde {\Phi }.$

Аналитическая компонента удовлетворяет уравнению Лапласа:

(7)
${{\bar {\Phi }}_{{{\xi \xi }}}} + {{\bar {\Phi }}_{{\vartheta \vartheta }}} + {{\bar {\Phi }}_{{{\zeta \zeta }}}} = 0,$

с известным решением

(8)
$\bar {\Phi }(\xi ,\vartheta ,\zeta ,\tau ) = \sum\limits_{k,l} {{{{\bar {\varphi }}}_{{k,l}}}(\tau )\exp \left( {\left| k \right|{\zeta }} \right)} {{\Theta }_{{k,l}}},$

здесь $k$ и $l$ – волновые числа в направлениях $x$ и $y,$ $\left| k \right| = {{\left( {{{k}^{2}} + {{l}^{2}}} \right)}^{{1{\text{/}}2}}},$ ${{\bar {\varphi }}_{{k,l}}}$ – Фурье коэффициенты для поверхностного потенциала $\bar {\Phi }$at $\zeta = 0;$ ${{\Theta }_{{k,l}}}$ – базисные функции разложения Фурье. Нелинейная компонента потенциала удовлетворяет уравнению Пуассона

(9)
${{\tilde {\Phi }}_{{{\xi \xi }}}} + {{\tilde {\Phi }}_{{\vartheta \vartheta }}} + {{\tilde {\Phi }}_{{{\zeta \zeta }}}} = \Upsilon \left( \Phi \right),$

Уравнение (9) решается численно на растянутой по вертикали сетке в пространстве Фурье методом прогонки с граничными условиями:

(10)
$\begin{gathered} \varsigma = 0:\,\,\tilde {\Phi } = 0, \\ \varsigma \to - \infty :\,\,{{{\tilde {\Phi }}}_{\zeta }} \to 0. \\ \end{gathered} $

Преимущество разделения потенциала скорости на две компоненты состоит в том, что нелинейная компонента на 1–2 порядка меньше линейной, поэтому необходимая точность при решении уравнения (9) достигается при меньшем количестве итераций, чем требуется для полного потенциала скорости. В качестве начальных условий задаются поля возвышения и поверхностного потенциала, удовлетворяющие линейным уравнениям при заданном волновом спектре. Подробное описание модели дано в монографии [32].

Модель, описанная в упомянутых работах, в последнее время была усовершенствована путем небольших модификаций разностной схемы и схем параметризации притока энергии и диссипации. Адиабатическая версия последнего варианта модели была проверена путем моделирования движения волны Стокса почти предельной крутизны ${\text{АК}} = 0.40.$ Начальные условия для волны Стокса были рассчитаны по алгоритму, разработанному в работе [1].

На рис 1. показана эволюция амплитуд мод Стокса, рассчитанная по трехмерной модели. Как видно из рисунка, по крайней мере первые 15 мод (включая моды с амплитудой порядка ${{10}^{{ - 4}}}$) сохраняются с очень высокой точностью; моды с амплитудами ${{10}^{{ - 6}}} - {{10}^{{ - 4}}}$ сохраняют свою идентичность, и лишь моды с амплитудами $А < {{10}^{{ - 6}}}$ теряют упорядоченность. Скорость волны не считалась, но в предыдущих расчетах [41] было получено, что скорость совпадает с теоретической с пятью знаками. Заметим, что структура волны Стокса очень сложна: она состоит из бесконечного числа мод, каждая из которых движется со скоростью первой моды. Даже небольшие ошибки в формулировке задачи или кодах ведут к быстрому разрушению волны. Очевидно, что этого не происходит. Такая проверка модели является точной и не тривиальной. Заметим, что ни одна из известных моделей не подвергалась такой проверке.

Рис. 1.

Зависимость от времени $t$ (выраженном в периодах волны) амплитуд $А$ волн Стокса. Период приблизительно равен $2\pi .$

Первый цикл работ, основанных на модели (1–10), был посвящен исследованию адиабатических процессов в волновых полях, т.е. процессов, происходящих без притоков энергии и ее диссипации. В нелинейной гидродинамической системе с конечным числом степеней свободы всегда возникает поток энергии в сторону малых волновых чисел. Для предотвращения неустойчивости в этой области вводится селективный фильтр. Для поддержания состояния квазистационарности очень медленное изменение энергии Е (порядка $\left( {{{{10}}^{{ - 9}}} - {{{10}}^{{ - 7}}}} \right)Е$ за один шаг по времени) компенсируется интегральным притоком энергии. В таком режиме модель используется для исследования многих проблем волновой механики, в частности, проблемы статистики экстремальных волн.

Рис. 2.

Двухмерные волновые спектры в пространстве волновых чисел (правая ось соответствует k, левая – $l$). Верхняя секция – спектр, соответствующий начальным условиям. Нижняя секция – волновой спектр после интегрирования на 318 периодов волны пика. Спектры нормированы на максимальное значение невозмущенного спектра.

1. Было обнаружено, что статистика полной высоты волны (от подошвы до гребня) оказывается совершенно идентичной для линейной и нелинейной моделей (Chalikov, Babanin, 2016b). Нелинейность проявляется только в статистике высот волн и глубины подошвы, по отдельности. Этот факт до сих пор не получил объяснения.

2. Показано, что волновой спектр по мере своего развития приобретает сильную неоднородность: глубокие провалы в двухмерном спектре чередуются с хорошо выраженными пиками [41].

Интересно, что с увеличением разрешения в Фурье пространстве области спектра с низкой энергией расширяются. Можно предположить, что используемое разрешение было чрезмерным. Это утверждение выглядит слишком оптимистичным и малоправдоподобным. Тем не менее, его можно обосновать с помощью следующего рассуждения. Очевидно, многомодовая волновая динамика не похожа на многомасштабное турбулентное движение. Увеличение разрешения при моделировании турбулентности позволяет учитывать более мелкие структуры, но в разумных пределах не слишком значительно влияет на воспроизведение структуры потока в целом. В волновом поле с увеличением разрешения вводятся другие моды со своими фазовыми скоростями. Благодаря дисперсии, геометрия и эволюция поля с одинаковым спектром и удвоенным разрешением будут различными. Это означает, что численная волновая задача не удовлетворяет фундаментальному требованию сходимости при увеличении разрешения. Ситуацию может быть спасена, если существует предельное разрешение, при превышении которого некоторые моды не вносят своего вклада в движение из-за того, что их энергия слишком мала. Можно предположить, например, что очень близкие друг к другу моды не могут существовать по отдельности и быстро сливаются в одну, т.е., одна из мод сильно ослабляется, а вторая – увеличивается.

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

Если скорость ветра превышает фазовую скорость основных энергонесущих волн, энергия волн увеличивается, и волновой спектр смещается в стороны низких волновых чисел. Этот процесс многократно воспроизводился с помощью спектральных моделей с целью уточнения схем расчета скорости притока энергии от ветра и диссипации. Моделирование процесса развития волн с помощью фазоразрешающих моделей является намного более трудной задачей хотя бы потому, что расчет эволюции волн следует проводить, по крайней мере, на сотни периодов волны пика. Для существенно трехмерной проблемы, в которой на каждом шаге решаются уравнения Пуассона с итерациями, подобные расчеты осуществляются в течение десятков часов машинного времени. Для решения такой задачи к уравнениям (110) должны быть присоединены алгоритмы, описывающие приток энергии и импульса и диссипации энергии. Испытания и настройка таких схем для нестационарной и громоздкой задачи являются довольно трудоемким процессом. Надо признать, что, несмотря на огромное количество работ в области физики морских волн, состояние этой проблемы является далеко не удовлетворительным. Например, феномен диссипации волн в виде обрушения является объектом двухфазной гидродинамики. Общепризнанных надежных алгоритмов параметризации этого явления не существует. Проблема взаимодействия волн и ветра является более сложной, чем сама проблема моделирования волн, поскольку она требует проведения чрезвычайно сложных измерений или численного моделирования взаимодействия воздушного потока с многомодовым волновым полем [43, 44].

Обмен энергией и импульсом между водой и воздухом осуществляется неоднородным распределением динамического давления на поверхности ${{p}_{0}}$ в (3). Согласно общепринятому мнению, комплексные Фурье-компоненты давления ${{p}_{k}}$ и ${{p}_{{ - k}}}$ на поверхности связаны с компонентами возвышения ${{h}_{k}}$ и ${{h}_{{ - k}}}{\text{:}}$

(11)
${{p}_{k}} + i{{p}_{{ - k}}} = \frac{{{{\rho }_{a}}}}{{{{\rho }_{w}}}}\left( {{{\beta }_{k}} + i{{\beta }_{{ - k}}}} \right)({{h}_{k}} + i{{h}_{{ - k}}}),$
где ${{\beta }_{k}}$ и ${{\beta }_{{ - k}}}$ – реальный и мнимый коэффициенты $\beta $-функции, введенной Майлсом [45]. В настоящее время принято считать, что $\beta $‑ф-ункция зависит от единственного параметра $\Omega = {{\omega }_{k}}U\cos \psi $ = = ${U \mathord{\left/ {\vphantom {U {{{c}_{k}}\cos \psi }}} \right. \kern-0em} {{{c}_{k}}\cos \psi }}$ (где ${{\omega }_{k}}$ и U – безразмерные частота и скорость ветра; ${{c}_{k}}$ – фазовая скорость k-той моды; $\psi $ – угол между направлениями ветра и этой моды. Большинство схем, рекомендованных для вычисления этой функции, основаны на сильно упрощенных аналитических формулах (см., например, [46]). Линейное предположение (11) уже в значительной степени упрощено (например, не учтена зависимость коэффициента $\beta $ от амплитуды), поэтому использование дальнейших упрощений ничем не оправдано. Исследование функции $\beta $ должно основываться на одновременном воспроизведении полей возвышения $\eta $ и поверхностного давления. Такaя исключительно сложная работа уже проводилась в натурных условиях экспериментально с помощью датчика давления в воздухе, отслеживающего поверхность [47, 48]. Объем полученных таким образом данных не позволяет пока построить функцию $\beta $ в достаточно широком диапазоне. Второй метод исследования функции $\beta $ заключается в численном моделировании взаимодействия волн и ветра [49]. Этот метод позволяет получить очень большой объем данных, необходимых для установления функции $\beta ,$ но качество этих исследований, разумеется, зависит от полноты и точности формулировки задачи, особенно формулировки проблемы замыкания уравнений Рейнольдса. Первый вариант алгоритма расчета потока энергии к волнам был введен в прогностическую модель WAVEWATCH (см. [45, 46]). Следующая версия модели, описывающая взаимодействие пограничного слоя с волнами [43], дала результаты, заслуживающие наибольшего доверия, благодаря еще и тому, что они получены обработкой огромного объема данных, охватывающих широкий диапазон безразмерной частоты $\Omega .$ Функция $\beta ,$ полученная в этих расчетах, использована в моделировании развития волн [50].

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

Основным процессом диссипации волновой энергии является обрушение волн [51]. В той или иной форме обрушение волн учитывается во всех моделях волновых процессов, в частности, в прогностических волновых моделях типа WAVEWATCH. Поскольку сами волны в таких моделях отсутствуют, никаких локальных критериев обрушения не может быть сформулировано, поэтому диссипация волн представлена в спектральных моделях в искаженном виде. Обрушение волн происходит на сравнительно узких участках. В спектральных моделях спектр обрушения распространяется практически на весь спектральный диапазон, т.е. ослабляются все волны, в том числе – короткие. На самом деле обрушение понижает высоту и энергию главных волн. Противоречие возникает из-за того, что волны в спектральных моделях предполагаются линейными, тогда как на самом деле обрушивается нелинейная заостренная волна.

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

Механизм обрушения волн при реальном развитом спектре сильно отличается от механизма обрушения в волновом поле, представленном несколькими модами, как это предполагается в теоретических и лабораторных исследованиях. В известной степени, обрушение волны аналогично развитию экстремальной волны, которая тоже появляется внезапно, без заметной предыстории. В обоих явлениях нет ни малейших признаков модуляционной неустойчивости, выражающейся в росте линейной волны, отбирающей энергию у других мод. Ситуация с обрушением и возникновением экстремальных (“freak”) волн скорее обратная: энергия индивидуальной волны сохраняется, а колумнарная энергия концентрируется вокруг гребня волны, которая становится более заостренной и неустойчивой. Именно эти соображения послужили основой выбора в качестве критерия обрушения степень кривизны волнового профиля (см. [41, 50]).

Описание обрушения в прямом моделировании должно удовлетворять трем условиям: (1) предотвращать опрокидывание и следующую за ним вычислительную неустойчивость; (2) реалистично описывать потери локальной потенциальной и кинетической энергии; (3) сохранять объем. Параметризация обрушения может основываться на локальном высокоселективном диффузионном операторе с подходящим коэффициентом диффузии. Такой алгоритм не изменяет объем и понижает кинетическую и потенциальную энергию волн. Предполагается, что теряемый импульс переносится к горизонтальному потоку, а энергия – частично к турбулентности и течению [52]. Кроме того, энергия и импульс могут переноситься к другим волновым модам. Выбор параметров в схеме основывается на простых рассуждениях: волны должны максимально близко приблизиться к обрушению, но не обрушиться, а подвергнуться сглаживанию.

Алгоритмы притока энергии и диссипации были включены в модель (1–10), после чего модель была инициирована для длительного счета. Для расчетов использовалась сетка $2048 \times 1024$ узлов, число вертикальных уровней для расчета нелинейной компоненты трехмерного потенциала равнялось 10. Использовалась растянутая сетка $\Delta {{z}_{{i + 1}}} = \gamma \Delta {{z}_{i}},$ где $\gamma = 1.2$ – коэффициент растяжения. Глубина жидкости $H$ принималась равной $2\pi {{k}_{p}},$ где ${{k}_{p}}$ – текущее значение волнового числа пика спектра. В качестве начальных условий задавался спектр JONSWAP [53] с максимумом на волновом числе ${{k}_{p}} = 100.$ Начальные волны были почти однонаправленными: угловое разрешение в области энергосодержащих волн было принято пропорциональным ${{\left( {\cos \psi } \right)}^{{256}}}$ ($\psi $ – угол между направлением волновой моды и направлением ветра). Расчеты проводились на 1 200 000 шагов по времени с шагом $\Delta t = 0.005.$

Полная удельная энергия волнового поля $Е = {{E}_{p}} + {{E}_{k}}$ (${{E}_{p}}$ – потенциальная энергия, ${{E}_{k}}$ – кинетическая энергия) рассчитывалась по следующим формулам:

(12)
$\begin{gathered} {{E}_{p}} = 0.25\overline {{{\eta }^{2}}} ,\,\,\,\,{{E}_{k}} = 0.5\overline{\overline {\left( {\varphi _{x}^{2} + \varphi _{y}^{2} + \varphi _{z}^{2}} \right)}} , \\ E = {{E}_{p}} + {{E}_{k}}, \\ \end{gathered} $
где одна черта обозначает осреднение по координатам $\xi $ и $\vartheta ,$ а двойная черта – осреднение по всему объему. Уравнение эволюции интегральной энергии E имеет вид:
(13)
$\frac{{dЕ}}{{dt}} = \bar {I} + \overline {{{D}_{b}}} + \overline {{{D}_{t}}} + \bar {N},$
где $\bar {I}$ – скорость притока энергии к волнам от ветра, $\overline {{{D}_{b}}} $ – скорость диссипации энергии за счет обрушения волн, $\overline {{{D}_{t}}} $ – скорость диссипации за счет высокочастотной фильтрации, $\bar {N}$ – интегральный эффект нелинейных взаимодействий, описываемый правыми частями уравнений (2, 3), при условии, что поверхностное давление $p$ равно нулю. Дифференциальную форму компонент трансформации энергии можно вывести из уравнений (2–5), но это пока никому не удавалось. Более простой и свободный от ошибок аппроксимации метод состоит в вычислении прироста энергии различными группами вне основного цикла вычислений с применением фиктивного шага по времени. Например, значение притока энергии вычисляется по соотношению:
(14)
$\bar {\bar {I}} = \frac{1}{{\Delta t}}\left( {\overline{\overline {{{E}^{{t + \Delta t}}}}} - \overline{\overline {{{E}^{t}}}} } \right),$
где $\overline{\overline {{{E}^{{t + \Delta t}}}}} $ – интегральная энергия, полученная после одного шага по времени с правой частью уравнений (2, 3), содержащей только поверхностное давление.

Эволюция интегральных характеристик волнового поля показана на рис. 3a. Быстрое изменение всех характеристик при $t < 500$ объясняется приспособлением первичных линейных полей к нелинейности. Интегральный эффект нелинейных взаимодействий $\bar {\bar {I}}$ (прямая линия 1) оказался очень близок к нулю. Высокочастотная диссипация $\overline{\overline {{{D}_{t}}}} $ (кривая 2) значительно слабее, чем диссипация за счет разрушения $\overline{\overline {{{D}_{t}}}} $ (кривая 3). Величина $\overline{\overline {{{D}_{b}}}} $ обнаруживает сильные колебания, поскольку обрушение обладает свойством перемежаемости. Полная диссипация ${{\bar {\bar {D}}}_{b}} + {{\bar {\bar {D}}}_{t}}$ поглощает почти всю поступающую энергию, лишь малая часть которой расходуется на увеличение энергии. Баланс энергии $\bar {\bar {B}} = \bar {\bar {I}} + {{\bar {\bar {D}}}_{t}} + {{\bar {\bar {D}}}_{b}}$ (кривая 5) к концу интегрирования приближается к нулю.

Рис. 3.

a – Скорости изменения интегральной энергии волн, умноженные на 109 за счет: 1 – нелинейных взаимодействий; 2 – диссипации в высокочастотной части спектра; 3 – обрушения волн; 4 – притока энергии. Кривая 5 описывает баланс всех составляющих. Серые отрезки показывают мгновенные значения, сплошные кривые – результат сглаживания. б – зависимость интегральных характеристик от разгона $F$: 1 – взвешенная частота в максимуме спектра; 2 – частота максимума спектра; 3 – полная энергия волн, умноженная на 106; 4 – аппроксимация (17).

Дополнительные характеристики решения представлены на секции 3б. Кривая 1 описывает эволюцию взвешенной частоты пика ${{\omega }_{w}}{\text{:}}$

(15)
${{\omega }_{w}} = {{\left( {\frac{{\int {kSdkdl} }}{{\int {Sdkdl} }}} \right)}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}},$
где интеграл берется по всей Фурье области. Величина ${{\omega }_{w}}$ не зависит от деталей спектра, поэтому ее эволюция хорошо воспроизводит развитие (downshifting) спектра. Кривая 2 описывает эволюцию частоты максимума спектра. Ступенчатая форма этой кривой обнаруживает фундаментальное свойство: спектр смещается не монотонно, от волнового числа к ближайшему волновому числу, а скачками, в данном случае – на 3–4 волновых числа путем ослабления предыдущего максимума и роста нового. Очевидно, это тесно связано с дискретностью спектра, показанной на рис. 1. Любопытно, что такое явление наблюдалось со спектральной моделью [54].

Величину разгона $F$для периодической задачи можно оценить интегрированием по времени величины фазовой скорости ${{c}_{p}} = {{\left| k \right|}^{{ - 1/2}}}$ в пике спектра

(16)
$F = \int\limits_{{{t}_{0}}}^t {{{c}_{p}}dt} .$

Проведенные численные эксперименты по развитию волнового поля под действием постоянного ветра (за вычетом несообразностей, вводимых условием периодичности). соответствуют полевому эксперименту JONSWAP [53]. Обработка результатов эксперимента показали, что частота спектрального пика ${{\omega }_{p}}$ уменьшается как ${{F}^{{ - 1/3}}},$ а полная энергия растет с разгоном $F$линейно. Ни одна из этих зависимостей не может быть точной: первая из них сингулярна при $F = 0,$ и обе зависимости не учитывают приближение к режиму насыщения. В нашем случае более точной формулой, описывающей зависимость частоты ${{\omega }_{p}}$ от разгона, является

(17)
${{\omega }_{p}} = \frac{{75.6}}{{5.63 + {{F}^{{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-0em} 3}}}}}}.$

(кривая 4). Формула ${{\omega }_{p}} \sim {{F}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 3}} \right. \kern-0em} 3}}}}$ справедлива в более узком интервале F. Разумеется, коэффициенты в (17) зависят от скорости ветра. Линейная зависимость энергии E от разгона F вообще не подтверждается, поскольку модель воспроизводит гораздо более протяженный интервал разгона, чем тот, что наблюдался в эксперименте.

3. ЗАКЛЮЧЕНИЕ

Фазоразрешающие (прямые) трехмерные модели являются мощным средством исследования свойств реальных волн. В противоположность теоретическим исследованиям, основанным на упрощенных уравнениях, результаты, полученные на компьютере, могут рассматриваться почти как наблюдения в лаборатории или полевых условиях, с одним ценным исключением: численные данные практически лишены шума и помех, всегда имеющих место в реальном мире. Модель позволяет “наблюдать” практически любые характеристики волнового поля с высоким разрешением в пространстве и времени. Моделирование волн отставало от других разделов геофизической и, в особенности, технической гидродинамики. Численное моделирование можно уподобить совершенному лабораторному моделированию с широкими возможностями вариации параметров и возможности “наблюдения” любых свойств объекта. Поскольку математическое моделирование свободно от различий в масштабах, точность воспроизведения естественных процессов выше, чем точность физического моделирования. Эта ситуация типична для других разделов гидродинамики, особенно для технических приложений. Зачастую оказывалось, что численные методы обеспечивают более высокую точность, нежели эксперимент. Это относится также к волновым экспериментам в волновых каналах. Очень часто подобные эксперименты посвящены исследованиям нелинейных волновых процессов. Здесь хотелось бы сделать замечание. Предположим, что мы имеем двухмерную модель волн и располагаем доказательствами, что уравнения решаются точно. Мы пытаемся повторить физический эксперимент в волновом канале и получаем один из следующих результатов:

а) Довольно хорошее совпадение с его результатами. Если отмести извечное желание авторов представить свои результаты в лучшем виде, можно сделать вывод, что физический эксперимент проводился в соответствии с потенциальной теорией в пределах сопоставляемых параметров.

б) Совпадение получилось неудовлетворительным. Можно искать ошибки в кодах. Если таких ошибок нет, это означает, что имеется несоответствие модели и эксперимента. Если модель хорошо проверена, усовершенствовать ее трудно. Следовательно, эксперимент может быть и хорош, но он не соответствует модели. Проблема состоит в том, что исследователь никогда не знает точно характеристики волн, поскольку генерируется поверхность, а вторая переменная, потенциал скорости, остается неизвестной. При этом контроль над потенциалом, отсутствует. Взаимное приспособление “верной” поверхности к неизвестному потенциалу сопровождается появлением дополнительных и ненужных мод, которые ухудшают или вовсе уничтожают согласование. Неопределенность возрастает, когда исследуются многомодовые или распределенные по углу процессы. Относительно малые размеры волновых каналов препятствуют исследованию процессов на больших временных и пространственных масштабах. Поэтому для исследования экстремальных волн или процессов обрушения волн часто создают искусственно сильное взаимодействие волн со специально подобранными амплитудами и фазами.

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

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

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

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

Работа выполнена в рамках государственного задания РАН России (тема 0149-2019-0015). Раздел 2.2 выполнен при частичном финансировании грантом РФФИ (№ 18-05-01122). Автор благодарит О.И. Чаликову за помощь в оформлении статьи.

Список литературы

  1. Chalikov D., Sheinin D. Direct modeling of one-dimensional nonlinear potential waves // Nonlinear ocean waves. In: Advances in fluid mechanics (Perrie W. (ed)). 1997. V. 17. P. 207–258.

  2. Craig W., Sulem C. Numerical simulation of gravity waves // J. Comput. Phys. 1993. V. 108. P. 73–83.

  3. Harlow F.H., Welch E. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface // Phys. Fluids. 1965. V. 8. P. 2182–2189.

  4. Noh W.F., Woodward P. SLIC (simple line interface calculation). // In: Lecture Notes in Physics. Springer, New York. 1976. V. 59. P. 330–340.

  5. Hirt C.W., Nichols B.D. Volume of fluid method for the dynamics of free surface // J. Comput. Phys. 1981. V. 39. P. 201–225.

  6. Prosperetti A., Jacobs J.W. A numerical method for potential flow with free surface // J. Comp. Phys 1983. V. 51. P. 365–386.

  7. Miyata H. Finite-difference simulation of breaking waves // J. Comput. Phys. 1986. V. 5. P. 179–214.

  8. Iafrati A. Numerical Study of the Effects of the Breaking Intensity on Wave Breaking Flows. J. Fluid Mech. 2009. V. 622. P. 371–411.

  9. Zhao X., Liu B.-J., Liang S.-X., Sun Z.-C. Constrained Interpolation Profile (CIP) method and its application // J. Ship Mechanics. 2016. V. 20. P. 393–402.

  10. Young D.-L., Wu N.-J., Tsay T.-K. Method of Fundamental Solutions For Fully Nonlinear Water Waves // Advances in Numerical Simulation of Nonlinear Water Waves. 2010. P. 325–355.

  11. Ma Q.W., Yan S. Qale-FEM method and its application to the simulation of free responses of floating bodies and overturning waves // Advances in Numerical Simulation of Nonlinear Water Wave. 2010. P. 165–202.

  12. Beale J.T. A convergent boundary integral method for three-dimensional water waves // Math. Comput. 2001. V. 70. № 335. P. 977–1029.

  13. Xue M., Xu H., Liu Y., Yue D.K.P. Computations of fully nonlinear three-dimensional wave and wave–body interactions. I. Dynamics of steep three-dimensional waves // J. Fluid Mech. 2001. V. 438. P. 11–39.

  14. Grilli S., Guyenne P., Dias F. A fully nonlinear model for three-dimensional overturning waves over arbitrary bottom // Int. J. Num. Math. Fluids. 2001. V. 35. P. 829–867.

  15. Clamond D., Grue J. A fast method for fully nonlinear water wave dynamics // J. Fluid. Mech. 2001. V. 447. P. 337–355

  16. Clamond D., Fructus D., Grue J., Krisitiansen O. An efficient method for three-dimensional surface wave simulations. Part II: Generation and absorption // J. Comp Phys. 2005. V. 205. P. 686–705.

  17. Fructus D., Clamond D., Grue J., Kristiansen O. An efficient model for three-dimensional surface wave simulations. Part I: free space problems // J. Comp. Phys. 2005. V. 205. P. 665–685.

  18. Fochesato C., Dias F., Grilli S. Wave energy focusing in a three-dimensional numerical wave tank // Proc. R. Soc. 2006. V. A462. P. 2715–2735.

  19. Ducrozet G., Bonnefoy F., Le Touzé D., Ferrant P. 3D HOS simulations of extreme waves in open seas // Nat. Hazards Earth Syst. Sci. 2007. V. 7. P. 109–122. https://doi.org/10.5194/nhess-7-109-2007

  20. Ducrozet G., Bingham H.B., Engsig-Karup A.P., Bonnefoy F., Ferrant P.A. Comparative study of two fast nonlinear free-surface water wave models // Int. J. Numer. Meth. Fluids. V. 2012. V. 69. P. 1818–1834.

  21. Dommermuth D., Yue D. A high-order spectral method for the study of nonlinear gravity waves // J. Fluid Mech. 1987. V. 184. P. 267–288.

  22. West B., Brueckner K., Janda R., Milder D., Milton M.R. A new numerical method for surface hydrodynamics // J. Geophys. Res. 1987. V. 92. P. 11803–11824.

  23. Zakharov V.E. Stability of periodic waves of finite amplitude on the surface of deep water // J. Appl. Mech. Tech. Phys. 1968. V. 9. P. 190–194.

  24. Tanaka M., Dold J.W., Lewy M., Peregrine D.H. Instability and breaking of a solitary wave. J. Fluid Mech. 1987. V. 187. P. 235–248.

  25. Toffoli A., Onorato M., Bitner-Gregersen E., Monbaliu J. Development of a bimodal structure in ocean wave spectra // J. Geophys. Res. 2010. V. 115. № C3. https://doi.org/10.1029/2009JC005495

  26. Touboul J., Kharif C. Two-Dimensional Direct Numerical Simulations Of The Dynamics Of Rogue Waves Under Wind Action // Advances in Numerical Simulation of Nonlinear Water Waves. 2010. P. 43–74.

  27. Ducrozet G., Bonnefoy F., Le Touzé D., Ferrant P. HOS-ocean: Open-source solver for nonlinear waves in open ocean based on High-Order Spectral method. Comp. Phys. Comm. https://doi.org/10.1016/j.cpc.2016.02.017

  28. Engsig-Karup A.P., Bingham H.B., Lindberg O. An efficient flexible-order model for 3D nonlinear water waves // J. Comp. Phys. 2009. V. 228. P. 2100–2118.

  29. Dold J.W. An efficient surface-integral algorithm applied to unsteady gravity waves // J. Comp. Phys. 1992. V. 103. P. 90–115.

  30. Longuet-Higgins M.S., Tanaka M. On the crest instabilities of steep surface waves // J. Fluid Mech. 1997. V. 336. P. 51–68.

  31. Crapper G.D. An exact solution for progressive capillary waves of arbitrary amplitude // J. Fluid Mech. 1957. V. 96. P. 417–445.

  32. Chalikov D. Numerical modeling of sea waves, Springer. 2016. 330 p. https://doi.org/10.1007/978-3-319-32916-1

  33. Hasselmann K. On the non-linear energy transfer in a gravity wave spectrum, Part 1 // J. Fluid Mech. 1962. V. 12. P. 481–500.

  34. Benjamin T.B., Feir J.E. The disintegration of wave trains in deep water // J. Fluid Mech. 1967. V. 27. № 3. P. 417–430.

  35. Chalikov D. Simulation of Benjamin-Feir instability and its consequences // Phys. Fluid. 2007. V. 19. № 1. P. 016602–016615.

  36. Chalikov D., Babanin A.V. Simulation of wave breaking in one-dimensional spectral environment // J. Phys Oceanogr. 2012. V. 42. № 11. P. 745–1761.

  37. Johannessen T.B., Swan C. A numerical transient water waves – part 1. A numerical method of computation with comparison to 2-D laboratory data // Appl. Ocean. Res. 1997a. V. 19. P. 293–308.

  38. Johannessen T.B., Swan C. A laboratory study of the focusing of transient and directionally spread surface water waves // Proc. R. Soc. Lond. 1997b. V. A457. P. 971–1006.

  39. Johannessen T.B., Swan C. On the nonlinear dynamics of wave groups produced by the focusing of surface-water waves // Proc. R. Soc. Lond. 2003. V. A459. P. 1021–1052.

  40. Chalikov D., Babanin A.V. Nonlinear sharpening during superposition of surface waves // Ocean Dynamics 2016a. V. 66. № 8. P. 931–937.

  41. Chalikov D., Babanin A.V., Sanina E. Numerical modeling of three-dimensional fully nonlinear potential periodic waves // Ocean Dyn. 2014. V. 64. № 10. P. 1469–1486.

  42. Chalikov D., Babanin A.V. Simulation of one-dimensional evolution of wind waves in a deep water. Physics of Fluid. 2014. V. 26. № 9. P. 096697.

  43. Chalikov D., Rainchik S. Coupled numerical modelling of wind and waves and the theory of the wave boundary layer // Boundary-Layer Meteorol. 2010. V. 138. P. 1–41. https://doi.org/10.1007/s10546-010-9543-7

  44. Sullivan P.P., McWilliams J.C., Oatton E.G., Large-Eddy Simulation of Marine Atmospheric Boundary Layers above a Spectrum of Moving Waves // J. Atmos. Sci. 2014. V. 71. № 11. P. 4001–4027.

  45. Miles J.W. On the generation of surface waves by shear flows // J. Fluid Mech. 1957. V. 3. P. 185–204.

  46. Janssen P.A.E.M. Quasi-linear theory of wind-wave generation applied to wave forecasting. // J. Phys. Oceanogr. 1991. V. 21. № 11. P. 1631–1642.

  47. Donelan M.A., Babanin A.V., Young I.R., Banner M.L., McCormick C. Wave follower field measurements of the wind input spectral function. Part I. Measurements and calibrations // J. Atmos. Oceanic Tech. 2005. V. 22. P. 799–813.

  48. Donelan M.A., Babanin A.V., Young I.R., Banner M.L. Wave follower field measurements of the wind input spectral function. Part II. Parameterization of the wind input // J. Phys. Oceanogr. 2006. V. 36. P. 1672–1688.

  49. Chalikov D.V. Numerical simulation of the boundary layer above waves // Bound Layer Met. 1986. V. 34. № 1–2. P. 63–98.

  50. Chalikov D. Numerical modeling of surface wave development under the action of wind // Ocean Sci. 2018. V. 14. P. 453–470. https://doi.org/10.5194/os-14-453-2018

  51. Babanin A.V. Breaking and Dissipation of Ocean Surface Waves. Cambridge University Press. 2011. 463 c.

  52. Chalikov D., Belevich M. One–dimensional theory of the wave boundary layer // Boundary Layer Meteorol. 1993. V. 63. P. 65–96.

  53. Hasselmann K., Barnett R.P., Bouws E. et al. Measurements of wind-wave growth and swell decay during the Joint Sea Wave Project (JONSWAP) // Tsch. Hydrogh. Z. Suppl. 1973. V. A8. № 12. P. 1–95.

  54. Rogers W.E., Babanin A.V., Wang D.W. Observation-consistent input and whitecapping-dissipation in a model for wind-generated surface waves: description and simple calculations // J Atmos. Oceanic Tech. 2012. V. 29. № 9. P. 1329–1346.

Дополнительные материалы отсутствуют.