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

которой решение задачи о распаде разрыва не будет иметь физического смысла. Это приводит к остановке вычислительного процесса. Поэтому подавление нефизических осцилляций в методе WAF имеет первостепенное значение. В простейшем случае можно применить метод искусственной вязкости как в схеме (6.9), но эта процедура требует настройки алгоритма для каждой конкретной задачи. Существует более универсальный способ уменьшения амплитуды осцилляций, основанный на применении специальных параметров (ограничителей), зависящих от изменения приближённого решения.

111

ГЛАВА 8

Движение несжимаемой вязкой жидкости

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

8.1 Введение

Уравнения Навье-Стокса для несжимаемой жидкости имеют вид: 1) уравнение неразрывности

		div(u) = 0 ,
2) уравнение движения		(8.1)
ρ	∂u	+ ρ(u )u = − p + νΔu + fe (x,t) ,
	∂t
		x D, t≥0,

u(x,0) = g0 (x) + граничные условия.

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

					3	∂∂uxk = 0 ,
				∑
				k =1			k
	∂u	3	∂u		∂ p		3	∂2u
ρ	∂tn + ρ∑uk ∂ xn			= −			+ ν∑		2n + fe,n (x1	, x2	, x3 ,t) ,
					∂ x			∂x
		k =1	k			n	k =1		k
						112

n=1, 2, 3.

Здесь

	u			u1
u =			=			-
		v		u2
	w			u
				3

скорость частиц жидкости в точке x=(x1, x2, x3) в момент времени t, ρ - плотность жидкости, p=p(x, t) - давление, ν - коэффициент вязкости и fe(x, t) - внешние силы.

Запишем уравнения (8.1) в безразмерном виде. Для этого введём новые переменные u ,

p , x и t :

u = u uc , p = p ρuc2 , x = x l ,

t = t l / uc , fe = fe ρuc2 / l, ,

где l и uc представляют некоторые характерные параметры задачи. После этого уравнения (8.1) принимают вид

div(u ) = 0 ,

(8.2)

	∂u	+ (u )u = −p +	1	u + fe (x ,t ) ,
	∂t		Re

здесь Re=ρucl/ν - число Рейнольдса. В дальнейшем, мы будем использовать безразмерные уравнения. При этом знак « » мы будем опускать. Характер течения определяется величиной числа Рейнольдса. Можно выделить три основных режима течения:

1)Re>>1. В этом случае вязкие силы значительно меньше инерционных сил, что приводит к очень сложной картине течения из-за возникновения турбулентности.

2)Re=O(1). В этом режиме вязкие силы сравнимы с инерциальными.

3)Re<<1. Этот режим характерен для течения очень вязких жидкостей. В этом случае уравнения (8.2) можно заменить линейным стационарным уравнением.

Мы будем рассматривать методы решения задач для которых реализуются первый и второй режимы течения.

Одной из трудностей в решении уравнений (8.2) является то, что отсутствует эволюционное уравнение для давления. Для несжимаемой жидкости, давление не имеет своего обычного термодинамического смысла, а является следствием несжимаемости.

Обычно мы сталкиваемся с задачами о течении жидкости в некоторых замкнутых областях или с задачами об обтекании тел. В этом случае возникает вопрос о постановке

условий на границах области или поверхности тела. Из вида уравнений (8.2) следует, что

113

удобнее всего формулировать граничные условия для скорости. Разберём, например, постановку граничных условий в двумерном случае. Пусть n есть единичный вектор нормали к границе в точке (x, y), а k - вектор касательной в той же точке. Часто на практике задаются следующие условия:

1) для нормальной скорости -

(u,n ) = gn (x, y,t) ,

(8.3)

здесь ( , ) – скалярное произведение, 2) для касательной скорости -

(u,k ) = gk (x, y,t) ,

(8.4)

или условие свободного проскальзывания

∂( )

∂n = ( (u,k ),n ) = 0 , (8.5)

3)равенство нулю нормального и касательного напряжения на свободной поверхности (более подробно этот случай мы рассмотрим позже).

Мы рассмотрим два подхода к решению уравнений Навье-Стокса. Первый из них

основан на уравнениях (8.2), где конвективный член записывается в дивергентной форме, что даёт следующие уравнения u,k

