Elmer Magnetics 3D comparison with 2D Axisymmetric Lumpy

[kataja]

Hi,

note also that integrating the normal component of the variable "magnetic flux density e" over some surface gives you exactly the right line integral of vector potential over its boundary (see SaveData module and convective flux operator). Internally, the "e" field is of first order (or higher) even though its visualization in paraview is piecewise constant. So even for second order vector potential field the "magnetic flux density e" is optimal.

Cheers,
Juhani

[crobar]

Hi,

note also that integrating the normal component of the variable "magnetic flux density e" over some surface gives you exactly the right line integral of vector potential over its boundary (see SaveData module and convective flux operator). Internally, the "e" field is of first order (or higher) even though its visualization in paraview is piecewise constant. So even for second order vector potential field the "magnetic flux density e" is optimal.

Cheers,
Juhani

Yes, but sometimes it is much easier to integrate the closed loop line integral than the surface normal dot product integral. This is the actual reason for wanting the vector potential. Imagine a surface where all points on the boundary are not coplanar for instance, this is often the case for generator coils. In this case finding the correct surface and surface normal is not trivial.

[kataja]

Hi!

You are right, sometimes it is.

However, with the convective flux operator one doesn't need to define the normal. It is automatically deduced. Moreover, integrating over dim-1 (surfaces in 3D and curves in 2D) domains, where dim is the ambient space dimension, can be naturally done with elmer because dim-1 domains are handled in "boundary condition" subsections and they can be named etc. whereas there is next to no (up to my knowledge) framework for dealing with dim-2 objects in elmer.

So in terms of workflow, defining the flux computation surfaces (how complicated they may be) during geometry/mesh definition phase will lead in huge simplification in the total workflow, because in this phase, we usually utilize some GUI tool with which the flux surfaces can be accurately and reliably identified. Internally, we could perform line integrals over the boundaries of these identified surfaces, but that will not lead in any additional accuracy.

What I'm trying to say that even though generally computing the flux with line integrals might seem easy. In terms of clean and reliable workflow, it is irrelevant.

Cheers,
Juhani

[crobar]

What I'm trying to say that even though generally computing the flux with line integrals might seem easy. In terms of clean and reliable workflow, it is irrelevant.

Cheers,
Juhani

Thanks, I appreciate your explanation and advice, but I'm still not sure how to easily apply it for a typical generator coil topology. Would you consider having a look at my post here:

viewtopic.php?f=3&t=4311&p=15365&hilit= ... 5ff#p15365

Which has an image of a typical coil. And suggesting how you would approach this problem and what workflow you would choose to get the double integral over the 3D coil shape in Elmer using the surface integral method on the flux density you suggest rather than multple line integrals of the vector potential?

Even suggesting how you would do it for the coil shown in the first post of this thread (the tube wrapped around the smaller cylinder) would be really helpful. I would need to integrate the flux going through the tube normal to the hole in the tube, and do this also over its axial length. Can I do this without defining a large number of surfaces and boundaries in preprocessing?

[kataja]

I'm a bit confused over which surface the flux is to be calculated?

It seems, however, somehow sensible to only consider line integrals of the vector potential instead of surface integrals because only the boundary line and vector potential on it matter in terms of flux..

-Juhani

[crobar]

I'm a bit confused over which surface the flux is to be calculated?

It seems, however, somehow sensible to only consider line integrals of the vector potential instead of surface integrals because only the boundary line and vector potential on it matter in terms of flux..

-Juhani

Well, in the example in the first post of this thread, the tube represents a coil of wire with multiple turns wrapped concenticily around a little cylinder of iron (which is near the cylindrical magnet). If the coil was just a single turn, we would want the integral of the flux going through the middle of this loop of wire. In practice the coil is multiple turns packed together into a tube shape. The flux passing through every turn will be slightly different due to the different diameter, so a slightly different surface would really be needed for every turn, or at least an integral relflecting this. The integral also has to be repeated in the axial direction as the loops are also stacked on top of one another.

