''' ----------------------------------------------------------------------------- Semi elliptic crack in a rectangular plate modeled using continuum elements (C3D8). ----------------------------------------------------------------------------- ''' from abaqus import * import testUtils testUtils.setBackwardCompatibility() from abaqusConstants import * import part, material, section, assembly, step, interaction import regionToolset, displayGroupMdbToolset as dgm, mesh, load, job import inpReader #---------------------------------------------------------------------------- # Create a model Mdb() modelName = 'RectBlEllipticCrC3D8' myModel = mdb.Model(name=modelName) # Create a new viewport in which to display the model # and the results of the analysis. myViewport = session.Viewport(name=modelName) myViewport.makeCurrent() myViewport.maximize() #--------------------------------------------------------------------------- # Create a part # Create a sketch for the base feature mySketch = myModel.Sketch(name='barProfile',sheetSize=200.0) mySketch.sketchOptions.setValues(viewStyle=REGULAR) mySketch.setPrimaryObject(option=STANDALONE) mySketch.rectangle(point1=(0.0, 0.0), point2=(13.333, 26.667)) myBlock = myModel.Part(name='block', dimensionality=THREE_D, type=DEFORMABLE_BODY) myBlock.BaseSolidExtrude(sketch=mySketch, depth=1.6667) mySketch.unsetPrimaryObject() del myModel.sketches['barProfile'] # Create a partition for the crack line face1 = myBlock.faces.findAt((5,0,0.6),) edge1 = myBlock.edges.findAt((13.333,0,0.83335),) t = myBlock.MakeSketchTransform(sketchPlane=face1, sketchUpEdge=edge1, sketchPlaneSide=SIDE1, origin=(6.6665, 0.0, 0.83335)) mySketch = myModel.Sketch(name='barProfile', sheetSize=26.87, gridSpacing=0.67, transform=t) g = mySketch.geometry mySketch.setPrimaryObject(option=SUPERIMPOSE) myBlock.projectReferencesOntoSketch(sketch=mySketch, filter=COPLANAR_EDGES) mySketch.sketchOptions.setValues(gridOrigin=(-6.6665, 0.83335)) mySketch.EllipseByCenterPerimeter(center=(-6.6665, 0.83335), axisPoint1=(-2.4998, 0.83335), axisPoint2=(-6.6665, -0.16665)) f = myBlock.faces pickedFaces = f[face1.index:(face1.index+1)] myBlock.PartitionFaceBySketch(sketchUpEdge=edge1, faces=pickedFaces, sketch=mySketch) mySketch.unsetPrimaryObject() del myModel.sketches['barProfile'] # Create a shell tube for the inner partition (around the crack tip) # First create the sweep path mySketch1 = myModel.Sketch(name='sweepPath', sheetSize=50.0) g1 = mySketch1.geometry mySketch1.setPrimaryObject(option=STANDALONE) mySketch1.EllipseByCenterPerimeter(center=(0.0, 0.0), axisPoint1=(4.1667, 0.0), axisPoint2=(0.0, -1.0)) mySketch1.Line(point1=(0.0, 0.0), point2=(0.0, -5.0)) mySketch1.Line(point1=(0.0, 0.0), point2=(6.0, 0.0)) mySketch1.autoTrimCurve(curve=g1[3], parameter1=0.655947506427765) mySketch1.delete(objectList=(g1[5], g1[7])) mySketch1.unsetPrimaryObject() # Create a circle and sweep it along the curve created above mySketch2 = myModel.Sketch(name='innerCircleProfile', sheetSize=50.0, transform=(-1.46951586845991e-17, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.46951586845991e-17, 0.0, 2.55128364723585e-16, -1.0, 0.0)) g2 = mySketch2.geometry mySketch2.setPrimaryObject(option=SUPERIMPOSE) mySketch2.ObliqueConstructionLine(point1=(-23.57, 0.0), point2=(23.57, 0.0)) mySketch2.ObliqueConstructionLine(point1=(0.0, -23.57), point2=(0.0, 23.57)) mySketch2.CircleByCenterPerimeter(center=(0.0, 0.0), point1=(0.005, 0.0)) myInnerTube = myModel.Part(name='innerTube', dimensionality=THREE_D, type=DEFORMABLE_BODY) myInnerTube.BaseShellSweep(sketch=mySketch2, path=mySketch1) mySketch2.unsetPrimaryObject() myViewport.setValues(displayedObject=myInnerTube) del myModel.sketches['innerCircleProfile'] # Create a shell tube for the outer partition (around the inner tube) # For outer partition the sweep path is same as inner partition mySketch2 = myModel.Sketch(name='outerCircleProfile', sheetSize=50.0, transform=(-1.46951586845991e-17, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.46951586845991e-17, 0.0, -7.65385094170755e-16, -1.0, 0.0)) g3 = mySketch2.geometry mySketch2.setPrimaryObject(option=SUPERIMPOSE) mySketch2.ObliqueConstructionLine(point1=(-23.57, 0.0), point2=(23.57, 0.0)) mySketch2.ObliqueConstructionLine(point1=(0.0, -23.57), point2=(0.0, 23.57)) mySketch2.CircleByCenterPerimeter(center=(0.0, 0.0), point1=(0.2, 0.0)) myOuterTube = myModel.Part(name='outerTube', dimensionality=THREE_D, type=DEFORMABLE_BODY) myOuterTube.BaseShellSweep(sketch=mySketch2, path=mySketch1) mySketch2.unsetPrimaryObject() myViewport.setValues(displayedObject=myOuterTube) del myModel.sketches['outerCircleProfile'] del myModel.sketches['sweepPath'] #--------------------------------------------------------------------------- # Create an assembly # Create an instance of the bar myAssembly = myModel.rootAssembly myViewport.setValues(displayedObject=myAssembly) myAssembly.DatumCsysByDefault(CARTESIAN) myAssembly.Instance(name='myBlock-1', part=myBlock) myBlockInstance = myAssembly.instances['myBlock-1'] # Create an instance of the outer partition tube myAssembly.Instance(name='myOuterTube-1', part=myOuterTube, dependent=OFF) myOuterTubeInstance = myAssembly.instances['myOuterTube-1'] # Position the outer tube around the crack tip myOuterTubeInstance.rotateAboutAxis(axisPoint=(0.0, 0.0, 0.0), axisDirection=(13.333, 0.0, 0.0), angle=90.0) myOuterTubeInstance.translate(vector= (8.88178419700125e-16, 0.0, 1.6667)) # Create an instance of the inner partition tube myAssembly.Instance(name='myInnerTube-1', part=myInnerTube, dependent=OFF) myInnerTubeInstance = myAssembly.instances['myInnerTube-1'] # Position the inner tube around the crack tip myInnerTubeInstance.rotateAboutAxis(axisPoint=(0.0, 0.0, 0.0), axisDirection=(13.333, 0.0, 0.0), angle=90.0) myInnerTubeInstance.translate(vector= (8.88178419700125e-16, 0.0, 1.6667)) # Subtract the inner and outer tubes from the block myAssembly.PartFromBooleanCut(name='partitionedBlock', instanceToBeCut=myBlockInstance, cuttingInstances=(myInnerTubeInstance, myOuterTubeInstance, )) # Create an instance of the partitioned block myPtBlock = myModel.parts['partitionedBlock'] myAssembly.Instance(name='myPtBlock-1', part=myPtBlock, dependent=OFF) myPtBlInstance = myAssembly.instances['myPtBlock-1'] # Create a new assembly of the parts to be suppressed myAssembly1 = myModel.rootAssembly myAssembly1.suppressFeatures(('myInnerTube-1', 'myOuterTube-1', 'myBlock-1', )) # Create a set for the new bar instance cells = myPtBlock.cells myPtBlock.Set(cells=cells, name='All') #--------------------------------------------------------------------------- # Assign material properties myViewport.setValues(displayedObject=myPtBlock) # Create linear elastic material myModel.Material(name='LinearElastic') myModel.materials['LinearElastic'].Elastic(table=((30000000.0, 0.3), )) myModel.HomogeneousSolidSection(name='SolidHomogeneous', material='LinearElastic', thickness=1.0) region = myPtBlock.sets['All'] # Assign the above section to the part (myPtBlock.SectionAssignment(region=region, sectionName='SolidHomogeneous')) #--------------------------------------------------------------------------- # Create a step for applying a load myModel.StaticStep(name='ApplyLoad', previous='Initial', description='Apply the load') # Create a partition to make the entire part meshable c1 = myPtBlInstance.cells e1 = myPtBlInstance.edges.findAt((2.217181,0.,0.82003)) pickedEdges = (e1,) sPath = myPtBlInstance.edges.findAt((0.,13.3335,1.6667)) myAssembly.PartitionCellBySweepEdge(sweepPath=sPath, cells=c1, edges=pickedEdges) myViewport.setValues(displayedObject=myAssembly) #--------------------------------------------------------------------------- # Create interaction properties # Create a set for the top surface surf = myPtBlInstance.faces surf1 = myPtBlInstance.faces.findAt((6.6665,26.667,0.6)) surf2 = myPtBlInstance.faces.findAt((2.08335,26.667,0.9)) surfAll = (surf[surf1.index:(surf1.index+1)] + surf[surf2.index:(surf2.index+1)]) myAssembly.Surface(side1Faces=surfAll, name='topSurf') # Create a set for the X face faces1 = (myPtBlInstance.faces.findAt(((6.6665,0,0.8),), ((2.3,0,0.7),), ((2.179725,0.,0.81192),),)) myAssembly.Set(faces=faces1, name='xFace') # Create a set for Y face faces1 = (myPtBlInstance.faces.findAt(((0,13.3335,0.4),), ((0,13.3335,1.0),), ((0,0.1,0.541),), ((0,0.1,0.808),), ((0,0.02,0.641),), ((0,0.02,0.691),),)) myAssembly.Set(faces=faces1, name='yFace') # Create a set for the vertex to be fixed in Z v1 = myPtBlInstance.vertices.findAt(((0.,0.,1.6667),),) myAssembly.Set(vertices=v1, name='pt') # Create a set for the crack line edges = myPtBlInstance.edges edges1 = myPtBlInstance.edges.findAt((2.217181,0.,0.82003)) edges1 = edges[edges1.index:(edges1.index+1)] myAssembly.Set(edges=edges1, name='crackLine') # Create the contour integral definition for the crack crackFront = crackTip = myAssembly.sets['crackLine'] v1 = myPtBlInstance.vertices.findAt((0.,0.,0.6667)) v2 = myPtBlInstance.vertices.findAt((0.,26.667,0.6667)) myAssembly.engineeringFeatures.ContourIntegral(name='Crack', symmetric=ON, crackFront=crackFront, crackTip=crackTip, extensionDirectionMethod=CRACK_NORMAL, crackNormal=(v1, v2), midNodePosition=0.25, collapsedElementAtTip=SINGLE_NODE,crackTipSense=REVERSE) #--------------------------------------------------------------------------- # Create loads and boundary conditions # Fix the X face in the Y directions region = myAssembly.sets['xFace'] myModel.DisplacementBC(name='yFixed', createStepName='Initial', region=region, u2=0.0, fixed=OFF, distributionType=UNIFORM, localCsys=None) # Fix the Y face in the X direction region = myAssembly.sets['yFace'] myModel.DisplacementBC(name='xFixed', createStepName='Initial', region=region, u1=0.0, fixed=OFF, distributionType=UNIFORM, localCsys=None) # Fix a vertex in the Z direction region = myAssembly.sets['pt'] myModel.DisplacementBC(name='zFixed', createStepName='Initial', region=region, u3=0.0, fixed=OFF, distributionType=UNIFORM, localCsys=None) # Apply a pressure load on the top face region = myAssembly.surfaces['topSurf'] myModel.Pressure(name='prLoad', createStepName='ApplyLoad', region=region, distributionType=UNIFORM, magnitude=-1.0, amplitude=UNSET) #--------------------------------------------------------------------------- # Create a mesh # Seed all the edges pickedEdges1 = (myPtBlInstance.edges.findAt(((4.308121,0.141421,1.6667),), ((4.025279,0.141421,1.6667),), ((0.,0.141396,0.808146),), ((0.,0.141419,0.525276),))) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=6, constraint=FIXED) pickedEdges1 = (myPtBlInstance.edges.findAt(((4.1667,0.025,1.6667),), ((4.1417,0,1.6667),), ((4.1917,0,1.6667),), ((0,0.025,0.6667),), ((0.,0.,0.6417),), ((0.,0.,0.6917),))) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=1, constraint=FIXED) pickedEdges1 = (myPtBlInstance.edges.findAt(((0.,0.125,0.6667),), ((0.,0.,0.5417),), ((0.,0.,0.7917),), ((4.1667,0.125,1.6667),), ((4.0417,0.,1.6667),), ((4.2917,0.,1.6667),))) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5, constraint=FIXED) pickedEdges1 = (myPtBlInstance.edges.findAt(((2.217181,0.,0.82003),), ((1.090774,0.,0.696563),), ((2.372976,0.,0.642024),), ((1.089261,0.,0.701475),), ((2.06162,0.,0.999576),), ((1.087582,0.,0.706376),))) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=20, constraint=FIXED) pickedEdges1 = myPtBlInstance.edges.findAt(((1.98335,0.,1.6667),)) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=20, constraint=FIXED) pickedEdges1 = myPtBlInstance.edges.findAt(((0.,0.,0.68455),)) pickedEdges2 = (myPtBlInstance.edges.findAt(((0.,0.014928,0.6667),), ((0.,0.,0.651772),))) myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2, ratio=3.0, number=10) ''' pickedEdges1 = (myPtBlInstance.edges.findAt(((4.186988,0.,1.6667),), ((4.1667,0.030668,1.6667),))) pickedEdges2 = myPtBlInstance.edges.findAt(((4.141774,0.,1.6667),)) myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2, ratio=3.0, number=10) ''' pickedEdges1 = myPtBlInstance.edges.findAt(((4.141774,0.,1.6667),)) myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, ratio=3.0, number=10, constraint=FIXED) pickedEdges1 = myPtBlInstance.edges.findAt(((4.1667,0.024933,1.6667),)) myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, ratio=3.0, number=10, constraint=FIXED) pickedEdges2 = myPtBlInstance.edges.findAt(((4.181774,0.,1.6667),)) myAssembly.seedEdgeByBias(end2Edges=pickedEdges2, ratio=3.0, number=10, constraint=FIXED) pickedEdges1 = myPtBlInstance.edges.findAt(((0.,0.,1.2667),)) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5) pickedEdges1 = myPtBlInstance.edges.findAt(((0.,0.,0.23335),)) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5, constraint=FIXED) pickedEdges1 = (myPtBlInstance.edges.findAt(((13.333,0.,0.83335),), ((13.333,26.667,0.83335),))) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5) pickedEdges1 = (myPtBlInstance.edges.findAt(((8.84985,0.,1.6667),), ((8.74985,26.667,1.6667),))) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=15, constraint=FIXED) pickedEdges1 = (myPtBlInstance.edges.findAt(((6.6665,0.,0),), ((6.6665,26.667,0),))) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=25, constraint=FIXED) pickedEdges1 = myPtBlInstance.edges.findAt(((2.217181,26.667,0.82003),)) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=20, constraint=FIXED) pickedEdges1 = myPtBlInstance.edges.findAt(((2.08335,26.667,1.6667),)) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=20, constraint=FIXED) pickedEdges1 = myPtBlInstance.edges.findAt(((0.,26.667,1.1667),)) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5) pickedEdges1 = myPtBlInstance.edges.findAt(((0.,26.667,0.33335),)) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5, constraint=FIXED) pickedEdges1 = (myPtBlInstance.edges.findAt(((4.1667,13.4335,1.6667),), ((0.,13.4335,0.6667),), ((13.333,13.3335,1.6667),), ((0.,13.3335,1.6667),), ((0.,13.3335,0.),), ((13.333,13.3335,0.),))) myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=30, constraint=FIXED) # Assign meshing controls to the respective cells cells = myPtBlInstance.cells.findAt(((2.167847,0.,0.916656),), ((2.378009,0.,0.742006),),) myAssembly.setMeshControls(regions=cells, elemShape=HEX, technique=STRUCTURED) cells = myPtBlInstance.cells.findAt(((2.271117,0.,0.830836),), ((2.276242,0.,0.826576),),) myAssembly.setMeshControls(regions=cells, elemShape=WEDGE, technique=SWEEP) # All the remaining cells will be meshed with hex elements using # structured meshing. This is also the default elemType1 = mesh.ElemType(elemCode=C3D8, elemLibrary=STANDARD) elemType2 = mesh.ElemType(elemCode=C3D6, elemLibrary=STANDARD) elemType3 = mesh.ElemType(elemCode=C3D4, elemLibrary=STANDARD) cells1 = myPtBlInstance.cells pickedRegions =(cells1, ) myAssembly.setElementType(regions=pickedRegions, elemTypes=(elemType1, elemType2, elemType3)) partInstances = (myPtBlInstance, ) myAssembly.generateMesh(regions=partInstances) #--------------------------------------------------------------------------- # Request history output for the crack myModel.historyOutputRequests.changeKey(fromName='H-Output-1', toName='JInt') myModel.historyOutputRequests['JInt'].setValues(contourIntegral='Crack', numberOfContours=6) myModel.HistoryOutputRequest(name='StrInt', createStepName='ApplyLoad', contourIntegral='Crack', numberOfContours=6, contourType=K_FACTORS) myModel.HistoryOutputRequest(name='TStr', createStepName='ApplyLoad', contourIntegral='Crack', numberOfContours=6, contourType=T_STRESS) #--------------------------------------------------------------------------- # Create the job myJob = mdb.Job(name=modelName, model=modelName, description='Contour integral analysis') mdb.saveAs(pathName=modelName) #---------------------------------------------------------------------------