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

参数全部从外部读取

This commit is contained in:
n3040 2021-12-22 16:11:14 +08:00
parent 2251966b7e
commit 7f03fc2b9c
4 changed files with 450 additions and 238 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

163
core.py
View File

@ -7,6 +7,38 @@ gMSP = None
gCount = 1 gCount = 1
class Parameter:
h_g_sag: float # 地线弧垂
h_c_sag: float # 导线弧垂
h_whole: float # 杆塔全高
voltage_n: int # 工作电压分成多少份来计算
td: int # 雷暴日
insulator_c_len: float # 串子绝缘长度
string_c_len: float
string_g_len: float
gc_x: [float] # 导、地线水平坐标
ground_angels: [float] # 地面倾角，向下为正
h_arm: float # 导、地线垂直坐标
altitude: int # 海拔，单位米
max_i: float # 最大尝试电流，单位kA
para = Parameter()
def rg_line_function_factory(_rg, ground_angel): # 返回一个地面捕雷线的直线方程
y_d = _rg / math.cos(ground_angel) # y轴上的截距
# 利用公式y-y0=k(x-x0) 得到直线公式
y0 = y_d
x0 = 0
k = math.tan(math.pi - ground_angel)
def f(x):
return y0 + k * (x - x0)
return f
class Draw: class Draw:
def __init__(self): def __init__(self):
self._doc = ezdxf.new(dxfversion="R2010") self._doc = ezdxf.new(dxfversion="R2010")
@ -14,7 +46,20 @@ class Draw:
global gCAD global gCAD
gCAD = self gCAD = self
def draw(self, i_curt, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, color): def draw(
self,
i_curt,
u_ph,
rs_x,
rs_y,
rc_x,
rc_y,
rg_x,
rg_y,
rg_type,
ground_angel,
color,
):
doc = self._doc doc = self._doc
msp = doc.modelspace() msp = doc.modelspace()
global gMSP global gMSP
@ -33,8 +78,15 @@ class Draw:
(rc_x, rc_y), circle_intersection, dxfattribs={"color": color} (rc_x, rc_y), circle_intersection, dxfattribs={"color": color}
) # 地线 ) # 地线
if rg_type == "g": if rg_type == "g":
msp.add_line((0, rg), (2000, rg), dxfattribs={"color": color}) ground_angel_func = rg_line_function_factory(rg, ground_angel)
circle_line_section = solve_circle_line_intersection(rc, rg, rc_x, rc_y) msp.add_line(
(0, ground_angel_func(0)),
(2000, ground_angel_func(2000)),
dxfattribs={"color": color},
)
circle_line_section = solve_circle_line_intersection(
rc, rc_x, rc_y, ground_angel_func
)
if not circle_line_section: if not circle_line_section:
pass pass
else: else:
@ -94,40 +146,77 @@ def solve_circle_intersection(
# 圆与捕雷线交点 # 圆与捕雷线交点
def solve_circle_line_intersection(radius, rg, center_x, center_y): def solve_circle_line_intersection(
distance = distance_point_line(center_x, center_y, 0, rg, 0) # 捕雷线到暴露圆中点的距离 radius, center_x, center_y, ground_surface_func
): # 返回交点的x和y坐标
x0 = 0
y0 = ground_surface_func(x0)
x1 = 1
y1 = ground_surface_func(x1)
k = (y1 - y0) / (x1 - x0)
distance = distance_point_line(center_x, center_y, x0, y0, k) # 捕雷线到暴露圆中点的距离
if distance > radius: if distance > radius:
return [] return []
else: else:
r = (radius ** 2 - (rg - center_y) ** 2) ** 0.5 + center_x # r = (radius ** 2 - (rg - center_y) ** 2) ** 0.5 + center_x
return [r, rg] a = center_x
b = center_y
c = y0
d = x0
bb = -2 * a + 2 * c * k - 2 * d * (k ** 2) - 2 * b * k
aa = 1 + k ** 2
rr = radius
cc = (
a ** 2
+ c ** 2
- 2 * c * k * d
+ (k ** 2) * (d ** 2)
- 2 * b * (c - k * d)
+ b ** 2
- rr ** 2
)
_x = (-bb + (bb ** 2 - 4 * aa * cc) ** 0.5) / 2 / aa
_y = ground_surface_func(_x)
# 验算结果
equ = (center_x - _x) ** 2 + (center_y - _y) ** 2 - radius ** 2
assert abs(equ) < 1e-5
return [_x, _y]
def min_i(string_len, u_ph): def min_i(string_len, u_ph):
u_50 = 530 * string_len + 35 # 海拔修正
altitude = para.altitude
k_a = math.exp(altitude / 8150) # 气隙海拔修正
u_50 = 1 / k_a * (530 * string_len + 35) # 50045 上附录的公式，实际应该用负极性电压的公式
z_0 = 300 # 雷电波阻抗 z_0 = 300 # 雷电波阻抗
z_c = 251 # 导线波阻抗 z_c = 251 # 导线波阻抗
# 新版大手册公式 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)
return r return r
def thunder_density(i): # l雷电流幅值密度函数 def thunder_density(i): # l雷电流幅值密度函数
r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44) td = para.td
r = None
if td == 20:
r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44) # 雷暴日20d
if td == 40:
r = -(10 ** (-i / 88)) * math.log(10) * (-1 / 88) # 雷暴日40d
return r return r
def angel_density(angle): # 入射角密度函数 angle单位是弧度 def angel_density(angle): # 入射角密度函数 angle单位是弧度
r = 0.75 * abs((np.cos(angle - math.pi / 2) ** 3)) r = 0.75 * abs((np.cos(angle - math.pi / 2) ** 3)) # 新版大手册公式3-275
return r return r
def rs_fun(i): def rs_fun(i):
r = 10 * (i ** 0.65) r = 10 * (i ** 0.65) # 新版大手册公式3-271
return r return r
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) r = 1.63 * ((5.015 * (i ** 0.578) - 0.001 * u_ph) ** 1.125) # 新版大手册公式3-272
return r return r
@ -136,16 +225,16 @@ 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) 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) rg = 5.5 * (i_curt ** 0.65) # 新版大手册公式3-273
elif typ == "c": # 此时返回的是圆半径 elif typ == "c": # 此时返回的是圆半径
rg = rc_fun(i_curt, u_ph) rg = rc_fun(i_curt, u_ph)
return rg return rg
def intersection_angle( def intersection_angle(
rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_surface, rg_type rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_angel, rg_type
): # 暴露弧的角度 ): # 暴露弧的角度
rs = rs_fun(i_curt) rs = rs_fun(i_curt)
rc = rc_fun(i_curt, u_ph) rc = rc_fun(i_curt, u_ph)
@ -154,35 +243,38 @@ def intersection_angle(
rs, rc, rs_x, rs_y, rc_x, rc_y rs, rc, rs_x, rs_y, rc_x, rc_y
) # 两圆的交点 ) # 两圆的交点
circle_line_or_rg_intersection = None circle_line_or_rg_intersection = None
rg_line_func = rg_line_function_factory(rg, ground_angel)
if rg_type == "g": if rg_type == "g":
circle_line_or_rg_intersection = solve_circle_line_intersection( circle_line_or_rg_intersection = solve_circle_line_intersection(
rc, rg, rc_x, rc_y rc, rc_x, rc_y, rg_line_func
) # 暴露圆和补雷线的交点 ) # 暴露圆和补雷线的交点
if rg_type == "c": if rg_type == "c":
circle_line_or_rg_intersection = solve_circle_intersection( circle_line_or_rg_intersection = solve_circle_intersection(
rg, rc, rg_x, rg_y, rc_x, rc_y rg, rc, rg_x, rg_y, rc_x, rc_y
) # 两圆的交点 ) # 两圆的交点
if circle_line_or_rg_intersection: # TODO 应该是不存在落到地面线以下的情况，先把以下注释掉
(
circle_line_or_rg_intersection_x, # if circle_line_or_rg_intersection:
circle_line_or_rg_intersection_y, # (
) = circle_line_or_rg_intersection # circle_line_or_rg_intersection_x,
if ( # circle_line_or_rg_intersection_y,
ground_surface(circle_line_or_rg_intersection_x) # ) = circle_line_or_rg_intersection
> circle_line_or_rg_intersection_y # if (
): # 交点在地面线以下，就可以不积分 # ground_surface(rg, 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 # # 找到暴露弧和地面线的交点
) # circle_line_or_rg_intersection = circle_ground_surface_intersection(
# rc, rc_x, rc_y, ground_surface
# )
theta1 = None theta1 = None
np_circle_intersection = np.array(circle_intersection) np_circle_intersection = np.array(circle_intersection)
theta2_line = np_circle_intersection - np.array([rc_x, rc_y]) theta2_line = np_circle_intersection - np.array([rc_x, rc_y])
theta2 = math.atan(theta2_line[1] / theta2_line[0]) theta2 = math.atan(theta2_line[1] / theta2_line[0])
np_circle_line_or_rg_intersection = np.array(circle_line_or_rg_intersection) np_circle_line_or_rg_intersection = np.array(circle_line_or_rg_intersection)
if not circle_line_or_rg_intersection: if not circle_line_or_rg_intersection:
if rc_y - rc > rg: # rg在rc下面 if rc_y - rc > rg_line_func(rc_x): # rg在rc下面
# 捕捉线太低了，对高塔无保护，θ_1从-90°开始计算。 # 捕捉线太低了，对高塔无保护，θ_1从-90°开始计算，即从与地面垂直的角度开始就已经暴露了。
theta1 = -math.pi / 2 theta1 = -math.pi / 2
else: else:
theta1_line = np_circle_line_or_rg_intersection - np.array([rc_x, rc_y]) theta1_line = np_circle_line_or_rg_intersection - np.array([rc_x, rc_y])
@ -202,6 +294,7 @@ def func_calculus_pw(theta, max_w):
if segments < 2: # 最大最小太小，没有可以积分的 if segments < 2: # 最大最小太小，没有可以积分的
return 0 return 0
w_samples, d_w = np.linspace(0, max_w, segments, retstep=True) w_samples, d_w = np.linspace(0, max_w, segments, retstep=True)
# 童中宇 750KV信洛Ⅰ线雷电防护性能研究 公式 3-10
cal_w_np = abs(angel_density(w_samples)) * np.sin(theta - (w_samples - math.pi / 2)) cal_w_np = abs(angel_density(w_samples)) * np.sin(theta - (w_samples - math.pi / 2))
r_pw = np.sum((cal_w_np[:-1] + cal_w_np[1:])) / 2 * d_w r_pw = np.sum((cal_w_np[:-1] + cal_w_np[1:])) / 2 * d_w
return r_pw return r_pw
@ -216,7 +309,7 @@ def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y): # 对θ进行积分
# 求保护弧到直线的距离，判断是否相交 # 求保护弧到直线的距离，判断是否相交
d_to_rs = distance_point_line(rs_x, rs_y, line_x, line_y, k) d_to_rs = distance_point_line(rs_x, rs_y, line_x, line_y, k)
if d_to_rs < rs: # 相交 if d_to_rs < rs: # 相交
# 要用过直线上一点到暴露弧的切线 # 要用过这一点到保护弧的切线
new_k = tangent_line_k(line_x, line_y, rs_x, rs_y, rs, init_k=k) new_k = tangent_line_k(line_x, line_y, rs_x, rs_y, rs, init_k=k)
if new_k >= 0: if new_k >= 0:
max_w = math.atan(new_k) # 用于保护弧相切的角度 max_w = math.atan(new_k) # 用于保护弧相切的角度
@ -249,15 +342,16 @@ def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y): # 对θ进行积分
# ) # )
# gCAD.save() # gCAD.save()
pass pass
# 童中宇 750KV信洛Ⅰ线雷电防护性能研究 公式 3-10
r = rc / math.cos(theta) * func_calculus_pw(theta, max_w) r = rc / math.cos(theta) * func_calculus_pw(theta, max_w)
return r return r
def bd_area( def bd_area(
i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, ground_surface, rg_type i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, ground_angel, rg_type
): # 暴露弧的投影面积 ): # 暴露弧的投影面积
theta1, theta2 = intersection_angle( theta1, theta2 = intersection_angle(
rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_surface, rg_type rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_angel, rg_type
) # θ角度 ) # θ角度
theta_fineness = 0.01 theta_fineness = 0.01
rc = rc_fun(i_curt, u_ph) rc = rc_fun(i_curt, u_ph)
@ -364,6 +458,7 @@ def circle_ground_surface_intersection(radius, center_x, center_y, ground_surfac
# u_ph是相电压 # u_ph是相电压
# insulator_c_len绝缘子闪络距离 # insulator_c_len绝缘子闪络距离
def arc_possibility(rated_voltage, insulator_c_len): # 建弧率 def arc_possibility(rated_voltage, insulator_c_len): # 建弧率
# 50064 中附录给的公式
_e = rated_voltage / (3 ** 0.5) / insulator_c_len _e = rated_voltage / (3 ** 0.5) / 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

468
main.py
View File

@ -1,225 +1,309 @@
import math
import sys import sys
import tomli
from loguru import logger from loguru import logger
from core import * from core import *
import timeit import timeit
def egm(): def egm():
h_g_avr_sag = 11.67 * 2 / 3 if len(sys.argv) < 2:
h_c_avr_sag = 14.43 * 2 / 3 toml_file_path = "default.toml"
h_whole = 130 # 杆塔全高 # # logger.info('没指定计算文件!程序结束。')
voltage_n = 3 # 工作电压分成多少份来计算 # # sys.exit(0)
td = 20 # 雷暴日 # h_g_sag = 20 # 地线弧垂
insulator_c_len = 6.6 # 串子绝缘长度 # h_c_sag = 20 # 导线弧垂
string_c_len = 9.2 # h_whole = 106.1 # 杆塔全高
string_g_len = 0.5 # voltage_n = 3 # 工作电压分成多少份来计算
gc_x = [17.9, 17, 15, 17.0] # td = 20 # 雷暴日
# insulator_c_len = 6.98 # 串子绝缘长度
# 以后考虑地形角度,地面线 # string_c_len = 9.2
def ground_surface(x): # string_g_len = 0.63
return 0 # gc_x = [32.2 / 2, 32.2 / 2, 15, 17.0] # 导、地线水平坐标
# ground_angels = [40 / 180 * math.pi] # 地面倾角，向下为正
# h_arm = 34
# gc_y = [
# h_whole - string_g_len - h_g_sag * 2 / 3, # 地线对地平均高
# # h_whole - string_c_len - h_c_sag - 2.7, # 导线对地平均高
# h_whole - string_c_len - h_c_sag * 2 / 3 - h_arm, # 导线对地平均高
# # h_whole - string_c_len - h_c_sag - 35.7, # 导线对地平均高
# ]
else:
toml_file_path = sys.argv[1]
logger.info(f"读取文件{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.max_i = toml_parameter["max_i"]
toml_optional = toml_dict["optional"]
para.voltage_n = toml_optional["voltage_n"] # 工作电压分成多少份来计算
#########################################################
# 以上是需要设置的参数
h_whole = para.h_whole
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 = [ gc_y = [
h_whole - string_g_len - h_g_avr_sag, # 地线对地平均高 h_whole - string_g_len - h_g_sag * 2 / 3, # 地线对地平均高
# 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(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: # 双回路 if len(gc_y) > 2: # 双回路
phase_n = 3 # 边相导线数量 phase_n = 3 # 边相导线数量
else: else:
phase_n = 1 phase_n = 1
######################################################### # 地闪密度 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
rg_type = None td = para.td
# 以上是需要设置的参数
avr_n_sf = 0 # 考虑电压的影响计算的跳闸率
rg_x = None
rg_y = None
cad = Draw()
# 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
ng = func_ng(td) ng = func_ng(td)
n_sf_phases = np.zeros((phase_n, voltage_n)) # 计算每一相的跳闸率 avr_n_sf = 0 # 考虑电压的影响计算的跳闸率
if np.any(np.array(gc_y) < 0): ground_angels = para.ground_angels
logger.info("导线可能掉地面了，程序退出。") for ground_angel in ground_angels:
return 0 logger.info(f"地面倾角{ground_angel/math.pi*180:.3f}°")
for phase_conductor_foo in range(phase_n): rg_type = None
exposed_curve_shielded = False rg_x = None
rs_x = gc_x[phase_conductor_foo] rg_y = None
rs_y = gc_y[phase_conductor_foo] cad = Draw()
rc_x = gc_x[phase_conductor_foo + 1] voltage_n = para.voltage_n
rc_y = gc_y[phase_conductor_foo + 1] n_sf_phases = np.zeros((phase_n, voltage_n)) # 存储每一相的跳闸率
if phase_n == 1: if np.any(np.array(gc_y) < 0):
rg_type = "g" logger.info("导线可能掉地面下了，程序退出。")
if phase_n > 1: # 多回路 return 0
if phase_conductor_foo < 2: for phase_conductor_foo in range(phase_n):
rg_type = "c" exposed_curve_shielded = False
rg_x = gc_x[phase_conductor_foo + 2] rs_x = gc_x[phase_conductor_foo]
rg_y = gc_y[phase_conductor_foo + 2] rs_y = gc_y[phase_conductor_foo]
else: rc_x = gc_x[phase_conductor_foo + 1]
rc_y = gc_y[phase_conductor_foo + 1]
if phase_n == 1:
rg_type = "g" rg_type = "g"
# TODO 保护角公式可能有问题，后面改 if phase_n > 1: # 多回路
shield_angle = ( if phase_conductor_foo < 2:
math.atan((rc_x - rs_x) / ((rs_y - rc_y) + string_c_len)) * 180 / math.pi rg_type = "c" # 捕捉弧有下面一相导线的击距代替
) # 保护角 rg_x = gc_x[phase_conductor_foo + 2]
logger.info(f"保护角{shield_angle:.3f}°") rg_y = gc_y[phase_conductor_foo + 2]
logger.debug(f"最低相防护标识{rg_type}") else:
for u_bar in range(voltage_n): rg_type = "g"
u_ph = ( # TODO 保护角公式可能有问题，后面改
-math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732 shield_angle = (
) # 运行相电压 math.atan((rc_x - rs_x) / ((rs_y - rc_y) + string_c_len))
logger.info(f"计算第{phase_conductor_foo + 1}相，电压为{u_ph:.2f}kV") * 180
# 迭代法计算最大电流 / math.pi
i_max = 0 ) # 保护角
i_min = min_i(insulator_c_len, u_ph / 1.732) logger.info(f"保护角{shield_angle:.3f}°")
_min_i = i_min # 尝试的最小电流 logger.debug(f"最低相防护标识{rg_type}")
_max_i = 200 # 尝试的最大电流 for u_bar in range(voltage_n): # 计算不同工作电压下的跳闸率
# cad.draw(i_min, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 2) u_ph = (
for i_bar in np.linspace( math.sqrt(2)
_min_i, _max_i, int((_max_i - _min_i) / 0.1) * 750
): # 雷电流 * math.cos(2 * math.pi / voltage_n * u_bar)
# logger.info(f"尝试计算电流为{i_bar:.2f}") / 1.732
rs = rs_fun(i_bar) ) # 运行相电压
rc = rc_fun(i_bar, u_ph) logger.info(f"计算第{phase_conductor_foo + 1}相，电压为{u_ph:.2f}kV")
rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type) # 迭代法计算最大电流
####### i_max = 0
# cccCount += 1 insulator_c_len = para.insulator_c_len
# if cccCount % 30 == 0: i_min = min_i(insulator_c_len, u_ph / 1.732)
# import core _min_i = i_min # 尝试的最小电流
# _max_i = para.max_i # 尝试的最大电流
# core.gMSP.add_circle((0, h_gav), rs) # cad.draw(i_min, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 2)
# core.gMSP.add_circle( for i_bar in np.linspace(
# (dgc, h_cav), rc_fun(i_bar, -u_ph), dxfattribs={"color": 4} _min_i, _max_i, int((_max_i - _min_i) / 0.1)
# ) ): # 雷电流
# core.gMSP.add_circle((dgc, h_cav), rc) # logger.info(f"尝试计算电流为{i_bar:.2f}")
####### rs = rs_fun(i_bar)
rg_rc_circle_intersection = solve_circle_intersection( rc = rc_fun(i_bar, u_ph)
rs, rc, rs_x, rs_y, rc_x, rc_y rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type)
) rg_line_func = None
i_max = i_bar
if not rg_rc_circle_intersection: # if circle_intersection is []
logger.debug("保护弧和暴露弧无交点，检查设置参数。")
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
)
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:
# 暴露弧和捕捉弧无交点
if rg_type == "g": if rg_type == "g":
if rg > rc_y: rg_line_func = rg_line_function_factory(rg, ground_angel)
i_min = i_bar #######
logger.info(f"捕捉弧在暴露弧之上，设置最小电流为{i_min:.2f}") # 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
)
i_max = i_bar
if not rg_rc_circle_intersection: # if circle_intersection is []
logger.debug("保护弧和暴露弧无交点，检查设置参数。")
continue
circle_rc_line_or_rg_intersection = None
if rg_type == "g":
circle_rc_line_or_rg_intersection = (
solve_circle_line_intersection(rc, rc_x, rc_y, rg_line_func)
)
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:
# 暴露弧和捕捉弧无交点
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: else:
logger.info("暴露弧和捕捉弧无交点，检查设置参数。") logger.info("上面的导地线无法保护下面的导地线，检查设置参数。")
continue continue
else: min_distance_intersection = (
logger.info("暴露弧和捕捉弧无交点，检查设置参数。")
continue
min_distance_intersection = (
np.sum(
(
np.array(rg_rc_circle_intersection)
- np.array(circle_rc_line_or_rg_intersection)
)
** 2
)
** 0.5
) # 计算两圆交点和地面直线交点的最小距离
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
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.sum(
( (
np.array(circle_rs_line_or_rg_intersection) np.array(rg_rc_circle_intersection)
- np.array([rc_x, rc_y]) - np.array(circle_rc_line_or_rg_intersection)
) )
** 2 ** 2
) )
** 0.5 ** 0.5
) ) # 计算两圆交点和地面直线交点的最小距离
if distance > rc: if min_distance_intersection < 0.1:
logger.info("暴露弧已经完全被屏蔽") break # 已经找到了最大电流
exposed_curve_shielded = True # 判断是否以完全被保护
break if (
# if phase_conductor_foo == 2: rg_rc_circle_intersection[1]
cad.draw(i_min, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 2) < circle_rc_line_or_rg_intersection[1]
cad.draw(i_max, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 6) ):
cad.save_as(f"egm{phase_conductor_foo + 1}.dxf") circle_rs_line_or_rg_intersection = None
# 判断是否导线已经被完全保护 if rg_type == "g":
if abs(i_max - _max_i) < 1e-5: circle_rs_line_or_rg_intersection = (
logger.info("无法找到最大电流，可能是杆塔较高。") solve_circle_line_intersection(
logger.info(f"最大电流设置为自然界最大电流{i_max}kA") rs, rs_x, rs_y, rg_line_func
logger.info(f"最大电流为{i_max:.2f}") ) # 保护弧和捕雷弧的交点
logger.info(f"最小电流为{i_min:.2f}") )
if exposed_curve_shielded: if rg_type == "c":
logger.info("暴露弧已经完全被屏蔽,不会跳闸。") circle_rs_line_or_rg_intersection = (
continue solve_circle_intersection(
curt_fineness = 0.1 # 电流积分细度 rs, rg, rs_x, rs_y, rg_x, rg_y
if i_min > i_max or abs(i_min - i_max) < curt_fineness: )
logger.info("最大电流小于最小电流，没有暴露弧。") )
continue # 判断与保护弧的交点是否在暴露弧外面
# 开始积分 distance = (
curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份 np.sum(
i_curt_samples, d_curt = np.linspace( (
i_min, i_max, curt_segment_n + 1, retstep=True np.array(circle_rs_line_or_rg_intersection)
) - np.array([rc_x, rc_y])
bd_area_vec = np.vectorize(bd_area) )
cal_bd_np = ( ** 2
bd_area_vec( )
i_curt_samples, ** 0.5
)
if distance > rc:
logger.info("暴露弧已经完全被屏蔽")
exposed_curve_shielded = True
break
# if phase_conductor_foo == 2:
cad.draw(
i_min,
u_ph, u_ph,
rc_x,
rc_y,
rs_x, rs_x,
rs_y, rs_y,
rc_x,
rc_y,
rg_x, rg_x,
rg_y, rg_y,
ground_surface,
rg_type, rg_type,
ground_angel,
2,
) )
* thunder_density(i_curt_samples) cad.draw(
) i_max,
calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt u_ph,
# for i_curt in i_curt_samples[:-1]: rs_x,
# cal_bd_first = bd_area(i_curt, u_ph, dgc, h_gav, h_cav) rs_y,
# cal_bd_second = bd_area(i_curt + d_curt, u_ph, dgc, h_gav, h_cav) rc_x,
# cal_thunder_density_first = thunder_density(i_curt) rc_y,
# cal_thunder_density_second = thunder_density(i_curt + d_curt) rg_x,
# calculus += ( rg_y,
# ( rg_type,
# cal_bd_first * cal_thunder_density_first ground_angel,
# + cal_bd_second * cal_thunder_density_second 6,
# ) )
# / 2 cad.save_as(f"egm{phase_conductor_foo + 1}.dxf")
# * d_curt # 判断是否导线已经被完全保护
# ) if abs(i_max - _max_i) < 1e-5:
# if abs(calculus-0.05812740052770032)<1e-5: logger.info("无法找到最大电流，可能是杆塔较高。")
# abc=123 logger.info(f"最大电流设置为自然界最大电流{i_max}kA")
# pass logger.info(f"最大电流为{i_max:.2f}")
n_sf = 2 * ng / 10 * calculus * arc_possibility(750, insulator_c_len) logger.info(f"最小电流为{i_min:.2f}")
avr_n_sf += n_sf / voltage_n if exposed_curve_shielded:
n_sf_phases[phase_conductor_foo][u_bar] = n_sf logger.info("暴露弧已经完全被屏蔽,不会跳闸。")
logger.info(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}") continue
logger.info(f"跳闸率是{avr_n_sf:.6f}") curt_fineness = 0.1 # 电流积分细度
logger.info(f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=6)}") 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)
cal_bd_np = (
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(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 * arc_possibility(750, 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:.6}")
logger.info(f"跳闸率是{avr_n_sf:.6f}")
logger.info(
f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=6)}"
)
def speed(): def speed():

