Projekt okładki i stron tytułowych Przemysław Spiechowski

Ilustracja na okładce fluidsystems.pl

Wydawca Adam Filutowski

Koordynator ds. redakcji Adam Kowalski

Redaktor Maria Kasperska

Produkcja Mariola Grzywacka

Dział reklamy Magdalena Bystrzycka (magdalena.lewocka@pwn.pl)
Urszula Obrycka (urszula.obrycka@pwn.pl)

Dział marketingu marketing.pwn@pwn.pl

Skład wersji elektronicznej na zlecenie Wydawnictwo Naukowe PWN S.A.:
Michał Latusek

Książka, którą nabyłeś, jest dziełem twórcy i wydawcy. Prosimy, abyś przestrzegał praw, jakie im przysługują. Jej zawartość możesz udostępnić nieodpłatnie osobom bliskim lub osobiście znanym. Ale nie publikuj jej w internecie. Jeśli cytujesz jej fragmenty, nie zmieniaj ich treści i koniecznie zaznacz, czyje to dzieło. A kopiując jej część, rób to jedynie na użytek osobisty. Szanujmy cudzą własność i prawo.

Szanujmy cudzą własność i prawo.

Więcej na www.legalnakultura.pl
Polska Izba Książki

Copyright © by Wydawnictwo Naukowe PWN SA

Warszawa 2025

ISBN 978-83-01-23966-4
https://doi.org/10.53271/2024.133

eBook został przygotowany na podstawie wydania papierowego z 2025 r. (Wydanie I)

Warszawa 2025

Wydawnictwo Naukowe PWN SA

02-460 Warszawa, ul. Gottlieba Daimlera 2

tel. 22 69 54 321; faks 22 69 54 288

infolinia 801 33 33 88

e-mail: pwn@pwn.com.pl; reklama@pwn.pl

www.pwn.pl

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 [1] oraz [2].

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 [4]

Algorytmy symulacji statycznej opracowywane są w oparciu o:

• I prawo Kirchhoffa

[1.1]

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

[1.2]

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)

[1.3]

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 [3],[5]

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 [3]:

• 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

[1.4]

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

[1.5]

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

ρ c² = p

[1.6]

gdzie:

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

p ‒ ciśnienie gazu.

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

[1.7]

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

[1.8]

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

[1.9]

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

[1.10]

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

[1.11]

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ą)

[1.12]

gdzie:

R ‒ indywidualna stała gazowa, wyrażona w m²s‒²K‒¹,

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ć

[1.13]

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

[1.14]

Równanie energii

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

[1.15]

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:

[1.16]

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ć

[1.17]

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 [3], [5]

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

[1.18]

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

[1.19]

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

[1.20]

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

[1.21]

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

[1.22]

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ć

[1.23]

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 [6]

[1.24]

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ń

[1.25]

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, E.@F.N. Spon, London 1987.

4. Osiadacz A.J., Statyczna symulacja sieci gazowych, Fluid Systems, Warszawa 2001.

5. Osiadacz A.J., Symulacja i optymalizacja sieci gazowych, [w:] 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 ai,j będący liczbą rzeczywistą lub zespoloną, lub operatorem (np. różniczkowania), bądź wielomianem. Macierz zapisujemy w postaci tablicy

[2.1]

lub krócej jako A = [ai,j], i = 1, 2, ..., m; j = 1, 2, ..., n; Am×n. Element ai,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)

[2.2]

można macierz (2.1) zapisać jako

[2.3]

Jeśli wiersze macierzy potraktujemy jako n-wymiarowe wektory wierszowe postaci Ai = [ai₁ ai₂ ... ... ain] (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, dii ≠ 0 (i = 1, 2, …, n), występują tylko na głównej przekątnej, nazywamy macierzą diagonalną i oznaczamy symbolem diag(d₁₁, d₂₂, …, dnn)

[2.4]

Jeżeli w macierzy diagonalnej wszystkie dii ≠ 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

[2.5]

Wprowadzając symbol Kroneckera

[2.6]

macierz jednostkową zapiszemy w postaci I = [δij].

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

Macierz kwadratowa, której elementy spełniają warunek aij = aji 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 = [aij]mn i B = [bij]mn są równe wtedy i tylko wtedy, gdy ich odpowiednie elementy są sobie równe, tzn. A = B, gdy [aij] = [bij] dla i = 1, 2, …, m; j = 1, 2, …, n.

Sumą (różnicą) dwóch macierzy A = [aij]mn i B = [bij]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 cij = [aij] ± [bij] 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 = [aij] o wymiarach m × n nazywamy macierz o wymiarach m × n, którą otrzymano w rezultacie zamiany wierszy z kolumnami. Elementy [aij] macierzy A zastąpiono elementami [aji] macierzy AT.

Dla dowolnej macierzy prostokątnej

macierz transponowana ma postać:

Wektor wierszowy xT = (x₁, x₂, …, xn) 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. aij = aji 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

[2.7]

gdzie Aij jest macierzą o wymiarach mi × ni (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 pi = qi (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ść

[2.8]

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|:

[2.9]

W przypadku macierzy prostokątnej A = [aij]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 Mij elementu aij 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 = ai₁Ai₁ + ai₂Ai₂ + … + ainAin (1 ≤ i ≤ n)

[2.10]

det A = a₁j A₁j + a₂j A₂j + … + anj Anj (1 ≤ j ≤ n)

[2.11]

gdzie A₁k (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₂₂, …, · lnn, det(U) = u₁₁· u₂₂, ..., · unn,

• 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 aij (i, j = 1, 2, …, n) macierzy A możemy obliczyć dopełnienie algebraiczne Aij = (–1)i⁺jMij, gdzie Mij 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 aij 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₂, …, xn 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 xi = yi dla i = 1, 2, …, m.

Sumą/różnicą dwóch wektorów nazywamy wektor przestrzeni m-wymiarowej o współrzędnych odpowiednio xi + yi lub xi – yi (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₂ + … + ckxk nazywany jest kombinacją liniową wektorów x₁, x₂, ..., xk. Wektory te są nazywane liniowo niezależnymi, jeśli równość c₁x₁ + c₂x₂ + … + ckxk = 0 zachodzi tylko pod warunkiem, że c₁ = c₂ = ... = ck = 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 ∈ Rn można więc wyrazić jako

x = a₁y₁ + a₂y₂ + … + anyn

[2.12]

gdzie a₁, a₂, …, an 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

[2.13]

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

[2.14]

lub w postaci macierzowo-wektorowej

[2.15]

gdzie A ‒ λI jest macierzą charakterystyczną dla macierzy A = [aij]n×n.

Wyznacznikiem charakterystycznym macierzy A nazywamy wyznacznik

[2.16]

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 ki (i = 1, 2, …, n) są funkcją elementów macierzy A:

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

ki – suma minorów głównych stopnia trzeciego itd., kn = 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

[2.17]

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