Nonlinear Analysis

Chapter 12 - Nonlinear Analysis

CHAPTER 12 Important issues regarding nonlinear analysis of structures are described.Three types of nonlinearities are introduced with an emphasis on geometrical and material nonlinearities. Nonlinear formulations for one dimensional bars and beams are described, as well as generalization to multi-dimensional problems. The most recent solution techniques are presented. Applications are placed on beams, frames, plates and shells.

Nonlinear analysis

12.1

Introduction.............................................................................................................................................................. 12.1.1 General....................................... ...................................... ..................................... ................................... ... 12.1.2 A nonlinear geometrical problem................................................... ...................................... ...................... 12.1.3 Nonlinear material behaviour........................................................... .................................... ...................... 12.2 Stiffness relationship for beam with axial force........................................ ............................... ............................ 12.2.1 General....................................... ...................................... ..................................... ................................... ... 12.2.2 Comparison of alternative stiffness matrices for lateral deformations of a bem with axial force.............. 12.3 Formulations for nonlinear geometrical behaviour of bars and beams with axial and lateral deformation 12.3.1 General....................................... ...................................... ..................................... ................................... ... 12.3.2 Methods with updated coordinates...................................................... ...................................... ................. 12.3.3 Total Lagrangian formulation for a beam with axial and lateral deformation............................................ 12.3.4 Generalization................................................... ........................................ ..................................... ............. 12.4 Nonlinear material behaviour........................................................ .................................... .................................... 12.4.1 One dimensional case...................................... ...................................... ................................. .................... 12.4.2 Generalization................................................... ........................................ ..................................... ............. 12.4.3 Cyclic plasticity, shakedown and ratchetting............................................. .................................. ................ 12.5 Solution techniques.......................................... .................................... ................................. ................................... 12.5.1 General....................................... ...................................... ..................................... ................................... ... 12.5.2 Load increnmental methods....................................................... ................................... .............................. 12.5.3 Iterative methods........................................................... ..................................... ..................................... ... 12.5.4 Combined methods............................................................ .................................... ............................... ...... 12.5.5 Advanced solution procedures...................................................... .................................... ......................... 12.5.6 Direct integration methods........................................................ ..................................... ............................ 12.6 Applications............................................................................................................................................................. 12.6.1 General....................................... ...................................... ..................................... ................................... ... 12.6.2 Beams and frames.......................................... ...................................... ................................. ...................... 12.6.3 Plane stress, plates and shells........................................... ................................... ................................... .... 12.7 Analysis of accidental load effects............................. ................................. .................................... .................... 12.6.1 General....................................... ...................................... ..................................... ................................... ... 12.6.2 Fires and explosions............................................................ ..................................... ............................... ... 12.6.3 Ship impacts.............................................. ..................................... ...................................... ....................... Appendix A Solution of the differential equation of a beam with axial load................................... ........................ Appendix B General formulation for geometrically nonlinear behaviour................................................. ............. Appendix C Plasticity theory............................................ ................................. ................................ .......................... Reffrences ...................................................................................................................................................................

page 12.2 12.2 12.5 12.14 12.17 12.17 12.18 12.20 12.20 12.23 12.28 12.35 12.36 12.36 12.42 12.43 12.45 12.45 12.48 12.57 12.56 12.58 12.63 12.67 12.67 12.67 12.74 12.84 12.84 12.86 12.88 12.94 12.99 12.106 12.116

12.1


12 Nonlinear Analysis 12.1 Introduction 12.1.1 General Linear versus nonlinear nonlinear analysis

Structural analysis – including the finite element method – is based on the following principles:   

Equilibrium(expressed by stresses) Kinematic compatibility (expressed by strains) Stress-strain relationship

So far, the analysis has been based on the assumptions that  

Displacements are small The material is linear and elastic

When the displacements are small, the equilibrium equations can be established with reference to the initial configuration. Moreover, this implies that the strains are linear functions of displacement gradients (derivatives). The linear elastic stress-strain relationship corresponds to Hooke’s law. The relationship between load and displacement for structures with nonlinear behaviour may be as shown in Fig. 12.1. When the ultimate strength of structures that buckle and collapse is to be calculated, the assumptions about small displacements and linear material need to be modified. If the change of geometry is accounted for, when establishing the equilibrium equations and calculating the strains from displacements, a geometrical nonlinear behaviour is accounted for. Various examples are given in Fig. 12.1. In section 12.1.2, a quantitative example is completely washed out. Analogously, material nonlinear behaviour is associated with nonlinear stress-strain relationship. An example is given in Section 12.1.3. Finally, nonlinearity may be associated with the boundary condition, i.e. when a large displacement leads to contact. Boundary non-linearity occurs in most contact problems, in which two surfaces come into or out of contact. The displacements and stresses of the contacting bodies are usually not linearly dependent on the applied loads. This type of non-linearity may occur even if the material behavior is assumed linear and the displacement are infinitesimal, due to the fact that that the size of the contact area is usually not linearly dependent on the applied loads, i.e. doubling the applied loads does not necessarily produce double the displacement. If the effect of friction is included in the analysis, then slick-slip behaviour may occur in the contact area which adds a further non-linear complexity that is normally dependent on the loading history. 12.2


12 Nonlinear Analysis 12.1 Introduction 12.1.1 General Linear versus nonlinear nonlinear analysis

Structural analysis – including the finite element method – is based on the following principles:   

Equilibrium(expressed by stresses) Kinematic compatibility (expressed by strains) Stress-strain relationship

So far, the analysis has been based on the assumptions that  

Displacements are small The material is linear and elastic

When the displacements are small, the equilibrium equations can be established with reference to the initial configuration. Moreover, this implies that the strains are linear functions of displacement gradients (derivatives). The linear elastic stress-strain relationship corresponds to Hooke’s law. The relationship between load and displacement for structures with nonlinear behaviour may be as shown in Fig. 12.1. When the ultimate strength of structures that buckle and collapse is to be calculated, the assumptions about small displacements and linear material need to be modified. If the change of geometry is accounted for, when establishing the equilibrium equations and calculating the strains from displacements, a geometrical nonlinear behaviour is accounted for. Various examples are given in Fig. 12.1. In section 12.1.2, a quantitative example is completely washed out. Analogously, material nonlinear behaviour is associated with nonlinear stress-strain relationship. An example is given in Section 12.1.3. Finally, nonlinearity may be associated with the boundary condition, i.e. when a large displacement leads to contact. Boundary non-linearity occurs in most contact problems, in which two surfaces come into or out of contact. The displacements and stresses of the contacting bodies are usually not linearly dependent on the applied loads. This type of non-linearity may occur even if the material behavior is assumed linear and the displacement are infinitesimal, due to the fact that that the size of the contact area is usually not linearly dependent on the applied loads, i.e. doubling the applied loads does not necessarily produce double the displacement. If the effect of friction is included in the analysis, then slick-slip behaviour may occur in the contact area which adds a further non-linear complexity that is normally dependent on the loading history. 12.2


Fig. 12.1c shows a typical contact problem of a cylindrical roller on a flat plane. Initially the contact is at a single point, and then spreads as the load is increased. The increase in the contact area and the change in the contact pressure are not linearly proportional to the applied load. Another example is shown in Fig. 12.1c where the tip t ip of the cantilever comes into contact with a rigid surface.

a) Response of a thin plate/shell (e.g. due to water pressure or explosion pressure)

b) Response of a column subjected to lateral and axial compressive lo ad

c) Representation of contact Figure 12.1 Typical nonlinear geometrically behaviour. 12.3


Table 12.1 Comparison of linear and non-linear analysis [Becker, 2001]

12.4


Reasons for nonlinear nonlinear stress analysis

There are several areas where nonlinear stress analysis may be necessary (Moan et al., 2002): 



 

Direct use in design for ultimate and accidental collapse limit states. Modern structural design codes refer to truly ultimate failure modes and not only first yield and analogous modes. Use in the assessment of existing structures whose integrity may be in doubt due to (a) visible damage (crack, etc.) concern over corrosion or general ageing. The above will largely relate to the ultimate limit state because, in many cases, the serviceabil ity limit state will already have been exceeded and yet key question still remain such as: Is the structure safe? Should it be repaired and if so, how will any proposed strengthening work? Can it be kept in service for a little time longer? Use to help to establish the causes of a structural failure. Use in code development and research: (a) to help to establish simple ‘code-based’ methods of analysis and design, (b) to help understand basic structural behaviour and (c) to test the validity of proposed ‘material models’.

With the new generation of inexpensive yet powerful computers, solution cost is no longer the major obstacle it has been. However, the complexity of nonlinear stress analysis still remains to provide the ‘expert’ as well as the unwary novice with many headaches. Nonlinear analyses are applied in all the ways mentioned above. However, a significant increase in the use of nonlinear stress analyses in the assessment of existing structures is envisaged and eventually in the direct design of more routine structures.This will occur as hardware becomes cheaper and faster and software becomes more robust and user-friendly. It will simply become easier for an engineer to apply direct analysis rather than code based charts. However, problems will arise because the latter often include ‘fiddle factors’ relating to experience, uncertainty, etc. The advent of more computer-based analysis procedures will lead to the need for a ‘surrounding’, probably computer based, ‘code’ to incorporate the ‘partial factors’ including those factors (often now hidden) relating to the degree of uncertainty of the analysis. The analysis would have to be directly embedded in a statistical reliability framework.

12.1.2 A nonlinear geometrical problem

Geometrical nonlinearity may be illustrated by the bar system shown in Fig. 12.2a (Bergan and Syvertsen, 1978). See also Crisfield (1991).

12.5


a) Geometry

b) Deformation and equilibrium equilib rium for small displacements ( r)

c) Deformation and equilibrium for large displacements ( r) Figure 12.2 Two-bar systems Linear model

When r is small compared to h, the axial strain in the bar is ε =

r sin α o  / cos α o

=

r sin α o cos α o



(ε is positive when the bar shortens) and the axial force becomes: S = EAε =

EA 

sin α o cos α o ⋅ r

S is positive in compression. Equilibrium as referred to the initial geometry, gives R = 2 S sin α o

or

=

2 EA 

sin 2 α o cos α o ⋅ r

(12.1)

R = Kr

where K =

2 EA 

sin 2 α o cos α o

The stiffness K is constant, implying a linear relationship between force and displacement. 12.6


If the angle αo is small (αo << 1) sinαo ∼ αo, cosαo ∼ 1 i.e. R =

2 EAα o2 

r

(12.2)

Nonlinear model (large deformations)

The true axial shortening for finite (not small) value of r is

Δ=



cos α o

−



cos α

and the strain is (positive in compression) ε =

Δ  / cos α o

=1 −

cos α o cos α

Equilibrium for the deformed truss is:

⎛ ⎝

R = 2 S sin α = 2 EAε sin α = 2 EA sin α ⎜1 −

cos α o ⎞

⎟

(12.3)

cos α ⎠

By introducing h − r

sin α = 

+ (h − r ) 2

2



, cos α = 

2

+ (h − r ) 2

, cos α o =

 

2

+ h2

Eq. (12.3) may be written as R =

⎞ ⎛ ⎜ − 1⎟ ⎜⎜ ⎝ r ⎠ ⎝

2 EA ⎛ h 

 

2

+ ( h − r ) 2

−

⎞ ⎟ r 2 2 ⎟  + h ⎠ 

(12.4)

or R = K (r ) ⋅ r

The stiffness now depends upon the displacement r and the force-displacement relationship (12.4) is nonlinear . For small angles α and αo: h − r sin α ≈ α ≅ tgα = 

1

cos α ≈ 1 − α 2 cos α 0

1 ⎛ h − r ⎞

2

≈1 − ⎜ ⎟ 2 ⎝  ⎠ 2 1 2 1 ⎛ h ⎞ ≈1 − α 0 ≈1 − ⎜ ⎟ 2 2 ⎝  ⎠ 2

12.7


Introducing these expressions into Eq. (12.4) yields

⎛ h − r ⎞ ⎛ h − 1 r ⎞ ⎜ ⎟⎜ ⎟ 2 EA ⎝   ⎠ ⎝  2  ⎠ 2 EA ⎛ h r ⎞ ⎛ h 1 r ⎞ R = r ≈ ⎜ − ⎟ ⎜ − ⎟ r 2 2  ⎠   ⎝   ⎠ ⎝  1 ⎛ h ⎞ 1− ⎜ ⎟ 2 ⎝  ⎠ 2

when assuming that

By introducing

R

or

=

h



≈ α o , the following equation results:



2 EA

⎛ h ⎞ << 1 . ⎜ ⎟ ⎝  ⎠

⎛ ⎝

α o2 ⎜1 −

r ⎞ ⎛ r ⎞ ⎟ ⎜1 − ⎟ r h ⎠ ⎝ 2h ⎠

(12.5)

R = K (r )r

(12.6)

where K (r ) =

2 EA

=

2 EA





⎛ ⎝

α o2 ⎜1 −

α o2 +

r ⎞ ⎛ r ⎞ ⎟ ⎜1 − ⎟ h ⎠ ⎝ 2h ⎠

EA 

⎛ r ⎞ r α o2 ⎜ − 3 ⎟ = K o + K g ⎝ h ⎠ h

(12.7)

The first term of Eq. (12.7) is the linear stiffness term; see Eq. (12.2), while the second term is a correction due to nonlinear geometrical effects. The stiffness relationship is a third degree polynomial, plotted in Fig. 12.3a. Eq. (12.7) does not always give unique solutions for a given load. Fig. 12.3b shows that three equilibrium points (A, B, C) may correspond to a given load level. This also means that it would not be possible to follow the load-displacement curve by increasing the load. Actually when the point D is reached, the solution will jump to E, which is a stable equilibrium condition. This phenomenon is called snap-through.

a) Load-deflection

b) Possible equilibrium conditions

Figure 12.3 Load-deflection characteristics of two-bar system. 12.8


Eq. (12.5) expresses a formulation of equilibrium between external loads (R ext = R) and internal loads R int = K(r) r. The stiffness K(r) is denoted secant stiffness . The Eq. (12.6) can be solved, i.e. r can be determined for a given R, analytically. However, in general this is not possible and iterative methods need to be used. Then it is convenient to express the equilibrium condition, Eq. (12.5) on a differential form: dR =

d dr

( K (r )r )dr = K I dr

(12.8)

where K I (r ) =

d

( K (r ) r ) dr is denoted the tangent stiffness or incremental stiffness . The formulation (12.8) allows a solution by initial value problem or incremental methods. Such a method may be combined with or replaced by iterative methods. This topic is discussed in Section 12.5. The incremental stiffness K I( r) for the problem defined by Eqs (12.6 – 12.7) is

⎛ ⎛ r ⎞ ⎛ r ⎞ ⎞ ⎜ ⎜1 − h ⎟ ⎜1 − 2h ⎟ r ⎟ dr ⎝ ⎝  ⎠⎝ ⎠ ⎠ 2 2 EA 2 ⎛ r 3 ⎛ r ⎞ ⎞ = α o ⎜1 − 3 + ⎜ ⎟ ⎟ ⎜  h 2 ⎝ h ⎠ ⎟⎠ ⎝

K I ( r ) =

=

2 EA

α o ⋅

2 EA 

d

α o2 +

6 EA 

α o2

(12.9)

⎛ r − 1 ⎞ = K + K ⎜ ⎟ o G h ⎝ 2h ⎠ r

Here, K o is the linear (initial) stiffness, while K G represents the change in incremental stiffness due to change of geometry by deformation. It is sometimes called geometric stiffness. Stiffness concepts

The different stiffness concepts (K o, K G, K I, K) are shown in Fig. 12.4.

K o- linear (initial) stiffness K - secant stiffness K I - incremental (tangent) stiffness K G-geometrical stiffness (negative in the fig.)

a) At a stable point

b) At an unstable point Figure 12.4 Stiffness definitions. 12.9


As illustrated in Fig. 12.4b the incremental stiffness in an unstable point is zero, or dR = K I dr = ( K o

+ K G )dr = 0

This condition is analogous to the stability criterion in case of linearised buckling, which for a beam is: ( K − λ P Kσ ) r = 0 as addresses in Chapter 2.9.4. This means that buckling occurs when the “load factor” equals the elastic stiffness, K.

λP is large enough that λPK G

Further comments about the two-bar problem

Assume that the bars initially have an axial force, S0 (positive in compression). When an external load, R is applied, the total axial force in each bar is denoted by S. The axial force imposed due to the load, R is then S – S 0, and the following equation applies S − S0

= EAε = Aσ

(12.10)

And the following equation replaces Eq. (12.1) R = 2 S sin α

=

2 EA 

sin 2 α 0 ⋅ cos α 0 ⋅ r + 2S 0 sin α

(12.11)

where the first term is as derived before. In particular it is noted that sin α

≈ α ≈

h − r 

and

cos α

≈ cos α 0 ≈ 1.0

K(r) in Eq. (12.7) becomes K (r ) = K 0

+ K g + Ks

(12.12)

Where the additional term is Ks

= 2S0

h−r 

r

= 2S 0α 0 (1 − ) h

(10.12a)

K I in eq. (12.9) becomes K I = K 0 + K G + K σ

(12.13)

where 12.10


K σ

=−

2 S 0 

(12.13a)

These terms constitute  the linear stiffness, K 0  the initial displacement(geometric) stiffness, K G  the initial stress stiffness, K σ K o is familiar from small displacement structural analysis. K G reflects the effect of changing deplacement on the stiffness and K σ is the effect on the stiffness of the initial member forces. Initial means existing prior to further displacement. K 0 and K σ are essential to any linear stability analysis wherein the bucling load is found which gives an internal force (e.g. N ) such that the stiffness K 0 + K σ is singular. K σ is also called the geometric stiffness because it represents the change in the forces maintaining the structure (e.g. the bar) in equilibrium, which results from a rotation, with the internal forces (e.g. S0) remaining constant. It is interesting to consider the special care that h = 0, with a pretension force S 0

= −S 0

(Where S 0 is positive in tension). Noting that α 0 = h /  the governing equation, Eq. (12.5), with account of the geometric stiffness, Eq. (12.13 a) becomes

⎡ EA ( r ) 2 r ⎤ + 2S 0 r R = ⎢ ⎥ ⎣



⎦





which describes the force-displacement relationship for a pretensined cable structure. It is observed that the incremental stiffness has two terms – one initial deplacement (K G) and one initial stress term, namely: KG

=(

3 EA 

3

r ), K σ

=2

S 0 

R(N)

R, r

800

S 0

S 0

Formatted: Font:

(Tension)

= 104

600

a = 2l

Increasing linearity 400

S 0

= 102

a =5000 mm 6

E =200×10

2

N/ mm 2

200

A =0. 25 mm h =0 mm 0

10

20

30

40

50

r,mm

2

S0 =10

4

or 10

N

Fig.12.5 Load-displacement relationship for a pretension cable with lateral load. 12.11

10 pt


It is also interesting to establish K I directly by K I =

dR dr

for the total problem. Based

on the equilibrium equation R = 2 S

h − r 

and, S − S0

h 1 r = EAε ≅ EA ⎜⎛ − ⎟⎞ r ⎝ 2 ⎠

thedrivative

dS becomes : dr

h r = EA ⎡⎢ − ⎤⎥ dr ⎣ ⎦

dS

Then, the tangent stiffness is

⎡ 2S h − r ⎤ ⎥  ⎦ dr dr ⎢⎣ dS h − r d h− r =2 + 2S ( ) K I =

dR

dr

=

=

=

d

dr





r ⎞⎛ h − r ⎞ 1 ⎜  −  ⎟⎜  ⎟ − 2S 0  ⎝ ⎠⎝ ⎠

2 EA ⎛ h 

2 EA ⎛ h ⎞

2

⎜⎟ + ⎝ ⎠



2 EA ⎡ hr 

⎢− 2 ⎢⎣ 

r + ⎛⎜ ⎞⎟ ⎝⎠

= K 0 + KG + K σ K 0

=

2 EA 

=

K σ

= −2



S 0 

⎤ S 0 ⎥ −2  ⎥⎦ (12.14)

(12.15a)

⎡⎛ r ⎞2 r ⎤ ⎜  ⎟ ⎢⎜ h ⎟ − 2 h ⎥ ⎝ ⎠ ⎢⎣⎝ ⎠ ⎥⎦

(12.15b)

(12.15c)

2 EA ⎛ h ⎞

K G

2

2

Another problem is concerned with a cable with a pretension force N0 and its own weight and an initial geometry defined by w0.

12.12


N

u N

N 0

N 0 w0 (Initial geometry)

u Fig. 12.6. Load-displacement relationship for a catenary cable. Comments on strain expressions

In the present derivations the engineering strain, ε defined by ε = (  −  0 ) /  0 has been used. As discussed later other strain definition may be applied – much as the Green’s strain, ε G discussed subsequently. However, for the bar or cable discussed above, it follows by definition that εG

=

1



2

−  20 2 0

2

=

1



2

− 0  + 0 0

0

=

for structural materials where typically of the stress is the same, i.e. 

cos α 0

=

1 2

ε (ε

+ 1 +1) = ε +

1 2

ε2

≈ε

ε ≤ 0.003. Also, with small strains the length



cos α

Moreover, the cross-sector area remains the same. Also the engineering stress σ can be assumed to the same as the continued mechanics stress. Cauchy stress and 2nd Piola- Kirchoff stresses in local coordinates .

Generalization

The stiffness expressions (12.16, 12.18) can be generalized to systems with many degrees of freedom K (r )r=R implying K I (r)dr

(12.16a)

= ( Ko +

K G (r)) dr

=

dR

(12.16b)

where R and r are load and displacement vectors, respectively. Also, the various stiffness concepts indicated in Fig. 12.4, can also be generalized. 12.13


Eq. (12.16a) expresses equilibrium between external loads, R and internal (reaction) forces, Kr. The formulation (12.16a) with secant stiffness, however, is not very practical. The differential formulation (12.16b) may be written on a finite incremental form K I Δr = ΔR

(12.17a)

implying

Δr = K −1ΔR

(12.17b)

ΔR and Δr are corresponding increments in load and displacement, respectively. With a given condition (r, R ), K I can be calculated and the displacement increment Δr due to a load increment, ΔR can be calculated by Eq. (12.17b). 12.1.3 Nonlinear material behaviour

Sofar a linear (elastic) relationship between stress and strain has been assumed. Material tests of metals show that the linearity does not apply when the stress exceeds a level, σP proportionality limit . Above this level a nonlinear elasto-plastic condition prevails. See the two typical relationships in Fig. 12.7. In the example with the mild steel the stress reaches a plateau, yield stress level, until the stress again increases. This phenomenon is denoted hardening.

a) Mild steel

σY

b) High strength steel, aluminium Figure 12.7 Stress-strain curves for metals.

Unloading from a stress condition above σP takes place along a straight line parallel with the initial linear stress-strain relationship, as shown by the dashed line in Fig. 12.7. When the stress is zero, a residual plastic strain, εP remains (see Fig. 12.7a). The nominal or engineering strain is defined as the ratio of the change of length over a given gange length to the original length, l0 as follows ε e

=

l − l0 l0

=

l l0

−1

( 12.18)

12.14


Comment on nominal and true stress and strain

In uniaxial tests used to obtain the stress-strain curve for a given material, the test piece cross-sectional area and length change with the onset of plasticity. Assuming that the load remains constant on the test-piece, the stress will continue to rise, not fall, as the neck develops before final fracture. Therefore it becomes necessary to take account of the ‘current’ cross-sectional area and length, rather than refer to the original (nominal) dimensions of the test-piece. In order to obtain a better understanding of plastic flow behavior, ‘true’ stresses ans strains are defined. To take account of the change in the dimensions of the test-piece, the true stress, is defined as the force over the current (instantaneous) area, A, as follows: σ true

=

F A

σtrue

(12.19)

The true strain (also called natural or logarithmic strain), ε true , is defined as follows: ε true

l

dl

0

l

= ∫l

(12.20)

Where l is the current (instantaneous) length. Using the incompressibility (constant volume) assumption in the plastic flow of metals, i.e. Al = A0l0 , the following relationship between the nominal and current dimensions is obtained: A =

A0l0 l

=

A0

1 + ε 0

(12.21)

Substituting this value of A in equation in equation (12.19), the true stress can be expressed as follows: σ true

=

F (1 + ε 0 ) A0

= σ (1 + ε 0 )

(12.22)

Integrating the expression for the true strain in equation in equation (12.20) gives:

ε true

⎛l⎞ = ln ⎜ ⎟ = ln(1 + ε 0 ) ⎝ l0 ⎠

(12.23)

which is the total (logarithmic) strain between the original and current limits of length. Therefore, the true plastic strain component can be obtained by subtracting the elastic strain from the total strain, as follows:

