Belka w inżynierii. Część II: Nieliniowość geometryczna. Podstawy MES

Leszek Chodor, 26 października 2014
26.10.2025 – nadal naprawa po poważnej awarii portalu.
W przypadku nieczytelnych treści, proszę powiadomić: leszek@chodor.co

W ciągu ostatnich 24 godzin z artykułu korzystało 5 Czytelników

CZĘŚC II: NIELINIOWOŚĆ I IMPERFEKCJE geometryczne oraz podstawy MES

Wprowadzenie

Analiza nieliniowej ścieżki równowagi  konstrukcji  pod obciążeniem parametryzowanym

Przeprowadzono analizę ścieżki równowagi ściskanego pręta, obejmującą zakres przedkrytyczny, stan utraty stateczności Eulera oraz zakres pokrytyczny. Analiza zachowania pokrytycznego elementów ma istotne znaczenie praktyczne, ponieważ osiągnięcie obciążenia krytycznego nie oznacza jeszcze utraty zdolności konstrukcji do przenoszenia dalszych obciążeń. Po bifurkacyjnej utracie stateczności element przechodzi na pokrytyczną ścieżkę równowagi i nadal może przenosić wzrastające obciążenie, wykorzystując rezerwy nośności wynikające z redystrybucji sił wewnętrznych oraz nieliniowej pracy konstrukcji. Dopiero dalszy wzrost przemieszczeń i odkształceń prowadzi do osiągnięcia stanu granicznego nośności, najczęściej związanego z uplastycznieniem przekroju, powstaniem mechanizmu zniszczenia lub inną postacią utraty zdolności nośnej. Z tego względu analiza pokrytyczna stanowi niezbędne uzupełnienie klasycznej analizy stateczności, umożliwiając ocenę rzeczywistego zachowania konstrukcji po przekroczeniu obciążenia krytycznego.

W szczególności w konstrukcjach stalowych, aluminiowych oraz cienkościennych zakres pokrytyczny często decyduje o rzeczywistej rezerwie bezpieczeństwa elementu i nie powinien być pomijany podczas oceny jego pracy. Z tego względu współczesne metody obliczeniowe coraz częściej wykorzystują analizę nieliniową drugiego rzędu (GMNIA), umożliwiającą wyznaczenie pełnej ścieżki równowagi konstrukcji – od stanu przedkrytycznego, poprzez punkt bifurkacji, aż do osiągnięcia granicznej nośności. Takie podejście pozwala nie tylko wyznaczyć obciążenie krytyczne, lecz również określić rzeczywistą nośność elementu oraz jego zdolność do bezpiecznej pracy po utracie stateczności, co ma zasadnicze znaczenie dla racjonalnego i ekonomicznego projektowania konstrukcji.

Analizowana konstrukcja obciążana jest proporcjonalnie względem przyjętej konfiguracji obciążeń odniesienia \( \mathbf{F}_0 \). Dla dowolnego poziomu obciążenia \( \mathbf{F} \) zachodzi zależność

\[ \mathbf{F}=\Lambda\,\mathbf{F}_0, \tag{II.1}\label{II.1} \]

gdzie \( \Lambda \) oznacza mnożnik obciążenia. Przyjmuje się, że wszystkie schematy statyczne oraz dane liczbowe odnoszą się do konfiguracji obciążeń odniesienia \( \mathbf{F}_0 \). Przykładowo, jeżeli na schemacie statycznym podano siłę poziomą \(H=10\,\mathrm{kN}\), moment skupiony \(M=20\,\mathrm{kNm}\) lub obciążenie równomiernie rozłożone \(q=3\,\mathrm{kN/m}\), to cały układ tych obciążeń wraz z miejscami ich przyłożenia należy interpretować jako konfigurację odniesienia

\[ \mathbf{F}_0= \{10\,\mathrm{kN};\,20\,\mathrm{kNm};\,3\,\mathrm{kN/m}\}. \]

Aktualne wartości wszystkich obciążeń odpowiadające zadanemu mnożnikowi \( \Lambda \) wyznaczane są z zależności (\ref{II.1}).

Wartość krytyczna mnożnika obciążenia \( \Lambda_{cr} \), wyznaczana z liniowej analizy wyboczeniowej (LBA), odpowiada stanowi utraty stateczności konstrukcji. W literaturze oraz normach parametr ten oznaczany jest również symbolami \( \alpha_{cr} \) lub \( k_{cr} \). Wygodną miarą stopnia wykorzystania nośności wyboczeniowej jest bezwymiarowy parametr

\[ \bar{\Lambda}=\frac{\Lambda}{\Lambda_{cr}}, \tag{II.2}\label{II.2} \]

który osiąga wartość równą jedności w stanie krytycznym. Parametr \( \bar{\Lambda} \) stanowi zatem wygodną miarę stopnia wykorzystania nośności wyboczeniowej oraz zbliżenia konstrukcji do stanu krytycznego. Wszystkie zależności opisujące zachowanie elementów i konstrukcji formułowane są w funkcji parametru \( \bar{\Lambda} \). Dzięki temu uzyskane rozwiązania mają charakter uniwersalny i pozostają niezależne od bezwzględnych wartości obciążeń odniesienia, zachowując jedynie ich wzajemną konfigurację.

Uogólnienie teorii belki Bernoulliego–Eulera o ściskanie. Nieliniowości geometryczne

Klasyczna teoria zginania Bernoulliego–Eulera może zostać uogólniona na przypadek elementów poddanych jednoczesnemu zginaniu oraz ściskaniu osiowemu. W dalszej części rozważony zostanie szczególny przypadek konstrukcji złożonej z pojedynczego elementu prętowego ściskanego jedną siłą osiową \(F\). Nie ogranicza to ogólności rozważań dotyczących nieliniowości geometrycznych, natomiast pozwala przedstawić ich podstawowe własności w zwartej postaci analitycznej oraz ułatwia interpretację otrzymanych wyników.

Nieliniowa macierz sztywności belki-słupa (\ref{II.22}), wykorzystująca funkcje statecznościowe \( \Psi_i \), (\ref{II.23}) należy do klasycznych rezultatów teorii stateczności konstrukcji prętowych. Jej rozwój związany jest przede wszystkim z pracami Livesleya i Chandlera (1956) [1] oraz późniejszą systematyzacją teorii przeprowadzoną przez Przemienieckiego (1968) [2]. Od czasu opublikowania tych fundamentalnych prac macierz statecznościowa belki-słupa była wielokrotnie analizowana, weryfikowana oraz implementowana w programach wykorzystujących metodę sztywności i metodę elementów skończonych. Ponad sześćdziesięcioletnia historia jej stosowania w analizie drugiego rzędu konstrukcji prętowych stanowi najlepsze potwierdzenie zarówno jej poprawności teoretycznej, jak i wysokiej przydatności praktycznej.

Klasyczna teoria została rozszerzona o ścisły równoważny wektor obciążeń węzłowych oraz nowe elementy dyskretne umożliwiające modelowanie podatnych połączeń, podpór i więzów sprężystych. Rozszerzenia te zachowują pełną zgodność z klasycznym sformułowaniem Livesleya, jednocześnie znacząco zwiększając zakres jego praktycznych zastosowań. Dzięki temu możliwe jest modelowanie szerokiej klasy konstrukcji prętowych z uwzględnieniem nieliniowości geometrycznych oraz lokalnych podatności połączeń przy zachowaniu ścisłego opisu pracy elementu.

Równania równowagi konstrukcji nieliniowej

Równanie równowagi nieliniowego  układu konstrukcyjnego można zapisać w klasycznej postaci

\[ [\tilde{K}]\,\mathbf{q} = \mathbf{\tilde{P}}, \tag{II.3}\label{II.3} \]

gdzie:

\( [\tilde{K}] \) – zmodyfikowana globalna macierz sztywności układu, otrzymana poprzez złożenie macierzy sztywności poszczególnych elementów oraz uwzględnienie kinematycznych warunków brzegowych,

\( \mathbf{q} \) – wektor niewiadomych przemieszczeń i obrotów węzłowych,

\( \mathbf{\tilde{P}} \) – zmodyfikowany wektor równoważnych obciążeń węzłowych.

Rozwiązanie układu równań kanonicznych (\ref{II.3}) ma postać

\[ \mathbf{q} = [\tilde{K}]^{-1}\mathbf{\tilde{P}}, \tag{II.4}\label{II.4} \]

gdzie \( [\tilde{K}]^{-1} \) oznacza globalną macierz podatności układu.

Klasyczną liniową macierz sztywności elementu Bernoulliego–Eulera będziemy oznaczać symbolem \( [K]_B \), natomiast nieliniową macierz sztywności wykorzystującą funkcje statecznościowe symbolem \( [K]_L \) (Livesley).  Jeżeli nie będzie to prowadziło do niejednoznaczności, indeksy te mogą być pomijane.

Pomimo zastosowania nieliniowej macierzy sztywności elementu \( [K]_L \), zarówno przed uwzględnieniem warunków brzegowych (\ref{II.33}), jak i po ich uwzględnieniu (\ref{II.34}), równanie kanoniczne układu zachowuje postać liniową względem niewiadomego wektora przemieszczeń \( \mathbf{q} \). Wynika to z faktu, że nieliniowość geometryczna elementu została ujęta w funkcjach statecznościowych \( \Psi_i \), które dla ustalonej wartości siły osiowej \(N\) (lub równoważnie – odpowiadającego jej parametru stateczności \( \mu \)) stanowią znane współczynniki macierzy sztywności. Oznacza to, że dla każdej ustalonej wartości siły osiowej układ równań pozostaje liniowy i może być rozwiązany klasycznymi metodami algebry liniowej. Zmiana wartości siły osiowej powoduje jedynie zmianę funkcji statecznościowych \( \Psi_i \), a tym samym aktualizację współczynników macierzy sztywności elementu. Dzięki temu analiza geometrycznie nieliniowa może być prowadzona jako ciąg kolejnych rozwiązań liniowych odpowiadających kolejnym poziomom obciążenia.

Przedstawione sformułowanie stanowi jedną z najważniejszych zalet ścisłej macierzy sztywności belki-słupa. Pozwala ono zachować prostotę klasycznej metody sztywności i metody elementów skończonych, przy jednoczesnym uwzględnieniu wpływu nieliniowości geometrycznych na pracę konstrukcji. W rezultacie otrzymuje się model obliczeniowy łączący ścisłość opisu mechanicznego z wysoką efektywnością numeryczną, co ma szczególne znaczenie w analizie drugiego rzędu oraz podczas śledzenia ścieżki równowagi konstrukcji.

Model elementu w konstrukcji

Model elementu $[e]$ z węzłami (i), (j) w konstrukcji modelowanej  podatnymi elementami węzłowymi Ξ , nazywany dalej interaktywami  pokazano  na rys. II.1.

Przedstawiony na rys. II.1 model elementu belki-słupa stanowi podstawowy obiekt dalszych rozważań. W odróżnieniu od klasycznego ujęcia, w którym warunki brzegowe opisuje się za pomocą idealnych podpór, analizowany element połączony jest z otaczającą konstrukcją za pośrednictwem interaktywów $\mathbf{\Xi}$.

Interaktywem nazywa się uogólniony model opisujący mechaniczne współoddziaływanie analizowanego elementu z jego otoczeniem. W najprostszym przypadku liniowym interaktywy mogą być reprezentowane przez sztywności translacyjne i obrotowe odpowiadające podatnym podporom lub podatnym połączeniom. W ujęciu bardziej ogólnym interaktywy opisują dowolne elementy dyskretne lub modele zastępcze otaczającej konstrukcji, których oddziaływanie zostało sprowadzone do węzłów analizowanego elementu. Mogą one reprezentować między innymi podatność fundamentów, podatność połączeń konstrukcyjnych, wpływ sąsiednich elementów, a także zastępczy model całej pozostałej części konstrukcji.
Pojęcie interaktywu nie ogranicza się do modeli liniowo-sprężystych. W zależności od przyjętego modelu mechanicznego może on opisywać zjawiska nieliniowe, kontakt jednostronny, tarcie, tłumienie, degradację właściwości mechanicznych, procesy uszkodzenia oraz odnowy elementów konstrukcji. Dzięki temu interaktyw stanowi uniwersalny model współoddziaływania elementu z jego otoczeniem i tworzy wspólną podstawę opisu klasycznych zagadnień stateczności, nieliniowości geometrycznej oraz dalszych analiz niezawodności konstrukcji.

W dalszych rozważaniach analizowany będzie szczególny przypadek elementu belki-słupa poddanego działaniu stałej siły ściskającej $N$. Założenie to jest spełnione dla prostych układów prętowych obciążonych siłami skupionymi przyłożonymi w węzłach konstrukcji i umożliwia uzyskanie ścisłych rozwiązań analitycznych. Jednocześnie nie ogranicza ono ogólności wyprowadzanych zależności, ponieważ parametr bezwymiarowy $\bar{\Lambda}$ został wcześniej zdefiniowany dla dowolnej konfiguracji obciążeń odniesienia. Dla analizowanego elementu siła osiowa jest stała na całej długości pręta  $  N=\mathrm{const}, $ a odpowiadające jej obciążenie krytyczne określa klasyczna teoria Eulera

\[ F_{cr}=N_{cr}=\frac{\pi^2EI}{L^2}. \tag{II.5} \label{II.5} \]

W tym szczególnym przypadku parametr $\bar{\Lambda}$ przyjmuje prostą postać

\[ \bar{\Lambda} =\frac{F}{F_{cr}} =\frac{N}{N_{cr}}, \tag{II.6} \label{II.6} \]

co stanowi bezpośrednią konsekwencję zależności ogólnej (\ref{II.2}). Oznacza to, że parametr $\bar{\Lambda}$ może być interpretowany jako stosunek aktualnej siły ściskającej do klasycznego obciążenia krytycznego Eulera.

Parametr $\bar{\Lambda}$ pełni rolę podstawowego parametru teorii belki-słupa. Wszystkie zależności opisujące funkcje kształtu, macierze sztywności, równoważne wektory obciążeń oraz rozwiązania zagadnień nieliniowych wyrażane będą jako funkcje parametru $\bar{\Lambda}$. Dzięki temu otrzymane rozwiązania zachowują postać bezwymiarową i mogą być stosowane niezależnie od wymiarów geometrycznych elementu oraz jego właściwości materiałowych.

na rys rys. II.2a przedstawiono lokalny układ współrzędnych elementu oraz przyjęte dodatnie zwroty przemieszczeń i sił węzłowych. Wektor stopni swobody elementu przyjmuje postać:

$ \{q\}^{[e]}= [  u_1 ,\, w_1 ,\,  \varphi_1,\,  u_2,\, w_2,\, \varphi_2 ]^T $ , a odpowiadający mu wektor sił węzłowych $ \{F\}^{[e]}= [  N_1 ,\, V_1 ,\,  M_1,\,  N_2,\, V_2,\, M_2 ]^T $

Przyjęta konwencja znaków obowiązuje w całej dalszej części pracy.

Na rys. II.2 przedstawiono też  dwa podstawowe modele referencyjne wykorzystywane w dalszej części pracy do wyprowadzenia ścisłych funkcji kształtu, macierzy sztywności oraz równoważnych wektorów obciążeń elementu belki-słupa: 

• Pierwszy model (rys. II.2b) przedstawia belkę swobodnie podpartą, obciążoną równomiernie rozłożonym obciążeniem poprzecznym $q=-p$ oraz ściskaną stałą siłą osiową $F$. Prawy węzeł elementu połączony jest z otaczającą konstrukcją za pośrednictwem interaktywu translacyjnego, opisanego bezwymiarowym parametrem $\bar{C}_{\Delta}=C_{\Delta}L^{3}/EI$. Model ten stanowi podstawę wyprowadzenia ścisłego równoważnego wektora obciążeń oraz analizy wpływu ściskania osiowego na przemieszczenia i siły wewnętrzne elementu.
• Drugi model (rys. II.2c) przedstawia ściskany słup wspornikowy poddany działaniu poprzecznej siły skupionej $H$. W analogiczny sposób mogą być analizowane również inne rodzaje obciążeń poprzecznych, w tym obciążenia rozłożone. Podstawa słupa połączona jest z otaczającą konstrukcją za pośrednictwem interaktywu obrotowego, opisanego bezwymiarowym parametrem $\bar{C}_{\varphi}=C_{\varphi}L/EI$. Model ten wykorzystany zostanie do wyprowadzenia ścisłych zależności opisujących wpływ podatności obrotowej na zachowanie ściskanego elementu oraz do analizy efektów drugiego rzędu i utraty stateczności.

Oba modele stanowią szczególne przypadki ogólnego modelu elementu belki-słupa przedstawionego na rys. II.1. Zastosowanie interaktywów umożliwia ujęcie w jednolitym formalizmie zarówno klasycznych warunków podporowych, jak i bardziej złożonych modeli współoddziaływania elementu z otaczającą konstrukcją.

Ścisła macierz sztywności elementu belki-słupa

Nieliniową macierz sztywności ściskanego elementu prętowego [e], pokazanego na rys. II.1a, wyznaczymy przez rozwiązanie równania różniczkowego belki-słupa bez obciążenia poprzecznego $q_z=0$ (I.33), przedstawionego w I części artykułu, uzupełnionego o wyraz geometrycznej nieliniowości wynikający z działania osiowej siły ściskającej $N$.

\[ 0= EI \cfrac{d^4w}{dx^4} + N\cfrac{d^2w}{dx^2} \tag{II.7} \label{II.7} \]

Wprowadzamy bezwymiarową współrzędną osi pręta $ \xi=\cfrac{x}{L}$  oraz parametr ściskania $ \mu^2=\cfrac{NL^2}{EI}$

Ponieważ dla rozważanego przypadku zachodzi $\bar{\Lambda} = \cfrac{N}{N_{cr}}$ oraz $ N_{cr} =\cfrac{\pi^2EI}{L^2}$, to parametr ściskania można zapisać w postaci $\mu^2 = \pi^2\bar{\Lambda}$

czyli

\[ \mu = \pi\sqrt{\bar{\Lambda}}. \tag{II.8} \label{II.8}\]

W granicznym stanie krytycznym zachodzi $\bar{\Lambda}=1$, a zatem $\mu=\pi$. Po podstawieniu powyższych zależności do równania (\ref{II.7}) otrzymujemy jego postać bezwymiarową

\[ \cfrac{d^4w}{d\xi^4} + \mu^2 \cfrac{d^2w}{d\xi^2} = 0 \tag{II.9} \label{II.9} \]

której rozwiązaniem jest funkcja ugięcia

\[ w(\xi) = C_1 + C_2\mu\xi + C_3\cos(\mu\xi) + C_4\sin(\mu\xi) \tag{II.10} \label{II.10} \]

oraz funkcja obrotów

\[ \cfrac{dw(\xi)}{d\xi} = \mu \left[ C_2 – C_3\sin(\mu\xi) + C_4\cos(\mu\xi) \right] \tag{II.11} \label{II.11} \]

Warunki brzegowe przyjmują postać

\[ w(0)=w_1, \qquad w(1)=w_2, \qquad \cfrac{dw(0)}{d\xi}=\varphi_1L, \qquad \cfrac{dw(1)}{d\xi}=\varphi_2L.\tag{II.12} \label{II.12} \]

i mogą zostać zapisane w postaci macierzowej

\[ [D] \begin{Bmatrix} C_1\\ C_2\\ C_3\\ C_4 \end{Bmatrix} = \begin{Bmatrix} w_1\\ \varphi_1L\\ w_2\\ \varphi_2L \end{Bmatrix} \tag{II.13} \label{II.13}  \]

gdzie

\[ [D] = \begin{bmatrix}
1 & 0 & 1 & 0\\
0 & \mu & 0 & \mu\\
1 & \mu & \cos\mu & \sin\mu\\
0 & \mu & -\mu\sin\mu & \mu\cos\mu
\end{bmatrix} \tag{II.14}  \label{II.14} \]

której wyznacznik wynosi

\[ \det[D] = \mu^2 (2\beta-\varepsilon). \tag{II.15} \label{II.15} \]

gdzie:
$ \begin{array}{ll}
\beta=1-\cos\mu, & \text{(II.13a)}\\
\varepsilon=\mu\sin\mu. & \text{(II.13b)}
\end{array} $

Po odwróceniu macierzy D otrzymamy 

\[ [D]^{-1}=
\cfrac{1}{\mu(2\beta-\varepsilon)} \begin{bmatrix}
\mu\gamma & \delta & -\mu\beta & \eta \\[2mm]
-\varepsilon & \beta & \varepsilon & \beta \\[2mm]
-\mu\beta & -\delta & \mu\beta & -\eta \\[2mm]
\varepsilon & \gamma & -\varepsilon & -\beta
\end{bmatrix}.\tag{II.16} \label{II.16}\]

gdzie
$\beta$ (II.13a}, $\varepsilon$ (II.13b), oraz:
$ \begin{array}{ll}
\gamma=1-\cos\mu-\mu\sin\mu, & \text{(II.16a)}\\
\delta=\mu\cos\mu-\sin\mu, & \text{(II.16b)}\\
\eta=\mu-\sin\mu. & \text{(II.16c)} \end{array}$

Siły przekrojowe wyznaczymy z zależności:

\[ M = EIw^{(II)} = -\cfrac{EI}{L^2} \cfrac{d^2w}{d\xi^2} = \mu^2 \cfrac{EI}{L^2} \left[ C_3\cos(\mu\xi) + C_4\sin(\mu\xi) \right] \tag{II.17} \label{II.17} \]

\[ V = M^{(I)} = -\cfrac{EI}{L^3} \cfrac{d^3w}{d\xi^3} = \mu^3 \cfrac{EI}{L^3} \left[ -C_3\sin(\mu\xi) + C_4\cos(\mu\xi) \right] \tag{II.18} \label{II.18} \]

Równania (\ref{II.17}) oraz (\ref{II.18}) określają rozkład momentów zginających i sił poprzecznych wzdłuż elementu w funkcji stałych całkowania $C_3$ i $C_4$. Stałe te nie są niezależne, lecz wynikają z warunków brzegowych (\ref{II.11}). 

Przy znajomości odwróconej macierzy $[D]$ zależność (\ref{II.11}) można zapisać w postaci

\[ \begin{Bmatrix} C_1\\ C_2\\ C_3\\ C_4 \end{Bmatrix} = [D]^{-1} \begin{Bmatrix} w_1\\ \varphi_1L\\ w_2\\ \varphi_2L \end{Bmatrix}. \tag{II.19} \label{II.19} \]

Podstawienie $\xi=0$ oraz $\xi=1$ do zależności (\ref{II.17})–(\ref{II.18}) prowadzi bezpośrednio do sił przekrojowych na końcach elementu:
$ M_1=M(\xi=0), \qquad V_1=V(\xi=0), \qquad M_2=M(\xi=1), \qquad  V_2=V(\xi=1). $

Wartości te odpowiadają siłom wewnętrznym działającym na przekroje końcowe elementu. Aby zachować zgodność z konwencją dodatnich sił węzłowych przyjętą na rys. II.1a, należy uwzględnić zmianę znaków siły poprzecznej i momentu na prawym końcu elementu. Ostatecznie wektor sił węzłowych części zginanej przyjmuje postać

\[   \{F_b\}\} ^{[e]}= \{ V(\xi=0) , ,\ M(\xi=0), ,\ -V(\xi=1), ,\ -M(\xi=1) \} \label{II.20} \tag{II.20} \]

Zmiana znaków wynika z przejścia od sił przekrojowych do odpowiadających im sił węzłowych i zapewnia symetrię macierzy sztywności wynikającą z twierdzenia Maxwella-Bettiego. Korzystając z trzeciego i czwartego wiersza macierzy odwrotnej (\ref{II.16}) wyrażamy stałe całkowania $C_3$ i $C_4$ przez stopnie swobody elementu. Po uporządkowaniu otrzymanych zależności względem przemieszczeń i obrotów węzłowych uzyskujemy równanie kanoniczne elementu

\[ \{F\}^{[e]} = [k]^{[e]} \{q\}^{[e]}. \tag{II.21} \label{II.21} \]

Po wyeliminowaniu stałych całkowania oraz przekształceniu powyższych zależności otrzymujemy ścisłą macierz sztywności geometrycznie nieliniowego elementu belki-słupa  

\[ [k]^{[e]} = \left[\begin{array}{cccccc}
\dfrac{EA}{L} & 0 & 0 & -\dfrac{EA}{L} & 0 & 0 \\
& \dfrac{12EI}{L^3}\Psi_1 & \dfrac{6EI}{L^2}\Psi_2 & 0 & -\dfrac{12EI}{L^3}\Psi_1 &\dfrac{6EI}{L^2}\Psi_2 \\
& & \dfrac{4EI}{L}\Psi_3 & 0 & -\dfrac{6EI}{L^2}\Psi_2 & \dfrac{2EI}{L}\Psi_4 \\
&&& \dfrac{EA}{L}&0&0\\
& \mathrm{SYM} & & & \dfrac{12EI}{L^3}\Psi_1 & -\dfrac{6EI}{L^2}\Psi_2\\
& & & & & \dfrac{4EI}{L}\Psi_3 \end{array} \right] \tag{II.22} \label{II.22}\]

Współczynniki \(\Psi_i\) są bezwymiarowymi funkcjami statecznościowymi opisującymi wpływ siły osiowej na sztywność zginania elementu i wynoszą:

\[ \begin{cases}
\Psi_1=\alpha_c\Psi_2\\
\Psi_2=\cfrac{\alpha^2}{3(1-\alpha_c)}\\
\Psi_3=\cfrac{3}{4}\Psi_2+\cfrac{\alpha_c}{4}\\
\Psi_4=\cfrac{3}{2}\Psi_2-\cfrac{\alpha_c}{2}
\end{cases}  \tag{II.23} \label{II.23} \]

\[ \alpha_c=\alpha\operatorname{ctg}\alpha,
\qquad \alpha=\cfrac{\mu}{2} =\cfrac{\pi}{2}\sqrt{\bar{\Lambda}}.
\tag{II.24}  \label{II.24} \]

Zachodzą następujące tożsamości funkcji statecznościowych:

\[\begin{array}{l}
4\Psi_1\Psi_3-3\Psi_2^{\,2}
=3(\alpha_c-1)\Psi_2^{\,2}+\alpha_c^{\,2}\Psi_2,\\
2\Psi_3-\Psi_4=\alpha_c,\\
2\Psi_3+\Psi_4=3\Psi_2,\\
4\Psi_3=3\Psi_2+\alpha_c.
\end{array} \tag{II.25}\label{II.25} \]

Dla rozważanego w artykule przypadku pojedynczego pręta zachodzi $ \bar{\Lambda}=\cfrac{N}{N_{cr}}$

W dalszej części przeprowadzona zostanie analiza asymptotyczna funkcji statecznościowych oraz wyznaczone zostaną współczynniki rozwinięć Taylora w funkcji parametru $\bar{\Lambda}$. Pozwoli to bezpośrednio porównać ścisłą teorię belki-słupa z kolejnymi przybliżeniami stosowanymi w analizie drugiego rzędu i teorii imperfekcyjnej.

Analiza asymptotyczna funkcji statecznościowych

W przypadku znajomości siły osiowej $N$ funkcje statecznościowe (\ref{II.24}) pozwalają wyznaczyć ścisłą sztywność elementu, a w konsekwencji również ścisłą sztywność geometrycznie nieliniową całego układu konstrukcyjnego. W klasycznych algorytmach MES siły osiowe są wielkościami poszukiwanymi, dlatego funkcje $\Psi_i$ wprowadzają nieliniowość i konieczność iteracyjnego rozwiązywania układu równań konstrukcji. W tym celu wygodnie jest rozwinąć funkcje statecznościowe w szereg Taylora względem parametru obciążenia. Klasyczne rozwinięcie jest podawane względem kwadratu parametru ściskania $\mu^2$ w postaci

\[  \begin{equation} \Psi_i \approx \sum_{j=0}^{4} C_{i,j} (\mu^2)^j \tag{II.26} \label{II.26} \end{equation} \]

gdzie:
$(i=1,2,3,4)$ – numer funkcji statecznościowej,
$(j=0,1,2,3,4)$ – numer wyrazu szeregu

Po pozostawieniu dwóch pierwszych wyrazów rozwinięcia w pobliżu $\alpha=0$ otrzymujem

\[  \begin{cases}
\Psi_1 = 1-\cfrac{1}{10}\mu^2,\\
\Psi_2 = 1-\cfrac{1}{60}\mu^2,\\
\Psi_3 = 1-\cfrac{1}{30}\mu^2,\\
\Psi_4 = 1+\cfrac{1}{60}\mu^2.
\end{cases} \tag{II.27} \label{II.27} \] 
Wyrażenia (\ref{II.27}) stanowią pierwsze niezerowe przybliżenie ścisłych funkcji statecznościowych. Po podstawieniu ich do macierzy (\ref{II.22}) otrzymuje się klasyczną macierz geometryczną elementu prętowego stosowaną w analizie drugiego rzędu.

W analizach stateczności konstrukcji bardziej użyteczne jest przedstawienie funkcji statecznościowych jako funkcji globalnego mnożnika obciążenia $\bar{\Lambda}=\cfrac{\Lambda}{\Lambda_{cr}}$dla którego w rozważanym przypadku zachodzi zależność  $ \mu^2=\pi^2\bar{\Lambda}.$

W zależności od liczby zachowanych członów rozwinięcia (\ref{II.26}) mówimy o stopniu analizy MES. Człony zerowego rzędu ($j=0$) odpowiadają analizie I rzędu. Człony pierwszego rzędu ($j=1$) odpowiadają klasycznej analizie II rzędu (LBA), często nazywanej analizą P–Δ. Uwzględnienie kolejnych członów rozwinięcia prowadzi do analiz wyższych rzędów.

Rozwinięcia funkcji statecznościowych bezpośrednio względem mnożnika obciążenia konstrukcji $\bar{\Lambda}$ można zapisać w postaci

\[ \Psi_i \approx \sum_{j=0}^{4} D_{i,j}\,\bar{\Lambda}^{j}. \tag{II.28} \label{II.28}  \]

Ponieważ zachodzi zależność $\mu^2=\pi^2\bar{\Lambda}$, to współczynniki obu rozwinięć są związane relacją:

\[ D_{i,j}=C_{i,j}\pi^{2j}. \tag{II.29} \label{II.29} \]

Współczynniki $C_{i,j}$ odnoszą się do rozwinięcia względem parametru ściskania $\mu^2$, natomiast współczynniki $D_{i,j}$ do rozwinięcia względem globalnego mnożnika obciążenia $\bar{\Lambda}$. Funkcje statecznościowe mogą być interpretowane bezpośrednio jako funkcje stopnia wykorzystania nośności wyboczeniowej konstrukcji. Parametr $\bar{\Lambda}$ posiada jednoznaczną interpretację fizyczną i jest naturalnym parametrem stosowanym w analizie statecznościowej, analizie imperfekcyjnej oraz procedurach oceny nośności konstrukcji.

Współczynniki rozwinięcia można wyznaczyć przez rozwinięcie funkcji (\ref{II.24}) w szereg Taylora w punkcie $\bar{\Lambda}=0$. Otrzymane szeregi stanowią lokalne przybliżenie ścisłych funkcji statecznościowych i pozwalają analizować wpływ kolejnych rzędów nieliniowości geometrycznej na sztywność elementu. 
W tab. II.1 zestawiono współczynniki rozwinięcia Taylora funkcji statecznościowych $\Psi_i$ względem stopnia wykorzystania nośności wyboczeniowej $\bar{\Lambda}$. Dla pojedynczego pręta rozwinięcia te są równoważne rozwinięciom względem parametrów $\mu^2=\pi^2\bar{\Lambda}$ oraz $\alpha=\cfrac{\pi}{2}\sqrt{\bar{\Lambda}}$.

Tab. II.1. Współczynniki rozwinięcia Taylora funkcji statecznościowych $\Psi_i$ względem potęg mnożnika obciażenie $\bar{\Lambda}$

\[ \begin{array}{|c|c|c|c|c|c|}
\hline\text{ Typ analizy \to}  & LA & LBA & SOA & TOA& FOA\\
\hline \text{ Rząd teorii\to} & 0 & 1 & 2 & 3 & 4\\
\hline \Psi_1 & 1 & -\dfrac{\pi^2}{10} & -\dfrac{\pi^4}{8400}
& -\dfrac{\pi^6}{75600} & -\dfrac{37\pi^8}{2328480000} \\
\hline \Psi_2 & 1 & -\dfrac{\pi^2}{60} & -\dfrac{\pi^4}{8400} & -\dfrac{\pi^6}{75600} & -\dfrac{37\pi^8}{2328480000}\\
\hline \Psi_3 & 1 & -\dfrac{\pi^2}{30} & -\dfrac{11\pi^4}{25200} & -\dfrac{\pi^6}{108000} & -\dfrac{509\pi^8}{2328480000} \\
\hline \Psi_4 & 1 & +\dfrac{\pi^2}{60} & +\dfrac{13\pi^4}{25200} & +\dfrac{11\pi^6}{756000} & +\dfrac{907\pi^8}{2328480000}\\
\hline \end{array} \]

Uwagi:

(1) Skróty typów analizy oraz ich typowe obszary zastosowania:

• LA (Linear Analysis) – klasyczna analiza liniowa (teoria zerowego rzędu), stosowana dla konstrukcji o małych przemieszczeniach, gdy wpływ nieliniowości geometrycznej można pominąć.
• LBA (Linear Buckling Analysis) – liniowa analiza stateczności (wyboczeniowa), służąca do wyznaczania obciążeń krytycznych oraz wykorzystywana jako podstawa analiz typu P–Δ.
• SOA (Second-Order Analysis) – analiza drugiego rzędu, stanowiąca minimalny poziom teorii wymagany do opisu zachowania konstrukcji cięgnowych oraz elementów o istotnych efektach geometrycznej nieliniowości.
• TOA (Third-Order Analysis) – analiza trzeciego rzędu, zalecana do modelowania tkanin technicznych oraz konstrukcji cięgnowo-membranowych, w których występują silne efekty geometrycznie nieliniowe.
• FOA (Fourth-Order Analysis) – analiza czwartego rzędu, przeznaczona dla konstrukcji szczególnie wrażliwych na nieliniowości geometryczne, zwłaszcza cienkich membran i bardzo wiotkich powłok.

(2) Rozwinięcie Taylora jest poprawne w otoczeniu punktu \(\bar{\Lambda}=0\). Dokładność przybliżenia maleje wraz ze wzrostem parametru \(\bar{\Lambda}\), zwłaszcza w pobliżu stanu rytycznego. W granicy \(\bar{\Lambda}\rightarrow 1\) zaleca się stosowanie pełnych funkcji stateczności \(\Psi_i\) lub ścisłej macierzy sztywności elementu.
(3) Przedstawiona interpretacja zachowuje ważność zarówno dla pojedynczego
elementu prętowego, jak i dla wieloelementowych modeli MES, pod warunkiem że
parametr \(\bar{\Lambda}\) wyznaczono na podstawie globalnego zagadnienia
stateczności konstrukcji.

Wniosek:

Rozwinięcie Taylora funkcji stateczności względem bezwymiarowego mnożnika obciążenia $\Lambda$ tworzy naturalny pomost pomiędzy klasyczną analizą liniową (LA), liniową analizą stateczności (LBA), analizą drugiego rzędu (SOA) oraz ścisłą teorią belki-słupa. Kolejne wyrazy szeregu można interpretować jako sukcesywne przybliżenia geometrycznie nieliniowej sztywności elementu, odpowiadające coraz wyższym rzędom teorii. Parametr $=bar \Lambda$  posiada przy tym jednoznaczną interpretację fizyczną jako stopień wykorzystania nośności wyboczeniowej konstrukcji. Dzięki temu stanowi on naturalną zmienną opisującą poziom nieliniowości geometrycznej oraz wygodny parametr zarówno w analizie stateczności, jak i w teorii imperfekcyjnej.

Macierz  sztywności elementu w  konstrukcji 

Ścisła macierz sztywności elementu belki-słupa (\ref{II.22}) została wyprowadzona dla elementu swobodnego i opisuje pełne właściwości mechaniczne pojedynczego elementu. Stanowi ona punkt wyjścia do analizy konstrukcji prętowych metodą sztywności oraz metodą elementów skończonych.

Po włączeniu elementu do modelu konstrukcji (rys. II.1) jego końce współpracują z sąsiednimi elementami oraz podporami tworzącymi otoczenie konstrukcyjne. W niniejszej pracy wpływ tego otoczenia opisano za pomocą interaktywów węzłowych $\Xi^{(i)}$ oraz $\Xi^{(j)}$, reprezentujących uogólnione właściwości połączeń elementu z konstrukcją. Interaktyw może modelować między innymi przeguby, więzy kinematyczne, podpory sprężyste, podatności translacyjne i obrotowe, a także bardziej złożone modele połączeń.

W rezultacie równanie elementu

\[\{F\}^{[e]}=[k]^{[e]}\{q\}^{[e]} \tag{II.30} \label{II.30} \]

ulega modyfikacji i przyjmuje postać

\[ \{\tilde{F}\}^{[e]}= [\tilde{k}]^{[e]} \{\tilde{q}\}^{[e]}. \tag{II.31} \label{II.31} \]

