0
$\begingroup$

I am using the Galerkin method (Discontinuous to be precise) to discretize in space the scalar linear wave equation and the explicit second order centered finite difference scheme to discretize in time, resulting in a semi-discrete system of the kind:

$M\ddot{U} + KU=0, \quad U(0)=U_0, \quad \dot{U}(0) = U_1,$

where $M, K$ are the mass and stiffness matrix of the Galerkin method. Now the Idea to study the stability is to observe that if the bilinear form $a$ from which $K$ originates has an $L^2$-orthonormal basis of eigen-functions on the discrete space $V_h$: $a(w_i, v)=\lambda_i (w_i, v), \quad \forall v \in V_h, \quad i = 1 \dots N, \quad \dim{V_h} =N,$

then after using the time integration scheme and expressing $u^k= \sum_{j=1}^N u_j^kw_j$:

$(u^{n+1}-2u^n+u^{n-1}, w_i) + \Delta t^2a(u^n,w_i)=0$

becomes, using the fact that $(w_j, w_i) = \delta_{i,j}:$

$u^{n+1}_i = (2-\lambda_i \Delta t ^2)u_i^n - u_i^{n-1}, \quad \forall i = 1 \dots N$

Now the point is basically solving this second order difference equation and asking that it doesn't oscillate or explode and this results in the CFL condition of the kind $\Delta t \le C h$ (after bounding the $\lambda_i$ with the mesh size $h$).

So being a second order difference equation there can be cases where the discriminat is positive, negative (two complex roots) or zero. Is this the right track to do these kind of stability analysis? Are these requirements all that is needed to have the CFL condition? Thanks a lot for clarifying.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

Yes. That's all there is to the stability condition.

Taking the material properties - shear modulus ($\mu$), bulk modulus ($\kappa$) and density ($\rho$) - into account, the global critical time step is evaluated as the minimum of the critical time step for each element ($\Delta t^e$)

$\Delta t^e = CFL * h^e / c_{\kappa}$

where CFL is the Courant-Friedrichs-Lewy number, $h^e$ is the characteristic length of the element, and $c_{\kappa}$ is the bulk wave speed which is computed as

$c_{\kappa}=\sqrt{\frac{\kappa+4\mu/3}{\rho}}$.

For additional details and numerical examples, I refer you to my papers on elastodynamics and wave propagation. Paper 1 and Paper 2.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks a lot for your answer, I'll definitely read it, but in this case I don't see how the equation $u^{n+1}-(2-\lambda \Delta t^2) u^n + u^{n-1}=0$ has both roots with absolute value less than one (to avoid explosion or oscillation). I've never found anywhere the explicit solution of the difference equation (except in the parabolic case), and I am a bit lost. I know the solution should be $A\rho_1^n + B \rho_2^n$ if two roots, $(A+Bn)\rho^n$ if only one root. Rewrite the difference equation as $u^{n+1}-2\beta u^n + u^{n-1}=0$ the roots are $\rho_{1,2}= \beta \pm \sqrt{\beta^2-1}$. $\endgroup$
    – lucmobz
    Jan 17, 2021 at 7:55
  • $\begingroup$ Now if $|\beta| \ge 1$, in the real roots case there are no solutions with absolute value of the roots less than one, if $\beta < 1$ the roots are complex and become $\rho_{1,2} = \beta \pm i \sqrt{1-\beta^2}$, taking the absolute value, $|\rho_{1,2}| = 1$. I am probably missing something here, or doing some stupid mistake. $\endgroup$
    – lucmobz
    Jan 17, 2021 at 8:02
  • $\begingroup$ You can follow the eigenvalue analysis or von Neumann stability analysis. Have a look at these notes for computing the stability condition using von Neumann stability analysis. math.ubc.ca/~peirce/M257_316_2012_Lecture_8.pdf The stability condition is $\Delta t < 2/\sqrt{\lambda}$. $\endgroup$
    – Chenna K
    Jan 17, 2021 at 13:05
  • $\begingroup$ If you want to do derive the stability condition using eigenvalues, then you have to use the condition $-1 \leq \rho_{1,2} \leq 1$. This, after some algebraic operators, yields the stability condition, $\Delta t < 2/\sqrt{\lambda}$. $\endgroup$
    – Chenna K
    Jan 17, 2021 at 13:19
  • $\begingroup$ Ah, ok 1 is included! $\Delta t < 2 / \sqrt{\lambda}$ means two complex roots $\rho_{1,2} = \beta \pm i\sqrt{1-\beta^2}$ and their absolute value is always 1 under that condition, I tought that one was excluded. Thanks a lot for clarifying and for your references! $\endgroup$
    – lucmobz
    Jan 17, 2021 at 15:19

Your Answer

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

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