Ćwiczenie 9 Zastosowanie metody różnic skończonych do obliczania
płyt zginanych (opracowali M. Kłos Z. Waszczyszyn) wersja dla
studentów
Podstawowe wzory metody różnic skończonych (MRS)
1. Wzory różnicowe dla funkcji jednej zmiennej
y =� f(�x)�
Zakładamy, że funkcja
jest ciągła tak, że możemy obliczyć jej pochodne
w węz
le i na osi zmiennej niezależnej x. Posługujemy się skróconymi oznaczeniami
widocznymi na rys.9.1.
y
2Dx
x
X Xi-1 Xi Xi+1
i-2
i-2 i-1 i i+1 i+2
Dx Dx Dx Dx
Rys.9.1
f(�x)�
Funkcję rozwijamy w szereg Taylora w otoczeniu węzła i :
Dx2 Dx3 Dx4
f(�x)� =� fi +� fi'Dx +� fi" +� fi'" +� fiIV +� ...,
(9.1)
2 3! 4!
gdzie:
dfi d2fi
fi' =� , fi" =�
itd.
x=�xi x=�xi
dx
dx2
1
i-1
i
i+1
i+2
f
f
f
f
)
x
(
f
i =� ... , i -� 2, i -�1, i, i +�1, ...
Przyjmujemy równą odległość węzłów , określoną przez
przyrost tak, że
Dx
xi-�2 =� xi -� 2Dx, xi-�1 =� xi -� Dx, xi+�1 =� xi +� Dx, xi+�2 =� xi +� 2Dx
itd.
i +�1, i -�1
Zgodnie ze wzorem (9.1) obliczamy wartości funkcji w węzłach :
Dx2 Dx3 Dx4
fi+�1 =� fi +� fi'Dx +� fi" +� fi'" +� fiIV +� ....
2 6 24
(9.2)
Dx2 Dx3 Dx4
fi-�1 =� fi -� fi'Dx +� fi" -� fi'" +� fiIV +� ....
2 6 24
Jeżeli odejmiemy stronami zależności (9.2) to otrzymamy wzór na 1-szą pochodną:
fi' =�
(9.3)
Stąd otrzymujemy wzór różnicowy
fi+�1 -� fi-�1 Dfi
fi' � =� (9.4)
2Dx 2Dx
fi' =� (�df / dx)�
W stosunku do ścisłego wzoru analitycznego wzór różnicowy (9.4)
x=�xi
0(�-� fi'"Dx2 / 6)�
ma błąd o rzędzie wynikającym z pierwszego z pomijanych wyrazów
szeregu Taylora. Należy dodać, że na ogół piszemy znak � zamiast = w (9.4).
Jeśli dodamy stronami zależność (9.2) to otrzymamy wzór różnicowy na drugą
pochodną
fi" =�
(9.5)
W podobny sposób możemy wyprowadzić wzory różnicowe na wyższe pochodne
fi'", fiIV itd. Istotą tych wzorów jest, że pochodne obliczamy za pomocą wartości
funkcji w sąsiednich węzłach.
2
2. Wzory różnicowe dla funkcji dwóch zmiennych
f(�x, y)� (�x, y)�
W przypadku funkcji dziedzina funkcji jest określona na płaszczyz
nie
i węzły otrzymujemy przez przecięcie linii siatki różnicowej. Na rys.9.2 pokazano
Dx�� Dy
siatkę o oczkach prostokątnych . W przypadku funkcji 2-ch zmiennych
k, i
obliczamy pochodne cząstkowe jako pochodne kierunkowe wzdłuż linii
ś�f ś�f 1
=� � (�fi+�1, k -� fi-�1,k )�,
xi ,yk i,k
ś�x ś�x 2Dx
(9.6)
ś�f ś�f 1
=� � (�fi, k+�1 -� fi,k-�1)�,
xi ,yk i,k
ś�y ś�y 2Dy
i,k+2
k+2
y
i-1,k+1 i,k+1 i+1,k+1
k+1
i-2,k i-1,k i,k i+2,k x
i+1,k
k
i,k-1
i-1,k-1 i+1,k-1 i+2,k-1
k-1
i,k-2
k-2
Dx Dx Dx Dx
Rys.9.2
Wzór (9.5) wykorzystujemy dla obliczenia 2-gich pochodnych cząstkowych.
ś�2 f 1
� (� fi+�1,k -� 2 fi,k +� fi-�1,k )�,
ś�x2 i,k Dx2
(9.7)
ś�2 f
�
ś�y2 i,k
Pochodną mieszaną liczymy za pomocą wzoru (9.8):
3
i-2
i-1
i
i+1
i+2
D
y
D
y
D
y
D
y
D
y
ć� �� ć�
ś�2f ś� ś�f ś� 1 1 1
=� �� �� �� (�fi,k+�1 -� fi,k-�1)��� =� (�fi+�1,k+�1 -� fi-�1,k+�1)� -�
� ��
i,k i,k
�� �� �� ��
ś�xś�y ś�x ś�y ś�x Dy Dx Dy
Ł� ł� Ł� ł�
(9.8)
1 1 1
-� (�fi+�1,k-�1 -� fi-�1,k-�1)� =� (�fi+�1,k+�1 +� fi-�1,k-�1 -� fi-�1,k+�1 -� fi+�1,k-�1)�.
Dx Dy DxDy
Wzory różnicowe można schematycznie przedstawić w postaci tzw. gwiazd
różnicowych. W gwiazdach piszemy wartości współczynników w miejscach
występowania wartości funkcji. Dalej piszemy gwiazdy w odniesieniu do siatki
l� =� Dx =� Dy
kwadratowej i ustalonej długości oczka siatki . Dla pochodnych rzędu 1 i
2 gwiazdy przyjmują postać:
1
ś�f 1 ś�f 1
i, k
- 1 1
=� �� =� ��
i,k i,k
i , k
ś�x 2l� ś�y 2l�
- 1
(9.9)
1
i , k
ś�2 f 1 ś�2 f 1 - 2
1 - 2 1
=� �� =� ��
ś�x2 i,k l�2 ś�y2 i,k l�2
i, k
1
- 1
1
ś�2f 1
i , k
=� ��
i,k
ś�x ś�y
l�2
W podobny sposób dochodzimy do
1 gwiazd dla laplasjanu i biilaplasjanu:
- 1
4
a b
1
1
2 - 8 2
1
1 - 4 1
�2 f =� ��
i, j
i , k
l�2
1
1 - 8 2 0 - 8 1
�2�2 f =� ��
i, j
1
l�4
2 - 8 2
1
(9.10)
3.Wzory różnicowe dla zginanej płyty
Przygotowane wzory różnicowe możemy zastosować do równania płyty
p(�x, y)�
�2�2w =� ,
D
(9.11)
które za pomocą gwiazdy (9.8b)
pi,k
�2�2w =� , (9.12)
i,k
D
piszemy dla każdego węzła siatki różnicowej na płaszczyz
nie środkowej płyty.
Pojawia się przy tym konieczność wprowadzenia węzłów fikcyjnych poza obszarem
płyty. Jeśli ograniczymy się do brzegów podpartych niepodatnie to będziemy mieli
dwa przypadki warunków brzegowych:
a)podparcie przegubowe, b) podparcie utwierdzone
w n ę t r z e
w n ę t r z e
b - 1 b b + 1
b - 1 b b + 1
p ł y t y
p ły t y
b - 1 b + 1
b
b - 1
b b + 1
5
b
b
W = 0
W = 0
b + 1
W
b
b - 1
b + 1
W
- 1
W
W
Rys. 9.3
Jeśli napiszemy warunek przegubowego podparcia
ś�2w
mx b =� 0�� =� 0
ś�n2 b
to po podstawieniu wzoru różnicowego (9.5) i uwzględnieniu warunku niepodatnego
w =� 0
podparcia otrzymujemy zgodnie z rys.9.3a
b
w =� -�w
a) . (9.13)
b+�1 b-�1
j�b =� 0
Tak samo z warunku b) otrzymujemy dla
wb+�1 =� wb-�1
b) (9.14)
pi,k
Pozostaje do rozwiązania dyskretna wartość obciążenia we wzorze (9.12).
Obliczamy ją jako wartość uśrednioną dzieląc wypadkową z obciążenia
przykładanego do oczka siatki różnicowej przez pole powierzchni siatki. Na rys.9.4
pokazano przypadek częściowego obciążenia oczka siatki obciążeniem
pi,k
równomiernie rozłożonym i obciążenie siłą skupioną. Obliczone wartości
wynoszą:
p0 P
pi,k =�
a) pi,k =� , oraz b) (9.15)
2 l�2
a) b)
p0
P
l� l�
l� l�
Rys.9.4
6
Po obliczeniu wartości ugięć w węzłach siatki różnicowej płaszczyzny środkowej
płyty możemy obliczyć w tych węzłach wartości sił przekrojowych. Gwiazdy
w
różnicowe rysujemy też dla sił wewnętrznych wyrażonych przez funkcję ugięcia
(i,k).
W odniesieniu do momentów obowiązują wzory:
ć� �� ć� ��
ś�2w ś�2w ś�2w ś�2w 1-� ś�2w
�� �� �� ��
mx =� -�D�� +� my =� -�D�� +� mxy =� -� D ,
(9.16)
ś�x2 ś�y2 ��, ś�y2 ś�x2 ��, 2 ś�xś�y
Ł� ł� Ł� ł�
dla których możemy narysować następujące gwiazdy różnicowe
-1
-�
i,k i,k
D -� D
mx i,k =� -� �� my i,k =� ��
-� 2+2
-�
l�2 -1 2+2 -1 l�2
-1
-�
(9.17)
1 - 1
(�1-� )�D
i , k
mxy i,k =� ��
8l�2
1
- 1
Rys. 9.5
Gwiazdy różnicowe dla sił poprzecznych i reakcjach w narożach można znalez
ć
w literaturze np. skrypt Z. Waszczyszyna i M. Radwańskiej i podręcznik Z.
Kączkowskiego.
Przykład liczbowy
Przykład 9.1: Płyta kwadratowa przegubowo podparta obciążona równomiernie.
Ze względu na symetrię pokazano na rys.9.6 odpowiednią numerację punktów.
Dla skrócenia zapisu przyjęto 1-indeksowe oznaczenie węzłów. Węzły zewnętrzne
7
oznaczono indeksami ujemnymi, aby uwzględnić warunek brzegowy 9.13. Zgodnie z
9.12 i 9.10.b piszemy równania dla kolejnych węzłów przy podziale boku na 4 odcinki
o długości l� = a /4, por. Rys. 9.6.
- 3 0 - 3 - 2 - 3 0 - 3
x
0 0 0 0 0 0 0
- 3 0 3 2 3 - 3
- 2 0 2 1 2 - 2
- 3 0 3 2 3 - 3
0 0 0
- 3 0 - 3 - 2 - 3 0 - 3
y
Rys. 9.6
pkt.1:
pkt.2.:
pkt.3:
Z warunku (9.13) podparcia przegubowego wynika:
w3 =� -�w
-�3
i po podstawieniu i uporządkowaniu otrzymujemy układ równań:
20w1 -� 32w2 +� 8w3 =� A
-� 8w1 +� 24w -�16w3 =� A
2
2w1 -�16w +� 20w3 =� A
2
8
4
p0a
gdzie: A =� .
256D
Rozwiązaniem tego równania są ugięcia:
w1 =� 1.03125A, w =� 0.75A, w3 =� 0.546875A
2
4
p0a
i po podstawieniu za A =� otrzymujemy
256D
4
p0a
.
w1 =� w =� 0.0040283
max
d
W porównaniu z tablica 8.1 (ćwiczenie 8) błąd wynosi 0,6%.
Maksymalny moment zginający dla =� 0.3 wynosi zgodnie z (9.16):
D D
2
(�mx )� =� [�(�2 +� 2)�w1 -� (�2 +� 2)�w ]� =� 0.73125A =� 0.045703p0a
,
1 2
l�2 l�2
0.04579p0a2
co w porównaniu z wartością daje błąd ok. 5.6%.
l� =� a /8
Dla siatki różnicowej , por. Rys. 9.7, otrzymujemy układ dziesięciu równań
wykorzystując symetrię układu. Przykładowo dla węzła1i6 otrzymujemy:
9
- 1 0 - 8 - 7 - 4
- 1 0 1 0 8 7 4
- 8 9 6 3 6
- 7 6 5 2 5 6
- 4 4 3 2 1 2 3 4
5 2 5
9 6 3 6
3 5 4 7
Rys.9.7
Pełny układ równań różnicowych możemy napisać w postaci macierzowej:
20 -� 32 4 0 8 0 0 0 0 0 w1 0.0625A
�� ł� �� ł� �� ł�
ę�-� 8 25 -� 8 -�1 -�16 6 0 0 0 0 ś� ę� ś� ę�0.0625Aś�
w2
ę� ś� ę� ś� ę� ś�
ę� ś� ę� ś� ę� ś�
1 -� 8 20 -� 8 4 -�16 4 0 2 0 w3 0.0625A
ę� ś� ę�
0 1 -� 8 19 0 4 -�16 2 0 0 w4 ś� ę�0.0625Aś�
ę� ś� ę� ś� ę� ś�
ę� ś� ę�
2 -�16 4 0 22 -�16 2 0 2 0 w5 ś� ę�0.0625Aś�
=�
ę� ś� ę� ś� ę� ś�
.
0 3 -� 8 2 -� 8 23 -� 8 3 -� 8 0 w6 ś� ę�0.0625Aś�
ę� ś� ę�
ę� ś� ę�
0 0 2 -� 8 1 8 20 -� 8 2 1 w7 ś� ę�0.0625Aś�
ę� ś� ę� ś� ę� ś�
0 0 0 1 0 3 -� 8 21 -� 8 -� 8ś� w8 ś� ę�0.0625Aś�
ę� ę�
ę� ś� ę�
0 0 2 0 2 -�16 4 -�16 20 2 w9 ś� ę�0.0625Aś�
ę� ś� ę� ś� ę� ś�
ę� 0 0 0 0 0 0 2 -�16 2 18 ś� ę�w10 ś� ę�0.0625Aś�
�� �� �� �� �� ��
Obliczone ugięcia w węzłach siatki wynoszą:
w =1.038018 A, w = 0.965236 A, w = 0.751953 A
1 2 3
w = 0.416468 A, w = 0.897733 A, w = 0.699773 A
4 5 6
w = 0.387912 A, w = 0.303741 A, w = 0.546421 A
7 8 9
10
w = 0.169650 A.
10
4
p0a
Błąd z jakim określono przy tej siatce w =� w1 =� 1.038018A =� 0.0040547
wynosi
max
D
(�mx )�
0.13% natomiast moment określono z dokładnością 0.987 %.
1
Na rysunku 9.8 przedstawiono wykresy ugięcia i momentów.
Rys.9.8
Literatura:
1. Z. Waszczyszyn,M. Radwańska:Ustroje powierzchniowe
2. Z.Kączkowski: Płyty obliczenia statyczne
11