CC BY
47
УДК 519.63+532.529.5
БО!: 10.14529/ттр 160203
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЛИЯНИЯ СХЕМНОЙ ВЯЗКОСТИ ЧИСЛЕННЫХ МЕТОДОВ НА ВЕЛИЧИНУ ИМПУЛЬСА, ПЕРЕДАВАЕМОГО УДАРНЫМИ ВОЛНАМИ
Ю.М. Ковалев, Е.Е. Пигасов
В данной работе приводится анализ влияния ширины фронта ударной волны на величину импульса, передаваемого неподвижной твердой поверхности. Профиль ударной волны был рассчитан модифицированным методом крупных частиц, который позволяет проводить расчеты распространения ударных волн без явного введения искусственной вязкости. Показано, что ширина «размазанного» фронта ударной волны, рассчитанная модифицированным методом крупных частиц, не зависит от интенсивности ударной волны. Подобная картина наблюдается, если в численных методах для расчета поведения ударных волн используется в явном виде квадратичная искусственная вязкость. Для проведения серийных расчетов по исследованию передачи импульса ударных волн твердой неподвижной стенке, в данной работе получено аналитическое решение для профиля давления в ударном переходе в случае квадратичной искусственной вязкости в переменных Эйлера. Для ударной волны с треугольным профилем было показано, что величина импульса, передаваемого твердой стенке, не зависит от ширины ударного перехода.
Ключевые слова: численный метод; математическая модель; законы сохранения; ударные волны; число Куранта.
Введение
Активное применение быстропротекающих процессов в различных областях науки и техники требует повышенного внимания к обеспечению безопасности при работе на создаваемых установках. Развитие математических моделей и увеличение быстродействия вычислительной техники позволяет просчитывать многие аварийные ситуации. Появление новых математических моделей, с одной стороны, связано с появлением новых сложных практических задач. С другой стороны, развитие вычислительной техники позволяет получать численные решения для все более сложных математических моделей механики сплошных сред.
Проблема определения импульса, передаваемого ударными волнами твердой стенке, заключается в том, что существующие численные методы размазывают фронт ударной волны за счет схемной вязкости разностного метода либо путем введения искусственной вязкости в уравнения сохранения импульса и энергии. Если ширина ударной волны для реального газа имеет порядок длины свободного пробега в газе, то ширина ударной волны, получаемая в расчетах, может быть 5-10 ячеек. Это значительно больше длины свободного пробега. Появляется проблема достоверности результатов расчетов по определению импульса, передаваемого ударными волнами твердой стенке, и адекватности их реальной ситуации.
Успешное решение многочисленных задач газовой динамики и аэродинамики методом крупных частиц [1] и его модификациями [2 тсI позволяет надеяться на то^ ч^то
применение метода крупных частиц позволит решить задачу о передаче импульса ударной волной, распространяющейся по газу, твердой стенке. Целью дяннои работы является анализ результатов расчетов взаимодействия ударных волн, распространяющихся в газах, с неподвижной твердой стенкой.
1. Анализ профиля ударной волны
Рассмотрим одномерный плоский случай течения газа. Система законов сохранения массы, импульса и энергии в переменных Эйлера имеет следующий вид
I + % = « «
дру друу др дЬ дх дх' дрЕ дрЕу дру ^ у2
-дг + -&■ + дх = 0- Е =е(3>
Система уравнений (1) - (3) замыкается уравнениями состояния
р
(k - 1)р-
(4)
В данной работе не рассматриваются более сложные уравнения состояния [5, 6]. Здесь p,v,p,e, E,x,t - плотность, скорость, давление, внутренняя и полная энер-
k
Уравнения: (1) - уравнение неразрывности; (2) - уравнение сохранения импульса; (3) - уравнение сохранения полной энергии.
Для определения структуры ударного перехода воспользуемся модификацией метода крупных частиц [2], которая позволяет проводить расчеты по распространению ударных волн без явного введения искусственной вязкости даже для газовзвесей [3, 4J при больших значениях числа Куранта. На рис. 1 приведена структура стацио-норных ударных волн для трех значении безразмерного ДсШЛбНИЯ p/p0 = 5,10,20 соответственно, dx - шаг по пространству. На рисунке видно, что ширина ударного перехода не зависит от интенсивности стационарной ударной волны. Следовательно, схемная вязкость модификации метода крупных частиц [2] работает таким образом, что ширина ударного перехода является величиной постоянной, не зависящей от ин-тбнсивности ударной волны.
Крестиками на рис. 1 обозначены точки перегиба кривых в зоне ударного перехода, которые совпадают с точностью до dx/2 с положением стационарных ударных волн в соответствующий момент времени, рассчитанным из условия Гюгонио [7, 8]. Подобная картина наблюдается, если в численных методах для расчета поведения ударных волн используется в явном виде квадратичная искусственная вязкость [8]. Для проведения серийных расчетов по исследованию передачи импульса ударных волн твердой неподвижной стенке, получим аналитическое решение для профиля давления в ударном переходе в случае квадратичной искусственной вязкости в пере-М6ННЫХ Эйлера. В случае переменных Лагранжа аналогичное решение было получено в монографии A.A. Самарского, Ю.П. Попова [8]. С этой целью выделим в явном
Рис. 1. Структура стационарной ударной волны (х/йх)
виде в уравнениях (2) и (3) выражения, связанные с квадратичной искусственной ВЯЗКОСТЬЮ. Получим следующие уравнения
дру друу др дд дЬ дх дх дх'
дрЕ + дрЕу + дру + дду ^ ду ду
дЬ дх дх дх ' ^ ^Рдхдх Коэффициент вязкости ^ считаем постоянным. Считаем, что ударная волна распространяется по газу с постоянной скоростью Б. Начальное состояние газа определяется известными параметрами РО'УО'РО, е0, которые расположены справа от ударной волны и простираются в положительном направлении оси х до Будем считать, что газ перед ударной волной покоится Уо = 0. Полагаем также, что на обеих бесконечностях имеет место установившееся течение
(Щ
х
= 0' (7)
где Щ - любой из параметров: р, у,р, е.
Будем искать решение задачи о структуре фронта ударной волны в виде
Щ(х, Ь) = Щ(в)' в = х - БЬ. (8)
Замена переменных (8) позволяет выразить производные по х и Ь через производную в
д б й д й
дЬ ¿в дх ¿в
и свести систему уравнений в частных производных (1). (5), (6) к системе обьткно-ве н н ых дифференциальных уравнений:
й
-(р(у - Б)) = 0' (10)
( - £ - % м
¿> / ( Б) ( + ¿Ру бюйю
¿8 \ \ 2 у у ¿8 ^^ ¿8 ¿8
Равенство (7) играет роль граничных условий и имеет следующий вид
df ds
= 0. (13)
s—>±оо
Интегрируя уравнения (10) - (12) от 8 = до некоторого текущего значения, получим систему алгебраических уравнений для определения структуры фронта ударной
ВОЛНЫ
р(у - Б) = -роБ, (14)
ру(у - Б) = -р - д - р0, (15)
р(у - Б) ^е + + РоеоБ = -(р + д)у, (16)
= ¿у_^у_ = р
¿8 (к - 1)р
Переходя в системе уравнений (14) - (17) от плотности к удельному объему V =
1/р и исключая из системы уравнений все неизвестные кроме удельного объема V и
д
DV (k + 1)D )
* Jk-Y) = 2(k - 1)V0(Vo - V)
3 V- ((k - 1)V , 2kpoVo2 \
v(k + 1) 0 + (k + 1)D2y
(18)
д
вязкость) равен нулю, то значение удельного объема за ударной волной должно быть равно
(k - 1) 2kpoVo
2
Vi = Vo +
o
(к + 1) " (к + 1)Б2' Подставляя данное выражение в равенство (18), получим
д <-ггщ-т ^ - V XV - V.). <»>
Если в уравнении (19) перейти к переменным Лагранжа для диссипативного члена и массовой скорости ударной волны, то уравнение (19) совпадет с уравнением для
диссипативного члена, полученным в монографии А.А. Самарского, Ю.П. Попова [8].
д
лучим дифференциальное выражение для определения профиля удельного объема
V (5) _
¿V = М+1(Ъ - V)(У - V.). (20)
ds У 2р
Интегрируя уравнение (20), получим выражение для профиля удельного объема
+ с)
лп \ (Vo + V) , (Vo - V) (k + 1) .
v(s) = —^— + —2—cos п/ o.. s + Ci)' (21)
где Ci - постоянная интегрирования, определяющая положение начала системы координат. Значениям Ci = и s = 0 соответствует точка перегиба профиля удельного объема V (0) = (Vo + V)/0, а значениям Ci = 0 и s = 0- начальное значение удельно-Vo
ширина ударного перехода
4 = \ШУ п (22)
Для известной толщины ударного перехода 4 из выражения (22) определяется значение коэффициента квадратичной вязкости
. = \(k + 1)( С23)
Следовательно, профиль удельного объема для известном толщины ударного пе~ рехода 4 и Ci = 0 имеет следующий вид
V+ + ^ cos ((24)
Зная профиль удельного объема в зоне ударного перехода, легко получить распределение скорости и давления по формулам
v = D(1 - poV(s)), p = po + DpoV(s). (25)
2. Постановка задачи
Рассмотрим задачу о передаче импульса твердой стенке ударной волной с треугольным профилем, состоящим из ударного скачка и примыкающей к нему волны разрежения. Параметры газа на фронте ударной волны связаны с параметрами газа перед фронтом соотношениями Гюгонио. В начальный момент времени t = 0 в области 0 < x < xB параметры газа за фронтом ударной волны в зоне волны разрежения задаются на основе изоэнтропических зависимостей для простой волны Римана с ли-неиным профилем скорости газа
Ч й. (2б)
Остальные величины определяются из условия адиабатичности течения уравнениями
C = Cb - k-^vJ 1 - CB = kPPB, (27)
pB
Р = Рв(^СВ^ , е = ев(^, Р ={к - 1)ре, (28)
где С - скорость звука, а индексом В обозначены газодинамические параметры в точке В.
При ^ > 0 на правой границе х = ха задано уеловие и = 0, на левой границе при х = 0 условие свободного протекания газа.
В данной работе были рассмотрены два случая передачи импульса при взаимодействии ударных волн со стенкой:
1. Ударный скачок располагается непосредственно у стенки (рис. 2);
2. Ударный переход имеет профиль. определенный формулами (24). (25), представлен на рис. 3.
О 0.2 0.4 Об 0.& 1 12 1.4 1.6 1.8 2 22 2.4 2.6 2.8 3 К, М
Рис. 2. Начальное положение ударной волны (случай 1)
2.S6 188 24 1S1 194 2М 258
Рис. 3. Профиль давления в начальный момент времени (случай 2)
Для первого случая значения ха и хв совпадают. На рис. 3 штриховой линией обозначено истинное расположение ударного скачка.
Величина импульса, передаваемая ударными волнами стенке, вычислялась по формуле
I = [ (р - Ро)^. (29)
.70
Расчеты импульса, передаваемого стенке ударными волнами, продолжались до тех пор, пока давление на стенке не становилось равным начальному значению. В результате проведенных расчетов было показано, что величина импульса, передаваемого неподвижной стенке, для случаев 1 и 2 различается менее 2%. Следовательно, можно сделать вывод о том, что величина импульса, передаваемого неподвижной стенке, не зависит от ширины и профиля ударного перехода.
Авторы, выражают свою благодарность профессору В.Ф. Куропатенко за полезные обсуждения и интерес к работе.
Литература
1. Белоцерковский, О.М. Метод крупных частиц в газовой динамике / О.М. Белоцерков-ский, Ю.М. Давыдов. - М.: Наука, 1982.
2. Гришин, Ю.А. Новые схемы метода крупных частиц и использование их для оптимизации газовоздушных трактов двигателей / Ю.А. Гришин // Математическое моделирование. - 2002. - Т. 14, № 8. - С. 51-55.
3. Грищенко, Д.С. Модификация метода крупных частиц для исследования течений газовзвесей / Д.С. Грищенко, Ю.М. Ковалев, Е.А. Ковалева // Вестник ЮУрГУ. Серия: Математическое моделирование и программирование. - 2015. - Т. 8, № 2. - С. 36-42.
4. Ковалев, Ю.М. Анализ некоторых модификаций метода крупных частиц на примере исследования течений газовзвесей / Ю.М. Ковалев, Е.А. Ковалева, Е.Е. Пигасов // Вестник ЮУрГУ. Серия: Математика. Механика. Физика. - 2015. - Т. 7, № 3. - С. 71-77.
5. Ковалев, Ю.М. Уравнения состояния и температуры ударного сжатия кристаллических ВВ / Ю.М. Ковалев // Физика горения и взрыва. - 1984. - Т. 20, № 2. - С. 102-107.
6. Моделирование взрыва шнурового заряда в пологе леса при отсутствии пожара / В.А. Антонов, A.M. Гришин, Ю.М. Ковалев, Л.Ю. Наймушина // Физика горения и взрыва. - 1993. - Т. 29, № 4. - С. 115-123.
7. Куропатенко, В.Ф. Модели механики сплошных сред / В.Ф. Куропатенко. - Челябинск: Челябинский государственный университет, 2007.
8. Самарский, A.A. Разностные схемы газовой динамики / A.A. Самарский, Ю.П. Попов. -М.: Наука, 1975.
Юрий Михайлович Ковалев, доктор физико-математических наук, профессор, заведующий кафедрой «Вычислительная механика сплошных сред>, ЮжноУральский государственный университет (г. Челябинск, Российская Федерация), yum_kov@mail.ru.
Егор Евгеньевич Пигасов, старший преподаватель кафедры «Вычислительная механика сплошных сред>, Южно-Уральский государственный университет (г. Челябинск, Российская Федерация), pigasovee@mail.ru.
Поступила в редакцию 30 марта 2016 г.
MSC 76M25 DOI: 10.14529/mmp160203
Mathematical Modelling of Influence of Circuit Viscosity of Numerical Methods on a Value of the Impulse Transferred by Shock Waves
Yu.M. Kovalev, South Ural State University, Chelyabinsk, Russian Federation, yum_kov@mail.ru,
E.E. Pigasov, South Ural State University, Chelyabinsk, Russian Federation, pigasovee@mail.ru
The analysis of influence of width of the front of a shock wave to the size of the impulse transferred by a motionless firm surface is provided in this work. The profile of a shock wave has been calculated by the modified method of large particles which allows to carry out calculations of distribution of shock waves without obvious introduction of artificial viscosity. It is shown that the width of the "smeared" front of a shock wave calculated by the modified method of large particles doesn't depend on intensity of the shock wave. The similar picture is observed if in numerical methods for calculation of behavior of shock waves square artificial viscosity is used in an explicit form. For carrying out serial calculations for research of transfer of an impulse of shock waves to a firm motionless wall, in this work the analytical solution for pressure profile in shock transition in case of square artificial viscosity in Euler's variables is received. For a shock wave with a triangular profile it has been shown that the size of the impulse transferred to a firm wall doesn't depend on the width of the shock transition.
Keywords: numerical method; mathematical model; conservation laws; shock waves; Courant's number.
References
1. Belotserkovsky O.M., Davydov Yu.M. Metod krupnykh chastits v gazovoy dinamike [Method of Large Particles in the Gas Dynamics]. Moscow, 1982. 392 p. (in Russian)
2. Grishin Yu.A. [The New Scheme of the Method and the Use of Large Particles them for Optimization Gas Engines Paths]. Mathematical Models and Computer Simulations, 2002, vol. 14, no. 8, pp. 51-55. (in Russian)
3. Grishchenko D.S., Kovalev Yu.M., Kovaleva E.A. Modification of Method of Large Particles for Research of Currents of Gas-Suspensions. Bulletin of the South Ural State University. Series: Mathematical Modelling, Programming and Computer Software, 2015, vol. 8, no. 2, pp. 36-42. DOI: 10.14529/mmpl50203 (in Russian)
4. Kovalev Yu.M., Kovaleva E.A., Pigasov E.E. The Analysis of Some Modifications of the Large-Particle Method on the Basis of Research of Gas-Suspension Currents. Bulletin of the South Ural State University. Series: Mathematics. Mechanics. Physics, 2015, vol. 7, no. 3, pp. 71-77. (in Russian)
5. Kovalev Yu.M. Equations of State and Shock Compression Temperature of Crystal Explosives. Combustion, Explosion, and Shock Waves, 1984, vol. 20, no. 2. pp. 219-223. DOI: 10.1007/BF00751596
6. Antonov V.A., Grishin A.M., Kovalev Yu M., Naimushina L.Yu. Modeling Primer Cord Detonation in a Forest Canopy without a Fire. Combustion, Explosion, and Shock Waves, 1993, vol. 29, no. 4. pp. 527-534. DOI: 10.1007/BF00782981
7. Kuropatenko V.F. Modeli mekhaniki sploshnikh sred [Models of Continuum Mechanics]. Chelyabinsk, 2007. 302 p. (in Russian)
8. Samarskiy A.A., Popov Yu.P. Raznostniye skhemy gazovoy dinamiki [The Difference Schemes of Gas Dynamics]. Moscow, 1975. 352 p. (in Russian)
Received March 30, 2016
