Leszek Chodor, 26 października 2014
20.10. 2025 scalenie kilku odrębnych artykułów,
26.10. 2025- nadal naprawa po poważnej awarii portalu.
W przypadku nieczytelnych treści, proszę powiadomić: leszek@chodor.co
W ciągu ostatnich 24 godzin z artykułu korzystało 5 Czytelników
W literaturze można znaleźć różne definicje rodzajów zginania: czyste, proste, poprzeczne, ukośne. Definicje te uporządkowano w sposób spójny. Belka jest elementem prętowym w najprostszym przypadku zginanym, a w ogólności również rozciąganym – ściskanym, ścinanym i skręcanym swobodnie oraz nieswobodnie. Element skończony jest jednowymiarowy, w którym wymiar dominujący (długość pręta L) jest o rząd mniejszy od wymiarów przekroju poprzecznego (h i b). W zależności od rodzaju wytężenia i innych założeń definiujemy kilka elementów belkowych (prętowych): Bernoulli, Timoshenko, Własowa. Uwzględnienie naprężeń stycznych w elemencie Timoschenko pozwala za belki traktować średnio-krępe pręty zginane o smukłości $4 \le L/h < 10 (8)$.
Nowoczesne programy inżynierskie mają zaimplementowany element Timoszenko-Własow. Element jest używany nie tylko w programach naukowych, jak Abaqus, Ansys, COMSOL Multiphysics, ale także inżyniwrskich: SAP2000. OpenSees, SolidWorks, ETABS, a nawet Mathcad, LTBEAM i Consteel dla belek krępych.
W artykule przedstawiono podstawy teorii Bernoulli-Euler oraz Timoshenko -Ehrenfest. Pracę opatrzono przykładami rachunkowymi, dotyczącymi nieliniowości geometrycznych pryzmatycznej ściskanej i zginanej belki Bernoulli oraz rozwiązania statycznego belki Timoshenko na podłożu sprężystym.
Belka i zginanie
Belka
W konstrukcjach budowlanych belka odgrywa kluczową rolę w zapewnieniu wytrzymałości i stateczności budynków oraz innych obiektów. W klasycznej definicji belki, przyjmuje się, że belka jest prętem , którego długość $L$ jest znacznie większa od największego wymiaru poprzecznego wysokości $h$ lub szerokości $b$:
$$ \begin{equation} \lambda_b=\cfrac{L}{max [h, b]} \le \lambda_{b, lim} \label{1} \end{equation} $$
gdzie smukłość graniczną $\lambda_{b,lim}$ w normach i zaleceniach projektowania różnych rodzaj w konstrukcji zaleca się przyjmować następująco:
10 w klasycznej teorii sprężystości;
8 w technicznej teorii belek;
3 w konstrukcjach żelbetowych;
4 w konstrukcjach kompozytowych.
Powyżej podane graniczne smukłości dotyczą analiz prowadzonych klasyczną, techniczną teorią belkową Euler -Bernoulli.
W pracy [1] napodsrawie symulacji numerycznych pokazano że. smukłości graniczne $\lambda_{b.lim} przy zastosowaniu różnych teorii belkowych, różniących się stopniem uwzględnienia odkształcalności postaciowej, smukłości graniczne wynoszą:
4 dla teorii belki Euler -Bernoulli,
0,9 dla teorii belki Timoshenko -Ehrenfest,
0,3 dla teorii belki Reddy – Bickford
0,1 przy zastosowaniu teorii wyższego rzędu (7. lub 5. rzędu).
Własne badania wskazują, że graniczne smukłości belek, przy których uzyskuje się dokładność do ok 3%, zwykle wystraczającą w zastosowaniach inżynierskich, wynoszą:
8 dla teorii belki Euler – Bernoulli,
4 dla teorii belki Timoshenko – Ehrenfest,
2 dla teorii belki Reddy – Bickford,
1 przy zastosowaniu teorii wyższego rzędu (7. lub 5. rzędu).
Graniczne smukłości w niewielkim stopni zależą od materiału belki, a przede wszystkim zależą od typu teorii belkowej. Poniżej tych smukłości zaleca się w praktyce stosować elementy dwuwymiarowe ( tarczowe, powłokowe).
W niniejszym artykule przedstawiono dwie teorie stosowane w technice, a mianowicie: Euler -Bernoulli oraz Timoshenko-Ehrenfest.
Zginanie
Zginanie jest definiowane następująco:
- wg słownika języka polskiego [2] : nadawanie czemuś trwałej krzywizny, kabłąkowatego, półkolistego kształtu,
- wg innego słownika języka polskiego [3]: jeśli coś zostało zgięte lub jeśli zgięło się, to gnąc się zmieniło chwilowo lub na stałe swój kształt,
- wg [4] stan deformacji, przy którym prosty w stanie niezdeformowanym pręt, po deformacji jest zakrzywiony (wykazuje różną od zera krzywiznę). Zginanie jest dominującym sposobem pracy prętów, nazywanych belkami.
W artykule przedstawiono podstawowe informacje o zginaniu prętów, które są wykorzystywane w praktycznym projektowaniu belek wykonanych z różnych materiałów konstrukcyjnych, w tym stali, betonu drewna:
- belki stalowe w klasycznym podejściu normowym
- belki żelbetowe w podejściu normowym,
- belki drewniane w ujęciu normowym.
Każda z tych belek ma swoje unikalne właściwości i zastosowania.
Definicje rodzajów zginania pręta pryzmatycznego wynikają z rozwiązań teorii sprężystości, np [5] i są związane z dokładnością wypełnienia warunków statycznych, kinematycznych oraz fizycznych tej teorii.
Definicje przytaczamy dla przypadku zginania pręta pryzmatycznego, to znaczy takiego, którego podstawy (ścianki poprzeczne) leżą w płaszczyznach równoległych do siebie, a ścianki boczne są powierzchniami, zbudowanymi z pęku prostych , łączących odpowiadające punkty obu podstaw. W ogólności pręt pryzmatyczny nie musi mieć prostoliniowej osi lub stałego przekroju. Przypadek pręta zakrzywionego lub pręta o zmiennej sztywności po długości należy rozpatrywać odrębnie, bo występują tam dodatkowe zjawiska nie uwzględnione w przypadku elementarnym.
Zginanie czyste
Przypadek czystego zginania zdefiniowano na rys. 1. Charakteryzuje się on tym, że ścianki czołowe pręta obciążone są bryłami obciążeń, a nie siłami przekrojowymi o gęstości:
$$ \begin{equation} q_{vz} = k \cdot z \label{2} \end{equation} $$
gdzie k jest dowolną stałą, z współrzędną pionową punktów przekroju o wartości z=0 w osi pręta.
Przy takim obciążeniu uzyskuje się ścisłe rozwiązanie teorii sprężystości z wykorzystaniem tylko podstawowych równań teorii; statycznych (Naviera), geometrycznych (Cauchyego) oraz fizycznych (Hooka).
Założenie płaskich przekrojów Bernoulli można wyprowadzić z rozwiązania ścisłego jako wniosek W tym sensie założenie Bernoulli nie jest pierwotne (nie jest założeniem w sensie matematycznym), a jedynie wnioskiem z rozwiązania ścisłego [6].
W wyniku rozwiązania ścisłego zagadnienia „czystego” zginania uzyskujemy również związki na stałą matematyczną k występującą w ($\ref{2}$), a mianowicie:
$$ \begin{equation} k = \cfrac{E}{\rho}= E\cdot w” \left( = \cfrac {M}{I_y} \right )\label{3} \end{equation} $$
gdzie:
E- moduł Younga materiału,
$w$ ugięcie osi, $w”$ – druga pochodna funkcji ugięcia.
$\rho$ – promień krzywizny zginanego pręta
W przypadku czystego zginania promień zginania $\rho$ oraz $w”$ są stałe po długości pręta.
W nawiasie okrągłym w ($\ref{3}$) podano zależność na stałą k w przypadku prostego zginania momentem zginającym M przyłożonym z przeciwnymi znakami do końców pręta , którego przekrój ma moment bezwładności $I_y$
Odwrotność promienia krzywizny $\rho(x)$ (rys.4b) jest nazywana krzywizną ugięcia $\kappa$ i jest związana z ugięciem w(x) znaną zależnością z geometrii różniczkowej [7]
$$ \begin{equation} \kappa \stackrel{\rm def}{=} \cfrac{1}{\rho(x)} = \cfrac{ | w”(x) |}{\left [ 1+ (w'(x))^2 \right ]^{3/2}}\label{4} \end{equation} $$
Znak bezwzględnej wartości użyto w celu podkreślenia, że krzywizna jest nieujemna..

Rys.1. Czyste zginanie (opracowano na podstawi[8]
W przypadku czystego zginania ściśle określone naprężenia we wszystkich przekrojach pręta są takie same (nie zależą od x oraz y) i wynoszą:
$$ \begin{equation} \sigma_x= – k \cdot z \label{5} \end{equation} $$
a pozostałe składowe macierzy naprężeń są zerowe : $\sigma_y=\sigma_z= \tau_{xy}-\tau_{xz}= \tau_{yz}=0$
Ściśle wyznaczone odkształcenia wynoszą:
$$ \begin{equation} \varepsilon_x= – \cfrac{k}{E} \cdot z \quad ; \quad \varepsilon_y = \varepsilon_z= \nu \cdot \varepsilon_x \label{6} \end{equation} $$
gdzie $\nu$ – współczynnik Poissone’a materiału
Wszystkie odkształcenia kątowe (postaciowe) są zerowe: $\gamma_{xy}=\gamma_{yz}=\gamma_{xz} =0$
Zginanie proste
Proste zginanie jest takim przypadkiem obciążenie pręta, że w poprzecznym przekroju pręta działa moment zginający tylko w jednej z głównych płaszczyzn [8]. W żadnym przekroju pręta nie działają siły poprzeczne (tnące). Oznacza to taki stan pręta, w którym momenty zginające są stałe po długości pręta, to znaczy obciążenie zewnętrzne działa tylko na podstawy (końce pręta) i może być sprowadzone do zewnętrznych momentów zginających o takiej samej wartości, lecz przeciwnych zwrotach (kierunkach).
Przypadek prostego zginania zdefiniowano na rys.2.

Rys.2 Proste zginanie (opracowano na podstawie [8]
Na podstawy pręta działa moment zginający M, który jest statycznie równoważny, obciążeniom powierzchniowym $q_{vz}$, przedstawionym na rys.1. Obciążenia te związane są zależnością
$$ \begin{equation} q_{vx} = \cfrac{M}{I_y}\cdot z \label{7} \end{equation} $$
Po podstawieniu wyrażenia na stałą k $ (\ref{3})^3$ (w nawiasie) do ($\ref{5}$) oraz ($\ref{6}$) uzyskujemy wyrażenia na naprężenia oraz odkształcenia w postaci
$$ \begin{equation} \sigma_x= – \cfrac {M}{I_y} \cdot z \label{8} \end{equation} $$
$$ \begin{equation} \varepsilon_x = – \cfrac{M}{EI_y} \cdot z \quad ; \quad \varepsilon_y = \varepsilon_z= \nu \cdot \varepsilon_x\label{9} \end{equation} $$
Pozostałe składowe macierzy naprężeń oraz odkształceń są zerowe : $\sigma_y=\sigma_z= \tau_{xy}-\tau_{xz}= \tau_{yz}=0$, $\gamma{xy} =\gamma_{yz}=\gamma_{xz} =0$
Poprzeczne zginanie
W konstrukcjach budowlanych praktycznie nie mamy do czynienia z czystym lub prostym zginaniem, lecz wszystkie przypadki zginania to zginanie poprzecze, to znaczy współdziałanie zginania i ścinania pręta.
Poprzeczne zginanie (rys.3) jest to działanie takiego obciążenie na nieswobodny (podparty) pręt pryzmatyczny, o przekroju poprzecznym symetrycznym względem osi z, które jest przyłożone symetrycznie względem osi z, na powierzchni A pobocznicy i ściankach poprzecznych i ma gęstość $q_v$ = – q(x). Znak minus oznacza, że obciążenie ma znak przeciwny do kierunku osi z (jest skierowane pionowo w dół).
Znalezienie ścisłego rozwiązanie zagadnienia poprzecznego zginania jest w ogólnym przypadku bardzo trudne, ponieważ zależy od kształtu przekroju poprzecznego,. W przypadku najprostszego przekroju prostokątnego nie są spełnione równania nierozdzielności odkształceń. Znane rozwiązania techniczne (wytrzymałości materiałów) są tylko pewnym przybliżeniem rozwiązania dokładnego. Z zasady de Saint Venanta ( o lokalnym efekcie działania obciążeń i statyczne równoważności) wynika, że rozwiązania czystego i prostego zginania mogą być stosowane do innych obciążeń pręta, w tym w analizie zginania poprzecznego. Należy jednak pamiętać, że nie będą to rozwiązania ścisłe, a tylko ich aproksymacja
Na czole belki poprzecznie zginanej z rys.3 pokazano siły przekrojowe wyznaczone zgodnie z techniczną teorią zginania.

Rys.3. Poprzeczne zginanie [8] (uzupełnione o obciążenie siłę osiową w celu wskazania, że belka często oprócz zginania i ścinania jest też rozciągana/ściskana
Belka klasyczna Bernoulli-Euler jest opracowana na bazie teorii czydstego lub prostego zginania. Natomiast teoria poprzecznego zginania jest uogólnieniem teorii Bernoulli-Euler w spsob pokazany w rozdziale.
Ukośne zginanie
Ukośne zginanie jest przypadkiem obciążenia pręta jednocześnie zginaniem w płaszczyźnie (x,z) (na rys.3 $M_y$) jak i w płaszczyźnie (x-,y) momentem $M_z$. Zgodnie z zasadą superpozycji w stanie sprężystym rozwiązanie w złożonym stanie jest sumą rozwiązań od poszczególnych prostych obciążeń.
Obciążenie momentem $M_x$ w płaszczyźnie (y,z) nazywa się skręcaniem i jest jakościowo innym od zginania przypadkiem obciążenia i wytężenia, który omówiono w artykule Skręcanie pretów. Prety cienkościenne.
Belka Bernoulli-Euler
Rys historyczny i założenia klasycznej teorii belkowej
Teoria belki Bernoulli powstała ok 1750 roku dzięki wkładowi Leonarda Euler i Daniela Bernoulli. Bernoulli podał wzór na energię odkształcenia przy zginaniu belki, z którego Euler wyprowadził i rozwiązał równanie różniczkowe, przy czym skorzystali również osiągnięciach Jacoba Bernoulli. Wcześniej problem zginania belki podejmowali Leonardo da Vinci i Galileusz .
Dwa kluczowe założenia teorii belki Euler -Bernoulli to:
1) Materiał jest liniowo sprężysty zgodnie z prawem Hooke’a
2) Przekrój płaski przed odkształceniem, pozostaje płaski i prostopadły do odkształconej osi pręta po jego obciążeniu.
Założenie płaskich przekrojów jest nazywane hipotezą Naviera lub zasadą zesztywnienia.
Założenie płaskich przekrojów jest również stosowane do zagadnień dużych sił ściskających (zagadnienia niestateczności giętnej ), jak i dużych obrotów (zagadnienia giętno-skrętne). Założenie płaskich przekrojów i wynikający z niego rozkład odkształceń zginania po wysokości przekroju zilustrowano na rys. 4 .

Rys. 4 Belka Bernoulli-Euler : a) belka przed odkształceniem, b) belka po odkształceniu, c) rozkład odkształceń po wysokości przekroju zginanego
W modelu zakłada się, że na odkształcenia przekroju $\varepsilon$ mają wpływ jedynie naprężenia normalne wywołane momentem zginającym $M$, a pomija się naprężenie styczne. Konsekwencją takiego uproszczenia jest zasada płaskich przekrojów, a także liniowy rozkład odkształceń po wysokości przekroju wg zależności:
$$ \begin{equation} \varepsilon (x,z) =\Theta (x) \cdot z \\ \Theta= \cfrac{\partial w(x)}{\partial x} \label{10} \end{equation} $$
gdzie:
$\Theta(x)$ – kąt ugięcia osi belki w punkcie o współrzędnej x,
$ w(x)$ – linia ugięcia osi obojętnej belki O’
Ugięcie $w(x)$ jest przemieszczeniem osi belki w kierunku osi z (na rys 4, ugięcie jest ujemne).
Opis modelu belki
Na rys. 5 zebrano ilustracje do wyprowadzenia równań rządzących zagadnieniem poprzecznego zginania belki Bernoulli-Euler .

Rys. 5 Równowaga elementu dx belki: a) równowaga elementu dx zginanego , b) równowaga elemntu dx przy rozwarstwianiu, c) dodatnie ugięcie $w$ i kąt obrotu przekroju $\Theta$ (opracowano na podstawie [9]
Rys. 5b nie jest ilustracją do oryginalnych analiz teorii Bernoulli-Euler , ale uzupełnia model oryginalny o analizę rozwarstwiania belki na skutek działania sił poprzecznych i w tym sensie jest elementem zmodyfikowanej teorii zginania belki Bernoulli-Euler .
Równania równowagi
Równania równowagi belki uzyskuje się, rozpatrując równowagę nieskończenie małego elementu belki dx , pokazanego na rys. 5a. Z warunku sumy rzutów na oś pionową ( $\sum Z=0$) oraz sumy momentów ($\sum M=0$) otrzymujemy
$$ \begin{equation} q_z =\cfrac{dV}{dx} \quad ; \quad V= \cfrac{dM}{dx} +(…) \label{11} \end{equation} $$
W $(\ref{11})$ pominięto człony wyższego rzędu (pochodne po $dx^2$ i wyżej) , oznaczone (…).
Moment zginający $M = M_y $ uzyskuje się poprzez Całkowanie naprężeń osiowych $\sigma= \sigma_x$ po całym przekroju poprzecznym $A$ :
$$ \begin{equation} M= \int \limits_A – \sigma \cdot z dA \label{12} \end{equation} $$
gdzie pojawia się znak minus, ponieważ naprężenia ściskające (ujemne) są w dodatnim obszarze osi z .
Ugięcia $w(x)$ i momenty zginające $M(x)$ są dodatnie gdy będą miały takie same zwroty i w związku z tym w technice umówiono się znakować w i M na podstawie wykresów w taki sposób, że znak momentu zginającego lub ugięcia jest dodatni, jeśli ich wykresy leżą po stronie spodu pręta. Powyższa umowa o znakowaniu prowadzi do pojęcia „spodu” belki, który na wykresach oznaczany jest najczęściej podkreśleniem linią przerywaną tej strony belki.
Równania fizyczne
W przypadku czystego zginania belki liniowo sprężystej obowiązuje jednoosiowe prawo Hooke’a wiążące naprężenia $\sigma=\sigma_x$ z odkształceniami $\varepsilon = \varepsilon_x$ w płaskim stanie naprężenia:
$$ \begin{equation} \sigma = E \varepsilon \label{13} \end{equation} $$
W przypadku belki, w której odkształcenia w kierunku osi y są ograniczone ( np. podciąg w płycie stropowej) z warunku zerowania się odkształceń poprzecznych $\varepsilon_y$ w dwuwymiarowym prawie Hooke’a otrzymujemy:
$$ \begin{equation} \varepsilon_y = \cfrac{\sigma_y}{E} – \nu \cdot \cfrac{\sigma_x}{E} =0 \to \sigma_y = \nu \cdot \sigma_x \label{14} \end{equation} $$
i w tym przypadku bardziej odpowiednie jest stosowanie prawa fizycznego w postaci
$$ \begin{equation} \varepsilon_x =\cfrac{\sigma_x}{E} – \nu \cfrac{\sigma_y}{E} \to \sigma_x = \cfrac {E}{1-\nu^2} \cdot \varepsilon_x \label{15} \end{equation} $$
czyli jak dla płaskiego stanu odkształceń.
Wszystkie poniższe wyprowadzenia i wyniki oparte są na prawie materialnym ($\ref{15}$). Prawo Hook’a w postaci ($\ref{15}$) można w prosty sposób wprowadzić zastępując prosty moduł Younga E przez moduł zastępczy $\overline E$:
$$ \begin{equation} \overline E= \cfrac {E}{1-\nu^2}\label{16} \end{equation} $$
Równania geometryczne
Odkształcenia zginające prowadzą do skracania (powyżej osi obojętnej) lub wydłużania (poniżej osi) włókien w przekroju poprzecznym. Z rys 5d, odczytujemy, zę odkształcenie włókien $\varepsilon = \varepsilon_x$ wynosi
$$ \begin{equation}\varepsilon = \cfrac {du}{dx} \label{17} \end{equation} $$
Jednocześnie przemieszczenie osiowe u jest powiązane z obrotem przekroju poprzecznego. Rozważmy nieskończenie mały obrót w kierunku przeciwnym do ruchu wskazówek zegara $d\Theta$ nieskończenie małego elementu belki $dx$ na rys. 5d . Przy okazji, należy zauważyć, że $d\Theta$ jest równe krzywiźnie $\kappa =1/\rho$. Z założenia płaskich przekrojów wynika, że odkształcenie każdego włókna w przekroju poprzecznym zmienia długość proporcjonalnie do swojej odległości od osi obojętnej. Stopień skrócenia lub wydłużenia zależy od obrotu przekroju poprzecznego. Rozważania geometryczne przeprowadzone na rys 5d pokazują, że skrócenie lub wydłużenie, tj. przemieszczenie osiowe, każdego nieskończenie krótkiego włókna wynosi
$$ \begin{equation} du = – d \Theta \cdot z \label{18} \end{equation} $$
Obrót $\Theta$ jest powiązany z ugięciem, jak pokazano na rys. rys 5c. Z zależności geometrycznych i założenia małych obrotów $\Theta$ mamy:
$$ \begin{equation} tg \Theta = \cfrac{dw}{dx} \approx \Theta \to \Theta = arctg \left( \cfrac{dw}{dx} \right ) \label{19} \end{equation} $$
Połączenie przedstawionych wyżej zależności daje następujące równanie kinematyczne dla belki
$$ \begin{equation} \varepsilon = \cfrac{du}{dx} = -\cfrac{d \Theta}{dx} \cdot z = – \cfrac{d^2 w}{dx^2}\label{20} \end{equation} $$
Uzyskaliśmy przybliżenie dokładnej krzywizny belki ($\ref{4}$) stosowane w technicznej teorii zginania:
$$ \begin{equation} \kappa = \cfrac{d \Theta}{dx} \approx \cfrac{w”}{1+w”}\label{21} \end{equation} $$
co wynika z założenie małych przemieszczeń, w tym $dw/dx$.
Podstawowe zależności technicznej teorii belek
W technicznej teorii zginania porzecznego stosuje się zależności uzyskane problemu czystego zginania ( mimo, że porzeczne zginanie jest innym przypadkiem niż zginanie czyste ) w postaci :
- zlinearyzowanej (uzyskanej po odrzuceniu członów wyższego rzędu) zależności krzywizny od ugięcia belki
$$ \begin{equation} \kappa = \cfrac{1}{\rho(x)} \approx | w”(x)|\label{22} \end{equation} $$
- zlinearyzowaną zależność krzywizny belki od momentu zginającego M(x)
$$ \begin{equation} \kappa \approx \cfrac{| M(x)|}{EI_y}\label{23} \end{equation} $$
- i stąd zależność różniczkową na linii ugięcia belki
$$ \begin{equation} w”(x) = – \cfrac{ M(x)}{EI_y}\label{24} \end{equation} $$
- wyrażenia na naprężenia normalne przy zginaniu z połączenia ($\ref{3}$) i ($\ref{5}$), czyli ($\ref{8}) w postaci
$$ \begin{equation} \sigma_x = \cfrac{ M(x)}{EI_y}\cdot z \label{25} \end{equation} $$
a po wprowadzeniu pojęcia wskaźnika wytrzymałości przekroju
$$ \begin{equation} W_y \stackrel{\rm def}{=} \cfrac{I_y}{max \, z} \label{26} \end{equation} $$
gdzie: max z = maksymalna odległość włókna przekroju od osi obojętnej
Z ($\ref{25}$) mamy wzór na maksymalne naprężenie normalne w przekroju belki
$$ \begin{equation} max \, \sigma_x = \cfrac{ M(x)}{W_y} \label{27} \end{equation} $$
stanowiący jeden z najczęściej stosowanych wzorów wytrzymałości materiałów.
Uogólnienie teorii belki Bernoulli- Euler o naprężenia styczne. Wzór Żurawskiego
Teoria belki Euler -Bernoulli nie modeluje naprężeń stycznych w przekroju. Teoria zakłada, że płaskie przekroje pozostają płaskie i prostopadłe do osi obojętnej. Innymi słowy, jedynym występującym odkształceniem jest osiowe skrócenie lub wydłużenie włókien w przekroju poprzecznym. Brak ścinania wystąpi w przypadku zginania czystego lub prostego. i na takie przypadki została opracowana teoria Bernoulli-Euler .
Natomiast w przypadku zginania poprzecznego siła poprzeczna, wywołująca ścinanie wystąpi praktycznie w każdym przypadku (oprócz zginania momentem zginającym stałym na długości belki, bo siła poprzeczna jest pochodną momentu zginającego). Tę anomalię teorii Bernoulli poprawia rozwiązanie belki Timoshenko.
W tym rozdziale poprzestaniemy na analizie naprężeń ścinających $\tau$ „sił rozwarstwiający $q_s$” w dowolnym punkcie przekroju poprzecznego. Rozważmy nieskończenie krótki element belki , pokazany na na rys. 5b. przecięty na wysokości z siła rozwarstwiąjąca $q_s$ w tym miejscu.
Siła rozwarstwiająca $q_s$, to siła na jednostkę długości belki, która zapewnia równowagę z naprężeniami normalnymi $\sigma$. Naprężenia $\sigma$ są większe po jednej stronie belki od naprężeń po drugiej ze względu na przyorst wartości momentu zginającego $dM$:
$$ \begin{equation} q_s \cdot dx = \int\limits_A d \sigma dA = \int\limits _{A_s} \cfrac{dM}{I_y }\cdot z dA \label{28} \end{equation} $$
gdzie: $A_s$ jest pole przekroju poprzecznego powyżej przecięcia
Biorąc pod uwagę, że $V = \cfrac{dM}{dx}$ , otrzymujemy
$$ \begin{equation} q_s =\cfrac{V} {I_y} \cdot \overline S_y \label{29} \end{equation} $$
gdzie moment statyczny odciętej części przekroju wynosi
$$ \begin{equation} \overline S_y =\int \limits_{A_s} z dA\label{30} \end{equation} $$
W praktyce w celu wyznaczenia $\overline S_y$ przekrój poprzeczny jest często dyskretyzowany na kilka części o polu powierzchni $A_i$ i odległości $z_i$ od osi obojętnej do środka ciężkości części. Następnie wylicza się
$$ \begin{equation} \overline S_y =\sum \limits_{i=1}^N A_i \cdot z_i \label{31} \end{equation} $$
przy czym każdy z N części przekroju należy do $A_s$
Poziome naprężenia styczne, rozwarstwiające odcięty przekrój wynoszą
$$ \begin{equation} \tau = \cfrac{V\cdot \overline S_y}{I_y \cdot t}\label{32} \end{equation} $$
gdzie t jest szerokością przekroju w miejscu przecięcia.
Należy zauważyć, że równania ($\ref{29}$) i ($\ref{32}$) są tylko przybliżeniami uzyskanymi w drodze wielokrotnych uproszczeń . Największą dokładność oszacowania uzyskuje się dla belki pryzmatycznej o prostej osi z równomiernym przekrojem poprzecznym i wówczas, gdy siła ścinająca działa równolegle do osi z i krawędzie przekroju są równoległe do tej osi, a przekrój jest wysoki i cienki . Takie warunki spełnia prostokąt, ale już w trójkącie lub przekroju kołowym. Timoshenko i Goodier [10] stwierdzili, że dokładność wzoru Żurawskiego ($\ref{32}$) jest całkiem dobra nawet dla przekrojów kołowych i trójkątnych.
Dodatkowe, pomocnicze wskazówki do wyznaczeni naprężeń stycznych w przekroju belki wynikają z analiz równań na warunki brzegowe na pobocznicy pręta praz ze ścisłych warunków równoważności wewnątrz przekroju. Wynika z nich, że :
- wektor naprężeń stycznych $\\vec{tau}$ musi być styczny do konturu przekroju,
- suma momentów od napreżeń stycznych $\tau_{xy}$ i $\tau_{xz}$ względem środka ciężkości musi być równa zeru, skąd wynika, że dla przekroju symetrycznego względem osi pionowej $z$ rozkład tych naprężeń względem tej osi musi być również symetryczny.
- zwrot naprężenia $\tau_{xy}$ i $\tau_{xz}$ powinien być zgodny ze zwrotem siły poprzecznej $V_z(x)$ .
Równanie różniczkowe zginania poprzecznego belki Bernoulli
Na podstawie analiz z poprzednich rozdziałów można już wyprowadzić równanie różniczkowe rządzące zagadnieniem zginania belki Bernoulli następujący sposób:
$$ \begin{equation} q_z = \cfrac{V}{dx } = \cfrac{d^2 M}{dx^2} = – \cfrac{d^2}{dx^2} \int \sigma \cdot z dA = – \cfrac{d^2}{dx^2} \int E \cdot \varepsilon \cdot z dA = – \cfrac{d^2}{dx^2} \int_A E \cdot \cfrac{d^2 w}{dx^2 }\cdot z^2 dA = EI_y \cfrac{d^4 w }{dx^4} \label{33} \end{equation} $$
gdzie wprowadzono definicję momentu bezładności przekroju
$$ \begin{equation}I_y \stackrel{\rm def}{=} \int_A z^2 dA \label{34} \end{equation} $$
W dobie komputeryzacji rozwiązywanie równania różniczkowego dla zginania belki jest rzadko stosowane w praktyce inżynierskiej, warto przestudiować jego rozwiązanie w prostych przypadkach odniesienia. W szczególności rozwiązanie równania różniczkowego stanowi punkt wyjścia do wyboru funkcji kształtu w metodzie elementów skończonych. Funkcje kształtu są często przybliżone, podczas gdy rozwiązanie równania różniczkowego ujawnia dokładny kształt, gdy element ulega odkształceniu. Ogólne rozwiązanie równania różniczkowego ujawnia, czy funkcje kształtu elementów skończonych są dokładne, czy nie. W przypadku belek ogólne rozwiązanie równania różniczkowego uzyskuje się przez całkowanie czterokrotne, zakładając, że obciążenie $q_z$ jest równomiernie rozłożone wzdłuż belki w postaci paraboli czwartego stopnia
$$ \begin{equation} w(x) = \cfrac{1}{24} \cfrac{q_z}{EI_y}\cdot x^4 +C_1 x^3 +c_2 x^2 +C_3 x +C_4\label{36} \end{equation} $$
Stałe całkowania wyznacza się z warunków brzegowych po dwa na każdym końcu belki lub w innych miejscach pręta, przy czym mogą one dotyczyć ugięcia w lub jego pochodnych:
$ \cfrac {dw}{dx}=\Theta$,
$ EI_y \cfrac{d^2w}{dx^2}=M$ ,
$EI_y \cfrac{d^3 w}{dx^3} = V$.
Uogólnienie teorii belki Bernoulli- Euler o ściskanie. Nieliniowości geometryczne
Teoria zginania Bernoulli-Euler może być rozszerzona na przypadek ściskania belki. Na rys.6b pokazano belkę z obciążeniem poprzecznym i ściskającym. Model elementu belkowego [e] zginanego i ściskanego przedstawiono na rys. 6 a. Na rys, 6c pokazano słup wspornikowy zginany i ściskany.

Rys. 6 Belka-słup Bernoulli zginana poprzecznie i ściskana: a) model elementu [e]; b) belka-slup poziomy, c) słup-belka pionowy
Ścisła macierz sztywności elementu belkowego ściskanego
Nieliniową macierz sztywności ściskanego elementu prętowego [e], pokazanego na rys.6a, wyznaczymy w drodze rozwiązania równania różniczkowego ($\ref{33}$) belki-słupa bez obciążenia poprzecznego $q_z = 0 $, ale uzupełnionego o wyraz wynikający ze ściskania belki w postaci:
$$ \begin{equation} 0 = EI_y \cfrac{d^4 w }{dx^4} + P \cdot \cfrac{dw}{dx}\label{37} \end{equation} $$
Wprowadzamy zmienne:
$\xi=\cfrac{x}{L}$ – bezwymiarowa współrzędna osi belki,
$ \lambda^2=\cfrac{PL^2}{EI} =\pi^2 \alpha_{cr}$,
czyli $\lambda=\pi \sqrt{\alpha_{cr}}$ — miara obciążenia ściskającego,
gdzie:
$\alpha_{cr}= \cfrac{P}{P_{cr}}$,
$P_{cr}=\pi^2 EI/L^2$ – obciążenie krytyczne Euler
W granicznym stanie krytycznym $\alpha_{cr}=1$ i stopień wytężenia siła ściskającą $\lambda=\pi$.
Równanie pręta ($\ref{37}$) można teraz przepisać do postaci
$$\begin{equation} \cfrac{d^4 w}{d \xi^2}+ \lambda^2 \cfrac{d^2 w}{d \xi^2}=0 \label {38} \end {equation} $$
Rozwiązanie równania ($\ref{38}$) jest funkcją:
$$\begin{equation} w(\xi) = C_1 + C_2 \lambda \xi + C_3 \cos{( \lambda \xi)} + C_4 \sin{( \lambda \xi)} \label {39} \end {equation} $$
Po wykonaniu różniczkowania funkcji ugięcia uzyskamy funkcję obrotów
$$ \begin{equation} \cfrac{ d w ( \xi )} { d \xi } = \lambda [ C_2 + C_3 \sin{ ( \lambda \xi) } + C_4 \cos { ( \lambda \xi )} ] \label {40} \end {equation} $$
Warunki brzegowe wynoszą:
$$\begin{equation} w(0)= w_1 , \quad w (1)=w_2, \quad \cfrac{d w(0)}{d \xi} = \varphi_1 L, \quad \cfrac{d w(1)}{d \xi} = \varphi_2 L \label {41} \end {equation} $$
Warunki brzegowe ($\ref{41}$) można zapisać w formie macierzowej $ [D] \cdot \begin{vmatrix} C1 \\ C_2 \\ C_3 \\ C_4 \end{vmatrix} = \begin {vmatrix} w_1 \\ \varphi_1 L \\ w_2 \\ \varphi_2 L \end {vmatrix}$, gdzie [D]= $\begin {bmatrix} 1 & 0 &1 &0 \\ 0 & \lambda & 0 & \lambda \\ 1 & \lambda & \cos {\lambda} & {\sin \lambda} \\ 0 & \lambda & -\lambda \sin{ \lambda} & \lambda \cos{ \lambda} \end {bmatrix} $
Po odwróceniu macierzy [D] , otrzymamy
$ [D]^{-1} =\cfrac {1} {\Theta} \begin {bmatrix}
\lambda \beta & \gamma & \lambda (1- \cos{ \lambda}) & – ( \lambda – \sin{ \lambda} ) \\
\lambda \sin{ \lambda} & 1- \cos{\lambda} & -\lambda \sin{\lambda} & 1- \cos { \lambda} \\
\lambda (1- \cos{ \lambda}) & -\gamma & -\lambda(1-\cos{ \lambda}) & \lambda -\sin {\lambda} \\
– \lambda \sin{\lambda} & \beta & \lambda \sin { \lambda} & -(1-\cos{ \lambda})\\
\end {bmatrix}$
gdzie:
$\Theta = 2 \lambda (1-\cos {\lambda})- \lambda^2 \sin{\lambda}$,
$ \beta= 1- \cos{\lambda } -\sin{\lambda }$,
$\gamma=\lambda\cos{\lambda }-\sin{\lambda }$.
Stałe całkowania wynoszą
$\begin {vmatrix} C_1 \\ C_2 \\ C_3 \\ C_4 \end {vmatrix} = [D]^{-1} \cdot \begin {vmatrix} w_1 \\ \varphi_1 L \\ w_2 \\ \varphi_2 L \end {vmatrix} =
\cfrac{1}{\Theta} \begin {vmatrix}
\lambda (\Sigma w -\varphi_2 L)+ \lambda(\varphi_1 L -\Sigma w ) \cos{\lambda}- (\lambda^2 w_1 -\Delta \varphi L) \sin{\lambda} \\
\Sigma \varphi L \cdot (1-\cos{\lambda}) -\lambda \Delta w \sin {\lambda} \\
\lambda(\varphi_2L -\Delta w ) -\lambda(\varphi_1 L -\Delta w) \cos{\lambda} -\Delta \varphi L \sin {\lambda} \\
-\Delta \varphi L \cdot (1-\cos{\lambda}) -\lambda (\varphi_1 L -\Delta w ) \sin {\lambda} \\
\end {vmatrix}$
Siły przekrojowe uzyskamy z zależności:
$M = EI w^{(ii)} = – \cfrac{EI}{L^2} \cfrac{dw^2}{d\xi}= \lambda^2 \cfrac{EI}{L^2} (C_3 \cos{\lambda \xi} +C_4 \sin{\lambda \xi})$,
$ V = M^{(i)} = – \cfrac{EI}{L^3} \cfrac{dw^3}{d \xi} = \lambda^3 \cfrac{EI}{L^3} ( – C_3 \sin{\lambda \xi} +C_4 \cos{\lambda \xi})$
Wektor sił węzłowych wynosi
$ \begin {vmatrix} M_1 \\ T_1 \\ M_2 \\ T_2 \end {vmatrix} = \begin {vmatrix} M(\xi=0) \\ T(0) \\ M(1)\\ T(1) \end {vmatrix} =
\lambda^2 \cfrac{EI}{L^2} \cdot \begin {bmatrix} 1 & 0 \\ 0 & \cfrac{\lambda}{L} \\ \cos{\lambda} & \sin {\lambda} \\ \cfrac {- \lambda}{ L \sin{\lambda}} & \cfrac{ \lambda}{ L \cos{\lambda}} \end {bmatrix} \cdot \begin {vmatrix} C_3\\ C_4 \end {vmatrix}$
Po przekształceniach powyższych wyrażeń, uzyskamy macierz sztywności nieliniowego geometrycznie elementu ściskanego [e] w postaci:
$$\begin{equation} [k]^{[e]} = \left [ \begin {array} {cccccc}
\cfrac{EA}{L} & 0& 0 & \cfrac{-EA}{L} &0 &0 \\
& \cfrac{12 EI}{L^3} \Lambda_1 & \cfrac{6EI}{L^2} \Lambda_2 & 0 & \cfrac{-12 EI}{L^3} \Lambda_1 & \cfrac{6EI}{L^2} \Lambda_2 \\
& & \cfrac{4EI}{L} \Lambda_3 & 0 & \cfrac{-6EI}{L^2} \Lambda _2& \cfrac{2EI}{L} \Lambda_4 \\
& & & \cfrac{EA}{L} & 0 & 0 & \\
& SYM & & & \cfrac{12EI}{L^3} \Lambda_1 & \cfrac{-6 EI}{L^2} \Lambda_2 \\
& & & & & \cfrac{4 EI}{L} \Lambda_3 \\
\end {array} \right ] \label {42} \end {equation} $$
Macierz ($\ref{42}$) zawiera konwencjonalne współczynniki pręta zginanego zmodyfikowane współczynnikami $\Lambda_i \quad (i=1,..,4)$, które definiują ścisłą, geometrycznie nieliniową macierz sztywności pręta ściskanego, zależną od wielkości siły osiowej N w elemencie, wywołanej siłami P i wynoszą:
$$\begin{equation} \begin {cases}
\Lambda_1= & \alpha_c \cdot \Lambda_2 \\
\Lambda_2= & \cfrac{\alpha^2}{3(1-\alpha_c) } \\
\Lambda_3= & \cfrac{3}{4} \Lambda_2+\cfrac {\alpha_c}{4} \\
\Lambda_4= & \cfrac{3}{2} \Lambda_2 -\cfrac{\alpha_c}{2} \\
\end{cases} \label {43} \end {equation} $$
gdzie:
$\alpha_c = \alpha \cdot ctg{\alpha}$,
$\alpha = \cfrac{\lambda}{2}= \cfrac{\pi}{2} \sqrt{\cfrac{N}{N_{cr}}}$
W przypadku znajomości siły osiowej w elemencie, współczynniki ($\ref{43}$), pozwalają wyznaczyć ścisłą sztywność elementu [e] i w rezultacie nieliniową sztywność systemu konstrukcyjnego. W prostych przypadkach metoda jest użyteczna do obliczeń symbolicznych oraz porównawczych.
W konwencjonalnych algorytmach MES, siły osiowe są poszukiwaną zmienną, więc współczynniki $\Lambda_i$ wywołują nieliniowość i potrzebę iteracyjnego rozwiązywania systemu konstrukcyjnego. Do tego celu użyteczne jest rozwinięcie funkcji $\Lambda_i \quad (i=1,2,3,4)$ w szereg Taylora podług potęg :
$\lambda^2 = \cfrac{NL^2}{EI}$
$$\begin{equation} \Lambda_i \approx \sum \limits_{i=0}^{i=4} C_{1,i} \cdot \left( \lambda^2 \right)^i \label {44} \end{equation} $$
Po pozostawieniu dwóch pierwszych wyrazów rozwinięcia w pobliżu $\alpha = 0$, otrzymamy:
$$\begin{equation} \begin {cases}
\Lambda_1=1- \cfrac {1}{53} \lambda^2\\
\Lambda_2= 1- \cfrac {1}{75} \lambda^2\\
\Lambda_3= 1- \cfrac {1}{66} \lambda^2\\
\Lambda_4= 1+ \cfrac {1}{75} \lambda^2\\
\end{cases} \label {45} \end {equation} $$
Wyrażenia ($\ref{45}$) definiują konwencjonalną macierz geometryczną elementu prętowego.
W zależności od stopnia zachowanych członów w rozwinięciu ($\ref{44}$) mówimy o stopniu analizy MES. Człony liniowe (i=0) nie generują macierzy geometrycznej (jest to stopień 1-szy analizy MES) ; człony (i=1) są uwzględnione w analizie II rzędu, czyli liniowej analizie wyboczeniowej LBA), czasami nazwaną anlizą P-Δ, a współczynniki rozwinięcia, zawierające określają tak zwaną macierz geometryczną ($\ref{45}$). Człony rozwinięcia wyższych rzędów są uwzględniane w analizie trzeciego rzędu (i=2), czwartego (i=4) itd. Współczynniki rozwnięcia ścisłej macierzy geometrycvznej nieliniowości zestawiono w tab.1.
Tab. 1 Współczynniki $C_{1, i}$ rozwinięcia funkcji $\Lambda_i$ w szereg Taylora podług potęg $\lambda^2 = \cfrac{NL^2}{EI}$
W analizach wpływu stopnia nieliniowości teorii na efekty, bardziej użyteczne jest rozwinięcie funkcji szereg podług potęg współczynnika $\alpha = \cfrac{\pi}{2} \sqrt{\cfrac{N}{N_{cr}}}$.
$$\begin{equation} \Lambda_i \approx \sum \limits_{i=0}^{i=4} C_{2,i} \cdot \alpha^i \label {46} \end{equation} $$
Współczynniki rozwnięcia ($\rf{46}$) zestawiono w tab.2.
Tab. 2 Współczynniki $C_{2, i}$ rozwinięcia funkcji w szereg Taylora podług potęg $\alpha$
Ugięcie belki ściskanej przed i w momencie utraty stateczności Eulera
Z zależności ($\ref{38}$) po podstawieniu uzyskanych wyrażeń na stałe całkowania można wyznaczyć funkcję ugięcia pręta pomiędzy węzłami w zależności od przemieszczeń węzłowych:
$$ \begin{equation} w(\xi)=\cfrac{1}{\Theta}
\begin {vmatrix}
1- \cos{\lambda} + \cos{ \lambda \xi} -\cos{[ \lambda(1-\xi)] } – \lambda \sin { [\lambda(1-\xi)] } \\
– \lambda \xi – \lambda (1-\xi) \cos{\lambda}+\xi\cos{ [ \lambda(1-\xi)] } +\sin{\lambda} +\sin{[ \lambda\xi] } \sin{[ \lambda(1-\xi)] } \\
1- \cos{\lambda} -\cos{\lambda \xi} +\cos{[ \lambda(1-\xi)] } -\lambda \xi \sin \lambda \\
– \lambda (1- \xi ) – \lambda \xi \cos{\lambda} +\lambda \cos{\lambda \xi } +\sin{\lambda} – \sin{\lambda \xi } -\sin { [\lambda(1-\xi ) ]} \\
\end {vmatrix}
\cdot \begin {bmatrix} \lambda w_1 & \varphi_1 L & \lambda w_2 & \varphi_2 L\end {bmatrix} \label {47} \end {equation}$$
Z formuły ($\ref{47}$) wynika, że pręt idealnie prosty, ugina się pod wpływem osiowej siły ściskającej w każdym przypadku, gdy węzły doznają przemieszczeń liniowych lub obrotów wstępnych, czyli w każdej konstrukcji rzeczywistej.
W stanie granicznym krytycznym, dla $\lambda= \pi$ mamy:
$$\begin{equation} w(\xi)=\cfrac{1}{2 \pi} \left [ \Sigma w +\Sigma \varphi L (\xi-1/2) + \cfrac{ \pi}{2} [ \Sigma \varphi L -2 \Delta w ] \cos{ ( \pi \xi)}- \Delta \varphi L sin (\pi \xi) \right ] \label {48} \end {equation} $$
Maksymalne ugięcie pręta, czyli strzałka wygięcia wystąpi w miejscu pręta o współrzędnej $\xi_{sup}$
$$\begin{equation} \xi_{sup} = \pm \cfrac{1}{\lambda} \cdot arcos {( \pm \Gamma)} \label {49} \end {equation} $$
gdzie $ \Gamma = \cfrac {1} { \sqrt{ (C_3/C_4)^2 +1}} $.
Wyrażenie uzyskano w drodze porównania do zera drugiej pochodnej funkcji $\ref{37}$), a w rezultacie uzyskano cztery pierwiastki (przy wszystkich kombinacjach znaków). Po odrzuceniu pierwiastków ujemnych jako fizycznie niemożliwych i po podstawieniu dodatniej współrzędnej ($\ref{49}$) do ($\ref{47}$) otrzymamy strzałkę ugięcia pręta w zależności od warunków brzegowych wyrażonych stałymi $C_3 , C_4$.
W przypadku symetrycznych warunków brzegowych strzałka ugięcia wystąpi w środku pręta i wyniesie
$$\begin{equation} f=w (\xi =0,5 ) = \cfrac{1}{2} \left ( \Sigma w -\cfrac { tg(\lambda/4) \cdot \Delta \varphi L }{\lambda}\ \right) \label {50} \end {equation} $$
W granicznym stanie krytycznym: $ \lambda=\pi$ i strzałka ugięcia $(\ref{50}$) osiągnie wartość:
$$\begin{equation} f_{cr} = \cfrac{\Sigma w}{2} – \cfrac{\Delta \varphi L}{2 \pi} \label {51} \end {equation} $$
Dla zerowej siły osiowej $N=0 $ jest $\lambda=0$, więc $f_0 = \lim \limits_{\lambda \to 0}= (\Sigma w) /2 -(\Delta \varphi L) / 8 $, co jest zgodne z rezultatem teorii I rzędu.
Wyrażenia $(\ref{48})$ i $(\ref{51})$ wskazują, że w granicznym stanie krytycznym (Eulera):
- Funkcja ugięcia m postać $(\ref{48})$, to znaczy jest mieszaną funkcji liniowej z sinusoidą oraz kosinusoidą, można więc ją sprowadzić do szeregu Fouriera
- Strzałka ugięcia pręta jest skończona i niezerowa.
W przykładzie 1 rozważono trzy przypadki szczególne belek-słupów.
Energia potencjalna w naprężeniach i przemieszczeniach
W artykule przedstawimy rozwiązanie klasycznej belki Bernoulli-Euler , obciążonej rozłożonym obciążeniem porzecznym $q(x)$ (rys. 3), ale także ściskanej siłą $N(x)$, oraz poddaną działaniu poprzecznych obciążeń skupionych $[ V, H, M](x_k)$ zlokalizowanym w punkcie o współrzędnej $x_k$. Indeks k oznacza w zależności od kontekstu miejsce przyłożenia podpory skupionej lub obciążeń skupionych o różnej naturze
Powszechnie znaną zależność na gęstość energii potencjalnej ciała sprężystego zapiszmy w postaci [8]:
$$ \begin{equation} \Phi = \cfrac{1}{2E} \{ \sigma_x^2 + \sigma_y^2+ \sigma_z^2 -2\nu (\sigma_x \sigma_y + \sigma_x \sigma_z +\sigma_y \sigma_z)+ 2(1+\nu)(\tau_{xy}^2+\tau_{xz}^2+ \tau_{yz}^2) \} \label{52} \end{equation} $$
gdzie w celu zwiększenia czytelności – zastosowano symbole klasyczne (x,y,z) ( a nie zwykle stosowany zapis wskaźnikowy).
W przypadku zginania poprzecznego wzór ten znacznie upraszcza się. Ponieważ wpływ naprężeń $\sigma_y$, $\tau_{xz}$,$\tau_{yz}$ jest pomijalnie mały, więc przyjmuje się, że różne od zera są jedynie naprężenia $\sigma_x$ , $\tau_{xy}$.
W podejściu klasycznym Bernoulli-Euler pomija się również wpływ naprężeń $\tau_{xy}$ na przemieszczenia i energię potencjalną. Wpływ naprężeń stycznych uwzględnimy w rozdziale dotyczącym belki Timoshenko.
Przy poczynionych założeniach upraszczających energia potencjalna belki zginanej wynosi:
$$ \begin{equation} \Pi_\sigma= \iiint \limits_V{\Phi \cdot dV}=\cfrac {1} {2E} \iiint \limits _V \sigma^2_{xy} \cdot dV\label{53} \end{equation} $$
Energię $(\ref{53})$ wyrazimy przez siły przekrojowe, a następnie przez przemieszczenia. Zachodzi:
$$ \begin{equation} \sigma_x= -\cfrac {N}{A} + \cfrac {M(x)}{I_y} \cdot z \label{54} \end{equation} $$
gdzie: A – pole przekroju pręta, Iy – moment bezwładności przekroju względem osi głównej y, z-odległość (pionowa) punktu przekroju od osi głównej. Dalej moment bezwładności oznaczamy bez indeksu.
W przypadku przeważającego zginania przy obliczaniu minimum energii układu wpływ rozciągającej siły osiowej N jest pomijalnie mały. Wpływ siły ściskającej N uwzględniamy w dalszej części opracowania.
Całkę w równaniu $(\ref{53})$ przekształcimy tak, by otrzymać energię w przemieszczeniach.
Ponieważ zachodzi związek pomiędzy momentem zginającym a krzywizną
$$ \begin{equation} M(x) = – EIw”(x) \label{55} \end{equation} $$
więc mamy:
$$ \begin{equation}\Pi_\sigma= \iiint \limits_V {\cfrac {1}{2E}}\left [ – \cfrac {EIw”(x)} {I}z \right ]^2 dV=\cfrac {E} {2} \iiint \limits_V [w”(x) \}^2 \cdot z^2 dV =
\cfrac {E} {2} \int \limits _0 \limits^L [w”(x)]^2 dx \cdot \iint \limits_A z^2 dA = \cfrac {EI_y} {2} \int \limits_0 \limits^L [w”(x) \}^2 dx \label{56} \end{equation} $$
Praca zewnętrznych sił poprzecznych
Praca zewnętrznych sił poprzecznych q(x) i momentów zginających m(x) rozłożonych na długości belki wynosi
$$ \begin{equation} L_q = \int \limits _0 \limits^L [q(x)w(x) +m(x) w'(x)] dx\label{57} \end{equation} $$
Obciążenia skupione zlokalizowane w punkcie $x_k$ możemy zapisać następująco:
$$ \begin{equation} Q_k=Q\cdot \delta (x_k) \quad ; \quad M_k=M\cdot \delta (x_k) \label{58} \end{equation} $$
gdzie delta (symbol) Diraca zdefiniowany jesr następująco
$$ \begin{equation} \delta (x_k)=\begin {cases}
1 ,& \text {dla} x = x_k \\
0 ,& \text {dla} x \ne x_k \\
\end {cases} \label{59} \end{equation} $$
Po podstawieniu $(\ref{58})$ do $(\ref{57})$ mamy:
$$ \begin{equation} L_q = \int \limits_0 \limits^L [q(x)w(x) +m(x) w'(x)] dx+\sum \limits _k Q_k \cdot w(x_k) + M_k \cdot w'(x_k) \label{60} \end{equation} $$
Efekt drugiego rzedu P-Δ. Praca siły ściskającej N
Praca siły ściskającej N(x) jest wykonywana na przemieszczeniu poziomym belki, wywołanym jej ugięciem, przy tym w zagadnieniu belki pomijamy skrócenie belki pod wpływem przekrojowej siły N. W celu oszacowania tego przemieszczenia należy przeanalizować zagadnienie II rzędu zilustrowane na rys.7.

Rys. 7 Przemieszczenie poziome Δ(dx) przekroju belki wywołane jej ugięciem w(x) – Efekt drugiego rzedu P-Δ przy ściskaniu
Analizujemy odcinek dx belki. Ugięcie powoduje przemieszczenie punktów końcowych do położenia dx’. Pozioma odległość Δ(dx) pomiędzy punktem przed i po ugięciu na której pracuje siła N wynosi:
$$ \begin{equation} \Delta(dx)= dx=dx’=dx- \sqrt{(dx)^2-(\cfrac {d}{dx} dx)^2} =dx \left [ 1-\sqrt{1-(w’)^2} \right ] \label{61} \end{equation} $$
Po rozwinięciu w szereg MacLaurina ostatniego pierwiastka w ($\ref{61}$) i pozostawieniu pierwszego nieliniowego wyrazu szeregu, otrzymamy aproksymację, która jest podstawą teorii drugiego rzędu całej mechaniki konstrukcji:
$$ \begin{equation} \Delta(dx) \approx dx \left [ 1 – \left ( 1 – \cfrac{1}{2} {w’}^2 + \cfrac{1}{10} {w’}^4 – … \right) \right ] \approx \cfrac{1}{2}\left \{ w'(x) \right \}^2 \label{62} \end{equation} $$
W rezultacie praca siły siły podłużnej, ściskającej, wynosi:
$$ \begin{equation} L_N = \int \limits_0^L N \cdot \Delta(dx) dx = \int \limits_0^L \cfrac{N}{2} \{w'(x) \} ^2 \label{63} \end{equation} $$
Praca sił odporu
Siły odporu qz(x) są reakcją podłoża r(x) na nacisk belki. Z definicji zależą one bezpośrednio od przemieszczenia w(x) zgodnie z zależnością ($\ref{64}$), gdzie znak (-) wynika ze zwrotu sił odporu, które są przeciwnie skierowane do ugięcia.
$$ \begin{equation} q_z (x) = – r(x) = – c(x) \cdot w(x) \label{64} \end{equation} $$
Na rys. 5 pokazano odpór pionowy v podłoża o stałej sprężystości cv (małe c) oraz odpór obrotowy m podłoża o stałej sprężystości cm. Analogicznie można zdefiniować odpór poziomy h podłoża o stałej sprężystości ch.
Reakcje V, H lub M skupionych w punkcie xk podpór sprężystych o stałych sprężystości odpowiednio CV, CH i CM (duże C) otrzymamy z definicji ($\ref{64}$) i zapisu w funkcją Diraca (\ref{59}$) Na przykład dla sprężystej podpory pionowej i obrotowej mamy (rys.9):
$$ \begin{equation} q_k= – {C_k}^V \cdot w(x) \cdot \delta(x_k) \quad ; \quad M_k=-{C_k}^M \cdot w^{’}(x) \cdot \delta(x_k) \label{65} \end{equation} $$
Nadając przyrosty przemieszczenia w nadajemy przyrosty siłom odporu. Zapisując pracę wirtualną sił odporu w postaci:
$$ \begin{equation} \delta L_c= -\int \limits _0 \limits ^L – r(x) \cdot \delta w(x) dx = -\int \limits _0^L c \cdot w(x) \delta w(x) dx\label{66} \end{equation} $$
widzimy, że nie ma możliwości wyłączenia znaku wariacji przed całkę. Będziemy to mogli zrobić dopiero wtedy, gdy ostatnią całkę zapiszemy w postaci:
$$ \begin{equation} \delta L_c= -\int \limits _0^L c \cdot \cfrac {1}{2} \cdot \delta [w^2(x)] dx = – \delta \int \limits_0 ^L \cfrac {c}{2} \cdot w^2(x)dx \label{67} \end{equation} $$
Teraz możemy opuścić znak wariancji (uprościć) i otrzymamy wyrażenie na pracę sił odporu podłoża sprężystego
$$ \begin{equation} \delta L_c= \int \limits_0 ^L \cfrac {c}{2} \cdot w^2(x) dx \label{68} \end{equation} $$
Stosowne wyrażenia na pracę skupionych podpór sprężystych, otrzymamy z ($\ref{68}$) po podstawieniu ($\ref{64}$).
Funkcjonał Lagrange’a belki Bernoulli
Ostatecznie funkcjonał Lagrange’a dla belki Bernoulli ściskanej siłą osiową N, na podłożu sprężystym i na podporach sprężystych zapiszemy w postaci:
$$ \begin{equation} \Pi_B = \int \limits_0^L \left [ \cfrac {EI_y} {2} [w”(x)]^2 – \cfrac {N(x)} {2}[w'(x)]^2 +\cfrac {c^v (x)}{2}[w(x)]^2 – q(x) w(x) \right ] dx +\sum \limits k \left[ \cfrac{{C_k}^v } {2} [w(x_k)]^2+\cfrac{{C_k}^M} {2} [w'(x_k)]^2 – V_k w(x_k) – M_kw’ (x_k) \right ] \label{69} \end{equation} $$
Poszukiwanie punktu stacjonarnego (minimalizację) funkcjonału ($\ref{69}$) możemy przeprowadzić dowolną metodą, na przykład metodą Ritza, w sposób pokazany niżej dla belki Timoshenko.
Belka Timoshenko
Belka Timoshenko została zdefiniowana ok 1956 roku w ramach teorii Timoshenko-Ehrenfest. Timoshenko odszedł od założenia prostopadłości przekroju odkształconego do osi belki na rzecz niewielkich odchyleń pochodzących od odkształceń postaciowych powierzchni przekroju poprzecznego, wywołanych naprężeniem stycznym od sił poprzecznych. Stan wytężenia i odkształcenia belki Timoshenko jest zwykle używany w analizie prętów o dużych wymiarach poprzecznych lub o przekrojach złożonych (np. kratownice, słupy złożone). Belka Timoshenko jest skuteczna również dla belek klasycznych, jest więc uogólnieniem belki Bernoulli.
Belka Timoshenko uwzględnia korektę pierwszego rzędu teorii zginania poprzecznego Bernoulli.
Klasyczne założenie kinematyczne zesztywnienia jest zmienione na „płaski przekrój pozostaje płaski, ale nie jest już prostopadły do zdeformowanej osi obojętnej belki”. Dla prostego elementu belkowego założenie to zilustrowano na rys. 10.
Zastosowanie belki Timoshenko w praktyce eliminuje kryterium podziału belek i tarcz ( p. definicja -inżynierska-teoria-zginania/#Belka”>belki).
Współczynnik β efektywności przekroju na ścinanie
Na rys. 9 przedstawiono ugięcia zastępczego elementu belki o długości dx obciążonego na czole obciążeniem $qz= \tau$ Element w rzeczywistości odkształca się mniej więcej tak, jak pokazano na rysunku 9c, czyli nieliniowo z paczeniem przekroju czołowego. W teorii Timoshenko-Ehrenfest przyjmuje się założenie upraszczające, że przekrój pozostaje płaski, ale obraca się o kąt $\gamma_v$ jak na rys. 9d.

Rys.9: Ilustracja do współczynnika efektywności przekroju [11]
W celu oszacowania średniego odkształcenie ścinającego$\gamma_v (rys. 9d)$, przyjętego w definicji belki Timoshenko – porównamy pracę wszystkich odkształcających się włókien (rys. 9a,b ) z pracą całego przekroju siłą poprzeczną V jak następuje:
$$ \begin{equation} \int\limits_A \cfrac{1}{2} \cdot dw \cdot \tau dA = \cfrac {1}{2} dw_v \cdot V \label{70} \end{equation} $$
Podstawiając do ($\ref{70}$) związek $dw = \gamma \cdot dx $, wynikający z rys. 9b oraz wyrażenie, wynikające z fizycznego prawa ścinania
$\tau= G \cdot \gamma$,
gdzie moduł Kirchoffa związany jest z dwoma innymi stałymi teorii sprężystości, związkiem
$$ \begin{equation} G=\cfrac{E}{2\cdot (1+\nu)} \label{71} \end{equation} $$
otrzymujemy:
$$ \begin{equation} \int \limits_A \cfrac{1}{2} \left(\cfrac{\tau}{G} \right) \cdot dx \cdot \tau \cdot dA = \cfrac{1}{2} \left ( \cfrac{\tau_v }{G} \right ) \cdot dx \cdot V \label{72} \end{equation} $$
Należy zachować ostrożność przy interpretacji $\tau_v$. NIE jest to po prostu średnie naprężenie ścinające, uzyskane przez równomierne rozłożenie siły poprzecznej (ścinającej) V na całym przekroju poprzecznym. Naprężnie $\tau_v$ rozkłada się tylko na części przekroju, nazywanym polem przekroju na ścinanie lub polem ścinania. $A_v$ .
$$ \begin{equation} A_v = \beta \cdot A \label{73} \end{equation} $$
gdzie $\beta$ jest współczynnikiem efektywności przekroju na ścinanie
$$ \begin{equation} \int \limits_A \cfrac{1}{2}\cdot \cfrac{1}{G}\cdot dx \cdot \left ( \cfrac{V \cdot \cdot \overline S_y}{I_y \cdot t}\right) ^2 dA= \cfrac{1}{2} \cdot \cfrac{V}{\beta A}\cdot \cfrac{1}{G} \cdot dx \cdot V \label{74} \end{equation} $$
5/6 ,& \text {dla przekroju prostokątnego }\\
9/10 ,& \text {dla przekroju kołowego }\\
1/2 ,& \text {dla rury okrągłej }\\
∼ \cfrac{A_{||}}{A} & \text {dla przekrojów cienkościennych }\\
\end {cases} \label{76} \end{equation} $$
W przypadku przekrojów cienkościennych, złożonych z cienkich ścianek prostokątnych $A_{||}$ jest polem przekroju ścianek równoległych do siły poprzecznej (ścinającej). Na przykład dla dwuteownika ścinanego siła poprzeczną pionową jest to pole przekroju środnika ; dla dwuteownika ścinanego siła poziomą jest to suma pół przekroju pasów; dla rusy prostokątnej jest to suma pól ścianek równoległych do siły poprzecznej.
Obrót przekroju belki na skutek zginania i ścinania
W celu oszacowania obrotu przekroju belki poprzecznie zginanej zastosujemy równanie prac wirtualnych
$$ \begin{equation} \delta W_{ext} = \delta W_{int} \label{77} \end{equation} $$
gdzie praca sił zewnętrznych $F_{ext}$ na przemieszczeniach $\Delta_{ext}$:
$$ \begin{equation} \delta W_{ext} = \sum \limits_{i=1 }^{N_f} F_{ext,i } \cdot \Delta_{ext,i} \label{78} \end{equation} $$
gdzie:
$N_f$ – liczba wszystkich sił zewnętrznych
Praca sił wewnętrznych wynosi:
$$ \begin{equation} \delta W_{int } = \sum \limits_{e=1}^{N_e} \delta W_{N,e} + \delta W_{M,e} + \delta W_{V,e} \label{79} \end{equation} $$
gdzie: e – numer elementu (pręta itd) ; $N_e$ – liczba wszystkich elementów
Praca sił wewnętrznych w elemencie wywołana siła osiową (z pominięciem efektów 2-rzędu) $\delta W_N$ wynosi
$$ \begin{equation} \delta W_N = \delta N \cdot \cfrac{N \cdot L}{EA}\label{80} \end{equation} $$
Praca sił wewnętrznych w elemencie wywołana zginaniem $\delta W_M$ wynosi
$$ \begin{equation} \delta W_M = \int \limits_0^L \delta M \cdot \cfrac{M \cdot }{EI} dx \label{81} \end{equation} $$
Praca sił wewnętrznych w elemencie wywołana ścinaniem $\delta W_V$ wynosi
$$ \begin{equation} \delta W_V = \int \limits_0^L \delta V \cdot \cfrac{V \cdot }{G A_v } dx \label{82} \end{equation} $$
gdzie aproksymację rzeczywistego przemieszczenia ścinania przyjęto po analizie przedstawionej w rozdziale Współczynnik efektywności przekroju na ścinanie.
Wnioskiem z wyrażenia ($\ref{77}$) do ($\ref{82}$) jest podstawowa hipoteza Timoshenko, że:
$$ \begin{equation} \varphi = \Theta +\gamma_v \label{83} \end{equation} $$
gdzie:
$\varphi$ – całkowity kąt obrotu w węźle belki ( obrotowy stopień swobody węzła belki Timoshenko),
$\Theta = \cfrac{dw}{dx}$, – kąt obrotu wg teorii Bernoulli- Euler ($\ref{10}$) ,
$\gamma_v =\cfrac{\tau_v}{G} = \cfrac{V}{G\cdot A_v}$ – kąt obrotu wywołany ścinaniem.
(Niżej pokażemy, że wskutek przyjętej w praktyce umowy znakowania sił poprzecznych, kąt $\gamma_v$ jest ujemny).
Ponieważ $\varphi \neq \Theta$, więc przekroje poprzeczne pozostają na swoim miejscu, ale obrót przekroju poprzecznego nie jest już równy obrotowi osi belki ($\ref{10}$), czyli nie jest prostopadły do niej.
Powyższe rozważania zilustrowano na rys. 10 .

Rys.10. Dwuwęzłowy element Timoshenko: a) przemieszczenia końców elementu, w tym $\varphi=\Theta + \gamma $ – kąt ugięcia b) kąt ugięcia $\Theta$ od zginania, c) kąt ugięcia $\gamma$ od ścinania (zmodyfikowane i uzupełnione [12]
Pole przemieszczeń i smukłość ścinania
Pole przemieszczeń liniowych węzła belki Timoshenko możemy zapisać w postaci:
$$ \begin{equation} u(x,z)=u(x)+z \cdot \varphi \quad \wedge \quad w(x,z)=w(x) \label{84} \end{equation} $$
Zachodzi [12]:
$$ \begin{equation} \varphi =\Theta \, – \gamma_v \quad \wedge \quad \gamma_v = \cfrac {\tau} {G}=\cfrac {V}{GA_v} = \cfrac{V}{S_v} \label{85} \end{equation} $$
gdzie:
$V$ siła poprzeczna
$\gamma_x = \gamma_v (x)$ jest „obrotem ścinania”, uśrednionym po przekroju,
$G$ – moduł Kirchoffa,
$A_v$ efektywne polem przekroju na ścinanie,
$S_v = GA_v$ sztywność postaciowa(na ścinanie) przekroju belki.
Ujemny znak $\gamma_v$ w ($\ref{85}$} wprowadzono ze względu na umowę znakowania sił poprzecznych w technice, zgodnie z którą dodatnie siły tnące wywołują ujemny obrót ścinania (rys. 11).
W przypadku przekrojów złożonych (wielogałęziowych) sztywność postaciową $S_v$ wyznacza się na drodze porównania przemieszczeń konstrukcji rzeczywistej oraz pręta zastępczego. Dla najczęstszych przypadków prętów o pasach równoległych wyrażenia na sztywność postaciową podano w tab.3.
Tab.3. Podatność na ścinanie prętów złożonych [13]

Felippa (2001) [12] macierz sztywności elementu skończonego belki Timoshenko uzyskał poprzez modyfikację macierzy sztywności standardowego elementu belkowego Bernoulli poprzez wprowadzenie bezwymiarowego współczynnika $\Psi$, charakteryzującego „smukłość ścinania” . Wyrażenie na współczynnik $\Psi$ jest wnioskiem z powyższych analiz:
$$ \begin{equation} \Psi = \cfrac {12\cdot EI_y} {GA_v L^2}= \cfrac {12 \cdot EI_y} {S_v\cdot L^2} =12 \cdot k \label{86} \end{equation} $$
gdzie
$$ \begin{equation} k = \cfrac{EI_y}{S_v \cdot L^2} \label{87} \end{equation} $$
jest współczynnikiem podatności na ścinanie wprowadzonym w pracy Pałkowski (2009) [13].
Jeśli $\Psi \to 0$, to model Timoshenko redukuje się do klasycznego modelu Bernoulli-Euler .
Energia potencjalna belki Timoshenko
Energia potencjalna belki Timoshenko wynosi:
$$ \begin{equation} \Pi_\sigma= \cfrac {1}{2} \left \{ \int_0^L EI \cdot \kappa^2(x) +S_v \cdot \gamma^2(x) \right \} dx \label{88} \end{equation} $$
gdzie:
krzywizna $ \kappa =w”(x)$ .
Kąt ścinania wyznaczymy z wykorzystaniem zależności różniczkowej pomiędzy siłami przekrojowymi, jak następuje:
$$ \begin{equation} V(x)=\cfrac {dM(x)} {dx} \quad \wedge \quad M(x) = – EI_y w”(x) \to V(x)= – EI_y w^{3′}(x) \to \gamma(x)= \cfrac {V(x)}{GA_v} =-\cfrac {EI_y} {GA_v} w^{3′}\label{89} \end{equation} $$
Stąd otrzymujemy:
$$ \begin{equation} \Phi_{\sigma} =\cfrac {EI_y} {2} \int \limits_0 \limits ^L \left ( [w”(x)]^2 + \cfrac {EI_y} {GA_z} [w^{3′}(x)]^2 \right )dx \label{90} \end{equation} $$
Praca siły ściskającej i praca sił zewnętrznych
Ponieważ przy wyprowadzaniu formuł na pracę siły ściskającej oraz pracę sił zewnętrznych, w tym odporu podpor sprężystych nie wykorzystywaliśmy hipotezy kinematycznej Bernoulli, więc zależności słuszne są również dla belki Timoshenko.
Funkcjonał Lagrange’a belki Timoshenko
Ostatecznie funkcjonał Lagrange’a dla belki Timoshenko ściskanej siła osiową N , na podłożu sprężystym i na podporach sprężystych zapiszemy w postaci ($\ref{91}$):
$$ \begin{equation} \Pi_T = \int \limits_0\limits^L \left [ \cfrac {EI_y} {2} \left ( [w”(x)]^2 + \underline { \cfrac {EI_y} {GA_v} [w^{3′}(x)]^2 } \right ) – \cfrac {N(x)} {2}[w'(x)]^2 +\cfrac {c^v (x)}{2}[w(x)]^2-q(x)w(x) \right ] dx \\ +\sum \limits k \left[ \cfrac{{C_k}^v } {2} [w(x_k)]^2+\cfrac{{C_k}^M} {2} [w'(x_k)]^2- V_k w(x_k)- M_kw’ (x_k) \right ] \label{91} \end{equation} $$
Funkcjonał $\Pi_T$ dla belki Timoshenko różni się od funkcjonału belki Bernoulli ($\ref{69}$) jedynie składnikiem zaznaczonym podkreśleniem w formule ($\ref{91}$).
Minimalizacja funkcjonału Lagranga metodą Ritza
Układ kanoniczny równań Lagranga-Ritza dla belek
W metodzie Ritza poszukuje się funkcji realizującej minimum funkcjonału w klasie wielomianów:
$$ \begin{equation} w(x)=\varphi_0+\sum \limits_{i=1} \limits ^{n} a_i \cdot \varphi_i (x) \label{92} \end{equation} $$
gdzie funkcje bazowe φi(x) mogą być dowolnie dobrane, ale tak, by były dopuszczalne, czyli spełniały kinematyczne warunki brzegowe.
Uwaga: oznaczenie funkcji bazowych $\varphi_i (x)$ jest podobne do oznaczenia stopnia swobody obrotu. węzła $\varphi$. Należy zachować ostrożność i interpretować symbol zależnie od kontekstu.
Podstawiając ($\ref{92}$) do ($\ref{91}$), otrzymamy wyrażenie na $\Pi_T$ , w którym funkcjonał energii jest teraz funkcją współczynników ai. Funkcja ta osiąga minimum przy warunku:
$$ \begin{equation} \cfrac {\partial \Pi_T} {\partial a_j}=0 \quad , \quad (j=1,2,…n) \label{93} \end{equation} $$
Różniczkując ($\ref{91}$) (z aproksymacją ($\ref{92}$) podług $a_i$ i przyrównując pochodną do zera zgodnie z ($\ref{93}$) po niewielkich przekształceniach otrzymamy układ równań liniowych
$$ \begin{equation} \sum \limits_{i=1} \limits^n a_i \delta_{ij}= \Delta_{Fj} \quad , \quad (j=1,2,…n) \label{94} \end{equation} $$
gdzie współczynniki $\delta_{ij}$ oblicza się z wyrażeń:
$$ \begin{equation}\delta_{ij} \stackrel{\rm def}{=} \sum \limits_{i=1} \limits^n \int \limits_0 \limits ^L \left [ EI_y \left (\varphi^{„}_{i}(x) \varphi^{„}_{j}(x) + \cfrac {EI_y} {GA_z} \varphi^{3′}_{i}(x) \varphi^{3′}_{j}(x)\right) – \\ – N \varphi^{’}_{i}(x) \varphi^{’}_{j}(x) + c^V \varphi_{i}(x) \varphi_{j}(x)\right] dx + \sum \limits_k [C^V_k \varphi_{i}(x_k) \varphi_{j}(x_k)+C^M_k \varphi^{’}_{i}(x_k) \varphi^{’}_{j}(x_k)] \label{95} \end{equation} $$
a wyrazy wolne z wyrażeń:
$$ \begin{equation} \Delta_{Fj} \stackrel{\rm def}{=} \int \limits_0^L q(x) \varphi_{j}(x)dx+\sum \limits_k V_k \varphi_{j}(x_k)+\sum \limits_k M_k \cdot \varphi^{’}_{j}(x_k) \label{96} \end{equation} $$
Współczynniki metody $a_i$ są wartościami przemieszczeń węzłowych ęzłów belki $ (przemieszczeń liniowych $w_i$, $u_i$ oraz kątowych $\varphi_i$ zestawionych dla wszytskich węzłów belki)
Przy tych oznaczeniach otrzymujemy następujący układ równań liniowych , z których wyznaczamy współczynniki $a_i$ – przemieszczenie węzłowe belki:
$$ \begin{equation} \sum \limits_{i=1}^n a_i \cdot \delta_{ij} = \Delta_{Fj} \quad , \quad (j=1,,2, …, n) \label{97} \end{equation} $$
który nazywać będziemy kanonicznymi równaniami Lagranga-Ritza dla belek
Układ kanoniczny ($\ref{97}$)ma analogiczną strukturę do klasycznych wyrażeń mechaniki konstrukcji w metodzie sił lub metodzie przemieszczeń. W kolejnym rozdziale układ kanoniczny równań zapiszemy w postaci macierzowej i pokażemy, że zachowane są ogólne zależności metody elementów skończonych, w której zastosowano funkcje kształtu w postaci wielomianów Ritza.
Macierzowe sformułowanie zagadnienia belki Timoshenko
Równanie kanoniczne
Sformułowanie macierzowe zagadnienia belki Timoshenko polega na zastosowaniu standardowego algorytmu metody MES z korektami macierzy sztywności i innych zmiennych, wynikającą z teorii belki Timoshenko przedstawionej w poprzednich rozdziałach. Zagadnienie belki Bernoulli jest zagadnieniem szczególnym (zdegenerowanym) zagadnienia belki Timoshenko.
Równanie kanoniczne ($\ref{94}$) w zapisie macierzowym przyjmuje postać
$$ \begin{equation} [K] \cdot |a| =|F| \label{98} \end{equation} $$
gdzie:
$[K]$ – macierz sztywności konstrukcji o wymiarze (mxm),
$|a|$ – poszukiwane przemieszczenia węzłowe o wymiarze m
$|F|$ – wektor obciążenia o wymiarze m,
gdzie:
$m = w \cdot sq $ – wymiar macierzy i wektorów
w- całkowita liczba węzłów konstrukcji,
q – liczba stopni swobody w węźle. W przypadku rozpatrywanej belki q=2.
Ponadto w tym rozdziale stosuje się następujące oznaczenia:
$N^{[e]}$ – liczba elementów w konstrukcji,
[e] (1, .., N^{[e]}$- numer elementu,
(i) ( 1, …, w) – numer węzła w konstrukcji,
$ r,s (1, …, m)$ – numery wierszy i kolumn globalnej macierzy sztywności.
Poniżej opisano algorytm postępowania przy rozwiązaniu belki metodą Lagrange’a – Ritza, uwypuklając elementy merytoryczne nad automatyzację obliczeń. Algorytm jest analogiczny do metody elementów skończonych MES [14], [15]. Technika metody MES jest bardziej rozbudowana (np. w algorytmie siły osiowe w belce-słupie są elementem rozwiązania, więc liczba stopni swobody w węźle wynosi q=3 , a w niniejszej prezentacji przyjęto siły osiowe jako znane.
Dyskretyzacja belki
Obliczenia prowadzi się na konstrukcji podzielonej na elementy skończone, np. w sposób pokazany na rys,12, gdzie belka podzielona jest na 8 elementów skończonych zawiera w = 9 węzłów, to znaczy
Macierz sztywności belki i elementu
Macierz sztywności belki uzyskujemy drogą superpozycji – „sumowania” po elementach skończonych, według zasady:
$$ \begin{equation} [K]=[K_{rs}]=[\sum [k_{ij}]^{(e)} ] \label{99} \end{equation}$$
gdzie sumowanie dotyczy elementów $[e]$ zbiegających się w węźle (i)
Ścisła macierz sztywności $k^{[e]}$ elementu [e] belki z uwzględnieniem ściskania przy założeniu rozdzielenia stanu zgięciowego od ściskania i bez uwzględnienia sztywności podłużnej, ale z uwzględnieniem sztywności postaciowej jest złożona z podmacierzy $[k_{ij}]^{[e]}$ mających postać:
$$\begin{equation} [k_{ij}]^{[e]} = \cfrac{EI}{1+ 12 \Psi } \begin{bmatrix}
\cfrac{12}{L^3} \Lambda_1 & \cfrac{6}{L^2} \Lambda_2 & \cfrac{- 12}{L^3} \Lambda_1 & \cfrac{6}{L^2} \Lambda_2\\
\cfrac{6}{L^2} \Lambda_2 & \cfrac{4 \cdot (1+3\Psi )}{L} \Lambda_3 & \cfrac{ – 6}{L^2} \Lambda_2 & \cfrac{2 (1-6 \Psi )}{L} \Lambda_4\\
\cfrac{- 12}{L^3} \Lambda_1 & \cfrac{- 6}{L^2} \Lambda_2 & \cfrac{ 12}{L^3} \Lambda_1 & \cfrac{6}{L^2} \Lambda_2\\
\cfrac{6}{L^2} \Lambda_2 & \cfrac{2 (1- 6 \Psi )}{L}\Lambda_4 & \cfrac{-6} {L^2} \Lambda_2 & \cfrac{4 (1+3 \Psi}{L} \Lambda_3\\
\label {100}\end{bmatrix} \end{equation}$$
gdzie $L=L_{[e]}$, $EI = EI_{y, [e]}$, $\Psi=\Psi{[e]}$ – długość L , sztywność giętna $EI_y$ oraz funkcja $\Psi$ (${86}$) dla elementu [e].
Macierz ($\ref{100}$) powstała w wyniku uogólnienia konwencjonalnej macierzy sztywności elementu belkowego o:
- sztywność postaciową (element Timoshenko) współczynnikami $\Psi$ ($\ref{86}$),
- działanie siły osiowej ściskającej współczynnikami $Λ_i$, współczynnikami $\Lambda$ ($\ref{45$).
Sposób składania (inaczej „sumowania” lub „agregacji”) macierzy systemu $[K]$ z macierzy elementów $ [k_{ij}]^{[e]}$ pokazano w przykładzie rachunkowym 2.
Wektor równoważników węzłowych
Dla obciążenia równomiernie rozłożonego pomiędzy węzłami obciążenie węzłowe jest równe reakcjom belki utwierdzono-utwierdzonej. W przypadku podłoża sprężystego kinematyczne warunki brzegowe zastąpiono statycznymi warunkami brzegowymi. Można wykazać [12], że równoważniki węzłowe dla elementu belki Timoshenko [e] o długości $L$ , spoczywającego na podłożu sprężystym o stałej sprężystości $c_i$ wynoszą:
$$ \begin{equation} | F_i | = ( – c_i) \cdot L^2 \begin{bmatrix}
\cfrac {13 + 294 \Psi +1680 \Psi^2}{L} & \cfrac {11 + 231 \Psi + 1260 \Psi^2}{6} & \cfrac{3 \cdot ( 3 + 84 \Psi + 560 \Psi^2) } {2 \cdot L } & \cfrac{ 13 + 378 \Psi + 2520 \Psi^2} {12}\\
\cfrac {11 + 231 \Psi + 1260 \Psi^2}{6} & \cfrac{ 1 + 21 \Psi + 126 \Psi^2} {3} \cdot L & \cfrac { 13 + 378 \Psi + 2520 \Psi^2} {12} & \cfrac{ 1 + 28 \Psi + 168 \Psi^2} {4} \cdot L \\
\cfrac {3 \cdot ( 3 + 84 \Psi + 560 \Psi^2) } {2\cdot L} & \cfrac { 13 + 378 \Psi + 2520 \Psi^2} {12} & \cfrac{ 13 + 294 \Psi + 1680\Psi^2} {L} & \cfrac{11 + 231 \Psi + 11260 \Psi^2}{6} &\\
\cfrac {13 \cdot 378 \Psi + 2520 \Psi^2 } {12} & \cfrac { 1 + 28 \Psi + 168 \Psi^2} {4}\cdot L & \cfrac{ 11 + 231 \Psi + 1280\Psi^2} {6} & \cfrac{1 + 21 \Psi + 126 \Psi^2}{3} \cdot L &\\
\end{bmatrix}
\cdot \begin{vmatrix} w_1 \\ \varphi_1 \\ w_2 \\ \varphi_2 \end{vmatrix}
\label {101} \end{equation}$$
Kinematyczne warunki brzegowe. Modyfikacja macierzy sztywności
Jeśli przemieszczenia punktu węzłowego jest znane , np. $w_k = \delta$, to w celu zachowania wymiaru macierzy sztywności – kinematyczne warunki brzegowe można uwzględnić w następujący sposób:
1) przenieść iloczyny odpowiednich wyrazów macierzy sztywności i tego przemieszczenia na stronę wektora równoważników węzłowych:
$$ \begin{equation} \overline F_r = -k_{rk}\cdot \delta + F_r \label{102} \end{equation}$$
gdzie:
$F_r$ – r-ty wyraz pierwotnego wektora równoważników węzłowych,
$\overline F_r $ zmodyfikowany wyraz wektora równoważników węzłowych
a jednocześnie:
2) w wierszu macierzy , odpowiadającemu przemieszczeniu należy wyzerować wszystkie wyrazy oprócz wyrazu na głównej przekątnej, który przyjmujmy równy 1.
Jednym zdaniem do równania kanonicznego wpisujemy równość $ w_k = \delta $.
Wyrazy zawierające niewiadome przemieszczenia, występujące w wektorze obciążenia przenosi się na stronę macierzy sztywności.
Zmodyfikowane równanie kanoniczne
Po operacjach kondensacji macierzy sztywności, otrzymuje się zmodyfikowane równanie równanie kanoniczne
$$ \begin{equation} [ \widetilde K ] \cdot |u|=|\overline F| \label{103} \end{equation}$$
Przemieszczenia węzłowe
Przemieszczenia węzłowe uzyskuje się w drodze rozwiązania liniowego układu równań )$\ref{103}$).
Siły przekrojowe
Siły przekrojowe obliczono dla każdego elementu (e), wykorzystując zasadę superpozycji tzn. wyznaczając siły przywęzłowe z sumy sił od znanych już przemieszczeń |u|^{(e)} węzłowych (krok odwrotny w celu określenia sił od przemieszczeń) i sił wywołanych obciążeniami elementowych |F^0|^{(e)}:
$$ \begin{equation} |F|^{[e]} = [K]^{[e]} \cdot |u|^{(e)} +|F^0|^{[e]}\label {104} \end{equation}$$
czyli dla indeksów w numeracji lokalnej:
$$ \begin{equation}
\begin{vmatrix} V_1 \\ M_1 \\ \hline V_2 \\ M_2\\ \end{vmatrix}^{(e)} =
\begin{bmatrix} k_{11} & | & k_{12} \\ \hline k_{21} & | & k_{22} \\ \end{bmatrix}^{(e)}
\cdot \begin{vmatrix} w_1 \\ \varphi_1 \\ \hline w_2 \\ \varphi _2\\ \end{vmatrix}^{(e)}
+ \begin{vmatrix} V_1^0 \\ M_1^0 \\ \hline V_2^0 \\ M_2^0 \\ \end{vmatrix}^{(e)}
\label{105} \end{equation}$$
W rezultacie uzyskamy wektory sił przywęzłowych | V [kN], M [kNm] | dla poszczególnych elementów skończonych.
Inne elementy belkowe
Element Własowa
Najbardziej ogólnym, klasycznym elementem belkowym jest element Własowa zdefiniowany w roku 1959. Własow wprowadził dodatkowy stopień swobody — paczenie i siłę przekrojową bimoment. Teoria prętów cienkościennych pozwala wyjaśnić zjawisko zwichrzenia oraz utraty stateczności giętno-skrętnej. W klasycznym sformułowaniu Własowa nie uwzględniano odkształceń postaciowych.
Element Timoshenko – Własow
Uogólnionym elementem belkowym (prętowym), krótko elementem prętowym, jest element cienkościenny Własowa z uwzględnieniem efektu Timoshenko. Element Timoshenko oraz element Bernoulli są przypadkami szczególnymi elementu prętowego.
Element Timoshenko-Własow opisany opisany zgodnie z teorią Rankine–Timoshenko–Własowa umożliwia uwzględnienie:
- .ginania z wuzględniem ścinania (efektu Yimoshenki),
- skręcania swobodnego i nieswobodnego z paczeniem (efekt Własowa),
- Wyboczenie giętno-skrętne (zwichrzenie) cienkościennych profili otwartych.
- „Coupling efekt”, czyli sprzężenie zginania, skręcania i ścinania,
Specjalne elementy belkowe
Specjalne elementy belkowe, w tym belka Timoshenko i belka Własowa, są stosowane w różnych zastosowaniach inżynieryjnych, aby zwiększyć analizę konstrukcji.
Stosowane są również:
- element Reddy-Bickford przedstawione w pracy [16] – element Timoshenko wyższych rzędów
- belkowy element hybrydowy , stosowany w sytuacji bardzo smukłych belek, gdzie naprężenia osiowe są znacznie większe od naprężeń zginania,
- elementy mieszane, stosowane sytuacjach, gdy siła osiowa lub siły poprzeczne są traktowane jako niezależny stopień swobody.
Przykłady rachunkowe
Przykład 1 [ Belka Bernoulli nieliniowa geometrycznie ]
Korzystając z wyników podanych w rozdziału Ugięcie belki ściskanej w stanie utraty stateczności Euler przeprowadzić analizę nieliniowych geometrycznie: belki wolnopodpartej, słupa wspornikowego oraz belki z imperfekcjami osi w kształcie sinusoidy. Belkę pokazano an rys. 6b, słup wspornikowy na rys,. 6d).
Nieliniowa geometrycznie belka wolnodparta
Warunki brzegowe
Warunki brzegowe dla belki wolnopodpartej , pokazanego na rys. 6b są następujące:
$ u_1=0, \quad w_1=0, \quad V _2= – C_{\Delta} \cdot u_2 $
gdzie: $C_\Delta$ jest stała sprężystości podparcia w węźle prawym (2).
Stopnie swobody węzłów belki ; $\begin{vmatrix} \varphi_1 \\u_2 \\ w_2 \\ \varphi_2 \end {vmatrix}$
Zmodyfikowana macierz sztywności
$ [ \widetilde K ] = \cfrac{EI}{L} \cdot \begin{bmatrix}
4 \Lambda_3 & 0 & \cfrac {-6 \Lambda_2} {L} & 2 \Lambda_4 \\
& \cfrac {EA}{EI} & 0 & 0 \\
&& \cfrac{12 \Lambda_1 + c_{\Delta}} {L^2} & \cfrac { – 6\Lambda_2}{L} \\
& SYM & \cfrac{ -6 \Lambda_2}{L} & 4 \Lambda_3 \\
\end{bmatrix}$
Odwrotność macierzy zmodyfikowanej
$[\widetilde K]^{-1} = \cfrac{EI}{L} \cdot \begin{bmatrix}
k_{11} & 0 & 3 L \cdot \Lambda_2 & k_{14} \\
0 & \cfrac {\Theta} {EA} & 0 & 0 \\
& & L ( 2 \Lambda_3+\Lambda_4 ) & 3 L \cdot \Lambda_2 \\
& SYM & & k_{11} \\
\end{bmatrix}$
gdzie:
$c_\Delta= C_\Delta \cdot \cfrac{L^3}{EI}$,
$k_{11}= \cfrac{-9\Lambda_2^2 +(c_\Delta + 12 \Lambda_1 ) \Lambda_3}{2 \Lambda_3 -\Lambda_4}$,
$k_{14}=\cfrac{18 \Lambda_2^2 – (c_\Delta + 12 \Lambda_1 ) \Lambda_4}{4 \Lambda_3 -2 \Lambda_4}$,
$\Theta=EI\cdot [36 \Lambda_2^2 -(12\Lambda_1 +c_\Delta) \cdot (2 \Lambda_3 +\Lambda_4)]$.
Siły węzłowe od obciążenia rozłożonego $q$
$ [F]_q = \cfrac{qL}{2} [ L/6 \quad 0 \quad 1 \quad – L/6 ]$
Z formuły ($\ref {50} $) można wyznaczyć ścisłą strzałkę ugięcia pręta (z dokładnością do odstępstwa o d symetrycznych warunków brzegowych).
$ f_N =w_N (\xi=0,5) = \cfrac{qL^4}{48EI} \cdot \left ( \cfrac{12 (2 \Lambda_3 + \Lambda_4)}{36 \Lambda_2^2 +(12 c_\Delta \Lambda_1) \cdot (2 \Lambda_3 +\Lambda_4)} +\cfrac{\alpha \cdot tg(\alpha/2)}{2 \Lambda_3 +\Lambda_4 } \right)$
Strzałka ugięcia $f_q$ pochodząca od sił poprzecznych jest całką szczególną zagadnienia, która powinna być dodana do powyższej. W przypadku obciążenia równomiernie rozłożonego całka szczególna wynosi:
Ponieważ
$w_q (\xi) = \cfrac {qL^4}{24 EI}$, to
$f_q = w(\xi=0,5) = \cfrac{q L^4}{384 EI}$
Nieliniowy geometrycznie słup wspornikowy
Ściskany pręt wspornikowy obciążony jest poprzecznie siłą H w wierzchołku (1) lub obciążeniem równomiernie rozłożonym h po wysokości. Pręt jest sprężyście utwierdzony w węźle dolnym (1) ze stałą sprężystości $c_\varphi$. (rys.6c)
Warunki brzegowe
$ u_2=0, \quad w_2=0, \quad M_2 = – C^{\varphi} \cdot \varphi_2$
Macierz zmodyfikowana
W przypadku konstrukcji jednoelementowej macierz ($ \ref{42}$) jest macierzą konstrukcji i można ją bezpośrednio modyfikować w celu uwzględnienia warunków brzegowych . Modyfikację macierzy sztywności dokonamy poprzez wykreślenie wierszy i kolumn, odpowiadających odebranym stopniom swobody i bezpośrednie wpisanie warunku sprężystej podpory do równania kanonicznego. W rezultacie otrzymamy macierz zmodyfikowaną $[\widetilde K]$:
$ [ \widetilde{K]} = \cfrac{EI}{L} \cdot
\begin{bmatrix}
\cfrac{EA}{EI} & 0 & 0 & 0 \\
& \cfrac{12 }{L^2} \Lambda_1 & \cfrac{6}{L} \Lambda_2 & \cfrac{6}{L} \Lambda_2 \\
& & 4 \Lambda_3 & 2 \Lambda_4 \\
& SYM & & 4 \Lambda_3+ c_{\varphi} \\
\end{bmatrix}
\begin{matrix}
u_2 \\ w_2 \\ \varphi_2 \\ \varphi_1
\end {matrix} $
Odwrócona macierz zmodyfikowana
Poprzez odwrócenie macierzy $[ \widetilde K]$ uzyskujemy macierz podatności nieliniowego geometrycznie wspornika:
$ [ \widetilde{K]}^{-1} = \cfrac{L}{\Theta}
\begin{bmatrix}
\cfrac{\Theta}{EA} & 0 & 0 & 0 \\
& k_{22} & k_{24}+\Lambda_2 \cdot c_{\varphi}\cdot L/2 & k_{24} \\
& & 3 \Lambda_2^2 -\Lambda_1(4 \Lambda_3+c_{\varphi}) & -3\Lambda_2^2 +2 \Lambda_1 \Lambda_4 \\
& SYM & & 3 \Lambda_2^2 -4 ] \Lambda_1 \Lambda_3 \\
\end{bmatrix}
\begin{matrix}
u_2 \\ w_2 \\ \varphi_2 \\ \varphi_1
\end {matrix} $
gdzie:
$c_{\varphi}=C_{\varphi} \cdot L/EI$,
$k_{24}=\Lambda_2L(2\Lambda_3-\Lambda_4)$,
$k_{22}=-\cfrac{L^3}{3}\cdot(4\Lambda_3^2-\Lambda_4^2+\Lambda_3c_{\varphi} ) $,
$\Theta = – 4EI \cdot [(2\Lambda_3 -\Lambda_4) \cdot \Lambda_1 (2 \Lambda_3 + \Lambda_4) – 3\Lambda_2^2] + c_{\varphi} \cdot L \cdot ( 3 \Lambda_2^2 – 4 \Lambda_1 \Lambda_3) $.
Przemieszczenia i siły przekrojowe
Przemieszczenia węzłowe uzyskuje się poprzez przemnożenie macierzy podatności $ [ \widetilde{K]}^{-1}$ przez wektor równoważników węzłowych, który wynosi:
od poziomego obciążenia rozłożonego h
$ [F]_h = \cfrac{hL}{2} \cdot [\quad 0 \quad 1 \quad L/6 \quad -L/6 \quad ] $
od obciążenia skupionego H
$ [F]_H = [ \quad 0 \quad H \quad 0 \quad 0 \quad ] $
Moment zginający w utwierdzeniu wspornika otrzymuje się po dodaniu do odpowiedniej współrzędnej wektora siły przywęzłowej, uzyskanej z przemnożenia wektora przemieszczeń węzłowych przez ostatni wiersz pierwotnej macierzy sztywności.
Belka krzywoliniowa w kształcie sinusoidy
Równanie osi elementu
Oś elementu jest opisana ogólną zależnością parametryczną
$$\begin{equation} x=x(s) \quad y=y(s) \label {110} \end {equation} $$
Macierz podatności elementu
Macierz podatności (odwrotność macierzy sztywności) elementu opisanego równaniem paramteryczym )$\ref{110}$) , otrzymana po przprowadzeniu algorytmu opisanego w tekście wynosi
$$\begin{equation} [ K]^{-1} =\begin{bmatrix}
\int \limits_0^L \left ( \cfrac{(dx/dy)^2}{EA}+\cfrac{y^2}{EI} \right ) ds & \int \limits_0^L \left ( \cfrac{(dx/dy) (dy/ds)} {EA}- \cfrac{xy}{EI} \right ) ds & \int \limits_0^L \cfrac{y}{EI} ds \\
& \int \limits_0^L \left ( \cfrac{(dy/ds)^2} {EA}- \cfrac{x^2}{EI} \right ) ds & – \int \limits_0^L \cfrac{x}{EI} ds \\
& SYM & \int \limits_0^L \cfrac{1}{EI} ds \\
\end{bmatrix} \label {111} \end {equation} $$
gdzie: $A(s)$ i $I(s)$ – pole przekroju i moment bezwładności sa znanymi funkcjami współrzędnej bieżącej s
Wyrażenie ($\ref{111}$) jest zgodne z podanym w pracy [17].
Macierz podtności w standardowym układzie współrzędnych
Kształt elementu wstępnie wygiętego w łuk sinusoidy opisany w standardowym układzie współrzędnych (x,y)
$$\begin{equation} y= \eta \sin{(\omega x)} \label {112} \end {equation} $$
gdzie: $\eta$ – amplituda, $\omega= \cfrac{ n \ pi x}{L}$ – częstotliwość dla $n=1,2,3,…$ – liczba naturalna oznaczająca numer postaci własnej elementu,
Równanie krzywej ($\ref{112}$) można zapisać w postaci parametrycznej:
$$\begin{equation} y= \eta \sin {\alpha} ; \quad x =\cfrac{\alpha}{\omega} \label {113} \end {equation}$$
gdzie parametr $\alpha \in [ \quad \alpha_1= -\omega L/2 ; \quad \alpha_2= \omega L/2 \quad ]$,
Poszczególne wyrazy macierzy podatności $( \ref{111}$) wynoszą:
$k_{11}= \int \limits_0^L \left ( \cfrac{(dx/dy)^2}{EA}+\cfrac{y^2}{EI} \right ) ds=\int \limits_{\alpha_1}^{\alpha_2} \left ( \cfrac{(dx/dy)^2}{EA}+\cfrac{y^2}{EI} \right ) d \alpha= \cfrac{ \eta ^2 [ n \pi – \sin {(n \pi) ] } }{2EI}$,
$k_{12}= \int \limits_{\alpha_1}^{\alpha_2} \left ( \cfrac{(dx/dy) (dy/ds)} {EA}- \cfrac{xy}{EI} \right ) d \alpha = \cfrac{ \eta L [ n \pi \cdot \cos {(n \pi /2) } – 2 \sin{(n \pi / 2)]} }{n \pi \cdot EI }$,
$k_{22}= \int \limits_{\alpha_1}^{\alpha_2} \left ( \cfrac{(dy/ds)^2} {EA}- \cfrac{x^2}{EI} \right ) d \alpha = \cfrac{L}{12} \left ( \cfrac{6 \eta^2 [ n \pi +\sin{(n \pi) ] }}{EA}+\cfrac{L^2 \cdot n \pi}{EI}\right )$,
$k_{13}= \int \limits_{\alpha_1}^{\alpha_2} \cfrac{y}{EI} d \alpha =0$,
$k_{23}= \int \limits_{\alpha_1}^{\alpha_2} \cfrac{-x}{EI} d \alpha =0$,
$k_{33}= \int \limits_{\alpha_1}^{\alpha_2} \cfrac{1}{EI} d \alpha = \cfrac{n \pi}{EI}$,
Macierz sztywności w układzie (x,y)
Macierz sztywności uzyskamy przez odwrócenie macierzy podatności opisanej wyżej . Stosunkowo proste wyrażenie uzyskamy dla kształtu pręta proporcjonalnego do pierwszej postaci własnej, czyli dla n=1:
$$\begin{equation} [K] = \cfrac{1}{\Theta} \begin{bmatrix}
2 \pi^3 EI ( EA L^2/ \eta +6EI \cdot \eta & 48 \pi L \cdot EA \cdot EI & 0 \\
& 12 \pi^3 EA \cdot EI \cdot \eta & 0 \\
& SYM & EI \Theta / \pi \\
\end{bmatrix} \label {114} \end {equation} $$
gdzie: $\Theta= \eta [ EA L^2 (\pi^2 -96) +6EI \pi^4 \eta^2 ]$
Macierz sztywności ($\ref{114}) jest przydatna do badania wpływu imperfekcji elementu belki – słupa. Badanie tych zależności nie jest przedmiotem artykułu.
Przyklad 2 [ Belka Timoshenko na sprężystym podłożu ]
Rozwiązać belkę pokazaną na rys. 12 z uwzględnieniem wpływu sztywności ścinania na przemieszczenia i siły przekrojowe (belka Timoshenko).
Belka jest kratownicą Pratt’a (z wykratowaniem typu N) – wg tab. 3 poz. a) i danych:
- wysokość kratownicy h=720 mm,
- kąt nachylenia krzyżulców α=45º,
- pasy HEA 160, krzyżulce RK 40x40x3, słupki RK 50x50x3.
Zastosować macierzowe sformułowanie zagadnienia belki. Wykonać wykresy linii ugięcia oraz momentów zginających.
Dyskretyzacja oraz sztywność belki-kratownicy
Przyjęto podział belki na $N^{e} = 8$ elementów skończonych.
Na rys.12 węzły elementów wskazano wierzchołami trójkątów, w którego wnętrzu podano numer globalny węzła.
Liczba węzłów $w = 9$
W każdym węźle nietrywialne są dwa stopnie swobody $(w, \varphi)$ =(ugięcie, obrót) $q=2$
Macierze i wektory mają rozmiar
$m = w x q =9\cdot 2 = 18$
Długości elementów wynoszą:
[(1)-(2)], [(2)-(3)], [(7)-(8)], [(8)-(9)] 0,75 m,
a pozostałych elementów 1,00 m.
Belka ma stałą sztywność giętną oraz postaciową wszystkich elementów.
Pasy wykonano z profilu HEA160.
Moment bezwładności pasa ściskanego i rozciąganego wynosi $ I_c= I_t = 1670 \, cm^4$,
Pole przekroju $A_c = A_t =38,8 \, cm^2$.
Odległość pasa od osi obojętnej przekroju
$ z_c = z_t = h/2 = 72/2=36 \, cm$,
$A_k = 4,34 \, cm^2$ (krzyżulce RK 40x40x3),
$A_s = 5,54 \, cm^2$ (słupki RK 50x50x3).
Moment bezwładności na zginanie przekroju kratownicy, traktowanej jako belka o przekroju złożonym z pasów ( przekroju dwupunktowym) wynosi
$ I =I_c + A_c \cdot z_c^2 + I_t + A_t \cdot z_t^2 = 2\cdot (1670 + 38,8 \cdot 36^2)= 103909,6 \, cm^4 $
Moduł Younga stali wynosi E=210 GPa. W przykładzie nie uwzględniamy redukcji modułu Younga ($\ref {16}$), bowiem belka-kratownica ma swobodę przemieszczeń bocznych.
Sztywność giętna belki wynosi $EI= 210 GPa \cdot 103909,6 \, cm^4 \approx 2,182 \cdot 10^5 \, kNm^2 $
Sztywność postaciową belki-kratownicy obliczono z zależności podanej w poz a) tab.3 Po przeprowadzeniu stosownych przeliczeń, przyjęto $S_v = 25233,8 kN$
Macierze sztywności elementów
Z zależności ($\ref{100}$) uzyskano:
Elementy [(1)-(2)] i [(2)-(3)]
Elementy [(1)-(2)] i [(2)-(3)] mają długość L=0,75 m i są ściskane siłą N= 100 kN.
Współczynniki Λi, modyfikujące sztywność ($\ref{102}$) lub ($\ref{103}$) w tym przypadku wynoszą:
$\Lambda_i \approx 1,0$ (i=1,2,3,4).
Parametr $\Psi$, uwzględniający sztywność postaciową elementu [r] zdefiniowany wyrażeniem ($\ref{85}$) wynosi:
$ \Psi_{[(1)-(2)]}= \cfrac {2,182 \cdot 10^5} {0,75^2 \cdot 25233,81}=15,373 $
Dla elementów [e] =[(1)-(2)] oraz [(2)-(3)] macierze sztywności są takie same: $ k^{[(1)-(2)]} = k^{[(2)-(3)]}$.
Pierwszy element tych macierzy wynosi
$ k_{1,1} ^{[(1)-(2)]}= \cfrac {EI} {1+12 \Psi_{[(1)-(2)]}} \cdot \cfrac {55} {L_1^3} \cdot \Lambda_\Psi_{[(1)-(2)]} =\cfrac{2,182 \cdot 10^5} {1+12 \cdot 15,373} \cdot \cfrac{55} {0,75^2} \cdot 1,0=33462,825 $ kN/m
a pełna macierz sztywności :
$ [k]^{[(1)-(2)], [(2)-(3)]}= \small {{\left \lbrack \matrix { 33462,825 & 12548,829 & -33462,825 & 12548,829 \cr 12548,829 & 295650,171 & -12548,829 & -286242,279 \cr -33462,825 & -12548,829 & 33462,825 &-12548,829 \cr 12548,829 & -286242,279 & -12548,829 & 295650,171} \right \rbrack} } = \left | \matrix { {[k_{54}]} & {[k_{55}]} \cr {[k_{64}]} & {[k_{65}]} } \right |$
W celu przygotowania do „sumowania MES macierz sztywności elementu składamy z czterech bloków. Bloki macierzy sztywności elementu [(1)-(2)] i [(2)-(3)] oznaczamy bez nadkreślenia np.: $[k_{54}]$. Natomiast różniące się od nich bloki macierzy innych elementów będziemy oznaczać jednym lub dwoma nadkreśleniami, $[\overline k_{54}]$ lub $[ {\overset = k }_{54}]$.
Elementy [(3)-(4)] do [(6)-(7)]
Elementy nie są ściskane P(=N)=0 → Λi =0 (i=1,2,3,4). Macierz sztywności rozpatrywanych elementów wynosi więc:
$ [k]^{[(3)-(4)]…[(6)-(7)]}=\small{ {\left \lbrack \matrix { 24992,963 & 12496,481 & -24992,963 & 12496,481 \cr 12496,481 & 224458,401 & -12496,481 & -211961,919 \cr -24992,963 & -12496,481 & 24992,963 &-12496,481 \cr 12496,481 & -211961,919 & -12496,481 & 224458,401} \right \rbrack }} =\left | \matrix { {[\overline k_{54}]} & {[ \overline k_{55}]} \cr {[\overline k_{64}]} & {[ \overline k_{65}]} } \right |$
Elementy [(7)-(8)] i [(8)-(9)]
Macierz sztywności przy uwzględnieniu N=0 oraz długości L=L1=0,75 m , wynosi:
$ [k]^{[(7)-(8)], [(8)-(9)]}= \small {\left \lbrack \matrix { 33463,688 & 12548,883 & -33463,688 & 12548,883 \cr 12548,883 & 295652,711 & -12548,883 & -286241,049 \cr -33463,688 & -12548,883 & 33463,688 &-12548,883 \cr 12548,883 & -286241,049 & -12548,883 & 295652,711} \right \rbrack} = \left | \matrix { {[{\overset = k }_{54}]} & {[ {\overset = k }_{55}]} \cr {[{\overset = k }_{64}]} & {[ {\overset = k }_{65}]} } \right |$
Macierz sztywności układu
Macierz sztywności układu [K] ma wymiar (9×2)x(9×2)=18×18.
Mamy na przykład:
$ K_{54} = [k_{54}] = \small {\left \lbrack \matrix { +33462,825 & +12548,829 \cr +12548,829 & +295650,171} \right \rbrack} $
$ K_{55} = [k_{55}] = \small {\left \lbrack \matrix { -33462,825 & +12548,829 \cr -12548,829 & -286241,279} \right \rbrack} $
$K_{56}=K_{57}=K_{58}=K_{59}=0$
$ K_{65} = [k_{65}]+[k_{54}] = \small {\left \lbrack \matrix { +33462,825 & -12548,829 \cr -12548,829 & +295650,171} \right \rbrack +\left \lbrack \matrix { +33462,825 & +12548,829 \cr +12548,829 & +295650,171} \right \rbrack}$
$K_{69}=[k_{65}]+[\overline k_{54}]$, itd.
W tab.P2-1. zestawiono zagregowaną („zsumowaną”) macierz sztywności całego układu.
Tab.P2-1. Macierz sztywności systemu

