Next: Feladat: upwind módszer kontinuitási Up: Gáz és plazma dinamika Previous: Feladat: FTCS módszer a Tartalomjegyzék

Stabilitás

Az elõzõ fejezetben felírt (6) diszkretizált egyenlet nem mûködik megfelelõen mint az a 2. feladat megoldásából kiderült. Ennek az az oka, hogy a kontinuitási egyenletre az FTCS módszer numerikusan instabil. Ezt könnyen láthatjuk, ha elvégezzük a Neumann féle stabilitási vizsgálatot. Bontsuk fel a kezdeti feltételt Fourier komponensekre, és vizsgáljuk meg, hogy az egyes komponensek amplitúdója hogyan változik az FTCS módszer alkalmazása során. Mivel a (6) egyenlet lineáris $\rho$ -ban, valamint a határfeltételeket periodikusnak választottuk, a stabilitás analízis egzakt eredményt fog adni. A hullámszámú Fourier komponens a

$\begin{displaymath} \rho^n_j = G_n \exp(i k x_j)~~~~~~~~~~j=1,\ldots,N \end{displaymath}$

(8)

alakban írható fel, ahol

az amplitúdó, és $i=\sqrt{-1}$ az imaginárius egység. Ha behelyettesítjük a (8) formulát az (6) egyenletbe, akkor a $G=G_{n+1}/G_n$ amplitúdó növekedési faktorra a következõ megoldást kapjuk:

$\begin{displaymath} \vert G\vert=\left\vert 1-v\frac{\Delta t}{2\Delta x} \lef... ...vert 1-v\frac{\Delta t}{\Delta x}i\sin(k\Delta x)\right\vert>1 \end{displaymath}$

(9)

Mint látható tetszõlegesen kis $\Delta t$ idõlépés mellett is $\vert G\vert>1$ , azaz a Fourier komponens amplitúdója minden lépésben nõ. Tehát az FTCS módszer a kontinuitási egyenletre nézve instabil.

Vizsgáljuk meg, hogyan javítható a diszkretizált egyenlet viselkedése. Az egyik megoldás az, ha a térbeli deriváltat féloldalasra változtatjuk:

$\begin{displaymath} \rho_j^{n+1}=\rho_j^n-v \Delta t \frac{\rho^n_{j}-\rho^n_{j-1}} {\Delta x} \end{displaymath}$

(10)

Természetesen nem mindegy, hogy a bal vagy a jobb oldali differenciát használjuk. Mivel itt

, azaz az áramlás balról jobbra halad, logikus, hogy a bal oldali $\rho^n_{j}-\rho^n_{j-1}$ differenciát használjuk. Ezt áramlásfelöli (angolul ,,upwind'') deriváltnak nevezzük.

Végezzük el a stabilitás vizsgálatot a (10) módszerre vonatkozóan. A (8) Fourier komponsnek az upwind módszer (10) képletébe való helyettesítésével az erõsítési faktorra a

$\begin{displaymath} G=1-v\frac{\Delta t}{\Delta x} \left(1-e^{-ik\Delta x}\right) = 1-C + C e^{-ik\Delta x} \end{displaymath}$

(11)

megoldást kapjuk, ahol

$\begin{displaymath} C=v\frac{\Delta t}{\Delta x} \end{displaymath}$

(12)

a Courant szám, ami azt adja meg, hogy $\Delta t$ idõ alatt az információ hány rácstávolságnyira terjed. A (11) egyenlet szerint az erõsítési faktor a komplex számsíkon egy

sugarú körön helyezkedik el, melynek középpontja a valós tengelyen

-nél van amint ez a 2. ábrán látható.

**Ábra 2.:** A G komplex erõsítési faktor az upwind módszerre. G a folytonos vonallal rajzolt körön helyezkedik el. A stabil tartományt a szaggatott kör jelöli.
$\begin{figure}\centerline{\epsfig{file=upwind-stab.eps,width=6cm}} \end{figure}$

A stabilitás feltétele

$\begin{displaymath} 0 \le C \le 1 \end{displaymath}$

(13)

azaz a Courant szám nem lehet 1-nél nagyobb. A

, azaz $\Delta x=v\Delta t$ választás rendkívül érdekes eredményt ad, ugyanis ekkor $G=\exp(ikx)$ , azaz a Fourier komponens értéke

$\begin{displaymath} \rho^n_j = (G)^n \rho^0_j = e^{i k x_j -i k n \Delta x} = e^{i k (x_j - v t_n)} \end{displaymath}$

(14)

ami éppen az egzakt megoldása a kontinuitási egyenletnek. Mivel minden Fourier komponensre egzakt a diszkrét megoldás, általában is igaz, hogy az upwind módszer a konstans sebességû kontinuitási egyenletre

Courant szám mellett egzakt megoldást ad. Ez másképpen is belátható:

esetén az upwind módszer minden idõlépésben éppen egy rácstávolságnyival viszi arrébb a megoldást. Valójában persze ez nem túl hasznos, hiszen a konstans sebességû kontinuitási egyenlet megoldásához nincs szükség számítógépes modellre. Reális, általánosabb egyenletre is jellemzõ eredményt a

Courant számmal végzett szimulációk mutatnak.

Az is látható, hogy ha elõjele negatív lenne, akkor $\Delta t$ -tõl függetlenül instabillá válik a módszer. Éppen ezért az általános upwind módszer elõjelének függvényében választja a bal vagy jobboldali differencia formulát:

$\begin{displaymath} \rho_j^{n+1}=\left\{ \begin{array}{ll} \rho_j^n-C (\rho^n_... ...ho^n_{j+1}-\rho^n_{j}) & \mbox{ha $v\le0$} \end{array}\right. \end{displaymath}$

(15)

Az így definiált upwind módszer feltételesen stabil $\vert C\vert<1$ -re, most már

elõjelétõl függetlenül.

Alfejezetek

Feladat: upwind módszer kontinuitási egyenletre

Next: Feladat: upwind módszer kontinuitási Up: Gáz és plazma dinamika Previous: Feladat: FTCS módszer a Tartalomjegyzék

Gabor Toth 2000-09-04