3.2 The FE approach

Back

Many engineers consider frame analysis settled for good because all problems

of frame analysis are solved. No one expects new results, more complex

structures are analyzed with finite elements, and the benefit of 1-D models

seems questionable. But the main advantage of 1-D analysis lies in the clear

and descriptive representation of structures because results are immediately

presented in integral form. This gives many engineers the false impression that

• 1-D elements are exact

• 1-D elements are simple

But the more effects that must be considered in the analysis, the more baroque

1-D analysis becomes. It is far easier to program a 2-D finite element than to

consider all the intricacies when a 2-D structure is reduced to a series of 1-D

elements.

Degrees of freedom

In a typical analysis, the x-axis follows the long dimension of the element and

the coordinates y and z point in orthogonal directions. Displacements, and

rotations, forces and moments respectively are defined analogously; see Fig.

3.3. Simple formulations for 1-D elements are only obtained if the description

concentrates on the axis. Preferably resultant internal actions refer to

the centroids of the cross sections. But different construction stages possibly

mean changes in the location of the centroids, and often the inclination of

the neutral axis changes as well. At every stage the axes can differ in length,

and the normal force and shear forces can point in different directions. Then

the quantities must either be transformed, or one works with average values.

Occasionally when transferring data from the CAD-model to the structural

3.2 The FE approach 271

Fig. 3.3. Internal actions in a beam

model, new axes or surfaces must be generated that will not coincide with the

CAD-model.

To avoid these difficulties it is a good strategy to allow from the start for

eccentricities in the position of the centroid and the shear center with respect

to the axes. And it is also a good approach to refer the resultant internal

actions to the user-provided y and z axes, and not to the centroidal axes of

the cross section, because the position of the principal axes within a tapered

beam may rotate.

In such a beam (see Fig. 3.4), the internal actions can be referred to at

least three different coordinate systems. The simplest choice are the main axes

shown in Fig. 3.4 a. This is the classical treatment, where the normal force

and shear forces follow the dashed line, indicating the longitudinal axis. These

components N and V are also the components referred to when the stresses

are calculated and the reinforcement is designed.

But it would be a more elegant and also more efficient approach if the

longitudinal force and the transverse forces of second-order beam theory which

refer to the axes of the undeformed system, were used (see Fig. 3.4 b and c).

Superposition of the stresses when the centroid changes position is only

easy in case (c) because in the other two cases the position of the centroid

must be known, and in case (a) the inclination of the longitudinal axis must be

known as well. But the superposition of internal actions to calculate resulting

stresses only makes sense if the cross section does not change. Therefore the

longitudinal and transverse forces are better referred to the centroid of the

cross section, as in Fig. 3.4 b, because if an arbitrary point of reference is

chosen as in Fig. 3.4 c, it is much harder to interpret the results. Add to

this that if the internal actions are to be plotted, then besides the bending

moment the transverse forces, the longitudinal force, and the rotation of the

beam axis according to second-order theory must also be known.

At each end of the element lies one node whose displacements are given

in terms of the global coordinate system. It is helpful to transform these

272 3 Frames

Fig. 3.4. Possible

coordinate systems

deformations into the local element coordinate system. In accordance with

Fig. 3.4 b, we choose as internal reference nodes the centroids of the cross

sections at both ends of the element. Each of these nodes can have up to

seven degrees of freedom:

1. the longitudinal displacement ux at the centroid of the cross section

2. local transverse deformations uy and uz as well as rotations ϕy and ϕz

3. the twist ϕx and the angle of twist per unit length ϕ_

x

The displacements and rotations of the two nodes are determined by the global

displacements (uX, uY , uZ) and rotations (ϕX, ϕY , ϕZ) if the position of the

local nodes with respect to the global nodes are known (see Fig. 3.5). While

the rotations are identical for the local and the associated global node, the

displacements must be transformed by taking the eccentricities into account.

The longitudinal displacements at the centroid xc, for example, are

ux0(xc) = ux(xi) + ϕy Δz − ϕz Δy xi = node. (3.1)

To consider effects due to warping, matching transition conditions must be

