Barotropic model

class barotropy.barotropic_model.BarotropicModel(fluid_name, p_in: float, T_in: float = None, rho_in: float = None, p_out: float = None, p_min: float = None, p_max: float = None, efficiency: float = 1.0, mixture_ratio: float = None, process_type: str = None, calculation_type: str = None, blending_onset: float = None, blending_width: float = None, HEOS_solver: str = 'hybr', HEOS_tolerance: float = 1e-06, HEOS_max_iter: int = 100, HEOS_print_convergence: bool = False, ODE_solver: str = 'LSODA', ODE_tolerance: float = 1e-08, polynomial_degree: int = 8, polynomial_format: str = 'horner', polynomial_variables: list = ['density', 'viscosity', 'speed_sound', 'void_fraction', 'vapor_quality'], output_dir: str = 'barotropic_model')

Bases: object

Coordinates the simulation and polynomial fitting for a barotropic process.

Parameters:

fluid_namestr or list of str

The name(s) of the fluid(s) for the barotropic model.

• Specify a single string (e.g., ‘co2’) to use the single-component model.

• Specify a list of two strings (e.g., [‘water’, ‘nitrogen’]) to use the two-component model.

T_infloat

Inlet temperature of the fluid in Kelvin.

p_infloat

Inlet pressure of the fluid in Pascals.

rho_infloat

Inlet density of the fluid in kilogram per cubic meter.

Note

Applicable only to the one-component model.

p_min: float

Minimum pressure for the barotropic process integration in Pascals.

p_maxfloat

Maximum pressure for the barotropic process integration in Pascals.

efficiencyfloat

The efficiency of the polytropic process, dimensionless.

mixture_ratiofloat

Mass ratio of the first to the second fluid in the mixture.

Note

Applicable only to the two-component model.

process_typestr, optional

The type of polytropic process that the fluid experiences. Must be ‘expansion’ or ‘compression’.

calculation_typestr, optional

The type of fluid property calculation for the one-component model. Options include:

• equilibrium: Computes equilibrium properties only.

• metastable: Computes metastable properties only.

• blending: Computes both equilibrium and metastable properties, blends them, and returns the blended properties.

Note

Applicable only to the one-component model.

blending_onsetfloat, optional

The onset of blending in the process, typically a value between 0 and 1. Required when calculation_type is blending.

Note

Applicable only to the one-component model.

blending_widthfloat, optional

The width of the blending region, typically a value between 0 and 1. Required when calculation_type is blending.

Note

Applicable only to the one-component model.

HEOS_solverstr, optional

The solver algorithm used to compute the metastable states. Valid options:

Solver name	Description
hybr	Powell’s hybrid trust-region algorithm
lm	Levenberg–Marquardt algorithm

See Scipy root() for more info. Recommended solvers: Both hybr and lm work well for the tested cases.

Note

Applicable only to the one-component model.

HEOS_tolerancefloat, optional

The tolerance for the HEOS solver.

Note

Applicable only to the one-component model.

HEOS_max_iterint, optional

The maximum number of iterations for the HEOS solver.

Note

Applicable only to the one-component model.

HEOS_print_convergencebool, optional

If True, prints convergence information for the HEOS solver.

Note

Applicable only to the one-component model.

ODE_solverstr, optional

The solver to use for the ODE integration. Valid options:

Solver name	Description
RK23	Explicit Runge-Kutta method of order 3(2)
RK45	Explicit Runge-Kutta method of order 5(4)
DOP853	Explicit Runge-Kutta method of order 8
Radau	Implicit Runge-Kutta method of the Radau IIA family of order 5
BDF	Implicit multi-step variable-order (1 to 5) method based on a backward differentiation formula for the derivative approximation
LSODA	Adams/BDF method with automatic stiffness detection and switching

See Scipy solver_ivp() for more info. Recommended solvers: BDF, LSODA, or Radau for stiff problems or RK45 for non-stiff problems with smooth blending.

ODE_tolerancefloat, optional

The relative and absolute tolerance for the ODE solver.

polynomial_degreeint

Degree of the polynomials to fit.

Note

