[serles]

Hello everybody!

I tried to implement a simple solver for Helmholtz equation (wave equation in frequency domain) of the form:

laplace(p) + k^2 * p = 0

where p is the acoustic pressure and k the complex wave number. Since OpenFOAM doesn't support complex numbers I decomposed the equations in two (introducing p = p_Re + i*p_Im and same for k) and implemented them like this:

solve
                (
                        fvm::laplacian(p_Re)
                        + fvm::Sp(k_Re_sq, p_Re)
                        == k_Im_sq * p_Im
                );

            solve
                (
                        fvm::laplacian(p_Im)
                        + fvm::Sp(k_Re_sq, p_Im)
                        == - k_Im_sq * p_Re
                );

where Im and Re stand for imaginary and real.

As one can see the problem has to be solved iteratively: first for Re(p), then for the Im(p) using the result from the previous equation and afterwards repeat from the beginning until both solutions converge.

Set-up: I tested the code on 1D benchmark with the following conditions:

• wave number is constant for the hole domain
• BC of the acoustic source: Dirichlet with Re(p)=const and Im(p)=const
• BC of the opposite end: Dirichlet with Re(p)=Im(p)=0
• other boundaries are defined as empty (since it's 1D)

Result: The solution doesn't converge for several test cases!
It converges for small frequencies and small wave numbers, but not for higher values. The pressure starts to grow with each iteration step. If it converges, I get good results which perfectly match analytical solution or results from other software.

Possible solutions: I've tried all the numerical tweaks I could find/use in OF:

• different linear solvers (btw GAMG gives an error in almost all of the test cases)
• different interpolation/laplacian/snGrade schemes
• resolution of the domain from 100 cells to 10^6
• relaxation of the equations with different relaxation factors

I've even tried the block matrix solver from the openfoam-extend which is suitable for such coupled problems - no improvement.
Next thing I want to try is the Finite Element solver from foam-extend. I suppose this will bring problems when I start to couple acoustics with fluid dynamics...

But I truly cannot believe that Finite Volume Method is not capable to solve it. I'm not a specialist in CFD, so may be I don't see some simple mistake
Does somebody have any suggestions what else I can try?

I'm open for any help or discussion.

[mkraposhin]

Hi, maybe you must try to present p as vector instead of 2 scalars?

[serles]

Hi, thank you for the input!

But I don't see a way how to introduce such a vector in given equations.
Given the Helmholtz eqn. introduce complex p and k:

Δp + k^2 ⋅ p = 0
p = p_Re + i ⋅ p_Im
k^2 = K_Re + i ⋅ K_Im

Gives two coupled eqns:

Δp_Re + K_Re ⋅ p_Re = K_Im ⋅ p_Im
Δp_Im + K_Re ⋅ p_Im = - K_Im ⋅ p_Re

where the second eqn. is purely imaginary.
If I introduce a vector p_V = (p_Re; p_Im) what am I going to do with the RHS?

If I understood right, something similar, what you mean, the block matrix from the extend project is doing. It incorporates p_Re and p_Im as a vector and solves the equations implicitly. But in my case it doesn't help:

Time = 1

BiCGStab:  Solving for blockT, Initial residual = (0.018707 0.0155012), Final residual = (2.11629e-24 5.10652e-24), No Iterations 1
ExecutionTime = 0.75 s  ClockTime = 0 s

BiCGStab:  Solving for blockT, Initial residual = (0.017386 0.0209848), Final residual = (4.9824e-24 5.08883e-24), No Iterations 1
ExecutionTime = 0.77 s  ClockTime = 0 s

It converges perfectly during a single iteration but the residuals are large again when the new solution is introduced to the solver.

[mkraposhin]

Hi,

ok, now i understand what you are doing.

I would suggest three improvements for weakly coupled iterative solution scheme
1) Use relaxation
2) Use hybrid implicit / explicit source term.
3) Update all sources only once before solving equations

for example:

