CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Particle tracking algorithms

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   July 17, 2013, 06:28
Default Particle tracking algorithms
  #1
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
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, C_D^* 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?
flotus1 is offline   Reply With Quote

Old   July 17, 2013, 11:04
Default
  #2
agd
Senior Member
 
Join Date: Jul 2009
Posts: 192
Rep Power: 8
agd is on a distinguished road
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.
agd is offline   Reply With Quote

Old   July 17, 2013, 11:55
Default
  #3
Member
 
Reza
Join Date: Mar 2009
Location: Appleton, WI
Posts: 99
Rep Power: 8
triple_r is on a distinguished road
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
triple_r is offline   Reply With Quote

Old   July 17, 2013, 12:15
Default
  #4
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,645
Rep Power: 23
FMDenaro will become famous soon enough
I suppose that for very small dimensions you can treat directly the Lagrangian equation since inertia is disregardable...
FMDenaro is offline   Reply With Quote

Old   July 17, 2013, 13:39
Default
  #5
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
Thanks for your answers.

Quote:
Originally Posted by agd View Post
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 View Post
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 View Post
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?
flotus1 is offline   Reply With Quote

Old   July 17, 2013, 14:10
Default
  #6
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,645
Rep Power: 23
FMDenaro will become famous soon enough
I know something is treated in this paper and its references
http://arxiv.org/pdf/1004.5598.pdf
FMDenaro is offline   Reply With Quote

Old   July 17, 2013, 14:45
Default
  #7
Member
 
Reza
Join Date: Mar 2009
Location: Appleton, WI
Posts: 99
Rep Power: 8
triple_r is on a distinguished road
Quote:
Originally Posted by flotus1 View Post

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.
triple_r is offline   Reply With Quote

Old   July 18, 2013, 03:40
Default
  #8
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
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...
flotus1 is offline   Reply With Quote

Old   July 18, 2013, 04:45
Default
  #9
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,645
Rep Power: 23
FMDenaro will become famous soon enough
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....
FMDenaro is offline   Reply With Quote

Old   July 18, 2013, 06:23
Default
  #10
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
Quote:
Originally Posted by FMDenaro View Post
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.
flotus1 is offline   Reply With Quote

Old   July 18, 2013, 11:06
Default
  #11
Member
 
Reza
Join Date: Mar 2009
Location: Appleton, WI
Posts: 99
Rep Power: 8
triple_r is on a distinguished road
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.
triple_r is offline   Reply With Quote

Old   July 24, 2013, 02:53
Default
  #12
New Member
 
Pavel
Join Date: Aug 2010
Location: St-Petersburg
Posts: 10
Rep Power: 6
SuperPahan is on a distinguished road
Send a message via Skype™ to SuperPahan
How can I apply the ring cone in fluent like in cfx with indication outer and inner radii of the particle injection ring?
SuperPahan is offline   Reply With Quote

Old   July 24, 2013, 03:09
Default
  #13
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
Please open a new thread in the Fluent section of the forum.
Hijacking other threads is not very polite.
flotus1 is offline   Reply With Quote

Old   July 30, 2013, 08:12
Default
  #14
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
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 is offline   Reply With Quote

Old   July 31, 2013, 05:47
Default
  #15
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
I made a boo-boo in the implementation of the implicit method
Now the results are reasonable.

But still I have no clue how S_0 in the above formulas can be derived. If someone could shed light on this topic, I would be very pleased.
flotus1 is offline   Reply With Quote

Old   August 1, 2013, 15:17
Default
  #16
Member
 
Reza
Join Date: Mar 2009
Location: Appleton, WI
Posts: 99
Rep Power: 8
triple_r is on a distinguished road
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:

P(x, t|x_0)

that says what is the probability of a particle being at location x at time t, granted it was at location x_0 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.
triple_r is offline   Reply With Quote

Old   August 1, 2013, 16:47
Default
  #17
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
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.
flotus1 is offline   Reply With Quote

Old   August 18, 2013, 10:56
Default
  #18
Senior Member
 
Join Date: Aug 2011
Posts: 251
Rep Power: 6
leflix is on a distinguished road
Quote:
Originally Posted by flotus1 View Post
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.
leflix is offline   Reply With Quote

Old   August 18, 2013, 11:52
Default
  #19
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
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.
flotus1 is offline   Reply With Quote

Old   August 18, 2013, 12:38
Default
  #20
Senior Member
 
Join Date: Aug 2011
Posts: 251
Rep Power: 6
leflix is on a distinguished road
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.
leflix is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
dispersion model with lagragian particle tracking model for incompressible flows eelcovv OpenFOAM Running, Solving & CFD 48 January 31, 2015 11:10
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


All times are GMT -4. The time now is 02:54.