2D Frame System Solver

Use html5 supporting browser

New System

The larger of total system width and height[mm]

x[mm] y[mm] x constraint y constraint rotation constraint

x-Force[N] y-Force[N] B.Moment[Nmm] Add Node

Young's modulus[MPa] Area[mm2] Iz[mm4]

First Node Second Node Dist. load[kN/m] Add Bar

Member List Solve

Direct Stiffness Method for 2D Frames

The elements that make up a frame structure are capable of carrying shear forces and bending moments in addition to the axial forces. Also, in addition to the translational degrees of freedom at the two nodes of an element, frame members have rotational degrees of freedom. The numbering convention and the positive directions of these degrees of freedom are shown below.

Figure 1: Translational and rotational degrees of freedom of a frame member

The first step in the derivation of the element stiffness matrix is to describe the flexural displacement \(v(x)\) and axial displacement \(u(x)\) of the frame members in polynomial form.

\[v(x) = c_1+c_2x+c_3x^2+c_4x^3\]

\[u(x) = c_5+c_6x\]

In the Euler-Bernoulli beam theory the equation for the bending moment is \(M(x)=EI\displaystyle\frac{d^2v}{dx^2}\). Also, since in the finite element method loads are applied at the nodes of an element, \(M(x)\) must vary linearly between the nodes (the contributions of the distributed loads to the bending moments are added after the system equations are solved). Therefore \(v(x)\) is described as a third order polynomial so that its second derivative varies linearly between the element nodes. The description of the axial displacement as a first order polynomial follows from the fact that the axial force carried by a member is assumed to be constant along the member length. This leads to constant axial strain which has the equation \(\varepsilon_x=\displaystyle\frac{du}{dx}\). After the application of the boundary conditions \(u(x)\vert_{x=0}=u_1\), \(v(x)\vert_{x=0}=v_1\), \(v'(x)\vert_{x=0}=\theta_1\), \(u(x)\vert_{x=L}=u_2\), \(v(x)\vert_{x=L}=v_2\), \(v'(x)\vert_{x=L}=\theta_2\), we obtain the following expression for \(v(x)\) and \(u(x)\):

\[v(x) = \Big(1-\frac{3x^2}{L^2}+\frac{2x^3}{L^3}\Big)v_1+\Big(x-\frac{2x^2}{L}+\frac{x^3}{L^2}\Big)\theta_1+\Big(\frac{3x^2}{L^2}-\frac{2x^3}{L^3}\Big)v_2+\Big(\frac{x^3}{L^2}-\frac{x^2}{L}\Big)\theta_2\]

\[u(x) = u_1+\frac{(u_2-u_1)}{L}x\]

The strain energy in a member can be computed as \(U=\displaystyle\frac{1}{2}\displaystyle\int_V\sigma_x\varepsilon_xdV\) where \(V\) is the total volume of the member. The strain energy can be computed as the superposition of the strain energies related to flexure \((U_f)\) and axial loading \((U_a)\). In case of axial loading stress, strain and resulting strain energy are computed as follows:

\[\varepsilon_x = \frac{du}{dx}=\frac{(u_2-u_1)}{L},\quad \sigma_x=E\varepsilon_x=E\frac{du}{dx}\]

\[U_a = \frac{1}{2}\int_VE\Big(\frac{du}{dx}\Big)^2dV=\frac{1}{2}E\int_0^L\Big(\frac{du}{dx}\Big)^2\int_AdAdx=\frac{1}{2}EA\int_0^L\frac{(u_2-u_1)^2}{L^2}dx\]

\[U_a = \frac{EA}{2L}(u_2-u_1)^2\]

The stress, strain and strain energy associated with the flexural loading can be computed as follows:

\[\varepsilon_x = -y\frac{d^2v}{dx^2}=-y(2c_3+6c_4x),\quad \sigma_x=E\varepsilon_x=-yE\frac{d^2v}{dx^2}\]

\[U_f = \frac{1}{2}E\int_Vy^2\Big(\frac{d^2v}{dx^2}\Big)^2dV=\frac{1}{2}E\int_0^L\Big(\frac{d^2v}{dx^2}\Big)^2\int_Ay^2dAdx=\frac{1}{2}EI\int_0^L\Big(\frac{d^2v}{dx^2}\Big)^2dx\]

At this point it is convenient to describe \(v(x)\) as a combination of 4 shape functions

\[N_1=1-\frac{3x^2}{L^2}+\frac{2x^3}{L^3},\quad N_2=x-\frac{2x^2}{L}+\frac{x^3}{L^2}\]

\[N_3=\frac{3x^2}{L^2}-\frac{2x^3}{L^3},\quad N_4=\frac{x^3}{L^2}-\frac{x^2}{L}\]

\[v(x) = N_1v_1+N_2\theta_1+N_3v_2+N_4\theta_2\]

\[U_f =\frac{1}{2}EI\int_0^L\Big(\frac{d^2N_1}{dx^2}v_1+\frac{d^2N_2}{dx^2}\theta_1+\frac{d^2N_3}{dx^2}v_2+\frac{d^2N_4}{dx^2}\theta_2\Big)^2dx\]

The total strain energy of the frame member can be computed using the superposition of \(U_a\) and \(U_f\) as follows:

\[U =\frac{EA}{2L}(u_2-u_1)^2+\frac{1}{2}EI\int_0^L\Big(\frac{d^2N_1}{dx^2}v_1+\frac{d^2N_2}{dx^2}\theta_1+\frac{d^2N_3}{dx^2}v_2+\frac{d^2N_4}{dx^2}\theta_2\Big)^2dx\]

Once the total strain energy of a member is known, the forces and moments in the local coordinates acting at its nodes can be computed using Castigliano’s first theorem such that:

\[N_1 =\frac{\partial U}{\partial u_1}=\frac{EA}{L}(u_1-u_2)\]

\[V_1 =\frac{\partial U}{\partial v_1}=EI\int_0^L(N_1^{''}N_1^{''}v_1+N_2^{''}N_1^{''}\theta_1+N_3^{''}N_1^{''}v_2+N_4^{''}N_1^{''}\theta_2)dx\]

\[M_1 =\frac{\partial U}{\partial \theta_1}=EI\int_0^L(N_1^{''}N_2^{''}v_1+N_2^{''}N_2^{''}\theta_1+N_3^{''}N_2^{''}v_2+N_4^{''}N_2^{''}\theta_2)dx\]

\[N_2 =\frac{\partial U}{\partial u_2}=\frac{EA}{L}(u_2-u_1)\]

\[V_2 =\frac{\partial U}{\partial v_2}=EI\int_0^L(N_1^{''}N_3^{''}v_1+N_2^{''}N_3^{''}\theta_1+N_3^{''}N_3^{''}v_2+N_4^{''}N_3^{''}\theta_2)dx\]

\[M_2 =\frac{\partial U}{\partial \theta_2}=EI\int_0^L(N_1^{''}N_4^{''}v_1+N_2^{''}N_4^{''}\theta_1+N_3^{''}N_4^{''}v_2+N_4^{''}N_4^{''}\theta_2)dx\]

The above equations can be written in matrix form:

