From db2788f116d69361046b6849e89e3d0ed088152e Mon Sep 17 00:00:00 2001 From: facat Date: Sun, 12 Sep 2021 16:55:11 +0800 Subject: [PATCH] =?UTF-8?q?=E5=8A=A0=E5=85=A5=E4=BA=86=E8=AE=A1=E7=AE=97?= =?UTF-8?q?=E6=9C=80=E5=A4=A7=E5=85=A5=E5=B0=84=E8=A7=92=E5=BA=A6=E7=9A=84?= =?UTF-8?q?=E5=85=AC=E5=BC=8F?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- main.py | 130 +++++++++++++++++++++++++++++++++++++++----------------- 1 file changed, 90 insertions(+), 40 deletions(-) diff --git a/main.py b/main.py index d563fba..023fb70 100644 --- a/main.py +++ b/main.py @@ -2,30 +2,37 @@ import math import ezdxf import numpy as np +gCAD = None +gMSP = None + 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): + 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) + 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) + 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), (200, rg)) + msp.add_line((0, rg), (200, 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) # 导线和圆的交点 + # 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 @@ -34,12 +41,6 @@ class Draw: # 圆交点 def solve_circle_intersection(rs, rc, hgav, hcav, dgc): - # x = Symbol('x', real=True) - # y = Symbol('y', real=True) - # equ = [ - # x ** 2 + (y - hgav) ** 2 - rs ** 2, - # (x - dgc) ** 2 + (y - hcav) ** 2 - rc ** 2, - # ] # 用牛顿法求解 x = 300 y = 300 @@ -55,22 +56,12 @@ def solve_circle_intersection(rs, rc, hgav, hcav, dgc): if np.all(np.abs(X_set) < 1e-5): return [x, y] return [] - # list_set = list(X_set) - # solve_set = nonlinsolve(equ, [x, y]) - # print(ask(Q.real(solve_set))) - # list_set = list(solve_set) - # pprint(list_set) - # if not np.all(np.isreal(list_set)): - # return [] - # for value in list_set: - # if value[0] > 0 and value[1] > 1: - # return value - # return [] # 圆与地面线交点 -def solve_circle_line_intersection(rc, rg, hcav, dgc): - r = (rc ** 2 - (rg - hcav) ** 2) ** 0.5 + dgc +def solve_circle_line_intersection(rc, rg, h_cav, dgc): + # TODO: 需要考虑地面捕雷线与暴露弧完全没交点的情况 + r = (rc ** 2 - (rg - h_cav) ** 2) ** 0.5 + dgc return [r, rg] @@ -110,27 +101,62 @@ def rg_fun(i, h_cav): return rg -def intersection_angel(dgc, h_gav, h_cav, i_curt, u_ph): # 暴露弧的角度 +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) + 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) 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]) - if theta1 < 0: - # print(f"θ_1角度为负数{theta1:.4f},人为设置为0") - theta1 = 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 bd_area(i_curt, u_ph, dgc, h_gav, h_cav): # 暴露弧的投影面积 - theta1, theta2 = intersection_angel(dgc, h_gav, h_cav, i_curt, u_ph) + theta1, theta2 = intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph) rc = rc_fun(i_curt, u_ph) + rs = rs_fun(i_curt) + rg = rg_fun(i_curt, h_cav) + # 求暴露弧上一点的切线 + line_x = math.cos(theta1) * rc + dgc + line_y = math.sin(theta1) * rc + h_cav + max_w = 0 # 入射角 + if theta1 < 0: + max_w = theta1 + math.pi / 2 + k = math.tan(max_w) + # 求保护弧到直线的距离，判断是否相交 + 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) + max_w = math.atan(new_k) # 用于保护弧相切的角度 + 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) + ) + gMSP.add_line((0, rg), (1000, rg)) + gCAD.save() + pass + + # k=tangent_line_k(point_x, point_y, dgc, h_cav,rc) # 暂时不考虑雷电入射角的影响 r = (math.cos(theta1) - math.cos(theta2)) * rc return r @@ -140,11 +166,35 @@ def bd_area(i_curt, u_ph, dgc, h_gav, h_cav): # 暴露弧的投影面积 # 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 = init_k + 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 + ) + d_k = -fk / d_fk + k += d_k + if abs(d_k) < 1e-5: + dd = distance_point_line(center_x, center_y, line_x, line_y, k) + if abs(dd - radius) < 1e-5: + return k + return None + + def egm(): u_ph = 750 / 1.732 # 运行相电压 h_cav = 160 # 导线对地平均高 - h_gav = h_cav + 9.5 + 2.2 - dgc = 2 # 导地线水平距离 + h_gav = h_cav + 9.5 + 2.7 + dgc = -2 # 导地线水平距离 # 迭代法计算最大电流 i_max = 0 _min_i = 20 # 尝试的最小电流 @@ -176,8 +226,8 @@ def egm(): break i_min = min_i(6.78, 750 / 1.732) cad = Draw() - cad.draw(i_min, u_ph, h_gav, h_cav, dgc) - cad.draw(i_max, u_ph, h_gav, h_cav, dgc) + cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2) + cad.draw(i_max, u_ph, h_gav, h_cav, dgc, 6) cad.save() if abs(i_max - _max_i) < 1e-5: print("无法找到最大电流，可能是杆塔较高。") @@ -189,7 +239,7 @@ def egm(): print("最大电流小于最小电流，没有暴露弧，程序结束。") return # 开始积分 - curt_fineness = 0.001 # 电流积分细度 + curt_fineness = 0.1 # 电流积分细度 curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份 calculus = 0 i_curt_samples, d_curt = np.linspace(i_min, i_max, curt_segment_n + 1, retstep=True)