邱瑩宇，郭 然

(昆明理工大學(xué) 建筑工程學(xué)院，昆明 650504)

頁巖是沉積巖中的一種，具有十分明顯的層理構(gòu)造[1]。頁巖氣是指以吸附或游離且有時(shí)還帶有流體相的形態(tài)賦存于泥頁巖中的非常規(guī)天然氣[2]。頁巖氣在中國資源儲(chǔ)量豐富，地域分布廣泛。頁巖氣開采能緩解我國常規(guī)油氣產(chǎn)量不足、煤化石燃料引起環(huán)境污染等問題，已成為中國綠色能源開發(fā)的重要領(lǐng)域[3]。中國頁巖氣埋藏深，賦存條件差，自然豐度低[2]，高效開發(fā)理論與產(chǎn)能評(píng)價(jià)處于起步階段，因此，高效開采面臨更多的困難和挑戰(zhàn)[4]。而無論是在用水力壓裂手段[1]開采頁巖氣時(shí)，還是處理頁巖中含有流體的大量微小孔洞時(shí)，不可避免地會(huì)出現(xiàn)流固耦合的數(shù)值模擬問題。因此，帶有流體的固體單元的研究具有重要意義。

卞學(xué)鐄[5]在1964年創(chuàng)建了假設(shè)應(yīng)力有限元方法，即應(yīng)力雜交元。最初的應(yīng)力雜交元是通過修正余能泛函建立的，假定單元中存在滿足平衡方程的應(yīng)力場，單元邊界滿足力的平衡條件，并且通過對(duì)界面節(jié)點(diǎn)位移進(jìn)行插值獲得邊界位移。1988年下半年，Accorsi根據(jù)最小勢(shì)能原理推導(dǎo)了一個(gè)新單元[6-7]。2016年，黃永霞等[8]進(jìn)行了顆粒增強(qiáng)復(fù)合材料有效模量的Voronoi單元有限元法分析。2017年，Zhang等[9]使用Voronoi單元法提出了含孔洞單元的修正余能泛函形式，在少量單元的情況下，獲得了在內(nèi)部壓力作用下的應(yīng)力分布。2019年，韓寧等[10]模擬了顆粒夾雜復(fù)合材料的力學(xué)性能,分析了網(wǎng)格劃分對(duì)復(fù)合材料力學(xué)分析精度的影響。同年，Han等[11]提出了一種新的單元，其裂紋從內(nèi)孔的邊緣開始，并且在單元內(nèi)和單元間擴(kuò)展。他們采用的Voronoi單元網(wǎng)格劃分方法在處理形狀簡單和尺寸較小的微結(jié)構(gòu)方面具有優(yōu)勢(shì)，但當(dāng)孔體積分?jǐn)?shù)相對(duì)較大或彼此之間的距離非常接近時(shí)(如圖1所示)，將很難進(jìn)行網(wǎng)格化。

圖1 頁巖中流體微結(jié)構(gòu)的分布情況(藍(lán)色代表流體，綠色代表頁巖)Fig. 1 Distribution of fluid microstructure in shale (blue represents fluid, and green represents shale)

針對(duì)大量含流體孔洞的頁巖問題，筆者提出了一種采用均勻網(wǎng)格并滿足流固界面壓力平衡的新的流固耦合單元fluid-structure coupling stress hybrid element(FCSHE)。該單元可以進(jìn)一步推廣，解決一般性的流固耦合問題。

通過對(duì)簡單四單元模型FCSHE與MARC模型的應(yīng)力分布情況以及等效彈性模量的對(duì)比驗(yàn)證，證明單元的有效性后，建立多個(gè)模型，分別通過改變形狀、體分比、空間分布位置和半徑，研究了等效彈性模量變化的影響因素[12-16]。

1 含流體相的應(yīng)力雜交單元(FCSHE)

雜交元是基于最小余能原理的有限元模型，其中應(yīng)力函數(shù)為場變量，余能泛函如式(1)所示[5]。

(1)

為了滿足表面力的邊界條件

(2)

(3)

式中：n+和n-分別是相鄰單元邊的法向向量，使用拉格朗日乘子法引入2個(gè)約束，以獲得修正余能泛函[2]

(4)

式中：?Ve是所有單元的邊界，?Ωte是給定面力的邊界，u是單元邊上的位移。

如圖2所示，在流體和固體的界面處，滿足表面力的平衡條件

nfs·(σs-σf)=0，

(5)

式中：σs是固體側(cè)應(yīng)力，σf是流體側(cè)應(yīng)力，nfs是流體側(cè)指向固體側(cè)的交界處的法向矢量。

圖2 流固交界處的典型單元Fig. 2 Typical element at fluid-solid junction

引入該約束條件得到修正余能泛函

(6)

式中：nes是元素的外部固體邊界上的外部法線向量；nef是流體邊界上的是從流體側(cè)指向固體側(cè)的界面處的法線向量；σ是在單元域中滿足平衡的應(yīng)力場；ues是在單元固體邊界上的位移；uef是在流體邊界上的位移。

在固體部分Ωes中，將自平衡應(yīng)力場引入單元中[5]：

σs=Psβs。

(7)

式中：Ps為固體部分的插值矩陣，βs為固體部分的待求系數(shù)列陣。

在流體部分Ωef中，考慮液體壓強(qiáng)為一個(gè)常數(shù)，此時(shí)，σx=σy，τxy=0，所以

(8)

(9)

式中：x和y是構(gòu)造出的應(yīng)力函數(shù)的自變量；若為二階函數(shù)，

(10)

注意，在靜態(tài)分析中，頁巖孔中流體的壓力場通常是恒定的，因此以下所有模型均采用方程式(8)。

引入位移插值功能后，邊界位移可以用節(jié)點(diǎn)位移表示：

單元固體側(cè)的外部邊界?Ωes

ues=Lqes，

(11)

單元流體側(cè)的外部邊界?Ωef

uef=Lqef，

(12)

流固交界面?Ωi

ui=Lqi。

(13)

式中：qes是固體界面的節(jié)點(diǎn)位移；qef是流體界面的節(jié)點(diǎn)位移；qi是流體和固體交界面的節(jié)點(diǎn)位移；L是插值函數(shù)。

將方程(7)(8)和(11)～(13)代入方程(6)，得到離散系統(tǒng)修正余能泛函的弱形式。根據(jù)修正余能泛函的平衡條件：

(14)

(15)

可以得到運(yùn)動(dòng)關(guān)系的弱表達(dá)形式

(16)

其中，

(17)

(18)

(19)

(20)

(21)

(22)

式中：Ωs是固體區(qū)域，Ωf是流體區(qū)域，Ωis是流固交界面上固體部分的邊界，Ωif是流固交界面上流體部分的邊界，nis是流固交界面上固體部分的法向向量，nif是流固交界面上流體部分的法向向量，Lis是流固交界面上固體部分的插值函數(shù)，Lif是流固交界面上流體部分的插值函數(shù)。

從上面可以看出，應(yīng)力參數(shù)在單元中的位移表達(dá)式為

β=H-1GQ，

(23)

其中，

(24)

(25)

(26)

(27)

節(jié)點(diǎn)位移qes、qef和qi關(guān)于系統(tǒng)Πsf的總能量的一階偏導(dǎo)可以表示為

(28)

(29)

(30)

隨即可以得到每個(gè)單元中的運(yùn)動(dòng)關(guān)系弱表達(dá)形式

(31)

(32)

其中，

(33)

將方程式(23)和(32)組合在一起，得到表達(dá)式(34)來求解廣義位移：

(34)

得到求解廣義位移的方程組

(35)

其中，

(36)

方程(35)的剛度矩陣為

(37)

式中：Ke表示單元的剛度，e表示所在單元。

2 算 例

2.1 驗(yàn)證模型

在此算例中，建立了一個(gè)帶孔的正方形板模型，分析了其在平面應(yīng)力狀態(tài)下的線彈性靜力學(xué)行為。正方形板邊長為0.40 mm，在其中心有一個(gè)半徑為0.05 mm的孔，并充滿了水(不考慮流動(dòng))。由于MARC軟件中沒有液體單元，因此通過在孔中施加壓力(從FCSHE獲取)來表示含水孔。FCSHE模型由4個(gè)單元組成，MARC模型由11 482個(gè)單元組成。邊界條件如圖3(a)所示，水平位移為0.15 μm。在模型仿真分析中使用的材料屬性如下：

固體：彈性模量Es=72.0 GPa，泊松比νs= 0.33。

流體：彈性模量Ef=2.18 GPa，泊松比νf= 0.50。

