Hebraud-Lequeux (HL) Model: Equations and Physics

Reference: P. Hebraud and F. Lequeux, “Mode-Coupling Theory for the Pasty Rheology of Soft Glassy Materials,” Phys. Rev. Lett. 81, 2934 (1998). DOI: 10.1103/PhysRevLett.81.2934

1. Physical Picture

The HL model is a mean-field elastoplastic model for yield-stress fluids and soft glassy materials (foams, emulsions, dense suspensions). The material is divided into mesoscopic blocks, each carrying a local shear stress sigma. The central quantity is the probability density P(sigma, t) of finding a block at stress sigma at time t.

Each block evolves via four competing processes:

  1. Elastic loading: External shear rate gamma_dot loads stress at rate G_0 * gamma_dot (G_0 is the local elastic modulus, often normalized to 1).

  2. Mechanical noise (diffusion): Plastic events elsewhere create random stress kicks, modeled as Gaussian noise with diffusion coefficient D(t).

  3. Plastic yielding: When |sigma| > sigma_c, the block yields at rate 1/tau.

  4. Stress renewal: Yielded blocks are reinjected at sigma = 0 (delta function source).

2. The Fokker-Planck Equation (Core PDE)

The evolution of P(sigma, t) is governed by:

dP/dt = -G_0 * gamma_dot * dP/dsigma          [Elastic advection]
        + D(t) * d^2P/dsigma^2                 [Mechanical noise diffusion]
        - (1/tau) * P * Theta(|sigma| - sigma_c) [Plastic yielding sink]
        + Gamma(t) * delta(sigma)               [Renewal source at sigma=0]

where:

  • G_0: Local elastic modulus (often set to 1 in dimensionless form)

  • gamma_dot: Applied shear rate

  • sigma_c: Critical yield stress threshold

  • tau: Microscopic yielding timescale

  • Theta(x): Heaviside step function (1 if x > 0, else 0)

  • delta(sigma): Dirac delta function centered at sigma = 0

Normalization constraint:

integral from -inf to +inf of P(sigma, t) dsigma = 1   for all t

3. Self-Consistent Closure Relations

3a. Plastic activity (yielding rate):

Gamma(t) = (1/tau) * integral_{|sigma| > sigma_c} P(sigma, t) dsigma

This is the total rate of plastic rearrangements – the fraction of blocks with |sigma| > sigma_c, weighted by the attempt rate 1/tau.

3b. Diffusion coefficient (mechanical noise):

D(t) = alpha * Gamma(t)

where alpha is the dimensionless coupling parameter. This is the key self-consistent closure: the diffusion coefficient is proportional to the total plastic activity. More yielding -> more mechanical noise -> more yielding.

This closure makes the equation nonlinear and self-consistent: D depends on P through Gamma, and P depends on D.

3c. Conservation (mass balance):

The renewal source exactly balances the yielding sink:

Gamma(t) = (1/tau) * integral_{|sigma| > sigma_c} P(sigma, t) dsigma

ensuring that integral P dsigma = 1 is maintained.

4. Macroscopic Stress

The macroscopic shear stress is the first moment of the distribution:

sigma_macro(t) = integral from -inf to +inf of sigma * P(sigma, t) dsigma

The macroscopic viscosity is:

eta = sigma_macro / gamma_dot

5. Phase Behavior and Critical Point

The coupling parameter alpha controls the phase state:

alpha

Behavior

alpha < 1/2

Glassy/Jammed: Finite yield stress

alpha = 1/2

Critical point: Yield stress vanishes

alpha > 1/2

Fluid: No yield stress, Newtonian-like at low rates

Physical interpretation:

  • Small alpha: Weak mechanical noise coupling. Yielding events do not generate enough noise to trigger further events. The system jams and has a yield stress.

  • Large alpha: Strong coupling. Each yielding event efficiently triggers further events (positive feedback). The system flows easily.

6. Steady-State Flow Curve

In steady state (dP/dt = 0), the flow curve sigma(gamma_dot) is obtained by solving:

