diff --git a/biosteam/units/phase_equilibrium.py b/biosteam/units/phase_equilibrium.py index ea7b86aa..b148f70c 100644 --- a/biosteam/units/phase_equilibrium.py +++ b/biosteam/units/phase_equilibrium.py @@ -241,7 +241,6 @@ def _run(self, stacklevel=1, P=None, solvent=None, update=True, else: self.K = K_new self.IDs = IDs - class MultiStageEquilibrium(Unit): @@ -794,6 +793,7 @@ def hot_start(self): for i in partitions: i.T = T i.phi = phi + for j in i.outs: j.T = T else: P = self.P # TODO: Set up much better initial conditions @@ -867,10 +867,62 @@ def hot_start(self): for i in bottom_side_draws: for s in stages[i].splitters: s._run() for i in partitions: i.IDs = IDs - top_flow_rates = self._get_top_flow_rates(False) + top_flow_rates = self._get_top_flow_rates() return top_flow_rates - def get_vle_phase_ratios(self): + def get_energy_balance_temperatures(self): + # ENERGY BALANCE + # Hv1 + Cv1*(dT1) + Hl1 + Cl1*dT1 - Hv2 - Cv2*dT2 - Hl0 - Cl0 = Q1 + # (Cv1 + Cl1)dT1 - Cv2*dT2 - Cl0*dT0 = Q1 - Hv1 - Hl1 + Hv2 + Hl0 + stages = self.stages + for i in stages: + if i.partition.B is not None: + raise RuntimeError('cannot specify vapor fraction and solve for temperature') + N_stages = self.N_stages + a = np.zeros(N_stages) + b = a.copy() + c = a.copy() + d = a.copy() + Ts = np.zeros(N_stages) + stage_mid = stages[0] + stage_bottom = stages[1] + partition_mid = stage_mid.partition + partition_bottom = stage_bottom.partition + C_out = sum([i.C for i in partition_mid.outs]) + C_bottom = stage_bottom.outs[0].C + Ts[0] = partition_mid.T + b[0] = C_out + c[0] = -C_bottom + Q = partition_mid.Q or 0. + d[0] = Q + stage_mid.H_in - partition_mid.H_out + for i in range(1, N_stages-1): + stage_top = stage_mid + stage_mid = stage_bottom + stage_bottom = stages[i+1] + partition_mid = partition_bottom + partition_bottom = stage_bottom.partition + C_out = sum([i.C for i in partition_mid.outs]) + C_top = stage_top.outs[1].C + C_bottom = stage_bottom.outs[0].C + Ts[i] = partition_mid.T + a[i] = -C_top + b[i] = C_out + c[i] = -C_bottom + Q = partition_mid.Q or 0. + d[i] = Q + stage_mid.H_in - partition_mid.H_out + stage_top = stage_mid + stage_mid = stage_bottom + partition_mid = partition_bottom + Ts[-1] = partition_mid.T + C_out = sum([i.C for i in partition_mid.outs]) + C_top = stage_top.outs[1].C + a[-1] = -C_top + b[-1] = C_out + Q = partition_mid.Q or 0. + d[-1] = Q + stage_mid.H_in - partition_mid.H_out + return Ts + solve_TDMA(a, b, c, d) + + def get_energy_balance_phase_ratios(self): # ENERGY BALANCE # Intermediate stage: # Hv1*V1 + Hl1*L1 - Hv2*V2 - Hl0*L0 = Q1 @@ -913,8 +965,7 @@ def get_vle_phase_ratios(self): next_stage = stages[i - 1] # Going up vapor_next, liquid_next, *_ = next_stage.outs Q += liquid_next.H - try: adjacent_stages = [i for i in (next_stage, last_stage) if i] - except: breakpoint() + adjacent_stages = [i for i in (next_stage, last_stage) if i] other_feeds = [i for i in stage.ins if i.source not in adjacent_stages] Q += sum([i.H for i in other_feeds]) if vapor.isempty(): @@ -931,13 +982,21 @@ def get_vle_phase_ratios(self): last_stage = stage return phase_ratios - def _get_top_flow_rates(self, solve_vapor_fraction): + def update_energy_balance_vapor_fractions(self): + phase_ratios = self.get_energy_balance_phase_ratios() + for i, j in zip(self.stages, phase_ratios): + i.partition.phi = 1 if j == inf else j / (1 + j) + + def update_energy_balance_temperatures(self): + temperatures = self.get_energy_balance_temperatures() + for stage, T in zip(self.stages, temperatures): + partition = stage.partition + partition.T = T + for i in partition.outs: i.T = T + + def _get_top_flow_rates(self): stages = self.stages partitions = [i.partition for i in stages] - if solve_vapor_fraction: - phase_ratios = self.get_vle_phase_ratios() - for i, j in enumerate(phase_ratios): - stages[i].partition.phi = 1 if j == inf else j / (1 + j) phase_ratios = [] partition_coefficients = [] Ts = [] @@ -997,7 +1056,7 @@ def multistage_equilibrium_iter(self, top_flow_rates): P = self.P eq = 'vle' if self.multi_stream.phases[0] == 'g' else 'lle' if eq == 'vle': - for n, i in enumerate(stages): + for i in stages: mixer = i.mixer partition = i.partition mixer.outs[0].mix_from( @@ -1006,19 +1065,23 @@ def multistage_equilibrium_iter(self, top_flow_rates): partition._run(P=self.P, update=False, couple_energy_balance=False) T = partition.T - for i in partition.outs: i.T = T - new_top_flow_rates = self._get_top_flow_rates(True) + for i in (partition.outs + i.outs): i.T = T + self.update_energy_balance_vapor_fractions() + new_top_flow_rates = self._get_top_flow_rates() else: # LLE for i in stages: - i.mixer._run() + mixer = i.mixer partition = i.partition + mixer._run() + partition.T = mixer.outs[0].T partition._run(P=P, solvent=self.solvent_ID, update=False, couple_energy_balance=False) for i in stages: partition = i.partition T = partition.T - for i in partition.outs: i.T = T - new_top_flow_rates = self._get_top_flow_rates(False) + for j in (i.outs + partition.outs): j.T = T + self.update_energy_balance_temperatures() + new_top_flow_rates = self._get_top_flow_rates() mol = top_flow_rates.flatten() mol_new = new_top_flow_rates.flatten() mol_errors = abs(mol - mol_new) @@ -1074,6 +1137,7 @@ def solve_TDMA(a, b, c, d): # Tridiagonal matrix solver http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm http://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm) + """ n = d.shape[0] - 1 # number of equations minus 1 for i in range(n): @@ -1227,7 +1291,7 @@ def hot_start_bottom_flow_rates( d[i] += top_feed_flows[i] - top_flows[i + 1] * asplit[i + 1] return solve_LBDMA(a, b, d) -# @njit(cache=True) +@njit(cache=True) def flow_rates_for_multistage_equilibrium( phase_ratios, partition_coefficients,