
[Sponsors] 
July 17, 2013, 06:28 
Particle tracking algorithms

#1 
Senior Member
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,306
Rep Power: 22 
With my first naive approach to track the motion of particles in a fluid, I ran into the following problem.
Assuming a particle is only affected by drag forces, I write its equation of motion as: where Now solving this with any of the standard explicit schemes to obtain the particle trajectories (Euler, Heun, RungeKutta...) imposes a stability limit on the maximum timestep allowed. As the particles get smaller, increases and the largest stable timestep decreases. For particles smaller than a micrometer, the timestep gets unreasonably small for my application. How is this issue dealt with in general? 

July 17, 2013, 11:04 

#2 
Senior Member
Join Date: Jul 2009
Posts: 238
Rep Power: 11 
Can you treat the integration of the particle equation of motion as a series of substeps within one CFD timestep? This is equivalent to treating the flowfield as fixed during the integration of the particle equation during those substeps, but for very small particles (assuming the particle loading is small) the loss of coupling back to the flowfield should be negligible.


July 17, 2013, 11:55 

#3 
Senior Member
Reza
Join Date: Mar 2009
Location: Appleton, WI
Posts: 116
Rep Power: 9 
How about trying an implicit method? or a mix of implicit and explicit (implicit in time derivative, but explicit in drag coefficient for example). You might even be able to integrate it semianalytically if you know the form of the drag force and if you can assume that drag force doesn't change within a time step (from one location to the next).
For example, in the form that you have written the equation of motion, if C_d and u_fluid can be assumed to be constant in a single time step, then the solution will be (from time step n to time step n+1): u_p(n+1) = u_f(n) + (u_p(n)  u_f(n)) * exp(C_d(n)*dt) Another way is to use a higher order explicit method, or even better an adaptive high order method. Last edited by triple_r; July 17, 2013 at 12:24. Reason: forgot the minus sign :p 

July 17, 2013, 12:15 

#4 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,694
Rep Power: 33 
I suppose that for very small dimensions you can treat directly the Lagrangian equation since inertia is disregardable...


July 17, 2013, 13:39 

#5  
Senior Member
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,306
Rep Power: 22 
Thanks for your answers.
Quote:
Quote:
The method implemented so far is 4th order RungeKutta, and I would have a hard time explaining why my "improved" method has lower accuracy. Quote:
But the point is that I need to add a force term for Brownian motion on the right hand side later. And the only approach I found applied to CFD is the socalled EulerMaruyama method, which is based on the explicit Euler method and includes calculation of the drag force. Has anyone ever implemented particles with brownian motion without the EulerMaruyama scheme? 

July 17, 2013, 14:10 

#6 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,694
Rep Power: 33 
I know something is treated in this paper and its references
http://arxiv.org/pdf/1004.5598.pdf 

July 17, 2013, 14:45 

#7  
Senior Member
Reza
Join Date: Mar 2009
Location: Appleton, WI
Posts: 116
Rep Power: 9 
Quote:
Again, this is only if you want to use analytical integration. You don't have to though. Now that you are using 4th order RK and getting into stability problems, have you tried using a higher order method? 6th or 8th maybe? or even better, a variable order RK or adaptive RK? Take a look at Numerical Recipes if you want to implement the methods yourself and not use MATLAB for example. Regarding the Brownian motion and EulerMaruyama, that is going to be a half order SDE solver, but you can find higher order SDE solvers to use as well, like Milstein (first order) or there should even be RungeKutta variations for SDEs. 

July 18, 2013, 03:40 

#8 
Senior Member
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,306
Rep Power: 22 
Its about both, accuracy and stability.
For the value of the drag coefficient, I can make the assumption of Stokes flow around spherical particles with some kind of Cunningham correction factor for high Knudsen numbers. Thus the drag coefficient remains constant during the entire simulation. Thus the only variables are the fluid velocity (which I would like to treat as nonconstant during a time step) and the particle velocity. As far as I understood, ALL explicit methods have a rather small area of stability. The maximum stable step size is higher for higher order Methods, but not significantly. Am I right? Sorry for the basic question. How far can I get with a predictorcorrector method concerning accuracy and stability? Are these still among the explicit methods? I have lots of papers dealing with higher order integration of stochastic differential equations (also RK). But I havent seen any of these methods applied to CFD. Fluent for example seems to use EulerMaruyama according to the manual. I wonder why... 