0 = -G_0 * gamma_dot * dP/dsigma + D_ss * d^2P/dsigma^2
    - (1/tau) * P * Theta(|sigma| - sigma_c) + Gamma_ss * delta(sigma)

with D_ss = alpha * Gamma_ss determined self-consistently.

Known analytical results:

For alpha < 1/2 (Glassy phase, yield stress):

Near the yield point (gamma_dot -> 0+):

sigma(gamma_dot) = sigma_y + A * gamma_dot^(1/2) + ...

This is a Herschel-Bulkley flow curve with exponent n = 1/2 (equivalently, beta = 2 in sigma - sigma_y ~ gamma_dot^(1/beta)).

The yield stress sigma_y depends on alpha and sigma_c. It is determined implicitly by the steady-state equation at gamma_dot = 0+.

At large shear rates (gamma_dot -> infinity):

sigma(gamma_dot) ~ eta_inf * gamma_dot

The material behaves as a Newtonian fluid with effective viscosity eta_inf. This was rigorously proven in Olivier (2012), Sci. China Math. 55, 435-452.

For alpha = 1/2 (Critical point):

sigma(gamma_dot) ~ gamma_dot^(1/2)    (power-law fluid, no yield stress)

For alpha > 1/2 (Fluid phase):

sigma(gamma_dot) ~ eta_0 * gamma_dot   at low shear rates (Newtonian)
sigma(gamma_dot) ~ eta_inf * gamma_dot  at high shear rates (Newtonian)

with possible shear-thinning in between.

7. Dimensionless Form

Setting G_0 = 1 and sigma_c = 1, the model has only two dimensionless parameters:

  • alpha: coupling strength

  • gamma_dot * tau: Weissenberg number (Wi)

The dimensionless PDE becomes:

dP/dt = -gamma_dot * dP/dsigma + D * d^2P/dsigma^2
        - P * Theta(|sigma| - 1) + Gamma * delta(sigma)

with D = alpha * Gamma and Gamma = integral_{|sigma|>1} P dsigma (using tau = 1).

8. Aging Dynamics (No Shear, gamma_dot = 0)

In the absence of shear (Sollich, Olivier, Bresch 2016, arXiv:1611.06681):

  • The stress diffusion constant decays as D(t) ~ 1/t^2 during aging.

  • The cumulative memory integral (integral of D(t) dt) is finite.

  • The shear stress relaxation function decays only to a plateau (not to zero).

  • The system becomes progressively more elastic as it ages.

  • Frequency-dependent shear modulus: G’’(omega) decreases with system age, while relaxation times scale linearly with age.

9. Extensions and Variants

9a. Bocquet-Colin-Ajdari (BCA) model

Extension to include spatial heterogeneity and flow coupling. Adds a Stokes-like equation coupling the stress distribution to the local flow field.

9b. Thermal HL (Baron & Biroli, 2023, arXiv:2312.03627)

Generalization to include thermal fluctuations. The yielding rate becomes:

nu(sigma, sigma_c) = exp(-(sigma_c - |sigma|)/T)  for |sigma| < sigma_c
nu(sigma, sigma_c) = 1                             for |sigma| >= sigma_c

This introduces a temperature T and activated dynamics below sigma_c. Predicts MCT-like power-law divergence near T_c and Arrhenius below T_c.

9c. Multidimensional generalization (Olivier, 2012, HAL hal-01263132)

Extension to tensorial stress in multidimensional flow configurations.

9d. Generalized HL with different yielding rules

  • Different Herschel-Bulkley exponents beta can arise from non-Gaussian mechanical noise or modified yielding criteria (Lin & Wyart, 2018, PRX 8, 011005).

  • The standard HL Gaussian noise gives beta = 2 (i.e., n = 1/2).

  • Non-Gaussian (e.g., Levy) noise gives different exponents.

10. Benchmark Results for Numerical Verification

