Improvement of the P1-model

[Jonathan Gerardin]

Hello,
 
I'm starting a improvement of the P1-model. Currently, the radiative calculation doesn't take account of the scattering, and the P1-model is not efficiency in thin optically media. I will make an improvement of this model :
- addition of the scattering, possibility to have a linear anisotropic scattering.
- improvement of the P1-model in optically thin media thanks to the IDA (Improved Differential Approximation).
 
I have made the first improvement and compare the result with those of Tan and Hsu (Transient radiative transfer in three-dimensional homogeneous and non-homogeneous participating media, JQSRT 1992). They calculate the radiative heat flux in a box of 1m x 1m x 1m with a face at 86.28K, the other faces and the media are at 0K, the optical thickness is equal to 1 with an albedo of 0.1 (so absorption = 0.9 and scattering = 0.1). The scattering is isotropic.
I have put the results of the P1-model in 3 cases + 1 result in DOM (128 directions) by ploting the radiative heat flux along the center line of the box, which starts in the center of the "hot" face and finish in the center of the opposite face :
- P1 : absorption = 1 (black)
- P1 : absorption = 0.9 (red)
- P1 : absorption = 0.9 + scattering = 0.1 (blue)
-DOM : absorption = 1 (green)
 
The result in DOM are very close to those of Tan and Hsu. It's a little disturbing because we don't consider the scattering in DOM model.
For the P1 model, we show an improvement of the result when we take account of the scattering.
 
Now, i will start the add of the IDA. The IDA consist to launch ray from all the wall to measure the optically thickness between the walls and the elements. There is a weighting between the result of the P1-model and the radiosity of the wall through the optically thickness. The IDA will improve the results in optically thin media (when the radiosity from the wall is more accurate than the P1-model) and near the wall (in optically thick media, when you are near the wall, the optical thickness could be thin because of the small distance).

[Alexandre Douce]

Great job ! ;)
I understand that your are disappointed, but keep in mind that the DO model is better (and more CPU expensive) than the P-1 model. It is
likely that the improvement of P-1 is not enough to improve sufficiently
to make up the overall quality of the model, comparing to the DO model. I imagine that if you take into account the scattering in the
DO model, the results would be even closer to the curve of Tan and
Hsu.
 
In order to analyze your curves more precisely could you add the results of Tan and
Hsu?
 
 

[David Monfort]

As Alexandre wrote, great job ;-)
Keep us informed of your progress!
David

[Jonathan Gerardin]

Here is a better comparaison of the result. And you will see that the DOM128 is very close to the result of Tan and Hsu even if there is no diffusion in the DOM calculation.

[Jean Deschodt]

Hi Jonathan,
For the DO model, 128 directions, it is a very fine angular discretization. Do you try the coarser 32 directions parameter?
Nevertheless, what you did is a useful functionality for the radiative transfer model :D Do you plan to do the add scattering for the DO model too?
 

[Jonathan Gerardin]

Hi,
 
I haven't try the 32 directions parameter, but I think the 32 directions will stay better than the P1 model with scattering.
My work will only be on the P1, I know this method largely better than the DOM. But the modification I made in the P1 model to add the scattering was very tiny, so it's maybe easy to add it in the DOM model too.
 
To implement the IDA, I need the position of each element and each wall (and if it's a wall or a symmetry), how can I find that?
 
Thank you.

[Jean Deschodt]

I agree with you for the 32 directions, but how close from the experimental data should be interesting.

For geometrical data:
1) The center of gravity of cells of the mesh is stored in the array XYZCEN.
2) The center of gravity of intern faces and boundary faces are stored in arrays CDGFAC and CDGFBO.
3) The type of each boundary face is stored in the array ITYPFB (see usclim.f90 or usray2.f90). For a wall you have ITYPFB(IFAC) = IPAROI or IPARUG, for a symmetry you have ITYPFB(IFAC) = ISYMET

You find more informations in the § 5.2 in the user's guide

[Jonathan Gerardin]

Hello,
 
I've got a problem with raypun.f90. I think my source file is not good.
I tried a calculation without raypun.f90 in the SRC folder of my calculation. It works.
I tried the same calculation with raypun.f90 (without any modification inside) in the SRC folder of my calculation. It doesn't work. The calculation stops because of absorption coefficient which is not enough high. I'm ok with that because P1 doesn't work well in optically thin media (before put the raypun.f90 in the SRC folder, I can make a P1 calculation with very low absorption coefficient and there was no warning, so it wasn't normal).
So i put a high absorption coefficient. But it doesn't work too. And I don't understand why. (it's not the same error of absorption coefficient too low, no mention of that in the listing file)
Here is my error file, my listing file and my XML file. (no error of compile)
 
Thank you for your help!

[David Monfort]

Hmmm, I think you have a mismatch between the Code_Saturne version and the raypun.f90 version.

I assume you are using 2.0-rc1 (since this is what is written in the listing file). But the error message you get indicates that the solver used for the linear equation of the radiative transfert is the multigrid algorithm ; and this algorithm was activated by default in the 2.0-beta2 but not in the 2.0-rc1 as it could lead to numerical unstabilities... This is the only difference between both versions of the subroutine.

Can you test again with setting imgr1 = 0 (line 291 in beta2, or 278 in rc1) in raypun.f90 ?

[David Monfort]

However, I don't know what could cause the different behavior for optically-thin media ; but if I understand correctly, the current behavior is better than the previous one ?
Alexandre, could you check that point ?