FEA of Nonlinear Problems 2011 Ivanco

NONLINEAR FINITE ELEMENT ANALYSIS SCRIPT OF LECTURES

Doc. Ing. Vladimír Ivančo, PhD. Faculty of Mechanical Engineering, Technical University of Košice, Slovakia

HS Wismar, June 2011

2

NONLINEAR FINITE FINITE ELEMENT ANALYSIS ANALYSIS

CONTENTS 1.

Structural nonlinearities................... nonlinearities ......................................... ............................................ ..............................4 ........4 1.1 Introduction.......................................................................................4 1.2 1.3

2.

Types of structural nonlinearities.......................... nonlinearities................................................ ............................5 ......5 Concept of time curves ......................................... ............................................................... ............................5 ......5 Geometrically nonlinear finite element analysis .................................7 .................................7

2.1 2.2

Large displacement and small strain behavior..................... behavior ..................................7 .............7 Incremental - iterative solutions.................................................. solutions...................................................... ....12 12 2.2.1 Incremental method............................................. method...........................................................14 ..............14 2.2.2 Newton-Raphson Newton-Raphson method................................. method..................................................15 .................15 2.2.3 Modified Newton-Raphson method..................................16 2.2.4 Quasi-Newton methods....................................... methods.....................................................17 ..............17

2.3 2.4

Linear stability analysis ............................................ ..................................................................17 ......................17 Large displacement and large strain behavior.................................18 behavior.................................18 2.4.1 Total Lagrangian formulation ........................................ ........................................... ...19 19 2.4.2 Updated Lagrangian formulation ......................................19 ......................................19 Material nonlinearities ........................................... ................................................................. ............................ ......20 20

3. 3.1 3.2

Introduction.....................................................................................20 Nonlinear elasticity models.................................................. models.............................................................20 ...........20

3.3

Elastoplastic material model ......................................... ...........................................................22 ..................22 3.3.1 Yielding criterion ......................................... ..............................................................22 .....................22 3.3.2 Post yielding behavior.............................. behavior.................................................... ......................... ...24 24 3.3.3 Constitutional equations of elastoplastic material.............26 3.3.4 Integration of constitutive equations.................................28 3.3.4.1 Generalised trapezoidal rule................................................29 3.3.4.2 Generalised mid point rule..................................................29

3.3.5 Numerical procedures ......................................... .......................................................30 ..............30 Appendix........................................................................ Appendix.................................................. ............................................ ....................................32 ..............32 Nonlinear analyses analyses with COSMOS/M program............................... program..........................................33 ...........33

NONLINEAR FINITE ELEMENT ANALYSIS ANALYSIS

3

1.1 Properties of finite elements..................... elements ........................................... ........................... ..... 33 1.2 Nonlinear analysis setup...................................................... setup...................................................... 36 1.3 Time curves and time parameters................................... parameters........................................ ..... 42 1.4 Restart possibilities.......................................... possibilities.............................................................. .................... 47 Examples of of COSMOS/M program use use ........................................... ..................................................... .......... 48 1.5 Compressed slender beam .......................................... ................................................... ......... 48 1.6 Thick-walled pipe subjected to internal pressure ................ 51 1.7 Deformation of a thin-walled tank ...................................... ...................................... 55

4

NONLINEAR FINITE ELEMENT ANALYSIS

1. STRUCTURAL NONLINEARITIES 1.1 Introduction Solution of many engineering problems is based on linear approximations. In structural analyses, these approximations are represented by consideration that • displacements are small and can be neglected in equilibrium

equations, • the strain is proportional to the stress (linear Hookean material

model), • loads are conservative, independent on displacements, • supports of the structure remain unchanged during loading.

Consequently, in the Finite Element Analysis (FEA) the set of equations, describing the structural behavior is then linear K d = F ,

(1.1)

where K is the stiffness matrix of the structure, d is the nodal displacements vector and F is the external nodal force vector. Characteristics of solution of this linear problem is that • the displacements are proportional to the loads, • the stiffness of the structure is independent on the value of the

load level. In reality, behavior of structures is nonlinear, but divergences from linear response are usually small and may be neglected in most practical problems. On other hand, solution of many engineering problems needs abandonment of linear approximations. For example, displacements of slender structures (like crane towers, masts etc.) may be so large that changes of the structure shape (or configuration changes) cannot be neglected. Many materials behave nonlinearly or linear material model cannot be used if stress exceeds some value. Moreover, loads may change their orientations according to displacements and supports may change during loading. Consequently, structure behaves nonlinearly. If these phenomena are included in a FEA, the set of equilibrium equations becomes nonlinear and instead of set of linear equations (1.1) we obtain a set of nonlinear algebraic equations R (d )

= F .

(1.2)


5

1.2 Types of structural nonlinearities Structural nonlinearities can be specified as 1. Geometrical nonlinearities: The effect of large displacements on the overall geometric configuration of the structure. 2. Material nonlinearities : Material behavior is nonlinear. Possible material models are: a) nonlinear elastic, b) elastoplastic, c) viscoelastic, d) viscoplastic. 3. Boundary nonlinearities, i.e. displacement dependent boundary conditions. The most frequent boundary nonlinearities are encountered in contact problems. Consequences of nonlinear structural behavior that have to be recognized are: a) The principle of superposition cannot be applied. Thus, for example, the results of several load cases cannot be combined. Results of the nonlinear analysis cannot be scaled. b) Only one load case can be handled at a time. c) The sequence of application of loads (loading history) may be important. Especially, plastic deformations depend on a manner of loading. This is a reason for dividing loads into small increments in nonlinear FE analysis. d) The structural behavior can be markedly non-proportional to the applied load. e) The initial state of stress (e.g. residual stresses from heat treatment, welding, cold forming etc.) may be important.

1.3 Concept of time curves For nonlinear static analysis, the loads are applied in incremental steps using time curves. The “time” variable represents a pseudo time, which denotes the intensity of the applied loads at certain step. For nonlinear dynamic analysis and nonlinear static analysis with timedependent material properties,1 “time” represents the real time associated with the loads’ application. 1

i.e. analysis of creep and relaxation problems by use of viscoelastic or viscoplastic material models.

6


As an example, time curves of forces F 1 and F 2 loading simple beam are displayed in Figure 1.1. Values of forces at any time are defined as F 1 = λ 1 (t ) f 1

and

F 2 = λ 2 (t ) f 2

where f 1 and f 2 are input values of forces and λ 1 and λ 1 are load parameters that are functions of time t .

Figure 1.1: Example of time curves The choice of time step size depends on several factors such as the level of nonlinearities2 of the problem and the solution procedure. Generally, sufficiently small steps are necessary to simulate nonlinear response of a structure with satisfactory accuracy. On the other hand, large number of too small time steps uselessly increases consumption of CPU time. Computer programs are usually equipped with an adaptive automatic stepping algorithm to facilitate the analysis and to reduce the solution time demands.

2

Highly nonlinear problems need smaller load increments.

7


