О влиянии свойств разностных схем на математическое моделирование динамических процессов Текст научной статьи по специальности «Физика»

Вычислительные технологии Том 11, часть 2, Специальный выпуск, 2006

О ВЛИЯНИИ СВОЙСТВ РАЗНОСТНЫХ СХЕМ НА МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ДИНАМИЧЕСКИХ ПРОЦЕССОВ*

В.Ф. Куропатенко, И. А. Доровских, И. Р. Макеева Российский федеральный ядерный центр — ВНИИ технической физики им. акад. Е. И. Забабахина,

Снежинск, Россия e-mail: v. f . kuropatenko@vniitf . ru

The paper analyses the results of numerical experiments using several classical techniques and a new difference scheme that optimally combines the monotony and small discontinuity distraction properties.

Введение

В результате импульсных воздействий на вещество (электровзрыв фольги, детонация взрывчатого вещества, удар пластиной, поток электронов или ионов, рентгеновское, лазерное или другие виды излучений), как правило, возникают нестационарные ударные волны сжатия и разрежения. Картина взаимодействия этих волн может быть весьма сложной. Точность математического моделирования ударно-волновых процессов зависит как от применяемой физической модели, так и от свойств разностных уравнений (немонотонность, дистракция разрывов, производство энтропии, энтропийные следы и др.). Поскольку эти свойства не зависят от физической модели, их изучение можно проводить для такой модели, которая позволяет строить аналитические решения конкретных задач. При таком подходе сравнение результатов расчетов с аналитическим решением является эффективным способом оценки точности.

Основной целью данной работы является анализ точности расчетов разрушения вещества, По мнению авторов, результаты расчетов разрушения наиболее чувствительны к особенностям разностных уравнений. Именно поэтому в работе анализируются результаты расчетов по широко применяемым методам Неймана — Рихтмайера [1], Годунова [2], Лакса — Вендрофа [3] и по новой разностной схеме, оптимальным образом сочетающей свойства монотонности с малой дистракцией разрывов. Показано, что погрешности разностного метода существенно влияют на точность определения массы отколовшегося слоя.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 04-01-00050).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

1. Разностная схема

Рассмотрим разностную схему для расчета непрерывных и разрывных решений уравнений идеальной сплошной среды, в которой минимизация амплитуды колебаний и дис-тракцпп достигается путем выбора специальной формы разностных уравнений и применением условий Гни опт) — Ренкина для описания диссипации энергии в зоне ударного слоя, которым заменяется сильный разрыв при разностном счете. Уравнения траектории, неразрывности и движения аппроксимируются следующими разностными уравнениями:

С+1 = с + ти*,

