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

完善了双回路EGM模型代码。

This commit is contained in:
facat 2021-09-22 00:18:06 +08:00
parent dd44de030e
commit 392eeb0168
2 changed files with 229 additions and 142 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
core.py
View File

@ -25,7 +25,7 @@ class Draw:
rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
msp.add_circle((rs_x, rs_y), rs, dxfattribs={"color": color})
msp.add_line((0, 0), (rs_x, rs_y)) # 地线
msp.add_circle((rc_x, rc_y), rc, dxfattribs={"color": color+2})
msp.add_circle((rc_x, rc_y), rc, dxfattribs={"color": color + 2})
msp.add_line((rc_x, 0), (rc_x, rc_y)) # 导线
msp.add_line((rs_x, rs_y), (rc_x, rc_y))
# 角度线
@ -36,9 +36,12 @@ class Draw:
if rg_type == "g":
msp.add_line((0, rg), (2000, rg), dxfattribs={"color": color})
circle_line_section = solve_circle_line_intersection(rc, rg, rc_x, rc_y)
msp.add_line(
(rc_x, rc_y), circle_line_section, dxfattribs={"color": color}
) # 导线和圆的交点
if not circle_line_section:
pass
else:
msp.add_line(
(rc_x, rc_y), circle_line_section, dxfattribs={"color": color}
) # 导线和圆的交点
if rg_type == "c":
msp.add_circle((rg_x, rg_y), rg, dxfattribs={"color": color})
rg_rc_intersection = solve_circle_intersection(
@ -138,7 +141,7 @@ def rg_fun(i_curt, h_cav, u_ph, typ="g"):
def intersection_angle(
rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, rg_type
rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_surface, rg_type
): # 暴露弧的角度
rs = rs_fun(i_curt)
rc = rc_fun(i_curt, u_ph)
@ -155,6 +158,18 @@ def intersection_angle(
circle_line_or_rg_intersection = solve_circle_intersection(
rg, rc, rg_x, rg_y, rc_x, rc_y
) # 两圆的交点
(
circle_line_or_rg_intersection_x,
circle_line_or_rg_intersection_y,
) = circle_line_or_rg_intersection
if (
ground_surface(circle_line_or_rg_intersection_x)
> circle_line_or_rg_intersection_y
): # 交点在地面线以下，就可以不积分
# 找到暴露弧和地面线的交点
circle_line_or_rg_intersection = circle_ground_surface_intersection(
rc, rc_x, rc_y, ground_surface
)
np_circle_intersection = np.array(circle_intersection)
theta2_line = np_circle_intersection - np.array([rc_x, rc_y])
theta2 = math.atan(theta2_line[1] / theta2_line[0])
@ -224,9 +239,11 @@ def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y): # 对θ进行积分
return r
def bd_area(i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, rg_type): # 暴露弧的投影面积
def bd_area(
i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, ground_surface, rg_type
): # 暴露弧的投影面积
theta1, theta2 = intersection_angle(
rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, rg_type
rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_surface, rg_type
) # θ角度
theta_fineness = 0.01
rc = rc_fun(i_curt, u_ph)
@ -304,3 +321,19 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
def func_ng(td): # 地闪密度
return 0.023 * (td ** 1.3)
# 圆和地面线的交点，只去正x轴上的。
def circle_ground_surface_intersection(radius, center_x, center_y, ground_surface):
# 最笨的办法，一个个去试
x_series = np.linspace(0, radius, int(radius / 0.001)) + center_x
part_to_be_squared = radius ** 2 - (x_series - center_x) ** 2 # 有可能出现-0.00001的数值，只是一个数值稳定问题。
part_to_be_squared[
(part_to_be_squared < 0) & (abs(part_to_be_squared) < 1e-3)
] = 0 # 强制为0
y_series = center_y - part_to_be_squared ** 0.5
ground_surface_y = ground_surface(x_series)
equal_location = np.abs(ground_surface_y - y_series) < 0.5
r_x = x_series[equal_location][0]
r_y = ground_surface(r_x)
return r_x, r_y

324
main.py
View File

@ -1,158 +1,212 @@
import numpy as np
from core import *
import timeit
def egm():
avr_n_sf = 0 # 考虑电压的影响计算的跳闸率
voltage_n = 1 # 工作电压分成多少份来计算
ng = func_ng(20)
h_g_avr_sag = 11.67 * 2 / 3
h_c_avr_sag = 14.43 * 2 / 3
h_whole = 140 # 杆塔全高
voltage_n = 3 # 工作电压分成多少份来计算
td = 20 # 雷暴日
insulator_c_len = 6.8 # 串子绝缘长度
string_c_len = 9.2
string_g_len = 0.5
rc_x = -0.0 # 导地线水平距离
rs_x = 0
rg_x = 0
vertical_dgc = 2.7 # 导地线挂点垂直距离
h_g_avr_sag = 11.67 * 2 / 3
h_c_avr_sag = 14.43 * 2 / 3
rs_y = h_whole - string_g_len - h_g_avr_sag # 地线对地平均高
rc_y = h_whole - string_c_len - vertical_dgc - h_c_avr_sag # 导线对地平均高
rg_y = rc_y - 20
shield_angle = (
math.atan(rc_x / (vertical_dgc + string_c_len)) * 180 / math.pi
) # 保护角
print(f"保护角{shield_angle:.3f}°")
rg_type = "c"
for u_bar in range(voltage_n):
u_ph = (
math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732
) # 运行相电压
# 迭代法计算最大电流
i_max = 0
i_min = min_i(insulator_c_len, u_ph / 1.732)
_min_i = i_min # 尝试的最小电流
_max_i = 200 # 尝试的最大电流
cad = Draw()
# cad.draw(i_min, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 2)
for i_bar in np.linspace(_min_i, _max_i, int((_max_i - _min_i) / 0.1)): # 雷电流
print(f"尝试计算电流为{i_bar:.2f}")
rs = rs_fun(i_bar)
rc = rc_fun(i_bar, u_ph)
rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type)
#######
# cccCount += 1
# if cccCount % 30 == 0:
# import core
#
# core.gMSP.add_circle((0, h_gav), rs)
# core.gMSP.add_circle(
# (dgc, h_cav), rc_fun(i_bar, -u_ph), dxfattribs={"color": 4}
# )
# core.gMSP.add_circle((dgc, h_cav), rc)
#######
rg_rc_circle_intersection = solve_circle_intersection(
rs, rc, rs_x, rs_y, rc_x, rc_y
)
if not rg_rc_circle_intersection: # if circle_intersection is []
print("保护弧和暴露弧无交点，检查设置参数。程序退出。")
continue
circle_rc_line_or_rg_intersection = None
if rg_type == "g":
circle_rc_line_or_rg_intersection = solve_circle_line_intersection(
rc, rg, rc_x, rc_y
gc_x = [17.9, 16, 15, 16]
# 以后考虑地形角度,地面线
def ground_surface(x):
return 0
gc_y = [
h_whole - string_g_len - h_g_avr_sag, # 地线对地平均高
h_whole - string_c_len - h_c_avr_sag - 2.7, # 导线对地平均高
h_whole - string_c_len - h_c_avr_sag - 20, # 导线对地平均高
h_whole - string_c_len - h_c_avr_sag - 35.7, # 导线对地平均高
]
if len(gc_y) > 2: # 双回路
phase_n = 3 # 边相导线数量
else:
phase_n = 1
#########################################################
rg_type = None
# 以上是需要设置的参数
avr_n_sf = 0 # 考虑电压的影响计算的跳闸率
rg_x = None
rg_y = None
cad = Draw()
for phase_conductor in range(phase_n):
rs_x = gc_x[phase_conductor]
rs_y = gc_y[phase_conductor]
rc_x = gc_x[phase_conductor + 1]
rc_y = gc_y[phase_conductor + 1]
if phase_n == 1:
rg_type = "g"
if phase_n > 1: # 多回路
if phase_conductor < 2:
rg_type = "c"
rg_x = gc_x[phase_conductor + 2]
rg_y = gc_y[phase_conductor + 2]
else:
rg_type = "g"
# TODO 保护角公式可能有问题，后面改
shield_angle = (
math.atan(rc_x / ((rc_y - rs_y) + string_c_len)) * 180 / math.pi
) # 保护角
print(f"保护角{shield_angle:.3f}°")
print(f"最低相防护标识{rg_type}")
ng = func_ng(td)
for u_bar in range(voltage_n):
u_ph = (
math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732
) # 运行相电压
print(f"计算第{phase_conductor + 1}相，电压为{u_ph:.2f}kV")
# 迭代法计算最大电流
i_max = 0
i_min = min_i(insulator_c_len, u_ph / 1.732)
_min_i = i_min # 尝试的最小电流
_max_i = 200 # 尝试的最大电流
# cad.draw(i_min, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 2)
for i_bar in np.linspace(
_min_i, _max_i, int((_max_i - _min_i) / 0.1)
): # 雷电流
# print(f"尝试计算电流为{i_bar:.2f}")
rs = rs_fun(i_bar)
rc = rc_fun(i_bar, u_ph)
rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type)
#######
# cccCount += 1
# if cccCount % 30 == 0:
# import core
#
# core.gMSP.add_circle((0, h_gav), rs)
# core.gMSP.add_circle(
# (dgc, h_cav), rc_fun(i_bar, -u_ph), dxfattribs={"color": 4}
# )
# core.gMSP.add_circle((dgc, h_cav), rc)
#######
rg_rc_circle_intersection = solve_circle_intersection(
rs, rc, rs_x, rs_y, rc_x, rc_y
)
elif rg_type == "c":
circle_rc_line_or_rg_intersection = solve_circle_intersection(
rg, rc, rg_x, rg_y, rc_x, rc_y
)
if not circle_rc_line_or_rg_intersection:
continue
min_distance_intersection = (
np.sum(
(
np.array(rg_rc_circle_intersection)
- np.array(circle_rc_line_or_rg_intersection)
)
** 2
)
** 0.5
) # 计算两圆交点和地面直线交点的最小距离
i_max = i_bar
if min_distance_intersection < 0.1:
break
# 判断是否以完全被保护
if rg_rc_circle_intersection[1] < circle_rc_line_or_rg_intersection[1]:
circle_rs_line_or_rg_intersection = None
i_max = i_bar
if not rg_rc_circle_intersection: # if circle_intersection is []
print("保护弧和暴露弧无交点，检查设置参数。")
continue
circle_rc_line_or_rg_intersection = None
if rg_type == "g":
circle_rs_line_or_rg_intersection = solve_circle_line_intersection(
rs, rg, rs_x, rs_y
circle_rc_line_or_rg_intersection = solve_circle_line_intersection(
rc, rg, rc_x, rc_y
)
if rg_type == "c":
circle_rs_line_or_rg_intersection = solve_circle_intersection(
rs, rg, rs_x, rs_y, rg_x, rg_y
elif rg_type == "c":
circle_rc_line_or_rg_intersection = solve_circle_intersection(
rg, rc, rg_x, rg_y, rc_x, rc_y
)
# 判断与保护弧的交点是否在暴露弧外面
distance = (
if not circle_rc_line_or_rg_intersection:
# 暴露弧和捕捉弧无交点
if rg_type == "g":
if rg > rc_y:
i_min = i_bar
print(f"捕捉弧在暴露弧之上，设置最小电流为{i_min:.2f}")
else:
print("暴露弧和捕捉弧无交点，检查设置参数。")
continue
else:
print("暴露弧和捕捉弧无交点，检查设置参数。")
continue
min_distance_intersection = (
np.sum(
(
np.array(circle_rs_line_or_rg_intersection)
- np.array([rc_x, rc_y])
np.array(rg_rc_circle_intersection)
- np.array(circle_rc_line_or_rg_intersection)
)
** 2
)
** 0.5
)
if distance > rc:
print("暴露弧已经完全被屏蔽")
) # 计算两圆交点和地面直线交点的最小距离
if min_distance_intersection < 0.1:
break
cad.draw(i_min, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 2)
cad.draw(i_max, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 6)
cad.save()
# 判断是否导线已经被完全保护
if abs(i_max - _max_i) < 1e-5:
print("无法找到最大电流，可能是杆塔较高。")
print(f"最大电流设置为自然界最大电流{i_max}kA")
print(f"最大电流为{i_max:.2f}")
print(f"最小电流为{i_min:.2f}")
curt_fineness = 0.1 # 电流积分细度
if i_min > i_max or abs(i_min - i_max) < curt_fineness:
print("最大电流小于最小电流，没有暴露弧，程序结束。")
return
# 开始积分
curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份
i_curt_samples, d_curt = np.linspace(
i_min, i_max, curt_segment_n + 1, retstep=True
)
bd_area_vec = np.vectorize(bd_area)
cal_bd_np = bd_area_vec(
i_curt_samples, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, rg_type
) * thunder_density(i_curt_samples)
calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt
# for i_curt in i_curt_samples[:-1]:
# cal_bd_first = bd_area(i_curt, u_ph, dgc, h_gav, h_cav)
# cal_bd_second = bd_area(i_curt + d_curt, u_ph, dgc, h_gav, h_cav)
# cal_thunder_density_first = thunder_density(i_curt)
# cal_thunder_density_second = thunder_density(i_curt + d_curt)
# calculus += (
# (
# cal_bd_first * cal_thunder_density_first
# + cal_bd_second * cal_thunder_density_second
# )
# / 2
# * d_curt
# )
# if abs(calculus-0.05812740052770032)<1e-5:
# abc=123
# pass
n_sf = (
2 * ng / 10 * calculus
) # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
avr_n_sf += n_sf / voltage_n
print(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}")
print(f"跳闸率是{avr_n_sf:.6}")
# 判断是否以完全被保护
if rg_rc_circle_intersection[1] < circle_rc_line_or_rg_intersection[1]:
circle_rs_line_or_rg_intersection = None
if rg_type == "g":
circle_rs_line_or_rg_intersection = (
solve_circle_line_intersection(rs, rg, rs_x, rs_y)
)
if rg_type == "c":
circle_rs_line_or_rg_intersection = solve_circle_intersection(
rs, rg, rs_x, rs_y, rg_x, rg_y
)
# 判断与保护弧的交点是否在暴露弧外面
distance = (
np.sum(
(
np.array(circle_rs_line_or_rg_intersection)
- np.array([rc_x, rc_y])
)
** 2
)
** 0.5
)
if distance > rc:
print("暴露弧已经完全被屏蔽")
break
if phase_conductor == 1:
cad.draw(i_min, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 2)
cad.draw(i_max, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 6)
cad.save()
# 判断是否导线已经被完全保护
if abs(i_max - _max_i) < 1e-5:
print("无法找到最大电流，可能是杆塔较高。")
print(f"最大电流设置为自然界最大电流{i_max}kA")
print(f"最大电流为{i_max:.2f}")
print(f"最小电流为{i_min:.2f}")
curt_fineness = 0.1 # 电流积分细度
if i_min > i_max or abs(i_min - i_max) < curt_fineness:
print("最大电流小于最小电流，没有暴露弧。")
continue
# 开始积分
curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份
i_curt_samples, d_curt = np.linspace(
i_min, i_max, curt_segment_n + 1, retstep=True
)
bd_area_vec = np.vectorize(bd_area)
cal_bd_np = (
bd_area_vec(
i_curt_samples,
u_ph,
rc_x,
rc_y,
rs_x,
rs_y,
rg_x,
rg_y,
ground_surface,
rg_type,
)
* thunder_density(i_curt_samples)
)
calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt
# for i_curt in i_curt_samples[:-1]:
# cal_bd_first = bd_area(i_curt, u_ph, dgc, h_gav, h_cav)
# cal_bd_second = bd_area(i_curt + d_curt, u_ph, dgc, h_gav, h_cav)
# cal_thunder_density_first = thunder_density(i_curt)
# cal_thunder_density_second = thunder_density(i_curt + d_curt)
# calculus += (
# (
# cal_bd_first * cal_thunder_density_first
# + cal_bd_second * cal_thunder_density_second
# )
# / 2
# * d_curt
# )
# if abs(calculus-0.05812740052770032)<1e-5:
# abc=123
# pass
n_sf = (
2 * ng / 10 * calculus
) # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
avr_n_sf += n_sf / voltage_n
print(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}")
print(f"跳闸率是{avr_n_sf:.6f}")
def speed():