diff --git a/core.py b/core.py new file mode 100644 index 0000000..7d4b807 --- /dev/null +++ b/core.py @@ -0,0 +1,276 @@ +import math +import ezdxf +import numpy as np + +gCAD = None +gMSP = None +gCount = 1 + + +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, h_gav, h_cav, dgc, 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, h_cav) + msp.add_circle((0, h_gav), rs, dxfattribs={"color": color}) + msp.add_line((0, 0), (0, h_gav)) # 地线 + msp.add_circle((dgc, h_cav), rc, dxfattribs={"color": color}) + msp.add_line((dgc, 0), (dgc, h_cav)) # 导线 + msp.add_line((0, h_gav), (dgc, h_cav)) + msp.add_line((0, rg), (2000, rg), dxfattribs={"color": color}) + # 计算圆交点 + # circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc) + # msp.add_line((0, h_gav), circle_intersection) # 地线 + # msp.add_line((dgc, h_cav), circle_intersection) # 导线 + # circle_line_section = solve_circle_line_intersection(rc, rg, h_cav, dgc) + # msp.add_line((0, 0), circle_line_section) # 导线和圆的交点 + + def save(self): + doc = self._doc + doc.saveas("egm.dxf") + + +# 圆交点 +def solve_circle_intersection(rs, rc, hgav, hcav, dgc): + # 用牛顿法求解 + x = rc # 初始值 + y = rc # 初始值 + for bar in range(0, 10): + A = [[-2 * x, -2 * (y - hgav)], [-2 * (x - dgc), -2 * (y - hcav)]] + b = [ + x ** 2 + (y - hgav) ** 2 - rs ** 2, + (x - dgc) ** 2 + (y - hcav) ** 2 - rc ** 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, rg, center_x, center_y): + distance = distance_point_line(center_x, center_y, 0, rg, 0) # 捕雷线到暴露圆中点的距离 + if distance > radius: + return [] + else: + r = (radius ** 2 - (rg - center_y) ** 2) ** 0.5 + center_x + return [r, rg] + + +def min_i(string_len, u_ph): + u_50 = 530 * string_len + 35 + z_0 = 300 # 雷电波阻抗 + z_c = 251 # 导线波阻抗 + 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) + return r + + +def angel_density(angle): # 入射角密度函数 angle单位是弧度 + r = 0.75 * (math.cos(angle - math.pi / 2) ** 3) + return r + + +def rs_fun(i): + r = 10 * (i ** 0.65) + return r + + +def rc_fun(i, u_ph): + r = 1.63 * ((5.015 * (i ** 0.578) - 0.001 * u_ph) ** 1.125) + # r=14.7*(i**0.42) + return r + + +def rg_fun(i_curt, h_cav): + if h_cav < 40: + rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65) + else: + rg = 5.5 * (i_curt ** 0.65) + return rg + + +def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph): # 暴露弧的角度 + rs = rs_fun(i_curt) + rc = rc_fun(i_curt, u_ph) + rg = rg_fun(i_curt, h_cav) + circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc) # 两圆的交点 + circle_line_intersection = solve_circle_line_intersection( + rc, rg, dgc, h_cav + ) # 暴露圆和补雷线的交点 + np_circle_intersection = np.array(circle_intersection) + if not circle_intersection: + abc = 123 + theta2_line = np_circle_intersection - np.array([dgc, h_cav]) + theta2 = math.atan(theta2_line[1] / theta2_line[0]) + np_circle_line_intersection = np.array(circle_line_intersection) + theta1_line = np_circle_line_intersection - np.array([dgc, h_cav]) + theta1 = math.atan(theta1_line[1] / theta1_line[0]) + # 考虑雷电入射角度,所以theta1可以小于0,即计算从侧面击中的雷 + # if theta1 < 0: + # # print(f"θ_1角度为负数{theta1:.4f},人为设置为0") + # theta1 = 0 + return np.array([theta1, theta2]) + + +def distance_point_line(point_x, point_y, line_x, line_y, k): + d = abs(k * point_x - point_y - k * line_x + line_y) / ((k ** 2 + 1) ** 0.5) + return d + + +def fun_calculus_pw(theta, max_w): + w_fineness = 0.01 + r_pw = 0 + if int(max_w / w_fineness) < 0: + abc = 123 + pass + w_samples, d_w = np.linspace(0, max_w, int(max_w / w_fineness) + 1, retstep=True) + for cal_w in w_samples[: -1]: + r_pw += ( + ( + abs(angel_density(cal_w)) * math.sin(theta - cal_w + math.pi) + + abs(angel_density(cal_w + d_w)) + * math.sin(theta - cal_w - d_w + math.pi) + ) + / 2 + ) * d_w + return r_pw + + +def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): # 对θ进行积分 + max_w=0 + # 求暴露弧上一点的切线 + line_x = math.cos(theta) * rc + dgc + line_y = math.sin(theta) * rc + h_cav + k = math.tan(theta + math.pi / 2) # 入射角 + # 求保护弧到直线的距离,判断是否相交 + d_to_rs = distance_point_line(0, h_gav, line_x, line_y, k) + if d_to_rs < rs: # 相交 + # 要用过直线上一点到暴露弧的切线 + new_k = tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k) + + if not new_k: + a = 12 + tangent_line_k(line_x, line_y, 0, h_gav, 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 + if max_w < 0: + abc = 123 + tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k) + global 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 # 入射角 + 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 + r = rc / math.cos(theta) * fun_calculus_pw(theta, max_w) + return r + + +def bd_area(i_curt, u_ph, dgc, h_gav, h_cav): # 暴露弧的投影面积 + theta1, theta2 = intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph) # θ角度 + theta_fineness = 0.01 + rc = rc_fun(i_curt, u_ph) + rs = rs_fun(i_curt) + rg = rg_fun(i_curt, h_cav) + r_bd = 0 + theta_sample, d_theta = np.linspace( + theta1, theta2, int((theta2 - theta1) / theta_fineness) + 1, retstep=True + ) + 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 + for bar in range(0, 30): + 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 + # 把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] diff --git a/main.py b/main.py index 2e73b91..f4bed65 100644 --- a/main.py +++ b/main.py @@ -1,279 +1,15 @@ -import math -import ezdxf -import numpy as np - -gCAD = None -gMSP = None -gCount = 1 - - -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, h_gav, h_cav, dgc, 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, h_cav) - msp.add_circle((0, h_gav), rs, dxfattribs={"color": color}) - msp.add_line((0, 0), (0, h_gav)) # 地线 - msp.add_circle((dgc, h_cav), rc, dxfattribs={"color": color}) - msp.add_line((dgc, 0), (dgc, h_cav)) # 导线 - msp.add_line((0, h_gav), (dgc, h_cav)) - msp.add_line((0, rg), (2000, rg), dxfattribs={"color": color}) - # 计算圆交点 - # circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc) - # msp.add_line((0, h_gav), circle_intersection) # 地线 - # msp.add_line((dgc, h_cav), circle_intersection) # 导线 - # circle_line_section = solve_circle_line_intersection(rc, rg, h_cav, dgc) - # msp.add_line((0, 0), circle_line_section) # 导线和圆的交点 - - def save(self): - doc = self._doc - doc.saveas("egm.dxf") - - -# 圆交点 -def solve_circle_intersection(rs, rc, hgav, hcav, dgc): - # 用牛顿法求解 - x = rc # 初始值 - y = rc # 初始值 - for bar in range(0, 10): - A = [[-2 * x, -2 * (y - hgav)], [-2 * (x - dgc), -2 * (y - hcav)]] - b = [ - x ** 2 + (y - hgav) ** 2 - rs ** 2, - (x - dgc) ** 2 + (y - hcav) ** 2 - rc ** 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, rg, center_y, center_x): - # TODO: 需要考虑地面捕雷线与暴露弧完全没交点的情况 - r = (radius ** 2 - (rg - center_y) ** 2) ** 0.5 + center_x - return [r, rg] - - -def min_i(string_len, u_ph): - u_50 = 530 * string_len + 35 - z_0 = 300 # 雷电波阻抗 - z_c = 251 # 导线波阻抗 - 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) - return r - - -def angel_density(angle): # 入射角密度函数 angle单位是弧度 - r = 0.75 * (math.cos(angle - math.pi / 2) ** 3) - return r - - -def rs_fun(i): - r = 10 * (i ** 0.65) - return r - - -def rc_fun(i, u_ph): - r = 1.63 * ((5.015 * (i ** 0.578) - 0.001 * u_ph) ** 1.125) - return r - - -def rg_fun(i_curt, h_cav): - if h_cav < 40: - rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65) - else: - rg = 5.5 * (i_curt ** 0.65) - return rg - - -def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph): # 暴露弧的角度 - rs = rs_fun(i_curt) - rc = rc_fun(i_curt, u_ph) - rg = rg_fun(i_curt, h_cav) - circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc) # 两圆的交点 - circle_line_intersection = solve_circle_line_intersection( - rc, rg, h_cav, dgc - ) # 暴露圆和补雷线的交点 - np_circle_intersection = np.array(circle_intersection) - if not circle_intersection: - abc = 123 - theta2_line = np_circle_intersection - np.array([dgc, h_cav]) - theta2 = math.atan(theta2_line[1] / theta2_line[0]) - np_circle_line_intersection = np.array(circle_line_intersection) - theta1_line = np_circle_line_intersection - np.array([dgc, h_cav]) - theta1 = math.atan(theta1_line[1] / theta1_line[0]) - # 考虑雷电入射角度,所以theta1可以小于0,即计算从侧面击中的雷 - # if theta1 < 0: - # # print(f"θ_1角度为负数{theta1:.4f},人为设置为0") - # theta1 = 0 - return np.array([theta1, theta2]) - - -def distance_point_line(point_x, point_y, line_x, line_y, k): - d = abs(k * point_x - point_y - k * line_x + line_y) / ((k ** 2 + 1) ** 0.5) - return d - - -def fun_calculus_pw(theta, max_w): - w_fineness = 0.01 - r_pw = 0 - if int(max_w / w_fineness) < 0: - abc = 123 - pass - w_samples, d_w = np.linspace(0, max_w, int(max_w / w_fineness) + 1, retstep=True) - for cal_w in w_samples: - r_pw += ( - ( - abs(angel_density(cal_w)) * math.sin(theta - cal_w + math.pi) - + abs(angel_density(cal_w + d_w)) - * math.sin(theta - cal_w + math.pi - d_w) - ) - / 2 - ) * d_w - return r_pw - - -def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): # 对θ进行积分 - # 求暴露弧上一点的切线 - line_x = math.cos(theta) * rc + dgc - line_y = math.sin(theta) * rc + h_cav - k = math.tan(theta + math.pi / 2) # 入射角 - # 求保护弧到直线的距离,判断是否相交 - d_to_rs = distance_point_line(0, h_gav, line_x, line_y, k) - if d_to_rs < rs: # 相交 - # 要用过直线上一点到暴露弧的切线 - new_k = tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k) - - if not new_k: - a = 12 - tangent_line_k(line_x, line_y, 0, h_gav, 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 - if max_w < 0: - abc = 123 - tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k) - global gCount - gCount += 1 - if gCount % 1000 == 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, new_k * (-500 - line_x) + line_y), - (500, new_k * (500 - line_x) + line_y), - ) - gMSP.add_line((0, rg), (1000, rg)) - gCAD.save() - else: - max_w = theta + math.pi / 2 # 入射角 - r = rc / math.cos(theta) * fun_calculus_pw(theta, max_w) - return r - - -def bd_area(i_curt, u_ph, dgc, h_gav, h_cav): # 暴露弧的投影面积 - theta1, theta2 = intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph) - theta_fineness = 0.01 - rc = rc_fun(i_curt, u_ph) - rs = rs_fun(i_curt) - rg = rg_fun(i_curt, h_cav) - r_bd = 0 - theta_sample, d_theta = np.linspace( - theta1, theta2, int((theta2 - theta1) / theta_fineness) + 1, retstep=True - ) - 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 - - # r1=rc*(-math.cos(thyta2)+math.cos(thyta1)) - # 入射角密度函数积分 - # arrival_angle_fineness=0.0001 - # for calculus_arv_angle in np.linspace() - - -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是圆心 - # TODO:需要检验k值不存在的情况 - - 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 - for bar in range(0, 30): - 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 - # 把k转化成相应的角度,从x开始,逆时针为正 - k_angle = [] - for kk in k_candidate: - if kk == 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] +from core import * +import timeit def egm(): for u_bar in range(1): - u_ph = math.sqrt(2) * 750 * math.cos(2 * math.pi / 6 * 0) / 1.732 # 运行相电压 - h_whole = 150 # 杆塔全高 - h_gav = h_whole - 0.5 - 11.67 * 2 / 3 - h_cav = h_gav - 9.5 - 2.7 - (14.43 - 11.67) * 2 / 3 # 导线对地平均高 - dgc = 12 # 导地线水平距离 + u_ph = math.sqrt(1) * 750 * math.cos(2 * math.pi / 3 * 0) / 1.732 # 运行相电压 + h_whole = 140 # 杆塔全高 + string_len = 6.8 # 串子绝缘长度 + h_gav = h_whole - 0.5 - 11.67 * 2 / 3 # 地线对地平均高 + h_cav = h_gav - 9.2 - 2.7 - (14.43 - 11.67) * 2 / 3 # 导线对地平均高 + dgc = 0 # 导地线水平距离 # 迭代法计算最大电流 i_max = 0 _min_i = 20 # 尝试的最小电流 @@ -285,10 +21,13 @@ def egm(): rg = rg_fun(i_bar, h_cav) circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc) if not circle_intersection: # if circle_intersection is [] + # print("保护弧和暴露弧无交点,检查设置参数。程序退出。") continue circle_rc_line_intersection = solve_circle_line_intersection( - rc, rg, h_cav, dgc + rc, rg, dgc, h_cav ) + if not circle_rc_line_intersection: + continue min_distance_intersection = ( np.sum( ( @@ -304,7 +43,7 @@ def egm(): break if circle_intersection[1] < circle_rc_line_intersection[1]: circle_rs_line_intersection = solve_circle_line_intersection( - rs, rg, h_gav, 0 + rs, rg, 0, h_gav ) # 判断与保护弧的交点是否在暴露弧外面 distance = ( @@ -317,7 +56,7 @@ def egm(): if distance > rc: print("暴露弧已经完全被屏蔽") break - i_min = min_i(6.78, u_ph / 1.732) + i_min = min_i(string_len, u_ph / 1.732) cad = Draw() cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2) cad.draw(i_max, u_ph, h_gav, h_cav, dgc, 6) @@ -338,7 +77,7 @@ def egm(): i_curt_samples, d_curt = np.linspace( i_min, i_max, curt_segment_n + 1, retstep=True ) - for i_curt in i_curt_samples: + 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) @@ -351,13 +90,17 @@ def egm(): / 2 * d_curt ) - n_sf = 2 * 2.7 / 10 * calculus # 调整率 + n_sf = 2 * 2.7 / 10 * calculus # 跳闸率 print(f"跳闸率是{n_sf:.6}") - # draw(rs, rc, rg, h_gav, h_cav, dgc) + +def speed(): + a = 0 + for bar in range(100000000): + a += bar if __name__ == "__main__": - tangent_line_k(1, 0, 0, 0, 1) - egm() + run_time = timeit.timeit("egm()", globals=globals(), number=1) + print(f"运行时间:{run_time:.2f}s") print("Finished.")