I have a question on 2D static magnetic field simulation. My goal is to compute the magnetic field distribution around experimental setup in our precision physics experiment. The experimental setup is immersed in a uniform magnetic field pointing along a certain direction.
First, I am trying to solve a problem for which analytical solution exist so that I can compare Elmer results to the analytical solution. Depending on the success, I would like to move towards more complicated geometries.
All quantities are in SI units below.
Problem: There is circular surface (0.1 in radius) made of Iron of relative permeability 5000. It is placed in vacuum, which I model with another surrounding spherical surface of radius 1, with relative permeability of 1. Then I set the vector potential using MATC method such that there is uniform magnetic field of 1 Tesla in x direction.
The SIF file is as follows:
Code: Select all
Header
CHECK KEYWORDS Warn
Mesh DB "." "."
Include Path ""
Results Directory ""
End
Simulation
Max Output Level = 5
Coordinate System = Cartesian
Coordinate Mapping(3) = 1 2 3
Simulation Type = Steady state
Steady State Max Iterations = 1
Output Intervals = 1
Timestepping Method = BDF
BDF Order = 1
Solver Input File = case.sif
Post File = case.vtu
End
Constants
Gravity(4) = 0 -1 0 9.82
Stefan Boltzmann = 5.67e-08
Permittivity of Vacuum = 8.8542e-12
Boltzmann Constant = 1.3807e-23
Unit Charge = 1.602e-19
End
Body 1
Target Bodies(1) = 2
Name = "Body Property 1"
Equation = 1
Material = 2
End
Body 2
Target Bodies(1) = 1
Name = "Body Property 2"
Equation = 1
Material = 1
End
Solver 2
Equation = MgDyn2DPost
Procedure = "MagnetoDynamics2D" "BSolver"
Exec Solver = Before Saving
Stabilize = True
Bubbles = False
Lumped Mass Matrix = False
Optimize Bandwidth = True
Steady State Convergence Tolerance = 1.0e-5
Nonlinear System Convergence Tolerance = 1.0e-7
Nonlinear System Max Iterations = 20
Nonlinear System Newton After Iterations = 3
Nonlinear System Newton After Tolerance = 1.0e-3
Nonlinear System Relaxation Factor = 1
Linear System Solver = Iterative
Linear System Iterative Method = BiCGStab
Linear System Max Iterations = 500
Linear System Convergence Tolerance = 1.0e-10
BiCGstabl polynomial degree = 2
Linear System Preconditioning = ILU0
Linear System ILUT Tolerance = 1.0e-3
Linear System Abort Not Converged = False
Linear System Residual Output = 10
Linear System Precondition Recompute = 1
End
Solver 1
Equation = MgDyn2D
Procedure = "MagnetoDynamics2D" "MagnetoDynamics2D"
Variable = Potential
Exec Solver = Always
Stabilize = True
Bubbles = False
Lumped Mass Matrix = False
Optimize Bandwidth = True
Steady State Convergence Tolerance = 1.0e-5
Nonlinear System Convergence Tolerance = 1.0e-7
Nonlinear System Max Iterations = 20
Nonlinear System Newton After Iterations = 3
Nonlinear System Newton After Tolerance = 1.0e-3
Nonlinear System Relaxation Factor = 1
Linear System Solver = Iterative
Linear System Iterative Method = BiCGStab
Linear System Max Iterations = 500
Linear System Convergence Tolerance = 1.0e-10
BiCGstabl polynomial degree = 2
Linear System Preconditioning = ILU0
Linear System ILUT Tolerance = 1.0e-3
Linear System Abort Not Converged = False
Linear System Residual Output = 10
Linear System Precondition Recompute = 1
End
Equation 1
Name = "Coupled equations"
Active Solvers(2) = 2 1
End
Material 1
Name = "Air_vacuum"
Relative Permeability = 1.0
End
Material 2
Name = "Iron"
Relative Permeability = 1.0E6
End
Boundary Condition 1
Target Boundaries(1) = 2
Name = "Vector_pot_uniform_flux"
$B_applied = 1.0
Potential = Variable Coordinate 2
Real MATC "B_applied * tx"
End
The left figure shows the 2D profile for B-field strength (magnetic flux density). On the right I plot the components of B-field computed by Elmer: B = (B_x, B_y, B_z) along the x-axis. I used the "plot over line" filter in Paraview and saved the screenshot.
As one expects, B_y and B_z are identically zero. The B_x profile is correct in the sense that:
(i) it has the right values at the boundaries x=+1 and x=-1,
(ii) it has the correct symmetry about x=0, and
(iii) it is flat within the iron sphere/surface.
However, the analytical solution for this problem dictates that within the iron sphere/surface, B_x should have a value ~2.99 times the uniform field at the outer boundaries (1 Tesla). Instead what you see in the graph is that it is about 1.975 Tesla. I derived the analytical expression which, as far as I know, is exact because of the symmetry of the problem. It is analogous to the corresponding electrostatics problem from Jackson's Classical Electrodynamics textbook (3rd edition), Section 4.4, example associated with figure 4.6. The only change is that you have to replace dielectric constant by relative permeability in equations 4.53 to 4.55. Sections 5.10 and 5.11 also give the same result.
In the limit of relative permeability ----> infinity, the maximum B_x inside the iron sphere/surface is 3 times the field at the far-field boundary. I ran simulations with relative permeability from 1 all the way to 10^6. For relative permeability=1, I recover the uniform vacuum field distribution as expected. The high permeability limit however, does not agree with analytic solution.
I also tried the same version with high mesh density in the region of the iron sphere/surface which gives the same result, as shown in the attached figure.
Let me know if anyone has thoughts on the mismatch of peak B_x inside the iron sphere/surface, between the analytical vs simulated results. Or if you have other suggestions that I could try.
Thank you for reading this post.