When calculation_type is blending the degree of the polynomial in the blending region is set to 4 to achieve sufficient accuracy while preventing numerical round-off errors associated with single-precision arithmetic in CFD solvers.

polynomial_formatstr, optional

Type of polynomial representation (horner or standard).

polynomial_variableslist of str

A list of variable names to fit polynomials to, such as ‘density’, ‘viscosity’, ‘speed_sound’, ‘void_fraction’, ‘vapor_quality’.

output_dirstr

The directory where output will be saved.

Methods

export_expressions_cfx([output_dir])	Exports the polynomial expressions in a format suitable for Ansys CFX.
export_expressions_fluent([output_dir])	Exports the polynomial expressions in a format suitable for Ansys Fluent.
fit_polynomials()	Fits polynomials to the states using the PolynomialFitter class.
solve()	Solves the equations for the one-component or two-component barotropic model and stores the fluid properties.

export_expressions_cfx(output_dir=None)

Exports the polynomial expressions in a format suitable for Ansys CFX.

Parameters:

output_dirstr, optional

The directory where the expressions will be saved. It uses a default directory if not provided.

See also

ExpressionExporter

Class for exporting polynomial expressions for use in CFD software.

export_expressions_fluent(output_dir=None)

Exports the polynomial expressions in a format suitable for Ansys Fluent.

Parameters:

output_dirstr, optional

The directory where the expressions will be saved. It uses a default directory if not provided.

See also

ExpressionExporter

Class for exporting polynomial expressions for use in CFD software.

fit_polynomials()

Fits polynomials to the states using the PolynomialFitter class.

See also

PolynomialFitter

Class used for generating fitting polynomials.

solve()

Solves the equations for the one-component or two-component barotropic model and stores the fluid properties. The type of calculation performed is selected automatically depending on the number of fluid names defined when initializing the class.

See also

barotropic_model_one_component

Calculation of fluid properties for the one-component model.

barotropic_model_two_component

Calculation of fluid properties for the two-component model.

class barotropy.barotropic_model.ExpressionExporter(poly_fitter, poly_format='horner')

Bases: object

Exports expressions describing the barotropic model for use within CFD solvers.

Parameters:

poly_fitterPolynomialFitter

An instance of PolynomialFitter containing the polynomial coefficients and breakpoints.

poly_formatstr, optional

Type of polynomial representation (‘horner’ or ‘standard’).

Methods

export_expressions_cfx([output_dir])	Exports the barotropic model polynomials as expressions for use in CFX CFD software.
export_expressions_fluent([output_dir])	Exports the barotropic model polynomials as expressions for use in Fluent CFD software.
generate_expressions(var, unit)	Generates a piecewise polynomial expression with extrapolation safeguards for a given variable.
if_expression(expression_1, expression_2, ...)	Generates an IF expression string for piecewise functions.
polynomial_expression(coefficients)	Generates a polynomial expression string.

export_expressions_cfx(output_dir=None)

Exports the barotropic model polynomials as expressions for use in CFX CFD software.

Parameters:

output_dirstr, optional

The directory where the expressions will be saved. It uses a default directory if not provided.

export_expressions_fluent(output_dir=None)

Exports the barotropic model polynomials as expressions for use in Fluent CFD software.

Parameters:

output_dirstr, optional

The directory where the expressions will be saved. It uses a default directory if not provided.

generate_expressions(var, unit)

Generates a piecewise polynomial expression with extrapolation safeguards for a given variable.

Parameters:

varstr

The variable name for which to generate the expression.

unitstr

The unit of the variable.

Returns:

str

The piecewise polynomial expression with extrapolation safeguards.

if_expression(expression_1, expression_2, transition_pressure)

Generates an IF expression string for piecewise functions.

Parameters:

expression_1str

The expression for the condition being true.

expression_2str

The expression for the condition being false.

transition_pressurefloat

The transition pressure where the expressions change.

Returns:

str

The IF expression string.

polynomial_expression(coefficients)

Generates a polynomial expression string.

Parameters:

coefficientsarray-like

Coefficients of the polynomial.

Returns:

str

The polynomial expression string.

class barotropy.barotropic_model.PolynomialFitter(states, state_in, variables, degree, calculation_type, model_type, output_dir='barotropic_model')

Bases: object

