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

初步完成了双回路公式

This commit is contained in:
facat 2021-09-21 20:00:03 +08:00
parent 5a75df4542
commit dd44de030e
2 changed files with 140 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

128
core.py
View File

@ -15,26 +15,41 @@ class Draw:
global gCAD global gCAD
gCAD = self gCAD = self
def draw(self, i_curt, u_ph, h_gav, h_cav, dgc, color): 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 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, h_cav) rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
msp.add_circle((0, h_gav), rs, dxfattribs={"color": color}) msp.add_circle((rs_x, rs_y), rs, dxfattribs={"color": color})
msp.add_line((0, 0), (0, h_gav)) # 地线 msp.add_line((0, 0), (rs_x, rs_y)) # 地线
msp.add_circle((dgc, h_cav), rc, dxfattribs={"color": color}) msp.add_circle((rc_x, rc_y), rc, dxfattribs={"color": color+2})
msp.add_line((dgc, 0), (dgc, h_cav)) # 导线 msp.add_line((rc_x, 0), (rc_x, rc_y)) # 导线
msp.add_line((0, h_gav), (dgc, h_cav)) 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}) msp.add_line((0, rg), (2000, rg), dxfattribs={"color": color})
circle_line_section = solve_circle_line_intersection(rc, rg, rc_x, rc_y)
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
@ -42,15 +57,26 @@ class Draw:
# 圆交点 # 圆交点
def solve_circle_intersection(rs, rc, h_gav, h_cav, dgc): def solve_circle_intersection(
radius1,
radius2,
center_x1,
center_y1,
center_x2,
center_y2,
):
# 用牛顿法求解 # 用牛顿法求解
x = rc # 初始值 x = radius2 # 初始值
y = rc # 初始值 y = radius2 # 初始值
# TODO 考虑出现2个解的情况
for bar in range(0, 10): for bar in range(0, 10):
A = [[-2 * x, -2 * (y - h_gav)], [-2 * (x - dgc), -2 * (y - h_cav)]] A = [
[-2 * (x - center_x1), -2 * (y - center_y1)],
[-2 * (x - center_x2), -2 * (y - center_y2)],
]
b = [ b = [
x ** 2 + (y - h_gav) ** 2 - rs ** 2, (x - center_x1) ** 2 + (y - center_y1) ** 2 - radius1 ** 2,
(x - dgc) ** 2 + (y - h_cav) ** 2 - rc ** 2, (x - center_x2) ** 2 + (y - center_y2) ** 2 - radius2 ** 2,
] ]
X_set = np.linalg.solve(A, b) X_set = np.linalg.solve(A, b)
x += X_set[0] x += X_set[0]
@ -84,7 +110,7 @@ def thunder_density(i): # l雷电流幅值密度函数
def angel_density(angle): # 入射角密度函数 angle单位是弧度 def angel_density(angle): # 入射角密度函数 angle单位是弧度
r = 0.75 * (np.cos(angle - math.pi / 2) ** 3) r = 0.75 * abs((np.cos(angle - math.pi / 2) ** 3))
return r return r
@ -95,35 +121,50 @@ 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
def rg_fun(i_curt, h_cav): # typ 如果是g，代表捕雷线公式，c代表暴露弧公式
def rg_fun(i_curt, h_cav, u_ph, typ="g"):
rg = None
if typ == "g":
if h_cav < 40: if h_cav < 40:
rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65) rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65)
else: else:
rg = 5.5 * (i_curt ** 0.65) rg = 5.5 * (i_curt ** 0.65)
elif typ == "c": # 此时返回的是圆半径
rg = rc_fun(i_curt, u_ph)
return rg return rg
def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph): # 暴露弧的角度 def intersection_angle(
rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, 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, h_cav) rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc) # 两圆的交点 circle_intersection = solve_circle_intersection(
circle_line_intersection = solve_circle_line_intersection( rs, rc, rs_x, rs_y, rc_x, rc_y
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
) # 两圆的交点
np_circle_intersection = np.array(circle_intersection) np_circle_intersection = np.array(circle_intersection)
theta2_line = np_circle_intersection - np.array([dgc, h_cav]) theta2_line = np_circle_intersection - np.array([rc_x, rc_y])
theta2 = math.atan(theta2_line[1] / theta2_line[0]) theta2 = math.atan(theta2_line[1] / theta2_line[0])
np_circle_line_intersection = np.array(circle_line_intersection) np_circle_line_or_rg_intersection = np.array(circle_line_or_rg_intersection)
theta1_line = np_circle_line_intersection - np.array([dgc, h_cav]) theta1_line = np_circle_line_or_rg_intersection - np.array([rc_x, rc_y])
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
@ -137,17 +178,17 @@ def func_calculus_pw(theta, max_w):
return r_pw return r_pw
def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): # 对θ进行积分 def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y): # 对θ进行积分
max_w = 0 max_w = 0
# 求暴露弧上一点的切线 # 求暴露弧上一点的切线
line_x = math.cos(theta) * rc + dgc line_x = math.cos(theta) * rc + rc_x
line_y = math.sin(theta) * rc + h_cav line_y = math.sin(theta) * rc + rc_y
k = math.tan(theta + math.pi / 2) # 入射角 k = math.tan(theta + math.pi / 2) # 入射角
# 求保护弧到直线的距离，判断是否相交 # 求保护弧到直线的距离，判断是否相交
d_to_rs = distance_point_line(0, h_gav, line_x, line_y, k) d_to_rs = distance_point_line(rs_x, rs_y, line_x, line_y, k)
if d_to_rs < rs: # 相交 if d_to_rs < rs: # 相交
# 要用过直线上一点到暴露弧的切线 # 要用过直线上一点到暴露弧的切线
new_k = tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k) new_k = tangent_line_k(line_x, line_y, rs_x, rs_y, 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:
@ -183,20 +224,21 @@ def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): # 对θ进行积分
return r return r
def bd_area(i_curt, u_ph, dgc, h_gav, h_cav): # 暴露弧的投影面积 def bd_area(i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, 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, 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, h_cav) rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
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, dgc, h_cav, h_gav) 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 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 += (
@ -214,7 +256,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不存在
@ -248,10 +290,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:

73
main.py
View File

@ -6,20 +6,26 @@ import timeit
def egm(): def egm():
avr_n_sf = 0 # 考虑电压的影响计算的跳闸率 avr_n_sf = 0 # 考虑电压的影响计算的跳闸率
voltage_n = 3 # 工作电压分成多少份来计算 voltage_n = 1 # 工作电压分成多少份来计算
ng = func_ng(20) ng = func_ng(20)
h_whole = 140 # 杆塔全高 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
dgc = -0.0 # 导地线水平距离 rc_x = -0.0 # 导地线水平距离
rs_x = 0
rg_x = 0
vertical_dgc = 2.7 # 导地线挂点垂直距离 vertical_dgc = 2.7 # 导地线挂点垂直距离
h_g_avr_sag = 11.67 * 2 / 3 h_g_avr_sag = 11.67 * 2 / 3
h_c_avr_sag = 14.43 * 2 / 3 h_c_avr_sag = 14.43 * 2 / 3
h_gav = h_whole - string_g_len - h_g_avr_sag # 地线对地平均高 rs_y = h_whole - string_g_len - h_g_avr_sag # 地线对地平均高
h_cav = h_whole - string_c_len - vertical_dgc - h_c_avr_sag # 导线对地平均高 rc_y = h_whole - string_c_len - vertical_dgc - h_c_avr_sag # 导线对地平均高
shield_angle = math.atan(dgc / (vertical_dgc + string_c_len)) * 180 / math.pi rg_y = rc_y - 20
shield_angle = (
math.atan(rc_x / (vertical_dgc + string_c_len)) * 180 / math.pi
) # 保护角
print(f"保护角{shield_angle:.3f}°") print(f"保护角{shield_angle:.3f}°")
rg_type = "c"
for u_bar in range(voltage_n): for u_bar in range(voltage_n):
u_ph = ( u_ph = (
math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732 math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732
@ -30,12 +36,12 @@ def egm():
_min_i = i_min # 尝试的最小电流 _min_i = i_min # 尝试的最小电流
_max_i = 200 # 尝试的最大电流 _max_i = 200 # 尝试的最大电流
cad = Draw() cad = Draw()
cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2) # 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)): # 雷电流 for i_bar in np.linspace(_min_i, _max_i, int((_max_i - _min_i) / 0.1)): # 雷电流
# print(f"尝试计算电流为{i_bar:.2f}") print(f"尝试计算电流为{i_bar:.2f}")
rs = rs_fun(i_bar) rs = rs_fun(i_bar)
rc = rc_fun(i_bar, u_ph) rc = rc_fun(i_bar, u_ph)
rg = rg_fun(i_bar, h_cav) rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type)
####### #######
# cccCount += 1 # cccCount += 1
# if cccCount % 30 == 0: # if cccCount % 30 == 0:
@ -47,20 +53,28 @@ def egm():
# ) # )
# core.gMSP.add_circle((dgc, h_cav), rc) # core.gMSP.add_circle((dgc, h_cav), rc)
####### #######
circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc) rg_rc_circle_intersection = solve_circle_intersection(
if not circle_intersection: # if circle_intersection is [] rs, rc, rs_x, rs_y, rc_x, rc_y
# print("保护弧和暴露弧无交点，检查设置参数。程序退出。")
continue
circle_rc_line_intersection = solve_circle_line_intersection(
rc, rg, dgc, h_cav
) )
if not circle_rc_line_intersection: if not rg_rc_circle_intersection: # if circle_intersection is []
print("保护弧和暴露弧无交点，检查设置参数。程序退出。")
continue
circle_rc_line_or_rg_intersection = None
if rg_type == "g":
circle_rc_line_or_rg_intersection = solve_circle_line_intersection(
rc, rg, rc_x, rc_y
)
elif rg_type == "c":
circle_rc_line_or_rg_intersection = solve_circle_intersection(
rg, rc, rg_x, rg_y, rc_x, rc_y
)
if not circle_rc_line_or_rg_intersection:
continue continue
min_distance_intersection = ( min_distance_intersection = (
np.sum( np.sum(
( (
np.array(circle_intersection) np.array(rg_rc_circle_intersection)
- np.array(circle_rc_line_intersection) - np.array(circle_rc_line_or_rg_intersection)
) )
** 2 ** 2
) )
@ -69,14 +83,24 @@ def egm():
i_max = i_bar i_max = i_bar
if min_distance_intersection < 0.1: if min_distance_intersection < 0.1:
break break
if circle_intersection[1] < circle_rc_line_intersection[1]: # 判断是否以完全被保护
circle_rs_line_intersection = solve_circle_line_intersection( if rg_rc_circle_intersection[1] < circle_rc_line_or_rg_intersection[1]:
rs, rg, 0, h_gav 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 = ( distance = (
np.sum( np.sum(
(np.array(circle_rs_line_intersection) - np.array([dgc, h_cav])) (
np.array(circle_rs_line_or_rg_intersection)
- np.array([rc_x, rc_y])
)
** 2 ** 2
) )
** 0.5 ** 0.5
@ -84,8 +108,8 @@ def egm():
if distance > rc: if distance > rc:
print("暴露弧已经完全被屏蔽") print("暴露弧已经完全被屏蔽")
break break
cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2) 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, h_gav, h_cav, dgc, 6) cad.draw(i_max, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 6)
cad.save() cad.save()
# 判断是否导线已经被完全保护 # 判断是否导线已经被完全保护
if abs(i_max - _max_i) < 1e-5: if abs(i_max - _max_i) < 1e-5:
@ -99,13 +123,12 @@ def egm():
return return
# 开始积分 # 开始积分
curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份 curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份
calculus = 0
i_curt_samples, d_curt = np.linspace( i_curt_samples, d_curt = np.linspace(
i_min, i_max, curt_segment_n + 1, retstep=True i_min, i_max, curt_segment_n + 1, retstep=True
) )
bd_area_vec = np.vectorize(bd_area) bd_area_vec = np.vectorize(bd_area)
cal_bd_np = bd_area_vec( cal_bd_np = bd_area_vec(
i_curt_samples, u_ph, dgc, h_gav, h_cav i_curt_samples, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, rg_type
) * thunder_density(i_curt_samples) ) * thunder_density(i_curt_samples)
calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt
# for i_curt in i_curt_samples[:-1]: # for i_curt in i_curt_samples[:-1]: