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

Рисунок 8.1. Область устойчивости метода Рунге-Кутта 3-го порядка вблизи мнимой оси. Область устойчивости находится слева от кривой.

В качестве примера можно привести схему Рунге-Кутта 3-го порядка, которая задаётся следующей таблицей Бучера:

0
1/2	1/2
3/4	0	3/4
	2/9	1/3	4/9

Тогда процедура решения системы (8.21) строится следующим образом:

sk +1	− sk		1
m	m	=		(2d1	+ 3d2	+ 4d3 ),
			9
τk

где

d1 = Am smk , d2 = Am (smk + 12 τk d1 ), d3 = Am (smk + 34 τk d2 ).

Аналогичная процедура применяется затем к системе (8.22).

Суммируя всё изложенное выше, процедура вычисления одного шага по времени организуется следующим образом:

1)по известному в момент времени tk полю скоростей вычисляется шаг по времени τk,

2)во всех внутренних точках области вычисляется завихрённость в момент времени tk+1. Для этого используется либо схема (8.13), либо схема (8.21), (8.22),

3)функция тока в момент времени tk+1 определяется из решения уравнения (8.14) с граничными условиямми вида (8.16). Так как коэффициенты этого уравнения постоянны, то можно применить методы на основе FFT,

121

4) по известной функции тока вычисляется завихрённость в граничных узлах и поле скоростей в момент времени tk+1 по формулам (8.15).

______________________________________________________________________________

Пример 8.1 (движение жидкости в полости)

Рассмотрим движение жидкости в квадратной полости {0≤x≤1, 0≤y≤1}. Пусть в момент времени t=0 жидкость покоится:

u(x, y, 0) = v(x, y, 0)=0.

Нормальные к поверхности полости компоненты скорости равны нулю: u(0, y, t) = u(1, y, t)=0, v(x, 0, t) = v(x, 1, t)=0,

На верхней границе полости задаётся касательная скорость u(x, 1, t) = u0 sin(0.1π t),

что приводит к возникновению вихревого движения внутри полости. На остальных границах задаётся условие свободного проскальзывания:

∂v	(0, y,t) =	∂v	(0, y,t) = 0 ,	∂u	(x,0,t) = 0 .
∂x		∂x		∂ y

На рисунке 8.2 показано решение этой задачи в момент времени t=15. Течение жидкости происходит вдоль изолиний функции тока.

figure8_2.tif

Рисунок 8.2. Изолинии функции тока для течения в полости (Re=103, N=M=32).

______________________________________________________________________________________

122

8.3 Метод маркеров и ячеек.

Применение переменных функция тока-завихрённость бывает не всегда удобно для решения задач о течении вязкой несжимаемой жидкости. Так в течениях со свободными поверхностями, использование этих переменных сильно затрудняет постановку граничных условий. Использование в качестве неизвестных функций давления и компонент скорости позволяет избежать многих сложностей, особенно в трёхмерных задачах. Поэтому в дальнейшем мы будем работать с уравнениями (8.6). Одним из часто используемых методов является метод маркеров и ячеек (MAC, Marker and Cell). Первоначальная идея была предложена Харлоу и Уэлчем, которая затем получила дальнейшее развитие в различных модификациях метода MAC. В этом методе производится интегрирование уравнений движения несжимаемой вязкой жидкости с одновременным рассмотрением маркеров-частиц, перемещающихся со скоростью жидкости. Мы обсудим основные особенности этого подхода на примере решения двумерных задач в прямоугольной области.

Введём разностную сетку с шагами hx= xn+1-xn и hy= ym+1-ym. Параметры жидкости будут определяться в “смещённых” точках, как показано на рисунке 8.3.

Рисунок 8.3. Расчётная ячейка в методе MAC. В точках «•» вычисляется компонента скорости u; в точках «×» - компонента v; давление вычисляется в точках « ».

Предполагая, что давление в момент времени tk известно, на первом этапе вычисляется промежуточное поле скоростей на основе следующей аппроксимации уравнений (8.6):

123

unk,+m1/+1/2 2 − unk,m+1/ 2

+

(u2

)nk+1/ 2,m+1/ 2 −

(u2 )nk−1/ 2,m+1/ 2

+

τk

hx

(u v )nk,m+1 − (u v )nk,m

= −

pnk+1/ 2,m+1/ 2 − pnk−1/ 2,m+1/ 2

+

hy

hx

1

(

xunk,m+1/ 2 + yunk,m+1/ 2 )+ fe,1(xn , ym+1/ 2 ,tk )

Re

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

vnk++1/1/2,2 m − unk+1/ 2,m + (v2 )nk +1/ 2,m+1/ 2 −(v2 )nk +1/ 2,m−1/ 2 +

τk

hy

(u v )nk +1,m − (u v )nk,m

= −

pnk+1/ 2,m+1/ 2 − pnk+1/ 2,m−1/ 2

+

hx

hy

1

(

xvnk+1/ 2,m + yvnk+1/ 2,m )+ fe,2 (xn+1/ 2 , ym ,tk )

Re

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

где

(u2 )k

= unk+1,m+1/ 2 unk,m+1/ 2 , (v2 )k

