LMTD calculation in heat exchanger

LMTD stands for Logarithmic Mean Temperature Difference, which is a method of calculating the arithmetic mean temperature difference in a heat exchanger. The LMTD is calculated using the following formula:

Formula of LMTD

LMTD = (ΔT1-ΔT2) / ln (ΔT1/ΔT2)

Here, ΔT1 is the temperature difference between the hot and cold fluids on one side of the heat exchanger, and ΔT2 is the temperature difference on the other side. This formula is used to calculate the heat transfer surface area of a heat exchanger.

There are two types of heat exchangers such as parallel-flow and counter-flow heat exchangers that differ in the direction of fluid flow. In a parallel-flow heat exchanger, the hot and cold fluids enter at the same end, flow in the same direction, and leave at the same end. In contrast, in a counter-flow heat exchanger, the fluids enter at opposite ends, flow in opposite directions, and leave at opposite ends.

The counter-flow heat exchanger is more efficient than the parallel-flow heat exchanger because it allows for a greater temperature difference between the hot and cold fluids. This results in a higher rate of heat transfer and a smaller heat exchanger size.

For parallel flow in tubular heat exchangers

LMTD = ((Thi-Tci)-(Tho-Tco))/ln((Thi-Tci)/(Tho-Tco))

For counterflow in tubular heat exchangers
LMTD = ((Tho-Tci)-(Thi-Tco))/ln((Tho-Tci)/(Thi-Tco)) 

Where,
LMTD = log mean temperature difference (degC or degF)
Thi = inlet temperature of the hot fluid (degC or degF)
Tho = outlet temperature of the hot fluid (degC or degF)
Tci = inlet temperature of the cold fluid (degC or degF)
Tco = outlet temperature of the cold fluid (degC or degF)

Python code

LMTD calculation in the heat exchanger is coded in Python as follows.

import math
def lmtd(Thi, Tho, Tci, Tco):
    if min(Thi, Tho) > max(Tci, Tco):
        dt = (Thi-Tci) if (Thi-Tci) == (Tho-Tco) else ((Thi-Tci)-(Tho-Tco))/math.log((Thi-Tci)/(Tho-Tco))
    else:
        dt = (Thi-Tci) if (Tho-Tci) == (Thi-Tco) else ((Tho-Tci)-(Thi-Tco))/math.log((Tho-Tci)/(Thi-Tco))
    return dt
print("parallel flow case = ", lmtd(150, 120, 50, 80))
print("counter flow case = ", lmtd(180, 130, 100, 120))

Here, Thi, Tho, Tci and Tco are the temperatures of the hot and cold fluids at the inlet and outlet of the heat exchanger. The function lmtd calculates the LMTD using the formula (ΔT1-ΔT2)/ln(ΔT1/ΔT2) and returns the result.

parallel flow case =  65.48
counter flow case =  33.66

LHV calculation of mixed fuel gas (ASME PTC-22)

ASME Performance Test Codes (PTC 22) establishes directions and rules for the conduct and results reporting of thermal performance tests for open cycle gas turbine power plants and gas turbine engines. This performance test code provides explicit instruction on determining corrected power, heat rate (efficiency), exhaust flow, exhaust energy, and exhaust temperature. 

Guidance is also provided for designing testing requirements and programs to satisfy different goals such as absolute performance and comparative performance.

Introducing the fuel heating value table introduced here.
The procedure for calculating the mixed fuel gas LHV is written in Python as follows.

Run Python code below : https://www.mycompiler.io/view/41iScgW4CFz

dbHC = [# Formula, Molecular Weight [lb/lbmol], Standard Density [lb/1000 ft3], High Heating Value [Btu/lbm], Low Heating Value [Btu/lbm]
       ["CH4", 16.0425, 42.274, 23892.2, 21511.9],
       ["C2H6", 30.0690, 79.237, 22334.1, 20429.2],
       ["C3H8", 44.0956, 116.199, 21654.1, 19922.2],
       ["C4H10", 58.1222, 153.161, 21232.3, 19589.8],
       ["C4H10", 58.1222, 153.161, 21300.2, 19657.8],
       ["C5H12", 72.1488, 190.123, 21043.7, 19455.9],
       ["C5H12", 72.1488, 190.123, 21085, 19497.2],
       ["C6H14", 86.1754, 227.085, 20943.8, 19392.9],
       ["C7H16", 100.2019, 264.048, 20839.1, 19314.7],
       ["C8H18", 114.2285, 301.010, 20759.7, 19255.4],
       ["C9H20", 128.2551, 337.972, 20701, 19212.3],
       ["C10H22", 142.2817, 374.934, 20651.6, 19175.5],
       ["CO", 28.0101, 73.811, 4342.2, 4342.2],
       ["CO2", 44.0095, 115.972, 0.0, 0.0],
       ["H2S", 34.0809, 89.808, 7094.1, 6533.8],
       ["Air", 28.9651, 76.328, 0.0, 0.0],
       ["H2", 2.0159, 5.312, 61022.3, 51566.7],
       ["O2", 31.9988, 84.322, 0.0, 0.0],
       ["N2", 28.0134, 73.820, 0.0, 0.0],
       ["H2O", 18.0153, 47.473, 1059.8, 0.0],
       ["He", 4.0026, 10.547, 0.0, 0.0],
       ["Ar", 39.9480, 105.269, 0.0, 0.0],
       ]

