Skip to main content

Stress In 2D Elements

Strain Definitions

The normal definitions of strain used are as follows

εxx=uxεxx=vyεxx=wzγxy=uy+vxγyz=vz+wyγzx=wx+uz\begin{aligned} \varepsilon_{xx} &= \frac{\partial u}{\partial x} & \varepsilon_{xx} &= \frac{\partial v}{\partial y} & \varepsilon_{xx} &= \frac{\partial w}{\partial z}\\ \gamma_{xy} &= \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} & \gamma_{yz} &= \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} & \gamma_{zx} &= \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \end{aligned}

An alternative definition which fits more neatly in tensor form is

εxx=uxεxx=vyεxx=wzεxy=12(uy+vx)εyz=12(vz+wy)εzx=12(wx+uz)\begin{aligned} \varepsilon_{xx} &= \frac{\partial u}{\partial x} & \varepsilon_{xx} &= \frac{\partial v}{\partial y} & \varepsilon_{xx} &= \frac{\partial w}{\partial z}\\ \varepsilon_{xy} &= \frac{1}{2}\left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) & \varepsilon_{yz} &= \frac{1}{2}\left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right) & \varepsilon_{zx} &= \frac{1}{2}\left( \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \right) \end{aligned}

with the strain tensor defined as

ε=[εxxεxyεxzεyxεyyεyzεzxεzxεzz]\mathbf{\varepsilon} = \begin{bmatrix} \varepsilon_{xx} & \varepsilon_{xy} & \varepsilon_{xz} \\ \varepsilon_{yx} & \varepsilon_{yy} & \varepsilon_{yz} \\ \varepsilon_{zx} & \varepsilon_{zx} & \varepsilon_{zz} \\ \end{bmatrix}

The calculation of principal strains ε1,ε2,ε3\varepsilon_{1},\varepsilon_{2},\varepsilon_{3} follows from

E3(εxx+εyy+εzz)E2+(εxxεyy+εyyεzz+εzzεxxεxy2εyz2εzx2)E(εxxεyyεzz+2εxyεyzεzxεxxεyz2εyyεzx2εzzεxy2)=0\begin{aligned}E^{3} - &\left( \varepsilon_{xx} + \varepsilon_{yy} + \varepsilon_{zz} \right)E^{2} + \\ &\left( \varepsilon_{xx}\varepsilon_{yy} + \varepsilon_{yy}\varepsilon_{zz} + \varepsilon_{zz}\varepsilon_{xx} - {\varepsilon_{xy}}^{2} - {\varepsilon_{yz}}^{2} - {\varepsilon_{zx}}^{2} \right)E - \\ &\left( \varepsilon_{xx}\varepsilon_{yy}\varepsilon_{zz} + 2\varepsilon_{xy}\varepsilon_{yz}\varepsilon_{zx} - \varepsilon_{xx}{\varepsilon_{yz}}^{2} - \varepsilon_{yy}{\varepsilon_{zx}}^{2} - \varepsilon_{zz}{\varepsilon_{xy}}^{2} \right) = 0\end{aligned}

The maximum shear strain is calculated from the principal strain as

εmax shear=12(ε1ε3)\varepsilon_{\text{max shear}} = \frac{1}{2}\left( \varepsilon_{1} - \varepsilon_{3} \right)

or

γmax shear=(ε1ε3)\gamma_{\text{max shear}} = \left( \varepsilon_{1} - \varepsilon_{3} \right)

In a similar way to the definitions of average and von Mises stress a volumetric and effective strain can be calculated as

εav=13(εxx+εyy+εzz)εvM=12(1+ν)[(εxxεyy)2+(εyyεzz)2+(εzzεxx)2+6(εxy2+εyz2+εzx2)]\begin{aligned}\varepsilon_{av} &= \frac{1}{3}\left( \varepsilon_{xx} + \varepsilon_{yy} + \varepsilon_{zz} \right)\\ \varepsilon_{vM} &= \frac{1}{\sqrt{2}(1 + \nu)}\left\lbrack \left( \varepsilon_{xx} - \varepsilon_{yy} \right)^{2} + \left( \varepsilon_{yy} - \varepsilon_{zz} \right)^{2} + \left( \varepsilon_{zz} - \varepsilon_{xx} \right)^{2} + 6\left( {\varepsilon_{xy}}^{2} + {\varepsilon_{yz}}^{2} + {\varepsilon_{zx}}^{2} \right) \right\rbrack\end{aligned}

