Particle tracking algorithms

flotus1 · July 17, 2013, 06:28

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:

$\frac{\partial u_\text{particle}}{\partial t}=C_D^* \left( u_\text{fluid} - u_\text{particle} \right)$ where $\frac{\partial x_\text{particle}}{\partial t} = u_\text{particle}$

Now solving this with any of the standard explicit schemes to obtain the particle trajectories (Euler, Heun, Runge-Kutta...) 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?

agd · July 17, 2013, 11:04

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.

triple_r · July 17, 2013, 11:55

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 semi-analytically 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.

FMDenaro · July 17, 2013, 12:15

I suppose that for very small dimensions you can treat directly the Lagrangian equation since inertia is disregardable...

flotus1 · July 17, 2013, 13:39

Thanks for your answers.

Quote:

Originally Posted by agd

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.

Actually, I dont need the flow field to be coupled to the particles in a first step. The particle tracking is more like a post-processing operation on the constant flow field.

Quote:

Originally Posted by triple_r

You might even be able to integrate it semi-analytically 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.

Correct me if I am wrong: Assuming a constant drag force would render any method first order accurate? If not, how to obtain a suitable predictor for the drag force?
The method implemented so far is 4th order Runge-Kutta, and I would have a hard time explaining why my "improved" method has lower accuracy.

Quote:

Originally Posted by FMDenaro

I suppose that for very small dimensions you can treat directly the Lagrangian equation since inertia is disregardable...

You are right.
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 so-called Euler-Maruyama 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 Euler-Maruyama scheme?

FMDenaro · July 17, 2013, 14:10

I know something is treated in this paper and its references
http://arxiv.org/pdf/1004.5598.pdf

triple_r · July 17, 2013, 14:45

Quote:

Originally Posted by flotus1

Correct me if I am wrong: Assuming a constant drag force would render any method first order accurate? If not, how to obtain a suitable predictor for the drag force?
The method implemented so far is 4th order Runge-Kutta, and I would have a hard time explaining why my "improved" method has lower accuracy.

I thought the issue was stability and not accuracy, my bad. That was just an example for how to get an stable solution. You can make it more accurate if you assume a known form for the drag coefficient, for example assuming spherical particles and Stokes flow, then you can integrate the equation without fixing the drag coefficient. The only thing that is being fixed is going to be the fluid velocity in that time step. So this might be even considered first-order depending on how smooth the velocity field is. To get even better solutions you might want to use a predictor-corrector approach: taking a step assuming constant u_f, using the u_f in the new location and the u_f from last location to assume a functionality for fluid velocity (like linear or constant but equal to the average value) and integrate from the previous location once more. This is going to get harder and harder the more terms you add to the equation of motion.

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 Euler-Maruyama, 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 Runge-Kutta variations for SDEs.

flotus1 · July 18, 2013, 03:40

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 non-constant 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 predictor-corrector 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 Euler-Maruyama according to the manual. I wonder why...

FMDenaro · July 18, 2013, 04:45

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 time-step 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....

flotus1 · July 18, 2013, 06:23

Quote:

Originally Posted by FMDenaro

Anyway, I don't think you need an accuracy higher for the particles integration than for the fluid solution....

That is a good point! I will keep this in mind.

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.

triple_r · July 18, 2013, 11:06

If you are going to use E-M 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.

SuperPahan · July 24, 2013, 02:53

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

flotus1 · July 24, 2013, 03:09

Please open a new thread in the Fluent section of the forum.
Hijacking other threads is not very polite.

flotus1 · July 30, 2013, 08:12

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

$\frac{\partial u_p}{\partial t} = C_D^*(u_f-u_p) + n_b$

The numerical method described by Li and Ahmadi (Aerosol science and Technology 16:209-226 (1992)) models the Brownian force
$n_b = G \sqrt{\frac{\pi S_0}{\Delta t}}$
where G is a normally distributed random variable with zero mean and unit variance and
$S_0 = \frac{216 \nu k T}{\pi^2 \rho d^5 S^2 C_c}$
Is the factor containing all the physics of the brownian motion.
The factor in the square-root 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.

flotus1 · July 31, 2013, 05:47

I made a boo-boo in the implementation of the implicit method

Now the results are reasonable.

But still I have no clue how

in the above formulas can be derived. If someone could shed light on this topic, I would be very pleased.

triple_r · August 1, 2013, 15:17

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.

flotus1 · August 1, 2013, 16:47

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.

$\sigma^2_\text{position} = 2Dt$

With my limited math for engineers I can go one step further to conclude that the variance for the discretized velocity has to be

