DPM Brownian Motion -- Langevin Equation, Diffusion Coefficient, UDF for PM2.5

#+CATEGORY: cfd

#+TAGS: DPM Brownian Langevin diffusion PM2.5 nanoparticle stochastic UDF body-force

#+MATURITY: growing

DPM Brownian Motion :: Langevin Equation and UDF Implementation for PM2.5

Brownian motion dominates the dynamics of particles below approximately 1 micrometre. In a shisha simulation, PM2.5 particles (0.1--2.5 micrometres) experience random collisions with air molecules that produce displacement comparable to or larger than convective transport. Fluent's standard DPM does not include Brownian force by default; it must be added via DEFINE_DPM_BODY_FORCE. This page derives the model from the Langevin equation, gives the diffusion coefficient and force magnitude estimates for PM2.5, and provides the complete UDF.

Why Brownian Motion Matters for PM2.5

The Brownian displacement over time t scales as:

where D is the diffusion coefficient. For a spherical particle:

For PM2.5, Brownian displacement is small in absolute terms but becomes significant relative to the mesh cell size. If the cell size is 0.5 mm, then a 2.5 um particle's Brownian displacement (0.011 mm/s) is a non-trivial fraction; and for the sub-micron fraction (the majority of shisha aerosol by count), Brownian motion is essential.

The Langevin Equation

The full equation of motion for a particle under Brownian forcing:

The Brownian force is modelled as a Gaussian white noise:

where:

• \boldsymbol{\xi}(t) is a standard normal random vector (zero mean, unit variance)
• =k_B = 1.381 \times 10^{-23}= J/K (Boltzmann constant)
• T is the local fluid temperature [K]
• f is the drag coefficient (force per unit velocity)
• \Delta t is the particle time step

For Stokes drag, =f = 3 \pi \mu d_p=, so:

This is the force per unit particle mass is returned from DEFINE_DPM_BODY_FORCE:

Since =m_p = \rho_p \pi d_p^3 / 6=:

Key insight

The Brownian acceleration scales as d_p^{-3/2}. Halving the particle diameter increases the Brownian acceleration by 2^{3/2} \approx 2.83. This is why Brownian motion is negligible for large particles but dominates for sub-micron aerosols.

When to Enable Brownian Motion

The ratio of Brownian displacement to convective displacement over one cell transit is a Péclet-like number:

For the shisha simulation:

• Typical velocity in chamber: 0.1--0.5 m/s
• Cell size: 0.001--0.005 m
• D for 1 um particle: 2.3e-10 m^2/s
• =Pe_p = 0.3 * 0.003 / 2.3e-10 \approx 4 \times 10^6=

So in the main chamber, Brownian motion is negligible for 1 um particles. However, in the sub-micron fraction (0.1--0.5 um) and in low-velocity regions (near walls, in recirculation zones, at the outlet far field), Pe_p drops substantially and Brownian motion contributes to wall deposition and dispersion.

Recommendation for shisha PM2.5

Enable Brownian force if: 1. The particle size distribution has significant sub-micron count (< 0.5 um) 2. You care about near-wall deposition patterns 3. You want physically realistic dispersion rather than deterministic filaments

Disable Brownian force if: 1. All particles are > 2 um (mass distribution dominated by large particles) 2. You are only interested in bulk concentration (not deposition) 3. DRW stochastic tracking alone provides sufficient dispersion

Complete Brownian Force UDF

See cfd_dpm-udf for macro reference. The implementation below uses the Box-Muller transform to generate standard normal random numbers.

  #include "udf.h"
  #include "dpm.h"
  #include <math.h>

  /* Physical constants */
  #define KB_BROWNIAN  1.38065e-23   /* Boltzmann constant [J/K] */
  #define PI_BROWNIAN  3.14159265359

  DEFINE_DPM_BODY_FORCE(brownian_force_pm25, p, i)
  {
    /* 
     * Returns Brownian force per unit mass in direction i.
     * Model: F_B/m = xi * sqrt(2 kB T f / (dt * m_p))
     * For Stokes drag: f = 3 pi mu d_p
     */
    
    /* Get local cell properties */
    Thread *t = P_CELL_THREAD(p);
    cell_t c  = P_CELL(p);
    
    /* Fluid temperature [K] */
    real T = C_T(c, t);
    
    /* Fluid dynamic viscosity [Pa s] */
    /* Use laminar viscosity for Brownian force; 
       turbulent viscosity acts through DRW separately */
    real mu = C_MU_L(c, t);
    
    /* Particle diameter [m] */
    real d = P_DIAM(p);
    
    /* Particle mass [kg] */
    real mp = P_MASS(p);
    
    /* Particle time step [s] */
    real dt = CURRENT_TIMESTEP;
    if (dt <= 0.0) dt = 1e-6;
    
    /* Guard: if particle has effectively zero mass, 
       return zero to avoid singularity */
    if (mp < 1e-30) return 0.0;
    
    /* Drag coefficient for Stokes regime */
    real f_drag = 3.0 * PI_BROWNIAN * mu * d;
    
    /* Generate standard normal random number (Box-Muller transform) */
    real u1 = ((real)rand() + 1.0) / ((real)RAND_MAX + 2.0);
    real u2 = ((real)rand() + 1.0) / ((real)RAND_MAX + 2.0);
    
    /* Clamp u1 away from zero to avoid log(0) */
    if (u1 < 1e-20) u1 = 1e-20;
    
    real xi = sqrt(-2.0 * log(u1)) * cos(2.0 * PI_BROWNIAN * u2);
    
    /* Brownian force magnitude per unit mass [m/s^2] */
    real sigma_sq = 2.0 * KB_BROWNIAN * T * f_drag / (dt * mp);
    
    /* Clamp to avoid sqrt of negative due to floating-point */
    if (sigma_sq < 0.0) sigma_sq = 0.0;
    
    real force = xi * sqrt(sigma_sq);
    
    return force;
  }

