From 257d5bb23b9ec6050adf66ae1bc7e2017d31158d Mon Sep 17 00:00:00 2001 From: facat Date: Wed, 22 Sep 2021 11:22:13 +0800 Subject: [PATCH] =?UTF-8?q?=E4=BF=AE=E5=A4=8D=E5=87=A0=E4=B8=AA=E8=BE=B9?= =?UTF-8?q?=E7=95=8C=E6=9D=A1=E4=BB=B6bug?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- core.py | 59 +++++++++++++++++++++++++++++++++++++-------------------- main.py | 11 ++++++++--- 2 files changed, 46 insertions(+), 24 deletions(-) diff --git a/core.py b/core.py index 2532f56..ae0395f 100644 --- a/core.py +++ b/core.py @@ -158,24 +158,31 @@ def intersection_angle( circle_line_or_rg_intersection = solve_circle_intersection( rg, rc, rg_x, rg_y, rc_x, rc_y ) # 两圆的交点 - ( - circle_line_or_rg_intersection_x, - circle_line_or_rg_intersection_y, - ) = circle_line_or_rg_intersection - if ( - ground_surface(circle_line_or_rg_intersection_x) - > circle_line_or_rg_intersection_y - ): # 交点在地面线以下，就可以不积分 - # 找到暴露弧和地面线的交点 - circle_line_or_rg_intersection = circle_ground_surface_intersection( - rc, rc_x, rc_y, ground_surface - ) + if circle_line_or_rg_intersection: + ( + circle_line_or_rg_intersection_x, + circle_line_or_rg_intersection_y, + ) = circle_line_or_rg_intersection + if ( + ground_surface(circle_line_or_rg_intersection_x) + > circle_line_or_rg_intersection_y + ): # 交点在地面线以下，就可以不积分 + # 找到暴露弧和地面线的交点 + circle_line_or_rg_intersection = circle_ground_surface_intersection( + rc, rc_x, rc_y, ground_surface + ) + theta1 = None np_circle_intersection = np.array(circle_intersection) theta2_line = np_circle_intersection - np.array([rc_x, rc_y]) theta2 = math.atan(theta2_line[1] / theta2_line[0]) np_circle_line_or_rg_intersection = np.array(circle_line_or_rg_intersection) - theta1_line = np_circle_line_or_rg_intersection - np.array([rc_x, rc_y]) - theta1 = math.atan(theta1_line[1] / theta1_line[0]) + if not circle_line_or_rg_intersection: + if rc_y - rc > rg: # rg在rc下面 + # 捕捉线太低了，对高塔无保护，θ_1从-90°开始计算。 + theta1 = -math.pi / 2 + else: + theta1_line = np_circle_line_or_rg_intersection - np.array([rc_x, rc_y]) + theta1 = math.atan(theta1_line[1] / theta1_line[0]) return np.array([theta1, theta2]) @@ -187,7 +194,10 @@ def distance_point_line(point_x, point_y, line_x, line_y, k) -> float: def func_calculus_pw(theta, max_w): w_fineness = 0.01 - w_samples, d_w = np.linspace(0, max_w, int(max_w / w_fineness), retstep=True) + segments = int(max_w / w_fineness) + if segments < 2: # 最大最小太小，没有可以积分的 + return 0 + w_samples, d_w = np.linspace(0, max_w, segments, retstep=True) cal_w_np = abs(angel_density(w_samples)) * np.sin(theta - (w_samples - math.pi / 2)) r_pw = np.sum((cal_w_np[:-1] + cal_w_np[1:])) / 2 * d_w return r_pw @@ -282,7 +292,8 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0): for ind, k_cdi in enumerate(list(k_candidate)): k = k_candidate[ind] k_candidate[ind] = None - for bar in range(0, 30): + max_iteration = 30 + for bar in range(0, max_iteration): fk = (k * center_x - center_y - k * line_x + line_y) ** 2 - ( radius ** 2 ) * (k ** 2 + 1) @@ -304,13 +315,17 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0): if abs(dd - radius) < 1: k_candidate[ind] = k break + # 解决数值稳定性 + if bar == max_iteration - 1: + if abs(math.atan(k)) * 180 / math.pi > 89: + k_candidate[ind] = k # 把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: @@ -327,7 +342,9 @@ def func_ng(td): # 地闪密度 def circle_ground_surface_intersection(radius, center_x, center_y, ground_surface): # 最笨的办法，一个个去试 x_series = np.linspace(0, radius, int(radius / 0.001)) + center_x - part_to_be_squared = radius ** 2 - (x_series - center_x) ** 2 # 有可能出现-0.00001的数值，只是一个数值稳定问题。 + part_to_be_squared = ( + radius ** 2 - (x_series - center_x) ** 2 + ) # 有可能出现-0.00001的数值，只是一个数值稳定问题。 part_to_be_squared[ (part_to_be_squared < 0) & (abs(part_to_be_squared) < 1e-3) ] = 0 # 强制为0 diff --git a/main.py b/main.py index 9b6f6c8..3644709 100644 --- a/main.py +++ b/main.py @@ -1,3 +1,5 @@ +import numpy as np + from core import * import timeit @@ -5,13 +7,13 @@ import timeit def egm(): h_g_avr_sag = 11.67 * 2 / 3 h_c_avr_sag = 14.43 * 2 / 3 - h_whole = 140 # 杆塔全高 + h_whole = 260 # 杆塔全高 voltage_n = 3 # 工作电压分成多少份来计算 td = 20 # 雷暴日 insulator_c_len = 6.8 # 串子绝缘长度 string_c_len = 9.2 string_g_len = 0.5 - gc_x = [17.9, 16, 15, 16] + gc_x = [17.9, 17, 15, 17] # 以后考虑地形角度,地面线 def ground_surface(x): @@ -34,6 +36,7 @@ def egm(): rg_x = None rg_y = None cad = Draw() + 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] @@ -149,7 +152,7 @@ def egm(): if distance > rc: print("暴露弧已经完全被屏蔽") break - if phase_conductor == 1: + 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() @@ -205,8 +208,10 @@ def egm(): 2 * ng / 10 * calculus ) # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13 avr_n_sf += n_sf / voltage_n + n_sf_phases[phase_conductor][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)}") def speed():