Getting Started

Thermofluids makes modeling thermofluid systems easier. The essential element is the statepoint. statepoint objects work like physical statepoints: they have two independent variables, many dependent variables, and sometimes extensive variables. You can use them like this:

sp = Statepoint(fluid='water',  # define a fluid
                t=25.0,         # define one independent variable
                p=1.0,          # define a second independent variable
                m=0.1)          # define extensive variable

print(sp.rho, sp.h, sp.visc)


>>> (<Quantity(997.047039, 'kilogram / meter ** 3')>, <Quantity(104.918893, 'kilojoule / kilogram')>, <Quantity(0.000890022674, 'kilogram / meter / second')>)

Representing statepoints this way makes it easy to model thermodynamic cycles. The module includes some common thermodynamic unit processes to manipulate statepoints in physically meaningful ways. For example, we can model a refrigeration cycle:

from thermofluids import Statepoint, Quantity
from thermofluids.ProcessEquipment import Compressor


# Select a refrigerant
REFR = 'R134A'

# Find the saturation pressures of the refrigerant
ssp = Statepoint(fluid=REFR, t=-10, x=1.).p
scp = Statepoint(fluid=REFR, t=25., x=1.).p

# Define a statepoint representing the state between the evaporator and compressor
evap_to_comp = Statepoint(fluid=REFR, t=-5., p=ssp, m=0.1)

# Use the built in compressor object to calculate the statepoint at the compressor outlet
comp = Compressor(inlet=evap_to_comp,
                  eta=0.7,
                  Pout=scp)
comp_to_cond = comp.get_outlet()
print(f'Compressor work: {comp.work.to("kW")}')

# Derive the statepoint representing the condensed state between the condenser and expansion valve
cond_to_exp = Statepoint(fluid=REFR, t=20., p=scp, m=evap_to_comp.m)
q_cond = cond_to_exp.m * (comp_to_cond.h - cond_to_exp.h)
print(f'Condenser heat: {q_cond.to("kW")}')

# Derive the statepoint representing the isenthalpically expanded refrigerant
exp_to_evap = Statepoint(fluid=REFR, p=ssp, h=cond_to_exp.h, m=cond_to_exp.m)
q_evap = evap_to_comp.m * (evap_to_comp.h - exp_to_evap.h)
print(f'Evaporator heat: {q_evap.to("kW")}')
cop = q_evap / comp.work
print(f'Coefficient of performance: {cop.to("")}')


>>> Compressor work: 3.637 kW
>>> Condenser heat: 20.58 kW
>>> Evaporator heat: 16.95 kW
>>> Coefficient of performance: 4.659

Installation

Thermofluids is installable using standard python tools such as pip install python. See https://pypi.org/project/thermofluids/

Acknowledgement of CoolProp and Pint

Statepoints are, at a basic level, a nice looking wrapper around CoolProp property lookups and Pint Quantity for unit management. We’re grateful to those projects!

Code Docs

Essential Functionality

class thermofluids.Base.Statepoint.Statepoint

A statepoint instance represents a thermofluid statepoint. It fully characterizes a fluid at a point. Specify and two independent intensive properties and the rest will be automatically calculated. Some extensive properties are also logged, such as mass flow rate. This class uses CoolProp for the backend property lookups.

Do not change the properties of a statepoint once it has been fully defined- create a new instance instead.

Parameters:

  • name – Name of the statepoint for your own reference. Any string

  • fluid – Any CoolProp fluid. Can be a simple option like “water” or “R404A”, an incompressible fluid like “INCOMP::MMA”, or a Refprop fluid, if you have access to that library, like “REFPROP-MIX:METHANOL&WATER”

  • t – Temperature in Celcius

  • p – Pressure in bar

  • x – Quality in fraction 0-1. All liquid is 0, all vapor is 1.

  • h – Enthalpy in kJ/kg. Reference points are those from CoolProp/Refprop

  • s – Entropy kJ/kg-C

  • c – Concentration (used for INCOMP and REFPROP-MIX mixtures) kg/kg

  • m – Mass flow rate kg/s

  • twb – Wet bulb temperature C

  • tdp – Dew point temperature C

  • rh – Relative humidity -

  • omega – Humidity ratio kg water/kg dry air

  • source – Another statepoint instance to import properties from

  • source_vars – A tuple of strings of the properties you’d like to import from source

