Metody numeryczne - Metoda różnic skończonych

Treść

Rozwiązać problem brzegowy z danymi jak poniżej

𝑦′′⁡(𝑥)=2⁢𝑥𝑥∈[𝑎.𝑏]𝑎=3𝑏=7𝑦⁡(𝑎)=4𝑦′⁡(𝑏)=2

Dzieląc przedział na 𝑛 =4 elementy

Rozwiązanie

Długość pojedynczego elementu:

ℎ=𝑏−𝑎𝑛=1

Położenie poszczególnych węzłów

𝑥𝑖=𝑥0+𝑖⋅ℎ

Zamieniamy równanie 𝑦″⁡(𝑥) =2⁢𝑥 korzystając ze wzorów różnicowych:

1ℎ2⁢(𝑖−1−2⋅𝑖+𝑖𝑝⁢1)=2⁢𝑥𝑖−1−2⋅𝑖+𝑖𝑝⁢1=ℎ2⋅2⁢𝑥

Dla każdego węzła 𝑖 =1,2,3,4 zapisujemy równanie różnicowe zgodnie z powyższym wzorem, pomocniczo zapisujemy wartośći x dla poszczególnych węzłów:

𝑖=1𝑖0−2⋅𝑖1+𝑖2=ℎ2⋅2⁢𝑥1𝑥1=𝑥0+ℎ=4𝑖=2𝑖1−2⋅𝑖2+𝑖3=ℎ2⋅2⁢𝑥2𝑥2=𝑥0+2⁢ℎ=5𝑖=3𝑖2−2⋅𝑖3+𝑖4=ℎ2⋅2⁢𝑥3𝑥3=𝑥0+3⁢ℎ=6𝑖=4𝑖3−2⋅𝑖4+𝑖5=ℎ2⋅2⁢𝑥4𝑥4=𝑥0+4⁢ℎ=7

Uwzględniamy warunki brzegowe, również zamieniając je na wzory różnicowe

𝑦⁡(𝑎)=4𝑦0=4𝑦′⁡(𝑏)=212⁢ℎ⁢(−𝑖3+𝑖5)=2−𝑖3+𝑖5=4⁢ℎ

Zapisujemy powyższe równania w postaci macierzowej 𝐴 ⋅𝑦 =𝐵

⎡⎢⎣1−2100001−2100001−2100001−21100000000−101⎤⎥⎦ ⋅⎡⎢⎣𝑦0𝑦1𝑦2𝑦3𝑦4𝑦5⎤⎥⎦ =⎡⎢⎣ℎ2⋅2⋅𝑥1ℎ2⋅2⋅𝑥2ℎ2⋅2⋅𝑥3ℎ2⋅2⋅𝑥444⋅ℎ⎤⎥⎦ =⎡⎢⎣810121444⎤⎥⎦

I ostatecznie rozwiązanie (w tym przypadku uzyskane w programie Mathcad):

𝑦 =𝐴−1⁢𝐵 =⎡⎢⎣4−31−58−75−80−71⎤⎥⎦

W praktyce interesuje nas jedynie pierwsze 5 wyników, ponieważ ostatni znajduje się de facto poza przedziałem badanym 𝑖 =5 − >𝑥 =8

Rozwiązanie analitycznego tego równania uzyskane klasycznymi metodami rozwiązywania równań różniczkowych:

𝑦𝑎⁢𝑛⁡(𝑥) =𝑥33 −47⁢𝑥 +136

W tabeli zestawiono wartości wraz z błędami dla poszczególnych punktów

i𝑥𝑦𝑦𝑎 Bład bezwzgledny  Bład wzgledny 034400,00%14−31−30,6670,3331,09%25−58−57,3330,6671,16%36−75−7411,35%47−80−78,6671,3331,69%

Ddybyśmy podzielili przedział na więcej elementów dokładność była by oczywiście większa, natomiast w takim przypadku konieczne byłoby użycie programów obliczeniowych typu Matlab

Treść

Rozwiązać podaną belkę korzystając z MRS, 𝐸⁢𝐼 =8000⁢𝑘⁢𝑁⁢𝑚2

Rozwiązanie

Zapisujemy równanie różnicowe dla punktów i=1,2,3,4

𝑖=11⁢𝑣−1−4⁢𝑣0+6⁢𝑣1−4⁢𝑣2+1⁢𝑣3=𝑏𝑖=21⁢𝑣0−4⁢𝑣1+6⁢𝑣2−4⁢𝑣3+1⁢𝑣4=𝑏𝑖=31⁢𝑣1−4⁢𝑣2+6⁢𝑣3−4⁢𝑣4+1⁢𝑣5=𝑏𝑖=41⁢𝑣2−4⁢𝑣3+6⁢𝑣4−4⁢𝑣5+1⁢𝑣6=𝑏

