Metody numeryczne w symulacji sieci gazowych - ebook

PRZEDMOWA

Książka jest przeznaczona dla studentów, doktorantów i pracowników naukowych zajmujących się modelowaniem matematycznym i symulacją komputerową sieci gazowych, a także ogólnie sieci płynowych. Książka powstała na podstawie moich notatek (rozszerzonych o badania własne) do wykładów _Metody numeryczne_ oraz _Komputerowa symulacja sieci gazowych_ prowadzonych dla studentów Specjalności – Systemy Ciepłownicze i Gazownicze na Wydziale Inżynierii Środowiska w Politechnice Warszawskiej. Do pełnego zrozumienia materiału zawartego w książce czytelnik powinien posiadać wiedzę z matematyki przynajmniej w zakresie podanym w pracach oraz .

Książka składa się z 11 rozdziałów. W rozdziale 1 (Wiadomości wstępne) Czytelnik zostaje wprowadzony w problematykę symulacji sieci gazowych. Algorytmy i metody numeryczne omawiane w następnych rozdziałach są bowiem niezbędnymi narzędziami umożliwiającymi konstruowanie efektywnych algorytmów komputerowej symulacji sieci gazowych. Rozdział drugi zawiera podstawowe informacje dotyczące rachunku macierzowego pomocne w zrozumieniu treści pozostałych rozdziałów. Rozdział trzeci dotyczy wybranych elementów teorii grafów. W przypadku komputerowej symulacji sieci gazowych o dowolnych strukturach istotnym jest taki zapis struktury sieci, aby był jednoznaczny, umożliwiał łatwą modyfikację i pozwalał na efektywne prowadzenie obliczeń. Te kryteria spełnia graf. Rozdział czwarty to podstawowe metody rozwiązywania układów równań algebraicznych liniowych. Omówione zostały zarówno metody bezpośrednie, jak również iteracyjne. W rozdziale piątym scharakteryzowane zostały metody rozwiązywania układów równań liniowych, ale z rozrzedzonymi macierzami. Rozrzedzona struktura macierzy współczynników wynika ze struktury symulowanej sieci. Rozdział szósty to metody rozwiązywania równań i układów równań nieliniowych. Rozdziały czwarty, piąty i szósty to podstawowa wiedza wykorzystywana w algorytmach symulacji stanów ustalonych sieci gazowych. Rozdział siódmy i rozdział ósmy to odpowiednio interpolacja i aproksymacja funkcji. Wiedza niezbędna przy przygotowywaniu danych do symulacji dynamicznej, jak również do dokładnej analizy rezultatów symulacji. W rozdziale dziewiątym przedstawiono podstawowe algorytmy całkowania numerycznego – bardzo istotny element algorytmów symulacji dynamicznej. Rozdziały dziesiąty i jedenasty zawierają podstawową wiedzę dotyczącą numerycznych metod rozwiązywania równań różniczkowych zwyczajnych i cząstkowych konieczną do budowania efektywnych algorytmów symulacji dynamicznej.

Dziękuję Panu prof. dr Yuriyovi Povstenko oraz Panu prof. dr hab. inż. Waldemarowi Kamratowi recenzentom mojej książki. Ich bardzo wnikliwe recenzje pozwoliły uniknąć wielu pomyłek.

Andrzej J. Osiadacz

Warszawa, 25.01.2025 r.ROZDZIAŁ 1.
WIADOMOŚCI WSTĘPNE

Prawidłowa eksploatacja istniejących systemów gazowniczych, jak również ich rozbudowa nie może być realizowana bez pomocy programów obliczeniowych. Prawidłowa eksploatacja sieci to przede wszystkim realizacja dostaw gazu odbiorcom zgodnie z podpisanymi umowami, przy jednoczesnej minimalizacji szeroko rozumianych kosztów eksploatacji. Z kolei prawidłowa rozbudowa sieci to wybór wariantu, który spełni warunki przyjętego kryterium kosztów eksploatacji oraz budowy.

Oprogramowanie wykorzystujące modele matematyczne elementów sieci powinno zatem stanowić podstawowe narzędzie pracy służb odpowiedzialnych za transport gazu.

Model matematyczny jest układem równań, który w określony sposób i z określoną dokładnością reprezentuje rzeczywisty system. Model matematyczny jest oczywiście uproszczeniem systemu rzeczywistego, ponieważ zawiera tylko te aspekty systemu, które odpowiadają wymaganiom stawianym modelowi.

Modele można klasyfikować zgodnie z założeniami przyjętymi przy ich tworzeniu. Jeżeli system ma zmienne, które zmieniają się w sposób ciągły w czasie, to należy używać modeli dynamicznych opisywanych równaniami ogólnie różniczkowymi. Jeżeli zmienne nie są funkcją czasu, wtedy należy stosować model statyczny opisywany równaniami algebraicznymi. Wybór modelu zależy od konkretnej sytuacji. Na przykład, w sieci gazowej niskociśnieniowej zmiany ciśnienia i przepływu są bardzo szybkie. Są one pomijane w większości rozpatrywanych przypadków, stąd stosowanie modeli stanu ustalonego ‒ równań algebraicznych. W przypadku sieci wysokiego ciśnienia dynamika przepływającego gazu jest znacznie wolniejsza ze względu na ogromne ilości gazu zmagazynowane w rurociągach. Pominięcie jej prowadziłoby do dużych błędów w opisie zjawiska. Stąd konieczność stosowania modeli matematycznych w postaci równań różniczkowych.

Głównym celem stosowania programów symulacyjnych jest badanie zachowania się sieci gazowych w określonych warunkach. Rozwiązanie określonego układu równań algebraicznych lub różniczkowych opisujących sieć dla założonych parametrów obciążenia i zasilania sieci da nam odpowiedź na pytania dotyczące wartości ciśnienia lub przepływu w wybranych punktach lub fragmentach sieci. Dzięki algorytmom symulacji możemy znacznie zmniejszyć nakłady na wyposażenie pomiarowe systemów gazowniczych, a także, zakładając poprawność modeli, uzyskiwać pełną i dokładną informację o systemie. Efektywne zarządzanie takim systemem wymaga zastosowania zaawansowanych metod numerycznych, algorytmów obliczeniowych o dużym stopniu złożoności i odpowiedniej jakości środków technicznych.

W zależności od wykorzystywanych modeli matematycznych opisujących przepływ gazu w gazociągu symulację dzielimy na statyczną i dynamiczną.

Symulacja statyczna sieci gazowej

Algorytmy symulacji statycznej opracowywane są w oparciu o:

• I prawo Kirchhoffa

gdzie:

‒ wektor obciążenia w węzłach,

‒ wektor przepływu w łukach grafu sieci,

‒ zredukowana macierz incydencji węzłów i łuków,

_n_ ‒ liczba węzłów,

_m_ ‒ liczba łuków,

_n_₁ ‒ liczba węzłów wyróżnionych (źródła);

• II prawo Kirchhoffa

gdzie:

‒ wektor spadków ciśnienia w łukach,

‒ macierz incydencji oczek podstawowych i łuków;

k – liczba niezależnych oczek;

• równanie przepływu, które jest nieliniowym równaniem algebraicznym _Q_ = _f_ (Δ_p_), a dla całej sieci jest postaci

Q = Ψ(ΔP)

gdzie:

Ψ(ΔP) ‒ wektor funkcji spadku ciśnienia w łukach,

‒ wektor spadku ciśnienia w łukach.

Danymi wejściowymi do programów symulacji statycznej są:

• struktura topologiczna sieci,

• wymiary geometryczne poszczególnych rurociągów,

• jakość wewnętrznej powierzchni rurociągów,

• wartości obciążenia oraz parametry zasilania w określonych węzłach sieci.

Rezultatem symulacji statycznej są wartości przepływu w łukach spełniające I prawo Kirchhoffa, wartości ciśnienia w węzłach zapewniające spadki ciśnienia w każdym łuku spełniające II prawo Kirchhoffa w każdym oczku, a relacje między spadkiem ciśnienia a przepływem w każdym łuku spełniają przyjęte równanie przepływu.

Symulacja dynamiczna sieci gazowej ,

Stany nieustalone w sieci występują w sytuacji jej niezbilansowania. Ilości gazu dostarczane do sieci w określonym przedziale czasu są różne od ilości gazu odbieranych z sieci. Mamy wtedy do czynienia ze stanem sczerpywania sieci – gdy ilości dostarczane są mniejsze od odbieranych, lub ze stanem ładowania w przypadku odwrotnym. Wartość ciśnienia w dowolnym punkcie sieci jest funkcją drogi _x_ (przebytej przez strumień gazu) oraz czasu _t_. Nieustalony przepływ gazu może być opisany układem równań różniczkowych zwyczajnych lub cząstkowych typu hiperbolicznego lub parabolicznego. Model matematyczny nieustalonego przepływu gazu wyprowadzany jest przy założeniu pewnych uproszczeń, z których najważniejszym jest założenie jednowymiarowości przepływu. Przyjęcie przepływu gazu jako jednowymiarowy jest zasadne, jeśli długość przewodu jest wielokrotnie dłuższa od jego średnicy, szybkość zmian takich parametrów jak ciśnienie, gęstość, temperatura czy prędkość w kierunku normalnym do kierunku linii prądu jest pomijalna w stosunku do szybkości zmian tych parametrów wzdłuż linii prądu. Oznacza to, że parametry przepływu gazu w dowolnym przekroju poprzecznym rurociągu traktowane są jako stałe, a przepływ jako jednorodny. Poza tym przyjmuje się, że promień krzywizny gazociągu jest duży w porównaniu z jego średnicą ‒ przekrój poprzeczny gazociągu jest stały (przedziałami stały), rurociąg jest sztywny (nie odkształca się pod wpływem wewnętrznego ciśnienia panującego w rurociągu) oraz kształt profili prędkości i temperatury (w przypadku przepływu nieizotermicznego) są w przybliżeniu stałe wzdłuż gazociągu.

Model matematyczny nieustalonego przepływu gazu

Do tworzenia modelu matematycznego nieustalonego przepływu gazu w rurociągu stosuje się następujące równania :

• równanie ciągłości,

• równanie ruchu,

• równanie energii,

• równanie stanu.

Równanie ciągłości

Równanie ciągłości jest zasadą zachowania materii w odniesieniu do zjawiska przepływu płynu. Masa, która jest miarą ilości materii, nie może powstać ani zniknąć, stąd zmiana masy w objętości kontrolnej może być spowodowana tylko różnicą wpływu i wypływu do i z objętości. Równanie (1.4) jest prawem zachowania masy w postaci różniczkowej

gdzie: ρ ‒ gęstość gazu (kg/m³), _w_ ‒ prędkość strumienia gazu (m/s²).

Mnożąc obie strony równania przez pole przekroju poprzecznego _A_ oraz uwzględniając, że _M =_ ρ ∙ _w_ ∙ _A_, gdzie _M_ ‒ przepływ masowy, otrzymamy

Zakładając izotermiczny proces przepływu gazu, równanie stanu możemy zapisać

ρ _c_² = _p_

gdzie:

_c_ ‒ izotermiczna prędkość dźwięku w gazie,

_p_ ‒ ciśnienie gazu.

Podstawiając równanie (1.4) do (1.6), otrzymujemy

Równanie ruchu

Równanie ruchu (II prawo dynamiki Newtona) mówi, że algebraiczna suma sił działających na element płynu o określonej masie w określonej chwili jest równa zmianie pędu tego elementu w tej chwili

