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

核心代码独立出来

This commit is contained in:
facat 2021-09-20 20:51:09 +08:00
parent 1cc8070c34
commit 2b898e41db
2 changed files with 299 additions and 280 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

276
core.py Normal file
View File

@ -0,0 +1,276 @@
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, h_gav, h_cav, dgc, 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, h_cav)
msp.add_circle((0, h_gav), rs, dxfattribs={"color": color})
msp.add_line((0, 0), (0, h_gav)) # 地线
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), (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) # 地线
# 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):
doc = self._doc
doc.saveas("egm.dxf")
# 圆交点
def solve_circle_intersection(rs, rc, hgav, hcav, dgc):
# 用牛顿法求解
x = rc # 初始值
y = rc # 初始值
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 []
# 圆与地面线交点
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 * (math.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)
# r=14.7*(i**0.42)
return r
def rg_fun(i_curt, h_cav):
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)
return rg
def intersection_angle(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, dgc, h_cav
) # 暴露圆和补雷线的交点
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)
theta1_line = np_circle_line_intersection - np.array([dgc, h_cav])
theta1 = math.atan(theta1_line[1] / theta1_line[0])
# 考虑雷电入射角度，所以theta1可以小于0，即计算从侧面击中的雷
# if theta1 < 0:
# # print(f"θ_1角度为负数{theta1:.4f},人为设置为0")
# theta1 = 0
return np.array([theta1, theta2])
def distance_point_line(point_x, point_y, line_x, line_y, k):
d = abs(k * point_x - point_y - k * line_x + line_y) / ((k ** 2 + 1) ** 0.5)
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[: -1]:
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 - d_w + math.pi)
)
/ 2
) * d_w
return r_pw
def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): # 对θ进行积分
max_w=0
# 求暴露弧上一点的切线
line_x = math.cos(theta) * rc + dgc
line_y = math.sin(theta) * rc + h_cav
k = math.tan(theta + math.pi / 2) # 入射角
# 求保护弧到直线的距离，判断是否相交
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)
global 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 # 入射角
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) * 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)
r_bd = 0
theta_sample, d_theta = np.linspace(
theta1, theta2, int((theta2 - theta1) / theta_fineness) + 1, retstep=True
)
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
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
)
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
# 把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]

303
main.py
View File

@ -1,279 +1,15 @@
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, h_gav, h_cav, dgc, 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, h_cav)
msp.add_circle((0, h_gav), rs, dxfattribs={"color": color})
msp.add_line((0, 0), (0, h_gav)) # 地线
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), (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) # 地线
# 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):
doc = self._doc
doc.saveas("egm.dxf")
# 圆交点
def solve_circle_intersection(rs, rc, hgav, hcav, dgc):
# 用牛顿法求解
x = rc # 初始值
y = rc # 初始值
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 []
# 圆与地面线交点
def solve_circle_line_intersection(radius, rg, center_y, center_x):
# TODO: 需要考虑地面捕雷线与暴露弧完全没交点的情况
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 * (math.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
def rg_fun(i_curt, h_cav):
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)
return rg
def intersection_angle(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)
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)
theta1_line = np_circle_line_intersection - np.array([dgc, h_cav])
theta1 = math.atan(theta1_line[1] / theta1_line[0])
# 考虑雷电入射角度，所以theta1可以小于0，即计算从侧面击中的雷
# if theta1 < 0:
# # print(f"θ_1角度为负数{theta1:.4f},人为设置为0")
# theta1 = 0
return np.array([theta1, theta2])
def distance_point_line(point_x, point_y, line_x, line_y, k):
d = abs(k * point_x - point_y - k * line_x + line_y) / ((k ** 2 + 1) ** 0.5)
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 calculus_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
k = math.tan(theta + math.pi / 2) # 入射角
# 求保护弧到直线的距离，判断是否相交
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)
global gCount
gCount += 1
if gCount % 1000 == 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, new_k * (-500 - line_x) + line_y),
(500, new_k * (500 - line_x) + line_y),
)
gMSP.add_line((0, rg), (1000, rg))
gCAD.save()
else:
max_w = theta + math.pi / 2 # 入射角
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)
r_bd = 0
theta_sample, d_theta = np.linspace(
theta1, theta2, int((theta2 - theta1) / theta_fineness) + 1, retstep=True
)
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
# r1=rc*(-math.cos(thyta2)+math.cos(thyta1))
# 入射角密度函数积分
# arrival_angle_fineness=0.0001
# for calculus_arv_angle in np.linspace()
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值不存在的情况
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
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
)
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
# 把k转化成相应的角度，从x开始，逆时针为正
k_angle = []
for kk in k_candidate:
if kk == 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]
from core import *
import timeit
def egm():
for u_bar in range(1):
u_ph = math.sqrt(2) * 750 * math.cos(2 * math.pi / 6 * 0) / 1.732 # 运行相电压
h_whole = 150 # 杆塔全高
h_gav = h_whole - 0.5 - 11.67 * 2 / 3
h_cav = h_gav - 9.5 - 2.7 - (14.43 - 11.67) * 2 / 3 # 导线对地平均高
dgc = 12 # 导地线水平距离
u_ph = math.sqrt(1) * 750 * math.cos(2 * math.pi / 3 * 0) / 1.732 # 运行相电压
h_whole = 140 # 杆塔全高
string_len = 6.8 # 串子绝缘长度
h_gav = h_whole - 0.5 - 11.67 * 2 / 3 # 地线对地平均高
h_cav = h_gav - 9.2 - 2.7 - (14.43 - 11.67) * 2 / 3 # 导线对地平均高
dgc = 0 # 导地线水平距离
# 迭代法计算最大电流
i_max = 0
_min_i = 20 # 尝试的最小电流
@ -285,10 +21,13 @@ def egm():
rg = rg_fun(i_bar, h_cav)
circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
if not circle_intersection: # if circle_intersection is []
# print("保护弧和暴露弧无交点，检查设置参数。程序退出。")
continue
circle_rc_line_intersection = solve_circle_line_intersection(
rc, rg, h_cav, dgc
rc, rg, dgc, h_cav
)
if not circle_rc_line_intersection:
continue
min_distance_intersection = (
np.sum(
(
@ -304,7 +43,7 @@ def egm():
break
if circle_intersection[1] < circle_rc_line_intersection[1]:
circle_rs_line_intersection = solve_circle_line_intersection(
rs, rg, h_gav, 0
rs, rg, 0, h_gav
)
# 判断与保护弧的交点是否在暴露弧外面
distance = (
@ -317,7 +56,7 @@ def egm():
if distance > rc:
print("暴露弧已经完全被屏蔽")
break
i_min = min_i(6.78, u_ph / 1.732)
i_min = min_i(string_len, 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)
@ -338,7 +77,7 @@ def egm():
i_curt_samples, d_curt = np.linspace(
i_min, i_max, curt_segment_n + 1, retstep=True
)
for i_curt in i_curt_samples:
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)
@ -351,13 +90,17 @@ def egm():
/ 2
* d_curt
)
n_sf = 2 * 2.7 / 10 * calculus # 调整率
n_sf = 2 * 2.7 / 10 * calculus # 跳闸率
print(f"跳闸率是{n_sf:.6}")
# draw(rs, rc, rg, h_gav, h_cav, dgc)
def speed():
a = 0
for bar in range(100000000):
a += bar
if __name__ == "__main__":
tangent_line_k(1, 0, 0, 0, 1)
egm()
run_time = timeit.timeit("egm()", globals=globals(), number=1)
print(f"运行时间:{run_time:.2f}s")
print("Finished.")