Skip to content

Commit

Permalink
finally get no/full condenser/ 0 or inf reflux working
Browse files Browse the repository at this point in the history
  • Loading branch information
cortespea committed Dec 17, 2023
1 parent 5422404 commit afb95f7
Show file tree
Hide file tree
Showing 2 changed files with 80 additions and 18 deletions.
42 changes: 42 additions & 0 deletions biosteam/units/distillation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2259,6 +2259,48 @@ class MESHDistillation(MultiStageEquilibrium, new_graphics=False):
Total purchase cost USD 1.16e+05
Utility cost USD/hr 91.3
Simulate distillation column with a full condenser, 5 stages, a 0.673 reflux ratio,
2.57 boilup ratio, and feed at stage 2:
>>> import biosteam as bst
>>> bst.settings.set_thermo(['Water', 'Ethanol'], cache=True)
>>> feed = bst.Stream('feed', Ethanol=80, Water=100, T=80.215 + 273.15)
>>> D1 = bst.MESHDistillation(None, N_stages=5, ins=[feed], feed_stages=[2],
... outs=['vapor', 'liquid', 'distillate'],
... reflux=0.673, boilup=2.57,
... LHK=('Ethanol', 'Water'),
... full_condenser=True,
... )
>>> D1.simulate()
>>> vapor, liquid, distillate = D1.outs
>>> distillate.imol['Ethanol'] / feed.imol['Ethanol']
0.81
>>> distillate.imol['Ethanol'] / distillate.F_mol
0.70
>>> D1.results()
Distillation Units
Electricity Power kW 0.874
Cost USD/hr 0.0683
Low pressure steam Duty kJ/hr 1.28e+07
Flow kmol/hr 330
Cost USD/hr 78.5
Design Theoretical stages 5
Actual stages 6
Height ft 22.9
Diameter ft 3.82
Wall thickness in 0.312
Weight lb 4e+03
Purchase cost Trays USD 7.58e+03
Tower USD 3.62e+04
Platform and ladders USD 9.8e+03
Condenser - Double pipe USD 5.19e+03
Pump - Pump USD 4.35e+03
Pump - Motor USD 387
Reboiler - Floating head USD 2.3e+04
Total purchase cost USD 8.65e+04
Utility cost USD/hr 78.6
Notes
-----
The convergence algorithm decouples the equilibrium relationships,
Expand Down
56 changes: 38 additions & 18 deletions biosteam/units/phase_equilibrium.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -161,6 +161,7 @@ def _run(self, stacklevel=1, P=None, solvent=None, update=True,
else:
ms = self.feed.copy()
ms.phases = self.phases
if ms.isempty(): return
top, bottom = ms
partition_data = self.partition_data
if partition_data:
Expand Down Expand Up @@ -369,8 +370,9 @@ class MultiStageEquilibrium(Unit):
>>> vapor.imol['AceticAcid'] / feed.imol['AceticAcid']
0.74
# Simulate distillation column with 5 stages, a 0.673 reflux ratio,
# 2.57 boilup ratio, and feed at stage 2:
Simulate distillation column with 5 stages, a 0.673 reflux ratio,
2.57 boilup ratio, and feed at stage 2:
>>> import biosteam as bst
>>> bst.settings.set_thermo(['Water', 'Ethanol'], cache=True)
>>> feed = bst.Stream('feed', Ethanol=80, Water=100, T=80.215 + 273.15)
Expand All @@ -386,6 +388,24 @@ class MultiStageEquilibrium(Unit):
>>> vapor.imol['Ethanol'] / vapor.F_mol
0.69
Simulate the same distillation column with a full condenser, 5 stages, a 0.673 reflux ratio,
2.57 boilup ratio, and feed at stage 2:
>>> import biosteam as bst
>>> bst.settings.set_thermo(['Water', 'Ethanol'], cache=True)
>>> feed = bst.Stream('feed', Ethanol=80, Water=100, T=80.215 + 273.15)
>>> MSE = bst.MultiStageEquilibrium(N_stages=5, ins=[feed], feed_stages=[2],
... outs=['vapor', 'liquid', 'distillate'],
... stage_specifications={0: ('Reflux', float('inf')), -1: ('Boilup', 2.57)},
... bottom_side_draws={0: 0.673 / (1 + 0.673)}
... )
>>> MSE.simulate()
>>> vapor, liquid, distillate = MSE.outs
>>> distillate.imol['Ethanol'] / feed.imol['Ethanol']
0.81
>>> distillate.imol['Ethanol'] / distillate.F_mol
0.70
"""
_N_ins = 2
_N_outs = 2
Expand Down Expand Up @@ -567,8 +587,8 @@ def update(self, top_flow_rates):
bottom_flow = sum([i.mol for i in stage.ins]) - s_top.mol
mask = bottom_flow < 0.
if mask.any():
has_infeasible_flow = (bottom_flow[mask] < flow_tol[mask]).any() and i not in infeasible_checks
if has_infeasible_flow:
has_infeasible_flow = (bottom_flow[mask] < flow_tol[mask]).any()
if i not in infeasible_checks and has_infeasible_flow:
infeasible_checks.add(i)
infeasible_index, = np.where(mask[index])
# TODO: Find algebraic solution to keeping top flow rates within feasible region.
Expand All @@ -577,6 +597,7 @@ def update(self, top_flow_rates):
top_flow_rates[i, infeasible_index] += infeasible_flow
break
else:
has_infeasible_flow = False
bottom_flow[mask] = 0.
s_bottom.mol[:] = bottom_flow
if stage.bottom_split: stage.splitters[-1]._run()
Expand Down Expand Up @@ -857,7 +878,7 @@ def get_vle_phase_ratios(self):
Q = partition.Q or 0
vapor, liquid, *_ = stage.outs
Q -= liquid.H
if last_stage:
if last_stage and B_last:
Q += vapor_last.h * liquid_last.F_mol * B_last
if i != 0:
next_stage = stages[i - 1] # Going up
Expand All @@ -866,10 +887,10 @@ def get_vle_phase_ratios(self):
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 liquid.isempty():
B = B_max
elif vapor.isempty():
if vapor.isempty():
B = B_min
elif liquid.isempty():
B = B_max
else:
B = Q / (vapor.h * liquid.F_mol)
if B > B_max: B = B_max
Expand Down Expand Up @@ -908,7 +929,7 @@ def _get_top_flow_rates(self, inner_vle_loop):
almost_one = 1. - 1e-16
almost_zero = 1e-16
for i in range(N_stages):
partition = stages[i].partition
partition = partitions[i]
phi = partition.phi
if phi is not None and (partition.B is not None or almost_zero < phi < almost_one):
index.append(i)
Expand Down Expand Up @@ -966,7 +987,7 @@ def multistage_equilibrium_iter(self, top_flow_rates):
mixer.ins, conserve_phases=True, energy_balance=False,
)
partition._run(P=self.P, update=False,
couple_energy_balance=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)
Expand All @@ -976,8 +997,8 @@ def multistage_equilibrium_iter(self, top_flow_rates):
i.partition._run(P=P, solvent=self.solvent_ID, update=False,
couple_energy_balance=False)
new_top_flow_rates = self._get_top_flow_rates(False)
mol = np.hstack([top_flow_rates[0], top_flow_rates[-1]])
mol_new = np.hstack([new_top_flow_rates[0], new_top_flow_rates[-1]])
mol = top_flow_rates.flatten()
mol_new = new_top_flow_rates.flatten()
mol_errors = abs(mol - mol_new)
if mol_errors.any():
mol_error = mol_errors.max()
Expand Down Expand Up @@ -1201,7 +1222,6 @@ def flow_rates_for_multistage_equilibrium(
inf_mask = phase_ratios >= 1e32
ok_mask = ~zero_mask & ~inf_mask
phase_ratios = np.expand_dims(phase_ratios, -1)
asplit = np.expand_dims(asplit, -1)
safe = ok_mask.all()
if safe:
component_ratios = 1. / (phase_ratios * partition_coefficients)
Expand All @@ -1217,9 +1237,9 @@ def flow_rates_for_multistage_equilibrium(
for i in inf_index:
component_ratios[i] = 0.
b = 1. + component_ratios
c = bsplit[1:]
c = asplit[1:]
d = feed_flows.copy()
a = asplit * component_ratios
a = np.expand_dims(bsplit, -1) * component_ratios
if safe:
return solve_TDMA(a, b, c, d)
else:
Expand All @@ -1235,9 +1255,10 @@ def flow_rates_for_multistage_equilibrium(
ok_index, = np.nonzero(ok_mask)
inf_index, = np.nonzero(inf_mask)
zero_index, = np.nonzero(zero_mask)
special_index, = np.nonzero(inf_mask & (ai == inf))
special_index, = np.nonzero(inf_mask & (ai == -inf))
special = bsplit[i]
for j in inf_index: m[j] = 0
for j in special_index: m[j] = asplit[j]
for j in special_index: m[j] = special
for j in zero_index: m[j] = inf
for j in ok_index: m[j] = ai[j] / bi[j]
b[inext] = b[inext] - m * c[i]
Expand All @@ -1247,7 +1268,6 @@ def flow_rates_for_multistage_equilibrium(

for i in range(n-1, -1, -1):
b[i] = (d[i] - c[i] * b[i+1]) / b[i]

return b

def get_neighbors(index, all_index):
Expand Down

0 comments on commit afb95f7

Please sign in to comment.