MAXSURF | MOSES | SACS | OpenWindPower - Geometric Stiffness

Introduction

SACS performs P-Delta analyses using a geometric stiffness matrix which is derived from the potential energy function of each element due to axial, bending and shear effects. In this post, I will show how the geometric stiffness matrix is derived and applied to the elastic stiffness matrix for the structural analysis. This derivation is based upon the formulation presented in Concepts and Applications of Finite Element Analysis (Cook 1989).

Geometric Stiffness Matrix Derivation

In the case of the elastic stiffness matrix, the displacement of an element is dependent upon the force and stiffness of the structure. However as an element is displaced, the stiffness of the structure changes. The geometric stiffness matrix represents the change in stiffness of the structure as elements are displaced.

First let's consider an cable element with lateral loads at each end.

Cable Element

The formula for the relative deflection between the ends would be:

$F_i=\frac{T}{L}(v_i-v_j)$

And static equilibrium requires that:

$F_j=-F_i$

The resulting geometric stiffness matrix is:

$\begin{bmatrix} F_i\\ F_j\\ \end{bmatrix} =\frac{T}{L} \begin{bmatrix} 1 & -1\\ -1 & 1\\ \end{bmatrix} \begin{bmatrix} v_i\\ v_j\\ \end{bmatrix}$

Now if we consider a deformed beam element.

Beam Element

The deformed shape of the frame element, h(x), can be expressed as:

$h(x) = a_0+a_1x+a_2x^2+a3x^3$

The slope of the element is the derivative of the shape with respect to x:

$h'(x) = a_1+2a_2x+3a_3x^2$

Assuming small rotations the end displacements and rotations at each end of the beam can be solved as:

$h(0) = v_i = a_0$

$h'(0) = \theta_i = a_1$

$h(L) = v_j = a_0+a_1L+a_2L^2+a_3L^3$

$h'(L) = \theta_j = a_1+2a_2L+3a_3L^2$

solving for the variables $a_1$ through $a_3$ :

$a_0=v_i$

$a_1=\theta_i$

$a_2=3\frac{v_j-v_i}{L^2}-\frac{2\theta_i+\theta_j}{L}$

$a_3 = -2\frac{v_j-v_i}{L^3}+\frac{\theta_i+\theta_j}{L^2}$

The cubic deformation function can also be written as a weighted sum of cubic polynomials:

$h(x) = v_i b_2(x) + \theta_i b_3(x) + v_j b_5(x) + \theta_j b_6(x)$

$b_2(x) = 1-3(x/L)^2 + 2(x/L)^3$

$b_3(x) = x (1-x/L)^2$

$b_5(x) = 3 (x/L)^2 - 2(x/L)^3$

$b_6(x) = (x/L)^2(x/L-1)$

where:

$b_2(0) = 1$

$b'_3(0) = 1$

$b_5(L) = 1$

$b'_6(L) = 1$

The axial deformations can be expressed using the following equations:

$b_1(x) = (1-x/L)$

$b_4(x) = (x/L)$

The internal elastic strain energy of the beam element due to bending is:

$U_b = \int_0^L\frac{M^2(x)}{EI}dx$

Assuming small deformations the internal strain energy can be rewritten as:

$U_b =\frac{1}{2}\int_0^LM(x)h''(x)dx = \frac{1}{2} EI\int_0^L h''(x)^2dx$

A beam element will undergo axial deformation $u'(x)$ due axial loading and transverse rotation $h'(x)$ . The behavior can be expressed using the Pythagorean theorem:

$(dx)^2 = (dx-du)^2 + dh^2$

$(dx)^2 = dx^2 -2dudx + du^2 + dh^2$

$2dudx = du^2 + dh^2$

$\frac{du}{dx} = \frac{1}{2}\left(\frac{du}{dx}\right)^2 + \frac{1}{2}\left(\frac{dh}{dx}\right)^2$

Because axial strain is traditionally much smaller than the transverse rotations of a beam element:

$u'(x) = 1/2(h'(x))^2$

Because $T(x)$ does not increase linearly with $h(x)$ the internal strain energy force due to load $T$ can be expressed as:

$U_G = \int_0^L T(x) \frac{du(x)}{dx}dx = \frac{1}{2} \int_0^L T(x) \left(\frac{dh(x)}{dx}\right)^2dx$

If the axial load is constant over the length of the beam then the geometric strain energy is then:

$U_G = \frac{1}{2} T \int_0^L (h'(x))^2 dx$

If we then substitute the equation for $h'(x)$ into the equation we get:

$U_G = \frac{T}{30L} (-L\theta_i\theta_j - 3v_jL\theta_i - 3v_jL\theta_j + 3v_iL\theta_i + 3v_iL\theta_j + 18 v_j^2 - 36 v_jv_i + 18 v_i^2 + 2L^2\theta_i^2 + 2L^2\theta_j^2)$

Using Castigliano's theorem, the partial derivative of the geometric strain energy function with respect to a displacement coordinate yields the force of that displacement coordinate. The end forces for each coordinate can be expressed as:

$F_i = \frac{\partial U_G}{\partial v_i} = \frac{T}{30L} (36v_i + 3L\theta_i - 36 v_j + 3L \theta_j)$

$M_i = \frac{\partial U_G}{\partial \theta_i} = \frac{T}{30L} (3Lv_i + 4L^2\theta_i - 3L v_j - L^2 \theta_j)$

$F_j = \frac{\partial U_G}{\partial v_j} = \frac{T}{30L} (-36v_i - 3L\theta_i + 36 v_j - 3L \theta_j)$

$M_j = \frac{\partial U_G}{\partial \theta_j} = \frac{T}{30L} (3Lv_i - L^2\theta_i - 3L v_j + 4L^2 \theta_j)$

This expression can be rewritten in matrix form which gives us the geometric stiffness matrix:

Geometric Stiffness Matrix

The elastic stiffness solution can be expressed as:

Elastic Stiffness Matrix

So the resulting total stiffness solution is:

$F_T=F_E+F_G=[k_e+k_g]v = k_tv$

SACS Example

We will now consider twoP-Delta benchmark problems from the AISC 360 specification.

Case 1

Case 1 is a simply supported beam with a vertical load at one end and a transverse lateral along the length.

Case 1

The load P varies from 0 to 450 kips. Both Bernouli and Timoshenko beams are considered in this problem. The following results are compared with the published AISC values:

Axial Force	0	150	300	450
M_mid (kip-in) SD	235	270	316	380
M_mid (kip-in)	235	269	313	375
Delta_mid (in) SD	0.202	0.230	0.269	0.322
Delta_mid (in)	0.197	0.224	0.261	0.311

Note that the SD indicates shear deformation (Timoshenko beams).

The SACS results from the attached analyses are:

Axial Force	0	150	300	450
M_mid (kip-in) SD	235.2	269.64	315.43	379.26
M_mid (kip-in)	235.2	268.77	313.1	374.32
Delta_mid (in) SD	0.201	0.230	0.267	0.320
Delta_mid (in)	0.197	0.224	0.260	0.309

and the error (%) for each value is:

Axial Force	0	150	300	450
M_mid (kip-in) SD	-0.09	0.13	0.18	0.19
M_mid (kip-in)	-0.09	0.09	-0.03	0.18
Delta_mid (in) SD	0.42	0.17	0.58	0.58
Delta_mid (in)	0.17	0.08	0.51	0.59

Case 2

Case 2 is a cantilever column vertical load and lateral along at the tip:

Case 2

The load P varies from 0 to 200 kips. Again both Bernouli and Timoshenko beams are considered in this problem. The following results are compared with the published AISC values:

Axial Force	0	100	150	200
M_base (kip-in) SD	336	470	601	856
M_base (kip-in)	336	469	598	848
Delta_tip (in) SD	0.907	1.34	1.77	2.6
Delta_tip (in)	0.901	1.33	1.75	2.56

Note that the SD indicates shear deformation (Timoshenko beams).

The SACS results from the attached analyses are:

Axial Force	0	100	150	200
M_base (kip-in) SD	336	469.93	600.78	854.81
M_base (kip-in)	336	468.66	597.60	845.95
Delta_tip (in) SD	0.905	1.339	1.765	2.594
Delta_tip (in)	0.899	1.323	1.744	2.550

and the error (%) for each value is:

Axial Force	0	100	150	200
M_base (kip-in) SD	0.00	0.01	0.04	0.14
M_base (kip-in)	0.00	0.07	0.07	0.24
Delta_tip (in) SD	0.18	0.05	0.27	0.23
Delta_tip (in)	0.22	0.26	0.34	0.40

Conclusion

The derivation of the geometric stiffness matrix used in SACS was shown. The comparison with the benchmark values in AISC show that the implementation of the geometric stiffness matrix is with the 3% and 5% error for the moment and displacement respectively indicating that the SACS provides a valid P-Delta solution.