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?
