Прямое численное моделирование аттракторов внутренних волн стратифицированной жидкости в трапециедальной области с колеблющейся вертикальной стенкой Текст научной статьи по специальности «Физика»

Прямое численное моделирование

аттракторов внутренних волн стратифицированной жидкости в трапециедальной области с колеблющейся вертикальной стенкой

К. Брузе 1 < Christophe.brouzet&.ens-lvon.fr > Т. Доксуа 1 <thierry.dauxois(a),ens-lyon.fr > Е. Ерманюк 1,2 <ermanyuk(a>gm ail.com> С. Жубо 1 <sylvain.joubaud(a>.ens-lyon.fr> M. Крапошин 3,4 <os-cfd(a)yandex.ru> II. Сибгатуллин,4,5,6 <ilias_s(a)maH.ru> 1 laboratoire de Physique de I Ecole Normale Supérieure de lyon, Université de lyon, France ' Институт гидродинамики им. MA. Лаврентьева, Новосибирск, Россия 3Национальный исследовательский центр "Курчатовский институт",

Москва, Россия

4 Институт системного программирования РАН, Москва, Россия 5Механико-математический факультет и институт механики МГУ им. MB. Ломоносова, Москва, Россия 6 Институт океанологии им. П.П. Ширшова РАН, Москва, Россия

Список обозначений

s - абсолютная солёность, —

^ - скорость среды, м/с Р - плотность среды, мЗ/кг

Р - статическое давление среды. Па

Р - пьезометрическое давление. Па

- плотность солёной воды, кг/мЗ У - ускорение свободного падения, м/с2

г - радиус-вектор, м ^ - динамическая вязкость, Па-с v - кинематическая вязкость, м2/с

5 - коэффициент объёмного сужения жидкости за счёт солёности, ~

— число Шмидта, ~ МКО — Метод Конечного Объёма МСЭ — Метод Спектральных Элементов

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

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

1. Введение

Внутренние гравитационные волны в однородно стратифицированной жидкости при постоянной частоте плавучести могут распространяться под фиксированным углом к вертикали, который определяется только частотой вынуждающей силы [1]. При наличии наклонённой стенки волны могут фокусироваться и это приводит к образованию определённых путей, которые определяются лишь частотой. Эти пути были названы волновыми аттракторами. Необходимо отметить, что к аттракторам в фазовом пространстве из теории динамических систем данное понятие прямого отношения не имеет. Первые эксперименты, подтверждающие существование предсказанных ранее волновых аттракторов описаны в [1]. Существование подобных аттракторов и анализ их структуры, устойчивости и динамики важным для понимания процессов перемешивания в морях и озёрах, а также для анализа колебаний во вращающихся системах, таких как ядро Земли, поскольку инерционные волны также могут распространяться по аттракторам. Устойчиво стратифицированное состояние жидкости обладает частотой плавучести, или частотой Брента-Вяйсяля, которая характеризует малые колебания смещённого из исходного состояния элемента жидкости. Возмущения исходного состояния с постоянной частотой, меньшей частоты плавучести, приводят к внутренним волнам, распространяющимся вдоль прямых линий (в приближении постоянной частоты плавучести). Лучи могут отражаться от поверхностей, при этом после отражения лучи по прежнему

распространяются под углом, равным по модулю арккосинусу отношения вынуждающей частоты и частоты плавучести. За счёт такого отличия в характере отражения становится возможной фокусировка внутренних волн в замкнутой области при определенной геометрии. Замкнутые пути, возникающие в результате фокусировки, были названы волновыми аттракторами [1]. В частности, в работе [1] было показано, что такая фокусировка возможна, если одна из стенок области наклонена по отношению к вертикали. Также при формировании аттрактора в реальных условиях играет роль баланс фокусировки и вязкости [6]. Анализ двумерного около критично го отражения слабонелинейных внутренних гравитационных волн от наклоненной границы в однородностратифицированной жидкости приведен в [2,5]. По существу волновые аттракторы являются двумерными, что показано как теоретически, так и экспериментально [1-5].

Помимо теоретических и экспериментальных исследований волновых аттракторов, предпринимались попытки и их численного исследования: в работе [4] с помощью программного кода общей модели циркуляции (MIT général circulation model), использующей метод конечного объёма, исследовались волновые аттракторы, возникающие при гармоническом вертикальном колебании всей области. Проведено сравнение с экспериментом для устойчивого аттрактора. Известно, что методы конечного объёма обладают собственной численной вязкостью и их применение для анализа устойчивости сильно нелинейных режимов на больших временах требует осторожности. Помимо, этого, в [4] численные исследования проводились для числа Шмидта 100, а не 770, как в эксперименте, из-за невозможности разрешения диффузионных масштабов, а на нижней и правой стенке ставилось условие проскальзывания вместо прилипания. Таким образом численная модель не вполне соответствовала эксперименту. Тем не менее, при определённых параметрах удалось воспроизвести общий вид аттрактора. В работах [10,11] авторами проведены предварительные расчеты методом спектральных элементов (МСЭ).

В работе [7] для смежной задачи для инерционных волн в прямоугольной геометрии проведено прямое численное моделирование аттракторов вращающейся жидкости спектральным методом. Спектральный метод позволяет с высокой точностью рассчитывать нелинейные режимы и предпочтителен по сравнению с конечно объёмным для исследования нелинейных режимов, но его применение в непрямоугольной геометрии с разнородными граничными условиями затруднено.

Таким образом, при проведении и планировании расчетных исследований течений подобного рода мы сталкиваемся со следующей дилеммой:

• либо использовать метод конечного объёма, основными преимуществами которого являются простота реализации расчётной модели и устойчивость, но при этом внести в результаты расчётов ошибку, связанную с высокой численной диффузией самого метода;

• либо использовать метод спектральных элементов, порядок аппроксимации которого намного выше, чем в МКО (2-ой порядок), при этом потеряв гибкость и простоту создания расчётных моделей.

В этом случае хорошим решением может быть тестирование обоих методов на простой геометрии области моделирования с целью выявления характерных различий между МКО и МСЭ. А в следующих расчётных исследованиях использовать МКО и накопленный опыт для минимизации его недостатков, иными словами — численной диффузии. В качестве таких задач могут использоваться экспериментальные исследования, схожие с описываемыми в работах [3,7]. Эффективное применение МКО к задачам распространения внутренних волн малой амплитуды и движения тел в стратифицированных средах было реализовано в [13].

Другой принципиальной сложностью, например как в случае [3,7], является сопоставление с экспериментом расчётных данных. Это связано с некоторой неопределенностью постановки граничных и начальных условий, заданием физических констант — например, частоты плавучести, распределения плотности в начальный момент времени, числа Шмидта, амплитуды колебания поверхности исследуемого пространства и т. д. Более того, при интерпретации различных режимов течения с помощью расчётного анализа необходимо наличие критериев, показывающих границы этих режимов. Поиск таких критериев с помощью численного моделирования связан с обработкой большого количества данных.

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

С учётом указанных выше вызовов, исследование целесообразно проводить в следующем порядке:

1) Анализ условий экспериментального стенда для последующего численного моделирования;

2) Выбор математической модели;

3) Поиск приближений и упрощений математической модели, сформулированной в п. 2) для реализации численной модели;

4) Реализация численной модели течения методом конечного объёма (МКО) и методом спектральных элементов (МСЭ);

5) Сравнение МКО и МСЭ для выбранных расчётных случаев;

6) Исследование динамических характеристик моделей выбранных режимов матричными методами;

7) Обобщённый анализ полученных результатов.

2. Физическая постановка задачи и её математическое описание

Рассматривается двумерное течение однофазной несжимаемой вязкой жидкости — раствора соли в поле силы тяжести. В замкнутом объёме жидкости с заданной стратификацией в начальный момент времени, в движение приводится одна из вертикальных стенок (правая) — см. рис. 1-3. Вертикальная стенка «слева» движется в горизонтальном направлении по заданному закону, стенки снизу и справа — неподвижны, горизонтальная граница сверху — поверхность раздела, на которой задаётся условие «проскальзывания».

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

Рисунок 1. Принципиальная схема экспериментальной установки. Жирным пунктиром показана область моделирования

со

*

к I

ги aj

2 н и

QJ

к

CL ГО

S I

2 £

СЬ X

О m

QJ ч: о

tt с

v A

Поверхность раздела (непроницаемая стенка с проскальзыванием)

DH(V,t)=acos(jTV/H)cos((D0t) \

Неподвижная стенка

Н

Рисунок 2. Принципиальная схема расчётной области и физических внешних границ

0,0012

0,001

1 0,0008 X

О)

