О слабой сходимости на разрывных решениях TVD-схемы Хартена второго порядка аппроксимации

Вычислительные технологии

Том 2, № 5, 1997

О СЛАБОЙ СХОДИМОСТИ НА РАЗРЫВНЫХ РЕЩЕНИЯХ TVD-СХЕМЫ ХАРТЕНА ВТОРОГО ПОРЯДКА АППРОКСИМАЦИИ *

В. В. ОСТАПЕНКО Институт гидродинамики им. М. А. Лаврентьева СО РАН Новосибирск, Россия

The notion of a weak convergence of the difference solution to the exact generalized solution of the approximated quasilinear hyperbolic systems of conservation laws is introdused. It is shown that the TVD Harten scheme (wich has the second order approximation on smooth solutions) has only the first order of weak convergence in calculating the problem of the breaking of a dam with formations of discontinuity and depression waves. The obtained first order of weak convergence results in the similar decrease in the order of pointwise strong convergence on smooth parts of the exact solution behind the discontinuous wave front.

1. Введение

В настоящее время широкое распространение получили разностные схемы повышенной точности для сквозного расчета обобщенных [1, 2] (слабых [3]) решений гиперболических систем законов сохранения (в частности, системы законов сохранения газовой динамики). Однако в большинстве работ, посвященных построению таких схем (см., например, [4-11]), под точностью схемы понимается порядок аппроксимации на гладких решениях, что не гарантирует соответствующего повышения точности на разрывных решениях.

С другой стороны, при построении и анализе разностных схем, предназначенных для расчета обобщенных решений, применяются методы, использующие слабую аппроксимацию дифференциальных уравнений в различных интегральных нормах. Для эллиптических уравнений теория таких методов была развита в [12]. Для гиперболических систем законов сохранения различные вопросы, связанные со слабой аппроксимацией и сходимостью разностных схем, изучались в работах [13-19]. В частности, в [17] были получены необходимые и достаточные условия слабой аппроксимации (в том числе и с повышенным порядком) для однородных разностных схем, а в [18] — достаточные условия слабой аппроксимации для схем, заданных на неоднородной сетке.

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований гранты №96-01-01546, 96-01-01547 © В. В. Остапенко, 1997.

является порядок их слабой сходимости в соответствующей интегральной норме. В настоящей работе предлагается метод, позволяющий эффективно оценивать порядок слабой сходимости разностной схемы на обобщенных решениях аппроксимируемой гиперболической системы законов сохранения. Этот метод (идея которого принадлежит С. К. Годунову и В. С. Рябенькому) основан на экспериментальной проверке скорости сходимости первых интегралов от получаемого разностного решения, которые берутся по областям, содержащим различные особенности аппроксимируемого точного решения. При помощи указанного метода проведен анализ точности ТУБ-схемы Хартена [8, 9], имеющей второй порядок аппроксимации на гладких решениях. Показано, что при расчете задачи о разрушении плотины (т. е. задачи о распаде начального разрыва глубин в теории мелкой воды [19, 20]) с образованием прерывной волны и волны понижения, эта схема имеет лишь первый порядок слабой сходимости, который приводит к соответствующему снижению порядка сильной сходимости в гладких частях точного решения за фронтом прерывной волны. Отметим, что этот результат опровергает достаточно широко распространенное мнение о том, что монотонные разностные схемы повышенного порядка аппроксимации на гладких решениях сохраняют соответствующий повышенный порядок локальной сходимости в гладких частях разрывных решений квазилинейных гиперболических систем законов сохранения.

2. Определение слабой сходимости и описание метода, позволяющего оценивать ее порядок

Рассмотрим задачу Коши для гиперболической системы квазилинейных законов сохранения [2, 3, 20]

щ + /

где и(г,х) — искомая, а ио(х) — заданная кусочно-непрерывные вектор-функции, зависящие от т компонент; f (и) — заданная гладкая функция потока, зависящая от т компо-

нент. Будем считать, что задача (1) имеет единственное устойчивое обобщенное решение и(г,х). Численное решение этой задачи будем искать по явной двухслойной по времени и симметричной по пространству консервативной разностной схеме

г’П+1 = 0 — ^га

