DARPA SUBOFF潜艇建模（二）

Python求解围壳型值

---

一、围壳方程

二、代码

import numpy as np
import matplotlib.pyplot as plt
fig = plt.figure()
fig.add_subplot(projection="3d")

forebody_len = 0.325521  # 围壳前体长
parallel_middle_body_len = 0.200521  # 围壳中段长
afterbody_len = 0.682292  # 围壳后段长
total_len = forebody_len + parallel_middle_body_len + afterbody_len
deltx = 0.01  # 剖面间隔，等间隔的剖面其实并不好，下边没有应用这个
delty = 0.001  # 每个剖面高度方向点间隔
z_max = 0.109375  # 俯视围壳截面最大宽度
y_max = 1.507813  # 围壳主体最高平面，也是围壳cap下边缘
location_start = 3.032986  # 围壳开始位置


# forebody
def z_f(x_f):
    D = 3.072000 * (x_f - location_start)
    A = 2 * D * (D - 1) ** 4
    B = 1 / 3 * D ** 2 * (D - 1) ** 3
    C = 1 - (D - 1) ** 4 * (4 * D + 1)
    z_1 = z_max * (2.094759 * A + 0.2071781 * B + C) ** (1 / 2)
    return z_1


# parallel
def z_p(x_p):
    return z_max


# afterbody
def z_a(x_a):
    E = (location_start + total_len - x_a) / afterbody_len
    z_3 = z_max * (2.238361 * (E * (E - 1) ** 4) +
                   3.106529 * (E ** 2 * (E - 1) ** 3) + (1 - (E - 1) ** 4 * (4 * E + 1)))
    return z_3


def draw(x):
    zset = []
    yset = []
    z_2set = []
    for i in range(len(x)):
        if location_start <= x[i] < location_start + forebody_len:
            if x[i] == location_start:
                z = 0
            else:
                z = z_f(x[i])
        elif location_start + forebody_len <= x[i] < location_start + forebody_len + parallel_middle_body_len:
            z = z_p(x[i])
        elif location_start + forebody_len + parallel_middle_body_len <= x[i] <= location_start + total_len:
            if abs(x[i] - (location_start + total_len)) < 0.001:
                z = 0
            else:
                z = z_a(x[i])
        try:
            y = np.arange(y_max, y_max + z / 2, delty)
            y = np.append(y, [y_max + z / 2], axis=0)
        except:
            y = np.array([y_max])
        z_2 = (z ** 2 - (2 * (y - y_max)) ** 2)
        z_2 = np.where(z_2 < 0, 0, z_2)
        z_2 = z_2 ** (1 / 2)
        z_2[0] = z
        plt.plot(np.full(y.shape, x[i]), z_2, y, "o-")
        zset.append(z)  # 存放了每个x位置处的围壳主体一半宽度
        yset.append(y)  # 存放了每个x位置处的围壳cap截面高度
        z_2set.append(z_2)  # 存放了每个x位置处的围壳cap截面一半宽度
    with open("sj.txt", "w") as tt:
        tt.write("StartLoft\n")
        for i in range(len(z_2set)):
            tt.write("StartCurve\n")
            for j in range(len(z_2set[i])):
                tt.write(f"{round(x[i]*1000,3)}\t{round(z_2set[i][j]*1000,3)}\t{round(yset[i][j]*1000,3)}\n")
            tt.write("EndCurve\n")
        tt.write("EndLoft\n")
        tt.write("End\n")


x = [location_start, 3.039, 3.051, 3.070, 3.095, 3.119, 3.156, 3.259,
      3.599, 3.748, 3.831, 3.908, 3.995, 4.069, 4.155, 4.223, location_start + total_len]
# 下面是按等间隔剖面生成的型值，效果不好不建议使用
# x = np.arange(location_start, location_start + total_len, deltx)
# x = np.append(x, [location_start + total_len], axis=0)

draw(x)
plt.show()

程序生成的点如下图所示：

生成的数据文件“sj.txt”位于程序运行同目录下，里面包含了围壳Cap所有剖面的点，可以通过catia安装目录下的GSD_PointSplineLoftFromExcel.xls中的脚本导入catia，导入时由于开始和结束位置剖面都仅有1个点，可以先删掉，再导入剩余曲线，然后在catia里手工建这两个点。导入的曲线如下图所示：

手工建立首尾两点后，以中纵剖线为引导线建立多截面曲面，再将首尾的洞填充完整，对称后即获得围壳Cap。

Cap所有截面最下方的点构成了围壳主体轮廓，将其拉伸，再用艇体截断，就获得了围壳主体。