from abaqus import mdb from abaqusConstants import OFF import os as _sw_os _sw_before_jobs = {str(n): mdb.jobs[n] for n in mdb.jobs.keys()} _sw_written_jobs = [] def _sw_record_write(_job, *args, **kwargs): _answer = _job.writeInput(*args, **kwargs) _name = str(getattr(_job, 'name', '')) if not _name or _name not in mdb.jobs.keys(): _names = [str(n) for n in mdb.jobs.keys() if mdb.jobs[n] is _job] if len(_names) != 1: raise RuntimeError('Cannot identify the job which wrote the INP') _name = _names[0] _sw_written_jobs.append((_name, _sw_os.path.abspath(_name + '.inp'))) return _answer from abaqus import * from abaqusConstants import * import mesh import math base = 'BatteryOrphanImpact' modelName = base index = 1 while modelName in mdb.models.keys() or modelName + '_Job' in mdb.jobs.keys(): modelName = base + '_%02d' % index index += 1 model = mdb.Model(name=modelName, modelType=STANDARD_EXPLICIT) assembly = model.rootAssembly assembly.DatumCsysByDefault(CARTESIAN) P = dict(impactX=-300.0, impactY=-260.0, radius=15.0, mass=0.01, energy=150000.0, duration=0.006, friction=0.2, slip=0.005) speed = math.sqrt(2.0 * P['energy'] / P['mass']) materials = (('Steel', 7.85e-09, 210000.0, 0.3, ((400.0, 0.0), (560.0, 0.02), (650.0, 0.06), (735.0, 0.12), (800.0, 0.2), (840.0, 0.28))), ('CoolingAl', 2.73e-09, 69000.0, 0.33, ((50.0, 0.0), (90.0, 0.02), (130.0, 0.05), (160.0, 0.1), (180.0, 0.16), (190.0, 0.24))), ('HousingAl', 2.73e-09, 69000.0, 0.33, ((180.0, 0.0), (195.0, 0.005), (200.0, 0.012), (201.0, 0.04), (202.0, 0.07))), ('FrameAl', 2.7e-09, 70000.0, 0.33, ((200.0, 0.0), (240.0, 0.04), (280.0, 0.15))), ('PVC', 1.4e-09, 30.0, 0.35, ()), ('Glue', 1.5e-09, 20.0, 0.35, ()), ('Core', 1.9e-09, 50.0, 0.2, ()), ('Foam', 2e-10, 3.0, 0.05, ())) for (name, density, modulus, poisson, plastic) in materials: material = model.Material(name=name) material.Density(table=((density,),)) material.Elastic(table=((modulus, poisson),)) if plastic: material.Plastic(table=plastic) model.HomogeneousSolidSection(name=name + '_Section', material=name) def axis(a, b, divisions): return [a + (b - a) * float(i) / divisions for i in range(divisions + 1)] def refined_axis(a, b, center): return axis(a, center - 40.0, int(math.ceil((center - 40.0 - a) / 28.0)))[:-1] + axis(center - 40.0, center + 40.0, 20)[:-1] + axis(center + 40.0, b, int(math.ceil((b - center - 40.0) / 28.0))) def hex_part(name, xs, ys, zs, material, composite=False): part = model.Part(name=name, dimensionality=THREE_D, type=DEFORMABLE_BODY) (nx, ny, nz) = (len(xs), len(ys), len(zs)) for z in zs: for y in ys: for x in xs: part.Node(coordinates=(x, y, z)) nodes = part.nodes sides = dict(Bottom=[], Top=[], YMin=[], XMax=[], YMax=[], XMin=[]) coreLabels = [] housingLabels = [] label = 0 for k in range(nz - 1): for j in range(ny - 1): for i in range(nx - 1): n0 = (k * ny + j) * nx + i ids = (n0, n0 + 1, n0 + 1 + nx, n0 + nx, n0 + nx * ny, n0 + nx * ny + 1, n0 + nx * ny + 1 + nx, n0 + nx * ny + nx) part.Element(nodes=tuple((nodes[q] for q in ids)), elemShape=HEX8) label += 1 for (condition, side) in ((k == 0, 'Bottom'), (k == nz - 2, 'Top'), (j == 0, 'YMin'), (i == nx - 2, 'XMax'), (j == ny - 2, 'YMax'), (i == 0, 'XMin')): if condition: sides[side].append(label) if composite: if 0 < i < nx - 2 and 0 < j < ny - 2 and (0 < k < nz - 2): coreLabels.append(label) else: housingLabels.append(label) part.setElementType(regions=(part.elements[:],), elemTypes=(mesh.ElemType(elemCode=C3D8R, elemLibrary=EXPLICIT),)) if composite: groups = (('Core', coreLabels), ('HousingAl', housingLabels)) else: groups = ((material, list(range(1, label + 1))),) for (mat, labels) in groups: region = part.Set(name=mat + '_Elements', elements=part.elements.sequenceFromLabels(labels=tuple(labels))) part.SectionAssignment(region=region, sectionName=mat + '_Section') for (side, labels) in sides.items(): part.Set(name=side, elements=part.elements.sequenceFromLabels(labels=tuple(labels))) if not len(part.elements): raise RuntimeError('Empty mesh: ' + name) return part def place(part, name, shift=(0.0, 0.0, 0.0)): inst = assembly.Instance(name=name, part=part, dependent=ON) if shift != (0.0, 0.0, 0.0): assembly.translate(instanceList=(name,), vector=shift) return inst faceKeys = dict(Bottom='face1Elements', Top='face2Elements', YMin='face3Elements', XMax='face4Elements', YMax='face5Elements', XMin='face6Elements') def surface(inst, side, name): return assembly.Surface(name=name, **{faceKeys[side]: inst.sets[side].elements}) def bond(name, main, secondary): options = dict(name=name, positionToleranceMethod=SPECIFIED, positionTolerance=0.001, adjust=OFF, tieRotations=OFF) try: model.Tie(main=main, secondary=secondary, **options) except TypeError: model.Tie(master=main, slave=secondary, **options) xs = refined_axis(-800.0, 800.0, P['impactX']) ys = refined_axis(-550.0, 550.0, P['impactY']) items = {} for (name, bottom, top, divisions, material) in (('PVC', -1.0, 0.0, 1, 'PVC'), ('Guard', 0.0, 2.0, 3, 'Steel'), ('Foam', 2.0, 7.0, 3, 'Foam'), ('Cooling', 7.0, 9.0, 3, 'CoolingAl'), ('Glue', 9.0, 10.0, 1, 'Glue'), ('Cover', 130.0, 131.5, 2, 'FrameAl')): part = hex_part(name, xs, ys, axis(bottom, top, divisions), material) items[name] = place(part, name) bond('PVC_Guard', surface(items['Guard'], 'Bottom', 'GuardBottom'), surface(items['PVC'], 'Top', 'PVCTop')) bond('Glue_Cooling', surface(items['Cooling'], 'Top', 'CoolingTop'), surface(items['Glue'], 'Bottom', 'GlueBottom')) glueTop = surface(items['Glue'], 'Top', 'GlueTop') for (name, xx, yy) in (('LeftRail', axis(-820.0, -800.0, 2), axis(-550.0, 550.0, 40)), ('RightRail', axis(800.0, 820.0, 2), axis(-550.0, 550.0, 40)), ('FrontRail', axis(-820.0, 820.0, 60), axis(-570.0, -550.0, 2)), ('RearRail', axis(-820.0, 820.0, 60), axis(550.0, 570.0, 2))): part = hex_part(name, xx, yy, axis(0.0, 131.5, 8), 'FrameAl') items[name] = place(part, name) railSides = (('LeftRail', 'XMax', 'XMin'), ('RightRail', 'XMin', 'XMax'), ('FrontRail', 'YMax', 'YMin'), ('RearRail', 'YMin', 'YMax')) for (railName, innerSide, plateSide) in railSides: railSurface = surface(items[railName], innerSide, railName + '_Inner') for plateName in ('Guard', 'Cooling', 'Cover'): bond(plateName + '_' + railName, railSurface, surface(items[plateName], plateSide, plateName + '_' + plateSide)) for (endName, endSide, sideSide) in (('FrontRail', 'YMax', 'YMin'), ('RearRail', 'YMin', 'YMax')): endSurface = surface(items[endName], endSide, endName + '_CornerMain') for sideName in ('LeftRail', 'RightRail'): bond(endName + '_' + sideName, endSurface, surface(items[sideName], sideSide, sideName + '_' + sideSide)) cx = [-55.0] + axis(-54.0, 54.0, 10) + [55.0] cy = [-240.0] + axis(-239.0, 239.0, 24) + [240.0] cz = [0.0] + axis(1.0, 99.0, 8) + [100.0] cell = hex_part('Cell', cx, cy, cz, 'HousingAl', composite=True) for (row, y) in enumerate((-260.0, 260.0)): for column in range(12): name = 'Cell_%02d_%02d' % (row + 1, column + 1) inst = place(cell, name, (-660.0 + 120.0 * column, y, 10.0)) bond(name + '_Glue', glueTop, surface(inst, 'Bottom', name + '_Bottom')) assembly.regenerate() for railName in ('FrontRail', 'RearRail'): inst = items[railName] labels = [] for node in inst.nodes: (x, y, z) = node.coordinates if abs(z - 131.5) < 1e-06: if min((abs(x - q) for q in (-600.0, -200.0, 200.0, 600.0))) < 22.0: labels.append(node.label) if not labels: raise RuntimeError('Empty mounting set: ' + railName) region = assembly.Set(name=railName + '_Mounts', nodes=inst.nodes.sequenceFromLabels(labels=tuple(labels))) model.EncastreBC(name=railName + '_Fixed', createStepName='Initial', region=region) head = model.Part(name='Head', dimensionality=THREE_D, type=DISCRETE_RIGID_SURFACE) r = P['radius'] (nt, np) = (20, 48) points = [(0.0, 0.0, r)] for j in range(1, nt): theta = math.pi * j / nt for i in range(np): phi = 2.0 * math.pi * i / np points.append((r * math.sin(theta) * math.cos(phi), r * math.sin(theta) * math.sin(phi), r * math.cos(theta))) points.append((0.0, 0.0, -r)) for point in points: head.Node(coordinates=point) triangles = [] for i in range(np): triangles.append((0, 1 + i, 1 + (i + 1) % np)) for j in range(nt - 2): for i in range(np): q = 1 + j * np + i qn = 1 + j * np + (i + 1) % np triangles.extend(((q, q + np, qn + np), (q, qn + np, qn))) last = len(points) - 1 for i in range(np): q = 1 + (nt - 2) * np + i qn = 1 + (nt - 2) * np + (i + 1) % np triangles.append((q, last, qn)) nodes = head.nodes for ids in triangles: head.Element(nodes=tuple((nodes[q] for q in ids)), elemShape=TRI3) head.setElementType(regions=(head.elements[:],), elemTypes=(mesh.ElemType(elemCode=R3D3, elemLibrary=EXPLICIT),)) rp = head.ReferencePoint(point=(0.0, 0.0, 0.0)) head.Set(name='RP', referencePoints=(head.referencePoints[rp.id],)) head.Set(name='Body', elements=head.elements[:]) moment = 0.4 * P['mass'] * r * r head.engineeringFeatures.PointMassInertia(name='Mass', region=head.sets['RP'], mass=P['mass'], i11=moment, i22=moment, i33=moment) headInst = place(head, 'Head', (P['impactX'], P['impactY'], -r - 1.05)) model.RigidBody(name='RigidHead', refPointRegion=headInst.sets['RP'], bodyRegion=headInst.sets['Body']) model.DisplacementBC(name='Guide', createStepName='Initial', region=headInst.sets['RP'], u1=0.0, u2=0.0, u3=UNSET, ur1=0.0, ur2=0.0, ur3=0.0) model.Velocity(name='InitialVelocity', region=headInst.sets['RP'], velocity1=0.0, velocity2=0.0, velocity3=speed, omega=0.0) model.ExplicitDynamicsStep(name='Impact', previous='Initial', timePeriod=P['duration'], nlgeom=ON) prop = model.ContactProperty('Contact') prop.NormalBehavior(pressureOverclosure=HARD, allowSeparation=ON) prop.TangentialBehavior(formulation=PENALTY, table=((P['friction'],),), maximumElasticSlip=FRACTION, fraction=P['slip']) contact = model.ContactExp(name='GeneralContact', createStepName='Initial') contact.includedPairs.setValuesInStep(stepName='Initial', useAllstar=ON) contact.contactPropertyAssignments.appendInStep(stepName='Initial', assignments=((GLOBAL, SELF, 'Contact'),)) model.fieldOutputRequests['F-Output-1'].setValues(variables=('S', 'PEEQ', 'U', 'V', 'A'), numIntervals=100) model.historyOutputRequests['H-Output-1'].setValues(variables=PRESELECT, numIntervals=500) model.HistoryOutputRequest(name='HeadMotion', createStepName='Impact', variables=('U3', 'V3', 'A3'), region=headInst.sets['RP'], numIntervals=1000) assembly.regenerate() job = mdb.Job(name=modelName + '_Job', model=modelName, type=ANALYSIS, numCpus=4, numDomains=4) _sw_record_write(job, consistencyChecking=ON) print('Input written: ' + job.name + '.inp') print('SOUTHWEITION_JOB_PREPARED') if not _sw_written_jobs: _sw_candidates = [n for n in mdb.jobs.keys() if n not in _sw_before_jobs or mdb.jobs[n] is not _sw_before_jobs[n]] if len(_sw_candidates) != 1: raise RuntimeError('Cannot identify this script job: call writeInput() on exactly the intended job') _sw_record_write(mdb.jobs[_sw_candidates[0]], consistencyChecking=OFF) _sw_name, _sw_inp = _sw_written_jobs[-1] if not _sw_os.path.isfile(_sw_inp): raise RuntimeError('This script did not produce its INP: ' + _sw_inp) print('SOUTHWEITION_CWD ' + _sw_os.path.dirname(_sw_inp)) print('SOUTHWEITION_JOB ' + _sw_name) print('SOUTHWEITION_MODEL ' + str(getattr(mdb.jobs[_sw_name], 'model', ''))) print('SOUTHWEITION_INP ' + _sw_inp) print('SOUTHWEITION_JOB_PREPARED')