Math-Net.Ru
Общероссийский математический портал
Н. Я. Моисеев, Т. А. Мухамадиева, Метод Ньютона для решения задачи
о распаде произвольного рарыва в средах с уравнениями состояния общего
вида, Ж. вычисл. матем. и матем. физ., 2008, том 48, номер 6, 1102–1110
Использование Общероссийского математического портала Math-Net.Ru подразумевает, что вы прочитали и согласны с пользовательским соглашением
http://www.mathnet.ru/rus/agreement
Параметры загрузки:
IP: 185.201.44.167
13 декабря 2024 г., 12:54:13
ЖУРНАЛ
ВЫЧИСЛИТЕЛЬНОЙ
МАТЕМАТИКИ
И МАТЕМАТИЧЕСКОЙ
ФИЗИКИ,
2008, том 48, M 6, с.
1102-1110
УДК 519.634
МЕТОД Н Ь Ю Т О Н А ДЛЯ Р Е Ш Е Н И Я З А Д А Ч И
О РАСПАДЕ ПРОИЗВОЛЬНОГО Р А З Р Ы В А В СРЕДАХ
С УРАВНЕНИЯМИ СОСТОЯНИЯ ОБЩЕГО ВИДА
© 2008 г. Н. Я. Моисеев, Т. А. Мухамадиева
(456770 Снежинок, Челябинская обл.,
ул. Васильева, 13, а.я. 245, ФГУП РФЯЦ-ВНИИТФ им. ак. Е.И. Забабахина)
e-mail:
[email protected]
Поступила в редакцию 19.10.2007 г.
Переработанный вариант 12.12.2007 г.
Предлагается подход к решению задачи о распаде произвольного разрыва в средах с нормаль­
ными уравнениями состояния, основанный на методе Ньютона. Для эффективного вычисле­
ния интегралов Римана используется кубическая аппроксимация изэнтропы, обеспечившая
по сравнению с методом Симпсона более высокую точность, скорость сходимости и эконо­
мичность. Возможности подхода демонстрируются на примерах решения задач для сред, под­
чиняющихся уравнению состояния Ми-Грюнайзена, для которого получены в явном виде ал­
гебраическое уравнение изэнтропы и некоторые точные решения для конфигураций с волна­
ми разрежения. Библ. 16. Фиг. 3. Табл. 1.
Ключевые слова: уравнение газовой динамики, задача о распаде произвольного разрыва,
численный метод Ньютона, уравнения состояния Ми-Грюнайзена.
1. В В Е Д Е Н И Е
Задача Римана о распаде произвольного разрыва для уравнений газовой динамики хорошо из­
вестна и подробно описана в научной литературе, например в [1], [2]. Решения задачи широко ис­
пользуются в численных методиках и во многих из них являются массовыми операциями (см. [3]).
Поэтому особое внимание уделяется эффективности ее решения, как точного (см. [4]-[8]), так и
приближенного (см. [9], [10]).
Согласно универсальному подходу из [4], уравнения состояния (УPC) записываются в двухпараметрическом виде, что приводит к вычислению этих параметров из решения систем уравне­
ний, соответствующих реализовавшимся конфигурациям распада разрыва. З а т е м вычисляются
все термодинамические величины. Системы решаются методом Ньютона в конфигурациях с
ударными волнами (УВ). Интегрирование адиабаты Пуассона сводится к решению системы диф­
ференциальных уравнений для параметрических функций. В [5] нахождение термодинамических
величин сводится к решению систем уравнений методом Н ь ю т о н а в конфигурациях с УВ. В кон­
фигурациях с волнами разрежения (BP) решаются д и ф ф е р е н ц и а л ь н ы е уравнения методом Рунге-Кутты. В [6] изложен подход к решению задачи в о б щ е м виде без подробного описания кон­
кретных алгоритмов. В [8], [10] приводится подробное решение задачи в средах с простыми У PC
методом Ньютона и методом обратной параболической интерполяции (см. [11]) соответственно.
Метод обратной параболической интерполяции обеспечивает кубическую скорость сходимости.
Анализ алгоритмов из [4]-[6] показывает, что они являются достаточно трудоемкими, чтобы
их можно было применять как массовые операции в численных методиках. Так, если сравнивать
эффективность решения задачи в средах с идеальным газом методами из [5] и [10], то решения в
[5] находятся за 12-19 итераций, а в [10] - за 2 - 3 . П о э т о м у эффективность решения задачи о рас­
паде произвольного разрыва остается актуальной проблемой.
В данной работе рассмотрен новый численный метод решения задачи о распаде произвольно­
го разрыва в средах с нормальными УРС (см. [1]) методом Ньютона. Метод является обобщени­
ем приближенного метода (см. [7]), основанного на подходах работ [8], [10], и заключается в сле­
дующем. На каждой итерации исходные У Р С локально аппроксимируются двучленными, для ко­
торых решение задачи известно (см. [3], [8], [10]) и в ы р а ж а е т с я через элементарные функции.
По этим решениям определяются новые состояния сред в точках на исходных ветвях (и, /^-диа­
грамм. Задача о распаде разрыва снова решается т а к ж е , как и на предыдущих итерациях, с на­
чальными данными, которые соответствуют найденным новым состояниям. Возможности под-
хода демонстрируются на примерах решения задач в средах, подчиняющихся У Р С М и - Г р ю н а й ­
зена. Для этого УРС выписаны в явном виде алгебраическое уравнение изэнтропы и н е к о т о р ы е
точные решения для конфигураций с BP. Результаты численных расчетов модельных задач о
распаде произвольного разрыва показали эффективность рассмотренного метода, согласуются
с точными решениями и решениями, полученными по алгоритму из [7]. Метод реализован в про­
грамме, написанной на языке программирования C++ в среде WINDOWS для решения задач на
персональном компьютере в интерактивном режиме.
2. З А Д А Ч А О Р А С П А Д Е П Р О И З В О Л Ь Н О Г О Р А З Р Ы В А
Пусть в начальный момент времени t = 0 среда в пространстве разделена на две части - левую
и правую - плоскостью х = 0. Среда слева от этой плоскости имеет состояние р и р справа состояние р , и р . Среды подчиняются УРС, записанным в форме p = р\(р, E), р = /? (р, £ ) , ко­
торые удовлетворяют условиям Бете-Вейля [1], [12], являются выпуклыми и называются нор­
мальными. Здесь р, и,р,Еэто удельная плотность, скорость, давление и внутренняя энергия
соответственно. Среды разделены перегородкой, которая в момент времени t = 0 мгновенно уби­
рается. После этого среды приходят в соприкосновение и начинают взаимодействовать между
собой. Требуется определить состояние сред в момент t > 0.
В такой постановке задача известна как задача Римана о распаде произвольного разрыва.
Предполагаем, что для давлений выполнено условие р < р . Поэтому можно рассмотреть только
основные конфигурации - это конфигурации из двух УВ, двух BP, одной УВ и BP. Решение за­
дачи наглядно представляется в плоскости состояний сред и находится в точке пересечения вет­
вей (и, /?)-диаграмм (см. [2]). Поскольку задача о распаде разрыва хорошо известна, то м ы не бу­
дем приводить вывод основных уравнений, а выпишем сразу необходимые из них в следующем
общем виде:
0
ь
2
ъ
ъ
2
ь
2
х
л
^ - Р *
U - Ll +
Р
n
U - И2
о
-
{
(J,
2
РГ - Р,
о
P
~PÏ
J
= U,
v
Ü;
=
<
(D
i
1
P<Pf
Здесь / = 1, 2, £/, Р - скорость и давление на контактном разрыве (KP) в зонах постоянного тече­
ния, V = 1/р - удельный объем, величины p * , Е* , Vf - плотность, удельная внутренняя энергия
и удельный объем за фронтом УВ соответственно. Система уравнений (1) в зависимости от реа­
лизовавшейся конфигурации дополняется уравнениями адиабат Гюгонио
£ f - £ , - 0 . 5 ( V , . - V * ) ( / > , + P*) = 0,
P* =
(2)
p.(V?,E*),
и/или адиабат Пуассона
d
-b = -р(р,Е),
j = 1,2.
(3)
dv
Из (1)-(3) следует, что в конфигурации из двух У В требуется решить систему нелинейных ал­
гебраических уравнений, из двух BP - систему интегральных уравнений, из У В и BP - смешан­
ную систему уравнений. Для определения конфигурации проводим анализ поведения функции
F(P) (см. [3], [8]):
1) если и - и > U = F(p ), то Р > р > р и реализуется конфигурация из двух УВ;
J
х
2
2
2
2
2) если U у = F(p ) <и -и <
{
{
х
U , то р < Р <р и реализуется конфигурация из УВ и BP;
2
2
х
2
3) если U = F(p ) < щ < и < U то Р < р < р и реализуется конфигурация из двух BP;
0
0
2
h
{
2
4) если щ - и < (/ , то реализуется область вакуума со значениями р = 0, с = 0.
2
0
В (1) вычитая второе уравнение из первого, получаем уравнение для давления в виде
F(P)^f (P)
{
+ f (P)
2
= щ-и .
2
(4)
1104
МОИСЕЕВ, МУХАМАДИЕВА
Здесь
p-Pi
а
a-,
Л(Р) =
P f - P,'
(5)
Pi
fedp,
P<p,
И з уравнения (4) определяем давление на контактном разрыве (KP), а затем вычисляем осталь­
ные величины.
3. М Е Т О Д Н Ь Ю Т О Н А ДЛЯ РЕШЕНИЯ З А Д А Ч И О Р А С П А Д Е П Р О И З В О Л Ь Н О Г О
Р А З Р Ы В А В СРЕДАХ С П Р О С Т Ы М И У Р А В Н Е Н И Я М И С О С Т О Я Н И Я
Предположим, ч т о из уравнений (2), (3) можно получить уравнения адиабат Гюгонио и Пуас­
сона в аналитической форме. УРС с такими свойствами в дальнейшем будем называть простыми.
О б ы ч н а я процедура решения нелинейной системы уравнений методом Ньютона, к а к правило,
состоит в линеаризации всех уравнений и последующей организации итерационного процесса по
нахождению решений системы линейных уравнений. М ы рассмотрим подход к р е ш е н и ю систем
уравнений (1), который отличается от подходов, в [4]-[6]. Для решения этих систем уравнений
оказывается достаточно линеаризации только исходных УРС. Алгоритм решения в э т о м случае
строится следующим образом. Проводим локальную линеаризацию исходных УРС. Вследствие
такой линеаризации среды будут подчиняться двучленным УРС вида
р = (у-1)р£ + Со(р-ро).
(6)
Здесь у - показатель адиабаты, р , с - параметры вещества, вычисленные в результате линеари­
зации. В этом случае уравнение (4) упрощается и его решение находится по алгоритму из [8] с уче­
том модификации в [10] методом обратной парабалической интерполяции. Подставив найденные
значения Р, p f , р * давления и плотностей в уравнение (4), проверяем выполнение условия
0
0
\F(p)-(u -u )\<z\u -u \,
l
2
l
2
0 < 8 < ^ 1.
(7)
Если условие выполняется, то считаем, что приближенное решение находится в окрестности
точного решения уравнения (4), заканчиваем итерационный процесс и переходим к счету осталь­
ных величин. Если условие (7) не выполняется, т о поступаем следующим образом. И з получен­
ного приближенного решения берем плотности pf , р * в средах слева и справа от KP, вычисляем
давления pf , pf из исходных уравнений (2), (3) адиабат Гюгонио или Пуассона и скорости uf ,
uf , из уравнений (1). Очевидно, что величины ( p f , pf , uf ) , (uf , pf , p * ), найденные на этом
шаге, вычисляются без итераций и определяют новые состояния сред в точках, л е ж а щ и х на ис­
ходных ветвях (и, /?)-диаграмм. Далее повторяем описанный выше процесс получения очередно­
го приближенного решения, взяв в качестве начальных данных величины ( p f , pf , uf ) , (uf , pf ,
p * ). Процесс повторяется до тех пор, пока неравенство (7) не будет выполнено.
Геометрическая интерпретация (см. фиг. 1, 2) описанного алгоритма означает, ч т о в плоско­
сти переменных (р, р) кривую УРС аппроксимируем (заменяем) касательной, а адиабаты Гюго­
нио или Пуассона и ветви состояний исходных У Р С на (и, /?)-диаграмме аппроксимируем соответ­
ствующими адиабатами и ветвями для двучленных УРС (6).
Поскольку имеется взаимно однозначное соответствие между УРС, адиабатами и ветвями, т о
аппроксимирующие адиабаты и ветви по построению будут касаться исходных адиабат и ветвей
соответственно. Следовательно, описанный метод решения уравнения (4) соответствует методу
Ньютона, в котором вместо касательных используются адиабаты, соответствующие У Р С (6). Т о
ж е самое относится и к ветвям (и, /?)-диаграмм в плоскости состояний. Поэтому скорость сходимо­
сти к решению должна быть выше, чем в обычном методе с использованием только касательных.
Приближенные решения можно находить также и на основе модифицированного уравнения
из [10], которое отличается от уравнения (6) наличием свободного безразмерного параметра пе-
Фиг. 1. Две BP: (а) - плоскость (р, р) и (б) - плоскость (w, р).
Фиг. 2. Две УВ в плоскости (и,р).
ред вторым слагаемым. Это позволит сократить число итераций на начальной стадии поиска ре­
шения. Однако на заключительной стадии поиска лучше переходить на использование уравне­
ния (6), поскольку алгоритм выбора значений свободных параметров не обеспечивает квадра­
тичной скорости сходимости.
4. У Р А В Н Е Н И Е С О С Т О Я Н И Я М И - Г Р Ю Н А И З Е Н А
Возможности описанного подхода к решению задачи о распаде разрыва рассмотрим на при­
мере сред, подчиняющихся уравнению состояния Ми-Грюнайзена
ET
— E — Еу,
т
x
= —I
Еу
x
:— +
ил ц - 1
5 =
5
ï-, к = 1,2,
которое запишем в виде
p = ( _i)p£
Y
+
j
^£i|( -i)(ô-l)
Y
+
IJ(l-8^)|.
(8)
Здесь p с - параметры вещества, у, р - параметры уравнения состояния, которые будем считать
постоянными. У Р С хорошо известен и широко применяется для описания поведения многих ве­
ществ. Мы не будем подробно останавливаться на его свойствах, с которыми м о ж н о познако­
миться, например, в [13], [14], а отметим следующее. Скорость звука вычисляется по формуле
h
к
1/2
С
' Р + Рк
у-—— + Р - У ^ - !
Рк =
Рк к
Уравнение адиабаты Гюгонио в виде явной зависимости Р = Р(р) просто получить из у р а в н е ­
ний (2) и (8). Следовательно, решение задачи о распаде разрыва по описанному выше алгоритму
в конфигурациях с УВ не представляет каких-либо трудностей не только для У Р С М и - Г р ю н а й ­
зена, но и для У Р С такого типа. Сложнее обстоит вопрос с решением задачи в конфигурациях с
BP, поскольку необходимо вычислять вдоль изэнтропы интеграл Римана (5), который не всегда
можно выразить через элементарные функции. Более того, для многих нормальных У Р С часто
неизвестно и само уравнение изэнтропы в аналитическом виде. Поэтому нам представляется
целесообразным выписать это уравнение для УРС Ми-Грюнайзена в удобном для д а л ь н е й ш е ­
го исследования виде. Подставив в уравнение (3) вместо Е выражение из (8), получим следую­
щее дифференциальное уравнение:
V
aW
=
M
" ^ - Y ^ + (Y-^)P*v?v" .
(9)
Проинтегрировав уравнение (9) методом вариации произвольной постоянной и выбрав эту по­
стоянную, согласованную с двучленным уравнением состояния, получим уравнение изэнтропы в
явном аналитическом виде:
l
p =
7
y
>i
i
-o(S)p -p -p b (l-6 -' ).
k
k
(10)
Здесь ü(S) - энтропийная функция, которую можно вычислять по формуле
n
о(5) = Y
^ +Y
•
P
Р
Скорость звука вдоль изэнтропы можно вычислять по формулам
(П)
7
1
1
/о^р'-'-^-'-цб' - ) = y ^
m-y)Ejf.
(12)
А/
Рк
v
P
P
П о л у ч е н н ы е уравнения изэнтропы и энтропийной функции для УРС М и - Г р ю н а й з е н а п о з в о ­
л я ю т существенно упростить решение задачи о распаде разрыва и получить частные в ы р а ж е ­
ния для параметров у, р, при которых интеграл Римана записывается в квадратурах. Если у = р ,
т о ф о р м у л ы ( 10)—( 12) переходят в формулы для У Р С (6), который о б ы ч н о н а з ы в а ю т согласо­
ванным.
с =
+
5. Т О Ч Н Ы Е Р Е Ш Е Н И Я
Рассмотрим значения параметров у, р, при которых интеграл Римана (5) можно записать в
квадратурах. Преобразовав интеграл (5) к виду
/(p)=f^p
J p
p
(13)
Pf
на отрезке [ р * , р,] и подставив под знаком интеграла (13) выражение для скорости звука из (12),
получим
Р,
8,{y i)n
m
b - (a
Pf
y
l
Здесь а = c(S)p ~
k
+
b^) db.
ôf
2
- ус^/р, b - (\ip )/p
k
k
= q , 8
г
= р//р*. Согласно теореме Чебышёва, интеграл
ti
jV" (а + bx ydx может быть выражен через элементарные функции, если одно из чисел (т +
р + (т+ 1)Аг, р целое. Здесь р,т,прациональные числа. В нашем случае m - 0 . 5 ( у - 3), п - р - у,
р = 0.5. Следовательно, если параметры p, у связаны одним из следующих соотношений:
(2*+1)у-1
2Jfcy-l
.
,
0
то интегралы будут табличными (см. [15], [16]) и будут выражаться через элементарные функ­
ции. Решения задачи о распаде разрыва с такими параметрами у, Ц будем называть точными и
будем их использовать для контроля точности и эффективности описанного метода. Так, если у
= 3 + 2к,к = 0, 1,
и р = у + 1, т о получим интеграл вида
к
lL
/ = Jх [а +
bx] dx.
ôf
Для к = 1 параметры p, у из первого условия связаны соотношением р = у + 0 . 5 ( у - 1). В этом слу­
чае скорости uf , uf можно вычислять из следующих точных уравнений:
и*
К ^ + ад^-^ + М Г Л
= м^ + ^
/ = 1,2.
6. П Р И Б Л И Ж Е Н Н Ы Й М Е Т О Д В Ы Ч И С Л Е Н И Я ИНТЕГРАЛОВ Р И М А Н А
В общем случае значение интеграла Римана (5) находят каким-либо численным методом (тра­
пеций, прямоугольников, Симпсона и т.п.) либо аппроксимируют подынтегральную функцию
другой функцией, интеграл от которой выражается в квадратурах. Рассмотрим второй подход
приближенного вычисления интегралов. Аппроксимируем исходную изэнтропу (10) кубической
параболой вида
3
2
fe p + b i p + b p + fe = р
0
2
3
(14)
9
которая совпадает с изэнтропой на концах интервала интегрирования вместе со своими первыми
производными. К о э ф ф и ц и е н т ы кубической параболы находятся из решения соответствующей
системы уравнений и имеют достаточно простой вид:
2
Ь =
P2
p
с1
( ~i
2
(p 2
P
l
с
\+ г\
) Vp2-Pi
2
У
2
1
7
Ъ =
(p 2
2
2
J
P l
c
+
\
2
C
Рг-Р\\
2
+ 3 ( р + Pi)!
V 3
Pl^2 + Р 2 ^ 1
у
b
2
Г
,
2
)L
=cf-3ft pf-2^^!,
Ъ =
0
p -b p\-b p]-b p .
ъ
l
Ç)
x
2
{
Построив параболу, можем получить выражение для скорости звука
2
с = 7 з Ь р + 2/?!р + ô ,
0
2
которое будем использовать в (13) вместо (12). Тогда интеграл (13) преобразуется к виду
„
P
.
t
f
73/7 p
0
2
+ 2fe p + fe
1
/(р)«
2
dp
p
J
pf
и точно выражается через элементарные функции.
Покажем корректность рассмотренной аппроксимации. Целесообразность аппроксимации
кубическим полиномом обсуждается в [9] ввиду того, что полином может быть не монотонным
на отрезке, на котором осуществляется аппроксимация. Покажем, что в нашем случае кубиче­
ский полином (14) будет монотонным на интервале [pf , p j . Напомним, что мы рассматриваем
выпуклые уравнения состояния.
Проведем хорду через точки ( p f , pf ), (р , /? ), которые лежат на изэнтропе, и восстановим
перпендикуляр из середины этой хорды. Перейдем от исходной системы координат к системе ко2
.
i
2
к
ординат, в которой координатные оси (pj , р ) совпадают с хордой и перпендикуляром соответ­
ственно. В этой системе координат построенная кубическая парабола будет представляться в ви­
де полинома третьей степени, к о т о р ы й обращается в ноль в двух точках, соответствующих точх
кам ( р р ) , (р , р ) в исходной системе координат, и имеет производные разных знаков в этих
точках. Очевидно, что эти точки соответствуют действительным корням полинома в новой си­
стеме координат. Поскольку полином третьей степени может иметь один и два мнимых корня
или три действительных, то наш случай соответствует трем действительным корням. Следова­
тельно, поведение полинома внутри интервала [ p f , p j зависит от положения третьего корня.
Если корень лежит между двумя отмеченными корнями, то это соответствует наличию макси­
мума и минимума внутри интервала. В этом случае в силу построения полинома знак производ­
ной на одном из концов интервала будет совпадать со знаком заданной производной, а на другом
конце будет противоположен знаку исходной производной. Это противоречит условию постро­
ения полинома. Поэтому третий к о р е н ь будет лежать за пределами рассматриваемого интерва­
ла. Таким образом, доказано, что кубическая парабола (14) будет выпуклой на интервале [pf ,
р,]. Следовательно, производная будет непрерывно изменяться от первой точки ко второй без
смены знака.
ь
х
2
2
7. Р Е Ш Е Н И Е З А Д А Ч И О Р А С П А Д Е Р А З Р Ы В А В СРЕДАХ
СО С Л О Ж Н Ы М И У Р А В Н Е Н И Я М И СОСТОЯНИЯ
Рассмотрим адаптацию описанного алгоритма к решению задачи о распаде разрыва в средах,
подчиняющихся УРС, для которых нельзя выписать адиабаты Гюгонио или Пуассона в явном
виде, но которые можно аппроксимировать двучленными УРС. Такие УРС будем называть
сложными. Пусть реализовалась конфигурация с УВ. Тогда, решив задачу для двучленных У Р С
(6), давление находим из решения системы (2) и далее по описанному алгоритму. В случае кон­
фигураций с BP поступаем следующим образом. Решив задачу для двучленных УРС (6), по най­
денным значениям плотности и внутренней энергии вычисляем давление pf и скорость звука с*
из исходного УРС. Затем находим п а р а м е т р ы кубической параболы (14), которая проходит че­
рез точки (p Pi) и (pf , pf ) в плоскости переменных (р, р) и имеет касательные, наклоны кото­
рых равны квадратам скоростей звука c с* в этих точках соответственно. Предполагаем, что
эта парабола аппроксимирует изэнтропу в окрестности точки (р,-,/?;), в которой давления P , jP ,
вычисленные из исходного УРС и уравнения кубической параболы соответственно, удовлетво­
ряют условию
h
h
es
|^е.-^ср|<ер,-,
cp
0 < е « 1 .
Решение модельных задач показало, ч т о такие окрестности существуют. Взяв любое значение
плотности pf из этой окрестности и вычислив давление pf , продолжим решение задачи о рас­
паде разрыва согласно описанному в ы ш е алгоритму. Естественно, что в этом случае решение бу­
дет приближенным.
Замечание. Описанный подход к решению задачи о распаде разрыва позволяет наметить пути решения
и в средах, подчиняющихся ненормальным УРС. Среди таких УРС особый интерес представляют УРС, для
которых вторая производная
2
д р(У,
S)
2
дУ
является знакопеременной величиной. Если предположить, что условия существования разрыва выполне­
ны для таких УРС (см. [1]), то решение может быть следующим. Развиваем область определения УРС на
подобласти, в которых УРС будет выпуклым, т.е. вторая производная знакопостоянна. Тогда решение за­
дачи о распаде разрыва, по всей видимости, может быть получено по описанному методу путем последо­
вательного прохождения этих подобластей. Однако детальное исследование решения выходит за рамки
данной работы, поэтому мы ограничимся только таким замечанием.
8. Р Е З У Л Ь Т А Т Ы Р А С Ч Е Т О В Т Е С Т О В Ы Х З А Д А Ч
Верификация предложенного метода проводилась на модельных задачах путем сравнения с
точными решениями и решениями, полученными по программе из [7]. Анализ результатов рас­
четов показывает, что численные решения совпадают с точными решениями и решениями, по­
лученными по методу из [7], до 8—10 знаков. Численные решения задач без учета аналитического
уравнения изэнтропы (10) также согласуются с точными решениями и совпадают до двух знаков
h
0.2
0.1
0.05
0.01
0.005
0.001
Давление на KP,
метод Симпсона
Скорость на KP,
метод Симпсона
22.100332436991
1.10581609692354
-
Давление на KP,
метод кубической
интерполяции
Скорость на KP,
метод кубической
интерполяции
22.100332436991
22.10032
22.100330
22.100332432
22.10033243698
22.100332436993
22.10033243697
1.10581609692354
1.1058162
1.1058162
1.10581609695
1.1058160969232
1.1058160962329
1.10581609692357
-
1.1054
1.1057
1.105812
1.105815
1.10581606
22.11
22.103
22.1004
22.10036
22.100333
после запятой. Здесь представлены результаты расчетов задачи, в которой реализовалась кон­
фигурация из двух BP. Заданы следующие состояния сред:
слева от разрыва
справа
pi=4
р =6
2
и =0
м =6
{
2
Р\ = 46
р = 221.25
Е = 3.83(3)
Е = 14.125
р =2
р =3
{
2
0
2
0
с =2
с =3
{)
{)
у=3
у=3
р =4
р =4
Среды подчиняются УРС Ми-Грюнайзена. Задача имеет точное решение. На решениях этой
задачи проведено сравнение эффективности метода в зависимости от точности вычисления ин­
тегралов Римана методами Симпсона и кубической интерполяции изэнтропы. Число интерва­
лов, на которое разбивался отрезок интегрирования, находилось по формуле N = [\p - Pn\/h] + 1,
/ = 1,2. Здесь h - шаг интегрирования, р , p - плотности в средах до и после решения задачи о
распаде разрыва с двучленными УРС, [ ] - целая часть числа. В таблице в зависимости от шагов
интегрирования приведены результаты расчетов давлений, скоростей на KP при вычислении ин­
тегралов Римана методами Симпсона и кубической интерполяции изэнтропы. В первой строке
приведены точные решения.
i2
п
i2
Из анализа результатов, представленных в таблице, следует, ч т о точность вычисления инте­
гралов и скорость сходимости выше при кубической интерполяции изэнтропы, чем при вычис­
лении интегралов методом Симпсона.
На фиг. 3 представлены графики (1.3) и (2.4) исходных и аппроксимирующих изэнтроп соот­
ветственно, а также результаты расчетов той ж е задачи, но с измененной скоростью и = 13. Рас­
считанные значения давлений в точном (точки 5 на графиках) и приближенном (точки 6) реше­
ниях равны 0.69475 и -43.197 соответственно. Видно, что в этой задаче приближенные решения
дают совершенно неправильный результат.
2
250
200
150
100
50
0
-50
-100
0
1
2
3
Фиг. 3.
4
5
6
К сожалению, авторам не удалось найти данных по решению задачи о распаде разрыва мето­
дами в [4], [6], чтобы можно было сравнить с этими методами э ф ф е к т и в н о с т ь представленного
метода. В [5] приведена информация о том, что для решения задач в средах с У Р С Ми-Грюнай­
зена требуется до 25 итераций. Однако не приводятся постановки таких задач.
9. З А К Л Ю Ч Е Н И Е
Рассмотрен подход к решению задачи о распаде произвольного р а з р ы в а в средах, подчиняю­
щихся нормальным УРС, методом Ньютона на основе решений простейших задач с двучленны­
ми уравнениями состояния, которые локально аппроксимируют исходные У Р С . Для сред, подчи­
няющихся УРС Ми-Грюнайзена, выписаны в явном виде алгебраическое уравнение изэнтропы
и некоторые точные решения для конфигураций с BP. Предложен алгоритм приближенного вы­
числения интегралов Римана, обладающий более высокой точностью и скоростью сходимости,
чем алгоритм в методе Симпсона.
СПИСОК ЛИТЕРАТУРЫ
1. Рождественский БЛ., Яненко H.H. Системы квазилинейных уравнений. М.: Наука, 1968.
2. Овсянников Л.В. Лекции по основам газовой динамики. М.: Наука, 1981.
3. Годунов С.К., Забродин A.B., Иванов М.Я. и др. Численное решение многомерных задач газовой ди­
намики. М.: Наука, 1976.
4. Аладыкин Г.Б., Годунов С.К., Киреева ИЛ., ПлинерЛ.А. Решение одномерных задач газовой динами­
ки в подвижных сетках. М.: Наука, 1970.
5. Шустов Ю.М. Расчет распада разрыва для произвольных уравнений состояния // Числ. методы механ.
сплошной среды. Новосибирск: ВЦ СО АН СССР, 1978. Т. 9. № 4. С. 131-138.
6. Куропатенко В.Ф., Коваленко Г.В., Кузнецов В.И. и др. Комплекс программ "ВОЛНА" и неоднород­
ный разностный метод для расчета неустановившихся движений сжимаемых сплошных сред // ВАНТ.
Сер. Методики и программы числ. решения задач матем. физ. 1989. Вып. 2. С. 19-25.
7. Куликовский А.Г., Погорелое Н.В., Семёнов А.Ю. Математические вопросы численного решения ги­
перболических систем уравнений. М.: Физматлит, 2001.
8. Прокопов Г.П. Расчет распада разрыва для пористых сред и сплошных сред с двучленными уравнени­
ями состояния // ВАНТ. Сер. Методики и программы числ. решения задач матем. физ. 1982. Вып. 2(10).
С. 32-40.
9. Прокопов Г.П. О приближенных реализациях метода Годунова: Препринт № 15. М.: ИПМ матем.
РАН, 2007.
10. Кобзева Т.А., Моисеев Н.Я. Метод неопределенных коэффициентов для решения задачи о распаде
произвольного разрыва // ВАНТ. Сер. Методики и программы числ. решения задач матем. физ. 2003.
№ 1. С. 3-9.
11. КоллатцЛ. Функциональный анализ и вычислительная математика. М.: Мир, 1969.
12. Weyl H. Shock waves in arbitrary fluids // Communs Pure and Appl. Math. 1949. № 2. P. 103-122.
13. Забабахин Е.И. Некоторые вопросы газодинамики взрыва. Снежинск: РФЯЦ-ВНИИТФ, 1997.
14. Кобылкин И.Ф., Селиванов В.В., Соловьев B.C., Сысоев H.H. Ударные и детонационные волны. Ме­
тоды исследования. М.: Физматлит, 2004.
15. Двайт Г.Б. Таблицы интегралов и другие математические формулы. М.: Наука, 1978.
16. Бронштейн H.H., Семендяев К.А. Справочник по математике для инженеров и учащихся втузов. М.:
Наука, 1986.
