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}


\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?


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?

  • 1
    $\begingroup$ At first look this question seems to be off-topic, but it asks about the difference between an analytical solution and a numerical one. $\endgroup$
    – nicoguaro
    Apr 11, 2018 at 4:09
  • $\begingroup$ How many element did you use in your finite element model? Does this error changes when you refine your discretization? $\endgroup$
    – nicoguaro
    Apr 11, 2018 at 4:10
  • $\begingroup$ 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. $\endgroup$
    – user27430
    Apr 11, 2018 at 4:16
  • 2
    $\begingroup$ Would the stress be defined by $\sigma = E \epsilon$ instead of your expression? $\endgroup$
    – nicoguaro
    Apr 11, 2018 at 4:27
  • $\begingroup$ 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}$ $\endgroup$
    – user27430
    Apr 11, 2018 at 4:32

1 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}


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.