本征邊界積分方程法模擬含流體粒子固體的性能

周吉成， 和東宏， 馬 杭

(1. 上海大學上海市應用數學和力學研究所， 上海200072; 2. 上海大學理學院， 上海200444)

自然界中的許多材料， 如生物組織、飽和巖土、膠體材料、發(fā)泡塑料以及燒結陶瓷等通常都是流體與固體材料的復合體. 在這類含流體夾雜的彈性介質中， 流體粒子與基體之間存在復雜的相互作用， 其形狀、體積分數和內壓等都會對材料的有效力學性能產生影響. 顯然， 這類材料的力學行為不同于含有固體粒子的材料. 近年來流體粒子對材料宏觀有效性能影響規(guī)律的研究受到力學及復合材料研究者的關注[1-3].

對于夾雜粒子和非均質體的研究， 從1957 年Eshelby[4]創(chuàng)造性地提出等效夾雜理論后獲得快速發(fā)展， 出現了大量對夾雜物和非均質體的研究[5-6]， 以及基于Eshelby 等效夾雜和本征應變解的各種理論分析[7-8]和數值計算[9]的研究. 在各種物理問題中， 本征應變解可適用于非協(xié)調熱應變、相變、塑性應變以及殘余應力的固有應變[10]等問題. 通過對Eshelby等效夾雜物的替換， 本征應變還可用來處理固體基體中的異質體如夾雜粒子、孔洞[11]等問題. Eshelby 等效夾雜理論在材料科學、力學、物理、工程等問題中均有廣泛的應用. 近年來， 通過將邊界積分方程與Eshelby 等效夾雜理論相結合， 已有研究者提出了本征應變邊界積分方程的計算模型， 并成功地對含有大量固體粒子的固體進行了二維和三維數值模擬[12-13].

本工作將Eshelby 本征應變和等效夾雜理論與邊界積分方程相結合， 提出了含大量流體粒子固體的本征應變邊界積分方程的計算模型及迭代算法， 拓寬了模型的應用范圍. 通過數值算例， 驗證了計算模型的正確性和計算方法的可行性， 討論了材料參數如流體粒子分布、孔隙率及壓縮系數等對含流體粒子固體材料整體性能的影響.

1 計算模型

1.1 本征應變形式的邊界積分方程

在本模型中， Ω 為求解域即固體基體， 其外邊界為Γ; 基體中流體粒子總數NI所在的區(qū)域為 ΩI(I =1，2···，NI)， 其邊界即粒子與基體的界面 ΓI(ΓI= ΩI∩Ω). 對于問題的求解域Ω，位移和應力邊界積分方程[9，12-13]分別為

式中:

為式(2)域積分的自由項， Ωε為ΩI中ε無窮小域， 其邊界為Γε， xl為場點 x 和源點 y 之間距離的粒子后， 區(qū)域ΩI中的本征應變.

在上述邊界積分方程中， 本征應變及流體粒子的狀態(tài)都是未知的， 需要在數值求解過程中逐步加以確定. 每個粒子的本征應變取決于外加應力(或者應變)、粒子的幾何形狀、粒子和基體的材料參數等因素. 需要指出的是， 邊界積分方程(1)和(2)描述的是彈性均勻介質的應力場和位移場， 要描述含有非均質流體粒子的固體， 必須借助Eshelby 張量進行夾雜物的等效替換.

1.2 Eshelby 張量和等效夾雜理論

依據Eshelby 理論， 考慮二維全空間中的單個粒子， 在粒子域ΩI上存在的本征應變ε0kl滿足

式中: δjl為 Kronecker 符號; μ 和 ν 分別為基體的剪切模量和泊松比. 需要說明的是， 式(5)的適用條件是本征應變在域ΩI內均勻分布， 并利用以下積分恒等式[12-13]，

使得區(qū)域型積分轉化為邊界型積分. 假設流體為線性非黏性可壓縮流體， 其本構關系可表示為

其中外加應變、拘束應變和本征應變三者應變偏量的表達式為

將式(8)和(9)改寫為

式中:

將式(4)代入式(11)， 便得到流體粒子的本征應變與外加應變的關系，

將式(13)改寫成矩陣形式， 即

1.3 局部Eshelby 矩陣

