Matlab Pdetool Lezione 5  

 Corren'  Indo*e    

Formulazione  nel  caso  di  corren'  indo*e  

EJB

JJH

BE

σ==⋅∇

+=×∇∂∂−=×∇

eddy

seddy

0

tSonda  con  corrente  imposta  Js  

Pia*o  condu*ore  sede  di  corren'  indo*e  Jeddy  

Formulazione  

gaugedicondizione00

0

=⇒=⋅∇

+∇−∂∂−=×∇

∇−∂∂−=

⇓∂∂×−∇=×∇

×∇=⇒=⋅∇

ϕ

ϕσσ

ϕ

A

JAH

AE

AE

ABB

st

t

t

Il potenziale vettore A è univocamente determinato dalla gauge di Coulomb

Equazione  risolvente  

∇× 1µ0

∇×A = −σ ∂A∂t

+ Js

∇×∇×A = −∇2A +∇(∇⋅A)A = AMe

jωt

⇓

−∇2A = − jωµ0σA + µ0Js

Caso  Assial-­‐simmetrico  

ΦΦΦΦ =⇒= iAiJ ˆAˆJs

Φ=+∂∂

∂∂−

∂∂

∂∂− J1)1()1( 00 µ

ρσωµ

ρρρρPhijPhi

zzPhi

Φ= APhi ρx

y

Pdetool  

Configurazione  Geometrica  

Sonda  

Pia*o  condu*ore  

Definizione  Parametri  Sonda-­‐Pia*o  

Sonda

Parametri geometrici Rmin Rmax Zmin Zmax Nspire Corrente Lift off sonda-piatto Conducibilità

Rmin_disk=[0]; Rmax_disk=[0.04]; Zmin_disk=[-0.0015/2]; Zmax_disk=[0.0015/2]; sigma_disk=1/(0.027e-6); %%%%%%%%%%%%%%%%%%%%%%%%%%%%% liftoff=0.5e-3; Rmin_transmitter=[0.005]; Rmax_transmitter=[0.005+1e-3]; Zmin_transmitter=[liftoff+0.0015/2]; Zmax_transmitter=[liftoff+4.5e-3+0.0015/2]; sigma_transmitter=0; I_transmitter_pp=1;%corrente unitaria I_transmitter=[I_transmitter_pp/(2*sqrt(2))]; % corrente rms N_transmitter=[16]; areat=(Rmax_transmitter-Rmin_transmitter)*(Zmax_transmitter-Zmin_transmitter); Jt=(I_transmitter*N_transmitter)/areat; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% Rmin_grid=[0]; Rmax_grid=[5e-2]; Zmin_grid=[-5e-2]; Zmax_grid=[5e-2]; sigma_grid=0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Creazione  Geometria  e  Mesh   [pde_fig,ax]=pdeinit; pdetool('appl_cb',1); set(ax,'DataAspectRatio',[1 1 1]); set(ax,'PlotBoxAspectRatio',[1 1 1]); set(ax,'XLim',[Rmin_grid-0.5*(Rmax_grid-Rmin_grid) Rmax_grid+0.5*(Rmax_grid-Rmin_grid)]); set(ax,'YLim',[Zmin_grid-0.5*(Zmax_grid-Zmin_grid) Zmax_grid+0.5*(Zmax_grid-Zmin_grid)]); set(ax,'XTickMode','auto'); set(ax,'YTickMode','auto'); pderect([0 Rmax_grid Zmin_grid Zmax_grid],'R1'); pderect([Rmin_disk Rmax_disk Zmin_disk Zmax_disk],'R2') pderect([Rmin_transmitter Rmax_transmitter Zmin_transmitter Zmax_transmitter], 'R3') set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','R1+R2+R3'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Mesh generation: 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'); [p,e,t]=refinemesh(dl,p,e,t,'regular'); pdeplot(p,e,t)

0 0.01 0.02 0.03 0.04 0.05-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05

Condizioni  al  contorno  %%%condizioni al contorno pdetool('changemode',0) pdesetbd(11,... 'dir',... 1,... '1',... '0') pdesetbd(10,... 'dir',... 1,... '1',... '0') pdesetbd(9,... 'dir',... 1,... '1',... '0') pdesetbd(6,... 'dir',... 1,... '1',... '0') pdesetbd(2,... 'dir',... 1,... '1',... '0') pdesetbd(1,... 'dir',... 1,... '1',... '0') % PDE coefficients: h=findobj(get(pde_fig,'Children'),'flat','Tag','PDEBoundMenu'); bl=get(findobj(get(h,'Children'),'flat','Tag','PDEBoundMode'),'UserData');

Definizione  parametri  regioni  box=find(t(4,:)==1); %indici delle regioni sonda=find(t(4,:)==2); piatto=find(t(4,:)==3); Js=zeros(size(t,2),1); %densità di corrente imposta Js(sonda)=Jt; sigma(box)=0; sigma(piatto)=sigma_disk; sigma(sonda)=0; figure pdeplot(p,e,t(:,box)) %%%verifica regioni title('mesh box’); axis equal figure pdeplot(p,e,t(:,piatto)) title('mesh piatto'); axis equal figure pdeplot(p,e,t(:,sonda)) title('mesh sonda'); axis equal

Definendo le regioni del dominio totale è possibile evitare l’utilizzo dei cicli for sui triangoli della mesh

Equazione  Risolvente  

Phi=assempde(bl,p,e,t,c,a,f) %assembla e risolve  

-div(c grad(Phi))+a Phi=f

Φ=+∂∂

∂∂−

∂∂

∂∂− J1)1()1( 00 µ

ρσωµ

ρρρρPhijPhi

zzPhi

c=1/r a=jωµ0σ/r f=µ0JF

Soluzione    

rc=pdeintrp(p,t,p(1,:)');% r area=abs(pdetrg(p,t));%area di ogni triangolo mu0=4*pi*1e-7; f=20*1e3; % freq aa=j*2*pi.*f.*mu0.*sigma./rc; ff=mu0.*Js; Phi = assempde(bl,p,e,t,1./rc,aa,ff','0'); figure %plot linee di flusso subplot(1,2,1) pdegplot(dl) hold on pdecont(p,t,real(Phi)'*2*pi,10) title('Real') subplot(1,2,2) pdegplot(dl) hold on pdecont(p,t,imag(Phi)'*2*pi,10) title('Imag‘)

Φ=+∂∂

∂∂−

∂∂

∂∂− J1)1()1( 00 µ

ρσωµ

ρρρρPhijPhi

zzPhi

0 0.02 0.04 0.06-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05Real

0 0.02 0.04 0.06-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05Imag

0 0.02 0.04 0.06-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05Real

0 0.02 0.04 0.06-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05Imag

Linee  di  flusso  

Conducibilità piastra=0 Conducibilità piastra alluminio=37 MS/m

Inserendo il piatto conduttore il flusso ha anche una componente immaginaria; la tensione indotta presenta uno sfasamento che dipende dalla frequenza di eccitazione e dalla conducibilità del piatto.

Esercitazione  Laboratorio  

Configurazione  geometrica  

Sonda  trasmi*ente  

Sonda  ricevente  

Pia*o  condu*ore  

Geometria  Rmin_disk=[0]; Rmax_disk=[0.04]; Zmin_disk=[-0.0015/2]; Zmax_disk=[0.0015/2]; sigma_disk=1/(0.027e-6); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% liftoff=0.5e-3; Rmin_transmitter=[0.005]; Rmax_transmitter=[0.005+1e-3]; Zmin_transmitter=[liftoff+0.0015/2]; Zmax_transmitter=[liftoff+4.5e-3+0.0015/2]; sigma_transmitter=0; I_transmitter_pp=1;%corrente unitaria I_transmitter=[I_transmitter_pp/(2*sqrt(2))]; % corrente rms N_transmitter=[16]; areat=(Rmax_transmitter-Rmin_transmitter)*(Zmax_transmitter-Zmin_transmitter); Jt=(I_transmitter*N_transmitter)/areat; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Rmin_receiver=[0.005]; Rmax_receiver=[0.005+1e-3]; Zmin_receiver=[-0.0015/2-liftoff-4.5e-3]; Zmax_receiver=[-0.0015/2-liftoff]; sigma_receiver=0; N_receiver=[16]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Rmin_grid=[0]; Rmax_grid=[5e-2]; Zmin_grid=[-5e-2]; Zmax_grid=[5e-2]; sigma_grid=0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  

 

Sonda

Parametri geometrici Rmin Rmax Zmin Zmax Nspire Corrente Lift off sonda-piatto Conducibilità

Creazione  geometria  e  mesh  

0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05

[pde_fig,ax]=pdeinit; pdetool('appl_cb',1); set(ax,'DataAspectRatio',[1 1 1]); set(ax,'PlotBoxAspectRatio',[1 1 1]); set(ax,'XLim',[Rmin_grid-0.5*(Rmax_grid-Rmin_grid) Rmax_grid+0.5*(Rmax_grid-Rmin_grid)]); set(ax,'YLim',[Zmin_grid-0.5*(Zmax_grid-Zmin_grid) Zmax_grid+0.5*(Zmax_grid-Zmin_grid)]); set(ax,'XTickMode','auto'); set(ax,'YTickMode','auto'); pderect([0 Rmax_grid Zmin_grid Zmax_grid],'R1'); pderect([Rmin_disk Rmax_disk Zmin_disk Zmax_disk],'R2') pderect([Rmin_transmitter Rmax_transmitter Zmin_transmitter Zmax_transmitter], 'R3') pderect([Rmin_receiver Rmax_receiver Zmin_receiver Zmax_receiver], 'R4') set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','R1+R2+R3+R4'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Mesh generation: 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'); [p,e,t]=refinemesh(dl,p,e,t,'regular'); pdeplot(p,e,t)

Condizioni  al  contorno  % Boundary conditions: pdetool('changemode',0) pdesetbd(9,... 'dir',... 1,... '1',... '0') pdesetbd(8,... 'dir',... 1,... '1',... '0') pdesetbd(7,... 'dir',... 1,... '1',... '0') pdesetbd(6,... 'dir',... 1,... '1',... '0') pdesetbd(2,... 'dir',... 1,... '1',... '0') pdesetbd(1,... 'dir',... 1,... '1',... '0') % PDE coefficients: h=findobj(get(pde_fig,'Children'),'flat','Tag','PDEBoundMenu'); bl=get(findobj(get(h,'Children'),'flat','Tag','PDEBoundMode'),'UserData');

Definizione  parametri  regioni  

box=find(t(4,:)==1); %indici delle regioni piatto=find(t(4,:)==2); ric=find(t(4,:)==3); trasm=find(t(4,:)==4); Js=zeros(size(t,2),1); Js(trasm)=Jt; sigma(box)=0; sigma(piatto)=sigma_disk; sigma(trasm)=0; sigma(ric)=0; rc=pdeintrp(p,t,p(1,:)');%baricentro triangoli area=abs(pdetrg(p,t));%area di ogni triangolo mu0=4*pi*1e-7;

Definendo le regioni del dominio totale è possibile evitare l’utilizzo dei cicli for sui triangoli della mesh

Calcolo  della  tensione  sulla  bobina  ricevente  

f=[0.5 1 2 5 10 20]*1e3; % freq for kf=1:length(f) %calcolo su più frequenze aa=j*2*pi.*f(kf).*mu0.*sigma./rc; ff=mu0.*Js; Phi = assempde(bl,p,e,t,1./rc,aa,ff','0'); omega(kf)=2*pi*f(kf); PhiT=pdeintrp(p,t,Phi); ee=2*pi*(Jt.*PhiT(ric).*area(ric)); energiamutua=sum(ee); Vreceiver(kf)=energiamutua*j*omega(kf); figure subplot(1,2,1) pdegplot(dl) hold on pdecont(p,t,real(Phi)'*2*pi,10) title('Real') subplot(1,2,2) pdegplot(dl) hold on pdecont(p,t,imag(Phi)'*2*pi,10) title('Imag') end

12121

22121

21

2

1

iMjV

dVJAii

MV

ω=

= ∫

Tensione indotta sulla bobina ricevente (2)

Plot  flusso  campo  magne'co  

0 0.05 0.1-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05Real

0 0.05 0.1-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05Imag

0 0.05 0.1-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05Real

0 0.05 0.1-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05Imag

Conducibilità piastra=0 Conducibilità piastra alluminio=37 MS/m

Inserendo il piatto conduttore il flusso ha anche una componente immaginaria; la tensione indotta presenta uno sfasamento che dipende dalla frequenza di eccitazione e dalla conducibilità del piatto.

Flusso  a  diverse  frequenze  Sigma=58 MS/m conducibilità rame

1 KHz 10 KHz

100 KHz

All’aumentare della frequenza diminuisce lo spessore di penetrazione d e il campo viene schermato quasi completamente dalla piastra conduttrice

0

2σωµ

=δ

Calcolo  conducibilità   Effettuando misure sperimentali è possibile risalire alla conducibilità dell’oggetto in esame, valutando la variazione della fase della tensione indotta a varie frequenze e inserendo nel calcolo numerico una conducibilità di tentativo vicina a quella del tipo di materiale considerato.

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2

x 104

-150

-100

-50

0

50

f (Hz)

V ph

ase

(deg

)

Numerico sigma=37 MS/mNumerico sigma=50 MS/mNumerico sigma=20 MS/mSperimentale

figure, plot(f,angle(Vreceiver)*180/pi,f,angle(Vreceiver2)*180/pi,f,angle(Vreceiver3)*180/pi,f,exp,'or') xlabel('f (Hz)') ylabel('V phase (deg)') legend('Numerico sigma=37 MS/m','Numerico sigma=50 MS/m','Numerico sigma=20 MS/m','Sperimentale') axis([0 22000 -180 80])  

NB: le vreceiver vanno calcolate singolarmente al variare di sigma_disk