Author: CheGuide

Shortcut Distillation

Shortcut Distillation

Distillation is the most common unit operation for separating liquid mixtures into valuable and high purity products. The normal procedure for solving a typical multi-component distillation problem is to solve the MESH (Material balance, Equilibrium, Summation and Heat) balance equations stage-by-stage.

Though computer programs are normally available for the rigorous solution of the MESH equations, short-cut methods are still useful in the preliminary design work, and as an aid in defining problems for computer simulation.

This article describes a widely used shortcut distillation method commonly referred to as the Fenske-Underwood-Gilliland (FUG) method.

  1. The components that have their distillate and bottoms fractional recoveries specified are called key elements. The most volatile of the keys is called the Light Key (LK) and the least volatile is called the heavy key (HK). The other components are called non-keys(NK). First step is to choose LK and HK and fix their distribution in top and bottom product.
  2. Next step is to estimate overall top and bottom flow rate compositions assuming only light non-key components in top and heavy non-key component in bottoms as first attempt.
  3. Estimate Dew point TDew for top composition and Bubble Point TBubble for bottom composition.
  4. Estimate relative volatility with respect to the HK of all components at the TDew (top) and TBubble (bottom).
    αi = Ki / KHK

    Calculate average volatilities for all components.

    αi,Avg = (αi, Top * αi, Bottom)0.5
  5. Recalculate overall top and bottom flow rate compositions on the basis of the Hengsteback and Geddes equation.
    log(di/bi) = A + B*log(αi)

    Recovery of heavy key in bottoms is defined as bHK/fHK, recovery of light key in distillate is defined as dLK/fLK. Constants A & B are defined as:

    A = log((1-bHK/fHK)/(bHK/fHK))
    B = log[(dLK/fLK)/(1-dLK/fLK)*(bHK/fHK)/(1-bHK/fHK)]/log(αLK)

    Recovery of ith component in distillate is defined as:

    (di/fi) = (10AiB)/(1 + 10AiB)

    Recovery of ith component in bottom is defined as:

    (bi/fi) = 1 - (di/fi)

    After calculating top and bottom compositions, go back to Step (3) and update TDew and TBubble and recalculate relative volatilites.

  6. Estimate minimum number of stages, Nm using Fenske equation.
    Nm = log[(xLK/xHK)d.(xHK/xLK)b]/log(αLK)

    where, (xLK/xHK)d is the ratio of mol fraction of light key and heavy key in distillate and
    (xHK/xLK)b is the ratio of mol fraction of heavy key and light key in bottoms.

  7. Estimate minimum reflux ratio, Rm using Underwood’s equation.
    Rm + 1 = Σ(αi.xi,d)/(αi - θ)

    where Underwood constant, θ is calculated as the root of following equation:

    1 - q = Σ(αi.xi,f)/(αi - θ)

    where, 1 < θ < αLK and q is the quality of feed.

  8. Select an operating reflux ratio, R as a factor of minimum reflux ratio, typically a factor of 1.2 is considered. Number of stages, N is calculated based on Gilliland’s correlation.
    X = (R - Rm)/(R + 1)
    Y = (N - Nm)/(N + 1)
    Y = 1 - exp[(1 + 54.4*X)/(11 + 117.2*X)*(X - 1)/X0.5]
  9. Feed stage location is estimated based on Kirkbride’s equation.
    m/p = [(B/D)*(xHK/xLK)F*(xLK,B/xHK,D)²]0.206
    N = m + p

    where, m is number of stages above feed plate and n is number of stages below feed plate.

Web based calculation available at


Spreadsheet for Shortcut Multi-Component Distillation based on FUG method


Vapor Liquid Vertical Separator

Vapor Liquid Vertical Separator

Vapor Liquid separators are one of the most common types of process equipment. This article provides step-by-step procedure for two phase vertical vapor liquid separator design.

Vapor liquid separation occurs in three stages.

  • The first stage, primary separation uses an inlet diverter so that momentum of the liquid entrained in the vapor causes the largest droplets to impinge on the diverter and then drop by gravity.
  • The next phase, secondary separation is gravity separation of smaller droplets as the vapor flows through the disengagement area.
  • The final stage is mist elimination where the smallest droplets are formed which will separate by gravity.

For secondary separation, maximum allowable vapor velocity is calculated by using the Souders-Brown equation.

V = K ((ρL - ρV)/ ρV)0.5

where, ρL, ρV are liquid and vapor densities. K is calculated as following.

K = 0.35 - 0.01(P-100)/100

where P is in psig, K is in feet/sec. K value is multiplied by a factor of 0.6 for Glycol or Amine solution, 0.7 for compressor suction scrubber, 0.7 for expander inlet separators and 0.5 when there is no mist eliminator.

Vessel diameter is calculated based on vapor velocity V.

D = (4*Q / (π*V))0.5

where, Q is the vapor volumetric flowrate. For separator with mesh pad, width of support ring is added to diameter calculated above. It varies from 0.3 ft to 1 ft depending on diameter of mesh pad. Clearance of 4 inch is provided around the mesh pad for installation which increases the diameter by 0.66 ft.

Separator Height

Vertical Vapor Liquid Separator

h1 : Low Liquid Level to Bottom Tangent Line

For operating pressure less than 1000 psig (70 barg), it is taken as 6 inch (150 mm) else it is kept 12 inch (300 mm).

h2 : Level Range