All the inputs are also available as attributes. Additional dependent attributes include:

  • cp: heat capacity at constant pressure kJ/kg-C

  • cv: heat capacity at constant volume kJ/kg-C

  • rho: density kg/m^3

  • u: internal energy kJ/kg

  • gamma: heat capacity ratio

  • k: thermal conductivity kW/m-C

  • tlsat: saturation temperature of liquid aka bubble point C (Refprop only)

  • tvsat: saturation temperature of vapor aka dew point C (Refprop only)

  • cvap: equilibrium concentration of vapor for a liquid kg/kg (Refprop only)

  • ss: speed of sound m/s

  • visc: Dynamic viscosity (mu) kg/m-s^2

  • pr: Prandtl number -

  • sigma: surface tension N/m

  • beta: coefficient of thermal expansion C^-1

  • t_freeze: freezing temperature C

Statepoints can be added or divided (see methods below).

Examples:

  • sp = Statepoint(fluid='water', t=25., p=1.)     # water at room conditions

  • sp = Statepoint(fluid='INCOMP::MMA', t=-20., p=20., c=0.25, m=0.2)  # flow of methanol and water

  • sp = Statepoint(fluid='REFPROP-MIX:METHANOL&WATER', x=0.8, p=0.25, c=0.95, m=1.5)

  • sp = Statepoint(source=sp, source_vars=('fluid', 'm', 't'), x=0)    # import everything but quality from another statepoint

Warnings:

You can set Statepoint.allow_warnings = False to suppress warnings

__add__(other)

Adding two statepoints mixes them isenthalpically, conserving mass and energy. It assumes the higher pressure fluid is expanded to the lower pressure The function returns another statepoint instance with fluid, mass, enthalpy, and pressure (and optionally concentration) set

__div__(other)

Diving a statepoint by an integer divides the mass flow rate by the integer leaving all other properties equal :param other: an integer :return: new statepoint

Fluids

The Fluids submodule supplies functions to perform common fluids calculations including non-dimensional characterization and pressure drops in various situations.

thermofluids.Fluids.Drag.f_drag(shape, density, velocity, area=None)
Parameters:

  • shape – Shape of the object <str>

  • density – Density of the fluid <Quantity, kg/m^3>

  • velocity – Velocity of the fluid relative to the thing getting dragged <Quantity, m/s>

  • area – Optional, frontal area of the object <Quantity, m^2>

Returns:

Pressure on the object (without area specified) or force on the object (with area specified)

thermofluids.Fluids.Drag.settling_vel(statepoint, d_solid, rho_solid, d_h_ext, re='auto', detailed=False)