當彈性固體中含有大量流體粒子時， 粒子間存在相互作用， 其大小取決于粒子間距離. 因此， 通過外加應力計算當前粒子的應變時， 除了要考慮遠場荷載， 還需要考慮其他粒子對當前粒子的作用. 準確地說， 這種相互作用產生了兩個因素需要考慮. 首先， 多粒子條件下本征應變不再是常數. 但是計算實踐證明， 當粒子間無量綱距離滿足一定條件時， 本征應變的常數假定能夠近似成立[12-13]. 按本征應變的常數假定， 本模型中本征應變和外加應變都在粒子的幾何中心處進行計算. 第二個需要考慮的因素是， 迭代計算的收斂性受到粒子間相互作用的影響. 為了考慮這一影響， 特別是處理近距離粒子間較強的相互作用， 將計算區(qū)域的粒子分成兩組， 即近場群和遠場群.

如圖1 所示， 以當前粒子I 為中心， 將虛線所圍成的圓內粒子定義為近場群粒子， 其他為遠場群粒子. 假設全空間中近場群粒子數為N， 暫時忽略遠場群粒子的作用， 利用式(2)， 當前粒子的應力邊界積分方程可以表示成

圖1 全空間中當前粒子I 的近場群定義Fig.1 Definition of the near-field group for the current pore I in full space

運用Eshelby 張量和等效夾雜理論以及式(11)， 仿照1.2 節(jié)對單個粒子的推導過程， 將近場群粒子進行離散， 得到本征應變邊界積分方程的矩陣形式

式中: {ε0}={ε01，ε02，··· ，ε0N}T表示近場群粒子的本征應變向量; {ε}={ε1，ε2，··· ，εN}T表示近場群粒子的外加應變向量; 矩陣[S]的表達式為

在矩陣[S]中， 主對角線元素Skk反映了粒子本身的本征應變和外加應變的關系， 非主對角線元素Sjk(j ／=k)表示近場群中k 粒子的外加應變對于j 粒子本征應變的影響. 矩陣[S]的結構與方程(14)中的矩陣S 相似， 其表達式參見式(13). 需要指出的是， 當粒子隨機分布時， 每個粒子的近場群都是不同的. 假設當前粒子的編號為I， 每個粒子的本征應變向量ε0I都可通過其近場群粒子群的矩陣表達式和外加應變向量 {εI}={ε1， ε2，··· ，εN}TI計算得到，

式中: TI是對矩陣[S]求逆和縮并得到的， 稱為局部Eshelby 矩陣. 這樣， 基于本征應變邊界積分方程和局部Eshelby 矩陣的算法， 固體中每個流體粒子的本征應變由以下三部分組成: 第一部分在外加載荷的作用下產生; 第二部分由近場群粒子之間的相互作用而產生， 通過式(18)計算得到; 第三部分由遠場群粒子的弱相互作用產生， 通過迭代來計算[12-13].

2 數值算例

在數值算例中， 所有數值計算結果都通過本征應變邊界積分方程法和子域邊界積分方程法[14]獲得. 利用邊界點法對邊界積分方程進行離散時， 對于奇異邊界單元采用二次移動單元來離散， 對于非奇異邊界單元采用常數單元來離散[15]. 對于粒子的邊界或界面， 本征應變邊界積分方程法和子域邊界積分方程法分別用20 和40 節(jié)點來離散. 在所有的算例中， 泊松比取0.3， 楊氏模量取1. 圓形流體粒子的半徑用r0表示， 橢圓形粒子的短軸與長軸之比定義為縱橫比.

2.1 全空間中的單個圓形流體粒子

對于遠場均布單位壓力(p0=1)下全空間中的單個圓形流體粒子， 流體粒子的內部壓力p、固體基體徑向位移ur和徑向應力σr的理論解[14]為

式(19)～(21)中: E 表示固體基體的彈性模量; K 表示流體粒子的壓縮系數; r0為流體粒子半徑; r 為全空間中粒子外部一點到流體粒子中心的距離.

表1 列出了分別采用本征應變邊界積分方程法和子域邊界積分方程法對單個圓形流體粒子的數值計算結果以及與理論解的對比， 表中負值表示壓縮. 在表2 和3 中分別比較了3 種方法計算得到的基體內一點(x1=2r0， x2=r0)的徑向位移和徑向應力. 從表1～3 可以看出， 兩個數值程序的結果與理論解的對比均具有很好的一致性， 雖然本征應變邊界積分方程法使用了較少的節(jié)點， 但其計算結果比子域邊界積分方程法更接近理論解， 顯示了計算模型的正確性和可行性. 從理論上來說， 由于在本征應變邊界積分方程法中， 邊界未知量不直接出現在系統(tǒng)矩陣中， 因此較少的節(jié)點數不會明顯降低其計算精度.

表1 單個圓形流體粒子內部壓力的對比Table 1 Comparison of inner pressures in a single circular fluid-filled pore