bool converged = false;
while (converged)
{

    fvScalarMatrix ReEqn
    (
         fvm::laplacian(p_Re)
          + fvm::SuSp(k_Re_sq, p_Re)
          == k_Im_sq * p_Im
    );

    fvScalarMatrix ImEqn
    (
         fvm::laplacian(p_Im)
         + fvm::SuSp(k_Re_sq, p_Im)
         == - k_Im_sq * p_Re
    );

    ReEqn.relax();
    ImEqn.relax();

    ReEqn.solve();
    ImEqn.solve();

    //converged = ....
}

Oh, I'm sorry, i read your first message again and found, that you tried relaxation already. But nevertheless, i think that relaxation will stablize convergence process.
Maybe you need to play with relaxation coefficients, from 0.5 down to 0.05

I hope, this will help

[mkraposhin]

Hi,

I found that solution doesn't converge for cases with non-zero imaginary part of wave number,

you can see example and solver in the attachment

maybe you need different B.C. or treatment for the source due to Im(k^2)?

[serles]

Hi,
thank you for the interest on the problem. I've looked at your solver and it is exactly what I was trying to do. As I found out solving the Helmholtz equation is not such a trivial task: Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods. Not all the solvers are suggested and special preconditioners are needed.


Thus, I persued the approch with the Coupled Block Matrix solver and found out that I initialized the block matrix in a wrong way. In the first try I did it like in foam-extend-3.2/applications/solvers/coupled/blockCoupledScalarTransportFoam. They insert the implicit terms in the diagonal parts of the block matrix but not in the upper and lower parts, which I don't understand why.

After some time I found this thread: http://www.cfd-online.com/Forums/ope...-set-pdes.html, where Tron has done exactly what I needed. It looks now like this:

 
// separate 1st eqn in two parts and store them in separate matrices
// 2nd matrix serves only to insert the coupled terms into the 1st one
fvScalarMatrix Ap_ReM
(
    fvm::laplacian(Ap_Re) + fvm::Sp(k_Re_sq, Ap_Re)
);

Ap_ReM.relax();

fvScalarMatrix Ap_ImM2
(
    -fvm::Sp(k_Im_sq, Ap_Im)
);

Ap_ImM2.relax();

// separate 2nd eqn in two parts and store them in separate matrices
fvScalarMatrix Ap_ImM
(
    fvm::laplacian(Ap_Im) + fvm::Sp(k_Re_sq, Ap_Im)
);

Ap_ImM.relax();

fvScalarMatrix Ap_ReM2
(
    fvm::Sp(k_Im_sq, Ap_Re)
);

Ap_ReM2.relax();

// Prepare block system
BlockLduMatrix<vector2> blockM(mesh);
BlockLduMatrix<vector2> blockM2(mesh);

// Grab block diagonal and set it to zero
Field<tensor2>& d = blockM.diag().asSquare();
d = tensor2::zero;
Field<tensor2>& d2 = blockM2.diag().asSquare();
d2 = tensor2::zero;

// Grab linear off-diagonal and set it to zero
Field<tensor2>& u = blockM.upper().asSquare();
Field<tensor2>& l = blockM.lower().asSquare();
u = tensor2::zero;
l = tensor2::zero;
Field<tensor2>& u2 = blockM2.upper().asSquare();
Field<tensor2>& l2 = blockM2.lower().asSquare();
u2 = tensor2::zero;
l2 = tensor2::zero;

vector2Field& blockX = blockT.internalField();
vector2Field& blockX2 = blockT2.internalField();
vector2Field blockB(mesh.nCells(), vector2::zero);
vector2Field blockB2(mesh.nCells(), vector2::zero);

//- Insert equations into block Matrix
blockMatrixTools::insertEquation(0, Ap_ReM, blockM, blockX, blockB);
blockMatrixTools::insertEquation(1, Ap_ImM, blockM, blockX, blockB);
blockMatrixTools::insertEquation(0, Ap_ReM2, blockM2, blockX2, blockB2);
blockMatrixTools::insertEquation(1, Ap_ImM2, blockM2, blockX2, blockB2);

