# Discretization used for steady linear elasticity - followup to previous question

This is a follow-up of my previous scicomp question (https://scicomp.stackexchange.com/posts/28863/edit). I figured I'd start a new thread on this as the question is a bit different from my previous question.

I am solving the linear elasticity equations, using Hooke's Law and the strain-displacement relationship.

My discretization is as follows using a finite volume formulation:

$$\nabla \cdot \boldsymbol{\sigma} = \nabla \cdot (\mu \boldsymbol{\nabla u} + \mu(\boldsymbol{\nabla u})^T + \lambda\boldsymbol{I}tr(\boldsymbol{\nabla u})) = 0$$

In integral form, we have:

$$\iint_A [\mu \boldsymbol{\nabla u} + \mu(\boldsymbol{\nabla u})^T + \lambda\boldsymbol{I}tr(\boldsymbol{\nabla u})] \cdot dA = 0$$

I evaluate this surface integral using: $$\iint_A [\mu \boldsymbol{\nabla u} + \mu(\boldsymbol{\nabla u})^T + \lambda\boldsymbol{I}tr(\boldsymbol{\nabla u})] \cdot dA = \sum\limits_{k=1}^N |A_k|n_k\cdot[\mu (\boldsymbol{\nabla u})_k + \mu(\boldsymbol{\nabla u})_k^T + \lambda\boldsymbol{I}tr(\boldsymbol{\nabla u})_k]$$

I set the above equation to $f(\boldsymbol{u})$ for simplicity. I can expand this term using a Taylor series expansion of the form: $$f(\boldsymbol{u}) = f(0) + \frac{df}{d\boldsymbol{u}}\boldsymbol{u} = f(0) + \boldsymbol{Ju}$$

I compute the jacobian, $\boldsymbol{J}$, using a forward differencing scheme:

$$J = \frac{d\boldsymbol{f}}{d\boldsymbol{u}} = \frac{f(\boldsymbol{u}+h)-f(\boldsymbol{u})}{h}$$

So this gives us a linear system of the form: $$\boldsymbol{Ju} = \boldsymbol{b} = -f(0)$$

Since this is a finite volume formulation, my mesh consists of cells. I compute the cell gradients using the discrete green gauss method:

$$\boldsymbol{\nabla} u = -\frac{1}{V}\sum\limits_{j=1}^N \hat{u}_jA_j\boldsymbol{n_j}$$

Where $\hat{u}=\frac{u_L+u_R}{2}$, i.e., the face displacement is obtained by averaging the neighboring cells' displacement.

I compute face gradients using a simple average and applying a correction term: $$\nabla u_{face} = \frac{1}{2}(\nabla u_L + \nabla u_R) - \underbrace{\hat{\boldsymbol{d_{lr}}} [\frac{1}{2}(\nabla u_L + \nabla u_R)\cdot\hat{\boldsymbol{d_{lr}}} ] + (u_L - u_R)\frac{\hat{\boldsymbol{d_{lr}}}}{|\boldsymbol{d_{lr}}|}}_{\text{Correction term}}$$

So this is the gist of my formulation.

I tried to model a cantilever problem, and it ended in a disaster (see previous question). I have tested out this formulation on the classical 1-D compression problem and was able to match analytical results.

Does anyone see anything wrong with this formulation?

• Why have you linearized a system that indeed is linear? Do Lame coeff. depend on $\vec{u}$?Wouldn't it be easier with a finite element formulation? Finite element was developed for solving elasticity problems. In fact this problem leads to find the kernel $Ku=0$ of a symmetric positive definite operator $K$. Anyways... When you compute the volume average gradient of $\vec{u}$, if $\vec{n}$ is the outer normal, the minus sign is wrong (see Divergence Th.). If you are using a cartesian mesh (even more convenient if you are modeling a cantilever), there is no need to correct the gradient. – HBR Feb 23 '18 at 7:11
• @HBR I don't have a good reason for linearizing the system other than that I was being lazy. I could form the Jacobian matrix manually, as opposed to the linearizing&using finite difference. Maybe I will try that instead, but I doubt it'll change my solution much? Lame coefficients aren't dependent on $\vec{u}$ in my code. I am developing this code to run along with another FVM code, so I would prefer to use FVM. My mistake on the minus sign (I have it correct in my code). I am currently testing my code on a cartesian mesh, but when I remove the correction term, my solution is very different. – David Feb 25 '18 at 0:13
• The correction term should be 0 for a cartesian mesh, if you use a centred scheme. Is it? – HBR Feb 25 '18 at 10:45
• @HBR Hmmmm, I am definitely not getting 0 for my mesh. I'm currently running this on 4 columns of cells. The cells are all cubes and each face is aligned with one cartesian axis. – David Feb 25 '18 at 20:50
• Why don't you remove the correction part and compute the gradient with the simple average?. Try this. With this, the numerics for this type of equation will be correct. Let's see the results. – HBR Feb 26 '18 at 7:50