formulated at corner points and at rigid joints. The choice of components is

arbitrary, insofar as also higher derivatives of the displacements could be used

as nodal values, although these would be difficult to control.

Along the beam axis these displacements are interpolated as follows:

1. Linear interpolation of the axial displacement ux

2. Coupled interpolation of the displacement uy and the rotation ϕz

3. Coupled interpolation of the displacement uz and the rotation ϕy

4. Coupled interpolation of the twist ϕx and the angular twist per unit length

ϕ_ (linear interpolation of ϕx if warping torsion is neglected)

The coupled interpolation can be done with cubic splines, in which case the

rotations and the warping effects respectively are simply the derivatives of the

3.2 The FE approach 273

Fig. 3.5. Position of the

centroidal axis

displacements and the rotations. Then it is convenient to choose the shape

functions of a beam

H1 = (1 − ξ2) (1+2ξ) H2 = l (1 − ξ2) ξ (3.2)

H3 = ξ2(3 − 2 ξ) H4 = −l (1 − ξ) ξ2 ξ = x/l (3.3)

so that

uz(x) =

_

i

Hi(x) · uzi ϕy =

_

i

H

_

i(x) · uzi (3.4)

uy(x) =

_

i

Hi(x) · uyi ϕz =

_

i

H

_

i(x) · uyi . (3.5)

Here uyi and uzi are nodal displacements in the local coordinate system at

the two ends of the element.

But the coupling can also be done by considering the shearing strains

according to Timoshenko:

θy = Vz

GAz

, ϕy = uz,x +θy , (3.6)

θz = Vy

GAy

, ϕz = uy,x +θz . (3.7)

The displacements within a cross section are then obtained by a product

approach,

uj =

_7

i=1

Nij(y, z) · ui(x) (3.8)

or if we substitute for ui(x) the interpolating functions

uj =

_2

k=1

_7

i=1

Nij(y, z) · Hik(x) · uik . (3.9)

In general the following functions are used to model the displacements and

higher-order terms within the cross section:

274 3 Frames

ux(y, z) = ux0 + ϕy · (z − zs) − ϕz · (y − ys)

+ Uw · (θx − ϕ

_

y ϕz + ϕy ϕ

_

z) + Uy · θy + Uz · θz + Uw2 · θt2 ,

(3.10)

uy(y, z) = uy0 − ϕx · (z − zm) − 1

2

(ϕ2

x + ϕ2z

) · (y − ym) , (3.11)

uz(y, z) = uz0 − ϕx · (y − ym) − 1

2

(ϕ2

x + ϕ2

y) · (z − zm) . (3.12)

Here ux0, uy0 and uz0 are the displacements of the centroid, ys, zs are the

cross-sectional coordinates of the centroid, and ym, zm are the coordinates of

the shear center.

The first three terms of the longitudinal displacement are determined

by the requirement that the cross section remain flat, in agreement with

Bernoulli’s hypothesis. Next the unit warping functions Uw, Uy, Uz are added,

and finally terms (Uw2, θt2), which take into account deformations caused by

shear forces and secondary torsional moments. The latter add complex patterns

of warping, which in general are not easy to investigate.

Displacements orthogonal to the axis are simple translations and rotations

(rigid-body motions), that is, it is assumed that the cross section keeps

its shape. To allow for changes in the cross section would require either a

sophisticated 3-D model of flat shell elements—think of a double T beam—or

the single terms must be decoupled.

The strains are the derivatives of the displacements. If the higher-order

terms coming from second-order theory are neglected, the stresses become

σx = E εx = E u,x = E [u,x +ϕy,x (z − zs)

ϕz,x (y − ys) − zs,x ϕy + ys,x ϕz +

_

(Ui,x · θi + Ui · θi,x ) x] ,

(3.13)

τxy = Gγxy = G[ux,y +uy,x ] = G[(uy0,x −ϕz)

+

_

(Ui,y · θi) − (z − zm) ϑ,x +zm,x ϑ] , (3.14)

τxz = Gγxz = G[ux,z +uz,x ] = G[(uz0,x −ϕy)

+

_

(Ui,z · θi) − (y − ym) ϑ,x −ym,x ϑ] , (3.15)

σy = σz = τyz = 0. (3.16)

Note that the last three stress components vanish, because the movements

orthogonal to the axis are assumed to be only rigid-body motions.

When the derivatives are calculated, the product rule ensures that terms

appear which are absent in a prismatic member, because in such members the

position of the centroidal axis and the axis of the shear center do not change.

In a tapered beam these terms should not be neglected, because the element

will benefit from their inclusion.

3.2 The FE approach 275

With regard to the shear stresses two approaches are possible. If Timoshenko’s

and Mindlin’s approach is applied, the terms in the first bracket

effect a constant shear distribution, and the cross section remains flat.

A shear distribution which instead satisfies the equilibrium conditions locally

as well (differential equation) is determined by the warping of the cross

section. But this approach can also be used for the classical beam (if the first

bracket becomes zero). The shear stresses then simply depend on the unit

warping deformations, which must be scaled in such a way that equilibrium

is satisfied with regard to the total shear force in the cross section.

In frame analysis the focus is on the resulting stresses, that is, the so-called

internal actions

N =

_

A

σx dA = EA(u,x −zs,x ϕy + ys,x ϕz) + EAz ϕy,x −EAy ϕy,x

(3.17)

My =

_

A

z σx dA = EAz (u,x −zs,y ϕy + ys,x ϕz) + EAzz ϕy,x −EAyzϕy,x

(3.18)

Mz =

_

A

y σx dA = EAy (u,x −zs,x ϕy + ys,x ϕz) + EAyz ϕy,x −________EAyyϕy,x

(3.19)

Vy =

_

A

τxy dA , Vz =

_

A

τxz dA (3.20)

Mt =

_

A

[(y − ym) τxz − (z − zm) τxy] dA (3.21)

where

EA =

_

E dA EAzz =

_

E z2 dA (3.22)

EAy =

_

E ydA EAyy =

_

E y2 dA (3.23)

EAz =

_

E z dA EAyz =

_

E yz dA. (3.24)

The centroid of the cross section is the point at which the area integrals EAy

and EAz vanish. In addition, the unit warping functions Ui are scaled in such

a way that their contributions to the first three integrals vanish. The crosssectional

centrifugal moment EAyz should be retained in any case, because a

transformation of the cross section to the principal axes is neither appropriate

nor always possible in the presence of shear deformations.

Cross sectional values

The cross-sectional values can be easily determined if the cross section is

polygonal (see Fig. 3.6)

276 3 Frames

Fig. 3.6. Cross section

A =

_n

i=0

1

2

(zi+1 + zi) (yi+1 − yi) (3.25)

Iy =

_n

i=0

1

12

(zi+1 + zi) (yi+1 − yi) (z2

i+1 + z2

i ) . (3.26)

Sometimes a correction term in the form of an effective width is added if

Bernoulli’s hypothesis seems too crude an approximation. But these correction

terms depend on what is to be corrected, that is, they are different if crosssectional

values and thereby the stiffnesses are to be modified, or if they

modify terms which are relevant for the design. In such situations one must

be careful to group the correct position of the centroidal axes with the correct

cross-sectional values. For example it is a conventional approach in the design

of prestressed concrete to derive the internal actions for a cross section, which

is different from the approach that is later used for design purposes.

Because shear stresses do not appear in the constitutive equations of an

Euler–Bernoulli beam, they are usually introduced by formulating an equilibrium

condition:

τ,s = σ,x, τ= V Z

I b

. (3.27)

But unfortunately this formula has many drawbacks:

• The shear force V is only correct if the normal force is constant and the

beam has a constant cross section (no tapered beam).

• The section modulus Z is only correct if the cross section is simply connected.

• Shear stresses must especially not be constant across the width of the

beam.

• Eventually I must be replaced by Swain’s formula for skew bending.

Mainly the second remark indicates the basic problem namely that the force

method is not very appropriate for computer programs.

The natural choice would be a displacement method in which the primary

variable is the warping of the cross section. The constitutive equation for this

