Introduction to Shape Optimization
Finite element (FE) simulations are continuing to evolve and becoming a more and more useful engineering tool. It is easy these days to create very accurate FE models of most physical products in general load environments. One of the newer areas in which FE simulations have become useful is in shape optimizing or topology optimization of a component for some specific application. Recall that topology optimization is a general search for the optimal topology of the component, which can be very different from the original starting geometry and can include holes, etc. Topology optimization typically is based on a linear FE analysis. Shape optimization is in some sense a simpler problem since the topology of the part is not changed, only some of the dimensions are changed in order to minimize an objective function (for example the max stress or weight).
In this example, I will show how you can use MCalibration to simultaneously find the material parameters and the optimal shape of a component!
If you take a square 2D specimen and shear it horizontally such that the shear distance is equal to half the height, then the sides of the specimen do not stay straight. Here are 2 images that show the initial shape and the final shape.
(a) Undeformed shape.
(b) Deformed shape.
(1) What undeformed shapes should the left and right sides of the specimen have in order to become straight after the shearing deformation? (2) What stiffness should the material have in order for the final shearing force to become 10 N?
Solution Step 1: Create FE Model Generation Script
I first created a complete FE model of the specimen geometry, boundary conditions, and materials in Abaqus CAE. To facilitate the geometry optimization I used cubic splines with 4 control points along the left and right sides. Later on, I will have MCalibration search for the optimal locations of these control points.
After finishing the FE model, I quit Abaqus CAE. This creates a file called abaqus.rpy which contains the Python commands that were used to create the model. I renamed this file to Create_FE_Model_Script.py and opened it with a text editor. Most of the commands in this file have now been set up, but there are a few changes that are needed. The first change is to change the definition of the spline control points from being fixed values, to variables that are read in from an external file. There are the Python commands that I used for that:
# read in the geometry parameters
fo = open("parameters.inp", "r")
lines = fo.readlines()
fo.close
a1 = float(lines[6][2:])
a2 = float(lines[7][2:])
a3 = float(lines[8][2:])
a4 = float(lines[9][2:])
b1 = float(lines[10][2:])
b2 = float(lines[11][2:])
b3 = float(lines[12][2:])
b4 = float(lines[13][2:])
pp = lines[16].split(',')
mu = float(pp[0])
kappa = float(pp[1])
s.Spline(points=((0.0, 0.0), (a1, 2.0), (a2, 4.0), (a3, 6.0), (a4, 8.0), (0.0, 10.0)))
s.Spline(points=((10.0, 0.0), (10.0+b1, 2.0), (10.0+b2, 4.0), (10.0+b3, 6.0), (10.0+b4, 8.0), (10.0, 10.0)))
I’m reading in the parameters from a file called parameters.inp. That file will later be automatically created from within MCalibration as a material model template file. In the code above I also obtain the shear modulus and bulk modulus from the files exported from MCalibration. Later, in the Python file I have the following command to instruct Abaqus CAE to create a Neo-Hookean material model with the desired properties (which again are obtained from MCalibration).
mdb.models['Model-1'].materials['mat'].Hyperelastic(materialType=ISOTROPIC,
testData=OFF, type=NEO_HOOKE, volumetricResponse=VOLUMETRIC_DATA, table=((
mu, kappa), ))
Solution Step 2: Create Data Extraction Script
The next step is to create a short Python script that will be used to extract the key results from the FE model once it has been run. The script that I created is listed below. It first opens the Abaqus odb-file, and then creates Python lists for the three nodes sets of interest (left, right, top). If then extracts the coordinates of all nodes along the left and right node sets, and the reaction force along the top node set. Finally, it saves the results in a very specific way in order to facilitate the communication with MCalibration. Since the results file generated by this script will be read by a General External Solver load case in MCalibration, the results file need to have 3 columns corresponding to: (1) generalized time, (2) generalized displacement, (3) generalized force. In this case I selected to generate a results file with 4 rows of data. The first row is just 0, 0, 0. The second row is “1, 1, [error in left side predictions]”. The third row is “2, 2, [error in right side predictions]”, and the forth row is “3, 3, [error in force predictions]”. Using this approach I am able to see how the left, right, and force errors evolve during the optimization.
from odbAccess import *
from time import *
import sys
odb = openOdb(path='Shear_Model.odb')
# get the node sets
nodes_top = []
ns = odb.rootAssembly.instances['SPECIMEN'].nodeSets['TOP']
for i in range(len(ns.nodes)):
nodes_top.append(ns.nodes[i].label)
nodes_left = []
ns = odb.rootAssembly.instances['SPECIMEN'].nodeSets['LEFT']
for i in range(len(ns.nodes)):
nodes_left.append(ns.nodes[i].label)
nodes_right = []
ns = odb.rootAssembly.instances['SPECIMEN'].nodeSets['RIGHT']
for i in range(len(ns.nodes)):
nodes_right.append(ns.nodes[i].label)
# extract the history variables
step = odb.steps['theStep']
errL = 1.0
for nn in nodes_left:
name = 'Node SPECIMEN.' + str(nn)
LX = step.historyRegions[name].historyOutputs['COOR1'].data[-1][-1]
LY = step.historyRegions[name].historyOutputs['COOR2'].data[-1][-1]
errL += abs(LY - 2.0 * LX)
errR = 1.0
for nn in nodes_right:
name = 'Node SPECIMEN.' + str(nn)
RX = step.historyRegions[name].historyOutputs['COOR1'].data[-1][-1] - 10.0
RY = step.historyRegions[name].historyOutputs['COOR2'].data[-1][-1]
errR += abs(RY - 2.0 * RX)
RFtot = 0.0
for nn in nodes_top:
name = 'Node SPECIMEN.' + str(nn)
RFtot += step.historyRegions[name].historyOutputs['RF1'].data[-1][-1]
errF = 1.0 + abs(RFtot - 10.0)
# save the results
outf = open("data.txt", 'w')
outf.write('0.0, 0.0, 1.0\n')
outf.write('1.0, 1.0, ' + str(errL) + '\n')
outf.write('2.0, 2.0, ' + str(errR) + '\n')
outf.write('3.0, 3.0, ' + str(errF) + '\n')
outf.close()
odb.close()
Solution Step 3: Create the MCalibration File
First, setup the material model template in MCalibration. In this case I used an Abaqus-Template material with the definition shown below. I define the variables [a1, a2, a3, a4] which correspond to the horizontal positions of the control points of cubic spline defining the left side of the specimen. I also define the variables [b1, b2, b3, b4] which correspond to the horizontal positions of the control points of the cubic spline defining the right side of the specimen. Finally, I define the Neo-Hookean shear modulus as a variable. Note that MCalibration will export these variables into a file that will be used by the model generation script.
Now we can set up MCalibration to use the model generation and data extraction scripts. The setup that I used is shown below. Note that I left the data extraction command to be empty since that is performed by the Setup_and_Run.py script.
The target_data.txt file has the following values:
0.0, 0.0, 1.0
1.0, 1.0, 1.0
2.0, 2.0, 1.0
3.0, 3.0, 1.0
Solution Step 4: Run the Shape Optimization
Everything is now ready to run. Here is the results from Run Once in MCalibration.
The final results after running the optimization are shown below.
The figure to the left shows the optimized undeformed shape, and the figure to the right shows the deformed shape. As perhaps expected, if the initial shape is slightly narrower in the middle then the deformed (sheared) shape gets straight edges. I have not showed it here, but the applied force reaches the target 10 N.
That is is! This example illustrated how you can use MCalibration to not only calibrate a material model, but also perform a shape optimization given any objective function you are interested in.