где

/0+1/2 = /0+1,---0), 1

П— 1

о = у(ьп^к), го = 0, Ьп = ^2,Т*, = Тп/к,

г=0

тп — шаг схемы по времени на п-м временном слое гп, к — постоянный шаг схемы по пространству, / — непрерывная функция численного потока (конкретный ее вид будет приведен ниже).

Временной шаг Тп в схеме (2) будем определять из условия устойчивости Куранта

Тп < Тптах = г  к/ тах |аг(^П)|, г,0

где г = 0.5 — коэффициента запаса; аг(и), г = 1, ... , т — собственные значения матрицы Якоби А(и) = /и системы (1). При этом, если точное решение и(г, х) необходимо вычислить

на момент времени Т > 0, то временной шаг тп будем выбирать согласно следующему правилу:

' ттх, если Ьп + т”ах < Т,

Т - Ьп, если Ьп + т”ах > Т.

Получаемое таким образом разностное решение У/(Т,]Н) = у^, где ЬN = Т, путем линейной интерполяции доопределим до кусочно линейной функции

ун

заданной при Ух е И. Для введения понятия слабой сходимости разностного решения У/(Т,х) к точному решению и(Т,х) аппроксимируемой задачи (1) зафиксируем число а е И, и зададим следующие интегралы:

X X

иа

а а

Определение. Разностное решение у/(Т,х) слабо сходится с г-м (0 < г < 2) порядком на отрезке [а,х\ С И к точному решению и(Т,х), если

У£

где С — ограниченный векторный функционал, не зависящий от Н.

Ограничение г < 2, входящее в условие данного определения, связано с тем, что интеграл У/ (с учетом линейной интерполяции (4)) вычисляется по формуле трапеций, которая сама имеет лишь второй порядок точности на гладких функциях. Поскольку в данной работе изучаются ТУБ-схемы второго порядка аппроксимации на гладких решениях, то указанное ограничение является вполне допустимым. Причем, если отрезок интегрирования [а, х\ не содержит особенностей точного решения, то на нем, очевидно, порядок слабой сходимости будет совпадать с порядком локальной сходимости на гладких решениях, то есть будет вторым. Если же отрезок [а, х\ содержит внутри себя разрывы точного решения, то порядок слабой сходимости на нем, в общем случае, будет ниже второго, что связано с тем, что ТУБ-схемы, как показано в [17], не имеют второго порядка слабой аппроксимации на разрывных функциях. При этом корректность определения на таких отрезках следует из того, что интеграл У/ удовлетворяет интегро-интерполяционной формуле

УГЛ

где /п = ^Л+1/2 + Л-1/2^ /2, которая имеет место в силу линейности интерполяции (4) и консервативности схемы (2) и благодаря которой, как известно [1-3], предельные разрывные решения схемы (2) являются обобщенными решениями системы (1). Отметим, что, если численное интегрирование разностного решения проводить не на основе формулы трапеций (4), (5), а при помощи квадратурной формулы более высокого порядка точности на гладких функциях (например, формулы Симпсона), то тогда соотношение (7) будет нарушено и в результате определение слабой сходимости станет некорректным на разрывных решениях.

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

Предположим, что схема (2) удовлетворяет разностному аналогу энтропийного неравенства [15, 21], при помощи которого [3, 20] определяются устойчивые слабые решения системы (1) (как показано в [8], приводимая ниже ТУО-схема этому условию удовлетворяет). Тогда предельное разрывное решение разностной задачи Коши (2) будет устойчивым слабым решением дифференциальной задачи Коши (1). В этом случае для приближенного определения порядка слабой сходимости (применяя известное правило Рунге) достаточно провести 3 расчета задачи (2) с достаточно малыми пространственными шагами Н1 = Н, Н2 = Н/2 и Н3 = Н/4 . Для каждого из них выполнено соотношение (6), то есть

У£ - иа = СНГ + о(Нг), г = 1, 2, 3.

Отсюда следуют формулы

Ы = К - У£+11 = \СЩ - Н+!) + о

Предположив, что \С\ > £ > 0 (где е ^ 1 — число, определяемое машинной точностью

проведения численных расчетов) и взяв отношение равенств (8) с точностью до о(Нг), получим

^ = К-Щ = 1-

щ - Н3

2

Далее будут приведены результаты вычисления погрешностей (8) и порядка точности г для численного расчета по ТУБ-схеме [8] задачи разрушения плотины в русле непостоянной глубины.

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

В качестве конкретной разностной схемы (2) выберем одну из ТУБ-схем, предложенных Хартеном в [8]. Численный поток (3) в этой схеме вычислятся по следующим формулам, при записи которых для кратности опущен индекс временного слоя п:

Л+1/2 = 2

где

^ т

! (у3 ) + ! (у3 + 1) - Д X/ вк+1/2^-) +1/2

к=1

в)+1/2 = Я (.+1/2 + 1к]+1/^ а)+1/2 - (9к) + 9к)+1),

