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

Схема С. К. Годунова выгодно отличается от других схем простотой расчёта граничных условий. Например, если в точке x=0 задано условие u(0,t)=f(t), то скорость в граничном узле вычисляется как u0k = f (tk ) , и нам ещё необходимо определить давление в этом узле. Состояние в точке x=0 определяется волной движущейся в направлении убывания координаты x и зависит от значений u и p в ячейке, примыкающей к границе. Тогда в течение

времени tk≤t≤tk+τ условие сохранения инварианта S примет вид

p0k −cρu0k = p1/k 2 −cρu1/k 2 ,

откуда следует формула для вычисления давления в граничном узле. Другие типы граничных условий вычисляются аналогично.

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

________________________________________________________________________________

Пример 11.7 (схема С. К. Годунова первого порядка)

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

Этот результат показывает, что амплитуда приближённого решения значительно отличается от точного решения. Это объясняется тем, что схема С. К. Годунова при r=2/3 приводит к диссипации сеточного решения. К тому же для r=2/3, групповая скорость cm>c

61

(смотри рисунок 11.25) и при t=1 приближённое решение включает в себя отражённую от жёсткой стенки волну. В точном же решении отражение происходит при t>1.

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

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

________________________________________________________________________________

11.3.7. Разностные схемы для решения многомерных задач.

Построение разностных схем для волновых уравнений с двумя и тремя пространственными переменными очевидным образом более сложная задача, однако, используемые принципы идентичны применяемым в одномерном случае. Рассмотрим вначале двумерное уравнение колебаний в прямоугольной области

2				2			2
∂ u			2	∂ u			∂ u
		= c				+				,	(11.67)
	2				2			2
∂t				∂x			∂y

0≤x≤lx, 0≤y≤ly , 0≤t≤T

62

с некоторыми начальными и граничными условиями. Чтобы построить разностную схему для этого уравнения, поступим так же, как в одномерном случае. С этой целью построим разностную сетку с узлами (xn,ym,tk) и шаги этой сетки равны xn+1-xn=hx, ym+1-ym=hy и tk+1- tk=τ. После этого входящие в уравнение (11.67) частные производные заменяются соответствующими конечно-разностными соотношениями. В результате мы получим схему

k+1

k

k−1

k

k

k

k

k

k

un,m −2un,m + un,m

2

un

+1,m −2un,m + un−1,m

un,m+1 −2un,m + un,m−1

= c

+

, (11.68)

2

2

2

τ

h

h

x

y

n=1, …, N-1; m=1,…, M-1; k=1, 2, … ,

которая является непосредственным обобщением схемы (11.52). Исследование аппроксимации производится, так же как и в одномерном случае. В результате мы получим, что схема (11.68) аппроксимирует уравнение (11.67) с точностью до членов O(τ2+(hx)2+(hy)2). Процедура исследования устойчивости аналогична процедуре, которая использовалась в одномерном случае. Только в двумерном случае решение задаётся в виде

unk,m = λk exp(iα1n + iα2m).																			(11.69)
Подстановка этого решения в схему (11.68) даёт уравнение
λ −2λ + 1 = −4λc2τ2		sin2			(		1	α1 )−4λc2τ2				sin2	(		1	α2 ).
							2								2
h2											h2
x											y
Решив это уравнение мы найдём λ(α1,α2). Тогда из критерия фон Неймана следует, что
2			1						1							hxhy
		2
c	τ					+					≤1 или τ ≤							.	(11.70)
				2						2			2				2
									h
		h		x									c hx				+ hy
									y

В частности, для квадратной сетки hx=hy=h это условие принимает вид

r = chτ ≤ 12

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

Перейдём к построению разностных схем для волновых задач, сформулированных в виде системы дифференциальных уравнений. В качестве примера возьмём двумерное уравнение акустики

	∂s	+ A	∂s	+ B	∂s	= o ,	(11.71)
	∂t		∂x		∂y
			63

				0≤x≤lx, 0≤y≤ly , 0≤t≤T,
где
											0		0
			0 0 1/ ρ				0 0
u
		, A =	0 0		0		, B = 0 0 1/					ρ
s = v													, o = 0 .
			2							2
p			ρc 0		0			0	ρc		0			0

По аналогии с одномерным случаем, построение схемы Лакса-Вендроффа для уравнения (11.71) осуществляется на основе разложения

s(x,y,t + τ) = s(x,y,t) + τ	∂s	(x,y,t) +	1	τ2	∂2s	(x,y,t) + ... .
	∂t		2		∂t2