24
plot_data.py
View File

@ -40,24 +40,24 @@ data_130m塔高_不同接地电阻 = np.array(
] ]
) )
category_names_130m塔高_不同地线保护角 = ["-1", "-3", "-6"] category_names_150m塔高_不同地线保护角 = ["-1", "-3", "-6"]
data_130m塔高_不同地线保护角 = np.array( data_150m塔高_不同地线保护角 = np.array(
[ [
[0.000170, 0.103132], [0.000440, 0.155414],
[0.000168, 0.079659], [0.000433, 0.128227],
[0.000167, 0.055598], [0.000429, 0.099819],
] ]
) )
category_names_66m串长_不同塔高 = ['100', '120', '140'] category_names_66m串长_不同塔高 = ["100", "120", "140"]
data_66m串长_不同塔高 = np.array( data_66m串长_不同塔高 = np.array(
[ [
[0.000053 , 0.023285], [0.000053, 0.023285],
[0.000139 , 0.083229], [0.000139, 0.083229],
[0.000470 , 0.145586], [0.000470, 0.145586],
] ]
) )
@ -65,8 +65,8 @@ category_names_68m串长_不同塔高 = [100, 120, 140]
data_68m串长_不同塔高 = np.array( data_68m串长_不同塔高 = np.array(
[ [
[0.000039 , 0.019094], [0.000039, 0.019094],
[0.000098 , 0.073033], [0.000098, 0.073033],
[0.000287 , 0.130923], [0.000287, 0.130923],
] ]
) )

33
server.py Normal file
View File

@ -0,0 +1,33 @@
from fastapi import FastAPI
import uvicorn
from pydantic import BaseModel
class CalculationParameter(BaseModel):
loop: str # single or double 回路数
insulator_c_len: float # 绝缘长度
string_c_len: float # 导线串长
string_g_len: float # 地线串长
td: int # 雷暴日
h_g_avr_sag: float # 地线平均弧垂
h_c_avr_sag: float # 导线平均弧垂
h_whole: float # 杆塔全高
gc_x: tuple[float] # 导、地线水平坐标
ground_angels: tuple[float] # 地面倾角，向下为正
fastapi_app = FastAPI()
@fastapi_app.post("/calculation")
async def calculation(cp: CalculationParameter):
print(cp)
return {"Hello": "World"}
def server_start():
pass
if __name__ == "__main__":
uvicorn.run("server:fastapi_app", host="127.0.0.1", port=5000, log_level="info")