孫芳錦, 梁 爽， 張大明

(1.遼寧工程技術(shù)大學(xué) 建筑工程學(xué)院, 遼寧 阜新 123000; 2 遼寧工程技術(shù)大學(xué) 經(jīng)濟(jì)管理學(xué)院, 遼寧 阜新 123000)

物質(zhì)點(diǎn)法在風(fēng)與膜結(jié)構(gòu)耦合作用中的應(yīng)用研究

孫芳錦1, 梁 爽1， 張大明2

(1.遼寧工程技術(shù)大學(xué) 建筑工程學(xué)院, 遼寧 阜新 123000; 2 遼寧工程技術(shù)大學(xué) 經(jīng)濟(jì)管理學(xué)院, 遼寧 阜新 123000)

為克服常用數(shù)值方法在風(fēng)與膜結(jié)構(gòu)的耦合計(jì)算中網(wǎng)格易畸變、耗費(fèi)機(jī)時(shí)大的缺點(diǎn)，提出采用物質(zhì)點(diǎn)法計(jì)算風(fēng)與膜結(jié)構(gòu)的耦合作用。結(jié)構(gòu)域采用物質(zhì)點(diǎn)法描述，流體域采用浸入邊界法描述。結(jié)構(gòu)對(duì)流體的作用通過對(duì)結(jié)構(gòu)最邊緣的網(wǎng)格節(jié)點(diǎn)施加邊界條件來實(shí)現(xiàn)，流體對(duì)結(jié)構(gòu)的作用通過在結(jié)構(gòu)表面施加壓力和剪力實(shí)現(xiàn)。分別給出了結(jié)構(gòu)域、流體域和流固耦合計(jì)算的求解步驟。將該算法應(yīng)用于二維膜結(jié)構(gòu)的耦合計(jì)算，得到了結(jié)構(gòu)響應(yīng)和流體壓力場(chǎng)等結(jié)果，與采用有限元法計(jì)算結(jié)果進(jìn)行了對(duì)比，結(jié)果符合良好；同時(shí)與有限元法對(duì)比了計(jì)算效率，結(jié)果表明物質(zhì)點(diǎn)法計(jì)算效率要高于有限元法。結(jié)果表明物質(zhì)點(diǎn)法計(jì)算風(fēng)與膜結(jié)構(gòu)的流固耦合作用中可以得到準(zhǔn)確的結(jié)果，且計(jì)算效率得到較大提升。

膜結(jié)構(gòu)，流固耦合，物質(zhì)點(diǎn)法，計(jì)算效率

膜結(jié)構(gòu)柔性大、質(zhì)量輕等特點(diǎn)導(dǎo)致了風(fēng)荷載是其控制荷載。由于膜結(jié)構(gòu)較強(qiáng)的幾何非線性和材料非線性，風(fēng)與膜結(jié)構(gòu)的流固耦合作用是其抗風(fēng)設(shè)計(jì)中需要重點(diǎn)解決也是尚未很好解決的問題。數(shù)值模擬是研究流固耦合問題的有效方法之一，數(shù)值求解風(fēng)與膜結(jié)構(gòu)流固耦合問題的關(guān)鍵問題之一是避免流體和結(jié)構(gòu)交界面的網(wǎng)格過度變形，目前計(jì)算風(fēng)與膜結(jié)構(gòu)的流固耦合作用最常用的數(shù)值方法是有限元法、歐拉方法等, 這些方法在處理流固耦問題時(shí)常遇到網(wǎng)格畸變、界面追蹤和對(duì)流項(xiàng)處理等困難[1]，且網(wǎng)格劃分或重新生成也帶來了巨大的計(jì)算耗費(fèi)，也易導(dǎo)致計(jì)算精度下降或計(jì)算意外終止。因此研究高效簡(jiǎn)便的數(shù)值方法求解風(fēng)與膜結(jié)構(gòu)的耦合作用一直是計(jì)算風(fēng)工程領(lǐng)域被廣泛關(guān)注的熱點(diǎn)問題。

