We found a solution of the completely linearized problem under the action of the external loading. Let us denote the corresponding displacement and stress field in a generic beam of a frame as $$ u_0(s),\; v_0(s),\;N_0,\;T_0, \; M_0(s);\qquad \text{(solution of the linear problem)} $$ We call it the fundamental solution. In the buckling analysis we study the existence of other solutions close to this fundamental solution. Since this solution is found by a linear analysis, it is linear in the applied loading. In particular, if we denote by $\bar N_0$ the axial stress in the bar for a reference distributioon of the external loadings, if all the loading are multiplied by a factor $\lambda$, the axial stress in the beam will be $N_0=\lambda\,\bar N_0$, where $\lambda$ is called load multiplier
We apply an Euler buckling analysis, where we suppose that the fundamental configuration is with negligible displacements ($u_0(s),v_0(s)\ll 1$) but possibly large internal stresses ($N_0, T_0, M_0\simeq 1$). This is possible for stiff frames ($ES, EI\gg 1$).
Hence, we introduce the effect of the geometrical nonlinearities and we perform a new linearization around this fundamental solution with negligible deflection but large stress.
We use the moderate rotations approximation (see the first notebook): $$ \begin{cases} u'(s)&=(1+\epsilon(s))\cos\theta(s) - 1\simeq \epsilon(s)-\dfrac{\theta(s)^2}{2}\\ v'(s)&=(1+\epsilon(s))\sin\theta(s) \simeq \theta(s)+o(\theta^2) \end{cases} $$
Within this approximation the deformations are given by $$ \begin{cases} \epsilon(s) &=& u'(s) + \dfrac{v'(s)^2}{2}\\ k(s) &=& v''(s) \end{cases} $$ Here, a non-constant transverse deflection induces an axial stretching.
With this approximation the elastic energy reads as $$ \mathcal{E}(u,v)=\int_0^L\dfrac{EI}{2} v''(s)^2 + \dfrac{ES}{2} \left(u'(s)+ \dfrac{v'(s)^2}{2}\right)^2\,ds $$
The energy above in non-linear (4th order) in $u$ and $v$. This means that the equilibrium equation (which is obtained after taking the derivative of the energy) will be of the 3rd order. We introduce now some approximations that allows us to obtain a quadratic version of the energy (linearized model around the prestressed state).
Expanding the last term and neglecting the contribution in $v'^4$ we get: $$ \mathcal{E}(u,v)=\int_0^L\dfrac{EI}{2} v''(s)^2 + \dfrac{ES}{2} u'(s)^2+ \dfrac{ES}{2}u'(s){v'(s)^2}\,ds $$ We further approximate the last cubic term by replacing ${ES}u'(s)$ with the normal stress of the fundamental solution $N_0$. Hence we finally have:
$$ \mathcal{E}(u,v)=\int_0^L\dfrac{EI}{2} v''(s)^2 + \dfrac{ES}{2} u'(s)^2+ N_0{v'(s)^2}\,ds $$In this new problem, there are no other external loading than the prestresses $N_0$. Hence, the total potential energy coincides with the sum of the elastic energies of each beam segment. The equilibrium configuration will then be the stationary points of this energy among all admissible variations.
Let us introduce the discretisation for the displacement in a beam segment described in the first notebook (linear in $u$ and cubic in $v$). Let us denote by $U^{(e)}$ the element dof vector, having as components the displacements and rotations at the two nodes of the element.
Using a linear interpolation for the axial displacement, the axial stress $N_0$ will be constant in each beam segment.
The discretised elastic energy of an element will then read as $$ \mathcal{E}(U^{(e)})=\dfrac{1}{2}K^{(e)}U^{(e)}\cdot U^{(e)}-\dfrac{\lambda}{2}G^{(e)}U^{(e)}\cdot U^{(e)} $$ where $K^{(e)}$ is the elastic stiffness matrix of the element (including bending and extension contribution) defined in the first notebook, $\lambda$ is the load multiplier of the fundamental solution, and $G^{(e)}$ is a new geometric stiffness.
Exercice. Show that $$ \boxed{G^{(e)}_{ij}=-\bar N_0^{(e)}\int_0^LS_{1i}'(s)S_{1j}'(s) ds},\quad i,j=1,\ldots,6 $$ where $S$ is the shape function matrix introduced in the first notebook and $\bar N_0^{(e)}$ the value of the axial stress in the fundamental solution in the element $e$ for a unit value of the load multiplier.
Remark. In buckling problems it is costumary to introduce the negative sign to get a positive matrix for compressive loadings.
We can obtain the energy of the whole frame by summing over the element and using the assembling procedure already used for the elastic stiffness matrix $K$. Hence, if we denote by $U$ the global nodal displacement vector for the whole frame, the total energy of the frame will read as
$$ \mathcal{E}(U)=\dfrac{1}{2}(K-\lambda G)U\cdot U $$We assume that the prestresses are the only loading acting on the system. Hence the total potential energy coincides with $\mathcal{E}$ above and the equilibrium of the frame is then given by the stationarity condition of this new energy:
$$ \dfrac{\partial\mathcal{E}}{\partial U}=(K-\lambda G)U = 0 $$To account for blocked degrees of freedom associated to Dirichlet boundary conditions, we can use one of the two techniques presented at the end of the first notebook. For example we can set to 0 the corresponding rows and columns in the $K$ and $G$ matrices, and set 1 on the diagonal.
The system above admits always the solution $U=0$. This concides with the fundamental solution.
The buckling problem consists in looking for the value of the load multiplier $\lambda$ for which the system $$ (K-\lambda G)U =0 $$ admits non null solutions. We will call these values of the load multiplier $\lambda_c^{(i)}$ the buckling loads and the corresponding non null solutions $U^{(i)}_c$ the buckling modes. This problem is a generalized eigenvalue problem that can be solved numerically with the techniques introduced in the dedicated notebook.
In summary, to find the buckling loads and modes one can use the following practical recipe.
Consider a given frame with a given loading distribution.
Solve the linear problem $K U =F$ of the notebook 2 under the action of the given loadings. This will give the fundamental solution.
Extract the normal stress in each bar $\bar N^{(e)}$ for the fundamental solution, that we can collect in a vector having as length the number of elements.
Assemble the geometric stiffness $G$ and modify it to account for the Dirichlet Boundary conditions.
Solve the eigenvalue problem $$ (K-\lambda \,G)U =0 $$ which gives the buckling load multipliers $\lambda_c^{(i)}$ and the buckling modes $U^{(i)}_c$.
We can conclude from this analysis that the fundamental solution is locally unique (and stable, indeed) for load multipliers smaller than the first critical load $\lambda_c^{(0)}$. At $\lambda=\lambda_c^{(0)}$ a new solution departing from the fundamental solution with the shape associated with the first buckling mode $U^{(0)}_c$ is possible. We can also show, that for $\lambda>\lambda_c^{(0)}$ the fundamental solution will be unstable.
In [ ]: