Metody niejawne jednokrokowe

Do tej grupy metod rozwiązania równań różniczkowych zwyczajnych zalicza się nie-jawną metodę Eulera oraz niejawną metodę trapezową. Obydwie metody można łatwo wyprowadzić, wykorzystując ogólny wzór (6.12).

Metodę niejawną Eulera otrzymujemy jako wynik przybliżenia pola trapezu krzywoli-niowego polem prostokąta o wymiarach h ⋅ f(xj+1, yj+1) (rys. 4.7), czyli

) , ( ))

( ,

( _{1 1}

1

+

⋅ +

∫

⁺ ≈ ^{j j}

x

x

y x f h dx x y x f

j

j

. (6.37)

Rys. 6.7. Zastąpienie pola pod krzywą f(x, y) polem prostokąta – niejawna metoda Eulera

Równanie (6.12) przyjmie więc postać

yj+1 = yj + h ⋅ f(xj+1, yj+1). (6.38)

112 6. Rozwiązywanie równań różniczkowych zwyczajnych

Romuald Szymkiewicz – Metody numeryczne w inżynierii wodnej Jeśli w równaniu (6.12) całkę zastąpimy polem trapezu prostoliniowego (rys. 6.8)

)) , ( ) , ( 2( )) ( ,

( 1 1

1

+

+ +

∫

⁺ ≈ ^{j j j j}

x

x

y x f y x h f dx x y x f

j

j

, (6.39)

to otrzymamy niejawną metodę trapezową

)) , ( ) , (

2( ^{1 1}

1 + +

+ = j+ j j + j j

j h f x y f x y

y

y . (6.40)

Jak widzimy, formuły metody niejawnej Eulera oraz niejawnej metody trapezów prowadzą do algebraicznych równań nieliniowych. Zatem w każdym kroku całkowania poszukiwaną wartość yj+1 otrzymamy, rozwiązując nieliniowe równanie algebraiczne dowolną metodą opisaną w rozdziale 2. Najczęściej stosowany wariant metody iteracji prostej, nazywany algorytmem „predyktor-korektor”, przedstawiono w punkcie 6.2.4.

Obydwie wymienione metody należą do najczęściej stosowanych metod rozwią zywa-nia układów równań różniczkowych zwyczajnych.

Rys. 6.8. Zastąpienie pola pod krzywą f(x, y) polem trapezu – niejawna metoda trapezowa

„

Przykład 6.2

Czas opróżniania zbiornika retencyjnego

Obliczyć czas konieczny do przygotowania zbiornika retencyjnego na przyjęcie fali wezbraniowej, to znaczy czas, w którym zwierciadło wody w zbiorniku H(t) z począ tkowe-go poziomu Hp (rys. 6.2.1) obniży się do poziomu pożądanego Hk.

Obliczenie należy wykonać, przyjmując następujące założenia i dane:

⎯ pocz^ątkowy poziom wody w zbiorniku wynosi Hp = 15 m npp, zaś poziom, do którego należy obniżyć zwierciadło wody, wynosi Hk = 10 m npp;

⎯ dopływ wody do zbiornika nie zmienia si^ę w czasie i wynosi q = 3 m³/s;

⎯ odpływ wody odbywa si^ę przez upust denny przy jego całkowitym otwarciu. Jego po-wierzchnia wynosi A = 2,50 m², zaś współczynnik wydatku ϕ = 0,65;

⎯ poziom wody na stanowisku dolnym jest stały i wynosi Hd = 5 m npp;

⎯ powierzchni^ę zbiornika na poziomie zwierciadła wody, będącą funkcją napełnienia zbiornika F(H), definiuje tabela 6.2.1 i rys. 6.2.2.

Hp H_k

H_d poziom porównawczy

Q q

Rys. 6.2.1. Schemat zbiornika retencyjnego

Tabela 6.2.1 Zależność powierzchni zbiornika od jego napełnienia

Hi [m npp] 5,00 6,00 7,40 8,30 10,40 13,25 16,00 Fi [m²] 1000 20000 40000 60000 120000 240000 400000

Dla pośrednich wartości H wartość F(H) oblicza się, interpolując liniowo pomiędzy punk-tami. Zatem F(H) zdefiniowana jest następująco:

( 1,2,3,4,5,6,7)

. dla

dla ) (

dla )

(

6 6

1 1

1

1 1

=

⎪⎪

⎩

⎪⎪⎨

⎧

≥

<

<

− − + −

≤

= ₊

+

+ i

H H F

H H H H

H H H

F F F

H H F

H

F _{i i i}

i i

i i

i (6.2.1)

F[x10 m ]^{4 2} 2 6 10 14 18 22 26 30 34 38 H

[m npm]

16 14 12 10 8 6

Rys. 6.2.2. Zależność powierzchni zbiornika od jego napełnienia

114 6. Rozwiązywanie równań różniczkowych zwyczajnych

Romuald Szymkiewicz – Metody numeryczne w inżynierii wodnej Dla zbiornika przedstawionego na rys. 6.2.1 można napisać następujące różniczkowe równanie retencji:

))

