Skip to content
1 contributor

Users who have contributed to this file

115 lines (87 sloc) 4.63 KB
// This example shows the resolution of a 3D magnetodynamic problem of magnetic induction using the
// so called a-v formulation. An aluminium tube is surrounded by 5 turns of a thick copper wire (1cm radius).
// Because of the alternating voltage (50 Hz) applied to the ends of the wire, eddy currents appear
// in the conducting aluminium tube. In this example the skin effect in the thick wire can also be observed.
// In order to remove the singularity (that comes from the magnetic equations) in the generated algebraic matrix
// a gauge condition is used. For that a spanning tree is created on the whole mesh, starting the growth on the
// region where the magnetic potential vector field 'a' will be constrained (here the domain boundary).
// This example was adapted from, and validated with an example developed for the GetDP
// software (Patrick Dular and Christophe Geuzaine).
#include "sparselizardbase.h"
using namespace mathop;
void sparselizard(void)
// The domain regions as defined in 'inductionheating.geo':
int coil = 1, tube = 2, air = 3, coilskin = 4, tubeskin = 5, vin = 6, vout = 7, domainboundary = 8;
// Load the mesh file. It can include curved elements.
mesh mymesh("inductionheating.msh");
// Define extra regions for convenience:
int conductor = regionunion({coil,tube});
int wholedomain = regionunion({conductor,air});
domainboundary = regionunion({domainboundary,vin,vout});
parameter mu, sigma;
// Define the magnetic permeability mu [H/m] everywhere (the materials here are not magnetic):
mu|wholedomain = 4*getpi()*1e-7;
// Conductivity of the copper coil and aluminium tube [S/m]:
sigma|coil = 6e7;
sigma|tube = 3e7;
// Set the working frequency to 50 Hz:
// Define a spanning tree to gauge the magnetic vector potential (otherwise the matrix to invert is singular).
// Start growing the tree from the regions with constrained potential vector (here the domain boundary):
spanningtree spantree({domainboundary});
// Use nodal shape functions 'h1' for the electric scalar potential 'v'.
// Use edge shape functions 'hcurl' for the magnetic vector potential 'a'.
// A spanning tree has to be provided to field 'a' for gauging.
// Since the solution has a component in phase with the electric actuation
// and a quadrature component we need 2 harmonics at 50Hz
// (harmonic 1 is DC, 2 is sine at 50Hz and 3 cosine at 50Hz).
field a("hcurl", {2,3}, spantree), v("h1", {2,3});
// Gauge the vector potential field on the whole volume:
// Select adapted interpolation orders for field a and v (default is 1):
a.setorder(wholedomain, 0);
// Put a magnetic wall (i.e. set field a to 0) all around the domain (no magnetic flux can cross it):
// Also ground v on face 'vout':
// Set v to 1V on face 'vin' for the in-phase component and to 0 for the quadrature component:
v.harmonic(2).setconstraint(vin, 1);
// Define the weak magnetodynamic formulation:
formulation magdyn;
// The strong form of the magnetodynamic a-v formulation is
// curl( 1/mu * curl(a) ) + sigma * (dt(a) + grad(v)) = js, with b = curl(a) and e = -dt(a) - grad(v)
// Magnetic equation:
magdyn += integral(wholedomain, 1/mu* curl(dof(a)) * curl(tf(a)) );
magdyn += integral(conductor, sigma*dt(dof(a))*tf(a) + sigma* grad(dof(v))*tf(a) );
// Electric equation:
magdyn += integral(conductor, sigma*grad(dof(v))*grad(tf(v)) + sigma*dt(dof(a))*grad(tf(v)) );
// Generate the algebraic matrices A and vector b of the Ax = b problem:
// Get the x solution vector:
vec sol = solve(magdyn.A(), magdyn.b());
// Update the fields with the solution that has just been computed:
a.setdata(wholedomain, sol);
v.setdata(conductor, sol);
// Write the magnetic induction field b = curl(a) [T], electric field e = -dt(a) - grad(v) [V/m] and current density j [A/m^2]:
expression e = -dt(a) - grad(v);
curl(a).write(wholedomain, "b.pos");
e.write(conductor, "e.pos");
(sigma*e).write(conductor, "j.pos");
v.write(conductor, "v.pos");
// Code validation line. Can be removed:
double minj = norm(sigma*(-2*getpi()*50*a.harmonic(2) - grad(v.harmonic(3)))).min(tube,4)[0];
std::cout << (minj < 216432 && minj > 216430);
int main(void)
return 0;
You can’t perform that action at this time.