Разностные схемы

Рисунок 11.18. Решение уравнений (11.54) с граничными условиями (11.53) в момент времени t=1 (l=1, c=1, ρ=1); - точное решение, - приближённое решение

(h=0.05, τ=0.0333, r=2/3).

Это решение подтверждает результаты проведённого выше теоретического анализа: амплитуда волны уменьшается (схема (11.56) диссипативна) и происходит запаздывание волны, так как скорость распространения всех гармоник меньше чем групповая скорость, определяемая системой уравнений (11.54) (смотри рисунок 11.17).

________________________________________________________________________________

11.3.4. Характеристическая форма уравнений акустики.

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

∂∂ut + ρ1 ∂∂xp = 0 , .

∂∂pt + ρc2 ∂∂ux = 0.

Умножим первое уравнение на сρ, а затем получим два новых уравнения путём сложения этих уравнений и вычитания первого уравнения из второго:

51

	∂		(p + cρ u) + c			∂		(p + cρ u) = 0 ,
∂t					∂x
							,		(11.59).
	∂					∂
				(p −cρ u)−c				(p −cρ u) = 0.
	∂t					∂x

В результате мы получим характеристическую форму уравнений (11.54). Она представляет собой два уравнения переноса для величин R=p+сρ u и S=p-сρ u, которые называются инвариантами Римана. Решение исходной системы может быть легко выражено через эти

инварианты:

u = 2c1ρ(R −S), . p = 21 (R + S).

Решение этой системы можно записать в виде R(x,t)=R(x-ct) и S(x,t)=S(x+ct). Решение такого вида говорит о том, что инвариант R сохраняется вдоль характеристики dx/dt=c, то есть R=const при x-ct=const, и, соответственно, инвариант S сохраняется вдоль характеристики dx/dt= -c: S=const при x+ct=const.

Используя эти свойства инвариантов Римана, построим решение задачи о распаде разрыва. Пусть начальное условие для задачи (11.4) имеет следующий вид:

u		, x ≤ 0
l			,
u(x, 0) =		, x > 0
u	r
		, x ≤ 0 .
p
	l		.
p(x, 0) =		, x > 0
p
	r

При t>0, эта конфигурация распадается на две волны, распространяющиеся влево и вправо от точки x=0 (рисунок 11.19).

Обозначим решение в точке x=0 как u0 и p0. Тогда, учитывая свойства инвариантов

Римана, в точке x=0 можно записать

p0 +cρ u0 = R(−ct) = pl +cρ ul ,

p0 −cρ u0 = S(ct) = pr −cρ ur .

В результате, мы приходим к решению

p0	=	1	(pl + pr ) +	1 cρ(ul −ur ),	(11.60)
		2		2

u0 = 21 (ul + ur )+ 2c1ρ(pl − pr ).

52

Рисунок 11.19. Задача о распаде разрыва; - начальное состояние; - состояние в некоторый момент времени t>0.

11.3.5. Недиссипативная схема Ф. Роу.

Схема, предложенная Ф. Роу, основывается на характеристической форме (11.59). Так как эта форма представляет собой систему из двух уравнений переноса, то имеет смысл рассмотреть сначала построение схемы для уравнения

∂∂ut + c ∂∂ux = 0 .

Как обычно, введём разностную сетку с узлами (xn,tk) и шагами h и τ. Основная особенность этой схемы состоит в методе аппроксимации производной по времени: приближение для ∂u/∂t задаётся как полусумма разности вперёд по времени в точке xn и разности назад по времени в точке xn-1. Производная по x аппроксимируется разностью назад относительно точки xn. Шаблон этой схемы имеет вид, показанный на рис. 11.20, и разностная схема записывается как

1	uk+1		−uk	uk	−uk−1			uk	−uk
		n	n +	n−1	n−1		+ c	n	n−1	= 0 ,	(11.61)
2			τ		τ				h

n=1, …, N-1; k=1, 2, … .

53

Рисунок 11.20. Шаблон схемы (11.61).

В отличии от рассмотренных нами схем, схема (11.61) является трёхслойной (разностное уравнение связывает три момента времени).

Разностная схема для одномерного уравнения акустики строится на основе характеристической формы (11.59). Для аппроксимации первого уравнения применяется схема (11.61). Аппроксимация второго уравнения производится с помощью аналогичной схемы, учитывающей, что это уравнение описывает волну распространяющуюся влево. В результате, мы получим следующую разностную схему:

1

Rk+1

−Rk

Rk

−Rk−1

Rk −Rk

n

n +

n−1

n−1

+ c

n

n−1

= 0 ,

2

τ

τ

h

1

Sk+1

−Sk

Sk

−Sk−1

Sk

−Sk

n

n +

n+1

n+1

−c

n+1

n

= 0 ,

(11.62)

2

τ

τ

h

Rk

= pk

+ cρ uk

,

n

n

−cρ

n

,

Sk

= pk

uk

,

n

n

n

unk = 0.5(Rnk −Snk )/cρ, .

pnk = 0.5(Rnk + Snk ).

Для анализа свойств этой схемы, достаточно рассмотреть схему (11.61), так как разностные уравнения в (11.62) непосредственно не связаны друг с другом. Подставляя точное решение в (11.11) и используя выражения (3.14), получим

∂u

2

∂u

2

1

(xn ,tk ) +

1

τ

∂ u

(xn

,tk ) +

(xn−1,tk ) −

1

τ

∂ u

+

2

∂t

2

∂t

2

∂t

2

∂t

2

(xn−1,tk )

∂u

1

2

c

(xn

,tk ) −

h

∂ u

(xn

+O (h

2

+ τ

2

) = 0.

∂x

2

∂x

2

,tk )

54

Так как мы проводим анализ аппроксимации в точке (xn,tk), то нам ещё необходимо разложить величины определённые в точке (xn-1,tk) в ряд Тейлора в точке xn по степеням h. Тогда мы приходим к следующему дифференциальному приближению:

∂u

(x

,t

) + c ∂u

(x

,t ) =

∂t

n

k

∂x

n

k

1

∂2u

1

∂2u

2

2

h

(xn ,tk ) +

ch

∂x2 (xn ,tk ) +O (h

+ τ

).

2

∂x∂t

2

Из уравнения переноса следует, что ∂2u/∂x∂t=-c∂2u/∂x2. Тогда два первых слагаемых в правой части предыдущего выражения сокращаются и, следовательно, схема (11.12) имеет второй порядок аппроксимации по h и τ. Анализ устойчивости схемы (11.12) оставляется читателю в качестве упражнения.

Для анализа диссипации и дисперсии применим подход, рассмотренный в пункте

11.3.3. В результате получим следующее дисперсионное соотношение: sin(ωτ − 12 kh) = (1 −2r)sin(21 kh).

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

							1
cm =	dω	=	c		−	(1 −2r)cos(			kh)
								2
				1									.
						2	2				1
	dk		2r			1 −(1 −2r)	sin			(		kh)
											2

На рисунке 11.21 показана зависимость безразмерного параметра a=сm/c от параметра Куранта.

55

Рисунок 11.21. Зависимость параметра a=сm/c от параметра Куранта для схемы Ф. Роу; - безразмерная групповая скорость волны cg/c=1 для дифференциального уравнения; - параметр a при kh=π/2.

Рассмотрим теперь расчёт граничных условий. Пусть, например, задано граничное условие при x=0: u(0,t)=f(t). Вычисление этого условия не представляет особого труда:

uk+1	= f (t	k+1	) . Для вычисления давления в граничной точке мы используем схему (11.12).
0

Так как состояние в точке x=0 определяется волной распространяющейся влево, то следует

взять второе уравнение при n=0:

21τ (p0k+1 −cρ u0k+1 −S0k + S1k −S1k−1 )− hc (S1k −S0k ) = 0 .

Тогда второе граничное условие вычисляется по формуле