w którym chwilowy odpływ Q(t) jest funkcją różnicy poziomów wody w zbiorniku i poni-żej zapory:

Do rozwiązania zastosujmy niejawny schemat Eulera (6.38), który ma postać:

1

1 +

+ = j+Δ ′j

j H tH

H , (6.2.4)

gdzie: j − indeks punktu obliczeniowego, Δt − krok całkowania w czasie.

Ponieważ zastosowany schemat całkowania jest niejawny, otrzymane równanie jest nieliniowe. Zastosujmy więc metodę iteracji prostej. Załóżmy jednocześnie, że pierwsze przybliżenie obliczymy jawną metodą Eulera (6.15):

(

¹²

)

Tak obliczony wynik korygujemy, stosując (6.2.5) w postaci:

(

^{( )1 12}

)

Proces kontynuujemy do chwili, gdy różnica wyników w dwóch kolejnych iteracjach spełni warunek:

gdzie: ε− przyjęta dokładność rozwiązania.

Zatem odpowiedź na pytanie o czas opróżnienia zbiornika otrzymujemy, rozwiązując rów-nanie różniczkowe zwyczajne (6.2.2) z warunkiem początkowym:

H(t = 0) = H_p. (6.2.9)

Nie interesuje nas postać funkcji H(t), lecz czas t, po którym zwierciadło wody w zbiorniku osiągnie wartość H_k. Równanie różniczkowe należy więc całkować z krokiem Δt do chwili spełnienia warunku

H(t) ≤ Hk (6.2.10)

Algorytm obliczeń można zapisać w następującej ogólnej postaci:

1) podstaw j = 1, H_j = H_p;

2) oblicz pierwsze przybliżenie H⁽_j⁰₊⁾₁ według (6.2.6);

3) oblicz następne przybliżenie H⁽_j¹₊⁾₁ według (6.2.7);

4) sprawdź warunek (6.2.8):

⎯ jeśli jest spełniony:

podstaw: H_j+1 = H⁽_j¹₊⁾₁ na j podstaw j + 1 i przejdź do punktu 5,

⎯ jeśli nie jest spełniony:

podstaw: H⁽_j⁰₊⁾₁ = H⁽_j¹₊⁾₁ i przejdź do punktu 3;

5) sprawdź warunek (6.2.10);

⎯ jeśli nie jest spełniony: przejdź do punktu 2;

⎯ jeśli jest spełniony, zakończ obliczenia, ponieważ t = (j – 1) ⋅ Δt jest poszukiwanym czasem, po którym zwierciadło wody osiągnie poziom H_k.

Obliczenia wykonane według powyższego algorytmu dla przyjętych danych wykazały, że przygotowanie zbiornika na przyjęcie fali wezbraniowej wymaga ok. 17 godzin i 36 minut. Wynik ten uzyskano, przyjmując dokładność rozwiązania równania nieliniowego ε = 0,0005 m. Krok całkowania Δt zmieniano w kolejnych wariantach zadania w granicach od Δt = 180 s do Δt = 720 s. Zmiany te nie miały istotnego wpływu na wynik ostateczny. „

„

Przykład 6.3

Obliczenie krzywej spiętrzenia

Obliczenie krzywej spiętrzenia polega na wyznaczeniu profilu zwierciadła wody w ka-nale, który wywołuje budowla piętrząca przy danym natężeniu przepływu (rys. 6.3.1).

W tym celu należy scałkować równanie różniczkowe zwyczajne opisujące ruch ustalony wolnozmienny w kanale z warunkiem początkowym odpowiadającym spiętrzeniu w prze-kroju posadowienia budowli. W układzie współrzędnych, jak na rys. 6.3.1, równanie to można zapisać w następującej postaci (Szymkiewicz, 2000):

dx S

dE = , (6.3.1) gdzie: E(x) – całkowita energia mechaniczna liczona względem przyjętego poziomu odniesienia,

S – spadek linii energii.

Energię mechaniczną E definiuje wyrażenie

2 2

2gA h Q

E_{= +} α _{, (6.3.2)}

zaś spadek linii energii przy wyrażeniu oporów ruchu przez równanie Manninga określa wzór

2 3 4

2 2

A R

Q

S = n , (6.3.3)