= ε true − ε elastic = ln(1 + ε 0 ) −

σ 0

(12.24) (1 + ε 0 ) E If the unaxial stress-strain curve in Fig. 12.7, which is derived from the load-extension readings from a uniaxial tension test, is re-plotted with true-stress vs. true strain, as ε plastic

12.15


shown in Fig. 12.8, it can be seen that the curve continues rising beyond the point where necking appears. The curve clearly indicates a ‘strain hardening’ effect, i.e. the material becomes harder as the strain is increased. It shouldbe noted that the true stress-strain curve is strictly only valid up to the onset of necking, since the formation of neck gives rise to a complex state of stress which is no longer uniaxial, i.e. the stress is not simply the force divided by the cross-sectional area.

Fig. 12.8. Nominal and true stress and strain curves. The elastic strain is expressed by Hooke’s law: ε e = σ / E

Figure 12.9 Definitions of material properties. To deal with nonlinear material properties additional terminology needs to be introduced. The stress in point A, may be expressed as: σ = E S ε

(12.25)

where ES is the secant modulus , which depends upon the stress (strain) level. When loading is introduced at A, the change of stress, Δσ can be obtained from

Δσ = E T Δε

(12.26) 12.16


where ET is the tangent modulus . By unloading from A Hookes law applies:

Δσ = E Δε (unloading)

(12.27)

Effect of nonlinear material behaviour on structural behavior

The effect of elasto-plastic behaviour may be illustrated by considering buckling of a simply supported column with length, . The critical buckling stress is then (according to F. Engesser, 1889) obtained by replacing E with the tangent modulus, ET in the expression for the Euler buckling stress: σ Cr =

π 2 E T λ 2

, λ =  / I / A

(12.28)

ET depends upon the stress σ = σCr (see Fig. 12.7). Eq. (12.28), therefore, could be solved by iteration by first assuming σCr , calculating ET corresponding to σ = σCr , calculating a new σCr by Eq. (12.28); use new σCr to calculate ET etc. This example is to illustrate that nonlinear material behaviour (plasticity) affects the behaviour. A more refined approach is needed to calculate the ultimate strength accurately. Moreover, the effect of initial imperfections and residual stresses need to be accounted for. Other non-linear material problems

Material non-linearities are classified into three categories:  



Time-independent behaviour such as the elastic-plastic behavior of metals in which the structure is loaded past the yield point as briefly described above. Time-dependent behaviour such as creep of metals at high temperatures in which the effect of variation of stress/strain with time is of interest and a power law stress-strain relationship is often used. Viscoelastic/viscoplastic behaviour in which both the effects of plasticity and creep are exhibited. Here the stress is dependent on the strain rate and the material behaviour can be represented by a combination of a spring and a dashpot.

12.2 Stiffness relationship for beam with axial force 12.2.1 General

In Section 12.1.2 the stiffness relationship for a two-bar problem was established. This section deals with element stiffness relationships (S = kv) which account for the effect of change of geometry. When such element relationships have been established, the relationship for an entire structural system can be established by the direct stiffness method 12.17


S i = k i v i ⇒ R = Kr Only elasto-plastic behaviour will be pursued herein.

In Section 2.9.3 the stiffness matrices for trusses and beams with axial forces have been derived based on assumed polynominal shape functions and the principle of vertical work. The result was written as S = (k 0 + k σ ) v

(12.29)

(when axial forces in tension are defined as positive) where k 0 – stiffness matrix corresponding to small displacements. k σ - initial stress (geometric) stiffness matrix.

For beams with an axial force it is actually possible to exactly establish the stiffness relationship. This approach is based on the solution of the govering differential equation, as illusrated in Appendix A. It is interesting to compare the approxiomate stiffness matrix obtained in Section 2.9.3 with the exact one obtained in Appendix A.

12.2.2. Comparison of alternative stiffness matrices for lateral deformations of a beam with axial force

In Section 2.9.4 the stiffness relationship involving lateral deformation and rotations for a beam with axial force. The exact stiffness relationship is derived based on the differential equation in Appendix A. It is observed that the exact stiffness matrix, k is the same as for a beam without axial force, but with each stiffness term multiplied by a φi-function, which depends upon axial force. The analytical expressions for these φi-functions are given by Eq. (A.16), and they are plotted as a function of ρ =

P P E

=

P 2

π 2 EI

⎛  ⎜ β = ⎜ ⎝ 2

P EI

=

π 2

⎞ ⎟ ⎠

P⎟

(12.30)

in Fig. 12.10 . Note that PE (Euler load) does not depend on the boundary condition of the beam.

12.18


Figure 12.10 Stability (Livesley) functions. It is observed that all axial force.

φi-functions are equal to 1.0 for β = ρ = 0, i.e. a case with zero

For compressive axial force (positive ρ), all φi’ s (and the stiffness), except decrease with increasing axial force. For tensile axial force, all φi’ s, except increase with the magnitude of P.

φ4, φ4,

The exact stiffness matrix thus established based on linear elstic material behaviour is applied in computer programs like USFOS. The stiffness matrix obtained in Section 2.9.4 is approximate because the interpolation functions do not satisfy the differential equation. However, accurate results may be obtained by using several finite elements. Now, typical terms in the stiffness matrix k ’ NL in Eq. (A.19) may be compared to the corresponding term of k ’L = k ’ – k ’σ . For instance, k ' NL (1.1)

=

12 EI 12 EI 12 EI 12 P ≈ − = − 2 φ ρ [ ] 1 5 3 3 3 π    

(12.31a)

(by visual consideration of the diagram in Fig. 12.10. A formal linearization may be made by considering analytical expressions for the φi-functions.) k ' L (1.1)

=

12 EI 

3

−

6P 5

(12.31b)

which shows that a linearization of k’ NL(1.1) with respect to ρ yields a result very close to k L' (1.1) . The same applies to the other terms. 12.19


12.3 Formulations for non-linear geometrical behaviour of bars and beams 12.3.1. General

Numerical solutions of problems involving geometric non-linearity (GNL) usually attempt to replace the continuous non-linear displacement by a series of linearised increments of displacements. The theories used in the analysis of GNL problems using the FE method are not familiar to most engineers and, unfortunelately, suffer from a multitude og jargon terms – usually terms named after the scientists and mathematicians who first invented them, e.g. Lagrange, Green, Cauchy, Almansi, Piola-Kirchoff, etc. In this section some of features of GNL problems are discussed in order to provide a brief introduction to this type of non-linearity. Further information on all aspects of GNL analysis using FE can be found in the textbooks by Crisfield [1991, 1997], and Zienkiewicz and Taylor [1991]. Classifications of GNL problems

Most engineers associate GNL problems with ‘large’ displacements effects such as buckling. However, GNL problems can also involve small displacements (in fact, GNL theory often predicts less displacement than the corresponding linear loaddisplacement theory). In order to classify GNL problems, it is best to focus on whether strain is assumed small or large, as this has a direct influence on which non-linear definitions of stresses and strains are used. (a)

Small strain GNL problems

These problems are associated with small or large rotations. Problems in which small rotations occur include shallow struts, shells and arches deflected by a transverse load, clamped circular plates under transverse point loads and shallow spherical caps. Examples of large rotation problems include a fishing rod bent under the weight of a heavy fish, buckling of an imperfect Euler strut, and a deep arch. (b)

Large strain GNL problems

These are the most complex of all GNL problems and are usually associated with metal forming and manufacturing processes, such as deep drawing of drink cans, forging, extrusion and rolling. With large strains, it is also important to model material non-linearity such as plasticity. An exception is rubber which can undergo very large strains, of the order of unity, but remains elastic. This type of behaviour is called ‘hyperelastic’ or ‘non-linear elastic’ behaviour. The constitutive equations for rubber can be derived from the expressions for the potential energy density. Definitions of stresses and strains in GNL problems

The conventional definition of ‘engineering strain’ may not be adequate when dealing with GNL problems because it measures the change in length over the original (undeformed) length. A more suitable definition is one which takes into account the new length, such as the so-called ‘logarithmic strain’. Two strain definitions which 12.20


have been widely used in GNL problems are called the ‘Green strain’ and ‘Almansi strain’. These strains are based on the ‘square’ of the length, and have been shown to be very effective in dealing with a wide rang e of GNL problems, including large strain problems As in the strain measures, the ‘engineering’ or ‘nominal’ stress, defined as the force divided by the original undeformed area, may be inappropriate for use in GNL problems in which the cross-sectional area may exhibit large changes. Instead, as in material non-linearity, a ‘truestress’ (also called ‘Cauchy stress’) can be defined as the force divided by the ‘current’ cross-section area, rather than the original area. Another feature of GNL problems is that the relationships between stresses and strains have to be carefully defined. In conventional elasticity equations, stresses are linked to strains through the constitutive law, i.e. Hooke’s law. In GNL problems, in addition to the constitutive equations, stresses are usually associated with the corresponding strains using the ‘virtual work theorem’ or ‘total potential energy theorem’. This is important in large strain problems, where such stresses are called ‘work-conjugates’ to the corresponding strains. Conservative and non-conservative (follower) loads

A conservative load is that which always applies in a fixed direction regardless of the deformation of the body. A typical example is a gravitational load, which aways applies vertically. A non-conservative (follower) load is one which changes its direction during the deformation, i.e. it follows the deformation of the body, e.g. an internal pressure in a vessel changes its position and direction as the vessel deforms, in order to remain perpendicular to the surface. Figure 12.11 shows a schematic representation of conservative and non-conservative loads.

Fig. 12.11. Conservative and non-conservative loads.

12.21


Formulations

Formulation of geometrical non-linear problems requires choice of reference systems for describing the structures geometry and deformations. The most common modes of describing the deformations of solids and fluids are the Eulerian and the Lagrangian approaches:  Eulerian description of motion is also denoted spatial description , because it refers to what happens at a certain place in space. In the description the current coordinates x1, x2, x3 and time t are the independent variables. This description is especially suited for hydrodynamics. 

Lagrangian description of motion refers to what happens at a material particle. Hence, this description is also denoted material description . Independent variables are the initial coordinates x1, x2, x3 and time t.

The Lagrangian description of motion is commonly preferred in structural mechanics since the initial configuration is usually known. This formulation will be adopted in the present work. Adopting a Lagrangian description, the geometry at a certain load condition is to refer it to  the initial geometry/fixed global coordinate system  assume that the loading takes place in a stepwise manner and refer each element to a local coordinate system that follow the structure during deformation. This is called a corotational system . The coordinates are updated . The purpose is to show how incremental relations like

ΔR =

K I Δr

can be established for various types of structures. However, since the assemble of the global stiffness is straightforwardwhen the element relationship S = kv is known, the focus will be on the element relationship. It is interesting to note that the total Lagrange formulation for the two-bar problem corresponds to Eq. 12.13. K ITL

= K 0 + K G + Kσ

(12.31c)

The updated Lagrange formulation corresponds to K I UL

= K 0 + K σ

(12.31d)

where S (positive in tension) is the actual S at any load step. A simplified version of the latter reference system is described in Section 12.3.2. The Section 12.3.3 the use of the first type of reference system will be demonstrated. 12.22


12.3.2 Method with updated coordinates General

In this procedure the element stiffness relationships are first determined in a local system and are then transformed into a fixed, global coordinate system before the global stiffness relation is assembled. The location of the local coordinate systems would have to be updated when deformations cause the geometry to change. In this method the non-linear geometrical effects are accounted for by continuously changing the transformation matrices. It is assumed here that each element behave in a linear manner when it is referred to a corotational coordinates. This implies that small deformations on the local level are assumed. Element relationships

Fig. 12.12 shows an arbitrary plane element in its initial and deformed condition. The element is connected to global nodes a and b and a local axis defined by a straight line through these two points. The nodal forces (S) (in the local system) are defined in the conventional manner. The deformations of the element are completely described by the axial elongation, u and rotation ωa and ω b, which are referred to the chord ab (Fig. 12.12b). By means of these three deformation parameters the force vector S can be uniquely obtained from EA ⎤ ⎡ 0 − 0 ⎢ ⎥  ⎢ 6 EI 6 EI ⎥ ⎡ S 1 ⎤ ⎢− 2 − 2 0 ⎥ ⎢ S ⎥ ⎢   ⎥ ⎢ 2 ⎥ ⎢ 4 EI 2 EI ⎥ 0 ⎢ S 3 ⎥ ⎢  ⎥  S=⎢ ⎥ = ⎢ EA ⎥ ⎢ S 4 ⎥ ⎢ 0 0 ⎥  ⎥ ⎢ S 5 ⎥ ⎢ ⎢ ⎥ ⎢ 6 EI 6 EI 0 ⎥ ⎢⎣ S 6 ⎥⎦ ⎢  2 2 ⎥  ⎢ 2 EI 4 EI ⎥ 0 ⎢ ⎥  ⎣  ⎦

a) Nodal forces

⎡ω a ⎤ ~~ ⎢ω ⎥ = k b ⎢ ⎥ v ⎢⎣ u ⎥⎦

(12.32)

b) Deformations

Figure 12.12 Forces and deformations in a local coordinate system. 12.23


The vector ~ v describes the deformations relative to the local coordinate system xy. In addition the element will be subject to rigid body motions (translation and rotation) that correspond to the location of the local xy-system relative to the global x y -system, and are defined by the global coordinates of nodes a and b. Rigid body motions do not imply internal forces in the element. Hence, the relationship (12.32) does not depend on rigid body motions and is valid for any position of the local system. The reason ~ why k has dimensions 6 by 3 is that the three components of the rigid body motion is removed from Eq. (12.32). The assumptions of small displacements implies that u, ωa and ω b are small. The next step is to express the global nodal forces, S when the global displacement v are known. S and v are defined in Fig. 12.13a and b.

a) Nodal forces Figure 12.13 positive.

b) Displacements

Global nodal forces and displacements. All angles in for figure are

The relationship between the nodal forces defined in the global and local coordinate system is given by S = T T S

(12.33)

where

⎡ cosθ ⎢− sin θ ⎢ ⎢ 0 T=⎢ ⎢ ⎢ ⎢ ⎢⎣

sin θ

0

cosθ 0 0

0

1

0

cosθ

sin θ

− sin θ

cosθ

0

0

⎤ ⎥ ⎥ ⎥ ⎥ 0⎥ 0⎥ ⎥ 1⎥⎦

and θ is the angle between the x - and x-axes (Fig. 12.12). By combining Eqs (12.32, 12.33) gives ~ S = T T S = T T k ~ v = k * ~ v

(12.34) 12.24


~ k * = T T k can be expressed explicitly by carrying out the matrix multiplication by hand. If the local deformations ~ v are expressed by the global displacements v Eq. (12.34) can be transformed into a relationship between S and v . This can be achieved as follows. The coordinates of nodal points a and b in deformed condition, are:

= y a = xb = y b =

x a

x a 0 +

v1

y a 0 +

v2

+ y b 0 +

v4

xb 0

(12.35)

v5

where x a 0 etc. are the coordinates of the initial condition. The chord length ab in deformed condition is obtained as 

= ( xb − xa )2 + ( yb − y a )2

(12.36)

Relative to the initial condition the rotation of the chord axis ab (Fig. 10.15) is ψ = θ − θ 0

⎛ yb 0 − y a 0 ⎞ ⎟⎟  0 ⎝ ⎠ ⎛ y − y a ⎞ θ = arcsin ⎜ b ⎟ ⎝  ⎠ θ 0 = arcsin ⎜⎜

where

 0 is

(12.37)

the element length in the initial condition. To determined θ0 and θ in a

unique manner, the sign of cosθ0 and cosθ needs to be checked. The local displacements are then obtained (from Fig. 12.12b, 12.13b) as:

+ ψ ⎤ ⎥ ⎥ ⎢ + ψ ⎥ ⎥⎦ ⎢⎣  −  0 ⎥⎦

⎡ω a ⎤ ⎡v3 ⎢ ⎥ ⎢ ~ v = ω b = v6 ⎢ ⎢⎣ u

(12.38)

Eqs (12.34, 12.38) make it possible to find the element forces S for any known displacements, v . This relation for an individual element can then be used to express the relationship for the structural systems in the following manner. Structural system relationships

The relationship between the displacement vector for element No. i and the global displacement vector r is expressed by the kinematic relation: 12.25


vi

= air

(12.39)

Suppose that the structure is subjected to given external forces, corresponding to a vector R . The principle of virtual displacements yields

∑ (~v ) S i T

= ~r T R

i

i

where the (~) indicates virtual displacements. compatible ~i v

Since virtual displacements are

= a i ~r

which yields ~ r T

∑ (a ) S i T

i

= ~r T R

i

Since ~ r can be arbitrary, the equilibrium condition becomes

∑ (a )

i T

= R

Si

(12.40)

i

The left-hand side of Eq. (12.40) is a function of the displacement vector, through Eqs (12.34, 12.38, 12.39). Eq. (12.40) is a non-linear equilibrium equation. The left-hand side expresses the sum of internal forces from each element in the structure. To find the displacements r that fulfill Eq. (12.40) requires iteration. If r’ is a displacement vector that does not fulfill Eq. (12.40), the corresponding element forces S i do not satisfy this equation. Hence, this displacement condition implies a set of unbalanced forces or residual forces , given by: R r =

∑ (a )

i T

S ′i

− R

(12.41)

R r represents the additional forces required to fulfill the global equilibrium in this case (with r’). R + R r =

∑ (a ) i

T

S ′i

(12.42)

i

As indicated in Section 12.1.2 and discussed more in detail later, an incremental stiffness relation Eq. (2.103), will be applied. In the following it will be shown how this relationship can be established. In a local coordinate system the incremental stiffness relation for an element can be obtained as shown in Section 2.9.4 and may be written as: (when a tensile axial force is considered positive) (12.43) dS = ( k 0 + k σ ) dv = k I d v 12.26


where S and v is defined as in Fig. 12.14 and the S vector correspondingly.

Figure 12.14 Local degrees of freedom for an element. k 0 and k σ are the elastic and geometrical stiffness matrices, respectively as defined by Eqs (2.110) and (2.111). However, now the axial force is considered positive if in tension (this explains the +sign before k σ in Eq, (12.43)). This force can be calculated from the axial strain by: u (12.44) P = EA 0

By transformation to a global system d S = T T d S = T T k I d v = T T k I d v = k I d v

(12.45)

where k I = T T k I T S and v are defined in Fig. 12.12 and T in Eq. (12.33). Kinematic relationship and static equilibrium for a system of elements can be expressed as: d v i

= a i d r,

d R =

∑ (a )d S i

i

i

The incremental system stiffness then becomes d R =

∑ (a ) d S = ∑ (a ) i T

i T

i

i

k I i a i d r

i

or d R = K I d r

(12.46)

where K I =

∑ (a )

i T

k I i a i

i

Eqs (12.40, 12.46) form the basis for obtaining solutions of the geometrically nonlinear problem, by the methods described later in Section 12.4. The method described in this section is relatively simple in use since most of the matrices are well known from linear analysis. 12.27


Two practical issues should be noted when implementing this method in a computer program. The stiffness matrices k 0 and k σ could be calculated based on initial element length  0 , rather than the updated length,  , without loss of accuracy. In this way the element stiffness matrices can be calculated once and stored and retrieved when needed. However, the axial force, P in k σ need to be updated. The second issue is that the coordinate updating, Eq. (12.35), element length, Eq. (12.36) and rotations, Eq. (12.36) need to be done accurately to avoid that round off errors reduce/destroy the accuracy and cause convergence problems in the iterative solution. This formulation presented in this section is often called “Updated Langrange formulation”. However, this is not quite true, since corotated coordinate system in general will be a curvilinear one at not a cartesian one as in this formulation. The formulation is, therefore, an approximately updated Lagrangian formulation. Alternatively it may be called a method based on corotated coordinates or corotating coordinates.

12.3.3

Total Lagrangian formulation for a beam with axial and lateral deformations In the previous section a local coordinate system was assumed to follow the rigid body displacements as the structure deforms. The deformations are referred to these updated coordinate systems. Contrary to this, the total Lagrange formulation is based on a fixed coordinate system. When calculating displacements and strains, the expressions used in previous sections for the strain then need to be modified.

In this section a beam undergoing axial and lateral deformations will be used to illustrate the total Lagrangian formulation. Classical linear beam theory will serve as basis for this. This is because shear deformation contribute little to the behaviour of beams that have large deformatio ns. Fig. 12.15 shows a beam element undergoing large displacements. The displacements of a plane beam are uniquely described by the displacement u and w of the neutral axis. For a case with small displacements the strain at a location (x, z) is ε x ( x, z ) = u , x − z ⋅ w, xx ( x)

(12.47)

If the displacements and rotations are large (but the strains are small), the strain due to displacement in the x-direction remains unchanged. However, the large deflection, w(x) causes an additional axial strain. This effect can be estimated as follows:

a) Element with large displacement

b) Geometric locations for a differential element of length dx

Figure 12.15 Kinematics for beam with large displacements. 12.28


Let a small lateral displacement w = w(x) take place. Thus each differential length dx is changed to a new length, ds, where ds > dx because the distance between supports is not allowed to change. From Fig. 12.14b.

ds

= (1 + w,

)

2 1/ 2 x

⎛ w, x2 ⎞ ⎟⎟ dx dx ≈ ⎜⎜1 + 2 ⎝ ⎠

(12.48)

where the latter approximation comes from the first two terms of the binomial expansion. The approximation is valid is w, x2 << 1 , which restricts this development to small rotations. Axial membrane strain, therefore

=

ε m

ds − dx dx

≈

w, x2

2

εm in

the bar due the large deflection, w is

(12.49)

In the linear theory of elasticity we ignore terms of order w, x2 . But here we seek the consequences of retaining the more important of the higher-order terms that linear theory neglects. We are taking a physical approach to formulating these terms. They may also be obtained by a systematic procedure of linearization, as will be touched upon in Appendix C. This means that a beam with axial displacement of the neutral axis u = u(x) and lateral displacement w = w(x) has a strain of ε x

= u, x +

1 2

w, x2 − z ⋅ w, xx

(12.50)

It is noted that if the axial strains are not small, an additional terms needs to be included in Eq. (12.50). The strain then becomes E xx

= u, x − z ⋅ w, xx +

1 2

u , x2 +

1 2

w, x2

(12.51)

This is the so-called Green’s strains as presented in Appendix B. For metal structures, however, the additional term 12 u, x2 is negligible. Another issue is the effect of an initial lateral deflection, w ( x) of the beam. A beam with initial lateral deflection may also be considered a shallow arch. In this case the length of an arch with a projected length along the x-axis equal to dx would be ds = 1 + 12 w , x2 dx , according to Eq. (12.40). When an additional lateral deformation,

)

w(x) occurs, the resulting arch length becomes

[

ds* = 1 + (w , x + w, x )

]

2 1/ 2

⎡ ⎣

dx ≈ ⎢1 +

1 2

(w , x + w, x )2 ⎤⎥ dx ⎦

(12.52)

and the membrane strain due to lateral direction becomes 12.29


=

ε m

ds * − ds ds

= w , x w, x +

1 2

w, x2

(12.53)

For metal beams with initial lateral deflection, w ( x ) the Green strain becomes E xx

= u, x − z ⋅ w, xx + w , x w, x +

1 2

w, x2

(12.54)

Alternatively Exx may be written as

E xx

= u, x − zw, xx + (w ,x + w,x )2 − w ,x2

Stresses

Stresses may be defined with reference to the deformed structure or its initial configuration. True stresses referred to the deformed configuration are denoted σij (Eulerian or Cauchy stresses). Stresses referred to the initial configuration are denoted Sij (2 nd Piola-Kirchhoff stresses). The latter stress is consistent with the Green’s strain that also refers to the initial configuration. If the Green’s strain is applied the stresses should therefore be Piola-Kirchoff stresses. For the one-dimensional case the Piola-Kirchoff stress is related to the Eulerian stress (defined in the deformed configuration) as: S xx

=

∂ X ⎛ ∂u ⎞ σ σ xx = ⎜1 − ⎟ xx ∂ x ⎝ ∂ x ⎠

where X and x are coordinates of the initial and deformed configuration, respectively. In metal structures the strains are, but the rotation may be large and nonlinear geometric effects. In such cases S xx and σxx are for all practical purposes, considered to be equal. Virtual work (Principle of virtual displacements)

Assume that a beam structure in an equilibrium position is given a virtual axial and lateral displacement, δu and δw, which are kinematically consistent with the boundary conditions. By neglecting volume forces the equation of virtual work may according to Appendix B.2 the equation be expressed a

∫ S δ E = ∫ q xx

V

xx



u

∫

δ u ( x)dx + q w δ w( x ) dx

(12.55)