為克服上述傳統(tǒng)方法的缺點(diǎn)，本文采用物質(zhì)點(diǎn)法[2]求解風(fēng)與膜結(jié)構(gòu)的耦合作用。物質(zhì)點(diǎn)法是一種無網(wǎng)格方法，物質(zhì)點(diǎn)法的最大優(yōu)勢(shì)是非常適合處理像膜結(jié)構(gòu)這樣經(jīng)歷大變形的結(jié)構(gòu)，它可以避免有限元法存在的網(wǎng)格畸變以及糾正網(wǎng)格畸變重新生成網(wǎng)格的巨額計(jì)算資源等問題。物質(zhì)點(diǎn)法是一種將質(zhì)點(diǎn)和網(wǎng)格聯(lián)系起來的一種無網(wǎng)格方法，它采用互無關(guān)聯(lián)的拉格朗日物質(zhì)點(diǎn)離散結(jié)構(gòu)或流體，在這些物質(zhì)點(diǎn)上求解連續(xù)方程的所有變量，并采用規(guī)則的歐拉背景網(wǎng)格計(jì)算空間導(dǎo)數(shù)和動(dòng)量方程, 從而實(shí)現(xiàn)了質(zhì)點(diǎn)間的相互作用與聯(lián)系, 并避免了網(wǎng)格畸變問題。Sulsky 等[3]首先采用該方法分析了剛體的大轉(zhuǎn)動(dòng)和彈性體的碰撞問題, 展現(xiàn)了物質(zhì)點(diǎn)法潛在的優(yōu)勢(shì)。 此后, MPM 受到廣泛關(guān)注并得到了持續(xù)地發(fā)展。物質(zhì)點(diǎn)法的應(yīng)用廣泛，包括爆炸沖擊問題[4-5]，地質(zhì)學(xué)[6-7]，裂縫擴(kuò)展分析[8]，多相流分析[9]等。目前國(guó)內(nèi)外尚沒有將物質(zhì)點(diǎn)法應(yīng)用于風(fēng)與膜結(jié)構(gòu)流固耦合分析中的應(yīng)用。

本文采用物質(zhì)點(diǎn)法研究風(fēng)與膜結(jié)構(gòu)的流固耦合作用。結(jié)構(gòu)域采用物質(zhì)點(diǎn)法描述并給出求解步驟，考慮與流體域的數(shù)據(jù)交換，采用浸入邊界法對(duì)流體域進(jìn)行描述。浸入邊界表面采用非結(jié)構(gòu)三角形網(wǎng)格離散，這些表面網(wǎng)格節(jié)點(diǎn)構(gòu)成了一系列拉格朗日控制點(diǎn)，用以描述膜結(jié)構(gòu)的運(yùn)動(dòng)。給出了物質(zhì)點(diǎn)法求解流固耦合作用的步驟。最后將該方法應(yīng)用于二維膜結(jié)構(gòu)與風(fēng)的流固耦合分析中。

1 控制方程和邊界條件

1.1 流體控制方程

為便于與物質(zhì)點(diǎn)法描述的膜結(jié)構(gòu)的數(shù)據(jù)交換，這里采用曲線坐標(biāo)系下的浸入邊界法描述空氣流體，因曲線坐標(biāo)系下的浸入邊界法可以較好的處理邊界附近的流體域計(jì)算，也便于結(jié)構(gòu)域的邊界處理。

空氣流體的控制方程為N-S方程，在曲線坐標(biāo)下采用浸入邊界法可以表示為:

(1)

其中:

1.2 結(jié)構(gòu)控制方程

這里對(duì)經(jīng)歷大變形的膜結(jié)構(gòu)采用物質(zhì)點(diǎn)法進(jìn)行描述。浸入邊界表面采用非結(jié)構(gòu)三角形網(wǎng)格離散，這些表面網(wǎng)格節(jié)點(diǎn)構(gòu)成了一系列拉格朗日控制點(diǎn)，用以描述膜結(jié)構(gòu)的運(yùn)動(dòng)。膜結(jié)構(gòu)的控制方程可以表示為:

(2)

(3)

其中:ρs是膜結(jié)構(gòu)密度，σs是Cauchy應(yīng)力張量，bs是指定體積力，vs=dus/dt是結(jié)構(gòu)速度，us是結(jié)點(diǎn)位移，us=ds-Ds,ds是當(dāng)前(已變形)位形上的物質(zhì)點(diǎn)坐標(biāo)，Ds是初始(未變形)位形上的物質(zhì)點(diǎn)坐標(biāo)。

對(duì)于膜結(jié)構(gòu)，考慮膜材的材料非線性，其本構(gòu)關(guān)系可用超彈性模型描述，則σs可以通過如下方式求得：

(4)

(5)

其中:c1和c2是材料特性常數(shù)，k是體積模量，

J1=E1(E3)-1/3,
J2=E2(E3)-2/3,J3=(E3)1/2

(6)

(7)

1.3 邊界條件

在所有的浸入表面流體速度應(yīng)滿足

vf=vΩ, 在Ω上

(8)

其中:vΩ是浸入表面的速度。

結(jié)構(gòu)表面的表面牽引力應(yīng)滿足。

σf·nΩ=pΩ,在Ω上

(9)

其中:σf是流體應(yīng)力張量，pΩ是牽引力向量，nΩ是表面Ω的法線向量。

2 求解方法

2.1 結(jié)構(gòu)域的求解

結(jié)構(gòu)方程采用物質(zhì)點(diǎn)法進(jìn)行求解。這里的物質(zhì)點(diǎn)法采用無網(wǎng)格體系進(jìn)行，物質(zhì)點(diǎn)法認(rèn)為結(jié)構(gòu)域是由任意數(shù)量的物質(zhì)點(diǎn)構(gòu)成，物質(zhì)點(diǎn)可以分為內(nèi)部點(diǎn)和外部點(diǎn)(針對(duì)結(jié)構(gòu)表面而言)。

對(duì)于膜結(jié)構(gòu)，物質(zhì)點(diǎn)的質(zhì)量可以定義為:

(10)

其中:h0是膜結(jié)構(gòu)的厚度，ΔSi是第i個(gè)單元的表面積，Ni是離散膜結(jié)構(gòu)的單元數(shù)量，這里的變量在每一時(shí)間步進(jìn)行更新以滿足式(3)。在每一時(shí)間步里，將物質(zhì)點(diǎn)的信息向背景計(jì)算網(wǎng)格進(jìn)行內(nèi)插值計(jì)算求解結(jié)構(gòu)方程，該背景網(wǎng)格選定在Cartesian坐標(biāo)下。

在每一時(shí)間步采用物質(zhì)點(diǎn)法求解結(jié)構(gòu)域的過程可以歸納為：

(1) 對(duì)于每個(gè)物質(zhì)點(diǎn)對(duì)其進(jìn)行背景網(wǎng)格上的映射計(jì)算，

(11)

(2) 將動(dòng)量從物質(zhì)點(diǎn)向包含這些質(zhì)點(diǎn)的單元節(jié)點(diǎn)進(jìn)行映射，

(12)

(3) 求得在背景網(wǎng)格上每個(gè)點(diǎn)的內(nèi)力向量，對(duì)于膜結(jié)構(gòu)采用三角形單元離散，則有

(13)

其中:

(4) 根據(jù)動(dòng)態(tài)邊界條件(9)求得背景網(wǎng)格的外力向量，

(14)

并更新節(jié)點(diǎn)速度，

(15)

(5) 將當(dāng)前速度再映射回給物質(zhì)點(diǎn)，并計(jì)算物質(zhì)點(diǎn)當(dāng)前的位移和位置，

(16)

(17)

(18)

(6) 將物質(zhì)點(diǎn)的位移映射回給背景網(wǎng)格，并計(jì)算時(shí)間步Δt內(nèi)的變形梯度矩陣，

(19)

(20)

(7) 計(jì)算從t=0到t+tn+1的變形梯度矩陣、Green-Lagrangian應(yīng)力張量以及應(yīng)力，

(21)

(22)

(23)

可以看出，物質(zhì)點(diǎn)法將在每一時(shí)間步內(nèi)定義結(jié)構(gòu)新的速度和位置，求得的數(shù)值將傳遞給流體域進(jìn)行計(jì)算。

2.2 流體域的求解

流體域的求解主要采用文獻(xiàn)[10]的方法，該方法將壓力法和虛擬壓縮法相結(jié)合，采用四階Runga-Kutta方法計(jì)算求得對(duì)角線壓力算子，再進(jìn)行方程求解。受篇幅所限，詳細(xì)步驟可參考文獻(xiàn)[10]。

2.3 流固耦合算法

(4) 根據(jù)物質(zhì)點(diǎn)新的位置和壓力，求解結(jié)構(gòu)周圍新的流體方程，主要是估計(jì)浸入邊界節(jié)點(diǎn)的速度和壓力。這里采用的主要數(shù)值方法是將浸入邊界看成是零厚度單元進(jìn)行處理，根據(jù)結(jié)果表面已知解和流體節(jié)點(diǎn)，沿已經(jīng)定義好的結(jié)構(gòu)法線方向重新構(gòu)造解，將邊界條件應(yīng)用于浸入邊界附近的節(jié)點(diǎn)上。為方便任意復(fù)雜浸入邊界解的重構(gòu)，采用與結(jié)構(gòu)附近類似的三角形非結(jié)構(gòu)網(wǎng)格離散浸入邊界。

3 算例分析

這里采用上述方法對(duì)一膜結(jié)構(gòu)與風(fēng)的流固耦合作用進(jìn)行了研究。該膜結(jié)構(gòu)由受重力作用的彈性膜結(jié)構(gòu)屋蓋和兩個(gè)用以支撐膜屋蓋的豎向剛性墻構(gòu)成。膜結(jié)構(gòu)的計(jì)算簡(jiǎn)圖如圖1所示。膜結(jié)構(gòu)的厚度為0.01 m，密度為ρ=1.0×103kg/m3，彈性模量E=1.0×109N/m3，泊松系數(shù)ν=0，空氣流體的密度ρ=1.25 kg/m3，粘性系數(shù)μ=0.1 Ns/m2，入口風(fēng)速為13.75 m/s，特征長(zhǎng)度，即膜結(jié)構(gòu)的長(zhǎng)度l=10 m，因此雷諾數(shù)Re=1 700。這里未考慮湍流影響。

圖1 膜結(jié)構(gòu)計(jì)算簡(jiǎn)圖Fig.1 Computation sketch for membran structure

圖2 膜結(jié)構(gòu)表面的物質(zhì)點(diǎn)示意Fig.2 Material points sketch on the membrane structure

物質(zhì)點(diǎn)方法中，膜結(jié)構(gòu)表面由所有獨(dú)立的物質(zhì)點(diǎn)共同定義的。采用物質(zhì)點(diǎn)法模擬流固耦合問題的基本思路與其他類型的物質(zhì)點(diǎn)法模擬是一致的，區(qū)別是在流固耦合問題中，一些物質(zhì)點(diǎn)是膜材的，余下的物質(zhì)點(diǎn)是流體的。膜材的初試質(zhì)量是通過輸入膜結(jié)構(gòu)密度、膜材厚度和用戶指定點(diǎn)間距確定的。流體對(duì)結(jié)構(gòu)的作用或結(jié)構(gòu)對(duì)流體的作用是通過求解每個(gè)網(wǎng)格節(jié)點(diǎn)的動(dòng)量方程確定的。膜結(jié)構(gòu)表面的物質(zhì)點(diǎn)示意如圖2所示。流體域的離散使用了19 486/44 818個(gè)三角形單元/節(jié)點(diǎn)，時(shí)間步Δt=1.0×10-4s，膜結(jié)構(gòu)時(shí)間步Δt=5×10-4s，經(jīng)過試算，這種網(wǎng)格尺寸是獨(dú)立的，即對(duì)最后的計(jì)算結(jié)果無影響。計(jì)算中共使用了5次子迭代。

圖3分別給出了膜結(jié)構(gòu)中點(diǎn)豎向位移時(shí)程和膜結(jié)構(gòu)后方流體豎向速度時(shí)程，從圖中可以看出，經(jīng)過初試階段后膜結(jié)構(gòu)的響應(yīng)趨于規(guī)律性的變化，這是由于經(jīng)歷初試階段后，漩渦開始以接近膜結(jié)構(gòu)基本頻率從膜結(jié)構(gòu)前緣脫落，這也導(dǎo)致膜結(jié)構(gòu)響應(yīng)和空氣漩渦間存在較強(qiáng)的流固耦合作用。

圖3 膜結(jié)構(gòu)中點(diǎn)豎向位移和流體豎向速度時(shí)程Fig.3 Time histories of structure vertical displacement and fluid vertical velocity

