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.
Convection-diffusion equation and/or advection equation: incorrect solutions
Posted Jul 19, 2013, 12:41 p.m. EDT 4 Replies
Please login with a confirmed email address before reporting spam
Hi,
I have been working with the convection-diffusion module of the COMSOL Multiphysics module (version 4.3a).
I want to study how a charge density (a cloud of electrons) drifts horizontally and diffuses in all directions.
My question divides into three parts:
1.) 2D convection and diffusion
I started with only 2 dimensions and I defined a block of Silicon of 1 mm height and 2 mm length. The max time is set to 15 ns. In the convection-diffusion module, I set f (source) = 0, d_a (damping) = 1, and "beta" (the drift velocity) to 1.4E5 m/s in the x direction and 0 in y. (This comes from supposing a voltage of 2000 V, and a electron mobility of 1400 cm^2/(V*s).)
The diffusion coefficient "c" I set to isotropic and I varied the value.
The initial conditions for the density "u" I put to a 2D Gaussian shape,
u = exp(-((sqrt((x-DelX)^2+(y-DelY)^2)))^2/radius^2),
where DelX = 0.01, DelY = 0.05, and radius = 0.01
"u" is defined as a "Space charged density" (C/m^3).
The mesh was set to "physics-controlled" and "extremely fine".
When I set the diffusion coefficient c = 3.6E-1 [m^2/s], I see the following:
At time = 0 ns, I get a nice round blob, as expected
(see file "2D_ConvDiff_ExtrFineMesh_Circle_h0_2mm_0ns.jpg")
At time > 0 ns (for instance 3.5 ns and 7 ns), this blob moves to the right, slightly diffused in all directions
However, there seems to be a strange ripple in the tail after the blob...
(see file "2D_ConvDiff_ExtrFineMesh_Circle_h0_2mm_3p5ns.jpg")
So, my question is: why do I get these ripples?
2.) 2D addvection
When I set the diffusion coefficient to 0 (or even a very small value like 3.6E-3 [m^2/s]), the solution seems so "explode".
Unrealistic large values for w appear in the geometry and the results does not make sense anymore.
As a matter of fact, with the diffusion, c, set to 0, the equation is actually equally to an "advection equation", where I expect the density shape to move horizontally from left to right without diffusion.
However, the solution always seems to "explode" into huge values for u (of the order of 1E18, whereas the maximum should actually be 1.0).
So, my question is: how come with diffusion set to 0, the calculation collapses in 2D?
3.) 1D addvection
Next, I tried to do the same (with diffusion set to 0) in 1 dimension.
The line-length is 0.2 mm, the drift velocity (beta) is 1.4E5 m/s and the maximum time is set to 1.5 ns.
The initial density "u" I set to a Gaussian shape:
u = exp(-(x-DelX)^2/radius^2),
where DelX = 0.04 mm (to make sure the Gauss is well within the bounds of the line) and radius = 0.01 mm.
Now, I get a result similar to case 1) (however, now I have diffusion set to 0).
At time = 0 ns, I get a nice Gaussian shape
(see file "1D_DriftOnly_ExtrFineMesh_h0_0p2mm_0ns.jpg")
At time > 0 ns (for instance 3.5 ns and 7 ns), the curve moves to the right.
However, again, there appear strange ripples after the tail of the curve
(see file "1D_DriftOnly_ExtrFineMesh_h0_0p2mm_3p5ns.jpg")
So, my question is: why do I get these ripples again?
Best regards,
Machiel Kolstein
I have been working with the convection-diffusion module of the COMSOL Multiphysics module (version 4.3a).
I want to study how a charge density (a cloud of electrons) drifts horizontally and diffuses in all directions.
My question divides into three parts:
1.) 2D convection and diffusion
I started with only 2 dimensions and I defined a block of Silicon of 1 mm height and 2 mm length. The max time is set to 15 ns. In the convection-diffusion module, I set f (source) = 0, d_a (damping) = 1, and "beta" (the drift velocity) to 1.4E5 m/s in the x direction and 0 in y. (This comes from supposing a voltage of 2000 V, and a electron mobility of 1400 cm^2/(V*s).)
The diffusion coefficient "c" I set to isotropic and I varied the value.
The initial conditions for the density "u" I put to a 2D Gaussian shape,
u = exp(-((sqrt((x-DelX)^2+(y-DelY)^2)))^2/radius^2),
where DelX = 0.01, DelY = 0.05, and radius = 0.01
"u" is defined as a "Space charged density" (C/m^3).
The mesh was set to "physics-controlled" and "extremely fine".
When I set the diffusion coefficient c = 3.6E-1 [m^2/s], I see the following:
At time = 0 ns, I get a nice round blob, as expected
(see file "2D_ConvDiff_ExtrFineMesh_Circle_h0_2mm_0ns.jpg")
At time > 0 ns (for instance 3.5 ns and 7 ns), this blob moves to the right, slightly diffused in all directions
However, there seems to be a strange ripple in the tail after the blob...
(see file "2D_ConvDiff_ExtrFineMesh_Circle_h0_2mm_3p5ns.jpg")
So, my question is: why do I get these ripples?
2.) 2D addvection
When I set the diffusion coefficient to 0 (or even a very small value like 3.6E-3 [m^2/s]), the solution seems so "explode".
Unrealistic large values for w appear in the geometry and the results does not make sense anymore.
As a matter of fact, with the diffusion, c, set to 0, the equation is actually equally to an "advection equation", where I expect the density shape to move horizontally from left to right without diffusion.
However, the solution always seems to "explode" into huge values for u (of the order of 1E18, whereas the maximum should actually be 1.0).
So, my question is: how come with diffusion set to 0, the calculation collapses in 2D?
3.) 1D addvection
Next, I tried to do the same (with diffusion set to 0) in 1 dimension.
The line-length is 0.2 mm, the drift velocity (beta) is 1.4E5 m/s and the maximum time is set to 1.5 ns.
The initial density "u" I set to a Gaussian shape:
u = exp(-(x-DelX)^2/radius^2),
where DelX = 0.04 mm (to make sure the Gauss is well within the bounds of the line) and radius = 0.01 mm.
Now, I get a result similar to case 1) (however, now I have diffusion set to 0).
At time = 0 ns, I get a nice Gaussian shape
(see file "1D_DriftOnly_ExtrFineMesh_h0_0p2mm_0ns.jpg")
At time > 0 ns (for instance 3.5 ns and 7 ns), the curve moves to the right.
However, again, there appear strange ripples after the tail of the curve
(see file "1D_DriftOnly_ExtrFineMesh_h0_0p2mm_3p5ns.jpg")
So, my question is: why do I get these ripples again?
Best regards,
Machiel Kolstein
Attachments:
4 Replies Last Post Sep 12, 2013, 9:40 a.m. EDT
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi,
It is not surprising that you have oscillations because your problem has a very high Peclet number and in fact with diffusion set to zero it is purely a hyperbolic equation best solved via lagrangian or particle based or any suitable methods. What you experience is a well-known numerical issue and have been addressed extensively in literatures (e.g. cece.ucf.edu/people/yeh/briefcv/index.asp).
Within COMSOL itself you could find schemes to minimise problems with high Peclet numbers. For example, you can use their transport in dilute species module which offers consistent and inconsistent stabilization (if you don't have this module, you can also find the provision in heat transfer module or laminar flow module). They also give an example of a 2D heat transfer solution with high Peclet number. Of course, you must be aware that use of such schemes would result in artificial dispersion. If you don't have any, then you are stuck with putting an artificial diffusion yourself as you have done and ignore the oscillations if it is not that critical.
Alternatively, you could try adaptive meshing or particle tracing module (considering that yours is purely advective) in COMSOL.
Suresh
It is not surprising that you have oscillations because your problem has a very high Peclet number and in fact with diffusion set to zero it is purely a hyperbolic equation best solved via lagrangian or particle based or any suitable methods. What you experience is a well-known numerical issue and have been addressed extensively in literatures (e.g. cece.ucf.edu/people/yeh/briefcv/index.asp).
Within COMSOL itself you could find schemes to minimise problems with high Peclet numbers. For example, you can use their transport in dilute species module which offers consistent and inconsistent stabilization (if you don't have this module, you can also find the provision in heat transfer module or laminar flow module). They also give an example of a 2D heat transfer solution with high Peclet number. Of course, you must be aware that use of such schemes would result in artificial dispersion. If you don't have any, then you are stuck with putting an artificial diffusion yourself as you have done and ignore the oscillations if it is not that critical.
Alternatively, you could try adaptive meshing or particle tracing module (considering that yours is purely advective) in COMSOL.
Suresh
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi,
Thanks for your fast and interesting answer. Unfortunately, the coming 3 weeks I will be out of office, so I only will be able to study your suggestions in more detail after that.
In any case, what I understand is that in the case of a purely advection equation, the convection+diffusion module is not a good option, because the Peclet number goes to infinity (although in 1D, it works a lot better than for 2D). I would have to check whether I can just use the ordinary PDF module in COMSOL for a pure advection equation, or whether COMSOL still implicitly would use a computational model where the Peclet number causes problems.
About the adaptive messing, are you referring to the "user-defined meshing"? I suppose that this would mean define an extremely finer mess in the places where the density is highest, right?
Thanks again and best regards,
Machiel Kolstein
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi,
I am not fully proficient with COMSOL, so I may be wrong, but in my opinion even if you were to use PDEs to construct your own convection-diffusion equation but without numerical treatment or appropriate choice of shape functions, you would end up with the same problem, if you use classical Galerkin FEM.
If I may suggest that you go through an example from COMSOL's model library (www.comsol.com/showroom/documentation/model/72/) in which they teach you how to use a PDE and weak form contribution to introduce streamline diffusion and discontinuous elements to stabilize a stationary convective only problem. Note however that no crosswind is addressed in this example. This will help you if you do not have any modules that I mentioned earlier.
I suggest that you look for other discussions in the forum which seem to address similar problem as yours. Also I faintly remember some publications in COMSOL conference in which some authors implemented their own stabilization techniques.
Adaptive meshing is an inbuilt feature.
Finally, the appropriateness of chosen methods to 2D and 3D must be examined before investing your time in coding.
Hope this helps.
Suresh
I am not fully proficient with COMSOL, so I may be wrong, but in my opinion even if you were to use PDEs to construct your own convection-diffusion equation but without numerical treatment or appropriate choice of shape functions, you would end up with the same problem, if you use classical Galerkin FEM.
If I may suggest that you go through an example from COMSOL's model library (www.comsol.com/showroom/documentation/model/72/) in which they teach you how to use a PDE and weak form contribution to introduce streamline diffusion and discontinuous elements to stabilize a stationary convective only problem. Note however that no crosswind is addressed in this example. This will help you if you do not have any modules that I mentioned earlier.
I suggest that you look for other discussions in the forum which seem to address similar problem as yours. Also I faintly remember some publications in COMSOL conference in which some authors implemented their own stabilization techniques.
Adaptive meshing is an inbuilt feature.
Finally, the appropriateness of chosen methods to 2D and 3D must be examined before investing your time in coding.
Hope this helps.
Suresh
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi Suresh,
Your replies and those of some other people have suggested various possible solutions and I am still investigating them. One question about the Pectlet number though.
The Peclet number is defined as the drift velocity times the grid (mesh) size divided by the diffusion rate. However, even when defining the Mesh as "extremely fine", the solution still shows the same fluctuations and values smaller than zero.
Is there any way, as a user, to set the Mesh size even smaller in COMSOL?
Cheers,
Machiel
Your replies and those of some other people have suggested various possible solutions and I am still investigating them. One question about the Pectlet number though.
The Peclet number is defined as the drift velocity times the grid (mesh) size divided by the diffusion rate. However, even when defining the Mesh as "extremely fine", the solution still shows the same fluctuations and values smaller than zero.
Is there any way, as a user, to set the Mesh size even smaller in COMSOL?
Cheers,
Machiel