2. GEOMETRICALLY NONLINEAR FINITE ELEMENT ANALYSIS 2.1 Large displacement and small strain behavior To examine geometrically nonlinear behavior we will start with an example. We assume large displacement, but small rotation and, what is the most important, small strain. The structure is very simple – only one bar truss as is shown in Figure 2.1. At the beginning, when the force P is zero, the axial force N in the bar is zero too and bar has its initial length L0.

Figure 2.1: Example of nonlinear structure – single bar truss

Using the free body diagram shown in Figure 2.1 the equilibrium equation is N sin α − P = 0

and after substituting sin α = ( h + u ) L N

h+u L

− P = 0.

(2.1)

Assume that material is linearly elastic with Young’s modulus E. The assumption of small strains means here that changes of the bar cross sectional area A can be neglected. Then axial force in the bar is

8


N = E A0 ε

(2.2)

where A0 is the initial cross sectional area and ε is the engineering strain defined as ε =

L − L0 L0

.

(2.3)

As lengths are given as L0 =

a2 + h2

and

L =

a 2 + (h + u )2

(2.4)

the expression for strain is getting rather complicated. We can overcame this problem by introducing Green’s strain defined as L2 − L20 ε G = 2 L20

(2.5)

which for our problem becomes 2

1  u  ε G = +   . L0 L0 2  L0  h

u

(2.6)

Use of this new measure of strain is possible because we can define strain arbitrarily. The only condition is that the strain measure must be objective, which means that is have to be independent on choice of coordinate system and insensitive to a rigid body movement. From equations (2.3) and (2.5), it follows that L2 − L20 L − L0 L + L0 ε G = = = 2 2 L0 2 L0 L0

 L + L0 1 L + L0 1  ε = ε  + 1 − 1 =  2 L0 2  L0 



L − L0 2 L0  1   = ε + 2  L0 L0  

or 1 2 ε . 2 Noting that the constitution equation was measured as ε G = ε +

σ =

N A0

= E ε

the same constitutive equation when using Green’s strain should be

(2.7)

(2.8)

9


σ = E

ε ε G

ε G = E

ε

1 2 ε + ε 2

ε G =

E

1 1 + ε 2

ε G .

(2.9)

This means that we should use value E ∗ =

E

1 + (1 / 2)ε

instead of E in the constitution equation. Fortunately, we can ignore this complication now because for small engineering strain is the difference between engineering and Green’s strain negligible. For example, consider that ε = 0,002 (e.g. mild steel yields at about this value), then ε G = 0,002 + 0,5 ⋅0,0022 = 0,002002. This means that difference is only 0,1% i.e. a value that can be usually neglected. Assuming that strain is small, we can write σ ≈ E ε G and according to equation (2.6) N =

E A0  1  h u + u2  .  2 2  L0 

(2.10)

Substituting (2.10) to equilibrium equation (2.1) and assuming that for small strain is L ≈ L0 gives the equilibrium equation

)

(

E A0 3 u + 3 h u 2 + 2 h 2u = P . 3 2 L0

2.11)

Obviously, the equation is nonlinear with respect to displacement u. That means that relation between load P and displacement u is represented not by a straight line as it is when changes of configuration are neglected but by a curve. This nonlinear characteristic for E = 2,1⋅105 MPa, A0 = 100 mm2, a = 200 mm and h = 20 mm is shown in Figure 2.2. P [N]

2500

2000

1500

1000

500

0 -10

-5

0

5

10

u [mm] -500

-1000

Figure 2.2: Geometrically nonlinear behavior of a single bar truss

10


There is another possibility to obtain equation of equilibrium (2.1) or (2.11). From principle of virtual displacements, it follows that when the structure is in equilibrium, virtual works of internal and external forces are equal for every kinematic admissible set of virtual displacements. For our structure with one

degree of freedom, only one virtual displacement virtual displacements has a form

δ u is

possible and principle of

∫ σ δε G dV = P δ u

(2.12)

V

where δε G is virtual strain corresponding to virtual displacement strain can be expressed from equation (2.6) as dε d δε G = G δ u = du du

 h u 1  u 2  h+u  +    δ u = 2 δ u.  L0 L0 2  L0   L0  

δ u .

The virtual

(2.13)

It is assumed in principle of virtual displacements that virtual displacement is infinitesimal and hence the stress σ = ( N / A) remains unchanged. Noting that σ and δε G are constant over the whole volume V in this case and assuming that changes of the volume can be neglected due to small strain, i.e. V ≈ V 0 = A0 L0 , equation (2.12) becomes N h + u δ u A L0 = P δ u A0 L20

and from this equation it follows that N

h+u

= P .

L0

This is the same equation as the equation of equilibrium (2.1). After substituting for N from (2.10) the equation (2.11) will be received again. Utilization of principle of virtual displacements (PVD) is a convenient way to obtain conditions of equilibrium for complex structures. For general threedimensional case we have three components of displacement u, v, w and six components of Green’s strain 2 2 2 ∂u 1  ∂u   ∂v   ∂w   ε x = +   +   +    ∂ x 2  ∂ x   ∂ x   ∂ x    

,

2 2 2 ∂v 1  ∂u   ∂v   ∂w   +   +   +    ε y = ∂ y 2  ∂ y   ∂ y   ∂ y    

,

2 2 2 ∂w 1  ∂u   ∂v   ∂w   ε x = +   +   +    , ∂ z 2  ∂ z   ∂ z   ∂ z    

(2.14)

11


γ xy =

∂u ∂v ∂u ∂u ∂v ∂v ∂w ∂w , + + + + ∂ y ∂ x ∂ x ∂ y ∂ x ∂ y ∂ x ∂ y

γ xz =

∂u ∂w ∂u ∂u ∂v ∂v ∂w ∂w + + + + , ∂ z ∂ ∂ ∂ z ∂ ∂ z ∂ x ∂ z

γ yz =

∂v ∂w ∂u ∂u ∂v ∂v ∂w ∂w . + + + + ∂ z ∂ y ∂ y ∂ z ∂ y ∂ z ∂ y ∂ z

In finite element method are displacement interpolated within the finite elements as u=

∑ N i ui ,

∑ N i vi ,

v=

i

w=

i

∑ N i wi ,

(2.15)

i

where ui, vi, wi are nodal displacements and N i are shape functions. Substituting these equations into expressions of Green’s strain components, we obtain ε = (B L

+

1 B ) d . 2 N

(2.16)

In matrix equation (2.16)  ε x   ε   y   ε z  ε =  , γ yz   γ xz    γ xy 

(2.17)

and d is matrix of nodal displacements. Matrix B L is the usual small displacement matrix and matrix B N reflects the fact that Green’s strain is a nonlinear function of displacements. Elements of this matrix are linear functions of nodal displacements d . It might be shown that virtual strain corresponding to the virtual nodal displacements δ d is δ ε = (B L + B N ) δ d = B δ d .

(2.18)

According to the principle of virtual displacements, virtual work of internal forces must be equal to virtual work of external forces if the structure is in equilibrium. This is represented by the equation

∫

δ ε T σ dV = δ d T F

V

where F is matrix of nodal forces. We suppose linear relation between stress and strain components, hence σ = Dε

where D is matrix of material elastic constants.

(2.19)

12


Substituting (2.18) into (2.19) gives

∫