= vnk+1/ 2,m+1 vnk+1/ 2,m ,

n+1/ 2,m+1/ 2

n+1/ 2,m+1/ 2

(u v )k

=

1

(unk,m+1/ 2 + unk,m−1/ 2 )(vnk+1/ 2,m + vnk−1/ 2,m ),

4

n,m

unk,m+1/ 2

1

unk+1,m+1/ 2 −

2unk,m+1/ 2

+ unk−1,m+1/ 2

=

,

h2

x

v

k

v

k

− 2v

k

+ v

k

x

n+3 / 2,m

n+1/ 2,m

n−1/ 2,m

n+1/ 2,m

unk,m+1/ 2

1

unk,m+3 / 2 − 2unk,m+1/ 2 + unk,m−1/ 2

y

=

.

hy2

−

2vk

+ vk

vk

vk

n+1/ 2,m

n+1/ 2,m+1

n+1/ 2,m

n+1/ 2,m−1

(8.25)

(8.26)

Соотношения (8.25) и (8.26) необходимо ещё дополнить граничными условиями, записанными в разностной форме. Рассмотрим, например, расчёт граничных условий, заданных на вертикальной границе x=0. Тогда вектор нормали n=(1,0), вектор касательной k=(0,1) и условия (8.3)-(8.5) принимают вид

u(0, y) = gn ( y,t) ,	(8.27)
v(0, y) = gk ( y,t) ,	(8.28)
или
∂v (0, y) = 0 .	(8.29)
∂ x

Разностная форма условия (8.27) записывается очень просто:

u0,k +m1/+1/2 2 = gn ( ym+1/ 2 ,tk +1) , m=0, … , M-1.

124

Условие (8.28) или (8.29) используется для определения скорости v в точке (1/2, m) через расчёт величин

(u v )0,k m и xv1/k 2,m .

Эти величины включают в себя мнимое значение v-1/2,m. Для исключения этого значения в случае условия (8.28) можно положить

12 (v1/k 2,m + v−k1/ 2,m )= gk ( ym ,tk ) v−k1/ 2,m = 2gk ( ym ,tk ) − v1/k 2,m ,

а в случае условия (8.29) –

vk	− vk
1/ 2,m	−1/ 2,m	= 0 v−k1/ 2,m = v1/k 2,m .
	hx

Расчёт условий на других границах производится аналогичным образом.

Исследование устойчивости схем (8.25), (8.26) проводится так же как и для схемы (8.13). В результате мы получим следующее условие

							4(γ3 + γ4 )													1			,
τk	≤ min																,
			4(γ		3	+ γ		4	)2	+ (γ		+ γ		2	)2			2		(γ3 + γ4	)
											1
где
	γ1 =	1		max | unk,m | ,							γ2	=		1		max | vnk,m | ,
		hx n,m											hy n,m
				γ3 =			1			,	γ4	=			1					.
							Re h2						Re h2
									x								y

Вычисленные промежуточные значения скорости могут не удовлетворять условию неразрывности, записанному для каждой расчётной ячейки, то есть

	k +1/ 2		uk +1/ 2	− uk +1/ 2		vk +1/ 2	− vk +1/ 2
			n+1,m+1/ 2	n,m+1/ 2		n+1/ 2,m+1	n+1/ 2,m		,
	dn+1/ 2,m+1/ 2	=			+			≠ 0
			hx			hy

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

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

	ϕk +1	− ϕk +1
unk,+m1+1/ 2 = unk,+m1/+1/2 2 −	n+1/ 2,m+1/ 2	n−1/ 2,m+1/ 2	,
		hx

(8.30)

125

			ϕk +1		− ϕk +1
vnk++1/1 2,m = vnk++1/1/2,2 m −			n+1/ 2,m+1/ 2		n+1/ 2,m−1/ 2	.
					hy
Тогда
k +1		k +1/ 2	−
dn+1/ 2,m+1/ 2	= dn+1/ 2,m+1/ 2
		ϕnk++13/ 2,m+1/ 2 − 2ϕnk++1/1 2,m+1/ 2 + ϕnk−+1/1 2,m+1/ 2					−.
				h2
				x
		ϕnk++1/1 2,m+3/ 2 − 2ϕnk++1/1 2,m+1/ 2 + ϕnk++1/1 2,m−1/ 2
				hy2
Из требования
		dnk++1/1 2,m+1/ 2 = 0
следует разностное уравнение Пуассона для вычисления функции ϕ:
		k +1		k +1/ 2			(8.31)
		P(ϕn+1/ 2,m+1/ 2 )		= dn+1/ 2,m+1/ 2 .

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

( hϕ,n ) = 0

во всех точках границы, здесь h - разностная форма оператора .

Давление в момент времени tk+1 вычисляется по известным значениям функции ϕ. Так соотношение (8.25) можно представить в виде

k +1/ 2	k		∂ p k	.
un,m+1/ 2	= un,m+1/ 2	+ τk (...)− τk
			∂x n,m+1/ 2

Корректировку (8.30) можно трактовать как введение некоторого дополнительного поля давления, то есть

unk,+m1+1/ 2