(гп+1\а _ (гп+1\а т/га+1 _ т гп V г+1 /_V г /

4+0.5 - 4+0.5 (г«+1)а _ ' ^

тУп

ТТП+1 _ ТТП+1 _ »+0.5 ( и* _ р*\

иг+0.5 ~ иг+0.5 „ _ „ \Гг+1 Гг)' 'г+1 'г

Здесь У — удельный объем; и — скорость; Р — давление; г — пространственная координата; т — шаг сетки по времени а = 1 в случае плоской симметрии, а = 2 в случае цилиндрической симметрии, а = 3 в случае сферической симметрии; г — индекс узла сетки по г; п — индекс узла сетки по

В случае и*+0 5 — и*_0 5 < 0 предполагается, что из точки тг распространяется сеточная ударная волна, величины перед фронтом которой определяются следующим образом:

(Р, У, Е)_ = (Р, У, Е)г+0.5, если Р^ > РП+0.5,

(Р,У,Е)_ = (Р, У, Е)г_0.5, если Р^ < Рг+о.5-

Скорость за фронтом ударной волны полагается равной скорости на противоположном конце промежутка [тг_0.5, тг+0.5]:

и+ = Ц_0.5, если Р™_0.5 > РГ+0.5,

и+ = Ц+0.5, если Р™_0.5 < РГ+0.5. Из условий на сильном разрыве находится значение давления за фронтом ударной волны Р+. Значения Р+, и+ принимаются [4] в качестве вспомогательных значений

Р* = Р+, и* = и+.

Удельная внутренняя энергия Е™^ на ударных волнах находится из уравнения

Еп+1 — Е п 1 Уп+1 — Уп

Ег+0.5 Ег+0.5 _ 1 (. и*\ 4+0.5 4+0.5

; - _ 2{г+1 + г' ;

1 ип+1 — и п

+ о.б - и:+1 - и:) .

иг* , Рг*

ными уравнениями

т

2К'

_ (гп )а_1ип _ (гп )а_1ип ^ ^

и*— и™ {Р™+о.ь РР- о.б))

Г г/ ^ К

< =

(а?+0.5 )2 + (°П—0.5 )2

где а — жесткость вещества;

1,

1

при а = 1,

Ф-

^ \ „1+0.5 + Г 1-0.5,

(гп )2 | гп гп | (гп )2 ('г+0.5) + ' ¿+0.5' г—0.5 + ('г-0.5) ,

^ = ^ п при а = 2, (3)

4/3п при а = 3.

На волне разрежения новое значение энергии Е™^ определяется интегрированием уравнения изэнтропы [5, 6]

уп+1

Еп+1 - еп + J р(V, Е)дУ = 0.

V п

Значение давления на (п + 1)-м шаге определяется по уравнению состояния.

Для минимизации немонотонности в окрестности слабых разрывов уравнения для вспомогательных величин корректируются с учетом членов третьего порядка малости, В случае ударной волны и в случае волны разрежения описанная разностная схема

та

устойчива при К = < 1- Эффективная дистракция [7], определенная для стационарной

п

ударной волны произвольной амплитуды в идеальном газе, равна

В

э = 2(1 -К) I л/Уо + л/Уг 7+1 \у/У0-у/\

2. Тестирование разностных схем

Для тестирования разностных методов были рассчитаны задачи, имеющие аналитические решения.

Задача 1. Автомодельная сферическая сходящаяся ударная волна. Начальные данные и граничные условия выбираются следующим образом: в области 0 < г < 1 при £ = 0 находится покоящийся холодный идеальный газ с 7 = 5/3, Е0 = 0,р0 = 1,Р0 = 0, и0 = 0. Условие в центре симметрии I] = 0 при г = 0, Условие на правой границе, которая является траекторией макрочастицы, берется из автомодельного решения Л,Д. Ландау и К.П. Станюковича [8] и аппроксимируется табличной зависимостью ] от

При £ = 0 начальные данные задавались па равномерной сетке по радиусу с числом точек N. На рис, 1 приведены результаты расчетов этой задачи, аналитического решения и расчетов по программе "Волна'Г[9] при N = 80.

£ ипг £ ипг £ ЦЦд р

0 -1 0.328273 -1.04392 0.568207 -0.876229

0.02 -1.004474 0.36 -1.0395 0.613521 -0.77159

0.102828 -1.02181 0.39144 -1.03309 0.632854 -0.71116

0.16 -1.0323 0.42 -1.0227 0.650751 -0.64335

0.20701 -1.0391 0.438247 -1.01479 0.667652 -0.56491

0.241516 -1.04277 0.475136 -0.991499 0.675862 -0.51971

0.286642 -1.04513 0.50547 -0.964469 0.683804 -0.46939

0.31 -1.045 0.540474 -0.92164

Р

-0.8

-1

-1.2

-1.4

-1.6

-1.:

и

4 Г " "1 : 2

ц

г

0.34 0.36 0.38 0.4 0.42 0.44 0.46

4

г

3

Рис. 1. Задача 1. Зависимости давления Р (г) и скорости и (г) на момент Ь = 0.4: кривые 1 — фрагмент аналитического решения, кривые 2 — результат данной работы, соответственно кривые 3 и 4 — расчеты по программе "Волна" с выделением точного фронта и без него [9].

Задача 2. Коллапс сферического слоя идеального газа [10]. В момент Ь = 0 фокусируется слой холодного идеального газа с 7 = 5/3, Начальные данныв! в области 0 < г < 1 и0 = — 1, р0 = 1, Е0 = 0, Р0 = 0. Задача имеет особенность в точке Ь = 0, г = 0, где производная дУ/дЬ бесконечна. Это следует из того, что при и = const и ¿Ы = г2рйг закон

*2С/° п тт п

сохранения массы принимает вид -—— =-и при г —> 0 и и о < --> — сю, В конечной

