First upload

.gitignore

@@ -0,0 +1,2 @@
+*.asv
+*robosan*

drying_example.m

@@ -0,0 +1,73 @@
+close all; clearvars; clc
+
+% Input data
+p_amb = 101325; % ambient pressure in Pascal
+T_amb = 22;     % temperature in Celsius of ambient air
+phi_amb = 0.6;  % relative humidity (0.05 - 1.00) of ambient air
+T_prh = 72;     % temperature in Celsius of preheated air
+
+% Water vapor pressure in Pascal
+vaporPressure = @(T) 0.61121 * exp((18.678 - T/234.5) .* T./(257.14 + T)) * 1000;
+pv_amb = vaporPressure(T_amb); % water vapor pressure of ambient air
+pv_prh = vaporPressure(T_prh); % water vapor pressure of preheated air
+
+% Air humidity in kg water per kg dry air
+Y = 18.01/28.96 * phi_amb * pv_amb ./ (p_amb - phi_amb * pv_amb);
+
+% Relative humidity of preheated air
+phi_prh = p_amb / pv_prh * Y / (18.01/28.96 + Y);
+
+% Air and water vapor constants
+Cpg = 1000; % air specific heat at constant pressure in J / (kg K)
+Cpv = 1860; % water vapor specific heat at constant pressure in J / (kg K)
+Cpl = 4200; % liquid water specific heat at constant pressure in J / (kg K)
+delta_hv_0 = 2500900; % water vaporization enthalpy at 0°C in J/kg
+
+% Enthalpy of preheated air
+h = Cpg * T_prh + Y * (delta_hv_0 + Cpv*T_prh);
+
+% Saturation coordinates
+Y_star_fun = @(T) 18.01/28.96 * vaporPressure(T) ./ (p_amb - vaporPressure(T));
+h_star_fun = @(T) Cpg * T + Y_star_fun(T) * (delta_hv_0 + Cpv*T);
+adb_line_fun = @(T) (h_star_fun(T) - h) / (Y_star_fun(T) - Y) - Cpl*T;
+
+T_star = fzero(adb_line_fun,T_prh);
+Y_star = Y_star_fun(T_star);
+h_star = h_star_fun(T_star);
+
+% Drying rate
+alpha = 20; % heat transfer coefficient in W / (m2 K)
+delta_hv = 2430300; % vaporization enthalpy at T* in J/kg
+Mdot = alpha * (T_prh - T_star) / delta_hv * 3600;  % drying rate in kg / (m2 h)
+
+%% Plot of the Mollier diagram
+temp = [-20, 80];
+mollier([phi_amb, phi_prh, 1.0], [40, 80, 120]*1000, temp, p_amb)
+
+% Add annotations
+hold on
+ax = plot([Y, Y]*1000, [temp(1), T_amb], 'k--');
+label(ax,'Y', 'location','left')
+
+ax = plot([Y_star, Y_star]*1000, [temp(1), T_star], 'k--');
+label(ax,'Y^*', 'location','left')
+
+ax = plot([0, Y]*1000, [T_amb, T_amb], 'k');
+label(ax,'T_{amb}', 'location','top')
+
+ax = plot([0, Y]*1000, [T_prh, T_prh], 'k');
+label(ax,'T_{preheat}', 'location','top')
+
+ax = plot([Y, Y_star]*1000, [T_prh, T_star], 'k', 'LineWidth',0.1);
+label(ax,'adiabatic saturation line', 'location','left','slope')
+
+ax = plot([0, Y_star]*1000, [T_star, T_star], 'k');
+label(ax,'T^*', 'location','top')
+
+quiver(0,T_amb,Y*1000,0,'off','LineWidth',2)
+quiver(Y*1000,T_amb,0,T_prh-T_amb,'off','LineWidth',2)
+quiver(Y*1000,T_prh,(Y_star-Y)*1000,T_star-T_prh,'off','LineWidth',2);
+quiver(Y_star*1000,T_star,-Y_star*1000,0,'off','LineWidth',2)
+
+hold off
+xlabel('Humidity, g H_2O / kg dry air')

mollier.m

