How can I implement second order derivatives of shape functions of a 3D elements?

Question

I am developing an Abaqus UEL with 3D 8 nodes brick elements and I need second order derivatives of the shape functions, I have already mapped the first order derivatives from the element coordinates to the physical coordinates through: $\frac{\partial N}{\partial \textbf{X} } = J^{-1} \frac{\partial N}{\partial \textbf{x}} $ where $N$ are the shape functions, $X$ the physical coordinates, $x$ the local coordinates and the Jacobian matrix is $J_{ij}=\frac{\partial N}{\partial x_i} X_j$.

How can I compute and implement the second order derivatives of the shape functions with respect to global coordinates?

Answer 1

You want to compute the following (let me use $\nabla$ instead of your $\frac{\partial}{\partial\textbf{X}}$ for derivatives in real space and $\hat\nabla$ instead of $\frac{\partial}{\partial\textbf{x}}$ for derivatives on the reference cell) $$ (\nabla^2 N)_{ij} = \nabla_i(\nabla_j N) = [J^{-1}]_{ik} \hat\nabla_k ([J^{-1}]_{jl} \hat\nabla_l N). $$ By the product rule, this is $$ (\nabla^2 N)_{ij} = [J^{-1}]_{ik} [J^{-1}]_{jl} (\hat\nabla^2 N)_{kl} + [J^{-1}]_{ik} \hat\nabla_k[J^{-1}]_{jl} \hat\nabla_l \hat N. $$ You know $J^{-1}$ here from computing the first derivatives, and computing the derivatives $\hat\nabla^2 N$ is easy because it happens on the reference cell. The difficulty is the derivative of the inverse of the Jacobian: $$ \hat\nabla [J^{-1}]. $$ This has a simple solution: We know $$ \hat\nabla [JJ^{-1}]=0 $$ and so $$ J \hat\nabla [J^{-1}] + [\hat\nabla J] [J^{-1}] = 0. $$ which you can solve for $$ \hat\nabla [J^{-1}] = -[J^{-1}] [\hat\nabla J] [J^{-1}] $$ where $\hat\nabla J$ is simply the second derivative of the transformation, and is as easily computed as $J$ was in the first place.

Answer 2

While the answer of Wolfgang Bangerth is mathematically sound, it should be noted that doing this for shape functions of an 8 node brick would result in a 0 tensor. The shape functions are all linear, and hence, the second order derivatives are all 0. You would probably need to use higher order elements like a 27 node brick (better than the 20 node serendipity elements IMO) to have second derivatives of the shape functions.