From a153e69eb76de39c9f95b1b7849fd6fdf1dfa8c5 Mon Sep 17 00:00:00 2001 From: dmy Date: Mon, 2 Mar 2026 18:18:46 +0800 Subject: [PATCH] =?UTF-8?q?=E6=8F=90=E4=BA=A4=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Makefile | 2 +- README.md | 160 ++++++++++++++++++++++ animation.py | 4 +- article.toml | 6 +- main-batch.py | 365 ++++++++++++++++++++++++++++++++++++++++++++++++++ main.py | 4 +- 6 files changed, 533 insertions(+), 8 deletions(-) create mode 100644 README.md create mode 100644 main-batch.py diff --git a/Makefile b/Makefile index bdf9a18..d0d0757 100644 --- a/Makefile +++ b/Makefile @@ -11,4 +11,4 @@ build: .PHONY:clean clean: - rm -fr build \ No newline at end of file + rm -fr build diff --git a/README.md b/README.md new file mode 100644 index 0000000..da739e8 --- /dev/null +++ b/README.md @@ -0,0 +1,160 @@ +# EGM - 输电线路绕击跳闸率计算程序 + +基于电气几何模型（Electro-Geometric Model, EGM）的架空输电线路雷电防护性能计算工具，用于评估输电线路的绕击跳闸率。 + +## 功能特点 + +- 支持单回和双回输电线路的绕击跳闸率计算 +- 考虑工作电压对雷电击距的影响 +- 支持地面倾角参数 +- 支持自定义地闪密度和雷电流概率密度曲线 +- 支持交流/直流线路计算 +- 输出CAD图形（DXF格式）可视化击距模型 +- 提供动画演示模式（可选项） + +## 安装 + +### 环境要求 + +- Python >= 3.12 + +### 安装依赖 + +```bash +# 使用 uv（推荐） +uv sync + +# 或使用 pip +pip install -r requirements.txt +``` + +### 依赖包 + +- ezdxf - DXF文件生成 +- loguru - 日志记录 +- matplotlib - 数据可视化和动画 +- numpy - 数值计算 +- tomli - TOML配置文件解析 + +## 使用方法 + +### 基本使用 + +```bash +# 使用默认配置文件 +python main.py + +# 指定配置文件 +python main.py <配置文件路径>.toml +``` + +### 批量计算 + +批量计算不同保护角下的跳闸率： + +```bash +python main-batch.py <配置文件路径>.toml +``` + +结果将输出到 `r.txt` 文件中。 + +### 打包为可执行文件 + +使用 Makefile 打包： + +```bash +make +``` + +生成的可执行文件位于 `dist/Lightening.exe` + +## 配置文件格式 + +配置文件使用 TOML 格式，包含三个主要部分： + +### [parameter] - 基本参数 + +```toml +[parameter] +rated_voltage = 750 # 额定电压等级 (kV) +h_c_sag = 14.43 # 导线弧垂 (m) +h_g_sag = 11.67 # 地线弧垂 (m) +insulator_c_len = 7.02 # 导线串子绝缘长度 (m) +string_c_len = 9.2 # 导线串长 (m) +string_g_len = 0.5 # 地线串长 (m) +h_arm = [150, 130] # 导、地线挂点垂直距离 (m)，第一个值为地线挂点高度 +gc_x = [17.9, 17] # 导、地线水平坐标 (m) +ground_angels = [0] # 地面倾角 (°)，向下为正，支持多个角度 +altitude = 1000 # 海拔高度 (m) +td = 20 # 雷暴日 (d) +``` + +### [advance] - 高级参数 + +```toml +[advance] +ng = -1 # 地闪密度 (次/(km²·a))，大于0时使用此值，否则通过雷暴日计算 +Ip_a = -1 # 雷电流概率密度曲线系数a，大于0时使用此值 +Ip_b = -1 # 雷电流概率密度曲线系数b，大于0时使用此值 +``` + +注意：当 `ng` > 0 时，不会通过雷暴日计算地闪密度；当 `Ip_a` 和 `Ip_b` > 0 时，不会使用默认雷暴日对应的概率密度。 + +### [optional] - 可选参数 + +```toml +[optional] +voltage_n = 3 # 计算时电压分成多少份（考虑电压波动影响） +max_i = 200 # 最大尝试雷电流 (kA) +``` + +## 计算原理 + +### 击距模型 + +程序使用基于电气几何模型的方法计算绕击跳闸率，主要涉及以下击距公式： + +1. **地线击距**：$r_s = 10 \times I^{0.65}$ +2. **导线击距**：$r_c = 1.63 \times (5.015 \times I^{0.578} - 0.001 \times U_{ph})^{1.125}$ +3. **地面击距**： + - $h_{av} < 40m$: $r_g = (3.6 + 1.7 \ln(43 - h_{av})) \times I^{0.65}$ + - $h_{av} \ge 40m$: $r_g = 5.5 \times I^{0.65}$ + +### 地闪密度计算 + +根据 Q/GDW 11452-2015 导则： + +$N_g = 0.023 \times T_d^{1.3}$ + +### 跳闸率计算 + +通过积分暴露弧面积与雷电流概率密度的乘积得到最终跳闸率，并考虑建弧率。 + +## 输出结果 + +1. **控制台日志**：显示详细的计算过程和中间结果 +2. **DXF文件**：egm1.dxf、egm2.dxf 等，可视化击距模型 +3. **跳闸率结果**：单位为 次/(100km·a) + +## 项目结构 + +``` +EGM/ +├── main.py # 主程序入口 +├── main-batch.py # 批量计算程序 +├── core.py # 核心计算模块 +├── animation.py # 动画演示模块 +├── article.toml # 示例配置文件 +├── Makefile # 构建脚本 +├── pyproject.toml # 项目配置 +├── README.md # 说明文档 +└── CSharp/ # C# 版本实现 +``` + +## 技术支持 + +程序基于新版大手册公式和 Q/GDW 11452-2015《架空输电线路防雷导则》实现。 + +## 许可证 + +请参考项目许可证文件。 diff --git a/animation.py b/animation.py index 5f56906..b2bd9fc 100644 --- a/animation.py +++ b/animation.py @@ -29,8 +29,8 @@ class Animation: return wrapTheFunction - def disable(self, _disable): - self._disable = _disable + def enable(self, _enable): + self._disable = not _enable @switch_decorator def init_fig(self): diff --git a/article.toml b/article.toml index f2a94e5..9a4a476 100644 --- a/article.toml +++ b/article.toml @@ -3,13 +3,13 @@ title = "绕击跳闸率计算文件" rated_voltage = 750 # 额定电压等级 h_c_sag = 14.43 # 导线弧垂 h_g_sag = 11.67 # 地线弧垂 -insulator_c_len = 7.0 # 导线串子绝缘长度 +insulator_c_len = 7.02 # 导线串子绝缘长度 string_c_len = 9.2 # 导线串长 string_g_len = 0.5 # 地线串长 -h_arm = [100, 80] # 导、地线挂点垂直距离,计算的是中相 +h_arm = [150, 130] # 导、地线挂点垂直距离,计算的是中相 gc_x = [17.9, 17] # 导、地线水平坐标,计算的是中相 ground_angels = [0] # 地面倾角，向下为正，单位° -altitude = 1500 # 海拔，单位米 +altitude = 1000 # 海拔，单位米 td = 20 # 雷暴日 [advance] # ng=29.6 #地闪密度 ！！注意！！ 如果地闪密度大于0，则不会通过雷暴日计算地闪密度。填-1则忽略该项数据。 diff --git a/main-batch.py b/main-batch.py new file mode 100644 index 0000000..00b3c28 --- /dev/null +++ b/main-batch.py @@ -0,0 +1,365 @@ +import math +import os.path +import sys +import tomli +from loguru import logger +from core import * +import timeit + + +# 打印参数 +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: + 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"] + para.h_g_sag = toml_parameter["h_g_sag"] # 地线弧垂 + para.h_c_sag = toml_parameter["h_c_sag"] # 导线弧垂 + # para.h_whole = toml_parameter["h_whole"] # 杆塔全高 + para.td = toml_parameter["td"] # 雷暴日 + para.insulator_c_len = toml_parameter["insulator_c_len"] # 串子绝缘长度 + para.string_c_len = toml_parameter["string_c_len"] + para.string_g_len = toml_parameter["string_g_len"] + para.gc_x = toml_parameter["gc_x"] # 导、地线水平坐标 + para.ground_angels = [ + angel / 180 * math.pi for angel in toml_parameter["ground_angels"] + ] # 地面倾角，向下为正 + 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"ZC27103B-2000m.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) + h_whole = para.h_arm[0] # 塔全高 + string_g_len = para.string_g_len + string_c_len = para.string_c_len + h_g_sag = para.h_g_sag + h_c_sag = para.h_c_sag + gc_x = para.gc_x + shield_angle = [0, 10, 11, 12, 13, 14, 15, 16, 17, 18] + with open('r.txt','w',encoding='utf-8') as sa_result_file: + for sa in shield_angle: + gc_x[0] = gc_x[1] + math.tan(sa / 180 * math.pi) * (10 + 6 - 0.5) + h_arm = para.h_arm + gc_y = [ + h_whole - string_g_len - h_g_sag * 2 / 3, # 地线对地平均高 + ] + if len(h_arm) > 1: + for hoo in h_arm[1:]: + gc_y.append(hoo - string_c_len - h_c_sag * 2 / 3) # 导线平均高 + if len(gc_y) > 2: # 双回路 + phase_n = 3 # 边相导线数量 + else: + phase_n = 1 + # 地闪密度 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13 + td = para.td + ng = func_ng(td) + avr_n_sf = 0 # 考虑电压的影响计算的跳闸率 + ground_angels = para.ground_angels + for ground_angel in ground_angels: + logger.info(f"地面倾角{ground_angel/math.pi*180:.3f}°") + rg_type = None + rg_x = None + rg_y = None + cad = Draw() + voltage_n = para.voltage_n + n_sf_phases = np.zeros((phase_n, voltage_n)) # 存储每一相的跳闸率 + if np.any(np.array(gc_y) < 0): + logger.info("导线可能掉地面下了，程序退出。") + return 0 + for phase_conductor_foo in range(phase_n): + exposed_curve_shielded = False + rs_x = gc_x[phase_conductor_foo] + rs_y = gc_y[phase_conductor_foo] + rc_x = gc_x[phase_conductor_foo + 1] + rc_y = gc_y[phase_conductor_foo + 1] + if phase_n == 1: + rg_type = "g" + if phase_n > 1: # 多回路 + if phase_conductor_foo < 2: + rg_type = "c" # 捕捉弧有下面一相导线的击距代替 + rg_x = gc_x[phase_conductor_foo + 2] + rg_y = gc_y[phase_conductor_foo + 2] + else: + rg_type = "g" + # TODO 保护角公式可能有问题，后面改 + shield_angle_at_avg_height = ( + math.atan( + (rc_x - rs_x) + / ( + (h_arm[0] - string_g_len - h_arm[phase_conductor_foo + 1]) + + string_c_len + ) + ) + * 180 + / math.pi + ) # 挂点处保护角 + logger.info(f"挂点处保护角{shield_angle_at_avg_height:.3f}°") + logger.debug(f"最低相防护标识{rg_type}") + rated_voltage = para.rated_voltage + for u_bar in range(voltage_n): # 计算不同工作电压下的跳闸率 + # u_ph = ( + # math.sqrt(2) + # * 750 + # * math.cos(2 * math.pi / voltage_n * u_bar) + # / 1.732 + # ) # 运行相电压 + u_ph = rated_voltage + logger.info(f"计算第{phase_conductor_foo + 1}相，电压为{u_ph:.2f}kV") + # 迭代法计算最大电流 + i_max = 0 + insulator_c_len = para.insulator_c_len + # i_min = min_i(insulator_c_len, u_ph / 1.732) + i_min = min_i(insulator_c_len, u_ph) + _min_i = i_min # 尝试的最小电流 + _max_i = para.max_i # 尝试的最大电流 + # 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) + ): # 雷电流 + # logger.info(f"尝试计算电流为{i_bar:.2f}") + rs = rs_fun(i_bar) + rc = rc_fun(i_bar, u_ph) + rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type) + rg_line_func = None + if rg_type == "g": + rg_line_func = rg_line_function_factory(rg, ground_angel) + ####### + # 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) + ####### + rg_rc_circle_intersection = solve_circle_intersection( + rs, rc, rs_x, rs_y, rc_x, rc_y + ) + i_max = i_bar + if ( + not rg_rc_circle_intersection + ): # if circle_intersection is [] + logger.debug("保护弧和暴露弧无交点，检查设置参数。") + continue + circle_rc_line_or_rg_intersection = None + if rg_type == "g": + circle_rc_line_or_rg_intersection = ( + solve_circle_line_intersection( + rc, rc_x, rc_y, rg_line_func + ) + ) + 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: + # 暴露弧和捕捉弧无交点 + if rg_type == "g": + if rg_line_func(rc_x) > rc_y: + i_min = i_bar # 用于后面判断最小和最大电流是否相等，相等意味着暴露弧一直被屏蔽 + logger.info(f"捕捉面在暴露弧之上，设置最小电流为{i_min:.2f}") + else: + logger.info("暴露弧和地面捕捉弧无交点，检查设置参数。") + continue + else: + logger.info("上面的导地线无法保护下面的导地线，检查设置参数。") + continue + min_distance_intersection = ( + np.sum( + ( + np.array(rg_rc_circle_intersection) + - np.array(circle_rc_line_or_rg_intersection) + ) + ** 2 + ) + ** 0.5 + ) # 计算两圆交点和地面直线交点的最小距离 + if min_distance_intersection < 0.1: + break # 已经找到了最大电流 + # 判断是否以完全被保护 + 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, rs_x, rs_y, rg_line_func + ) # 保护弧和捕雷弧的交点 + ) + 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_or_rg_intersection) + - np.array([rc_x, rc_y]) + ) + ** 2 + ) + ** 0.5 + ) + if distance > rc: + logger.info(f"电流为{i_bar}kV时，暴露弧已经完全被屏蔽") + exposed_curve_shielded = True + break + # if phase_conductor_foo == 2: + cad.draw( + i_min, + u_ph, + rs_x, + rs_y, + rc_x, + rc_y, + rg_x, + rg_y, + rg_type, + ground_angel, + 2, + ) + cad.draw( + i_max, + u_ph, + rs_x, + rs_y, + rc_x, + rc_y, + rg_x, + rg_y, + rg_type, + ground_angel, + 6, + ) + cad.save_as(f"egm{phase_conductor_foo + 1}.dxf") + # 判断是否导线已经被完全保护 + if abs(i_max - _max_i) < 1e-5: + logger.info("无法找到最大电流，可能是杆塔较高。") + logger.info(f"最大电流设置为自然界最大电流{i_max}kA") + logger.info(f"最大电流为{i_max:.2f}") + logger.info(f"最小电流为{i_min:.2f}") + if exposed_curve_shielded: + logger.info("暴露弧已经完全被屏蔽,不会跳闸。") + continue + curt_fineness = 0.1 # 电流积分细度 + if i_min > i_max or abs(i_min - i_max) < curt_fineness: + logger.info("最大电流小于等于最小电流，没有暴露弧。") + continue + # 开始积分 + curt_segment_n = int((i_max - i_min) / curt_fineness) # 分成多少份 + i_curt_samples, d_curt = np.linspace( + 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, + u_ph, + rc_x, + rc_y, + rs_x, + rs_y, + rg_x, + rg_y, + ground_angel, + rg_type, + ) * 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]: + # cal_bd_first = bd_area(i_curt, u_ph, dgc, h_gav, h_cav) + # cal_bd_second = bd_area(i_curt + d_curt, u_ph, dgc, h_gav, h_cav) + # cal_thunder_density_first = thunder_density(i_curt) + # cal_thunder_density_second = thunder_density(i_curt + d_curt) + # calculus += ( + # ( + # cal_bd_first * cal_thunder_density_first + # + cal_bd_second * cal_thunder_density_second + # ) + # / 2 + # * d_curt + # ) + # if abs(calculus-0.05812740052770032)<1e-5: + # abc=123 + # pass + rated_voltage = para.rated_voltage + n_sf = ( + 2 + * ng + / 10 + * calculus + * arc_possibility(rated_voltage, insulator_c_len) + ) + avr_n_sf += n_sf / voltage_n + n_sf_phases[phase_conductor_foo][u_bar] = n_sf + logger.info(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.16f}次/(100km·a)") + logger.info(f"线路跳闸率是{avr_n_sf:.16f}次/(100km·a)") + logger.info( + f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=16)}次/(100km·a)" + ) + sa_result_file.write(f'{sa},{avr_n_sf}\n') + + +def speed(): + a = 0 + for bar in range(100000000): + a += bar + + +if __name__ == "__main__": + logger.remove() + logger.add(sys.stderr, level="DEBUG") + run_time = timeit.timeit("egm()", globals=globals(), number=1) + logger.info(f"运行时间:{run_time:.2f}s") + logger.info("Finished.") diff --git a/main.py b/main.py index 157586a..ae3a3e3 100644 --- a/main.py +++ b/main.py @@ -61,7 +61,7 @@ def read_parameter(toml_file_path): def egm(): if len(sys.argv) < 2: - toml_file_path = r"内自500kV-ZCK上相.toml" + toml_file_path = r"D:/code/EGM/历史/平乾750kV.toml" else: toml_file_path = sys.argv[1] if not os.path.exists(toml_file_path): @@ -96,7 +96,7 @@ def egm(): ground_angels = para.ground_angels # 初始化动画 animate = Animation() - animate.disable(False) + animate.enable(False) # animate.show() for ground_angel in ground_angels: logger.info(f"地面倾角{ground_angel/math.pi*180:.3f}°")