It is noted that Exx and Sxx are expressed by the displacements u and w according to Eq. (12.54) and Hooke’s law for an elastic beam. By applying the approach in Appendix B the incremental form of the virtual work can be written as 12.30


∫ ΔS

xx

∫

δ E xx dV + S xx Δ δ E xx dV +

V

V

∫ Δ S

xx

Δ δ E xx dV

V

⎡ ⎤ = ∫ Δ qu δ u dx + ∫ Δ q w δ w dx − ⎢ ∫ S xx δ E xx dV − ∫ qu δ u dx − ∫ q w δ w dx ⎥ ⎣⎢V ⎦⎥

(12.56)

If the configuration Cn is in equilibrium the parenthesis on the right hand side of Eq. (12.56) will vanish according to Eq. (12.55). However, due to approximations in the solution procedure equilibrium will not be generally satisfied. Hence, the terms in the parenthesis on the right h and side will serve as equilibriu m correction terms. Finite element model

The incremental stiffness relation for an element can then be established based on Eq. (12.56) and choice of interpolation functions for the displacements u =[u, w]T. This implies choice of nodes and degrees of freedom in each node. One option is shown in Fig. 12.16, based on: - quadratic polynomial for u, with the parameter vu = (u1, u2, u3)T - cubic polynomial for w, with the parameters vw = (w1, θ1, w2, θ2)T The initial lateral deflection, w may be described by the same type of polynomial as w. In matrix notation the displacement may be written as u ( x) = N u v u w( x) = B w v w w ( x)

(12.57a-c)

= Nwvw

z

Figure 12.16 Stiffener elements. If the matrix notation of Eqs (12.57a-c) are introduced, Exx in Eq. (12.54) may be written as E xx

= N u , x v u − z ⋅ N w , xx v w + v T w N T w , x N w , x v w +

1 2

v T w N T w , x N w , x v w

(12.58a)

or 12.31


E xx

1

= Bv + v T w N T w , x N w , x v w +

2

v T w N T w , x N w , x v w

(12.58b)

where B = [N u , x , N w , xx ] v

⎡v ⎤ =⎢ u⎥ ⎣v w ⎦

(12.58c)

The first variation in strain, δExx due to a variation in displacement δv is δ E xx (Note: δ E xx

= B δ v + v T w N T w' x N w , x δ v w + = E xx (v + δ v ) − E xx ( v) ,

1 2

v T w N T w , x N w , x δ v w

(12.59)

which is most easily performed on the scalar

expression (12.59) first. The high order term (δw)2 is neglected.) From one configuration to another, which is close to the first one, the change in strain variation is

Δδ E xx = Δv T w N T w , x N w , x δ v w

(12.60)

The stress increment between the two configurations becomes

ΔS xx = DΔ E xx

(12.61)

where D is the incremental material stiffness coefficient, elastic or elasto-plastic. See Section 12.4. By recognizing the integrand of the first integrand of Eq. (12.56) is a scalar, it may be written as:

∫ ΔS

xx

∫

T δ E xx dV = δ E xx D Δ E xx dV

V

(12.62)

V

When introducing the strain expression (12.58 b-c) into Eq. (12.62), the expression is simplified if the following notation is introduced: v wt

= vw + vw

(12.63)

Then, by using Eq. (12.62)

12.32


∫ ΔS

xx

∫(

)

δ E xx dV = δ v T B T D B dV Δv

V

V

+ δ v

T

∫ (B D v T

T wt

)

N T w , x N w , x dV Δv w

V

+ δ v T w ∫ ( N T w , x N w , x v wt D B dV )Δv

(12.64)

V

+ δ v T w ∫ ( N T w , x N w , x v wt D v T wt N T w , x N w , x dV )Δv w = δ v T k 00 Δv + δ v T k 0 L Δv w + δ v T w k T 0 L Δ v + δ v w k LL Δv w By introducing

vw

⎡v ⎤ = [0 I ]⎢ u ⎥ = Hv ⎣v w ⎦

(12.65)

Eq. (12.64) may be written as

∫ ΔS δ E xx

xx

(

dV = δ v T k 00

+ k 0 L H + H T k T 0 L + H T k LL H )Δv

V

= δ v

(12.66)

k 1Δ v

T

The second term on the left hand side of Eq. (12.56) becomes

∫ Δ δ E S xx

V

xx

dV = δ v T w

∫ (N

)

, N w , x S xx dV Δv w

T w x

V

= δ v k Lσ Δv w = δ v T H T k Lσ HΔv = δ v T k 2 Δv T w

(12.67)

The third term on the left-hand side of Eq. (12.56) is neglected as a higher order term. The load term of Eq. (12.56) may be formulated in matrix notation in a similar manner as for linear analysis. The resulting virtual work during an increment may be written as T

∫ Δqu δ u dx + ∫ Δq w δ w dx = δ v ΔS

(12.68)

When neglecting the third term on the left hand side and the parenthesis on the right hand side, Eq. (12.56) hence becomes δ v T [(k 1

+ k 2 ) Δv − ΔS ]= 0

(12.69)

Before completing the derivation the nodal set of parameters, v may be rearranged from v

= [u1 , u 2 , u 3 , w1 ,θ 1 , w2 ,θ 2 ]T

(12.70a) 12.33


to v ' = [(u , w,θ )1 ; (u , w,θ ) 2 ; u 3 ] T

(12.70b)

The connection between v and v’ is given by the matrix G: v = Gv'

(12.71)

Such a rearrangement involves a matrix with elements which are 0 or 1. Introducing

= GT k1G; k '2 = GT k 2 G k I ' = k '1 + k '2 ΔS = GT ΔS k1'

(12.72)

Eq. (12.69) turns into δ v 'T ( k I ' Δ v '− ΔS ' ) =0

(12.73)

Since the virtual displacement δv’ is arbitrary, the principle of virtual work gives k I ' Δv ' − ΔS ' = 0

(12.74)

This is the stiffness relation on incremental form for the beam element. The (9x9) symmetric matrix is k I ' is called incremental element stiffness matrix. k '1 is denoted initial displacement matrix (which includes the linear small displacement stiffness matrix) or large displacement matrix and the other contribution k ' 2 is known as initial stress matrix or geometric matrix.

In a similar manner the incremental stiffness relation for a beam element based on the updated Lagrangian formulation can be obtained. In this case the resulting incremental element stiffness matrix may be written as: k ' = k '0 + k 'σ

(12.75)

where the stiffness matrix k '0 is that due to small deflection theory, namely the first term on the left hand side in Eq. (12.75), and k σ ' is the (linearized) geometric stiffness matrix, resulting from second integral on the left hand side. In linear analysis of plane structures the behaviour in plane (axial) behaviour and out of plane or lateral (beam bending) is uncoupled. In geometrically nonlinear behaviour, axial and lateral behaviour is coupled. This can be seen for a Total Lagrangian Formulation from the fact that mean strain depends upon u,x as well as 0.5 w, x2 .

12.34


Comments

Various finite element models may be established based on the choice of interpolation polynomials for u and w. Since the axial strain would vary significantly over the element undergoing large deflections, the displacement u should be approximated by a quadratic or higher degree polynomial, while w would be interpolated by a cubic polynomial. Moreover, this choice of displacements reduces the effect of “selfstraining”, a phenomenon resembling “shear looking” for membrane elements. When the interpolation polynomials are given, the stiffness matrix is calculated by means of the expressions (12.73, 12.72, 12.67, 12.64-12.66). Internal degrees of freedom, e.g. associated with one of the d.o.f. associated with the displacement u are normally eliminated by static condensation, see Chapter 9. Finally, the global incremental stiffness relation is obtained by the direct stiffness method.

12.3.4 Generalization

In the previous section nonlinear geometrical formulations for one dimensional bars/beams are treated. The one dimensional formulation presented above can be generalized to two dimensional problems, i.e. thin plates subjected to in-plane and lateral loads. This generalization involves 



linearized problem in terms of a stiffness and geometric stiffness matrix that can be used for buckling analysis of plates subjected to in-plane forces per unit length: Nx = h⋅σx, Ny = h⋅σy and Nxy = h⋅τxy, where h is the plate thickness. Updated or Total Lagrangian Formulation for geometrically nonlinear problems.

Linear buckling analysis

Analogous with the expression (2.106) for a beam the virtual work done by in-plane forces Nx. Ny and Nxy acting on a plate is W ext

~, w, + N w ~ ~ (~ )) = ∫ ( N x w x x y , y w, y + N xy w, x w, y + w, y w, x dA ~, w ~, ⎡ N x N xy ⎤ ⎡ w, x ⎤ dA w ⎥⎢ ⎥ x y ⎢ ⎣ N xy N y ⎦ ⎣ w, y ⎦

= ∫[

]

(12.76)

Tensile forces are considered to be positive. By introducing the interpolation polynomial w = Nv

12.35


⎡ w, x ⎤ ⎡ N, x ⎤ ⎢ w, ⎥ = ⎢ N, ⎥ v = Bσ v ⎣ y ⎦ ⎣ y ⎦

(12.77)

the external work done is written as Wext = v T

⎡ ∫ BσT G Bσ dA⎤ v ⎣ ⎦

(12.78)

where

⎡ N = ⎢ x ⎣ N xy

G

N xy ⎤ N y

⎥ ⎦

(12.78a)

By equating internal and external work it is found that the element stiffness matrix k

= k '0 + k σ '

(12.79)

where k 0’ is the stiffness matrix for amall displacement theory as given in Chapter 7. k σ'

= ∫ BT σ G Bσ dA

(12.79 a-b)

Geometrical nonlinear effects can be formulated for two-dimensional structures, such as plates and shells according to the same principles as for bars and beams. Reference is made to textbooks such as Crisfield (1991), Hughes (2000), NAFEMS (1992), Taylor and Zienkiewicz (2000).

12.4 Nonlinear material behaviour 12.4.1 One dimensional case

In Section 12.1.3 typical one-dimensional stress-strain relationships for metals were shown, and it was pointed out that a nonlinear relationship occurs when the stress exceeds the proportionality limit. A material is called nonlinear if stresses and strains are related by a straindependent matrix rather than a matrix of constants. Thus the computational difficulty is that equilibrium equations must be written using material properties that depend on strains, but strains are not known in advance. Plastic flow is often a cause of material nonlinearity. The present section deals with formulation of elastic-plastic problems by considering the one dimensional case. Assume that yielding has already occurred; then a strain increment dε (between points A and B) takes place (Fig. 12.17a). This strain increment can be regarded as 12.36


composed of an elastic contribution dεe and a plastic contribution dε p, so that dε = dεe + dε p. The corresponding stress increment dσ can be written in various ways d σ = Ed ε e

= E (d ε − d ε p ); d σ = E t d ε ;

and d σ = H ' d ε p

(12.80)

where H’ is called the plastic tangent modulus as given by ∂σ / ∂ε p . (Note that H(ε p ).) Substitution of the first and third of Eqs. (12.80) into the second yields H ' =

1 1 E t

−

1 E

⎛ ⎝

or E t = E ⎜1 −

⎞ ⎟ E + H ' ⎠ E

σ=

(12.81)

where Et is the tangent modulus. When written in this form, the expression for Et is similar to a more general expression used for multiaxial states of stress. If E is finite and Et = 0, then H’ = 0, and the material is called “elastic-perfectly plastic”.

a)

Stress-strain plot in uniaxial stress, idealized as two straight lines, where σy is the stress at first onset of yielding.

b) Kinematic and isotropic hardening rules.

Figure 12.17 characteristic features of one-dimensional stress-strain relationships. A summary of elastic-plastic action in uniaxial stress is as follows: 1) The yield criterion states that yielding begins when |σ| reaches σY, where in practice σY is usually taken as the tensile yield strength. Subsequent plastic derformation may alter the stress needed to produce renewed or continued yielding; this stress exceeds the initial yield strength σY if Et > 0. 2) A hardening rule , which describes how the yield criterion is changes by the history of plastic flow. For example, imagine that the material first has been loaded to point B and then unloading occurs from point B to point C in Fig. 12.17a. With reloading from point C, response will be elastic until σ > σB, when renewed yielding occurs. Assume then that unloading occurs from point B and progresses into a reversed loading as shown in Fig. 12.17b. If the yielding is assumed to occur at σ | | =σB the hardening is said to be isotropic. However, for common metals, such a rule is in conflict with the observed behaviour that yielding reappears at a stress of approximate magnitude σB - 2σY when loading is reversed. Accordingly, a better match to observed behaviour is provided by the “kinematic hardening” rule, which (for uniaxial stress) says that a total elastic range of 2 σY is preserved. 12.37


3) A flow rule can be written in multidimensional problems. It leads to a relation between stress increments d and strain increments d . In uniaxial stress this relation is simply dσ = Et dε, which describes the increment of stress produced by an increment of strain. Note, however, that if the material has yet to yield or is unloading, then dσ = E d ε (e.g., in Fig. 12.17a, complete unloading from point B leads to point C and a permanent strain ε p). The discussion in the foregoing paragraph does not require that post-elastic response be idealized as a straight line. In other words, Et, need not be constant. Bar structures with vriable cross-section area

As a simple application of one-dimensional plasticity, imagine that a tapered bar is to be loaded by an axial force P (Fig. 12.18). Material properties are those depicted in Fig. 12.17. The bar is modeled by two-d.o.f. bar elements, each of constant cross section. For elastic conditions, the element stiffness matrix is given by Eq.(2.48), where E = dσ/dε when |σ| < σY. Upon yielding, the stress-strain relation becomes Et = dσ/dε. Accordingly, letting Eep represent the “elastic-plastic” stiffness, the element tangent-striffness matrix is expressed by: k T

=

AEep L

⎡ 1 −1⎤ ⎢ ⎥ ⎣ −1 1 ⎦

(12.82)

where Eep = E if the yield criterion is not exceeded or if unloading is taking place, and Eep = Et if plastic flow is involved. In numerical solutions, material may take the transition from elastic to plastic within an iterative cycle of the solution process. For example, imagine that dε spans εD to εA in Fig. 12.17a. The problem of “rounding the corner” can be addressed by combining E and Et according to the fraction m of the total step dε that is elastic. Thus let E ep

= mE + (1 − m) E t

a) Geometry and model

where m =

ε Y − ε D ε A

− ε D

(12.83)

b) Load-displacement incrementation

Figure 12.18 Tapered bar subjected to axial loading

12.38


Alternatively, by substituting stresses for strains and using the fictitious stress σ* =EεA, we can write m in terms of stresses, as m = ( σY - σD)/(σ* - σD). Refinements of this scheme are possible. As another simple elasto-plastic problem consider the single degree of freedom elastoplastic problem shown in Figure 12.19 . Each bar is constructed of an identical material which has bilinear elastic perfectly plastic behavior. Fig. 12.17a with Et = 0. Each bar is of cross-sectional area A and elastic modulus E . The yield stress values for the material in the bars are σ Y(1) , σ Y(2) and σ Y(3) where σ Y(1) < σ Y(2) < σ Y(3) . When the force f is applied all behaviour is initially elastic. As the load is increased the stress in element 1, reaches the yield and can carry no further increase of stress. Hence, the internal resisting force emanating from element 1 from this stage onward is p (1)

= σ Y (1) A

(12.84)

and element 1 loses its stiffness for an increasing load. As the load is increased the additional load must be balanced by additional internal resisting forces in elements 2 and 3. Eventually, the stress in element 2 reaches the yield stress and this element also yields, loses its stiffness and can carry no further increase of load. Thus, all extra load must be taken by element 3. Finally, element 3 yields and loses its stiffness and the structure can take no further load. The load at which this element yields is f ult = A( σ Y(1)

+ σ Y(2) + σ Y (3) )

(12.85)

Thus, the source of nonlinearity in the elastoplastic problem occurs in the evaluation of the stress in the elements and hence in the internal resisting force. The stress is effectively a nonlinear function of the displacements. The load-displacement curve for this problem is shown in Figure 12.19 and consists of a piecewise linear curve with the end of each segment signaling the onset of yield of a new element until failure of the whole structure. The slope of the load-displacement curve at any stage og the analysis is called the ‘tangential stiffness’ K T. When elastic KT = 3EA / 

(12.86)

As the first bar yields, KT = 2 EA / 

(12.87)

When the second bar yields KT = EA / 

(12.88) 12.39


and finally, after the last bar yields K T = 0

a) Three-bar

b) Equilibrium

c)Load-displacement curve

Fig. 12.19 Three-bar elastoplastic problem Solution procedures for elasto-plastic problems will be discussed in Section 12.5.

Beam structures

The deformation of a beam subjected to axial force and bending is described by assuming that plane sections remain plane. This implies that the strain is given by Eq. (12.54). The corresponding stresses need to be calculated in each layer (co-ordinate z referred to the centroid (neutral axis)) at a longitudinal location, x, by using the incremental expression: d σ E t d ε . In Section 12.4.2 a more accurate expression for beams with thin-walled cross-sections will be given for the incremental stress-strain relationship. Even fo Elasto-plastic material behaviour plane section remain plane. Hence the strain generally can be written as ε ε m

z κ

(12.89)

where εm and κ is the membrane strain and curvature, respectively. In Section 12.4.2 εm and κ have been described by their interpolation functions in terms of the coordinate x. The volume integrals that express the stiffness matrices, Eqs. (12.64, 12.67), can be carried out as follows ~ Dε dv ε v

~ [ε m

~ ) D (ε z κ m

zκ ) dA]dx

(12.90)

x A

For an elastic beam with D=E and z is defined with reference to the centroid the integral (12.90) may be written as ~ Dε dV ε v

~ ε dx EAε m m

~κ dA) dx ( Ez 2κ

x

~ ε dx EAε m x

x A

~κ dx EI κ

(12.91 )

x

12.40


~ κ dA 0 , when z is referred to the centroid. This is because for instance zε m A

For a beam with Elasto-plastic material behaviour the bending stress depends upon the strain e.g. as indicated in Fig. 12.20. Eq. (12.90) then needs to be integrated numerically and would include coupling terms of for instance the form: ~ κ dV z D( z )ε m v

σ

∫ σ dA = 0 ∫ σ ⋅ zdA = M

N.A. H ε

elastic a) stress-strain curve

(no axial force) elasto-plastic

b) stress distribution in beam secti on for various level of bending moment (Note: the bending strain is ε b = − z ⋅ H , i.e. varying linearly over the cross-section)

Formatted: Font:

10 pt

Formatted: English

(U.S.)

Formatted: English

(U.S.)

Formatted: English

(U.S.)

Formatted: Font:

10 pt

Formatted: English

(U.S.)

Formatted: English

(U.S.)

Formatted: English

(U.S.)

Formatted: English

(U.S.)

Fig. 12.20 Stress distribution in a beam cross-section as a function of bending moment (no axial force) It should be noted that the integral (12.90) also would include coupling terms of z is not defined with refrence to the centroid. The integration over the cross-section must be carried out with check of the stress level to decide whether: D = E (elastic) D = Et (elasto-plastic), Eqs. (12.80, 12.81) and needs to be carried by numerical integration by either standard Gaussian or by Lobatto integration (Appendix, Crisfield (1991)). In the formulation described above the stress and strain is described over the height of the beam- as a formulartion of z. Moreover, the relevant material properties in different sections along each beam element need pursued. However, if the behaviour of beams and frames under lateral (and axial) loading is observed, it is often seen that failure occurs by the formation of plastic hinges . Simulated by this observation more efficient methods can be developed, based on the assumption that plastic hinges develop at predefined locations. Moreover, the elasto-plastic behaviour in a hing is handled by using axial force, N and bending moment, M instead of pursuing the stresses (and strains) in each fiber level of the corss-section. Future comments on the use of the plastic hinges and generalized forces (M,N) are provided in Section 12.6.2.

12.4.2 Generalization 12.41


In the same manner as for one-dimensional case elasto-plastic behaviour of metals in multidimensional stress state is characterized by   

An initial yield condition , i.e. the state of stress for which plastic deformation first occurs. A hardening rule which describes the modification of the yield condition due to strain hardening during plastic flow. A flow rule which allows the determination of plastic strain increments at each point in the load history.

It is assumed that the material is isotropic, which implied that the stiffness properties are independent of orientation at a point. A review of elasto-plastic theory for multi-dimensional stress states based on isotropic hardening is given in Appendix C. It is shown that the relationship between stress and strain increments may be written as follows ep = Dijkl d ε kl

(12.92a)

ep Dijkl

= E ijkl − β sij s kl

(12.92b)

E ijkl

= 2Gε ije +

d σ ij

where

β =

σ

=

2ν G 1 − 2ν

9G 2 σ 2 ( H '+3G ) 3 2

sij sij

=

e δ ij ε kk

(12.92c)

3 2

(12.92d)

σ ijσ ij

− 12 (σ kk ) 2

(12.92e)

where the Einstein’s summation convention is applied i.e. Dep ijkl dε kl

n

n

= ∑∑ Depijkldε kl ;

n

s ij s ij

k =1 l =1

n

= ∑∑ sij2 ; σ

kk =

i =1 j =1

n

∑σ k =1

kk

;

σ m = 13 σkk .

(12.92f)

sij denotes the deviatoric, or reduced stress tensor sij

where

= σ ij − 13 δ ijσ kk

(12.92g)

δij is equal to 1 if i=j and equal to 0 if i± j.

It is seen that Eq. (12.92a ) may be simplified to the following form if one dimensional plane stress condition is considered with dεxx ≠ 0 and dεyy = dεzz = 0 and ν= 0.3

12.42


⎛ E(1 − ν ) 9G 2 ⎞ − dσ xx = ⎜ ⎟ dεxx + ν − ν + (1 )(1 2 ) (H ' 3G) ⎝ ⎠ 2 ⎛ 1.33E ⎞ ≈ ⎜1.33 − ⎟ dε xx H '+ 1.16E ⎠ ⎝

(12.93)

It is seen that Eq. (12.93) resembles the second of the expressions (12.81). In dealing with thin-walled metal structures, a plane stress condition is more relevant. The incremental stress-strain relationships for one- and two-dimensional conditions can be obtained from Eq. (12.92a-d). The resulting expressions are shown in Appendix C.

12.4.3 Cyclic Plasticity, Shakedown and Ratchetting

When some materials are subjected to uniaxial tension beyond yield, then unloaded and reloaded in compression, it is found that the yield stress in compression is less than the equivalent value in tension. This effect is called the Bauschinger effect and occurs because of the permanent strains and residual stresses remaining after the first yield point is reached. These residual stresses add to the reversed stresses in compression loading thus lowering the second yield point. Figure 12.21(a) shows a typical ‘hysteresis loop’ which forms when reversed loading is applied to a metal, where the Bauschinger effect is exhibited by the fact that the yield stress, σ Y ( c ) in compression is less than σ Y (t ) . In isotropic hardening it is assumed that during cyclic loading in which the load changes from tensile to compressive, the yield point and the effects of work hardening are the same in tension and compression (i.e. no Bauschinger effect), whereas in kinematic hardening the yield point in compression is lower. Referring to Figure 12.21(a), the first yield occurs at point A and then the material hardens up to point B. Upon unloading from point B, the material follows a straight line from B to C and the second yield occurs at point C. With continued load cycles between fixed limits, a stable hysteresis loop may be observed. Cyclic hardening or softening can also be observed in some metals. Under fully reversed constant amplitude stress-controlled experiment, it is observed that the strain amplitude either decreases or increases with each cycle. Similarly, under a straincontrolled loading, the stress amplitude either decreases or increases with each cycle. It is often convenient to represent this type of cyclic hardening or softening by a ‘backbone’ or ‘cyclic’ stress-strain curve, which is drawn through the tips of the stable hysteresis loops, as shown in Figure 12.21(b). The accumulation of plastic strains in cyclic loading is particularly important in ‘low cycle fatigue’ problems, i.e. when the number of load cycles is usually less than 10,000 cycles. The accumulation of strains, which eventually leads to failure, is usually caused by a mechanical or a thermal cyclic load, or a combination of both.

12.43


σ y

σ y

(t)

(c)

(a) Hysteresis loop

(c)

Cyclic stress-strain curve

(d)

Elastic Shakedown

Figure 12.21: Cyclic and reversed loading

Three important phenomena can be observed in cyclic loading when the amplitude is kept constant:

(i)

Elastic Shakedown This occurs when the plastic strain in the cycle is relatively small, i.e. the total strain is less than twice the yield strain (the strain when the stress reaches the yield stress). Referring to Figure 12.21(c), yield first occurs at point A and strain hardening causes the stress to rise to point B. When unloading occurs, the behaviour is linear elastic from B to C, with no further yielding. Therefore, in subsequent cycles, provided the applied load does not go below point C, the material behaves elastically, i.e. it moves up and down the line BC with no further development of plastic strains. Thus the material has ‘shaken down’ to a stabilizing condition, i.e. the structure is assumed to have ‘settled down’ to an elastic state.

12.44


(ii)

Ratchetting Depending on the load level, in some loading situations where a constant amplitude of stress is imposed on the material a stable hysteresis loop may not be reached. Instead, plastic strains keep on accumulating incrementally with each cycle, leading to eventual failure. This mechanism is called ‘ratchetting’, also known as ‘cyclic creep’, and can occur due to a cyclic thermal loading under a constant mechanical load. This occurs due to a cyclic thermal loading under a constant mechanical load. This phenomenon is often observed in materials which exhibit a difference in the yield stress between load cycles of tension and compression, such as cast iron and most composites.

(iii)

Alternating plasticity This phenomenon occurs in some cyclic load situation, where the behaviour can settle down to a state where the plastic strains in ech cycle are equal and opposite, and there is a progressive increase in total strain. It should be noted that in real-life applications, alternating plasticity is of practical concern since a limited amount of incremental material damage occurs in each cycle of reversed plasticity. Such damage can be correlated to the equivalent plastic strain.

12.5

Solution techniques

12.5.1

General

Characteristic features of static non-linear response

Two different types of structural non-linearities have primarily been described in this chapter, namely geometrical and material non-linearities. A third non-linearity is associated with boundary conditions, e.g. so-called contact problems. Characteristic features of various types of nonlinear response are illustrated in Fig. 12.22. The response diagrams illustrate three “monotonic” types of response: linear, hardening, and softening. Symbols F and L identify failure and limit points, respectively. A response such as in (a) is characteristic of glass and certain high strength composite materials. A response such as in (b) is typical of cable, netted and pneumatic (inflatable) structures, which may be collectively called tensile structures . The stiffening effect comes from geometry “adaptation” to the applied loads. Some flat-plate assemblies also display this behaviour initially. A response such as in Fig. 12.22c is more common for structural materials than the previous two. An almost linear regime is followed by a softening regime that may occur slowly or suddenly. Alternative “softening flavours” are given in d – g) 12.45


The diagrams, d – g illustrate a “combination of basic flavours” that can complicate the response as well as the task of the analyst. Here B and T denote bifurcation and turning points, respectively. The snap-through response (d) combines softening with hardening following the second limit point. The response branch between the two limit points has a negative stiffness and is therefore unstable. (If the structures is subject to a prescribed constant load, the structure “takes off” dynamically when the first limit point is reached). A response of this type is typical of slightly curved structures such as shallow arches or shells. The snap-back response (e) is an exaggerated snap-through, in which the response curve “turns back” in itself with the consequent appearance of turning points. The equilibrium between the two turning points may be stable and consequently physically realisable. This type of response is exhibited by trussed-dome, folded and thin-shell structures in which “moving arch” effects occur following the first limit point; for example cylindrical shells with free edges and supported by end diaphragms. In all previous diagrams the response was a unique curve. The presence of bifurcation (popularly known as “buckling” by structural engineers) points as in (f) and (g) introduces more features. At such points more than one response path is possible. The structure takes the path that is dynamically preferred (in the sense of having a lower energy) over the others. Bifurcation points may occur in any sufficiently thin structure that experiences compressive stresses. Bifurcation, limit and turning points may occur in many combinations as illustrated in (g). One striking example of such a complicated response is provided by thin cylindrical shells under axial compression.

12.46


(a)

(b)

Example problem N

(c)

N

N

Formatted: English

(U.S.)

Formatted: English

(U.S.)

Formatted: English

(U.S.)

Formatted: English

(U.S.)

Formatted: English

(U.S.)

Formatted: English

(U.S.)

Formatted: English

(U.S.)

Formatted: English

(U.S.)

k s

(d)

(e) Figure 12.22

Basic response patterns:

(f)

(g)

Characteristic features of nonlinear response: (a) Linear until brittle failure, (b) Stiffening or hardening, (c) Softening.

More complex response patterns:(d) snap-through, (e) snap-back, (f)&(g) bifurcation combined with limit points and snap-back

While in linear analysis the solution always is unique, this may no longer be the case in non-linear problems. Thus the solution achieved may not necessarily be the solution sought. For instance, for a given load R in the case shown in Fig. 12.22d) there may be three different displacement states.

