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 34 Czytelników
CZĘŚC II: NIELINIOWOŚĆ geometryczna i podstawy MES
Wprowadzenie
W teorii imperfekcyjnej wprowadzono mnożnik obciążenia konstrukcji $\Lambda$, definiowany 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 \]
Wartość krytyczna mnożnika $\Lambda_{cr}$ odpowiada stanowi utraty stateczności konstrukcji. Wygodną miarą stopnia wykorzystania nośności wyboczeniowej jest bezwymiarowa wielkość
\[ \bar{\Lambda}=\frac{\Lambda}{\Lambda_{cr}} \tag{II.1} \label{II.1} \]
która osiąga wartość równą jedności w stanie krytycznym.
Uogólnienie teorii belki Bernoulli-Euler o ściskanie. Nieliniowości geometryczne
Teoria zginania Bernoulli-Euler może zostać rozszerzona na przypadek elementów poddanych jednocześnie 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 postaci analitycznej.
Model elementu $[e]$ przedstawiono na rys. II.1. Siła przekrojowa jest stała na całej długości elementu i wynosi $ N=F$. Obciążenie krytyczne konstrukcji sprowadza się w tym przypadku do klasycznego obciążenia krytycznego Eulera
\[ F_{cr}=N_{cr} = \frac{\pi^2EI}{L^2} \tag{II.2} \label{II.2} \]
Stopień wykorzystania nośności wyboczeniowej (\ref{II.1}) można więc w tym przypadku zapisać w postaci
\[ \bar{\Lambda} = \frac{F}{F_{cr}} = \frac{N}{N_{cr}}\tag{II.3} \label{II.3} \]
Dalsze rozważania prowadzone są dla tego szczególnego przypadku, przy zachowaniu notacji zgodnej z ogólną teorią stateczności oraz imperfekcji konstrukcji.
Na rys. II.1a przedstawiono model elementu belki-słupa $[e]$, na rys. II.1b belkę obciążoną poprzecznie i ściskaną, natomiast na rys. II.1c słup wspornikowy zginany i ściskany. Układy takie prowadzą do zagadnień geometrycznie nieliniowych związanych z efektem drugiego rzędu oraz utratą stateczności.
Rys. II.1 . Belka-słup Bernoulli zginana poprzecznie i ściskana: a) model elementu [e]; b) belka-słup poziomy; c) słup-belka pionowy
Ścisła macierz sztywności elementu belkowi-słupa
Nieliniową macierz sztywności ściskanego elementu prętowego [e], pokazanego na rys. 6a, wyznaczymy przez rozwiązanie równania różniczkowego belki-słupa bez obciążenia poprzecznego $q_z=0$ (I.33), uzupełnionego o wyraz geometrycznej nieliniowości wynikający z działania osiowej siły ściskającej $P$.
\[ 0= EI_y\cfrac{d^4w}{dx^4} + N\cfrac{d^2w}{dx^2} \tag{II.4} \label{II.4} \]
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} = \frac{N}{N_{cr}}$ oraz $ N_{cr} =\frac{\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.5} \label{II.5}\]
W granicznym stanie krytycznym, dla $\bar{\Lambda}=1$ oraz $\mu=\pi$, otrzymujemy p o podstawieniu powyższych zależności do równania (II.4) otrzymujemy jego postać bezwymiarową
\[ \cfrac{d^4w}{d\xi^4} + \mu^2 \cfrac{d^2w}{d\xi^2} = 0 \tag{II.6} \label{II.6} \]
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.7} \label{II.7} \]
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.8} \label{II.8} \]
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. \]
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.9} \label{II.9} \]
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.10} \label{II.10} \]
Po odwróceniu macierzy $[D]$ otrzymujemy
\[ [D]^{-1} = \cfrac{1}{\Theta} \begin{bmatrix}
\mu\beta & \gamma & \mu(1-\cos\mu) & -(\mu-\sin\mu)\\
\mu\sin\mu & 1-\cos\mu & -\mu\sin\mu & 1-\cos\mu\\
\mu(1-\cos\mu) & -\gamma & -\mu(1-\cos\mu) & \mu-\sin\mu\\
-\mu\sin\mu & \beta & \mu\sin\mu & -(1-\cos\mu)
\end{bmatrix}
\tag{II.11} \label{II.11} \]
gdzie:
$ \Theta = 2\mu(1-\cos\mu) – \mu^2\sin\mu $,
$\beta = 1-\cos\mu-\sin\mu$,
$\gamma = \mu\cos\mu-\sin\mu $.
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.12} \label{II.12} \]
\[ 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.13} \label{II.13} \]
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.14} \label{II.14}\]
Współczynniki $\Psi_i$ są funkcjami statecznościowymi opisującymi wpływ siły osiowej na sztywność 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.15} \label{II.15} \]
gdzie:
$ \alpha_c = \alpha\operatorname{ctg}\alpha$
$\alpha = \cfrac{\mu}{2} = \cfrac{\pi}{2}\sqrt{\bar{\Lambda}}$
Dla rozważanego przypadku pojedynczego pręta zachodzi bowiem $ \bar{\Lambda} = \frac{N}{N_{cr}}.$
W przypadku znajomości siły osiowej $N$ funkcje statecznościowe (II.15 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 ściskanie $\mu^2$ w postaci
\[ \Psi_i \approx \sum_{j=0}^{4} C_{i,j} (\mu^2)^j \tag{II.16} \label{II.16} \]
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$ otrzymujemy
\[ \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.17} \label{II.17} \]
Wyrażenia (II.17) definiują klasyczną macierz geometryczną elementu prętowego stosowaną w analizie II rzędu.
W analizach stateczności konstrukcji bardziej użyteczne jest przedstawienie funkcji statecznościowych jako funkcji globalnego mnożnika obciążenia $\ bar{\Lambda}=\frac{\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 (II.16) 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–Δ. Człony wyższych rzędów prowadzą odpowiednio do analiz III, IV i kolejnych 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.18} \label{II.18} \]
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.19} \label{II.19} \]
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.
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=\frac{\pi}{2}\sqrt{\bar{\Lambda}}.$
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.20}\label{II.20} \]
Funkcje statecznościowe mogą być więc 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.
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=\frac{\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 j & 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) Pozostawienie wyłącznie wyrazu zerowego odpowiada klasycznej analizie I rzędu,
(2) Uwzględnienie wyrazu liniowego prowadzi do klasycznej macierzy geometrycznej stosowanej w analizie II rzędu (LBA, P–Δ),
(3) Kolejne wyrazy szeregu tworzą systematyczne przybliżenia ścisłych funkcji statecznościowych,
(4) Rozwinięcie jest poprawne w otoczeniu punktu $\bar{\Lambda}=0$,
(5) Dokładność rozwinięcia maleje wraz ze wzrostem parametru $\bar{\Lambda}$, szczególnie w pobliżu stanu krytycznego,
(6) W granicy $\bar{\Lambda}\rightarrow 1$ zaleca się stosowanie pełnych funkcji $\Psi_i$ lub ścisłej macierzy sztywności,
(7) Klasyczna macierz geometryczna jest jedynie pierwszym niezerowym wyrazem rozwinięcia ścisłej macierzy sztywności belki-słupa,
(8) Rozwinięcie Taylora tworzy formalny pomost pomiędzy analizą P–Δ a ścisłą teorią geometrycznie nieliniowej sztywności elementu,
(9) Ta sama interpretacja zachowuje ważność zarówno dla pojedynczego elementu, jak i dla wieloprętowych modeli MES, jeżeli parametr $\bar{\Lambda}$ pochodzi z globalnego zagadnienia stateczności konstrukcji.
Wniosek:
Rozwinięcie funkcji statecznościowych względem parametru $\bar{\Lambda}$ pozwala bezpośrednio powiązać kolejne rzędy teorii geometrycznie nieliniowej z klasycznymi procedurami MES stosowanymi w analizie stateczności konstrukcji. Parametr $\bar{\Lambda}$ posiada przy tym jednoznaczną interpretację fizyczną zarówno dla pojedynczego elementu, jak i całej konstrukcji, dzięki czemu stanowi naturalną zmienną do porównywania analiz I rzędu, analiz P–Δ oraz wyższych rzędów przybliżenia ścisłej macierzy sztywności.
Ugięcie belki ściskanej przed i w momencie utraty stateczności Eulera
Z zależności (\ref{II.7}), po podstawieniu uzyskanych wyrażeń na stałe całkowania, można wyznaczyć funkcję ugięcia pręta pomiędzy węzłami w zależności od przemieszczeń węzłowych
\[ w(\xi)= \cfrac{1}{\Theta} \begin{vmatrix}
1-\cos\mu+\cos(\mu\xi)-\cos[\mu(1-\xi)]-\mu\sin[\mu(1-\xi)]\\
-\mu\xi-\mu(1-\xi)\cos\mu+\xi\cos[\mu(1-\xi)]+\sin\mu+\sin(\mu\xi)\sin[\mu(1-\xi)]\\
1-\cos\mu-\cos(\mu\xi)+\cos[\mu(1-\xi)]-\mu\xi\sin\mu\\
-\mu(1-\xi)-\mu\xi\cos\mu+\mu\cos(\mu\xi)+\sin\mu-\sin(\mu\xi)-\sin[\mu(1-\xi)]
\end{vmatrix} \cdot
\begin{bmatrix}
\mu w_1 &
\varphi_1L &
\mu w_2 &
\varphi_2L
\end{bmatrix}
\tag{II.21}
\label{II.21} \]
Z zależności (\ref{II.21}) wynika, że pręt idealnie prosty ugina się pod wpływem osiowej siły ściskającej w każdym przypadku, gdy węzły doznają przemieszczeń liniowych lub obrotów wstępnych, czyli w każdej konstrukcji rzeczywistej. W granicznym stanie krytycznym, dla $\bar{\Lambda}=1$ oraz $\mu=\pi$, otrzymujemy, otrzymujemy
\[ w(\xi)= \cfrac{1}{2\pi} \left[ \Sigma w +\Sigma\varphi\,L(\xi-\tfrac{1}{2}) +\cfrac{\pi}{2} \left( \Sigma\varphi\,L-2\Delta w \right) \cos(\pi\xi) -\Delta\varphi\,L\sin(\pi\xi) \right] \tag{II.22}\label{II.22} \]
Maksymalne ugięcie pręta, czyli strzałka wygięcia, wystąpi w punkcie o współrzędnej
\[ \xi_{sup} = \pm \cfrac{1}{\mu} \operatorname{arccos}(\pm\Gamma) \tag{II.23} \label{II.23} \]
gdzie $ \Gamma=\cfrac{1}{\sqrt{(C_3/C_4)^2+1}} $.
Wyrażenie uzyskano przez przyrównanie do zera pierwszej pochodnej funkcji (\ref{II.21}). W rezultacie otrzymuje się cztery pierwiastki odpowiadające wszystkim kombinacjom znaków. Po odrzuceniu pierwiastków ujemnych jako fizycznie niemożliwych i po podstawieniu dodatniej współrzędnej (\ref{II.23}) do wzoru (\ref{II.21}) otrzymujemy strzałkę ugięcia pręta w zależności od warunków brzegowych wyrażonych stałymi całkowania $C_3$ i $C_4$. W przypadku symetrycznych warunków brzegowych strzałka ugięcia wystąpi w środku pręta i wyniesie
\[ f=w(\xi=0.5) = \cfrac{1}{2} \left[ \Sigma w – \cfrac{\operatorname{tg}(\mu/4)\,\Delta\varphi\,L}{\mu} \right] \tag{II.24} \label{II.24} \]
W granicznym stanie krytycznym, dla $ \bar{\Lambda}=1, \qquad \mu=\pi,$ strzałka ugięcia (\ref{II.24}) osiąga wartość
\[ f_{cr} = \cfrac{\Sigma w}{2} – \cfrac{\Delta\varphi\,L}{2\pi} \tag{II.25} \label{II.25} \]
Dla zerowej siły osiowej $N=0$ zachodzi $\mu=0$, zatem
\[ f_0= \lim_{\mu\to 0} (\ref{II.24}) = \frac{\Sigma w}{2} – \frac{\Delta\varphi\,L}{8}, \tag{II.26} \label{II.26} \]
co jest zgodne z rezultatem teorii I rzędu.
Wyrażenia (\ref{II.22}) i (\ref{II.25}) wskazują, że w granicznym stanie krytycznym Eulera:
(1) funkcja ugięcia ma postać opisaną zależnością (\ref{II.22}), czyli jest kombinacją funkcji liniowej, sinusoidy i kosinusoidy, a więc może zostać przedstawiona w postaci szeregu Fouriera,
(2) strzałka ugięcia pręta pozostaje wielkością skończoną i niezerową.
W przykładzie 1 rozważono trzy przypadki szczególne belek-słupów.
W dalszej części przechodzimy od rozwiązania równania różniczkowego belki-słupa do sformułowania energetycznego. Podejście takie stanowi podstawę metody Ritza oraz współczesnej metody elementów skończonych.
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 [1]:
\[ \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.27} \label{II.27} \]
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.28} \label{II.28} \]
Energię (\ref{II.28}) wyrazimy przez siły przekrojowe, a następnie przez przemieszczenia. Zachodzi
\[ \sigma_x = -\cfrac{N}{A} + \cfrac{M(x)}{I_y}\,z \tag{II.29} \label{II.29} \]
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.28}) 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.30} \label{II.30} \]
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.31} \label{II.31} \]
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.32} \label{II.32} \]
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.33} \label{II.33} \]
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.34} \label{II.34} \]
Własność (\ref{II.34}) 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.35} \label{II.35} \]
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.36}\label{II.36} \]
Po rozwinięciu w szereg Maclaurina ostatniego pierwiastka w równaniu (\ref{II.36}) 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.37} \label{II.37} \]
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.38} \label{II.38} \]
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.39} \label{II.39} \]
gdzie znak minus wynika z faktu, że siły odporu są przeciwnie skierowane do ugięcia.
Na rys. II.3 pokazano 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.39}) z zapisem obciążeń skupionych za pomocą funkcji Diraca (\ref{II.33}). 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.40} \label{II.40} \]
Rys. II.4 Reakcje od punktowej podpory sprężystej pionowej i obrotowej
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.41} \label{II.41} \]
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.42} \label{II.42} \]
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.43} \label{II.43} \]
Analogiczne wyrażenia dla skupionych podpór sprężystych otrzymujemy po podstawieniu zależności (\ref{II.40}) do wzoru (\ref{II.43}).
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 $N$, wynosi:
\[ \Delta(dx) = dx-dx’ = dx-\sqrt{(dx)^2-\left(\cfrac{d w}{d x}dx\right)^2} = dx\left[1-\sqrt{1-(w’)^2}\right] \tag{II.44}\label{II.44}
Po rozwinięciu w szereg Maclaurina ostatniego pierwiastka w ($\ref{61}$) i pozostawieniu pierwszego nieliniowego wyrazu szeregu otrzymujemy aproksymację, która jest 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.45}\label{II.45}
W rezultacie praca siły podłużnej ściskającej wynosi:
\[ L_N = \int\limits_0^L N\cdot\Delta(dx) = \int\limits_0^L \cfrac{N}{2}[w'(x)]^2\,dx \tag{II.46}\label{II.46}
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_B= \int\limits_0^L \left[ \cfrac{EI_y}{2}[w”(x)]^2 – \cfrac{N(x)}{2}[w'(x)]^2 + \cfrac{c^v(x)}{2}[w(x)]^2 – q(x)w(x) \right] dx + \sum_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.46}\label{II.47}
Poszukiwanie punktu stacjonarnego funkcjonału (\ref{II.47}) można przeprowadzić dowolną metodą wariacyjną, na przykład metodą Ritza, również w postaci jej implementacji w metodzie elementów skończonych. Przykład 2 przedstawia takie rozwiązanie dla belki Timoshenko.
Uogólnienie teorii belki na elementy krzywoliniowe
Funkcjonał (\ref{II.47}) 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ą.
Macierz podatności elementu z imperfekcją geometryczną
Przedstawione wcześniej sformułowanie energetyczne odnosi się do elementu o osi początkowo prostoliniowej. W praktyce inżynierskiej założenie takie stanowi jedynie idealizację, ponieważ rzeczywiste elementy konstrukcyjne zawsze wykazują pewne odchylenia od geometrii nominalnej.
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ą na podstawie energii odkształcenia oraz twierdzenia Castigliano.
Rys. II.5 Belka krzywoliniowa, b) sinusoidalna stanowiąca imperefekcję geometryczną (opracowano na podstawie pracy [2]
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.
Przykład 1 [ Nieliniowa geometrycznie belka wolnodparta]
Przeprowadzić analizę nieliniowej geometrycznie: belki wolnopodpartej, przedstawionej na rys. II.1.b z imperfekcjami osi w kształcie sinusoidy.
Warunki brzegowe
Warunki brzegowe dla belki wolnopodpartej , pokazanego na rys. 6b są następujące:
$ u_1=0, \quad w_1=0, \quad V _2= – C_{\Delta} \cdot u_2 $
gdzie: $C_\Delta$ jest stała sprężystości podparcia w węźle prawym (2).
Stopnie swobody węzłów belki ; $\begin{vmatrix} \varphi_1 \\u_2 \\ w_2 \\ \varphi_2 \end {vmatrix}$
Zmodyfikowana macierz sztywności
$ [ \widetilde K ] = \cfrac{EI}{L} \cdot \begin{bmatrix}
4 \Lambda_3 & 0 & \cfrac {-6 \Lambda_2} {L} & 2 \Lambda_4 \\
& \cfrac {EA}{EI} & 0 & 0 \\
&& \cfrac{12 \Lambda_1 + c_{\Delta}} {L^2} & \cfrac { – 6\Lambda_2}{L} \\
& SYM & \cfrac{ -6 \Lambda_2}{L} & 4 \Lambda_3 \\
\end{bmatrix}$
Odwrotność macierzy zmodyfikowanej
$[\widetilde K]^{-1} = \cfrac{EI}{L} \cdot \begin{bmatrix}
k_{11} & 0 & 3 L \cdot \Lambda_2 & k_{14} \\
0 & \cfrac {\Theta} {EA} & 0 & 0 \\
& & L ( 2 \Lambda_3+\Lambda_4 ) & 3 L \cdot \Lambda_2 \\
& SYM & & k_{11} \\
\end{bmatrix}$
gdzie:
$c_\Delta= C_\Delta \cdot \cfrac{L^3}{EI}$,
$k_{11}= \cfrac{-9\Lambda_2^2 +(c_\Delta + 12 \Lambda_1 ) \Lambda_3}{2 \Lambda_3 -\Lambda_4}$,
$k_{14}=\cfrac{18 \Lambda_2^2 – (c_\Delta + 12 \Lambda_1 ) \Lambda_4}{4 \Lambda_3 -2 \Lambda_4}$,
$\Theta=EI\cdot [36 \Lambda_2^2 -(12\Lambda_1 +c_\Delta) \cdot (2 \Lambda_3 +\Lambda_4)]$.
Siły węzłowe od obciążenia rozłożonego $q$
$ [F]_q = \cfrac{qL}{2} [ L/6 \quad 0 \quad 1 \quad – L/6 ]$
Z formuły ($\ref {50} $) można wyznaczyć ścisłą strzałkę ugięcia pręta (z dokładnością do odstępstwa o d symetrycznych warunków brzegowych).
$ f_N =w_N (\xi=0,5) = \cfrac{qL^4}{48EI} \cdot \left ( \cfrac{12 (2 \Lambda_3 + \Lambda_4)}{36 \Lambda_2^2 +(12 c_\Delta \Lambda_1) \cdot (2 \Lambda_3 +\Lambda_4)} +\cfrac{\alpha \cdot tg(\alpha/2)}{2 \Lambda_3 +\Lambda_4 } \right)$
Strzałka ugięcia $f_q$ pochodząca od sił poprzecznych jest całką szczególną zagadnienia, która powinna być dodana do powyższej. W przypadku obciążenia równomiernie rozłożonego całka szczególna wynosi:
Ponieważ
$w_q (\xi) = \cfrac {qL^4}{24 EI}$, to
$f_q = w(\xi=0,5) = \cfrac{q L^4}{384 EI}$
Przykład 2 [ Nieliniowy geometrycznie słup wspornikowy]
Przeprowadzić analizę nieliniowej geometrycznie: słupa wspornikowego, przedstawionego na rys. II.1.c z imperfekcjami osi w kształcie sinusoidy.
Ściskany pręt wspornikowy obciążony jest poprzecznie siłą H w wierzchołku (1) lub obciążeniem równomiernie rozłożonym h po wysokości. Pręt jest sprężyście utwierdzony w węźle dolnym (1) ze stałą sprężystości $c_\varphi$. (rys.II.1c)
Warunki brzegowe
$ u_2=0, \quad w_2=0, \quad M_2 = – C^{\varphi} \cdot \varphi_2$
Macierz zmodyfikowana
W przypadku konstrukcji jednoelementowej macierz ($ \ref{42}$) jest macierzą konstrukcji i można ją bezpośrednio modyfikować w celu uwzględnienia warunków brzegowych . Modyfikację macierzy sztywności dokonamy poprzez wykreślenie wierszy i kolumn, odpowiadających odebranym stopniom swobody i bezpośrednie wpisanie warunku sprężystej podpory do równania kanonicznego. W rezultacie otrzymamy macierz zmodyfikowaną $[\widetilde K]$:
$ [ \widetilde{K]} = \cfrac{EI}{L} \cdot
\begin{bmatrix}
\cfrac{EA}{EI} & 0 & 0 & 0 \\
& \cfrac{12 }{L^2} \Lambda_1 & \cfrac{6}{L} \Lambda_2 & \cfrac{6}{L} \Lambda_2 \\
& & 4 \Lambda_3 & 2 \Lambda_4 \\
& SYM & & 4 \Lambda_3+ c_{\varphi} \\
\end{bmatrix}
\begin{matrix}
u_2 \\ w_2 \\ \varphi_2 \\ \varphi_1
\end {matrix} $
Odwrócona macierz zmodyfikowana
Poprzez odwrócenie macierzy $[ \widetilde K]$ uzyskujemy macierz podatności nieliniowego geometrycznie wspornika:
$ [ \widetilde{K]}^{-1} = \cfrac{L}{\Theta}
\begin{bmatrix}
\cfrac{\Theta}{EA} & 0 & 0 & 0 \\
& k_{22} & k_{24}+\Lambda_2 \cdot c_{\varphi}\cdot L/2 & k_{24} \\
& & 3 \Lambda_2^2 -\Lambda_1(4 \Lambda_3+c_{\varphi}) & -3\Lambda_2^2 +2 \Lambda_1 \Lambda_4 \\
& SYM & & 3 \Lambda_2^2 -4 ] \Lambda_1 \Lambda_3 \\
\end{bmatrix}
\begin{matrix}
u_2 \\ w_2 \\ \varphi_2 \\ \varphi_1
\end {matrix} $
gdzie:
$c_{\varphi}=C_{\varphi} \cdot L/EI$,
$k_{24}=\Lambda_2L(2\Lambda_3-\Lambda_4)$,
$k_{22}=-\cfrac{L^3}{3}\cdot(4\Lambda_3^2-\Lambda_4^2+\Lambda_3c_{\varphi} ) $,
$\Theta = – 4EI \cdot [(2\Lambda_3 -\Lambda_4) \cdot \Lambda_1 (2 \Lambda_3 + \Lambda_4) – 3\Lambda_2^2] + c_{\varphi} \cdot L \cdot ( 3 \Lambda_2^2 – 4 \Lambda_1 \Lambda_3) $.
Przemieszczenia i siły przekrojowe
Przemieszczenia węzłowe uzyskuje się poprzez przemnożenie macierzy podatności $ [ \widetilde{K]}^{-1}$ przez wektor równoważników węzłowych, który wynosi:
od poziomego obciążenia rozłożonego h
$ [F]_h = \cfrac{hL}{2} \cdot [\quad 0 \quad 1 \quad L/6 \quad -L/6 \quad ] $
od obciążenia skupionego H
$ [F]_H = [ \quad 0 \quad H \quad 0 \quad 0 \quad ] $
Moment zginający w utwierdzeniu wspornika otrzymuje się po dodaniu do odpowiedniej współrzędnej wektora siły przywęzłowej, uzyskanej z przemnożenia wektora przemieszczeń węzłowych przez ostatni wiersz pierwotnej macierzy sztywności.
Przykład 3 [ Belka krzywoliniowa w kształcie sinusoidy]
Równanie osi elementu
Oś elementu jest opisana ogólną zależnością parametryczną
$$\begin{equation} x=x(s) \quad y=y(s) \label {110} \end {equation} $$
Macierz podatności elementu
Macierz podatności (odwrotność macierzy sztywności) elementu opisanego równaniem paramteryczym )$\ref{110}$) , otrzymana po przprowadzeniu algorytmu opisanego w tekście wynosi
$$\begin{equation} [ K]^{-1} =\begin{bmatrix}
\int \limits_0^L \left ( \cfrac{(dx/dy)^2}{EA}+\cfrac{y^2}{EI} \right ) ds & \int \limits_0^L \left ( \cfrac{(dx/dy) (dy/ds)} {EA}- \cfrac{xy}{EI} \right ) ds & \int \limits_0^L \cfrac{y}{EI} ds \\
& \int \limits_0^L \left ( \cfrac{(dy/ds)^2} {EA}- \cfrac{x^2}{EI} \right ) ds & – \int \limits_0^L \cfrac{x}{EI} ds \\
& SYM & \int \limits_0^L \cfrac{1}{EI} ds \\
\end{bmatrix} \label {111} \end {equation} $$
gdzie: $A(s)$ i $I(s)$ – pole przekroju i moment bezwładności sa znanymi funkcjami współrzędnej bieżącej s
Wyrażenie ($\ref{111}$) jest zgodne z podanym w pracy [3].
Macierz podtności w standardowym układzie współrzędnych
Kształt elementu wstępnie wygiętego w łuk sinusoidy opisany w standardowym układzie współrzędnych (x,y)
$$\begin{equation} y= \eta \sin{(\omega x)} \label {112} \end {equation} $$
gdzie: $\eta$ – amplituda, $\omega= \cfrac{ n \ pi x}{L}$ – częstotliwość dla $n=1,2,3,…$ – liczba naturalna oznaczająca numer postaci własnej elementu,
Równanie krzywej ($\ref{112}$) można zapisać w postaci parametrycznej:
$$\begin{equation} y= \eta \sin {\alpha} ; \quad x =\cfrac{\alpha}{\omega} \label {113} \end {equation}$$
gdzie parametr $\alpha \in [ \quad \alpha_1= -\omega L/2 ; \quad \alpha_2= \omega L/2 \quad ]$,
Poszczególne wyrazy macierzy podatności $( \ref{111}$) wynoszą:
$k_{11}= \int \limits_0^L \left ( \cfrac{(dx/dy)^2}{EA}+\cfrac{y^2}{EI} \right ) ds=\int \limits_{\alpha_1}^{\alpha_2} \left ( \cfrac{(dx/dy)^2}{EA}+\cfrac{y^2}{EI} \right ) d \alpha= \cfrac{ \eta ^2 [ n \pi – \sin {(n \pi) ] } }{2EI}$,
$k_{12}= \int \limits_{\alpha_1}^{\alpha_2} \left ( \cfrac{(dx/dy) (dy/ds)} {EA}- \cfrac{xy}{EI} \right ) d \alpha = \cfrac{ \eta L [ n \pi \cdot \cos {(n \pi /2) } – 2 \sin{(n \pi / 2)]} }{n \pi \cdot EI }$,
$k_{22}= \int \limits_{\alpha_1}^{\alpha_2} \left ( \cfrac{(dy/ds)^2} {EA}- \cfrac{x^2}{EI} \right ) d \alpha = \cfrac{L}{12} \left ( \cfrac{6 \eta^2 [ n \pi +\sin{(n \pi) ] }}{EA}+\cfrac{L^2 \cdot n \pi}{EI}\right )$,
$k_{13}= \int \limits_{\alpha_1}^{\alpha_2} \cfrac{y}{EI} d \alpha =0$,
$k_{23}= \int \limits_{\alpha_1}^{\alpha_2} \cfrac{-x}{EI} d \alpha =0$,
$k_{33}= \int \limits_{\alpha_1}^{\alpha_2} \cfrac{1}{EI} d \alpha = \cfrac{n \pi}{EI}$,
Macierz sztywności w układzie (x,y)
Macierz sztywności uzyskamy przez odwrócenie macierzy podatności opisanej wyżej . Stosunkowo proste wyrażenie uzyskamy dla kształtu pręta proporcjonalnego do pierwszej postaci własnej, czyli dla n=1:
$$\begin{equation} [K] = \cfrac{1}{\Theta} \begin{bmatrix}
2 \pi^3 EI ( EA L^2/ \eta +6EI \cdot \eta & 48 \pi L \cdot EA \cdot EI & 0 \\
& 12 \pi^3 EA \cdot EI \cdot \eta & 0 \\
& SYM & EI \Theta / \pi \\
\end{bmatrix} \label {114} \end {equation} $$
gdzie: $\Theta= \eta [ EA L^2 (\pi^2 -96) +6EI \pi^4 \eta^2 ]$
Macierz sztywności ($\ref{114}) jest przydatna do badania wpływu imperfekcji elementu belki – słupa. Badanie tych zależności nie jest przedmiotem artykułu.
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ą na podstawie energii odkształcenia oraz twierdzenia Castigliano.
Rozważmy element o osi środkowej opisanej funkcją sinusoidalną
\[ y(x)=e_0\sin\left(\frac{\pi x}{L}\right) \tag{II.A.1}\label{II.A.1}\]
gdzie $e_0$ oznacza amplitudę imperfekcji geometrycznej, natomiast $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.
\[ \frac{dy}{dx} = \frac{\pi e_0}{L} \cos\left(\frac{\pi x}{L}\right) \tag{II.A.2}\label{II.A.2} \]
\[ ds= \sqrt{ 1+ \left( \frac{\pi e_0}{L} \right)^2 \cos^2 \left( \frac{\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( \frac{\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=U_N+U_M \tag{II.A.5}\label{II.A.5} \]
\[ U_N= \int_0^L \frac{N^2}{2EA} \,ds \tag{II.A.6}\label{II.A.6} \]
\[ U_M= \int_0^L \frac{M^2(x)}{2EI} \,ds \tag{II.A.7}\label{II.A.7} \]
Po uwzględnieniu zależności (\ref{II.A.3}) i (\ref{II.A.4}) otrzymujemy
\[ U_N= \frac{N^2}{2EA} \int_0^L \sqrt{ 1+ \left( \frac{\pi e_0}{L} \right)^2 \cos^2 \left( \frac{\pi x}{L} \right) } \,dx \tag{II.A.8}\label{II.A.8} \]
\[ U_M= \frac{1}{2EI} \int_0^L \left[ M+Vx- Ne_0 \sin\left( \frac{\pi x}{L} \right) \right]^2 \sqrt{ 1+ \left( \frac{\pi e_0}{L} \right)^2 \cos^2 \left( \frac{\pi x}{L} \right) } \,dx \tag{II.A.9}\label{II.A.9}\]
Wzory (\ref{II.A.8}) i (\ref{II.A.9}) stanowią pełne energetyczne sformułowanie elementu imperfekcyjnego o osi sinusoidalnej. Na ich 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
\[ \{P\}= \begin{Bmatrix} N\\ V\\ M \end{Bmatrix} \tag{II.A.10}\label{II.A.10} \]
natomiast odpowiadający mu wektor przemieszczeń końcowych
\[ \{q\}= \begin{Bmatrix} u\\ v\\ \varphi \end{Bmatrix} \tag{II.A.11}\label{II.A.11} \]
Zgodnie z twierdzeniem Castigliano zachodzi zależność
\[ \{q\} = [F] \{P\} \tag{II.A.12}\label{II.A.12}\]
gdzie $[F]$ jest macierzą podatności elementu, której poszczególne skłądowe wyznacza się poprzez drugie pochodne energii odkształcenia względem sił uogólnionych
\[ f_{ij} = \frac{\partial^2 U} {\partial P_i\,\partial P_j} \tag{II.A.13}\label{II.A.13} \]
Zgodnie z zależnością (\ref{II.A.14}) poszczególne współczynniki macierzy podatności można wyznaczyć przez różniczkowanie energii odkształcenia względem odpowiednich sił uogólnionych. Ponieważ energia osiowa (\ref{II.A.8}) zależy wyłącznie od siły osiowej $N$, natomiast energia zginania (\ref{II.A.9}) zawiera wszystkie trzy siły uogólnione, współczynniki macierzy podatności można zapisać w postaci
\[ [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.14}\label{II.A.14} \]
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.8})–(\ref{II.A.9}) do definicji (\ref{II.A.13}) 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} = \frac{1}{EA} \int_0^L \sqrt{ 1+k^2\cos^2\left(\frac{\pi x}{L}\right) } \,dx + \frac{e_0^2}{EI} \int_0^L \sin^2\left(\frac{\pi x}{L}\right) \sqrt{ 1+k^2\cos^2\left(\frac{\pi x}{L}\right)
} \,dx =\\
\frac{2L}{\pi EA} \sqrt{1+k^2}\;E(m)+ \frac{2Le_0^2}{\pi EI} \sqrt{1+k^2} \left[ E(m) – \frac{E(m)-(1-m)K(m)}{m} \right]\tag{II.A.15}\label{II.A.15} $
$f_{12} = – \frac{e_0}{EI} \int_0^L x \sin\left( \frac{\pi x}{L} \right) \sqrt{ 1+ \left( \frac{\pi e_0}{L} \right)^2 \cos^2\left( \frac{\pi x}{L} \right) } \,dx\tag{II.A.16}\label{II.A.16}$
$ f_{13} = – \frac{e_0}{EI} \int_0^L \sin\left( \frac{\pi x}{L} \right) \sqrt{ 1+ \left( \frac{\pi e_0}{L} \right)^2 \cos^2\left( \frac{\pi x}{L} \right) } \,dx= – \frac{e_0L}{\pi EI} \left[ \sqrt{1+k^2} + \frac{\operatorname{arsinh}(k)}{k} \right] \tag{II.A.17}\label{II.A.17} $
$ f_{22} = \frac{1}{EI} \int_0^L x^2 \sqrt{ 1+ \left( \frac{\pi e_0}{L} \right)^2 \cos^2\left( \frac{\pi x}{L} \right) } \,dx =\\
\frac{L^3}{\pi^3EI} \int_0^\pi \theta^2 \sqrt{ 1+k^2\cos^2\theta }\,d\theta \tag{II.A.18}\label{II.A.18} $
$ f_{23} = \frac{1}{EI} \int_0^L x \sqrt{ 1+ \left( \frac{\pi e_0}{L} \right)^2 \cos^2\left( \frac{\pi x}{L} \right) }
\,dx = \frac{L}{2}\,f_{33} = \frac{L^2}{\pi EI} \sqrt{1+k^2}\;E(m) \tag{II.A.19}\label{II.A.19} $
$ f_{33} = \frac{1}{EI} \int_0^L \sqrt{ 1+k^2\cos^2\left(\frac{\pi x}{L}\right) } \,dx = \frac{2L}{\pi EI} \sqrt{1+k^2}\;E(m) \tag{II.A.20}\label{II.A.20} $
gdzie wprowadzono oznaczenia:
$ k=\frac{\pi e_0}{L}$,
$m=\frac{k^2}{1+k^2}$,
$K(m) = \int_0^{\pi/2} \frac{d\theta} {\sqrt{1-m\sin^2\theta}}$ oznacza zupełną całkę eliptyczną pierwszego rodzaju,
$E(m) = \int_0^{\pi/2} \sqrt{1-m\sin^2\theta}\,d\theta $ oznacza zupełną całkę eliptyczną drugiego rodzaju.
Współczynniki (\ref{II.A.15})–(\ref{II.A.20}) stanowią ścisłe przedstawienie macierzy podatności elementu imperfekcyjnego. Dla współczynników $f_{11}$, $f_{23}$ i $f_{33}$ uzyskano postać zamkniętą wyrażoną za pomocą funkcji eliptycznych, natomiast współczynniki $f_{12}$ i $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 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}$ wyrażone za pomocą funkcji elementarnych i całek eliptycznych. 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ółczynniki te można zapisać odpowiednio w postaci
$ f_{12} = – \frac{e_0}{EI} \int_0^L x \sin\left( \frac{\pi x}{L} \right) \sqrt{ 1+ \left( \frac{\pi e_0}{L} \right)^2 \cos^2\left( \frac{\pi x}{L} \right) } \,dx$ oraz
$ f_{22} = \frac{1}{EI} \int_0^L x^2 \sqrt{ 1+ \left( \frac{\pi e_0}{L} \right)^2 \cos^2\left( \frac{\pi x}{L} \right) } \,dx. $
Struktura całek oraz możliwość ich przybliżonego przedstawienia w postaci rozwinięć względem małego parametru imperfekcji
\[ \varepsilon=\frac{e_0}{L}. \tag{II.A.21}\label{II.A.21} \]
Podejście takie 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 bezwymiarową współrzędną
$ \xi=\frac{x}{L}, \qquad 0\le\xi\le1,$
oraz parametr
$ k=\pi\varepsilon. $
Współczynniki pozostające w postaci całkowej przyjmują wówczas postać
\[ f_{12} = -\frac{e_0L^2}{EI} \int_0^1 \xi \sin(\pi\xi) \sqrt{1+k^2\cos^2(\pi\xi) } \,d\xi \tag{II.A.24}\label{II.A.24}\]
\[ f_{22} = \frac{L^3}{EI} \int_0^1 \xi^2 \sqrt{ 1+k^2\cos^2(\pi\xi) } \,d\xi. \tag{II.A.25}\label{II.A.25} \]
Dla niewielkich imperfekcji geometrycznych zachodzi
$ k\ ll1, \tag{II.A.26}\label{II.A.26}
$$
co umożliwia rozwinięcie funkcji podpierwiastkowej w szereg Taylora
$$
\sqrt{1+z}
=
1+\frac{z}{2}
-\frac{z^2}{8}
+O(z^3).
\tag{II.A.27}\label{II.A.27}
$$
Po podstawieniu
$$
z=
k^2\cos^2(\pi\xi)
\tag{II.A.28}\label{II.A.28}
$$
otrzymujemy
$$
\sqrt{
1+k^2\cos^2(\pi\xi)
}
=
1
+
\frac{k^2}{2}\cos^2(\pi\xi)
–
\frac{k^4}{8}\cos^4(\pi\xi)
+
O(k^6).
\tag{II.A.29}\label{II.A.29}
$$
Podstawiając rozwinięcie (\ref{II.A.29}) do wzoru (\ref{II.A.24}) otrzymujemy
$$
f_{12}
=
-\frac{e_0L^2}{EI}
\int_0^1
\xi\sin(\pi\xi)\,d\xi
–
\frac{e_0L^2k^2}{2EI}
\int_0^1
\xi\sin(\pi\xi)\cos^2(\pi\xi)\,d\xi
+
O(k^4).
\tag{II.A.30}\label{II.A.30}
$$
Analogicznie dla współczynnika $f_{22}$
$$
f_{22}
=
\frac{L^3}{EI}
\int_0^1
\xi^2\,d\xi
+
\frac{L^3k^2}{2EI}
\int_0^1
\xi^2\cos^2(\pi\xi)\,d\xi
–
\frac{L^3k^4}{8EI}
\int_0^1
\xi^2\cos^4(\pi\xi)\,d\xi
+
O(k^6).
\tag{II.A.31}\label{II.A.31}
$$
Otrzymane całki zawierają już wyłącznie wielomiany oraz funkcje trygonometryczne i mogą zostać obliczone analitycznie metodami elementarnymi. Dzięki temu współczynniki $f_{12}$ oraz $f_{22}$ można przedstawić w postaci rozwinięć asymptotycznych względem parametru imperfekcji $\varepsilon=e_0/L$, zachowując jednocześnie pełną interpretację fizyczną kolejnych wyrazów szeregu.
Otrzymane całki zawierają już wyłącznie funkcje wielomianowe i trygonometryczne, dzięki czemu mogą zostać obliczone analitycznie metodami elementarnymi. W rezultacie współczynniki $f_{12}$ oraz $f_{22}$ można przedstawić w postaci szeregów potęgowych względem parametru imperfekcji $\varepsilon=e_0/L$, analogicznie do rozwinięć stosowanych wcześniej przy analizie rozwiązania Livesleya.
Macierz sztywności elementu
Macierz sztywności wyznacza się przez odwrócenie macierzy podatności
\[ [K]=[F]^{-1} \tag{II.A.23}\label{II.A.23} \]
Uzyskana w ten sposób macierz może być bezpośrednio wykorzystana w metodzie elementów skończonych. Klasyczna belka Bernoulliego-Eulera stanowi przypadek szczególny odpowiadający granicy $e_0 \rightarrow 0$, dla której oś środkowa staje się linią prostą
Po podstawieniu zależności (\ref{II.A.8})–(\ref{II.A.9}) 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. Są one bezpośrednią konsekwencją początkowej krzywizny osi środkowej.
Literatura
- Piechnik, St., Wytrzymałość materiałów dla wydziałów budowlanych. PWN, Warszawa-Kraków, 1980
- Livesley R.K. (1975), Matrix Methods of Structural Analysis, 2nd Edition, Pergamon Press, Oxford, Chapter 3.6.
- Livesley R. K. (1975), Matrix methods of structural analysis, Pergamon Press [http://books.google.com/books?id=tYsoAQAAMAAJ ]
________________________________
Related Baza wiedzy
