When No Chart Fits: Finite Differences

Charts need a uniform starting concentration and simple boundaries. The explicit finite-difference (Schmidt) method handles any starting profile with one spreadsheet formula, provided you respect its stability limit.

Learning objectives

Write the explicit finite-difference equation for an interior node, and the special forms for insulated and convective surfaces.
Choose $Δx$ and $Δt$ so that $M = \frac{Δx^{2}}{D\,Δt}$ satisfies the stability limit.
Build the march in a spreadsheet, and check convergence by refining the grid.
Recognize unstable and non-monotonic results, and explain where they come from.

Background

Real problems often start from a measured, non-uniform profile. Some have a surface condition that changes, or a diffusivity that varies with concentration, and no chart covers those. The fix is to cut the slab into nodes $Δx$ apart and march forward in time steps $Δt$ . A mass balance on each node gives:

c_{n}(t+Δt) = \dfrac{c_{n+1} + (M-2)c_{n} + c_{n-1}}{M} \qquad M = \dfrac{Δx^{2}}{D\,Δt}

With $M = 2$ (the Schmidt method), each new value is just the average of its two neighbours: $c_{n}(t+Δt) = \frac{c_{n+1} + c_{n-1}}{2}$ . In a spreadsheet that is one formula filled right and down. At an insulated face or a plane of symmetry, use a mirror node. A convective surface needs its own node equation, with $N = \frac{k_{c}Δx}{KD}$ , and a stricter limit of $M \ge 2N + 2$ .

The method is stable only for $M \ge 2$ . That is the same as Perry’s limit $\frac{DΔt}{Δx^{2}} \le \frac{1}{2}$ for the explicit (forward Euler) method in Sec. 3, “Parabolic Equations in One Dimension” (pp. 3-46 to 3-47). Below that limit the weight on $c_{n}$ turns negative, and the solution oscillates and can go negative. Implicit and Crank-Nicolson methods remove the limit at the cost of solving simultaneous equations.

Always check a numerical answer: halve $Δx$ , cut $Δt$ by four to keep $M$ the same, and see whether the answer moves. The simulator’s exact series reference shows the error directly.

Simulator

TRANSPORT PHENOMENA | INTERACTIVE SIMULATOR

Finite-Difference Diffusion: Schmidt Stepper

Explicit finite differences | Slab | Stability, monotonicity and convergence

Diffusion profile

Half-slab: x = 0 is the symmetry plane, x = x₁ is the exposed face.

Δx

Δt

M limit

D in Perry units

Plot units

Use Draw initial profile to change the starting concentration with a pointer, or use the sliders for a uniform initial profile.

Current time0 s

X = Dt / x₁²0

FD centre / sealed face1

Error vs exact0%

All values are within the plot axis.

Centre / sealed-face history

Reference: 80-term analytical slab series, uniform initial concentration.

Model settings

Slab geometry

Exposed-face condition

Typed initial profile (comma-separated node values)

Enter Nₛ + 1 values (4 to 40 slices). The slice count follows the length of the list; Δx = x₁/Nₛ. Any consistent concentration unit works; plots switch to dimensional values.

Diffusivity D (m²/s), logarithmic slider

Half-thickness x₁ (m)

Number of slices Nₛ

Schmidt parameter M

Uniform initial c₀ (mol/m³)

Fixed surface cₛ (mol/m³)

Surface node in row 0

Worked solutions differ on this choice. It changes the first few rows of the node table, not the long-time behaviour.

Playback speed (steps/s)

Model assumptions: constant D, one-dimensional slab, no reaction, constant bath and partition factor, and boundary conditions imposed at t ≥ 0. Y = (c − c_eq) / (c₀ − c_eq), where c_eq = cₛ or cₐ. When c₀ = c_eq, Y is undefined; dimensional concentration remains valid.

Governing equations

Fick’s second law: Perry’s 9th ed., Eq. 5-193b (p. 5-47) with no bulk velocity and no reaction ∂c/∂t = D ∂²c/∂x²

Geankoplis, Schmidt explicit interior node cᵢʲ⁺¹ = [cᵢ₊₁ʲ + (M − 2)cᵢʲ + cᵢ₋₁ʲ] / M

Reflected symmetry / sealed boundary c₀ʲ⁺¹ = c₀ʲ + 2(c₁ʲ − c₀ʲ)/M

Convective surface, solid-side units (cₐ = K c_bulk) c𝒇ʲ⁺¹ = [2N cₐ + (M − 2N − 2)c𝒇ʲ + 2c𝒇₋₁ʲ] / M

Grid and dimensionless groups

Δx = a/Nₛ; M = Δx²/(DΔt); Δt = Δx²/(MD)

     X = Dt/a²; n = x/a; N = k_cΔx/(KD); Bi = k_ca/(KD)

Analytical reference: plane-slab series

Y = Σ Cₖ exp(−λₖ²X) cos(λₖn)

     λₖ tan λₖ = Bi; Cₖ = 4 sin λₖ/(2λₖ + sin 2λₖ)

     Fixed face: λₖ = (k − ½)π, k = 1, 2, …

Monotonicity and stability. M ≥ 2N + 2 guarantees monotonicity (non-negative coefficients); violating it is not always divergence. For a fixed surface N = 0, so the criterion is M ≥ 2. Below M = 2 the interior and centre coefficients (M − 2)/M are negative and the shortest-wavelength error is amplified every step, so a fixed-surface run oscillates and eventually grows without bound. With a convective surface, 2 ≤ M < 2N + 2 makes only the surface coefficient negative; the run can stay bounded and still converge.

The uniform-initial-condition formula for Cₖ is used only for uniform profiles. Drawn profiles use numerical Fourier projection onto the same eigenfunctions.

Spreadsheet node table

The last 15 time steps, with row 0 retained until it scrolls out. Columns 0 to Nₛ span the sealed / centre node through the exposed surface. Scroll horizontally to view all nodes.

Self-checks

Benchmark convention: uniform c₀ = 1, equilibrium concentration = 0, x₁ = 0.01 m, D = 10⁻⁹ m²/s, 10 slices, X = 0.5. Expected values come from an independent implementation; every check starts afresh and runs the real solver. None is a homework answer.

Open this simulator on its own page.

Try this

Load the “unstable” preset and step slowly. What happens to neighbouring nodes from one step to the next, and why?
Compare the default grid with “fine-grid” at the same $X$ . How much does the centre value change, and is the coarse grid good enough?
Switch on Draw initial profile and sketch a profile that is high near the surface and low at the centre. What happens at the centre at first, and why can a chart not show it?

Credits and sources

Perry’s Chemical Engineers’ Handbook, 9th ed.: Sec. 3, “Parabolic Equations in One Dimension”, pp. 3-46 to 3-47; Eq. 5-193b, p. 5-47.
Secondary: Geankoplis, 3rd ed., Sec. 7.7 (numerical methods for unsteady diffusion) and Sec. 5.4 (heat analog).
Simulator: Dr. Ray’s Learning Center. Drafted with AI assistance, then checked against the exact slab series.

Course: Transport Phenomena. Previous: Multidimensional Bodies: the Product Rule.