From 392eeb01688b5c3bdefa06b046066fb8ffbdf1e5 Mon Sep 17 00:00:00 2001 From: facat Date: Wed, 22 Sep 2021 00:18:06 +0800 Subject: [PATCH] =?UTF-8?q?=E5=AE=8C=E5=96=84=E4=BA=86=E5=8F=8C=E5=9B=9E?= =?UTF-8?q?=E8=B7=AFEGM=E6=A8=A1=E5=9E=8B=E4=BB=A3=E7=A0=81=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- core.py | 47 ++++++-- main.py | 324 +++++++++++++++++++++++++++++++++----------------------- 2 files changed, 229 insertions(+), 142 deletions(-) diff --git a/core.py b/core.py index 7631be5..2532f56 100644 --- a/core.py +++ b/core.py @@ -25,7 +25,7 @@ class Draw: rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type) msp.add_circle((rs_x, rs_y), rs, dxfattribs={"color": color}) msp.add_line((0, 0), (rs_x, rs_y)) # 地线 - msp.add_circle((rc_x, rc_y), rc, dxfattribs={"color": color+2}) + msp.add_circle((rc_x, rc_y), rc, dxfattribs={"color": color + 2}) msp.add_line((rc_x, 0), (rc_x, rc_y)) # 导线 msp.add_line((rs_x, rs_y), (rc_x, rc_y)) # 角度线 @@ -36,9 +36,12 @@ class Draw: if rg_type == "g": 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 not circle_line_section: + pass + else: + 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( @@ -138,7 +141,7 @@ def rg_fun(i_curt, h_cav, u_ph, typ="g"): def intersection_angle( - rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, rg_type + rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_surface, rg_type ): # 暴露弧的角度 rs = rs_fun(i_curt) rc = rc_fun(i_curt, u_ph) @@ -155,6 +158,18 @@ 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 + ) 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]) @@ -224,9 +239,11 @@ def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y): # 对θ进行积分 return r -def bd_area(i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, rg_type): # 暴露弧的投影面积 +def bd_area( + i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, ground_surface, rg_type +): # 暴露弧的投影面积 theta1, theta2 = intersection_angle( - rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, rg_type + rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_surface, rg_type ) # θ角度 theta_fineness = 0.01 rc = rc_fun(i_curt, u_ph) @@ -304,3 +321,19 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0): def func_ng(td): # 地闪密度 return 0.023 * (td ** 1.3) + + +# 圆和地面线的交点，只去正x轴上的。 +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[ + (part_to_be_squared < 0) & (abs(part_to_be_squared) < 1e-3) + ] = 0 # 强制为0 + y_series = center_y - part_to_be_squared ** 0.5 + ground_surface_y = ground_surface(x_series) + equal_location = np.abs(ground_surface_y - y_series) < 0.5 + r_x = x_series[equal_location][0] + r_y = ground_surface(r_x) + return r_x, r_y diff --git a/main.py b/main.py index 5d958c5..9b6f6c8 100644 --- a/main.py +++ b/main.py @@ -1,158 +1,212 @@ -import numpy as np - from core import * import timeit def egm(): - avr_n_sf = 0 # 考虑电压的影响计算的跳闸率 - voltage_n = 1 # 工作电压分成多少份来计算 - ng = func_ng(20) + h_g_avr_sag = 11.67 * 2 / 3 + h_c_avr_sag = 14.43 * 2 / 3 h_whole = 140 # 杆塔全高 + voltage_n = 3 # 工作电压分成多少份来计算 + td = 20 # 雷暴日 insulator_c_len = 6.8 # 串子绝缘长度 string_c_len = 9.2 string_g_len = 0.5 - rc_x = -0.0 # 导地线水平距离 - rs_x = 0 - rg_x = 0 - vertical_dgc = 2.7 # 导地线挂点垂直距离 - h_g_avr_sag = 11.67 * 2 / 3 - h_c_avr_sag = 14.43 * 2 / 3 - rs_y = h_whole - string_g_len - h_g_avr_sag # 地线对地平均高 - rc_y = h_whole - string_c_len - vertical_dgc - h_c_avr_sag # 导线对地平均高 - rg_y = rc_y - 20 - shield_angle = ( - math.atan(rc_x / (vertical_dgc + string_c_len)) * 180 / math.pi - ) # 保护角 - print(f"保护角{shield_angle:.3f}°") - rg_type = "c" - for u_bar in range(voltage_n): - u_ph = ( - math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732 - ) # 运行相电压 - # 迭代法计算最大电流 - i_max = 0 - i_min = min_i(insulator_c_len, u_ph / 1.732) - _min_i = i_min # 尝试的最小电流 - _max_i = 200 # 尝试的最大电流 - cad = Draw() - # 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)): # 雷电流 - print(f"尝试计算电流为{i_bar:.2f}") - rs = rs_fun(i_bar) - rc = rc_fun(i_bar, u_ph) - rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type) - ####### - # cccCount += 1 - # if cccCount % 30 == 0: - # import core - # - # core.gMSP.add_circle((0, h_gav), rs) - # core.gMSP.add_circle( - # (dgc, h_cav), rc_fun(i_bar, -u_ph), dxfattribs={"color": 4} - # ) - # core.gMSP.add_circle((dgc, h_cav), rc) - ####### - rg_rc_circle_intersection = solve_circle_intersection( - rs, rc, rs_x, rs_y, rc_x, rc_y - ) - 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 + gc_x = [17.9, 16, 15, 16] + + # 以后考虑地形角度,地面线 + def ground_surface(x): + return 0 + + gc_y = [ + h_whole - string_g_len - h_g_avr_sag, # 地线对地平均高 + h_whole - string_c_len - h_c_avr_sag - 2.7, # 导线对地平均高 + h_whole - string_c_len - h_c_avr_sag - 20, # 导线对地平均高 + h_whole - string_c_len - h_c_avr_sag - 35.7, # 导线对地平均高 + ] + if len(gc_y) > 2: # 双回路 + phase_n = 3 # 边相导线数量 + else: + phase_n = 1 + ######################################################### + rg_type = None + # 以上是需要设置的参数 + avr_n_sf = 0 # 考虑电压的影响计算的跳闸率 + rg_x = None + rg_y = None + cad = Draw() + 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 phase_n == 1: + rg_type = "g" + if phase_n > 1: # 多回路 + if phase_conductor < 2: + rg_type = "c" + rg_x = gc_x[phase_conductor + 2] + rg_y = gc_y[phase_conductor + 2] + else: + rg_type = "g" + # TODO 保护角公式可能有问题，后面改 + shield_angle = ( + math.atan(rc_x / ((rc_y - rs_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") + # 迭代法计算最大电流 + i_max = 0 + i_min = min_i(insulator_c_len, u_ph / 1.732) + _min_i = i_min # 尝试的最小电流 + _max_i = 200 # 尝试的最大电流 + # 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) + ): # 雷电流 + # print(f"尝试计算电流为{i_bar:.2f}") + rs = rs_fun(i_bar) + rc = rc_fun(i_bar, u_ph) + rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type) + ####### + # cccCount += 1 + # if cccCount % 30 == 0: + # import core + # + # core.gMSP.add_circle((0, h_gav), rs) + # core.gMSP.add_circle( + # (dgc, h_cav), rc_fun(i_bar, -u_ph), dxfattribs={"color": 4} + # ) + # core.gMSP.add_circle((dgc, h_cav), rc) + ####### + rg_rc_circle_intersection = solve_circle_intersection( + rs, rc, rs_x, rs_y, 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 - min_distance_intersection = ( - np.sum( - ( - np.array(rg_rc_circle_intersection) - - np.array(circle_rc_line_or_rg_intersection) - ) - ** 2 - ) - ** 0.5 - ) # 计算两圆交点和地面直线交点的最小距离 - i_max = i_bar - if min_distance_intersection < 0.1: - break - # 判断是否以完全被保护 - if rg_rc_circle_intersection[1] < circle_rc_line_or_rg_intersection[1]: - circle_rs_line_or_rg_intersection = None + i_max = i_bar + if not rg_rc_circle_intersection: # if circle_intersection is [] + print("保护弧和暴露弧无交点，检查设置参数。") + continue + circle_rc_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 + circle_rc_line_or_rg_intersection = solve_circle_line_intersection( + rc, rg, rc_x, rc_y ) - if rg_type == "c": - circle_rs_line_or_rg_intersection = solve_circle_intersection( - rs, rg, rs_x, rs_y, rg_x, rg_y + elif rg_type == "c": + circle_rc_line_or_rg_intersection = solve_circle_intersection( + rg, rc, rg_x, rg_y, rc_x, rc_y ) - # 判断与保护弧的交点是否在暴露弧外面 - distance = ( + if not circle_rc_line_or_rg_intersection: + # 暴露弧和捕捉弧无交点 + if rg_type == "g": + if rg > rc_y: + i_min = i_bar + print(f"捕捉弧在暴露弧之上，设置最小电流为{i_min:.2f}") + else: + print("暴露弧和捕捉弧无交点，检查设置参数。") + continue + else: + print("暴露弧和捕捉弧无交点，检查设置参数。") + continue + min_distance_intersection = ( np.sum( ( - np.array(circle_rs_line_or_rg_intersection) - - np.array([rc_x, rc_y]) + np.array(rg_rc_circle_intersection) + - np.array(circle_rc_line_or_rg_intersection) ) ** 2 ) ** 0.5 - ) - if distance > rc: - print("暴露弧已经完全被屏蔽") + ) # 计算两圆交点和地面直线交点的最小距离 + if min_distance_intersection < 0.1: break - 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 abs(i_max - _max_i) < 1e-5: - print("无法找到最大电流，可能是杆塔较高。") - print(f"最大电流设置为自然界最大电流{i_max}kA") - print(f"最大电流为{i_max:.2f}") - print(f"最小电流为{i_min:.2f}") - curt_fineness = 0.1 # 电流积分细度 - if i_min > i_max or abs(i_min - i_max) < curt_fineness: - print("最大电流小于最小电流，没有暴露弧，程序结束。") - return - # 开始积分 - curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份 - i_curt_samples, d_curt = np.linspace( - i_min, i_max, curt_segment_n + 1, retstep=True - ) - bd_area_vec = np.vectorize(bd_area) - cal_bd_np = bd_area_vec( - i_curt_samples, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, rg_type - ) * thunder_density(i_curt_samples) - calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt - # 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) - # cal_thunder_density_second = thunder_density(i_curt + d_curt) - # calculus += ( - # ( - # cal_bd_first * cal_thunder_density_first - # + cal_bd_second * cal_thunder_density_second - # ) - # / 2 - # * d_curt - # ) - # 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 - avr_n_sf += n_sf / voltage_n - print(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}") - print(f"跳闸率是{avr_n_sf:.6}") + # 判断是否以完全被保护 + if rg_rc_circle_intersection[1] < circle_rc_line_or_rg_intersection[1]: + 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 = ( + np.sum( + ( + np.array(circle_rs_line_or_rg_intersection) + - np.array([rc_x, rc_y]) + ) + ** 2 + ) + ** 0.5 + ) + if distance > rc: + print("暴露弧已经完全被屏蔽") + break + if phase_conductor == 1: + 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 abs(i_max - _max_i) < 1e-5: + print("无法找到最大电流，可能是杆塔较高。") + print(f"最大电流设置为自然界最大电流{i_max}kA") + print(f"最大电流为{i_max:.2f}") + print(f"最小电流为{i_min:.2f}") + curt_fineness = 0.1 # 电流积分细度 + if i_min > i_max or abs(i_min - i_max) < curt_fineness: + print("最大电流小于最小电流，没有暴露弧。") + continue + # 开始积分 + curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份 + i_curt_samples, d_curt = np.linspace( + i_min, i_max, curt_segment_n + 1, retstep=True + ) + bd_area_vec = np.vectorize(bd_area) + cal_bd_np = ( + bd_area_vec( + i_curt_samples, + u_ph, + rc_x, + rc_y, + rs_x, + rs_y, + rg_x, + rg_y, + ground_surface, + rg_type, + ) + * thunder_density(i_curt_samples) + ) + calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt + # 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) + # cal_thunder_density_second = thunder_density(i_curt + d_curt) + # calculus += ( + # ( + # cal_bd_first * cal_thunder_density_first + # + cal_bd_second * cal_thunder_density_second + # ) + # / 2 + # * d_curt + # ) + # 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 + avr_n_sf += n_sf / voltage_n + print(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}") + print(f"跳闸率是{avr_n_sf:.6f}") def speed():