//- Add off-diagonal terms and remove from block source for the diagonal
forAll(d, i)
{
    d[i](0,1) = d2[i](1,1);
    d[i](1,0) = d2[i](0,0);
    blockB[i][0] -= blockB2[i][1];
    blockB[i][1] -= blockB2[i][0];
}

//- Add off-diagonal terms and remove from block source upper and lower
forAll(u, j)
{
    u[j](0,1) = u2[j](1,1);
    u[j](1,0) = u2[j](0,0);
}

forAll(l, j)
{
    l[j](0,1) = l2[j](1,1);
    l[j](1,0) = l2[j](0,0);
}

//- Block coupled solver call for Matrix1 with inserted terms
BlockSolverPerformance<vector2> solverPerf =
    BlockLduSolver<vector2>::New
    (
        blockT.name(),
        blockM,
        mesh.solutionDict().solver(blockT.name())
    )->solve(blockX, blockB);

solverPerf.print();

// Retrieve solution
blockMatrixTools::blockRetrieve(0, Ap_Re.internalField(), blockX);
blockMatrixTools::blockRetrieve(1, Ap_Im.internalField(), blockX);

Ap_Re.correctBoundaryConditions();
Ap_Im.correctBoundaryConditions();

Anyway, my solver with this implementation works nicely for the 1D case now, but I have still some convergence problems with 2D and 3D for the desired wave numbers. There are two linear solvers I tried out: BiCGStab and GMRES. These are my issues:

1. It doesn't matter what under-relaxation coefficients for the matrices (Ap_ReM, Ap_ImM, …) I use with both solvers. The residuals don't change at all.
2. BiCGStab converges fast but it's unstable (what an irony). Residuals may go down to 1e-5 and then shoot to the sky.
3. GMRES is much more stable but it converges very slow. It's better with higher value of nDirections but it needs ages to converge.

Does someone have an idea about these issues or do you have a tip to try out? The problem about under-relaxation concerns me mainly.


In the attached file you can find the solver (you need foam-extend-3.2 to compile it) and the mentioned test cases, if you want to try it out and get it to converge (faster) =))

[mkraposhin]

Hi,
i have several colleagues, who are working in the same or close area. I shall try to interest them

[through_horizon]

Hi,
which preconditioner you're using with BiCGStab? You could try using GAMG (in segregated case) or AMG with F-cycle in coupled case.

[serles]

Originally Posted by through_horizon

Hi,
which preconditioner you're using with BiCGStab? You could try using GAMG (in segregated case) or AMG with F-cycle in coupled case.

Hi,
I've tried all the available Preconditioners with BiCGStab, Cholesky looked for as the most sufficient.
I've tried GAMG as well, but it couldn't accomplish any iteration. I also found some works telling that Multigrid-solvers are not suitable for wave equations (e.g. here chapter 2.5), which seems logical to me since the "wave-like" structure get lost on a coarser grid. I might try the Multigrid Solver from foam-extend, but from the considerations above I'm doubting it'll be any better.

[Ali Blues]

Dear Sergey,
I was wondering if you were able to resolve your issue with solving the system. I also need to solve the Helmholtz equation, so just tired out the solver you attached for the 2D case. If I use the GMRES solver it just gets stuck at the first iteration (Time=1). Same goes for BiCGStab.
Thanks
Ali

[serles]

Hi Ali,
may be it's just very slow... =) How large is your mesh or do you use one of my test cases? What value of nDirections do you use with GMRES solver?

[Ali Blues]

Hi Sergey,
Thanks a lot for the reply, and sorry for the late response. Well I was looking at your 2d test case. In can't get simpler than this.