12.47


Non-linear equations

The resultant of internal forces can be expressed as (a i ) T S

R int

i

i

and the total equilibrium can be expressed as R int

R

Hence, the equations that need to be solved are formulated in terms of a total and an incremental equation of equilibrium

∑ (a ) i

T

Si

= R

(12.94a)

i

K I (r )d r

= d R

(12.94b)

For a given external load, R the displacement vector r is sought. Various techniques for solving these non-linear problems exist. Herein three types of methods will be briefly described, namely:   

incremental or stepwise procedures iterative procedures combined methods

The basic principles of these methods will first explained in Sections 12.5.2 – 12.5.4 with reference to problems where the R-r relationship is monotonously increasing with r. In Section 12.5.5 a more general approach to deal with particular problems associated with limit points, will be described.

Physical insight into the nature of the structural problem is essential in addition to pure mathematical knowledge to obtain a successful method of solution. Often a combination of different methods is used to achieve optimal efficiency/accuracy.

12.5.2 Load incremental methods

Incremental methods provide a solution of the non-linear problem by a stepwise application of the external loading . For each step the displacement increment, Δr is determined by Eq. (12.94b). The total displacement is obtained by adding displacement increments. The incremental stiffness matrix, K I is calculated based on the known displacement and stress condition before a new load increment is applied, and is kept constant during the increment. This method is also called Euler-Cauchy method . For load increment No. (m+1) is may be expressed as 12.48


ΔR m +1 = R m+1 − R m Δrm +1 = K I (rm )−1 ΔR m+1 rm +1 = rm + Δrm +1

(12.95)

with the initial condition r0 = 0. In this way the load may be incremented up to the desired level. The method is illustrated for a single degree of freedom in Fig. 12.23.

Figure 12.23 Euler-Cauchy incrementing. It is noted that the method does not include fulfillment of the total equilibrium equation, Eq. (12.94a). For this reason total equilibrium will not be fulfilled. This is illustrated by the deviation between approximate and true K (r)r = R in Fig. 12.23. The accuracy may be increased by reducing the load increment. Also, the load increment should be adjusted according to the degree of non-linearity. Computer programs commonly include procedures for automatical cho ice of load increment.

Example 1

Incremental solution of two-bar problem in Section 12.1.2

In Section 12.1.2 the exact solution to a two-bar problem (for small h/ ratio) was presented. Also, the incremental stiffness, K I was given by explicit formula Eq. (12.9). In practice, analytical solutions can rarely be provided. Rather, numerical solutions need to be used.

12.49


The problem in Section 12.1.2 may be solved by an ”updated Lagrange”, corotational formulation based on the element for the bar in Section 2.9.4, i.e. with a stiffness relationship given by Eq. (2.103). Eq. (2.103) is given for the element in Fig. 2.16 and needs to be transformed to a single d.o.f., r 2 refer to the global degrees of freedom, shown in Fig. 2.16. By establishing the stiffness matrix k and introducing the boundary condition (zero displacement) in the support nodes, a two d.o.f. system (r 1 and r 2) results. Due to symmetry only r 2 is needed to describe the problem in Fig. 12.2. The resulting stiffness matrix for the symmetric system in Fig. 2.16 is R2

2(

EA /c

s

P

2

/c

2 c ) r 2

where P is the axial (compression) force in the member and c = cos θ s = sin θ and θ is the angle between the element axis and the horizontal axis. The trigonometric functions may be approximated by sin

θ =θ−

cos

θ =1−

tg

θ =θ+

θ3

−

6

θ2 2

+

θ3 3

+

θ5 120

θ4 24

−

−

2θ5 15

+

If it is assumed that

2

θ

is negligible relative to 1.0, the following approximations

should be retained when these functions are used in products and sums

sin

θ = θ,

cos θ = 1,

θ=θ

tg

However, when differences between trigonometric functions are calculated, two terms in the expansions need to be retained. In the following, this will only apply to the expression for cos θ. For the system in Fig. 12.2 the angle θ =

α is

small and the incremental stiffness

relationship may be written as:

ΔR = 2(

EA 

α2 −

P 

) Δr

Global equilibrium is described by R

2 Ps

2 Pα 12.50


A solution can be obtained by the following incremental approach: After computational computational cycle (j-1) the following information is available

Displacement: r j Geometry: α j

j 1

Δr i

1

h 1

Axial force: P j

r j

1

h

α 0

j 1

1

i 1

Δ Ri /( 2 sin α j

1)

Next step is calculated calculated as follows: follows:

EA

Displacement increment

ΔR j = 2[

Update total displacement

r j = r j-1 + Δr j

Updated geometry

α j

Updated axial force

P j

h

r j h

=



α 2j−1 −

P j−1 

]Δrj

= K I Δr j

α 0

R j−1 + ΔR j 2α j

The incrementation should be initiated by the geometry corresponding to α 0 and the axial P0 = 0. The load may be incremented by conveniently expressed as a fraction Δλ j of the ultimate capacity of the present system, i.e.

Δ R j

3 3 EAα 0 Δλ j 8

The load incremental factor Δλ is chosen as a number between 0 and 1. The results are shown in Fig. 12.24. Ex.1 based on E= 210 Gpa, A = 0.0001 m2 and Δλ = 0.5 and 0.1. α0 = 0.02. The load incremental factor Δλ

12.51


Fig. 12.24 Load-displacement relations obtained by load incrementation. Incrementation with with equilibrium corrections

A simple improvement of the Euler-Cauchy method can be achieved by an equilibrium correction. Consider the condition after the m’th step – the total load is R m and calculated displacement is rm. However, Eq. (12.94a) is not fulfilled. The unbalanced or residual forces are then given by R r

= ∑ (a i ) T S i (rm ) − R m = R int (rm ) − R m

(12.96)

i

The unbalanced forces may be accounted for by adding them to the next load increment, when rm+1 is to be calculated. This means that the external loads are reduced so that global equilibrium is restored. This principle of equilibrium correction is illustrated in Fig. 12.25 for a single d.o.f. system. Formally the method may be expressed as follows ΔR m+1

= R m+1 - R m

R eq = R m - R int (rm ) Δrm+1

= K I (rm ) -1 ΔR m+1 - K I (r (rm )-1 (R int (rm ) - R m ) = K I (r (rm )-1 ⎡⎣ΔR m+1 + R eq ⎤⎦

rm+1 = rm + Δrm+1

(12.97)

Figure12.25 Euler-Cauchy procedure with equilibrium correction. The additional effort required in this modified Euler-Cauchy method consists in calculating the internal force vector R int int(rm). Euler-Cauchy’s method is based on a single point, rm by calculation of K I, and is denoted one-step method. In principle, multi-step methods are envisaged, but in practice it would n ot be optimal to g o beyond a two-step method. metho d. Example 1- continued 12.52


Equilibrium corrections may be used in the Example 1 in the following manner: The simple load incrementation procedure used in Example 1 may be improved by equilibrium corrections. The unbalanced force, R r r for the two-bar problem may be calculated according to Eq. (12.96). The total load after increment j is then estimated by: j

R j

= ∑ ΔR j i =1

R int int may be calculated as follows Rint (r j )

2 P j sin α j

2 P jα j

when the displacement r j and axial force, P j in each bar is calculated as follows: j

r j

= ∑ Δri ≅ (α 0 − α j ) i

P j

= EA EA ⋅ ≈

EA 2

Δ  / cos

α j

EA

=

(α 02 − α 2j ) =

 / cos

EA 2

αj

(2α 0

( / cos α 0 −  / cos α j ) r j 

−(

r j 

)2 )

Hence R int (rj ) = EA (2α 0

r j 

−(

rj 

r j ) 2 )(α 0 − ) 

Δr i and Δ Ri are incremental r j

= ∑ Δr i = ∑

ΔR i K I (r (ri −1 )

Since K I (r i 1 ) varies there is no simple explicit relationship between r j and R j

Δ Ri

The unbalanced force after increment j is

Rr EA 2α 0

r j 

(

r j 

)2

α 0

r j 

R j

where j

r j

i 1

Δr i ,

j

R j

i 1

Δ Ri

12.53


The geometry and axial force, and hence, the incremental stiffness to be used to calculate the next displacement increment, is the same as described in Example 1 in Section 12.5.2.

Example 2 Incremental solution of elasto-plastic bar problem (after Cook et al., 1988) The purpose is to calculate the load-displacement relationship for the bar in Fig. 12.18a, by using an incremental approach with small but not infinitesimal strains, so that dε becomes Δε. A numerical representation of the stress-strain relation must be stored, so that σ, E and Et can be obtained for any ε. The algorithm outlined below requires that we also store, and update after each computational cycle, the nodal displacements, r element strains ε, and element stresses σ. With two-d.o.f. bar elements, σ and ε are constant over each element length L. 1. For the first computational cycle (j = 1), assume Eep = E for all elements. Apply the first load increment, ΔR 1 . 2. Using the current strains, determine the current Eep in each element. Use Eq. (12.82) to obtain k I for each element i. Obtain K I ( j

1)

the

Δr j

tangent

ΔR j

stiffness

for Δr j . From

matrix,

K I ( j

1)

(a i ) T k it ( j

1)

.

Solve

Δ r j obtain current strain increments Δε ji for each

element (i). 3. Optional. If any elements make the elastic-to-plastic transition, use Eq. (12.83) to revise Eep for each such element, and go back to step 2. Without changing the applied load ΔR j repeat steps 2 and 3 until convergence, which may be defined as

Δε being less than a prescribed fraction of the accumulated total ε in every element. These operations represent secant-stiffness iterations within one of the load steps of the tangent-stiffness procedure. 4. Update: r j

r j

1

Δr j

and for each element, ε ji and σ ji

σ ji

ε ji

1

Δσ ji where Δσ j

Δε ji ( E ep ) j Δ

j

.

For the first cycle (j=1), initial values (subscript j – 1) of displacement, strain, and stress are typically all zero if one starts from the unloaded configuration, but are nonzero if one starts from a state in which plastic action impends. 5. Apply the next load increment and return to step 2. 6. Stop when Δ R j reaches the total applied load. Three cycles of the foregoing algorithm are depicted in Fig. 12.18b. Each cycle produces a line segment whose slope corresponds to the current stiffness. Drift from the exact path can be reduced by using smaller load increments, by exercising step 3 previously discussed, and by using “corrective loads,” which are discussed below. 12.54


Step 3 can be avoided by using load increments ΔR j that bring a single element to the verge of yielding as each load increment is added. This is easily accomplished by scaling the incremental tangent-stiffness solutions. 12.5.3 Iterative methods

The most frequently used iterative method for solving non-linear structural problems is the Newton-Raphson method. The Newton-Raphson algorithm to solve x for the problem: f(x) = 0 is

x n +1

= xn −

f (x n ) f '(x n )

where f '(x n ) is the derivative of f(x) with respect to x, at x = xn. f ( x n ) / tgθ f ( x n ) / f ' ( x n )

f(x) f(xn) θ

f ( xn ) tgθ

=

xn f ( xn ) f '( xn )

Fig. 12.26 Newton-Raphson algorithm This approach can be generalized to solve Eq. (12.92a): In this case K I(r) represents the generalisation of the ∂f / ∂x in Newton’s method for a single d.o.f. See also discussion in Section 12.1.2. Eq. (12.92a) is solved by the iteration formula rn +1 − rn

= Δrn +1 = K I−1 (rn )

(R − Rint )

(12.98)

or rn +1

= rn − K I−1 (rn )

(Rint

− R)

The basic principle for this iteration is illustrated in Fig. 12.27 system.

for a single d.o.f.

This method requires that K I is established and that Δrn+1 is solved from 12.55


R − R int

= K I(n) Δrn +1

(12.98a)

in each iterative step. This is time-consuming. By updating K I less frequently reduced efforts are needed. Since this approach implies only a limited loss of rate of convergence, such modified Newton-Raphson iteration is beneficial. Fig. 12.28 illustrates two alternative for modified Newton-Raphson methods, one with no updating of K I and one method where K I is updated after the first iteration.

Fig. 12.27 Newton-Raphson iteration.

a) No updating of K I

b) K I updated after 1. iteration

Figure 12.28 Modified Newton-Raphson methods for single d.o.f. The iteration is stopped when the accuracy is acceptable. The convergence criterion may be based on the change of displacement from one iteration to the next. The convergence criterion may be expressed by || rn +1

− rn || < ε

(12.99)

where || ⋅ || is a vectornorm and ε is a small, positive number, say, of a magnitude, of the order 10-2 –10-4. The vectornorm is a measure of the size of the vector. There are different vectornorms that may be applied. One alternative is the modified Euclidean norm defined by:

12.56


|| r ||=

1

N

∑ (r / r ) N

2

k

ref

(12.100 )

k =1

where N is the number of components in the vector r and r ref is a reference size, e.g. max(ri ) . N

12.5.4 Combined methods

Incremental and iterative methods are often combined. The external load is applied in increments and in each increment equilibrium is achieved by iteration. Fig. 12.29 illustrated a combination of Euler-Cauchy incrementation and a modified Newton-Raphson iteration.

Figure 12.29 Combined incremental and iterative solution procedures

The procedure is carried out by applying loading according to Eq. (12.95) followed by iteration at each load level by using Eq. 12.98). Commonly a modified NewtonRaphson method is used keeping the gradient K I constant during several iteration cycles. Iteration is stopped according to the criterion ( 12.99). As long as the load curve are monotonically increasing with displacement, the methods described are efficient. However, if there is an extremal point in the loaddisplacement curves (as e.g. in Fig. 12.3) special procedures need to be adopted to achieve acceptable results.

12.57


12.5.5 Advanced solution procedures General

The solution procedure described so far are a combination of incremental load coupled with full or modified Newton-Raphson iterations. Because the plastic flow rules are incremental in nature elasto-plastic problems should strictly be solved using small incremental steps. For, no matter how accurately flow rules and keeping on the yield surface may be satisfied within an increment, the solution is only in equilibrium at the end of each increment after equilibrium iterations. However, often acceptable solutions can be obtained with large steps. Although incremental-iterative techniques provide the basis for most nonlinear finite element computer programs, additional sophistications are required to produce effective, robust solution algorithms. An extensive of more refined methods are discussed e.g. in Chapter 9 of Crisfield (1991). In this section a brief review of such methods is given. For instance, severe difficulties are encountered when load incrementation methods are used to pass a limit point, L (Fig. 12.22c),i.e. when the target stiffness becomes zero. Using incrementation in displacement instead of load can solve this problem. This approach will be effective for problems characterised by Fig. 12.22(c-d) for which the load is uniquely determined by the displacement. However, displacement incrementation will fail at turning points (T) (“snap-back point”) e.g. in Fig. 12.22e). In the present section emphasis will be placed on arc-length techniques for solving these problems. Prior to their introduction, analysts either used artificial springs, switched from load to displacement control or abandoning equilibrium iteration in the close vicinity of the limit point. In relation to structural analysis, the arc-length method was originally introduced by Riks [1972] and Wempner [1971] with later modifications being made by a number of authors. Before describing such methods, one may ask why we need to pursue the response beyond a limit point (L) in Fig. 12.22c). After all, the limit point represents the ultimate strength. There are several reasons: i)

In many cases it may be important to know not just the collapse load, but whether or not this collapse is of a “ductile” or “brittle”nature.

ii)

The structure with the characteristic displayed in Fig. 12.22 may represent a component in structure. The ultimate behaviour of a redundant structure consisting of such components, would depend upon the post-ultimate beyond limit point, L) behaviour of the component.

Method

As a starting point the global equilibrium equation is written as: 12.58


g(r, λ )

R int (r )

λ R ref

0

(12.101 )

where R ref is a fixed external load vector and the scalar λ is a load level parameter. Equation (12.101) defines a state of “ proportional loading ” in which the loading pattern is kept fixed. Non-proportional loading will be briefly mentioned later in this section. The essence of the arc-length method is that the solution is viewed as the discovery of a single equilibrium path in a space defined by the nodal variables, r and the loading parameter, λ. Development of the solution requires a combined  incremental (also called predictor )  iterative (also called corrector ) approach.

Many of the materials (and possibly loadings) of interest will have path-dependent response. For these reasons, it is essential to limit the increment size. The increment size is limited by moving a given distance along the tangent line to the current solution point and then searching for equilibrium in the plane that passes through the point thus obtained and that is orthogonal to the same tangent line (Fig. 12.30c). In Fig. 12.30c the arc-length control strategies in the solution of nonlinear equations are illustrated and compared with load and displacement control. For instance if load incrementation is applied, the iterations are carried out to correct the displacements. When the arc-length method is applied the itereations are carried out with respect to both the load and displavements.

Fig.12.30 Geometric representation of different control strategies of non-linear solution methods for single d.o.f. An increment is made along a tangential path, SP. Correction to reach g = R int - λR ref = 0 is obtained by iteration controlled by the (hyper)plane c = 0 a) load control, b) state control, c) arc-length control

The arc length is formulated as an additional variable involving both the load and displacement. The increment in the load-displacement space can be described by a displacement vector Δr and a load increment parameter Δλ , such that ΔR = Δλ R ref . This formulation results in an additional equation to be solved. The advantage of the 12.59


extra equation is that the solution matrix never becomes ‘singular’ even at the limit points. Therefore, the solution matrix is re-assembled with N+1 variables, where N is the total number of the variables (degrees of freedom) of the system. However, the disadvantage is that, in some FE formulations, the solution matrix becomes unsymmetric, which may incur an increase in computing time and/or computer storage, particularly for very large problems. First the increment (predictor) from the “First point” is made along the tangent. Then, this solution is corrected iteratively to reach the “Second point” and so on. Several methods exist to obtain the arc length, for example by making the iteration path follow a plane perpendicular to the tangent of the load-displacement curve, as shown in Figure 12.31. Alternatively, instead of a normal plane, more sophisticated paths such as spherical or cylindrical planes can be followed, and the solution matrix can be manipulated to become symmetric (see, for example, Crisfield [1991]).

Load factor λ iteration

increment

Figure 12.31: Schematic representation of the arc-length technique. A geometrical interpretation of the incremental iterative approaches by Riks-Wempner and Ramm is sketched in Fig. 12.32. While in Ramm’s method the iterative corrector is orthogonal to the current tangential plane during the iteration, it is orthogonal to the incremental vector ( Δr0 , Δλ 0 ) in the Riks-Wempner methods.

a) Riks-Wempner’s method

b) Ramm’s method

Fig. 12.32 Arc-length control methods (Crisfield, 1991) 12.60


An alternative iterative method is so-called orthogonal trajectory iterations (Fried, 1984). The first step in this method can be illustrated by reference to Fig. 12.30. The first iteration is then assumed to be orthogonal to the vector S’P’ instead of SP. The resulting iterative solution will appear as shown in Fig. 12.33. Haugen (1994) found that this method was more efficient than the normal plane iterations.

Figure 12.33 Arc-length method with orthogonal trajectory iterations.

Automatic incrementation

To achieve computational efficiency the load increment Δλ should be chosen depending upon the degree of nonlinearity of the problem. Methods have been established based on the curvature of the nonlinear path (den Heijer and Rheinboldt, 1981) or the so-called current stiffness parameter (Bergan et al, 1978): S

i p

Δr1T ΔR 1 Δλ i2 = Δλ12 ΔriT ΔR i

(12.102)

S pi refers to increment No.i.

The initial value of S pi ( S p1 ) is 1.0. For stiffening system it will increase. For softening system it will decrease.

If S pi changes sign the sign of the increment should be

changed. Numerical experiments show that nearly the same number of iterations were requested to restore equilibrium when the increments were chosen according to the approach of Bergan et al.( 1978).

12.61


Ramm (1981) proposed another approach for estimating the necessary increment Δλ (load incrementation) or Δλ (for arc-length method). The new arc-length,  n is obtained by

Δ n

Δ 0

I d I 0

1/ 2

(12.103)

where Δ 0 is the “old” arc-length, and Id and I0 are the desired number of iterations (given as input) and the number of iterations when the old arc-length was used. This approach requires a suitable estimate of the initial arc-length. An alternative tactic is to apply load incrementation for early increments and switch to arc-length control once a limit point is approached. The current stiffness parameter can be used to decide the switch from load incrementation (or displacement control) to the arc-length method. An alternative indicator of when the limit point is approached is the check of negative values on the diagonal of the incremental stiffness matrix, i.e. negative pivot elements in the solution algorithm. In particular the current stiffness parameter may be used to control the solution strategy at limit points or bifurcation points. Alternative changes may be made when the current stiffness is below a limit value, namely - the sign of the incrementation is changed - iteration may be suppressed and a simple incrementation may be used. Iterations are then resumed when S pi increases beyond a specific limit (see Fig. 12.34).

Fig. 12.34 Possible choice of solution algorithm for a problem with limit point

Non-proportional loading

The solution procedures in this chapter have been based on the equilibrium relationship of (12.101) which implies a single loading (or displacing) vector, R ref , is proportionally scaled via λ. For many practical structural problems, this loading 12.62


regime is too restrictive. For example, we often wish to apply the “dead load” or “self-weight” and then monotonically increase the environmental load. Even more general load conditions may be required. Fortunately, many such loading regimes can be applied by means of a series of loading sequences involving two loading vectors, one that will be scaled (the previous R ref ) and one that will be fixed ( ( R ref ) . The external loading can then be represented by R R ref

λ R ref

(12.104)

so that the out-of-balance force vector becomes g

R int

R ref

λ R ref

(12.105)

12.5.6 Direct integration methods General

Up to now the methods for directly solving the statistic nonlinear equation have been based on incrementation of loads or displacements. Possibly combined with iterative methods. These are often considered standard methods for solving nonlinear problems (e.g. in ABAQUS). An alternative approach is to use so-called finite difference methods for direct integration of the dynamic equation of motion : M r (t) + Cr (t) + Kr(t) = R(t)

(12.106)

to solve the static problem : K r = R . Nonlinear structural effects make K a function of r, K (r ) .This means that the loading R is increased (artificially) or as a function of time. The loading time needs to be sufficiently long so that the inertia and demping forces do not have an effect on the behaviour on the static problem that is to be solved. A finite difference approximation is used when the time derivatives of (12.106) r and r ) are replaced by differences of displacement (r) at various instants of time. ( The direct integration methods are alternatives to modal methods, and they can be used to successfully treat both geometric and material non-linearities. The finite difference methods are called explicit if the displacements at the new time step, t + Δt , can be obtained by the displacements, velocities and accelerations of previous time steps. r (t

+ Δt ) = f {r (t ), r (t ), r(t ), r(t − Δt ), r(t − Δt ), r(t − Δt ),...}

(12.107)

or

ri +1

= f {ri , ri , ri , ri−1 , ri−1 , ri−1 ,...} 12.63


This is as opposed to the implicit finite difference formulations where the displacements at the new time step, t + Δt , are expressed by the velocities and accelerations at the new time step, in addition to the historical information at previous time steps. ri +1

= f {ri +1 , ri +1 , ri , ri , ri ,...}

(12.108)

Many of the implicit methods are unconditionally stable and the restrictions on the time step size are only due to requirements of accuracy. Explicit methods, on the other hand, are only stable for very short time steps. Central difference method

To illustrate this approach, one of the explicit solution methods, the central difference method is described in the following. The central difference method is based on the assumption that the displacements at the new time step, t + Δt , and the previous time step, t − Δt ,can be found by Taylor series expansion. t2 t 3 Δ Δ   ri +1 = r0 (t ) + Δt ri + ri + ri + ...

2

ri −1

= ri − Δt ri +

Δt 2 2

 ri

6

−

Δt 3 6

 ri

+ ...

(with r0 (t ) = ri )

(12.109)

(12.110)

The terms with time steps to the power of three and higher are neglected. Subtracting Eq. (12.110) for Eq. (12.109) yields : ri +1 − ri −1

= 2Δt ri

(12.111)

Adding Eq. (12.109-110) yields : ri +1 + ri −1

= 2r + Δt 2ri

(12.112)

Rearranging Eq. (12.111-112), the velocities and accelerations at the current time step can be expressed as: ri

=

 ri

=

1

{ri +1 − ri −1}

(12.113)

{ri +1 − 2ri (t ) + ri−1}

(12.114)

2Δt

1

Δt 2

Finally inserting Eqs. (12.113-114) into the dynamic equation of motion Eq. (12.106) gives:

12.64


⎧ 1 M + 1 C⎫ r = R (t ) − K ⎨ 2 ⎬ i +1 i 2Δt ⎭ ⎩ Δt

ri (t ) +

1

Δt 2

M {2ri

− ri −1} +

1 2Δt

C ri −1

(12.115) If the mass matrix, M, and the damping matrix, C, are diagonal, the equations will be uncoupled, and the displacements at the next time step, t + Δt , can be optained without solving simultaneous equations. The characteristic features of Eq. (12.115) are best illustrated by an example. Let us consider a system with three global directions of freedom. The mass matrix, M and damping matrix Care assumed to be diagonal. Eq. (12.115) may then be written as:

⎧ ⎡ M 11 0 0 ⎤ ⎡C11 0 0 ⎤ ⎫ ⎡ r1( i +1) ⎤ ⎡ R1( i) ⎤ ⎡ K11 K12 1 ⎢ ⎪ 1 ⎢ 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎥⎪ ⎢ 0 C22 0 ⎬ ⎢ r2( i +1) ⎥ = ⎢ R2( i) ⎥ − K21 K22 M 22 0 + ⎨ 2⎢ ⎥ ⎢ ⎥ ⎢ ⎪ Δt ⎢ 0 0 M ⎥ 2Δt ⎢ 0 0 C ⎥ ⎪ ⎢ r ⎥ ⎢ R ⎥ ⎢ K K 33 ⎦ 33 ⎦ ⎭ ⎣ 3( i +1) ⎦ ⎣ ⎩ ⎣ ⎣ 3( i) ⎦ ⎣ 31 32 ⎡ M 11 0 0 ⎤ ⎧ ⎡ r1( i ) ⎤ ⎡ r1( i −1) ⎤ ⎫ ⎡C11 0 0 ⎤ ⎡ r1( i −1) ⎤ ⎪ ⎢ ⎥ ⎢ ⎪ 1 ⎢ 1 ⎢ ⎥ ⎥ ⎥ ⎥⎢ + 2 ⎢ 0 M 22 0 ⎥ ⎨2 ⎢ r2( i ) ⎥ − ⎢r2( i −1) ⎥ ⎬ + 0 C22 0 ⎢r 2( i −1) ⎥ ⎢ ⎥ 2Δt Δt ⎢⎣ 0 0 M 33 ⎥⎦ ⎪⎩ ⎢⎣ r3( i) ⎥⎦ ⎢⎣ r3( i −1) ⎥⎦ ⎪⎭ ⎢⎣ 0 0 C33 ⎥⎦ ⎢⎣ r 3( i −1) ⎥⎦

⎤ ⎡ r 1( i) ⎤ ⎥⎢ ⎥ K32 ⎢r 2( i) ⎥ ⎥ ⎥ K33 ⎥⎦ ⎢⎣ r 3( i) ⎦ K13

(12.116)

The first equation in Eq. (12.121) is explicitty written as :

⎧ 1 M + 1 C ⎫ r = R (t ) − K r − K r − K r ⎨ 2 11 11 ⎬ 1( i +1) 1 11 1( i) 12 2( i) 13 3( i) 2Δt ⎭ ⎩ Δt +

1

Δt 2

{

M 11 2r1( i )

− r1( i −1) } +

1

(12.117)

C11 r 1( i −1) 2Δt

This shows that r i ( i +1) can be directly, explicity determined by the response at time t. There is no coupling between displacements, r j (i +1) at the time t + Δt . Because the expressions for the displacements are explicitly given, there is no need to invert the tangent stiffness matrix at every time step. The explicit method also has the advantage of drastically reducing the need for computer memory capacity. The stiffness forces, or internal force vector, can be found by summation of element contributions. The global stiffness vector, K, need not to be stored in the computers core memory. As already mentioned, Eq. (12.115) is conditionally stable and requires that

Δt <

2 ω max

(12.118)

12.65


where ω max is the highest natural frequency of det(K − ω 2 M ) = 0

(12.119)

The maximum frequency of Eq. (12.119) is bounded by the maximum frequency of the constituent unassembled and unsupported elements. When finding the maximum natural frequency of an element, one will see that the time step, Δt , must be short enough that information does not propagate across more than one element per time step. The maximum allowable time step will therefore be limited by a characteristic length, λ e , of the element and the acoustic wave speed, c.

Δt <

λ e c

(12.120)

Higher order elements yield higher maximum frequencies and should be avoided when doing explicit integration. Many alternative methods exist. See. e.g. Solution of static problems

The explicitmethod is very well suited to treat dynamic problems. As indicated above the method can also be used to solve static problems. It is obvious that the periode of the loading, or the amount of time for the loading to reach it’s maximum value, must be much larger than the largest eigenperiod to avoid dynamic effects as determined from the lowest eigenfrequency found for Eq. (12.119). The response of the structure is also dependent on the magnitude of the loading, not only on the period, and this complicates the picture. In addition, failures due to collapse or cracking of parts of the structure will cause vibrations. These events will not be captured by a traditional static analysis. All effects taken into account; if the time of the loading to reach its maximum level is conservatively chosen to be 30 times the longest eigenperiod of the system, the dynamic effects have shown to be negligible. Another problem with explicit analyses is that post-collapse behaviour cannot be traced if the loading is given as applied forces. In many cases this can be avoided by switching to displacement control. If displacement control is not possible or desirable, implicit solution procedures using arc length solution methods can be used. An advantage with the explicit solution procedure is that it is very easy to use. The user of an explicit finite element program is left with the difficulty of applying loads sufficiently slowly to avoid dynamic phenomena and sufficiently fast to avoid too large computional efforts times. In static analyses, and even in dynamic analyses, the computational time van be considerably reduced by changing mass densities in elements. The time step will be governed by the smallest element in the model. Artificially increasing the mass of small elements will reduce the acoustic wave speed and hence allow longer time steps. Similarly very large elements can be given mass reduction and hence be less affected by inertia forces. Systematic increase and reduction of element masses can be 12.66


performed to improve computational efficiency, but the details in these methods will not be elaborated on.

12.6

Applications

12.6.1 General

In practical structural analysis it is good practice that a nonlinear analysis is preceeded by a linear analysis, that provides a basis for planning t he nonlinear analysis. The linear analysis may be e.g. a stress analysis or buckling analysis. The stress analysis provides information about which nonlinearities which may be important, i.e. where stresses exceed yield limit and implies that elasto-plastic behaviour should be considered. In case of compressive loading, buckling may be an important issue. A linear buckling analysis based on elastic and geometric stiffness matrices, then provides buckling loads and buckling modes. Knowledge about the elastic buckling load (“Euler load”) can, together with code formulations for ultimate strength yield information about the ultimate capacity to be expected. The buckling modes are crucial for the introduction of initial imperfections for the nonlinear analysis. Usually, a linear combination of linear buckling modes are used to describe the imperfection pattern, and the maximum value is scaled to some given imperfect ion tolerance level.

12.6.2 Beams and frames General computer programs

General purpose finite element programs for nonlinear analysis, such as ABAQUS, LSDYNA etc. contain several options for beam models. Special purpose programs

The program system USFOS was initially developed at NTNU/SINTEF to deal with nonlinear analysis of frames and trussworks. Later it has been extended also to cover plates and shells. The basic for USFOS is summarized in this section. More details may e.g. be found in Skallerud and Amdahl (2002). USFOS may be used to assess the ultimate global capacity of space frame structures and to document the residual strength of such structures. The formulation of beam models allows the use of very coarse finite element modelling of the structure, but still obtain results good accuracy. Local flexibility of tubular joint is included through a simplified, but very efficient formulation. The formulation gives very good results compared with shell analysis of the joint, but requires no special modelling of the joint 12.67


geometry. The flexibility characteristics are calculated form the information already present in the design model o f the structure. The program is based on an updated Lagrange (incremental-iterative) procedure. It uses a (nonlinear) Green strain formulation (Section 12.3.3) . Thus, the USFOS beam element is valid for large displacements, but restricted to moderate strains. The influence of axial force on the bending stiffness of the element is introduced by the nonlinear terms in the Green strain formulation. The tangent and secant stiffness matrices are then obtain by introducing interpolation functions for the element displacements. The shape functions in USFOS are taken as the exact solution of the 4th order differential equation for a beam column, i.e. the stability or Livesly functions (Section 12.2 and Appendix A). With these shape functions all integration in the element stiffness expressions is carried out analytically, giving closed form solutions for the nonlinear elastic stiffness matrix. In addition to accounting for the coupling between axial displacements and lateral rotations on the element level, large displacement effects are also taken into account by updating the local reference co-ordinate system at each step. Material nonlinearities are modelled by plastic hinges at element ends or at element midspan. The control of element plastification is a control of stress resultants (forces and moments) against the total plastic capacity of the cross section, instead of local stresses. Possible unloading of plastic hinges into the elastic range is checked at each load step. Elasto-plastic element properties are established based on the flow theory. The plasticity formulation is based on interaction formulae for plastic capacity of stress resultants.

Fig. 12.35 USFOS formulation: One finite element per structural element.

12.68


Hinges may occur at beam-ends and at midspan. In the latter case the original element is subdivided into two sub elements. The extra nodal point is introduced automatically and eliminated by static condensation before adding into the global stiffness matrix. The difference between yield hinge formulations have been investigated (Hellan et al. (1994 -1995), and show significant differences. Elastic-perfectly plastic models over predict the column buckling capacity of ideally straight tubulars. For yield-hinge models incorporating first fiber yield, gradual plastification and strain hardening, the buckling capacity depends on the plasticity parameters given to model the transition from elastic to plastic behaviour. In general, these formulations may be slightly conservative for stocky columns, and slightly un-conservative for slender columns. For both these formulations, geometrical imperfections and residual stresses for a member can be accounted for by introducing equivalent initial imperfections in the element formulation to ensure exact fit to any given column curve, e.g. Eurocode 3, which is based on extensive test results. However, particular assessment of the imperfection mode for systems analysis is required. By comparing various imperfection patterns, it was found that imperfection in the direction of the global load yields the lowest capacity. This is discussed in more detail in Hellan et al. (1994). Joints can in principle be modeled by shell elements. Systematic FE studies of linearly elastic joints have resulted in parametric formula for the joint flexibility (Ultiguide, 1999). Under extreme loads, the nonlinear deformations of the joint and failure characteristics can influence the disposition of forces and the overall structural response. Failure of tubular joints generally involves some combination of the following local and global modes : 

Local plastic deformation (yield) of the chord around the brace intersection  Cracking in the chord at the weld toe (and propagation to severance)  Local buckling in compression areas of the chord  Ovalisation of the chord cross-section  Beam shear failure across a gap K joint chord  Beam bending of the chord especially for T/Y action The specific response depends upon the type of joint (T/Y, X, K; simple, stiffened, grouted; etc), the loading (Axial – tension/compression; bending – in-plane bending/out-of-plane bending, etc) and the joint geometry parameters ( β , γ , etc ). Fig.12.36b illustrates typical load-deformation responses for axially loaded joints as seen in isolated component tests performed with idealized boundary conditions.

12.69


a) member behaviour

b) Joint behaviour

c) Joint model

Figure 12.36: Typical behaviour of members and joints in jacket platforms Several methods are in use to model joint behaviour in non-linear analyses as summarized in Ultiguide (1999). One method is to model the joint behaviour by a nonlinear spring between the brace end and the center of the chord. Spring properties can then determined from test data, FE analyses or parametric formulae. Alternatively, the springs can be replaced by a beam-column element with the ultimate capacities determined from the capacity equations for the joint (Fig. 12.36c) given in appropriate codes. The joint behaviour may be modeled by a plastic potential, with interaction between the axial force, in-plane bending and out-of-plane bending. It is, however, an open question whether the interaction in the post-ultimate range is adequately described in this way. Formulations have also been published that accounts for braceto-brace interaction by adding ‘beam’ elements between the brace ends. It is important that the nonlinear analysis method can predict failure in accordance with recognized failure criteria. Hence, besides modeling of stiffness and capacity, yield characteristics, post ultimate behaviour and ductility limits (if applicable) should be represented, including local failure modes such as local denting, local buckling, joint overload, joint fracture, etc. Local buckling may be accounted for by various approaches. One alternative is to modify the stress-strain relationship by a drop in the stress when the critical strain εcrit, is reached (Fig. 12.37a). Another alternative is to simulate local buckling by using a model for a dented tube (Taby and Moan, 1985). This is done by assuming an initial dent and accounting for the growth of the buckle as the member deforms (Fig. 12.37b). The same model is applied in the analysis of structures with dents imposed by e.g. accidental loads. General shell FEM formulations may be able to capture local buckling through the shell modeling, but initial geometric imperfections are generally needed to initiate the local buckle correctly. Obviously, a global analysis based on a model which accommodates local buckling modes by shell elements will be very time consuming.

a) Equivalent stress-strain

b) Dented member model

Figures 12.37: Models for damaged members 12.70


An automatic load algorithm is implemented. At each step the occurrence of new plastic hinges is examined. The step length is scaled to make the stress resultants satisfy the yield criterion exactly. Equilibrium iterations may be performed including correction to the yield surface due to drift-off of stress resultants during finite increments. Possible unloading of plastic hinges into the elastic range is checked during the load incrimination step. If unloading occurs, the element stiffness is modified and the system stiffness recalculated before the step is completed. For a detailed description of theory and implementation, reference is made to the theory manual of the program (Søreide et al., 1994). Moan et al (1985) compared results obtained by the USFOS formulations with those obtained by “conventional” a FE formulation (Engseth, 1986). Ultimate limit state checks of structures have traditionally been carried out by linear global analyses to calculate member forces; and only include inelastic and second order effects in the determination of the ultimate capacity of the individual components. The strength of materials formulations used in codes for component limit states have been extensively validated against test results. This methodology hence, focuses on the first failure of a structural component and not the overall collapse of the structure, which obviously is of main concern in view of the failure consequences. Tools for nonlinear analysis that includes second order geometrical and plasticity effects, p rovides a means to account for possible re-distribution of the forces and subsequent component failures until system collapse, and, hence, a more realistic estimate of the ultimate capacity. Initially, such methods were developed for calculating the residual strength of systems with damage (progressive limit state checks of NPD). More recently, such methods are also applied for re-assessment of ageing structures to determine the consequences of fatigue induced fracture of members in connection with inspection planning. Examples are presented in dr.ing.thesis by Ø. Hellan (1995) and Hellan et al (1994).

Fixed platform analyses are carried out by modeling the pile-soil behaviour by equivalent linear or nonlinear concentrated springs or, distributed springs along the piles, or continuum (finite element) model that represent stiffness and foundation, capacity, appropriately using the material properties in the different soil layers. Most of the existing pile-soil models are empirical or semi-empirical based on the acquisition of data from a limited number of large diameter pile tests from early 1970’s until now. Although these empirical models have provided practical tools for the designer of the offshore pile-soil foundation, they are associated with some uncertainties, see e.g. Horsnell and Toolan (1996) and Nadim and Dahlberg (1996). Soils exhibit nonlinear behaviour, even at low load levels, which needs to be accounted for. The response calculations are carried out as integrated pile-soilstructure analyses. A typical finite element model of the platform in Fig. 12.38a is shown in Fig. 12.38c. 12.71


(261) (363) (463)

(455) (456)

a)

b)

c)

d)

Figure 12.38: Global ultimate behaviour of a North Sea jacket The effects of global seabed scour and local scour in granular soils, and the partial loss of soil-pile contact in cohesive soils should be accounted for. When modeling the individual piles in a pile group, nonlinear soil P-y and T-z curves have to be adjusted to account for pile group effects. Cyclic loads cause deterioration of the lateral bearing capacity as indicated in Fig. 12.38. The soil capacity and the nonlinear P-δ characteristics given in most codes represent the fully degraded properties of the soil, based on the cyclic (hysteretic) behaviour. The capacity under monotonic, static loading can be significantly higher, as shown in Fig. 12.38. High loading rates (compared to laboratory test loading) can also contribute to increase the capacity up to 40% for lateral loading and 50% for axial loading for wave loading. It is noted that this strength increase only applies to the dynamic (variable) part of the soil reaction.

Analysis strategy and solution procedure

To facilitate the interpretation and verification of ultimate strength predictions a systematic approach to the non-linear analysis is necessary. For each scenario a 12.72


sequence of analyses of increasing complexity and increasing realism should be performed, beginning with the structure alone; then introducing, for example, soilstructure interaction and non-linear joint characteristics in turn, as appropriate. Also, a static analysis is done before a dynamic one, and a monotonous analysis before a cyclic one. At each stage the results should be examined and the influence of the modeling changes examined. Various options for solution of the non-linear problem are used. For static problems an incremental iterative, proportional loading algorithm is applied. Load increment size is adjusted when plastic hinges are formed, predefined maximum displacement increment is exceeded and when load limit points are passed. In the latter case arc length iterations is used as one option. Even though automatic solution controls are generally available, manual intervention may be required. Sudden changes in stiffness, rapid unloading and alternative equilibrium conditions in the non-linear region may cause numerical instability and present challenging problems for all analysis software; small step sizes may be required to coax the solution. Apparent solution difficulties should be investigated in terms of the non-linear events being predicted. Failure of an analysis does not necessarily mean that the ultimate strength of the system has been reached. Solution should be continued until a clearly defined peak or sustained limit load has been reached and the post-ultimate response characteristics determined. Case studies

As demonstrated e.g. by Moan, Azadi and Hellan (1997) the choice of pile-soil model can affect the load distribution in the structures and, hence, the failure mode and corresponding ultimate strength. The most important issue is, of course, that a pure linear pile-soil model would not represent a possible soil failure and, hence overestimate the system strength if the pile-soil is the critical part of the system. For the jacket in Fig. 12.37a with plugged piles the pile-foundation is not critical. Yet the difference in jacket failure mode when using a linear instead of a nonlinear model, results in an ultimate load which is about 15% smaller for the former case. Determination of the global ultimate capacity by monotonically increasing wave loading (pushover analysis) has become a well-established approach. As mentioned above, systems analysis is particularly of interest to demonstrate robustness in connection with structures that are damaged due to accidental loads of fatigue fracture. It is found that the nonlinear approach especially yields more realistic predictions of the ultimate strength of a damaged system than a conventional approach based on a linear global analysis and an ultimate strength assessment of the structural component. To illustrate the effect of damage the jacket in Fig. 12.38 is considered. The ultimate capacity is normalized with respect to F100. For this case a linear spring model is applied for the pile-soil. Damage in terms of removal of individual braces is considered as indicated in Fig. 12.38b. Table 12.2 shows that failure of the braces 261 and 463 for broad side loading does not reduce the ultimate strength, and, most importantly, for all cases with a single brace failure the reserve capacity is at least 12.73


0.7 x 2.73 = 2.08 times the 100 year characteristic load, while the normal total safety factor for design checks of components of offshore structures is about 1.5. Table 12.2: Residual strength of damaged North Sea jacket. Linear pile-soil behaviour. Wave height incremenation Loading and damage condition

Broad side loading End-on loading Brace 261 263 Brace Brace 463 Brace 455 Brace 456 Ultimate strength Fult / FH100 2.73 Residual strength Fult (d) / Fult 1.0 Note: Fult

2.73

2.73

2.89

2.89

0.76 1.0 0.91 0.85 – ultimate strength, Fult(d) – ultimate strength of damaged platform

Conclusions

For a realistic evaluation of the structural safety of offshore structures with complex three-dimensional and stochastic sea loads and complex structural geometries, the global ultimate strength needs to be calculated, yet, considering relevant local failure modes of members and joints. The platforms have to be analyzed in intact and damage condition. The efficiency of the calculation methods is then detrimental, especially in a design context, when repeated analyses are required. Members are modeled as beam elements and joints by spring elements. It is shown how local failure modes, especially local buckling of members and joints, are accommodated by a phenomenological approach. The sensitivity of the global ultimate strength to damage of individual members is examined. 12.6.3 Plane stress and bending of plates and shells Formulations

The one-dimensional formulations in terms of strain, Exx, stress, Sxx or their incremental form, need to be generalized for multidimensional problems. The basis for this generalization is given in appendices. To a large extent this generalization corresponds to the generalization of the linear theory from one – to e.g. twodimensional problems. For instance, the expressions for strains and stresses when geometrical nonlinearities are considered, can be established by introducing “correction terms” on the linear strain expression in the same way as for Exx as discussed in Section 12.3.3. However, the generalization of elasto-plastic theory from one-to e.g. two-dimensional problems include the following issues:  yield criterion (von Mises)  hardening rule  flow rule 12.74


The von Mises yield criterion is well known. The generalization of isotropic hardening rule is relatively straight-forward (see Appendix C), other hardening rules are more complex and are not described in the Appendix C. The generalization of the flow rule is the most important issue in the generalization. The main starting point is Eq. (C.15) which follows from Drucker’s postulate. Eqs. (C.19) shows that the plastic strain increment is proportional to the equivalent stress increment, d σ and the deviatoric stress, sij. Eq. (C.30) expresses the stress increment as a function of total strain increment. It is this relationship that is used in multidimensional elasto-plastic finite element formulations.

