Direct simulation Monte Carlo

<h2 id="dsmc-algorithm">DSMC Algorithm</h2>
<p>The direct simulation Monte Carlo algorithm is like <a href="/facts/Molecular_dynamics/yR3DNt6U">molecular dynamics</a> in that the state of
the system is given by the positions and velocities of the
particles, 
  
    
      
        {
        
          
            r
          
          
            i
          
        
        ,
        
          
            
              v
            
          
          
            i
          
        
        }
      
    
    {\displaystyle \{\mathbf {r} _{i},{\textbf {v}}_{i}\}}
  
, for 
  
    
      
        i
        =
        1
        ,
        …
        ,
        N
      
    
    {\displaystyle i=1,\ldots ,N}
  
. 
Unlike molecular dynamics, each particle in a DSMC simulation represents 
  
    
      
        
          F
          
            N
          
        
      
    
    {\displaystyle F_{N}}
  
 molecules in
the physical system that have roughly the same position and velocity. 
This allows DSMC to rescale length and time for the modeling of macroscopic systems (e.g., <a href="/facts/Atmospheric_entry/p7TMk9sa">atmospheric entry</a>). 
Specifically, the system volume is 
  
    
      
        V
        =
        (
        N
        
          F
          
            N
          
        
        )
        
          /
        
        n
      
    
    {\displaystyle V=(NF_{N})/n}
  
, where 
  
    
      
        n
      
    
    {\displaystyle n}
  
 is the number
density and each collision between simulation particles represents 
  
    
      
        
          F
          
            N
          
        
      
    
    {\displaystyle F_{N}}
  
 collisions
among molecules in the physical system.
As a rule of thumb there should be 20 or more particles per cubic mean free path
for accurate results.
</p><p>The evolution of the system is integrated in time steps, 
  
    
      
        τ
      
    
    {\displaystyle \tau }
  
, which are
typically on the order of the mean collision time for a particle. 
At each time step all the particles are moved and then a random set of pairs collide.
In the absence of external fields (e.g., gravity) the particles move ballistically as

r
          
          
            i
          
        
        (
        t
        +
        τ
        )
        =
        
          
            r
          
          
            i
          
        
        (
        t
        )
        +
        
          
            v
          
          
            i
          
        
        (
        t
        )
        τ
      
    
    {\displaystyle \mathbf {r} _{i}(t+\tau )=\mathbf {r} _{i}(t)+\mathbf {v} _{i}(t)\tau }
  
.
Any particle that reaches a boundary or a surface has its position and velocity reset accordingly
(e.g., <a href="/facts/Periodic_boundary_conditions/8hdzjIAD">periodic boundary conditions</a>).
After all the particles have moved, they are sorted into cells and some are randomly selected to collide.
based on probabilities and collision rates obtained from the <a href="/facts/Kinetic_theory_of_gases/WUaJuMHU">kinetic theory of gases</a>.
After the velocities of all colliding particles have been reset, statistical sampling is performed and then
the process is repeated for the next time step.
</p>
<h3>Collisions</h3>
<p>On each timestep the particles are sorted into spatial cells and only particles in the same cell
are allowed to collide. Typically the dimension of a cell is no larger than a mean free path.
All pairs of particles in a cell are candidate collision partners, regardless of their actual trajectories. 
</p><p>The details of how collisions are calculated in DSMC depend on the molecular interaction model;
here we take the <a href="/facts/Hard_spheres/91Gfe05N">hard spheres</a> model, which is the simplest.
In the hard spheres model, the collision probability for the pair of particles, 
  
    
      
        i
      
    
    {\displaystyle i}
  
 and 
  
    
      
        j
      
    
    {\displaystyle j}
  
, is
proportional to their relative speed,

P
          
            
              c
              o
              l
              l
            
          
        
        [
        i
        ,
        j
        ]
        =
        
          
            
              
                |
              
              
                
                  v
                
                
                  i
                
              
              −
              
                
                  v
                
                
                  j
                
              
              
                |
              
            
            
              
                ∑
                
                  m
                  =
                  1
                
                
                  
                    N
                    
                      
                        c
                      
                    
                  
                
              
              
                ∑
                
                  n
                  =
                  1
                
                
                  m
                  −
                  1
                
              
              
                |
              
              
                
                  v
                
                
                  m
                
              
              −
              
                
                  v
                
                
                  n
                
              
              
                |
              
            
          
        
      
    
    {\displaystyle P_{\mathrm {coll} }[i,j]={{|\mathbf {v} _{i}-\mathbf {v} _{j}|} \over {\sum _{m=1}^{N_{\mathrm {c} }}\sum _{n=1}^{m-1}|\mathbf {v} _{m}-\mathbf {v} _{n}|}}}

