diff --git a/docs/theory/examples/cartesian_linear-elasticity.py b/docs/theory/examples/cartesian_linear-elasticity.py new file mode 100644 index 0000000..2606cb6 --- /dev/null +++ b/docs/theory/examples/cartesian_linear-elasticity.py @@ -0,0 +1,238 @@ + +import numpy as np + +np.set_printoptions(linewidth=200) + +# ================================================================================================== +# identity tensors +# ================================================================================================== + +II = np.zeros((3,3,3,3)) +I4 = np.zeros((3,3,3,3)) +I4rt = np.zeros((3,3,3,3)) +Is = np.zeros((3,3,3,3)) + +for i in range(3): + for j in range(3): + for k in range(3): + for l in range(3): + if ( i == j and k == l ): + II[i,j,k,l] = 1. + +for i in range(3): + for j in range(3): + for k in range(3): + for l in range(3): + if ( i == l and j == k ): + I4[i,j,k,l] = 1. + +for i in range(3): + for j in range(3): + for k in range(3): + for l in range(3): + if ( i == k and j == l ): + I4rt[i,j,k,l] = 1. + +I4s = ( I4 + I4rt ) / 2. +I4d = ( I4s - II/3. ) + +# ================================================================================================== +# elasticity tensor +# ================================================================================================== + +# bulk and shear modulus +K = 0.8333333333333333 +G = 0.3846153846153846 + +# elasticity tensor (3d) +C4 = K * II + 2. * G * I4d + +# ================================================================================================== +# mesh definition +# ================================================================================================== + +# number of elements +nx = 10 +ny = 10 + +# mesh dimensions +nelem = nx * ny # number of elements +nnode = (nx+1) * (ny+1) # number of nodes +nne = 4 # number of nodes per element +nd = 2 # number of dimensions +ndof = nnode * nd # total number of degrees of freedom + +# coordinates and connectivity: zero-initialize +coor = np.zeros((nnode,nd ), dtype='float') +conn = np.zeros((nelem,nne), dtype='int' ) + +# coordinates: set +# - grid of points +x,y = np.meshgrid(np.linspace(0,1,nx+1), np.linspace(0,1,ny+1)) +# - store from grid of points +coor[:,0] = x.ravel() +coor[:,1] = y.ravel() + +# connectivity: set +# - grid of node numbers +inode = np.arange(nnode).reshape(ny+1, nx+1) +# - store from grid of node numbers +conn[:,0] = inode[ :-1, :-1].ravel() +conn[:,1] = inode[ :-1,1: ].ravel() +conn[:,2] = inode[1: ,1: ].ravel() +conn[:,3] = inode[1: , :-1].ravel() + +# DOF-numbers per node +dofs = np.arange(nnode*nd).reshape(nnode,nd) + +# ================================================================================================== +# quadrature definition +# ================================================================================================== + +# integration point coordinates (local element coordinates) +Xi = np.array([ + [-1./np.sqrt(3.), -1./np.sqrt(3.)], + [+1./np.sqrt(3.), -1./np.sqrt(3.)], + [+1./np.sqrt(3.), +1./np.sqrt(3.)], + [-1./np.sqrt(3.), +1./np.sqrt(3.)], +]) + +# integration point weights +W = np.array([ + [1.], + [1.], + [1.], + [1.], +]) + +# ================================================================================================== +# assemble stiffness matrix +# ================================================================================================== + +# allocate matrix +K = np.zeros((ndof, ndof)) + +# loop over nodes +for e in conn: + + # - nodal coordinates + xe = coor[e,:] + + # - allocate element stiffness matrix + Ke = np.zeros((nne*nd, nne*nd)) + + # - numerical quadrature + for w, xi in zip(W, Xi): + + # -- shape function gradients (w.r.t. the local element coordinates) + dNdxi = np.array([ + [-.25*(1.-xi[1]), -.25*(1.-xi[0])], + [+.25*(1.-xi[1]), -.25*(1.+xi[0])], + [+.25*(1.+xi[1]), +.25*(1.+xi[0])], + [-.25*(1.+xi[1]), +.25*(1.-xi[0])], + ]) + + # -- Jacobian + J = np.zeros((nd, nd)) + + for m in range(nne): + for i in range(nd): + for j in range(nd): + J[i,j] += dNdxi[m,i] * xe[m,j] + + Jinv = np.linalg.inv(J) + Jdet = np.linalg.det(J) + + # -- shape function gradients (w.r.t. the global coordinates) + dNdx = np.zeros((nne,nd)) + + for m in range(nne): + for i in range(nd): + for j in range(nd): + dNdx[m,i] += Jinv[i,j] * dNdxi[m,j] + + # -- assemble to element stiffness matrix + for m in range(nne): + for n in range(nne): + for i in range(nd): + for j in range(nd): + for k in range(nd): + for l in range(nd): + Ke[m*nd+j, n*nd+k] += dNdx[m,i] * C4[i,j,k,l] * dNdx[n,l] * w * Jdet + + # - assemble to global stiffness matrix + iie = dofs[e,:].ravel() + + K[np.ix_(iie,iie)] += Ke + +# ================================================================================================== +# partition and solve +# ================================================================================================== + +# zero-initialize displacements and forces +f = np.zeros((ndof)) +u = np.zeros((ndof)) + +# fixed displacements +# - zero-initialize +iip = np.empty((0),dtype='int' ) +up = np.empty((0),dtype='float') +# - 'y = 0' : 'u_y = 0' +idof = dofs[inode[0,:], 1] +iip = np.hstack(( iip, idof )) +up = np.hstack(( up , 0.0 * np.ones(idof.shape) )) +# - 'x = 0' : 'u_x = 0' +idof = dofs[inode[:,0], 0] +iip = np.hstack(( iip, idof )) +up = np.hstack(( up , 0.0 * np.ones(idof.shape) )) +# - 'x = 1' : 'u_x = 0.1' +idof = dofs[inode[:,-1], 0] +iip = np.hstack(( iip, idof )) +up = np.hstack(( up , 0.1 * np.ones(idof.shape) )) + +# free displacements +iiu = np.setdiff1d(dofs.ravel(), iip) + +# partition +# - stiffness matrix +Kuu = K[np.ix_(iiu, iiu)] +Kup = K[np.ix_(iiu, iip)] +Kpu = K[np.ix_(iip, iiu)] +Kpp = K[np.ix_(iip, iip)] +# - external force +fu = f[iiu] + +# solve +uu = np.linalg.solve(Kuu, fu - Kup.dot(up)) +fp = Kpu.dot(uu) + Kpp.dot(up) + +# assemble +u[iiu] = uu +u[iip] = up +f[iip] = fp + +# convert to nodal displacement and forces +disp = u[dofs] +fext = f[dofs] + +# ================================================================================================== +# plot +# ================================================================================================== + +print(fext[inode[:,-1],0].sum()) + +import matplotlib.pyplot as plt + +fig, axes = plt.subplots(ncols=2, figsize=(16,8)) + +axes[0].scatter(coor[:,0], coor[:,1]) +axes[0].quiver (coor[:,0], coor[:,1], disp[:,0], disp[:,1]) + +axes[1].scatter((coor+disp)[:,0], (coor+disp)[:,1]) +axes[1].quiver ((coor+disp)[:,0], (coor+disp)[:,1], fext[:,0], fext[:,1], scale=.1) + +for axis in axes: + axis.set_xlim([-0.4, 1.5]) + axis.set_ylim([-0.4, 1.5]) + +plt.show() diff --git a/docs/theory/examples/cylindrical_linear-elasticity.py b/docs/theory/examples/cartesian_linear-elasticity_B-matrix.py similarity index 52% copy from docs/theory/examples/cylindrical_linear-elasticity.py copy to docs/theory/examples/cartesian_linear-elasticity_B-matrix.py index d788089..46f290e 100644 --- a/docs/theory/examples/cylindrical_linear-elasticity.py +++ b/docs/theory/examples/cartesian_linear-elasticity_B-matrix.py @@ -1,212 +1,195 @@ import numpy as np np.set_printoptions(linewidth=200) # ================================================================================================== # tensor products # ================================================================================================== def ddot43(A,B): nd = A.shape[0] C = np.zeros((nd,nd,nd)) for i in range(nd): for j in range(nd): for k in range(nd): for l in range(nd): for m in range(nd): C[i,j,m] += A[i,j,k,l] * B[l,k,m] return C # -------------------------------------------------------------------------------------------------- def ddot33(A,B): nd = A.shape[0] C = np.zeros((nd,nd)) for i in range(nd): for j in range(nd): for k in range(nd): for l in range(nd): C[i,l] += A[i,j,k] * B[k,j,l] return C +# -------------------------------------------------------------------------------------------------- + +def transpose3(A): + + nd = A.shape[0] + + C = np.zeros((nd,nd,nd)) + + for i in range(nd): + for j in range(nd): + for k in range(nd): + C[k,j,i] = A[i,j,k] + + return C + # ================================================================================================== # identity tensors # ================================================================================================== II = np.zeros((3,3,3,3)) I4 = np.zeros((3,3,3,3)) I4rt = np.zeros((3,3,3,3)) Is = np.zeros((3,3,3,3)) for i in range(3): for j in range(3): for k in range(3): for l in range(3): if ( i == j and k == l ): - II[i,i,j,l] += 1. + II[i,j,k,l] = 1. for i in range(3): for j in range(3): for k in range(3): for l in range(3): if ( i == l and j == k ): - I4[i,i,j,l] += 1. + I4[i,j,k,l] = 1. for i in range(3): for j in range(3): for k in range(3): for l in range(3): if ( i == k and j == l ): - I4rt[i,i,j,l] += 1. + I4rt[i,j,k,l] = 1. I4s = ( I4 + I4rt ) / 2. I4d = ( I4s - II/3. ) # ================================================================================================== # elasticity tensor # ================================================================================================== -K = 1. -G = 1. +K = 0.8333333333333333 +G = 0.3846153846153846 C4 = K * II + 2. * G * I4d # ================================================================================================== # mesh definition # ================================================================================================== nne = 4 ndim = 2 xe = np.array([ [0., 0.], [1., 0.], [1., 1.], [0., 1.], ]) Xi = np.array([ [-1./np.sqrt(3.), -1./np.sqrt(3.)], [+1./np.sqrt(3.), -1./np.sqrt(3.)], [+1./np.sqrt(3.), +1./np.sqrt(3.)], [-1./np.sqrt(3.), +1./np.sqrt(3.)], ]) W = np.array([ [1.], [1.], [1.], [1.], ]) Ke = np.zeros((nne*ndim, nne*ndim)) for w, xi in zip(W, Xi): - N = np.array([ - [.25 * (1.-xi[0]) * (1.-xi[1])], - [.25 * (1.+xi[0]) * (1.-xi[1])], - [.25 * (1.+xi[0]) * (1.+xi[1])], - [.25 * (1.-xi[0]) * (1.+xi[1])], - ]) - dNdxi = np.array([ - [-.25 * (1.-xi[1]), -.25 * (1.-xi[0])], - [+.25 * (1.-xi[1]), -.25 * (1.+xi[0])], - [+.25 * (1.+xi[1]), +.25 * (1.+xi[0])], - [-.25 * (1.+xi[1]), +.25 * (1.-xi[0])], + [-.25*(1.-xi[1]), -.25*(1.-xi[0])], + [+.25*(1.-xi[1]), -.25*(1.+xi[0])], + [+.25*(1.+xi[1]), +.25*(1.+xi[0])], + [-.25*(1.+xi[1]), +.25*(1.-xi[0])], ]) J = np.zeros((ndim, ndim)) - for n in range(nne): + for m in range(nne): for i in range(ndim): for j in range(ndim): - J[i,j] += dNdxi[n,i] * xe[n,j] + J[i,j] += dNdxi[m,i] * xe[m,j] Jinv = np.linalg.inv(J) Jdet = np.linalg.det(J) dNdx = np.zeros((nne,ndim)) - for n in range(nne): + for m in range(nne): for i in range(ndim): for j in range(ndim): - dNdx[n,i] += Jinv[i,j] * dNdxi[n,j] - - xk = np.zeros((ndim)) - - for n in range(nne): - for i in range(ndim): - xk[i] += N[n] * xe[n,i] + dNdx[m,i] += Jinv[i,j] * dNdxi[m,j] - rk = xk[0] + for m in range(nne): - for n in range(nne): + Bm = np.zeros((3,3,3)) - Bn = np.zeros((3,3,3)) + Bm[0,0,0] = dNdx[m,0] + Bm[0,1,1] = dNdx[m,0] - Bn[0,0,0] = dNdx[n,0] # B(n, r r r ) - Bn[0,1,0] = -1./rk * N[n] # B(n, r \theta r ) - Bn[0,1,1] = dNdx[n,0] # B(n, r \theta \theta ) - Bn[0,2,2] = dNdx[n,0] # B(n, r z z ) + Bm[1,0,0] = dNdx[m,1] + Bm[1,1,1] = dNdx[m,1] - Bn[1,0,0] = 0. # B(n, \theta r r ) - axisymmetric - Bn[1,1,0] = +1./rk * N[n] # B(n, \theta \theta r ) - Bn[1,1,1] = 0. # B(n, \theta \theta \theta ) - axisymmetric - Bn[1,2,2] = 0. # B(n, \theta z z ) - axisymmetric + for n in range(nne): - Bn[2,0,0] = dNdx[n,1] # B(n, z r r ) - Bn[2,1,1] = dNdx[n,1] # B(n, z \theta \theta ) - Bn[2,2,2] = dNdx[n,1] # B(n, z z z ) + Bn = np.zeros((3,3,3)) - for m in range(nne): + Bn[0,0,0] = dNdx[n,0] + Bn[0,1,1] = dNdx[n,0] - Bm = np.zeros((3,3,3)) + Bn[1,0,0] = dNdx[n,1] + Bn[1,1,1] = dNdx[n,1] - Bm[0,0,0] = dNdx[m,0] # B(n, r r r ) - Bm[0,1,0] = -1./rk * N[n] # B(n, r \theta r ) - Bm[0,1,1] = dNdx[m,0] # B(n, r \theta \theta ) - Bm[0,2,2] = dNdx[m,0] # B(n, r z z ) + Kmn = ddot33(transpose3(Bm),ddot43(C4,Bn)) - Bm[1,0,0] = 0. # B(n, \theta r r ) - axisymmetric - Bm[1,1,0] = +1./rk * N[n] # B(n, \theta \theta r ) - Bm[1,1,1] = 0. # B(n, \theta \theta \theta ) - axisymmetric - Bm[1,2,2] = 0. # B(n, \theta z z ) - axisymmetric + iim = np.array([ m*ndim+0, m*ndim+1 ]) + iin = np.array([ n*ndim+0, n*ndim+1 ]) - Bm[2,0,0] = dNdx[m,1] # B(n, z r r ) - Bm[2,1,1] = dNdx[m,1] # B(n, z \theta \theta ) - Bm[2,2,2] = dNdx[m,1] # B(n, z z z ) - - K = ddot33(Bn,ddot43(C4,Bm)) - - iin = np.array([ 2*n+0, 2*n+1 ]) - iim = np.array([ 2*m+0, 2*m+1 ]) - - Ke[np.ix_(iin,iim)] += K[np.ix_([0,2], [0,2])] - - # apply the integration point volume - Ke *= w * Jdet * rk + Ke[np.ix_(iim,iin)] += Kmn[np.ix_([0,1], [0,1])] * w * Jdet print(Ke) +# print(np.linalg.inv(Ke)) + # Ke[] diff --git a/docs/theory/examples/cylindrical_linear-elasticity.py b/docs/theory/examples/cylindrical_linear-elasticity.py index d788089..7cbe7bc 100644 --- a/docs/theory/examples/cylindrical_linear-elasticity.py +++ b/docs/theory/examples/cylindrical_linear-elasticity.py @@ -1,212 +1,224 @@ import numpy as np np.set_printoptions(linewidth=200) # ================================================================================================== # tensor products # ================================================================================================== def ddot43(A,B): nd = A.shape[0] C = np.zeros((nd,nd,nd)) for i in range(nd): for j in range(nd): for k in range(nd): for l in range(nd): for m in range(nd): C[i,j,m] += A[i,j,k,l] * B[l,k,m] return C # -------------------------------------------------------------------------------------------------- def ddot33(A,B): nd = A.shape[0] C = np.zeros((nd,nd)) for i in range(nd): for j in range(nd): for k in range(nd): for l in range(nd): C[i,l] += A[i,j,k] * B[k,j,l] return C +# -------------------------------------------------------------------------------------------------- + +def transpose3(A): + + nd = A.shape[0] + + C = np.zeros((nd,nd,nd)) + + for i in range(nd): + for j in range(nd): + for k in range(nd): + C[k,j,i] = A[i,j,k] + + return C + # ================================================================================================== # identity tensors # ================================================================================================== II = np.zeros((3,3,3,3)) I4 = np.zeros((3,3,3,3)) I4rt = np.zeros((3,3,3,3)) Is = np.zeros((3,3,3,3)) for i in range(3): for j in range(3): for k in range(3): for l in range(3): if ( i == j and k == l ): - II[i,i,j,l] += 1. + II[i,j,k,l] = 1. for i in range(3): for j in range(3): for k in range(3): for l in range(3): if ( i == l and j == k ): - I4[i,i,j,l] += 1. + I4[i,j,k,l] = 1. for i in range(3): for j in range(3): for k in range(3): for l in range(3): if ( i == k and j == l ): - I4rt[i,i,j,l] += 1. + I4rt[i,j,k,l] = 1. I4s = ( I4 + I4rt ) / 2. I4d = ( I4s - II/3. ) # ================================================================================================== # elasticity tensor # ================================================================================================== -K = 1. -G = 1. +K = 0.8333333333333333 +G = 0.3846153846153846 C4 = K * II + 2. * G * I4d # ================================================================================================== # mesh definition # ================================================================================================== nne = 4 ndim = 2 xe = np.array([ [0., 0.], [1., 0.], [1., 1.], [0., 1.], ]) Xi = np.array([ [-1./np.sqrt(3.), -1./np.sqrt(3.)], [+1./np.sqrt(3.), -1./np.sqrt(3.)], [+1./np.sqrt(3.), +1./np.sqrt(3.)], [-1./np.sqrt(3.), +1./np.sqrt(3.)], ]) W = np.array([ [1.], [1.], [1.], [1.], ]) Ke = np.zeros((nne*ndim, nne*ndim)) for w, xi in zip(W, Xi): N = np.array([ [.25 * (1.-xi[0]) * (1.-xi[1])], [.25 * (1.+xi[0]) * (1.-xi[1])], [.25 * (1.+xi[0]) * (1.+xi[1])], [.25 * (1.-xi[0]) * (1.+xi[1])], ]) dNdxi = np.array([ - [-.25 * (1.-xi[1]), -.25 * (1.-xi[0])], - [+.25 * (1.-xi[1]), -.25 * (1.+xi[0])], - [+.25 * (1.+xi[1]), +.25 * (1.+xi[0])], - [-.25 * (1.+xi[1]), +.25 * (1.-xi[0])], + [-.25*(1.-xi[1]), -.25*(1.-xi[0])], + [+.25*(1.-xi[1]), -.25*(1.+xi[0])], + [+.25*(1.+xi[1]), +.25*(1.+xi[0])], + [-.25*(1.+xi[1]), +.25*(1.-xi[0])], ]) J = np.zeros((ndim, ndim)) - for n in range(nne): + for m in range(nne): for i in range(ndim): for j in range(ndim): - J[i,j] += dNdxi[n,i] * xe[n,j] + J[i,j] += dNdxi[m,i] * xe[m,j] Jinv = np.linalg.inv(J) Jdet = np.linalg.det(J) dNdx = np.zeros((nne,ndim)) - for n in range(nne): + for m in range(nne): for i in range(ndim): for j in range(ndim): - dNdx[n,i] += Jinv[i,j] * dNdxi[n,j] + dNdx[m,i] += Jinv[i,j] * dNdxi[m,j] xk = np.zeros((ndim)) for n in range(nne): for i in range(ndim): xk[i] += N[n] * xe[n,i] rk = xk[0] - for n in range(nne): - - Bn = np.zeros((3,3,3)) + for m in range(nne): - Bn[0,0,0] = dNdx[n,0] # B(n, r r r ) - Bn[0,1,0] = -1./rk * N[n] # B(n, r \theta r ) - Bn[0,1,1] = dNdx[n,0] # B(n, r \theta \theta ) - Bn[0,2,2] = dNdx[n,0] # B(n, r z z ) + Bm = np.zeros((3,3,3)) - Bn[1,0,0] = 0. # B(n, \theta r r ) - axisymmetric - Bn[1,1,0] = +1./rk * N[n] # B(n, \theta \theta r ) - Bn[1,1,1] = 0. # B(n, \theta \theta \theta ) - axisymmetric - Bn[1,2,2] = 0. # B(n, \theta z z ) - axisymmetric + Bm[0,0,0] = dNdx[m,0] # B(m, r r r ) + Bm[0,1,0] = -1./rk * N[m] # B(m, r \theta r ) + Bm[0,1,1] = dNdx[m,0] # B(m, r \theta \theta ) + Bm[0,2,2] = dNdx[m,0] # B(m, r z z ) - Bn[2,0,0] = dNdx[n,1] # B(n, z r r ) - Bn[2,1,1] = dNdx[n,1] # B(n, z \theta \theta ) - Bn[2,2,2] = dNdx[n,1] # B(n, z z z ) + Bm[1,0,0] = 0. # B(m, \theta r r ) - axisymmetric + Bm[1,1,0] = +1./rk * N[m] # B(m, \theta \theta r ) + Bm[1,1,1] = 0. # B(m, \theta \theta \theta ) - axisymmetric + Bm[1,2,2] = 0. # B(m, \theta z z ) - axisymmetric - for m in range(nne): + Bm[2,0,0] = dNdx[m,1] # B(m, z r r ) + Bm[2,1,1] = dNdx[m,1] # B(m, z \theta \theta ) + Bm[2,2,2] = dNdx[m,1] # B(m, z z z ) - Bm = np.zeros((3,3,3)) + for n in range(nne): - Bm[0,0,0] = dNdx[m,0] # B(n, r r r ) - Bm[0,1,0] = -1./rk * N[n] # B(n, r \theta r ) - Bm[0,1,1] = dNdx[m,0] # B(n, r \theta \theta ) - Bm[0,2,2] = dNdx[m,0] # B(n, r z z ) + Bn = np.zeros((3,3,3)) - Bm[1,0,0] = 0. # B(n, \theta r r ) - axisymmetric - Bm[1,1,0] = +1./rk * N[n] # B(n, \theta \theta r ) - Bm[1,1,1] = 0. # B(n, \theta \theta \theta ) - axisymmetric - Bm[1,2,2] = 0. # B(n, \theta z z ) - axisymmetric + Bn[0,0,0] = dNdx[n,0] # B(n, r r r ) + Bn[0,1,0] = -1./rk * N[n] # B(n, r \theta r ) + Bn[0,1,1] = dNdx[n,0] # B(n, r \theta \theta ) + Bn[0,2,2] = dNdx[n,0] # B(n, r z z ) - Bm[2,0,0] = dNdx[m,1] # B(n, z r r ) - Bm[2,1,1] = dNdx[m,1] # B(n, z \theta \theta ) - Bm[2,2,2] = dNdx[m,1] # B(n, z z z ) + Bn[1,0,0] = 0. # B(n, \theta r r ) - axisymmetric + Bn[1,1,0] = +1./rk * N[n] # B(n, \theta \theta r ) + Bn[1,1,1] = 0. # B(n, \theta \theta \theta ) - axisymmetric + Bn[1,2,2] = 0. # B(n, \theta z z ) - axisymmetric - K = ddot33(Bn,ddot43(C4,Bm)) + Bn[2,0,0] = dNdx[n,1] # B(n, z r r ) + Bn[2,1,1] = dNdx[n,1] # B(n, z \theta \theta ) + Bn[2,2,2] = dNdx[n,1] # B(n, z z z ) - iin = np.array([ 2*n+0, 2*n+1 ]) - iim = np.array([ 2*m+0, 2*m+1 ]) + Kmn = ddot33(transpose3(Bm),ddot43(C4,Bn)) - Ke[np.ix_(iin,iim)] += K[np.ix_([0,2], [0,2])] + iim = np.array([ m*ndim+0, m*ndim+1 ]) + iin = np.array([ n*ndim+0, n*ndim+1 ]) - # apply the integration point volume - Ke *= w * Jdet * rk + Ke[np.ix_(iim,iin)] += Kmn[np.ix_([0,2], [0,2])] * w * Jdet * rk print(Ke) # Ke[] diff --git a/docs/theory/readme.tex b/docs/theory/readme.tex index 93b2c99..373e410 100644 --- a/docs/theory/readme.tex +++ b/docs/theory/readme.tex @@ -1,1094 +1,1228 @@ %!TEX program = xelatex \documentclass[times,namecite]{goose-article} \usepackage{lipsum,verbatim,mdframed,array} \title{% The Finite Element Method in Continuum Mechanics } \author{T.W.J.~de~Geus} \contact{% $^*$Contact: % \href{http://goosefem.geus.me}{goosefem.geus.me}% \hspace{1mm}--\hspace{1mm} % \href{mailto:tom@geus.me}{tom@geus.me} % \hspace{1mm}--\hspace{1mm} % \href{http://www.geus.me}{www.geus.me}% } \hypersetup{pdfauthor={T.W.J. de Geus}} \header{% The Finite Element Method in Continuum Mechanics, T.W.J.\ de Geus } \begin{document} \maketitle \begin{abstract} The theory of the Finite Element Method in Continuum Mechanics is discussed in a compact way. This discussion is by no means comprehensive and one is invited to dive in more complete textbooks, but also to get one's hands dirty by implementing some (simple) examples. \end{abstract} \setcounter{tocdepth}{2} \tableofcontents \section{Statics} \subsection{The basic concept} The discussion starts by considering a static, solid mechanics, problem. Loosely speaking the goal is to find a deformation map, $\vec{x} = \varphi(\vec{X},t)$, that maps a body $\Omega_0$ to a deformed state $\Omega$ that satisfies equilibrium and the boundary conditions applied on $\Gamma$. This is illustrated in Fig.~\ref{fig:problem}, whereby it is emphasized that the body can be subjected to two kinds of boundary conditions: \begin{itemize} \item \emph{Essential} or \emph{Dirichlet} boundary conditions on $\Gamma_p$, whereby the displacements are prescribed. \item \emph{Natural} or \emph{Neumann} boundary conditions on $\Gamma_u$, whereby the tractions are prescribed. Note that traction-free natural boundary conditions are also perfectly acceptable. \end{itemize} \begin{figure}[htp] \centering \includegraphics[width=.5\textwidth]{figures/problem.pdf} \caption{Generic problem definition} \label{fig:problem} \end{figure} In practice, one is not looking for the map $\vec{x} = \varphi(\vec{X},t)$, but for the deformation gradient $\bm{F}$, or in fact for a displacement field $\vec{u} = \vec{x} - \vec{X}$. To make things a bit more explicit, the deformation gradient is defined as follows: \begin{equation} \vec{x} = \bm{F} \cdot \vec{X} \end{equation} hence \begin{equation} \bm{F} = \frac{\partial \varphi}{\partial \vec{X}} = \big( \vec{\nabla}_0 \, \vec{x} \big)^T = \bm{I} + \big( \vec{\nabla}_0 \, \vec{u} \big)^T \end{equation} \subsection{Momentum balance} The momentum balance reads \begin{equation} \label{eq:static:momentum} \vec{\nabla} \cdot \bm{\sigma}(\vec{x}) = \vec{0} \qquad \vec{x} \in \Omega \end{equation} where $\bm{\sigma}$ is the Cauchy stress which depends on the new position $\vec{x}$ and thus on the displacement $\vec{u}$. It has been assumed that all actions are instantaneous (no inertia) and, for simplicity, that there are no body forces. Loosely speaking the interpretation of this equation is that \emph{the sum of all forces vanishes} everywhere in the domain $\Omega$. The crux of the Finite Element Method is that this non-linear differential equation is solved in a weak sense. I.e. \begin{equation} \label{eq:static:int} \int\limits_\Omega \vec{\phi}(\vec{x}) \cdot \big[\, \vec{\nabla} \cdot \bm{\sigma}(\vec{x}) \,\big] \; \mathrm{d}\Omega = 0 \qquad \forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d \end{equation} where $\vec{\phi}$ are test functions. For reasons that become obvious below, integration by parts is applied, which results in \begin{equation} \label{eq:static:weak} \int\limits_\Omega \big[\, \vec{\nabla} \vec{\phi}(\vec{x}) \,\big] : \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega = \int\limits_\Omega \vec{\nabla} \cdot \big[\, \vec{\phi}(\vec{x}) \cdot \bm{\sigma}(\vec{x}) \,\big] \; \mathrm{d}\Omega \qquad \forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d \end{equation} In going from Eq.~\eqref{eq:static:int} to Eq.~\eqref{eq:static:weak}, use has been made of the following chain rule \begin{equation} \vec{\nabla} \cdot \big[\, \vec{\phi} \cdot \bm{\sigma}^T \,\big] = \big[\, \vec{\nabla} \vec{\phi} \,\big] : \bm{\sigma}^T + \vec{\phi} \cdot \big[\, \vec{\nabla} \cdot \bm{\sigma} \,\big] \end{equation} together with the symmetry of the Cauchy stress \begin{equation} \bm{\sigma} = \bm{\sigma}^T \end{equation} The right-hand side of this equation can be reduced to an area integral by employing Gauss's divergence theorem: \begin{equation} \int\limits_\Omega \vec{\nabla} \cdot \vec{a}(\vec{x}) \; \mathrm{d}\Omega = \int\limits_\Gamma \vec{n}(\vec{x}) \cdot \vec{a}(\vec{x}) \; \mathrm{d}\Gamma \end{equation} where $\vec{n}$ is the normal along the surface $\Gamma$. The result reads \begin{equation} \label{eq:static:weak:final} \int\limits_\Omega \big[\, \vec{\nabla} \vec{\phi}(\vec{x}) \,\big] : \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega = \int\limits_\Gamma \vec{\phi}(\vec{x}) \cdot \underbrace{ \vec{n}(\vec{x}) \cdot \bm{\sigma}(\vec{x}) }_{ \vec{t}(\vec{x}) } \; \mathrm{d}\Gamma \qquad \forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d \end{equation} \subsection{Discretization} The problem is now discretized using $n$ \emph{nodes} that are connected through \emph{elements}, which define the discretized domain $\Omega^h_0$. Shape functions $N_i(\vec{x})$ are used to extrapolate the nodal quantities throughout the domain $\Omega^h_0$ (and $\Omega^h$). For example the displacement field $\vec{u}(\vec{x})$ \begin{equation} \vec{u}(\vec{x},t) \approx \vec{u}^h(\vec{x},t) = - \sum_{i=1}^{n} N_i (\vec{x}) \; \vec{u}_i (t) - = - \underline{N}^\mathsf{T} (\vec{x}) \; \underline{\vec{x}} (t) + \sum_{m=1}^{n_\mathrm{nodes}} N_m (\vec{x}) \; \vec{u}_m (t) \end{equation} whereby we will not further consider the time dependence (denoted with $t$) for now. Following standard Galerkin the test functions are interpolated in the same way, i.e.\ \begin{equation} \label{eq:discretization} \vec{\phi}(\vec{x}) \approx \vec{\phi}^h(\vec{x}) = - \sum_{i=1}^{n} N_i (\vec{x}) \; \vec{\phi}_i - = - \underline{N}^\mathsf{T} (\vec{x}) \; \underline{\vec{\phi}} + \sum_{m=1}^{n_\mathrm{nodes}} N_m (\vec{x}) \; \vec{\phi}_m \end{equation} Applied to the problem sketch from Fig.~\ref{fig:problem}, a discretization might look like Fig.~\ref{fig:problem:discretized}. The nodes are clearly marked as circles. The lines connecting the nodes clearly mark the elements which are in this case three-node linear triangles. \begin{figure}[htp] \centering \includegraphics[width=.25\textwidth]{figures/problem-discretized.pdf} \caption{Example of a discretization of the problem in Fig.~\ref{fig:problem} using Finite Elements (linear triangles in this case).} \label{fig:problem:discretized} \end{figure} Applied to the balance equation the following is obtained \begin{equation} - \underline{\vec{\phi}}^\mathsf{T} \cdot + \vec{\phi}_m \cdot \int\limits_{\Omega^h} - \big[\, \vec{\nabla} \underline{N}(\vec{x}) \,\big] + \big[\, \vec{\nabla} N_m (\vec{x}) \,\big] \cdot \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega = - \underline{\vec{\phi}}^\mathsf{T} \cdot + \vec{\phi}_m \cdot \int\limits_{\Gamma^h} - \underline{N}(\vec{x}) \cdot + N_m (\vec{x}) \cdot \vec{t}(\vec{x}) \; \mathrm{d}\Gamma \qquad - \forall \; \underline{\vec{\phi}} \in \mathbb{R}^d_n + \forall \; \vec{\phi}_m \in \mathbb{R}^d_n \end{equation} -from which the dependency on $\underline{\vec{\phi}}$ can be dropped: +from which the dependency on $\vec{\phi}_m$ can be dropped: \begin{equation} \int\limits_{\Omega^h} - \big[\, \vec{\nabla} \underline{N}(\vec{x}) \,\big] + \big[\, \vec{\nabla} N_m (\vec{x}) \,\big] \cdot \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega = \int\limits_{\Gamma^h} - \underline{N}(\vec{x}) \cdot + N_m (\vec{x}) \cdot \vec{t}(\vec{x}) \; \mathrm{d}\Gamma \end{equation} This corresponds to a (non-linear) set of nodal balance equations: \begin{equation} - \underline{\vec{f}}^\mathrm{int}(\vec{x}) + \vec{f}_m^\mathrm{int}(\vec{x}) = - \underline{\vec{f}}^\mathrm{ext}(\vec{x}) + \vec{f}_m^\mathrm{ext}(\vec{x}) \end{equation} with: \begin{itemize} \item \emph{Internal forces} \begin{equation} - \underline{\vec{f}}^\mathrm{int}(\vec{x}) + \vec{f}_m^\mathrm{int}(\vec{x}) = \int\limits_{\Omega^h} - \big[\, \vec{\nabla} \underline{N}(\vec{x}) \,\big] + \big[\, \vec{\nabla} N_m (\vec{x}) \,\big] \cdot \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega \end{equation} \item \emph{External forces} \begin{equation} - \underline{\vec{f}}^\mathrm{ext}(\vec{x}) + \vec{f}_m^\mathrm{ext}(\vec{x}) = \int\limits_{\Gamma^h} - \underline{N}(\vec{x}) \cdot + N_m (\vec{x}) \cdot \vec{t}(\vec{x}) \; \mathrm{d}\Gamma \end{equation} Note that this term is zero in the interior of the domain, i.e. in $\Omega^h \bigcap \Gamma^h$, while it can be zero or non-zero in $\Gamma^h$ depending on the problem details. Most often it following completely from the boundary conditions of the boundary value problem that we set out to solve. \end{itemize} \subsection{Iterative solution -- small strain} A commonly used strategy to solve the non-linear system is the iterative Newton-Raphson scheme (see Appendix~\ref{sec:newton-raphson}). The idea is thereby to formulate an initial guess for the solution, determine possible residual forces, and use these forces to come to a better guess for the solution. This is continued until the solution has been found, i.e. when the residual vanishes. This solution technique is discussed here in the context of small deformations. Assuming the deformations (and therefore rotations) to be small leads to the assumption that $\Omega = \Omega_0$, and thus that $\nabla = \nabla_0$. In this context one typically uses the linear strain tensor \begin{equation} \bm{\varepsilon} = \tfrac{1}{2} \left[ \nabla_0 \vec{u} + \big[\, \nabla_0 \vec{u} \,\big]^T \right] = \mathbb{I}_s : \big[\, \nabla_0 \vec{u} \,\big] \end{equation} and some (non-linear) relationship between it and the stress as follows \begin{equation} \bm{\sigma} = \bm{\sigma} \big( \bm{\varepsilon} \big) \end{equation} To further simplify the discussion, the boundary tractions are assumed to be some displacement independent quantity, which is \textit{a-priori} known for all relevant boundary nodes. The nodal equilibrium equations now reads \begin{equation} \label{eq:static:eq-non-lin} - \underline{\vec{r}}(\vec{X}, \vec{u}) + \vec{r}_m(\vec{X}, \vec{u}) = - \underline{\vec{f}}^\mathrm{ext}(\vec{X}) + \vec{f}_m^\mathrm{ext}(\vec{X}) - - \underline{\vec{f}}^\mathrm{int}(\vec{X}, \vec{u}) + \vec{f}_m^\mathrm{int}(\vec{X}, \vec{u}) = \underline{\vec{0}} \end{equation} with \begin{equation} - \underline{\vec{f}}^\mathrm{int}(\vec{X}, \vec{u}) + \vec{f}_m^\mathrm{int}(\vec{X}, \vec{u}) = \int\limits_{\Omega^h_0} - \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big] + \big[\, \vec{\nabla}_0 N_m(\vec{X}) \,\big] \cdot \bm{\sigma}(\vec{X}, \vec{u}) \; \mathrm{d}\Omega_0 \end{equation} The displacement field $\vec{u}$ is now iteratively updated starting from some initial guess, i.e. \begin{equation} \vec{u}_{(i+1)} = \vec{u}_{(i)} + \delta \vec{u} \end{equation} The update, $\delta \vec{u}$, follows from the linearization of Eq.~\eqref{eq:static:eq-non-lin}. The results of which is \begin{equation} \int\limits_{\Omega^h_0} - \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big] + \big[\, \vec{\nabla}_0 N_m(\vec{X}) \,\big] \cdot \mathbb{K}\big(\vec{X},\vec{u}_{(i)}\big) \cdot - \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big]^\mathsf{T} \; + \big[\, \vec{\nabla}_0 N_n(\vec{X}) \,\big] \; \mathrm{d}\Omega_0 - \cdot \delta \underline{\vec{u}} + \cdot \delta \vec{u}_n = - \underline{\vec{f}}^\mathrm{ext}(\vec{X}) + \vec{f}_m^\mathrm{ext}(\vec{X}) - \int\limits_{\Omega^h_0} - \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big] + \big[\, \vec{\nabla}_0 N_m(\vec{X}) \,\big] \cdot \bm{\sigma}\big(\vec{X},\vec{u}_{(i)}\big) \; \mathrm{d}\Omega_0 \end{equation} where, at each position $\vec{X}$, \begin{equation} \mathbb{K}\big(\vec{u}_{(i)}\big) = \left. \frac{\partial \bm{\sigma}}{\partial \bm{\varepsilon}} \right|_{\vec{u}_{(i)}} : \mathbb{I}_s \end{equation} is the \emph{constitutive tangent operator}. In a shorter notation, the iterative update reads: \begin{equation} - \underline{\underline{\bm{K}}}_{(i)} \cdot \delta \underline{\vec{u}} + \bm{K}_{mn,(i)} \cdot \delta \vec{u}_m = - \underline{\vec{f}}^\mathrm{ext} + \vec{f}_m^\mathrm{ext} - - \underline{\vec{f}}^\mathrm{int}_{(i)} + \vec{f}_{m,(i)}^\mathrm{int} \end{equation} with \begin{equation} - \underline{\underline{\bm{K}}}_{(i)} + \bm{K}_{mn,(i)} = \int\limits_{\Omega^h_0} - \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big] + \big[\, \vec{\nabla}_0 N_m(\vec{X}) \,\big] \cdot \mathbb{K}\big(\vec{X}, \vec{u}_{(i)}\big) \cdot - \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big]^\mathsf{T} \; + \big[\, \vec{\nabla}_0 N_m(\vec{X}) \,\big] \; \mathrm{d}\Omega_0 \end{equation} and \begin{equation} - \underline{\vec{f}}^\mathrm{int}_{(i)} + \vec{f}_{m,(i)}^\mathrm{int} = \int\limits_{\Omega^h_0} - \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big] + \big[\, \vec{\nabla}_0 N_m(\vec{X}) \,\big] \cdot \bm{\sigma}\big(\vec{X}, \vec{u}_{(i)}\big) \; \mathrm{d}\Omega_0 \end{equation} \section{Dynamics} \subsection{Momentum balance} The next step is to add inertia to the balance of momentum in Eq.~\eqref{eq:static:momentum}: \begin{equation} \rho\, \vec{a} = \vec{\nabla} \cdot \bm{\sigma}(\vec{x}) \qquad \vec{x} \in \Omega \end{equation} where $\rho$ is the mass density, and second time derivative of the position $\vec{x}$ is the acceleration $\vec{a} = \ddot{\vec{u}} = \ddot{\vec{x}}$. Note that this function is the continuum equivalent of $\vec{f} = m \vec{a}$. Like before, this equation is solved in a weak sense \begin{equation} \int\limits_\Omega \rho(\vec{x})\; \vec{\phi}(\vec{x}) \cdot \vec{a} \; \mathrm{d}\Omega = \int\limits_\Omega \vec{\phi}(\vec{x}) \cdot \Big[\, \vec{\nabla} \cdot \bm{\sigma}(\vec{x}) \,\Big] \; \mathrm{d}\Omega \qquad \forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d \end{equation} Integration by parts results in \begin{equation} \int\limits_\Omega \rho(\vec{x})\; \vec{\phi}(\vec{x}) \cdot \vec{a} \; \mathrm{d}\Omega = \int\limits_\Gamma \vec{\phi}(\vec{x}) \cdot \vec{t}(\vec{x}) \; \mathrm{d}\Gamma - \int\limits_\Omega \big[\, \vec{\nabla} \vec{\phi}(\vec{x}) \,\big] : \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega \qquad \forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d \end{equation} which is discretized as before: \begin{equation} - \underline{\vec{\phi}}^\mathsf{T} \cdot + \vec{\phi}_m \cdot \int\limits_\Omega - \rho(\vec{x})\; \underline{N}(\vec{x})\; \underline{N}^\mathsf{T}(\vec{x}) \; + \rho(\vec{x})\; N_m(\vec{x})\; N_n(\vec{x}) \; \mathrm{d}\Omega \; - \underline{\vec{a}} + \vec{a}_n = - \underline{\vec{\phi}}^\mathsf{T} \cdot + \vec{\phi}_m \cdot \int\limits_\Gamma - \underline{N}(\vec{x})\; \vec{t}(\vec{x}) \; + N_m(\vec{x})\; \vec{t}(\vec{x}) \; \mathrm{d}\Gamma - - \underline{\vec{\phi}}^\mathsf{T} \cdot + \vec{\phi}_m \cdot \int\limits_\Omega - \big[\, \vec{\nabla} \underline{N}(\vec{x}) \,\big] + \big[\, \vec{\nabla} N_m(\vec{x}) \,\big] : \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega \qquad - \forall \; \underline{\vec{\phi}} \in \mathbb{R}^d_n + \forall \; \vec{\phi}_m \in \mathbb{R}^d_n \end{equation} This is independent of the test functions, hence: \begin{equation} \label{eq:dynamics:system} \underbrace{ \int\limits_\Omega - \rho(\vec{x})\; \underline{N}(\vec{x})\; \underline{N}^\mathsf{T}(\vec{x}) \; + \rho(\vec{x})\; N_m(\vec{x})\; N_n(\vec{x}) \; \mathrm{d}\Omega - }_{\underline{\underline{M}}(\vec{x})} \; - \underline{\vec{a}} + }_{\displaystyle + M_{mn}(\vec{x}) + } \; + \vec{a}_n = \underbrace{ \int\limits_\Gamma - \underline{N}(\vec{x})\; \vec{t}(\vec{x}) \; + N_m(\vec{x})\; \vec{t}(\vec{x}) \; \mathrm{d}\Gamma - }_{\underline{\vec{f}}^\mathrm{ext}(\vec{x})} + }_{\displaystyle + \vec{f}_m^\mathrm{ext}(\vec{x}) + } - \underbrace{ \int\limits_\Omega - \big[\, \vec{\nabla} \underline{N}(\vec{x}) \,\big] + \big[\, \vec{\nabla} N_m(\vec{x}) \,\big] : \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega - }_{\underline{\vec{f}}^\mathrm{int}(\vec{x})} + }_{\displaystyle + \vec{f}_m^\mathrm{int}(\vec{x}) + } \end{equation} This is denoted as follows \begin{equation} - \underline{\underline{M}}(\vec{x})\; \underline{\vec{a}} + M_{mn}(\vec{x})\; \vec{a}_n = - \underline{\vec{f}}^\mathrm{ext}(\vec{x}) + \vec{f}_m^\mathrm{ext}(\vec{x}) - - \underline{\vec{f}}^\mathrm{int}(\vec{x}) + \vec{f}_m^\mathrm{int}(\vec{x}) \end{equation} -Where $\underline{\underline{M}}(\vec{x})$ is the \emph{mass matrix}, $\underline{\vec{f}}^\mathrm{ext}(\vec{x})$ are the \emph{external forces}, and $\underline{\vec{f}}^\mathrm{int}(\vec{x})$ are the \emph{internal forces}. +Where $M_{mn}(\vec{x})$ is the \emph{mass matrix}, $\vec{f}_m^\mathrm{ext}(\vec{x})$ are the \emph{external forces}, and $\vec{f}_m^\mathrm{int}(\vec{x})$ are the \emph{internal forces}. \subsection{Time discretization} \label{eq:dynamics:time} To solve the second order differential equation in time, one typically also discretizes time with some finite difference based scheme. For this we also need to distinguish the velocity field $\vec{v} = \dot{\vec{u}} = \dot{\vec{x}}$. \subsubsection{Verlet} In the case that the forces are velocity independent one can use the Velocity Verlet scheme, in which energy is conserved, i.e.\ it is reversible in time. The scheme reads are follows: \begin{enumerate} \item Compute the velocity on a dummy time grid: \begin{equation} \vec{v}_{n+1/2} = \vec{v}_{n-1/2} + \Delta_t \; \vec{a}_n \end{equation} Note that this involves solving Eq.~\eqref{eq:dynamics:system} for $\vec{a}_n$. From this it is directly obvious why the forces need to be velocity independent. \item Update the positions \begin{equation} \vec{u}_{n+1} = \vec{u}_n + \Delta_t \;\vec{a} \vec{v}_{n + 1/2} \end{equation} \end{enumerate} \subsubsection{Velocity Verlet} In the case that the forces are dependent on the velocity (i.e.\ when there is damping) the previous scheme cannot be used anymore. An additional approximation has to be made in this case. Note that this is not very problematic as energy is in no case conserved. The resulting scheme reads: \begin{enumerate} \item Compute the position at $t_{n+1} = t_{n} + \Delta_t$: \begin{equation} \vec{u}_{n+1} = \vec{u}_{n} + \Delta_t \vec{v}_{n} + \tfrac{1}{2} \Delta_t^2 \vec{a}_{n} \end{equation} \item Estimate the velocity at $t_{n+1} = t_{n} + \Delta_t$ (involves solving Eq.~\eqref{eq:dynamics:system}): \begin{equation} \hat{\vec{v}}_{n+1} = \vec{v}_{n} + \tfrac{1}{2} \Delta_t \Big[\, \vec{a}_{n} + \vec{a} ( \vec{u}_{n+1} , \vec{v}_{n} + \Delta_t \vec{a}_{n} , t_{n+1} ) \, \Big] \end{equation} \item Correct $\hat{\vec{v}}_{n+1}$ (involves solving Eq.~\eqref{eq:dynamics:system}): \begin{equation} \vec{v}_{n+1} = \vec{v}_{n} + \tfrac{1}{2} \Delta_t \Big[\, \vec{a}_{n} + \vec{a} ( \vec{u}_{n+1} , \hat{\vec{v}}_{n+1} , t_{n+1} ) \, \Big] \end{equation} \end{enumerate} \subsection{Approximations} For problems where the local volume is conversed (either weakly through slightly compressible elasticity, or strongly by adding an incompressibility constraint) it makes sense to assume the mass matrix constant, as any change of volume results in an equivalent change of the density. In that case \begin{equation} \int\limits_{\Omega} \rho(\vec{x}) \;\mathrm{d}\Omega = \int\limits_{\Omega_0} \rho(\vec{X}) \;\mathrm{d}\Omega_0 \end{equation} Which results in: \begin{equation} - \underline{\underline{M}}(\vec{x}) + M_{mn}(\vec{x}) = \int\limits_{\Omega_0} - \rho(\vec{X})\; \underline{N}(\vec{X})\; \underline{N}^\mathsf{T}(\vec{X}) \; + \rho(\vec{X})\; N_m(\vec{X})\; N_n(\vec{X}) \; \mathrm{d}\Omega_0 = \mathrm{constant} \end{equation} To enhance computational efficiency, it may be a good option concentrate the mass to point masses on the nodes. This has to strong advantage that the mass matrix becomes diagonal, known as the \emph{lumped mass matrix}. Consequently, instead of solving a linear system one just has to solve fully independent equations. \appendix \section{Nomenclature} \begin{itemize} \item Vector \begin{equation} \vec{u} = \sum\limits_i u_i \vec{e}_i \end{equation} for example for Cartesian coordinates \begin{equation} \vec{u} = u_x \vec{e}_x + u_y \vec{e}_y + u_z \vec{e}_z \end{equation} % \item Second-order tensor \begin{equation} \bm{A} = \sum\limits_i \sum\limits_j A_{ij} \vec{e}_i \vec{e}_j \end{equation} % \item Fourth-order tensor \begin{equation} \mathbb{A} = \sum\limits_i \sum\limits_j \sum\limits_k \sum\limits_l A_{ijkl} \vec{e}_i \vec{e}_j \vec{e}_k \vec{e}_l \end{equation} % \item Divergence \begin{equation} \vec{\nabla} \cdot \bm{\sigma} = \frac{ \partial \sigma_{ij} }{ \partial x_i } \end{equation} % \item Double tensor contraction \begin{equation} C = \bm{A} : \bm{B} = A_{ij} B_{ji} \end{equation} \end{itemize} \section{Shape functions} In the Finite Element Method a geometry is discretized using nodes. The nodes are grouped in elements which define the domain $\Omega^h$. The crux of the method is that nodal quantities, for example $\vec{u}_i$, are interpolated throughout the discretized domain $\Omega^h$ using shape functions $N_i (\vec{x})$. The shape functions are polynomials that spanned by the nodes, such that they constitute a partition of unity and satisfy $N_i (\vec{x}_j) = \delta_{ij}$, i.e. it is one in the node $i$, and zero in all other nodes. This implies that, while each shape function is globally supported, $N_i (\vec{x}) \neq 0$ only in the elements containing node $i$. For all other elements it is zero (even in between the nodes). This again implies that $\partial N_i / \partial \vec{x}$ is discontinuous across element boundaries. For a one-dimensional problem comprising four linear elements and five nodes the shape functions are sketched in Fig.~\ref{fig:shape:1d-domain}. To emphasize the global support of each shape function the shape function of node 2 is isolated in Fig.~\ref{fig:shape:1d-domain-2}. As can be seen, node 2 is only non-zero in elements 1 and 2, while it is zero everywhere else. To evaluate for example $f_2$, the internal force one node 2, the integration can be limited to these two elements: \begin{equation} f_2 = \int\limits_{\Omega^{(1)}} \frac{\mathrm{d} N^{(1)}_2}{\mathrm{d} x} \; \sigma (x) \; \mathrm{d}\Omega^{(1)} + \int\limits_{\Omega^{(2)}} \frac{\mathrm{d} N^{(2)}_2}{\mathrm{d} x} \; \sigma (x) \; \mathrm{d}\Omega^{(2)} \end{equation} \begin{figure}[htp] \centering \captionsetup[subfigure]{justification=centering} \begin{minipage}[t]{.5\textwidth} \centering \includegraphics[width=1.\textwidth]{figures/shape-functions-1d.pdf} \subcaption{All shape functions $N_i(x)$, each with a different color.} \label{fig:shape:1d-domain} \end{minipage} \\ \begin{minipage}[t]{.5\textwidth} \centering \includegraphics[width=1.\textwidth]{figures/shape-functions-1d-node-2.pdf} \subcaption{Isolated shape function $N_2(x)$.} \label{fig:shape:1d-domain-2} \end{minipage} \caption{(a) Shape functions for a one-dimensional domain discretized using four two node, linear, element. (b) Isolated shape functions for node 2. The node numbers are displayed using a unique color, while the element numbers are shown in black in a different font, in between the nodes.} \label{fig:shape:1d-domain:fig} \end{figure} By now it should be clear that the above allows us assemble $\underline{f}$ element-by-element. For the example in Fig.~\ref{fig:shape:1d-domain:fig} this is graphically shown in Fig.~\ref{fig:shape:1d-domain:assembly}. \begin{figure}[htp] \centering \captionsetup[subfigure]{justification=centering} \begin{minipage}[t]{.5\textwidth} \centering \includegraphics[width=1.\textwidth]{figures/shape-functions-1d-element-0.pdf} \end{minipage} \\ $+$ \\ \begin{minipage}[t]{.5\textwidth} \centering \includegraphics[width=1.\textwidth]{figures/shape-functions-1d-element-1.pdf} \end{minipage} \\ $+$ \\ \begin{minipage}[t]{.5\textwidth} \centering \includegraphics[width=1.\textwidth]{figures/shape-functions-1d-element-2.pdf} \end{minipage} \\ $+$ \\ \begin{minipage}[t]{.5\textwidth} \centering \includegraphics[width=1.\textwidth]{figures/shape-functions-1d-element-3.pdf} \end{minipage} \caption{Element-by-element assembly of the problem in Fig.~\ref{fig:shape:1d-domain:fig}.} \label{fig:shape:1d-domain:assembly} \end{figure} where the indices show that the \emph{shape functions} are evaluated compared to some generic element definition. \section{Iso-parametric transformation and quadrature} A very important concept in the Finite Element Method is the iso-parametric transformation. It allows an arbitrarily shaped element with volume $\Omega^e$ to be mapped onto a generic \emph{iso-parametric element} of constant volume $Q$, see Fig.~\ref{fig:isoparametric}. Using this mapping it is easy to perform interpolation and numerical quadrature using a generic implementation. \begin{figure}[htp] \centering \includegraphics[width=.5\textwidth]{figures/isoparametric-transform.pdf} \caption{Iso-parameteric transformation for a four node, bi-linear, quadrilateral element.} \label{fig:isoparametric} \end{figure} The mapping between the generic domain $Q$ and the physical domain $\Omega^e$ is as follows \begin{equation} \label{eq:isoparametric:coordinate} - \vec{x} ( \vec{\xi} ) = \big[\, \underline{N}^{e} \,\big]^\mathsf{T} \underline{x}^e + \vec{x} ( \vec{\xi} ) = N_m^{e} \; \underline{x}_m^e \end{equation} where the column $\underline{x}^e$ contains the position vectors of the element nodes in the actual coordinates (i.e.\ spanning $\Omega^e$). In order to perform the quadrature on $Q$ we must also map the gradient operator: \begin{equation} \vec{\nabla}_{\xi}\, = \vec{e}_i \frac{\partial}{\partial \xi_i} = \vec{e}_i \frac{\partial x_j(\vec{\xi})}{\partial \xi_i} \frac{\partial}{\partial x_j} = \vec{e}_i \frac{\partial x_j(\vec{\xi})}{\partial \xi_i} \vec{e}_j \cdot \vec{e}_k \frac{\partial}{\partial x_k} = \big[\, \vec{\nabla}_{\xi}\, \vec{x}(\vec{\xi}) \,\big] \cdot \vec{\nabla} = \bm{J}(\vec{\xi}) \cdot \vec{\nabla} \end{equation} i.e. \begin{equation} J_{ij} = \frac{\partial x_i}{\partial \xi_j} \end{equation} from which the inverse relationship trivially follows \begin{equation} \vec{\nabla} = \bm{J}^{-1}(\vec{\xi}) \cdot \vec{\nabla}_{\xi}\, \end{equation} where $\bm{J}$ is the Jacobian which can be obtained from the gradient of the shape functions with respect to the iso-parametric coordinates and the nodal coordinates in the actual coordinates through Eq.~\eqref{eq:isoparametric:coordinate}: \begin{equation} \bm{J}(\vec{\xi}) = \vec{\nabla}_{\xi}\, \vec{x}(\vec{\xi}) = - \big[\, \vec{\nabla}_{\xi}\, \underline{N}^{e} \,\big]^\mathsf{T} \; \underline{x}^e + \vec{\nabla}_{\xi}\, \big[\, N_m^{e} \; \vec{x}_m^e \,\big] + = + \big[\, \vec{\nabla}_{\xi}\, N_m^{e} \,\big] \; \vec{x}_m^e \end{equation} The gradient of the shape functions with respect to the actual coordinates can now be computed though \begin{equation} - \vec{\nabla} \underline{N}^{e} + \vec{\nabla} N_m^{e} = - \bm{J}^{-1}(\vec{\xi}) \cdot \big[\, \vec{\nabla}_{\xi}\, \underline{N}^{e} \,\big] + \bm{J}^{-1}(\vec{\xi}) \cdot \big[\, \vec{\nabla}_{\xi}\, N_m^{e} \,\big] \end{equation} -Furthermore, the mapping between the real and the master volume is equal to + +The change of variables should also be accounted for in the integration: +\begin{equation} + \int_\Omega (\ldots) \; \mathrm{d} \Omega + = + \iiint (\ldots) \; \mathrm{d} x \, \mathrm{d} y \, \mathrm{d} z + = + \iiint (\ldots) \; \left| \frac{\partial \left( x, y, z \right)}{\partial \left( \xi, \eta, \mu \right)} \right| \mathrm{d} \xi \, \mathrm{d} \eta \, \mathrm{d} \mu +\end{equation} +where the determinant of the Jacobian +\begin{equation} + \frac{\partial \left( x, y, z \right)}{\partial \left( \xi, \eta, \mu \right)} + = + \left| + \begin{array}{ccc} + \displaystyle \frac{\partial x}{\partial \xi } & + \displaystyle \frac{\partial x}{\partial \eta} & + \displaystyle \frac{\partial x}{\partial \mu } + \\ + \displaystyle \frac{\partial y}{\partial \xi } & + \displaystyle \frac{\partial y}{\partial \eta} & + \displaystyle \frac{\partial y}{\partial \mu } + \\ + \displaystyle \frac{\partial z}{\partial \xi } & + \displaystyle \frac{\partial z}{\partial \eta} & + \displaystyle \frac{\partial z}{\partial \mu } + \end{array} + \right| +\end{equation} +We can now make use of the Jacobian that we already obtained +\begin{equation} + \int_\Omega (\ldots) \; \mathrm{d} \Omega + = + \iiint (\ldots) \det \left( \bm{J}(\vec{\xi}) \right) \; \mathrm{d} \xi \, \mathrm{d} \eta \, \mathrm{d} \mu + = + \int_Q (\ldots) \det \left( \bm{J}(\vec{\xi}) \right) \; \mathrm{d} Q +\end{equation} +Note that another way to obtain this result is using \begin{align} \mathrm{d} \Omega &= \mathrm{d} \vec{x}_0 \times \mathrm{d} \vec{x}_1 \cdot \mathrm{d} \vec{x}_2 = \left[ \mathrm{d} \vec{x}_0 \cdot \bm{J}(\vec{\xi}) \right] \times \left[ \mathrm{d} \vec{x}_1 \cdot \bm{J}(\vec{\xi}) \right] \cdot \left[ \mathrm{d} \vec{x}_2 \cdot \bm{J}(\vec{\xi}) \right] = \det \big( \bm{J}(\vec{\xi}) \big)\, \mathrm{d} \vec{\xi}_0 \times \mathrm{d} \vec{\xi}_1 \cdot \mathrm{d} \vec{\xi}_2 \\ &= \det \big( \bm{J}(\vec{\xi}) \big)\, \mathrm{d} Q \end{align} Numerical quadrature can be formulated (exactly) on the master element. It corresponds to taking the weighted sum of the integrand evaluated at specific \emph{quadrature points} (or \emph{integration points}). The practical implications, for example for the internal force, are as follows: \begin{align} - \underline{\vec{f}^e} + \vec{f}_m^e &= \int\limits_{\Omega^e} - \big[\, \vec{\nabla} \underline{N} \,\big] + \big[\, \vec{\nabla} N_m \,\big] \cdot \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega = \int\limits_{Q} - \big[\, \vec{\nabla} \underline{N} \,\big] + \big[\, \vec{\nabla} N_m \,\big] \cdot \bm{\sigma}(\vec{x}) \; \det \big( \bm{J}(\vec{\xi}) \big) \; \mathrm{d}Q \\ &= \sum_{k}^{n_k} w_k - \big[\, \vec{\nabla} \underline{N} \,\big]_{\vec{\xi} = \vec{\xi}_k} + \big[\, \vec{\nabla} N_m \,\big]_{\vec{\xi} = \vec{\xi}_k} \cdot \bm{\sigma}\big(\vec{x}(\vec{\xi}_k)\big) \; \det \big( \bm{J}(\vec{\xi}_k) \big) \; \end{align} \subsection*{Total Lagrange} To obtain \begin{equation} \vec{x}(\vec{\xi}),\, \vec{\nabla}_0,\,\, \text{and}\, \int\limits_{\Omega_0} . \;\mathrm{d}\Omega \end{equation} -simply replace $\underline{x}^e$ with $\underline{X}^e$ in Eq.~\eqref{eq:isoparametric:coordinate}. For this reason the same element implementation can be used in small strain and finite strain (total Lagrange and updated Lagrange), by providing either $\underline{X}^e$ or $\underline{X}^e + \underline{u}^e$ as input. +simply replace $\vec{x}_m^e$ with $\vec{X}_m^e$ in Eq.~\eqref{eq:isoparametric:coordinate}. For this reason the same element implementation can be used in small strain and finite strain (total Lagrange and updated Lagrange), by providing either $\vec{X}_m^e$ or $\vec{X}_m^e + \vec{u}_m^e$ as input. \section{Newton-Raphson in one dimension} \label{sec:newton-raphson} The objective is to find $x$ such that \begin{equation} r(x) = 0 \end{equation} An initial guess is made for $x$ which is (hopefully) iteratively improved. This iterative value is denoted using $x_{(i)}$. It is updated by making use of the following Taylor expansion \begin{equation} r \big( x_{(i+1)} \big) = r \big( x_{(i)} \big) + \left. \frac{\mathrm{d} r}{\mathrm{d} x} \right|_{x = x_{(i)}} \delta x + \mathcal{O} \big( \delta x^2 \big) \approx 0 \end{equation} where \begin{equation} \delta x = x_{(i+1)} - x_{(i)} \end{equation} The value of $\delta x$ can be obtained by neglecting higher order terms: \begin{equation} r \big( x_{(i)} \big) + \left. \frac{\mathrm{d} r}{\mathrm{d} x} \right|_{x = x_{(i)}} \delta x = 0 \end{equation} From which it follows that \begin{equation} \delta x = - \left[ \left. \frac{\mathrm{d} r}{\mathrm{d} x} \right|_{x = x_{(i)}} \right]^{-1} r \big( x_{(i)} \big) \end{equation} Thereafter, the improved `guess' of the solution is \begin{equation} x_{(i+1)} = x_{(i)} + \delta x \end{equation} These steps are then repeated, until the solution has been reached within a certain accuracy $\epsilon$: \begin{equation} \left| r \big( x_{(i+1)} \big) \right| < \epsilon \end{equation} The iterative scheme is well understood from the illustration in Fig.~\ref{fig:newton-raphson}. \begin{figure}[htp] \centering \includegraphics[width=.3\textwidth]{figures/newton-raphson.pdf} \caption{Schematic of the Newton-Raphson iterations in a one-dimensional problem.} \label{fig:newton-raphson} \end{figure} +\section{B-matrix} + +\subsection{Discretization} + +\eqref{eq:static:weak:final} + +\begin{equation} + \int\limits_\Omega + \big[\, \vec{\nabla} \vec{\phi}(\vec{x}) \,\big] : \bm{\sigma}(\vec{x}) \; + \mathrm{d}\Omega + = + \int\limits_\Gamma + \vec{\phi}(\vec{x}) \cdot + \vec{t}(\vec{x}) \; + \mathrm{d}\Gamma + \qquad + \forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d +\end{equation} + +The discretization from Eq.~\eqref{eq:discretization} can be written more precisely as + +\begin{equation} + \phi_x(\vec{x}) \vec{e}_x + + \phi_y(\vec{x}) \vec{e}_y + + \phi_z(\vec{x}) \vec{e}_z + \approx + \sum_{m=1}^{n_\mathrm{nodes}} + N_m (\vec{x}) \; \phi_{m,x} \vec{e}_x + + N_m (\vec{x}) \; \phi_{m,y} \vec{e}_y + + N_m (\vec{x}) \; \phi_{m,z} \vec{e}_z +\end{equation} + + +\begin{align} + & + \vec{\nabla} \vec{\phi}(\vec{x}) + \approx + \sum_{m=1}^{n_\mathrm{nodes}} + \bigg( + \vec{e}_x \frac{\partial}{\partial x} + + \vec{e}_y \frac{\partial}{\partial y} + + \vec{e}_z \frac{\partial}{\partial z} + \bigg) + \bigg( + N_m (\vec{x}) \; \phi_{m,x} \vec{e}_x + + N_m (\vec{x}) \; \phi_{m,y} \vec{e}_y + + N_m (\vec{x}) \; \phi_{m,z} \vec{e}_z + \bigg) + \\ + & + = + \phi_{m,x} \frac{\partial N_m}{\partial x} \vec{e}_x \vec{e}_x + + \phi_{m,y} \frac{\partial N_m}{\partial x} \vec{e}_x \vec{e}_y + + \phi_{m,z} \frac{\partial N_m}{\partial x} \vec{e}_x \vec{e}_z + \\ + & + + + \phi_{m,x} \frac{\partial N_m}{\partial y} \vec{e}_y \vec{e}_x + + \phi_{m,y} \frac{\partial N_m}{\partial y} \vec{e}_y \vec{e}_y + + \phi_{m,z} \frac{\partial N_m}{\partial y} \vec{e}_y \vec{e}_z + \\ + & + + + \phi_{m,x} \frac{\partial N_m}{\partial z} \vec{e}_z \vec{e}_x + + \phi_{m,y} \frac{\partial N_m}{\partial z} \vec{e}_z \vec{e}_y + + \phi_{m,z} \frac{\partial N_m}{\partial z} \vec{e}_z \vec{e}_z +\end{align} + +Which we can denote in short as +\begin{equation} + \vec{\nabla} \vec{\phi} (\vec{x}) + \approx + \sum_{m=1}^{n_\mathrm{nodes}} + \mathcal{B}_m (\vec{x}) \cdot \vec{\phi}_m +\end{equation} + + +\begin{alignat}{3} + \mathcal{B}_{m,x x x} &= \frac{\partial N_m}{\partial x} & \qquad + \mathcal{B}_{m,y x x} &= \frac{\partial N_m}{\partial y} & \qquad + \mathcal{B}_{m,z x x} &= \frac{\partial N_m}{\partial z} + \nonumber + \\ + \mathcal{B}_{m,x y y} &= \frac{\partial N_m}{\partial x} & \qquad + \mathcal{B}_{m,y y y} &= \frac{\partial N_m}{\partial y} & \qquad + \mathcal{B}_{m,z y y} &= \frac{\partial N_m}{\partial z} + \\ + \mathcal{B}_{m,x z z} &= \frac{\partial N_m}{\partial x} & \qquad + \mathcal{B}_{m,y z z} &= \frac{\partial N_m}{\partial y} & \qquad + \mathcal{B}_{m,z z z} &= \frac{\partial N_m}{\partial z} +\end{alignat} + +Note that similarly +\begin{equation} + \vec{\nabla} \vec{u} (\vec{x}) + \approx + \sum_{m=1}^{n_\mathrm{nodes}} + \mathcal{B}_m (\vec{x}) \cdot \vec{u}_m +\end{equation} + + +\begin{equation} + \int\limits_{\Omega^h} + \left( \mathcal{B}_m \cdot \vec{\phi}_m \right) : \bm{\sigma}(\vec{x}) \; + \mathrm{d}\Omega + = + \int\limits_{\Gamma^h} + N_m (\vec{x}) \; \vec{\phi}_m \cdot \vec{t}(\vec{x}) \; + \mathrm{d}\Gamma + \qquad + \forall \; \vec{\phi}_m \in \mathbb{R}^d_n +\end{equation} + + +\begin{equation} + \vec{\phi}_m \cdot + \int\limits_{\Omega^h} + \mathcal{B}_m^T : \bm{\sigma}(\vec{x}) \; + \mathrm{d}\Omega + = + \vec{\phi}_m \cdot + \int\limits_{\Gamma^h} + N_m(\vec{x}) \; \vec{t}(\vec{x}) \; + \mathrm{d}\Gamma + \qquad + \forall \; \vec{\phi}_m \in \mathbb{R}^d_n +\end{equation} + + +\begin{equation} + \int\limits_{\Omega^h} + \mathcal{B}_m^T : \bm{\sigma}(\vec{x}) \; + \mathrm{d}\Omega + = + \int\limits_{\Gamma^h} + N_m(\vec{x}) \; \vec{t}(\vec{x}) \; + \mathrm{d}\Gamma +\end{equation} + +\subsection{Linear elasticity} + +\begin{equation} + \int\limits_{\Omega_0^h} + \tilde{\mathcal{B}}_m^T : \Big( \mathbb{C} : \bm{\varepsilon}(\vec{x}) \Big) \; + \mathrm{d}\Omega + = + \int\limits_{\Gamma^h} + N_m(\vec{x}) \; \vec{t}(\vec{x}) \; + \mathrm{d}\Gamma +\end{equation} + +\begin{equation} + \int\limits_{\Omega_0^h} + \tilde{\mathcal{B}}_m^T : \Big( \mathbb{C} : \vec{\nabla}\vec{u}(\vec{x}) \Big) \; + \mathrm{d}\Omega + = + \int\limits_{\Gamma^h} + N_m(\vec{x}) \; \vec{t}(\vec{x}) \; + \mathrm{d}\Gamma +\end{equation} + +\begin{equation} + \int\limits_{\Omega_0^h} + \tilde{\mathcal{B}}_m^T : \Big( \mathbb{C} : ( \tilde{\mathcal{B}}_n \cdot \vec{u}_n) \Big) \; + \mathrm{d}\Omega + = + \int\limits_{\Gamma^h} + N_m(\vec{x}) \; \vec{t}(\vec{x}) \; + \mathrm{d}\Gamma +\end{equation} + +\begin{equation} + \int\limits_{\Omega_0^h} + \tilde{\mathcal{B}}_m^T : \Big( \mathbb{C} : \tilde{\mathcal{B}}_n \Big) \; + \mathrm{d}\Omega + \cdot \vec{u}_n + = + \int\limits_{\Gamma^h} + N_m(\vec{x}) \; \vec{t}(\vec{x}) \; + \mathrm{d}\Gamma +\end{equation} \section{Cylindrical coordinates} \subsection{??} \begin{figure}[htp] \centering \includegraphics[width=.2\textwidth]{figures/cylindrical_coordinates} \caption{Caption here} \label{fig:} \end{figure} \begin{equation} x = r \cos \theta \qquad y = r \sin \theta \qquad z = z \end{equation} \begin{equation} r = \sqrt{ x^2 + y^2 } \qquad \theta = \mathrm{arctan}\left( y/x \right) \qquad z = z \end{equation} \begin{align} \vec{e}_r &= \cos \theta \; \vec{e}_x + \sin \theta \; \vec{e}_y \nonumber \\ \vec{e}_\theta &= - \sin \theta \; \vec{e}_x + \cos \theta \; \vec{e}_y \nonumber \\ \vec{e}_r &= \vec{e}_z \end{align} \begin{equation} \vec{\nabla} = \vec{e}_r \frac{\partial}{\partial r} + \frac{1}{r} \vec{e}_\theta \frac{\partial}{\partial \theta} + \vec{e}_z \frac{\partial}{\partial z} \end{equation} \begin{alignat}{3} \frac{\partial \vec{e}_r}{\partial r } &= 0 &\qquad \frac{\partial \vec{e}_\theta}{\partial r } &= 0 &\qquad \frac{\partial \vec{e}_z}{\partial r } &= 0 \nonumber \\ \frac{\partial \vec{e}_r}{\partial \theta} &= \vec{e}_\theta &\qquad \frac{\partial \vec{e}_\theta}{\partial \theta} &= - \vec{e}_r &\qquad \frac{\partial \vec{e}_z}{\partial \theta} &= 0 \nonumber \\ \frac{\partial \vec{e}_r}{\partial z } &= 0 &\qquad \frac{\partial \vec{e}_\theta}{\partial z } &= 0 &\qquad \frac{\partial \vec{e}_z}{\partial z } &= 0 \end{alignat} We start from the weak form in Eq.~\eqref{eq:static:weak:final}, that has been obtained strictly without specifying a coordinate system. Our point of departure is thus \begin{equation} \int\limits_\Omega \big[\, \vec{\nabla} \vec{\phi}(\vec{x}) \,\big] : \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega = \int\limits_\Gamma \vec{\phi}(\vec{x}) \cdot \vec{t}(\vec{x}) \; \mathrm{d}\Gamma \qquad \forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d \end{equation} The next step is to discretize using Eq.~\eqref{eq:discretization} thereby we have to be a be careful. We begin by making Eq.~\eqref{eq:discretization} more specific: \begin{equation} \phi_r (\vec{x}) \vec{e}_r + \phi_\theta (\vec{x}) \vec{e}_\theta + \phi_z (\vec{x}) \vec{e}_z \approx \sum\limits_{i=1}^{n} \phi_{i,r} N_i (\vec{x}) \vec{e}_r + \phi_{i,\theta} N_i (\vec{x}) \vec{e}_\theta + \phi_{i,z} N_i (\vec{x}) \vec{e}_z \end{equation} We can also evaluate the discretized gradient \begin{align} &\vec{\nabla} \vec{\phi} (\vec{x}) \approx \sum\limits_{i=1}^{n} \bigg( \vec{e}_r \frac{\partial}{\partial r} + \frac{1}{r} \vec{e}_\theta \frac{\partial}{\partial \theta} + \vec{e}_z \frac{\partial}{\partial z} \bigg) \bigg( \phi_{i,r} N_i (\vec{x}) \vec{e}_r + \phi_{i,\theta} N_i (\vec{x}) \vec{e}_\theta + \phi_{i,z} N_i (\vec{x}) \vec{e}_z \bigg) \\ &= \sum\limits_{i=1}^{n} \phi_{i,r } \frac{\partial N_i}{\partial r} \vec{e}_r \vec{e}_r + \phi_{i,\theta} \frac{\partial N_i}{\partial r} \vec{e}_r \vec{e}_\theta + \phi_{i,z } \frac{\partial N_i}{\partial r} \vec{e}_r \vec{e}_z \nonumber \\ &\quad + \phi_{i,r } \frac{1}{r} \frac{\partial N_i }{\partial \theta} \vec{e}_\theta \vec{e}_r + \phi_{i,r } \frac{1}{r} N_i (\vec{x}) \vec{e}_\theta \vec{e}_\theta + \phi_{i,\theta} \frac{1}{r} \frac{\partial N_i }{\partial \theta} \vec{e}_\theta \vec{e}_\theta - \phi_{i,r} \frac{1}{r} N_i (\vec{x}) \vec{e}_r \vec{e}_\theta + \phi_{i,z } \frac{1}{r} \frac{\partial N_i }{\partial \theta} \vec{e}_\theta \vec{e}_z \nonumber \\ &\quad + \phi_{i,r } \frac{\partial N_i }{\partial z} \vec{e}_z \vec{e}_r + \phi_{i,\theta} \frac{\partial N_i }{\partial z} \vec{e}_z \vec{e}_\theta + \phi_{i,z } \frac{\partial N_i }{\partial z} \vec{e}_z \vec{e}_z \end{align} Which we can denote in short as \begin{equation} \vec{\nabla} \vec{\phi} (\vec{x}) \approx \sum\limits_{i=1}^{n} \mathcal{B}_i (\vec{x}) \cdot \vec{\phi}_i = - \underline{\mathcal{B}}^\mathsf{T} \cdot \underline{\vec{\phi}} + \underline{\mathcal{B}}^\mathsf{T} \cdot \vec{\phi}_m \end{equation} where the non-zero components of the third order $\mathcal{B}_i$ read \begin{alignat}{3} \mathcal{B}_{i,r r r } &= \frac{\partial N_i}{\partial r} & \qquad \mathcal{B}_{i,\theta r r } &= \frac{1}{r} \frac{\partial N_i}{\partial \theta} & \qquad \mathcal{B}_{i,z r r } &= \frac{\partial N_i}{\partial z} \nonumber \\ \mathcal{B}_{i,r \theta r } &= - \frac{1}{r} N_i (\vec{x}) & \qquad \mathcal{B}_{i,\theta \theta r } &= \frac{1}{r} N_i (\vec{x}) & \qquad \nonumber \\ \mathcal{B}_{i,r \theta \theta} &= \frac{\partial N_i}{\partial r} & \qquad \mathcal{B}_{i,\theta \theta \theta} &= \frac{1}{r} \frac{\partial N_i}{\partial \theta} & \qquad \mathcal{B}_{i,z \theta \theta} &= \frac{\partial N_i}{\partial z} \\ \mathcal{B}_{i,r z z } &= \frac{\partial N_i}{\partial r} & \qquad \mathcal{B}_{i,\theta z z } &= \frac{1}{r} \frac{\partial N_i}{\partial \theta} & \qquad \mathcal{B}_{i,z z z } &= \frac{\partial N_i}{\partial z} \end{alignat} %\begin{align} % \mathcal{B}_{i,r r r } &= \frac{\partial N_i}{\partial r} \\ % %\mathcal{B}_{i,r r \theta} &= 0 \\ % %\mathcal{B}_{i,r r z } &= 0 \\ % \mathcal{B}_{i,r \theta r } &= - \frac{1}{r} N_i (\vec{x}) \\ % \mathcal{B}_{i,r \theta \theta} &= \frac{\partial N_i}{\partial r} \\ % %\mathcal{B}_{i,r \theta z } &= 0 \\ % %\mathcal{B}_{i,r z r } &= 0 \\ % %\mathcal{B}_{i,r z \theta} &= 0 \\ % \mathcal{B}_{i,r z z } &= \frac{\partial N_i}{\partial z} \\ % \mathcal{B}_{i,\theta r r } &= \frac{1}{r} \frac{\partial N_i}{\partial \theta} \\ % %\mathcal{B}_{i,\theta r \theta} &= 0 \\ % %\mathcal{B}_{i,\theta r z } &= 0 \\ % \mathcal{B}_{i,\theta \theta r } &= \frac{1}{r} N_i (\vec{x}) \\ % \mathcal{B}_{i,\theta \theta \theta} &= \frac{1}{r} \frac{\partial N_i}{\partial \theta} \\ % %\mathcal{B}_{i,\theta \theta z } &= 0 \\ % %\mathcal{B}_{i,\theta z r } &= 0 \\ % %\mathcal{B}_{i,\theta z \theta} &= 0 \\ % \mathcal{B}_{i,\theta z z } &= \frac{1}{r} \frac{\partial N_i}{\partial \theta} \\ % \mathcal{B}_{i,z r r } &= \frac{\partial N_i}{\partial z} \\ % %\mathcal{B}_{i,z r \theta} &= 0 \\ % %\mathcal{B}_{i,z r z } &= 0 \\ % %\mathcal{B}_{i,z \theta r } &= 0 \\ % \mathcal{B}_{i,z \theta \theta} &= \frac{\partial N_i}{\partial z} \\ % %\mathcal{B}_{i,z \theta z } &= 0 \\ % %\mathcal{B}_{i,z z r } &= 0 \\ % %\mathcal{B}_{i,z z \theta} &= 0 \\ % \mathcal{B}_{i,z z z } &= \frac{\partial N_i}{\partial z} \\ %\end{align} -Note that similarly -\begin{equation} - \vec{\nabla} \vec{u} (\vec{x}) - \approx - \sum\limits_{i=1}^{n} - \mathcal{B}_i (\vec{x}) \cdot \vec{u}_i - = - \underline{\mathcal{B}}^\mathsf{T} \cdot \underline{\vec{u}} -\end{equation} -\begin{equation} - \int\limits_\Omega - \left( \underline{\mathcal{B}}^\mathsf{T} \cdot \underline{\vec{\phi}} \right) : \bm{\sigma}(\vec{x}) \; - \mathrm{d}\Omega - = - \int\limits_\Gamma - \underline{N}^\mathsf{T}(\vec{x}) \; \underline{\vec{\phi}} \cdot \vec{t}(\vec{x}) \; - \mathrm{d}\Gamma - \qquad - \forall \; \underline{\vec{\phi}} \in \mathbb{R}^d -\end{equation} - - -\begin{equation} - \underline{\vec{\phi}}^\mathsf{T} \cdot - \int\limits_\Omega - \underline{\mathcal{B}} : \bm{\sigma}(\vec{x}) \; - \mathrm{d}\Omega - = - \underline{\vec{\phi}}^\mathsf{T} \cdot - \int\limits_\Gamma - \underline{N}(\vec{x}) \vec{t}(\vec{x}) \; - \mathrm{d}\Gamma - \qquad - \forall \; \underline{\vec{\phi}} \in \mathbb{R}^d -\end{equation} - - -\begin{equation} - \int\limits_\Omega - \underline{\mathcal{B}} : \bm{\sigma}(\vec{x}) \; - \mathrm{d}\Omega - = - \int\limits_\Gamma - \underline{N}(\vec{x}) \vec{t}(\vec{x}) \; - \mathrm{d}\Gamma -\end{equation} - -\begin{equation} - \int\limits_\Omega - \underline{\mathcal{B}} : \left( \mathbb{C} : \bm{\varepsilon}(\vec{x}) \right) \; - \mathrm{d}\Omega - = - \int\limits_\Gamma - \underline{N}(\vec{x}) \vec{t}(\vec{x}) \; - \mathrm{d}\Gamma -\end{equation} - -\begin{equation} - \int\limits_\Omega - \underline{\mathcal{B}} : \left( \mathbb{C} : \vec{\nabla}\vec{u}(\vec{x}) \right) \; - \mathrm{d}\Omega - = - \int\limits_\Gamma - \underline{N}(\vec{x}) \vec{t}(\vec{x}) \; - \mathrm{d}\Gamma -\end{equation} - -\begin{equation} - \int\limits_\Omega - \underline{\mathcal{B}} : \left( \mathbb{C} :( \underline{\mathcal{B}}^\mathsf{T} \cdot \underline{\vec{u}}) \right) \; - \mathrm{d}\Omega - = - \int\limits_\Gamma - \underline{N}(\vec{x}) \vec{t}(\vec{x}) \; - \mathrm{d}\Gamma -\end{equation} - -\begin{equation} - \int\limits_\Omega - \underline{\mathcal{B}} : \left( \mathbb{C} : \underline{\mathcal{B}}^\mathsf{T} \right) \; - \mathrm{d}\Omega - \cdot \underline{\vec{u}} - = - \int\limits_\Gamma - \underline{N}(\vec{x}) \vec{t}(\vec{x}) \; - \mathrm{d}\Gamma -\end{equation} - \begin{equation} \vec{\nabla}_\xi = \bm{J} \cdot \vec{\nabla} \end{equation} \begin{equation} \vec{e}_\xi \frac{\partial}{\partial \xi } + \vec{e}_\eta \frac{\partial}{\partial \eta} + \vec{e}_\mu \frac{\partial}{\partial \mu} = \bm{J} \cdot \left( \vec{e}_r \frac{\partial}{\partial r } + \frac{1}{r} \vec{e}_\theta \frac{\partial}{\partial \theta} + \vec{e}_z \frac{\partial}{\partial z } \right) \end{equation} \begin{align} \bm{J} &= \vec{e}_\xi \frac{\partial r }{\partial \xi } \vec{e}_r + \vec{e}_\xi r \frac{\partial \theta}{\partial \xi } \vec{e}_\theta + \vec{e}_\xi \frac{\partial z }{\partial \xi } \vec{e}_z \nonumber \\ &+ \vec{e}_\eta \frac{\partial r }{\partial \eta} \vec{e}_r + \vec{e}_\eta r \frac{\partial \theta}{\partial \eta} \vec{e}_\theta + \vec{e}_\eta \frac{\partial z }{\partial \eta} \vec{e}_z \nonumber \\ &+ \vec{e}_\mu \frac{\partial r }{\partial \mu } \vec{e}_r + \vec{e}_\mu r \frac{\partial \theta}{\partial \mu } \vec{e}_\theta + \vec{e}_\mu \frac{\partial z }{\partial \mu } \vec{e}_z \end{align} \subsection{Axisymmetric} % \bibliography{example_refs} \end{document}