Чтобы пронаблюдать поведение частицы по градиенту плотности в стандартном SPH-методе и в GSPH, был рассчитан тест на подавление возмущений. Выбирается двухмерная рассчетная область [-Lx, Lx] x [-Ly, Ly]. Здесь Lx = Ly = π/2. Сначала частицы расположены в решетке со смещением Δх/2, а перепад плотности между верхним и нижним слоями составляет 1:2. К верхнему и нижнему слоям добавляется небольшое смещение положения ξ (ξx, ξy), соответственно:
(19)
(20)
Здесь х и у – первоначальное положение частиц (т.е. решетка). Амплитуда и волновое число смещения положения, А, устанавливаются равными 0.01 π и 2 соответственно. ξ добавляется к исходному положению, чтобы переместить частицу к возмущенному положению. Предполагается, что по всей области вычисления давление равномерно, а скорость звука в верхнем слое устанавливается равной 1. Время прохождения звука, tsc,в вертикальном направлении в верхнем слое становится равным π/2. На рис. 2 показаны снимки результатов, полученных в стандартном SPH-методе и GSPH в различные моменты времени.
Рисунок 2. Тест на подавление возмущений в стандартном SPH-методе (левая колонка) и GSPH (правая колонка). Отдельные снимки показывают положение частиц при t = 0.0, 1.3, 2.7 и 4.0tsc,,соответственно, сверху вниз. Результаты стандартного SPH-метода показывают, что первоначальное возмущение полностью подавлено в 4tsc вследствие отталкивания частиц по градиенту плотности, хотя давление при этом остается постоянным. А результаты GSPH, напротив, показывают, что возмущение сохранилось. Первоначальный перепад плотностей верхнего и нижнего слоев составляет 1:2.
Тест не предполагает какого-либо движения частиц, поскольку во всей области вычисления давление равномерно. Однако результаты, полученные при использовании стандартного SPH-метода, показывают, что первоначальная граница разрыва плотности сглаживается. А при использовании GSPH разрыв плотности достаточно хорошо сохраняет свою форму. Становится ясно, что отталкивание частиц вследствие несогласованности стандартного SPH-метода гасит возмущение. Сила отталкивания действует в нормальном направлении по отношению к разрыву плотности, создавая поверхностное натяжение (Price, 2008). Мы изменяли кривизну первоначального возмущения и утверждаем, что подавление зависит от кривизны.
Лагранжиан
Еще один способ получить уравнения SPH-метода – использование Лагранжиана (например, Price и Monaghan, 2004, и ссылки данной работы). Лагранжиан, L, для жидкости выглядит следующим образом:
(21)
где u – удельная внутренняя энергия. Функция Лагранжа для стандартного SPH-метода:
(22)
Для данной функции Лагранжа уравнение Эйлера-Лагранжа дает уравнение движения для стандартного SPH-метода. Однако, функция Лагранжа уже использует аппроксимацию частиц (уравнение 22), поэтому уравнение импульса, полученное из функции Лагранжа, по-прежнему сохраняет несогласованность в неравномерном распределении частиц.
Взаимосвязь между этой функцией Лагранжа и точной функцией Лагранжа для жидкости показана в I02, где получена точная функция Лагранжа для системы частиц, а затем аппроксимированный Лагранжиан:
(23)
который имеет 2-ой порядок точности. Новая функция Лагранжа очень похожа на ту, которая используется в стандартном SPH-методе, единственное отличие заключается в члене, обозначающем удельную внутреннюю энергию. Удельная внутренняя энергия выглядит как сглаженная еще раз по сравнению со стандартным SPH-методом, но данная форма в виде второго члена функции Лагранжа точно такая же, как и соответствующий член в функции Лагранжа для реальной жидкости (см. уравнения (29) и (41) в I02). Уравнение импульса, полученное из уравнения (23), абсолютно такое же, как и уравнение, полученное при свертке ядра.
Для интегрирования уравнения (17) необходимо знать функциональные формы плотности и давления. Линейная или кубическая сплайн-интерполяция была использована в I02 в качестве функции плотности около частиц i и j, но требуется дальнейшее уточнение для более точного учета поля плотности. Для определения давления и скорости между частицами i и j использовалось решение задачи Римана (RPS или РЗР). Поэтому данный метод именуется «метод SPH Годунова». Как и в обычном сеточном методе Годунова, никакой явной диссипации (например, искусственной вязкости) не требуется благодаря применению решения задачи Римана (РЗР).
Заметили, что использование решения задачи Римана в GSPH напрямую не связано с отсутствием НКГ или с проблемой согласованности. Появление нефизической силы вследствие несогласованности уравнения импульса в стандартном SPH-методе было устранено в новом уравнении импульса в GSPH, полученном в результате свертки ядра или из новой функции Лагранжа. Решение задачи Римана используется для описания ударных волн, поскольку создает небольшое, но достаточное рассеяние в окрестности ударных волн. Чтобы проверить это, мы провели моделирование НКГ с помощью простейшего GSPH, предложенного Cha и Whitworth (2003a). В простейшем GSPH используется то же уравнение импульса, что и в стандартном SPH-методе, однако вместо искусственной вязкости применяется решение задачи Римана. Простейший GSPH также показывает отсутствие НКГ по градиенту плотности.
ТЕСТИРОВАНИЕ
Было проведено два вида тестирования. Первый – это традиционное моделирование НКГ в двух слоях со сдвигом скоростей, а второй – тест с каплей. Все тесты проводились при помощи двухмерного кода GSPH второго порядка[1], в который заведено адиабатическое уравнение состояния. Отношение удельных теплоемкостей, γ, равно 5/3 во всех моделях.
4.1 НКГ в двух слоях (ρu: ρl = 1: 2)
Имеется два слоя различной плотности при изначально равновесном давлении. Равновесное давление устанавливается равным 2.5, а отношение плотностей между верхним и нижним слоями полагается равным 1:2. Эти два слоя движутся в противоположных по отношению друг к другу направлениях с числами Маха в верхнем и нижнем слоях 0.22 и 0.3 соответственно. Расчетная область - [0, 1/3] x [- 1/6, 1/6]. Размера расчетной области меньше, чем у A07, с целью экономии времени вычисления. По осям х и y были заданы, соответственно, периодические и зеркальные граничные условия. Общее количество частиц внутри расчетной области составляет ~ 105, а первоначальная конфигурация распределения частиц – решетка (A07).
Первоначальное возмущение скорости в направлении у определяется как
(24)
где Ao – амплитуда возмущения, которая устанавливается как 1/40 от первоначального сдвига скоростей. Здесь λ – это длина волны первоначального возмущения, которая устанавливается равной 1/6. Таким образом, в области вычисления ожидаются два вихря. Первоначальное возмущение задается только в тонком слое (|y| < 0.05) около первоначального контактного разрыва.
При таком первоначальном сдвиге скоростей и перепаде плотностей шкала времени для НКГ определяется выражением
(25)
Здесь ρu и ρl – плотности верхнего и нижнего слоев соответственно, а 
– разница скоростей между двумя слоями, τкг в условных единицах есть 0.43.
На рисунке 3 представлены снимки, сделанные в различные моменты развития, t = 0.5, 1.0, 1.5 и 2.0 τкг. При t = 0.5τкг начальный контактный разрыв колеблется вследствие первоначального возмущения и сдвига скоростей. В последующих снимках видны хорошо закрученные вихри, развивающиеся вокруг разрыва. На снимке, сделанном в момент времени t = 2.0τкг, видно искажение вихрей, и вокруг первоначального контактного разрыва ожидается образование смешанного слоя перемешивания. В конечном итоге слой перемешивания остановит развитие НКГ. В отличие от стандартного SPH-метода, GSPH намного меньше зависит от нефизической силы по градиенту плотности, поэтому с его помощью можно описать НКГ в слоях различной плотности. Заметим, что в данном моделировании нет никакой дополнительной явной диссипации, такой как искусственная вязкость (или искусственная теплопроводность).
На рисунке 4 показано распределение давления при t = 1.0 и 2.0 τкг, На нем видно, что помех значительно меньше, чем в расчетах по стандартному SPH-методу, в котором наблюдаются всплески давления поперек контактного разрыва (см. рис. 6 в Price (2008)). Заметим, что Price (2008) получил аналогичную карту распределения давления также и с использованием искусственной теплопроводности.
4.2 НКГ в двух слоях ((ρu: ρl = 1: 10)
Моделирование НКГ, представленное в предыдущем разделе, было проведено еще, на сей раз с другим перепадом плотностей. Перепад плотностей намного больше того, который использовался в предыдущем моделировании, и составляет 1:10. Начальные числа Маха в верхнем и нижнем слоях равны 0.2 и 0.63 соответственно. Общее количество частиц, используемых в расчете, составляет ~ 105. Первоначальное возмущение такое же, как и в предыдущем расчете.
На рисунке 5 представлены результаты. Два снимка сделаны при t = 1.0 и 1.25τкг соответственно. На более ранней стадии (до 1.0τкг) ситуация очень сходна со случаем меньшего перепада плотностей, который рассматривался в предыдущем разделе. Однако на более поздней стадии завихрения не закручены, а растянуты.
Причина такого растяжения не понятна, но мы полагаем, что это может происходить вследствие худшего, по сравнению с предыдущим моделированием, разрешением верхнего слоя (Price 2008). GSPH (как и стандартный SPH-метод) является лагранжевым методом, поэтому численное разрешение зависит от численной плотности частиц. При аналогичном общем количестве частиц более высокий перепад плотностей между двумя слоями приводит в конечном итоге к худшему разрешению слоя с более низкой плотностью. Другой возможной причиной растяжения вихрей является начальное давление. В данном моделировании в качестве равновесного давления мы использовали значение 2.5, но если выбрать иное значение, то результат может измениться. Тем не менее, мы хотим особо отметить, что НКГ может иметь место в случае большой разницы плотностей.