# TODO update docstring Fits polynomials to the thermodynamic properties of a fluid across various states.

polynomial_degreeint

Degree of the polynomials to fit.

Note

When calculation_type is blending the degree of the polynomial in the blending region is set to 4 to achieve sufficient accuracy while preventing numerical round-off errors associated with single-precision arithmetic in CFD solvers.

polynomial_formatstr, optional

Type of polynomial representation (horner or standard). Default is ‘horner’.

polynomial_variableslist of str

A list of variable names to fit polynomials to, such as ‘density’, ‘viscosity’, ‘speed_sound’, ‘void_fraction’, ‘vapor_quality’.

output_dirstr

The directory where output will be saved.

Methods

evaluate_polynomial(p, variable)	Evaluates the polynomials for a given variable at a specified pressure values.
fit_polynomials()	Fits polynomials to the data obtained from the ODE solution based on the specified calculation type.
plot_phase_diagram(fluid, var_x, var_y[, ...])	Plots the barotropic process in the phase diagram of the specified fluid.
plot_polynomial(var[, p_eval, showfig, ...])	Plots the polynomials for the specified variable.
plot_polynomial_and_error(var[, num_points, ...])	Plots the barotropic polynomials and original data points from the ODE solution to illustrate the quality of the fit.

evaluate_polynomial(p, variable)

Evaluates the polynomials for a given variable at a specified pressure values.

The function evaluates the polynomial values at a given physical pressure p for a variable (e.g., density). It automatically determines the correct branch of the piecewise polynomial based on the breakpoints, and applies safeguards for pressures outside the limits to ensure sensible values and continuity.

The evaluation follows a piecewise definition depending on the range of \(\hat{p}=p/p_{\text{in}}\):

\[\begin{split}\phi(\hat{p}) = \begin{cases} a_1 \, e^\left(\frac{\hat{p} - \hat{p}_{\text{out}}}{a_2}\right) & \text{ for } & \hat{p} < \hat{p}_{\text{out}} \\ \sum_{i=0}^{d} b_{i,1} \, \hat{p}^i & \text{ for } & \hat{p}_{\text{ out}} \leq \hat{p} \leq \hat{p}_{1} \\ & \;\; \vdots & \\ \sum_{i=0}^{d} b_{i,n} \, \hat{p}^i & \text{ for } & \hat{p}_{n} \leq \hat{p} \leq \hat{p}_{\text{in}} \\ c_1 + c_2 \left(\hat{p} - \hat{p}_{\text{in}}\right) & \text{ for } & \hat{p} > \hat{p}_{\text{in}} \end{cases}\end{split}\]

where:

• \(\phi(p)\) is the variable being fitted (e.g., density, viscosity).

• \(a_1\) and \(a_2\) are constants of the exponential function used to extrapolate the properties below the outlet pressure. The numerical values are determined to match the polynomial value and its first derivative at the endpoint \(\hat{p}=\hat{p}_{\text{out}}\).

• \(b_{i,j}\) is the \(i\)-th polynomial coefficient of the \(j\)-th polynomial segment, where \(d\) is the degree of the polynomials, \(n\) is the number of breakpoints, and \(n+1\) is the number of polynomial segments. Additionally, \([\hat{p}_{\text{out}},\, \hat{p}_{1},\,\ldots,\, \hat{p}_{n},\,\hat{p}_{\text{in}}]\) are normalized pressures at the breakpoints.

• \(c_1\) and \(c_2\) are constants of the linear function used to extrapolate the properties above the inlet pressure. The numerical values are determined to match the polynomial value and its first derivative at the endpoint \(\hat{p}=\hat{p}_{\text{in}}\).

Note

The safeguards for pressures outside the specified limits help to avoid numerical issues by applying exponential decay for pressures below \(p_{\text{out}}\) and linear extrapolation for pressures above \(p_{\text{in}}\). These measures prevent properties from becoming negative or excessively large, which can occur with very high or low (even negative) pressures during internal iterations of CFD solvers.

Parameters:

parray-like

The physical pressure values at which to evaluate the polynomials. Given in Pascals (Pa).

variablestr

The variable for which to evaluate the polynomial, such as ‘density’, ‘viscosity’, ‘speed_of_sound’, or ‘void_fraction’.