表1 膜結(jié)構(gòu)響應(yīng)統(tǒng)計(jì)值比較

為了說明本文方法的計(jì)算效率，將本文的物質(zhì)點(diǎn)法的計(jì)算效率與采用有限元法的計(jì)算效率進(jìn)行了對(duì)比，如表2所示。

表1中給出了采用本文方法與采用有限元法計(jì)算得到的結(jié)構(gòu)響應(yīng)參數(shù)的比較，表1 中給出了采用本文方法與采用有限元法計(jì)算得到的結(jié)構(gòu)響應(yīng)參數(shù)的比較，其中有限元法的流體網(wǎng)格為26 957個(gè)，計(jì)算殘差小于等于10-4時(shí)認(rèn)為計(jì)算收斂。從表中可以看出，采用本文的物質(zhì)點(diǎn)法計(jì)算得到的膜結(jié)構(gòu)響應(yīng)與有限元法計(jì)算結(jié)果非常接近，證明了本文方法的正確性。

表2 不同方法計(jì)算效率對(duì)比

從表2中可以看出，采用本文的物質(zhì)點(diǎn)法達(dá)到收斂于指定值所需迭代次數(shù)要遠(yuǎn)小于有限元法，計(jì)算效率提高約42%；在經(jīng)歷同樣的迭代次數(shù)后，物質(zhì)點(diǎn)法的計(jì)算殘差要小于有限元法計(jì)算殘差，說明采用本文物質(zhì)點(diǎn)法的計(jì)算效率要高于一般的有限元法，這主要是由于物質(zhì)點(diǎn)法對(duì)于流固交界面處的變形處理方法較簡(jiǎn)單，因?yàn)橹恍杼幚砟そY(jié)構(gòu)表面的物質(zhì)點(diǎn)，而無需像有限元方法那樣處理大量的網(wǎng)格變形問題，這樣就節(jié)省了較多的計(jì)算資源，簡(jiǎn)化了問題。

圖4和圖5分別給出了瞬時(shí)壓力場(chǎng)分布和速度場(chǎng)分布。從圖中可以看出，漩渦沿著膜結(jié)構(gòu)的前緣脫落，并形成了典型的渦街。同時(shí)可以看出當(dāng)漩渦經(jīng)過尾緣時(shí)并沒有在結(jié)構(gòu)壁面的下游形成回流。值得注意的是，本例中的漩渦脫離現(xiàn)象主要是由膜結(jié)構(gòu)的振動(dòng)引起的，因?yàn)槟そY(jié)構(gòu)屋蓋表面的壓力變化過小，如果沒有結(jié)構(gòu)的振動(dòng)是不會(huì)形成漩渦的。

圖4 不同時(shí)刻瞬時(shí)壓力圖Fig.4 Instantaneous pressure at different time

圖5 膜結(jié)構(gòu)周圍速度場(chǎng) (t=55)Fig.5 Velocity field around membrane structure

4 結(jié) 論

本文采用物質(zhì)點(diǎn)法求解風(fēng)與膜結(jié)構(gòu)的耦合作用，采用物質(zhì)點(diǎn)法描述膜結(jié)構(gòu)并求解，為適應(yīng)與物質(zhì)點(diǎn)法在交界面處的數(shù)據(jù)交換，流體域采用浸入邊界法進(jìn)行描述，分別給出了物質(zhì)點(diǎn)法求解結(jié)構(gòu)方程和流固耦合作用的求解步驟，并將該方法應(yīng)用于二維膜結(jié)構(gòu)的耦合分析中，得到了較理想的結(jié)果，包括結(jié)構(gòu)響應(yīng)和流場(chǎng)壓力、速度等分布；計(jì)算對(duì)比發(fā)現(xiàn)物質(zhì)點(diǎn)法的計(jì)算效率要高于有限元法。計(jì)算證實(shí)了物質(zhì)點(diǎn)法在風(fēng)與膜結(jié)構(gòu)耦合分析中的準(zhǔn)確性和高效性，物質(zhì)點(diǎn)法最大的優(yōu)點(diǎn)在于避免了有限元法中對(duì)于交界面處網(wǎng)格的繁瑣計(jì)算，大大提高了準(zhǔn)確性和計(jì)算效率。本文目前只是將其應(yīng)用于二維算例中，后續(xù)研究應(yīng)進(jìn)一步將其擴(kuò)展到三維膜結(jié)構(gòu)的耦合計(jì)算中。另外計(jì)算中發(fā)現(xiàn)在計(jì)算膜結(jié)構(gòu)內(nèi)力和速度時(shí)有時(shí)會(huì)出現(xiàn)計(jì)算不穩(wěn)定性，因此如何提高計(jì)算過程中的穩(wěn)定性也是值得進(jìn)一步研究的問題。

