diff --git a/PharmaPy/Drying_Model.py b/PharmaPy/Drying_Model.py index 47e847e..a55fc13 100644 --- a/PharmaPy/Drying_Model.py +++ b/PharmaPy/Drying_Model.py @@ -222,8 +222,6 @@ def get_y_equilib(self, temp_cond, x_liq, p_gas): p_sat = self.Liquid_1.AntoineEquation(temp=temp_cond) # [Pa] gamma = self.Liquid_1.getActivityCoeff(mole_frac=x_liq_mole_frac) # [-] - # p_gas_total = np.sum(x_liq_mole_frac * p_sat[:, self.idx_volatiles], - # axis=1) p_partial = (gamma * x_liq_mole_frac * p_sat[:, self.idx_volatiles]).T # [Pa] y_equil = p_partial / p_gas # [-] @@ -252,11 +250,9 @@ def get_drying_rate(self, x_liq, temp_cond, y_gas, p_gas): y_gas_mole_frac = self.Vapor_1.frac_to_frac(mass_frac=y_gas) y_equil = self.get_y_equilib(temp_cond, x_liq, p_gas) - y_volat = y_gas_mole_frac[:, self.idx_volatiles].T # * p_gas - dry_volatiles = self.k_y * self.a_V * (y_equil - y_volat).T - - # dry_volatiles = self.k_y * (y_equil - y_volat).T - dry_rates = np.zeros_like(y_gas_mole_frac) + y_volat = y_gas_mole_frac[:, self.idx_volatiles].T # [-] + dry_volatiles = self.k_y * self.a_V * (y_equil - y_volat).T # [mol/m**3/s] + dry_rates = np.zeros_like(y_gas_mole_frac) # [mol/m**3/s] dry_rates[:, self.idx_volatiles] = dry_volatiles dry_rates[dry_rates < 0] = 0 @@ -337,7 +333,7 @@ def unit_model(self, time, states, sw=None): sat_red = np.clip(sat_red, 0, 1) # [-] k_ra = (1 - sat_red)**2 * (1 - sat_red**1.4) # [-] vel_gas = self.k_perm * k_ra * self.dPg_dz / visc_gas # [m/s] - + # ---------- Drying rate term mw_avg_gas = np.dot(y_gas, self.Vapor_1.mw) # legacy MW surrogate [g/mol] rho_gas = self.pres_gas / gas_ct / temp_gas * mw_avg_gas / 1000 # [kg/m**3] @@ -512,19 +508,19 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq, When ``return_terms`` is True, returns the diagnostic terms ``(convec_term, drying, heat_cond, heat_loss_emp)`` instead. ``convec_term`` is the raw ``u_gas * dTg_dz`` diagnostic - [kg*K/m**3/s]. ``drying`` preserves the legacy diagnostic - denominator on this branch; PR #124 owns the condensed latent-heat - diagnostic correction. + [kg*K/m**3/s]. ``drying`` is the condensed-temperature latent + contribution [K/s]; ``heat_cond`` and ``heat_loss_emp`` are + gas-temperature-rate contributions [K/s]. Notes ----- - This branch fixes the gas material-balance holdup denominator. The - legacy latent-heat factor in ``drying_terms`` is left unchanged here; - PR #124 owns that scoped correction. The existing gas-convection - discretization is also preserved; ``dTg_dz`` [kg*K/m**4] and - ``conv_term`` [J*kg/m**6/s] are annotated as implemented so that the - remaining dimensional debt is explicit without expanding this branch's - behavior change. + ``latent_heat`` is requested on a mass basis [J/kg], so multiplying by + ``dry_rate`` gives the full latent power [J/m**3/s]. The existing + gas-convection discretization is preserved in this branch; ``dTg_dz`` + [kg*K/m**4] and ``conv_term`` [J*kg/m**6/s] are annotated as + implemented so that the remaining dimensional debt is explicit without + expanding the #24 behavior change. The legacy gas molecular-weight + surrogate is documented as issue #28 rather than changed here. """ mw_avg_gas = np.dot(y_gas, self.Vapor_1.mw) # legacy MW surrogate [g/mol] @@ -547,10 +543,14 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq, cpl_mix = self.Liquid_1.getCp(temp=temp_sol, mass_frac=xliq_extended, basis='mass') # [J/kg/K] - temp_wb = 22+273 # [K] - sensible_heat = cpg_mix * (temp_gas - temp_wb) * dry_rate.sum(axis=1) # [J/m**3/s] - - heat_transf = self.h_T_j * self.a_V * (temp_gas - temp_sol) # [J/m**3/s] + temp_wb = 22 + 273 # [K] + sensible_heat = ( + cpg_mix * (temp_gas - temp_wb) * dry_rate.sum(axis=1) + ) # [J/m**3/s] + + heat_transf = ( + self.h_T_j * self.a_V * (temp_gas - temp_sol) + ) # [J/m**3/s] heat_loss = np.zeros_like(temp_gas) # [J/m**3/s] fluxes_Tg = high_resolution_fvm(temp_gas, boundary_cond=temp_gas_inputs) # [K] @@ -567,7 +567,7 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq, dens_liq = self.rho_liq # [kg/m**3] heat_loss_cond = np.zeros_like(temp_sol) # [J/m**3/s] drying_terms = ( - dry_rate[:, self.idx_volatiles] * latent_heat * 2 + dry_rate[:, self.idx_volatiles] * latent_heat ).sum(axis=1) # [J/m**3/s] denom_cond = self.rho_sol * (1 - self.porosity) * self.cp_sol + \ self.porosity * satur * cpl_mix * dens_liq # [J/m**3/K] @@ -576,7 +576,7 @@ def energy_balance(self, time, temp_gas, temp_sol, satur, y_gas, x_liq, if return_terms: self.convec_term = u_gas * dTg_dz # [kg*K/m**3/s] - self.drying = drying_terms/ denom_gas # [K/s] + self.drying = drying_terms / denom_cond # [K/s] self.heat_cond = heat_transf/ denom_gas # [K/s] self.heat_loss_emp = heat_loss/ denom_gas # [K/s] @@ -621,6 +621,10 @@ def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True, distributed and uniform single-node cake initial conditions. Cake permeability is computed as ``1 / (alpha * rho_sol * (1 - porosity))`` [m**2]. + Liquid density [kg/m**3] computed from the initial liquid composition + is used for the irreducible-saturation estimate. During integration, + ``unit_model`` recomputes ``self.rho_liq`` from the current liquid + state before material and energy balances are evaluated. """ p_atm=101325 # [Pa] @@ -640,20 +644,14 @@ def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True, states_in_dict = dict(zip(self.names_states_in, len_states_in)) self.states_in_dict = {'Inlet': states_in_dict} - # Molar fractions - # y_gas_init = np.tile(self.Vapor_1.mole_frac, (self.num_nodes,1)) - # x_liq_init = self.CakePhase.Liquid_1.mole_frac[:, idx_volatiles] - y_gas_init = self.Vapor_1.mass_frac # [-] x_liq_init = self.CakePhase.Liquid_1.mass_frac # [-] satur_init = self.CakePhase.saturation # [-] - # Temperatures temp_cond_init = self.CakePhase.Solid_1.temp # [K] temp_gas_init = self.Vapor_1.temp # [K] - z_cake = self.CakePhase.z_external # [m], for drying_script_inyoung - # z_cake = self.CakePhase.z_external # This line for 2MSMPR_Filter.py + z_cake = self.CakePhase.z_external # [m] if x_liq_init.ndim == 1: x_liq_init = x_liq_init[idx_volatiles] @@ -681,9 +679,6 @@ def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True, temp_cond_init = np.ones_like(satur_init) * temp_cond_init if isinstance(temp_gas_init, float): temp_gas_init = np.ones_like(satur_init) * temp_gas_init - # temp_gas_init = np.tile(temp_gas_init, (self.num_nodes,1)) - # temp_cond_init = np.tile(temp_cond_init, (self.num_nodes,1)) - states_stacked = np.column_stack( (satur_init, y_gas_init, x_liq_init, temp_gas_init, temp_cond_init)) @@ -691,23 +686,11 @@ def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True, interp_obj = CubicSpline(z_cake, states_stacked) states_prev = interp_obj(self.z_centers) - # states_init = states_prev - # Merge states and interpolate in the grid nodes - # states_prev = np.column_stack((satur_init, y_gas_init, x_liq_init, - # temp_gas_init, temp_cond_init)) - - # states_init = define_initial_state(state=states_stacked, z_after=self.z_centers, - # z_before=z_cake, indexed_state=True) - - #z_cake = self.cake_height * self.CakePhase.z_external - - # Physical properties alpha = self.CakePhase.alpha # [m/kg] rho_sol = self.Solid_1.getDensity() # [kg/m**3] porosity = self.CakePhase.porosity # [-] xliq = states_prev[:, num_comp + 1: num_comp + 1 + self.num_volatiles] - # xliq = states_init[:, num_comp + 1: num_comp + 1 + self.num_volatiles] xliq_init = np.zeros((self.num_nodes, num_comp)) xliq_init[:, self.idx_volatiles] = xliq @@ -717,20 +700,13 @@ def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True, mass_frac=xliq_init) # [N/m] self.k_perm = 1 / alpha / rho_sol / (1 - porosity) # [m**2] - # self.rho_liq = rho_liq self.rho_sol = rho_sol self.porosity = porosity self.cp_sol = self.Solid_1.getCp() # [J/kg/K] - # Mass transfer moments = self.Solid_1.getMoments(mom_num=[0, 1, 2, 3, 4]) - sauter_diam = moments[1] / moments[0] # [m] - # self.a_V = 6 / sauter_diam # [m**2/m**3] self.a_V = moments[2] * (1 - porosity) / moments[3] # [m**2/m**3] - # Gas pressure - # deltaP_media = deltaP*self.resist_medium / \ - # (alpha*rho_sol*self.cake_height + self.resist_medium) deltaP_media = deltaP*self.resist_medium / \ (alpha*rho_sol*(1 - porosity)*self.cake_height + @@ -743,12 +719,10 @@ def solve_unit(self, deltaP, runtime=None, time_grid=None, any_event=True, self.pres_gas = np.linspace(p_top, p_top - deltaP, num=self.num_nodes) # [Pa] - # Irreducible saturation - x_csd = self.Solid_1.x_distrib #* 1e-6 - csd = self.Solid_1.distrib# * 1e6 - mom_zero = self.Solid_1.moments[0] + x_csd = self.Solid_1.x_distrib # [m] + csd = self.Solid_1.distrib # [-] + mom_zero = self.Solid_1.moments[0] # [-] - # rholiq_mass = np.mean(rho_liq[0] * self.Liquid_1.mw_av[0]) # [kg/m**3] self.s_inf = get_sat_inf(x_csd, csd, deltaP, porosity, self.cake_height, mom_zero, (np.mean(surf_tens), rho_liq[0])) # [-] diff --git a/tests/test_drying_energy_rate_basis.py b/tests/test_drying_energy_rate_basis.py index 14e5e91..c66b1eb 100644 --- a/tests/test_drying_energy_rate_basis.py +++ b/tests/test_drying_energy_rate_basis.py @@ -117,12 +117,11 @@ def test_unit_model_uses_mass_drying_rate_in_energy_terms(monkeypatch): 2475.0 / 700.0, ]) # [K/s] - # The mass-rate latent powers are [256000, 338000, 128000] [J/m**3/s]. Issue #24 - # separately tracks the existing factor of two applied to those powers. + # The mass-rate latent powers are [256000, 338000, 128000] [J/m**3/s]. expected_condensed_temperature_rate = np.array([ - -512000.0 / 900000.0, - -676000.0 / 900000.0, -256000.0 / 900000.0, + -338000.0 / 900000.0, + -128000.0 / 900000.0, ]) # [K/s] np.testing.assert_allclose( diff --git a/tests/test_drying_latent_heat_factor.py b/tests/test_drying_latent_heat_factor.py new file mode 100644 index 0000000..50352ff --- /dev/null +++ b/tests/test_drying_latent_heat_factor.py @@ -0,0 +1,102 @@ +"""Regression coverage for Drying condensed-phase latent heat.""" + +from types import SimpleNamespace + +import numpy as np +import pytest + +pytest.importorskip("assimulo") +import PharmaPy.Drying_Model as drying_module +from PharmaPy.Drying_Model import Drying + +pytestmark = pytest.mark.assimulo + + +def test_condensed_energy_uses_single_latent_heat_factor(monkeypatch): + """Mass drying rates multiplied by latent heat already give full power.""" + dryer = Drying(number_nodes=4, supercrit_names=["nitrogen"]) + dryer.idx_volatiles = np.array([0, 2]) # component indices [-] + dryer.idx_supercrit = np.array([1]) # component indices [-] + dryer.porosity = 0.5 # [-] + dryer.rho_sol = 1000.0 # [kg/m**3] + dryer.cp_sol = 1000.0 # [J/kg/K] + dryer.rho_liq = np.full(4, 800.0) # [kg/m**3] + dryer.h_T_j = 0.0 # [W/m**2/K] + dryer.a_V = 1.0 # [m**2/m**3] + dryer.h_T_loss = 0.0 # [W/m**2/K] + dryer.cake_height = 1.0 # [m] + dryer.T_ambient = 298.15 # [K] + dryer.dz = np.ones(4) # [m] + + dryer.Vapor_1 = SimpleNamespace( + mw=np.array([28.0, 28.0, 28.0]), # [g/mol] + getCp=lambda temp, mass_frac, basis: np.full(4, 1000.0), # [J/kg/K] + getHeatVaporization=lambda temp, basis: np.array([2.0e6, 1.0e6]), # [J/kg] + ) + dryer.Liquid_1 = SimpleNamespace( + getCp=lambda temp, mass_frac, basis: np.full(4, 2000.0), # [J/kg/K] + ) + + monkeypatch.setattr( + drying_module, + "high_resolution_fvm", + lambda values, boundary_cond: np.zeros(values.size + 1), # [K] + ) + + temp_gas = np.full(4, 295.0) # [K], matches energy_balance wet-bulb default + temp_sol = np.array([299.0, 304.0, 309.0, 314.0]) # [K] + satur = np.full(4, 0.5) # [-] + y_gas = np.array([ + [0.10, 0.80, 0.10], + [0.20, 0.70, 0.10], + [0.30, 0.60, 0.10], + [0.40, 0.50, 0.10], + ]) # gas mass fractions [-] + x_liq = np.array([ + [0.0, 1.0], + [0.0, 1.0], + [0.0, 1.0], + [0.0, 1.0], + ]) # liquid volatile mass fractions [-] + dry_rate = np.array([ + [0.036, 0.0, 0.184], + [0.054, 0.0, 0.230], + [0.018, 0.0, 0.092], + [0.027, 0.0, 0.046], + ]) # [kg/m**3/s] + + _, dTcond_dt = dryer.energy_balance( + time=0.0, + temp_gas=temp_gas, + temp_sol=temp_sol, + satur=satur, + y_gas=y_gas, + x_liq=x_liq, + u_gas=np.zeros(4), # [m/s] + rho_gas=np.ones(4), # [kg/m**3] + dry_rate=dry_rate, + inputs={"temp": 295.0}, # [K] + ) + + expected_condensed_temperature_rate = np.array([ + -256000.0 / 900000.0, + -338000.0 / 900000.0, + -128000.0 / 900000.0, + -100000.0 / 900000.0, + ]) # [K/s] + np.testing.assert_allclose(dTcond_dt, expected_condensed_temperature_rate) + + _, drying_term, _, _ = dryer.energy_balance( + time=0.0, + temp_gas=temp_gas, + temp_sol=temp_sol, + satur=satur, + y_gas=y_gas, + x_liq=x_liq, + u_gas=np.zeros(4), # [m/s] + rho_gas=np.ones(4), # [kg/m**3] + dry_rate=dry_rate, + inputs={"temp": 295.0}, # [K] + return_terms=True, + ) + np.testing.assert_allclose(drying_term, -expected_condensed_temperature_rate)