Stress In 2D Elements Strain Definitions The normal definitions of strain used are as follows
ε x x = ∂ u ∂ x ε x x = ∂ v ∂ y ε x x = ∂ w ∂ z γ x y = ∂ u ∂ y + ∂ v ∂ x γ y z = ∂ v ∂ z + ∂ w ∂ y γ z x = ∂ w ∂ x + ∂ u ∂ z \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} ε x x γ x y = ∂ x ∂ u = ∂ y ∂ u + ∂ x ∂ v ε x x γ y z = ∂ y ∂ v = ∂ z ∂ v + ∂ y ∂ w ε x x γ z x = ∂ z ∂ w = ∂ x ∂ w + ∂ z ∂ u An alternative definition which fits more neatly in tensor form is
ε x x = ∂ u ∂ x ε x x = ∂ v ∂ y ε x x = ∂ w ∂ z ε x y = 1 2 ( ∂ u ∂ y + ∂ v ∂ x ) ε y z = 1 2 ( ∂ v ∂ z + ∂ w ∂ y ) ε z x = 1 2 ( ∂ w ∂ x + ∂ u ∂ z ) \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} ε x x ε x y = ∂ x ∂ u = 2 1 ( ∂ y ∂ u + ∂ x ∂ v ) ε x x ε y z = ∂ y ∂ v = 2 1 ( ∂ z ∂ v + ∂ y ∂ w ) ε x x ε z x = ∂ z ∂ w = 2 1 ( ∂ x ∂ w + ∂ z ∂ u ) with the strain tensor defined as
ε = [ ε x x ε x y ε x z ε y x ε y y ε y z ε z x ε z x ε z z ] \mathbf{\varepsilon} = \begin{bmatrix} \varepsilon_{xx} & \varepsilon_{xy} & \varepsilon_{xz} \\ \varepsilon_{yx} & \varepsilon_{yy} & \varepsilon_{yz} \\ \varepsilon_{zx} & \varepsilon_{zx} & \varepsilon_{zz} \\ \end{bmatrix} ε = ⎣ ⎢ ⎡ ε x x ε y x ε z x ε x y ε y y ε z x ε x z ε y z ε z z ⎦ ⎥ ⎤ The calculation of principal strains
ε 1 , ε 2 , ε 3 \varepsilon_{1},\varepsilon_{2},\varepsilon_{3} ε 1 , ε 2 , ε 3 follows from
E 3 − ( ε x x + ε y y + ε z z ) E 2 + ( ε x x ε y y + ε y y ε z z + ε z z ε x x − ε x y 2 − ε y z 2 − ε z x 2 ) E − ( ε x x ε y y ε z z + 2 ε x y ε y z ε z x − ε x x ε y z 2 − ε y y ε z x 2 − ε z z ε x y 2 ) = 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} E 3 − ( ε x x + ε y y + ε z z ) E 2 + ( ε x x ε y y + ε y y ε z z + ε z z ε x x − ε x y 2 − ε y z 2 − ε z x 2 ) E − ( ε x x ε y y ε z z + 2 ε x y ε y z ε z x − ε x x ε y z 2 − ε y y ε z x 2 − ε z z ε x y 2 ) = 0 The maximum shear strain is calculated from the principal strain as
ε max shear = 1 2 ( ε 1 − ε 3 ) \varepsilon_{\text{max shear}} = \frac{1}{2}\left( \varepsilon_{1} - \varepsilon_{3} \right) ε max shear = 2 1 ( ε 1 − ε 3 ) or
γ max shear = ( ε 1 − ε 3 ) \gamma_{\text{max shear}} = \left( \varepsilon_{1} - \varepsilon_{3} \right) γ max shear = ( ε 1 − ε 3 ) In a similar way to the definitions of average and von Mises stress a
volumetric and effective strain can be calculated as
ε a v = 1 3 ( ε x x + ε y y + ε z z ) ε v M = 1 2 ( 1 + ν ) [ ( ε x x − ε y y ) 2 + ( ε y y − ε z z ) 2 + ( ε z z − ε x x ) 2 + 6 ( ε x y 2 + ε y z 2 + ε z x 2 ) ] \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} ε a v ε v M = 3 1 ( ε x x + ε y y + ε z z ) = 2 ( 1 + ν ) 1 [ ( ε x x − ε y y ) 2 + ( ε y y − ε z z ) 2 + ( ε z z − ε x x ) 2 + 6 ( ε x y 2 + ε y z 2 + ε z x 2 ) ] Stress Definitions Stress can be considered as a tensor quantity whose components can be
represented in matrix form as
σ = [ σ x x σ x y σ x z σ y x σ y y σ y z σ z x σ z x σ z z ] \boldsymbol{\sigma} = \begin{bmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{zx} & \sigma_{zx} & \sigma_{zz} \\ \end{bmatrix} σ = ⎣ ⎢ ⎡ σ x x σ y x σ z x σ x y σ y y σ z x σ x z σ y z σ z z ⎦ ⎥ ⎤ where each term corresponds to a force per unit area. The following
notation for the stress components is common
σ x = σ x x σ y = σ y y σ z = σ z z τ x y = σ x y τ y z = σ y z τ z x = σ z x \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} σ x τ x y = σ x x = σ x y σ y τ y z = σ y y = σ y z σ z τ z x = σ z z = σ z x The principal stresses σ 1 , σ 2 , σ 3 \sigma_{1},\sigma_{2},\sigma_{3} σ 1 , σ 2 , σ 3 are calculated
as the roots ( S ) (S) ( S ) of the cubic
S 3 − I 1 S 2 + I 2 S − I 3 = 0 S^{3} - I_{1}S^{2} + I_{2}S - I_{3} = 0 S 3 − I 1 S 2 + I 2 S − I 3 = 0 where the terms I 1 , I 2 , I 3 I_{1},I_{2},I_{3} I 1 , I 2 , I 3 are stress invariants defined as
I 1 = σ x + σ y + σ z I 2 = σ x σ y + σ y σ z + σ z σ x − τ x y 2 − τ y z 2 − τ z x 2 I 3 = σ x σ y σ z + 2 τ x y τ y z τ z x − σ x τ y z 2 − σ y τ z x 2 − σ z τ x y 2 \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} I 1 I 2 I 3 = σ x + σ y + σ z = σ x σ y + σ y σ z + σ z σ x − τ x y 2 − τ y z 2 − τ z x 2 = σ x σ y σ z + 2 τ x y τ y z τ z x − σ x τ y z 2 − σ y τ z x 2 − σ z τ x y 2 Alternatively the principal stress equation can be written
S 3 − ( σ x + σ y + σ z ) S 2 + ( σ x σ y + σ y σ z + σ z σ x − τ x y 2 − τ y z 2 − τ z x 2 ) S − ( σ x σ y σ z + 2 τ x y τ y z τ z x − σ x τ y z 2 − σ y τ z x 2 − σ z τ x y 2 ) = 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} S 3 − ( σ x + σ y + σ z ) S 2 + ( σ x σ y + σ y σ z + σ z σ x − τ x y 2 − τ y z 2 − τ z x 2 ) S − ( σ x σ y σ z + 2 τ x y τ y z τ z x − σ x τ y z 2 − σ y τ z x 2 − σ z τ x y 2 ) = 0 The maximum shear stress is calculated from the principal stress as
τ max = 1 2 ( σ 1 − σ 3 ) {\tau_{\max}=\frac{1}{2}\left( \sigma_{1} - \sigma_{3} \right)} τ m a x = 2 1 ( σ 1 − σ 3 ) Two other stress measures that are used are the average or hydrostatic
stress and the von Mises stress; these are defined as
σ = 1 3 ( σ z + σ y + σ z ) σ v M = 1 2 [ ( σ x − σ y ) 2 + ( σ y − σ z ) 2 + ( σ z − σ x ) 2 + 6 ( τ x y 2 + τ y z 2 + τ z x 2 ) ] \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} σ σ v M = 3 1 ( σ z + σ y + σ z ) = 2 1 [ ( σ x − σ y ) 2 + ( σ y − σ z ) 2 + ( σ z − σ x ) 2 + 6 ( τ x y 2 + τ y z 2 + τ z x 2 ) ] 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} 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} C for example for and elastic isotropic
material the material matrix is
C = E 1 − ν 2 [ 1 ν ν 1 1 − ν 2 ] \mathbf{C} = \frac{E}{1 - \nu_{2}}\begin{bmatrix} 1 & \nu & \\ \nu & 1 & \\ & & \frac{1 - \nu}{2} \\ \end{bmatrix} C = 1 − ν 2 E ⎣ ⎢ ⎡ 1 ν ν 1 2 1 − ν ⎦ ⎥ ⎤ Thus the strains are
ε = B u \boldsymbol{\varepsilon} = \mathbf{Bu} ε = B u and the stresses are
σ = C ε = C B u \boldsymbol{\sigma} = \mathbf{C}\boldsymbol{\varepsilon} = \mathbf{CBu} σ = C ε = C B u 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.
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 = a f = a + b x + c y f = a + b x + c y + d x y \begin{aligned} f &= a\\ f &= a + bx + cy\\ f &= a + bx + cy + dxy\end{aligned} f f f = a = a + b x + c y = a + b x + c y + d x y 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
f 0 = a + b x 0 + c y 0 f 1 = a + b x 1 + c y 1 f 2 = a + b x 2 + c y 2 \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} f 0 f 1 f 2 = a + b x 0 + c y 0 = a + b x 1 + c y 1 = a + b x 2 + c y 2 This can then be used to calculate the coefficients a , b , c a,b,c a , b , c
c = ( f 0 − f 1 ) ( x 0 − x 2 ) − ( f 0 − f 2 ) ( x 0 − x 1 ) ( y 0 − y 1 ) ( x 0 − x 2 ) − ( y 0 − y 2 ) ( x 0 − x 1 ) b = ( f 0 − f 1 ) − c ( y 0 − x 1 ) ( x 0 − x 1 ) a = f 0 − b x 0 − c y 0 \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} c b a = ( y 0 − y 1 ) ( x 0 − x 2 ) − ( y 0 − y 2 ) ( x 0 − x 1 ) ( f 0 − f 1 ) ( x 0 − x 2 ) − ( f 0 − f 2 ) ( x 0 − x 1 ) = ( x 0 − x 1 ) ( f 0 − f 1 ) − c ( y 0 − x 1 ) = f 0 − b x 0 − c y 0 For 4 Gauss points a similar approach can be used, but in this case the
locations of the Gauss points are at
r , s = ± 1 3 r,s = \pm \frac{1}{\sqrt{3}} r , s = ± 3 1 so the equations can be written in the form
f 0 = a − b x − c y + d x y f 1 = a + b x − c y + d x y f 2 = a + b x + c y + d x y f 3 = a − b x + c y − d x y \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} f 0 f 1 f 2 f 3 = a − b x − c y + d x y = a + b x − c y + d x y = a + b x + c y + d x y = a − b x + c y − d x y This can then be solved for the coefficients a , b , c , d a,b,c,d a , b , c , d
a = f 0 + f 1 + f 2 + f 3 4 c = − f 0 − f 1 + f 2 + f 3 4 y b = − f 0 + f 1 + f 2 − f 3 4 x d = f 0 − f 1 + f 2 − f 3 4 x y \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} a b = 4 f 0 + f 1 + f 2 + f 3 = 4 x − f 0 + f 1 + f 2 − f 3 c d = 4 y − f 0 − f 1 + f 2 + f 3 = 4 x y f 0 − f 1 + f 2 − f 3 Once these are established the polynomial functions can be used to
establish the values at any position on the element.
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 = a 0 + a 1 x + a 2 y + a 3 x y + a 4 x 2 + a 5 y 2 f = a 0 + a 1 x + a 2 y + a 3 x y + a 4 x 2 + a 5 y 2 + a 6 x 2 y + a 7 x y 2 \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} f f = a 0 + a 1 x + a 2 y + a 3 x y + a 4 x 2 + a 5 y 2 = a 0 + a 1 x + a 2 y + a 3 x y + a 4 x 2 + a 5 y 2 + a 6 x 2 y + a 7 x y 2 The square of the error for any point is then
e 2 = ( p ( x , y ) − f ) 2 e^{2} = \left( p(x,y) - f \right)_{2} e 2 = ( p ( x , y ) − f ) 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
[ n ∑ x ∑ y ∑ x y ∑ x 2 ∑ y 2 ∑ x 2 y ∑ x y 2 ∑ x 2 ∑ x y ∑ x 2 y ∑ x 3 ∑ x y 2 ∑ x 3 y ∑ x 2 y 2 ∑ y 2 ∑ x y 2 ∑ x 2 y ∑ y 3 ∑ x 2 y 2 ∑ x y 3 ∑ x 2 y 2 ∑ x 3 y ∑ x y 3 ∑ x 3 y 2 ∑ x 2 y 3 ∑ x 4 ∑ x 2 y 2 ∑ x 4 y ∑ x 3 y 2 ∑ y 4 ∑ x 2 y 3 ∑ x y 4 ∑ x 4 y 2 ∑ x 3 y 3 ∑ x 2 y 4 ] { a 0 a 1 a 2 a 3 a 4 a 5 a 6 a 7 } = { ∑ f ∑ f x ∑ f y ∑ f x y ∑ f x 2 ∑ f y 2 ∑ f x 2 y ∑ f x y 2 } \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} ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ n ∑ x ∑ x 2 ∑ y ∑ x y ∑ y 2 ∑ x y ∑ x 2 y ∑ x y 2 ∑ x 2 y 2 ∑ x 2 ∑ x 3 ∑ x 2 y ∑ x 3 y ∑ x 4 ∑ y 2 ∑ x y 2 ∑ y 3 ∑ x y 3 ∑ x 2 y 2 ∑ y 4 ∑ x 2 y ∑ x 3 y ∑ x 2 y 2 ∑ x 3 y 2 ∑ x 4 y ∑ x 2 y 3 ∑ x 4 y 2 ∑ x y 2 ∑ x 2 y 2 ∑ x y 3 ∑ x 2 y 3 ∑ x 3 y 2 ∑ x y 4 ∑ x 3 y 3 ∑ x 2 y 4 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎧ a 0 a 1 a 2 a 3 a 4 a 5 a 6 a 7 ⎭ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎫ = ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎧ ∑ f ∑ f x ∑ f y ∑ f x y ∑ f x 2 ∑ f y 2 ∑ f x 2 y ∑ f x y 2 ⎭ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎫
The 6 node version is the same except that the a 6 a_{6} a 6 and a 7 a_{7} a 7 terms
are ignored. This can then be solved for the coefficients.