基于Python的结构优化？ 50
浏览：1170
abaqus二次开发，源程序由我提供，但是需要进行改变，将单材料的优化程序改为多材料的优化程序，基于Python语言，源代码写的很清楚，改动不会太大，但是会附带后处理的问题解决，详情细聊。
import math,customKernel
from abaqus import getInput,getInputs
from odbAccess import openOdb
## Function of formatting Abaqus model for stiffness optimisation
def fmtMdb(Mdb):
mdl = Mdb.models['Model-1']
part = mdl.parts['Part-1']
# Build sections and assign solid section
mdl.Material('Material01').Elastic(((1.0, 0.3), ))
mdl.HomogeneousSolidSection('sldSec','Material01')
mdl.Material('Material02').Elastic(((0.001**3, 0.3), ))
mdl.HomogeneousSolidSection('voidSec','Material02')
part.SectionAssignment(part.Set('ss',part.elements),'sldSec')
# Define output request
mdl.FieldOutputRequest('SEDensity','Step-1',variables=('ELEDEN', ))
mdl.HistoryOutputRequest('ExtWork','Step-1',variables=('ALLWK', ))
## Function of running FEA for raw sensitivities and objective function
def FEA(Iter,Mdb,Xe,Ae):
Mdb.Job('Design_Job'+str(Iter),'Model-1').submit()
Mdb.jobs['Design_Job'+str(Iter)].waitForCompletion()
opdb = openOdb('Design_Job'+str(Iter)+'.odb')
seng = opdb.steps['Step-1'].frames[-1].fieldOutputs['ESEDEN'].values
for en in seng: Ae[en.elementLabel]=en.data/Xe[en.elementLabel]
obj=opdb.steps['Step-1'].historyRegions['Assembly ASSEMBLY'].historyOutputs['ALLWK'].data[-1][1]
opdb.close()
return obj
## Function of preparing filter map (Fm={elm1:[[el1,el2,...],[wf1,wf2,...]],...})
def preFlt(Rmin,Elmts,Nds,Fm):
# Calculate element centre coordinates
c0 = {}
for el in Elmts:
nds = el.connectivity
c0[el.label]=[sum([Nds[nd].coordinates[i]/len(nds) for nd in nds]) for i in range(3)]
# Weighting factors
for el in Elmts:
Fm[el.label] = [[],[]]
for em in Elmts:
dis=math.sqrt(sum([(c0[el.label][i]-c0[em.label][i])**2 for i in range(3)]))
if dis<Rmin:
Fm[el.label][0].append(em.label)
Fm[el.label][1].append(Rmin - dis)
sm = sum(Fm[el.label][1])
for i in range(len(Fm[el.label][0])): Fm[el.label][1][i] /= sm
## Function of filtering sensitivities
def fltAe(Ae,Fm):
raw = Ae.copy()
for el in Fm.keys():
Ae[el] = 0.0
for i in range(len(Fm[el][0])): Ae[el]+=raw[Fm[el][0][i]]*Fm[el][1][i]
## Function of optimality update for design variables and Abaqus model
def BESO(Vf,Xe,Ae,Part,Elmts):
lo, hi = min(Ae.values()), max(Ae.values())
tv = Vf*len(Elmts)
while (hi-lo)/hi > 1.0e-5:
th = (lo+hi)/2.0
for key in Xe.keys(): Xe[key] = 1.0 if Ae[key]>th else 0.001
if sum(Xe.values())-tv>0: lo = th
else: hi = th
# Label elements as solid or void
vlb, slb = [], []
for el in Elmts:
if Xe[el.label] == 1.0: slb.append(el.label)
else: vlb.append(el.label)
# Assign solid and void elements to each section
Part.SectionAssignment(Part.SetFromElementLabels('ss',slb),'sldSec')
Part.SectionAssignment(Part.SetFromElementLabels('vs',vlb),'voidSec')
## ====== MAIN PROGRAM ======
if __name__ == '__main__':
# Set parameters and inputs
pars = (('VolFrac:','0.5'), ('Rmin:', '1'), ('ER:', '0.02'))
vf,rmin,ert = [float(k) if k!=None else 0 for k in getInputs(pars,dialogTitle='Parameters')]
if vf<=0 or rmin<0 or ert<=0: sys.exit()
mddb = openMdb(getInput('Input CAE file:',default='Test.cae'))
# Design initialization
fmtMdb(mddb)
part = mddb.models['Model-1'].parts['Part-1']
elmts, nds = part.elements, part.nodes
oh, vh = [], []
xe, ae, oae, fm = {}, {}, {}, {}
for el in elmts: xe[el.label] = 1.0
if rmin>0: preFlt(rmin,elmts,nds,fm)
# Optimisation iteration
change, iter, obj = 1, -1, 0
while change > 0.001:
iter += 1
# Run FEA
oh.append(FEA(iter,mddb,xe,ae))
# Process sensitivities
if rmin>0: fltAe(ae,fm)
if iter > 0: ae=dict([(k,(ae[k]+oae[k])/2.0) for k in ae.keys()])
oae = ae.copy()
# BESO optimisation
vh.append(sum(xe.values())/len(xe))
nv = max(vf,vh[-1]*(1.0-ert))
BESO(nv,xe,ae,part,elmts)
if iter>10: change=math.fabs((sum(oh[iter-4:iter+1])-sum(oh[iter-9:iter-4]))/sum(oh[iter-9:iter-4]))
# Save results
mddb.customData.History = {'vol':vh,'obj':oh}
mddb.saveAs('Final_design.cae')