3"

§ 0,0006 (С

=с

н 0,0004

с;

с

< 0,0002 0

50 100 150 200 250 Физическое время, с

300

350

400

Рисунок 3. Схема «вывода» амплитуды смещения левой вертикальной стенки на

максимальное значение.

Для описания расчётной модели используются следующие уравнения сохранения:

уравнение сохранения массы

|£+v-(p&)=o

(1)

уравнение сохранения импульса

дрй дг

+ У-[ри®и)-У-о =-Ур+рд (2)

уравнение переноса массы растворённой соли

дг

(3)

Слагаемое в правой части уравнения (3) позволяет учесть диффузионные эффекты, которые могут возникать на длительных временах моделирования.

С

Для учёта влияния массовой концентрации растворённой соли на плотность жидкости используется выражение:

Р(5) = Ро+(|г)(5"5о) (4)

Где в — массовая концентрация солёности жидкости — отношение массы соли в элементарном объёме к общей массе среды.

Жидкость является Ньютоновской и касательные напряжения ° подчиняются линейному закону:

а = ^уи + (УиУ) (5)

3. Описание расчетной модели

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

а) пренебречь изменениями плотности в уравнении неразрывности (1);

б) пренебречь изменениями плотности в левой части баланса импульса (2).

Введём следующие величины:

Пьезометрическое давление Р ~Р ~Р 9'г Коэффициент объёмного сужения жидкости за счёт солёности

R Ifép

ls~v Us" PUs

Кинематическая плотность среды

Тогда система уравнений (1)-(4) может быть переписана в виде:

У-[> = 0 (6)

DU dt

V -(и О и) - V -(v(ví7 + (VE/У)) = - V-^— -(g ■ г)

) Vp

(7)

^- + V-(sÜ) = V~Vs (8)

sí se

Кинематическая плотность вычисляется как ^ Ps(s so) g

Ps

данном исследовании принимается, что величина s , также как и кинематическая вязкость v и число Шмидта считаются заданными для каждого расчётного случая и не меняются.

4. Численная модель с использованием метода конечного объема

Для численной реализации математической модели, описываемой уравнениями (6)-(8) используется метод конечного объёма второго порядка точности по времени и пространству и схема с расщеплением переменных (алгоритм PISO). Использование метода конечного объёма предполагает, что балансные соотношения импульса (7), массы соли (8) и объёма системы (6) должны быть сформулированы в интегральном виде.

Поскольку балансные соотношения составлены для подвижной контрольной массы, необходимо выполнить переход к контрольному объёму. Для этого используется транспортная теорема Рейнольдса (теорема дифференцирования интеграла по переменному объёму) - для любой экстенсивной величины 4х характеризуемой интенсивной величиной ^ , содержащейся в контрольной

массе CM(t) движущейся со скоростью ^" , изменение может быть вычислено как сумма изменения Ч* , содержащейся в подвижном объёме

^WO , охватывающем контрольную массу CM{t) в момент времени [

\\?ÜrdS = ii>(tf0— Üb)dS

плюс относительный поток этой величины r ^ °

через границы контрольного объёма при его движении относительно

контрольной массы:

dW d л ^ Л Л

-=— j рц/dV =— j рц/dV + j ptyiu ■ ds} ^-g-j

dt dt ,.,( , , dt „-( , , „,-( , ,

Или, переписав первое слагаемое правой части с учётом движения границ контрольного объёма ^'' :

d¥ d , г друг _ - _

-=— j j -+V-( PVUjdV + j pV(u-dS)

dt dt ,:M(t) lT(,) dt ,)

В соответствии с (9) или (10) возможно два подхода:

1) Положить контрольный объём постоянным в пространстве (метод

Эйлера) - . Этот подход удобен либо для случаев, в которых

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

2) В тех же случаях, когда скорость движения контрольной массы может быть направлена нормально по отношению к поверхности контрольного объёма, но при этом невозможно выделить выделить поверхности с известными параметрами среды на границе системы, целесообразно принять поверхность контрольного объёма движущейся вместе с поверхностью системы (совмещённый подход Лагранжа-Эйлера).

В соответствии с физической постановкой задачи, в данной работе используется второй подход, поскольку левая стенка исследуемого объёма движется в нормальном направлении. В этом случае уравнения (6)-(8) приводятся к виду (11)-(13):

J Ua-dS = 0

dCV{t)

j- J UadV+ J t7a(t7r-dS) =

Gl CV(t) dCV(t)

J (v(Vua+(VuJ)\ds- J ^+(g-T)pkds

dCV(t) dCV(t)

4 J sdV+ J s(Ur-dS) = J f-(vs-ds)

al CVit) dCVit) dCVit) ^C

(id

(12)

(13)

5. Процедура численного интегрирования уравнений математической модели

В соответствии с выбранным численным методом, был реализован алгоритм интегрирования уравнений (11)-(13). На каждом шаге по времени выполнялись следующие действия:

1) Вычисление нового положения границ расчётной области, обновление геометрии сеточных линий

2) Прогноз поля солёности по имеющимся полям объёмных потоков (полю скорости с предыдущего шага)

3) Обновление поля кинематической плотности

4) Прогноз поля скорости жидкости

5) Решение уравнения для давления

6) Обновление поля объёмных потоков, обновление поля скорости

7) Обновление поля солёности в соответствии с новым полем скорости

8) Переход к новому шагу по времени

При реализации численной модели применялись следующие численные «трюки»:

1) Учёт небаланса импульса на первом шаге по времени

Поскольку в начальный момент времени жидкость предполагается находящейся в состоянии покоя, необходимо выполнение следующего условия (см. (12)):

Это условие можно обеспечить двумя способами: либо найти такое распределение пьезометрического давления, которое в начальный момент времени будет компенсировать столб жидкости, либо же вычесть этот столб перед началом интегрирования:

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

2) Учёт градиента солёности в уравнении для давления

Чтобы получить уравнение для давления, воспользуемся известной процедурой PISO, представив уравнение сохранения импульса (12) в полудискретном виде

дС

У(ли -Я(с/)) = - | V——Л'- | £ ■ - £ ■ г О = 0) О = о)

Р

3-01 Р,

с1У

Где А — вклад в диагональ линеаризованного уравнения сохранения импульса Н — сумма вкладов от недиагональных элементов и источников

Перейдя к средне интегральным значениям и сократив на объём V, получаем:

*

(Н(иа)) (д.гУ9к-д.г(1=0)Урк(1 = 0))

и„

А

А

А

После интерполяции на грани ячеек и скалярного умножения на площади граней, получаем выражение для объёмных потоков вещества:

/ * \

Ф, =