July 18, 2013, 04:45 

#9 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 2,694
Rep Power: 33 
To tell the truth, I have not clear what you need.... I can see these problems:
1) unknown drag for the particle 2) unresolved fluid velocity 3) accuracy in the interpolation of the velocity particle on the Eulerian grid. Problems 2) and 3) appear when you are not solving in DNS approach, otherwise you don't need any adding of Brownian motion and have to face only with 1). In DNS you have all scales resolved and the timestep for the integration is constrained by physical assumptions to be quite small If you are not in a DNS you solve either for URANS or LES. Such cases, resolving different velocities, require different approach as function of the inertia of the particle. Actually, transport of particles in URANS is quite questionable in my opinion... Anyway, I don't think you need an accuracy higher for the particles integration than for the fluid solution.... 

July 18, 2013, 06:23 

#10  
Senior Member
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,306
Rep Power: 22 
Quote:
The flow field is computed with a DNS (well actually the flow is steady and laminar due to low Reynolds number at high Knudsen numbers). Interpolation of the fluid velocity at arbitrary points should not be a mayor problem since the grid is cartesian. I dont know if I got your point right. I dont need the "Brownian" motion to mimic the movement of particles caused by unresolved eddies, but I want to model the real Brownian motion of the small particles. So what I need (or what would be favorite method) is an integration algorithm that can deal with "large" timesteps, at least large compared to the stability constraint of the explicit methods. And it should be at least second order accurate. And of course the integration of the stochastic part should be dealt with. But for this last point, I think I have to read a bit more about SDEs. Anyway, thanks for your valuable input so far. 

July 18, 2013, 11:06 

#11 
Senior Member
Reza
Join Date: Mar 2009
Location: Appleton, WI
Posts: 116
Rep Power: 9 
If you are going to use EM for SDE when you add Brownian motion, then I don't see the point in using a second order method for other terms.
I myself have only used EM, and my reason was its being easy to program, so I didn't have any reasons other than my laziness :p However, you might want to use higher order solvers if you want your over all integration to be of second order or higher. 

July 24, 2013, 02:53 

#12 
New Member

How can I apply the ring cone in fluent like in cfx with indication outer and inner radii of the particle injection ring?


July 24, 2013, 03:09 

#13 
Senior Member
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,306
Rep Power: 22 
Please open a new thread in the Fluent section of the forum.
Hijacking other threads is not very polite. 

July 30, 2013, 08:12 

#14 
Senior Member
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,306
Rep Power: 22 
I have an additional question regarding the approximation of a Brownian motion with the Euler Maruyama method.
The particle equation of motion in one dimension reads The numerical method described by Li and Ahmadi (Aerosol science and Technology 16:209226 (1992)) models the Brownian force where G is a normally distributed random variable with zero mean and unit variance and Is the factor containing all the physics of the brownian motion. The factor in the squareroot of the brownian force is obviously chosen to produce the correct diffusion properties for the displacement of the particles. Now my problem is: how exactly can this factor be determined? This step appears to me as crucial for the implementation of any numerical scheme modeling brownian motion this way. For example: if I use the factor given above with an explicit Euler integration to obtain the particle velocity, everything works fine. But the same factor used with an implicit euler scheme does not give the correct diffusion coefficient for the displacement. So how can the factor (the variance of the brownian force) be determined for arbitrary integration schemes of the equation of motion? It would already help if i knew how the factor is determined for the Euler Maruyama scheme. Starting from there, I should be able to figure out the procedure for other schemes. 

July 31, 2013, 05:47 

#15 
Senior Member
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,306
Rep Power: 22 

August 1, 2013, 15:17 

#16 
Senior Member
Reza
Join Date: Mar 2009
Location: Appleton, WI
Posts: 116
Rep Power: 9 
I didn't have access to paper you mentioned, so I don't know how they have derived, but my guess is it comes from Einstein's equation for the PDF for Brownian motion of a particle:
that says what is the probability of a particle being at location x at time t, granted it was at location at time 0. He derives a differential equation for P, and then solves it. Later on he gives a relationship between the friction factor and the Brownian diffusion coefficient. My guess is this equation is from that relationship when they substitute the friction factor for a specific shape (for example friction factor for a sphere in Stokes flow). If I remember the title correctly, it was "on the theory of Brownian motion" or something similar. Because it is by Einstein, even translated versions should be downloadable for free, just google for it. 

August 1, 2013, 16:47 

#17 
Senior Member
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,306
Rep Power: 22 
They dont derive the formula in the paper, thats my problem... However, here is a link to the paper
I know the factor is chosen to yield the correct variance for the displacement of a spherical particle with a Stokes law of friction corrected by a Cunningham factor for large Knudsen numbers. You are right with this assumption, apart from the Cunningham correction that is Einsteins formula for the diffusion. With my limited math for engineers I can go one step further to conclude that the variance for the discretized velocity has to be since the sum of variances of the velocity yields the variance of the displacement (at least for noncorrelated random processes) For example i have n steps of size for the velocity, each with a variance . When summed up to get the displacement, the variance is . So far, so good It should be straightforward to carry out the same step to get the variance for the acceleration, from which the factor for the stochastic term can be derived. But from the discretized form of the acceleration, I am unable to determine the variance due to a lack of mathematic skills. I left my notes at the office, so I will post the exact formula here tomorrow. Anyway, thanks for your interest in my problem so far. 

August 18, 2013, 10:56 

#18  
Senior Member
Join Date: Aug 2011
Posts: 251
Rep Power: 8 
Quote:
you should rather use: m dup/dt = Fdrag + mg dx/dt=up check that here the derivative is d and not partial derivative (round d)... In this case you use a one way coupling between fluid and particles. If the fluid velocity uf is computed accounting for the presence of particles then you will get two way coupling. But I'm not sure about what you are trying to do. Indeed you can also compute streaklines. In this case particles are seen as tracers. Then you just solve dx/dt = uf Secondly using explicit or implicit methods to solve your ordinary differential equations does not matter. Explicit method are nice for this problem. 2nd order Runge Kutta is ok if you are looking for accuracy if not you can even use 1rst order Euler method. It's enough. In every cases you will need to evaluate uf at the particle location. Bilinear interpolation is requested. 

August 18, 2013, 11:52 

#19 
Senior Member
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,306
Rep Power: 22 
You seem to have further knowledge about this topic.
Would you mind commenting on the other 16 posts after the initial one? Like I said, I assume the particle is only affected by drag forces (and Brownian forces) because gravity is negligible in my case. 

August 18, 2013, 12:38 

#20 
Senior Member
Join Date: Aug 2011
Posts: 251
Rep Power: 8 
what is your exact question? because I'm not sure what is your 16th post.
Unless you are in microgravity or zero condition or if m is very small, I do not see how you can neglect mg compare to Fdrag.... In any case the equation should be m du/dt = Fdrag and not du/dt=Fdrag In the case where m is very small it would give du/dt =Fdrag/m which could be quite high...pay attention to this. 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
dispersion model with lagragian particle tracking model for incompressible flows  eelcovv  OpenFOAM Running, Solving & CFD  52  July 24, 2016 19:59 
Ubuntu 12.10 + openfoam2.2.0 ==> paraview error message  peteryuan  OpenFOAM Installation  6  August 18, 2013 18:00 
Particle tracking prob, urgent.  sakurabogoda  CFX  1  March 11, 2013 22:11 
Blood Damage Modelling via Particle Tracking in a Centrifugal Heart Pump  scatman  CFX  5  May 5, 2011 07:23 
Particle Tracking for ion  Jun  CFX  2  August 31, 2010 08:19 