gdzie: h(x) – rzędna zwierciadła wody w punkcie x, α – współczynnik de Saint-Venanta,

116 6. Rozwiązywanie równań różniczkowych zwyczajnych

Romuald Szymkiewicz – Metody numeryczne w inżynierii wodnej Q – natężenie przepływu,

g – przyspieszenie grawitacyjne,

A – powierzchnia przekroju czynnego kanału, R – promień hydrauliczny przekroju czynnego, n – współczynnik szorstkości kanału wg Manninga.

x x = 0

Rys. 6.3.1. Profil zwierciadła wody wywołany budowlą piętrzącą

Równanie (6.3.1) rozwiązujemy niejawną metodą trapezową (6.40) )

gdzie: j – indeks przekroju obliczeniowego, Δx – krok całkowania.

Podstawiając do powyższego wzoru wyrażenia (6.3.2) i (6.3.3), otrzymujemy

⎟⎟

to tym samym w równaniu (6.3.5) wszystkie zmienne z indeksem j są znane. Równanie to jest więc równaniem z jedną niewiadomą hj+1, reprezentującą rzędną zwierciadła wody w przekroju 2 odległym od zapory o Δx. Równanie jest jednak nieliniowe, gdy^ż Aj+1 = A(hj+1) oraz Rj+1 = R(hj+1). Rozwiązując je jedną z metod rozwiązywania równań algebraicznych nieliniowych, otrzymujemy przybliżoną wartość h2. Następnie identyczny tok postępowania powtarzamy dla następnego odcinka kanału. Obliczenia prowadzimy do miejsca, w którym spełniony jest następujący warunek:

m

gdzie: Y_n – głębokość normalna w kanale odpowiadająca natężeniu przepływu Q.

Powyższy warunek określa koniec krzywej spiętrzenia.

Tabela 6.3.1 Obliczony profil zwierciadła wody w kanale

Indeks

przekroju km Rzędna dna r [m npp]

Rzędna zw.

wody h [ m npp]

Głębokość y [m]

1 0,000 5,000 8,500 3,500 2 0,250 5,063 8,501 3,438 3 0,500 5,125 8,501 3,376 4 0,750 5,188 8,502 3,314 5 1,000 5,250 8,503 3,253 6 1,250 5,313 8,504 3,191 7 1,500 5,375 8,504 3,129 8 1,750 5,438 8,505 3,068 9 2,000 5,500 8,506 3,006 10 2,250 5,563 8,508 2,945 11 2,500 5,625 8,509 2,884 12 2,750 5,688 8,510 2,823 13 3,000 5,750 8,512 2,762 14 3,250 5,813 8,513 2,701 15 3,500 5,875 8,515 2,640 16 3,750 5,938 8,517 2,579 17 4,000 6,000 8,519 2,519 18 4,250 6,063 8,521 2,459 19 4,500 6,125 8,524 2,399 20 4,750 6,188 8,527 2,339 21 5,000 6,250 8,530 2,280 22 5,250 6,313 8,533 2,221 23 5,500 6,375 8,537 2,162 24 5,750 6,438 8,541 2,104 25 6,000 6,500 8,546 2,046 26 6,250 6,563 8,551 1,988 27 6,500 6,625 8,557 1,932 28 6,750 6,688 8,564 1,876 29 7,000 6,750 8,571 1,821 30 7,250 6,813 8,579 1,767 31 7,500 6,875 8,588 1,713 32 7,750 6,938 8,598 1,661 33 8,000 7,000 8,610 1,610 34 8,250 7,063 8,623 1,560 35 8,500 7,125 8,637 1,512

118 6. Rozwiązywanie równań różniczkowych zwyczajnych

Romuald Szymkiewicz – Metody numeryczne w inżynierii wodnej Opisany algorytm rozwiązania zastosowano do wyznaczenia krzywej spiętrzenia w ka-nale trapezowym o następujących parametrach:

⎯ szeroko^ść kanału B = 5 m,

⎯ nachylenie skarp kanału M = 1,5,

⎯ spadek dna kanału s = 0,00025,

⎯ współczynnik szorstko^ści wg Manninga n = 0,030,

⎯ gł^ębokość normalna Yn = 1,5 m,

⎯ nat^ężenie przepływu Q = 3 m³/s,

⎯ rz^ędna zwierciadła wody po podpiętrzeniu w przekroju budowli: hp = 8,5 m.

Wyniki obliczeń, uzyskane przy kroku całkowania Δx = 250 m oraz dokładno^ści roz-wiązania równania nieliniowego (6.3.5) metodą siecznych ε = 0,0005, przedstawiono

w tabeli 6.3.1. „

W dokumencie Metody numeryczne w inżynierii wodnej (Stron 112-119)