三相变压器建模及仿真研究 beta=20;
gama=2*0.016867; C=0.2227;
% 处理初始剩磁 if (begin==1) Hzero=I0; Bzero=B0; if(DI1<0)
sgn=1;
B=(alfa*B0-(gama*I0-pi/2))/(alfa*(atan(beta*(I0+C))+pi/2))*(atan(beta*(I+C))+pi/2)+1/alfa
*(gama*I-pi/2);
L=(alfa*B0-(gama*I0-pi/2))/(alfa*(atan(beta*(I0+C))+pi/2))*beta/(1+(beta*(I+C))*(beta*(I+
C)))+gama/alfa;
else
sgn=-1;
B=(alfa*B0-(gama*I0+pi/2))/(alfa*(atan(beta*(I0-C))-pi/2))*(atan(beta*(I-C))-pi/2)+1/alfa*(
gama*I+pi/2);
L=(alfa*B0-(gama*I0+pi/2))/(alfa*(atan(beta*(I0-C))-pi/2))*beta/(1+(beta*(I-C))*(beta*(I-C )))+gama/alfa; end;
% 结束处理起始磁通 else
% 处理回环时的状态 if(DI0*DI1<0) I0=I; B0=B;
if(DI1>0)
B=(alfa*B0-(gama*I0+pi/2))/(alfa*(atan(beta*(I0-C))-pi/2))*(atan(beta*(I-C))-pi/2)+1
/alfa*(gama*I+pi/2);
L=(alfa*B0-(gama*I0+pi/2))/(alfa*(atan(beta*(I0-C))-pi/2))*beta/(1+(beta*(I-C))*(bet
a*(I-C)))+gama/alfa;
else
B=(alfa*B0-(gama*I0-pi/2))/(alfa*(atan(beta*(I0+C))+pi/2))*(atan(beta*(I+C))+pi/2)+1/alf
a*(gama*I-pi/2);
L=(alfa*B0-(gama*I0-pi/2))/(alfa*(atan(beta*(I0+C))+pi/2))*beta/(1+(beta*(I+C))*(beta*(I
+C)))+gama/alfa;
30
三相变压器建模及仿真研究 end else
if(DI1<0)
B=(alfa*B0-(gama*I0-pi/2))/(alfa*(atan(beta*(I0+C))+pi/2))*(atan(beta*(I+C))+pi/2)+1
/alfa*(gama*I-pi/2);
L=(alfa*B0-(gama*I0-pi/2))/(alfa*(atan(beta*(I0+C))+pi/2))*beta/(1+(beta*(I+C))*(bet
a*(I+C)))+gama/alfa;
else
B=(alfa*B0-(gama*I0+pi/2))/(alfa*(atan(beta*(I0-C))-pi/2))*(atan(beta*(I-C))-pi/2)+1/a
lfa*(gama*I+pi/2);
L=(alfa*B0-(gama*I0+pi/2))/(alfa*(atan(beta*(I0-C))-pi/2))*beta/(1+(beta*(I-C))*(beta
*(I-C)))+gama/alfa;
end; end;
%结束处理回环时状态 end;
%文件名：rk4fixed.m
%功能：四阶定步长龙格－库塔法。
function [T,X]=rk4fixed(Fcn,Tspan,X0,N,K)
% Matlab implementation of a fixed-step RK4 algorithm. h = (Tspan(2)-Tspan(1))/N; halfh = 0.5*h; neqs=size(X0);
X=zeros(neqs(1),N); T=zeros(1,N); X(:,1)=X0; T(1)=Tspan(1); Td = Tspan(1); Xd = X0; for i=2:N,
RK1 = feval(Fcn,Td,Xd,K); Thalf = Td + halfh;
Xtemp = Xd + halfh*RK1;
RK2 = feval(Fcn,Thalf,Xtemp,K); Xtemp = Xd + halfh*RK2;
RK3 = feval(Fcn,Thalf,Xtemp,K);
31
三相变压器建模及仿真研究 Tfull = Td + h;
Xtemp = Xd + h*RK3;
RK4 = feval(Fcn,Tfull,Xtemp,K);
X(:,i) = Xd + h*(RK1+2.0*(RK2+RK3)+RK4)/6; T(i) = Tfull; Xd = X(:,i); Td = T(i); end
X=X';T=T'; End
%文件：trastate.m
% 目的：得到系统微分方程 dIm=inv(L)*(U-r*IM),此文件适合于仿真涌流过程
function dIm = transtatea(t,Im,K) R1=0.7925; L1=0.07336; Um=100000;
dIm=(inv([(L1+K(1)),0,0;0,(L1+K(2)),0;0,0,(L1+K(3))]))*(Um*[sin(314*t);sin(314*t+pi*2/3);sin(314*t+pi*4/3)]-[R1,0,0;0,R1,0;0,0,R1]*Im);
附录B.在Ynd11接线方式下两段反正切函数拟和极限磁滞回环的程序
% 文件名：transa.m
% 目的：仿真电流互感器的暂态过程
%%%%%%%%%% 初始化%%%%%%%%%%%% clear; %%清除变量 %磁滞回线起始点 B0a=0; I0a=0; B0b=0; I0b=0; B0c=0; I0c=0; B0d=0; I0d=0;
%变量：磁通 Ba=0;
32
三相变压器建模及仿真研究 Bb=0; Bc=0; Bd=0;
K=[1242.496;1242.496;1242.496;0]; %变量：励磁电感 options=[];%RK法参数设定，此处为默认值 deltaT=0.0001; %采样时间间隔 Tmax=0.4; %仿真时间
tspan = [0,deltaT]; %RK法时间宽度 yzero = [0;0;0;0]; %RK法初始值
%励磁电流变化率，象征磁滞回线的改换方向的点 DI0a=0; DI0b=0; DI0c=0; DI0d=0;
%%%%%%%%%%%%%%%%仿真开始，是个循环过程%%%%%%%%%%5 for i=1:Tmax/deltaT
% 四阶RK解微分方程,注意此为定步长RK [t,num_y] =rk4fixed(@transtatea,tspan,yzero,2,K); % 更新参数，为下次RK做准备
tspan=[t(length(t),1),t(length(t),1)+deltaT]; tt(i)=[t(length(t),1)]; yzero=num_y(2,:); yzero=yzero';
Im(i,:)=num_y(2,:); %取得励磁电流变化率
DI1a=diff(num_y(:,1))./diff(t); DI1b=diff(num_y(:,2))./diff(t); DI1c=diff(num_y(:,3))./diff(t); DI1d=diff(num_y(:,4))./diff(t); %取得当前状态的励磁电感值
[Ka,B0a,I0a,Ba]=getLm(Ba,Im(i,1),B0a,I0a,DI1a,DI0a,i); [Kb,B0b,I0b,Bb]=getLm(Bb,Im(i,2),B0b,I0b,DI1b,DI0b,i); [Kc,B0c,I0c,Bc]=getLm(Bc,Im(i,3),B0c,I0c,DI1c,DI0c,i); [Kd,B0d,I0d,Bd]=getLm(Bd,Im(i,4),B0d,I0d,DI1d,DI0d,i); %为计算励磁电感更新参数 DI0a=DI1a; DI0b=DI1b; DI0c=DI1c; DI0d=DI1d;
%输出中间变量，电感和磁通 K=[Ka;Kb;Kc;0]; end;
33
三相变压器建模及仿真研究
%文件名：getLm.m
功能:按磁滞回线计算动态电感 可以使用§1中的getLm.m程序。
%文件名：rk4fixed.m
功能：四阶定步长龙格－库塔法。 可以使用§1中的rk4fixed.m程序。
%文件：trastate.m
% 目的：得到系统微分方程 dIm=inv(L)*(U-r*IM),此文件适合于仿真涌流过程
function dIm = transtatea(t,Im,K) R1=0.7925; L1=0.07336; Um=100000;
dIm=(inv([(L1+K(1)),0,0,L1;0,(L1+K(2)),0,L1;0,0,(L1+K(3)),L1;L1,L1,L1,6*L1]))*(Um*[sin(314*t);sin(314*t+pi*2/3);sin(314*t+pi*4/3);0]-[R1,0,0,R1;0,R1,0,R1;0,0,R1,R1;R1,R1,R1,6*R1]*Im);
附录C.在Yny0接线方式下两段反正切函数拟和极限磁滞回环的程序
%文件名：transa.m
目的：仿真电流互感器的暂态过程 可以使用§1中的transa.m程序。
%文件名：getLm.m
功能:按磁滞回线计算动态电感 可以使用§1中的getLm.m程序。
34