# Discussion Forum

## How to solve the problem of negative concentration

 Topics: Fluid Flow, 3.5a
RSS feed   |   Turn on email notifications   |   32 Replies   Last post: June 20, 2015 3:05pm UTC

Yuanhai Su

June 6, 2012 10:27am UTC

How to solve the problem of negative concentration

I am a totally new user on comsol. Now I try to simulate the mixing in the microchannel without reactions. I use the laminar flow and Transport of dilluted species. There is a similar example from comsol model library. However, I got negative concentration and the larger value beyond physic in simulation results. Somebodies say this problem can be solved by seting up logarithmic transform of concentration. But how to do that in the model of Transport of diluted specices? Please help me. Thank you so much.

Ivar Kjelberg

June 6, 2012 8:05pm UTC in response to Yuanhai Su

Re: How to solve the problem of negative concentration

Hi

I do not know enough about your model, but I have noticed that if this happens in a time dependent solver case, it often comes from an initial condition of a concentration = 0 on the domain and a value of C>>0 at a boundary. Then the mesh along this boundary is often rather coarse and the diffusivity low, such that the time stepping times mesh size does not link correctly to the diffusivity and the transient gradient along the C>0 boundary is not well resolved. This then results in concentration oscillations and mostly negative values.

If so, you can avoid this by ramping up your concentration with a smoothened step function (Definition Functions) starting at or closer to "0" at t=0 and to use a denser mesh along these abrupt gradient regions

But your issue might be another case too, this is just one guess and one issue I have often observed

--
Good luck
Ivar

Yuanhai Su

June 7, 2012 11:19am UTC in response to Ivar Kjelberg

Re: How to solve the problem of negative concentration

Ivar, thank you so much! My model is stationary, and it seems simple. I agree with your opinion. Maybe it is the problem of mesh. If the mesh is not enough fine, the numerical error is easy to occur (it is right or not?). Next week I will change another new computer with much larger memory. Another thing let me feel inexplicable. I use the segregated strategy for solving velocity and concentration fields under the same meshes. The velocity field had accurate results. But why could be the concentration field right? Certainly, the PDE for laminar flow is relatively simple, and even it has analyse solution. However, the resolution of concentration field must base on the velocity field. Is the the reason that the velocity field is more easy to get the right simulation results?

Ivar Kjelberg

June 7, 2012 11:46am UTC in response to Yuanhai Su

Re: How to solve the problem of negative concentration

Hi
Indeed my reark is more applicable to time series analyses
Mesh is important and a too coarse mesh is a common issue, specially with little RAM.

Normally a laminar flow needs quite a dense mesh, but the concentration case could even need a finer mesh, it depends on the diffusivity

Using segegated solvers is one way to minimise RAM, but if you have a werey weak coupling you could consider even solving on.y the CFD part first and then separately the concetration by using the velocity field from the CFD case.

But this might not be valid for your model, it all depends, "many ways to Rome"

--
Good luck
Ivar

Suresh Seetharam

June 8, 2012 9:11am UTC in response to Yuanhai Su

Re: How to solve the problem of negative concentration

Hi,

For stationary you have (ignore my tensorial mistakes),

where c is the primary variable defined in comsol.

In comsol model builder, define a variable

c1=log(c)

Therefore,

c = exp(c1) and

Now go to "transport in dilute species" and click on convection and diffusion equation.

In the velocity field, u, feed
velocity in x = exp(c1)*ux
velocity in y = exp(c1)*uy

ux and uy are the velocities you get from your flow solution.

Now, go to the diffusion field, and feed:

D*exp(c1)

where D is the effective diff coeff.

Finally, if you are solving with initial conc = 0, then put a very small number say 10^-20 so that the variable c1 does not have a problem (log(0)).

This is roughly the derivation, but I have never tried this. If there is a mistake let me know.

You could switch on adaptive meshing with a reasonable initial meshing as a starting point.

Suresh

Yuanhai Su

June 8, 2012 10:13am UTC in response to Ivar Kjelberg

Re: How to solve the problem of negative concentration

Thank you so much! The CFD part was first simulated, and then the concentration field was calculated based on this velocity field. I used the segregated strategy. Now, I am out of my work for vacation. After returning, I will try to simulate it with more finer meshes.

Yuanhai Su

June 8, 2012 10:14am UTC in response to Suresh Seetharam

Re: How to solve the problem of negative concentration

Suresh, Thank you! I will try it like that after returning to work.

