[zfaraday]

Dear foamers,

Lately I've been doing some tests in a case I'm working in. I'm using chtMultiRegionSimpleFoam and fvDOM for the radiation. I was trying to find the proper amount of rays to be used in my case. After some tests I found out that the Qr at the external patches of my domain behaved in a strange way depending on the number of angular divisions. At the begining I didn't find the patern of the variation of the radiative heat fluxes and it seemed that it was a random fact. Besides that, in that case it wasn't easy to predict the Qr behavior on the external patches because the maximum effect of radiation takes place in the center of the domain. However, in the end I started to suspect that it was related to the number of nTheta angles used in the "constant/fluid_Region/radiationProperties" dictionary.

So as to have a better view and understanding of what it was going on, I took (again) a simple test case I recently used in order to found out that thermalBaffle1D and externalWallHeatFluxTemperature BC's had some bugs in their codes. The case is attached and the explanation of its characteristics is discussed here. Besides what is said in the link, I added an option to generate a 3D version of the case in order to check the nTheta bad behavior. 3D version may be activated in the Allrun script. To be able to run the case swak4foam must be installed, otherwise you must coment out the functions subdictionary in the controlDict file. I used these functions in order to check the energy balance within the domain, it computes the integral of the Qr field over the whole areas.

In order to see the diferences I also attached an Excel file with the results shown by these functions in all the cases I simulated. The .xls file contains some sheets with the test name of the corresponding sheet. The name has a structure such as "test_xy" where x is the number of nPhi and y the number of nTheta defined in the radiationProperties file. In the last sheet a summary table is also added with the variation in percentage of the radiative power in each patch for every test case.

In the last table it is easy to see the patern I mentioned previously: an even number of nTheta angles seems to give physical results according to the conditions of the case, whereas an odd number for nTheta gives totally meaningless values for the radiative power.

Has anyone expericed it before? Why does this diference exist between using an odd or an even number of nTheta angles? What is wrong with the fvDOM for the radiation calculation? Can I really assume that using an even number of nTheta I am getting reasonable results?

Many thanks in advance!

Best regards,

Alex

Note: The .xls file was created in LibreOffice as a .ods and I converted to .xls because of the uploader restrictions. If anyone has problems with the .xls file, let me know and I will upload the original .ods file compressed as .zip.

[zfaraday]

Has anyone experienced this behaviour when using fvDOM for radiation calculation? Does anyone have any idea regarding this issue?

[Astrodan]

I can't help you with your problem, but since I'm working with the fvDOM Model too, I'm at least interested.

If the angular resolution, especially the number if nTheta angles matters, you could try an even simpler test case (compare this thread, modified to 3D). Maybe if I have time I will try that. However, I'm only using a custom solver radiationFoam (you can find on vesion in the linked post), and so far only check if the computed radiation intensitites at the surfaces are correct.

And I came here from your other thread (link) and just want to note that I by now read a whole lot of the fvDOM source code and they don't require the field G to exist anywhere.

[zfaraday]

Thanks Timm for your words.

As I point out in the last lines of my other thread, the fact of defining G in fvDOM doesn't make a diference in the results. However, it's confusing that almost all the tutorials where radiation is computed by means of the fvDOM have G defined in the 0/ directory...

The test case I proposed is also quite simple, but I could also try the even simpler one you proposed if I have time...

Best regards,

Alex

[tniemi]

Hi

I commented this issue in mantis http://www.openfoam.org/mantisbt/view.php?id=1740.

Did my comments make any sense to you? (I probably had an error in my comment and should have talked about fixedWalls instead of floor and ceil, but anyway the problems are caused by rays that are parallel to some walls.)

[zfaraday]

Hi Timo,

I was going to say something in the report but I couldn't after it was closed.

I appreciate your comments and yes, what you say makes sense to me. What I didn't understand much is when you say

One benefit of having (at least some of the) rays align with axes is that then the discretization error is probably smaller with structured grids.

I think it doesn't make much sense since it is exactly the case in the test case I used and, in my opinion, the results are wrong when nTheta are odd.

We can keep talking about this matter in this thread in order to put things in its right place since I want to be sure that everything is going well when I use fvDOM! And I'm still not an expert in this matter...

Thanks for the insight!

[tniemi]

I think it doesn't make much sense since it is exactly the case in the test case I used and, in my opinion, the results are wrong when nTheta are odd.

Yes, you are right, it doesn't make sense in this case, because it is entirely wall driven. However, with absorption/emission and different geometry odd nTheta may work ok and I was just thinking a reason why somebody would want to have it odd. (by discretization error I meant the error that happens inside the domain, when the rays are "convected" from one cell to another. This error is small if the rays are normal to cell faces) But I guess it is a bit dangerous and using even number is probably better in practice at least when small number of rays is used. (if you have time, you could test with a large, odd nTheta. The results should improve)