def hclhv(Comp):

    n = len(dbHC)
    sumComp = sum(Comp)
    for i in range(n):
        Comp[i] = Comp[i] / sumComp

    sumLHV = 0
    for i in range(n):
        LHV = dbHC[i][4] * dbHC[i][1] / ((10.7316*519.67) / 14.696)
        sumLHV = sumLHV + Comp[i] * LHV
    sumComp = sum(Comp)

    return sumLHV / sumComp

Comp = [90, 4, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0]
print("Mixed Fuel Gas LHV [Btu/Scf] = ", hclhv(Comp))

When run the code, you get the results below.

Mixed Fuel Gas LHV [Btu/Scf] =  959.622

DOE motor challenge program

Motor Challenge is a partnership program between the U.S. Department of Energy and the nation’s industries. The program is committed to increasing the use of industrial energy-efficient electric motor systems and related technologies.

This program is wholly funded by the U.S. Department of Energy and is dedicated to helping industry increase its competitive edge, while conserving the nation’s energy resources and enhancing environmental quality.

Over half of all electrical energy consumed in the United States is used by electric motors. Improving the efficiency of electric motors and the equipment they drive can save energy, reduce operating costs, and improve our nation’s productivity.

Energy efficiency should be a major consideration when you purchase or rewind a motor. The annual energy cost of running a motor is usually many times greater than its initial purchase price. For example, even at the relatively low energy rate of $0.04/kWh, a typical 20-horsepower (hp) continuously running motor uses almost $6,000 worth of electricity annually, about six times its initial purchase price.

Select a new energy-efficient motor under any of the following conditions:

• The motor is less than 40 hp.

• An energy-efficient motor is recommended according to Table 3.

• The cost of the rewind exceeds 65% of the price of a new motor.

• The motor was rewound before 1980 year

Survey your motors. Gather nameplate information and obtain field measurements (voltage, amperage, power factor, operating speed) under typical operating conditions. Initially focus on motors that exceed minimum size and operating duration criteria. Typical selection criteria include:

• Three-phase NEMA design B motor

• Non-specialty motor

• 10 to 600 hp

• At least 2000 hours per year of operation

• Constant load (not intermittent, cyclic, or fluctuating

• Older or rewound standard efficiency motors

• Easy access

• Readable nameplate.

BHP = (ρ*(Q/3600)*H)/(102*η)

BHP : Pump Power [Kw]
ρ  : Fluid Density [kg/m3]
Q   : Flow Rate [m3/hr]
H   : Total Diff' Head [m]
η  : Pump Efficiency [%]

Design Condition
ρ = 1,126 kg/m3 
Q = 173.0 m3/hr
P1 = 0.6 kgf/cm2g
P2 = 14.4 kgf/cm2g
H = (14.4-0.6)*10/(1126/1000) = 122.6 m
η = 66.5%
BHP = (1126*(173/3600)*122.6)/(102*66.5/100) = 97.8 kW

Actual Condition
ρ = 1,126 kg/m3
Q = 84.4 m3/hr
P1 = 0.6 kgf/cm2g
P2 = 14.6 kgf/cm2g
H = (14.6-0.6)*10/(1126/1000) = 124.3 m
n = 53.0%
BHP = (1126*(84.4/3600)*124.3)/(102*53.0/100) = 60.7 kW

Load Factor = Pump Operation BHP / Pump Rated BHP
            = 60.7 kW /97.8 kW = 62%
Pump Imbalance = [(Actual Q * Actual H) / (Design Q * Design H) -1] * 100%
               = ((84.4*124.3)/(173*122.6)-1)) = -51%
Flow Ratio = Actual Q / Design Q

           = 84.4 m3/hr / 173 m3/hr = 49%
Motor Load = Operation P / Rated P

           = 63.8 kW / 110 kW = 58%

If the calculation results of the four parameters calculated above meet the conditions below, energy engineer should consider an investment project to save energy consumption.

Pump Imbalance: Review below -20%
Pressure Ratio: Review when exceeding 130%
Flow Ratio: Review when below 70%
Motor Load: Review when below 40%

Crane differential pressure or flow rate for compressible fluid (Imperial)

This Python code calculates differential pressure or flow rate according to given pipe & fitting information for compressible fluid.

Formula of DP and flow rate

Calculate differential pressure in pipes based on information on the following pipes/accessory equipment

w    : Flow rate [lb/hr]
△P  : Delta Pressure [psia]
L    : Length [ft]
D    : Pipe Diameter [ft]
d    : Pipe Diameter [in]
ρ    : Density [lb/ft3]
f    : Friction factor

Here, Friction factor (f) must be calculated using separate charts and programs.

Calculation formulas based on CRANE BOOK

Kpipe = f * l/D
K     = Kpipe + Kfitting

compressible fluid
△P = (K * W^2) / (1891^2 * d^4 * ρ  * Y^2)
w = 1891 * d^2 * √(△P * ρ / K) * Y

liquid fluid
△P = (K * W^2) / (1891^2 * d^4 * ρ)
w = 1891 * d^2 * √(△P * ρ / K)

Python code

Run Python code below : https://www.mycompiler.io/view/6fwILhy1i1P

import math

def compressibledp(W, l, d, r, Y, f, Kf):
    D = d * 12
    Kp = f*l/D
    K = Kp + Kf
    dp = (K*pow(W, 2))/(pow(1891, 2)*pow(d, 4)*r*pow(Y, 2))
    return dp

