CC BY
75
УДК 533.6.011.51 +533.6.011.72+532.529.5
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ РАСПРОСТРАНЕНИЯ ЦИЛИНДРИЧЕСКИХ УДАРНЫХ ВОЛН В ГЕТЕРОГЕННЫХ СРЕДАХ
Ю.М. Ковалев, Е.С. Шестаковская
Проведено численное исследование распространения цилиндрической ударной волны и ее взаимодействие с прилегающим слоем гетерогенной двухфазной среды. Для решения уравнений сохранения двухфазной среды была применена модификация метода Неймана - Рихтмайера, адаптированная для расчетов течений гетерогенных сред. Сложности расчета течений, содержащих разрывы типа ударных волн, связаны с необходимостью постоянного пересчета расчетных сеток в области высоких градиентов давления. Совокупность задач: построение математической модели двухфазной среды, модификация численного метода расчета, алгоритмы перестройки расчетных сеток представлены в виде законченной компьютерной технологии, позволяющей эффективно проводить расчеты многомерных задач по распространению сильных разрывов в гетерогенных средах.
Ключевые слова: численное исследование, математическая модель, гетерогенная среда, расчетная сетка, ударная волна.
Введение
Исследование процессов, связанных с изучением распространения цилиндрических ударных волн в гетерогенных слоях, представляет практический интерес по ряду причин: оценка последствий взрывов болидов в атмосфере Земли, анализ эффективности воздействия ударных волн с тепловыми и химическими неоднородностями [1, 2], разработка систем гашения ударных волн [3] и т. д.
Сложность рассматриваемых в работе процессов связана с необходимостью решения комплекса задач: построение математической модели двухфазной среды, модификация численных методов расчета, адаптированных для расчета гетерогенных сред, разработка алгоритмов перестройки расчетных сеток, представленных в виде законченных компьютерных технологий, позволяющих эффективно проводить расчеты многомерных задач по распространению сильных разрывов в гетерогенных средах. Решение всего этого комплекса задач будет представлено ниже на примере распространения цилиндрической ударной волны (УВ) и взаимодействия ее с гетерогенным слоем.
1. Постановка задачи и математическая модель
В настоящей работе проведено численное исследование взрыва цилиндрического заряда взрывчатого вещества (ВВ) (рис. 1) над плоской поверхностью (поверхностью Земли) и над плоской поверхностью с прилегающим к ней слоем гетерогенной среды (поверхность Земли, покрытая лесом; городской массив и т. д.) с различной объемной долей конденсированной фазы (к-фазы).
Пусть цилиндрический заряд ВВ расположен вдоль оси Oz на определенной высоте h над поверхностью Земли (см. рис. 1). Вблизи поверхности Земли располагается слой двухфазной гетерогенной среды высотой И. Задача является симметричной относительно плоскостей хОу и yOz, поэтому будем рассматривать только область х > 0.
Законы сохранения, описывающие динамику распространения УВ в двухфазной среде с неподвижной, неде-формируемой к-фазой и инвариантные относительно преобразования Галилея [4], в лагранжевых массовых переменных имеют следующий вид
1 d а _ди а dt дх ду’
(1)
(2)
(3)
dE dа
(4)
А ’ Р
Здесь р - плотность газа; р - его давление; а _ 1/р - удельный объем; Е - удельная внутренняя энергия газа; у - показатель адиабаты; и и V - х и у компоненты вектора скорости соответственно (V _ {и,V} ); f; и ^ - х и у компоненты силы межфазного взаимодействия; а - объемная доля к-фазы.
Система уравнений (1)-(4) решалась численно модифицированным методом Неймана - Рихт-майера [5]. На нижней и левой границах расчетной области (в силу симметрии задачи) задавались граничные условия для жесткой стенки.
В расчетах моделировался взрыв цилиндрического заряда диаметром В = 3 см, находящегося на высоте h = 1 м от поверхности земли. Вблизи поверхности земли задавался слой двухфазной гетерогенной среды высотой И = 0,5 м, так что взрыв заряда происходил над гетерогенным слоем. Диаметр частиц к-фазы выбирался равным d = 0,1 мм. При моделировании заряда ВВ были использованы данные работы [6]. Расчеты проводились для различных объемных концентраций а частиц к-фазы в гетерогенном слое.
2. Численный метод
Классическим методом решения уравнений механики сплошной среды для разрывных течений является разностный метод Неймана - Рихтмайера [5]. Данный метод является методом сквозного счета и для дистракции (сглаживания) сильных разрывов в нем используется псевдовязкость ю (искусственная вязкость) как добавка к давлению, а переменные Лагранжа позволяют легко описывать контактные разрывы.
В данной работе проведено обобщение подхода, изложенного в работе [5], на случай двухфазной гетерогенной среды.
Запишем систему уравнений в разностном виде. Полагаем, что плотность, давление, внутренняя энергия, сила сопротивления со стороны к-фазы и псевдовязкость относятся к центрам ячеек, а координаты и скорости - к узлам сетки. При интегрировании по времени используем схему с перешагиванием, когда скорости считаются на полушаге по времени ^+1/2, а остальные величины - в моменты времени {'.
Используя формулу Грина для интегрального определения частных производных, запишем уравнения движения (2) следующим образом
Для построения разностной схемы уравнения (5) и (6) интегрируются по некоторой площадке, состоящей из частей ячеек, окружающих рассматриваемый узел сетки, причем давление и плотность считаются постоянными в каждой ячейке. При этом получаются уравнения сохранения импульса для каждой указанной ячейки.
Пусть в расчетной области построена четырехугольная сетка. В некоторый момент времени Г
п-1/2 п-1/2 п п гп гп _п п п т?п
известны величины и , V , х , у , /х , Jy во всех узлах сетки и величины а , ю , р , Е - в центрах четырехугольников, образующих сетку. Получим разностные соотношения для определения этих величин внутри и на границе расчетной области в моменты времени t = Г + х/2 - для скоростей и ^ = {' + х - для остальных величин, где т - шаг по времени.
(5)
(6)
Внутренний узел сетки. Любой внутренний узел сетки окружен четырьмя четырехугольниками (рис. 2). При интегрировании уравнений (5), (6) для данного узла в качестве площадки интегрирования S выберем четырехугольник, ограниченный линиями, проходящими через середины общих сторон и отсекающими от соседних с узлом ячеек четверть площади.
Будем считать, что сумма р + ю постоянна внутри каждой ячейки. Заменив производные по времени разностными отношениями, получим для определения скоростей в момент ^+12 в точке {хо,уо}:
Рис. 2. Расчетная схема
и+1/2 _.и-1/2
40
= и,
0
2Ат0 ;=1
£(РП +<)(У+1 -УП)+"Frm,
4
П+1/2 _ . 1-1/2 . Т V"' / п . ,^п\( п п\ .
уо = уо + 2Ат ^( Р1 + Юг)(- х,)+
■тК
0 і=1
п—1/2 У ,
(7)
(8)
где при i = 4 индекс i = 1 полагаем равным единице. Масса т0, связанная с узлом, определяется по формуле
1 4
Ат0 ={р^ = - £р^”, (9)
5 4 /=1
где - площадь /-го четырехугольника. Величина т0 остается неизменной в ходе расчетов, но
меняется при перестройке сетки, рассмотренной далее.
В уравнениях (7), (8) введены обозначения Fx и Fy для компонент удельной силы межфазного взаимодействия
К
п—1/2
( f У
Ух
п—1/2
У
л
р
= —г сг, а*/| и,
с, ^(иоп—1/2 )2+(V;-1/2)
С,^(и0-1'2 )2 + (<-1'2 )2 • V
п—1/2
о ,
= ——с, а*/| и,
п—1/2 0 .
(10)
(11)
Добавки (10) и (11) к разностным уравнениям движения (7), (8) учитывают наличие конденсированной фазы с объемной долей а. Из уравнений (3) получаем новые координаты узлов сетки
(12)
и+1
х0-- = х0 + ти0
п+1/2 п+1 п . _„п+1/2
, Уо = Уо + "^о .
Граничные узлы. На границах расчетной области используется условие на жесткой стенки. Вдоль жестких стенок возможно только скольжение вещества. В численном алгоритме [5] это достигается за счет введения фиктивных ячеек.
Уравнение энергии и уравнение состояния (5), как и в [5], записываем в разностном виде и
п+1 г7п+1
разрешаем относительно р и Е
2
рп+1 =1
Еп+1 = -
I / I I ,~.п+1\/—И+1\
-р а + ^р + ю + ю Да —а ^
У + 1 п+1 п
а — а
У—1
^п+1
(у — 1)р
п+1
(13)
(14)
х
Р
о
о
Для определения плотности р”+ в новый момент времени воспользуемся тем, что полная масса газа, заключенного в ячейке, не меняется в ходе расчетов Sn
и+1 п
р =р
£
п+1
(15)
где £ - площадь рассматриваемой счетной ячейки (четырехугольника).
Выбор шага х интегрирования уравнений по времени осуществляется исходя из условия устойчивости Куранта [5], где в качестве линейного масштаба расчетной ячейки выбирается d шт -наименьшая диагональ четырехугольника.
В двумерных расчетах с использованием лагранжевых координат сетка может сильно искажаться. Если не проводить перестройку сетки, расчеты на искаженной сетке могут дать очень большую погрешность, что, как правило, приводит к неустойчивости решения.
В данной работе используется алгоритм перестройки сетки, предложенный в работе [7]. Перестройка сетки осуществляется на всей расчетной области. Пересчет термодинамических параметров среды со старой сетки на новую сетку основан на применении законов сохранения массы, импульса и энергии, для вычисления скоростей в узлах новой сетки используется закон сохранения импульса [7]. Реализация на практике описанного в [7] метода перестройки сетки показала, что он обеспечивает выполнение законов сохранения массы, импульса и энергии с точностью не хуже 1 % даже при решении задач, требующих выполнения нескольких десятков процедур перестройки. В качестве критерия перестройки сетки в наших расчетах выступало резкое уменьшение временного шага. Тестирование разработанной программы расчетов было проведено на решении задач, предложенных в работе [5], и получено удовлетворительное совпадение результатов.
Рассмотрим взаимодействие УВ с гетерогенным слоем при а = 0,001 (рис. 3).
1
32.00
30.00
28.00 26.00
24.00
22.00 20.00 18.00 16.00
14.00
12.00 10.00 8.00 6.00
4.00
2.00 0.00
I
10.00
9.00
8.00
7.00
6.00
5.00
4.00
3.00
2.00 1.00
0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50
t=0.65 тс t=1.5 тс
Рис. 3. Эволюция поля давления при взрыве шнурового заряда над гетерогенным слоем а = 0,001
Полученные результаты показывают, что, входя в гетерогенный слой, УВ испытывает торможение. Так как объемная доля частиц мала, УВ частично отражается от поверхности слоя, проходит через него и отражается от плоской поверхности. Со временем регулярное отражение переходит в нерегулярное (маховское). При этом ножка Маха формируется на поверхности Земли. Фронт УВ деформируется в гетерогенной области, однако его форма остается близкой к цилиндрической.
При большей объемной доле к-фазы а = 0,01 (рис. 4) наблюдается более существенное отражение УВ от поверхности гетерогенного слоя, и формирование ножки Маха происходит уже на поверхности гетерогенного слоя. В пристеночном гетерогенном слое образуется область повы-
шенного давления. Величина давления в этой области уменьшается существенно медленнее, нежели происходит эволюция УВ в чистом воздухе, тем самым увеличивая передаваемый на поверхность Земли импульс и область разрушений.
п,
t=0.65 тс t=1.5 тс
Рис. 4. Эволюция поля давления при взрыве шнурового заряда над гетерогенным слоем а = 0,01
6.00
5.00
4.00
3.00
2.00
1.00
Заключение
В работе представлена комплексная технология, позволяющая проводить численное исследование динамики распространения цилиндрических УВ в областях, заполненных гетерогенными смесями. Проведение расчетов в лагранжевых переменных позволяет выделять сильные и контактные разрывы, что очень важно при анализе динамики изменения структуры УВ, а пересчет расчетной сетки дает возможность увеличивать расчетную область без потери точности.
Результаты проведенных расчетов показывают, что структура ударной волны в гетерогенной среде становится существенно сложнее, чем в чистом газе. Происходит запирание УВ внутри слоя и увеличение передаваемого УВ поверхности Земли импульса, что приводит к значительному увеличению масштабов разрушения.
Литература
1. Гришин, А.М. Об усилении ударных волн при их взаимодействии с фронтом лесного пожара /А.М. Гришин, Ю.М. Ковалев //ДАН СССР. - 1990. - Т. 312, № 1. - С. 50-54.
2. Ковалев, Ю.М. Ослабление воздушных ударных волн системой решеток / Ю.М. Ковалев, А.Ю. Черемохов // Вопросы атомной науки и техники. Сер. «Мат. моделирование физ. процессов». - 1997. - Вып. 3. - С. 39-43.
3. Ковалев, Ю.М. Взаимодействие плоской ударной волны с нагретым слоем газа /Ю.М. Ковалев, А.Ю.Черемохов //Доклады РАН. - 1999. - Т. 367, № 6. - С. 769-771.
4. Ковалев, Ю.М. Анализ инвариантности относительно преобразования Галилея некоторых моделей математических многокомпонентных сред / Ю.М. Ковалев, В. Ф. Куропатенко // Вестник ЮУрГУ. Сер. «Математическое моделирование и программирование». - 2012. -Вып. 13, № 27 (286). - С. 69-73.
5. Софронов, И.Д. Методика расчета нестационарных двумерных задач газовой динамики в лагранжевых координатах / И.Д. Софронов, Н.А. Дмитриев, Л.В. Дмитриева, Е.В. Малиновская // Теоретические основы и конструирование численных алгоритмов задач математической физики. - М. : Наука, 1979. - С. 175-201.
6. Моделирование взрыва шнурового заряда в пологе леса при отсутствии пожара / В.А. Антонов, А.М. Гришин, Ю.М. Ковалев, Л.Ю. Наймушина // ФГВ. - 1993. - Т. 29, № 4, - С. 115-123.
7. Яловец А.П. Расчет течений среды при воздействии интенсивных потоков заряженных частиц / А.П. Яловец // Журнал прикладной механики и технической физики. - 1997. - № 1. -С. 51-66.
Работа выполнена при поддержке РФФИ грант №13 - 01 - 00072.
Ковалев Юрий Михайлович, д-р физ.-мат. наук, профессор, заведующий кафедрой вычислительной механики сплошных сред, Южно-Уральский государственный университет (г. Челябинск); [email protected].
Шестаковская Елена Сергеевна, канд. физ.-мат. наук, доцент кафедры вычислительной механики сплошных сред, Южно-Уральский государственный университет (г. Челябинск); [email protected].
Bulletin of the South Ural State University Series “Computer Technologies, Automatic Control, Radio Electronics”
2013, vol. 13, no. 3, pp. 102-108
NUMERICAL ANALYSIS OF CYLINDRICAL SHOCK WAVES IN HETEROGENEOUS ENVIRONMENTS
Y.M. Kovalev, South Ural State University, Chelyabinsk, Russian Federation, yum_kov@mail. ru,
E.S. Shestakovskaya, South Ural State University, Chelyabinsk, Russian Federation, [email protected]
In this paper, a numerical study of propagation of cylindrical shock wave and its interaction with the adjacent layer of a heterogeneous two-phase medium is held. To solve the conservation equations of two-phase medium was used a modification of the method of Neumann-Richtmyer, adapt to calculations for heterogeneous environments. The complexity of the calculation of flows containing discontinuities such as shock waves dues to the necessity of translation computational grids in the high pressure gradients. Set of problems: the construction of a mathematical model of two-phase medium, the modification of the numerical method of calculation algorithms for restructuring computational grids are presented in the form of a complete computer technology, which effectively carry out calculations for the dissemination of multi-dimensional problems of strong discontinuities in heterogeneous environments.
Keywords: numerical analysis, mathematical model, heterogeneous medium, the computational grid, the shock wave.
References
1. Grishin A.M., Kovalev Yu.M. About the Strengthening Shock Waves through their Interaction with the Front of a Forest Fire [Ob usilenii udarnykh voln pri ih vzaimodeystvii s frontom lesnogo pozhara] Doklady akademii nauk [Academy of Sciences of the USSR Doklads], 1990, vol. 312, no.1, pp.50-54.
2. Kovalev Yu.M., Cheremokhov A.Yu. The Weakening of the System of Air Shock Wave Gratings [Oslablenie vozdushnykh udarnykh voln sistemoy reshetok] Voprosy atomnoy nauki i tekhniki. Seriya “Matematicheskoe modelirovanie fizicheskikh protsessov” [Problems of Atomic Science and Technology. Series «Mathematical modelling of physical processes», 1997, no. 3, pp. 39-43.
3. Kovalev Yu.M., Cheremokhov A.Yu. The Interaction of a Plane Shock Wave with a Hot Gas layer [Vzaimodeystvie ploskoy udarnoy volny s nagretym sloem gaza] Doklady akademii nauk [Academy of Sciences of the USSR Doklads], 1999, vol. 44, no. 8, pp. 582-584.
4. Kovalev Yu.M., Kuropatenko V.F. Analysis of the Invariance under the Galilean Transformation of Some Mathematical Models of Multi-media [Analiz invariantnosti otnositel'no preobrazovaniya Gali-leya Nnekotorykh Modeley matematicheskikh mnogokomponentnykh sred] Bulletin of the South Ural
State University. Series “Mathematical modelling, Programming & Computer Software”, 2012, no. 27, pp.69-73. (in Russian)
5. Sofronov I.D., Dmitriev N.A., Dmitrieva L.V., Malinovskaya Е.В. Method of Calculation of Non-stationary Two-dimensional Gas Dynamics in Lagrangian Coordinates [Metodika raschyeta nestat-sionarnykh dvumernykh zadach gazovoy dinamiki v lagranzhevykh koordinatakh] V kn.: Teoreticheskie osnovy i konstruirovanie chislennykh algoritmov zadach matematicheskoy fiziki [In a book : Theoretical foundations and construction of numerical algorithms for problems of mathematical physics]. Moscow, Nauka, 1979, pp. 175-201.
6. Antonov V.A., Grishin А.М., Kovalev Yu.M., Naimushina L.Yu. An Explosion Modeling of the Flexing Charge within Wood’s Curtains in Absenting Fire [Modelirovanie vzryva shnurovogo zaryada v pologe lesa pri otsutstvii pozhara] Fizika goreniya i vzryva [J. The Physics of Combustion and Explosion], 1993, vol. 29, no. 4, pp. 115-123.
7. Yalovets A.P. Calculation of Fluid Flow under the Influence of Intense Beams of Charged Particles [Raschet techeniy sredy pri vozdeystvii intensivnykh potokov zaryazhennykh chastits ] Zhurnal prikladnoy vtkhaniki i fiziki [J. of Applied Mechanics and Technical Physics], 1997, no. 1, pp. 51-66.
Поступила в редакцию 9 июня 2013 г.
