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

Question

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?

Answer 1

It is a matter of boundary conditions on the longitudinal faces. As you noted, the axial stress and strain for a linear isotropic material will satisfy the following relation:

\begin{equation} \sigma_{xx} = \lambda(e_x + e_y + e_z) + 2\mu e_x. \end{equation}

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

\begin{equation} \sigma_{xx} = \frac{\mu(3\lambda+2\mu)}{\lambda+\mu}e_x = Ee_x \end{equation}