Elasto-plastic analysis of stresses and strains in a corner of a large frame composed of box members (Moan and Nordsve, 1979)

The purpose of this example is to study the behaviour of a corner in a large frame (Fig. 12.38a). The analysis has to relevance to serviceability and ultimate limit state dimensioning criteria. The structure analyzed is a plane idealization of the structure shown in Fig. 12.39c. Both an unstiffed configuration (No. I) and a stiffened one (No. II) are considered. In configuration No. II only the effective part of traversal panels and flanges of the actual panel are included, see Figs. 12.39c-d. Buckling stiffeners are smeared and contribute to the effective plate thickness. The plate material is assumed to be ideally elasto plastic with yield stress σ Y = 360N/mm2 . The loads are introduced according to simple beam theory on the boundary AA, see Fig. 12.39d. The level of the load is expressed by the nominal stress level according to beam theory in section BB. The external moment and shear forces at AA are such that the maximum equivalent stress in the outer fiber in BB is η ⋅ σ Y . The average shear stress in BB is 0.18 η ⋅ σ Y and the average normal stress in CC 0.58 η ⋅ σ Y . η is a usage factor. Two load conditions are considered, namely η = 0.63 (Condition A) and η = 0.84 (Condition B). These usage factors correspond to typical values allowed for the operational and extreme load condition for marine steel structures.

12.75


a) Gravity platform

c)Corner of large frame

b) Deck structure and columns

d) Plane stress symmetric model. Configuration I: unstiffed Configuration II: stiffened

e) Mesh of constant strain elements and (possible) bar elements, and uniaxial model of material. Fig. 12.39. Idealization of a corner in a large frame composed of box members.

12.76


The finite element solution is obtained with an idealization of the panel based on 288 constant strain plane stress elements. Out-of-plane displacements are not considered. Stiffeners are represented by bar elements, totaling 77. The resulting mesh is shown in Fig. 12.39e. Because the nonlinear behaviour is very local in the present case, it is computationally advantageous to use the so-called substructure approach. The two substructures used, denoted by No.1 and 2, are shown in Fig. 12.39e. A selection of the results obtained in this study is displayed in Fig. 12.40 and Table 12.3. Fig. 12.40 shows how the yield zone develops in the two configurations. It is observed that the initial extent of the yield zone for Configuration I can be estimated fairly well by using the yield condition in conjunction with stresses calculated according to the elasticity theory. However, the gradual local softening of the structure due to plastification influences the further development of the yield zone.

Fig. 12.40. Yield zones in substructure No.2, see Fig. 12.39e. In a plastic zone, strain is obviously a better measure of strength than stress. Table 12.3 displays the maximum strains found. For both configurations considered the maximum strain is far below fracture strain, which may be of the order of 10ε Y . This is because the plastic flow is contained in a fairly large elastic medium. The location HS1 and HS2 are the most highly strained zones in Configuration I and II, respectively. In particular the shear straining at HS2 in Configuration II should be noted. The residual equivalent stress when unloading from maximum loading is displayed in Table 12.3.

12.77


Table 12.3: Characteristic stress and strain responses in corner of a large frame.

Elastic analysis Configuration

I II

Elastic stress concentration factor: σ max / σ ref 1.91 1.78

Elasto-plastic analysis referring to max. loading (η = 0.84) Max. equivalent Max. residual equivalent stress plastic strain: after unloading : p ∗ ε max σ max p / ε Y ε max

∗ / σ Y σ max

1.6 2.1

0.62 0.52

σ max

maximum pointwise equivalent stress

σ ref = ησ Y

maximum equivalent stress in Section BB, see Fig. 12.39d.

∗ σ max

maximum pointwise equivalent residual stress.

σ Y

yield stress

ε Y

strain in uniaxial specimen at the yield point.

p ε max

maximum pointwise accumulated equivalent plastic strain.

The above case study may serve as a basis for a discussion of the adequacy of linear elastic stress analysis and possible other, improved methods of analysis to represent actual limit states for a steel structure. The traditional way of assessing the strength of a structure is by ensuring that the stresses obtained in a linear elastic analysis are less than the yield stress. Commonly the analysis is performed by means of a beam or a frame model, which represents the gross stress variation. If the stress analysis is based on a finite element method, the very local stresses can be calculated. Therefore, due to the high elastic stress concentration present in most steel structures, the stresses obtained in a linear elastic analysis may exceed the yield level even or operational loads. In the above case study the maximum equivalent stress for Configuration I and II at Load Level A (η = 0.63) is, respectively, 1.12σ Y and 1.20σ Y . This fact does not necessary imply an inadequate design. This is because ductile collapse strength depends on the average stress over a region and not the point wise maximum value – given that sufficient ductility allows a 1 redistribution of stresses to comply with the stress distribution at ultimate collapse. However, fracture strength under single overload or repeated large loads (low cycle fatigue) depends on the local response. Even in that case the stresses obtained in a linear analysis may not be particularly relevant as a measure of strength. A better measure of strength is the maximum plastic strain, as calculated in a elasto-plastic analysis. Nor does a linear elastic stress analysis reflect the fact that the local yielding caused during overloads influences the subsequent growth of a crack under relatively small amplitude (high cycle fatigue conditions). The possible local yielding may blunt the 1

A ductile collapse of s tructures with in-plane compressive or shear stresses are associated with loss of stiffness due to out-of-plane displacements, and an accurate representation may require incorporation of plasticity as well as large deflection effects. 12.78


crack tips and/or change the residual stress pattern. Regarding the possible change of a residual stress field, it should be noted that current (high cycle) fatigue design rules have calibrated the strength data according to the assumption that there will be a tensile yield zone adjacent to the weld, see e.g. Gurney and Maddox (1972).

Large deflection, elasto-plastic behaviour of plate-strip subjected to lateral loading.

The purpose of the present example is to investigate the structural behaviour of a plate subjected to local lateral load such as a wheel load on a helideck, or on the deck of a ro-ro carrier etc. The actual design criterion may be stated in terms of an ultimate or serviceability limit state condition. A typical stiffened panel is shown in Fig. 12.41. If an ultimate strength criterion is used for the stiffeners and girders the corresponding strength can be relatively accurately estimated by a consideration of plastic-hinge theory. If the design criterion is formulated in terms of a serviceability requirement allowing no yielding; the structural analysis may be properly accomplished according to linear elastic beam theory. Therefore, the attention here is focused on the load-carrying behaviour of the plating. A plate-strip of unit width with material and geometric properties as shown in Fig. 12.42 is considered. The material is assumed to be mild steel with no strain hardening. The geometries correspond to slenderness ratios, b / t of 40 and 60.

Fig. 12.41 Stiffened panel with categorization of regions according to structural behaviour. The plate-strip is considered to be simply supported on the stiffeners. Two configurations of horizontal restraint are analyzed, namely complete constraint against movement and an elastic restraint corresponding to the membrane flexibility of the adjacent spans, see Fig. 12.42. The lateral load is uniform with intensity q over the central 200 mm of the loaded span. The extreme load intensity q = 2.0 N / mm 2 may refer to an accidental, crash landing of a helicopter of such a magnitude that the wheel legs collapse. q = 1.2 N / mm may be a more normal load during the landing of a helicopter. The plate-strips are modeled by beam elements of equal length, b/2, with due consideration of the symmetry of the problem. The thickness of the plate is divided by the layers.

12.79


Fig.12.42 Geometrical and material data for plate-strip. Typical load-deflection curves obtained by a finite element large deflection, elasto plastic analysis are displayed in Fig. 12.43. They are compared with solutions obtained according to small deflection theory combined with linear elastic and rigid plastic representation of the material behaviour. Due to the development of membrane action in a real plate-strip undergoing large deflections, a prediction by a small displacement theory will underestimate the stiffness and strength significantly. According to this theory the first yielding occurs at q = 0.35 N / mm 2 , and plastic collapse at q = 0.53 N / mm 2 . The more precise structural representation predicts initial yielding

to take place at about q = 0.54 N / mm 2 (for b / t

= 40;

σ Y = 320 N / mm 2 ) and no

danger of collapse even at q = 2.0 N / mm 2 .

Fig. 12.43 Load-deflection curves for a plate-strip of thickness t = 10 mm; b = 400 mm , and yield stress σ Y = 320 N /mm 2 .

12.80


Table 12.4 summarizes characteristic response values of various plate-strip configurations as obtained by a large deflection, elasto-plastic analysis. It is noted that ultimate capacity is not explicitly given in Table 12.4, because ductile collapse mechanisms are not developed for any configuration at a load intensity less than 2.0 N / mm2 . However, for the case b / t = 60 and σ Y = 240 N / mm 2 yielding over the complete cross-section is developed at q = 1.8 N / mm 2 . By sufficient additional deflection a stable equilibrium can be maintained even at q = 2.0 N / mm 2 . The maximum strain is reported because it may be a measure for ultimate capacity as defined by fracture. (Another point here is that for the two cases exhibiting the maximum strains strain hardening will occur and limit both maximum deflections and strains.) Maximum total and permanent deflection are given because they are relevant as an SLS criterion. It is observed that a load implying a plastic collapse mechanism according to small deflections theory barely seems to lead to yielding and permanent deflections when the true finite deflections are taken into account. It is particularly noted that the membrane boundary condition is an important parameter as regards total deflection. However, the permanent deflection appears to be less dependent upon the stiffness of the membrane spring used in this case study (see Fig. 12.42).

Table 12.4 Responses of plate-strips of different configurations subjected load, see Fig. 12.42. Case Response quantities Maximum Permanent Membrane central central Geometric boundary Material Load deflection, deflection parameter conditions; property σ Y q w p wmax b/t spring

stiffness, k: 40 240 1.2 ∞ 40 240 2.0 ∞ 40 320 1.2 ∞ 40 320 2.0 ∞ 40 E/40 320 1.2 40 E/40 320 2.0 60 240 1.2 ∞ 60 240 2.0 ∞ 1) The numerical solution procedure is not able case

wmax / t

w p / t

to a lateral

Maximum strain, ε max ε max / ε Y

0.95 0.64 5.7 1.36 1.14 15.6 0.81 0.35 3.1 1.14 0.76 5.5 1.10 0.41 3.8 1.49 0.79 6.1 1.48 0.90 6.6 2.20 -1) 17.5 to follow the unloading path for this

It is generally expected that the plate-strip analysis yields conservative results. The inherent conservativism should be explored by a proper two-dimensional analysis before it is applied in design. The above study has the following design implications. If sufficient membrane restraints can be mobilized by the surrounding structure, the ultimate ductile strength of a laterally loaded plate can be very high. Therefore, the total or permanent deflections may be governing in the dimensioning. A reasonable serviceability requirement on maximum allowable permanent deflection could be of the same order 12.81


of magnitude as the deflections caused during fabrication processes. If the structural component can be subjected to in-plane compressive loads after the lateral load, the possible reduction strength due to the permanent deflections must be observed.

Stiffened plates

The stiffened panel is a fundamental structural component in ship hulls. Due to the simplicity in fabrication and excellent strength to weight ratio, it is also widely used in civil engineering, bridge, aerospace, offshore and other engineering fields design of ship structures. Fig. 12.44 shows a typical stiffened plate structure.

b

Figure 12.44 A stiffened panel The main purpose of the plates is to transfer the hydrostatic loads (the difference between external and internal pressure) to the stiffeners, which again, through beam action, transfer the loads to the transverse girders. These are parts of the transverse frames of the hull girder. From the vertical girders the loads are transferred to the heavy longitudinal girders or as membrane stresses in the side of the hull. The stiffened panel in ships is generally subjected to combined in-plane and lateral pressure loads. In plane loads include biaxial compression/tension, biaxial in-plane bending and edge shear, which are mainly induced by overall hull girder bending and / or torsion of the vessel. Lateral pressure loads are due to water pressure and/or cargo. These loads components are not always applied simultaneously, but more than one normally exist and interact . Theoretically, the primary modes of overall failure for a stiffened panel can be categorized into the following five groups, namely:  Excessive plasticity of the plate after overall (grillage) buckling occurs in the elastic regime 12.82


 

 

Yielding along plate-stiffener intersection after local buckling or collapse of plating between stiffeners occurs Column or beam-column type collapse of the plate-stiffener combination as representative of the stiffened panel after local buckling of collapse of plating between stiffeners occurs Local buckling of stiffener web after local buckling or collapse of plating between stiffeners occur Tripping of stiffener after local buckling or collapse of plating between stiffeners occur

For a stiffened panel, the plate flange of a panel stiffener is usually not fully effective because plate buckling results in a non-uniform stress distribut ion (see Fig. 12.45). The rationale of effective width of buckled plating is such that the stress distribution of the stiffener and its associated effective plating can be regarded as uniform and equal to the peak stress σmax (Fig. 12.45). It would be incorrect to use the panel overall average stress, which is given as σmean. The effects of plate buckling are usually taken into account in design codes by means of the effective width approach. As illustrated in Fig. 12.45 the regions close to the stiffener are considered to be effective in carrying load while the regions remote from the stiffeners are considered to be fully ineffective in resisting compression. The ultimate load is generally considered to be reached when the maximum stress reaches the yield stress. Many different effective width equations have been proposed. Faulkner simple and useful expression is: be 2 1 b

=

β

−

β 2

, where β is the plate slenderness parameter and is defined as: β =

b

σ Y

t p

E

.

Real stress distribution max

Effective stress distribution b be mean

b

Fig. 12.44 Stress distribution in a plate in post-buckling regime

12.83


12.7 Analysis of accidental load effects 12.7.1 General

As mentioned above, the ALS check is a survival check of the structural system which is damaged due to accidental actions or abnormal strength. Accidental actions are caused by human errors or technical faults, and include fires and explosions, ship impacts, dropped objects, unintended distribution of variable deck loads and ballast, change of intended pressure difference. Most notable in this connection is of course accidental loads such as ship impacts, fires and explosions which should not occur, but do so because of operational errors and omissions. The accidental actions and abnormal conditions of structural strength are supposed to be determined by risk analysis, see e.g.Vinnem (1999), by accounting for relevant factors that affect the accidental loads. In particular, risk reduction can be achieved by reducing the probability of initiating event; leakage and ignition (that can cause fire or explosion), ship impact, etc. or by reducing the consequences of hazards. Passive or active measures can be used to control the magnitude of the accidental event and, thereby, its consequences. For instance, the fire action is limited by sprinkler/inert gas system or by fire walls. Fenders can be used to reduce the damage due to collisions. ALS checks apply to all relevant failure modes as indicated in Table 12.5. Account of accidental loads in conjunction with the design of the structure, equipment as well as safety systems is a crucial safety measure, to prevent accidents to escalate. Typical situations where direct design may affect the layout and scantlings are indicted in Table 12.6. Table 12.5 Safety criteria

Limit states Ultimate (ULS)

Fatigue (FLS)

Remarks Different for bottom – supported, buoyant, ---. Component design check Component design check depending on residual system strength after fatigue failure System design check

- Ultimate capacity1) of damaged structure (due to fabrication defects or accidental loads) or operational error Capacity to resist “rigid body” instability or total structural failure

Accidental collapse (ALS)

1)

Description - Overall “rigid body” stability - Ultimate strength of structure, mooring or possible foundation - Failure of joint

12.84


Table 12.6 Examples on accidental actions for relevant failure modes of platforms Structural concept

Failure mode

Fixed platform

Structural failure Structural failure Instability

Floating platform

Tension leg platform

Relevant accidental action

All All • Collision, dropped object, unintended pressure …, unintended balla st that initiate flooding Mooring system •Collision on platform Strength • (Abnormal strength) Structural failure All Mooring - slack • Accidental actions that initiate flooding system - strengt • Collision on platform • Dropped object on tether • (Abnormal strength)

Table 12.7 Design implications of accidental loads

Direct design of - structure (to avoid progressive structural failure or flooding) - equipment & protective barrier (to avoid damage and escalation of accident) Load

Fire Explosion Ship impact

Dropped object

Structure Columns /deck (if not protected) Topside (if not protected) Waterline structure (…. subdivision) (if not protected) Deck Buoyancy elements

Equipment

Passive protection system

Exposed equip. (if not protected) Exposed equip. (if not protected) Possibly exposed risers, … (if not protected) Equipment on deck, risers and subsea (if not protected)

Fire barriers Blast / Fire barriers Possible fender systems

Impact protection

Different subsystems, like:  loads-carrying structure & mooring system  process equipment  evacuation and escape system are designed according to criteria given for the particular subsystems. For instance, safety criteria for structural design are given in terms of ULS, FLS and ALS with specific target levels, and, hence, implies a certain residual risk level. ALS is carried out by checking the system strength after e.g. the effect of accidental actions with annual exceedance probability of 10-4 (as determined by risk analysis). A complete identification of a “Design Accidental Event” should also include an estimate of the probability of occurrence. For each physical phenomenon (fire, explosions, collisions, ..) there is normally a continuous spectrum of accidental events. 12.85


A finite number of events has to be selected by judgement. These events represent different action intensities at different probabilities. The characteristic accidental load on different components of a given installation, could be determined as follows:  establish exceedance diagram for the load on each component -4  allocate a certain portion of the reference exceedance probability (10 ) to each component  determine the characteristic load for each component from the relevant load exceedance diagram and reference probability. Iin view of the uncertainties associated with the probabilistic analysis, more pragmatic approach would normally suffice. Yet, significant analysis efforts are involved in identifying the relevant design scenarios for the different types of accidental loads. For each design accident scenario the damage imposed on the offshore installation needs to be estimated, followed by an analysis of the residual ultimate strength of the damaged structure in order to demonstrate survival of the installation. To estimate damage , i.e. permanent deformation, rupture etc of parts of the structure, nonlinear material and geometrical structural behaviour need to be accounted for. While in general nonlinear finite element methods need to be applied, simplified methods, e.g. based on plastic mechanisms, are developed and calibrated using more refined methods, to limit the computational effort required. In the following sections the estimation of damage due to accidental loads will be exemplified. 12.7.2 Fires and explosions

The dominant fire and explosion events are associated with hydrocarbon leakage from flanges, valves, equipment seals, nozzles etc. Commonly the effect of 40 - 60 scenarios need to be analyzed. This means that location and magnitude e.g. of relevant hydrocarbon leaks, likelihood of ignition, as well as combustion and temperature development(in a fire) and pressure-time development (for an explosion) needs to be estimated, followed by a structural assessment of the potential damage. The fire thermal flux may be calculated on the basis of the type of hydrocarbons, release rate, combustion, time and location of ignition, ventilation and structural geometry, using simplified conservative semi-empirical formulae or analytical/numerical models of the combustion process. The heat flux may be determined by empirical, phenomenological or numerical method (SCI,1993;BEFETS,1998). Typical thermal loading in hydrocarbon fire scenarios may be 200- 300 kW/m2 for a 15 min – 2 hours period. The structural effect is primarily due to the reduced strength with increasing temperature. In case of explosion scenarios the analysis of leaks is followed by a gas dispersion and possible formation of gas clouds, ignition, combustion and development of overpressure. Tools such as FLACS, PROEXP, or AutoReGas are available for this effort. Typical overpressures for topsides of North Sea platforms is 0.2-0.6 barg, with a duration of 0.1-0.5 s. The damage due to explosion should be determined with due account of the nonlinear and dynamic character of the action effects. Simple, conservative single degree of 12.86


freedom models may be applied. In particular cases where simplified methods have not been calibrated, nonlinear time domain analyses based on numerical methods like the finite element method should be applied. A recent overview of such methods may be found in Czujko (2001). For instance the behaviour of the topside structure of the 6-legged North Sea jacket shown in Fig. 12.46 under blast loading, has been studied. Fig. 12.47 shows the failure mode of the stiffened lower deck. The analysis was carried out by LS-DYNA using strain-based rupture criteria (Czujko, 2001). As indicated in this figure the final failure is a rupture. Fire and explosion events that result from the same scenario of released combustibles and ignition should be assumed to occur at the same time, i.e. to be fully dependent. The fire and blast analyses should be performed by taking into account the effects of one on the other. The damage done to the fire protection by an explosion preceding the fire, should be considered.

Fig. 12.46 Layout of topside of 6-legged North Sea jacket

Fig. 12.47 Failure mode of lower deck in topside structure in Fig. 12.46.

12.87


12.7.3 Ship impacts

Ship impacts on fixed platforms could cuase reduction of structural strength and possible progressive structural failure. For buoyant structures the impact damage can lead to flooding and, hence, loss of buoyancy. The measure of damage in this connection is the maximum indentation implying loss of watertightness. However, in case of large damage, reduction of structural strength is also of concern for floating structures. Ship collision loads are characterised by a kinetic energy, described by the mass of the ship, including hydrodynamic added mass and the speed of the ship at the instant of impact. If the collision is non-central, i.e. the contact force does not go through the centre of gravity of the platform (installation) and the ship, a part of the kinetic energy may remain as kinetic energy after the impact. The remainder of the kinetic energy has to be dissipated as strain energy in the installation and, possibly, in the vessel. Generally this involves large plastic strains and significant structural damage to either the installation or the ship or both. The most probable impact locations and impact geometry should be established based on the dimensions and geometry of the offshore structure and vessel, and should account for tidal changes, operational sea-state and motions of the vessel and structure which has free modes of behaviour. Impact scenarios should be established representing bow, stern and side impacts on the structure as appropriate The collision problem comprises both internal mechanics related large, inelastic deformations at the point of contact as well as global hull bending of struck vessel and interaction with the surrounding fluid (added mass, viscous forces etc.). A fully integrated analysis is fairly demanding. It is, therefore, often found convenient to split the problem into two uncoupled analyses, namely, the external collision mechanics dealing with global inertia forces and hydrodynamic effects, and internal mechanics dealing with the energy dissipation and distribution of damage in the two structures. Only the latter issue is pursued herein. This involves estimating how the energy is shared among the installation and the ship. The structural response of the ramming ship and installation can formally be represented as load-deformation relation ships as illustrated in Fig. 12.48. The strain energy dissipated by the ship and installation equals the total area under the load-deformation curves. The total energy dissipation may be expressed by:

Es

w s, max

= E s,s + E = ∫0 s, i

R s dw s

w i, max

+ ∫0

R i dw i

where Ri and Rs are the resistance of installation and ship, respectively; and dwi and dws are the deformation of installation and ship. As the load level is not known a priori an incremental procedure is generally needed. It is customary to establish the load-deformation relationships for the ship and the installation independently of each other assuming the other object infinitely rigid. Fig. 12.48b shows approximate resistance-indentation for ships. This approach may imply severe limitations, because both structures will inevitably dissipate some energy regardless of their relative strength. Care should therefore be exercised that the load-deformation curves calculated are representative for the true, interactive nature of the contact between the two structures. 12.88

Chapter 12 - Nonlinear Analysis 50

R s

R i

dws

Ship

FPSO

10

Stern end Stern corner D = 10 m = 1.5 m D

) N M ( e30 c r o f t c20 a p m I

Es,i

Es,s

40

Broad side D = 10 m = 1.5 m

dwi

D

D

Bow

0 0

1

2 Indentation (m)

3

4

a) Impact force-intendation

b) Impact force- intendation for supply vessel impact on rigid cylinrical column Fig. 12.48 Energy absorption based on force-intendation relationship

Based on the relative energy absorption capabilities of the installation and ship, the design of the installation different design principles may be distinguished; namely: strength design; ductility design , or shared-energy design. As indicated in Fig.12.49 the distribution depends upon the relative strength of the two structures. Strength design implies that the installation is strong enough to resist the collision force with minor deformation, so that the ship is forced to deform and dissipate the major part of the energy. Ductility design implies that the installation undergoes large, plastic deformations and dissipates the major part of the collision energy. Shared energy design implies that both the installation and ship contribute significantly to the energy dissipation.

n o i t a p i s s i

d y g r e n E

Ductile design

Shared-energy design

Strength design

ship installation

Relative strength - installation/ship

Fig 12.49 Ship impact design principle based on reløative energy sharing between ship and installation From the calculation point of view strength design or ductility design is favourable. In this case the response of the «soft» structure can be calculated on the basis of simple considerations of the geometry of the «rigid» structure. For instance Fig. 12.48b can be used when the ship is soft while the platform is rigid to carry out a strength design of the platform in the case of supply vessel impact. In shared energy design both the magnitude and distribution of the collision force depends upon the deformation of both structures. The analysis has to be carried out incrementally on the basis of the current deformation field, contact area and force distribution over the contact area. It is the current weaker structure that is forced to deform most, whereas the damage of the 12.89