Lewa strona równania (1.8) może być wyrażona następująco:

W przypadku transportu gazu rurociągiem wyróżnia się trzy składowe siły działającej na element gazu

gdzie:

_F_₁ ‒ siła od ciężaru gazu znajdującego się wewnątrz powierzchni kontrolnej,

_F_₂ ‒ siła od oporu hydraulicznego,

_F_₃ ‒ siła od ciśnienia gazu,

α – kąt nachylenia osi rurociągu do poziomu,

λ – współczynnik oporu hydraulicznego Darcy’ego-Weisbacha,

_D_ ‒ średnica rurociągu.

Po przekształceniach otrzymujemy

Równanie (1.11) jest ogólną postacią równania ruchu Newtona dla jednowymiarowego izotermicznego przepływu gazu.

W przedstawionych powyżej równaniach ciągłości oraz ruchu występują cztery zmienne: ciśnienie, gęstość, przepływ oraz temperatura, które są funkcjami czasu i zmiennej przestrzennej. Określenie wartości tych zmiennych będzie możliwe po dodaniu dodatkowych dwóch równań. Tymi równaniami są równanie stanu oraz równanie energii. Należy przy tym zauważyć, że jeśli przyjmiemy przemianę izotermiczną i odpowiadające jej równanie stanu uwzględnimy w równaniu ciągłości, to równanie (1.7) wraz z równaniem (1.11) tworzą model matematyczny nieustalonego izotermicznego przepływu gazu w gazociągu.

Równanie stanu

Równanie stanu opisuje wzajemną zależność pomiędzy parametrami stanu gazu (ciśnieniem, gęstością, temperaturą)

gdzie:

_R_ ‒ indywidualna stała gazowa, wyrażona w m²s‒2K‒1,

_Z_ ‒ współczynnik ściśliwości, zależny od wartości krytycznych temperatury i ciśnienia oraz ich wartości rzeczywistych,

_T_ – temperatura gazu.

W przemianie izotermicznej (_T_ = const) zależność pomiędzy ciśnieniem a gęstością przyjmuje postać

Przekształcając powyższe równanie oraz korzystając z definicji izotermicznej prędkości dźwięku w gazie, otrzymujemy

Równanie energii

Pierwsze prawo termodynamiki w odniesieniu do układu znajdującego się wewnątrz objętości kontrolnej jest postaci

gdzie:

_Q_ ‒ ilość ciepła dostarczonego do układu,

_W_ – praca wykonana przez układ,

Δ_E_ – zmiana energii układu.

Energia związana z masą układu wyrażana jest następująco:

gdzie:

_U_ ‒ energia wewnętrzna związana z oddziaływaniami molekularnymi w układzie,

(1/2)_mw_² ‒ energia kinetyczna układu,

_mgz_ ‒ energia potencjalna związana z położeniem układu w polu grawitacyjnym.

Po przekształceniach i pewnych uproszczeniach wynikających z definicji warunków, w jakich odbywa się przepływ, równanie energii w postaci różniczkowej przyjmie postać

gdzie:

_u_ ‒ energia wewnętrzna właściwa, tzn. wartość na jednostkę masy gazu;

‒ człon określający ilość ciepła wywiązującą się w płynie wewnątrz powierzchni kontrolnej w jednostce czasu;

‒ człon dotyczący lokalnej zmiany energii odpowiednio – wewnętrznej, kinetycznej i potencjalnej w czasie;

‒ człon opisujący strumień energii przez powierzchnię kontrolną;

‒ człon dotyczący zmiany strumienia energii wskutek konwekcji;

‒ człon reprezentujący strumień energii związany z pracą naprężeń normalnych.

W przypadku przepływu adiabatycznego oraz izotermicznego człon wyrażający ilość ciepła wywiązującego się w płynie na jednostkę czasu δ_Q /dt_ przyjmuje odpowiednio wartości:

• δ_Q /dt_ = 0 ‒ przepływ adiabatyczny,

• δ_Q /dt_ ≠ 0, _T_ = const ‒ przepływ izotermiczny.

Przepływ adiabatyczny dotyczy procesów szybkich, w których zjawiska związane z wymianą ciepła są pomijalne. Przepływy izotermiczne (temperatura gazu wzdłuż gazociągu nie ulega zmianie, _T_ = const) dotyczą wolnych zmian przepływu, w których zjawiska związane z wymianą ciepła pomiędzy gazem a otaczającym gazociąg gruntem nie są uwzględnione. Równaniami opisującymi nieustalony, izotermiczny przepływ gazu w gazociągu jest równanie ciągłości (1.7) oraz równanie ruchu (1.11). W przypadku przepływu nieizotermicznego (_T_ ≠ const) następuje wymiana ciepła między gazem a otoczeniem. W takim przypadku należy uzupełnić układ równań (1.5) i (1.11) o równanie energii (1.17) oraz równanie stanu (1.12).

Modele nieizotermiczne w odróżnieniu od izotermicznych pozwalają uwzględniać zmiany temperatury gazu w sieci. Stabilizacji temperatury w strumieniu służy dostatecznie duża długość gazociągu, przy niezmiennych warunkach termicznych panujących w gruncie. Jednak dla małych odległości od stacji przetłoczonych, gdzie następuje wzrost temperatury sprężanego gazu, zmiany temperatury gazu _T_(_x_) są wyraźne.

Przedstawione powyżej modele matematyczne nieustalonego nieizotermicznego i izotermicznego przepływu gazu wymagają dużych nakładów obliczeniowych. Należy pamiętać, że nie w każdym przypadku istnieje potrzeba korzystania z pełnego modelu, aby zapewnić oczekiwaną dokładność obliczeń. W praktyce szukamy modeli, które są kompromisem między dokładnością opisu zjawiska symulowanego a niezbędnym nakładem obliczeniowym wymaganym do jego rozwiązania.