where 
  
    
      
        
          N
          
            
              c
            
          
        
      
    
    {\displaystyle N_{\mathrm {c} }}
  
 is the number of particles in the cell and the summations are over particles within the cell.
Because of the double sum in the denominator it can be computationally expensive to use this collision probability directly.
Instead, the following <a href="/facts/Rejection_sampling/6kj9Bahq">rejection sampling</a> scheme can be used to select collision pairs:
</p>
<ol><li>A pair of candidate particles, 
  
    
      
        i
      
    
    {\displaystyle i}
  
 and 
  
    
      
        j
      
    
    {\displaystyle j}
  
, is chosen at random and their relative speed, 
  
    
      
        
          v
          
            
              r
            
          
        
        =
        
          |
        
        
          
            v
          
          
            i
          
        
        −
        
          
            v
          
          
            j
          
        
        
          |
        
      
    
    {\displaystyle v_{\mathrm {r} }=|\mathbf {v} _{i}-\mathbf {v} _{j}|}
  
, is computed.</li>
<li>The pair is accepted as collision partners if 
  
    
      
        
          v
          
            
              r
            
          
        
        >
        
          v
          
            
              r
            
          
          
            
              m
              a
              x
            
          
        
        ℜ
      
    
    {\displaystyle v_{\mathrm {r} }>v_{\mathrm {r} }^{\mathrm {max} }\Re }
  
, where 
  
    
      
        
          v
          
            
              r
            
          
          
            
              m
              a
              x
            
          
        
      
    
    {\displaystyle v_{\mathrm {r} }^{\mathrm {max} }}
  
 is the maximum relative speed in the cell and 
  
    
      
        ℜ
      
    
    {\displaystyle \Re }
  
 is a <a href="/facts/Continuous_uniform_distribution/XbnlVljT"> uniform deviate</a> in [0, 1).</li>
<li>If the pair is accepted, the collision is processed; the velocities of the particles are reset but positions are unchanged.</li>
<li>After the collision is processed or if the pair is rejected, return to step 1.</li></ol>
<p>This procedure is correct even if the value
of 
  
    
      
        
          v
          
            
              r
            
          
          
            
              m
              a
              x
            
          
        
      
    
    {\displaystyle v_{\mathrm {r} }^{\mathrm {max} }}
  
 is  overestimated, although it is less efficient
in the sense that more candidates are rejected.
</p><p>After the collision pair is chosen, their post-collision velocities,

v
          
          
            i
          
          
            ∗
          
        
      
    
    {\displaystyle \mathbf {v} _{i}^{*}}
  
 and 
  
    
      
        
          
            v
          
          
            j
          
          
            ∗
          
        
      
    
    {\displaystyle \mathbf {v} _{j}^{*}}
  
, are evaluated.
Writing the relative velocity in terms of <a href="/facts/Spherical_coordinate_system/5cVnNof3"> spherical angles</a>, 
  
    
      
        θ
      
    
    {\displaystyle \theta }
  
 and 
  
    
      
        ϕ
      
    
    {\displaystyle \phi }

v
          
          
            
              r
            
          
          
            ∗
          
        
        =
        
          v
          
            
              r
            
          
        
        [
        (
        sin
        ⁡
        θ
        cos
        ⁡
        ϕ
        )
        
          
            
              
                x
              
              ^
            
          
        
        +
        (
        sin
        ⁡
        θ
        sin
        ⁡
        ϕ
        )
        
          
            
              
                y
              
              ^
            
          
        
        +
        cos
        ⁡
        θ
        
        
          
            
              
                z
              
              ^
            
          
        
        ]
      
    
    {\displaystyle \mathbf {v} _{\mathrm {r} }^{*}=v_{\mathrm {r} }[(\sin \theta \cos \phi ){\hat {\mathbf {x} }}+(\sin \theta \sin \phi ){\hat {\mathbf {y} }}+\cos \theta \,{\hat {\mathbf {z} }}]}

these angles are selected by a Monte Carlo process with distributions given by the collision model.
For the hard spheres model these angles are uniformly distributed over the <a href="/facts/Unit_sphere/5UeoTcug">unit sphere</a>. 
The azimuthal angle is uniformly distributed between 0 and 
  
    
      
        2
        π
      
    
    {\displaystyle 2\pi }
  
, so it is selected as 
  
    
      
        ϕ
        =
        2
        π
        
          ℜ
          
            1
          
        
      
    
    {\displaystyle \phi =2\pi \Re _{1}}

where 
  
    
      
        
          ℜ
          
            1
          
        
      
    
    {\displaystyle \Re _{1}}
  
 is a <a href="/facts/Continuous_uniform_distribution/XbnlVljT"> uniform deviate</a> in [0, 1).
The polar angle is distributed according to the probability density,

P
          
            θ
          
        
        (
        θ
        )
        
        d
        θ
        =
        
          
            
              
                1
                2
              
            
          
        
        sin
        ⁡
        θ
        
        d
        θ
      
    
    {\displaystyle P_{\theta }(\theta )\,d\theta ={\textstyle {\frac {1}{2}}}\sin \theta \,d\theta }

Using the change of variable 
  
    
      
        q
        =
        −
        cos
        ⁡
        θ
      
    
    {\displaystyle q=-\cos \theta }
  
, we have 
  
    
      
        
          P
          
            q
          
        
        (
        q
        )
        
        d
        q
        =
        (
        
          
            
              
                1
                2
              
            
          
        
        )
        
        d
        q
      
    
    {\displaystyle P_{q}(q)\,dq=({\textstyle {\frac {1}{2}}})\,dq}
  
 so

−
        cos
        ⁡
        θ
        =
        q
         
        
          a
          n
          d
        
         
        sin
        ⁡
        θ
        =
        
          
            1
            −
            
              q
              
                2
              
            
          
        
         
        
          w
          h
          e
          r
          e
        
         
        q
        =
        2
        
          ℜ
          
            2
          
        
        −
        1
      
    
    {\displaystyle -\cos \theta =q~\mathrm {and} ~\sin \theta ={\sqrt {1-q^{2}}}~\mathrm {where} ~q=2\Re _{2}-1}

The post-collision velocities are set as

v
          
          
            i
          
          
            ∗
          
        
        =
        
          
            v
          
          
            
              c
              m
            
          
          
            ∗
          
        
        +
        
          
            1
            2
          
        
        
          
            v
          
          
            
              r
            
          
          
            ∗
          
        
        
        
          
            v
          
          
            j
          
          
            ∗
          
        
        =
        
          
            v
          
          
            
              c
              m
            
          
          
            ∗
          
        
        −
        
          
            1
            2
          
        
        
          
            v
          
          
            
              r
            
          
          
            ∗
          
        
      
    
    {\displaystyle \mathbf {v} _{i}^{*}=\mathbf {v} _{\mathrm {cm} }^{*}+{1 \over 2}\mathbf {v} _{\mathrm {r} }^{*}\qquad \mathbf {v} _{j}^{*}=\mathbf {v} _{\mathrm {cm} }^{*}-{1 \over 2}\mathbf {v} _{\mathrm {r} }^{*}}

Note that by conservation of linear momentum and energy the center of mass velocity
and the relative speed are unchanged in a collision. That is,

v
          
          
            
              c
              m
            
          
        
        =
        
          
            1
            2
          
        
        (
        
          
            v
          
          
            i
          
        
        +
        
          
            v
          
          
            j
          
        
        )
        =
        
          
            1
            2
          
        
        (
        
          
            v
          
          
            i
          
          
            ∗
          
        
        +
        
          
            v
          
          
            j
          
          
            ∗
          
        
        )
        =
        
          
            v
          
          
            
              c
              m
            
          
          
            ∗
          
        
      
    
    {\displaystyle \mathbf {v} _{\mathrm {cm} }={1 \over 2}(\mathbf {v} _{i}+\mathbf {v} _{j})={1 \over 2}(\mathbf {v} _{i}^{*}+\mathbf {v} _{j}^{*})=\mathbf {v} _{\mathrm {cm} }^{*}}

and

