【高速入水】【CalcModel】模型的体积、质量、质心位置和转动惯量计算

模型的体积、质量、质心位置和转动惯量计算

函数描述

计算模型的体积、质量、质心位置和转动惯量。

参数

• $\pi$：圆周率。
• $L_{\text{mm}}$：模型总长度（mm）。
• $\text{Top0}$：顶点位置的无量纲量。
• $R_{\text{nmm}}$：空化器半径（mm）。
• $N_{\text{con}}$：模型节数。
• $\text{Part}$：每节的无量纲长度。
• $L_f$：前体长度（mm）。
• $\text{ConeLen}$：每节的长度（mm）。
• $\text{BaseDiam}$：每节的末端直径（mm）。
• $L_h$：腔体长度（mm）。
• $D_{Lh}$：腔体左直径（mm）。
• $D_{Rh}$：腔体右直径（mm）。
• $\rho_f$：前体材料密度（kg/m³）。
• $\rho_a$：后体材料密度（kg/m³）。
• $\rho_h$：腔体材料密度（kg/m³）。

返回值

• $TM_{\text{kg}}$：模型总质量（kg）。
• $X_c$：质心位置（与模型长度相关）。
• $I_c$：质心处的转动惯量（kg*mm²）。

计算公式

1. 初始化变量

• $\text{SumLen}$：每节的累积长度（mm）。
• $R$：每节的半径（mm）。
• $Q_{\text{fw}}$：前体各节的体积（mm³）。
• $Q_{\text{aw}}$：后体各节的体积（mm³）。
• $M_{\text{fw}}$：前体各节的质量（kg）。
• $M_{\text{aw}}$：后体各节的质量（kg）。
• $X_{\text{ci}}$：各节的质心位置（mm）。

2. 计算顶点位置和空化器体积

$\text{Top} = \text{Top0} \times L_{\text{mm}}$
$Q_n = \pi \times R_{\text{nmm}}^2 \times \text{Top} / 3$

3. 计算每节的长度和半径

$\text{SumLen}[i] = \text{Part}[i] \times L_{\text{mm}}$
$R[i] = \text{BaseDiam}[i] / 2$

4. 找到前体末端的锥体编号

$L_f > \text{SumLen}[i]$
$N_f = i$
$\text{NPf} = 2 \times i + 4$
$\text{NPa} = 2 \times N_{\text{con}} - \text{NPf} + 8$

5. 计算腔体参数

$R_{Lh} = D_{Lh} / 2$
$R_{Rh} = D_{Rh} / 2$
$L_{hx} = L_{\text{mm}} - L_h$
$R_{fh} = R_{Lh} + (L_f - L_{hx}) \times (R_{Rh} - R_{Lh}) / L_h$

6. 计算每节的体积

6.1 前体部分

$L_f \leq \text{ConeLen}[0]$
$R_f = R_{\text{nmm}} + L_f \times (R[0] - R_{\text{nmm}}) / \text{ConeLen}[0]$
$Q_{\text{fw}}[0] = \pi \times L_f \times (R_f^2 + R_{\text{nmm}}^2 + R_f \times R_{\text{nmm}}) / 3$
$Q_{\text{aw}}[0] = \pi \times (\text{ConeLen}[0] - L_f) \times (R[0]^2 + R_f^2 + R[0] \times R_f) / 3$

6.2 后体部分

$L_f > \text{SumLen}[j]$
$R_f = R[j-1] + (L_f - \text{SumLen}[j-1]) \times (R[j] - R[j-1]) / \text{ConeLen}[j]$
$Q_{\text{fw}}[j] = \pi \times (L_f - \text{SumLen}[j-1]) \times (R_f^2 + R[j-1]^2 + R_f \times R[j-1]) / 3$
$Q_{\text{aw}}[j] = \pi \times (\text{SumLen}[j] - L_f) \times (R[j]^2 + R_f^2 + R_f \times R[j]) / 3$

7. 计算每节的质量

$M_{\text{fw}}[i] = \rho_f \times Q_{\text{fw}}[i] / 1000$
$M_{\text{aw}}[i] = \rho_a \times Q_{\text{aw}}[i] / 1000$

8. 计算总体积

$Q_f = \sum Q_{\text{fw}}$
$Q_a = \sum Q_{\text{aw}}$
$Q = Q_f + Q_a + Q_n$

9. 计算腔体体积

$Q_h = \pi \times L_h \times (R_{Lh}^2 + R_{Lh} \times R_{Rh} + R_{Rh}^2) / 3$