由于沒有流固耦合材料的解析解，所以采用與MARC結(jié)果作對(duì)比的形式來驗(yàn)證有效性。圖3中的3幅圖分別是工況圖(a)、由FCSHE和MARC得到的模型(b)和(c)，圖4～6為FCSHE和MARC分別得到的σx、σy和τxy應(yīng)力云圖。一條線上的路徑和結(jié)果顯示在圖7中。比較結(jié)果發(fā)現(xiàn)，盡管存在一些細(xì)微的差異，應(yīng)力的總體分布、應(yīng)力集中的位置以及應(yīng)力場值的趨勢(shì)都非常相似，而且應(yīng)力集中出現(xiàn)在流固交界面的上下兩側(cè)，不同于固體材料的左右兩側(cè)。

圖3 工況及網(wǎng)格劃分圖Fig. 3 Boundary conditions and mesh diagram

圖4 FCSHE(a)和MARC(b)的應(yīng)力在x方向上的應(yīng)力分量云圖Fig. 4 Contour plots of stress x of FCSHE (a) and MARC (b)

圖5 FCSHE(a)和MARC(b)的應(yīng)力在y方向上的應(yīng)力分量云圖Fig. 5 Contour plots of stress y of FCSHE (a) and MARC (b)

圖6 FCSHE(a)和MARC(b)的切應(yīng)力的應(yīng)力云圖Fig. 6 Contour plots of shear stress xy of FCSHE (a) and MARC (b)

圖7 應(yīng)力路徑圖Fig. 7 Diagram of stress on the path

在FCSHE方法中，流體與固體之間的界面上節(jié)點(diǎn)很少，并且每2個(gè)節(jié)點(diǎn)之間的位移采用線性插值法，因此界面上的位移場相對(duì)簡單。在MARC中，界面處節(jié)點(diǎn)很多，因此應(yīng)力場會(huì)更加豐富。這是兩個(gè)模型的計(jì)算結(jié)果不同的主要原因。

經(jīng)計(jì)算可知，F(xiàn)CSHE法得出的等效彈性模量為62.84 GPa，MARC中算出的等效彈性模量為62.60 GPa，兩者相對(duì)誤差為0.385%，證明了FCSHE法的有效性。

2.2 等效彈性模量的影響因素

2.2.1 體分比對(duì)等效彈性模量的影響

在8 mm×12 mm的長方形上建立模型，其上分布有10個(gè)含水孔。邊界條件為：上下兩側(cè)約束y方向位移，左側(cè)約束x方向位移，右側(cè)給定一個(gè)0. 3 μm的方向向右的位移。模型一共100個(gè)單元，外邊界含121個(gè)節(jié)點(diǎn)，流固交界面的邊界包含4個(gè)節(jié)點(diǎn)。含水孔體分比w分別為4%、8%、12%。在模擬分析中所用材料的彈性模量和泊松比與模型驗(yàn)證中所用材料相同。

圖8為不同體分比的單元模型示意圖，圖9為等效彈性模量隨體分比的變化情況。很容易看出，材料的等效彈性模量隨孔的體分比增大而減小，且基本呈線性變化。通過對(duì)體分比的研究得知，可以通過減小含流體孔洞材料的體分比提高材料的等效彈性模量。

圖8 不同體分比單元模型Fig. 8 Element model with different volume fractions

圖9 等效彈性模量隨體分比變化Fig. 9 The variation of equivalent elastic modulus with volume fraction

2.2.2 半徑對(duì)等效彈性模量的影響

8 mm×12 mm的長方形上分布有若干圓形含水孔。邊界條件為：上下兩側(cè)約束y方向位移，左側(cè)約束x方向位移，右側(cè)給定一個(gè)0.3 μm的方向向右的位移。單元模型一共64個(gè)單元，外邊界含81個(gè)節(jié)點(diǎn)，流固交界面上有4個(gè)節(jié)點(diǎn)。單元最大邊數(shù)為6，最大內(nèi)邊數(shù)為7。保持孔體分比都為7%，孔的半徑r為0.4,0.5,0.6 mm(單元模型如圖10所示)。

圖10 不同半徑單元模型Fig. 10 Element model with a different radius

在孔的體分比不變的情況下，半徑與孔的數(shù)量成反比。圖11顯示的是等效彈性模量隨半徑的變化情況，由圖可知，半徑增大時(shí)，等效彈性模量也逐漸增大，不呈線性。

