371 lines
16 KiB
Python
371 lines
16 KiB
Python
import math
|
||
import os.path
|
||
import sys
|
||
import time
|
||
import tomli
|
||
from loguru import logger
|
||
from core import *
|
||
import timeit
|
||
from animation import Animation
|
||
|
||
|
||
# 打印参数
|
||
def parameter_display(para_dis: Parameter):
|
||
logger.info(f"额定电压 kV {para_dis.rated_voltage}")
|
||
logger.info(f"导线弧垂 m {para_dis.h_c_sag}")
|
||
logger.info(f"地线弧垂 m {para_dis.h_g_sag}")
|
||
logger.info(f"全塔高 m {para_dis.h_arm[0]}")
|
||
logger.info(f"串绝缘距离 m {para_dis.insulator_c_len}")
|
||
logger.info(f"导线串长 m {para_dis.string_c_len}")
|
||
logger.info(f"地线串长 m {para_dis.string_g_len}")
|
||
logger.info(f"挂点垂直坐标 m {para_dis.h_arm}")
|
||
logger.info(f"挂点水平坐标 m {para_dis.gc_x}")
|
||
logger.info(f"地面倾角 ° {[an * 180 / math.pi for an in para_dis.ground_angels]}")
|
||
logger.info(f"海拔高度 m {para_dis.altitude}")
|
||
if para_dis.ng > 0:
|
||
logger.info("不采用雷暴日计算地闪密度和雷电流密度")
|
||
logger.info(f"地闪密度 次/(每平方公里·每年) {para_dis.ng}")
|
||
logger.info(f"概率密度曲线系数a {para_dis.Ip_a}")
|
||
logger.info(f"概率密度曲线系数b {para_dis.Ip_b}")
|
||
pass
|
||
else:
|
||
logger.info(f"雷暴日 d {para_dis.td}")
|
||
|
||
|
||
def read_parameter(toml_file_path):
|
||
with open(toml_file_path, "rb") as toml_fs:
|
||
toml_dict = tomli.load(toml_fs)
|
||
toml_parameter = toml_dict["parameter"]
|
||
para.h_g_sag = toml_parameter["h_g_sag"] # 地线弧垂
|
||
para.h_c_sag = toml_parameter["h_c_sag"] # 导线弧垂
|
||
# para.h_whole = toml_parameter["h_whole"] # 杆塔全高
|
||
para.td = toml_parameter["td"] # 雷暴日
|
||
para.insulator_c_len = toml_parameter["insulator_c_len"] # 串子绝缘长度
|
||
para.string_c_len = toml_parameter["string_c_len"]
|
||
para.string_g_len = toml_parameter["string_g_len"]
|
||
para.gc_x = toml_parameter["gc_x"] # 导、地线水平坐标
|
||
para.ground_angels = [
|
||
angel / 180 * math.pi for angel in toml_parameter["ground_angels"]
|
||
] # 地面倾角，向下为正
|
||
para.h_arm = toml_parameter["h_arm"]
|
||
para.altitude = toml_parameter["altitude"]
|
||
para.rated_voltage = toml_parameter["rated_voltage"]
|
||
toml_advance = toml_dict["advance"]
|
||
para.ng = toml_advance["ng"] # 地闪密度
|
||
para.Ip_a = toml_advance["Ip_a"] # 概率密度曲线系数a
|
||
para.Ip_b = toml_advance["Ip_b"] # 概率密度曲线系数b
|
||
toml_optional = toml_dict["optional"]
|
||
para.voltage_n = toml_optional["voltage_n"] # 工作电压分成多少份来计算
|
||
para.max_i = toml_optional["max_i"]
|
||
|
||
|
||
def egm():
|
||
if len(sys.argv) < 2:
|
||
toml_file_path = r"内自500kV-ZCK上相.toml"
|
||
else:
|
||
toml_file_path = sys.argv[1]
|
||
if not os.path.exists(toml_file_path):
|
||
logger.info(f"无法找到数据文件{toml_file_path}，程序退出。")
|
||
sys.exit(0)
|
||
logger.info(f"读取文件{toml_file_path}")
|
||
read_parameter(toml_file_path)
|
||
#########################################################
|
||
# 以上是需要设置的参数
|
||
parameter_display(para)
|
||
h_whole = para.h_arm[0] # 挂点高
|
||
string_g_len = para.string_g_len
|
||
string_c_len = para.string_c_len
|
||
h_g_sag = para.h_g_sag
|
||
h_c_sag = para.h_c_sag
|
||
gc_x = para.gc_x
|
||
h_arm = para.h_arm
|
||
gc_y = [
|
||
h_whole - string_g_len - h_g_sag * 2 / 3, # 地线对地平均高
|
||
]
|
||
if len(h_arm) > 1:
|
||
for hoo in h_arm[1:]:
|
||
gc_y.append(hoo - string_c_len - h_c_sag * 2 / 3) # 导线平均高
|
||
if len(gc_y) > 2: # 双回路
|
||
phase_n = 3 # 边相导线数量
|
||
else:
|
||
phase_n = 1
|
||
# 地闪密度 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
|
||
td = para.td
|
||
ng = func_ng(td)
|
||
avr_n_sf = 0 # 考虑电压的影响计算的跳闸率
|
||
ground_angels = para.ground_angels
|
||
# 初始化动画
|
||
animate = Animation()
|
||
animate.disable(False)
|
||
# animate.show()
|
||
for ground_angel in ground_angels:
|
||
logger.info(f"地面倾角{ground_angel/math.pi*180:.3f}°")
|
||
rg_type = None
|
||
rg_x = None
|
||
rg_y = None
|
||
voltage_n = para.voltage_n
|
||
n_sf_phases = np.zeros((phase_n, voltage_n)) # 存储每一相的跳闸率
|
||
if np.any(np.array(gc_y) < 0):
|
||
logger.info("导线可能掉地面下了，程序退出。")
|
||
return 0
|
||
for phase_conductor_foo in range(phase_n):
|
||
exposed_curve_shielded = False
|
||
rs_x = gc_x[phase_conductor_foo]
|
||
rs_y = gc_y[phase_conductor_foo]
|
||
rc_x = gc_x[phase_conductor_foo + 1]
|
||
rc_y = gc_y[phase_conductor_foo + 1]
|
||
if phase_n == 1:
|
||
rg_type = "g"
|
||
if phase_n > 1: # 多回路
|
||
if phase_conductor_foo < 2:
|
||
rg_type = "c" # 捕捉弧由下面一相导线的击距代替
|
||
rg_x = gc_x[phase_conductor_foo + 2]
|
||
rg_y = gc_y[phase_conductor_foo + 2]
|
||
else:
|
||
rg_type = "g"
|
||
# TODO 保护角公式可能有问题，后面改
|
||
shield_angle_at_avg_height = (
|
||
math.atan(
|
||
(rc_x - rs_x)
|
||
/ (
|
||
(h_arm[0] - string_g_len - h_arm[phase_conductor_foo + 1])
|
||
+ string_c_len
|
||
)
|
||
)
|
||
* 180
|
||
/ math.pi
|
||
) # 挂点处保护角
|
||
logger.info(f"挂点处保护角{shield_angle_at_avg_height:.3f}°")
|
||
logger.debug(f"最低相防护标识{rg_type}")
|
||
rated_voltage = para.rated_voltage
|
||
for u_bar in range(voltage_n): # 计算不同工作电压下的跳闸率
|
||
# TODO 需要区分交、直流
|
||
# u_ph = (
|
||
# math.sqrt(2)
|
||
# * rated_voltage
|
||
# * math.cos(2 * math.pi / voltage_n * u_bar)
|
||
# / 1.732
|
||
# ) # 运行相电压
|
||
u_ph = rated_voltage / 1.732
|
||
logger.info(f"计算第{phase_conductor_foo + 1}相，电压为{u_ph:.2f}kV")
|
||
# 迭代法计算最大电流
|
||
i_max = 0
|
||
insulator_c_len = para.insulator_c_len
|
||
# i_min = min_i(insulator_c_len, u_ph / 1.732)
|
||
# TODO 需要考虑交、直流
|
||
i_min = min_i(insulator_c_len, u_ph)
|
||
_min_i = i_min # 尝试的最小电流
|
||
_max_i = para.max_i # 尝试的最大电流
|
||
# 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) / 1)
|
||
): # 雷电流
|
||
logger.info(f"尝试计算电流为{i_bar:.2f}")
|
||
rs = rs_fun(i_bar)
|
||
animate.add_rs(rs, rs_x, rs_y)
|
||
rc = rc_fun(i_bar, u_ph)
|
||
animate.add_rc(rc, rc_x, rc_y)
|
||
rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type)
|
||
rg_line_func = None
|
||
if rg_type == "g":
|
||
rg_line_func = rg_line_function_factory(rg, ground_angel)
|
||
animate.add_rg_line(rg_line_func)
|
||
rs_rc_circle_intersection = solve_circle_intersection(
|
||
rs, rc, rs_x, rs_y, rc_x, rc_y
|
||
)
|
||
i_max = i_bar
|
||
if not rs_rc_circle_intersection: # if circle_intersection is []
|
||
logger.debug("保护弧和暴露弧无交点，检查设置参数。")
|
||
continue
|
||
circle_rc_or_rg_line_intersection = None
|
||
if rg_type == "g":
|
||
circle_rc_or_rg_line_intersection = (
|
||
solve_circle_line_intersection(rc, rc_x, rc_y, rg_line_func)
|
||
)
|
||
elif rg_type == "c":
|
||
circle_rc_or_rg_line_intersection = solve_circle_intersection(
|
||
rg, rc, rg_x, rg_y, rc_x, rc_y
|
||
)
|
||
if not circle_rc_or_rg_line_intersection:
|
||
# 暴露弧和捕捉弧无交点
|
||
if rg_type == "g":
|
||
if rg_line_func(rc_x) > rc_y:
|
||
i_min = i_bar # 用于后面判断最小和最大电流是否相等，相等意味着暴露弧一直被屏蔽
|
||
logger.info(
|
||
f"捕捉面在暴露弧之上，设置最小电流为{i_min:.2f}"
|
||
)
|
||
else:
|
||
logger.info("暴露弧和地面捕捉弧无交点，检查设置参数。")
|
||
continue
|
||
else:
|
||
logger.info(
|
||
"上面的导地线无法保护下面的导地线，检查设置参数。"
|
||
)
|
||
continue
|
||
animate.add_expose_area(
|
||
rc_x,
|
||
rc_y,
|
||
*rs_rc_circle_intersection,
|
||
*circle_rc_or_rg_line_intersection,
|
||
)
|
||
cad = Draw()
|
||
cad.draw(
|
||
i_min,
|
||
u_ph,
|
||
rs_x,
|
||
rs_y,
|
||
rc_x,
|
||
rc_y,
|
||
rg_x,
|
||
rg_y,
|
||
rg_type,
|
||
ground_angel,
|
||
2,
|
||
) # 最小电流时
|
||
cad.draw(
|
||
i_max,
|
||
u_ph,
|
||
rs_x,
|
||
rs_y,
|
||
rc_x,
|
||
rc_y,
|
||
rg_x,
|
||
rg_y,
|
||
rg_type,
|
||
ground_angel,
|
||
6,
|
||
) # 最大电流时
|
||
cad.save_as(f"egm{phase_conductor_foo + 1}.dxf")
|
||
min_distance_intersection = (
|
||
np.sum(
|
||
(
|
||
np.array(rs_rc_circle_intersection)
|
||
- np.array(circle_rc_or_rg_line_intersection)
|
||
)
|
||
** 2
|
||
)
|
||
** 0.5
|
||
) # 计算两圆交点和地面直线交点的最小距离
|
||
if min_distance_intersection < 0.1:
|
||
break # 已经找到了最大电流
|
||
# 判断是否以完全被保护
|
||
if (
|
||
rs_rc_circle_intersection[1]
|
||
< circle_rc_or_rg_line_intersection[1]
|
||
):
|
||
circle_rs_line_or_rg_intersection = None
|
||
if rg_type == "g":
|
||
circle_rs_line_or_rg_intersection = (
|
||
solve_circle_line_intersection(
|
||
rs, rs_x, rs_y, rg_line_func
|
||
) # 保护弧和捕雷弧的交点
|
||
)
|
||
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:
|
||
logger.info(f"电流为{i_bar}kV时，暴露弧已经完全被屏蔽")
|
||
exposed_curve_shielded = True
|
||
break
|
||
animate.pause()
|
||
# 判断是否导线已经被完全保护
|
||
if abs(i_max - _max_i) < 1e-5:
|
||
logger.info("无法找到最大电流，可能是杆塔较高。")
|
||
logger.info(f"最大电流设置为自然界最大电流{i_max}kA")
|
||
logger.info(f"最大电流为{i_max:.2f}")
|
||
logger.info(f"最小电流为{i_min:.2f}")
|
||
# if exposed_curve_shielded:
|
||
# logger.info("暴露弧已经完全被屏蔽,不会跳闸。")
|
||
# continue
|
||
curt_fineness = 0.1 # 电流积分细度
|
||
if i_min > i_max or abs(i_min - i_max) < curt_fineness:
|
||
logger.info("最大电流小于等于最小电流，没有暴露弧。")
|
||
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)
|
||
td = para.td
|
||
ip_a = para.Ip_a
|
||
ip_b = para.Ip_b
|
||
bd_area_vec_result = bd_area_vec(
|
||
i_curt_samples,
|
||
u_ph,
|
||
rc_x,
|
||
rc_y,
|
||
rs_x,
|
||
rs_y,
|
||
rg_x,
|
||
rg_y,
|
||
ground_angel,
|
||
rg_type,
|
||
)
|
||
thunder_density_result = thunder_density(
|
||
i_curt_samples, td, ip_a, ip_b
|
||
) # 雷电流幅值密度函数
|
||
cal_bd_np = bd_area_vec_result * thunder_density_result
|
||
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
|
||
rated_voltage = para.rated_voltage
|
||
n_sf = (
|
||
2
|
||
* ng
|
||
/ 10
|
||
* calculus
|
||
* arc_possibility(rated_voltage, insulator_c_len)
|
||
)
|
||
avr_n_sf += n_sf / voltage_n
|
||
n_sf_phases[phase_conductor_foo][u_bar] = n_sf
|
||
logger.info(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.16f}次/(100km·a)")
|
||
logger.info(f"线路跳闸率是{avr_n_sf:.16f}次/(100km·a)")
|
||
logger.info(
|
||
f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=16)}次/(100km·a)"
|
||
)
|
||
|
||
|
||
def speed():
|
||
a = 0
|
||
for bar in range(100000000):
|
||
a += bar
|
||
|
||
|
||
if __name__ == "__main__":
|
||
logger.remove()
|
||
logger.add(sys.stderr, level="DEBUG")
|
||
egm()
|
||
# run_time = timeit.timeit("egm()", globals=globals(), number=1)
|
||
# logger.info(f"运行时间:{run_time:.2f}s")
|
||
# input('enter any key to exit.')
|
||
logger.info("Finished.")
|