From ef60e1474b667097cd990ecfc7634efecda9540e Mon Sep 17 00:00:00 2001 From: facat Date: Mon, 13 Sep 2021 01:34:21 +0800 Subject: [PATCH] =?UTF-8?q?1.=E5=A4=84=E7=90=86=E4=BA=8690=C2=B0=EF=BC=8Ck?= =?UTF-8?q?=E4=B8=8D=E5=AD=98=E5=9C=A8=E7=9A=84=E6=83=85=E5=86=B5=E3=80=82?= =?UTF-8?q?=202.=E9=87=8D=E5=91=BD=E5=90=8D=E4=BA=86=E4=B8=80=E4=BA=9B?= =?UTF-8?q?=E5=87=BD=E6=95=B0=E3=80=82=203.=E5=A2=9E=E5=8A=A0=E4=BA=86?= =?UTF-8?q?=E4=B8=80=E4=BA=9B=E4=B8=8D=E7=BB=95=E5=87=BB=E7=9A=84=E5=88=A4?= =?UTF-8?q?=E6=96=AD=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- main.py | 137 ++++++++++++++++++++++++++++++++++---------------------- 1 file changed, 83 insertions(+), 54 deletions(-) diff --git a/main.py b/main.py index 6b293d3..e2ceab1 100644 --- a/main.py +++ b/main.py @@ -4,6 +4,7 @@ import numpy as np gCAD = None gMSP = None +gCount = 1 class Draw: @@ -59,9 +60,9 @@ def solve_circle_intersection(rs, rc, hgav, hcav, dgc): # 圆与地面线交点 -def solve_circle_line_intersection(rc, rg, h_cav, dgc): +def solve_circle_line_intersection(radius, rg, center_y, center_x): # TODO: 需要考虑地面捕雷线与暴露弧完全没交点的情况 - r = (rc ** 2 - (rg - h_cav) ** 2) ** 0.5 + dgc + r = (radius ** 2 - (rg - center_y) ** 2) ** 0.5 + center_x return [r, rg] @@ -93,11 +94,11 @@ def rc_fun(i, u_ph): return r -def rg_fun(i, h_cav): +def rg_fun(i_curt, h_cav): if h_cav < 40: rg = (3.6 + 1.7 ** math.log(43 - h_cav)) ** 0.65 else: - rg = 5.5 * (i ** 0.65) + rg = 5.5 * (i_curt ** 0.65) return rg @@ -148,17 +149,17 @@ def fun_calculus_pw(theta, max_w): return r_pw -def cal_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): +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 - max_w = theta + math.pi / 2 # 入射角 - k = math.tan(max_w) + 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) @@ -169,18 +170,21 @@ def cal_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): 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, k * (-500 - line_x) + line_y), (500, k * (500 - line_x) + line_y) - # ) - # gMSP.add_line((0, rg), (1000, rg)) - # gCAD.save() - tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k) - - pass + 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 @@ -195,11 +199,11 @@ def bd_area(i_curt, u_ph, dgc, h_gav, h_cav): # 暴露弧的投影面积 theta_sample, d_theta = np.linspace( theta1, theta2, int((theta2 - theta1) / theta_fineness) + 1, retstep=True ) - for cal_theta in theta_sample[:-1]: + for calculus_theta in theta_sample[:-1]: r_bd += ( ( - cal_bd(cal_theta, rc, rs, rg, dgc, h_cav, h_gav) - + cal_bd(cal_theta + d_theta, rc, rs, rg, dgc, h_cav, h_gav) + 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 @@ -217,31 +221,44 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0): # 牛顿法求解k # f(k)=(k*x1-y1-k*x0+y0)**2-R**2*(k**2+1) x1,y1是圆心 # TODO:需要检验k值不存在的情况 - k_candidate = [-100, 100] - 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 - ) - 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: - k_candidate[ind] = k + k_candidate = [-100, 100] + if abs(center_y - line_y) < 1e-4 and abs(line_x - center_x - radius) < 1e-4: + # 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-4: + dd = distance_point_line(center_x, center_y, line_x, line_y, k) + if abs(dd - radius) < 1e-5: + 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: @@ -253,9 +270,9 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0): def egm(): for u_bar in range(1): u_ph = math.sqrt(2) * 750 * math.cos(2 * math.pi / 6 * 0) / 1.732 # 运行相电压 - h_cav = 90 # 导线对地平均高 - h_gav = h_cav + 9.5 + 2.7 - dgc = 2.9 # 导地线水平距离 + h_gav = 30 + h_cav = h_gav - 9.5 - 2.7 - 5 # 导线对地平均高 + dgc = -2.9 # 导地线水平距离 # 迭代法计算最大电流 i_max = 0 _min_i = 20 # 尝试的最小电流 @@ -263,23 +280,20 @@ def egm(): for i_bar in np.linspace(_min_i, _max_i, int((_max_i - _min_i) / 0.1)): # 雷电流 print(f"尝试计算电流为{i_bar:.2f}") rs = rs_fun(i_bar) - # if not np.isreal(rs): - # continue rc = rc_fun(i_bar, u_ph) - # if not np.isreal(rc): - # continue rg = rg_fun(i_bar, h_cav) - # if not np.isreal(rg): - # continue circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc) if not circle_intersection: # if circle_intersection is [] continue - circle_line_intersection = solve_circle_line_intersection( + circle_rc_line_intersection = solve_circle_line_intersection( rc, rg, h_cav, dgc ) min_distance_intersection = ( np.sum( - (np.array(circle_intersection) - np.array(circle_line_intersection)) + ( + np.array(circle_intersection) + - np.array(circle_rc_line_intersection) + ) ** 2 ) ** 0.5 @@ -287,14 +301,29 @@ def egm(): i_max = i_bar if min_distance_intersection < 0.1: break + if circle_intersection[1] < circle_rc_line_intersection[1]: + circle_rs_line_intersection = solve_circle_line_intersection( + rs, rg, h_gav, 0 + ) + # 判断与保护弧的交点是否在暴露弧外面 + distance = ( + np.sum( + (np.array(circle_rs_line_intersection) - np.array([dgc, h_cav])) + ** 2 + ) + ** 0.5 + ) + if distance > rc: + print("暴露弧已经完全被屏蔽") + break i_min = min_i(6.78, 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) cad.save() + # 判断是否导线已经被完全保护 if abs(i_max - _max_i) < 1e-5: print("无法找到最大电流，可能是杆塔较高。") - i_max = 300 print(f"最大电流设置为自然界最大电流{i_max}kA") print(f"最大电流为{i_max:.2f}") print(f"最小电流为{i_min:.2f}") @@ -328,6 +357,6 @@ def egm(): if __name__ == "__main__": - thunder_density(2) + tangent_line_k(1, 0, 0, 0, 1) egm() print("Finished.")