Compare commits
No commits in common. "392eeb01688b5c3bdefa06b046066fb8ffbdf1e5" and "5a75df45429cf8bc644e3a6847ffda2320ee8988" have entirely different histories.
392eeb0168
...
5a75df4542
173
core.py
173
core.py
|
|
@ -15,44 +15,26 @@ 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, h_gav, h_cav, dgc, color):
|
||||||
doc = self._doc
|
doc = self._doc
|
||||||
msp = doc.modelspace()
|
msp = doc.modelspace()
|
||||||
global gMSP
|
global gMSP
|
||||||
gMSP = msp
|
gMSP = msp
|
||||||
rs = rs_fun(i_curt)
|
rs = rs_fun(i_curt)
|
||||||
rc = rc_fun(i_curt, u_ph)
|
rc = rc_fun(i_curt, u_ph)
|
||||||
rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
|
rg = rg_fun(i_curt, h_cav)
|
||||||
msp.add_circle((rs_x, rs_y), rs, dxfattribs={"color": color})
|
msp.add_circle((0, h_gav), rs, dxfattribs={"color": color})
|
||||||
msp.add_line((0, 0), (rs_x, rs_y)) # 地线
|
msp.add_line((0, 0), (0, h_gav)) # 地线
|
||||||
msp.add_circle((rc_x, rc_y), rc, dxfattribs={"color": color + 2})
|
msp.add_circle((dgc, h_cav), rc, dxfattribs={"color": color})
|
||||||
msp.add_line((rc_x, 0), (rc_x, rc_y)) # 导线
|
msp.add_line((dgc, 0), (dgc, h_cav)) # 导线
|
||||||
msp.add_line((rs_x, rs_y), (rc_x, rc_y))
|
msp.add_line((0, h_gav), (dgc, h_cav))
|
||||||
# 角度线
|
msp.add_line((0, rg), (2000, rg), dxfattribs={"color": color})
|
||||||
circle_intersection = solve_circle_intersection(rs, rc, rs_x, rs_y, rc_x, rc_y)
|
|
||||||
msp.add_line(
|
|
||||||
(rc_x, rc_y), circle_intersection, dxfattribs={"color": color}
|
|
||||||
) # 地线
|
|
||||||
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)
|
|
||||||
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(
|
|
||||||
rg, rc, rg_x, rg_y, rc_x, rc_y
|
|
||||||
)
|
|
||||||
msp.add_line(
|
|
||||||
(rc_x, rc_y), rg_rc_intersection, dxfattribs={"color": color}
|
|
||||||
) # 圆和圆的交点
|
|
||||||
# 计算圆交点
|
# 计算圆交点
|
||||||
|
# circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
|
||||||
|
# msp.add_line((0, h_gav), circle_intersection) # 地线
|
||||||
# msp.add_line((dgc, h_cav), circle_intersection) # 导线
|
# msp.add_line((dgc, h_cav), circle_intersection) # 导线
|
||||||
|
# circle_line_section = solve_circle_line_intersection(rc, rg, h_cav, dgc)
|
||||||
|
# msp.add_line((0, 0), circle_line_section) # 导线和圆的交点
|
||||||
|
|
||||||
def save(self):
|
def save(self):
|
||||||
doc = self._doc
|
doc = self._doc
|
||||||
|
|
@ -60,26 +42,15 @@ class Draw:
|
||||||
|
|
||||||
|
|
||||||
# 圆交点
|
# 圆交点
|
||||||
def solve_circle_intersection(
|
def solve_circle_intersection(rs, rc, h_gav, h_cav, dgc):
|
||||||
radius1,
|
|
||||||
radius2,
|
|
||||||
center_x1,
|
|
||||||
center_y1,
|
|
||||||
center_x2,
|
|
||||||
center_y2,
|
|
||||||
):
|
|
||||||
# 用牛顿法求解
|
# 用牛顿法求解
|
||||||
x = radius2 # 初始值
|
x = rc # 初始值
|
||||||
y = radius2 # 初始值
|
y = rc # 初始值
|
||||||
# TODO 考虑出现2个解的情况
|
|
||||||
for bar in range(0, 10):
|
for bar in range(0, 10):
|
||||||
A = [
|
A = [[-2 * x, -2 * (y - h_gav)], [-2 * (x - dgc), -2 * (y - h_cav)]]
|
||||||
[-2 * (x - center_x1), -2 * (y - center_y1)],
|
|
||||||
[-2 * (x - center_x2), -2 * (y - center_y2)],
|
|
||||||
]
|
|
||||||
b = [
|
b = [
|
||||||
(x - center_x1) ** 2 + (y - center_y1) ** 2 - radius1 ** 2,
|
x ** 2 + (y - h_gav) ** 2 - rs ** 2,
|
||||||
(x - center_x2) ** 2 + (y - center_y2) ** 2 - radius2 ** 2,
|
(x - dgc) ** 2 + (y - h_cav) ** 2 - rc ** 2,
|
||||||
]
|
]
|
||||||
X_set = np.linalg.solve(A, b)
|
X_set = np.linalg.solve(A, b)
|
||||||
x += X_set[0]
|
x += X_set[0]
|
||||||
|
|
@ -113,7 +84,7 @@ def thunder_density(i): # l雷电流幅值密度函数
|
||||||
|
|
||||||
|
|
||||||
def angel_density(angle): # 入射角密度函数 angle单位是弧度
|
def angel_density(angle): # 入射角密度函数 angle单位是弧度
|
||||||
r = 0.75 * abs((np.cos(angle - math.pi / 2) ** 3))
|
r = 0.75 * (np.cos(angle - math.pi / 2) ** 3)
|
||||||
return r
|
return r
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -124,62 +95,35 @@ 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)
|
r = 1.63 * ((5.015 * (i ** 0.578) - 0.001 * u_ph) ** 1.125)
|
||||||
|
# r=14.7*(i**0.42)
|
||||||
return r
|
return r
|
||||||
|
|
||||||
|
|
||||||
# typ 如果是g，代表捕雷线公式，c代表暴露弧公式
|
def rg_fun(i_curt, h_cav):
|
||||||
def rg_fun(i_curt, h_cav, u_ph, typ="g"):
|
if h_cav < 40:
|
||||||
rg = None
|
rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65)
|
||||||
if typ == "g":
|
else:
|
||||||
if h_cav < 40:
|
rg = 5.5 * (i_curt ** 0.65)
|
||||||
rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65)
|
|
||||||
else:
|
|
||||||
rg = 5.5 * (i_curt ** 0.65)
|
|
||||||
elif typ == "c": # 此时返回的是圆半径
|
|
||||||
rg = rc_fun(i_curt, u_ph)
|
|
||||||
return rg
|
return rg
|
||||||
|
|
||||||
|
|
||||||
def intersection_angle(
|
def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph): # 暴露弧的角度
|
||||||
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)
|
rs = rs_fun(i_curt)
|
||||||
rc = rc_fun(i_curt, u_ph)
|
rc = rc_fun(i_curt, u_ph)
|
||||||
rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
|
rg = rg_fun(i_curt, h_cav)
|
||||||
circle_intersection = solve_circle_intersection(
|
circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc) # 两圆的交点
|
||||||
rs, rc, rs_x, rs_y, rc_x, rc_y
|
circle_line_intersection = solve_circle_line_intersection(
|
||||||
) # 两圆的交点
|
rc, rg, dgc, h_cav
|
||||||
circle_line_or_rg_intersection = None
|
) # 暴露圆和补雷线的交点
|
||||||
if rg_type == "g":
|
|
||||||
circle_line_or_rg_intersection = solve_circle_line_intersection(
|
|
||||||
rc, rg, rc_x, rc_y
|
|
||||||
) # 暴露圆和补雷线的交点
|
|
||||||
if rg_type == "c":
|
|
||||||
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)
|
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([dgc, h_cav])
|
||||||
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_intersection = np.array(circle_line_intersection)
|
||||||
theta1_line = np_circle_line_or_rg_intersection - np.array([rc_x, rc_y])
|
theta1_line = np_circle_line_intersection - np.array([dgc, h_cav])
|
||||||
theta1 = math.atan(theta1_line[1] / theta1_line[0])
|
theta1 = math.atan(theta1_line[1] / theta1_line[0])
|
||||||
return np.array([theta1, theta2])
|
return np.array([theta1, theta2])
|
||||||
|
|
||||||
|
|
||||||
# 点到直线的距离
|
|
||||||
def distance_point_line(point_x, point_y, line_x, line_y, k) -> float:
|
def distance_point_line(point_x, point_y, line_x, line_y, k) -> float:
|
||||||
d = abs(k * point_x - point_y - k * line_x + line_y) / ((k ** 2 + 1) ** 0.5)
|
d = abs(k * point_x - point_y - k * line_x + line_y) / ((k ** 2 + 1) ** 0.5)
|
||||||
return d
|
return d
|
||||||
|
|
@ -193,17 +137,17 @@ def func_calculus_pw(theta, max_w):
|
||||||
return r_pw
|
return r_pw
|
||||||
|
|
||||||
|
|
||||||
def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y): # 对θ进行积分
|
def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): # 对θ进行积分
|
||||||
max_w = 0
|
max_w = 0
|
||||||
# 求暴露弧上一点的切线
|
# 求暴露弧上一点的切线
|
||||||
line_x = math.cos(theta) * rc + rc_x
|
line_x = math.cos(theta) * rc + dgc
|
||||||
line_y = math.sin(theta) * rc + rc_y
|
line_y = math.sin(theta) * rc + h_cav
|
||||||
k = math.tan(theta + math.pi / 2) # 入射角
|
k = math.tan(theta + math.pi / 2) # 入射角
|
||||||
# 求保护弧到直线的距离，判断是否相交
|
# 求保护弧到直线的距离，判断是否相交
|
||||||
d_to_rs = distance_point_line(rs_x, rs_y, line_x, line_y, k)
|
d_to_rs = distance_point_line(0, h_gav, 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, 0, h_gav, rs, init_k=k)
|
||||||
if new_k >= 0:
|
if new_k >= 0:
|
||||||
max_w = math.atan(new_k) # 用于保护弧相切的角度
|
max_w = math.atan(new_k) # 用于保护弧相切的角度
|
||||||
elif new_k < 0:
|
elif new_k < 0:
|
||||||
|
|
@ -239,23 +183,20 @@ def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y): # 对θ进行积分
|
||||||
return r
|
return r
|
||||||
|
|
||||||
|
|
||||||
def bd_area(
|
def bd_area(i_curt, u_ph, dgc, h_gav, h_cav): # 暴露弧的投影面积
|
||||||
i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, ground_surface, rg_type
|
theta1, theta2 = intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph) # θ角度
|
||||||
): # 暴露弧的投影面积
|
|
||||||
theta1, theta2 = intersection_angle(
|
|
||||||
rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_surface, rg_type
|
|
||||||
) # θ角度
|
|
||||||
theta_fineness = 0.01
|
theta_fineness = 0.01
|
||||||
rc = rc_fun(i_curt, u_ph)
|
rc = rc_fun(i_curt, u_ph)
|
||||||
rs = rs_fun(i_curt)
|
rs = rs_fun(i_curt)
|
||||||
rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
|
rg = rg_fun(i_curt, h_cav)
|
||||||
|
r_bd = 0
|
||||||
theta_sample, d_theta = np.linspace(
|
theta_sample, d_theta = np.linspace(
|
||||||
theta1, theta2, int((theta2 - theta1) / theta_fineness), retstep=True
|
theta1, theta2, int((theta2 - theta1) / theta_fineness), retstep=True
|
||||||
)
|
)
|
||||||
if len(theta_sample) < 2:
|
if len(theta_sample) < 2:
|
||||||
return 0
|
return 0
|
||||||
vec_calculus_bd = np.vectorize(calculus_bd)
|
vec_calculus_bd = np.vectorize(calculus_bd)
|
||||||
calculus_bd_np = vec_calculus_bd(theta_sample, rc, rs, rg, rc_x, rc_y, rs_x, rs_y)
|
calculus_bd_np = vec_calculus_bd(theta_sample, rc, rs, rg, dgc, h_cav, h_gav)
|
||||||
r_bd = np.sum(calculus_bd_np[:-1] + calculus_bd_np[1:]) / 2 * d_theta
|
r_bd = np.sum(calculus_bd_np[:-1] + calculus_bd_np[1:]) / 2 * d_theta
|
||||||
# for calculus_theta in theta_sample[:-1]:
|
# for calculus_theta in theta_sample[:-1]:
|
||||||
# r_bd += (
|
# r_bd += (
|
||||||
|
|
@ -273,7 +214,7 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
|
||||||
# 直线方程为 y-y0=k(x-x0)，x0和y0为经过直线的任意一点
|
# 直线方程为 y-y0=k(x-x0)，x0和y0为经过直线的任意一点
|
||||||
# 牛顿法求解k
|
# 牛顿法求解k
|
||||||
# f(k)=(k*x1-y1-k*x0+y0)**2-R**2*(k**2+1) x1,y1是圆心
|
# f(k)=(k*x1-y1-k*x0+y0)**2-R**2*(k**2+1) x1,y1是圆心
|
||||||
# 已考虑两个解的判别
|
|
||||||
k_candidate = [-100, 100]
|
k_candidate = [-100, 100]
|
||||||
if abs(center_y - line_y) < 1 and abs(line_x - center_x - radius) < 1:
|
if abs(center_y - line_y) < 1 and abs(line_x - center_x - radius) < 1:
|
||||||
# k不存在
|
# k不存在
|
||||||
|
|
@ -307,10 +248,10 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
|
||||||
# 把k转化成相应的角度，从x开始，逆时针为正
|
# 把k转化成相应的角度，从x开始，逆时针为正
|
||||||
k_angle = []
|
k_angle = []
|
||||||
for kk in k_candidate:
|
for kk in k_candidate:
|
||||||
# if kk is None:
|
if kk is None:
|
||||||
# abc = 123
|
abc = 123
|
||||||
# # tangent_line_k(line_x, line_y, center_x, center_y, radius)
|
# tangent_line_k(line_x, line_y, center_x, center_y, radius)
|
||||||
# pass
|
pass
|
||||||
if kk >= 0:
|
if kk >= 0:
|
||||||
k_angle.append(math.atan(kk))
|
k_angle.append(math.atan(kk))
|
||||||
if kk < 0:
|
if kk < 0:
|
||||||
|
|
@ -321,19 +262,3 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
|
||||||
|
|
||||||
def func_ng(td): # 地闪密度
|
def func_ng(td): # 地闪密度
|
||||||
return 0.023 * (td ** 1.3)
|
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
|
|
||||||
|
|
|
||||||
311
main.py
311
main.py
|
|
@ -1,212 +1,135 @@
|
||||||
|
import numpy as np
|
||||||
|
|
||||||
from core import *
|
from core import *
|
||||||
import timeit
|
import timeit
|
||||||
|
|
||||||
|
|
||||||
def egm():
|
def egm():
|
||||||
h_g_avr_sag = 11.67 * 2 / 3
|
avr_n_sf = 0 # 考虑电压的影响计算的跳闸率
|
||||||
h_c_avr_sag = 14.43 * 2 / 3
|
|
||||||
h_whole = 140 # 杆塔全高
|
|
||||||
voltage_n = 3 # 工作电压分成多少份来计算
|
voltage_n = 3 # 工作电压分成多少份来计算
|
||||||
td = 20 # 雷暴日
|
ng = func_ng(20)
|
||||||
|
h_whole = 140 # 杆塔全高
|
||||||
insulator_c_len = 6.8 # 串子绝缘长度
|
insulator_c_len = 6.8 # 串子绝缘长度
|
||||||
string_c_len = 9.2
|
string_c_len = 9.2
|
||||||
string_g_len = 0.5
|
string_g_len = 0.5
|
||||||
gc_x = [17.9, 16, 15, 16]
|
dgc = -0.0 # 导地线水平距离
|
||||||
|
vertical_dgc = 2.7 # 导地线挂点垂直距离
|
||||||
# 以后考虑地形角度,地面线
|
h_g_avr_sag = 11.67 * 2 / 3
|
||||||
def ground_surface(x):
|
h_c_avr_sag = 14.43 * 2 / 3
|
||||||
return 0
|
h_gav = h_whole - string_g_len - h_g_avr_sag # 地线对地平均高
|
||||||
|
h_cav = h_whole - string_c_len - vertical_dgc - h_c_avr_sag # 导线对地平均高
|
||||||
gc_y = [
|
shield_angle = math.atan(dgc / (vertical_dgc + string_c_len)) * 180 / math.pi
|
||||||
h_whole - string_g_len - h_g_avr_sag, # 地线对地平均高
|
print(f"保护角{shield_angle:.3f}°")
|
||||||
h_whole - string_c_len - h_c_avr_sag - 2.7, # 导线对地平均高
|
for u_bar in range(voltage_n):
|
||||||
h_whole - string_c_len - h_c_avr_sag - 20, # 导线对地平均高
|
u_ph = (
|
||||||
h_whole - string_c_len - h_c_avr_sag - 35.7, # 导线对地平均高
|
math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732
|
||||||
]
|
) # 运行相电压
|
||||||
if len(gc_y) > 2: # 双回路
|
# 迭代法计算最大电流
|
||||||
phase_n = 3 # 边相导线数量
|
i_max = 0
|
||||||
else:
|
i_min = min_i(insulator_c_len, u_ph / 1.732)
|
||||||
phase_n = 1
|
_min_i = i_min # 尝试的最小电流
|
||||||
#########################################################
|
_max_i = 200 # 尝试的最大电流
|
||||||
rg_type = None
|
cad = Draw()
|
||||||
# 以上是需要设置的参数
|
cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2)
|
||||||
avr_n_sf = 0 # 考虑电压的影响计算的跳闸率
|
for i_bar in np.linspace(_min_i, _max_i, int((_max_i - _min_i) / 0.1)): # 雷电流
|
||||||
rg_x = None
|
# print(f"尝试计算电流为{i_bar:.2f}")
|
||||||
rg_y = None
|
rs = rs_fun(i_bar)
|
||||||
cad = Draw()
|
rc = rc_fun(i_bar, u_ph)
|
||||||
for phase_conductor in range(phase_n):
|
rg = rg_fun(i_bar, h_cav)
|
||||||
rs_x = gc_x[phase_conductor]
|
#######
|
||||||
rs_y = gc_y[phase_conductor]
|
# cccCount += 1
|
||||||
rc_x = gc_x[phase_conductor + 1]
|
# if cccCount % 30 == 0:
|
||||||
rc_y = gc_y[phase_conductor + 1]
|
# import core
|
||||||
if phase_n == 1:
|
#
|
||||||
rg_type = "g"
|
# core.gMSP.add_circle((0, h_gav), rs)
|
||||||
if phase_n > 1: # 多回路
|
# core.gMSP.add_circle(
|
||||||
if phase_conductor < 2:
|
# (dgc, h_cav), rc_fun(i_bar, -u_ph), dxfattribs={"color": 4}
|
||||||
rg_type = "c"
|
# )
|
||||||
rg_x = gc_x[phase_conductor + 2]
|
# core.gMSP.add_circle((dgc, h_cav), rc)
|
||||||
rg_y = gc_y[phase_conductor + 2]
|
#######
|
||||||
else:
|
circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
|
||||||
rg_type = "g"
|
if not circle_intersection: # if circle_intersection is []
|
||||||
# TODO 保护角公式可能有问题，后面改
|
# print("保护弧和暴露弧无交点，检查设置参数。程序退出。")
|
||||||
shield_angle = (
|
continue
|
||||||
math.atan(rc_x / ((rc_y - rs_y) + string_c_len)) * 180 / math.pi
|
circle_rc_line_intersection = solve_circle_line_intersection(
|
||||||
) # 保护角
|
rc, rg, dgc, h_cav
|
||||||
print(f"保护角{shield_angle:.3f}°")
|
)
|
||||||
print(f"最低相防护标识{rg_type}")
|
if not circle_rc_line_intersection:
|
||||||
ng = func_ng(td)
|
continue
|
||||||
for u_bar in range(voltage_n):
|
min_distance_intersection = (
|
||||||
u_ph = (
|
np.sum(
|
||||||
math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732
|
(
|
||||||
) # 运行相电压
|
np.array(circle_intersection)
|
||||||
print(f"计算第{phase_conductor + 1}相，电压为{u_ph:.2f}kV")
|
- np.array(circle_rc_line_intersection)
|
||||||
# 迭代法计算最大电流
|
)
|
||||||
i_max = 0
|
** 2
|
||||||
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
|
|
||||||
)
|
)
|
||||||
i_max = i_bar
|
** 0.5
|
||||||
if not rg_rc_circle_intersection: # if circle_intersection is []
|
) # 计算两圆交点和地面直线交点的最小距离
|
||||||
print("保护弧和暴露弧无交点，检查设置参数。")
|
i_max = i_bar
|
||||||
continue
|
if min_distance_intersection < 0.1:
|
||||||
circle_rc_line_or_rg_intersection = None
|
break
|
||||||
if rg_type == "g":
|
if circle_intersection[1] < circle_rc_line_intersection[1]:
|
||||||
circle_rc_line_or_rg_intersection = solve_circle_line_intersection(
|
circle_rs_line_intersection = solve_circle_line_intersection(
|
||||||
rc, rg, rc_x, rc_y
|
rs, rg, 0, h_gav
|
||||||
)
|
)
|
||||||
elif rg_type == "c":
|
# 判断与保护弧的交点是否在暴露弧外面
|
||||||
circle_rc_line_or_rg_intersection = solve_circle_intersection(
|
distance = (
|
||||||
rg, rc, rg_x, rg_y, rc_x, rc_y
|
|
||||||
)
|
|
||||||
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.sum(
|
||||||
(
|
(np.array(circle_rs_line_intersection) - np.array([dgc, h_cav]))
|
||||||
np.array(rg_rc_circle_intersection)
|
|
||||||
- np.array(circle_rc_line_or_rg_intersection)
|
|
||||||
)
|
|
||||||
** 2
|
** 2
|
||||||
)
|
)
|
||||||
** 0.5
|
** 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.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)
|
if distance > rc:
|
||||||
)
|
print("暴露弧已经完全被屏蔽")
|
||||||
calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt
|
break
|
||||||
# for i_curt in i_curt_samples[:-1]:
|
cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2)
|
||||||
# cal_bd_first = bd_area(i_curt, u_ph, dgc, h_gav, h_cav)
|
cad.draw(i_max, u_ph, h_gav, h_cav, dgc, 6)
|
||||||
# cal_bd_second = bd_area(i_curt + d_curt, u_ph, dgc, h_gav, h_cav)
|
cad.save()
|
||||||
# cal_thunder_density_first = thunder_density(i_curt)
|
# 判断是否导线已经被完全保护
|
||||||
# cal_thunder_density_second = thunder_density(i_curt + d_curt)
|
if abs(i_max - _max_i) < 1e-5:
|
||||||
# calculus += (
|
print("无法找到最大电流，可能是杆塔较高。")
|
||||||
# (
|
print(f"最大电流设置为自然界最大电流{i_max}kA")
|
||||||
# cal_bd_first * cal_thunder_density_first
|
print(f"最大电流为{i_max:.2f}")
|
||||||
# + cal_bd_second * cal_thunder_density_second
|
print(f"最小电流为{i_min:.2f}")
|
||||||
# )
|
curt_fineness = 0.1 # 电流积分细度
|
||||||
# / 2
|
if i_min > i_max or abs(i_min - i_max) < curt_fineness:
|
||||||
# * d_curt
|
print("最大电流小于最小电流，没有暴露弧，程序结束。")
|
||||||
# )
|
return
|
||||||
# if abs(calculus-0.05812740052770032)<1e-5:
|
# 开始积分
|
||||||
# abc=123
|
curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份
|
||||||
# pass
|
calculus = 0
|
||||||
n_sf = (
|
i_curt_samples, d_curt = np.linspace(
|
||||||
2 * ng / 10 * calculus
|
i_min, i_max, curt_segment_n + 1, retstep=True
|
||||||
) # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
|
)
|
||||||
avr_n_sf += n_sf / voltage_n
|
bd_area_vec = np.vectorize(bd_area)
|
||||||
print(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}")
|
cal_bd_np = bd_area_vec(
|
||||||
print(f"跳闸率是{avr_n_sf:.6f}")
|
i_curt_samples, u_ph, dgc, h_gav, h_cav
|
||||||
|
) * 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}")
|
||||||
|
|
||||||
|
|
||||||
def speed():
|
def speed():
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue