main.cpp 5.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144
  1. // This example shows the time-simulation of a 2D axisymmetric high-temperature superconductor (HTS):
  2. // a superconducting tube is used to shield its inside from an externally applied magnetic field
  3. // that increases linearly over time. For the simulation the so called a-v magnetodynamic formulation
  4. // is used. The v part can be disregarded in the configuration of this example.
  5. //
  6. // The tube height is 210 mm, its thickness is 1 mm and its inner radius is 10 mm.
  7. //
  8. // This example is based on and validated with the experimental and simulation work done by K. Hogan in his thesis
  9. // "Passive magnetic shielding with bulk high temperature superconductors under a non-uniform magnetic field".
  10. //
  11. // Details on the equations for HTS can also be found in "Numerical simulation of the magnetization of
  12. // high-temperature superconductors: a 3D finite element method using a single time-step iteration" (G. Lousberg)
  13. #include "sparselizardbase.h"
  14. using namespace mathop;
  15. void sparselizard(void)
  16. {
  17. // The physical region numbers as defined in the mesh file:
  18. int tube = 1, air = 2, domainskin = 3;
  19. setaxisymmetry();
  20. mesh mymesh("tubeinair.msh");
  21. // Define the whole domain for convenience:
  22. int wholedomain = regionunion({tube, air});
  23. // Magnetic permeability of the air:
  24. double mu = 4.0*getpi()*1e-7;
  25. // In this axisymmetric problem the magnetic vector potential can
  26. // be expressed as a vector with only a non-zero z component:
  27. field az("h1"), x("x");
  28. expression a = array3x1(0,0,az);
  29. // Put a magnetic wall (i.e. set field a to 0) all
  30. // around the domain (no magnetic flux can cross it):
  31. az.setconstraint(domainskin);
  32. // The magnetic vector potential a is rewritten as the sum of a known
  33. // source term and an unknown part to be determined: atotal = a + asource.
  34. // We want to apply an external induction field b of 1 mT in
  35. // the y direction that linearly increases over time:
  36. expression bsource = array3x1(0,1e-3*t(),0);
  37. // Since b = curl(a) (in cylindrical coordinates) a valid choice for a is:
  38. expression asource = array3x1(0,0,-1e-3*0.5*t()*x);
  39. // And dt(a)/dt is:
  40. expression dtasource = array3x1(0,0,-1e-3*0.5*x);
  41. // The strong form of the magnetodynamic a-v formulation is
  42. //
  43. // curl( 1/mu * curl(a) ) + sigma * (dt(a) + grad(v)) = js, with b = curl(a) and e = -dt(a) - grad(v)
  44. //
  45. // Here grad(v) is zero because of axisymmetry and y direction bsource and the electric field becomes:
  46. //
  47. expression e = -dt(a)-dtasource;
  48. // The conductivity of the high temperature superconductor is modeled using
  49. // a power law relation between the electric field and the current density:
  50. //
  51. // J = Jc/Ec^(1/n) * norm(E)^((1-n)/n) * E
  52. // |_______________________________|
  53. // sigma(E)
  54. //
  55. // where the critical current density Jc [A/m^2] and critical electric field Ec [V/m] are supposed
  56. // constant here and the J(E) relation is isotropic. An extra 1e-11 is added to never divide by 0.
  57. //
  58. double n = 30.0, Ec = 1e-4, Jc = 1e8;
  59. expression sigma = Jc/( pow(Ec,1.0/n) ) * pow( norm(e) + 1e-11, (1.0-n)/n );
  60. // Define the weak magnetodynamic formulation:
  61. formulation magdyn;
  62. magdyn += integral(wholedomain, 1/mu*( curl(dof(a)) + bsource ) * curl(tf(a)) );
  63. // Add an extra odd integration degree for convergence:
  64. magdyn += integral(tube, sigma*( dt(dof(a)) + dtasource )*tf(a), +1 );
  65. // Start the implicit Euler time resolution with an all zero solution and time derivative:
  66. vec initsol(magdyn), initdtsol(magdyn);
  67. impliciteuler eul(magdyn, initsol, initdtsol);
  68. // Tolerance on the nonlinear iteration:
  69. eul.settolerance(1e-4);
  70. // Use a 1 second timestep:
  71. double timestep = 1.0;
  72. // Run from 0 to 150 sec with a fixed-point nonlinear iteration at every timestep:
  73. std::vector<std::vector<vec>> sols = eul.runnonlinear(0.0, timestep, 150.0);
  74. // Field a is available at every timestep in sols[0] and its time derivative in sols[1]:
  75. std::vector<vec> asol = sols[0];
  76. std::vector<vec> dtasol = sols[1];
  77. // Loop over all time solutions:
  78. std::vector<double> bcenter(asol.size());
  79. for (int i = 0; i < asol.size(); i++)
  80. {
  81. // Set variable t() to the time of the current timestep:
  82. settime(timestep*i);
  83. // Make the data of vector asol[i] available in field az:
  84. az.setdata(wholedomain, asol[i]);
  85. // Write a and b = curl(a) to disk:
  86. norm(a).write(wholedomain, "norma" + std::to_string(i+1000) + ".vtu");
  87. norm(curl(a) + bsource).write(wholedomain, "normbtotal" + std::to_string(i+1000) + ".vtu");
  88. // Output the b induction field [T] at the tube center to assess the shielding effectiveness.
  89. // Interpolate at a x coordinate slightly away from 0 to avoid NaN issues:
  90. bcenter[i] = norm(curl(a) + bsource).interpolate(wholedomain, {1e-10,0,0})[0];
  91. std::cout << "b source " << timestep*i << " mT --> b tube center " << bcenter[i]*1e3 << " mT" << std::endl;
  92. }
  93. // Combine all timesteps in a ParaView .pvd file for convenience:
  94. std::vector<double> timestepvalues(151);
  95. for (int i = 0; i < 151; i++)
  96. timestepvalues[i] = timestep*i;
  97. grouptimesteps("normbtotal.pvd", "normbtotal", 1000, timestepvalues);
  98. // Write to file the field at the tube center for all timesteps:
  99. writevector("bcenter.csv", bcenter);
  100. // Code validation line. Can be removed:
  101. std::cout << (bcenter[asol.size()-1] < 0.03746 && bcenter[asol.size()-1] > 0.03744);
  102. }
  103. int main(void)
  104. {
  105. SlepcInitialize(0,{},0,0);
  106. sparselizard();
  107. SlepcFinalize();
  108. return 0;
  109. }