John H.

June 8, 2012 2:08pm UTC in response to Suresh Seetharam

Re: How to solve the problem of negative concentration

In comsol model builder, define a variable

c1=log(c)

[...]

Now, go to the diffusion field, and feed:

D*exp(c1)

where D is the effective diff coeff.

[...]

This is roughly the derivation, but I have never tried this. [...]

I've tried this once — some time ago when I had a numerically unstable transport problem with negative concentrations popping up all over the place. I derived the formulas the same way you did (and for the record, see nothing wrong in your derivation), All of this because of the promising advice given in that Comsol Knowledge Base article, which reasons that, if you take the logarithm, the concentrations simply cannot become negative anymore.

Well, this is obviously true. However, at least in my case, the logarithmic transformation made everything much worse. Sure, the concentrations were no longer negative. But they would fluctuate much more wildly than before.

In the end, this isn't surprising at all. As one can see above, transforming the equation in terms of c to something in terms of log(c) introduces those pesky exponentials in diffusion coefficients, and, usually, in source terms as well (since they typically include a dependence on concentration). So except for some pathological cases, this doesn't look like it would solve the stability problem at all.

Also, I couldn't find much literature on "logarithmic transformation" or something (but maybe I didn't search for the right keywords), but the few things I did find here and there pretty much confirmed my observation.

Yuanhai Su

June 8, 2012 3:03pm UTC in response to John H.

Re: How to solve the problem of negative concentration

Thank you. You means that it is not a good idea to solve this problem by the method of logarithmic transformation? Except for choosing finer meshes and the method of lograrithmic transformation, it seems that there is not anything other strategy to solve this kind of problem.

In comsol model builder, define a variable

c1=log(c)

[...]

Now, go to the diffusion field, and feed:

D*exp(c1)

where D is the effective diff coeff.

[...]

This is roughly the derivation, but I have never tried this. [...]

I've tried this once — some time ago when I had a numerically unstable transport problem with negative concentrations popping up all over the place. I derived the formulas the same way you did (and for the record, see nothing wrong in your derivation), All of this because of the promising advice given in that Comsol Knowledge Base article, which reasons that, if you take the logarithm, the concentrations simply cannot become negative anymore.

Well, this is obviously true. However, at least in my case, the logarithmic transformation made everything much worse. Sure, the concentrations were no longer negative. But they would fluctuate much more wildly than before.

In the end, this isn't surprising at all. As one can see above, transforming the equation in terms of c to something in terms of log(c) introduces those pesky exponentials in diffusion coefficients, and, usually, in source terms as well (since they typically include a dependence on concentration). So except for some pathological cases, this doesn't look like it would solve the stability problem at all.

Also, I couldn't find much literature on "logarithmic transformation" or something (but maybe I didn't search for the right keywords), but the few things I did find here and there pretty much confirmed my observation.

John H.

June 8, 2012 3:43pm UTC in response to Yuanhai Su

Re: How to solve the problem of negative concentration

You means that it is not a good idea to solve this problem by the method of logarithmic transformation? Except for choosing finer meshes and the method of lograrithmic transformation, it seems that there is not anything other strategy to solve this kind of problem.

Yes, I think the logarithmic transformation is unlikely to help. So use it only as a last, desperate resort. The meshing is much more important, probably the most important setting. Other settings that, in my experience are worth playing around with, are: the tuning parameters for the various stabilization techniques offered by the Transport of Diluted Species interface, the (relative/absolute) tolerances, and (manual) scaling of dependent variables — roughly in this order of importance.

Yuanhai Su

June 8, 2012 8:00pm UTC in response to John H.

Re: How to solve the problem of negative concentration

Thank you. Various stabilization technique? Which parameters are related to this technique? I know that if the the (relative/absolute) tolerance is larger, it is more easy to get convergence. The default value of tolerance is 0.001. Sometimes I can not get the convergence with the default value, but after changing it to 0.01, I got the convergence results. How does the tolerance actually affect the simulation results?

Yuanhai Su

June 15, 2012 2:15pm UTC in response to Suresh Seetharam

Re: How to solve the problem of negative concentration

Hello.
I return to work. i had tried your suggestion. But it still did not pass the problem of negative concentration. It showed that 'Attempt to evaluate real logarithm of negative number.
- Function: log
'

How can I use comsol physcis well? So far I do not get the right results. It seems that they is very few parameters that can be controlled. Just mesh and solvers can be controlled............ Maybe I am wrong.....