$\sigma^2_\text{velocity} = 2D \Delta t$

since the sum of variances of the velocity yields the variance of the displacement (at least for non-correlated random processes)
For example i have n steps of size $\Delta t=\frac{t}{n}$ for the velocity, each with a variance $2D \Delta t$ . When summed up to get the displacement, the variance is $n 2 D \Delta t = 2Dt$ .
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.

leflix · August 18, 2013, 10:56

Quote:

Originally Posted by flotus1

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:

$\frac{\partial u_\text{particle}}{\partial t}=C_D^* \left( u_\text{fluid} - u_\text{particle} \right)$ where $\frac{\partial x_\text{particle}}{\partial t} = u_\text{particle}$

first of all I'm doubtfull about your equation. Why don't you account for gravity force too?
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.

flotus1 · August 18, 2013, 11:52

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.

leflix · August 18, 2013, 12:38

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.

July 17, 2013, 06:28	Particle tracking algorithms	#1
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,399 Rep Power: 46	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: $\frac{\partial u_\text{particle}}{\partial t}=C_D^* \left( u_\text{fluid} - u_\text{particle} \right)$ where $\frac{\partial x_\text{particle}}{\partial t} = u_\text{particle}$ Now solving this with any of the standard explicit schemes to obtain the particle trajectories (Euler, Heun, Runge-Kutta...) 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:55		#3
triple_r Senior Member Reza Join Date: Mar 2009 Location: Appleton, WI Posts: 116 Rep Power: 17	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 semi-analytically 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*

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
dispersion model with lagragian particle tracking model for incompressible flows	eelcovv	OpenFOAM Running, Solving & CFD	54	April 10, 2018 09:36
Blood Damage Modelling via Particle Tracking in a Centrifugal Heart Pump	scatman	CFX	7	January 8, 2018 00: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 21:11
Particle Tracking for ion	Jun	CFX	2	August 31, 2010 08:19

July 17, 2013, 11:04		#2
agd Senior Member Join Date: Jul 2009 Posts: 354 Rep Power: 18	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, 12:15		#4
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	I suppose that for very small dimensions you can treat directly the Lagrangian equation since inertia is disregardable...

July 17, 2013, 14:10		#6
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	I know something is treated in this paper and its references http://arxiv.org/pdf/1004.5598.pdf

July 18, 2013, 03:40		#8
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,399 Rep Power: 46	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 non-constant 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 predictor-corrector 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 Euler-Maruyama according to the manual. I wonder why...

July 18, 2013, 04:45		#9
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	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 time-step 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, 11:06		#11
triple_r Senior Member Reza Join Date: Mar 2009 Location: Appleton, WI Posts: 116 Rep Power: 17	If you are going to use E-M 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
SuperPahan New Member Pavel Join Date: Aug 2010 Location: St-Petersburg Posts: 10 Rep Power: 15	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
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,399 Rep Power: 46	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
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,399 Rep Power: 46	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 $\frac{\partial u_p}{\partial t} = C_D^*(u_f-u_p) + n_b$ The numerical method described by Li and Ahmadi (Aerosol science and Technology 16:209-226 (1992)) models the Brownian force $n_b = G \sqrt{\frac{\pi S_0}{\Delta t}}$ where G is a normally distributed random variable with zero mean and unit variance and $S_0 = \frac{216 \nu k T}{\pi^2 \rho d^5 S^2 C_c}$ Is the factor containing all the physics of the brownian motion. The factor in the square-root 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
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,399 Rep Power: 46	I made a boo-boo in the implementation of the implicit method Now the results are reasonable. But still I have no clue how in the above formulas can be derived. If someone could shed light on this topic, I would be very pleased.

August 1, 2013, 15:17		#16
triple_r Senior Member Reza Join Date: Mar 2009 Location: Appleton, WI Posts: 116 Rep Power: 17	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
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,399 Rep Power: 46	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. $\sigma^2_\text{position} = 2Dt$ With my limited math for engineers I can go one step further to conclude that the variance for the discretized velocity has to be $\sigma^2_\text{velocity} = 2D \Delta t$ since the sum of variances of the velocity yields the variance of the displacement (at least for non-correlated random processes) For example i have n steps of size $\Delta t=\frac{t}{n}$ for the velocity, each with a variance $2D \Delta t$ . When summed up to get the displacement, the variance is $n 2 D \Delta t = 2Dt$ . 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, 11:52		#19
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,399 Rep Power: 46	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
leflix Senior Member Join Date: Aug 2011 Posts: 272 Rep Power: 15	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.