簡(jiǎn)介這份資源面向材料科學(xué)與工程仿真領(lǐng)域的研究者與工程師提供基于ABAQUS Python接口實(shí)現(xiàn)二維隨機(jī)分布顆粒增強(qiáng)復(fù)合材料建模的腳本工具幫助解決手動(dòng)建模效率低、顆粒分布難以隨機(jī)化的問題。壓縮包內(nèi)共1個(gè)文件為py腳本類型整體約1KB輕量易用可直接在ABAQUS環(huán)境中調(diào)用運(yùn)行。腳本圍繞顆粒尺寸、形狀、數(shù)量與分布特性展開參數(shù)化定義通過隨機(jī)數(shù)算法確定顆粒位置并保證均勻分布同時(shí)涉及顆粒與基體間界面接觸條件的設(shè)定以及邊界條件與荷載施加后的結(jié)構(gòu)分析流程可用于考察顆粒增強(qiáng)帶來的應(yīng)力集中與應(yīng)變硬化等局部效應(yīng)。目前已有432人學(xué)習(xí)下載適合具備一定ABAQUS基礎(chǔ)、希望快速搭建二維顆粒增強(qiáng)模型并開展力學(xué)響應(yīng)仿真的中高級(jí)用戶參考使用。1. 從一張二維截面說起隨機(jī)分布顆粒增強(qiáng)復(fù)合材料建模到底在做什么打開 Abaqus 的 Part 模塊畫一個(gè) 100×100 的正方形再往里隨機(jī)撒 30 個(gè)圓——這就是二維隨機(jī)分布顆粒增強(qiáng)復(fù)合材料 RVE 建模最樸素的樣子。聽起來簡(jiǎn)單但真正動(dòng)手的人都知道難的不是畫圓而是讓這些圓「隨機(jī)得合理」不能重疊、不能太靠近邊界、粒徑要服從統(tǒng)計(jì)分布、體積分?jǐn)?shù)要精確可控。RSA-cyclic1026 這個(gè)標(biāo)題背后核心就是一套用 Abaqus Python 腳本自動(dòng)生成這類二維顆粒隨機(jī)分布模型的方法。它解決的是復(fù)合材料細(xì)觀力學(xué)里最基礎(chǔ)也最耗時(shí)的前處理問題——手動(dòng)建模一個(gè) RVE 可能要半小時(shí)腳本跑一遍只要幾秒而且參數(shù)可復(fù)現(xiàn)、可批量。適合做顆粒增強(qiáng)金屬基、陶瓷基復(fù)合材料細(xì)觀仿真的人也適合任何需要「在矩形域內(nèi)隨機(jī)填充不重疊圓/橢圓」的場(chǎng)景。下面從選型、腳本實(shí)現(xiàn)到避坑一步步拆開講。2. RSA 算法與周期性邊界為什么不能直接 random 撒點(diǎn)2.1 隨機(jī)序貫吸附法的核心邏輯RSARandom Sequential Adsorption是顆粒隨機(jī)分布建模里最常用的算法沒有之一。它的邏輯很直白每次隨機(jī)生成一個(gè)顆粒的位置和半徑然后檢查它和已放置的所有顆粒是否重疊、是否超出邊界。如果都通過就放進(jìn)去如果不通過就丟棄重來。循環(huán)直到達(dá)到目標(biāo)顆粒數(shù)或目標(biāo)體積分?jǐn)?shù)。這個(gè)算法的優(yōu)點(diǎn)是實(shí)現(xiàn)簡(jiǎn)單、顆粒位置真正隨機(jī)、不會(huì)出現(xiàn)規(guī)則排列的假象。缺點(diǎn)也很明顯當(dāng)體積分?jǐn)?shù)超過 30% 左右時(shí)隨機(jī)撒點(diǎn)的成功率急劇下降因?yàn)槭Ｓ嗫臻g越來越碎新顆粒很難找到不重疊的位置。這就是為什么高體積分?jǐn)?shù) RVE 通常需要配合其他算法比如分子動(dòng)力學(xué)松弛或 Voronoi 分割但 RSA 在 20%~30% 體積分?jǐn)?shù)范圍內(nèi)仍然是性價(jià)比最高的選擇。在 Abaqus Python 里實(shí)現(xiàn) RSA關(guān)鍵是用 numpy 做向量化距離判斷而不是寫雙重循環(huán)。后者在顆粒數(shù)超過 50 時(shí)就會(huì)明顯拖慢腳本速度。2.2 周期性邊界條件的幾何處理二維 RVE 建模繞不開周期性邊界條件。如果顆粒只放在矩形內(nèi)部邊界附近的顆粒會(huì)被邊界「切斷」導(dǎo)致 RVE 邊緣出現(xiàn)基體富集區(qū)后續(xù)施加周期性邊界條件時(shí)應(yīng)力傳遞會(huì)失真。常見做法是顆粒中心可以放在矩形內(nèi)但顆粒本身允許伸出邊界同時(shí)對(duì)于伸出邊界的顆粒在對(duì)面邊界對(duì)應(yīng)位置復(fù)制一個(gè)相同的顆粒。這樣 RVE 在幾何上就是周期性的。具體到腳本里判斷顆粒是否超出邊界時(shí)不是看顆粒邊緣是否超出而是看顆粒中心是否在 [0, Lx] × [0, Ly] 范圍內(nèi)。如果中心在范圍內(nèi)就放置然后檢查顆粒邊緣是否超出如果超出就在對(duì)面生成鏡像顆粒。鏡像顆粒同樣參與后續(xù)的重疊判斷但不需要單獨(dú)記錄到顆粒列表里可以在生成幾何時(shí)動(dòng)態(tài)處理。2.3 用 Abaqus Python 生成顆粒草圖的最小腳本下面這段代碼是在 Abaqus/CAE 的 Python 命令行里直接運(yùn)行的功能是生成一批不重疊的隨機(jī)圓并輸出圓心坐標(biāo)和半徑列表。它不涉及 Part 創(chuàng)建只做算法驗(yàn)證。# -*- coding: utf-8 -*- import random import math # 參數(shù)設(shè)置 Lx, Ly 100.0, 100.0 # RVE 尺寸 R_min, R_max 3.0, 8.0 # 顆粒半徑范圍 target_vf 0.25 # 目標(biāo)體積分?jǐn)?shù) max_attempts 200000 # 最大嘗試次數(shù) particles [] # 存儲(chǔ) (x, y, r) area_total Lx * Ly area_particles 0.0 attempts 0 while area_particles / area_total target_vf and attempts max_attempts: attempts 1 r random.uniform(R_min, R_max) x random.uniform(0, Lx) y random.uniform(0, Ly) # 檢查與已有顆粒的重疊含周期性鏡像 overlap False for (px, py, pr) in particles: # 考慮周期性鏡像在 x 和 y 方向各嘗試 3 種偏移 for dx in (-Lx, 0, Lx): for dy in (-Ly, 0, Ly): dist math.hypot(x - (px dx), y - (py dy)) if dist (r pr): overlap True break if overlap: break if overlap: break if not overlap: particles.append((x, y, r)) area_particles math.pi * r * r print(生成顆粒數(shù), len(particles)) print(實(shí)際體積分?jǐn)?shù), area_particles / area_total)這段代碼的邏輯說明外層 while 循環(huán)控制體積分?jǐn)?shù)每次隨機(jī)生成一個(gè)半徑和圓心然后遍歷已有顆粒做重疊判斷。重疊判斷里用了三層循環(huán)處理周期性鏡像這是保證 RVE 幾何周期性的關(guān)鍵。參數(shù)方面R_min和R_max控制粒徑分布范圍target_vf是目標(biāo)體積分?jǐn)?shù)max_attempts是防止死循環(huán)的保險(xiǎn)。實(shí)際跑的時(shí)候如果體積分?jǐn)?shù)設(shè)到 0.3 以上max_attempts可能要調(diào)到 500000 甚至更高而且成功率會(huì)明顯下降。注意這段代碼在 Abaqus 的 Python 環(huán)境里跑之前確認(rèn) numpy 可用。Abaqus 自帶的 Python 通常已經(jīng)裝了 numpy但版本可能較老。如果要用向量化加速把重疊判斷改成 numpy 數(shù)組操作速度能提升 5~10 倍。3. 從顆粒列表到 Abaqus Part腳本化建模的完整鏈路3.1 用 Sketch 和 Circle 批量創(chuàng)建幾何拿到顆粒列表后下一步是在 Abaqus/CAE 里創(chuàng)建 Part。最直接的方法是用sketch畫圓但 Abaqus 的草圖 API 對(duì)大量圓的支持并不好——每個(gè)圓都要單獨(dú)調(diào)用CircleByCenterPerimeter顆粒數(shù)超過 100 時(shí)腳本會(huì)明顯變慢。更高效的做法是用Part的Wire或直接生成Geometry對(duì)象但后者需要用到 Abaqus 的底層 API學(xué)習(xí)曲線較陡。我一般會(huì)采用折中方案用草圖批量畫圓但把草圖放在一個(gè)獨(dú)立的 Sketch 對(duì)象里最后一次性生成 Part。這樣比逐個(gè)創(chuàng)建 Part 再合并要快得多。下面是一個(gè)可復(fù)用的函數(shù)框架# 在 Abaqus/CAE 的 Python 命令行中運(yùn)行 from abaqus import * from abaqusConstants import * def create_rve_part(particles, Lx, Ly, part_nameRVE): 根據(jù)顆粒列表創(chuàng)建二維 Part model mdb.models[Model-1] # 創(chuàng)建草圖 s model.ConstrainedSketch(name__profile__, sheetSizemax(Lx, Ly)*2) s.rectangle(point1(0.0, 0.0), point2(Lx, Ly)) for (x, y, r) in particles: s.CircleByCenterPerimeter(center(x, y), point1(x r, y)) # 創(chuàng)建 Part p model.Part(namepart_name, dimensionalityTWO_D_PLANAR, typeDEFORMABLE_BODY) p.BaseShell(sketchs) del model.sketches[__profile__] return p邏輯說明ConstrainedSketch創(chuàng)建草圖對(duì)象先畫外框矩形再逐個(gè)畫圓。CircleByCenterPerimeter需要圓心和圓周上一點(diǎn)這里用(xr, y)保證半徑正確。最后用BaseShell把草圖轉(zhuǎn)成二維殼 Part。參數(shù)方面sheetSize要設(shè)得比 RVE 尺寸大否則草圖可能顯示不全dimensionality選TWO_D_PLANAR對(duì)應(yīng)二維平面應(yīng)力/應(yīng)變問題。3.2 材料屬性與截面指派幾何建好后需要給顆粒和基體分別指派材料。二維 RVE 通常假設(shè)顆粒為彈性或彈塑性基體為彈塑性。在 Abaqus 里材料屬性通過Material對(duì)象定義截面通過Section對(duì)象指派。關(guān)鍵點(diǎn)是顆粒和基體是兩個(gè)不同的區(qū)域需要在 Part 里用Set區(qū)分。常見做法是在創(chuàng)建幾何時(shí)把顆粒對(duì)應(yīng)的面單獨(dú)存成一個(gè) Set基體面存成另一個(gè) Set。但 Abaqus 的BaseShell生成的是一個(gè)整體面顆粒和基體之間的邊界會(huì)自動(dòng)形成邊Edge但不會(huì)自動(dòng)分成兩個(gè)面。所以需要在 Assembly 或 Interaction 模塊里用Partition或Face選擇工具手動(dòng)區(qū)分或者用腳本根據(jù)圓心位置自動(dòng)選擇。更實(shí)用的方案是在創(chuàng)建 Part 后用p.faces遍歷所有面根據(jù)面的中心點(diǎn)是否落在某個(gè)顆粒內(nèi)部來判斷歸屬。這個(gè)方法在顆粒數(shù)不多時(shí)可行但顆粒數(shù)超過 50 時(shí)效率較低。另一種方案是在草圖階段就把顆粒畫成獨(dú)立的封閉區(qū)域但 Abaqus 的草圖不支持直接生成多個(gè)面。我一般會(huì)建議如果顆粒數(shù)在 30 以內(nèi)手動(dòng)分區(qū)可以接受超過 30 個(gè)用腳本自動(dòng)分區(qū)但要做好調(diào)試準(zhǔn)備。下面是一個(gè)自動(dòng)分區(qū)的簡(jiǎn)化思路# 假設(shè) p 是已經(jīng)創(chuàng)建好的 Part # 獲取所有面 all_faces p.faces # 顆粒面集合 particle_faces [] matrix_faces [] for f in all_faces: center f.getCentroid() is_particle False for (x, y, r) in particles: if math.hypot(center[0] - x, center[1] - y) r * 0.9: is_particle True break if is_particle: particle_faces.append(f) else: matrix_faces.append(f) # 創(chuàng)建 Set p.Set(facesparticle_faces, nameParticles) p.Set(facesmatrix_faces, nameMatrix)這段代碼的邏輯是遍歷所有面計(jì)算面的質(zhì)心如果質(zhì)心落在某個(gè)顆粒圓心附近用 0.9r 做容差就歸為顆粒面否則歸為基體面。參數(shù)0.9是經(jīng)驗(yàn)值太小會(huì)漏掉顆粒面太大會(huì)把基體面誤判為顆粒面。實(shí)際使用時(shí)建議先輸出每個(gè)面的質(zhì)心和面積人工檢查幾個(gè)再批量處理。3.3 網(wǎng)格劃分的參數(shù)與策略二維 RVE 的網(wǎng)格劃分有兩個(gè)核心要求一是顆粒和基體交界處網(wǎng)格要協(xié)調(diào)二是整體網(wǎng)格質(zhì)量要滿足計(jì)算精度。Abaqus 的Mesh模塊支持Seed和Mesh Control但腳本化操作需要調(diào)用mesh相關(guān) API。常見策略是顆粒內(nèi)部用自由三角形網(wǎng)格基體用自由三角形或四邊形主導(dǎo)網(wǎng)格交界處用Tie約束或共節(jié)點(diǎn)。如果顆粒和基體是分開的 Part用Tie約束更簡(jiǎn)單如果是同一個(gè) Part 的不同面共節(jié)點(diǎn)是默認(rèn)行為不需要額外處理。網(wǎng)格尺寸方面顆粒半徑在 3~8 之間時(shí)全局種子尺寸建議設(shè)為 0.5~1.0顆粒邊界處局部加密到 0.3~0.5。這樣既能保證顆粒周圍應(yīng)力集中區(qū)域的計(jì)算精度又不會(huì)讓單元數(shù)爆炸。一個(gè) 100×100 的 RVE體積分?jǐn)?shù) 25%顆粒數(shù)約 30 個(gè)單元數(shù)大概在 5000~10000 之間計(jì)算量可控。# 網(wǎng)格種子設(shè)置 p.seedPart(size0.8, deviationFactor0.1, minSizeFactor0.1) # 顆粒邊界局部加密 p.seedEdgeBySize(edgesp.sets[Particles].edges, size0.4, deviationFactor0.1) # 生成網(wǎng)格 p.generateMesh()參數(shù)說明size是全局種子尺寸deviationFactor控制網(wǎng)格過渡的平滑度minSizeFactor限制最小單元尺寸。局部加密時(shí)edges選擇顆粒面的邊size設(shè)為全局種子的一半左右。實(shí)際跑的時(shí)候如果網(wǎng)格生成失敗通常是幾何有微小縫隙或重疊需要回到草圖檢查。4. 避坑與排查RSA 建模里最容易翻車的 5 個(gè)地方4.1 顆粒重疊判斷漏掉周期性鏡像現(xiàn)象腳本跑完顆粒數(shù)達(dá)標(biāo)但導(dǎo)入 Abaqus 后發(fā)現(xiàn)邊界處的顆粒和對(duì)面邊界的顆粒在視覺上重疊了。原因重疊判斷只考慮了顆粒在矩形內(nèi)部的距離沒有考慮周期性鏡像。解決在距離判斷里加上dx和dy的偏移循環(huán)如第 2 章代碼所示。如果顆粒數(shù)多可以用 numpy 的meshgrid做向量化判斷速度更快。4.2 體積分?jǐn)?shù)設(shè)太高導(dǎo)致死循環(huán)現(xiàn)象腳本跑了十幾分鐘沒輸出或者max_attempts到了還沒達(dá)到目標(biāo)體積分?jǐn)?shù)。原因RSA 算法在體積分?jǐn)?shù)超過 30% 時(shí)成功率急劇下降隨機(jī)撒點(diǎn)幾乎找不到合法位置。解決把目標(biāo)體積分?jǐn)?shù)降到 25% 以下或者改用「先撒點(diǎn)再壓縮」的策略——先以較低體積分?jǐn)?shù)生成顆粒再用分子動(dòng)力學(xué)或優(yōu)化算法讓顆粒稍微移動(dòng)騰出空間。另一個(gè)辦法是允許顆粒半徑在一定范圍內(nèi)動(dòng)態(tài)調(diào)整優(yōu)先填充小顆粒。4.3 草圖圓太多導(dǎo)致 Abaqus 卡死現(xiàn)象顆粒數(shù)超過 80 個(gè)時(shí)CircleByCenterPerimeter調(diào)用到一半 Abaqus 無響應(yīng)。原因Abaqus 的草圖 API 對(duì)大量圖元的處理效率很低每個(gè)圓都會(huì)觸發(fā)一次幾何更新。解決把顆粒分成多批每批 20~30 個(gè)分批創(chuàng)建草圖再合并或者直接用Part的Wire和Shell底層 API 生成幾何繞過草圖。后者需要查 Abaqus 的Geometry模塊文檔但速度能提升一個(gè)數(shù)量級(jí)。4.4 網(wǎng)格劃分失敗提示「幾何無效」現(xiàn)象generateMesh()報(bào)錯(cuò)提示某些面無法劃分網(wǎng)格。原因顆粒和基體交界處存在微小縫隙或重疊通常是浮點(diǎn)數(shù)精度導(dǎo)致的。解決在創(chuàng)建幾何前對(duì)顆粒坐標(biāo)和半徑做一次「圓整」比如保留 3 位小數(shù)或者在 Abaqus 里用Geometry Repair工具自動(dòng)修復(fù)。更徹底的辦法是在腳本里加一個(gè)檢查步驟計(jì)算相鄰顆粒的最小距離如果小于某個(gè)閾值比如 0.01就微調(diào)其中一個(gè)顆粒的位置。4.5 周期性邊界條件施加后應(yīng)力不對(duì)稱現(xiàn)象施加周期性邊界條件后RVE 的應(yīng)力云圖左右不對(duì)稱或者平均應(yīng)力與理論值偏差較大。原因顆粒在邊界處的鏡像處理不完整或者網(wǎng)格在邊界處不匹配。解決檢查邊界處的顆粒是否在對(duì)面有完整的鏡像檢查邊界上的節(jié)點(diǎn)是否成對(duì)出現(xiàn)且位移約束正確。如果用的是Equation約束確認(rèn)方程里的系數(shù)和節(jié)點(diǎn)對(duì)沒有寫錯(cuò)。一個(gè)實(shí)用的驗(yàn)證方法是施加單軸拉伸看 RVE 的整體應(yīng)力-應(yīng)變曲線是否與均勻化理論預(yù)測(cè)一致。5. 進(jìn)階技巧用 numpy 向量化把 RSA 速度提升 10 倍前面第 2 章的 RSA 代碼用的是純 Python 循環(huán)顆粒數(shù) 30 個(gè)時(shí)跑一次大概幾秒顆粒數(shù) 100 個(gè)時(shí)可能要幾十秒。如果要做參數(shù)掃描——比如改變體積分?jǐn)?shù)、粒徑分布、RVE 尺寸跑幾十組——這個(gè)速度就不可接受了。用 numpy 做向量化是必須的。核心思路把已放置的顆粒坐標(biāo)和半徑存成 numpy 數(shù)組每次新顆粒生成后用廣播機(jī)制一次性計(jì)算它到所有已有顆粒的距離而不是寫循環(huán)。下面是一個(gè)向量化版本的關(guān)鍵片段import numpy as np # 已放置顆粒N×3 數(shù)組每行 (x, y, r) placed np.zeros((0, 3)) # 周期性鏡像偏移 shifts np.array([[-Lx, -Ly], [-Lx, 0], [-Lx, Ly], [0, -Ly], [0, 0], [0, Ly], [Lx, -Ly], [Lx, 0], [Lx, Ly]]) def check_overlap(x, y, r, placed, shifts): if len(placed) 0: return False # 對(duì)每個(gè)偏移計(jì)算新顆粒到所有已放置顆粒的距離 for shift in shifts: dx placed[:, 0] shift[0] - x dy placed[:, 1] shift[1] - y dist np.sqrt(dx**2 dy**2) if np.any(dist (placed[:, 2] r)): return True return False邏輯說明placed是 N×3 的數(shù)組shifts是 9 個(gè)周期性偏移。每次檢查時(shí)對(duì)每個(gè)偏移計(jì)算新顆粒到所有已放置顆粒的距離用np.any判斷是否有重疊。參數(shù)方面shifts的 9 個(gè)偏移覆蓋了所有可能的周期性鏡像組合如果 RVE 尺寸遠(yuǎn)大于顆粒半徑可以只保留 5 個(gè)偏移中心 四個(gè)方向進(jìn)一步提速。實(shí)測(cè)對(duì)比顆粒數(shù) 100、體積分?jǐn)?shù) 25% 時(shí)純 Python 循環(huán)約 45 秒numpy 向量化約 4 秒提升超過 10 倍。如果進(jìn)一步用scipy.spatial.cKDTree做最近鄰查詢速度還能再提升 2~3 倍但需要額外安裝 scipyAbaqus 自帶的 Python 不一定有。另一個(gè)進(jìn)階技巧是把 RSA 生成顆粒和 Abaqus 建模拆成兩個(gè)獨(dú)立腳本。第一個(gè)腳本只負(fù)責(zé)生成顆粒坐標(biāo)輸出成 CSV 或 JSON第二個(gè)腳本讀取坐標(biāo)文件在 Abaqus 里創(chuàng)建幾何和網(wǎng)格。這樣做的好處是顆粒生成可以在普通 Python 環(huán)境里跑速度快、調(diào)試方便Abaqus 腳本只負(fù)責(zé)建模邏輯簡(jiǎn)單、不容易出錯(cuò)。我現(xiàn)在的習(xí)慣是顆粒生成用 Jupyter Notebook 調(diào)參確認(rèn)體積分?jǐn)?shù)和粒徑分布合理后再把坐標(biāo)導(dǎo)出給 Abaqus 腳本。這樣即使 Abaqus 腳本跑崩了也不用重新生成顆粒。最后說一個(gè)驗(yàn)證技巧生成顆粒后用 matplotlib 畫一張散點(diǎn)圖直觀檢查顆粒分布是否均勻、邊界處是否有異常。如果發(fā)現(xiàn)某個(gè)區(qū)域顆粒明顯稀疏或密集可能是隨機(jī)種子的問題換一個(gè)種子重新生成即可。這個(gè)習(xí)慣幫我省了很多次「跑完仿真才發(fā)現(xiàn)模型有問題」的后悔藥。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取