W = 220      # Flow rate (W, lb/hr)
l = 30.48    # Length (L, ft)
d = 1.969    # Pipe Diameter (d, in)
r = 0.0774   # Density (ρ, lb/ft3)
Y = 1.0      # Net Expansion Factor (Y)
f = 0.005    # Friction factor (f)
Kf = 500     # Resistance coefficient of fittings

dp = compressibledp(W, l, d, r, Y, f, Kf)
print("delta pressure of pipe = ", dp, "psia")

def compressibleflow(dp, l, d, r, Y, f, Kf):
    D = d * 12
    Kp = f * l / D
    K = Kp + Kf
    W = 1891 * pow(d, 2) * math.sqrt(dp * r / K) * Y
    return W

dp = 5.933   # Delta Pressure (△P, psia)
m = 30.48    # Length (L, ft)
d = 1.969    # Pipe Diameter (d, in)
r = 0.0774   # Density (ρ, lb/ft3)
Y = 1        # Net Expansion Factor (Y)
f = 0.005    # Friction factor (f)
Kf = 500     # Resistance coefficient of fittings

flow = compressibleflow(dp, l, d, r, Y, f, Kf)
print("mass flow of pipe = ", flow, "lb/hr")

When run the code, you will receive the following results.

delta pressure of pipe =  5.82 psia
mass flow of pipe =  222 lb/hr

Gas mean pressure in long pipe

Equation for mean pressure for steady flow in a gas pipe at pressures below 1000 psia and temperatures above 60°F, the change in compressibility factor z with Pressure is approximately linear. Therefore, the compressibility factor z can be expressed as a constant, so the average pressure of the gas fluid can be calculated as follows.

Pm = 2/3 * (P1^3 - P2^3)/(P1^2 - P2^2)

def pmean(p1, p2):

    return 2/3*(pow(p1, 3) - pow(p2, 3))/( pow(p1, 2) - pow(p2, 2))

print("Gas mean pressure of 120 psia and 100 psia = ", pmean(120, 100))

When run the code, you will receive the following results.

Gas mean pressure of 120 psia and 100 psia =  111.429

How to Draw Heat & Cool Compositive Curve in HYSYS

How to Draw Heat & Cool Compositive Curve in HYSYS

1. Complete the LNG exchange Object.

2. Select the [Tools]-[Utilities] menu.

3. Add the Compositive Curve Utility menu in the right window.

4. Confirm that Compositive Curve Utility is added to the left window and double-click it.
    Then select the LNG Excange that was previously completed.

5. Check the Pinch Temperature in the Pinch Result tab.

6. Check Compositive Curve in the Plot tab.

  ** If necessary, copy it to the clipboard and add it to the report.

Heat loss method for heater efficiency

The heat loss method disused in ASME PTC 4.1 may be used for evaluating the efficiency of steam generators. For quick estimates of oil and natural gas fired steam generators, however, the following equations may be used. These equations were arrived at by the author after performing several calculations.

1. Natural gas 

Efficiency (%, HHV) = 89.4 - (0.001123 + 0.0195 * EA)*(Tg - Ta)
Efficiency (%, LHV) = 99.0 - (0.001244 + 0.0216 * EA)*(Tg - Ta)

Tg = exit gas temperatures, degF
Ta = reference air temperatures, degF
EA = excess air, EA = K * 21/(21-O2), K = 0.98 for natural gas 

2. Fuel Oil

Efficiency (%, HHV) = 92.9-(0.001298 + 0.0195 * EA) * (Tg-Ta) 
Efficiency (%, LHV) = 99.0-(0.001383 + 0.0203 * EA) * (Tg-Ta) 

Tg = exit gas temperatures, degF
Ta = reference air temperatures, degF
EA = excess air, EA = K * 21/(21-O2), K = 1.00 for fuel oil

Example-1:

A natural gas fired boiler with 15% excess air has an exit gas temperature of 380 degF, ambient = 90 degF. Determine efficiency on LHV basis. 

EA = 1.15
Efficiency (%, LHV) = 99.0 - (0.001244 + 0.0216 * EA)*(Tg - Ta)
Efficiency (%, LHV) = 99.0 - (0.001244 + 0.0216 * 1.15)*(380 - 90) = 91.4%

Example-2:

A natural gas fired boiler with 3% excess oxygen has an exit gas temperature of 350 degF, ambient = 85 degF. Determine efficiency on LHV basis. 

EA = K * 21/(21-O2) = 0.98 * 21/(21-3) = 1.143
Efficiency (%, LHV) = 99.0 - (0.001244 + 0.0216 * EA)*(Tg - Ta)

Efficiency (%, LHV) = 99.0 - (0.001244 + 0.0216 * 1.143)*(350 - 85) = 92.1%

Bernoulli Equation

Bernoulli Equation (SI unit)

Bernoulli's law is an important principle used in fluid mechanics, describing the relationship between the velocity and pressure of a fluid. According to this law, as speed increases, pressure decreases, and as speed decreases, pressure increases. This principle is expressed by the following formula:

Energy basis   : P1/ρ1 + u1^2/2 + z1 = P2/ρ2 + u2^2/2 + z2 [Pa, kg/m3, m/s, m]
Pressure basis : P1 + ρ1*u1^2/2 + ρ1*g*z1 = P1 + ρ2*u2^2/2 + ρ2*g*z2 [Pa, kg/m3, m/s, m]
Head basis     : P1/(ρ1*g) + u1^2/2g + z1 = P2/(ρ2*g) + u2^2/2g + z2 [Pa, kg/m3, m/s, m]
 
