# -*- coding: utf-8 -*- # Run inside Abaqus/CAE: File -> Run Script. # Units: mm, N, s, tonne, MPa. # Reference: Kumar et al., Engineering Structures 206 (2020) 110125. # Slab-11: 1000x1000x100 mm; TNT 0.85 kg; standoff 100 mm. # REDUCED APPROXIMATION, NOT an exact reproduction: # Paper: HJC, 1.5-5 mm mesh, 1,685,625 concrete elements/quarter. # Here: CDP, graded Cartesian C3D8R mesh, 40,000 concrete elements. # No HJC SDV_DAMAGE, concrete erosion, or fabricated contour results. # CDP softening, steel curve interpolation and support elasticity # are modeling assumptions requiring calibration. from abaqus import mdb, session from abaqusConstants import * import regionToolset import material import section import interaction import step import load import mesh import job import os import math import tempfile base = 'Slab11_CDP_Light' model_name = base counter = 1 while model_name in mdb.models or model_name in mdb.jobs: model_name = '%s_%02d' % (base, counter) counter += 1 # Quarter slab: x,y=0..500; bottom z=0; loaded face z=100. # Opposite full-model supports are at x=+/-400..500. # Figure 4 is interpreted as 150 mm overall depth: # slab 100 mm + support 50 mm. Support width is 100 mm. nodes = [] groups = {} nsets = {} faces = {} element_count = [0] def node(x, y, z): label = len(nodes) + 1 nodes.append((label, float(x), float(y), float(z))) return label def element(group, conn): element_count[0] += 1 label = element_count[0] groups.setdefault(group, []).append((label, tuple(conn))) return label def linspace(a, b, count): return [a + (b-a)*i/float(count) for i in range(count+1)] def brick_grid(xs, ys, zs, group): nx, ny, nz = len(xs)-1, len(ys)-1, len(zs)-1 ids = {} for k, z in enumerate(zs): for j, y in enumerate(ys): for i, x in enumerate(xs): ids[i, j, k] = node(x, y, z) bottom, top = [], [] for k in range(nz): for j in range(ny): for i in range(nx): conn = [ids[i,j,k], ids[i+1,j,k], ids[i+1,j+1,k], ids[i,j+1,k], ids[i,j,k+1], ids[i+1,j,k+1], ids[i+1,j+1,k+1], ids[i,j+1,k+1]] lab = element(group, conn) if k == 0: bottom.append(lab) if k == nz-1: top.append(lab) return bottom, top # Cell widths increase smoothly from 5 to 15 mm. xy = [0.0] for i in range(50): xy.append(xy[-1] + 5.0 + 10.0*i/49.0) xy[-1] = 500.0 faces['SLAB_BOTTOM'], faces['SLAB_TOP'] = brick_grid( xy, xy, linspace(0.0, 100.0, 16), 'CONCRETE') _, faces['SUPPORT_TOP'] = brick_grid( linspace(400.0, 500.0, 5), linspace(0.0, 500.0, 25), linspace(-50.0, 0.0, 4), 'SUPPORT') # Single mid-depth reinforcement mat: diameter 10, spacing 100. # Bar ends at 490 mm: assumed 10 mm end cover. # Bars lying on a symmetry plane receive half cross-sectional area. # Orthogonal bars use independent nodes and are embedded in concrete. for direction in (0, 1): for station in (0.0, 100.0, 200.0, 300.0, 400.0): chain = [] for q in linspace(0.0, 490.0, 49): x, y = (q, station) if direction == 0 else (station, q) chain.append(node(x, y, 50.0)) group = 'REBAR_HALF' if station == 0.0 else 'REBAR_FULL' for a, b in zip(chain[:-1], chain[1:]): element(group, (a, b)) nsets['XSYM'] = [n[0] for n in nodes if abs(n[1]) < 1.e-8] nsets['YSYM'] = [n[0] for n in nodes if abs(n[2]) < 1.e-8] nsets['FIXED'] = [n[0] for n in nodes if abs(n[3]+50.0) < 1.e-8] nsets['CENTER_BACK'] = [ n[0] for n in nodes if abs(n[1])+abs(n[2])+abs(n[3]) < 1.e-8] # Importing a generated orphan mesh avoids tens of thousands of # individual CAE node/element API calls. handle, mesh_file = tempfile.mkstemp(suffix='.inp', prefix='slab11_') os.close(handle) def write_ids(stream, labels): for start in range(0, len(labels), 16): stream.write(', '.join(str(v) for v in labels[start:start+16])+'\n') with open(mesh_file, 'w') as stream: stream.write('*HEADING\nReduced Slab-11 mesh\n*PART, NAME=BODY\n*NODE\n') for n in nodes: stream.write('%d, %.12g, %.12g, %.12g\n' % n) for group in ('CONCRETE', 'SUPPORT', 'REBAR_FULL', 'REBAR_HALF'): kind = 'T3D2' if group.startswith('REBAR') else 'C3D8R' stream.write('*ELEMENT, TYPE=%s, ELSET=%s\n' % (kind, group)) for label, conn in groups[group]: stream.write(', '.join(str(v) for v in (label,)+conn)+'\n') stream.write('*ELSET, ELSET=REBAR\nREBAR_FULL, REBAR_HALF\n') for name, labels in nsets.items(): stream.write('*NSET, NSET=%s\n' % name) write_ids(stream, labels) stream.write('*END PART\n*ASSEMBLY, NAME=ASSEMBLY\n') stream.write('*INSTANCE, NAME=BODY-1, PART=BODY\n*END INSTANCE\n') stream.write('*END ASSEMBLY\n') m = mdb.ModelFromInputFile(name=model_name, inputFileName=mesh_file) os.remove(mesh_file) p = m.parts['BODY'] # Concrete density and shear modulus follow Table 2. # Poisson ratio 0.20 and all CDP parameters below are approximations. nu = 0.20 young = 2.0*14648.0*(1.0+nu) concrete = m.Material(name='Concrete_CDP_Approx') concrete.Density(table=((2.430e-9,),)) concrete.Elastic(table=((young, nu),)) concrete.ConcreteDamagedPlasticity( table=((35.0, 0.1, 1.16, 0.667, 0.0),)) cdp = concrete.concreteDamagedPlasticity cdp.ConcreteCompressionHardening(table=( (18.0, 0.0), (35.0, 0.0005), (47.0, 0.0015), (35.0, 0.004), (12.0, 0.010), (2.0, 0.025))) cdp.ConcreteCompressionDamage(table=( (0.0, 0.0), (0.0, 0.0005), (0.0, 0.0015), (0.15, 0.004), (0.65, 0.010), (0.90, 0.025))) # Stress versus cracking displacement gives mesh-regularized softening. cdp.ConcreteTensionStiffening(type=DISPLACEMENT, table=( (3.80, 0.0), (2.0, 0.02), (0.50, 0.08), (0.038, 0.15))) cdp.ConcreteTensionDamage(type=DISPLACEMENT, table=( (0.0, 0.0), (0.45, 0.02), (0.85, 0.08), (0.98, 0.15))) steel = m.Material(name='RebarSteel') steel.Density(table=((7.85e-9,),)) steel.Elastic(table=((200000.0, 0.30),)) # Yield 609 MPa and ultimate engineering stress 745 MPa from paper. # Intermediate engineering curve points are approximate readings. eng_curve = ((609.0, 609.0/200000.0), (645.0, 0.025), (700.0, 0.060), (735.0, 0.100), (745.0, 0.120)) plastic = [(609.0, 0.0)] for stress, strain in eng_curve[1:]: true_stress = stress*(1.0+strain) plastic.append((true_stress, math.log(1.0+strain)-true_stress/200000.0)) steel.Plastic(table=tuple(plastic)) steel.DuctileDamageInitiation(table=((0.13, 0.333, 0.025),)) steel.ductileDamageInitiation.DamageEvolution( type=ENERGY, table=((3582.0,),)) support = m.Material(name='SupportElastic_Assumed') support.Density(table=((7.85e-9,),)) support.Elastic(table=((200000.0, 0.30),)) for name, mat in (('ConcreteSection', 'Concrete_CDP_Approx'), ('SupportSection', 'SupportElastic_Assumed')): m.HomogeneousSolidSection(name=name, material=mat, thickness=None) area = math.pi*10.0**2/4.0 for name, fraction_area in (('BarFull', 1.0), ('BarHalf', 0.5)): m.TrussSection(name=name, material='RebarSteel', area=area*fraction_area) for group, sec in (('CONCRETE', 'ConcreteSection'), ('SUPPORT', 'SupportSection'), ('REBAR_FULL', 'BarFull'), ('REBAR_HALF', 'BarHalf')): p.SectionAssignment(region=p.sets[group], sectionName=sec) for name, labels in faces.items(): elems = p.elements.sequenceFromLabels(labels=tuple(labels)) if name == 'SLAB_BOTTOM': p.Surface(name=name, face1Elements=elems) else: p.Surface(name=name, face2Elements=elems) a = m.rootAssembly a.regenerate() inst = a.instances['BODY-1'] m.EmbeddedRegion(name='EmbeddedRebar', embeddedRegion=inst.sets['REBAR'], hostRegion=inst.sets['CONCRETE']) for name, axis in (('XSYM', 1), ('YSYM', 2)): args = {'u%d' % axis: 0.0} m.DisplacementBC(name=name, createStepName='Initial', region=inst.sets[name], **args) m.DisplacementBC(name='SupportBase', createStepName='Initial', region=inst.sets['FIXED'], u1=0.0, u2=0.0, u3=0.0) m.ExplicitDynamicsStep(name='Blast', previous='Initial', timePeriod=0.008, nlgeom=ON) contact = m.ContactProperty('FrictionlessHard') contact.NormalBehavior(pressureOverclosure=HARD, allowSeparation=ON) # Positional surfaces accommodate master/slave -> main/secondary renaming. m.SurfaceToSurfaceContactExp( 'SupportContact', 'Initial', inst.surfaces['SUPPORT_TOP'], inst.surfaces['SLAB_BOTTOM'], FINITE, 'FrictionlessHard', mechanicalConstraint=KINEMATIC) source = a.ReferencePoint(point=(0.0, 0.0, 200.0)) source_region = regionToolset.Region( referencePoints=(a.referencePoints[source.id],)) m.IncidentWaveProperty(name='TNT085', definition=AIR_BLAST, massTNT=0.85, massFactor=1.0, lengthFactor=0.001, timeFactor=1.0, pressureFactor=1.e6) # TNT mass is the FULL charge: do not divide by four. # standoffPoint is unused for CONWEP; supplied for API compatibility. m.IncidentWave(name='AirBlast', createStepName='Blast', sourcePoint=source_region, standoffPoint=source_region, surface=inst.surfaces['SLAB_TOP'], interactionProperty='TNT085', definition=CONWEP, detonationTime=0.0, magnitudeFactor=1.0) for key in list(m.fieldOutputRequests.keys()): del m.fieldOutputRequests[key] for key in list(m.historyOutputRequests.keys()): del m.historyOutputRequests[key] m.FieldOutputRequest(name='Mechanical', createStepName='Blast', variables=('U', 'V', 'S', 'LE'), numIntervals=40) m.FieldOutputRequest(name='ConcreteDamage', createStepName='Blast', region=inst.sets['CONCRETE'], variables=('DAMAGET', 'DAMAGEC', 'PEEQ'), numIntervals=40) m.HistoryOutputRequest(name='Energy', createStepName='Blast', variables=('ALLKE', 'ALLIE', 'ALLAE', 'ALLVD', 'ETOTAL'), numIntervals=160) m.HistoryOutputRequest(name='CenterResponse', createStepName='Blast', region=inst.sets['CENTER_BACK'], variables=('U3', 'V3'), numIntervals=160) # No mass scaling. Check energy balance and mesh sensitivity after solving. # Post-solve in CAE: plot DAMAGET/DAMAGEC, range 0..1; inspect both faces. # These fields are NOT the paper's HJC SDV_DAMAGE or crack geometry. job = mdb.Job(name=model_name, model=model_name, type=ANALYSIS, numCpus=4, numDomains=4, memory=80, memoryUnits=PERCENTAGE, explicitPrecision=DOUBLE, nodalOutputPrecision=SINGLE) job.writeInput(consistencyChecking=ON) mdb.saveAs(pathName=os.path.abspath(model_name+'.cae')) if len(session.viewports): viewport = session.viewports[list(session.viewports.keys())[0]] viewport.setValues(displayedObject=a) viewport.assemblyDisplay.setValues(mesh=ON) viewport.view.setValues(session.views['Iso']) viewport.view.fitView() print('Created: %s; concrete elements: %d' % (model_name, len(groups['CONCRETE']))) print('Input written. Submit the existing job using CAE Job Manager.') 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')