FEA TipsFinite Element Theory & MATLAB Programs
Contents/Section 25.1

Section 25.1 of Chapter 25: Hexahedral isoparametric finite element. Linear case

Hexahedral 8-node isoparametric finite element for linear analysis

01

Element and hypotheses

StrainSmall
DisplacementSmall
MaterialLinear elastic

The hexahedral 8-node isoparametric finite element is the equivalent of the 2D quadrilateral element presented in Section 9.4. For more details, refer to [1]. This finite element has eight nodes and three DOFs per node:

Fel=keluel\htmlClass{fea-math-bold-italic}{\boldsymbol{F}}_{el}=\htmlClass{fea-math-bold-italic}{\boldsymbol{k}}_{el}\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}_{el}
Fel={X1Y1Z1X2Y2Z2X3Y3Z3X4Y4Z4⋯X8Y8Z8}T\htmlClass{fea-math-bold-italic}{\boldsymbol{F}}_{el}=\begin{Bmatrix}X_1&Y_1&Z_1&X_2&Y_2&Z_2&X_3&Y_3&Z_3&X_4&Y_4&Z_4&\cdots&X_8&Y_8&Z_8\end{Bmatrix}^{T}
uel={u1v1w1u2v2w2u3v3w3u4v4w4⋯u8v8w8}T\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}_{el}=\begin{Bmatrix}u_1&v_1&w_1&u_2&v_2&w_2&u_3&v_3&w_3&u_4&v_4&w_4&\cdots&u_8&v_8&w_8\end{Bmatrix}^{T}

Fel\htmlClass{fea-math-bold-italic}{\boldsymbol{F}}_{el} denotes the nodal forces, uel\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}_{el} the nodal displacements, and kel\htmlClass{fea-math-bold-italic}{\boldsymbol{k}}_{el} the stiffness matrix of the finite element.

The theory is briefly presented, following Section 9.4.

02

Natural coordinates and interpolation

HEX8 natural and Cartesian coordinates
Figure 1. Natural coordinates s,t,rs,t,r and Cartesian coordinates x,y,zx,y,z; HEX8 node numbering.

Interpolation functions:

{u(s,t,r)v(s,t,r)w(s,t,r)}=∑i=18[Ni(s,t,r)000Ni(s,t,r)000Ni(s,t,r)]{uiviwi}=[NuNvNw]uel=Nuel\begin{Bmatrix}u(s,t,r)\\v(s,t,r)\\w(s,t,r)\end{Bmatrix}=\sum_{i=1}^{8}\begin{bmatrix}N_i(s,t,r)&0&0\\0&N_i(s,t,r)&0\\0&0&N_i(s,t,r)\end{bmatrix}\begin{Bmatrix}u_i\\v_i\\w_i\end{Bmatrix}=\begin{bmatrix}\htmlClass{fea-math-bold-italic}{\boldsymbol{N}}_u\\\htmlClass{fea-math-bold-italic}{\boldsymbol{N}}_v\\\htmlClass{fea-math-bold-italic}{\boldsymbol{N}}_w\end{bmatrix}\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}_{el}=\htmlClass{fea-math-bold-italic}{\boldsymbol{N}}\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}_{el}

The same interpolation functions are also used for the Cartesian coordinates:

{x(s,t,r)y(s,t,r)z(s,t,r)}=∑i=18[Ni(s,t,r)000Ni(s,t,r)000Ni(s,t,r)]{xiyizi}=[NuNvNw]xel=Nxel\begin{Bmatrix}x(s,t,r)\\y(s,t,r)\\z(s,t,r)\end{Bmatrix}=\sum_{i=1}^{8}\begin{bmatrix}N_i(s,t,r)&0&0\\0&N_i(s,t,r)&0\\0&0&N_i(s,t,r)\end{bmatrix}\begin{Bmatrix}x_i\\y_i\\z_i\end{Bmatrix}=\begin{bmatrix}\htmlClass{fea-math-bold-italic}{\boldsymbol{N}}_u\\\htmlClass{fea-math-bold-italic}{\boldsymbol{N}}_v\\\htmlClass{fea-math-bold-italic}{\boldsymbol{N}}_w\end{bmatrix}\htmlClass{fea-math-bold-italic}{\boldsymbol{x}}_{el}=\htmlClass{fea-math-bold-italic}{\boldsymbol{N}}\htmlClass{fea-math-bold-italic}{\boldsymbol{x}}_{el}

xel\htmlClass{fea-math-bold-italic}{\boldsymbol{x}}_{el} contains the nodal coordinates:

xel={x1y1z1x2y2z2x3y3z3x4y4z4⋯x8y8z8}T\htmlClass{fea-math-bold-italic}{\boldsymbol{x}}_{el}=\begin{Bmatrix}x_1&y_1&z_1&x_2&y_2&z_2&x_3&y_3&z_3&x_4&y_4&z_4&\cdots&x_8&y_8&z_8\end{Bmatrix}^{T}

The interpolation functions are:

N1(s,t,r)=(1+s)(1+t)(1+r)8N2(s,t,r)=(1−s)(1+t)(1+r)8N3(s,t,r)=(1−s)(1−t)(1+r)8N4(s,t,r)=(1+s)(1−t)(1+r)8N5(s,t,r)=(1+s)(1+t)(1−r)8N6(s,t,r)=(1−s)(1+t)(1−r)8N7(s,t,r)=(1−s)(1−t)(1−r)8N8(s,t,r)=(1+s)(1−t)(1−r)8\begin{aligned}N_{1}(s,t,r)&=\frac{(1+s)(1+t)(1+r)}{8}&\qquad N_{2}(s,t,r)&=\frac{(1-s)(1+t)(1+r)}{8}\\[5pt]N_{3}(s,t,r)&=\frac{(1-s)(1-t)(1+r)}{8}&\qquad N_{4}(s,t,r)&=\frac{(1+s)(1-t)(1+r)}{8}\\[5pt]N_{5}(s,t,r)&=\frac{(1+s)(1+t)(1-r)}{8}&\qquad N_{6}(s,t,r)&=\frac{(1-s)(1+t)(1-r)}{8}\\[5pt]N_{7}(s,t,r)&=\frac{(1-s)(1-t)(1-r)}{8}&\qquad N_{8}(s,t,r)&=\frac{(1+s)(1-t)(1-r)}{8}\end{aligned}
s,t,r∈[−1,1]s,t,r\in[-1,1]

03

Strains and Jacobian

Strains:

ε={εxεyεzγxyγyzγzx}={∂u∂x∂v∂y∂w∂z∂u∂y+∂v∂x∂v∂z+∂w∂y∂w∂x+∂u∂z}\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}=\begin{Bmatrix}\varepsilon_x\\\varepsilon_y\\\varepsilon_z\\\gamma_{xy}\\\gamma_{yz}\\\gamma_{zx}\end{Bmatrix}=\begin{Bmatrix}\dfrac{\partial u}{\partial x}\\\dfrac{\partial v}{\partial y}\\\dfrac{\partial w}{\partial z}\\\dfrac{\partial u}{\partial y}+\dfrac{\partial v}{\partial x}\\\dfrac{\partial v}{\partial z}+\dfrac{\partial w}{\partial y}\\\dfrac{\partial w}{\partial x}+\dfrac{\partial u}{\partial z}\end{Bmatrix}
ε={∂u∂x00∂u∂y0∂u∂z}+{0∂v∂y0∂v∂x∂v∂z0}+{00∂w∂z0∂w∂y∂w∂x}\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}=\begin{Bmatrix}\dfrac{\partial u}{\partial x}\\0\\0\\\dfrac{\partial u}{\partial y}\\0\\\dfrac{\partial u}{\partial z}\end{Bmatrix}+\begin{Bmatrix}0\\\dfrac{\partial v}{\partial y}\\0\\\dfrac{\partial v}{\partial x}\\\dfrac{\partial v}{\partial z}\\0\end{Bmatrix}+\begin{Bmatrix}0\\0\\\dfrac{\partial w}{\partial z}\\0\\\dfrac{\partial w}{\partial y}\\\dfrac{\partial w}{\partial x}\end{Bmatrix}
ε=[100000000010000001]{∂u∂x∂u∂y∂u∂z}+[000010000100001000]{∂v∂x∂v∂y∂v∂z}+[000000001000010100]{∂w∂x∂w∂y∂w∂z}\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}=\begin{bmatrix}1&0&0\\0&0&0\\0&0&0\\0&1&0\\0&0&0\\0&0&1\end{bmatrix}\begin{Bmatrix}\dfrac{\partial u}{\partial x}\\\dfrac{\partial u}{\partial y}\\\dfrac{\partial u}{\partial z}\end{Bmatrix}+\begin{bmatrix}0&0&0\\0&1&0\\0&0&0\\1&0&0\\0&0&1\\0&0&0\end{bmatrix}\begin{Bmatrix}\dfrac{\partial v}{\partial x}\\\dfrac{\partial v}{\partial y}\\\dfrac{\partial v}{\partial z}\end{Bmatrix}+\begin{bmatrix}0&0&0\\0&0&0\\0&0&1\\0&0&0\\0&1&0\\1&0&0\end{bmatrix}\begin{Bmatrix}\dfrac{\partial w}{\partial x}\\\dfrac{\partial w}{\partial y}\\\dfrac{\partial w}{\partial z}\end{Bmatrix}

With:

{∂u∂x∂u∂y∂u∂z}=J−1{∂u∂s∂u∂t∂u∂r}{∂v∂x∂v∂y∂v∂z}=J−1{∂v∂s∂v∂t∂v∂r}{∂w∂x∂w∂y∂w∂z}=J−1{∂w∂s∂w∂t∂w∂r}\begin{Bmatrix}\dfrac{\partial u}{\partial x}\\\dfrac{\partial u}{\partial y}\\\dfrac{\partial u}{\partial z}\end{Bmatrix}=\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}^{-1}\begin{Bmatrix}\dfrac{\partial u}{\partial s}\\\dfrac{\partial u}{\partial t}\\\dfrac{\partial u}{\partial r}\end{Bmatrix}\qquad \begin{Bmatrix}\dfrac{\partial v}{\partial x}\\\dfrac{\partial v}{\partial y}\\\dfrac{\partial v}{\partial z}\end{Bmatrix}=\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}^{-1}\begin{Bmatrix}\dfrac{\partial v}{\partial s}\\\dfrac{\partial v}{\partial t}\\\dfrac{\partial v}{\partial r}\end{Bmatrix}\qquad \begin{Bmatrix}\dfrac{\partial w}{\partial x}\\\dfrac{\partial w}{\partial y}\\\dfrac{\partial w}{\partial z}\end{Bmatrix}=\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}^{-1}\begin{Bmatrix}\dfrac{\partial w}{\partial s}\\\dfrac{\partial w}{\partial t}\\\dfrac{\partial w}{\partial r}\end{Bmatrix}

J\htmlClass{fea-math-bold-italic}{\boldsymbol{J}} is the Jacobian:

J(s,t,r)=[∂x∂s∂y∂s∂z∂s∂x∂t∂y∂t∂z∂t∂x∂r∂y∂r∂z∂r]\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}(s,t,r)=\begin{bmatrix}\dfrac{\partial x}{\partial s}&\dfrac{\partial y}{\partial s}&\dfrac{\partial z}{\partial s}\\\dfrac{\partial x}{\partial t}&\dfrac{\partial y}{\partial t}&\dfrac{\partial z}{\partial t}\\\dfrac{\partial x}{\partial r}&\dfrac{\partial y}{\partial r}&\dfrac{\partial z}{\partial r}\end{bmatrix}

Using the interpolation functions:

J(s,t,r)=[∂x∂s∂y∂s∂z∂s∂x∂t∂y∂t∂z∂t∂x∂r∂y∂r∂z∂r]=H(s,t,r)[x1y1z1x2y2z2x3y3z3⋮⋮⋮x8y8z8]\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}(s,t,r)=\begin{bmatrix}\dfrac{\partial x}{\partial s}&\dfrac{\partial y}{\partial s}&\dfrac{\partial z}{\partial s}\\\dfrac{\partial x}{\partial t}&\dfrac{\partial y}{\partial t}&\dfrac{\partial z}{\partial t}\\\dfrac{\partial x}{\partial r}&\dfrac{\partial y}{\partial r}&\dfrac{\partial z}{\partial r}\end{bmatrix}=\htmlClass{fea-math-bold-italic}{\boldsymbol{H}}(s,t,r)\begin{bmatrix}x_1&y_1&z_1\\x_2&y_2&z_2\\x_3&y_3&z_3\\\vdots&\vdots&\vdots\\x_8&y_8&z_8\end{bmatrix}

With:

H(s,t,r)=18[(1+t)(1+r)(1+s)(1+r)(1+s)(1+t)−(1+t)(1+r)(1−s)(1+r)(1−s)(1+t)−(1−t)(1+r)−(1−s)(1+r)(1−s)(1−t)(1−t)(1+r)−(1+s)(1+r)(1+s)(1−t)(1+t)(1−r)(1+s)(1−r)−(1+s)(1+t)−(1+t)(1−r)(1−s)(1−r)−(1−s)(1+t)−(1−t)(1−r)−(1−s)(1−r)−(1−s)(1−t)(1−t)(1−r)−(1+s)(1−r)−(1+s)(1−t)]T\htmlClass{fea-math-bold-italic}{\boldsymbol{H}}(s,t,r)=\frac18\begin{bmatrix}(1+t)(1+r)&(1+s)(1+r)&(1+s)(1+t)\\-(1+t)(1+r)&(1-s)(1+r)&(1-s)(1+t)\\-(1-t)(1+r)&-(1-s)(1+r)&(1-s)(1-t)\\(1-t)(1+r)&-(1+s)(1+r)&(1+s)(1-t)\\(1+t)(1-r)&(1+s)(1-r)&-(1+s)(1+t)\\-(1+t)(1-r)&(1-s)(1-r)&-(1-s)(1+t)\\-(1-t)(1-r)&-(1-s)(1-r)&-(1-s)(1-t)\\(1-t)(1-r)&-(1+s)(1-r)&-(1+s)(1-t)\end{bmatrix}^{T}

04

Strain-displacement matrix

Therefore:

ε={εxεyεzγxyγyzγzx}=[100000000010000001]J−1{∂u∂s∂u∂t∂u∂r}+[000010000100001000]J−1{∂v∂s∂v∂t∂v∂r}+[000000001000010100]J−1{∂w∂s∂w∂t∂w∂r}\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}=\begin{Bmatrix}\varepsilon_x\\\varepsilon_y\\\varepsilon_z\\\gamma_{xy}\\\gamma_{yz}\\\gamma_{zx}\end{Bmatrix}=\begin{bmatrix}1&0&0\\0&0&0\\0&0&0\\0&1&0\\0&0&0\\0&0&1\end{bmatrix}\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}^{-1}\begin{Bmatrix}\dfrac{\partial u}{\partial s}\\\dfrac{\partial u}{\partial t}\\\dfrac{\partial u}{\partial r}\end{Bmatrix}+\begin{bmatrix}0&0&0\\0&1&0\\0&0&0\\1&0&0\\0&0&1\\0&0&0\end{bmatrix}\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}^{-1}\begin{Bmatrix}\dfrac{\partial v}{\partial s}\\\dfrac{\partial v}{\partial t}\\\dfrac{\partial v}{\partial r}\end{Bmatrix}+\begin{bmatrix}0&0&0\\0&0&0\\0&0&1\\0&0&0\\0&1&0\\1&0&0\end{bmatrix}\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}^{-1}\begin{Bmatrix}\dfrac{\partial w}{\partial s}\\\dfrac{\partial w}{\partial t}\\\dfrac{\partial w}{\partial r}\end{Bmatrix}

Using the interpolation functions again:

{∂u∂s∂u∂t∂u∂r}=H(s,t,r){u1u2⋮u8}{∂v∂s∂v∂t∂v∂r}=H(s,t,r){v1v2⋮v8}{∂w∂s∂w∂t∂w∂r}=H(s,t,r){w1w2⋮w8}\begin{Bmatrix}\dfrac{\partial u}{\partial s}\\\dfrac{\partial u}{\partial t}\\\dfrac{\partial u}{\partial r}\end{Bmatrix}=\htmlClass{fea-math-bold-italic}{\boldsymbol{H}}(s,t,r)\begin{Bmatrix}u_1\\u_2\\\vdots\\u_8\end{Bmatrix}\qquad \begin{Bmatrix}\dfrac{\partial v}{\partial s}\\\dfrac{\partial v}{\partial t}\\\dfrac{\partial v}{\partial r}\end{Bmatrix}=\htmlClass{fea-math-bold-italic}{\boldsymbol{H}}(s,t,r)\begin{Bmatrix}v_1\\v_2\\\vdots\\v_8\end{Bmatrix}\qquad \begin{Bmatrix}\dfrac{\partial w}{\partial s}\\\dfrac{\partial w}{\partial t}\\\dfrac{\partial w}{\partial r}\end{Bmatrix}=\htmlClass{fea-math-bold-italic}{\boldsymbol{H}}(s,t,r)\begin{Bmatrix}w_1\\w_2\\\vdots\\w_8\end{Bmatrix}

