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.
Determining magnetic body force density - Possible bug in Comsol 3.5a
Posted Nov 24, 2009, 3:21 p.m. EST 4 Replies
Please login with a confirmed email address before reporting spam
Hi
I have attached a very simple problem that gives me erroneous results when compared to analytical solutions. I have looked at the COMSOL documentation and examples for this and it is insufficiently explained and could be a serious bug.
The problem involves determining the magnetic body force density given by F=mu0*(Magnetization_vector DOT GRAD) Hvector
COMSOL deals with calculating total force densities (using stress tensors and virtual work) but I want to calculate a spatial distribution of the force density. There is one example where they solve for the body force density in the example "Magnetic drug targeting in cancer therapy" . However when I try to replicate a simple enough analytical problem I get an incorrect result.
For an infinitely long cylinder of mu_r>1 in a uniform transverse field will have a uniform magnetic field inside because the demagnetizing factors are equal in all directions and are equal to 1/2. (Hinside=Houtside-(1/2)*Minside) .
This can be easily modelled in COMSOL in 2D and I have attached a file BodyForce.mph which outlines a 2D circle (representing an infinitely long cylinder in the z direction) in a uniform x directed magnetic field. The volume body force density in the x direction is given by the expression FMx. Plotting FMx vs x from x=-1 to 1 gives you a non-zero result - which is incorrect.
If I reduce the model to a 1 region problem like in BodyForce_oneregion.mph. The results are better but there is still a large enough non-zero FMx when it should be on the order of the precision of the machine since it should be 0.
Any help on this matter will be greatly appreciated and will help many people trying to solve spatial electromagnetic force calculations as there is only one example that outlines this.
Best
Shahriar
I have attached a very simple problem that gives me erroneous results when compared to analytical solutions. I have looked at the COMSOL documentation and examples for this and it is insufficiently explained and could be a serious bug.
The problem involves determining the magnetic body force density given by F=mu0*(Magnetization_vector DOT GRAD) Hvector
COMSOL deals with calculating total force densities (using stress tensors and virtual work) but I want to calculate a spatial distribution of the force density. There is one example where they solve for the body force density in the example "Magnetic drug targeting in cancer therapy" . However when I try to replicate a simple enough analytical problem I get an incorrect result.
For an infinitely long cylinder of mu_r>1 in a uniform transverse field will have a uniform magnetic field inside because the demagnetizing factors are equal in all directions and are equal to 1/2. (Hinside=Houtside-(1/2)*Minside) .
This can be easily modelled in COMSOL in 2D and I have attached a file BodyForce.mph which outlines a 2D circle (representing an infinitely long cylinder in the z direction) in a uniform x directed magnetic field. The volume body force density in the x direction is given by the expression FMx. Plotting FMx vs x from x=-1 to 1 gives you a non-zero result - which is incorrect.
If I reduce the model to a 1 region problem like in BodyForce_oneregion.mph. The results are better but there is still a large enough non-zero FMx when it should be on the order of the precision of the machine since it should be 0.
Any help on this matter will be greatly appreciated and will help many people trying to solve spatial electromagnetic force calculations as there is only one example that outlines this.
Best
Shahriar
Attachments:
4 Replies Last Post Dec 8, 2009, 6:07 p.m. EST
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi
first to minimise confusion, I would use either the plot parameters contour Az, or use the streamline but with "plot type" "magnitude controlled", otherwise the streamlines are very assymetric
Then I beleive it's correct, in the sens that you are observing FEM precision limitations, you need a very very fine mesh along the boundary, but even then you cannot avoid quite some gradients change over an element, and you have a "mu" of only 100, with soft iron >2000 it's even worse!.
jsut try to plot the different variables used for the FMx calculations along the same line
The magnetic flux density along X is a nice parabole, with some border artefacts realted directly to the meshing
It is like stress concetration in structural analysis. I would use adaptive meshing, even if the L2 norm used by COMSOL is not very good to "just" increase the mesh at stress concentration (or magnetic boundaries), it add quite some nodes all over.
If you want better explanations, I would suggest to submit the case to Comsol support, if so pls report back its of general interest, there are many dicussions around force calculations and their precision
Good luck
Ivar
first to minimise confusion, I would use either the plot parameters contour Az, or use the streamline but with "plot type" "magnitude controlled", otherwise the streamlines are very assymetric
Then I beleive it's correct, in the sens that you are observing FEM precision limitations, you need a very very fine mesh along the boundary, but even then you cannot avoid quite some gradients change over an element, and you have a "mu" of only 100, with soft iron >2000 it's even worse!.
jsut try to plot the different variables used for the FMx calculations along the same line
The magnetic flux density along X is a nice parabole, with some border artefacts realted directly to the meshing
It is like stress concetration in structural analysis. I would use adaptive meshing, even if the L2 norm used by COMSOL is not very good to "just" increase the mesh at stress concentration (or magnetic boundaries), it add quite some nodes all over.
If you want better explanations, I would suggest to submit the case to Comsol support, if so pls report back its of general interest, there are many dicussions around force calculations and their precision
Good luck
Ivar
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi
Thanks for the response. I have actually made the mesh as fine as possible and have even used adaptive mesh refinement but it makes no difference. I have increased the boundary layer mesh as well and it makes no difference.
The example that I have posted should analytically give a 0 force spatially. I have done simple things like d(Hx_emqa,x) and the results give non-zero jumps at the boundaries and non-zero (higher order of magnitude than machine precision) for regions within the cylinder. The jumps at the boundary could be because of surface stress but those should be impulses and should not propagate into the subdomain of interest. Within the subdomain of interest the d(Hx_emqa,x) should be analytically 0 or 10^-15 in COMSOL.
Whats interesting is that COMSOL gives well behaved (or better) answers when it is a one region problem than when it is a two region problem, which means it could be the very nature of FEM.
Nonetheless, the result COMSOL gives is incorrect and it can affect several people trying to do force calculations or any calculation that involves a spatial derivative. I am working on a complicated model that requires calculating the magnetic body force density and inserting that into the navier stokes equation. If COMSOL cannot compute this simple case I dont think I can have confidence in any result it gives for more complicated problems that I work on.
I have submitted it to COMSOL support. I hope they will respond.
Best
Shahriar
Thanks for the response. I have actually made the mesh as fine as possible and have even used adaptive mesh refinement but it makes no difference. I have increased the boundary layer mesh as well and it makes no difference.
The example that I have posted should analytically give a 0 force spatially. I have done simple things like d(Hx_emqa,x) and the results give non-zero jumps at the boundaries and non-zero (higher order of magnitude than machine precision) for regions within the cylinder. The jumps at the boundary could be because of surface stress but those should be impulses and should not propagate into the subdomain of interest. Within the subdomain of interest the d(Hx_emqa,x) should be analytically 0 or 10^-15 in COMSOL.
Whats interesting is that COMSOL gives well behaved (or better) answers when it is a one region problem than when it is a two region problem, which means it could be the very nature of FEM.
Nonetheless, the result COMSOL gives is incorrect and it can affect several people trying to do force calculations or any calculation that involves a spatial derivative. I am working on a complicated model that requires calculating the magnetic body force density and inserting that into the navier stokes equation. If COMSOL cannot compute this simple case I dont think I can have confidence in any result it gives for more complicated problems that I work on.
I have submitted it to COMSOL support. I hope they will respond.
Best
Shahriar
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi
Well we will see, but I'm still not 100% convinced that the field should be absolutely flat, and you cannot avoid steps at intefarces, such as step you see around changes of ur or E (stress for structural)
I have notived often that our analytical solutions are already approximations (that we have forgotten about) and that often we see more subtile effects with detailed FEM analysis. In youre case it could well be border effects due to the round surface, perhaps one should try a square central part.
Ivar
Well we will see, but I'm still not 100% convinced that the field should be absolutely flat, and you cannot avoid steps at intefarces, such as step you see around changes of ur or E (stress for structural)
I have notived often that our analytical solutions are already approximations (that we have forgotten about) and that often we see more subtile effects with detailed FEM analysis. In youre case it could well be border effects due to the round surface, perhaps one should try a square central part.
Ivar
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
This is the response I got from COMSOL Support
"With numerical modelling, you cannot expect to get zero result at the
limit of machine precision. It is the shape of the mesh elements at the
interface that decide how large the difference is. In this particular case,
you also use a second order derivative in the expression of the body force.
Because you are using second order elements (The Az-field is expressed with
second order polynomials), the second order derivative is a constant and
not very accurate. When you need high accuracy of higher order derivatives,
it is better to increase the element order than refining the mesh.
When I increased the element order to five in your model, the resulting
force was in the order of 1e-7. To decide whether this is a high value or a
low value, you should disturb your model so you get a non-zero force. Then
you compare that value with the "zero" value you get. As an example, I came
across another setup where the "zero" force was in the order of 1e3 N with
a virtual work calculation, which is a quite large force. When disturbing
that model the resulting force was in the order of 1e7 N, and it is this
relative comparison that you should look at.
To increase the shape order, select all domains, and open the subdomain
settings. The go to the element page, and select Lagrange - Quintic."
After changing to Quintic elements, I find that the problem works in the 2D case.
When I move to simulating a 3D problem I get errors. I think its just inherent to the FEM
"With numerical modelling, you cannot expect to get zero result at the
limit of machine precision. It is the shape of the mesh elements at the
interface that decide how large the difference is. In this particular case,
you also use a second order derivative in the expression of the body force.
Because you are using second order elements (The Az-field is expressed with
second order polynomials), the second order derivative is a constant and
not very accurate. When you need high accuracy of higher order derivatives,
it is better to increase the element order than refining the mesh.
When I increased the element order to five in your model, the resulting
force was in the order of 1e-7. To decide whether this is a high value or a
low value, you should disturb your model so you get a non-zero force. Then
you compare that value with the "zero" value you get. As an example, I came
across another setup where the "zero" force was in the order of 1e3 N with
a virtual work calculation, which is a quite large force. When disturbing
that model the resulting force was in the order of 1e7 N, and it is this
relative comparison that you should look at.
To increase the shape order, select all domains, and open the subdomain
settings. The go to the element page, and select Lagrange - Quintic."
After changing to Quintic elements, I find that the problem works in the 2D case.
When I move to simulating a 3D problem I get errors. I think its just inherent to the FEM