3.2 The FE approach 277

model is either the principle of minimum potential energy, or equivalently the

set of equations that formulates the equilibrium conditions:

τxy = G(w,y −z θx,x ) (3.28)

τxz = G(w,z +y θx,x ) (3.29)

GΔw = G(w,xx +w,yy ) = −σx,x (3.30)

τxy ny + τxz nz = 0 on the edge . (3.31)

This problem can either be solved in one step or divided into four subproblems:

• The primary torsion problem: dθ/dx = warping ; σx = 0

• The two shear problems: dθ/dx = 0; σx = from the shear force V

• The secondary torsion problem: dθ/dx = 0 ; σx = from the warping moment

In the case of thin-walled cross sections where the stresses in thickness direction

are nearly constant, the problem is easily solvable. Otherwise the Poisson

equation must be solved numerically.

The result of such an advanced analysis is a precise shear stress distribution,

which facilitates the design. By employing an equivalence principle, the

internal strain energy can also be used to calculate the shear deformations,

and the shear stress distribution also defines the unit warping functions Ui.

Stiffness matrix

If the interaction between the mixed terms of the stresses and the unit warping

functions are neglected, the strain energy of the element can be written

Πi =

1

2

_

E ε2 dV =

1

2

_

(EA

_

(u

_

0)2 − 2 u

_

0 [ϕy z

_

s

ϕz y

_

s]

_

+ EA

_

ϕ2

y (z

_

s)2 + ϕ2z

(y

_

s)2 − 2 ϕy z

_

s ϕz y

_

s

_

+ EIy (ϕ

_

y)2 + EIz (ϕ

_

z)2 − 2 EIyz ϕ

_

y ϕ

_

z) dx . (3.32)

This integral yields a stiffness matrix that not only takes into account the

normal displacements and the bending stiffness, but also the stiffening effect

of a tapered beam.

If the shape functions satisfy the differential equations, the stiffness matrix

is exact. If this is not the case, as for example in the presence of an elastic

foundation or if the beam is tapered or second-order theory effects must be

considered, the beam must be subdivided into small elements to improve the

approximation.

In the case of a tapered beam (see Fig. 3.7) the results obtained with different

elements are displayed in the following table.

278 3 Frames

Fig. 3.7. Tapered beam

w(mm) Ne (kN) Nm (kN) Mye (kN m) Mym (kN m)

Inclined axis

* 1 element 0.397 −80.5 −78 −73.58 31.91

* 8 elements 0.208 −46.3 −43.8 −94.87 19.17

** 1 element 0.172 −39.8 −37.3 −93.65 22

** 8 elements 0.206 −45.8 −43.3 −95.02 19.14

Horizontal axis

1 element 0.168 −37.9 −37.9 −93.01 22.52

8 elements 0.204 −44.2 −44.2 −94.85 19.1

Here * indicates that the shape of 1/EI was interpolated, while ** indicates

that the shape of EI was interpolated and e = end and m = middle refer to

the position of the cross section.

Shear stresses and shear deformations

Shear deformations are neglected in classical beam theory. It would not be

correct to simply add strain energy terms that account for the equivalent shear

cross sections, because these add stiffness to the system but not flexibility.

This can only be done in a formulation based on the complementary energy

approach.

The simplest approach is it to reduce the bending stiffness in such a way

that the deformations become the same. But this technique is limited to prismatic

members, and would only achieve the intended effect for one particular

load case. Also an extension to tapered beams is not so easy.

Timoshenko proposed to do the coupling of the rotations and the derivative

of the deflection with a Lagrange multiplier. Hence the equations

M = −EI w

__

, V= −EI w

___

, (3.33)

are replaced by the two equations

3.2 The FE approach 279

M = −EI ϕ,x, V= −GAθ = GA(ϕ − w,x ) . (3.34)

If the rotations and displacements are interpolated with linear polynomials,

the bending moment is constant and the shear force is linear. This leads to

two new problems that need attention:

• For large values of GA locking becomes a problem. A possible remedy is

the introduction of a so-called Kirchhoff mode, which guarantees that at

