Skip to content

Commit e52dd76

Browse files
DanRileycwhanse
DanRiley
authored and
cwhanse
committed
Added initial Boyle/Coello soiling model (#26)
* Added initial Boyle/Coello soiling model * Corrected error in syntax defs of hsu_soiling * Changes to Boyle/Coello soiling model per comments * Updated reference for Coello/Boyle soiling model
1 parent 85efebf commit e52dd76

File tree

1 file changed

+355
-0
lines changed

1 file changed

+355
-0
lines changed

‎pvl_soiling_hsu.m

Lines changed: 355 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,355 @@
1+
function SR = pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, varargin)
2+
% PVL_SOILING_HSU Calculates soiling rate over time given particulate and rain data
3+
%
4+
% Syntax
5+
% [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10)
6+
% [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, ModelType)
7+
% [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, ModelType, RainAccPeriod)
8+
% [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, ModelType, RainAccPeriod, LUC, WindSpeed, Temperature)
9+
%
10+
% Description
11+
% [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10)
12+
% Uses the time, rainfall amounts, tilt angle, particulate matter
13+
% concentrations, and rain cleaning threshold to predict soiling rate
14+
% using the fixed settling model. Uses hourly rain accumulation to
15+
% determine cleaning.
16+
% [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, ModelType)
17+
% Uses the time, rainfall amounts, tilt angle, particulate matter
18+
% concentrations, and rain cleaning threshold to predict soiling rate
19+
% using the model type selected by ModelType. Uses hourly rain
20+
% accumulation to determine cleaning.
21+
% [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, ModelType, RainAccPeriod)
22+
% Uses the time, rainfall amounts, tilt angle, particulate matter
23+
% concentrations, and rain cleaning threshold to predict soiling rate
24+
% using the model type selected by ModelType. The user specifies the
25+
% rain accumulation period (in hours) over which to determine whether
26+
% sufficient rain has fallen to clean the module.
27+
% [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, ModelType, RainAccPeriod, LUC, WindSpeed, Temperature)
28+
% Uses the parameters of the fixed models along with Land Use Category
29+
% (LUC), wind speed, and ambient temperature to predict soiling ratio
30+
% using the variable deposition velocity model.
31+
%
32+
% Input Parameters:
33+
% Time is a struct with elements as described in pvl_maketimestruct.
34+
% Time values for the soiling function do not need to be
35+
% regularly spaced, although large gaps in timing are discouraged.
36+
%
37+
% Rain is a vector of rainfall in mm of the same length as Time, where
38+
% each element of Rain correlates with the element of Time.
39+
% Rainfall values should be in mm of rainfall. Programmatically, rain
40+
% is accumulated over a given time period, and cleaning is applied
41+
% immediately after a time period where the cleaning threshold is
42+
% reached.
43+
%
44+
% RainThresh is a scalar for the amount of rain, in mm, in an accumulation
45+
% period needed to clean the PV modules. In periods where the
46+
% accumulated rain meets or exceeds RainThresh, the panels are assumed
47+
% to be cleaned immediately after the accumulation period
48+
% [1] suggests a good RainThresh could be 1mm, but the time period is
49+
% not specified. Rain accumulation period length can be adjusted in the
50+
% optional input RainAccPeriod.
51+
%
52+
% Tilt is a scalar or vector for the tilt of the PV panels. Changing tilt
53+
% angles (e.g. in tracking cases) can be accomodated, and tilt angles
54+
% are correlated with the entries in Time.
55+
%
56+
% PM2_5 is the concentration of airborne particulate matter (PM) with diameter less
57+
% than 2.5 microns, in g/m^3. Historical US concentration data sets may be
58+
% found at https://aws.epa.gov/aqsweb/airdata/download_files.html. With
59+
% file descriptions (for hourly data sets) at
60+
% https://aqs.epa.gov/aqsweb/airdata/FileFormats.html#_hourly_data_files
61+
%
62+
% PM10 is the concentration of airborne particulate matter (PM) with diameter less than
63+
% 10 microns, in g/m^3. Historical US PM data sets found as described
64+
% above.
65+
%
66+
% ModelType is an optional input to the function to determine the
67+
% the model type to be used in the soiling model, see [1]. A
68+
% value of "1" indicates that the Variable Deposition Velocity model
69+
% shall be used, a value of "2" indicates that the Fixed Settling
70+
% Velocity model shall be used, and a value of "3" indicates that the
71+
% Fixed Deposition Velocity model shall be used. [1] indicates that the
72+
% Fixed Settling Velocity model performs best under a wide range of
73+
% conditions, and thus "2" is the default ModelType if ModelType is omitted.
74+
% Validation efforts by Sandia National Laboratories
75+
% confirm these findings. If an incorrect ModelType is provided, the
76+
% Fixed Settling Velocity (type 2) will be used (with a warning).
77+
%
78+
% RainAccPeriod is an optional input that specifies the period, in hours,
79+
% over which to accumulate rainfall totals before checking against the
80+
% rain cleaning threshold. For example, if the rain threshold is
81+
% 0.5 mm per hour, then RainThresh should be 0.5 and RainAccPeriod
82+
% should be 1. If the threshold is 1 mm per hour, then the values
83+
% should be 1 and 1, respectively. The minimum RainAccPeriod is 1
84+
% hour. The default value is 1, indicating hourly rain accumulation.
85+
% Accumulation periods exceeding 24 (daily accumulation) are not
86+
% recommended.
87+
%
88+
% LUC is an optional input to the function, but it is required for the
89+
% Variable Deposition Model. LUC is the Land Use Category as specified
90+
% in Table 19.2 of [2]. LUC must be a numeric scalar with value 1, 4,
91+
% 6, 8, or 10, corresponding to land with evergreen trees, deciduous
92+
% trees, grass, desert, or shrubs with interrupted woodlands. If
93+
% omitted, the default value of 8 (desert) is used.
94+
%
95+
% WindSpeed is an optional input to the function, but is required for the
96+
% Variable Deposition Model. WindSpeed is a scalar or vector value with
97+
% the same number of elements as Time, and must be in meters per
98+
% second. If WindSpeed is omitted, the value of 2 m/s is used as
99+
% default.
100+
%
101+
% Temperature is an optional input to the function, but is required for
102+
% the Variable Deposition Model. Temperature is a scalar or vector
103+
% value with the same number of Elements as Time and must be in degrees
104+
% C. If Temperature is omitted, the value of 12 C is used as default.
105+
%
106+
%
107+
% Output Parameters:
108+
% SR = The soiling ratio (SR) of a tilted PV panel, this is a number
109+
% between 0 and 1. SR is a time series where each element of SR
110+
% correlates with the accumulated soiling and rain cleaning at the times
111+
% specified in Time.
112+
%
113+
% Background and recognition:
114+
% This soiling model was developed by Liza Boyle and Merissa Coello at
115+
% Humboldt State University (HSU). Many thanks to Liza and Merissa for
116+
% model development, initial code, and subsequent code and performance
117+
% reviews.
118+
%
119+
% References
120+
% [1] M. Coello and L. Boyle, "Simple Model For Predicting Time Series
121+
% Soiling of Photovoltaic Panels," in IEEE Journal of Photovoltaics.
122+
% doi: 10.1109/JPHOTOV.2019.2919628
123+
% [2] Atmospheric Chemistry and Physics: From Air Pollution to Climate
124+
% Change. J. Seinfeld and S. Pandis. Wiley and Sons 2001.
125+
%
126+
% See also PVL_MAKETIMESTRUCT
127+
128+
p=inputParser;
129+
130+
p.addRequired('Time', @isstruct);
131+
p.addRequired('Rain', @(x) all(isnumeric(x) & all(x>=0) & isvector(x)));
132+
p.addRequired('RainThresh', @(x) all(isnumeric(x) & isscalar(x) & all(x>=0)));
133+
p.addRequired('Tilt', @(x) all(isnumeric(x) & isvector(x)));
134+
p.addRequired('PM2_5', @(x) all(isnumeric(x) & all(x>=0) & isvector(x)))
135+
p.addRequired('PM10', @(x) all(isnumeric(x) & all(x>=0) & isvector(x)))
136+
p.addOptional('ModelType', 2, @(x) all(isscalar(x)));
137+
p.addOptional('RainAccPeriod', 1, @(x) all(isscalar(x) & x>=1 & isnumeric(x)));
138+
p.addOptional('LUC', 8, @(x) all(isscalar(x)))
139+
p.addOptional('WindSpeed', 2, @(x) all(isnumeric(x) & all(x>=0) & isvector(x)));
140+
p.addOptional('Temperature', 12, @(x) all(isnumeric(x) & all(x>=0) & isvector(x)));
141+
p.parse(Time, Rain, RainThresh, Tilt, PM2_5, PM10, varargin{:});
142+
143+
RainAccPeriod = p.Results.RainAccPeriod;
144+
ModelType = p.Results.ModelType;
145+
LUC = p.Results.LUC;
146+
WindSpeed = p.Results.WindSpeed;
147+
Temperature = p.Results.Temperature;
148+
149+
% Create a datenum from the time structure to handle duration easily
150+
TimeAsDatenum = datenum(Time.year, Time.month, Time.day, Time.hour, ...
151+
Time.minute, Time.second);
152+
153+
% Following section accumulates the rain in each accumulation period
154+
% To change this to a more frequent accumulation rate, create a
155+
% new array based on TimeAsDatenum where the unit is the given unit of
156+
% accumulation (e.g. for hourly accumulation you might use
157+
% TimeAsDatenum*24)
158+
RainAccAsDatenum = floor(TimeAsDatenum * 24 ./ RainAccPeriod);
159+
% Determine the number of unique days, the first instance of that day, and
160+
% create an index of which entries are part of that day
161+
[RainAccTimes, UnqRainAccFrstVal, UnqRainAccIndx] = unique(RainAccAsDatenum);
162+
% Accumulate rain readings according to the unique index of days
163+
RainAtAccTimes = accumarray(UnqRainAccIndx, Rain);
164+
% Create an accumulated rain vector, then populate the last entry of each
165+
% date with the rain that fell on that day
166+
AccumRain = zeros(size(Rain));
167+
AccumRain(UnqRainAccFrstVal(2:end)-1) = RainAtAccTimes(1:end-1);
168+
% Rain accumulated on the last day is placed at the last time (does not
169+
% affect soiling calculation)
170+
AccumRain(end) = RainAtAccTimes(end);
171+
172+
switch ModelType
173+
case 1 % Variable Deposition Velocity
174+
vd = depo_veloc(Temperature, WindSpeed, LUC);
175+
case 2 % Fixed Settling Velocity in m/s
176+
vd(1,2) = 0.004;
177+
vd(1,1) = 0.0009;
178+
case 3 % Fixed Deposition Velcoity in m/s
179+
vd(1,2) = 0.0917;
180+
vd(1,1) = 0.0015;
181+
otherwise
182+
warning('Unknown Model Type specified, using Fixed Setting Velocity model')
183+
vd(1,2) = 0.004;
184+
vd(1,1) = 0.0009;
185+
end
186+
187+
188+
% Corrects PM10 measurement since PM2.5 is included in measurement
189+
PMConcentration=nan(numel(TimeAsDatenum),2); % pre-allocate with NAN
190+
191+
PMConcentration(:,1) = PM2_5; % fill PM2.5 data in column 1
192+
% Since PM2.5 particulate is counted in PM10 measurements, and we need only
193+
% the particular from 2.5 micron to 10 micron, subtract the 2.5 micron
194+
% particulate from the 10 micron particulate
195+
PMConcentration(:,2) = PM10 - PM2_5; % fill in PM2.5-PM10 data in column 2
196+
197+
% For any instances where 2.5 micron particulate was more thatn 10 micron,
198+
% set the amount of 2.5 to 10 micron particulate to 0
199+
PMConcentration(PM10 - PM2_5 < 0 , 2) = 0;
200+
201+
% EPA concentrations reported in micrograms/m^3, we need it in g/m^3
202+
PMConcentration = PMConcentration .* 1e-6;
203+
204+
% Calculate the mass deposition rate of each particulate size
205+
F=PMConcentration .* vd; % g * m^-2 * s^-1, by particulate size
206+
HorizontalTotalMassRate = F(:,1) + F(:,2); % g * m^-2 * s^-1, total
207+
208+
% Mass depositition rate on a tilted surface. In the event of a changing
209+
% tilt (e.g. tracking system) the mass rate is adjusted by the tilt angle
210+
% at each time step, and thus a linear interpolation of changing tilt
211+
% angles are used (when using trapezoidal integration).
212+
TiltedMassRate = HorizontalTotalMassRate .* cosd(Tilt);
213+
214+
% Mass on the tilted surface if no rain occurs. Uses trapezoidal
215+
% integration and determines the amount of mass deposited on a tilted
216+
% module.
217+
TiltedMassNoRain = cumtrapz(TimeAsDatenum*86400, TiltedMassRate);
218+
219+
% Create a Tilted soiling mass that will account for rain cleaning
220+
TiltedMass = TiltedMassNoRain;
221+
% For each rain accumulation period, check to see if enough rain fell to meet or
222+
% exceed the cleaning threshold. If the threshold was met or exceeded,
223+
% then subtract from each subsuequent time step the mass of future time
224+
% steps. This essentially "resets" the cumulative soiling accumulation when
225+
% the rain meets or exceeds the stated thresshold. Note that this means
226+
% that the mass will never be reported as zero because the cleaning is
227+
% applied immediately after a time step, and new mass accumulates between
228+
% the cleaning and the first time reported after the cleaning.
229+
% Note that this could be done more clearly by going through each time step
230+
% and accumulating or clearing based on daily rainfall, but this is much
231+
% faster, and can utilize trapezoidal soiling integration.
232+
for cntr1 = 1:numel(RainAtAccTimes)-1
233+
if RainAtAccTimes(cntr1) >= RainThresh
234+
TiltedMass(UnqRainAccFrstVal(cntr1+1):end) = TiltedMass(UnqRainAccFrstVal(cntr1+1):end)-TiltedMass(UnqRainAccFrstVal(cntr1+1)-1);
235+
end
236+
end
237+
238+
% Soiling rate is the % transmittance reduction using Heagzy equation
239+
SoilingRate = 34.37 * erf(0.17 .* (TiltedMass.^0.8473));
240+
241+
% Soiling Ratio is the % soil transmittance (0.95 = 5% soiling loss)
242+
SR = (100 - SoilingRate) ./ 100;
243+
244+
end
245+
246+
% This function creates deposition velocities for the variable deposition
247+
% model.
248+
function vd = depo_veloc(T, WindSpeed, LUC)
249+
250+
% convert temperature into Kelvin
251+
T = T + 273.15;
252+
253+
% save wind data
254+
u = WindSpeed;
255+
256+
g=9.81; %gravity in m/s^2
257+
Na=6.022*10^23; %avagadros number
258+
R=8.314; %Universal gas consant in m3Pa/Kmol
259+
k=1.38*10^-23; %Boltzmann's constant in m^2kg/sK
260+
P=101300; %pressure in Pa
261+
rhoair= 1.2041; %density of air in kg/m3
262+
z0=1;
263+
rhop=1500; %Assume density of particle in kg/m^3
264+
265+
switch LUC
266+
case {1, 4}
267+
gamma = 0.56;
268+
case {6, 8, 10}
269+
gamma = 0.54;
270+
otherwise
271+
warning('Unknown Land Use Category (LUC), assuming LUC 8.');
272+
gamma = 0.54;
273+
end
274+
275+
276+
%Diameter of particle in um
277+
Dpum=[2.5,10];
278+
Dpm=Dpum*10^-6; %Diameter of particle in m
279+
280+
%Calculations
281+
mu=1.8*10^-5.*(T./298).^0.85; %viscosity of air in kg/m s
282+
nu=mu/rhoair;
283+
lambda1=2*mu./(P.*(8.*0.0288./(pi.*R.*T)).^(1/2)); %mean free path
284+
285+
Cc = 1+2.*[lambda1 lambda1]./Dpm.*(1.257+0.4.*exp(-1.1.*Dpm./([lambda1 lambda1].*2))); %slip correction coefficient
286+
287+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
288+
%Calculate vs
289+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
290+
vs = rhop.*Dpm.^2 .* (g .* Cc ./(mu.*18)); %particle settling velocity
291+
292+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
293+
%Calculate rb
294+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
295+
ustar = NaN(size(u)); % pre-allocate ustar
296+
%Equation 11.66 in Ramaswami (and 16.67 and Sienfeld &Pandis)
297+
ustar(u > 0) = 0.4 * u(u > 0) ./ log(10/z0);
298+
ustar(u<=0) = 0.001;
299+
300+
D=k*T.*(Cc./(3*pi*mu*Dpm));
301+
302+
Sc=nu./D;
303+
%gamma=.56; %for urban
304+
%alpha=1.5; %for urban
305+
EB=Sc.^(-1 * gamma);
306+
St = vs .* (ustar .^ 2) ./ (g .* nu);
307+
308+
EIM=10.^(-3./St); %For smooth surfaces
309+
%EIM=((St)./(0.82+St)).^2;
310+
311+
R1=exp(-St.^(1/2)); %percentage of particles that stick
312+
313+
rb=1./(3*(EB+EIM).*ustar.*R1);
314+
315+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
316+
%Calculate ra
317+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
318+
a=[-0.096,-0.037,-0.002,0,0.004,0.035];
319+
b=[0.029,0.029,0.018,0,-0.018,-0.036];
320+
321+
% For wind speeds <= 3, use a = -0.037 and b = 0.029
322+
% For wind speeds >3 and <=5, use a = -.002, b = 0.018
323+
% For wind speeds > 5, use a = 0, b = 0
324+
avals = a(2) .* ones(size(u));
325+
avals(u>3) = a(3);
326+
avals(u>5) = a(4);
327+
328+
bvals = b(2) .* ones(size(u));
329+
bvals(u>3) = b(3);
330+
bvals(u>5) = b(4);
331+
332+
L = 1 ./ (avals + bvals .* log(z0));
333+
334+
zeta0=z0./L;
335+
zeta=10./L;
336+
eta = ((1-15.*zeta).^(1./4));
337+
eta0 = ((1-15.*zeta0).^(1./4));
338+
339+
ra = NaN(size(zeta)); % Preallocate memory
340+
ra(zeta == 0) = (1 ./ (0.4 .* ustar(zeta == 0))) .* log(10 ./ z0);
341+
ra(zeta > 0) = (1./(0.4.*ustar(zeta > 0))).*(log(10./z0) + ...
342+
4.7.*(zeta(zeta > 0)-zeta0(zeta > 0)));
343+
ra(zeta < 0) = (1 ./ (0.4 .* ustar(zeta < 0))) .* (log(10 ./ z0) + ...
344+
log((eta0(zeta < 0).^2 + 1) .* (eta0(zeta < 0)+1).^2 ./ ...
345+
((eta(zeta < 0).^2 + 1) .* (eta(zeta < 0)+1).^2)) + ...
346+
2 .* (atan(eta(zeta < 0))-atan(eta0(zeta < 0))));
347+
348+
349+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
350+
%Calculate vd and mass flux
351+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
352+
353+
vd=1./(ra+rb)+vs;
354+
355+
end

0 commit comments

Comments
 (0)