提交。

This commit is contained in:
dmy
2026-03-02 18:18:46 +08:00
parent 7a5bb05f58
commit a153e69eb7
6 changed files with 533 additions and 8 deletions

View File

@@ -11,4 +11,4 @@ build:
.PHONY:clean
clean:
rm -fr build
rm -fr build

160
README.md Normal file
View File

@@ -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《架空输电线路防雷导则》实现。
## 许可证
请参考项目许可证文件。

View File

@@ -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):

View File

@@ -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则忽略该项数据。

365
main-batch.py Normal file
View File

@@ -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
# 地闪密度 利用QGDW 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.")

View File

@@ -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}°")