Takie modele otrzymuje się poprzez odrzucenie niektórych członów z równań dokładnych, w oparciu o ocenę ich ilościowego wpływu na wynik obliczeń. Oznacza to, że wykorzystywany w procesie symulacji model nieustalonego przepływu gazu musi być odpowiedni do jasno zdefiniowanych granicznych wartości wybranych parametrów eksploatacji gazociągu. Analiza tych parametrów jest warunkiem koniecznym odpowiedniego wyboru modelu.

Przykładowe uproszczone modele nieustalonego przepływu gazu ,

Dokładny model matematyczny nieustalonego, izotermicznego przepływu gazu w gazociągu poziomym ma postać

Zakładając, że mamy do czynienia z gazociągiem płaskim, a zmiany prędkości przepływającego gazu są powolne, równanie ruchu uprościć można do postaci

W przypadku szybkich zmian prędkości przepływu gazu, spowodowanych na przykład dużymi zmianami obciążenia gazociągu lub nagłymi zmianami parametrów dostawy gazu, w równaniu ruchu koniecznym staje się uwzględnienie członu inercyjnego. Wówczas

Uzupełniając równanie (1.19) o równanie ciągłości w procesie izotermicznym (1.7), otrzymamy model nieustalonego izotermicznego przepływu gazu opisany układem dwóch równań różniczkowych cząstkowych typu parabolicznego

Jeśli z kolei przyjmiemy równanie ruchu (1.20), to otrzymamy

Układ (1.22) jest układem dwóch równań różniczkowych cząstkowych typu hiperbolicznego.

Dokładny model matematyczny nieustalonego, nieizotermicznego przepływu gazu w gazociągu poziomym ma postać

Przy założeniu, że wymiana ciepła może nastąpić jedynie jako przenikanie przez ściankę gazociągu lub jako efekt podłużnego przewodnictwa cieplnego (wszelkie możliwe reakcje chemiczne, dysocjację, jonizację, pochłanianie lub emisję promieniowania zaniedbujemy, przy czym gaz traktujemy jako ośrodek jednorodny o stałej przewodności cieplnej), z przeprowadzonej analizy wartości poszczególnych członów dokładnego równania energii dla określonych warunków eksploatacyjnych i geometrii gazociągu wynika, że równanie energii można uprościć do postaci

Oznacza to, że uproszczona postać modelu nieustalonego nieizotermicznego przepływu gazu w rurociągu poziomym może być opisana poniższym układem równań

gdzie:

‒ lokalna zmiana gęstości w czasie,

‒ strumień masy przypadający na jednostkę powierzchni,

‒ zmiana ciśnienia gazu na długości,

‒ siły bezwładności w układzie,

‒ zmiana energii kinetycznej (zmiana pola prędkości),

‒ składowa sił grawitacji,

‒ siła od oporu hydraulicznego,

‒ zmiana ilości ciepła w strumieniu gazu na długości,

‒ rozprężanie gazu na długości,

‒ zmiana energii potencjalnej na długości,

‒ zmiana ilości ciepła w jednostce objętości gazu w czasie,

‒ przenikanie ciepła przez ściankę gazociągu.

Bardzo istotnym zagadnieniem jest odpowiedź na pytanie, w jakich warunkach eksploatacji gazociągu/systemu przesyłowego należy stosować model hiperboliczny, a w jakich model paraboliczny. Wiadomo, że model hiperboliczny ze względu na swoją strukturę dokładniej opisuje szybko zachodzące zmiany parametrów gazu spowodowane nagle zmieniającymi się warunkami brzegowymi. Oznacza to, że istnieje pewna graniczna częstotliwość zmian warunków brzegowych, powyżej której zaleca się ze względu na dokładność opisu zjawiska stosowanie modelu hiperbolicznego. Należy jednocześnie pamiętać, że to, jaka będzie reakcja gazociągu na zakłócenie, zależy nie tylko od jakości zakłócenia, ale także od pojemności akumulacyjnej gazociągu, czy też od średniego ciśnienia w gazociągu. Sformułowanie jasnych kryteriów doboru modeli do warunków pracy gazociągu jest bardzo istotne, jako że model paraboliczny jest łatwiejszy do numerycznego rozwiązania, istnieje wiele efektywnych obliczeniowo numerycznych algorytmów, a także, co jest nie bez znaczenia, znacznie łatwiej sformułować algorytm symulacji sieci o dowolnej strukturze z elementami nierurowymi (m.in. tłocznie, sprężarki, zawory, reduktory, stacje gazowe) o różnych wariantach pracy tych elementów. Modele paraboliczne aczkolwiek są znacznie łatwiejsze do numerycznego rozwiązywania i wykorzystywania do tworzenia programów symulacyjnych pozwalających na obliczanie sieci o dowolnych strukturach i dowolnej liczbie elementów nierurowych, jednak w pewnych warunkach mniej dokładnie opisują zmiany zachodzące w sieci w stosunku do modeli hiperbolicznych. Te ostatnie z kolei są bardziej pracochłonne obliczeniowo, a także trudniej je wykorzystać do stworzenia efektywnego symulatora sieci o dowolnej strukturze.

Algorytmy symulacji nieustalonego przepływu gazu w sieci o dowolnej strukturze

Algorytm symulacji dynamicznej oblicza wartości ciśnienia lub przepływu w wybranych punktach (węzłach) sieci w określonych chwilach. W porównaniu z algorytmem stanów ustalonych liczba danych wejściowych niezbędnych do wykonania obliczeń jest znacznie większa, ponieważ zarówno parametry obciążenia (odbiór gazu z sieci), jak również parametry zasilania systemu (parametry źródeł) są funkcją czasu. Rośnie w sposób znaczący w porównaniu z symulacją statyczną czas obliczeń, ponieważ:

• problem polega na rozwiązaniu układu równań różniczkowych cząstkowych lub zwyczajnych zamiast algebraicznych,

• symulacja dynamiczna dotyczy określonego przedziału czasu, a więc wielokrotnego rozwiązywania układu równań w tym przedziale czasu.

Wnioski

Dzięki algorytmom symulacji możemy znacznie zmniejszyć nakłady na wyposażenie pomiarowe systemów gazowniczych, a także, zakładając poprawność modeli, uzyskiwać pełną i dokładną informację o systemie. Algorytmy symulacji powinny umożliwiać realizację następujących funkcji:

• Przeprowadzanie symulacji sieci w całości lub wybranych jej fragmentów w trybie off-line w celu analizy reakcji sieci na określone wymuszenia lub zmiany jej konfiguracji.

• Symulacja on-line na bazie danych dostarczonych przez system telemetryczny.

• Ocena jakości pracy systemu telemetrycznego i urządzeń pomiarowych.

• Symulacja w celu monitorowania pracy sieci.

• Śledzenie zmian akumulacji gazu w systemie.

• Wydawanie zapewnień dostawy gazu przyszłym odbiorcom.

• Sprawdzenie jakości funkcjonowania sieci podczas prowadzenia remontu wybranej jej części.

• Określenie masy gazu traconej w systemie na skutek nieszczelności.

• Ocena możliwości przepustowych sieci.

• Śledzenie procesu mieszania się gazów o różnych składach.

• Lokalizacja nieszczelności.

Literatura do rozdziału pierwszego

1. Kołodziej W., _Analiza matematyczna_, Państwowe Wydawnictwo Naukowe, Warszawa 2020.

2. Krysicki W., Włodarski L., _Analiza matematyczna w zadaniach, Część I, Część II,_ Państwowe Wydawnictwo Naukowe, Warszawa 1999.

3. Osiadacz A.J., _Simulation and Analysis of Gas Networks_, [email protected]. Spon, London 1987.

4. Osiadacz A.J., _Statyczna symulacja sieci gazowych_, Fluid Systems, Warszawa 2001.

5. Osiadacz A.J., _Symulacja i optymalizacja sieci gazowych_, _Vademecum Gazownika_, Tom IV, Rozdział 6, SITPNIG, Kraków 2012.

6. Osiadacz A.J., Chaczykowski M._, Comparison of isothermal and non-isothermal pipeline flow models_, Chemical Engineering Journal, vol. 81, no 1‒3, 2001, s. 41‒51.ROZDZIAŁ 2.
PODSTAWY RACHUNKU MACIERZOWEGO

2.1. MACIERZE ORAZ DZIAŁANIA NA MACIERZACH

_Macierzą_ A (dokładniej macierzą dwuskładnikową) nazywamy funkcję dwóch zmiennych, która każdej parze liczb naturalnych (_i_,_j_) (_i_ = 1, 2, ..., _m_; _j_ = 1, 2, ..., _n_) przyporządkowuje dokładnie jeden element _a_i,j będący liczbą rzeczywistą lub zespoloną, lub operatorem (np. różniczkowania), bądź wielomianem. Macierz zapisujemy w postaci tablicy

lub krócej jako A = , _i_ = 1, 2, ..., _m_; _j_ = 1, 2, ..., _n_; Am×n. Element _a_i,j znajduje się w _i_-tym wierszu oraz w _j_-tej kolumnie macierzy.

Macierz prostokątną, której liczba kolumn jest równa jeden (_n_ = 1), nazywamy _macierzą kolumnową_, __ którą można traktować jako _m_-wymiarowy wektor kolumnowy. Macierz prostokątną, której liczba wierszy wynosi jeden (_m_ = 1), nazywamy _macierzą wierszową_, którą można traktować jako _n_-wymiarowy wektor wierszowy. Traktując kolumny macierzy jako _m_-wymiarowe wektory

( _j_ = 1, 2, ..., _n_)

można macierz (2.1) zapisać jako

Jeśli wiersze macierzy potraktujemy jako _n_-wymiarowe wektory wierszowe postaci Ai = (_i_ = 1, 2, …, _m_), to macierz (2.1) zapiszemy jako

_Wymiarem macierzy_ (2.1) nazywamy parę uporządkowaną zawierającą liczbę wierszy _m_ i kolumn _n_, którą oznaczamy przez _m_ × _n_.

Macierz, której liczba wierszy _m_ jest równa liczbie kolumn _n_, nazywamy macierzą _kwadratową._ W przeciwnym wypadku mamy do czynienia z macierzą _prostokątną._

Macierz kwadratową, w której różne od zera elementy, _d_ii ≠ 0 (_i_ = 1, 2, …, _n_), występują tylko na głównej przekątnej, nazywamy _macierzą diagonalną_ i oznaczamy symbolem diag(_d_₁₁, _d_₂₂, …, _d_nn)

Jeżeli w macierzy diagonalnej wszystkie _d_ii ≠ 0 są sobie równe, to taką macierz nazywamy macierzą _skalarną_. Szczególnym przypadkiem macierzy D jest _macierz jednostkowa_ In stopnia _n_, określona wzorem

Wprowadzając symbol Kroneckera

macierz jednostkową zapiszemy w postaci I = .

Macierz, której wszystkie elementy są równe zero, nazywamy _macierzą zerową._

Macierz kwadratowa, której elementy spełniają warunek _a_ij = _a_ji dla _i_ = 1, 2, ..., _m_; _j_ = 1, 2, ..., _n_, nazywamy _macierzą symetryczną._

Macierz kwadratową nazywamy trójkątną dolną ‒ L oraz górną ‒ U, jeśli wszystkie elementy położone poniżej/powyżej głównej przekątnej są równe zero.

Macierz trójkątna dolna stopnia _n_ jest postaci