Н [7,

\

А

ху ро а

[д-гУРк-д-г{1 = 0)УРк{1 = 0) А

И

Представив уравнение неразрывности (11) в дискретном виде через потоки

Ф г

= (и.

получаем уравнение для давлениям

I

(.Н(иа)> (д-гУрк-д-г(1=0)Урк(1=0)

А

А

Н

* \

1

\

А

Из решения уравнения для давления получаем величину поправки к объёмным потокам, необходимую для точного выполнения уравнения неразрывности (11). Получившаяся схема предполагает хранение поля скорости в центрах контрольных объёмов и потоков, обеспечивающих консервативность — в центрах граней.

3) Граничное условие для давления на внешних границах расчётной области

На границах расчётной области, в которых скорость среды задана постоянной (в данном случае — стенка), следует задать такое граничное условие, которое обеспечит нулевую коррекцию к потокам вещества

/ £- + (д.г)рк^~ / (д-г(г=0))рк(г = 0)<К = 0

асу(()

асу(()

4) Движение сеточных линий

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

Для учёта деформации расчётной области в каждый момент времени решается уравнение Лапласа для смещений граней расчётной области:

/ УтиУ-иьду = о (14)

6. Описание расчётной области и граничных условий

При разработке расчётной модели и расчётной области (постановке задачи) преследовалась цель максимально возможного соответствия условиям физического эксперимента (см. рис. 1). На рисунке 1 представлена схема экспериментальной установки, на рисунке 2 — соответствующая ей расчётная область. Принято, что поверхность раздела фаз можно считать непроницаемой стенкой, на которой реализуется условие проскальзывания скорости среды. Правая боковая и нижняя стенки непроницаемы и неподвижны, левая боковая стенка движется по синусоидальному закону (см. рис. 2). Следующим важным предположением является характер выхода на установившийся режим колебаний вертикальной стенки — как следует из рисунка 3, весьма вероятно, что максимальная амплитуда колебаний стенки на начальных этапах процесса (до 200 секунд) изменялась от 0 до номинального значения. Подтвердить это предположение на основе экспериментальных невозможно, поскольку датчики ускорения на стенке отсутствовали. Косвенным подтверждением может быть сравнение режимов с постепенным выходом на номинальную амплитуду и мгновенным приложением смещений в нулевой момент времени. В соответствии с этим, в расчётах амплитуда колебаний стенки изменялась от 0 до номинального значения по ступенчатому закону, показанному на рисунке 3. Данный закон был подобран на основе экспериментальных данных о динамике горизонтальной составляющей скорости в точке отбора (см. рис. 4).

Рис. 4 Экспериментальные данных о динамике горизонтальной составляющей скорости в точке отбора, по горизонтали отложены периоды колебаний левой

стенки.

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

Граничные и начальные условия

Граничные и начальные условия задавались в соответствии с постановкой задачи.

Пьезометрическое Давление:

На левой, нижней и правой стенках и верхней стенках — условие непроницаемости.

Скорость:

На левой, нижней и верхней стенках — условие прилипания. На верхней — условие проскальзывания.

Солёность:

На всех стенках условие непротекания (нулевой градиент)

Скорость смещения расчётной области:

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

горизонтальной 129

Uh (v, t )=a cos (Jt v/H)(— ш)sin(w01)

/

на правой вертикальной стенке — постоянное нулевое значение на верхней и нижней — нулевой градиент

Начальные условия

Скорость — жидкость находится в покое (все компоненты равны 0) Пьезометрическое давление — слои жидкости находятся в равновесии — равно давлению на верхней границе

Скорость смещения сеточных линий — все компоненты поля скорости смещения равны нулю

Солёность — распределена по заданному линейному закону от SO на верхней границе до максимального значения на нижней границе

У 1 s dz Po dz

Вместо абсолютной солёности воспользуемся относительной, которая меняется от линейно от 0 (верхняя граница) до 1 (нижняя граница). Задавшись плотностью жидкости на верхней плоскости, скажем 1000кг/мЗ, вычисляем плотность на нижней поверхности. Вычисляем значение

коэффициента s . N = 1.059 рад/с

g=9.81

р0 = ЮОО

pnN2/g = 114.32018 Соответствующий градиент плотности и

A z = 0.3 А р=34.296

ß, = —гг = —34.296 Х103 ls A sA Р°

Угловая частота 0.623 рад/с

Амплитуда колебаний 0.15см (1.5мм). Как показали дальнейшие расчётные исследования, величина 0.15см слишком большая для получения устойчивого аттрактора на промежутке времени 0-1500с, поэтому расчётным путём была подобрана амплитуда 0.105см

Положение точки отбора (от нижнего левого угла): +28.2см, +20см

7. Описание расчетной модели методом спектральных элементов.

При описании расчетной модели с использованием метода спектральных элементов [9,12] удобнее пользоваться дифференциальной формой законов сохранения. Уравнения сохранения импульса и диффузии соли в приближении Буссинеска представляются в виде:

^ + укУку) = -V— + V АУ + вд (15)

ОЬ ) рт

(д + укук3) = ЛЛв' (16)

сЦу(г>) = 0 (17)

Поскольку левая граница совершает гармонические колебания вида

жь(0, = а соз(7ту/Н) соз(и;о£),

граничные условия в силу малости амплитуды можно записать

и{0,2/, Ь) = аио соб(7гу/Н) V = 0

Остальные параметры аналогичны параметрам при описании численной модели методом конечного объема. Спектрально-элементный подход сочетает свойства конечно-элементного и спектрального методов. Расчетная область разбивается на элементы, в которых уравнения сохранения дискретизируются с помощью вариационного подхода, по сути являющимся Галеркинским методом. Для разнородных граничных условий удобно использовать базисные функции Лежандра при пространственной дискретизации элементов Гауса-Лобатто [9,12].

При расчетах методом спектральных элементов визуально профиль скорости аналогичен профилю скорости, полученному методом конечного объема:

а = 0.15 ст,Ы = 1.059газ/з,а;о = 0.623гай/в

Рисунок 5. Поле модуля скорости расчетов с помощью метода спектральных

элементов

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

Рисунок б. Поле модуля скорости и значение горизонтальной составляющей скорости

х = 28.2cm, у =20ст.

п mnwp '

Рисунок 7. Экспериментальное распределение поля модуля скорости, полученное авторами [3] (слева) и соответствующее колебание горнзонтапъной компоненты скорости в точке х = 28.2см,у = 20см(справа). По горизонтапи отложены периоды

вынужденных колебаний.

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

В соответствии с поставленной задачей выполнено тестирование расчётной модели с использованием метода конечного объёма и метода спектральных элементов. Список режимов показан в таблице 1. Динамика горизонтальной скорости как наиболее показательный параметр для сравнения с экспериментом для базового случая показана на рис. 8, результаты тестирования сеточной сходимости — на рис. 9 и 10. На рис. 11 показано поведение аттрактора (горизонтальной скорости жидкости в выбранной точке) без учёта условий эксперимента. На рисунках 12-13 — распределение вертикальной и соответственно горизонтальной компоненты скорости вдоль линии Х=28.2см, на рис. 14-15 — поле модуля скорости. На рисунке 16 показаны точки времени, для которых выполнялось сравнялось МКО и МСЭ. Результаты сравнения МКО и МСЭ представлены на рис. 17. На рис. 18 показано изменение поля горизонтальной компоненты скорости среды во времени в соответствии с рис. 12-13.

Таблица 1. Список режимов для кросс-тестирования МКО и МСЭ

№№ Описание

1 Базовый расчётный случай (соответствует описанию в разделе 0), сетка 225x150

2 Как случай №1, но с сеткой 450x300

3 Как случай №2, но с шагом по времени 2.5Е-3

4 Как случай №3, но с сеткой 900x600

5 Как случай №1, но с отрицательной амплитудой (-0.105см)

Для этого случая выводится распределение горизонтальной и вертикальной компонент скорости на линии у=20см от нижнего края

6 Как случай №1, но с постоянной амплитудой (0.105см)

7 Как случай №1, но с постоянной амплитудой (-0.105см)

Все расчёты проводились с использованием вычислительных мощностей кластера платформы ишНиВ ИСП РАН (http://www.unihub.ru) и суперкомпьютера Ломоносов.

Рабочая обллсть

Случай 1

СИ

:

Сб. сен]. 13 1E:D9:40 2014

Бремя, с

Рисунок 8. Динамика горизонтальной составляющей скорости в точке (28.2, 20) для

базового случая (случай 1)

Рабочая обллсть

О.ООЭ-1

0.002-O.OOi,

iii

illii

| J

J .и, i-11 1 * N . | I

П1Ш1

Ш111Ш1

lift

- ■■!

'■r

IHiHIIIliK

ill:

2»

CS. cphi 1E:D9:40 20X4

400

Бремя,С

Рисунок 9. Сравнение динамики горизонтальной составляющей скорости среды в точке (28.2, 20) для трёх различных сеточных разрешений (225x150, 450x300, 900x600), время с 200с по 600с

)&&##<P""" "!"5(*+%%$''++() ###**+%(("

■350 460 470 464 490 500 510

Бремя, с

Рисунок 10. Сравнение динамики горизонтальной составляюгцей скорости среды в точке (28.2, 20) для трёх различных сеточных разрешений (225x150, 450x300, 900x600), время с 450с по 550с

Случай 7 I I

Время, с

Рисунок 11. Динамика горизонтальной составляюгцей скорости в точке (28.2, 20) для случая 7 (без постепенного вывода скорости движения подвижной стенки)

Расстояние от стенни слева, м

Рисунок 12. Профиль горизонтальной и вертикальной скоростей на линии у=20см для момента 501.5с (максимум горизонтальной скорости)

Рисунок 13. Профиль горизонтальной и вертикальной скоростей на линии у=20см для момента 506.5с (минимум горизонтальной скорости)

U Magnitude 0.008244 "0.008 -0,006 [0.004 »0.002

1,655е-6

Рисунок 14. Визуализация мгновенного поля модуля скорости момента 501.5с (максимум горизонтальной скорости)

Рисунок 15.

U Magnitude D.008244 '0.008 -0.006 10.004 f0.002

1.655е-6

Визуализация мгновенного поля модуля скорости момента 506.5с (минимум горизонтальной скорости)

Рабочая область

2

л" 3

Случай 1 Случай 2 Случай 4

сетка 225x150 сетка 450x300 сетка ЭООхбОСЬ

450 460

Сб. сент. 13 18:09:40 2014

480 490 500 510

Время, с

520 530 540 550

Рисунок 16. Сравнение сеточного разрешения и схема выбора точек для кросс-тестирования МКО и МСЭ

з: I-

о и

со о

^ о.

Яг °

о *

1— и

-Л

У/-

1=501.5с

1=506.5с

Рисунок 17. Сравнение МКО и МСЭ для различных моментов времени

— Ux-500 Ux-500,5—- Ux-501 — Ux-501,5 — Ux-502 — Ux-502,5

0 0, )5 0 1 0, 15 III / 2 / / / 0 2/ / 0 3 0,

Рисунок 18. Изменение поля скорости продольной скорости среды в зависимости от

времени

9. Выводы

Двумерные численные расчеты хорошо воспроизводят экспериментальные данные при учете изменения линейного профиля солености у свободной поверхности. При этом амплитуда колебаний в численных расчетах увеличена за счет потерь энергии волнопродуктора в трехмерном контейнере. При проведении эксперимента предполагалось создание линейного профиля, но добиться идеально линейной стратификации у верхней границы не представляется возможным. При учете изменения профиля солености у верхней границы форма аттрактора в расчете и эксперименте совпадает. Несмотря на то, что общий вид аттрактора хорошо воспроизводится как методом спектральных элементов, так и методом конечных объемов, в профилях скоростей имеются отличия у левой границы. Это требует дальнейшего исследования, поскольку такие отличия могут влиять на нелинейную динамику волн при развитии неустойчивостей. В расчётах использовались вычислительные мощности платформы UniHUB ИСП РАН (http://www.unihub.ru) и суперкомпьютерных кластеров МГУ им. М.В.Ломоносова. Работы выполнялись при финансовой поддержке Министерства образования и науки Российской Федерации (уникальный идентификатор соглашения RFMEFI60714X0090) и гранта РФФИ 15-01-06363.

Список источников

[1]. Maas, L. R. М. & Lam, F.-Р. АЛ Geometrie focusing of internal waves. J. Fluid Mech, 1995,. 300, 1^1

[2]. Dauxois, Thierry; Young, W.// Journal of Fluid Mechanics, 1999, vol. 390, Issue 01, p.271-295

[3]. Scolan, H., Ermanyuk, E., Dauxois, T.// 2013, Physical Review Letters, 110, 234501 138

[4]. Grisouard, N., Staquet, С., Pairaud, I.//2008, Journal of Fluid Mechanics, 614, 1

[5]. Hazewinkel, J., van Breevoort, P., Dalziel, S.~B., Maas, L.~R.~M.// 2008, Journal of Fluid Mechanics, 598, 373

[6]. Frans-Peter A. Lam, Leo R.M. Maas. Internal wave focusing revisited; a reanalysis and new theoretical links // Fluid Dynamics Research 40 (2008) 95 - 122.

[7]. Mercier, Matthieu J.; Gamier, Nicolas В.; Dauxois, Thierry Reflection and diffraction of internal waves analyzed with the Hilbert transform Physics of Fluids, Volume 20, Issue 8, pp. 086601-086601-10 (2008).

[8]. Jouve, L., Ogilvie, G.~I.// Journal of Fluid Mechanics, 2014, 745,223.

[9]. Fischer P., RonquistE. Spectral element methods for large scale parallel Navier— Stokes calculations, Computer Methods in Applied Mechanics and Engineering, vol. 116, issue 1-4, pp. 69-76, 1994

[10]. Сибгатуллин И.Н. Моделирование волнового аттрактора в стратифицированной жидкости. Международная конференция "Mode Conversion, Coherent Structures and Turbulence", Space Research Institute of RAS, Moscow, Russia, 2014

[11]. Brouzet C., Dauxois Т., Ermanyuk E., Kraposhin M., Sibgatullin I. Modelling of Wave Attractors in Stratified Fluids, 5-я международная научная школа молодых ученых «ВОЛНЫ И ВИХРИ В СЛОЖНЫХ СРЕДАХ», Москва, 2014

[12]. Fischer P.F. An overlapping schwarz method for spectral element solution of the incompressible Navier-Stokes equations. J. Comput. Phys. 133 (1), 84-101, 1997.

[13]. Загуменный Я.В., Чашечкнн Ю.Д. Расчет течений непрерывно стратифицированной жидкости с использованием открытых вычислительных пакетов на базе технологической платформы UniHUB. Труды Института системного программирования РАН, том 24, 2013 г. стр. 87-106.

Direct numerical simulation of internal gravity wave attractor in trapezoidal domain with oscillating vertical wall

C. Brouzet1 < Christophe.brouzet(a),ens-lyon.fr> T. Dauxois1 < thierry.dauxois(ciiens-lyon.fr>

E. Ermanvuk1" <ermanyuk(a)gmail.com> S. Joubaud1 <sylvain.joubaud(a),ens-lyon.fr> M. Kraposhin 3,4 <os-cfd(q)yandex.ru> I. Sibgatullin'4'5'6 <ilias_s(a),mail.ru> 1 Laboratoire de Physique de l'École Normale Supérieure de Lyon, Université de Lyon, France ' Lavrentyev Institute of Hydrodynamics, Novosibirsk, Russia 3 P National Research Center "Kurchatow institute", Moscow, Russia 4 Institute of system programming of Russian Academy of Sciences, Moscow, Russia 3 Faculty of Mechanics and Mathematics, and Institute of Mechanics of Moscow

State University, Russia 6 Shirshov Oceanology Institute of Russian Academy of Sciences, Moscow, Russia

Abstract. Direct numerical simulation of internal gravity waves focusing and developement of a wave attractor was performed with the help of two different numerical approaches. Mathematical formulation corresponds to experiments on exitations of inner waves in a trapezoidal container with salt solutions through forced oscillations of the left boundary. It was shown that numerical simulations reproduce the experiments after taking into account ther imperfection of linear salinity profile near the upper boundary. The amplitudes of resulting oscillations in both numerical simulations are increased as compared to the experiments due to loss of energy of the 3D wave generator in the experiments. Despite the fact that the general shape of the attractor is reproduced by both method, there are differences in velocity profiles near the left boundary. This fact requires further investigations since this discrepancy may in influence nonlinear dynamics of developing instabilities.

Keywords: attractor, internal waves, gravity waves, DNS, direct numerical simulation.

References

[1]. Maas L. R. M. & Lam F.-P. A.// Geometric focusing of internal waves. J. Fluid Mech, 1995,. 300, 1^1

[2]. Dauxois Thierry, Young W.// Journal of Fluid Mechanics, 1999, vol. 390, Issue 01, p.271-295

[3]. Scolan H., Ermanyuk E., Dauxois T.// 2013, Physical Review Letters, 110, 234501

[4]. Grisouard N., Staquet C., Pairaud I.// 2008, Journal of Fluid Mechanics, 614, 1 140

[5]. Hazewinkel J., van Breevoort P., Dalziel S.B., Maas L.~R.~M.// 2008, Journal of Fluid Mechanics, 598, 373

[6]. Frans-Peter A. Lam, Leo R.M. Maas. Internal wave focusing revisited; a reanalysis and new theoretical links // Fluid Dynamics Research 40 (2008) 95 - 122.

[7]. Mercier Matthieu J., Gamier Nicolas B., Dauxois Thierry Reflection and diffraction of internal waves analyzed with the Hilbert transform Physics of Fluids, Volume 20, Issue 8, pp. 086601-086601-10 (2008).

[8]. Jouve L., Ogilvie GUI Journal of Fluid Mechanics, 2014, 745, 223.

[9]. Fischer P., RonquistE. Spectral element methods for large scale parallel Navier— Stokes calculations, Computer Methods in Applied Mechanics and Engineering, vol. 116, issue 1-4, pp. 69-76, 1994

[10]. Sibgatullin I.N. Modeling of wave attractor in strafitied fluid. International conference "Mode Conversion, Coherent Structures and Turbulence", Space Research Institute of RAS, Moscow, Russia, 2014

[11]. Brouzet C., Dauxois T., Ermanyuk E., Kraposhin M., Sibgatullin I. Modelling of Wave Attractors in Stratified Fluids, 5-th international school of young scientists «Waves and vorticesin complex media», Moscow, 2014

[12]. Fischer P.F. An overlapping schwarz method for spectral element solution of the incompressible Navier-Stokes equations. J. Comput. Phys. 133 (1), 84-101, 1997.

[13]. Zagumennyi Ia.V. , Chashechkin Yu.D. Calculations of continuously stratified fluid flows using open source computational packages based on the technological platform UniHUB. Papers of Institute of System Programming of Russian Academy of Sciences, vol 24, 2013, pp. 87-106.