import math import ezdxf import numpy as np gCAD = None gMSP = None gCount = 1 class Parameter: h_g_sag: float # 地线弧垂 h_c_sag: 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 rated_voltage: float # 额定电压 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") self._doc.layers.add("EGM", color=2) 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, ground_angel, color, ): doc = self._doc msp = doc.modelspace() global gMSP gMSP = msp rs = rs_fun(i_curt) rc = rc_fun(i_curt, u_ph) 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_line((rc_x, 0), (rc_x, rc_y)) # 导线 msp.add_line((rs_x, rs_y), (rc_x, rc_y)) # 角度线 circle_intersection = solve_circle_intersection(rs, rc, rs_x, rs_y, rc_x, rc_y) msp.add_line( (rc_x, rc_y), circle_intersection, dxfattribs={"color": color} ) # 地线 if rg_type == "g": 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: 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( rg, rc, rg_x, rg_y, rc_x, rc_y ) 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) # 导线 def save(self): doc = self._doc doc.saveas("egm.dxf") def save_as(self, file_name): doc = self._doc doc.saveas(file_name) # 圆交点 def solve_circle_intersection( radius1, radius2, center_x1, center_y1, center_x2, center_y2, ): # 用牛顿法求解 x = radius2 + center_x2 # 初始值 y = radius2 + center_y2 # 初始值 # TODO 考虑出现2个解的情况 for bar 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_set = np.linalg.solve(A, b) x += X_set[0] y += X_set[1] if np.all(np.abs(X_set) < 1e-5): return [x, y] return [] # 圆与捕雷线交点 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 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): # 海拔修正 altitude = para.altitude k_a = math.exp((altitude-1000) / 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雷电流幅值密度函数 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)) # 新版大手册公式3-275 return r def rs_fun(i): 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 return r # typ 如果是g，代表捕雷线公式，c代表暴露弧公式 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 else: 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_angel, rg_type ): # 暴露弧的角度 rs = rs_fun(i_curt) rc = rc_fun(i_curt, u_ph) rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type) circle_intersection = solve_circle_intersection( 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, 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 ) # 两圆的交点 # 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_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]) theta1 = math.atan(theta1_line[1] / theta1_line[0]) return np.array([theta1, theta2]) # 点到直线的距离 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) return d def func_calculus_pw(theta, max_w): w_fineness = 0.01 segments = int(max_w / w_fineness) 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 def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y): # 对θ进行积分 max_w = 0 # 求暴露弧上一点的切线 line_x = math.cos(theta) * rc + rc_x line_y = math.sin(theta) * rc + rc_y k = math.tan(theta + math.pi / 2) # 入射角 # 求保护弧到直线的距离，判断是否相交 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) # 用于保护弧相切的角度 elif new_k < 0: max_w = math.atan(new_k) + math.pi # TODO to be removed # global gCount # gCount = gCount+1 # if gCount % 100 == 0: # gMSP.add_circle((0, h_gav), rs) # gMSP.add_circle((dgc, h_cav), rc) # gMSP.add_line((dgc, h_cav), (line_x, line_y)) # gMSP.add_line( # (-500, new_k * (-500 - line_x) + line_y), # (500, new_k * (500 - line_x) + line_y), # ) # gCAD.save() # pass else: max_w = theta + math.pi / 2 # 入射角 # TODO to be removed if gCount % 200 == 0: # # intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph) # gMSP.add_circle((0, h_gav), rs) # gMSP.add_circle((dgc, h_cav), rc) # gMSP.add_line((dgc, h_cav), (line_x, line_y)) # gMSP.add_line( # (-500, k * (-500 - line_x) + line_y), # (500, k * (500 - line_x) + line_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_angel, rg_type ): # 暴露弧的投影面积 theta1, theta2 = intersection_angle( 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) rs = rs_fun(i_curt) rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type) 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) calculus_bd_np = vec_calculus_bd(theta_sample, rc, rs, rg, rc_x, rc_y, rs_x, rs_y) r_bd = np.sum(calculus_bd_np[:-1] + calculus_bd_np[1:]) / 2 * d_theta # for calculus_theta in theta_sample[:-1]: # r_bd += ( # ( # calculus_bd(calculus_theta, rc, rs, rg, dgc, h_cav, h_gav) # + calculus_bd(calculus_theta + d_theta, rc, rs, rg, dgc, h_cav, h_gav) # ) # / 2 # * d_theta # ) return r_bd def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0): # 直线方程为 y-y0=k(x-x0)，x0和y0为经过直线的任意一点 # 牛顿法求解k # f(k)=(k*x1-y1-k*x0+y0)**2-R**2*(k**2+1) x1,y1是圆心 # 已考虑两个解的判别 k_candidate = [-100, 100] if abs(center_y - line_y) < 1 and abs(line_x - center_x - radius) < 1: # k不存在 k_candidate = [99999999, 99999999] else: for ind, k_cdi in enumerate(list(k_candidate)): k = k_candidate[ind] k_candidate[ind] = None 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) d_fk = ( 2 * (k * center_x - center_y - k * line_x + line_y) * (center_x - line_x) - 2 * (radius ** 2) * k ) if abs(d_fk) < 1e-5 and abs(line_x - center_x - radius) < 1e-5: # k不存在，角度为90°，k取一个很大的正数 k_candidate[ind] = 99999999999999 break d_k = -fk / d_fk k += d_k if abs(d_k) < 1e-3: dd = distance_point_line(center_x, center_y, line_x, line_y, k) 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 >= 0: k_angle.append(math.atan(kk)) if kk < 0: k_angle.append(math.pi + math.atan(kk)) # 返回相对x轴最大的角度k return np.array(k_candidate)[np.max(k_angle) == k_angle].tolist()[-1] 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 # 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