v
          
            
              r
            
          
        
        =
        
          |
        
        
          
            v
          
          
            i
          
        
        −
        
          
            v
          
          
            j
          
        
        
          |
        
        =
        
          |
        
        
          
            v
          
          
            i
          
          
            ∗
          
        
        −
        
          
            v
          
          
            j
          
          
            ∗
          
        
        
          |
        
        =
        
          v
          
            
              r
            
          
          
            ∗
          
        
      
    
    {\displaystyle v_{\mathrm {r} }=|\mathbf {v} _{i}-\mathbf {v} _{j}|=|\mathbf {v} _{i}^{*}-\mathbf {v} _{j}^{*}|=v_{\mathrm {r} }^{*}}

This process is repeated for every pair of colliding particles.
</p><p>From the collision frequency, 
  
    
      
        
          f
          
            
              c
              o
              l
              l
            
          
        
      
    
    {\displaystyle f_{\mathrm {coll} }}
  
, given by kinetic theory the total
number of hard sphere collisions in a cell during a time 
  
    
      
        τ
      
    
    {\displaystyle \tau }
  
 is

M
          
            
              c
              o
              l
              l
            
          
        
        =
        
          
            1
            2
          
        
        (
        
          N
          
            
              c
            
          
        
        −
        1
        )
        
          F
          
            N
          
        
        
          f
          
            
              c
              o
              l
              l
            
          
        
        τ
        =
        
          
            
              
                N
                
                  
                    c
                  
                
              
              (
              
                N
                
                  
                    c
                  
                
              
              −
              1
              )
              
                F
                
                  N
                
              
              π
              
                d
                
                  2
                
              
              ⟨
              
                v
                
                  
                    r
                  
                
              
              ⟩
              τ
            
            
              2
              
                V
                
                  
                    c
                  
                
              
            
          
        
      
    
    {\displaystyle M_{\mathrm {coll} }={1 \over 2}(N_{\mathrm {c} }-1)F_{N}f_{\mathrm {coll} }\tau ={{N_{\mathrm {c} }(N_{\mathrm {c} }-1)F_{N}\pi d^{2}\langle v_{\mathrm {r} }\rangle \tau } \over {2V_{\mathrm {c} }}}}

where 
  
    
      
        d
      
    
    {\displaystyle d}
  
 is the particle diameter and 
  
    
      
        
          V
          
            
              c
            
          
        
      
    
    {\displaystyle V_{\mathrm {c} }}
  
 is the volume of the cell. 
Since collision candidates go through a rejection sampling procedure 
the ratio of total accepted to total candidates for hard sphere particles is

M
                
                  
                    c
                    o
                    l
                    l
                  
                
              
            
            
              
                M
                
                  
                    c
                    a
                    n
                    d
                  
                
              
            
          
        
        =
        
          
            
              ⟨
              
                v
                
                  
                    r
                  
                
              
              ⟩
            
            
              
                v
                
                  
                    r
                  
                
                
                  max
                
              
            
          
        
      
    
    {\displaystyle {{M_{\mathrm {coll} }} \over {M_{\mathrm {cand} }}}={{\langle v_{\mathrm {r} }\rangle } \over {v_{\mathrm {r} }^{\max }}}}

The number of collision candidates selected in a cell over a time step 
  
    
      
        τ
      
    
    {\displaystyle \tau }
  
 is

M
          
            
              c
              a
              n
              d
            
          
        
        =
        
          
            
              
                N
                
                  
                    c
                  
                
              
              (
              
                N
                
                  
                    c
                  
                
              
              −
              1
              )
              
                F
                
                  N
                
              
              π
              
                d
                
                  2
                
              
              
                v
                
                  
                    r
                  
                
                
                  max
                
              
              τ
            
            
              2
              
                V
                
                  
                    c
                  
                
              
            
          
        
      
    
    {\displaystyle M_{\mathrm {cand} }={{N_{\mathrm {c} }(N_{\mathrm {c} }-1)F_{N}\pi d^{2}v_{\mathrm {r} }^{\max }\tau } \over {2V_{\mathrm {c} }}}}

This approach for determining the number of collisions is known as the <i>No-Time-Counter</i> (NTC) method.
If 
  
    
      
        
          v
          
            
              r
            
          
          
            max
          
        
      
    
    {\displaystyle v_{\mathrm {r} }^{\max }}
  
 is set excessively high then the algorithm processes the same number of collisions (on average)
