考虑了切线的k为负的情况
This commit is contained in:
parent
db2788f116
commit
b98d2534ab
289
main.py
289
main.py
|
|
@ -26,7 +26,7 @@ class Draw:
|
|||
msp.add_circle((dgc, h_cav), rc, dxfattribs={"color": color})
|
||||
msp.add_line((dgc, 0), (dgc, h_cav)) # 导线
|
||||
msp.add_line((0, h_gav), (dgc, h_cav))
|
||||
msp.add_line((0, rg), (200, rg), dxfattribs={"color": color})
|
||||
msp.add_line((0, rg), (2000, rg), dxfattribs={"color": color})
|
||||
# 计算圆交点
|
||||
# circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
|
||||
# msp.add_line((0, h_gav), circle_intersection) # 地线
|
||||
|
|
@ -42,8 +42,8 @@ class Draw:
|
|||
# 圆交点
|
||||
def solve_circle_intersection(rs, rc, hgav, hcav, dgc):
|
||||
# 用牛顿法求解
|
||||
x = 300
|
||||
y = 300
|
||||
x = rc # 初始值
|
||||
y = rc # 初始值
|
||||
for bar in range(0, 10):
|
||||
A = [[-2 * x, -2 * (y - hgav)], [-2 * (x - dgc), -2 * (y - hcav)]]
|
||||
b = [
|
||||
|
|
@ -79,7 +79,7 @@ def thunder_density(i): # l雷电流幅值密度函数
|
|||
|
||||
|
||||
def angel_density(angle): # 入射角密度函数 angle单位是弧度
|
||||
r = 0.75 * (math.cos(angle) ** 3)
|
||||
r = 0.75 * (math.cos(angle - math.pi / 2) ** 3)
|
||||
return r
|
||||
|
||||
|
||||
|
|
@ -110,6 +110,8 @@ def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph): # 暴露弧的角度
|
|||
rc, rg, h_cav, dgc
|
||||
) # 暴露圆和补雷线的交点
|
||||
np_circle_intersection = np.array(circle_intersection)
|
||||
if not circle_intersection:
|
||||
abc = 123
|
||||
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)
|
||||
|
|
@ -127,39 +129,83 @@ def distance_point_line(point_x, point_y, line_x, line_y, k):
|
|||
return d
|
||||
|
||||
|
||||
def fun_calculus_pw(theta, max_w):
|
||||
w_fineness = 0.01
|
||||
r_pw = 0
|
||||
if int(max_w / w_fineness) < 0:
|
||||
abc = 123
|
||||
pass
|
||||
w_samples, d_w = np.linspace(0, max_w, int(max_w / w_fineness) + 1, retstep=True)
|
||||
for cal_w in w_samples:
|
||||
r_pw += (
|
||||
(
|
||||
abs(angel_density(cal_w)) * math.sin(theta - cal_w + math.pi)
|
||||
+ abs(angel_density(cal_w + d_w))
|
||||
* math.sin(theta - cal_w + math.pi - d_w)
|
||||
)
|
||||
/ 2
|
||||
) * d_w
|
||||
return r_pw
|
||||
|
||||
|
||||
def cal_bd(theta, rc, rs, rg, dgc, h_cav, h_gav):
|
||||
# 求暴露弧上一点的切线
|
||||
line_x = math.cos(theta) * rc + dgc
|
||||
line_y = math.sin(theta) * rc + h_cav
|
||||
max_w = theta + math.pi / 2 # 入射角
|
||||
k = math.tan(max_w)
|
||||
# 求保护弧到直线的距离，判断是否相交
|
||||
d_to_rs = distance_point_line(0, h_gav, line_x, line_y, k)
|
||||
if d_to_rs < rs: # 相交
|
||||
# 要用过直线上一点到暴露弧的切线
|
||||
new_k = tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k)
|
||||
if not new_k:
|
||||
a = 12
|
||||
tangent_line_k(line_x, line_y, 0, h_gav, 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
|
||||
if max_w < 0:
|
||||
abc = 123
|
||||
tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k)
|
||||
# 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)
|
||||
# )
|
||||
# gMSP.add_line((0, rg), (1000, rg))
|
||||
# gCAD.save()
|
||||
tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k)
|
||||
|
||||
pass
|
||||
r = rc / math.cos(theta) * fun_calculus_pw(theta, max_w)
|
||||
return r
|
||||
|
||||
|
||||
def bd_area(i_curt, u_ph, dgc, h_gav, h_cav): # 暴露弧的投影面积
|
||||
theta1, theta2 = intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph)
|
||||
theta_fineness = 0.01
|
||||
rc = rc_fun(i_curt, u_ph)
|
||||
rs = rs_fun(i_curt)
|
||||
rg = rg_fun(i_curt, h_cav)
|
||||
# 求暴露弧上一点的切线
|
||||
line_x = math.cos(theta1) * rc + dgc
|
||||
line_y = math.sin(theta1) * rc + h_cav
|
||||
max_w = 0 # 入射角
|
||||
if theta1 < 0:
|
||||
max_w = theta1 + math.pi / 2
|
||||
k = math.tan(max_w)
|
||||
# 求保护弧到直线的距离，判断是否相交
|
||||
d_to_rs = distance_point_line(0, h_gav, line_x, line_y, k)
|
||||
if d_to_rs < rs: # 相交
|
||||
# 要用过直线上一点到暴露弧的切线
|
||||
new_k = tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k)
|
||||
max_w = math.atan(new_k) # 用于保护弧相切的角度
|
||||
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)
|
||||
r_bd = 0
|
||||
theta_sample, d_theta = np.linspace(
|
||||
theta1, theta2, int((theta2 - theta1) / theta_fineness) + 1, retstep=True
|
||||
)
|
||||
for cal_theta in theta_sample[:-1]:
|
||||
r_bd += (
|
||||
(
|
||||
cal_bd(cal_theta, rc, rs, rg, dgc, h_cav, h_gav)
|
||||
+ cal_bd(cal_theta + d_theta, rc, rs, rg, dgc, h_cav, h_gav)
|
||||
)
|
||||
gMSP.add_line((0, rg), (1000, rg))
|
||||
gCAD.save()
|
||||
pass
|
||||
/ 2
|
||||
* d_theta
|
||||
)
|
||||
return r_bd
|
||||
|
||||
# k=tangent_line_k(point_x, point_y, dgc, h_cav,rc)
|
||||
# 暂时不考虑雷电入射角的影响
|
||||
r = (math.cos(theta1) - math.cos(theta2)) * rc
|
||||
return r
|
||||
# r1=rc*(-math.cos(thyta2)+math.cos(thyta1))
|
||||
# 入射角密度函数积分
|
||||
# arrival_angle_fineness=0.0001
|
||||
|
|
@ -170,94 +216,113 @@ 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是圆心
|
||||
# TODO:应该找到两个角度值后再比较
|
||||
k = init_k
|
||||
for bar in range(0, 30):
|
||||
fk = (k * center_x - center_y - k * line_x + line_y) ** 2 - (radius ** 2) * (
|
||||
k ** 2 + 1
|
||||
)
|
||||
# TODO:需要检验k值不存在的情况
|
||||
k_candidate = [-100, 100]
|
||||
for ind, k_cdi in enumerate(list(k_candidate)):
|
||||
k = k_candidate[ind]
|
||||
k_candidate[ind] = None
|
||||
for bar in range(0, 30):
|
||||
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
|
||||
)
|
||||
d_k = -fk / d_fk
|
||||
k += d_k
|
||||
if abs(d_k) < 1e-5:
|
||||
dd = distance_point_line(center_x, center_y, line_x, line_y, k)
|
||||
if abs(dd - radius) < 1e-5:
|
||||
return k
|
||||
return None
|
||||
d_fk = (
|
||||
2
|
||||
* (k * center_x - center_y - k * line_x + line_y)
|
||||
* (center_x - line_x)
|
||||
- 2 * (radius ** 2) * k
|
||||
)
|
||||
d_k = -fk / d_fk
|
||||
k += d_k
|
||||
if abs(d_k) < 1e-5:
|
||||
dd = distance_point_line(center_x, center_y, line_x, line_y, k)
|
||||
if abs(dd - radius) < 1e-5:
|
||||
k_candidate[ind] = k
|
||||
break
|
||||
# 把k转化成相应的角度，从x开始，逆时针为正
|
||||
k_angle = []
|
||||
for kk in k_candidate:
|
||||
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 egm():
|
||||
u_ph = 750 / 1.732 # 运行相电压
|
||||
h_cav = 160 # 导线对地平均高
|
||||
h_gav = h_cav + 9.5 + 2.7
|
||||
dgc = -2 # 导地线水平距离
|
||||
# 迭代法计算最大电流
|
||||
i_max = 0
|
||||
_min_i = 20 # 尝试的最小电流
|
||||
_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
|
||||
for u_bar in range(1):
|
||||
u_ph = math.sqrt(2) * 750 * math.cos(2 * math.pi / 6 * 0) / 1.732 # 运行相电压
|
||||
h_cav = 90 # 导线对地平均高
|
||||
h_gav = h_cav + 9.5 + 2.7
|
||||
dgc = 2.9 # 导地线水平距离
|
||||
# 迭代法计算最大电流
|
||||
i_max = 0
|
||||
_min_i = 20 # 尝试的最小电流
|
||||
_max_i = 300 # 尝试的最大电流
|
||||
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)
|
||||
# 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
|
||||
)
|
||||
** 0.5
|
||||
) # 计算两圆交点和地面直线交点的最小距离
|
||||
i_max = i_bar
|
||||
if min_distance_intersection < 0.1:
|
||||
break
|
||||
i_min = min_i(6.78, 750 / 1.732)
|
||||
cad = Draw()
|
||||
cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2)
|
||||
cad.draw(i_max, u_ph, h_gav, h_cav, dgc, 6)
|
||||
cad.save()
|
||||
if abs(i_max - _max_i) < 1e-5:
|
||||
print("无法找到最大电流，可能是杆塔较高。")
|
||||
i_max = 300
|
||||
print(f"最大电流设置为自然界最大电流{i_max}kA")
|
||||
print(f"最大电流为{i_max:.2f}")
|
||||
print(f"最小电流为{i_min:.2f}")
|
||||
if i_min > i_max:
|
||||
print("最大电流小于最小电流，没有暴露弧，程序结束。")
|
||||
return
|
||||
# 开始积分
|
||||
curt_fineness = 0.1 # 电流积分细度
|
||||
curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份
|
||||
calculus = 0
|
||||
i_curt_samples, d_curt = np.linspace(i_min, i_max, curt_segment_n + 1, retstep=True)
|
||||
for i_curt in i_curt_samples:
|
||||
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
|
||||
min_distance_intersection = (
|
||||
np.sum(
|
||||
(np.array(circle_intersection) - np.array(circle_line_intersection))
|
||||
** 2
|
||||
)
|
||||
** 0.5
|
||||
) # 计算两圆交点和地面直线交点的最小距离
|
||||
i_max = i_bar
|
||||
if min_distance_intersection < 0.1:
|
||||
break
|
||||
i_min = min_i(6.78, u_ph / 1.732)
|
||||
cad = Draw()
|
||||
cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2)
|
||||
cad.draw(i_max, u_ph, h_gav, h_cav, dgc, 6)
|
||||
cad.save()
|
||||
if abs(i_max - _max_i) < 1e-5:
|
||||
print("无法找到最大电流，可能是杆塔较高。")
|
||||
i_max = 300
|
||||
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("最大电流小于最小电流，没有暴露弧，程序结束。")
|
||||
return
|
||||
# 开始积分
|
||||
curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份
|
||||
calculus = 0
|
||||
i_curt_samples, d_curt = np.linspace(
|
||||
i_min, i_max, curt_segment_n + 1, retstep=True
|
||||
)
|
||||
n_sf = 2 * 2.7 / 10 * calculus # 调整率
|
||||
print(f"跳闸率是{n_sf:.6}")
|
||||
for i_curt in i_curt_samples:
|
||||
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
|
||||
)
|
||||
n_sf = 2 * 2.7 / 10 * calculus # 调整率
|
||||
print(f"跳闸率是{n_sf:.6}")
|
||||
|
||||
# draw(rs, rc, rg, h_gav, h_cav, dgc)
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue