DARPA SUBOFF潜艇建模（一）

根据suboff几何特征报告，运用Python求suboff模型型值点

---

一、艇体建模

艇体为一轴对称的回转体，根据回转体母线方程，求得母线型值点，以x轴为回转轴，建立回转体。根据suboff几何特征报告，艇体分为Bow、Parallel middle body、Afterbody、Afterbody Cap四个部分，其方程分别为：

将其用代码表达，为：

import numpy as np
import matplotlib.pyplot as plt

Forebody_Length = 3.333333  # 前体长
Parallel_Middle_Body_Length = 7.3125  # 平行中体长
Afterbody_Length = 3.645833  # 后体长
Aft_Perpendicular_loc = 13.979167  # 尾柱位置
Total_Body_Length = 14.291666  # 总长
Maximum_Body_Diameter = 1.666667  # （未使用的参数）
Lambda = 24  # 缩尺比
R_max = 5 / 6  # 回转体最大半径


def equation(forebody_length, parallel_middle_body_length, aft_perpendicular_loc, total_body_length):
    # 首部方程
    x_bow = np.arange(0, forebody_length, 0.01)
    R_bow = R_max * (1.126395101 * x_bow * (0.3 * x_bow - 1) ** 4 + 0.442874707 * x_bow ** 2 *
                     (0.3 * x_bow - 1) ** 3 + 1 - (0.3 * x_bow - 1) ** 4 * (1.2 * x_bow + 1)) ** (1 / 2.1)
    plt.plot(x_bow, R_bow)

    # 平行中体方程
    x_parallel = np.arange(forebody_length, forebody_length + parallel_middle_body_length, 0.01)
    R_parallel = np.full(x_parallel.shape, R_max)
    plt.plot(x_parallel, R_parallel)

    # 尾部方程
    r_h, K_o, K_1 = 0.1175, 10, 44.6244
    x_afterbody = np.arange(forebody_length + parallel_middle_body_length, aft_perpendicular_loc, 0.01)
    eplison = (aft_perpendicular_loc - x_afterbody) / forebody_length
    R_afterbody = R_max * (r_h ** 2 + r_h * K_o * eplison ** 2 +
                           (20 - 20 * r_h ** 2 - 4 * r_h * K_o - 1 / 3 * K_1) * eplison ** 3 +
                           (-45 + 45 * r_h ** 2 + 6 * r_h * K_o + K_1) * eplison ** 4 +
                           (36 - 36 * r_h ** 2 - 4 * r_h * K_o - K_1) * eplison ** 5 +
                           (-10 + 10 * r_h ** 2 + r_h * K_o + 1 / 3 * K_1) * eplison ** 6) ** 0.5
    plt.plot(x_afterbody, R_afterbody)

    # 尾部cap方程
    x_afterbody_cap = np.arange(aft_perpendicular_loc, total_body_length, 0.01)
    R_afterbody_cap = 0.1175 * R_max * (1 - (3.2 * x_afterbody_cap - 44.733333) ** 2) ** 0.5
    plt.plot(x_afterbody_cap, R_afterbody_cap)
    plt.show()

    # 导出数据文件
    tmp_1 = np.concatenate((x_bow.reshape(-1, 1), R_bow.reshape(-1, 1)), axis=1)
    tmp_2 = np.concatenate((x_parallel.reshape(-1, 1), R_parallel.reshape(-1, 1)), axis=1)
    tmp_3 = np.concatenate((x_afterbody.reshape(-1, 1), R_afterbody.reshape(-1, 1)), axis=1)
    tmp_4 = np.concatenate((x_afterbody_cap.reshape(-1, 1), R_afterbody_cap.reshape(-1, 1)), axis=1)
    tmp = np.concatenate((tmp_1, tmp_2, tmp_3, tmp_4), axis=0)
    np.savetxt('艇体数据文件.txt', tmp, fmt='%.3f', delimiter=',')

equation(Forebody_Length, Parallel_Middle_Body_Length, Aft_Perpendicular_loc, Total_Body_Length)

运行上述代码，生成回转体母线型值点，如下图所示。

同时，在程序运行路径下生成型值点文件，打开文件，复制所有内容，同时在autocad里点击“样条曲线拟合”命令，将上述复制的内容粘贴至autocad命令框里，生成样条曲线，保存为dwg文件，导入catia，进而根据该样条曲线建立回转体。

二、尾部附体建模

代码如下（示例）：
尾部附体截面为如下所示的翼型，其中横坐标为无因次x位置。

截面方程为：

可按以下过程求解各截面型值点：
（1）先确定无因次x位置序列，在下边步骤正是求解在这些无因次x位置处的截面宽度。
（2）求尾部附体与艇体的交线（斜截面），在每个x位置处不断加减截面所处的高度，使得截面的点逐渐逼近艇体。
（3）以固定间隔求上部完整剖面，尾部附体最高位置为艇体最大半径处。
（4）根据求得的各截面，建立一个尾部附体，另外三个通过第一个旋转而来。
代码实现过程如下：

import matplotlib.pyplot as plt
import numpy as np


class SternAppendage:
    # 译自FORTRAN源码，变量保留了原来的名字
    def __init__(self):
        self.XXI = [0, 0.005, 0.0125, 0.025, 0.05, 0.075, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95,
                    1.0]  # 无因次x位置
        self.RH = 0.1175
        self.AK0 = 10
        self.AK1 = 44.6244
        self.NP = 19
        self.RMAX = 5 / 6
        self.SRS = ''
        self.Z = ''
        self.ITR = ''
        self.RHAS = ''
        self.RHA = ''
        self.X = ''
        self.HH = 13.146284
        self.CY = ''
        self.XI = ''
        self.RBSMAX = ''
        self.RR = ''
        self.RO = ''
        self.DELR = ''
        self.eplison = ''
        self.spacing = ''
        self.spacingX = ''

        self.PointTmp = []
        self.PointTmp_neg = []
        self.USet = []

    def half_sec(self):
        # 求交线（斜剖面）
        for J in range(self.NP):
            self.XI = self.XXI[J]
            self.RR = 0.075
            self.DELR = 0.025
            self.RR += self.DELR
            self.ITR = 0
            while True:
                self.CY = -0.466308 * self.RR + 0.88859
                self.X = (self.XI - 1.0) * self.CY + self.HH

                eplison = (13.979167 - self.X) / 3.333333
                self.RHA = self.RMAX * (self.RH ** 2 + self.RH * self.AK0 * eplison ** 2 +
                                        (20 - 20 * self.RH ** 2 - 4 * self.RH * self.AK0 - 1 / 3 * self.AK1) * eplison ** 3 +
                                        (-45 + 45 * self.RH ** 2 + 6 * self.RH * self.AK0 + self.AK1) * eplison ** 4 +
                                        (36 - 36 * self.RH ** 2 - 4 * self.RH * self.AK0 - self.AK1) * eplison ** 5 +
                                        (-10 + 10 * self.RH ** 2 + self.RH * self.AK0 + 1 / 3 * self.AK1) * eplison ** 6) ** 0.5
                self.RHAS = self.RHA ** 2

                # 定义尾部附体
                self.Z = 0.2969 * self.XI ** 0.5 - 0.126 * self.XI - 0.3516 * self.XI ** 2 + 0.2852 * self.XI ** 3 - 0.1045 * self.XI ** 4
                self.Z = self.CY * self.Z
                self.SRS = self.RR ** 2 + self.Z ** 2
                # 附体在艇体表面
                if abs(self.SRS - self.RHAS) <= 0.00001:
                    self.PointTmp.append(f'{self.X * 1000},{self.Z * 1000},{self.RR * 1000}')
                    #if self.XI != 0:剖面为左右对称的，如果想将两边曲线组合为一条曲线，则不加入重复点
                    self.PointTmp_neg.append(f'{self.X * 1000},{-self.Z * 1000},{self.RR * 1000}')
                    if self.XI == 0:
                        self.RBSMAX = self.RR
                    break
                self.ITR += 1
                if self.ITR > 20:
                    self.DELR *= 0.5  # 半径增量减半

                if self.SRS > self.RHAS:
                    self.RR -= self.DELR
                if self.SRS <= self.RHAS:
                    self.RR += self.DELR

        self.PointTmp.reverse()
        self.USet.append(self.PointTmp.copy())
        self.USet.append(self.PointTmp_neg.copy())

        self.PointTmp.clear()
        self.PointTmp_neg.clear()

    def sec(self):
        self.DELR = 0.05  # 半径增量
        for I in range(self.NP):
            self.RO = self.RR
            self.RR = self.RBSMAX + I * self.DELR
            if self.RR > self.RMAX:
                self.RR = self.RMAX
            if self.RR == self.RO:
                return
            self.CY = -0.466308 * self.RR + 0.88859
            for J in range(self.NP):
                self.XI = self.XXI[J]
                self.X = (self.XI - 1.0) * self.CY + self.HH
                self.Z = 0.2969 * self.XI ** 0.5 - 0.126 * self.XI - \
                         0.3516 * self.XI ** 2 + 0.2852 * self.XI ** 3 - 0.1045 * self.XI ** 4
                self.Z = self.CY * self.Z

                self.PointTmp.append(f'{self.X * 1000},{self.Z * 1000},{self.RR * 1000}')
                #if self.XI != 0:剖面为左右对称的，如果想将两边曲线组合为一条曲线，则不加入重复点
                self.PointTmp_neg.append(f'{self.X * 1000},{-self.Z * 1000},{self.RR * 1000}')

            self.PointTmp.reverse()
            self.USet.append(self.PointTmp.copy())
            self.USet.append(self.PointTmp_neg.copy())

            self.PointTmp.clear()
            self.PointTmp_neg.clear()


ste = SternAppendage()
ste.half_sec()
ste.sec()
# ste.Uset里存放了各截面型值点，绘制成图像观察
fig = plt.figure()
ax1 = fig.add_subplot(111, projection='3d')
for i in range(len(ste.USet)):
    point = []
    for j in range(len(ste.USet[i])):
        tmp = ste.USet[i][j].split(',')
        tmp = [float(k) for k in tmp]
        point.append(tmp)
    point = np.array(point)
    ax1.plot(point[:, 0], point[:, 1], point[:, 2], 'o-')
plt.show()

运行程序，得各剖面形状为：

剖面为对称的，一般仅需其一半数据，在上述代码中注释掉所有包含“self.PointTmp_neg”的语句既可。结果如下：

Uset里存放了各截面型值点，将其导出进行建模既可。

若想省去在各个软件间的来回操作，在上述代码里添加对catia的操作，直接根据点建立样条曲线。大致流程为：
（1）若对catia的COM接口不熟悉，可先在catia里录制宏，宏的内容为随便画一条样条曲线。
（2）导入操作COM接口的库，并与catia建立连接，如下：

import win32com.client
catia = win32com.client.Dispatch('CATIA.application')

（3）查看宏的代码，为VBA语言，将其复制到Python编辑环境里修改，一般来说需要修改的是，删去参数定义语句、将索引由圆括号改为方括号、将for循环等改为Python格式，VBA里的Sub类似于Python里的def。
（3）找到其中画样条曲线的命令，按照实际需要对其进行修改使用。