PDF
Brownian Motion
Introduction
Particle tracing offers an attractive alternative to continuum-based numerical methods, such as the finite element method, for modeling species transport in strongly convecting flows because the particle tracing method is insensitive to the magnitude of the Péclet number. In most real systems, species transport has both a convective and diffusive component. Particle tracing can be used to solve purely convective motion, purely diffusive motion and anything in-between. Thus using a particle-based approach the entire spectrum of Péclet numbers can be handled without encountering the numerical instabilities associated with the continuum approach.
In this example the agreement between the continuum and particle-based numerical methods is verified in the case of purely diffusive motion.
Model Definition
The diffusion equation is solved in two different ways. First, the species concentration is computed using the Transport of Diluted Species interface, which uses a continuum model in which the concentration is discretized using a finite element mesh in the modeling domain. The equation governing the evolution of the concentration c (SI unit: mol/m3) in a stagnant background fluid (u = 0) is:
where the diffusion coefficient D (SI unit: m2/s) is defined as
where
•
kB = 1.3806488 × 10−23 J/K is Boltzmann’s constant,
•
T (SI unit: K) is the absolute fluid temperature,
•
μ (SI unit: Pa·s) is the fluid viscosity, and
•
rp (SI unit: m) is the particle radius.
The initial concentration is given by a Dirac delta function at the origin,
That is, the initial concentration is infinitely large at the origin and zero everywhere else, such that the surface integral over any region containing the origin is unity. Because an initial condition that is infinitely large at a point is impractical to model, the initial concentration is instead given a very large finite value in a small area surrounding the origin.
The model geometry consists of two concentric circles as shown in Figure 1. The initial concentration diffuses from the origin radially outward in all directions. After 100 seconds, some of the initial concentration has diffused from the inner circular domain to the outer domain. The transmission probability for diffusion from the inner domain to the outer domain is defined as:
where I denotes the inner domain and O the outer domain.
Figure 1: Model geometry. The concentration is initially a delta function at (0,0) and the particles are all released at (0,0) with an initial velocity of zero.
Diffusion can also be modeled using a particle-based approach. The combination of the Brownian Force and the Drag Force results in diffusion of particles from regions of high to low number density. The equations of motion are:
This is the Stokes drag law, which is appropriate when the relative Reynolds number of the particles in the fluid is small. Because the fluid is stagnant and the particle speeds are low, the Stokes drag law is applicable in this example. The Brownian force is given by
•
mp (SI unit: kg) is the particle mass,
•
dp (SI unit: m) is the particle diameter,
•
τp (SI unit: s) is the particle velocity response time,
•
v (SI unit: m/s) is the velocity of the particle,
•
u (SI unit: m/s) is the fluid velocity, which in this example is set to zero representing a stagnant background fluid,
•
Δt (SI unit: s) is the size of the time step taken by the solver, and
•
ζ (dimensionless) is a vector of independent, normally distributed random numbers with zero mean and unit standard deviation.
As explained in Ref. 1, independent values are chosen for all components of ζ. A different value of ζ is created for each particle, at each time step for each component of the Brownian force. The Brownian force leads to spreading of particles from regions of high particle density to low density.
Initially, 5000 particles are released from the origin with initial velocity components all zero. The transmission probability from the inner to the outer domain is computed by counting the number of particles in the outer domain and dividing it by the total number of particles.
Results and Discussion
The computed transmission probability for the two different methods is shown in Table 1. Since the Brownian force uses pseudorandom number generators, the problem is solved five times, each time with a different seed. In each case the transmission probability is slightly different but all cases agree with the result from solving the diffusion equation.
Table 1: table of results for the computed transmission probability.
method
transmission probability
Diffusion Equation
0.23645
Particle Tracing 1
0.24460
Particle Tracing 2
0.23080
Particle Tracing 3
0.24860
Particle Tracing 4
0.23760
Particle Tracing 5
0.23840
Figure 2 plots the location of the particles at the final solution time for 4 different runs.
Figure 2: Plot of the particle location after 100 seconds. For each run, different random numbers were generated.
It is clear from these results that diffusive processes can be modeled using a particle-based approach. Furthermore, if a significant nonzero background fluid velocity were applied, the particle based approach would remain numerically stable.
Reference
1. M. Kim and A.L. Zydney, “Effect of Electrostatic, Hydrodynamic, and Brownian Forces on Particle Trajectories and Sieving in Normal Flow Filtration,” J. Colloid Interface Sci., vol. 269, pp. 425–431, 2004.
Application Library path: Particle_Tracing_Module/Verification_Examples/brownian_motion
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D.
2
In the Select Physics tree, select Chemical Species Transport > Transport of Diluted Species (tds).
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies > Time Dependent.
6
Click  Done.
Geometry 1
Start by adding some definitions for the geometry and physical properties of the background fluid.
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
In the table, enter the following settings:
Name
Expression
Value
Description
router
0.0005
5E-4
Outer radius
rinner
0.00025
2.5E-4
Inner radius
rp
1E-7[m]
1E-7 m
Particle radius
T
300[K]
300 K
Temperature
eta
2E-5[Pa*s]
2E-5 Pa·s
Fluid viscosity
D
k_B_const*T/(6*pi*eta*rp)
1.0987E-10 m²/s
Diffusivity
ds
1
1
Input to random number generator
Geometry 1
Circle 1 (c1)
1
In the Model Builder window, expand the Component 1 (comp1) > Geometry 1 node.
2
Right-click Geometry 1 and choose Circle.
3
In the Settings window for Circle, locate the Size and Shape section.
4
In the Radius text field, type router.
Circle 2 (c2)
1
In the Geometry toolbar, click  Circle.
2
In the Settings window for Circle, locate the Size and Shape section.
3
In the Radius text field, type rinner.
Point 1 (pt1)
1
In the Geometry toolbar, click  Point.
2
In the Settings window for Point, click  Build All Objects.
3
Click the  Zoom Extents button in the Graphics toolbar. The geometry should look like Figure 1.
Definitions
Construct an expression for the initial concentration, which is a smoothed delta function.
Variables 1
1
In the Model Builder window, expand the Component 1 (comp1) > Definitions node.
2
Right-click Definitions and choose Local Variables.
3
In the Settings window for Variables, locate the Variables section.
4
In the table, enter the following settings:
Name
Expression
Unit
Description
smooth
2E-7
 