a macierz trójkątna górna

PODSTAWOWE DZIAŁANIA NA MACIERZACH

Dwie macierze A = mn i B = mn są równe wtedy i tylko wtedy, gdy ich odpowiednie elementy są sobie równe, tzn. A = B, gdy = dla _i_ = 1, 2, …, _m_; _j_ = 1, 2, …, _n_.

_Sumą_ (_różnicą_) _dwóch macierzy_ A = mn i B = mn o jednakowych wymiarach jest macierz C o elementach równych sumom (różnicom) odpowiednich elementów macierzy A i B, tzn. C = A ± B, gdy _c_ij = ± dla _i_ = 1, 2, …, _m_; _j_ = 1, 2, …, _n._

Z powyższej definicji wynikają następujące własności:

(A + B) + C = A + (B + C )

A + B = B + A

A + 0mn = A

_Iloczyn_ A × _k_ lub _k_ × A macierzy A przez liczbę _k_ jest macierzą o elementach równych odpowiednim elementom macierzy A pomnożonym przez _k._

Niech będą dane dwie macierze

w których liczba kolumn macierzy A jest równa liczbie wierszy macierzy B.

_Iloczynem macierzy_ A i B nazywamy macierz , której elementy określone są zależnością dla _i_ = 1, 2, …, _m_; _k_ = 1, 2, …, _q_.

Na przykład iloczynem macierzy przez jest macierz .

Mnożenie macierzy na ogół nie jest przemienne, czyli AB ≠ BA.

W związku z nieprzemiennością mnożenia macierzy wyróżnia się mnożenie prawostronne macierzy A przez macierz B, czyli C = AB, oraz mnożenie lewostronne, czyli D = BA.

Macierze A i B spełniające warunek AB = BA nazywamy _macierzami przemiennymi_.

Ogólne własności mnożenia macierzy są następujące:

Iloczyny macierzy trójkątnych dolnych są także macierzami trójkątnymi dolnymi, a iloczyny macierzy trójkątnych górnych są macierzami trójkątnymi górnymi.

Transpozycja macierzy

_Macierzą transponowaną_ AT macierzy A = o wymiarach _m_ × _n_ nazywamy macierz o wymiarach _m_ × _n_, którą otrzymano w rezultacie zamiany wierszy z kolumnami. Elementy macierzy A zastąpiono elementami macierzy AT.

Dla dowolnej macierzy prostokątnej

macierz transponowana ma postać:

Wektor wierszowy xT = (_x_₁, _x_₂, …, _x_n) jest transpozycją pewnego wektora kolumnowego. Można wykazać, że:

Jeżeli A jest macierzą kwadratową, dla której A = AT, to macierz A jest macierzą symetryczną, tzn. _a_ij = _a_ji dla wszystkich _i_ oraz wszystkich _j_. W macierzy symetrycznej elementy symetryczne względem głównej przekątnej są równe.

Iloczyn dwóch macierzy symetrycznych A i B jest symetryczny tylko pod warunkiem, że AB = BA.

Macierz symetryczna nazywa się _dodatnio określoną_, jeśli związana z nią forma kwadratowa xTAx spełnia warunek xTAx > 0 dla każdego rzeczywistego _x_ ≠ 0.

_Macierzą ortogonalną_ nazywamy macierz kwadratową An×n spełniającą równość AAT = AAT= In.

Macierze blokowe

Dowolną macierz A można przedstawić jako _macierz blokową_, zbudowaną z pewnej liczby macierzy o mniejszych wymiarach

gdzie Aij jest macierzą o wymiarach _m_i × _n_i (_i_ = 1, 2, ..., _q_ ≤ _m_; _j_ = 1, 2, ..., _p_ ≤ _n_).

W przypadku kiedy macierze na przekątnej Aii są kwadratowe, macierz A również musi być kwadratowa, a _p_i = _q_i (_i_ = 1, 2, …, _n_). Dodawanie i mnożenie takich macierzy blokowych wykonuje się, traktując, że bloki są liczbami.

Na przykład, dla C = AB istnieje zależność

_Macierzą blokowo-przekqtniową_ nazywamy macierz, którą można zapisać w postaci

A = diag(A₁₁, A₂₂, …, Ann)

przy czym macierze Aii muszą być kwadratowe.

Analogicznie definiuje się _macierz blokowo-trojkątną_.

Wyznaczniki

_Wyznacznik macierzy kwadratowej_ A stopnia _n_ oznacza się symbolem det(A) lub | A|:

W przypadku macierzy prostokątnej A = m×n skreślamy pewna liczbę wierszy i kolumn, tak aby z elementów nieskreślonych macierzy prostokątnej utworzyć macierz kwadratową. Wyznacznik tak otrzymanej macierzy nazywamy _minorem macierzy_ A.

Najwyższy spośród stopni minorów macierzy A różnych od zera nazywamy _rzędem macierzy_ A i oznaczamy _r_(A). Macierz stopnia _m_ × _n_ ma rząd spełniający nierówność 0 ≤ _r_(A) ≤ minor(_m_,_n_).

Rząd macierzy jest równy maksymalnej liczbie liniowo niezależnych wierszy lub kolumn. Rząd kolumnowy macierzy A to liczba liniowo niezależnych kolumn macierzy A_._ Rząd wierszowy macierzy A to liczba liniowo niezależnych wierszy macierzy A_._

Zawsze rząd kolumnowy i wierszowy macierzy A są sobie równe i będziemy je nazywać rzędem macierzy A (oznaczenie _r_(A)).

W szczególnym przypadku, jeśli _r = m = n_, macierz A jest _macierzą nieosobliwą_.

PRZYKŁAD

Określmy rząd macierzy

Tworzymy minory stopnia 3

Minor stopnia 2 wynosi

Oznacza to, że rząd macierzy wynosi 2.

Minorem _M_ij elementu _a_ij wyznacznika det A stopnia _n_ nazywamy wyznacznik stopnia _n_ ‒ 1 powstały przez usunięcie z wyznacznika det A _i_-tego wiersza oraz _j_-tej kolumny.

Wartość liczbową wyznacznika detA stopnia _n_ można obliczyć, stosując poniższe wzory:

det A = _a_i1_A_i1 + _a_i2_A_i2 _+_ … + _a_in_A_in (1 ≤ _i_ ≤ _n_)

det A = _a_1j _A_1j + _a_2j _A_2j _+_ … + _a_nj _A_nj (1 ≤ _j_ ≤ _n_)

gdzie A1k (_k =_ 1, 2, ..., _n_) oznacza wyznacznik stopnia _n_ ‒ l, który powstaje przez skreślenie z macierzy A pierwszego wiersza i _k-_tej kolumny.

Własności wyznaczników:

• wyznacznik macierzy trójkątnej jest równy iloczynowi elementów głównej przekątnej det(L) = _l_₁₁ · _l_₂₂, …, · _l_nn, det(U) = _u_₁₁· _u_₂₂, ..., · _u_nn,

• transponowanie macierzy kwadratowej nie zmienia wyznacznika tej macierzy detAT = detA,

• zamiana dwóch wierszy (kolumn) macierzy kwadratowej zmienia wartość wyznacznika tej macierzy na przeciwną,

• jeśli w macierzy dwa wiersze (dwie kolumny) są identyczne , to wyznacznik tej macierzy wynosi zero,

• jeżeli macierz B powstaje z macierzy A przez pomnożenie wszystkich elementów pewnego wiersza (kolumny) przez liczbę α, to detB = αdetA,

• dodanie do wiersza (kolumny) wielokrotności innego wiersza (kolumny) nie zmienia wyznacznika tej macierzy,

• wyznacznik iloczynu dwu macierzy kwadratowych jest równy iloczynowi wyznaczników tych macierzy det(A · B) = detA detB.

OBLICZANIE WYZNACZNIKÓW

1. Macierz drugiego stopnia

2. Macierz trzeciego stopnia:

• wzór Laplace’a

• schemat Sarrusa

• reguła Chio

przy założeniu, że _a_₁₁ ≠ 0

ODWRACANIE MACIERZY

_Macierzą odwrotną_ A–1 macierzy kwadratowej stopnia _n_ nazywamy macierz spełniającą następujący warunek: AA‒1 = A‒1A = I. W przypadku iloczynu macierzy obowiązuje zależność (AB)–1 = A–1B–1.

Macierz odwrotna A–1 istnieje wtedy i tylko wtedy, gdy A jest macierzą _nieosobliwą_, tzn. jeśli det(A) ≠ 0. W przeciwnym wypadku macierz określa się jako _macierz osobliwą_.

Macierz odwrotną obliczamy zgodnie z zależnością

gdzie AD – _macierz dołączona_.

Dla każdego elementu _a_ij (_i_, _j_ = 1, 2, …, _n_) macierzy A możemy obliczyć dopełnienie algebraiczne _A_ij = (–1)i+j_M_ij, gdzie _M_ij jest minorem, który powstał przez skreślenie _i_-tego wiersza oraz _j_-tej kolumny w macierzy A.

Aby wyznaczyć macierz dołączoną AD macierzy A, należy:

• wyznaczyć macierz dopełnień algebraicznych dla danej macierzy kwadratowej A, tzn. zastąpić każdy element _a_ij macierzy A odpowiadającym mu dopełnieniem algebraicznym,

• wyznaczyć macierz transponowaną macierzy dopełnień algebraicznych.

PRZYKŁADY

Aby wyznaczyć macierz odwrotną macierzy , obliczamy: detA = (3 · 1 · 2 + 2 · 1 · 1 + 2 · 0 · 0) – (1 · 1 · 0 + 0 · 1 · 3 + 2 · 2 · 2) = 0, detA = 0, co oznacza, że macierz A–1 nie istnieje.

Aby wyznaczyć macierz odwrotną macierzy , sprawdzamy, że detB = 8, co oznacza, że B–1 istnieje. Obliczamy macierz dopełnień algebraicznych dla macierzy B.

Ostatecznie

Otrzymane rozwiązanie jest poprawne ponieważ BB–1 = I.

PRZESTRZEŃ LINIOWA WEKTOROWA

_Wektorem_ x w przestrzeni _n_-wymiarowej nazywamy uporządkowany zbiór składający się z _n_ liczb rzeczywistych lub zespolonych: . Liczby _x_₁, _x_₂, …, _x_n nazywamy współrzędnymi wektora.

Przykłady wektorów:

• wektor na płaszczyźnie jest wektorem przestrzeni dwuwymiarowej,

• w macierzy A o wymiarze _m_ × _n_ każdy jej wiersz jest wektorem _n_-wymiarowym, a każda kolumna wektorem _m_-wymiarowym.

Dwa wektory x oraz y w przestrzeni _m_-wymiarowej, , są równe, jeśli _x_i = _y_i dla _i_ = 1, 2, …, _m_.

Sumą/różnicą dwóch wektorów nazywamy wektor przestrzeni _m_-wymiarowej o współrzędnych odpowiednio _x_i + _y_i lub _x_i – _y_i (_i_ = 1, 2, …, _m_).

W przypadku dodawania wektorów słuszne są następujące równości:

x + y = y + x __ (przemienność)

(x + y) + z = x + ( y + z) (łączność)

Iloczynem wektora przez liczbę α nazywamy wektor .

Dowolny zbiór wektorów _n_-wymiarowych, w którym są określone operacje dodawania wektorów oraz mnożenia wektora przez liczbę, przy czym rezultaty tych działań należą do tego zbioru, nazywamy _liniową przestrzenią wektorową_ i oznaczamy przez En.