Wektor równoważników węzłowych
Wektor równoważników węzłowych $$\ref{101}$) dla danych zadania ( po uwzględnieniu obciążenia międzywęzłowego q oraz reakcji od podłoża ze stałą sprężystości c – wektor równoważników węzłowych) zestawiono w tab .P2-2a i 2b
Tab. P2- 2a Wektor równoważników węzłowych od odporu podłoża oraz od podpory sprężystej
Tab. P2-2b Wektor równoważników węzłowych od obciążeń czynnych rozłożonych i skupionych
Uwaga: Funkcję kształtu przy obliczaniu wektora równoważników węzłowych elementów belki leżących na sprężystym podłożu, przyjęto jak dla elementu poddanego jedynie zginaniu, choć element ten jest również ściskany. Przybliżenie takie jest wystarczająco dokładne w obliczeniach praktycznych.
Po uwzględnieniu wartości stałych sprężystości c i C oraz obciążenia P i M oraz faktu, że q<0, P<0 (ujemne), otrzymamy wektor równoważników węzłowych zestawiony w tab.P1-2c.
TabP2-2c. Wektor równoważników węzłowych po przeniesieniu wyrazów od sprężystych podpór
Kinematyczne warunki brzegowe. Modyfikacja macierzy sztywności
Postępując wg algorytmu ($\ref{105}$) wprowadzimy warunek
$w_3=0$.
Wyrazy zawierające niewiadome przemieszczenia, występujące w wektorze obciążenia przenosi się na stronę macierzy sztywności.
Zmodyfikowane równanie kanoniczne ($\ref{103}$) dla prezentowanego przykładu podano w tab.P2-3.
Tab.P2-3. Zmodyfikowana macierz sztywności
Rozwiązanie równania kanonicznego
Przemieszczenia
W wyniku rozwiązania zmodyfikowanego, kanonicznego układu równań(Tab.3), otrzymano następujące wartości przemieszczeń ($w$ [m] $\varphi $ [rad] :
$ |u| = \small {\left | \matrix { w_1= 0,00000 \cr \varphi_1 – 3,80841 \cr w_2= 0,00000 \cr \varphi_2= -3,80871 \cr w_3=0,00000 \cr \varphi_3 =-3,80964 \cr w_4=-9,21016 \cr \varphi_4 = -4,21969 \cr w_5=-1,62009 \cr \varphi_5 = -4,17215 \cr w_6= -20,65211 \cr \varphi_6=-3,94199 \cr w_7= -22,51866 \cr \varphi_7 = -3,80417 \cr w_8= -2,85051 \cr \varphi_8= -4,19084 \cr w_9= -34,68487 \cr \varphi_9= -4,31973} \right |\cdot 10^{-4} }$
Wykresy linii ugięcia pokazano na rys. P2-1. Wykresy uzyskano poprzez połączenie liniami prostymi wartości węzłowych lub przywęzłowych. Nie uwzględniono rzeczywistej, najczęściej nieliniowej zmienności wykresu pomiędzy węzłami. To samo dotyczy wykresu momentów zginających belkę (rys. P2-2).
Siły przekrojowe
Siły przekrojowe obliczono, wykonując krok odwrotny wg zależności ($\ref{104}$) i ($\ref{105}$), uzyskując dla poszczególnych elementów skończonych:
$ |S|= [V_1 | M_1 | V_2 | M_2]^{(e)}=$
$\small [ -10,000 | 0,000 | 625,91 | -7,97]^{[(1)-(2)]} $
$\small [ -625,91 | 7,97| 1455,25 | -15,94]^{[(2)-(3)]}$
$\small [ 15,99 | 15,94 | -9,99 | -2,96]^{[(3)-(4)]}$
$\small [ 9,99 | 2,96 | -3,99 | 4,03]^{(4-5)}$
$\small [ 3,99 | 4,03 | 2,02 | 5,02]^{(5-6)}$
$\small [ -2,02 | -5,02 | 8,02 | 0,00]^{[(6)-(7)]}$
$\small [ 10,00 | 15,00 | -10,00 | -7,50]^{[(7)-(8)]}$
$\small [ 10,00 | 7,50 | -10,00 | 0,00]^{[(8)-(9)]}$
Wykresy momentów zginających pokazano na rys. P2-2.
Literatura
- Rodrigo F., N., Comparison of Beam Theories for Cross Sectional Deflection Parameter analysis of Euler, Timoshenko and Reddy beam theories, Bachelor Report, Eindhoven University of Technology, Department of Mechanical Engineering (Energy Technology), Eindhoven, July 1, 2023
- Szymczyk, M. (1981). Słownik języka polskiego. PWN
- Bańko, M. (Ed.). (2000). Inny słownik języka polskiego. PWN
- Wikipedia. (2013). Zginanie,[ https://pl.wikipedia.org/w/index.php?title=Zginanie&oldid=37876761]
- Timoshenko S., Goodier J.N., Teoria sprężystości, Arkady Warszawa 1962
- Birger, I. A., & Panowko, J. G. (1968). Procnost, ustojcivost, kolebanija, Spravocnik (Vol. 1). Mosinostroejnije
- Tablice matematyczne,, Mizerski W (red.), Wydawnictwo Adamantan, Warszawa 2002
- Piechnik, St., Wytrzymałość materiałów dla wydziałów budowlanych. PWN, Warszawa-Kraków, 1980
- Haukaas T., Euler-Bernoulli Beams, Lectures, The University of British Columbia, Vancouver, [ terje.civil.ubc.ca]
- Timoshenko, S. P., Goodier, J. N. (1969). Theory of Elasticity. McGraw-Hill
- Haukaas T., Timoshenko Beams, Lectures, The University of British Columbia, Vancouver, [ terje.civil.ubc.ca]
- Felippa, C. A., Introduction to Finite Element (ASEN 5007) . Part 13. Lecture for Curse Fall 2001, [http://www.colorado.edu/engineering/CAS/courses.d/Structures.d/IAST.Lect03.d/IAST.Lect03.pdf ]
- Pałkowski S., (2009), Konstrukcje stalowe: wybrane zagadnienia obliczania i projektowania. Wydawnictwo Naukowe PWN, Warszawa
- Zienkiewicz O.C., Taylor R., L., The Finite Element Method, Vol I. Basic Formulation and Linear Problems, Fourth Ed. McGraw -Hill Boook Company, 1989
- Norri D., de Vries G., An Introduction to Finite Element Method, Dep. of Mechanical Engineering , The University pf Calgary, Calgary, 1978
- Heyliger P., R., Reddy J.,N. , A higher order beam finite element for bending and vibration problems, Journal of Sound and Vibration, 126(2):309–326, 1988
- Livesley R. K. (1975), Matrix methods of structural analysis, Pergamon Press [http://books.google.com/books?id=tYsoAQAAMAAJ ]
________________________________













