HPLC Simulator

A packed column is not solid all the way through. Liquid fills the gaps between the particles and the small pores inside them. These spaces determine how much liquid the column holds and how quickly an unretained compound can pass through it.

Porosity means the fraction of a volume that liquid can occupy. Only pores accessible to the mobile phase are counted. Interparticle porosity, ε_e, is the fraction of the empty-column volume between particles. Intraparticle porosity, ε_i, is the fraction of each particle’s volume accessible to mobile phase. Total porosity, ε_T, is the fraction of the column volume accessible to mobile phase:²

ε_T = ε_e + ε_i(1 − ε_e)

For a column with internal diameter d and length L, the cross-sectional area O, mobile-phase hold-up volume V₀, and hold-up time t₀ are:

O = πd²/4
V₀ = O L ε_T
t₀ = V₀/F

F is volumetric flow rate. Hold-up time is the time required for an unretained solute to pass through the column; any post-column tubing delay is calculated separately. Use consistent volume and time units when evaluating these equations.

The same flow rate can be expressed as several linear velocities:

u_s = F/O
u_e = u_s/ε_e
u = u_s/ε_T = L/t₀

Superficial velocity, u_s, is the speed the liquid would have if the column were an empty tube.² Interstitial velocity, u_e, refers to the space between particles. Chromatographic velocity, u, accounts for all accessible mobile-phase volume and describes an unretained solute.

Viscosity describes how resistant a liquid is to flowing. A more viscous liquid requires more pressure to push through the column. The liquid’s temperature and the proportions of water and organic solvent both affect its viscosity. “Eluent” is another name for the mobile phase.

Viscosity, η, determines resistance to flow and influences solute diffusion. The simulator uses equations fitted to measurements of acetonitrile–water and methanol–water mixtures. Acetonitrile or methanol is solvent B; water is solvent A. Let φ be the volume fraction of solvent B, from 0 to 1, and let T_K = T + 273.15, where T is temperature in °C and T_K is absolute temperature in kelvin. Both expressions return viscosity in cP.

Acetonitrile–water³

ln(η) = φ(−3.476 + 726/T_K)
+ (1 − φ)(−5.414 + 1566/T_K)
+ φ(1 − φ)(−1.762 + 929/T_K)

Methanol–water⁴

ln(η) = φ(−4.597 + 1211/T_K)
+ (1 − φ)(−5.961 + 1736/T_K)
+ φ(1 − φ)(−6.215 + 2809/T_K)

Viscosity is the exponential of the corresponding expression. Mixing water with an organic solvent can produce a viscosity that is quite different from the average of the two liquids. These equations account for that behavior.

In a gradient run, solvent composition changes over time. The simulator uses the starting composition to calculate column efficiency and spreading in the tubing. The optional viscosity and backpressure plots follow the changing liquid composition entering the column.

The pump must generate pressure to move liquid through the packed column and narrow tubing. This pressure is called backpressure. Faster flow, a longer column, smaller particles, or more viscous liquid generally increase it.

The pressure needed to push liquid through the column, ΔP_col, is estimated from its resistance to flow. The equation combines Darcy’s law with the Kozeny–Carman relationship for a packed bed:^5,6,7

ΔP_col = 180η L u_s(1 − ε_e)² / (ε_e³d_p²)

d_p is particle diameter. In SI units, viscosity is in Pa·s, lengths in meters, and velocity in m/s; the result is in pascals. The simulator converts this to bar (1 bar = 100,000 Pa). At fixed flow and column dimensions, smaller particles increase pressure approximately as 1/d_p².

Tubing between the column and detector adds resistance too. For a straight, round tube, the Hagen–Poiseuille equation estimates the pressure needed:⁵

ΔP_tube = 8η L_tube F/(πr⁴)
ΔP_total = ΔP_col + ΔP_tube

Here ΔP_tube is the tubing pressure drop, ΔP_total is total backpressure, L_tube is tubing length, and r is the tubing’s internal radius. The symbols η, L, u_s, ε_e, and F retain their definitions above: viscosity, column length, superficial velocity, interparticle porosity, and volumetric flow rate. The strong fourth-power dependence explains why narrow tubing can add substantial resistance. The calculation assumes smooth flow in layers (laminar flow). It does not add separate resistance from fittings or the small liquid chamber inside a detector. For gradients, the pressure overlay scales the initial pressure by the ratio of current inlet viscosity to initial viscosity; it does not separately calculate viscosity at every position inside the column.

Even without flow, dissolved molecules spread out through random motion. This is diffusion. Inside a column, diffusion and differences in the paths molecules take cause a narrow sample band to spread into a wider peak. Column efficiency describes how well the column limits that spreading.

An analyte is a compound being separated. Its diffusion coefficient measures how quickly it spreads through the liquid. The simulator estimates one representative value for the sample using the Wilke–Chang relationship:⁸

D_e = 7.4 × 10⁻⁸T_K√(χ M) / (η V_m^0.6)

D_e is the representative analyte diffusion coefficient in cm²/s, η is in cP, M is the solvent molecular weight in g/mol, and V_m is the analyte molar volume in cm³/mol. Molar volume is the volume occupied by one mole of a substance. The simulator averages the selected compounds’ molar volumes, giving each compound equal weight regardless of concentration. It uses that average to estimate a single diffusion coefficient for the sample.

χ = 2.6(1 − φ) + 1.9φ
M = 18(1 − φ) + M_org φ

The solvent association parameter χ is an empirical correction for interactions between solvent molecules. Its value changes with the fraction of water and organic solvent. The organic modifier’s molecular weight, M_org, is 41 for acetonitrile and 32 for methanol. Compound molar volumes are supplied in the library or entered for custom compounds.

Reduced velocity, ν, compares transport by flowing liquid with spreading by diffusion. It has no units:⁹

ν = u_e d_p/D_e

Here u_e, d_p, and D_e must use consistent length units. Reduced plate height follows the reduced van Deemter equation:¹⁰

h = A + B/ν + C ν
H = h d_p
N = L/H

A, B, and C control three sources of spreading: molecules taking different paths around particles (eddy dispersion), diffusion along the column (longitudinal diffusion), and the time needed to move between the mobile and stationary phases (resistance to mass transfer). h is the dimensionless reduced plate height; H is the height equivalent to a theoretical plate (HETP); N is the theoretical plate count. A “theoretical plate” is a way of describing separation efficiency; there are no physical plates inside the column. Smaller plate height means more theoretical plates in a given column and narrower calculated peaks. At very low flow, diffusion has more time to spread a band; at very high flow, molecules have less time to move between phases. The equation describes this trade-off. All compounds share this efficiency estimate, which uses the initial composition in gradient mode.

In an isocratic run, the percentage of solvent B stays constant. Compounds that spend more time interacting with the stationary phase move more slowly and appear later in the chromatogram. Retention time is the time from injection to the peak’s center at the detector.

The dimensionless retention factor, k′, is the ratio of time a compound spends retained by the stationary phase to time it spends in the mobile phase. The simulator estimates how solvent B changes retention using a straight-line relationship between the logarithm of retention factor and the fraction of solvent B:¹⁴

log₁₀(k′) = log₁₀(k′_w) − S φ

k′_w is the estimated retention factor in pure water, and S describes how strongly retention changes as solvent B is added. Temperature also affects retention. The model describes its effect with straight-line fits to measured data.¹¹ In the notation used by the custom-compound controls:

log₁₀(k′_w) = m_w T + b_w
−S = m_φ T + b_φ
log₁₀(k′) = m_w T + b_w + (m_φ T + b_φ)φ

T is in °C here, not kelvin. The coefficients m_w and b_w are the slope and intercept of log₁₀(k′_w) versus temperature. The coefficients m_φ and b_φ are the slope and intercept of −S versus temperature. The slopes have units of °C⁻¹; the intercepts are dimensionless. All four depend on compound, stationary phase, and organic solvent. The last two coefficients describe −S, so they are added in the final expression. These fitted relationships help explain trends. Predictions far from the conditions used to measure the coefficients may differ from a real experiment.

