计算磁滞回线Jiles-Atherton模型matlab程序；可用于各向异性磁性材料的磁滞回线计算.zip

clear all clc fprintf('\n\nDemonstration of Jiles-Atherton model - three hysteresis loops.'); Ms=1.6e6; a=1100; alpha=1.6e-3; k=400; c=0.2; % Parameters of Jiles-Atherton model fprintf('Calculation for parameters: \na=%f(A/m), k=%f(A/m), c=%f, Ms=%e(A/m), alpha=%e \n\n',a,k,c,Ms,alpha); H=[0:0.01:1 1:-0.01:-1 -1:0.01:1]'; % magnetizing field H - column vector H=[H.*1000 H.*3000 H.*6000]; fprintf('Calculations started. Expected time of calculations: below 2 minutes.\n'); Bmodel = func_ja(a,k,c,Ms,alpha,H); fprintf('Calculations finished.\n\n'); plot(H, Bmodel,'-o'); xlabel('H (A/m)');ylabel('B (T)');grid; function BsimT = func_ja(a,k,c,Ms,alpha,HmeasT) % INPUT: % a - quantifies domain density, A/m (scalar) % k - quantifies average energy required to break pinning side, A/m (scalar) % c - magnetization reversibility, 0..1 (scalar) % Ms - saturation magnetization, A/m (scalar) % alpha - Bloch coefficient (scalar) % H - magnetizing field, A/m (column vector) % M0 - initial value of magnetization for ODE, A/m (scalar) % OUTPUT: % Hw - set of output values of magnetizing field H, A/m (vector) % Bw - set of output values of flux density B, T (vector) if ((a<0) || (k<0) || (c<0) || (c>1) || (alpha<0)) fprintf('\n\n*** ERROR in JAn_loops: unphysical parameters.\n\n'); BsimT=zeros(size(HmeasT)); return end % check if parameters are physical BsimT=zeros(size(HmeasT)); % Prepare output variable InitialCurve=0; if sum(HmeasT(1,:))==0 InitialCurve=1; end M0=0; % sample demagnetized at the beginning for lT=1:size(HmeasT,2) H=HmeasT(:,lT); if InitialCurve==0 H=[0; H]; end % check if there is initial magnetization curve in H [Hw,Bw] = func_jaloop(a,k,c,Ms,alpha,H,M0); if InitialCurve==0 Hw(1)=[]; Bw(1)=[]; H(1)=[]; Bw=Bw+linspace(0,Bw(1)-Bw(numel(Bw)),numel(Bw))'; % Close hysteresis loop Bw=Bw-(max(Bw)+min(Bw))./2 ; % Symmetrize Bw end BsimT(:,lT)=Bw; % Remember the result end end function [Hw,Bw] = func_jaloop(a,k,c,Ms,alpha,H,M0) mi0=4.*pi.*1e-7; ip=1; ik=2; Hw=zeros(size(H)); Hw(1)=H(1); Mw=zeros(size(H)); Mw(1)=M0; while ik<numel(H) if H(ik)>H(ip) if H(ik+1)>=H(ik) ik=ik+1; else [Hi, Mi]=func_jasolver(a,k,c,Ms,alpha,H(ip),H(ik),M0); Hi(isnan(Hi))=1; Mi(isnan(Mi))=1; if numel(Hi)>2 Mw(ip+1:ik)=interp1(Hi,Mi,H(ip+1:ik),'pchip'); else Mw(ip+1)=Mi(end); end Hw(ip+1:ik)=H(ip+1:ik); M0=Mi(end); ip=ik; ik=ip+1; end end if H(ik)<H(ip) if H(ik+1)<=H(ik) ik=ik+1; else [Hi, Mi]=func_jasolver(a,k,c,Ms,alpha,H(ip),H(ik),M0); Hi(isnan(Hi))=1; Mi(isnan(Mi))=1; if numel(Hi)>2 Mw(ip+1:ik)=interp1(Hi,Mi,H(ip+1:ik),'pchip'); else Mw(ip+1)=Mi(end); end Hw(ip+1:ik)=H(ip+1:ik); M0=Mi(end); ip=ik; ik=ip+1; end end if H(ik)==H(ip) ik=ik+1; end end [Hi, Mi]=func_jasolver(a,k,c,Ms,alpha,H(ip),H(ik),M0); Hi(isnan(Hi))=1; Mi(isnan(Mi))=1; if numel(Hi)>2 Mw(ip+1:ik)=interp1(Hi,Mi,H(ip+1:ik),'pchip'); else Mw(ip+1)=Mi(end); end Hw(ip+1:ik)=H(ip+1:ik); M0=Mi(end); Bw=Mw.*mi0+Hw.*mi0; end function [H,M] = func_jasolver(a,k,c,Ms,alpha,Hstart,Hend,M0) options=odeset('RelTol',1e-4,'AbsTol',1e-6,'MaxStep',abs(Hend-Hstart)./10,'InitialStep',(Hend-Hstart)./10); dMdH_=@(H,M) dMdH(a,k,c,Ms,alpha,M,H,Hstart,Hend); [H,M] = ode45(dMdH_,[Hstart Hend],M0,options); end