I've read that in modeling structures problems, the finite element method (FEM) is typically used. I am unfamiliar with FEM, but I am wondering, in particular, if using FEM, as opposed to finite difference/volume, alleviates some of the numerical stability issues that may arise due to the large structural constants?

I am currently trying to solve the transient linear elastic equations (see https://en.wikipedia.org/wiki/Linear_elasticity#Direct_tensor_form) using FVM. For ease of demonstration, we can consider the 1-D form (just the classical second order hyperbolic wave equation) of this equation in the absence of shear and assuming constant density:

$\frac{\partial^2 u}{\partial t^2} = \frac{2\mu + \lambda}{\rho} * \frac{\partial^2 u}{\partial x^2}$.

Where $u$ is the x-displacement, $\rho$ is the mass density and $\mu$ and $\lambda$ are Lame constants (see http://scienceworld.wolfram.com/physics/LameConstants.html). I solved this by turning the second derivative equation into two first order equations: Let $$\frac{\partial^2 u}{\partial t^2} = \frac{dz}{dt}$$ where $$z = \frac{du}{dt}$$ or in vector form:

$$\frac{d\boldsymbol{q}}{t} = \boldsymbol{g}$$ where $\boldsymbol{q} = [u \ \ \ \ \ q]$ and $g = [z \ \ \ \ \ \frac{2\mu+\lambda}{\rho}* \frac{\partial^2u}{\partial x^2}=f $]

In preliminary modeling, I used Explicit Euler and set my constants = 1. The code "appeared" to be stable. Then I set the constants to be representative of steel and some other materials of interest, where the Lame constants fall between O(10^8)-O(10^11). This severely limited my the timestep size that I could use for my code.

Due to this, I decided to try to implement an implicit scheme. To do this, I tried to formulate my problem using the Jacobian matrix as such $$ \frac{d\boldsymbol{q}}{dt} = A\boldsymbol{q} $$ where A is Jacobian matrix $$ A = \frac{d\boldsymbol{g}}{d\boldsymbol{q}} = \begin{bmatrix} 0 & I \\ \frac{df}{du} & 0 \end{bmatrix} $$

Then I found the eigenvalues of this Jacobian matrix and found eigenvalues with a POSITIVE REAL root, indicating that this problem can't be stable (regardless of the large order of magnitude of the structural constants) at least with most conventional time integration schemes.

Earlier, when I stated that my explicit Euler scheme "appeared" to be stabled. I realized that it was mildly unstable. I had used very small structural constants such that the eigenvalues hugged the imaginary axis and the real roots magnitudes are smaller, but still positive.

Going back to my first question: Does FEM alleviate some of these mentioned issues? Second question: Does anyone have insight into solving such a problem where the jacobian matrix's eigenvalues have positive real roots? Third question: Any recommended schemes for tackling linear elasticity with FV that won't be severely limited by stability with EXPLICIT schemes?

  • $\begingroup$ Your one dimensional equation can be rewritten with the constants being of magnitude 1. $\endgroup$ – nicoguaro May 23 '17 at 5:09
  • $\begingroup$ Could you please elaborate on this? Like normalizing/nondinensionalizing? $\endgroup$ – David May 23 '17 at 5:25
  • $\begingroup$ Yes, exactly that is what I was mentioning $\endgroup$ – nicoguaro May 23 '17 at 5:27
  • $\begingroup$ Hmm how would you perform the nondimensionalization? I assumed that's what you meant but I don't understand how this helps with the stability issue because when you nondimensionize the constant by another constant(s), you would just be introducing other constants such that all the constants multiplied together would still be the same magnitude as before? $\endgroup$ – David May 23 '17 at 5:37
  • $\begingroup$ That's true.... But in this case, just make $(2\mu+\lambda)/\rho=1$. That won't make much difference but then you don't have the "big parameters" as part of your problem. $\endgroup$ – nicoguaro May 23 '17 at 5:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.