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

Количество шагов по времени k, необходимое для того, чтобы ошибка начального приближения уменьшилась в 10p раз, можно оценить из соотношения α=10-p. Тогда

					p ln(10)
k = −
					−2 sin2 (π 2N )						)
			ln 1
				(
и для больших значений N
	2					2		p				2
k ≈			p ln(10) N					≈		N			.
	π	2							2

Для расчёта одного шага по времени по схеме (11.110) требуется O(N2) операций. Тогда общее число операций будет O(pN4/2)=O((Ns)2).

Можно построить более эффективный метод, используя схему переменных направлений (11.49). Для уравнения (11.109) она записывается следующим образом:

	wn,m −vnk,m		a		k	1
		=		( xwn,m +	yvn,m )−		f (xn ,ym ) ,	(11.117)
	τ		2			2

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

w0,m = g1(ym ) , wN ,m = g3(ym ) , m=0, …, M,

	vnk,+m1 −wn,m		a		k+1	1
		=		( xwn,m +	yvn,m )−		f (xn ,ym ) ,
	τ		2			2

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

vn0,m = g(xn ,ym ),

vnk,0+1 = g2(xn ) , vnk,+M1 = g4(xn ) , n=0, …, N.

Поведение ошибки (11.112) для этой схемы также описывается выражением (11.115), только

λ	=	(1 −β1	sin2	(πr	2N ))(1 −β2 sin2 (πs 2N )) .
r,s		(1 + β1	sin2	(πr	2N ))(1 + β2 sin2 (πs 2N ))

Минимальное значение α достигается при

τ =	hxhy		≈	Nhxhy	(11.118)
	2a sin(π2N )cos(π2N )			aπ

и для больших значений N равно

91

					h			γ
	π hx				y
α = 1 −				+			≈ 1 −		.
								N
	N h				h
			y		x

Количество шагов по времени k, необходимое для того, чтобы ошибка начального приближения уменьшилась в 10p раз, оценивается как

p ln(10)
k ≈		N .
	γ

Тогда общее число операций будет O((Ns)3/2), что примерно соответствует вычислительным затратам итерационного метода SOR.

Дальнейшее ускорение сходимости можно достигнуть, проводя вычисления со специальным образом подобранным набором переменных шагов по времени {τk}. Однако построение такого набора возможно только для задач допускающих разделение переменных, а в этом случае методы на основе FFT являются более эффективными. Поэтому метод установления обычно применяется для решения задач с переменными коэффициентами в областях сложной формы, а также для некоторых нелинейных задач. Естественно, что для таких задач невозможно получить аналитические выражения для оптимального шага по времени, но всё же полученные нами соотношения можно использовать как руководство для выбора τ.

В качестве примера рассмотрим построение разностной схемы для уравнения (11.85) при a(u)>0 и граничными условиями Дирихле. По аналогии с (11.98) схему переменных направлений в этом случае можно записать в следующем виде:

wn,m −vnk,m

=

ank+1/2,m

(wn+1,m −wn,m )−

ank−1/2,m

(wn,m −wn−1,m )+

τ

k

2h2

2h2

x

x

,

ank,m+1/2

ank,m−1/2

k

k

k

k

1

(vn,m+1

−vn,m )−

(vn,m −vn,m−1 )−

2 f (xn ,ym )

2h2

2h2

y

y

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

vnk,+m1 −wn,m

=

ank+1/ 2,m

(wn+1,m −wn,m )−

ank−1/ 2,m

(wn,m −wn−1,m )+

τ

k

2h2

2h2

x

x

,

ank,m+1/2

ank,m−1/ 2

k+1

k+1

k+1

k+1

1

(vn,m+1

−vn,m )−

(vn,m −vn,m−1 )−

2 f (xn ,ym )

2h2

2h2

y

y

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

vn0,m = g(xn ,ym ),

92

v0,k+m1 = g1(ym ) , vNk+,m1 = g3(ym ) ,

m=0, …, M,

vnk,0+1 = g2(xn ) , vnk,+M1 = g4(xn ) ,

n=0, …, N.

где

ak+ = a (1 (vk+

n 1/ 2,m 2 n 1,m

ak + = a (1 (vk +

n,m 1/2 2 n,m 1

+ vnk,m )), . + vnk,m )).

Здесь для приближения нелинейных членов мы использовали явную аппроксимацию, поэтому вычисления по этой схеме производятся также как и для линейной схемы (11.117). Для выбора шага по времени используем соотношение (11.118), учитывая, что теперь параметр a зависит от решения задачи:

	τk =			(N + M )hxhy				.
		(	n,m	( n,m )	n,m	( n,m ))
		π	maxa vk		+ mina vk

________________________________________________________________________________

Пример 11.10 (метод установления)

Рассмотрим уравнение (11.2) при a(u)=exp(-u) в квадрате {0≤x≤1,0≤y≤1} с

граничными условиями
u(0,y)=0, u(1,y)=ln(1+y),	(11.119)

u(x,0)=0, u(x,1)=ln(1+x).

Правая часть уравнения

2(x2 +y2 ) f (x,y) = − (1+xy)3 .

Зададим разностную сетку с N=M=10. При начальном условии

vn(0),m = 0 ,

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

и ε=10-5 процесс установления сходится за 14 шагов. Приближённое решение показано на рисунке 11.29. Точное решение этой задачи – u(x,y)=ln(1+xy). Ошибка приближённого решения

93

e

=

v(14) −u(e)

2

≈ 3 10−4 .

a

u(e)

2

Рисунок 11.29. Приближённое решение уравнения (11.85) с граничными условиями (11.119),

метод установления (N=M=10, ε=10-5).

________________________________________________________________________________

94

ГЛАВА 6

Движение сжимаемой жидкости (газа)

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

Трудности аналитического исследования определяются сложной математической постановкой: в общем случае – это краевые задачи для нестационарных нелинейных уравнений, описывающие вихревые течения сжимаемой жидкости (газа). Течение жидкости (газа) часто характеризуется внутренними разрывами, такими как, например, ударные волны, на которых требуется ставить специальные условия. Эти условия описываются известными уравнениями Ренкина-Гюгонио, но их применение на практике затруднено тем, что поверхности разрывов перемещаются по жидкости, и это движение не известно заранее, а определяется из самих дифференциальных уравнений и условий на разрыве. Поэтому, как правило, лишь численное моделирование и довольно сложные физические эксперименты позволяют определить необходимые характеристики течения.

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

6.1. Схема Лакса-Вендроффа для одномерных задач.

Сначала запишем уравнения гидродинамики (смотри параграф 2.5) в дивергентной форме, то есть в виде законов сохранения. Введём новые зависимые переменные ρu и

ρ(e+1/2u2). Тогда уравнения гидродинамики в одномерном случае примут вид

95

∂s	+	∂ f (s)	= o ,	(6.1)
∂t		∂ x

где s и f(s) - векторы определённые равенствами

s1				ρ
s = s		=		ρu			,
2					2
						)