Macierz $ [\tilde{k}]^{[e]} $ nazywana będzie dalej  zmodyfikowaną macierzą sztywności elementu. Symbol tyldy oznacza, że została ona otrzymana ze ścisłej macierzy sztywności (\ref{II.22}) po uwzględnieniu rzeczywistych warunków współpracy elementu z otoczeniem konstrukcyjnym. Ogólne określenie modyfikacja macierzy sztywności, obejmuje wszystkie przekształcenia prowadzące od ścisłej macierzy elementu do modelu odpowiadającego rzeczywistym warunkom pracy konstrukcji. Modyfikacja może polegać na zmianie warunków brzegowych, redukcji lub kondensacji stopni swobody, dołączeniu elementów podatnych lub sprężystych, wprowadzeniu przegubów, a także rozbudowie modelu o dodatkowe stopnie swobody. Określenia redukcja i kondensacja odnoszą się jedynie do wybranych operacji algebraicznych, natomiast termin <i>modyfikacja</i> obejmuje wszystkie przekształcenia wykonywane podczas budowy modelu elementu konstrukcyjnego.

Należy podkreślić, że ścisła macierz sztywności (\ref{II.22}) opisuje element swobodny, a zatem zawiera stopnie swobody odpowiadające ruchom bryły sztywnej. W konsekwencji jest macierzą osobliwą i nie posiada macierzy odwrotnej, co uniemożliwia bezpośrednie wyznaczenie macierzy podatności elementu. Podstawowym celem dalszych rozważań jest takie zmodyfikowanie ścisłej macierzy sztywności, aby usunąć jej osobliwość i otrzymać nieosobliwą macierz uniwersalnego elementu konstrukcyjnego $[\tilde{k}]^{[e]},$  dla której istnieje macierz odwrotna  $[\tilde{k}]^{-1},$

Macierz $[\tilde{k}]^{-1},$ nazywana jest dalej uniwersalną macierzą podatności elementu konstrukcyjnego Wyznaczenie tej macierzy stanowi zasadniczy cel niniejszego rozdziału, ponieważ umożliwia analityczne wyprowadzenie zależności opisujących przemieszczenia, obroty oraz siły wewnętrzne dla szerokiej klasy elementów konstrukcyjnych. Wszystkie szczególne przypadki elementów otrzymuje się następnie przez odpowiedni dobór parametrów interaktywów węzłowych $\Xi^{(i)}$ i $\Xi^{(j)}$, bez konieczności ponownego wyznaczania macierzy odwrotnej.

Przykładowo, po uwzględnieniu sztywności translacyjnych i obrotowych interaktywów węzłowych $\Xi^{(i)}$ oraz $\Xi^{(j)}$, macierz sztywności interaktywów można zapisać w postaci

\[[k_{\Xi}]^{[e]}= \operatorname{diag} \left\{ 0,\, C_{\Delta}^{(i)},\, C_{\varphi}^{(i)},\, 0,\, C_{\Delta}^{(j)},\, C_{\varphi}^{(j)} \right\}. \tag{II.32} \label{II.32} \]

Uniwersalną zmodyfikowaną macierz sztywności elementu konstrukcyjnego otrzymujemy przez dodanie macierzy interaktywów do ścisłej macierzy sztywności elementu

\[[\tilde{k}]^{[e]} = [k]^{[e]} + [k_{\Xi}]^{[e]}= \left[ \begin{array}{cccccc} \dfrac{EA}{L} & 0 & 0 &-\dfrac{EA}{L} & 0 & 0 \\
0 & \dfrac{12EI}{L^3}\Psi_1+C_{\Delta}^{(i)} & \dfrac{6EI}{L^2}\Psi_2 & 0 & -\dfrac{12EI}{L^3}\Psi_1 & \dfrac{6EI}{L^2}\Psi_2 \\
0 & \dfrac{6EI}{L^2}\Psi_2 & \dfrac{4EI}{L}\Psi_3+C_{\varphi}^{(i)} & 0 &-\dfrac{6EI}{L^2}\Psi_2 &\dfrac{2EI}{L}\Psi_4 \\
-\dfrac{EA}{L} & 0 & 0 & \dfrac{EA}{L} & 0 & 0 \\
0 & -\dfrac{12EI}{L^3}\Psi_1 & -\dfrac{6EI}{L^2}\Psi_2 & 0 & \dfrac{12EI}{L^3}\Psi_1+C_{\Delta}^{(j)} & -\dfrac{6EI}{L^2}\Psi_2 \\
0 & \dfrac{6EI}{L^2}\Psi_2 & \dfrac{2EI}{L}\Psi_4 & 0 & -\dfrac{6EI}{L^2}\Psi_2 & \dfrac{4EI}{L}\Psi_3+C_{\varphi}^{(j)}
\end{array} \right]. \tag{II.33} \label{II.33} \]

gdzie $C_{\Delta}^{(i)}$ i $C_{\Delta}^{(j)}$ oznaczają sztywności translacyjne, natomiast $C_{\varphi}^{(i)}$ i $C_{\varphi}^{(j)}$ sztywności obrotowe interaktywów węzłowych $\Xi^{(i)}$ oraz $\Xi^{(j)}$, opisujących współpracę elementu z otoczeniem konstrukcyjnym.

Macierz podatności elementu w konstrukcji

Nieliniowa macierz sztywności elementu przed modyfikacją (\ref{II.33}) i po modyfikacji (\ref{II.38})  nie zaburza  liniowości układu kanonicznego konstrukcji (\re{55}), bowiem nieliniowośc jest zawarta w funkcjach statecznosciowych $\Psi_i$, które są traktowane jako parametry (dla ustalonej siły osiowej w elemencie $N$)

Macierz podatności elementu uzykuje się z odwrócenia (\ref{38}), a po wykorzystaniu tożsamości Livesleya oraz wprowadzeniu pomocniczych oznaczeń można ją zapisać w zwartej postaci:

\[\tilde{\mathbf K}^{-1}= \begin{bmatrix} \dfrac{EA+C_{nj}L}{D_N} &0&0&\dfrac{EA}{D_N}&0&0\\[2mm]
0&\dfrac{1}{C_{di}+C_{dj}} &-\dfrac{NC_1}{D_{C1}} &0 &\dfrac{1}{C_{di}+C_{dj}}&-\dfrac{NC_2}{D_{C1}}\\[2mm]
0 &-\dfrac{NC_1}{D_{C1}} &\dfrac{L\,NC_5}{D_{C2}} &0 &\dfrac{NC_3}{D_{C1}} &\dfrac{L\,NC_6}{D_{C2}} \\[2mm]
\dfrac{EA}{D_N}&0&0&\dfrac{EA+C_{ni}L}{D_N}&0&0\\[2mm]
0 &\dfrac{1}{C_{di}+C_{dj}} &\dfrac{NC_3}{D_{C1}} &0 &\dfrac{1}{C_{di}+C_{dj}} &\dfrac{NC_4}{D_{C1}} \\[2mm]
0 &-\dfrac{NC_2}{D_{C1}} &\dfrac{L\,NC_6}{D_{C2}} &0 &\dfrac{NC_4}{D_{C1}} &\dfrac{L\,NC_5}{D_{C2}}
\end{bmatrix} \tag{II.34} \label{II.34}  \]

gdzie:
$ D_N=(C_{ni}+C_{nj})EA+C_{ni}C_{nj}L, $
$ D_{C1}=24(ac-1)\,ac\,(C_{di}+C_{dj})EI^{2}\Psi_{2},$
$ D_{C2}=12(ac-1)\,ac\,EI\,\Psi_{2}, $
$ NC_1 = C_{dj}L^{2}\left(2ac\,EI+C_{fj}L\right), $
$ NC_2 = C_{dj}L^{2}\left(2ac\,EI+C_{fi}L\right), $
$NC_3 = C_{di}L^{2}\left(2ac\,EI+C_{fj}L\right), $
$ NC_4 = C_{di}L^{2}\left(2ac\,EI+C_{fi}L\right), $
$ NC_5 = ac^{2}+3(ac-1)\Psi_{2}, $
$ NC_6 = ac^{2}+3(1-ac)\Psi_{2}. $

W powyższych zależnościach parametr \(ac=\alpha\cot\alpha\) jest funkcją stateczności zdefiniowaną we wzorze (\ref{II.24}) $ \alpha=\frac{\pi}{2}\sqrt{\bar{\Lambda}},$. Funkcje stateczności \(\Psi_1,\Psi_2,\Psi_3,\Psi_4\) (\ref{II.23}) wyrażone są również przez ten  parametr.

Ścisłe funkcje kształtu nieliniowej belki-słupa

Po rozwiązaniu równania różniczkowego belki-słupa, nałożeniu warunków węzłowych (\ref{II.11}), wyznaczeniu stałych całkowania \(C_i\) oraz podstawieniu ich do rozwiązania ogólnego (\ref{II.10}), funkcję ugięcia można zapisać w postaci interpolacyjnej

\[ w(\xi)= \begin{bmatrix} N_1(\xi) & N_2(\xi) & N_3(\xi) & N_4(\xi) \end{bmatrix} \begin{Bmatrix} w_1\\ \varphi_1\\ w_2\\ \varphi_2 \end{Bmatrix}, \tag{II.35} \label{II.35} \]

gdzie \(N_i(\xi)\) są ścisłymi funkcjami kształtu Livesleya odpowiadającymi klasycznym stopniom swobody dwuwęzłowego elementu belkowego. Parametr

\[ \mu=L\sqrt{\cfrac{N}{EI}}=\pi\sqrt{\bar{\Lambda}} \tag{II.36} \label{II.36} \]

jest bezwymiarowym parametrem stateczności elementu, natomiast \(\bar{\Lambda}=N/N_E\) oznacza znormalizowaną siłę osiową odniesioną do siły krytycznej Eulera. Funkcje kształtu mają postać

\[ N_1(\xi)= \cfrac{ 1-\cos\mu+\cos(\mu\xi)-\cos[\mu(1-\xi)] -\mu(1-\xi)\sin\mu }{\Theta}, \tag{II.37} \label{II.37} \]
\[ N_2(\xi)= L\, \cfrac{ \mu\xi+\mu(1-\xi)\cos\mu -\mu\cos[\mu(1-\xi)] -\sin\mu+\sin(\mu\xi)+\sin[\mu(1-\xi)] }{\Theta}, \tag{II.38} \label{II.38} \]
\[ N_3(\xi)= \cfrac{ 1-\cos\mu-\cos(\mu\xi)+\cos[\mu(1-\xi)] -\mu\xi\sin\mu }{\Theta}, \tag{II.39} \label{II.39} \]
\[ N_4(\xi)= L\, \cfrac{ -\mu(1-\xi)-\mu\xi\cos\mu +\mu\cos(\mu\xi) +\sin\mu-\sin(\mu\xi)-\sin[\mu(1-\xi)] }{\Theta}, \tag{II.40} \label{II.40} \]

gdzie:
\[ \Theta= 2(1-\cos\mu)-\mu\sin\mu. \tag{II.41} \label{II.41} \]

W granicy \(\mu\rightarrow0\) funkcje kształtu (\ref{II.35}) przechodzą dokładnie w klasyczne wielomianowe funkcje Hermite’a stosowane w liniowej teorii belek. Oznacza to, że ścisły element belki-słupa stanowi naturalne uogólnienie klasycznego elementu Bernoulliego na przypadek prętów obciążonych siłą osiową.

Ścisły równoważny wektor obciążeń węzłowych nieliniowej belki-słupa

Rozkład przemieszczeń opisany ścisłymi funkcjami kształtu (\ref{II.35}) umożliwia bezpośrednie wyznaczenie ścisłego równoważnego wektora obciążeń węzłowych od dowolnego obciążenia międzywęzłowego działającego na element. Zgodnie z zasadą prac wirtualnych wektor sił węzłowych otrzymuje się przez rzutowanie rozkładu obciążenia na funkcje kształtu elementu

\[ \{P\}^{[e]}= \int_0^L [\mathbf N]^T(x)\, p(x)\, dx, \tag{II.42} \label{II.42} \]

gdzie \(p(x)\) oznacza dowolne obciążenie poprzeczne przypadające na jednostkę długości elementu, natomiast

\[ [\mathbf N]= \begin{bmatrix} N_1(\xi) & N_2(\xi) & N_3(\xi) & N_4(\xi) \end{bmatrix} \tag{II.43} \label{II.43} \]

jest macierzą ścisłych funkcji kształtu odpowiadających klasycznym stopniom swobody elementu belkowego. Po przejściu do współrzędnej bezwymiarowej \(\xi=x/L\), przy czym \(dx=L\,d\xi\), zależność (\ref{II.42}) przyjmuje postać

\[ \{P\}^{[e]} = L \int_0^1 \{\mathbf N\}^{T}(\xi)\, p(\xi)\, d\xi. \tag{II.44} \label{II.44} \]

W szczególnym przypadku równomiernie rozłożonego obciążenia \(p(x)=q=\mathrm{const}\) otrzymuje się ścisły równoważny wektor obciążeń węzłowych

\[ \{P_q\}^{[e]}_{NL} = \frac{qL}{2} \left\{ 1\,;\, L\left( \frac{1}{\mu^{2}} -\frac{\cot(\mu/2)}{2\mu} \right)\,;\, 1\,;\, -L\left( \frac{1}{\mu^{2}} -\frac{\cot(\mu/2)}{2\mu} \right) \right\}^{T}. \tag{II.45}\label{II.45}
\]

W granicy \(\mu\rightarrow0\) ścisła zależność przechodzi dokładnie w klasyczny równoważny wektor obciążeń elementu Bernoulliego 

\[ \{P_q\}^{[e]}_{B} = \lim_{\mu\rightarrow0}\{P_q\}^{[e]}_{NL} = \frac{qL}{2} \left\{ 1\,;\, \frac{L}{6}\,;\, 1\,;\, -\frac{L}{6} \right\}^{T}. \tag{II.46}\label{II.46} \]

Otrzymany wynik potwierdza, że ścisły równoważny wektor obciążeń stanowi naturalne uogólnienie klasycznego wektora obciążeń elementu Bernoulliego na przypadek belki-słupa obciążonej siłą osiową. Różnica pomiędzy ścisłym i klasycznym wektorem obciążeń wynosi

\[ \Delta\{P_q\}^{[e]} = \{P_q\}^{[e]}_{NL} – \{P_q\}^{[e]}_{B} = \frac{qL^{2}}{2} \left( \frac{1}{\mu^{2}} – \frac{\cot(\mu/2)}{2\mu} – \frac{1}{6} \right) \left\{ 0\,;\, 1\,;\, 0\,;\, -1 \right\}^{T}.\tag{II.47}\label{II.47} \]

Ponieważ

\[\lim_{\mu\rightarrow0} \left( \frac{1}{\mu^{2}} – \frac{\cot(\mu/2)}{2\mu} \right) = \frac{1}{6}, \]

to zależność (\ref{II.47}) spełnia warunek zgodności z klasyczną teorią Bernoulliego, tj.

\[ \lim_{\mu\rightarrow0} \Delta\{P_q\}^{[e]} = \mathbf{0}. \tag{II.48}\label{II.48} \]

Oznacza to, że ścisły równoważny wektor obciążeń stanowi spójne uogólnienie klasycznego wektora obciążeń elementu belkowego i w sposób ciągły przechodzi do rozwiązania liniowego wraz z zanikiem wpływu siły
osiowej.
Zależność (\ref{II.46}) pokazuje, że wpływ siły osiowej nie zmienia składowych poprzecznych równoważnego wektora obciążeń, lecz wyłącznie modyfikuje momenty węzłowe. Wielkość tej korekty jest opisana funkcją

\[ \frac{\cot(\mu/2)}{2\mu} -\frac{1}{\mu^{2}} -\frac{1}{6}, \tag{II.49}\label{II.49} \]

która zanika dla \(\mu\rightarrow0\), dzięki czemu ścisły wektor obciążeń przechodzi w sposób ciągły w klasyczny wektor obciążeń elementu Bernoulliego.

Konsekwencje zastosowania ścisłego wektora obciążeń węzłowych

Wyprowadzenie ścisłego równoważnego wektora obciążeń węzłowych zakończyło budowę ścisłego elementu belki-słupa. W przeciwieństwie do klasycznego elementu Bernoulliego-Eulera zmianie ulega bowiem nie tylko macierz sztywności wynikająca z zastosowania ścisłych funkcji kształtu, lecz również równoważny wektor obciążeń węzłowych. Oba te elementy stanowią nierozłączną część sformułowania metody elementów skończonych i powinny być wyprowadzone z tego samego modelu mechanicznego. Zastosowanie ścisłych funkcji kształtu wyłącznie do budowy macierzy sztywności, przy jednoczesnym pozostawieniu klasycznego wektora obciążeń Bernoulliego, prowadzi do wewnętrznej niespójności modelu elementu. 
Konsekwencje tej niespójności zostały wykazane przez autora już w roku 2014. Wykazano wówczas, że zastosowanie ścisłej macierzy podatności Livesleya \( [K]_{L}^{-1} \) wraz z klasycznym równoważnym wektorem obciążeń Bernoulliego \( \mathbf{F}_B \) prowadzi do rozwiązania różniącego się od ścisłego rozwiązania równania belki-słupa dokładnie o klasyczne rozwiązanie Bernoulliego. Wynik ten można przedstawić w postaci wunikająceej z (\ref{II.3})

\[[K]_{L}^{-1}\mathbf{F}_B+\mathbf{q}_B = \mathbf{q}_{L}, \tag{II.50}\label{II.50} \]

gdzie \( \mathbf{q}_B \) oznacza klasyczny wektor przemieszczeń elementu Bernoulliego, natomiast \( \mathbf{q}_{L} \) jest ścisłym rozwiązaniem elementu belki-słupa.

Po przekształceniu powyższej zależności otrzymuje się

\[ [K]_{L}^{-1}\Delta\mathbf{F} = \mathbf{q}_B, \tag{II.51}\label{II.51} \]

gdzie

\[ \Delta\mathbf{F} = \mathbf{F}_{L}-\mathbf{F}_B. \tag{II.52}\label{II.52} \]

Równanie (\ref{II.51}) ma fundamentalne znaczenie interpretacyjne. Pokazuje ono jednoznacznie, że różnica pomiędzy rozwiązaniem klasycznym i ścisłym nie wynika z zastosowania nieliniowej macierzy sztywności. Przyczyną tej różnicy jest wykorzystanie klasycznego równoważnego wektora obciążeń, który nie jest zgodny z nieliniowymi funkcjami kształtu wykorzystanymi do budowy macierzy sztywności. Oznacza to, że klasyczny element Livesleya nie stanowi jeszcze rozwiązania całkowicie ścisłego. Ścisła jest jedynie macierz sztywności, natomiast równoważny wektor obciążeń pozostaje nadal oparty na klasycznych funkcjach kształtu Bernoulliego. Dopiero jednoczesne zastosowanie ścisłej macierzy sztywności oraz ścisłego równoważnego wektora obciążeń prowadzi do pełnego, wzajemnie zgodnego sformułowania elementu belki-słupa.

Dla obciążenia równomiernie rozłożonego różnica równoważnych wektorów obciążeń opisana jest zależnością (\ref{II.47}). Wynika z niej, że wpływ siły osiowej nie zmienia równoważnych sił poprzecznych, lecz wyłącznie równoważne momenty węzłowe. Zmiana ta powoduje odpowiednią korektę przemieszczeń węzłowych, pola przemieszczeń oraz rozkładu sił wewnętrznych, prowadząc do pełnej zgodności rozwiązania metody elementów skończonych ze ścisłym rozwiązaniem równania różniczkowego belki-słupa.

Przedstawione zależności prowadzą do istotnego wniosku. Wprowadzenie ścisłego równoważnego wektora obciążeń powoduje, że rozwiązanie szczególne równania różniczkowego zostaje automatycznie uwzględnione w równaniu kanonicznym metody elementów skończonych. W konsekwencji całkowite rozwiązanie elementu otrzymuje się bez konieczności dodawania klasycznego rozwiązania Bernoulliego \( \mathbf{q}_B \), ponieważ jego wpływ został już zawarty w ścisłym równoważnym wektorze obciążeń.

Ugięcie pręta nieliniowego

Ścisłe funkcje kształtu (\ref{II.35}) umożliwiają wyznaczenie pola przemieszczeń w dowolnym punkcie elementu wyłącznie na podstawie współrzędnych węzłowych. Po wprowadzeniu wielkości symetrycznych i różnicowych

\[ \Sigma w=w_1+w_2,\qquad \Delta w=w_2-w_1, \]

\[\Sigma\Phi=L(\varphi_1+\varphi_2),\qquad \Delta\Phi=L(\varphi_2-\varphi_1), \]

oraz korzystając z definicji (\ref{II.41}) funkcję ugięcia elementu nieliniowego można zapisać w zwartej postaci

\[ \begin{aligned} w(\xi)=& \cfrac{\Sigma w}{2} -\cfrac{\Delta w}{2\Theta} \left[ 2\cos(\mu\xi) -2\cos\!\bigl(\mu(1-\xi)\bigr) +\mu(2\xi-1)\sin\mu \right] \\[2mm] &-\cfrac{\Delta\Phi}{\mu} \csc\!\left(\cfrac{\mu}{2}\right) \sin\!\left(\cfrac{\mu\xi}{2}\right) \sin\!\left(\cfrac{\mu(1-\xi)}{2}\right) \\[2mm] & -\cfrac{\Sigma\Phi}{2\Theta} \left[ -2\xi +(2\xi-1)\cos\mu -\cos(\mu\xi)
+\cos\!\bigl(\mu(1-\xi)\bigr) +1 \right]. \end{aligned} \tag{II.53} \label{II.53} \]

Zależność (\ref{II.53}) stanowi ścisłe rozwiązanie jednorodnego równania belki-słupa wyrażone za pomocą klasycznych współrzędnych uogólnionych elementu. W przeciwieństwie do klasycznej interpolacji Hermite’a uwzględnia ona w sposób ścisły wpływ siły osiowej poprzez parametr \(\mu\), dzięki czemu zachowuje ważność również w zakresie geometrycznej nieliniowości.

W granicy \(\mu\rightarrow0\) zależność (\ref{II.53}) przechodzi dokładnie w klasyczną interpolację Hermite’a stosowaną w liniowej teorii belek Bernoulliego, co potwierdza zgodność ścisłego rozwiązania z klasycznym elementem belkowym.

Po rozwiązaniu układu równań elementu otrzymuje się współrzędne uogólnione opisujące część jednorodną rozwiązania. Podstawiając je do zależności interpolacyjnej (\ref{II.61}) oraz wyznaczając ugięcie w środku rozpiętości otrzymuje się bezwymiarową funkcję ugięcia $ Y=\frac{f}{f_h}$, gdzie \(f_h =5pL^4/(384EI)\) oznacza ugięcie belki liniowej swobodnie podpartej poddanej równomiernie rozłożonemu obciążeniu \(p\).

Zgodnie z interpretacją przedstawioną w punkcie  Interpretacja rozwiązania elementu belki-słupa  pole przemieszczeń elementu stanowi sumę rozwiązania jednorodnego oraz rozwiązania szczególnego zgodnie z zależnością (\ref{II.61}). W analizowanym przykładzie rozwiązanie szczególne nie jest wyznaczane jawnie. Jego wpływ został uwzględniony poprzez równoważny wektor obciążeń węzłowych \(\mathbf{P}_q\), wyznaczony zgodnie z zasadą prac wirtualnych. Rozwiązanie układu równań elementu prowadzi zatem bezpośrednio do całkowitego ugięcia belki, przy zachowaniu ścisłej zgodności z rozwiązaniem równania różniczkowego.

Po podstawieniu zależności Livesleya oraz wykorzystaniu definicji funkcji stateczności otrzymuje się następującą zależność na bezwymiarowe ugięcie pręta

$Y=\frac{f}{f_h}=\ldots$

Otrzymana zależność stanowi ścisłe rozwiązanie analizowanego zagadnienia i będzie stanowiła podstawę dalszej analizy oraz porównania z rozwiązaniem uzyskanym metodą Ritza.

Strzałka ugięcia elementu nieliniowego

Wzór (\ref{II.61}) umożliwia wyznaczenie ugięcia elementu w dowolnym punkcie jego długości. Szczególne znaczenie praktyczne ma maksymalne ugięcie pręta, zwane dalej strzałką ugięcia. W ogólnym przypadku współrzędną punktu ekstremalnego wyznacza się z warunku zaniku pierwszej pochodnej funkcji ugięcia względem współrzędnej bezwymiarowej $ \cfrac{dw(\xi)}{d\xi}=0.$

W praktycznych zastosowaniach inżynierskich współrzędną punktu maksymalnego ugięcia najczęściej wyznacza się numerycznie lub bezpośrednio z wykresu linii ugięcia. W dalszej części pracy rozpatrywane są głównie przypadki symetrycznych warunków brzegowych oraz symetrycznego obciążenia, dla których punkt maksymalnego ugięcia znajduje się w środku elementu, tj. $ \xi=\cfrac12.$ Po podstawieniu tej współrzędnej do zależności (\ref{II.61}) otrzymuje się prosty wzór na strzałkę ugięcia

\[ f= w\!\left(\cfrac12\right) = \cfrac12 \left[ \Sigma w – \cfrac{\tan(\mu/4)}{\mu}\, \Delta\Phi \right], \tag{II.54} \label{II.54} \] gdzie

\[ \Sigma w=w_1+w_2, \qquad \Delta\Phi=L(\varphi_2-\varphi_1). \]

Powyższą zależność otrzymuje się bezpośrednio przez podstawienie współrzędnej \(\xi=\cfrac12\) do funkcji ugięcia (\ref{II.61}).

Dla granicznego stanu bifurkacji Eulera

\[ \bar{\Lambda}=1, \qquad \mu=\pi, \]

zależność (\ref{II.54}) upraszcza się do postaci

\[ f_{cr} = \cfrac12 \left( \Sigma w – \cfrac{\Delta\Phi}{\pi} \right), \tag{II.55} \label{II.55} \]

opisującej strzałkę ugięcia elementu w punkcie bifurkacji. Zależność ta stanowi szczególny przypadek ogólnego rozwiązania (\ref{II.54}) i odpowiada osiągnięciu przez element obciążenia krytycznego Eulera.

Równanie (\ref{II.55}) określa geometrię postaci wyboczeniowej, natomiast jej amplituda zależy od przemieszczeń i obrotów węzłów wyznaczanych z globalnego układu równań konstrukcji. W pobliżu obciążenia krytycznego mogą one osiągać bardzo duże wartości wskutek utraty odwracalności macierzy sztywności. W konsekwencji amplituda postaci wyboczeniowej zależy od rozwiązania całego układu konstrukcyjnego, imperfekcji geometrycznych oraz sposobu obciążenia.

Dla zerowej siły osiowej $ N=0, \qquad \mu=0$  i zależność (\ref{II.54}) przechodzi w granicy do postaci

\[ f_0= \lim_{\mu\rightarrow0}f = \cfrac12 \left( \Sigma w – \cfrac{\Delta\Phi}{4} \right) = \cfrac12 \left( \Sigma w – \cfrac{L\,\Delta\varphi}{4} \right), \tag{II.56} \label{II.56} \]

gdzie $ \Delta\varphi=\varphi_2-\varphi_1. $. Otrzymana zależność jest zgodna z klasycznym wzorem na strzałkę ugięcia elementu Hermite’a, co stanowi dodatkową weryfikację poprawności wyprowadzonego rozwiązania. W przykładach II.P1–II.P4 przedstawiono zastosowanie otrzymanych zależności do wybranych zagadnień belek-słupów o różnych warunkach brzegowych i sposobach obciążenia.

Na rys.  II.3 zaprezentowano wyniki uzyskane dla przykładu II.P1 dotyczącego belki-słupa z  prawą podporą sprężystą. 

Na rys.  II.3 zaprezentowano wyniki uzyskane dla przykładu II.P2 dotyczącego słupa utwierdzonego sprężyście i obciążonego poprzecznie obciążeniem rozłożonym $h$. 

 Z rys. II.2 i II.3 wynika, że odpowiedź układu jest wynikiem silnego współdziałania dwóch zjawisk: geometrycznej nieliniowości wywołanej siłą osiową oraz podatności podpory sprężystej. W przeciwieństwie do klasycznej teorii wyboczenia, w której stan krytyczny utożsamiany jest z osiągnięciem siły Eulera, w analizowanym układzie podatność podpory może prowadzić do gwałtownego wzrostu ugięcia przy znacznie mniejszych wartościach obciążenia. Oznacza to, że zachowanie konstrukcji jest determinowane nie tylko przez stateczność pręta, lecz również przez rzeczywistą sztywność warunków podparcia.
Podatność podpory i efekt wyboczeniowy nie są zjawiskami niezależnymi, lecz pozostają ze sobą w silnym sprzężeniu. Geometryczna nieliniowość powoduje wzrost wrażliwości układu na podatność podpory, natomiast skończona sztywność podpory przyspiesza rozwój efektów drugiego rzędu.

W rezultacie odpowiedź konstrukcji nie może być interpretowana wyłącznie w kategoriach klasycznej siły krytycznej Eulera, lecz wymaga jednoczesnego uwzględnienia obu mechanizmów. .
Wraz ze wzrostem bezwymiarowej sztywności podpory wpływ jej podatności stopniowo maleje, a przebiegi szybko zbliżają się do rozwiązania odpowiadającego podporze sztywnej. Oznacza to, że dla odpowiednio dużych wartości $\bar C_\Delta$ dalsze zwiększanie sztywności podpory wywiera już niewielki wpływ na ugięcie belki, natomiast o zachowaniu układu decydują przede wszystkim efekty geometrycznej nieliniowości związane ze zbliżaniem się do klasycznego stanu krytycznego Eulera. Uzyskane wyniki mają również istotne znaczenie praktyczne. Pokazują one, że projektowanie elementów ściskanych z podatnymi podporami (czyli w praktyce  każdej  konstrukcji ) nie powinno opierać się wyłącznie na klasycznej analizie wyboczeniowej. W wielu przypadkach o zachowaniu konstrukcji decyduje bowiem wzajemne oddziaływanie podatności podpór i efektów drugiego rzędu, które może prowadzić do znacznego zwiększenia przemieszczeń jeszcze przed osiągnięciem klasycznej siły krytycznej Eulera.
W rozważanym przypadku  asymptota Eulera  została zachowana tylko dlatego, że swoboda obrotu na podporze (1) lub (2) spowodowała, że w stanie krytycznym różnica kątów obrotu rośnie nieskończenie. 

Wnioski

Wyrażenia (\ref{II.61}) i (\ref{II.55}) wskazują, że w granicznym stanie krytycznym Eulera:

(1) W punkcie bifurkacyjnym Eulera funkcja ugięcia jest kombinacją funkcji liniowych oraz trygonometrycznych. Oznacza to, że postać wyboczeniowa należy do przestrzeni funkcji generowanej przez składniki stałe, liniowe oraz harmoniczne. Własność ta ma istotne znaczenie przy budowie funkcji aproksymacyjnych w metodzie Ritza i metodzie elementów skończonych, ponieważ uzasadnia stosowanie baz wielomianowych, trygonometrycznych oraz ich kombinacji do opisu zachowania konstrukcji w pobliżu punktu bifurkacji.

(2) Strzałka ugięcia pręta pozostaje wielkością skończoną i niezerową o ile tylko skończone są węzłowe przemieszczenia i kąty obrotu.
Równanie (\ref{II.55}) prowadzi do istotnego wniosku dotyczącego interpretacji klasycznej teorii stateczności. W punkcie bifurkacyjnym Eulera [3], odpowiadającym wartości (\bar{\Lambda}=1), strzałka ugięcia pozostaje wielkością skończoną i jest jednoznacznie określona przez przemieszczenia oraz obroty węzłowe. Oznacza to, że samo osiągnięcie obciążenia krytycznego nie prowadzi do matematycznej osobliwości funkcji ugięcia ani do nieskończonych przemieszczeń konstrukcji. Obciążenie Eulera wyznacza punkt bifurkacji, w którym układ traci jednoznaczność rozwiązania i pojawia się możliwość przejścia na nową gałąź równowagi, nie stanowi jednak dowodu utraty zdolności przenoszenia obciążenia przez pręt. Konstrukcja może nadal pozostawać w stanie równowagi i przenosić obciążenia osiowe przy zmienionej konfiguracji geometrycznej.|

(3) Wynik pozostaje w pełnej zgodności z teorią nośności nadkrytycznej konstrukcji oraz z licznymi obserwacjami eksperymentalnymi. Badania prowadzone od początku XX wieku wykazały, że wiele układów prętowych, płytowych i powłokowych zachowuje zdolność przenoszenia obciążeń również po osiągnięciu obciążenia krytycznego Eulera. Szczególne znaczenie miały prace Heinza Wagnera dotyczące pracy nadkrytycznej cienkościennych elementów konstrukcyjnych oraz późniejsze uogólnienia teorii bifurkacji i stateczności opracowane przez Koitera. W teoriach tych obciążenie krytyczne nie jest utożsamiane ze zniszczeniem konstrukcji, lecz z utratą stabilności dotychczasowej gałęzi równowagi i pojawieniem się nowych możliwych stanów równowagi. Heinz Wagner — już w 1931 r. [4] wykazał, że po wyboczeniu cienkie środniki i płyty mogą nadal przenosić obciążenie dzięki pracy w stanie nadkrytycznym (teoria pola rozciąganego).  Theodore von Kármán i Hsue-Shen Tsien [5] pokazali, że po osiągnięciu obciążenia krytycznego możliwe są dalsze gałęzie równowagi, a zachowanie konstrukcji zależy od charakteru bifurkacji. Warner T. Koiter [6]  w pracy doktorskiej z 1945 r. stworzył ogólną teorię stateczności i zachowania nadkrytycznego konstrukcji sprężystych.

W punkcie bifurkacyjnym Eulera, odpowiadającym wartości $\bar{\Lambda}=1$, strzałka ugięcia pozostaje wielkością skończoną i jest jednoznacznie określona przez skończone przemieszczenia oraz obroty węzłowe. Przemieszczenie i obroty węzłowe pozostają skończone w przypadku, gdy zapewniają to kinematyczne warunki brzegowe, np dla belki utwierdzono- utwierdzonej.  W przykładach z rys II.2 i II.3  wyodrębnione elementy konstrukcji mają swobodne węzły, co skutkuje ujawnieniem się asymptot krytycznych., cow konstrukcjach rzeczywistych będzie osłabione  , a nawet wyeliminowane poprzez przeskok na pokrytyczną  gałąź ścieżki równowagi. Wynik ten pozostaje w pełnej zgodności z teorią nośności nadkrytycznej konstrukcji oraz z licznymi obserwacjami eksperymentalnymi. Badania prowadzone od początku XX wieku wykazały, że wiele układów prętowych, płytowych i powłokowych zachowuje zdolność przenoszenia obciążeń również po osiągnięciu obciążenia krytycznego Eulera. Szczególne znaczenie miały prace Heinza Wagnera dotyczące pracy nadkrytycznej cienkościennych elementów konstrukcyjnych oraz późniejsze uogólnienia teorii bifurkacji i stateczności opracowane przez Koitera. W teoriach tych obciążenie krytyczne nie jest utożsamiane ze zniszczeniem konstrukcji, lecz z utratą stabilności dotychczasowej gałęzi równowagi i pojawieniem się nowych możliwych stanów równowagi.

Imperfekcje przechyłowe i ich równoważniki 

Rzeczywiste imperfekcje geometryczne przechyłowe  w praktyce perojektowej zastępuje się  imperfekcjami zastępczymi (równoważnymi).  

Ścisła macierz sztywności oraz odpowiadająca jej macierz podatności, wyprowadzone w poprzednich rozdziałach, umożliwiają ilościową ocenę tej równoważności bez stosowania dodatkowych założeń upraszczających. Dzięki temu możliwe jest określenie zakresu stosowalności klasycznych modeli zastępczych oraz ocena wpływu kolejnych przybliżeń teorii na dokładność odwzorowania rzeczywistych imperfekcji przechyłowych.

Globalne imperfekcje przechyłowe należą do podstawowych źródeł geometrycznej nieliniowości konstrukcji prętowych. Powstają wskutek niedokładności wykonania i montażu, odchyłek geometrycznych, przemieszczeń podpór oraz innych czynników powodujących odchylenie osi konstrukcji od położenia projektowego. W obecności siły osiowej imperfekcje te wywołują dodatkowe momenty drugiego rzędu, których wartość rośnie wraz ze wzrostem obciążenia ściskającego. W praktyce inżynierskiej rzeczywiste imperfekcje przechyłowe zastępowane są najczęściej równoważnymi obciążeniami poziomymi. Dla przechyłu konstrukcji

\[ \Delta \stackrel{\mathrm{def}}{=} \phi L, \tag{II.57}\label{II.57} \]

gdzie $\phi$ oznacza kąt przechyłu, który będziemy oznaczali jako $1/n_L$. Podstawowa wartość podzielnika (faktora) impefekcij wynosi  $N_{L0}= 200$, ale w praktyce może dochodzić do $n_L=500$

W praktyce analizy konstrukcji imperfekcję przechyłową zastępuje się   przez równoważną  siłę poziomą stowarzyszoną z kazym obciązęniem grawiacyjnym $F $ (skupionym, rozłożonym lub innwj natury prowadzącym do „ściskania”  konstrukcji.     

\[ H_{eq} =  F \phi = \cfrac{F}{n_L} \tag{II.58}\label{II.58} \]

Imperfekcja zastępcza (\ref{II.58}) wynika  z warunku równoważności momentów $ H \cdot L = F \cdot \Delta.$.  Przyjmuje się , że tak zdefiniowana siła wywołuje ten sam globalny moment drugiego rzędu $N\Delta$ , a metoda równoważnych obciążeń została powszechnie przyjęta w projektowaniu konstrukcji stalowych, żelbetowych oraz innych układów prętowych. Przyjęta równoważność ma jednak charakter wyłącznie statyczny. Rzeczywista imperfekcja przechyłowa stanowi wymuszenie geometryczne, natomiast jej równoważnik obciążeniowy jest zewnętrznym oddziaływaniem mechanicznym. Oba modele mają zatem odmienny charakter fizyczny i mogą prowadzić do różnych odpowiedzi konstrukcji, zwłaszcza w analizie geometrycznie nieliniowej.

W przypadku konstrukcji żelbetowych często stosuje się zastępczy (równoważny) mimośród  działania obciążeń  $e_{eq}$, w tym mimośró niezamierzzony integrujący imperfekcje elementu. 

Kryterium przemieszczeniowe

Kryterium przemieszczeniowe  oznacza porównanie przemieszczeń (liniowych , obrotów, itp) w określonym punkcie konstrukcji i można je zapisać w prostej postaci dla punktu przyłożenia siły $F$:

\[ \delta (F_{eq}) =\Delta \tag{II.59}\label{II.59} \]

gdzie:
$\Delta$ – imperfekcją przechyłowa (\ref{II.57}) ,
$\delta (F_{eq})$ – przmieszczenie wywołane zastępczym wymuszeniem, a najcześciej zastęczą siła poziomą $H_{eq}$ lub zastępczym miośrodem $e_{eq}$ 

Kryteria równoważności statycznej i energetycznej

w poprzednim punkcie jako kryterium równoważności przyjęto zgodność odpowiedzi konstrukcji, rozumianą jako zgodność przemieszczeń uogólnionych. Obejmuje ona przemieszczenia liniowe, obroty oraz ugięcia charakterystycznych punktów konstrukcji. Kryterium to posiada bezpośrednią interpretację fizyczną i umożliwia jednoznaczną ocenę, czy rzeczywista imperfekcja przechyłowa oraz odpowiadające jej obciążenie zastępcze wywołują taki sam stan deformacji konstrukcji.

Kolejnymi kryteriami oceny są:

• zgodność sił przekrojowych lubodpowiadających im naprężeń (kryterium statyczne), 
• zgodność energii odkształcenia zgromadzonej w konstrukcji (kryterium energetyczne).

Kryterium statyczne przyjmuje się najczęściej jako warunek równości momentów zginajacych wywołanych  imperfekcją zastępczą i  wymuszeniem imperefekcjąprzmiszczeniopwą   $\Delta$  \ref{II.49$)

Równoważnośc imperfekcji w stanie krytycznym konstrukcji

Analiza stateczności konstrukcji posiadającej rzeczywistą imperfekcję geometryczną oraz analiza konstrukcji idealnej z równoważnym obciążeniem poziomym prowadzą na ogół do różnych wartości obciążenia krytycznego

 \ [ \Lambda_{cr,\Delta}\neq\Lambda_{cr,H}. \]

Oznacza to, że klasyczne równoważniki obciążeniowe nie zachowują równoważności statecznościowej. Już sam ten fakt wskazuje, że równoważność oparta wyłącznie na warunku \(HL=N\Delta\) ma charakter ograniczony i wymaga dalszej oceny w zakresie odpowiedzi przedkrytycznej konstrukcji.  Jakość modeli zastępczych istotnie zalreży od poziomu obciążęnia konstrukcji. oraz sztywności  konstrukcji podpierajacych. 

Zastosowanie kryterium przmieszczeniowego

Równoważność siły poziomej $H_{eq}$

W literaturze trudno znaleźć ścisłe wyprowadzenie zależności  $H=N/nLH=N/n_LH=N/nL$​ z ogólnych zasad mechaniki. Jest ona uzasadniane przede wszystkim prosyotą w zastoswaniach praktycznych i prowadzi do poprawnych wyników inżynierskich, jednak nie stanowi konsekwencji zasady minimum energii potencjalnej ani innego globalnego kryterium mechaniki konstrukcji.

Poniżej sprawdzimy dobroć dopasowania imperfekcji zastępczych do wynikających z definicji (\ref{II.57} \]  a także  proponuje się alternatywne kryterium równoważności oparte na energii odkształcenia. W przeciwieństwie do klasycznego kryterium momentowego ma ono charakter globalny i wynika bezpośrednio z zasad wariacyjnych mechaniki. Z punktu widzenia praktyki projektowej otrzymane poniżej  wyniki prowadzą do istotnego wniosku, że  wartość imperfekcji przechyłowej jest zależna nie tylko od rodzaju konstrukcji, ale także od poziomu jej obciążenia, sztywności układu oraz konstuklckcji podporowoej , co nie zapewnia stałej jakości odwzorowania rzeczywistego zachowania konstrukcji.

Wraz ze wzrostem względnego poziomu obciążenia \(\bar{\Lambda}\) szybko narasta wpływ nieliniowości geometrycznej, wskutek czego zgodność pomiędzy rozwiązaniem liniowym i ścisłym ulega systematycznemu pogorszeniu.
Stała wartość imperfekcji przechylwej nie może być traktowana  jako uniwersalny równoważnik mechaniczny rzeczywistej imperfekcji przechyłowej. Jest ona jedynie praktycznym parametrem projektowym, którego dokładność odwzorowania zależy od stanu pracy konstrukcji. W konsekwencji ta sama wartość imperfekcji może bardzo dobrze opisywać zachowanie konstrukcji przy niewielkich obciążeniach, natomiast w zakresie obciążeń zbliżonych do stanu krytycznego prowadzić do coraz większych odchyleń od rozwiązania ścisłego.

Na rys. II.4 przedstawiono wyniki analizy geometrycznie nieliniowego wspornika obciążonego poziomą siłą skupioną \(H\), przyłożoną do głowicy słupa. W przeciwieństwie do przykładu II.2, w którym analizowano odpowiedź konstrukcji na dowolną siłę poziomą, przedstawiona analiza służy do wyznaczenia zastępczej siły poziomej równoważnej zadanej imperfekcji geometrycznej słupa. Jako kryterium równoważności (\ref{II.59}) przyjęto zgodność poziomego przemieszczenia $\delta(H_{eq}$ głowicy wspornika wywołanego zastępczą siłą poziomą \(H_{eq}= \) \cfrac{N}{n_L}.$ z normową imperfekcją przechyłową $ \Delta=\cfrac{L}{n_L}.$

Formuła analityczna na funkcję przedstawiona  na wykresie rys. II.4 jest następująca:

\[ \bar{f}_{H}= 1 – \cfrac{12}{\pi^3\bar{\Lambda}^{3/2}} \left( \pi\sqrt{\bar{\Lambda}} + \cfrac{2\bar{C}_{\varphi} \cdot s} {\alpha \cdot s -\bar{C}_{\varphi} \cot c} \right), \tag{II.60}\label{II.60} \]

gdzie:
$\alpha=\cfrac{\pi}{2}\sqrt{\bar{\Lambda}}$
$c=\cos\alpha$$,
$s= \sin\alpha.$

Uwzględniając, że: liniowe przemieszczenie głowicy wspornika od siły wynosi $H_{eq}$  wynosi $ \delta_0=\cfrac{H_{eq}L^3}{3EI},$;  rzeczywiste przemieszczenie opisuje zależność  $ \delta=\bar f\,\delta_0,$ ; siła krytyczna wspornika ma postać $ N_{cr}=\cfrac{\pi^2EI}{4L^2}$, to po podstawieniu powyższych zależności otrzymuje się zwartą postać współczynnika zgodności

\[ \eta_{H=\Delta} = \cfrac{\delta(H_{eq})}{\Delta} = \cfrac{\pi^2}{12}\, \bar f\, \bar{\Lambda} \approx 0.823\,\bar f\,\bar{\Lambda}. \tag{II.61}\label{II.61}\]

Współczynnik \(\bar f\) odczytuje się z wykresu rys. II.4 dla zadanego poziomu wytężenia \(\bar{\Lambda}\) oraz bezwymiarowej sztywności podpory \(\bar C_{\varphi}\). Współczynnik (\ref{II.61}) stanowi bezwymiarową miarę zgodności pomiędzy przemieszczeniem wywołanym zastępczą siłą poziomą \(H_{eq}=N/n_L\) a odpowiadającą jej normową imperfekcją przechyłową \(\Delta=L/n_L\). Charakterystyczną cechą otrzymanej zależności jest całkowite wyeliminowanie współczynnika \(n_L\), co oznacza, że jakość równoważności obu modeli nie zależy od przyjętej wartości imperfekcji normowej, lecz wyłącznie od poziomu obciążenia \(\bar{\Lambda}\) oraz sztywności podpory \(\bar C_{\varphi}\).

Warunek pełnej równoważności obu modeli spełniony jest dla $ \eta_{H=\Delta}=1$. Dla \(\eta_{H=\Delta}<1\) zastępcza siła pozioma wywołuje przemieszczenie mniejsze od odpowiadającej imperfekcji przechyłowej, natomiast dla \(\eta_{H=\Delta}>1\) powoduje przemieszczenie większe od odpowiadającej jej imperfekcji.

W tabeli poniżej przedstawiono jakość dopasowania dla wspornika z podporą o sztywności \(\bar C_{\varphi}=100\). 

\[ \begin{array}{c|c|c}
\bar{\Lambda} & \bar{f} & \eta_{H=\Delta} \\
\hline 0.20 & 2.2905 & 0.3768\\
0.4304 & 2.8336 & 1.0000\\
0.50 & 3.0876 & 1.2696\\
0.90 & 13.4096 & 9.9259 \end{array} \]

Pełna zgodność obu modeli występuje jedynie dla \(\bar{\Lambda}\approx 0,430\). Dla mniejszych wartości \(\bar{\Lambda}\) model zastępczej siły poziomej prowadzi do niedoszacowania efektu imperfekcji przechyłowej, natomiast dla większych poziomów obciążenia powoduje jego coraz większe przeszacowanie.

Równoważność  mimośrodu $e_{eq}$

Na rys. II.5 przedstawiono wykres ugięcia względnego $\bar f={\delta(e_q)}{\delta_0} dla wspornika obciążonego w wierzchołku siłą $F=N$ na mimośrodzie $e_{eq}$, czyli równoważnie momentem zginającym $M= N\cdot e_{eq}$ 

gdzie wartość odniesienia $ \delta_0= \cfrac{M \cdot L^2}{2 EI}= \cfrac{N \cdot e_{eq} \cdot L^2}{2EI}$ jest klasycznym ugięciem liniowego wspornika.

Formuła analityczna przedstawiona na wykresie to:

\[ \bar{f}_{e} = \cfrac{-16\pi^{2}\bar{\Lambda}\left[\pi\sqrt{\bar{\Lambda}}\,s\sqrt{1-s^{2}}-\bar{C}_{\varphi}\left(1-2s^{2}\right)\right]-64\,\bar{C}_{\varphi}\,s\left(\pi\sqrt{\bar{\Lambda}}-8s\right)}{\pi^{4}\bar{\Lambda}^{2}\left[\pi\sqrt{\bar{\Lambda}}\,s\sqrt{1-s^{2}}- \bar{C}_{\varphi}\left(1-2s^{2}\right) \right] }, \tag{II.62}\label{II.62} \]

gdzie: $ s=\sin\!\left(\cfrac{\pi\sqrt{\bar{\Lambda}}}{4}\right), $

Wykres rys.  II.5  sporządzono dla znormalizowanego wygiecia wspornika$\bar f$.  Wartość rzeczywista ugięcia wyniesie natomiast:

\[ \delta=  bar f  \cdot \delta_0 =  \bar f  \cdot \cfrac{N \cdot e_{eq} \cdot L^2}{2EI} \] 

Mimośród równoważny działania siły $N$ w wierzchołku słupa zapiszemy w postaci bezwymiarowej $ \bar e=\cfrac{e_{eq}}{\Delta}, $, gdzie $ \Delta =L/n_L $ – impefekcja przechyłowa (\ref{II.57})

Siłę  osiową we wsporniku przy pozimie obciazenie $\bar \Lambda$ można wyznaczyć z wzoru Eulera  $ N=\bar \Lambda N_{cr} = \bar\Lambda  \cdot \cfrac{\pi^2EI}{4L^2}$

Z powyższych zależności uzyskujemy wyrażenie na współczynnik równoważności  kryterium $\delta_{eq} = \Delta$:

\[ \eta_{[\delta({eq}) =\Delta} = \cfrac{\delta}{\Delta} =\cfrac{\pi^2}{8}  \bar f  \cdot \bar \Lambda  \cdot \bar {e_{eq}}  \approx 1,233 \bar f  \cdot \bar \Lambda  \cdot \bar {e_{eq}} \tag{II.63}\label{II.63} \]

Pierwotna miara projektowej imperfekcji oraz praktyczny sposób jej realizacji

Przeprowadzona analiza zastosowania kryterium przmieszczeniowego wykazała, że zarówno klasyczna imperfekcja przechyłowa $\Delta=\cfrac{L}{n_L},$ jak i odpowiadające jej modele zastępcze w postaci siły poziomej $ H_{eq}=\cfrac{N}{n_L},$
lub mimośrodu $ e_{eq},$ nie stanowią uniwersalnych modeli równoważnych. Otrzymane wyniki wskazują, że stopień zgodności pomiędzy poszczególnymi modelami zależy od poziomu obciążenia, sztywności konstrukcji oraz warunków podparcia. Oznacza to, że geometryczna imperfekcja przechyłowa nie może być traktowana jako jednoznaczny mechaniczny odpowiednik rzeczywistego zachowania konstrukcji. Imperfekcja przechyłowa należy do wielkości najtrudniejszych do jednoznacznego zdefiniowania i zrealizowania w praktyce. Rzeczywista konstrukcja nigdy nie posiada jednej, dokładnie określonej wartości przechyłu początkowego. Na jej zachowanie wpływają jednocześnie odchyłki wykonawcze, krzywizny początkowe elementów, mimośrody montażowe, odkształcenia podpór oraz przemieszczenia powstające podczas realizacji i eksploatacji obiektu. Z tego względu normowa imperfekcja przechyłowa stanowi jedynie uproszczony model projektowy, a nie bezpośrednio mierzalną cechę konstrukcji.

Proponujemy uniwesalną definicję imperferkcji przchyłowych i sposób jej realizacji zwanego dalej metodą statyczną.
Zachowuje się normową interpretację współczynnika

\[ n_L, \tag{II.64}\label{II.64} \]

zwanego dalej  faktorem imperfekcji, który pozostaje podstawową miarą projektowej imperfekcji konstrukcji. Wartość tego współczynnika nadal ustala się zgodnie z obowiązującymi wymaganiami normowymi. Zmianie ulega natomiast sposób realizacji imperfekcji w modelu obliczeniowym.

Zamiast bezpośredniego wprowadzania geometrycznej imperfekcji
\[ \Delta=\cfrac{L}{n_L}, \]  proponuje się jako podstawowy model obliczeniowy stosowanie równoważnego obciążenia zastępczego 

\[ H_d=\cfrac{F_d}{n_L}, \tag{II.65}\label{II.65} \]

gdzie \(F_d\) oznacza projektowe obciążenie grawitacyjne (pionowe) działające na konstrukcję Wielkość \(F_d\) odpowiada rzeczywistemu schematowi obciążenia i może reprezentować siłę skupioną, obciążenie rozłożone, oddziaływanie termiczne lub inne oddziaływanie wyznaczone z odpowiedniej kombinacji obciążeń dla analizowanego stanu granicznego (SGN, SGU, FAT itp.). W takim ujęciu faktor imperfekcji \(n_L\) zachowuje interpretację normową, jednak jego rola ulega uogólnieniu. Staje się on bezwymiarowym współczynnikiem określającym poziom projektowej imperfekcji poprzez odpowiednie skalowanie rzeczywistego obciążenia odniesienia. Obciążeniem tym jest zawsze konfiguracja obciążeń właściwa dla analizowanego stanu granicznego. Przykładowo, dla stanu granicznego nośności przyjmuje się $F_0=F_d$,,  natomiast dla innych stanów granicznych obciążeniem odniesienia jest odpowiednia wartość wynikająca z obowiązujących kombinacji oddziaływań.

W takim ujęciu współczynnik \(n_L\) pozostaje wielkością pierwotną, natomiast geometryczna imperfekcja przechyłowa staje się wielkością wtórną, wyznaczaną z przyjętego kryterium równoważności. Kryterium tym może być zgodność przemieszczeń, energii odkształcenia, pracy sił zewnętrznych lub innego efektu mechanicznego uznanego za reprezentatywny dla analizowanego zagadnienia. Oznacza to odwrócenie klasycznego sposobu interpretacji imperfekcji projektowej. Nie narzuca się z góry wartości geometrycznego przechyłu, lecz wyznacza się jego wartość równoważną odpowiadającą rzeczywistemu oddziaływaniu projektowemu. Proponowane podejście zachowuje pełną zgodność z normową interpretacją współczynnika \(n_L\), a jednocześnie lepiej odzwierciedla rzeczywisty sposób pracy konstrukcji, w której stan naprężenia i przemieszczeń jest wywoływany przez układ działających obciążeń, a nie przez samą imperfekcję geometryczną. Dzięki temu projektowa realizacja imperfekcji staje się bardziej jednoznaczna, łatwiejsza do zastosowania oraz umożliwia wykorzystanie różnych kryteriów równoważności mechanicznej w zależności od analizowanego zagadnienia. W konsekwencji proponuje się zmianę punktu odniesienia całej analizy.

Rozwiązaniem pierwotnym staje się odpowiedź konstrukcji wywołana projektowym obciążeniem zastępczym
\[ H_d=\cfrac{F_d}{n_L}, \] natomiast wszystkie pozostałe modele imperfekcji, w szczególności imperfekcja przechyłowa  $ \Delta=\cfrac{L}{n_L}$ oraz równoważny mimośród $ e_{eq}$, traktowane są jako modele wtórne. Ich zadaniem nie jest już definiowanie imperfekcji, lecz możliwie wierne odtworzenie odpowiedzi konstrukcji uzyskanej dla modelu podstawowego. W dalszej części pracy wszystkie proponowane modele zastępcze będą oceniane poprzez porównanie ich odpowiedzi z rozwiązaniem referencyjnym otrzymanym dla obciążenia \(H_d\)

Pierwotną miarą imperfekcji przychyłowej jest faktor imperfekcji $n_L$. Faktor imperfekcji  jest mnożnikiem konfiguracji odniesienia obciążęń w stanie granicznym $F_d$ i wyznacza równoważne poziome siły imperekcji $H_d= \cfrac{F_d}{n_L} $  stowarzyszone z każdym grawitacyjnycm (pionowym) obciazęniem $F_d$

Równoważność energetyczna

Funkcjonał całkowitej energii potencjalnej można zapisać w postaci

\[  \Pi = U-W, \tag{II.66}\label{II.66} \]

gdzie \(U\) oznacza energię odkształcenia układu, natomiast \(W\) jest pracą sił zewnętrznych. Funkcja przemieszczeń wewnątrz elementu, opisana ścisłymi funkcjami kształtu Livesleya (\ref{II.35}), prowadzi do funkcjonału zawierającego całki postaci $ \int_0^L EI\,N_i”N_j”\,dx,$ oraz  $  \int_0^L N\,N_i’N_j’\,dx.$ Całki te zostały już wyznaczone podczas wyprowadzania dokładnej macierzy sztywności elementu. Ich wyniki zostały ujęte w funkcjach stabilności \(\Psi_1,\Psi_2,\Psi_3,\Psi_4\) zdefiniowanych zależnościami (\ref{II.23})–(\ref{II.24}), na podstawie których zbudowano dokładną macierz sztywności geometrycznej elementu. W konsekwencji funkcjonał całkowitej energii potencjalnej można zapisać w zwartej postaci kwadratowej

\[ \Pi= \cfrac12\, \mathbf{q}^{\mathrm T} \mathbf{K}(\Psi_i) \mathbf{q} – \mathbf{q}^{\mathrm T} \mathbf{P}, \tag{II.67}\label{II.67} \] 

gdzie:  $ \mathbf{q} = \begin{Bmatrix} w_1 & \varphi_1 & w_2 & \varphi_2 \end{Bmatrix}^{\mathrm T},$,
$ \mathbf{K}(\Psi_i)$ oznacza dokładną macierz sztywności elementu zależną od funkcji stabilności.(\ref{II.22}).

Po uwzględnieniu warunków brzegowych oraz eliminacji więzów, analogicznie jak w przykładzie II.P2, liczba niewiadomych ulega dalszemu zmniejszeniu. Po podstawieniu rozwiązania

\[ \mathbf{q} = \mathbf{\tilde K}^{-1}\mathbf{\tilde P}, \tag{II.68}\label{II.68} \] 

gdzie: $ \mathbf{\tilde K}^{-1}$ jest macierzą potanosći konstrukcji , uzyskaną poprzez odwróćenie zmodyfikowanej macierzy sztywności konstrukcji. 
$\mathbf{\tilde P} odpowiadajacy zmodyfikowany wektor róęnoeąznikó węzłowych (p. przyklad II.P1 i II.P2) 

funkcjonał przyjmuje wyjątkowo prostą postać

\[ \Pi = -\cfrac{1}{2}\, \tilde{\mathbf P}^{\,\mathrm T}\, \tilde{\mathbf K}^{-1}\, \tilde{\mathbf P}. \tag{II.69}\label{II.69} \]

gdzie::
$\Pi$ – całkowita energia potencjalna układu,
$\tilde{\mathbf P}$ – zmodyfikowany  wektor obciążeń uogólnionych,
$\tilde{\mathbf K}$ – zmodyfikowana  macierz sztywności układu,
$\tilde{\mathbf K}^{-1}\) – odwrotność zredukowanej macierzy sztywności. (macierz podatnosći)

Wzór (II.58) stanowi ogólną postać energii potencjalnej po rozwiązaniu układu równań równowagi i wyeliminowaniu przemieszczeń uogólnionych. Jest on punktem wyjścia do wyznaczania energii odkształcenia, współczynnika równoważności energetycznej \(\eta_E\) oraz równoważnego mimośrodu \(\bar e\).

Oznacza to, że całkowita energia potencjalna elementu może zostać wyrażona wyłącznie za pomocą elementów odwrotnej macierzy sztywności, bez konieczności ponownego wykonywania całkowania. W szczególności  porównanie modeli z równoważną siłą poziomą oraz z równoważnym mimośrodem można sprowadzić do analizy odpowiednich funkcjonałów energii, co stanowi podstawę sformułowania kryterium równoważności energetycznej.

Rownoważność  zastępczego mimośrodu

Równowazność przemieszczeniowa 

Na rys. II.6. przedstawiono przebieg współczynnika równoważności względnego mimośrodu $bar e=e/\Delta$  w odniesieniu do referencyjnej imperfekcji w postaci zastępczej siły poziomej $(H=N/n_L\$

Współczynnik wyznaczono z kryterium równoważności przemieszczeniowej polegającego na zrównaniu wychylenia wspornika wywołanego działaniem osiowej siły ściskającej $N$ przyłożonej z mimośrodem $e$ z wychyleniem wywołanym działaniem referencyjnej siły poziomej $H$.

Przedstawione przebiegi odpowiadają poziomowi obciążenia $ \bar{\Lambda}=\cfrac{\Lambda}{\Lambda_{cr}}=0{,}5$  odpowiadającemu połowie obciążenia krytycznego Eulera. Wybór ten nie ma jednak charakteru ograniczającego. Analiza numeryczna wykazała bowiem, że współczynnik równoważności $\eta_{(e=H)}$ wykazuje niewielką wrażliwość na poziom obciążenia osiowego w całym analizowanym zakresie wartości $\bar{\Lambda}$. Oznacza to, że zmiana poziomu siły ściskającej wpływa na przebieg krzywych jedynie nieznacznie, dlatego wykres dla $\bar{\Lambda}=0{,}5$można traktować jako reprezentatywny dla pozostałych poziomów obciążenia. wrażliwośc wspóczynnika równowązności na poziom obciązenia.

Współczynnik \(\eta_{(e=H)}\) określa stopień równoważności obu modeli imperfekcji pod względem wywoływanego przemieszczenia oraz umożliwia zastąpienie modelu z zastępczą siłą poziomą równoważnym modelem z mimośrodem. Zależność opisującą tę równoważność zapisano w postaci

\[ \eta_{(e=H)} = \cfrac{H_{\mathrm{eq}}}{H} = \cfrac{1}{1+\cfrac{3}{2}\,\bar e\,\eta_{(H=M)}}. \tag{II.70}\label{II.70} \]

gdzie
\[ \eta_{(H=M)} = \cfrac{\Psi_{2}\left(\bar{C}_{\varphi}+4\Psi_{3}-2\Psi_{4}\right)} {\Psi_{3}\left(\bar{C}_{\varphi}+4\Psi_{3}\right)-\Psi_{4}^{2}} \tag{II.71}\label{II.71} \]

jest współczynnikiem równoważności przemieszczeniowej schematu z momentem zginającym \(M=Ne\) względem referencyjnego schematu z zastępczą siłą poziomą \(H\). Szczegółowe wyprowadzenie zależności (\ref{II.70})–(\ref{II.71}}) przedstawiono w przykładzie II.P3.  .

Przedstawione wyniki wskazują, że współczynnik równoważności maleje monotonicznie wraz ze wzrostem względnego mimośrodu \(\bar{e}\). Tym samym wzrost mimośrodu powoduje systematyczną utratę równoważności względem referencyjnej imperfekcji w postaci zastępczej siły poziomej. W granicy $ \bar{e}\rightarrow\infty$ współczynnik $\eta_{(e=H)}$ dąży do zera, co oznacza, że dla dużych mimośrodów ich odwzorowanie za pomocą równoważnej siły poziomej staje się praktycznie niemożliwe.

Wniosek częściowy. Przeprowadzona analiza prowadzi do wniosku, że równoważność imperfekcji mimośrodowej względem referencyjnej imperfekcji siłowej jest determinowana przede wszystkim przez geometrię imperfekcji, natomiast wpływ poziomu obciążenia osiowego jest drugorzędny. Oznacza to, że geometryczny mimośród nie stanowi uniwersalnej postaci imperfekcji równoważnej zastępczej sile poziomej, nawet przy przyjęciu jedynie kryterium zgodności przemieszczeń. W praktyce projektowej mimośród wymaga zatem odpowiedniej kalibracji względem referencyjnej imperfekcji siłowej i nie powinien być stosowany jako jej bezpośredni zamiennik. W kolejnym punkcie zostanie wykazane, że zastosowanie kryterium statycznego, opartego na porównaniu momentów zginających w podstawie pręta, prowadzi do jeszcze bardziej rygorystycznej oceny równoważności obu modeli.

Równoważność statyczna – momentów podporowych

Na rys. II.7.  wodróżnieniu od poprzedniego podrozdziału, w którym jako kryterium równoważności przyjęto zgodność przemieszczeń, jako podstawę porównania jest zgodność efektów statycznych wyrażonych momentem zginającym w utwierdzeniu. $ M_e=N\,e.$
Jak wykazano w przykładzie II.P3, współczynnik równoważności wynikający z warunku statycznego można zapisać w postaci

\[  \eta_{(e=H)} = \cfrac{M_e}{M_H} = \cfrac{HL/N}{e} = \cfrac{\pi}{2}\, \cfrac{\sqrt{\bar{\Lambda}}} {\sin\!\left(\dfrac{\pi}{2}\sqrt{\bar{\Lambda}}\right)} \tag{II.72}\label{II.72} \]

gdzie:
$e$ – równoważny mimośród obciążenia osiowego,
$H$ – równoważna siła pozioma,
$N$ – siła osiowa ściskająca,
$ (M_e = Ne$ – moment zginający wywołany równoważnym mimośrodem $e$,
$M_H = HL$ – moment zginający wywołany równoważną siłą poziomą $H$,
$\bar{\Lambda}$ – bezwymiarowym parametrem statecznościowym zdefiniowanym zależnością (\ref{II.4}).
Wartość współczynnika $\eta_{(e=H)}$ określa, ile razy moment od modelu mimośrodowego jest większy od momentu odpowiadającego równoważnej sile poziomej.  W granicy $ bar{\Lambda}\rightarrow 0)$ współczynnik dąży do jedności, natomiast wraz ze wzrostem obciążenia  monotonicznie rośnie.

W przeciwieństwie do kryterium przemieszczeniowego współczynnik \(\eta_{[e=H)}\) zależy wyłącznie od poziomu obciążenia osiowego, opisanego parametrem \(\bar{\Lambda}\). Oznacza to, że mimośród równoważny nie jest wielkością stałą, lecz zmienia się wraz ze stanem pracy konstrukcji. Ścisła równoważność pomiędzy obciążeniem mimośrodowym i zastępczą siłą poziomą zachodzi jedynie dla określonej wartości parametru \(\bar{\Lambda}\), a więc jedynie lokalnie w przestrzeni stanów pracy konstrukcji. Na szczególną uwagę zasługuje niemal liniowy charakter otrzymanej zależności. Wraz ze wzrostem obciązenia rośnie współczynnik równoważności 

Wniosek częściowy Kryterium statyczne prowadzi do jeszcze bardziej rygorystycznego wniosku niż kryterium przemieszczeniowe. Równoważność mimośrodu względem referencyjnej imperfekcji w postaci zastępczej siły poziomej zależy bezpośrednio od poziomu obciążenia osiowego. Tym samym wartość mimośrodu równoważnego nie jest wyłącznie cechą geometryczną konstrukcji, lecz zależy również od aktualnego stanu jej pracy.<p>
Wynik ten ma istotne konsekwencje mechaniczne. Imperfekcja geometryczna jest własnością konstrukcji istniejącą przed przyłożeniem obciążenia i nie powinna zależeć od wartości tego obciążenia. Tymczasem równoważny mimośród wyznaczony z warunku statycznego jest funkcją parametru \(\bar{\Lambda}\). Oznacza to, że nie może być interpretowany jako rzeczywista imperfekcja geometryczna konstrukcji, lecz jedynie jako parametr zastępczy odpowiadający określonemu stanowi obciążenia.
W praktyce projektowania, zwłaszcza konstrukcji żelbetowych, wpływ imperfekcji geometrycznych jest często reprezentowany przez zastępczy mimośród działania siły ściskającej $ M=N\,e$.
Przeprowadzona analiza wskazuje jednak, że takie podejście nie zachowuje ścisłej równoważności mechanicznej nawet w odniesieniu do podstawowego kryterium statycznego. Uzyskany mimośród jest bowiem zależny od poziomu obciążenia osiowego, a więc nie może być traktowany jako jednoznaczny odpowiednik rzeczywistej imperfekcji geometrycznej.
Łączna analiza kryterium przemieszczeniowego oraz statycznego prowadzi do ogólnego wniosku, że zastępcza siła pozioma stanowi bardziej fundamentalny model referencyjny imperfekcji niż mimośród działania siły osiowej. Mimośród może być wyznaczany jedynie jako parametr równoważny względem tej imperfekcji i tylko dla ściśle określonego kryterium oraz stanu pracy konstrukcji. Nie istnieje zatem pojedyncza wartość mimośrodu zdolna do zachowania równoważności mechanicznej w całym zakresie obciążeń i dla wszystkich kryteriów oceny odpowiedzi konstrukcji.

Rownoważność energetyczna zastępczego mimośrodu

Dla modelu z równoważną siłą poziomą \(H\) oraz modelu z równoważnym mimośrodem \(e\) warunek równoważności energetycznej można zapisać w postaci

\[ \Pi_H=\Pi_e. \tag{II.73}\label{II.73} \]

Po podstawieniu odpowiednich wektorów obciążenia oraz wykorzystaniu zależności (\ref{II.69}) funkcjonały całkowitej energii potencjalnej obu modeli można wyrazić wyłącznie za pomocą elementów odwrotnej macierzy sztywności. Ponieważ macierz ta została wyprowadzona analitycznie z wykorzystaniem ścisłych funkcji kształtu Livesleya, wszystkie całki energetyczne zostały już uwzględnione podczas wyznaczania funkcji statecznościowych
\(\Psi_1\), \(\Psi_2\), \(\Psi_3\) oraz \(\Psi_4\). W rezultacie otrzymuje się

\[ \cfrac{\Pi_e}{\Pi_H}= \cfrac{ 3e^{\,2}N^{2} \left( -3\Psi_2^{\,2} +\bar C_{\varphi}\Psi_1 +4\Psi_1\Psi_3 \right) } { H^{2}L^{2} \left[ \Psi_3\left(\bar C_{\varphi}+4\Psi_3\right) -\Psi_4^{\,2} \right] }.
\tag{II.74}\label{II.74} \]

Po wprowadzeniu bezwymiarowego mimośrodu

\[ \bar e=\cfrac{e}{\Delta}, \qquad \Delta=\cfrac{L}{n_L}, \qquad H=\cfrac{N}{n_L}, \tag{II.75}\label{II.75} \]

otrzymuje się współczynnik równoważności energetycznej

\[ \eta_E = \cfrac{\Pi_e}{\Pi_H} = \cfrac{ 3\bar e^{\,2} \left( -3\Psi_2^{\,2} +\bar C_{\varphi}\Psi_1 +4\Psi_1\Psi_3 \right) } { \Psi_3\left(\bar C_{\varphi}+4\Psi_3\right) -\Psi_4^{\,2}}. \tag{II.76}\label{II.76} \]

Wartość \(\eta_E=1\) odpowiada pełnej równoważności energetycznej obu modeli, natomiast dla \(\eta_E<1\) energia odkształcenia modelu z mimośrodem jest mniejsza od energii modelu odniesienia, zaś dla \(\eta_E>1\) – większa. Uwzględniając warunek (\ref{II.73}), otrzymuje się ścisły wzór określający bezwymiarowy mimośród równoważny

\[ \bar e_{eq}= \sqrt{ \cfrac{ \Psi_3\left(\bar C_{\varphi}+4\Psi_3\right) -\Psi_4^{\,2} } { 3\left( -3\Psi_2^{\,2} +\bar C_{\varphi}\Psi_1 +4\Psi_1\Psi_3 \right) } }. \tag{II.77}\label{II.77} \]

Równanie (\ref{II.77}) stanowi ścisłe rozwiązanie problemu równoważności energetycznej i jednoznacznie określa wartość zastępczego mimośrodu jako funkcję poziomu obciążenia \(\bar\Lambda\) oraz bezwymiarowej  podatności obrotowej podpory \(\bar C_{\varphi}\). Ponieważ funkcje \(\Psi_1\)-\(\Psi_4\) zależą wyłącznie od parametru \(\bar\Lambda\), równanie to definiuje pełną charakterystykę \(\bar e_{eq}=\bar e_{eq} \bar\Lambda,\bar C_{\varphi})\).

Na rys. II.7 przedstawiono rodziny charakterystyk współczynnika równoważności energetycznej \(\eta_E\) w funkcji bezwymiarowego mimośrodu \(\bar e_{eq}\). Parametrami wykresu są bezwymiarowy mnożnik obciążenia
\(\bar\Lambda=\Lambda/\Lambda_{cr}\) oraz bezwymiarowa podatność obrotowa podpory \(\bar C_{\varphi}\). Linie ciągłe odpowiadają \(\bar C_{\varphi}=8\), natomiast linie przerywane \(\bar C_{\varphi}=25\).

Pozioma prosta \(\eta_E=1\) wyznacza warunek pełnej równoważności energetycznej. Punkty przecięcia tej prostej z poszczególnymi krzywymi określają wartości zastępczego mimośrodu \(\bar e_{eq}\), dla których oba modele magazynują identyczną energię odkształcenia. Są to rozwiązania równania (\ref{II.77}), oznaczone na wykresie odpowiednimi wartościami. </p> <p> Z przedstawionych charakterystyk wynika, że wraz ze wzrostem parametru \(\bar\Lambda\), odpowiadającym zbliżaniu się konstrukcji do stanu krytycznego, wartość równoważnego mimośrodu maleje. Oznacza to, że narastające efekty drugiego rzędu powodują wzrost energii  dkształcenia
układu już przy coraz mniejszych mimośrodach. Jednocześnie zwiększenie sztywności obrotowej podpory prowadzi do wzrostu wartości \(\bar e_{eq}\), ponieważ sztywniejszy układ wymaga większego mimośrodu, aby zgromadzić tę samą energię odkształcenia. Otrzymane charakterystyki stanowią praktyczny nomogram umożliwiający szybkie wyznaczenie energetycznie równoważnego mimośrodu bez konieczności każdorazowego rozwiązywania równania (\ref{II.77}).

Szczególnie interesujące są dwie graniczne postacie rozwiązania (\ref{II.77}). Dla przypadku idealnie sztywnego zamocowania \(\bar C_{\varphi}\rightarrow\infty\) otrzymuje się

\[ \lim_{\bar C_{\varphi}\rightarrow\infty}\bar e_{eq} = \sqrt{\cfrac{\Psi_3}{3\Psi_1}}. \tag{II.78}\label{II.78} \]

Oznacza to, że dla dostatecznie dużej sztywności obrotowej wartość równoważnego mimośrodu zależy wyłącznie od funkcji statecznościowych, a tym samym wyłącznie od poziomu obciążenia
\(\bar{\Lambda}\).

Natomiast w granicy zaniku obciążenia osiowego \(\bar{\Lambda}\rightarrow0\), gdy \(\Psi_1=\Psi_2=\Psi_3=\Psi_4=1\), otrzymuje się prostą zależność

\[ \lim_{\bar{\Lambda}\rightarrow0}\bar e_{eq} = \cfrac{1}{\sqrt3} \sqrt{\cfrac{\bar C_{\varphi}+3} {\bar C_{\varphi}+1}}. \tag{II.79}\label{II.79} \]

Zależność (\ref{II.79}) pokazuje, że nawet przy pomijalnym wpływie siły osiowej wartość równoważnego mimośrodu zależy od podatności obrotowej podpory. W granicy idealnego utwierdzenia otrzymuje się \(\bar e_{eq}=1/\sqrt3\), natomiast dla podpór podatnych wartość ta  wzrasta zgodnie z równaniem (\ref{II.79}).

Po podstawieniu ścisłych funkcji stateczności Livesleya (\ref{II.17})–(\ref{II.20}) współczynnik równoważności energetycznej (\ref{II.76}) można zapisać wyłącznie jako funkcję parametrów
\(\bar{\Lambda}\), \(\bar C_{\varphi}\) oraz $\bar e$:/

\[ \eta_E= \cfrac{ \bar e^{\,2}\, \bar{\Lambda}\pi^2 \left[ \pi\sqrt{\bar{\Lambda}} \cos\!\left(\cfrac{\pi\sqrt{\bar{\Lambda}}}{2}\right) + 2\bar C_{\varphi} \sin\!\left(\cfrac{\pi\sqrt{\bar{\Lambda}}}{2}\right) \right]
} { 2\left(4\bar C_{\varphi} +\pi^2\bar{\Lambda}\right) \sin\!\left(\cfrac{\pi\sqrt{\bar{\Lambda}}}{2}\right) – 4\pi\bar C_{\varphi}\sqrt{\bar{\Lambda}} \cos\!\left(\cfrac{\pi\sqrt{\bar{\Lambda}}}{2}\right)}. \tag{II.80}\label{II.80}\]

Równanie (\ref{II.80}) stanowi bezpośrednią podstawę do sporządzenia charakterystyk przedstawionych na rys. II.7. Punkty przecięcia krzywych z poziomą prostą \(\eta_E=1\) wyznaczają wartości
bezwymiarowego mimośrodu równoważnego spełniające warunek pełnej równoważności energetycznej.

Wniosek o zastępczym mimośrodzie

Wniosek dla praktyki projektowej i rozwoju metod normowych.Zastępowanie imperfekcji przechyłowych konstrukcji geometrycznym mimośrodem działania siły osiowej nie zapewnia zachowania ścisłej równoważności mechanicznej. Nie istnieje jedna wartość mimośrodu zdolna do jednoczesnego odwzorowania przemieszczeń i sił wewnętrznych w całym zakresie pracy konstrukcji. Mimośród działania siły osiowej nie powinien być traktowany jako uniwersalny model zastępczy imperfekcji geometrycznych, lecz jedynie jako parametr pomocniczy, którego wartość wymaga kalibracji względem jednoznacznie zdefiniowanego modelu referencyjnego. Wyniki przeprowadzonych analiz wskazują, że rolę takiego modelu w sposób bardziej konsekwentny spełnia zastępcza siła pozioma, stanowiąca bezpośrednią przyczynę efektów drugiego rzędu i niezależna od aktualnego stanu obciążenia konstrukcji. Pokazane wyniki wyniki stanowią przesłankę do ponownej oceny metod projektowych wykorzystujących pojedynczy mimośród zastępczy jako reprezentację imperfekcji geometrycznych. W szczególności wskazują na potrzebę opracowania metod projektowania, w których równoważność imperfekcji będzie oparta na kryteriach mechanicznych, a nie wyłącznie na geometrycznej interpretacji mimośrodu działania siły ściskającej.

Analiza pręta w konstrukcji

Zgodnie z klasyfikacją przyjętą w normie [7] węzły konstrukcji klasyfikowane zsą e względu na sztywność jako: sztywne, podatne (półsztywne) oraz przegubowe. Granice pomiędzy poszczególnymi klasami wyznaczane są na podstawie początkowej sztywności obrotowej $S_{j,ini}$ j-tego węzła. Odpowiadające obszary klasyfikacji przedstawiono na rys. II.7.

W klasyfikacji normowej początkową sztywność obrotową węzła $S_{j,ini}$   ustala się na podstawie  sztywności zginania $ EI_b/L_b$, rozpatrywanego elementu konstrukcji $b= [e] $  gdzie $EI_b$ oznacza sztywność zginania elementu , a $L_b$ jego długość . Oznaczenia normowe odpowiadają oznczeniom $\bar C_{\varphi}$ przyjętym w artykule:

$ S_{j,ini}= \bar C_{\varphi}$, $ EI_b=EI,\qquad L_b=L.$

Współczynnik $k_b$ występujący na rys. II.7  zależy od rodzaju układu konstrukcyjnego. Zgodnie z normą dla ram stężonych  przyjmuje się $k_b=8$, natomiast dla ram niestężonych $k_b=25$. W tym drugim przypadku wymagane jest ponadto spełnienie warunku $ \cfrac{K_b}{K_c}\ge 0,1$,  gdzie $K_b$ oraz $K_c$ oznaczają odpowiednio sztywności belek i słupów analizowanej kondygnacji konstrukcji budynku.

Analizę nieliniowe pręta w konstrukcji przedstawiono na rys. II.8., Pręt jest  utwierdzony  sprężyście w stopie fundamentowej (lub w konstrukcji niższej kondygnacji) za pomocą sprężyny obrotowej o sztywności \(C_{1,\varphi}\). Górny koniec pręta połączony jest z konstrukcją stropu (lub przekrycia) poprzez sprężynę obrotową o sztywności \(C_{2,\varphi}\) oraz sprężynę poziomą o sztywności \(C_{2,\Delta}\). Sprężyny te stanowią model zastępczy wpływu pozostałej części konstrukcji na analizowany pręt. Odpowiadające im sztywności normowe wynoszą: $\bar C_{1,\varphi}=S_{1,\mathrm{ini}}$, $\bar C_{2,\varphi}=S_{2,\mathrm{ini}}$.

W analizie parametrycznej rozpatrzono wartości sztywności skrętnych $\bar C_{\varphi}=8,25$ – graniczne   pomiędzy wężłęm podatnym, a sztywnym$

Sztywność pozioma górnego węzła \(\bar C_{2,\Delta}\) zależy od globalnej sztywności poziomej całej konstrukcji, a w szczególności od liczby i sztywności słupów, belek oraz układu stężeń. Jej jednoznaczne wyznaczenie wymaga analizy całego układu konstrukcyjnego i wykracza poza zakres niniejszego opracowania. Z tego względu parametr \(\bar C_{2,\Delta}\) traktowany jest jako niezależny parametr analizy, którego wpływ na zachowanie pręta badany jest metodą analizy parametrycznej. Do obliczeń przyjęto reprezentatywne wartości

$\bar C_{2,\Delta}= 1,8,25,50$, obejmujące zakres od całkowitego braku bardo podatnego podparcia poziomego do praktycznie nieprzesuwnego węzła.

Takie ujęcie pozwala w sposób uniwersalny zastąpić wpływ całej otaczającej, rzeczywistej konstrukcji równoważnymi sztywnościami sprężystymi podpór. Dzięki temu możliwe jest zbadanie mechanizmu współpracy pręta z konstrukcją w postaci ogólnej, niezależnie od szczegółowego rozwiązania konstrukcyjnego. 

Wniosek:

Z analizy przedstawionej na rys  rys II.8 wynika, że : 

(1) Elementarny model pręta ze sprężystymi podporami ujawnia podstawowe mechanizmy odpowiedzialne za redystrybucję energii i powstawanie rezerwy stateczności konstrukcji,
(2) Dla małych sztywności podpory energia wyboczenia koncentruje się w analizowanym pręcie, dlatego pierwsza osobliwość jest silna. Wraz ze wzrostem sztywności podpory coraz większa część energii jest przejmowana przez układ podparcia. Mechanizm utraty stateczności pierwszej postaci staje się coraz słabszy, a jego wpływ na przemieszczenia maleje. W granicy sztywności dążącej do nieskończoności odpowiedź przechodzi do innego schematu statycznego, więc pierwsza osobliwość może przestać być obserwowalna.
(2) Wraz ze wzrostem bezwymiarowej sztywności podpory maleje intensywność pierwszej osobliwości odpowiedzi. Dla odpowiednio dużych wartości parametrów sztywności jej wpływ staje się praktycznie niezauważalny, co oznacza stopniowe wygaszanie mechanizmu utraty stateczności związanego z pierwszą postacią wyboczenia.

Konstrukcja z wymuszonymi przemieszczeniami

W poprzednich rozdziałach rozpatrywano zagadnienia, w których wymuszeniem konstrukcji były wyłącznie obciążenia zewnętrzne. Obecnie poddamy anlizie sytuację, w której  część przemieszczeń konstrukcji jest z góry określona i stanowi dane zadania. Do tej grupy należą między innymi imperfekcje przechyłowe, zadane przemieszczenia podpór, wymuszone obroty węzłów, błędy montażowe, osiadania podpór oraz inne wymuszenia kinematyczne. W klasycznej metodzie przemieszczeń wymuszenia takie najczęściej zastępowane są równoważnymi obciążeniami zewnętrznymi. Podejście to upraszcza analizę konstrukcji, jednak równoważność pomiędzy wymuszeniem kinematycznym a odpowiadającym mu układem obciążeń nie jest oczywista w przypadku analizy geometrycznie nieliniowej. Zachowanie zgodności momentów lub reakcji nie gwarantuje zgodności przemieszczeń, obrotów ani efektów drugiego rzędu. W pracy p rzyjęto odmienne podejście. Punktem wyjścia analizy jest rzeczywista konstrukcja z zadanymi wymuszeniami kinematycznymi, a nie model zastępczy z równoważnymi obciążeniami. Rozwiązanie otrzymane dla rzeczywistych wymuszeń będzie traktowane jako rozwiązanie odniesienia, względem którego oceniana będzie dokładność klasycznych procedur projektowych oraz modeli zastępczych opartych na równoważnych obciążeniach.

Niech globalny wektor stopni swobody konstrukcji zostanie podzielony na dwie części

\[ \{q\}= \begin{Bmatrix} q_n\\ q_w \end{Bmatrix}, \tag{II.81}\label{II.81} \]

gdzie $q_n$ oznacza wektor nieznanych przemieszczeń i obrotów, natomiast $q_w$ jest wektorem zadanych wymuszeń kinematycznych.

Po odpowiednim uporządkowaniu stopni swobody globalne równanie równowagi przyjmuje postać

\[ \begin{bmatrix}K_{nn} & K_{nw}\\ K_{wn} & K_{ww} \end{bmatrix} \begin{Bmatrix} q_n\\ q_w \end{Bmatrix}= \begin{Bmatrix} P_n\\ P_w
\end{Bmatrix}. \tag{II.82}\label{II.82} \]

Ponieważ wektor wymuszeń kinematycznych $q_w$ jest znany, jego wpływ można przenieść na prawą stronę równania, otrzymując 

\[ K_{nn}q_n=P_n-K_{nw}q_w. \tag{II.83}\label{II.83} \]

Po odwróceniu macierzy $K_{nn}$ otrzymuje się

\[ q_n= K_{nn}^{-1} \left( P_n-K_{nw}q_w \right), \tag{II.84}\label{II.84}\]

natomiast reakcje odpowiadające zadanym wymuszeniom wyznacza się z zależności

\[ R_w= K_{wn}q_n+ K_{ww}q_w – P_w.\tag{II.85}\label{II.85} \]

Przedstawione równania stanowią ogólną metodologię analizy konstrukcji z wymuszonymi przemieszczeniami i obrotami. Obejmują one wszystkie przypadki wymuszeń kinematycznych, niezależnie od ich fizycznego pochodzenia, takie jak imperfekcje przechyłowe, zadane przemieszczenia podpór, wymuszone obroty węzłów, mimośrodowe przyłożenie sił osiowych czy osiadania podpór.

Analizujemy dokładność zastępowania rzeczywistych imperfekcji przechyłowych równoważnymi obciążeniami poziomymi oraz wpływ kolejnych przybliżeń teorii konstrukcji na zachowanie tej równoważności. Ścisła teoria belki-słupa stanowi rozwiązanie odniesienia umożliwiające ilościową ocenę błędów modeli zastępczych oraz procedur obliczeniowych wykorzystywanych w praktyce projektowej.

Funkcjonał Lagrange’a belki-słupa  Bernoulli

Energia potencjalna belki-słupa w naprężeniach i przemieszczeniach

Przedstawimy rozwiązanie klasycznej belki Bernoulli-Euler obciążonej rozłożonym obciążeniem poprzecznym $q(x)$, a także ściskanej siłą osiową $N(x)$ oraz poddanej działaniu skupionych obciążeń poprzecznych $[V,H,M](x_k)$ zlokalizowanych w punkcie o współrzędnej $x_k$. Indeks $k$ oznacza, w zależności od kontekstu, miejsce przyłożenia podpory skupionej lub obciążenia skupionego.

Powszechnie znaną zależność na gęstość energii potencjalnej ciała sprężystego zapiszmy w postaci [9]:

\[ \Phi= \cfrac{1}{2E} \left\{ \sigma_x^2+\sigma_y^2+\sigma_z^2 -2\nu(\sigma_x\sigma_y+\sigma_x\sigma_z+\sigma_y\sigma_z) +2(1+\nu)(\tau_{xy}^2+\tau_{xz}^2+\tau_{yz}^2) \right\} \tag{II.94} \label{II.94} \]

gdzie w celu zwiększenia czytelności zastosowano symbole klasyczne $(x,y,z)$ zamiast zwykle stosowanego zapisu wskaźnikowego.

W przypadku zginania poprzecznego wzór ten znacznie upraszcza się. Ponieważ wpływ naprężeń $\sigma_y$, $\tau_{xz}$ oraz $\tau_{yz}$ jest pomijalnie mały, przyjmuje się, że różne od zera są jedynie naprężenia $\sigma_x$ oraz $\tau_{xy}$.

W podejściu klasycznym Bernoulli-Euler pomija się również wpływ naprężeń $\tau_{xy}$ na przemieszczenia i energię potencjalną. Wpływ naprężeń stycznych uwzględnimy w rozdziale dotyczącym belki Timoshenko. Przy poczynionych założeniach upraszczających energia potencjalna belki zginanej wynosi

\[ \Pi_\sigma = \iiint\limits_V \Phi\,dV = \cfrac{1}{2E} \iiint\limits_V \sigma_x^2\,dV \tag{II.95} \label{II.95} \]

Energię (\ref{II.95}) wyrazimy przez siły przekrojowe, a następnie przez przemieszczenia. Zachodzi

\[ \sigma_x = -\cfrac{N}{A} + \cfrac{M(x)}{I_y}\,z \tag{II.96} \label{II.96} \]

gdzie: $A$ – pole przekroju pręta, $I_y$ – moment bezwładności przekroju względem osi głównej $y$, $z$ – odległość punktu przekroju od osi głównej. Dalej moment bezwładności oznaczamy bez indeksu.

W przypadku przeważającego zginania wpływ energii od siły osiowej $N$ na minimum funkcjonału jest zwykle niewielki. Jednak wraz ze wzrostem stopnia wykorzystania nośności wyboczeniowej $\bar{\Lambda}$ wpływ siły ściskającej staje się dominujący i prowadzi do efektów drugiego rzędu oraz utraty stateczności. Zjawiska te zostaną uwzględnione w dalszej części opracowania.

Całkę w równaniu (\ref{II.95}) przekształcimy tak, aby otrzymać energię wyrażoną przez przemieszczenia.

Ponieważ zachodzi związek pomiędzy momentem zginającym a krzywizną

\[ M(x) = -EI\,w”(x) \tag{II.97} \label{II.97} \]

więc mamy

  \[ \Pi_\sigma = \iiint\limits_V \cfrac{1}{2E} \left[ -\cfrac{EI\,w”(x)}{I}\,z \right]^2 dV = \cfrac{E}{2} \iiint\limits_V [w”(x)]^2 z^2\,dV = \]

  \[ \cfrac{E}{2} \int\limits_0^L [w”(x)]^2\,dx \cdot \iint\limits_A z^2\,dA = \cfrac{EI_y}{2} \int\limits_0^L [w”(x)]^2\,dx \tag{II.98} \label{II.98} \]

Praca zewnętrznych sił poprzecznych

Praca zewnętrznych sił poprzecznych $q(x)$ oraz momentów zginających $m(x)$ rozłożonych na długości belki wynosi

\[L_q= \int\limits_0^L \left[ q(x)\,w(x) + m(x)\,w'(x) \right] dx \tag{II.99} \label{II.99} \]

Obciążenia skupione zlokalizowane w punkcie $x_k$ można zapisać za pomocą funkcji Diraca w postaci

\[ Q_k(x)=Q_k\,\delta(x-x_k), \qquad M_k(x)=M_k\,\delta(x-x_k) \tag{II.100} \label{II.100} \]

gdzie $\delta(x-x_k)$ jest funkcją Diraca reprezentującą siłę lub moment skupiony przyłożony w punkcie $x_k$. Formalnie funkcja Diraca  dla dowolnej funkcji ciągłej $f(x)$. pełnia zależność

\[ \int\limits_{-\infty}^{+\infty} f(x)\, \delta(x-x_k)\,dx = f(x_k) \tag{II.101} \label{II.101} \]

Własność (\ref{II.101}) powoduje, że całkowanie obciążenia skupionego po długości pręta prowadzi do wartości funkcji przemieszczeń obliczonej dokładnie w punkcie przyłożenia siły lub momentu.

Po uwzględnieniu obciążeń skupionych równanie pracy sił zewnętrznych przyjmuje postać

  \[ L_q= \int\limits_0^L \left[ q(x)\,w(x) + m(x)\,w'(x) \right] dx + \sum_k \left[ Q_k\,w(x_k) + M_k\,w'(x_k) \right] \tag{II.102} \label{II.102} \]

Pierwsza całka opisuje pracę obciążeń ciągłych rozłożonych na długości elementu, natomiast suma uwzględnia pracę sił skupionych i momentów skupionych przyłożonych w punktach dyskretnych.

W dalszej części opracowania zapis z wykorzystaniem funkcji Diraca będzie wykorzystywany do reprezentacji obciążeń skupionych, podpór sprężystych oraz innych oddziaływań skoncentrowanych w pojedynczych punktach konstrukcji.

Efekt drugiego rzędu P-Δ. Praca siły ściskającej N

W przypadku belki obciążonej siłą osiową $N$ pojawia się dodatkowy składnik energii związany z geometryczną zmianą położenia osi pręta. Jest to tzw. efekt drugiego rzędu $P-\Delta$, wynikający z pracy wykonywanej przez siłę osiową podczas ugięcia elementu. W celu wyznaczenia tej pracy w niniejszym artykule rozważmy elementarny odcinek belki pokazany na rys. II.2.

Analizujemy odcinek $dx$ belki. Ugięcie powoduje przemieszczenie punktów końcowych do położenia $dx’$. Pozioma odległość $\Delta(dx)$ pomiędzy punktem przed i po ugięciu, na której pracuje siła osiowa $N$, wynosi

\[ \Delta(dx) = dx-dx’ = dx-\sqrt{(dx)^2-\left(\cfrac{dw}{dx}dx\right)^2} = dx\left[1-\sqrt{1-(w’)^2}\right] \tag{II.103}\label{II.103} \]

Po rozwinięciu w szereg Maclaurina ostatniego pierwiastka w równaniu ( \ref{II.103}) i pozostawieniu pierwszego nieliniowego wyrazu szeregu otrzymujemy aproksymację stanowiącą podstawę teorii drugiego rzędu całej mechaniki konstrukcji

\[ \Delta(dx) \approx dx \left[ 1- \left( 1-\cfrac{1}{2}(w’)^2 +\cfrac{1}{8}(w’)^4 -\ldots \right) \right] \approx \cfrac{1}{2}(w’)^2\,dx \tag{II.104} \label{II.104} \]

W rezultacie praca osiowej siły ściskającej wynosi

  \[ L_N = \int\limits_0^L N\,\Delta(dx) = \int\limits_0^L \cfrac{N}{2} [w'(x)]^2\,dx \tag{II.105} \label{II.105} \]

Praca sił odporu podłoża i podpór sprężystych

Siły odporu $q_z(x)$ są reakcją podłoża $r(x)$ na nacisk belki. Z definicji zależą one bezpośrednio od przemieszczenia $w(x)$ zgodnie z zależnością

\[q_z(x) = -r(x) = -c(x)\,w(x) \tag{II.106} \label{II.106} \]

gdzie znak minus wynika z faktu, że siły odporu są przeciwnie skierowane do ugięcia.

Na rys. II.4pokazano odpór pionowy $v$ podłoża o stałej sprężystości $c^v$ (małe $c$) oraz odpór obrotowy $m$ podłoża o stałej sprężystości $c^m$. Analogicznie można zdefiniować odpór poziomy $h$ podłoża o stałej sprężystości $c^h$.

Reakcje skupionych podpór sprężystych otrzymujemy przez połączenie zależności (\ref{II.106}) z zapisem obciążeń skupionych za pomocą funkcji Diraca (\ref{II.100}). Dla podpory pionowej i obrotowej otrzymujemy

\[ q_k = -{C_k}^V\,w(x)\,\delta(x-x_k), \qquad m_k = -{C_k}^M\,w'(x)\,\delta(x-x_k) \tag{II.107} \label{II.107} \]

Nadając przyrosty przemieszczenia $w(x)$, nadajemy jednocześnie przyrosty siłom odporu. Zapisując pracę wirtualną sił odporu w postaci

\[ \delta L_c = -\int\limits_0^L [-r(x)] \,\delta w(x)\,dx = -\int\limits_0^L c\,w(x)\,\delta w(x)\,dx \tag{II.108} \label{II.108} \]

widzimy, że nie ma możliwości wyłączenia znaku wariacji przed całkę. Można to zrobić dopiero po zapisaniu całki w postaci

\[ \delta L_c = -\int\limits_0^L c\cdot\cfrac{1}{2} \, \delta[w^2(x)] \,dx = -\delta \int\limits_0^L \cfrac{c}{2} \,w^2(x)\,dx \tag{II.109} \label{II.109} \]

Po opuszczeniu znaku wariacji otrzymujemy wyrażenie na energię potencjalną odkształcenia podłoża sprężystego

 \[ L_c = \int\limits_0^L \cfrac{c}{2} \,w^2(x)\,dx \tag{II.110} \label{II.110} \]

Analogiczne wyrażenia dla skupionych podpór sprężystych otrzymujemy po podstawieniu zależności (\ref{II.107}) do wzoru (\ref{II.110}).\]

Funkcjonał Lagrange’a belki Bernoulli

Ostatecznie funkcjonał Lagrange’a dla belki Bernoulli ściskanej siłą osiową $N$, na podłożu sprężystym i podporach sprężystych zapiszemy w postaci:

  \[ \Pi= \int\limits_0^L \left[ \cfrac{EI_y}{2}[w”(x)]^2 – \cfrac{N(x)}{2}[w'(x)]^2 + \cfrac{c^v(x)}{2}[w(x)]^2 – q(x)w(x) \right] dx + \sum_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_k w'(x_k) \right] \tag{II.111}\label{II.111} \]

Poszukiwanie punktu stacjonarnego funkcjonału (\ref{II.102}) można przeprowadzić dowolną metodą wariacyjną, na przykład metodą Ritza, również w postaci jej implementacji w metodzie elementów skończonych.

Minimalizacja funkcjonału Lagranga metodą Ritza

W klasycznej metodzie Ritza poszukuje się funkcji realizującej minimum funkcjonału (\ref{II.111}) w przyjętej przestrzeni funkcji aproksymujących. Najczęściej przyjmuje się aproksymację w postaci kombinacji liniowej funkcji bazowych

\[ w(x)=\varphi_0(x)+\sum_{i=1}^{n}a_i\,\varphi_i(x), \tag{II.84} \label{II.84} \]

gdzie \(a_i\) są nieznanymi współczynnikami aproksymacji (stopniami swobody metody Ritza), natomiast \(\varphi_i(x)\) oznaczają funkcje bazowe.

W przedstawionym ujęciu klasyczna metoda Ritza zostanie rozszerzona na układy złożone z elementów ciągłych oraz dyskretnych, takich jak podpory sprężyste, sprężyny translacyjne i obrotowe. W takich zagadnieniach całkowite pole przemieszczeń może zawierać, oprócz części opisującej odkształcenie elementu ciągłego, również składnik odpowiadający ruchowi elementów dyskretnych konstrukcji. W konsekwencji oprócz współczynników aproksymacji \(a_i\) w funkcjonale mogą występować dodatkowe niewiadome, opisujące przemieszczenia lub obroty podpór sprężystych. Wielkości te nie stanowią współczynników rozwinięcia funkcji bazowych, lecz są dodatkowymi stopniami swobody konstrukcji wyznaczanymi równocześnie z współczynnikami metody Ritza na podstawie zasady minimum energii potencjalnej.

Należy zauważyć, że oznaczenie $\varphi_i(x)$ funkcji bazowych jest podobne do oznaczenia kąta obrotu przekroju $\varphi$. Znaczenie symbolu wynika każdorazowo z kontekstu.

Dobór funkcji bazowych

Jednym z najistotniejszych etapów metody Ritza jest dobór funkcji bazowych, od którego zależą zarówno dokładność aproksymacji, jak i szybkość zbieżności rozwiązania. Funkcje te mogą być wybierane w zasadzie dowolnie, pod warunkiem spełnienia dwóch podstawowych wymagań:

• kinematycznej dopuszczalności, tzn. spełnienia wszystkich istotnych (geometrycznych) warunków brzegowych analizowanego zagadnienia,
• liniowej niezależności, tzn. aby żadnej z funkcji zbioru nie można było przedstawić jako kombinacji liniowej pozostałych funkcji bazowych.

Należy podkreślić, że wymaganie kinematycznej dopuszczalności odnosi się wyłącznie do geometrycznych więzów elementów ciągłych.
Warunki wynikające z obecności elementów dyskretnych, takich jak podpory sprężyste lub sprężyny obrotowe, nie muszą być uwzględniane przez funkcje bazowe. Wpływ tych elementów jest opisywany przez odpowiednie składniki funkcjonału energii, natomiast odpowiadające im przemieszczenia lub obroty stanowią dodatkowe niewiadome wariacyjne.
Takie ujęcie pozwala zachować tę samą bazę funkcji dla różnych konfiguracji podpór, a zmiana modelu konstrukcji sprowadza się jedynie do odpowiedniej modyfikacji funkcjonału całkowitej energii potencjalnej.

Liniowa niezależność funkcji bazowych oznacza, że żadnej z funkcji \(\varphi_i(x)\) nie można przedstawić jako kombinacji liniowej pozostałych funkcji zbioru. Warunek ten zapisuje się w postaci

\[ \sum_{i=1}^{n}c_i\varphi_i(x)=0, \tag{II.113} \label{II.113}, \]

przy czym równość ta zachodzi dla każdego punktu rozważanego przedziału wtedy i tylko wtedy, gdy

\[ c_1=c_2=\ldots=c_n=0. \tag{II.114} \label{II.114} \]

Warunek liniowej niezależności gwarantuje, że każda funkcja bazowa wnosi do aproksymacji nową informację, a otrzymany układ równań Lagrange’a-Ritza posiada jednoznaczne rozwiązanie. W praktyce odpowiedni dobór funkcji modulujących ma zasadniczy wpływ nie tylko na dokładność rozwiązania, lecz również na szybkość zbieżności metody oraz uwarunkowanie numeryczne macierzy sztywności.

Klasyczna metoda Ritza

W klasycznej metodzie Ritza najwygodniej jest rozdzielić oba wymagania( kinemtycznej dopuszczalności i linioerj niezlażnosci) , konstruując funkcje bazowe jako iloczyn funkcji generatorowej oraz funkcji modulujących:

\[ \varphi_i(\xi)=\phi(\xi)\,\psi_i(\xi),\qquad i=1,2,\ldots,n. \tag{II.115} \label{II.115} \]

Funkcja generatorowa $\phi(\xi)$ jest wyznaczana w taki sposób, aby spełniała wszystkie zadane warunki kinematyczne na brzegach obszaru. Odpowiada ona wyłącznie za zapewnienie kinematycznej dopuszczalności całej aproksymacji. Natomiast funkcje modulujące $\psi_i(\xi)$ służą do wzbogacania przestrzeni aproksymacyjnej i powinny być wzajemnie liniowo niezależne oraz dostatecznie gładkie. Dzięki takiej konstrukcji każda funkcja bazowa automatycznie spełnia wymagane warunki brzegowe, niezależnie od wyboru funkcji modulujących.

Metoda ta ma charakter całkowicie ogólny i może być stosowana do zagadnień jedno-, dwu- i trójwymiarowych. W zależności od rodzaju analizowanego problemu funkcjami modulującymi mogą być wielomiany, funkcje trygonometryczne, funkcje ortogonalne (Legendre’a, Czebyszewa, Hermite’a), funkcje sklejane (B-splajny), funkcje falowe (wavelets) lub inne układy funkcji zapewniające dobrą zbieżność aproksymacji. O wyborze konkretnego zbioru decydują charakter rozwiązania, oczekiwana dokładność oraz własności numeryczne otrzymywanego układu równań.

Funkcję aproksymującą (\ref{II.84}) można zapisać w postaci

\[ w(x)=\varphi_0(x)+\sum_{i=1}^{n}a_i\varphi_i(x). \tag{II.116} \label{II.116} \]

Funkcja $ \varphi_0(x)$  służy do uwzględnienia niejednorodnych warunków brzegowych, natomiast funkcje bazowe $ \varphi_i(x)$  spełniają odpowiadające im jednorodne warunki kinematyczne. W przypadku tylko jednorodnych warunków brzegowych można przyjąć $\varphi_0(x)=0,$, co oznacza, że aproksymacja jest budowana wyłącznie z funkcji bazowych spełniających wymagania kinematycznej dopuszczalności.

Uwagi dotyczące zakresu stosowalności klasycznej metody Ritza

Klasyczna postać metody Ritza zakłada, że wszystkie warunki kinematyczne konstrukcji są znane przed rozpoczęciem procesu aproksymacji.
W takim przypadku funkcja \(\varphi_0(x)\), uwzględniająca niejednorodne warunki brzegowe, jest funkcją zadaną, natomiast niewiadomymi pozostają wyłącznie współczynniki aproksymacji $a_i$.
Takie sformułowanie jest wystarczające dla większości klasycznych zagadnień mechaniki konstrukcji, w których przemieszczenia i obroty podpór są znane lub wynikają z idealnych więzów geometrycznych. 
W praktyce inżynierskiej często spotyka się jednak konstrukcje współpracujące z elementami dyskretnymi, takimi jak podpory sprężyste, sprężyny translacyjne i obrotowe, podatne węzły ram, łożyska, a także zastępcze modele fragmentów konstrukcji o skończonej podatności. W takich przypadkach przemieszczenia lub obroty punktów podparcia nie są znane a priori, lecz stanowią część poszukiwanego rozwiązania.
Powstaje wówczas naturalne pytanie, czy klasyczna metoda Ritza może zostać rozszerzona w taki sposób, aby oprócz współczynników aproksymacji wyznaczać również przemieszczenia i obroty elementów dyskretnych. Innymi słowy, czy niewiadomymi wariacyjnymi mogą być nie tylko współczynniki rozwinięcia funkcji bazowych, lecz również wybrane stopnie swobody konstrukcji wynikające z jej rzeczywistej podatności.
W kolejnych punktach zostanie przedstawione takie uogólnienie klasycznej metody Ritza. Zachowuje ono wszystkie podstawowe założenia metody wariacyjnej, rozszerzając jedynie zbiór niewiadomych o dodatkowe stopnie swobody opisujące współpracę elementów ciągłych z elementami dyskretnymi. Dzięki temu możliwe staje się jednolite modelowanie belek i ram współpracujących z podporami sprężystymi, sprężynami obrotowymi oraz innymi elementami podatnymi bez konieczności zmiany przyjętej bazy funkcji aproksymacyjnych.

Rozszerzona metoda Ritza

W przypadkach konstrukcji ciągłych współpracujących z elementami dyskretnymi przemieszczenia lub obroty punktów podparcia nie są znane przed rozwiązaniem zadania, lecz stanowią część poszukiwanego rozwiązania.  Prowadzi to do naturalnego rozszerzenia klasycznej metody Ritza.

Zachowując niezmienioną postać funkcjonału całkowitej energii potencjalnej oraz klasyczne założenia wariacyjne, rozszerza się jedynie zbiór niewiadomych. Oprócz współczynników aproksymacji $a_i$ do zbioru zmiennych wariacyjnych wprowadza się dodatkowe stopnie swobody opisujące przemieszczenia lub obroty elementów dyskretnych konstrukcji. W ogólnym przypadku funkcję aproksymującą można zapisać w postaci

\[ w(x)=\varphi_0 (x,\mathbf{q}_d) +  \sum_{i=1}^{n}a_i\,\varphi_i(x), \tag{II.117} \label{II.117} \]

gdzie:
$\mathbf{q}_d= \{ q_1,q_2,\ldots,q_m\}$ –  wektor dodatkowych stopni swobody elementów dyskretnych,
$\varphi_0(x,\mathbf{q}_d)$  funkcja opisujaca opisująca wpływ dodatkowuyych stopni swobody na pole przemieszczeń.

W odróżnieniu od klasycznej metody Ritza funkcja $ \varphi_0(x,\mathbf{q}_d)$ nie jest funkcją zadaną, lecz zależy od niewiadomych parametrów konstrukcji i jest wyznaczana równocześnie ze współczynnikami aproksymacji $a_i$ z warunku stacjonarności funkcjonału energii.  W rezultacie poszukiwanymi niewiadomymi stają się zarówno

współczynniki aproksymacji $a_1,a_2,\ldots,a_n,$
jak również dodatkowe stopnie swobody konstrukcji  $ q_1,q_2,\ldots,q_m.$

Oznacza to, że przemieszczenia i obroty elementów dyskretnych są wyznaczane w tym samym procesie wariacyjnym co współczynniki aproksymacji, bez konieczności ich wcześniejszego zadawania lub wyznaczania z odrębnych zależności statycznych.

Układ kanoniczny równań Lagranga-Ritza 

Podstawiając aproksymację (\ref{II.115}) do funkcjonału (\ref{II.111}), otrzymujemy funkcjonał całkowitej energii potencjalnej $\Pi$, który staje się funkcją zmiennych wariacyjnych.

Warunek konieczny istnienia minimum funkcjonału polega na zaniku jego pochodnych cząstkowych względem wszystkich zmiennych wariacyjnych. Otrzymujemy zatem układ równań

\[ \cfrac{\partial\Pi}{\partial \eta_k}=0, \qquad k=1,2,\ldots, r, \tag{II.118} \label{II.118} \]

gdzie:
$ r = n + m.$ – liczba zmiennych wariayjnych $\eta_k$

Równania (\ref{II.118}) stanowią kanoniczny układ równań Lagranga-Ritza (uogólniony na dla układu ciągło-dyskretny).

Po wykonaniu różniczkowania funkcjonału zgodnie z warunkiem stacjonarności (\ref{II.118}) otrzymuje się układ liniowych równań algebraicznych

\[ \sum_{l=1}^{\,n+m}\delta_{kl}\, \eta_l = \Delta_{Fk}, \qquad k=1,2,\ldots,n+m, \tag{II.119} \label{II.119} \]

gdzie:
$ \eta_l= \begin{cases} a_l, & l=1,2,\ldots,n,\\[2mm] q_{ \, l-n}, & l=n+1,n+2,\ldots,n+m, \end{cases}$ – zbiór wszystkich niewiadomych wariacyjnych,
$\eta _l$ – $l$-ta zmienna wariacyjna obejmująca zarówno współczynniki aproksymacji, jak i dodatkowe stopnie swobody elementów dyskretnych,
$a_i$ – $i$ -ty współczynnik aproksymacji metody Ritza,
$q_j$ – $j$ -ty dodatkowy stopień swobody elementu dyskretnego (np. przemieszczenie podpory sprężystej lub obrót zamocowania podatnego),
$\delta_{kl}$ – współczynnik układu równań wynikający z różniczkowania funkcjonału całkowitej energii potencjalnej,
$\Delta_{Fk}$ – $k$-ta składowa wektora obciążeń uogólnionych,
$n$ – liczba współczynników aproksymacji,
$m$ – liczba dodatkowych stopni swobody elementów dyskretnych,
$r=n+m$ – całkowita liczba zmiennych wariacyjnych,
$k,l$ – indeksy numerujące równania oraz odpowiadające im zmienne wariacyjne.

W przypadku klasycznej metody Ritza nie występują dodatkowe stopnie swobody elementów dyskretnych (\(m=0\)). Wówczas układ (\ref{II.119}) redukuje się do klasycznego układu równań Lagranga-Ritza.

Współczynniki układu równań określone są zależnością

\[ \delta_{kl} \stackrel{\rm def}{=} \cfrac{\partial^2\Pi} {\partial\eta_k\,\partial\eta_l}, \qquad k,l=1,2,\ldots,r, \tag{II.120} \label{II.120} \]

W zależności od rodzaju zmiennych wariacyjnych współczynniki
\(\delta_{kl}\) mogą odpowiadać zarówno klasycznym współczynnikom metody Ritza, jak również sztywnościom elementów dyskretnych oraz współczynnikom sprzężenia pomiędzy elementami ciągłymi i dyskretnymi.

Wyrazy wolne wynoszą:

\[ \Delta_{Fk} \stackrel{\rm def}{=} – \left. \cfrac{\partial\Pi}{\partial\eta_k} \right|_{\eta=0}, \qquad k=1,2,\ldots,r. \tag{II.121} \label{II.121} \]

W celu ujednolicenia zapisu klasycznych oraz rozszerzonych równań Lagranga-Ritza wprowadza się funkcje

\[ N_k(x) \stackrel{\rm def}{=} \cfrac{\partial w(x)}{\partial\eta_k}, \qquad k=1,2,\ldots,r. \tag{II.122} \label{II.122}\]

Funkcje \(N_k(x)\) opisują wpływ poszczególnych zmiennych wariacyjnych na pole przemieszczeń i stanowią uogólnienie klasycznych funkcji bazowych metody Ritza.

Dla współczynników aproksymacji $\eta_k=a_i $  otrzymuje się
$N_k(x)=\varphi_i(x),$

natomiast dla dodatkowych stopni swobody elementów dyskretnych $ \eta_k=q_j$ otrzymuje się
$ N_k(x) = \cfrac{\partial\varphi_0(x,\mathbf q_d)} {\partial q_j}. $

Oznacza to, że funkcje odpowiadające elementom ciągłym oraz dyskretnym są traktowane w jednakowy sposób i różnią się jedynie interpretacją mechaniczną.

Korzystając z definicji (\ref{II.125}) współczynniki układu równań można zapisać w jednolitej postaci

\[ \delta_{kl} = \int_0^L \left[ EI_y \left( N_k”N_l” + \cfrac{EI_y}{GA_z} N_k”’N_l”’ \right) – N\,N_k’N_l’ + c^V N_kN_l \right]dx + \sum_p \left[ C_p^V N_k(x_p)N_l(x_p) + C_p^M
N_k'(x_p)N_l'(x_p) \right]. \tag{II.123} \label{II.123} \]

\[ \Delta_{Fk} = \int_0^L q(x)N_k(x)\,dx + \sum_p V_pN_k(x_p) + \sum_p M_pN_k'(x_p). \tag{II.124} \label{II.124}\]

W przypadku klasycznej metody Ritza, gdy jedynymi zmiennymi wariacyjnymi są współczynniki aproksymacji $ \eta_i\equiv a_i, $ zależności (\ref{II.120}) i (\ref{II.121}) przyjmują znaną postać

\[ \delta_{ij} = \int_{0}^{L} \left[ EI_y \left( \varphi_i”\varphi_j” + \cfrac{EI_y}{GA_z} \varphi_i”’\varphi_j”’ \right) – N\varphi_i’\varphi_j’ + c^V\varphi_i\varphi_j \right]dx +
\sum_k \left[ C_k^V\varphi_i(x_k)\varphi_j(x_k) + C_k^M \varphi_i'(x_k)\varphi_j'(x_k) \right], \tag{II.125} \label{II.125} \]

\[ \Delta_{Fj} = \int_0^L q(x)\varphi_j(x)\,dx + \sum_kV_k\varphi_j(x_k) + \sum_kM_k\varphi_j'(x_k). \tag{II.126} \label{II.126} \]

Zmienne wariacyjne $\eta_i$ stanowią uogólnione współrzędne przyjętej aproksymacji. Obejmują one zarówno współczynniki aproksymacji metody Ritza \(a_i\), opisujące odkształcenie elementów ciągłych, jak również dodatkowe stopnie swobody $q_j$, opisujące przemieszczenia lub obroty elementów dyskretnych konstrukcji. W szczególnym przypadku funkcji kształtu stosowanych w metodzie elementów skończonych zmienne te mogą być interpretowane jako przemieszczenia i obroty węzłów elementu, tj. przemieszczenia podłużne \(u_i\), poprzeczne \(w_i\) oraz obroty $\varphi_i$.
Równania (\ref{II.119}) stanowią kanoniczny układ równań Lagranga-Ritza dla układów ciągło-dyskretnych. W przypadku klasycznej metody Ritza, gdy nie występują dodatkowe stopnie swobody elementów dyskretnych (\(m=0\)), redukują się one do klasycznego układu równań Lagranga-Ritza.

Struktura otrzymanego układu jest analogiczna do równań stosowanych w metodzie przemieszczeń oraz metodzie elementów skończonych, co wskazuje na ich wspólne podstawy energetyczne wynikające z zasady minimum całkowitej energii potencjalnej. Rozwiązanie układu pozwala wyznaczyć wszystkie zmienne wariacyjne $\eta_i$, a więc zarówno współczynniki aproksymacji metody Ritza, jak i dodatkowe stopnie swobody elementów dyskretnych.
Stanowi to teoretyczne uzasadnienie tych metod obliczeniowych oraz wyjaśnia ich wspólne podstawy energetyczne. Z układu (\ref{II.119}) wyznacza się współczynniki Ritza $\eta_i$, które najczęściej mają znaczenie fizyczne przemieszczeń węzłowych elementów skończonych. 

Macierzowe ujęcie równań Lagranga-Ritza

Układ równań (\ref{II.119}) można przedstawić w zwartej postaci macierzowej

\[ [\delta]\{\eta\}=\{\Delta_F\}, \tag{II.127} \label{II.127} \]

gdzie

\[ [\delta]= \begin{bmatrix} \delta_{11} & \delta_{12} & \cdots & \delta_{1r}\\ \delta_{21} & \delta_{22} & \cdots & \delta_{2r}\\ \vdots & \vdots & \ddots & \vdots\\ \delta_{r1} & \delta_{r2} & \cdots & \delta_{rr} \end{bmatrix}, \tag{II.128} \label{II.128} \]

\[ \{\eta\}= \begin{Bmatrix} \eta_1\\ \eta_2\\ \vdots\\ \eta_r \end{Bmatrix}, \qquad r=n+m, \tag{II.129} \label{II.129} \]

\[ \{\Delta_F\}= \begin{Bmatrix} \Delta_{F1}\\ \Delta_{F2}\\ \vdots\\ \Delta_{Fr} \end{Bmatrix}. \tag{II.130} \label{II.130} \]

Macierz \([\delta]\) jest uogólnioną macierzą sztywności metody Ritza. Jej elementy opisują wzajemne sprzężenie wszystkich zmiennych wariacyjnych poprzez całkowitą energię potencjalną układu. Wektor \(\{\eta\}\) zawiera wszystkie niewiadome wariacyjne, obejmujące zarówno współczynniki aproksymacji, jak i dodatkowe stopnie swobody elementów dyskretnych, natomiast wektor $ \{\Delta_F\}$ reprezentuje odpowiadające im obciążenia uogólnione.

Jeżeli funkcjonał całkowitej energii potencjalnej jest dodatnio określony, macierz \([\delta]\) jest symetryczna i nieosobliwa, dzięki czemu rozwiązanie układu jest jednoznaczne. Zmienne wariacyjne wyznacza się z zależności

\[ \{\eta\} = [\delta]^{-1} \{\Delta_F\}. \tag{II.131} \label{II.131} \]

Po wyznaczeniu wektora $ \{\eta\}\$ pole przemieszczeń konstrukcji wyznacza się z przyjętej funkcji aproksymującej

\[ w(x) = \varphi_0(x,\mathbf q_d) + \sum_{i=1}^{n} a_i\varphi_i(x), \tag{II.132} \label{II.132} \]

przy czym współczynniki \(a_i\) oraz dodatkowe stopnie swobody \(\mathbf q_d\) stanowią odpowiednie składowe wektora \(\{\eta\}\).

Tak sformułowany układ równań jest formalnie identyczny z układami równań równowagi stosowanymi w metodzie przemieszczeń oraz metodzie elementów skończonych

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

Rolę macierzy sztywności konstrukcji \([K]\) pełni tutaj macierz \([\delta]\), rolę wektora niewiadomych \(\{q\}\) pełni wektor zmiennych wariacyjnych \(\{\eta\}\), natomiast rolę wektora obciążeń \(\{P\}\) pełni wektor obciążeń uogólnionych \(\{\Delta_F\}\).

Oznacza to, że metoda elementów skończonych może być interpretowana jako szczególny przypadek metody Ritza, w którym funkcje aproksymujące przyjmują postać funkcji kształtu elementów skończonych, natomiast niewiadome wariacyjne odpowiadają stopniom swobody węzłów elementów.

Interpretacja macierzy współczynników metody Ritza

Macierz współczynników \([\delta]\), zdefiniowana zależnością (\ref{II.120}), jest hesjanem funkcjonału całkowitej energii potencjalnej względem zmiennych wariacyjnych \(\eta_i\). Jej elementy opisują wpływ poszczególnych składników energii potencjalnej na sprzężenie pomiędzy zmiennymi wariacyjnymi. W zależności od rodzaju analizowanego zagadnienia mogą one odpowiadać energii odkształcenia elementów ciągłych, energii elementów dyskretnych oraz wzajemnemu sprzężeniu obu tych grup elementów.

W przypadku klasycznej liniowej teorii sprężystości jedynymi zmiennymi wariacyjnymi są współczynniki aproksymacji $ \eta_i \equiv a_i,$

a elementy macierzy współczynników mają postać

\[ \delta_{ij} = EI \int_0^L \phi_i”(x)\phi_j”(x)\,dx. \tag{II.134} \label{II.134} \]

Oznacza to, że macierz współczynników pokrywa się z macierzą sztywności sprężystej

\[ [\delta]=[K_e], \tag{II.135} \label{II.135} \]

gdzie \([K_e]\) jest macierzą sztywności odpowiadającą energii odkształcenia przy zginaniu belki.

W ogólnym przypadku funkcjonał całkowitej energii potencjalnej może zawierać oprócz energii odkształcenia przy zginaniu również składniki odpowiadające energii geometrycznej, energii sprężystego podłoża oraz energii elementów dyskretnych, takich jak podpory lub zamocowania sprężyste. Wówczas elementy macierzy współczynników można przedstawić w postaci

\[ \delta_{kl} = k_{e,kl} + k_{g,kl} + k_{c,kl} + k_{C,kl}, \tag{II.136} \label{II.136} \]

gdzie

\[ k_{e,kl} = EI \int_0^L N_k”(x)N_l”(x)\,dx, \tag{II.137} \label{II.137} \]

\[ k_{g,kl} = – N \int_0^L N_k'(x)N_l'(x)\,dx, \tag{II.138} \label{II.138} \]

\[ k_{c,kl} = \int_0^L C(x)\, N_k(x)N_l(x)\,dx, \tag{II.139} \label{II.139} \]

\[ k_{C,kl} = \sum_{p=1}^{n_s} C_p\, N_k(x_p)N_l(x_p), \tag{II.140} \label{II.140} \]

przy czym

\[ N_k(x) = \cfrac{\partial w(x)}{\partial\eta_k} \]

oznacza funkcję odpowiadającą \(k\)-tej zmiennej wariacyjnej.

Dla współczynników aproksymacji metody Ritza zachodzi

\[ N_k(x)=\phi_k(x), \]

natomiast dla dodatkowych stopni swobody elementów dyskretnych

\[ N_k(x) = \cfrac{\partial\varphi_0(x,\mathbf q_d)} {\partial q_k}. \]

Po zsumowaniu wszystkich składników otrzymujemy

\[ [\delta] = [K_e] + [K_g] + [K_c] + [K_C]. \tag{II.141} \label{II.141} \]

gdzie:
$K_e]$ – macierz sztywności sprężystej,
$ [K_g]$ – macierz sztywności geometrycznej,
$ [K_c]$ – macierz sztywności sprężystego podłoża,
$ [K_C]$ – macierz sztywności elementów dyskretnych (podpór sprężystych, zamocowań podatnych, sprężyn translacyjnych i obrotowych).

Układ równań Lagranga-Ritza przyjmuje zatem postać

\[ ([K_e]+[K_g]+[K_c]+[K_C]) \{\eta\} = \{\Delta_F\}. \tag{II.142} \label{II.142} \]

W przypadku zagadnień stateczności, gdy siłę osiową zapisuje się w postaci

\[ N=\Lambda N_0, \]

otrzymuje się

\[ [K_g] = -\Lambda[K_{g0}], \tag{II.143} \label{II.143} \]

a równanie równowagi przyjmuje postać

\[ ([K_e]+[K_c]+[K_C]-\Lambda[K_{g0}]) \{\eta\} = \{\Delta_F\}. \tag{II.144} \label{II.144} \]

Dla klasycznego zagadnienia wyboczenia pomija się obciążenia poprzeczne, czyli

\[ \{\Delta_F\} = \{0\}, \]

co prowadzi do jednorodnego układu równań

\[ ([K_e]+[K_c]+[K_C]-\Lambda[K_{g0}]) \{\eta\} = \{0\}. \tag{II.145} \label{II.145} \]

Warunek istnienia rozwiązania niezerowego ma postać

\[\det\left([K_e]+[K_c]+[K_C]-\Lambda[K_{g0}]\right)=0. \tag{II.146} \label{II.146} \]

Jest to równanie własne stateczności konstrukcji. Wartości własne $\Lambda$ wyznaczają mnożniki obciążenia krytycznego, natomiast odpowiadające im wektory \(\{\eta\}\) określają postacie wyboczenia wraz z odpowiadającymi im dodatkowymi stopniami swobody elementów dyskretnych.

Przedstawiona interpretacja pokazuje, że każdy składnik funkcjonału całkowitej energii potencjalnej generuje odpowiadającą mu macierz współczynników. Macierz całkowita jest sumą wkładów wynikających z energii odkształcenia elementów ciągłych, energii geometrycznej, energii sprężystego podłoża oraz energii elementów dyskretnych. W konsekwencji klasyczna metoda Ritza, jej rozszerzona postać oraz metoda elementów skończonych posiadają wspólne podstawy energetyczne, różniąc się przede wszystkim sposobem doboru funkcji aproksymujących oraz interpretacją zmiennych wariacyjnych.

Metoda elementów skończonych (MES) jako lokalna implementacja metody Ritza

Przedstawiona w poprzednim rozdziale interpretacja energetyczna metody Ritza prowadzi w sposób naturalny do metody elementów skończonych (MES). Obie metody opierają się na tym samym funkcjonale całkowitej energii potencjalnej oraz identycznym warunku jego stacjonarności. Różnica pomiędzy nimi nie dotyczy podstaw teoretycznych, lecz sposobu aproksymacji pola przemieszczeń.

W metodzie Ritza pole przemieszczeń opisuje się za pomocą funkcji globalnych zdefiniowanych na całym obszarze konstrukcji, natomiast w metodzie elementów skończonych konstrukcję dzieli się na skończoną liczbę elementów, na których stosowane są funkcje lokalne. W konsekwencji współczynniki wariacyjne $\{\eta\}$ zostają zastąpione przez wektory stopni swobody poszczególnych elementów, a funkcje aproksymujące przez funkcje kształtu elementów skończonych.

Ponieważ funkcjonał całkowitej energii potencjalnej pozostaje niezmieniony, wszystkie zależności wyprowadzone wcześniej dla metody Ritza zachowują swoją ważność również w metodzie elementów skończonych. W szczególności interpretacja macierzy współczynników przedstawiona równaniami (\ref{II.134})–(\ref{II.141}) prowadzi bezpośrednio do macierzy sztywności elementu skończonego. Poszczególne składniki energii odkształcenia, energii geometrycznej, energii sprężystego podłoża oraz energii elementów dyskretnych generują odpowiadające im składowe macierzy sztywności elementu, analogicznie jak miało to miejsce w metodzie Ritza.

Równanie równowagi pojedynczego elementu ma zatem identyczną strukturę jak równanie (\ref{II.142}). Różnica polega jedynie na tym, że odnosi się ono do lokalnych stopni swobody danego elementu. Następnie wszystkie macierze elementarne są składane zgodnie z połączeniami węzłów, tworząc globalną macierz sztywności całej konstrukcji. Proces ten, określany jako montaż macierzy, stanowi podstawową operację metody elementów skończonych i nie zmienia energetycznej interpretacji układu równań.

Analogicznie zagadnienia stateczności prowadzą do globalnego równania własnego, którego postać jest zgodna z równaniami (\ref{II.143})–(\ref{II.146}). Jedyną różnicą jest fakt, że poszczególne macierze powstają jako suma wkładów wszystkich elementów skończonych. Wartości własne wyznaczają mnożniki obciążenia krytycznego, natomiast odpowiadające im wektory własne opisują globalne postacie utraty stateczności konstrukcji.

Przedstawione porównanie prowadzi do istotnego wniosku, że metoda elementów skończonych nie stanowi odrębnej teorii mechaniki konstrukcji. Jest ona numeryczną implementacją metody Ritza, wykorzystującą lokalne funkcje aproksymujące zamiast funkcji globalnych. Obie metody opierają się na tej samej zasadzie minimum całkowitej energii potencjalnej, prowadzą do identycznych równań równowagi oraz posiadają jednakową interpretację energetyczną macierzy sztywności.

Ścisłe funkcje kształtu belki-słupa wyprowadzone wcześniej według Livesleya mogą być interpretowane jako dokładne funkcje kształtu elementu skończonego. Oznacza to, że uzyskane wcześniej ścisłe macierze sztywności stanowią jednocześnie dokładne macierze elementarne MES dla analizowanego elementu belkowo-słupowego, tworząc bezpośredni pomost pomiędzy klasyczną metodą Ritza a współczesną metodą elementów skończonych.

Agregacja  globalnej macierzy sztywności

Po wyznaczeniu macierzy sztywności poszczególnych elementów skończonych następuje ich połączenie w jeden globalny układ równań opisujący zachowanie całej konstrukcji. Operacja ta nosi nazwę agregacji macierzy sztywności i stanowi podstawowy etap metody elementów skończonych.  Każdy element skończony posiada własny lokalny układ stopni swobody oraz odpowiadającą mu macierz sztywności

\[ [K]^{(e)}, \]

wyznaczoną zgodnie z zasadą minimum całkowitej energii potencjalnej przedstawioną w poprzednim rozdziale. Macierz ta opisuje jedynie zachowanie pojedynczego elementu i nie uwzględnia jego połączenia z pozostałymi elementami konstrukcji. Po podziale konstrukcji na elementy skończone wszystkie węzły otrzymują numerację globalną. Dzięki temu każdy lokalny stopień swobody może zostać jednoznacznie przyporządkowany odpowiedniemu stopniowi swobody całej konstrukcji. Proces montażu polega na przeniesieniu współczynników macierzy elementarnych do odpowiednich miejsc globalnej macierzy sztywności oraz dodaniu współczynników odnoszących się do tych samych stopni swobody.

Jeżeli dwa lub więcej elementów posiada wspólny węzeł, odpowiadające mu współczynniki sztywności nie są zastępowane, lecz sumowane. Wynika to bezpośrednio z addytywności całkowitej energii potencjalnej układu. Energia całej konstrukcji jest bowiem sumą energii wszystkich elementów

\[ \Pi=\sum_{e=1}^{n_e}\Pi^{(e)}. \tag{II.147}\label{II.147} \]

Po dwukrotnym różniczkowaniu funkcjonału względem globalnych stopni swobody otrzymuje się globalną macierz sztywności

\[ [K]=\sum_{e=1}^{n_e}[K]^{(e)}. \tag{II.148}\label{II.148} \]

Zapis (\ref{II.148}) ma charakter symboliczny. W rzeczywistości sumowaniu podlegają jedynie te współczynniki macierzy elementarnych, które odpowiadają jednakowym globalnym stopniom swobody. Operację tę realizuje się z wykorzystaniem numeracji węzłów lub macierzy lokalizacji elementów.

Przykładowo, jeżeli dwa sąsiednie elementy posiadają wspólny węzeł, to współczynniki odpowiadające temu węzłowi zostają dodane do tych samych pozycji globalnej macierzy sztywności. W rezultacie sztywność wspólnego węzła stanowi sumę wkładów wszystkich elementów dochodzących do tego węzła.

Otrzymana w ten sposób globalna macierz sztywności zachowuje wszystkie właściwości macierzy elementarnych. Jest ona macierzą symetryczną, a dla liniowych zagadnień sprężystych również dodatnio określoną po uwzględnieniu warunków brzegowych. Ponadto zachowuje interpretację energetyczną przedstawioną wcześniej dla metody Ritza, gdyż każdy jej współczynnik odpowiada drugiej pochodnej całkowitej energii potencjalnej względem odpowiednich globalnych stopni swobody.

Po zakończeniu procesu agregacji oraz uwzględnieniu warunków brzegowych otrzymuje się globalny układ równań równowagi konstrukcji

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

który posiada identyczną strukturę jak równanie (\ref{II.142}) wyprowadzone wcześniej metodą Ritza. Oznacza to, że różnica pomiędzy obiema metodami nie dotyczy postaci równań równowagi, lecz jedynie sposobu ich konstruowania. W metodzie Ritza globalna macierz sztywności powstaje bezpośrednio z funkcji aproksymujących zdefiniowanych na całej konstrukcji, natomiast w metodzie elementów skończonych jest ona budowana poprzez agregację macierzy elementarnych odpowiadających poszczególnym fragmentom konstrukcji.

Przedstawiony proces posiada również prostą interpretację fizyczną. Każdy element wnosi do układu własną sztywność wynikającą z magazynowanej energii odkształcenia. W miejscach połączeń elementów wkłady te sumują się, dzięki czemu globalna macierz sztywności odzwierciedla jednocześnie lokalne właściwości wszystkich elementów oraz wzajemne oddziaływanie pomiędzy nimi. Z matematycznego punktu widzenia montaż macierzy jest więc bezpośrednią konsekwencją addytywności funkcjonału całkowitej energii potencjalnej, natomiast z fizycznego punktu widzenia stanowi zapis współpracy wszystkich elementów konstrukcji w przenoszeniu obciążeń.

Współczesne elementy węzłowe

Klasyczna metoda elementów skończonych zakładała, że połączenia pomiędzy elementami konstrukcji mają charakter idealny. Węzły traktowano jako połączenia całkowicie sztywne lub idealnie przegubowe, natomiast ich zadaniem było jedynie zapewnienie zgodności przemieszczeń oraz przekazanie sił pomiędzy sąsiednimi elementami.
Rzeczywiste właściwości mechaniczne połączeń elementów prętowych są bardziej złożóne i  w istotny sposób wpływają  na pracę całej konstrukcji. W konsekwencji współczesne metody obliczeniowe nie modeluja już idealnych węzłów na rzecz modeli uwzględniających ich rzeczywistą podatność.

Podstawową ideą współczesnego modelowania jest zastąpienie rzeczywistego połączenia układem elementów węzłowych opisujących poszczególne mechanizmy odkształcenia. Każdy z takich elementów reprezentuje określoną właściwość mechaniczną połączenia i posiada własną charakterystykę sztywności, która może być liniowa lub nieliniowa. Do najczęściej stosowanych elementów węzłowych należą elemenety węzłowe : podatności obrotowej,  translacyjnej, osiowej, skrętnej; mimośrodowego połączenia (offset); sztywnego ramienia (Rigid Link); jednostronnego kontaktu; przenoszący wyłącznie ściskanie, rozciąganie; element kontaktowy pomiędzy częściami konstrukcji; elementy współpracujące z  prętami Timoschenko lub Własowa.

Elementy te nie opisują zachowania całego pręta ani płyty lub tarczy, lecz lokalne właściwości mechaniczne połączenia. W zależności od stopnia szczegółowości modelu mogą one reprezentować pojedynczy mechanizm odkształcenia lub grupę współpracujących mechanizmów. Szczególne znaczenie posiada element podatności obrotowej, który umożliwia modelowanie połączeń półsztywnych. i przgubów plastycznych lub uogólnionych Zastępuje on klasyczne założenie idealnego przegubu lub sztywnego zamocowania zależnością pomiędzy momentem zginającym i względnym obrotem końców łączonych elementów. Pozwala to znacznie dokładniej odwzorować rzeczywistą pracę połączeń stalowych i żelbetowych.

Coraz większe znaczenie mają również elementy opisujące lokalną podatność translacyjną. Znajdują one zastosowanie przy modelowaniu podatnych podpór, połączeń z przekładkami elastomerowymi, łożysk mostowych oraz lokalnych deformacji elementów konstrukcyjnych. W nowoczesnych programach do analizy połączeń stalowych połączenie nie jest już reprezentowane pojedynczym przegubem ani jedną sprężyną. Zastępuje się je układem współpracujących elementów podatnych odpowiadających poszczególnym składnikom konstrukcyjnym. Każdy z tych elementów opisuje odrębny mechanizm odkształcenia i posiada własną charakterystykę sztywności.

Do najczęściej modelowanych składników połączeń należą: blachy czołowe, blachy węzłowe, śruby, spoiny, żebra usztywniające, środniki i półki kształtowników, strefy docisku pomiędzy elementami. Elementy te nie opisują zachowania całego pręta, płyty ani tarczy, lecz lokalne właściwości mechaniczne samego połączenia. W zależności od przyjętego poziomu szczegółowości modelu mogą reprezentować pojedynczy mechanizm odkształcenia lub zespół współpracujących mechanizmów zachodzących w obrębie węzła. Najważniejszym z nich jest element podatności obrotowej, umożliwiający modelowanie połączeń półsztywnych. W modelach nieliniowych może on również reprezentować przeguby plastyczne oraz uogólnione przeguby nieliniowe. Zastępuje on klasyczne założenie idealnego przegubu lub idealnego utwierdzenia zależnością pomiędzy momentem zginającym a względnym obrotem końców łączonych elementów. Pozwala to znacznie dokładniej odwzorować rzeczywistą pracę połączeń stalowych, żelbetowych i zespolonych, a także prawidłowo określić redystrybucję sił wewnętrznych w konstrukcji. Coraz większe znaczenie mają również elementy opisujące lokalną podatność translacyjną. Stosowane są one do modelowania podatnych podpór, połączeń z przekładkami elastomerowymi, łożysk mostowych, styków prefabrykatów, lokalnych deformacji fundamentów oraz innych połączeń, w których występują przemieszczenia względne pomiędzy łączonymi elementami.

W bardziej zaawansowanych modelach wykorzystywane są również elementy kontaktowe. Pozwalają one opisać zjawiska jednostronnego kontaktu, odrywania, zamykania szczelin, tarcia oraz selektywnego przenoszenia obciążeń, np. wyłącznie ściskania lub wyłącznie rozciągania. Elementy tego typu są obecnie powszechnie stosowane przy modelowaniu połączeń śrubowych, styków prefabrykatów, fundamentów oraz kontaktu pomiędzy współpracującymi częściami konstrukcji.

W nowoczesnych programach do analizy połączeń konstrukcyjnych pojedynczy przegub lub pojedyncza sprężyna zostały zastąpione układem wielu współpracujących elementów podatnych. Każdy z nich odpowiada określonemu składnikowi połączenia oraz opisuje jeden mechanizm jego odkształcenia. Podejście to umożliwia wierne odwzorowanie rzeczywistej pracy węzła bez konieczności modelowania jego pełnej geometrii za pomocą elementów bryłowych

W konstrukcjach żelbetowych podobne podejście stosowane jest do modelowania podatności połączeń monolitycznych i prefabrykowanych. Poszczególne elementy węzłowe mogą reprezentować podatność obrotową połączenia, podatność zakotwienia zbrojenia, podatność styków prefabrykatów, podatność przekładek elastomerowych lub lokalne odkształcenia stref przypodporowych. Dzięki temu ten sam formalizm elementów węzłowych może być stosowany niezależnie od materiału konstrukcyjnego.

Każdy z wymienionych składników może być traktowany jako lokalny element podatny o własnej charakterystyce sztywności. W rezultacie cały węzeł konstrukcyjny staje się układem współpracujących elementów mechanicznych, których wspólne działanie decyduje o sztywności, nośności i trwałości połączenia.

Takie podejście stanowi naturalne rozwinięcie klasycznej metody elementów skończonych. Modelowane są już nie tylko elementy ciągłe konstrukcji, lecz również lokalne mechanizmy odkształcenia występujące w samym połączeniu. Pozwala to analizować rzeczywistą pracę węzłów z dokładnością niemożliwą do uzyskania przy wykorzystaniu modeli opartych wyłącznie na idealnych przegubach i utwierdzeniach.

Z punktu widzenia rozszerzonej metody Ritza wszystkie elementy węzłowe podlegają identycznej interpretacji energetycznej. Każdy z nich wnosi do funkcjonału całkowitej energii potencjalnej własny składnik energii sprężystej odpowiadający określonemu mechanizmowi odkształcenia. Po wykonaniu operacji wariacyjnych składniki te generują odpowiednie przyrosty globalnej macierzy współczynników. Oznacza to, że elementy ciągłe, elementy dyskretne oraz elementy węzłowe różnią się wyłącznie sposobem opisu energii, natomiast podlegają temu samemu formalizmowi matematycznemu. Przyjęcie takiego ujęcia otwiera możliwość systematycznego rozszerzania modeli obliczeniowych poprzez wprowadzanie kolejnych elementów węzłowych opisujących nowe mechanizmy fizyczne. W dalszej części pracy koncepcja ta zostanie wykorzystana do modelowania podatności ścinania zgodnie z teorią Timoshenki oraz podatności deplanacji i transformacji bimomentów w teorii cienkościennych prętów Własowa.
Każdy z wymienionych składników może być traktowany jako lokalny element podatny o własnej sztywności. W rezultacie cały węzeł staje się układem współpracujących elementów mechanicznych, których wspólne działanie decyduje o sztywności oraz nośności połączenia. Takie podejście stanowi rozwinięcie klasycznej metody elementów skończonych. Modelowane są już nie tylko elementy ciągłe konstrukcji, lecz również lokalne mechanizmy odkształcenia występujące w samym połączeniu. Współczesne programy obliczeniowe umożliwiają dzięki temu analizę rzeczywistej pracy węzłów z uwzględnieniem podatności ich poszczególnych składników. 
Z punktu widzenia rozszerzonej metody Ritza wszystkie elementy węzłowe mogą być interpretowane w jednakowy sposób. Każdy z nich wnosi do funkcjonału całkowitej energii potencjalnej własny składnik energii sprężystej odpowiadający określonemu mechanizmowi odkształcenia. Po wykonaniu operacji wariacyjnych składniki te generują odpowiednie przyrosty globalnej macierzy współczynników. Oznacza to, że elementy ciągłe, elementy dyskretne oraz elementy węzłowe podlegają temu samemu formalizmowi matematycznemu i różnią się jedynie postacią funkcjonału opisującego ich zachowanie.

Takie ujęcie otwiera możliwość dalszego rozszerzania modeli obliczeniowych poprzez wprowadzanie nowych elementów węzłowych opisujących kolejne mechanizmy fizyczne. W następnych rozdziałach zostanie ono wykorzystane do modelowania podatności ścinania zgodnie z teorią Timoshenki oraz podatności deplanacji i transformacji bimomentów w teorii cienkościennych prętów Własowa.

Dyskretny element węzłowy z podatnością obrotową

Dyskretny element węzłowy z podatnością obrotową stanowi najprostszy przykład elementu węzłowego stosowanego we współczesnej metodzie elementów skończonych. Opisuje on lokalną podatność połączenia pomiędzy dwoma elementami prętowymi, zastępując klasyczne założenie idealnego przegubu lub całkowicie sztywnego zamocowania. W praktyce inżynierskiej wiele połączeń wykazuje zachowanie nieliniowe wynikające z uplastycznienia materiału, poślizgu elementów łączących, powstawania szczelin, degradacji sztywności lub innych zjawisk fizycznych. W takich przypadkach liniowa zależność konstytutywna zostaje zastąpiona ogólną funkcją opisującą charakterystykę moment–obrót połączenia.

Rozważmy połączenie dwóch prętów, którego jedynymi stopniami swobody są obroty końców elementów. Zakłada się, że zależność pomiędzy względnym obrotem końców elementu \(\Delta\varphi\) a momentem zginającym \(M\) przenoszonym przez połączenie opisana jest ogólną funkcją konstytutywną

\[M=f(\Delta\varphi), \tag{II.150} \label{II.150} \]

gdzie
\[ \Delta\varphi=\varphi_2-\varphi_1. \tag{II.151} \label{II.151} \]

oznacza względny obrót końców elementu.

Energia odkształcenia zgromadzona w elemencie wynika z pracy wykonanej podczas wzrostu względnego obrotu i może zostać zapisana w postaci

\[ \Pi= \int_{0}^{\Delta\varphi} M(\theta)\,d\theta, \tag{II.152} \label{II.152} \]

gdzie zmienna \(\theta\) oznacza bieżącą wartość względnego obrotu podczas procesu obciążania.

Siły uogólnione elementu wyznacza się z pierwszej pochodnej funkcjonału energii

\[ \mathbf{Q} = \cfrac{\partial\Pi} {\partial\mathbf{q}}, \tag{II.153} \label{II.153} \]

gdzie
\[ \mathbf{q} = \begin{Bmatrix} \varphi_1\\ \varphi_2 \end{Bmatrix}. \tag{II.154} \label{II.154} \]

Po wykonaniu różniczkowania otrzymuje się

\[ \mathbf{Q} = M(\Delta\varphi) \begin{Bmatrix} -1\\ 1 \end{Bmatrix}. \tag{II.155} \label{II.155} \]

W analizie nieliniowej szczególne znaczenie posiada macierz sztywności stycznej elementu, wyznaczana z drugiej pochodnej funkcjonału energii lub równoważnie z pochodnej wektora sił uogólnionych względem stopni swobody

\[ \mathbf{K}_t = \cfrac{\partial\mathbf{Q}} {\partial\mathbf{q}}. \tag{II.156} \label{II.156} \] 

Po uwzględnieniu zależności

\[ C_t= \cfrac{dM} {d(\Delta\varphi)}, \tag{II.157} \label{II.157} \]

otrzymuje się

\[ \mathbf{K}_t= C_t \begin{bmatrix} 1&-1\\ -1&1 \end{bmatrix} \tag{II.158} \label{II.158} \]

gdzie \(C_t\) oznacza chwilową (styczną) sztywność obrotową połączenia.  Element ten posiada jedynie dwa stopnie swobody odpowiadające obrotom końców połączenia i nie jest związany z żadną długością geometryczną. Opisuje wyłącznie lokalną podatność węzła, dlatego zaliczany jest do grupy dyskretnych elementów węzłowych. Po zmontowaniu z globalną macierzą sztywności wnosi jedynie własny wkład do odpowiednich wierszy i kolumn macierzy globalnej zgodnie z klasyczną procedurą agregacji stosowaną w metodzie elementów skończonych. Dzięki temu może zostać umieszczony pomiędzy dowolnymi elementami belkowymi, ramowymi lub płytowymi bez jakiejkolwiek modyfikacji podstawowego algorytmu metody elementów skończonych.

Liniowy element podatności obrotowej stanowi szczególny przypadek przedstawionego modelu w którym 

\[ C_t=C_{\varphi}=\mathrm{const}, \tag{II.159} \label{II.159} \]

Odpowiedni dobór funkcji \(M=f(\Delta\varphi)\) umożliwia modelowanie różnych typów elementów węzłowych, między innymi połączeń półsztywnych o charakterystyce nieliniowej, przegubów plastycznych, przegubów z umocnieniem materiału, przegubów z degradacją sztywności oraz połączeń wyznaczonych bezpośrednio na podstawie badań doświadczalnych.
W kolejnych przykładach przedstawione zostaną najczęściej stosowane modele funkcji \(M=f(\Delta\varphi)\), poczynając od klasycznego modelu sprężysto-plastycznego, poprzez modele z umocnieniem materiału, aż do charakterystyk nieliniowych wykorzystywanych do opisu rzeczywistych połączeń konstrukcyjnych.

Węzeł – przegub nośnościowy 

Szczególnym przypadkiem przedstawionego modelu jest element węzłowy o stałej nośności momentowej</b>, zwany dalej przegubem nośnościowym. Po osiągnięciu granicznego momentu nośności \(M_R\) element zachowuje zdolność przenoszenia stałego momentu, natomiast jego sztywność styczna zanika. Model ten nie opisuje wyłącznie klasycznego przegubu plastycznego, lecz dowolny mechanizm osiągnięcia stanu granicznego połączenia, niezależnie od materiału konstrukcyjnego.

Zależność konstytutywna elementu przyjmuje postać

\[ M=M_R=\mathrm{const}, \tag{II.160} \label{II.160} \]

Energia odkształcenia elementu wynosi

\[ \Pi = \int_0^{\Delta\varphi} M_R\,d\theta = M_R\,\Delta\varphi. \tag{II.161} \label{II.161} \]

Wektor sił węzłowych elementu jest równy

\[ \mathbf{Q} = M_R \begin{Bmatrix} -1\\ 1 \end{Bmatrix}. \tag{II.162} \label{II.162} \]

Ponieważ przenoszony moment nie zależy od dalszego wzrostu względnego obrotu, chwilowa sztywność styczna elementu wynosi

\[ C_t = \cfrac{dM}{d(\Delta\varphi)} = 0. \tag{II.163} \label{II.163} \]

Po podstawieniu zależności (II.124) do równania (II.119) otrzymuje się macierz sztywności stycznej przegubu nośnościowego

\[  \mathbf{K}_t = \mathbf{0}.  \tag{II..164} \label{II..164} \]

Oznacza to, że po osiągnięciu granicznej nośności momentowej element nadal przenosi stały moment \(M_R\), lecz nie zwiększa już swojej sztywności. Dalszy wzrost względnego obrotu powoduje jedynie rozwój odkształceń związanych z osiągniętym mechanizmem granicznym, natomiast element przestaje wnosić wkład do globalnej macierzy sztywności konstrukcji.

Przedstawiony model ma charakter uniwersalny i obejmuje znacznie szerszą klasę zjawisk niż klasyczny przegub plastyczny. Może on opisywać między innymi:

• uplastycznienie przekroju stalowego,
• powstanie przegubu plastycznego w konstrukcjach żelbetowych,
• osiągnięcie nośności połączenia śrubowego lub spawanego,
• zniszczenie połączeń drewnianych,
• osiągnięcie nośności połączeń prefabrykowanych,
• dowolny lokalny mechanizm utraty nośności połączenia.

Z tego względu określenie   przegub nośnościowy  wydaje się bardziej właściwe niż tradycyjnie stosowane pojęcie  przegubu plastycznego,  które odnosi się wyłącznie do jednego z możliwych mechanizmów osiągnięcia stanu granicznego. W praktyce odpowiedni dobór funkcji \(M=f(\Delta\varphi)\) umożliwia płynne przejście od modelu liniowo sprężystego, poprzez modele nieliniowe, aż do modelu  przegubu nośnościowego oraz modeli uwzględniających degradację lub umocnienie materiału.

Uogólnienie teorii belki na elementy krzywoliniowe

Funkcjonał (\ref{II.111}) stanowi kompletne energetyczne sformułowanie belki Bernoulli-Eulera uwzględniające zginanie, ściskanie osiowe, efekt drugiego rzędu P–Δ, podłoże sprężyste oraz podpory sprężyste.
Wyprowadzenie to opiera się jednak na założeniu, że oś elementu w stanie początkowym jest linią prostą. W praktyce inżynierskiej założenie takie stanowi jedynie idealizację, ponieważ rzeczywiste elementy konstrukcyjne zawsze wykazują pewne odchylenia od geometrii nominalnej.

Macierz podatności elementu z imperfekcją geometryczną

Ogólna procedura wyznaczania macierzy podatności elementów krzywoliniowych została przedstawiona przez Livesleya [10]. W najogólniejszym przypadku geometria elementu opisywana jest funkcjami parametrycznymi osi środkowej, natomiast współczynniki macierzy podatności wyznaczane są na podstawie energii odkształcenia oraz twierdzenia Castigliano.

W teorii imperfekcyjnej pełne sformułowanie ogólne nie jest jednak konieczne. Znacznie bardziej interesujące jest prześledzenie tej samej procedury dla geometrii odpowiadającej początkowemu wygięciu pręta. Pozwala to bezpośrednio ocenić wpływ imperfekcji geometrycznej na podatność, sztywność oraz stateczność elementu.

W dalszych rozważaniach ograniczymy się do imperfekcji opisanej funkcją sinusoidalną zgodną z pierwszą postacią wyboczeniową pręta przegubowo podpartego. Wybór taki nie oznacza, że rzeczywiste imperfekcje mają przebieg sinusoidalny. Jest to model obliczeniowy odpowiadający składowej imperfekcji wywierającej największy wpływ na utratę stateczności.
Rzeczywiste imperfekcje geometryczne powstają w procesie produkcji, transportu, montażu oraz eksploatacji konstrukcji. Mają one charakter losowy i mogą przyjmować praktycznie dowolny przebieg geometryczny. W ogólnym przypadku imperfekcja rzeczywista stanowi superpozycję wielu składowych o różnych amplitudach i długościach fal.

W analizie stateczności największe znaczenie posiada jednak składowa zgodna z pierwszą postacią wyboczeniową. To właśnie ona ulega najsilniejszemu wzmocnieniu pod wpływem ściskania i odpowiada kierunkowi najszybszego obniżania obciążenia krytycznego. Z tego względu w teorii imperfekcyjnej oraz w normach projektowych stosuje się imperfekcje zastępcze, których kształt nawiązuje do odpowiednich postaci wyboczeniowych konstrukcji.
Podejście takie stanowi świadome uproszczenie prowadzące do modeli konserwatywnych. Nie ma ono na celu wiernego odwzorowania rzeczywistej geometrii elementu, lecz uchwycenie dominującego mechanizmu odpowiedzialnego za utratę stateczności. W tym sensie przyjmowane w obliczeniach kształty imperfekcji są związane z postaciami własnymi układu.

Pełne wyprowadzenie  macierzy podatności krzywoliniowego elementu belki Bernoulliego z impefekcją w kształcie sinusoidy podano w  Dodatku II.A. 

Przykłady 

Przykład II.P1 [ Nieliniowa geometrycznie belka wolnodparta. Rozwiązanie ścisłe] 

Przeprowadzić analizę nieliniową geometrycznie belki wolnopodpartej ze sprężystą podporą pionową w węźle (2) (rys. II.1b), ściskanej osiowo , dla której warunki podporowe postać

\[ u_1=0, \qquad w_1=0, \qquad V_2=-C_{\Delta}w_2, \tag{II.P1.1}\]

gdzie \(C_{\Delta}\) oznacza sztywność pionowej podpory sprężystej w węźle (2).

Warunek $ u_1=0$  eliminuje przemieszczenie podłużne w podporze lewej, natomiast warunek $ w_1=0$ eliminuje przemieszczenie pionowe. Zależność $ V_2=-C_{\Delta}w_2$ opisuje reakcję podpory sprężystej zgodnie z prawem Hooke’a.  Wprowadza się bezwymiarowy parametr sztywności podpory

\[ \bar C_{\Delta}=\cfrac{C_{\Delta}L^3}{EI}. \tag{II.P1.2} \label{II.P1.2}\]

Dla \(\bar C_{\Delta}=0\) otrzymujemy klasyczną belkę wolnopodpartą, natomiast dla \(\bar C_{\Delta}\rightarrow\infty\) podpora zachowuje się jak podpora nieprzemieszczalna w kierunku pionowym. 

Zmodyfikowana macierz sztywności

Ponieważ analizowany system konstrukcyjny składa się z pojedynczego elementu belkowego, globalna macierz sztywności układu jest identyczna z macierzą sztywności elementu

\[ [K]=[k]^{(e)}, \tag{II.P1.3}\]

gdzie \([k]^{(e)}\) określono równaniem (\ref{II.22}).

Uwzględnienie warunków brzegowych można przeprowadzić poprzez  eliminację stopni swobody $u_1$ oraz $w_1$ ospowiadajacych podporom stałym. Warunek sprężystego podparcia prowadzi do zwiększenia współczynnika sztywności odpowiadającego przemieszczeniu $w_2$ o dodatkowy składnik wynikający ze sztywności podpory. W konsekwencji wektor niewiadomych przemieszczeń przyjmuje postać

\[ \{q\} = \begin{Bmatrix} \varphi_1 & u_2 & w_2 & \varphi_2 \end{Bmatrix}^{T}, \tag{II.P1.4}\]

W konsekwencji zmodyfikowana macierz sztywności układu wynosi</p>

\[ [\tilde K]=\cfrac{EI}{L^3} \begin{bmatrix}
4L^2\Psi_3 & 0 & -6L\Psi_2 & 2L^2\Psi_4 \\[3mm]
0 & \dfrac{EA\,L^2}{EI} & 0 & 0 \\[3mm]
& \mathrm{SYM} & 12\Psi_1+\bar C_{\Delta} & -6L\Psi_2 \\[3mm]
& & & 4L^2\Psi_3 \end{bmatrix}. \tag{II.P1.5} \label{II.P1.5}\]

gdzie:
$\bar C_{\Delta}$ jest bezwymiarową sztywnością podpory sprężystej zdefiniowaną równaniem (\ref{II.P1.2}), 
$\Psi_i , \quad (i=1,\ldots,4)$  – funkcje statecznościowe określone zależnościami  (\ref{II.23}) – (\ref{II.24})

Zmodyfikowana macierz podatności

W celu wyznaczenia przemieszczeń przywęzłowych od zadanych obciążeń odwraca się zmodyfikowaną macierz sztywności $[\widetilde K]$. Po symbolicznym wykonaniu operacji odwrócenia oraz wyłączeniu wspólnego czynnika otrzymuje się zwartą postać zmodyfikowanej macierzy podatności

\[ [\widetilde K]^{-1}= \cfrac{L}{EI\,(2\Psi_3-\Psi_4)\,\Theta} \begin{bmatrix}
-\dfrac{EI(2\Psi_3-\Psi_4)\Theta}{EA} & 0 & 0 & 0\\[2mm]
0 & k_{22}L^2 & k_{23}L & k_{23}L\\[2mm]
0 & k_{23}L & k_{33} & k_{34}\\[2mm]
0 & k_{23}L & k_{34} & k_{33}
\end{bmatrix},
\tag{II.P1.6} \label{II.P1.6} \]

gdzie:
$\Theta= 36\Psi_2^2- (2\bar C_{\Delta}+24\Psi_1)\Psi_3- (\bar C_{\Delta}+12\Psi_1)\Psi_4, \tag{II.P1.7}$
$ k_{22}=\Psi_4^{\,2}-4\Psi_3^{\,2}, \tag{II.P1.8}$,
$ k_{23}=3\Psi_2(2\Psi_3-\Psi_4), \tag{II.P1.9}$,
$k_{33}=9\Psi_2^{\,2}-(\bar C_{\Delta}+12\Psi_1)\Psi_3, \tag{II.P1.10}$
$ k_{34}=\cfrac{1}{2}(\bar C_{\Delta}+12\Psi_1)\Psi_4-9\Psi_2^{\,2}. \tag{II.P1.11}$

Współczynniki $k_{22}$–$k_{34}$ zależą wyłącznie od funkcji statecznościowych $\Psi_i$ oraz bezwymiarowej sztywności podpory sprężystej $\bar C_{\Delta}$. Takie przedstawienie macierzy podatności umożliwia wyodrębnienie wspólnego mianownika, dzięki czemu wszystkie elementy macierzy przyjmują postać prostych wielomianów względem funkcji statecznościowych i parametrów układu, co znacznie upraszcza zarówno dalsze przekształcenia analityczne, jak i implementację komputerową.

Zmodyfikowane równoważniki węzłowe 

W przykładzie wykorzystano ścisły równoważny wektor obciążeń elementu (\ref{II.46}). Po uwzględnieniu warunków brzegowych oraz redukcji układu do aktywnych stopni swobody zmodyfikowanej macierzy podatności (\ref{II.P1.6}) otrzymuje się zredukowany wektor równoważników obciążenia

\[ \{\widetilde P_p\}= \{ 0 \, , \, \dfrac{Lp}{2} \, , \, Lp\left(\dfrac{L}{\mu}-\dfrac{L}{2}\cot\dfrac{\mu}{2}\right) \, , \, \dfrac{L^2p}{2\mu} \left(\mu\cot\dfrac{\mu}{2}-2\right) \}^T, \tag{II.P1.12} \label{II.P1.12} \]

zgodnym z kolejnością przyjętą w zmodyfikowanej macierzy sztywności (\ref{II.P1.6}). Zerowa pierwsza składowa wynika z braku osiowego obciążenia węzłowego odpowiadającego przemieszczeniu $u_2$. Przyjęto dodatni zwrot obciążenia poprzecznego zgodny z dodatnim kierunkiem osi lokalnej elementu; dla oznaczeń stosowanych na rysunku zachodzi zależność $p = -q$.

Strzałka ugięcia

Zgodnie z interpretacją przedstawioną w (\ref{II.61}) całkowita strzałka ugięcia belki stanowi sumę rozwiązania szczególnego oraz rozwiązania jednorodnego 

\[ f=f_q+f_h, \tag{II.P1.13} \label{II.P1.13} \]

gdzie rozwiązanie jednorodne $f_h$ wynika z przemieszczeń i obrotów przywęzłowych wyznaczonych z macierzy podatności (\ref{II.P1.6}) i nieliniowego wektora równoważników obciążenia (\ref{II.P1.12}).

Rozwiązanie szczególne odpowiada bezpośredniemu działaniu równomiernie rozłożonego obciążenia poprzecznego i ma posta

\[ w_q(\xi)= \cfrac{qL^4}{24EI}
\left( \xi-2\xi^3+\xi^4 \right). \tag{II.P1.14} \label{II.P1.14} \]

Po podstawieniu $\xi=0.5$ otrzymuje się klasyczną strzałkę ugięcia belki swobodnie podpartej

\[f_q= w_q(0.5)= \cfrac{5qL^4}{384EI} \equiv f_0. \tag{II.P1.15} \label{II.P1.15} \]

Rozwiązanie jednorodne $f_h$ wyznaczono symbolicznie z macierzy  podatności (\ref{II.P1.6}), oraz wektora obciążeń (\ref{II.P1.12}), traktując funkcje statecznościowe $\Psi_i$ jako niezależne współczynniki algebraiczne. Po wykonaniu przekształceń symbolicznych otrzymano

\[ f_h= \cfrac{24\,qL^4}{384EI} \left[ \cfrac{2\Psi_3+\Psi_4} {(\bar C_\Delta+12\Psi_1)(2\Psi_3+\Psi_4)-36\Psi_2^{\,2}} + \cfrac{\left(2-\mu\cot\cfrac{\mu}{2}\right)\tan\cfrac{\mu}{4}} {(2\Psi_3-\Psi_4)\mu^2} \right]. \tag{II.P1.16} \label{II.P1.16} \]

Po znormalizowaniu całkowitej strzałki ugięcia względem wartości odniesienia $f_0$ otrzymuje się ścisłą zależność

\[ \bar f = \cfrac{f}{f_0} =  \cfrac{96}{5} \left[ \cfrac{2\Psi_3+\Psi_4} {(\bar C_\Delta+12\Psi_1)(2\Psi_3+\Psi_4)-36\Psi_2^{\,2}} + \cfrac{\left(2-\mu\cot\cfrac{\mu}{2}\right)\tan\cfrac{\mu}{4}}
{(2\Psi_3-\Psi_4)\mu^2} \right]. \tag{II.P1.17} \label{II.P1.17} \]

gdzie $f_0$ określono zależnością (\ref{II.P1.15}), natomiast funkcje statecznościowe $\Psi_i$, parametr ściskania $\mu$ oraz bezwymiarową sztywność podpory $\bar C_\Delta$ zdefiniowano wcześniej.

Zależność (\ref{II.P1.17}) stanowi ścisłe rozwiązanie analizowanego zagadnienia i będzie w dalszej części pracy wykorzystywana jako rozwiązanie odniesienia przy ocenie dokładności metod przybliżonych, w szczególności metody Ritza.

Na rys. II.2  w tekście przedstawiono przebiegi znormalizowanej strzałki ugięcia $\bar f$ w funkcji bezwymiarowego współczynnika obciążenia osiowego $\bar\Lambda$ dla różnych wartości bezwymiarowej sztywności podpory sprężystej $\bar C_\Delta$. Wartość $\bar f=1$ odpowiada klasycznej strzałce ugięcia belki wywołanej wyłącznie obciążeniem poprzecznym. Wraz ze wzrostem siły osiowej rośnie udział rozwiązania jednorodnego, powodując nieliniowy wzrost całkowitego ugięcia. Zbliżanie się do obciążenia krytycznego prowadzi do gwałtownego wzrostu przemieszczeń, będącego konsekwencją sprzężenia efektów geometrycznej nieliniowości i podatności podpory sprężystej.
Współczynnik 

\[ \Theta= 36\Psi_2^{\,2}- (2\bar C_\Delta+24\Psi_1)\Psi_3- (\bar C_\Delta+12\Psi_1)\Psi_4 \tag{II.P1.18} \label{II.P1.18} \]

jest wyznacznikiem bezwymiarowej części zmodyfikowanej macierzy sztywności (\ref{II.P1.5}). Warunek

\[ \Theta=0  \tag{II.P1.19} \label{II.P1.19} \]

oznacza utratę odwracalności macierzy podatności (\ref{II.P1.7}) i wyznacza stan krytyczny analizowanego układu. Jest to jednocześnie równanie stateczności belki z podporą sprężystą, określające zależność pomiędzy bezwymiarowym obciążeniem osiowym $\bar\Lambda$ a bezwymiarową sztywnością podpory $\bar C_\Delta$.
Współczynnik $\Theta$ stanowi wyznacznik bezwymiarowej części zmodyfikowanej macierzy sztywności.

Asymptota nie wynika z nieograniczonego wzrostu przemieszczeń węzłowych, lecz z nieograniczonego wzrostu różnicy obrotów końcowych elementu. W miarę zbliżania się do obciążenia krytycznego część jednorodna rozwiązania zaczyna dominować nad rozwiązaniem szczególnym, powodując gwałtowny wzrost całkowitej strzałki ugięcia. Wzrost ten zachodzi jednak dopiero w bardzo wąskim otoczeniu punktu krytycznego. Dla przykładu, przy $a=1{,}569$ otrzymano już $\bar f\approx362$, podczas gdy odpowiadająca wartość współczynnika obciążenia wynosi $\bar\Lambda\approx0{,}9977$, a więc różni się od wartości krytycznej jedynie o około $0{,}23\%$.

Z punktu widzenia praktyki inżynierskiej asymptota ma przede wszystkim znaczenie teoretyczne. W rzeczywistych konstrukcjach stan graniczny wyznaczany jest zwykle przez uplastycznienie materiału, imperfekcje geometryczne, mimośrody obciążenia lub inne zjawiska występujące przed osiągnięciem idealnego obciążenia krytycznego Eulera. W bezpośrednim otoczeniu punktu krytycznego nawet niewielkie odchylenia parametrów materiałowych, geometrycznych lub podporowych powodują bardzo duże zmiany obliczonego ugięcia, prowadząc do gwałtownego wzrostu niepewności rozwiązania. W konsekwencji matematyczna asymptota wyznacza granicę stosowalności idealnego modelu liniowo-sprężystego, a nie rzeczywistą nośność konstrukcji.
Przedstawione rozwiązanie pozwala również na ilościową ocenę wpływu podatności podpory na zachowanie układu. Wzrost bezwymiarowej sztywności podpory $\bar C_\Delta$ powoduje zmniejszenie ugięć oraz przesunięcie punktu krytycznego w kierunku większych wartości współczynnika obciążenia $\bar\Lambda$. Oznacza to jednoczesny wzrost sztywności i nośności statecznościowej konstrukcji. Otrzymane wyniki wykazują, że nawet pojedyncza podpora sprężysta może istotnie ograniczyć efekty drugiego rzędu oraz znacząco zwiększyć odporność układu na utratę stateczności.

Przykład II.P2 [ Nieliniowy geometrycznie słup wspornikowy. Rozwiązanie ścisłe]

Przeprowadzić analizę geometrycznie nieliniową słupa wspornikowego przedstawionego na rys. II.1c. Podstawa słupa jest sprężyście utwierdzona za pomocą podpory obrotowej o sztywności $C^\varphi$. Rozpatruje się dwa przypadki obciążenia poprzecznego: siłę skupioną \(H\) przyłożoną w węźle górnym (1) oraz obciążenie równomiernie rozłożone $h$ na całej wysokości słupa.

Wprowadza się bezwymiarowy parametr sztywności obrotowej podpory

\[ \bar C_\varphi=\cfrac{C^\varphi L}{EI}. \tag{II.P2.1} \]

Dla $ \bar C_\varphi=0$ otrzymuje się podporę przegubową, natomiast granica $\bar C_\varphi\rightarrow\infty$ odpowiada klasycznemu utwierdzeniu sztywnemu. Parametr $\bar C_\varphi$ stanowi zatem miarę stopnia zamocowania obrotowego podstawy słupa.

Ponieważ analizowany układ konstrukcyjny składa się z pojedynczego elementu belki-słupa, globalna macierz sztywności jest identyczna z macierzą sztywności elementu

\[ [K]=[k]^{(e)}. \tag{II.P2.2} \]

gdzie macierz \([k]^{(e)}\) określono równaniem (\ref{II.22}).

Warunki brzegowe i zmodyfikowana macierz sztywności

Warunki podporowe dla słupa wspornikowego mają postać

\[ u_2 = 0,\qquad w_2 = 0, \qquad M_2= – C^\varphi\varphi_2. \tag{II.P2.3} \]

Warunki $u_2=0$ oraz $ w_2=0$ eliminują odpowiednio przemieszczenie podłużne i poprzeczne podstawy słupa. Natomiast zależność $ M_2=-C^\varphi\varphi_2$ opisuje reakcję sprężystej podpory obrotowej zgodnie z prawem Hooke’a. Uwzględnienie warunków podporowych prowadzi do eliminacji odpowiednich stopni swobody oraz modyfikacji równania równowagi węzła 2 przez dołączenie sztywności obrotowej podpory. Dla wektora niewiadomych przemieszczeń

\[ \{q\}= \begin{Bmatrix} u_1 & w_1 & \varphi_1 & \varphi_2 \end{Bmatrix}^{T}, \tag{II.P2.4} \]

Po uwzględnieniu warunków brzgwowych otrzymuje się zmodyfikowaną macierz sztywności układu</p>

\[ [\widetilde K] = \cfrac{EI}{L} \begin{bmatrix}
\dfrac{EA}{EI} & 0 & 0 & 0 \\[2mm]
0 & \dfrac{12}{L^2}\Psi_1 & \dfrac{6}{L}\Psi_2 & \dfrac{6}{L}\Psi_2 \\[2mm]
0 & \dfrac{6}{L}\Psi_2 & 4\Psi_3 & 2\Psi_4 \\[2mm]
0 & \dfrac{6}{L}\Psi_2 & 2\Psi_4 & 4\Psi_3+\bar C_\varphi \end{bmatrix}. \tag{II.P2.5} \]

Współczynnik $\bar C_\varphi$ występuje wyłącznie w elemencie odpowiadającym obrotowi podpory, co wynika z bezpośredniego uwzględnienia sprężystego zamocowania obrotowego w równaniu równowagi węzła 2.</p>

Odwrotność macierzy zmodyfikowanej

Po symbolicznym odwróceniu zmodyfikowanej macierzy sztywności oraz uporządkowaniu otrzymanych wyrażeń względem funkcji statecznościowych otrzymuje się zwartą postać macierzy podatności geometrycznie nieliniowego słupa wspornikowego:

\[ [\widetilde K]^{-1} = \cfrac{L}{EI\,\Theta} \begin{bmatrix}
\dfrac{EI\,\Theta}{EA} & 0 & 0 & 0 \\[3mm]
& k_{22} & k_{24}+\dfrac{L}{2}\Psi_2\bar C_\varphi & k_{24} \\[3mm]
& & 3\Psi_2^2-\Psi_1\left(4\Psi_3+\bar C_\varphi\right) & 2\Psi_1\Psi_4-3\Psi_2^2 \\[3mm]
& \mathrm{SYM} & & 3\Psi_2^2-4\Psi_1\Psi_3
\end{bmatrix}. \tag{II.P2.6} \]

gdzie:
$ k_{22} = -\cfrac{L^2}{3} \left( 4\Psi_3^2-\Psi_4^2+\bar C_\varphi\Psi_3 \right), \tag{II.P2.7} $
$ k_{24} = L\Psi_2(2\Psi_3-\Psi_4)\tag{II.P2.8} $

\[ \Theta= \bar C_\varphi \left( 3\Psi_2^{\,2}-4\Psi_1\Psi_3 \right) – 4\left(2\Psi_3-\Psi_4\right) \cdot \left [\Psi_1\left(2\Psi_3+\Psi_4\right)-3\Psi_2^{\,2} \right]. \tag{II.P2.9} \]

Współczynnik \(\Theta\) jest wspólnym mianownikiem wszystkich elementów macierzy podatności oraz jednocześnie wyznacznikiem bezwymiarowej części zmodyfikowanej macierzy sztywności.

Warunek $ \Theta=0 $ określa stan krytyczny słupa wspornikowego. W tym punkcie macierz sztywności układu staje się osobliwa, a rozwiązanie zagadnienia traci jednoznaczność, co odpowiada utracie stateczności konstrukcji

Równoważniki węzłowe (wektor obciążeń zastępczych)

Rozpatrzono trzy podstawowe przypadki obciążenia poprzecznegowspornika:

• $H$ siła skupiona przyłożona do górnego końca słupa (1),
• $h$ obciążenie równomiernie rozłożone , działające na całej jego wysokości (1)-(2).
• $M$ moment skupiony przyłożony do górnego końca słupa(1) ,, w dalszych przykladach interetowany jako $M=N \cdot e$ , czyli pochodzący od siły pinowej N na mmośrodzie $e$ 

Równoważniki węzłowe (wektr obciążeń, zastępczych), odpowiadający stopniom swobody ( II.P2.4)  wynoszą

 $H$:
$\{F\}_H= \begin{bmatrix} 0 \, ; \,  & H\, , \, & 0\, , \, & 0 \end{bmatrix}^{T}. \tag{II.P2.10}$
$M$:
$ \{F\}_M= \begin{bmatrix} 0\, ; \, & 0\, , \, & M\, , \, & 0 \end{bmatrix}^{T}.\tag{II.P2.11}$
$h$
dla równomiernie rozłożonego obciążenia poprzecznego \(h\) nieliniowy wektor obciążenia wyznacza się z zależności (II.36) w postaci

\[ \{F\}_h= \begin{bmatrix} 0 & \dfrac{hL}{2} &  -\dfrac{hL^{2}\left(\mu\cot\left(\dfrac{\mu}{2}\right)-2\right)}{2\mu^{2}} & \dfrac{hL^{2}\left(\mu\cot\left(\dfrac{\mu}{2}\right)-2\right)}{2\mu^{2}} \end{bmatrix}^{T}. \tag{II.P2.12}\]

Wychylenie wierzchołka wspornika

Znajomość macierzy podatności (II.P2.6) umożliwia wyznaczenie przemieszczeń przywęzłowych dla dowolnego przypadku obciążenia. Przemieszczenia wyznacza się z zależności

\[ \{q\}=[\widetilde K]^{-1}\{ F \}. \tag{II.P2.13} \label{II.P2.13}\]

Wyznaczone z zależności (II.P2.11) przemieszczenie poziome $w$ wierzchołka wspornika zostało znormalizowane względem klasycznego ugięcia wspornika wywołanego odpowiednim obciążeniem ($H$, $M$ lub $h$):

\[ \bar w=\cfrac{w_2}{w_{2,0}}. \tag{II.P2.14} \]

Dla poszczególnych przypadków obciążenia, przed podstawieniem funkcji statecznościowych , uzyskano następujące zależności na znormalizowane wychylenie pozziome wierzchołka wspornika$w$ od poszczególnych wynuszeń: 

• od siła poziomej $H$

\[ \bar{w}_H=  \cfrac{\Psi_{3}\left(\bar{C}_{\varphi}+4\Psi_{3}\right)-\Psi_{4}^{\,2}} {D_H}. \tag{II.P2.15}\label{II.P2.15}\]

\[D_H= \bar{C}_{\varphi}\left(4\Psi_{1}\Psi_{3}-3\Psi_{2}^{\,2}\right) + 4\left(2\Psi_{3}-\Psi_{4}\right)\cdot \left[ \Psi_{1}\left(2\Psi_{3}+\Psi_{4}\right)-3\Psi_{2}^{\,2} \right].\]

Wielkość jest znormalizowana względem klasycznego ugięcia

\[ w_{2,0}=\cfrac{HL^3}{3EI} \tag{II.P2.16}\label{II.P2.16}\]

• od momentu zginającego $M$

\[ \bar{w}_M= 1+ \cfrac{\Psi_{2}\left(\bar{C}_{\varphi}+4\Psi_{3}-2\Psi_{4}\right)} {D_M}. \tag{II.P2.17}\label{II.P2.17}\]

\[ D_M=D_H.\]

Wielkość jest znormalizowana względem klasycznego ugięcia

\[ w_{2,0}=\cfrac{ML^2}{2EI} \tag{II.P2.18}\label{II.P2.18}\]

• od obciążenie równomiernie rozłożonego $h$

\[ \cfrac{ 2\left[ -2\mu^{2}\left(\Psi_{3}\left(\bar C_{\varphi}+4\Psi_{3}\right)-\Psi_{4}^{\,2}\right) -3\bar C_{\varphi}\Psi_{2}\,\mu\cot\!\left(\cfrac{\mu}{2}\right) +6\bar C_{\varphi}\Psi_{2}
\right] } {D_h}. \tag{II.P2.19}\label{II.P2.19} \]

\[ D_h=3\mu^2D_H. \]

Wielkość znormalizowano względem klasycznego ugięcia

\[w_{2,0}=\cfrac{hL^4}{8EI}. \tag{II.P2.20}\label{II.P2.20}\]

Rozwiązanie można rozdzielić na dwie części o jednoznacznej interpretacji mechanicznej. Mianownik opisuje właściwości geometryczne i statecznościowe układu, natomiast licznik odzwierciedla wpływ rodzaju obciążenia na odpowiedź konstrukcji. Dzięki temu wszystkie analizowane przypadki wykazują tę samą osobliwość odpowiadającą utracie stateczności, różniąc się jedynie wartością przemieszczenia przed osiągnięciem stanu krytycznego.

Po podstawieniu definicji funkcji statecznościowych $\Psi_i $  otrzymuje się zwarte postacie i wykresy przedstawione w rozdziale Imperfekcje_przechylowe i ich_rownoważniki

Podobieństwo  statecznościowe 

Występuje zgodność mianowników

\[D_H=D_M=\cfrac{D_h}{3\mu^2}. \tag{II.P2.21}\label{II.P2.21} \]

Wzystkie trzy rozwiązania mają identyczną strukturę algebraiczną. Różnice pomiędzy rozwiązaniami wynikają wyłącznie z postaci liczników, natomiast mianownik opisuje wspólne właściwości geometryczne i statecznościowe układu. Oznacza to, że niezależnie od rodzaju obciążenia ($H$, $M$ lub $h$) utrata stateczności następuje przy spełnieniu tego samego warunku osobliwego. Wyjaśnia to duże podobieństwo przebiegów przedstawionych na rys. II.4–II.6, które różnią się jedynie wpływem rodzaju obciążenia na wartość ugięcia, zachowując tę samą osobliwość odpowiadającą utracie stateczności. O położeniu punktu osobliwego decydują wyłącznie parametry geometryczne i sztywnościowe elementu, opisane funkcjami statecznościowymi $\Psi_i$ oraz bezwymiarową sztywnością obrotową podpory $\bar C_\varphi$. Wyzerowanie mianownika odpowiada zanikowi wyznacznika macierzy sztywności, a tym samym osiągnięciu stanu krytycznego utraty stateczności.

Wpływ rodzaju obciążenie

Lliczniki poszczególnych rozwiązań zależą od rodzaju przyłożonego obciążenia i opisują sposób jego oddziaływania na układ. Określają one wartość przemieszczenia dla danego przypadku obciążenia, nie wpływając jednak na położenie punktu osobliwego rozwiązania. Zmiana schematu obciążenia prowadzi zatem do zmiany postaci licznika, natomiast wspólny mianownik pozostaje niezmienny.

Szczególowa analiza  rozwiązań  

Na rys. II.4 i II.5  przedstawiono przebieg bezwymiarowych wychyleń końca słupa wspornikowego wyznaczonych z zależności (II.P2.12) oraz (II.P2.13). 

Dla małych wartości parametru obciążenia \(\bar\Lambda\) wpływ siły osiowej na wychylenie wspornika jest niewielki, a odpowiedź układu pozostaje zbliżona do rozwiązania liniowego. Wraz ze wzrostem obciążenia ściskającego obserwuje się stopniowe zmniejszanie efektywnej sztywności konstrukcji, czego efektem jest coraz szybszy wzrost wychyleń. W pobliżu obciążenia krytycznego tempo przyrostu przemieszczeń gwałtownie rośnie, a krzywe przyjmują charakter asymptotyczny. Odpowiada to utracie stateczności układu określonej warunkiem \(\Theta=0\). Dla każdej wartości parametru \(\bar C_\varphi\) punkt krytyczny występuje przy innej wartości \(\bar\Lambda_{cr}\), co ilustrują pionowe linie przerywane. Zwiększenie sztywności obrotowej podpory powoduje przesunięcie punktu krytycznego w kierunku większych wartości parametru \(\bar\Lambda\), a tym samym wzrost nośności statecznościowej wspornika. Jednocześnie maleją wychylenia odpowiadające temu samemu poziomowi obciążenia osiowego, co świadczy o korzystnym wpływie częściowego lub pełnego utwierdzenia podstawy.

Porównanie obu pęków krzywych pokazuje ponadto, że jakościowy charakter odpowiedzi konstrukcji jest taki sam dla obciążenia skupionego $H$ i obciążenia rozłożonego $h$ oraz momentu $M$ . Dominujący wpływ na zachowanie układu wywierają parametr stateczności \(\bar\Lambda\) oraz sztywność obrotowa podpory \(\bar C_\varphi\).

Momemnt skupiony $M$, obciążenie równomiernie rozłożone \(h\) nie jest statycznie równoważne sile skupionej \(H\) przyłożonej na końcu wspornika, ze względu na osdmienny przebieg sił przkrojowych po długości pręta. 

Otrzymany wynik ma istotne znaczenie praktyczne dla metod wykorzystujących imperfekcje zastępcze. W projektowaniu konstrukcji często przyjmuje się wartości imperfekcji opracowane dla prętów ściskanych o schemacie belki swobodnie podpartej i stosuje je bez dodatkowej analizy również do słupów wspornikowych. Przeprowadzona analiza wskazuje jednak, że takie postępowanie nie posiada ścisłego uzasadnienia teoretycznego. Odpowiedź wspornika na obciążenie poprzeczne zależy nie tylko od wielkości obciążenia, lecz również od jego rozkładu, parametrów statecznościowych oraz warunków zamocowania. W konsekwencji nie istnieje uniwersalna imperfekcja zastępcza zapewniająca równoważność wszystkich schematów obciążenia. Oznacza to, że bezpośrednie przenoszenie wartości imperfekcji zastępczych wyznaczonych dla prętów swobodnie podpartych na słupy wspornikowe może prowadzić do błędnej oceny efektów drugiego rzędu oraz nośności statecznościowej. Problem ten staje się szczególnie istotny dla konstrukcji o dużej smukłości oraz dla układów o częściowo podatnym zamocowaniu.

W dalszej części pracy zagadnienie to zostanie przeanalizowane szczegółowo na podstawie ścisłych rozwiązań równań stateczności. Zostanie pokazane, że dla wsporników klasyczne procedury imperfekcji zastępczych mogą prowadzić do wyników istotnie odbiegających od rzeczywistej odpowiedzi konstrukcji.

Kąt obrotu w podstawie

Kąt obrotu w podstaie wspornika wynosi (dla $h$ i $H$ odpowiednio):

\[ \varphi_{2,H} = \cfrac{HL^2}{EI} \, \cfrac{ \Psi_2(2\Psi_3-\Psi_4) } {\Theta}. \tag{II.P2.22} \]

\[ \varphi_{2,h} = -\cfrac{hL^3}{6EI} \, \cfrac{ 3\Psi_2(\Psi_2-2\Psi_3+\Psi_4) -\Psi_1(2\Psi_3+\Psi_4) } {\Theta}. \tag{II.P2.23} \]

\[ \cfrac{\varphi_{ 2,h}}{\varphi_{2,H}} = -\cfrac{hL}{6H} \,\cfrac{ 3\Psi_2(\Psi_2-2\Psi_3+\Psi_4) -\Psi_1(2\Psi_3+\Psi_4) } { \Psi_2(2\Psi_3-\Psi_4) }.\tag{II.P2.24} \]

Iloraz kątów obrotu wywołąnych $h$ i $H$  (II.P2.22) wskazuje, że stosunek obrotów podstawy słupa wywołanych obciążeniem rozłożonym i skupionym nie zależy bezpośrednio od sztywności obrotowej podpory \(\bar C_\varphi\). Parametr ten wpływa wprawdzie na bezwzględne wartości obrotów poprzez współczynnik \(\Theta\), jednak jego wpływ znika po utworzeniu ilorazu. Oznacza to, że względna skuteczność obu schematów obciążenia w wywoływaniu obrotu podstawy jest określona wyłącznie przez funkcje statecznościowe \(\Psi_i\), a więc przez poziom obciążenia osiowego \(\bar\Lambda\).

Wynik ten jest szczególnie interesujący w kontekście metod imperfekcji zastępczych. Pokazuje on bowiem, że nawet dla wspornika z podatnym utwierdzeniem nie istnieje stały współczynnik pozwalający zastąpić obciążenie rozłożone równoważną siłą skupioną. Współczynnik taki zależy od aktualnego stanu pracy konstrukcji i zmienia się wraz ze wzrostem siły ściskającej.

Siły przekrojowe w utwierdzeniu 

Po wyznaczeniu przemieszczeń przywęzłowych siły przekrojowe można obliczyć z wykorzystaniem pierwotnej macierzy sztywności elementu (\ref{II.22}). Wektor sił przywęzłowych otrzymuje się z zależności

\[ \{P\} = [k]^{[e]}\{q\}. \tag{II.P2.25} \label{II.P2.25}\]

gdzie wektor przemieszczeń przywęzłowych $\{q \}$ wyznaczono  z zależności (\ref{II.P2.13}), a  $[k]^{[e]}$  jest pierwotną macierzą sztywności elementu [e]. 

W analizowanym przypadku wspornika szczególne znaczenie posiada moment zginający w utwierdzeniu \(M_2\), odpowiadający reakcji momentowej w podstawie słupa. Wartość ta uwzględnia zarówno wpływ obciążenia poprzecznego, jak i nieliniowych efektów drugiego rzędu wywołanych siłą osiową.
Ścisłe zależności opisujące moment utwierdzenia dla obciążenia skupionego $H$ oraz dla obciążenia równomiernie rozłożonego $h$ umożliwiają ocenę wpływu parametrów statecznościowych $ \Psi_i$, sztywności obrotowej podpory $\bar C_\varphi$ oraz obciążenia osiowego $\bar\Lambda$ na rozwój sił przekrojowych w  węźle podporowym (2) . Po wykonaniu przypisanych działań uzyskano:

\[ M_{2,H} = -\cfrac{H\,\bar C_\varphi\,\Psi_2(2\Psi_3-\Psi_4)} {\Theta}. \tag{II.P2.26} \]

\[ M_{2,h} = \cfrac{hL^2\,\bar C_\varphi} {6\,\Theta} \left[ 3\Psi_2(\Psi_2-2\Psi_3+\Psi_4) – \Psi_1(2\Psi_3+\Psi_4) \right]. \tag{II.P2.27} \]

\[ \cfrac{M_{2,h}}{M_{2,H}} = -\cfrac{hL}{6H} \, \cfrac{ 3\Psi_2(\Psi_2-2\Psi_3+\Psi_4) – \Psi_1(2\Psi_3+\Psi_4) } { \Psi_2(2\Psi_3-\Psi_4) }. \tag{II.P2.28} \]

Ponieważ moment utwierdzenia jest proporcjonalny do obrotu podstawy słupa, stosunek momentów wywołanych obciążeniem rozłożonym i skupionym jest identyczny jak stosunek odpowiadających im obrotów. Oznacza to, że brak ścisłej równoważności pomiędzy obciążeniem rozłożonym i skupionym dotyczy nie tylko przemieszczeń, lecz również sił przekrojowych. W szczególności nie istnieje uniwersalny współczynnik pozwalający zastąpić obciążenie rozłożone równoważną siłą skupioną lub równoważną imperfekcją geometryczną niezależnie od poziomu obciążenia osiowego. Współczynnik taki zależy od funkcji statecznościowych $\Psi_i$, a więc od aktualnego stanu pracy konstrukcji.

Przykład II.P3. [ Równoważniki imperfekcji przechyłowych wspornika ]

W przykładzie porównano trzy niezależne kryteria wyznaczania zastępczego mimośrodu równoważnego imperfekcji przechyłowej wspornika z przykładu II.P2 . Rozpatrzono trzy kryteria równoważności:

• kryterium przemieszczeniowego, polegającego na zrównaniu wychylenia wierzchołka wspornika,
• kryterium statycznego, polegającego na zrównaniu momentów zginających w podstawie wspornika,
• kryterium energetyczne – oparte na zrównaniu całkowitej energii odkształcenia obu modeli.

We wszystkich przypadkach porównywano referencyjny model imperfekcji reprezentowany przez zastępczą siłę poziomą \(H\) z modelem mimośrodowego przyłożenia osiowej siły ściskającej \(N\), dla którego moment zastępczy wynosi  $ M=N\,e.$

Pierwsze dwa kryteria wyprowadzono bezpośrednio na podstawie ścisłego rozwiązania przedstawionego w przykładzie II.P2. Natomiast kryterium energetyczne stanowi bezpośrednie zastosowanie ogólnej teorii
równoważności energetycznej przedstawionej w rozdziale „Równoważność energetyczna zastępczego mimośrodu”.

Równoważnośc przemieszczeniowa 

Obciążenie elementu osiową siłą ściskającą \(N\) przyłożoną z mimośrodem równoważnym \(e_{eq}\) jest statycznie równoważne układowi osiowej siły \(N\) oraz momentu zginająceg

\[ M=N\,e_{eq}. \tag{II.P3.1}\label{II.P3.1} \]

Oznacza to, że do dalszej analizy można bezpośrednio wykorzystać rozwiązanie dla momentu skupionego otrzymane w przykładzie II.P2. Porównując zależności (\ref{II.P2.15}) oraz (\ref{II.P2.17}) zauważa się, że oba rozwiązania posiadają identyczny mianownik \(D_H=D_M\), zgodnie z zależnością (\ref{II.P2.21}). Wynika stąd, że oba modele imperfekcji osiągają stan krytyczny przy tej samej wartości parametru stateczności, natomiast różnią się jedynie licznikami opisującymi wpływ rodzaju imperfekcji na wartość przemieszczenia.

Warunek równoważności przemieszczeniowej wymaga, aby oba modele imperfekcji wywoływały jednakowe wychylenie poziome wierzchołka wspornika

\[ \bar w_H=\bar w_M. \tag{II.P3.2}\label{II.P3.2} \]

Po podstawieniu zależności (\ref{II.P2.15}) oraz (\ref{II.P2.17}) otrzymuje si

\[ 1+ \cfrac{\Psi_{3}\left(\bar C_{\varphi}+4\Psi_{3}\right)-\Psi_{4}^{\,2}} {D_H} = 1+ \cfrac{\Psi_{2}\left(\bar C_{\varphi}+4\Psi_{3}-2\Psi_{4}\right)} {D_M}. \tag{II.P3.3}\label{II.P3.3} \]

Składnik równy jedności odpowiada wspólnej części szczególnej rozwiązania liniowego pierwszego rzędu, dlatego występuje po obu stronach równania i ulega redukcji. Ponadto, zgodnie z zależnością (\ref{II.P2.21}), zachodzi \(D_H=D_M\), dzięki czemu porównaniu podlegają wyłącznie części nieliniowe obu rozwiązań. Wprowadzając współczynnik równoważności przemieszczeniowej odpowiedzi wywołanej momentem zginającym względem odpowiedzi wywołanej referencyjną siłą poziomą, otrzymuje się zależnoś

\[ \Psi_{2}\left(\bar C_{\varphi}+4\Psi_{3}-2\Psi_{4}\right) = \eta_{(H=M)} \left[ \Psi_{3}\left(\bar C_{\varphi}+4\Psi_{3}\right)-\Psi_{4}^{\,2} \right]. \]

Stąd współczynnik równoważności przemieszczeniowej przyjmuje posta

\[ \eta_{(H=M)} = \cfrac{ \Psi_{2}\left(\bar C_{\varphi}+4\Psi_{3}-2\Psi_{4}\right) } { \Psi_{3}\left(\bar C_{\varphi}+4\Psi_{3}\right)-\Psi_{4}^{\,2} }. \tag{II.P3.4}\label{II.P3.4} \]

W praktyce projektowej imperfekcja geometryczna opisywana jest jednak mimośrodem równoważnym \(e_{eq}\). Przyjmując referencyjną imperfekcję w postaci zastępczej siły poziomej

\[ H=\cfrac{N}{n_L}. \tag{II.P3.5}\label{II.P3.5} \]

oraz wprowadzając bezwymiarowy mimośród

\[ \bar e=\cfrac{e_{eq}}{\Delta}, \qquad n_L=\cfrac{L}{\Delta}. \tag{II.P3.6}\label{II.P3.6} \]

otrzymuje się zależność

\[ \cfrac{M}{H} = \bar e\,L. \tag{II.P3.7}\label{II.P3.7} \]

Uwzględniając wielkości normalizujące zastosowane w przykładzie II.P2,

\[ w_{H,0}=\cfrac{HL^{3}}{3EI}, \qquad w_{M,0}=\cfrac{ML^{2}}{2EI}, \]

otrzymuje si

\[ \cfrac{w_{M,0}}{w_{H,0}} = \cfrac{3}{2}\,\bar e. \tag{II.P3.8}\label{II.P3.8} \]

Ostatecznie współczynnik równoważności mimośrodu względem referencyjnej siły poziomej przyjmuje postać

\[ \eta_{(e=H)} = \cfrac{1} {1+\dfrac{3}{2}\,\bar e\,\eta_{(H=M)}}. \tag{II.P3.9}\label{II.P3.9} \]

gdzie współczynnik \(\eta_{(H=M)}\) określono zależnością (\ref{II.P3.4}).

Zależność (\ref{II.P3.9}) stanowi ścisłe rozwiązanie przyjętego kryterium równoważności przemieszczeniowej i będzie podstawą dalszej analizy wpływu mimośrodu, smukłości elementu oraz podatności obrotowej podpory na równoważność modeli imperfekcji.

Równoważnośc statyczna – moementów u podstawy

Oprócz równoważności przemieszczeniowej można rozpatrywać również równoważność statyczną modeli imperfekcji. Kryterium to polega na porównaniu momentów zginających w przekroju utwierdzenia wywołanych referencyjną siłą poziomą oraz mimośrodem równoważnym. Moment u podstawy jest bowiem podstawowym parametrem wykorzystywanym podczas wymiarowania przekroju i połączenia z fundamentem.

Dla referencyjnej imperfekcji w postaci zastępczej siły poziomej moment u podstawy wynosi 

\[ M_H=HL. \tag{II.P3.10}\label{II.P3.10} \]

Natomiast dla osiowej siły ściskającej \(N\) przyłożonej z mimośrodem równoważnym \(e_{eq}\)

\[ M_e=N\,e_{eq}. \tag{II.P3.11}\label{II.P3.11} \]

Po podstawieniu zależności (\ref{II.P3.5}) oraz (\ref{II.P3.6}) otrzymuje się

\[ \cfrac{M_e}{M_H} = \bar e. \tag{II.P3.12}\label{II.P3.12} \]

Wprowadzając współczynnik równoważności statycznej mimośrodu względem referencyjnej siły poziomej zdefiniowany jako stosunek momentów zginających w przekroju utwierdzenia

\[ \eta_{(e=H)} = \cfrac{M_e}{M_H}, \tag{II.P3.13}\label{II.P3.13} \]

oraz uwzględniając zależność (\ref{II.P3.12}), otrzymuje się

\[ \eta_{(e=H)} = \bar e. \tag{II.P3.14}\label{II.P3.14} \]

Równoważnośc energetyczna i  porównanie kryteriów równoważności

Kryterium róęnowżności energetycznej wyczerpująco przedtsawiono w tekscie . 

Przeprowadzona analiza wykazała, że wartość zastępczego mimośrodu zależy od przyjętego kryterium równoważności. W kryterium przemieszczeniowym zachowana zostaje zgodność wychylenia wierzchołka wspornika, natomiast w kryterium statycznym zapewniona jest zgodność momentów zginających w przekroju utwierdzenia. Kryterium energetyczne prowadzi natomiast do zgodności całkowitej energii odkształcenia obu modeli, uwzględniając jednocześnie wpływ rozkładu sił wewnętrznych na całej długości elementu.

W przeciwieństwie do dwóch pierwszych kryteriów, równoważność energetyczna nie opiera się na porównaniu pojedynczej wielkości kinematycznej lub statycznej, lecz na porównaniu globalnej odpowiedzi układu. Z tego względu może być traktowana jako najbardziej ogólne kryterium wyznaczania zastępczego mimośrodu. Najważniejsze zależności uzyskane dla poszczególnych kryteriów zestawiono w tabeli II.P3.1.

Tab. II.P3.1 Porównanie  kryteruów róenoważnosci

\[ \begin{array}{|c|c|c|c|}
\hline \textbf{Kryterium} & \textbf{Warunek} & \textbf{Współczynnik} & \textbf{Zależność wyznaczająca } \bar e_{eq} \\
\hline \text{Przemieszczeniowe} & w_H=w_e & \eta_w & \eta_{(e=H)} = \cfrac{1} {1+\dfrac{3}{2}\,\bar e\,\eta_{(H=M)}} \;(\ref{II.P3.9}) \\
\hline \text{Statyczne} & M_H=M_e & \eta_M & \eta_{(e=H)} = \bar e \;(\ref{II.P3.14}) \\ 
\hline \text{Energetyczne} & \Pi_H=\Pi_e & \eta_E & \bar e_{eq} = \sqrt{ \cfrac{ \Psi_3\!\left(\bar C_{\varphi}+4\Psi_3\right)-\Psi_4^{\,2} } { 3\!\left( -3\Psi_2^{\,2} +\bar C_{\varphi}\Psi_1 +4\Psi_1\Psi_3 \right)}
} \;(\ref{II.77}) \\
\hline \end{array} \]

Wnioski

(1) Przedstawiona analiza wykazała, że wartość zastępczego mimośrodu zależy od przyjętego kryterium równoważności. Kryterium przemieszczeniowe zapewnia zgodność odpowiedzi kinematycznej elementu, kryterium statyczne prowadzi do zgodności momentów zginających w przekroju utwierdzenia, natomiast kryterium energetyczne zapewnia zgodność całkowitej energii odkształcenia obu modeli.
(2) W przeciwieństwie do kryteriów lokalnych, opartych na porównaniu pojedynczej wielkości kinematycznej lub statycznej, kryterium energetyczne uwzględnia globalną odpowiedź konstrukcji wynikającą z rozkładu sił wewnętrznych na całej długości elementu. Dzięki wykorzystaniu ścisłych funkcji kształtu Livesleya oraz funkcji stateczności \(\Psi_i\) możliwe było wyprowadzenie zamkniętego wzoru na równoważny mimośród bez stosowania aproksymacji charakterystycznych dla klasycznej metody Ritza.
(3) Otrzymane zależności stanowią podstawę do wyznaczania zastępczych mimośrodów odpowiadających różnym modelom imperfekcji geometrycznych oraz mogą być wykorzystane zarówno w analizie pojedynczych elementów prętowych, jak i w geometrycznie nieliniowych modelach MES.

Przykład II.P4 [ Analiza belki-słupa metodą Ritza ]

Zastosować metodę Ritza do analizy belki przedstawionej na rys. II.1b, dla której w przykładzie rys. II.1b, dla której w przykładzie II.P1 uzyskano rozwiązanie ścisłe. Obliczenia przeprowadzićz wykorzystaniem funkcjonału Lagrange’a (\ref{II.111}). Belka-słup jest ściskana osiowo siłą (N) oraz podparta sprężyście w kierunku pionowym podporą zlokalizowaną w węźle (2). Sztywność podpory opisuje bezwymiarowy parametr
$ \bar C_\Delta=\frac{L^3}{EI},C_\Delta,$ , gdzie $$C_\Delta$ oznacza sztywność pionową podpory sprężystej.

W celu zachowania ogólności uwzględniono jednoczesne działanie siły osiowej $N$ oraz równomiernie rozłożonego obciążenia poprzecznego $q$. Pozwala to na równoczesne wyznaczenie ugięcia belki oraz ocenę wpływu siły osiowej na odpowiedź konstrukcji. Szczególną uwagę poświęcono wpływowi doboru funkcji bazowych oraz liczby wyrazów aproksymacji na dokładność otrzymywanych rozwiązań.

Dobór funkcji bazowych

Funkcję aproksymującą ugięcie przyjęto w postaci określonej zależnością (\ref{II.116})

\[ w(\xi)=a_1\varphi_1(\xi)+a_2\varphi_2(\xi)+a_3\varphi_3(\xi)+a_4\varphi_4(\xi),
\tag{II.P2.1}\label{II.P2.1} \]

gdzie (\xi=x/L) jest bezwymiarową współrzędną wzdłuż długości belki, (a_i;(i=1,\ldots,4)) są współczynnikami aproksymacji, natomiast (\varphi_i(\xi)) oznaczają funkcje bazowe.

Funkcje bazowe wyznaczono z funkcji generatorowej przyjętej w postaci wielomianu czwartego stopnia

\[ \varphi(\xi)=a+b\xi+c\xi^2+d\xi^3+e\xi^4, \tag{II.P2.2}\label{II.P2.2} \]

której współczynniki wyznaczono z warunków kinematycznych analizowanego zagadnienia. W rozpatrywanym przypadku jedynym geometrycznym warunkiem brzegowym jest nieprzemieszczalność podpory w węźle (1). Pozostałym stopniom swobody nie przypisano geometrycznych warunków brzegowych, dzięki czemu ich wartości zostaną wyznaczone z warunku stacjonarności funkcjonału Lagrange’a. Warunki kinematyczne można zatem zapisać w postaci

\ [ w(0)=0,\qquad w'(0)\neq0,\qquad w(1)\neq0,\qquad w'(1)\neq0. \tag{II.P2.3}\label{II.P2.3}]

Warunek (w(1)\neq0) wynika z obecności podpory sprężystej. W przeciwieństwie do podpory nieprzemieszczalnej nie narzuca ona zerowego przemieszczenia pionowego, lecz jedynie wytwarza reakcję proporcjonalną do przemieszczenia. Wartość (w(1)) nie jest zatem znana a priori, lecz zostanie wyznaczona z warunku stacjonarności funkcjonału Lagrange’a.

Uwzględnienie geometrycznego warunku brzegowego (\ref{II.P2.3}) prowadzi do zależności

\[ a=0, \tag{II.P2.4}\label{II.P2.4} \]

a funkcja generatorowa przyjmuje postać

\[ \varphi(\xi)=b\xi+c\xi^2+d\xi^3+e\xi^4. \tag{II.P2.5}\label{II.P2.5} \]

Funkcje bazowe (\varphi_i(\xi)) wyznacza się przez kolejne przyjęcie jednostkowej wartości współczynników (b), (c), (d) oraz (e), przy jednoczesnym wyzerowaniu pozostałych współczynników. W rezultacie otrzymuje się cztery liniowo niezależne funkcje bazowe, które zostaną wykorzystane do wyznaczenia współczynników metody Ritza   w posataci  funkcji aproksymującej (\refl{II.P2.1}).

Wyznaczenie współczynników Ritza  

Przykład jest ilustracją do   rozdziału  Układ kanoniczny metody Ritza w  ujeciu klasycznym (bez elementów dyskretnych)

Ponieważ funkcje bazowe przyjęto w postaci wielomianów, ich pierwsze i drugie pochodne wyznacza się przez bezpośrednie różniczkowanie. Pochodne te zostaną następnie wykorzystane do wyznaczenia współczynników \(\delta_{ij}\) zgodnie z zależnością (\ref{II.125}). Otrzymuje się

\[ \begin{aligned} \varphi_1′(\xi)&=\frac{d}{d\xi}\left(\xi\right)=1, & \varphi_2′(\xi)&=\frac{d}{d\xi}\left(\xi^2\right)=2\xi, & \varphi_3′(\xi)&=\frac{d}{d\xi}\left(\xi^3\right)=3\xi^2, & \varphi_4′(\xi)&=\frac{d}{d\xi}\left(\xi^4\right)=4\xi^3,\\[2mm]
\varphi_1”(\xi)&=\frac{d^2}{d\xi^2}\left(\xi\right)=0, & \varphi_2”(\xi)&=\frac{d^2}{d\xi^2}\left(\xi^2\right)=2, & \varphi_3”(\xi)&=\frac{d^2}{d\xi^2}\left(\xi^3\right)=6\xi, & \varphi_4”(\xi)&=\frac{d^2}{d\xi^2}\left(\xi^4\right)=12\xi^2. \end{aligned} \tag{II.P2.9}\label{II.P2.9} \]

Po podstawieniu pochodnych funkcji bazowych do zależności (\ref{II.P2.7}) i (\ref{II.P2.8}) otrzymuje się pierwszą i drugą pochodną funkcji aproksymującej w postaci

\[w'(\xi)=\frac{dw}{d\xi} =a_1\varphi_1′(\xi)+a_2\varphi_2′(\xi)+a_3\varphi_3′(\xi)+a_4\varphi_4′(\xi) =a_1+2a_2\xi+3a_3\xi^2+4a_4\xi^3. \tag{II.P2.7}\label{II.P2.7} \]

\[ w”(\xi)=\frac{d^2w}{d\xi^2} =a_1\varphi_1”(\xi)+a_2\varphi_2”(\xi)+a_3\varphi_3”(\xi)+a_4\varphi_4”(\xi) =2a_2+6a_3\xi+12a_4\xi^2. \tag{II.P2.8}\label{II.P2.8}\]

Przemieszczenie pionowe węzła (2), w którym zlokalizowana jest podpora sprężysta, jest niezbędne do wyznaczenia współczynników macierzy \([K_C]\) zgodnie z zależnością (\ref{II.140}). Podstawiając do funkcji aproksymującej (\ref{II.P2.1}) wartości funkcji bazowych w punkcie \(\xi=1\) $  \varphi_1(1)=1,\qquad \varphi_2(1)=1,\qquad \varphi_3(1)=1,\qquad \varphi_4(1)=1, $ otrzymuje się

\[w(1)=a_1+a_2+a_3+a_4. \tag{II.P2.10}\label{II.P2.10}\]

Wyznaczone pochodne stanowią podstawę do obliczenia współczynników \(\delta_{ij}\) zgodnie z zależnością (\ref{II.125}) oraz składowych wektora obciążeń uogólnionych \(\Delta_{Fi}\) zgodnie z zależnością (\ref{II.126}). Korzystając z zależności (\ref{II.137}) wyznacza się współczynniki macierzy sztywności sprężystej \([K_e]\), odpowiadającej energii odkształcenia przy zginaniu. Po wykonaniu całkowania otrzymuje się

Współczynniki macierzy Ritza wyznacza się z zależności (\ref{II.125}) oraz (\ref{II.126}). Po podstawieniu wyznaczonych funkcji bazowych współczynniki macierzy Ritza zostały wyznaczone symbolicznie z wykorzystaniem programu Mathematica. Pozwoliło to na analityczne obliczenie całek występujących w równaniach oraz automatyczne uproszczenie otrzymanych wyrażeń.

Korzystając z zależności (\ref{II.137}) wyznacza się współczynniki macierzy sztywności sprężystej \([K_e]\), odpowiadającej energii odkształcenia przy zginaniu. Po podstawieniu przyjętych funkcji bazowych oraz wykonaniu całkowania otrzymuje się

\[[K_e] = \frac{EI}{L^3} \begin{bmatrix}
0 & 0 & 0 & 0\\[2mm]
0 & 4 & 6 & 8\\[2mm]
0 & 6 & 12 & 18\\[2mm]
0 & 8 & 18 & \dfrac{144}{5}
\end{bmatrix}. \tag{II.P2.11}\label{II.P2.11} \]

Korzystając z zależności (\ref{II.118}) wyznacza się współczynniki macierzy sztywności geometrycznej \([K_g]\), odpowiadającej energii potencjalnej siły osiowej. Wprowadzając bezwymiarowy parametr obciążenia

\[N=\bar{\Lambda}\frac{EI}{L^2}, \]

po wykonaniu całkowania otrzymuje się

\[ [K_g]=-\frac{\bar{\Lambda}EI}{L^3}
\begin{bmatrix}1 & 1 & 1 & 1\\[2mm]
1 & \dfrac{4}{3} & \dfrac{3}{2} & \dfrac{8}{5}\\[2mm]
1 & \dfrac{3}{2} & \dfrac{9}{5} & 2\\[2mm]
1 & \dfrac{8}{5} & 2 & \dfrac{16}{7}
\end{bmatrix}. \tag{II.P2.12}\label{II.P2.12}\]

Korzystając z zależności (\ref{II.140}) wyznacza się współczynniki macierzy sztywności elementu dyskretnego \([K_C]\), odpowiadającej energii odkształcenia podpory sprężystej. Uwzględniając zależność (\ref{II.P2.10}) oraz bezwymiarowy parametr sztywności podpory

\[C_\Delta=\bar C_\Delta\frac{EI}{L^3},\]

otrzymuje się

\[[K_C]= \frac{\bar C_\Delta EI}{L^3}
\begin{bmatrix}1 & 1 & 1 & 1\\[2mm]
1 & 1 & 1 & 1\\[2mm]1 & 1 & 1 & 1\\[2mm]
1 & 1 & 1 & 1\end{bmatrix}.
\tag{II.P2.13}\label{II.P2.13}\]

Korzystając z zależności (\ref{II.141}) całkowitą macierz współczynników metody Ritza otrzymuje się przez zsumowanie macierzy sztywności sprężystej, geometrycznej oraz macierzy elementu dyskretnego. Otrzymuje się

\[ [\delta] = [K_e] + [K_g] + [K_C]. \tag{II.P2.14}\label{II.P2.14} \]

Korzystając z zależności (\ref{II.126}) wyznacza się składowe wektora obciążeń uogólnionych \(\{\Delta_F\}\), odpowiadającego pracy równomiernie rozłożonego obciążenia poprzecznego. W analizowanym przykładzie przyjęto dodatni zwrot obciążenia \(q\) zgodny z dodatnim zwrotem osi \(w\). W dalszej części, podczas porównania z rozwiązaniem ścisłym przykładu II.P1, zostanie przyjęta zależność

\[ q=-p, \]

gdzie \(p\) oznacza intensywność obciążenia stosowaną w rozwiązaniu ścisłym. Po podstawieniu funkcji bazowych oraz wykonaniu całkowania otrzymuje się

Korzystając z zależności (\ref{II.126}) wyznacza się wektor obciążeń uogólnionych. Dla równomiernie rozłożonego obciążenia poprzecznego \(q=\mathrm{const}\) otrzymuje się

\[ \{\Delta_F\}
= – qL \begin{Bmatrix} \dfrac12\\[2mm] \dfrac13\\[2mm] \dfrac14\\[2mm] \dfrac15 \end{Bmatrix}. \tag{II.P2.14}\label{II.P2.14} \]

Po wyznaczeniu współczynników aproksymacji oblicza się linię ugięcia belki z zależności (\ref{II.P2.1}), a następnie strzałkę ugięcia w środku rozpiętości

\[ f=w\!\left(\frac{L}{2}\right). \tag{II.P2.16}\label{II.P2.16}\]

W celu porównania z rozwiązaniem ścisłym wprowadza się bezwymiarową strzałkę ugięcia

\[ \bar f=\frac{fEI}{qL^4}. \tag{II.P2.17}\label{II.P2.17}\]

Na rysunku II.P2.1 przedstawiono porównanie rozwiązania otrzymanego metodą Ritza z rozwiązaniem ścisłym przykładu II.P1. Wraz ze wzrostem liczby funkcji bazowych dokładność aproksymacji systematycznie wzrasta, a dla przyjętej czterowyrazowej aproksymacji uzyskano bardzo dobrą zgodność z rozwiązaniem ścisłym.

Wyznaczenie obciążenia krytycznego

Po zsumowaniu macierzy (\ref{II.P2.11})–(\ref{II.P2.13}) zgodnie z zależnością (\ref{II.141}) otrzymuje się macierz współczynników metody Ritza (\ref{II.P2.15}). Wyznacznik tej macierzy ma postać

\[ \det[\delta] = \frac{EI^{4}}{10500\,L^{12}} (\bar{\Lambda}-60) (\bar{\Lambda}-\bar{C}_{\Delta}) \left(\bar{\Lambda}^{2}-180\bar{\Lambda}+1680\right). \tag{II.P2.16}\label{II.P2.16} \]

Ponieważ czynnik \(\dfrac{EI^{4}}{10500\,L^{12}}\) jest różny od zera, warunek utraty stateczności (\ref{II.146}) prowadzi do równania charakterystycznego

\[ (\bar{\Lambda}-60) (\bar{\Lambda}-\bar{C}_{\Delta}) \left(\bar{\Lambda}^{2}-180\bar{\Lambda}+1680\right)=0. \tag{II.P2.17}\label{II.P2.17} \]

Rozwiązaniem równania charakterystycznego są cztery wartości własne

\[ \bar{\Lambda}_1=60, \] \[ \bar{\Lambda}_2=2\left(45-\sqrt{1605}\right)\approx 9.875, \] \[ \bar{\Lambda}_3=2\left(45+\sqrt{1605}\right)\approx 170.125, \]\[ \bar{\Lambda}_4=\bar{C}_{\Delta}. \tag{II.P2.18}\label{II.P2.18} \]

O obciążeniu krytycznym decyduje najmniejsza dodatnia wartość własna. Dla analizowanego układu zależy ona od sztywności podpory sprężystej \(\bar{C}_{\Delta}\). Jeżeli \(\bar{C}_{\Delta}<9.875\), to wartością krytyczną jest \(\bar{\Lambda}=\bar{C}_{\Delta}\). W przeciwnym przypadku obciążenie krytyczne osiąga wartość

\[\bar{\Lambda}_{cr}=2\left(45-\sqrt{1605}\right)\approx 9.875. \]

Uzyskane rozwiązanie stanowi czterowyrazową aproksymację metody Ritza, której dokładność zostanie oceniona przez porównanie z rozwiązaniem ścisłym przedstawionym w przykładzie II.P1.

Przykład II.P5 [ Belka krzywoliniowa w kształcie sinusoidy]

W Dodatku II.A  przedstawiono ogólną procedurę Livesleya wyznaczania macierzy podatności elementów prętowych o dowolnej osi krzywoliniowej. W szczególności wyprowadzono ogólną postać energii odkształcenia (\ref{II.A.6}), definicję macierzy podatności wynikającą z twierdzenia Castigliano (\ref{II.A.7}) – (\ref{II.A.9}) oraz jawne zależności dla elementu o osi sinusoidalnej. W niniejszym przykładzie pokazano zastosowanie tej procedury do wyznaczenia macierzy podatności i odpowiadającej jej macierzy sztywności elementu belki Bernoulliego o zadanej geometrii początkowej.

Rozważyć element belki Bernoulliego o długości \(L\), którego oś środkowa posiada początkową imperfekcję w postaci pojedynczej lub wielokrotnej półfali sinusoidalnej przedstawionej na  na rys. II. 5  z amplitudą   

\[ f_0=\cfrac{L}{n_L}, \]

gdzie:
$L$ – długość elementu w rzucie,
$n_L$ – współczynnikimperfekcji łukowej.  W praktyce obliczeniowej przyjmuje się zwykle \(n_L\approx 200\), co odpowiada początkowej amplitudzie wygięcia rzędu $(L/200)$.

Równanie osi elementu

Oś środkową elementu opisuje zależność

\[ y(x)=f_0\sin(\omega x), \tag{II.P5.1} \label{II.P5.1} \]

gdzie:
$ \omega=\cfrac{n \pi}{L},$
$ n=1,2,3,\ldots$ – numer postaci własnej opisującej kształt imperfekcji.

Równanie (\ref{II.P5.1}) jest szczególnym przypadkiem ogólnej geometrii elementu opisanej zależnościami (II.A.1)–(II.A.3), dla których w Dodatku II.A wyprowadzono pełne zależności dla elementu imperfekcyjnego.

Do obliczeń wygodnie jest wprowadzić parametr

\[  x=\cfrac{\alpha}{\omega}, \qquad y=f_0\sin\alpha, \tag{II.P5.2} \label{II.P5.2} \]

gdzie $ \alpha\in \left[ -\cfrac{n\pi}{2}, \; \cfrac{n\pi}{2} \right]. $

Macierz podatności

Zgodnie z procedurą przedstawioną w Dodatku II.A energia odkształcenia elementu opisana jest zależnością (II.A.6), natomiast macierz podatności wyznacza się z twierdzenia Castigliano zgodnie z równaniami (\ref{II.A.7}) – (\ref{II.A.9}).

Po podstawieniu równania osi (\ref{II.P5.1}) do ogólnych zależności Dodatku II.A oraz wykonaniu odpowiednich całkowań otrzymuje się współczynniki macierzy podatności.

\[ k_{11} = \cfrac{f_0^2\left(n\pi-\sin(n\pi)\right)} {2EI}, \tag{II.P5.3} \label{II.P5.3} \]

\[ k_{12} = \cfrac{f_0L \left[ n\pi\cos\left(\cfrac{n\pi}{2}\right) – 2\sin\left(\cfrac{n\pi}{2}\right) \right]} {n\pi EI}, \tag{II.P5.4} \label{II.P5.4} \]

\[ k_{22} = \cfrac{L}{12} \left[ \cfrac{6f_0^2\left(n\pi+\sin(n\pi)\right)} {EA} + \cfrac{L^2n\pi} {EI} \right], \tag{II.P5.5} \label{II.P5.5} \]

\[ k_{13} = k_{23}=0, \tag{II.P5.6} \label{II.P5.6} \]

\[ k_{33} = \cfrac{n\pi}{EI}. \tag{II.P5.7} \label{II.P5.7} \]

W rezultacie macierz podatności przyjmuje postać
\[  [F] = \begin{bmatrix} k_{11}&k_{12}&0\\ k_{12}&k_{22}&0\\
0&0&k_{33} \end{bmatrix}.\tag{II.P5.8} \label{II.P5.8} \]

Otrzymana macierz stanowi szczególny przypadek ogólnej macierzy podatności elementu imperfekcyjnego wyprowadzonej w Dodatku II.A. W przeciwieństwie do rozwiązania ogólnego, współczynniki macierzy zostały tutaj wyznaczone w postaci zamkniętej dla przyjętego przebiegu sinusoidalnego.

Macierz sztywności

Zgodnie z ogólną zależnością (\ref{II.A.53})Macierz sztywności wyznacza się przez odwrócenie macierzy podatności

\[ [K]=[F]^{-1}, \]

Dla pierwszej postaci imperfekcji (\(n=1\)) oraz przy założeniu stałych parametrów przekroju \(EA=\mathrm{const}\) i \(EI=\mathrm{const}\), odwrócenie macierzy podatności (\ref{II.P5.8}) prowadzi do zwartej postaci macierzy sztywności

\[ [K]= \cfrac{1}{\Theta} \begin{bmatrix} \dfrac{2\pi^{3}EI\left(EAL^{2}+6EI\,f_{0}^{2}\right)}{f_{0}^{2}} & \dfrac{48\pi EAEI}{f_{0}} & 0 \\[3mm]
& \dfrac{12\pi^{3}EAEI}{L} & 0 \\
& SYM & \dfrac{EI\,\Theta}{\pi} \end{bmatrix}, \tag{II.P5.9} \label{II.P5.9} \]

gdzie

\[ \Theta= \pi^{4}EAL^{2} – 96EAL + 6\pi^{4}EI\,f_{0}^{2}. \tag{II.P5.10} \label{II.P5.10} \]

Przedstawiona macierz stanowi szczególny przypadek ogólnej macierzy sztywności elementu imperfekcyjnego Livesleya omówionego opracowanego w Dodatku II.A dla osi elementu opisanej pojedynczą półfalą sinusoidalną. Zwarta postać zależności (\ref{II.P5.9})–(\ref{II.P5.10}) umożliwia bezpośrednie wykorzystanie elementu w analizie belek i belek-słupów z imperfekcjami geometrycznymi oraz stanowi punkt wyjścia do badania wpływu amplitudy imperfekcji na sztywność i stateczność konstrukcji. Analiza tych zagadnień wykracza jednak poza zakres niniejszego opracowania.

Dodatek II.A. Element imperfekcyjny Livesleya

Ogólna procedura wyznaczania macierzy podatności elementów krzywoliniowych została przedstawiona przez Livesleya [Livesley R.K. (1975), Matrix Methods of Structural Analysis, 2nd Edition, Pergamon Press, Oxford, rozdz. 3.6]. W najogólniejszym przypadku geometria elementu opisywana jest funkcjami parametrycznymi osi środkowej, natomiast współczynniki macierzy podatności wyznaczane są z energii odkształcenia przy wykorzystaniu twierdzenia Castigliano.

Rozważmy element o osi środkowej opisanej funkcją sinusoidalną

\[ y(x)=e_0\sin\left(\cfrac{\pi x}{L}\right) \tag{II.A.1}\label{II.A.1}\]

gdzie $e_0$ oznacza amplitudę imperfekcji geometrycznej, a $L$ długość elementu. Funkcja ta nie opisuje rzeczywistego przebiegu imperfekcji, lecz stanowi model zastępczy odpowiadający pierwszej postaci wyboczeniowej pręta przegubowo podpartego. Przyjęcie pojedynczej półfali sinusoidalnej pozwala uchwycić dominujący mechanizm utraty stateczności przy zachowaniu stosunkowo prostej postaci analitycznej.

W metodzie Livesleya wszystkie wielkości geometryczne odnoszone są do rzeczywistej osi elementu. W przeciwieństwie do klasycznej teorii belki prostej długość elementu, kierunki lokalnych osi oraz energia odkształcenia wyznaczane są względem osi krzywoliniowej.

\[  \cfrac{dy}{dx} = \cfrac{\pi e_0}{L} \cos\left(\cfrac{\pi x}{L}\right) \tag{II.A.2}\label{II.A.2} \]

\[ ds= \sqrt{ 1+ \left( \cfrac{\pi e_0}{L} \right)^2 \cos^2 \left( \cfrac{\pi x}{L} \right) } \,dx \tag{II.A.3}\label{II.A.3}\]

Zależność (\ref{II.A.3}) określa elementarną długość rzeczywistej osi środkowej i stanowi punkt wyjścia do energetycznego opisu elementu imperfekcyjnego.

Energia odkształcenia

Rozważmy element obciążony siłą osiową $N$, siłą poprzeczną $V$ oraz momentem końcowym $M$. W dowolnym punkcie osi moment zginający wynosi

\[  M(x)= M+Vx- Ne_0 \sin\left( \cfrac{\pi x}{L} \right) \tag{II.A.4}\label{II.A.4} \]

Całkowita energia odkształcenia elementu jest sumą energii odkształceń osiowych i energii zginania

\[ U =  \int_0^L \cfrac{N^2}{2EA} \,ds + \int_0^L \cfrac{M^2(x)}{2EI} \,ds \tag{II.A.5}\label{II.A.5} \]

Po uwzględnieniu zależności (\ref{II.A.3}) i (\ref{II.A.4}) otrzymujemy

\[ U = \int_0^L \left \{ \cfrac{N^2}{2EA}  + \cfrac{1}{2EI} \left[ M+Vx- Ne_0 \sin\left( \cfrac{\pi x}{L} \right) \right]^2 \right\} \sqrt{ 1+ \left( \cfrac{\pi e_0}{L} \right)^2 \cos^2 \left( \cfrac{\pi x}{L} \right) } \,d\, x \tag{II.A.6}\label{II.A.6}\]

Energia oskształcenia   (\ref{II.A.6})  stanowi pełne energetyczne sformułowanie elementu imperfekcyjnego o osi sinusoidalnej. Na jej podstawie można wyznaczyć współczynniki macierzy podatności bez wprowadzania dodatkowych uproszczeń geometrycznych.

Macierz podatności elementu

Wektor sił uogólnionych zapisujemy w postaci $ \mathbf{P} = |N, \,  V, \,  M|^T $, 

natomiast odpowiadający mu wektor przemieszczeń węzłowych  $\mathbf{q} = |u, \, v , \, \varphi |^T $

Zgodnie z twierdzeniem Castigliano zachodzi:

\[ \mathbf{q}  = [F] \, \mathbf{P} \tag{II.A.7}\label{II.A.7}\]

gdzie $[F]$ jest macierzą podatności elementu,  której poszczególne składowe  wyznacza się poprzez drugie pochodne energii odkształcenia względem sił uogólnionych

\[  f_{ij} = \cfrac{\partial^2 U} {\partial P_i\,\partial P_j} \tag{II.A.8}\label{II.A.8} \] 

Energia odkształcenia osiowego zależy wyłącznie od siły osiowej $N$, natomiast energia odkształcenia zginania zawiera wszystkie trzy siły uogólnione, więc macierz podatności przyjmie  postać

\[  [F] = \begin{bmatrix} f_{11} & f_{12} & f_{13} \\
f_{12} & f_{22} & f_{23} \\
f_{13} & f_{23} & f_{33} \end{bmatrix}
= \begin{bmatrix}
\dfrac{\partial^2 U}{\partial N^2} & \dfrac{\partial^2 U}{\partial N\,\partial V} & \dfrac{\partial^2 U}{\partial N\,\partial M}\\[4mm]
\dfrac{\partial^2 U}{\partial N\,\partial V} & \dfrac{\partial^2 U}{\partial V^2} & \dfrac{\partial^2 U}{\partial V\,\partial M} \\[4mm]
\dfrac{\partial^2 U}{\partial N\,\partial M} & \dfrac{\partial^2 U}{\partial V\,\partial M} & \dfrac{\partial^2 U}{\partial M^2}\end{bmatrix}
\tag{II.A.9}\label{II.A.9} \]Po obliczeniu odpowiednich całek otrzymuje się pełną macierz podatności elementu imperfekcyjnego.

Jawna postać współczynników macierzy podatności

Po podstawieniu zależności (\ref{II.A.6}) do definicji (\ref{II.A.8}) otrzymujemy układ całek zależnych wyłącznie od parametrów geometrycznych elementu oraz amplitudy imperfekcji $e_0$. W przeciwieństwie do klasycznej belki prostoliniowej pojawiają się dodatkowe wyrazy sprzęgające oddziaływania osiowe i zginające, będące bezpośrednią konsekwencją początkowej krzywizny osi środkowej. Dla części współczynników możliwe jest uzyskanie postaci zamkniętej wyrażonej za pomocą funkcji eliptycznych i funkcji elementarnych, natomiast pozostałe zachowują postać całkową.

Ze względu na symetrię macierzy podatności wynikającą z twierdzenia Maxwella–Bettiego wystarczy wyznaczyć sześć niezależnych współczynników.

$ f_{11} = \cfrac{1}{EA} \int_0^L \sqrt{ 1+k^2\cos^2\left(\cfrac{\pi x}{L}\right) } \,dx + \cfrac{e_0^2}{EI} \int_0^L \sin^2\left(\cfrac{\pi x}{L}\right) \sqrt{ 1+k^2\cos^2\left(\cfrac{\pi x}{L}\right)} \,dx =\\
\cfrac{2L}{\pi EA} \sqrt{1+k^2}\;E(m)+ \cfrac{2Le_0^2}{\pi EI} \sqrt{1+k^2} \left[ E(m) – \cfrac{E(m)-(1-m)K(m)}{m} \right]\tag{II.A.10}\label{II.A.10} $

$f_{12} = – \cfrac{e_0}{EI} \int_0^L x \sin\left( \cfrac{\pi x}{L} \right) \sqrt{ 1+ \left( \cfrac{\pi e_0}{L} \right)^2 \cos^2\left( \cfrac{\pi x}{L} \right) } \,dx\tag{II.A.11}\label{II.A.11}$

$ f_{13} = – \cfrac{e_0}{EI} \int_0^L \sin\left( \cfrac{\pi x}{L} \right) \sqrt{ 1+ \left( \cfrac{\pi e_0}{L} \right)^2 \cos^2\left( \cfrac{\pi x}{L} \right) } \,dx=  – \cfrac{e_0L}{\pi EI} \left[ \sqrt{1+k^2} + \cfrac{\operatorname{arsinh}(k)}{k} \right] \tag{II.A.12}\label{II.A.12} $

\[ f_{22} = \cfrac{1}{EI} \int_0^L x^2 \sqrt{ 1+ \left( \cfrac{\pi e_0}{L} \right)^2 \cos^2\left( \cfrac{\pi x}{L} \right) } \,dx = \cfrac{L^3}{\pi^3EI} \int_0^\pi \theta^2 \sqrt{ 1+k^2\cos^2\theta
} \,d\theta. \tag{II.A.13}\label{II.A.13} \]

$ f_{23} = \cfrac{1}{EI} \int_0^L x \sqrt{ 1+ \left( \cfrac{\pi e_0}{L} \right)^2 \cos^2\left( \cfrac{\pi x}{L} \right) } \, dx = \cfrac{L}{2}\,f_{33} = \cfrac{L^2}{\pi EI} \sqrt{1+k^2}\, E(m) \tag{II.A.14}\label{II.A.14} $

$ f_{33} = \cfrac{1}{EI} \int_0^L \sqrt{ 1+k^2\cos^2\left(\cfrac{\pi x}{L}\right) } \,dx = \cfrac{2L}{\pi EI} \sqrt{1+k^2}\;E(m) \tag{II.A.15}\label{II.A.15} $

gdzie wprowadzono oznaczenia:

• pomocnicze parametery imperferkcji 

\[ k = \cfrac{\pi e_0}{L} = \pi\varepsilon. \tag{II.A.16}\label{II.A.16} \]

\[ m=\cfrac{k^2}{1+k^2} \tag{II.A.17}\label{II.A.17} \]

• zupełna całka eliptyczna pierwszego rodzaju, 

\[ K(m) = \int_0^{\pi/2} \cfrac{d\theta} {\sqrt{1-m\sin^2\theta}}\tag{II.A.18}\label{II.A.18} \] 

• zupełna całka eliptyczna drugiego rodzaju.

\[ E(m) = \int_0^{\pi/2} \sqrt{1-m\sin^2\theta}\,d\theta \tag{II.A.19}\label{II.A.19} \] 

Współczynniki (\ref{II.A.10})–(\ref{II.A.15}) stanowią ścisłe przedstawienie macierzy podatności elementu imperfekcyjnego. Dla współczynników $f_{11}, f_{13}, f_{23}, f_{33}$  uzyskano postacie zamknięte, natomiast współczynniki $f_{12}, f_{22}$​ pozostają w postaci całkowej.

W granicy $e_0 \rightarrow 0$  otrzymujemy $ds \rightarrow dx $,  a współczynniki macierzy podatności redukują się do klasycznych zależności teorii Bernoulliego-Eulera.

Analiza asymptotyczna współczynników pozostających w postaci całkowej

W poprzednim punkcie uzyskano ścisłe postacie współczynników $f_{11}$, $f_{13}$, $f_{23}$ oraz $f_{33}$.  Współczynniki $f_{11}$, $f_{23}$ i $f_{33}$ wyrażają się za pomocą zupełnych całek eliptycznych, natomiast współczynnik $f_{13}$ można zapisać za pomocą funkcji elementarnych.
Dla współczynników $f_{12}$ oraz $f_{22}$ nie udało się natomiast otrzymać analogicznych postaci zamkniętych. Wynika to z obecności dodatkowych czynników wagowych $x$ oraz $x^2$, które naruszają strukturę klasycznych całek eliptycznych pojawiających się w rozwiązaniu Livesleya.

Współczynnikiwpostaci odpowiednio (\ref{II.A.11}) i (\ref{II.A.13})  można rozwinąć względem małego parametru imperfekcji zdefiniowanego następująco:

\[ \varepsilon=\cfrac{e_0}{L} \tag{II.A.20}\label{II.A.20}\]

 Rozwinięcie współczynników pozostających w postaci całkowej względem małego parametru jest szczególnie uzasadnione z punktu widzenia zastosowań inżynierskich, ponieważ w praktyce amplituda imperfekcji stanowi zwykle niewielki ułamek długości elementu. Pozwala to zastąpić złożone całki szeregiem kolejnych poprawek do rozwiązania klasycznej belki Bernoulliego–Eulera oraz ocenić wpływ imperfekcji geometrycznych na poszczególne składniki macierzy podatności.

Wprowadzamy ponadto bezwymiarową współrzędną $ \xi=\cfrac{x}{L}, \qquad 0\le\xi \le 1,$

Współczynniki pozostające w postaci całkowej przyjmują wówczas postać

\[  f_{12} = -\cfrac{e_0L^2}{EI} \int_0^1 \xi \sin(\pi\xi) \sqrt{ 1+k^2\cos^2(\pi\xi) } \,d\xi \tag{II.A.21}\label{II.A.21} \]

\[ f_{22} = \cfrac{L^3}{EI} \int_0^1 \xi^2 \sqrt{ 1+k^2\cos^2(\pi\xi) } \,d\xi. \tag{II.A.22}\label{II.A.22} \]

Dla niewielkich imperfekcji geometrycznych zachodzi

\[  k=\pi\varepsilon\ll 1. \tag{II.A.23}\label{II.A.23} \]

Rozwijając funkcję podpierwiastkową w szereg Taylora względem parametru $k$ otrzymujemy

\[ \sqrt{ 1+k^2\cos^2(\pi\xi) } = 1 + \cfrac{k^2}{2}\cos^2(\pi\xi) – \cfrac{k^4}{8}\cos^4(\pi\xi) +\ldots \tag{II.A.24}\label{II.A.24} \]

Podstawiając rozwinięcie (\ref{II.A.24}) do wzoru (\ref{II.A.21}) otrzymujemy

\[  f_{12} = -\cfrac{e_0L^2}{EI} \int_0^1 \xi\sin(\pi\xi)\,d\xi – \cfrac{e_0L^2k^2}{2EI} \int_0^1 \xi\sin(\pi\xi)\cos^2(\pi\xi)\,d\xi +\ldots \tag{II.A.25}\label{II.A.25} \]

\[  f_{22} = \cfrac{L^3}{EI} \int_0^1 \xi^2\,d\xi + \cfrac{L^3k^2}{2EI} \int_0^1 \xi^2\cos^2(\pi\xi)\,d\xi – \cfrac{L^3k^4}{8EI} \int_0^1 \xi^2\cos^4(\pi\xi)\,d\xi +\ldots \tag{II.A.26}\label{II.A.26} \]

Pierwsze całki występujące w rozwinięciach mają postać

\[  \int_0^1 \xi\sin(\pi\xi)\,d\xi = \cfrac{1}{\pi} \tag{II.A.27} \label{II.A.27}  \] 

Po podstawieniu zależności (II.A.27)–(II.A.28) do wzorów (II.A.25)–(II.A.26) otrzymujemy wiodące wyrazy rozwinięć asymptotycznych

\[  \int_0^1 \xi^2\,d\xi = \cfrac{1}{3}. \tag{II.A.28}\label{II.A.28} \]

\[  f_{12} = -\cfrac{e_0L^2}{\pi EI} + O(k^2) \tag{II.A.29}\label{II.A.29}\]

\[  f_{22} = \cfrac{L^3}{3EI} + O(k^2). \tag{II.A.30}\label{II.A.30} \]

 Otrzymane całki zawierają już wyłącznie funkcje elementarne i mogą zostać obliczone analitycznie. Pierwszy wyraz rozwinięcia współczynnika $f_{22}$ odpowiada klasycznej podatności belki Bernoulliego–Eulera, natomiast kolejne wyrazy opisują wpływ początkowej krzywizny osi elementu. Analogicznie współczynnik $f_{12}$ stanowi miarę sprzężenia pomiędzy oddziaływaniami osiowymi i poprzecznymi wywołanego obecnością imperfekcji geometrycznej.

Rozwinięcia współczynników macierzy podatności względem małego parametru imperfekcji

W celu dalszej analizy wpływu imperfekcji geometrycznych na właściwości elementu rozwijamy pozostałe oprócz $f_{12}$ (\ref{II.A.29}) oraz $f_{22}$ (\ref{II.A.30}) współczynniki macierzy podatności względem małego parametru (\ref{II.A.20}) .
Struktura rozwinięć wynika z symetrii zagadnienia względem zmiany znaku imperfekcji \(e_0\rightarrow -e_0\). Współczynniki diagonalne oraz współczynnik \(f_{23}\) są funkcjami parzystymi amplitudy imperfekcji, natomiast współczynniki sprzęgające \(f_{12}\) oraz \(f_{13}\) są funkcjami nieparzystymi. Dla współczynników diagonalnych otrzymujemy rozwinięcia zawierające wyłącznie parzyste potęgi parametru imperfekcji. 

Składowa $f_{11}$

\[  f_{11} = f_{11}^{(0)} + f_{11}^{(2)}\varepsilon^2 + f_{11}^{(4)}\varepsilon^4 +\ldots = \cfrac{L}{EA} + \cfrac{L^3}{2EI}\varepsilon^2 + \cfrac{\pi^2L^3}{16EI}\varepsilon^4 + O(\varepsilon^6).  \tag{II.A.31}\label{II.A.31} \]

Przy wyprowadzeniu składowej $f_{11}$ skorzystano  z rozwinięcia  (\ref{II.A.24}) oraz z tożsamosci matematycznych:

$ \sqrt{1+\pi^2\varepsilon^2\cos^2(\pi\xi)} = 1 + \cfrac{\pi^2\varepsilon^2}{2}\cos^2(\pi\xi) – \cfrac{\pi^4\varepsilon^4}{8}\cos^4(\pi\xi) + O(\varepsilon^6)$

$ \int_0^1\sin^2(\pi\xi)\,d\xi = \cfrac12$

$ \int_0^1\sin^2(\pi\xi)\cos^2(\pi\xi)\,d\xi = \cfrac18 $,

Składowa $f_{12}$

\[ f_{12} = -\cfrac{e_0L^2}{EI} \left[ \int_0^1 \xi\sin(\pi\xi)\,d\xi + \cfrac{k^2}{2} \int_0^1 \xi\sin(\pi\xi)\cos^2(\pi\xi)\,d\xi \right] + O(k^4) = -\cfrac{L^3}{\pi EI}\varepsilon -\cfrac{5\pi L^3}{18EI}\varepsilon^3 + O(\varepsilon^5).   \tag{II.A.32}\label{II.A.32} \]

Przy wyprowadzeniu składowej $f_{12}$ skorzystano  z rozwinięcia  (\ref{II.A.24}) podstawienia (\ref{II.A.16}) oraz z własnosći całki:

$\int_0^1 \xi\sin(\pi\xi) \cos^2(\pi\xi)\,d\xi = \cfrac{5}{9\pi}$ 

Składowa $f_{13}$

Współczynniki sprzęgające oddziaływania osiowe i zginające przyjuje postać

\[ f_{13} = f_{13}^{(1)}\varepsilon + f_{13}^{(3)}\varepsilon^3 +\ldots \tag{II.A.33}\label{II.A.33} \]

\[ f_{13} = -\cfrac{\varepsilon L^2}{EI} \int_0^1 \sin(\pi\xi)\,d\xi – \cfrac{\pi^2\varepsilon^3L^2}{2EI} \int_0^1 \sin(\pi\xi)\cos^2(\pi\xi)\,d\xi + O(\varepsilon^5) =-\cfrac{2L^2}{\pi  I}\varepsilon -\cfrac{\pi L^2}{3EI}\varepsilon^3 + O(\varepsilon^5). \tag{II.A.34}\label{II.A.34} \]

Składowa $f_{22}$

\[f_{22} = \cfrac{L^3}{EI} \left[ \int_0^1 \xi^2\,d\xi + \cfrac{k^2}{2} \int_0^1 \xi^2\cos^2(\pi\xi)\,d\xi – \cfrac{k^4}{8} \int_0^1 \xi^2\cos^4(\pi\xi)\,d\xi \right] + O(k^6) = \cfrac{L^3}{3EI} + \cfrac{L^3}{EI} \left( \cfrac{\pi^2}{12} + \cfrac18 \right)\varepsilon^2 – \cfrac{L^3}{EI} \left( \cfrac{\pi^4}{64} + \cfrac{15\pi^2}{512} \right)\varepsilon^4 + O(\varepsilon^6).\tag{II.A.35}\label{II.A.35} \]

przy czym wykorzystano zależnosći: $\int_0^1 \xi^2\cos^2(\pi\xi)\,d\xi = \cfrac16+\cfrac{1}{4\pi^2}$ , $ \int_0^1 \xi^2\cos^4(\pi\xi)\,d\xi = \cfrac18+\cfrac{15}{64\pi^2}$. 

Zatem

\[ f_{22}^{(0)} = \cfrac{L^3}{3EI}, \qquad f_{22}^{(2)} = \cfrac{L^3}{EI} \left( \cfrac{\pi^2}{12} + \cfrac18 \right), \qquad f_{22}^{(4)} = -\cfrac{L^3}{EI} \left( \cfrac{\pi^4}{64} + \cfrac{15\pi^2}{512} \right). \tag{II.A.36}\label{II.A.36} \]

Składowa $f_{23}$

współczynnik sprzęgający zginanie i obrót przyjmuje postać

\[ f_{23} = f_{23}^{(0)} + f_{23}^{(2)}\varepsilon^2 + f_{23}^{(4)}\varepsilon^4 +\ldots \tag{II.A.37}\label{II.A.37} \]

Składowa $f_{33}$

\[ f_{33} = f_{33}^{(0)} + f_{33}^{(2)}\varepsilon^2 + f_{33}^{(4)}\varepsilon^4 +\ldots \tag{II.A.38}\label{II.A.38} \]

Dla współczynników \(f_{23}\) i \(f_{33}\) wygodniej jest wykorzystać wcześniej uzyskane postacie zawierające zupełne całki eliptyczne i rozwinąć funkcje \(K(m)\) (\ref{II.A.18}) oraz \(E(m)\) (\ref{II.A.19}) względem małego parametru \(m\).

\[ K(m) = \cfrac{\pi}{2} \left( 1+\cfrac{m}{4} +\cfrac{9m^2}{64} +O(m^3) \right), \qquad E(m) = \cfrac{\pi}{2} \left( 1-\cfrac{m}{4} -\cfrac{3m^2}{64} +O(m^3) \right). \tag{II.A.39}\label{II.A.39} \]

Po podstawieniu powyższych rozwinięć do zależności (\ref{II.A.14})–(\ref{II.A.15}) otrzymujemy

\[ f_{23} = \cfrac{L^2}{2EI} +\cfrac{\pi^2L^2}{8EI}\varepsilon^2 -\cfrac{3\pi^4L^2}{128EI}\varepsilon^4 +O(\varepsilon^6). \tag{II.A.40}\label{II.A.40} \]

\[ f_{33} = \cfrac{L}{2EI} -\cfrac{\pi^2L}{8EI}\varepsilon^2 +\cfrac{9\pi^4L}{128EI}\varepsilon^4 + O(\varepsilon^6). \tag{II.A.41}\label{II.A.41} \]

Otrzymane wyniki pokazują, że wpływ imperfekcji geometrycznej pojawia się już w drugim rzędzie rozwinięcia. Współczynniki diagonalne pozostają dodatnie, podobnie jak współczynnik sprzęgający \(f_{23}\). Wszystkie poprawki zależą od parzystych potęg parametru \(\varepsilon\), co jest konsekwencją symetrii zagadnienia względem zmiany zwrotu początkowej imperfekcji osi elementu.

Wyznacznik $D_F$ macierzy  podatności $\mathbf{F}$

Wyznacznik macierzy podatności

\[ D_F=\det[F]. \tag{II.A.42}\label{II.A.42} \]

decyduje o odwracalności macierzy podatności oraz o właściwościach odpowiadającej jej macierzy sztywności. Dla rozpatrywanej macierzy symetrycznej  (\ref{II.A.9}) wyznacznik można zapisać w postaci

\[ D_F= f_{11}f_{22}f_{33} + 2f_{12}f_{13}f_{23} – f_{11}f_{23}^{\,2} – f_{22}f_{13}^{\,2} – f_{33}f_{12}^{\,2}. \tag{II.A.43}\label{II.A.43} \]

Ponieważ współczynniki \(f_{12}\) oraz \(f_{13}\) są funkcjami nieparzystymi parametru imperfekcji, natomiast pozostałe współczynniki są funkcjami parzystymi, wyznacznik również jest funkcją parzystą

\[ D(\varepsilon)=D(-\varepsilon). \tag{II.A.44}\label{II.A.44} \]

Można więc przedstawić go w postaci szeregu

\[ D_F= D_F^{(0)} + D_F^{(2)}\varepsilon^2 + D_F^{(4)}\varepsilon^4 +\ldots \tag{II.A.45}\label{II.A.45} \]

W zerowym przybliżeniu odpowiadającym idealnie prostemu prętowi otrzymujemy

\[ D_F^{(0)} = f_{11}^{(0)} \left( f_{22}^{(0)}f_{33}^{(0)} – \left(f_{23}^{(0)}\right)^2 \right). \tag{II.A.46}\label{II.A.46} \]

Po podstawieniu wartości zerowego rzędu

\[ f_{11}^{(0)}=\cfrac{L}{EA}, \qquad f_{22}^{(0)}=\cfrac{L^3}{3EI}, \qquad f_{23}^{(0)}=\cfrac{L^2}{2EI}, \qquad f_{33}^{(0)}=\cfrac{L}{EI} \]

otrzymujemy

\[ D_F^{(0)} = \cfrac{L}{EA} \left( \cfrac{L^4}{3(EI)^2} – \cfrac{L^4}{4(EI)^2} \right) = \cfrac{L^5}{12\,EA\,(EI)^2}. \tag{II.A.47}\label{II.A.47} \]

Wyznacznik zerowego rzędu jest dodatni, co potwierdza dodatnią określoność macierzy podatności odpowiadającej klasycznej belce Bernoulliego–Eulera. Oznacza to również istnienie macierzy sztywności będącej odwrotnością macierzy podatności.

Pierwsza poprawka związana z imperfekcją pojawia się w rzędzie \(\varepsilon^2\). Jej wartość zależy od wszystkich współczynników drugiego rzędu oraz od iloczynów współczynników nieparzystych \(f_{12}\) i \(f_{13}\). W szczególności człony

\[ -f_{22}f_{13}^{\,2} \qquad\text{oraz}\qquad -f_{33}f_{12}^{\,2} \tag{II.A.48}\label{II.A.48} \]

powodują zmniejszenie wyznacznika, a więc prowadzą do obniżenia efektywnej sztywności układu. Jest to bezpośrednia konsekwencja pojawienia się sprzężeń pomiędzy oddziaływaniami osiowymi i zginającymi wywołanych przez początkową krzywiznę osi elementu.

Dla typowych wartości imperfekcji spotykanych w praktyce inżynierskiej, rzędu \(e_0/L \approx 1/200\), parametr rozwinięcia wynosi \(\varepsilon \approx 5\cdot10^{-3}\). Oznacza to, że składniki proporcjonalne do \(\varepsilon^4\) są około cztery rzędy wielkości mniejsze od składników drugiego rzędu. W konsekwencji w większości zastosowań wystarczające jest ograniczenie analizy do wyrazów proporcjonalnych do \(\varepsilon^2\)

W celu wyznaczenia kolejnych wyrazów rozwinięcia podstawiamy rozwinięcia współczynników macierzy podatności do wzoru (\ref{II.A.43})  i grupujemy składniki względem kolejnych potęg  parametru

\[ D_F= D_F^{(0)} + D_F^{(2)}\varepsilon^2 + D_F^{(4)}\varepsilon^4 + O(\varepsilon^6). \tag{II.A.49}\label{II.A.49} \]

Współczynnik drugiego rzędu otrzymujemy przez zebranie wszystkich składników proporcjonalnych do \(\varepsilon^2\)

\[ \begin{aligned} D_F^{(2)}={}& f_{11}^{(2)}f_{22}^{(0)}f_{33}^{(0)} + f_{11}^{(0)}f_{22}^{(2)}f_{33}^{(0)} + f_{11}^{(0)}f_{22}^{(0)}f_{33}^{(2)} \\
& + 2f_{12}^{(1)}f_{13}^{(1)}f_{23}^{(0)} – 2f_{11}^{(0)}f_{23}^{(0)}f_{23}^{(2)} \\
& – f_{22}^{(0)}\!\left(f_{13}^{(1)}\right)^2 – f_{33}^{(0)}\!\left(f_{12}^{(1)}\right)^2 – f_{23}^{(0)2}f_{11}^{(2)}.
\end{aligned} \tag{II.A.50}\label{II.A.50} \]

W analogiczny sposób współczynnik czwartego rzędu przyjmuje postać

\[ D_F^{(4)} = D_{F,a}^{(4)} + D_{F,b}^{(4)} + D_{F,c}^{(4)}. \tag{II.A.51}\label{II.A.51} \]

gdzie pierwszy składnik obejmuje iloczyny współczynników czwartego rzędu

\[ \begin{aligned} D^{(4)}_{(a)}={}& f_{11}^{(4)}f_{22}^{(0)}f_{33}^{(0)} + f_{11}^{(0)}f_{22}^{(4)}f_{33}^{(0)} + f_{11}^{(0)}f_{22}^{(0)}f_{33}^{(4)} \\
& + f_{11}^{(2)}f_{22}^{(2)}f_{33}^{(0)} + f_{11}^{(2)}f_{22}^{(0)}f_{33}^{(2)} + f_{11}^{(0)}f_{22}^{(2)}f_{33}^{(2)}. \end{aligned} \tag{II.A.52}\label{II.A.52} \]

Wyznacznik $D_K$ macierzy sztywności $\mathbf{K}$

Znajomość rozwinięcia wyznacznika macierzy podatności umożliwia bezpośrednie wyznaczenie rozwinięcia wyznacznika macierzy sztywności. Ponieważ macierz sztywności jest odwrotnością macierzy podatności

\[ [K]=[F]^{-1}. \tag{II.A.53}\label{II.A.53} \]

Zachodzi zależność

\[ D_K = \det[K] = \cfrac{1}{\det[F]} = \cfrac{1}{D_F}. \tag{II.A.54}\label{II.A.54} \]

Podstawiając rozwinięcie asymptotyczne wyznacznika podatności (\ref{II.A.45}) i rozwijając odwrotność w szereg względem parametru \(\varepsilon\), otrzymujemy

\[ D_K = \cfrac{1}{D_F^{(0)}} \left[ 1 – \cfrac{D_F^{(2)}}{D_F^{(0)}}\varepsilon^2 + \left( \cfrac{\left(D_F^{(2)}\right)^2}{\left(D_F^{(0)}\right)^2} – \cfrac{D_F^{(4)}}{D_F^{(0)}} \right)\varepsilon^4 + O(\varepsilon^6) \right].\tag{II.A.55}\label{II.A.55} \]

Otrzymana zależność opisuje wpływ amplitudy imperfekcji na wyznacznik macierzy sztywności. W praktyce inżynierskiej dominujące znaczenie ma poprawka drugiego rzędu, natomiast składniki wyższych rzędów stają się istotne dopiero przy dużych imperfekcjach lub analizach wymagających podwyższonej dokładności.

Macierz sztywności elementu 

Macierz sztywności wyznacza się przez odwrócenie macierzy podatności (\ref{II.A.53}). Ponieważ współczynniki macierzy podatności zostały wyznaczone w postaci ścisłej, możliwe jest również uzyskanie ścisłej postaci odpowiadającej macierzy sztywności.

W granicy $ e_0 \rightarrow 0$  oś środkowa elementu przechodzi w linię prostą, a macierz sztywności redukuje się do klasycznej macierzy elementu belkowego Bernoulliego–Eulera. Dla niezerowej amplitudy imperfekcji w macierzy sztywności pojawiają się dodatkowe wyrazy sprzęgające oddziaływania osiowe i zginające. Są one bezpośrednią konsekwencją początkowej krzywizny osi środkowej i stanowią podstawową różnicę pomiędzy elementem imperfekcyjnym Livesleya a klasycznym elementem prostoliniowym.

Wnioski ogólne

Przeprowadzona analiza wykazała, że dla elementu imperfekcyjnego o osi opisanej pojedynczą półfalą sinusoidalną większość współczynników macierzy podatności można przedstawić w postaci zamkniętej. Współczynniki $f_{11}$, $f_{13}$, $f_{23}$ oraz $f_{33}$ wyrażają się za pomocą funkcji elementarnych oraz zupełnych całek eliptycznych pierwszego i drugiego rodzaju. Jedynie współczynniki $f_{12}$ oraz $f_{22}$ nie prowadzą bezpośrednio do standardowych funkcji specjalnych. Przyczyną jest obecność dodatkowych czynników wagowych $x$ oraz $x^2$, które zmieniają strukturę odpowiednich całek i uniemożliwiają ich bezpośrednią redukcję do klasycznych całek eliptycznych.

Wprowadzenie bezwymiarowej współrzędnej $\xi=x/L$ oraz małego parametru imperfekcji $\varepsilon=e_0/L$ pozwala jednak na skuteczne rozwinięcie tych współczynników w szeregi asymptotyczne. Otrzymane rozwinięcia pokazują, że pierwszy wyraz współczynnika $f_{22}$ odpowiada klasycznej podatności belki Bernoulliego–Eulera, natomiast kolejne wyrazy reprezentują poprawki wynikające z początkowej krzywizny osi elementu. Współczynnik $f_{12}$ ma charakter sprzęgający i zanika wraz z amplitudą imperfekcji. Jego obecność stanowi bezpośrednią konsekwencję początkowej krzywizny osi środkowej i opisuje wzajemne oddziaływanie deformacji osiowych oraz poprzecznych.

Otrzymane wyniki potwierdzają, że klasyczna teoria belki prostoliniowej stanowi graniczny przypadek rozwiązania Livesleya odpowiadający przejściu $ \varepsilon=\cfrac{e_0}{L}\rightarrow0$ . W tym sensie model imperfekcyjny można interpretować jako naturalne rozszerzenie teorii Bernoulliego–Eulera na elementy posiadające geometrycznie niedoskonałą oś środkową.

Znaczenie dla współczynnika redukcyjnego (wvboczeniowego) Perry’ego

Przeprowadzona analiza pokazuje, że współczynnik redukcyjny Perry’ego stanowi w istocie skondensowany opis skutków geometrycznych wynikających z początkowej imperfekcji osi środkowej pręta. W klasycznej teorii wyboczenia wpływ imperfekcji sprowadzany jest do pojedynczego parametru zastępczego, podczas gdy model Livesleya pozwala prześledzić mechanizm jego powstawania na poziomie macierzy podatności elementu. Iimperfekcja geometryczna prowadzi do dwóch odmiennych efektów:
1) pojawienie się dodatkowych sprzężeń pomiędzy oddziaływaniami osiowymi i poprzecznymi, reprezentowanych przez współczynniki mieszane macierzy podatności.
2) modyfikacja samych podatności głównych, a więc efektywnej sztywności elementu.

W teorii Perry’ego oba efekty zostają zastąpione pojedynczym współczynnikiem redukcyjnym zależnym od smukłości i przyjętej amplitudy imperfekcji. Takie podejście jest uzasadnione dopóty, dopóki dominujący wpływ na zachowanie pręta wywiera pierwsza postać wyboczeniowa, a amplituda imperfekcji pozostaje niewielka w porównaniu z długością elementu. Analiza Livesleya pokazuje jednak, że wraz ze wzrostem imperfekcji coraz większe znaczenie uzyskują wyrazy wyższych rzędów rozwinięcia asymptotycznego. Współczynnik Perry’ego nie rozróżnia źródeł tych efektów i sprowadza je do jednej liczby redukcyjnej. W rezultacie dokładność teorii Perry’ego maleje wraz ze wzrostem amplitudy imperfekcji oraz wtedy, gdy rzeczywisty kształt osi odbiega od pojedynczej półfali sinusoidalnej.

Szczególnie istotne jest to, że teoria Perry’ego zachowuje wyłącznie informację o maksymalnej amplitudzie imperfekcji, natomiast nie uwzględnia szczegółowej geometrii jej przebiegu. Tymczasem analiza energetyczna pokazuje, że współczynniki macierzy podatności zależą od całego rozkładu krzywizny osi środkowej. Dwa pręty posiadające tę samą wartość maksymalnego odchylenia mogą zatem wykazywać różne właściwości podatnościowe i różną nośność graniczną.

Oznacza to, że teoria Perry’ego jest najbardziej dokładna dla niewielkich imperfekcji oraz dla przypadków zdominowanych przez pierwszą postać wyboczeniową. W miarę wzrostu amplitudy imperfekcji lub pojawiania się bardziej złożonych kształtów początkowej osi środkowej przewagę uzyskują modele energetyczne i imperfekcyjne, w których geometria elementu jest uwzględniana bezpośrednio. Z tego punktu widzenia współczynnik Perry’ego można interpretować jako pierwszy wyraz rozwinięcia pełnego modelu imperfekcyjnego. Model Livesleya stanowi natomiast jego naturalne uogólnienie, zachowujące pełną informację o geometrii elementu i pozwalające ocenić zakres stosowalności klasycznej teorii wyboczeniowej.  

Literatura

1. R. K. Livesley, D. B. Chandler, <i>Stability Functions for Structural Frameworks</i>, Manchester University Press, Manchester, 1956
2. Przemieniecki J. S., <i>Theory of Matrix Structural Analysis</i>, McGraw–Hill Book Company, New York, 1968
3. Euler (1744), Methodus Inveniendi Lineas Curvas Maximi Minive Proprietate Gaudentes. Lausanne and Geneva: Marc-Michel Bousquet
4. Wagner, H. (1931), Flat Sheet Metal Girders with Very Thin Metal Webs. Part I: General Theories and Assumptions. NACA Technical Memorandum TM-604
5. von Karman, T., Tsien, H. S. (1941), The Buckling of Thin Cylindrical Shells under Axial Compression, Journal of the Aeronautical Sciences, 8(8), 303–312
6. Koiter, W. T. (1945), Over de Stabiliteit van het Elastisch Evenwicht, Doctoral Dissertation, Delft University of Technology
7. PN-EN-1993-8, Projektowanie konstrukcji stalowych. Część 1-8: Projektowanie węzłów i połączeń
8. PN-EN-1993-8, Projektowanie konstrukcji stalowych, Węzły i połączenia
9. Piechnik, St., Wytrzymałość materiałów dla wydziałów budowlanych. PWN, Warszawa-Kraków, 1980
10. Livesley R.K. (1975), Matrix Methods of Structural Analysis, 2nd Edition, Pergamon Press, Oxford

________________________________

Related Baza wiedzy

• Belki stalowe w podejściu normowym
• Stalowe dźwigary kratowe (+)