Stress and Strain Calculation

The stress and strain for a shell element can be written in vector notation. Each component is a stress or strain feature of the element.

The generalized strain ε can be written as:(1)
{ ε } = { e x , e y , e x y , k x , k y , k x y }
Where,
e ij MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyzamaaBa aaleaacaWGPbGaamOAaaqabaaaaa@38E9@
Membrane strain
χ i j
Bending strain or curvature
The generalized stress can be written as:(2)
{ } = { N x , N y , N x y , M x , M y , M x y }
Where,
N x = t / 2 t / 2 σ x d z M x = t / 2 t / 2 σ x z d z
N y = t / 2 t / 2 σ y d z M y = t / 2 t / 2 σ y z d z
N x y = t / 2 t / 2 σ x y d z M x y = t / 2 t / 2 σ x y z d z
N y z = t / 2 t / 2 σ y z d z N x z = t / 2 t / 2 σ x z d z

Isotropic Linear Elastic Stress Calculation

The stress for an isotropic linear elastic shell for each time increment is computed using:(3)
{ e l ( t + Δ t ) } = { ( t ) } + L { ε ˙ } Δ t
Where,(4)
L = [ L m 0 0 L b ]
(5)
L m = [ E t 1 v 2 v E t 1 v 2 0 v E t 1 v 2 E t 1 v 2 0 0 0 E t 1 + v ]
(6)
L b = [ E t 3 12 ( 1 v 2 ) v E t 3 12 ( 1 v 2 ) 0 v E t 3 12 ( 1 v 2 ) E t 3 12 ( 1 v 2 ) 0 0 0 E t 3 12 ( 1 + v ) ]
E MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamiDaaaa@36EF@
Young's or Elastic modulus
υ
Poisson's ratio
t MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamiDaaaa@36EF@
Shell thickness

Isotropic Linear Elastic-Plastic Stress Calculation

An incremental step-by-step method is usually used to resolve the nonlinear problems due to elasto-plastic material behavior. The problem is presented by the resolution of the following equation:(7)
σ ˙ = C : ( ε ˙ ε ˙ p )
(8)
f ( σ , σ y ) = σ e q 2 σ y 2 = 0 ; σ ˙ y = H ε ˙ p
(9)
f ˙ = 0
and(10)
ε ˙ p = f σ λ ˙
f is the yield surface function for plasticity for associative hardening. The equivalent stress σ e q may be expressed in form:(11)
σ e q 2 = { σ } t [ A ] { σ }

with { σ } = { σ x x σ y y σ x y } and [ A ] = [ 1 1 2 0 1 2 1 0 0 0 3 ] for von Mises criteria.

The normality law (Stress and Strain Calculation, 式 10) for associated plasticity is written as:(12)
{ ε ˙ p } = f ( σ ) λ ˙ = 2 [ A ] { σ } λ ˙ = ε ˙ p σ y [ A ] { σ }
Where,
ε ˙ p
Equivalent plastic deformation
Stress and Strain Calculation, 式 7 is written in an incremental form:(13)
{ σ } n + 1 = { σ } n + { d σ } = { σ } n + [ C ] ( { d ε } { d ε p } ) = { σ } d ε p σ y n + 1 [ C ] [ A ] { σ } n + 1
Where, { σ * } represents stress components obtained by an elastic increment and [ C ] MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaWaamWaaeaaca WGdbaacaGLBbGaayzxaaaaaa@38B0@ the elastic matrix in plane stress. The equations in Stress and Strain Calculation, 式 7 to 式 13 lead to obtain the nonlinear equation:(14)
f ( d ε p ) = 0

that can be resolved by an iterative algorithm as Newton-Raphson method.

To determine the elastic-plastic state of a shell element, a number of steps have to be performed to check for yielding and defining a plasticity relationship. Stress-strain and force-displacement curves for a particular ductile material are shown in 図 1.


図 1. Material Curve
The steps involved in the stress calculation are:
  1. Strain calculation at integration point z
    The overall strain on an element due to both membrane and bending forces is:(15)
    ε x = e x z χ x
    (16)
    ε y = e y z χ y
    (17)
    ε x y = e x y z χ x y
    (18)
    { ε } = { ε x , ε y , ε x y }
  2. Elastic stress calculation
    The stress is defined as:(19)
    { σ } = { σ x , σ y , σ x y }
    It is calculated using explicit time integration and the strain rate:(20)
    { σ e l ( t + Δ t ) } = { σ ( t ) } + L { ε ˙ } Δ t
    L = [ E 1 v 2 v E 1 v 2 0 v E 1 v 2 E 1 v 2 0 0 0 E 1 + v ]
    The two shear stresses acting across the thickness of the element are calculated by:(21)
    σ y z e l ( t + Δ t ) = σ y z ( t ) + α E 1 + v e ˙ y z Δ t σ x z e l ( t + Δ t ) = σ x z ( t ) + α E 1 + v e ˙ x z Δ t

    Where, α is the shear factor. Default is Reissner's value of 5/6.

  3. von Mises yield criterion
    The von Mises yield criterion for shell elements is defined as:(22)
    σ v m 2 = σ x 2 + σ y 2 σ x σ y + 3 σ x y 2
    For type 2 simple elastic-plastic material, the yield stress is calculated using:(23)
    σ y i e l d ( t ) = a + b ε p n ( t )
  4. Plasticity Check
    The element's state of stress must be checked to see if it has yielded. These values are compared with the von Mises and Yield stresses calculated in the previous step. If the von Mises stress is greater than the yield stress, then the material will be said to be in the plastic range of the stress-strain curve.


    図 2. Plasticity Check
  5. Compute plastically admissible stresses

    If the state of stress of the element is in the plastic region, there are two different analyses that can be used as described in the next paragraph. The scheme used is defined in the shell property set, card 2 of the input.

  6. Compute thickness change
    The necking of the shells undergoing large strains in hardening phase can be taken into account by computing normal strain ε z z in an incremental process. The incompressibility hypothesis in plasticity gives:(24)
    d ε z z p = ( d ε x x p + d ε y y p )
    Where, the components of membrane strain d ε x x p and d ε y y p are computed by 式 12 as:(25)
    { d ε x x p d ε y y p } = d ε p σ y [ A ] 2 x 2 { σ x x σ y y } MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake GabaaVqmaaceaabiqaaWledaGacaabceqabaaVqeaacaWGKbGaeqyT du2aa0baaSqaaiaadIhacaWG4baabaGaamiCaaaaaOqaaiaadsgacq aH1oqzdaqhaaWcbaGaamyEaiaadMhaaeaacaWGWbaaaaaakiaaw2ha aaGaay5EaaGaeyypa0ZaaSaaaeaacaWGKbGaeqyTdu2aaWbaaSqabe aacaWGWbaaaaGcbaGaeq4Wdm3aaSbaaSqaaiaadMhaaeqaaaaakmaa dmaabaGaamyqaaGaay5waiaaw2faamaaBaaaleaacaaIYaGaamiEai aaikdaaeqaaOWaaiqaaeGabaaVqmaaciaaeiqabeaa8craaiabeo8a ZnaaBaaaleaacaWG4bGaamiEaaqabaaakeaacqaHdpWCdaWgaaWcba GaamyEaiaadMhaaeqaaaaakiaaw2haaaGaay5Eaaaaaa@61E0@
  7. The plan stress condition d σ z z = 0 allows to resolve for d ε z z :(26)
    d ε z z = υ 1 υ ( d ε x x + d ε y y ) + 1 2 υ 1 υ d ε z z p

Plastically Admissible Stresses

