Схема предиктор-корректор, сохраняющая гидравлический скачок Текст научной статьи по специальности «Физика»

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

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

СХЕМА ПРЕДИКТОР-КОРРЕКТОР, СОХРАНЯЮЩАЯ ГИДРАВЛИЧЕСКИЙ СКАЧОК*

Ю.И. Шокин, Г. С. Хакимзянов Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: shokin@ict.nsc.ru, khak@ict.nsc.ru

An explicit predictor-corrector difference scheme of the second order of approximation for solution of the one-dimensional nonlinear shallow water equations is considered, which yields nonoscillatory profiles of numerical solutions under the special choice of approximating viscosity. It is shown that the considered scheme preserves a stationary hydraulic jump for calculations using a uniform grid, and it preserves a moving jump when adaptive grids are used.

Введение

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

В настоящей работе рассматривается явная схема предиктор-корректор [5] для уравнений мелкой воды. В экспериментах [8] на одномерных задачах газовой динамики установлено, что в окрестности фронта ударной волны схема предиктор-корректор генерирует небольшие осцилляции численных профилей и дает приемлемую зону размазывания скачка. Схема имеет свободный параметр в, от которого зависят многие ее свойства. Например, в случае уравнения переноса с постоянным коэффициентом а

ut + aux = 0

* Работа выполнена при финансовой поддержке Программы интеграционных фундаментальных исследований СО РАН (гранты № 2006.113 и № 2006.2.12) и Программы поддержки ведущих научных школ РФ (грант № НШ-9886.2006.9).

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

явная схема предиктор-корректор состоит из шага предиктор

1

//+1/2 2 ^+1 + ^

+ а

т

¿+1/2

£П _ £П

■Ъ+1

к

(1)

аппроксимирующего уравнение

и шага корректор

Щ_+ ¿¿+1/2

/ + а/; /

* X*

- Л-

¿-1/2

0,

(2)

/

т

*

¿+1/2

т К

шаг по времени; К — шаг равномерной сетки с узлами x¿ = jh^, т¿* (1 + 0™+1/2)> ^/+1/2 — параметр схемы — сеточная функция, которая может меняться от

узла к узлу и от одного временного слоя к другому. При в = 1/Сг — 1 схема предиктор-корректор совпадает с противопоточной схемой, которая при условии

Сг = |а|ж < 1, ж = т/К = сопв^

(3)

устойчива, сохраняет монотонность численного решения, но является схемой лишь первого порядка аппроксимации. При в = 1/Сг2 — 1 схема предиктор-корректор переходит в схему Лакса, Если в = —1, то получаем неустойчивую схему второго порядка аппроксимации по К, При в = 0 схема предиктор-корректор переходит в схему Лакса — Вендроффа, которая имеет второй порядок аппроксимации по т и К, но не относится к классу схем, сохраняющих монотонность численного решения. Наконец, если в = О(К), то схема (1), (2) также имеет второй порядок аппроксимации и может при специальном выборе параметра в

В статье [9] на основе анализа дисперсии дифференциального приближения [10] схемы (1),

в

в

¿+1/2

0 их^'+1/2 < ^,¿ + 1/2-8

в0 (1 — £¿+1/2) при их^'+1/2 > ^,¿ + 1/2-8

их^+1/2 ^,¿+1/2-8 > 0, их^+1/2 ^,¿+1/2-8 > 0,

(4)

при ^,¿+1/2 ^,¿+1/2-8

-8 < 0.

Здесь 5 = ^па в0 — некоторая положительная постоянная; = (^¿+1 — ц)/К;

^¿+1/2 = их^+1/2-8/ихз-+1/2, здесь и ниже верхний индекс п опущен. При в0 = 1/Сг — 1 и условии (3) схема предиктор-корректор (1), (2) будет сохранять монотонность численного решения. Аналогичным свойством обладает и схема, аппроксимирующая нелинейное скалярное уравнение,

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

0

*

0

в

0

1. Схема предиктор-корректор на равномерной сетке для уравнений мелкой воды

Рассмотрим схему предиктор-корректор на равномерной сетке для системы уравнений мелкой воды с одной пространственной переменной

ди д ^ ^

Здесь

Ни ) ; £(и) = ( Ни2 + #2/2) ; С = ( ЯК*

где £ — время; и(х, £) — скорость; Н = п + К — полная глубина; п(х, £) — возвышение поверхности над невозмущенным уровнем; у = —К(х) — функция, задающая дно бассейна. Уравнения (5), дополненные начальными и краевыми условиями, решаются при £ > 0 в области О = (0, I).

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

(5):

дд

где А = М/ди — матрица Якоби, Отличительной чертой предлагаемой схемы является то, что численные потоки через границы ячеек вычисляются на основе уравнения для вектора потоков £

дд

+ = (7)

А

На шаге предиктор аппроксимируется система (7)

-^-+ (тг2?Л(ЛР - £С));+1/2 = О, (8)

а на шаге корректор — система (5)

„ГС+1 „п **

% - % + 5+1/2 ~ 5-1/2 = _ (9)

т К

при этом

А+1/2 = = (^¿3+1/2,

— ип01п 4- Ип Оч1п

и3 из + 1 + Н 3+1/2 2и3+1/2

Ли+1/2 0 \ / 1 + 0)+1/2 0

п

3+1/2 0 \ I 1 + ^3+1/2

Л™ , „ = I , vn

3+1/2 = | , ^3+1/2 -

0 Л2,]+1/2 ) V 0 1 + 3+1/2

%+1/2 = (и3-+1 + и^) /2, Лк (к =1, 2) — собственные значения матрицы Ап+1/2,

(А1,2)я-1/2 = ^+1/2 Т (Ю)

0

1

4 I 1 V л л ! -1 1

Гп _ 4 - 2 ~ 1

1-1-1 /9 --, чо I I 1 /V.- I

— (Л2_Л1,Ч_Л1 ^^ 4

элементы 0^+1/2 диагональной матрицы Р задаются по формуле, аналогичной (4)

0 при

0к+1/2 = \ 0к .+1/2(1 - Ск+1/2) ПРИ

9.?+1/2 <

9.?+1/2 >

^5+1/2 ^+1/2-^ > 0, ^?+1/2 &1/2-в* > 0 (И)

в0.+1/2 ПРИ .1/2 .1/2-** < 0

В КОТОрОЙ в к ^ЩП ЛП;.+1/^ 9^+1 /2 = Ж (в5С5+1/2 = .+1/2-**Рк,3'+1/2

компоненты вектора Р^+1/2 = )™+1/2> их,^+1/2 = (Чж — %)/Ь,

&о,з+1/2 = т;--1) 1/2 = |Л^+1/2|ж,

о \ _ М^я-О - М^-О

/о — ттп. 7_ 1 — пи 1 —

4+1/2 ^ Щ+1/2кх,]+1/2 ) ' ^ ^ ) ' 2Л

Видно, что для реализации шага корректор (9) требуются не только потоки Г^+1/2, но и предикторные величины #*, определенные в целых узлах и для вычисления которых требуются свои разностные уравнения. Они получаются [11] на основе аппроксимации в целых узлах недивергентной системы (6),

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

тт( Г #1, если х < хь(¿), , . Г С/1, если х < хь(£), , л

Н(х, ¿) = < тт м(х,£) = < „ (12)

[ #2, если х > хь(¿), [ и2, если х > хь(¿),

где Н1 > Н2, хь (¿) = х0 + ^ — скорость движения скачка,

= + <13>

На скачке выполняется соотношение

W (#2 — #1 ) = #2^2 — #1^1. (14)

Скачок является стационарным, если W = 0, В этом случае

#2^2 = #1^1. (15)

Лемма 1. Схема (8), (9), (11) сохраняет стационарный скачок (12).

Доказательство, Пусть = и(х., ¿п), = # (х., ¿п), W = 0 и н а п-м слое по времени скачок располагается между узлами х.0 и х.0+1; т, е, х.0 < х0 < х^0+1. Тогда Р™+1/2 = 0 при ] = Из формул (10) и (13) следует, что

= Ъ + ^ = й' (16)

п

0.0 -0.5 -1.0

100

200

и

1.0

0.5

0.0

100

200

Рис. 1. Графики свободной поверхности (а) и скорости (б). Сетка равномерная, И = 1, И2 = 0.01.

Т-е- Л^0+1/2

0, Учитывая равенство (15), получаем, что (ЛР)™+1/2 = 0, Тогда из уравнений шага предиктор (8) следует, что при всех ] потоки вычисляются по одной и той же формуле Ц+1/2 = (+ Р") /2, вследствие чего будем иметь

"+1 з

и

1 9 (К,3+1/2 + ^-1/2) '

В силу равенства

получаем

х,3+1/2 = Л3+1/2их,3+1/2

(17)

п+1 = п _

о з 2

(^ЛР )"+1/2 + П-

3-1/2

Отсюда с учетом равенства (ЛР)™±1/2 = 0 будет следовать, что и^1 = и", т. е. гидравлический скачок остается на месте и не размазывается рассматриваемой схемой.

Отметим, что указанное в лемме 1 свойство схемы остается справедливым и при других аппроксимациях матрицы Л, удовлетворяющих условию (17), например, когда в качестве Л"+1/2 берется матрица Еое [12],

Разумеется, что схема на равномерной сетке не сохраняет движущийся скачок. Для иллюстрации свойств схемы на нестационарных решениях рассмотрим задачу о разрушении плотины в случае плоского горизонтального дна К(х) = 1, Предполагается, что в начальный момент времени жидкость покоится, а полная глубина имеет постоянное значение Н1 слева от некоторой точки хо и значение Н2 справа от нее, 0 < х0 <1. На рис, 1 представлены результаты численного решения (штриховые линии) на равномерной сетке с числом узлов N = 100 в момент времени £ = 50 при I = 200, х0 = 100, Видно, что схема дает неосциллирующие профили, близкие к точным [13] (сплошные линии) и хорошую разрешимость сильного и слабых разрывов.

0

0

2. Схема предиктор-корректор на подвижной сетке

Координаты узлов подвижной неравномерной сетки обозначим через х". Схема строится на равномерной неподвижной сетке (^состоящей из узлов д3- = ]К {] = 0,

с" подвижной сетки равенствами х" = х(дз,

К = 1/^, которые связаны с узлами х" подвижной сетки равен ствами х" = х(д3-, £п), где

х = х(д, £) — некоторое невырожденное преобразование координат, которое в каждый момент времени £ взаимно-однозначно отображает единичный отрезок < = [0, 1] на О, причем х(0, £) = 0 х(1, £) = I.

Использование повой сетки < фактически означает переход в новую систему координат (д, £) [14]. Запишем в этой системе уравнения мелкой воды (5)

(/VЬ + ^ - х^)я = С=( Д) . (18)

Здесь £) = и(х(д, £), £). Аналогами уравнений (6) и (7) будут в данном случае следующие уравнения:

Ъ + хг£) V, = Ь + -^Л (Л - V, = ^

Схема предиктор-корректор состоит из шага предиктор

^•+1/2-^(^+1 + ^) /1

+ (-7гХ>Л(ЛР-£С))П =0 (19)

Т/2 V/ / 5+1/2

и шага корректор

(^)П+1 - (/у)" (^7+1/2 - (х^)П+1/^ - (^-1/2 - (х^)П-1/2)

где

+ ^--= С*, (20)

V™. ^"

д п _ дп _ п с -рп _ пп п п _ 5++_5_

Л1+1/2 - Л1+1/2 Х^+1/2С, 1/2 - 5+1/2 5,5+1/21 ^,5+1/2 - ^ )

г» _( о \ / о \ _ у?+1 + у?

5+1/2 - I ттп ип , 5 I I ' 5+1/2 - 9 )

\ Пз+1/2'\,з + 1/2 ) \ пз Я?,5 ) 2

„п | ™п „,п

„п _ Н,з "I" „ _ ^_, ,

3+1/2 — 9 ) Н,з — _ 1

Т

хп хп Тп + 7п I 1 х„- ______/9 Т О •_

тп _ „п _ Х5'+1 5 тп _ 5+1/2 ^ 3 — 1/2 , ^

5+1/2 - 5 + 1/2 ~ ^ ' 3 ~ 2 '

Функции ^/2 вычисляются по той же формуле (11), что использовалась в схеме предиктор-корректор на равномерной сетке. Только теперь учитывается подвижность сетки:

3+1/2 = |5п,5+1/2| (1 - Сгк,з+1/^ рП,3+1 /2' Сгк,5+1/2

ж|5п,5+1/2|

к,3-\-1/'А\ ~ ■ К,ЗТ1-! ±} г К,3 + 1/Л - ' К,31-1/* тп

^5+1/2

Здесь р" 5+1/2 — компоненты век тора Р"+1/^ 5П5+1/2 — диагональные элементы матрицы

лП+1/2; ^ = 8§п 5П, 5+1 /2 ? к = 12

Рассмотрим движущийся гидравлический скачок (12), (13) и предположим, что в начальный момент времени £ = 0 от располагался между узлами х0о и х°0+1, т.е. х0о < х0 < х0о+1, Покажем, что схема предиктор-корректор точно описывает этот скачок, если на каждом слое по времени выполняется равенство

х"50+1/2 = Ж. (23)

а б в

Рис. 2. Графики свободной поверхности (а), скорости (б) и траекторий узлов адаптивной сетки (в). И = 1, И = 0.01.

Для справедливости условия (23) достаточно, например, выполнения при всех п равенств х— = х—+1 = т. е. достаточно, чтобы два ближайших к разрыву узла подвижной сетки двигались с той же скоростью, что и разрыв.

Лемма 2. При условии (23) схема предиктор-корректор на подвижной сетке сохраняет движущийся скачок (12), (13).

Доказательство. В силу равенств (14), (16) и условия (23) получаем, что (ЛР)™+1/2 = 0, В других узлах равенство (ЛР)"+1/2 = 0 также выполняется, поскольку Р™+1/2 = 0 при у = j0. Поэтому из уравнения шага предиктор (19) будет следовать, что ^?+1/2 = (-+1 + —О/2- Тогда на шаге корректор (20) получим

(^)Г1 = (»)? -\ (ъ+1/а + +

Используя равенства

^¿+1/2 = ^+1/2^+1/2 = )"+1/2 = (ЯЛР)-+1/2 + (xtVq—+1/2 = (х^q —+1/2 ,

а также формулу дифференцирования произведения двух функций

(а^+1/2 ~ (ж^'-1/2 = ед^ + \ [(хъ%+1/2 + ,

запишем равенство (^)™+1 = (— + т(xtУ^-) V™, которое в силу тождества /™+1 = — + т(xt)п - означает, что = V™, т.е. на подвижной сетке, удовлетворяющей условию (23), схема предиктор-корректор не размазывает движущийся скачок (12), (13).

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

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

[1| Марчук Г.И., Яненко H.H. Применение метода расщепления (дробных шагов) для решения задач математической физики / / Некоторые вопросы вычислительной и прикладной математики. Новосибирск: Наука, 1966. С. 5-22.

[2] Тушева Л.А., Шокин Ю.И., Яненко H.H. О построении разностных схем повышенного порядка аппроксимации на основе дифференциальных следствий // Некоторые проблемы вычислительной и прикладной математики. Новосибирск: Наука, 1975. С. 184-191.

[3] Яненко H.H., Шокин К).И. О групповой классификации разностных схем для системы уравнений газовой динамики // Тр. МИЛИ СССР. 1973. Т. 122. С. 85-96.

[4] ГОДУНОВ С.К. Разностный метод численного расчета разрывных решений уравнений гидродинамики // Математический сборник. 1959. Т. 47, вып. 3. С. 271-306.

[5] Яненко H.H., Шокин К).И. Об аппроксимационной вязкости разностных схем // Докл. АН СССР. 1968. Т. 182, № 2. С. 280-281.

[6] Яненко H.H., Яушев И.К. Об одной абсолютно устойчивой схеме интегрирования уравнений гидродинамики // Тр. МИЛИ СССР. 1966. Т. 74. С. 141-146.

[7] Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001.

[8] Яушев И.К. О численном расчете нестационарных течений газа в одномерном приближении в каналах со скачком площади сечения // Изв. СО АН СССР. Сер. Техн. науки. 1967. № 8, вып. 2. С. 39-48.

[9] Сергеева Ю.В., Хакимзянов Г.С. Об использовании дифференциального приближения при построении монотонных схем // Вычисл. технологии. 2004. Т. 9, спецвыпуск. С. 139-149.

[10] Шокин К).И., Яненко H.H. Метод дифференциального приближения. Применение к газовой динамике. Новосибирск: Наука, 1985.

[11] Шокин К).И., Сергеева Ю.В., Хакимзянов Г.С. О монотонизации явной схемы предиктор-корректор // Вест. КазНУ. Сер. Математика, механика, информатика. 2005. № 2. С. 103-114.

[12] Roe P.L. Approximate Riemann problem solvers, parameter vectors, and difference schemes // J. Comp. Phys. 1981. Vol. 43, N 2. P. 357-372.

[13J Численное моделирование течений жидкости с поверхностными волнами /P.C. Хакимзянов, Ю.И. Шокин, В.Б. Барахнин, Н.Ю. Шокина. Новосибирск: Изд-во СО РАН, 2001.

[14] Яненко H.H., Шокин Ю.И., Урусов А.И. О разностных схемах в произвольной криволинейной системе координат // Докл. АН СССР. 1978. Т. 242, № 3. С. 552-555.

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