diff --git a/main.py b/main.py index 9a0f578..d563fba 100644 --- a/main.py +++ b/main.py @@ -3,6 +3,35 @@ import ezdxf import numpy as np +class Draw: + def __init__(self): + self._doc = ezdxf.new(dxfversion="R2010") + self._doc.layers.add("EGM", color=2) + + def draw(self, i_curt, u_ph, h_gav, h_cav, dgc): + doc = self._doc + msp = doc.modelspace() + 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) + 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) # 导线和圆的交点 + + def save(self): + doc = self._doc + doc.saveas("egm.dxf") + + # 圆交点 def solve_circle_intersection(rs, rc, hgav, hcav, dgc): # x = Symbol('x', real=True) @@ -45,26 +74,6 @@ def solve_circle_line_intersection(rc, rg, hcav, 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 # 雷电波阻抗 @@ -74,7 +83,7 @@ def min_i(string_len, u_ph): def thunder_density(i): # l雷电流幅值密度函数 - r = -10 ** (-i / 44) * math.log(10) * (-1 / 44) + r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44) return r @@ -113,10 +122,14 @@ def intersection_angel(dgc, h_gav, h_cav, i_curt, u_ph): # 暴露弧的角度 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]) + if theta1 < 0: + # print(f"θ_1角度为负数{theta1:.4f},人为设置为0") + theta1 = 0 return np.array([theta1, theta2]) -def bd_area(i_curt, u_ph, theta1, theta2): # 暴露弧的投影面积 +def bd_area(i_curt, u_ph, dgc, h_gav, h_cav): # 暴露弧的投影面积 + theta1, theta2 = intersection_angel(dgc, h_gav, h_cav, i_curt, u_ph) rc = rc_fun(i_curt, u_ph) # 暂时不考虑雷电入射角的影响 r = (math.cos(theta1) - math.cos(theta2)) * rc @@ -129,12 +142,12 @@ def bd_area(i_curt, u_ph, theta1, theta2): # 暴露弧的投影面积 def egm(): u_ph = 750 / 1.732 # 运行相电压 - h_cav = 60 # 导线对地平均高 - h_gav = h_cav + 9.5 + 7.2 - dgc = 2 + h_cav = 160 # 导线对地平均高 + h_gav = h_cav + 9.5 + 2.2 + dgc = 2 # 导地线水平距离 # 迭代法计算最大电流 i_max = 0 - _min_i = 30 # 尝试的最小电流 + _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}") @@ -158,27 +171,33 @@ def egm(): ) ** 0.5 ) # 计算两圆交点和地面直线交点的最小距离 - if min_distance_intersection < 0.01: - i_max = i_bar - draw(rs, rc, rg, h_gav, h_cav, dgc) + i_max = i_bar + if min_distance_intersection < 0.1: break - print(f"最大电流为{i_max:.2f}") i_min = min_i(6.78, 750 / 1.732) + cad = Draw() + cad.draw(i_min, u_ph, h_gav, h_cav, dgc) + cad.draw(i_max, u_ph, h_gav, h_cav, dgc) + 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) - d_curt = (i_max - i_min) / curt_segment_n + curt_fineness = 0.001 # 电流积分细度 + curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份 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) + 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 @@ -187,8 +206,8 @@ def egm(): / 2 * d_curt ) - n_sf=2*2.7/10*calculus - print(f'跳闸率是{n_sf}') + n_sf = 2 * 2.7 / 10 * calculus # 调整率 + print(f"跳闸率是{n_sf:.6}") # draw(rs, rc, rg, h_gav, h_cav, dgc)