I want to perform a hygrothermal simulation of a wall on which every day a shadow is cast for some time by a building on the other side of the street. WUFI does not offer an option to allow for such a shadow, but I could simply use a self-created *.KLI file by converting the measured radiation myself and allowing for the sadow in this process. But, how is the conversion of the radiation data done?
First you need to determine the radiation incident on the surface of your building element from the measured data describing the radiation on a horizontal surface. For this purpose it is necessary to determine the position of the sun in the sky at the time of the measurement.
Position of the Sun:
Let J be the number of the day in the year (1 .. 365 or 366). Then compute the auxiliary quantity x:
x = 0.9856° * J - 2.72°
and the equation of time Z (in minutes):
Z = -7.66*sin(x) - 9.87*sin( 2*x + 24.99° + 3.83°*sin(x) ) [min].
The equation of time describes the variable difference in time between the actual culmination of the sun and noon. Because of the ellipticity of the Earth's orbit and the obliquity of the Earth's axis the sun wanders with slightly irregular speed across the sky. During the course of the year there are thus times where it reaches culmination earlier than a fictitious sun with constant speed (the so-called 'mean' sun) and times where it reaches culmination later.
The local meridian is the great circle that rises from the horizon due north, passes through the point directly above the observer and crosses the horizon again due south. The instant at which the sun crosses the local meridian on its daily path from east to west is also the instant where its position is due south and where it reaches its daily greatest height.
When the apparent sun (i.e. the actually observed sun) crosses the meridian it is 12 noon local apparent solar time (LAT); when the mean sun crosses the meridian it is 12 noon local mean time (LMT). The equation of time is therefore the difference between LAT and LMT (Z = LAT - LMT).
Furthermore, since the place where the measurements were taken is usually not located on the reference meridian of the time zone (15° East for the Central European Time Zone, CET), the difference between local mean time and zone time must be allowed for, which is 4 minutes for 1° difference in geographical longitude L and one hour for 15° difference. If the measurement was timed in Central European Summer Time CEST, convert to CET first by subtracting one hour (CET = CEST - 1h).
In this way you can now compute the corresponding local apparent time LAT from the known measurement time (in CET):
LAT = CET - (15°-L)/(15°/h) + Z/(60 min/h) [h]
and thus determine the position of the sun: at 12 noon LAT the sun is exactly on the meridian, before noon it stands at an appropriate distance to the east of the meridian, after noon, an appropriate distance to the west.
The distance between the sun and the meridian is measured by the hour angle omega:
omega = (LAT - 12h) * 15°/h.
The hour angle omega is reckoned perpendicular to the meridian; it is negative before noon, zero at noon and positive after noon; it increases steadily by 15° per hour.
The hour angle gives the distance of the sun from the meridian; the declination delta, i.e. the distance of the sun from the celestial equator, then fixes the position of the sun completely. The declination varies between -23°26' at winter solstice, 0° at the equinoxes, and 23°26' at the summer solstice. Since its change during one day is very small, it suffices to compute it once for the day J under consideration:
sin(delta) = 0.3978 * sin( x - 77.51° + 1.92° * sin(x) ),
cos(delta) = sqrt(1 - sin(delta)^2)
where x is the auxiliary quantity introduced above.
The last step is the transformation from the coordinate system determined by omega and delta into the more familiar coordinates altitude gamma and azimuth psi (=compass direction). The geographical latitude phi of the measurement location is needed for this.
sin(gamma) = cos(delta)*cos(omega)*cos(phi)+sin(delta)*sin(phi)
cos(gamma) = sqrt(1 - sin(gamma)^2)
if cos(gamma)=0 then psi = 0
else begin
sin(psi) = cos(delta)*sin(omega)/cos(gamma)
cos(psi) = (cos(delta)*cos(omega)*sin(phi)-sin(delta)*cos(phi))/cos(gamma)
psi = atn2(sin(psi), cos(psi))
end
This formula uses atn2(A,B), the arctangent function for two arguments A and B, which is provided by many programming languages, and which gives the arctangent of A/B in the correct quadrant. If this function is not available to you, you can use the ordinary arctangent and then explicitly determine the correct quadrant (i.e. you compute psi=atn(A/B), and in the case B<0 you add 180° if psi<0 or subtract 180° if psi>0. If B=0 and A<0 then psi=-90°, if B=0 and A>0 then psi=+90°).
The azimuth psi is counted from south=0°, positive towards the west and negative towards the east.
Examples for Munich (48.13°N, 11.58°E):
CET Altitude Azimuth Declination
(J=1):
1 Jan. 2001 09:00 6.436° -44.614° -22.987°
1 Jan. 2001 12:00 18.836° -4.209° -22.977°
1 Jan. 2001 16:00 3.441° 49.590° -22.962°
1 Jan. 2001 16:25 0.434° 54.316° -22.961°
(J=79):
20 Mar. 2001 07:00 6.498° -82.661° -0.126°
20 Mar. 2001 12:21 41.851° -0.041° -0.038°
20 Mar. 2001 16:00 22.726° 62.242° +0.022°
(J=172):
21 Jun. 2001 08:00 34.501° -87.522° +23.437°
21 Jun. 2001 12:00 65.126° -8.437° +23.437°
21 Jun. 2001 18:00 19.771° 103.475° +23.436°
These values were computed with an astronomical ephemeris program. Of course, the simplified method described above cannot reproduce these data exactly, in particular for low altitudes of the sun (1 Jan. 16:25), since it does not allow for atmospheric refraction. On the other hand, the comparison allows you to assess the overall accuracy of this simple method. Your results should agree with these exact positions within a few tenths of a degree. The declinations have been included as well for testing purposes.
Converting the Radiation Data:
We assume that your input data are measured hourly values of the global (I_glob) and the diffuse radiation (I_diff) on a horizontal surface.
The radiation incident on the measuring or the component surface is split up into a direct and a diffuse component. The direct component is received directly from the sun and is therefore a directed quantity that depends on the position of the sun. The direct radiation vertically incident on a surface which is facing the sun is the direct normal radiation I_dir_normal. The direct radiation I_dir obliquely incident on a horizontal measuring surface depends on the solar altitude gamma:
I_dir = I_dir_normal * sin(gamma).
Since I_dir can be computed as the difference between the measured values of global and diffuse radiation and gamma can be determined from the measurement location and time by the method given above, the corresponding direct normal radiation is
I_dir_normal = (I_glob - I_diff) / sin(gamma).
The angle of incidence eta, i.e. the angle that the direct normal radiation makes with the normal to the component surface which is tilted by the angle beta and oriented in the direction alpha, is
cos(eta) = sin(gamma)*cos(beta) + cos(gamma)*sin(beta)*cos(alpha-psi)
eta: Angle of incidence (vertical=0°)
gamma: Altitude of the sun
psi: Azimuth of the sun (south=0°, positive towards west, negative towards east)
beta: Tilt of the component surface (vertical wall=90°)
alpha: Azimuth of the normal to the component surface (south=0°, west positive).
The direct radiation incident on the component surface is therefore:
I_dir_in = I_dir_normal * cos(eta) = (I_glob - I_diffus) * cos(eta) / sin(gmma).
The diffuse component consists of the radiation scattered by the air ("blue sky") and the clouds which comes from all directions and can approximately be treated as isotropic. Diffuse radiation is measured by blocking the direct radiation with a shadow ring around the solarimeter. The measurement gives I_diff, the diffuse radiation incident on the horizontal measuring surface from the entire sky hemisphere. A component surface with arbitrary tilt and orientation receives the same diffuse radiation (since it is isotropic), but for non-horizontal surfaces the fact has to be allowed for that the sky covers a smaller part of its field of view and the total amount of incident diffuse radiation is reduced proportionately (a vertical wall sees sky only in the upper half of its field of view):
I_diff_in = I_diff * ( cos(beta/2) )^2.
Additionally, you may add the global radiation reflected from the ground:
I_refl_in = rho * I_glob * ( sin(beta/2) )^2,
where rho is the short-wave albedo of the ground and the reflection is assumed to be isotropic. In the current version, WUFI ignores the reflected component of the radiation.
The total radiation incident on the surface of the building component is the sum of the components:
I_in = I_dir_in + I_diff_in + I_refl_in.
You may now modify or supplement this conversion method according to your needs. For example, you can allow for shadows by setting the direct radiation to zero at times where the sun is behind the obstacle, and by reducing at all times the diffuse radiation in proportion to the reduction of the field of view caused by the obstacle. On the other hand, at times where the sun illuminates the facing side of the obstacle, it may be necessary to add some reflected radiation.
Hint: if the radiation data to be converted have been averaged over some longer interval (e.g. one hour), please note the following:
It is advisable to compute the solar positions for the middle of the measuring interval, i.e. the averaged data measured between 9h and 10h should be converted using the solar position computed for 9:30h.
If the sun has risen or set during such a measuring intervall (which is easy to check for, using the solar altitude), the solar position must be computed for the middle of the visibility interval, not for the middle of the measuring interval.
Independent of the duration of the measuring interval, radiation data obtained at very low solar altitudes should not be used, since under these circumstances the direct normal radiation must be calculated from very small and unreliable values obtained for the direct radiation at grazing angles of incidence.
Details on these conversion methods can be found in:
VDI 3789 Umweltmeteorologie, Blatt 2: Wechselwirkungen zwischen Atmosphäre und Oberflächen; Berechnung der kurz- und der langwelligen Strahlung.
In addition to data on global and diffuse radiation, the weather file IBP1991.WET included with WUFI contains radiation data obtained with a west-facing solarimeter which you can use to test your conversion routines.