I also attach my simulation as following:

Hi,

For stationary you have (ignore my tensorial mistakes),

where c is the primary variable defined in comsol.

In comsol model builder, define a variable

c1=log(c)

Therefore,

c = exp(c1) and

Now go to "transport in dilute species" and click on convection and diffusion equation.

In the velocity field, u, feed
velocity in x = exp(c1)*ux
velocity in y = exp(c1)*uy

ux and uy are the velocities you get from your flow solution.

Now, go to the diffusion field, and feed:

D*exp(c1)

where D is the effective diff coeff.

Finally, if you are solving with initial conc = 0, then put a very small number say 10^-20 so that the variable c1 does not have a problem (log(0)).

This is roughly the derivation, but I have never tried this. If there is a mistake let me know.

You could switch on adaptive meshing with a reasonable initial meshing as a starting point.

Suresh

Tero Hietanen

June 20, 2012 8:39am UTC in response to Yuanhai Su

Re: How to solve the problem of negative concentration

Hi,

Here is one trick to solve negative concentrations.

I made variable "C" which is defined by "if(c<0,eps,c)" under Definitions => Variables. Instead of concentration "c" I used variable "C". Now variable "C" is never below zero. Eps is very small number defined is Comsol.

Best regard

Tero

Yuanhai Su

June 23, 2012 9:56pm UTC in response to Tero Hietanen

Re: How to solve the problem of negative concentration

Thank you! I just returned to work yesterday from a vacation. I will try to do it like that and tell you the results.

I try to use the boundary layer mesh. The appropriate velocity field could be obtained, but the concentration field was not right (some of the values much exceeded the expectable largest value, and some were much negative). Now my computer memory is enough. However, I do not get the proper results so far.........

Yuanhai Su

June 26, 2012 12:15pm UTC in response to Tero Hietanen

Re: How to solve the problem of negative concentration

Thank you. Your method seems available. But if another condition that the concentration value exceeds the expectable largest value also occurs, how to set the variable C? That is, if c<0, c=eps; if c>1, c=1. Please help me. Thank you again.

Hi,

Here is one trick to solve negative concentrations.

I made variable "C" which is defined by "if(c<0,eps,c)" under Definitions => Variables. Instead of concentration "c" I used variable "C". Now variable "C" is never below zero. Eps is very small number defined is Comsol.

Best regard

Tero

Tero Hietanen

June 26, 2012 12:35pm UTC in response to Yuanhai Su

Re: How to solve the problem of negative concentration

Hi,

Try this: if(c<0,eps,if(c>1,1,c)).

Best regards

Tero

Yuanhai Su

June 26, 2012 1:10pm UTC in response to Tero Hietanen

Re: How to solve the problem of negative concentration

Dear Tero,

It is so good. You are good on using comsol physics. I should learn from experienced experts on this field like you! Thank you so much!

Hi,

Try this: if(c<0,eps,if(c>1,1,c)).

Best regards

Tero

Yuanhai Su

June 27, 2012 7:39pm UTC in response to Tero Hietanen

Re: How to solve the problem of negative concentration

Hello Tero,

I think about this issue again. I use the method that you mentioned, all positive values of concentration could be obtained. However, maybe it just related to the range of the values that you want to get. For example, by this method we can use C to instead of c, but if you want to show the values of C after finishing the simulation, it is the same as the previous ones without instead of c by C.
Now I doubt this method.........We can continue to discuss this issue. Thank you!

Hi,

Try this: if(c<0,eps,if(c>1,1,c)).

Best regards

Tero

Tero Hietanen

June 28, 2012 5:32am UTC in response to Yuanhai Su

Re: How to solve the problem of negative concentration

Hi,

I not sure did I understand you correctly...? This is just modeling trick how to get over problems during solving. Like Ivar said other things (mesh etc) are still main solution to this problems (if you want to get accurate results). When other things are ok, then difference between values C and c should be minimal.

Best regards

Tero

Yuanhai Su

June 28, 2012 7:49am UTC in response to Tero Hietanen

Re: How to solve the problem of negative concentration

Hello, good morning. For my model, the concentration should be between 0-1 mol/l. However, some values are over 1, the largest values are about 1.03 (3%of the expectable largest value), and some values are under 0, about -0.03 for the lowest values. The two inlet zones (one inlet concentration is set as 1, another inlet concentration is set as 0) are special, which produce most of these special values. Even I let the meshes finer (after enough meshes), these special values almost do not change. Certainly, when I use coarse meshes, the non-physical values are more obvious. Maybe it is different to totally eliminate these special values, and they are related to calculation error.