δ d T BT σ d V = δ d T F

(2.20)

V

for any kinematic admissible set of virtual displacements δ d . Then

∫

T

B σ d V = F

.

(2.21)

V

The last equation is a matrix representation of a set of nonlinear algebraic equations for unknown nodal displacements d . R (d )

= F .

(2.22)

2.2 Incremental-iterative solutions We have seen that assumption of large displacements leads to nonlinear equation of equilibrium (2.1) or (2.11) for a simple bar truss example. Generally, in finite element analysis we have a set of nonlinear equations (2.22). Let us start with the bar-truss example. The equation of equilibrium (2.1) or (2.11) can be written in a form R (u ) = P

(2.23)

where R (u ) = N

h+u L0

=

)

(

E A 3 u + 3 h u 2 + 2 h 2u . 3 2 L0

(2.24)

represents a component of internal force. The basic step to solve the nonlinear equation (2.24) is a linear approximation for small increment of force and corresponding increment of displacement. Assume that for a prescribed value of force P we managed to find (e.g. by error and trial method) a displacement u satisfying the equation (2.23). Internal force R (u + d u ) for new external force P 1 ( P + d P ) can be approximated by the linear function  d R   du u d   u

R (u + d u ) = R (u ) + 

and approximate condition of equilibrium is  d R   d u = P + d P .  d u  u

R (u ) + 

Assuming equation (2.23) gives  d R    d u = d P u d   u

(2.25)

13


or K T ( u ) d u = d P

(2.26)

where  d R    d u  u

K T = 

(2.27)

is called the tangent stiffness. For the particular case of the bar-truss, tangent stiffness can be easily found as K T =

 u + h  d N d  u + h    N +   . L u d u  L0  d  0 

Using the equation (2.10) gives d N EA u + h = d u L0 L0 from which K T = K 0 + K u + K σ

(2.28)

where K 0 =

E A  h 

2

  L0  L0 

E A  u

 u  K u = 2 +   L0  h  h   K σ =

N L0

is the linear stiffness 2

 h       L0 

2

is the initial displacement stiffness

is the initial stress stiffness.

The linear stiffness, which is independent on displacement, is familiar from small displacement structural analysis. The initial displacement stiffness reflects the effect of displacement on stiffness3. The initial stress stiffness reflects the fact that there is an axial force in the bar prior to load increment. In like manner, we can precede in a general case described by the equation (2.21) or (2.22) and derive R (d + d d )

= R (d ) + K T d d

and K T dd =

d F

(2.29)

where 3

e.g. it can be seen from the diagram in Figure 2.2 that for compressive load P stiffness decreases and for tensional force P increases.

14


K T

∂R ∂d

=

is the tangent stiffness matrix. We can also find out that K T

= K 0 + K u + K σ .

(2.30)

where K T is linear stiffness matrix, K u is initial displacement stiffness matrix and K σ is initial stress stiffness matrix. Introduction of tangent stiffness matrix is crucial for solution of nonlinear equations (2.22). The most widely used methods are briefly introduced in the following text: 2.2.1 Incremental method

The load is divided into a set of small increments ∆F i . Increments of displacements ∆d i are calculated from the set of linear simultaneous equations K T (i −1)

∆d i = ∆F i

.31)

and an updated solution is obtained as d i

= d i −1 + ∆d i .

Figure 2.3: Incremental method

32)

15


The procedure is shown in Figure 2.3. It is obvious that solution error, i.e. difference from exact solution gradually cumulates. To reduce error, large number of small incremental steps has to be done that is inefficient. On the other hand, division of loading process into sufficiently small increments is necessary to model load path dependent behavior of a structure. Dependence of response on a manner of loading, not only of final values of loads is typical for problems with plastic deformation and with friction. In these problems, incremental method is usually combined with one of following methods. 2.2.2 Newton-Raphson method

Suppose that initial displacements d 0 are known. The first guess of nodal displacements for load F is calculated by solving set of linear algebraic equations K T ( 0) d 1

= F

(2.33)

where K T ( 0)

= K T (d 0 )

is tangent stiffness matrix calculated for initial displacements. As the displacements d 1 are most probably not accurate, the equilibrium equation (2.22) is not satisfied and R (d 1 )

≠ F

that means there are unbalanced (or residual) nodal forces r 1

= R (d 1 ) − F .

(2.34)

By computing new tangential stiffness matrix K T (1)

= K T (d 1 )

and solving new set of algebraic linear equations K T (1)

∆d 1 = r 1

(2.35)

we will obtain an improved solution d 2

= d 1 + ∆d 1 .

(2.36)

If r 2 = R (d 2 ) − F ≠ 0 the procedure is repeated until the sufficiently accurate solution is obtained. The iterations are schematically shown in Figure 2.4. This method, known as Newton-Raphson 4 method (NR) is often combined with incremental method as displayed in Figure 2.5.

4

Joseph Raphson (1648-1715) was an English mathematician, a Fellow of the Royal Society of London and friend of Newton.

16


Figure 2.4: Standard Newton-Raphson (NR) method.

Figure 2.5: Combination of Newton-Raphson and incremental methods. 2.2.3 Modified Newton-Raphson method

The standard Newton-Raphson method, although effective in many cases, needs the solution of the set of linear equations (2.35) which is time demanding for large systems. Modified Newton-Raphson method (MNR) differs from standard NR algorithm in that the stiffness matrix is only updated occasionally. In the example shown in Figure 2.6, the tangential stiffness matrix is formed and decomposed at the beginning and used throughout the iterations. Advantage of the method is in saving computer time, because factorization of the tangent stiffness matrix is


17

performed only once for the load increment. On the other hand, number of iterations needed is usually larger.

Figure 2.6: Modified Newton-Raphson(MNR) method.

2.2.4 Quasi-Newton methods

There exist many other methods for solution of the set of nonlinear algebraic equations, so called quasi-Newton methods. The most popular among them is Broyden – Fletcher – Goldfarb – Shanno (BFGS) method.

2.3 Linear stability analysis Theoretically, below a certain critical load a structure is in position of stable equilibrium, whilst above that load the equilibrium may be unstable. Unstable equilibrium means that though the structure is in equilibrium, any arbitrary small disturbance will cause loss of this equilibrium. In many practical problems, the displacements are small for load less than critical and behavior of the structure can be considered as a linear function of applied load. The typical example is Euler strut buckling, Figure 2.7. For axial force N that is less than critical, the strut is in stabile equilibrium. This equilibrium is possible if a lateral load P then deflects the strut as well. If the lateral load is removed, the strut will return to its straight shape.

18


If the force N is greater than critical, the strut can remain (theoretically) straight but its equilibrium is unstable, any small lateral load will cause deflection increasing until the collapse5.

Figure 2.7: Buckling of a strut

For load less than critical small longitudinal (in plane) and lateral displacements allow the initial displacement stiffness matrix K u to be ignored. The equilibrium equation can be written as [K 0 + λ K σ (N )] d = P

(2.37)

The elastic critical (buckling) load is given by the lowest value of load parameter λ for which d ≠ 0 when the lateral load P = 0 . Physically this means that equilibrium is possible with very small lateral displacements in the absence of any lateral load. In mathematical sense, we have to solve the eigenvalue problem [K 0 + λ K σ (N )] d = 0

