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

Рассмотрим бесконечную пластину толщины l. Пусть в момент времени t=0 она имеет температуру u0, а температура окружающей среды поддерживается равной нулю. Тогда распределение температуры в пластине описывается уравнением (11.1) с условиями

u(x,0)=u0 (начальное условие), (11.29) u(0,t)=u(l,t)=0 (граничные условия).

Точное решение этой задачи имеет вид

	4u0	∞	exp(−κπ2(2m + 1)2t l2 )			π(2m	+ 1)x
u(x,t) =		∑			sin			.	(11.30)
	π m=0		2m + 1			l

На рисунке 11.2 показано распределение температуры в момент времени T=0.1 при κ=1, l=1 и u0=1. Для оценки ошибки приближённого решения введём вектор приближённого

решения в момент времени tk

vak = (u0k ,...,uNk )T

и точное решение, вычисленное в узлах разностной сетки в момент времени tk

vek = (u(x0,tk ),...,u(xN ,tk ))T .

Тогда относительная ошибка приближённого решения определяется как

e(vak ) =	|| vek − vak ||2 .
	|| vek ||2

Для решения, показанного на рисунке 11.8 ошибка e(vak ) = 0.0347 .

Рисунок 11.8. Решение уравнения (11.19) с условиями (11.29); - точное решение (11.30),- приближённое решение ( явная схема (11.20), h=0.1, τ=0.005).

21

На рисунке 11.9 показано распределение температуры в момент времени T=0.1, полученное с помощью схемы (11.21). При этом значение τ выбиралось из условия (11.22), а шаг по времени выбирался из условия β = κτ/h2 = 1.66666 . Для этого решения e(vak ) = 0.0136 , что демонстрирует эффективность применения неявной схемы для решения такого типа задач.

Рисунок 11.3. Решение уравнения (11.19) с условиями (11.29); - точное решение (11.30),

- приближённое решение ( схема (11.21), h=0.1, τ=1.66666 10-2, γ=0.45).

________________________________________________________________________________

11.2.2. Построение разностных схем для одномерного уравнения теплопроводности с переменными коэффициентами.

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

ρ(x)c(x) ∂∂ut = ∂∂x

Введём обозначение для теплового потока

q =

a(x

a(x)

	∂u
)		, 0≤t≤T, 0≤x≤l .	(11.31)
	∂x

∂∂ux ,

22

тогда предыдущее уравнение примет вид
ρ(x)c(x)	∂u	=	∂q .	(11.32)
	∂t		∂x
Введём две сетки: Sm с узлами (xn, tk) (основная сетка)				и Sa с узлами (xn+1/2, tk)

(вспомогательная сетка), как показано на рисунке 11.10. Далее, температура будет определяться в узлах основной сетки, а потоки в узлах вспомогательной сетки.

Рисунок 11.10. Основная сетка ; вспомогательная сетка .

Проинтегрируем уравнение (11.32) по разностной ячейке Sa:

xn +1/2 tk +1

∂u

tk +1

xn +1/2

∂q

∫ ∫

ρ(x)c(x)

=

dt dx

dx dt .

∂t

∫ ∫

∂x

xn−1/2 tk

tk xn−1/2

Вычисление интегралов даёт

xn +1/2

tk +1

∫

ρ(x)c(x) u(x,t

k+1

) −u(x,t )

dx =

q(x

n+1/2

,t) −q(x

n−1/2

,t) dt .

(11.33)

(

k )

∫ (

)

xn−1/2

tk

Это выражение представляет собой закон сохранения энергии, записанный для разностной ячейки. Такая форма записи является более универсальной, так как, например, если функция a(x) разрывна в некоторых точках отрезка 0≤x≤l, то уравнение (11.31) в этих точках смысла не имеет. Для того чтобы построить разностную схему, необходимо сделать некоторые предположения о распределении u(x,tk) на отрезке [xn-1/2, xn+1/2] и q(xn-1/2,t) на отрезке [tk, tk+1]. Пусть u(x,tk) на каждом отрезке [xn-1/2, xn+1/2] принимает постоянное значение, которое мы

обозначим unk , а значение потока q(xn+1/2,t) постоянно на каждом интервале [tk, tk+1] и равно

q(xn+1/2,t) ≈ γq(xn+1/2,tk+1) + (1 − γ)q(xn+1/2,tk ) = γqnk++1/1 2 + (1 − γ)qnk+1/2 .

Введём обозначения

23

1 xn +1/2

(ρc)n = h ∫ ρ(x)c(x) dx , an+1/2 = a(xn+1/2 ),

xn−1/2

тогда выражение (11.33) принимает вид

(ρc)n (unk+1 −unk )h = (γ(qnk++1/1 2 −qnk−+1/1 2 )+ (1 − γ)(qnk+1/2 −qnk−1/2 ))τ .

Нам осталось выразить тепловой поток через температуру. Аппроксимируя выражение для теплового потока центральной разностью относительно точки xn+1/2, получим

qnk+1/2 = an+1/2	uk	−uk
	n+1	n	.
		h

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

k

+1

k

k+1

k+1

k+1

k+1

un

−un

un+1

−un

un

−un−1

(ρc)

=

γ a

−a

n−1/2

+

2

2

n

τ

n+1/2

h

h

(11.34)

k

k

k

k

un+1 −un

un

−un−1

(1 − γ) a

−a

.

n+1/2

2

n−1/ 2

2

h

h

Легко показать, что когда ρ, c и a постоянны, эта схема переходит в схему (11.21).

В случае уравнения с переменными коэффициентами аналитическое исследование аппроксимации не всегда возможно. Допустим, что a(x) является дважды непрерывно дифференцируемой функцией. Тогда можно показать, что, в общем случае, разностная схема (11.34) аппроксимирует дифференциальное уравнение (11.31) с точностью до членов O(τ+h2), а если положить γ=1/2 , то с точностью до членов O(τ2+h2). В случае кусочнонепрерывной функции a(x), о порядке аппроксимации нельзя ничего сказать заранее, так как он будет зависеть от решения конкретной задачи.

Для анализа устойчивости применим принцип замороженных коэффициентов. Заменим переменные параметры (ρc)n, an-1/2 и an+1/2 в схеме (11.34) на некоторые константы c1, c2 и c3, соответственно. Далее исследуем полученную схему на устойчивость. В результате мы прийдём к условию (11.23) где κ=(c2+c3)/(2c1). Учитывая, что значения коэффициентов зависят от номера сеточного узла, мы окончательно получим следующее условие устойчивости схемы (11.34):

	τ			an−1/2 + an+1/2			1
										1
β =	h2	max		(ρc)		≤	1 −2γ	, 0	≤ γ <		2 .
		n
				n

(11.35)

β - любое>0, 12 ≤ γ ≤1

24

Для неявной схемы (γ≠0) вычисления производятся также как и для схемы (11.21), только коэффициенты матрицы A из (11.28) определяются как

pn = −γ aρn−1/ 2τ2 ,

( c)n h

sn = 1 + γ (an−1/ 2 + an+1/ 2 )τ , ,

(ρc)n h2

rn = −γ aρn+1/2τ2 ,

( c)n h

n=1, …, N-1.

Для нелинейной задачи

∂		∂		∂u
	(ρ(u)c(u) u) =		a(u)		, 0	≤ t ≤ T, 0	≤ x ≤ l ,
∂t
		∂x		∂x

разностная схема строится аналогичным образом. Только для данной задачи следует использовать явную аппроксимацию с переменным шагом по времени, так как неявная схема требует решения системы нелинейных уравнений на каждом временном шаге, что может потребовать неоправданно большого объёма вычислений. Учитывая это замечание, преобразуем схему (11.34) к следующему виду:

ρkck

uk+1

−uk

uk

−uk

n

n = ak

n+1

n

−ak

h2

n n

τk

n+1/2

n−1/ 2

где

ρk

= ρ(uk ) ,

n

n

ck

= c(uk ) ,

.

n

n

ak

= a

1

(uk

+ uk ) .

(2

n+1/2

n+1

n )

uk	−uk	,
n	n−1
	h2

Условие устойчивости для данной схемы следует из условия (11.35) при γ=0:

τk				k		k
			an−1/ 2 + an+1/ 2
	2	max						≤1 .
					k	k
h		n			ρ c
					n	n

________________________________________________________________________________

Пример 11.2 (распределение температуры в неоднородной среде при воздействии теплового потока).

Рассмотрим задачу о воздействии теплового потока на границу неоднородного полупространства. Тогда при x>0 распределение температуры будет описываться уравнением (11.13). Начальные и граничные условия зададим в виде

25

u(x, 0) = 0 , x>0 ,

a(0)	∂u (0,t)	= −f , t > 0 ,	(11.36)
	∂x	0
	lim u(x,t) = 0 .
	x→∞
Пусть зависимость теплофизических параметров от координаты имеет следующий вид:
ρ(x)c(x) = ρ0c0(1 + x) , a(x) = a0(1 + x).			(11.37)

Такая зависимость характерна, например, для грунтов.

Для численного решения рассматриваемой задачи применим явную схему (11.34), где

(ρc)n = ρ0c0(1 + xn ) , an+1/2 = a0(1 + xn + 0.5h) .

Легко проверить, что в данном случае условие устойчивости (11.35) совпадает с условием устойчивости для однородного уравнения: β≤1/2.

Применяя преобразование Лапласа по времени, можно получить точное решение

рассматриваемой задачи:			f0 κK0 ((1 + x) s /κ)
	u	(x,s) =			,
			s3/2K1 ( s /κ)

где κ=a0/(ρ0c0), s – переменная преобразования Лапласа, K0 и K1 - модифицированные функции Бесселя второго рода (функции Макдональда) порядка 0 и 1, соответственно. Обращение этого выражения даёт точное решение задачи u(x,t).

На рисунке 11.11 показано распределение температуры в момент времени T=0.5, полученное с помощью схемы (11.16), вместе с точным решением при ρ0c0=1, a0=1, f0=1.

26

Рисунок 11.11. Решение уравнения (11.31) с условиями (11.37) и (11.38); - точное

решение, - приближённое решение (явная схема (11.34), h=0.25, τ=0.0278).

________________________________________________________________________________

11.2.3. Метод расщепления по пространственным переменным

Метод расщепления является одним из важнейших инструментов для построения разностных схем для многомерных задач математической физики. Мы обсудим основные идей этого метода на примере двумерного уравнения теплопроводности, но аналогичные рассуждения справедливы и для других типов задач, например, для уравнений акустики и газовой динамики.

Рассмотрим следующую задачу

∂u		2			2
		∂ u		+	∂ u			(11.38)
	= κ		2			2	,
∂t		∂x			∂y

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

u(x,y, 0) = g(x,y) ,

b u(0,y,t) +b ∂u

(0,y,t) = f

(y,t) , b u(l

,y,t) +b ∂u

(l

,y,t) = f

(y,t) ,

1

2 ∂x

1

3

x

4 ∂x x

2

b u(x, 0,t) +b ∂u

(x, 0,t) = f

(x,t) , b u(x,l

,t) +b

∂u

(x,l

,t) = f

(x,t) .

5

6 ∂y

3

7

y

8 ∂y

y

4

Запишем дифференциальное уравнение этой задачи в операторном виде:

27

∂u = Du = D u + D u = κ ∂2u				+ κ ∂2u .
∂t	x	y	∂x2	∂y2

Пусть решение в некоторый момент времени tk известно. Тогда решение в момент времени tk+τ можно выразить через это известное решение, используя разложение в ряд Тейлора:

u(x,y,t + τ) = u(x,y,t		) + τ ∂u	(x,y,t ) +O(τ2 ) =
k	k	∂t	k	(11.39)

(I + τD)u(x,y,tk ) +O(τ2 ),

здесь I – единичный оператор (Iu=u). Расщепим исходную задачу (11.38) на две одномерные задачи

∂v = D v

,

(11.40)

∂t

x

