egm/main.py

import math
import ezdxf
import numpy as np

gCAD = None
gMSP = None


class Draw:
    def __init__(self):
        self._doc = ezdxf.new(dxfversion="R2010")
        self._doc.layers.add("EGM", color=2)
        global gCAD
        gCAD = self

    def draw(self, i_curt, u_ph, h_gav, h_cav, dgc, 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), (200, rg), 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
        doc.saveas("egm.dxf")


# 圆交点
def solve_circle_intersection(rs, rc, hgav, hcav, dgc):
    # 用牛顿法求解
    x = 300
    y = 300
    for bar in range(0, 10):
        A = [[-2 * x, -2 * (y - hgav)], [-2 * (x - dgc), -2 * (y - hcav)]]
        b = [
            x ** 2 + (y - hgav) ** 2 - rs ** 2,
            (x - dgc) ** 2 + (y - hcav) ** 2 - rc ** 2,
        ]
        X_set = np.linalg.solve(A, b)
        x += X_set[0]
        y += X_set[1]
        if np.all(np.abs(X_set) < 1e-5):
            return [x, y]
    return []


# 圆与地面线交点
def solve_circle_line_intersection(rc, rg, h_cav, dgc):
    # TODO: 需要考虑地面捕雷线与暴露弧完全没交点的情况
    r = (rc ** 2 - (rg - h_cav) ** 2) ** 0.5 + dgc
    return [r, rg]


def min_i(string_len, u_ph):
    u_50 = 530 * string_len + 35
    z_0 = 300  # 雷电波阻抗
    z_c = 251  # 导线波阻抗
    r = (u_50 + 2 * z_0 / (2 * z_0 + z_c) * u_ph) * (2 * z_0 + z_c) / (z_0 * z_c)
    return r


def thunder_density(i):  # l雷电流幅值密度函数
    r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44)
    return r


def angel_density(angle):  # 入射角密度函数 angle单位是弧度
    r = 0.75 * (math.cos(angle) ** 3)
    return r


def rs_fun(i):
    r = 10 * (i ** 0.65)
    return r


def rc_fun(i, u_ph):
    r = 1.63 * ((5.015 * (i ** 0.578) - 0.001 * u_ph) ** 1.125)
    return r


def rg_fun(i, h_cav):
    if h_cav < 40:
        rg = (3.6 + 1.7 ** math.log(43 - h_cav)) ** 0.65
    else:
        rg = 5.5 * (i ** 0.65)
    return rg


def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph):  # 暴露弧的角度
    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, h_cav, dgc
    )  # 暴露圆和补雷线的交点
    np_circle_intersection = np.array(circle_intersection)
    theta2_line = np_circle_intersection - np.array([dgc, h_cav])
    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])
    theta1 = math.atan(theta1_line[1] / theta1_line[0])
    # 考虑雷电入射角度，所以theta1可以小于0，即计算从侧面击中的雷
    # if theta1 < 0:
    #     # print(f"θ_1角度为负数{theta1:.4f},人为设置为0")
    #     theta1 = 0
    return np.array([theta1, theta2])


def distance_point_line(point_x, point_y, line_x, line_y, k):
    d = abs(k * point_x - point_y - k * line_x + line_y) / ((k ** 2 + 1) ** 0.5)
    return d


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)
    rc = rc_fun(i_curt, u_ph)
    rs = rs_fun(i_curt)
    rg = rg_fun(i_curt, h_cav)
    # 求暴露弧上一点的切线
    line_x = math.cos(theta1) * rc + dgc
    line_y = math.sin(theta1) * rc + h_cav
    max_w = 0  # 入射角
    if theta1 < 0:
        max_w = theta1 + math.pi / 2
        k = math.tan(max_w)
        # 求保护弧到直线的距离，判断是否相交
        d_to_rs = distance_point_line(0, h_gav, 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)
            max_w = math.atan(new_k)  # 用于保护弧相切的角度
            intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph)
            gMSP.add_circle((0, h_gav), rs)
            gMSP.add_circle((dgc, h_cav), rc)
            gMSP.add_line((dgc, h_cav), (line_x, line_y))
            gMSP.add_line(
                (-500, k * (-500 - line_x) + line_y), (500, k * (500 - line_x) + line_y)
            )
            gMSP.add_line((0, rg), (1000, rg))
            gCAD.save()
            pass

    # k=tangent_line_k(point_x, point_y, dgc, h_cav,rc)
    # 暂时不考虑雷电入射角的影响
    r = (math.cos(theta1) - math.cos(theta2)) * rc
    return r
    # r1=rc*(-math.cos(thyta2)+math.cos(thyta1))
    # 入射角密度函数积分
    # arrival_angle_fineness=0.0001
    # for calculus_arv_angle in np.linspace()


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是圆心
    # TODO:应该找到两个角度值后再比较
    k = init_k
    for bar in range(0, 30):
        fk = (k * center_x - center_y - k * line_x + line_y) ** 2 - (radius ** 2) * (
            k ** 2 + 1
        )

        d_fk = (
            2 * (k * center_x - center_y - k * line_x + line_y) * (center_x - line_x)
            - 2 * (radius ** 2) * k
        )
        d_k = -fk / d_fk
        k += d_k
        if abs(d_k) < 1e-5:
            dd = distance_point_line(center_x, center_y, line_x, line_y, k)
            if abs(dd - radius) < 1e-5:
                return k
    return None


def egm():
    u_ph = 750 / 1.732  # 运行相电压
    h_cav = 160  # 导线对地平均高
    h_gav = h_cav + 9.5 + 2.7
    dgc = -2  # 导地线水平距离
    # 迭代法计算最大电流
    i_max = 0
    _min_i = 20  # 尝试的最小电流
    _max_i = 80  # 尝试的最大电流
    for i_bar in np.linspace(_min_i, _max_i, int((_max_i - _min_i) / 0.01)):  # 雷电流
        print(f"尝试计算电流为{i_bar:.2f}")
        rs = rs_fun(i_bar)
        if not np.isreal(rs):
            continue
        rc = rc_fun(i_bar, u_ph)
        if not np.isreal(rc):
            continue
        rg = rg_fun(i_bar, h_cav)
        if not np.isreal(rg):
            continue
        circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
        if not circle_intersection:  # if circle_intersection is []
            continue
        circle_line_intersection = solve_circle_line_intersection(rc, rg, h_cav, dgc)
        min_distance_intersection = (
            np.sum(
                (np.array(circle_intersection) - np.array(circle_line_intersection))
                ** 2
            )
            ** 0.5
        )  # 计算两圆交点和地面直线交点的最小距离
        i_max = i_bar
        if min_distance_intersection < 0.1:
            break
    i_min = min_i(6.78, 750 / 1.732)
    cad = Draw()
    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.save()
    if abs(i_max - _max_i) < 1e-5:
        print("无法找到最大电流，可能是杆塔较高。")
        i_max = 300
        print(f"最大电流设置为自然界最大电流{i_max}kA")
    print(f"最大电流为{i_max:.2f}")
    print(f"最小电流为{i_min:.2f}")
    if i_min > i_max:
        print("最大电流小于最小电流，没有暴露弧，程序结束。")
        return
    # 开始积分
    curt_fineness = 0.1  # 电流积分细度
    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)
    for i_curt in i_curt_samples:
        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
        )
    n_sf = 2 * 2.7 / 10 * calculus  # 调整率
    print(f"跳闸率是{n_sf:.6}")

    # draw(rs, rc, rg, h_gav, h_cav, dgc)


if __name__ == "__main__":
    thunder_density(2)
    egm()
    print("Finished.")