From interface to equilibrium
A closed diffusion couple with composition-dependent mobility.
Concentration profile
Dashed vertical line: original interface. Values are cell averages.
Final local diffusivity
Zero diffusivities are omitted from the logarithmic axis.
Eight independent runs over ±200 K, clipped to 300–2500 K.
Simulator documentation
This simulator explores one-dimensional cation diffusion between a source layer and a surrounding matrix. You can compare concentration profiles, local diffusivities, time evolution and temperature sensitivity.
The calculations use a phenomenological, composition-dependent form of Fick’s second law. The model is intended for exploring diffusion behaviour; quantitative predictions require appropriate material parameters and numerical convergence checks.
1. Defining the diffusion couple
The domain extends from x = 0 to x = L. Initially, the source occupies 0 ≤ x < a, where a is the source length. The remaining domain initially contains only the matrix.
- Temperature (K): The uniform annealing temperature, between 300 and 2500 K. Temperature remains constant during each individual simulation.
- Annealing time (h): The duration of diffusion. Enter zero to inspect the initial state. The solver converts hours to seconds internally.
- Source length (µm): The initial thickness of the cation-containing layer. It must be positive and smaller than the total length.
- Total length (µm): The length of the closed simulation domain. Lengths are converted to metres for the calculations.
- Grid cells: The number of equal-width computational cells, from 20 to 400. More cells resolve sharper spatial features but generally require more time steps.
Both ends are impermeable: no species enters or leaves the domain. The source is therefore a finite reservoir; its concentration is allowed to decrease during diffusion.
2. Source composition and matrix balance
Add up to four distinct cations. Each source concentration is expressed as a percentage on the same composition basis. The sum of the entered cation concentrations must not exceed 100%.
Cᵢ(x, 0) = 0 outside the source
Cmatrix(x, t) = 100 − Σᵢ Cᵢ(x, t)
For example, a source containing 25% Ni and 15% Co initially contains 60% matrix. Outside the source, the initial matrix concentration is 100%.
Concentrations are stored as cell averages. If the source interface crosses a cell, that cell receives the corresponding fractional source loading. This preserves the entered source inventory without rounding the source length to a grid boundary.
Cation parameters
- Cation preset: Selecting an element fills its D₀, Q and β fields. These values remain editable.
- D₀ (m²/s): The Arrhenius pre-exponential factor. Setting it to zero makes that cation immobile.
- Q (eV): The activation energy used in the thermal mobility factor.
- Source concentration (%): The initial concentration of that cation inside the source layer.
- β: A dimensionless empirical parameter controlling how mobility changes with the cation’s own local concentration.
- Colour: The colour used to identify the cation in the plots. It does not affect the calculation.
Preset values are illustrative, unreferenced estimates. An element name alone does not determine its diffusivity. Use parameters appropriate to the host material, phase, composition and temperature range being studied.
Matrix reference parameters
The matrix name labels the balance species. Its D₀ and Q define an Arrhenius reference diffusivity:
This reference appears in the diffusivity plot and exported data. It does not drive the cation profiles. The matrix concentration is calculated by composition balance, rather than by solving a separate matrix diffusion equation.
3. Composition-dependent diffusivity
For each cation i, the solver evaluates a local diffusion coefficient from temperature, its own concentration and the optional blocking factor:
kB = 8.617333262145 × 10⁻⁵ eV/K
Thermal activation
For positive Q, increasing temperature increases the Arrhenius diffusivity. At a fixed temperature and D₀, increasing Q decreases it. Q is entered in eV, so the solver uses Boltzmann’s constant in eV/K.
Concentration sensitivity
- β > 0: The self-concentration factor increases with cation concentration.
- β < 0: The self-concentration factor decreases with cation concentration.
- β = 0: The self-concentration factor equals one. Diffusivity is spatially constant at fixed temperature only when blocking is also disabled.
β is an empirical mobility parameter. The exponential factor is not a thermodynamic factor derived from an activity or free-energy model. When blocking is enabled, the combined diffusivity also depends on the matrix fraction.
Phenomenological site blocking
Blocking disabled: φ = 1
With blocking enabled, decreasing the matrix fraction reduces cation mobility. For example, a matrix concentration of 60% gives φ = 0.60. Below 1% matrix, φ remains fixed at 0.01.
Matrix fraction is not vacancy concentration. This factor is a simplified description of composition-dependent blocking, not a calculated vacancy population. Its 1% floor means blocking alone never completely stops diffusion.
4. Transport equation and boundary conditions
Each cation evolves according to Fick’s second law in conservative form:
∂Cᵢ/∂t = −∂Jᵢ/∂x = ∂/∂x [Dᵢ ∂Cᵢ/∂x]
Jᵢ(0, t) = Jᵢ(L, t) = 0
The negative sign makes each diffusive flux point down that species’ concentration gradient. Since D may vary with position and composition, the solver retains it inside the spatial derivative.
Here C is expressed in percentage units. The numerical flux therefore has units of percentage points × m/s; it is not directly a molar flux. Converting it to mol/(m²·s) would require a specified concentration basis and an appropriate molar or site density.
5. Conservative numerical method
The domain is divided into N equal-width finite-volume cells. Positions reported in plots and downloads are the cell centres:
xⱼ = (j + ½) Δx, for j = 0, …, N − 1
At every internal face, diffusivity is the harmonic mean of the adjacent cell values. If either adjacent diffusivity is zero, the face diffusivity is zero.
Jᵢ,ⱼ₊½ = −Dᵢ,ⱼ₊½ (Cᵢ,ⱼ₊₁ − Cᵢ,ⱼ) / Δx
Cᵢ,ⱼ,new = Cᵢ,ⱼ + (Δt / Δx) (Jᵢ,ⱼ₋½ − Jᵢ,ⱼ₊½)
A face flux is shared by its two neighbouring cells: the amount leaving one enters the other. Exterior face fluxes are exactly zero. Consequently, each cation’s integrated concentration is conserved to floating-point precision.
Adaptive explicit time stepping
The solver recalculates local diffusivities before each step. It chooses the smallest of the stability limit, the time-resolution limit and the time remaining to the next output:
Dmax is the largest current cation diffusivity. When all cation diffusivities are zero, the stability restriction is unnecessary. A zero-duration run returns the initial state without taking steps.
Positive-duration history runs record 26 states: the initial state and 25 equally spaced output times, including the requested final time. Steps are shortened to reach those output times.
Each individual run has a limit of 500,000 steps. Exceeding this limit produces an error, not a result falsely labelled with the requested duration. Reducing grid resolution, annealing time or diffusivity can reduce the required number of steps.
Conservation is not an accuracy estimate
The inventory-error indicator measures numerical conservation. A very small value does not establish that the profile is sufficiently resolved. Compare results on progressively finer grids and check that features relevant to your interpretation change negligibly.
6. Understanding the result views
- Concentration: Final cation concentrations and matrix balance versus depth. The dashed vertical line marks the original source interface; it is not a moving marker.
- Diffusivity: Final local cation diffusivities and the matrix reference. The vertical coordinate is log₁₀[D / (m²/s)]. Zero values cannot appear on this logarithmic scale.
- 3D · Time: Concentration versus depth and elapsed annealing time. Select one cation at a time to inspect its surface.
- 3D · Temperature: Eight independent simulations spanning ±200 K around the selected temperature, restricted to 300–2500 K. Each run uses the same initial composition and annealing duration. This is a temperature comparison, not a heating ramp.
7. Automatic updates, cancellation and downloads
- Automatic updates: Parameter edits start a new calculation after a short pause. Previous results are cleared so they are not mistaken for results from the new inputs.
- Run simulation: Starts immediately, skipping the automatic-update delay. It also allows a cancelled calculation to be restarted without changing a parameter.
- Cancel: Stops the active background calculation. Cancelling a temperature sweep preserves an already completed base run.
- Download data: Exports the completed base run’s parameters, actual duration, cell-centre positions, concentrations and final diffusivities as tab-separated text. The download does not contain the full time history or temperature sweep.
Calculations run in a background worker to keep the interface responsive. The concentration and diffusivity plots work offline. Interactive 3D plots require an internet connection to load Plotly.
8. Physical scope and limitations
The cation equations interact through the optional matrix-dependent blocking factor, but they do not constitute a complete multicomponent transport model. The solver does not include cross-diffusion coefficients, electrostatic fields, charge compensation, reactions, phase changes or an evolving vacancy population.
The matrix is an algebraic balance species. If a parameter set produces a local cation sum above 100% beyond numerical tolerance, the simulation is rejected. Clipping the matrix to zero would conceal an inconsistent composition and would not fix the underlying transport model.
Why no Kirkendall marker is shown
The earlier marker calculation used an arbitrary scaling factor. It has been removed because it did not provide a physically justified displacement. Predicting Kirkendall motion requires a consistent intrinsic-flux and lattice-velocity formulation with an appropriate volume or site-density basis. The present concentration-balance model does not supply that formulation.
9. Corrections from version 5.4
- Replaced mass-losing endpoint updates with conservative boundary fluxes.
- Removed silent truncation of runs at the computational step limit.
- Corrected time overshoot and inaccurate history timestamps.
- Recalculated final diffusivities at every cell, including boundary cells.
- Added validation for geometry, numerical inputs and source composition.
- Removed hidden matrix clipping and the unsupported marker displacement.
- Clarified the empirical blocking factor, its floor and the preset limitations.
- Replaced Tailwind with embedded CSS and 2D chart dependencies with native SVG.
Diffusion_main.html · NitaD · Univ Paris-Saclay · September 2026
Illustrative model · No-flux boundaries · Cell-centered finite volumes