tk ≤ t ≤ tk + τ, 0 ≤ x ≤ lx , 0 ≤ y ≤ ly,

v(x,y,tk ) = u(x,y,tk ),

b v(0,y,t) +b

∂v

(0,y,t) = f (y,t) , b v(l

,y,t) +b

∂v

(l

,y,t) = f

(y,t) ,

1

2 ∂x

1

3

x

4 ∂x

x

2

и

∂w

= D w ,

(11.41)

∂t

y

tk ≤ t ≤ tk + τ, 0 ≤ x ≤ lx , 0 ≤ y ≤ ly,

w(x,y,tk ) = v(x,y,tk + τ),

b w(x, 0,t) +b

∂w

(x, 0,t) = f

(x,t) , b w(x,l

,t) +b ∂w

(x,l

,t) = f

(x,t).

5

6 ∂y

3

7

y

8 ∂y

y

4

Последовательное решение этих задач даёт решение w(x,y,tk+τ) и, естественно, это решение

должно совпадать с u(x,y,tk+τ). По аналогии с (11.39) можно показать, что v(x,y,tk + τ) = (I + τDx )v(x,y,tk ) +O(τ2 )

и

w(x,y,tk + τ) = (I + τDy )w(x,y,tk ) +O(τ2 ) .

Тогда из (11.39),(11.40) и (11.41) следует, что

28

w(x,y,tk + τ) = (I + τDy )v(x,y,tk + τ) +O(τ2 ) =

(I + τDy )(I + τDx )v(x,y,tk ) +O(τ2 ) = . (I + τ(Dx + Dy ))u(x,y,tk ) +O(τ2 ) =

(I + τD)u(x,y,tk ) +O(τ2 ) = u(x,y,tk + τ) +O(τ2 )

Таким образом, мы показали, что решение, полученное после последовательного решения задач (11.40) и (11.41) совпадает с точностью до членов O(τ2) с решением задачи (11.38). Это даёт возможность заменить многомерную задачу последовательностью более простых одномерных задач.

11.2.4. Разностные схемы для решения многомерного уравнения теплопроводности.

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

∂u		2			2
		∂ u		+	∂ u			+ fs (x,y,t),	(11.42)
	= κ		2			2
∂t		∂x			∂y

0 ≤ t ≤ T, 0 ≤ x ≤ lx , 0 ≤ y ≤ ly, u(x,y, 0) = g(x,y) , +граничные условия.

Чтобы применить метод конечных разностей в этой ситуации, поступим совершнно так же, как и в одномерном случае. С этой целью построим множество равноотстоящих точек xn (n=0,…,N) на отрезке 0≤x≤lx с x0=0, xN=lx, xn+1-xn=h, а также множество равноотстоящих точек ym (m=0,…,M) на отрезке 0≤y≤ly с y0=0, yM=ly, ym+1-ym=h. После этого область, где ищется решение, покрывается прямоугольной пространственной сеткой: через каждую точку xn проводятся прямые параллельные оси y, а через каждую точку ym – прямые параллельные оси x (рисунок 11.12).

29

Рисунок 11.12. Двумерная пространственная сетка.

В дополнение введём множество точек по времени tk (k=0,1,…) с tk+1-tk=τ. Тогда типичная точка разностной сетки имеет координаты вида (xn, ym, tk). После этого входящие в уравнение (11.25) частные производные заменяются соответствующими конечноразностными аппроксимациями.

Для более компактной записи разностных схем, примем следующие обозначения для разностных операторов аппроксимирующих вторые производные по x и y

	xunk,m =		1			(unk+1,m −2unk,m + unk−1,m ),
			h2
	yunk,m =	1				(unk,m+1 −2unk,m + unk,m−1 ).
		h2
Тогда двумерный аналог схемы (11.21) записывается как
	unk,+m1 −unk,m			= κγ( xunk,+m1 + yunk,+m1 )+
							(11.43)
	τ					,

κ(1 − γ)( xunk,m + yunk,m )+ fs (xn ,ym ,tk )

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

Шаблон, схематически представляющий эту схему, показан на рис. 11.13.

30