200 lines
6.4 KiB
Python
200 lines
6.4 KiB
Python
import math
|
|
import ezdxf
|
|
import numpy as np
|
|
|
|
|
|
# 圆交点
|
|
def solve_circle_intersection(rs, rc, hgav, hcav, dgc):
|
|
# x = Symbol('x', real=True)
|
|
# y = Symbol('y', real=True)
|
|
# equ = [
|
|
# x ** 2 + (y - hgav) ** 2 - rs ** 2,
|
|
# (x - dgc) ** 2 + (y - hcav) ** 2 - rc ** 2,
|
|
# ]
|
|
# 用牛顿法求解
|
|
x = 300
|
|
y = 300
|
|
for bar in range(0, 10):
|
|
A = [[-2 * x, -2 * (y - hgav)], [-2 * (x - dgc), -2 * (y - hcav)]]
|
|
b = [
|
|
x ** 2 + (y - hgav) ** 2 - rs ** 2,
|
|
(x - dgc) ** 2 + (y - hcav) ** 2 - rc ** 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 []
|
|
# list_set = list(X_set)
|
|
# solve_set = nonlinsolve(equ, [x, y])
|
|
# print(ask(Q.real(solve_set)))
|
|
# list_set = list(solve_set)
|
|
# pprint(list_set)
|
|
# if not np.all(np.isreal(list_set)):
|
|
# return []
|
|
# for value in list_set:
|
|
# if value[0] > 0 and value[1] > 1:
|
|
# return value
|
|
# return []
|
|
|
|
|
|
# 圆与地面线交点
|
|
def solve_circle_line_intersection(rc, rg, hcav, dgc):
|
|
r = (rc ** 2 - (rg - hcav) ** 2) ** 0.5 + dgc
|
|
return [r, rg]
|
|
|
|
|
|
def draw(rs, rc, rg, h_gav, h_cav, dgc):
|
|
doc = ezdxf.new(dxfversion="R2010")
|
|
doc.layers.add("EGM", color=2)
|
|
msp = doc.modelspace()
|
|
msp.add_circle((0, h_gav), rs)
|
|
msp.add_line((0, 0), (0, h_gav)) # 地线
|
|
msp.add_circle((dgc, h_cav), rc)
|
|
msp.add_line((dgc, 0), (dgc, h_cav)) # 导线
|
|
msp.add_line((0, h_gav), (dgc, h_cav))
|
|
msp.add_line((0, rg), (200, rg))
|
|
# 计算圆交点
|
|
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) # 导线
|
|
circle_line_section = solve_circle_line_intersection(rc, rg, h_cav, dgc)
|
|
msp.add_line((0, 0), circle_line_section) # 导线和圆的交点
|
|
doc.saveas("egm.dxf")
|
|
solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
|
|
|
|
|
|
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 * (math.cos(angle) ** 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
|
|
|
|
|
|
def rg_fun(i, h_cav):
|
|
if h_cav < 40:
|
|
rg = (3.6 + 1.7 ** math.log(43 - h_cav)) ** 0.65
|
|
else:
|
|
rg = 5.5 * (i ** 0.65)
|
|
return rg
|
|
|
|
|
|
def intersection_angel(dgc, h_gav, h_cav, i_curt, u_ph): # 暴露弧的角度
|
|
rs = rs_fun(i_curt)
|
|
rc = rc_fun(i_curt, u_ph)
|
|
rg = rg_fun(i_curt, h_cav)
|
|
circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
|
|
circle_line_intersection = solve_circle_line_intersection(rc, rg, h_cav, dgc)
|
|
np_circle_intersection = np.array(circle_intersection)
|
|
theta2_line = np_circle_intersection - np.array([dgc, h_cav])
|
|
theta2 = math.atan(theta2_line[1] / theta2_line[0])
|
|
np_circle_line_intersection = np.array(circle_line_intersection)
|
|
theta1_line = np_circle_line_intersection - np.array([dgc, h_cav])
|
|
theta1 = math.atan(theta1_line[1] / theta1_line[0])
|
|
return np.array([theta1, theta2])
|
|
|
|
|
|
def bd_area(i_curt, u_ph, theta1, theta2): # 暴露弧的投影面积
|
|
rc = rc_fun(i_curt, u_ph)
|
|
# 暂时不考虑雷电入射角的影响
|
|
r = (math.cos(theta1) - math.cos(theta2)) * rc
|
|
return r
|
|
# r1=rc*(-math.cos(thyta2)+math.cos(thyta1))
|
|
# 入射角密度函数积分
|
|
# arrival_angle_fineness=0.0001
|
|
# for calculus_arv_angle in np.linspace()
|
|
|
|
|
|
def egm():
|
|
u_ph = 750 / 1.732 # 运行相电压
|
|
h_cav = 60 # 导线对地平均高
|
|
h_gav = h_cav + 9.5 + 7.2
|
|
dgc = 2
|
|
# 迭代法计算最大电流
|
|
i_max = 0
|
|
_min_i = 30 # 尝试的最小电流
|
|
_max_i = 80 # 尝试的最大电流
|
|
for i_bar in np.linspace(_min_i, _max_i, int((_max_i - _min_i) / 0.01)): # 雷电流
|
|
print(f"尝试计算电流为{i_bar:.2f}")
|
|
rs = rs_fun(i_bar)
|
|
if not np.isreal(rs):
|
|
continue
|
|
rc = rc_fun(i_bar, u_ph)
|
|
if not np.isreal(rc):
|
|
continue
|
|
rg = rg_fun(i_bar, h_cav)
|
|
if not np.isreal(rg):
|
|
continue
|
|
circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
|
|
if not circle_intersection: # if circle_intersection is []
|
|
continue
|
|
circle_line_intersection = solve_circle_line_intersection(rc, rg, h_cav, dgc)
|
|
min_distance_intersection = (
|
|
np.sum(
|
|
(np.array(circle_intersection) - np.array(circle_line_intersection))
|
|
** 2
|
|
)
|
|
** 0.5
|
|
) # 计算两圆交点和地面直线交点的最小距离
|
|
if min_distance_intersection < 0.01:
|
|
i_max = i_bar
|
|
draw(rs, rc, rg, h_gav, h_cav, dgc)
|
|
break
|
|
print(f"最大电流为{i_max:.2f}")
|
|
i_min = min_i(6.78, 750 / 1.732)
|
|
print(f"最小电流为{i_min:.2f}")
|
|
# 开始积分
|
|
curt_fineness = 0.1 # 电流积分细度
|
|
curt_segment_n = int((i_max - i_min) / curt_fineness)
|
|
d_curt = (i_max - i_min) / curt_segment_n
|
|
calculus = 0
|
|
for curt in np.linspace(i_min, i_max, curt_segment_n):
|
|
cal_thyta_first = intersection_angel(dgc, h_gav, h_cav, curt, u_ph)
|
|
cal_bd_first = bd_area(curt, u_ph, cal_thyta_first[0], cal_thyta_first[1])
|
|
cal_thyta_second = intersection_angel(dgc, h_gav, h_cav, curt + d_curt, u_ph)
|
|
cal_bd_second = bd_area(
|
|
curt + d_curt, u_ph, cal_thyta_second[0], cal_thyta_second[1]
|
|
)
|
|
cal_thunder_density_first = thunder_density(curt)
|
|
cal_thunder_density_second = thunder_density(curt + d_curt)
|
|
calculus += (
|
|
(
|
|
cal_bd_first * cal_thunder_density_first
|
|
+ cal_bd_second * cal_thunder_density_second
|
|
)
|
|
/ 2
|
|
* d_curt
|
|
)
|
|
n_sf=2*2.7/10*calculus
|
|
print(f'跳闸率是{n_sf}')
|
|
|
|
# draw(rs, rc, rg, h_gav, h_cav, dgc)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
thunder_density(2)
|
|
egm()
|
|
print("Finished.")
|