This complexity is why the vector potential is normally used. I know FEMM for instance (which is 2D only) computes this flux linkage calculation by integrating the vector potential over the coil cross section, (the rectangular cross-section of the tube in this case) dividing by the area and then multiplying by the number of turns in the coil. I want to do the same in 3D (where it's a dot product integral) and have already written the code to do it. However, I'd obviously take the easier workflow if it exists.

Does that make things clearer at all?

I am reasonably sure that the commercial solvers must do it using the vector potential as in at least one you can define 'virtual' coils which don't actually affect the meshing and just represent regions where the integration will take place in post-procesing. Obviosly these virtual coils can't then be carrying a current. I suppose this could still be doing an integration with the flux density but they don't mesh any special features to achieve this as far as I can tell.

[raback]

Hi

Inspired by the discussion I started to work in baby steps towards being more easily to do some line integrals. I've looked today on the SaveLine with the intent to 1) allow saving of also edge element solution on the line and 2) extending the lines to lines in cylindrical coordinates i.e. arbitrary circular paths. If this would work in SaveLine it would be quite little additional work to intagrate over the line.

Now the question is whether this would be usefull to your purposes? Is a circular closed path something that could be applied? Currently the circle is defined by one point, the normal of the cylinder axis, and a radius i.e. 7 parameters.

It could be possible to fit the circles with some heuristics since the coils certainly are not always simple. Also for this reason there is the CoilSolver which tries to automatically find the center of the coil and then generates a potential the gradient of which is always pointing to the direction of the current.

-Peter

[crobar]

Hi
Now the question is whether this would be usefull to your purposes? Is a circular closed path something that could be applied? Currently the circle is defined by one point, the normal of the cylinder axis, and a radius i.e. 7 parameters.
-Peter

Well, not exactly. The circular coil example is only a very trivial case I chose in order to do a comparison. In real life problems I deal with the coil shape is rarely a perfect circle, and can be quite complex. I only chose the example in this thread because it allowed a direct comparison with an axisymmetric finite element solution, not because it is a realistic problem.

In my own attempts to calculate the integrals I have just chosen a 'trajectory' which starts and ends at the same point. I then also created a 2D mesh of a surface representing a coil cross-section (it is reasonable to assume the coil cross-section shape stays constant in most cases).

I then aligned the normal of this surface to the tangent of the trajectory and integrated the dot product on each triangle of this surface, and repeat this moving along the closed loop summing up the results. This is still experimental but gave reasonable results at first glace on the example problem in this thread.

Actually, for me, the main thing is not sorting out doing the line integrals in Elmer really (although a rigorous built-in way to do them would be nice), but rather making sure we can trust the vector potential results.

Richard

[macek]

I've made a similar comparsion between FEMM, Elmer2D and Elmer3D.

The geometry consits of an aluminum strand, stainless steel shell and a coil. All materials have magnetic permeability mu=1. Frequency is 15 Hz.

Meshes for Elmer3D and 2D. Mesh for 2D was much finer.

Results on a line in the radial direction, where y is at the coil centre, matches for all solvers. For 3D case I've taken lines at each degree and then averaged them. I've use point data (e.g. 'magnetic flux density re')


But there is a difference between Elmer 3D solver and the 2D solvers (FEMM, Elmer), but only in the strand (in the conductive medium). Plots at different radii confirm this.



Plots for r=0.5R, r=1.0R, r=1.2R, r=1.4R, where R is the strand radius.

First I thought that there is a problem in the mesh (bottom strand was meshed with triangles and the extruded, resulting elements were prisms), but then I created a mesh with tetrahedrons in the strand (shown above) and got the same results. Also I run the simulation with lower tolerance and still got the same results.

Results matches, if I took element data (e.g. 'magnetic flux density re e'), but only at the strand border. So I guess the results matches, because the data is from the air not from the strand.

The following two pictures show, that the problem is most probably at the shell-strand contact.

And we can also see two spikes for the r=R figure above.
Then I compared the results, where I've set the shell's electric conductivity to zero (but assigned different material than for the air) and got the same problem.



How to overcome this problem? Is this some kind of element numbering issue?

This is my AV solver code:

Code: Select all

Solver $solver_AV
  Equation = "MgDyn"
  Procedure = "MagnetoDynamics" "WhitneyAVHarmonicSolver"
  Variable = AV[AV re:1 AV im:1]

  Angular Frequency = $ omega

  $ jfix_tol = 5.0e-9
  $ jfix_max_iter = 1000
  $ lin_tol = 1.0e-6
  $ lin_max_iter = 8000
  $ nonlin_tol = 1.0e-6
  $ nonlin_max_iter = 1

  Exec Solver = Always  
 
  Linear System Timing = True
  Linear System Timing Cumulative = True
  Solver Timing = True
  Solver Timing Cumulative = True
  Fix Input Current Density = Logical True

  Stabilize = True
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = $ nonlin_tol

  JFix: Linear System Convergence Tolerance = $ jfix_tol
  JFix: Linear System Max Iterations = $ jfix_max_iter
  JFix: Linear System Abort Not Converged = False
  JFix: Linear System Residual Output = 10
  JFix: Linear System Solver = Iterative
  JFix: Linear System Iterative Method = BiCGStabl
  JFix: BiCGstabl polynomial degree = 4
  JFix: Linear System Preconditioning = Diagonal !ILU2 only for nonparallel run !!!
  JFix: Linear System ILUT Tolerance = 1.0e-4
  JFix: Linear System Precondition Recompute = 1

  Linear System Symmetric = Logical True
  Linear System Convergence Tolerance = $ lin_tol
  Linear System Max Iterations = $ lin_max_iter
  Linear System Abort Not Converged = False
  Linear System Residual Output = 1
  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGstabl
  BiCGstabl polynomial degree = 4
  Linear System Preconditioning = Diagonal ! Only None or Diagonal preconditioning
  Linear System Precondition Recompute = 1
  
  Nonlinear System Convergence Tolerance = $ nonlin_tol
  Nonlinear System Max Iterations = $ nonlin_max_iter 
  Newton-Raphson Iteration = Logical True
End

And code for the B-field solver

Code: Select all

Solver $solver_MagField

  Equation = "MgDynPost"
  Procedure = "MagnetoDynamics" "MagnetoDynamicsCalcFields"
  Potential Variable = String "AV"

  Exec Solver = Always


  Angular Frequency = $omega
  Show Angular Frequency = Logical True

  Linear System Timing = True
  Linear System Timing Cumulative = True
  Solver Timing = True
  Solver Timing Cumulative = True

  Calculate Magnetic Field Strength = Logical True
  Calculate Electric Field = Logical True
  Calculate Current Density = Logical True
  Calculate JxB = Logical True

  Stabilize = True
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-7

  Linear System Symmetric = Logical True
  Linear System Convergence Tolerance = 1.0e-9
  Linear System Max Iterations = 10
  Linear System Abort Not Converged = False
  Linear System Residual Output = 1
  Linear System Solver = Iterative
  Linear System Iterative Method = CG
  BiCGstabl polynomial degree = 2
  Linear System Preconditioning = Diagonal
  Linear System Precondition Recompute = 1

  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 10
  Nonlinear System Relaxation Factor = 1

End

[kataja]

It seems that the mesh includes hexahedral and prism elements so using "Use piola transform = logical true" inside Solver $solver_AV -- End block might help. This will enrich the edge element basis to a "more complete" one.