表2 固體基體中一點(x1 =2r0， x2 =r0)的徑向位移對比Table 2 Comparison of radial displacements at the point (x1 =2r0， x2 =r0) in solid media

表3 固體基體中一點(x1 =2r0， x2 =r0)的徑向應力對比Table 3 Comparison of radial stresses at the point (x1 =2r0， x2 =r0) in solid media

2.2 粒子附近處基體的應力

圖2(a)和(b)分別比較了全空間遠場均布雙軸單位壓縮荷載情況下， 2 和4 個圓形流體粒子附近基體的應力. 從圖2 中可以看出， 本征應變邊界積分方程法和子域邊界積分方程法的結果相互吻合， 驗證了本工作提出的計算模型的正確性和可行性.

圖2 全空間2 和4 個圓形粒子雙向受壓時粒子附近的應力Fig.2 Dimensionless stresses at solid media adjacent to the two and four fluid-filled pores under uniform biaxial compression in full space

2.3 含單個流體粒子的代表性體積單元的整體性能

在單軸壓縮荷載條件下， 分別采用本征應變邊界積分方程法和子域邊界積分方程法對含有單個流體粒子(NI= 1)的正方形代表性體積單元(representative volume element， RVE)的整體性能進行計算和比較. 在兩種數值方法中， 代表性體積單元的邊界都用200 個節(jié)點進行離散.

圖3(a)和(b)分別給出了含有單個圓形流體粒子的代表性體積單元的無量綱彈性模量E*/E 及泊松比ν*隨孔隙率的變化. 孔隙率定義為流體粒子與代表性體積單元的面積之比.圖4(a)和(b)分別給出了含有單個橢圓流體粒子的代表性體積單元的無量綱彈性模量E*/E 及泊松比ν*隨孔隙率的變化， 其中E1*和ν1*分別表示加載方向和橢圓長軸都平行于x1方向時的整體性能， E2*和ν2*分別表示加載方向和橢圓短軸都平行于x2方向時的整體性能. 顯然， 含有橢圓粒子時， 代表性體積單元的整體性能表現出各向異性. 圖4 和5 表明兩種計算方法得到的結果吻合良好.

圖3 RVE 中含單個圓形粒子時彈性模量和泊松比與孔隙率的關系Fig.3 Overall elastic modulus and Poisson’s ratio of RVE with a single circular pore as a function of porosity ratio

圖4 RVE 中含單個橢圓粒子時彈性模量和泊松比與孔隙率的關系Fig.4 Overall elastic modulus and Poisson’s ratio of RVE with an elliptical pore as a function of porosity ratio

2.4 含多個流體粒子的代表性體積單元的整體性能

本工作計算分析了代表性體積單元中含有多個規(guī)則分布的流體粒子時的整體性能， 粒子的分布如圖5(a)所示. 根據當前流體粒子位置的不同， 近場群的定義也有所不同(見圖5(b)).圖5(b)中編號1 的情況最常見， 其近場群包含9 個粒子. 而編號2 和3 的粒子分別靠近邊界或角點， 近場群包含粒子數相對較少. 多粒子條件下代表性體積單元邊界的離散、孔隙率的定義與2.3 節(jié)相同.

在保持孔隙率不變的條件下， 粒子數目對代表性體積單元整體性能的影響如圖6 所示. 計算結果表明， 當流體粒子的數目增加到足夠大， 即大約400 時， 若繼續(xù)增加粒子的數目， 計算得到的代表性體積單元的彈性模量E*/E 和泊松比ν*均保持了相對穩(wěn)定. 根據這一結果， 當討論多粒子代表性體積單元的整體性能時， 選取的粒子數NI=1 089.

在單軸壓縮條件下， 含多個圓形流體粒子的代表性體積單元的彈性模量和泊松比隨孔隙率的變化如圖7 所示. 在圖7(b)中， 泊松比隨孔隙率的變化表現出不同的趨勢， 這一現象與流體的性質有關. 靜態(tài)流體內部不存在切應力和切應變， 不存在對于形狀改變的抵抗力， 而體積改變則與壓縮系數有關. 隨著壓縮系數的減小， 流體粒子的力學性能將逐漸接近空洞的力學性能， 所以泊松比隨孔隙率的增加而增大. 當壓縮系數較大時， 隨著孔隙率的增加， 粒子間的相互作用也隨之增大， 使得泊松比反而減小. 這一整體性能的變化趨勢與Kachanov 等[2]和Huang等[14]的計算結果也是一致的.