Activation in Fluent

1. Compile the UDF (Models > Discrete Phase > Interaction > drag user-defined not needed; just compile) 2. = define/models/dpm/non-interactive/enable = (if not already on) 3. Models > Discrete Phase > Physical Models:

• Enable Body Forces
• Under Body Forces, set x/y/z Body Force to User-Defined and select brownian_force_pm25

4. Ensure Stochastic Tracking is enabled (Number of Tries \geq 10)

Mesh and Time-Step Requirements

Brownian motion imposes two additional mesh constraints beyond the standard DPM requirements:

Cell size vs. Brownian displacement

The cell size must resolve the Brownian displacement over one time step:

For =D = 2.3 \times 10^{-10}= m^2/s (1 um particle) and =\Delta t_p = 10^{-4}= s: =\sqrt{2 D \Delta t_p} = 0.21= nm. This is far below any practical mesh, so the mesh constraint from Brownian motion alone is not active.

Time step constraint

The Brownian force is a white-noise approximation; it is only valid when the correlation time of the molecular collisions (approximately 10^{-10} s for air at STP) is much shorter than the particle time step. This is always satisfied for CFD time steps.

The more active constraint is numerical: the Brownian acceleration must not cause the particle to traverse more than one cell per time step:

For PM2.5: a_{Brownian} \sim 0.1 m/s^2, =\Delta x = 0.001= m, so \Delta t_p < 0.01 s. Standard Fluent DPM sub-timesteps are much smaller than this.

Coupling with DRW (Stochastic Tracking)

There are two separate stochastic mechanisms acting on PM2.5 particles:

These are independent and additive. DRW represents dispersion by resolved + unresolved turbulence eddies; Brownian represents dispersion by molecular bombardment. Both should be enabled for PM2.5.

Important: do not double-count

Do NOT enable Brownian Motion as a separate Fluent model (there is a toggle in newer versions) and the UDF simultaneously. Pick one implementation. The UDF is preferred because: 1. It matches the Langevin equation exactly 2. You can validate the force magnitude analytically 3. It is portable across Fluent versions

Number of tries

For the combined DRW + Brownian model, statistical convergence requires more tries than DRW alone:

Validation of Brownian Implementation

Test case: quiescent air

1. Set up a 2D box (0.01 x 0.01 m) with zero velocity 2. Inject 10000 parcels of 1 um particles at the center 3. Enable only Brownian force (disable DRW, gravity, drag coupling) 4. Run for 10 seconds 5. Compute mean-square displacement from particle positions 6. Compare to =\langle r^2 \rangle = 2 \cdot 2 \cdot D \cdot t= (2D in a 2D plane)

Expected D for 1 um particle: 2.3 \times 10^{-10} m^2/s Expected \langle r^2 \rangle after 10 s: 9.2 \times 10^{-9} m^2 Expected RMS displacement: 0.096 mm

If the simulation matches this within 10--20% (statistical variation), the Brownian UDF is correct.

Test case: compare to analytical diffusion profile

In quiescent air, particles from a point source spread as a Gaussian:

Bin the particle positions at time t into radial shells, compute the histogram, and compare to this analytical Gaussian. The fit should be excellent for > 5000 parcels.

PM2.5-Specific Considerations for Shisha

For the shisha simulation, the particle size distribution is typically log-normal with:

• Geometric mean diameter: d_g \sim 0.2--0.5 um (sub-micron peak)
• Geometric standard deviation: \sigma_g \sim 1.5--2.0

The number distribution is dominated by sub-micron particles (where Brownian is strongest), while the mass distribution is dominated by particles near the PM2.5 boundary (where Brownian is weaker but still present).

Implications for the simulation

1. Number concentration predictions require Brownian motion for accuracy 2. Mass concentration predictions may be acceptable with DRW alone, but adding Brownian improves dispersion near walls 3. Deposition on walls is strongly affected by Brownian diffusion for d_p < 0.5 um -- without Brownian, you will underpredict deposition in the bowl chamber and hose 4. Temperature gradients: thermophoretic force (see cfd_dpm-physical-models) may dominate over Brownian near hot surfaces. Both should be enabled simultaneously.

Force ranking for PM2.5 in shisha

See Also