Therefore:

ε={εxεyεzγxyγyzγzx}=[100000000010000001]⏟Pub{u1u2⋮u8}+[000010000100001000]⏟Pvb{v1v2⋮v8}+[000000001000010100]⏟Pwb{w1w2⋮w8}\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}=\begin{Bmatrix}\varepsilon_x\\\varepsilon_y\\\varepsilon_z\\\gamma_{xy}\\\gamma_{yz}\\\gamma_{zx}\end{Bmatrix}=\underbrace{\begin{bmatrix}1&0&0\\0&0&0\\0&0&0\\0&1&0\\0&0&0\\0&0&1\end{bmatrix}}_{\htmlClass{fea-math-bold-italic}{\boldsymbol{P}}_u}\htmlClass{fea-math-bold-italic}{\boldsymbol{b}}\begin{Bmatrix}u_1\\u_2\\\vdots\\u_8\end{Bmatrix}+\underbrace{\begin{bmatrix}0&0&0\\0&1&0\\0&0&0\\1&0&0\\0&0&1\\0&0&0\end{bmatrix}}_{\htmlClass{fea-math-bold-italic}{\boldsymbol{P}}_v}\htmlClass{fea-math-bold-italic}{\boldsymbol{b}}\begin{Bmatrix}v_1\\v_2\\\vdots\\v_8\end{Bmatrix}+\underbrace{\begin{bmatrix}0&0&0\\0&0&0\\0&0&1\\0&0&0\\0&1&0\\1&0&0\end{bmatrix}}_{\htmlClass{fea-math-bold-italic}{\boldsymbol{P}}_w}\htmlClass{fea-math-bold-italic}{\boldsymbol{b}}\begin{Bmatrix}w_1\\w_2\\\vdots\\w_8\end{Bmatrix}

Where:

b=J−1H=[b11b12b13b14b15b16b17b18b21b22b23b24b25b26b27b28b31b32b33b34b35b36b37b38]\htmlClass{fea-math-bold-italic}{\boldsymbol{b}}=\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}^{-1}\htmlClass{fea-math-bold-italic}{\boldsymbol{H}}=\begin{bmatrix}b_{11}&b_{12}&b_{13}&b_{14}&b_{15}&b_{16}&b_{17}&b_{18}\\b_{21}&b_{22}&b_{23}&b_{24}&b_{25}&b_{26}&b_{27}&b_{28}\\b_{31}&b_{32}&b_{33}&b_{34}&b_{35}&b_{36}&b_{37}&b_{38}\end{bmatrix}

It follows that B\htmlClass{fea-math-bold-italic}{\boldsymbol{B}} is a 6×246\times24 matrix:

ε=Pub{u1u2⋮u8}+Pvb{v1v2⋮v8}+Pwb{w1w2⋮w8}=Buel\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}=\htmlClass{fea-math-bold-italic}{\boldsymbol{P}}_u\htmlClass{fea-math-bold-italic}{\boldsymbol{b}}\begin{Bmatrix}u_1\\u_2\\\vdots\\u_8\end{Bmatrix}+\htmlClass{fea-math-bold-italic}{\boldsymbol{P}}_v\htmlClass{fea-math-bold-italic}{\boldsymbol{b}}\begin{Bmatrix}v_1\\v_2\\\vdots\\v_8\end{Bmatrix}+\htmlClass{fea-math-bold-italic}{\boldsymbol{P}}_w\htmlClass{fea-math-bold-italic}{\boldsymbol{b}}\begin{Bmatrix}w_1\\w_2\\\vdots\\w_8\end{Bmatrix}=\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}_{el}

B\htmlClass{fea-math-bold-italic}{\boldsymbol{B}} matrix:

B=[b1100b1200b1300b1400b1500b1600b1700b18000b2100b2200b2300b2400b2500b2600b2700b28000b3100b3200b3300b3400b3500b3600b3700b38b21b110b22b120b23b130b24b140b25b150b26b160b27b170b28b1800b31b210b32b220b33b230b34b240b35b250b36b260b37b270b38b28b310b11b320b12b330b13b340b14b350b15b360b16b370b17b380b18]\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}=\begin{bmatrix}b_{11}&0&0&b_{12}&0&0&b_{13}&0&0&b_{14}&0&0&b_{15}&0&0&b_{16}&0&0&b_{17}&0&0&b_{18}&0&0\\0&b_{21}&0&0&b_{22}&0&0&b_{23}&0&0&b_{24}&0&0&b_{25}&0&0&b_{26}&0&0&b_{27}&0&0&b_{28}&0\\0&0&b_{31}&0&0&b_{32}&0&0&b_{33}&0&0&b_{34}&0&0&b_{35}&0&0&b_{36}&0&0&b_{37}&0&0&b_{38}\\b_{21}&b_{11}&0&b_{22}&b_{12}&0&b_{23}&b_{13}&0&b_{24}&b_{14}&0&b_{25}&b_{15}&0&b_{26}&b_{16}&0&b_{27}&b_{17}&0&b_{28}&b_{18}&0\\0&b_{31}&b_{21}&0&b_{32}&b_{22}&0&b_{33}&b_{23}&0&b_{34}&b_{24}&0&b_{35}&b_{25}&0&b_{36}&b_{26}&0&b_{37}&b_{27}&0&b_{38}&b_{28}\\b_{31}&0&b_{11}&b_{32}&0&b_{12}&b_{33}&0&b_{13}&b_{34}&0&b_{14}&b_{35}&0&b_{15}&b_{36}&0&b_{16}&b_{37}&0&b_{17}&b_{38}&0&b_{18}\end{bmatrix}

05

Constitutive law and stiffness matrix

The stiffness matrix of the finite element is:

kel=∫VelBTDB dV\htmlClass{fea-math-bold-italic}{\boldsymbol{k}}_{el}=\int_{V_{el}}\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}^{T}\htmlClass{fea-math-bold-italic}{\boldsymbol{D}}\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}\,dV

With EE denoting Young’s modulus and ν\nu Poisson’s ratio:

D=E(1−ν)(1+ν)(1−2ν)[1ν1−νν1−ν000ν1−ν1ν1−ν000ν1−νν1−ν10000001−2ν2(1−ν)0000001−2ν2(1−ν)0000001−2ν2(1−ν)]\htmlClass{fea-math-bold-italic}{\boldsymbol{D}}=\frac{E(1-\nu)}{(1+\nu)(1-2\nu)}\begin{bmatrix}1&\dfrac{\nu}{1-\nu}&\dfrac{\nu}{1-\nu}&0&0&0\\\dfrac{\nu}{1-\nu}&1&\dfrac{\nu}{1-\nu}&0&0&0\\\dfrac{\nu}{1-\nu}&\dfrac{\nu}{1-\nu}&1&0&0&0\\0&0&0&\dfrac{1-2\nu}{2(1-\nu)}&0&0\\0&0&0&0&\dfrac{1-2\nu}{2(1-\nu)}&0\\0&0&0&0&0&\dfrac{1-2\nu}{2(1-\nu)}\end{bmatrix}

Because B\htmlClass{fea-math-bold-italic}{\boldsymbol{B}} is a function of s,t,rs,t,r, it is more convenient to work in natural coordinates. Consequently:

kel=∫−11∫−11∫−11B(s,t,r)TDB(s,t,r)det ⁣(J(s,t,r)) ds dt dr\htmlClass{fea-math-bold-italic}{\boldsymbol{k}}_{el}=\int_{-1}^{1}\int_{-1}^{1}\int_{-1}^{1}\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}(s,t,r)^{T}\htmlClass{fea-math-bold-italic}{\boldsymbol{D}}\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}(s,t,r)\mathop{\mathit{det}}\!\bigl(\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}(s,t,r)\bigr)\,ds\,dt\,dr

det ⁣(J(s,t,r))\mathop{\mathit{det}}\!\bigl(\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}(s,t,r)\bigr) is the determinant of the Jacobian J\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}.

06

Gauss quadrature

For second-order Gauss integration:

kel≃∑i=12∑j=12∑k=12E ⁣((−1)i33,(−1)j33,(−1)k33)\htmlClass{fea-math-bold-italic}{\boldsymbol{k}}_{el}\simeq\sum_{i=1}^{2}\sum_{j=1}^{2}\sum_{k=1}^{2}\htmlClass{fea-math-bold-italic}{\boldsymbol{E}}\!\left((-1)^i\frac{\sqrt{3}}{3},(-1)^j\frac{\sqrt{3}}{3},(-1)^k\frac{\sqrt{3}}{3}\right)

Here E\htmlClass{fea-math-bold-italic}{\boldsymbol{E}} is the integrand:

E(s,t,r)=B(s,t,r)TDB(s,t,r)det ⁣(J(s,t,r))\htmlClass{fea-math-bold-italic}{\boldsymbol{E}}(s,t,r)=\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}(s,t,r)^T\htmlClass{fea-math-bold-italic}{\boldsymbol{D}}\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}(s,t,r)\mathop{\mathit{det}}\!\bigl(\htmlClass{fea-math-bold-italic}{\boldsymbol{J}}(s,t,r)\bigr)

07

MATLAB program

A simple MATLAB program can be downloaded below. Run main.m.

MMATLAB — HEX8 linear analysisSeven MATLAB files and README; default displacement scale 1 and stress scale 0.Download ZIP

08

Example

Cantilever beam: L=200 mmL=200\,\mathit{mm}, cross-section 50 mm×40 mm50\,\mathit{mm}\times40\,\mathit{mm}, E=2×105 MPaE=2\times10^5\,\mathit{MPa}, ν=0.3\nu=0.3, F=8000 NF=8000\,\mathit{N}.

Mesh: 60×16×12=11 52060\times16\times12=11\,520 HEX8 elements, 13 481 nodes and 40 443 equations.

Cantilever mesh, constraints and loading
Figure 2. Cantilever mesh, constraints and loading. All nodes at the clamped end are fixed in the xx-direction. The nodes at the clamped end with z=h/2z=h/2 (on the horizontal axis of symmetry) are also fixed in the yy- and zz-directions.

Vertical displacement ww:

QuantityMATLABANSYS
wmax≈0w_{\mathop{\mathit{max}}}\approx00.0002311 mm0.0002311\,\mathit{mm}0.000176 mm0.000176\,\mathit{mm}
∣w∣max|w|_{\mathop{\mathit{max}}}0.41216 mm0.41216\,\mathit{mm}0.413512 mm0.413512\,\mathit{mm}

The tip displacement is negative: wtip=−0.41216 mmw_{\mathit{tip}}=-0.41216\,\mathit{mm} in MATLAB and −0.413512 mm-0.413512\,\mathit{mm} in ANSYS.

MATLAB — vertical displacement w [mm]

ANSYS — vertical displacement w [mm]

Figure 3. Vertical displacement ww [mm]: MATLAB (left) and ANSYS (right). The deformation display scale is 120.

Axial stress σx\sigma_x:

QuantityMATLABANSYS
σx,min\sigma_{x,\mathop{\mathit{min}}}−119.242 MPa−118.092 MPa
σx,max\sigma_{x,\mathop{\mathit{max}}}119.242 MPa118.093 MPa

MATLAB — axial stress σx [MPa]

ANSYS — axial stress σx [MPa]

Figure 4. Axial stress σx\sigma_x [MPa]: MATLAB (left) and ANSYS (right). The deformation display scale is 120.

Classical beam verification:

σmax=MW=8000⋅20050⋅402/6=120 MPa.\sigma_{\mathop{\mathit{max}}}=\frac{M}{W}=\frac{8000\cdot200}{50\cdot40^2/6}=120\,\mathit{MPa}.

Taking into account both bending and shear deformation:

∣w∣max=FL33EI+FLkGA,|w|_{\mathop{\mathit{max}}}=\frac{FL^3}{3EI}+\frac{FL}{kGA},

with:

I=50⋅40312,k=56,G=E2(1+ν),I=\frac{50\cdot40^3}{12},\qquad k=\frac56,\qquad G=\frac{E}{2(1+\nu)},

which gives approximately:

∣w∣max=0.4125 mm.|w|_{\mathop{\mathit{max}}}=0.4125\,\mathit{mm}.

The ANSYS macro can be downloaded below.

AANSYS — Cantilever exampleCantilever comparison macro.Download MAC

09

References

  1. Carlos A. Felippa, Introduction to Finite Element Methods, University of Colorado, Boulder, Colorado.