This function calculates the settling velocity for a solid particle in a flowing liquid. If the fluid is moving slower than the settling velocity, the particle will tend to settle out. If its moving faster, it will tend to remain entrained. There are three regimes:

  1. Stokes flow (Re<10: force balance is Stokes drag vs gravity

  2. Laminar flow (10<Re<3000): force balance is hydrodynamic drag vs gravity

  3. Turbulent flow (Re>3000): force balance is hydrodynamic drag vs gravity

Parameters:

  • statepoint – statepoint instance of the fluid containing the solids

  • d_solid – diameter of the solid particle

  • rho_solid – density of the solid particle

  • d_h_ext – hydraulic diameter of the flow area

  • re – Reynolds number, if known. If left default=’auto’, it will be calculated from the settling velocity

Returns:

settling velocity [m/s]

From http://nullege.com/codes/show/src@l@i@libtsl-HEAD@mr_darcy.py Edited by Luke Erickson to include transition region

thermofluids.Fluids.Moody.frictionfactor(Re, relrough)

Returns the Darcy friction factor based on reynolds # and relative roughness.

thermofluids.Fluids.Re.Re_ext(statepoint, vel, d1, d2=None, shape=None)

Reynolds number for external flow :param statepoint: statepoint of fluid flowing over object :param vel: velocty m/s :param d1: dimension 1 [m] :param d2: dimension 2 [m] :param shape: not implemented :return:

thermofluids.Fluids.PressureDrop.pressure_drop(statepoint, roughness, L, d1, d2=None, v=None, shape='round', props='REFPROP', units='bar', detailed=False)
Parameters:

  • statepoint

  • roughness

  • L

  • d1

  • d2

  • shape

  • props

Returns:

thermofluids.Fluids.PressureDrop.pressure_drop_corrugated_plates(statepoint, length, width, spacing, chevron_angle, detailed=False, units='bar')

Calcualtes the pressure drop through a channel of corrugated plates, such as in a brazed plate heat exchanger Based on http://achp.sourceforge.net/ACHPComponents/PlateHeatExchanger.html#single-phase-flow Seems to match Alfa Lavals pressure drop calculations pretty well :param statepoint: a statepoint instance of the fluid flowing through the hx :param length: plate length (in direction of flow) :param width: plate width (perpendicular to direction of flow) :param spacing: spacing between the plates [m] :param chevron_angle: Angle between flow direction and chevron lines in radians. pi/2 is perpendicular to flow, 0 is parallel to flow :return: pressure drop in specified units

thermofluids.Fluids.PressureDrop.pressure_drop_fitting(statepoint, fitting_type, d1, roughness, units='bar')

Calculates the pressure drop through an elbow Based on https://neutrium.net/fluid_flow/pressure-loss-from-fittings-equivalent-length-method/ Secondary source for L/d values: http://www.metropumps.com/ResourcesFrictionLossData.pdf :param statepoint: statepoint instance of fluid flowing through pipe :param fitting_type: ‘elbow’, ‘long elbow’, ‘45 elbow’ :param d1: Hydraulic diameter [m] :param roughness: pipe roughness [m] :param units: return units [str] :return: pressure drop [return units]

Heat Transfer

The Heat Transfer module supplies functions to perform common heat transfer calculations including non-dimensional characterization, effective temperature differences, two-phase flow regime prediction, and heat transfer coefficients.

thermofluids.HeatTransfer.lmtd.lmtd(t_hot_in, t_hot_out, t_cold_in, t_cold_out)

Log Mean Temperature Difference :param t_hot_in: C :param t_hot_out: C :param t_cold_in: C :param t_cold_out: C :return: log mean temp difference C

thermofluids.HeatTransfer.Nusselt.Nu(Re, Pr, f=None, mode=None)

Nusselt number for forced convection :param Re: Reynolds number :param Pr: Prandtl number :param f: Darcy friction factor :param mode: ‘heating’ or ‘cooling’ :return: Nusselt number

thermofluids.HeatTransfer.Nusselt.Nu_natural(Ra, Pr, shape)

Nusselt number for natural convection References: [1] Cengel Heat and Mass Transfer 3rd Ed p511 :param Ra: Rayleigh number :param Pr: Prandtl number :param shape: can be ‘vertical plate’, ‘inclined plate’, ‘horizontal plate’, ‘vertical cylinder’, ‘horizontal cylinder’, ‘sphere’ :return: Nusselt number <Quantity, ->

thermofluids.HeatTransfer.Ra.Ra(statepoint, gr=None, L=None, t_wall=None)

Rayleigh number for natural convection Please specify either the Grashof number or the characteristic length and wall temperature :param statepoint: A statepoint instance of the fluid being heated/cooled :param gr: Grashof number <Quantity [-]> :param L: Characteristic length <Quantity [m]> :param t_wall: Wall temperature <Quantity [C]> :return:

This file contains functions to calculate the heat transfer for convective boiling flow in a tube. It is based on the Shah correlation as described in Carey ‘Liquid-Vapor Phase-Change Phenomena’ p616.

thermofluids.HeatTransfer.TwoPhase.Bo(q_flux, G, h_v)

Boiling number - the ratio of flow rate due to boiling vs flow rate through pipe :param q_flux: heat flux [kW/m^2] :param G: mass flux [kg/m^2-s] :param h_v: enthalpy of vaporization [kJ/kg] :return: boiling number [-]

thermofluids.HeatTransfer.TwoPhase.Co(x, rho_v, rho_l)

Convection number - Some representation of volumetric liquid/gas ratio (higher is more liquid) :param x: quality [kg/kg] :param rho_v: vapor density [kg/m^3] :param rho_l: liquid density [kg/m^3] :return: Convenction number [-]

thermofluids.HeatTransfer.TwoPhase.Fr_le(G, rho_l, d)

Froude number - the ratio of flow inertia to gravity field :param G: mass flux [kg/m^2-s] :param rho_l: liquid density [kg/m^3] :param d: hydraulic diameter [m] :return: Froude number [-]

thermofluids.HeatTransfer.TwoPhase.crit_heat_flux(h_v, rho_l, rho_v, sigma)

Critical heat flux for boiling heat transfer according to: https://en.wikipedia.org/wiki/Critical_heat_flux#Correlations_for_critical_heat_flux :param h_v: enthalpy of vaporization [kJ/kg] :param rho_l: liquid density [kg/m^3] :param rho_v: vapor density [kg/m^3] :param sigma: interfacial tension [N/m] :return: critical heat transfer [kW/m^2]

thermofluids.HeatTransfer.TwoPhase.heat_flux_nucleate_boiling(mu_l, h_v, rho_l, rho_v, sigma, cp_l, T_s, T_sat, Csf, n, Pr_l)

heat flux for nucleate boiling based on Rohsenow correlation https://en.wikipedia.org/wiki/Nucleate_boiling#Boiling_Heat_transfer_correlations :param mu_l: Liquid viscosity [kg/m-s] :param h_v: enthalpy of vaporization [kJ/kg] :param rho_l: liquid density [kg/m^3] :param rho_v: vapor density [kg/m^3] :param sigma: interfacial tension [m^2/s] :param cp_l: liquid specific heat capacity [kJ/kg-C] :param T_s: surface temperature [C] :param T_sat: saturation temperature [C] :param Csf: surface fluid combination :param Pr_l: liquid prandtl number :param n: parameter depending on surface fluid combo, typically 1.0 or 1.7 :return: heat flux [kW/m^2]

thermofluids.HeatTransfer.TwoPhase.htc_conv_boil(x, rho_l, rho_v, G, h_v, d, q_flux, mu_l, k_l, pr_l)

Heat transfer coefficient of boiling flow through a tube based on Shah correlation as described in Carey ‘Liquid-Vapor Phase-Change Phenomena’ p616. :param x: quality [-] :param rho_l: liquid density [kg/m^3] :param rho_v: vapor density [kg/m^3] :param G: mass flux [kg/m^2-s] :param h_v: enthalpy of vaporization [kJ/kg] :param d: hydraulic diameter [m] :param q_flux: Boiling heat flux [kW/m^2] :param mu_l: liquid kinematic viscosity [kg/m-s] :param k_l: liquid thermal conductivity [kW/m-C] :param pr_l: liquid Prandtl number [-] :return: Heat transfer coefficient [kW/m^2-C]

Special Properties

Special properties are a handful of properties we’ve needed over time that aren’t available in CoolProp.

thermofluids.Properties.BrineProps.get_P_liq_sat(fluid, c, t)

Calculates the pressure at which a fluid is a saturated liquid, aka bubble point :param fluid: fluid name- must be a legal CoolProp->Refprop input :param c: concentration [wt] :param t: temperature [C] :return: bubble point pressure [bar]

thermofluids.Properties.BrineProps.get_P_vap_sat(fluid, c, t)

Calculates the pressure at which a fluid is a saturated liquid, aka dew point :param fluid: fluid name- must be a legal CoolProp->Refprop input :param c: concentration [wt] :param t: temperature [C] :return: dew point pressure [bar]

thermofluids.Properties.BrineProps.get_T_liq_sat(fluid, c, p)

Calculates the temperature at which a fluid mixture is a saturated liquid, aka bubble point :param fluid: fluid name- must be a legal CoolProp->Refprop input :param c: concentration [wt] :param p: pressure [bar] :return: bubble point temperature [C]

thermofluids.Properties.BrineProps.get_T_vap_sat(fluid, c, p)

Calculates the temperature at which a fluid mixture is a saturated vapor, aka dew point :param fluid: fluid name- must be a legal CoolProp->Refprop input :param c: concentration [wt] :param p: pressure [bar] :return: dew point temperature [C]

thermofluids.Properties.BrineProps.get_c_freeze(fluid, T)

calculates the the concentration at which a brine will freeze :param fluid: [str] fluid name :param T: [float] temperature in C :return: [float] concentration

Process Equipment

The Process Equipment supplies classes and functions that represent unit functions that take and produce thermodynamic statepoints. The performance/outputs are mostly theoretically derived as opposed to empirically derived from specific implementations of these types of process equipment.

class thermofluids.ProcessEquipment.Column.Column(P_col=None, inlet=None, outlet_vap=None, outlet_liq_c=None, fluid=None)

This class models distillation processes. There are three performance models:

  1. ME_balance - This model just does mass and energy balances to determine the inputs and outputs. Use this model when you need high speed calculations and you don’t need any info about column design

  2. mccabe - This model does tray-wise mass and (approx.) energy balances using the McCabe Thiele model Use it to get more detailed information about the internal performance of the column with minimal computation time.

  3. mccabe_diabatic - This model does tray-wise mass and (approx.) energy balances using the McCabe Thiele model. The operational line for this model is parallel to the VLE line to model a diabatic process. It solves the model from the top tray downward such that inlet feed and outlet vapor are defined. Use this model as a detailed, functional analysis of a diabatic column.

  4. mccabe_diabatic_bottomup - This model is the same as mccabe_diabatic except that it attempts to solve it from the bottom up such that inlet feed and outlet liquid are defined. DON’T use this model… this function is currently buggy (go with mccabe_diabatic instead)

  5. FullDetail - This model does a complete mass energy balance on every tray. Use it when you need the most internal performance detail, require a column design, and have some time to let it run. This function is currently buggy (go with mccabe_diabatic instead)

ME_balance()

This is the simplest column model. It just does mass and energy balances to determine how much heat is required to get the desired top and bottom products from the inlet stream.

You can use the following results from this model: Column.bottoms Column.Q_input

get_diabatic_curves(shift, plotfile=None, printmode=True)

This calculates the coordinates for the vapor-liquid equilibrium line and operating line for a given column. The operating line is that which minimized exergy destruction. In order to find this operating line, this function finds the shift value that will move the VLE line at a 45deg angle southeast such that it satisfies the inlet liquid and outlet vapor concentrations. That shift value is used in diabatic_mccabe during the McCabe-Thiele calculations. (See “Diabatic Column Optimization Compared to Isoforce Columns” by Sauar et al)

You can use the following results from this function:

column.x_sat <– saturated liquid concentrations for temperature range column.y_sat <– saturated vapor concentrations for temperature range column.x_oper <– operating liquid concentrations for temperature range column.y_oper <– operating vapor concentrations for temperature range column.oper_shift <– value used to make oper. line from VLE line (x_VLE + oper_shift, y_VLE - oper_shift)

Written by Fiona Majeau for Rebound Technology July 1 2016

mccabe(plotfile=None)

This is a graphical way to do an approximate mass and energy balance on each tray For more info: https://en.wikipedia.org/wiki/McCabe%E2%80%93Thiele_method

You can use the following results from this model: Column.n_trays Produces a png picture file with the McCabe Theile diagram

mccabe_diabatic(shift=None, plotfile=None, printmode=True)

This is a graphical way to do an approximate mass and energy balance on each tray for a diabatic column. Creates McCabe-Thiele diagram using an minimum exergy destruction operational line that sets the liquid outlet state (bottoms) given Statepoints for the inlet liquid and outlet vapor. Using COE, the necessary Q_sens is calculated for each tray. (See “Diabatic Column Optimization Compared to Isoforce Columns” by Sauar et al)

You can use the following results from this model:

Column.num_trays <– number of trays Column.Q_feed_to_satliq <– Q to heat inlet feed to saturated liquid (if subcooled) Column.Q_sens <– array of Q_sens to each tray (1st val is top tray) Column.Q_trays_total <– total Q going to trays (sum of Q_sens) Column.Q_total <– total Q going to column (Q_trays_total + Q_feed_to_satliq) Column.tray_temps <– array of temperature for each tray (1st val is top tray) Column.tray_concs <– array of concentration arrays for each tray (1st val is top tray) Column.tray_masses <– array of mass flow arrays for each tray (1st val is top tray) Column.bottoms <– Statepoint of liquid outlet

As well as the outputs from get_diabatic_curves (called within)

column.x_sat <– saturated liquid concentrations for temperature range column.y_sat <– saturated vapor concentrations for temperature range column.x_oper <– operating liquid concentrations for temperature range column.y_oper <– operating vapor concentrations for temperature range column.oper_shift <– value used to make oper. line from VLE line (x_VLE + oper_shift, y_VLE - oper_shift)

Produces a png picture file with the McCabe-Theile diagram (includes VLE line and operating line)

Written by Fiona Majeau for Rebound Technology June 30 2016 Rectifying section added by Luke Erickson for Rebound Technology Nov 30 2016

off_design(heating_stream, n_trays, h_tray, d_col, htc, plotfile=None)

This model calculates the performance given both an inlet stream and a heating stream. Most of the other models are set up to do design, meaning they output column specifications or heat requirements This one, in contrast, just outputs performance based on a specified design parameters and inputs.

To do:

-It currently is set up for only a square column heating from the perimeter. We should make other options available -Improve the solver speed. I think the best ways to do this would be to make the conc/mass convergence criteria adaptive and the heat transfer step size adaptive. A bigger but maybe better change would be to do something more complex than damped direct substitution. -LMTD for the top tray is currently not accounting for the fact that the Q/T curve is significantly non-linear because that section has both sensible and latent heating. This should be fixed. -Set up plotting again. -Make outputs more useful.

Parameters:

  • heating_stream – A statepoint instance representing the heating fluid coming into the bottom of the column

  • n_trays – The number of trays [-]

  • h_tray – The height of a tray [m]

  • d_col – Width of the column [m]

  • htc – Heat transfer coefficient [kW/m^2-C]. I typically assume this is equal the HTC of the heating fluid.

  • plotfile – [string] name of png file to output a mccabe thiele style diagram. Plotting is not fully enabled right now

Returns:

Currently just prints all the outputs

class thermofluids.ProcessEquipment.Compressor.Compressor(inlet=None, inject=None, eta=None, Pout=None, output_name=None, work=None)

Models a Compressor or Pump

Parameters:

  • inlet – [Statepoint] compressor inlet statepoint

  • eta – [Quantity] or [float] icentropic efficiency of turbo equipment

  • Pout – [Quantity] or [float] pressure at outlet of compressor or pump

Returns:

[Statepoint] compressor outlet statepoint

class thermofluids.ProcessEquipment.Turbine.Turbine(inlet=None, eta=None, Pout=None, output_name=None, work=None)

Models a turbine

Parameters:

  • inlet – [Statepoint] turbine inlet stream

  • eta – [Quantity] or [float] icentropic efficiency of turbo equipment

  • Pout – [Quantity] or [float] pressure at outlet of turbine

Returns:

[Statepoint] turbine outlet stream

Manipulations of definition of flow coefficient, Cv=Q*sqrt(SG/dP), to calculate valve performance. https://en.wikipedia.org/wiki/Flow_coefficient

thermofluids.ProcessEquipment.Valve.get_valve_C_v(statepoint, dP)

Calculate the flow coefficient of a valve from flow rate and pressure drop

Parameters:

  • statepoint

  • dP

Returns:

thermofluids.ProcessEquipment.Valve.get_valve_dP(statepoint, C_v, units='bar')

Calculate pressure drop across a valve from flow coefficient and flow rate

Parameters:

  • statepoint

  • C_v

  • units

Returns:

thermofluids.ProcessEquipment.Valve.get_valve_flow_rate(statepoint, dP, C_v, units='kg*s^-1', debug=False)

Calculate flow rate through a valve from pressure drop and flow coefficient

Parameters:

  • statepoint

  • dP – Pressure drop across the valve

  • C_v – Flow coefficient

  • units

  • debug

Returns: