facat
/
potential
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

基本完成逐次镜像法 

Signed-off-by: dmy@lab <dmy@lab.com>
This commit is contained in:
dmy@lab 2015-10-08 14:01:55 +08:00
parent 3a73dfd84e
commit 2ac1e54f69
2 changed files with 163 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

114
GradualMirror.m Normal file
View File

@ -0,0 +1,114 @@
clc
clear
%% 逐次镜像法
% [1]. 付宾兰, 邵.A., 高压输电线路分裂导线表面和周围电场的计算. 电网技术, 1984(Z1): 第83-91页.
% [2]. 韦钢与李海峰, 同杆并架多回线下方的电场强度和感应电压. 中国电力, 1999(03): 第39-42页.
% 采用文献的结果： 蒋兴良, 胡.A.舒.A., Analysis of Conductors＇ Surface Electric Field of UHVDC Transmission Lines Based on Optimized Charge Simulation Method. 高电压技术, 2008(12): p. 2547-2551.
%%
%设置几个参数
semi_lineDistance=450;%分裂间距
semi_lineCount=6;%分裂数
ConductorX=[-11000,11000];%导线间距
ConductorY=[22000,22000,];%导线距地高度
subconductorR=16.8;%子导线半径
phaseN=2;%相数，单回三相
%%
eslong=8.854187817*10^-12*1000;
%设置电压
% Volt_=[1100/sqrt(3);1100/sqrt(3)*exp(1j*4/3*pi);1100/sqrt(3)*exp(1j*2/3*pi);];
Volt_=[800;-800;];
Volt=[];
for vLoop=1:phaseN
Volt=[Volt;Volt_(vLoop)*ones(semi_lineCount,1);];
end
%按分裂数和分裂导线间距布置单相线路导线
%用极坐标
arc=2*pi/semi_lineCount;
%子导线中心到导线中心的距离
R=semi_lineDistance/2/sin(arc/2);
%计算导线互电位和自电位系数
subconductorPos=[];
for phaseLoop=1:phaseN
for sC=1:semi_lineCount
subconductorPos=[subconductorPos;ConductorX(phaseLoop)+1j*ConductorY(phaseLoop)+exp(1j*((sC-1)*arc+arc/2))*R];%移动到子导线中心
%同时计算匹配点的位置
% matchPos=[matchPos;simulationChargePos/r1*subconductorR+ConductorX(phaseLoop)+1j*ConductorY(phaseLoop)+exp(1j*((sC-1)*arc+arc/2))*R];
end
end
mirrorSubconductorPos=conj(subconductorPos);%获得子导线镜像
H=diag(imag(subconductorPos));
r=eye(length(H))*subconductorR;
matSubconductor=repmat(subconductorPos,1,length(subconductorPos));
conductor2conductorDistance=abs(matSubconductor-conj(matSubconductor'));
% conductor2conductorDistance=conductor2conductorDistance-diag(diag(conductor2conductorDistance));
conductor2MirrorDistance=abs(matSubconductor-repmat(conj(mirrorSubconductorPos'),length(subconductorPos),1));
% conductor2MirrorDistance=conductor2MirrorDistance-diag(diag(conductor2MirrorDistance));
P1=1/2/pi/eslong*log(2*H./r);
P1(isnan(P1))=0;
P2=1/2/pi/eslong*log(conductor2MirrorDistance./conductor2conductorDistance);
P2(isinf(P2))=0;
Pij=P1+P2;
%求电荷
QRI=Pij\Volt;
% 计算镜像电荷
%只计算同极性和同极性对地镜像的镜像
innerMirrorPos=[];
innerMirrorQ=[];%内部镜像的电荷
for phaseLoop=1:phaseN
for sCOuter=1:semi_lineCount
for sCInner=1:semi_lineCount
if sCInner==sCOuter
innerMirrorPos=[innerMirrorPos;subconductorPos((phaseLoop-1)*semi_lineCount+sCOuter)];%先预留一个位置
innerMirrorQ=[innerMirrorQ;sum(QRI(1+(phaseLoop-1)*semi_lineCount:phaseLoop*semi_lineCount))];
continue
end
innerMirrorPos=[innerMirrorPos;subconductorPos(sCOuter+(phaseLoop-1)*semi_lineCount)+subconductorR^2/abs(subconductorPos(sCOuter+(phaseLoop-1)*semi_lineCount)-subconductorPos(sCInner+(phaseLoop-1)*semi_lineCount))*(subconductorPos(sCInner+(phaseLoop-1)*semi_lineCount)-subconductorPos(sCOuter+(phaseLoop-1)*semi_lineCount))./(abs(subconductorPos(sCInner+(phaseLoop-1)*semi_lineCount)-subconductorPos(sCOuter+(phaseLoop-1)*semi_lineCount)))];
innerMirrorQ=[innerMirrorQ;-QRI(sCInner+(phaseLoop-1)*semi_lineCount)];
end
end
end
%选检验导线上一个角度
vrfRelA=linspace(0,2*pi,200)';%vrf=verify
%计算检验点相对于子导线的位置
vrfRelPos=exp(1j*vrfRelA)*subconductorR;
%移动坐标，使验证的子导线中心和实际子导线中心重合。
vrfPos=[];
for phaseLoop=1:phaseN
for sC=1:semi_lineCount
vrfPos=[vrfPos;exp(1j*((sC-1)*arc+arc/2))*R+ConductorX(phaseLoop)+1j*ConductorY(phaseLoop)+vrfRelPos];
end
end
%计算这一点的电位系数
matVrfPos=repmat(vrfPos,1,length(innerMirrorPos));
vrf2ConductorDistance=abs(matVrfPos-repmat(conj(innerMirrorPos'),length(vrfPos),1));
vrf2MirrorDistance=abs(matVrfPos-repmat(conj(conj(innerMirrorPos')),length(vrfPos),1));
Pij=1/2/pi/eslong*log(vrf2MirrorDistance./vrf2ConductorDistance);
%计算电压
V=Pij*innerMirrorQ;
Vvalidation=[];
for phaseLoop=1:phaseN
Vvalidation=[Vvalidation;Volt_(phaseLoop)*ones(semi_lineCount*200,1);];
end
error=abs((V-Vvalidation)./Vvalidation);
% error(isinf(error))=0;
error=sum(error)/length(Vvalidation)
%以下是验证部分
display('Finished.');
%计算场强
ABCy=imag(repmat(innerMirrorPos,1,length(vrfPos)));
ABCx=real(repmat(innerMirrorPos,1,length(vrfPos)));
y=imag(conj(matVrfPos'));
x=real(conj(matVrfPos'));
ERy=sum( ( (ABCy-y)./( (ABCy-y).^2+(ABCx-x).^2 )-(ABCy+y)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(real(innerMirrorQ),1,length(vrfPos))./2/pi/eslong,1 );
EIy=sum( ( (ABCy-y)./( (ABCy-y).^2+(ABCx-x).^2 )-(ABCy+y)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(imag(innerMirrorQ),1,length(vrfPos))./2/pi/eslong,1 );
ERx=sum( ( (ABCx-x)./( (ABCy-y).^2+(ABCx-x).^2 )-(ABCx-x)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(real(innerMirrorQ),1,length(vrfPos))./2/pi/eslong,1 );
EIx=sum( ( (ABCx-x)./( (ABCy-y).^2+(ABCx-x).^2 )-(ABCx-x)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(imag(innerMirrorQ),1,length(vrfPos))./2/pi/eslong,1 );
E2=sqrt(ERy.^2+EIy.^2+ERx.^2+EIx.^2+((ERy.^2-EIy.^2+ERx.^2-EIx.^2).^2+4*(ERy.*EIy+ERx.*EIx).^2).^.5);
E3=sqrt(ERy.^2+EIy.^2+ERx.^2+EIx.^2);
Emat=1/pi/2./eslong.*repmat(conj(innerMirrorQ'),length(vrfPos),1)./(vrf2ConductorDistance.^2).*(matVrfPos-repmat(conj(innerMirrorPos'),length(vrfPos),1))./vrf2ConductorDistance;
E=sum(Emat,2);
max(sqrt(2)*abs(E));
scatter(real(innerMirrorPos(1:length(innerMirrorPos)/1)),imag(innerMirrorPos(1:length(innerMirrorPos)/1)),[],'r');
axis equal
hold on;
scatter(real(vrfPos),imag(vrfPos),[],'k');

83
SelfAdaptSimulation.m
View File

@ -3,20 +3,22 @@ clear
%% 自适应模拟电荷法 %% 自适应模拟电荷法
% [1]. 任巍巍, 孙.A.宗.A., 一种较准确的分裂导线表面场强计算方法. 电网技术, 2006(04): 第92-96页. % [1]. 任巍巍, 孙.A.宗.A., 一种较准确的分裂导线表面场强计算方法. 电网技术, 2006(04): 第92-96页.
% [2]. 陈习文, 特高压直流输电线路电磁环境的研究, 2012, 北京交通大学. % [2]. 陈习文, 特高压直流输电线路电磁环境的研究, 2012, 北京交通大学.
% 采用文献的结果： 蒋兴良, 胡.A.舒.A., Analysis of Conductors＇ Surface Electric Field of UHVDC Transmission Lines Based on Optimized Charge Simulation Method. 高电压技术, 2008(12): p. 2547-2551.
%% %%
%设置几个参数 %设置几个参数
semi_lineDistance=457;%分裂间距 semi_lineDistance=450;%分裂间距
semi_lineCount=4;%分裂数 semi_lineCount=6;%分裂数
ConductorX=[-13720,0,13720];%导线距地高度 ConductorX=[-11000,11000];%导线间距
ConductorY=[20830,20830,20830];%导线间距 ConductorY=[22000,22000,];%导线距地高度
CSM_N=80;%每一个子导线的模拟电荷数 CSM_N=40;%每一个子导线的模拟电荷数
subconductorR=29.95;%子导线半径 subconductorR=16.8;%子导线半径
phaseN=3;%相数，单回三相 phaseN=2;%相数，单回三相
%% %%
%设置电压 %设置电压
Volt_=[765/sqrt(3);765/sqrt(3)*exp(1j*4/3*pi);765/sqrt(3)*exp(1j*2/3*pi);]; % Volt_=[1100/sqrt(3);1100/sqrt(3)*exp(1j*4/3*pi);1100/sqrt(3)*exp(1j*2/3*pi);];
Volt_=[800;-800;];
Volt=[]; Volt=[];
for vLoop=1:length(Volt_) for vLoop=1:phaseN
Volt=[Volt;Volt_(vLoop)*ones(CSM_N*semi_lineCount,1);]; Volt=[Volt;Volt_(vLoop)*ones(CSM_N*semi_lineCount,1);];
end end
%按分裂数和分裂导线间距布置单相线路导线 %按分裂数和分裂导线间距布置单相线路导线
@ -26,46 +28,47 @@ CSM_arc=2*pi/CSM_N;
%子导线中心到导线中心的距离 %子导线中心到导线中心的距离
R=semi_lineDistance/2/sin(arc/2); R=semi_lineDistance/2/sin(arc/2);
%计算模拟电荷的位置 %计算模拟电荷的位置
r1=20; r1=5;
error=10000; error=10000;
step=1/10; step=1/10;
maxLoop=round((subconductorR-r1)/step); maxLoop=round((subconductorR-r1)/step);
for Loop=1:maxLoop; for Loop=1:maxLoop;
simulationChargePos=ones(CSM_N,1); simulationChargePos=ones(CSM_N,1);
simulationChargeABCPos=[]; simulationChargeABCPos=[];
matchPos=[];
for I=1:CSM_N for I=1:CSM_N
simulationChargePos(I)=exp(1j*((I-1)*CSM_arc+CSM_arc/2))*r1;%逆时针转一个角度 simulationChargePos(I)=exp(1j*((I-1)*CSM_arc+CSM_arc/2))*r1;%逆时针转一个角度
end end
for phaseLoop=1:phaseN for phaseLoop=1:phaseN
for sC=1:semi_lineCount for sC=1:semi_lineCount
simulationChargeABCPos=[simulationChargeABCPos;simulationChargePos+ConductorX(phaseLoop)+1j*ConductorY(phaseLoop)+exp(1j*((sC-1)*arc+arc/2))*R];%移动到子导线中心 simulationChargeABCPos=[simulationChargeABCPos;simulationChargePos+ConductorX(phaseLoop)+1j*ConductorY(phaseLoop)+exp(1j*((sC-1)*arc+arc/2))*R];%移动到子导线中心
%同时计算匹配点的位置
matchPos=[matchPos;simulationChargePos/r1*subconductorR+ConductorX(phaseLoop)+1j*ConductorY(phaseLoop)+exp(1j*((sC-1)*arc+arc/2))*R];
end end
end end
% simulationChargeAPos=simulationChargePos+ConductorX(1)+1j*ConductorY(1); % simulationChargeAPos=simulationChargePos+ConductorX(1)+1j*ConductorY(1);
% simulationChargeBPos=simulationChargePos+ConductorX(2)+1j*ConductorY(2); % simulationChargeBPos=simulationChargePos+ConductorX(2)+1j*ConductorY(2);
% simulationChargePos=simulationChargeABCPos; % simulationChargePos=simulationChargeABCPos;
%计算电位系数 %计算电位系数
H=diag(imag(simulationChargeABCPos)); % H=diag(imag(simulationChargeABCPos));
r=subconductorR*eye(length(imag(simulationChargeABCPos)));%导线自几何均距 % r=subconductorR*eye(length(imag(simulationChargeABCPos)));%导线自几何均距
%导线与导线的距离 %导线与导线的距离
matSimulationChargePos=repmat(simulationChargeABCPos,1,length(simulationChargeABCPos)); matSimulationChargePos=repmat(simulationChargeABCPos,1,length(simulationChargeABCPos));
conductor2conductorDistance=matSimulationChargePos-conj(matSimulationChargePos'); matMatchPos=repmat(conj(matchPos'),length(simulationChargeABCPos),1);
conductor2conductorDistance=abs(conductor2conductorDistance-diag(diag(conductor2conductorDistance))); CMS2MatchPointDistance=abs(matSimulationChargePos-matMatchPos);
% conductor2conductorDistance=abs(conductor2conductorDistance-diag(diag(conductor2conductorDistance)));
matMirrorChargePos=conj(matSimulationChargePos);%虚部取负号 matMirrorChargePos=conj(matSimulationChargePos);%虚部取负号
conductor2MirrorDistance=matSimulationChargePos-conj(matMirrorChargePos'); mirrorCharge2MatchPointDistance=abs(matMirrorChargePos-matMatchPos);
conductor2MirrorDistance=abs(conductor2MirrorDistance-diag(diag(conductor2MirrorDistance))); % conductor2MirrorDistance=abs(conductor2MirrorDistance-diag(diag(conductor2MirrorDistance)));
eslong=8.854187817*10; eslong=8.854187817*10^-12*1000;
P1=1/2/pi/eslong*log(2*H./r); % P1=1/2/pi/eslong*log(2*H./r);
P1(isnan(P1))=0; % P1(isnan(P1))=0;
P2=1/2/pi/eslong*log(conductor2MirrorDistance./conductor2conductorDistance); P2=1/2/pi/eslong*log(mirrorCharge2MatchPointDistance./CMS2MatchPointDistance);
P2(isnan(P2))=0; % P2(isnan(P2))=0;
P=P1+P2; % P=P1+P2;
P=P2;
%求电荷 %求电荷
QRI=P\Volt; QRI=P\Volt;
%以下是验证部分
if error<0.0001
break;
end
%选检验导线上一个角度 %选检验导线上一个角度
vrfRelA=linspace(0,2*pi,200)';%vrf=verify vrfRelA=linspace(0,2*pi,200)';%vrf=verify
%计算检验点相对于子导线的位置 %计算检验点相对于子导线的位置
@ -88,8 +91,14 @@ for Loop=1:maxLoop;
for phaseLoop=1:phaseN for phaseLoop=1:phaseN
Vvalidation=[Vvalidation;Volt_(phaseLoop)*ones(semi_lineCount*200,1);]; Vvalidation=[Vvalidation;Volt_(phaseLoop)*ones(semi_lineCount*200,1);];
end end
error=sum(abs((V-Vvalidation)./Vvalidation))/length(Vvalidation); error=abs((V-Vvalidation)./Vvalidation);
r1=r1+step; % error(isinf(error))=0;
error=sum(error)/length(Vvalidation)
%以下是验证部分
if error<0.001
break;
end
r1=r1+1*step;
end end
display('Finished.'); display('Finished.');
if Loop<maxLoop if Loop<maxLoop
@ -101,12 +110,18 @@ ABCy=imag(repmat(simulationChargeABCPos,1,length(vrfPos)));
ABCx=real(repmat(simulationChargeABCPos,1,length(vrfPos))); ABCx=real(repmat(simulationChargeABCPos,1,length(vrfPos)));
y=imag(conj(matVrfPos')); y=imag(conj(matVrfPos'));
x=real(conj(matVrfPos')); x=real(conj(matVrfPos'));
ERy=sum( ( (ABCy-y)./( (ABCy-y).^2+(ABCx-x).^2 )+(ABCy+y)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(real(QRI),1,length(vrfPos))./2/pi/eslong,1 ); ERy=sum( ( (ABCy-y)./( (ABCy-y).^2+(ABCx-x).^2 )-(ABCy+y)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(real(QRI),1,length(vrfPos))./2/pi/eslong,1 );
EIy=sum( ( (ABCy-y)./( (ABCy-y).^2+(ABCx-x).^2 )+(ABCy+y)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(imag(QRI),1,length(vrfPos))./2/pi/eslong,1 ); EIy=sum( ( (ABCy-y)./( (ABCy-y).^2+(ABCx-x).^2 )-(ABCy+y)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(imag(QRI),1,length(vrfPos))./2/pi/eslong,1 );
ERx=sum( ( (ABCx-x)./( (ABCy-y).^2+(ABCx-x).^2 )+(ABCx-x)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(real(QRI),1,length(vrfPos))./2/pi/eslong,1 ); ERx=sum( ( (ABCx-x)./( (ABCy-y).^2+(ABCx-x).^2 )-(ABCx-x)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(real(QRI),1,length(vrfPos))./2/pi/eslong,1 );
EIx=sum( ( (ABCx-x)./( (ABCy-y).^2+(ABCx-x).^2 )+(ABCx-x)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(imag(QRI),1,length(vrfPos))./2/pi/eslong,1 ); EIx=sum( ( (ABCx-x)./( (ABCy-y).^2+(ABCx-x).^2 )-(ABCx-x)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(imag(QRI),1,length(vrfPos))./2/pi/eslong,1 );
E=sqrt(ERy.^2+EIy.^2+ERx.^2+EIx.^2); E2=sqrt(ERy.^2+EIy.^2+ERx.^2+EIx.^2+((ERy.^2-EIy.^2+ERx.^2-EIx.^2).^2+4*(ERy.*EIy+ERx.*EIx).^2).^.5);
max(E); E3=sqrt(ERy.^2+EIy.^2+ERx.^2+EIx.^2);
Emat=1/pi/2./eslong.*repmat(conj(QRI'),length(vrfPos),1)./(vrf2ConductorDistance.^2).*(matVrfPos-repmat(conj(simulationChargeABCPos'),length(vrfPos),1))./vrf2ConductorDistance;
E=sum(Emat,2);
max(sqrt(2)*abs(E));
scatter(real(simulationChargeABCPos(1:length(simulationChargeABCPos)/1)),imag(simulationChargeABCPos(1:length(simulationChargeABCPos)/1)),[],'r'); scatter(real(simulationChargeABCPos(1:length(simulationChargeABCPos)/1)),imag(simulationChargeABCPos(1:length(simulationChargeABCPos)/1)),[],'r');
axis equal axis equal