[Abaqus]能幹什麼之Xfem
你可能會想到
這是個複雜的問題，官方是這樣說的：是一種解決斷裂力學問題的新的有限元方法，其理論最早亍 1999 年，由美國西北大學的教授 Belyschko 和 Black 首次提出，主要是採用獨立亍網格剖分的思想解決有限元中的裂紋擴展問題，在保留傳統有限元所有優點的同時，並不需要對結構內部存在的裂紋等缺陷進行網格劃分。
等等這個重要嗎？
這個當然不重要！！作為工科生我們只關心兩件事
第一，這個玩意怎麼搞，有個啥用
第二，這個玩意怎麼用
但是，作為一名搞科研的人，還是要把原理大致闡述一下，就是ABAQUS在6.7版本引入XFEM這種演算法可以很好地描述，斷裂的特性，應力、裂紋擴展的情況，在機械、岩石有著廣闊的應用，那麼作為一個工科生，我可以負責人的講，前方高能預警
要是用還是要干這兩檔子事，
第一，選擇模型中可能出現裂紋的區域，將其單元設置為具有擴展有限元性質的富集
單元。第二，選擇合適的破壞準則，使得單元在達到條件時發生破壞，裂紋得以擴展。
首先，在材料中添加Maxps Damage ，這個是傳說中的損傷係數，用來說明破壞
準則，Damage Evolution是損傷演化參數 主要用來損傷的發展情況
Damage Stabilization Cohesive 是損傷穩定性係數，改善收斂用的
第二，檔子事就是告訴電腦損傷區在哪？在哪，自己定義
當然， ABAQUS 允許在裂縫發展區域中插入一個初始裂縫，這個初始裂縫可以在 Assembly里進行組裝。
其實單元內部的位移函數（形函數）可以是任意形狀的，但大多數的計算軟體都採用了多項式戒者揑值多項式作為手段來描述單元內部的位移場，這是因為採用這種方法更加便亍在編程中進行處理。但是這種方法的缺點就是，由亍形函數的連續性，導致單元內部不可能存在間斷。直到 Belytschko 提出採用水平集函數作為手段這種擴充形函數能夠描述單元內位移場在裂縫兩邊的跳躍性，
但是由於裂縫存在於單元內部，其擴展獨立不其他單元，使得計算變得高效。但是這種方法也存在一些問題，XFEM 採用的形函數模式會導致其求解方程很容易接近線形相關，極大的增加了收斂難度，因此導致 XFEM 方法一直沒有辦法得到很好的推廣。 其實 ABAQUS 在集成 XFEM 方法時做了大量的簡化，目的都是減小求解的難度， ABAQUS 的幫助在介紹 XFEM 的時候其實說：
The following limitations exist with an enriched feature:
幾個限制：
An enriched element cannot be intersected by more than one crack.
不能模擬交叉裂縫
A crack is not allowed to turn more than 90° in one increment during an analysis.
裂縫擴展方向不能偏轉超過90
Only asymptotic crack-tip fields in an isotropic elastic material are considered for a stationary crack.
對於制定裂縫的方法來講，裂尖漸近位移場計算只能使用各項異性材料
Adaptive remeshing is not supported.
自適應網格也不能用
Composite solid elements are not supported.
複合材料單元也不能用
此外，計算過程中很容易發現裂縫是丌能停留在單元內部的，這說明了ABAQUS 放棄了單元內部對裂尖的描述。同時，由亍 ABAQUS 在計算 XFEM 的損傷時採用的是基亍能量釋放率的 paris 法則，雖然這是基亍彈塑性斷裂力學的經典手段，但由亍承認了裂尖位置的塑性效應，使得在模擬損傷時也只能對低周疲勞能有比較好的近似。
總之，儘管有如此之多的限制，但 XFEM 方法依然為有限元領域注入了一劑新鮮的空氣，它所提供的是一種廣義有限元的實現過程，為有限元理論的發展提供了一種新的思維方式。
這是極官方的評價
考慮大家基礎比較薄弱，下面較大家如何添加裂縫:
第一步，做一個Part
第四步，編輯裂縫
第七步，為了保證能夠看清楚裂縫擴展的過程，要把求解增量步變小
第七步，定義損傷穩定係數
第八步，選取輸出場變數
第九步，畫網格
咱也有法子，你把下邊的代碼丟到最下面的那個命令框里，就好了。
# === Modules ===from part import *from material import *from section import *from assembly import *from step import *from regionToolset import *from interaction import *from load import *from mesh import *from job import *from sketch import *from visualization import *from connectorBehavior import *#---------------------------------------------------------------------------------# === Parameters ===modelName = "modeII_2d_lefm"length = 3.0 # plate half-lengthheight = 6.0 # plate heightwidth = 1.0 # plate widthcrack_len = 1.5 # crack lengthcrack_y = 0.05 # crack offset in y-directionYM = 2.10e+11 # Young"s modulusNU = 0.3 # Poisson"s ratioMAXPS = 8.44e7 # Damage initiationDTOL = 0.05 # Damage toleranceGI = 42200.0 # Fracture energyETA = 1.0 # Power-law exponent#---------------------------------------------------------------------------------# === Create model ===Mdb()viewportName = session.Viewport(name=modelName)viewportName.makeCurrent()viewportName.maximize()plateModel = mdb.Model(name=modelName)del mdb.models["Model-1"]#---------------------------------------------------------------------------------# === Create parts ===plateSketch = plateModel.ConstrainedSketch(name="plateProfile",sheetSize=height)plateSketch.setPrimaryObject(option=STANDALONE)plateSketch.sketchOptions.setValues(decimalPlaces=3)plateSketch.Line(point1=(0.0, -height/2.0), point2=(0.0, height/2.0))plateSketch.Line(point1=(0.0, height/2.0), point2=(length, height/2.0))plateSketch.Line(point1=(length, height/2.0), point2=(length, -height/2.0))plateSketch.Line(point1=(length, -height/2.0), point2=(0.0, -height/2.0))platePart = plateModel.Part(dimensionality=TWO_D_PLANAR, name="plate", type=DEFORMABLE_BODY)platePart.BaseShell(sketch=plateModel.sketches["plateProfile"])plateSketch.unsetPrimaryObject()del plateSketch# Part for crack geometrycrackSketch = plateModel.ConstrainedSketch(name="crackProfile",sheetSize=height)crackSketch.Line(point1=(0.0, crack_y), point2=(crack_len, crack_y))crackPart = plateModel.Part(dimensionality=TWO_D_PLANAR, name="crack", type=DEFORMABLE_BODY)crackPart.BaseWire(sketch=plateModel.sketches["crackProfile"])crackSketch.unsetPrimaryObject()del crackSketch#---------------------------------------------------------------------------------# === Create geometry sets ===platePart.Set(faces=platePart.faces[:], name="All")e1 = platePart.edges.findAt(( (length/2.0, -height/2.0, 0.0), ))platePart.Set(edges=e1, name="bottom")e1 = platePart.edges.findAt(( (length/2.0, +height/2.0, 0.0), ))platePart.Set(edges=e1, name="top")#---------------------------------------------------------------------------------# === Define material and section properties ===plateMatl = plateModel.Material(name="elas")plateMatl.Elastic(table=((YM, NU), ))plateMatl.MaxpsDamageInitiation(table=((MAXPS, ), ), tolerance=DTOL)plateModel.HomogeneousSolidSection(material="elas", name="solid", thickness=width)#---------------------------------------------------------------------------------# === Assign section and orientation ===#platePart.MaterialOrientation(fieldName="", localCsys=None, orientationType=GLOBAL, region=# platePart.sets["All"], stackDirection=STACK_3)platePart.SectionAssignment(region=platePart.sets["All"], sectionName="solid")#---------------------------------------------------------------------------------# === Assign mesh controls and mesh plate ===# Element typeplatePart.setElementType(elemTypes=(ElemType(elemCode=CPE4, elemLibrary=STANDARD), ElemType(elemCode=CPE3, elemLibrary=STANDARD)), regions=platePart.sets["All"])# Mesh techniqueplatePart.setMeshControls(elemShape=QUAD, regions=platePart.faces[:], technique=STRUCTURED)# Seed meshex = 30; ey = 60;e1 = platePart.edges.findAt(( (length/2.0, -height/2.0, 0.0), ))platePart.seedEdgeByNumber(edges=e1, number=ex)e1 = platePart.edges.findAt(( (length, 0.0, 0.0), ))platePart.seedEdgeByNumber(edges=e1, number=ey)platePart.generateMesh()# === End part ===#---------------------------------------------------------------------------------# === Assemble ===plateModel.rootAssembly.DatumCsysByDefault(CARTESIAN)plateModel.rootAssembly.Instance(dependent=ON, name="plate_1", part=platePart)plateModel.rootAssembly.Instance(dependent=ON, name="crack_1", part=crackPart)# Reference points for displacement applicationrp_db = plateModel.rootAssembly.ReferencePoint(point=(length, -height/2.0, 0.0))plateModel.rootAssembly.features.changeKey(fromName="RP-1", toName="db")#---------------------------------------------------------------------------------# === Create assembly sets ===v1 = (plateModel.rootAssembly.referencePoints[rp_db.id], )plateModel.rootAssembly.Set(name="bdisp", referencePoints=v1)# === End assembly ===#---------------------------------------------------------------------------------# === Create constraint equation ===plateModel.Equation(name="ce_bot", terms=((1.0, "plate_1.bottom", 1), (-1.0, "bdisp", 1)))#---------------------------------------------------------------------------------# === Create step ===plateModel.StaticStep(timePeriod=0.483,initialInc=0.01, maxInc=0.05, maxNumInc=10000, minInc=1e-012, name="Static", nlgeom=ON, previous="Initial")plateModel.steps["Static"].control.setValues( allowPropagation=OFF, resetDefaultValues=OFF, timeIncrementation=(4.0, 8.0, 9.0, 16.0, 10.0, 4.0, 12.0, 20.0, 6.0, 3.0, 50.0))plateModel.fieldOutputRequests["F-Output-1"].setValues( variables=("S", "LE", "U", "RF", "PHILSM", "STATUSXFEM"))plateModel.HistoryOutputRequest(createStepName="Static", name="H-Output-2", rebar=EXCLUDE, region= mdb.models[modelName].rootAssembly.sets["bdisp"], sectionPoints= DEFAULT, variables=("U1", "RF1"))#---------------------------------------------------------------------------------# === Apply boundary conditions ===plateModel.TabularAmplitude(name="Amp-1", timeSpan=STEP, smooth=SOLVER_DEFAULT, data=((0.0, 0.0), (1, 1.0)))plateModel.DisplacementBC( createStepName= "Static", distributionType=UNIFORM, fieldName="", fixed=OFF, localCsys=None , name="rp",amplitude="Amp-1", region=plateModel.rootAssembly.sets["bdisp"], u1=0.01, u2= UNSET, ur3=UNSET)plateModel.DisplacementBC(amplitude=UNSET, createStepName= "Static", distributionType=UNIFORM, fieldName="", fixed=OFF, localCsys=None , name="bot", region=plateModel.rootAssembly.instances["plate_1"].sets["bottom"] , u1=UNSET, u2=0, ur3=UNSET)plateModel.DisplacementBC( createStepName= "Static", distributionType=UNIFORM, fieldName="", fixed=OFF, localCsys=None , name="top",amplitude="Amp-1", region=plateModel.rootAssembly.instances["plate_1"].sets["top"] , u1=-0.01, u2=0, ur3=UNSET)#---------------------------------------------------------------------------------# === Define enrichment and initial crack ===edges = plateModel.rootAssembly.instances["crack_1"].edgese1 = edges.findAt(( (crack_len/2.0, crack_y, 0.0), ))crackLocation = Region(edges=e1)plateModel.ContactProperty("contact")plateModel.rootAssembly.engineeringFeatures.XFEMCrack( name="enr1", crackDomain=plateModel.rootAssembly.instances["plate_1"].sets["All"] , interactionProperty="contact", crackLocation=crackLocation) plateModel.interactionProperties["contact"].FractureCriterion( initTable=((GI, GI, GI, 1.0, 1.0, 1.0), ), mixedModeBehavior=POWER, tolerance=0.2, viscosity=5e-05)plateModel.rootAssembly.engineeringFeatures.cracks["enr1"].setValues(interactionProperty="contact")#---------------------------------------------------------------------------------# === Create Job === mdb.Job(model=modelName, name="modeII_lefm_xfem_cpe4", description="Crack propagation in a plate under mixed-mode loading (XFEM)")#---------------------------------------------------------------------------------
搞定。。。
推薦閱讀：
※精通（或者說學好？）abaqus可以做什麼工作？
※如何修改abaqus並行計算的默認設置？
※abaqus批處理？
※如何使用Abaqus輸入隨時間變化的材料屬性，是否需要編寫用戶程序？