uk +

n,m 1/ 2

Тогда

k +1/ 2		∂ϕ k +1						=
= un,m+1/ 2		−	∂x
						n,m+1/ 2
+ τk (...)− τk					∂ p k
					∂x n,m+1/ 2
∂ϕ		= −τk				∂ p
−						∂x
∂ x

	∂ p k +1			.
+	∂ x
			n,m+1/ 2

p = ϕ τk

и давление в момент времени tk+1 вычисляется как

	ϕk +1
pnk++1/1 2,m+1/ 2 = pnk+1/ 2,m+1/ 2 +	n+1/ 2,m+1/ 2	,
	τk

126

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

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

	d xp	= u (xp , yp ,t ),	d yp		= v(xp , yp ,t )	(8.32)
	d t		d t
		p=1, …, Np ,
		xp (0 ), yp (0 )	- заданы,

где (xp, yp) - координаты маркера в декартовой системе координат.

Для получения разностной формы уравнений (8.32) используется линейная интерполяция скорости как по пространству так и по времени для каждой разностной ячейки. Пусть, например, маркер с номером p находится в момент времени tk в ячейке

{xn≤x≤xn+1, ym≤y≤ym+1} и значения скорости в узлах (n, m+1/2), (n+1, m+1/2), (n+1/2, m) и (n+1/2, m+1) известны как в момент времени tk, так и в момент времени tk+1. Тогда положение этого маркера в момент времени tk+1 вычисляется по формулам

p

p

2h

(

n,m+1/ 2

n,m+1/ 2 )(

n+1

p )

xk

+1

= xk +

τk

uk +1

+ uk

x

− xk

+

x

,

(

n+1,m+1/ 2

n+1,m+1/ 2

)(

p

n )

uk +1

+ uk

xk

− x

p

p

τk

(

n+1/ 2,m

n+1/ 2,m )(

m+1

p )

yk

+1

= yk +

2hy

vk +1

+ vk

y

− yk

+

.

( n+1/ 2,m+1

n+1/ 2,m+1 )(

p

m )

vk +1

+ vk

yk

− y

Положение жидкости определяется распределением маркеров по ячейкам разностной сетки. При этом выделяется три типа ячеек: пустая ячейка (E), которая не содержит маркеров; поверхностная ячейка (S), которая содержит частицы-маркеры и граничит с пустыми ячейками и ячейки жидкости (F) (смотри рисунок 8.4).

127

Рисунок 8.4. Частицы-маркеры определяют ячейки занятые жидкостью.

Давление в пустых ячейках равно pвнеш, которое для простоты мы будем считать равным нулю, так как мы всегда можем этого добиться путём замены p→p-pвнеш. Тогда граничным условием на свободной поверхности будет отсутствие нормального и касательного напряжений. Если n=(nx, ny) есть вектор нормали к свободной поверхности в некоторой точке, то эти условия имеют вид

2

∂u

2

∂u

∂v

∂v

2

= 0 ,

p −

nx

+

+

nxny +

ny

∂ x

∂ y

∂x

∂ y

Re

(8.33)

1

∂v

∂u

∂u

∂v

2

2

2

−

nxny +

+

(nx − ny

)

= 0 .

∂ y

∂x

∂ y

Re

∂x

Расчёт давления и компонент скорости в поверхностных ячейках осуществляется специальным образом и основывается на граничных условиях (8.33) и условии неразрывности. При этом выделяется четыре типа поверхностных ячеек, а именно:

1)ячейки граничащие с одной пустой ячейкой,

2)ячейки граничащие с двумя пустыми ячейками,

3)ячейки граничащие с тремя пустыми ячейками, и

4)ячейки окружённые пустыми ячейками и представляющие собой изолированные капли. Наличие свободной поверхности вносит свои коррективы и в вычисление потенциальной

функции ϕ. Для ячеек занятых жидкостью справедливо уравнение (8.31) с граничным

128

условием ϕ=0 в поверхностных ячейках. При этом область определения уравнений (8.31) задаётся положением частиц жидкости, а значит зависит от времени. Если эти уравнения дополнить соотношением ϕ=0 в пустых ячейках, то мы получим систему уравнений с переменными коэффициентами и с постоянной областью определения, которая совпадает с областью определения задачи. Для решения полученных уравнений обычно применяются итерационные методы. Так как поле скоростей меняется незначительно в течении одного шага по времени, то в качестве хорошего начального условия для вычисления ϕ в момент времени tk+1 можно взять значение этой функции в момент времени tk.

Суммируя всё изложенное выше, расчёт одного шага по времени для течениё со свободными границами производится следующим образом:

1)по известному в момент времени tk полю скоростей определяется шаг по времени τk,

2)по известным в момент времени tk координатам маркеров выделяются ячейки жидкости

(F) и поверхностные ячейки (S),

3)вычисляются промежуточные скорости в ячейках F,

4)решается уравнение Пуассона для определения функции ϕ,

5)в ячейках F вычисляются скорости и давление в момент времени tk+1,

6)вычисляются скорости и давление в ячейках S,

7)определяются новые положения маркеров в момент времени tk+1.

129