научная статья по теме TVD-СХЕМА ДЛЯ РАСЧЕТА ВОЛНОВЫХ ТЕЧЕНИЙ В ОТКРЫТЫХ РУСЛАХ Математика

Текст научной статьи на тему «TVD-СХЕМА ДЛЯ РАСЧЕТА ВОЛНОВЫХ ТЕЧЕНИЙ В ОТКРЫТЫХ РУСЛАХ»

ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2008, том 48, < 12, с. 2212-2224

УДК 519.634

TVD-CXEMA ДЛЯ РАСЧЕТА ВОЛНОВЫХ ТЕЧЕНИЙ В ОТКРЫТЫХ РУСЛАХ1

© 2008 г. М. В. Бунтина, В. В. Остапенко

(630090 Новосибирск, пр-т акад. Лаврентьева, 15, Ин-т гидродинамики СО РАН) e-mail: ostapenko@hydro.nsc.ru; ostapenko_w@ngs.ru Поступила в редакцию 09.11.2007 г.

Для уравнений первого приближения теории мелкой воды (уравнений Сен-Венана) разработана TVD-схема, предназначенная для сквозного расчета течений с прерывными волнами в открытых руслах. Предложенная схема использует специальную аппроксимацию недивергентной формы записи уравнения полного импульса, в которой отсутствуют интегралы, связанные с определением силы давления в поперечном сечении русла и силы реакции стенок русла. В стандартных консервативных разностных схемах на вычисление этих интегралов расходуется основная часть машинного времени. Тестовые расчеты показали, что предложенная схема передает соотношения на разрывах с точностью, достаточной для численного моделирования процесса распространения реальных прерывных волн. В качестве примера построены все качественно различные решения задачи о разрушении плотины в трапецеидальном русле, имеющем область сужения в нижнем бьефе. Библ. 23. Фиг. 8.

Ключевые слова: уравнения мелкой воды, прерывные волны, открытое русло, численная ТУБ-схема.

1. ВВЕДЕНИЕ

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

В настоящей работе для численного моделирования прерывных волн в открытых руслах на основе метода, предложенного в [7], разработана специальная модификация тУБ-схемы из [13], которая в случае призматического русла имеет второй порядок аппроксимации на гладких решениях. В этой схеме, в отличие от схемы, рассмотренной в [7], ширина размазывания фронта прерывной волны не зависит от ее интенсивности и составляет 3, 4 ячейки разностной сетки. Построенная схема является консервативной в случае прямоугольного русла и неконсервативной в случае русел другой формы. Несмотря на это тестовые расчеты показали, что данная схема передает условия Гюгонио на фронте прерывной волны, распространяющейся по треугольному,

1)Работа выполнена при финансовой поддержке гранта Президента Российской Федерации по государственной поддержке ведущих научных школ (< НШ-5873.2006.1), РФФИ (код проекта 04-01-00253) и в рамках Проекта фундаментальных исследований Президиума РАН < 16.2.

2212

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

2. СИСТЕМА УРАВНЕНИИ СЕН-ВЕНАНА

Дифференциальные уравнения теории мелкой воды (уравнения Сен-Венана) без учета влияния трения и уклона дна имеют вид (см. [1], [2])

fflf + qx = 0, (2.1)

qt + (q v + gP)x = gR, (2.2)

H

ю = J b (x,%) , (2.3)

0

HH

P = Jffl( x,%) d\ = J( H -b (x,%) d\, (2.4)

00 HH

R = Jfflx(x,%)d^ = J(H - %)bx(x,%)d^, (2.5)

00

где H = H(t, x) - уровень потока, отсчитываемый от дна, q = q(t, x) - расход воды в поперечном сечении потока, b = b(x, H) - ширина русла на уровне H, ю = ffl(x, H) - площадь поперечного сечения потока при уровне H, v = v(t, x) = q/ю - средняя по сечению скорость потока, pgP = pgP(x, H) -сила гидростатического давления воды в поперечном сечении потока, pgR = pgR(x, H) - сила реакции стенок русла, вызванная его непризматичностью, g - ускорение силы тяжести, p - плотность воды.

В случае призматического русла, для которого

bx = 0 ^ю = ю(HP = P(H), R = 0, (2.6)

система уравнений (2.1)-(2.5) с учетом замены площади поперечного сечения потока ю на плотность газа p эквивалентна уравнениям изоэнтропической газовой динамики (см. [3]) с уравнением состояния

p (p) = P( H (p)), (2.7)

где H(p) = H(ю) - функция, обратная к ю(Д). При степенной функции для ширины русла

„а + 1 а + 2

b(H) = Ha ^ ю(H) = У— ^ P(H) = --^----v ' v ' а +1 v ' (а +1 )(а + 2)

уравнение состояния (2.7) имеет вид

.у . 1 а+2 p = Ap , A = -, Y = --.

а--а----------а + 1

(а + 1 )а +1 (а + 2)

Отсюда следует, что при возрастании а от 0 до значение у убывает от 2 до 1. В частности, при а = 0, что соответствует прямоугольному руслу, получим классическую аналогию с газом, для которого у = 2, а при а = 1/2 - аналогию с одноатомным газом, для которого у = 5/3.

Уравнения (2.1) и (2.2), входящие в систему (2.1)-(2.5), представляют собой дифференциальную форму записи законов сохранения массы и полного импульса и поэтому дивергентная аппроксимация их левых частей является необходимым условием консервативности разностной схемы (см. [14, [15]), предназначенной для сквозного расчета прерывных волн. В случае реальных русел, для которых функция b = b(x, H) задается по результатам натурных измерений, это

приводит к необходимости численного определения интегралов (2.4) и (2.5), на что в конкретных расчетах уходит основная часть машинного времени. В то же время, использовав формулу

н

Рх = ю Нх +

^хД = ю нх + Я,

(2.8)

уравнение для полного импульса (2.2) можно переписать в следующей недивергентной форме:

Ъ + (^) х + £ю Нх = 0,

которая не содержит интегралов (2.4) и (2.5).

Запишем систему (2.1), (2.9) в векторной форме

и + f (и) х + Я ю ф (и) х = 0,

где

Применяя формулу

ю^ I ъ \ _ I 0

и = д, <•< и) = иу • Ф (и) = IН

н

юх = Ь(х, н)Нх + |Ьх(х,£)д^,

(2.9)

(2.10)

(2.11)

из которой следует, что

1

| ЬЛ

Н = т ю х -

V 0 >

векторное уравнение (2.10) можно представить в виде

иг + А (и )и х = Е (и),

где

А (и) =

2 2 -V с - V 2 V у

, Е (и) =

н

-| Ьхд^

здесь с = -У^ю/Ь - скорость распространения малых возмущений. Матрица А(и) имеет собственные значения

= V - с, А2 = V + с

и соответствующие им левые

I1 = к( А* 2, -1), I2 = К(-А,1,1)

и правые

г1 = к( )т, Г2 = к( )т

(2.12)

(2.13)

(2.14)

(2.15)

собственные векторы, в которых к = (2с)-1/2. Из (2.13) следует, что при Н > 0 система уравнений Сен-Венана является гиперболической при выполнении естественного условия Ь(Н) > 0 при Н > 0. Поскольку системы левых (2.14) и правых (2.15) собственных векторов являются биортогональ-ными, т.е. 1 Т' = 8у, то имеют место матричные соотношения

А = ^ Л 5 ^ А2 = ^ Л2

(2.16)

0

0

где

5 =

ч! у

, 5 1 = (г1, г2), Л =

Г А 0 л

V 0 А2У

(2.17)

которые используются ниже при построении ТУБ-схемы, аппроксимирующей систему (2.10), (2.11).

3. МОДИФИЦИРОВАННАЯ ТУБ-СХЕМА

Аппроксимируем систему уравнений (2.10), (2.11) явной двухслойной по времени разностной схемой

п +1 п „п „п

и 1 - и 1 * 1 + 1- * 1 -1 7лп птП

—-- + -2±1-—— + 1 = W1,

тп 2Л 6 1 1

(3.1)

которой

Й. = ю Ф-+1 - Ф- - 1 = "1 = ю1 2 Л =

Н1 +1 -Н

ю ;

1 -1

(3.2)

есть разностный оператор, дифференциальное представление (см. [16]) которого в случае призматического русла (2.6) является дивергентным с точностью до членов 0(Л4);

W1 = № 1+1/2-№

1 + 1/2 ~ " - - 1/2 Л

(3.3)

есть дивергентный разностный оператор, построенный путем ТУБ-модификации (см. [13]) схемы Лакса-Вендроффа (см. [14]). При записи схемы (3.1) использованы следующие сокращенные обозначения для сеточных функций:

ип = и(¿п, х1), х1 = ]Л, Хп = ^^ Тт,

т =1

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

а Л

а Л

/|л|т|^|тч , т\ т,

тахтах((А^^, |А2|^ ) тах(|V ^ | + с^ )

(3.4)

где а е (0, 1) - коэффициент запаса.

Для построения сеточной функции ю ^, входящей в разностный оператор (3.2), воспользуемся методом, предложенным в [7]. Рассмотрим симметричный разностный оператор

О 1 = ю1

Н1 + 1-Н1--2Л

(3.5)

который со вторым порядком аппроксимирует дифференциальный оператор юНх. В случае призматического русла (2.6), для которого с учетом (2.8) имеем юНх = Рх, дифференциальное представление (см. [16]) оператора (3.5) относительно центральной точки его шаблона х1 = х имеет вид

О = Рх + | ю Нххх + о (Л4).

(3.6)

Формальное повышение порядка аппроксимации этого оператора, связанное с аппроксимацией недивергентного оператора третьего порядка юНххх, приводит к появлению искусственной

п

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

юНххх = (юНхх - юхНх)х + юххНх

и симметричную разностную аппроксимацию

О * = ю + 1 - 22 ю + ю _ Н+1 - н - 1 ' = н2 2Н

дифференциального оператора юххНх, строим разностный операт

Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.

Показать целиком