# -*- coding: utf-8 -*- # Run in Abaqus/CAE using File > Run Script. # Units: m, kg, s, Pa. # Solid slab reconstruction, not a validated exact reproduction. from abaqus import * from abaqusConstants import * import regionToolset import mesh import math TNT_MASS = 20.0 MESH_SIZE = 0.10 LX, LY = 6.0, 3.0 NSTOREY, STOREY_H = 3, 3.0 SLAB_T, COL_B = 0.10, 0.25 BEAM_B, BEAM_H = 0.25, 0.50 COVER = 0.04 STANDOFF = 2.0 ROOF_Z = NSTOREY * STOREY_H + SLAB_T SOURCE_X, SOURCE_Y = 1.5, 1.5 BASE_NAME = 'SolidSlab_%gkg' % TNT_MASS NAME = BASE_NAME counter = 1 while NAME in mdb.models.keys() or NAME in mdb.jobs.keys(): NAME = BASE_NAME + '_%d' % counter counter += 1 model = mdb.Model(name=NAME) assembly = model.rootAssembly assembly.DatumCsysByDefault(CARTESIAN) # Concrete: strength exponents interpreted as MPa. concrete = model.Material(name='Concrete') concrete.Density(table=((2400.0,),)) concrete.Elastic(table=((24.0e9, 0.18),)) concrete.ConcreteDamagedPlasticity( table=((30.0, 1.0, 1.16, 0.666, 0.10),)) cdp = concrete.concreteDamagedPlasticity cdp.ConcreteCompressionHardening(table=( (20.0e6, 0.0), (15.0e6, 0.0011), (12.0e6, 0.004))) cdp.ConcreteCompressionDamage(table=( (0.0, 0.0), (0.20, 0.0011), (0.24, 0.004))) cdp.ConcreteTensionStiffening(type=STRAIN, table=( (3.30e6, 0.0), (3.20e6, 0.003), (3.10e6, 0.005), (3.00e6, 0.007), (2.95e6, 0.010))) cdp.ConcreteTensionDamage(type=STRAIN, table=( (0.0, 0.0), (0.15, 0.003), (0.21, 0.005), (0.28, 0.007), (0.30, 0.010))) steel = model.Material(name='Steel') steel.Density(table=((7850.0,),)) steel.Elastic(table=((203.0e9, 0.33),)) steel.Plastic(hardening=JOHNSON_COOK, table=( (304.33e6, 422.007e6, 0.345, 0.87, 1800.0, 293.0),)) steel.plastic.RateDependent( type=JOHNSON_COOK, table=((0.0156, 1.0),)) # Assumed reference strain rate: 1/s. # No invented steel fracture-evolution parameter is supplied. model.HomogeneousSolidSection( name='ConcreteSection', material='Concrete') for diameter in (8, 10, 12, 16): model.TrussSection( name='D%d' % diameter, material='Steel', area=math.pi * (diameter / 1000.0)**2 / 4.0) concrete_instances = [] def add_box(name, dimensions, position): sketch = model.ConstrainedSketch(name='BoxSketch', sheetSize=30.0) sketch.rectangle( point1=(0.0, 0.0), point2=(dimensions[0], dimensions[1])) part = model.Part( name=name, dimensionality=THREE_D, type=DEFORMABLE_BODY) part.BaseSolidExtrude(sketch=sketch, depth=dimensions[2]) del model.sketches['BoxSketch'] instance = assembly.Instance( name=name + '_i', part=part, dependent=ON) assembly.translate(instanceList=(instance.name,), vector=position) concrete_instances.append(instance) for ix in range(3): for iy in range(2): add_box( 'Column_%d_%d' % (ix, iy), (COL_B, COL_B, NSTOREY * STOREY_H), (3.0 * ix - COL_B / 2.0, 3.0 * iy - COL_B / 2.0, 0.0)) for level in range(1, NSTOREY + 1): z = level * STOREY_H add_box('Slab_%d' % level, (LX, LY, SLAB_T), (0.0, 0.0, z)) for iy in range(2): add_box( 'BeamX_%d_%d' % (level, iy), (LX, BEAM_B, BEAM_H), (0.0, iy * 3.0 - BEAM_B / 2.0, z - BEAM_H)) for ix in range(3): add_box( 'BeamY_%d_%d' % (level, ix), (BEAM_B, LY, BEAM_H), (ix * 3.0 - BEAM_B / 2.0, 0.0, z - BEAM_H)) rc_instance = assembly.InstanceFromBooleanMerge( name='ConcreteBody', instances=tuple(concrete_instances), keepIntersections=ON, originalInstances=SUPPRESS, domain=GEOMETRY) rc = model.parts['ConcreteBody'] rc.Set(name='AllConcrete', cells=rc.cells[:]) rc.SectionAssignment( region=rc.sets['AllConcrete'], sectionName='ConcreteSection') # Shared topology replaces separate concrete component ties. rc.seedPart(size=MESH_SIZE, deviationFactor=0.1, minSizeFactor=0.1) rc.setMeshControls( regions=rc.cells[:], elemShape=TET, technique=FREE) rc.setElementType( regions=(rc.cells[:],), elemTypes=(mesh.ElemType( elemCode=C3D10M, elemLibrary=EXPLICIT),)) rc.generateMesh() if not len(rc.elements): raise RuntimeError('Concrete mesh is empty.') # Keep each diameter in a separate part to prevent section misassignment. bar_lines = dict((d, []) for d in (8, 10, 12, 16)) def add_bar(diameter, start, end): bar_lines[diameter].append((tuple(start), tuple(end))) def add_stirrup(points): for j in range(4): add_bar(8, points[j], points[(j + 1) % 4]) offset = COL_B / 2.0 - COVER for x in (0.0, 3.0, 6.0): for y in (0.0, 3.0): for dx in (-offset, offset): for dy in (-offset, offset): add_bar(16, (x+dx, y+dy, 0.02), (x+dx, y+dy, NSTOREY * STOREY_H)) for j in range(45): z = 0.10 + j * 0.20 add_stirrup(( (x-offset, y-offset, z), (x+offset, y-offset, z), (x+offset, y+offset, z), (x-offset, y+offset, z))) def beam_point(axis, centre, u, v, z): if axis == 0: return (u, centre + v, z) return (centre + v, u, z) for level in range(1, NSTOREY + 1): z = level * STOREY_H zb, zt = z - BEAM_H + COVER, z - COVER beam_axes = [(0, y, LX) for y in (0.0, 3.0)] beam_axes += [(1, x, LY) for x in (0.0, 3.0, 6.0)] for axis, centre, length in beam_axes: half = BEAM_B / 2.0 - COVER for v in (-half, 0.0, half): add_bar(16, beam_point(axis, centre, COVER, v, zb), beam_point(axis, centre, length-COVER, v, zb)) for v in (-half, half): add_bar(12, beam_point(axis, centre, COVER, v, zt), beam_point(axis, centre, length-COVER, v, zt)) for j in range(int(round(length / 0.20))): u = 0.10 + j * 0.20 add_stirrup(( beam_point(axis, centre, u, -half, zb), beam_point(axis, centre, u, half, zb), beam_point(axis, centre, u, half, zt), beam_point(axis, centre, u, -half, zt))) for j in range(int(round(LY / 0.20))): y = 0.10 + j * 0.20 add_bar(10, (COVER, y, z+COVER), (LX-COVER, y, z+COVER)) for j in range(int(round(LX / 0.20))): x = 0.10 + j * 0.20 add_bar(10, (x, COVER, z+COVER+0.010), (x, LY-COVER, z+COVER+0.010)) reinforcement_instances = [] for diameter in (8, 10, 12, 16): part = model.Part( name='Rebar%d' % diameter, dimensionality=THREE_D, type=DEFORMABLE_BODY) part.WirePolyLine( points=tuple(bar_lines[diameter]), mergeType=IMPRINT, meshable=ON) part.Set(name='AllBars', edges=part.edges[:]) part.SectionAssignment( region=part.sets['AllBars'], sectionName='D%d' % diameter) part.seedPart( size=MESH_SIZE, deviationFactor=0.1, minSizeFactor=0.1) part.setElementType( regions=(part.edges[:],), elemTypes=(mesh.ElemType(elemCode=T3D2, elemLibrary=EXPLICIT),)) part.generateMesh() if not len(part.elements): raise RuntimeError('Empty reinforcement mesh: D%d' % diameter) instance = assembly.Instance( name='Rebar%d_i' % diameter, part=part, dependent=ON) reinforcement_instances.append(instance) assembly.regenerate() for instance in reinforcement_instances: model.EmbeddedRegion( name='Embedded_' + instance.name, embeddedRegion=regionToolset.Region(edges=instance.edges[:]), hostRegion=regionToolset.Region(cells=rc_instance.cells[:]), weightFactorTolerance=1.e-6, absoluteTolerance=0.0, fractionalTolerance=0.05, toleranceMethod=BOTH) tol = 1.e-6 base_faces = rc_instance.faces.getByBoundingBox( xMin=-1.0, yMin=-1.0, zMin=-tol, xMax=LX+1.0, yMax=LY+1.0, zMax=tol) if not len(base_faces): raise RuntimeError('Column base faces not found.') model.EncastreBC( name='FixedBases', createStepName='Initial', region=regionToolset.Region(faces=base_faces)) model.ExplicitDynamicsStep( name='Blast', previous='Initial', timePeriod=0.08) # Correct keyword is numIntervals, not numberInterval. model.fieldOutputRequests['F-Output-1'].setValues( variables=('S', 'U', 'V', 'PEEQ', 'DAMAGET', 'DAMAGEC'), numIntervals=20) model.HistoryOutputRequest( name='Energy', createStepName='Blast', variables=('ALLIE', 'ALLKE', 'ALLAE', 'ETOTAL'), numIntervals=200) roof_faces = rc_instance.faces.getByBoundingBox( xMin=-tol, yMin=-tol, zMin=ROOF_Z-tol, xMax=LX+tol, yMax=LY+tol, zMax=ROOF_Z+tol) if not len(roof_faces): raise RuntimeError('Roof faces not found.') roof_surface = assembly.Surface(name='Roof', side1Faces=roof_faces) rp = assembly.ReferencePoint( point=(SOURCE_X, SOURCE_Y, ROOF_Z+STANDOFF)) source = assembly.Set( name='BlastSource', referencePoints=(assembly.referencePoints[rp.id],)) model.IncidentWaveProperty( name='BlastProperty', definition=SURFACE_BLAST, massTNT=TNT_MASS, massFactor=1.0, lengthFactor=1.0, timeFactor=1.0, pressureFactor=1.0) model.IncidentWave( name='CONWEP', createStepName='Blast', definition=CONWEP, sourcePoint=source, surface=roof_surface, interactionProperty='BlastProperty', detonationTime=0.0, magnitudeFactor=1.0) job = mdb.Job( name=NAME, model=NAME, type=ANALYSIS, description='Solid slab reconstruction with documented assumptions', memory=90, memoryUnits=PERCENTAGE, explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, multiprocessingMode=DEFAULT, numCpus=4, numDomains=4) job.writeInput() print('Model, mesh and job created: ' + NAME) print('Input written; job has not been submitted.') print('SOUTHWEITION_JOB_PREPARED') import os as _sw_os try: print('SOUTHWEITION_CWD ' + _sw_os.getcwd()) def _sw_job_model(_job): _model = getattr(_job, 'model', '') try: return str(_model.name) except Exception: return str(_model or '') _sw_models = list(mdb.models.keys()) _sw_latest_model = _sw_models[-1] if _sw_models else '' _sw_name = '' for _cand in reversed(list(mdb.jobs.keys())): if _sw_latest_model and _sw_job_model(mdb.jobs[_cand]) == _sw_latest_model: _sw_name = str(_cand) break if not _sw_name: _sw_names = list(mdb.jobs.keys()) if _sw_names: _sw_name = str(_sw_names[-1]) if _sw_name: try: mdb.jobs[_sw_name].setValues(numCpus=4, numDomains=4) except Exception: pass mdb.jobs[_sw_name].writeInput(consistencyChecking=OFF) print('SOUTHWEITION_JOB ' + str(_sw_name)) print('SOUTHWEITION_MODEL ' + str(_sw_latest_model)) print('SOUTHWEITION_INP ' + _sw_os.path.abspath(str(_sw_name) + '.inp')) print('SOUTHWEITION_JOB_PREPARED') except Exception: print('SOUTHWEITION_JOB_PREPARED')