I am trying to solve the in-homogeneous Helmholtz equation with damping and forcing using finite element method (FEM) with FEniCSx. The equation is;

$c^2\nabla^2p - c^2\omega i d \nabla^2p + \omega^2 p = f$

where the $d$ is the damping parameter, the $i$ is complex number, the f is the forcing, the $\omega$ is the eigenvalue and the $p$ is the eigenvector. In matrix form;

$(\mathbf{A} + \omega \mathbf{B} + \omega^2 \mathbf{C})\mathbf{p} = \mathbf{D}\mathbf{p} $

where $\mathbf{B}$ is the damping matrix and $\mathbf{D}$ is the forcing matrix. For 1D case, the boundary conditions are Neumann boundary at both ends.

I am aiming to damp the system with $\mathbf{B}$ and did try 4 different cases;

  • $\mathbf{A} + \omega^2 \mathbf{C} = 0 $ (no damping and no forcing)
  • $\mathbf{A} + \omega \mathbf{B} + \omega^2 \mathbf{C} = 0$ (damping and no forcing)
  • $(\mathbf{A}-\mathbf{D}) + \omega^2 \mathbf{C} = 0 $ (no damping and forcing)
  • $(\mathbf{A}-\mathbf{D}) + \omega \mathbf{B} + \omega^2 \mathbf{C} = 0 $ (damping and forcing)

When I run this cases, I get very weird behaviour. Here is the list of the eigenvalues for the cases;

  • $\mathbf{A} + \omega^2 \mathbf{C} = 0 $ gives $\omega=1069.2399-0j$
  • $\mathbf{A} + \omega \mathbf{B} + \omega^2 \mathbf{C} = 0$ gives $\omega=1139.6967-0.2463j$
  • $(\mathbf{A}-\mathbf{D}) + \omega^2 \mathbf{C} = 0 $ gives $\omega=1074.748+3.2431j$
  • $(\mathbf{A}-\mathbf{D}) + \omega \mathbf{B} + \omega^2 \mathbf{C} = 0 $ gives $\omega=1146.0758+3.5129j$

The complex part of the eigenvalue is representing the growth rate of the system. As I impose a damping, the last case's complex part (+3.5129j) should be lower than the 3. case's complex part (3.2431j) according to the law of superposition but it is not.

Would anyone can enlight me about what I am doing wrong, I am struggling with this for more than a week.

Thank you!

  • $\begingroup$ Usually, the forced Helmholtz equation is not an eigenvalue problem. Could you elaborate? $\endgroup$
    – nicoguaro
    Commented Nov 30, 2023 at 0:25
  • $\begingroup$ Sure, the forcing is a function of an eigenvalue $\omega$. Hence I find the matrix D iteratively and I use a method called fixed point iteration. The fundamental paper about this problem stated in Section IV.B in hal.science/file/index/docid/908192/filename/paper.pdf $\endgroup$ Commented Nov 30, 2023 at 10:29
  • $\begingroup$ Which eigenvalues are being shown here? The one with the largest magnitude? Largest real part? $\endgroup$
    – whpowell96
    Commented Nov 30, 2023 at 13:26
  • $\begingroup$ These eigenvalues belong to 1st axial mode of the 1m tube. $\endgroup$ Commented Nov 30, 2023 at 14:32
  • 1
    $\begingroup$ Even though you are approximating the solution via a sequence of standard eigenvalue problems, the problem you are attempting to solve is a quadrativ eigenvalue problem and you are trying to udnerstand the effect of perturbing one of the matrices in the problem via a damping matrix. I don't see any easy way of analyzing the effect of this perturbation. $\endgroup$
    – whpowell96
    Commented Nov 30, 2023 at 20:37


Your Answer

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