Ćwiczenie 7a: Równanie falowe dla struny

Ćwiczenie 7a: Równanie falowe dla struny ^∗

Równanie struny (przy jednostkowej prędkości rozchodzenia się drgań) ma postać

∂ ² u

∂t ² = ∂ ² u

∂x ² , (1)

gdzie u(x, t) oznacza wychylenie struny w punkcie x w chwili t. Jednoznaczne rozwiązania równania (1) można podać, jeśli zadane zostaną warunki początko- we na wychylenie

u(x, t = 0) = u 0 (x) (2)

i prędkość w chwili początkowej t = 0:

v(x, t = 0) ≡ ∂u(x, t = 0)

∂t |

= v 0 (x). (3)

Struna jest sztywno zamocowana na końcach: u(x = 0, t) = u(x = 1, t) = 0.

Rozwiązanie równania falowego można przedstawić w rozwinięciu na drgania własne

u(x, t) =

∑

c

sin(ω

x) cos(ω

t) +

∑

s

sin(ω

x) sin(ω

t) (4)

z ω

= nπ. Współczynniki rozwinięcia c

i s

wyznaczamy z warunku począt- kowego korzystając z ortogonalności drgań normalnych

c

= 2

∫ 1 0

u 0 (x) sin(ω

x)dx (5)

s

= 2 ω

∫ 1 0

v ₀ (x) sin(ω

x)dx (6)

Dla warunków początkowych:

u 0 (x) = u(x, t = 0) = sin(πx) − 1

2 sin(2πx) (7)

oraz

v 0 (x) = v(x, t = 0) = 0 (8)

∗

Laboratorium z inżynierskich metod numerycznych, Wydział Fizyki i Informatyki Sto-

sowanej AGH 2014/2015. Bartłomiej Szafran (bszafran@agh.edu.pl), Krzysztof Kolasiński

(kolasinski@fis.agh.edu.pl), Elżbieta Wach (Elzbieta.Wach@fis.agh.edu.pl), Dariusz Żebrowski

(Dariusz.Zebrowski@fis.agh.edu.pl)

rozwiązanie analityczne ma postać:

u(x, t) = cos(πt) sin(πx) − 1

2 cos(2πt) sin(2πx). (9) Z tego wynikają analityczne wzory na prędkość i przyspieszenie:

v(x, t) = −π sin(πt) sin(πx) + π sin(2πt) sin(2πx), (10) a(x, t) = −π ² cos(πt) sin(πx) + 2π ² cos(2πt) sin(2πx). (11)

Zadanie 1: Schemat Verleta

Równanie (1) rozwiążemy metodą różnic skończonych na siatce N = 101 punktów w przedziale x ∈ [0, 1], z krokiem czasowym ∆t = 1/200. Przyjmiemy warunki początkowe i brzegowe odpowiadające rozwiązaniu (9). Użyjemy pręd- kościowego schematu Verleta. Schemat Verleta dany jest przez trzy równania:

v(x, t + ∆t/2) = v(x, t) + ∆t

2 a(x, t) (12)

u(x, t + ∆t) = u(x, t) + ∆tv(x, t + ∆t/2) (13) v(x, t + ∆t) = v(x, t + ∆t/2) + ∆t

2 a(x, t + ∆t), (14)

gdzie przyspieszenie a =

wyliczamy z prawej strony równania falowego (1) jako

a(x, t) = u(x + ∆x, t) + u(x − ∆x, t) − 2u(x, t)

∆x ² . (15)

Kroki algorytmu:

1. ustalenie warunku początkowego/brzegowego na wychylenie u(x, t = 0), prędkość v(x, t = 0), przyspieszenie a(x, t = 0): wzory (7),(8) oraz (11);

2. obliczenie prędkości pośredniej v(x, t + ∆t/2) dla całej struny: (12);

3. znalezienie wychylenia struny u(x, t + ∆t) w nowej chwili czasowej: (13);

4. obliczenie przyspieszenia a(x, t + ∆t) (15);

5. wyliczenie nowego rozkładu prędkości v(x, t + ∆t) wzdłuż struny: (14).

6. przejście do kolejnej chwili czasowej → krok 2.

a) Rozwiązać równanie (1) dla warunków początkowych

u 0 (x) = exp( −100(x − 0.5) ² ), (16) oraz

v ₀ (x) = 0 (17)

