# 1D analytical solution vs FEM solution for a bar under compression

I simulated the compression problem in ANSYS and compared to the analytical solution and found some discrepancies.

The classical solution to the 1-D compression problem is:

\begin{align} u(x) = Cx \end{align}

subjecting to the following BCs:

\begin{align} u(x=0) = 0 \\ \frac{du}{dx}(x=L) = C \end{align}

This solution is derived for 1-D in the absence of any shear phenomena. For a 10 m long rod subject to a compressive load of $\sigma$ at $x=L$, we can use the stress relationship to determine $\frac{du}{dx}$. In 1-D, we can write the following:

\begin{align} \sigma = (2\mu+\lambda)\frac{du}{dx} \end{align}

Thus,

\begin{align} \frac{du}{dx} = \frac{\sigma}{2\mu+\lambda} =C \end{align}

For the purpose of comparing this analytical solution to ANSYS, I will give some numbers:

\begin{align} \mu = 7.6923 \times 10^{10} \text{ Pa}\\ \lambda = 1.153 \times10^{11} \text{ Pa}\\ \sigma = 10^6 \text{ Pa} \\ L = 10 \text{ m} \end{align}

The analytical solution gives us

\begin{align} u(x=L) = 3.7143\times 10^{-5} \text{ m} \end{align}

When I simulate an analogous problem in ANSYS of a 10m rod with a 1m by 1m cross section subjected to the same load, I obtain a maximum displacement of $u(x=L) = 5 \times 10^{-5}\text{ m}$. The full 3-D equations are solved in ANSYS, so we have shear effects, displacements in other 2 directions, and intercomponent coupling among the displacement components. Are these the reasons why the solution is so different?

EDIT

In response to KNL's comment: So I understand that using Hooke's law, and the stress-strain relationship for small displacements, we get the following in 3-D: \begin{align} \boldsymbol{\sigma}= \begin{bmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{zx} & \sigma_{zy} & \sigma_{zz} \\ \end{bmatrix} = \mu \begin{bmatrix} \frac{du}{dx} & \frac{du}{dy} & \frac{du}{dz} \\ \frac{dv}{dx} & \frac{dv}{dy} & \frac{dv}{dz} \\ \frac{dw}{dx} & \frac{dw}{dy} & \frac{dw}{dz} \\ \end{bmatrix} + \mu \begin{bmatrix} \frac{du}{dx} & \frac{du}{dy} & \frac{du}{dz} \\ \frac{dv}{dx} & \frac{dv}{dy} & \frac{dv}{dz} \\ \frac{dw}{dx} & \frac{dw}{dy} & \frac{dw}{dz} \\ \end{bmatrix}^T + \lambda (\frac{du}{dx} + \frac{dv}{dy} + \frac{dw}{dz}) \boldsymbol{I} \end{align}

In 1-D, and absent of shear terms, we get: \begin{align} \sigma_{xx} = (2\mu+\lambda)\frac{du}{dx} \end{align} Which is the relationship that I originally used, but I don't understand what is wrong?

• At first look this question seems to be off-topic, but it asks about the difference between an analytical solution and a numerical one. – nicoguaro Apr 11 '18 at 4:09
• How many element did you use in your finite element model? Does this error changes when you refine your discretization? – nicoguaro Apr 11 '18 at 4:10
• The error does not change when I refine the discretization. I started with a coarse mesh of only about 16 cells and expanded to 4000 cells. There was almost no difference in the error at the loading end. I also tried to the FEM model with a slimmer rod of 0.1m x 0.1m cross section and this made no difference either. – user27430 Apr 11 '18 at 4:16
• Would the stress be defined by $\sigma = E \epsilon$ instead of your expression? – nicoguaro Apr 11 '18 at 4:27
• Hmmm, that's exactly it. But using simplifying the 3-D Hooke's law, I thought $\sigma = (2\mu+\lambda)\frac{du}{dx}$, whereas $\sigma = E\epsilon$ gives $\sigma = E\frac{du}{dx}$ – user27430 Apr 11 '18 at 4:32

$$\sigma_{xx} = \lambda(e_x + e_y + e_z) + 2\mu e_x.$$
If you assume zero transverse strains (for example in the case of sliding contact boundary conditions), you will obtain the result that you did, $\sigma_{xx}=(2\mu+\lambda)e_x$.
However, in the more commonly encountered case of free boundaries, the Poisson effect comes into play. Just substitute $e_y = e_z = -\nu e_x = -\frac{\lambda}{2(\lambda+\mu)}e_x$ in the above relation, to obtain
$$\sigma_{xx} = \frac{\mu(3\lambda+2\mu)}{\lambda+\mu}e_x = Ee_x$$