Hi,

I not sure did I understand you correctly...? This is just modeling trick how to get over problems during solving. Like Ivar said other things (mesh etc) are still main solution to this problems (if you want to get accurate results). When other things are ok, then difference between values C and c should be minimal.

Best regards

Tero

Ivar Kjelberg

June 28, 2012 8:15am UTC in response to Yuanhai Su

Re: How to solve the problem of negative concentration

Hi

often in time series models you have steep gradients at t=0 i.e. concentration at a boundary, nothing inside the domain.
Then I have noticed that using a "boundary mesh" on the c=1 (or T=cnst) boundary can help, as you do not need always the very fine mesh in the full domain (but this is very model dependent) this works fine for me with most of my diffusion driven models, such as HT and Chemistry models

--
Good luck
Ivar

Yuanhai Su

June 28, 2012 10:37am UTC in response to Ivar Kjelberg

Re: How to solve the problem of negative concentration

Hello Ivar,
Thank you so much.

My model is stationary model with laminar flow and transport of diluted specie models. To my model I used the boundary layer mesh for the laminar flow, and got physical solution (I had tested the dependence of simulation result on the mesh number, and found it did not change after putting enough meshes and it seemed right from physical meaning). However, when I tried to use the boundary layer mesh for the transport, it got much large negative concentration (for example, some points reached -0.07 mol/l, 7%of the expectable largest value; certainly, I also tested the effect of mesh number, these results were also not feasible). Therefore, I used boundary layer meshes for laminar, and used meshes without boundary layer mesh for transport. These two parts were separated to solve, and the transport part was based on the velocity field obtained from the laminar flow solution. With the later method, I got the concentration field. However, the concentration in the two inlet parts were not so good from physical opinion (larger than 1, and smaller than 0, -3%-3% and most of them were between -2%-2% of the expectable largest value; certainly, they were not far from the expectable values, and the main geometry part (main channnel) could get proper results).

Hi

often in time series models you have steep gradients at t=0 i.e. concentration at a boundary, nothing inside the domain.
Then I have noticed that using a "boundary mesh" on the c=1 (or T=cnst) boundary can help, as you do not need always the very fine mesh in the full domain (but this is very model dependent) this works fine for me with most of my diffusion driven models, such as HT and Chemistry models

--
Good luck
Ivar

Ivar Kjelberg

June 28, 2012 10:41am UTC in response to Yuanhai Su

Re: How to solve the problem of negative concentration

Hi

indeed my remark was more for time dependent cases with high gradients due to initialcondition steps.

The next on the list, in my view, is the numerical stabilisation parameters, but I'm not use to tweak these for such models so I have no real idea where to start. but its worth to read carefully the "stabilization" chapter in the COMSOl pdf documents

--
Good luck
Ivar

Thanh Lem

December 12, 2012 1:20am UTC in response to Tero Hietanen

Re: How to solve the problem of negative concentration

Hi,

Try this: if(c<0,eps,if(c>1,1,c)).

Best regards

Tero

Thank you so much, Tero. I appreciate your helps.

Naveen Aerpula

December 9, 2013 10:34am UTC in response to Tero Hietanen

Re: How to solve the problem of negative concentration

if(c<0,eps,if(c>1,1,c))

Dear Tero,
First i would like to say thank you because ,i am doing my MS project "modeling of electrosorption" using comsol version 4.3 at Hindustan Unilever R&D,Bangalore and i got the same problem of getting negative values or larger values which are not expected.I searched for resolving this problem in various literatures and available comsol materials but i did not find anywhere and finally i found randomly here(discussion forum) and i tried it and its working well.But i have a few doubts about defining this condition.
(1)How much small number is this eps? or is it constant small number?
(2)what is the condition if c lies between 0 and 1 i.e (0=<c=<1) (here =< represents less than or equall)?

could you please give me the combined condition for this(for all three conditions).

Thanking you,

A.Naveen

Tero Hietanen

December 9, 2013 7:31pm UTC in response to Naveen Aerpula

Re: How to solve the problem of negative concentration

Hi,

Eps is coming from Comsol and it is about 2.2e-16. I have used it instead of zero.

If c is between 0 and 1 it is keeping original value of c. For example if:

c is -1 => return value is 2.2e-16 (eps),
c is 2 => return value is 1,
c is 0.5 => return valus is 0.5,
etc...

