Dear people from stackexchange physics.

I have been using SPH for quite a while to simulate free boundary flow, and just recently we tried to include some kind of (simple) chemical reaction in our code. We are still doing some tests, but one curious thing that we found was that we needed much smaller time steps than we usually used to be necessary.

I have never dealed with chemical reactions, but I don know that some of them fall on the category of 'Siff Equations', and as far as the wikipedia description goes, it seems to be my case, but I would like to ask if anyone tried anything similar and also found that their equations were stiff.

I never studied stiff equations, I just knew that they existed and recently read a bit more about them, but that's it. Last but not least, do anyone know how to check if a set of SPH equations is stiff or not? The wikipedia recomendation seemed at least complicated to calculate for SPH.


closed as off-topic by Kyle Kanos, Brandon Enright, jinawee, user10851, John Rennie Jan 10 '14 at 10:25

  • This question does not appear to be about physics within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 4
    $\begingroup$ This question appears to be off-topic because it is about programming/computational methods and not physics. $\endgroup$ – Kyle Kanos Jan 9 '14 at 16:02
  • 3
    $\begingroup$ You may want to try asking this at scicomp.stackexchange.com. $\endgroup$ – Kyle Kanos Jan 9 '14 at 16:03
  • 1
    $\begingroup$ As Kyle Kanos points out, this is really a computational question. For non-reacting flows the time step limitation of your numerical method comes from the speed of sound in the medium. Reaction rates are typically much faster than sound waves and hence you will end up with a "stiff" set of equations. One approach to dealing with this is operator splitting whereby you use an explicit step for the flow followed by an implicit one for the reaction. $\endgroup$ – SimpleLikeAnEgg Jan 9 '14 at 20:56

Stiff equations tend to use "implicit" methods. I am actually working on an implicit simulation (rigid body constraints) which is infinitely stiff! Implicit methods are basically guess and check: try to find a set of forces that satisfies future constraints. This is a well-developed field, multigrid conjugate-gradient based methods are among the state of the art.

So how would you use an implicit method? Take a pendulum for which we have a single anchor-spring-mass as an example. The spring is stiff, almost a rigid rod, and you don't care about the high-frequency shortening and lengthening vibrations of the spring. You do care about modeling the overall swinging motion, however. "Normal" (explicit) solvers need a time-step that is small enough to capture the spring's vibrations. For an implicit scheme, the step is set by the dynamics we care about: the swing period of the pendulum, which is much longer than the vibrational period of the spring. Thus implicit methods, although more expensive per step, allow us to gloss over high-frequency degrees of freedom that we don't care about (as long as they these DOF's are nearly linear, as is true in our case with small spring vibrations).

For the pendulum case it is easy to create a new model bases on an ideal rigid rod, eliminating the annoying high-frequency DOF. For more complex cases this is difficult to do, however.

Stiff equations in chemistry typically arise when reaction time-constants are small. In a first order decay reaction: dC/dt = -k*C, the time constant is 1/k. For non-linear complex reactions, you can look at small perturbations of concentration and see how quickly they grow/decay.

Another way is for diffusivity to be high. The time constant is ~ neighbor distance^2/(diffusivity). It shrinks rapidly with decreased particle size.

Chemical stiffness is an issue if these time-constants are smaller than the SPH time-step needed. In this case you can decouple the simulations: take an SPH step and take several smaller diffusion and/or reaction sub-steps, or use an implicit model for the diffusion/reaction. Implicit methods for diffusion work well since error tends to blurs out, as long as total mass is conserved.

If the only stiff aspect is reaction kinetics, here's a simple way to use "sub-steps": for each SPH step calculate the diffusion fluxes. Assume a constant flux into your particle over the entire SPH timestep. For your particle, calculate reaction kinetics using much smaller sub-steps until you found your concentrations at the end of the SPH timestep. This way you capture the fast kinetics without any extra expensive SPH calculations.

SPH simulations themselves get stiff when fluids become incompressible and/or resolution gets higher. The step size ~ (particle size)/(speed of sound), the latter depends on the "spring constants" in the pressure terms and how many neighbors particles have, as well as the mass of the particles. Both higher speed of sound and lower fluid velocities make a more incompressible fluid, and both both cause "stiffness" in the sense that they increase how many steps are needed to for the dynamics you are interested in happen.


Not the answer you're looking for? Browse other questions tagged or ask your own question.