Belka w inżynierii: Część III: Teoria Timoshenko i metoda elementów skończonych

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 [1]

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

Podstawienie do lewej strony ($\ref{72})$ wyrażenia na naprężenie ścinające, określone wzorem Żurawskiego ($\ref{32}$) i przyjętej zależności $\tau_v=\cfrac{V}{A_v} = \cfrac{V} {\beta A}$ prowadzi do równania:

$$ \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} $$

Rozwiązując ($\ref {74}$) ze względu na $\beta$ uzyskujemy

$$ \begin{equation} \beta = \cfrac{I_y^2}{A \cdot \int \limits_A \left(\cfrac{V}{t}\right)^2 dA} \label{75} \end{equation} $$

W przypadku przekrojów belek o prostym kształcie po wykonaniu całkowania przypisanego w ($\ref{75}$), otrzymano:

$$ \begin{equation} \beta =\begin {cases}
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 [2]

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 [2]:

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

Rys.11 Umowa znakowania V-$\gamma$

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 [3]

Felippa (2001) [2] 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) [3].

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 a_i. 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 [4], [5]. 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

Rys.12 Przykład belki Timoshenko. Podział na elementy skończone MES

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ć [2], ż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 [6] – 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 [7].

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{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\phi_1 – 3310151336182.58w_2 + 845723833.31\phi_2 \\ -847552471.02w_1 – 126996250.17\phi_1 – 845723833.31w_2 + 126995195.48\phi_2 \\ (-3310151336182.58w_1 – 845723833.31\phi_1 – 9035679573.68w_2 + 847552471.02\phi_2) + (-9035679573.68w_2 – 847552471.02\phi_2 – 3310151336182.58w_3 + 845723833.31\phi_3) \\ (845723833.31w_1 + 2102324.63126995195.48\phi_1 + 847552471.02w_2 – 126996250.17\phi_2) + (-847552471.02w_2 – 126996250.17\phi_2 – 845723833.31w_3 + 126995195.48\phi_3) \\ -3310151336182.58w_2 – 845723833.31\phi_2 – 9035679573.68w_3 + 847552471.02\phi_3 \\ 845723833.31w_2 + 2102324.63126995195.48\phi_2 + 847552471.02w_3 – 126996250.17\phi_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\phi_1 – 3310151336182.58w_2 + 845723833.31\phi_2 + (-10) \\ -847552471.02w_1 – 126996250.17\phi_1 – 845723833.31w_2 + 126995195.48\phi_2 + (0) \\ (-3310151336182.58w_1 – 845723833.31\phi_1 – 9035679573.68w_2 + 847552471.02\phi_2) + (-9035679573.68w_2 – 847552471.02\phi_2 – 3310151336182.58w_3 + 845723833.31\phi_3) + (0) \\ (845723833.31w_1 + 2102324.63126995195.48\phi_1 + 847552471.02w_2 – 126996250.17\phi_2) + (-847552471.02w_2 – 126996250.17\phi_2 – 845723833.31w_3 + 126995195.48\phi_3) + (0) \\ -3310151336182.58w_2 – 845723833.31\phi_2 – 9035679573.68w_3 + 847552471.02\phi_3 + (-3) \\ 845723833.31w_2 + 2102324.63126995195.48\phi_2 + 847552471.02w_3 – 126996250.17\phi_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

Literatura

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 ]

________________________________

12 cze 26

By : Leszek Chodor

Comments : Off

O autorze

* dr inż. Leszek Chodor. Architekt i Inżynier Konstruktor; Rzeczoznawca budowlany. Autor wielu projektów budowli, w tym nagrodzonych w konkursach krajowych i zagranicznych, a między innymi: projektu wykonawczego konstrukcji budynku głównego Centrum "Manufaktura" w Łodzi, projektu budowlanego konstrukcji budynku PSE w Konstancinie Bielawa, projektów konstrukcji "Cersanit" ( Starachowice, Wałbrzych, Nowograd Wołyński-Ukraina), projektu konstrukcji hali widowiskowo-sportowej Arena Szczecin Autor kilkudziesięciu prac naukowych z zakresu teorii konstrukcji budowlanych, architektury oraz platformy BIM w projektowaniu.

Belka w inżynierii: Część III: Teoria Timoshenko i metoda elementów skończonych