Smoothing distance
xd
x[1/m]
 
x-coordinate
yd
y[1/m]
 
y-coordinate
c0
1
 
Peak initial concentration
c_init
2*c0*(1-flc2hs(xd^2+yd^2-smooth^2,5e-11))
 
Initial concentration
Define a pair of Integration component couplings so that the fraction of the concentration that diffuses from the inner domain to the outer domain can be computed.
Integration 1 (intop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
Select Domain 1 only.
Integration 2 (intop2)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
Select Domain 2 only.
Transport of Diluted Species (tds)
1
In the Model Builder window, under Component 1 (comp1) click Transport of Diluted Species (tds).
2
In the Settings window for Transport of Diluted Species, locate the Transport Mechanisms section.
3
Clear the Convection checkbox.
Concentration 1
1
In the Physics toolbar, click  Boundaries and choose Concentration.
2
Select Boundaries 1, 2, 5, and 8 only.
3
In the Settings window for Concentration, locate the Concentration section.
4
Select the Species c checkbox.
Fluid 1
1
In the Model Builder window, click Fluid 1.
2
In the Settings window for Fluid, locate the Diffusion section.
3
In the Dc text field, type D.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the c text field, type c_init.
Mesh 1
Use a very fine mesh at the point (0, 0).
1
In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2
In the Settings window for Mesh, locate the Physics-Controlled Mesh section.
3
From the Element size list, choose Extra fine.
Size 1
1
In the Mesh toolbar, click  Sizing and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Point.
4
Select Point 5 only.
5
Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section.
7
Select the Maximum element size checkbox. In the associated text field, type 6.7E-5*0.01.
Free Triangular 1
1
In the Mesh toolbar, click  Free Triangular.
2
In the Settings window for Free Triangular, click  Build All.
Study 1
Step 1: Time Dependent
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type 0 100.
4
From the Tolerance list, choose User controlled.
5
In the Relative tolerance text field, type 1E-4.
6
In the Study toolbar, click  Compute.
Results
Use the Global Evaluation feature to compute the fraction of the total concentration that diffused from the inner domain to the outer domain.
Global Evaluation 1
1
In the Results toolbar, click  Global Evaluation.
2
In the Settings window for Global Evaluation, locate the Data section.
3
From the Time selection list, choose Last.
4
Locate the Expressions section. In the table, enter the following settings:
Expression
Unit
Description
intop1(c)/(intop1(c)+intop2(c))
1
 
5
Clicknext to  Evaluate, then choose New Table.
Component 1 (comp1)
Now solve the same problem using a particle-based approach.
Add Physics
1
In the Home toolbar, click  Add Physics to open the Add Physics window.
2
Go to the Add Physics window.
3
In the tree, select Fluid Flow > Particle Tracing > Particle Tracing for Fluid Flow (fpt).
4
Find the Physics interfaces in study subsection. In the table, clear the Solve checkbox for Study 1.
5
Click the Add to Component 1 button in the window toolbar.
6
In the Home toolbar, click  Add Physics to close the Add Physics window.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Physics interfaces in study subsection. In the table, clear the Solve checkbox for Transport of Diluted Species (tds).
4
Find the Studies subsection. In the Select Study tree, select General Studies > Time Dependent.
5
Click the Add Study button in the window toolbar.
6
In the Home toolbar, click  Add Study to close the Add Study window.
Root
1
Click the  Show More Options button in the Model Builder toolbar.
2
In the Show More Options dialog, in the tree, select the checkbox for the node Physics > Advanced Physics Options.
3
Click OK.
Particle Tracing for Fluid Flow (fpt)
1
In the Model Builder window, under Component 1 (comp1) click Particle Tracing for Fluid Flow (fpt).
2
In the Settings window for Particle Tracing for Fluid Flow, click to expand the Advanced Settings section.
3
From the Wall accuracy order list, choose 1.
4
From the Arguments for random number generation list, choose User defined.
Particle Properties 1
1
In the Model Builder window, under Component 1 (comp1) > Particle Tracing for Fluid Flow (fpt) click Particle Properties 1.
2
In the Settings window for Particle Properties, locate the Particle Properties section.
3
In the dp text field, type 2*rp.
4
From the ρp list, choose User defined.
Drag Force 1
1
In the Physics toolbar, click  Domains and choose Drag Force.
2
In the Settings window for Drag Force, locate the Domain Selection section.
3
From the Selection list, choose All domains.
4
Locate the Drag Force section. From the μ list, choose User defined. In the associated text field, type eta.
Brownian Force 1
1
In the Physics toolbar, click  Domains and choose Brownian Force.
2
In the Settings window for Brownian Force, locate the Domain Selection section.
3
From the Selection list, choose All domains.
4
Locate the Model Input section. In the T text field, type T.
5
Locate the Brownian Force section. From the μ list, choose User defined. In the associated text field, type eta.
6
Locate the Advanced Settings section. In the i text field, type ds.
Release from Grid 1
Release 5000 particles at the origin with an initial velocity of zero.
1
In the Physics toolbar, click  Global and choose Release from Grid.
2
In the Settings window for Release from Grid, locate the Initial Velocity section.
3
From the Initial velocity list, choose Constant speed, spherical.
4
In the v0 text field, type 0.
5
In the Nvel text field, type 5000.
Particle Counter 1
1
In the Physics toolbar, click  Domains and choose Particle Counter.
2
Select Domain 1 only.
3
In the Settings window for Particle Counter, locate the Particle Counter section.
4
From the Release feature list, choose Release from Grid 1.
The Brownian force depends on the time step taken by the solver. The default tolerances are very strict for the particle tracing interfaces. When including forces with random components such as the Brownian force the tolerances need to be relaxed, otherwise the solver will take very small time steps and the model will take a long time to solve.
Study 2
Step 1: Time Dependent
1
In the Model Builder window, under Study 2 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type 0 100.
4
From the Tolerance list, choose User controlled.
5
In the Relative tolerance text field, type 1E-2.
6
In the Study toolbar, click  Compute.
Results
Global Evaluation 2
1
In the Results toolbar, click  Global Evaluation.
2
In the Settings window for Global Evaluation, locate the Data section.
3
From the Dataset list, choose Study 2/Solution 2 (sol2).
4
From the Time selection list, choose Last.
5
Click Replace Expression in the upper-right corner of the Expressions section. From the menu, choose Component 1 (comp1) > Particle Tracing for Fluid Flow > Particle Counter 1 > fpt.pcnt1.alpha - Transmission probability - 1.
6
Click  Evaluate.
Root
Finally, add another study and solve the same problem 4 times, with different random numbers generated for each run. A Parametric Sweep over the parameter ds is used to create unique random numbers for each run.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Physics interfaces in study subsection. In the table, clear the Solve checkbox for Transport of Diluted Species (tds).
4
Find the Studies subsection. In the Select Study tree, select General Studies > Time Dependent.
5
Click the Add Study button in the window toolbar.
6
In the Home toolbar, click  Add Study to close the Add Study window.
Study 3
Step 1: Time Dependent
1
In the Settings window for Time Dependent, locate the Study Settings section.
2
In the Output times text field, type 0 100.
3
From the Tolerance list, choose User controlled.
4
In the Relative tolerance text field, type 1E-2.
Parametric Sweep
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
Click  Add.
4
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
ds (Input to random number generator)
2 3 4 5
 
5
In the Study toolbar, click  Compute.
Results
Global Evaluation 3
1
In the Model Builder window, expand the Results > Datasets node.
2
Right-click Results > Derived Values and choose Global Evaluation.
3
In the Settings window for Global Evaluation, locate the Data section.
4
From the Dataset list, choose Study 3/Parametric Solutions 1 (sol4).
5
From the Time selection list, choose Last.
6
From the Table columns list, choose Time.
7
Click Replace Expression in the upper-right corner of the Expressions section. From the menu, choose Component 1 (comp1) > Particle Tracing for Fluid Flow > Particle Counter 1 > fpt.pcnt1.alpha - Transmission probability - 1.
8
Click the arrow next to the Evaluate button and select New table.
Particle Trajectories (fpt) 1
1
In the Model Builder window, under Results click Particle Trajectories (fpt) 1.
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Parameter value (ds) list, choose 2.
4
In the Particle Trajectories (fpt) 1 toolbar, click  Plot.
5
Click the  Zoom Extents button in the Graphics toolbar.
View the particle trajectories for other parameter values by selecting different options from the Parameter value (ds) list. Note that the distributions all look similar but that the particle positions are distinct in each plot. Four such plots are shown in Figure 2.