дЬ р0г дЬ

окрестности точки Ь = 0, г = 0 производиая дУ/дЬ также не мала. Это ставит любой разностный метод в очень жесткие рамки, поскольку шаги сетки нельзя взять равными нулю Аг = 0 и ДЬ = 0, Тем не менее расчеты проведены, их результаты представлены на рис, 2, В начальный момент сетка выбиралась равномерно по г, число точек N = 200,

Р

25

20

15

10

/2

*

3 /

80 70 60 50 40 30 20 10 0

: 1

^_^

V "г 2 /.....

Р

4

0 0.1 0.2 0.3

0.4

0.1 0.2 0.3

Р

5

Г

0

0

Рис. 2. Задача 2. Зависимости давления Р(г) и плотности р(г) на момент Ь = 0.6: кривые 1 — аналитическое решение, кривые 2 — результат данной работы, кривые 3 — расчет по программе "Волна" без выделения разрыва [9], кривые 4 — результат из [10].

Задача 3. О точечном сферическом взрыве в идеальном газе ("седовский взрыв" [11]). Начальные данные: при £ = 0 в области 0 < г < 0.1 заданы величины р =1, Е = 107, Р = 0.4 ■ 107, и = 0. В области 0.1 < г < 1.4 р = 1, Е = 0, Р = 0, и = 0. Уравнение состояния: идеальный газ с 7 = 1.4. Граничные условия: при £ > 0, г = 0, и = 0, г = 1.4, Р = 0 0 < г < 0.1

точек N = 5 Сетка в области 0.1 < г < 1.4 равномерная по радиусу с числом точек N1 = 200.

На рис. 3 приведены результаты расчета в момент выхода фронта ударной волны в точном решении на радиус г = 1 (£ = 3.5). Результаты приведены в безразмерном виде. Безразмерные функции Р(г),р(г) получаются путем их деления соответственно на значения Р+, р+ в точном решении на фронте ударной волны в момент £ = 3.5. На рис. 4 приводится фрагмент решения в момент £ = 3.5.

Р

0.8

0.6

0.4

0.2

/ - 2 У

-^-А

0.5

1.5

0.8

0.6

0.4

0.2

1^г

1 2 Ж^.........

у

0.5

1.5

Р

1

1

г

г

0

0

0

1

0

1

Рис. 3. Задача 3. Зависимость давления Р(г) и плотности р(г) на момент Ь = 3.5: кривые 1 аналитическое решение, кривые 2 — результат данной работы.

Рис. 4. Задача 3. Зависимость давления Р (г) и скорости и (г) на момент Ь = 3.5. Фрагмент решения в окрестности фронта ударной волны: кривые 1 — аналитическое решение, кривые 2 — результат данной работы.

Рис. 5. Задача 4. Зависимости давления Р и скорости и от координаты х за фронтом ударной волны. Фрагмент решения: кривые 1 — точное решение, кривые 2 — метод Неймана — Рихтмай-ера [1], кривые 3 — метод Лакса — Вендрофа [3], х х х — результаты данной работы.

Задача 4. Распространение стационарной ударной волны. При £ = 0 в области с координатами [0, 14] находится идеальный газ с уравнением состояния Р = (7 - 1)рЕ и следующими начальными параметрами: 7 = 4/3, р0 = 1, Е0 = 0, Р0 = 0, и0 = 0. На левой

и=3

следующие параметры за фронтом ударной волны: р1 = 7, Е1 = 4.5, Р1 = 10.5, и = 3, скорость распространения ударной волны Ш = 3.5. Сетка по пространству при £ = 0 задана равномерной с Дж = 0.2 (число то чек N = 70).

На рис. 5 приведены фрагменты решения в малой окрестности за фронтом ударной волны, полученного по трем разностным схемам, и результаты аналитического решения. Видно, что предложенная разностная схема дает значительно меньшие осцилляции за фронтом ударной волны.

Задача 5. Распространение волны разрежения по постоянному фону. В области с координатами [0,14] находится покоящийся идеальный газ с начальными параметрами: 7 = 2,р0 = 4.5, Е0 = 1.125, Р0 = 5.0625, и0 = 0. На левой границе системы задана постоянная скорость и = —1. Таким образом, в точке ж = 0, £ = 0 находится произвольный разрыв. При £ > 0 он распадается и в газ распространяется волна разрежения,

Р

Р

1.8

1.6

1.4

1.2

-0.2 -0.1 0 0.1 0.2 0.3

0.8

А, Ф .

* * 3 Л Л ^ « 0

-0.2 -0.1 0 0.1 0.2 0.3

Р

и

-0.2 -0.1 0 0.1 0.2 0.3

-0.2 -0.1 0 0.1 0.2 0.3

Рис. 6. Задача 5. Зависимость давления Р и скорости и от безразмерной координаты £ = х/с0Ь на волне разрежения. Фрагмент решения в окрестности слабого разрыва: — точное решение, ххх — результаты данной работы, □ — схема Годунова, ♦ — схема Неймана — Рихтмайера, ооо - схема Лакса — Вендрофа.

1

%

%

%

Ь=0

с Ах = 0.2 ^ = 70).

На рис. 6 приведены фрагменты решений в окрестности слабого разрыва типа 3, отделяющего область постоянного течения возле поршня от центрированной волны разрежения, полученные по четырем известным разностным схемам, и результаты аналитического решения. Видно, что предлагаемая разностная схема дает меньшие осцилляции по сравнению со схемами Неймана — Рихтмайера и Лакса — Вендрофа, а по сравнению со схемой Годунова [2] значительно меньше "размазывает" слабый разрыв.

3. Расчет откола

Немонотонность и дистракция разностного метода существенно проявляются при описании таких тонких эффектов, как места зарождения трещин или фазовых переходов. При расчете взаимодействия ударных волн и волн разрежения друг с другом или с контактными границами происходит необратимое накопление погрешностей, порождаемых осцилляционными свойствами и дистракцией разностного метода. Для оценки влияния

погрешностей, вносимых свойствами разностных схем, на расчет массы отколовшегося вещества в [12] построено аналитическое решение задачи о выходе ударной волны на свободную поверхность конденсированного вещества с последующим взаимодействием двух волн разрежения и образованием откола. Аналитическое решение удовлетворяет следующим требованиям:

1, Энтропия после выхода ударной волны на свободную поверхность постоянна во всей системе,

2, Масса откола составляет примерно треть массы всей системы для исключения влияния граничных эффектов при численном счете,

3, Откол происходит только в одной точке.

Задача о выходе нестационарной ударной волны на свободную поверхность вещества численно решалась в постановке, соответствующей аналитическому решению [12], В области с координатами [0, 2] находится конденсированное вещество, отбываемое двучленным уравнением состояния P = (y — 1)рЕ + C0,(р — p0k) с параметрами y = 3, p0k = 2.7, C0k = 3. При t = 0 заданы следующие параметры: р0 = 2.7, P0 = 0, U0 = 0, E0 = 0, Единицы измерения величин таковы: dim P — ГПа, dim E — кДж/г, dim р — г/см3, dim C — км/с. На левой границе системы задано переменное граничное условие U (t) в виде

' 1.5 при 0 < t < 0.11267067,

(/5 11267067\

1~\~т—) при °-11267067 <1 - °-4i7o7oi3>

1.368417417 при t > 0.41707013.

P=0

сетка с числом интервалов N = 70, Считалось, что откольное разрушение происходит в момент, когда давление достигнет предела прочности Pa = —1, В аналитическом решении получена зависимость минимального давления от массовой лагранжевой координаты

minP = —8.1

3.427051 (——— - 0.6632743^1 + 1 \ 64.753667 )

Рис. 7. Зависимость РШш(т): кривые 1 — аналитическое решение, 2 — расчет по схеме Годунова, 3 — расчет по схеме Лакса — Вендрофа, 4 — расчет по схеме Неймана — Рихтмайера, 5 — расчет по недивергентной схеме Куропатенко [4], 6 — результат данной работы.

Трещина, которая образуется в момент £ = 0.65033 (0,065033 мкс), откалывает от первоначальной плиты с массой 5,4 г/см2 массу Дтотк=1,84565 г/см2, Лагранжева координата этой точки т = 3.55435 г/см2.

Задача рассчитана по описанной разностной схеме и разностным схемам Неймана — Рихтмайера, Лакса — Вендрофа, Годунова и недивергентной схеме Куропатенко, На рис, 7 приведены зависимости Ршш от массовой координаты, соответствующие расчетам по этим методам и по изложенной выше разностной схеме. Зависимость РШш(т), полученная при расчете по изложенной в работе разностной схеме, в которой минимизированы осцил-ляционные и дистракционные свойства, удовлетворительно согласуется с аналитическим решением.

Заключение

Для моделирования процесса откола использованы широко известные методы [1-6]. Результаты расчетов стационарной ударной волны и волны разрежения демонстрируют (см, рис, 5 и 6), что методы Неймана — Рихтмайера и Лакса — Вендрофа дают осциллирующее численное решение как за фронтом ударной волны, так и за слабым разрывом, на котором

дР_

dx

ж+0

дР_

dx

> 0.

x-0

кк|)п. in как за ударной волной, заметно отличается от ну-

После выхода ударной волны на свободную границу от нее отражается волна разрежения, Но при этом все осцилляции, следующие за ударной волной, также последовательно отражаются от свободной границы. Суперпозиция осцилляций в области Р < 0 может привести к появлению разрушающего напряжения там, где оно отсутствует в точном решении (ложный откол). Метод Годунова дает монотонные п':

дР

так и за уединенным слабым разрывом. Однако, если ——

дж х-0

ля, как это имеет место в задаче об отколе, значение разрушающего напряжения может быть не достигнуто вообще. Иными словами, для моделирования отколов сильная дис-тракция слабых разрывов так же вредна, как и осцилляции. Это подтверждают результаты расчетов, приведенные на рис, 7, Видно, что при расчете по монотонному методу Годунова минимальное давление не достигает давления, при котором происходит разрушение, а при расчетах по разностным схемам Неймана — Рихтмайера, Лакса — Вендрофа и недивергентной схеме Куропатенко разрушение наступает из-за колебаний раньше, чем в аналитическом решении.

Список литературы

[1] Neumann J., Richtmayer R. A method for the numerical calculation of hydro dynamical shocks // J. Appl. Phys. 1950. Vol. 21, N 3. P. 232-237.

[2] годунов С.К. Разностный метод счета разрывных решений уравнений газодинамики // Мат. сб. 1959. № 47(89), вып. 3. С. 271-306.

[3] Lax P.D., Wendroff В. System of Conservation Laws // Сотр. Pure Appl. Math. 1960. N 13. P. 217.

[4] Куропатенко В.Ф. Метод расчета ударных волн // Докл. АН СССР. 1960. Вып. 3, № 4. С. 771.

[5] Куропатенко В.Ф. Математическое моделирование неустановившихся движений сред с равновесными фазовыми переходами // Вопр. атомной науки и техники. Сер. Методы и программы расчета задач мат. физики. М., 1979. Вып. 4(6). С. 3-12.

[6] Куропатенко В.Ф. Явный безусловно устойчивый разностный метод расчета течений жидкости // Числ. методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ, ИТПМ. 1984. Т. 15, № 4. С. 84-92.

[7] Куропатенко В.Ф., Макеева И.Р. Исследование дистракции разрывов в методах расчета ударных волн // Мат. моделирование. 2006. Т. 18, № 3. С. 120-128.

[8] Ландау Л.Д., Лифшиц Е.М. Гидродинамика. М.: Наука, 1988.

[9] Куропатенко В.Ф., Еськова Т.Е., Коваленко Г.В. и др. Комплекс программ "Волна" и неоднородный разностный метод для расчета неустановившихся движений сжимаемых сплошных сред // ВАНТ. Сер. Мат. моделирование физ. процессов. 1989. Вып. 2. С. 9-26.

[10] Noh W.F. Errors for calculations of strong shocks using an artificial viscosity and an artificial heat flux // J. of Сотр. Phys. 1986. Vol. 72. P. 78-120.

[11] Коробейников В.П., Мельникова H.С., Рязанов Е.В. Теория точечного взрыва. М.: Физматгиз, 1961.

[12] Куропатенко В.Ф., Макеева И.Р. О точности расчета откольного разрушения // Хим. физика. 2002. Т. 21. № 9. С. 72-78.

Поступила в редакцию 4 апреля 2006 г.