other may remain virtually unchanged during an incremental step. The relative strength of the two structures may vary both over the contact area as well as over time. Recent advances in computers and algorithms have made nonlinear finite element analysis (NLFEM) a viable tool for assessing collisions. There are generally two methodologies available: implicit analysis and explicit analysis. Implicit methodologies require solution of equation systems. This places demands on the equation solver and the computer capacity especially in terms of memory resources. Explicit systems do not require equation solving. Equilibrium is solved at element level. However, to maintain stability, very small time steps are needed. . Explicit methodologies based computer codes include ABAQUS/Explicit, DYTRAN, LSDYNA, PAM-CRASH and RADIOSS, and implicit methodologies based codes include ABAQUS/Standard, ANSYS, MARC and NASTRAN.

The internal accident mechanics involve yielding, crushing, tearing or fracture. Any non-linear FEM mesh for simulating the internal accident mechanics needs to be fine enough so as to capture such highly non-linear characteristics. It is found that a particularly fine mesh is required in order to obtain accurate results for components deforming by axial crushing. Higher order elements generally provides better accuracy and allow a less finer mesh, but they require more computational effort. The importance of mesh fineness or element types has been studied by many investigators (e.g. Amdal & Kavlie 1992, ). It is observed that a very large number of elements is required in order to obtain accurate results for components deformed by axial crushing forces. Accounting for realistic size and boundary conditions of FE models is also crucial. Analytical formulae derived for evaluating structural damage characteristics (e.g., failure patterns) may be used to determine relevant mesh size. For instance, it is recommended that more than eight (rectangular plate-shell type) finite elements are necessary to capture the structural crushing pattern within a half length of one structural fold, see Figure 12.50. Available analytical formulations for predicting the length of structural fold are suggested e.g. by Amdahl & Kavlie (1992). It is cautioned, however, that these formulae were derived for different crushing patterns to different structural geometries.

H

H H

H

Figure 12.50: A thin-walled structure crushed under predominantly axial compressive loads and cut at its mid-section (Paik & Thayamballi 2002) An improved modeling technique for fracture behavior may be to introduce a double set of nodes such that the elements are allowed to separate once the critical stress is attained (Amdahl & Stornes 2001). 12.90


The critical strain for fracture depends heavily on the stress-strain measure as well as the mesh size. Various options exist for the stress-strain relationship. Most often engineering stress-strain relationship or true stress-strain relationship is used. The true stress-strain relationship can model the physical process more accurately than the engineering stress-strain relationship, but it is more complex as the change of the element volume needs to be involved. Experience obtained thus far indicates that the difference between these two approaches can be neglected up to ultimate stress as long as the load–displacement relation and the associated energy dissipation are concerned. The dependency of fracture strain on element size has been studied e.g. by Simonsen & Lauridsen (2000). By comparing tensile tests with numerical simulations they concluded that a fairly small element mesh is required to capture the features of the tensile test. An empirical relationship between fracture strain and element size has been derived by Lehmann et al (2001), based on thickness measurements of structural elements from actual collisions. A major challenge in NLFEM analysis is prediction of ductile crack initiation and propagation. This problem is not yet solved. Crack initiation and propagation should be based on fracture mechanics analysis, using the J-integral or Crack Tip Opening Displacement method rather than simple strain considerations. A difficulty in this connection is that the strain depends upon the element mesh. The simplest approach to the problem is to remove elements once the critical strain is attained. This is fairly easily done in an explicit code because there is no need to assemble and invert the effective system stiffness matrix. However, deleting elements disregards the fact the large stresses can be maintained parallel to the cracks. An improved modelling is to introduce a double set of nodes such that the elements are allowed to separate once the critical stress is attained. A drawback with a double set of nodes is that the potential location of cracks needs to be defined prior to analysis. Paik et al.(2003) compared FE analyses with test results obtained by The Association of Structural Improvement of Shipbuilding Industry of Japan (ASIS 1993). One of the collision test models is a double side structure model made of mild steel. It was impacted from outer side shell by a 82.32 kN weight (striking bow) freely fallen from a height of 4.8m. The weight struck the double hull model at a speed of 9.7 m/s. In FE modelling, the element size should be fine enough so that deformation patterns be properly captured in the analysis. It is desirable that the shape of the element is rectangle and the aspect ratio of an element is near 1.0. While the deformation patterns of steel plates under axial compression at the ultimate limit state have a sinusoidal shape, ship collisions and grounding cause more complex deformation patterns involving folding and tearing as well as localized yielding. To investigate the effect of mesh fineness, two meshes are considered: one is a coarse mesh usually applied for an ultimate strength analysis and the other is a fine mesh more suitable for a collision analysis. For an ultimate strength analysis, five elements between a stiffener spacing may be enough to capture the collapse pattern of the plating between stiffeners, i.e., with one element size of 80mm for the collision test model considered. However, for a collision analysis involving structural folding, a finer mesh is required. As previously noted in Section 3.5.2, at least eight elements are needed within a half length of one structural fold, see Figure 11. In the present benchmark study, the fold length was estimated by 12.91


H = 0.983 b 2 / 3 t 1 / 3 where b = width of plating between stringers, t = plate thickness, H = a half fold length. With b=2000mm and t=7mm for the test model, a half length of one fold is H = 298.5 mm. Therefore, it is recommended that one element size must be smaller than 298.5/8=37.3 mm so that at least 11 elements are necessary between a stiffener spacing (i.e., 400mm). In the present benchmark study, 13 elements were used for modelling of plating between stiffeners, which corresponds to the element size of 400/13=30.77mm. Fracture strain as obtained by the tensile coupon test was used for both types of element sizes which are relatively large. The strain rate sensitivity effect on the material yield stress of mild steel was accounted for using the Cowper-Symonds formula (with C = 40.4 and q = 5 ) , while the effect of strain rate on fracture strain was not considered (since LS-DYNA does not deal with it). Two types of material models are considered: model I for accounting for the strainhardening effect but neglecting the neck effect and model II accounting for both strain-hardening and necking effects. Figure 12.51 shows the two types of finite element modelling for the collision test model, namely a coarse mesh for the ultimate strength analysis and a fine mesh for the collision analysis. Figure 12.52 shows the deformation patterns obtained from the two types of finite element sizes.

Figure 12.51: Two types of finite element models with meshes: coarse mesh for the ultimate strength analysis (upper ) and fine mesh for the collision analysis (lower )

12.92


(a) Coarse mesh, penetration of 500mm

(b) Fine mesh, penetration of 500mm

(c) Test, penetration of 800 mm (ASIS 1993) Figure 52: Deformation patterns obtained by the mesh type with material model II (a) coarse for the ultimate strength analysis, (b) fine for the collision analysis and (c) by the test Figure 12.53 compares the force-penetration curves (left) and the absorbed energy penetration curves (right). It is evident that the fine mesh for the collision analysis captures the folding mechanism more properly than the coarse mesh for the ultimate strength analysis. In this example, the difference between material models I and II is negligible because rupture associated with necking is not a dominant failure mode. 1.6 4

① Experiment ② Coarse mesh with Material model Ⅰ ③ Coarse mesh with Material model Ⅱ ④ Fine mesh with Material model Ⅰ ⑤ Fine mesh with Material model Ⅱ

3

) N M ( 2 e c r o F

④

1.2

①

⑤

0.4

0

0 200

400

Penetration (mm)

600

③ ②

) J M ( y 0.8 g r e n E

1

0

① Experiment ② Coarse mesh with Material modelⅠ ③ Coarse mesh with Material model Ⅱ ④ Fine mesh with Material modelⅠ ⑤ Fine mesh with Material model Ⅱ

②③

⑤ ④ ①

0

200

400

600

Penetration (mm)

Figure 1253: Collision force versus penetration (left ) and absorbed energy versus penetration (right ) 12.93


Appendix A Solution of the differential equation of a beam with a axial load A.1. Differential equation

The differential equation for a beam with axial force is (Eq. 2.105)

∂ 2 ⎛ ∂ 2 w ⎞ ∂2w ⎜ EI ⎟ + P 2 = q( x) ∂ x 2 ⎜⎝ ∂ x 2 ⎠⎟ ∂ x

(A.1)

For constant EI, this equation may be written as

∂ 4 w 2 ∂ 2 w q( x) 2 P + k 2 = ; k = EI ∂ x 4 ∂ x EI

(A.2)

The general solution of Eq. (A.2) is composed of a homogenous and particular solution. The homogenous solution is obtained by solving with respect to w(x) for q(x) ≡ 0. x wh ( x) = C 1 sin 2 β 

x

+ C 2 cos 2 β + C 3 x + C 4 

wh ( x) = C 1 sinh (kx ) + C 2 cosh (kx) + C 3 x + C 4

for compressive P (A.3) for tensile P

(A.4)

where β =



P

2 EI

(A.5)

and Ci are constants to be determined by boundary conditions. The particular solution satisfies Eq. (A.1) and would, hence, depend upon q(x). For example, with q ( x) = α 1

+ α 2 x

the particular solution becomes w p ( x) =

1 ⎛ 1

1 ⎞ ⎜ α 1 x 2 + α 2 x 3 ⎟ + Ax + B 6 P ⎝ 2 ⎠

(A.6)

The constants A and B can be chosen arbitrarily, but will be influenced by the choice of C3 and C4 in the homogenous solution, since the boundary conditions should be satisfied by the total solution: w( x) = wh ( x) + w p ( x)

(A.7) 12.94


Sofar, the exact solution for a problem with distributed load, q(x) has been established. Now the solution for a beam with end moments and shear forces will be established, see Fig. A.1. This case will serve as basis for obtaining the stiffness relationship for a beam with axial force

Figure A.1 Beam with axial force, with no lateral load. The solution for this problem is obtained in another manner than by solving Eq. (A.2) directly. The moment at location C may be expressed by moment equilibrium with respect to C of the left hand part of the beam: M ( x) = − M A

+ P( w − w A ) + Qx

(A.8)

The moment-curvature relationship is given by Eq. (2.77) as

∂2w M = − EI 2 ∂ x

(A.9)

By combining Eqs (A.4, A.5)

∂2w − EI 2 = − M A + P( w − w A ) + Qx ∂ x

(A.10)

∂2w P 1 ( Pw + M A − Qx) w + = ∂ x 2 EI EI A

(A.11)

or

Eq. (A.11) is a 2nd order differential equation as opposed to the 4th order one in Eq. (A.2). The homogeneous solution of Eq. (A.7) for a positive P is: wh

= C 1 sin kx + C 2 cos kx

; k 2

=

P EI

(A.12)

The solution for a negative P is similar, with hyperbolic functions instead of trigonometric ones. The total solution (for P>0) may be written as: 12.95


w( x) = C 1 sin kx + C 2 cos kx +

A.2

1 P

( Pw A + M A

− Qx)

(A.13)

Stiffness matrix with stability function

The exact solution for a beam with an axial force (and no lateral load) is

⎛ x ⎞ ⎛ x ⎞ 1 w( x) = C 1 sin ⎜ 2 β ⎟ + C 2 cos⎜ 2 β ⎟ + ( Pw A + M A ⎝  ⎠ ⎝  ⎠ P w ( x) = C 1 ,

where β =

2 β 



⎛ x ⎞ cos⎜ 2 β ⎟ − C 2 ⎝  ⎠

2 β 

⎛ x ⎞ Q sin ⎜ 2 β ⎟ + ⎝  ⎠ P

− Qx)

(A.14)

P

2 EI

This solution applies for an arbitrary beam as shown in Fig. A.2, with transverse displacement and rotations of the ends (wA, θA) and (wB, θB), end moments MA, MB while the shear force is constant and equal to Q (due to no lateral load, q(x)). Based on the solution (A.14), kinematic boundary conditions and equilibrium equations, a relationship between forces MA, M B and Q and displacements wA, θA, w B and θB can be established as shown in as follows.

Figure A.2: Beam with axial force. The boundary conditions are w(0) = C2 + w A w '(0) = C1

2β 

−

+

MA

Q

= −θA

P

P

= wA

w( ) = C1 sin 2β + C2 cos 2β + w A + w '( ) = C1

2β 

cos 2β − C 2

2β 

MA

sin 2β −

P Q P

−

Q



P

= wB

= −θB

The first two equations result in: 12.96


C 1

M Q ⎞  = ⎛ ⎜ − θ A ⎟ ; C 2 = − A P ⎝ P ⎠ 2 β

The last two equations then become: M A

(1 − c) −

P M A

2 β s −

Q ⎛ s ⎞ ⎜⎜1 − ⎟⎟ + w A P ⎝ 2 β ⎠

Q

P P where s = sin 2 β , c

(1 − c) − θ A c + θ A

− w B − θ A 

s

2 β

=0

=0

= cos 2 β

Solving these equations with respect to MA and Q yields:

⎡ 2 β (1 − c) ( w − w ) + ( s − 2 β c)θ + (2 β − s)θ ⎤ B A A B ⎥ 4 β 1 − c − β s ⎢⎣  ⎦ 1

P

M A

=

Q=

⎡ 2β s (w − w ) + (θ + θ )(1 − c)⎤ B A A B ⎥⎦ 2(1 − c −β s) ⎢⎣  P

(A.15a)

(A.15b)

Moment equilibrium requires MA

+ M B − P(w B − w A ) − Q = 0

which yields MB

= −M A + P(w B − w A ) + Q P 1 ⎡ 2β (1 − c)(w − w ) + (2β − s)θ + (s − 2β = A B A 4β 1 − c −β s ⎢⎣ 

⎤ ⎦

c)θ B ⎥

The expressions (A.15a-c) may be simplified by using the so-called defined by:

φ1 = β

cot gβ

1 3

β 1 − φ1

1

3

4

4

φ2 =

;

β=



P

2

EI

(A.15c)

φi-functions,

2

φ3 = φ1 +

φ2

1

3

2

2

φ4 = − φ1 +

(A.16)

φ2

φ5 = φ1φ2 By using these functions Eq. (A.15a-c) is transformed into: 12.97


MA

=

6EI

MB

=

6EI

Q=

12EI







2

2

3

φ2 (w B − w A ) +

4EI

φ2 (w B − w A ) +

4EI

φ5 (w B − w A ) +





6EI 

2

φ3θA +

2EI

φ 4 θA +

4EI

φ 2 θA −





6EI 

2

φ4 θB φ3 θB

(A.17)

φ2 θB

By means of Eq. (A.17) the stiffness matrix for an element with 4 d.o.f. can be established. The relations (A.17) may then be used to establish the stiffness relationship S = kv, with Si and vi defined as shown in Fig. A.3b, by observing that:

= −Q, S 3 = M A , S 5 = Q, S 6 = M B v 2 = w A , v3 = θ A , v5 = w B , v6 = θ B

S 2

(A.18)

By combining these relations with Eq. (A.18) the resulting equations are written on matrix form as:

⎡ S 2 ⎤ ⎡ 6φ 5 − 3φ 2 − 6φ 5 − 3φ 2 ⎤ ⎢ S ⎥ ⎢− 3φ 2 2φ 3φ ⎥ 2  φ 4 2 EI 3 2 3 2 ⎢ ⎥= ⎢ ⎥ ⎢ S 5 ⎥  3 ⎢ − 6φ 5 3φ 2 6φ 5 3φ 2 ⎥ ⎢ ⎥ ⎢ ⎥ 2 2 ⎣ S 6 ⎦ ⎣− 3φ 2  φ 4 3φ 2 2 φ 3 ⎦

⎡v 2 ⎤ ⎢v ⎥ ⎢ 3⎥ ⎢ v5 ⎥ ⎢ ⎥ ⎣v 6 ⎦

(A.19)

or S,

, = k NL v,

a) 6 d.o.f. model

b) 4 d.o.f. model

Figure A.3 Beam element with axial force.

12.98


Appendix B General formulations for geometrically nonlinear behaviour B.1 Continuum mechanics General

The purpose of the present appendix is to briefly define some fundamental concepts of the mechanics of continua, restricted to three-dimensional Euclidean spaces. Coordinates are means of identifying or labeling material particles. The simultaneous position of the set of particles comprising the body is called the configuration of the body. The motion of the body is a continuous sequence of configurations in space and time.

Figure B.1 Deformation of a body. Two configurations of a body in motion are shown in Fig. B.1. The initial configuration is denoted C0. P0 indicates the initial location of a material particle. The current configuration is denoted C, and P marks the present location of the material particle. The original configuration is described by the coordinates xi. To describe the motion og the body two choices of coordinates are preferred: 

the same rectangular Cartesian coordinates are used to describe both the original reference and deformed configurations (e.g. xi in Fig. B.1).



the frame of reference is distorted in such a way that the coordinates of Xi of a particle have the same numerical values xi as the reference configurations. Such coordinates that follow deformations become survilinear (called convective coordinates).

or

For the present purpose the first choice is most appropriate and is hence, used. If the reference configuration within this choice is the actual initial configuration at t0 = 0, and the independent variables are the coordinates xi and time t, then the description of motion is called material description or Lagrangian description . The alternative is spatial description whose independent variables are the present position xi occupied by the particle at time t and the present time t. The spatial description, 12.99


also called the Eulerian description focuses attention on fixed points of space instead of on a given points in the body. The description is most used in fluid mechanics. We shall limit ourselves here to material or Lagrange description. However, it is necessary to relate the coordinates xi of the initial configuration with the coordinates Xi of the displaced body. The relation between the two coordinates is given by the displacement as illustrated in Fig. B.2. In component form the relation is: X i

= xi + u i

(B.1)

where ui = ui(xi, t)

Figure B.2 Definition of displacement. Strains

Alternative expressions for strain and stress appropriate for large deformation problems are envisaged (e.g. Crisfield, 1991). They include the Almansi and Green strains. The Green (-Lagranage) strain is commonly applied. It is defined by the equations

E ij

∂u ∂u ∂u ⎞ 1 ⎛ ∂u = ⎜⎜ i + j + k k ⎟⎟ 2 ⎝ ∂ x j ∂ xi ∂ xi ∂ x j ⎠ =

1 2

(u

i ' j

(B.2)

+ u j 'i + u k 'i u k ' j ) 3

when adopting the Einstein summation convention that u k u k

= ∑ u k 2 . i =1

The Almansi strain is defined with reference to the deformed geometry. By defining ui by u1 = u, u2 = v and u3 = w, Eqs (B.2) may be written fully as:

12.100


E xx

= u, x +

E yy

= v, y +

E zz

= w, z +

1 2 1 2 1

(u,

2 x

+ v, x2 + w, x2 )

(u,

2 y

+ v, y2 + w, y2 )

(u,

2 z

+ v, z2 + w, z2 )

(B.3)

2 E xy = u, y + v, x + (u , x u , y

+ v, x +v, y + w, x w, y ) E yz = v, z + w, y + (u , y u , z + v, y + v, z + w, y w, z ) E zx = w, x + u , z + (u , z u , x + v, z + v, x + w, z w, x ) The initial terms in Eqs (B.2) are the customary engineering definitions of normal and shear strain (εx = u,x, etc.). The added terms, in parentheses, become significant if displacement gradients are not small. Green-Lagrange strains are zero for a rigid-body rotation of any magnitude. In Eqs (B.2), all displacement derivatives are computed in the original coordinate system, regardless of how large a rigid-body rotation may be supposed on the deformations. This is the “total Lagrangian” approach, in which all displacements are measured in a reference frame that is stationary rather than attached to the deforming structure. The stationary coordinates may also be called “material coordinates” and may be denoted in some papers by uppercase labels X, Y and Z (or xi). Green normal strains correspond to defining the strain of a line segment by the equation E =

1 ⎡⎛ ds * ⎞

⎤ 1 − ⎜ ⎟ ⎢ ⎥ 2 ⎣⎢⎝ ds ⎠ ⎥⎦ 2

(B.4)

where ds and ds* are respectively the initial and final lengths of the line segment. If ds ≈ ds*, Eq. (B.3) reduces to the usual small-strain approximation, E = E (ds*-ds)/ds. Stress

Consider the body in deformed state. In Fig. B.3 Δf is the force acting on the surface element Δa and n is the unit normal vector at Δa with outwards direction.

Figure B.3 Surface force on deformed body.

12.101


Assume that the ration Δf /Δa tends to a definite limit at Δa tends to zero. The force t per unit area acting on the surface in the limit condition is denoted traction

Δf d f = Δa →0 Δa da

t = lim

(B.5)

The equilibrium equation at the surface of the body gives the definition of the Eulerian (Cauchy) stress tensor , σij σ ij ni

= t j

(B.6)

In Eq. (B.6) t j and ni are, respectively, components of t and n. σij refers to the current configuration which is a natural physical concept. However, stresses must be related to strains. In the Lagrangian description of motion the strains Eij are referred to the original position of the particles. Therefore, it is convenient to define stresses to the initial area. The stress tensor should be energy conjugate with the strains Eij. Since Eij is symmetric, the stress tensor should also be symmetric. These stresses, Sij are denoted Piola-Kirchhoff stresses, and can be expressed by the Cauchy stress, σlk as follows S ji

=

∂xi ∂x j ∂Xk ∂X

σk 

/ det |

∂xi

/

∂X j

|

(B.7)



where det | ∂xi/∂X j | is the determinant of the deformation matrix {∂xi/∂X j}. For normal strains this determinant is equal to unity. Equilibrium equations

It can be shown (e.g. Crisfield, 1991) that equilibrium equations for a body with volume forces can be formulated with reference to the initial geometry as

∂u i ⎞⎤ ∂ ⎡ ⎛ ⎟⎥ + ρ F = 0 ⎢ S jk ⎜⎜ δ ik + ∂ x j ⎣ ⎝ ∂ xk ⎠⎟⎦ 0 0i

(B.8)

Virtual work

Eq. (B.8) expresses the equilibrium of a body in the direction i. The body is acted upon by surface tractions and body forces. The surface of the body under consideration may be thought of as consisting of two areas. On area S1 the surface tractions are specified, while the displacements are specified on the area S2. By giving the body from its equilibrium position a virtual displacement δui, which is kinematically consistent with the boundary conditions, the following equation may be derived:

12.102


1

∫

2 V

⎛ ⎜ ⎝

⎛ ∂u i ⎞ ⎛ ∂u j ⎞ ⎛ ∂u k ∂u k ⎞ ⎞⎟ ⎟ ⎜ ⎟ ⎟ ⎜ ⎜ ∂ x j ⎟ + δ ⎜⎝ ∂ xi ⎠⎟ + δ ⎜ ∂ xi ∂ x j ⎟ ⎟ dV ⎝ ⎠ ⎝ ⎠ ⎠

S ij ⎜ δ ⎜

(B.9)

− ∫ T i δ ui dS − ∫ ρ F i δ u i dV = 0 S1

V

The index zero has been omitted here as everything is referred to the original configuration. T i denotes the prescribed surface tractions. By introducing the Green’s strain tensor of Eq. (B.2) in Eq. (B.9) the virtual work equation may be written as

∫ S δ E dV − ∫ T δ u dS − ∫ ρ F δ u d V = 0 ij

ij

V

i

S1

i

i

i

(B.10)

V

B.2 Principle of virtual displacement on incremental form Total Lagrange formulation

In this section an incremental form of the virtual work equation will be shown. Fig. B.4 shows the spatial coordinate system xi and the body in three different configurations. C0 is the initial configuration. Cn is some deformed configuration and Cn+1 is a configuration close to Cn.

Figure B.4 Initial and deformed configuration of a body. The virtual work principle (Eq. (B.10)) in configuration Cn+1 reads

∫ S

*

ij

∫

∫

δ E ij* dV = T i * δ u i* dA + F i * δ u i* dV

v

S

(B.11)

V

In Eq. (B.11) superscript * denotes quantities in configuration Cn+1. The superindex  on the stress S ij denotes stress due to external loads. The relations between these quantities in configurations Cn+1 and Cn are S ij* = S ij + ΔS ij E ij* = E ij u k *,i

+ Δ E ij

= u k ,i + Δu k ,i

(B.12a-e)

T i * = T i + ΔT i F i * = F i + ΔF i 12.103


The Δ denotes increment when moving from Cn to Cn+1. The virtual displacment fields for configurations Cn and Cn+1 are assumed to be identical so that δ u i* = δ u i . The virtual work principle for configuration Cn is expressed as in Eq. (B.10). Introducing Eqs (B.12) into Eq. (B.11) yields an incremental form of the virtual work principle as

∫ ΔS δ E dV + ∫ S δ Δ E dV + ∫ ΔS δ Δ E dV 

ij



ij



ij

V

ij

V

ij

ij

V

⎛ ⎞ = ∫ ΔT iδ u i dS + ∫ ΔF iδ u i dV − ⎜⎜ ∫ S ij δ E ij dV − ∫ T i δ u i dS − ∫ F i δ u i dV ⎟⎟ S V S V ⎝ V ⎠

(B.13)



If the body in configuration Cn is in equilibrium, the parenthesis on the right hand side of Eq. (B.13) will vanish according to Eq. (B.10). However, due to approximations in the solution procedure, configuration Cn will generally not be in equilibrium before onset of further loading. The terms in the parenthesis will then act as an equilibrium correcting term. The variation of Green’s strain terror, Eq. (B.4), is δ E ij

= 12

δ u i , j

+ δ u j ,i + u k ,i δ u k , j + δ u k ,i u k , j )

