Конечно-разностная процедура случайной выборки численного решения уравнений механики сплошной среды Текст научной статьи по специальности «Математика»

УДК 519+531

КОНЕЧНО-РАЗНОСТНАЯ ПРОЦЕДУРА СЛУЧАЙНОЙ ВЫБОРКИ ЧИСЛЕННОГО РЕШЕНИЯ УРАВНЕНИЙ МЕХАНИКИ СПЛОШНОЙ СРЕДЫ

Н.А. Скоркин, В.Ф. Степанов, М.А. Шистирикова

Предлагается метод TWS численного решения задач механики сплошной среды, основанный на использовании эйлеровских конечно-разностных сеток. Слежение за движением контактных разрывов по эйлеровской сетке осуществляется с помощью процедуры случайной выборки. Счёт шага по времени разбивается на два этапа. На первом этапе по схеме Неймана определяются величины в центрах ячеек. На втором этапе расчёт переноса этих величин через границу ячеек осуществляется по схеме случайной выборки [1-3]. Приводятся результаты тестирования метода.

На примере линейного уравнения переноса продемонстрируем идею случайной выборки при дискретном представлении дифференциального уравнения переноса. Линейное уравнение переноса для функции и(х, 0 имеет вид

ди ди п п ...

— + а— = 0, а = сопв^а>0. (1)

д1дх

Простая явная схема для этого уравнения записывается следующим образом

«Г'^цм+О-йчГ. (2)

Ад:

Здесь /- номер узла на эйлеровской сетке, п - номер временного слоя, Ах = х1 , А/ - шаг счёта по времени /. Если посмотреть на уравнение (2), то можно увидеть, что значение функции и(*,/) на (и + 1) временном слое м,”+1 есть средневзвешенное значение известных величин и" и

и”_,. Это приводит к мысли: для вычисления м"+1 использовать методы теории вероятностей. Из

теории вероятностей известно, что для случайной дискретной величины £, принимающей с вероятностью Pj при проведении N испытаний значения в _/ -м испытании, математическое _ _ м

ожидание £ равно £ = ^ . Если с функцией и = м(х, ,?п+1) связать случайную величину

М

н?+1, которая может принимать только два значения: и"с вероятностью \-% и и"_х с вероятно-

,1

стью х, то математическое ожидание случайной величины и будет полностью совпадать с выражением (2). Отсюда для линейного одномерного уравнения переноса (1) получается следующий вероятностный алгоритм:

с вероятностью р = ^

и" с вероятностью р=1~Х> где р - вероятность реализации события. Алгоритм будет полностью определён, если указать способ моделирования случайной величины. Пусть £ - случайное действительное число, равномерно распределённое в интервале (0,1). Тогда вычислительный алгоритм определим так:

м„+1 = |им, если £<х, (4)

если %>х-

Применение вероятностной схемы типа (4) к системе нелинейных уравнений механики сплошной среды, её развитие и доработка осуществлены авторами работ [2, 3]. Перейдём к изложению сказанного. Расчёт временного шага At разбивается на два этапа: на лагранжевом этапе

и"+1 =

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

Р 1/ > £ 1/ , р. 1/ и др. Здесь и - скорость, р - плотность, р - давление, е - удельная внут-

/2 /2 /2

ренняя энергия. Индексом I обозначены величины в узлах эйлеровской сетки, а индексом / +1/2 - в середине эйлеровых ячеек. После определения величин первого этапа счёта шага (величины с крышечкой д) разыгрывается случайное число равномерно распределенное на промежутке (0, 1). Вероятностный алгоритм для второго этапа (численное решение уравнений переноса) для величин, определяемых в центре ячейки, имеет вид:

а) для им/2 > 0 б) для ы,+1/2 <0

уП+1 /+1/2

1/2, если £”<ДГ-^-.

Ах

Л+1/2, если

Ах

уП+1

/+1/2

I +3/2, если (1-£")>(1+А?-^-),

Ах

Л+1/2, если (1-£")>(1+А/-^-).

Ал:

(5)

Для расчёта переноса скоростей - величин, определяемых в узлах разностной сетки, апробированным [2] является следующий алгоритм:

а) для й. > 0 б) для и. < 0

ип+1 =

