Finite element stiffness matrices

Constitutive equations

Constitutive equations relate deformation to force. For an infinitesimal element constitutive equations relate stress to strain, σ=σ(ϵ). For a solid with no pre-stress, i.e., σ(0)=0, the linear truncation of the MacLaurin series of σ(ϵ) yields σ=C ϵ, with C a rank-4 tensor in general.

Hookean elasticity (linear, isotropic)

The material is characterized by Young modulus Y, Poisson ratio ν, and C reduces to a rank-2 tensor (matrix) H

H=Simplify[(Y/(1-(ν^2))){{1,ν,0},{ν,1,0},{0,0,((1-ν)/2)}}];
MatrixForm[H]

Hexahedron element

This computation follows http://what-when-how.com/the-finite-element-method/fem-for-3d-solids-finite-element-method-part-2/

Define the 8 nodes of hexahedral element in the physical coordinate system X=(x,y,z)

nodesX = Table[{(Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])},{i,8}]

Define the node coordinates in the local, transformed coordinate system Ξ=(ξ,η,ζ)

nodesΞ = Table[2IntegerDigits[i,2,3]-1,{i,0,7}]

Define the form functions

Ξ={ξ,η,ζ};
Nform[ξ_,η_,ζ_] = Table[Product[(1/2)(1+nodesΞ[[i,j]]Ξ[[j]]),{j,3}],{i,8}]

Define the nodal displacements

δ=Table[{(Subscript[u, i]),(Subscript[v, i]),(Subscript[z, i])},{i,8}]

Define the element displacement interpolation

U[ξ_,η_,ζ_] = Sum[ Nform[ξ,η,ζ][[i]] δ[[i]],{i,8}]

Also define the displacement interpolation in terms of a single vector argument, and check the element node numbering

U[Ξ_]:=U[Ξ[[1]],Ξ[[2]],Ξ[[3]]]

Table[U[nodesΞ[[j]]],{j,8}]

Define the form function matrix

Nblock[Ni_]:=IdentityMatrix[3] Ni;
NformM[ξ_,η_,ζ_]=Transpose[Flatten[Map[Nblock[#]&,Nform[ξ,η,ζ]],1]];
Dimensions[NformM[ξ,η,ζ]]

Define the gradient in transformed coordinates

gradC[f_]:={D[f,ξ],D[f,η],D[f,ζ]}

Compute the Jacobian matrix

J[ξP_,ηP_,ζP_]:=gradC[Nform[ξ,η,ζ]] . nodesX /. {ξ->ξP,η->ηP,ζ->ζP};

MatrixForm[J[ξ,η,ζ]]

Define the gradient in physical coordinates. Note that this function should be evaluated for numeric arguments (ξP,ηP,ζP) to avoid costly analytical evaluation of the Jacobian inverse.

gradP[f_,ξP_,ηP_,ζP_]:=Inverse[J[ξP,ηP,ζP]] . {D[f,ξ],D[f,η],D[f,ζ]}

Define the strain operator

L[Ni_,ξP_,ηP_,ζP_] := Outer[Times,{1,1,1},gradP[Ni,ξP,ηP,ζP]]

{{

?DiagonalMatrix

Outer[Times,{1,1,1},{a,b,c}]

MatrixForm[%]

Finite element stiffness matrix