Stress Definitions

Stress can be considered as a tensor quantity whose components can be represented in matrix form as

σ=[σxxσxyσxzσyxσyyσyzσzxσzxσzz]\boldsymbol{\sigma} = \begin{bmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{zx} & \sigma_{zx} & \sigma_{zz} \\ \end{bmatrix}

where each term corresponds to a force per unit area. The following notation for the stress components is common

σx=σxxσy=σyyσz=σzzτxy=σxyτyz=σyzτzx=σzx\begin{aligned}\sigma_{x} &= \sigma_{xx}&\sigma_{y} &= \sigma_{yy}&\sigma_{z} &= \sigma_{zz}\\ \tau_{xy} &= \sigma_{xy}&\tau_{yz} &= \sigma_{yz}&\tau_{zx} &= \sigma_{zx}\end{aligned}

The principal stresses σ1,σ2,σ3\sigma_{1},\sigma_{2},\sigma_{3} are calculated as the roots (S)(S) of the cubic

S3I1S2+I2SI3=0S^{3} - I_{1}S^{2} + I_{2}S - I_{3} = 0

where the terms I1,I2,I3I_{1},I_{2},I_{3} are stress invariants defined as

I1=σx+σy+σzI2=σxσy+σyσz+σzσxτxy2τyz2τzx2I3=σxσyσz+2τxyτyzτzxσxτyz2σyτzx2σzτxy2\begin{aligned}I_{1} &= \sigma_{x} + \sigma_{y} + \sigma_{z}\\ I_{2} &= \sigma_{x}\sigma_{y} + \sigma_{y}\sigma_{z} + \sigma_{z}\sigma_{x} - {\tau_{xy}}^{2} - {\tau_{yz}}^{2} - {\tau_{zx}}^{2}\\ I_{3} &= \sigma_{x}\sigma_{y}\sigma_{z} + 2\tau_{xy}\tau_{yz}\tau_{zx} - \sigma_{x}{\tau_{yz}}^{2} - \sigma_{y}{\tau_{zx}}^{2} - \sigma_{z}{\tau_{xy}}^{2}\end{aligned}

Alternatively the principal stress equation can be written

S3(σx+σy+σz)S2+(σxσy+σyσz+σzσxτxy2τyz2τzx2)S(σxσyσz+2τxyτyzτzxσxτyz2σyτzx2σzτxy2)=0\begin{aligned}S^{3} - &\left( \sigma_{x} + \sigma_{y} + \sigma_{z} \right)S^{2} + \\ &\left( \sigma_{x}\sigma_{y} + \sigma_{y}\sigma_{z} + \sigma_{z}\sigma_{x} - {\tau_{xy}}^{2} - {\tau_{yz}}^{2} - {\tau_{zx}}^{2} \right)S - \\ &\left( \sigma_{x}\sigma_{y}\sigma_{z} + 2\tau_{xy}\tau_{yz}\tau_{zx} - \sigma_{x}{\tau_{yz}}^{2} - \sigma_{y}{\tau_{zx}}^{2} - \sigma_{z}{\tau_{xy}}^{2} \right) = 0\end{aligned}

The maximum shear stress is calculated from the principal stress as

τmax=12(σ1σ3){\tau_{\max}=\frac{1}{2}\left( \sigma_{1} - \sigma_{3} \right)}

Two other stress measures that are used are the average or hydrostatic stress and the von Mises stress; these are defined as

σ=13(σz+σy+σz)σvM=12[(σxσy)2+(σyσz)2+(σzσx)2+6(τxy2+τyz2+τzx2)]\begin{aligned}\sigma &= \frac{1}{3}( \sigma_z + \sigma_y + \sigma_z )\\ \sigma_{vM} &= \frac{1}{\sqrt{2}}\left\lbrack \left( \sigma_{x} - \sigma_{y} \right)^{2} + \left( \sigma_{y} - \sigma_{z} \right)^{2} + \left( \sigma_{z} - \sigma_{x} \right)^{2} + 6\left( {\tau_{xy}}^{2} + {\tau_{yz}}^{2} + {\tau_{zx}}^{2} \right) \right\rbrack\end{aligned}

Stress in 2D elements

The stress in 2D elements is calculated via the strains. The strains are calculated from the displacements using the strain displacement relationship B\mathbf{B} (see 2D elements). Using the interpolation functions these can be calculated at any point in the element. Once the strains are calculated the stress can be calculated using the material elastic matrix C\mathbf{C} for example for and elastic isotropic material the material matrix is

C=E1ν2[1νν11ν2]\mathbf{C} = \frac{E}{1 - \nu_{2}}\begin{bmatrix} 1 & \nu & \\ \nu & 1 & \\ & & \frac{1 - \nu}{2} \\ \end{bmatrix}

Thus the strains are

ε=Bu\boldsymbol{\varepsilon} = \mathbf{Bu}

and the stresses are

σ=Cε=CBu\boldsymbol{\sigma} = \mathbf{C}\boldsymbol{\varepsilon} = \mathbf{CBu}

This can be used to evaluate the stress at any point in the element. However the stress is based on the strain which in turn is based on the displacement gradients in the element. Thus some of the strain terms in an element that has a parabolic displacement field are linear. It has been found that the best stress results are obtained by evaluating the stress at particular points (the points used for the element integration) and extrapolating the results to the nodes.

In order to have good stress results the mesh will have to be finer that the mesh required for the displacement solution and the stress results are likely to be influenced by high displacement gradients in the element.

Direct extrapolation of results

In the case of direct extrapolation a function is chosen to represent the variation of stress over the element based on the number of Gauss points. In practice this is used when there are 1, 3 or 4 Gauss points. The corresponding polynomial functions are

f=af=a+bx+cyf=a+bx+cy+dxy\begin{aligned} f &= a\\ f &= a + bx + cy\\ f &= a + bx + cy + dxy\end{aligned}

For 1 Gauss point the values are assumed to be constant over the whole element so the Gauss point values are simply copied to the nodes.

For 3 Gauss points the values are the Gauss points are known so the following set of equations can be set up

f0=a+bx0+cy0f1=a+bx1+cy1f2=a+bx2+cy2\begin{aligned}f_{0} &= a + bx_{0} + cy_{0}\\ f_{1} &= a + bx_{1} + cy_{1}\\ f_{2} &= a + bx_{2} + cy_{2}\end{aligned}

This can then be used to calculate the coefficients a,b,ca,b,c

c=(f0f1)(x0x2)(f0f2)(x0x1)(y0y1)(x0x2)(y0y2)(x0x1)b=(f0f1)c(y0x1)(x0x1)a=f0bx0cy0\begin{aligned}c &= \frac{\left( f_{0} - f_{1} \right)\left( x_{0} - x_{2} \right) - \left( f_{0} - f_{2} \right)\left( x_{0} - x_{1} \right)}{\left( y_{0} - y_{1} \right)\left( x_{0} - x_{2} \right) - \left( y_{0} - y_{2} \right)\left( x_{0} - x_{1} \right)}\\ b &= \frac{\left( f_{0} - f_{1} \right) - c\left( y_{0} - x_{1} \right)}{\left( x_{0} - x_{1} \right)}\\ a &= f_{0} - bx_{0} - cy_{0}\end{aligned}

For 4 Gauss points a similar approach can be used, but in this case the locations of the Gauss points are at

r,s=±13r,s = \pm \frac{1}{\sqrt{3}}

so the equations can be written in the form

f0=abxcy+dxyf1=a+bxcy+dxyf2=a+bx+cy+dxyf3=abx+cydxy\begin{aligned}f_{0} &= a - bx - cy + dxy\\ f_{1} &= a + bx - cy + dxy\\ f_{2} &= a + bx + cy + dxy\\ f_{3} &= a - bx + cy - dxy\end{aligned}

This can then be solved for the coefficients a,b,c,da,b,c,d

a=f0+f1+f2+f34c=f0f1+f2+f34yb=f0+f1+f2f34xd=f0f1+f2f34xy\begin{aligned}a &= \frac{f_{0} + f_{1} + f_{2} + f_3}{4}&c &= \frac{-f_{0} - f_{1} + f_{2} + f_3}{4y}&\\ b &= \frac{- f_{0} + f_{1} + f_{2} - f_3}{4x}&d &= \frac{f_{0} - f_{1} + f_{2} - f_3}{4xy}\end{aligned}

Once these are established the polynomial functions can be used to establish the values at any position on the element.

Least squares extrapolation of results

In this case a function chosen to fit through the points would imply a higher order polynomial than the one used to interpolate the geometry, so a least squares approach is used to find the polynomial to map the stresses from the Gauss points to the nodes. The interpolation functions used for 6 node and 8 node elements are respectively

f=a0+a1x+a2y+a3xy+a4x2+a5y2f=a0+a1x+a2y+a3xy+a4x2+a5y2+a6x2y+a7xy2\begin{aligned}f &= a_{0} + a_{1}x + a_{2}y + a_{3}xy + a_{4}x^{2} + a_{5}y^{2}\\ f &= a_{0} + a_{1}x + a_{2}y + a_{3}xy + a_{4}x^{2} + a_{5}y^{2} + a_{6}x^{2}y + a_{7}xy^{2}\end{aligned}

The square of the error for any point is then

e2=(p(x,y)f)2e^{2} = \left( p(x,y) - f \right)_{2}

This is summed over all the Gauss points and then the derivatives with respect to the coefficients are set to zero (selecting the coefficients that minimise the error). This leads to the matrix equation for 8 node elements

[nxyxyx2y2x2yxy2x2xyx2yx3xy2x3yx2y2y2xy2x2yy3x2y2xy3x2y2x3yxy3x3y2x2y3x4x2y2x4yx3y2y4x2y3xy4x4y2x3y3x2y4]{a0a1a2a3a4a5a6a7}={ffxfyfxyfx2fy2fx2yfxy2}\begin{bmatrix} n & \sum_{}^{}x & \sum_{}^{}y & \sum_{}^{}{xy} & \sum_{}^{}x^{2} & \sum_{}^{}y^{2} & \sum_{}^{}{x^{2}y} & \sum_{}^{}{xy^{2}} \\ & \sum_{}^{}x^{2} & \sum_{}^{}{xy} & \sum_{}^{}{x^{2}y} & \sum_{}^{}x^{3} & \sum_{}^{}{xy^{2}} & \sum_{}^{}{x^{3}y} & \sum_{}^{}{x^{2}y^{2}} \\ & & \sum_{}^{}y^{2} & \sum_{}^{}{xy^{2}} & \sum_{}^{}{x^{2}y} & \sum_{}^{}y^{3} & \sum_{}^{}{x^{2}y^{2}} & \sum_{}^{}{xy^{3}} \\ & & & \sum_{}^{}{x^{2}y^{2}} & \sum_{}^{}{x^{3}y} & \sum_{}^{}{xy^{3}} & \sum_{}^{}{x^{3}y^{2}} & \sum_{}^{}{x^{2}y^{3}} \\ & & & & \sum_{}^{}x^{4} & \sum_{}^{}{x^{2}y^{2}} & \sum_{}^{}{x^{4}y} & \sum_{}^{}{x^{3}y^{2}} \\ & & & & & \sum_{}^{}y^{4} & \sum_{}^{}{x^{2}y^{3}} & \sum_{}^{}{xy^{4}} \\ & & & & & & \sum_{}^{}{x^{4}y^{2}} & \sum_{}^{}{x^{3}y^{3}} \\ & & & & & & & \sum_{}^{}{x^{2}y^{4}} \\ \end{bmatrix}\begin{Bmatrix} a_{0} \\ a_{1} \\ a_{2} \\ a_{3} \\ a_{4} \\ a_{5} \\ a_{6} \\ a_{7} \\ \end{Bmatrix} = \begin{Bmatrix} \sum_{}^{}f \\ \sum_{}^{}fx \\ \sum_{}^{}fy \\ \sum_{}^{}{fxy} \\ \sum_{}^{}{fx^{2}} \\ \sum_{}^{}{fy^{2}} \\ \sum_{}^{}{fx^{2}y} \\ \sum_{}^{}{fxy^{2}} \\ \end{Bmatrix}

The 6 node version is the same except that the a6a_{6} and a7a_{7} terms are ignored. This can then be solved for the coefficients.