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

---------------------□ □-------------------------

Метод сумісної апроксимації для побудування різничних схем підвищеного порядку точності розширено для випадку квазілінійних гіперболічних рівнянь та їх систем. Наведено результати побудування різничних схем аж до одинадцятого порядку точності за часом та простором для одновимірного рівняння Бюргерса та системи рівнянь одновимірної динаміки нев’язкого сти-скаємого газу. Результати тестових розрахунків підтверджують теоретичні висновки

Ключові слова: метод сумісної апроксимації, кінцево-різничні схеми, високий порядок точності

□------------------------------------□

Метод совместной аппроксимации для построения разностных схем повышенного порядка точности расширен на случай квазилинейных гиперболических уравнений и их систем. Приведены результаты построения разностных схем вплоть до одиннадцатого порядка точности по времени и пространству для одномерного уравнения Бюргерса и системы уравнений одномерной динамики невязкого сжимаемого газа. Результаты тестовых расчетов подтверждают теоретические выводы

Ключевые слова: метод совместной аппроксимации, конечно-разностные схемы, высокий порядок точности ---------------------□ □-------------------------

1. Введение

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

2. Анализ исследований и публикаций

Ранее автором был предложен новый метод построения компактных разностных схем повышенного порядка точности - метод совместной аппроксимации (МСА) [1, 2, 3]. В основе этого метода лежит идея повышения порядка аппроксимации не за счет раздельного повышения порядка малости составляющих невязки, возникающих при аппроксимации производных по раз-

з.............................................

УДК 519.6

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

В. Л. Бучарски й

Кандидат технических наук, доцент Кафедра двигателестроения Днепропетровский национальный университет

им. О. Гончара

пр. Гагарина, 72, г. Днепропетровск, Украина,

49000

Е-mail: bucharsky@mail.ru

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

При построении разностных схем в работах [1, 2, 3] использовался метод неопределенных коэффициентов. И если в случае линейных уравнений проблем не наблюдалось, то при построении разностной схемы для квазилинейного уравнения Бюргерса (см. [2]) возникли сложности, связанные с переменностью неопределенных коэффициентов и экспоненциальным ростом их числа в зависимости от порядка аппроксимации. Это не позволило предложить надежный подход для построения разностных схем МСА в нелинейном случае.

3. Цель настоящей работы

Цель настоящей работы состоит в разработке подхода для применения МСА к построению разностных

© В. Л

схем наперед заданного порядка аппроксимации по пространству и времени для решения квазилинейных гиперболических уравнений и их систем в одномерном случае

4. Методика построения разностных схем МСА в квазилинейном случае

Рассмотрим в качестве исходного следующее одномерное квазилинейное уравнение, записанное в консервативной форме:

.k

д

df + М) = 0.

dt дх

(1)

df f"+1 - fn

dt

--I

dk+1f

T k=i (k +1)! dtk

Эта запись следует из разложения разности ■П"1 -^¡п) в ряд Тейлора в окрестности точки (^^) . Тогда исходное уравнение принимает вид