[1] 劉欣.無網(wǎng)格方法[M].北京: 科學(xué)出版社, 2011.

[2] Sulsky D, Schreyer H L. Axisymmetric form of the material point method with applications to upsetting and taylor impact problems[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139 (1-4): 409-429.

[3] Sulsky D, Chen Z, Schreyer H L. A particle method for history-dependent materials[J]. Computer Methods in Applied Mechanics and Engineering, 1994,118(1-2): 179-196.

[4] Hu W, Chen Z, Model-based simulation of the synergistic effects of blast and fragmentation on a concrete wall using the MPM[J]. Int. J. Impact Eng,2006,32 (12):2066-2096.

[5] Lian Y P, Zhang X, Zhou X, et al., Numerical simulation of explosively driven metal by material point method[J]. Int. J. Impact Eng, 2011, 38 :237-245.

[6] Zhang H W, Wang K P, Chen Z, Material point method for dynamic analysis of saturated porous media under external contact/impact of solid bodies[J]. Comput. Meth. Appl. Mech. Eng, 2009, 198: 1456-1472.

[7] Andersen S, Andersen L. Modelling of landslides with the material-point method[J].Comput. Geosci,2010, 14: 137-147.

[8] Nairn J A. Material point method calculations with explicit cracks[J]. Comput. Model Eng. Sci.,2003, 4: 649-663.

[9] Zhang D Z, Zou Q, VanderHeyden W B, et al, Material point method applied to multiphase flows[J] J. Comput. Phys,2008, 227: 3159-3173.

[10] Sotiropoulos F, Constantinescu G. Pressure-based residual smoothing operators for multistage pseudo-compressibility algorithms[J]. Journal of Computational Physics,1997,133(1):129-145.

Application of material point method in wind-membrane structure interaction

SUN Fang-jin1, LIANG Shuang1, ZHANG Da-ming2

(1. College of Architecture and civil engineering, Liaoning Technical University, Fuxin 123000, China;2. College of Economy & Management, Liaoning Technical University, Fuxin 123000, China)

Due to the disadvantages of mesh distortion and large computer time consuming in commonly used numerical methods, material point method (MPM) was proposed to deal with the wind-membrane structure interaction. The solid was described by MPM, whereas the fluid was described using immersed boundary method. The effect of the solid on fluid was taken into account through imposing boundary conditions on the immediate vicinity mesh nodes. The effect of the fluid on solid was taken into account by imposing pressure and shear stress on the solid surface. And the procedures for solving the solid equations, fluid equations and the fluid-structure interaction were presented respectively. The algorithm was applied to a 2D membrane structure interacting with wind. The results of structural responses and fluid pressure field were obtained, which are in good agreement with those obtained using finite element method (FEM). The computing efficiency was compared between MPM and FEM, where MPM shows better efficiency. It is shown that ideal results can be obtained in wind-structure interaction analysis by employing MPM, with great improvement in computing efficiency.

membrane structure; fluid-structure interaction; material point method; computing efficiency

國(guó)家自然科學(xué)基金(51108345);同濟(jì)大學(xué)土木工程防災(zāi)國(guó)家重點(diǎn)實(shí)驗(yàn)室開放基金(SLDRCE-MB-04);遼寧省教育廳基金一般項(xiàng)目(L2013134);江蘇省工程力學(xué)分析重點(diǎn)實(shí)驗(yàn)室(東南大學(xué))開放課題基金項(xiàng)目資助

2013-11-01 修改稿收到日期：2014-04-21

孫芳錦 女，博士，副教授，1981年生

TU973

A

10.13465/j.cnki.jvs.2015.09.018