10. 计算模型质量

$L_f \neq 0$
$M_n = \rho_f \times Q_n / 1000$
$M_f = \rho_f \times Q_f / 1000$
$M_a = \rho_a \times Q_a / 1000$
$M_h = \rho_h \times Q_h / 1000$

$TM_g = M_n + M_f + M_a$

$L_f \leq L_{hx}$
$TM_g = M_n + M_f + M_a - Q_h \times \rho_a / 1000$

11. 计算模型平均密度

$\rho = 1000 \times TM_g / Q$

12. 计算质心位置

$X_c = (X_{\text{cf}} \times (M_f + M_n) + X_{\text{ca}} \times M_a) / (M_f + M_n + M_a)$

13. 考虑腔体

$L_f \leq L_{hx}$
$X_c = (X_c \times (M_f + M_a + M_n) - X_{\text{ch}} \times (Q_h \times \rho_a / 1000)) / (M_n + M_f + M_a - (Q_h \times \rho_a / 1000))$

14. 计算转动惯量

$\beta = 180^\circ$
$I_{\text{of}} = 0$

$I_{\text{of}} = R_{\text{nmm}}^2 \times \text{Top}^3 \times (0.2 \times R_{\text{nmm}}^2 / \text{Top}^2 + 2 / 15)$

$I_{\text{of}} = \sum \left( \frac{1}{2} \times \text{BaseDiam}[i] \right)^2 \times \left( \left( \frac{1}{2} \times \text{BaseDiam}[i] \right)^2 \times \text{ConeLen}[i] + 4 \times (\text{Part}[i] \times L_{\text{mm}})^3 - (\text{Part}[i-1] \times L_{\text{mm}})^3 \right) / 3$

$I_{\text{of}} = I_{\text{of}} \times \pi \times \rho_f / 4 \times 10^{12}$

$I_{\text{oa}} = I_{\text{of}}$

$I_{\text{oa}} = \sum \left( \frac{1}{2} \times \text{BaseDiam}[i] \right)^2 \times \left( \left( \frac{1}{2} \times \text{BaseDiam}[i] \right)^2 \times \text{ConeLen}[i] + 4 \times (\text{Part}[i] \times L_{\text{mm}})^3 - (\text{Part}[i-1] \times L_{\text{mm}})^3 \right) / 3$

$I_{\text{oa}} = I_{\text{oa}} \times \pi \times \rho_a / 4 \times 10^{12}$

$I_o = I_{\text{of}} + I_{\text{oa}}$

$L_f = 0$
$I_o = I_{\text{oa}}$

$L_f = L_{\text{mm}}$
$I_o = I_{\text{of}}$

$I_c = I_o - TM_{\text{kg}} \times (X_c \times L_{\text{mm}})^2$

Python代码

