ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2012, том 52, № 4, с. 672-695
УДК 519.633
МОНОТОННЫЕ КОМПАКТНЫЕ СХЕМЫ БЕГУЩЕГО СЧЕТА ДЛЯ СИСТЕМ УРАВНЕНИЙ ГИПЕРБОЛИЧЕСКОГО ТИПА1)
© 2012 г. М. Н. Михайловская*, Б. В. Рогов**
(*141700Долгопрудный Московской обл., Институтский пер., 9, МФТИ (гос. ун-т);
**125047Москва, Миусская пл., 4, ИПМим. М.В. Келдыша РАН) e-mail: rogov@post.ru, m.n.mikhailovskaya@gmail.com Поступила в редакцию 22.06.2011 г.
Для квазилинейных уравнений гиперболического типа представлены консервативные абсолютно устойчивые компактные схемы, монотонные в широком диапазоне значений локального числа Куранта. Они имеют четвертый порядок аппроксимации по пространственной координате на компактном шаблоне и нечетный (первый или третий) порядок аппроксимации по времени. Схемы экономичны и решаются методом бегущего счета. Приводится детальное исследование скорости сходимости предложенных схем при сгущении разностной сетки для различных порядков гладкости решения. Возможности схем продемонстрированы на примере решений известных одномерных тестовых задач для уравнений газовой динамики. Библ. 31. Фиг. 17. Табл. 8.
Ключевые слова: квазилинейные уравнения гиперболического типа, компактные разностные схемы, монотонность, бегущий счет.
1. ВВЕДЕНИЕ
Уравнения и системы уравнений гиперболического типа составляют значительную часть математических моделей, используемых для решения разнообразных прикладных задач (см. [1]). Среди них широкий класс задач связан с нахождением обобщенных решений, которые описываются кусочно-гладкими функциями и могут иметь скачкообразный характер. Важнейшим свойством, которым должны обладать разностные схемы сквозного счета для численного решения данного класса задач, является их монотонность (см. [1]). Для вычислительной практики представляют интерес не только монотонные схемы, но и более широкий класс диссипативных схем (см. [2], [3]), обладающих аппроксимационной вязкостью. Последний класс схем допускает немонотонность численного решения вблизи скачков искомых функций, но эта немонотонность, как правило, пространственно локализована и невелика по амплитуде. Другое важнейшее требование, предъявляемое к схемам сквозного счета, — наличие консервативности (см. [4]).
В настоящее время среди схем высокого (выше первого) порядка аппроксимации большую популярность получили компактные схемы (см. [5]—[8]). При построении таких схем используются дифференциальные следствия исходных уравнений, что позволяет при разностной аппроксимации производных использовать небольшое число точек шаблона (компактный шаблон). Порядок компактной аппроксимации больше или равен числу точек шаблона (см. [5], [6]). Обычно при построении компактных схем производные от искомых функций также рассматриваются в качестве искомых переменных (см. [9]). Этот прием позволяет легко сформулировать граничные условия с высокой точностью на компактном разностном шаблоне. Компактность шаблона обеспечивает экономичность неявных компактных схем: разностные уравнения решаются либо прогонкой (см. [5], [6], [9]), либо бегущим счетом (см. [10]).
Компактные схемы обладают свойством консервативности (см. [5], [6], [11]). Схемы с нецен-трированными пространственными аппроксимациями (см. [5], [6]), ориентированными против потока, являются устойчивыми и диссипативными при подходящей аппроксимации производных по времени (см. [5], [6]). Однако для точного воспроизведения структуры скачков искомых функций с помощью этих схем нечетного порядка аппроксимации по пространственным пере-
Работа выполнена при финансовой поддержке РФФИ (коды проектов 09-01-00728, 11-01-00504) и Роснауки (госконтракт 02.740.11.0615 и П594).
менным приходится прибегать к монотонизаторам в виде ограничителей потоков (см. [12], [13]). Другим подходом для построения неосциллирующих компактных схем является применение процедур ENO и WENO для расчета потоков через границы разностных ячеек (см. [7], [8]). Однако этот подход полностью не устраняет немонотонность численного решения вблизи скачков (см. [14]). Кроме того, процедуры ENO и WENO достаточно трудоемки. Альтернативным подходом является построение монотонных бикомпактных (двухточечных компактных) схем с симметричными компактными аппроксимациями четного порядка по пространственным переменным и аппроксимациями по времени нечетного порядка (см. [10], [15]).
В [10] для линейного уравнения переноса были предложены бикомпактные разностные схемы, имеющие на двухточечном шаблоне четвертый порядок аппроксимации по пространственной координате. Из них схема первого порядка аппроксимации по времени является монотонной. На ее основе в [10] построена гибридная нелинейная монотонная схема, имеющая на гладких решениях третий порядок точности по времени. Обе схемы являются абсолютно устойчивыми и экономичными. Расчет по ним ведется по явным формулам бегущего счета.
В настоящей работе схемы из [10] обобщены на случай системы квазилинейных уравнений гиперболического типа, описывающих законы сохранения. Обобщенные схемы являются консервативными, абсолютно устойчивыми, а также монотонными в широком диапазоне значений локального числа Куранта. Расчет по схемам осуществляется по явным формулам бегущего счета, что говорит о простоте реализации предложенных схем. В построенных схемах не использованы искусственная вязкость и какие-либо ограничители потоков. Исследован вопрос о скорости сходимости предложенных схем при сгущении разностной сетки для различных порядков гладкости решения. Приведены результаты расчетов, демонстрирующие возможности предложенных схем на примере решения известных одномерных тестовых задач газовой динамики (см. [16], [17]).
2. КОМПАКТНАЯ СХЕМА ПЕРВОГО ПОРЯДКА АППРОКСИМАЦИИ ПО ВРЕМЕНИ
Рассмотрим сначала смешанную задачу Коши (см. [2]) для одного квазилинейного гиперболического уравнения, записанного в дивергентном виде:
ut + f (u)x = 0, a (u) = df(u)/du > 0,
u(x, 0) = u0(x), x > 0; u(0, t) = ^(t), t > 0.
Разностную схему для начально-краевой задачи построим с помощью метода прямых и интегро-интерполяционного метода (см. [18]). Для этого введем неравномерную сетку {xy, j> 0} на интервале [0, да). На временном слое t = const проинтегрируем уравнение (1) на отрезке xj< x < +1, используя формулу Симпсона (см. [2]):
J F(x)dx = j(Fj+i + 4Fj+i/2 + Fj) + 0(h5+i), hj+i = xy+i - Xj, (2)
xj
для F(x) = u. В результате получим дифференциальное по t и разностное по x уравнение
6 ((i + 4uj+i/2 + uj)t =-h (/+i - fj) fj = f /), j > a (3)
выполняющееся с точностью до O(h4) на точном решении уравнения (1). Здесь и в дальнейшем индекс при пространственном шаге h опущен.
Вычислим интегральное среднее от функцииf(u) на отрезке xj < x < xj +1:
xj+i
w^ , =1
j+i-h J f (u)dx (4)
к
х;
двумя различными способами с точностью до O(h4): по формуле Симпсона (2) и по формуле Эйлера—Маклорена (см. [2]):
х;+1 2
| F(x)йх = к((+1 + р)- к-2; - р';) + 0(к5) (5)
х;
для F(x) = f(u)/h. В результате получим соответственно следующие формулы:
Щ„ = f + VJttn+fj)+O(h'), f¡'+1/2 = /((+,/2), (6)
=1 ( + ,2 C/j- j Oh<). (7)
Исключим из формулы (7) значения производнойf в целых узлах с помощью уравнения (1), а затем приравняем друг другу выражения (6) и (7). В результате получим дифференциальное по t и разностное по x уравнение
h - uj)t =-((+1 - 2/+1/2 + /), ^ ° (8)
выполняющееся с точностью до O(h4) на точном решении уравнения (1).
Приведенные выше выкладки для скалярного уравнения (1) можно полностью повторить и в случае системы уравнений
u, + f(u) x = 0, (9)
где u(x, t) — искомая вектор-функция c m компонентами, f(u) — заданная вектор-функция размерности m. В результате можно получить системы обыкновенных дифференциальных уравнений (ОДУ), аналогичные (3), (8).
Подчеркнем, что именно добавление в дифференциально-разностную схему дополнительного уравнения (8) отличает данную схему от известной трехточечной симметричной (центрированной) компактной схемы четвертого порядка аппроксимации, в которую входит только уравнение (3) (см. [5]). Нетрудно увидеть, что уравнение (8) представляет собой разностную аппроксимацию дифференциального следствия уравнения (1):
(и + fx) x =
Важно отметить, что дифференциально-разностная схема (3), (8) может быть решена методом бегущего счета. Пусть функция uk(t) известна (например, при k = 0 она известна из граничного условия (1)), тогда из системы двух ОДУ (3), (8) при j = k можно определить две функции uk + 1/2(t), uk + 1(t) при заданных начальных данных, которые определяются из начального условия (1). Это означает, что uk + 1/2(t) являются вспомогательными функциями при расчете функций uk + 1(t) в целых узлах. Фактически схема (3), (8) является двухточечной компактной (т.е. бикомпактной в терминологии [19], [20]) схемой. На неравномерной сетке, состоящей из целых узлов, схема сохраняет четвертый порядок точности.
Введем, как в [20], первообразную функцию v от f(u):
vx = / (u). (10)
Интегрируя это уравнение на отрезке Xj < x < Xj +1 и используя формулу (6), получим связь между значением uJ +1/2 в полуцелом узле и значениями функций u и v в целых узлах:
h(+1 - Vj) = W+1 = 6/ И+1) + 4/(lij+i/2) + f (Uj)] + °(й4). (И)
Далее, исключая из схемы (3), (8) функцию uj + 1/2(t) с помощью формулы (11), можно получить бикомпактную схему для определения функций uy(t), vj(t) в целых узлах сетки. Однако, в отличие от схемы (3), (8), эту схему намного сложнее решать, за исключением случая, когда функция f(u) является линейной:
f (и) = au, a = const > 0. (12)
В этом случае схема (3), (8) после исключения из нее функции uJ + 1/2(t) с помощью формулы (11) переходит в бикомпактную дифференциально-разностную схему из [10], [20] для определения функций uy(t), Vj(t) в целых узлах. Рассматривая функции uj + 1/2(t) в бикомпактной схеме (3), (8) как вспомогательные функции, которые подобны функциям vj(t) в бикомпактных схемах из [10], [20], перепишем систему ОДУ (3), (8) в следующем виде:
6(Uj+1 + 4Wj+1 + Uj) =-h[f (Uj+1)- f (Uj)] j > 0 (13)
4 Ы - uj) =-[ (j)-2/ (j) + / () j * 0, (14)
где Wj(t) = Uj_ 1/2
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.