but the simulation is inefficient because many candidates are rejected.
</p><p>An alternative, more accurate and time-efficient algorithm is <i>Majorant Frequency </i> (MF) method, proposed by Mikhail Ivanov and Sergey Rogasinsky in 1988.<a class="footnote-ref" id="fnref:7" href="#fn:7"><sup>7</sup></a>
</p>

<h2 id="external-links">External links</h2>
<ul><li><a href="http://gab.com.au/">Direct Simulation Monte Carlo Method: Visual Simulation Programs created by GA Bird</a>.</li>
<li><a href="http://www.simba.us/misc/dsmc/dsmca.html">DSMC Demo Applet</a> by Greg Khanlarov</li>
<li><a href="http://homepage.univie.ac.at/franz.vesely/cp_tut/nol2h/new/c8hd_s4dsm.html">Course material on DSMC</a> (part of Computational Physics tutorial by Franz J. Vesely, University of Vienna)</li>
<li><a href="https://web.archive.org/web/20110717092648/https://www.ipam.ucla.edu/schedule.aspx?pc=kttut">Course material on DSMC and recent developments</a> (given at IPAM UCLA by Lorenzo Pareschi, University of Ferrara)</li></ul>

<h2 id="references">References</h2>

<ol>
<li id="fn:1"><p>Bird, G. A (1963). "Approach to Translational Equilibrium in a Rigid Sphere Gas". Physics of Fluids. 6 (10): 1518–1519. Bibcode:1963PhFl....6.1518B. doi:10.1063/1.1710976. <a href="/wiki/Bibcode_(identifier)" target="_blank">/wiki/Bibcode_(identifier)</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></p></li>
<li id="fn:2"><p>G. A. Bird, Molecular Gas Dynamics, Clarendon Press, Oxford (1976)[page needed] <a href="/wiki/Wikipedia:Citing_sources" target="_blank">/wiki/Wikipedia:Citing_sources</a> <a href="#fnref:2" class="footnote-back-ref">↩</a></p></li>
<li id="fn:3"><p>G. A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Clarendon Press, Oxford (1994)[page needed] <a href="/wiki/Wikipedia:Citing_sources" target="_blank">/wiki/Wikipedia:Citing_sources</a> <a href="#fnref:3" class="footnote-back-ref">↩</a></p></li>
<li id="fn:4"><p>Tsien, Hsue-Shen (1946). "Superaerodynamics, Mechanics of Rarefied Gases". Journal of the Aeronautical Sciences. 13 (12): 653–64. doi:10.2514/8.11476. <a href="/wiki/Doi_(identifier)" target="_blank">/wiki/Doi_(identifier)</a> <a href="#fnref:4" class="footnote-back-ref">↩</a></p></li>
<li id="fn:5"><p>M. N. Macrossan, 
'Scaling Parameters for Hypersonic Flow: Correlation of Sphere Drag Data'. In: M. S. Ivanov and A. K. Rebrov, 25th International Symposium on Rarefied Gas Dynamics, Siberian Division of the Russian Academy of Sciences, p.759 (2007). <a href="http://espace.library.uq.edu.au/view.php?pid=UQ:7959" target="_blank">http://espace.library.uq.edu.au/view.php?pid=UQ:7959</a> <a href="#fnref:5" class="footnote-back-ref">↩</a></p></li>
<li id="fn:6"><p>Roohi, E.; Stefanov, S. (2016). "Collision partner selection schemes in DSMC: From micro/nano flows to hypersonic flows". Physics Reports. 656 (1): 1–38. Bibcode:2016PhR...656....1R. doi:10.1016/j.physrep.2016.08.002. <a href="/wiki/Bibcode_(identifier)" target="_blank">/wiki/Bibcode_(identifier)</a> <a href="#fnref:6" class="footnote-back-ref">↩</a></p></li>
<li id="fn:7"><p>Ivanov, M. S.; Rogasinsky, S. V. (1988). "Analysis of numerical techniques of the direct simulation Monte Carlo method in the rarefied gas dynamics". Soviet Journal of Numerical Analysis and Mathematical Modelling. 3: 453–465. <a href="#fnref:7" class="footnote-back-ref">↩</a></p></li>
</ol>

Direct simulation Monte Carlo open-in-new

Direct simulation Monte Carlo