Wektor y = _c_₁x₁ + _c_₂x₂ + … + _c_kxk nazywany jest _kombinacją liniową wektorów_ x₁, x₂, ..., xk. Wektory te są nazywane _liniowo niezależnymi_, jeśli równość _c_₁x₁ + _c_₂x₂ + … + _c_kxk = 0 zachodzi tylko pod warunkiem, że _c_₁ = _c_₂ = ... = _c_k _=_ 0. W przeciwnym wypadku wektory te określa się jako _liniowo zależne_.

_Wymiarem liniowej przestrzeni wektorowej_ nazywamy maksymalną liczbę wektorów liniowo niezależnych. Dowolny zbiór _n_ wektorów liniowo niezależnych _n_-wymiarowej przestrzeni tworzy _bazę_ tej przestrzeni.

Każdy wektor x ∈ _R_n można więc wyrazić jako

x = _a_₁y₁ + _a_₂y₂ + … + _a_nyn

gdzie _a_₁, _a_₂, …, _a_n nazywa się _współrzędnymi_ tego wektora względem bazy y₁, y₂, ..., yn.

W przestrzeni liniowej rzeczywistej jest określony iloczyn skalarny, jeśli każdej parze wektorów x oraz y należących do tej przestrzeni przyporządkowana jest liczba rzeczywista oznaczona przez (x, y). Liniowa przestrzeń rzeczywista, w której zdefiniowany jest iloczyn skalarny, nazywa się _przestrzenią euklidesową_.

W przestrzeni euklidesowej aksjomaty iloczynu skalarnego zdefiniowane są następująco:

(x, x) > 0, jeśli x ≠ 0,

(x, y) = ( y, x)

(α x, y) = (x, α y) = α (x, y) α – liczba rzeczywista

(x + y, z) = (x, z) + ( y, z), (x, y + z) = (x, y) + (x, z)

Iloczyn skalarny wektorów x oraz y jest liczbą określoną wzorem

a w zapisie macierzowym

_Długością_ (_modułem_) wektora w przestrzeni euklidesowej nazywamy liczbę

2.2. Wartości własne i wektory własne

Jeśli dla danej liczby λ i wektora x ≠ 0 spełniona jest równość

Ax = λx __

to __ λ nazywa się _wartością własną_ macierzy A, a x ‒ wektorem własnym, odpowiadającym wartości λ. Zadanie polega na wyznaczeniu takich wartości parametru λ, przy których powyższe równanie ma rozwiązanie niezerowe.

Równanie (2.13) można zapisać w postaci układu jednorodnych równań algebraicznych liniowych

lub w postaci macierzowo-wektorowej

gdzie A ‒ λI jest _macierzą charakterystyczną_ dla macierzy A = n×n.

Wyznacznikiem charakterystycznym macierzy A nazywamy wyznacznik

Równanie (2.16) nazywane jest _równaniem charakterystycznym_ macierzy A, natomiast pierwiastki równania charakterystycznego nazywane są _wartościami własnymi_ macierzy A. Zbiór wartości własnych (λ₁, λ₂, …, λn) __ macierzy kwadratowej A nazywamy _widmem_ tej macierzy. Wielomian charakterystyczny ma postać

Współczynniki _k_i (_i_ = 1, 2, …, _n_) są funkcją elementów macierzy A:

– suma minorów głównych stopnia drugiego macierzy A,

_k_i – suma minorów głównych stopnia trzeciego itd., _k_n = det A.

W przypadku przestrzeni _R_₂ macierz A ‒ λI przyjmuje postać

czyli:

det(A ‒ λI) = (_a_₁₁ ‒ λ)(_a_₂₂ – λ) – _a_₁₂_a_₂₁ = λ² – (_a_₁₁ + _a_₂₂)λ + (_a_₁₁_a_₂₂ – _a_₁₂_a_₂₁) = 0

tzn. w tym przypadku λ jest pierwiastkiem wielomianu drugiego stopnia.

PRZYKŁAD

Niech

czyli

Wielomian charakterystyczny jest postaci

_D_(λ) = det (A ‒ λI)= (1 ‒ λ)(2 ‒ λ)(3 ‒ λ)

co oznacza, że λ₁ = 1, λ₂ = 2, λ₃ = 3.

Niech

x – wektor własny odpowiada λ₁,

y_ _– wektor własny odpowiada λ₂,

z – wektor własny odpowiada λ₃.

Aby znaleźć x, rozwiązujemy układ równań

czyli

Dla _x_₃ = 0, przyjmując _x_₁ = 1, otrzymujemy _x_₂ = ‒2, tzn.

Wartości wektora y obliczamy z układu równań

czyli

Ponieważ _y_₃ = 0, _y_₁ = 0, stąd _y_₂ =1, tzn.

Po rozwiązaniu poniższego układu równań otrzymujemy

TWIERDZENIE CAYLEYA-HAMILTONA

Ważnym twierdzeniem w teorii macierzy jest twierdzenie Cayleya-Hamiltona, które mówi, że każda kwadratowa macierz A spełnia swoje równanie charakterystyczne, tzn. |A –λI| = 0. Jeżeli _D_(λ) jest wielomianem charakterystycznym macierzy A = |_a_|n×n, to wówczas _D_(A) = 0.

Jeśli wielomian charakterystyczny macierzy A jest postaci

to wówczas jest spełnione równanie macierzowe

PRZYKŁAD

Niech

Wtedy

Twierdzenie Cayleye’a-Hamiltona można wykorzystać do znalezienia macierzy odwrotnej. Mnożąc obustronnie równanie (2.17) przez A–1, otrzymujemy

Ostatecznie

PRZYKŁAD

Aby znaleźć macierz odwrotną do macierzy

obliczamy

Ostatecznie