# -*- coding: utf-8 -*- from abaqus import mdb from abaqusConstants import * import math import os # SI units: m, kg, s, N, Pa. # Paper: j.aej.2019.11.004; approximate six-story Case 3. # Direct INP construction avoids geometry/API unpacking failures. # The paper depicts a building, not stairs. # Assumptions: story height 3.3 m, CDP curves, simplified blast pulse. # Simplifications: no gravity preload, prescribed column removal, # slab/beam reinforcement or stirrups; perfect casing/core bond. # These assumptions prevent claiming an exact contour reproduction. stem = 'BlastFrame_Case3' name = stem idx = 1 while name in mdb.models or name in mdb.jobs: name = stem + '_%03d' % idx idx += 1 H = 3.3 NF = 6 TS = 0.15 TEND = 0.2 GRID = 1.10 VERTICAL_GRID = 0.45 xs = [0.0, 6.6, 10.8, 17.4] ys = [0.0, 8.4, 16.8, 21.0, 29.6, 38.0] layout = [ ['C2', 'C3', 'C4', 'C1'], ['C7', 'C7', 'C8', 'C6'], ['C7', None, 'C8', 'C6'], ['C7', None, 'C7', 'C5'], ['C7', 'C7', 'C7', 'C6'], ['C2', 'C4', 'C3', 'C2']] dims = {'C1': [0.4, 0.5], 'C2': [0.5, 0.5], 'C3': [0.5, 0.7], 'C4': [0.5, 0.9], 'C5': [0.5, 1.0], 'C6': [0.5, 1.1], 'C7': [0.5, 1.5], 'C8': [0.5, 1.8]} counts = {'C1': 12, 'C2': 12, 'C3': 14, 'C4': 16, 'C5': 18, 'C6': 20, 'C7': 30, 'C8': 34} columns = [] vx = [] vy = [] for j in range(len(ys)): for i in range(len(xs)): typ = layout[j][i] if typ is None: continue b = dims[typ][0] d = dims[typ][1] thickness = 0.0 if i == 0: b = 0.4 if typ == 'C2' else 1.2 d = 0.4 thickness = 0.0046 if typ == 'C2' else 0.0065 c = dict(x=xs[i], y=ys[j], typ=typ, t=thickness, x0=xs[i]-b/2, x1=xs[i]+b/2, y0=ys[j]-d/2, y1=ys[j]+d/2) columns.append(c) vx.extend([c['x0'], c['x'], c['x1']]) vy.extend([c['y0'], c['y'], c['y1']]) for x in xs: vx.extend([x-0.15, x+0.15]) for y in ys: vy.extend([y-0.4, y+0.4]) def make_axis(values, size): values = sorted(set(round(float(v), 6) for v in values)) result = [values[0]] for q in range(len(values)-1): lo = values[q] hi = values[q+1] count = max(1, int(math.ceil((hi-lo)/size))) for k in range(1, count+1): result.append(round(lo+(hi-lo)*k/count, 6)) return result X = make_axis(vx, GRID) Y = make_axis(vy, GRID) Z = [0.0] for story in range(NF): top = (story+1)*H Z.extend(make_axis([story*H, top-0.70, top-0.37, top-TS, top], VERTICAL_GRID)[1:]) node_ids = {} nodes = [] solids = [] skins = [] rebars = [] fixed = [] monitor = {} loads = {} skinsets = {'T46': [], 'T65': []} barsets = dict((key, []) for key in counts) def get_node(i, j, k): key = (i, j, k) if key not in node_ids: label = len(nodes)+1 node_ids[key] = label nodes.append([X[i], Y[j], Z[k]]) if k == 0: fixed.append(label) return node_ids[key] # Conforming grid: slabs, beams and columns share nodes. plan = {} for i in range(len(X)-1): x = (X[i]+X[i+1])/2 for j in range(len(Y)-1): y = (Y[j]+Y[j+1])/2 owner = -1 for q in range(len(columns)): c = columns[q] if c['x0'] < x < c['x1'] and c['y0'] < y < c['y1']: owner = q break db = any(abs(x-v) < 0.150001 for v in xs) hb = any(abs(y-v) < 0.400001 for v in ys) plan[(i, j)] = [owner, db, hb] for k in range(len(Z)-1): zm = (Z[k]+Z[k+1])/2 story = min(NF-1, int(zm/H)) depth = (story+1)*H-zm for i in range(len(X)-1): for j in range(len(Y)-1): info = plan[(i, j)] if not (info[0] >= 0 or depth < TS+1e-7 or (info[1] and depth < 0.70) or (info[2] and depth < 0.37)): continue points = [(i,j,k), (i+1,j,k), (i+1,j+1,k), (i,j+1,k), (i,j,k+1), (i+1,j,k+1), (i+1,j+1,k+1), (i,j+1,k+1)] solids.append([get_node(*pt) for pt in points]) if info[0] >= 0: c = columns[info[0]] if c['t'] > 0 and abs(X[i]-c['x0']) < 1e-6: key = 'BL_%d_%d' % (info[0], story) if key not in loads: loads[key] = dict(labels=[], col=info[0], story=story) loads[key]['labels'].append(len(solids)) xi = dict((round(X[i],6), i) for i in range(len(X))) yi = dict((round(Y[i],6), i) for i in range(len(Y))) for c in columns: i0 = xi[round(c['x0'],6)] i1 = xi[round(c['x1'],6)] j0 = yi[round(c['y0'],6)] j1 = yi[round(c['y1'],6)] for k in range(len(Z)-1): if c['t'] > 0: edges = [] for j in range(j0, j1): edges.extend([[(i0,j),(i0,j+1)], [(i1,j+1),(i1,j)]]) for i in range(i0, i1): edges.extend([[(i+1,j0),(i,j0)], [(i,j1),(i+1,j1)]]) for edge in edges: u = edge[0] v = edge[1] skins.append([get_node(u[0],u[1],k), get_node(v[0],v[1],k), get_node(v[0],v[1],k+1), get_node(u[0],u[1],k+1)]) key = 'T46' if c['t'] < 0.005 else 'T65' skinsets[key].append(len(solids)+len(skins)) for pt in [(i0,j0), (i1,j0), (i1,j1), (i0,j1)]: rebars.append([c['typ'], get_node(pt[0],pt[1],k), get_node(pt[0],pt[1],k+1)]) if abs(c['x']) < 1e-6 and abs(c['y']-16.8) < 1e-6: for level in range(1, NF+1): kz = min(range(len(Z)), key=lambda q: abs(Z[q]-level*H)) monitor[level] = get_node(i0,j0,kz) path = os.path.abspath(name+'_source.inp') with open(path, 'w') as f: def out(text): f.write(text+'\n') def write_set(kind, key, labels): out('*'+kind+', '+kind+'='+key) for start in range(0, len(labels), 16): out(', '.join(str(v) for v in labels[start:start+16])) out('*HEADING\nApproximate Case 3; conforming orphan mesh') out('*PART, NAME=FRAME\n*NODE') for q in range(len(nodes)): out('%d, %.8g, %.8g, %.8g' % ( q+1, nodes[q][0], nodes[q][1], nodes[q][2])) for spec in [['C3D8R','CONCRETE',solids,0], ['S4R','CASING',skins,len(solids)]]: out('*ELEMENT, TYPE='+spec[0]+', ELSET='+spec[1]) for q in range(len(spec[2])): out(', '.join(str(v) for v in [spec[3]+q+1]+spec[2][q])) out('*ELEMENT, TYPE=T3D2, ELSET=REINFORCEMENT') offset = len(solids)+len(skins) for q in range(len(rebars)): row = rebars[q] label = offset+q+1 out('%d, %d, %d' % (label,row[1],row[2])) barsets[row[0]].append(label) write_set('NSET','FIXED',fixed) for level in sorted(monitor): write_set('NSET','MON_%d' % level,[monitor[level]]) for key in sorted(loads): write_set('ELSET',key,loads[key]['labels']) out('*SOLID SECTION, ELSET=CONCRETE, MATERIAL=CONC\n,') for key in sorted(skinsets): write_set('ELSET',key,skinsets[key]) t = 0.0046 if key == 'T46' else 0.0065 out('*SHELL SECTION, ELSET='+key+', MATERIAL=STEEL') out('%.8g, 5' % t) for key in sorted(barsets): if barsets[key]: write_set('ELSET','RB_'+key,barsets[key]) out('*SOLID SECTION, ELSET=RB_'+key+', MATERIAL=STEEL') out('%.10g' % (counts[key]*math.pi*0.016**2/16.0)) out('*END PART\n*ASSEMBLY, NAME=ASSEMBLY') out('*INSTANCE, NAME=FRAME-1, PART=FRAME\n*END INSTANCE\n*END ASSEMBLY') out('*MATERIAL, NAME=CONC\n*DENSITY\n2400.\n*ELASTIC\n3.0E10, 0.2') out('*CONCRETE DAMAGED PLASTICITY\n35., 0.1, 1.16, 0.667, 0.') out('*CONCRETE COMPRESSION HARDENING') out('18.E6,0.\n36.E6,0.0003\n45.E6,0.0005') out('30.E6,0.0016\n9.E6,0.0037\n1.5E6,0.008') out('*CONCRETE COMPRESSION DAMAGE') out('0.,0.\n0.,0.0003\n0.,0.0005\n0.35,0.0016\n0.75,0.0037\n0.95,0.008') out('*CONCRETE TENSION STIFFENING, TYPE=DISPLACEMENT') out('6.E6,0.\n3.E6,0.00005\n0.6E6,0.0002\n0.06E6,0.0005') out('*CONCRETE TENSION DAMAGE, TYPE=DISPLACEMENT') out('0.,0.\n0.45,0.00005\n0.85,0.0002\n0.98,0.0005') out('*MATERIAL, NAME=STEEL\n*DENSITY\n7850.\n*ELASTIC\n2.E11,0.3') out('*PLASTIC\n630.E6,0.\n750.E6,0.01\n960.E6,0.04\n1050.E6,0.08') for key in sorted(loads): data = loads[key] c = columns[data['col']] z = (data['story']+0.5)*H r = math.sqrt((c['x']+5)**2+(c['y']-19)**2+(z-1)**2) scaled = r/10.0 ps = max(0.,1772/scaled**3-114/scaled**2+108/scaled)*1000 data['peak'] = 2*ps*(7*101325+4*ps)/(7*101325+ps)*(5/r)**2 ta = r/1800. td = 0.008+0.0003*r out('*AMPLITUDE, NAME=A_'+key+', TIME=STEP TIME') out('0.,0., %.9g,0., %.9g,1., %.9g,0.' % (ta,ta+1e-5,ta+td)) out('%.9g,0.' % TEND) out('*BOUNDARY\nFRAME-1.FIXED,1,3') # Rotational constraints apply only to casing nodes at the base. rotation_nodes = sorted(set(v for row in skins for v in row if abs(nodes[v-1][2]) < 1e-8)) for label in rotation_nodes: out('FRAME-1.%d,4,6' % label) out('*STEP, NAME=Blast, NLGEOM=YES\n*DYNAMIC, EXPLICIT\n, %.9g' % TEND) for key in sorted(loads): out('*DLOAD, AMPLITUDE=A_'+key) out('FRAME-1.%s, P6, %.9g' % (key,loads[key]['peak'])) out('*OUTPUT, FIELD, NUMBER INTERVAL=40\n*NODE OUTPUT\nU,V') out('*ELEMENT OUTPUT, ELSET=FRAME-1.CONCRETE\nS,PEEQ,DAMAGET,DAMAGEC') out('*ELEMENT OUTPUT, ELSET=FRAME-1.CASING\nS,PEEQ') out('*OUTPUT, HISTORY, TIME INTERVAL=0.001') out('*ENERGY OUTPUT\nALLIE,ALLKE,ALLAE,ALLVD,ETOTAL') for level in sorted(monitor): out('*NODE OUTPUT, NSET=FRAME-1.MON_%d\nU1,U2,U3' % level) out('*END STEP') model = mdb.ModelFromInputFile(name=name, inputFileName=path) job = mdb.Job(name=name, model=name, type=ANALYSIS, numCpus=4, numDomains=4, parallelizationMethodExplicit=DOMAIN, explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE) job.writeInput(consistencyChecking=ON) print('Model and mesh imported: '+name) print('Job input written: '+name+'.inp') 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')