From 7f03fc2b9c653c1f7a05df0ef3205bb361d0998d Mon Sep 17 00:00:00 2001 From: n3040 Date: Wed, 22 Dec 2021 16:11:14 +0800 Subject: [PATCH] =?UTF-8?q?=E5=8F=82=E6=95=B0=E5=85=A8=E9=83=A8=E4=BB=8E?= =?UTF-8?q?=E5=A4=96=E9=83=A8=E8=AF=BB=E5=8F=96?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- core.py | 163 ++++++++++++++---- main.py | 468 ++++++++++++++++++++++++++++++--------------------- plot_data.py | 24 +-- server.py | 33 ++++ 4 files changed, 450 insertions(+), 238 deletions(-) create mode 100644 server.py diff --git a/core.py b/core.py index 55402b7..ccc41a6 100644 --- a/core.py +++ b/core.py @@ -7,6 +7,38 @@ gMSP = None gCount = 1 +class Parameter: + h_g_sag: float # 地线弧垂 + h_c_sag: float # 导线弧垂 + h_whole: float # 杆塔全高 + voltage_n: int # 工作电压分成多少份来计算 + td: int # 雷暴日 + insulator_c_len: float # 串子绝缘长度 + string_c_len: float + string_g_len: float + gc_x: [float] # 导、地线水平坐标 + ground_angels: [float] # 地面倾角，向下为正 + h_arm: float # 导、地线垂直坐标 + altitude: int # 海拔，单位米 + max_i: float # 最大尝试电流，单位kA + + +para = Parameter() + + +def rg_line_function_factory(_rg, ground_angel): # 返回一个地面捕雷线的直线方程 + y_d = _rg / math.cos(ground_angel) # y轴上的截距 + # 利用公式y-y0=k(x-x0) 得到直线公式 + y0 = y_d + x0 = 0 + k = math.tan(math.pi - ground_angel) + + def f(x): + return y0 + k * (x - x0) + + return f + + class Draw: def __init__(self): self._doc = ezdxf.new(dxfversion="R2010") @@ -14,7 +46,20 @@ class Draw: global gCAD gCAD = self - def draw(self, i_curt, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, color): + def draw( + self, + i_curt, + u_ph, + rs_x, + rs_y, + rc_x, + rc_y, + rg_x, + rg_y, + rg_type, + ground_angel, + color, + ): doc = self._doc msp = doc.modelspace() global gMSP @@ -33,8 +78,15 @@ class Draw: (rc_x, rc_y), circle_intersection, dxfattribs={"color": color} ) # 地线 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) + ground_angel_func = rg_line_function_factory(rg, ground_angel) + msp.add_line( + (0, ground_angel_func(0)), + (2000, ground_angel_func(2000)), + dxfattribs={"color": color}, + ) + circle_line_section = solve_circle_line_intersection( + rc, rc_x, rc_y, ground_angel_func + ) if not circle_line_section: pass else: @@ -94,40 +146,77 @@ def solve_circle_intersection( # 圆与捕雷线交点 -def solve_circle_line_intersection(radius, rg, center_x, center_y): - distance = distance_point_line(center_x, center_y, 0, rg, 0) # 捕雷线到暴露圆中点的距离 +def solve_circle_line_intersection( + radius, center_x, center_y, ground_surface_func +): # 返回交点的x和y坐标 + x0 = 0 + y0 = ground_surface_func(x0) + x1 = 1 + y1 = ground_surface_func(x1) + k = (y1 - y0) / (x1 - x0) + distance = distance_point_line(center_x, center_y, x0, y0, k) # 捕雷线到暴露圆中点的距离 if distance > radius: return [] else: - r = (radius ** 2 - (rg - center_y) ** 2) ** 0.5 + center_x - return [r, rg] + # r = (radius ** 2 - (rg - center_y) ** 2) ** 0.5 + center_x + a = center_x + b = center_y + c = y0 + d = x0 + bb = -2 * a + 2 * c * k - 2 * d * (k ** 2) - 2 * b * k + aa = 1 + k ** 2 + rr = radius + cc = ( + a ** 2 + + c ** 2 + - 2 * c * k * d + + (k ** 2) * (d ** 2) + - 2 * b * (c - k * d) + + b ** 2 + - rr ** 2 + ) + _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 + assert abs(equ) < 1e-5 + return [_x, _y] def min_i(string_len, u_ph): - u_50 = 530 * string_len + 35 + # 海拔修正 + altitude = para.altitude + k_a = math.exp(altitude / 8150) # 气隙海拔修正 + u_50 = 1 / k_a * (530 * string_len + 35) # 50045 上附录的公式，实际应该用负极性电压的公式 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) return r def thunder_density(i): # l雷电流幅值密度函数 - r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44) + td = para.td + r = None + if td == 20: + r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44) # 雷暴日20d + if td == 40: + r = -(10 ** (-i / 88)) * math.log(10) * (-1 / 88) # 雷暴日40d return r def angel_density(angle): # 入射角密度函数 angle单位是弧度 - r = 0.75 * abs((np.cos(angle - math.pi / 2) ** 3)) + r = 0.75 * abs((np.cos(angle - math.pi / 2) ** 3)) # 新版大手册公式3-275 return r def rs_fun(i): - r = 10 * (i ** 0.65) + 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) + r = 1.63 * ((5.015 * (i ** 0.578) - 0.001 * u_ph) ** 1.125) # 新版大手册公式3-272 return r @@ -136,16 +225,16 @@ 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) + rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65) # 新版大手册公式3-273 else: - rg = 5.5 * (i_curt ** 0.65) + rg = 5.5 * (i_curt ** 0.65) # 新版大手册公式3-273 elif typ == "c": # 此时返回的是圆半径 rg = rc_fun(i_curt, u_ph) return rg def intersection_angle( - rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_surface, rg_type + rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_angel, rg_type ): # 暴露弧的角度 rs = rs_fun(i_curt) rc = rc_fun(i_curt, u_ph) @@ -154,35 +243,38 @@ def intersection_angle( rs, rc, rs_x, rs_y, rc_x, rc_y ) # 两圆的交点 circle_line_or_rg_intersection = None + rg_line_func = rg_line_function_factory(rg, ground_angel) if rg_type == "g": circle_line_or_rg_intersection = solve_circle_line_intersection( - rc, rg, rc_x, rc_y + rc, rc_x, rc_y, rg_line_func ) # 暴露圆和补雷线的交点 if rg_type == "c": circle_line_or_rg_intersection = solve_circle_intersection( rg, rc, rg_x, rg_y, rc_x, rc_y ) # 两圆的交点 - 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 - ) + # TODO 应该是不存在落到地面线以下的情况，先把以下注释掉 + + # 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(rg, 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) if not circle_line_or_rg_intersection: - if rc_y - rc > rg: # rg在rc下面 - # 捕捉线太低了，对高塔无保护，θ_1从-90°开始计算。 + if rc_y - rc > rg_line_func(rc_x): # rg在rc下面 + # 捕捉线太低了，对高塔无保护，θ_1从-90°开始计算，即从与地面垂直的角度开始就已经暴露了。 theta1 = -math.pi / 2 else: theta1_line = np_circle_line_or_rg_intersection - np.array([rc_x, rc_y]) @@ -202,6 +294,7 @@ def func_calculus_pw(theta, max_w): if segments < 2: # 最大最小太小，没有可以积分的 return 0 w_samples, d_w = np.linspace(0, max_w, segments, retstep=True) + # 童中宇 750KV信洛Ⅰ线雷电防护性能研究 公式 3-10 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 @@ -216,7 +309,7 @@ def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y): # 对θ进行积分 # 求保护弧到直线的距离，判断是否相交 d_to_rs = distance_point_line(rs_x, rs_y, line_x, line_y, k) if d_to_rs < rs: # 相交 - # 要用过直线上一点到暴露弧的切线 + # 要用过这一点到保护弧的切线 new_k = tangent_line_k(line_x, line_y, rs_x, rs_y, rs, init_k=k) if new_k >= 0: max_w = math.atan(new_k) # 用于保护弧相切的角度 @@ -249,15 +342,16 @@ def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y): # 对θ进行积分 # ) # gCAD.save() pass + # 童中宇 750KV信洛Ⅰ线雷电防护性能研究 公式 3-10 r = rc / math.cos(theta) * func_calculus_pw(theta, max_w) return r def bd_area( - i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, ground_surface, rg_type + i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, ground_angel, rg_type ): # 暴露弧的投影面积 theta1, theta2 = intersection_angle( - rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_surface, rg_type + rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_angel, rg_type ) # θ角度 theta_fineness = 0.01 rc = rc_fun(i_curt, u_ph) @@ -364,6 +458,7 @@ def circle_ground_surface_intersection(radius, center_x, center_y, ground_surfac # u_ph是相电压 # 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 return r diff --git a/main.py b/main.py index 65c1c17..0d463f6 100644 --- a/main.py +++ b/main.py @@ -1,225 +1,309 @@ +import math import sys - +import tomli from loguru import logger from core import * import timeit def egm(): - h_g_avr_sag = 11.67 * 2 / 3 - h_c_avr_sag = 14.43 * 2 / 3 - h_whole = 130 # 杆塔全高 - voltage_n = 3 # 工作电压分成多少份来计算 - td = 20 # 雷暴日 - insulator_c_len = 6.6 # 串子绝缘长度 - string_c_len = 9.2 - string_g_len = 0.5 - gc_x = [17.9, 17, 15, 17.0] - - # 以后考虑地形角度,地面线 - def ground_surface(x): - return 0 - + if len(sys.argv) < 2: + toml_file_path = "default.toml" + # # logger.info('没指定计算文件!程序结束。') + # # sys.exit(0) + # h_g_sag = 20 # 地线弧垂 + # h_c_sag = 20 # 导线弧垂 + # h_whole = 106.1 # 杆塔全高 + # voltage_n = 3 # 工作电压分成多少份来计算 + # td = 20 # 雷暴日 + # insulator_c_len = 6.98 # 串子绝缘长度 + # string_c_len = 9.2 + # string_g_len = 0.63 + # gc_x = [32.2 / 2, 32.2 / 2, 15, 17.0] # 导、地线水平坐标 + # ground_angels = [40 / 180 * math.pi] # 地面倾角，向下为正 + # h_arm = 34 + # gc_y = [ + # h_whole - string_g_len - h_g_sag * 2 / 3, # 地线对地平均高 + # # h_whole - string_c_len - h_c_sag - 2.7, # 导线对地平均高 + # h_whole - string_c_len - h_c_sag * 2 / 3 - h_arm, # 导线对地平均高 + # # h_whole - string_c_len - h_c_sag - 35.7, # 导线对地平均高 + # ] + else: + toml_file_path = sys.argv[1] + logger.info(f"读取文件{toml_file_path}") + with open(toml_file_path, "rb") as toml_fs: + toml_dict = tomli.load(toml_fs) + toml_parameter = toml_dict["parameter"] + para.h_g_sag = toml_parameter["h_g_sag"] # 地线弧垂 + para.h_c_sag = toml_parameter["h_c_sag"] # 导线弧垂 + para.h_whole = toml_parameter["h_whole"] # 杆塔全高 + para.td = toml_parameter["td"] # 雷暴日 + para.insulator_c_len = toml_parameter["insulator_c_len"] # 串子绝缘长度 + para.string_c_len = toml_parameter["string_c_len"] + para.string_g_len = toml_parameter["string_g_len"] + para.gc_x = toml_parameter["gc_x"] # 导、地线水平坐标 + para.ground_angels = [ + angel / 180 * math.pi for angel in toml_parameter["ground_angels"] + ] # 地面倾角，向下为正 + para.h_arm = toml_parameter["h_arm"] + para.altitude = toml_parameter["altitude"] + para.max_i = toml_parameter["max_i"] + toml_optional = toml_dict["optional"] + para.voltage_n = toml_optional["voltage_n"] # 工作电压分成多少份来计算 + ######################################################### + # 以上是需要设置的参数 + h_whole = para.h_whole + string_g_len = para.string_g_len + string_c_len = para.string_c_len + h_g_sag = para.h_g_sag + h_c_sag = para.h_c_sag + gc_x = para.gc_x + h_arm = para.h_arm 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, # 导线对地平均高 + h_whole - string_g_len - h_g_sag * 2 / 3, # 地线对地平均高 ] + if len(h_arm) > 1: + for hoo in h_arm[1:]: + gc_y.append(hoo - string_c_len - h_c_sag * 2 / 3) 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() - # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13 + # 地闪密度 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13 + td = para.td ng = func_ng(td) - n_sf_phases = np.zeros((phase_n, voltage_n)) # 计算每一相的跳闸率 - if np.any(np.array(gc_y) < 0): - logger.info("导线可能掉地面了，程序退出。") - 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_foo < 2: - rg_type = "c" - rg_x = gc_x[phase_conductor_foo + 2] - rg_y = gc_y[phase_conductor_foo + 2] - else: + avr_n_sf = 0 # 考虑电压的影响计算的跳闸率 + ground_angels = para.ground_angels + for ground_angel in ground_angels: + logger.info(f"地面倾角{ground_angel/math.pi*180:.3f}°") + rg_type = None + rg_x = None + rg_y = None + cad = Draw() + voltage_n = para.voltage_n + n_sf_phases = np.zeros((phase_n, voltage_n)) # 存储每一相的跳闸率 + if np.any(np.array(gc_y) < 0): + logger.info("导线可能掉地面下了，程序退出。") + 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" - # TODO 保护角公式可能有问题，后面改 - shield_angle = ( - math.atan((rc_x - rs_x) / ((rs_y - rc_y) + string_c_len)) * 180 / math.pi - ) # 保护角 - logger.info(f"保护角{shield_angle:.3f}°") - logger.debug(f"最低相防护标识{rg_type}") - for u_bar in range(voltage_n): - u_ph = ( - -math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732 - ) # 运行相电压 - logger.info(f"计算第{phase_conductor_foo + 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) - ): # 雷电流 - # logger.info(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 - ) - i_max = i_bar - if not rg_rc_circle_intersection: # if circle_intersection is [] - logger.debug("保护弧和暴露弧无交点，检查设置参数。") - 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 - ) - 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: - # 暴露弧和捕捉弧无交点 + if phase_n > 1: # 多回路 + if phase_conductor_foo < 2: + rg_type = "c" # 捕捉弧有下面一相导线的击距代替 + 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 - rs_x) / ((rs_y - rc_y) + string_c_len)) + * 180 + / math.pi + ) # 保护角 + logger.info(f"保护角{shield_angle:.3f}°") + logger.debug(f"最低相防护标识{rg_type}") + for u_bar in range(voltage_n): # 计算不同工作电压下的跳闸率 + u_ph = ( + math.sqrt(2) + * 750 + * math.cos(2 * math.pi / voltage_n * u_bar) + / 1.732 + ) # 运行相电压 + logger.info(f"计算第{phase_conductor_foo + 1}相，电压为{u_ph:.2f}kV") + # 迭代法计算最大电流 + i_max = 0 + insulator_c_len = para.insulator_c_len + i_min = min_i(insulator_c_len, u_ph / 1.732) + _min_i = i_min # 尝试的最小电流 + _max_i = para.max_i # 尝试的最大电流 + # 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) + ): # 雷电流 + # logger.info(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) + rg_line_func = None if rg_type == "g": - if rg > rc_y: - i_min = i_bar - logger.info(f"捕捉弧在暴露弧之上，设置最小电流为{i_min:.2f}") + rg_line_func = rg_line_function_factory(rg, ground_angel) + ####### + # 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 + ) + i_max = i_bar + if not rg_rc_circle_intersection: # if circle_intersection is [] + logger.debug("保护弧和暴露弧无交点，检查设置参数。") + continue + circle_rc_line_or_rg_intersection = None + if rg_type == "g": + circle_rc_line_or_rg_intersection = ( + solve_circle_line_intersection(rc, rc_x, rc_y, rg_line_func) + ) + 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: + # 暴露弧和捕捉弧无交点 + if rg_type == "g": + if rg_line_func(rc_x) > rc_y: + i_min = i_bar # 用于后面判断最小和最大电流是否相等，相等意味着暴露弧一直被屏蔽 + logger.info(f"捕捉面在暴露弧之上，设置最小电流为{i_min:.2f}") + else: + logger.info("暴露弧和地面捕捉弧无交点，检查设置参数。") + continue else: - logger.info("暴露弧和捕捉弧无交点，检查设置参数。") - continue - else: - logger.info("暴露弧和捕捉弧无交点，检查设置参数。") - continue - min_distance_intersection = ( - np.sum( - ( - np.array(rg_rc_circle_intersection) - - np.array(circle_rc_line_or_rg_intersection) - ) - ** 2 - ) - ** 0.5 - ) # 计算两圆交点和地面直线交点的最小距离 - 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 - 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 = ( + logger.info("上面的导地线无法保护下面的导地线，检查设置参数。") + 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: - logger.info("暴露弧已经完全被屏蔽") - exposed_curve_shielded = True - break - # 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.save_as(f"egm{phase_conductor_foo + 1}.dxf") - # 判断是否导线已经被完全保护 - if abs(i_max - _max_i) < 1e-5: - logger.info("无法找到最大电流，可能是杆塔较高。") - logger.info(f"最大电流设置为自然界最大电流{i_max}kA") - logger.info(f"最大电流为{i_max:.2f}") - logger.info(f"最小电流为{i_min:.2f}") - if exposed_curve_shielded: - logger.info("暴露弧已经完全被屏蔽,不会跳闸。") - continue - curt_fineness = 0.1 # 电流积分细度 - if i_min > i_max or abs(i_min - i_max) < curt_fineness: - logger.info("最大电流小于最小电流，没有暴露弧。") - 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, + ) # 计算两圆交点和地面直线交点的最小距离 + 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 + if rg_type == "g": + circle_rs_line_or_rg_intersection = ( + solve_circle_line_intersection( + rs, rs_x, rs_y, rg_line_func + ) # 保护弧和捕雷弧的交点 + ) + 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: + logger.info("暴露弧已经完全被屏蔽") + exposed_curve_shielded = True + break + # if phase_conductor_foo == 2: + cad.draw( + i_min, u_ph, - rc_x, - rc_y, rs_x, rs_y, + rc_x, + rc_y, rg_x, rg_y, - ground_surface, rg_type, + ground_angel, + 2, ) - * 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 * arc_possibility(750, insulator_c_len) - avr_n_sf += n_sf / voltage_n - n_sf_phases[phase_conductor_foo][u_bar] = n_sf - logger.info(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}") - logger.info(f"跳闸率是{avr_n_sf:.6f}") - logger.info(f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=6)}") + cad.draw( + i_max, + u_ph, + rs_x, + rs_y, + rc_x, + rc_y, + rg_x, + rg_y, + rg_type, + ground_angel, + 6, + ) + cad.save_as(f"egm{phase_conductor_foo + 1}.dxf") + # 判断是否导线已经被完全保护 + if abs(i_max - _max_i) < 1e-5: + logger.info("无法找到最大电流，可能是杆塔较高。") + logger.info(f"最大电流设置为自然界最大电流{i_max}kA") + logger.info(f"最大电流为{i_max:.2f}") + logger.info(f"最小电流为{i_min:.2f}") + if exposed_curve_shielded: + logger.info("暴露弧已经完全被屏蔽,不会跳闸。") + continue + curt_fineness = 0.1 # 电流积分细度 + if i_min > i_max or abs(i_min - i_max) < curt_fineness: + logger.info("最大电流小于等于最小电流，没有暴露弧。") + 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_angel, + 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 * arc_possibility(750, insulator_c_len) + avr_n_sf += n_sf / voltage_n + n_sf_phases[phase_conductor_foo][u_bar] = n_sf + logger.info(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}") + logger.info(f"跳闸率是{avr_n_sf:.6f}") + logger.info( + f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=6)}" + ) def speed(): diff --git a/plot_data.py b/plot_data.py index 6a63be6..d687110 100644 --- a/plot_data.py +++ b/plot_data.py @@ -40,24 +40,24 @@ data_130m塔高_不同接地电阻 = np.array( ] ) -category_names_130m塔高_不同地线保护角 = ["-1", "-3", "-6"] +category_names_150m塔高_不同地线保护角 = ["-1", "-3", "-6"] -data_130m塔高_不同地线保护角 = np.array( +data_150m塔高_不同地线保护角 = np.array( [ - [0.000170, 0.103132], - [0.000168, 0.079659], - [0.000167, 0.055598], + [0.000440, 0.155414], + [0.000433, 0.128227], + [0.000429, 0.099819], ] ) -category_names_66m串长_不同塔高 = ['100', '120', '140'] +category_names_66m串长_不同塔高 = ["100", "120", "140"] data_66m串长_不同塔高 = np.array( [ - [0.000053 , 0.023285], - [0.000139 , 0.083229], - [0.000470 , 0.145586], + [0.000053, 0.023285], + [0.000139, 0.083229], + [0.000470, 0.145586], ] ) @@ -65,8 +65,8 @@ category_names_68m串长_不同塔高 = [100, 120, 140] data_68m串长_不同塔高 = np.array( [ - [0.000039 , 0.019094], - [0.000098 , 0.073033], - [0.000287 , 0.130923], + [0.000039, 0.019094], + [0.000098, 0.073033], + [0.000287, 0.130923], ] ) diff --git a/server.py b/server.py new file mode 100644 index 0000000..dfc18c9 --- /dev/null +++ b/server.py @@ -0,0 +1,33 @@ +from fastapi import FastAPI +import uvicorn +from pydantic import BaseModel + + +class CalculationParameter(BaseModel): + loop: str # single or double 回路数 + insulator_c_len: float # 绝缘长度 + string_c_len: float # 导线串长 + string_g_len: float # 地线串长 + td: int # 雷暴日 + h_g_avr_sag: float # 地线平均弧垂 + h_c_avr_sag: float # 导线平均弧垂 + h_whole: float # 杆塔全高 + gc_x: tuple[float] # 导、地线水平坐标 + ground_angels: tuple[float] # 地面倾角，向下为正 + + +fastapi_app = FastAPI() + + +@fastapi_app.post("/calculation") +async def calculation(cp: CalculationParameter): + print(cp) + return {"Hello": "World"} + + +def server_start(): + pass + + +if __name__ == "__main__": + uvicorn.run("server:fastapi_app", host="127.0.0.1", port=5000, log_level="info")