Note: This discussion is about an older version of the COMSOL Multiphysics® software. The information provided may be out of date.
Discussion Closed This discussion was created more than 6 months ago and has been closed. To start a new discussion with a link back to this one, click here.
Numerical instability when (conditionally) switching a PDE on
Posted Apr 3, 2012, 5:08 a.m. EDT Materials, Modeling Tools & Definitions, Parameters, Variables, & Functions, Studies & Solvers Version 4.2a 6 Replies
Please login with a confirmed email address before reporting spam
Dear colleagues,
I have a challenge for you if you would be so kind :) .
I am trying to model char combustion, for which I use a PDE function in the form like this:
dm_char/dt = -P_O2*k_eff*da_comb, in which P_O2 is the partial oxygen pressure, k_eff a rate function dependent on the amount of char left, a kinetic expression and a diffusive expression and da_comb is a switch defined as:
if(m_VOL_l_kg<0.02*Volatiles_init_wtkg_zone,1,0)
In every element there is water, volatiles and char and I want the char combustion only to start when the amount of volatiles is less then 2% left. So I keep the right hand side of the PDE zero until the condition is met. So when m_VOL<2%, the PDE activates and is dependent on the mentioned variables. This seems to work. The problem is that I want to model the energy that is released by this reaction and that does not work. It doesn't converge.
I think the problem is an instability because of the switching on of the PDE. Please view the attached pictures.
In the pictures you will see the char conversion flux (kg/s) vs time for certain x positions in the reactor. (It is a porous bed of wood particles 1 m long through which hot air flows). The red line is the da_comb switch (scaled). The instability (green) can be clearly seen. I also included an overview picture of the PDE, maybe it is useful.
The gradient is quite steep, but more elements and smaller timesteps have not worked. Any tips, pointers, hints?
I have not included the model since it is quite large and it takes a lot of time to calculate.
I find it interesting in the graph that the instability starts BEFORE the switch is turned on. Why is that?
I would be very grateful for any help you can give on this.
Thanks and kind regards,
Ray.
PS.
I tried flc1hs, which is supposed to be some kind of damped heavyside function, but that didn't work.
I am now trying to set initial flux to 0 instead of a small starting number (1e-6), didn't work
I am going to try changing da_comb to if(m_VOL_l_kg>0.02*Volatiles_init_wtkg_zone,0,1) (reversed)
I am going to try not using da_comb at all and see what happens.
I have a challenge for you if you would be so kind :) .
I am trying to model char combustion, for which I use a PDE function in the form like this:
dm_char/dt = -P_O2*k_eff*da_comb, in which P_O2 is the partial oxygen pressure, k_eff a rate function dependent on the amount of char left, a kinetic expression and a diffusive expression and da_comb is a switch defined as:
if(m_VOL_l_kg<0.02*Volatiles_init_wtkg_zone,1,0)
In every element there is water, volatiles and char and I want the char combustion only to start when the amount of volatiles is less then 2% left. So I keep the right hand side of the PDE zero until the condition is met. So when m_VOL<2%, the PDE activates and is dependent on the mentioned variables. This seems to work. The problem is that I want to model the energy that is released by this reaction and that does not work. It doesn't converge.
I think the problem is an instability because of the switching on of the PDE. Please view the attached pictures.
In the pictures you will see the char conversion flux (kg/s) vs time for certain x positions in the reactor. (It is a porous bed of wood particles 1 m long through which hot air flows). The red line is the da_comb switch (scaled). The instability (green) can be clearly seen. I also included an overview picture of the PDE, maybe it is useful.
The gradient is quite steep, but more elements and smaller timesteps have not worked. Any tips, pointers, hints?
I have not included the model since it is quite large and it takes a lot of time to calculate.
I find it interesting in the graph that the instability starts BEFORE the switch is turned on. Why is that?
I would be very grateful for any help you can give on this.
Thanks and kind regards,
Ray.
PS.
I tried flc1hs, which is supposed to be some kind of damped heavyside function, but that didn't work.
I am now trying to set initial flux to 0 instead of a small starting number (1e-6), didn't work
I am going to try changing da_comb to if(m_VOL_l_kg>0.02*Volatiles_init_wtkg_zone,0,1) (reversed)
I am going to try not using da_comb at all and see what happens.
Attachments:
6 Replies Last Post Apr 16, 2012, 12:41 p.m. EDT
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
I included 2 more pictures. The first shows the error that I get when include the energy of combustion. I don't know why it says that it has found an undefined value in the stiffness matrix, but it just won't converge.
The second picture shows the amount of char left in a certain element. The amount decreases nicely as one would expect, but at the top small instabilities can be seen. These have the same origin as the instabilities in the flux plot, but where they come from...?
Please let me know if any of you have any suggestions or ideas.
Many thanks and kind regards,
Ray
PS.
I am currently simulating with a different initial condition and reversely defined da_comb and also one with no PDE switch used at all (so it will start from the beginning, which is physically wrong, but numerically interesting).
PPS.
Both methods did not converge.
The second picture shows the amount of char left in a certain element. The amount decreases nicely as one would expect, but at the top small instabilities can be seen. These have the same origin as the instabilities in the flux plot, but where they come from...?
Please let me know if any of you have any suggestions or ideas.
Many thanks and kind regards,
Ray
PS.
I am currently simulating with a different initial condition and reversely defined da_comb and also one with no PDE switch used at all (so it will start from the beginning, which is physically wrong, but numerically interesting).
PPS.
Both methods did not converge.
Attachments:
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi Raymond
Your work looks quite interesting. I have combustion model where I model gas combustion on grate. Grate part is modeled only as inlet "surface" which is providing gas for combustion. Due to this combustion has no effect to grate drying/devolatilization. Where did you found data for these PDE:s? I would like to make similar 1D model for grate part.
To your problem. Have you tried to different type of solver. I have used segregated solver where flow, temperature (energy) and concentration is calculated in own segregated steps. This helps sometimes...
Best regards
Tero
Your work looks quite interesting. I have combustion model where I model gas combustion on grate. Grate part is modeled only as inlet "surface" which is providing gas for combustion. Due to this combustion has no effect to grate drying/devolatilization. Where did you found data for these PDE:s? I would like to make similar 1D model for grate part.
To your problem. Have you tried to different type of solver. I have used segregated solver where flow, temperature (energy) and concentration is calculated in own segregated steps. This helps sometimes...
Best regards
Tero
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi Tero,
thanks for your interest. What PDE's are you interested in exactly? I use Arrhenius type equations (first order in amount of volatiles, water left) and got the data from various papers. There are a lot of papers on these subjects.
For drying I am using functions from a paper called: ''Measurements and particle resolved modelling of heat-up and drying of a packed bed.''
For devolatilization there are a lot of papers by Colomba Di Blasi and others.
For combustion of the solid char I am using relations used in ''Computational fluid dynamic simulation of a solid biomass combustor modelling approaches.''.
But these last functions, to describe the char combustion, are giving me a really hard time. I can't get them to work in COMSOL.
I haven't modeled the gas reactions yet, but if I do it, it will be done using a chemical equilibrium model (using matlab). I already modeled the devolatilization this way.
To your problem. Have you tried to different type of solver. I have used segregated solver where flow, temperature (energy) and concentration is calculated in own segregated steps. This helps sometimes...
Yeah, I have tried all timestepping methods and manny solvers, but to no avail. I tried using a segregated solving method, but it didn't work. I am not entirely sure how to set the options for it... Will it still give an accurate result if it works? Are you experienced with it?
Kind regards,
Ray
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi Ray
Thank you for these. There were familiar papers which I already had but also new ones.
So far I have tried to avoid making own PDE:s and used only ready modules of Comsol (not very handy with PDE:s). But it looks like that I have learn how to make and use them...
My model is stationary where I have to advance by steps to get it working (first flow, then reactions, last energy). Comsol seems to need very good start values or else solving doesn't work. For example I have Arrhennius type reaction rate for CO which is highly temperature dependent. To get this working I must have solution from mixing base reaction. I also have to remove temperature parameter from one segregated step to another step (with concentration parameters). I think this doesn't help you much since your model is time dependent.
Sorry, I don't have experience what is accuracy difference between segregated or coupled solver.
Best regards
Tero
Thank you for these. There were familiar papers which I already had but also new ones.
So far I have tried to avoid making own PDE:s and used only ready modules of Comsol (not very handy with PDE:s). But it looks like that I have learn how to make and use them...
My model is stationary where I have to advance by steps to get it working (first flow, then reactions, last energy). Comsol seems to need very good start values or else solving doesn't work. For example I have Arrhennius type reaction rate for CO which is highly temperature dependent. To get this working I must have solution from mixing base reaction. I also have to remove temperature parameter from one segregated step to another step (with concentration parameters). I think this doesn't help you much since your model is time dependent.
Sorry, I don't have experience what is accuracy difference between segregated or coupled solver.
Best regards
Tero
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
I am now sure that the oscillations happen because of the actual switching. Before I was thinking it could have been the combustion expression, but now I am sure that is not the problem.
To recapitulate: I have a PDE (general form) and I want to activate it (start using its expression) when the mass that is left (calculated by another PDE) falls below a certain value. This gives nasty numerical oscillations when the switching happens.
I tried to solve this by:
1) smaller timesteps (0.1 s, instead of 10 or 5)
2) more elements (3000 instead of 400, only makes oscillations smaller, does not remove them)
3) initial flux=0 does not work, initial flux = -1e-6 or -1e-8 solves but does not remove problems, or initial flux is taken the same as PDE function but with initial mass or that but with da_comb included, but these do not solve at all.
4) not use da_comb at all, but this does not work because the initial conditions of the simulation are not suitable for the PDE function, it has to happen later when the conditions are ok.
5) free BDF timestepping, intermediate, strict
6) quintic PDE shape function, hermite, discontinuous lagrange and setting all model functions to cubic.
7) make mu_gas, rho_gas, Cp_gas and k_gas constant as well as Cp_solid, k_solid
8) I tried making as many terms in the PDE constant, like this:
da_comb*-4800[Pa]*(1/((1/(1e-7*m_char_zone))+((1/0.08)*(400)*(R/(12*A_S_zone)))))
But this only seems to make it worse.
9) turned consistent stabilization on and off under the heat transfer modules
10) changed RelTol to 1e-6 and 1e-10
11) changes AbsTol to 1e-3 en 1e-6
12) used a smoothed step function to replace the switch, but this just makes the oscillations nice and round, but does not remove them.
13) Tried to make function dependent on temperature (instead of mass left) but this did not work either.
PLEASE, any suggestions? At all?
Thanks for the trouble
Ray.
PS. I included the model showing this behaviour (solves really quick).
The function of interest (da_comb) is under 'GlobalZoneVariables' at the bottom. The first graph shows the oscillation.
To recapitulate: I have a PDE (general form) and I want to activate it (start using its expression) when the mass that is left (calculated by another PDE) falls below a certain value. This gives nasty numerical oscillations when the switching happens.
I tried to solve this by:
1) smaller timesteps (0.1 s, instead of 10 or 5)
2) more elements (3000 instead of 400, only makes oscillations smaller, does not remove them)
3) initial flux=0 does not work, initial flux = -1e-6 or -1e-8 solves but does not remove problems, or initial flux is taken the same as PDE function but with initial mass or that but with da_comb included, but these do not solve at all.
4) not use da_comb at all, but this does not work because the initial conditions of the simulation are not suitable for the PDE function, it has to happen later when the conditions are ok.
5) free BDF timestepping, intermediate, strict
6) quintic PDE shape function, hermite, discontinuous lagrange and setting all model functions to cubic.
7) make mu_gas, rho_gas, Cp_gas and k_gas constant as well as Cp_solid, k_solid
8) I tried making as many terms in the PDE constant, like this:
da_comb*-4800[Pa]*(1/((1/(1e-7*m_char_zone))+((1/0.08)*(400)*(R/(12*A_S_zone)))))
But this only seems to make it worse.
9) turned consistent stabilization on and off under the heat transfer modules
10) changed RelTol to 1e-6 and 1e-10
11) changes AbsTol to 1e-3 en 1e-6
12) used a smoothed step function to replace the switch, but this just makes the oscillations nice and round, but does not remove them.
13) Tried to make function dependent on temperature (instead of mass left) but this did not work either.
PLEASE, any suggestions? At all?
Thanks for the trouble
Ray.
PS. I included the model showing this behaviour (solves really quick).
The function of interest (da_comb) is under 'GlobalZoneVariables' at the bottom. The first graph shows the oscillation.
Attachments:
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi Ray
Sorry, I cannot help you. I don't have license for Chemical reaction module.
Best regards
Tero
Sorry, I cannot help you. I don't have license for Chemical reaction module.
Best regards
Tero