diff --git a/GradualMirror.m b/GradualMirror.m new file mode 100644 index 0000000..e1506b7 --- /dev/null +++ b/GradualMirror.m @@ -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'); \ No newline at end of file diff --git a/SelfAdaptSimulation.m b/SelfAdaptSimulation.m index b1c0979..94d7431 100644 --- a/SelfAdaptSimulation.m +++ b/SelfAdaptSimulation.m @@ -3,20 +3,22 @@ clear %% 自适应模拟电荷法 % [1]. 任巍巍, 孙.A.宗.A., 一种较准确的分裂导线表面场强计算方法. 电网技术, 2006(04): 第92-96页. % [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_lineCount=4;%分裂数 -ConductorX=[-13720,0,13720];%导线距地高度 -ConductorY=[20830,20830,20830];%导线间距 -CSM_N=80;%每一个子导线的模拟电荷数 -subconductorR=29.95;%子导线半径 -phaseN=3;%相数，单回三相 +semi_lineDistance=450;%分裂间距 +semi_lineCount=6;%分裂数 +ConductorX=[-11000,11000];%导线间距 +ConductorY=[22000,22000,];%导线距地高度 +CSM_N=40;%每一个子导线的模拟电荷数 +subconductorR=16.8;%子导线半径 +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=[]; -for vLoop=1:length(Volt_) +for vLoop=1:phaseN Volt=[Volt;Volt_(vLoop)*ones(CSM_N*semi_lineCount,1);]; end %按分裂数和分裂导线间距布置单相线路导线 @@ -26,46 +28,47 @@ CSM_arc=2*pi/CSM_N; %子导线中心到导线中心的距离 R=semi_lineDistance/2/sin(arc/2); %计算模拟电荷的位置 -r1=20; +r1=5; error=10000; step=1/10; maxLoop=round((subconductorR-r1)/step); for Loop=1:maxLoop; simulationChargePos=ones(CSM_N,1); simulationChargeABCPos=[]; + matchPos=[]; for I=1:CSM_N simulationChargePos(I)=exp(1j*((I-1)*CSM_arc+CSM_arc/2))*r1;%逆时针转一个角度 end for phaseLoop=1:phaseN for sC=1:semi_lineCount 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 % simulationChargeAPos=simulationChargePos+ConductorX(1)+1j*ConductorY(1); % simulationChargeBPos=simulationChargePos+ConductorX(2)+1j*ConductorY(2); % simulationChargePos=simulationChargeABCPos; %计算电位系数 - H=diag(imag(simulationChargeABCPos)); - r=subconductorR*eye(length(imag(simulationChargeABCPos)));%导线自几何均距 +% H=diag(imag(simulationChargeABCPos)); +% r=subconductorR*eye(length(imag(simulationChargeABCPos)));%导线自几何均距 %导线与导线的距离 matSimulationChargePos=repmat(simulationChargeABCPos,1,length(simulationChargeABCPos)); - conductor2conductorDistance=matSimulationChargePos-conj(matSimulationChargePos'); - conductor2conductorDistance=abs(conductor2conductorDistance-diag(diag(conductor2conductorDistance))); + matMatchPos=repmat(conj(matchPos'),length(simulationChargeABCPos),1); + CMS2MatchPointDistance=abs(matSimulationChargePos-matMatchPos); +% conductor2conductorDistance=abs(conductor2conductorDistance-diag(diag(conductor2conductorDistance))); matMirrorChargePos=conj(matSimulationChargePos);%虚部取负号 - conductor2MirrorDistance=matSimulationChargePos-conj(matMirrorChargePos'); - conductor2MirrorDistance=abs(conductor2MirrorDistance-diag(diag(conductor2MirrorDistance))); - eslong=8.854187817*10; - P1=1/2/pi/eslong*log(2*H./r); - P1(isnan(P1))=0; - P2=1/2/pi/eslong*log(conductor2MirrorDistance./conductor2conductorDistance); - P2(isnan(P2))=0; - P=P1+P2; + mirrorCharge2MatchPointDistance=abs(matMirrorChargePos-matMatchPos); +% conductor2MirrorDistance=abs(conductor2MirrorDistance-diag(diag(conductor2MirrorDistance))); + eslong=8.854187817*10^-12*1000; +% P1=1/2/pi/eslong*log(2*H./r); +% P1(isnan(P1))=0; + P2=1/2/pi/eslong*log(mirrorCharge2MatchPointDistance./CMS2MatchPointDistance); +% P2(isnan(P2))=0; +% P=P1+P2; + P=P2; %求电荷 QRI=P\Volt; - %以下是验证部分 - if error<0.0001 - break; - end %选检验导线上一个角度 vrfRelA=linspace(0,2*pi,200)';%vrf=verify %计算检验点相对于子导线的位置 @@ -88,8 +91,14 @@ for Loop=1:maxLoop; for phaseLoop=1:phaseN Vvalidation=[Vvalidation;Volt_(phaseLoop)*ones(semi_lineCount*200,1);]; end - error=sum(abs((V-Vvalidation)./Vvalidation))/length(Vvalidation); - r1=r1+step; + error=abs((V-Vvalidation)./Vvalidation); +% error(isinf(error))=0; + error=sum(error)/length(Vvalidation) + %以下是验证部分 + if error<0.001 + break; + end + r1=r1+1*step; end display('Finished.'); if Loop