s3	ρ(e + 12 u

(6.2)

f1 f (s) = f2f3

ρu

s2

ρu

2

+ p

=

2

/ s1 + p

=

s2

.

u (ρe + 1

2 ρu

2

+ p )

(s3 + p )/ s1

s2

Мы будем рассматривать движение идеального газа с уравнением состояния

e = ρ(γp−1) .

Тогда давление и температура определяются как

p = ρe(γ −1) = (γ −1)	s		−	s22		,	T =	p	.
				2s
		3						ρR
				1

где R - газовая постоянная. Поэтому f(s) действительно является функцией от s. Система (6.1) называется дивергентной формой уравнений гидродинамики. Физическая интерпретация этого определения состоит в следующем: если проинтегрировать уравнения (6.1) по некоторому объёму, то мы получим законы сохранения величины s для этого объёма.

Введём разностную сетку с узлами (xn, tk), где xn=nh, n=0, 1, … и tk+1=tk+τk. Для аппроксимации уравнений (6.1) удобнее использовать двухшаговую схему ЛаксаВендроффа, аналогичную схеме (5.7) и (5.8). Сначала в центрах прямоугольных ячеек на плоскости (x,t) (рисунок 5.3) вычисляются промежуточные значения:

snk++1/1/22 = 12 (snk+1 + snk )− 2τkh ( fnk+1 − fnk ),

n=1, …, N-1,

где

fnk = f (snk ).

Окончательное решение вычисляется затем по схеме

snk +1 = snk −	τk ( fnk++1/1/22	− fnk−+1/1/22 ).	(6.3)
	h

Так же как и в случае уравнений акустики эта схема имеет второй порядок аппроксимации по h и τ.

96

Непосредственное исследование устойчивости схемы (6.3) невозможно, из-за нелинейности системы (6.1) которую эта схема аппроксимирует. Поэтому применим принцип замороженных коэффициентов. Сначала линеаризуем систему (6.1). Для этого представим решение в виде s=s0+w, где s0 – постоянный вектор вида

	ρ0
	ρ0u0		,
s0 =
ρ0	(e0 + 12 u02	)

а вектор w – малое возмущение. Тогда f(s)=f(s0+w)≈f(s0)+A(s0)w+…,

где

A( s0 ) = {aij ( s0 )} = ∂∂sfi ( s0 ) .

j

Подставляя это разложение в (6.1), мы получим линейную систему уравнений с постоянной матрицей A0=A(s0) вида

∂w	+ A ∂w	= o .	(6.4)
∂t	0 ∂x

Тогда схема (6.3) переходит в схему (5.6) с матрицей A=A0, а для анализа этой схемы уже можно применить критерий фон Неймана. Следуя процедуре, рассмотренной в параграфе 5.3, можно получить следующее выражение для собственных значений оператора перехода

λj (α) =1 − i ωhj τk sin(α) − ωhj τk 2 (1 − cos(α) ) ,

j=1, 2, 3,

где ωj есть собственные значения матрицы A0, которые равны u0, u0±с0, здесь с0=(γp0/ρ0)1/2 – скорость звука. Из критерия фон Неймана следует, что схема для линеаризованного уравнения (6.4) будет устойчива если

τ	k	max | ω	j	|		τk max (| u0 |,| u0 + c0 |,| u0 − c0 | )
		j			=		≤1.
		h				h

Тогда шаг по времени выбирается из условия

τk ≤ | u0 h| +c0 .

Это условие можно использовать для вывода условия устойчивости схемы (6.3). Учитывая, что при течении газа скорость среды и скорость звука являются функциями координаты и времени, мы приходим к условию устойчивости схемы (6.3):

97

τk=rτmax,

h

1/ 2

τk ≤

= τmax ,

cnk

= (γpnk / ρnk )

.

max (| unk | +cnk )

n

или в терминах компонент вектора s

τk

≤

h

.

(s

)k

γ

(s

)k

1

(s

)k 2 1/ 2

max

2

n

+

(γ −1)

3

n −

2

2

n

n

(s

)k

(s

)k

(s

)k

1

n

1

n

1

n

Величина

max (| unk | +cnk ),

(6.5)

n

оценивает наибольшую скорость волны в заданной области в момент времени t=tk. Данная оценка получена из анализа устойчивости линеаризованной системы (6.4), поэтому в некоторых случаях величина (6.5) будет давать заниженное значение наибольшей скорости волны. Тогда вычисления с шагом по времени τk близким к τmax могут быть неустойчивыми. Поэтому шаг по времени задаётся как

(6.6)

где r – параметр Куранта (0<r<1) и на практике часто выбирается r=0.9. Для более надёжного достижения устойчивости можно использовать следующий подход. Расчёт первых нескольких временных шагов производить с малым параметром Куранта. Это позволит течению сформироваться и оценка (6.5) станет более надёжной. Тогда последующие вычисления можно производить с большим параметром Куранта.

Для расчёта граничных условий вводятся фиктивные слои, что позволяет не нарушать единообразия вычислений для граничных ячеек. Мы рассмотрим два вида граничных условий: твёрдая стенка (или ось симметрии) и открытая граница расчётной области. В качестве примера рассмотрим расчёт этих условий на левой границе расчётной области. Пусть граничным точкам отвечают индексы (0, k), а соответствующим фиктивным узлам – индексы (-1, k). Тогда граничное условие на твёрдой стенке

u(0,t)=0, t>0

запишется как

0.5(u−k1 + u1k )= u0k			= 0 ,		ρ−k1 = ρ1k , p−k1 = p1k
или
s1	k				s1	k
			=
s2					−s2 .
s					s
3		−1			3
						1
			98

Далее расчитываются величины

s−k1/+1/22 и s0k +1

по приведённым выше формулам. Через открытые границы газ (жидкость) может втекать или вытекать из области, и здесь должны быть обеспечены условия непрерывности движения. Пусть, например, газ (жидкость) втекает в область с левой стороны. Тогда в фиктивных узлах задаются параметры набегающего потока (помечены знаком «∞»):

ρ k

u

p −1

ρ

=u .p ∞

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

s1 k

s2

s3 −1

s1 k

=s2 .s3 0

________________________________________________________________________________

Пример 6.1 (схема Лакса-Вендроффа)

В качестве тестового примера рассмотрим задачу об ударной трубе (задача Сода).

Пусть бесконечная труба, заполненная воздухом (γ=1.4), в точке x=0 разделена перегородкой. При t=0 состояние газа описывается кусочно-постоянной функцией:

ρ			(ρL		0 pL ), x ≤ 0
						.
u			=	(ρR	0 pR ), x > 0
p		t=0

В момент времени t=0 перегородка мгновенно убирается и эта конфигурация распадается на две волны распространяющиеся влево и вправо от точки x=0. Например, если

	ρL			1
		=		0		и
uL
				1
	pL

	ρR			0.125
		=		0		,	(6.7)
uR
				0.1
	pR

то влево распространяется волна разрежения, а вправо – контактный разрыв и ударная волна. На рисунках 6.1-6.4 показано приближённое решение уравнений (6.1) с начальными условиями (6.7), полученное с помощью схемы Лакса-Вендроффа, а также точное решение этой задачи.

99

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

Рисунок 6.1. Пространственное распределение плотности в момент времени t=0.5037 ; -

точное решение, - приближённое решение (схема Лакса-Вендроффа, h=0.02).

________________________________________________________________________________

6.2. Метод искусственной вязкости.

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

∂s	+	∂ f (s)	=	∂		μB(s)	∂s	,	(6.8)
∂t		∂ x
				∂ x			∂x

где B(s) – квадратная матрица; μ - коэффициент искусственной вязкости. Матрица B(s) должна быть подобрана таким образом, чтобы решение s(x,t) системы (6.8) обладало достаточной гладкостью и при μ→0 приближалось к решению исходной системы (6.1). Аппроксимируя систему (6.8) соответствующей разностной схемой, мы получим разностные уравнения с диссипативными членами, которые и определяют вязкостные эффекты схемы. Следует отметить, что введение искусственной вязкости имеет смысл только для схем

100