雷浚同学的ABAQUS的PYTHON二次开发——短线程网壳
短程线网壳的划分方式
短程线网壳的划分方式有多种：弦分法、弧分法、等弧再分法、椭圆投影法等等。弦分法的做法最为简单，即将20面体的基本三角形各边等分若干点，作划分线平行于该三角形的边，形成三角形网格，再将各点投影到外接球面上，连接球面上各点，得到短程线型球面网格。但得到的网格的质量不高，但在此次建模过程中采用弦分法的同时对网格质量进行了优化。
球形短线程网壳
短程线网壳网点的坐标计算
根据弦分法的划分规则，短程线网壳网点坐标计算时以球正二十面体的一个球面三角形为母体，进而根据所需的划分频率进行下级划分，因此在计算坐标时以一个球面三角形和划分出来的下级网点为基本单元，将这个基本单元进行在球内以圆心旋转复制即可获得一个完整的球的短程线网壳网点坐标。
但实际上，工程上并不需要一个完整球的网壳。通常情况下，根据场地的跨度和最大高度需要建立一个半球或者更小弧度的曲面。因此短程线网壳的更具体化参数化建模要求为输入跨度s和高度h以及用户所需要的划分频率即可自动建立网壳模。因此问题变为给出球的任意一个截面体，在上面进行合理的短程线网格划分。然而按照要求单纯截断短程线球体仅仅只能得出由边界杆件不完整的网壳，如图。
截断球体产生的不完整杆件
但是我们是否可以截取短程线网壳的完整某一部分进行拉伸调整变成适应于所要求曲面的短程线网壳？答案是可以的，同时我们也可以对截取部分进行限制性选择。 经过观察和比较，截取部分由两个部分组成：有5个大球面三角形组成的上部和多个小球面三角形组成的下部。
短线程网壳上部
短线程网壳下部
定义上部的划分频率为f1和下部的划分频率为f2（对应于网格层数，左图6层、右图3层），由于上部和下部的网格大小需要大致相同，因此对应f1：f2的值即可大致代表为上部和下部在截面体所占角度的比值α1：α2。
参数示意图
所以根据大致分配的角度进行上部和下部坐标计算，同时注意上部和下部连接节点的坐标协调即可完成短程线网格坐标计算。 最后在py脚本文件使用“p.WirePolyLine(points=point1,point2),meshable=ON)”语句依次生成连接网点的杆件就完成参数化建模。
自动优化网壳质量
然而在生成的短程线网格由于上下部占的角度为估算杆件的均匀度并不是高，因此对α1：α2的值进行适当调整。引入系数对角度比值进行变动，每一次变动记录网格杆件的均方差，在有限次变动后，选取均方差最小时的系数作为调整系数对角度比值进行最后的修正，再生成短程线网格。另外计算调整系数的时间在生成杆件的时间很小。可以从比较图中可以看出调整过后网格更加均匀。
同跨高网壳对比图
该短线程网壳的参数化建模用的弦分法，该划分方法得出网格质量比不上弧分法以及其他方法的。因此该参数化建模的优化方向是尝试用其他划分网格方法。 同时在了解工程实际后，进行短程程网格划分的实体并不是仅有球曲面，拥有其他异性构件，该参数化建模的拓展方向是建立对其他异性构件进行短线程网格划分的参数化方法。
附：代码与说明图
参数说明图1
参数说明图2
python
# -*- coding: utf-8 -*-
import math
import copy
import numpy as np
from abaqus import *
from abaqusConstants import *
from caeModules import *
from driverUtils import executeOnCaeStartup
"""简单坐标运算"""
def plus(x,y):
return (x[0] + y[0], x[1] + y[1], x[2] + y[2])
def cut(x,y):
return (x[0] - y[0], x[1] - y[1], x[2] - y[2])
def double(x,m):
return (x[0]*m, x[1]*m, x[2]*m)
def distance(x,y):
ai = x[0] - y[0]
bi = x[1] - y[1]
ci = x[2] - y[2]
d = math.pow((math.pow(ai,2)+math.pow(bi,2)+math.pow(ci,2)), 0.5)
return d
# 点映射到球面上的点坐标
def xyz_coordinates(x, radius):
if x[0] == 0:
if x[1] > 0:
xy_angle = float(math.pi) / 2
if x[1] < 0:
xy_angle = float(3 * math.pi) / 2
else:
xy_angle = math.atan(x[1] / x[0])
l_xy = math.pow((math.pow(x[0], 2) + math.pow(x[1], 2)), 0.5)
if l_xy == 0:
if x[2] >= 0:
z_angle = pi / 2
elif x[2] < 0:
z_angle = -1 * pi / 2
else:
z_angle = math.atan(x[2] / l_xy)
if x[0] >= 0:
xi = 1 * radius * math.cos(z_angle) * math.fabs(math.cos(xy_angle))
else:
xi = -1 * radius * math.cos(z_angle) * math.fabs(math.cos(xy_angle))
if x[1] >= 0:
yi = 1 * radius * math.cos(z_angle) * math.fabs(math.sin(xy_angle))
else:
yi = -1 * radius * math.cos(z_angle) * math.fabs(math.sin(xy_angle))
zi = radius * math.sin(z_angle)
return (xi, yi, zi)
"""粗略寻找较好的调整系数"""
#运用列举法以求杆长长度的方差最小时的调整系数
def search_best_adc(f1,f2,s,h):
Kn_number = 5 # 圆分数
pi = math.pi
angle = float(2 * pi) / Kn_number
radius = float(h) / 2 + (math.pow(s, 2) / 8) / h # 半径
z_angle = math.acos((radius - h) / radius)
coefficient = float(f1) / (f1 + f2) # 上下部分角度的分配系数
ad_list = [-0.02,-0.01,0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12]
t = radius*radius
for m in range(len(ad_list)):
adjustment_coefficient = ad_list[m]
Distances = []
last_coefficient = coefficient + adjustment_coefficient # 最终分配系数
z_angle_1 = last_coefficient * z_angle
"""上半部分坐标计算"""
coordinates0 = [ (math.sin(z_angle_1)*math.cos(angle*m),math.sin(z_angle_1)*math.sin(angle*m),math.cos(z_angle_1)) for m in range(Kn_number) ]
coordinates1 = [ ]
m = 0
for c in coordinates0:
a = double(cut(c , (0,0,1)), 1.0/f1)
coordinates = [ plus((0,0,1),double(a,m)) for m in range(f1+1) ]
coordinates1.append(coordinates)
m = m + 1
coordinates2 = [[(0,0,radius)],]
coordinates = [ ]
m=1
while m <= f1:
i=0
while i < Kn_number:
if i+1 < Kn_number :
a = double(cut(coordinates1[i][m],coordinates1[i+1][m]),1.0/m)
else:
a = double(cut(coordinates1[i][m],coordinates1[0][m]),1.0/m)
for j in range(m):
coordinates.append(cut(coordinates1[i][m], double(a,j)))
i = i + 1
coordinates2.append(coordinates)
coordinates = []
m = m + 1
for m in range(f1+1):
if m != 0:
for q in range(Kn_number*m):
coordinates2[m][q] = xyz_coordinates(coordinates2[m][q],radius)
"""上半部分杆长的计算"""
#同一层杆件的生成
m = 1
while m <= f1:
for q in range(Kn_number*m):
if q+1 < Kn_number*m:
Distances.append(distance(coordinates2[m][q],coordinates2[m][q+1]))
else:
Distances.append(distance(coordinates2[m][q],coordinates2[m][0]))
m = m + 1
#相邻层的杆件生成
coordinates3 = copy.deepcopy(coordinates2) #列表的深层复制
m = f1
while m > 0 :
q = Kn_number*m -1
while q >= 0:
if q % m == 0:
n = q / m * (m-1)
Distances.append(distance(coordinates2[m][q],coordinates2[m-1][int(n)]))
del coordinates3[m][q]
q = q-1
for q in range(Kn_number*(m-1)):
if q+1 < Kn_number*(m-1):
Distances.append(distance(coordinates3[m][q],coordinates2[m-1][q+1]))
Distances.append(distance(coordinates3[m][q],coordinates2[m-1][q]))
else:
Distances.append(distance(coordinates3[m][q],coordinates2[m-1][0]))
Distances.append(distance(coordinates3[m][q],coordinates2[m-1][q]))
m = m - 1
"""下半部分坐标计算"""
c1 = []
for i in range(5*f1):
if i % f1 ==0:
c1.append(double(coordinates2[f1][i], 1.0/radius))
#下界限坐标
c2 = []
for i in range(5):
c2.append( (math.sin(z_angle)*math.cos(float(f2)/f1*pi/5+i*2*pi/5), math.sin(z_angle)*math.sin(float(f2)/f1*pi/5+i*2*pi/5), math.cos(z_angle)) )
# 下界限辅助计算坐标
c21 = []
for i in range(5):
c21.append( (math.sin(z_angle)*math.cos(float(f2)/f1*pi/5+(float(f1-f2))/f1*2*pi/5+2*i*pi/5), math.sin(z_angle)*math.sin(float(f2)/f1*pi/5+(float(f1-f2))/f1*2*pi/5+2*i*pi/5), math.cos(z_angle)) )
#计算主要部分的节点坐标
m = 0
n = 0
c3 = []
while n < 5:
a = double(cut(c2[n],c1[m]), 1.0 / f2)
c0 = [plus(c1[m], double(a, j)) for j in range(f2 + 1)]
c3.append(c0)
if abs(m-n)>2:
break
if m > n :
n = n + 1
continue
if m == n:
m = m + 1
c2 [n] = c21[n]
if m == 5:
m = 0
#计算下半部分所用网点的坐标
c4 = []
c00 = []
for i in range(f2 + 1):
for m in range(10):
if m % 2 ==0 and i != f1:
a = double(cut(c3[m + 1][i], c3[m][i]), 1.0 / (f1 - i))
for v in range(int(f1) - i):
c00.append(plus(c3[m][i], double(a, v)))
elif m %2 != 0 and i!=0 :
if m < 9:
a = double(cut(c3[m + 1][i], c3[m][i]), 1.0 / i)
else:
a = double(cut(c3[0][i], c3[m][i]), 1.0 / i)
for v in range(i):
c00.append(plus(c3[m][i], double(a, v)))
c4.append(c00)
c00 = []
#点投影到球面上
for i in range(f2 +1):
for m in range(5*f1):
c4[i][m] = xyz_coordinates(c4[i][m],radius)
#!!!!!!上半部分和下半部分相连处的节点统一，使用上半部分的低层坐标
for i in range(5*f1):
c4[0][i] = coordinates2[f1][i]
"""下半部分的杆长计算"""
#同一层杆件的生成
m = 1
while m <= f2:
for q in range(5* f1):
if q + 1 < 5 * f1:
Distances.append(distance(c4[m][q], c4[m][q + 1]))
else:
Distances.append(distance(c4[m][q], c4[m][0]))
m = m+1
#不同层杆件的生成
m = 0
while m < f2:
for q in range(5 * f1):
if q > 0:
Distances.append(distance(c4[m][q], c4[m + 1][q]))
Distances.append(distance(c4[m][q], c4[m + 1][q-1]))
else:
Distances.append(distance(c4[m][q], c4[m + 1][q]))
Distances.append(distance(c4[m][q], c4[m + 1][-1]))
m = m + 1
"""方差的计算"""
distance_var = np.var(Distances)
#不同调整系数方差进行比较，取方差最小的时的调整系数
if distance_var <= t :
best_adc = adjustment_coefficient
t = distance_var
return best_adc
class Project_dcx():
def __init__(self, part_name, f1, f2 , s, h):
"""
参数设置
part_name = 'duanchengxian'
f1 = 5 #上半部分层数
f2 = 3 #下半部分层数，f2 <= f1/2
s = 200 #跨度
h = 120 #矢高
Kn_number = 5 # 圆分数
angle = float(2*pi)/Kn_number
radius = float(h)/2 + (math.pow(s,2)/8)/h #半径
z_angle = math.acos((radius-h)/radius)
coefficient = float(f1)/(f1+f2) #上下部分角度的分配系数
adjustment_coefficient # 修正系数，用于调整整体杆长方差
last_coefficient = coefficient + adjustment_coefficient #最终分配系数
z_angle_1 = coefficient*z_angle
"""
Kn_number = 5
angle = float(2*pi)/Kn_number
radius = float(h)/2 + (math.pow(s,2)/8)/h
z_angle = math.acos((radius-h)/radius)
coefficient = float(f1)/(f1+f2)
adjustment_coefficient = search_best_adc(f1,f2,s,h)
last_coefficient = coefficient + adjustment_coefficient
z_angle_1 = last_coefficient * z_angle
mdb.models['Model-1'].Part(name=part_name,dimensionality=THREE_D,type=DEFORMABLE_BODY)
p = mdb.models['Model-1'].parts[part_name]
self.f1 = f1
self.f2 = f2
self.p = p
self.radius = radius
self.z_angle_1 = z_angle_1
self.z_angle = z_angle
self.angle = angle
self.Kn_number = Kn_number
def part(self):
pi = math.pi
"""上半部分坐标计算"""
#底层坐标计算
coordinates0 = [ (math.sin(self.z_angle_1)*math.cos(self.angle*m),math.sin(self.z_angle_1)*math.sin(self.angle*m),cos(self.z_angle_1)) for m in range(self.Kn_number) ]
coordinates1 = [ ]
m = 0
#计算棱边上坐标
for c in coordinates0:
a = double(cut(c , (0,0,1)), 1.0/self.f1)
coordinates = [ plus((0,0,1),double(a,m)) for m in range(self.f1+1) ]
coordinates1.append(coordinates)
m = m + 1
coordinates2 = [[(0,0,self.radius)],] #手动补上顶点坐标
coordinates = [ ]
m=1
#计算各层坐标
while m <= self.f1:
i=0
while i < self.Kn_number:
if i+1 < self.Kn_number :
a = double(cut(coordinates1[i][m],coordinates1[i+1][m]),1.0/m)
else:
a = double(cut(coordinates1[i][m],coordinates1[0][m]),1.0/m)
for j in range(m):
coordinates.append(cut(coordinates1[i][m], double(a,j)))
i = i + 1
coordinates2.append(coordinates)
coordinates = []
m = m + 1
#直角坐标转换为球面坐标
for m in range(self.f1+1):
if m != 0:
for q in range(self.Kn_number*m):
coordinates2[m][q] = xyz_coordinates(coordinates2[m][q],self.radius)
"""上半部分单元生成"""
#同一层点的连接
m = 1
while m <= self.f1:
for q in range(self.Kn_number*m):
if q+1 < self.Kn_number*m:
self.p.WirePolyLine(points=(coordinates2[m][q],coordinates2[m][q+1]),meshable=ON)
else:
self.p.WirePolyLine(points=(coordinates2[m][q],coordinates2[m][0]),meshable=ON)
m = m + 1
#相邻层点的连接
coordinates3 = copy.deepcopy(coordinates2)
m = self.f1
while m > 0 :
q = self.Kn_number*m -1
while q >= 0:
if q % m == 0:
self.p.WirePolyLine(points=(coordinates2[m][q],coordinates2[m-1][q/(m)*(m-1)]),meshable=ON)
del coordinates3[m][q]
q = q-1
for q in range(self.Kn_number*(m-1)):
if q+1 < self.Kn_number*(m-1):
self.p.WirePolyLine(points=(coordinates3[m][q],coordinates2[m-1][q+1]),meshable=ON)