Where,
ρ = fluid density (kg/m3)
g = acceleration due to gravity (9.8 m/s)
P1 = pressure at elevation-1 (Pa)
u1 = velocity at elevation-1 (m/s)
h1 = height of elevation-1 (m)
P2 = pressure at elevation-2 (Pa)
u2 = velocity at elevation-2 (m/s)
h2 = height at elevation-2 (m)

Example of Bernoulli Equation (SI unit)

Water at 25 degC enters a pipe at the velocity of 1 m/s and a pressure of 101.3 kPa. The pipe has a constant diameter and friction loss is negligible. A change in the pipe’s elevation changes the downstream pressure in the pipe to 2.1 atm. Most nearly, what is the elevation change?

P1/(ρ1*g) + u1^2/2g + z1 = P2/(ρ2*g) + u2^2/2g + z2
u1 = u2, ρ1 = ρ2
P1/(ρ1*g) + z1 = P2/(ρ2*g) + z2
101.3 kPa = 101,300 Pa
2.1 atm = 212,730 Pa
z2 - z1 = (P1 - P2)/ρ*g = (101,300 Pa - 212,730 Pa)/(997 kg/m3 * 9.81 m/sec) = -11.4 m

Bernoulli Equation (Imperial unit)

Energy basis   : P1*gc/ρ1 + u1^2/2 + g*z1 = P2*gc/ρ2 + u2^2/2 + g*z2 [lbf/ft2, lbm/ft3, ft/s, ft, 32.2 lbm-ft/lbf-sec2, 32.2 ft/sec2]
Pressure basis : P1 + ρ1*u1^2/2gc + ρ1*g*z1/gc = P1 + ρ2*u2^2/2gc + ρ2*g*z2/gc [lbf/ft2, lbm/ft3, ft/s, ft, 32.2 lbm-ft/lbf-sec2, 32.2 ft/sec2]
Head basis     : P1*gc/(ρ1*g) + u1^2/2g + z1 = P2*gc/(ρ2*g) + u2^2/2g + z2 [lbf/ft2, lbm/ft3, ft/s, ft, 32.2 lbm-ft/lbf-sec2, 32.2 ft/sec2]

Where,
ρ = fluid density (lbm/ft3)
g = acceleration due to gravity (32.2 ft/sec2)
P1 = pressure at elevation-1 (lbf/ft2, = psia * 144 in2/ft2)
u1 = velocity at elevation-1 (ft/s)
h1 = height of elevation-1 (ft)
P2 = pressure at elevation-2 (lbf/ft2, = psia * 144 in2/ft2)
u2 = velocity at elevation-2 (ft/s)
h2 = height at elevation-2 (ft)

Example of Bernoulli Equation (Imperial unit)

The pump supplies water at a flow rate of 4.0 ft3/sec from an open reservoir through a horizontal 8 inch pipe. the head loss from the reservoir to the suction of the pump is 2 ft-lbf/lbm. and the discharge is at 90 psi in to 6 inch pipe (5.761 inch). The pump has an efficiency of 65%. The head that must be delivered by the pump to water is?

P1*gc/(ρ1*g) + u1^2/2g + z1 + hpump = P2*gc/(ρ2*g) + u2^2/2g + z2 + hf
ρ1 = ρ2 = 62.4 lbm/ft3
u1 = 0
z1 = 0
z2 = 0
hf = 0
P2 - P1 = 90 psi = 90 lbf/in2 = 12,960 lbf/ft2
6 inch pipe diameter = 3.14*(5.761 inch * (1 ft / 12 inch)/2)^2 = 0.181 ft2
u2 = 4.0 ft3/sec / 0.181 ft2 = 22.1 ft/sec
gc = 32.2 lbm-ft/lbf-sec2, 
g = 32.2 ft/sec2
hpump = (P2-P1)*gc/(ρ*g) + u2^2/2g = 12,960 lbf/ft2 / 62.4 lbm/ft3 + 22.1^2 / (2 * 32.2) = 215 ft

Darcy friction factor f calculation

Darcy friction factor f as a function of Reynolds number Re and pipe relative roughness ε / d, fitting the data of experimental studies of turbulent flow in smooth and rough pipes.

The equation can be used to (alliteratively) solve for the Darcy–Weisbach friction factor f.

For a conduit flowing completely full of fluid at Reynolds numbers greater than 4000, it is expressed as:

1/√(f) = -2*log((e/d)/3.7 + 2.51/(Nre*√(f))

As can see from the formula, simple calculation is not possible, so user should use iterative calculations or a solver to find the value of f that satisfies the formula.

Pipe Material Commercial Steel, e = 0.05
Pipe Internal Diameter (d, mm) = 25
Reynold's Number (Re) = 100000
Relative Roughness (ε/D)  0.002
Friction factor (f) 0.023

Run Python code below : https://www.mycompiler.io/view/4rbUNlvs1QN

import math
import numpy as np

def frictionfactor(piping, d, Nre):

    piping = "Commercial Steel"
    if piping == "Commercial Steel": e = 0.05
    if piping == "Cast Iron": e = 0.26
    if piping == "Galvanized Iron": e = 0.15
    if piping == "Asphalted Cast Iron": e = 0.12
    if piping == "Drawn Tubing": e = 0.0015

    if (Nre > 1 and Nre <= 2000):
        
        f_result = 64 / Nre

    elif Nre > 4000:

        delta_result = 1
        delta = 1

        for f in np.arange(0.001, 0.1, 0.001):
            lhs = 1/math.sqrt(f)
            rhs = -2*math.log10((e/d)/3.7 + 2.51/(Nre*math.sqrt(f)))
            delta = abs(lhs - rhs)
            if (delta < delta_result):
                delta_result = delta
                f_result = f
    else:
        pass
    
    return f_result

f = frictionfactor("Commercial Steel", 25, 10000)
print("darcy friction factor = ", f)

When run the code, you will receive the following results.

darcy friction factor =  0.034

Crane differential pressure or flow rate for compressible fluid (SI)

This Python code calculates differential pressure or flow rate according to given pipe & fitting information for compressible fluid.

Calculate differential pressure in pipes based on information on the following pipes/accessory equipment

Flow rate (W) [kg/hr]
Delta Pressure (△P) [kgf/cm2]
Length (L) [m]
Pipe Diameter (d) [mm]
Density (ρ) [kg/m3]
Net Expansion Factor (Y)
Friction factor (f)

Here, Net Expansion Factor (Y) and Friction factor (f) must be calculated using separate charts and programs.

Calculation formulas based on CRANE BOOK

Kpipe = f * l/(d/1000)
K = Kpipe + Kfitting

compressible fluid
△P = (K * W^2) / (1.2646^2 * 0.981 * d^4 * ρ * Y^2)
w = 1.2646 * d^2 * √(△P * 0.981 * ρ / K) * Y

liquid fluid
△P = (K * W^2) / (1.2646^2 * 0.981 * d^4 * ρ)
w = 1.2646 * d^2 * √(△P * 0.981 * ρ / K)

Run Python code below : https://www.mycompiler.io/view/2HpJLHzhEeQ

import math

def compressibledp(W, l, d, r, Y, f, Kf):

    Kp = f*l/(d/1000)

    K = Kp + Kf

    dp = (K*pow(W, 2))/(pow(1.2646, 2)*0.981*pow(d, 4)*r*pow(Y, 2))

    return dp

W = 100   # Flow rate (W, kg/hr)

l = 100   # Length (L, m)

d = 50    # Pipe Diameter (d, mm)

r = 1.24  # Density (ρ, kg/m3)

Y = 1.0   # Net Expansion Factor (Y)

f = 0.005 # Friction factor (f)

Kf = 500  # Resistance coefficient of fittings

dp = compressibledp(W, l, d, r, Y, f, Kf)

print("delta pressure of pipe = ", dp, "kgf/cm2")

def compressibleflow(dp, l, d, r, Y, f, Kf):

    

    Kp = f*l/(d/1000)

    K = Kp + Kf

    W = 1.2646 * pow(d, 2) * math.sqrt(dp * 0.981 * r / K) * Y

    return W

dp = 0.42  # Delta Pressure (△P ,kgf/cm2)

m = 100    # Length (L, m)

d = 50     # Pipe Diameter (d, mm)

r = 1.24   # Density (ρ, kg/m3)

Y = 1      # Net Expansion Factor (Y)

f = 0.005  # Friction factor (f)

Kf = 500   # Resistance coefficient of fittings

flow = compressibleflow(dp, l, d, r, Y, f, Kf)

print("mass flow of pipe = ", flow, "kg/hr")

When run the code, you will receive the following results.

delta pressure of pipe =  0.42 kgf/cm2
mass flow of pipe =  100 kg/hr

Equation of state

Energy engineers frequently calculate the equation of state for gases. The Python example below is an example of calculating the value of an unknown variable using the ideal equation of state and the Van der waals equation of state.

Ideal equation of state

      P = 0.082*T/V

Van Der Waals equation of state

      P = (0.082 * T) / (V - b) - a / V^2

      a = (27 * 0.082^2 * Tc^2) / (64 * Pc)
      b = (0.082 * Tc) / (8 * Pc)

Python code of equation of state 

Run Python code below : https://www.mycompiler.io/view/6RK7hWqB76n

import numpy as np

def eos(p, v, t, pc, tc, EOS, PVT):
 
    delta = 99999; 
    delta_result = 99999; 
 
    if EOS == "Ideal":

        if (PVT == 'P'):
            p = 0.082 * t / v
            result = p
        elif PVT == 'V':
            v = 0.082 * t / p
            result = v
        elif PVT == 'T':
            t = p * v / 0.082
            result = t

    elif EOS == "VanderWaals":

        a = 27*(pow(0.082,2)*pow(tc,2))/(64*pc)
        b = 0.082*tc/(8*pc)

        if PVT == 'P':
            p = (0.082*t)/(v-b)-a/pow(v,2)
            result = p
     
        elif PVT == 'V':
            v = 0.082 * t / p
            v_start = v*1-(v*0.5)
            v_end = v*1+(v*0.5)
            v_interval = v*1/2000

            Iteration = np.arange(v_start, v_end, v_interval)

            for v_temp in Iteration:
             
                p1 = p
                p2 = (0.082*t)/(v_temp-b)-a/pow(v_temp,2)
              
                delta = abs(p1 - p2)

                if (delta < delta_result):
                    delta_result = delta
                    v_result = v_temp                
          
            v_result = round(v_result*1000)/1000
            result = v_result

        elif PVT == 'T':
            t = (p*1+a/pow(v,2))*(v-b)/0.082
            result = t
    else:
        pass

    return result

P = 10          # Pressure (P, atm) = 10
V = 2.403       # Volume/mole (V, L/gmol) = 2.403
T = 298.15      # Temperature (T, K) = 298.15
Pc = 45.80      # Critical Pressure (Pc, atm) = 45.80
Tc = 190.7      # Critical Temperature (Tc, K) = 190.7

print(eos(P, V, T, Pc, Tc, 'Ideal', 'P'))
print(eos(P, V, T, Pc, Tc, 'Ideal', 'V'))
print(eos(P, V, T, Pc, Tc, 'Ideal', 'T'))

print(eos(P, V, T, Pc, Tc, 'VanderWaals', 'P'))
print(eos(P, V, T, Pc, Tc, 'VanderWaals', 'V'))
print(eos(P, V, T, Pc, Tc, 'VanderWaals', 'T'))

When run the code, you will receive the following results.

10.174
2.445
293.049
9.968
2.395

299.072

Unit converter

For energy engineers, unit conversion is the most basic of basics. In the field, there are gauge units and absolute units, so minor mistakes can lead to very different calculation results. Also, there are times when you want to verify your unit conversion.

Python code of unit converter

The Python code below is an example of temperature and pressure engineering unit conversion.
dbUnit = [
         ['C', 'Temperature', 1.0, 0.0],
         ['K', 'Temperature', 1.0, 273.15],
         ['F', 'Temperature', 1.8, 32.0],
         ['R', 'Temperature', 1.8, 491.67],
         ['kPa', 'Pressure',  1.0, 0.0],
         ['MPa', 'Pressure',  0.001, 0.0],
         ['bar', 'Pressure',  0.01, 0.0],
         ['N/m2', 'Pressure', 1000.0, 0.0],
         ['atm', 'Pressure', 0.009869233, 0.0],
         ['at', 'Pressure', 0.01019716, 0.0],
         ['kg/cm2', 'Pressure', 0.01019716, 0.0],
         ['psia', 'Pressure', 0.1450377, 0.0],
         ['lbf/ft2', 'Pressure', 20.88543, 0.0],
         ['torr', 'Pressure', 7.500615, 0.0],
         ['kPag', 'Pressure', 1.0, -101.325],
         ['MPag', 'Pressure', 0.001, -0.101325],
         ['bar_g', 'Pressure', 0.01, -1.01325],
         ['N/m2_g', 'Pressure', 1000.0, -101325.0],
         ['atm_g', 'Pressure', 0.009869, -1.0],
         ['ate', 'Pressure', 0.01019716, -1.033227],
         ['kg/cm2_g', 'Pressure', 0.01019716, -1.033227],
         ['psig', 'Pressure', 0.1450377, -14.696],
         ['lbf/ft2_g', 'Pressure', 20.88543, -2116.216],
         ['torr_g', 'Pressure', 7.500615, -760.0],
         ]

Run Python code below : https://www.mycompiler.io/view/8cFCBv5qrIk

def loadUnit(search):
    for i in range(len(dbUnit)):
        if dbUnit[i][0] == search:
            UnitGroup = dbUnit[i][1]
            UnitScale = float(dbUnit[i][2])
            UnitOffset = float(dbUnit[i][3])

    return UnitGroup, UnitScale , UnitOffset

def unitconv(number, Unit_old, Unit_new):
    
    Group_old, Scale_old, Offset_old = loadUnit(Unit_old)
    Group_new, Scale_new, Offset_new = loadUnit(Unit_new)

    if Scale_old == 0 or Group_old != Group_new: 
        return None

    return ((number - Offset_old) / Scale_old) * Scale_new + Offset_new

print(unitconv(50, "C", "R")) # 50 C = 581.67 R
print(unitconv(10, "bar_g", "kPa")) # 10 bar_g = 1101.325 kPA

Enthalpy chart for Air, N2, O2, H2O, CO, CO2, SO2

Enthalpy chart used when performing enthalpy calculations to calculate the energy loss of flue gas discharged from the stack. Air, N2, O2, H2O, CO, CO2, SO2 enthalpy curves (Btu/lb) are presented as a function of temperature (degF).

Air Preheater in fired heater

Air Preheat System is usually applied to increase a fired heater’s efficiency, and the economics of air preheating should be compared with other forms of flue gas heat recovery. Air preheat systems become more profitable with increasing fuel costs, with increasing process inlet temperature (i.e., higher stack flue gas temperature), and with increasing fired duty.

Economic analysis of preheater

The economic analysis should account for the APH system’s capital costs, operating costs, maintenance costs, fuel savings and the value (if any) of increased capacity. In the case of a system retrofit, the economic analysis should include the cost of incremental heater downtime for the preheat system installation. In addition to economics, the system’s impact on the heater’s operations and maintenance should also be considered. Compared to natural draft systems, air preheat systems typically provide the following operational advantages:

a. Reduced fuel consumption.
b. Improved control of combustion air flow.
c. Reduced oil burner fouling.
d. Better flame pattern control.
e. More complete combustion of difficult fuels. 
Air preheat systems typically have the following operational disadvantages (vs. natural draft systems):
a. Increased radiant section operating temperatures (coil, film, supports, etc.).
b. Increased potential for corrosion of flue gas wetted components downstream of the preheat exchanger, from sulfuric acid condensation.
c. Formation of acid mists, resulting in stack plume, if fuel sulfur content is high.
d. Increased maintenance requirements for mechanical equipment.
e. Increased nitrogen oxide concentration in the flue gas.
f. Reduced stack effluent velocity and dispersion of the flue gases.

In all applications, the use of an air preheat system will increase both the heater’s firebox temperature and radiant flux rate(s). Because of these hotter operating conditions, a thorough review of the heater’s mechanical and process design under APH operations should be performed on all retrofit applications. The hotter firebox temperatures could result in overheated tube supports, guides, tubes, and/or unacceptably high process film temperatures.

In some cases, an air preheat system may provide an increase in fired heater capacity or duty. For example, when a fired heater’s operation is limited by a large flame envelope or poor flame shape (flame impingement on tubes) or by inadequate draft (flue gas removal limitations), the addition of an air preheat system may increase the heater’s capacity.

Based on the flue gas and air flow through the system, the three system types are:

a. Balanced Draft APH Systems (most common type).
b. Forced Draft APH Systems. 
c. Induced Draft APH Systems.

The common “balanced draft” system has both a forced draft (FD) fan and an induced draft (ID) fan. The system is balanced because the combustion air charge, provided by the forced draft fan, is balanced by the flue gas removal of the induced draft fan. In most applications, the forced draft fan is controlled by a “duty controller” that is reset by the heater’s O2 analyzer and the induced draft fan is controlled by an arch pressure controller.

In comparison, the simpler “forced draft” system has only a forced draft fan to provide the heater’s combustion air requirements. All flue gases are removed by stack draft. Because of the low draft generation capabilities of a stack, the exchanger’s flue gas side pressure drop must be kept very low, thus increasing the size and cost of the APH exchanger.

The third and last designation based on fluid flow design is the “induced draft” system, which has only an induced draft fan to remove flue gases from the heater and maintain the appropriate system draft. Combustion air flow is induced by the sub-atmospheric pressure of the heater. In this application, the exchanger must be carefully designed to minimize the combustion air pressure drop while providing the necessary heat transfer.

A typical balanced draft APH system, employing a direct exchanger, is illustrated in figures below.

Air Preheat System Using Regenerative, Recuperative, or Heat Pipe Air Preheater

Air Preheat System Using an Indirect Closed Air Preheater and Mechanical Circulation

Nucleate and film boiling

Boiling of liquid occurs as nucleate or as film boiling. Figure below illustrates a typical flux curve for water. In the region 1-2 the liquid is being heated by natural convection; in 2-3 the nucleate pool boiling occurs with bubbles forming at active sites on the heat transfer surface, natural convection currents set up. 

Q/A varies as Δt. Where n is 3-4, and the peak flux is at point 3 corresponding to the critical Δt for nucleate boiling; in 3-4 film boiling begins; and at 4-5-6 film boiling occurs. In film boiling heat us transferred by conduction and radiation through a film in the heating surface.

Note that the rate of effective heat transfer decreases beyond point 3 and it is for this reason that essentially all process heating-boiling equipment is designed to operate to the left of point 3.

Compressor polytropic efficiency

Polytropic efficiency is a measure of the efficiency of a compressor. It is defined as the ratio of the actual work done by the compressor to the work that would be done if the compression process were isentropic (reversible and adiabatic). The polytropic efficiency of a compressor is always less than its isentropic efficiency.

Formula of Compressor polytropic efficiency

The polytropic efficiency of a compressor can be calculated using the following formula:

Po/Pi = (To/Ti)^(n/(n-1))
n/(n-1) = log(Po/Pi)/log(To/Ti)
η = (n/(n-1))/(k/(k-1))
η = log(Po/Pi)/log(To/Ti)/(k/(k-1))

where η is the polytropic efficiency, Pi and Ti are the inlet pressure and temperature, Po and To are the outlet pressure and temperature, and k is the specific heat ratio of the gas.

η  : polytropic efficiency
Po : outlet pressure (psig)
Ti : inlet temperature (degF)
To : outlet temperature (degF)
n  : compression ratio
k  : specific heat ratio

Python code of Compressor polytropic efficiency

The following Python code is an example of calculating the polytropic efficiency of a compressor when gas introduced at 30 psig, 140 degF is compressed to 125 psig, 350 degF. (specific heat ratio k = 1.243).

Run Python code below : https://www.mycompiler.io/view/JMIFkEueEHW

import math
def compressorpolytropiceff(Pi, Po, Ti, To, k):

Pi : inlet Pressure (psig)

    Po = Po + 14.696
    Pi = Pi + 14.696
    To = To + 460
    Ti = Ti + 460
    if k == 1 or Pi == 0 or Ti == 0: return -1
    return math.log(Po/Pi)/math.log(To/Ti) / (k/(k-1)) * 100

polyeff = compressorpolytropiceff(30, 125, 140, 350, 1.243)
print("compressor polytropic efficiency = ", polyeff)

When run the code, you will receive the following results.

compressor polytropic efficiency =  74.2%

Convert MMSCFH (Million Standard Cubic Feet) to Nm3/hr

MMSCF (Million Standard Cubic Feet) is a imperial unit that generally refers to the volume of natural gas when the temperature and pressure are 60 ℉ (15.56 ℃) and 14.696 psi (101.325 kPa), respectively.

Therefore, when want to find the mass of natural gas, you can find the mass of natural gas by multiplying the volume under these conditions by the standard density.

Calculation example) First, if the volume is converted to normal conditions (0 ℃, 101.325 kPa) based on Korea Gas Corporation,

1 SCF = 1 ft3 (60℉, 14.696 psia) = 0.02832 m3 (0℉, 14.696 psia)
= 0.02679 m3 (0℃, 101.325 kPa) = 0.02679 Nm3

Therefore, 1 MMSCFH = 0.02679 * 10^6 Nm3 = 26,790 Nm3/hr

Fluegas enthalpy calculation

Flue gas emitted from the stack of fired heaters are the main target for recovery of wasted heat. Energy engineers should calculate the amount of flue gas recoverable energy and decide whether to pursue investment projects. At this time, what energy engineers need is the flue gas enthalpy calculation library.

If know the composition of nitrogen, oxygen, argon, carbon dioxide, and water at atmospheric pressure and the discharge temperature, it is a simple procedure to calculate the enthalpy of each component in the mixed component state using the mixing rule.

Python code of flue gas enthalpy

The following Python code calculates the fluegas entalpy emitted at 500 degC with compositions of 70.0 mol% N2, 15.0 mole% O2, 5.0 mole% Ar, 5.0 mole% CO2, and 5.0 mole% H2O.

def fluegas(temp_given, temp_base, mol_n2, mol_o2, mol_arg, mol_co2, mol_h2o):
    if temp_given < 0 or temp_base < 0 or temp_base > temp_given: return -1
    coeff_a = [6.921525928892, 6.905598931009, 4.964694915254, 8.44682086499, 7.967187267909]
    coeff_b = [0.000264324135, 0.00156659667, 0, 0.00575849121, 0.000876904142]
    coeff_c = [0.000000433542006, -0.000000453431059, 0, -0.00000215927489, 0.000000575161121]
    coeff_d = [-1.17644032E-10, 5.37208504E-11, 0, 3.05898491E-10, -1.47239368E-10]
    mw = [28.0134, 31.9988, 39.948, 44.01, 18.0152]
    wt = [0, 0, 0, 0, 0] 
    temp_given = temp_given * 1.8 + 32
    temp_base = temp_base * 1.8 + 32
    sum_mol = mol_n2 + mol_o2 + mol_arg + mol_co2 + mol_h2o
    mol_n2 = mol_n2 / sum_mol
    mol_o2 = mol_o2 / sum_mol
    mol_arg = mol_arg / sum_mol
    mol_co2 = mol_co2 / sum_mol
    mol_h2o = mol_h2o / sum_mol
    sum_mw = mw[0] * mol_n2 + mw[1] * mol_o2 + mw[2] * mol_arg + mw[3] * mol_co2 + mw[4] * mol_h2o
    wt[0] = mw[0] * mol_n2 / sum_mw
    wt[1] = mw[1] * mol_o2 / sum_mw
    wt[2] = mw[2] * mol_arg / sum_mw
    wt[3] = mw[3] * mol_co2 / sum_mw
    wt[4] = mw[4] * mol_h2o / sum_mw
    sum_h = 0
    for comp in range(5):
        each_ha = coeff_a[comp] * (temp_given - temp_base)
        each_hb = coeff_b[comp] / 2 * (pow(temp_given, 2) - pow(temp_base, 2))
        each_hc = coeff_c[comp] / 3 * (pow(temp_given, 3) - pow(temp_base, 3))
        each_hd = coeff_d[comp] / 4 * (pow(temp_given, 4) - pow(temp_base, 4))
        each_h = (each_ha + each_hb + each_hc + each_hd) / mw[comp] * wt[comp] * 0.556
        sum_h = sum_h + each_h
    return sum_h
