diff --git a/core.py b/core.py index ae0395f..db93b48 100644 --- a/core.py +++ b/core.py @@ -1,6 +1,5 @@ import math import ezdxf -import numba import numpy as np gCAD = None @@ -47,9 +46,10 @@ class Draw: 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} - ) # 圆和圆的交点 + if rg_rc_intersection: + msp.add_line( + (rc_x, rc_y), rg_rc_intersection, dxfattribs={"color": color} + ) # 圆和圆的交点 # 计算圆交点 # msp.add_line((dgc, h_cav), circle_intersection) # 导线 @@ -58,6 +58,10 @@ class Draw: doc = self._doc doc.saveas("egm.dxf") + def saveas(self, file_name): + doc = self._doc + doc.saveas(file_name) + # 圆交点 def solve_circle_intersection( @@ -69,8 +73,8 @@ def solve_circle_intersection( center_y2, ): # 用牛顿法求解 - x = radius2 # 初始值 - y = radius2 # 初始值 + x = radius2 + center_x2 # 初始值 + y = radius2 + center_y2 # 初始值 # TODO 考虑出现2个解的情况 for bar in range(0, 10): A = [ @@ -89,7 +93,7 @@ def solve_circle_intersection( 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: @@ -259,9 +263,10 @@ def bd_area( rc = rc_fun(i_curt, u_ph) rs = rs_fun(i_curt) rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type) - theta_sample, d_theta = np.linspace( - theta1, theta2, int((theta2 - theta1) / theta_fineness), retstep=True - ) + theta_segments = int((theta2 - theta1) / theta_fineness) + if theta_segments < 2: + return 0 + theta_sample, d_theta = np.linspace(theta1, theta2, theta_segments, retstep=True) if len(theta_sample) < 2: return 0 vec_calculus_bd = np.vectorize(calculus_bd) @@ -322,10 +327,10 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0): # 把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 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: @@ -354,3 +359,11 @@ def circle_ground_surface_intersection(radius, center_x, center_y, ground_surfac r_x = x_series[equal_location][0] r_y = ground_surface(r_x) return r_x, r_y + + +# u_ph是相电压 +# insulator_c_len绝缘子闪络距离 +def arc_possibility(rated_voltage, insulator_c_len): # 建弧率 + _e = rated_voltage / (3 ** 0.5) / insulator_c_len + r = (4.5 * (_e ** 0.75) - 14) * 1e-2 + return r diff --git a/main.py b/main.py index 3644709..fd4e261 100644 --- a/main.py +++ b/main.py @@ -7,7 +7,7 @@ import timeit def egm(): h_g_avr_sag = 11.67 * 2 / 3 h_c_avr_sag = 14.43 * 2 / 3 - h_whole = 260 # 杆塔全高 + h_whole = 40 # 杆塔全高 voltage_n = 3 # 工作电压分成多少份来计算 td = 20 # 雷暴日 insulator_c_len = 6.8 # 串子绝缘长度 @@ -36,33 +36,38 @@ def egm(): rg_x = None rg_y = None cad = Draw() + # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13 + ng = func_ng(td) n_sf_phases = np.zeros((phase_n, voltage_n)) # 计算每一相的跳闸率 - for phase_conductor in range(phase_n): - rs_x = gc_x[phase_conductor] - rs_y = gc_y[phase_conductor] - rc_x = gc_x[phase_conductor + 1] - rc_y = gc_y[phase_conductor + 1] + if np.any(np.array(gc_y) < 0): + print("导线可能掉地面了，程序退出。") + return 0 + for phase_conductor_foo in range(phase_n): + exposed_curve_shielded = False + rs_x = gc_x[phase_conductor_foo] + rs_y = gc_y[phase_conductor_foo] + rc_x = gc_x[phase_conductor_foo + 1] + rc_y = gc_y[phase_conductor_foo + 1] if phase_n == 1: rg_type = "g" if phase_n > 1: # 多回路 - if phase_conductor < 2: + if phase_conductor_foo < 2: rg_type = "c" - rg_x = gc_x[phase_conductor + 2] - rg_y = gc_y[phase_conductor + 2] + rg_x = gc_x[phase_conductor_foo + 2] + rg_y = gc_y[phase_conductor_foo + 2] else: rg_type = "g" # TODO 保护角公式可能有问题，后面改 shield_angle = ( - math.atan(rc_x / ((rc_y - rs_y) + string_c_len)) * 180 / math.pi + math.atan((rc_x - rs_x) / ((rs_y - rc_y) + string_c_len)) * 180 / math.pi ) # 保护角 print(f"保护角{shield_angle:.3f}°") print(f"最低相防护标识{rg_type}") - ng = func_ng(td) for u_bar in range(voltage_n): u_ph = ( math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732 ) # 运行相电压 - print(f"计算第{phase_conductor + 1}相，电压为{u_ph:.2f}kV") + print(f"计算第{phase_conductor_foo + 1}相，电压为{u_ph:.2f}kV") # 迭代法计算最大电流 i_max = 0 i_min = min_i(insulator_c_len, u_ph / 1.732) @@ -151,17 +156,21 @@ def egm(): ) if distance > rc: print("暴露弧已经完全被屏蔽") + exposed_curve_shielded = True break - if phase_conductor == 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, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 6) - cad.save() + # if phase_conductor_foo == 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, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 6) + cad.saveas(f"egm{phase_conductor_foo+1}.dxf") # 判断是否导线已经被完全保护 if abs(i_max - _max_i) < 1e-5: print("无法找到最大电流，可能是杆塔较高。") print(f"最大电流设置为自然界最大电流{i_max}kA") print(f"最大电流为{i_max:.2f}") print(f"最小电流为{i_min:.2f}") + if exposed_curve_shielded: + print("暴露弧已经完全被屏蔽,不会跳闸。") + continue curt_fineness = 0.1 # 电流积分细度 if i_min > i_max or abs(i_min - i_max) < curt_fineness: print("最大电流小于最小电流，没有暴露弧。") @@ -204,11 +213,9 @@ def egm(): # if abs(calculus-0.05812740052770032)<1e-5: # abc=123 # pass - n_sf = ( - 2 * ng / 10 * calculus - ) # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13 + n_sf = (2 * ng / 10 * calculus) * arc_possibility(750, insulator_c_len) avr_n_sf += n_sf / voltage_n - n_sf_phases[phase_conductor][u_bar] = n_sf + n_sf_phases[phase_conductor_foo][u_bar] = n_sf print(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}") print(f"跳闸率是{avr_n_sf:.6f}") print(f"不同相跳闸率是{np.mean(n_sf_phases,axis=1)}")