если тах(й.

1 1 1 Ах

Щ, если тах(ы._ ,и.)—<4п, 1 1 ' Ах

А/

м/+ь если (1+тт(м.,м., ,))—<(1-#")> ' Ал

Щ, если (1+тт(й.,и. ))—>(1-#”).

1 1+1 Ал

(6)

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

Рис. 1. Профиль скорости частиц среды ударной волной

На рис. 1 представлено решение задачи о движении сильной ударной волны по идеальному газу, имеющему уравнение состояния

р = (у-\)ре.

•0.3 -0.2 -0.1 0.0 0.1

Рис. 2. Профиль давления

Начальная плотность среды (плотность перед фронтом сильного разрыва) = 2,67, скорость движения ударной волны D = 1. На рисунке представлен профиль скорости частиц U среды за сильным разрывом. Решение, приведённое на рис. 1 в виде графика-ступеньки, есть точное решение (Exact solution), а решение, полученное по вероятностной схеме (Probabilistic scheme), представлено на рисунке кривой, помеченной маркерами-квадратиками. Усматривается вполне удовлетворительное согласование решений.

На рис. 2 на некоторый момент времени приведены результаты расчёта задачи о распаде произвольного разрыва по апробированной схеме Неймана (Newman’s scheme) и по вероятностной (Probabilistic scheme). В процессе распада произвольного разрыва в идеальном газе появляются ударная волна, волна разрежения и контактный разрыв, что отражено на рисунке. рВ начальный момент времени состояние газа определяется следующими параметрами:

при х<0 р = 1; p = \;U = 0; у = 1,4; при дс > 0 р = 0,125;/? = 0,1; U = 0, у = 1,4.

Здесь х - координата, р - плотность, р

- давление, U - скорость, у - показатель изоэнтропы. Задача решалась при следующих параметрах разностной сетки: коэффициент линейной искусственной вязкости Q = 1; коэффициент квадратичной вязкости Q2 = 0,5; число Куранта v = 0,8. Из рисунка 2 следует, что вероятностная схема даёт решение, согласующееся с решением, полученным по апробированной схеме Неймана.

Другое приложение предлагаемой здесь вычислительной схемы связано с решением стационарных задачи упругопластичности.

На рис. 3 в графическом виде приведены результаты решения задачи о стационарном распределении напряжений в толстом сферическом сосуде. Толстая сферическая оболочка нагружена внутренним давлением [4]. Внутренний радиус оболочки^ = 10 м, внешний R2 - 20 м.

Давление внутри полости Р} - 0,1 ГПа, на внешней поверхности Р2 = 0. При пределе текучести У0 = 0,0865 ГПа материал переходит в пластическое состояние в области между внутренней

поверхностью и сферой радиуса Rq =14,49 м. Вне этой области остается упругим. Согласно [4] распределение напряжений в пластической области

R{ <, R < Rq задается уравнениями

р” - 2Y0 la(R/Л,)-ij,

„вв _ „ГГ , V

Р =Р +Y0>

P = -prr-2Y0/3.

В упругой области Rq<R<R2

^=p(l-(i?2/i?)3),

^=#(1 + 0.5(Л2/Л)3),

P = -qR*/(Rl-Rо3), q = 2Y0 It^Rq/R^)- P\,

P = -P.

■8 Р . Р, GPa

0.06 -]

0.04 -

0.02 -

0.00 -

-0.02 -

-0.04 -

-0.06 -

R,m

Рис. 3. Распределение напряжений по радиусу сферического сосуда

Уравнение состояния бралось в виде P = K(p/pQ-1). Начальная плотность р0 = 2 г/см3, модуль Юнга К = 33,333 ГПа, коэффициент Пуассона v = 1/3.

Сферически симметричная (одномерная) задача о распределении напряжений в толстой оболочке решалась численно в цилиндрических координатах, т.е. в двухмерной постановке. Для подавления инерционных эффектов через промежутки времени 1 мс скорости полагались равными нулю. Распределения давления р и азимутального напряжения рвв в оболочке, полученные в расчете по вероятностной схеме (прерывистая кривая TWS) на момент 60 мс (момент установления напряжений), сравниваются с точным (сплошная кривая Exact solution) решением, см. рис. 3. Из рисунка следует, что имеет место быть вполне удовлетворительное согласование сравниваемых параметров. Надо отметить, что результат удивителен. Дело в том, что искалось стационарное решение по вероятностному, существенно нестационарному по определению, алгоритму TWS.

Однако описанный выше численный метод TWS не является достаточно надежным при расчете положения границ различных веществ. Дело в том, что для расчета перемещения термодинамических параметров задачи, проводимого по схеме случайной выборки (5), и для расчета перемещения кинематических параметров задачи по схеме (6), используются разные скорости. Это обусловлено тем, что разностной ячейке приписана одна плотность рм/2 > °ДН0 давление Рм/2 ’ °Дна энергия sj+U2 , в то время как ячейка (одномерная) имеет две скорости и, и им. И это иногда в силу накопления вычислительных ошибок (например, в результате округления чисел в компьютерном коде) приводит к рассогласованию движения поля скоростей и поля термодинамических параметров по сетке.

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

Приступим к изложению полученного в результате метода SPH&MK, на примере одномерной газовой динамики. Будем предполагать, что в центре разностной сетки определяются все параметры задачи. Расчёт временного шага At разбивается на два этапа: на лагранжевом этапе численно по ХРЯ-методу решаются уравнения газовой динамики без конвективных слагаемых. На втором этапе решаются уравнения переноса по схеме (5).

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

Здесь где г= дсу-дсу - расстояние от центра текущей /-й ячейки (частицы) до соседней /-й

ячейки, и = г/И, к - область определения сглаживающей функции IV. Разностные уравнения, например, для уравнения импульсов на первом этапе счета шага А/ в форме БРН-аппроксимации для задач гидродинамики в одномерном случае имеют вид [1-4]

W{r,h) = —у 7(2-i>)3, лЬґ 4

1--и2+-іД 0<< 1;

2 4

-(2-of, 1<и<2;

0.

2<о.

С7)

В этих уравнениях т] есть масса у'-й ячейки. После определения по схеме (7) параметров

