Determination of the domain of nonlinearity in a Neo-Hook solid model (Finite elements)

Question

For a FEM simulation of a Neo-Hook solid model, how do we know we are in the "regime" of nonlinearity of the solid? In other words, how do I know the hyperelastic material law is really used, and not just a needless additional computational cost because it can be modeled with a linear elastic Hook law. Following these Wikipedia figures for example, does this mean the nonlinearity really happens when we see a nonlinear relationship between $\lambda_{11}$ and $P_{11}$?

For my case/simulation, I chose the element with the biggest stress values, and I calculated the Lagrange Green strain tensor $\pmb{E}$ there, I computed the eigendecomposition to get the biggest eigenvalue $E_{11}$ and its associated eigenvector $\pmb{v}_{11}$. I then compute the stretch as : $\lambda_{11} = \sqrt{1+2E_{11}}$. For the principal stress I compute : $P_{11} = \pmb{v}_{11} \cdot (\pmb{P}\pmb{v}_{11})$ where $\pmb{P}$ is the first Piola-Kirchoff stress. Should I expect to see a nonlinear relationship here? Is my definition of these quantities correct? In my case, I get something linear:

I used $W(\pmb{F}) = \frac{\mu_s}{2}\, (|\pmb{F}|^2 - 3 - 2\, \log(J)) $ with $\mu = 1923076$.

Answer 1

Since your $\lambda_{11}$ is close to unity, it means that your response will be close to linear, because your strain is small. You have not provided enough strain for the non-linearity to kick-in.

I would plot over a larger range of $\lambda_{11}$, say from 1 to 5. That should be enough to see the non-linearity.

Also, before doing FEM, I would derive a non-linear relationship between $\lambda_{11}$ and $P_{11}$, solve the resulting equations with some non-linear solver and plot the results. And then compare with FEM.

Edit: Also try plotting an exponential constitutive relation such as the Veronda Westmann model.

Answer 2

In the FEM simulation, the part is loaded using forces, so I don't control the range of $\lambda_{11}$. Specifically, it's a 2D beam under a bending. So it gives something like this :

As suggested by @NNN I derived the relationship $P_{11}(\lambda_{11})$ and I have $P_{11} =(\lambda_{11}^2 - 1) \mu / \lambda_{11} $, which is indeed not very "nonlinear", unless $\mu$ is very low. In the comparison between FEA (the above plot) and the analytical expression, it does seem to match. The problem is that I probably "sampled" the stress and its principal component the wrong way, since using the principal direction only gives stretches $\lambda \in [.999, 1.0003]$ whether the $yy$ component for example (as shown in the paraview contour) reaches values much higher : If I plot the $P_{yy} = f(\sqrt{1+2*E_{yy}})$ for the two corner elements, I get :

Could you please indicate what's the right way to compute the principal directions for the principal stretch ?

Thank you.