Radial return
Iplas=2
When the shell plane stress plasticity flag is set to 0 on card 1 of the shell property type definition, a radial return plasticity analysis is performed. Thus, Step 5 of the stress computation is:
The hardening parameter is calculated using the material stress-strain curve:(27)
σ y i e l d ( t ) = a + b ε p n ( t ) ε ˙ p Δ t = σ v m σ y i e l d E
Where, ε p is the plastic strain rate.
The plastic strain, or hardening parameter, is found by explicit time integration:(28)
ε p ( t + Δ t ) = ε p ( t ) + ε ˙ p Δ t
Finally, the plastic stress is found by the method of radial return. In case of plane stress this method is approximated because it cannot verify simultaneously the plane stress condition and the flow rule. The following return gives a plane stress state:(29)
σ i j p a = σ y i e l d σ v m σ i j e l
Iterative algorithm
Iplas=1
If flag 1 is used on card 1 of the shell property type definition, an incremental method is used. Step 5 is performed using the incremental method described by Mendelson. 1 It has been extended to plane stress situations. This method is more computationally expensive, but provides high accuracy on stress distribution, especially when one is interested in residual stress or elastic return. This method is also recommended when variable thickness is being used. After some calculations, the plastic stresses are defined as:(30)
σ x p a = σ x e l + σ y e l 2 ( 1 + Δ r 1 v ) + σ x e l σ y e l 2 ( 1 + 3 Δ r 1 + v )
(31)
σ y p a = σ x e l + σ y e l 2 ( 1 + Δ r 1 v ) σ x e l σ y e l 2 ( 1 + 3 Δ r 1 + v )
(32)
σ x y p a = σ x y e l 1 + 3 Δ r 1 + v
Where,(33)
Δ r = E Δ ε p 2 σ y i e l d ( t + Δ t )
The value of Δ ε P must be computed to determine the state of plastic stress. This is done by an iterative method. To calculate the value of Δ ε P , the von Mises yield criterion for the case of plane stress is introduced:(34)
σ x 2 + σ y 2 σ x σ y + 3 σ x y 2 = σ y i e l d 2 ( t + Δ t )
and the values of σ x , σ y , σ xy and σ yield are replaced by their expression as a function of Δ ε P (Hourglass Modes, 式 1 to Hourglass Modes, 式 4), with for example:(35)
σ y i e l d ( t + Δ t ) = a + b ε p n ( t + Δ t )
and:(36)
ε p ( t + Δ t ) = ε p ( t ) + Δ ε p
The nonlinear equation 式 34 is solved iteratively for Δ ε P by Newton's method using three iterations. This is sufficient to obtain Δ ε P accurately.

Plastic Plane Stress with Hill's Criterion

In the case of Hill's orthotropic criterion, the equivalent stress is given by:(37)
σ e q 2 = A 1 σ x x 2 + A 2 σ y y 2 A 3 σ x x σ y y + A 12 σ x y 2

with [ A ] = [ A 1 A 3 2 0 A 2 0 s y m A 12 ]

式 13 is then written as:(38)
[ B ] { σ } n + 1 = { σ }
[ B ] = [ 1 + ( A 1 A 3 2 υ ) d R ( A 2 υ A 3 2 ) d R 0 ( A 1 υ A 3 2 ) d R 1 + ( A 2 A 3 2 υ ) d R 0 0 0 1 + 1 υ 2 A 12 d R ] ; d R = d ε p σ y E 1 υ 2
Changing the stress variables to { σ ¯ } :(39)
{ σ ¯ } = [ Q ] { σ }

with:

[ Q ] = [ 1 A 1 A 2 + C 2 ( A 2 υ A 3 2 ) 0 A 1 A 2 + C 2 ( A 1 υ A 3 2 ) 1 0 0 0 1 ] ; C = ( 1 υ 2 ) ( A 1 A 2 ) 2 + [ A 3 ( A 1 + A 2 ) υ ] 2

The matrix [ B ¯ ] = [ Q ] [ B ] [ Q ] 1 is diagonal:(40)
[ B ¯ ] = [ 1 + d R ( A 2 A 3 υ 2 C J Q ) 0 0 0 1 + d R ( A 1 A 3 υ 2 + C J Q ) 0 0 0 1 + d R . A 12 ( 1 υ 2 ) ]
Where, J Q = 1 + ( A 1 A 2 + C ) 2 4 ( A 1 υ A 3 2 ) ( A 2 υ A 3 2 ) is the Jacobian of [ Q MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyyaaaa@36DC@ ]. 式 39 is now written as:(41)
[ B ¯ ] { σ ¯ } n 1 = { σ ¯ }
This will enable to give explicitly the expression of the yield surface 式 14:(42)
f n + 1 = ( σ e q n + 1 ) 2 ( σ y n + 1 ) 2 = { σ ¯ } t [ B ¯ ] t [ A ¯ ] [ B ¯ ] 1 { σ ¯ } ( σ y n + 1 ) 2 = A ¯ 11 B ¯ 1 2 σ ¯ x x 2 + A ¯ 22 B ¯ 2 2 σ ¯ y y 2 + A 12 B ¯ 3 2 σ ¯ x y 2 + A ¯ 12 B ¯ 1 B ¯ 2 σ ¯ x x σ ¯ y y ( σ y n + 1 ) 2

With [ A ¯ ] 2 × 2 = [ Q ] 2 × 2 t [ A ] 2 × 2 [ Q ] 2 × 2 1 .

The derivative of f n + 1 is carried out in order to use the Newton-Raphson method:(43)
Δ ε n + 1 p = Δ ε n p f n + 1 f n + 1
1 Mendelson A., 「Plasticity: Theory and Application」, MacMillan Co., New York, 1968.