div(u) = 0 ,

(8.6)

∂u

+

∂ f

x

(u)

+

∂ fy (u)

+

∂ f

z

(u)

= −p +

1

u + fe (x,t) ,

∂t

∂ y

Re

∂ x

∂ z

Здесь

	u2
fx (u ) =	u v ,

u w

v u		w u
fy (u ) = v2	, fz (u ) =	w v .
			2
v w		w

Если течение не содержит свободных поверхностей, то можно исключить давление, взяв операцию rot от уравнения движения (8.2). Учитывая, что rot(grad(p))=0 можно получить

	∂ω	+ (u )ω − (ω )u =	1	ω + rot ( fe (x,t) ),	(8.7)
	∂t		Re

ω(x,0) = rot ( g0 (x) ), x D, t≥0,

где ω=rot(u) – завихрённость. Левая часть уравнения описывает перенос завихрённости вместе с жидкостью. Определим скорость через функцию тока ψ следующим образом:

114

u = rot (ψ) .

(8.8)

Так как div(rot(ψ))=0, то условие неразрывности выполняется. Для того чтобы замкнуть наши уравнения нам необходимо соотношение между завихрённостью и функцией тока. Для этого возьмём операцию rot от (8.8) и получим

rot (u ) = ω = rot (rot (ψ)) = − ψ + grad(div(ψ)).

Функция тока определяется из уравнения (8.8) неединственным образом. Поэтому нам необходимо некоторое дополнительное условие. Удобнее всего положить div(ψ)=0. Тогда связь между ω и ψ выражается уравнением Пуассона

ψ = −ω.

(8.9)

Таким образом, мы получили замкнутый набор уравнений для описания движения несжимаемой жидкости: эволюционное уравнение для завихрённости, уравнение Пуассона (8.9) для функции тока из которой определяется поле скоростей согласно уравнению (8.8).

Граничные условия для ω и ψ следуют из граничных условий для u. Уравнения (8.7)-(8.9) называются уравнениями Навье-Стокса в переменных функция тока-завихрённость. Эта форма особенно удобна в случае двумерных задач, когда

	u(x, y,t)							fe,1	(x, y,t)
u =				и	fe	=		fe,2
	v(x, y,t)								(x, y,t) .
		0							0

Тогда нетрудно получить, что

0

0

ω = 0

и ψ = 0

.

ω

ψ

Система уравнений (8.7)-(8.9) при этом значительно упрощается и принимает вид

∂ω

+ u

∂ω

+ v

∂ω

=

1

∂2ω

+

∂2ω

+ fs (x, y,t) , (8.10)

∂t

∂x

∂ y

∂ x2

Re

∂ x2

ω(x, y,0)

∂g0,2

∂g0,1

, (x, y) D, t≥0,

=

−

∂x

∂ y

∂2ψ

+

∂2ψ

= −ω(x, y,t) ,

(8.11)

∂ x2

∂ y2

u

=

∂ψ ∂ y

(8.12)

,

v

−∂ψ ∂ x

где

fs (x, y,t) = ∂∂f2,ye − ∂∂f1,xe .

115

Обычно граничные условия задаются на u и v. Тогда из уравнения (8.8) мы получим граничные условия для ψ, которые используются при решении уравнения (8.9).

8.2 Разностные схемы в переменных функция тока-завихрённость.

Перейдём к построению разностных схем для уравнений (8.10) - (8.12) в прямоугольной области {0 ≤ x ≤ lx , 0 ≤ y ≤ ly}. Предположим, что fe=const. Граничные условия примем в следующем виде:

u									0		,
v						= g			( y,t)
					(0, y )		1
u						g		2	(x,t)		,
						=
v					( x,0 )				0

которые, согласно (8.12), преобразуются в

ψ

0

,

∂ψ ∂ x

=

−g

( y,t)

ψ

(0, y )

1

0

=

,

∂ψ ∂ y

( x,0 )

g2

(x,t)

u

0

,

=

v

(lx , y )

g3

( y,t)

u

g

4 (x,t)

,

=

v

(x,ly )

0

ψ

=

0

,

∂ψ ∂ x

(lx , y )

−g3 ( y,t)

ψ

=

0

.

∂ψ ∂ y

(x,ly )

g4 (x,t)

Введём разностную сетку с узлами (xn, ym, tk) и шагами hx=lx/N, hy=ly/M и τk. Для построения разностной схемы для уравнения (8.10) применим метод расщепления. Производная по времени аппроксимируется разностью вперёд, первые производные по пространству – центральной разностью и вторые производные – обычным трёхточечным разностным отношением. В результате получим

ωk +1/ 2

− ωk

ωk

− ωk

ωk

− 2ωk

+ ωk

n,m

n,m

= −unk,m

n+1,m

n−1,m

+

n+1,m

n,m

n−1,m

,

τk

2hx

Re h2

x

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

(8.13)

ωk +1

− ωk +1/ 2

ωk +1/ 2

− ωk +1/ 2

ωk +1/ 2

− 2ωk +1/ 2 + ωk +1/ 2

n,m

n,m

= −vnk,m

n,m+1

n,m−1

+

n,m+1

n,m

n,m−1

,

τk

2hy

Re hy2

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

Уравнение (8.11) представляет собой уравнение Пуассона с постоянными коэффициентами. Поэтому для аппроксимации этого уравнения можно использовать схему (7.17) с a=1 и b=0,

116

левую часть которой мы обозначим как P(un, m). Тогда разностная схема для определения ψ в момент времени tk+1 имеет вид

P(ψnk,+m1 )= −ωnk,+m1 ,

(8.14)

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

Аппроксимируя уравнения (8.12), мы получим формулы для вычисления скоростей в момент времени tk+1:

		ψk +1		− ψk +1
unk,+m1	=		n,m+1	n,m−1	,
			2hy
				(8.15)
			ψk +1	− ψk +1
vnk,+m1	= −		n+1,m	n−1,m		,
			2hx
n=1, … , N-1 ; m=1, … , M-1.

Нам осталось сформулировать граничные условия для уравнений (8.10) и (8.11). Для примера рассмотрим граничные условия при x=0. Первое условие в разностной форме имеет простой вид:

ψ0,k +m1 = 0 , m=0, … , M.

(8.16)

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

ψk +1 − ψk +1

1,m 2hx −1,m = −g1( ym ,tk +1) .

Как обычно, запишем уравнения (8.14) при n=0:

1	(ψ1,k +m1 − 2ψ0,k +m1 + ψ−k1,+1m )+	1	(ψ0,k +m1+1 − 2ψ0,k +m1 + ψ0,k +m1−1 )= −ω0,k +m1 .
h2		h2
x		y

Далее, выражая фиктивное значение из предыдущего соотношения и учитывая условие (8.16), окончательно получим

ω0,k +m1 = −	2	(ψ1,k +m1 + hx g1( ym ,tk +1) ),	(8.17)
	h2
	x

m=1, … , M-1.

Это соотношение определяет граничное условие для завихрённости при x=0, так как ψ1, m в момент времени tk+1 нам известно из решения уравнения (8.14). Вычисление граничных условий на других границах осуществляется совершенно аналогичным способом.

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

117

«замороженных» коэффициентов, мы должны сначала исследовать на устойчивость схему (8.13) для уравнения

∂ω

= −c

∂ω

+

1 ∂2ω

,

∂t

∂ x

Re

∂x2

где c=const. Из критерия фон Неймана можно получить ограничение на шаг по времени

2Re h2

1

τ ≤ min

x

,

Re h2

.

4 + c2 Re2 hx2

2

x

Распространение этого результата на схему (8.13) приводит с следующему условию

τk ≤ min (τ1,τ2 ),

(8.18)

где

2 Re h2

1

τ1

≤ min

x

,

Re hx2

,

+

(

max (uk

)2

)

Re2 h2

2

4

n,m

n,m

x

2Re hy2

1

2

τ2

≤ min

,

Re hy

.

+

(

max (vk

)2

)

Re2 h2

2

4

n,m

n,m

y

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

(8.17), что приводит к крайне неэффективному алгоритму для нахождения (ωk+1,ψk+1). Решение этой проблемы требует изменения процедуры интегрирования по времени. Схему расщепления можно представить в следующем виде

ds0,m

= 0 ,

dsN ,m

= 0 ,

dt

dt

dsn,m

= −unk,m

sn+1,m − sn−1,m

+

sn+1,m − 2sn,m + sn−1,m

,

dt

2hx

Re h2

x

tk ≤ t ≤ tk+τk ; n=1, … , N-1 ; m=1, … , M-1,

sn,m (tk ) = ωnk,m ,

dqn,0

= 0 ,

dqn,M

= 0 ,

dt

dt

dqn,m

= −vnk,m

qn,m+1 − qn,m−1

+

qn,m+1 − 2qn,m

+ sn,m−1

,

dt

2hy

Re hy2

tk ≤ t ≤ tk+τk ; n=1, … , N-1 ; m=1, … , M-1,

qn,m (tk ) = sn,m (tk + τk ).

118

(8.19)

(8.20)

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

ωk +1 = q

(t

k

+ τ

k

).

n,m n,m

Схемы (8.19), (8.20) можно записать в матричном виде

dsm

= Am sm , tk ≤ t ≤ tk+τk,

(8.21)

dt

sm(tk) – задано,

m=1, … , M-1,

dqn

= B q

n

, t

k

≤ t ≤ t

+τ ,

(8.22)

dt

n

k

k

qn(tk) – задано,

n=1, … , N-1,

где

sm = {s0,m , ... ,

sN ,m }T ,

qn = {qn,0 , ... , qn,M }T .

и Am, Bn – трёхдиагональные матрицы. Устойчивость того или иного метода решения систем обыкновенных дифференциальных уравнений вида ds/dt=As зависит от собственных значений матрицы A, а именно: параметр p=τkλj(A) должен принадлежать некоторой области на комплексной плоскости для всех j=1 , … , J, где J - число уравнений. Форма области зависит от метода. Оценим собственные значения матриц Am и Bn. Эти матрицы не являются постоянными, поэтому на основе принципа «замороженных» коэффициентов мы будем работать с некоторой постоянной матрицей A0, которая имеет такую же структуру как матрицы Am и Bn. Нас интересует случай больших чисел Рейнольдса. Поэтому в пределе

Re→∞ матрица A0 будет иметь вид

				0				0					0				0
		1		c	0 −c						c		1	0 −1
A0																			,
	=				.	.	.			=				.	.	.
		2h				c	0				2h				1 0		−1
								−c
				0				0					0				0

где h принимает значение либо hx, либо hy и c=const. Собственные значения этой матрицы равны

λj ( A0 ) = −i hc sin (πjL ), j=1 , … , L,

119

то есть все собственные значения принадлежат отрезку [0, -i c/h], здесь i – мнимая единица.

Возвращаясь к матрицам Am и Bn, можно заключить, что при Re→∞ собственные значения матрицы Am в момент времени tk принадлежат отрезку

i

k

i

k

λj ( Am )

−

max

un,m

,

max

un,m

,

(8.23)

hx n,m

hx n,m

j=1 , … , N

и аналогично для матрицы Bn

i

k

i

k

λj ( Bn )

−

max

vn,m

,

max

vn,m

,

(8.24)

hy n,m

hy n,m

j=1 , … , M.

Эти результаты говорят о том, что для интегрирования систем (8.21) и (8.22) необходимо применять метод, область устойчивости которого включает мнимую ось. Наиболее подходящими являются методы Рунге-Кутта порядка выше второго (смотри ОНВ, рисунок 6.4). Возьмём трёхстадийный метод третьего порядка. На рисунке 8.1 показана часть области устойчивости этого метода вблизи мнимой оси. Схема Рунге-Кутта 3-го порядка будет устойчива для систем (8.21) и (8.22), если параметр p для собственных значений (8.23) и (8.24) будет принадлежать отрезку [-i a, i a ] , где a ≈1.6. Учитывая, что при конечном значении числа Рейнольдса члены со второй пространственной производной также влияют на устойчивость, мы приходим к условию

	τk ≤ min (τ1,τ2 ),
где
τ1 =	1.6 min(hx ,hy )									,
	max (maxn,m		unk,m		,maxn,m		vnk,m		)

τ2 = 12 Re(min(hx ,hy ) )2 .

которое накладывает более слабое ограничение на шаг по времени, чем условие (8.18), при том же шаге по пространству.

120