ground_vegetation_radi results

[ericvmueller]

I've been revisiting the ground_vegetation_radi verification case to try to understand inaccuracies we might have with radiation attenuation by particles/vegetation. We give the expected solution for attenuation of radiation from vegetation between a hot ceiling to cold floor as:

where kappa is the absorption coefficient (15.46 1/m) and delta is the thickness of the vegetation layer (5cm)

The first thing worth noting is that the guide/expected solution assumes a hot wall temperature of 993 K, but the actual case is run at 923 K, so both the boundary fuel and the particle method are even further off than the figure suggests (the expected value would be ~8.7 kW/m2).

The second thing worth noting is that we state in the guide that the particle prediction is off because we only use a single layer of particles (we actually use two). I wanted to assess the idea that the prediction should be improved with increasing resolution. I ran the particle case with a series of resolutions between 5 cm and 0.3125 cm. I also increased the hot wall temperature to 1273 K. The typical assumption is that we should do okay with ~10 or more grid cells over the attenuation length scale ...which is ~6.5 cm here.

As we can see, even at high resolutions we are pretty far off. So I have two questions:

1. Are we convinced the analytical solution provided fits the scenario? (there isn't a derivation or reference in the guide)
2. Can we do any better than this?

The vegetation_absorb case tests a similar concept and we appear to do okay but these have lower attenuation coefficients.

[ericvmueller]

So I explored increasing the resolution a bit more. For the boundary model I set a fixed cell size (0.00125m) to get 40 layers. For the particle model I have gone all the way down to 0.5 mm (resolving the vegetation layer with 100 particle layers).

The boundary model appears to be converging to the solution in the guide if we add enough resolution. This is not the case for the particles.

I cracked open Siegel and Howell, and to me the case resembles this problem:

where the net heat flux (in the positive x direction) at x=0 is:

In this case we only want the downward flux at x=0 so we set T1 = 0, T2 = 1273 K, and T in the participating medium (vegetation) is 293 K, and we reverse the sign of q. Someone can check my math but this gives:

If the upper surface does not contribute (T2 = 0) we just have emission from an optically thin medium, and the solution matches that of 9.2/9.3 in the verification guide ...and on the other hand if kappa=0 we just see the heat flux from the hot wall ...so I think I'm on the right track.

This means the expected value is ~45 kW/m2. This is closer to the particle solution at high resolution though we are still off by 10%

Can anyone point me to the source for the equation in the FDS guide?

[ericvmueller]

Well... I found my 10% error for the particles at high resolution. The default emissivity is 0.9 which I wasn't taking into account in my analytical solution. Doh.

I'm testing a scheme to locally refine the solution when we have poor resolution with respect to kappa so hopefully we will be much less grid sensitive (which might also help move us in the direction of favoring particles over the boundary model in all scenarios)

[shostikk]

This is very good test case. I agree with the analytical solution. This is basically the same as radiation_plane_layer tests in radiation folder. I did, however, make some mods to the input file to

• reduce time to avoid heating of particles
• measure only incident heat flux instead of gauge
  File is enclosed, so you can compare and check we are doing the same.

The only way (that I can think of) to improve the RTE solution in cases with very few cells (like two in this case), is to abandon the first-order upwind. Instead, I can try to implement first-order exponential approximation of cell-face intensities (exp(-kx) \approx 1-kx, now exp(-kx) \approx 1 there).

It is easier to test this with KAPPA0 than with particles. So, I will start this by adding the pure absorption versions of the plane_layer test cases, and extending them towards the very poor spatial resolution.
ground_vegetation_radi.fds.txt

[ericvmueller]

Sounds good - I made similar modifications/simplifications to my test case.

I have tried a few different ideas, some of which I think are similar to what you are describing (basically storing an estimation at the cell center and the downwind faces, so that the downwind values can be used for the next cell). I can get a much better solution along the shortest path of attenuation (z-direction here) but something is still off as the angle is more oblique. Here's an example of one of my modifications showing the directional intensity at the lower gauge for all different angles of inclination (ie. path lengths through the attenuating medium)

Anyway maybe we can take the discussion offline and combine efforts.

[ericvmueller]

Also, something that came up in this case when I was looking at devices to measure directional intensity is this bit of code was picking up the wrong radiation angle when I was trying to measure at ORIENTATION = 0.707,0.0,0.707

fds/Source/radi.f90

Lines 3410 to 3426 in 85b5026

	ALLOCATE(COSINE_ARRAY(1:NRA))
	ALLOCATE(NEAREST_RADIATION_ANGLE(N_ORIENTATION_VECTOR))
	ALLOCATE(VIEW_ANGLE_AREA(N_ORIENTATION_VECTOR))
	VIEW_ANGLE_AREA = 0._EB
	DO IO=1,N_ORIENTATION_VECTOR
	DO N=1,NRA
	COSINE_ARRAY(N) = ORIENTATION_VECTOR(1,IO)*DLX(N) + &
	ORIENTATION_VECTOR(2,IO)*DLY(N) + &
	ORIENTATION_VECTOR(3,IO)*DLZ(N)
	IF (-(ORIENTATION_VECTOR(1,IO)*DLANG(1,N) + &
	ORIENTATION_VECTOR(2,IO)*DLANG(2,N) + &
	ORIENTATION_VECTOR(3,IO)*DLANG(3,N)) > ORIENTATION_VIEW_ANGLE(IO)) &
	VIEW_ANGLE_AREA(IO) = VIEW_ANGLE_AREA(IO) - COSINE_ARRAY(N)
	ENDDO
	NEAREST_RADIATION_ANGLE(IO) = MINLOC(COSINE_ARRAY,DIM=1)
	VIEW_ANGLE_AREA(IO) = PI/VIEW_ANGLE_AREA(IO)
	ENDDO

I think it should be something more like:

         COSINE_ARRAY(N) = ORIENTATION_VECTOR(1,IO)*DLANG(1,N) + &
                           ORIENTATION_VECTOR(2,IO)*DLANG(2,N) + &
                           ORIENTATION_VECTOR(3,IO)*DLANG(3,N)
         IF (-COSINE_ARRAY(N) > ORIENTATION_VIEW_ANGLE(IO)) &
            VIEW_ANGLE_AREA(IO) = VIEW_ANGLE_AREA(IO) - (ORIENTATION_VECTOR(1,IO)*DLX(N) + &
                                  ORIENTATION_VECTOR(2,IO)*DLY(N) + &
                                  ORIENTATION_VECTOR(3,IO)*DLZ(N))  

What do you think?

[mcgratta]

I'll take a look at the orientation/angle thing.

[mcgratta]

I removed an unused loop that was a leftover from the past when we used to have multiple orientations allowed for a solid particle. This loop was doing nothing and does not address your issue. I just wanted to get it out of the way.

[mcgratta]

DLANG is an array that was introduced this year. I'm not sure exactly what it does and how it's different than DLX etc. I have to look into it more. I'm adding Jason to the issue because I think he added it.

[ericvmueller]

Its just the unit vector in the direction of the radiation angle and DLX etc are integrated over the solid angle. How this is causing the difference when comparing to the device orientation vector I haven't fully understood yet.

[mcgratta]

I cleaned up some of the variable names related to radiation and oriented particles. Should not change functionality.