Fig. 1. left to right physical model of curved crease model created by Sam; bottom view of model crushed in experiment taken by Sam; one unit of curved crease model undeformed in Abaqus simulation UI; one unit of curved crease model crushed in simulation
Energy-absorbing materials are used in a variety of applications ranging from crash panels in cars to the soles of shoes. Currently, the industry standard for these types of materials are foams of different types and specific types of cellular structures like honeycomb sandwich panels. In this project, we try to create origami and kirigami structures (folded and cut from sheets of material) that have energy-absorbing properties comparable to the industry standard. One advantage of foldables is that these structures are tuneable as opposed to monolithic structures. Because 2D processes can be used to make these structures, they are also more efficient to manufacture as compared to the industry standard.
These foldable cellular structures are defined by several different parameters depending on the specific structure. To optimize for these different parameters, it is useful to simulate on these and relatively compare the effect of each parameter on the energy absorption. We use Abaqus for this purpose. Using Abaqus' Python API, we are able to run several simulations very quickly while ensuring that the boundary conditions and loads remain accurate and consistent throughout the models.
This particular application is a collaboration with Toyota Motor Corporation with the aim of improving energy absorbing capability per unit mass over existing technologies used in vehicle bodies. Typically, aluminum honeycombs are used because they are easily available and provide reasonable performance. Theoretical work has shown (Gattas and You, 2015) that folded cellular materials (particularly those with curved creases) can greatly exceed the performance of honeycombs, but experimental validation only showed modest improvement, mostly due to inadequate manufacturing methods. Our work will rely on a method developed by Sam which can accurately reproduce curved crease structures and make use of advanced composite materials. Our goal will be to experimentally produce cellular materials with much greater energy absorption per mass than aluminum honeycombs.
My python script takes in several inputs to run a simulation of an applied force. In order, they are the step file, the number of increments, the amount of force applied per increment, the thickness of the shell section, the name of the output data file, and the boundary condition clearance. The Abaqus scripting manual can be found here but this manual is mostly a repository of commands and does not explain how to set up a simulation. This document explains the set up of an Abaqus script better but only sets up a very basic simulation. The documentation below details the commands used specifically in the simulation I set up which constrains the bottom of any model and applies a force to the top nodes of the model. A periodic boundary condition is also applied tying the front to the back of the model in this specific simulation. This simulation attempts to model the real-life experiment of applying a force across the top of the model as the model sits on top of a stationary surface. The experimental setup that Sam uses is shown below.
Fig. 2. Experimental setup used to crush models (taken from Sam's documentation on energy absorption)
The code for setting up the simulation can be found here while the code for graphing the results can be found here.
As mentioned before, we used curved-crease Miura-Ori folds for this energy absorption application. After experimenting with Abaqus' internal geometry creation functions, I decided the best way to simulate in Abaqus with origami and kirigami in particular was to import already-constructed geometries into Abaqus. For this purpose, models were generated using a python script written by Sam which can be accessed here. This script generates STEP files with the model information based on two different parameters: radius of cylinder section (R), and the fold angle between the rule lines and the plane of the curved crease (gamma). For initial simulation, radius values of 5, 7, 9, 11, 13, and 15 cm and gamma values of 30, 45 and 60 degrees were used.
Integrating the STEP file into Abaqus involves first opening the STEP file and then creating a part from the geometry. In Abaqus, parts are individual components that are not held in the same coordinate space. If there are multiple components moving in relation to one another in the simulation, each component must be imported or created as separate parts that exist in their own coordinate space. These parts are placed in relation to each other in a later step. In this specific case, for the purpose of importing STEP files, if the model has multiple different faces, the input stitch tolerance must be small enough to stitch those faces into one part. Otherwise, the model will be imported into Abaqus as several parts where each part is a different face of your original model. In this step, geometries can also be set to 2-dimensional or 3-dimensional. For this simulation, we set the dimensionality to "THREE-D". Below are a list of python commands to perform the operations described.
Command used to import the STEP file: mdb.openStep('C:\Users\Soma\Documents\ea_sim' + filename, scaleFromFile=OFF) where filename is the name of your STEP file and is set in command line.
Command used to create part from geometry: mdb.models['Model-1'].PartFromGeometryFile(combine=True, dimensionality=THREE_D, geometryFile=mdb.acis, name='simple-cc-test', type=DEFORMABLE_BODY, stitchTolerance=1e-6).
After the creation of parts, the materials and material properties you plan to use must be defined. These material properties can fall under a variety of categories including General (density, depvar, regularization, etc.), Mechanical (elasticity, plasticity, damping, viscosity, etc.), Thermal (conductivity, heat generation, specific heat, etc.), Electrical/Magnetic (electrical conductivity, dielectric, piezoelectric, magnetic permeability), and Other (acoustic medium, mass diffusion, etc.). Because we are mainly performing structural simulations, we only used elasticity and plasticity fields for our material properties. In elasticity, a Young's modulus and Poisson's ratio are required. For the plasticity field, a table of values is required which specifies different stress and strain values. Plasticity is extrapolated from that data. Python commands are described below. At this point in the simulation, I have already set model = mdb.models['Model-1']
Command used to create a new material: model.Material(name='aluminum')
Command used to set elasticity: model.materials['aluminum'].Elastic(table=((70000000000.0, 0.3), ))
Command used to set plasticity: model.materials['aluminum'].Plastic(table=((300000000.0,0), (300000000.0, 0.06), ))
For our purposes of simulating folded structures, we decided to use a homogenous shell section to simulate the sheet of material. My script allows for the model to be set to any input thickness. For most of our simulation tests, we used a thickness of 100 microns. In the creation of this shell section, you also specify what pre-defined material to set the section to (in our case, we used aluminum for all our samples to keep consistency). After this section is created, each part must be assigned to the section it belongs to. In our case, we had only one part and only one section but it is possible to assign multi-part models to different sections with different material properties. Python commands are described below. At this point in the simulation, I have set simplepart = model.parts['simple-cc-test']
Command used to create a new homogeneous shell section: model.HomogeneousShellSection(idealization=NO_IDEALIZATION, integrationRule=SIMPSON, material='aluminum', name='shell', nodalThicknessField='', numIntPts=5, poissonDefinition=DEFAULT, preIntegrate=OFF, temperature=GRADIENT, thickness=thickness, thicknessField='', thicknessModulus=None, thicknessType=UNIFORM, useDensity=OFF) where thickness is set by command line.
Command used to assign a part to a section (in this situation a shell section): simplepart.SectionAssignment(offset=0.0, offsetField='', offsetType=MIDDLE_SURFACE, region=simplepart.sets['body'], sectionName='shell', thicknessAssignment=FROM_SECTION).
After section creation and assignment, each part must be placed in the assembly. For this step of the simulation, a coordinate system must first be defined so that the part can be placed in a geometric space in relation to other parts. In Abaqus, placing a part in this way is called creating an instance. Python commands for this are described below. I have set root = mdb.models['Model-1'].rootAssembly
Command used to create new coordinate system: root.DatumCsysByDefault(CARTESIAN).
Command used to create instance of part in coordinate system: root.Instance(dependent=OFF, name='simple-cc-test-1', part=simplepart).
After the part has been created and placed in an assembly, we can begin building the simulation to apply to the model. The first step in this is creating steps to apply forces at. For this portion of the simulation, I created a new step for each increment of force I wanted to apply. The number of increments and value of increment was taken from command line arguments. Steps were created in a loop based on these numbers. The python command for step creation and the loop used in my simulation are below.
Command used for step creation: model.StaticStep(description='apply force '+str(i), name='force-'+str(i), previous='Initial').
Loop used for creating step increments: for i in range(0,numsteps+1):
if(i == 0): .
model.StaticStep(description='apply force '+str(i), name='force-'+str(i), previous='Initial')
else:
model.StaticStep(description='apply force '+str(i), name='force-'+str(i), previous='force-'+str(i-1))
#adding backwards steps
for i in range(0, numsteps+1, 1):
if(i == 0):
model.StaticStep(description='release force '+str(-i), name='rforce-'+str(i), previous='force-'+str(numsteps))
else:
model.StaticStep(description='release force '+str(-i), name='rforce-'+str(i), previous='rforce-'+str(i-1))
After these steps are created, we can now set up the specific conditions of the simulation. I start by applying boundary conditions constraining the bottom of the model to stay stationary through the simulation. For this task, we prescribe the displacement of the bottom curves of the model to 0. To apply this condition through all of the steps, I use a for loop. The python command for boundary conditions constraining the bottom are described below. This same command with different inputs is used to constrain the front and back edges in the y-direction. It is important to note that these are not the only boundary conditions we are applying. We will apply periodic boundary conditions after meshing the part instance.
Command used for boundary condition creation: model.DisplacementBC(amplitude=UNSET, createStepName='force-'+str(i), distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None, name='constrainbottom-'+str(i), region=root.sets['bottomedge'], u1=0.0, u2=0.0, u3=0.0, ur1=UNSET, ur2=UNSET, ur3=UNSET) where i is used to denote which step we're applying a boundary condition in.
To observe crushing dynamics and self-collisions in the simulation, we decided to apply general contact properties. In Abaqus, this is an interaction that makes surfaces crush on impact with other surfaces in the model. A general contact interaction must first be created. Within this interaction, different interaction properties like tangential contact or normal contact can be defined. In this simulation, we chose to simulate normal or "hard" contact. The python commands defining this are described below.
Command used for first defining a contact property: mdb.models['Model-1'].ContactProperty('normal')
Command used for defining the properties of this contact (we chose 'HARD' in this simulation): mdb.models['Model-1'].interactionProperties['normal'].NormalBehavior(allowSeparation=ON, constraintEnforcementMethod=DEFAULT, pressureOverclosure=HARD).
Command used for first defining a type of interaction: mdb.models['Model-1'].ContactStd(createStepName='Initial', name='general-contact').
Command used for first defining what parts are in the interation: mdb.models['Model-1'].interactions['general-contact'].includedPairs.setValuesInStep(stepName='Initial', useAllstar=ON).
Command used for defining contact properties of part: mdb.models['Model-1'].interactions['general-contact'].contactPropertyAssignments.appendInStep(assignments=((GLOBAL, SELF, 'normal'), ), stepName='Initial')
Command used to assign surface thickness properties to general contact interation: mdb.models['Model-1'].interactions['general-contact'].surfaceThicknessAssignments.appendInStep(assignments=((mdb.models['Model-1'].rootAssembly.surfaces['all-faces'], ORIGINAL, 1.0), ), stepName='Initial')
Before we can apply loads and periodic boundary conditions, we must mesh the part instance. This is because the loads and boundary conditions are applied at mesh points rather than points defined on the initial model. To generate a mesh, the part instance must first be seeded. The size parameter in this command decides how dense the mesh is. It is important to consider this parameter especially when the scale of the model is changed. For this simulation, I always use size = 0.00059 because the models are around 5 mm (0.005 m) tall and 20 mm long. The python command for seeding and meshing are defined below.
Command to seed: root.seedPartInstance(deviationFactor=0.1, minSizeFactor=0.1, regions=( instance, ), size= 0.00059)
Command to mesh: root.generateMesh(regions=( instance, ))
After applying contact interactions and meshing, we apply loads at each step. My simulation code takes in the input force from command line, finds the number of nodes through the top of the model and finds how much force to apply per node to evenly distribute the load across the whole model. The python command for applying loads is below.
Command to apply loads: model.ConcentratedForce(cf3=-forceincrem, createStepName='force-'+str(i), distributionType=UNIFORM, field='', localCsys=None, name='Load-'+str(i), region=root.sets['topnodes'])
to be updated
The final step in simulation setup is job creation and submission. The python commands to do so are described below.
Command to create job: mdb.Job(atTime=None, contactPrint=OFF, description='apply force', echoPrint=OFF, explicitPrecision=SINGLE, getMemoryFromAnalysis=True, historyPrint=OFF, memory=90, memoryUnits=PERCENTAGE, model='Model-1', modelPrint=OFF, multiprocessingMode=DEFAULT, name='putforce', nodalOutputPrecision=SINGLE, numCpus=1, numGPUs=0, queue=None, resultsFormat=ODB, scratch='', type=ANALYSIS, userSubroutine='', waitHours=0, waitMinutes=0)
Command to submit job: mdb.jobs['putforce'].submit(consistencyChecking=OFF)
Command to wait for job completion before extracting values: mdb.jobs['putforce'].waitForCompletion()
to be updated
We decided to compare different thicknesses of material first. The model used for this comparison had an angle of 30 degrees and a curvature of 5 cm.
Fig. 1. left: Applied force in Newtons vs. Displacement in microns over different thicknesses; right: undeformed 30 degree, 5 cm curved-crease Miura-Ori fold
Fig. 2. Deformed models in simulation; left to right: 100 microns thick, 150 microns thick, 200 microns thick; top to bottom: after 10 steps of applying force; after 10 steps of releasing force
We also compared different radii of curvature. In the graphs shown below, each model has a fold angle of 45 degrees and we iterate over the following radii of curvature: 5 cm, 9 cm and 13 cm.
Below is the graph comparing these different radii and the deformed results from these simulations. It is worth noting that to get comparable graphs, we must do different simulations with different force increments. Currently, the displacement of one simulation eclipses that of the other simulation so we only see the other simulations as lines.
Fig. 3. From left to right: Effective stress vs. displacement over curvatures of 5 cm, 9 cm, and 13 cm; Applied force in Newtons vs. displacement in microns over curvatures of 5 cm, 9 cm and 13 cm.
Fig. 4. Deformed models applying and releasing force
Finally, we compared different fold angles. As mentioned before, the fold angles used for comparison are 30 degrees, 45 degrees and 60 degrees. The models are below.
Fig. 5. Undeformed models; left to right: 30 degrees, 45 degrees, and 60 degrees in fold angle
Below are the graphs and deformed models associated.
Fig. 6. Applied force in Newtons vs. displacement in microns over fold angles of 30 degrees, 45 degrees, and 60 degrees.
Fig. 7. Deformed models applying and releasing force
Implement contact in simulation so that crumpling is taken into account.
Work with corrugated models.
Programatically get deformed images from Abaqus output database.
In deformed images, get same scale factor across all simulations so visualizations are comparable
Try to piece together trends based on curvature and fold angle from data.