10a. Flow curve shape

  • For alpha = 0.3, tau = 1, sigma_c = 1:

    • Yield stress sigma_y > 0 (finite, since alpha < 0.5)

    • Near yield: sigma ~ sigma_y + A * sqrt(gamma_dot)

    • High gamma_dot: sigma ~ eta_inf * gamma_dot (Newtonian)

10b. Phase transition

  • At alpha = 0.5: yield stress vanishes continuously

  • For alpha just below 0.5: very small yield stress

10c. Probability distribution shape

  • Steady state at finite gamma_dot: P(sigma) is asymmetric (skewed in flow direction)

  • P(sigma) = 0 is NOT reached; distribution has finite support up to edges

  • Center bin (sigma = 0) has a peak from the renewal source

10d. Macroscopic stress

  • sigma_macro = integral(sigma * P * dsigma) is always positive for gamma_dot > 0

  • For gamma_dot = 0 from rest: sigma_macro = 0

10e. Conservation check

  • integral(P * dsigma) must remain exactly 1 at all times

  • Activity Gamma must match the yielded fraction / tau

10f. Large shear rate asymptote

  • sigma / gamma_dot -> constant (Newtonian viscosity) as gamma_dot -> infinity

  • This is proven rigorously in the mathematical literature

11. Numerical Implementation Notes

The PDE is solved on a finite grid sigma in [-sigma_max, sigma_max] with n_bins points. The existing RheoJAX implementation uses:

  1. Finite Volume / Finite Difference discretization

  2. Upwind scheme for the advection term (elastic loading)

  3. Central differences for the diffusion term

  4. Explicit Euler time stepping with CFL-based sub-stepping

  5. Delta function approximated as mass injection into the center bin

  6. Boundary conditions: P = 0 at grid edges

  7. Positivity enforcement and renormalization at each step

Key numerical parameters:

  • Grid: sigma in [-5sigma_c, 5sigma_c], 501 bins (default)

  • Sub-stepping: 25 fixed sub-steps per outer time step for CFL stability

  • Max scan steps: 20000 (flow curve), 500 (creep, due to XLA compilation cost)

References

  1. P. Hebraud, F. Lequeux, “Mode-Coupling Theory for the Pasty Rheology of Soft Glassy Materials,” Phys. Rev. Lett. 81, 2934 (1998). https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.81.2934

  2. P. Sollich, “Rheological constitutive equation for a model of soft glassy materials,” Phys. Rev. E 58, 738 (1998). arXiv:cond-mat/9712001. https://arxiv.org/abs/cond-mat/9712001

  3. P. Sollich, J. Olivier, D. Bresch, “Aging and linear response in the Hebraud-Lequeux model for amorphous rheology,” arXiv:1611.06681 (2016). https://arxiv.org/abs/1611.06681

  4. A. Gimenez et al., “Stability and Numerical Analysis of the Hebraud-Lequeux Model for Suspensions,” Discrete Dyn. Nat. Soc. 2011, 415921 (2011). https://onlinelibrary.wiley.com/doi/10.1155/2011/415921

  5. J. Olivier, “Large shear rate behavior for the Hebraud-Lequeux model,” Sci. China Math. 55, 435-452 (2012). https://link.springer.com/article/10.1007/s11425-011-4350-2

  6. Asymptotic analysis in flow curves for a model of soft glassy fluids, Z. Angew. Math. Phys. 61, 445-466 (2010). https://link.springer.com/article/10.1007/s00033-009-0022-2

  7. J. Lin, M. Wyart, “Microscopic processes controlling the Herschel-Bulkley exponent,” Phys. Rev. X 8, 011005 (2018). arXiv:1708.00516.

  8. J.W. Baron, G. Biroli, “Mean-Field Analysis of the Glassy Dynamics of an Elastoplastic Model of Super-Cooled Liquids,” arXiv:2312.03627 (2023). https://arxiv.org/html/2312.03627

  9. L. Bocquet, A. Colin, A. Ajdari, “Kinetic Theory of Plastic Flow in Soft Glassy Materials,” Phys. Rev. Lett. 103, 036001 (2009).