Matlab Pdetool

Lezione 3Magnetostatica

Ing. Flavio Calvano

• Modello magnetostatica;• Utilizzo del pdetool in modalità grafica;• Esercitazione di laboratorioCalcolo auto-induttanza;Calcolo mutua-induttanza;

Confronto dati numerici-sperimentali.

Magnetostatica

Equazioni integrali magnetostatica

Formulazione differenziale

coulombdigauge0

''

0

S0S0

=⋅∇↓

×∇=×∇=⇒ϕ∇+=µ=×∇×∇⇒µ=×∇

×∇=⇒=⋅∇

A

AABAAJAJB

ABB

Il potenziale vettore A è univocamente determinato dalla gauge di coulomb

Caso 2d assial-simmetrico

( )

2

0 0 02

0

22 2 2 2

02

1 1 1 0

2 2ˆ ˆ ˆ

ˆ ˆ

z

con

A A A AA A A z

J AA

A J

ρ ρ ϕ ρρ ϕ

ϕϕ ϕ

µ µ µ

µ

ρ ϕρ ρ ϕ ρ ρ ϕ

ϕ ϕ

µρ

Φ Φ

∇× ∇× = ∇×∇× = ∇∇ ⋅ − ∇ ∇ ⋅ =

∇ = −

∂ ∂ ∇ = ∆ − − + ∆ − − + ∆ ∂ ∂

= ⇒ =

−∆ + =

A A A A A

A J

A

J A

J

Formulazione

x

y

Equazione risolvente

( )

( )

2

2

0 02 2

2

2 2 2

0

1 ( )

1 1( )

1 1 1 1 ( )

1 1( ) ( )

A A Az

A AA J A A J

z zA AA A AA A A

A

Jz z

ϕ

ϕ ϕϕ ϕ ϕ

ϕ ϕ

ϕ

ρρ ρ ρ

µ ρ ρ µρ ρ ρ ρ ρ ρ

ρρ ρρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ

ρ

µρ ρ ρ ρ

Φ Φ

Φ Φ

Φ Φ ΦΦ Φ Φ

Φ

∂ ∂ ∂∆ = +

∂ ∂ ∂

∂ ∂ ∂ ∂−∆ + = ⇒ − + − =

∂ ∂ ∂ ∂

∂ ∂ ∂∂ ∂ ∂ ∂ ∂ ∂− = − − = + − − = − +

∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂

⇓Ψ =

∂ ∂ ∂ ∂− Ψ − Ψ =

∂ ∂ ∂ ∂

Pdetool

Creazione geometria

xmax=2.2e-02; dr=2.5e-03xmin= xmax-dr;dh=4.0e-02;pderect([xmin xmax -dh dh],'R1');pderect([0 10*xmax 10*(-dh) 10*dh],'R2');

Condizioni al contornopdesetbd(8,...'dir',...1,...'1',...'0')pdesetbd(5,...'dir',...1,...'1',...'0')pdesetbd(4,...'dir',...1,...'1',...'0')pdesetbd(3,...'dir',...1,...'1',...'0')

Mesh

R2=box aria ->t(4,:)=1R1=solenoide ->t(4,:)=2

Dall’indice della quarta riga della matrice t è possibile riconoscere i domini e imporre il termine noto J

0 0.05 0.1 0.15 0.2 0.25-0.25

-0.2

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

0.25

set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','R2+R1');gd=get(findobj(get(pde_fig,'Children'),'flat','Tag','PDEMeshMenu'),'UserData');dl=decsg(gd);[p,e,t]=initmesh(dl,'Hmax',5e-1,'init','off');[p,e,t]=refinemesh(dl,p,e,t,'regular');h=findobj(get(pde_fig,'Children'),'flat','Tag','PDEBoundMenu');bl=get(findobj(get(h,'Children'),'flat','Tag','PDEBoundMode'),'UserData');

Grafico Regioni separate

box=find(t(4,:)==1);solenoide=find(t(4,:)==2);pdeplot(p,e,t(:,solenoide))

pdeplot(p,e,t(:,box))

0 0.02 0.04 0.06 0.08 0.1 0.12

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

0 0.01 0.02 0.03 0.04 0.05 0.06

-0.08

-0.06

-0.04

-0.02

0

0.02

0.04

0.06

0.08

0.1

Definizione del forzamento

Nspire=500;dr=2.0e-03;dh=8.0e-02;I=1;J0=Nspire*I/(dh*dr); %nell’area occupata dalla bobina

Calcolo Soluzione

%coordinata ρ dei baricentri dei triangoli