(2.38)

where λ is the eigenvalue and d is the corresponding eigenvector. It should be noted that due to assumptions accepted the solution represents itself only an estimation of the upper bound of the structure load capacity.

2.4 Large displacement and large strain behavior When strain is large, it is inadmissible to neglect shape and volume changes of a structure. For example, in the simple bar example we have to introduce current cross sectional A instead of initial A0 and current length L instead of initial length L0 in the equations (2.10) and (2.11). Accordingly, integration in the equation (2.19) expressing the principle of virtual displacements has to be taken over the current volume. This brings problems, as the current volume is unknown, because it depends on displacements that are unknown too and must be calculated first. To solve this problem, it is necessary to introduce a transformation so that integrals are taken over known volume. Two possible ways are briefly described bellow:

5

In reality, unstable equilibrium is due to initial imperfections (e.g. eccentricity of force N , initial curvature of the strut etc.) impossible, but estimation of critical load may be useful in many cases.

19


2.4.1 Total Lagrangian formulation

In a Total Lagrangian (TL) formulation all integrals are calculated with respect to the initial undeformed configuration of the structure

∫ δ ε G σ P d V = δ d T

T F

(2.39)

V 0

where V 0 is the initial volume. Due to transformation, new measure for stress so called second Piola-Kirchhoff stress tensor σ σ P has to be introduced with Green’s strain tensor ε εG . 2.4.2 Updated Lagrangian formulation

In an Updated Lagrangian (UL) formulation, a known deformed configuration i is taken as an initial state for subsequent configuration ( i+1) and this is continually updated as the calculation proceeds

∫ δ

(i +1) T ( i +1) ε σ A

T C d V = δ d F

(2.40)

V 0

In the left side of the equation (2.40), σ σC is Cauchy stress tensor and ε ε A is Almansi (i+1) (i+1) ε σ strain tensor respectively. Notation ε A and σC means that the strain and stress are in configuration (i+1). Integration is done over volume V i that is in current configuration i. Use of different measures for stress and strain in TL and UL formulation follows condition that virtual work of internal forces must be the same irrespective of the volume over which is integration taken 6.

6

That means that stress and strain measures must be work conjugate.

20


3. MATERIAL NONLINEARITIES 3.1 Introduction Linear elastic FE analysis is based on linear constitutive stress-strain equations σ = D ε

(3.1)

in which the terms of material matrix D are expressed as functions of constant values of modulus of elasticity and Poisson’s ratio. The constant D matrix leads to a constant stiffness matrix K , which is for strain-displacement relationship ε = B d

(3.2)

given by K

=

∫

T B DB

d V

(3.3)

V

Departure from linear elasticity implies that the linear elastic constitutive equations are no longer valid, as the material matrix is no longer constant. The non-constant material matrix D represents nonlinear constitutive equations corresponding to the adopted nonlinear material model. Consequently, the conditions of equilibrium derived in FEM from principle of virtual displacements are nonlinear like equations (2.21) and (2.22). Solution of these equations is based on the same methods as in geometrically nonlinear case. Usually it is necessary to divide load into increments and perform equilibrium iterations (e.g. by MNR or NR method) for each increment. Moreover, for each load increment there must be performed stress iterations, as the material matrix is function of strain. The strain is unknown a priori and will be computed only. Material nonlinearities are often combined with geometrical and/or boundary nonlinearities.

3.2 Nonlinear elasticity models Nonlinear elastic behavior of materials can be formulated in several ways. The simplest is total formulation, where the stress and strains are defined in terms of the secant modulus of elasticity E s, see Figure 3.1, σ = E s (ε ) ε .

(3.4)

In hypo-elastic formulation, the relationship between the increments of stress and strain are defined by the tangential modulus of elasticity E t d ε = E t (ε ) d σ

(3.5)

21


The nonlinear elastic material law can also be formulated in terms of hyperelastic formulation, which assumes the existence of strain energy density function U and the corresponding complementary energy density function U ∗ such that d U ∗ d U σ = and ε = d σ d ε The hyperelastic material model is usually used for rubber-like materials7.

(3.6)

Figure 3.1: Nonlinear elasticity model. Material models for multiaxial states of stress are usually based on generalization of one-dimensional concepts. For example, in a hyperelastic formulation components of stress tensor are computed as σ =

∂U ∂ε

(3.7)

that means σ x =

∂U , ∂ε x

σ y =

∂U , ∂ε y

τ xy =

∂U ∂γ xy

etc.

Figure 3.2: Strain energy density functions U and U ∗ 7

An example is the Mooney-Rivlin material model used for modelling rubber-like materials.

22


For any nonlinear elastic material model, it is possible to define relation between stress and strain increments increments as d σ = DT d ε

(3.8)

where matrix DT is function of strains ε ε. Consequently, a set of equilibrium equations we receive in FEM is nonlinear and must be solved by use of any method (e.g. NR) described above.

3.3 Elastoplastic material model 3.3.1 Yielding criterion criterion

