facat
/
egm
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

加了一些TODO

This commit is contained in:
n3040 2024-01-03 14:28:28 +08:00
parent 9b852235f1
commit 5dc1613d2a
1 changed files with 43 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

72
core.py
View File

@ -1,6 +1,7 @@
import math import math
import ezdxf import ezdxf
import numpy as np import numpy as np
from typing import List
gCAD = None gCAD = None
gMSP = None gMSP = None
@ -15,8 +16,8 @@ class Parameter:
insulator_c_len: float # 串子绝缘长度 insulator_c_len: float # 串子绝缘长度
string_c_len: float string_c_len: float
string_g_len: float string_g_len: float
gc_x: [float] # 导、地线水平坐标 gc_x: List[float] # 导、地线水平坐标
ground_angels: [float] # 地面倾角，向下为正 ground_angels: List[float] # 地面倾角，向下为正
h_arm: float # 导、地线垂直坐标 h_arm: float # 导、地线垂直坐标
altitude: int # 海拔，单位米 altitude: int # 海拔，单位米
max_i: float # 最大尝试电流，单位kA max_i: float # 最大尝试电流，单位kA
@ -131,14 +132,14 @@ def solve_circle_intersection(
x = radius2 + center_x2 # 初始值 x = radius2 + center_x2 # 初始值
y = radius2 + center_y2 # 初始值 y = radius2 + center_y2 # 初始值
# TODO 考虑出现2个解的情况 # TODO 考虑出现2个解的情况
for bar in range(0, 10): for _ in range(0, 10):
A = [ A = [
[-2 * (x - center_x1), -2 * (y - center_y1)], [-2 * (x - center_x1), -2 * (y - center_y1)],
[-2 * (x - center_x2), -2 * (y - center_y2)], [-2 * (x - center_x2), -2 * (y - center_y2)],
] ]
b = [ b = [
(x - center_x1) ** 2 + (y - center_y1) ** 2 - radius1 ** 2, (x - center_x1) ** 2 + (y - center_y1) ** 2 - radius1**2,
(x - center_x2) ** 2 + (y - center_y2) ** 2 - radius2 ** 2, (x - center_x2) ** 2 + (y - center_y2) ** 2 - radius2**2,
] ]
X_set = np.linalg.solve(A, b) X_set = np.linalg.solve(A, b)
x += X_set[0] x += X_set[0]
@ -166,22 +167,22 @@ def solve_circle_line_intersection(
b = center_y b = center_y
c = y0 c = y0
d = x0 d = x0
bb = -2 * a + 2 * c * k - 2 * d * (k ** 2) - 2 * b * k bb = -2 * a + 2 * c * k - 2 * d * (k**2) - 2 * b * k
aa = 1 + k ** 2 aa = 1 + k**2
rr = radius rr = radius
cc = ( cc = (
a ** 2 a**2
+ c ** 2 + c**2
- 2 * c * k * d - 2 * c * k * d
+ (k ** 2) * (d ** 2) + (k**2) * (d**2)
- 2 * b * (c - k * d) - 2 * b * (c - k * d)
+ b ** 2 + b**2
- rr ** 2 - rr**2
) )
_x = (-bb + (bb ** 2 - 4 * aa * cc) ** 0.5) / 2 / aa _x = (-bb + (bb**2 - 4 * aa * cc) ** 0.5) / 2 / aa
_y = ground_surface_func(_x) _y = ground_surface_func(_x)
# 验算结果 # 验算结果
equ = (center_x - _x) ** 2 + (center_y - _y) ** 2 - radius ** 2 equ = (center_x - _x) ** 2 + (center_y - _y) ** 2 - radius**2
assert abs(equ) < 1e-5 assert abs(equ) < 1e-5
return [_x, _y] return [_x, _y]
@ -189,12 +190,20 @@ def solve_circle_line_intersection(
def min_i(string_len, u_ph): def min_i(string_len, u_ph):
# 海拔修正 # 海拔修正
altitude = para.altitude altitude = para.altitude
k_a = math.exp((altitude - 1000) / 8150) # 气隙海拔修正 if altitude > 1000:
k_a = math.exp((altitude - 1000) / 8150) # 气隙海拔修正
else:
k_a = 1
u_50 = 1 / k_a * (530 * string_len + 35) # 50045 上附录的公式，实际应该用负极性电压的公式 u_50 = 1 / k_a * (530 * string_len + 35) # 50045 上附录的公式，实际应该用负极性电压的公式
# u_50 = 1 / k_a * (533 * string_len + 132) # 串放电路径 1000m海拔
# u_50 = 1 / k_a * (477 * string_len + 99) # 串放电路径 2000m海拔
# u_50 = 615 * string_len # 导线对塔身放电 1000m海拔
# u_50= 263.32647401+533.90081562*string_len
z_0 = 300 # 雷电波阻抗 z_0 = 300 # 雷电波阻抗
z_c = 251 # 导线波阻抗 z_c = 251 # 导线波阻抗
# 新版大手册公式 3-277 # 新版大手册公式 3-277
r = (u_50 + 2 * z_0 / (2 * z_0 + z_c) * u_ph) * (2 * z_0 + z_c) / (z_0 * z_c) r = (u_50 + 2 * z_0 / (2 * z_0 + z_c) * u_ph) * (2 * z_0 + z_c) / (z_0 * z_c)
# r = 2 * (u_50 - u_ph) / z_c
return r return r
@ -207,12 +216,15 @@ def thunder_density(i, td, ip_a, ip_b): # 雷电流幅值密度函数
r = -( r = -(
-ip_b / ip_a / ((1 + (i / ip_a) ** ip_b) ** 2) * ((i / ip_a) ** (ip_b - 1)) -ip_b / ip_a / ((1 + (i / ip_a) ** ip_b) ** 2) * ((i / ip_a) ** (ip_b - 1))
) )
return r
else: else:
if td == 20: if td == 20:
r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44) # 雷暴日20d r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44) # 雷暴日20d
return r
if td == 40: if td == 40:
r = -(10 ** (-i / 88)) * math.log(10) * (-1 / 88) # 雷暴日40d r = -(10 ** (-i / 88)) * math.log(10) * (-1 / 88) # 雷暴日40d
return r return r
raise Exception("检查雷电参数!")
def angel_density(angle): # 入射角密度函数 angle单位是弧度 def angel_density(angle): # 入射角密度函数 angle单位是弧度
@ -221,12 +233,12 @@ def angel_density(angle): # 入射角密度函数 angle单位是弧度
def rs_fun(i): def rs_fun(i):
r = 10 * (i ** 0.65) # 新版大手册公式3-271 r = 10 * (i**0.65) # 新版大手册公式3-271
return r return r
def rc_fun(i, u_ph): def rc_fun(i, u_ph):
r = 1.63 * ((5.015 * (i ** 0.578) - 0.001 * u_ph) ** 1.125) # 新版大手册公式3-272 r = 1.63 * ((5.015 * (i**0.578) - 0.001 * u_ph * 1) ** 1.125) # 新版大手册公式3-272
return r return r
@ -235,9 +247,9 @@ def rg_fun(i_curt, h_cav, u_ph, typ="g"):
rg = None rg = None
if typ == "g": if typ == "g":
if h_cav < 40: if h_cav < 40:
rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65) # 新版大手册公式3-273 rg = (3.6 + 1.7 * math.log(43 - h_cav)) * (i_curt**0.65) # 新版大手册公式3-273
else: else:
rg = 5.5 * (i_curt ** 0.65) # 新版大手册公式3-273 rg = 5.5 * (i_curt**0.65) # 新版大手册公式3-273
elif typ == "c": # 此时返回的是圆半径 elif typ == "c": # 此时返回的是圆半径
rg = rc_fun(i_curt, u_ph) rg = rc_fun(i_curt, u_ph)
return rg return rg
@ -294,7 +306,7 @@ def intersection_angle(
# 点到直线的距离 # 点到直线的距离
def distance_point_line(point_x, point_y, line_x, line_y, k) -> float: 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) d = abs(k * point_x - point_y - k * line_x + line_y) / ((k**2 + 1) ** 0.5)
return d return d
@ -404,14 +416,14 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
max_iteration = 30 max_iteration = 30
for bar in range(0, max_iteration): for bar in range(0, max_iteration):
fk = (k * center_x - center_y - k * line_x + line_y) ** 2 - ( fk = (k * center_x - center_y - k * line_x + line_y) ** 2 - (
radius ** 2 radius**2
) * (k ** 2 + 1) ) * (k**2 + 1)
d_fk = ( d_fk = (
2 2
* (k * center_x - center_y - k * line_x + line_y) * (k * center_x - center_y - k * line_x + line_y)
* (center_x - line_x) * (center_x - line_x)
- 2 * (radius ** 2) * k - 2 * (radius**2) * k
) )
if abs(d_fk) < 1e-5 and abs(line_x - center_x - radius) < 1e-5: if abs(d_fk) < 1e-5 and abs(line_x - center_x - radius) < 1e-5:
# k不存在，角度为90°，k取一个很大的正数 # k不存在，角度为90°，k取一个很大的正数
@ -447,7 +459,7 @@ def func_ng(td): # 地闪密度，通过雷暴日计算
if para.ng > 0: if para.ng > 0:
r = para.ng r = para.ng
else: else:
r = 0.023 * (td ** 1.3) r = 0.023 * (td**1.3)
return r return r
@ -456,12 +468,12 @@ def circle_ground_surface_intersection(radius, center_x, center_y, ground_surfac
# 最笨的办法，一个个去试 # 最笨的办法，一个个去试
x_series = np.linspace(0, radius, int(radius / 0.001)) + center_x x_series = np.linspace(0, radius, int(radius / 0.001)) + center_x
part_to_be_squared = ( part_to_be_squared = (
radius ** 2 - (x_series - center_x) ** 2 radius**2 - (x_series - center_x) ** 2
) # 有可能出现-0.00001的数值，只是一个数值稳定问题。 ) # 有可能出现-0.00001的数值，只是一个数值稳定问题。
part_to_be_squared[ part_to_be_squared[
(part_to_be_squared < 0) & (abs(part_to_be_squared) < 1e-3) (part_to_be_squared < 0) & (abs(part_to_be_squared) < 1e-3)
] = 0 # 强制为0 ] = 0 # 强制为0
y_series = center_y - part_to_be_squared ** 0.5 y_series = center_y - part_to_be_squared**0.5
ground_surface_y = ground_surface(x_series) ground_surface_y = ground_surface(x_series)
equal_location = np.abs(ground_surface_y - y_series) < 0.5 equal_location = np.abs(ground_surface_y - y_series) < 0.5
r_x = x_series[equal_location][0] r_x = x_series[equal_location][0]
@ -473,6 +485,8 @@ def circle_ground_surface_intersection(radius, center_x, center_y, ground_surfac
# insulator_c_len绝缘子闪络距离 # insulator_c_len绝缘子闪络距离
def arc_possibility(rated_voltage, insulator_c_len): # 建弧率 def arc_possibility(rated_voltage, insulator_c_len): # 建弧率
# 50064 中附录给的公式 # 50064 中附录给的公式
_e = rated_voltage / (3 ** 0.5) / insulator_c_len # TODO 需要区分交直流
r = (4.5 * (_e ** 0.75) - 14) * 1e-2 # _e = rated_voltage / (3**0.5) / insulator_c_len #交流
_e = abs(rated_voltage) / (1) / insulator_c_len # 直流
r = (4.5 * (_e**0.75) - 14) * 1e-2
return r return r