rc=pdeintrp(p,t,p(1,:)');

Nt=size(t,2); Numero dei triangoli

J= zeros(Nt,1);

J(solenoide)=J0;

J definita solo in R2 (t(4,:)=2)corrispondente all’area della bobina, nulla all’esterno

Phi = assempde(bl,p,e,t,c,'0',f);

-div∙c(grad(Phi))+aPhi=f

Con:mu0=4.e-7*pi;c= 1./(rc*mu0)f= J;

Plot linee di flusso

∫∫∫γ

Φπ=⋅=⋅×∇=⋅= rA2dlndSdSSS

AAnBΦ pdecont(p,t,Phi*2*pi,10)

S

Esempio 2

pderect([xmin xmax ymax-dh ymax],'R1');pderect([xmin xmax ymax-2*dh-h ymax-1*dh-h],'R2');pderect([xmin xmax ymax-3*dh-2*h ymax-2*dh-2*h],'R3');pderect([xmin xmax ymax-4*dh-3*h ymax-3*dh-3*h],'R4');pderect([xmin xmax ymax-5*dh-4*h ymax-4*dh-4*h],'R5');pderect([xmin xmax ymax-6*dh-5*h ymax-5*dh-5*h],'R6');pderect([xmin xmax ymax-7*dh-6*h ymax-6*dh-6*h],'R7');pderect([xmin xmax ymax-8*dh-7*h ymax-7*dh-7*h],'R8');pderect([xmin xmax ymax-9*dh-8*h ymax-8*dh-8*h],'R9');pderect([xmin xmax ymax-10*dh-9*h ymax-9*dh-9*h],'R10');pderect([0 25*xmax -15*xmax 15*xmax],'R11');

xmax=1.e-2;dr=1.e-3;xmin=xmax-dr;dh=1.e-3;ymax=2.e-2;h=3.e-3;

Condizioni al contornopredettolo('changemode',0)pdesetbd(5,...'dir',...1,...'1',...'0')pdesetbd(6,...'dir',...1,...'1',...'0')pdesetbd(3,...'dir',...1,...'1',...'0')pdesetbd(2,...'dir',...1,...'1',...'0')

Assemblaggio termine noto

%coordinata ρ dei baricentri dei triangoli

rc=pdeintrp(p,t,p(1,:)');

Nt=size(t,2); %Numero dei triangoli

J= zeros(1,Nt);

spire=find(t(4,:)>1);J0=I/(dh*dr);J(spire)=J0;

J uguale a zero in R11 (t(4,:)=1)corrispondente al box aria , mentre è diversa da zero negli avvolgimenti che hanno tutte la stessa corrente

Phi= assempde(bl,p,e,t,c,'0',f);

Con:µυ0=4.ε−7∗πι;χ= 1./(ρχ∗µυ0);φ= ϑ;

Flusso

Linee di flusso

figurepdegplot(dl)hold onpdecont(p,t,Phi,20)axis equal

Esercitazione laboratorio

Definizione Geometria

dh=4.4e-02;dr=1.0e-03;xmin=0.03-dr/2;xmax=0.03+dr/2;xmax2=0.032+dr/2;xmin2=0.032-dr/2;

pderect([xmin xmax -dh/2 dh/2],'R1');pderect([xmin2 xmax2 -dh/2 dh/2],'R2');

pderect([0 20*xmax 10*(-xmax) 10*xmax],'R3');

set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String',‘R3+R1+R2');

t(4,:)=1 t(4,:)=2 t(4,:)=3

Nspire=500;I1=0.0023; %dati laboratorioI2=0.0014;J01=Nspire*I1/(dh*dr);J02=Nspire*I2/(dh*dr);

mu0=4.e-7*pi;

J1=zeros(1,size(t,2));J2=zeros(1,size(t,2));

%coordinata ρ del baricentro dei triangolirc=pdeintrp(p,t,p(1,:)');area=abs(pdetrg(p,t));%area di ogni triangolo

spira1=find( t(4,:)==2);spira2=find( t(4,:)==3);

J1(spira1)=J01;J2(spira2)=J02;

Phi1 = assempde(bl,p,e,t,1/(rc.*mu0),’0’,J1,'0');Phi2 = assempde(bl,p,e,t,1/(rc.*mu0),’0’,J2,'0');

Definizione forzamento

Calcolo induttanze

Calcolo AUTO-INDUTTANZA

spira1=find(t(4,:)==2);

Phi1 = assempde(bl,p,e,t,1./(rc.*mu0),’0’,J1,'0');PhiT1=pdeintrp(p,t,Phi1);

area=abs(pdetrg(p,t));%area di ogni triangolo

Energia11=2*pi*sum(J1(spira1)’.*PhiT1(spira1).*area(spira1));

L1=Energia11/I1^2

Calcolo MUTUA-INDUTTANZA

( )2 2

21 2 1 2 2 1 2 2 21 2 1 2

1 2

V S

M dV J A d dzi i i i

π ρ ρ= ⋅ =∫ ∫J A

Phi1 = assempde(bl,p,e,t,1/(rc.*mu0),’0’,J1,'0');PhiT1=pdeintrp(p,t,Phi1);

area=abs(pdetrg(p,t));%area di ogni triangolo

spira2=find(t(4,:)==3);

Energia12=2*pi*sum(J2(spira2).*PhiT1(spira2).*area(spira2));

M12=Energia12/(I1*I2)

Confronto dati numerici-sperimentali

M21energia = 11.8 mHM21misurata=14.8 mH

L11energia = 12 .2 mHL11misurata=13.5 mH

L22energia= 13.4 mHL22misurata= 16 .0 mH

M12energia = 11.8 mHM12misurata=14.6 mH