根据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)找到其中画样条曲线的命令,按照实际需要对其进行修改使用。