Experiments indicate that linear elastic model is acceptable only within a limited range of stress. As an example, the stress-strain curve from tension test of steel specimen is shown in Figure 3.3. Until the yield stress represented by point A (in the given case σ y = 280 MPa) the deformations are elastic and stress-strain relation may be described as σ = E ε . When the stress level exceeds the yield stress, an elastoplastic constitutive law governs the relationship between increments of stress and strain. Due to lack of information,8 approximate stress-strain curves are usually used in analysis. Bilinear approximation defined by yield stress, modulus of elasticity E and tangential modulus E T T is shown in Figure 3.4. If E T T = 0, material model is E T elastic-perfectly plastic. If E T ≠ 0 material model assumes strain hardening. 450 400 350 300 ] a P 250 M [ 200 σ σ

A

150 100 50 0 0.00

0.02

0.04

0.06

0.08

0.10

0.12

ε ε

Figure 3.3: Typical stress-strain curve for mild steel It should be noted that curves in Figures 3.3 and 3.4 are for tensile behavior. It is usually assumed that similar curves for compressive behavior are applicable if there has been no history of plastic deformation. 8

In a design process, the real material curve is usually unknown, only basic values like yield stress etc. are available. Moreover, the material properties slightly differ by different supplies.

23


Figure 3.4: Elastoplastic model with linear strain hardening The indication of yielding under multiaxial conditions in metals is obtained from experiments usually conducted on cylindrical samples subjected to combined axial load and torque. Experiments suggest that there is i s no significant difference in behavior of metals in tension or compression and no volume change associated with yielding and no effect of mean stress level on yielding can be assumed. In a mathematical description, onset of yielding may be represented by a scalar function termed the yield function F. The yield function is written in a form, which leads to the conditions F < 0

for elastic

and

F = 0

for plastic deformation.

(3.9)

In engineering practice, two following conditions for yielding are most frequently used: ● Von Mises yield criterion F = (σ 1 − σ 1 )2 + (σ 2 − σ 3 ) 2 + (σ 3 − σ 1 )2 − 2σ y = 0

(3.10)

where σ 1, σ 2 and σ 3 are principal stresses. Thus, yield occurs when the effective stress σ eff eff reaches the yield stress value σ y σ eff =

1 (σ 1 − σ 1 )2 + (σ 2 − σ 3 )2 + (σ 3 − σ 1 )2 = σ y . 2

(3.11)

● Tresca yield criterion F = [(σ 1 − σ 2 ) 2 − σ y2 ] [(σ 2 − σ 3 ) 2 − σ y2 ] [(σ 3 − σ 1 ) 2 − σ y2 ] = 0 .

(3.12)

The largest difference between these two classical yield criteria is about 15% for the pure shear stress state. For other stress states is the difference less. Hence, both criteria are frequently frequently considered as as equivalent in engineering engineering practice. Any yield condition that is function of stress tensor components σ σ and material parameters κ κ

24


F (σ , κ ) = 0

(3.13)

defines a yield surface in principal stress space, see Figure 3.5. Stress points that lie inside the yield surface are associated with elastic stress states whereas those that lie on the surface represent plastic stress states. No stress point can be outside the yield surface.

Figure 3.5: Yield surface 3.3.2 Post yielding behavior

The fundamental assumption in describing post-yielding behavior is the decomposition of the total strain increment into an elastic (recoverable) part and a plastic (irreversible) part. For For uniaxial stress state is, according to Figure 3.6 d σ = E T , d ε (3.14) d σ = E d ε e and plastic strain increment is then d ε p = d ε − d ε e =

E − E T E E T

d σ .

Figure 3.2: Decomposition of the total strain increment

(3.15)

25


By analogy, in multiaxial stress state the total strains we decompose into elastic and plastic parts too e ε = ε

+ ε p

(3.16)

In multiaxial cases, subsequent loading after first yield produces further plastic deformation that can result in a modification of the shape and/or position of the yield surface. For a perfectly plastic material, the yield surface remains unchanged during plastic deformation. For a strain hardening material, plastic deformation produces a change in shape and position of the yield surface. This means that initial yield surface is gradually replaced by the subsequent yield surfaces. A modified yield function is adopted which has a form such as F (σ , ε p , K ) = 0 .

(3.17)

This yield function depends on the stresses but also the plastic strains and a hardening parameter K. The way in which the plastic strains modify the yield function is defined by hardening rules:

Figure 3.7: Isotropic hardening. 1. An isotropic hardening law implies that the yield surface increases in size but maintains its original shape under loading conditions. Schematic representation of isotropic hardening for uniaxial and biaxial stress state is shown in Figure 3.7. 2. In kinematic hardening , the original yield surface is translated to a new position in stress space with no change of its shape and size as shown in Figure 3.8. Kinematic hardening has paramount importance in modelling cyclic behavior.

26


3. The combination of the two principal hardening laws leads to a mixed hardening law, where the initial yield surface both expands and translates as a consequence of plastic flow.

Figure 3.8: Kinematic hardening. 3.3.3 Constitutional equations of elastoplastic material

The yield criterion says whether plastic deformation will occur but says nothing about the plastic behavior of a material after onset of plastic deformations. This is defined by so-called flow rule in which is the rate and the direction of plastic strains is related to the stress state and the stress rate. This relation can be expressed as d ε p ij = dλ

∂Q ∂σ ij

(3.17)

d ε p = dλ

∂Q ∂σ

(3.18)

or in matrix form as

where dλ is a scalar value (to be determined) and Q is a scalar valued function of stress components called plastic potential . For metals, the so-called associated flow rule, in which the plastic potential surface coincides with the yield surface, i.e. Q = F

can be adapted to model plastic flow. For some other materials, non-associated flow rule in which Q ≠ F has to be used to model plastic flow adequately. In the following text we will deal with associated flow rule d ε p = dλ

∂ F ∂σ

(3.19)

27


Consider a uniaxial stress state first. The plastic behavior of material is described as dσ = E T d ε

(3.20)

where E T is constant for a bilinear material as obvious from the equations (3.14) and (3.15). In a multiaxial stress state, we can formulate a similar constitutive equation d σ = DT d ε

(3.21)

where tangential material matrix DT can be derived from known stress tensor σ σ, strain tensor ε ε and constitutive matrix D from equation (3.1) in following way: The first step is strain decomposition into elastic d ε e and plastic part dε p d ε = d ε e + d ε p . From constitutive law it follows that

(3.22)

d σ = D d ε e hence d σ = D d ε − d ε p .

(3.24)

From associated plastic flow rule, it follows that d ε p = dλ

∂ F = a d λ ∂σ

(3.25)

where a

 ∂ F  = . σ ∂  

(3.26)

Using equations (3.24) and (3.25) we obtain d σ = D (d ε − a d λ ) .

(3.27)

The stress point must lie in yield surface 9 and hence the following consistency conditions must be fulfilled T

T

∂ F   ∂ F  p d F =    d σ +  p  d ε = 0  ∂σ   ∂ε 

or with respect to equations (3.25) and (3.25) d F = a

T

 ∂ F  d σ +  p   ∂ε 

T a d λ = 0 .

After substituting from equation (3.27) we obtain 9

note that yield surface may change according to hardening rule

(3.29)

28


T  ∂ F  T d F = a D (d ε − a d λ ) +  p  a d λ = 0  ∂ε 

or d F = a T D (d ε − a d λ ) − A d λ = 0

(3.30)

where scalar quantity A is defined as  ∂ F  A = −   ∂ε p 

T a.

(3.31)

Now, we can derive parameter dλ from equation (3.29) d λ =

T a D d ε

(3.32)

A + a T D a

and substituting this expression for d λ into equation (3.28) we finally obtain T  D aa D   d ε . d σ =  D − T    A + a D a 

(3.33)

When compare the last equation with equation (3.21) we can see that DT = D −

Note that material matrix symmetric as well.

D

T Daa D

A + a T D a

is symmetric, i.e.

D

.

T

= D , hence matrix

(3.34) DT

is

3.3.4 Integration of constitutive equations

We have derived that for infinitesimal increments of stress and strain it holds d σ = DT d ε . In FE analysis we need to work with finite increments ∆σ and ∆ε for which is the relation above approximate only, so if we use relation ∆σ = DT ∆ε

for large increments of stress and strain, an error occurs as stress σ + ∆σ in subsequent step will not satisfy constitutive law and consistency condition. Hence, we need to integrate over the increment of pseudo-time ∆σ = ∫ d σ ∆t

where, according to equation (3.24) d σ = D d ε − d ε p .

29


It is important to note that if plastic flow is present, the DT changes during increment ∆t and as a result, the ratio between total and plastic strain changes too. To obtain correct results, various stress increment integration schemes that differs in the degree of approximation have been developed. Frequently used are the following schemes: 3.3.4.1 Generalized trapezoidal rule

Consider that we know stress n. Then at step n+1 σ n +1

σ n ,

total strain

ε n

and plastic strain

= σ n + ∆σ n +1 = σ n + D ∆ε n +1 − ∆ε p n +1

p ε n

at time step (3.35)

∆ε p n +1 = ∆λ [(1 − α )an + α an +1 ]

(3.36)

F n +1 = 0

(3.37)

3.3.4.2 Generalized mid point rule σ n +1

= σ n + ∆σ n +1 = σ n + D ∆ε n +1 − ∆ε p n +1

(3.35)

∆ε p n +1 = ∆λ an +α

(3.36)

F n +α = 0

(3.37)

In both rules, α is a parameter ranging from 0 to 1. For α = 0 we obtain explicit forward Euler integration scheme. Advantage of this algorithm is in its simplicity; disadvantage is that it is conditionally stable only. That means that step increment has to be smaller than some critical value to avoid instability of the solution. For α = 1 we obtain implicit backward Euler integration scheme ∆σ n +1 = D ∆ε n +1 − ∆ε p n +1 ∆ε p n +1 = ∆λ an +1 F n +1 = F (σ n +1 , ε n +1 ) = 0

30


It is obvious that in difference with forward scheme, we deal with values defined at the end of the increment, which are unknown at start of it. Hence, the procedure is of an iterative nature. This means that at beginning of the increment, the trial stress is estimated by assuming elastic deformation and computed values are then checked whether consistency condition and constitutional equation are satisfied. If not, the process is repeated with improved values until the conditions are satisfied. 3.3.5 Numerical procedures

The tangential material matrix DT is used to form a tangential stiffness matrix K T. When the tangential stiffness matrix is defined, the displacement increment is obtained for a known load increment K T ∆d

= ∆F

(3.22)

As load and displacement increments are final, not infinitesimal, displacements obtained by solution of this set of linear algebraic equation will be approximate only. That means, conditions of equilibrium of internal and external nodal forces will not be satisfied and iterative process is necessary. Any of methods mentioned above may be used. The problem that arises now is the fundamental problem in computational elastoplasticity - not only equilibrium equations but also constitutive equations of material must be satisfied. That means that within the each equilibrium iteration step check of stress state and iterations to find elastic and plastic part of strains at every integration point must be included. The iteration process continues until both, equilibrium conditions and constitutive equations are satisfied simultaneously. The converged solution at the end of load increment is then used at the start of new load increment.


31

REFERENCES [1]

Hinton, E: Introduction to Nonlinear Finite Element Analysis. NAFEMS, Glasgow 1992

[2]

Crisfield, M. A.: Non-linear Finite Element Analysis of Solids and Structures. John Wiley & Sons 1991

[3]

Bittnar, Z., Šejnoha, J.: Numerické metody mechaniky (Numerical methods of mechanics). ČVUT, Praha 1992

[4]

Okrouhlík, M., Höshl, C., Plešek, J., Pták, S., Nadrchal, J.: Mechanika poddajných těles, numerická matematika a superpočítače (Mechanic of solids, numerical mathematics and supercomputers). Czech Academy of Science, Prague 1997.

[5]

Okrouhlík, M.: Implementation of Nonlinear Continuum Mechanics in Finite Element Codes. Institute of Thermodynamics, Prague 1995.

[6]

Hinton, E., Ezatt, M., H.: Fundamental Tests for Two and Three Dimensional, Small Strain, Elastoplastic Finite Element Analysis. NAFEMS, Glasgow 1987.

[7]

Electronic documentation of program COSMOS/M, version 2.95. SRAC, Los Angeles 2005, www,cosmosm.com.

Falzon, B., G., Hitchings, D.: An Introduction to Modeling Buckling and Collapse. NAFEMS, Glasgow 2006 [9] Eurocode 3: Design of steel structures - Part 1-1: General rules and rules for buildings. CEN, Brussels 2005. [8]

32


APPENDIX

33


A 1. NONLINEAR ANALYSES WITH COSMOS/M PROGRAM Properties of finite elements

For nonlinear analyses it is necessary to compute tangential stiffness matrices of individual elements. This possibility is available only for some elements according to Table 1. Table 1: Properties of finite elements Element Name TRUSS2D TRUSS3D BEAM2D BEAM3D SPRING

Element Description

Plane truss Space truss Plane beam Space beam Axial an/or torsional spring 4 to 8 node plane stress, plane PLANE2D strain, axisymmetric element 3 to 6 node plane stress, plane TRIANG strain, axisymmetric triangular element 3-node triangular thin shell SHELL3 4-node quadrilateral thin shell SHELL4 SHELL3T 3-node triangular thick shell 6-node triangular thin shell SHELL6 6-node triangular thick shell SHELL6 SHELL4T 4-node quadrilateral thick shell SHELL3L 3-node triangular composite shell 4-node quadrilateral composite SHELL4L shell 8 to 20 node hexahedron SOLID TETRA4 4-node tetrahedron TETRA10 10-node tetrahedron GAP Contact element with friction MASS Concentrated mass

GNL Le Ne Pl Ve He

Geometric nonlinearities Linear elastic material Nonlinearly elastic material Elastic-plastic material Viscoelastic material Hyper elastic material

GNL

Material model Le Ne Pl Ve He

• • • • •

• • • • •

• • • • •

• • • •

• • • •

•

•

•

•

•

•

•

•

•

•

•

•

• • • • • • •

• • •

• • •

• • •

• • •

• • •

• •

•

•

•

•

•

•

• • •

• • •

• • •

• • •

• • •

• • •

34


Nonlinear properties of individual elements should be defined in element options. For example, steps needed for specification of SOLID element properties are shown in Figure 4 and subsequent figures.

Figure 4: Definition of finite element group

Figure 5: Element specification

Figure 6: Option 1 of the SOLID element


Figure 7: Definition of material model

Figure 8: Definition of geometric nonlinearity

35

36


Figure 9: Listing of element groups – element group 1 defines a linear element SOLID; element group 2 defines geometric and material nonlinear element.

Nonlinear analysis setup

Parameters of nonlinear analysis are specified from MENU ANALYSIS / NONLINEAR, see Figure 10. Individual items in this menu have following use: 1. Solution Control – definition of method of solution of system of nonlinear equations. Possible values are in Figure 11. 2. Integration Options – definition of method of solution for dynamic problems, see Figure 12. 3. Auto Step Option – if is set “ON”, the program automatically selects optimal size of time increments according to convergence in previous time steps, see Figure 13. 4. Base Motion Parameter and Damping Coefficient – applicable for dynamic problems only. 5. Print Options – specification what will be written in output file name of problem .OUT. The advisable is to oppress printing of large files that are difficult to read. Recommended is to use setup according to Figure 14. 6. Plot Options – specification for which steps will program store results. Normally, if nothing is specified, results are stored for the last successfully accomplished step only. If there is need to store more results to plot and list, these steps have to be specified as shown Figure 15.


37

7. Response Options – specification for which nodes will program store results to plot nodal values as displacements and reactions, see Figure 16. This specification is useful if we do not need results for all nodes. If a plot option is specified for all steps, this option is not necessary.

Figure 10: Nonlinear Analysis Menu

38


Figure 11: Solution Control menu

Figure 12: Integration method for solution of dynamic problems specification

Figure 13: Auto step definition – minimum acceptable time increment is set to 10-6, maximum allowed time increment is 1,5 and maximum number of trials to select proper time increment is limited to 5.


39

Figure 14: Output file specification

Figure 15: Specification that results will be stored for steps 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 and then for steps 12, 14, 16, 18 and 20

40


Figure 16: Specification of nodes to plot graphs

8. Nonlinear Analysis Options – detailed specification of analysis. Default setup, which is appropriate for most nonlinear static problems is in Figure 17. 9. Run Nonlinear Analysis – start of computations. For other items of the NONLINEAR menu refer to the program COSMOS/M documentation.


Figure 17: Specification of nonlinear analysis – default values of parameters

41

42


Time curves and time parameters

At least one time curve has to be specified for nonlinear analysis. Time curves can be defined from menu Loads BC , see Figure 18 and subsequent figures.

or from icons in left panel

Figure 18: Start of definition of time curves

Figure 19: Time curve definition – 1 st step


43

Figure 20: Time curve definition – 2 nd step (in this example, definition of the time curve number 1 is prepared)

Figure 21: Time curve definition – 3 rd step –coordinates of points will start with first point 10

Figure 22: Definition of points of the time curve 10

start point number higher than 1 we use in a case that curve is defined and we want to add some points

44


To plot a time curve, commands from menu DISPLAY can be selected as shown in figure and subsequent figures.

Figure 23: Commands to prepare a plot of a time curve


Figure 24: Commands to plot a time curve

Figure 25: Graph of the time curve number 1 as defined in Figure 22

45

46


Note that after definition, the time curve became active. This means that boundary conditions defined subsequently will be associated with it. To deactivate a time curve or to switch activation to another one, the commands according to Figure 26 should be used.

Figure 26: Making a time curve active (or not active)


47

Restart possibilities

Nonlinear analyses are time demanding. This is a reason that there is a possibility to continue in computation after stop of program. For this so called restart of computation the parameter RESTART, see Figure 27 have to be set to 1. If set to 0, computation will start from beginning.

Figure 27: Restart options

48


A 2. EXAMPLES OF COSMOS/M USE Compressed slender beam

Problem description: Slender elastic beam is compressed by force F and laterally loaded by a small force P . The forces increase according to time curves11 TC 1 and TC 2. Force P = 200 N is associated with the curve TC 1 and force F = 100000 N with the curve TC 212. Dimensions of the beam are L = 1600 mm, a = 80 mm and t = 3 mm. Modulus of elasticity is E = 2,1⋅105 MPa.

Fig. 1: Compressed beam 1.8

TC 1

1.6

TC 2

1.4 1.2 1 0.8 0.6 0.4 0.2 0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Fig. 2: Time curves

11

The time is now only a parameter controlling the loading, because the problem is static. Association with a time curve means that at every time the value of force is multiplied by instantaneous value of time curve ordinate.

12

49


Modelling hints: 1. The BEAM2D element group is used to model beam. As the problem is geometrically non-linear, finer finite element mesh must be used than for a linear problem. Geometric nonlinearities are taken into account if option 6 for the element group BEAM2D is set to 1 (large displacement formulation). 2. It is useful to perform linear and buckling analyses before the nonlinear analysis. For comparison, the Euler critical force is given by F cr =

π 2 E I z

(2 L) 2

.

For the dimensions and modulus of elasticity given above is F cr = 185088,7 N. 3. Load should be divided into small increments and NR, or MNR iterations of equilibrium should be used for each incremental step. 4. Finer load increments are necessary when F → F cr . 5. Restart from the last successfully accomplished load increment can be performed. Commands: ( NOTE: C* means a comment only – row with C* is not processed by the program ) C* 1.

Geometry creation

VIEW,0,0,1,0 PT,1,0,0,0 CREXTR,1,1,1,Y,1600 CREXTR,2,2,1,Y,1600 C* 2. Defining material properties MPROP,1,EX,2.1E5 C* 3.

Defining elemnt group and beam section

EGROUP,1,BEAM2D,0,0,0,0,0,1,0,0 BMSECDEF,1,1,4,1,9,80,80,3,3,0,0,0,0,0 C* 4.

Creating FE mesh

M_CR,1,2,1,2,5,1 NMERGE,1,12,1,1,0,0,0 NCOMPRESS,1,12 C* 5.

Defining boundary conditions

C*

5.1 prescribed displacements

DPT,1,UX,0,3,2, DPT,1,UY,0,1,1,

50


C*

5.2 defining force P

CURDEF,TIME,1,1,0,0,.1,1,10,1 FPT,2,FX,200,2,1 C* 5.3 force F CURDEF,TIME,2,1,0,0,10,10,10 FPT,3,FY,-1E5,3,1 C* 6. Adjusting nonlinear analysis A_NONLINEAR,S,1,1,20,0.001,0,N,0,0,1E+010,0.001,0.01,0,1,0,0 NL_PLOT,1,100,1,0 C* 7. specification of load increments TIMES,0,1.7,.1 C* 7. nonlinear analysis R_NONLIN C* restart from the last step RESTART,1 TIMES,1.7,1.8,0.02 R_NONLIN C* RESTART,1 TIMES,1.8,1.82,0.005 R_NONLIN

Results: Typical nonlinear response is obvious from the graph (Figure 3) showing beam deflection at the middle of span versus load parameter (time). 25

20

] 15 m m [ X U

10

5

0 0

0,2

0,4

0,6

0,8

1

1,2

1,4

TIME

Fig. 3: Deflection versus time

1,6

1,8

2

51


Thick-walled pipe subjected to internal pressure

Problem description: The long cylindrical pressure vessel is subjected to internal pressure. Both ends are assumed fixed. Compute stress distribution during pressure test and during normal operation. Assume small plastic deformations according to bilinear material model if yield stress is σ Y = 260 MPa, modulus of elasticity is E = 2,1.105 MPa and tangent modulus is E T = 100 MPa.

Fig. 4: Thick walled pipe

300

250

200 s s e r 150 t s

100

50

0 0

0.005

0.01

0.015

0.02

0.025

0.03

strain

Fig. 5: Bilinear material mode Time curve reflects load history. During test is internal pressure increased to 1,2 × nominal value, which is 120 MPa.

52


Fig. 6: Load history – time curve Modelling hints: 1. Due to axial symmetry, only a part of the cylinder longitudinal section may be modelled. 2. 8 node PLANE2D elements are used with option 5 set to “Von Mises kinematic hardening”. 3. Load incrementation with NR iterations should be used. 4. Restart is used to change load (time) increment.

Fig. 7: FE mesh and boundary conditions Commands: C* GEOMETRY CREATION VIEW,0,0,1,0 PT,1,45,0,0 CREXTR,1,1,1,X,35


53

SFEXTR,1,1,1,Y,10 C* DEFINING ELEMENT GROUP EGROUP,1,PLANE2D,0,2,1,0,2,0,0,0 C* DEFINING MATERIAL PROPERTIES MPROP,1,EX,2.1E5 MPROP,1,NUXY,.3 MPROP,1,ETAN,100 MPROP,1,SIGYLD,260 C* MESHING M_SF,1,1,1,8,8,2,1,1 C* DEFINING BOUNDARY CONDITIONS DCR,2,UY,0,1,1, CURDEF,TIME,1,1,0,0,1,1.2,1.1,1.2,2,0,3,1 PCR,3,120,3,1,120,4 C* ADJUSTING NONLINEAR ANALYSIS NL_CONTROL,0,1 NL_PRINT,0,0,0,0,0,0,0 NL_PLOT,1,50,1,0 C* TEST SIMULATON TIMES,0,1,.2 R_NONLIN C* RESTART,1 TIMES,1,1.1,.1 R_NONLIN C* UNLOADING TIMES,1.1,2,0.09 R_NONLIN C* C* NOMINAL PRESSURE LOADING TIMES,2,3,0.2 R_NONLIN C*

NOTE: C* means comment only – row with C* is not processed by the program

54


Results: Note that 1. Linear analysis gives non-realistic values of stress (maximal von Mises stress is about 304 MPa). 2. Plastic deformations occur during test loading. 3. After test, residual stress is present in the cylinder. This redistribution of stress is useful and desired. 4. Due to stress redistribution, only stress within elastic range is present during nominal loading after the test.


55

Deformation of a thin-walled tank

Problem description: Thin-walled pressure vessel – tank is subjected to internal pressure. Compute stress distribution according to linear and nonlinear static analysis. In the nonlinear analysis, both geometrical and material nonlinearities should be taken into account. Material is elastoplastic, the same as in previous example (Figure 5). Nominal value of internal pressure is 0,6 MPa. Symmetric half of the tank is shown in Figure 8. Diameter of the manhole cover (thickness of 25 mm) is 600 mm.

Fig. 8: Thin-walled tank Modeling hints: 1. SHELL4T elements with options for large displacement, small plastic strain, von Mises yield criterion and kinematic hardening are used. 2. Due to symmetry, only one eighth of the tank may be modelled. 3. Stresses are calculated in element co-ordinate systems. To receive compatible co-ordinate systems curves are reoriented. 4. Element co-ordinate system of SHELL4T element is defined according to Figure 9. The element x-axis goes from the first node to the second. The y-axis lies in the plane defined by the first three nodes perpendicular to the x-axis toward the fourth node. The z -axis completes a right-hand Cartesian system. When meshing surfaces, element x-axes are parallel to the first surface curve marked by asterisks in Figure 10. FE mesh is shown in Figure 11. 5. The automatic time stepping and restart options are used for nonlinear FE analysis.

56


z

TOP FACE y

x

3 2 BOTTOM FACE 1

Fig. 9: Element co-ordinate system (ECS)

Fig. 10: Orientation of shell elements

Fig. 11: Finite element mesh


57

Figure 28: Example of inappropriate orientation of elements – orientation of elements on surface 1 differs from orientation of other elements 13

Figure 29: Orientation of surfaces identification

13

top faces of shell elements are red colour, bottom faces are blue

58


Commands: C* C* 1. DEFINING GEOMETRY C* C* 1.1. creating points and curves C* VIEW,0,0,1,0 PT,1,0,0,0 PT,2,2900,0,0 PLANE,Z,0,1 CRPCIRCLE,1,2,1,2900,-80,1 PT,4,0,1700,0 C* CREXTR,4,4,1,X,2450 C* 1.1.1 intersecting curves CRINTCC,1,2,2,1,2,5E-005 C* 1.1.2 deleting useless curves CRDEL,4,2,2 C* 1.1.3 creating fillet etc. CRFILLET,4,1,3,400,1,0,1E-006 PT,10,-20,300,0 CREXTR,10,10,1,X,100 CRINTCC,1,5,5,1,0,5E-005 CRDEL,1,1,1 CRDEL,5,5,1 CREXTR,12,12,1,Y,-300 C*

1.1.4

C*

changing orientation of curves 3 and 7 in order to receive desired ECS orientation

CRREPAR,7,3,4 VIEW,-1,2,3,0 CRCOMPRESS,1,7 C*


C* 1.2 creating surfaces C* SFSWEEP,1,4,1,X,90,2 C* C* 2. FE MESH DEFINITION C* 2.1 defining element type and real constant sets EGROUP,1,SHELL4T,1,0,0,1,2,1,0,0 RCONST,1,1,1,6,6,0,0,0,0,0 RCONST,1,2,1,6,25,0,0,0,0,0 C* 2.2 defining material properties MPROP,1,EX,2.1E5,NUXY,.3,SIGYLD,260 MPROP,1,ETAN,100 C* 2.3 meshing ACTSET,RC,1 M_SF,5,6,1,4,10,6,1.8,1 VIEW,-1.5,2,2,0 ACTSET,RC,2 MA_NUSF,7,8,1,4,4,6,0 MASFCH,7,8,1,Q,4,1,0.4,1 ACTSET,RC,1 M_SF,3,4,1,4,3,6,1,1 M_SF,1,2,1,4,10,6,2,1 NMERGE,1,400,1,1,0,0,0 NCOMPRESS,1,400 NCOMPRESS,1,330 ECOMPRESS,1,315 C* 2.4 checking real constants assignment ACTECLR,1,RC,1 SHADE,1,4,1 EPLOT,1,ELMAX,1 ACTECLR,0 C* C* 3. DEFINING BOUNDARY CONDITIONS C*

59

60


C* 3.1 displacements SELPIC,CR,20,16,12,6,0 DCR,1,SY,0,CRMAX,1, INITSEL,CR,1,1 SELPIC,CR,4,3,2,1,0 DCR,1,SZ,0,CRMAX,1, INITSEL,ALL,1,1 DCR,8,SX,0,10,2, CLS,1 C* 3.2 pressure CURDEF,TIME,1,1,0,0,10,10 PSF,1,0.6,8,1,0.6,0.6,4 C* C* 4. ANALYSIS C* 4.1 preliminary linear static analysis C* R_STAT C* C* 4.2 nonlinear analysis adjustment NL_AUTOSTEP,1,.001,.1,5 NL_PLOT,1,20,1,0 NL_PRINT,0,0,0,0,0,0,0 A_NONLINEAR,S,1,1,20,0.001,0,N,0,1,1E+010,0.001,0.01,0,1,0,0 C* 4.3 running nonlinear analysis TIMES,0,1,.1 R_NONLIN C* C* 4.4 continuing in NL analysis RESTART,1 TIMES,1,1.5,.1 R_NONLIN C* C* 5. POSTPROCESSING C* plotting stress, displacements and xy graphs. C*

NOTE: C* means comment only – row with C* is not processed by the program


61

Results:

1. Von Mises stress distribution for time step 16 ( t = 1,5) on top faces of shell elements is in . Von Mises stress distribution on bottom faces is in . 2. There are large differences between results of linear and nonlinear analysis. In linear analysis, the von Mises stress is above the yield stress 260 MPa. Difference between linear and nonlinear analysis results is a consequence of different solution methods. While in the linear analysis is stress considered as elastic and geometry of the tank is considered unchanged, in the nonlinear analysis changes of the tank shape are taken into account together with material nonlinearities. 3. Onset of plastic deformation is obvious from the graph in Fig. 12.

Figure 30: Von Mises stress distribution on top faces at t = 1,5.

FEA of Nonlinear Problems 2011 Ivanco

Recommend Documents