Belka w inżynierii: Część III: Teoria Timoshenko. Podstawy MES

Część III: Belka w inżynierii. Teoria Timoshenko oraz podstawy MES

Belka Timoshenko

Belka Timoshenko została zdefiniowana ok. 1956 roku w ramach teorii Timoshenko-Ehrenfest [1]. Timoshenko lojalnie odnotował kluczowy wkład Ehrenfesta w przypisach swoich pierwszych publikacji. 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ężeniami stycznymi 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 Bernoulliego.

Belka Timoshenko uwzględnia korektę pierwszego rzędu teorii zginania poprzecznego Bernoulliego.
Klasyczne założenie kinematyczne zesztywnienia zostało 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. III.10.

Zastosowanie belki Timoshenko w praktyce eliminuje kryterium podziału belek i tarcz (zob. definicja belki).

Współczynnik $\beta$ efektywności przekroju na ścinanie

Na rys. III.1 przedstawiono ugięcie zastępczego elementu belki o długości $dx$, obciążonego na czole naprężeniem stycznym $q_z=\tau$. Element w rzeczywistości odkształca się mniej więcej tak, jak pokazano na rys. III.1c, 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 pokazano na rys. III.1d.

Rys. III.1. Ilustracja do współczynnika efektywności przekroju [2]

W celu oszacowania średniego odkształcenia ścinającego $\gamma_v$ (rys. III.1d), przyjętego w definicji belki Timoshenko, porównamy pracę wszystkich odkształcających się włókien (rys. III.1a,b) z pracą całego przekroju pod działaniem siły poprzecznej $V$:

\[ \int\limits_A \frac{1}{2}\,dw\,\tau\,dA = \frac{1}{2}\,dw_v\,V \tag{III.1} \label{III.1} \]

Podstawiając do ($\ref{III.1}$) zależność $dw=\gamma\,dx$, wynikającą z rys. III.1b, oraz związek wynikający z prawa Hooke’a dla ścinania $\tau=G\gamma$,

gdzie moduł Kirchoffa określa zależność

\[ G=\frac{E}{2(1+\nu)} \tag{III.2} \label{III.2} \]

otrzymujemy:

\[ \int\limits_A \frac12 \left(\frac{\tau}{G}\right) dx\,\tau\,dA = \frac12 \left(\frac{\tau_v}{G}\right) dx\,V \tag{III.3} \label{III.3} \]

Należy zachować ostrożność przy interpretacji $\tau_v$. Nie jest to zwykłe średnie naprężenie ścinające otrzymane przez równomierne rozłożenie siły poprzecznej $V$ na całym przekroju. Naprężenie $\tau_v$ odnosi się do efektywnego pola ścinania przekroju $A_v$, zdefiniowanego zależnością

\[ A_v=\beta A \tag{III.4} \label{III.4} \]

gdzie $\beta$ jest współczynnikiem efektywności przekroju na ścinanie.

Podstawienie do lewej strony równania ($\ref{III.3}$) wzoru Żurawskiego na naprężenie ścinające (I.32) oraz zależności $\tau_v=\dfrac{V}{A_v}=\dfrac{V}{\beta A}$ prowadzi do równania

\[ \int\limits_A \frac12 \frac1G dx \left( \frac{V\,\overline{S}_y}{I_y\,t} \right)^2 dA = \frac12 \frac{V}{\beta A} \frac1G dx\,V \tag{III.5} \label{III.5} \]

Rozwiązując równanie ($\ref{III.5}$) względem $\beta$, otrzymujemy

\[ \beta= \frac{I_y^2} {A\displaystyle\int\limits_A \left( \frac{\overline S_y}{t} \right)^2dA} \tag{III.6} \label{III.6} \]

Po wykonaniu całkowania dla najczęściej spotykanych przekrojów otrzymuje się:

\[ \beta= \begin{cases} \dfrac56, & \text{dla przekroju prostokątnego},\\[2mm] \dfrac9{10}, & \text{dla przekroju kołowego},\\[2mm] \dfrac12, & \text{dla rury okrągłej},\\[2mm] \approx\dfrac{A_{\parallel}}{A}, & \text{dla przekrojów cienkościennych}. \end{cases} \tag{III.7} \label{III.7} \]

W przypadku przekrojów cienkościennych, złożonych z cienkich ścianek prostokątnych, $A_{\parallel}$ jest polem przekroju ścianek równoległych do kierunku działania siły poprzecznej. Przykładowo, dla dwuteownika ścinanego siłą pionową jest to pole przekroju środnika, dla dwuteownika ścinanego siłą poziomą – suma pól półek, natomiast dla rury prostokątnej – suma pól ścianek równoległych do kierunku działania siły poprzecznej.

Obrót przekroju belki na skutek zginania i ścinania

W celu wyznaczenia obrotu przekroju belki zginanej poprzecznie zastosujemy zasadę prac wirtualnych

\[ \delta W_{ext}=\delta W_{int} \tag{III.8} \label{III.8} \]

gdzie praca sił zewnętrznych wynosi

\[ \delta W_{ext} = \sum_{i=1}^{n_F} F_{ext,i}\,\Delta_{ext,i} \tag{III.9} \label{III.9} \]

gdzie:
$n_F$ – liczba wszystkich sił zewnętrznych.

Praca sił wewnętrznych ma postać

\[ \delta W_{int} = \sum_{[e]=1}^{N_{[e]}} \left( \delta W_{N,[e]} + \delta W_{M,[e]} + \delta W_{V,[e]} \right) \tag{III.10} \label{III.10} \]

gdzie:
$[e]$ – numer elementu,
$N_{[e]}$ – liczba elementów konstrukcji.

Praca sił wewnętrznych wywołana siłą osiową (z pominięciem efektów II rzędu) wynosi

\[ \delta W_{N,[e]} = \delta N \frac{NL}{EA} \tag{III.11} \label{III.11} \]

Praca sił wewnętrznych wywołana zginaniem wynosi

\[ \delta W_{M,[e]} = \int_0^L \delta M \frac{M}{EI} dx \tag{III.12} \label{III.12} \]

Natomiast praca sił wewnętrznych wywołana ścinaniem jest równa

\[ \delta W_{V,[e]} = \int_0^L \delta V \frac{V}{GA_v} dx \tag{III.13} \label{III.13} \]

Aproksymację rzeczywistego przemieszczenia ścinania przyjęto zgodnie z analizą przedstawioną w rozdziale Współczynnik efektywności przekroju na ścinanie.

Wnioskiem wynikającym z zależności ($\ref{III.8}$)–($\ref{III.13}$) jest podstawowa hipoteza teorii Timoshenko

\[ \varphi=\Theta+\gamma_v \tag{III.14} \label{III.14} \]

gdzie:
$\varphi$ – całkowity kąt obrotu przekroju (obrotowy stopień swobody elementu Timoshenko),
$\Theta=\dfrac{dw}{dx}$ – kąt obrotu osi belki według teorii Bernoulliego-Eulera (I.10),
$\gamma_v=\dfrac{\tau_v}{G}=\dfrac{V}{GA_v}$ – kąt obrotu wywołany odkształceniem postaciowym.

(Niżej zostanie pokazane, że zgodnie z przyjętą w mechanice budowli konwencją znaków dodatnia siła poprzeczna powoduje ujemny obrót ścinania $\gamma_v$.)

Ponieważ $\varphi\neq\Theta$, przekroje poprzeczne pozostają płaskie, lecz ich obrót nie jest już równy obrotowi osi belki (I.10). Oznacza to, że przekrój poprzeczny nie pozostaje prostopadły do odkształconej osi belki.

Powyższe zależności zilustrowano na rys. III.2.

Rys. III.2. Dwuwęzłowy element Timoshenko: a) przemieszczenia końców elementu, w tym $\varphi=\Theta+\gamma_v$; b) obrót $\Theta$ wywołany zginaniem; c) obrót $\gamma_v$ wywołany ścinaniem (na podstawie: [3]).

Pole przemieszczeń i smukłość ścinania

Pole przemieszczeń liniowych belki Timoshenko można zapisać w postaci

\[ u(x,z)=u(x)+z\varphi, \qquad w(x,z)=w(x) \tag{III.15} \label{III.15} \]

Zachodzi ponadto [3]

\[ \varphi=\Theta-\gamma_v, \qquad \gamma_v=\frac{\tau_v}{G} =\frac{V}{GA_v} =\frac{V}{S_v} \tag{III.16} \label{III.16} \]

gdzie:
$V$ – siła poprzeczna,
$\gamma_v=\gamma_v(x)$ – uśredniony po przekroju kąt ścinania,
$G$ – moduł Kirchoffa,
$A_v$ – efektywne pole przekroju na ścinanie,
$S_v=GA_v$ – sztywność postaciowa przekroju.

Ujemny znak przy $\gamma_v$ w równaniu ($\ref{III.16}$) wynika z przyjętej w mechanice budowli konwencji znaków, zgodnie z którą dodatnia siła poprzeczna wywołuje ujemny obrót ścinania (rys. III.3).

Rys. III.3. Umowa znakowania siły poprzecznej $V$ i kąta ścinania $\gamma_v$.

W przypadku przekrojów złożonych (wielogałęziowych) sztywność postaciową $S_v$ wyznacza się przez porównanie przemieszczeń konstrukcji rzeczywistej i pręta zastępczego. Dla najczęściej spotykanych prętów o pasach równoległych odpowiednie zależności zestawiono w tab. III.1.

Tab. III.1. Podatność na ścinanie prętów złożonych [4]

Funkcjonał Lagrange’a dla belki Timoshenko

W klasycznej teorii Timoshenko obrót przekroju $\varphi(x)$ traktowany jest jako dodatkowa zmienna wariacyjna, niezależna od przemieszczenia poprzecznego $w(x)$. Podejście takie jest wygodne przy wyprowadzaniu równań różniczkowych oraz funkcji kształtu elementów skończonych, nie jest jednak jedynym możliwym sformułowaniem zagadnienia.

Podobnie jak w funkcjonałach mieszanych, np. Hu-Washizu lub Hellingera-Reissnera, liczba niezależnych pól może zostać zwiększona w celu uproszczenia wyprowadzenia równań, jednak zmienne pomocnicze nie są wielkościami podstawowymi. Po wykorzystaniu związków konstytutywnych i równań równowagi mogą zostać wyeliminowane, prowadząc do równoważnego funkcjonału zapisanego wyłącznie za pomocą jednego pola przemieszczeń.

W niniejszej pracy przyjęto właśnie takie podejście. Obrót przekroju nie jest traktowany jako niezależna niewiadoma, lecz zostaje wyeliminowany z wykorzystaniem zależności pomiędzy momentem zginającym, siłą poprzeczną i odkształceniem postaciowym. Otrzymany funkcjonał zachowuje jednopolową postać charakterystyczną dla teorii Bernoulliego, natomiast wpływ odkształceń postaciowych uwzględniony jest przez dodatkowy składnik energii. Takie sformułowanie stanowi naturalny punkt wyjścia do wyprowadzenia uściślonych macierzy sztywności metodą energetyczną.

Energia potencjalna belki zginanej i ścinanej Timoshenko

Energia potencjalna belki Timoshenko, po uwzględnieniu odkształceń postaciowych wywołanych siłami poprzecznymi, stanowi uogólnienie energii potencjalnej belki Bernoulliego-Eulera (II.47). W przeciwieństwie do teorii Bernoulliego energia odkształcenia składa się z dwóch części: energii zginania oraz energii ścinania, co prowadzi do zależności

\[ \Pi_{\sigma,T}=\frac{1}{2}\int_0^L\left(EI_y\kappa^2(x)+S_v\gamma^2(x)\right)dx, \tag{III.17} \label{III.17} \]

gdzie:
$\kappa(x)=w”(x) $ krzywizna osi belki,
$ \gamma(x)$ uśredniony kąt odkształcenia postaciowego.

W celu zachowania jednopolowego opisu przemieszczenia elementu – kąt ścinania nie jest traktowany jako niezależna zmienna wariacyjna, lecz zostaje wyznaczony z zależności równowagi pomiędzy momentem zginającym, siłą poprzeczną oraz odkształceniem postaciowym otrzymuje się 

W celu zachowania jednopolowego opisu przemieszczenia elementu kąt ścinania nie jest traktowany jako niezależna zmienna wariacyjna, lecz zostaje wyznaczony z zależności równowagi pomiędzy momentem zginającym, siłą poprzeczną oraz odkształceniem postaciowym. Korzystając z zależności

$ M(x)=-EI_y\,w”(x),$  oraz równania równowagi $ V(x)=\frac{dM(x)}{dx},$  otrzymuje się $  V(x)=-EI_y\,w”'(x),$

a następnie

\[ \gamma(x)=\frac{V(x)}{GA_v} = -\frac{EI_y}{GA_v}\,w”'(x). \tag{III.18} \label{III.18} \]

Po podstawieniu zależności (\ref{III.18}) do funkcjonału (\ref{III.17}) energia potencjalna belki Timoshenko zostaje wyrażona wyłącznie przez funkcję ugięcia oraz jej pochodne

\[ \Pi_{\sigma,T} = \frac{EI_y}{2} \int_0^L \left( [w”(x)]^2 + \frac{EI_y}{GA_v}[w”'(x)]^2 \right) dx. \tag{III.19} \label{III.19} \]

Otrzymany funkcjonał zachowuje jednopolową postać charakterystyczną dla teorii Bernoulliego (II.47), natomiast wpływ odkształceń postaciowych uwzględniony jest przez dodatkowy składnik energii zawierający trzecią pochodną funkcji ugięcia.

Praca siły ściskającej , sił poprzecznych i podpór sprężystych

Ponieważ  przy wyprowadzaniu formuł dla belki Bernoulliego  (II.60) na pracę siły ściskającej (II.54) , pracę  sił poprzecznych  (II.51) czy odporu podpór sprężystych (II.59) NIE  wykorzystywaliśmy  hipotezy kinematycznej Bernoulli, więc zależności słuszne są również dla belki Timoshenko.

Praca siły ściskającej, obciążeń poprzecznych oraz podpór sprężystych zależy wyłącznie od przemieszczeń elementu oraz układu obciążeń zewnętrznych. Przy wyprowadzaniu odpowiednich składników funkcjonału Lagrange’a dla belki Bernoulliego (II.60) nie wykorzystywano hipotezy płaskich przekrojów ani zależności kinematycznych charakterystycznych dla teorii Bernoulliego. Z tego względu wyrażenia opisujące pracę siły osiowej, pracę obciążeń poprzecznych oraz energię podpór i podłoża sprężystego pozostają niezmienione również w teorii Timoshenko.

Oznacza to, że przejście od teorii Bernoulliego do teorii Timoshenko wymaga jedynie modyfikacji części funkcjonału opisującej energię odkształcenia elementu przez dodanie energii ścinania. Wszystkie pozostałe składniki funkcjonału zachowują swoją postać, co w dalszej części pracy umożliwi bezpośrednie uściślenie klasycznych macierzy sztywności metodą energetyczną, bez konieczności ich ponownego wyprowadzania od podstaw.

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:

\[  \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 ] \tag{III.20} \label{III.20} \]

Funkcjonał $\Pi_T$ dla belki Timoshenko różni się od funkcjonału belki Bernoulliego (II.60) jedynie dodatkowym składnikiem energii odkształceń postaciowych, zaznaczonym podkreśleniem we wzorze (\ref{III.20}). Wszystkie pozostałe składniki funkcjonału, opisujące energię zginania, pracę siły osiowej, obciążeń poprzecznych, podłoża sprężystego oraz podpór sprężystych, pozostają niezmienione. 

Wynik ten ma istotne znaczenie metodologiczne. Oznacza on, że przejście od teorii Bernoulliego do teorii Timoshenko nie wymaga ponownego wyprowadzania całego funkcjonału Lagrange’a. Wystarczy jedynie uzupełnić energię odkształcenia o składnik opisujący wpływ ścinania, zachowując niezmienioną postać wszystkich pozostałych wyrażeń.
Własność ta stanowi podstawę energetycznej metody konstruowania uściślonych macierzy sztywności przedstawionej w następnym rozdziale. Zamiast wyprowadzać od podstaw ścisłą macierz elementu Timoshenko z rozwiązania równań różniczkowych lub klasyczny element MES z dodatkowymi polami interpolacji, wystarczy odpowiednio uściślić klasyczną macierz Bernoulliego poprzez uwzględnienie dodatkowej energii odkształceń postaciowych. Takie podejście prowadzi do prostych i efektywnych obliczeniowo macierzy sztywności, zachowujących strukturę klasycznego elementu Bernoulliego, a jednocześnie uwzględniających wpływ ścinania z dokładnością wystarczającą dla zdecydowanej większości zastosowań inżynierskich.

Macierz sztywności elementu Timoshenko 

Jak wykazano w poprzednim rozdziale, funkcjonał Lagrange’a belki Timoshenko różni się od funkcjonału belki Bernoulliego jedynie dodatkowym składnikiem energii odkształceń postaciowych. Oznacza to, że przejście od teorii Bernoulliego do teorii Timoshenko nie wymaga ponownego wyprowadzania całej macierzy sztywności od podstaw. Wystarczy odpowiednio uściślić klasyczną macierz elementu Bernoulliego, zachowując jej strukturę, liczbę stopni swobody oraz podstawowe własności numeryczne.

Podstawą modyfikacji jest bezwymiarowy parametr smukłości ścinania $\Xi$,  wynikający wprost z analiz prowadzonych w poprzedniej części pracy i który można zapisać w postaci

\[  \Xi = \cfrac {12\cdot EI_y} {GA_v L^2}= \cfrac {12 \cdot EI_y} {S_v\cdot L^2} = 12 \cdot k \tag{III.21} \label{III.21} \]

Jeśli  $\Xi \to 0$,  to model Timoshenko redukuje się do klasycznego modelu Bernoulli-Eulera. Parametr $\Xi$ jest miarą wpływu odkształceń postaciowych na pracę elementu. Dla elementów smukłych zachodzi $\Xi\ll1,$. Natomiast dla elementów krótkich lub o małej sztywności postaciowej jego wartość wzrasta, zwiększając wpływ ścinania na odpowiedź konstrukcji. 
Współczynnik $\Xi$ pozostaje w ściślym ziwąku z innym inżynirskim prarametrem podatności  na ścinanie  $k$  wprowadzonym przez Pałkowskiego (2009) [5] i który jest cząsto stosowany przez inżynierów w meijsce podatności postaciowej pręta $S_v$ zgodnie z zależnością

\[  k = \cfrac{EI_y}{S_v \cdot L^2} \tag{III.22} \label{III.22} \]

Wyrażenia na podatności $S_v$  dla  często stosowanych belek kratownicowych lub słupów wielogałęziowych  zestawiono w Tab. III.1 . Zależą one od typu skratowania kąta nachyklenia krzyżulcó $\alpha$ oraz powierzchni przekroju pretów $A_k, A_s, $  oraz kilku innych paramterów zgodnie z zestawionymi zależnosciami.

Efektywne pole przekroju ścinania wyraża zależność (\ref{III.4}$  (A_v=\beta A)$, a współczynnik efektywności ścinania $\beta$ (\ref{III.7}).

Macierz sztywności elementu Timoshenko uzyskuje się poprzez bezpośrednią modyfikację klasycznej macierzy belki Eulera–Bernoulliego o z funkcjami statecznościowymi  $\Psi_i$, (II.24)  z użyciem smukłości ścinania $\Xi$ (\ref{21), i po zdefiniowaniu współczynnika redukcji sztywności zginania 

\[ \eta=\frac{1}{1+\Xi}. \tag{III.23} \label{III.23} \]

Element belki Timoshenko uzyskamy przez uściślenie elementu Bernoulliego, co prowadzi do praktycznej postaci elementu Timoshenko bez porzeby wyprowadzania ścisłeego-analitycznego elementu  Timoshenki (exact Timoshenko beam element z rozwiązania równań róniczkowych, ani też klasycznego elementu MES z  interpolacji funkcji kształtu.

W dalszej części pracy wykorzystywane będą trzy postacie macierzy elementu Timoshenki:

pełna – z korekcją η (\ref{III.23}) i składnikami  Ξ (\ref{III.21})
uproszczona – z korekcją części zginanej przez η 
praktyczna (Felippy) – z globalnym współczynnikiem  η dla całej macierzy.

Takie podejście umożliwia uzyskanie macierzy sztywności stosowane już w prajycznych algorytmach z ich niewielikimi modyfikacjami rachunkowymi , ale z zachoweaniem stryktury oraz liczby stopni swobody., ale także uzyskanie wyników wystraczająco dokładnych dla praktyki in żynierskiej.

Pełna uściślona macierz Timoshenko

Pełną uściśloną macierz otrzymuje się przez modyfikację współczynników klasycznej macierzy Bernoulliego zawierającej funkcje stateczności $\Psi_i$ (II.24). Modyfikacja obejmuje dwie operacje:

• redukcję wszystkich współczynników odpowiadających części zginanej przez współczynnik $\eta$, przy niezmienionej części osiowej.
• wprowadzenie  dodatkowych składników zależnych od parametru $\Xi$ w wyrazach odpowiadających sztywności obrotowej przekrojów. Modyfikacji nie podlegają same funkcje stateczności $\Psi_i$, lecz współczynniki macierzy, w których funkcje te występują. Odpowiednie podstawienia mają postać:

\[ \Psi_1 \longrightarrow \eta\Psi_1 = \frac{\Psi_1}{1+\Xi}, \tag{III.24} \label{III.24} \]

\[ \Psi_2 \longrightarrow \eta\Psi_2 = \frac{\Psi_2}{1+\Xi}, \tag{III.25} \label{III.25} \]

\[ 4\Psi_3 \longrightarrow \eta\left(4\Psi_3+\Xi\right) = \frac{4\Psi_3+\Xi}{1+\Xi}, \tag{III.26} \label{III.26} \]

\[ 2\Psi_4 \longrightarrow \eta\left(2\Psi_4-\Xi\right) = \frac{2\Psi_4-\Xi}{1+\Xi}. \tag{III.27} \label{III.27} \]

Po podstawieniu powyższych zależności do macierzy Bernoulliego (II.23) otrzymuje się macierz sztywnosći elementu Timoshenki  w postaci:

\[ [k]^{[e]}_{T}= \left[ \begin{array}{cccccc} \dfrac{EA}{L} & 0 & 0 & -\dfrac{EA}{L} & 0 & 0 \\[2mm]
& \dfrac{12EI_y}{L^3}\eta\Psi_1 & \dfrac{6EI_y}{L^2}\eta\Psi_2 & 0 & -\dfrac{12EI_y}{L^3}\eta\Psi_1 & \dfrac{6EI_y}{L^2}\eta\Psi_2 \\[2mm]
& & \dfrac{EI_y}{L}\eta\left(4\Psi_3+\Xi\right) & 0 & -\dfrac{6EI_y}{L^2}\eta\Psi_2 & \dfrac{EI_y}{L}\eta\left(2\Psi_4-\Xi\right) \\[2mm]
&&& \dfrac{EA}{L} & 0 & 0 \\[2mm] & \mathrm{SYM} & & & \dfrac{12EI_y}{L^3}\eta\Psi_1 & -\dfrac{6EI_y}{L^2}\eta\Psi_2 \\[2mm]
& & & & & \dfrac{EI_y}{L}\eta\left(4\Psi_3+\Xi\right) \end{array} \right] \tag{III.28}  \label{III.28}  \]

Zmianie ulegają wyłącznie współczynniki związane ze zginaniem i obrotem przekrojów. Część osiowa pozostaje identyczna jak w teorii Bernoulliego. Otrzymaną macierz można zatem interpretować jako uściślenie klasycznego elementu Bernoulliego poprzez uwzględnienie energii odkształceń postaciowych.

W granicznym przypadku $ \Xi\rightarrow0, \qquad \eta\rightarrow1, $ achodzi: $\eta\Psi_i\rightarrow\Psi_i, $, $ $\eta(4\Psi_3+\Xi)\rightarrow4\Psi_3,$ $\eta(2\Psi_4-\Xi)\rightarrow2\Psi_4, $ co oznacza, że pełna macierz Timoshenko przechodzi dokładnie w klasyczną macierz elementu Bernoulliego.

Uproszczona macierz Timoshenko

W praktycznych obliczeniach inżynierskich bardzo często analizowane są elementy smukłe, dla których wpływ odkształceń postaciowych jest niewielki. W takim przypadku parametr $Xi=\frac{12EI_y}{\kappa GA L^2}$ osiąga małe wartości $ \Xi\ll1,$m, co oznacza, że wyrazy korekcyjne występujące wyłącznie w członach opisujących sztywność obrotową mają wpływ drugiego rzędu.
Dla małych wartości parametru $\Xi$ można rozwinąć odpowiednie współczynniki w szereg Taylora $ \frac{4\Psi_3+\Xi}{1+\Xi}= 4\Psi_3+\left(1-4\Psi_3\right)\Xi+O(\Xi^2),$ $\frac{2\Psi_4-\Xi}{1+\Xi} = 2\Psi_4-\left(2\Psi_4+1\right)\Xi+O(\Xi^2). $.   

Ponieważ zarówno $\Xi$, jak i wyrazy rzędu $\Xi^2$ są niewielkie, ich wpływ na globalną odpowiedź elementu jest zazwyczaj znacznie mniejszy niż wpływ samej redukcji sztywności zginania opisanej współczynnikiem $\Xi$ (\ref {III,17}). Wówczas  pełną macierz elementu Timoshenki można z dużą dokładnością zastąpić uproszczoną postacią otrzymaną przez pominięcie składników korekcyjnych $\pm\Xi$ występujących w wyrazach obrotowych. Wówczas cała część odpowiadająca zginaniu ulega jedynie redukcji przez współczynnik $ \eta$ (\ref{III.23})  Macierz sztywności przyjmie postać

\[ [k]^{[e]}_{T,\,\mathrm{sm}} = [k]^{[e]}_{N} + \eta\,[k]^{[e]}_{M}, \tag{III.29}  \label{III.29}  \]

gdzie

\[ [k]^{[e]}_{N}= \frac{EA}{L} \left[ \begin{array}{cccccc} 1&0&0&-1&0&0
\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ -1&0&0&1&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0 \end{array} \right] \tag{III.30}  \label{III.30}  \]

jest macierzą sztywności osiowej, natomiast

\[ [k]^{[e]}_{M} = \frac{EI_y}{L} \left[ \begin{array}{cccccc} 0 & 0 & 0 & 0 & 0 & 0 \\[2mm]
0 & \dfrac{12\Psi_1}{L^2} & \dfrac{6\Psi_2}{L} & 0 & -\dfrac{12\Psi_1}{L^2} & \dfrac{6\Psi_2}{L} \\[2mm]
0 & \dfrac{6\Psi_2}{L} & 4\Psi_3 & 0 & -\dfrac{6\Psi_2}{L} & 2\Psi_4 \\[2mm]
0 & 0 & 0 & 0 & 0 & 0 \\[2mm]
0 & -\dfrac{12\Psi_1}{L^2} & -\dfrac{6\Psi_2}{L} & 0 & \dfrac{12\Psi_1}{L^2} & -\dfrac{6\Psi_2}{L} \\[2mm]
0 & \dfrac{6\Psi_2}{L} & 2\Psi_4 & 0 & -\dfrac{6\Psi_2}{L} & 4\Psi_3 \end{array} \right]. \tag{III.31}  \label{III.31}  \]

oznacza część zginaną klasycznej macierzy Bernoulliego zawierającą funkcje stateczności $\Psi_i$.

Taka postać zachowuje strukturę elementu Bernoulliego, a wpływ odkształceń postaciowych uwzględniany jest wyłącznie przez globalny współczynnik redukcyjny $\eta$. Przybliżenie to jest bardzo dokładne dla elementów smukłych i znacznie upraszcza implementację numeryczną, gdyż nie wymaga modyfikacji poszczególnych współczynników macierzy.

Uproszczenie nie jest stosowane dla  dla belek krótkich, głębokich lub wykonanych z materiałów o stosunkowo małym module ścinania i wóczas stosuje ć pełną macierz Timoshenki z zachowaniem dodatkowych składników $\pm\Xi$. W tych przypadkach odkształcenia postaciowe mają istotny wpływ na rozkład sił wewnętrznych, ugięcia oraz wartości własne zagadnienia stateczności.

Praktyczna macierz Timoshenko

W praktycznych implementacjach elementów Timoshenki często stosuje się jeszcze dalsze uproszczenie Felippa (2001)[6]  przedstawił taką postać w sposób szczególnie przejrzysty. Przyjął mianowicie, że z wystarczającą dla praktyki dokładnością współczynnik redukyjny $\eta$ (\ref{III.23}) można zastosować dla pełnej macierzy sztysności Bernullniego-Eulera (II.23) (bez jej podziału na część pohodząca od zginania i od rozciągania/ściskania).  Takie uproszczenie Felippa przedstawił dla klasycznej (liniowej) macierzy sztywności elementu belkowego. Może ono zostać w naturalny sposób rozszerzone również na uściślone macierze sztywności zawierające funkcje stateczności $\Psi_i$, ponieważ sprowadza się jedynie do wprowadzenia globalnego współczynnika redukcyjnego $\eta$   

Formalnie błąd tego przybliżenia jest pierwszego rzędu względem parametru $\Xi$. W praktyce jego zakres stosowalności jest jednak znacznie szerszy, ponieważ w typowych konstrukcjach budowlanych dominują elementy smukłe, dla których wpływ ścinania pozostaje niewielki.
Wobec tego z praktycznego punktu widzenia uproszczenie to ma znacznie szerszy zakres stosowalności niż wynikałoby jedynie z analizy matematycznej. W typowych konstrukcjach budowlanych dominują elementy smukłe, dla których sztywność postaciowa jest wielokrotnie większa od sztywności zginania, a zatem parametr $\Xi$ przyjmuje niewielkie wartości. Dotyczy to w szczególności belek walcowanych, dwuteowników, ceowników, profili zamkniętych.
Szczególnie duże uzasadnienie praktyczne uproszczenie to znajduje w modelowaniu prętów wielogałęziowych i kratownicowych za pomocą pręta zastępczego. Już samo zastąpienie rzeczywistej konstrukcji prętem ciągłym stanowi przybliżenie, którego wpływ na wyniki jest zazwyczaj większy niż błąd wynikający z pominięcia składników $\pm\Xi$. W takich zastosowaniach dominującym źródłem niedokładności jest idealizacja modelu konstrukcyjnego, a nie sposób uwzględnienia odkształceń postaciowych w macierzy elementu.

Z tego względu uproszczona macierz Felippy stanowi bardzo korzystny kompromis pomiędzy dokładnością a prostotą implementacji. Zachowuje identyczną strukturę klasycznej macierzy Bernoulliego, nie wymaga zmian algorytmów numerycznych  ani liczby stopni swobody, a uzyskiwana dokładność jest w zdecydowanej większości analiz inżynierskich w pełni wystarczająca.

Praktyczna macierz Felippy stanowi zatem kompromis pomiędzy dokładnością i prostotą implementacji. Zachowuje identyczną strukturę klasycznej macierzy Bernoulliego, nie wymaga zmian algorytmów montażu ani liczby stopni swobody elementu, a uzyskiwana dokładność jest wystarczająca dla zdecydowanej większości zastosowań inżynierskich. Pełną uściśloną macierz Timoshenko należy natomiast stosować przede wszystkim w analizie elementów krótkich, głębokich, wykonanych z materiałów o małym module ścinania oraz w zagadnieniach wymagających wysokiej dokładności wyznaczania ugięć, częstości drgań lub obciążeń krytycznych.

Pełna macierz Timoshenko powinna być natomiast stosowana przede wszystkim w analizie elementów krótkich, głębokich, wykonanych z materiałów o małym module ścinania oraz w zagadnieniach wymagających wysokiej dokładności wyznaczania ugięć, częstości drgań lub obciążeń krytycznych.

Układ kanoniczny równań Lagranga-Ritza dla belki Timoshenko

Układ kanoniczny rónań  Lagrangea-Ritza jest analogiczne jak dla belki Bernoulliego-Eulera jest uogólnieniem  zagadnienia przedstawionego w  części II artykulu. Minimalizacja funkjonału meyodą Ritza   wyrazęniem (II.64) w postaci analiztycznej i w (II.67) w zapisie macierzowym.

Macierzowe sformułowanie zagadnienia belki dla belki Timoshenko

Równanie kanoniczne

Sformułowanie macierzowe zagadnienia belki Timoshenko jedt uogólnieniem  standardowego ujęcia macierzowego dl abelki Bernoulliego. Uohólnienie polega na zastosowaniu macierzy sztywności elemntu Timoshenko  i innych elementów terorii belki Timoshenko przedstawionej w poprzednich rozdziałach. Zagadnienie belki Bernoulli jest zagadnieniem szczególnym zagadnienia belki Timoshenko.

Równanie kanoniczne (II.67)  w zapisie macierzowym zachowuje postać 

\[ [\delta] \{a\} =\{\Delta_F \} \tag{II.67a} \label{III.67a}\]

a zaisie  dostosowanym do metody elmentów skończonych postać 

\[ [K] \{q\}= \{P\} \tag{II.72} \label{II.72} \]

gdzie:

$[K]$ – macierz sztywności konstrukcji o wymiarze (mxm),
$\{ q\}$ – poszukiwane przemieszczenia węzłowe o wymiarze m
$\{P\}$ – wektor obciążenia o wymiarze m,
$m = w  \cdot sq $ – wymiar macierzy i wektorów
$w$- całkowita liczba węzłów konstrukcji,
$s$ –  liczba stopni swobody w węźle. W przypadku rozpatrywanej belki $s=2$.

W tym rozdziale stosuje się róenież 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 [7], [8].
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.

Macierz sztywności elementu i konstrukcji 

Dyskretyzacja konstrukcji

Obliczenia numeryczne prowadzi się na  konstrukcji podzielonej na elementy skończone 

Macierz sztywności elementu [e] blki Timoshenko  zdefiniowano wyrażeniem (\ref{60}).

Macierz sztywnosci konstrukcji 

Macierz sztywności  konstrukcji uzyskujemy drogą superpozycji – „sumowania” po elementach skończonych, według zasady:

\[  [K]=[K_{rs}]=[\sum [k_{ij}]^{[e]} ] \tag{III.99} \label{III.99} \]

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 \Xi } \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\Xi )}{L} \Lambda_3 & \cfrac{ –  6}{L^2} \Lambda_2 & \cfrac{2 (1-6 \Xi )}{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 \Xi )}{L}\Lambda_4 & \cfrac{-6} {L^2} \Lambda_2 & \cfrac{4 (1+3 \Xi}{L} \Lambda_3\\
 \label{III.100}\end{bmatrix} \]

gdzie $L=L_{[e]}$, $EI = EI_{y, [e]}$, $\Xi=\Xi{[e]}$ – długość L , sztywność giętna $EI_y$ oraz funkcja $\Xi$ (${III.21}$) 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 $\Xi$ ($\ref{III.21}$),
• 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 1.

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 [9] – 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 

Przyklad 1[ 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 $\Xi$,  uwzględniający sztywność postaciową elementu [r]  zdefiniowany  wyrażeniem ($\ref{III.16}$) wynosi:

$ \Xi_{[(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{2,182 \cdot 10^5} {1+12 \cdot 15,373} \cdot \cfrac{55} {0,75^3} \cdot 1,0 \approx 33462,825 \text{ 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=L₁=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

$$K = \begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|} \hline & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 & 16 & 17 & 18 \\ \hline 1 & 33462.83 & 12548.83 & -33462.83 & 12548.83 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 2 & 12548.83 & 295650.17 & -12548.83 & 286242.28 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 3 & -33462.83 & -12548.83 & 66925.65 & 0 & -33462.83 & 12548.83 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 4 & 12548.83 & 286242.28 & 0 & 591300.34 & -12548.83 & -286242.28 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 5 & 0 & 0 & -33462.83 & -12548.83 & 58455.79 & -52.35 & -24992.96 & 12496.48 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 6 & 0 & 0 & 12548.83 & 286242.28 & -52.35 & 520108.57 & -12496.48 & -211961.92 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 7 & 0 & 0 & 0 & 0 & -24992.96 & -12496.48 & 49985.93 & 0 & -24992.96 & 12496.48 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 8 & 0 & 0 & 0 & 0 & 12496.48 & -211961.92 & 0 & 448916.8 & -12496.48 & -211961.92 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 9 & 0 & 0 & 0 & 0 & 0 & 0 & -24992.96 & -12496.48 & 49985.93 & 0 & -24992.96 & 12496.48 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 10 & 0 & 0 & 0 & 0 & 0 & 0 & 12496.48 & -211961.92 & 0 & 448916.8 & -12496.48 & 211961.92 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 11 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -24992.96 & -12496.48 & 49985.93 & 0 & -24992.96 & 12496.48 & 0 & 0 & 0 & 0 \\ \hline 12 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 12496.48 & -211961.92 & 0 & 448916.8 & -12496.48 & 211961.92 & 0 & 0 & 0 & 0 \\ \hline 13 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -24992.96 & -12496.48 & 58456.65 & 52.4 & -33463.69 & 12548.88 & 0 & 0 \\ \hline 14 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 12496.48 & 211961.92 & 52.4 & 520111.11 & -12548.88 & 286241.05 & 0 & 0 \\ \hline 15 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -33463.69 & -12548.88 & 66927.38 & 0 & -33463.69 & 12548.88 \\ \hline 16 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 12548.88 & 286241.05 & 0 & 591305.42 & -12548.88 & -286241.05 \\ \hline 17 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -33463.69 & -12548.88 & 33463.69 & -12548.88 \\ \hline 18 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 12548.88 & 286241.05 & -12548.88 & 295652.71 \\ \hline \end{array}$$

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

$$|F| = \begin{pmatrix} V_1^0 \\ M_1^0 \\ V_2^0 \\ M_2^0 \\ V_3^0 \\ M_3^0 \\ V_4^0 \\ M_4^0 \\ V_5^0 \\ M_5^0 \\ V_6^0 \\ M_6^0 \\ V_7^0 \\ M_7^0 \\ V_8^0 \\ M_8^0 \\ V_9^0 \\ M_9^0 \end{pmatrix} = \begin{bmatrix} -9035679573.68w_1 – 847552471.02\Xi_1 – 3310151336182.58w_2 + 845723833.31\Xi_2 \\ -847552471.02w_1 – 126996250.17\Xi_1 – 845723833.31w_2 + 126995195.48\Xi_2 \\ (-3310151336182.58w_1 – 845723833.31\Xi_1 – 9035679573.68w_2 + 847552471.02\Xi_2) + (-9035679573.68w_2 – 847552471.02\Xi_2 – 3310151336182.58w_3 + 845723833.31\Xi_3) \\ (845723833.31w_1 + 2102324.63126995195.48\Xi_1 + 847552471.02w_2 – 126996250.17\Xi_2) + (-847552471.02w_2 – 126996250.17\Xi_2 – 845723833.31w_3 + 126995195.48\Xi_3) \\ -3310151336182.58w_2 – 845723833.31\Xi_2 – 9035679573.68w_3 + 847552471.02\Xi_3 \\ 845723833.31w_2 + 2102324.63126995195.48\Xi_2 + 847552471.02w_3 – 126996250.17\Xi_3 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ -c_{w7} \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}$$

Tab. P2-2b   Wektor równoważników węzłowych od obciążeń czynnych rozłożonych i skupionych

$$+ \, q \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \frac{1}{2} l_2 \\ \frac{1}{12} l_2^2 \\ \frac{1}{2} l_2 + \frac{1}{2} l_2 = l_2 \\ \frac{1}{12} l_2^2 – \frac{1}{12} l_2^2 = 0 \\ l_2 \\ 0 \\ l_2 \\ 0 \\ \frac{1}{2} l_2 \\ -\frac{1}{12} l_2^2 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} P \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ M \\ 0 \\ 0 \\ P \\ 0 \end{pmatrix} \begin{matrix} 1 \\ 1 \\ 2 \\ 2 \\ 3 \\ 3 \\ 4 \\ 4 \\ 5 \\ 5 \\ 6 \\ 6 \\ 7 \\ 7 \\ 8 \\ 8 \\ 9 \\ 9 \end{matrix}$$

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

Oto zapis wektora sił węzłowych $|F|$ z obrazka image_732db5.png w formacie MathJax, uwzględniający dodane wartości liczbowe:

$$|F| = \begin{bmatrix} -9035679573.68w_1 – 847552471.02\Xi_1 – 3310151336182.58w_2 + 845723833.31\Xi_2 + (-10) \\ -847552471.02w_1 – 126996250.17\Xi_1 – 845723833.31w_2 + 126995195.48\Xi_2 + (0) \\ (-3310151336182.58w_1 – 845723833.31\Xi_1 – 9035679573.68w_2 + 847552471.02\Xi_2) + (-9035679573.68w_2 – 847552471.02\Xi_2 – 3310151336182.58w_3 + 845723833.31\Xi_3) + (0) \\ (845723833.31w_1 + 2102324.63126995195.48\Xi_1 + 847552471.02w_2 – 126996250.17\Xi_2) + (-847552471.02w_2 – 126996250.17\Xi_2 – 845723833.31w_3 + 126995195.48\Xi_3) + (0) \\ -3310151336182.58w_2 – 845723833.31\Xi_2 – 9035679573.68w_3 + 847552471.02\Xi_3 + (-3) \\ 845723833.31w_2 + 2102324.63126995195.48\Xi_2 + 847552471.02w_3 – 126996250.17\Xi_3 + (-0.5) \\ -6 \\ 0 \\ -6 \\ 0 \\ -6 \\ 0 \\ (-8000w_7) + (-3) \\ 0.5 + 15 \\ 0 \\ 0 \\ -10 \\ 0 \end{bmatrix}$$

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

Oto macierz sztywności $\tilde{K}$ (prawdopodobnie po uwzględnieniu warunków brzegowych lub transformacji) z załączonego obrazu image_74457e.png, przygotowana w formacie MathJax:

$$\tilde{K} = \begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|} \hline & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 & 16 & 17 & 18 \\ \hline 1 & 35713036.51 & 47565019.85 & 51302719.75 & 45711284.48 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 2 & 47565019.85 & 27291900.34 & 45711284.48 & 27281437.76 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 3 & 51302719.75 & 45711284.48 & 71426073.01 & 0 & 0 & 45711284.48 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 4 & 45711284.48 & 27281437.76 & 0 & 54583800.68 & 0 & 27281437.76 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 5 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 6 & 0 & 0 & 45711284.48 & 27281437.76 & 0 & 27516358.74 & -12496.48 & -211961.92 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 7 & 0 & 0 & 0 & 0 & 0 & -12496.48 & 49985.93 & 0 & -24992.96 & 12496.48 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 8 & 0 & 0 & 0 & 0 & 0 & -211961.92 & 0 & 448916.8 & -12496.48 & -211961.92 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 9 & 0 & 0 & 0 & 0 & 0 & 0 & -24992.96 & -12496.48 & 49985.93 & 0 & -24992.96 & 12496.48 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 10 & 0 & 0 & 0 & 0 & 0 & 0 & 12496.48 & -211961.92 & 0 & 448916.8 & -12496.48 & -211961.92 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 11 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -24992.96 & -12496.48 & 49985.93 & 0 & -24992.96 & 12496.48 & 0 & 0 & 0 & 0 \\ \hline 12 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 12496.48 & -211961.92 & 0 & 448916.8 & -12496.48 & -211961.92 & 0 & 0 & 0 & 0 \\ \hline 13 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -24992.96 & -12496.48 & 66456.65 & 52.4 & -33463.69 & 12548.88 & 0 & 0 \\ \hline 14 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 12496.48 & -211961.92 & 52.4 & 520111.11 & -12548.88 & -286241.05 & 0 & 0 \\ \hline 15 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -33463.69 & -12548.88 & 66927.38 & 0 & -33463.69 & 12548.88 \\ \hline 16 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 12548.88 & -286241.05 & 0 & 591305.42 & -12548.88 & -286241.05 \\ \hline 17 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -33463.69 & -12548.88 & 33463.69 & -12548.88 \\ \hline 18 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 12548.88 & -286241.05 & -12548.88 & 295652.71 \\ \hline \hline \end{array}$$

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).

Rys. P2-1. Linia ugięcia belki-kratownicy z przykładu

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.

Rys. P2-2 Momenty zginające belki-kratownicy z przykładu

Przykład Podstawy MES

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ć [10], ż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ą:

\[ | F_i | = ( – c_i) \cdot L^2 \begin{bmatrix}
\cfrac {13 + 294 \Xi +1680 \Xi^2}{L} & \cfrac {11 + 231 \Xi + 1260 \Xi^2}{6} & \cfrac{3 \cdot ( 3 + 84 \Xi + 560 \Xi^2) } {2 \cdot L } & \cfrac{ 13 + 378 \Xi + 2520 \Xi^2} {12}\\
\cfrac {11 + 231 \Xi + 1260 \Xi^2}{6} & \cfrac{ 1 + 21 \Xi + 126 \Xi^2} {3} \cdot L  & \cfrac { 13 + 378 \Xi + 2520 \Xi^2} {12} & \cfrac{ 1 + 28 \Xi + 168 \Xi^2} {4} \cdot L \\
\cfrac {3 \cdot ( 3 + 84 \Xi + 560 \Xi^2) } {2\cdot L} & \cfrac { 13 + 378 \Xi + 2520 \Xi^2} {12} & \cfrac{ 13 + 294 \Xi + 1680\Xi^2} {L}  & \cfrac{11 + 231 \Xi + 11260 \Xi^2}{6} &\\
\cfrac {13 \cdot 378 \Xi + 2520  \Xi^2 } {12} & \cfrac { 1 + 28 \Xi + 168 \Xi^2} {4}\cdot L  & \cfrac{ 11 + 231 \Xi + 1280\Xi^2} {6}  & \cfrac{1 + 21 \Xi + 126 \Xi^2}{3} \cdot L &\\
\end{bmatrix}
\cdot \begin{vmatrix} w_1 \\ \varphi_1 \\ w_2 \\ \varphi_2 \end{vmatrix}
\tag{101} \label{III.101} \]

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:

\[  \overline F_r = -k_{rk}\cdot \delta + F_r \tag{102} \label{III.102} \]

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

\[    [ \widetilde K ]  \cdot |u|=|\overline F| \tag{103} \label{III.103} \]

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]}:

\[ |F|^{[e]} = [K]^{[e]} \cdot |u|^{[e]} +|F^0|^{[e]} \tag{104} \label{III.104} \]

czyli dla indeksów w numeracji lokalnej:

\[ \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]}
\tag{105} \label{III.105} \] 

W rezultacie uzyskamy wektory sił przywęzłowych | V  [kN], M [kNm] | dla poszczególnych elementów skończonych.

Przykład 2 [ 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 {\Xi} {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}$,
$\Xi=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}{\Xi}
\begin{bmatrix}
\cfrac{\Xi}{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} ) $,
$\Xi = – 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{III.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{III.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 [11].

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{III.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{III.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}{\Xi} \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 \Xi / \pi \\
\end{bmatrix} \label{III.114} \end {equation} $$

gdzie: $\Xi= \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.

Literatura

1. Timoshenko, S. P. (1956). Strength of Materials: Advanced Theory and Problems (Cz. 2, 3. wyd.). D. Van Nostrand Company
2. Haukaas T., Timoshenko Beams, Lectures, The University of British Columbia, Vancouver, terje.civil.ubc.ca
3. Felippa, C. A., Introduction to Finite Element (ASEN 5007), Part 13, University of Colorado, 2001
4. Pałkowski S. (2009), Konstrukcje stalowe. Wybrane zagadnienia obliczania i projektowania, PWN, Warszawa
5. Pałkowski S., (2009), Konstrukcje stalowe: wybrane zagadnienia obliczania i projektowania. Wydawnictwo Naukowe PWN, Warszawa
6. 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 ]
7. Zienkiewicz O.C., Taylor R., L., The Finite Element Method, Vol I. Basic Formulation and Linear Problems, Fourth Ed. McGraw -Hill Boook Company, 1989
8. Norri D., de Vries G., An Introduction to Finite Element Method, Dep. of Mechanical Engineering , The University pf Calgary, Calgary, 1978
9. 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
10. 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 ]
11. Livesley R. K. (1975), Matrix methods of structural analysis, Pergamon Press [http://books.google.com/books?id=tYsoAQAAMAAJ ]

________________________________