load "iovtk" //load paraview lib verbosity=0; ////////////////////////////////////////////////////////////////////////////////////////////// // BENCHMARK PROBLEM ///////////////////////////////////////////////////////////////////////////////////////////// //--Geometry parameters---------------------------------- //mold real rmold = 25.255e-03; // mold ext radius real hmold = 34.74e-03; // mold height //air real rair = 200e-03; // cap height real hair = 100e-03; //spires real rSpireExt = 6e-03; //spire external radius real rSpireInt = 3.05e-03; //spire internal radius real airGap = 12.5e-03; real deltaY1 = 8.72e-03; // vertical air gap between top-mid spires 8.72e-03 real deltaY2 = 8.72e-03; // vertical air gap between mid-bot spires //14.25e-03 //--Mesh parameters--// real Lcmold = 40; real Lcskin = 40; real Lcair = 10; real Lcspire = 40; // Electromagnetic------------------------------------- // real Pw =53.04455934993124; // real Qs = 918103.2985457362; // sigma = 9 real Power = 99; real IrmsMax=85; real Irms=(Power/100)*IrmsMax; //////////////////////////////////////////////////////////////// //----------------- Materials parameters----------------------- //////////////////////////////////////////////////////////////// real mu0 = 4.*pi*1.e-7; //Vacuum magnetic permeability //--Mold-- real muMoldRelative = 1800; //relative permeability mu/mu0 real muMoldReal= muMoldRelative*mu0; real sigmaMold = 9.0e06; //--Coils-- real muSpireRelative =1; //relative permeability mu/mu0 real muSpireReal= muSpireRelative*mu0; real sigmaSpire = 58.7e06; //--Air-- real muAirRelative =1; //relative permeability mu/mu0 real muAirReal= muAirRelative*mu0; real sigmaAir = 0.0; //////////////////////////////////////////////////////////////////////////////////////////// // GEOMETRY & MESH START // ////////////////////////////////////////////////////////////////////////////////////////// //mold ext border topMold(t = 0,rmold){x=t;y=hmold/2;label=400;} border rightMold(t = hmold/2,-hmold/2){x=rmold;y=t;label=401;} border botMold(t = rmold,0){x=t;y=-hmold/2;label=402;} border leftMold(t = -hmold/2,hmold/2){x=0;y=t;label=403;} //spires real spireExtArea = pi*(rSpireExt)^2; real spireIntArea = pi*(rSpireInt)^2; real spireArea = spireExtArea-spireIntArea; //spire mid border spireMidExt(t=0, 2*pi){x=rmold+airGap+rSpireExt*cos(t); y=rSpireExt*sin(t); label=1001;} border spireMidInt(t=0, 2*pi){x=rmold+airGap+rSpireInt*cos(t); y=rSpireInt*sin(t); label=1002;} //spire top border spireTopExt(t=0, 2*pi){x=rmold+airGap+rSpireExt*cos(t); y=deltaY1+rSpireExt+rSpireExt*sin(t); label=1003;} border spireTopInt(t=0, 2*pi){x=rmold+airGap+rSpireInt*cos(t); y=deltaY1+rSpireExt+rSpireInt*sin(t); label=1004;} //spire bot border spireBotExt(t=0, 2*pi){x=rmold+airGap+rSpireExt*cos(t); y=-deltaY2-rSpireExt-rSpireExt*sin(t); label=1005;} border spireBotInt(t=0, 2*pi){x=rmold+airGap+rSpireInt*cos(t); y=-deltaY2-rSpireExt-rSpireInt*sin(t); label=1006;} // //air (electromagnetic geometry) border topAir(t = 0,rair){x=t;y=hmold/2+hair;label=600;} border rightAir(t = hmold/2+hair,-hmold/2-hair){x=rair;y=t;label=601;} border botAir(t = rair,0){x=t;y=-hmold/2-hair;label=602;} border botLeftAir(t =-hmold/2-hair,-hmold/2 ){x=0;y=t;label=603;} border topLeftAir(t =hmold/2+hair,+hmold/2 ){x=0;y=t;label=604;} real J0 = Irms/spireArea;//Current density //---------------------------------------------------------------------- // Defined edge for BC conditions (A=0) // int[int] potentialBC = [604,600,601,602,603,403]; int[int] potentialBC = [600,601,602]; plot(topMold(-Lcskin)+rightMold(-Lcskin)+botMold(-Lcskin)+leftMold(-Lcskin) +topLeftAir(Lcair)+topAir(-Lcair)+rightAir(-Lcair)+botAir(-Lcair)+botLeftAir(-Lcair) +spireMidExt(Lcspire)+spireMidInt(-Lcspire) +spireTopExt(Lcspire)+spireTopInt(-Lcspire) +spireBotExt(Lcspire)+spireBotInt(+Lcspire), wait=true); // Electromagnetic mesh creation mesh Th1 = buildmesh(topMold(-Lcskin)+rightMold(-Lcskin)+botMold(-Lcskin)+leftMold(-Lcskin) +topLeftAir(Lcair)+topAir(-Lcair)+rightAir(-Lcair)+botAir(-Lcair)+botLeftAir(-Lcair) +spireMidExt(Lcspire)+spireMidInt(-Lcspire) +spireTopExt(Lcspire)+spireTopInt(-Lcspire) +spireBotExt(Lcspire)+spireBotInt(+Lcspire) ); plot(Th1,wait=true); //Set up materials regions int moldRegion = Th1(rmold-1e-03,0).region; int spireMidRegion = Th1(rmold+airGap+1e-03,0).region; int spireTopRegion = Th1(rmold+airGap+1e-03,deltaY1+1e-03).region; int spireBotRegion = Th1(rmold+airGap+1e-03,-deltaY2-1e-03).region; int airRegion = Th1(rair-1e-03,0).region; // Set material properties to regions // Real magnetic permeability func Mu =(muSpireReal)*(region==spireMidRegion)+ (muSpireReal)*(region==spireTopRegion)+(muSpireReal)*(region==spireBotRegion) +(muMoldReal)*(region==moldRegion)+(muAirReal)*(region==airRegion); // relative magnetic permeability // func Mu = (muSpireRelative)*(region==spireMidRegion)+ // (muSpireRelative)*(region==spireTopRegion)+(muSpireRelative)*(region==spireBotRegion) // +(muMoldRelative)*(region==moldRegion)+(muAirRelative)*(region==airRegion); // electrical conductivity func Sigma = (sigmaSpire)*(region==spireMidRegion)+ (sigmaSpire)*(region==spireTopRegion)+(sigmaSpire)*(region==spireBotRegion) +(sigmaMold)*(region==moldRegion)+(sigmaAir)*(region==airRegion); //---------------------------------------------------------------------------------// real f = 52000; //frequency real w =2*pi*f; // angular frequency complex iomega = 1i*w; //Fespace func Pk = P2; fespace Vh(Th1, Pk); Vh A,v; // current source Vh Js = 0 + J0*((Th1(x,y).region == spireMidRegion)+(Th1(x,y).region == spireTopRegion) +(Th1(x,y).region == spireBotRegion)); // zoom in region of interest plot(Js,bb=[[0, -hmold/2-hair/4],[rmold+airGap+2*rSpireExt,hmold/2+hair/4 ]],nbiso=50,fill=true,value=true,wait=true); plot(Js, nbiso=50, fill=true, value=true, cmm="Source current density J [A/m^2]",wait=true); //Macro for vectorial operator macro Curl(A) [dy(A), -dx(A)] // func Nu = 1./Mu; // formulation cartesian coordonates solve maxwellTimeHarmoCartesian(A,v) =int2d(Th1)(Nu *Curl(A)' * Curl(v)) +int2d(Th1)(iomega*Sigma * A * v) - int2d(Th1)(Js * v) + on(potentialBC, A=0) ; //formulation 2 cylindrical coordinate (r,z) (ref Handbook V.rudnev et.al.) // solve maxwellTimeHarmoCylindrical(A,v) // =int2d(Th1)(Nu*(dx(A)*dx(v)+dy(A)*dy(v))*x) // +int2d(Th1)(iomega*Sigma * A * v) // -int2d(Th1)(Js * x * v) // + on(potentialBC, A=0) // ; //formulation 3 (ref Handbook V.rudnev et.al.) // solve maxwellTimeHarmoCylindrical(A,v) // =int2d(Th1)(Nu*(dx(A)*dx(v)+dy(A)*dy(v))) // +int2d(Th1)(iomega*Sigma * A * v) // -int2d(Th1)(Js * v) // + on(potentialBC, A=0) // ; // Real part Az Vh Areal = real(A); // Imaginary part Az Vh Aim = imag(A); // Modulus A Vh Amod; Amod = abs(A); //Magnetic induction Vh Bx, By; Bx = -dy(Amod); By = dx(Amod); // cylindrical Vh B = sqrt(Bx^2 + By^2); //Magnetic field // Vh Hx, Hy; // Hx = Nu * Bx; // Hy = Nu * By; // Vh H = sqrt(Hx^2 + Hy^2); // Induced Current // Vh J = dx(Hy) - dy(Hx); // Eddy current Vh J; J = abs(-iomega*Sigma*A); ////////////////////////////// PLOT ///////////////////////// // plot region numbers cout<< "-----------SIMULATION PARAMETERS-----------"<