Returns:

var_valuesndarray

The evaluated values of the variable at the given physical pressures.

fit_polynomials()

Fits polynomials to the data obtained from the ODE solution based on the specified calculation type.

plot_phase_diagram(fluid, var_x, var_y, savefig=True, showfig=True, output_dir=None, plot_spinodal_line=True, plot_quality_isolines=True)

Plots the barotropic process in the phase diagram of the specified fluid.

Note

This function is applicable to single-component systems

Parameters:

fluidFluid

Fluid object used to plot the phase diagram.

var_xstr

Variable for the x-axis (e.g., “s” for entropy).

var_ystr

Variable for the y-axis (e.g., “T” for temperature).

savefigbool, optional

If True, saves the figure to the specified directory. Default is True.

showfigbool, optional

If True, displays the figure after plotting. Default is True.

output_dirstr, optional

Directory where the plot will be saved if savefig is True. If not specified, uses the default output directory of the class.

Returns:

figmatplotlib.figure.Figure

The figure object containing the plotted phase diagram.

plot_polynomial(var, p_eval=None, showfig=True, savefig=False, output_dir=None)

Plots the polynomials for the specified variable. If the pressure values for plotting are not provided, a suitable range is generated internally.

Parameters:

varstr

The variable for which to plot the polynomial.

p_evalarray-like, optional

The pressure values at which to plot the polynomial. If not provided, automatic values are generated.

showfigbool, optional

If True, displays the plot.

savefigbool, optional

If True, saves the plot to a file.

output_dirstr, optional

The directory where the plot will be saved if savefig is True. If not specified, uses the default output directory of the class.

Returns:

figmatplotlib.figure.Figure

The figure object containing the plotted data.

plot_polynomial_and_error(var, num_points=500, showfig=True, savefig=False, output_dir=None)

Plots the barotropic polynomials and original data points from the ODE solution to illustrate the quality of the fit. Additionally, plots the relative error between the polynomial and the data.

Parameters:

varstr

The variable to plot (e.g., ‘density’, ‘viscosity’, etc.).

num_pointsint, optional

Number of plot points for each of the polynomial segments.

showfigbool, optional

If True, displays the plot.

savefigbool, optional

If True, saves the plot to a file.

output_dirstr, optional

The directory where the plot will be saved if savefig is True. If not specified, uses the default output directory of the class.

Returns:

figmatplotlib.figure.Figure

The figure object containing the plotted data and error.

barotropy.barotropic_model.barotropic_model_one_component(fluid_name, p_in, rho_in, p_out, efficiency, process_type=None, calculation_type=None, blending_onset=None, blending_width=None, HEOS_solver='hybr', HEOS_tolerance=1e-06, HEOS_max_iter=100, HEOS_print_convergence=False, ODE_solver='lsoda', ODE_tolerance=1e-08)

Simulates a polytropic process for a given fluid from an inlet state to a specified outlet pressure.

This function integrates the differential equation describing the polytropic process, which may involve phase changes, using the relationship:

\[\frac{\text{d}h}{\text{d}p} = \frac{\eta_p}{\rho}\]

where \(h\) is the specific enthalpy, \(p\) is the pressure, \(\eta_p\) is the process efficiency, and \(\rho\) is the density of the fluid.

The evaluation of fluid properties varies based on the calculation_type parameter:

• calculation_type='equilibrium': Computes properties assuming thermodynamic equilibrium using the CoolProp property solver.

• calculation_type='metastable': Computes properties based on the Helmholtz energy equations of state, iterating on density and temperature.

• calculation_type='blending': Blends equilibrium and metastable properties according to the fluid vapor quality.

In single-phase regions, all calculation types yield identical results. The differences arise in the evaluation within the two-phase region, where the handling of equilibrium and metastable fluid properties differs.

Note

The choice of ODE solver can significantly impact the accuracy and stability of the simulation. The BDF and LSODA solvers are recommended for general use, particularly for stiff problems involving equilibrium phase change or narrow blending widths, while RK45 may be suitable for non-stiff problems with no phase change or smooth phase transitions.

Parameters:

fluid_namestr