one point on the axis the derivative M_ equals the shear force V. This

condition relates u to ϕ, which enables one to modify the shape functions

accordingly. If Hughes’ approach is adopted [121] then

M = −EI ϕ,x (3.35)

V = −GAθ = GA[ϕi + ϕj

2

− wj − wi

L

] . (3.36)

• In addition, the element can no longer model simple problems correctly, as

for example a cantilever beam that carries a single force at the free end. The

rotations and the internal actions are correct, but the deflection is not. To

overcome this problem the beam must be subdivided into many elements,

which is hardly a sound approach from the standpoint of the user. In this

case it helps to supplement the rotations with a nonconforming quadratic

function ϕm so that—considering the Kirchhoff-condition—it follows

ϕ = ϕi (1 − ξ) + ϕj ξ + ϕm (4 ξ (1 − ξ)) , (3.37)

M = −EI

L

[(ϕj − ϕi) + 1.5 ϕm (8 ξ − 4)] , (3.38)

V = −GAθ = −GA[ϕi + ϕj

2

+ ϕm − wj − wi

L

] . (3.39)

The function produces exactly the missing linear variation of the bending

moment. The shear force stays at its maximum value, which requires a

correction factor 1.5 for the bending moment.

But even such an advanced model as a Timoshenko beam cannot accommodate

all effects. In numerical tests the results depended for some crosssectional

shapes on the orientation of the coordinate system, because if the

shear force vector as well as the vector of the shearing strains is transformed,

the tensor of the inverse shear areas is obtained

_

θy

θx

_

=

_

1/GAy 1/GAyz

1/GAyz 1/GAz

__

Vy

Vz

_

. (3.40)

On the one hand the introduction of a mixed shear area 1/Ayz removes the

inconsistency in the results; on the other hand, this additional area cannot

simply be added to the Timoshenko beam as an additional stiffness, because

normally it is infinite, and it would lead to a coupling of the bending moments

about the two axes.

There is no easy way to account for this matrix in the classical beam

element. There are only two possibilities:

280 3 Frames

Fig. 3.8. Shear stresses

• The first is to use a flexibility matrix. This may be established via numerical

integration of the differential equations. Extra shear deformations

are simply added and eventually included when this matrix is inverted to

derive the stiffness matrix.

• The other solution based on the Timoshenko beam operates with the inverse

of the above matrix, which yields the following term for the potential:

[Vy, Vz]

1

1 − a

_

Ay −b

b Ay

_ _

Vy

Vz

_

a = Ay Az

Ayz

b = Ay Ay

Ayz

. (3.41)

In the Timoshenko model the shear stresses are constant within the cross

section, while they have a parabolic shape in rectangular cross sections. The

shear areas for the deformations and for the maximum stresses are therefore

in general not identical:

θ = V

GA

(A = 0.833 · b · h for a rectangle) ,

τ = V

A

(A = 0.666 · b · h for a rectangle) .

One aspect deserves our particular attention: normally the prismatic bar has

a constant cross section and in calculating the shear stresses it is assumed

that they vanish at the edge of the cross section. But in the case of a tapered

beam this is not correct; see Fig. 3.8.

Shear stresses appear at the edge of a tapered beam (see Fig. 3.9). In a

plane orthogonal to the reference axis the distribution of the shear stresses is

asymmetric, (Fig. 3.10 a), while in a plane orthogonal to the centroidal axis

the shear stress distribution is more harmonic (Fig. 3.10 b).

In the cross section x = 1.0 the transverse force is 40 kN, hence the shear

force is 39.95 kN. If the cross section were rectangular with a height of 93.75

cm, the shear stress would be 63.92 kN/m2. Because the bending moment in

this cross section is 52.6 kN m, the shear force M/d · tan α is reduced by 5.6

kN, so that the real value is 54.9 kN/m2. The FE shear stress is 53 kN/m2 at

the center and 30 kN/m2 at the edge. Hence this reduction of the shear force

comes close to the maximum value, but the distribution over the cross section

3.2 The FE approach 281

Fig. 3.9. Shear stresses τxy