przy sztywnych warunkach brzegowych, tj. u(x = 0, t) = u(x = 1, t) = 0 (jak wszędzie powyżej) dla t ∈ [0, 4]. Narysować mapę u(x, t) ( 10 pkt) [na osi poziomej przyjąć t, na pionowej x].

b) Wykonać obliczenia z podpunktu a) z innymi warunkami brzegowymi: tym razem końce struny nie są przytwierdzone na stałe, tylko mogą się swobodnie poruszać:

(x = 0, t) =

(x = 1, t) = 0, dla t ∈ [0, 4]. Luźne warunki brzegowe realizujemy przypisując po każdej chwili czasowej u(x = 0, t) :=

u(x = ∆x, t) oraz u(x = 1, t) := u(x = 1 −∆x, t) - wówczas pochodna będzie znikać.

Narysować wychylenie początku i końca struny w funkcji czasu, tj. u(x = 0, t) oraz u(x = 1, t) (10 pkt).

c) Wracamy do sztywnych warunków brzegowych. Zbadajmy ruch identyczne- go pakietu w przypadku dwóch różnych ośrodków. Modyfikujemy równanie struny:

∂ ² u

∂t ² = 1 ρ(x)

∂ ² u

∂x ² , (18)

gdzie ρ(x) jest gęstością ośrodka. Definiujemy:

ρ(x) =

{ 1 dla x < x ₀ ρ ₀ dla x ­ x 0

(19) (x ₀ to granica ośrodków). Niejednostkowa gęstość powoduje zmianę we wzo- rze na przyspieszenie - teraz mamy:

a(x, t) = u(x + ∆x, t) + u(x − ∆x, t) − 2u(x, t)

∆x ² · 1

ρ(x) . (20)

Należy przeprowadzić identyczne obliczenia, jak w podpunkcie a) - zamie-

niamy tylko wzór na przyspieszenie. Przyjąć x 0 = 0.75, ρ 0 = 3.0. Narysować

mapę u(x, t) (10 pkt).

Zadanie 3: Drgania tłumione

Wracamy do równania (1) - tj. pomijamy gęstość ośrodka. Wprowadzimy tłumienie drgań proporcjonalne do prędkości struny

∂ ² u

∂t ² = ∂ ² u

∂x ² − 2β ∂u

∂t , (21)

gdzie β jest współczynnikiem tłumienia. Przyspieszenie potrzebne do schematu Verleta dane jest teraz przez

a(x, t) = u(x + ∆x, t) + u(x − ∆x, t) − 2u(x, t)

∆x ² − 2β u(x, t) − u(x, t − ∆t)

∆t .

(22) Zmiana w programie polega więc znowu tylko na podstawieniu nowego wzoru na przyspieszenie. Narysować trzy mapy u(x, t): dla β = 0.2, 1.0 i 3.0 przy warun- kach początkowych (16), (17) oraz sztywnych warunkach brzegowych (10 pkt).

Zadanie 4: Drgania wymuszone

Do poprzedniego zadania dodajemy przyspieszenie związane z zewnętrzną siłą a

(x, t).

∂ ² u

∂t ² = ∂ ² u

∂x ² − 2β ∂u

∂t + a

(x, t), (23)

Siłę przykładamy punktowo na środku struny:

a

(x, t) = cos(ωt) exp( −10 ⁵ (x − x 0 ) ² ) (24) z x ₀ = 0.5 (czynnik −10 ⁵ w eksponencie gwarantuje wygaszenie wartości dla x różnych od x ₀ ). W chwili początkowej struna spoczywa (v ₀ (x) = 0) w rów- nowadze (u 0 (x) = 0). Wstawimy słabe tłumienie β = 1 i częstość wymuszenia ω = π/2. Narysować mapę u(x, t) dla t ∈ [0, 10] ( 10 pkt). Widzimy, że drga- nia osiągają stan ustalony po pewnym czasie. Jaki jest okres drgań w stanie ustalonym?

Zadanie 5: Rezonanse

Energia struny w chwili t dana jest wyrażeniem

dla t 1 = 16, t 2 = 20 w przypadku z zadania 4. Narysować ⟨E⟩ w funkcji

ω ∈ (0, 10π) ( 10 pkt). Pętla zewnętrzna powinna przebiegać po wartościach

ω; w środku pętla po czasie t, a w niej kroki algorytmu.