t_R,col = t₀(1 + k′)
t_R = t_R,col + t_tube

Here t_R,col is the time from injection to the column outlet, t_tube is the post-column tubing transit time, and t_R is the retention time at the detector. Increasing k′ makes retention longer. The selected-compound position overlay follows z = L min[1, t/(t₀(1 + k′))]; where z is distance from the column inlet and t is elapsed time. The position reaches the column outlet before any post-column delay.

In a gradient run, the percentage of solvent B changes during the separation. A programmed change does not reach the column immediately: liquid must first travel through tubing and mixing spaces. The volume between the point where solvents are combined and the column inlet is called dwell volume.

Between two rows in the gradient table, the programmed percentage of solvent B changes at a constant rate. Before time zero the first composition is used; after the final row the last composition is held. The composition arriving at the column differs from this program because of non-mixing and mixing volumes upstream.

t_delay = V_nonmix/F
φ_in(t) = φ_program(t − t_delay)

Here V_nonmix is the upstream non-mixing volume, t_delay is its transport delay, φ_program(t) is the programmed solvent-B fraction, and φ_in(t) is the delayed fraction entering the mixer. Mixing volume acts like a small stirred tank filled with the starting mixture. New liquid mixes with liquid already in the tank, so a sharp change becomes more gradual. For a short time increment Δt, a volume FΔt enters and leaves the mixer. The new composition is calculated by keeping track of how much solvent B enters and leaves:¹³

φ_mix,new = φ_mix,old + α(φ_in − φ_mix,old)
α = min[1, FΔt/V_mix]

The mixing-chamber volume is V_mix. The fractions φ_mix,old and φ_mix,new describe the chamber before and after a time step. The replacement fraction α is the pumped volume divided by mixer volume, capped at 1. The resulting fraction φ_mix(t) is the column-inlet composition as a function of time.

If mixing volume is zero, the delayed incoming composition is used directly. The update is also limited to complete replacement when a time step pumps at least one mixer volume. Mixing rounds and delays changes in composition, while non-mixing volume simply shifts the program in time. Their sum is the nominal dwell volume.

The calculation horizon is the time span over which the simulator calculates the gradient and follows compounds through the column. The numerical mixer time step is this horizon divided by 999. These small calculation steps are separate from the detector samples used to draw the chromatogram. The solvent-B overlay shows the resulting composition at the column inlet.

During a gradient, a compound can speed up as a stronger mobile phase reaches it. A single, fixed retention factor is therefore not enough. The simulator follows the compound through the column, updating its speed as the surrounding solvent changes.

During a gradient, a compound’s retention factor changes as the solvent surrounding it changes. The simulator divides migration into small steps and recalculates k′ using the same composition and temperature relationship as in isocratic mode.

The solvent at a compound’s position entered the column earlier. To account for this travel time, the model uses a corrected time, θ: elapsed time minus the time an unretained solute would take to reach that position. For a uniform column, the equation below adds up the compound’s progress. The compound leaves when the accumulated fraction of the column traveled reaches 1:^14,15

∫₀^{t_R,col − t₀} dθ / [t₀ k′(φ_mix(θ), T)] = 1

At each step the model evaluates the inlet composition at the compound’s solvent-time coordinate. With x = z/L as fractional column position (0 at the inlet and 1 at the outlet), Δx as its increment, Δt as the elapsed-time increment, and Δθ as the solvent-time integration step, migration is accumulated as:

Δx = Δθ/(t₀ k′)
Δt = Δθ + t₀Δx

The simulation tracks elapsed time, retention factor, and position until x reaches 1. If a step carries the compound past the outlet, the simulator estimates where within that step it crossed the outlet. It then adds the travel time through the post-column tubing. The step Δθ is the calculation horizon divided by 1000. In automatic-duration mode, the horizon is extended as needed, up to the acquisition limit. A compound that has not exited within the available horizon is reported as not eluted.

The optional plots of the selected compound’s retention factor and position use the recorded progress through the column, estimating values between calculation steps. The position is reported in millimeters along the column.

