← Назад Далее →

2.5. Частные решения линейных разностных схем. Диссипативные и дисперсионные поверхности

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

$$u_j^n=D\cdot \exp \left\{i[\omega \cdot n\tau -kjh]\right\}=A\cdot q^n\cdot \exp \{-i\left(kjh\right)\};(2.38)$$

здесь $q=\exp \left\{i\left(\omega \tau \right)\right\}$ - множитель переходана следующий временной слой.

Подставим (2.38) в разностный оператор схемы «уголок» (2.26). Получаем:

$$\begin{array}{l}q=r\cdot e^{ikh}+(1-r)\cdot ,\Rightarrow \\ \Rightarrow \left|q(r,kh)\right|=\sqrt{{[(1-r)+r\cdot \cos \left(kh\right)]}^2+{\left[\mathrm{s}\mathrm{i}\mathrm{n}\left(kh\right)\right]}^2}\end{array}(2.39)$$

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

Модуль множителя перехода $\left|q\right|$ является функцией от числа Куранта - Фридрихса-Леви $r$ и приведенного волнового числа $kh$, изменяющегося в пределах от $-\pi$: $(-\pi \le kh\le \pi )$. Он показывает, как изменяется амплитуда бегущей волны при переходе на следующий слой. Если модуль перехода меньше или равен единице, то амплитуда волны не возрастает. Если для заданного $r$ амплитуды всех волн при $(-\pi \le kh\le \pi )$ не возрастают, то говорят, что данное число $r$ принадлежит области устойчивости дискретного оператора. Опираясь на (2.39) можно показать, что область устойчивости разностной схемы «уголок» представляет собой единичный отрезок $r\in [0,1]$. Если в области устойчивости модуль перехода любой волны равен единице, то дискретный оператор называется бездиссипативным, в противном случае имеет место численная диссипация.

Диссипативные свойства дискретного оператора для всех приведенных волновых чисел наиболее наглядно можно представить в виде двумерной поверхности $\left|q\right|=f(r,kh)$. Для схемы «уголок» диссипативная поверхность изображена на рис. 6а, где по оси $z$ отложены значения модуля перехода, по оси $x$ - значения числа $r\in [0,1]$, по оси $y$ - значения приведенного волнового числа $kh\in [-\pi ,\pi ]$. Видно, что модуль перехода стремится к единице для длинноволновых гармоник, для которых $kh\approx 0$, и равен единице на границе области устойчивости, при $r=1,(-\pi \le kh\le \pi )$.

а)	б)
Рис. 6. Диссипативная и дисперсионная поверхности схемы «уголок»

Из (2.39) также следует дисперсионное соотношение:

$$q=e^{i\omega \tau }=r\cdot e^{ikh}+(1-r)\cdot ,\Rightarrow \omega =-\frac{i}{\tau }\arg [r\cdot e^{ikh}+(1-r)]$$

откуда получается выражение для приведенной фазовой скорости бегущих волн:

$$\begin{array}{l}\gamma (r,kh)=\frac{\omega }{kc}=-\frac{i}{\tau kc}\arg [r\cdot e^{ikh}+(1-r)]=\\ =-\frac{i}{r\cdot kh}\arg [r\cdot e^{ikh}+(1-r)]\end{array}$$

Приведенная фазовая скорость бегущих волн, описываемых исходным дифференциальным оператором, не зависит от номера гармоники и равна единице. Как ранее уже отмечалось, зависимость фазовой скорости от длины волны называется дисперсией. В разностном случае приведенная фазовая скорость может принимать значения как меньшие, так и большие единицы. В первом случае дисперсия называется «нормальной», во втором - «аномальной». В случае нормальной дисперсии гармоника отстает от «эталонной» гармоники, двигающейся с фазовой скоростью $c$, в случае аномальной дисперсии - опережает ее.

Дисперсионные свойства дискретного оператора представляются его «дисперсионной поверхностью». На рис. 6в приведена дисперсионная поверхность разностного оператора «уголок». Здесь по оси $z$ отложены значения приведенной фазовой скорости $\gamma$, по оси $x$- значения числа $r\in [0,1]$, по оси $y$ - значения приведенного волнового числа $kh\in [-\pi ,\pi ]$. Дисперсия стремится к нулю, что соответствует стремлению $\gamma$ к единице, для длинноволновых гармоник, для которых $kh\approx 0$ и на границе области устойчивости, при $r=1,(-\pi \le kh\le \pi )$.