Используя уравнение (11.71) и полученное на его основе соотношение

∂2s		2 ∂2s		+ (AB + BA)	∂2s	2 ∂2s
∂t2	= A					+ B		,
			∂x2		∂x ∂y		∂y2

мы придем к двумерному аналогу выражения (11.5)

∂s

∂s

s(x,y,t + τ) = s(x,y,t) −τ A

(x,y,t) + B

(x,y,t)

+

∂x

∂y

.

2

2

2

1

2

2

∂ s

∂ s

2 ∂ s

τ

A

(x,y,t) + (AB + BA)

(x,y,t) + B

(x,y,t) .

2

2

∂x

∂x ∂y

∂y

2

Замена производных по x и y на соответствующие разностные отношения второго порядка точности даёт схему Лакса-Вендроффа для системы (11.71):

skn+,m1

τ22 A2

τ2

8hxhy

= s

k

k

+ B

k

k

+

−τ A

sn+1,m −sn−1,m

sn,m+1

−sn,m−1

k

n,m

2h

2h

x

y

k

k

k

k

k

k

sn+1,m −2sn,m

+ sn−1,m

2 sn,m+1 −2sn,m

+ sn,m−1

+ B

+

(11.72)

h

2

h

2

x

y

(AB + BA)(sk + + −sk − + −sk + − + sk − − ).

n 1,m 1 n 1,m 1 n 1,m 1 n 1,m 1

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

Из процедуры построения следует, что эта схема аппроксимирует уравнение (11.72) со вторым порядком по h и τ. Устойчивость схемы определяется условием (11.70).

Для вычисления граничных условий можно использовать метод фиктивных областей, вводя дополнительные (фиктивные) узлы с координатами (x-1,ym,tk), (xN+1,ym,tk), (xn,y-1,tk) и

64

(xn,yM+1,tk). Процедура вычисления решения в граничных точках аналогична рассмотренной нами процедуре для одномерных уравнений (пункт 11.3.3).

В пункте 11.3.5 мы рассмотрели одномерную схему Ф. Роу, которая основана на характеристической форме волнового уравнения. В многомерном случае такая форма имеет несколько более сложный вид. Рассмотрим систему волновых уравнений вида

∂s + ∑Ap	∂s = o ,
3
∂t p=1	∂xp

где (x1,x2,x3)=(x,y,z) и Ap – некоторые постоянные матрицы размером J×J. Умножим эту систему на lj, левый собственный вектор матрицы A1, которому соответствуют собственные значения λj. В результате получим систему уравнений

∂(lj s)			∂(lj s)		3	∂s
		+ λj			+ ∑lj Ap		= o , j=1, …, J,
∂t						∂xp
			∂x1 p=2

которая описывает распространение волн в направлении оси x1. Для волн, распространяющихся в других направлениях, получаются аналогичные системы. Построим характеристическую форму двумерных уравнений акустики (11.71), учитывая, что A1=A и A2=B. Собственные значения матрицы A определяются как λ=±c,0. Собственный вектор соответствующий λ3=0 не даёт ничего полезного, а два других вектора l1,2=(ρc,0,±1). Тогда для волн, распространяющихся в направлении оси ±x, получим следующие характеристические уравнения:

	∂	(p + cρ u) + c	∂	(p + cρ u) = −c2ρ	∂v	,	(11.73)
	∂t		∂x		∂y

∂∂t (p −cρ u)−c ∂∂x (p −cρ u) = −c2ρ ∂∂yv .

Аналогично, для волн, распространяющихся в направлении оси ±y, характеристические уравнения имеют вид

		∂		(p +cρ v) +c	∂		(p + cρ v) = −c2ρ	∂u	,	(11.74)
	∂t				∂y			∂x
		∂		(p −cρ v)−c	∂		(p −cρ v) = −c2ρ	∂u .
		∂t			∂y
								∂x

Для обобщения схемы (11.62) на случай двумерных уравнений требуется две сетки (рисунок 11.27). Для аппроксимации уравнений (11.73) введём сетку с узлами (xn,ym+1/2,tk), где

xn=n hx, n=0, …, N,

65

ym+1/2=(m+1/2) hy, m=0, …, M-1,

(11.75)

tk=kτ, k=0, 1, … ,

а для аппроксимации уравнений (11.74) - сетку с узлами (xn+1/2,ym,tk), где

xn+1/2=(n+1/2) hx, n=0, …, N-1,

ym=m hy, m=0, …, M,

(11.76)

tk=kτ, k=0, 1, … .

Такие сетки называются смещёнными. Они возникают, когда различные величины или группы величин вычисляются в различных узлах.

Рисунок 11.27. Пространственная сетка для аппроксимации уравнений (11.73) и (11.74); в узлах «•» определяются u и p; в узлах «×» определяются v и p.

Рассмотрим сначала построение разностной схемы для уравнений (11.73). Эти уравнения, при каждом значении y, представляют собой два уравнения переноса с отличной от нуля правой частью. Поэтому левая часть уравнений (11.73) аппроксимируется так же как в схеме (11.62) на сетке (11.75), а правая часть – на сетке (11.76). Тогда разностная схема примет следующий вид

66

Rk+1

−Rk

Rk

−Rk−1

Rk

−Rk

n,m+1/ 2

n,m+1/2 + n−1,m+1/ 2

n−1,m+1/2

+ c

n,m+1/2

n−1,m+1/ 2 =

2τ

2τ

hx

vk

−vk

−ρc2

n−1/2,m+1

n−1/2,m

,

hy

Sk+1

−Sk

Sk

−Sk−1

−c

Sk

−Sk

n,m+1/2

n,m+1/2 +

n+1,m+1/2

n+1,m+1/2

n+1,m+1/2

n,m+1/2 =

2τ

2τ

hx

,

vk

−vk

−ρc2

n+1/2,m+1

n+1/2,m

,

hy

Rk

= pk

+ cρ uk

,

n,m+1/2

n,m+1/ 2

n,m+1/2

Sk

= pk

−cρ uk

,

n,m+1/2

n,m+1/ 2

n,m+1/2

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

Выражая из этих соотношений неизвестные величины в момент времени tk+1, мы получим расчётные формулы. Окончательно, скорость и давление вычисляются как

uk+1+ = 1 (Rk+1+ −Sk+1+ ), n,m 1/2 2cρ n,m 1/2 n,m 1/2

pk+1+ = 1 (Rk+1+ + Sk+1+ ). n,m 1/2 2 n,m 1/2 n,m 1/2

Разностная схема для уравнений (11.74) строится аналогичным образом:

k+1

−Q

k

Q

k

−Q

k−1

Q

k

−Q

k

Q

+1/2,m +

+1/ 2,m

+1/2,m−1 =

n+1/ 2,m

n

n+1/ 2,m−1

n+1/2,m−1

+ c

n

n

2τ

2τ

hy

uk

−uk

−ρc2

n+1,m−1/2

n,m−1/2

,

hx

k+1

−T

k

T

k

−T

k−1

T

k

−T

k

T

+1/ 2,m +

n+1/2,m

n

n+1/2,m+1

n+1/2,m+1

−c

n+1/2,m+1

n+1/2,m =

2τ

2τ

hy

,

uk

−uk

−ρc2

n+1,m+1/2

n,m+1/2

,

hx

Qk

= pk

+ cρ vk

,

n+1/2,m

n+1/ 2,m

n+1/2,m

T k

= pk

−cρ vk

,

n+1/ 2,m

n+1/2,m

n+1/2,m

n=1, …, N-1; m=1,…, M-1; k=1, 2, … ,

и

67

vnk++1/1 2,m = 2c1ρ(Qnk++1/1 2,m −Tnk++1/1 2,m ),

pnk,+m1+1/2 = 21 (Qnk,+m1+1/ 2 +Tnk,m+1+1/2 ).

Из построения этой схемы следует, что она аппроксимирует двумерные уравнения акустики со вторым порядком по h и τ. Исследование устойчивости этой схемы приводит к очень громоздкому уравнению четвёртого порядка для собственных значений оператора перехода. Анализ этого уравнения даёт необходимое условие устойчивости рассматриваемой схемы:

	1		1
cτ		+			≤1 ,
			h
h
	x		y

а в случае квадратной сетки hx=hy=h -

chτ ≤ 21 .

Для построения разностных схем для многомерных задач можно использовать метод расщепления. Рассмотрим этот подход на примере схемы С. К. Годунова. Следуя процедуре представленной в пункте 11.2.3, вычисление решения в момент времени tk+1 по известному решению в момент времени tk производится в два этапа. На первом этапе вычислим промежуточные значения, решая одномерные задачи по x с использованием схемы (11.66):

	k+1/2			k		k	k
	sn+1/2,m+1/2			−sn+1/2,m+1/2	+ A	sn+3/2,m+1/2	−sn−1/ 2,m+1/ 2	−
				τ		2hx			(11.77)
	1		A2 (snk +3/ 2,m+1/ 2 −2snk +1/2,m+1/2 + snk −1/2,m+1/2 ) = 0,
		2chx

