diff --git a/core.py b/core.py index 7ad9a0d..31d6119 100644 --- a/core.py +++ b/core.py @@ -1,6 +1,7 @@ import math import ezdxf import numpy as np +from typing import List gCAD = None gMSP = None @@ -15,8 +16,8 @@ class Parameter: insulator_c_len: float # 串子绝缘长度 string_c_len: float string_g_len: float - gc_x: [float] # 导、地线水平坐标 - ground_angels: [float] # 地面倾角，向下为正 + gc_x: List[float] # 导、地线水平坐标 + ground_angels: List[float] # 地面倾角，向下为正 h_arm: float # 导、地线垂直坐标 altitude: int # 海拔，单位米 max_i: float # 最大尝试电流，单位kA @@ -131,14 +132,14 @@ def solve_circle_intersection( x = radius2 + center_x2 # 初始值 y = radius2 + center_y2 # 初始值 # TODO 考虑出现2个解的情况 - for bar in range(0, 10): + for _ in range(0, 10): A = [ [-2 * (x - center_x1), -2 * (y - center_y1)], [-2 * (x - center_x2), -2 * (y - center_y2)], ] b = [ - (x - center_x1) ** 2 + (y - center_y1) ** 2 - radius1 ** 2, - (x - center_x2) ** 2 + (y - center_y2) ** 2 - radius2 ** 2, + (x - center_x1) ** 2 + (y - center_y1) ** 2 - radius1**2, + (x - center_x2) ** 2 + (y - center_y2) ** 2 - radius2**2, ] X_set = np.linalg.solve(A, b) x += X_set[0] @@ -166,22 +167,22 @@ def solve_circle_line_intersection( b = center_y c = y0 d = x0 - bb = -2 * a + 2 * c * k - 2 * d * (k ** 2) - 2 * b * k - aa = 1 + k ** 2 + bb = -2 * a + 2 * c * k - 2 * d * (k**2) - 2 * b * k + aa = 1 + k**2 rr = radius cc = ( - a ** 2 - + c ** 2 + a**2 + + c**2 - 2 * c * k * d - + (k ** 2) * (d ** 2) + + (k**2) * (d**2) - 2 * b * (c - k * d) - + b ** 2 - - rr ** 2 + + b**2 + - rr**2 ) - _x = (-bb + (bb ** 2 - 4 * aa * cc) ** 0.5) / 2 / aa + _x = (-bb + (bb**2 - 4 * aa * cc) ** 0.5) / 2 / aa _y = ground_surface_func(_x) # 验算结果 - equ = (center_x - _x) ** 2 + (center_y - _y) ** 2 - radius ** 2 + equ = (center_x - _x) ** 2 + (center_y - _y) ** 2 - radius**2 assert abs(equ) < 1e-5 return [_x, _y] @@ -189,12 +190,20 @@ def solve_circle_line_intersection( def min_i(string_len, u_ph): # 海拔修正 altitude = para.altitude - k_a = math.exp((altitude - 1000) / 8150) # 气隙海拔修正 + if altitude > 1000: + k_a = math.exp((altitude - 1000) / 8150) # 气隙海拔修正 + else: + k_a = 1 u_50 = 1 / k_a * (530 * string_len + 35) # 50045 上附录的公式，实际应该用负极性电压的公式 + # u_50 = 1 / k_a * (533 * string_len + 132) # 串放电路径 1000m海拔 + # u_50 = 1 / k_a * (477 * string_len + 99) # 串放电路径 2000m海拔 + # u_50 = 615 * string_len # 导线对塔身放电 1000m海拔 + # u_50= 263.32647401+533.90081562*string_len z_0 = 300 # 雷电波阻抗 z_c = 251 # 导线波阻抗 # 新版大手册公式 3-277 r = (u_50 + 2 * z_0 / (2 * z_0 + z_c) * u_ph) * (2 * z_0 + z_c) / (z_0 * z_c) + # r = 2 * (u_50 - u_ph) / z_c return r @@ -207,12 +216,15 @@ def thunder_density(i, td, ip_a, ip_b): # 雷电流幅值密度函数 r = -( -ip_b / ip_a / ((1 + (i / ip_a) ** ip_b) ** 2) * ((i / ip_a) ** (ip_b - 1)) ) + return r else: if td == 20: r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44) # 雷暴日20d + return r if td == 40: r = -(10 ** (-i / 88)) * math.log(10) * (-1 / 88) # 雷暴日40d - return r + return r + raise Exception("检查雷电参数!") def angel_density(angle): # 入射角密度函数 angle单位是弧度 @@ -221,12 +233,12 @@ def angel_density(angle): # 入射角密度函数 angle单位是弧度 def rs_fun(i): - r = 10 * (i ** 0.65) # 新版大手册公式3-271 + r = 10 * (i**0.65) # 新版大手册公式3-271 return r def rc_fun(i, u_ph): - r = 1.63 * ((5.015 * (i ** 0.578) - 0.001 * u_ph) ** 1.125) # 新版大手册公式3-272 + r = 1.63 * ((5.015 * (i**0.578) - 0.001 * u_ph * 1) ** 1.125) # 新版大手册公式3-272 return r @@ -235,9 +247,9 @@ def rg_fun(i_curt, h_cav, u_ph, typ="g"): rg = None if typ == "g": if h_cav < 40: - rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65) # 新版大手册公式3-273 + rg = (3.6 + 1.7 * math.log(43 - h_cav)) * (i_curt**0.65) # 新版大手册公式3-273 else: - rg = 5.5 * (i_curt ** 0.65) # 新版大手册公式3-273 + rg = 5.5 * (i_curt**0.65) # 新版大手册公式3-273 elif typ == "c": # 此时返回的是圆半径 rg = rc_fun(i_curt, u_ph) return rg @@ -294,7 +306,7 @@ def intersection_angle( # 点到直线的距离 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 @@ -404,14 +416,14 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0): 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) + radius**2 + ) * (k**2 + 1) d_fk = ( 2 * (k * center_x - center_y - k * line_x + line_y) * (center_x - line_x) - - 2 * (radius ** 2) * k + - 2 * (radius**2) * k ) if abs(d_fk) < 1e-5 and abs(line_x - center_x - radius) < 1e-5: # k不存在，角度为90°，k取一个很大的正数 @@ -447,7 +459,7 @@ def func_ng(td): # 地闪密度，通过雷暴日计算 if para.ng > 0: r = para.ng else: - r = 0.023 * (td ** 1.3) + r = 0.023 * (td**1.3) return r @@ -456,12 +468,12 @@ def circle_ground_surface_intersection(radius, center_x, center_y, ground_surfac # 最笨的办法，一个个去试 x_series = np.linspace(0, radius, int(radius / 0.001)) + center_x part_to_be_squared = ( - radius ** 2 - (x_series - center_x) ** 2 + 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 + 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] @@ -473,6 +485,8 @@ def circle_ground_surface_intersection(radius, center_x, center_y, ground_surfac # insulator_c_len绝缘子闪络距离 def arc_possibility(rated_voltage, insulator_c_len): # 建弧率 # 50064 中附录给的公式 - _e = rated_voltage / (3 ** 0.5) / insulator_c_len - r = (4.5 * (_e ** 0.75) - 14) * 1e-2 + # TODO 需要区分交直流 + # _e = rated_voltage / (3**0.5) / insulator_c_len #交流 + _e = abs(rated_voltage) / (1) / insulator_c_len # 直流 + r = (4.5 * (_e**0.75) - 14) * 1e-2 return r