I want to be sure that everything is going well when I use fvDOM! And I'm still not an expert in this matter...

I'm not an expert either and I don't have huge amounts of experience with fvDOM. However, I have used it in some cases and the results have been reasonable.

[Astrodan]

Hi Alex,
since I already tested the 2D case I mentioned above, it was quite straightforward for me to run the 3D simulation. So here are my general results:

Theory:
The geometry is a cube with a length of 1m. If one surface radiates energy (left) AT T=100K, and all others only absorb (T=0K), the view factors can be calculated giving a view factors of approx. 0.2 for every other surface (compare: equation / calculator).
Thus we have:

• Emitted radiation:
• View factors
• Incident radiation on every surface:

A case for this is attached (cubeFvDOM.zip), written for OpenFOAM 2.3.x for use with the radiationFoam solver (compare this thread). It includes an angleTest bash script, which runs the calculations I have attached. However, note that these are 400 simulations, so I strongly advise to parallise it in some way (e.g. fork the for-loop, setup different cases with for loops from 1-5, 6-10, ..).



Performance:
performance.jpg

Weird enough it seems that the number of phi-angles as nearly no influence on the number of iteratiosn until convergence. However (not surprising) the execution time increases both with nPhi and nTheta. Note that the actual time might vary quite a lot depending on your system. Also I was working on my computer while the simulation ran, and though I tried to run the simulation on cores that were nott busy otherwise this might attribute to variations.
For some simulations my script faild to store the execution time (didn't bother figuring out why), so I replayed the values with NaN (empty patches in the figure).

Results:
results.jpg

The solution seems to converge to the physical result quite fast, with nPhi = nTheta = 5 maybe sufficient. It appears, that the phi-convergence is achieved faster, although (looking in the source code) nPhi produces twice the number oif rays as nTheta does, so no surpise here either. Still, the variations in the solution due to the number of theta-rays seems to be more significant, and oscillations around the true solution are quite apparent, which might be related to your problem. Anyways, for this case I wouldn't consider the solution "unphysical".

The computed results are quite close to the ones expected (see above), though I have not estimated the error that occured during to the CFD calculations. Still, when run with many rays one can still see the effects of the "ray-resolution" on the right wall of the domain (checker-board like).

Evaluation:
Also attached is my MATLAB Script (cubeRadiationTest.zip) including the raw data generated. You can use this to produce above figures or use your own data. Colours are merely to be able to see the 3D surface. I tried to have the colorbar match the paraview colorbar, which sadly did not work (might be due to linux/matlab combination, I had some problems before).


I hope this helps and is an interesting base case for anyone who wants to use this solver. Maybe it even helps with an estimation of how many rays (nPhi, nTheta) are reasonable.

[zfaraday]

Hi Timm,

Thanks for your detailed analysis about the number of rays in fvDOM. It is a matter of big interest for me at this time and it resolves my doubts. In the graphs you posted its quite obvious that the solution depends on the number of nTheta used but not (at least it's not such a strong dependency) on nPhi, something that I stated before in my first post. Thanks again for developing such a big analysis using such an amount of rays!

According to the graphs, it seems that the fact that nTheta is even or odd affects the value of the solution, being more important the difference when the values for both are low and barely having importance when values are large, say more than 10 for each parameter. Besides that, You said that a good value would be nPhi=nTheta=5 and it leads me to a question: comparing to the analytical result which set of values gives less error on the solution: nPhi=nTheta=4 or nPhi=nTheta=5?

By the way, I just realized that the solution seems not to depend on nPhi value at all, which means that the value of the solution would be the same for nPhi=3, nTheta=5 as for nPhi=5, nTheta=5. Having noticed that, Why don't you suggest, for instance, nPhi=3, nTheta=5 instead of nPhi=nTheta=5? The results must be very close one to another but one consumes more than the other.

Again, thanks for your great contribution!

Best regards,

Alex

[Astrodan]

Hi Alex,

there are some good questions, which would require some more analysis of the data. Unfortunately I don't have time right now. Maybe you can check some of the things yourself with the raw data attachted in my previous post somewhere.

To answer at least one question, nPhi=3 pr nPhi=5: I did a 2D simulation of this case before, too, and plottet the relative error of the solution obtained (compared to analytical solutions) over the number of phi-rays. There you could still see some difference up to nPhi=8. At nPhi = 5 the error was somewhere bewteen 1-2%, which sounds acceptable, while nPhi = 3 gave 8-10% error.
Maybe it wold be a good idea to evaluate the data here for the relative error, too.

So, while it looks fine here, I think a different plot would be better to analyse it. Furthermore, you should always consider that the additional computational power between nPhi =3 and nPhi=5 is not that much and should be manageable with current PC systems. Also we have a very simple test case here. For a slightly more complex case too few rays might be quite bad, sinmply because some geometric angles can not be reproduced if the ray resolution is sufficiently fine.