Fig. 3.10. Shear stresses: a) τxy in a vertical plane, b) τnq in a plane orthogonal

to the centroidal axis

is only an approximation, so it is questionable whether the calculated principal

tensile stresses or the maximum von Mises stresses are always correct.

Influence of the shear center, warping torsion

If the influence of warping torsion is to be considered as well (and assuming

that the shear center axis changes its location) so many terms contribute to

the strain energy integral that no complete analysis seems to have been done

yet. In the simpler case that the reference axis coincides with the rotational

axis, use could be made of the following formula for the strain energy:

Πi2 =

1

2

_ l

0

(ECM (ϑ

__)2 + GI (ϑ

_)2

+ N [2 ϑ

_

zm v

_

m + 2ϑ

_

ym w

_

m + (v

_

m)2 + (w

_

m)2 + im (ϑ

_)2]

+ My [−2ϑw

__

m + rMy (ϑ

_)2] +Mt [2ϑv

__

m + rMz (ϑ

_)2]

+ Mb [rMw (ϑ

_)2] +Mt [v

_

m w

__

m

− v

__

m w

_

m]) dx . (3.42)

Of course these terms are to be supplemented with additional terms that

depend on the load.

282 3 Frames

With regard to the torsional moment, two components are now to be

considered. The total moment can be split into a Saint-Venant part and a

part that stems from secondary torsion:

Mt = Mtv +Mt2 = GIT ϑ

_ −ECM ϑ

___

. (3.43)

Evidently a cubic polynomial is now needed to model the twist of the longitudinal

axis. This is best achieved by introducing two additional degrees of

freedom for the warping of the cross section. Then the same standard thirddegree

Hermite polynomials are also used to model the twist. As in the case

of the shear force, one obtains in each element a constant and a stepped distribution

of the secondary torsional moment. These formulas are also valid if

no warping occurs (CM = 0).

With regard to the shearing deformations due to warping, the same observations

are true as for the shear deformations. They can either be incorporated

directly, or one can use an approach similar to the approach used in a Timoshenko

beam.

The general beam element

In the most general case—if arbitrary loads are to be applied, if the cross section

can have any shape and support conditions can be rather complicated—

use can be made of transfer matrices. The transfer matrix for a beam with

constant values of EI, GAs and a constant load p is

⎢⎢⎣

w

w_

M

V

⎥⎥⎦

=

⎢⎢⎢⎢⎣

1 −x x2

2 EI

x3

6 EI

− x

GAs

0 1 − x

EI

− x2

2 EI

0 0 −1 −x

0 0 0 −1

⎥⎥⎥⎥⎦

⎢⎢⎣

w0

w_

0

M0

V0

⎥⎥⎦

+

⎢⎢⎢⎢⎢⎣

x4

24 EI + x2

2GAs

− x3

6 EI

−x2

2

−x

⎥⎥⎥⎥⎥⎦

p .

(3.44)

From an abstract point of view these equations can be interpreted as the

result of a direct or numerical integration of the differential equations. Numerical

integration is best done with Runge–Kutta methods. In simple cases

the exact solution is obtained after just one step and the additional effort is

minimal. In the worst case a stiff system of differential equations is obtained,

and then the beam must be subdivided into many elements to control the numerical

solution. The advantage of transfer matrices over variational methods

is that by using transfer matrices a flexibility matrix can be derived, which

corresponds to the principle of minimum complementary energy. This is an

important property when it comes to a correct assessment of the shear deformations.

An additional advantage is that the differential equation is easier to

program.

3.2 The FE approach 283

Resistance factor design

There are other areas where a generalization of the current models opens up

new possibilities. Some building codes allow to work in cross sections with

plastic limit loads. Here too we have in principle the possibility of using finite

elements and to apply elasto-plastic laws, or the whole analysis can be done

in the cross section. In the case of thin-walled cross sections a complete description

of the interaction of all resultant stresses is feasible [143]. If a general

nonlinear procedure is employed, then the analysis must be done iteratively

due to the changing stiffnesses, and the need to handle the interaction in the

partially plasticized zones.