(f,

fn+1 - fn ,dg(f) f,

- + — 1

dk+1f

T dx k=i (k +1)! dtk

= 0.

(2)

Отметим, что все производные определены в точке

(tn,xi ) .

Далее, в соответствии с идеей МСА, выразим производные по времени через производные по пространству, используя Г-форму дифференциального представления (2). В результате имеем следующие выражения для производных по времени от искомой функции:

dk+1f

dk

= —"Цт |(-A)k — I, k = 1, ~.

dtk+1 dxk ’ dx 1

(3)

Тождества (3) в общем случае несложно доказать методом математической индукции. Подставляя (3) в полудискретную форму (2), получим

Г-£ + dgf + i*

t dx ^=1 (k +1)! dxk

(-A)k dx = 0.

(4)

В (4) все производные берутся по пространственной координате на п - м временном слое. Отметим, что уравнение (4) полностью эквивалентно уравнению (1). Погрешность появится при учете конечного числа слагаемых в сумме в (4). Как видно, для получения разностной схемы порядка точности М необходимо учесть (М -1) слагаемых в сумме в уравнении (4).

Из вида (4) следует простая двухшаговая методика построения разностной схемы заданного порядка точности М для решения уравнения (1).

Шаг 1. Строится разностная аппроксимация про-

изводнои

здесь Г = Г (^х) - искомая функция независимых переменных ^х, g = g (Г) - поток величины Г. Отметим, что в случае гиперболичности уравнения (1) существует действительная функция А(Г) = — , называемая якобианом

Эг

системы.

В области {0 < х < 1,0 < t < Т} построим равномерную сетку ю4 =^ = Ш|1 = 1,^;^ = пт|п = 1,^}, где т,Ь -шаги дискретизации по времени и пространству соответственно. В дальнейшем будем полагать Гп = Г(tn,xi). Ранее было обнаружено, что при применении МСА использование неявного шаблона не приводит к безусловной устойчивости даже в линейном случае [1]. Это обусловлено введением в структуру разностной схемы исходного гиперболического уравнения, область зависимости решения которого ограничена. Следовательно, будем проводить все рассуждения для явной двухслойной разностной схемы. Рассмотрим уравнение (1) в точке (^^). Запишем его в полудискретной форме,

ЭГ

для чего представим производную — в следующей форме:

dx

порядка точности M

Ф, = l1m (gn ) =

dg(f)

dx

-0(hM).

(5)

Здесь LIM - некоторый разносный оператор к - го поряд-

Эк

ка, аппроксимирующий производную —| с порядком

Эх

точности М. В качестве его могут быть использованы любые разностные операторы - от простейших центрально-разностных до широко распространенных в последние годы монотонных WENO аппроксимаций.

Шаг 2.Значение искомой функции на новом временном слое определяется по следующей зависимости:

f,n+1 = f,n - T ■ ф,-I

LM-k ((-A)k-1 Ф,).

(6)

Задача построения операторов LM тривиальна -можно, например, воспользоваться методом неопределенных коэффициентов или неким другим подходом (см., например,[5, 6]).

По методике (5, 6) несложно построить разностные схемы для систем квазилинеИных уравнении. В

dg

этом случае якобиан A(f) = -j будет матрицей. Отме-

df

тим существенную простоту предлагаемого подхода по сравнению с предложенным в [7, 8]. В цитируемой работе для построения алгоритма необходимо вычислять и манипулировать со старшими производны-dkg

ми A(k)(f) = djT , которые в случае систем уравнений

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

5. Пример построения конечно-разностных схем

В работе приведены результаты построения

разностных схем нечетного порядка точности

dk

M = {3,5,7,9,11} . Для аппроксимации операторов —^

dxk

были использованы простейшие центрально-разностные операторы L1^ по причинам, изложенным в [3]. При этом для аппроксимации на первом шаге dk

первой производной —г использовался кососим-dxk

k

k

T

k

T

метричный оператор , для аппроксимации

Э2к11

нечетных производных д 2к11 на втором шаге - кососимметричные операторы L2c|(M_2k), а для аппрокси-

Э2к

мации четных производных э 2к на втором шаге -

симметричные операторы L2ckM_2k+1). При этом, для схемы, например, пятого порядка точности:

ф,=ц* & +°(ь6)■

I- ( а ї)=^ <л'ф')+°<ь<>■■

& ( а’ ї )=^ (а’фі)+°<ь‘ >

& ( аз і )=^ (а’Ф')+° <ьг )■• I?( а< і )=^<а>')+° <ьг).

Выражение для операторов Lc(j) несложно получить, используя классические методики построения конечно-разностных аппроксимаций (см., например, [5, 6]). Например,

^ (ц )=-^-2 + 4-1 - 4+1 + %2

ь

Для подавления высокочастотных осцилляций в правую часть (6) в настоящей работе вводилась диссипативный член порядка М

V, = (^Ь,“ (їі" )).

Величина параметра а в последнем выражении определялась хорошо известным minmod-лимитером потоков [9], который обеспечивает невозрастание полной вариации решения. Естественно, использование монотонизирующих процедур, настроенных на конкретную схему, привело бы к лучшим результатам, однако такие исследования выходят за рамки настоящей работы.

Анализ устойчивости полученных разностных схем спектральными методами и методом дифференциального приближения показал, как и ожидалось, их условную устойчивость при выполнении условия Куранта-Фридрихса-Леви

с|Хш

1< 1.

Ь

здесь Хтах - максимальная скорость распространения возмущений в решаемой задаче.

6. Результаты тестовых расчетов

Для практической оценки проведенных построений был решен ряд общепринятых тестовых задач для

линейных и квазилинейных гиперболических уравнений и их систем. В тестовых задачах на проверку порядка аппроксимации оценка погрешности решения є проводилась в конечномерных аналогах норм в іі-і,- С .

і.Иї = ЬХІ4 11и4 =*ь£1 и.1 -

и, и = шахи,

II НІС і=1,к 1 11

здесь и - некоторая сеточная функция, N - число узлов сетки по пространству. Порядок аппроксимации определялся по экстраполяции Ричардсона (см., например, [9]). Во всех задачах полагалось о = ^ = 0.5.

Квазилинейное скалярное уравнение. Была решена следующая задача для уравнения Бюргерса:

Э£ _Э_

Эt Эх

= °-

ї (°-х) = °.5 + sin (пх), х є[-1-1)

с периодическими граничными условиями. Расчет проводился до момента времени t = °.^Л - когда решение все еще гладкое. Результаты расчетов- представленные на рис. 1- показали- что все предложенные разностные схемы подтвердили заявленный порядок аппроксимации.

Рис. 1. Зависимость С-нормы погрешности численного решения |ЩС от числа точек в расчетной области N для квазилинейного скалярного уравнения

Здесь и далее к - оценка порядка аппроксимации, вычисленная по результатам расчетов. Отметим, что введение искусственной диссипации несколько снижает порядок точности. Решение этой задачи подтвердило работоспособность предложенных схем при знакопеременном якобиане.

Система квазилинейных уравнений. В качестве системы уравнений были рассмотрены уравнения Эйлера динамики невязкого сжимаемого газа в одномерной постановке:

м+м>=°

Эt Эх ■

V / \ ( ри )

ї= ри ■ g= ри2 + р

1Е) ч и (Е+р))

(7)

=1

=1

здесь р,и,Е,р - плотность, скорость, полная энергия и «Точное» решение этой задачи взято из [8]. Зна-

давление газа соответственно. Замыкалась система (10) чение плотности газа для t = 1.8 представлено на

уравнением состояния идеального газа рис. 3.

Р = (у-1)

Е -

ри

2

2

с показателем адиабаты у = 1,4 .

Для оценки скорости сходимости система (10) решалась со следующими начальными данными

р = 1+0^т (пх); и = 1; р = 1

и с периодическими граничными условиями (задача распространения энтропийных волн [7]). Известное точное решение этой задачи

р = 1+0^т (п(х - ^); и = 1; р = 1

позволило оценить скорость сходимости предложенных разностных схем. Приведенные на рис. 2 зависимости L1 -нормы погрешности определения плотности газа от шага дискретизации в логарифмических координатах подтверждают заявленный порядок аппроксимации и в случае систем квазилинейных уравнений.

Рис. 2. Зависимость L1 -нормы погрешности численного решения от числа точек в расчетной области N для системы квазилинейных уравнений

В следующей серии расчетов были оценены способности предложенных разностных схем отрабатывать разрывные решения. Решение классических задач о распаде разрыва (задачи Лакса и Сода) не показало сколько-нибудь значимых преимуществ предлагаемых схем перед хорошими TVD схемами второго порядка.

Следует ожидать большей выгоды от применения разностных схем высокого порядка точности при решении задач, в которых присутствуют сильные и слабые возмущения. Для проверки этого была решена задача о взаимодействии сильной ударной волны ( Мж = 3) с энтропийной волной [10]. В расчетной области х е[-5,5] задавались следующие начальные условия

(р-и-р) =

(3.857143- 2.629369- 1°.333333)-х <-4 (1 + °^іп (5х)-°-1)-х > -4.

Рис. 3. Зависимость плотности газа р от пространственной координаты х

В результате серии расчетов для каждой разностной схемы определялись минимальное число ячеек сетки КШ1П - при котором отклонение полученного решения от «точного» в норме І1 - не превышает °-°1 (за исключением зоны головной ударной волны) и время ТШП - необходимое для проведения расчета на такой сетке.

Эти данные для предложенных разностных схем приведены на рис. 4.

Рис. 4. Зависимости времени расчета Т™ (белые прямоугольники) и потребного числа точек расчетной области (черные прямоугольники) от порядка точности схемы

Здесь по оси абсцисс указан порядок точности схемы, темные прямоугольники отражают число (в тысячах) ячеек N1:““ , а белые - затраченное на расчет относительное время.

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

Е

7. Выводы

В заключение можно отметить следующее:

• в работе приведена и обоснована методика построения двухслойных двухшаговых разностных схем метода совместной аппроксимации высокого порядка точности для решения квазилинейных гиперболических уравнений в частных производных и их систем в одномерном случае;

• с помощью предложенной методики получены разностные схемы М = {3,5,7,9,11} порядков аппрокси-

мации по времени и пространству для решения одномерных квазилинейных гиперболических уравнений в частных производных и их систем;

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

В дальнейшем предполагается обобщить предложенный подход на случай нескольких пространственных переменных.

Литература

1. Бучарский, В.Л. Метод совместной аппроксимации построения разностных схем для решения уравнений в частных про-

изводных / В.Л. Бучарский // Техническая механика. - 2007. - № 1. - с. 50 - 57.

2. Бучарский В.Л. Разностная схема метода совместной аппроксимации для решения квазилинейных гиперболических

уравнений / В.Л. Бучарский // Проблеми обчислювальної механіки і міцності конструкцій. - Дніпропетровськ: ДНУ, 2008 . - с. 7 - 15.

3. Бучарский В.Л., Калинчук Е.М. Симметричные разностные схемы метода совместной аппроксимации для решения

линейного уравнения переноса / В.Л. Бучарский, Е.М Калинчук // Математичні машини і системи. - 2011. - №4. -с. 161-165.

4. Shokin Y.I. The Method of Differential Approximation / Shokin Y.I. - Springer-Verlag Berlin And Heidelberg Gmbh & Co.

- 1983. - 224р.

5. Самарский А.А. Теория разностных схем. Учебное пособие / Самарский А.А. - М.: Главная редакция физико-математи-

ческой литературы, 1980. - 616 с.

6. Tannehill J. C. Computational Fluid Mechanics and Heat Transfer / J.C. Tannehill; R.H. Pletcher; D.A. Anderson. - Hemisp-

here Pub, 1997. - 792p.

7. Qiu J., Shu C.-W. Finite-difference WENO schemes with Lax-Wendroff-type time discretizations / J. Qiu, C.-W. Shu // SIAM

J.Sci.Comput. - 2003. - v.24. - №6. - p.2185-2198.

8. Qiu J. Hermite weno schemes with Lax-Wendroff type time discretizations for Hamilton-Jacobi equations/ J. Qiu // Journal of

Computational Mathematics. - 2007. - v.25 - p.131-144.

9. Encyclopedia of Computational Mechanics Volume 1 Fundamentals /Editors Erwin Stein, Rene de Borst, Thomas J. R. Hughes -

WILEY, 2004 - 798 p.

10. Shu C.-W., Osher S. Efficient implementation of essentially non-oscillatory shock capturing schemes / C.-W. Shu, S.Osher // J. Comp. Phys. - 1988. - v.77. - p.439-471.

3