圖5 RVE 中含有多個規(guī)則分布的圓形流體粒子示意圖及3 種不同情況下近場群的定義Fig.5 Schematics of the RVE with multiple circular fluid-filled pores in regular distribution and the definitions of the near-field group in three different cases

圖6 RVE 中含有多個圓形流體粒子時整體彈性模量和泊松比隨粒子數NI 的變化Fig.6 Overall elastic modulus and Poisson’s ratio of RVE with multiple circular pores as a function of pore number NI

2.5 粒子隨機分布對RVE 整體性能的影響

在實際問題中， 更常見的是固體中流體粒子的參數如位置、形狀、大小和壓縮系數等隨機分布的情況， 這些參數可通過Fortran 語言中的偽隨機函數random 來生成， 圖8(a)中示意地給出了多粒子代表性體積單元中粒子位置隨機分布的一種情況. 對于橢圓粒子， 其傾角θ 如圖8(b)所示. 本工作計算了兩種不同隨機粒子分布時的代表性體積單元的體積模量， 并與流體粒子規(guī)則分布的代表性體積單元的體積模量進行了比較. 第一種為壓縮系數在10-9～100范圍內隨機變化且流體粒子的位置隨機變化; 在此基礎上， 第二種隨機分布增加了粒子大小以及橢圓粒子傾角的隨機變化， 計算中橢圓粒子的縱橫比保持不變(0.6).

圖7 雙向受壓時含多個圓形流體粒子的RVE 的彈性模量和泊松比與孔隙率的關系Fig.7 Overall elastic modulus and Poisson’s ratio of RVE with multiple circular pores as a function of porosity ratio under uniaxial compression

圖8 RVE 中含有多個位置隨機分布的流體粒子示意圖及橢圓粒子傾角θ的定義Fig.8 Schematics of the RVE with multiple fluid-filled pores in random distribution by position and the definition of tilting angle θ of elliptic pore

在雙軸均勻受壓荷載下， 本工作計算了兩種不同隨機粒子分布時的代表性體積單元的體積模量隨孔隙率的變化， 并與粒子規(guī)則分布時的結果進行了比較， 如圖9 所示. 結果表明， 兩種隨機分布時的體積模量分別與壓縮系數K =8.75×10-2和4.55×10-2的規(guī)則分布時的代表性體積單元的體積模量相吻合. 由此可知， 在一定條件下， 隨機分布流體粒子固體的整體性能可以借助于具有規(guī)則分布粒子的固體來進行研究. 然而這一做法的效果是有限的， 因為粒子參數與固體力學性能的關系比較復雜， 需要借助數值計算的手段進行研究， 特別是粒子隨機分布的情況. 本工作提出的計算模型為這一研究提供了有效的工具.

2.6 計算效率

圖10 給出了本征應變邊界積分方程法和子域邊界積分方程法CPU 時間的對比. 可以看出， 隨著粒子數的增加， 本征應變邊界積分方程法計算模型的效率遠高于子域邊界積分方程法， 因為在子域邊界積分方程法中， 系統(tǒng)矩陣的維數隨著粒子數的增加而呈幾何級數增長， 當粒子數超過百位數時， 由于計算機內存的限制， 程序已無法運行. 而在本征應變邊界積分方程法計算模型中， 系統(tǒng)矩陣的維數與粒子數無關， 只有局部Eshelby 矩陣的數量與粒子數成正比，計算量隨著粒子數的增加按算術級數增加， 因此本工作提出的計算模型具有極高的計算效率.

圖9 RVE 中含有多個流體粒子時隨機分布和規(guī)則分布的體積模量隨孔隙率變化的對比Fig.9 Comparison of overall bulk modulus of RVE with multiple pores between random and regular distributions as a function of porosity ratio

圖10 兩種計算模型(本征應變法和子域法)CPU 時間隨粒子數NI 變化的對比Fig.10 Comparison of CPU times for the two computing procedures，the Eigenstrain and the Subdomain， as a function of total pore numbers NI

3 結束語

本工作將Eshelby 本征應變和等效夾雜理論引入邊界積分方程中， 提出了含大量流體粒子固體的本征應變邊界積分方程的計算模型及迭代算法. 通過全空間中近場群的定義和局部Eshelby 矩陣的引入， 解決了粒子間相互作用的問題， 保證了迭代計算的收斂性. 通過數值算例的比較， 不僅驗證了所提出計算模型的正確性和計算方法的可行性， 也證明計算模型具有很高的計算效率， 為含有大量流體粒子固體的數值模擬提供了新的有效的研究手段.