Deterministic Elastic Rough-Surface Contact Calculator (FFT-BEM)

Surface, material, and load

Use a square power-of-two grid. Start at 32 × 32, then repeat at finer resolution to assess grid sensitivity.

Elastic contact result

Pressure is compressive and nonnegative. Area and peak pressure are discretization-sensitive.

Real contact fraction—

Preparing surface preview…

Contact pressure, p(x,y)

Contact pressure (GPa)

———

Zero pressure to the calculated peak pressure after simulation.

Equations used

SI units are used internally. Heights are positive toward the opposing flat; pressure and load are positive in compression. The zero-height datum is the map mean.

Generated self-affine fractal surface

The generated option uses seeded Fourier spectral synthesis rather than spatial Gaussian smoothing. Random complex Fourier coefficients have a radial power spectral density (PSD) with a flat long-wavelength roll-off, a self-affine band, and a short-wavelength cutoff. An inverse 2-D FFT produces z(x,y), which is mean-centered and scaled to the selected Rq.

q_r = 2π/λ_r; q_s = 2π/λ_s; C(q) = C₀ for q ≤ q_r; C(q) = C₀(q/q_r)^−2(1+H) for q_r < q ≤ q_s; C(q) = 0 for q > q_s

z = R_q[IFFT₂(√C(q) ξ(q)) − mean]/RMS, with ξ(−q) = ξ(q)^*

Here H is the dimensionless Hurst exponent; λ_s is the short-wavelength cutoff in m and must satisfy λ_s ≥ 2Δ; λ_r is the long-wavelength roll-off in m and satisfies λ_s ≤ λ_r ≤ NΔ; q, q_r, and q_s are radial spatial angular frequencies in rad/m; C is the isotropic height PSD; C₀ is an arbitrary positive spectral normalization whose magnitude cancels when the synthesized map is scaled to R_q, so it is not a calibrated material property or user input; ξ is seeded complex Gaussian noise with Hermitian symmetry; the asterisk denotes complex conjugation; IFFT₂ is the inverse two-dimensional discrete Fourier transform; and RMS is the unscaled generated root-mean-square height.

Elastic contact

1/E^* = (1 − ν₁²)/E₁ + (1 − ν₂²)/E₂

u_ij = Σ_kl K_ij,klp_kl, where K = (1/πE^*) ∬_cell(1/r)dA

g_ij = d − z_ij + u_ij; p_ij ≥ 0; g_ij ≥ 0; p_ijg_ij = 0

p_c = max(1 Pa, p_max × 10⁻⁸); N_c = count(p_ij > p_c)

W_calc = Σp_ijΔ²; A₀ = N²Δ²; A_r = N_cΔ²; A_r/A₀ = N_c/N²; p̄ = W_calc/A_r

δ = max(z) − d; k_n = W_calc/δ

ε_W = |W_calc − W|/W; ε_g = max(max_pij>0|g_ij|, max_pij=0max(0, −g_ij))

tolerance_m = ε × max(R_q, 0.05Δ) × 0.05; converged when ε_W ≤ ε and ε_g ≤ 20 × tolerance_m

Symbols and units: E^*, E₁, and E₂ are the contact reduced and body Young’s moduli in Pa; ν₁ and ν₂ are dimensionless Poisson ratios; i and j index observation cells; k and l index source cells; z_ij is mean-centered surface height in m; u_ij is elastic normal displacement in m; d is mean-plane separation in m; g_ij is local gap in m; p_ij and p_kl are cell-average compressive pressures in Pa; p_max is peak cell pressure in Pa; p_c is the numerical contact-count threshold in Pa; r is the in-plane source-to-observation distance in m; dA is differential cell area in m²; K_ij,kl is rectangular-cell half-space compliance in m/Pa; K₀₀ is the self-cell compliance in m/Pa; Δ is square-cell spacing in m; N is the number of cells per side; N_c is the number of cells above p_c; A₀ and A_r are nominal and discretized real areas in m²; W is the target compressive load in N; W_calc is the integrated calculated load in N; p̄ is mean contact pressure in Pa; δ is peak indentation in m; k_n is secant normal stiffness in N/m; R_q is RMS height in m; ε is the selected dimensionless tolerance; ε_W is relative load residual; ε_g and tolerance_m are dimensional gap residual and tolerance in m.

Numerical method

The rectangular-cell influence kernel is zero-padded and convolved with pressure by a radix-2 two-dimensional FFT. Projected iterations enforce nonnegative pressure and unilateral contact; an inner step requires both ε_g ≤ tolerance_m and the maximum pressure update multiplied by the self-cell compliance K₀₀ to be ≤ tolerance_m. Separation is bisected until integrated pressure matches the selected load. The final convergence flag uses the load and gap criteria shown above. Residuals test numerical balance only.

Assumptions and limits

The model assumes frictionless normal contact between homogeneous isotropic linear-elastic half-spaces, small strain, a single-valued topography, square uniform sampling, and no adhesion. It does not include yielding, hardness caps, plastic residual shape, coatings, tangential traction, friction, viscoelasticity, subsurface stress, debris, or thermal effects. Keep contact away from map edges and repeat with a finer grid and larger sampled area. Do not use after yielding or for safety-critical qualification without specialist review.

References

Stanley and Kato (1997), doi:10.1115/1.2833523.
Polonsky and Keer (1999), doi:10.1016/S0043-1648(99)00113-1.
Love, half-space elasticity treatment.

How to use

Select Generated surface for a seeded self-affine fractal topography or Measured grid for pasted/uploaded heights.
Choose a square power-of-two grid and enter pixel spacing Δ in µm. In measured mode, provide exactly N rows of N finite heights in µm; comma, semicolon, or whitespace separators are accepted.
For a generated surface, set RMS roughness Rq, Hurst exponent H, short-wavelength cutoff λs, long-wavelength roll-off λr, and an integer seed. Use 2Δ ≤ λs ≤ λr ≤ NΔ.
Select Generate / preview surface. Inspect the displayed topography and its height colorbar in µm before proceeding. Changing any surface setting disables simulation until a fresh preview is generated.
Enter E₁, ν₁, E₂, and ν₂ for the two bodies. Moduli must be positive and −1 < ν < 0.5.
Enter total compressive normal load W in N and select a numerical tolerance, then select Run simulation. Use Reset example to restore the 32 × 32 fractal steel example without automatically running contact.
Read the canonical real contact fraction, integrated load, area, peak/mean pressure, indentation, secant stiffness, convergence residuals, and pressure map. Use the pressure colorbar in GPa to interpret the heat map.
Correct any red input or non-convergence message before interpreting results. A small residual does not establish physical validity.
Repeat at finer spacing/grid size and a larger physical patch. Reject results that change materially or whose pressure reaches the map boundary.
Do not use this elastic model when pressure can cause yielding, adhesion/coatings matter, the contact is conformal, tangential loading dominates, or the result is a safety-critical design decision.

Interpretation

Real area counts cells whose pressure exceeds max(1 Pa, 10⁻⁸ of peak pressure); it changes with mesh resolution.
Peak pressure is a cell average and is especially mesh-sensitive.
The reported stiffness is W_calc/δ, a secant screening value rather than the tangent slope dW/dδ.
Generated surfaces are reproducible for a fixed seed and settings, but one realization is not a statistical uncertainty study.