From 94c4878b1ba21f3eb59180ba320ac1af813cfd3b Mon Sep 17 00:00:00 2001 From: facat Date: Sat, 11 Sep 2021 09:29:04 +0800 Subject: [PATCH] =?UTF-8?q?=E5=88=9D=E7=89=88=E5=AE=8C=E6=88=90=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- main.py | 199 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 199 insertions(+) create mode 100644 main.py diff --git a/main.py b/main.py new file mode 100644 index 0000000..9a0f578 --- /dev/null +++ b/main.py @@ -0,0 +1,199 @@ +import math +import ezdxf +import numpy as np + + +# 圆交点 +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 + 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 [] + # 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 + return [r, rg] + + +def draw(rs, rc, rg, h_gav, h_cav, dgc): + doc = ezdxf.new(dxfversion="R2010") + doc.layers.add("EGM", color=2) + msp = doc.modelspace() + msp.add_circle((0, h_gav), rs) + msp.add_line((0, 0), (0, h_gav)) # 地线 + msp.add_circle((dgc, h_cav), rc) + msp.add_line((dgc, 0), (dgc, h_cav)) # 导线 + msp.add_line((0, h_gav), (dgc, h_cav)) + msp.add_line((0, rg), (200, rg)) + # 计算圆交点 + 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) # 导线和圆的交点 + doc.saveas("egm.dxf") + solve_circle_intersection(rs, rc, h_gav, h_cav, dgc) + + +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) ** 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, 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) + return rg + + +def intersection_angel(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) + 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]) + return np.array([theta1, theta2]) + + +def bd_area(i_curt, u_ph, theta1, theta2): # 暴露弧的投影面积 + rc = rc_fun(i_curt, u_ph) + # 暂时不考虑雷电入射角的影响 + r = (math.cos(theta1) - math.cos(theta2)) * rc + return r + # r1=rc*(-math.cos(thyta2)+math.cos(thyta1)) + # 入射角密度函数积分 + # arrival_angle_fineness=0.0001 + # for calculus_arv_angle in np.linspace() + + +def egm(): + u_ph = 750 / 1.732 # 运行相电压 + h_cav = 60 # 导线对地平均高 + h_gav = h_cav + 9.5 + 7.2 + dgc = 2 + # 迭代法计算最大电流 + i_max = 0 + _min_i = 30 # 尝试的最小电流 + _max_i = 80 # 尝试的最大电流 + for i_bar in np.linspace(_min_i, _max_i, int((_max_i - _min_i) / 0.01)): # 雷电流 + 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(rc, rg, h_cav, dgc) + min_distance_intersection = ( + np.sum( + (np.array(circle_intersection) - np.array(circle_line_intersection)) + ** 2 + ) + ** 0.5 + ) # 计算两圆交点和地面直线交点的最小距离 + if min_distance_intersection < 0.01: + i_max = i_bar + draw(rs, rc, rg, h_gav, h_cav, dgc) + break + print(f"最大电流为{i_max:.2f}") + i_min = min_i(6.78, 750 / 1.732) + print(f"最小电流为{i_min:.2f}") + # 开始积分 + curt_fineness = 0.1 # 电流积分细度 + curt_segment_n = int((i_max - i_min) / curt_fineness) + d_curt = (i_max - i_min) / curt_segment_n + calculus = 0 + for curt in np.linspace(i_min, i_max, curt_segment_n): + cal_thyta_first = intersection_angel(dgc, h_gav, h_cav, curt, u_ph) + cal_bd_first = bd_area(curt, u_ph, cal_thyta_first[0], cal_thyta_first[1]) + cal_thyta_second = intersection_angel(dgc, h_gav, h_cav, curt + d_curt, u_ph) + cal_bd_second = bd_area( + curt + d_curt, u_ph, cal_thyta_second[0], cal_thyta_second[1] + ) + cal_thunder_density_first = thunder_density(curt) + cal_thunder_density_second = thunder_density(curt + d_curt) + calculus += ( + ( + cal_bd_first * cal_thunder_density_first + + cal_bd_second * cal_thunder_density_second + ) + / 2 + * d_curt + ) + n_sf=2*2.7/10*calculus + print(f'跳闸率是{n_sf}') + + # draw(rs, rc, rg, h_gav, h_cav, dgc) + + +if __name__ == "__main__": + thunder_density(2) + egm() + print("Finished.")