Frame analysis is usually based on the following assumptions:

• linear elastic distribution of the normal stresses and the stresses caused by

torsional warping

• linear elastic distribution of the shear stresses due to shear forces and

torsional moments

• a plane strain state plus warping

• a yield condition

To calculate the linear distribution of shear stresses in the cross section advanced

methods—based either on a displacement or a force method—are necessary

to determine the shear flow.

In the first step, given a certain combination of loads, a linear equivalent

stress distribution can be calculated, and these stresses can then be related

to the yield stresses.

In the second step one might “somehow” reduce these stresses. Besides the

special case where the shear force is taken care of separately, other techniques

are available.

In the third step, the stresses can be transformed by numerical integration

into resultant internal actions, and these can then be used to calculate effective

nonlinear stiffnesses. By an iterative procedure based on the stiffnesses, the

equilibrium of the internal actions ultimately is maintained [141].

To control the yield condition it suffices to check just the normal stress

σx and the shear stresses τxy and τxz, because the shear stresses τyz and the

stresses σy and σz will only be noticeable near the loaded regions, and they

are probably not influenced by changes in the stiffness of the cross section.

For an investigation into the local limit load of such regions, 1-D analysis is

not appropriate. Rather this would require experimental tests or an elaborate

FE analysis with flat shell elements.

With regard to the intended reduction of the stresses, in principle three

methods are conceivable. Common to all three methods is that each stress

point is independent of the neighboring points. Plastic strains that might lead

to additional stresses because the movement in y or z direction is impeded

are therefore ignored, and perhaps this provides an additional safety factor

with respect to experimental results. Formally this is in agreement with other

284 3 Frames

engineering approaches, as for example the Winkler model. We can choose

between three methods:

• Prandtl solution In adopting Prandtl’s yield conditions, we calculate

the plastic strains, which are orthogonal to the yield surface. To this end

[72], [258], first the onset of the development of a plastic zone is calculated

by a uniform reduction according to the first method, and then for the

remaining plastic strain increments an elasto-plastic elasticity matrix is

calculated by considering the elasticity matrix C:

σ =

_

C − q · C · q_

q_ · C · q

_

· ε , q = ∂F

σ

. (3.45)

• Isotropic reduction Shear and normal stresses are relaxed in the same

ratio, so that the equivalent stress just reaches the yield stress:

σ =

_

fy

σv, elastic

_

· σv, elastic, τ =

_

fy

σv, elastic

_

· τv, elastic . (3.46)

• Shear stresses Shear stresses are absorbed fully and therefore normal

stresses are reduced. This is the usual strategy in manual calculations. It

remains unsatisfactory in the presence of strong shear stresses, because it

can lead to situations where an increase in the curvature has no effect on

the system

τ = min

_

√fy

3

, τ elastic

_

, σ= min

#

f2

y

− 3 τ 2, σelastic

$

. (3.47)

With one of these methods, the resulting internal actions can be calculated

by numerical integration.

If these internal actions exceed the actual internal actions, the structure is

safe. If we do an elastic–plastic analysis, we have only to check the slenderness

ratio (or a similar criterion) to pass the design check. Otherwise an iterative

analysis must follow to allow for a redistribution of the forces.

If plastic zones develop in a cross section, the stiffnesses change. To assess

the rearrangement of internal forces, these changes must be considered in the

analysis. In bending-dominated problems, the equations can easily be modified

if the following equation is either solved for the plastic curvature κ0 or used

to calculate a secant stiffness:

_

My

Mz

_

=

_

E Iy E Iyz

E Iyz E Iz

_ _

κy

κz

_

+

_

κy0

κz0

_

. (3.48)

With regard to shear strains, it is important to decide whether a rearrangement

of the shear strains is possible or desirable. In some cases a reduction in

the bending stiffness alone will lead to a reduction in the shear strains. But

in general a static analysis must incorporate the shear deformations, because

3.2 The FE approach 285

Fig. 3.11. Nonlinear analysis of a beam

only then can the shear stiffnesses be reduced. But unlike the situation in

bending problems, there is no simple rule to determine nonlinear shear deformations.

Nevertheless one can simply reduce the shear stiffness according to

the ratio between internal and external shear force, but note that during the

iteration the stiffness can increase again if the shear forces decrease.

Nonlinear analysis

In nonlinear analysis the stiffness of a structure depends on the actual strains

and stresses within the structure. By introducing an equivalent secant modulus

the nonlinear problem can be reduced to an elastic problem but then an

iterative analysis must be performed to find the strains which are compatible

with the internal actions. In the final step of such an iterative analysis it must

be checked whether the internal actions are compatible with the iteratively

determined stiffnesses. We perform this check for the beam in Fig. 3.11, grade

C 20 concrete and grade S 500 reinforcement steel (Eurocode), [180].

If the origin of the system of coordinates does not coincide with the elastic

centroid the matrix which describes the relation between the strains and the

internal actions is fully populated

Nx

My

Mz

⎦ =

EA EAz −EAy

EAz EAzz −EAyz

−EAy −EAyz EAyy

u_

−w__

v__

⎦ . (3.49)

For a definition of the terms EA,EAz, etc. see (3.22). A nonlinear analysis

with an FE program rendered a value of 4.70 cm2 for the reinforcement and

286 3 Frames

Fig. 3.12. σ ε diagram

for the strains the values

εu = −2.0_ εl = 11.12_ εm =

11.12 − 2.0

2

= 4.56_ (3.50)

so that the curvature of the cross sections is

− w

__ =

2.0_+ 11.12_

0.6m

= 21.87 · 10−3 m−1 . (3.51)

The nonlinear analysis was based on the following stress-strain law for the

concrete (see Fig. 3.12)

σc = ε · (1 − ε

4

) · α · fcd α = 0.85, fcd =

20

1.5

= 13.3MN/m2 (3.52)

so that the secant stiffness becomes

E = σc

ε

= (1 − ε

4

) · α · fcd . (3.53)

Both the concrete (b = 30 cm)

EAc :=

_ zu

z=0

(1 − ε

4

) · α · fcd · b dz = 234.5MN (3.54)

and the steel

EAs := σ________s

εs

· As = 20.65MN (3.55)

contribute to the longitudinal stiffness

EA :=

_

E dA = EAc + EAs = 234.5 + 20.65 = 255.15MN. (3.56)

In the same sense the statical moment is the sum of two parts

3.2 The FE approach 287

Fig. 3.13. An FE analysis of a standard frame: a) linear analysis; b) nonlinear

analysis; c) nonlinear analysis + axial displacements; d) plus tension stiffening

EAz :=

_

E z dA = EAz,c + EAz,s = −58.37 + 5.16 = −53.21MN

(3.57)

and the moment of inertia as well

EAzz :=

_

E z2 dA = EAzz,c + EAzz,s = 14.69 + 1.29 = 15.98MNm

(3.58)

so that

288 3 Frames

Fig. 3.14. An FE analysis of standard frames yields the exact values

Nx

My

Mz

⎦ =

255.15 −53.21 −EAy

−53.21 15.98 0

−EAy 0 EAyy

4.56 · 10−3

21.87 · 10−3

0

⎦ =

0

107

0

⎦ . (3.59)

The result agrees with the internal actions, N = 0,My = 107kNm,Mz = 0.

The constant strain εm of the x-axis causes a horizontal displacement

ux =

_ l

0

εm dx = 4.56_· 1.0m = 4.56mm. (3.60)

Often the horizontal displacements are neglected but their influence on the

structural reaction can be considerable and sometimes they produce quite

surprising results.

Today it is standard that FE programs allow to consider shear deformations

in frame analysis. It is hoped that in the future it will be possible to

include the effects of nonlinear axial strains and tension stiffening as well in

the analysis. The influence of these different effects

• linear analysis

• nonlinear analysis

• nonlinear analysis + axial displacements

• nonlinear analysis + axial displacements + tension stiffening

on the bending moment distribution of a concrete frame is depicted in Fig.

3.13.

3.3 Finite elements and the slope deflection method 289

3.15. Any deflection

curve can be split into two

parts