(Note that δ E ij = E ij (u i

(B.14)

+ δ u i ) − E ij (u i ) .)

The increment of Green’s strain tensor when moving from configuration Cn to Cn+1 is

ΔEij =

1 2

( Δu

i,j

+ Δu j,i + Δuk ,i uk ,j + uk,i Δ uk ,j + Δ uk ,i Δ uk ,j )

(B.15)

The variation of the increment of Green’s strain tensor reads δ Δ E ij

= 12

δ u k , j Δu k ,i

+ δ u k ,i Δu k , j )

(B.16)

The constitutive relation is assumed to have the following form

ΔS ij = C ijk Δ E k 





(B.17)

where C ijk  is the material tensor. The material modeling will be discussed in Section 10.4. Eq. (B.17) may also be written as a vector-matrix relation if the components ΔSij (Δ Ekl ) are collected in vectors, and Cijkl in a matrix. Updated Lagrangian formulation

In the updated Lagrangian description the displacement increments are referred to the preceding configuration Cn. The large displacements and rotations are taken into account by updating the geometry. Hence, the nonlinear terms may be neglected in Eqs (B.14 – B.15). They are, however, the only terms contributing to the variation of the strain increment, see Eq. (B.17). 12.104


The discussion has revealed that the equations resulting from the virtual work principle on incremental form as expressed in Eq. (B.13) are nonlinear in Δui. Linearization is obtained by neglecting second order terms. This is done by omitting the third integral on the left hand side of Eq. (B.13) as well as by using the linearized strains

= 12 (δ u i , j + δ u j ,i ) Δ E ij = 12 ( Δ u i, j + Δu j ,i )

δ E ij

(B.18a-b)

Finally, the linearized form of the virtual work principle in incremental form can be written as

∫C

ijk

ΔEk δEijdV + 

V

= ∫ ΔTi δui dS + S

∫ S δΔE dV ij

ij

V

⎛ Δ F δ u dV − ⎜ ∫ SijδEijdV − ∫ Tiδ i i ∫V S ⎝V

ui dS

−

∫

V

⎞

Fiδui dV ⎟

⎠

(B.19) where Sij are the total stresses within the structure.

12.105


Appendix C Plasticity theory In the same manner as for the one dimensional case, elasto-plastic behaviour of metals in a multiaxial stress state can be characterized by   

An initial yield condition , i.e. the state of stress for which plastic deformation first occurs. A hardening rule which describes the modification of the yield condition due to strain hardening during plastic flow. A flow rule which allows the determination of plastic strain increments at each point in the load history.

It is assumed that the material is isotropic, which implies that the stiffness properties are independent of orientation at a point. In this section a brief review of the plasticity theory will be given. More details may be found in textbooks like Crisfield(1991), Chapter 6. C.1 The von Mises initial yield condition

The yield condition of a material defines the limit of purely elastic behaviour under any combination of stresses. Many initial yield conditions have been proposed. Experiments indicate that the von Mises yield condition best represents material behaviour of most metals. Another advantage of the von Mises criterion is its simple continuous function of stress components which makes it especially attractive to numerical analysis. The mathematical expression for the von Mises initial yield surface reads f = σ − σ Y

= 0

(C.1)

where f is the loading function and σY is the initial uniaxial yield stress of the material. The equivalent stress σ is given by σ

=

3 2

sij sij

=

3 2

σ ijσ ij

− 12 (σ kk ) 2

(C.2)

where the Einstein’s summation convention is applied i.e. n

sij s ij

n

= ∑∑ sij2 ; σ i =1 j =1

kk =

n

∑σ k =1

kk

;

σ m = 13 σ kk .

sij denotes the deviatoric, or reduced stress tensor sij

= σ ij − 13 δ ijσ kk

(C.3)

where

⎧1 δ ij = ⎨ ⎩0

for

i= j

for

i≠ j

12.106


C.2 Hardening rule

By using the von Mises yield criterion, the yield criterion in connection with hardening may be written as:

)

σ Y = H ε p

(C.4)

where the equivalent plastic strain is ε p

ε

p

=

∫ d ε

p

(C.5)

0

with the equivalent plastic strain increment d ε p =

(2 / 3)d ε pij d ε pij

(C.6)

(again using Einstein’s summation convention). It is also convenient to introduce the plastic work, W p: ε p ij

W p

=

ε p ij

∫ σ d ε = ∫ σ d ε ij

0

p ij

p

(C.7)

0

The hardening may be described by different models  

isotropic kinematic or generally anisotropic manner as illustrated for a one-dimensional and two-dimensional condition in Fig. C.1.

The choice between models especially matters if the load condition/stress is reversed or cyclic.

a) One-dimensional

b) Two-dimensional

Figure C.1 Hardening models 12.107


Experiments with metals show that a phenomenon denoted Bauschinger effect occurs. This means that the material yields a lower stress level when the loading is reversed than during the initial loading. As shown e.g. for the one-dimensional case in Fig. C.1a, this feature is not captured by the isotropic hardening. In the present section only isotropic hardening is considered. The yield condition is written as f = σ − H (ε p ) = 0

(C.8)

where σ is given by Eq. (C.2). C.3 The flow theory of plasticity

The relation ship between stress and (plastic) strain may be obtained mainly by the two major different plasticity theories, namely the deformation theory and the flow theory. Experiments show that the flow theory is the better one when treating problems with general loading paths, e.g. reversed and cyclic loading. Therefore, the flow theory of plasticity is described herein. The flow theory of plasticity yields an incremental relationship between Cauchy stresses σij and true strains εij. For small strains the same constitutive law can be used between Piola-Kirchhoff stresses and Green’ s strains. It is a fundamental assumption that the total strains εij may be decomposed into elastic components ε ije and plastic components ε ij p ε ij

= ε ije + ε ij p

(C.9)

The elastic components of strain are related to the stresses by Hooke’s law.

εeij = Cijklσ kl = 9K1 δ ijσ kk + 2G1 s ij =

1+ ν E

ν σ ij − δijσ kk E

(C.10a)

or inversely

σij = Eijklεekl = 2Gε eij +

3νK 1+ ν

δ ijε ekk

(C.10b)

where the compression modulus K, modulus of elasticity, E and shear modulus, G are such that:

3(1 − 2ν ) K = 2(1 + ν )G = E

(C.10c)

The plastic mode is assumed to be incompressible, implying ε ii p

=0 12.108


As was expected in Eq. (C.1) for initial yielding it is assumed that in general there exists a loading function f. Setting f equal to zero defines a yield surface that bounds the elastic range. The mathematical expression for the yield surface may be written f σ ij , ε ij p , κ

)= 0

(C.11)

for all σij and ε ij p components. Eq. (C.11) defines the initial as well as subsequent yield surfaces depending upon the values of the variables. κ is a hardening parameter. Different values of f define different stress states f < 0 ; elastic state of stress f = 0 ; plastic state of stress

(C.12a-c)

f > 0 inadmissib le

Taking the differential of f from a plastic state gives df =

∂ f ∂σ ij

d σ ij

+

∂ f ∂ε ij p

d ε ij p

+

∂ f d κ ∂κ

(C.13)

again by invoking Einstein’s summation convention. Then, three different loading conditions can be defined:

∂ f ∂σ ij ∂ f ∂σ ij ∂ f ∂σ ij

d σ ij

< 0 ; f = 0

during unloading

d σ ij

= 0 ; f = 0

during naytral loading

d σ ij

> 0 ; f = 0

during loading

(C.14a-c)

The final assumption is that the so-called Drucker’s postulate is valid. This implies that the yield surface is convex. Moreover, it means that the plastic strain increment vector is directed along the outward normal to the yield surface. The normality implication may be expressed as d ε ij p

= d λ

∂ f ∂σ ij

(C.15)

where dλ is a nonnegative constant. Eq. (C.15) indicates that the loading function may be taken as a plastic potential . According to Eq. (C.14b), the function value of f remains unchanged, like zero, from one plastic state to another. 12.109


Since,

)

κ = κ ε ij p

(C.16)

the function f is dependent upon two sets of variables, the stresses and plastic strains

)

f = f σ ij , ε ij p

(C.17)

Hence, Eq. (C.13) may be written as df =

∂ f ∂σ ij

d σ ij

+

∂ f ∂ε ij p

p

d ε ij

= 0

(C.18)

By substituting Eq. (C.18) into Eq. (C.15) yields

∂ f d σ ∂σ ij ij d λ = − ∂ f ∂ f ∂ε pmn ∂σ mn

(C.19)

From Eqs (C.2, C.8)

∂ f ∂σ 3 sij = = ∂σ ij ∂σ ij 2 σ

(C.20)

and

∂ f ∂ = p p ∂ε ij ∂ε ij

(σ − H (ε )) = − p

∂ ∂ε

p ij

( H (ε ))

1 ∂ H ∂ε p ∂W p = − p = − H ' σ ij p p σ ∂ε ∂W ∂ε ij

p

(C.21)

where the plastic work, W p (Eq. (C.7)) has been utilized. By noting

∂ f ∂σ = ∂σ ij ∂σ ij and

σijs ij = σijσij − 13 σ ijδ ijσ m = σ ijσ ij − 13 σ2m = 23 σ2

12.110


as well as using Eqs (C.20 – C.21) to rewrite the demoninator of Eq. (C.19), dλ may be written as:

d λ = −

∂σ ∂σ ij

d σ ij

σ ⎞ ⎛ s ⎞ ⎛ ⎜⎜ − H ' ij ⎟⎟ ⎜⎜ 3 ij ⎟⎟ σ ⎠ ⎝ 2 σ ⎠ ⎝

=

d σ H '

(C.22)

Now, the derivative of the hardening parameter H with respect to plastic strain, ε p can be obtained from uniaxial test. This is obtained from Eq. (10.94) by replacing σ and ε p with the equivalent quantities σ and ε p for multidimensional cases. By combining Eqs (C.15, C.20, C.22) the following expressions for plastic strain increments result: dε pij

=

sij

3

2 H 'σ

dσ

(C.23)

These expressions are called the Brandtl-Reuss equations for an isotropically hardening material, obeying the von Mises yield criterion. It is seen from Eq. (C.23) that the relative increments of yielding is determined by the relative magnitude of total stresses. The stress increments only relate to the scalar d σ . It may also be convenient to express the plastic strain increments by the stress increments. The differential d σ in Eq. (C.25) is then expressed as follows: d σ

=

∂σ ∂σ kl

d σ kl

=

3 2σ

s kl d σ kl

(C.24)

Then by inserting in Eq. (C.23) dε pij

= α ⋅ sijs kldσ kl

(C.25)

where α =

9 4σ 2 H ' (σ )

Besides the relation between plastic strain and (equivalent) stress increment it is also of interest in finite element formulations to have the relationship between plastic strain increments and the total strain increments, i.e. including elastic increments. By reformulating Eq. (C.19) and introducing Hooke’s law:

12.111


dλ

⎛ ∂f ∂f ⎞ ∂f ⎜ ∂ε p ∂σ ⎟ = − ∂σ dσij ⎝ mn mn ⎠ ij ∂f = − (E ijkldε ekl ) ∂σij = −

∂f ∂σij

E ijkl

(C.26)

⎛ ∂f ⎞ ⎜ dε kl − dλ ∂σ ⎟ ⎝ kl ⎠

Introducing the following relations: E ijkl s kl E ijkl

= 2Gsij

= E klij

(C.27)

and applying Eqs (C.20, C.21) gives d λ =

3Gs kl d ε kl σ ( H '+3G )

(C.28)

By Eqs (C.15, C.20) the plastic strain increment now becomes d ε ij p

=

9Gsij s kl 2σ 2 ( H '+3G )

d ε kl

(C.29)

This is the relationship between plastic and total strain. Finally, the incremental stress-strain relationship than gets the following form: d σ ij

= E ijkl d ε kle

( Hooke' s law)

= E ijkl (d ε kl − d ε pkl )

(C.30)

ep d ε kl = Dijkl

where ep ijkl

D

= E ijkl −

9G 2 s ij s kl σ 2 ( H '+3G )

(C.31)

It is noted that Dijkl is symmetric Dijkl

= Dklij

because both terms in the expression ( C.31) are symmetric.

12.112


C.4 Incremental stress-strain relationships for oriented bodies made of thin plates Beam

The following assumption are made: Strains: γ xy γ zx

= d γ xy = 0 = d γ zx = 0

; γ yz

= d γ yz = 0

(C.32)

Stresses: σ xy σ yy

= d σ xy = 0 = d σ yy = 0

; σ yz ; σ zz

= d σ yz = 0 = d σ zz = 0

(C.33)

By taking into consideration these assumptions, the following stress-strain relation can be found

D ⎛ d σ xx ⎞ ⎛ ⎜ ⎟ ⎜ xxxx ⎜ d σ yy ⎟ = ⎜ D yyxx ⎜ d σ ⎟ ⎜ D ⎝ xx ⎠ ⎝ zzxx

D xxyy D xxzz ⎞ ⎛ d ε xx ⎞

⎟⎜

⎟

D yyyy D yyzz ⎟ ⎜ d ε yy ⎟

(C.34)

⎟

D zzzz ⎠ ⎜⎝ d ε zz ⎠⎟

D zzyy

or d σ ij

ep d ε kl = Dijkl

(C.35)

or in matrix notation dσ = Ddε

(C.36)

Eq. (C.33) and Eq. (C.34) gives, when symmetry of stiffness matrix is taken into consideration d σ xx

= ( D xxxx − −

D xxyy 2 D yyyy D zzzz − D yyzz

D xxzz 2 D yyyy D zzzz − D yyzz

( D

( D

D zzzz − D zzxx D yyzz

yyxx

D yyyy − D yyxx D yyzz

zzxx

) ) d ε

xx

)

The deviatoric stress components expressed with total stresses become s xx = 23 σ xx ; s yy = s zz = − 13 σ xx

(C.36a)

(C.37)

12.113


Eq- (C.31) and Eq. (C.37) now give the following expressions for the stiffness components of Eq. (C.34) D xxxx

2 9G 2s xx

= E xxxx −

E(1 − ν )

=

−

2 4G 2 σxx

σ 2 (H '+ 3G) (1 + ν )(1 − 2ν ) σ2 (H '+ 3G) 2 Eν 2G 2 σxx − D xxyy = D yyxx = (1 + ν )(1 − 2ν ) σ2 (H '+ 3G) D xxzz = D zzxx = D xxyy D yyyy

= D zzzz =

E(1 − ν ) (1 + ν )(1 − 2ν )

(C.38)

2 G 2 σxx

−

σ2 (H '+ 3G) 2 Eν G 2 σxx − 2 D yyzz = D zzyy = (1 + ν )(1 − 2ν ) σ (H '+ 3G) By combination of Eq. (C.36) and Eq. (C.38) the stress-strain relation can be reduced to an expression of the form d σ xx

= Dd ε xx

(C.39)

Plate

As is generally done for thin plate theory, the following case will be considered. Strains: γ zx

= d γ zx = 0 ;

γ yz

= d γ yz = 0

;

σ yz = dσ yz =

(C.40)

Stresses:

σ zx = dσzx = 0 σ zz = d σzz = 0

0

(C.41)

The stress-strain relation can then be reduced to a 4x4 expression

⎛ d σ xx ⎞ ⎜ ⎟ ⎜ d σ yy ⎟ ⎜ d σ ⎟ ⎜⎜ xy ⎟⎟ ⎝ d σ zz ⎠

⎛ D xxxx ⎜ ⎜ D yyxx ⎜ D ⎜ xyxx ⎜ D zzxx ⎝

D xxyy D xxxy D xxzz ⎞ D yyyy D yyxy D xyyy D xyxy D zzyy

D zzxy

⎟ D yyzz ⎟ D xyzz ⎟ ⎟ D zzzz ⎠⎟

⎛ d ε xx ⎞ ⎜ ⎟ ⎜ d ε yy ⎟ ⎜ d ε ⎟ ⎜⎜ xy ⎟⎟ ⎝ d ε zz ⎠

(C.42)

or

dσij

= Dijkldε kl + D ijzzdε zz

(C.43)

and d σ zz

= D zzkl d ε kl + D zzzz d ε zz

(C.44)

This equation system may be contracted to the following expression 12.114


d σ ij

D D ⎞ ⎛ = ⎜⎜ Dijkl − ijzz zzkl ⎟⎟ d ε kl D zzzz ⎠ ⎝

(C.45)

It is seen that the system (10.111) is symmetric. The deviatoric stress components for this two-dimensional case becomes s xx

= 13 ( 2σxx − σ yy )

;

s yy

= 13 ( 2 σyy − σ xx )

s xy

= σ xy

;

s zz

= − 13 ( σ xx + σ yy ) =− σ m

(C.46)

In Eq. (10.108) the stiffness matrix elements are D xxxx D xxyy

= E xxxx −

2 9G 2s xx

σ 2 (H '+ 3G)

=

Eν

= D yyxx =

(1 + ν )(1 − 2ν )

−

9G 2 s2xx

σ2 (H '+ 3G)

9G 2s xx s yy

−

σ2 (H '+ 3G) Eν 9G 2 sxx σ m D xxzz = D zzxx = + (1 + ν )(1 − 2ν ) σ2 (H '+ 3G) 9G 2s xx σxy D xxxy = D xyxx = − 2 σ (H '+ 3G) D yyyy D yyzz D yyxy

=

(1 + ν )(1 − 2ν )

E(1 − ν )

E(1 − ν ) (1 + ν )(1 − 2ν)

= D zzyy =

−

9G 2s2yy

σ2 (H '+ 3G)

Eν (1 + ν )(1 − 2ν )

= D xyyy = −

+

9G 2s yyσm

σ2 (H '+ 3G)

9G 2s yy σ xy

σ 2 (H '+ 3G) 9G 2 σm σ xy D zzxy = D xyzz = 2 σ (H '+ 3G) 9G 2 σ2xy D xyxy = G − 2 σ (H '+ 3G) 2 E(1 − ν ) 9G 2σm D zzzz = − (1 + ν )(1 − 2ν ) σ2 (H '+ 3G)

(C.47)

For the equivalent stress σ the following expression comes out of Eq. (C.2) σ

=

2 σ xx

+ σ yy2 − σ xxσ yy + 3σ xy2

(C.48)

and the equivalent plastic strain increment from Eq. ( C.6 ) becomes d ε p

=

4 3

(d ε

p 2 xx

p p + d ε yy + d ε xx p d ε yy ) + 13 d γ xy p 2

2

(C.49) 12.115


References Amdahl, J. and Kavlie, D. (1992). “Experimental and numerical simulation of double hull stranding”, DNV-MIT Workshop on Mechanics of Ship Collision and Grounding, Høvik, September.

Amdahl, J. and Stornes, A. (2001). “Energy dissipation in aluminium high-speed vessels during collision and grounding”, Proceedings of ICCGS’01, Copenhagen, 203219. ASIS (1993). The Conference on “Prediction Methodology of Tanker Structural Failure & Consequential Oil Spill”, Association of Structural Improvement of Shipbuilding Industry of Japan, III.39-III.47. Bell, K. (1994). “Matrix Statics”, Tapir, Trondheim (in Norwegian). Becker, A.A. (2001). “Understanding Non-linear Finite ElementAnalysis through Illustrative Benchmarks”, NAFEMS, Glasgow. BEFETS: “Blast and Fire Engineering for Topside Systems”, Phase 2, SCI Publication No. 253, Ascot, UK, 1998. Bergan, P.G. and Syvertsen, T.G. (1978). ”Buckling of Columns and Frames”, Tapir, Trondheim (in Norwegian). Bergan, P.G., Horrigmoe, G., Kråkeland, B. & Søreide, T.H., (1978), “Solution techniques for non-linear finite element problems”, Int. J. Num. Meth. Engngn., 12 , 1677-1696. Cook, Robert D., Malkus, David S. and Plesha, Michael E. (1989) “Concepts and Applications of Finite Element Analysis” John Wiley & Sons, ISBN 0-471-84788-7. Crisfield, M.A. (1991). “Non-linear Finite Element Analysis of Solids and Structures”, Vol. 1, J. Wiley & Sons, Chichester. Cullen, The Hon. Lord, 1990. ‘The Public Inquiry into the Piper Alpha Disaster’, London: HMSO. Czujko, Jerzy (Ed.) (2001). “Design of Offshore Facilities to Resist Gas Explosion Hazard.” Engineering handbook. First edition. CorrOcean ASA, ISBN 82-996080. Den Heijer, C. & Rheinboldt, W.C. (1981), “On step-length algorithms for a class of continuation methods”, SIAM J. Num. Analysis, 18 , 925-948. Emami Azadi, M., Moan, T. and Hellan, Ø., (1997). “Nonlinear Dynamic vs. static Analysis of Jacket Systems for Ultimate Limit State Check”, proc. Int. Conf. on Advanced in Marine Structures, DERA, Dunfirmline. Engseth, A.G. (1985). “Finite element collapse analysis of tubular steel offshore structures”, Div. of Marine Structures, Report UR-85-46, Norwegian Institute of Technology, Trondheim, Norway.

12.116


ISO 13819, (1994), ‘Petroleum and Natural Gas Industries – Offshore Structures – Part 1: General Requirements’, (1994), ‘Part 2: Fixed Steel Structures’, (2001), Draft,. Int. Standardization Organization, London. Fried, I. (1984) ”Orthogonal trajectory accession on the nonlinear equilibrium curve” J. Comp. Methods Appl. Mech. Engng. Vol.47, pp 283-297. Haugen,B. (1994) ” Buckling and stability problems for thin shell structures using high performance finite elements” Ph.D. thesis, Univ. of Colorado, Boulder. Hellan, Ø. (1995) “Nonlinear pushover and cyclic analyses in ultimate state design and reassessment of tubular steel Offshore Structures”, Dr.ing.thesis, Report MTA 1995: 108, Department of Marine Structures, Norwegian Institute of Technology, Norway. Hellan, Ø., Moan, T. and Drange, S.O. (1994) “Use of Nonlinear Pushover Analysis in Ultimate Limit State Design and Integrity Assessment of Jacket Structures” Proc. Behaviour of Offshore Structures, BOSS, Boston, Massachusetts, vol. 3, pp. 323-345. Horsnell, M.R. and Toolan, F.E., (1996), ‘Risk of Foundation Failure of Offshore Jacket Piles’, Paper No. 28 th Offshore Technology Conf. ., Houston, 381-392. Moan,T and Nordsve, N.T. (1979) “Numerical Prediction of Ultimate Behaviour of Marine steel Structures”, Int. Symp. on Marine Technology, NTH, Trondheim. Moan, T. and Amdahl, J. (2001). "Risk Analysis of FPSOs, with Emphasis on Collision Risk", Report RD 2001-12, American Bureau of Shipping, Houston (restricted). Moan, T., Amdahl, J. Engseth, A.G., and Granli, T. (1985), , “Collapse behaviour of trusswork steel platforms, “ Proc. Int. Conference on Behaviour of Offshore Structures, BOSS, Delft, Holland. Moan, T., Amdahl, Hellan, Ø. (2002), “Numerical Analysis for Ultimate and Accidental Limit State Design and Requlification of Offshore Platforms“ invited papers, Fifth World Congress on Computational Mechanics, Viena, Austria, ISBN 39501554-0-6, http://www.wccm.tuwien.ac.at NAFEMS (1992), “Introduction to Nonlinear Finite Element Analysis”, E.Hinton (ed.) NAFEMS, Glasgow. Nadim, F. and Dahlberg, R. (1996). ‘Numerical Modelling of Cyclic Pile Capacity in Clay’, Proc. OTC Conf., 1, pp. 347-356. Lehmann, E., Egge, E.D., Scharrer, M. and Zhang, L. (2001). ”Calculation of collisions with the aid of linear FE models”, Proceedings of PRADS’01, Shanghai, 1293-1300. Paik, J.K. and Thayamballi, A.K. (2002). “Ultimate limit state design of steel-plated structures”, John Wiley & Sons, Chichester, U.K. Paik, J.K. et al. (2003). “Report of ISSC Committee V.3, Collision and Grounding”, Proc. 15th ISSC, SanDiego, Published by Elsvier, Amesterdam. 12.117

Nonlinear Analysis

Recommend Documents