diff --git a/core.py b/core.py index 06a728c..edf8691 100644 --- a/core.py +++ b/core.py @@ -21,6 +21,9 @@ class Parameter: altitude: int # 海拔，单位米 max_i: float # 最大尝试电流，单位kA rated_voltage: float # 额定电压 + ng: float # 地闪密度 次/(每平方公里·每年) + Ip_a: float # 概率密度曲线系数a + Ip_b: float # 概率密度曲线系数b para = Parameter() @@ -186,7 +189,7 @@ def solve_circle_line_intersection( def min_i(string_len, u_ph): # 海拔修正 altitude = para.altitude - k_a = math.exp((altitude-1000) / 8150) # 气隙海拔修正 + k_a = math.exp((altitude - 1000) / 8150) # 气隙海拔修正 u_50 = 1 / k_a * (530 * string_len + 35) # 50045 上附录的公式，实际应该用负极性电压的公式 z_0 = 300 # 雷电波阻抗 z_c = 251 # 导线波阻抗 @@ -195,13 +198,20 @@ def min_i(string_len, u_ph): return r -def thunder_density(i): # l雷电流幅值密度函数 - td = para.td +def thunder_density(i, td, ip_a, ip_b): # 雷电流幅值密度函数 + # td = para.td r = None - if td == 20: - r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44) # 雷暴日20d - if td == 40: - r = -(10 ** (-i / 88)) * math.log(10) * (-1 / 88) # 雷暴日40d + # ip_a = para.Ip_a + # ip_b = para.Ip_b + if ip_a > 0 and ip_b > 0: + r = -( + -ip_b / ip_a / ((1 + (i / ip_a) ** ip_b) ** 2) * ((i / ip_a) ** (ip_b - 1)) + ) + else: + if td == 20: + r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44) # 雷暴日20d + if td == 40: + r = -(10 ** (-i / 88)) * math.log(10) * (-1 / 88) # 雷暴日40d return r @@ -433,7 +443,7 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0): return np.array(k_candidate)[np.max(k_angle) == k_angle].tolist()[-1] -def func_ng(td): # 地闪密度 +def func_ng(td): # 地闪密度，通过雷暴日计算 return 0.023 * (td ** 1.3) diff --git a/main.py b/main.py index d1c2f09..581f03a 100644 --- a/main.py +++ b/main.py @@ -8,30 +8,29 @@ import timeit # 打印参数 -def parameter_display(para: Parameter): - logger.info(f"额定电压 kV {para.rated_voltage}") - logger.info(f"导线弧垂 m {para.h_c_sag}") - logger.info(f"地线弧垂 m {para.h_g_sag}") - logger.info(f"全塔高 m {para.h_arm[0]}") - logger.info(f"串绝缘距离 m {para.insulator_c_len}") - logger.info(f"导线串长 m {para.string_c_len}") - logger.info(f"地线串长 m {para.string_g_len}") - logger.info(f"挂点垂直坐标 m {para.h_arm}") - logger.info(f"挂点水平坐标 m {para.gc_x}") - logger.info(f"地面倾角 ° {[an*180/math.pi for an in para.ground_angels]}") - logger.info(f"海拔高度 m {para.altitude}") - logger.info(f"雷暴日 d {para.td}") - - -def egm(): - if len(sys.argv) < 2: - toml_file_path = r"article.toml" +def parameter_display(para_dis: Parameter): + logger.info(f"额定电压 kV {para_dis.rated_voltage}") + logger.info(f"导线弧垂 m {para_dis.h_c_sag}") + logger.info(f"地线弧垂 m {para_dis.h_g_sag}") + logger.info(f"全塔高 m {para_dis.h_arm[0]}") + logger.info(f"串绝缘距离 m {para_dis.insulator_c_len}") + logger.info(f"导线串长 m {para_dis.string_c_len}") + logger.info(f"地线串长 m {para_dis.string_g_len}") + logger.info(f"挂点垂直坐标 m {para_dis.h_arm}") + logger.info(f"挂点水平坐标 m {para_dis.gc_x}") + logger.info(f"地面倾角 ° {[an * 180 / math.pi for an in para_dis.ground_angels]}") + logger.info(f"海拔高度 m {para_dis.altitude}") + if para_dis.ng > 0: + logger.info("不采用雷暴日计算地闪密度和雷电流密度") + logger.info(f"地闪密度 次/(每平方公里·每年) {para_dis.ng}") + logger.info(f"概率密度曲线系数a {para_dis.Ip_a}") + logger.info(f"概率密度曲线系数b {para_dis.Ip_b}") + pass else: - toml_file_path = sys.argv[1] - if not os.path.exists(toml_file_path): - logger.info(f'无法找到数据文件{toml_file_path}，程序退出。') - sys.exit(0) - logger.info(f"读取文件{toml_file_path}") + logger.info(f"雷暴日 d {para_dis.td}") + + +def read_parameter(toml_file_path): with open(toml_file_path, "rb") as toml_fs: toml_dict = tomli.load(toml_fs) toml_parameter = toml_dict["parameter"] @@ -49,9 +48,25 @@ def egm(): para.h_arm = toml_parameter["h_arm"] para.altitude = toml_parameter["altitude"] para.rated_voltage = toml_parameter["rated_voltage"] + toml_advance = toml_dict["advance"] + para.ng = toml_advance["ng"] # 地闪密度 + para.Ip_a = toml_advance["Ip_a"] # 概率密度曲线系数a + para.Ip_b = toml_advance["Ip_b"] # 概率密度曲线系数b toml_optional = toml_dict["optional"] para.voltage_n = toml_optional["voltage_n"] # 工作电压分成多少份来计算 para.max_i = toml_optional["max_i"] + + +def egm(): + if len(sys.argv) < 2: + toml_file_path = r"default.toml" + else: + toml_file_path = sys.argv[1] + if not os.path.exists(toml_file_path): + logger.info(f"无法找到数据文件{toml_file_path}，程序退出。") + sys.exit(0) + logger.info(f"读取文件{toml_file_path}") + read_parameter(toml_file_path) ######################################################### # 以上是需要设置的参数 parameter_display(para) @@ -223,7 +238,7 @@ def egm(): ** 0.5 ) if distance > rc: - logger.info("暴露弧已经完全被屏蔽") + logger.info(f"电流为{i_bar}kV时，暴露弧已经完全被屏蔽") exposed_curve_shielded = True break # if phase_conductor_foo == 2: @@ -273,6 +288,9 @@ def egm(): i_min, i_max, curt_segment_n + 1, retstep=True ) bd_area_vec = np.vectorize(bd_area) + td = para.td + ip_a = para.Ip_a + ip_b = para.Ip_b cal_bd_np = ( bd_area_vec( i_curt_samples, @@ -286,7 +304,7 @@ def egm(): ground_angel, rg_type, ) - * thunder_density(i_curt_samples) + * thunder_density(i_curt_samples, td, ip_a, ip_b) ) calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt # for i_curt in i_curt_samples[:-1]: diff --git a/metadata.yml b/metadata.yml index 0338e6a..6ea90f2 100644 --- a/metadata.yml +++ b/metadata.yml @@ -1,4 +1,4 @@ -Version: 1.1.0 +Version: 1.2.0 #CompanyName: My Imaginary Company #FileDescription: Simple App #InternalName: Simple App diff --git a/untest.py b/untest.py new file mode 100644 index 0000000..8e7d4f3 --- /dev/null +++ b/untest.py @@ -0,0 +1,29 @@ +import unittest +from core import thunder_density +from main import read_parameter +import numpy as np +import sympy + + +class Testing(unittest.TestCase): + @classmethod + def setUpClass(cls): + # read_parameter('default.toml') + pass + + def test_thunder_density(self): + i = sympy.symbols("i") + a = np.random.random() + b = np.random.random() + p = 1-1 / (1 + (i / a) ** b) + d_p = sympy.diff(p, i) + random_i = np.random.randint(10, 100) + v_from_thunder_density = thunder_density(random_i, 0, a, b) + v_from_diff = d_p.evalf(subs={i: random_i}) + self.assertTrue( + abs(v_from_thunder_density - v_from_diff) < 1e-5, "与自动微分结果不一致" + ) # add assertion here + + +if __name__ == "__main__": + unittest.main()