370 lines
13 KiB
Python
370 lines
13 KiB
Python
import math
|
||
import ezdxf
|
||
import numpy as np
|
||
|
||
gCAD = None
|
||
gMSP = None
|
||
gCount = 1
|
||
|
||
|
||
class Draw:
|
||
def __init__(self):
|
||
self._doc = ezdxf.new(dxfversion="R2010")
|
||
self._doc.layers.add("EGM", color=2)
|
||
global gCAD
|
||
gCAD = self
|
||
|
||
def draw(self, i_curt, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, color):
|
||
doc = self._doc
|
||
msp = doc.modelspace()
|
||
global gMSP
|
||
gMSP = msp
|
||
rs = rs_fun(i_curt)
|
||
rc = rc_fun(i_curt, u_ph)
|
||
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_line((rc_x, 0), (rc_x, rc_y)) # 导线
|
||
msp.add_line((rs_x, rs_y), (rc_x, rc_y))
|
||
# 角度线
|
||
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
|
||
)
|
||
if rg_rc_intersection:
|
||
msp.add_line(
|
||
(rc_x, rc_y), rg_rc_intersection, dxfattribs={"color": color}
|
||
) # 圆和圆的交点
|
||
# 计算圆交点
|
||
|
||
# msp.add_line((dgc, h_cav), circle_intersection) # 导线
|
||
|
||
def save(self):
|
||
doc = self._doc
|
||
doc.saveas("egm.dxf")
|
||
|
||
def saveas(self, file_name):
|
||
doc = self._doc
|
||
doc.saveas(file_name)
|
||
|
||
|
||
# 圆交点
|
||
def solve_circle_intersection(
|
||
radius1,
|
||
radius2,
|
||
center_x1,
|
||
center_y1,
|
||
center_x2,
|
||
center_y2,
|
||
):
|
||
# 用牛顿法求解
|
||
x = radius2 + center_x2 # 初始值
|
||
y = radius2 + center_y2 # 初始值
|
||
# TODO 考虑出现2个解的情况
|
||
for bar in range(0, 10):
|
||
A = [
|
||
[-2 * (x - center_x1), -2 * (y - center_y1)],
|
||
[-2 * (x - center_x2), -2 * (y - center_y2)],
|
||
]
|
||
b = [
|
||
(x - center_x1) ** 2 + (y - center_y1) ** 2 - radius1 ** 2,
|
||
(x - center_x2) ** 2 + (y - center_y2) ** 2 - radius2 ** 2,
|
||
]
|
||
X_set = np.linalg.solve(A, b)
|
||
x += X_set[0]
|
||
y += X_set[1]
|
||
if np.all(np.abs(X_set) < 1e-5):
|
||
return [x, y]
|
||
return []
|
||
|
||
|
||
# 圆与捕雷线交点
|
||
def solve_circle_line_intersection(radius, rg, center_x, center_y):
|
||
distance = distance_point_line(center_x, center_y, 0, rg, 0) # 捕雷线到暴露圆中点的距离
|
||
if distance > radius:
|
||
return []
|
||
else:
|
||
r = (radius ** 2 - (rg - center_y) ** 2) ** 0.5 + center_x
|
||
return [r, rg]
|
||
|
||
|
||
def min_i(string_len, u_ph):
|
||
u_50 = 530 * string_len + 35
|
||
z_0 = 300 # 雷电波阻抗
|
||
z_c = 251 # 导线波阻抗
|
||
r = (u_50 + 2 * z_0 / (2 * z_0 + z_c) * u_ph) * (2 * z_0 + z_c) / (z_0 * z_c)
|
||
return r
|
||
|
||
|
||
def thunder_density(i): # l雷电流幅值密度函数
|
||
r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44)
|
||
return r
|
||
|
||
|
||
def angel_density(angle): # 入射角密度函数 angle单位是弧度
|
||
r = 0.75 * abs((np.cos(angle - math.pi / 2) ** 3))
|
||
return r
|
||
|
||
|
||
def rs_fun(i):
|
||
r = 10 * (i ** 0.65)
|
||
return r
|
||
|
||
|
||
def rc_fun(i, u_ph):
|
||
r = 1.63 * ((5.015 * (i ** 0.578) - 0.001 * u_ph) ** 1.125)
|
||
return r
|
||
|
||
|
||
# typ 如果是g，代表捕雷线公式，c代表暴露弧公式
|
||
def rg_fun(i_curt, h_cav, u_ph, typ="g"):
|
||
rg = None
|
||
if typ == "g":
|
||
if h_cav < 40:
|
||
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
|
||
|
||
|
||
def intersection_angle(
|
||
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)
|
||
rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
|
||
circle_intersection = solve_circle_intersection(
|
||
rs, rc, rs_x, rs_y, rc_x, rc_y
|
||
) # 两圆的交点
|
||
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
|
||
) # 两圆的交点
|
||
if circle_line_or_rg_intersection:
|
||
(
|
||
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
|
||
)
|
||
theta1 = None
|
||
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])
|
||
np_circle_line_or_rg_intersection = np.array(circle_line_or_rg_intersection)
|
||
if not circle_line_or_rg_intersection:
|
||
if rc_y - rc > rg: # rg在rc下面
|
||
# 捕捉线太低了，对高塔无保护，θ_1从-90°开始计算。
|
||
theta1 = -math.pi / 2
|
||
else:
|
||
theta1_line = np_circle_line_or_rg_intersection - np.array([rc_x, rc_y])
|
||
theta1 = math.atan(theta1_line[1] / theta1_line[0])
|
||
return np.array([theta1, theta2])
|
||
|
||
|
||
# 点到直线的距离
|
||
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)
|
||
return d
|
||
|
||
|
||
def func_calculus_pw(theta, max_w):
|
||
w_fineness = 0.01
|
||
segments = int(max_w / w_fineness)
|
||
if segments < 2: # 最大最小太小，没有可以积分的
|
||
return 0
|
||
w_samples, d_w = np.linspace(0, max_w, segments, retstep=True)
|
||
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
|
||
return r_pw
|
||
|
||
|
||
def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y): # 对θ进行积分
|
||
max_w = 0
|
||
# 求暴露弧上一点的切线
|
||
line_x = math.cos(theta) * rc + rc_x
|
||
line_y = math.sin(theta) * rc + rc_y
|
||
k = math.tan(theta + math.pi / 2) # 入射角
|
||
# 求保护弧到直线的距离，判断是否相交
|
||
d_to_rs = distance_point_line(rs_x, rs_y, line_x, line_y, k)
|
||
if d_to_rs < rs: # 相交
|
||
# 要用过直线上一点到暴露弧的切线
|
||
new_k = tangent_line_k(line_x, line_y, rs_x, rs_y, rs, init_k=k)
|
||
if new_k >= 0:
|
||
max_w = math.atan(new_k) # 用于保护弧相切的角度
|
||
elif new_k < 0:
|
||
max_w = math.atan(new_k) + math.pi
|
||
# TODO to be removed
|
||
# global gCount
|
||
# gCount = gCount+1
|
||
# if gCount % 100 == 0:
|
||
# gMSP.add_circle((0, h_gav), rs)
|
||
# gMSP.add_circle((dgc, h_cav), rc)
|
||
# gMSP.add_line((dgc, h_cav), (line_x, line_y))
|
||
# gMSP.add_line(
|
||
# (-500, new_k * (-500 - line_x) + line_y),
|
||
# (500, new_k * (500 - line_x) + line_y),
|
||
# )
|
||
# gCAD.save()
|
||
# pass
|
||
else:
|
||
max_w = theta + math.pi / 2 # 入射角
|
||
# TODO to be removed
|
||
if gCount % 200 == 0:
|
||
# # intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph)
|
||
# gMSP.add_circle((0, h_gav), rs)
|
||
# gMSP.add_circle((dgc, h_cav), rc)
|
||
# gMSP.add_line((dgc, h_cav), (line_x, line_y))
|
||
# gMSP.add_line(
|
||
# (-500, k * (-500 - line_x) + line_y),
|
||
# (500, k * (500 - line_x) + line_y),
|
||
# )
|
||
# gCAD.save()
|
||
pass
|
||
r = rc / math.cos(theta) * func_calculus_pw(theta, max_w)
|
||
return r
|
||
|
||
|
||
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, ground_surface, rg_type
|
||
) # θ角度
|
||
theta_fineness = 0.01
|
||
rc = rc_fun(i_curt, u_ph)
|
||
rs = rs_fun(i_curt)
|
||
rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
|
||
theta_segments = int((theta2 - theta1) / theta_fineness)
|
||
if theta_segments < 2:
|
||
return 0
|
||
theta_sample, d_theta = np.linspace(theta1, theta2, theta_segments, retstep=True)
|
||
if len(theta_sample) < 2:
|
||
return 0
|
||
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)
|
||
r_bd = np.sum(calculus_bd_np[:-1] + calculus_bd_np[1:]) / 2 * d_theta
|
||
# for calculus_theta in theta_sample[:-1]:
|
||
# r_bd += (
|
||
# (
|
||
# calculus_bd(calculus_theta, rc, rs, rg, dgc, h_cav, h_gav)
|
||
# + calculus_bd(calculus_theta + d_theta, rc, rs, rg, dgc, h_cav, h_gav)
|
||
# )
|
||
# / 2
|
||
# * d_theta
|
||
# )
|
||
return r_bd
|
||
|
||
|
||
def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
|
||
# 直线方程为 y-y0=k(x-x0)，x0和y0为经过直线的任意一点
|
||
# 牛顿法求解k
|
||
# f(k)=(k*x1-y1-k*x0+y0)**2-R**2*(k**2+1) x1,y1是圆心
|
||
# 已考虑两个解的判别
|
||
k_candidate = [-100, 100]
|
||
if abs(center_y - line_y) < 1 and abs(line_x - center_x - radius) < 1:
|
||
# k不存在
|
||
k_candidate = [99999999, 99999999]
|
||
else:
|
||
for ind, k_cdi in enumerate(list(k_candidate)):
|
||
k = k_candidate[ind]
|
||
k_candidate[ind] = None
|
||
max_iteration = 30
|
||
for bar in range(0, max_iteration):
|
||
fk = (k * center_x - center_y - k * line_x + line_y) ** 2 - (
|
||
radius ** 2
|
||
) * (k ** 2 + 1)
|
||
|
||
d_fk = (
|
||
2
|
||
* (k * center_x - center_y - k * line_x + line_y)
|
||
* (center_x - line_x)
|
||
- 2 * (radius ** 2) * k
|
||
)
|
||
if abs(d_fk) < 1e-5 and abs(line_x - center_x - radius) < 1e-5:
|
||
# k不存在，角度为90°，k取一个很大的正数
|
||
k_candidate[ind] = 99999999999999
|
||
break
|
||
d_k = -fk / d_fk
|
||
k += d_k
|
||
if abs(d_k) < 1e-3:
|
||
dd = distance_point_line(center_x, center_y, line_x, line_y, k)
|
||
if abs(dd - radius) < 1:
|
||
k_candidate[ind] = k
|
||
break
|
||
# 解决数值稳定性
|
||
if bar == max_iteration - 1:
|
||
if abs(math.atan(k)) * 180 / math.pi > 89:
|
||
k_candidate[ind] = k
|
||
# 把k转化成相应的角度，从x开始，逆时针为正
|
||
k_angle = []
|
||
for kk in k_candidate:
|
||
# if kk is None:
|
||
# abc = 123
|
||
# tangent_line_k(line_x, line_y, center_x, center_y, radius)
|
||
# pass
|
||
if kk >= 0:
|
||
k_angle.append(math.atan(kk))
|
||
if kk < 0:
|
||
k_angle.append(math.pi + math.atan(kk))
|
||
# 返回相对x轴最大的角度k
|
||
return np.array(k_candidate)[np.max(k_angle) == k_angle].tolist()[-1]
|
||
|
||
|
||
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
|
||
|
||
|
||
# u_ph是相电压
|
||
# insulator_c_len绝缘子闪络距离
|
||
def arc_possibility(rated_voltage, insulator_c_len): # 建弧率
|
||
_e = rated_voltage / (3 ** 0.5) / insulator_c_len
|
||
r = (4.5 * (_e ** 0.75) - 14) * 1e-2
|
||
return r
|