CC BY
52
УДК 532.59; 532.595, 533
ПРИМЕНЕНИЕ ПОДХОДА ЛАГРАНЖА К РЕШЕНИЮ ОДНОМЕРНОЙ ЗАДАЧИ РАСПРОСТРАНЕНИЯ УДАРНЫХ ВОЛН В ГАЗЕ
© 2011 В.В. Никонов, В.Г. Шахов
Самарский государственный аэрокосмический университет
Поступила в редакцию 06.04.2011
Рассматривается моделирование одномерной задачи распространения ударных волн в газе с помощью метода, использующего подход Лагранжа к перемещению частиц сплошной среды. В основе метода лежит решение Годунова для задачи о распаде разрыва. Результаты численного решения для двух тестовых случаев сравниваются с точным решением и с результатами, полученными Хартеном. Отмечается, что у предлагаемого метода отсутствует численная диффузия на ударных волнах. Ключевые слова: газ, одномерная волна, численный метод, метод Годунова, подход Лагранжа, шаг по времени, точное решение, начальные условия, численная диффузия
Актуальной проблемой при численном моделировании сверхзвуковых течений газа является корректное воспроизведение скачков уплотнения. В настоящее время для расчета таких течений широко используются методы, построенные на так называемых ТУБ схемах [1, 2], использующих подход Эйлера к рассмотрению движения сжимаемой сплошной среды. Одним из недостатков данных методов является большая численная вязкость, приводящая к размазыванию скачков уплотнения (фронтов ударных волн). Для решения этой проблемы Годунов [3] и др. предлагали использовать подвижные сетки, что существенно усложняет задачу, особенно при переходе к рассмотрению двумерных и трехмерных областей течения. В данной работе, как и в [4], предлагается метод, использующий подход Лагранжа к рассмотрению движения сжимаемой сплошной среды и фиксированную однородную сетку. Как будет показано ниже на тестовых задачах, у предлагаемого метода отсутствует численная диффузия на ударных волнах.
1. ОСНОВНЫЕ ПОЛОЖЕНИЯ И ИДЕЯ МЕТОДА
Уравнения Эйлера, описывающие одномерную задачу движения сжимаемой сплошной среды, в размерных переменных имеют вид [1]
др др ди — + и — + р—= 0
д! д х дх '
Никонов Валерий Владимирович, кандидат технических наук, научный сотрудник лаборатории "Аэродинамика". E-mail: v_nikonov@mail.ru
Шахов Валентин Гаврилович, кандидат технических наук, профессор, заведующий кафедрой аэрогидродинамики. E-mail: shakhov@ssau.ru
ди ди 1 д р
— + и — +--- = 0
д! д х р д х '
дв дв р ди л
— + и — +——= 0 (1)
д! д х р д х ' ^ '
Здесь и - скорость, р - плотность, в - внутренняя удельная энергия, р - давление, х -координата, ! - время.
Система уравнений (1) замыкается с помощью уравнения состояния для идеального газа р = (к-1) рв, (2)
где к - показатель адиабаты.
С.К. Годуновым [3] был предложен метод решения данных уравнений, заключающийся в подстановке в явную по времени центральную конечно-разностную схему решения задачи о распаде разрыва. В данной работе предлагается метод, использующий подход Лагранжа для рассмотрения движения жидкости. Сначала для каждой пары ячеек решается задача о распаде разрыва предложенным в [3] методом, в результате получаем значения "больших" переменных: И - плотность, и - скорость, Е - внутренняя удельная энергия, Р - давление, Б - скорость распространения ударной волны или волны разрежения, - справа и слева от границы ячеек. Величины справа и слева от границы разрыва будем обозначать индексами И и Ь. При этом процессы конвекции и "акустики" рассматриваются отдельно, т.е. производится расщепление по физическим процессам, так как местная "акустическая" скорость может в разы отличаться от конвекционной скорости течения. Если шаг по времени для расчета процесса конвекции оказывается меньше шага по времени для процесса "акустики", то конвекция не рассчитывается до первого расчёта этапа "акустики", а за-
тем выполняется столько расчетов этапа конвекции, сколько было пропущено.
2. СХЕМА МЕТОДА ДЛЯ ЭТАПА"АКУСТИКИ"
то для всех ячеек сетки, для которых выполняется условие x1 < xJ < x2, решение в следующий
момент времени к будет тривиальным P = к-1 + P к-1)
1 Rj Ri т 1 Ri+ 1 / ,
На этапе "акустики" в численном методе аналогично работе [4] рассматриваются две волны, бегущие влево _и_ вправо, причём левая переносит величины ={рь, иь,Яь, Бь) с местной "акустической" скоростью Сь = -иь , а правая - W2={PR, ик,Як, Бк) со скоростью Ск = DR-UR.
Предварительно для каждой пары ячеек 1 и 1+1 с набором данных (р^и^р^} и {ры,иы, pi+1,8i+1} решается задача о распаде разрыва предложенным Годуновым в [3] методом, в результате получаются значения больших переменных {р, и,.ь,Бь^ь, которые записываются в ячейки сетки следующим образом
С! = Р , иь1+11 = и , ^ = ^ ,
'Ь - , (3)
они составляют вектор ^ волны, бегущей в отрицательном направлении оси х (влево); и
р.к-1=I5, и.1-1=и, ..к-1=..,
Е.к-1 = Е., С.к-1 = DR-U , (4)
которые, в свою очередь, составляют вектор w2 волны бегущей в положительном направлении оси х (вправо). В случае волны разрежения переменой Б присваивается скорость более медленной характеристики волны разрежения. Например, для правой волны разрежения
DR = и + (ц+1-Ц), (5)
Ек 1 = Е С к 1 = D -и
Ь1+1 ^Ь ' Ь1+1 ^т ^ '
где
Ci+1 =
К-
■4+1
(6)
Рассмотрим волну, бегущую в положительном направлении оси х (вправо). Ячейки рассматриваются попарно. Координаты ячеек после перемещения на этапе акустики для "правой" волны определятся как
= + С к-1 + Х1 = Х1 + CRi д 1с ,
= + С к-1 +
Х2 = Xi+1 + CRi+1 Д 1с
(7)
Решения задачи на этапе "акустики" ищется в следующем виде:
1) Если выполняются условия
Р 11 ^ ^1+1 Rl
< е
1
С к-1 - С к-1
CRi+1 CRi
< е1, (8)
и к = 1(и 1-1 + и к-1)
и ^ Ri ^ и Ri+1 >
RRk=^к"+^к^1),
ERk = ^к-1+ERk:1l)
(9)
Здесь е1 = 10-5 - допуск, связанный с погрешностью вычисления "больших" величин в итерационном методе Годунова. 2) Если выполняется условие
(CRk+l1-CRk-1 )< е1, (10)
тогда
2.1) если
(Х2-Х1) > е2, (11)
где е2 = 10 8 - допуск, связанный погрешностью машинной арифметики, то имеем "слабую" ударную волну, которая не обгоняет решение из впереди стоящей ячейки, поэтому решение определится как
р^к=р*к-1, ^к=, RRk=RRk-l ^ = ERk-1, если Х1 < xj < (х1 + х2)/2 , р к = р к-1 и к = и к-1 те к = те к-1
Rj +1 , +1 , +1 ,
ERk = ERk+-l1, если (х1 + х2)/2<х <х2; (12)
2.2) если
(Х2-Х1) < е2, (13)
то решение определится как
р к = р к-1 и к = и к-1 те к = те к-1
PRj = ^ , ^ = ^ , ^ = ,
ERk = ERk-1, если (^^ -URk-1) < 0 и
(х^Ь/2) < < (х1 + Ь/2), (14)
так как имеем ударную волну, и
р к = р к-1 и к = и к-1 те к = те к-1
^ +1 , +1 , Rj +1 ,
ERk = , если (URk+-l1 - URk-1) > 0
и (х2-ь/2) < Xj< (Х2 + Ь/2) , (15)
так как имеем волну разрежения.
3) Если (CRk+l1
а) если ную волну
-С к-1 )> е 1 — с
1, то
/т т к-1 т т к-1 \ а) если I - UR1 ] < е1, то имеем удар-
Р к : _ Р к -1 _ ГЯ1 , Ияк - Ияк-1, Яяк - ^ -1
Е к к-1 - ЕЯ1 , если (х1 -Ь/2) <х, < (х1 + Ь/2),
Р к : _ Р к -1 ГЯ1+1 , Ияк - Ияк-1, ..к - Я к -1 ■ ЯЯ1+1
Е к ЕЯ, - Е к-1 _ Я1+1 , если (х1 + Ь/2) < х, <х2. (16)
Данный вариант возможен в силу следующего отношения "акустической" скорости ударной волны к скорости звука невозмущенного течения перед её фронтом, которое будет иметь вид к -1
СяГ _ бя - И _ рп
С
к -1 К1+1
Я.
/(к + 1)Р + к-1
2к рп
4к
(17)
или
Б. - И рп (к + 1) к + 1 - (к -1) ^ р Яя к-1
С11 ЯЯ ^ к (к + 1) -(к-1) 1 ЯЯ V 2к ,(18)
где величины с индексом II относятся зоне перед фронтом ударной волны.
Из графического анализа данного отношения следует, что для слабых ударных волн оно будет меньше единицы, а потом неограниченно возрастать.
В частности, тестовый случай ударной волны, рассмотренный Рое [1], относится к вариантУ а). / к-1 к-А
Ь) Если I иш+1 - ИЯ1 I > е1, то имеем волну разрежения и решение для х1 < xJ < х2 определяется из условия сохранения инвариантов Римана
и.к = и.к-1
и к-1 - и к-1
(xJ - х1)
Е к -
..к -
х^ (х,- х1) ^
(19)
( Т, к л:
Я
Я1
+ Я
Я1+1
1 л
к Лк-1
^Я1+1
/2
Ряк - (к- 1)ЯякЕяк.
При выполнении условий 1) - 3), если оказывается, что в ячейку сетки на данном подшаге уже было записано решение, то оно заменяется "новым", когда давление "нового" решения больше давления решения уже имеющегося в ячейке сетки. Кроме того, при записи решения в ячейки сетки также проверяется условие, чтобы распространение данных решений с местной "акусти-
ческой" скоростью не обгоняло (не затирало) решение впереди идущей ударной волны, если таковая существует.
Для волны, бегущей в отрицательном направлении оси х (влево) и переносящей величины
w1 ={РЬ, Иь ,ЯЬ, Еь }, условия и выражения, записанные выше, получаются аналогичным образом.
После перемещения левой и правой волн в каждой ячейке сетки мы имеем наборы величин
= иДЯД ЕД и = {РЯк, ИЯк,ЯЯк, ЕЯк} , для которых снова решается задача о распаде разрыва предложенным в [3] методом. Причём в начальных данных задачи значения ^ считаются расположенными слева от границы разрыва, а значения Wl - справа. В результате снова получаем значения "больших" переменных
{Рк*, Ик*,Яьк\ Еьк*Дяк\ ЕЯк*} .
Так как давление и конвективная скорость не испытывают скачка в распаде разрыва, то их значения после подшага "акустики" определятся как
рк - рк*, ик - Ик*.
(20)
Значения для плотности и энергии выбираются исходя из физического смысла решения. Если Р1к* > рк-1, то выбирается решение для ударной волны. В случае если имеем две ударные волны или пограничный случай, когда ударные волны отсутствуют, то решение находится следующим образом
рк=ЯЯк*, вк=ЕЯк*, если Я^*>ЯЯк*>рк-1 или ЯЯк*>рк-1>Я^*, рк=^Г, ек=Е1Г, если ЯЯк*>Я^*>рк-1 или Як*>рк-1>ЯЯк*. (21)
Если Р^* < рк-1, то выбирается решение для волны разрежения. В случае если имеем две волны разряжения или пограничный случай, когда волны разряжения отсутствуют, то решение находится следующим образом
рк=ЯЯк*, вк=ЕЯк*, если рк-^Я/Ж^* или Я^*>рк-1>ЯЯк*, рк =RLk*, ек=Е1Г, если рк-1>Я^*>ЯЯк* или ЯЯк*>рк-1>Я^*. (22)
3. СХЕМА МЕТОДА ДЛЯ ЭТАПА КОНВЕКЦИИ
В предлагаемом методе на этапе конвекции ячейки также рассматриваются попарно. Координаты ячеек после перемещения на этапе конвекции определяются как
х1 - х1 + и1
к-1
Д ^ , х2 - х1+1 + ик+-; д 1и . (23)
Решения задачи на этапе конвекции отыскиваются в следующем виде:
1) если выполняются условия
Рм-Ц <
е , |р1+1 "Р^ < е1, |в1+1 < е, (24)
то для всех ячеек сетки, для которых выполняется условие х1 < xJ < х2, решение в следующий момент времени к будет тривиальным
3) Если
(и1+1-и1) > е,, (34)
имеем волну разрежения и решение для х1 < xJ < х2 определяется из условия сохранения инвариантов Римана
к к-1 . и1+1 и1 / \ ^ = и +——Мх-х1)
Рк=^2(рк-1+рк-1), ик = ,2(ик-1+ик-1)
1
рк = 2(рк-1 +Р1+1
к-1
1
2
), е к=- (ек-1+ег-'). (25)
При записи решения (25) в ячейки сетки отдельно проверяется условие, чтобы распространение данного решения со скоростью конвекции не обгоняло (не затирало) решение впереди идущей ударной волны, если таковая существует. Для этого в начале этапа конвекции определяются положения границ всех ударных волн в следующий момент времени.
2) Если выполняется условие
(и1+1-и1) < е1 , (26)
то имеем ударную волну а) если
(х2-х1) > е2, (27)
то решение определится как
Рк=Рк-1 к , и = к-1 к к-1 :и1 , Р: =Р1 ,
е к =ек-1, если х1 < х: < (х1 + х2)/2 ,
Рк=Рк+-1 к , и = к-1 к к-1 : и1+1 , Р: = Р1+1
Р к = е к-1 ^ _ +1 , если (х1 + х2)/2 <х: < х2 ;
Ь)если
(х2-х1) < е2
(28) (29)
то
- если
рк+-;
-к-1
Р1
< е1, сначала определяется граница разрыва
х = х1 + Ь/2 + (ик-1 + ик-11)д 1С, (30) затем в ячейку j слева от границы разрыва записывается решение
Рк = Рк-
к к-1 и: = и: ,
„к к-1 Р: =Р1
е к =ек-1, (31)
: 1 : 1 : а в ячейку j+1 справа от границы разрыва запи сываются величины
к-1
к-1
Р: = Р1+1 , и: = и1+1 , Р: =Р
1+1 , е: е1+1 , (32)
- если ся как
Рк+-;
Р1
к-1
> е1, решение определит-
Рк=Рк-1
к к-1 и: = и:
если
к к-1 Р: =Р1
„к к-1 е: = е.
„к-1 „ (^-Ь/2) < х: < (х1 + Ь/2),
Р1 > Р1+1 к „к-1 к
к-1
Р: = Р1+1 , и: = и1+1 ,
J
, к-1
,к-1
Р: Р1+1 , е: е1+1 ,
если Рк-1 < Рк+-1 и (х2-Ь/2) <х: <(х2 + Ь/2). (33)
ек =
Рк =
к-1 , Л^М ~ Уек 1
■(х,"х1)
(35)
( ± ( ек
к-1
Р1
ч*"у
, ~к-1 + РМ
( е к ^к -1 ек-1
Чь 1+1 У
/2
Рк = (к -1)Рк ек.
При выполнении условий 1)-3), если оказывается, что в ячейку сетки на данном подшаге уже было записано решение, то оно заменяется "новым", когда давление "нового" решения больше давления решения уже имеющегося в ячейке сетки.
4. ТЕСТИРОВАНИЕ МЕТОДА
Предлагаемый метод тестировался на двух широко известных задачах, отличающихся начальными условиями (НУ).
Первая тестовая задача была в свое время предложена Рое [1] и имеет следующие НУ:
Р0 = 1, и0 = 0, Рк = 1, если х: < 0,
Р0 = 0.1, и0 = 0, Р0 = 0.125 , если х: > 0. (36)
Область моделирования принималась равной х е [-4.5,5.5], сетка содержала 100 ячеек, шаг сетки составлял Ь = 0.1. Эти данные аналогичны параметрам моделирования задачи в работе [2].
Шаг по времени для процесса "акустики" выбирался согласно критерию Куранта-Фридри-ха-Леви
= к ь л = ко"с,
(37)
где С = Б - и - "акустическая" скорость распространения ударной волны; кс = 1, т.е. ударная волна на этапе "акустики" перемещалась на однуячейку сетки.
Шаг по времени для процесса конвекции выбирался по аналогичному правилу
Ь
Л, = к,, —, И'
и и
(38)
2
где и - конвективная скорость распростране- Этап "акустики" выполнялся через 6 шагов по
ния ударной волны; ки = 2 , т.е. ударная волна времени, т.е. д ^ = 6 д 1, этап конвекции - через
на этапе конвекции перемещалась на две ячейки 10 шагов д 1;и = 10 д 1.
сетки. Результаты расчетов после 110 маршевых
Маршевый шаг равнялся д 1 = 0.02020929 . шагов по времени показаны на рис. 1-3 в сравне-
Рис. 1. Распределение плотности в задаче о распаде разрыва (36) точное решение Годунова [3], О - решение Хартена [2], - предлагаемый метод
и
1
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
тлттгкп ГТТ)1(Г] кГц г ТГГ>г-.
С
^ ПОШИТО — да 10 №№
-4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5 5.5
Рис. 2. Распределение скорости в задаче о распаде разрыва (36): --точное решение Годунова [3], о - решение Хартена [2], - предлагаемый метод
^ ЩЩ XI) иш \\и иш с
с
■ШЛД1 Ши 1 Цл1)
-4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5 5.5
Рис. 3. Распределение давления в задаче о распаде разрыва (36) --точное решение Годунова [3], о - решение Хартена [2], - предлагаемый метод
нии с точным решением Годунова [3] и данными работы Хартена [2], полученными им для схемы ИЬТ1С.
Графики показывают отсутствие численной диффузии на ударных волнах, в отличие от метода Хартена. Правда следует отметить, что из-за накопления ошибки округления ударная волна запаздывает в показанный момент времени на однуячейку сетки.
НУ второй тестовой задачи [2] задавались следующим образом
р0 = 3.52773 , и0 = 0.69888 , = 0.445 ,
если xJ < 0,
р0 = 0.571, и0 = 0, р0 = 0.5 ,
если х^ > 0. (39)
Область моделирования принималась равной х е [-8,6], сетка содержала 140 ячеек, шаг сетки составлял И = 0.1. Эти данные были также аналогичны параметрам моделирования задачи в работе [2]. Шаг по времени для процесса "акустики" выбирался согласно критерию (37), где кс = 2, т.е. ударная волна на этапе "акустики" перемещалась на две ячейки сетки. Шаг по времени для процесса конвекции выбирался по правилу (38), где ки = 3 , т.е. ударная волна на
этапе конвекции перемещалась на три ячейки сетки. Маршевый шаг равнялся д 1 = 0.2 , этапы "акустики" и конвекции выполнялись на каждом шаге по времени, т.е. д 1;с = д 1;и = д 1.
Результаты расчетов в момент времени 1 = 2.0 показаны на рисунках 4-6 в сравнении с данными работы Хартена [2], полученными им для схемы иЬТ1С.
Графики показывают отсутствие численной диффузии на ударных волнах, в отличие от метода Хартена. Как и в первом примере, из-за накопления ошибки округления задний фронт ударной волны запаздывает в показанный момент времени на одну ячейку сетки.
В заключение можно сделать следующие выводы:
1) Преимуществом предлагаемого метода по сравнению с методом Хартена является отсутствие численной вязкости (диффузии) на ударных волнах.
2) С течением времени ударные волны или волны разрежения могут распространяться несколько медленнее или быстрее, чем в точном решении. Это происходит из-за округления положения фронтов волн с точностью до ячейки сетки вследствие применения фиксированной расчетной сетки.
1.4
1.2
0.8
0.6
0.4
0.2
ЛЛВДИ с
3 3
(ВШТОШШМ с ) Э дашвз
'*™?Й11гттМ"ШТ>П (ТГШШ Н1+Ш+1Л1 мнЦЧНГГ
Ш ЩЦШЛ!) СШЛ
-8 -6
X
Рис. 4. Распределение плотности в задаче о распаде разрыва (39) точное решение Годунова [3], О - решение Хартена [2], - предлагаемый метод
1.8 1.6 1.4 1.2 1
0.8 0.6 0.4 0.2 0 -0.2
с ■>
т
0 э
(ДШШШШ
о
и
-6
-4
-2
Рис. 5. Распределение скорости в задаче о распаде разрыва (39) точное решение Годунова [3], О - решение Хартена [2], + - предлагаемый метод
Р
ц)1 IflilKljII Ifflllt! i)i (MB (ш щ?1 HJki til □11Щ1(ВЩ)11 ШСЩ) с
-1
J
3.5 3 2.5
1.5
0.5
X
Рис. 6. Распределение давления в задаче о распаде разрыва (39) точное решение Годунова [3], О - решение Хартена [2], - предлагаемый метод
СПИСОК ЛИТЕРАТУРЫ
Roe P.L. Approximate Riemann solvers, parameter vectors, and difference schemes // J. of Comp. Phys. V. 135. 1997. P. 250-258.
Harten A. High resolution schemes for hyperbolic conservation laws // J. of Comp. Phys. v. 49. 1983. P. 357-393.
3.
4.
Численное решение многомерных задач газовой динамики / С.К. Годунов, А.В. Забродин, МЯ. Иванов, А.Н. Крайко, Г.П. Прокопов. М.: Наука, 1976. 400 с. Никонов В.В., Шахов В.Г. Применение подхода Лаг-ранжа к решению одномерной задачи распространения волн в газе в рамках применимости адиабатического закона // Известия Самарского научного центра РАН. 2009. Т. 11. № 3. С. 33-37.
LAGRANGE APPROACH APPLYING TO THE SOLUTION OF ONE-DIMENSIONAL SHOCK WAVES PROPAGATION PROBLEM IN AIR
© 2011 V.V. Nikonov, V.G. Shakhov
Samara State Aerospace University
One-dimensional problem simulation of shock waves propagation in gas using Lagrange approach to fluid particles movement is considered. The method uses Godunov's solution for Riemann problem as basis. The numerical solution results for two test cases are compared with exact one and Harten's data. It is shown, that proposed method doesn't have numerical diffusion for shock waves.
Key words: gas, one-dimensional wave, numerical method, Godunov method, Lagrange approach, time step, exact solution, initial conditions, numerical diffusion.
Valeriy Nikonov, Candidate of Technics, Research Fellow at the Aerodynamics Laboratory. Е-maü: v_nikonov@mail.ru Valentin Shakhov, Candidate of Technics, Professor, Head at the Aero-Hydrodynamics Department. E-mail: shakhov@ssau.ru
