Creating Helmholtz Solver: no convergence
Hello everybody!
I tried to implement a simple solver for Helmholtz equation (wave equation in frequency domain) of the form: Quote:
Code:
solve 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:
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:
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 :confused: Does somebody have any suggestions what else I can try? I'm open for any help or discussion. |
Hi, maybe you must try to present p as vector instead of 2 scalars?
|
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: HTML Code:
Δp + k^2 ⋅ p = 0 HTML Code:
Δp_Re + K_Re ⋅ p_Re = K_Im ⋅ p_Im 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: Code:
Time = 1 |
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: Code:
bool converged = false; Maybe you need to play with relaxation coefficients, from 0.5 down to 0.05 I hope, this will help |
1 Attachment(s)
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)? |
1 Attachment(s)
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: Code:
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) =)) |
Hi,
i have several colleagues, who are working in the same or close area. I shall try to interest them |
Hi,
which preconditioner you're using with BiCGStab? You could try using GAMG (in segregated case) or AMG with F-cycle in coupled case. |
Quote:
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. |
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 |
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? |
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 |
Quote:
Quote:
Quote:
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. |
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 |
Great interest in this project
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. |
All times are GMT -4. The time now is 18:36. |