def CalcModel(self):
    # PI, Lmm, Top0, Rnmm, Ncon, Part, Lf, ConeLen, BaseDiam, Lh, DLh, DRh, Rhof, Rhoa, Rhoh
    """
    计算模型的体积、质量、质心位置和转动惯量。

    参数:
    PI (float): 圆周率。
    Lmm (float): 模型总长度（mm）。
    Top0 (float): 顶点位置的无量纲量。
    Rnmm (float): 空化器半径（mm）。
    Ncon (int): 模型节数。
    Part (list of float): 每节的无量纲长度。
    Lf (float): 前体长度（mm）。
    ConeLen (list of float): 每节的长度（mm）。
    BaseDiam (list of float): 每节的末端直径（mm）。
    Lh (float): 腔体长度（mm）。
    DLh (float): 腔体左直径（mm）。
    DRh (float): 腔体右直径（mm）。
    Rhof (float): 前体材料密度（kg/m^3）。
    Rhoa (float): 后体材料密度（kg/m^3）。
    Rhoh (float): 腔体材料密度（kg/m^3）。

    返回:
    TMkg (float): 模型总质量（kg）。
    Xc (float): 质心位置（与模型长度相关）。
    Ic (float): 转动惯量（与质心相关，kg*mm^2）。
    """
    # 取回变量
    # 点集分辨率
    DPI = self.DPI
    # 模型的节数
    Ncon = self.Ncon
    # 圆周率，瞎注释啊哈哈哈哈
    PI = self.PI
    # 顶点的无量纲位置
    Top0 = self.Top0
    # 模型总长度
    Lmm = self.Lmm
    # 空化器半径，单位毫米
    Rnmm = self.Rnmm
    # 美节的无量纲长度
    Part = self.Part
    # 模型前体（Forebody）的长度
    Lf = self.Lf
    # 每节长度，单位毫米
    ConeLen = self.ConeLen
    # 每节末端直径
    BaseDiam = self.BaseDiam
    # 模型腔体长度，单位毫米
    Lh = self.Lh
    # 模型腔体左端直径
    DLh = self.DLh
    # 模型腔体右端直径
    DRh = self.DRh
    # 模型前体（Forebody）部分的材料密度 (float): 前体材料密度（kg/m^3）。
    Rhof = self.Rhof
    # Rhoa (float): 后体材料密度（kg/m^3）。
    Rhoa = self.Rhoa
    # Rhoh (float): 腔体材料密度（kg/m^3）。
    Rhoh = self.Rhoh
    # 空化器直径
    Dnmm = self.Dnmm
    # 空化器锥角
    Beta = self.Beta

    # 初始化变量
    # mm，这里就是将无量纲各节末端长度转换为mm量纲，并非求和。
    SumLen = np.zeros(DPI)
    # 每一节末端半径实际长度
    R = np.zeros(DPI)
    # 存储模型前体（Forebody）各节的体积
    Qfw: ndarray = np.zeros(DPI)
    # 存储模型后体（Aftbody）各节的体积
    Qaw = np.zeros(DPI)
    # 用于存储模型前体（Forebody）各节的质量
    Mfw = np.zeros(DPI)
    # 存储模型后体（Aftbody）各节的质量
    Maw = np.zeros(DPI)
    #  存储模型各节的质心位置
    Xci = np.zeros(DPI)

    # 计算
    Top = Top0 * Lmm  # 顶点位置（mm）
    Qn = PI * Rnmm ** 2 * Top / 3.0  # 空化器体积（mm^3）

    # 计算每节的长度和半径
    for i in range(Ncon):
        SumLen[i] = Part[i] * Lmm  # 每节的长度（mm）
        R[i] = BaseDiam[i] / 2.0  # 每节的半径（mm）

    # 找到前体末端的锥体编号
    i = 0
    # Lf 是前体长度
    while Lf > SumLen[i]:
        i += 1
    # 前体末端锥体编号
    Nf = i
    NPf = 2 * i + 4  # 前体多边形节点个数
    NPa = 2 * Ncon - NPf + 8  # 后体多边形节点个数

    RLh = DLh / 2.0  # 腔体左半径（mm）
    RRh = DRh / 2.0  # 腔体右半径（mm）
    Lhx = Lmm - Lh  # 无腔体部分长度（mm）
    if Lf > Lhx:
        Rfh = RLh + (Lf - Lhx) * (RRh - RLh) / Lh  # 前体端空心半径（mm）

    # 计算每节的体积
    if Lf <= ConeLen[0]:
        Rf = Rnmm + Lf * (R[0] - Rnmm) / ConeLen[0]
        Qfw[0] = PI * Lf * (Rf ** 2 + Rnmm ** 2 + Rf * Rnmm) / 3.0
        Qaw[0] = PI * (ConeLen[0] - Lf) * (R[0] ** 2 + Rf ** 2 + R[0] * Rf) / 3.0
        for i in range(1, Ncon):
            Qfw[i] = 0.0
            Qaw[i] = PI * ConeLen[i] * (R[i] ** 2 + R[i - 1] ** 2 + R[i] * R[i - 1]) / 3.0
    else:
        for j in range(1, Ncon):
            if Lf <= SumLen[j]:
                Rf = R[j - 1] + (Lf - SumLen[j - 1]) * (R[j] - R[j - 1]) / ConeLen[j]
                Qfw[0] = PI * ConeLen[0] * (R[0] ** 2 + Rnmm ** 2 + R[0] * Rnmm) / 3.0
                Qaw[0] = 0.0
                for i in range(1, j):
                    Qfw[i] = PI * ConeLen[i] * (R[i] ** 2 + R[i - 1] ** 2 + R[i] * R[i - 1]) / 3.0
                    Qaw[i] = 0.0
                Qfw[j] = PI * (Lf - SumLen[j - 1]) * (Rf ** 2 + R[j - 1] ** 2 + Rf * R[j - 1]) / 3.0
                Qaw[j] = PI * (SumLen[j] - Lf) * (R[j] ** 2 + Rf ** 2 + Rf * R[j]) / 3.0
                for i in range(j + 1, Ncon):
                    Qfw[i] = 0.0
                    Qaw[i] = PI * ConeLen[i] * (R[i] ** 2 + R[i - 1] ** 2 + R[i] * R[i - 1]) / 3.0
                break

    # 计算每节的质量
    for i in range(Ncon):
        Mfw[i] = Rhof * Qfw[i] / 1e3
        Maw[i] = Rhoa * Qaw[i] / 1e3

    Qf = sum(Qfw)  # 前体体积
    Qa = sum(Qaw)  # 后体体积
    Q = Qf + Qa + Qn  # 模型总体积

    # 计算腔体体积
    Qh = PI * Lh * (RLh ** 2 + RLh * RRh + RRh ** 2) / 3.0  # mm^3

    # 计算模型质量
    if Lf != 0.0:
        Mn = Rhof * Qn / 1e3
        Mf = Rhof * Qf / 1e3
        Ma = Rhoa * Qa / 1e3
        Mh = Rhoh * Qh / 1e3
    else:
        Mn = Rhoa * Qn / 1e3
        Mf = 0.0
        Ma = Rhoa * Qa / 1e3
        Mh = Rhoh * Qh / 1e3

    TMg = Mn + Mf + Ma  # 不考虑腔体的质量
    if Lf <= Lhx:
        if Rhoh == 0.0:
            TMg = Mn + Mf + Ma - Qh * Rhoa / 1e3
        else:
            TMg = Mn + Mf + Ma - Qh * Rhoa / 1e3 + Mh
    else:
        Qhf = PI * (Lf - Lhx) * (RLh ** 2 + Rfh * RLh + Rfh ** 2) / 3.0  # 前体部分空心体积
        Qha = PI * (Lmm - Lf) * (Rfh ** 2 + Rfh * RRh + RRh ** 2) / 3.0  # 后体部分空心体积
        if Rhoh == 0.0:
            TMg = Mn + Mf + Ma - (Qha * Rhoa / 1e3) - (Qhf * Rhof / 1e3)
        else:
            TMg = Mn + Mf + Ma - (Qha * Rhoa / 1e3) - (Qhf * Rhof / 1e3) + Mh

    # 计算模型平均密度
    Rho = 1e3 * TMg / Q

    # 确定中段截面
    Imid = 0
    for i in range(1, Ncon):
        if BaseDiam[i] > BaseDiam[i - 1]:
            Imid = i
    # Xp 在 CalcModel 中代表了模型的中段截面的相对位置，不知道是否共享了啊
    if Imid == Ncon:
        Xp = 1.0 - 4.0 * Q / (PI * Lmm * BaseDiam[Ncon] ** 2)
    elif BaseDiam[Imid] > BaseDiam[Ncon] and 1.0 - Part[Imid] < 0.02:
        Lu = 1.0 / Part[Imid]
        Xp = (0.67 - 0.75 * Lu + 0.2 * Lu ** 2) / (1.0 - Lu + 0.25 * Lu ** 2)
    elif BaseDiam[Imid] == BaseDiam[Ncon]:
        Lu = 1.0 / Part[Imid]
        work = 0.5 * (BaseDiam[0] - Dnmm) / ConeLen[0]
        La = Lmm * Part[Imid] / BaseDiam[Imid]
        Xp = (0.73 + 0.67 * La * (Lu ** 2 - 1.0)) / (Lu * (1.57 + 1.33 * La * (Lu ** 2 - 1.0)))

    # 计算前体的质心
    if Nf > 1:
        if BaseDiam[0] == Dnmm:
            Xci[0] = ConeLen[0] / 2.0
        else:
            hi = R[0] * ConeLen[0] / (R[0] - Rnmm)
            xc1 = SumLen[0] - hi / 4.0
            xc2 = - (hi - ConeLen[0]) / 4.0
            Q1 = PI * R[0] ** 2 * hi / 3.0
            Q2 = PI * Rnmm ** 2 * (hi - ConeLen[0]) / 3.0
            Xci[0] = (Q1 * xc1 - Q2 * xc2) / (Q1 - Q2)

        for i in range(1, Nf - 1):
            if BaseDiam[i] == BaseDiam[i - 1]:
                Xci[i] = SumLen[i] - ConeLen[i] / 2.0
            else:
                hi = R[i] * ConeLen[i] / (R[i] - R[i - 1])
                xc1 = SumLen[i] - hi / 4.0
                xc2 = SumLen[i] - ConeLen[i] - (hi - ConeLen[i]) / 4.0
                Q1 = PI * R[i] ** 2 * hi / 3.0
                Q2 = PI * R[i - 1] ** 2 * (hi - ConeLen[i]) / 3.0
                Xci[i] = (Q1 * xc1 - Q2 * xc2) / (Q1 - Q2)

        if Rf * 2.0 == BaseDiam[Nf - 1]:
            Xci[Nf] = Lf - (Lf - SumLen[Nf - 1]) / 2.0
        else:
            hi = Rf * (Lf - SumLen[Nf - 1]) / (Rf - R[Nf - 1])
            xc1 = Lf - hi / 4.0
            xc2 = Lf - (Lf - SumLen[Nf - 1]) - (hi - (Lf - SumLen[Nf - 1])) / 4.0
            Q1 = PI * Rf ** 2 * hi / 3.0
            Q2 = PI * R[Nf - 1] ** 2 * (hi - (Lf - SumLen[Nf - 1])) / 3.0
            Xci[Nf] = (Q1 * xc1 - Q2 * xc2) / (Q1 - Q2)

    else:
        if Rf * 2.0 == Dnmm:
            Xci[Nf] = Lf / 2.0
        else:
            hi = Rf * Lf / (Rf - Rnmm)
            xc1 = Lf - hi / 4.0
            xc2 = - (hi - Lf) / 4.0
            Q1 = PI * Rf ** 2 * hi / 3.0
            Q2 = PI * Rnmm ** 2 * (hi - Lf) / 3.0
            Xci[Nf] = (Q1 * xc1 - Q2 * xc2) / (Q1 - Q2)

    work1 = - Mn * Top / 4.0
    for i in range(Nf):
        work1 += Mfw[i] * Xci[i]
    Xcf = work1 / (Mf + Mn)

    # 计算后体的质心
    if BaseDiam[Nf] == Rf * 2.0:
        Xci[Nf] = SumLen[Nf] - (SumLen[Nf] - Lf) / 2.0
    else:
        hi = R[Nf] * (SumLen[Nf] - Lf) / (R[Nf] - Rf)
        xc1 = SumLen[Nf] - hi / 4.0
        xc2 = SumLen[Nf] - (SumLen[Nf] - Lf) - (hi - (SumLen[Nf] - Lf)) / 4.0
        Q1 = PI * R[Nf] ** 2 * hi / 3.0
        Q2 = PI * Rf ** 2 * (hi - (SumLen[Nf] - Lf)) / 3.0
        Xci[Nf] = (Q1 * xc1 - Q2 * xc2) / (Q1 - Q2)

    for i in range(Nf + 1, Ncon):
        if BaseDiam[i] == BaseDiam[i - 1]:
            Xci[i] = SumLen[i] - ConeLen[i] / 2.0
        else:
            hi = R[i] * ConeLen[i] / (R[i] - R[i - 1])
            xc1 = SumLen[i] - hi / 4.0
            xc2 = SumLen[i] - ConeLen[i] - (hi - ConeLen[i]) / 4.0
            Q1 = PI * R[i] ** 2 * hi / 3.0
            Q2 = PI * R[i - 1] ** 2 * (hi - ConeLen[i]) / 3.0
            Xci[i] = (Q1 * xc1 - Q2 * xc2) / (Q1 - Q2)

    work1 = 0.0
    for i in range(Nf, Ncon):
        work1 += Maw[i] * Xci[i]
    Xca = work1 / Ma

    # 计算腔体的质心
    if RLh == RRh:
        Xch = Lmm - Lh / 2.0
    else:
        hh = RRh * Lh / (RRh - RLh)
        xc1h = Lmm - hh / 4.0
        xc2h = Lmm - Lh - (hh - Lh) / 4.0
        Q1h = PI * RRh ** 2 * hh / 3.0
        Q2h = PI * RLh ** 2 * (hh - Lh) / 3.0
        Xch = (Q1h * xc1h - Q2h * xc2h) / (Q1h - Q2h)

    # 计算模型质心
    Xc = (Xcf * (Mf + Mn) + Xca * Ma) / (Mf + Mn + Ma)

    if Lf == 0.0:
        Xc = Xca
    elif Lf == Lmm:
        Xc = Xcf

    # 考虑腔体
    if (Lh > 0.0 and RLh >= 0.0 and RRh > 0.0) or (Lh > 0.0 and RLh > 0.0 and RRh >= 0.0):
        if Lf <= Lhx:
            Xc = (Xc * (Mf + Ma + Mn) - Xch * (Qh * Rhoa / 1e3)) / (Mn + Mf + Ma - (Qh * Rhoa / 1e3))
        elif Lf == Lmm:
            Xc = (Xcf * (Mf + Mn) - Xch * (Qh * Rhof / 1e3)) / (Mn + Mf - (Qh * Rhof / 1e3))
        else:
            if DLh == DRh:
                Xcha = Lmm - (Lmm - Lf) / 2.0
                Xchf = Lf - (Lf - Lhx) / 2.0
            else:
                hhf = Rfh * (Lf - Lhx) / (Rfh - RLh)
                xc1hf = Lf - hhf / 4.0
                xc2hf = Lf - (Lf - Lhx) - (hhf - (Lf - Lhx)) / 4.0
                Q1hf = PI * Rfh ** 2 * hhf / 3.0
                Q2hf = PI * RLh ** 2 * (hhf - (Lf - Lhx)) / 3.0
                Xchf = (Q1hf * xc1hf - Q2hf * xc2hf) / (Q1hf - Q2hf)

                hha = RRh * (Lmm - Lf) / (RRh - Rfh)
                xc1ha = Lmm - hha / 4.0
                xc2ha = Lmm - (Lmm - Lf) - (hha - (Lmm - Lf)) / 4.0
                Q1ha = PI * RRh ** 2 * hha / 3.0
                Q2ha = PI * Rfh ** 2 * (hha - (Lmm - Lf)) / 3.0
                Xcha: float = (Q1ha * xc1ha - Q2ha * xc2ha) / (Q1ha - Q2ha)

            Xc = (Xcf * (Mf + Mn) + Xca * Ma - Xcha * (Qha * Rhoa / 1e3) - Xchf * (Qhf * Rhof / 1e3)) / (
                        Mf + Ma + Mn - (Qha * Rhoa / 1e3) - (Qhf * Rhof / 1e3))

        if Rhoh != 0.0:
            if Lhx >= Lf:
                Xc = (Xc * (Mf + Ma + Mn - (Qh * Rhoa / 1e3)) + Xch * Mh) / (Mn + Mh + Mf + Ma - (Qh * Rhoa / 1e3))
            else:
                Xc = (Xc * (Mf + Mn + Ma - (Qha * Rhoa / 1e3) - (Qhf * Rhof / 1e3)) + Xch * Mh) / (
                            Mh + Mf + Mn + Ma - (Qha * Rhoa / 1e3) - (Qhf * Rhof / 1e3))

    Xc = Xc / Lmm  # 质心位置与模型长度相关

    # 计算转动惯量
    if Beta == 180.0:
        Iof = 0.0
    else:
        Iof = Rnmm ** 2 * Top ** 3 * (0.2 * Rnmm ** 2 / Top ** 2 + 2.0 / 15.0)

    # 前体
    # 如果不是第一个
    if Nf > 0:
        if BaseDiam[0] == Dnmm:
            Iof += Rnmm ** 2 * ConeLen[0] * (Rnmm ** 2 + 4.0 * ConeLen[0] ** 2 / 3.0)
        else:
            work2 = 0.5 * (BaseDiam[0] - Dnmm) / ConeLen[0]
            work3 = Rnmm
            work4 = work3 ** 2 * Lmm ** 3 * Part[0] ** 3 / 3.0 + work3 * work2 * Lmm ** 4 * Part[
                0] ** 4 / 2.0 + work2 ** 2 * Lmm ** 5 * Part[0] ** 5 / 5.0
            Iof += 0.2 * ((work3 + work2 * ConeLen[0]) ** 5 - work3 ** 5) / work2 + 4.0 * work4

        for i in range(1, Nf - 1):
            if BaseDiam[i] == BaseDiam[i - 1]:
                Iof += (0.5 * BaseDiam[i]) ** 2 * (
                            (0.5 * BaseDiam[i]) ** 2 * ConeLen[i] + 4.0 * (Part[i] * Lmm) ** 3 - (
                                Part[i - 1] * Lmm) ** 3) / 3.0
            else:
                work2 = 0.5 * (BaseDiam[i] - BaseDiam[i - 1]) / ConeLen[i]
                work3 = 0.5 * BaseDiam[i - 1] - Part[i - 1] * Lmm * work2
                work4 = work3 ** 2 * Lmm ** 3 * (
                            Part[i] ** 3 - Part[i - 1] ** 3) / 3.0 + work3 * work2 * Lmm ** 4 * (
                                    Part[i] ** 4 - Part[i - 1] ** 4) / 2.0 + work2 ** 2 * Lmm ** 5 * (
                                    Part[i] ** 5 - Part[i - 1] ** 5) / 5.0
                Iof += 0.2 * ((work3 + work2 * Part[i] * Lmm) ** 5 - (
                            work3 + work2 * Part[i - 1] * Lmm) ** 5) / work2 + 4.0 * work4

        if Rf * 2.0 == BaseDiam[Nf - 1]:
            Iof += (0.5 * 2.0 * Rf) ** 2 * (
                        (0.5 * 2.0 * Rf) ** 2 * (Lf - SumLen[Nf - 1]) + 4.0 * (Lf ** 3 - SumLen[Nf - 1] ** 3) / 3.0)
        else:
            work2 = 0.5 * (2.0 * Rf - BaseDiam[Nf - 1]) / (Lf - SumLen[Nf - 1])
            work3 = 0.5 * BaseDiam[Nf - 1] - SumLen[Nf - 1] * work2
            work4 = work3 ** 2 * (Lf ** 3 - SumLen[Nf - 1] ** 3) / 3.0 + work3 * work2 * (
                        Lf ** 4 - SumLen[Nf - 1] ** 4) / 2.0 + work2 ** 2 * (Lf ** 5 - SumLen[Nf - 1] ** 5) / 5.0
            Iof += 0.2 * ((work3 + work2 * Lf) ** 5 - (work3 + work2 * SumLen[Nf - 1]) ** 5) / work2 + 4.0 * work4

    else:
        if 2.0 * Rf == Dnmm:
            Iof += Rnmm ** 2 * Lf * (Rnmm ** 2 + 4.0 * Lf ** 2 / 3.0)
        else:
            work2 = 0.5 * (2.0 * Rf - Dnmm) / Lf
            work3 = Rnmm
            work4 = work3 ** 2 * Lf ** 3 / 3.0 + work3 * work2 * Lf ** 4 / 2.0 + work2 ** 2 * Lf ** 5 / 5.0
            Iof += 0.2 * ((work3 + work2 * Lf) ** 5 - work3 ** 5) / work2 + 4.0 * work4
    # 单位是kg*m^2
    Iof = Iof * PI * Rhof / 4e12

    # 后体
    if Lf == 0.0:
        if Beta == 180.0:
            Ioa = 0.0
        else:
            Ioa = Rnmm ** 2 * Top ** 3 * (0.2 * Rnmm ** 2 / Top ** 2 + 2.0 / 15.0)
    else:
        Ioa = 0.0

    if BaseDiam[Nf] == Rf * 2.0:
        Ioa += (0.5 * BaseDiam[Nf]) ** 2 * (
                    (0.5 * BaseDiam[Nf]) ** 2 * (SumLen[Nf] - Lf) + 4.0 * (SumLen[Nf] ** 3 - Lf ** 3) / 3.0)
    else:
        work2 = 0.5 * (BaseDiam[Nf] - 2.0 * Rf) / (SumLen[Nf] - Lf)
        work3 = 0.5 * (2.0 * Rf) - Lf * work2
        work4 = work3 ** 2 * (SumLen[Nf] ** 3 - Lf ** 3) / 3.0 + work3 * work2 * (
                    SumLen[Nf] ** 4 - Lf ** 4) / 2.0 + work2 ** 2 * (SumLen[Nf] ** 5 - Lf ** 5) / 5.0
        Ioa += 0.2 * ((work3 + work2 * SumLen[Nf]) ** 5 - (work3 + work2 * Lf) ** 5) / work2 + 4.0 * work4

    for i in range(Nf + 1, Ncon):
        if BaseDiam[i] == BaseDiam[i - 1]:
            Ioa += (0.5 * BaseDiam[i]) ** 2 * (
                        (0.5 * BaseDiam[i]) ** 2 * ConeLen[i] + 4.0 * (Part[i] * Lmm) ** 3 - (
                            Part[i - 1] * Lmm) ** 3) / 3.0
        else:
            work2 = 0.5 * (BaseDiam[i] - BaseDiam[i - 1]) / ConeLen[i]
            work3 = 0.5 * BaseDiam[i - 1] - Part[i - 1] * Lmm * work2
            work4 = work3 ** 2 * Lmm ** 3 * (Part[i] ** 3 - Part[i - 1] ** 3) / 3.0 + work3 * work2 * Lmm ** 4 * (
                        Part[i] ** 4 - Part[i - 1] ** 4) / 2.0 + work2 ** 2 * Lmm ** 5 * (
                                Part[i] ** 5 - Part[i - 1] ** 5) / 5.0
            Ioa += 0.2 * ((work3 + work2 * Part[i] * Lmm) ** 5 - (
                        work3 + work2 * Part[i - 1] * Lmm) ** 5) / work2 + 4.0 * work4

    Ioa = Ioa * PI * Rhoa / 4e12

    # 完整模型
    Io = Iof + Ioa

    if Lf == 0.0:
        Io = Ioa
    elif Lf == Lmm:
        Io = Iof

    # 考虑腔体
    if (Lh > 0.0 and RLh >= 0.0 and RRh > 0.0) or (Lh > 0.0 and RLh > 0.0 and RRh >= 0.0):
        if RRh == RLh:
            Iohw = (0.5 * DLh) ** 2 * ((0.5 * DRh) ** 2 * Lh + 4.0 * (Lmm ** 3 - Lhx ** 3) / 3.0)
        else:
            work2 = 0.5 * (DRh - DLh) / Lh
            work3 = 0.5 * DLh - Lhx * work2
            work4 = work3 ** 2 * (Lmm ** 3 - Lhx ** 3) / 3.0 + work3 * work2 * (
                        Lmm ** 4 - Lhx ** 4) / 2.0 + work2 ** 2 * (Lmm ** 5 - Lhx ** 5) / 5.0
            Iohw = 0.2 * ((work3 + work2 * Lmm) ** 5 - (work3 + work2 * Lhx) ** 5) / work2 + 4.0 * work4

        if Lf <= Lhx:
            Ioh = Iohw * PI * Rhoa / 4e12
        else:
            if RRh == RLh:
                Iohf = (0.5 * DLh) ** 2 * ((0.5 * DLh) ** 2 * (Lf - Lhx) + 4.0 * (Lf ** 3 - Lhx ** 3) / 3.0)
            else:
                work2 = 0.5 * (2.0 * Rfh - DLh) / (Lf - Lhx)
                work3 = 0.5 * DLh - Lhx * work2
                work4 = work3 ** 2 * (Lf ** 3 - Lhx ** 3) / 3.0 + work3 * work2 * (
                            Lf ** 4 - Lhx ** 4) / 2.0 + work2 ** 2 * (Lf ** 5 - Lhx ** 5) / 5.0
                Iohf = 0.2 * ((work3 + work2 * Lf) ** 5 - (work3 + work2 * Lhx) ** 5) / work2 + 4.0 * work4
            Iohf = Iohf * PI * Rhof / 4e12

            if RRh == RLh:
                Ioha = (0.5 * DLh) ** 2 * ((0.5 * DLh) ** 2 * (Lmm - Lf) + 4.0 * (Lmm ** 3 - Lf ** 3) / 3.0)
            else:
                work2 = 0.5 * (DRh - 2.0 * Rfh) / (Lmm - Lf)
                work3 = Rfh - Lf * work2
                work4 = work3 ** 2 * (Lmm ** 3 - Lf ** 3) / 3.0 + work3 * work2 * (
                            Lmm ** 4 - Lf ** 4) / 2.0 + work2 ** 2 * (Lmm ** 5 - Lf ** 5) / 5.0
                Ioha = 0.2 * ((work3 + work2 * Lmm) ** 5 - (work3 + work2 * Lf) ** 5) / work2 + 4.0 * work4
            Ioha = Ioha * PI * Rhoa / 4e12
            Ioh = Iohf + Ioha

        if Rhoh == 0.0:
            Io = Io - Ioh
        else:
            Io = Io - Ioh + Iohw * PI * Rhoh / 4e12

    # 计算总质量
    TMkg = TMg / 1e3  # kg

    # 计算质心处的转动惯量
    Ic = Io - TMkg * (Xc * Lmm) ** 2  # kg*mm^2
    ###################
    # Imid 好像定义了
    # 公共变量Imid
    self.Imid = Imid

    # 总质量
    self.TMkg = TMkg
    # 模型质心位置
    self.Xc = Xc
    # Xp 在 CalcModel 中代表了模型的中段截面的相对位置
    # 不知道是否共享了啊
    self.Xp = Xp
    # 模型质心转动惯量
    self.Ic = Ic
    # 模型总质量
    self.TMg = TMg
    # 模型平均密度
    self.Rho = Rho
    # 模型的总转动惯量
    self.Io = Io
    # 前体末端锥体编号
    self.Nf = Nf
    # 前体多边形节点个数
    self.NPf = NPf
    # 后体多边形节点个数
    self.NPa = NPa