After leaving the column, compounds still have to travel to the detector. Connecting tubing makes peaks arrive later and can spread them out. Longer or wider tubing holds more liquid, while very narrow tubing can add substantial backpressure.

A tube of length L_tube and internal radius r adds a volume and transport time:

V_tube = πr²L_tube
t_tube = V_tube/F
u_tube = F/(πr²)

Here V_tube is internal tubing volume, t_tube is tubing transit time, and u_tube is mean flow velocity in the tube.

In smooth flow through a round tube, liquid at the center moves faster than liquid near the wall. Molecules also diffuse between these regions. The model estimates how these effects spread a sample band along the tube.¹² It then converts spreading in distance into spreading in arrival time. Variance is the square of a standard deviation:

σ_z² = 2D_e L_tube/u_tube + r²L_tube u_tube/(24D_e)
σ_tube² = σ_z²/u_tube²

The quantity σ_z is the standard deviation of axial position caused by tubing dispersion, and σ_tube is the corresponding standard deviation in arrival time. Their squares are spatial variance (length²) and time variance (time²), respectively. Use consistent length and time units. This is a simplified laminar-tube dispersion estimate, rather than a detailed model of fittings, bends, or a detector cell. Tubing pressure is included separately in total backpressure. With zero tubing length, the tubing delay, dispersion, and pressure contributions are all zero.

A real sample does not leave the column all at once. Its molecules arrive over a range of times, producing a peak with measurable width. In this simulator, a larger standard deviation, σ, means a wider peak.

Each compound is represented by a Gaussian peak with time standard deviation σ. Spreading from different sources is combined by adding their variances (the squares of their standard deviations). In isocratic mode:

σ² = t_R²/N + τ² + σ_tube² + (V_inj/F)²/12

From left to right, the terms account for spreading in the column, the detector’s response time, spreading in the tubing, and the time taken to inject the sample. Injection is treated as delivery at a constant concentration for a fixed time—a rectangular pulse.¹² τ is the detector time constant, and V_inj is injection volume; V_inj/F is injection duration. The current implementation uses detected retention time, including tubing delay, in the first term. The detector time constant describes how quickly the detector responds. Here its effect is represented as extra peak width rather than a detailed simulation of the detector electronics.

In gradient mode, the width depends on retention factor near the column exit, k′_e:¹⁶

σ² = [t₀(1 + k′_e)]²/N + τ² + σ_tube²

This approximation can produce much narrower peaks than using total gradient retention time in the column term. Injection-volume broadening is omitted in gradient mode. A compound that is strongly retained in the starting solvent can collect in a narrow band at the column entrance before moving onward. This is called focusing. Omitting injection broadening is a useful approximation in that case, but it can make peaks too narrow for compounds that are weakly retained at the start. A gradient can also make the back of a band move faster than its front, squeezing the band narrower. The model does not include a separate correction for this band-compression effect.¹⁷

The table reports σ in seconds. For a Gaussian peak, the width measured halfway up the peak (full width at half maximum) is approximately 2.355σ; the conventional baseline width is approximately 4σ. The model produces symmetric peaks. It does not simulate a stretched trailing edge (tailing), a stretched leading edge (fronting), or changes in peak shape caused by injecting more sample than the column can handle (overload).

Once the arrival time and width are known, the simulator draws a smooth, bell-shaped peak for each compound. Peak height depends on how much compound was injected and how widely it has spread. The individual peaks add together wherever they overlap.

The simulated detector measures concentration: it gives the same signal for every compound at the same concentration in moles per liter. This is a simplifying assumption. Real detectors generally respond differently to different compounds: a UV detector gives a stronger signal for compounds that absorb more strongly at the selected wavelength, while a mass spectrometer’s response depends on ionization and detection efficiency. Readily ionized or already ionic compounds can give stronger responses, depending on measurement conditions.

For compound i, injected amount n_i is sample concentration c_sample,i multiplied by injection volume V_inj. The subscript i identifies a compound in the sample:

n_i = c_sample,i V_inj

Each compound’s bell-shaped peak is described by a Gaussian equation:¹

