# -*- coding: utf-8 -*- """ Created on Tue Aug 1 14:20:24 2023 @author: rsp1n17 """ # Creation of push-out assembly from abaqus import * from abaqusConstants import * import section import regionToolset import displayGroupMdbToolset as dgm import part import material import assembly import step import interaction import load import mesh import optimization import job import sketch import visualization import xyPlot import displayGroupOdbToolset as dgo import connectorBehavior import numpy import csv def PlaneByThreePoints(part,x1, y1, z1, x2, y2, z2, x3, y3, z3): p=pushoutModel.parts[part] myPlane=p.DatumPlaneByThreePoints(point1=(x1,y1,z1), point2=(x2,y2,z2), point3=(x3,y3,z3)) myID=myPlane.id return myID def CreateSetEdge(x,y,z,part,set_name): edge = () p = pushoutModel.parts[part] e = p.edges myEdge = e.findAt((x,y,z),) edge = edge + (e[myEdge.index:myEdge.index+1], ) p.Set(edges=edge, name=set_name) return myEdge def CreateStudShankExtrude(part, id_plane, edge, x1,y1,z1, d1, h1, h2): p=pushoutModel.parts[part] r=0.5*d1 e, d = p.edges, p.datums t = p.MakeSketchTransform(sketchPlane=d[id_plane], sketchUpEdge=edge, sketchPlaneSide=SIDE1, sketchOrientation=RIGHT, origin=(x1, y1, z1)) s = pushoutModel.ConstrainedSketch(name='StudShank', sheetSize=600, transform=t) g, v, d1, c = s.geometry, s.vertices, s.dimensions, s.constraints s.setPrimaryObject(option=SUPERIMPOSE) p.projectReferencesOntoSketch(sketch=s, filter=COPLANAR_EDGES) studLine1=s.Line(point1=(0,(300+r)), point2=(0,(300-r))) studArc1=s.Arc3Points(point1=(0,(300+r)), point2=(0,(300-r)), point3=(r,300)) studLine2=s.Line(point1=(0,(550+r)), point2=(0,(550-r))) studArc2=s.Arc3Points(point1=(0,(550+r)), point2=(0,(550-r)), point3=(r,550)) p.SolidExtrude(sketchPlane=d[id_plane], sketchUpEdge=edge, sketchPlaneSide=SIDE1, sketchOrientation=RIGHT, sketch=s, flipExtrudeDirection=OFF, depth=(h1-h2)) s.unsetPrimaryObject() def CreateStudHeadExtrude(part, id_plane, edge, x1,y1,z1,d2,h2): p=pushoutModel.parts[part] r=0.5*d2 e, d = p.edges, p.datums t = p.MakeSketchTransform(sketchPlane=d[id_plane], sketchUpEdge=edge, sketchPlaneSide=SIDE1, sketchOrientation=RIGHT, origin=(x1, y1, z1)) s = pushoutModel.ConstrainedSketch(name='StudHead', sheetSize=600, transform=t) g, v, d1, c = s.geometry, s.vertices, s.dimensions, s.constraints s.setPrimaryObject(option=SUPERIMPOSE) p.projectReferencesOntoSketch(sketch=s, filter=COPLANAR_EDGES) studLine1=s.Line(point1=(0,(300+r)), point2=(0,(300-r))) studArc1=s.Arc3Points(point1=(0,(300+r)), point2=(0,(300-r)), point3=(r,300)) studLine2=s.Line(point1=(0,(550+r)), point2=(0,(550-r))) studArc2=s.Arc3Points(point1=(0,(550+r)), point2=(0,(550-r)), point3=(r,550)) p.SolidExtrude(sketchPlane=d[id_plane], sketchUpEdge=edge, sketchPlaneSide=SIDE1, sketchOrientation=RIGHT, sketch=s, flipExtrudeDirection=OFF, depth=h2) s.unsetPrimaryObject() def CreateStudShankCutExtrude(part, id_plane, edge, x1,y1,z1, d1, h1, h2): p=pushoutModel.parts[part] r=0.5*d1 e, d = p.edges, p.datums t = p.MakeSketchTransform(sketchPlane=d[id_plane], sketchUpEdge=edge, sketchPlaneSide=SIDE2, sketchOrientation=RIGHT, origin=(x1, y1, z1)) s = pushoutModel.ConstrainedSketch(name='StudShank', sheetSize=600, transform=t) g, v, d1, c = s.geometry, s.vertices, s.dimensions, s.constraints s.setPrimaryObject(option=SUPERIMPOSE) p.projectReferencesOntoSketch(sketch=s, filter=COPLANAR_EDGES) studLine1=s.Line(point1=(0,(250+r)), point2=(0,(250-r))) studArc1=s.Arc3Points(point1=(0,(250+r)), point2=(0,(250-r)), point3=(r,250)) studLine2=s.Line(point1=(0,(500+r)), point2=(0,(500-r))) studArc2=s.Arc3Points(point1=(0,(500+r)), point2=(0,(500-r)), point3=(r,500)) p.CutExtrude(sketchPlane=d[id_plane], sketchUpEdge=edge, sketchPlaneSide=SIDE2, sketchOrientation=RIGHT, sketch=s, flipExtrudeDirection=OFF, depth=(h1-h2)) s.unsetPrimaryObject() def CreateStudHeadCutExtrude(part, id_plane, edge, x1,y1,z1,d2,h2): p=pushoutModel.parts[part] r=0.5*d2 e, d = p.edges, p.datums t = p.MakeSketchTransform(sketchPlane=d[id_plane], sketchUpEdge=edge, sketchPlaneSide=SIDE2, sketchOrientation=RIGHT, origin=(x1, y1, z1)) s = pushoutModel.ConstrainedSketch(name='StudHead', sheetSize=600, transform=t) g, v, d1, c = s.geometry, s.vertices, s.dimensions, s.constraints s.setPrimaryObject(option=SUPERIMPOSE) p.projectReferencesOntoSketch(sketch=s, filter=COPLANAR_EDGES) studLine1=s.Line(point1=(0,(250+r)), point2=(0,(250-r))) studArc1=s.Arc3Points(point1=(0,(250+r)), point2=(0,(250-r)), point3=(r,250)) studLine2=s.Line(point1=(0,(500+r)), point2=(0,(500-r))) studArc2=s.Arc3Points(point1=(0,(500+r)), point2=(0,(500-r)), point3=(r,500)) p.CutExtrude(sketchPlane=d[id_plane], sketchUpEdge=edge, sketchPlaneSide=SIDE2, sketchOrientation=RIGHT, sketch=s, flipExtrudeDirection=OFF, depth=h2) s.unsetPrimaryObject() #Define parameters d1=[16,19,22,25] stud_height=[60,75,100,150,200,250] C=100 session.viewports['Viewport: 1'].setValues(displayedObject=None) #------------------------------------------------------------ #Create the model for j in range(len(stud_height)): h1=stud_height[j] if h1==60 or h1==75 or h1==100: hc=200 else: hc=h1+50 for i in range(len(d1)): d=d1[i] mdb.Model(name='PA_C{}_d{}_h{}'.format(C,d,h1)) pushoutModel=mdb.models['PA_C{}_d{}_h{}'.format(C,d,h1)] if d==16: d2=32 h2=8 if d==19: d2=32 h2=10 if d==22: d2=35 h2=10 if d==25: d2=41 h2=12 #------------------------------------------------------------ #Create the I-Beam and studs part #sketch one quarter of the I-beam profile beamProfileSketch=pushoutModel.ConstrainedSketch(name='I-Beam Profile', sheetSize=130) Line1=beamProfileSketch.Line(point1=(0,0), point2=(0,130)) Line2=beamProfileSketch.Line(point1=(0, 130), point2=(130, 130)) Line3=beamProfileSketch.Line(point1=(130, 130), point2=(130, 115)) Line4=beamProfileSketch.Line(point1=(130, 115), point2=(5, 115)) Line5=beamProfileSketch.Line(point1=(5, 115), point2=(5, 0)) Line6=beamProfileSketch.Line(point1=(5, 0), point2=(0, 0)) #Create the beam by extruding the sketch beamPart=pushoutModel.Part(name='I-Beam', dimensionality=THREE_D, type=DEFORMABLE_BODY) beamPart.BaseSolidExtrude(sketch=beamProfileSketch, depth=700) #Creat the stud shanks by extruding on the flange studSketchRightEdge=CreateSetEdge(130,130,350,'I-Beam', 'flangeRightEdge') studShankSketchPlane=PlaneByThreePoints('I-Beam', 0, 130, 700, 130,130,700, 0, 130, 0) CreateStudShankExtrude('I-Beam', studShankSketchPlane, studSketchRightEdge, 0,130,700, d, h1,h2) #Create the stud heads studHeadSketchPlane=PlaneByThreePoints('I-Beam', 0,130+h1-h2,700,9.5,130+h1-h2,700,0,130+h1-h2,0) CreateStudHeadExtrude('I-Beam', studHeadSketchPlane, studSketchRightEdge, 0,130+h1-h2,700,d2,h2) #Partition the studs from the beam beamPart.PartitionCellByDatumPlane(cells=beamPart.cells, datumPlane=beamPart.datums[studShankSketchPlane]) #------------------------------------------------------------ #Create the slab part #Sketch the slab profile slabProfileSketch=pushoutModel.ConstrainedSketch(name='Slab profile', sheetSize=700) Line8=slabProfileSketch.Line(point1=(0,30), point2=(0,650)) Line9=slabProfileSketch.Line(point1=(0,650), point2=(300,650)) Line10=slabProfileSketch.Line(point1=(300,650), point2=(300,0)) Line11=slabProfileSketch.Line(point1=(300,0), point2=(100,0)) Line12=slabProfileSketch.Line(point1=(100,0), point2=(100,30)) Line13=slabProfileSketch.Line(point1=(100,30), point2=(0,30)) #Create the slab by extruding the sketch slabPart=pushoutModel.Part(name='Slab', dimensionality=THREE_D, type=DEFORMABLE_BODY) slabPart.BaseSolidExtrude(sketch=slabProfileSketch, depth=hc) #Cut extrude the stud shanks slabExtrudeRightEdge=CreateSetEdge(300,325,hc,'Slab', 'slabRightEdge') slabShankCutPlane=PlaneByThreePoints('Slab', 0,0,hc,0,650,hc,300,0,hc) CreateStudShankCutExtrude('Slab', slabShankCutPlane, slabExtrudeRightEdge, 0,0,hc,d,h1,h2) #Cut extrude the stud heads slabHeadCutPlane=PlaneByThreePoints('Slab', 0,30,hc-h1+h2, 0,650,hc-h1+h2,300,0,hc-h1+h2) CreateStudHeadCutExtrude('Slab', slabHeadCutPlane, slabExtrudeRightEdge, 0,0,hc-h1+h2,d2,h2) #Partition the slab slabPart.PartitionCellByPlaneThreePoints(cells=slabPart.cells, point1=(0,0,hc-30), point2=(100,100,hc-30), point3=(0,100,hc-30)) if h1>60: slabPart.PartitionCellByPlaneThreePoints(cells=slabPart.cells, point1=(0,0,hc-60), point2=(100,100,hc-60), point3=(0,100,hc-60)) if h1>100: slabPart.PartitionCellByPlaneThreePoints(cells=slabPart.cells, point1=(0,0,hc-90), point2=(100,100,hc-90), point3=(0,100,hc-90)) if h1>200: slabPart.PartitionCellByPlaneThreePoints(cells=slabPart.cells, point1=(0,0,hc-120), point2=(100,100,hc-120), point3=(0,100,hc-120)) #------------------------------------------------------------ #Create the rebar part #Sketch the wire profile rebarWireSketch=pushoutModel.ConstrainedSketch(name='Rebar wire profile', sheetSize=700) Line14=rebarWireSketch.Line(point1=(0,65), point2=(285,65)) Line15=rebarWireSketch.Line(point1=(0,165), point2=(285,165)) Line16=rebarWireSketch.Line(point1=(0,315), point2=(285,315)) Line17=rebarWireSketch.Line(point1=(0,465), point2=(285,465)) Line18=rebarWireSketch.Line(point1=(0,615), point2=(285,615)) Line19=rebarWireSketch.Line(point1=(90,45), point2=(90,635)) Line20=rebarWireSketch.Line(point1=(270,45), point2=(270,635)) #Create the rebar from the wire sketch rebarPart=pushoutModel.Part(name='Rebar', dimensionality=THREE_D, type=DEFORMABLE_BODY) rebarPart.BaseWire(sketch=rebarWireSketch) #------------------------------------------------------------ #Create the rigid surface RigidSurfaceSketch=pushoutModel.ConstrainedSketch(name='Rigid surface profile', sheetSize=hc) RigidSurfaceLine=RigidSurfaceSketch.Line(point1=(0,0), point2=(200,0)) RigidSurfacePart=pushoutModel.Part(name='Rigid Surface', dimensionality=THREE_D, type=ANALYTIC_RIGID_SURFACE) RigidSurfacePart.AnalyticRigidSurfExtrude(sketch=RigidSurfaceSketch, depth=hc) RigidSurfacePart.ReferencePoint(point=(0.5*hc,0,0)) #------------------------------------------------------------ #Create the materials #Beam material beamMaterial=pushoutModel.Material(name='Beam material') #### read the data f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/BeamElastic.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]))) elasticity=tuple(Data) f.close() f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/BeamPlastic.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]))) plasticity=tuple(Data) f.close() ####define material beamMaterial.Elastic(table=elasticity) beamMaterial.Plastic(table=plasticity) #Stud material studMaterial=pushoutModel.Material(name='Stud material') #### read the data f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/AusteniticStudElastic.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]))) elasticity=tuple(Data) f.close() f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/AusteniticStudPlastic.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]))) plasticity=tuple(Data) f.close() f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/AusteniticStudInitiation.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]),float(data[2]))) initiation=tuple(Data) f.close() f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/AusteniticStudEvolution.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]),)) evolution=tuple(Data) f.close() studMaterial.Elastic(table=elasticity) studMaterial.Plastic(table=plasticity) studMaterial.DuctileDamageInitiation(table=initiation) studMaterial.ductileDamageInitiation.DamageEvolution(type=DISPLACEMENT, table=evolution, softening=TABULAR) #Concrete material slabMaterial=pushoutModel.Material(name='Concrete') #### read the data f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/C100Elastic.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]))) elasticity=tuple(Data) f.close() f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/CDP.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]), float(data[2]), float(data[3]), float(data[4]))) plasticity=tuple(Data) f.close() f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/C100CH.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]))) CompressionHardening=tuple(Data) f.close() f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/C100TS.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]))) TensionStiffening=tuple(Data) f.close() f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/C100CDamage.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]))) CDamage=tuple(Data) f.close() f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/C100TDamage.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]))) TDamage=tuple(Data) f.close() slabMaterial.Elastic(table=elasticity) slabMaterial.ConcreteDamagedPlasticity(table=plasticity) slabMaterial.concreteDamagedPlasticity.ConcreteCompressionHardening(table=CompressionHardening) slabMaterial.concreteDamagedPlasticity.ConcreteTensionStiffening(type=DISPLACEMENT, table=TensionStiffening) slabMaterial.concreteDamagedPlasticity.ConcreteCompressionDamage(table=CDamage) slabMaterial.concreteDamagedPlasticity.ConcreteTensionDamage(type=DISPLACEMENT, table=TDamage) #Rebar material rebarMaterial=pushoutModel.Material(name='Rebar material') #### read the data f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/RebarElastic.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]))) elasticity=tuple(Data) f.close() f=open('/mainfs/home/rsp1n17/ParametricMarch24/Material Data/RebarPlastic.txt', 'r') Data=[] for line in f: data=line.split() Data.append((float(data[0]),float(data[1]))) plasticity=tuple(Data) f.close() rebarMaterial.Elastic(table=elasticity) rebarMaterial.Plastic(table=plasticity) #------------------------------------------------------------ #Create sections beamSection=pushoutModel.HomogeneousSolidSection(name='I-Beam section', material='Beam material') studSection=pushoutModel.HomogeneousSolidSection(name='Stud section', material='Stud material') slabSection=pushoutModel.HomogeneousSolidSection(name='Slab section', material='Concrete') rebarSection=pushoutModel.TrussSection(name='Rebar section', material='Rebar material', area=78.5) #Assign sections to parts slabRegion=(slabPart.cells,) slabPart.SectionAssignment(region=slabRegion, sectionName='Slab section') rebarRegion=(rebarPart.edges,) rebarPart.SectionAssignment(region=rebarRegion, sectionName='Rebar section') studRegion=(beamPart.cells.findAt(((0,130+h1-h2,150),),((0,130+h1-h2,400),),),) beamPart.SectionAssignment(region=studRegion, sectionName='Stud section') beamRegion=(beamPart.cells.findAt(((0,0,0),),),) beamPart.SectionAssignment(region=beamRegion, sectionName='I-Beam section') #------------------------------------------------------------ #Create the assembly assembly=pushoutModel.rootAssembly beamInstance=assembly.Instance(name='Beam Instance', part=beamPart, dependent=OFF) slabInstance=assembly.Instance(name='Slab Instance', part=slabPart, dependent=OFF) rebarInstance1=assembly.Instance(name='Rebar Instance 1', part=rebarPart, dependent=OFF) rebarInstance2=assembly.Instance(name='Rebar Instance 2', part=rebarPart, dependent=OFF) rigidSurfaceInstance=assembly.Instance(name='Rigid Surface Instance', part=RigidSurfacePart, dependent=OFF) #Position rebar rebarInstance1.translate(vector=(0,0,20)) rebarInstance2.translate(vector=(0,0,180)) #Position beam beamFace=beamInstance.faces.findAt((2.5,5,700),) slabTopFace=slabInstance.faces.findAt((100,650,100),) beamFlangeFace=beamInstance.faces.findAt((10,130,hc),) slabFlangeFace=slabInstance.faces.findAt((10,300,hc),) assembly.FaceToFace(movablePlane=beamFace, fixedPlane=slabTopFace, flip=OFF, clearance=150) assembly.FaceToFace(movablePlane=beamFlangeFace, fixedPlane=slabFlangeFace, flip=ON, clearance=0) #Position Rigid surface rigidSurfaceInstance.translate(vector=(100,0,0.5*hc)) #------------------------------------------------------------- #Partition the assembly #Create the faces to be used for partitions def PartitionByFace(x,y,z, instance): face=assembly.instances[instance].faces.findAt((x,y,z),) beamCells=beamInstance.cells slabCells=slabInstance.cells allCells=beamCells+slabCells assembly.PartitionCellByExtendFace(cells=allCells, extendFace=face) #SlabTopFace PartitionByFace(75,650,0.5*hc, 'Slab Instance') #SlabRecessBottomFace PartitionByFace(50,30,0.5*hc, 'Slab Instance') #SlabRecessSideFace PartitionByFace(100,15,0.5*hc, 'Slab Instance') #FlangeBottomEdgeFace PartitionByFace(65,100,hc+7.5, 'Beam Instance') #WebOutsideFace PartitionByFace(5,350,hc+20, 'Beam Instance') #FlangeOutsideEdgeFace PartitionByFace(130,350,hc+7.5, 'Beam Instance') #FlangeBottomFace PartitionByFace(50,300,hc+15, 'Beam Instance') #StudHeadTop PartitionByFace(5,250,hc-h1, 'Beam Instance') #StudHeadBottom PartitionByFace(1+0.5*d,250,hc-h1+h2, 'Beam Instance') #StudHeadOuter PartitionByFace(0.5*d2,250,hc-h1+0.5*h2, 'Beam Instance') PartitionByFace(0.5*d2,500,hc-h1+0.5*h2, 'Beam Instance') #StudShankOuter PartitionByFace(0.5*d,250,hc-0.5*h1, 'Beam Instance') PartitionByFace(0.5*d,500,hc-0.5*h1, 'Beam Instance') #Partition by datum planes def PartitionByPlaneThreePoints(x1,y1,z1,x2,y2,z2,x3,y3,z3): beamCells=beamInstance.cells slabCells=slabInstance.cells allCells=beamCells+slabCells assembly.PartitionCellByPlaneThreePoints(cells=allCells, point1=(x1,y1,z1), point2=(x2,y2,z2), point3=(x3,y3,z3)) PartitionByPlaneThreePoints(0,250,100,10,250,100,5,250,150) PartitionByPlaneThreePoints(0,500,100,10,500,100,5,500,150) PartitionByPlaneThreePoints(0,375,100,10,375,100,5,375,150) #PartitionByPlaneThreePoints(0,0,hc-0.5*h1,hc,0,hc-0.5*h1,hc,600,hc-0.5*h1) # #Slab15mmpartitions # PartitionByPlaneThreePoints(0,0,hc-30,100,100,hc-30,100,600,hc-30) # if h1>60: # PartitionByPlaneThreePoints(0,0,hc-60,100,100,hc-60,100,600,hc-60) # if h1>100: # PartitionByPlaneThreePoints(0,0,hc-90,100,100,hc-90,100,600,hc-90) # if h1>200: # PartitionByPlaneThreePoints(0,0,hc-120,100,100,hc-120,100,600,hc-120) #---------------------------------------------------------------- #Create surfaces for BCs and contact #Beam and slab contact surface #Slab contact surfaces SlabFlangeFaces=slabInstance.faces.getByBoundingBox(xMin=0, yMin=100, zMin=hc, xMax=130, yMax=650, zMax=hc+10) SlabFlangeSurface=assembly.Surface(side2Faces=SlabFlangeFaces, name='Slab flange surface') SlabStudHeadFacesBottom=slabInstance.faces.getByBoundingBox(xMin=0, yMin=250-0.5*d2, zMin=hc-h1, xMax=0.5*d2, yMax=250+0.5*d2, zMax=hc-h1+h2) SlabStudHeadSurfaceBottom=assembly.Surface(side2Faces=SlabStudHeadFacesBottom, name='Slab stud head bottom surface') SlabStudHeadFacesTop=slabInstance.faces.getByBoundingBox(xMin=0, yMin=500-0.5*d2, zMin=hc-h1, xMax=0.5*d2, yMax=500+0.5*d2, zMax=hc-h1+h2) SlabStudHeadSurfaceTop=assembly.Surface(side2Faces=SlabStudHeadFacesTop, name='Slab stud head top surface') SlabStudShankFacesBottom=slabInstance.faces.getByBoundingBox(xMin=0, yMin=250-0.5*d, zMin=hc-h1+h2, xMax=0.5*d, yMax=250+0.5*d, zMax=hc) SlabStudShankSurfaceBottom=assembly.Surface(side2Faces=SlabStudShankFacesBottom, name='Slab stud shank bottom surface') SlabStudShankFacesTop=slabInstance.faces.getByBoundingBox(xMin=0, yMin=500-0.5*d, zMin=hc-h1+h2, xMax=0.5*d, yMax=500+0.5*d, zMax=hc) SlabStudShankSurfaceTop=assembly.Surface(side2Faces=SlabStudShankFacesTop, name='Slab stud shank top surface') #Beam contact surfaces BeamContactSurfaceFaces=beamInstance.faces.getByBoundingBox(xMin=0, yMin=100, zMin=0, xMax=130, yMax=650, zMax=hc) BeamContactSurface1=assembly.Surface(side2Faces=BeamContactSurfaceFaces, name='Beam surface') #Remove symmetry surface BeamSymmetryFaces=beamInstance.faces.getByBoundingBox(xMin=-10, yMin=0, zMin=0, xMax=0, yMax=800, zMax=600) BeamSymmetrySurface=assembly.Surface(side2Faces=BeamSymmetryFaces, name='Beam symmetry surface') BeamContactSurface2=assembly.SurfaceByBoolean(name='BeamContactSurface', surfaces=(BeamContactSurface1, BeamSymmetrySurface), operation=DIFFERENCE) SlabContactSurface=assembly.SurfaceByBoolean(name='SlabContactSurface', surfaces=(SlabFlangeSurface, SlabStudHeadSurfaceBottom, SlabStudHeadSurfaceTop, SlabStudShankSurfaceBottom, SlabStudShankSurfaceTop)) #YZ Plane symmetry surfaces SlabSymmetryFaces=slabInstance.faces.getByBoundingBox(xMin=-10, yMin=0, zMin=0, xMax=0, yMax=800, zMax=hc) SlabSymmetrySurface=assembly.Surface(side2Faces=SlabSymmetryFaces, name='Slab symmetry surface') #XY Plane symmetry surfaces WebSymmetryFaces=beamInstance.faces.getByBoundingBox(xMin=0, yMin=0, zMin=hc+130, xMax=15, yMax=800, zMax=hc+140) WebSymmetrySurface=assembly.Surface(side2Faces=WebSymmetryFaces, name='Web symmetry surface') #Top of flange for displacement RP FlangeTop=beamInstance.faces.getByBoundingBox(xMin=0, yMin=800, zMin=hc, xMax=130, yMax=810, zMax=hc+130) FlangeTopSurface=assembly.Surface(side2Faces=FlangeTop, name='Flange top surface') #Bottom of slab for BCs SlabBottom=slabInstance.faces.getByBoundingBox(xMin=100, yMin=0, zMin=0, xMax=300, yMax=10, zMax=hc) SlabBottomSurface=assembly.Surface(side2Faces=SlabBottom, name='Slab bottom surface') #--------------------------------------------------------- #Create step pushoutModel.StaticStep(name='DisplacementApplied', previous='Initial', nlgeom=ON, initialInc=0.003, minInc=1E-15, maxNumInc=1000) #---------------------------------------------------- #Apply BCs and contact #Base fixed BC r1 = assembly.instances['Rigid Surface Instance'].referencePoints refPoints1=(r1[2], ) BaseFixedRegion = regionToolset.Region(referencePoints=refPoints1) pushoutModel.DisplacementBC(name='BaseFixed', createStepName='Initial', region=BaseFixedRegion, u1=0, u2=0, u3=0, ur1=0, ur2=0, ur3=0) #Applied displacement------------------------------------------------------------------------------------------------------------------------ DisplacementRP=assembly.ReferencePoint(instanceName='DisplacementRP', point=(0,800,330)) DisplacementRPRegion=regionToolset.Region(referencePoints=(assembly.referencePoints[DisplacementRP.id],)) TieRegion=regionToolset.Region(faces=FlangeTop) pushoutModel.RigidBody(name='DisplacementRP', refPointRegion=DisplacementRPRegion, tieRegion=TieRegion) pushoutModel.DisplacementBC(name='Displacement', createStepName='DisplacementApplied', region=TieRegion, u1=0, u2=-12,u3=0, ur1=0, ur2=0, ur3=0) #Symmetry on the web XYSymmetryRegion=regionToolset.Region(faces=(WebSymmetryFaces)) pushoutModel.ZsymmBC(name='WebSymmetry', createStepName='Initial', region=XYSymmetryRegion) #Symmetry on the beam and slab YZSymmetryRegion=regionToolset.Region(faces=(BeamSymmetryFaces + SlabSymmetryFaces)) pushoutModel.XsymmBC(name='CrossSectionSymmetry', createStepName='Initial', region=YZSymmetryRegion) #---------------------------------------------------------------------------------------------------------------------------------------------- #Steel to concrete contact #Create interaction property HardContactFriction=pushoutModel.ContactProperty(name='Hard Contact') HardContactFriction.NormalBehavior(pressureOverclosure=HARD, allowSeparation=ON, contactStiffness=DEFAULT, contactStiffnessScaleFactor=1.0, clearanceAtZeroContactPressure=0.0, constraintEnforcementMethod=AUGMENTED_LAGRANGE) HardContactFriction.TangentialBehavior(formulation=PENALTY, table=((0.4,),), fraction=0.005) #Create contact interaction pushoutModel.SurfaceToSurfaceContactStd(name='SteelConcreteContact', createStepName='Initial', main=SlabContactSurface, secondary=BeamContactSurface2, sliding=SMALL, interactionProperty='Hard Contact', adjustMethod=OVERCLOSED) #Base rigid surface to concrete contact #Create interaction property BaseFriction=pushoutModel.ContactProperty(name='Base friction') BaseFriction.NormalBehavior(pressureOverclosure=HARD, allowSeparation=ON, contactStiffness=DEFAULT, contactStiffnessScaleFactor=1.0, clearanceAtZeroContactPressure=0.0, constraintEnforcementMethod=AUGMENTED_LAGRANGE) BaseFriction.TangentialBehavior(formulation=PENALTY, table=((0.45,),), fraction=0.005) rigidSurfaceTop= rigidSurfaceInstance.faces.getByBoundingBox(xMin=0, yMin=0, zMin=0, xMax=300, yMax=10, zMax=hc) RigidSurfaceTop=assembly.Surface(side1Faces=rigidSurfaceTop, name='Rigid Surface Top') pushoutModel.SurfaceToSurfaceContactStd(name='SlabBaseContact', createStepName='Initial', main=RigidSurfaceTop, secondary=SlabBottomSurface, sliding=FINITE, interactionProperty='Base friction') #Rebar embedded in concrete Rebar1=rebarInstance1.edges[:] Rebar2=rebarInstance2.edges[:] Rebar=regionToolset.Region(edges=Rebar1+Rebar2) ConcreteHost=regionToolset.Region(cells=slabInstance.cells[:]) pushoutModel.EmbeddedRegion(name='RebarEmbedded', embeddedRegion=Rebar,hostRegion=ConcreteHost) #------------------------------------------------------ #Create mesh #Slab and beam mesh def MeshParts(InstanceName_str,Size): partInstances =(assembly.instances[InstanceName_str], ) assembly.seedPartInstance(regions=partInstances, size=Size, deviationFactor=0.1, minSizeFactor=0.1) for i in assembly.instances.keys(): if i=='Beam Instance': MeshParts(i,15) elif i=='Slab Instance': MeshParts(i,15) else: MeshParts(i,30) #Seed stud edges Stud1Edges=beamInstance.edges.getByBoundingBox(xMin=0, yMin=250-0.5*d2, zMin=hc-h1+h2, xMax=30, yMax=250+0.5*d2, zMax=hc) Stud2Edges=beamInstance.edges.getByBoundingBox(xMin=0, yMin=500-0.5*d2, zMin=hc-h1+h2, xMax=30, yMax=500+0.5*d2, zMax=hc) StudEdges=Stud1Edges + Stud2Edges assembly.seedEdgeBySize(edges=StudEdges, size=1.5, constraint=FINER) for iname in assembly.instances.keys(): partInstances =(assembly.instances[iname], ) assembly.generateMesh(regions=partInstances) elemType1 = mesh.ElemType(elemCode=C3D8R, elemLibrary=STANDARD, kinematicSplit=AVERAGE_STRAIN, secondOrderAccuracy=OFF, hourglassControl=DEFAULT, distortionControl=DEFAULT, viscosity=0.02, elemDeletion=ON, maxDegradation=1.0) pickedRegions = beamInstance.cells[:] + slabInstance.cells[:] AllCells=assembly.Set(name='All cells', cells=pickedRegions) assembly.setElementType(regions=AllCells, elemTypes=(elemType1,)) elemType2 = mesh.ElemType(elemCode=T3D2) rebarRegions=rebarInstance1.edges[:] + rebarInstance2.edges[:] rebarSet=assembly.Set(name='Rebar', edges=rebarRegions) assembly.setElementType(regions=rebarSet, elemTypes=(elemType2,)) #----------------------------------------------------- #Field output requests pushoutModel.fieldOutputRequests['F-Output-1'].setValues(variables=('S','PE','PEEQ','PEMAG','LE','U','RF','DAMAGEC','DAMAGET','SDEG', 'CSTRESS', 'CFORCE')) #History output requests pushoutModel.historyOutputRequests['H-Output-1'].setValues(region=DisplacementRPRegion, variables=('RF2', 'RF3')) pushoutModel.historyOutputRequests.changeKey(fromName='H-Output-1', toName='Load') SlipPointsSlab=slabInstance.vertices.findAt(((130,250,hc),)) SlipPointsSlabSet=assembly.Set(name='Slab displacement', vertices=SlipPointsSlab) SlipPointsSlabRegion=assembly.sets['Slab displacement'] SlipPointBeam=beamInstance.vertices.findAt(((130,250,hc+15),)) SlipPointBeamSet=assembly.Set(name='Beam displacement', vertices=SlipPointBeam) SlipPointBeamRegion=assembly.sets['Beam displacement'] pushoutModel.HistoryOutputRequest(name='SlabDisplacement', createStepName='DisplacementApplied', region=SlipPointsSlabRegion, variables=('U2',)) pushoutModel.HistoryOutputRequest(name='BeamDisplacement', createStepName='DisplacementApplied', region=SlipPointBeamRegion, variables=('U2',)) SlabFrictionFaceRegion=regionToolset.Region(faces=SlabFlangeFaces) #pushoutModel.HistoryOutputRequest(name='Slab Friction', createStepName='DisplacementApplied', region=SlabFrictionFaceRegion, variables=('CSTRESS',)) #------------------------------------------------------- #Create job mdb.Job(name='PA_C{}_d{}_h{}'.format(C,d,h1), model='PA_C{}_d{}_h{}'.format(C,d,h1), description='', type=ANALYSIS, atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', scratch='', resultsFormat=ODB, multiprocessingMode=DEFAULT, numCpus=4, numDomains=4, numGPUs=0) #Run the job mdb.jobs['PA_C{}_d{}_h{}'.format(C,d,h1)].submit(consistencyChecking=OFF) #Do not return control till job is finished running mdb.jobs['PA_C{}_d{}_h{}'.format(C,d,h1)].waitForCompletion() #-------------------------------------------------------- #POSTPROCESSING import odbAccess import visualization ODB_path='PA_C{}_d{}_h{}.odb'.format(C,d,h1) odb_object=session.openOdb(name=ODB_path) #ODB get data keyarray=session.odbData['PA_C{}_d{}_h{}.odb'.format(C,d,h1)].historyVariables.keys() BeamDisplacementOutputVariableName=[] SlabDisplacementOutputVariableName=[] LoadOutputVariableName=[] for x in keyarray: if (x.find('U2')>-1) and (x.find('Beam')>-1): BeamDisplacementOutputVariableName.append(x) for x in keyarray: if (x.find('U2'))>-1 and (x.find('Slab')>-1): SlabDisplacementOutputVariableName.append(x) for x in keyarray: if (x.find('RF2')>-1): LoadOutputVariableName.append(x) BeamDisplacementData=session.XYDataFromHistory(name='Beam Displacement', odb=odb_object, outputVariableName=BeamDisplacementOutputVariableName[0], steps=('DisplacementApplied',),) SlabDisplacementData=session.XYDataFromHistory(name='Slab Displacement', odb=odb_object, outputVariableName=SlabDisplacementOutputVariableName[0], steps=('DisplacementApplied',),) RFData=session.XYDataFromHistory(name='Load', odb=odb_object, outputVariableName=LoadOutputVariableName[0], steps=('DisplacementApplied',),) Slip=-1*(BeamDisplacementData-SlabDisplacementData) Load=-0.004*RFData r=open('PA_C{}_d{}_h{}_LoadSlipData.txt'.format(C,d,h1),'w') ###pathData is the location of txt file and ModelName is the abaqus model slip=[] load=[] for i in range(len(Slip)): a=Slip[i] b=a[1] slip.append(b) c=Load[i] e=c[1] load.append(e) r.write(repr(slip[i])+'\t'+repr(load[i])+'\n') r.close() P_u=max(load) P_Rk=0.9*P_u idx_P_u=load.index(P_u) L=len(load) Postpeak=load[idx_P_u:L] idx_2=min(range(len(Postpeak)), key=lambda i: abs(Postpeak[i]-P_Rk)) idx_d_u=idx_P_u+idx_2-1 d_u=slip[idx_d_u] r2=open('PA_C{}_d{}_h{}_PRk_du.txt'.format(C,d,h1),'w') r2.write(repr(P_Rk)+'\t'+repr(d_u)) r2.close() XYData=session.XYDataFromFile(name='LoadSlipData', fileName='PA_C{}_d{}_h{}_LoadSlipData.txt'.format(C,d,h1)) LoadSlipCurve=session.Curve(xyData=XYData) if 'LoadSlipPlot' in session.xyPlots.keys(): del session.xyPlots['LoadSlipPlot'] LoadSlipPlot=session.XYPlot(name='LoadSlipPlot') chartName=LoadSlipPlot.charts.keys()[0] chart=LoadSlipPlot.charts[chartName] chart.setValues(curvesToPlot=(LoadSlipCurve,),) LoadSlipViewport=session.Viewport(name='LoadSlip') LoadSlipViewport.setValues(displayedObject=LoadSlipPlot) StressViewport=session.Viewport(name='Stress') StressViewport.setValues(displayedObject=odb_object, width=200, height=120) StressViewport.odbDisplay.setPrimaryVariable(variableLabel='S', refinement=(INVARIANT, 'Mises'), outputPosition=INTEGRATION_POINT) StressViewport.odbDisplay.display.setValues(plotState=(CONTOURS_ON_DEF,)) StressViewport.view.setValues(session.views['Left']) StressViewport.view.setValues(width=2000,height=1500, projection=PARALLEL) session.printToFile(fileName='PA_C{}_d{}_h{}_Stress'.format(C,d,h1), format=PNG, canvasObjects=(StressViewport,))