Я хочу выполнить параметры c исследования оболочки цилиндра с использованием abaqus и python. elemRaduis и elemHeight - это параметры, которые я собираюсь изменить. Но, к сожалению, я не могу получить желаемые значения, деформацию для начала, из файла odb. Я хочу получить усреднение деформации в специальных элементах (хранящихся в списке), сохраненных в списке.
Заранее спасибо
'''
from abaqus import *
from abaqusConstants import *
import regionToolset
import numpy as np
#Define the function used to calculate in Aabaqus
def calc_cylinder(radius, height, thickness, eModul, poisson, pressure, elemHeight, elemRadius, elemType):
#---------------------------------------------
#Create the part
import sketch
import part
#Sketch the topview of the quarter cylinder
cylinderSketch = cylinder.ConstrainedSketch(name='Cylinder SKetch',
sheetSize = max(radius, height) *1.1)
cylinderSketch.ArcByCenterEnds(center=(0.0, 0.0), point1=(radius, 0.0), point2=(0.0, radius))
#From the created sketch the Part() method ist used to create a 3D part out of the sketch
cylinderPart = cylinder.Part(name='Cylinder Part',
dimensionality = THREE_D,
type = DEFORMABLE_BODY)
#extrude the topview sketch/part by the height of the cylinder
cylinderPart.BaseShellExtrude(depth = height,
sketch = cylinderSketch)
#---------------------------------------------
#Create the Material
import material
cylinderMaterial = cylinder.Material(name='material1')
cylinderMaterial.Elastic(table=((eModul, poisson), ))
#---------------------------------------------
#Create and assign the section
import section
cylinderSection = cylinder.HomogeneousShellSection(name= 'Cylinder Section',
material = 'material1',
thicknessType = UNIFORM,
thickness = thickness)
#Assign the section th the surface
cylinderCenterPonit = (np.cos(45)*radius, np.sin(45)*radius, height/2)
cylinderFace = cylinderPart.faces.findAt((cylinderCenterPonit,))
cylinderRegion = (cylinderFace, )
cylinderPart.SectionAssignment(region=cylinderRegion,
sectionName = 'Cylinder Section',
offset = 0.0,
offsetType = MIDDLE_SURFACE,
offsetField = '')
#---------------------------------------------
#Create assembly
import assembly
cylinderAssembly = cylinder.rootAssembly
cylinderInstance = cylinderAssembly.Instance(name='Cylinder Instance',
part = cylinderPart,
dependent = ON)
#--------------------------------------------
#Create the steps
import step
cylinder.StaticStep(name='Load Step', #or 'Step-1
previous='Initial',
description='Apply pressure on the inside of the shell ',
nlgeom=OFF)
#---------------------------------------------
#Create the field output request
#---------------------------------------------
#Apply Boundary Conditions
#Clamped Edge on bottom of cylinder
clampedEdge = cylinderInstance.edges.findAt(((np.cos(45)*radius, np.sin(45)*radius, 0), ))
clampedRegion = regionToolset.Region(edges=clampedEdge)
cylinder.EncastreBC(name = 'Encaster Edge',
createStepName = 'Initial',
region = clampedRegion)
#Edge on the x Axis
xAxisEdge = cylinderInstance.edges.findAt(((radius, 0, height/2), ))
xAxisRegion = regionToolset.Region(edges=xAxisEdge)
cylinder.DisplacementBC(name = 'Support on x-Axis',
createStepName = 'Initial',
region = xAxisRegion,
u1 = UNSET, u2 = SET, u3 = UNSET,
ur1 = SET, ur2 = UNSET, ur3 = SET,
amplitude = UNSET,
distributionType = UNIFORM,
fieldName = '',
localCsys = None)
#Edge on the y Axis
yAxisEdge = cylinderInstance.edges.findAt(((0, radius, height/2), ))
yAxisRegion = regionToolset.Region(edges=yAxisEdge)
cylinder.DisplacementBC(name = 'Support on y-Axis',
createStepName = 'Initial',
region = yAxisRegion,
u1 = SET, u2 = UNSET, u3 = UNSET,
ur1 = UNSET, ur2 = SET, ur3 = SET,
amplitude = UNSET,
distributionType = UNIFORM,
fieldName = '',
localCsys = None)
#---------------------------------------------
#Apply Loads
shellFace = cylinderInstance.faces.findAt(((np.cos(45)*radius, np.sin(45)*radius, height/2), ))
pressureSurface = regionToolset.Region(side2Faces = shellFace)
cylinder.Pressure(name = 'Apply Pressure',
createStepName = 'Load Step',
region = pressureSurface,
magnitude = pressure)
#---------------------------------------------
#Create the mesh
import mesh
cylinderMeshRegion = cylinderRegion
cylinderElemType = mesh.ElemType(elemCode=SymbolicConstant(elemType), elemLibrary=STANDARD)
cylinderPart.setElementType(regions = cylinderMeshRegion,
elemTypes = (cylinderElemType, ))
#Seed the edges by number to perform the parametric studies
meshEdgesVertical = cylinderPart.edges.findAt(((radius, 0.0, height/2),),
((0.0, radius, height /2),))
meshEdgesHorizontal = cylinderPart.edges.findAt(((np.cos(45)*radius, np.sin(45)*radius, 0.0),),
((np.cos(45)*radius, np.sin(45)*radius, height),))
cylinderPart.seedEdgeByNumber (edges = meshEdgesHorizontal, number = elemRadius)
cylinderPart.seedEdgeByNumber (edges = meshEdgesVertical, number = elemHeight)
cylinderPart.generateMesh()
#---------------------------------------------
#Stundet version cant handle more than 1000 nodes
numNodes = len(cylinderPart.nodes)
#If the mesh has more than 1000 nodes, return an error
if numNodes > 1000:
output = 'over 1000 nodes'
#If less than 1000 nodes, create josb and calculate them
else:
#---------------------------------------------
#Create Jobs
import job
#Create the jobs
job_name = ('Cylinder_' + elemType + '_' + str(elemRadius) + '_' + str(elemHeight))
mdb.Job(name = job_name,
model = 'Cylinder',
type = ANALYSIS,
description = 'Linear Analysis')
#run the job
mdb.jobs[job_name].submit()
#wait for the job to complete
mdb.jobs[job_name].waitForCompletion()
#---------------------------------------------
#Find needed displacemnt, moment etc
odb = session.openOdb(name=job_name + '.odb')
#open step
odb_step = odb.steps['Load Step']
#open first frame, not zeroth frame, for geom. nl analysis probably different frame
odb_frame = odb_step.frames[-1]
#open needed variable
odb_output = odb_frame.fieldOutputs['U']
#Get the elemen numbers nin an array
start = int(round(elemRadius / 2))
stop = int(round((elemRadius + 1) * (elemHeight + 1)))
step = int(round(elemRadius + 1))
elem_numbers = range(start, stop, step)
#leere Liste erezugen fuer die Verformung
deformation_u = list()
#durch die Elemente schleifen und die Verformungen in die liste schreiben
for i in elem_numbers:
odb_knot = odb_output.values
deformation_u.append(odb_knot.data)
odb_return = deformation_u
# Close Database
odb.close()
return rtrn
#Create Viewport in Abaqus and display nothing
session.viewports['Viewport: 1'].setValues(displayedObject=None)
#Create the model
mdb.models.changeKey(fromName='Model-1', toName= 'Cylinder' )
cylinder = mdb.models['Cylinder']
#Define parameters for the calculation
radius = 4
height = 8
thickness = 0.2
eModul = 34000000
poisson = 0.2
pressure = 10
elemRadius = [16] #[1, 2, 4, 8, 16, 32, 64]
elemHeight = [16] #[1, 2, 4, 8, 16, 32, 64]
elemType = ['S4'] #['S4', 'S4R', 'S8R']
results = list()
#Loop through the parameters and call the calculation function
for elem in elemType:
for heig in elemHeight:
for rad in elemRadius:
output = calc_cylinder(radius, height, thickness, eModul, poisson, pressure, heig, rad, elem)
#Define Str for results
result = elem + '\t' + str(heig) + '\t' + str(rad) + '\t' + str(list(output))
results.append(result)
'''