It is calculated based on holdup time required and slug capacity of the separator. Holdup volume is calculated based on liquid flowrate multiplied by holdup (residence) time.

h2 = 4*(VHoldup + VSlug)/(π.D²)

Minimum height is taken as 18 inch (450 mm) for proper level control.

h3: High Liquid Level to Inlet Nozzle

For separator with half pipe at inlet nozzle it is taken as maximum of 0.3D or 12 inch (300 mm).

h4 : Diameter of Inlet Nozzle

This height is taken based on outside diameter of inlet nozzle.

h5 : Inlet nozzle to Mist Eliminator/ Top Tangent Line

For separator without mist eliminator this height is taken as maximum of 0.9D or 36 inch (900 mm). For separator with mist eliminator this height is taken as maximum of 0.45D or 24 inch (600 mm)

h6 : Height of Mist Eliminator

Height of mist eliminator is typically 6 inch (150 mm)

h7 : Mist Eliminator to Top Tangent Line

This height is taken as maximum of 0.15D or 12 inch (300 mm).

All heights are summed up to calculate overall length of separator. L/D ratio is typically kept in the range of 3 to 5. If it is out of this range, new diameter is selected and all above steps are repeated till desired L/D ratio is obtained.


Spreadsheet for design sizing of Vapor Liquid Vertical Separator


Pump Sizing Calculation

Pump Sizing Calculation

Pump is a most common equipment used in a chemical plant to transfer fluid from one location to another. This article shows how to do pump sizing calculation to determine differential head required to be generated by pump based on suction and discharge conditions.

Pump Sizing Calculation

Suction Pressure

Pressure at pump suction is calculated as following

PSuct. = P1 + Pstatic - ΔPEquipment - ΔPfriction


  • P1 is pressure at liquid surface in suction vessel.
  • Pstatic is pressure due to height of liquid level above pump suction.
  • ΔPEquipment is pressure drop in an equipment at pump suction like strainers, filters etc.
  • ΔPfriction is pressure drop due to suction pipe and fittings.


Pstatic (psi) = h1.(SG)/2.31
Pstatic (bar) = h1.ρ.g/100000

where, h1 is height of liquid above pump suction, SG is specific gravity of liquid, ρ is liquid density (kg/m³) and g is gravitational constant (9.81).


ΔPfriction = ΔPPipe + ΔPFittings

ΔPPipe is pressure drop in a pipe due to single phase fluid flow. ΔPFittings is pressure drop due to pipe fittings, which can be calculated based on 2-K & 3-K method.

Discharge Pressure

Pressure at pump discharge is calculated as following

PDisch. = P2 + Pstatic + ΔPEquipment + ΔPfriction


  • P2 is pressure at liquid surface in discharge vessel.
  • Pstatic is pressure due to height of liquid level above pump suction.
  • ΔPEquipment is pressure drop in equipment at pump discharge like heat exchangers, control valve, flowmeter, valves etc.
  • ΔPfriction is pressure drop due to suction pipe and fittings.


Pstatic (psi) = h2.(SG)/2.31
Pstatic (bar) = h2.ρ.g/100000

where, h2 is height of liquid above pump suction at which liquid is to be discharged.

ΔPfriction is calculated in similar way as mentioned above for discharge piping.

Differential Head

Differential head required to be generated by pump is calculated as following.

Head (ft) = (PDisch. - PSuct.)*2.31/SG
Head (m) = (PDisch. - PSuct.)*100000/(ρ.g)

Hydraulic Power

Hydraulic power is calculated as following.

Power (bhp) = Head(ft) * Flow(gpm) * SG / 3960
Power (kW) = Head(m) * Flow(m³/h) * SG * g / 3600

NPSH Available

Net Positive Suction Head (NPSH) available is calculated as following.

NPSH Avail.(ft) = (PSuct. - PVapor)* 2.31 / SG
NPSH Avail.(m) = (PSuct. - PVapor)* 100000 / (ρ.g)

where, PVapor is vapor pressure of liquid at suction vessel temperature.

System Head Curve

System head curve is a graphical representation of the relationship between flow and hydraulic losses in a given piping system. It is prepared by calculating differential head as specified above at different flow. The intersection of the pump manufacturer’s curve with system curve defines the operating point of the pump.

Pump System Curve

Pump curve and system curve can be fitted into a second order polynomial. Operating point is calculated by solving these equations for positive roots.

Web based calculation available at


Spreadsheet for Pump Sizing Calculation

Property Estimation Joback Method

Property Estimation Joback Method

The Joback method is a group contribution method which calculates thermophysical and transport properties as a function of the sum of group parameters. It uses a very simple and easy to assign group scheme.

This article shows how to calculate properties using Joback Method in an excel spreadsheet.

Estimate properties for 1-Butanol based on Joback Method

Structure of 1-Butanol consists of following groups

-CH3 Group : 1
-CH2 Group : 3
-OH (alcohol) Group : 1

Normal Boiling Point

TNBP (K) = 198 + ΣTb,i

where Tb,i is contribution due to each group. These values for individual group are available in literature and Wikipedia reference below. On combining values for each group normal boiling point comes out to be –

TNBP (K) = 383.10

Critical Temperature

Tc (K) = TNBP /(0.584 + 0.965*ΣTc,i - (ΣTc,i)²)
Tc (K) = 545.08

Critical Pressure

Pc (bar) = (0.113 + 0.0032*NA - ΣPc,i)-2
Pc (bar) = 43.86

where NA is number of atoms in the molecular structure.
Critical Volume

Vc (cm³/mol) = 17.5 + ΣVc,i
Vc (cm³/mol) = 278.50

Critical Compressibility

Zc = (Pc.Vc)/(R.Tc)
Zc = 0.2695

Acentric Factor
Lee-Kesler method can be used to estimate the acentric factor.

ω = α/β
θ = TNBP / Tc
α = -ln(Pc) - 5.92714 + 6.09648/θ + 1.28862.ln(θ) - 0.169347.θ6
β = 15.2518 - 15.6875/θ - 13.4721.ln(θ) + 0.43577.θ6
ω = 0.6602

Freezing Point

Tm (K) = 122.5 + ΣTm,i
Tm (K) = 195.66

Heat of Formation (Ideal Gas, 298 K)

Hformation (kJ/mol) = 68.29 + ΣHform,i
Hformation (kJ/mol) = -278.12

Gibbs Energy of Formation (Ideal Gas, 298 K)

Gformation (kJ/mol) = 53.88 + ΣGform,i
Gformation (kJ/mol) = -154.02

Heat of Vaporization (at Normal Boiling Point)
Reidel’s equation can be used to estimate a liquid’s heat of vaporization at its normal boiling point.

ΔHv (kJ/mol) = 1.092.R.TNBP .(ln(Pc) -1.013)/(0.930 - (TNBP / Tc))
ΔHv (kJ/mol) = 42.38

Heat of Fusion

ΔHfus (kJ/mol) = -0.88 + ΣHfus,i
ΔHfus (kJ/mol) = 10.2

Heat Capacity (Ideal Gas)

CP (J/mol.K) = A + B.T + C.T² + D.T³
A = Σai - 37.93
B = Σbi + 0.210
C = Σci - 3.91*10-4
D = Σdi + 2.06*10-7
CP (J/mol.K) = 110.96 (at 300 K)

Heat of Vaporization (@ Temperature T)
Watson equation can be used to estimate heat of vaporization at different temperature as following.

ΔHv,2 (kJ/mol) = Hv,1*((Tc - T2)/(Tc - T1))0.38
ΔHv,2 (kJ/mol) = 49.60 (at 300 K)

Liquid Viscosity

ηL (Pa.s) = Mw.exp(A/T + B)
A = Σηa - 597.82
B = Σηb - 11.202
ηL (Pa.s) = 1.936*10-3 (at 300 K)

where Mw is the molecular weight.
Liquid Density
Rackett equation can be used to estimate liquid density from critical properties as following –

ρL (gm/cm³) = Mw/( (R.Tc/Pc)*Zc^(1 + (1- T/Tc)2/7))
ρL (gm/cm³) = 0.756 (at 300 K)

Liquid Vapor Pressure
Lee-Kesler equation can be used to estimate liquid vapor pressure as following –

Pvp (bar) = Pc.exp(f0 + ω.f1)
f0 = 5.92714 - 6.09648/Tr - 1.28862*ln(Tr) + 0.169347.Tr6
f1 = 15.2518 - 15.6875/Tr - 13.4721*ln(Tr) + 0.43577.Tr6
Tr = T / Tc
Pvp (bar) = 0.018

Web based calculation available at


Spreadsheet for property estimation based on Joback Method


  1. Joback Method – Wikipedia
  2. Molecular Knowledge Systems
  3. Chemstation Physical Properties User Guide
Heat Exchanger Rating (Bell-Delaware Method)

Heat Exchanger Rating (Bell-Delaware Method)

Stream Analysis Heat Exchanger

In Bell Delaware method, the fluid flow in the shell is divided into a number of individual streams. Each of these streams introduces a correction factor which is used to correct heat transfer coefficient and pressure drop across the shell. This article gives step-by-step guidance on doing heat exchanger rating analysis based on Bell-Delware method.

Shell Side Heat Transfer Coefficient, hs

Cross flow area, Sm is the minimum flow area in one baffle space at the center of the tube bundle. It is calculated by following equation:

Sm = B[(Ds - DOTL) + (DOTL - Do)(PT - Do)/PT,eff ]

where, PT is tube pitch, B is central baffle spacing, DOTL is outer tube limit diameter, Ds is shell diameter and Do is tube outside diameter.

PT,eff = PT for 30° and 90° layouts 
PT,eff = 0.707*PT for 45° layout

Shell side cross flow mass velocity, GS is defined as:

GS = mS/Sm

where, mS is shell side mass flow rate. Shell side Reynolds number ReS is then calculated from

ReS = Do.GS / μS

where, μS is the shell side fluid dynamic viscosity at average bulk temperature.

Shell side Prandtl number PrS is calculated as following :

PrS = CP,SS / kS

where, CP,S is the shell side fluid specific heat and kS is the shell side fluid thermal conductivity.

Colburn j-factor for an ideal tube bank is defined as:
Colburn j-factor
where a1, a2, a3 and a4 are correlation constants listed below.
Correlation constants j factor
The ideal tube bank based coefficient is calculated from –
ideal heat transfer coefficient
where, μS,W is shell side fluid viscosity at wall temperature.

Correction factor for Baffle Window Flow, JC

The factor JC accounts for heat transfer in the baffle windows. It has a value of 1.0 for exchanger with no tubes in the windows.

JC = 0.55 + 0.72FC 
FC = 1 - 2FW 
FW = (θCTL - Sin(θCTL))/2π 
θCTL = 2cos-1(Ds(1 - 2*Bc/100)/DCTL) 

where, Bc is segemental baffle cut in %.

Correction factor for Baffle Leakage, JL

The correction factor JL considers the effects of the tube-to-baffle and shell-to-baffle leakage streams on heat transfer.

JL = 0.44(1-rS) + (1-0.44(1-rS))exp(-2.2rL)
rS = Ssb /(Ssb + Stb) 
rL = (Ssb + Stb)/ Sm 
Ssb = Ds*DSB(π - 0.5θDS) 
Stb = (π/4)((Do+LTB)2 - Do2)Nt(1-FW)
θDS = 2cos-1(1 - 2Bc/100)

where, Nt is number of tubes, DSB is diametral clearance between shell & baffle and LTB is diametral clearance between tube and baffle.

Correction factor for Bundle Bypass effects, JB

Bundle bypass correction factor JB accounts for the bundle bypass stream flowing in the gap between the outermost tubes and the shell. The number of effective rows crossed in one cross flow section, Ntcc between the baffle tips is provided by following equation.

Ntcc = (Ds/Pp)(1 - 2Bc/100) 
Pp = PT 30.5/2 for 30° layout 
Pp = PT / 20.5 for 45° layout 
Pp = PT for 90° layout

Ratio of sealing strips to tube rows rss is provided by

rss = Nss/ Ntcc

where Nss is number of sealing strips (pairs) in one baffle.

The bundle bypass flow area, Sb is defind as

Sb = B(Ds - DOTL - Do/2)

where, B is central baffle spacing. Correction factor JB is then calculated as following –

JB = exp(-Cj(Sb / Sm)(1 - (2rss)1/3)) for rss < 0.5 
JB = 1 for rss >= 0.5 
Cj = 1.35 for ReS < 100 
Cj = 1.25 for ReS >= 100

Correction factor for adverse temperature gradient, JR

The factor JR accounts for the decrease in the heat transfer coefficient with downstream distance in laminar flow.

Ntcw = (0.8/Pp)(Ds(Bc/100) - (Ds-(DOTL-Do))/2 )
NB = 1 + (int)(L - 2Ls - LBIn - LBOut)/B 
NC = (Ntcw + Ntcc)(1 + NB) 
JRL = (10/NC)0.18 
JR = 1, ReS > 100 
JR = JRL + (20-ReS)(JRL - 1)/80, ReS <= 100, ReS > 20 
JR = JRL, ReS <= 20

where, L is tube length, Ls is tubesheet thickness, LBIn is inlet baffle spacing and LBOut is outlet baffle spacing.

Correction factor for unequal baffle spacing, JS

n1 = 0.6, ReS >= 100 
n1 = 1/3, ReS < 100 
JS = ((NB-1)+(LBIn/B)1-n1 + (LBOut/B)1-n1)/((NB-1)+(LBIn/B) + (LBOut/B))

Shell side heat transfer coefficient is calculated as

hs = hIdeal(JC.JL.JB.JS.JR)

Shell Side Pressure Drop, ΔPs

Friction factor for ideal tube bank is calculated as following –
friction factor for ideal tube bank
where b1, b2, b3 and b4 are correlational constants listed below.
friction factor constants
Pressure drop for an ideal tube bank is calculated from

ΔPIdeal = 2f(GS²/ρS)(μSS,W)0.14 Ntcc

Correction factor for Baffle Leakage, RL

RL = exp(-1.33(1+rS)rLp) 
p = 0.8 - 0.15(1+rS)

Pressure drop for window section, ΔPW

Following terms are calculated as –

SWG = (Ds²/8)(θDS - Sin(θDS)) 
SWT = Nt.FW(πDo²/4) 
GW = mS/(Sm.SW)0.5 
DW = 4.SW /(π.Do.Nt.FW + θDS.Ds)

Pressure drop for laminar and turbulent flow is calculated.

ΔPW, Turb = NB.RL(2+0.6*Ntcw).GW²/(2.ρS)

pressure drop for laminar flow

ΔPW = ΔPW, Turb , ReS >= 100 
ΔPW = ΔPW,Laminar , ReS < 100

Correction factor for Bundle Bypass effect, RB

RB = exp(-Cr(Sb / Sm)(1 - (2rss)1/3)) for rss < 0.5 
RB = 1 for rss >= 0.5 
Cr = 4.5 for ReS < 100 
Cr = 3.7 for ReS >= 100

Correction factor for unequal baffle spacing, RS

n = 0.2, ReS >= 100 
n = 1.0, ReS < 100 
RS = 0.5((B/LBIn)2-n + (B/LBOut)2-n)

Pressure drop in Central Baffle spaces, ΔPC is defined as –

ΔPC = (NB - 1)ΔPIdeal.RL.RB

Pressure drop in entrance & exit baffle spaces, ΔPE is calculated as –

ΔPE = ΔPIdeal(1 + Ntcw/Ntcc).RB.RS

Shell side pressure drop is calculated as following –


Tube Side Heat Transfer Coefficient, ht

Reynold’s number and Prandtl number are calculated as following –

ReT = Di.v.ρtt 
PrT = Cp,tt/kt

where, Di is tube inside diameter, v is velocity, ρt is density, μt is viscosity, kt is thermal conductivity and Cp,t is specific heat for fluid on tube side.

For laminar flow, ReT < 2300, Sieder and Tate correlation is used for Nusselt’s nubmer.

Nu = 1.86(ReT.PrT.Di/Leff)1/3 
Leff = L - 2*Ls

For turbulent flow, ReT > 10,000, following equation developed by Petukhov-Kirillov can be used.

Nu = (f/2)ReT.PrT/(1.07+12.7(f/2)0.5(PrT2/3-1)) 
f = (1.58 ln(ReT) - 3.28)-2

For transient flow, Nusselt number can be interpolated from Nu Laminar & Nu Turbulent.

Heat transfer coefficient is calculated as following –

ht = Nu.(kt/Di)(μtt, w)0.14

Tube Side Pressure Drop, ΔPt

Tube side pressure drop is calculated by following equation –

ΔPt = (4.f.Leff.Np/Di + 4.Np)ρt.v²/2

where, Np is number of tube passes.

Overall Heat Transfer Coefficient, U

Resistance due to tube wall is calculated as following

Rtube = Do/(2.ln(Do/Di).ktube)

where, ktube is thermal conductivity of tube material. Overall clean heat transfer coefficient, UClean is calculated as per below equation

UClean = 1/(hS + Do/( + Rtube)

Overall dirty heat transfer coefficient, UDirty is calculated as per below expression

UDirty = 1/(1/UClean + fshell + ftube)

where, fshell & ftube are fouling factors for shell and tube side.

Heat transfer coefficient required, URequired is calculated as following

URequired = Q /(A x LMTDcorrected)

where, Q is heat duty, A is heat transfer area and LMTDcorrected is corrected logarithmic mean temperature difference.

Over Surface, % = (UClean/URequired - 1)*100 
Over Design,  % = (UDirty/URequired - 1)*100

Web based calculation available at


Spreadsheet for Heat Exchanger Rating based on Bell-Delaware Method


LMTD Correction Factor Charts

LMTD Correction Factor Charts

LMTD Correction factor charts

Heat transfer rate in the exchanger is represented by

q = U * A * F * LMTD

here F (< 1) is interpreted as a geometric correction factor, that when applied to the LMTD (Log Mean Temperature Difference) of a counter flow heat exchanger, provides the effective temperature difference of the heat exchanger under consideration.

It is a measure of the heat exchanger’s departure from the ideal behavior of a counter flow heat exchanger having the same terminal temperatures. The F-LMTD method is widely used in heat exchanger analysis, particularly for heat exchanger selection, (sizing problems) when as a result of the process requirements the temperatures are known and the size of the heat exchanger is required.

Log Mean Temperature Difference is defined as

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


ΔT1 = T1 - t2
ΔT1 = T2 - t1

T1, T2 are inlet and outlet temperature of Fluid 1;  t1, t2 are inlet and outlet temperature of Fluid 2.

Log Mean Temperature Difference Correction Factor F is dependent on temperature effectiveness P and heat capacity rate ratio R for a given flow arrangement. Temperature effectiveness P is different for each fluid of a two fluid exchanger.

For fluid 1, it is defined as the ratio of the temperature range of fluid 1 to the inlet temperature difference.

P1 = ( T2 - T1 ) / ( t1 - T1 )

Heat Capacity Ratio R is defined as

R1 = ( t1 - t2 ) / ( T2 - T1 )

N (Shell) – 2M (Tube) Pass Tema E

Following general equation is used for shell and tube heat exchanger having N shell passes and 2M tube passes per shell.

S = (R1² + 1)0.5 / (R1 - 1)
W = [(1 - P1.R1)/(1 - P1)]1/N
F = S.ln(W)/ ln[( 1 + W - S + S.W) /( 1 + W + S - S.W)]

For limiting case of R1 = 1,

W' = (N - N.P1)/( N - N.P1 + P1 )
F = 20.5 [(1 - W')/ W' ]/ln[( W'/(1-W') + 1/20.5)/( W'/(1-W') - 1/20.5)]

For plotting the correction factor charts P1 values are listed from 0.01 to 1 with increment of 0.01 and then F values are calculated for each P1 and R1 based on above equations.

For following type of exchangers F values depend on NTU along with P1 and R1, where NTU is defined as number of transfer units. Range of NTU values are listed from 0.01 to 32 and F value is calculated for each value. LMTD Correction factor, F is determined as following –

F = [ 1/(NTU1.(1 - R1))].ln[(1 - R1.P1)/(1 - P1)]  ---- (1)

For limiting case of R1 = 1,

F = P1 / [NTU1 . (1 - P1)] ---- (2)

Cross Flow Fluid 1 Unmixed Tema X

Following relation is used to calculate P1 using NTU.

K = 1 - exp(-NTU1)
P1 = [1 - exp(-K.R1)]/ R1

F factor is calculated as per equations (1) & (2).

Cross Flow Fluid 2 Unmixed Tema X

Following relation is used to calculate P1 using NTU.

K = 1 - exp(-R1.NTU1)
P1 = 1 - exp(-K/R1)

F factor is calculated as per equations (1) & (2).

Cross Flow Both Fluid Mixed Tema X

Following relation is used to calculate P1 using NTU.

K1 = 1 - exp(-NTU1)
K2 = 1 - exp(-R1.NTU1)
P1 = 1/[1/K1 + R1/K2 - 1/NTU1]

F factor is calculated as per equations (1) & (2).

In a similar way, LMTD correction charts can be prepared for different type of geometries based on relation between NTU, P & R.


Spreadsheet for LMTD Correction factor charts.


  • Journal of Heat Transfer, Vol. 125, June 2003 – A General Expression for the Determination of the Log Mean Temperature Correction Factor for Shell and Tube Heat Exchangers – Ahmad Fakheri
  • Handbook of Heat Transfer 3rd Edition – Chapter 17 – Heat Exchangers by R.K. Shah and D.P. Sekulic
Jacketed Vessel Heat Transfer (Half Pipe Coil)

Jacketed Vessel Heat Transfer (Half Pipe Coil)

Agitator equipped vessels with half pipe coil jackets are widely used in variety of process applications. This article shows how to calculate heat transfer in an agitated vessel provided with an external half pipe coil jacket.

Overall heat transfer coefficient, U is defined as

1/U = 1/hi + ffi + x/k + ffo + 1/ho


  • hi : film coefficient process side
  • ho : film coefficient coil side
  • ffi : fouling factor process side
  • ffo : fouling factor coil side
  • x : vessel wall thickness
  • k : vessel wall thermal conductivity

Process Side, hi

Process side film coefficient, hi depends upon type of impeller, Reynold’s number (Re) and Prandtl number (Pr).

Re = D².N.ρ / μ
Pr = Cp.μ/k

where, D is impellor diameter, N is impellor rpm, ρ is fluid density, μ is fluid viscosity, Cp is fluid specific heat and k is fluid thermal conductivity at bulk fluid temperature. hi is defined as following :

Nui = C.Rea. Prb. (μ/ μw)c. Gc 
hi.DT/k = Nui

where constants C, a, b & c are available in literature for different type of impellors. DT is vessel diameter. Gc is geometric correction factor for non-standard geometries. μ/ μw is viscosity correction factor due to difference in viscosities at bulk fluid and wall temperatures. These constants are available in references mentioned below.

Coil Side, ho

Pipe coils are made with a 180° central angle or a 120° central angle. Equivalent diameter (De) and Flow area (Ax) is defined as following.

180° Coil

De = (π/2).dci
Ax = (π/8).dci²

120° Coil

De = 0.708 dci
Ax = 0.154 dci²

where, dci is inner diameter of pipe.

Reynold’s and Prandtl number are calculated based on jacket fluid properties and velocities.

Re = De.v.ρ / μ
Pr = Cp.μ / k

where, ρ is coil fluid density, μ is coil fluid viscosity and k is coil fluid thermal conductivity. v is fluid velocity in coil.

Turbulent Flow

For Re > 10000

Nuc = 0.027 Re0.8 Pr0.33 (μ/ μw)0.14 (1 +3.5 De/Dc)

where Dc is the mean or centerline diameter of the coil. Coil outer diameter (Do) is determined as following.

180° Coil : Do = DT + 2(dci/2) + 2.x
120° Coil : Do = DT + 2(dci/4) + 2.x
Dc = ( Do + DT) / 2

Laminar Flow

For Re < 2100

Nuc = 1.86 [ Re.Pr.De/L ]0.33 (μ/ μw)0.14

where, L is coil length along the vessel.

Transient Flow

For 2100 < Re < 10000
Calculate NuLaminar based on Re = 2100 and NuTurbulent based on Re = 10,000.

NuTransient = NuLaminar + (NuTurbulent - NuLaminar).(Re - 2100)/(10000 - 2100)

ho is determined as following

ho.De/k = Nuc

ho and hi thus calculated are used to get value of overall heat transfer coefficient U.

Spreadsheet for Half Pipe Coil Agitated Vessel Heat Transfer


  • Heat Transfer Design Methods – John J. McKetta (1992)
  • Heat Transfer in Agitated Jacketed Vessels – Robert F. Dream, Chemical Engineering, January 1999
Estimating Binary Interaction Parameter for Peng Robinson EOS

Estimating Binary Interaction Parameter for Peng Robinson EOS

This article shows how to estimate binary interaction parameter (BIP) used in Peng Robinson (PR) Equation of State (EOS) from experimental data by regression in excel spreadsheet.

Determine binary interaction parameter used in Peng Robinson EOS for a mixture of acetone and chloroform from Tx-y experimental data available at 101.33 kPa.

VLE data acetone chloroform mixture

Obtain pure component properties like critical temperature (Tc), critical pressure (Pc) and accentric factor ω from literature. Experimental K1 is determined from x, y data.

K1Experimental = y1 / x1

Peng Robinson EOS parameters are calculated.

κi = 0.37464 + 1.54226ω - 0.26992ω²
αi = [ 1 + κi (1 - (T/Tc)0.5)]²
ai = 0.45724 (RTc)²α / Pc
bi = 0.07780 RTc / Pc

Mixture parameters are calculated next

a12= [(a1.a2)0.5(1 - k12)] = a21
a = a1.x1² + 2.a12.x1.x2 + a2.x2²
b = b1.x1 + b2.x2
A = aP/(RT)²
B = bP/RT

where k12 is Binary Interaction Parameter (BIP) to be estimated. For first run, use a guess value to start calculation.

Following cubic equation is solved to get ZL.

Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0

Above equation can be written as

Z³ + C2.Z² + C1.Z + C0 = 0

Solving Cubic Equation

Cubic equation is solved using following procedure. Calculate Q1, P1 & D.

Q1 = C2.C1/6 - C0/2 - C2³/27
P1 = C2²/9 - C1/3
D = Q1² - P1³

If D >= 0, then equation has only one real root provided by

Z1 = (Q1 + D0.5)1/3 + (Q1 - D0.5)1/3 - C2/3

If D < 0, then equation has 3 real roots, following parameters are calculated

t1 = Q1² / P1³
t2 = (1 - t1)0.5 / t10.5. Q1/abs(Q1)
θ = atan(t2)

Roots are calculated as following –

Z0 = 2.P10.5.cos(θ/3) - C2/3
Z1 = 2.P10.5.cos((θ + 2*Π)/3) - C2/3
Z2 = 2.P10.5.cos((θ + 4*Π)/3) - C2/3

Roots thus calculated are arranged in descending order, highest root gives ZV and lowest root gives ZL.


Based on ZL, liquid fugacity φiL is calculated for each data set.

Liquid fugacity phi using Peng Robinson EOS

In next step, Vapor phase fugacity is calculated. Mixture properties are estimated as following

a = a1.y1² + 2.a12.y1.y2 + a2.y2²
b = b1.y1 + b2.y2
A = aP/(RT)²
B = bP/RT

Cubic equation is solved as per method shown above to get ZV.

Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0

Based on ZV, vapor fugacity φiV is calculated for each data set.

Vapor Fugacity using Peng Robinson EOS

K1 is calculated as:

K1Calculated = φ1L1V

Difference of K1Experimental and K1Calculated is obtained for all data points. An objective function is defined as summation of all these differences. Click on Solver in Data Ribbon (Excel 2010) to open dialog box for Solver parameters and input data as shown below.

Solver for BIP Regression

Minimize the objective function by changing values of kij. Uncheck Make Unconstrained Variables Non-Negative, as k12 can take negative value. Click solve to start regression and new value k12 is calculated.

 k12 = -0.0605

Above value is then used to calculate y1 values and results are plotted to check the deviation.

Plot for Acetone Chloroform XY Data

Determine temperature dependent binary interaction parameter used in Peng Robinson EOS ( k12 = k1 + k2.T + k3/T ) for a mixture of acetone and n-Hexane from P x-y experimental data available at 318.15 °K.

VLE data for Acetone n-Hexane

Follow above steps to calculate pure component and mixture parameters. Change formula for kij to make it temperature dependent. As an initial guess put values for 3 parameters (k1,k2,k3) which gives a numeric difference value. In solver select these 3 variables to be varied for minimizing the difference value and then solve it. After doing regression binary interaction parameters are obtained and result is plotted as following.

k1 = -0.0218
k2 = 3.142x10-4
k3 = 10.155

Plot for Acetone n-Hexane XY data

Note : In case of non-convergence, try different initial guess values for kij and then let solver obtain optimized values.

Spreadsheet for Regressing BIP for Peng Robinson EOS

Simple Multicomponent Distillation

Simple Multicomponent Distillation

This article shows application of Thiele and Geddes method to solve simple multi-component distillation problem in excel spreadsheet based on θ method of convergence.

Antoine equation is used to calculate equilibrium constant Ki.

Ki = exp(Ai - Bi / (T + Ci))/P

where, Ai, Bi & Ci are Antoine equation constants for a component i, T & P are temperature and pressure at a stage.

Feed Analysis

To start with bubble point pressure (PBubble) and dew point pressure (PDew) are determined for feed mixture.

  • P < PDew, Mixture exists as superheated vapor.
  • P > PBubble, Mixture exists as sub-cooled liquid.
  • PDew < P < PBubble, mixture exist in vapor and liquid phase.

For last case, PT Flash calculation is performed to determine liquid (Lf) & vapor (Vf) fraction and molar composition.

Stage-wise Calculation

Liquid and Vapor flow across column are determined as:

L (Above Feed Stage) = R.D
V (Above Feed Stage) = (1+R).D
L (Below Feed Stage) = R.D + Lf
V (Below Feed Stage) = (1+R).D - Vf

Start of Iteration
A Linear temperature gradient is assumed across the column. Ki values are calculated for all component on each stage based on Antoine equation.

θ Method of Convergence
L/D ratios are calculated for all components from top to stage above feed stage.

Top Stage
Total Condenser

LD(1,i) = L(1)/D

Partial Condenser

LD(1,i) = L(1)/(K(1,i).D)

where 1 indicates top stage, i denotes the component and LD denotes L/D ratio.

Above Feed Stage

LD(j,i) = (LD(j-1,i) + 1).L(j)/(K(j,i).V(j))

where j represents the stage and i represents a single component.

Similarly L/B ratios are calculated from bottom upto feed stage

Bottom Reboiler

LB(bottom) = 1

L/B ratios from bottom upto feed tray is calculated as:

LB(j,i) = LB(j+1,i).K(j+1,i).V(j+1)/L(j+1)  + 1 

For Feed Stage

Vfb(i) = V(f).K(f,i)/( L(f).LB(f,i) )

B/D ratio for feed plate is calculated as:

bd(i) = ( LD(nf-1,i) + Lf.Xfi/(F.Zfi))/(Vfb(i) + Vf.Yfi/(F.Zfi))

where bd represent B/D ratio at feed plate.

Based on overall component material balance.

F.Xi = D.Xdi + B.Xbi
d(i) = D.Xdi
b(i) = B.Xbi
F.Xi = d(i) + b(i)
Σd(i) = D

A multiplier θ is defined as:

[b(i)/d(i)]corrected = θ[b(i)/d(i)]calculated

Above equation is rewritten as:

bd(i)co = θ.bd(i)ca

Function g(θ) is defined as:

d(i)co = FXi/(1+θ.bd(i)ca)
g(θ)  = Σd(i)co - D
      = Σi FXi/(1+θ.bd(i)ca) - D

Above equation is solved using Newton Raphson method.

g'(θ) = -Σi bd(i)ca.FXi/[1+θ.bd(i)ca

New estimate of θ is made as following:

θnew = θ - g(θ)/g'(θ)

Iteration is done till there is negligible change in value of θ. Converged value of θ is used to estimate d(i)co. Value of b(i)co is calculated as following:

b(i)co = θ.d(i)

Liquid mol fraction from top to stage above feed stage are calculated:

X(j,i) = LD(j,i)ca.d(i)co / Σi LD(j,i)ca.d(i)co

Liquid mol fraction from feed to bottom tray are calculated:

X(j,i) = LB(j,i)ca.b(i)co / Σi LB(j,i)ca.b(i)co

Bubble point temperature is calculated for each stage. With this first iteration completes. In next iteration new set of Ki values are calculated for each tray based on updated temperature profile and all of the above steps are repeated till there is no change in temperature profile and θ value becomes 1.

Spreadsheet for Simple Multicomponent Distillation


  • Fundamentals of Multicomponent Distillation – C.D. Holland
Shortcut Sizing for Air Cooled Heat Exchanger

Shortcut Sizing for Air Cooled Heat Exchanger

An air-cooled exchanger is used to cool fluids with ambient air. 1 inch OD tube is the most popular diameter and the most common fins are 1/2 inch or 5/8 inch high. Tube configuration used in this guideline : 1 inch OD tube, 5/8 inch fin height, 10 fins per inch and 2.5 inch triangular pitch.

APSF – External area in ft²/ft² of bundle face area.

Rows APSF Face Velocity, ft/min
3 80.4 700
4 107.2 660
5 134.0 625
6 160.8 600

APF – Total external area/ ft of fintube in ft²/ft ~ 5.58.

Face Velocity – Typical air face velocities (VFace ) used in design are tabulated above, these value result in optimum cost of exchanger.

Obtain Process Parameters
Obtain process duty (Q), hot process side inlet (T1) and outlet (T2) temperature. Select an overall heat transfer coefficient (U) from literature based on type of fluids. Select an air inlet temperature (t1) that is not exceeded during a certain percentage of time over the year (e.g. 95% of the time).

Calculate Air Density
Density at air inlet temperature and site elevation

ρo = 14.696 x 29 /(10.7316 x (t1 + 459.67))
ρAir / ρo = exp(-29 x z/ (1545 x (t1 + 459.67)))

t1: Air Inlet Temperature, °F
z : Elevation above sea level, feet

Calculate MTD
Assume an Air outlet temperature (t2) and calculate LMTD.

LMTD = ((T1-t2) -(T2-t1))/ ln((T1-t2)/(T2-t1))
R = (T1 - T2)/(t2 - t1)
S = (t2 - t1)/(T1 - t1)

LMTD Correction factor (F) is estimated based on following graphs

Cross Flow Single Pass

LMTD Correction Factor Cross Flow Single Pass

Cross Flow Two Pass

LMTD Correction Factor Cross Flow Multi Pass

Calculate Air Flowrate
Finned area is estimated using following equation :

AFinned = Q / (U * F * LMTD)

Bundle Face area is estimated using following equation :

AFace = AFinned / APSF

Air Flow is estimated using following relation :

VAir = AFace * VFace

Air Mass flowrate is estimated :

MAir = VAir * ρAir

Air temperature rise is calculated :

ΔT = Q / (MAir * CpAir )

Revised air outlet temperature is calculated :

t2 = t1 + ΔT

This temperature is again used in above steps to re-estimate air outlet temperature. Above steps are iterated till there is no change in air outlet temperature.

Air Cooler Dimension
Air cooler width is calculated :

W = AFace / LTube

Number of Tubes are calculated :

NTube = AFinned / (APF * LTube)

Number of Tubes per Row are calculated :

Nr Tube = NTube / No of Rows

LTube : Length of Tube

Air Side Pressure Drop
Air Side pressure drop is calculated as following :

ΔP Total = ΔPStatic + ΔPVelocity

Static Pressure Drop
Pressure drop across tube bundle.

ΔP Static = FP * No of Rows / DR
FP = 6*10-8 * ( GFace )1.825

DR : ρAir / ρAir at seal level and 70°F
GFace : Air face mass velocity in lb/h.ft² face area
ΔPStatic : Static pressure drop in inch of H2O

Velocity Pressure Drop
Typically 2 fans are used in air cooler. Fan area per fan (FAPF) is calculated as following :

FAPF = 0.4 * AFace / No. of Fans

Fan Diameter is calculated :

D = (4 * FAPF / π )0.5
ΔPVelocity = (ACFM / (4005 * (π* D2/4)) )2* DR

ACFM : Air flowrate in Actual Cubic Feet per Minute
ΔPVelocity : Velocity pressure drop in inch of H2O

Power Calculation
Break power for fan is calculated as following :

BHP = ΔP Total * ACFM / 6356 / ηFan

Motor power is calculated as following :

Power = BHP / ηMotor

Spreadsheet for Shortcut Air Cooled Heat Exchanger Design


  1. Process Heat Transfer: Principles and Applications,2007, Robert W Serth
  2. Handbook of Chemical Engineering Calculations, Nicholas P Chopey
  3. Rules of Thumb for Chemical Engineers, Carl R Branan
  4. GPSA, Engineering Databook, 12th Edition FPS