CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Main CFD Forum (https://www.cfd-online.com/Forums/main/)
-   -   Particle tracking algorithms (https://www.cfd-online.com/Forums/main/120910-particle-tracking-algorithms.html)

flotus1 July 17, 2013 07:28

Particle tracking algorithms
 
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?

agd July 17, 2013 12: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 12: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 13:15

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

flotus1 July 17, 2013 14:39

Thanks for your answers.

Quote:

Originally Posted by agd (Post 440343)
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 (Post 440354)
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 (Post 440356)
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 15:10

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

triple_r July 17, 2013 15:45

Quote:

Originally Posted by flotus1 (Post 440371)

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 04: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 05: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 07:23

Quote:

Originally Posted by FMDenaro (Post 440482)
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 12: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 03: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 04:09

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

flotus1 July 30, 2013 09: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 06:47

I made a boo-boo in the implementation of the implicit method :rolleyes:
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.

triple_r August 1, 2013 16: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:

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.

flotus1 August 1, 2013 17: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:D

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 11:56

Quote:

Originally Posted by flotus1 (Post 440270)
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 12: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 13: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.

flotus1 August 18, 2013 14:46

The question boiled down to: how is http://www.cfd-online.com/Forums/vbL...c9a3ef94-1.gif for the amplitude of the stochastic part of the particle equation of motion derived.

Dont worry, my equation is already divided by the particle mass and yes, the particle mass is very small.


All times are GMT -4. The time now is 13:43.