c_i(t) = n_i/[F σ_i√(2π)]
× exp[−(t − t_R,i)²/(2σ_i²)]

Here c_i(t) is compound i’s concentration at the detector at time t, t_R,i is its detected retention time, and σ_i is its peak’s time standard deviation. With amount in moles and flow in liters per second, this gives concentration in mol/L. The display converts to an appropriate concentration unit. Integrated peak area is n_i/F; for a fixed amount and flow, a narrower peak is taller. Overlapping peaks add together.

c_total(t) = ∑_i c_i(t) + c_offset + ε_noise(t)

Here c_total(t) is the total detector signal expressed as concentration, ∑_i means summation over all sample compounds, c_offset is the constant baseline offset, and ε_noise(t) is the zero-mean random detector-noise contribution at time t. The offset shifts the entire baseline. The noise is random variation around zero, with small deviations more common than large ones (a Gaussian distribution). A fresh noise pattern is generated whenever the chromatogram is recalculated. In the implementation’s numerical units, the noise standard deviation in µM is the Noise control value divided by 1000√(τ), with τ entered in seconds. Thus, increasing the detector time constant broadens peaks and reduces noise. Zooming, panning, and display-only changes reuse the stored samples and do not generate new noise.

Resolution describes how clearly two peaks are separated. Moving their centers farther apart increases resolution; making the peaks wider decreases it. Run duration determines how much of the separation is recorded, and sampling rate determines how frequently the detector records a value.

Resolution compares the separation in retention time with the widths of a pair of peaks:¹⁴

R_s = |t_R,2 − t_R,1| / [2(σ₁ + σ₂)]

Here R_s is dimensionless resolution; subscripts 1 and 2 identify the two compounds being compared. The times at the centers of their peaks are t_R,1 and t_R,2, and their peak standard deviations are σ₁ and σ₂. Vertical bars denote absolute value.

For each eluted compound, the results table compares it with the preceding and following compounds in retention order and reports the smaller resolution. This is the more difficult adjacent separation, which need not be the smallest retention-time gap when widths differ. With no neighboring eluted compound, resolution is undefined. Minimum resolution below the graph is the smallest adjacent-pair value among compounds whose retention times lie within the acquired run.

Automatic run duration uses the last-eluting compound:

t_{end,requested} = 1.05(t_R,last + 3σ_last)

Here t_{end,requested} is the requested time to stop recording, t_R,last is the last-eluting compound’s detected retention time, and σ_last is its peak standard deviation. The three-standard-deviation allowance includes most of the final peak’s trailing tail; the additional 5% leaves space afterward. Manual duration can instead be set in Time range & sampling.

At sampling rate f_s in samples per second, detector samples are taken at t_j = j/f_s. Here j is the integer sample index (0, 1, 2, …), and t_j is the time of that sample in seconds. The endpoint is rounded upward to a sampling interval. A maximum of 100,000 points is allowed, giving a maximum acquired time of 99,999/f_s seconds. If the requested run exceeds that limit, recording stops early and a warning is shown. A higher sampling rate gives more points per second but uses up the point limit sooner. Sampling covers the full acquired chromatogram; zooming changes only the viewed window.

Use a custom compound when you want to supply your own retention behavior. Random unknowns provide invented compounds for exercises, so students can explore a separation without relying on a familiar compound name.

Custom compounds use the same retention equations as library compounds. Enter a molar volume and four retention coefficients for each organic solvent. Molar volume affects the sample’s representative diffusion coefficient; the coefficients determine retention as composition and temperature change.

Random unknowns are synthetic teaching compounds, not measured chemical identities. Their retention coefficients are generated together so that they follow related trends, with some random variation. Separate sets of coefficients describe their behavior in acetonitrile and methanol. Molar volume is drawn between 90 and 190 cm³/mol, and concentration is chosen between 1 and 100 µM on a logarithmic scale, giving equal chance to the ranges 1–10 and 10–100 µM. Generated numerical parameters are rounded to three significant figures. Values explicitly entered for a custom compound retain the supplied precision.