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 import os import time tag = time.strftime('%Y%m%d_%H%M%S') modelName = 'BuriedPipe_' + tag jobName = 'BuriedPipeJob_' + tag model = mdb.Model(name=modelName) assembly = model.rootAssembly assembly.DatumCsysByDefault(CARTESIAN) D = 1219.0 THICKNESS = 26.4 LENGTH = 70.0 * D WIDTH = 11.0 * D HEIGHT = 9.0 * D COVER = 2000.0 YC = -COVER - 0.5 * D RO = 0.5 * D RM = RO - 0.5 * THICKNESS PRESSURE = 12.0 SLIP = 200.0 E_STEEL = 210000.0 NC = 24 RADIAL = [0.0, 0.035, 0.12, 0.3, 0.58, 1.0] NR = len(RADIAL) - 1 steel = model.Material(name='X80') steel.Density(table=((7.85e-09,),)) steel.Elastic(table=((E_STEEL, 0.3),)) steel.Plastic(table=((553.0, 0.0), (621.0, 0.03 - 621.0 / E_STEEL), (621.0, 0.3))) soil = model.Material(name='Soil') soil.Density(table=((1.9e-09,),)) soil.Elastic(table=((40.0, 0.35),)) soil.MohrCoulombPlasticity(table=((22.0, 5.0),)) soil.mohrCoulombPlasticity.MohrCoulombHardening(table=((0.04, 0.0), (0.04, 0.5))) model.HomogeneousSolidSection(name='SoilSection', material='Soil') model.HomogeneousShellSection(name='PipeSection', material='X80', thickness=THICKNESS, integrationRule=SIMPSON, numIntPts=5) def linspace(start, stop, count): values = [] for i in range(count + 1): values.append(start + (stop - start) * float(i) / count) return values def make_node_set(part, name, labels): values = list(labels) if len(values) == 0: raise RuntimeError('Empty node set: ' + name) selected = part.nodes.sequenceFromLabels(labels=values) part.Set(name=name, nodes=selected) def element_sequence(part, labels): return part.elements.sequenceFromLabels(labels=list(labels)) def add_element(part, indices, shape, label): selected = [] for index in indices: selected.append(part.nodes[index]) part.Element(nodes=selected, elemShape=shape, label=label) def assign_section(part, name, code): elements = part.elements[:] part.Set(name='ALL', elements=elements) elementType = mesh.ElemType(elemCode=code, elemLibrary=STANDARD) part.setElementType(regions=(elements,), elemTypes=(elementType,)) part.SectionAssignment(region=part.sets['ALL'], sectionName=name) corners = [(WIDTH / 2.0, 0.0), (-WIDTH / 2.0, 0.0), (-WIDTH / 2.0, -HEIGHT), (WIDTH / 2.0, -HEIGHT)] cornerAngles = [] for (x, y) in corners: cornerAngles.append(math.atan2(y - YC, x) % (2.0 * math.pi)) cornerAngles.sort() angles = [] for i in range(4): start = cornerAngles[i] stop = cornerAngles[(i + 1) % 4] if stop <= start: stop += 2.0 * math.pi for j in range(6): angles.append(start + (stop - start) * j / 6.0) crossSection = [] for fractionValue in RADIAL: ring = [] for theta in angles: cx = math.cos(theta) sy = math.sin(theta) tx = WIDTH / (2.0 * abs(cx)) if abs(cx) > 1e-12 else 1e+30 if sy > 1e-12: ty = -YC / sy elif sy < -1e-12: ty = (-HEIGHT - YC) / sy else: ty = 1e+30 outerRadius = tx if tx < ty else ty radius = RO + fractionValue * (outerRadius - RO) ring.append((radius * cx, YC + radius * sy)) crossSection.append(ring) distances = linspace(0.0, 4.0 * D, 12) distances.extend(linspace(4.0 * D, 12.0 * D, 8)[1:]) distances.extend(linspace(12.0 * D, LENGTH / 2.0, 12)[1:]) fractions = [] for distance in distances: fractions.append(distance / (LENGTH / 2.0)) NZ = len(fractions) - 1 soilInstances = [] def soil_index(k, radialIndex, angularIndex): return (k * (NR + 1) + radialIndex) * NC + angularIndex % NC for sideIndex in range(2): moving = sideIndex == 1 partName = 'MovingSoil' if moving else 'FixedSoil' part = model.Part(name=partName, dimensionality=THREE_D, type=DEFORMABLE_BODY) stations = list(fractions) if not moving: stations.reverse() (bottomNodes, sideNodes, endNodes) = ([], [], []) for (k, fractionValue) in enumerate(stations): for ir in range(NR + 1): for j in range(NC): (x, y) = crossSection[ir][j] faultZ = -(y - YC) endZ = LENGTH / 2.0 if moving else -LENGTH / 2.0 z = (1.0 - fractionValue) * faultZ + fractionValue * endZ label = soil_index(k, ir, j) + 1 part.Node(coordinates=(x, y, z), label=label) atBottom = abs(y + HEIGHT) < 0.001 if atBottom: bottomNodes.append(label) if abs(abs(x) - WIDTH / 2.0) < 0.001 and (not atBottom): sideNodes.append(label) if abs(fractionValue - 1.0) < 1e-10 and (not atBottom): endNodes.append(label) (holeElements, faultElements) = ([], []) for k in range(NZ): for ir in range(NR): for j in range(NC): indices = [soil_index(k, ir, j), soil_index(k, ir + 1, j), soil_index(k, ir + 1, j + 1), soil_index(k, ir, j + 1), soil_index(k + 1, ir, j), soil_index(k + 1, ir + 1, j), soil_index(k + 1, ir + 1, j + 1), soil_index(k + 1, ir, j + 1)] label = (k * NR + ir) * NC + j + 1 add_element(part, indices, HEX8, label) if ir == 0: holeElements.append(label) if moving and k == 0 or (not moving and k == NZ - 1): faultElements.append(label) assign_section(part, 'SoilSection', C3D8R) for (name, labels) in [('BOTTOM', bottomNodes), ('SIDE', sideNodes), ('END', endNodes)]: make_node_set(part, name, labels) part.Surface(name='HOLE', face6Elements=element_sequence(part, holeElements)) faultSequence = element_sequence(part, faultElements) if moving: part.Surface(name='FAULT', face1Elements=faultSequence) else: part.Surface(name='FAULT', face2Elements=faultSequence) instance = assembly.Instance(name=partName + '-1', part=part, dependent=ON) soilInstances.append(instance) part = model.Part(name='Pipe', dimensionality=THREE_D, type=DEFORMABLE_BODY) centralLength = 6.0 * D nFar = int(math.ceil((LENGTH / 2.0 - centralLength) / (0.3 * D))) nCenter = int(math.ceil(2.0 * centralLength / 100.0)) zValues = linspace(-LENGTH / 2.0, -centralLength, nFar) zValues.extend(linspace(-centralLength, centralLength, nCenter)[1:]) zValues.extend(linspace(centralLength, LENGTH / 2.0, nFar)[1:]) for (k, z) in enumerate(zValues): for (j, theta) in enumerate(angles): part.Node(coordinates=(RM * math.cos(theta), YC + RM * math.sin(theta), z), label=k * NC + j + 1) for k in range(len(zValues) - 1): for j in range(NC): nextJ = (j + 1) % NC indices = [k * NC + j, k * NC + nextJ, (k + 1) * NC + nextJ, (k + 1) * NC + j] add_element(part, indices, QUAD4, k * NC + j + 1) assign_section(part, 'PipeSection', S4R) part.Surface(name='OUTER', side1Elements=part.elements[:]) part.Surface(name='INNER', side2Elements=part.elements[:]) make_node_set(part, 'LEFT', list(range(1, NC + 1))) make_node_set(part, 'RIGHT', list(range((len(zValues) - 1) * NC + 1, len(zValues) * NC + 1))) pipe = assembly.Instance(name='Pipe-1', part=part, dependent=ON) assembly.regenerate() previous = 'Initial' for stepName in ['Gravity', 'InternalPressure', 'FaultSlip']: model.StaticStep(name=stepName, previous=previous, nlgeom=ON, timePeriod=1.0, initialInc=0.01, minInc=1e-07, maxInc=0.05, maxNumInc=2000, stabilizationMethod=DISSIPATED_ENERGY_FRACTION, stabilizationMagnitude=0.0002, adaptiveDampingRatio=0.05) previous = stepName for propertyName in ['PipeSoil', 'FaultFrictionless']: model.ContactProperty(propertyName) model.interactionProperties[propertyName].NormalBehavior(pressureOverclosure=HARD, allowSeparation=ON) model.interactionProperties['PipeSoil'].TangentialBehavior(formulation=PENALTY, directionality=ISOTROPIC, table=((0.6,),), maximumElasticSlip=ABSOLUTE_DISTANCE, absoluteDistance=0.5) for instance in soilInstances: model.SurfaceToSurfaceContactStd(name='Contact_' + instance.name, createStepName='Initial', main=instance.surfaces['HOLE'], secondary=pipe.surfaces['OUTER'], sliding=FINITE, interactionProperty='PipeSoil', thickness=ON, adjustMethod=NONE) model.SurfaceToSurfaceContactStd(name='FaultContact', createStepName='Initial', main=soilInstances[0].surfaces['FAULT'], secondary=soilInstances[1].surfaces['FAULT'], sliding=FINITE, interactionProperty='FaultFrictionless', adjustMethod=NONE) for (i, instance) in enumerate(soilInstances): prefix = 'Soil%d_' % i for (suffix, values) in [('BOTTOM', {'u1': 0.0, 'u2': 0.0, 'u3': 0.0}), ('SIDE', {'u1': 0.0}), ('END', {'u3': 0.0})]: bc = model.DisplacementBC(name=prefix + suffix, createStepName='Initial', region=instance.sets[suffix], **values) if i == 1: changes = {} for component in values.keys(): changes[component] = -SLIP if component == 'u3' else SLIP bc.setValuesInStep(stepName='FaultSlip', **changes) model.DisplacementBC(name='PipeLeft', createStepName='Initial', region=pipe.sets['LEFT'], u1=0.0, u2=0.0) model.DisplacementBC(name='PipeRight', createStepName='Initial', region=pipe.sets['RIGHT'], u1=0.0, u2=0.0, u3=0.0) model.boundaryConditions['PipeRight'].setValuesInStep(stepName='FaultSlip', u1=SLIP, u2=SLIP, u3=-SLIP) area = math.pi * (RO * RO - (RO - THICKNESS) ** 2) assembly.engineeringFeatures.SpringDashpotToGround(name='AxialSpring', region=pipe.sets['LEFT'], dof=3, springBehavior=ON, springStiffness=E_STEEL * area / (35.0 * D * NC), dashpotBehavior=OFF) model.Gravity(name='SelfWeight', createStepName='Gravity', comp2=-9810.0) model.Pressure(name='Pressure', createStepName='InternalPressure', region=pipe.surfaces['INNER'], magnitude=PRESSURE) force = PRESSURE * math.pi * (RO - THICKNESS) ** 2 / NC for (name, sign) in [('LEFT', -1.0), ('RIGHT', 1.0)]: model.ConcentratedForce(name='EndThrust_' + name, createStepName='InternalPressure', region=pipe.sets[name], cf3=sign * force) model.fieldOutputRequests['F-Output-1'].setValues(variables=('S', 'LE', 'PEEQ', 'U', 'RF', 'CSTRESS', 'CDISP'), numIntervals=20) model.HistoryOutputRequest(name='Energy', createStepName='Gravity', variables=('ALLIE', 'ALLSE', 'ALLPD', 'ALLSD', 'ALLWK'), numIntervals=40) assembly.regenerate() job = mdb.Job(name=jobName, model=modelName, numCpus=4, numDomains=4, memory=80, memoryUnits=PERCENTAGE, type=ANALYSIS) _sw_record_write(job, consistencyChecking=ON) inpPath = os.path.abspath(jobName + '.inp') with open(inpPath, 'r') as inputFile: inputText = inputFile.read().upper().replace(' ', '') if '*ELEMENT,TYPE=' not in inputText: raise RuntimeError('Input file contains no typed elements.') mdb.saveAs(pathName=os.path.abspath(modelName + '.cae')) elementCount = 0 for currentPart in model.parts.values(): elementCount += len(currentPart.elements) print('MODEL: ' + modelName) print('INP: ' + inpPath) print('ELEMENTS: %d' % elementCount) print('JOB CREATED; INPUT WRITTEN; NOT SUBMITTED.') 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')