Стремление к нулю дисперсии и диссипации длинноволновых гармоник является общим свойством всех дискретных операторов, аппроксимирующих исходный дифференциальный оператор. Характер касания дисперсионных и диссипативных поверхностей единичной плоскости определяется порядком аппроксимации разностного оператора. Тот факт, что при $r=1$все гармоники не диссипируют и не диспергируют, говорит о том, что разностный оператор при данном числе Куранта переносит начальный профиль без каких либо искажений.

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

Из рис. 6в видно, что дисперсия схемы «уголок» при $r<0.5$ является нормальной, а при $r>0.5$ - аномальной. При $r=0.5$дисперсия отсутствует.

Диссипативная и дисперсионная поверхности дают полное представление о свойствах разностной схемы, позволяют предсказывать поведение дискретного решения на разных начальных условиях. Рассмотрим модельную задачу (2.1) с периодическими начальными условиями на отрезке $x\in [0,2\pi ]$:

$$\begin{array}{l}\frac{\partial u}{\partial t}+c\cdot \frac{\partial u}{\partial x}=0,c=1;t\in [0,\mathrm{\infty });\\ u(x,0)=f\left(x\right);u(0,t)=u(1,t);x\in [0,2\pi ];\end{array}(2.40)$$

и двумя видами начальных условий - достаточно гладкими и разрывными.

В качестве гладкого начального условия возьмем «двойной гауссиан» с параметрами:

$$\begin{array}{l}f\left(x\right)=\exp [-{(x-x_1)}^2/{\mathrm{\Delta }}^2]+\exp [-{(x-x_2)}^2/{\mathrm{\Delta }}^2];\\ x_1=0.5\cdot \pi ,x_2=0.7\cdot \pi ,\mathrm{\Delta }=0.1\cdot \pi \end{array}(2.41)$$

в качестве разрывного - «ступеньку»:

$$f\left(x\right)=\{\begin{array}{l}1,\mathrm{i}\mathrm{f}0.2\cdot \pi \le x\le 0.4\cdot \pi \\ 0,\mathrm{i}\mathrm{f}(0\le x\le 0.2\cdot \pi )\mathrm{\&}(0.4\cdot \pi \le x\le 2\pi )\end{array}(2.42)$$

Аппроксимируем дифференциальную задачу (2.40) - (2.42) разностной схемой «уголок». Разобьём отрезок $x\in [0,2\pi ]$ на N=100 равных интервалов с шагом $h=0.2\pi$ и спроектируем на нее начальные данные (2.41) - (2.42):

$$\begin{array}{l}{f}_i^{\left(1\right)}=\exp \{-{[(i-1)h-x_1]}^2/{(0.1\cdot \pi )}^2\}+\\ +\exp \{-{[(i-1)h-x_2]}^2/{(0.1\cdot \pi )}^2\};\\ f_i^{\left(2\right)}=\{\begin{array}{c}1,\mathrm{i}\mathrm{f}11\le i\le 21\\ 0,\mathrm{i}\mathrm{f}(1\le i\le 11)\mathrm{\&}(21\le i\le 101)\end{array};i=1,...,101\end{array}(2.43)$$

Разложим эти сеточные функции в ряд Фурье:

$$f_n^{\left(m\right)}=\sum _{k=0}^N{a}_k^{\left(m\right)}\cdot \exp [-i\left(kh\right)\cdot n];m=1,2;n=1,...,N+1(2.44)$$

где коэффициенты разложения определяются по формуле:

$$a_k^{\left(m\right)}=\frac{1}{N+1}\cdot \sum _{n=1}^{N+1}{f}_n^{\left(m\right)}\cdot \exp [i\left(kh\right)\cdot n](2.45)$$

На рис. 7 представлены графики начальных данных (2.43), на рис. 8 - модули их коэффициентов разложения $a_k^{\left(m\right)}$.

Рис. 7

Рис. 8

В начальный момент времени все гармоники разложения (2.44) подогнаны друг к другу так, что в сумме составляют начальные данные. В последующие моменты времени каждая гармоника будет уменьшаться в амплитуде в соответствии с ее местом на диссипативной поверхности, и распространяться со своей фазовой скоростью, определяемой дисперсионной поверхностью, так что «пакет» гармоник начнет разъезжаться. Полная картина будет еще зависеть от числа Куранта-Фридрихса-Леви $r$.

Из рис. 6 видно, что наибольшая диссипация высоких гармоник имеет место при $r=0.5$. На рис. 9 приведен результаты расчета переноса начальных профилей (2.43), рассчитанные по схеме «уголок» при разных числах Куранта, до одного и того же момента времени. Сплошной линией показано аналитическое решение дифференциальной задачи, линией с маркерами - решение, рассчитанное по разностной схеме.

Рис. 9. Результаты расчета по схеме «уголок» (Simplest Upwind)

← Назад Далее →