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 78 Czytelników
CZĘŚC II: NIELINIOWOŚĆ I IMPERFEKCJE geometryczne oraz 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 podatnoś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ą z energii odkształcenia przy wykorzystaniu 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, 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.
\[ \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 = \int_0^L \frac{N^2}{2EA} \,ds + \int_0^L \frac{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 \{ \frac{N^2}{2EA} + \frac{1}{2EI} \left[ M+Vx- Ne_0 \sin\left( \frac{\pi x}{L} \right) \right]^2 \right\} \sqrt{ 1+ \left( \frac{\pi e_0}{L} \right)^2 \cos^2 \left( \frac{\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} = \frac{\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} = \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.10}\label{II.A.10} $
$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.11}\label{II.A.11}$
$ 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.12}\label{II.A.12} $
\[ 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.13}\label{II.A.13} \]
$ 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.14}\label{II.A.14} $
$ 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.15}\label{II.A.15} $
gdzie wprowadzono oznaczenia:
- pomocnicze parametery imperferkcji
\[ k = \frac{\pi e_0}{L} = \pi\varepsilon. \tag{II.A.16}\label{II.A.16} \]
\[ m=\frac{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} \frac{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=\frac{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=\frac{x}{L}, \qquad 0\le\xi \le 1,$
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.21}\label{II.A.21} \]
\[ f_{22} = \frac{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 + \frac{k^2}{2}\cos^2(\pi\xi) – \frac{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} = -\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 +\ldots \tag{II.A.25}\label{II.A.25} \]
\[ 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 +\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 = \frac{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 = \frac{1}{3}. \tag{II.A.28}\label{II.A.28} \]
\[ f_{12} = -\frac{e_0L^2}{\pi EI} + O(k^2) \tag{II.A.29}\label{II.A.29}\]
\[ f_{22} = \frac{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 = \frac{L}{EA} + \frac{L^3}{2EI}\varepsilon^2 + \frac{\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 + \frac{\pi^2\varepsilon^2}{2}\cos^2(\pi\xi) – \frac{\pi^4\varepsilon^4}{8}\cos^4(\pi\xi) + O(\varepsilon^6)$
$ \int_0^1\sin^2(\pi\xi)\,d\xi = \frac12$
$ \int_0^1\sin^2(\pi\xi)\cos^2(\pi\xi)\,d\xi = \frac18 $,
Składowa $f_{12}$
\[ f_{12} = -\frac{e_0L^2}{EI} \left[ \int_0^1 \xi\sin(\pi\xi)\,d\xi + \frac{k^2}{2} \int_0^1 \xi\sin(\pi\xi)\cos^2(\pi\xi)\,d\xi \right] + O(k^4) = -\frac{L^3}{\pi EI}\varepsilon -\frac{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łasności całki:
$\int_0^1 \xi\sin(\pi\xi) \cos^2(\pi\xi)\,d\xi = \frac{5}{9\pi}$
Składowa $f_{13}$
Współczynnik sprzęgający oddziaływania osiowe i zginające przyjmuje postać
\[ f_{13} = f_{13}^{(1)}\varepsilon + f_{13}^{(3)}\varepsilon^3 +\ldots = -\frac{\varepsilon L^2}{EI} \int_0^1 \sin(\pi\xi)\,d\xi – \frac{\pi^2\varepsilon^3L^2}{2EI} \int_0^1 \sin(\pi\xi)\cos^2(\pi\xi)\,d\xi + O(\varepsilon^5) =-\frac{2L^2}{\pi I}\varepsilon -\frac{\pi L^2}{3EI}\varepsilon^3 + O(\varepsilon^5). \tag{II.A.33}\label{II.A.33} \]
Składowa $f_{22}$
\[f_{22} = \frac{L^3}{EI} \left[ \int_0^1 \xi^2\,d\xi + \frac{k^2}{2} \int_0^1 \xi^2\cos^2(\pi\xi)\,d\xi – \frac{k^4}{8} \int_0^1 \xi^2\cos^4(\pi\xi)\,d\xi \right] + O(k^6) = \frac{L^3}{3EI} + \frac{L^3}{EI} \left( \frac{\pi^2}{12} + \frac18 \right)\varepsilon^2 – \frac{L^3}{EI} \left( \frac{\pi^4}{64} + \frac{15\pi^2}{512} \right)\varepsilon^4 + O(\varepsilon^6).\tag{II.A.34}\label{II.A.34} \]
Przy wyprowadzeniu składowej $f_{22}$ wykorzystano zależności:
$\int_0^1 \xi^2\cos^2(\pi\xi)\,d\xi = \frac16+\frac{1}{4\pi^2}$ , $ \int_0^1 \xi^2\cos^4(\pi\xi)\,d\xi = \frac18+\frac{15}{64\pi^2}$.
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 = \frac{L^2}{2EI} +\frac{\pi^2L^2}{8EI}\varepsilon^2 -\frac{3\pi^4L^2}{128EI}\varepsilon^4 +O(\varepsilon^6). \tag{II.A.35}\label{II.A.35} \]
Składowa $f_{33}$
\[ f_{33} = f_{33}^{(0)} + f_{33}^{(2)}\varepsilon^2 + f_{33}^{(4)}\varepsilon^4 +\ldots \= \frac{L}{2EI} -\frac{\pi^2L}{8EI}\varepsilon^2 +\frac{9\pi^4L}{128EI}\varepsilon^4 + O(\varepsilon^6). tag{II.A.36}\label{II.A.36} \]
Uwaga:
Dla współczynników \(f_{23}\) i \(f_{33}\) skorzytsano z wcześniej uzyskanych postaci 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) = \frac{\pi}{2} \left( 1+\frac{m}{4} +\frac{9m^2}{64} +O(m^3) \right)$
$ E(m) = \frac{\pi}{2} \left( 1-\frac{m}{4} -\frac{3m^2}{64} +O(m^3) \right)$
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.
Wyzncznik $D_F$ macierzy podatności $\mathbf{F}$
Wyznacznik macierzy podatności
\[ D_F=\det[F]. \tag{II.A.37}\label{II.A.37} \]
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.38}\label{II.A.38} \]
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.39}\label{II.A.39} \]
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.48}\label{II.A.48} \]
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.49}\label{II.A.49} \]
Po podstawieniu wartości zerowego rzędu
\[ f_{11}^{(0)}=\frac{L}{EA}, \qquad f_{22}^{(0)}=\frac{L^3}{3EI}, \qquad f_{23}^{(0)}=\frac{L^2}{2EI}, \qquad f_{33}^{(0)}=\frac{L}{EI} \]
otrzymujemy
\[ D_F^{(0)} = \frac{L}{EA} \left( \frac{L^4}{3(EI)^2} – \frac{L^4}{4(EI)^2} \right) = \frac{L^5}{12\,EA\,(EI)^2}. \tag{II.A.49}\label{II.A.49} \]
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.51}\label{II.A.51} \]
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.38}) 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.52}\label{II.A.52} \]
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.53}\label{II.A.53} \]
W analogiczny sposób współczynnik czwartego rzędu przyjmuje postać
\[ D_F^{(4)} = D_F^{(4)}_{(a)} + D_F^{(4)}_{(b)} + D_F^{(4)}_{(c)}, \tag{II.A.54}\label{II.A.54} \]
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.55}\label{II.A.55} \]
a kolejne ……
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.56}\label{II.A.56} \]
Zachodzi zależność
\[ D_K = \det[K] = \frac{1}{\det[F]} = \frac{1}{D_F}. \tag{II.A.57}\label{II.A.57} \]
Podstawiając rozwinięcie asymptotyczne wyznacznika podatności (\ref{II.A.48}) i rozwijając odwrotność w szereg względem parametru \(\varepsilon\), otrzymujemy
\[ D_K = \frac{1}{D_F^{(0)}} \left[ 1 – \frac{D_F^{(2)}}{D_F^{(0)}}\varepsilon^2 + \left( \frac{\left(D_F^{(2)}\right)^2}{\left(D_F^{(0)}\right)^2} – \frac{D_F^{(4)}}{D_F^{(0)}} \right)\varepsilon^4 + O(\varepsilon^6) \right].\tag{II.A.58}\label{II.A.58} \]
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.56}). 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
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=\frac{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 redukcyjneg (wyvoczenowego) 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
- 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
