From dd44de030edf790ae41af69276e0374548e4ce8d Mon Sep 17 00:00:00 2001 From: facat Date: Tue, 21 Sep 2021 20:00:03 +0800 Subject: [PATCH] =?UTF-8?q?=E5=88=9D=E6=AD=A5=E5=AE=8C=E6=88=90=E4=BA=86?= =?UTF-8?q?=E5=8F=8C=E5=9B=9E=E8=B7=AF=E5=85=AC=E5=BC=8F?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- core.py | 140 ++++++++++++++++++++++++++++++++++++-------------------- main.py | 75 +++++++++++++++++++----------- 2 files changed, 140 insertions(+), 75 deletions(-) diff --git a/core.py b/core.py index 8766001..7631be5 100644 --- a/core.py +++ b/core.py @@ -15,26 +15,41 @@ class Draw: global gCAD gCAD = self - def draw(self, i_curt, u_ph, h_gav, h_cav, dgc, color): + def draw(self, i_curt, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 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}) + rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type) + msp.add_circle((rs_x, rs_y), rs, dxfattribs={"color": color}) + msp.add_line((0, 0), (rs_x, rs_y)) # 地线 + msp.add_circle((rc_x, rc_y), rc, dxfattribs={"color": color+2}) + msp.add_line((rc_x, 0), (rc_x, rc_y)) # 导线 + msp.add_line((rs_x, rs_y), (rc_x, rc_y)) + # 角度线 + circle_intersection = solve_circle_intersection(rs, rc, rs_x, rs_y, rc_x, rc_y) + msp.add_line( + (rc_x, rc_y), circle_intersection, dxfattribs={"color": color} + ) # 地线 + if rg_type == "g": + msp.add_line((0, rg), (2000, rg), dxfattribs={"color": color}) + circle_line_section = solve_circle_line_intersection(rc, rg, rc_x, rc_y) + msp.add_line( + (rc_x, rc_y), circle_line_section, dxfattribs={"color": color} + ) # 导线和圆的交点 + if rg_type == "c": + msp.add_circle((rg_x, rg_y), rg, dxfattribs={"color": color}) + rg_rc_intersection = solve_circle_intersection( + rg, rc, rg_x, rg_y, rc_x, rc_y + ) + msp.add_line( + (rc_x, rc_y), rg_rc_intersection, 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 @@ -42,15 +57,26 @@ class Draw: # 圆交点 -def solve_circle_intersection(rs, rc, h_gav, h_cav, dgc): +def solve_circle_intersection( + radius1, + radius2, + center_x1, + center_y1, + center_x2, + center_y2, +): # 用牛顿法求解 - x = rc # 初始值 - y = rc # 初始值 + x = radius2 # 初始值 + y = radius2 # 初始值 + # TODO 考虑出现2个解的情况 for bar in range(0, 10): - A = [[-2 * x, -2 * (y - h_gav)], [-2 * (x - dgc), -2 * (y - h_cav)]] + A = [ + [-2 * (x - center_x1), -2 * (y - center_y1)], + [-2 * (x - center_x2), -2 * (y - center_y2)], + ] b = [ - x ** 2 + (y - h_gav) ** 2 - rs ** 2, - (x - dgc) ** 2 + (y - h_cav) ** 2 - rc ** 2, + (x - center_x1) ** 2 + (y - center_y1) ** 2 - radius1 ** 2, + (x - center_x2) ** 2 + (y - center_y2) ** 2 - radius2 ** 2, ] X_set = np.linalg.solve(A, b) x += X_set[0] @@ -84,7 +110,7 @@ def thunder_density(i): # l雷电流幅值密度函数 def angel_density(angle): # 入射角密度函数 angle单位是弧度 - r = 0.75 * (np.cos(angle - math.pi / 2) ** 3) + r = 0.75 * abs((np.cos(angle - math.pi / 2) ** 3)) return r @@ -95,35 +121,50 @@ def rs_fun(i): 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) +# typ 如果是g，代表捕雷线公式，c代表暴露弧公式 +def rg_fun(i_curt, h_cav, u_ph, typ="g"): + rg = None + if typ == "g": + 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) + elif typ == "c": # 此时返回的是圆半径 + rg = rc_fun(i_curt, u_ph) return rg -def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph): # 暴露弧的角度 +def intersection_angle( + rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, rg_type +): # 暴露弧的角度 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 - ) # 暴露圆和补雷线的交点 + rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type) + circle_intersection = solve_circle_intersection( + rs, rc, rs_x, rs_y, rc_x, rc_y + ) # 两圆的交点 + circle_line_or_rg_intersection = None + if rg_type == "g": + circle_line_or_rg_intersection = solve_circle_line_intersection( + rc, rg, rc_x, rc_y + ) # 暴露圆和补雷线的交点 + if rg_type == "c": + circle_line_or_rg_intersection = solve_circle_intersection( + rg, rc, rg_x, rg_y, rc_x, rc_y + ) # 两圆的交点 np_circle_intersection = np.array(circle_intersection) - theta2_line = np_circle_intersection - np.array([dgc, h_cav]) + theta2_line = np_circle_intersection - np.array([rc_x, rc_y]) 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]) + np_circle_line_or_rg_intersection = np.array(circle_line_or_rg_intersection) + theta1_line = np_circle_line_or_rg_intersection - np.array([rc_x, rc_y]) theta1 = math.atan(theta1_line[1] / theta1_line[0]) return np.array([theta1, theta2]) +# 点到直线的距离 def distance_point_line(point_x, point_y, line_x, line_y, k) -> float: d = abs(k * point_x - point_y - k * line_x + line_y) / ((k ** 2 + 1) ** 0.5) return d @@ -137,17 +178,17 @@ def func_calculus_pw(theta, max_w): return r_pw -def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): # 对θ进行积分 +def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y): # 对θ进行积分 max_w = 0 # 求暴露弧上一点的切线 - line_x = math.cos(theta) * rc + dgc - line_y = math.sin(theta) * rc + h_cav + line_x = math.cos(theta) * rc + rc_x + line_y = math.sin(theta) * rc + rc_y k = math.tan(theta + math.pi / 2) # 入射角 # 求保护弧到直线的距离，判断是否相交 - d_to_rs = distance_point_line(0, h_gav, line_x, line_y, k) + d_to_rs = distance_point_line(rs_x, rs_y, 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) + new_k = tangent_line_k(line_x, line_y, rs_x, rs_y, rs, init_k=k) if new_k >= 0: max_w = math.atan(new_k) # 用于保护弧相切的角度 elif new_k < 0: @@ -183,20 +224,21 @@ def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav): # 对θ进行积分 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) # θ角度 +def bd_area(i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, rg_type): # 暴露弧的投影面积 + theta1, theta2 = intersection_angle( + rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, rg_type + ) # θ角度 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 + rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type) theta_sample, d_theta = np.linspace( theta1, theta2, int((theta2 - theta1) / theta_fineness), retstep=True ) if len(theta_sample) < 2: return 0 vec_calculus_bd = np.vectorize(calculus_bd) - calculus_bd_np = vec_calculus_bd(theta_sample, rc, rs, rg, dgc, h_cav, h_gav) + calculus_bd_np = vec_calculus_bd(theta_sample, rc, rs, rg, rc_x, rc_y, rs_x, rs_y) r_bd = np.sum(calculus_bd_np[:-1] + calculus_bd_np[1:]) / 2 * d_theta # for calculus_theta in theta_sample[:-1]: # r_bd += ( @@ -214,7 +256,7 @@ 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不存在 @@ -248,10 +290,10 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0): # 把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 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: diff --git a/main.py b/main.py index 88862fc..5d958c5 100644 --- a/main.py +++ b/main.py @@ -6,20 +6,26 @@ import timeit def egm(): avr_n_sf = 0 # 考虑电压的影响计算的跳闸率 - voltage_n = 3 # 工作电压分成多少份来计算 + voltage_n = 1 # 工作电压分成多少份来计算 ng = func_ng(20) h_whole = 140 # 杆塔全高 insulator_c_len = 6.8 # 串子绝缘长度 string_c_len = 9.2 string_g_len = 0.5 - dgc = -0.0 # 导地线水平距离 + rc_x = -0.0 # 导地线水平距离 + rs_x = 0 + rg_x = 0 vertical_dgc = 2.7 # 导地线挂点垂直距离 h_g_avr_sag = 11.67 * 2 / 3 h_c_avr_sag = 14.43 * 2 / 3 - h_gav = h_whole - string_g_len - h_g_avr_sag # 地线对地平均高 - h_cav = h_whole - string_c_len - vertical_dgc - h_c_avr_sag # 导线对地平均高 - shield_angle = math.atan(dgc / (vertical_dgc + string_c_len)) * 180 / math.pi + rs_y = h_whole - string_g_len - h_g_avr_sag # 地线对地平均高 + rc_y = h_whole - string_c_len - vertical_dgc - h_c_avr_sag # 导线对地平均高 + rg_y = rc_y - 20 + shield_angle = ( + math.atan(rc_x / (vertical_dgc + string_c_len)) * 180 / math.pi + ) # 保护角 print(f"保护角{shield_angle:.3f}°") + rg_type = "c" for u_bar in range(voltage_n): u_ph = ( math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732 @@ -30,12 +36,12 @@ def egm(): _min_i = i_min # 尝试的最小电流 _max_i = 200 # 尝试的最大电流 cad = Draw() - cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2) + # cad.draw(i_min, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 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) + rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type) ####### # cccCount += 1 # if cccCount % 30 == 0: @@ -47,20 +53,28 @@ def egm(): # ) # core.gMSP.add_circle((dgc, h_cav), rc) ####### - 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, dgc, h_cav + rg_rc_circle_intersection = solve_circle_intersection( + rs, rc, rs_x, rs_y, rc_x, rc_y ) - if not circle_rc_line_intersection: + if not rg_rc_circle_intersection: # if circle_intersection is [] + print("保护弧和暴露弧无交点，检查设置参数。程序退出。") + continue + circle_rc_line_or_rg_intersection = None + if rg_type == "g": + circle_rc_line_or_rg_intersection = solve_circle_line_intersection( + rc, rg, rc_x, rc_y + ) + elif rg_type == "c": + circle_rc_line_or_rg_intersection = solve_circle_intersection( + rg, rc, rg_x, rg_y, rc_x, rc_y + ) + if not circle_rc_line_or_rg_intersection: continue min_distance_intersection = ( np.sum( ( - np.array(circle_intersection) - - np.array(circle_rc_line_intersection) + np.array(rg_rc_circle_intersection) + - np.array(circle_rc_line_or_rg_intersection) ) ** 2 ) @@ -69,14 +83,24 @@ 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, 0, h_gav - ) + # 判断是否以完全被保护 + if rg_rc_circle_intersection[1] < circle_rc_line_or_rg_intersection[1]: + circle_rs_line_or_rg_intersection = None + if rg_type == "g": + circle_rs_line_or_rg_intersection = solve_circle_line_intersection( + rs, rg, rs_x, rs_y + ) + if rg_type == "c": + circle_rs_line_or_rg_intersection = solve_circle_intersection( + rs, rg, rs_x, rs_y, rg_x, rg_y + ) # 判断与保护弧的交点是否在暴露弧外面 distance = ( np.sum( - (np.array(circle_rs_line_intersection) - np.array([dgc, h_cav])) + ( + np.array(circle_rs_line_or_rg_intersection) + - np.array([rc_x, rc_y]) + ) ** 2 ) ** 0.5 @@ -84,8 +108,8 @@ def egm(): if distance > rc: print("暴露弧已经完全被屏蔽") break - 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.draw(i_min, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 2) + cad.draw(i_max, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 6) cad.save() # 判断是否导线已经被完全保护 if abs(i_max - _max_i) < 1e-5: @@ -99,13 +123,12 @@ def egm(): return # 开始积分 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 ) bd_area_vec = np.vectorize(bd_area) cal_bd_np = bd_area_vec( - i_curt_samples, u_ph, dgc, h_gav, h_cav + i_curt_samples, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, rg_type ) * thunder_density(i_curt_samples) calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt # for i_curt in i_curt_samples[:-1]: