60
0
0
根据老学长的建模py脚本写了个机翼参数自动优化脚本
文章摘要
|
Chapter 1：从建模脚本开始
感谢老学长的贡献，感觉最难的一部分就是建模脚本化了（（（，赞美老学长.jpg
原始建模脚本：autowing2.2.py
介里系完整代码
from abaqus import *
from odbMaterial import *
from odbSection import *
from odbAccess import *
from abaqusConstants import *
import visualization
import mesh
from math import *
from struct import *
from string import *
import os
import os.path
import glob
import time
import part
import regionToolset
from abaqusConstants import *
import interaction
from caeModules import *
from te3 import *
from te2 import *
if org==1:
from te1 import *
else:
from te import *
from te4 import *
############ Parameters ######################
##part##
########### Materials ################
##LY-4##
DLY4 = 2850
ELY4 = (68000000000.0, 0.33)
##LY-9##
DLY9 = 2800
ELY9 = (68000000000.0,0.33)
##LY-12##
DLY12 = 2800
ELY12 = (68000000000.0,0.33)
##OCr15Ni7Mo2Al##
DOCr15 = 7800
EOCr15 = (199000000000,0.33)
############### Sketch ##########
dataName = str(profile3)
#dataName = path3+"wingprofile/"+profile3
print (dataName)
f = open(dataName, 'r')
data = f.readlines()
f.close()
aa0 = []
aa1 = []
for i in range(0, len(data)):
aa0.append(data[i].split()[0])
aa1.append(data[i].split()[1])
session.viewports['Viewport: 1'].setValues(displayedObject=None)
session.viewports['Viewport: 1'].setValues(displayedObject=None)
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__',
sheetSize=1000.0)
g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
s.setPrimaryObject(option=STANDALONE)
for i in range(0, len(data)):
s.Spot(point=(float(aa0[i]), float(aa1[i])))
bb = []
for i in range(0, len(data)):
bb.append((float(aa0[i]), float(aa1[i])))
bb.append((float(aa0[0]), float(aa1[0])))
s.Spline(points=(bb))
mdb.models['Model-1'].sketches.changeKey(fromName='__profile__',toName='Sketch-1')
##shuofang##xianchang###
object=[]
object.append(g[2])
for i in range(0,len(v)):
object.append(v[i])
s.scale(scaleValue=l, scaleCenter=(0.0, 0.0), objectList=(object))
s.unsetPrimaryObject()
###### Step Time ##########
st = 0.004
IT_target = 5.0e-8
##############################################
############ Materials ################
##
mdb.models['Model-1'].Material(name='LY-4')
mdb.models['Model-1'].materials['LY-4'].Density(table=((DLY4,),))
mdb.models['Model-1'].materials['LY-4'].Elastic(table=(ELY4,))
##
mdb.models['Model-1'].Material(name='LY-9')
mdb.models['Model-1'].materials['LY-9'].Density(table=((DLY9,),))
mdb.models['Model-1'].materials['LY-9'].Elastic(table=(ELY9,))
##
mdb.models['Model-1'].Material(name='LY-12')
mdb.models['Model-1'].materials['LY-12'].Density(table=((DLY12,),))
mdb.models['Model-1'].materials['LY-12'].Elastic(table=(ELY12,))
##
mdb.models['Model-1'].Material(name='OCr15Ni7Mo2Al')
mdb.models['Model-1'].materials['OCr15Ni7Mo2Al'].Density(table=((DOCr15,),))
mdb.models['Model-1'].materials['OCr15Ni7Mo2Al'].Elastic(table=(EOCr15,))
#################part#########
session.viewports['Viewport: 1'].setValues(displayedObject=None)
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__',
sheetSize=200.0)
g, v1, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
s.setPrimaryObject(option=STANDALONE)
s.sketchOptions.setValues(gridOrigin=(0.4997, 0.0154334306716919)) ###***
s.retrieveSketch(sketch=mdb.models['Model-1'].sketches['Sketch-1'])
session.viewports['Viewport: 1'].view.fitView()
p = mdb.models['Model-1'].Part(name='plate-1', dimensionality=THREE_D,
type=DEFORMABLE_BODY)
p = mdb.models['Model-1'].parts['plate-1']
p.BaseShellExtrude(sketch=s, depth=W)
s.unsetPrimaryObject()
p = mdb.models['Model-1'].parts['plate-1']
session.viewports['Viewport: 1'].setValues(displayedObject=p)
del mdb.models['Model-1'].sketches['__profile__']
c=p.faces
mdb.models['Model-1'].parts['plate-1'].Set(faces=c[0:len(c)], name='mengpi' )
################# partition
sum1 = 0.0
p = mdb.models['Model-1'].parts['plate-1']
for i in range(0, A1-1):
sum1 = sum1 + D1
p.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=sum1)
p = mdb.models['Model-1'].parts['plate-1']
p.DatumPlaneByPrincipalPlane(principalPlane=YZPLANE, offset=l1*l)
A2 = int((1.0 - l1) * l // D2)-1
A3=int((l1*l//D2))
sum2 = 0.0
sum3 = 0.0
d=p.datums
for i in range(0, A2-1):
sum2=sum2+D2
p.DatumPlaneByOffset(plane=d[2+A1],flip=SIDE1,offset=sum2)
for i in range(0, A3-1):
sum3=sum3+D2
p.DatumPlaneByOffset(plane=d[2+A1],flip=SIDE2,offset=sum3)
for i in range(0, A1+A2+A3-2):
p = mdb.models['Model-1'].parts['plate-1']
c, e, v, d = p.faces, p.edges, p.vertices, p.datums
p.PartitionFaceByDatumPlane(datumPlane=d[i+3], faces=c[0:len(c)])
###coveredges###leiban#fuban####
####leiban####
for i in range(0,A1-1):
str1=str('Partition '+'face'+str(-(i+1)))
p.CoverEdges(edgeList = p.getFeatureEdges(str1),tryAnalytical=True)
p = mdb.models['Model-1'].parts['plate-1']
e = p.edges
f=p.faces
e1=e.getByBoundingBox(-50.0, -50.0, -0.1, 50.0, 50.0, 0.1)
e2=e.getByBoundingBox(-50.0, -50.0, W-0.1, 50.0, 50.0, W+0.1)
p.CoverEdges(edgeList = e1+e2,tryAnalytical=True)
p.PartitionFaceByDatumPlane(datumPlane=d[A1+2], faces=c[0:len(c)])
p.PartitionFaceByDatumPlane(datumPlane=d[int(A1+2+0.6*A2)], faces=c[0:len(c)])
#setPTleiban#
fff=[]
for i in range(0,1):
ff=f.getByBoundingBox(-50.0, -50.0, i*D1-0.01, 50.0, 50.0, i*D1+0.01)
fff.append(ff)
for i in range(2,int(0.5*A1)):
ff = f.getByBoundingBox(-50.0, -50.0, i * D1 - 0.001, 50.0, 50.0, i * D1 + 0.001)
fff.append(ff)
for i in range(int(0.5*A1+1),A1+1):
ff = f.getByBoundingBox(-50.0, -50.0, i * D1 - 0.001, 50.0, 50.0, i * D1 + 0.001)
fff.append(ff)
mdb.models['Model-1'].parts['plate-1'].Set(faces=fff, name='PTleiban' )
#setJQleiban#
fff=[]
ff1=f.getByBoundingBox(-50.0, -50.0, 1*D1-0.0001, 50.0, 50.0, 1*D1+0.0001)
ff2=f.getByBoundingBox(-50.0, -50.0, int(0.5*A1)*D1-0.0001, 50.0, 50.0, int(0.5*A1)*D1+0.0001)
fff.append(ff1)
fff.append(ff2)
mdb.models['Model-1'].parts['plate-1'].Set(faces=fff, name='JQleiban' )
####fuban####
str2=str('Partition '+'face'+str(-(A1)))
e3=e.findAt(((l*l1,0,0),))
e4=e.findAt(((l*l1,0,W),))
p.CoverEdges(edgeList = p.getFeatureEdges(str2)+e3+e4,tryAnalytical=True)
str4=str('Cover '+ 'edges'+str(-(A1+1)))
mdb.models['Model-1'].parts['plate-1'].Set(faces=p.getFeatureFaces(str4), name='qianliangfuban')
e1=e.getByBoundingBox(l*l1+int(0.6*A2)*D2-0.001,0,-50,l*l1+int(0.6*A2)*D2+0.001,50,50)
e2=e.getByBoundingBox(l*l1+int(0.6*A2)*D2-0.001,-50,-50,l*l1+int(0.6*A2)*D2+0.001,0,50)
e3=e.findAt(((l*l1+int(0.6*A2)*D2,0,0),))
e4=e.findAt(((l*l1+int(0.6*A2)*D2,0,W),))
p.CoverEdges(edgeList = e1+e2+e3+e4,tryAnalytical=True)
str4=str('Cover '+ 'edges'+str(-(A1+2)))
mdb.models['Model-1'].parts['plate-1'].Set(faces=p.getFeatureFaces(str4), name='houliangfuban')
########Stringer#########
##changhang##
p = mdb.models['Model-1'].parts['plate-1']
e=p.edges
eee=[]
for i in range(0,int(A2*0.6)-1):
ee1=e.getByBoundingBox(l*l1+(i+1)*D2-0.001, -50.0, -50.0, l*l1+(i+1)*D2+0.001, 50.0, 50.0)
eee.append(ee1)
for i in range(int(A2*0.6),A2):
ee1 = e.getByBoundingBox(l*l1+(i+1)*D2-0.001,-50.0,-50.0,l*l1+(i+1)*D2+0.001,50.0,50.0)
eee.append(ee1)
p.Stringer(edges=eee,name='STchanghang1')
mdb.models['Model-1'].parts['plate-1'].Set(stringerEdges=(('STchanghang1',eee),), name='changhang1')
eee=[]
for i in range(0,A3):
ee2=e.getByBoundingBox(l*l1-(i+1)*D2-0.001, -50.0, -50.0, l*l1-(i+1)*D2+0.001, 50.0, 50.0)
eee.append(ee2)
p.Stringer(edges=eee,name='STchanghang2')
mdb.models['Model-1'].parts['plate-1'].Set(stringerEdges=(('STchanghang2',eee),), name='changhang2')
##liang##
p = mdb.models['Model-1'].parts['plate-1']
e=p.edges
eee=[]
ee1=e.getByBoundingBox(l*l1-0.001,0,-50,l*l1+0.001,50,50)
ee2=e.getByBoundingBox(l*l1-0.001,-50,-50,l*l1+0.001,0,50)
eee.append(ee1)
eee.append(ee2)
p.Stringer(edges=eee,name='STqianliang')
mdb.models['Model-1'].parts['plate-1'].Set(stringerEdges=(('STqianliang',eee),), name='qianliang')
##
eee=[]
ee1=e.getByBoundingBox(l*l1+int(0.6*A2)*D2-0.001,0,-50,l*l1+int(0.6*A2)*D2+0.001,50,50)
ee2=e.getByBoundingBox(l*l1+int(0.6*A2)*D2-0.001,-50,-50,l*l1+int(0.6*A2)*D2+0.001,0,50)
eee.append(ee1)
eee.append(ee2)
p.Stringer(edges=eee,name='SThouliang')
mdb.models['Model-1'].parts['plate-1'].Set(stringerEdges=(('SThouliang',eee),), name='houliang')
##JQleibanyuantiao##
p = mdb.models['Model-1'].parts['plate-1']
e=p.edges
eee=[]
ee1=e.getByBoundingBox(-50.0, 0, 1*D1-0.00001, 50.0, 50.0, 1*D1+0.00001)
ee2=e.getByBoundingBox(-50.0, -50.0, 1*D1-0.00001, 50.0, 0, 1*D1+0.00001)
ee3=e.getByBoundingBox(l1*l+0.0001, -50.0, 1*D1-0.00001, 50.0, 50.0, 1*D1+0.00001)
ee4=e.getByBoundingBox(-50.0, -50.0, 1*D1-0.00001, l1*l-0.0001, 50.0, 1*D1+0.00001)
eee.append(ee1)
eee.append(ee2)
eee.append(ee3)
eee.append(ee4)
p.Stringer(edges=eee,name='jiangqianglei1')
mdb.models['Model-1'].parts['plate-1'].Set(stringerEdges=(('jiangqianglei1',eee),), name='jiangqianglei1')
eee=[]
ee1=e.getByBoundingBox(-50.0, 0, int(0.5*A1)*D1-0.00001, 50.0, 50.0, int(0.5*A1)*D1+0.00001)
ee2=e.getByBoundingBox(-50.0, -50.0, int(0.5*A1)*D1-0.00001, 50.0, 0, int(0.5*A1)*D1+0.00001)
ee3=e.getByBoundingBox(l1*l+0.0001, -50.0, int(0.5*A1)*D1-0.00001, 50.0, 50.0, int(0.5*A1)*D1+0.00001)
ee4=e.getByBoundingBox(-50.0, -50.0, int(0.5*A1)*D1-0.00001, l1*l-0.0001, 50.0, int(0.5*A1)*D1+0.00001)
eee.append(ee1)
eee.append(ee2)
eee.append(ee3)
eee.append(ee4)
p.Stringer(edges=eee,name='jiangqianglei2')
mdb.models['Model-1'].parts['plate-1'].Set(stringerEdges=(('jiangqianglei2',eee),), name='jiangqianglei2')
#######poumian#######
if poumian==1:
mdb.models['Model-1'].LProfile(name='changhang', a=numA2, b=numB2, t1=numt1, t2=numt2)
mdb.models['Model-1'].RectangularProfile(name='liang', a=numA1, b=numB1)
mdb.models['Model-1'].RectangularProfile(name='leibanyuantiao', a=0.5*numA1, b=0.5*numB1)
elif poumian==2:
mdb.models['Model-1'].IProfile(name='liang', l=numL1, h=numH1, b1=numB1,b2=numB2, t1=numT1, t2=numT2, t3=numT3)
mdb.models['Model-1'].TProfile(name='changhang', b=numB3, h=numH2, l=numL2, tf=numTF, tw=numTW)
mdb.models['Model-1'].RectangularProfile(name='leibanyuantiao', a=0.5 * numB1, b=0.5 * numB1)
elif poumian==3:
mdb.models['Model-1'].TProfile(name='liang', b=numB4, h=numH3, l=numL3, tf=numTF2, tw=numTW2)
mdb.models['Model-1'].LProfile(name='changhang', a=numA2, b=numB2, t1=numt1, t2=numt2)
mdb.models['Model-1'].RectangularProfile(name='leibanyuantiao', a=0.5 * numB4, b=0.25 * numB4)
else:None
#######jiemian#######
###
mdb.models['Model-1'].BeamSection(name='liang', integration=DURING_ANALYSIS,
poissonRatio=0.0, profile='liang', material='OCr15Ni7Mo2Al',
temperatureVar=LINEAR, consistentMassMatrix=False)
mdb.models['Model-1'].HomogeneousShellSection(name='mp', preIntegrate=OFF,
material='LY-12', thicknessType=UNIFORM, thickness=0.002,
thicknessField='', idealization=NO_IDEALIZATION, poissonDefinition=DEFAULT,
thicknessModulus=None, temperature=GRADIENT, useDensity=OFF,
integrationRule=SIMPSON, numIntPts=5)
mdb.models['Model-1'].HomogeneousShellSection(name='liangfuban',
preIntegrate=OFF, material='LY-9', thicknessType=UNIFORM, thickness=numh1,
thicknessField='', idealization=NO_IDEALIZATION, poissonDefinition=DEFAULT,
thicknessModulus=None, temperature=GRADIENT, useDensity=OFF,
integrationRule=SIMPSON, numIntPts=5)
mdb.models['Model-1'].BeamSection(name='changhang',
integration=DURING_ANALYSIS, poissonRatio=0.0, profile='changhang',
material='LY-4', temperatureVar=LINEAR, consistentMassMatrix=False)
mdb.models['Model-1'].BeamSection(name='leibanyuantiao',
integration=DURING_ANALYSIS, poissonRatio=0.0, profile='leibanyuantiao',
material='LY-4', temperatureVar=LINEAR, consistentMassMatrix=False)
mdb.models['Model-1'].HomogeneousShellSection(name='PTleiban', preIntegrate=OFF,
material='LY-4', thicknessType=UNIFORM, thickness=numh2, thicknessField='',
idealization=NO_IDEALIZATION, poissonDefinition=DEFAULT,
thicknessModulus=None, temperature=GRADIENT, useDensity=OFF,
integrationRule=SIMPSON, numIntPts=5)
mdb.models['Model-1'].HomogeneousShellSection(name='JQleiban', preIntegrate=OFF,
material='LY-4', thicknessType=UNIFORM, thickness=2*numh2, thicknessField='',
idealization=NO_IDEALIZATION, poissonDefinition=DEFAULT,
thicknessModulus=None, temperature=GRADIENT, useDensity=OFF,
integrationRule=SIMPSON, numIntPts=5)
###########zhipai###
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['qianliangfuban']
p = mdb.models['Model-1'].parts['plate-1']
p.SectionAssignment(region=region, sectionName='liangfuban', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['houliangfuban']
p = mdb.models['Model-1'].parts['plate-1']
p.SectionAssignment(region=region, sectionName='liangfuban', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['PTleiban']
p = mdb.models['Model-1'].parts['plate-1']
p.SectionAssignment(region=region, sectionName='PTleiban', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['JQleiban']
p = mdb.models['Model-1'].parts['plate-1']
p.SectionAssignment(region=region, sectionName='JQleiban', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['mengpi']
p = mdb.models['Model-1'].parts['plate-1']
p.SectionAssignment(region=region, sectionName='mp', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['qianliang']
p = mdb.models['Model-1'].parts['plate-1']
p.SectionAssignment(region=region, sectionName='liang', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
if num==2:
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['houliang']
p = mdb.models['Model-1'].parts['plate-1']
p.SectionAssignment(region=region, sectionName='liang', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
elif num==1:
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['houliang']
p = mdb.models['Model-1'].parts['plate-1']
p.SectionAssignment(region=region, sectionName='changhang', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
else:
None
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['changhang1']
p = mdb.models['Model-1'].parts['plate-1']
p.SectionAssignment(region=region, sectionName='changhang', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['changhang2']
p = mdb.models['Model-1'].parts['plate-1']
p.SectionAssignment(region=region, sectionName='changhang', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['jiangqianglei1']
p = mdb.models['Model-1'].parts['plate-1']
p.SectionAssignment(region=region, sectionName='leibanyuantiao', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['jiangqianglei2']
p = mdb.models['Model-1'].parts['plate-1']
p.SectionAssignment(region=region, sectionName='leibanyuantiao', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
###ST fangxiang##
p = mdb.models['Model-1'].parts['plate-1']
region=p.sets['qianliang']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(1.0, 0.0,0.0))
p = mdb.models['Model-1'].parts['plate-1']
region=p.sets['houliang']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(1.0, 0.0,0.0))
if poumian == 1:
p = mdb.models['Model-1'].parts['plate-1']
region=p.sets['changhang1']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(-1.0, 0.0, 0.0))
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['houliang']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(-1.0, 0.0, 0.0))
p = mdb.models['Model-1'].parts['plate-1']
region=p.sets['changhang2']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(1.0, 0.0, 0.0))
elif poumian==2:
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['changhang1']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(1.0, 0.0, 0.0))
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['houliang']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(1.0, 0.0, 0.0))
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['changhang2']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(-1.0, 0.0, 0.0))
elif poumian == 3:
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['changhang1']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(-1.0, 0.0, 0.0))
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['houliang']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(-1.0, 0.0, 0.0))
p = mdb.models['Model-1'].parts['plate-1']
region = p.sets['changhang2']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(1.0, 0.0, 0.0))
else:None
p = mdb.models['Model-1'].parts['plate-1']
region=p.sets['jiangqianglei1']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(0.0, 0.0, -1.0))
p = mdb.models['Model-1'].parts['plate-1']
region=p.sets['jiangqianglei2']
p = mdb.models['Model-1'].parts['plate-1']
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(0.0, 0.0, -1.0))
############## Assembly #######################
a = mdb.models['Model-1'].rootAssembly
session.viewports['Viewport: 1'].setValues(displayedObject=a)
session.viewports['Viewport: 1'].assemblyDisplay.setValues(
optimizationTasks=OFF, geometricRestrictions=OFF, stopConditions=OFF)
a = mdb.models['Model-1'].rootAssembly
a.DatumCsysByDefault(CARTESIAN)
p = mdb.models['Model-1'].parts['plate-1']
a.Instance(name='plate-1-1', part=p, dependent=OFF)
############step#################
mdb.models['Model-1'].StaticStep(name='Step-1', previous='Initial',
stabilizationMagnitude=0.0002,
stabilizationMethod=DISSIPATED_ENERGY_FRACTION,
continueDampingFactors=False, adaptiveDampingRatio=0.05, initialInc=0.1,
minInc=1e-09, nlgeom=ON)
session.viewports['Viewport: 1'].assemblyDisplay.setValues(step='Step-1')
###########interactive#########
a = mdb.models['Model-1'].rootAssembly
session.viewports['Viewport: 1'].setValues(displayedObject=a)
session.viewports['Viewport: 1'].assemblyDisplay.setValues(interactions=ON, constraints=ON, connectors=ON, engineeringFeatures=ON)
a = mdb.models['Model-1'].rootAssembly
if num==2:
RP1=a.ReferencePoint(point=(0.5*l, 0.0, 1*D1))
RP2=a.ReferencePoint(point=(0.5*l, 0.0, int(0.5*A1)*D1))
RP3=a.ReferencePoint(point=(0.5*l, 0.0, int(0.3*A1)*D1))
elif num==1:
RP1 = a.ReferencePoint(point=(0.5 * l, 0.0, 1 * D1))
RP2 = a.ReferencePoint(point=(0.5 * l, 0.0, int(0.5 * A1) * D1))
RP3 = a.ReferencePoint(point=(0.5 * l, 0.0, int(0.3 * A1) * D1))
RP4=a.ReferencePoint(point=(l*l1+int(0.6*A2)*D2,0,0))
else:
None
a = mdb.models['Model-1'].rootAssembly
r1 = a.referencePoints
refPoints1=(r1[4], )
region1=regionToolset.Region(referencePoints=refPoints1)
a = mdb.models['Model-1'].rootAssembly
s1 = a.instances['plate-1-1'].faces
side1Faces1 = s1.getByBoundingBox(-50.0, -50.0, 1*D1-0.001, 50.0, 50.0, 1*D1+0.001)
region2=regionToolset.Region(side1Faces=side1Faces1)
mdb.models['Model-1'].Coupling(name='zaidan1', controlPoint=region1,
surface=region2, influenceRadius=WHOLE_SURFACE, couplingType=DISTRIBUTING,
weightingMethod=UNIFORM, localCsys=None, u1=ON, u2=ON, u3=ON, ur1=ON,
ur2=ON, ur3=ON)
region=regionToolset.Region(referencePoints=refPoints1)
mdb.models['Model-1'].rootAssembly.engineeringFeatures.PointMassInertia(
name='zaidan1', region=region, mass=mass*0.5*0.05, alpha=0.0, composite=0.0)
refPoints1=(r1[5], )
region1=regionToolset.Region(referencePoints=refPoints1)
a = mdb.models['Model-1'].rootAssembly
s1 = a.instances['plate-1-1'].faces
side1Faces1 = s1.getByBoundingBox(-50.0, -50.0, int(0.5*A1)*D1-0.001, 50.0, 50.0, int(0.5*A1)*D1+0.001)
region2=regionToolset.Region(side1Faces=side1Faces1)
mdb.models['Model-1'].Coupling(name='zaidan2', controlPoint=region1,
surface=region2, influenceRadius=WHOLE_SURFACE, couplingType=DISTRIBUTING,
weightingMethod=UNIFORM, localCsys=None, u1=ON, u2=ON, u3=ON, ur1=ON,
ur2=ON, ur3=ON)
region=regionToolset.Region(referencePoints=refPoints1)
mdb.models['Model-1'].rootAssembly.engineeringFeatures.PointMassInertia(
name='zaidan2', region=region, mass=mass*0.5*0.05, alpha=0.0, composite=0.0)
refPoints1=(r1[6], )
region1=regionToolset.Region(referencePoints=refPoints1)
a = mdb.models['Model-1'].rootAssembly
s1 = a.instances['plate-1-1'].faces
side1Faces1 = s1.getByBoundingBox(-50.0, -50.0, 1 * D1 - 0.001, 50.0, 50.0, 1 * D1 + 0.001)
side2Faces2 = s1.getByBoundingBox(-50.0, -50.0, int(0.5 * A1) * D1 - 0.001, 50.0, 50.0, int(0.5 * A1) * D1 + 0.001)
region2=regionToolset.Region(side1Faces=side1Faces1+side2Faces2)
mdb.models['Model-1'].Coupling(name='youxiang', controlPoint=region1,
surface=region2, influenceRadius=WHOLE_SURFACE, couplingType=DISTRIBUTING,
weightingMethod=UNIFORM, localCsys=None, u1=ON, u2=ON, u3=ON, ur1=ON,
ur2=ON, ur3=ON)
region=regionToolset.Region(referencePoints=refPoints1)
mdb.models['Model-1'].rootAssembly.engineeringFeatures.PointMassInertia(
name='youxiang', region=region, mass=mass*0.5*0.35, alpha=0.0, composite=0.0)
###########load###############
##BC##shuangliang
if num==2:
a = mdb.models['Model-1'].rootAssembly
e = a.instances['plate-1-1'].edges
e1 = e.findAt(((l*l1,0,0),))
e2 = e.findAt(((l*l1+int(0.6*A2)*D2,0,0),))
a.Set(edges=e1+e2, name='BC-1')
region = a.sets['BC-1']
mdb.models['Model-1'].EncastreBC(name='BC-1', createStepName='Initial',
region=region, localCsys=None)
##BC#####danliang
elif num==1:
a = mdb.models['Model-1'].rootAssembly
e = a.instances['plate-1-1'].edges
n = a.instances['plate-1-1'].nodes
e1 = e.findAt(((l*l1,0,0),))
e2 = e.findAt(((l*l1+int(0.6*A2)*D2,0,0),))
a.Set(edges=e1, name='BC-1')
a.Set(edges=e2, name='BC-12')
region = a.sets['BC-1']
mdb.models['Model-1'].EncastreBC(name='BC-1', createStepName='Initial',
region=region, localCsys=None)
a = mdb.models['Model-1'].rootAssembly
r1 = a.referencePoints
refPoints1=(r1[7], )
region = regionToolset.Region(referencePoints=refPoints1)
mdb.models['Model-1'].PinnedBC(name='BC-2', createStepName='Initial',
region=region, localCsys=None)
a = mdb.models['Model-1'].rootAssembly
r1 = a.referencePoints
refPoints1=(r1[7], )
region1=regionToolset.Region(referencePoints=refPoints1)
a = mdb.models['Model-1'].rootAssembly
region2 = a.sets['BC-12']
mdb.models['Model-1'].Coupling(name='houliangjiaojie', controlPoint=region1,
surface=region2, influenceRadius=WHOLE_SURFACE, couplingType=DISTRIBUTING,
weightingMethod=UNIFORM, localCsys=None, u1=ON, u2=ON, u3=ON, ur1=ON,
ur2=ON, ur3=ON)
else:
None
##load
mdb.models['Model-1'].Gravity(name='gravity', createStepName='Step-1',
comp2=9.8*overload, distributionType=UNIFORM, field='')
#############mesh###########
a = mdb.models['Model-1'].rootAssembly
partInstances =(a.instances['plate-1-1'], )
a.seedPartInstance(regions=partInstances, size=SIZE, deviationFactor=0.1,
minSizeFactor=0.1)
a = mdb.models['Model-1'].rootAssembly
partInstances =(a.instances['plate-1-1'], )
a.generateMesh(regions=partInstances)
##########job##########
mdb.Job(name='Job-wing', model='Model-1', description='', type=ANALYSIS,
atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90,
memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True,
explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF,
modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='',
scratch='', resultsFormat=ODB, multiprocessingMode=DEFAULT, numCpus=2,
numDomains=2, numGPUs=1)
mdb.jobs['Job-wing'].submit(consistencyChecking=OFF)
##########danliang########
Chapter 2：参数化求解与自动优化
一顿改改改，报错一顿修修修，电脑逐渐红温，眼皮逐渐打架，精神逐渐涣散......
参数化求解器：wing_solver.py
点击查看秘密情报
# -*- coding: utf-8 -*-
from abaqus import *
from abaqusConstants import *
from odbMaterial import *
from odbSection import *
from odbAccess import *
import regionToolset
from math import *
def run_wing_simulation(W=5.0, overload=2.0, aircraft_mass=6000.0,
D1=0.5, D2=0.15,
mp_t=0.002, numh1=0.005, numh2=0.005,
st_a=0.03, st_b=0.03, st_t1=0.005, st_t2=0.005,
bm_a=0.05, bm_b=0.05,
rebuild_geo=True):
# 11变量参数化（实际9，间距离散性变量不兼容）机翼求解器，支持基于拓扑判断的属性免新建模模式，提升优化效率。
job_name = 'Job-wing'
try:
if rebuild_geo:
# 模式一：完整建模 (当 D1 或 D2 发生变化时，重新构建几何与网格)
Mdb()
profile3 = r"E:\EXP\Fifth\wingprofile\NACA_2412.dat" # 懒得写个文件路径，直接写死了，反正不改翼型（（（（x
l = 1.5 # 弦长比例
l1 = 0.3 # 前梁位置
num = 1 # 单梁
SIZE = 0.05 # 全局网格尺寸
# 拓扑自适应划分
target_D1, target_D2 = D1, D2
A1_bays = max(2, int(round(W / target_D1)))
D1_real = W / A1_bays
A1 = A1_bays
L_rear = (1.0 - l1) * l
A2_bays = max(3, int(round(L_rear / target_D2)))
D2_rear = L_rear / A2_bays
A2 = A2_bays
L_front = l1 * l
A3_bays = max(2, int(round(L_front / target_D2)))
D2_front = L_front / A3_bays
A3 = A3_bays
# 材料定义
DLY4, DLY9, DLY12, DOCr15 = 2850, 2800, 2800, 7800
ELY4, ELY9, ELY12, EOCr15 = (680e9, 0.33), (680e9,0.33), (680e9,0.33), (199e9,0.33)
mdb.models['Model-1'].Material(name='LY-4').Density(table=((DLY4,),))
mdb.models['Model-1'].materials['LY-4'].Elastic(table=(ELY4,))
mdb.models['Model-1'].Material(name='LY-9').Density(table=((DLY9,),))
mdb.models['Model-1'].materials['LY-9'].Elastic(table=(ELY9,))
mdb.models['Model-1'].Material(name='LY-12').Density(table=((DLY12,),))
mdb.models['Model-1'].materials['LY-12'].Elastic(table=(ELY12,))
mdb.models['Model-1'].Material(name='OCr15Ni7Mo2Al').Density(table=((DOCr15,),))
mdb.models['Model-1'].materials['OCr15Ni7Mo2Al'].Elastic(table=(EOCr15,))
# 翼型草图绘制
f = open(str(profile3), 'r')
data = f.readlines()
f.close()
aa0, aa1, bb = [], [], []
for i in range(0, len(data)):
aa0.append(data[i].split()[0])
aa1.append(data[i].split()[1])
bb.append((float(aa0[i]), float(aa1[i])))
bb.append((float(aa0[0]), float(aa1[0])))
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=1000.0)
g, v = s.geometry, s.vertices
s.setPrimaryObject(option=STANDALONE)
for i in range(0, len(data)): s.Spot(point=(float(aa0[i]), float(aa1[i])))
s.Spline(points=(bb))
mdb.models['Model-1'].sketches.changeKey(fromName='__profile__',toName='Sketch-1')
object=[g[2]]
for i in range(0,len(v)): object.append(v[i])
s.scale(scaleValue=l, scaleCenter=(0.0, 0.0), objectList=(object))
s.unsetPrimaryObject()
# 三维拉伸
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
s.setPrimaryObject(option=STANDALONE)
s.retrieveSketch(sketch=mdb.models['Model-1'].sketches['Sketch-1'])
p = mdb.models['Model-1'].Part(name='plate-1', dimensionality=THREE_D, type=DEFORMABLE_BODY)
p.BaseShellExtrude(sketch=s, depth=W)
s.unsetPrimaryObject()
del mdb.models['Model-1'].sketches['__profile__']
c=p.faces
mdb.models['Model-1'].parts['plate-1'].Set(faces=c[0:len(c)], name='mengpi' )
# 生成切分基准面
sum1 = 0.0
for i in range(0, A1-1):
sum1 = sum1 + D1_real
p.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=sum1)
p.DatumPlaneByPrincipalPlane(principalPlane=YZPLANE, offset=l1*l)
sum2, sum3 = 0.0, 0.0
d=p.datums
for i in range(0, A2-1):
sum2 = sum2 + D2_rear
p.DatumPlaneByOffset(plane=d[2+A1],flip=SIDE1,offset=sum2)
for i in range(0, A3-1):
sum3 = sum3 + D2_front
p.DatumPlaneByOffset(plane=d[2+A1],flip=SIDE2,offset=sum3)
# 执行切分与蒙皮盖面
for i in range(0, A1+A2+A3-2):
c = p.faces
p.PartitionFaceByDatumPlane(datumPlane=d[i+3], faces=c[0:len(c)])
for i in range(0,A1-1):
str1=str('Partition '+'face'+str(-(i+1)))
p.CoverEdges(edgeList = p.getFeatureEdges(str1),tryAnalytical=True)
e, f = p.edges, p.faces
e1, e2 = e.getByBoundingBox(-50, -50, -0.1, 50, 50, 0.1), e.getByBoundingBox(-50, -50, W-0.1, 50, 50, W+0.1)
p.CoverEdges(edgeList = e1+e2,tryAnalytical=True)
p.PartitionFaceByDatumPlane(datumPlane=d[A1+2], faces=c[0:len(c)])
p.PartitionFaceByDatumPlane(datumPlane=d[int(A1+2+0.6*A2)], faces=c[0:len(c)])
# 创建面集合 (Sets)
fff = [f.getByBoundingBox(-50, -50, -0.01, 50, 50, 0.01)]
for i in range(2,int(0.5*A1)): fff.append(f.getByBoundingBox(-50, -50, i * D1_real - 0.001, 50, 50, i * D1_real + 0.001))
for i in range(int(0.5*A1+1),A1+1): fff.append(f.getByBoundingBox(-50, -50, i * D1_real - 0.001, 50, 50, i * D1_real + 0.001))
mdb.models['Model-1'].parts['plate-1'].Set(faces=fff, name='PTleiban' )
fff = [f.getByBoundingBox(-50, -50, 1*D1_real-0.0001, 50, 50, 1*D1_real+0.0001),
f.getByBoundingBox(-50, -50, int(0.5*A1)*D1_real-0.0001, 50, 50, int(0.5*A1)*D1_real+0.0001)]
mdb.models['Model-1'].parts['plate-1'].Set(faces=fff, name='JQleiban' )
str2=str('Partition '+'face'+str(-(A1)))
p.CoverEdges(edgeList = p.getFeatureEdges(str2)+e.findAt(((l*l1,0,0),))+e.findAt(((l*l1,0,W),)),tryAnalytical=True)
mdb.models['Model-1'].parts['plate-1'].Set(faces=p.getFeatureFaces(str('Cover edges'+str(-(A1+1)))), name='qianliangfuban')
e1, e2 = e.getByBoundingBox(l*l1+int(0.6*A2)*D2_rear-0.001,0,-50,l*l1+int(0.6*A2)*D2_rear+0.001,50,50), e.getByBoundingBox(l*l1+int(0.6*A2)*D2_rear-0.001,-50,-50,l*l1+int(0.6*A2)*D2_rear+0.001,0,50)
p.CoverEdges(edgeList = e1+e2+e.findAt(((l*l1+int(0.6*A2)*D2_rear,0,0),))+e.findAt(((l*l1+int(0.6*A2)*D2_rear,0,W),)),tryAnalytical=True)
mdb.models['Model-1'].parts['plate-1'].Set(faces=p.getFeatureFaces(str('Cover edges'+str(-(A1+2)))), name='houliangfuban')
# 创建梁及长桁单元特征 (Stringers)
eee=[]
for i in range(0,int(A2*0.6)-1): eee.append(e.getByBoundingBox(l*l1+(i+1)*D2_rear-0.001, -50, -50, l*l1+(i+1)*D2_rear+0.001, 50, 50))
for i in range(int(A2*0.6),A2-1): eee.append(e.getByBoundingBox(l*l1+(i+1)*D2_rear-0.001,-50,-50,l*l1+(i+1)*D2_rear+0.001,50,50))
p.Stringer(edges=eee,name='STchanghang1')
mdb.models['Model-1'].parts['plate-1'].Set(stringerEdges=(('STchanghang1',eee),), name='changhang1')
eee=[]
for i in range(0,A3-1): eee.append(e.getByBoundingBox(l*l1-(i+1)*D2_front-0.001, -50, -50, l*l1-(i+1)*D2_front+0.001, 50, 50))
p.Stringer(edges=eee,name='STchanghang2')
mdb.models['Model-1'].parts['plate-1'].Set(stringerEdges=(('STchanghang2',eee),), name='changhang2')
p.Stringer(edges=[e.getByBoundingBox(l*l1-0.001,0,-50,l*l1+0.001,50,50), e.getByBoundingBox(l*l1-0.001,-50,-50,l*l1+0.001,0,50)],name='STqianliang')
mdb.models['Model-1'].parts['plate-1'].Set(stringerEdges=(('STqianliang',p.stringers['STqianliang'].edges),), name='qianliang')
p.Stringer(edges=[e.getByBoundingBox(l*l1+int(0.6*A2)*D2_rear-0.001,0,-50,l*l1+int(0.6*A2)*D2_rear+0.001,50,50), e.getByBoundingBox(l*l1+int(0.6*A2)*D2_rear-0.001,-50,-50,l*l1+int(0.6*A2)*D2_rear+0.001,0,50)],name='SThouliang')
mdb.models['Model-1'].parts['plate-1'].Set(stringerEdges=(('SThouliang',p.stringers['SThouliang'].edges),), name='houliang')
eee=[e.getByBoundingBox(-50, 0, 1*D1_real-0.00001, 50, 50, 1*D1_real+0.00001), e.getByBoundingBox(-50, -50, 1*D1_real-0.00001, 50, 0, 1*D1_real+0.00001),
e.getByBoundingBox(l1*l+0.0001, -50, 1*D1_real-0.00001, 50, 50, 1*D1_real+0.00001), e.getByBoundingBox(-50, -50, 1*D1_real-0.00001, l1*l-0.0001, 50, 1*D1_real+0.00001)]
p.Stringer(edges=eee,name='jiangqianglei1')
mdb.models['Model-1'].parts['plate-1'].Set(stringerEdges=(('jiangqianglei1',eee),), name='jiangqianglei1')
eee=[e.getByBoundingBox(-50, 0, int(0.5*A1)*D1_real-0.00001, 50, 50, int(0.5*A1)*D1_real+0.00001), e.getByBoundingBox(-50, -50, int(0.5*A1)*D1_real-0.00001, 50, 0, int(0.5*A1)*D1_real+0.00001),
e.getByBoundingBox(l1*l+0.0001, -50, int(0.5*A1)*D1_real-0.00001, 50, 50, int(0.5*A1)*D1_real+0.00001), e.getByBoundingBox(-50, -50, int(0.5*A1)*D1_real-0.00001, l1*l-0.0001, 50, int(0.5*A1)*D1_real+0.00001)]
p.Stringer(edges=eee,name='jiangqianglei2')
mdb.models['Model-1'].parts['plate-1'].Set(stringerEdges=(('jiangqianglei2',eee),), name='jiangqianglei2')
# 定义截面属性 (使用传入的变量)
mdb.models['Model-1'].LProfile(name='changhang', a=st_a, b=st_b, t1=st_t1, t2=st_t2)
mdb.models['Model-1'].RectangularProfile(name='liang', a=bm_a, b=bm_b)
mdb.models['Model-1'].RectangularProfile(name='leibanyuantiao', a=0.5*bm_a, b=0.5*bm_b)
mdb.models['Model-1'].BeamSection(name='liang', integration=DURING_ANALYSIS, poissonRatio=0.0, profile='liang', material='OCr15Ni7Mo2Al', temperatureVar=LINEAR, consistentMassMatrix=False)
mdb.models['Model-1'].BeamSection(name='changhang', integration=DURING_ANALYSIS, poissonRatio=0.0, profile='changhang', material='LY-4', temperatureVar=LINEAR, consistentMassMatrix=False)
mdb.models['Model-1'].BeamSection(name='leibanyuantiao', integration=DURING_ANALYSIS, poissonRatio=0.0, profile='leibanyuantiao', material='LY-4', temperatureVar=LINEAR, consistentMassMatrix=False)
mdb.models['Model-1'].HomogeneousShellSection(name='mp', preIntegrate=OFF, material='LY-12', thicknessType=UNIFORM, thickness=mp_t, integrationRule=SIMPSON, numIntPts=5)
mdb.models['Model-1'].HomogeneousShellSection(name='liangfuban', preIntegrate=OFF, material='LY-9', thicknessType=UNIFORM, thickness=numh1, integrationRule=SIMPSON, numIntPts=5)
mdb.models['Model-1'].HomogeneousShellSection(name='PTleiban', preIntegrate=OFF, material='LY-4', thicknessType=UNIFORM, thickness=numh2, integrationRule=SIMPSON, numIntPts=5)
mdb.models['Model-1'].HomogeneousShellSection(name='JQleiban', preIntegrate=OFF, material='LY-4', thicknessType=UNIFORM, thickness=2*numh2, integrationRule=SIMPSON, numIntPts=5)
# 赋予属性与方向
for set_name, sec_name in [('qianliangfuban','liangfuban'), ('houliangfuban','liangfuban'), ('PTleiban','PTleiban'), ('JQleiban','JQleiban'), ('mengpi','mp'), ('qianliang','liang'), ('changhang1','changhang'), ('changhang2','changhang'), ('jiangqianglei1','leibanyuantiao'), ('jiangqianglei2','leibanyuantiao'), ('houliang','changhang')]:
p.SectionAssignment(region=p.sets[set_name], sectionName=sec_name, offset=0.0, offsetType=MIDDLE_SURFACE, thicknessAssignment=FROM_SECTION)
for set_name, vec in [('qianliang',(1.,0.,0.)), ('changhang1',(-1.,0.,0.)), ('houliang',(-1.,0.,0.)), ('changhang2',(1.,0.,0.)), ('jiangqianglei1',(0.,0.,-1.)), ('jiangqianglei2',(0.,0.,-1.))]:
p.assignBeamSectionOrientation(region=p.sets[set_name], method=N1_COSINES, n1=vec)
# 装配、分析步、相互作用与网格
a = mdb.models['Model-1'].rootAssembly
a.DatumCsysByDefault(CARTESIAN)
a.Instance(name='plate-1-1', part=p, dependent=OFF)
mdb.models['Model-1'].StaticStep(name='Step-1', previous='Initial', stabilizationMagnitude=0.0002, stabilizationMethod=DISSIPATED_ENERGY_FRACTION, adaptiveDampingRatio=0.05, initialInc=0.1, minInc=1e-09, nlgeom=ON)
a.ReferencePoint(point=(0.5 * l, 0.0, 1 * D1_real))
a.ReferencePoint(point=(0.5 * l, 0.0, int(0.5 * A1) * D1_real))
a.ReferencePoint(point=(0.5 * l, 0.0, int(0.3 * A1) * D1_real))
a.ReferencePoint(point=(l*l1+int(0.6*A2)*D2_rear,0,0))
s1, r1 = a.instances['plate-1-1'].faces, a.referencePoints
region1 = regionToolset.Region(referencePoints=(r1[4],))
region2 = regionToolset.Region(side1Faces=s1.getByBoundingBox(-50, -50, 1*D1_real-0.001, 50, 50, 1*D1_real+0.001))
mdb.models['Model-1'].Coupling(name='zaidan1', controlPoint=region1, surface=region2, influenceRadius=WHOLE_SURFACE, couplingType=DISTRIBUTING, weightingMethod=UNIFORM, u1=ON, u2=ON, u3=ON, ur1=ON, ur2=ON, ur3=ON)
mdb.models['Model-1'].rootAssembly.engineeringFeatures.PointMassInertia(name='zaidan1', region=region1, mass=aircraft_mass*0.5*0.05)
region1 = regionToolset.Region(referencePoints=(r1[5],))
region2 = regionToolset.Region(side1Faces=s1.getByBoundingBox(-50, -50, int(0.5*A1)*D1_real-0.001, 50, 50, int(0.5*A1)*D1_real+0.001))
mdb.models['Model-1'].Coupling(name='zaidan2', controlPoint=region1, surface=region2, influenceRadius=WHOLE_SURFACE, couplingType=DISTRIBUTING, weightingMethod=UNIFORM, u1=ON, u2=ON, u3=ON, ur1=ON, ur2=ON, ur3=ON)
mdb.models['Model-1'].rootAssembly.engineeringFeatures.PointMassInertia(name='zaidan2', region=region1, mass=aircraft_mass*0.5*0.05)
region1 = regionToolset.Region(referencePoints=(r1[6],))
region2 = regionToolset.Region(side1Faces=s1.getByBoundingBox(-50, -50, 1 * D1_real - 0.001, 50, 50, 1 * D1_real + 0.001) + s1.getByBoundingBox(-50, -50, int(0.5 * A1) * D1_real - 0.001, 50, 50, int(0.5 * A1) * D1_real + 0.001))
mdb.models['Model-1'].Coupling(name='youxiang', controlPoint=region1, surface=region2, influenceRadius=WHOLE_SURFACE, couplingType=DISTRIBUTING, weightingMethod=UNIFORM, u1=ON, u2=ON, u3=ON, ur1=ON, ur2=ON, ur3=ON)
mdb.models['Model-1'].rootAssembly.engineeringFeatures.PointMassInertia(name='youxiang', region=region1, mass=aircraft_mass*0.5*0.35)
e = a.instances['plate-1-1'].edges
a.Set(edges=e.findAt(((l*l1,0,0),)), name='BC-1')
a.Set(edges=e.findAt(((l*l1+int(0.6*A2)*D2_rear,0,0),)), name='BC-12')
mdb.models['Model-1'].EncastreBC(name='BC-1', createStepName='Initial', region=a.sets['BC-1'])
mdb.models['Model-1'].PinnedBC(name='BC-2', createStepName='Initial', region=regionToolset.Region(referencePoints=(r1[7],)))
mdb.models['Model-1'].Coupling(name='houliangjiaojie', controlPoint=regionToolset.Region(referencePoints=(r1[7],)), surface=a.sets['BC-12'], influenceRadius=WHOLE_SURFACE, couplingType=DISTRIBUTING, weightingMethod=UNIFORM, u1=ON, u2=ON, u3=ON, ur1=ON, ur2=ON, ur3=ON)
mdb.models['Model-1'].Gravity(name='gravity', createStepName='Step-1', comp2=9.8*overload, distributionType=UNIFORM)
a.seedPartInstance(regions=(a.instances['plate-1-1'], ), size=SIZE, deviationFactor=0.1, minSizeFactor=0.1)
a.generateMesh(regions=(a.instances['plate-1-1'], ))
# 创建分析任务
mdb.Job(name=job_name, model='Model-1', type=ANALYSIS, numCpus=2, numDomains=2)
else:
# 模式二：属性更新 (拓扑不变时，仅更新截面尺寸，无需重划网格)
m = mdb.models['Model-1']
# 1. 更新梁和长桁截面
m.profiles['changhang'].setValues(a=st_a, b=st_b, t1=st_t1, t2=st_t2)
m.profiles['liang'].setValues(a=bm_a, b=bm_b)
m.profiles['leibanyuantiao'].setValues(a=0.5*bm_a, b=0.5*bm_b)
# 2. 更新蒙皮和腹板厚度
m.sections['mp'].setValues(thickness=mp_t)
m.sections['liangfuban'].setValues(thickness=numh1)
m.sections['PTleiban'].setValues(thickness=numh2)
m.sections['JQleiban'].setValues(thickness=2*numh2)
# 3. 更新全局变量参数
m.loads['gravity'].setValues(comp2=9.8*overload)
m.rootAssembly.engineeringFeatures.inertias['zaidan1'].setValues(mass=aircraft_mass*0.5*0.05)
m.rootAssembly.engineeringFeatures.inertias['zaidan2'].setValues(mass=aircraft_mass*0.5*0.05)
m.rootAssembly.engineeringFeatures.inertias['youxiang'].setValues(mass=aircraft_mass*0.5*0.35)
# 任务提交与后处理
mdb.jobs[job_name].submit(consistencyChecking=OFF)
mdb.jobs[job_name].waitForCompletion()
# 提取质量和最大mises应力
mass_prop = mdb.models['Model-1'].parts['plate-1'].getMassProperties()
total_mass = mass_prop['mass']
odb = openOdb(path=job_name + '.odb')
stress_field = odb.steps.values()[-1].frames[-1].fieldOutputs['S']
max_mises = max([val.mises for val in stress_field.values if val.mises])
odb.close()
return total_mass, max_mises
# 报错返回
except Exception as e:
raise Exception("Abaqus 仿真引擎执行错误: " + str(e))
优化脚本：optimizer.py
Powered by SLSQP
点击查看能人智工优化方案
import wing_solver
try:
reload(wing_solver)
except NameError:
import importlib
importlib.reload(wing_solver)
import numpy as np
from scipy.optimize import minimize
import csv
import os
ALLOWABLE_STRESS = 300e6
W_span = 5.0
OVERLOAD = 2.0
AIRCRAFT_MASS = 6000.0
CSV_FILENAME = "optimizer_log.csv"
last_D1 = -1.0
last_D2 = -1.0
simulation_cache = {}
iteration_count = 0
def init_csv():
if os.path.exists(CSV_FILENAME):
try:
os.remove(CSV_FILENAME)
except:
pass
with open(CSV_FILENAME, mode='w', newline='') as f:
writer = csv.writer(f)
writer.writerow([
"迭代数", "模式", "肋间距", "长桁间距",
"蒙皮厚", "梁腹厚", "肋腹厚",
"长桁a", "长桁b", "长桁t1", "长桁t2",
"梁a", "梁b",
"质量(kg)", "最大应力(MPa)", "状态"
])
def run_abaqus_cached(x):
global iteration_count, last_D1, last_D2
D1, D2, mp_t, numh1, numh2, st_a, st_b, st_t1, st_t2, bm_a, bm_b = x
# 截断至8位小数作为缓存键值，避免浮点精度误差导致缓存失效
key = tuple(round(v, 8) for v in x)
if key in simulation_cache:
return simulation_cache[key]
iteration_count += 1
if abs(D1 - last_D1) > 1e-5 or abs(D2 - last_D2) > 1e-5:
rebuild_flag = True
mode_text = "[完整建模]"
else:
rebuild_flag = False
mode_text = "[属性更新]"
print("\n[Iter {:<3}] {} 迭代计算中...".format(iteration_count, mode_text))
try:
mass, max_stress = wing_solver.run_wing_simulation(
W=W_span, overload=OVERLOAD, aircraft_mass=AIRCRAFT_MASS,
D1=D1, D2=D2,
mp_t=mp_t, numh1=numh1, numh2=numh2,
st_a=st_a, st_b=st_b, st_t1=st_t1, st_t2=st_t2,
bm_a=bm_a, bm_b=bm_b,
rebuild_geo=rebuild_flag
)
stress_mpa = max_stress / 1e6
note = "安全" if max_stress <= ALLOWABLE_STRESS else "超标"
print(" -> 计算结果: 质量 {:.2f} kg | 应力 {:.2f} MPa | 状态: {}".format(mass, stress_mpa, note))
last_D1 = D1
last_D2 = D2
except Exception as e:
mass, max_stress = 1e6, 1e10
stress_mpa = max_stress / 1e6
note = "特征报错"
print(" -> 计算异常: 几何错误跳过 ({})".format(e))
last_D1 = -1.0
last_D2 = -1.0
simulation_cache[key] = (mass, max_stress)
try:
with open(CSV_FILENAME, mode='a', newline='') as f:
writer = csv.writer(f)
row = [iteration_count, mode_text, round(D1*1000,1), round(D2*1000,1),
round(mp_t*1000,2), round(numh1*1000,2), round(numh2*1000,2),
round(st_a*1000,1), round(st_b*1000,1), round(st_t1*1000,2), round(st_t2*1000,2),
round(bm_a*1000,1), round(bm_b*1000,1),
round(mass,2), round(stress_mpa,2), note]
writer.writerow(row)
except:
pass
return mass, max_stress
def objective(x):
mass, _ = run_abaqus_cached(x)
# 质量目标归一化处理，提高收敛性
return mass / 100.0
def stress_constraint(x):
_, max_stress = run_abaqus_cached(x)
# 应力归一化处理，防止量级差异过大，导致优化器难以收敛
return (ALLOWABLE_STRESS - max_stress) / 1e7
if __name__ == "__main__":
init_csv()
print("="*60)
print(" 机翼结构优化程序启动")
print(" 设定工况: 起飞重量 {}kg | 过载 {}g".format(AIRCRAFT_MASS, OVERLOAD))
print("="*60)
# 初始设计变量，设置较厚的板厚以确保初始构型满足应力约束
x0 = [
0.5, # D1: 肋间距 (初始值 0.5m)
0.15, # D2: 长桁间距 (初始值 0.15m)
0.006, # mp_t: 蒙皮厚度 (初始值 6mm)
0.006, # numh1: 梁腹板厚度 (初始值 6mm)
0.006, # numh2: 肋腹板厚度 (初始值 6mm)
0.03, # st_a: 长桁截面宽度/高度a (初始值 30mm)
0.03, # st_b: 长桁截面宽度/高度b (初始值 30mm)
0.005, # st_t1: 长桁截面厚度t1 (初始值 5mm)
0.005, # st_t2: 长桁截面厚度t2 (初始值 5mm)
0.05, # bm_a: 梁截面宽度/高度a (初始值 50mm)
0.05 # bm_b: 梁截面宽度/高度b (初始值 50mm)
]
bnds = (
(0.500, 0.500), # D1 约束：上下界相同，锁 0.5m，不改变拓扑
(0.150, 0.150), # D2 约束：上下界相同，锁 0.15m，不改变拓扑
(0.002, 0.010), # 蒙皮厚度范围：2mm ~ 10mm
(0.002, 0.020), # 梁腹板厚度范围：2mm ~ 20mm
(0.002, 0.020), # 肋腹板厚度范围：2mm ~ 20mm
(0.010, 0.060), # 长桁尺寸a范围：10mm ~ 60mm
(0.010, 0.060), # 长桁尺寸b范围：10mm ~ 60mm
(0.002, 0.010), # 长桁厚度t1范围：2mm ~ 10mm
(0.002, 0.010), # 长桁厚度t2范围：2mm ~ 10mm
(0.020, 0.100), # 梁截面尺寸a范围：20mm ~ 100mm
(0.020, 0.100) # 梁截面尺寸b范围：20mm ~ 100mm
)
cons = {'type': 'ineq', 'fun': stress_constraint}
res = minimize(
objective, x0, method='SLSQP', bounds=bnds, constraints=cons,
options={'disp': True, 'eps': 1e-4, 'maxiter': 50}
)
print("\n优化任务执行完毕。")
写在最后
这个格式怎么那么难调啊Ovo
不是，为什么会有半期之后直接开始考试月的啊QAQ，博主接下来一个月有四门考试，接下来六周大概率有六门考试......不说了，复习去了。只能之后再试试遗传算法了。
