Exercise Result Checks ¶ This page helps you check whether your own solution of an exercise is heading
in the right direction. It is not a set of solutions. For each exercise it gives
rounded target values or ranges, the trends and signs you should see, sanity
checks that a correct model passes, and pitfalls that often lead to wrong
results. Some exercises also show a reference figure.
How to use it:
Build your model with the data of the exercise first, without looking here.
Compare your results with the targets. Agreement within the stated rounding
or range means you are on track; small differences can come from grid size,
time step or solver tolerances, so refine before you worry.
If a result is far off, go through the sanity checks and pitfalls: units,
the sign convention of the boundary conditions, reading outlet values at the
boundary face rather than in the last cell, and converged solves.
For questions that ask you to explain or discuss, the checks name what a good
answer involves; the reasoning is yours to write.
First order reaction in a batch reactor ¶ With k = 1 s − 1 k = 1~\mathrm{s^{-1}} k = 1 s − 1 Δ t = 0.2 \Delta t = 0.2 Δ t = 0.2 c c c
Halving Δ t \Delta t Δ t
Try a few larger steps (k Δ t k\Delta t k Δ t c c c
With k Δ t = 0.2 k\Delta t = 0.2 k Δ t = 0.2 c c c
Halving Δ t \Delta t Δ t
Sanity check: with a very large step (for example k Δ t = 5 k\Delta t = 5 k Δ t = 5
A good answer places each numerical curve relative to the exact one (above or below) and gives the reason for the difference.
It states the order of accuracy of both schemes and supports it with a refinement of Δ t \Delta t Δ t
It compares the stability of the two schemes as k Δ t k\Delta t k Δ t
With a tight relative tolerance (around 10-8 ) the maximum error against c 0 e − k t c_0 e^{-kt} c 0 e − k t 0 ≤ t ≤ 5 0 \le t \le 5 0 ≤ t ≤ 5 10-7 or smaller.
Tightening rtol and atol should reduce the error; if it does not, check that you compare at the same times as the exact solution.
Pitfall: solve_ivp expects the right-hand side as a function of ( t , c ) (t, c) ( t , c )
Equilibrium consecutive batch reactions ¶ c A c_A c A c B c_B c B c C c_C c C c A , 0 c_{A,0} c A , 0
The sum c A + c B + c C c_A + c_B + c_C c A + c B + c C − 3 ^{-3} − 3
Your analytical solution must contain two negative exponential rates whose sum is − ( k 1 + k − 1 + k 2 ) -(k_1 + k_{-1} + k_2) − ( k 1 + k − 1 + k 2 ) k 1 k 2 k_1 k_2 k 1 k 2
With the stiff rate constants forward Euler blows up unless Δ t \Delta t Δ t 10-12 s.
With Δ t = 0.1 \Delta t = 0.1 Δ t = 0.1 Δ t \Delta t Δ t
The total concentration must be conserved to round-off.
With the stiff rate constants and a large step (for example 0.1 s) the solution must stay bounded and non-negative.
With the non-stiff constants your Newton-based step must reproduce the linear backward Euler result to solver tolerance, and Newton should converge in one or two iterations, since the residual is linear.
With the stiff constants and Δ t = 0.01 \Delta t = 0.01 Δ t = 0.01 t = 2 t = 2 t = 2 c C c_C c C c A c_A c A c B / c A c_B/c_A c B / c A
Pitfall: the Jacobian of the time-step residual contains the accumulation term I / Δ t \mathbf{I}/\Delta t I /Δ t
For the stiff set, implicit methods (Radau, BDF, LSODA) should agree on the final state at t = 2 t = 2 t = 2 c C c_C c C c B / c A c_B/c_A c B / c A k 1 / k − 1 k_1/k_{-1} k 1 / k − 1
Implicit methods should need at most a few thousand function evaluations.
An explicit method such as RK45 needs steps of order 10-13 s; do not run it over the full 2 s, try a very short interval and extrapolate the cost.
Steady-State 1D Fixed-Bed-Reactor Model: First-Order Exothermal Reaction ¶ The inlet total concentration is about 41 to 42 mol m− 3 ^{-3} − 3
The reactor reaches full conversion with an outlet temperature of 443 K in the adiabatic limit; the temperature rises monotonically, so there is no hot spot.
The reaction rate is not highest at the inlet: it peaks inside the bed, near z ≈ 0.2 z \approx 0.2 z ≈ 0.2 X A = 0.5 X_A = 0.5 X A = 0.5
Pitfall: Δ H r \Delta H_r Δ H r
With the molar fluxes as unknowns, conversion proceeds more slowly than with constant velocity: X A = 0.5 X_A = 0.5 X A = 0.5
The outlet velocity is about 6.0 m s− 1 ^{-1} − 1
Pitfall: the rate needs c A = ( F A / F t o t ) P / ( R T ) c_A = (F_A/F_\mathrm{tot})\,P/(RT) c A = ( F A / F tot ) P / ( RT ) F A F_A F A
Your equation for d v / d z \mathrm{d}v/\mathrm{d}z d v / d z
With no reaction and constant temperature your model must reduce to d c i / d z = 0 \mathrm{d}c_i/\mathrm{d}z = 0 d c i / d z = 0 v v v
Pitfall: the velocity belongs inside the derivative, d ( c i v ) / d z \mathrm{d}(c_i v)/\mathrm{d}z d ( c i v ) / d z c i d v / d z c_i\,\mathrm{d}v/\mathrm{d}z c i d v / d z
The results must agree with the molar-flux model of question 2: X A = 0.5 X_A = 0.5 X A = 0.5 − 1 ^{-1} − 1
Sanity check: ∑ i c i \sum_i c_i ∑ i c i P / ( R T ) P/(RT) P / ( RT )
A good answer computes the adiabatic temperature rise from Δ H r \Delta H_r Δ H r
It checks whether T − T i n T - T_\mathrm{in} T − T in
For the velocity, it separates the contribution of the mole-number change from that of thermal expansion, using the ideal gas law.
With about 250 steps your backward Euler profiles should lie close to the solve_ivp ones, with the same outlet temperature (443 K) and outlet velocity (about 6.0 m s− 1 ^{-1} − 1
Refining the step should move the position of X A = 0.5 X_A = 0.5 X A = 0.5 solve_ivp value; backward Euler smears the steep ignition region on coarse steps.
Pitfall: check that Newton converges in every step; use the previous step as initial guess.
Convection of one species using different schemes ¶ Reference case for the numbers on this page: v = 0.01 v = 0.01 v = 0.01 − 1 ^{-1} − 1 L = 1 L = 1 L = 1 c = 0 c = 0 c = 0 c i n = 1 c_\mathrm{in} = 1 c in = 1 N = 100 N = 100 N = 100 t = 50 t = 50 t = 50 x = 0.5 x = 0.5 x = 0.5
Upwind stays between 0 and 1 but smears the front; the L 1 L_1 L 1 N = 100 N = 100 N = 100
For a step profile the error decreases more slowly than first order on refinement (roughly with Δ x 1 / 2 \Delta x^{1/2} Δ x 1/2
Sanity check: until the front reaches the outlet, the amount in the domain must grow as v c i n t v\,c_\mathrm{in}\,t v c in t
Central differencing keeps the front steeper but oscillates behind it; with N = 100 N = 100 N = 100
Refining the grid reduces the L 1 L_1 L 1 N = 100 N = 100 N = 100
Pitfall: write the scheme in flux form so that the boundary cells stay conservative.
The minmod scheme gives the sharpest front of the three and stays between 0 and 1: no new extrema.
Its L 1 L_1 L 1 N = 100 N = 100 N = 100
Pitfall: the ratio r r r
Multicomponent convection-reaction with first-order chemical kinetics ¶ The CFL number v Δ t / Δ x v\,\Delta t/\Delta x v Δ t /Δ x
The largest eigenvalue magnitude of the reaction Jacobian is about 4.7 s− 1 ^{-1} − 1
With all rate constants set to zero you must recover pure convection of a step, bounded between 0 and 1.
Behind the front, c A + c B + c C c_A + c_B + c_C c A + c B + c C c A , i n c_{A,\mathrm{in}} c A , in
Nothing leaves the column before t = τ = L / v t = \tau = L/v t = τ = L / v
Sanity check: the steady outlet must equal the batch solution of exercise 1.2 at t = τ t = \tau t = τ
At v = 1.0 v = 1.0 v = 1.0 − 1 ^{-1} − 1 c A ≈ 0.22 c_A \approx 0.22 c A ≈ 0.22 c B ≈ 0.16 c_B \approx 0.16 c B ≈ 0.16 c C ≈ 0.62 c_C \approx 0.62 c C ≈ 0.62 v = 0.1 v = 0.1 v = 0.1 − 1 ^{-1} − 1
The maximum of c B c_B c B v v v
Behind the front at x = v t x = vt x = v t
At v = 0.3 v = 0.3 v = 0.3 − 1 ^{-1} − 1 t = L / v ≈ 3.3 t = L/v \approx 3.3 t = L / v ≈ 3.3
A tiny residual change (order 10-4 ) between snapshots after that can come from the limiter switching; it is not a mass leak.
Multicomponent convection-reaction with general chemical kinetics ¶ Run both methods on the same grid and time step; then the difference should be at round-off or solver-tolerance level, not at discretization level.
At v = 0.3 v = 0.3 v = 0.3 − 1 ^{-1} − 1 c C c_C c C
A good answer explains why the root-seeking step and the linear solve must agree for first-order kinetics.
The batch fixed point is ( c X , c Y ) = ( c A , c B / c A ) (c_X, c_Y) = (c_A, c_B/c_A) ( c X , c Y ) = ( c A , c B / c A ) c B = 1 + c A 2 c_B = 1 + c_A^2 c B = 1 + c A 2 c A = 1 c_A = 1 c A = 1 c B = 1.7 c_B = 1.7 c B = 1.7 c B = 3 c_B = 3 c B = 3
Sanity check: once the column is filled, the profile must equal the batch trajectory evaluated at t = x / v t = x/v t = x / v
To see oscillations in space, the residence time must exceed the oscillation period (several seconds); at short residence times both cases look alike.
Pitfall: make sure Newton converges in every time step, and keep concentrations non-negative.
Solving a one-component 1D diffusion-reaction equation ¶ Interior rows must be the three-point stencil with weights proportional to D / Δ x 2 D/\Delta x^2 D /Δ x 2
Check: a linear profile c = α + β x c = \alpha + \beta x c = α + β x
The boundary face lies half a cell from the first and last cell centres; using Δ x \Delta x Δ x Δ x / 2 \Delta x/2 Δ x /2
Check: a linear profile that takes the prescribed values at x = 0 x = 0 x = 0 x = L x = L x = L
Assemble the matrix as a sparse array from its three diagonals; it has N × N N \times N N × N 3 N − 2 3N-2 3 N − 2
If you compare with the pymrm operators: pymrm uses a second-order one-sided gradient at the boundary face, so only the first and last rows may differ from the simplest half-cell rows; the interior rows must agree to round-off.
Without reaction and with Dirichlet values on both sides the exact solution is a straight line, which your discretisation must reproduce to round-off (errors of order 10-14 ).
An error much larger than round-off here points to a wrong factor or sign in the boundary rows, not to discretisation error.
Outward normal: at x = 0 x = 0 x = 0 ∂ c / ∂ n = − d c / d x \partial c/\partial n = -\,dc/dx ∂ c / ∂ n = − d c / d x x = L x = L x = L ∂ c / ∂ n = + d c / d x \partial c/\partial n = +\,dc/dx ∂ c / ∂ n = + d c / d x
Check the limits: a = 0 a = 0 a = 0 b = 0 b = 0 b = 0
Without reaction the exact solution is still linear; its two coefficients follow from a 2 × 2 2\times 2 2 × 2
Agreement must again be at round-off level, for any choice of a a a b b b d d d
Pick a test with a ≠ 0 a \neq 0 a = 0
The reaction only adds − k -k − k
With k = 0 k = 0 k = 0
With zero flux on one side and c = 1 c = 1 c = 1
The error is no longer round-off: with N = 50 N = 50 N = 50 10-7 and 10-4 , growing with k k k
Choose k k k L k / D L\sqrt{k/D} L k / D
Refining the grid must reduce the error; if it does not, look for an error in the analytical coefficients or the boundary signs.
Sign check: as Δ t → ∞ \Delta t \to \infty Δ t → ∞
With k = 0 k = 0 k = 0
Evaluate the boundary vector at the new time level, like the other implicit terms.
A convenient test: uniform initial concentration, c = 0 c = 0 c = 0 e − k t e^{-kt} e − k t
At a fixed time the error should roughly halve when you halve Δ t \Delta t Δ t
The series converges slowly at very small times; use enough terms or compare at somewhat later times.
Multi-component 1D counter-diffusion with reaction ¶ Four balances, one per species, each with its own diffusivity; the reaction term has the same magnitude in all four, with a minus sign for A and B and a plus sign for C and D.
Check: adding the balances of A and C (or of B and D) must eliminate the reaction term.
Units: with k 1 k_1 k 1 m 3 k m o l − 1 s − 1 \mathrm{m^3\,kmol^{-1}\,s^{-1}} m 3 kmo l − 1 s − 1 k m o l m − 3 \mathrm{kmol\,m^{-3}} kmol m − 3
With Fick’s law the diffusion matrix is block diagonal per species: no species diffuses because of another species’ gradient.
A good answer shows that the boundary conditions enter only through a constant vector; count how many of its entries are non-zero for the given feed and explain why.
The reaction Jacobian couples only species within the same cell; its sparsity pattern is block diagonal with 4 × 4 4\times 4 4 × 4
Compare your analytical reaction Jacobian with a finite-difference one before relying on Newton.
Keep one array layout throughout: shape ( N , 4 ) (N, 4) ( N , 4 )
Create the numerical Jacobian for the field shape ( N , 4 ) (N, 4) ( N , 4 )
The products C and D peak at about 6.6 k m o l m − 3 \mathrm{kmol\,m^{-3}} kmol m − 3
A and B hardly coexist because the reaction is very fast; the peak position must agree with a simple flux-balance estimate for an instantaneous reaction.
Sanity check: because D A = D C D_A = D_C D A = D C c A + c C c_A + c_C c A + c C k m o l m − 3 \mathrm{kmol\,m^{-3}} kmol m − 3
The diffusion time L 2 / D L^2/D L 2 / D
Check: the long-time limit must coincide with your steady-state solution.
Check Newton convergence in every time step; with a very fast reaction the iterates can overshoot to negative concentrations.
Lowering k 1 k_1 k 1
A good answer relates the position of the reaction plane to the supply of A and B from both sides, and compares the computed positions with a flux-balance estimate.
Multi-component 1D counter-current convection with reaction ¶ One equation per species with its own signed velocity; the system is first order in x x x
Check: adding the balances of A and C (or of B and D) must remove the reaction term.
For each species, the off-diagonal entry in each row must sit on the upwind side; that side is not the same for all four species.
A good answer states at which end each species needs its condition and what the condition at the other end does to the matrix; count the non-zero entries of your boundary vector and explain that number.
The a a a b b b d d d
As Δ t → ∞ \Delta t \to \infty Δ t → ∞
The reaction part of the Jacobian couples only species within the same cell; compare it with a finite-difference Jacobian.
Everything that does not depend on the concentrations can be assembled once, outside the time and Newton loops.
Start each Newton solve from the previous time level and check convergence before accepting the step.
Create the numerical Jacobian for the field shape ( N , 4 ) (N, 4) ( N , 4 )
Make sure the velocity passed to the convective-flux operator has the right sign per species.
From an empty column the fronts need about one residence time L / v = 1 L/v = 1 L / v = 1
At steady state the conversion of A is about 99 %, and the profiles are mirror images: c A ( x ) = c B ( L − x ) c_A(x) = c_B(L - x) c A ( x ) = c B ( L − x )
Conservation check: at steady state c A + c C c_A + c_C c A + c C c B + c D c_B + c_D c B + c D k m o l m − 3 \mathrm{kmol\,m^{-3}} kmol m − 3
A direct steady-state Newton solve must agree with the end of the time march to round-off.
Test at least one high rate constant (for example k 1 = 50 k_1 = 50 k 1 = 50 c = 0 c = 0 c = 0 clip_approach.
Clipping must never change the converged profile, only the path to it.
A good answer relates the difficulty to the Damköhler number k 1 c A , i n L / v k_1 c_{A,\mathrm{in}} L / v k 1 c A , in L / v
Convection-dispersion-reaction in a 1D reactor model ¶ Pass the same boundary-condition tuple to the gradient and to the convective-flux constructor.
Mind the outward normal: at the inlet ∂ c / ∂ n = − ∂ c / ∂ x \partial c/\partial n = -\partial c/\partial x ∂ c / ∂ n = − ∂ c / ∂ x
Check: without reaction the steady state must be c = c 0 c = c_0 c = c 0
Without reaction the problem is linear, so the Jacobian of g \mathbf{g} g
As Δ t → ∞ \Delta t \to \infty Δ t → ∞
Run until t t t L / v L/v L / v L 2 / D a x L^2/D_{ax} L 2 / D a x D a x = 0.1 m 2 / s D_{ax} = 0.1\ \mathrm{m^2/s} D a x = 0.1 m 2 /s v = 1 v = 1 v = 1 L = 1 L = 1 L = 1
The accumulation term must use the solution of the previous time step, not the current Newton iterate.
For D a x = 0.1 m 2 / s D_{ax} = 0.1\ \mathrm{m^2/s} D a x = 0.1 m 2 /s v = 1 v = 1 v = 1 L = 1 L = 1 L = 1 k = 1 s − 1 k = 1\ \mathrm{s^{-1}} k = 1 s − 1 1 − e − 1 = 0.63 1 - e^{-1} = 0.63 1 − e − 1 = 0.63
With N = 100 N = 100 N = 100 10-3 , roughly halving each time you double N N N
A semi-infinite analytical solution ignores the outlet condition, so deviations near the outlet are expected; the finite-domain solution with both boundary conditions allows a sharper test.
Dispersion only: with v = 0 v = 0 v = 0
Convection only: a step front must be centred at x = v t x = vt x = v t N = 100 N = 100 N = 100
Outlet response to an inlet step: the mean residence time must be L / v L/v L / v
With α = 1 \alpha = 1 α = 1
Give the numerical Jacobian a field shape with a trailing component axis, ( N , 1 ) (N, 1) ( N , 1 )
For k = 1 s − 1 k = 1\ \mathrm{s^{-1}} k = 1 s − 1 α = 0.5 \alpha = 0.5 α = 0.5 α = 2 \alpha = 2 α = 2
Cross-check with an independent solver, for example a boundary-value solver for the steady equation: expect agreement to about 10-3 , limited by upwind numerical diffusion.
For α < 1 \alpha < 1 α < 1 k k k k = 4 s − 1 k = 4\ \mathrm{s^{-1}} k = 4 s − 1 α = 0.5 \alpha = 0.5 α = 0.5
Particle model: first order reaction ¶ The spherical geometry enters only through the divergence (face areas ∝ r 2 \propto r^2 ∝ r 2 ∝ r 3 \propto r^3 ∝ r 3
Check: the face at r = 0 r = 0 r = 0
Convert the surface flux into a rate per unit particle volume before comparing it with the intrinsic rate.
Evaluate the surface gradient at the boundary face, including the boundary-condition contribution, not from the last two cell values.
Sanity check: the effectiveness factor (apparent rate divided by k c s k c_s k c s ϕ \phi ϕ ϕ = 1 \phi = 1 ϕ = 1
Two limits: η → 1 \eta \to 1 η → 1 ϕ ≪ 1 \phi \ll 1 ϕ ≪ 1 η ≈ 3 / ϕ \eta \approx 3/\phi η ≈ 3/ ϕ ϕ ≫ 1 \phi \gg 1 ϕ ≫ 1 ϕ = R k / D \phi = R\sqrt{k/D} ϕ = R k / D
If your curve bends away from 3 / ϕ 3/\phi 3/ ϕ ϕ \phi ϕ
At ϕ = 100 \phi = 100 ϕ = 100 η \eta η
Refinement pays only when the reaction layer, about R / ϕ R/\phi R / ϕ ϕ = 1 \phi = 1 ϕ = 1
Strong geometric stretching with many cells can produce cells of zero width and NaN results; check the smallest cell width.
Weisz and Hicks model ¶ With β = 0 \beta = 0 β = 0
Normalise the Arrhenius factor so that it equals 1 at the surface temperature; the steady equations then contain only φ \varphi φ β \beta β γ \gamma γ
The time-dependent form needs one more group, comparing heat and mass diffusivities; it changes the transient but not the steady state.
Storing c ^ \hat c c ^ T ^ \hat T T ^
Sanity check: at steady state T ^ + β c ^ \hat T + \beta \hat c T ^ + β c ^ T ^ \hat T T ^ 1 + β 1 + \beta 1 + β
The time march must end on the steady solution; for φ = 0.4 \varphi = 0.4 φ = 0.4 β = 0.4 \beta = 0.4 β = 0.4 γ = 20 \gamma = 20 γ = 20 t ^ ≈ 1.5 \hat t \approx 1.5 t ^ ≈ 1.5
During ignition a large Δ t \Delta t Δ t
With β = 0 \beta = 0 β = 0 η \eta η φ \varphi φ
For γ = 20 \gamma = 20 γ = 20 β = 0.4 \beta = 0.4 β = 0.4 η \eta η φ = 0.4 \varphi = 0.4 φ = 0.4
For strongly exothermic cases (for example β = 0.8 \beta = 0.8 β = 0.8 φ \varphi φ η \eta η φ = 1 \varphi = 1 φ = 1
A good answer explains η > 1 \eta > 1 η > 1 γ β \gamma\beta γ β
Counter-Current Column Processes ¶ Reference case: m = 20 m = 20 m = 20 U g = U l = 1 U_g = U_l = 1 U g = U l = 1 k o v a = 1 k_\mathrm{ov}a = 1 k ov a = 1 − 1 ^{-1} − 1 L = 1 L = 1 L = 1 c g , i n = 0 c_{g,\mathrm{in}} = 0 c g , in = 0 c l , i n = 1 c_{l,\mathrm{in}} = 1 c l , in = 1 z = L z = L z = L z = 0 z = 0 z = 0
Gas concentration rises monotonically with z z z z = 0 z = 0 z = 0
The liquid enters at z = L z = L z = L z = 0 z = 0 z = 0
Overall balance: gas picked up, U g ( c g , o u t − c g , i n ) U_g(c_{g,\mathrm{out}} - c_{g,\mathrm{in}}) U g ( c g , out − c g , in ) U l ( c l , i n − c l , o u t ) U_l(c_{l,\mathrm{in}} - c_{l,\mathrm{out}}) U l ( c l , in − c l , out )
Your derivation should show that the driving force c g / m − c l c_g/m - c_l c g / m − c l S = m U g / U l S = m U_g / U_l S = m U g / U l
Limit check: for S = 1 S = 1 S = 1 z z z
With 100 cells the maximum difference between numerical and analytical profiles is of order 10-3 , and it should shrink roughly in proportion to the cell size (first-order upwinding).
Rearranged for the NTU, your analytical outlet concentrations must return the imposed NTU (the Kremser form) to many digits.
For every component, inlet flow plus production must equal outlet flow to round-off, provided the outlet flows are the face fluxes of the scheme and not the last cell values.
Setting the Langmuir constant K = 0 K = 0 K = 0
The non-linearity only matters where 1 + K c l 1 + K c_l 1 + K c l
Pitfall: once the source term is non-linear or couples components, its Jacobian must be rebuilt every Newton iteration and must include the cross-component derivatives; Newton should then converge in a handful of iterations.
Reactor Model for Heterogeneous Bubble Columns ¶ Nested hold-ups: the volume fractions per unit column volume come out at about 0.12 for the small bubbles and 0.78 for the slurry.
Sanity check: for an inert tracer the steady profiles must be flat, with both bubble phases at the inlet concentration and the slurry in equilibrium with them, c s = c i n / m c_s = c_\mathrm{in}/m c s = c in / m
Pitfall: a fed phase that also has dispersion needs a Danckwerts inlet (total flux equals the feed), not a fixed inlet value. The slurry has no throughflow and gets zero flux at both ends.
The large dispersion coefficients are numerical: doubling them should hardly change the results.
F ( t ) F(t) F ( t )
The mean residence time grows about linearly with the height: roughly 8 s at 5 m and 48 s at 30 m.
Independent check: the mean ∫ 0 ∞ ( 1 − F ) d t \int_0^\infty (1 - F)\,\mathrm{d}t ∫ 0 ∞ ( 1 − F ) d t
If your mean comes out about half the expected value, look at the inlet condition of the dispersed bubble phase.
At L = 10 L = 10 L = 10 k = 0.1 k = 0.1 k = 0.1 − 1 ^{-1} − 1 k = 1 k = 1 k = 1 − 1 ^{-1} − 1 k = 0 k = 0 k = 0
Conversion rises with k k k
A good discussion compares the rate of gas-slurry transfer with the rate of reaction and identifies which phase acts as a bypass.
Kunii and Levenspiel Model for a ‘Fine Particle’ Fluidized Bed ¶ With the given data u b r ≈ 0.45 u_{br} \approx 0.45 u b r ≈ 0.45 u b ≈ 0.54 u_b \approx 0.54 u b ≈ 0.54 δ ≈ 0.18 \delta \approx 0.18 δ ≈ 0.18
The convective velocities per phase are the superficial ones, δ u b \delta u_b δ u b ( 1 − δ ) u e (1-\delta)u_e ( 1 − δ ) u e u 0 u_0 u 0
Check the transfer term: exchange leaving the bubbles must appear with equal size and opposite sign in the emulsion balance.
Conversion at the heights 0.5, 1 and 2 m is about 0.68, 0.89 and 0.99.
The model is linear with constant coefficients, so an exact solution exists (a matrix exponential); your error against it should roughly halve when you double the number of cells with first-order upwinding.
The emulsion concentration collapses within about a millimetre of the inlet and stays far below the bubble concentration.
Pitfall: base X X X X X X
At a bed height of 1 m, X X X K r K_r K r u 0 u_0 u 0 d b d_b d b
Limit check: for K r → ∞ K_r \to \infty K r → ∞ K r K_r K r
A good answer distinguishes a reaction-limited and an exchange-limited regime and relates the effect of u 0 u_0 u 0 d b d_b d b K b e H / u b K_{be}H/u_b K b e H / u b
Computing Residence Time Distributions ¶ The routine should sum, over all phases, the convective plus dispersive flux at the outlet face; a stagnant phase contributes zero.
At a Neumann outlet the result must equal u u u u u u
Tracer fed minus tracer left, both from face fluxes, must equal the tracer hold-up to round-off, and F → 1 F \to 1 F → 1
Reference: L = 1 L = 1 L = 1 u = 1 u = 1 u = 1 P e = 20 \mathrm{Pe} = 20 Pe = 20 t = 1.2 t = 1.2 t = 1.2 F ≈ 0.76 F \approx 0.76 F ≈ 0.76
Pitfall: the boundary convention uses the outward normal, which points in − x -x − x F F F
The mean residence time ∫ 0 ∞ ( 1 − F ) d t \int_0^\infty (1-F)\,\mathrm{d}t ∫ 0 ∞ ( 1 − F ) d t L / u L/u L / u
The variance should follow the closed-vessel result: close to 2 / P e 2/\mathrm{Pe} 2/ Pe τ 2 \tau^2 τ 2
The erfc solution fits well at large Pe (deviation about 0.06 at Pe = 20) and worse as Pe decreases; at small Pe your curve must approach the CISTR curve.
Upwinding and backward Euler act like extra axial dispersion. At Pe of a few hundred a coarse grid or large time step gives a stable but far too wide F F F
Without exchange the mean residence time is the moving hold-up times L / u L/u L / u L / u L/u L / u
Slow exchange gives early breakthrough plus a long tail; fast exchange must approach a single dispersed phase with the total hold-up.
The variance should grow without bound as the exchange coefficient goes to zero (at fixed total hold-up), and fall towards the single-phase value as it grows.
For any exchange coefficient the mean residence time must equal total hold-up divided by total flow.
Without exchange F F F
Pitfall: the RTD is the total outlet flux over the total feed flux, not the outlet concentration of any single phase.
The Westerterp Wave-Model for Axial Dispersion in Packed Beds ¶ In dimensionless form the only group is 1 / B o = D a x / ( v L ) 1/Bo = D_{ax}/(vL) 1/ B o = D a x / ( vL ) P e = 10 Pe = 10 P e = 10 L / R = 100 L/R = 100 L / R = 100
The mean of E ( θ ) E(\theta) E ( θ )
Pitfall: first-order upwinding adds a numerical dispersion of about Δ ζ / 2 \Delta\zeta/2 Δ ζ /2 1 / B o 1/Bo 1/ B o
Mind the outward-normal sign of the dispersive term in the Danckwerts inlet condition.
Consistency of the given coefficients: ε 1 + ε 2 = 1 \varepsilon_1 + \varepsilon_2 = 1 ε 1 + ε 2 = 1 ε 1 v 1 + ε 2 v 2 = v \varepsilon_1 v_1 + \varepsilon_2 v_2 = v ε 1 v 1 + ε 2 v 2 = v L / v L/v L / v
In the long-tube limit the two-phase model must reproduce the Taylor-Aris coefficient v 2 R 2 / ( 48 D ) v^2R^2/(48D) v 2 R 2 / ( 48 D )
For P e = 10 Pe = 10 P e = 10 L / R = 100 L/R = 100 L / R = 100
Feed each phase in proportion to its flow ε i v i \varepsilon_i v_i ε i v i
Organise the results with 1 / B o 1/Bo 1/ B o N = ( L / R ) / ( 0.2131 P e ) N = (L/R)/(0.2131\,Pe) N = ( L / R ) / ( 0.2131 P e )
For N N N
At fixed L / R L/R L / R 1 / B o 1/Bo 1/ B o P e = 48 ≈ 7 Pe = \sqrt{48} \approx 7 P e = 48 ≈ 7
A good discussion looks at the earliest time tracer reaches the outlet in each model and relates it to the phase velocities v 1 v_1 v 1 v 2 v_2 v 2
Modeling a Desiccant Dryer ¶ Equilibrium loadings: about 0.16 kg/kg for the humid feed and about 0.012 kg/kg for the regeneration air.
The humidity and temperature fronts move far slower than the gas (v / ε v/\varepsilon v / ε Δ t = 0.05 \Delta t = 0.05 Δ t = 0.05
Pitfall: w g ∼ 1 0 − 2 w_g \sim 10^{-2} w g ∼ 1 0 − 2 T ∼ 3 ⋅ 1 0 2 T \sim 3\cdot10^{2} T ∼ 3 ⋅ 1 0 2 w g w_g w g
The accumulation Jacobian is block diagonal with one 2 × 2 2\times2 2 × 2
Starting from a bed at w g = 0.006 w_g = 0.006 w g = 0.006
The humidity breakthrough to the feed value is gradual and happens roughly between 230 and 280 s.
Water and enthalpy balances, closed with the face fluxes, must hold to round-off; the total water uptake should equal the equilibrium capacity, about 1.4 kg per m2 ^2 2
Grid convergence is slow with first-order upwinding; expect smeared fronts on 60 cells.
For regeneration the flow is reversed: negative velocity, inlet condition at x = L x = L x = L
The cyclic state is reached within one or two cycles; check this by comparing the bed state at the end of successive cycles.
At the cyclic state the mean humidity of the dried air is about 0.0019 kg/kg, so close to 90 % of the feed water is removed.
Over one cycle, water removed from the process air must equal water released to the regeneration air, and the net enthalpy stored must be zero.
Axial Convection with Radial Dispersion ¶ Your condition must follow from the angular symmetry of the problem and from the fact that no finite flux can leave a line of zero volume.
Implementation check: in the cylindrical finite-volume divergence the face at r = 0 r = 0 r = 0
Pitfall: pymrm boundary dictionaries use the outward normal, which points in the − r -r − r
You should end up with a linear system d c / d x = A c + b \mathrm{d}\mathbf{c}/\mathrm{d}x = \mathbf{A}\mathbf{c} + \mathbf{b} d c / d x = Ac + b A \mathbf{A} A
Sanity check: the least negative eigenvalue of A \mathbf{A} A v R 2 / D r a d vR^2/D_{rad} v R 2 / D r a d − λ 1 2 ≈ − 5.78 -\lambda_1^2 \approx -5.78 − λ 1 2 ≈ − 5.78 λ 1 = 2.405 \lambda_1 = 2.405 λ 1 = 2.405 J 0 J_0 J 0
Expect a very large ratio between the fastest and the slowest eigenvalue: the system is stiff.
Use a stiff (implicit) method and give it the constant matrix as the Jacobian; an explicit method needs a huge number of steps.
Near the inlet only a thin layer at the wall is depleted; far downstream the radial profile takes the shape of J 0 ( λ 1 r / R ) J_0(\lambda_1 r/R) J 0 ( λ 1 r / R )
In the fully developed region ⟨ c ⟩ \langle c\rangle ⟨ c ⟩ e − 5.8 ξ e^{-5.8\,\xi} e − 5.8 ξ ξ = D r a d x / ( v R 2 ) \xi = D_{rad}x/(vR^2) ξ = D r a d x / ( v R 2 ) ξ = 1.5 \xi = 1.5 ξ = 1.5 1.2 × 1 0 − 4 c i n 1.2\times10^{-4}\,c_{in} 1.2 × 1 0 − 4 c in
Take the gradient at the wall face (for example with compute_boundary_values), not a difference between the last cell centre and the wall value.
Balance check: v d ⟨ c ⟩ / d x = − 2 N w / R v\,\mathrm{d}\langle c\rangle/\mathrm{d}x = -2N_w/R v d ⟨ c ⟩ / d x = − 2 N w / R
Close to the inlet N w N_w N w c i n D r a d v / ( π x ) c_{in}\sqrt{D_{rad}v/(\pi x)} c in D r a d v / ( π x ) x − 1 / 2 x^{-1/2} x − 1/2
Sh starts high near the inlet and levels off to a constant for ξ ≳ 0.3 \xi \gtrsim 0.3 ξ ≳ 0.3 λ 1 2 ≈ 5.78 \lambda_1^2 \approx 5.78 λ 1 2 ≈ 5.78
At small ξ \xi ξ 2 / π ξ 2/\sqrt{\pi\xi} 2/ π ξ
Base Sh on the mean concentration and the wall face flux. Do not compare with the laminar-flow value 3.66, which belongs to a parabolic velocity profile.
Halving the radial cell size should reduce the error in Sh by about a factor four.
Diffusion-Reaction in a Cylindrical Pore ¶ An infinitely fast wall reaction is not a volumetric term: check that your PDE has no reaction term and that the reaction enters elsewhere.
Your equation must contain the cylindrical radial operator 1 r ∂ r ( r ∂ r c ) \frac{1}{r}\partial_r(r\,\partial_r c) r 1 ∂ r ( r ∂ r c )
Compare the diffusion time L 2 / D L^2/D L 2 / D
You need one condition on each of the four edges: pore mouth, dead end, axis and wall.
The wall condition should follow from a surface flux balance with rate constant k s k_s k s k s → ∞ k_s \to \infty k s → ∞
Pitfall: pymrm uses the outward normal, so the sign of a gradient term differs between the lower and upper boundary.
Notice the corner where the mouth meets the wall: the two conditions there are incompatible. Think about what that means for the gradient.
The matrix should have at most five nonzeros per row (a five-point stencil).
Only the mouth condition should give a nonzero constant vector b \mathbf{b} b
Implementation check: your 2D matrix should equal the Kronecker combination of the 1D axial and radial operators to round-off.
The mean concentration in the pore is about 0.065 c 0 0.065\,c_0 0.065 c 0 z = L / 2 z = L/2 z = L /2 1 0 − 3 c 0 10^{-3}\,c_0 1 0 − 3 c 0
Beyond the entrance region the area-averaged concentration decays exponentially with a decay length of about R / λ 1 ≈ 0.42 μ m R/\lambda_1 \approx 0.42~\mu\mathrm{m} R / λ 1 ≈ 0.42 μ m
Balance check with face values: the uptake through the mouth equals the consumption at the wall to round-off, and the flux through the dead end is zero.
Pitfall: if the total uptake at the mouth keeps growing by a similar amount at each grid halving, relate this to the corner of question 2, and use quantities inside the pore for your refinement study.
2D Membrane Fixed Bed Reactor Model ¶ The equilibrium conversion is about 0.50.
It should depend only on K = k f / k b K = k_f/k_b K = k f / k b
Sanity check: the net rate evaluated at your equilibrium composition must be zero.
Without axial dispersion each balance is first order in z z z
Only C C C + r +r + r
Limits your equations must reduce to: for k m = 0 k_m = 0 k m = 0 D e , r D_{e,r} D e , r 2 k m / R 2k_m/R 2 k m / R
A good answer identifies a Biot number k m R / D e , r k_m R/D_{e,r} k m R / D e , r
With k m = 0 k_m = 0 k m = 0 z z z
Balance check: the amount of C C C ⟨ c D ⟩ − ⟨ c C ⟩ \langle c_D\rangle - \langle c_C\rangle ⟨ c D ⟩ − ⟨ c C ⟩
Pitfall: take wall concentrations at the face, not in the last cell, and give the stiff integrator a sparse Jacobian.
The outlet conversion should be about 0.50, with all four outlet concentrations close to 5 m o l m − 3 5~\mathrm{mol\,m^{-3}} 5 mol m − 3
Equilibrium is practically reached about halfway down the tube; the radial profiles stay flat.
The conversion rises above equilibrium with k m k_m k m k m = 0.01 m s − 1 k_m = 0.01~\mathrm{m\,s^{-1}} k m = 0.01 m s − 1 k m ≈ 1 m s − 1 k_m \approx 1~\mathrm{m\,s^{-1}} k m ≈ 1 m s − 1
Upper bound: even a perfect sink for C C C D a / ( 1 + D a ) ≈ 0.83 Da/(1+Da) \approx 0.83 D a / ( 1 + D a ) ≈ 0.83 D a = k f c A , i n L / U 0 Da = k_f c_{A,in}L/U_0 D a = k f c A , in L / U 0
Explaining the saturation involves two resistances in series and the Biot number; describe how the wall value of c C c_C c C k m k_m k m
For the validity at k m = 100 m s − 1 k_m = 100~\mathrm{m\,s^{-1}} k m = 100 m s − 1 D e , r D_{e,r} D e , r
The conversion now saturates around 0.62, so the gain above equilibrium is less than half of that in question 5.
The core should stay close to equilibrium while C C C
A good discussion compares the radial dispersion distance D e , r L / U 0 \sqrt{D_{e,r}L/U_0} D e , r L / U 0
Steady-state 2D fixed bed reactor model: first order exothermal reaction ¶ You should have 2 n r 2n_r 2 n r z z z c A c_A c A T T T
Energy balance with the wall face flux: heat released = sensible heat leaving with the gas + heat conducted to the wall, to within quadrature error.
Limiting case: with an adiabatic wall the 2D model must reproduce the 1D adiabatic model in every radial cell.
Δ T a d = 150 K \Delta T_{ad} = 150~\mathrm{K} Δ T a d = 150 K T = 443 K T = 443~\mathrm{K} T = 443 K c A , i n = P / ( R T i n ) ≈ 41 m o l m − 3 c_{A,in} = P/(RT_{in}) \approx 41~\mathrm{mol\,m^{-3}} c A , in = P / ( R T in ) ≈ 41 mol m − 3
Check that ∑ i ν i C p , i \sum_i \nu_i C_{p,i} ∑ i ν i C p , i
The averaged conversion still shows an S-shaped ignition, reaching 0.5 at about z = 0.20 m z = 0.20~\mathrm{m} z = 0.20 m
The average temperature peaks below 443 K (in the 430 to 435 K range) and then falls to about 410 K at the outlet.
Pitfall: the area average equals the flow average only because v v v ρ C p \rho C_p ρ C p
The exit concentration profile is flat at essentially zero, but the exit temperature profile is not flat and not monotonic: its maximum lies off the axis and it drops to T w T_w T w
The hot spot is about 580 K, off the axis, around z ≈ 0.22 m z \approx 0.22~\mathrm{m} z ≈ 0.22 m
A good explanation compares the radial time scales of mass dispersion (R 2 / D e , r R^2/D_{e,r} R 2 / D e , r ρ C p R 2 / λ e , r \rho C_p R^2/\lambda_{e,r} ρ C p R 2 / λ e , r
Test your explanation: with D e , r D_{e,r} D e , r λ e , r / ( ρ C p ) \lambda_{e,r}/(\rho C_p) λ e , r / ( ρ C p )
A 2D Gas-Solid Fluidized Bed ¶ Only the bubble phase is convected; the cloud and emulsion balances are algebraic at every point and need no boundary conditions.
The bubble balance needs an inlet condition, a symmetry condition at the axis and a zero-flux wall, and no outlet condition.
Limit check: for species A A A K f K_f K f
Summed over A A A B B B C C C
The species balances must close to round-off when you use the axial face fluxes.
The inlet molar flow of A A A π R 2 ⟨ U 0 ⟩ c A , i n ≈ 7.9 m o l s − 1 \pi R^2\langle U_0\rangle c_{A,in} \approx 7.9~\mathrm{mol\,s^{-1}} π R 2 ⟨ U 0 ⟩ c A , in ≈ 7.9 mol s − 1
Pitfall: do not impose an outlet condition that the first-order model does not have; with upwinding the outlet face takes the value of the last cell.
Outlet c A c_A c A c C c_C c C
c B c_B c B r / R ≈ 0.7 r/R \approx 0.7 r / R ≈ 0.7
A good explanation uses the local residence time L / U 0 ( r ) L/U_0(r) L / U 0 ( r )
The emulsion, which holds most catalyst, should be depleted in A A A
X A X_A X A S B S_B S B
Pitfall: weight the outlet concentrations with the local velocity (cup mixing); a plain area average is wrong for a parabolic profile.
With first-order upwind in z z z n z n_z n z 10-3 in X A X_A X A
X A X_A X A S B S_B S B
The first-order decay constant of A A A K f ≈ 0.29 s − 1 K_f \approx 0.29~\mathrm{s^{-1}} K f ≈ 0.29 s − 1
Check: the ODE result must match the closed-form solution of a consecutive first-order reaction, and your 2D code with a uniform velocity must converge to it.
A 2D Bubble Column Model ¶ ⟨ ε G ⟩ = 0.13 \langle\varepsilon_G\rangle = 0.13 ⟨ ε G ⟩ = 0.13 ⟨ ε L ⟩ = 0.87 \langle\varepsilon_L\rangle = 0.87 ⟨ ε L ⟩ = 0.87 r r r
v L , 0 v_{L,0} v L , 0 Δ v G L \Delta v_{GL} Δ v G L
Check by quadrature that your profiles return U L = 1 c m s − 1 U_L = 1~\mathrm{cm\,s^{-1}} U L = 1 cm s − 1 U G = 10 c m s − 1 U_G = 10~\mathrm{cm\,s^{-1}} U G = 10 cm s − 1
An inert tracer fed at c = 1 c = 1 c = 1 c = 1 c = 1 c = 1
Where the axial flux profile changes (end zones), a radial flux from continuity is needed; without it these cells act as sources or sinks and the tracer test fails.
Use a flux (Danckwerts) inlet condition; with large end-zone dispersion a fixed inlet value lets extra material diffuse in.
Limit check: with radially uniform profiles and no end zones you must recover the 1D axial dispersion (Wehner-Wilhelm) conversion.
The liquid conversion is close to the ideal CSTR and slightly above it: about 0.52 at D a = 1 Da = 1 D a = 1 τ L = ⟨ ε L ⟩ H / U L \tau_L = \langle\varepsilon_L\rangle H/U_L τ L = ⟨ ε L ⟩ H / U L
Increasing the end-zone dispersion further should hardly change the result.
A good answer relates the behaviour to an axial Péclet number of the liquid and to the internal circulation.
The gas conversion is close to the ideal PFR and slightly below it: about 0.62 at D a = 1 Da = 1 D a = 1 τ G \tau_G τ G
The deficit below the PFR is largest at intermediate to high conversion and small at very high D a Da D a
A good answer splits the deficit into contributions (axial dispersion, mixed end zones, radial profiles) and compares the radial mixing time with the gas space time.
Open question, so no target values. Check that the uptake per column volume never exceeds k L a c L ∗ k_La\,c_L^{\ast} k L a c L ∗
With a nonlinear rate, scale the unknowns so that the Newton tolerance is meaningful, and confirm convergence with the residual.
The balance over both phases (gas absorbed = liquid consumed + liquid outflow) must close.
For oxygen in water with air, a saturation concentration of about 8 mg/L is a useful sanity check of your Henry coefficient.
Taylor Dispersion ¶ Write the model in dimensionless form: the RTD should depend only on P e = v ˉ R / D m Pe = \bar{v}R/D_m P e = v ˉ R / D m L / R L/R L / R κ = D m L / ( v ˉ R 2 ) \kappa = D_mL/(\bar{v}R^2) κ = D m L / ( v ˉ R 2 )
Check that the discrete volumetric flow ∑ j a j u j \sum_j a_j u_j ∑ j a j u j c = c i n c = c_{in} c = c in
Pitfall: first-order upwinding adds numerical axial dispersion of about u Δ z / 2 u\Delta z/2 u Δ z /2
Weight the outlet concentration with the local velocity (molar flow), not with area; F F F
For a closed vessel the mean residence time ∫ 0 ∞ ( 1 − F ) d θ \int_0^\infty(1-F)\,\mathrm{d}\theta ∫ 0 ∞ ( 1 − F ) d θ
Limit check: with negligible radial diffusion you must recover the laminar-flow RTD F = 1 − 1 / ( 4 θ 2 ) F = 1 - 1/(4\theta^2) F = 1 − 1/ ( 4 θ 2 ) θ ≥ 1 / 2 \theta \ge 1/2 θ ≥ 1/2 F = 0 F = 0 F = 0 θ = 1 / 2 \theta = 1/2 θ = 1/2
For large κ \kappa κ θ = 1 \theta = 1 θ = 1 κ \kappa κ θ = 0.5 \theta = 0.5 θ = 0.5
Compare the models on the same P e Pe P e L / R L/R L / R
Sanity check: in a long tube (κ ≫ 1 \kappa \gg 1 κ ≫ 1
A good answer relates the quality of each 1D model to κ \kappa κ 2 v ˉ 2\bar{v} 2 v ˉ
Ternary Diffusion with Maxwell-Stefan Equations ¶ Signs: H2 _2 2 N 1 < 0 N_1 < 0 N 1 < 0 2 _2 2 N 3 > 0 N_3 > 0 N 3 > 0 N 2 N_2 N 2 x 2 = 0.5 x_2 = 0.5 x 2 = 0.5
Magnitudes: N 1 N_1 N 1 − 1.8 ⋅ 1 0 − 2 -1.8 \cdot 10^{-2} − 1.8 ⋅ 1 0 − 2 N 2 N_2 N 2 8 ⋅ 1 0 − 3 8 \cdot 10^{-3} 8 ⋅ 1 0 − 3 N 3 N_3 N 3 9 ⋅ 1 0 − 3 m o l m − 2 s − 1 9 \cdot 10^{-3}~\mathrm{mol\,m^{-2}\,s^{-1}} 9 ⋅ 1 0 − 3 mol m − 2 s − 1 c t o t c_{tot} c t o t 39 m o l m − 3 39~\mathrm{mol\,m^{-3}} 39 mol m − 3
The N2 _2 2
Sanity check: set D 23 = D 13 D_{23} = D_{13} D 23 = D 13 2 _2 2 c t o t D 13 Δ x 3 / δ c_{tot} D_{13} \Delta x_3/\delta c t o t D 13 Δ x 3 / δ
The linearized fluxes have the same signs as in question 1; H2 _2 2 2 _2 2 2 _2 2
With a constant matrix the profiles are straight lines, so the N2 _2 2
Evaluate the matrix at the arithmetic mean composition of the two ends; in the binary limit D 23 = D 13 D_{23} = D_{13} D 23 = D 13 2 _2 2
Because N t o t = 0 N_{tot} = 0 N t o t = 0
The two eigenvalues of the Fick matrix [ D ] [D] [ D ] 10-5 to 1 0 − 4 m 2 s − 1 10^{-4}~\mathrm{m^2\,s^{-1}} 1 0 − 4 m 2 s − 1
Transform both the end compositions and the fluxes with the same eigenvector matrix, and transform back at the end.
Both ends are Dirichlet conditions (a = 0 a = 0 a = 0 b = 1 b = 1 b = 1
After convergence the flux must be the same on every face, and the profiles must coincide with the shooting profiles, including the N2 _2 2
Grid refinement: the fluxes converge with order about 2 towards the shooting fluxes of question 1.
The flux depends nonlinearly on the compositions, so use Newton with a sparse (block-tridiagonal) Jacobian and check that it converged.
Use A d = π d 2 / 4 A_d = \pi d^2/4 A d = π d 2 /4 d = 2.08 d = 2.08 d = 2.08 + z +z + z
Expected sequence for N2 _2 2 t = 0 t = 0 t = 0 N 2 = 0 N_2 = 0 N 2 = 0
Each species is conserved over the two bulbs to round-off, and both bulbs end at the volume-weighted mean composition (x 1 x_1 x 1
Pseudo-steady state is justified if the capillary volume is small compared with the bulbs; check this ratio.
Dehydrogenation of Ethanol ¶ Bootstrap: the fluxes follow the stoichiometry, N 2 = N 3 = − N 1 N_2 = N_3 = -N_1 N 2 = N 3 = − N 1 N t ≠ 0 N_t \neq 0 N t = 0
The right-hand sides of the three equations must add up to zero, so ∑ i x i = 1 \sum_i x_i = 1 ∑ i x i = 1
Limit check: for a binary mixture with a stagnant second species your equations must reduce to Fick’s law with the Stefan factor 1 / ( 1 − x 1 ) 1/(1 - x_1) 1/ ( 1 − x 1 )
The boundary conditions are the bulk composition at z = 0 z = 0 z = 0 N 1 = k r c t x 1 , δ N_1 = k_r c_t x_{1,\delta} N 1 = k r c t x 1 , δ z = δ z = \delta z = δ
N 1 N_1 N 1 0.86 m o l m − 2 s − 1 0.86~\mathrm{mol\,m^{-2}\,s^{-1}} 0.86 mol m − 2 s − 1 c t c_t c t 22 m o l m − 3 22~\mathrm{mol\,m^{-3}} 22 mol m − 3
The process is strongly film-limited: x 1 , δ x_{1,\delta} x 1 , δ
Upper bound: N 1 N_1 N 1 k r c t x 1 , 0 k_r c_t x_{1,0} k r c t x 1 , 0
Tighten the integration tolerance; the flux should not change in the first few digits.
The exact matrix method solves the same equations as question 2, so the fluxes and surface compositions must agree to the tolerance of your root finder.
Check your matrix-exponential profile by differentiating it numerically and comparing with the Maxwell-Stefan right-hand side.
The high-flux correction matrix must tend to the identity when the flux tends to zero; here it differs clearly from the identity.
Because k 13 = k 23 k_{13} = k_{23} k 13 = k 23 [ Φ ] [\Phi] [ Φ ]
The linearized flux is slightly below the exact one, by less than 1 %.
Keep the drift terms (N t ≠ 0 N_t \neq 0 N t = 0
Scale the unknowns (flux and mole fractions) to order one before calling the nonlinear solver.
The reaction conserves mass, so the mass-average velocity is zero in the film; in this frame the diffusion fluxes equal M i N i M_i N_i M i N i
Take M 1 = M 2 + M 3 M_1 = M_2 + M_3 M 1 = M 2 + M 3
The resulting N 1 N_1 N 1
Compare every method against an exact reference (question 2 or 3), and check that reference first against a case with a closed-form answer, for example all k i j k_{ij} k ij
A good answer relates the errors of the approximate methods to what each one assumes (profile shape, reference frame, composition at which the matrix is evaluated).
Show that your comparison can fail: solving with the wrong bootstrap (equimolar) must give a much larger error than any of the four methods.
Mass Transfer Limitations Using Maxwell-Stefan Equations ¶ Compare the largest reaction rate with the largest film transfer rate k a c A k\,a\,c_A k a c A a = 6 ( 1 − ε ) / d p a = 6(1-\varepsilon)/d_p a = 6 ( 1 − ε ) / d p 1.2 ⋅ 1 0 3 m 2 m − 3 1.2 \cdot 10^{3}~\mathrm{m^2\,m^{-3}} 1.2 ⋅ 1 0 3 m 2 m − 3
At the inlet (293 K) this ratio is of order 0.1; evaluate it also at higher temperatures up to the adiabatic limit (Δ T a d = 150 \Delta T_{ad} = 150 Δ T a d = 150
A good answer relates the conclusion to the activation energy: k r k_r k r
Bootstrap from the stoichiometry: N B = N C = − N A N_B = N_C = -N_A N B = N C = − N A N t = − N A N_t = -N_A N t = − N A
Dilute-A A A A A A 1 / k e f f = x B / k A B + x C / k A C 1/k_{eff} = x_B/k_{AB} + x_C/k_{AC} 1/ k e ff = x B / k A B + x C / k A C
Binary limit (k A B = k A C k_{AB} = k_{AC} k A B = k A C ( 1 + x ˉ A ) (1 + \bar{x}_A) ( 1 + x ˉ A )
The film equations must leave the total concentration unchanged across the film, and the rate basis (per gas, solid or bed volume) must be used consistently.
Expected outlet: conversion about 0.95 and temperature about 435 K, somewhat below the adiabatic 443 K; half conversion near z z z
Near the inlet c A s / c A b c_A^s/c_A^b c A s / c A b
At the surface B B B C C C
Sanity checks: the balance of A A A
With thousandfold coefficients the profiles must fall on the homogeneous ones, with half conversion near z z z
With the real coefficients the front moves downstream and the outlet conversion is no longer complete.
A good explanation treats the regions before and after ignition separately and relates each to the external Damköhler number and its temperature dependence.
In the binary limit your exact film must reproduce the exact log law with Stefan flow.
The exact solution must conserve the total concentration across the film, and agree with a direct numerical integration of the Maxwell-Stefan equations.
Your linearized relation of question 2 should be close to, but not identical with, the exact one for a test composition.
The exact and linearized films should give nearly the same reactor profiles; conversion differences of order 10-4 , far below 1 %.
Compare both films without grid error (or on the same grid), so that the discretisation error does not hide the difference.
A good discussion relates the size of the difference to how large the fluxes are compared with c t k i j c_t k_{ij} c t k ij
An effective diffusivity in the particle is not given; state the value you assume, since the results scale with it.
Check the particle model against the analytical effectiveness factor η = 3 ϕ 2 ( ϕ coth ϕ − 1 ) \eta = \frac{3}{\phi^2}(\phi\coth\phi - 1) η = ϕ 2 3 ( ϕ coth ϕ − 1 )
With a very large D e f f D_{eff} D e ff η → 1 \eta \to 1 η → 1
Adding internal resistance lowers the rate, so ignition moves downstream and the outlet conversion decreases.
Pressure-Velocity Coupling in a 2D Fixed-Bed Reactor ¶ In the adiabatic, constant-pressure reference conversion is complete and T = 443 T = 443 T = 443 T i n + 150 T_\mathrm{in} + 150 T in + 150 X = 1 X = 1 X = 1
The outlet interstitial velocity is about three times the inlet value, about 6.0 m/s.
Sanity check: v o u t / v i n v_\mathrm{out}/v_\mathrm{in} v out / v in ( 1 + X ) T o u t / T i n (1+X)\,T_\mathrm{out}/T_\mathrm{in} ( 1 + X ) T out / T in
Newton should converge in a handful of iterations, and a different initial guess (cold or adiabatic) should lead to the same state.
The total molar balance and the enthalpy balance (reaction heat = convected out + removed at the wall) should close to round-off.
With the solution’s parameters (h = 40 h = 40 h = 40 2 ^2 2
Pitfall: the wall Robin condition uses the outward normal (+ r +r + r T > T c T > T_c T > T c
With h = 0 h = 0 h = 0
The 1D reference has zero pressure drop by construction; compare the 2D pressure drop with a 1D route that also solves Darcy’s law. They should agree to within the axial discretisation error (the gap roughly halves when you halve Δ z \Delta z Δ z
Expect a pressure drop of order 10 kPa: about 12 kPa adiabatic and about 8 kPa with cooling.
Trend: wall cooling lowers the outlet temperature, the outlet velocity and the pressure drop.
Temperature is highest on the axis and falls towards the cooled wall; the axis-to-wall difference is largest near the hot spot (of order 20 K) and small near the inlet.
c A c_A c A
Part of the radial variation of every concentration is gas density (c t = p / ( R T ) c_t = p/(RT) c t = p / ( RT )
In every cell ∑ i c i \sum_i c_i ∑ i c i p / ( R T ) p/(RT) p / ( RT )
The mean outlet superficial velocity should match ( 1 + X ) N i n R T c u p / p o u t (1+X)\,N_\mathrm{in} R T_\mathrm{cup}/p_\mathrm{out} ( 1 + X ) N in R T cup / p out
The inlet pressure is above 1 bar (about 1.1 bar), so the inlet interstitial velocity is somewhat below 2 m/s.
Coupled Batch Reactor and Particle Model ¶ Sphere: use the spherical divergence and symmetry at r = 0 r = 0 r = 0 r = R r = R r = R D ∂ c s / ∂ r + k m c s = k m c f D\,\partial c_s/\partial r + k_m c_s = k_m c_f D ∂ c s / ∂ r + k m c s = k m c f
c f c_f c f c f c_f c f
r a p p r_\mathrm{app} r app
The exercise gives no numbers: choose them so that both the Thiele modulus and the Biot number k m R / D k_m R/D k m R / D
Steady particle at fixed c f c_f c f r a p p r_\mathrm{app} r app − ( 1 − ε b ) k o v c f -(1-\varepsilon_b)k_\mathrm{ov}c_f − ( 1 − ε b ) k ov c f 1 / k o v = 1 / ( η k ) + R / ( 3 k m ) 1/k_\mathrm{ov} = 1/(\eta k) + R/(3k_m) 1/ k ov = 1/ ( η k ) + R / ( 3 k m )
Without reaction, the uptake at fixed c f c_f c f Δ t \Delta t Δ t
The reactor alone with a prescribed first-order rate must give an exponential decay, and its discrete balance must close to round-off.
Make sure your checks can fail: flip the sign of the film term (k m k_m k m
Within a step, one of the two models has to use information from the previous time level; state which one and why.
c f c_f c f
Compare with a particle treated as quasi-steady (rate from the effectiveness factor) to judge whether the particle transient matters for your parameters.
Fluid plus particles plus the amount reacted must stay constant to round-off. If not, the rate you hand to the reactor is not the flux the particle received.
The splitting error is first order: halving Δ t \Delta t Δ t
For first-order kinetics the coupled problem is linear, so a Laplace transform in time gives an independent semi-analytical reference for c f ( t ) c_f(t) c f ( t )
Expect the lagged c f c_f c f
Check the Jacobian of your coupled step against finite differences before trusting the Newton iterations.
For linear kinetics Newton should converge in very few iterations; many iterations point to an inconsistent coupling term.
The implicit result must satisfy the same conservation check as the explicit one (question 4).
The Schur elimination must agree with a monolithic sparse solve of the same step to round-off.
At small Δ t \Delta t Δ t
Increase Δ t \Delta t Δ t
A fast-exchange case (fast diffusion, large k m k_m k m c f ≈ ε b c f , 0 e − ( 1 − ε b ) k t c_f \approx \varepsilon_b c_{f,0}\,e^{-(1-\varepsilon_b)kt} c f ≈ ε b c f , 0 e − ( 1 − ε b ) k t
Reactor-Particle Coupling ¶ The plug-flow outlet must converge to c b ( L ) / c i n = e − D a c_b(L)/c_\mathrm{in} = e^{-\mathrm{Da}} c b ( L ) / c in = e − Da D a = 2 \mathrm{Da} = 2 Da = 2
With first-order upwinding the outlet error should roughly halve when you double n z n_z n z
Pitfall: read the outlet on the outlet face, not in the last cell.
The effectiveness factor from the surface flux must approach η = 3 ( ϕ coth ϕ − 1 ) / ϕ 2 \eta = 3(\phi\coth\phi - 1)/\phi^2 η = 3 ( ϕ coth ϕ − 1 ) / ϕ 2 ϕ = 1 \phi = 1 ϕ = 1 ϕ ≈ 3.2 \phi \approx 3.2 ϕ ≈ 3.2
Use spherical geometry and read the flux on the outer face of the particle; the error should then drop by about a factor 4 per grid doubling.
The explicit iteration should converge (with plug flow it does), in of order ten cycles for D a = 2 \mathrm{Da} = 2 Da = 2
Check the converged state in the full coupled residual, not only by the change between cycles.
The outlet should approach exp ( − η D a ) \exp(-\eta\,\mathrm{Da}) exp ( − η Da )
For first-order kinetics the problem is linear: Newton should need one step plus one to confirm.
The outlet must equal the explicit result to within the stopping tolerance.
Sanity checks: the nested model equals a column with η k \eta k η k v ( c i n − c b ( L ) ) v\,(c_\mathrm{in} - c_b(L)) v ( c in − c b ( L )) r a p p r_\mathrm{app} r app
Both couplings must give the same conversion; only the cost differs.
Expect the number of explicit cycles to grow with η D a \eta\,\mathrm{Da} η Da
A good explanation relates the convergence of the explicit scheme to the strength of the particle-column exchange compared with transport, and considers what axial back-mixing does to how a rate error propagates.
2D Boundary Surface-Reaction Coupling ¶ An A-atom balance (bulk c A + 2 c B c_A + 2c_B c A + 2 c B θ A + 2 θ B \theta_A + 2\theta_B θ A + 2 θ B
Compare your Jacobian with finite differences at a perturbed state.
Near steady state on the 40 × 30 40 \times 30 40 × 30 c A c_A c A c B c_B c B
Pitfalls: read the outlet on the face, not the last cell; the wall Dirichlet value is an unknown, not a constant.
The weighted sum θ A + 2 θ B + θ ∗ \theta_A + 2\theta_B + \theta_\ast θ A + 2 θ B + θ ∗
The unweighted sum should not: it is off by − θ B -\theta_B − θ B
A drift away from one points to rates not all taken at the new time level, or to unconverged Newton steps.
A A A B B B
Along the wall the vacant fraction increases downstream as the wall gas concentration of A A A
Independent check: at steady state the coverages must satisfy r a d s = 2 r d i m = 2 r d e s r_{ads} = 2 r_{dim} = 2 r_{des} r a d s = 2 r d im = 2 r d es
Doubling k a d s k_{ads} k a d s c A c_A c A
At steady state c B , o u t = ( c A , i n − c A , o u t ) / 2 c_{B,\mathrm{out}} = (c_{A,\mathrm{in}} - c_{A,\mathrm{out}})/2 c B , out = ( c A , in − c A , out ) /2 B B B A A A
Repeat on a finer grid to show the change is larger than the discretisation error.
A good description looks at how the wall concentration of A A A
All bulk concentrations and coverages should stay non-negative in every step; the smallest values occur early, where the wall starts empty and the channel starts without B B B
If you see negative values, suspect a Newton overshoot or a non-monotone convection scheme.
A good answer covers where each source enters the discrete equations, at which concentration it is evaluated, and the units (per wall area or per volume, with the specific area 1 / H 1/H 1/ H
It considers the extra transport step to the wall; a numerical test in which the transverse diffusion is made very fast is a convincing way to show it.
Monolithic 2D Membrane Module ¶ One sparse solve should give a backward error at round-off, and all concentrations should lie between the two feed values.
The retentate is depleted along z z z z = 0 z = 0 z = 0
The retentate outlet cup-mixing concentration is about 0.46 on the base grid; the axial (upwind) error is below 1 %.
Cup-mixing retentate outlet about 0.46, permeate outlet about 0.73: the permeate leaves richer than the retentate outlet, which only counter-current flow allows.
Weight with the ring areas 2 π r Δ r 2\pi r\,\Delta r 2 π r Δ r
Pitfall: take outlets on the axial faces; last-cell values differ in the third decimal.
The net axial molar flow (retentate minus permeate) must be the same at every axial face, to round-off.
Retentate loss, membrane law at the face values, and permeate gain must all agree; a little over half of the retentate feed crosses the membrane in the base case.
Pitfall: the permeate side needs the area factor R p e r m / R r e t R_\mathrm{perm}/R_\mathrm{ret} R perm / R ret
The retentate outlet concentration falls monotonically with P P P
At small P P P P P P 1 − Q p e r m / Q r e t = 0.25 1 - Q_\mathrm{perm}/Q_\mathrm{ret} = 0.25 1 − Q perm / Q ret = 0.25
The 2D result should never lie below the 1D plug-flow result (radial diffusion only adds resistance).
A good answer shows what happens to a lagged (segregated) interface iteration as P P P P P P
Relate the amplification of an error in the membrane flux to P P P P → ∞ P \to \infty P → ∞
Reactor-Particle Coupling with Maxwell-Stefan Film Transfer ¶ Newton should converge in a few iterations to a backward error at round-off.
Reactor mole balances per species should close to round-off when you use the outlet face value; the sum of the film differences ∑ i ( c b , i − c g , i ) \sum_i (c_{b,i} - c_{g,i}) ∑ i ( c b , i − c g , i )
Limit check: with A A A
For A A A B B B C C C B B B C C C
The film drop of A A A c b , A / c g , A c_{b,A}/c_{g,A} c b , A / c g , A
The ratio of centre to boundary concentration of A A A
Outlet conversion about 0.86 on 40 axial cells, based on the outlet face value and c A , i n c_{A,\mathrm{in}} c A , in
The minimum of the state vector is positive, of order 0.01 to 0.1 mol/m3 ^3 3
Both residual blocks should be at round-off at convergence and decrease quadratically together during the Newton iterations.
The blocks have different units (mol/m3 ^3 3 3 ^3 3
Increasing D A B \mathcal{D}_{AB} D A B D A B \mathcal{D}_{AB} D A B
A good description separates the local effect on c b , A / c g , A c_{b,A}/c_{g,A} c b , A / c g , A c g , A c_{g,A} c g , A
A good explanation involves the net molar flux out of the particle (one mole in, two out) and how the Maxwell-Stefan equations couple the fluxes through the mole fractions.
Compare with a run using independent Fickian film coefficients; the conversion difference should exceed the grid error.
Check the total concentration across the film in both models.
