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

Это значит, что разностная схема (11.52) аппроксимирует дифференциальную задачу (11.51) со вторым порядком по h и τ. Далее, подставляя в разностную схему решение в виде (11.14), получим уравнение для собственных значений оператора перехода

λ2 −2λ + 1 = r2 (exp(iα) −2 + exp(−iα)) = −4r2 sin2(21 α) .

Решение этого уравнения имеет вид

и, согласно критерию фон Неймана, схема «крест» будет устойчива при r≤1.

________________________________________________________________________________

Пример 11.4 (схема «крест»)

Рассмотрим задачу (11.51) с нулевыми начальными условиями и граничными условиями

Решение, полученное с помощью схемы (11.52) показано на рисунке 11.13. Как показывают эти результаты, данная схема сохраняет амплитуду, но искажает форму волны. Более того, появляются высокочастотные осцилляции, отсутствующие в точном решении.

41

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

τ=0.0333, r=2/3).

________________________________________________________________________________

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

11.3.2. Диссипация и дисперсия сеточного волнового решения.

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

u(x,t)=u0exp(i(ωt-kx)),

где ω=2πν - круговая частота, k=2π/λ - волновое число. Если мы подставим это решение в какое-нибудь линейное волновое уравнение, то получим зависимость ω=ω(k), которое называется дисперсионным соотношением. Если ω принимает вещественные значения, то амплитуды гармоник не изменяются со временем. В противном случае, когда ω принимает

42

комплексные значения, волна будет затухать (диссипировать) как exp(-(Imω)t)=exp(-γt). Отношение

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

d

dk Re ω = cg

называется групповой скоростью. Она может интерпретироваться как скорость волнового пакета, состоящего из гармонических волн с близкими волновыми числами. Передача энергии в волне осуществляется с групповой скоростью. Если фазовая (групповая) скорость зависит от k, то гармоники с разными волновыми числами распространяются с разными скоростями. Это явление называется дисперсией. Для уравнения (11.51) нетрудно получить, что ω=ck и, следовательно, cf=cg=c. Это значит, что точное решение уравнения (11.51) представляет собой волну, которая распространяется без затухания и дисперсии. Другой пример: для линеаризованного уравнения Бюргерса

зависимость ω(k) имеет следующий вид

ω(k)=ck+iμk2.

Это соотношение описывает затухающую недиспергирующую волну. В частности, при μ=0 мы получаем уравнение переноса (11.1), для которого справедливо соотношение cf=cg=c и коэффициент затухания равен нулю.

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

Анализ диссипации и дисперсии сеточного решения мы будем проводить на решении unk = u exp(i(ωtk −kxn )).

43

Если подставить это решение в разностную схему, то мы получим некоторое соотношение вида ω=ω(k,τ,h). Сравнивая это соотношение с соответствующим дисперсионным соотношением для дифференциального уравнения, которое эта схема аппроксимирует, мы можем анализировать свойства данной схемы. В качестве примера рассмотрим свойства схемы «крест». Подставим приведённое выше решение в схему (11.52) и получим следующее

выражение exp(iωτ)-2+exp(-iωτ)=r2(exp(-ikh)-2+exp(ikh)) из которого следует, что sin(12 ωτ) = r sin(21 kh).

В данном случае ω принимает вещественные значения, и значит, сеточное волновое решение не затухает со временем. Выражая ω из предыдущего выражения, мы приходим к следующему дисперсионному соотношению:

ω = τ2 arcsin(r sin(21 kh)).

Тогда групповая скорость волны распространяющейся на сетке, определяется как

то есть, в отличии от уравнения (11.1), возникает дисперсия. Только когда r=1, групповая скорость приближённого решения совпадает с групповой скоростью точного решения.

Минимальная длина волны, которая может быть представлена на сетке, составляет λ=4h (kh=π/2). Понятно, что гармоники с такой длиной волны сильнее всего подвержены влиянию дисперсии. Для этого случая, зависимость безразмерного параметра a=сm/c от параметра Куранта показана на рис. 11.14.

44

Рисунок 11.14. Зависимость параметра a=сm/c от параметра Куранта для схемы «крест»;

- безразмерная групповая скорость волны cg/c=1 для дифференциального уравнения;

- параметр a при kh=π/2.

Зависимость, показанная на этом рисунке, представляет наихудшую ситуацию. Для длин волн много больших пространственного шага, дисперсия приближённого решения уменьшается.

11.3.3. Схема Лакса-Вендроффа.

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

Здесь u – скорость частиц среды, p – давление, ρ - плотность среды, и c – скорость звука.

Уравнения (11.54) дополняются начальными условиями

u(x, 0) = g1(x) , , p(x, 0) = g2(x)

и граничными условиями

α1u(0,t) + α2p(0,t) = f1(t) , . α3u(l,t) + α4p(l,t) = f2(t).

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

момент времени t+ t можно получить с помощью разложения функции s(x,t+ t) в ряд Тейлора в точке (x,t) по степеням t, то есть

Производные по t заменяются производными по x путём использования системы (11.4):

45

Рисунок 11.15. Вспомогательные узлы «×» для вычисления промежуточного решения.

Здесь

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

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

Выражения (11.57) и (11.58) называются двухшаговой схемой Лакса-Вендроффа. Легко проверить, что в частном случае, когда матрица A постоянна, эта схема переходит в схему

(11.56).

Перейдём к анализу устойчивости схемы (11.6). Только теперь мы имеем систему

уравнений, поэтому вместо (11.14) зададим решение в виде skn = λk (α)exp(iαn)s ,

где s* - некоторый постоянный вектор. Подстановка этого решения в разностную схему (11.56) даёт следующее выражение:

После некоторых преобразований, мы получим стандартную задачу на собственное значение

47

Решение характеристического уравнения det(B-λI)=0 даёт следующее выражение для

собственного значения оператора перехода

λ(α) = 1 −r2(1 −cos(α)) ± ir sin(α) .

Тогда из критерия фон Неймана следует, что разностная схема (11.6) будет устойчива при r≤1.

Анализ диссипативных и дисперсионных свойств разностных схем для системы

уравнений (11.4) мы будем проводить на решении

skn = s exp(i(ωtk −kxn )) .

Если подставить это решение в разностную схему (11.6),то мы получим систему уравнений

Из условия det(D)=0 следует дисперсионное соотношение

exp(iωτ) = 1 −r2 (1 −cos(kh))+ ir sin(kh)

и, следовательно, ω принимает комплексные значения, что говорит о диссипации сеточного решения. Введём обозначение ω=ω0+iγ, и приравнивая вещественные и мнимые части в

предыдущем выражении, получим

exp(−γτ)cos(ω0τ) = 1 −r2 (1 −cos(kh)), . exp(−γτ)sin(ω0τ) = r sin(kh).

Отсюда нетрудно получить коэффициент затухания

и групповую скорость волны на сетке

На рисунках 11.16 и 11.17 представлены зависимости параметра диссипации d=γτ и групповой скорости от параметра Куранта.

48

Рисунок 11.16. Зависимость параметра диссипации d=γτ от параметра Куранта для схемы Лакса-Вендроффа (kh=π/2).

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

49

Для аппроксимации граничных условий со вторым порядком по h и τ, можно использовать метод фиктивных областей. Пусть, например, задано граничное условие при x=0: u(0,t)=f(t). В разностном виде это граничное условие записывается очень просто:

= f (tk+1) . Для вычислений по схеме (11.56) нам необходимо также значение давления в граничном узле. Для определения этого значения введём фиктивный слой с координатами

(x-1,tk) и запишем второе уравнение (11.56) в виде

p0k+1 = p0k − rc2ρ(u1k −u−k 1 )+ 21 r2 (p1k −2p0k + p−k 1 ).

Для исключения фиктивных значений u−k1 и p−k1 используем очевидное приближение

0.5(u1k + u−k 1) = u0k = f (tk ) .

Тогда первое уравнение из (11.56) примет вид

2rρc (p1k − p−k 1 ) = u0k −u0k+1 = f (tk ) − f (tk+1) .

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

pk+1 = pk −ρc r (uk − f (t + ))+ r2 (pk − pk ).

0 0 1 k 1 1 0

В случае, когда на границе задано давление, аналогичная процедура применяется для определения скорости в граничной точке.

________________________________________________________________________________

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

Рассмотрим задачу (11.54) (которая эквивалентна задаче (11.51)) с нулевыми начальными условиями и граничными условиями (11.53). Приближённое решение, полученное с помощью схемы Лакса-Вендроффа (11.56) показано на рисунке 11.18.

50