备注
Go to the end to download the full example code.
参数研究¶
skewExpleple.py 来自 Using a script to perform a parametric study。
该脚本执行单元类型与倾斜角度的参数研究。 有关更多详细信息，请参阅 《Abaqus Benchmarks manual》中的问题 2.3.4。
在运行此脚本之前运行以下命令：
abaqus fetch job=viewer_tutorial
以下是 trakewExpleple.py 的源代码。参数研究完成以下工作：
打开模型数据库并创建引用存储在 Model-1 中的部件、装配体和部件实例的变量。
创建引用平面壳体零件实例中的四个面和九个顶点的变量。
通过修改基础特征草图中的角度尺寸来倾斜板。
定义四个面的逻辑角，并生成结构化网格。
对每个角度使用两种元素类型来运行一系列角度的分析。
计算壳体中心处的最大力矩和位移。
在以下各项的单独视口中显示 X - Y 绘图：
位移与倾斜角
最大弯矩与倾斜角度
最小弯矩与倾斜角
还绘制了理论结果。
import visualization
from abaqus import *
from abaqusConstants import *
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def getResults():
"""
Retrieve the displacement and calculate the minimum
and maximum bending moment at the center of plate.
"""
from visualization import ELEMENT_NODAL
# Open the output database.
odb = visualization.openOdb("skew.odb")
centerNSet = odb.rootAssembly.nodeSets["CENTER"]
frame = odb.steps["Step-1"].frames[-1]
# Retrieve Z-displacement at the center of the plate.
dispField = frame.fieldOutputs["U"]
dispSubField = dispField.getSubset(region=centerNSet)
disp = dispSubField.values[0].data[2]
# Average the contribution from each element to the moment,
# then calculate the minimum and maximum bending moment at
# the center of the plate using Mohr's circle.
momentField = frame.fieldOutputs["SM"]
momentSubField = momentField.getSubset(region=centerNSet, position=ELEMENT_NODAL)
m1, m2, m3 = 0, 0, 0
for value in momentSubField.values:
m1 = m1 + value.data[0]
m2 = m2 + value.data[1]
m3 = m3 + value.data[2]
numElements = len(momentSubField.values)
m1 = m1 / numElements
m2 = m2 / numElements
m3 = m3 / numElements
momentA = 0.5 * (abs(m1) + abs(m2))
momentB = sqrt(0.25 * (m1 - m2) ** 2 + m3**2)
maxMoment = momentA + momentB
minMoment = momentA - momentB
odb.close()
return disp, maxMoment, minMoment
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def createXYPlot(vpOrigin, vpName, plotName, data):
"""
Display curves of theoretical and computed results in
a new viewport.
"""
from visualization import USER_DEFINED
vp = session.Viewport(name=vpName, origin=vpOrigin, width=150, height=100)
xyPlot = session.XYPlot(plotName)
chart = xyPlot.charts.values()[0]
curveList = []
for elemName, xyValues in sorted(data.items()):
xyData = session.XYData(elemName, xyValues)
curve = session.Curve(xyData)
curveList.append(curve)
chart.setValues(curvesToPlot=curveList)
chart.axes1[0].axisData.setValues(useSystemTitle=False, title="Skew Angle")
chart.axes2[0].axisData.setValues(useSystemTitle=False, title=plotName)
vp.setValues(displayedObject=xyPlot)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def createModel():
"""
Create the skew example model, including material, step, load, bc, and job.
"""
import mesh
import part
import regionToolset
import step
# Create the Plate
m = mdb.models["Model-1"]
s = m.ConstrainedSketch(name="__profile__", sheetSize=5.0)
g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
s.sketchOptions.setValues(
sheetSize=5.0,
gridSpacing=0.1,
grid=ON,
gridFrequency=2,
constructionGeometry=ON,
dimensionTextHeight=0.1,
decimalPlaces=2,
)
s.setPrimaryObject(option=STANDALONE)
s.rectangle(point1=(0.0, 0.0), point2=(1.0, 1.0))
s.delete(objectList=(c[21], c[18], c[19], c[20]))
s.HorizontalConstraint(entity=g.findAt((0.5, 0.0)))
s.FixedConstraint(entity=v.findAt((0.0, 0.0)))
s.FixedConstraint(entity=v.findAt((1.0, 0.0)))
s.ParallelConstraint(entity1=g.findAt((0.0, 0.5)), entity2=g.findAt((1.0, 0.5)))
s.AngularDimension(line1=g.findAt((0.0, 0.5)), line2=g.findAt((0.5, 0.0)), textPoint=(0.2, 0.2), value=90.0)
p = m.Part(name="Plate", dimensionality=THREE_D, type=DEFORMABLE_BODY)
p.BaseShell(sketch=s)
s.unsetPrimaryObject()
vp = session.viewports["Viewport: 1"]
vp.setValues(displayedObject=p)
del mdb.models["Model-1"].sketches["__profile__"]
# Create the Steel material
m.Material("Steel")
m.materials["Steel"].Elastic(table=((30.0e6, 0.3),))
m.HomogeneousShellSection(
name="Shell",
preIntegrate=OFF,
material="Steel",
thickness=0.01,
poissonDefinition=DEFAULT,
temperature=GRADIENT,
integrationRule=SIMPSON,
numIntPts=5,
)
# Assign Steel to the plate
p = mdb.models["Model-1"].parts["Plate"]
region = (None, None, p.faces, None)
p.SectionAssignment(region=region, sectionName="Shell")
# Create the assembly
a = m.rootAssembly
vp.setValues(displayedObject=a)
a.DatumCsysByDefault(CARTESIAN)
a.Instance(name="Plate-1", part=p, dependent=OFF)
pi = a.instances["Plate-1"]
# Create the step
m.StaticStep(name="Step-1", previous="Initial", description="Apply pressure", timePeriod=1, initialInc=1)
vp.assemblyDisplay.setValues(step="Step-1")
m.fieldOutputRequests["F-Output-1"].setValues(frequency=1, variables=("U",))
m.FieldOutputRequest(name="F-Output-2", createStepName="Step-1", variables=("SF",), position=NODES)
del mdb.models["Model-1"].historyOutputRequests["H-Output-1"]
# Create the displacement BC
e = pi.edges
edges = e.findAt(
((0.25, 0.0, 0.0),),
((1.0, 0.25, 0.0),),
((0.75, 1.0, 0.0),),
((0.0, 0.75, 0.0),),
)
region = (None, edges, None, None)
m.DisplacementBC(name="Pinned", createStepName="Step-1", region=region, u1=0.0, u2=0.0, u3=0.0)
# Create the Pressure load
s1 = pi.faces
side1Faces1 = s1.findAt(
(
(0.333333333333333, 0.333333333333333, 0.0),
(0.0, 0.0, 1.0),
),
)
region = regionToolset.Region(side1Faces=side1Faces1)
m.Pressure(
name="Load-1", createStepName="Step-1", region=region, distributionType=UNIFORM, magnitude=1.0, amplitude=UNSET
)
# Partition the face
f1, e1 = pi.faces, pi.edges
faces = (f1.findAt(coordinates=(0.33333333333, 0.33333333333, 0.0)),)
pt1 = pi.InterestingPoint(edge=e1.findAt(coordinates=(0.0, 0.75, 0.0)), rule=MIDDLE)
pt2 = pi.InterestingPoint(edge=e1.findAt(coordinates=(1.0, 0.25, 0.0)), rule=MIDDLE)
a.PartitionFaceByShortestPath(faces=faces, point1=pt1, point2=pt2)
faces = (
f1.findAt(coordinates=(0.33333333333, 0.66666666667, 0.0)),
f1.findAt(coordinates=(0.66666666667, 0.33333333333, 0.0)),
)
pt1 = pi.InterestingPoint(edge=e1.findAt(coordinates=(0.75, 1.0, 0.0)), rule=MIDDLE)
pt2 = pi.InterestingPoint(edge=e1.findAt(coordinates=(0.25, 0.0, 0.0)), rule=MIDDLE)
a.PartitionFaceByShortestPath(faces=faces, point1=pt1, point2=pt2)
# Create the Geometry set CENTER
verts = pi.vertices.findAt(((0.5, 0.5, 0.0),))
a.Set(name="CENTER", vertices=verts)
# Create the mesh
a.seedPartInstance(regions=(pi,), size=0.25)
a.generateMesh(regions=(pi,))
# Create the job
mdb.Job(
name="skew",
model="Model-1",
type=ANALYSIS,
explicitPrecision=SINGLE,
description="",
userSubroutine="",
numCpus=1,
scratch="",
echoPrint=OFF,
modelPrint=OFF,
contactPrint=OFF,
historyPrint=OFF,
)
mdb.saveAs(pathName="skew.cae")
createModel()
import job
import mesh
import part
from mesh import S4, S8R, STANDARD, STRUCTURED
# Create a list of angle parameters and a list of
# element type parameters.
angles = [90, 80, 60, 40, 30]
elemTypeCodes = [S4, S8R]
# Open the model database.
openMdb("skew.cae")
model = mdb.models["Model-1"]
part = model.parts["Plate"]
feature = part.features["Shell planar-1"]
assembly = model.rootAssembly
instance = assembly.instances["Plate-1"]
job = mdb.jobs["skew"]
allFaces = instance.faces
regions = (allFaces[0], allFaces[1], allFaces[2], allFaces[3])
assembly.setMeshControls(regions=regions, technique=STRUCTURED)
face1 = allFaces.findAt(
(0.0, 0.0, 0.0),
)
face2 = allFaces.findAt(
(0.0, 1.0, 0.0),
)
face3 = allFaces.findAt(
(1.0, 1.0, 0.0),
)
face4 = allFaces.findAt(
(1.0, 0.0, 0.0),
)
allVertices = instance.vertices
v1 = allVertices.findAt(
(0.0, 0.0, 0.0),
)
v2 = allVertices.findAt(
(0.0, 0.5, 0.0),
)
v3 = allVertices.findAt(
(0.0, 1.0, 0.0),
)
v4 = allVertices.findAt(
(0.5, 1.0, 0.0),
)
v5 = allVertices.findAt(
(1.0, 1.0, 0.0),
)
v6 = allVertices.findAt(
(1.0, 0.5, 0.0),
)
v7 = allVertices.findAt(
(1.0, 0.0, 0.0),
)
v8 = allVertices.findAt(
(0.5, 0.0, 0.0),
)
v9 = allVertices.findAt(
(0.5, 0.5, 0.0),
)
# Create a copy of the feature sketch to modify.
tmpSketch = model.ConstrainedSketch("tmp", feature.sketch)
v, d = tmpSketch.vertices, tmpSketch.dimensions
# Create some dictionaries to hold results. Seed the
# dictionaries with the theoretical results.
dispData, maxMomentData, minMomentData = {}, {}, {}
dispData["Theoretical"] = ((90, -0.001478), (80, -0.001409), (60, -0.000932), (40, -0.000349), (30, -0.000148))
maxMomentData["Theoretical"] = ((90, 0.0479), (80, 0.0486), (60, 0.0425), (40, 0.0281), (30, 0.0191))
minMomentData["Theoretical"] = ((90, 0.0479), (80, 0.0448), (60, 0.0333), (40, 0.0180), (30, 0.0108))
# Loop over the parameters to perform the parameter study.
for elemCode in elemTypeCodes:
# Convert the element type codes to strings.
elemName = repr(elemCode)
dispData[elemName], maxMomentData[elemName], minMomentData[elemName] = [], [], []
# Set the element type.
elemType = mesh.ElemType(elemCode=elemCode, elemLibrary=STANDARD)
assembly.setElementType(regions=(instance.faces,), elemTypes=(elemType,))
for angle in angles:
# Skew the geometry and regenerate the mesh.
assembly.deleteMesh(regions=(instance,))
d[0].setValues(
value=angle,
)
feature.setValues(sketch=tmpSketch)
part.regenerate()
assembly.regenerate()
assembly.setLogicalCorners(region=face1, corners=(v1, v2, v9, v8))
assembly.setLogicalCorners(region=face2, corners=(v2, v3, v4, v9))
assembly.setLogicalCorners(region=face3, corners=(v9, v4, v5, v6))
assembly.setLogicalCorners(region=face4, corners=(v8, v9, v6, v7))
assembly.generateMesh(regions=(instance,))
# Run the job, then process the results.
job.submit()
job.waitForCompletion()
print("Completed job for %s at %s degrees" % (elemName, angle))
disp, maxMoment, minMoment = getResults()
dispData[elemName].append((angle, disp))
maxMomentData[elemName].append((angle, maxMoment))
minMomentData[elemName].append((angle, minMoment))
# Plot the results.
createXYPlot((10, 10), "Skew 1", "Displacement - 4x4 Mesh", dispData)
createXYPlot((160, 10), "Skew 2", "Max Moment - 4x4 Mesh", maxMomentData)
createXYPlot((310, 10), "Skew 3", "Min Moment - 4x4 Mesh", minMomentData)
由 Sphinx-Gallery 生成