Я(х) = {х2/(4£)+ е, \х\ < 2е; \х\, \х\ > 2е}, е = 0.1,

1]+1/2 = {(9j+1 9] )/а]+1/2, а]+1/2 = 0; 0, а]+1/2 = 0}

9) = .+1/2 тах I0, т1п

9.+1 = 2 [Я(^.+1/2) (у3+1/2) ] а]+1/2,

.+ 1/2 = Дак (vj+1/2), а)+1/2 = ^]+1/2(у.+1 - ), Щ+1/2 = Кк (vj+1/2), ^к+1/2 = ^ (vj+1/2), ^.+1/2 = (+ У.+1 )/2,

ак — к-е собственное значение матрицы Якоби ^, Кк и Ьк — соответствующие ему правый и левый собственные векторы. Применим эту схему для численного решения уравнений Сен-Венана теории мелкой воды [19], которые в случае горизонтального прямоугольного русла при отсутствии трения эквивалентны уравнениям политропной газовой динамики с показателем политропы 7 =2. Векторная форма записи (1) указанных уравнений имеет вид

Щ + f (п)х = 0,

и =( М, 1(п)=( 9 V

\я) \92/И + 9Н /2 у

где Н, 9 — глубина и расход потока, 9 — ускорение силы тяжести (в расчетах полагалось, что 9 = 10). На рис. 1-3 на момент времени Т = 0.2 показаны результаты расчета задачи Коши со следующими разрывными начальными данными:

Н

9

В этих расчетах полагалось, что

Н1 = Н = 0.1, Н2 = Н/2 = 0.05, Нз = Н/4 = 0.025.

На рис. 1 кружками показаны результаты расчета на основной сетке с шагом Н1, сплошной линией — результаты расчета на сетке с шагом Н3. На рис. 1, а пунктирной линией показано начальное положение уровня воды, задаваемое графиком функции (11). На рис. 2 приведены соответствующие данной задаче графики функций (8), (9). В интегралах (5), при помощи которых вычисляются эти функции, полагалось, что а = 12.

Из рис. 2 следует, что при х > 6, то есть когда отрезок интегрирования [а,х\ целиком лежит в гладкой части решения перед фронтом ударной волны, порядок слабой сходимости г(х) = 2. Если же х < 5 и отрезок интегрирования [а,х\ 7 содержит внутри себя

фронт прерывной волны, то тогда порядок слабой сходимости г(х) 1 и, следователь-

но, ТУБ-схема (2), (14) приближенно имеет лишь первый порядок точности при расчете задачи разрушения плотины (11), (12).

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

г0(х) = {Я(х), К(х) < 3; 3, К(х) > 3},

где

Я

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

Рис. 1. Сравнение точного и численного решений. Поясн. см. в тексте.

порядок локальной сходимости. При х > 6.5 порядок го(х) = 2, это объясняется тем, что в этой области решение определяется характеристиками, приходящими с той части оси Ох, на которой функция начальных данных (11) является гладкой. При х £ [4.3, 5.3] — го(х) 1.3 в связи с тем, что в этой области решение определяется как характеристикой,

приходящей с оси Ох, так и характеристикой, приходящей с фронта ударной волны. Наконец, при х £ [2, 3], то есть внутри центрированной волны разрежения, порядок локальной сходимости оказывается самым низким го (х) < 1, поскольку в этой области одна из характеристик, определяющих решение, приходит из точки х = 3.5 разрыва начальных данных

(11).

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

Рис. 3. Порядки сильной сходимости.

уже отмечалось, это связано с тем, что TVD-схема, в общем случае, приближенно имеет лишь первый порядок слабой сходимости (см. рис. 2, б). Последнее полностью согласуется с тем, что она (имея второй порядок аппроксимации на гладких решениях [8]) имеет лишь минимальный порядок o(1) слабой аппроксимации на разрывных функциях [17].

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

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

[1] Годунов С. К., РЯБЕНЬКИЙ В. С. Разностные схемы. Наука, М., 1973.

[2] РОЖДЕСТВЕНСКИЙ Б. Л., ЯНЕНКО Н. Н. Системы квазилинейных уравнений. Наука, М., 1978.

[3] Lax P. D. Hyperbolic systems of conservation laws and the mathematical theory of shock

waves. Soc. Industr. and Appl. Math., Philadelphia, 1972.

[4] ГОЛЬДИН В. Я., КАЛИТКИН Н. Н., ШИТОВА Т. В. Нелинейные разностные схемы для гиперболических уравнений. Журн. вычисл. матем. и матем. физ., 5, №5, 1965, 938-944.

[5] РУСАНОВ В. В. Разностные схемы третьего порядка точности для сквозного расчета разрывных решений. Докл. АН СССР, 180, №6, 1968, 1303-1305.

[6] КОЛГАН В. П. Применение принципа минимальных значений производной к построению конечно-разностных схем для расчета разрывных решений газовой динамики. Уч. зап. ЦАГИ, 3, №6, 1972, 68-78.

[7] ХОЛОДОВ А. С. О построении разностных схем повышенного порядка точности для

уравнений гиперболического типа. Журн. вычисл. матем. и матем. физ., 10, №6,

1980, 1601-1620.

[8] HARTEN A. High resolution schemes for hyperbolic conservation laws. J. Comp. Phys., 49, 1983, 357-393.

[9] HARTEN A. On a class of high resolution total-variation-stable finite-difference schemes. SIAM. J. Numer. Anal, 21, №1, 1984, 1-23.

[10] ТОЛСТЫХ А. И. Компактные разностные схемы и их применение в задачах аэрогидродинамики. Наука, М., 1990.

[11] ПИНЧУКОВ В. И. О построении монотонных схем типа предиктор—корректор произвольного порядка аппроксимации. Матем. моделирование, 3, №9, 1991, 95-103.

[12] САМАРСКИЙ А. А., ЛАЗАРЕВ Р. Д., МАКАРОВ В. Л. Разностные схемы для дифференциальных уравнений с обобщенными решениями. Высш. школа, М., 1987.

[13] Lax P., Wendroff B. Systems of conservation laws. Comm. Pure and Appl. Math., 13, 1960, 217-237.

[14] ГОДУНОВ С. К. Оценки невязок для приближенных решений простейший уравнений газовой динамики. Журн. вычисл. матем. и матем. физ., 1, №4, 1961, 622-637.

[15] HARTEN A., HYMAN J. M., Lax P. D. On finite-difference approximations and entropy conditions for shocks. Comm. Pure and Appl. Math., 29, 1976, 297-322.

[16] Roland J. Dl Perna. Finite difference schemes for conservation laws. Ibid., 25, 1982, 379-450.

[17] ОСТАПЕНКО В. В. Об аппроксимации законов сохранения разностными схемами сквозного счета. Журн. вычисл. матем. и матем. физ., 30, №9, 1990, 1405-1417.

[18] ОСТАПЕНКО В. В. Аппроксимация законов сохранения на неоднородной разностной сетке. Там же, 33, №11, 1993, 1663-1680.

[19] ВОЕВОДИН А. Ф., ШУГРИН С. М. Методы решения одномерных эволюционных систем. Наука, Новосибирск, 1993.

[20] STOKER J. J. Water waves. Interscience publ. Inc., N. Y., 1957.

[21] HARTEN A., Lax P. D. A random choice finite difference scheme for hyperbolic conservation lawes. SIAM, J. Numer. Anal., 18, No. 2, 1981, 289-315.

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