[952]Negative concentrations

[김태언]

There are several possible reasons for negative concentrations:

Numerical errors

The most common reason for negative concentrations is numerical noise: when the species concentration approaches zero, the numerical noise becomes significant in comparison to the concentration. If you are seeing negative concentrations of very small magnitude, numerical noise is probably the cause. In a pure convection-diffusion problem, the scale of the concentration does not matter, so in order to avoid the problem, you can add an arbitrary baseline concentration to keep the resulting overall concentration above zero. However, with a reaction term which depends on the concentration, the scale and origin do matter, which means that you need to think of other ways to keep the concentration strictly positive.

Discontinuous concentrations

Another common cause for slightly negative concentrations is a discontinuity in space or time, for instance in the initial condition. As an example, consider the one dimensional time dependent Convection and Diffusion equation dc/dt-0.01*div(grad(c)+u*c)=0, where the convection is in the positive x-direction (i.e. the direction vector u=1), with a uniformly zero initial condition, and boundary conditions which set the concentration at the end nodes to one and zero respectively. The physical interpretation of this PDE is an initially sharp, gradually diffusing front moving in the positive x-direction. However, for the default shape function (second order Lagrange), only continuous functions are admissible as FEM solutions, for which reason the discontinuous initial value is modified before the time-iterations can begin. This often results in a small dip in the solution for t=0, and in the above example, the concentration will locally be slightly negative at t=0, as shown in the figure below.

Solutions to the time-dependent Convection and diffusion equation at times t=0, 0.01, 0.1, 0.2, 0.3. This type of behaviour can also result in wildly oscillating solutions and convergence problems.

This problem can be avoided by using one of COMSOL Multiphysics's built in smoothed Heaviside functions, for instance flc2hs, in order to smooth out the initial discontinuity in a controlled manner. For example, in the Convection and Diffusion problem described above, you could, instead of the uniformly zero initial condition, use the initial condition (1-flc2hs(x-0.01,0.05)) in order to avoid negative values for concentration at t=0. For more information please refer to 'Variables and Expressions' in the COMSOL Multiphysics User's Guide.

Incorrect reaction term

Usually a significantly negative concentration (i.e. not noise around zero) indicates that the underlying mathematical model does not correctly describe the physics. In this case "fixing" the numerics does not take care of the problem. One potential cause is that you have a constant sink in your reaction term, which is an approximation that only works for large concentrations. When the concentration reaches zero, the reaction term continues to consume the species, finally resulting in negative concentration. In order to avoid this problem you need to make sure that your reaction rate is such that when the concentration of the species approaches zero, then so does the rate of its disappearance. This can be achieved for instance by using the following modifications for the source term: Q*(c>0), or Q*flc1hs(c-w/2,w), where w is a sufficiently small number.

Mesh resolution

One other thing that can be indicated by significantly negative concentrations is the lack of mesh resolution. The resulting convergence problems are often the underlying issue when negative concentrations are observed in high convection regimes (high Peclet number) and in those with large reaction terms or fast kinetics (high Damkohler number). In knowledgebase entry 103 a few tricks are presented which may alleviate these types of problems without extensively refining the mesh. It can also be useful to investigate whether the negative concentration problem gets better or worse with mesh refinement. If better, then you know in which direction to go. If worse, then the physics of the model probably need to be checked.

Logarithmic concentrations

A nice way of eliminating mesh resolution problems and negtive dips is to use the logarithm of the concentration and not the concentration itself as the dependent variable. The reason for this is that a linearly varying mesh can sometimes not capture the exponential behaviour of the concentration changes. In addition, modeling the logarithm of the concentration ensures that the real concentration can never become negative during the solution process. The two files below give an example of a stiff reacting system modeled with log(c).

Related Files

	catalytic_reaction_ignition.mph	339 KB
	catalytic_reaction_ignition.pdf	493 KB