• Author / Uploaded
• abhinavka123

• Categories
• Deformación (Mecánica)
• Tensión (Mecánica)
• Pandeo
• Elasticidad (Física)
• Sistema no lineal

Citation preview

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

Chapter 12 - Nonlinear Analysis

12 Nonlinear Analysis 12.1 Introduction 12.1.1 General Linear versus 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

Chapter 12 - Nonlinear Analysis

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 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 load

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

Chapter 12 - Nonlinear Analysis

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

12.4

Chapter 12 - Nonlinear Analysis

Reasons for 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 serviceability 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 codebased 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 computerbased, ‘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

Chapter 12 - Nonlinear Analysis

a) Geometry

b) Deformation and equilibrium 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 r = sin α o cos α o A / cos α o A (ε is positive when the bar shortens)

ε =

and the axial force becomes:

S = EAε =

EA sin α o cos α o ⋅ r A

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 A

(12.1)

R = Kr

where

K=

2 EA 2 sin α o cos α o A

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

Chapter 12 - Nonlinear Analysis

If the angle αo is small (αo 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

Chapter 12 - Nonlinear Analysis

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: kT =

AE ep ⎡ 1 −1⎤ ⎢ −1 1 ⎥ L ⎣ ⎦

(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 where m =

a) Geometry and model

εY − ε D εA −εD

(12.83)

b) Load-displacement incrementation

Figure 12.18 Tapered bar subjected to axial loading

12.38

Chapter 12 - Nonlinear Analysis

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’ KT. When elastic KT = 3EA / A

(12.86)

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

(12.87)

When the second bar yields KT = EA / A

(12.88) 12.39

Chapter 12 - Nonlinear Analysis

and finally, after the last bar yields KT = 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 = ∫ ∫ [ε m − z ⋅ κ ) D(ε m − zκ ) dA]dx v

(12.90)

xA

For an elastic beam with D=E and z is defined with reference to the centroid the integral (12.90) may be written as 2~ ~ ~ ∫ ε DεdV = ∫ EAε m ε m dx + ∫ ( ∫ Ez κ κdA) dx v

x

x A

= ∫ EAε~m εdx + ∫ EIκ~κdx x

(12.91 )

x

12.40

Chapter 12 - Nonlinear Analysis

This is because for instance ∫ zε~mκdA = 0 , when z is referred to the centroid. 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: ~ ∫ zD ( z )ε mκdV v

σ N.A. H

ε a) stress-strain curve

∫ σ dA = 0 ∫ σ ⋅ zdA = M

(no axial force) elastic elasto-plastic b) stress distribution in beam section 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

Chapter 12 - Nonlinear Analysis

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 dσ ij = Dijkl dε kl

(12.92a)

ep Dijkl = E ijkl − β sij s kl

(12.92b)

where

Eijkl = 2Gε ije +

β= σ =

2νG δ ij ε kke 1 − 2ν

(12.92c)

9G 2 σ ( H '+3G )

(12.92d)

2

3 2

sij sij =

3 2

σ ij σ ij − 12 (σ kk ) 2

(12.92e)

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

n

n

n

ep ep Dijkl dε kl = ∑∑ Dijkl dε kl ; sij sij = ∑∑ sij2 ; σ k =1 l =1

i =1 j =1

kk

=

n

∑σ k =1

kk

; σ m = 13 σ kk . (12.92f)

sij denotes the deviatoric, or reduced stress tensor

sij = σ ij − 13 δ ij σ

kk

(12.92g)

where δ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

Chapter 12 - Nonlinear Analysis

⎛ E(1 − ν) 9G 2 ⎞ dσ xx = ⎜ − ⎟ dε xx (H '+ 3G) ⎠ ⎝ (1 + ν)(1 − 2ν) ⎛ 1.33E 2 ⎞ ≈ ⎜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

Chapter 12 - Nonlinear Analysis

σy(t) σy(c)

(a) Hysteresis loop

Cyclic stress-strain curve

(c)

(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

Chapter 12 - Nonlinear Analysis

(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

Chapter 12 - Nonlinear Analysis

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

Chapter 12 - Nonlinear Analysis

(a)

(b)

(c)

Formatted: English (U.S.) Formatted: English (U.S.)

Example problem N

Formatted: English (U.S.)

N

ks

(d)

(e)

Formatted: English (U.S.)

N

(f)

(g)

Formatted: English (U.S.) Formatted: English (U.S.)

Figure 12.22 Characteristic features of nonlinear response: Basic response patterns:

Formatted: English (U.S.) Formatted: English (U.S.)

(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

Chapter 12 - Nonlinear Analysis

Non-linear equations The resultant of internal forces can be expressed as R int = ∑ (a i ) T S

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 )

Si = R

(12.94a)

K I (r ) dr = dR

(12.94b)

i T

i

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, KI 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

Chapter 12 - Nonlinear Analysis

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

(12.95)

rm +1 = rm + Δrm +1 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 choice 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/A ratio) was presented. Also, the incremental stiffness, KI 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

Chapter 12 - Nonlinear Analysis

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., r2 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 (r1 and r2) results. Due to symmetry only r2 is needed to describe the problem in Fig. 12.2. The resulting stiffness matrix for the symmetric system in Fig. 2.16 is R 2 = 2(

EA 2 P 2 s − c ) r2 A/c A/c

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

θ3 θ5 − − 6 120 θ2 θ4 cos θ = 1 − + − 2 24 θ3 2θ5 tg θ = θ + + + 3 15 sin θ = θ −

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 A A

Global equilibrium is described by R = 2 Ps = 2 Pα 12.50

Chapter 12 - Nonlinear Analysis

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

Displacement: r j −1 = ∑ Δri Geometry: α j −1 ≅

h − r j −1 h

α0

j −1

Axial force: Pj −1 = ∑ ΔRi /( 2 sin α j −1 ) i =1

Next step is calculated as follows:

Pj−1 EA 2 ]Δrj = K I Δrj α j−1 − A A

Displacement increment

ΔR j = 2[

Update total displacement

rj = rj-1 + Δr j

Updated geometry

αj =

Updated axial force

Pj =

h − rj h

α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 EAα 03 Δλ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 = 0.02. The load incremental factor Δλ = 0.5 and 0.1.

12.51

Chapter 12 - Nonlinear Analysis

Fig. 12.24 Load-displacement relations obtained by load incrementation. Incrementation 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 Rm 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 (rm )-1 (R int (rm ) - R m ) = K I (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 Rint(rm). Euler-Cauchy’s method is based on a single point, rm by calculation of KI, and is denoted one-step method. In principle, multi-step methods are envisaged, but in practice it would not be optimal to go beyond a two-step method. Example 1- continued 12.52

Chapter 12 - Nonlinear Analysis

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, Rr 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

Rint may be calculated as follows Rint (r j ) = 2 Pj sin α j ≈ 2 Pj α j

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

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

Pj = EA ⋅ ≈

EA ΔA (A / cos α 0 − A / cos α j ) = A / cos α j A / cos α j

rj rj EA 2 EA (α 0 − α 2j ) = (2α 0 − ( ) 2 ) 2 2 A A

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

rj

rj rj − ( ) 2 )(α 0 − ) A A A

Δri and ΔRi are incremental

rj = ∑ Δri = ∑

ΔR i K I (ri −1 )

Since K I (ri −1 ) varies there is no simple explicit relationship between rj and R j = ∑ ΔRi The unbalanced force after increment j is

where

rj r j ⎞⎛ rj ⎞ ⎛ Rr =EA⎜⎜ 2α 0 − ( ) 2 ⎟⎟⎜⎜ α 0 − ⎟⎟ − R j A A ⎠⎝ A⎠ ⎝

j

r j = ∑ Δri i =1

,

j

R j = ∑ ΔRi i =1

12.53

Chapter 12 - Nonlinear Analysis

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 kI for each element i. Obtain

the

tangent

stiffness

matrix,

K I ( j −1) = ∑ (a i ) T k it ( j −1) .

Solve

K I ( j −1) Δr j = ΔR j for Δr j . From Δ rj obtain current strain increments Δε for each i j

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, ε ij = ε ij −1 + Δε ij and σ ij = σ ij + Δσ ij where Δσ j = ( 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

Chapter 12 - Nonlinear Analysis

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 = x n −

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)

θ xn f ( xn ) f ( xn ) = tgθ f '( xn ) Fig. 12.26 Newton-Raphson algorithm

This approach can be generalized to solve Eq. (12.92a): In this case KI(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 − R int )

(12.98)

or

rn +1 = rn − K I−1 (rn ) (R int − R) The basic principle for this iteration is illustrated in Fig. 12.27 system.

for a single d.o.f.

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

Chapter 12 - Nonlinear Analysis

R − R int = K I(n) Δrn +1

(12.98a)

in each iterative step. This is time-consuming. By updating KI 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 KI and one method where KI is updated after the first iteration.

Fig. 12.27 Newton-Raphson iteration.

a) No updating of KI

b) KI 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

Chapter 12 - Nonlinear Analysis

|| r ||=

1 N

∑ (r N

k =1

k

/ rref

)

2

(12.100 )

where N is the number of components in the vector r and rref 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 KI 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

Chapter 12 - Nonlinear Analysis

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

Chapter 12 - Nonlinear Analysis

g (r, λ ) = R int (r ) − λR ref = 0

(12.101 )

where Rref 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 = Rint - λRref = 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

Chapter 12 - Nonlinear Analysis

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

Chapter 12 - Nonlinear Analysis

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 pi =

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

(12.102)

S ip refers to increment No.i.

The initial value of S ip ( S 1p ) is 1.0. For stiffening system it will increase. For softening system it will decrease. If S ip 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

Chapter 12 - Nonlinear Analysis

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

⎛I ⎞ (12.103) ΔA n = ΔA 0 ⎜⎜ d ⎟⎟ ⎝ I0 ⎠ where ΔA 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 ip 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, Rref, is proportionally scaled via λ. For many practical structural problems, this loading 12.62

Chapter 12 - Nonlinear Analysis

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 Rref) 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 :

Mr (t) + Cr (t) + Kr (t) = R (t)

(12.106)

to solve the static problem : Kr = 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 , ri ,  ri , ri −1 , ri −1 ,  ri −1 ,...} 12.63

Chapter 12 - Nonlinear Analysis

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 {ri +1 ,  ri +1 , ri , ri ,  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. ri +1 = r0 (t ) + Δt ri +

ri −1 = ri − Δt ri +

Δt 2 Δt 3   ri + ri + ... 2 6

Δt 2 Δt 3   ri − ri + ... 2 6

(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 ri

(12.111)

Adding Eq. (12.109-110) yields :

ri +1 + ri −1 = 2r + Δt 2ri

(12.112)

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

ri =

1 {ri +1 − ri −1} 2Δt

(12.113)

 ri =

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

(12.114)

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

12.64

Chapter 12 - Nonlinear Analysis

1 ⎫ 1 1 ⎧ 1 C ⎬ ri +1 = R i (t ) − K ri (t ) + 2 M {2ri − ri −1} + C ri −1 ⎨ 2 M+ 2 Δt ⎭ 2 Δt Δt ⎩ Δt (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:

⎧ ⎪ 1 ⎨ 2 ⎪ Δt ⎩

0 ⎤ ⎡ M 11 0 ⎡C11 0 0 ⎤ ⎫ ⎡ r1(i +1) ⎤ ⎡ R1(i ) ⎤ ⎡ K11 K12 K13 ⎤ ⎡ r1( i ) ⎤ 1 ⎢ ⎪⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ 0 M ⎥ 0 ⎥+ 0 C22 0 ⎥⎥ ⎬ ⎢ r2(i +1) ⎥ = ⎢ R2( i ) ⎥ − ⎢⎢ K 21 K 22 K 32 ⎥⎥ ⎢ r2( i ) ⎥ 22 ⎢ ⎢ 2Δt ⎢⎣ 0 ⎢⎣ 0 0 C33 ⎥⎦ ⎪⎭ ⎢⎣ r3( i +1) ⎥⎦ ⎢⎣ R3(i ) ⎥⎦ ⎢⎣ K 31 K 32 K 33 ⎥⎦ ⎢⎣ r3(i ) ⎥⎦ 0 M 33 ⎥⎦

0 ⎤ ⎧ ⎡ r1( i ) ⎤ ⎡ r1(i −1) ⎤ ⎫ ⎡ M 11 0 ⎡C11 0 0 ⎤ ⎡ r1(i −1) ⎤ ⎪ ⎢ ⎥ ⎢ 1 ⎢ ⎥⎪ 1 ⎢ ⎢ ⎥ ⎥ + 2 ⎢ 0 M 22 0 ⎥ ⎨ 2 ⎢ r2( i ) ⎥ − ⎢ r2( i −1) ⎥ ⎬ + 0 C22 0 ⎥⎥ ⎢ r2( i −1) ⎥ ⎢ Δt 2Δt ⎢⎣ 0 ⎢⎣ 0 0 C33 ⎥⎦ ⎢⎣ r3(i −1) ⎥⎦ 0 M 33 ⎥⎦ ⎪⎩ ⎢⎣ r3(i ) ⎥⎦ ⎢⎣ r3(i −1) ⎥⎦ ⎪⎭ (12.116) The first equation in Eq. (12.121) is explicitty written as : 1 ⎧ 1 ⎫ C11 ⎬ r1( i +1) = R1 (t ) − K11r1( i ) − K12 r2( i ) − K13 r3( i ) ⎨ 2 M 11 + Δ Δ 2 t t ⎩ ⎭ 1 1 + 2 M 11 {2r1(i ) − r1(i −1) } + C11r1(i −1) Δt 2 Δt

(12.117)

This shows that ri ( i +1) can be directly, explicity determined by the response at time t. There is no coupling between displacements, rj (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
0 inadmissible

Taking the differential of f from a plastic state gives

df =

∂f ∂f ∂f dσ ij + p dε ijp + dκ ∂σ ij ∂κ ∂ε ij

(C.13)

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

∂f dσ ij < 0 ; f = 0 during unloading ∂σ ij ∂f dσ ij = 0 ; f = 0 during naytral loading ∂σ ij

(C.14a-c)

∂f dσ ij > 0 ; f = 0 during loading ∂σ ij 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ε ijp = 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

Chapter 12 - Nonlinear Analysis

Since,

κ = κ (ε ijp )

(C.16)

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

(

f = f σ ij , ε ijp

)

(C.17)

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

∂f ∂f dσ ij + p dε ijp = 0 ∂σ ij ∂ε ij

(C.18)

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

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

(C.19)

From Eqs (C.2, C.8)

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

(C.20)

and ∂f ∂ ∂ = σ − H (ε p ) = − p H (ε p ) p p ∂ε ij ∂ε ij ∂ε ij

(

)

(

)

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

(C.21)

where the plastic work, Wp (Eq. (C.7)) has been utilized. By noting ∂f ∂σ = ∂σ ij ∂σ ij

and σijsij = σijσij − 13 σijδijσ m = σijσij − 13 σ2m = 23 σ 2

12.110

Chapter 12 - Nonlinear Analysis

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

dλ = −

∂σ dσ ij ∂σ ij

σ ⎛ ⎜⎜ − H ' ij σ ⎝

⎞ ⎛ 3 sij ⎞ ⎟⎟ ⎜⎜ ⎟⎟ ⎠⎝ 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εijp =

3 sij 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σ =

∂σ 3 dσ kl = s kl dσ kl ∂σ kl 2σ

(C.24)

Then by inserting in Eq. (C.23) dεijp = α ⋅ sijs kl dσ kl

(C.25)

where

α=

9 4σ H ' (σ ) 2

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

Chapter 12 - Nonlinear Analysis

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

(C.26)

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

= −

Introducing the following relations:

Eijkl s kl = 2Gsij

(C.27)

Eijkl = E klij 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ε ijp =

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)

(

= Eijkl dε kl − dε klp

)

(C.30)

ep = Dijkl dε kl

where ep Dijkl = Eijkl −

9G 2 sij 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

Chapter 12 - Nonlinear Analysis

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

The following assumption are made: Strains:

γ xy = dγ xy = 0 ; γ yz = dγ yz = 0

(C.32)

γ zx = dγ zx = 0 Stresses:

σ xy = dσ xy = 0 ; σ yz = dσ yz = 0

(C.33)

σ yy = dσ yy = 0 ; σ zz = dσ zz = 0

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

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

D xxyy D yyyy D zzyy

D xxzz ⎞ ⎛ dε xx ⎞ ⎟⎜ ⎟ D yyzz ⎟ ⎜ dε yy ⎟ D zzzz ⎟⎠ ⎜⎝ dε zz ⎟⎠

(C.34)

or ep dσ ij = Dijkl dε kl

(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

yyxx

D zzzz − D zzxx D yyzz )

D xxzz (Dzzxx D yyyy − D yyxx D yyzz ) ) dε xx 2 D yyyy D zzzz − D yyzz

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

Chapter 12 - Nonlinear Analysis

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

D xxxx = E xxxx −

9G 2s 2xx E(1 − ν) 4G 2 σ2xx = − σ 2 (H '+ 3G) (1 + ν)(1 − 2ν ) σ 2 (H '+ 3G)

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

D xxyy = D yyxx = D xxzz = D zzxx

D yyyy = D zzzz =

E(1 − ν) Gσ − (1 + ν)(1 − 2ν) σ 2 (H '+ 3G)

D yyzz = D zzyy =

Eν G 2 σ 2xx − 2 (1 + ν)(1 − 2ν) σ (H '+ 3G)

2

(C.38)

2 xx

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

(C.40)

Stresses: σ zx = dσ zx = 0

;

σ yz = dσ yz = 0

(C.41)

σ zz = dσ zz = 0

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 yyyy D xyyy D zzyy

D xxxy D yyxy D xyxy D zzxy

D xxzz ⎞ ⎟ D yyzz ⎟ D xyzz ⎟ ⎟ D zzzz ⎟⎠

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

(C.42)

or

dσij = Dijkl dε kl + Dijzz dε zz

(C.43)

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

(C.44)

and

This equation system may be contracted to the following expression 12.114

Chapter 12 - Nonlinear Analysis

Dijzz D zzkl ⎛ dσ ij = ⎜⎜ Dijkl − D zzzz ⎝

⎞ ⎟⎟ dε kl ⎠

(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 zz = − 13 ( σ xx + σ yy ) =− σ m

s xy = σ xy

(C.46)

In Eq. (10.108) the stiffness matrix elements are D xxxx = E xxxx −

9G 2s 2xx E(1 − ν) 9G 2s 2xx = − σ 2 (H '+ 3G) (1 + ν)(1 − 2ν ) σ 2 (H '+ 3G)

D xxyy = D yyxx =

9G 2s xx s yy Eν − 2 (1 + ν)(1 − 2ν) σ (H '+ 3G)

D xxzz = D zzxx =

Eν 9G 2s xx σ m + 2 (1 + ν)(1 − 2ν) σ (H '+ 3G)

D xxxy = D xyxx = − D yyyy

9G 2s xx σ xy σ 2 (H '+ 3G)

9G 2s 2yy E(1 − ν) − = (1 + ν)(1 − 2ν) σ 2 (H '+ 3G)

D yyzz = D zzyy =

9G 2s yy σ m Eν + 2 (1 + ν)(1 − 2ν) σ (H '+ 3G)

D yyxy = D xyyy = − D zzxy = D xyzz = D xyxy = G − D zzzz =

9G 2s yy σ xy σ 2 (H '+ 3G)

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

9G 2 σ2xy

(C.47)

σ 2 (H '+ 3G)

E(1 − ν) 9G 2 σ 2m − 2 (1 + ν)(1 − 2ν ) σ (H '+ 3G)

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

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

(C.48)

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

p

=

4 3

(dε

p2 xx

2

)

2

+ dε yyp + dε xxp dε yyp + 13 dγ xyp

(C.49) 12.115

Chapter 12 - Nonlinear Analysis

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

Chapter 12 - Nonlinear Analysis

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. 28th 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

Chapter 12 - Nonlinear Analysis

Ramm, E. (1982) “The Riks/Wempner approach – an extension of the displacement control method in non-linear analysis”, Non-linear Computational Mechanics, ed. E. Hinton et al., Pineridge Swansea, pp. 63-86. Ramm, E., (1981)“Strategies for tracing the non-linear response near limit-points, Non-linear Finite Element Analysis in structural mechanics”, ed. W. Wunderlich, Springer-Verlag, Berlin, pp 63-89. Riks, E. (1972) “The application of Newton’s method to the problem of elastic stability”, J. Appl. Mech. 39, 1060-6. SCI(1993) “Interim Guideance Notes for the Design and Protection of Topside Structures against Explosion and Fire”, Stell Construction Institute, Document SCI-P112/503, London. Skallerud, B. and Amdahl, J. (2002) “Nonlinear analysis of offshore structures”, Baldock: Research Studies Press, 2002 Simonsen, B.C. and Lauridsen, L.P. (2000). “Energy absorption and ductile fracture in metal sheets under lateral indentation by a sphere”, International Journal of Impact Engineering, 24, 1017-1039. Søreide, T.H., Amdahl, J. Granli, T. and Astrup, O.C. (1986), “Collapse analysis of framed offshore structures”, Paper OTC 5302 Offshore Technology Conference, OTC, Houston, Texas. Søreide, T.H., Amdahl, J., Eberg, E., Holmås, T. and Hellan, Ø. (1994), “USFOS – A computer program for progressive collapse analyses of steel offshore structures; Theory manual”, report STF71 F88038, 6th revision, SINTEF Structures and Concrete, Trondheim, Norway. Taby, J. and Moan, T. (1985) “Collapse and Residual Strength of Damaged Tubular Members”, Proc. 4th In. Conf. on the Behaviour of Offshore Structures, BOSS Conf., Elsevier, Amsterdam. Taylor, R. and Zienkiewicz, O.C. (2000), “The finite element method”,Vol. 1,2,3, 5th Edition Ultiguide (1999), “Best Practice Guidelines for Use of Nonlinear Analysis Methods in Documentation of Ultimate Limit States for Jacket Type Offshore Structures", DNVSINTEF-BOMEL, Oslo - Trondheim - London. Vinnem, J.E. (1999) “Offshore Risk Assessment”, Kluwer Academic Publishers, Doordrecht. WOAD: “Worldwide Offshore Accident Databank”, Det Norske Veritas, Oslo, 1996. Wempner, G.A. (1971), “Discrete approximations related to nonlinear theories of solids”, Int. J. Solids & Structs., 7, 1581-1599. Wolf, J.P. (1994). ‘Foundation Vibration Analysis using Simple Physical Models’, PTR Prentice Hall, Englewood Cliffs. 12.118