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.
Crazy oscillations in Darcy velocity profile
Posted 29 févr. 2012, 05:29 UTC−5 Fluid & Heat, Results & Visualization, Studies & Solvers Version 4.2a 17 Replies
Please login with a confirmed email address before reporting spam
Hello dear colleagues,
A small question this time: Does anyone know how to remove the crazy (!) oscillations that occur in my velocity profiles (and all other graphs dependent on the velocity) ?
I am using 1D Darcy, Transport of species in porous media and 2 heat transfer modules. Changing solvers, timesteps, boundary conditions and number and type of elements has no effect. If I take the gas viscosity (mu) constant instead of using a (linear, temperature dependent) function, it is less bad, but still there.
So, any ideas? (please view attached picture to see what I mean ;) )
Thank you so much for your time.
Kind regards,
Ray
Attachments:
Please login with a confirmed email address before reporting spam
Without seeing your model I can't say much about it, I also don't have experience with the kind of simulations you're doing.
Maybe you need to refine your mesh?
Check the convergence plots as well, see if you see something weird.
I would suggest to do a very basic simple simulation and see if you can reproduce the artifacts.
-Jaap
P.s, nice name.
Please login with a confirmed email address before reporting spam
What I tried is to get finer mesh, and see if oscillations have decreasing magnitude.
Please login with a confirmed email address before reporting spam
Why not trust a result because there is oscillations ? I would rather propose to try still to understand why, so you get around quicker next time.
Oscillations can come from several issues, such as (but you seem to have ruled out these) mesh density and time stepping, but also from numerical issues such as
underflow overflow, hence from the scaling of the variables. See also the KB on the numerical stabilisation (www.comsol.eu/support/knowledgebase/1104/),
yet another source. And sometimes one forget that a transient or time series solver adds in the inertial terms rho*d^2_u_/dt^2, so with non-coherent initial
values and some BC values one can hit the model like a structure with a hammer, it will resonate quite some time, if it even settles once ;)
Finally,what I have observed when fiddleing with simple PDEs, is that you sometimes get a spiky solution with extrema on the mesh nodes, and that follow the mesh
change, these can come from insufficient BCs so you have 2 solution of your equations and COMSOL is oscillating between both. In fact you get two solutions by
looking at the extremas (in linear discretization) for the same work ;) Probably, this is not your case here with the complex physics you have.
I can propose the "Jelly fish quiz" to illustrate the effects of time series and oscillations:
Start a new COMSOL in 2D, with SOLID,, with a time transient solver, default time steps. Then draw a square of default 1m^2, add a user material with E=1[GPa],
nu=0.3, rho = 2000[kg/m^2] i.e some material. You fix the lower Boundary #2, and apply a time increasing load of Fy=-10[MPa]*t[1/s] on the top Boundary #3. As we
have a compression stiffness of k=E*A/L = 1000[N/um] hence an expected final displacement of 10[mm]. Turn on the Results plot while solving and add a Point Probe
plot on one of the top vertex/points to follow the displacement. This solves and confirms our results, as expected.
Now, lets change the material properties to "jelly fish" E=100[Pa], rho = 1 and Boundary load to -1[Pa]*t[1/s], the load to stiffness relation and displacement
should be the same, what else ?
By observing carefully the probe plot (displacement versus time) we see clearly oscillations, this comes from the speed of sound in the material v=sqrt(E/rho) =
10[m/s] that we excite easily with our increasing boundary load. We can though remove this transient behaviour, by going back to the main solid physics
node, and under the tab: Structural Tansient behaviour, set it to Quasi static. Now we have no oscillations as we have removed the rho*d^2_u_/dt^2 term from the
equations (check the Solid Physics Equation tab).
Next step is to set back the inertial term, remove the *t[1/s] dependency on the boundary load and run the same two material simulations (start with the "Jelly
fish"), the oscillatory behaviour is now strongly enhanced.
il
THese simulatons are all "correct" but we observe different behaviours, and some we had probably forgotten about, as they are not that easy observed ouside the
FEM world
I have attached the 4.2a model hereby ;)
--
Good luck
Ivar
Attachments:
Please login with a confirmed email address before reporting spam
I have read the reference guide (again :P ) and I think I found the source of the problems (which I indeed to a certain degree managed to repeat in a simple model). The problem is that when the Peclet number is greater than 1, you are bound to get numerical instability. On page 354 of the reference guide a measure for instability is given:
Pe= (||v|| * h ) / (2 * c)
v is the convective velocity vector, h is the element size, c is the diffusion coefficient.
When Pe is greater than 1, and you have certain boundary conditions, you get instability.
Fortunately, they give the solution in terms of adding numerical diffusion which you can control under the 'stability' tab of many modules.
Unfortunately, such a tab is nonexistent in the Darcy flow module, the one that I am using of course, so,...I still don't really know how to get rid of this annoying instability.
They proposed to use more elements but it (my real simulation) is already slow with 250 elements, so using 1000's of elements is unworkable and I strongly suspect (based on testsimulations) that this will not work anyway.
My question to you all now is: any other ideas on how to solve this? And why isn't there any option for stabilization in Darcy / Transport of species in porous media ?
PS.
Please view the attached pictures and/or the corresponding attached models if you want to see it in actual action.
The first model is a very simplified model which does seem to show the same behaviour.
The second model is a stripped down version (to make it a lot faster to solve) of my actual model, which shows this behaviour very strongly.
Attachments:
Please login with a confirmed email address before reporting spam
Thanks for the feed back, at least we are all warned for next time ;)
--
Have fun COMSOLING
Ivar
Please login with a confirmed email address before reporting spam
no problem, at least now others with the same problem might find this topic and understand what is going wrong ;)
I personally don't have a solution yet however.
[To all:] Are there no further possibilities to tackle this problem? Any ideas at all? Should I take this up with support?
Thanks again everyone and kind regards,
Ray
Please login with a confirmed email address before reporting spam
This discussion enlightened me on stabilization techniques and verifying Pe number, thanks.
Regarding your model, I see that you sticked with the Darcy's module! I think you should read on multicomponent species transport and see why, usually, a problem is solved for N-1 species. If you solve for N species, like you do, the continuity equation in the Darcy's module becomes useless. The sum of the species equation gives the continuity equation, perform this calculation on paper. I persist in thinking that your problem is overdefined. Pressure and velocity should be calculated from the concentrations of the N species.
Please login with a confirmed email address before reporting spam
This type of convection/diffusion instability, usually occurs both for transient and steady state, and it goes away if you sufficiently refine the mesh. In your simple model refining the mesh does not help, and the oscillations nearly disappear after some time so they could be caused by a sudden change in loading in a transient analysis. If you try ramping up the C0,c2 inflow condition in the Species Transport and also the velocity condition in Darcy’s Laws with a fast acting step function, the solution becomes less oscillatory. I don’t know how much of that applies to your other detailed model. One more factor to investigate!
Please keep us posted if you resolve the issue.
Nagi Elabbasi
Veryst Engineering
Please login with a confirmed email address before reporting spam
and there is a common problem to most diffusion eqations in transient, with large gradients (typically initial conditions of the Dirac type. The solution an erfc(t*...) are very non linear at the transition, particularly for t<< and the meshes used often rather coarse, hence the integration over these first layers of the mesh is very wrong and can lead to large roundoff errors, hence "articifal numerical" oscillations
--
Good luck
Ivar
Please login with a confirmed email address before reporting spam
Hi Raymond,
This discussion enlightened me on stabilization techniques and verifying Pe number, thanks.
Regarding your model, I see that you sticked with the Darcy's module! I think you should read on multicomponent species transport and see why, usually, a problem is solved for N-1 species. If you solve for N species, like you do, the continuity equation in the Darcy's module becomes useless. The sum of the species equation gives the continuity equation, perform this calculation on paper. I persist in thinking that your problem is overdefined. Pressure and velocity should be calculated from the concentrations of the N species.
Hi Francois,
Thank you for your interest.
I am not sure if I am following you completely. I read the knowledge base article on this:
www.comsol.com/support/knowledgebase/1100/
, and I saw the explanation about the N-1 species solving method, but are you sure this applies to my situation also?
I have a diluted system of species (H2, CO, H2O, CO2, CH4, O2) diluted in N2 (so no concentrated system).
My idea was to use Species Transport to calculate the concentrations of every species in every element and to use Darcy to calculate the resulting changes in pressure and with that the velocity that I then feed back to Species Transport.
There is no other ready made function in COMSOL that I know of to calculate the velocity from differences in pressure (for 1D porous media). Using Darcy to do that does not necessarily overdefines it, does it?
PS.
I of course remember your comment in my other topic
www.comsol.com/community/forums/general/thread/25998/
and in that case I eventually found out that I made a mistake in one of my functions, since I corrected that, the integrals now give the right answer, it was not necessarily overdefined there. I thought it might help to know :)
Please login with a confirmed email address before reporting spam
In your simple model refining the mesh does not help, and the oscillations nearly disappear after some time so they could be caused by a sudden change in loading in a transient analysis. If you try ramping up the C0,c2 inflow condition in the Species Transport and also the velocity condition in Darcy’s Laws with a fast acting step function, the solution becomes less oscillatory. I don’t know how much of that applies to your other detailed model. One more factor to investigate!
Please keep us posted if you resolve the issue.
Dear Nagi,
Thank you for your input.
It seems you are right. I did more tests and must conclude that apparently the high Peclet number does not seem to be the problem in this particular case. I now used a stepfunction to increase the inflow velocity from 0 to v_in and this did indeed resolve the issue with the oscillations in the v plots.
However, I now noticed that the plot for the velocity is just a horizontal line (the v_in that I used as input), which means that apparently the inflow of mass (my drying reaction, input of H2O in gas stream) does not have any impact at all on the velocity. I dearly hope I did not miss anything fundamental or something, but it looks to me as if the increase in concentration/density does not reflect a change in pressure/velocity in this case, which seems wrong to me. Maybe the definitions of pressure in COMSOL are confusing me again.
[To all]:
The water particle source increases the total amount of particles in the gas flow so, since mass is added, I would expect the Darcy velocity dl.u to increase. The mass increase increases the pressure, but since it is an open system, the velocity should increase instead, right? I get really different results for different settings for Darcy flow so I want to calculate it myself, but I am confused.
Ideas gas law states: PV=nRT, which leads to: PM=RhoRT. But of course this is a flow system and not a closed system, so actually some kind of dynamic pressure should be used.
According to Bernoulli P1 + 1/2 * rho_1 * v_1^2 (no change in height) = P2 + 1/2 * rho_2 * v_2^2 , but Bernoulli can only be used in steady state right (which this isn't)? And only at constant T as well, right? I calculated the total (absolute) pressure before as P_total=SUM[C_i*R*T] but that seems wrong (way too high), is it wrong to do this? How should I do this?
If I set Darcy density to 'self defined' rho using rho_total=SUM[rho_i] I can't set the pressure and temperature. But if I set the density to 'Ideal Gas', I get the P and T option.
(see Darcy_Interface_self_defined_rho_vs_Ideal_Gas.jpg).
This leads to very different results.
(see rho_g_tot_vs_x.jpg, rho_g_tot is green, dl.rho is red).
Why is that different, according to (rewritten) Ideal Gas Law, rho_i=Ci*Mi, right?
So why the difference? How then does Darcy calculate the density? I have checked the manuals but can't find what I am looking for.
For the velocity this means that with the self defined rho, v=0.1 constantly, while with the Ideal Gas option, the velocity goes up to 2 m/s.
This could be right, but I have no idea to check if it is. Any help would again be greatly appreciated. Many thanks for reading.
Kind regards,
Raymond
Please login with a confirmed email address before reporting spam
1. rho_i = c_i*M_i is not "according to the ideal gas law", it's just a relation.
2. When I open your models, I can't access the equation view in the physics, even if my Options\Preferences... are set to show all/expand all. Is it also your case? If yes, you can't access the detailed equations and you can't see that Darcy's module calculates the density from ideal gas law in this way: rho = p/(RT) where R is the specific gas constant in J/(kg.K). Maybe it's a little bug and you should rebuild your models from scratch.
3. Regarding the difference between the ideal gas density calculated by Darcy's module and the density calculated manually by rho = sum(c_i*M_i) --> been there, done that! And I had weird results just like you. This is when we began to realize that solving the N species equations + the continuity equation made the problem overdefined.
In a way, you just proved it yourself: you used the sum of your species concentrations to calculate the gas density. From this density, you know that you can calculate the pressure with the ideal gas law. Now that you have the pressure, you just have to calculate the Darcy's velocity in a list a variable. The Darcy's velocity calculation is very simple, it's the same calculation as, for example, the heat flux that Comsol calculates when you plot a heat flux (q" = -k*d(T,x). You don't need a module to do that for you, and remember that the Darcy's module solves a form of the continuity equation, expressed with the pressure as the dependent variable, then it calculates the Darcy's velocity separately from the pressure.
An excellent reference regarding the sum of the N species equations is: Transport Phenomena 2nd Ed. by Bird, Stewart and Lightfoot, section 19.1: The equations of continuity for a multicomponent mixture.
Please login with a confirmed email address before reporting spam
I am not sure if my situation fits under this topic but I guess I'm having the same problem. I tried to solve the following dimensionless equations:
u1t=u1x+u2-k1*u1
u2t=c*u2xx-(k2+1)*u2+k3*u3+k1*u1
u3t=-u3x+k2*u2-k3*u3
(u is the concentration) with the boundary conditions:
u1(5,t)=0
u2(0,t)=0
u2(5,t)=0
u3(0,t)=1-heaviside(t-0.1)
and initial conditions:
u1(x,0)=0, u2(x,0)=0, u3(x,0)=0
They are describing material transport in a cell where advection and diffusion take place. I tried coefficient form but all I got is oscillations and negative concentration results. I already increased the diffusion coefficient, reduced the mesh size in order to make peclet number smaller than 1 and picked smoothed version of heaviside function but they didn't give any meaningful results. I also couldn't add stabilization part even though I chose it under the show and expand sections parts. Maybe it's because I'm using Comsol version 4.2 with classkit licence.
I am pretty much stuck and I don't know if I defined the coefficients in the wright way. I don't have any other idea to make this model work. Could someone take a quick look and see what could be wrong?
Thank you very much.
Attachments:
Please login with a confirmed email address before reporting spam
first I would rather use a step function rather than the "flsmhs()" on BC(1) that has two nice overshoots certainly playing you games
Try to plot your function, create an Definition analytical function and plot it, always best to check your equations ;)
Then when you have such a step on a boundary, you should use a finer mesh along this boundary, as your diffusion equations makes a rather steep gradient there and your are not resolving it correctly with your coarse time steps, this leads easily to numerical issues and "negative" overshoot values.
There is still something funny with your model, it crashes my COMSOL :)
--
Good luck
Ivar
Please login with a confirmed email address before reporting spam
I have access to Comsol via remote desktop application and it was crashing from time to time but I guessed that was normal. I don't know why it's doing the same thing to you. Sorry about that.
I played around the mesh size and timestep .Only the amount of oscillations changed and the results were not stable. I also changed flsmhs to triangle, rectangle, and flc1hs step functions and they didn't solve the problem, either. I am not sure how I should change size of transition zone or scale for the step functions.
Lastly, I used time discrete solver and it returned better results than time dependent solver. Is it right to use time discrete solver in my case? I didn't use any bdf or prev operator anywhere in my model.
I forgot to mention that reference guide lists some stabilization techniques for this kind of problem but even though I choose stabilization option under the show menu, I can't see them in my module. Are they not available in pde module or is it something specific to the version I'm using (version 4.2 with classkit licence)?
Thank you very much.
nihat
Please login with a confirmed email address before reporting spam
I would try out the change the step function to have a smoth transition and then steepen it to get closer to your initial step function (open up all tabs there you will see the step function transition time field)
There are a quite some discussion on stabilisation techniques around in the doc. I'm not that familiar with your type of imulations, and as I can get it runningI cant really say much more, but there are others out here with more experience in this field than me ;)
--
Good luck
Ivar
Please login with a confirmed email address before reporting spam
The time discrete solver is a good solver to use for transient problems even if you don’t use the prev or bdf operators.
Nagi Elabbasi
Veryst Engineering
Note that while COMSOL employees may participate in the discussion forum, COMSOL® software users who are on-subscription should submit their questions via the Support Center for a more comprehensive response from the Technical Support team.