Brownian Dynamics of a Spherical Particle
in a Fluid Flow with no External Force
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :

This example was inspired from the Matlab program available at the following MIT OpenCourseWare web page:

http://ocw.mit.edu/OcwWeb/Chemical-Engineering/10-34Fall-2005/CourseHome/index.htm

Integration up to t=1 for 5000 times (this computation takes several minutes)

<< Statistics`ContinuousDistributions`

x[0] = 0;

ξ  = 1; ksp = 1; Di = 1;

Δt = 0.001;

Vf = 1;

i = 1;

For[j = 1, 
      j < 5001, {i = 1, 
        While[i < 
            1001, {x[i] = 
              x[i - 1] + 
                Vf Δt + (Sqrt[2 Di Δt]) 
                  Random[NormalDistribution[0, 1]], i = i + 1}], 
        X[j] = x[i - 1], j++, 
        If[Mod[j, 100] ==  0, Print[j]]}]

Probability distribution function compared to solution of convection/diffusion equation

<< Statistics`DataManipulation`

t = BinCounts[Table[X[j], {j, 1, 4999}], {-5, 5, 0.1}]

midpoints = Table[i + 0.05, {i, -5, 5 - 0.1, 0.1}]

su = Sum[t[[i]], {i, 1, Length[t]}]*0.1

plt1 = 
      ListPlot[Table[ {midpoints[[i]], t[[i]]/su}, {i, 1, 
            Length[t]}], PlotJoined ->  True,  
        PlotStyle  ->  RGBColor[1, 0, 0]]

tf = 1;

σ = Sqrt[2 Di tf]

plt2 = 
      Plot[Exp[-(x - Vf tf)^2/(2*σ^2)]/(σ Sqrt[
                2 Pi]), {x, -20, 20}, 
        PlotStyle ->  RGBColor[0, 0, 1], 
        PlotRange ->  All]

Show[plt1, plt2, PlotRange ->  All]

Integration up to t=3  for 5000 times (this computation takes several minutes)

For[j = 1, 
      j < 5001, {i = 1, 
        While[i < 
            3001, {x[i] = 
              x[i - 1] + 
                Vf Δt + (Sqrt[2 Di Δt]) 
                  Random[NormalDistribution[0, 1]], i = i + 1}], 
        X[j] = x[i - 1], j++, 
        If[Mod[j, 100] ==  0, Print[j]]}]

Probability distribution function compared to solution of convection/diffusion equation

t = BinCounts[Table[X[j], {j, 1, 4999}], {-10, 10, 0.1}]

midpoints = Table[i + 0.05, {i, -10, 10 - 0.1, 0.1}]

su = Sum[t[[i]], {i, 1, Length[t]}]*0.1

plt1 = 
      ListPlot[Table[ {midpoints[[i]], t[[i]]/su}, {i, 1, 
            Length[t]}], PlotJoined ->  True,  
        PlotStyle  ->  RGBColor[1, 0, 0]]

tf = 3;

σ = Sqrt[2 Di tf]

plt2 = 
      Plot[Exp[-(x - Vf tf)^2/(2*σ^2)]/(σ Sqrt[
                2 Pi]), {x, -20, 20}, 
        PlotStyle ->  RGBColor[0, 0, 1], 
        PlotRange ->  All]

Show[plt1, plt2, PlotRange ->  All]