From 011db48f8bd9e81eea420e256d78980f9a5f7a0d Mon Sep 17 00:00:00 2001 From: facat Date: Tue, 21 Sep 2021 00:36:09 +0800 Subject: [PATCH] =?UTF-8?q?=E8=80=83=E8=99=91=E4=BA=86=E7=94=B5=E5=8E=8B?= =?UTF-8?q?=E7=9A=84=E5=BD=B1=E5=93=8D?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- core.py | 19 ++++++++++--------- main.py | 57 +++++++++++++++++++++++++++++++++++---------------------- 2 files changed, 45 insertions(+), 31 deletions(-) diff --git a/core.py b/core.py index 09fba3b..9996628 100644 --- a/core.py +++ b/core.py @@ -115,6 +115,7 @@ def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph): # 暴露弧的角度 rc, rg, dgc, h_cav ) # 暴露圆和补雷线的交点 np_circle_intersection = np.array(circle_intersection) + # TODO to be removed if not circle_intersection: abc = 123 theta2_line = np_circle_intersection - np.array([dgc, h_cav]) @@ -122,10 +123,6 @@ def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph): # 暴露弧的角度 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]) @@ -134,13 +131,14 @@ def distance_point_line(point_x, point_y, line_x, line_y, k): return d -def fun_calculus_pw(theta, max_w): +def func_calculus_pw(theta, max_w): w_fineness = 0.01 r_pw = 0 + # TODO to be removed 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) + w_samples, d_w = np.linspace(0, max_w, int(max_w / w_fineness), retstep=True) for cal_w in w_samples[:-1]: r_pw += ( ( @@ -164,7 +162,7 @@ def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): # 对θ进行积分 if d_to_rs < rs: # 相交 # 要用过直线上一点到暴露弧的切线 new_k = tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k) - + # TODO to be removed if not new_k: a = 12 tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k) @@ -172,9 +170,11 @@ def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): # 对θ进行积分 max_w = math.atan(new_k) # 用于保护弧相切的角度 elif new_k < 0: max_w = math.atan(new_k) + math.pi + # TODO to be removed if max_w < 0: abc = 123 tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k) + # TODO to be removed global gCount gCount += 1 if gCount % 100 == 0: @@ -189,6 +189,7 @@ def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): # 对θ进行积分 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) @@ -200,7 +201,7 @@ def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): # 对θ进行积分 # ) # gCAD.save() pass - r = rc / math.cos(theta) * fun_calculus_pw(theta, max_w) + r = rc / math.cos(theta) * func_calculus_pw(theta, max_w) return r @@ -212,7 +213,7 @@ def bd_area(i_curt, u_ph, dgc, h_gav, h_cav): # 暴露弧的投影面积 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 + theta1, theta2, int((theta2 - theta1) / theta_fineness), retstep=True ) for calculus_theta in theta_sample[:-1]: r_bd += ( diff --git a/main.py b/main.py index 4acba7e..b00ab9b 100644 --- a/main.py +++ b/main.py @@ -3,40 +3,49 @@ import timeit def egm(): - cccCount=0 - for u_bar in range(1): - u_ph = math.sqrt(1) * 750 * math.cos(2 * math.pi / 3 * 0) / 1.732 # 运行相电压 - h_whole = 140 # 杆塔全高 - insulator_c_len = 6.8 # 串子绝缘长度 - string_c_len = 9.2 - string_g_len = 0.5 - dgc = -2.9 # 导地线水平距离 - vertical_dgc = 2.7 # 导地线挂点垂直距离 - h_g_avr_sag = 11.67 * 2 / 3 - h_c_avr_sag = (14.43 - 11.67) * 2 / 3 - h_gav = h_whole - string_g_len - h_g_avr_sag # 地线对地平均高 - h_cav = h_gav - string_c_len - vertical_dgc - h_c_avr_sag # 导线对地平均高 - shield_angle = math.atan(dgc / (vertical_dgc + string_c_len)) * 180 / math.pi - print(f"保护角{shield_angle:.3f}°") + # TODO to be removed + cccCount = 0 + avr_n_sf = 0 # 考虑电压的影响 + voltage_n = 3 # 工作电压分成多少份来计算 + ng = func_ng(20) + h_whole = 140 # 杆塔全高 + insulator_c_len = 6.8 # 串子绝缘长度 + string_c_len = 9.2 + string_g_len = 0.5 + dgc = -0.9 # 导地线水平距离 + vertical_dgc = 2.7 # 导地线挂点垂直距离 + h_g_avr_sag = 11.67 * 2 / 3 + h_c_avr_sag = (14.43 - 11.67) * 2 / 3 + h_gav = h_whole - string_g_len - h_g_avr_sag # 地线对地平均高 + h_cav = h_gav - string_c_len - vertical_dgc - h_c_avr_sag # 导线对地平均高 + shield_angle = math.atan(dgc / (vertical_dgc + string_c_len)) * 180 / math.pi + print(f"保护角{shield_angle:.3f}°") + for u_bar in range(voltage_n): + u_ph = ( + math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732 + ) # 运行相电压 # 迭代法计算最大电流 i_max = 0 i_min = min_i(insulator_c_len, u_ph / 1.732) _min_i = i_min # 尝试的最小电流 _max_i = 200 # 尝试的最大电流 - #TODO remove it + # TODO remove it cad = Draw() cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2) for i_bar in np.linspace(_min_i, _max_i, int((_max_i - _min_i) / 0.1)): # 雷电流 - print(f"尝试计算电流为{i_bar:.2f}") + # print(f"尝试计算电流为{i_bar:.2f}") rs = rs_fun(i_bar) rc = rc_fun(i_bar, u_ph) rg = rg_fun(i_bar, h_cav) ####### - cccCount+=1 - if cccCount%30==0: + cccCount += 1 + if cccCount % 30 == 0: import core + core.gMSP.add_circle((0, h_gav), rs) - core.gMSP.add_circle((dgc, h_cav), rc_fun(i_bar, -u_ph),dxfattribs={"color": 4}) + core.gMSP.add_circle( + (dgc, h_cav), rc_fun(i_bar, -u_ph), dxfattribs={"color": 4} + ) core.gMSP.add_circle((dgc, h_cav), rc) ####### circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc) @@ -108,8 +117,12 @@ def egm(): / 2 * d_curt ) - n_sf = 2 * 2.7 / 10 * calculus # 跳闸率 - print(f"跳闸率是{n_sf:.6}") + n_sf = ( + 2 * ng / 10 * calculus + ) # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13 + avr_n_sf += n_sf / voltage_n + print(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}") + print(f"跳闸率是{avr_n_sf:.6}") def speed():