первого этапа счёта шага разыгрывается случайное число є (0,1). И алгоритм случайной выборки для второго этапа счета шага по времени для всех величин теперь имеет единый вид (5). Для иллюстрации возможностей предлагаемого метода рассмотрим задачу о распадении произвольного разрыва в идеальном газе. Начальные значения параметров газа слева и справа от поверхности контактного разрыва взяты из работы [6], там же приведены формулы для расчета точного решения на моменты времени />0. Расчетная область и число ячеек в ней в полном соответствии с цитированной работой [6]. Результаты расчетов для удельной внутренней энергии на момент времени 1-4 мс представлены на рисунке 4. При расчете распределения удельной внутренней энергии и давления по стандартному варианту БРН-метода на поверхности контактного разрыва возникает скачок энергии и давления, см. рис.4. При расчете по вероятностной схеме БРН&МК скачок отсутствует.

С /Е|

1.1

1.0

оя

0.8

0.7

0.6

0.5

0.4

В расчетах по вероятностной схеме SPH&MK по причине наличия разностной сетки поиск соседей осуществляется один раз в начальный момент времени. Перемещение частиц по сетке осуществляется по схеме со случайной выборкой (5). Вследствие этого расчет задачи по вероятностной схеме приблизительно в 10 раз быстрее по времени, чем расчет по SPH-методу. В расчетах использовалась только квадратичная вязкость, предложенная в работе [6],

9 = -80р,(ДсНум,)2. Численные методы, использующие SPH-идеологию, особенно эффективны в таких задачах, как расчет проникания тел в преграды, задачи о заглубленных взрывах, задачи астрофизики и т.п.

При разработке осушителей природного газа, работающих на принципе дросселирования давления и представляющих собой трубу переменного сечения, интерес представляет распределение параметров течения газа по длине трубы. Первоначально покоящийся газ в трубе переменного сечения, каковой является, например, сопло Лаваля, разгоняется до скорости звука в горловом сечении, потом расширяется в сверхзвуковой части сопла при заданном давлении на срезе сопла. Из всех уравнений газовой динамики только в уравнение неразрывности входит зависимость S(x), описывающая площадь сечения трубы как функцию координаты jc точки, расположенной на оси симметрии трубы

0.3

1 і I |_ I— 1

I I

г гз l_ _l>

' г rrri L

I I I I kill |\| I I |_ I і I Ц Гї\ I

0.0

0.2

0.4

0.6

0.8

1.0

XIX,

Рис. 4. Распределение внутренней энергии на конечный момент времени:

1 - точное решение; 2 - решение по методу ЭРН&МК; 3 - ЭРН, стандартный вариант

др дри S'(x)

— + -*-— + ри-dt дх S(x)

■ О.

Остальные уравнения имеют стандартный вид. Для расчетов были опробованы 4 схемы: выше изложенная TWS, Моретга, МакКормака и Лакса-Вендроффа [7]. Схемы МакКормака и Лак-са-Вендроффа дали схожие результаты, поэтому на рис. 5. они идентифицированы одним маркером «McCormack». Длина трубы бралась равной 3 м. Образующая сопла описывалась кривой у = х/2 + 1/2, 0 < х £ 3 м. В качестве газа рассматривался идеальный газ с уравнением состояния р = (у-1)ре,у = 1,4.

В начальный момент времени газ в трубе покоится, давление р-= 20 атм. На правом конце трубы поддерживалось давление р = 0,15р*. Граничные условия (на входе сопла и на срезе со-

пла) для скорости 1/(х) задавались в виде - =0. На рис. 5 приведены результаты расчетов

дх

распределения скорости V по длине трубы х. Результаты расчетов показали, что схемы МакКормака и Лакса-Вендроффа обладают сильным сглаживающим эффектом и при расчете задач на установление они непригодны.

Численный метод Моретти, использующий информацию о характеристических направлениях им,мс кой системы уравнений газовой

X, м

Рис. 5. Распределение скорости по длине сопла

динамики, наиболее точен. Схема дала результаты расчетов, близкие к тем, что получены по схеме Моретти. Таким образом, при проведении расчетов движения идеального газа в квазиодномерной постановке предпочтение следует отдать схеме Моретти.

Однако применение этой схемы в многомерном случае или, при использовании сложного по конструкции уравнения состояния, применение схемы TWS проще. Связано это с тем, что схема Моретти использует информацию о характеристических направлениях. А для этого необходимо уравнения движения записывать в специальном виде, т.е. уравнение состояния в вычислительном коде затруднительно представить в виде отдельного программного модуля.

Далее на рис. 6 приводится сравнение результатов расчетов, полученных по двухмерной программе TWS и по квазиодномерной схеме Моретти. Рассматривалось движение в трубе диамет-Тетрегйиге ром 57 мм, длиной 320 мм, на

оси которой на расстоянии 100 мм от левого конца располагался дроссель.

Дроссель - представлял тело, полученное вращением сегмента окружности радиуса 50 мм.

Высота сегмента 25 мм. Минимальный зазор между дросселем и поверхностью трубы получился равным 3 мм.

На рис. 6 представлены распределения температуры (шкала Цельсия). На этом рисунке маркер «TWS2D» соответствует двухмерным расчетам по разностной схеме 'РМЗ. Маркеру «МогеШ» соответствует расчетная кривая, полученная по квази одномерной программе, в которой реализована ко-нечно-разностная схема Моретти. Граничное условие на левом конце трубы такое же, как в выше

приведенном случае сопла Лаваля. На правом конце трубы поддерживалось давление р = 0,5р . Видно вполне удовлетворительное согласование сравниваемых параметров. Из этого следует, что двухмерная вычислительная программа пригодна для решения рассмотренных задач.

Рис. 6. Распределение температуры по длине трубы

Расчеты по вычислительной методике SPH&MK не проводились в силу неспецифичности задачи для данной численной процедуры.

Работа выполнена при финансовой поддержке фонда содействия развитию малых форм предприятий в научно-технической сфере, проект М 443бр/6842 от 29.06.2006 г.

Литература

1. Сод Г. Обзор некоторых разностных методов для консервативных нелинейных гиперболических систем// J. of Comp. Phys. - 1978. - V.27. - № 1. - P. 1-14.

2. Модификация метода Глимма к задачам проникания / С.Г. Андреев, В.В. Башуров, В.А. Свидинский, Н.А. Скоркин // ВАНТ, серия «Методики и программы численного решения задач математической физики». - 1985. - Вып. 3. - С. 80-88.

3. Скоркин Н.А. Об одной конечно-разностной схеме для задач механики сплошной среды// Доклады международной конференции «IV Забабахинские научные чтения». - Снежинок: РФЯЦ - ВНИИТФ, 1995. - С. 277-279.

4. Качанов JI.M. Основы теории пластичности. - М.: Наука, 1969. - С.103-105.

5. Libersky L. D., Petschek A.G. etc. High Strain Lagrangian Hydrodynamics (3D SPH code for Dynamic Material response) // J. of Comp. Phys. - 1993. - № 109. - P. 67-75.

6. Паршиков A.H. Применение решения задачи Римана в методе частиц// Журнал вычислительной математики и математической физики. - 1999. - Т. 396. - № 7. - С. 1216-1225.

7. Флетчер К.А. Вычислительные методы в динамике жидкостей. - М.: Мир, 1991, Т. 2. -С.172-175, С. 186-190.

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