Coverage for backend/idaes_service/solver/custom/water_tank_with_units.py: 74%
112 statements
« prev ^ index » next coverage.py v7.10.7, created at 2025-11-06 23:27 +0000
1# Water tank model from IDAES, with the units added. Ideally, we can remove this
2# once our tank units pr to IDAES is merged.
3#################################################################################
4# The Institute for the Design of Advanced Energy Systems Integrated Platform
5# Framework (IDAES IP) was produced under the DOE Institute for the
6# Design of Advanced Energy Systems (IDAES).
7#
8# Copyright (c) 2018-2024 by the software owners: The Regents of the
9# University of California, through Lawrence Berkeley National Laboratory,
10# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
11# University, West Virginia University Research Corporation, et al.
12# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md
13# for full copyright and license information.
14#################################################################################
15"""
16Watertank model
17The water tank has only one inlet and one outlet
19main assumptions:
211.- Heat loss is a variable given by the user (zero heat loss can be
22specified if adiabatic)
232.- Calculate pressure change due to gravity based on water level
24and contraction to downcomer
253.- Water level is either fixed for steady-state model or calculated for
26dynamic model
274.- Assume enthalpy_in = enthalpy_out + heat loss
295.- Subcooled water from economizer and saturated water from waterwall
30are mixed before entering the tank
32Created: November 04 2020
33"""
34# TODO: Missing docstrings
35# pylint: disable=missing-function-docstring
37# Import Pyomo libraries
38from pyomo.environ import value, Var, Reference, acos
39from pyomo.common.config import ConfigBlock, ConfigValue, In, Bool
41# Import IDAES cores
42from idaes.core import (
43 ControlVolume0DBlock,
44 declare_process_block_class,
45 MaterialBalanceType,
46 EnergyBalanceType,
47 MomentumBalanceType,
48 UnitModelBlockData,
49 useDefault,
50)
51from idaes.core.util.config import is_physical_parameter_block
52from idaes.core.solvers import get_solver
53from idaes.core.util.constants import Constants as const
55import idaes.logger as idaeslog
57__author__ = "Boiler Subsystem Team (J. Ma, D. Caballero, M. Zamarripa)"
58__version__ = "2.0.0"
60# Set up logger
61_log = idaeslog.getLogger(__name__)
64@declare_process_block_class("WaterTank")
65class WaterTankData(UnitModelBlockData):
66 """
67 Water Tank Unit Operation Class
68 """
70 CONFIG = UnitModelBlockData.CONFIG()
72 CONFIG.declare(
73 "tank_type",
74 ConfigValue(
75 default="simple_tank",
76 domain=In(
77 [
78 "simple_tank",
79 "rectangular_tank",
80 "vertical_cylindrical_tank",
81 "horizontal_cylindrical_tank",
82 ]
83 ),
84 description="Flag indicating the tank type",
85 doc="""Flag indicating the type of tank to be modeled, and
86then calculate the volume of the filled level consequently,
87**default** - simple_tank.
88**Valid values:** {
89**simple_tank** - use a general tank and provide the area,
90**rectangular_tank** - use a rectangular tank and provide the width and length,
91**vertical_cylindrical_tank** - use a vertical cylindrical tank
92and provide the diameter,
93**horizontal_cylindrical_tank** - use a horizontal cylindrical tank and
94provide the length and diameter.}""",
95 ),
96 )
97 CONFIG.declare(
98 "material_balance_type",
99 ConfigValue(
100 default=MaterialBalanceType.componentPhase,
101 domain=In(MaterialBalanceType),
102 description="Material balance construction flag",
103 doc="""Indicates what type of material balance should be constructed,
104**default** - MaterialBalanceType.componentPhase.
105**Valid values:** {
106**MaterialBalanceType.none** - exclude material balances,
107**MaterialBalanceType.componentPhase** - use phase component balances,
108**MaterialBalanceType.componentTotal** - use total component balances,
109**MaterialBalanceType.elementTotal** - use total element balances,
110**MaterialBalanceType.total** - use total material balance.}""",
111 ),
112 )
113 CONFIG.declare(
114 "energy_balance_type",
115 ConfigValue(
116 default=EnergyBalanceType.enthalpyTotal,
117 domain=In(EnergyBalanceType),
118 description="Energy balance construction flag",
119 doc="""Indicates what type of energy balance should be constructed,
120**default** - EnergyBalanceType.enthalpyTotal.
121**Valid values:** {
122**EnergyBalanceType.none** - exclude energy balances,
123**EnergyBalanceType.enthalpyTotal** - single ethalpy balance for material,
124**EnergyBalanceType.enthalpyPhase** - ethalpy balances for each phase,
125**EnergyBalanceType.energyTotal** - single energy balance for material,
126**EnergyBalanceType.energyPhase** - energy balances for each phase.}""",
127 ),
128 )
129 CONFIG.declare(
130 "momentum_balance_type",
131 ConfigValue(
132 default=MomentumBalanceType.pressureTotal,
133 domain=In(MomentumBalanceType),
134 description="Momentum balance construction flag",
135 doc="""Indicates what type of momentum balance should be constructed,
136**default** - MomentumBalanceType.pressureTotal.
137**Valid values:** {
138**MomentumBalanceType.none** - exclude momentum balances,
139**MomentumBalanceType.pressureTotal** - single pressure balance for material,
140**MomentumBalanceType.pressurePhase** - pressure balances for each phase,
141**MomentumBalanceType.momentumTotal** - single momentum balance for material,
142**MomentumBalanceType.momentumPhase** - momentum balances for each phase.}""",
143 ),
144 )
145 CONFIG.declare(
146 "has_heat_transfer",
147 ConfigValue(
148 default=False,
149 domain=Bool,
150 description="Heat transfer term construction flag",
151 doc="""Indicates whether terms for heat transfer should be constructed,
152**default** - False.
153**Valid values:** {
154**True** - include heat transfer terms,
155**False** - exclude heat transfer terms.}""",
156 ),
157 )
158 CONFIG.declare(
159 "has_pressure_change",
160 ConfigValue(
161 default=True,
162 domain=Bool,
163 description="Pressure change term construction flag",
164 doc="""Indicates whether terms for pressure change should be
165constructed,
166**default** - False.
167**Valid values:** {
168**True** - include pressure change terms,
169**False** - exclude pressure change terms.}""",
170 ),
171 )
172 CONFIG.declare(
173 "property_package",
174 ConfigValue(
175 default=useDefault,
176 domain=is_physical_parameter_block,
177 description="Property package to use for control volume",
178 doc="""Property parameter object used to define property calculations,
179**default** - useDefault.
180**Valid values:** {
181**useDefault** - use default package from parent model or flowsheet,
182**PhysicalParameterObject** - a PhysicalParameterBlock object.}""",
183 ),
184 )
185 CONFIG.declare(
186 "property_package_args",
187 ConfigBlock(
188 implicit=True,
189 description="Arguments to use for constructing property packages",
190 doc="""A ConfigBlock with arguments to be passed to a property block(s)
191and used when constructing these,
192**default** - None.
193**Valid values:** {
194see property package for documentation.}""",
195 ),
196 )
198 def build(self):
199 """
200 Begin building model (pre-DAE transformation).
201 Args:
202 None
203 Returns:
204 None
205 """
207 # Call UnitModel.build to setup dynamics
208 super().build()
210 # Build Control Volume
211 self.control_volume = ControlVolume0DBlock(
212 dynamic=self.config.dynamic,
213 has_holdup=self.config.has_holdup,
214 property_package=self.config.property_package,
215 property_package_args=self.config.property_package_args,
216 )
218 self.control_volume.add_geometry()
220 self.control_volume.add_state_blocks(has_phase_equilibrium=False)
222 self.control_volume.add_material_balances(
223 balance_type=self.config.material_balance_type
224 )
226 self.control_volume.add_energy_balances(
227 balance_type=self.config.energy_balance_type,
228 has_heat_transfer=self.config.has_heat_transfer,
229 )
231 self.control_volume.add_momentum_balances(
232 balance_type=self.config.momentum_balance_type, has_pressure_change=True
233 )
235 if self.config.has_holdup is True:
236 # enforce that holdup must start above zero. This fixes some problems with initialisation from steady state assumptions.
237 self.control_volume.material_holdup[0, :, :].setlb(0)
239 # Add Inlet and Outlet Ports
240 self.add_inlet_port()
241 self.add_outlet_port()
243 # Add object references
244 self.volume = Reference(self.control_volume.volume)
246 # Set references to balance terms at unit level
247 if ( 247 ↛ 253line 247 didn't jump to line 253 because the condition on line 247 was always true
248 self.config.has_heat_transfer is True
249 and self.config.energy_balance_type != EnergyBalanceType.none
250 ):
251 self.heat_duty = Reference(self.control_volume.heat)
253 if ( 253 ↛ 260line 253 didn't jump to line 260 because the condition on line 253 was always true
254 self.config.has_pressure_change is True
255 and self.config.momentum_balance_type != "none"
256 ):
257 self.deltaP = Reference(self.control_volume.deltaP)
259 # Set Unit Geometry and Holdup Volume
260 self._set_geometry()
262 # Construct performance equations
263 self._make_performance()
265 def _set_geometry(self):
266 """
267 Define the geometry of the unit as necessary
268 """
270 units = self.config.property_package.get_metadata().get_derived_units
272 if self.config.tank_type == "simple_tank": 272 ↛ 274line 272 didn't jump to line 274 because the condition on line 272 was never true
273 # Declare a variable for cross sectional area
274 self.tank_cross_sect_area = Var(
275 initialize=1.0,
276 doc="Cross-sectional area of the tank",
277 units=units("length") ** 2,
278 )
280 elif self.config.tank_type == "rectangular_tank": 280 ↛ 282line 280 didn't jump to line 282 because the condition on line 280 was never true
281 # Declare variables for width and length
282 self.tank_width = Var(
283 initialize=1.0, units=units("length"), doc="Width of the tank"
284 )
285 self.tank_length = Var(
286 initialize=1.0, units=units("length"), doc="Length of the tank"
287 )
289 elif ( 289 ↛ exitline 289 didn't return from function '_set_geometry' because the condition on line 289 was always true
290 self.config.tank_type == "horizontal_cylindrical_tank"
291 or "vertical_cylindrical_tank"
292 ):
293 # Declare a variable for diameter of the tank
294 self.tank_diameter = Var(
295 initialize=0.5, units=units("length"), doc="Inside diameter of the tank"
296 )
297 if self.config.tank_type == "horizontal_cylindrical_tank": 297 ↛ 299line 297 didn't jump to line 299 because the condition on line 297 was never true
298 # Declare a variable for length of the tank
299 self.tank_length = Var(
300 initialize=1, units=units("length"), doc="Length of the tank"
301 )
303 def _make_performance(self):
304 """
305 Define constraints which describe the behaviour of the unit model
306 """
307 units = self.config.property_package.get_metadata().get_derived_units
309 # Add performance variables
310 self.tank_level = Var(
311 self.flowsheet().time,
312 initialize=1.0,
313 doc="Water level from in the tank",
314 units=units("length"),
315 )
317 # Auxiliary expressions for volume
318 # Rectangular tank
319 if self.config.tank_type == "rectangular_tank": 319 ↛ 321line 319 didn't jump to line 321 because the condition on line 319 was never true
320 # Calculation of cross-sectional area of the rectangle
321 @self.Expression(doc="Cross-sectional area of the tank")
322 def tank_cross_sect_area(b):
323 return b.tank_width * b.tank_length
325 # Vertical cylindrical tank
326 elif self.config.tank_type == "vertical_cylindrical_tank": 326 ↛ 338line 326 didn't jump to line 338 because the condition on line 326 was always true
328 @self.Expression(doc="Radius of the tank")
329 def tank_radius(b):
330 return b.tank_diameter / 2
332 # Calculation of cross-sectional area of the vertical cylinder
333 @self.Expression(doc="Cross-sectional area of the tank")
334 def tank_cross_sect_area(b):
335 return const.pi * b.tank_radius**2
337 # Horizontal cylindrical tank
338 elif self.config.tank_type == "horizontal_cylindrical_tank":
339 # Calculation of area covered by the liquid level
340 # at one end of the tank
341 @self.Expression(doc="Radius of the tank")
342 def tank_radius(b):
343 return b.tank_diameter / 2
345 # Angle of the circular sector used to calculate the area
347 @self.Expression(
348 self.flowsheet().time,
349 doc="Angle of the circular" " sector of liquid level",
350 )
351 def alpha_tank(b, t):
352 return 2 * acos((b.tank_radius - b.tank_level[t]) / b.tank_radius)
354 @self.Expression(
355 self.flowsheet().time,
356 doc="Area covered by the liquid level" " at one end of the tank",
357 )
358 def tank_area(b, t):
359 return (
360 0.5 * b.alpha_tank[t] * b.tank_radius**2
361 - (b.tank_radius - b.tank_level[t])
362 * (2 * b.tank_radius * b.tank_level[t] - b.tank_level[t] ** 2)
363 ** 0.5
364 )
366 # Constraint for volume of the liquid in tank
367 @self.Constraint(self.flowsheet().time, doc="volume of liquid in the tank")
368 def volume_eqn(b, t):
369 if self.config.tank_type == "horizontal_cylindrical_tank": 369 ↛ 370line 369 didn't jump to line 370 because the condition on line 369 was never true
370 return b.volume[t] == b.tank_length * b.tank_area[t]
371 else:
372 return b.volume[t] == b.tank_level[t] * b.tank_cross_sect_area
374 # Pressure change equation due gravity
375 @self.Constraint(self.flowsheet().time, doc="pressure drop")
376 def pressure_change_eqn(b, t):
377 return (
378 b.deltaP[t]
379 == b.control_volume.properties_in[t].dens_mass_phase["Liq"]
380 * const.acceleration_gravity
381 * b.tank_level[t]
382 )
384 def set_initial_condition(self):
385 if self.config.dynamic is True:
386 self.control_volume.material_accumulation[:, :, :].value = 0
387 self.control_volume.energy_accumulation[:, :].value = 0
388 self.control_volume.material_accumulation[0, :, :].fix(0)
389 self.control_volume.energy_accumulation[0, :].fix(0)
391 def initialize_build(
392 blk, state_args=None, outlvl=idaeslog.NOTSET, solver=None, optarg=None
393 ):
394 """
395 Water tank initialization routine.
397 Keyword Arguments:
398 state_args : a dict of arguments to be passed to the property
399 package(s) for the control_volume of the model to
400 provide an initial state for initialization
401 (see documentation of the specific property package)
402 (default = None).
403 outlvl : sets output level of initialisation routine
405 * 0 = no output (default)
406 * 1 = return solver state for each step in routine
407 * 2 = return solver state for each step in subroutines
408 * 3 = include solver output information (tee=True)
410 optarg : solver options dictionary object (default=None, use
411 default solver options)
412 solver : str indicating which solver to use during
413 initialization (default = None, use default solver)
415 Returns:
416 None
417 """
418 init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit")
419 solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit")
421 # Create solver
422 opt = get_solver(solver, optarg)
424 init_log.info_low("Starting initialization...")
426 flags = blk.control_volume.initialize(
427 state_args=state_args, outlvl=outlvl, optarg=optarg, solver=solver
428 )
429 init_log.info_high("Initialization Step 1 Complete.")
431 # Fix outlet pressure
432 for t in blk.flowsheet().time:
433 blk.control_volume.properties_out[t].pressure.fix(
434 value(blk.control_volume.properties_in[t].pressure)
435 )
436 blk.pressure_change_eqn.deactivate()
438 # solve model
439 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
440 res = opt.solve(blk, tee=slc.tee)
441 init_log.info_high("Initialization Step 2 {}.".format(idaeslog.condition(res)))
443 # Unfix outlet pressure
444 for t in blk.flowsheet().time:
445 blk.control_volume.properties_out[t].pressure.unfix()
446 blk.pressure_change_eqn.activate()
448 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
449 res = opt.solve(blk, tee=slc.tee)
450 init_log.info_high("Initialization Step 3 {}.".format(idaeslog.condition(res)))
452 blk.control_volume.release_state(flags, outlvl)
453 init_log.info("Initialization Complete.")
455 def calculate_scaling_factors(self):
456 pass