3.4 Stiffness matrices

Back

To calculate the stiffness matrix of the bar in Fig. 3.18, the differential equation

that relates the axial load p to the axial displacement u(x),

− EAu

__(x) = p(x) (3.72)

and the general homogeneous solution of this equation

uh(x) = a0 + a1x (3.73)

3.4 Stiffness matrices 293

3.18. Unit displacedownwards

must both be known. By a proper choice of the coefficients a0 and a1, the

unit displacements

ϕ1(x) =

1 − x

l

ϕ1(0) = 1, ϕ1(l) = 0,

ϕ2(x) = x

l

ϕ2(0) = 0, ϕ2(l) = 1

(3.74)

are derived. Finally Green’s first identity is formulated,

G(u, ˆu) =

_ l

0

−EAu

__ˆudx + [Nˆu] l

0

_ _ _

δWe

_ l

0

EAu

_ ˆu

_

dx

_ _ _

δWi

= 0 N = EAu

_

(3.75)

to provide the definition of the strain energy product δWi, because the element

k ij of the stiffness matrix K is the strain energy product between the

unit displacements, ϕi and ϕj

k ij =

_ l

0

EAϕ

_

iϕ

_

j dx =

_ l

0

NiNj

EA

dx , (3.76)

so that

K = EA

l

_

1 −1

−1 1

_

(3.77)

or Ku = f +p, where the fi are the end forces, and the pi are the fixed-end

forces (×−1) resulting from the distributed load

p 1 =

_ l

0

p ϕ1 dx p 2 =

_ l

0

p ϕ2 dx . (3.78)

Fig.

displacements are plotted

ments of the bar. Horizontal

294 3 Frames

Fig. 3.19. Variable stiffness

EA = EA(x)

Tapered beam

In a tapered beam with a cross-sectional area such as

A(x) = A0 + A1x , (3.79)

the differential equation for the longitudinal displacement is

− EA(x)u

__(x) − EA

_(x)u

_(x) = p(x) . (3.80)

With Maple or Mathematica, the homogeneous solution

u(x) = c2 + c1

lnA(x)

A1

(3.81)

can be found, and thereby the unit displacements

ϕ1(x) =

lnA(x) − lnA(l)

lnA(0) − lnA(l), ϕ2(x) =

lnA(0) − lnA(x)

lnA(0) − lnA(l) . (3.82)

Upon substituting these into the strain energy product

k ij =

_ l

0

EA(x) ϕ

_

iϕ

_

j dx , (3.83)

the stiffness matrix is obtained:

K = k

_

1 −1

−1 1

_

k = A1 E

lnA(l) − lnA0

(lnA1 − lnA0)2 . (3.84)

According to the formula for the fixed-end forces

p1 =

_ l

0

p ϕ1 dx , p2 =

_ l

0

p ϕ2 dx (3.85)

a single force P will generate the following actions at the fixed ends

p1 = P · ϕ1(xP ), p2 = P · ϕ2(xP ), xP = location of P . (3.86)

3.4 Stiffness matrices 295

Fig. 3.20. Unit displacements of a

beam

The sum of the support reactions pi must equal P = p1 + p2, i.e., at each

point x the sum of ϕ1(x) and ϕ2(x) must be 1

ϕ1(x) + ϕ2(x) = 1 (100%) . (3.87)

If the cross-sectional area varies as

A(x) = A0 + A1x = 1+1 · x length l = 1 (3.88)

and if a single force P acts at the center of the shaft (see Fig. 3.19), then

because

ϕ1(0.5) = 0.415 ϕ2(0.5) = 0.585 (3.89)

about 41% of P will pull at the left end and about 59% will press on the other

end of the shaft.

Euler–Bernoulli beam

This technique for the derivation of stiffness matrices can be applied to all

possible deformations in a beam, such as the twist ϕ(x) of the axis, shear

deformations ws(x), or simply the deflection w(x) in which case the differential

equation is

EIwIV (x) = p(x) . (3.90)

Starting with the homogeneous solution

wh(x) = a0 + a1x + a2x2 + a3x3 , (3.91)

296 3 Frames

the four unit displacements (see Fig. 3.20) are found:

ϕ1(x) = 1 − 3x2

l2 +

2x3

l3

ϕ2(x) = −x +

2x2

l

− x3

l2

ϕ3(x) =

3x2

l2

− 2x3

l3

ϕ4(x) = x2

l

− x3

l2 .

(3.92)

Green’s first identity is

G(w, ˆ w) =

_ l

0

EIwIV ˆ wdx + [V ˆ w −M ˆ w

_]l

0

_ _ _

δWe

_ l

0

EIw

__ ˆ w

__

dx

_ _ _

δWi

= 0 (3.93)

(3.94)

and the element k ij of the stiffness matrix K is the strain energy product

between the unit displacements ϕi and ϕj

k ij =

_ l

0

EI ϕ

__

i ϕ

__

j dx , (3.95)

so that

K = EI

l3

⎢⎢⎣

12 −6l −12 −6l

−6l 4l2 6l 2l2

−12 6l 12 6l

−6l 2l2 6l 4l2

⎥⎥⎦

. (3.96)

The (negative) end-fixing forces p i (= equivalent nodal forces) of a distributed

load p are the scalar product of p and the unit displacements:

p i =

_ l

0

pϕi dx . (3.97)

Timoshenko beam

In the following it is assumed that the bending stiffness EI, the effective shear

cross section As and the shear modulus G are constant. The deformations of

the beam are described by the deflection w and the rotation θ (see Fig. 3.21).

The constitutive equations are

strains θ

_ − κ = 0 w

_ + θγ = 0 (3.98)

material law GAsγ − V = 0 EI κ −M = 0 (3.99)

equilibrium M

_ − V = 0 −V

_ = p (3.100)

3.4 Stiffness matrices 297

Fig. 3.21. Timoshenko beam

or

− EI θ

__ + GAs (w

_ + θ) = 0 (3.101)

−GAs (w

__ + θ

_) = p . (3.102)

The latter system can be read as the application of an operator −L (minus

because of second order) to the vector-valued function u = [w, θ]T. Green’s

first identity for this operator −L is

G(uu) =

_ l

0

Lu•ˆu dx +

_

V ˆ w +M ˆθ

_l

0

− a(uu) = 0 (3.103)

where

a(uu) =

_ l

0

[V ˆγ +M ˆκ] dx

=

_ l

0

[GAs(w

_ + θ) ( ˆ w

_ + ˆθ) + EI θ

_ ˆθ

_ ] dx (3.104)

is the strain energy product. By some simple algebra, it follows that the

homogeneous solutions of the system (3.101) and (3.102) must satisfy the

equations

wIV = 0 w

_ = −θ + EI

GAs

θ

__

. (3.105)

The homogeneous solutions of the first equation are

w(x) = c1 + c2 ξ + c3 ξ2 + c4 ξ3, ξ= x/l (3.106)

and the matching function θ(x) is found to be

θ(x) = c1 · 0 − c2

1

l

− 2 c3

ξ

l

− c4( η

2 l

+

3

l

ξ2), η=

12

l2

EI

GAs

. (3.107)

By a proper choice of the constants ci, the following nodal unit deflections for

the two ends of a beam are found [198]:

298 3 Frames

↓ w1(x) =

1

1 + η

[1 − 3 ξ2 + 2ξ3 + η (1 − ξ)] (3.108)

_ w2(x) = l

1 + η

[−ξ + 2ξ2 − ξ3 − η

2

(ξξ2)] (3.109)

↓ w3(x) =

1

1 + η

[3 ξ2 − 2 ξ3 + η ξ] (3.110)

_ w4(x) = l

1 + η

[ξ2 − ξ3 + η

2

(ξξ2)] . (3.111)

The corresponding rotations are

θ1(x) =

1

1 + η

[−6

l

ξ (ξ − 1)] (3.112)

θ2(x) =

1

1 + η

[1 − 4 ξ + 3ξ2 + (1 − ξ) η] (3.113)

θ3(x) =

1

1 + η

[−6

l

ξ (1 − ξ)] (3.114)

θ4(x) =

1

1 + η

[−ξ (2 − 3 ξη)] . (3.115)

The strain energy products, kij = a(ui,uj ), of the nodal unit deformations

u1 =

_

w1

θ1

_

u2 =

_

w2

θ2

_

u3 =

_

w3

θ3

_

u4 =

_

w4

θ4

_

(3.116)

constitute the stiffness matrix

K = EI

l3 (1 + η)

⎢⎢⎣

12 −6 l −12 −6 l

−6 l (4 + η) l2 6 l (2 − η) l2

−12 6 l 12 6 l

−6 l (2 − η) l2 6 l (4 + η) l2

⎥⎥⎦

. (3.117)

The nodal degrees of freedom ui have the same meaning as in an Euler–

Bernoulli beam (see Fig. 3.20).