main.cpp 5.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139
  1. // This code simulates the static magnetic field created by 5 permanent magnets (1cm x 1cm) placed side by side.
  2. // A steel region above the magnets is perturbing the magnetic field lines.
  3. //
  4. // The magnetic scalar potential is used to solve this problem. This is valid since there are no current sources.
  5. // The permanent magnets are treated as pre-magnetized pieces of non-magnetic material (since all magnetic domains are already fully oriented).
  6. #include "sparselizardbase.h"
  7. using namespace mathop;
  8. void sparselizard(void)
  9. {
  10. // The domain regions as defined in 'halbacharray.geo':
  11. int magnet1 = 1, magnet2 = 2, magnet3 = 3, magnet4 = 4, magnet5 = 5, steel = 6, air = 7, zeropotential = 8;
  12. // The mesh can be curved!
  13. mesh mymesh("halbacharray.msh");
  14. // Define new physical regions for convenience:
  15. int magnets = regionunion({magnet1, magnet2, magnet3, magnet4, magnet5});
  16. int wholedomain = regionunion({magnets, steel, air});
  17. // Nodal shape functions 'h1' for the magnetic scalar potential 'phi'.
  18. field phi("h1");
  19. // Use interpolation order 2 on the whole domain:
  20. phi.setorder(wholedomain, 2);
  21. // The magnetic scalar potential (just like the electric potential)
  22. // needs to be fixed at least at one node to have a well-posed problem.
  23. // Here it is forced to 0 at a permanent magnet corner point.
  24. phi.setconstraint(zeropotential);
  25. // Vacuum magnetic permeability [H/m]:
  26. double mu0 = 4*getpi()*1e-7;
  27. // Define the permeability in all regions:
  28. parameter mu;
  29. mu|air = mu0;
  30. mu|steel = 1000*mu0;
  31. mu|magnets = mu0;
  32. formulation magnetostatic;
  33. // In the absence of current sources the static Maxwell equations give:
  34. //
  35. // curl h = 0
  36. //
  37. // One can thus define a magnetic scalar potential 'phi' such that h = -grad(phi).
  38. //
  39. // Considering also that div(b) = 0 we get
  40. //
  41. // div(mu*(-grad(phi))) = 0
  42. //
  43. // with b = mu * h.
  44. //
  45. // In the permanent magnet region b = mu0 * (h + m),
  46. // i.e. the material is non-magnetic but it is pre-magnetized by the magnetization vector m [A/m].
  47. // We thus get:
  48. // div(mu0*(-grad(phi)) + mu0*m) = 0
  49. //
  50. // A permanent magnet with a [0, 800e3] A/m magnetization is considered.
  51. // The weak form corresponding to the above equations:
  52. magnetostatic += integral(wholedomain, -grad(dof(phi)) * mu * grad(tf(phi)) );
  53. // This is when all magnets are oriented in the y direction:
  54. magnetostatic += integral(magnets, array2x1(0, 800e3) * mu * grad(tf(phi)) );
  55. // This is in Halbach configuration (to maximise the magnetic field above the array):
  56. //magnetostatic += integral(magnet1, array2x1(-800e3, 0) * mu * grad(tf(phi)) );
  57. //magnetostatic += integral(magnet2, array2x1(0, -800e3) * mu * grad(tf(phi)) );
  58. //magnetostatic += integral(magnet3, array2x1(800e3, 0) * mu * grad(tf(phi)) );
  59. //magnetostatic += integral(magnet4, array2x1(0, 800e3) * mu * grad(tf(phi)) );
  60. //magnetostatic += integral(magnet5, array2x1(-800e3, 0) * mu * grad(tf(phi)) );
  61. magnetostatic.generate();
  62. // Solve the algebraic system Ax = b:
  63. vec sol = solve(magnetostatic.A(), magnetostatic.b());
  64. // Transfer the data from the solution vector to the phi field:
  65. phi.setdata(wholedomain, sol);
  66. // Write the magnetic scalar potential and the magnetic field with an order 2 interpolation.
  67. phi.write(wholedomain, "phi.pos", 2);
  68. norm(-grad(phi)).write(wholedomain, "hnorm.pos", 2);
  69. (-grad(phi)).write(wholedomain, "h.pos");
  70. // Evaluate the magnetic field 1.5cm above the center of the magnet array:
  71. std::vector<double> magfieldnorm = norm(-grad(phi)).interpolate(wholedomain, {0,0.02,0});
  72. std::cout << "Magnetic field 1.5cm above the array center: " << magfieldnorm[0] << " A/m" << std::endl;
  73. // Write 20 magnetic field lines starting on the top side of the magnet array:
  74. shape ln("line", -1, {-0.025,0.01,0, 0.025,0.01,0}, 20);
  75. (-grad(phi)).streamline(regionexclusion(wholedomain, magnets), "magneticfieldline.pos", ln.getcoords(), 0.01/5);
  76. // The MAGNETOSTATIC FORCE acting on the steel disk is computed below.
  77. // This field will hold the x and y component of the magnetic forces:
  78. field magforce("h1xy");
  79. // The magnetic force is projected on field 'magforce' on the steel region.
  80. // This is done with a formulation of the type dof*tf - force calculation = 0.
  81. formulation forceprojection;
  82. forceprojection += integral(steel, dof(magforce)*tf(magforce));
  83. // The force term integral must be performed on all elements in the steel region AND in the ELEMENT LAYER AROUND.
  84. // This is obtained by considering the elements on the whole domain 'wholedomain'. Since the force must only be
  85. // calculated on the steel disk, the test functions are only selected on the steel (this is done with tf(..., steel)).
  86. forceprojection += integral(wholedomain, - predefinedmagnetostaticforce(tf(magforce, steel), -grad(phi), mu));
  87. forceprojection.generate();
  88. vec solf = solve(forceprojection.A(), forceprojection.b());
  89. magforce.setdata(wholedomain, solf);
  90. std::cout << "Force acting on the steel disk (fx, fy) = (" << compx(magforce).integrate(steel, 5) << ", " << compy(magforce).integrate(steel, 5) << ") N per unit depth" << std::endl;
  91. // Code validation line. Can be removed.
  92. std::cout << (magfieldnorm[0] < 64963.8 && magfieldnorm[0] > 64963.6);
  93. }
  94. int main(void)
  95. {
  96. SlepcInitialize(0,{},0,0);
  97. sparselizard();
  98. SlepcFinalize();
  99. return 0;
  100. }