p0k+1 = cρ u0k+1 + (1 −2r)S0k −(1 −2r)S1k + S1k−1 .

Для вычисления граничных условий заданных на правой границе при x=l, применяется аналогичный подход. Только в этом случае, следует использовать уравнение для R инварианта при n=N.

________________________________________________________________________________

Пример 11.6 (схема Ф. Роу)

Как и прежде, рассмотрим задачу (11.54) с нулевыми начальными условиями и граничными условиями (11.53). Приближённое решение, полученное с помощью схемы (11.62), показано на рисунке 11.22.

56

un-1/2

Рисунок 11.22. Решение уравнений (11.54) с граничными условиями (11.53) в момент времени t=1 (l=1, c=1, ρ=1); - точное решение, - приближённое решение

(h=0.05, τ=0.0333, r=2/3).

________________________________________________________________________________

11.3.6. Метод С. К. Годунова.

До сих пор мы представляли приближённое решение как набор значений в узлах разностной сетки. Другой способ представления приближённого решения состоит в том, что распределение величин аппроксимируется ступенчатыми функциями, такими, что эти величины имеют постоянные значения на каждом интервале [xn,xn+1], как показано на рисунке 11.23. Значение некоторой величины на интервале [xn,xn+1] будет обозначаться как

и, для определённости, мы будем относить это значение к середине интервала. Далее мы рассмотрим построение разностных схем на основе такого представления приближённого решения.

Рисунок 11.23. Представление приближённого решения в виде кусочно-постоянной функции.

Для того чтобы построить основную схему, проинтегрируем уравнения (11.54) по площади разностной ячейки {xn≤x≤xn+1, tk≤t≤tk+1}:

∫ ∫			∂s dt dx + A∫ ∫				∂s dx dt = 0 .
xn	+1	tk +1		tk +1	xn	+1
			∂t
x	n	t		t	x	n	∂x
		k		k
				57

Вычисляя интегралы от производных, получим

xn +1	tk +1
∫ (s(x,tk+1) −s(x,tk ))dx + A∫ (s(xn+1,t) −s(xn ,t))dt = 0 .
xn	tk
Введём следующие обозначения:
xn +1
h1 ∫ s(x,tk )dx = snk +1/2 -среднее значение величины s на интервале [xn,xn+1],
xn
tk +1
∫ s(xn ,t)dt = fnk (s) = (fnk (u), fnk (p))T	- поток величины s через узел xn за время от tk до tk+1.
tk
С учётом этих обозначений, предыдущее соотношение принимает вид
(snk++11/2 −snk +1/ 2 )h + A(fnk+1(s) −fnk (s)) = o .		(11.63)

Самый простой способ вычисления потоков основывается на предположении, что величина

s(xn,t) постоянна на каждом интервале [tk,tk+1]

и мы обозначим её значение как

snk . Тогда

схема (11.63) записывается как

snk++11/2

−snk +1/ 2

+ A

snk

+1

−snk

= 0 ,

(11.64)

τ

h

или в развёрнутом виде

unk++1/1

2

−unk+1/2

+

pk

− pk

= 0 ,

n+1

n

τ

ρh

pnk++1/1

2 − pnk+1/ 2

uk

−uk

+ ρc2

n+1

n

= 0 ,

τ

h

n=1, …, N-1; k=0, 1, … .

Для того чтобы производить вычисления по схеме (11.64), нам необходимо иметь некоторый алгоритм вычисления значений величин в узлах xn,через средние значения этих величин в соседних ячейках. С. К. Годунов предложил использовать для этой цели решение задачи о распаде разрыва. Действительно, в некоторый момент времени tk приближённое решение имеет вид ступенчатой функции как показано на рисунке 11.23. При t>tk каждый разрыв распадается на две волны. Одна из этих волн движется в направлении возрастания координаты x, а вторая – в сторону её убывания. Между волнами, скажем в точке xn, u и p постоянны по времени до тех пор, пока волны от соседних разрывов не пройдут расстояние h, то есть в течении времени τ≤h/c (это и есть физическое обоснование устойчивости данной

58

схемы). Поэтому, при tk≤t≤tk+τ можно использовать решение задачи о распаде разрыва (11.60). Если обозначить

ul				k		ur				k
			un−1/2						un+1/2
						и
		=		k				=		k	,
p						p
	l						r
			pn−1/2						pn+1/ 2

то значения величин u и p в узле xn вычисляются по формулам

uk = 1 (uk− + uk+ )+ 1 (pk− − pk+ ),

n 2 n 1/2 n 1/2 2cρ n 1/2 n 1/2

(11.65)

pk = 1 (pk−1/2 + pk+1/2 )+ cρ(uk−1/ 2 −uk+1/2 ). n 2 n n 2 n n

Подставляя эти выражения в (11.64) мы получим схему С. К. Годунова для уравнений акустики:

unk++1/1

2 −unk+1/2

pnk

+3/2 − pnk−1/ 2

c

k

k

k

= −

+

(un+3 /2

−2un+1/2

+ un−1/ 2 ),

(11.66)

τ

2ρh

2h

pnk++1/1

2 − pnk+1/ 2

= −c2ρ

unk+3/ 2 −unk−1/2

+

c

(pnk+3/2 −2pnk+1/2 + pnk−1/ 2 ),

τ

2h

2h

n=1, … , N-1; k=0, 1, … .

Анализ аппроксимации проведём в точке (xn+1/2,tk). Подставляя точное решение в (11.66) и используя выражения (11.11), получим

∂u

(x

n+1/2

,t

) +

1 ∂p

(x

n+1/ 2

,t ) =

∂t

k

ρ ∂x

k

.

1

∂2u

1

∂2u

2

2

).

−

2

τ

∂t2

(xn+1/2,tk ) +

2 ch

∂x2 (xn+1/2,tk ) +O (h

+ τ

Приводится только первое уравнение из (11.66), так как для второго уравнения получается аналогичное выражение. Из уравнения (11.54) следует, что ∂2u/∂t2=c2∂2u/∂x2. Тогда предыдущее выражение преобразуется в

∂u	(x	n+1/ 2	,t ) +	1 ∂p	(x	n+1/2	,t ) =
∂t			k	ρ ∂x			k

1ch(1 −r) ∂2u (xn+1/2,tk ) +O (h2 + τ2 ).

2∂x2

иэто говорит о том, что схема (11.66) имеет только первый порядок по h и τ. Качественный анализ этого выражения показывает, что схема (11.66) диссипативна. Действительно, члены

вида μ∂2u/∂x2 представляют собой вязкие силы, поэтому с точностью до величин второго

59

порядка малости можно сказать, что схема (11.66) добавляет к исходному уравнению (11.54) «вязкость» 0.5ch(1-r)∂2u/∂x2. Такая «вязкость» (называемая аппроксимационной вязкостью) приводит к тому, что амплитуда волны уменьшается со временем при r<1.

Разностную схему (11.16) можно записать в виде

skn++11/ 2 = skn+1/ 2 − 2τh A(skn+3/ 2 −skn−1/2 )+ 2chτ A2 (sn+3 /2 −2skn+1/ 2 + skn−1/ 2 ).

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

λ(α) = 1 −2r sin2(α/2) ± ir sin(α).

Тогда из критерия фон Неймана следует, что схема (11.66) будет устойчива при r≤1. Коэффициент затухания и групповая скорость сеточного решения удовлетворяют

соотношениям

γτ = −21 ln (1 −r + r cos(kh))2 + r2 sin2(kh) ,

	c	=	dω0	= c	r + (1 −r)cos(kh)		.
			d k		(1 −r)2 + r2 + 2r(1 −r)cos(kh)
	m

На рисунках 11.24 и 11.25 представлены зависимости параметра диссипации d=γτ и групповой скорости от параметра Куранта.

Рисунок 11.24. Зависимость параметра диссипации d=γτ от параметра Куранта для схемы С. К. Годунова первого порядка (kh=π/2).

60