Q1) Thanks for the hint regarding nDirections. Indeed as mentioned in the Chalmers University's report by Klas Jareteg on block Coupled Calculations in OpenFOAM, the AMG solver is sensitive to both nDir and maxIter. I set it to nDir=10. But now I'm facing the issue you mentioned regarding convergence. So did you figure this out?

Q2) I noticed from your examples that you are or were working on sono-chemical reactors which is the topic I'm going to look into. From the solver description, I was wondering if you were using the paper of R. Jamshidi et al, Chemical Engr Jr 189-190 (2012) 364-375? since they also discuss eventually an iterative scheme for calculation the volume fraction and hence the k_m.

Q3) Regarding the boundary conditions, I was wondering why would you set the imaginary pressure to same amplitude as the real while it seems that, at first glance, it should be zero ?

Thanks,
Ali

[serles]

Originally Posted by Ali Blues

So did you figure this out?

Unfortunately not. It's sufficient for my problem at this point, since I work with stationary (in frequency domain) acoustic fields now. What I figured out is that, the problem has to do more with the preconditioner than with the solver. There are some but none has been implemented in OF yet. E.g. Analytic Incomplete LU (AILU) preconditioner for Krylov methods proposed by M. J. Gander 2000 or a Multigrid one proposed by Y. A. Erlangga 2006, whereby the latter looks much more complex. It sounds like we need to wait until somebody implements it in OF or do it ourselves =)

Originally Posted by Ali Blues

I was wondering if you were using the paper of R. Jamshidi

Exactly, I'm working with his models.

Originally Posted by Ali Blues

I was wondering why would you set the imaginary pressure to same amplitude as the real while it seems that, at first glance, it should be zero ?

Why should it be zero? The ratio of imaginary and real parts of the pressure contains only the information about the phase angle (phi = arctan(Ap_Im/Ap_Re)). Since we are not interested in phase, it doesn't matter what value Ap_Im or Ap_Re have. It's important that they give the desired pressure in common (Ap = sqrt(Ap_Im^2+Ap_Re^2)). I should mention that my background has not much to do with acousitcs, so, correct me if I'm wrong...
What I also learned about the BCs meanwhile is that one should use the BC of type "fixedGradient" for sonotrode/transducer/vibrating surface. This one is closer to the physics since it imposes (indirectly) a displacement of the surface, which is what actualy happens in the real world.

[Ali Blues]

Hi Sergey,
Thanks for the reply.

With regards to splitting the pressure amplitude at the source/BC between the real and imaginary part, it is true that the amplitude matters at the end but still slightly unsure about the end result being equivalent to case where the amplitude is entirely imposed on the real part, given the coupled nature of the problem. Perhaps since it is linear problem it wouldn't matter. My background is not acoustics either so ... The easiest way is to consider both and check whether we get the same solution. But the problem now is that we just don't have a stable algorithm to run the cases

Also a question regarding the imaginary and real wave numbers values you've included in the transportProperties dict. I just would like to double check I didn't make any errors in the steps to arrive at these values from a given volume fraction of bubbles (Beta), equations 1-8 in in the paper of Jamshidi. I noticed you always have positive values of the wave numbers for all values of beta, which is not the case for me and can be that both are +,-, or the two are of different signs. Say for f=25 kHz, and p0=1 bar=1e5 Pa, and a0=150 microns, and beta=0.1, I get kmSqr.real = -4.4e7, and kmSqr.imag =-8.3e6, for which km=6.7e3 as in figure 4.a in the paper. Could you please confirm this using your code, if possible.

Also, have you seen this recent ppt from openfoam workshop:
http://www.openfoamworkshop.org/asse...es/OFWP053.pdf

They've seem to include a new pre-conditioner for AMG solver, also they're working on new ones as well, so it seems your wish is coming true .

Best,
Ali

[anothr_acc]

Hi everyone.

1) Is work still going on on this?

2) Did you consider using a complex solver rather than splitting real and complex? I'm sure I read somewhere that a complex solver had been implemented.

I'll look around and post back if I find the solver.

Best regards,

Mark.