CC BY
0
2023 ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА Т. 10 (68). Вып. 4
МАТЕМАТИКА. МЕХАНИКА. АСТРОНОМИЯ
МЕХАНИКА
УДК 533 MSC 76-10
Применение гибридных римановских солверов на основе HLLC и HLL для моделирования течений с газодинамическими разрывами*
Г. В. Шоев
Институт теоретической и прикладной механики им. С. А. Христиановича СО РАН, Российская Федерация, 630090, Новосибирск, ул. Институтская, 4/1
Для цитирования: Шоев Г. В. Применение гибридных римановских солверов на основе HLLC и HLL для моделирования течений с газодинамическими разрывами // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2023. Т. 10 (68). Вып. 4. С. 775-791. https://doi.org/10.21638/spbu01.2023.413
Обсуждается применение гибридных приближенных римановских солверов на основе стандартных HLLC- и HLL-солверов. Рассмотрено три различных гибридных сол-вера. Первый гибридный солвер (rHLLC-HLL) использует взвешенную сумму HLLC и HLL так, что HLLC используется по нормали к ударной волне, а HLL вдоль нее. Второй гибридный солвер (HLLC-ADC) использует взвешенную сумму HLLC и HLL, применяя в качестве весов функцию от давления в центрах ячеек слева и справа от соседних граней. Третий гибридный солвер (HLLC-HLL) выполняет расчет невязких потоков по HLL внутри ударных волн и по HLLC в других областях течения. Грани внутри ударных волн определяются индикатором ударной волны, основанном на реконструированных значениях давления слева и справа от грани. Проведен ряд тестов, и показано, что гибридные солверы позволяют избавиться от карбункула, уменьшают осцилляции на ударных волнах. Точное разрешение контактных разрывов показали гибридные солверы rHLLC-HLL и HLLC-HLL.
Ключевые слова: конечно-объемный метод, приближенный римановский солвер, карбункул, римановский солвер HLLC, римановский солвер HLL.
1. Введение. Точное численное моделирование течений с газодинамическими разрывами является сложной и нерешенной задачей. Одной из проблем вычислительной аэродинамики является численный артефакт карбункул (carbuncle), который встречается в различных течениях, например на ударных волнах около затупленных тел [1, 2], или может косвенно проявляться при истечении из сопел [3]. При
* Работа выполнена в рамках государственного задания ИТПМ СО РАН. © Санкт-Петербургский государственный университет, 2023
обтекании затупленных тел карбункул характеризуется тем, что у ударной волны появляется излом вверх по потоку. Для преодоления карбункула разработаны различные методы. В частности, при обтекании осесимметричного затупленного тела в [4] вводилась дополнительная искусственная вязкость, значение которой увеличивается к оси. В результате, удается значительно улучшить качество численного решения. Недостатком такого подхода является то, что он применим только к осесимметричным течениям и, как отмечается в [4], не является решением проблемы карбункула. Применение искусственной вязкости получило развитие в [5-7], где предлагается более общий способ использования подавления карбункула, а именно использования искусственной вязкости в отдельных областях течений.
Другим способом является использование метода расчета невязких потоков, не приводящего к эффекту карбункула путем тестирования [8, 9] с последующим отбором наиболее подходящих солверов или применения солвера с большой численной диссипацией [10]. В последнем случае численное решение не содержит осцилляций на ударных волнах, но в областях без них также применяется солвер с большой численной диссипацией (например, на контактных разрывах). Для уменьшения численной диссипации и отсутствия карбункула необходима модификация римановского сол-вера (например, [11-13]). Некоторой сложностью здесь является огромный выбор солверов со своими достоинствами и недостатками (например, семейство AUSM-солверов). Для существующих газодинамических решателей дополнительный римановский солвер — это дополнительный код, который необходимо разработать, внедрить, верифицировать и поддерживать. В то же время существующие газодинамические решатели, как правило, имеют набор стандартных (фактически классических) солверов. Поэтому возникает естественная задача преодоления карбункула и других численных артефактов с использованием римановских солверов, ставших фактически классическими для современной вычислительной аэродинамики.
Целью настоящей статьи является реализация нескольких гибридных римановских солверов на основе солверов HLL [14] и HLLC [15, 16] для преодоления карбункула и уменьшения осцилляций, исходящих от ударных волн, а также сравнение и анализ результатов расчетов с использованием реализованных солверов в ряде стандартных случаев. Для реконструкции переменных слева и справа от грани используется стандартная MUSCL-реконструкция. Такая комбинация реконструкции и римановских солверов проста в реализации по сравнению с более современными методами [17-19]. В частности, в [19] используется схема ADER высокого порядка, требующая относительно сложной реализации: разложение в ряд Тэйлора по времени на грани с переходом от временных производных к функционалам пространственных производных и решением задачи Римана. В ADER необходима реконструкция высокого порядка, как правило, WENO, а не MUSCL. В качестве расчетного кода для реализации используется код HyCFS [20], разрабатываемый в ИТПМ СО РАН, и MUSCL в качестве базовой реконструкции.
Необходимо отдельно отметить, что проблема карбункула не всегда рассматривается как исключительно численный артефакт (см. обсуждения в [8, 21, 22]). Природа естественных возмущений на ударных волнах и их корректное моделирование требуют отдельного исследования, поэтому в настоящей работе будут рассмотрены только те случаи, где карбункул считается численным артефактом (т. е. не наблюдается в экспериментах).
Далее будут представлены основные уравнения, решаемые в настоящей работе, приближенные римановские солверы HLL и HLLC, реализованные в HyCFS, и их
гибриды, построенные на основе различных литературных источников. Затем показаны тестовые случаи моделирования течений с этими римановскими солверами. Основные выводы сформулированы в заключении.
2. Основные уравнения. В настоящей работе течение газа рассчитывается на основе уравнений Эйлера:
дQ | dFz | д¥х
dt dz dx
0,
(1)
где консервативные переменные Р и невязкие потоки ¥;;, Еж записываются как:
Q
( pw\ pu pE
\PJ
Fz
/ pww + p \ pwu
P[E +
V
pw
F F x
puw
puu + p p
\
P[E +
p
pu
/
|u|2
Здесь p — плотность; u = (w, и) — вектор скорости; E = CVT -\—---внутренняя
энергия, Cv — теплоемкость при постоянном объеме, T — температура, p — давление.
Для замыкания системы уравнений (1) используется уравнение состояния идеального газа:
П,
р = fJwT'
(2)
где ^ — универсальная газовая постоянная; Ш — молекулярный вес.
Уравнения (1) решаются численно с использованием конечно-объемного метода на структурированной сетке, состоящей из четырехугольных ячеек. Газодинамические параметры хранятся в центрах ячеек. Для численного решения уравнения (1) переписываются в обобщенных криволинейных координатах. Рис. 1 иллюстрирует преобразования от (х,у) координат в физической плоскости к координатам (£,£) вычислительной плоскости. Соответствующие уравнения переписываются следующим
w
образом:
dt
дС
или в полудискретной форме
д_
dt
(V Q) = -
|SC |F Л -(\SZ |F z ' ' Vk+1/2,i y Vfc-1/2,i
|sç|FЛ - (S|Fç)
J k,i+1/2 V J k,i+1/2
где использовано преобразование к обобщенным криволинейным координатам x = х(С), С = (С,С); J — якобиан преобразования; V — объем ячейки; k, i — индексы ячеек вдоль направлений Z и С- Целые индексы к и i соответствуют номерам ячеек, а полуцелые — номерам их граней. Невязкие потоки Fz и Fç связаны с потоками Fz и Fx следующими соотношениями:
iZ Fz + nXFx
Tz Fz + nXFx
(4)
где П = (п^, пХ) и п£ = (п|, пХ) — единичные нормали граней ячеек. В итоге невязкие потоки F^ и F£ имеют вид
( pwwn + pnZ \ puwn + pnX
р ( Е + - | Wr, PWn
\
( pwun + pnZ \ puun +pnX
p[E+P-P
PUn
\
(5)
/
где тп = тп^ + ипХ и ип = тп^ + ипХ.
Компоненты квазивектора Р определены в центрах ячеек, а невязкие потоки F^ и F£ должны быть вычислены на их гранях. Для вычисления невязких потоков F^ и F£ в настоящей работе рассмотрены приближенные римановские солве-ры ИЬЬС, ИЬЬ и их гибриды гИЬЬС-ИЬЬ, HLLC-ADC и ИЬЬС-ИЬЬ. Для использования римановского солвера необходимо определить состояния слева и справа от грани ячейки. Состояния слева и справа от грани вычисляются стандартной MUSCL-реконструкцией [23] первого порядка при моделировании течения около цилиндра и третьего порядка для всех остальных задач.
3. Приближенный римановский солвер HLL. Невязкие потоки F^ и F£ вычисляются по формуле
™ HLL
F j + 1/2
bpF l — bmF r + bmbp(Qn — Ql)
bp — bn
где индекс ] означает индекс ячейки вдоль направления £ (у = к) или £ (у = г); индексы Ь и К соответствуют параметрам слева и справа от грани ] + 1/2, а Ьр, Ьт, вычисляются по формулам
bm = min(bL, 0),
bL = min(vL — cL, vR — cr),
bp = max(bR, 0),
/ L R \
max( Vn + cL,vn + cr),
(6)
д
z
ç
z
ç
u
a
где с — скорость звука; Уп — скорость по нормали к грани 2 + 1/2 (уп = ™п, если 2 = к или уп = ип, если 2 = г).
4. Приближенный римановский солвер HLLC. Невязкие потоки Е^ и Е^ вычисляются по формуле
Е нььс
Е 3+1/2
Е ь, Г в,
г (и*- ь), Е (и* в),
если Ьь > 0, если Ья < 0, если Ьт > 0, иначе,
где квазивекторы невязких потоков Е(и*ь) и Е(И*я) определяются как:
Г (и *
ь) =
(Р™)Ь Ьт (ри)т Ьт
(Ет + Рт) Ьт
рт Ьт
Е (и *
я) =
(Р™)Гт Ьт (ри)т Ьт (Ет + Рт) Ьт пт Ь
РЯ Ьт
Здесь используются следующие обозначения:
(Р™)т = Рт™Ь + Пх (рт - Рь)оь (ри)т = ртиь + пх(рт - Рь)оь Ет = (рьЕь(Ьь-У^) -РЬУЬ + РтЬт)оь
т
Рь = аь оь аь = Рь (Ьь - ) 1
оь =
Ь - Ьт
(7)
(р™)т = р^я + Пх (Рт - Ря)оя, (ри)т = Ртив + Пх(Рт - Ря)оя,
Ет = (РЯ Ея(Ья - - РЯУЩ + РтЬт)ов,
Рв = ая оя, ая = Ря(Ья - уя)
оя =
Ья Ьт
где Уп определяется так же, как для ИЬЬ, Ьь и Ья — по выражениям (6), а Ьт и Рт выражаются как:
Ьт =
я
авУп ~ аьУ„ + Рь ~ Ря ая - а
Рт = Р ь +
а
Ьт
5. Гибридный солвер rHLLC-HLL с учетом направления ударной волны. Основная идея этого варианта гибридного солвера описана в [13, 24]. В общем виде невязкий поток через грань записывается следующим образом:
Е гИЬЬС-ИЬЬ _ ^ Е ИЬЬ / N , ^ Е N
Е 3+1/2 = а1 Е 3+1/2 (П1) + а2 Е 3+1/2 (п2 )
(8)
где коэффициенты «1 = |п • П1 \ и «2 = |п • П2|. Здесь п — вектор единичной нормали к грани ячейки. Вектор П1 в общем случае может определяться по-разному. В настоящей работе используется одно из определений [13]:
п1
Ад |Д<1Г
если |Д<| > е,
п+, иначе,
(9)
1
где q = (wr —wl, «д — ul)', | Aq| = \/(wr — wl)2 + (ur — ul)2; e — некоторое малое число, а n^ — единичный вектор, перпендикулярный к вектору n. Вектор П2 L ni определяется из уравнения ni • П2 = 0. Оба вектора ni и П2 являются единичными и направлены в одно полупространство с вектором n (рис. 2). Все рассмотренные векторы связаны соотношением: n = aini + а2n2. С учетом определения (9) выражение для невязкого потока (8) сводится к:
( ai FhlL/2(ni)+ a2 FH+L/^n), если |Aq| > e,
™ rHLLC-HLL I /1П\
F j+i/2 = ) Л (10) I F"HL-L/Cj (n), иначе.
Для трехмерного случая вектор n2 должен определяться по формуле из [13, 25]. Солвер rHLLC-HLL устроен так, что менее диссипативный HLLC использует направление по нормали к ударной волне, а более диссипативный HLL — направление вдоль нее. Детальное обсуждение этого вопроса дано в [13].
6. Гибридный солвер HLLC-ADC. Этот гибридный солвер, фактически предложенный в работе [12], рассчитывает невязкие потоки через грань ячейки следующим образом:
F HLLC-ADC = (1 ш) FHLL + ш FHLLC (11)
Ff =j + i/2 = (1 " ш) Ff =j+i/2 + ш Ff =j+i/2, (11)
где ш в двумерном случае вычисляется согласно
(PCR Pc'l \ а
ш = min (/г), /г = min —, — , / = /ь ... /4.
\PCL PCR)l
Здесь обозначения граней показаны на рис. 3, а pCL и pCR — давление в центре ячейки слева и справа соответственно. Константа a фактически регулирует вклад численной диссипации. В расчетах использовалось a = 3 [12].
Рис. 3. Шаблон для индикатора ударной волны в 2D-расчете.
7. Приближенный гибридный римановский солвер HLLC-HLL. Гибридный солвер объединяет реализации солверов HLL и HLLC. В ударной волне солвер использует HLL с большей численной вязкостью, чтобы избежать осцилляций.
fl Cell L fs f Cell R
f f
Вне ударной волны солвер переключается на менее диссипативный ИЬЬС, чтобы точнее рассчитывать потоки в областях без ударных волн (в частности, контактные разрывы). Для определения, лежит ли грань ] внутри ударной волны, используется индикатор ударной волны [11, 26], который работает следующим образом. Если на грани f или хотя бы на одной из соседних граней выполняется условие
\рп~рь\ Р_П±РЬ^
тт(рь,ря) 2
то считается, что грань находится внутри ударной волны. Здесь £р — некоторый управляющий параметр. В настоящей работе параметр ер подбирается индивидуально для каждого расчета так, чтобы неравенство (12) выполнялось только на ударных волнах, поскольку универсальный механизм определения параметра £р в настоящий момент автором не найден.
Индикатор ударной волны в двухмерном случае (рис. 3) определяется как:
SDf =max(sdf, , , , ) , (13)
где
1 если Ьд ~Рь\ > £ Рп + Рь М = { ' 1Ш11(рЬ,рп) р 2 ' , I = /, /ь /2, /з, /4-
0, иначе.
Необходимо заметить, что стандартные солверы ИЬЬ и ИЬЬС фактически выполняют вычисление потоков, не используя информацию о параметрах на соседних гранях (например, вычисление потоков вдоль направления £ прямо не зависит от параметров вдоль направления £, и наоборот). В этом смысле стандартные солверы ИЬЬ и ИЬЬС вычисляют потоки как в одномерном случае, независимо по каждому направлению. Гибридный солвер ИЬЬС-ИЬЬ учитывает информацию о давлении на соседних гранях, т. е. способ вычисления потоков вдоль направления £ зависит от распределения давления вдоль направления £. В этом заключается еще одно принципиальное различие стандартного и гибридного солверов.
Конечное выражение для невязких потоков:
Е ИЬЬС-ИЬЬ
* !=3+1/2
Г И=Ь+1/2, если SDf=3+1/2 = 1
Г И=Ь+С1/2, иначе.
8. Результаты численных расчетов. 8.1. Обтекание цилиндра. Одним из классических примеров течения, где встречается численный артефакт карбункула, является течение около цилиндра (рис. 4). Здесь и далее рассматривается газ с показателем адиабаты 7 = 1.4. На левой (входной) границе задается равномерный набегающий поток с = 6044 м/с, Тж = 219 К, рж = 3, 549 Па, что соответствует числу Маха Мж = 20. На цилиндре используется граничное условие непротекания (невязкая стенка). На выходных границах задается сверхзвуковой выход (снос всех параметров из расчетной области). В настоящей работе расчеты проводились на двух сетках: 60x281 и 180x281. Вдоль цилиндра выбрано нечетное число ячеек, чтобы центры среднего ряда ячеек лежали на линии торможения. На рис. 4
Рис.4. Изотермы около цилиндра в расчетах с HLL-солвером на сетках: о — 60x281, д — 180x281; с HLLC-солвером на сетках: Ь — 60x281, Ь — 180x281; с гHLLC-HLL-солвером на сетках: с — 60x281 и г — 180x281; с HLLC-ADC-солвером на сетках: ё — 60x281 и ] — 180x281; с HLLC-HLL-солвером на сетках: е — 60x281, к — 180x281, и f, I — индикатор ударных волн согласно (13) в расчетах с использованием HLLC-HLL-солвера.
видно, что HLL обеспечивает численное решение без осцилляций, исходящих от головного скачка уплотнения. Солвер HLLC дает карбункул на сетке 60x281, но при увеличении числа ячеек по нормали к поверхности цилиндра головной скачок приобретает корректную форму. Тем не менее на сетке 180x281 HLLC дает численный артефакт в поведении изотерм между скачком и поверхностью цилиндра вдоль линии торможения. Отметим, что все расчеты проведены с реконструкцией первого порядка, это на практике улучшает устойчивость численного решения. Гибридный солвер HLLC-HLL позволяет получить численное решение без карбункула и других численных артефактов на обеих сетках. В расчетах с солвером HLLC-HLL использовалось ep = 0.01, что позволило точно определить место положения скачка (рис. 4, f и l). На гранях в области SD использовался HLL, а на гранях вне этой области — HLLC. Гибридные варианты rHLLC-HLL и HLLC-ADC также позволяют избавиться от карбункула.
8.2. Двойное маховское отражение. Двойное маховское отражение также используется для тестирования римановских солверов (детальная классификация
г, м г, м
Рис. 5. Изотермы двойного маховского отражения в расчетах с использованием HLL (а), HLLC (Ь), гHLLC-HLL (с), HLLC-ADC (ё), HLLC-HLL (е) и индикатор ударных волн (/) согласно (13) с HLLC-HLL.
отражения ударных волн дана, например, в [27, 28]). В настоящей работе рассматривается отражение ударной волны с числом Маха 10, набегающей на клин с углом полураствора 30°. В расчетах использовались следующие параметры до (индекс 1) и после (индекс 2) ударной волны: и>1 =0, Т = 1 К, р1 = 1 Па, и>2 = 168, 2 м/с, Т2 = 20, 38 К, р2 = 116, 5 Па. Расчетная область представлена на рис. 5. На левой
границе задаются параметры, рассчитанные по соотношениям Ренкина — Гюгонио за нормальной ударной волной. На остальных границах условие непротекания (невязкая стенка). В начальный момент времени ударная волна расположена при z = 0. Расчет продолжается до времени t = 0.01 с. На рис. 5 представлены результаты расчетов. Солвер HLL позволяет получить численное решение с небольшим количеством осцилляций на ударных волнах. Солвер HLLC производит большее количество осцилляций на ударных волнах, а в поле течения наблюдаются небольшие предвестники карбункула на ножке Маха, отраженной от клина, и на падающей ударной волне на верхней границе. Солверы rHLLC-HLL и HLLC-ADC позволяют получить численное решение без карбункула. В то же время видно, что солвер rHLLC-HLL приводит к большему количеству осцилляций на ножке Маха, чем другие гибридные солверы. Гибридный солвер HLLC-HLL позволяет получить численное решение без больших осцилляций на ударных волнах. В расчетах использовалось ep = 0.0001. Индикатор ударной волны также показал хорошие результаты по определению положения ударных волн. Небольшая ошибка заметна на носике клина, где, согласно выражению (13), SDf = 1. Однако в области контактной поверхности индикатор ударной волны сработал безошибочно.
8.3. Дифракция ударной волны на угле 90°. В начальный момент времени в координате z = 0 задается ударная волна с числом Маха 5.09, которая далее распространяется за угол 90°. В расчетах использовались следующие параметры до (индекс 1) и после (индекс 2) ударной волны: wi =0, Ti = 1 K, pi = 1 Па, w2 = 83,1 м/с, T2 = 5, 98 K, p2 = 30,06 Па. Расчет продолжается до 0.0105 с. На левой границе задаются значения в соответствии с соотношениями Ренкина — Гюгонио, а на всех остальных границах используется граничное условие невязкой стенки. На рис. 6 показана расчетная область и изолинии плотности. Видно, что HLLC-солвер дает небольшие осцилляции на ударной волне при взаимодействии с верхней границей. Солверы HLL, rHLLC-HLL, HLLC-ADC и HLLC-HLL дают численное решение без сильных осцилляций на ударной волне; в HLLC-HLL-расчетах использовалось ep = 0.1, что привело к удовлетворительной работе индикатора ударной волны по выражению (13). Ложное срабатывание индикатора ударной волны произошло на угле, где формируется центрированная волна разрежения Прандт-ля — Майера. На самой ударной волне вблизи точки (0.1, 0.11) видно, что индикатор не сработал, но на поле плотности не наблюдается значительного увеличения осцилляций. Важно отметить, что индикатор ударной волны корректно сработал на ударной волне ниже по потоку за основной падающей ударной волной.
8.4. Контактный разрыв. Моделирование контактного разрыва производится в квадратной области 1x 1. В качестве начального условия используется разрывное решение в верхней и нижней половине расчетной области. В верхней половине задается число Маха 2, а в нижней — 1.1. В нижней половине плотность в 10 раз больше, чем в верхней. Давление одинаковое во всей расчетной области. В расчетах использовались следующие параметры (сверху (индекс t) и снизу (индекс b)) ударной волны: wt = 705, 5 м/с, Tt = 299, 3 K, pt = 101325 Па, wb = 122, 7 м/с, Ть = 29, 9 K, pb = 101325 Па. Аналогичные входные условия задаются на верхней и нижней половинах левой границы расчетной области. На правой границе сверхзвуковой выход, а на верхней и нижней границах условие непротекания. На рис. 7 показаны изолинии плотности, полученные в расчетах. Видно, что HLL- и HLLC-ADC-солверы плохо воспроизводят скачок плотности: контактный разрыв сильно
Рис. 6. Дифракция ударной волны на угле 90°. Изолинии плотности в расчетах с использованием HLL (а), HLLC (Ь), гHLLC-HLL (с), HLLC-ADC (ё), HLLC-HLL (е) и индикатора ударных волн (ё) согласно (13) с HLLC-HLL.
размазывается, и поле плотности в нижней половине расчетной области заметно неравномерное. Результаты расчетов с использованием солверов ИЬЬС- и ИЬЬС-ИЬЬ полностью совпадают и дают значительно меньшее размазывание контактного разрыва, сохраняя равномерность потока в обеих половинах расчетной области.
Рис. 7. Изолинии плотности на контактном разрыве в расчетах с использованием HLL (а), HLLC (b), rHLLC-HLL (с), HLLC-ADC (d), HLLC-HLL (e) и индикатора ударных волн (d) согласно выражению (13) с HLLC-HLL.
Также необходимо отметить, что индикатор ударной волны не работал, т. е. И^С-Н^ полностью перешел в И^С. В расчетах И^С-И^ использовался ер = 0.1. Гибридный вариант rHLLC-HLL хорошо сохраняет контактный разрыв.
9. Заключение. Проведено тестирование римановских солверов HLL, HLLC, rHLLC-HLL, HLLC-ADC и HLLC-HLL, реализованных в программном пакете HyCFS. Показано, что при численном решении уравнений Эйлера солверы HLL и HLLC-ADC позволяют получить численное решение без существенных осцилляций на ударных волнах, но плохо описывают контактные разрывы. Использование сол-вера HLLC может приводить к возникновениям осцилляций на ударных волнах и численного артефакта карбункула, но хорошо описывает контактный разрыв. Гибридный солвер HLLC-HLL объединяет достоинства HLL и HLLC: внутри ударных волн производится солвером HLL, а вне ударных волн — солвером HLLC. Это позволяет получать численные решения без значительных осцилляций на ударных волнах и точно описывать области течения с контактными разрывами. К недостаткам HLLC-HLL относится необходимость точного определения положения ударных волн, которое выполняется через использование индикатора ударных волн с одним управляющим параметром. Гибридный солвер rHLLC-HLL также допускает получить численное решение без карбункула и точно описывает контактные разрывы. Однако в случае двойного маховского отражения rHLLC-HLL показал большее количество осцилляций на ножке Маха.
Благодарности. Автор выражает благодарность Т. Ю. Шкредову за идею статьи, проф. А. Н. Кудрявцеву за плодотворную дискуссию и канд. физ.-мат. наук А. А. Шершневу за техническую поддержку при реализации гибридных солверов в HyCFS. Численные расчеты выполнялись с использованием ресурсов вычислительного центра Новосибирского государственного университета (НГУ) и ИТПМ СО РАН.
Литература
1. Туник Ю.В. Проблемы численного моделирования на основе некоторых модификаций схемы Годунова. Физико-химическая кинетика в газовой динамике 19 (1), 1—11 (2018). http://doi.Org/10.33257/PhChGD.19.1.701
2. Tunik Y. V. Problems of Numerical Simulation Based on Some Modifications of Godunov's Scheme. Fluid Dyn 57, S75-S83 (2022). https://doi.org/10.1134/S0015462822601334
3. Карпов А. В., Васильев Е. А. Численное моделирование истечения перерасширенной струи газа из короткого осесимметричного сопла. Вестник ВолГУ. Серия 1. Прикладная математика 9, 81-88 (2005).
4. MacCormack R. The Carbuncle CFD Problem. 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition (2011). https://doi.org/10.2514/6.2011-381
5. Rodionov A.V. Artificial viscosity to cure the carbuncle phenomenon: The three-dimensional case. Journal of Computational Physics 361, 50-55 (2018). https://doi.org/10.1016 /j.compfluid.2019.06.011
6. Родионов А. В. Искусственная вязкость для подавления численной неустойчивости типа «карбункул» в расчетах трехмерных задач. ВАНТ. Сер. Математическое моделирование физических процессов 3, 44-51 (2018).
7. Тагирова И. Ю., Родионов А. В. Применение искусственной вязкости для борьбы с «кар-бункул»-неустойчивостью в схемах типа Годунова. Математическое моделирование 27 (10), 47-64 (2015).
8. Исмагилов Д. Р., Сидельников Р. В. Особенности численного моделирования гиперзвукового обтекания простых тел. Вестник концерна ПВО «Алмаз — Антей» 2, 49-54 (2015). https://doi.org/10.38013/2542-0542-2015-2-49-54
9. Смирнова Н. С. Сравнение схем с расщеплением потока для численного решения уравнений Эйлера сжимаемого газа. Труды МФТИ. Математика 10 (1), 122-141 (2018).
10. Колесник Е. В., Смирнов Е. М., Смирновский А. А. Численное решение трехмерной задачи обтекания установленного на пластине цилиндрического тела сверхзвуковым потоком вязкого газа при M = 2,95. Математическое моделирование физических процессов 12 (2), 7-22 (2019). https://doi.org/10.18721/JPM.12201
11. Shen Z., Yan W., Yuan G. A robust HLLC-type Riemann solver for strong shock. Journal of Computational Physics 309, 185-206 (2016). https://doi.org/10.1016/j.jcp.2016.01.001
12. Simon S., Mandal J.C. A cure for numerical shock instability in HLLC Riemann solver using antidiffusion control. Computers and Fluids 174, 144-166 (2018). https://doi.org/10.1016 /j.compfluid.2018.07.001
13. Nishikawa H., Kitamura K. Very simple, carbuncle-free, boundary-layer-resolving, rotated-hybrid Riemann solvers. Journal of Computational Physics 227 (4), 2560-2581 (2008). https://doi.org/10.1016/j.jcp.2007.11.003
14. Harten A. High resolution schemes for hyperbolic conservation laws. Journal of Computational Physics 49 (3), 357-393 (1983). https://doi.org/10.1016/0021-9991(83)90136-5
15. Toro E. F., Spruce M., Speares W. Restoration of the contact surface in the HLL-Riemann solver. Shock Waves 4 (1), 25-34 (1994). https://doi.org/10.1007/BF01414629
16. Batten P., Leschziner M.A., Goldberg U. C. Average-State Jacobians and Implicit Methods for Compressible Viscous and Turbulent Flows. Journal of Computational Physics 137 (1), 38-78 (1997). https://doi.org/10.1006/jcph.1997.5793
17. Toro E. F., Titarev V. A. MUSTA fluxes for systems of conservation laws. Journal of Computational Physics 216, 403-429 (2006). https://doi.org/10.1016/j.jcp.2005.12.012
18. Titarev V.A., Romenski E., Toro E. F. MUSTA-type upwind fluxes for non-linear elasticity. International Journal for Numerical Methods in Engineering 73 (7), 897-926 (2008). https://doi.org/10.1002/nme.2096
19. Dematte R., Titarev V. A., Montecinos G.I., Toro E. F. ADER-methods for hyperbolic equations with a time-reconstruction solver for the generalized Riemann problem: the scalar case. Communications on Applied Mathematics and Computation 2, 369-402 (2020). https://doi.org/10.1007/s42967-019-00040-x
20. Shershnev A. A., Kudryavtsev A. N., Kashkovsky A. V., Shoev G. V., Borisov S. P., Shkre-dov T. Yu., Polevshchikov D.P., Korolev A. A., Khotyanovsky D. V., Kratova Yu. V. A Numerical Code for a Wide Range of Compressible Flows on Hybrid Computational Architectures. Supercomputing Frontiers and Innovations 9 (4), 85-99 (2022). https://doi.org/10.14529/jsfi220408
21. Тарнавский Г. А., Алиев А. В. Особенности аэродинамики высокоскоростного полета: компьютерное моделирование гиперзвукового обтекания головной части объекта. Вычислительные методы и программирование 9 (4), 371-394 (2008).
22. Родионов А. В. Разработка методов и программ для численного моделирования неравновесных сверхзвуковых течений в приложении к аэрокосмическим и астрофизическим задачам: дис. ... д-ра физ.-мат. наук. Саратов, Институт прикладной математики им. М. В. Келдыша РАН. (2019).
23. Leer B. van. Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method. Journal of Computational Physics 32 (1), 101-136 (1979). https://doi.org/10.1016/0021-9991(79)90145-1
24. Levy D. W., Powell K. G., van Leer B. Use of a rotated Riemann solver for the two-dimensional Euler equations. Journal of Computational Physics 106 (2), 201-214 (1993). https://doi.org/10.1016/S0021-9991(83)71103-4
25. Ren Yu-X. A robust shock-capturing scheme based on rotated Riemann solvers. Computers and Fluids 32 (10), 1379-1403 (2003). https://doi.org/10.1016/S0045-7930(02)00114-7
26. Sanders R., Morano E., Druguet M.-C. Multidimensional Dissipation for Upwind Schemes: Stability and Applications to Gas Dynamics. Journal of Computational Physics 145 (2), 511-537 (1998). https://doi.org/10.1006/jcph.1998.6047
27. Семенов А. Н., Березкина М. К., Красовская И. В. Классификация разновидностей отражения ударной волны от клина. Ч. 1. Границы и области существования различных типов отражения. Журнал технической физики 79 (4), 46-51 (2009).
28. Семенов А. Н., Березкина М. К., Красовская И. В. Классификация разновидностей отражения ударной волны от клина. Ч. 2. Экспериментальное и численное исследование разновидностей маховского отражения. Журнал технической физики 79 (4), 52-58 (2009).
Статья поступила в редакцию 7 марта 2023 г.;
доработана 2 мая 2023 г.; рекомендована к печати 18 мая 2023 г.
Контактная информация:
Шоев Георгий Валерьевич — канд. физ.-мат. наук, ст. науч. сотр.; shoev@itam.nsc.ru
Application of hybrid Riemann solvers based on HLLC and HLL for simulation of flows with gas-dynamic discontinuities*
G. V. Shoev
Khristianovich Institute of Theoretical and Applied Mechanics of the Siberian Branch of the Russian Academy of Sciences, 4/1, ul. Institutskaya, Novosibirsk, 630090, Russian Federation
For citation: Shoev G. V. Application of hybrid Riemann solvers based on HLLC and HLL for simulation of flows with gas-dynamic discontinuities. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2023, vol. 10(68), issue 4, pp. 775-791. https://doi.org/10.21638/spbu01.2023.413 (In Russian)
Application of hybrid approximate Riemannian solvers based on standard HLLC and HLL solvers is discussed. Three different hybrid solvers are considered. The first hybrid solver (rHLLC-HLL) uses a weighted sum of HLLC and HLL so that HLLC is applied in the direction normal to the shock wave while HLL is applied in the direction along it. The second hybrid solver (HLLC-ADC) uses the weighted sum of HLLC and HLL, applying as weights the pressure function at the left and right cell centers. The third hybrid solver (HLLC-HLL) computes inviscid fluxes by HLL inside shock waves, and by HLLC in other areas of the flow. Facets within shock waves are defined by a shock wave indicator based on the reconstructed left and right pressure values of the face. Several tests are performed, and it is shown that hybrid solvers prevent the carbuncle and reduce oscillations on shock waves.
Keywords: finite volume method, approximate Riemannian solver, carbuncle, Riemannian solver HLLC, Riemannian solver HLL.
References
1. Tunik Y. Problems of numerical modeling on the basis of some modifications of the Godunov's scheme. Physical-Chemical Kinetics in Gas Dynamics 19 (1), 1—11 (2018). http://doi.org/10.33257/PhChGD.19.L701 (In Russian)
2. Tunik Y. V. Problems of Numerical Simulation Based on Some Modifications of Godunov's Scheme. Fluid Dyn 57, S75-S83 (2022). https://doi.org/10.1134/S0015462822601334
3. Karpov A. V., Vasiliev E.A. Numerical simulation of the outflow of an overexpanded gas jet from a short axisymmetric nozzle. Vestnik VolGU. Series 1. Applied Mathematics 9, 81-88 (2005). (In Russian)
4. MacCormack R. The Carbuncle CFD Problem. 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition (2011). https://doi.org/10.2514/6.2011-381
5. Rodionov A.V. Artificial viscosity to cure the carbuncle phenomenon: The three-dimensional case. Journal of Computational Physics 361, 50-55 (2018). https://doi.org/10.1016 /j.compfluid.2019.06.011
6. Rodionov A.V. The artificial viscosity for suppressing the «carbuncle» type numerical instability in 3d problem simulations. QAST, ser. Mathematical modeling of physical processes 3, 44-51 (2018). (In Russian)
7. Tagirova I. Yu., Rodionov A. V. Application of the artificial viscosity for suppressing the carbuncle phenomenon in Godunov-type schemes. Mathematical models and computer simulations 27 (10), 47-64 (2015). https://doi.org/10.1134/S2070048216030091 (In Russian)
8. Ismagilov D. R., Sidelnikov R. V. Features of numerical simulation of hypersonic flow around simple bodies. Journal of «Almaz — Antey» Air and Space Defence Corporation 2, 49-54 (2015). https://doi.org/10.38013/2542-0542-2015-2-49-54 (In Russian)
9. Smirnova N. S. Comparison of flux splitting schemes for numerical solution of the compressible Euler equations. Proceedings of MIPT. Mathematics 10 (1), 122-141 (2018). (In Russian)
10. Kolesnik E. V., Smirnov E. M., Smirnovsky A. A. Numerical solution of a 3D problem on a supersonic viscous gas flow past a plate-cylindrical body junction at M = 2.95. St. Pe-
*The work was carried out within the framework of the state task of Institute of Theoretical and Applied Mechanics of the Siberian Branch of the Russian Academy Sciences.
tersburg Polytechnical State University Journal. Physics and Mathematics 12 (2), 7—22 (2019). https://doi.org/10.18721/JPM.12201 (In Russian)
11. Shen Z., Yan W., Yuan G. A robust HLLC-type Riemann solver for strong shock. Journal of Computational Physics 309, 185-206 (2016). https://doi.org/10.1016/j.jcp.2016.01.001
12. Simon S., Mandal J.C. A cure for numerical shock instability in HLLC Riemann solver using antidiffusion control. Computers and Fluids 174, 144-166 (2018). https://doi.org/10.1016 /j.compfluid.2018.07.001
13. Nishikawa H., Kitamura K. Very simple, carbuncle-free, boundary-layer-resolving, rotated-hybrid Riemann solvers. Journal of Computational Physics 227 (4), 2560-2581 (2008). https://doi.org/10.1016/j.jcp.2007.11.003
14. Harten A. High resolution schemes for hyperbolic conservation laws. Journal of Computational Physics 49 (3), 357-393 (1983). https://doi.org/10.1016/0021-9991(83)90136-5
15. Toro E. F., Spruce M., Speares W. Restoration of the contact surface in the HLL-Riemann solver. Shock Waves 4 (1), 25-34 (1994). https://doi.org/10.1007/BF01414629
16. Batten P., Leschziner M.A., Goldberg U. C. Average-State Jacobians and Implicit Methods for Compressible Viscous and Turbulent Flows. Journal of Computational Physics 137 (1), 38-78 (1997). https://doi.org/10.1006/jcph.1997.5793
17. Toro E. F., Titarev V. A. MUSTA fluxes for systems of conservation laws. Journal of Computational Physics 216, 403-429 (2006). https://doi.org/10.1016/jjcp.2005.12.012
18. Titarev V.A., Romenski E., Toro E. F. MUSTA-type upwind fluxes for non-linear elasticity. International Journal for Numerical Methods in Engineering 73 (7), 897-926 (2008). https://doi.org/10.1002/nme.2096
19. Dematte R., Titarev V. A., Montecinos G.I., Toro E. F. ADER-methods for hyperbolic equations with a time-reconstruction solver for the generalized Riemann problem: the scalar case. Communications on Applied Mathematics and Computation 2, 369-402 (2020). https://doi.org/10.1007/s42967-019-00040-x
20. Shershnev A. A., Kudryavtsev A. N., Kashkovsky A. V., Shoev G. V., Borisov S. P., Shkre-dov T. Yu., Polevshchikov D.P., Korolev A. A., Khotyanovsky D. V., Kratova Yu. V. A Numerical Code for a Wide Range of Compressible Flows on Hybrid Computational Architectures. Supercomputing Frontiers and Innovations 9 (4), 85-99 (2022). https://doi.org/10.14529/jsfi220408
21. Tarnavsky G.A., Aliev A.V. Specific features of high-speed flight aerodynamics: computer simulation of hypersonic flow around the head of an object. Numerical Methods and Programming 9 (4), 371-394 (2008). (In Russian)
22. Rodionov A. V. Development of methods and programs for numerical simulation of non-equilibrium supersonic flows in application to aerospace and astrophysical problems. Thesis for the degree of Doctor of phys- and mathem. Sciences. Saratov, Institut prikladnoi matematiki im. M. V. Keldysha RAN. (2019). (In Russian)
23. van Leer B. Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method. Journal of Computational Physics 32 (1), 101-136 (1979). https://doi.org/10.1016/0021-9991(79)90145-1
24. Levy D. W., Powell K. G., van Leer B. Use of a rotated Riemann solver for the two-dimensional Euler equations. Journal of Computational Physics 106 (2), 201-214 (1993). https://doi.org/10.1016/S0021-9991(83)71103-4
25. Ren Yu-X. A robust shock-capturing scheme based on rotated Riemann solvers. Computers and Fluids 32 (10), 1379-1403 (2003). https://doi.org/10.1016/S0045-7930(02)00114-7
26. Sanders R., Morano E., Druguet M.-C. Multidimensional Dissipation for Upwind Schemes: Stability and Applications to Gas Dynamics. Journal of Computational Physics 145 (2), 511-537 (1998). https://doi.org/10.1006/jcph.1998.6047
27. Semenov A. N., Berezkina M. K., Krasovskaya I. V. Classification of shock wave reflections from a wedge. Part 1: Boundaries and domains of existence for different types of reflections. Technical Physics 79 (4), 46-51 (2009). (In Russian) [Engl. trans.: Technical Physics 54 (4), 491-496 (2009). https://doi.org/10.1134/S1063784209040082].
28. Semenov A.N., Berezkina M.K., Krasovskaya I. V. Classification of shock wave reflections from a wedge. Part 2: Experimental and numerical simulations of different types of Mach reflections. Technical Physics 79 (4), 52-58 (2009). (In Russian) [Engl. trans.: 54 (4), 497-503 (2009). https://doi.org/10.1134/S1063784209040094].
Received: March 7, 2023 Revised: May 2, 2023 Accepted: May 18, 2023
АИ^ОГ'Б ¡п£огтаМоп:
Georgy V. ЭНоею — shoev@itam.nsc.ru
ХРОНИКА
18 мая 2023 г. на заседании секции теоретической механики им. проф. Н. Н. Поляхова в Доме ученых им. М. Горького (Санкт-Петербург) в рамках мини-симпозиума по актуальным проблемам механики робототехнических систем заслушано два доклада:
1. Доклад доктора физ.-мат. наук, профессора М. П. Юшкова и студентов С. О. Бондаренко и Т. Ю. Тепловой (СПбГУ) «Гашение колебаний одной робо-тотехнической системы».
Краткое содержание доклада:
Изучается движение в горизонтальной плоскости гибкой «руки» манипулятора с грузом. Требуется найти оптимальную управляющую силу, приложенную к основанию манипулятора, гасящую колебания элементов манипулятора при смещении за указанное время основания прибора на заданное расстояние. Гибкая «рука» манипулятора имитируется набором двух стержней со спиральными пружинами. Поставленная задача управления решается с помощью применения принципа максимума Понтрягина при минимизации функционала от квадрата искомой управляющей силы. Приводятся графики управляющей силы и изменения главных координат системы. Обсуждаются достоинства и недостатки использованного алгоритма решения задачи.
2. Доклад кандидата физ.-мат. наук, вед. науч. сотр. Л. А. Климиной (НИИ механики МГУ имени М. В. Ломоносова) «Робот-тримаран, управляемый одним внутренним маховиком».
Краткое содержание доклада:
Рассматривается робот, перемещающийся в воде за счет движения внутренних масс. Корпус аппарата состоит из платформы, частично погруженной в воду, и трех профилированных поплавков, жестко соединенных с платформой. На платформе расположен управляемый маховик с вертикальной осью вращения, а также мотор, приводящий его в движение, контроллер, система датчиков положения и ориентации объекта. Взаимодействие поплавков и платформы с жидкостью описано на основе гипотезы квазистационарного обтекания. Нестационарные эффекты учтены путем введения присоединенных масс. Построено управление маховиком, обеспечивающее существование установившегося режима движения, на котором средняя скорость центра масс робота имеет заданное направление. Исследована зависимость скорости движения от параметров модели. Собран лабораторный макет робота, проведена серия экспериментов.
