In this example you will use ABAQUS/CAE to create the model and store the model in a model database. The script opens the model database and performs a parametric study on the model. The example illustrates how you can use a combination of ABAQUS/CAE and the ABAQUS Scripting Interface to analyze a problem.
This example uses ABAQUS Scripting Interface commands to evaluate the sensitivity of the shell elements in ABAQUS to skew distortion when they are used as thin plates. Further details can be found in “Skew sensitivity of shell elements,” Section 2.3.4 of the ABAQUS Benchmarks Manual. The problem investigates the effects on the accuracy of the bending moment computed at the center of a shell using:
different shell formulations and
at different angles.
Figure 7–2 illustrates the basic geometry of the simply supported skew plate with a uniform distributed load.
The plate is loaded by a uniform pressure of 1.0 × 10–6 MPa applied over the entire surface. The edges of the plate are all simply supported. The analysis is performed for five different values of the skew angle,
: 90°, 80°, 60°, 40°, and 30°. The analysis is performed for two different quadrilateral elements: S4 and S8R.The example is divided into two scripts. The controlling script, skewExample.py, imports skewExampleUtils.py. Use the fetch utility to retrieve the scripts:
abaqus fetch job=skewExample abaqus fetch job=skewExampleUtils
The following topics are covered:
You should use ABAQUS/CAE to create your model and to save the resulting model database. You will then use scripting to parameterize your model, submit an analysis job, and operate on the results generated.
Start ABAQUS/CAE, and create a model database from the Start Session dialog box. By default, you are operating on a model named Model-1. The model should include the following:
Part
Create a three-dimensional planar shell part, and name it Plate. Use an approximate size of 5.0. All sides are 1.0 m long.
Material
Create a material, and name it Steel. The Young's modulus is 30 MPa, and the Poisson's ratio is 0.3.
Section
Create a homogeneous shell section that refers to the material called Steel. Name the section Shell. The plate thickness is 0.01 m. The length/thickness ratio is, thus, 100/1 so that the plate is thin in the sense that transverse shear deformation should not be significant. Assign the section to the plate.
Assembly
Create the assembly using a single instance of the part. ABAQUS/CAE names the part instance Plate-1.
Step
Create a static step and name it Step-1. Enter Apply pressure for the step Description. Accept the default time period of 1.0 and the default initial increment of 1.0.
Output database requests
Edit the default output database request for field output and select only U, Translations and rotations, and SF, Section forces and moments, for the whole model after every increment. Delete all requests for history output.
Boundary condition
Create a displacement boundary condition, and name it Pinned. The boundary condition pins the exterior edges of the plate.
Load
Create a pressure load, and name it Pressure. Apply the load to the face of the plate. Accept the default side of the plate and use a magnitude of 1.0. This positive pressure will result in a negative displacement in the 3-direction.
Set
Partition the plate into quarters. Create a set that contains the vertex at the center of the plate, and name the set CENTER.
Mesh
Create a 4 × 4 swept mesh of quadrilateral elements on the plate.
Keyword editor
You must use the Keyword Editor to request output of section forces and moments for the node at the center of the plate. When you edited the output database request to select the variable, SF, ABAQUS/CAE requested output of section forces and moments at element integration points. To modify this output request, you must add position=NODES to the OUTPUT REQUESTS block as follows:
*Element Output, position=NODES SF,
Job
Create a job, and name it skew. The job must refer to the model Model-1.
The parameterized script changes the skew angle of the plate and computes the maximum bending moment at the center for two different element types. The script changes the skew angle by modifying an angular dimension and selecting the vertices to move. You need to add the angular dimensions and determine the indices of the dimensions to modify and the vertices to move.
To add the angular dimensions:
Return to the Part module.
From the main menu bar, select Feature
Edit and select the plate to edit.
From the Edit Feature dialog box, select Edit Section Sketch.
From the Sketcher toolbox, select the angular dimension tool 
and dimension the angles at the lower corners of the plate as shown in Figure 7–3:
To determine the indices of the dimensions to modify and the vertices to move:
From the Sketcher toolbox, select the edit dimension tool 
.
Select the lower left angular dimension.
Select the upper left vertex to move.
Enter a dimension of 60.
Repeat the steps for the lower right dimension and the upper right vertex.
Exit the Sketcher tools, and exit the Sketcher.
From the Edit Feature dialog box, select OK.
Examine the replay file, abaqus.rpy. The last few lines of the replay file will contain the statements that modified the angular dimensions. The statements will look similar to the following:
s1.changeDimension(dimension=d[9], value=60.0,
vertexList=(v[2], ))
s1.changeDimension(dimension=d[10], value=60.0,
vertexList=(v[4], ))The example script, skewExample.py, contains similar statements that modify the angular dimensions of the plate. The indices of the angular dimensions and vertices in your model must be the same as the indices in the example script. If the indices are not the same, you must edit the example script and enter the correct indices.
tmpSketch.changeDimension(dimension=d[9], value=angle,
vertexList=(v[2], ))
tmpSketch.changeDimension(dimension=d[10], value=angle,
vertexList=(v[4], ))Save the model database, and name it skew. ABAQUS/CAE saves the model database in a file called skew.cae. The example script opens this model database and parameterizes the model it contains.
The following shows the contents of the script skewExample.py. The parametric study does the following:
Opens the model database and creates variables that refer to the part, the assembly, and the part instance stored in Model-1.
Creates variables that refer to the four faces and the nine vertices in the instance of the planar shell part.
Skews the plate by modifying the angular dimensions in the sketch of the base feature.
Defines the logical corners of the four faces, and generates a structured mesh.
Runs the analysis for a range of angles using two element types for each angle.
Calculates the maximum moment and displacement at the center of the shell.
Displays X–Y plots in separate viewports of the following:
Displacement versus skew angle
Maximum bending moment versus skew angle
Minimum bending moment versus skew angle
"""
skewExample.py
This script performs a parameter study of element type versus
skew angle. For more details, see Problem 2.3.4 in the
ABAQUS Benchmarks manual.
Before executing this script you must fetch the appropriate
files: abaqus fetch job=skewExample
abaqus fetch job=skewExampleUtils.py
"""
import part
import mesh
from mesh import S4, S8R, STANDARD, STRUCTURED
import job
from skewExampleUtils import getResults, createXYPlot
# 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.), )
face2 = allFaces.findAt((0.,1.,0.), )
face3 = allFaces.findAt((1.,1.,0.), )
face4 = allFaces.findAt((1.,0.,0.), )
allVertices = instance.vertices
v1 = allVertices.findAt((0.,0.,0.), )
v2 = allVertices.findAt((0.,.5,0.), )
v3 = allVertices.findAt((0.,1.,0.), )
v4 = allVertices.findAt((.5,1.,0.), )
v5 = allVertices.findAt((1.,1.,0.), )
v6 = allVertices.findAt((1.,.5,0.), )
v7 = allVertices.findAt((1.,0.,0.), )
v8 = allVertices.findAt((.5,0.,0.), )
v9 = allVertices.findAt((.5,.5,0.), )
# Create a copy of the feature sketch to modify.
tmpSketch = model.Sketch('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, -.001478), (80, -.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.
tmpSketch.changeDimension(dimension=d[9], value=angle,
vertexList=(v[2], ))
tmpSketch.changeDimension(dimension=d[10], value=angle,
vertexList=(v[4], ))
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)
The script imports two functions from skewExampleUtils. The functions do the following:
Retrieve the displacement and calculate the maximum bending moment at the center of the plate.
Display curves of theoretical and computed results in a new viewport.
"""
skewExampleUtils.py
Utilities for the scripting tutorial Skew Example.
"""
from abaqus import *
import visualization
from miscUtils import sorted
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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)
for elemName, xyValues in sorted(data.items()):
xyData = visualization.XYData(elemName, xyValues)
xyPlot.Curve(elemName, xyData)
xyPlot.setValues(curvesToPlot=sorted(data.keys()),
xAxisTitle='Skew Angle', xAxisTitleSource=USER_DEFINED,
yAxisTitle=plotName, yAxisTitleSource=USER_DEFINED)
vp.setValues(displayedObject=xyPlot)