\[\begin{split}\begin{bmatrix} N_1 \\ V_1 \\ M_1 \\N_2\\V_2\\M_2 \end{bmatrix} = \begin{bmatrix} \displaystyle\frac{EA}{L} &0&0&-\displaystyle\frac{EA}{L}&0&0\\ 0 &EI\int_0^L N_1^{''}N_1^{''}dx&EI\int_0^L N_2^{''}N_1^{''}dx & 0&EI\int_0^L N_3^{''}N_1^{''}dx&EI\int_0^L N_4^{''}N_1^{''}dx \\ 0 &EI\int_0^L N_1^{''}N_2^{''}dx&EI\int_0^L N_2^{''}N_2^{''}dx& 0 & EI\int_0^LN_3^{''}N_2^{''}dx & EI\int_0^L N_4^{''}N_2^{''}dx\\ -\displaystyle\frac{EA}{L}&0&0&\displaystyle\frac{EA}{L}&0&0\\0&EI\int_0^L N_1^{''}N_3^{''}dx &EI\int_0^L N_2^{''}N_3^{''}dx &0&EI\int_0^L N_3^{''}N_3^{''}dx&EI\int_0^L N_4^{''}N_3^{''}dx\\ 0&EI\int_0^L N_1^{''}N_4^{''}dx & EI\int_0^L N_2^{''}N_4^{''}dx & 0 & EI\int_0^L N_3^{''}N_4^{''}dx & EI\int_0^L N_4^{''}N_4^{''}dx\end{bmatrix}\begin{bmatrix} u_1\\v_1\\\theta_1\\u_2\\v_2\\\theta_2\end{bmatrix}\end{split}\]

After evaluating the integrals in the above matrix equation, the frame member stiffness matrix in local coordinates is found as

\[\begin{split}\mathbf{k^{'}} = \begin{bmatrix} \displaystyle\frac{EA}{L} &0&0& -\displaystyle\frac{EA}{L}&0&0 \\ 0&\displaystyle\frac{12EI}{L^3}&\displaystyle\frac{6EI}{L^2}&0&-\displaystyle\frac{12EI}{L^3}&\displaystyle\frac{6EI}{L^2}\\ 0&\displaystyle\frac{6EI}{L^2}&\displaystyle\frac{4EI}{L}&0&-\displaystyle\frac{6EI}{L^2}&\displaystyle\frac{2EI}{L}\\ -\displaystyle\frac{EA}{L} &0&0& \displaystyle\frac{EA}{L}&0&0\\ 0&-\displaystyle\frac{12EI}{L^3}&-\displaystyle\frac{6EI}{L^2}&0&\displaystyle\frac{12EI}{L^3}&-\displaystyle\frac{6EI}{L^2}\\ 0&\displaystyle\frac{6EI}{L^2}&\displaystyle\frac{2EI}{L}&0&-\displaystyle\frac{6EI}{L^2}&\displaystyle\frac{4EI}{L} \end{bmatrix}\end{split}\]

The transformation of the stiffness matrices into the global coordinate system and the assemblage of the global stiffness matrix can be done similar to 2 dimensional trusses.

In case of elements carrying distributed loading, the reaction forces that the distributed load would cause on a single beam element, are added to the load vectors of the element nodes. These reaction forces and moments are shown in Figure 2.

Figure 2: Reaction forces and moments of a frame member

In Figure 2 the nodes at the two ends of the frame element are denoted with \(N_1\) and \(N_2\). From elementary strength of materials we know that a distributed load \(q\) on a beam of length \(L\) clamped at both ends would cause end moments in magnitude \(\displaystyle\frac{qL^2}{12}\) and support shear forces in magnitude \(\displaystyle\frac{qL}{2}\). In Figure 2, the reaction forces and moments acting in the opposite directions of the degrees of freedom given in Figure 1, have negative sign. On the other hand the forces and moments acting on the nodes of the element because of the distributed load, act in the opposite directions of the reaction forces and moments. A nodal force or moment caused by \(q\) is only added to the system load vector if the node is not constrained in the direction of that force or moment.

After solving the system equations, the shear forces and bending moments are transformed to local coordinates. Finally according to the sign convention of Figure 1, \(qL/2\) and \(qL^2/12\) are added to the local forces and moments.

References

[1] Hibbeler R.C., Structural Analysis, 8th edition, ISBN:978-0132570534

[2] Hutton D.V., Fundamentals of Finite Element Analysis, ISBN:0072395362