n=1, …, N-2; m=0,…, M-1 .

Расчёт значений в ячейках, примыкающих к границам, которые имеют индексы (1/2,m+1/2) и (N-1/2,m+1/2), производятся с учётом граничных условий. На втором этапе решаются одномерные задачи по y, используя аналогичную схему, основанную на уже вычисленном промежуточном решении:

k+1				k+1/2			k+1/2	k+1/2
sn+1/2,m+1/2				−sn+1/2,m+1/ 2		+ B	sn+1/2,m+3/2 −sn+1/2,m−1/2		−
				τ				2hy		(11.78)
1
		2		k+1/2		k+1/2		k+1/2
		B	(sn+1/2,m+3/2		−2sn+1/2,m+1/2			+ sn+1/2,m−1/ 2 ) = 0,
	2chy

n=0, …, N-1; m=1,…, M-2,

68

+расчётные формулы для ячеек, примыкающих к границам.

Врезультате мы получим требуемое решение в момент времени tk+1.

Исследуем устойчивость этой схемы. Для этого зададим решение в виде

		k
	u
k
sn,m	= v
	p
		n+1/2,m+1/2

u
= v		exp(iα n + iα m).
			1	2
	p

Подстановка этого решения в схему (11.77) даёт следующее выражение для промежуточного решения

			k+1/2	λ1(α1)u
		u
	k+1/ 2
s					v
	n,m	= v		=				exp(iα n + iα m),
									1	2
		p
				λ1(α1)p
			n+1/ 2,m+1/ 2

где λ1(α1) - собственные значения оператора перехода для одномерной задачи (смотри пункт 11.2.3). Далее, подставляя полученное промежуточное решение в схему (11.78) мы получим

k+1 = u k+1 sn,m v

p n+1/2,m+1/2

	λ	(α )u
	1	1
	λ	(α )v
=				exp(iαn + iα m).
	2	2			1	2
λ1(α1)λ2(α2 )p

Тогда для устойчивости необходимо выполнение условий

|λ1(α1)|≤1,

|λ2(α2)|≤1,

|λ1(α1) λ1(α1)|=|λ1(α1)| |λ2(α2)|≤1,

для всех 0≤α1, α2≤2π. Таким образом, устойчивость двумерной схемы определяется устойчивостью одномерных схем. В свою очередь, одномерные схемы будут устойчивы если

cτ	≤1 и cτ	≤1,
h	h
x	y

и, следовательно, двумерная схема будет устойчива при условии

cτ	≤1 .
min(hx ,hy )

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

69

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

snk+,m1 = Dx (snk −1,m , snk ,m , snk +1,m , τ,hx )

(11.79)

представляет собой некоторую схему для решения системы

∂s	+ A		∂s		= o ,
∂t			∂x
а
snk+,m1 = Dy (snk ,m−1 , snk ,m , snk ,m+1 , τ,hy )						(11.80)
есть схема для решения системы
∂s	+ B		∂s		= o .
∂t			∂y

В этих обозначениях схема расщепления (11.77), (11.78) выглядит следующим образом:

sk+1/2 = D (sk − , sk , sk + , τ,h ),

n,m x n 1,m n,m n 1,m x

n=1, …, N-1; m=0,…, M,

snk+,m1 = Dy (snk+,m1/−12 , snk+,m1/ 2 , snk+,m1/+21 , τ,hy ),

(11.81)

n=0, …, N; m=1,…, M-1,

Построение более точной схемы расщепления требует, в общем случае, большего объёма вычислений. Наиболее эффективной представляется следующая процедура:

q(1)n,m = Dx (snk −1,m , snk	,m , snk +1,m ,	1	τ,hx ),
		2
n=1, …, N-1; m=0,…, M,
q(2)n,m = Dy (q(1)n,m−1 , q(1)n,m , q(1)n,m+1 , τ,hy ),				(11.82)

n=0, …, N; m=1,…, M-1,

skn+,m1 = Dx (q(2)n−1,m , q(2)n,m , q(2)n+1,m , 21 τ,hx ),

n=1, …, N-1; m=0,…, M.

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

70