圖11 等效彈性模量隨半徑變化Fig. 11 The variation of equivalent elastic modulus with radius

2.2.3 形狀對(duì)等效彈性模量的影響

8 mm×12 mm的長方形上分布有10個(gè)含水孔。邊界條件為：上下兩側(cè)約束y方向位移，左側(cè)約束x方向位移，右側(cè)給定一個(gè)0. 3 μm的方向向右的位移。單元模型一共64個(gè)單元，外邊界含81個(gè)節(jié)點(diǎn)，流固交界面上有4個(gè)節(jié)點(diǎn)。單元最大邊數(shù)為6，最大內(nèi)邊數(shù)為6。形狀改變方差控制了孔的長短軸，長短軸的取值符合正態(tài)分布式(38)。當(dāng)Svar為0時(shí)，孔的長短軸相等，孔呈圓形。當(dāng)Svar值增大時(shí)，保持其他控制參數(shù)不變，形狀改變方差Svar分別為0.00、0.10、0.15和0.20。形狀改變方差越大，孔越扁平，當(dāng)它為0時(shí)，孔為圓形。

(38)

式中：μ為平均數(shù)；σ為標(biāo)準(zhǔn)差，即Svar的算術(shù)平方根；a、f(a)為孔的長短軸。

圖12為不同方差單元模型圖，圖13為方差在0.00～0.20區(qū)間時(shí)等效彈性模量的變化情況。由圖13可知，在0.00～0.10區(qū)間時(shí)，等效彈性模量隨方差增大而增大。在0.10～0.20區(qū)間時(shí)，等效彈性模量隨方差增大而減小，且減小的幅度較大。如要得到符合正態(tài)分布的規(guī)律曲線，需要大量的樣本。

圖12 不同方差單元模型Fig. 12 Element model with different variances

圖13 等效彈性模量隨方差變化Fig. 13 The variation of equivalent elastic modulus with variance

2.2.4 空間分布對(duì)等效彈性模量的影響

在8 mm×12 mm的長方形上分布有3個(gè)圓形含水孔，孔半徑為0.5 mm。邊界條件為：上下兩側(cè)約束y方向位移，左側(cè)約束x方向位移，右側(cè)給定一個(gè)0.3 μm的方向向右的位移。單元模型一36個(gè)單元，外邊界含49個(gè)節(jié)點(diǎn)，流固交界面上有4個(gè)節(jié)點(diǎn)。單元最大邊數(shù)為6，最大內(nèi)邊數(shù)為3。保持其他控制參數(shù)不變，改變孔的空間分布位置，以3個(gè)孔的中心點(diǎn)連線與水平線的夾角為角度Sa，研究Sa分別為0°、30°、45°、60°時(shí)(如圖14所示)等效彈性模量的變化情況。由圖15可知，在0°到60°區(qū)間內(nèi)，等效彈性模量隨角度Sa增大而減小，并且?guī)缀醭示€性變化。

圖14 不同角度單元模型Fig. 14 Element model with different angles

圖15 等效彈性模量隨角度變化Fig. 15 The variation of equivalent elastic modulus with angle

3 結(jié) 論

基于普通應(yīng)力雜交元原理，推導(dǎo)了一種帶有新的流體相的固體單元(FCSHE)，在處理頁巖中大量孔的問題時(shí)，該單元可以大大減少計(jì)算量。它既包含流體又包含固體，適合作為流固耦合界面上的單元，可以很好地解決流固網(wǎng)格的過渡問題。建立了典型的四單元模型與普通商業(yè)有限元計(jì)算軟件MARC作對(duì)比，兩者之間應(yīng)力分布基本一致，等效彈性模量相對(duì)誤差在0.4%以內(nèi)，證明了其有效性。進(jìn)一步分別研究了形狀、體分比、空間分布位置和半徑變化時(shí)模型等效彈性模量的變化。結(jié)果顯示，不同形狀對(duì)等效彈性模量影響不大；同體分比時(shí)，隨含流體孔半徑增大，等效彈性模量隨之非線性增大；體分比增大時(shí)，等效彈性模量呈線性下降；當(dāng)角度逐漸增大時(shí)，等效彈性模量逐漸減小，并且下降幅度逐漸減小。