Best regards

Tero

Alex Alexian

April 10, 2014 6:27am UTC in response to Tero Hietanen

Re: How to solve the problem of negative concentration

Hi,

Here is one trick to solve negative concentrations.

I made variable "C" which is defined by "if(c<0,eps,c)" under Definitions => Variables. Instead of concentration "c" I used variable "C". Now variable "C" is never below zero. Eps is very small number defined is Comsol.

Best regard

Tero

Hi Tero
I try this for species transport in porous media but my finally concentrations are negative. my problem is unsteady.
Best regards

Tero Hietanen

April 11, 2014 10:30am UTC in response to Alex Alexian

Re: How to solve the problem of negative concentration

Hi,

Diffucult to say without seeing your model... have you check your model mesh? I have never tested this method with unsteady case.

Best regards

Tero

Sanjoy Nandi

June 12, 2014 11:11am UTC in response to Suresh Seetharam

Re: How to solve the problem of negative concentration

Hi Suresh,

Thanks for the analysis. Is there any way to solve with logarithmic concentration(or any alternative approach) for time dependent drift-diffusion equation. Physically why we can't use your approach time dependent eqn.?

Thanks again,
Sanjoy

Peng Lin

November 12, 2014 5:39pm UTC in response to Ivar Kjelberg

Re: How to solve the problem of negative concentration

Dear All,

I met the exactly same problem as you described.

I am working on 3 coupled PDEs, two from the Joule heating, one PDE form diffusion-drift equation.

The domain is a rectangle.

The boundary condition for the Diffusion-Drift is Zero flux on all edges.

The initial condition for the diffusion-drift is 1E27 [m^-3] in the left half domain and 0 for the right half.

The solver is time dependent with constant damping factor Newton's method and a MUMPS linear solver.

The issue:

Negative concentration keeps to present.

The approaches I have tried:

(1) C=if(c<0,eps,c). No problems caused by negative ion concentration. However, there is a time step within which the maximum iterations of (constant damping) newton's method has been reached. (Resolved by C=if(c<0,1E25,c))

(2) Ultra dense mesh+Boundary mesh layer. Problem "NaN or Inf found" in the linear MUMPS solver. And the speed is so slow that I believe segregated solver has to be implemented.

(3) As the diffusion drift equation is implemented in the PDE interface, there is no stabilization option, although the Joule heating equations are of available stabilization options.

(4) Put a background ion concentration for the right half domain. However, when time increases, there will be negative ion concentration in some regions of the left-half domain.

Can any one help? Many many thanks in advance.

Peng Lin

November 14, 2014 1:22am UTC in response to Peng Lin

Re: How to solve the problem of negative concentration

Up Up Up

Shiva Asapu

June 20, 2015 11:46am UTC in response to Tero Hietanen

Re: How to solve the problem of negative concentration

Hi,

Here is one trick to solve negative concentrations.

I made variable "C" which is defined by "if(c<0,eps,c)" under Definitions => Variables. Instead of concentration "c" I used variable "C". Now variable "C" is never below zero. Eps is very small number defined is Comsol.

Best regard

Tero

Hi,

What do you mean by - [Instead of concentration "c" I used variable "C"]?
Did you replace "c" with "C"? But c is the main variable.
I'm sorry, my doubt might be trivial since I'm new to COMSOL.

Thank You,
Shiva

Shiva Asapu

June 20, 2015 3:05pm UTC in response to Peng Lin

Re: How to solve the problem of negative concentration

(1) C=if(c<0,eps,c). No problems caused by negative ion concentration. However, there is a time step within which the maximum iterations of (constant damping) newton's method has been reached. (Resolved by C=if(c<0,1E25,c))

Hi,

I believe that I'm trying to solve the same set of PDEs that you are trying to solve as well. Can you please let me know how did you use 'C'?
It is clear that 'c' is the variable that we declare for the number of species in the drift-diffusion PDE. But after you define 'C' as 'if(c=0,eps,c)' or 'if(c<0,1E25,c)', how do we use it? I mean, are we supposed to replace 'c' with 'C' anywhere?

I have simply defined C=if(c=0,eps,c) as a local variable for all the domains where the species exist but still can find negative values of 'c' as and when the the solver solves.

And also which solver would you recommend? I tried both MUMPS and PARADISO. Unfortunately, each solver gives different results.

I'd really be grateful if you can help.

Sincerely,
Shiva

Rules and guidelines