The name of the fluid to be used in the simulation (e.g., ‘Water’).

p_infloat

Inlet pressure of the fluid in Pascal.

rho_infloat

Inlet density of the fluid in kilogram per cubic meter.

p_outfloat

Outlet pressure of the fluid in Pascals.

efficiencyfloat

The efficiency of the polytropic process, dimensionless.

calculation_typestr

The type of calculation to perform. Options include:

• equilibrium: Computes equilibrium properties only.

• metastable: Computes metastable properties only.

• blending: Computes both equilibrium and metastable properties, blends them, and returns the blended properties.

blending_onsetfloat, optional

The onset of blending in the process, typically a value between 0 and 1. Required when calculation_type is blending.

blending_widthfloat, optional

The width of the blending region, typically a value between 0 and 1. Required when calculation_type is blending.

HEOS_solverstr, optional

The solver algorithm used to compute the metastable states. Valid options:

Solver name	Description
hybr	Powell’s hybrid trust-region algorithm
lm	Levenberg–Marquardt algorithm

See Scipy root() for more info. Recommended solvers: Both hybr and lm work well for the tested cases.

HEOS_tolerancefloat, optional

The tolerance for the HEOS solver.

HEOS_max_iterint, optional

The maximum number of iterations for the HEOS solver.

HEOS_print_convergencebool, optional

If True, prints convergence information for the HEOS solver.

ODE_solverstr, optional

The solver to use for the ODE integration. Valid options:

Solver name	Description
RK23	Explicit Runge-Kutta method of order 3(2)
RK45	Explicit Runge-Kutta method of order 5(4)
DOP853	Explicit Runge-Kutta method of order 8
Radau	Implicit Runge-Kutta method of the Radau IIA family of order 5
BDF	Implicit multi-step variable-order (1 to 5) method based on a backward differentiation formula for the derivative approximation
LSODA	Adams/BDF method with automatic stiffness detection and switching

See Scipy solver_ivp() for more info. Recommended solvers: BDF, LSODA, or Radau for stiff problems or RK45 for non-stiff problems with smooth blending.

ODE_tolerancefloat, optional

The relative and absolute tolerance for the ODE solver.

Returns:

statesdictionary of arrays

A dictionary of Numpy arrays representing the properties of the fluid at each evaluation point.

solutionscipy.integrate.OdeResult

The result of the ODE integration containing information about the solver process.

barotropy.barotropic_model.barotropic_model_two_component(fluid_name_1, fluid_name_2, mixture_ratio, T_in, p_in, p_out, efficiency, process_type=None, ODE_solver='lsoda', ODE_tolerance=1e-08)

Simulates a polytropic process for a mixture of two different fluids.

TODO: add model equations and explanation

Parameters:

fluid_name_1str

The name of the first component of the mixture.

fluid_name_2str

The name of the second component of the mixture.

mixture_ratiofloat

Mass ratio of the first to the second fluid in the mixture.

T_infloat

Inlet temperature of the mixture in Kelvin.

p_infloat

Inlet pressure of the mixture in Pascals.

p_outfloat

Outlet pressure of the mixture in Pascals.

efficiencyfloat

The efficiency of the polytropic process, (between zero and one).

ODE_solverstr, optional

The solver to use for the ODE integration. Valid options:

Solver name	Description
RK23	Explicit Runge-Kutta method of order 3(2)
RK45	Explicit Runge-Kutta method of order 5(4)
DOP853	Explicit Runge-Kutta method of order 8
Radau	Implicit Runge-Kutta method of the Radau IIA family of order 5
BDF	Implicit multi-step variable-order (1 to 5) method based on a backward differentiation formula for the derivative approximation
LSODA	Adams/BDF method with automatic stiffness detection and switching

See Scipy solver_ivp() for more info. Recommended solvers: BDF, LSODA, or Radau for stiff problems or RK45 for non-stiff problems with smooth blending.

ODE_tolerancefloat, optional

The relative and absolute tolerance for the ODE solver.

Returns:

statesdictionary of arrays

A dictionary of Numpy arrays representing the properties of the fluid at each evaluation point.

solutionscipy.integrate.OdeResult

The result of the ODE integration containing information about the solver process.