@@ -0,0 +1,108 @@
+function mollier(varargin)
+% MOLLIER() without input data draws the Mollier chart for several curves.
+% By default, the fluids are water and air.
+%
+% REQUIRES THE 'LABEL' PACKAGE BY CHAD GREENE
+% https://www.mathworks.com/matlabcentral/fileexchange/47421-label
+%
+% MOLLIER([PHI],[ENT]) draws the Mollier chart with specified curves for
+% the relative humidity and enthalpy. PHI is a numeric array of relative
+% humidities (accept values between 0.01 and 1.00) and ENT is a numeric
+% array of specific entalpy in J/kg. The last PHI value will always be 1.0,
+% independtly from the user input.
+%
+% MOLLIER([PHI],[ENT],[TEM]) draws the Mollier chart as the previous case
+% with the prescribed range of temperature TEM, defined as [t_initial,
+% t_final] in Celsius degrees.
+%
+% MOLLIER([PHI],[ENT],[TEM],P) also prescribes the ambient pressure with
+% the scalar P in Pascal.
+%
+% MOLLIER([PHI],[ENT],[TEM],P,CPG,CPV,DELTA_HV_0) adds custom values for
+% the constant pressure spefici heat for the gas CPG and the vapor CPV in
+% J/(kg K), as well as the liquid vaporization enthalpy at 0°C in J/kg.
+
+% Check if the 'label' community package has been installed
+addons = matlab.addons.installedAddons;
+if ~any(strcmp("label", addons.Name))
+    error(['LABEL add-on not found. Please retrieve it from: ',...
+        'https://www.mathworks.com/matlabcentral/fileexchange/47421-label'])
+end
+
+narginchk(0,7)
+
+% Demo mode
+relHum = [0.05, 0.2, 0.6, 1.0]; % relative humidity array
+enthalpy = [40000, 80000, 120000]; % specific enthalpy array, J/kg
+temp = [-20,80]; % temperature array, Celsius
+pamb = 101325; % ambient pressure, Pa
+Cpg = 1000; % air specific heat at constant pressure, J / (kg K)
+Cpv = 1860; % water vapor specific heat at constant pressure, J / (kg K)
+delta_hv_0 = 2500900; % water vaporization enthalpy at 0°C, J/kg
+
+if nargin > 0
+    relHum = varargin{1};
+    enthalpy = varargin{2};
+
+    % Check limits for relative humidity
+    if min(relHum) < 0 || max(relHum) > 1
+        error('Relative humidity values are out of bounds [0, 1]')
+    end
+
+    % The last value of relative humidity must be 1.0
+    if relHum(end) < 1
+        relHum = [relHum, 1.0];
+    end
+end
+
+if nargin > 2
+    temp = varargin{3};
+end
+
+if nargin > 3
+    pamb = varargin{4};
+end
+
+if nargin > 4
+    Cpg = varargin{5};
+    Cpv = varargin{6};
+    delta_hv_0 = varargin{7};
+end
+
+% Water vapor pressure in Pascal
+T = temp(1):temp(2);
+pv_star = (0.61121 * exp((18.678 - T/234.5) .* T./(257.14 + T))) * 1000;
+
+% Check if vapor pressure gets larger than ambient pressure
+if any(pv_star > pamb)
+    warning(['Some values of the vapor pressure are larger than the ',...
+        'ambient pressure. Please check the temperature interval as ',...
+        'well as the ambient pressure, together with the gas properties.'])
+end
+
+% Air humidity in kg water per kg dry air
+figure
+xlim([0,40])
+hold on
+for phi = relHum
+    Y_phi = 18.01/28.96 * phi * pv_star ./ (pamb - phi * pv_star);
+
+    ax = plot(Y_phi*1000,T,'k');
+    label(ax,['\phi = ', num2str(phi,'%.2f')],'location','right','slope')
+end
+
+for h = enthalpy
+    Y_h = (h - Cpg*T) ./ (delta_hv_0 + Cpv*T);
+
+    % Clear data points beyond 100% relative humidity
+    Y_h (Y_h > Y_phi) = NaN;
+
+    ax = plot(Y_h*1000,T,'k--');
+    label(ax,['h = ', num2str(h/1000), ' kJ/kg'], 'location','center','slope')
+end
+
+hold off
+xlabel('Humidity, g vapor / kg dry gas')
+ylabel('Temperture, Celsius')
+
+end

mollier_script.m

@@ -0,0 +1,40 @@
+close all; clearvars; clc
+
+% Ambient pressure in Pascal
+p = 101325;
+
+% Input temperature array in Celsius
+T = -20:80;
+
+% Water vapor pressure in Pascal
+pv_star = (0.61121 * exp((18.678 - T/234.5) .* T./(257.14 + T))) * 1000;
+
+% Air humidity in kg water per kg dry air
+figure
+xlim([0,40])
+hold on
+for phi = [0.05, 0.2, 0.6, 1.0] % THE LAST MUST BE 1.0
+    Y_phi = 18.01/28.96 * phi * pv_star ./ (p - phi * pv_star);
+
+    ax = plot(Y_phi*1000,T,'k');
+    label(ax,['\phi = ', num2str(phi)],'location','right','slope')
+end
+
+% Enthalpy of moist air (inverse function)
+Cpg = 1000; % air specific heat at constant pressure in J / (kg K)
+Cpv = 1860; % water vapor specific heat at constant pressure in J / (kg K)
+delta_hv_0 = 2500900; % water vaporization enthalpy at 0°C in J/kg
+
+for h = [40000, 80000, 120000]
+    Y_h = (h - Cpg*T) ./ (delta_hv_0 + Cpv*T);
+
+    % Clear data points beyond 100% relative humidity
+    Y_h (Y_h > Y_phi) = NaN;
+
+    ax = plot(Y_h*1000,T,'k--');
+    label(ax,['h = ', num2str(h/1000), ' kJ/kg'], 'location','center','slope')
+end
+
+hold off
+xlabel('Y, g H_2O / kg dry air')
+ylabel('T, Celsius')