gdzie 𝑏 =ℎ4⋅𝑞𝐸⁢𝐼

Następnie zapisujemy warunki brzegowe, w x=0, i=0 utwierdzenie pełne, w x=4, i=4 swobodny koniec

𝑣0=0−1⁢𝑣−1+1⁢𝑣1=01⁢𝑣3−2⁢𝑣4+1⁢𝑣5=0−1⁢𝑣2+2⁢𝑣3−2⁢𝑣5+1⁢𝑣6=0

Podobnie jak w poprzednim przykładzie zapisujemy całość w postaci macierzowej i rozwiązujemy przy pomocy programu obliczeniowego

⎡⎢⎣1−46−4100001−46−4100001−46−4100001−46−4101000000−1010000000001−210000−120−21⎤⎥⎦ ⋅⎡⎢⎣𝑣−1𝑣0𝑣1𝑣2𝑣3𝑣4𝑣5𝑣6⎤⎥⎦ =⎡⎢⎣𝑏𝑏𝑏𝑏0000⎤⎥⎦ =⎡⎢⎣1.25⋅10−31.25⋅10−31.25⋅10−31.25⋅10−30000⎤⎥⎦
⎡⎢⎣𝑣−1𝑣0𝑣1𝑣2𝑣3𝑣4𝑣5𝑣6⎤⎥⎦ =𝐴−1⁢𝐵 =⎡⎢⎣5⋅10−305⋅10−31.562⋅10−22.875⋅10−24.25⋅10−25.625⋅10−27.062⋅10−2⎤⎥⎦

Mając wartości ugięcia możemy obliczyć np. wartość momentu gnącego w utwierdzeniu

𝑀0 =−𝐸⁢𝐼ℎ2⁢(𝑣−1−2⋅𝑣0+𝑣1) =−80⁢𝑘⁢𝑁⁢𝑀

Akurat ta belka jest bardzo typowa, możemy więc z łatwością znaleźć analityczne rozwiązanie i porównać je z naszym, numerycznym

𝑀𝑎 =−12⁢𝑞 ⋅𝐿2 =−80⁢𝑘⁢𝑁⁢𝑚

Jak widać rozwiązanie dla momentu zginającego jest identyczne, nieco bardziej interesująco ma się sprawa z ugięciem końca belki:

𝑣𝑎=𝑞⋅𝐿48⁢𝐸⁢𝐼=4⋅10−2𝑣4−𝑣𝑎𝑣𝑎=6.25%

Jak widać w tym przypadku błąd względny wynosi ponad 6%, podobnie jak poprzednio wyższe n sprawiłoby że ten błąd byłby mniejszy.

Bonus

W powyższym przykładzie pominęliśmy wartość siły tnącej w utwierdzeniu. Jak łatwo zauważyć stwarza ona pewien kłopot, ponieważ wzór:

𝑇0 =−𝐸⁢𝐼2⁢ℎ3⁢(−𝑣−2+2⋅𝑣−1−2⁢𝑣1+𝑣2)

zawiera wartość przemieszczenia 𝑣−2 której nie obliczyliśmy. Jeżeli jednak w zadaniu mielibyśmy za zadanie obliczyć wartość siły tnącej musielibyśmy dodatkow zapisać wzór różnicowy dla i=0, co zilustrowano poniżej

⎡⎢⎣1−46−41000001−46−41000001−46−41000001−46−41000001−46−410010000000−10100000000001−2100000−120−21⎤⎥⎦ ⋅⎡⎢⎣𝑣−2𝑣−1𝑣0𝑣1𝑣2𝑣3𝑣4𝑣5𝑣6⎤⎥⎦ =⎡⎢⎣𝑏𝑏𝑏𝑏𝑏0000⎤⎥⎦
⎡⎢⎣𝑣−2𝑣−1𝑣0𝑣1𝑣2𝑣3𝑣4𝑣5𝑣6⎤⎥⎦ : =𝐴−1⁢𝐵 =⎡⎢⎣2.562⋅10−25⋅10−305⋅10−31.562⋅10−22.875⋅10−24.25⋅10−25.625⋅10−27.062⋅10−2⎤⎥⎦

Oczywiście pozostałe wartości przemieszczenia nie zmieniły się, ale mając obliczone 𝑣−2 możemy zweryfikować wartość siły tnącej w utwierdzeniu:

𝑇𝑎=𝑞⋅𝐿=40⁢𝑘⁢𝑁𝑇0=−𝐸⁢𝐼2⁢ℎ3⁢(−𝑣−2+2⋅𝑣−1−2⁢𝑣1+𝑣2)=40⁢𝑘⁢𝑁

Jak widać wynik jest idelanie zgodny z wartością analityczną