Compare commits

..

2 Commits

Author SHA1 Message Date
n3040 5dc1613d2a 加了一些TODO 2024-01-03 14:28:28 +08:00
n3040 9b852235f1 准备进行陇东大跨越计算。 2022-07-15 12:47:43 +08:00
2 changed files with 47 additions and 33 deletions

26
core.py
View File

@ -1,6 +1,7 @@
import math import math
import ezdxf import ezdxf
import numpy as np import numpy as np
from typing import List
gCAD = None gCAD = None
gMSP = None gMSP = None
@ -15,8 +16,8 @@ class Parameter:
insulator_c_len: float # 串子绝缘长度 insulator_c_len: float # 串子绝缘长度
string_c_len: float string_c_len: float
string_g_len: float string_g_len: float
gc_x: [float] # 导、地线水平坐标 gc_x: List[float] # 导、地线水平坐标
ground_angels: [float] # 地面倾角,向下为正 ground_angels: List[float] # 地面倾角,向下为正
h_arm: float # 导、地线垂直坐标 h_arm: float # 导、地线垂直坐标
altitude: int # 海拔,单位米 altitude: int # 海拔,单位米
max_i: float # 最大尝试电流单位kA max_i: float # 最大尝试电流单位kA
@ -131,7 +132,7 @@ def solve_circle_intersection(
x = radius2 + center_x2 # 初始值 x = radius2 + center_x2 # 初始值
y = radius2 + center_y2 # 初始值 y = radius2 + center_y2 # 初始值
# TODO 考虑出现2个解的情况 # TODO 考虑出现2个解的情况
for bar in range(0, 10): for _ in range(0, 10):
A = [ A = [
[-2 * (x - center_x1), -2 * (y - center_y1)], [-2 * (x - center_x1), -2 * (y - center_y1)],
[-2 * (x - center_x2), -2 * (y - center_y2)], [-2 * (x - center_x2), -2 * (y - center_y2)],
@ -189,12 +190,20 @@ def solve_circle_line_intersection(
def min_i(string_len, u_ph): def min_i(string_len, u_ph):
# 海拔修正 # 海拔修正
altitude = para.altitude altitude = para.altitude
if altitude > 1000:
k_a = math.exp((altitude - 1000) / 8150) # 气隙海拔修正 k_a = math.exp((altitude - 1000) / 8150) # 气隙海拔修正
else:
k_a = 1
u_50 = 1 / k_a * (530 * string_len + 35) # 50045 上附录的公式,实际应该用负极性电压的公式 u_50 = 1 / k_a * (530 * string_len + 35) # 50045 上附录的公式,实际应该用负极性电压的公式
# u_50 = 1 / k_a * (533 * string_len + 132) # 串放电路径 1000m海拔
# u_50 = 1 / k_a * (477 * string_len + 99) # 串放电路径 2000m海拔
# u_50 = 615 * string_len # 导线对塔身放电 1000m海拔
# u_50= 263.32647401+533.90081562*string_len
z_0 = 300 # 雷电波阻抗 z_0 = 300 # 雷电波阻抗
z_c = 251 # 导线波阻抗 z_c = 251 # 导线波阻抗
# 新版大手册公式 3-277 # 新版大手册公式 3-277
r = (u_50 + 2 * z_0 / (2 * z_0 + z_c) * u_ph) * (2 * z_0 + z_c) / (z_0 * z_c) r = (u_50 + 2 * z_0 / (2 * z_0 + z_c) * u_ph) * (2 * z_0 + z_c) / (z_0 * z_c)
# r = 2 * (u_50 - u_ph) / z_c
return r return r
@ -207,12 +216,15 @@ def thunder_density(i, td, ip_a, ip_b): # 雷电流幅值密度函数
r = -( r = -(
-ip_b / ip_a / ((1 + (i / ip_a) ** ip_b) ** 2) * ((i / ip_a) ** (ip_b - 1)) -ip_b / ip_a / ((1 + (i / ip_a) ** ip_b) ** 2) * ((i / ip_a) ** (ip_b - 1))
) )
return r
else: else:
if td == 20: if td == 20:
r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44) # 雷暴日20d r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44) # 雷暴日20d
return r
if td == 40: if td == 40:
r = -(10 ** (-i / 88)) * math.log(10) * (-1 / 88) # 雷暴日40d r = -(10 ** (-i / 88)) * math.log(10) * (-1 / 88) # 雷暴日40d
return r return r
raise Exception("检查雷电参数!")
def angel_density(angle): # 入射角密度函数 angle单位是弧度 def angel_density(angle): # 入射角密度函数 angle单位是弧度
@ -226,7 +238,7 @@ def rs_fun(i):
def rc_fun(i, u_ph): def rc_fun(i, u_ph):
r = 1.63 * ((5.015 * (i ** 0.578) - 0.001 * u_ph) ** 1.125) # 新版大手册公式3-272 r = 1.63 * ((5.015 * (i**0.578) - 0.001 * u_ph * 1) ** 1.125) # 新版大手册公式3-272
return r return r
@ -235,7 +247,7 @@ def rg_fun(i_curt, h_cav, u_ph, typ="g"):
rg = None rg = None
if typ == "g": if typ == "g":
if h_cav < 40: if h_cav < 40:
rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65) # 新版大手册公式3-273 rg = (3.6 + 1.7 * math.log(43 - h_cav)) * (i_curt**0.65) # 新版大手册公式3-273
else: else:
rg = 5.5 * (i_curt**0.65) # 新版大手册公式3-273 rg = 5.5 * (i_curt**0.65) # 新版大手册公式3-273
elif typ == "c": # 此时返回的是圆半径 elif typ == "c": # 此时返回的是圆半径
@ -473,6 +485,8 @@ def circle_ground_surface_intersection(radius, center_x, center_y, ground_surfac
# insulator_c_len绝缘子闪络距离 # insulator_c_len绝缘子闪络距离
def arc_possibility(rated_voltage, insulator_c_len): # 建弧率 def arc_possibility(rated_voltage, insulator_c_len): # 建弧率
# 50064 中附录给的公式 # 50064 中附录给的公式
_e = rated_voltage / (3 ** 0.5) / insulator_c_len # TODO 需要区分交直流
# _e = rated_voltage / (3**0.5) / insulator_c_len #交流
_e = abs(rated_voltage) / (1) / insulator_c_len # 直流
r = (4.5 * (_e**0.75) - 14) * 1e-2 r = (4.5 * (_e**0.75) - 14) * 1e-2
return r return r

View File

@ -59,7 +59,7 @@ def read_parameter(toml_file_path):
def egm(): def egm():
if len(sys.argv) < 2: if len(sys.argv) < 2:
toml_file_path = r"default.toml" toml_file_path = r"article.toml"
else: else:
toml_file_path = sys.argv[1] toml_file_path = sys.argv[1]
if not os.path.exists(toml_file_path): if not os.path.exists(toml_file_path):
@ -333,10 +333,10 @@ def egm():
) )
avr_n_sf += n_sf / voltage_n avr_n_sf += n_sf / voltage_n
n_sf_phases[phase_conductor_foo][u_bar] = n_sf n_sf_phases[phase_conductor_foo][u_bar] = n_sf
logger.info(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.16f}") logger.info(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.16f}次/(km·a)")
logger.info(f"跳闸率是{avr_n_sf:.16f}") logger.info(f"线路跳闸率是{avr_n_sf:.16f}次/(km·a)")
logger.info( logger.info(
f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=16)}" f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=16)}次/(km·a)"
) )