enthalpy = fluegas(500,15,0.70, 0.15, 0.05, 0.05, 0.05) 
print("flue gas enthalpy = ", enthalpy, "kcal/kg")

When run the code, you will receive the following results.

flue gas enthalpy =  121.0 kcal/kg

Normal state and Standard state

When dealing with gas flow rates, normal condition, standard condition, and actual condition must be handled carefully. 

Definition of Normal state and Standard state

Normal state is 0℃, 1 atm, and Nm3/hr units are usually used.

Standard state is 15℃ (60℉), 1 atm.

Example of Normal state and Standard state

If the flow rate is 5,000 Nm3/hr at 1 bar and 0 ℃, What is the flow rate at 7 bar and 20 ℃?

Normal condition
Pn = 1 bar (=absolute pressure)
Tn = 0 ℃
Vn = 5000 Nm3/hr

Actual condition
Pa = 7 bar
Ta = 20 ℃
Va = ?

Pn*Vn/Tn = Pa*Va/Ta
Va = Pn/Pa*Ta/Tn*Vn
Va = 1/3*(273+50)/(273+0)*5000 = 1,972 m3/hr

Adiabatic compression power calculation

A compressor is equipment to raise the pressure of a compressible fluid. Work is required as an input to the compressor in this process. User can calculate theoretical energy consumption through the calculation below.

Formula of adiabatic compression

W = Pi*Vi*(k/(k-1))*[(Pe/Pi)^((k-1)/k)-1]

Where,
W : fluid or gas power (Watt)
Vi : inlet or suction flow (m3/s, actual condition)
Pe : exit or discharge pressure (N/m2)
ηc : compressor efficiency

Example of adiabatic compression

If an adiabatic compressor compresses 5,000 m3/hr of atmospheric air to 10 barg, how much energy will it consume (adiabatic efficiency η = 75%, ratio of specific heats k = 1.4)?

Pi : inlet or suction pressure (N/m2)

k : ratio of specific heats = cp/cv ()

Vi = 5,000 m3/hr = 1.389 m3/s

k = 1.4

Pi = 0 barg = 101,325 Pa
Pe = 10 barg = 1,101,325 Pa
η = 75%
W = Pi*Vi*(k/(k-1))*[(Pe/Pi)^((k-1)/k)-1]/η
W = 101,325*1.389*1.4/(1.4-1)*((1,101,325/101,325)^((1.4-1)/1.4)-1)/0.75 = 

641,832 W = 642 kW