多孔結(jié)構(gòu)多尺度隨機(jī)振動(dòng)分析的漸近均勻化-時(shí)域顯式法*

蘇 成，羅俊哲，許 秩

（1.華南理工大學(xué) 土木與交通學(xué)院，廣州 510640；2.華南理工大學(xué) 亞熱帶建筑科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，廣州 510640；3.人工智能與數(shù)字經(jīng)濟(jì)廣東省實(shí)驗(yàn)室，廣州 510330）

引言

多孔結(jié)構(gòu)是由含一定數(shù)量相互貫通或封閉孔洞的固體材料組成的結(jié)構(gòu)[1].由于具有高比強(qiáng)、高比剛度等優(yōu)點(diǎn)，多孔結(jié)構(gòu)在土木工程、機(jī)械工程和航天航空工程等領(lǐng)域得到了廣泛應(yīng)用[2-6].由于周期性多孔結(jié)構(gòu)更易制備，因而在工程中的應(yīng)用更為廣泛.在工程應(yīng)用中，多孔結(jié)構(gòu)通常承受隨機(jī)動(dòng)力荷載的作用，其隨機(jī)振動(dòng)及動(dòng)力可靠度分析對(duì)保障結(jié)構(gòu)安全具有重要的意義.

由于多孔結(jié)構(gòu)具有非均質(zhì)特性，因此采用單一尺度有限元法進(jìn)行結(jié)構(gòu)分析時(shí)，需要采用非常精細(xì)的有限元網(wǎng)格，導(dǎo)致計(jì)算量十分巨大.多尺度分析方法[7-8]通過(guò)建立不同尺度的模型，把單一尺度下的復(fù)雜問(wèn)題分解為不同尺度下的簡(jiǎn)單問(wèn)題，可在保證計(jì)算精度的前提下大幅提高整體計(jì)算效率，是解決多孔結(jié)構(gòu)計(jì)算問(wèn)題的一類(lèi)重要方法.在多尺度分析方法中，漸近均勻化法[9-11]是一種具有嚴(yán)格數(shù)學(xué)基礎(chǔ)的方法，當(dāng)引入的特征參數(shù)ε趨近于0 時(shí)，該方法的計(jì)算結(jié)果收斂于精確解答.

在隨機(jī)動(dòng)力荷載作用下多孔結(jié)構(gòu)隨機(jī)振動(dòng)問(wèn)題的研究尚不多見(jiàn).文獻(xiàn)[12]采用功率譜法研究了理想白噪聲平穩(wěn)隨機(jī)荷載作用下多孔彈性板的隨機(jī)振動(dòng)問(wèn)題.然而，對(duì)于非平穩(wěn)隨機(jī)振動(dòng)問(wèn)題，為了獲得結(jié)構(gòu)響應(yīng)的演化功率譜，功率譜法需要在大量的離散頻率點(diǎn)上進(jìn)行時(shí)域積分運(yùn)算，難以應(yīng)用于實(shí)際工程中的多孔結(jié)構(gòu)問(wèn)題.近年來(lái)提出的一類(lèi)非平穩(wěn)隨機(jī)振動(dòng)時(shí)域顯式法[13-16]，通過(guò)構(gòu)建結(jié)構(gòu)動(dòng)力響應(yīng)的時(shí)域顯式表達(dá)式，能夠在時(shí)域內(nèi)直接建立非平穩(wěn)響應(yīng)統(tǒng)計(jì)矩的顯式列式，實(shí)現(xiàn)任意自由度上的降維計(jì)算，已成功應(yīng)用于非均質(zhì)結(jié)構(gòu)非平穩(wěn)隨機(jī)響應(yīng)問(wèn)題[17].

本文綜合多尺度漸近均勻化法和隨機(jī)振動(dòng)時(shí)域顯式法的計(jì)算優(yōu)勢(shì)，實(shí)現(xiàn)了非平穩(wěn)隨機(jī)激勵(lì)下多孔結(jié)構(gòu)隨機(jī)振動(dòng)問(wèn)題的高效求解.首先采用多尺度漸近均勻化法推導(dǎo)并求解了多孔結(jié)構(gòu)動(dòng)力問(wèn)題的多尺度控制微分方程，建立了多孔結(jié)構(gòu)宏觀和細(xì)觀動(dòng)力響應(yīng)的時(shí)域顯式表達(dá)式；然后結(jié)合結(jié)構(gòu)隨機(jī)振動(dòng)時(shí)域顯式法，獲得了非平穩(wěn)隨機(jī)激勵(lì)下多孔結(jié)構(gòu)動(dòng)力響應(yīng)統(tǒng)計(jì)矩的演化規(guī)律.數(shù)值算例表明，本文所提出的方法具有理想的計(jì)算精度與計(jì)算效率.

1 多尺度動(dòng)力分析的漸近均勻化法

1.1 多尺度動(dòng)力學(xué)模型

周期性多孔結(jié)構(gòu)由周期性分布的帶孔單胞組成，如圖1所示.單胞的尺度比整個(gè)結(jié)構(gòu)的尺度小得多，即w×h?W×H.設(shè)兩者的比為ε量級(jí)，0<ε ?1，ε稱為特征參數(shù).在多尺度漸近均勻化分析方法中，需要用到宏觀坐標(biāo)(x,y)和細(xì)觀坐標(biāo)(ξ,η)，這兩組坐標(biāo)滿足以下關(guān)系式：

圖1 周期性多孔結(jié)構(gòu)和宏觀等效結(jié)構(gòu)Fig.1 The periodic porous structure and the macroscopic equivalent structure

假設(shè)場(chǎng)函數(shù)Φ(x,y)表示單一尺度坐標(biāo)下多孔結(jié)構(gòu)的材料場(chǎng)（如彈性模量、Poisson 比等）、荷載場(chǎng)（如體力）和結(jié)構(gòu)響應(yīng)場(chǎng)（如位移、應(yīng)力和應(yīng)變等）.對(duì)于周期性多孔結(jié)構(gòu)，場(chǎng)函數(shù)Φ(x,y)原則上應(yīng)包含周期性部分.將周期性部分改用細(xì)觀坐標(biāo)(ξ,η)表達(dá)，從而得到擴(kuò)充坐標(biāo)后的多尺度坐標(biāo)下多孔結(jié)構(gòu)的場(chǎng)函數(shù)Φε(x,y,ξ,η)，它是關(guān)于細(xì)觀坐標(biāo)(ξ,η)的周期函數(shù)，周期為單胞尺寸.顯然，Φ(x,y)和 Φε(x,y,ξ,η)滿足以下關(guān)系：

式中k1和k2為整數(shù)，w和h為單胞尺寸.

相應(yīng)地，Φ(x,y)的偏導(dǎo)數(shù)可以用Φε(x,y,ξ,η)的偏導(dǎo)數(shù)表示為

在單一尺度坐標(biāo)下，多孔結(jié)構(gòu)的控制微分方程可表示為

式中

按式(2)和式(3)的方式擴(kuò)充坐標(biāo)后，多尺度坐標(biāo)下的多孔結(jié)構(gòu)控制微分方程可表示為

式中

其中帶ε的函數(shù)分別為式(5)中各函數(shù)擴(kuò)充坐標(biāo)后得到的多尺度坐標(biāo)下的函數(shù)；?ε為關(guān)于(ξ,η)的微分算子矩陣.值得注意的是，由于周期性多孔結(jié)構(gòu)的材料屬性呈周期性變化，因此ρε(ξ,η)，cε(ξ,η)，Eε(ξ,η)和νε(ξ,η)與宏觀坐標(biāo)(x,y)無(wú)關(guān).

對(duì)多尺度坐標(biāo)下的位移向量進(jìn)行關(guān)于特征參數(shù)ε的多尺度漸近展開(kāi)，并保留前三項(xiàng)得[18]

式中u(0)，u(1)和u(2)分別為由位移的零階、一階和二階展開(kāi)函數(shù)組成的向量.

將式(8)代入式(6)，整理后得到多尺度坐標(biāo)下多孔結(jié)構(gòu)控制微分方程關(guān)于特征參數(shù)ε的展開(kāi)形式為

保留前三項(xiàng)并考慮到ε的任意性，上式成立的必要條件是

在式(8)中通常采用u(0)來(lái)反映位移隨宏觀坐標(biāo)(x,y)的變化，u(0)稱為宏觀位移向量.顯然，u(0)與細(xì)觀坐標(biāo)(ξ,η)無(wú)關(guān)，即u(0)=u(0)(x,y)，因此式(10)自動(dòng)成立.由式(11)和式(12)可以分別導(dǎo)出單胞控制微分方程和宏觀控制微分方程.

1.2 單胞控制微分方程

對(duì)u(1)進(jìn)行分離變量，將它表示為以下乘積形式：

式中X=X(ξ,η)，稱為單胞特征位移矩陣，其只與單胞內(nèi)細(xì)觀坐標(biāo)(ξ,η)有關(guān)，表示為

將式(13)代入式(11)，可得

由于?u(0)不恒等于0，所以上式成立的必要條件為

式(16)只與細(xì)觀坐標(biāo)(ξ,η)相關(guān)，用于求解單胞特征位移矩陣，稱為單胞控制微分方程.該方程也可以表達(dá)為如下形式：

式中Xi表示X的第i列；ei表示為

顯然，式(17)可以理解為在不同方向單位初應(yīng)變ei作用下的單胞靜力平衡微分方程，解之即可獲得Xi(i=1,2,3)，從而按式(14)構(gòu)造單胞特征位移矩陣X.值得注意的是，對(duì)于周期性多孔結(jié)構(gòu)，在求解式(17)時(shí)，只需要考慮一個(gè)典型單胞區(qū)域.對(duì)于圖1所示的具有雙對(duì)稱軸的單胞區(qū)域，考慮到周期性邊界條件的要求，在單位初應(yīng)變e1或e2作用下，需對(duì)單胞外圍邊界施以法向固定的邊界條件；在單位初應(yīng)變e3作用下，需對(duì)單胞外圍邊界施以切向固定的邊界條件[19].

顯然，式(17)可以理解為在不同方向單位初應(yīng)變ei作用下的單胞靜力平衡微分方程，解之即可獲得Xi(i=1,2,3)，從而按式(14)構(gòu)造單胞特征位移矩陣X.值得注意的是，對(duì)于周期性多孔結(jié)構(gòu)，在求解式(17)時(shí)，只需要考慮一個(gè)典型單胞區(qū)域.對(duì)于圖1所示的具有雙對(duì)稱軸的單胞區(qū)域，考慮到周期性邊界條件的要求，在單位初應(yīng)變e1或e2作用下，需對(duì)單胞外圍邊界施以法向固定的邊界條件；在單位初應(yīng)變e3作用下，需對(duì)單胞外圍邊界施以切向固定的邊界條件[19].

根據(jù)單胞特征位移矩陣和本構(gòu)方程進(jìn)一步定義單胞特征應(yīng)力矩陣為

該矩陣可用于計(jì)算下文的宏觀彈性矩陣.

1.3 宏觀控制微分方程

將式(13)代入式(12)，并考慮式(19)得

根據(jù)Fredholm 定理[20]，關(guān)于u(2)的橢圓型偏微分方程有且僅有唯一解的必要條件為方程的非齊次項(xiàng)在單胞內(nèi)積分為0，即

式中V為單胞的面積.

文獻(xiàn)[21]證明了以下等式：

將式(22)代入式(21)中，可得宏觀控制微分方程：

式中

顯然，式(23)可以理解為關(guān)于宏觀位移向量u(0)的宏觀運(yùn)動(dòng)微分方程，式中 ρH，cH，DH和fH分別為宏觀質(zhì)量密度矩陣、宏觀阻尼系數(shù)矩陣、宏觀彈性矩陣和宏觀體力向量.該方程的求解區(qū)域?yàn)閳D1所示的宏觀等效結(jié)構(gòu)區(qū)域，其邊界條件與實(shí)際問(wèn)題一致.式(24)～(27)稱為均勻化方程，用于求解各均勻化宏觀參數(shù)，其中宏觀彈 性矩陣DH依賴于1.2 小節(jié)獲得的單胞特征應(yīng)力矩陣 Υ.

1.4 回代分析

將式(13)代入式(8)，并略去 ε2項(xiàng)，可得單胞細(xì)觀位移向量uε的表達(dá)式為

式中單胞特征位移矩陣X可由單胞靜力平衡微分方程式(17)求得，而宏觀位移向量u(0)則由宏觀運(yùn)動(dòng)微分方程式(23)求得.

單胞細(xì)觀應(yīng)力向量 σε可根據(jù)幾何方程和本構(gòu)方程由uε求得，表達(dá)式為

將式(28)代入式(29)并略去 ε1項(xiàng)，可得

式中單胞特征應(yīng)力矩陣 Υ見(jiàn)式(19).

2 多尺度隨機(jī)振動(dòng)分析的時(shí)域顯式法

2.1 宏觀動(dòng)力響應(yīng)時(shí)域顯式表達(dá)式

首先進(jìn)行單胞分析.采用靜力有限元法求解式(17)所示的單胞靜力平衡微分方程，此時(shí)需要將圖1所示的單胞劃分為細(xì)觀有限元網(wǎng)格，由此求得單胞內(nèi)各細(xì)觀單元節(jié)點(diǎn)的特征位移矩陣和特征應(yīng)力矩陣，進(jìn)一步由式(26)獲得式(23)所需的宏觀彈性矩陣.

然后進(jìn)行宏觀分析.采用動(dòng)力有限元法求解式(23)所示的宏觀運(yùn)動(dòng)微分方程，此時(shí)需要將圖1所示的宏觀等效結(jié)構(gòu)劃分為宏觀有限元網(wǎng)格（宏觀單元的尺寸一般取為單胞尺寸），從而式(23)可以離散為以下形式：

式中M，C和K分別為宏觀質(zhì)量、阻尼和剛度矩陣；分別為宏觀節(jié)點(diǎn)位移、速度和加速度向量；F為非平穩(wěn)隨機(jī)激勵(lì)向量，L為其定位矩陣.

記時(shí)間步長(zhǎng)為 ?t，積分步數(shù)為n，把隨機(jī)激勵(lì)向量F離散為時(shí)刻0,t1,t2,···,tn處的荷載向量F0,F1,F2,···,Fn.假定結(jié)構(gòu)響應(yīng)具有零初值，采用Newmark-β數(shù)值積分格式求解式(31)，可以推導(dǎo)得到宏觀狀態(tài)向量的顯式表達(dá)式為

式中

其中

γ=0.5,β=0.25.

根據(jù)式(33)所揭示的系數(shù)矩陣之間的內(nèi)在關(guān)系，可以把各系數(shù)矩陣排列如表1所示.由表1 可見(jiàn)，僅第一列系數(shù)矩陣Ai,0和第二列系數(shù)矩陣Ai,1(i=1,2,···,n)需要計(jì)算和存儲(chǔ)，其計(jì)算量相當(dāng)于對(duì)宏觀等效結(jié)構(gòu)進(jìn)行2m次脈沖響應(yīng)分析的計(jì)算量，m為荷載向量F的分量個(gè)數(shù)[22].

表1 各時(shí)刻顯式表達(dá)式的系數(shù)矩陣Table 1 Coefficient matrices for explicit formulation at different instants

2.2 細(xì)觀動(dòng)力響應(yīng)時(shí)域顯式表達(dá)式

為了進(jìn)一步建立單胞細(xì)觀動(dòng)力響應(yīng)時(shí)域顯式表達(dá)式，需要采用式(28)和式(30)進(jìn)行回代分析.由該兩式可以得到單胞細(xì)觀單元節(jié)點(diǎn)位移向量及應(yīng)力向量與宏觀節(jié)點(diǎn)狀態(tài)向量之間的關(guān)系：

將式(32)代入式(35)得

式中Ai,j的表達(dá)式見(jiàn)式(33).

在多孔結(jié)構(gòu)的隨機(jī)振動(dòng)分析中，通常只需關(guān)注結(jié)構(gòu)的某些關(guān)鍵響應(yīng).假設(shè)r為所關(guān)注的關(guān)鍵響應(yīng)，如位移、應(yīng)力等，則由式(36)可直接得到ri=r(ti)的顯式表達(dá)式為

2.3 動(dòng)力響應(yīng)統(tǒng)計(jì)矩

前述工作利用多尺度漸近均勻化法建立了多孔結(jié)構(gòu)動(dòng)力響應(yīng)的顯式表達(dá)式，揭示了多孔結(jié)構(gòu)的物理演化規(guī)律.在此基礎(chǔ)上，可以進(jìn)一步研究動(dòng)力響應(yīng)統(tǒng)計(jì)矩的演化規(guī)律.根據(jù)一階矩和二階矩的運(yùn)算規(guī)則，由式(38)可得各時(shí)刻關(guān)鍵響應(yīng)ri的均值和方差如下：

式中E(F[i])為F[i]的均值向量，cov(F[i],F[i])為F[i]的協(xié)方差矩陣，分別表示為

式 中μF(t)和RFF(t,τ)分別為非平穩(wěn)隨機(jī)激勵(lì)向量F(t)的均值函數(shù)向量和互相關(guān)函數(shù)矩陣，i=1,2,···,n.

3 數(shù)值算例

如圖1所示，底部固支的周期性多孔有機(jī)玻璃方形板邊長(zhǎng)和厚度分別為W=H=100 mm 和t=1 mm，單胞邊長(zhǎng)為w=h=10 mm，單胞中心有一邊長(zhǎng)為a=b=5 mm 的方孔.有機(jī)玻璃方形板[23]的彈性模量E=5.3 GPa，Poisson 比μ=0.3，密度ρ = 1.18 × 10?6kg/mm3.選取Rayleigh 阻尼模型，阻尼比取為ξ= 0.002.

該周期性多孔有機(jī)玻璃方形板在上部受到零均值非平穩(wěn)隨機(jī)激勵(lì)f(t)的作用，f(t)取為均勻調(diào)制非平穩(wěn)隨機(jī)過(guò)程，即

式中g(shù)(t)為調(diào)制函數(shù)，用于反映隨機(jī)過(guò)程的非平穩(wěn)特性，取為

q(t)為零均值平穩(wěn)隨機(jī)過(guò)程，其相關(guān)函數(shù)取為

式中τ為時(shí)間間隔，λ=144 MPa2.于是，f(t)的相關(guān)函數(shù)可以表達(dá)為

根據(jù)式(45)所示的相關(guān)函數(shù)，可以采用Cholesky 分解法生成隨機(jī)激勵(lì)f(t)的樣本，其中一個(gè)樣本如圖2所示.

圖2 隨機(jī)激勵(lì)f(t)的一個(gè)樣本Fig.2 A sample of random excitation f(t)

分別采用基于漸近均勻化法的多尺度有限元法以及單一尺度有限元法建立多孔結(jié)構(gòu)動(dòng)力響應(yīng)時(shí)域顯式表達(dá)式.在多尺度有限元法中，首先將圖1所示的單胞劃分為300 個(gè)四節(jié)點(diǎn)單元，通過(guò)單胞分析獲得多孔材料宏觀彈性模量和宏觀Poisson 比分別為EH= 2.864 GPa 和μH= 0.1967；然后將整個(gè)宏觀等效方形板劃分為100 個(gè)四節(jié)點(diǎn)單元，并進(jìn)行宏觀分析.在單一尺度有限元法中，將整個(gè)多孔方形板劃分為30000 個(gè)四節(jié)點(diǎn)單元，每個(gè)單元的尺寸與多尺度有限元法中單胞分析的細(xì)觀單元尺寸一致.此外，在建立動(dòng)力響應(yīng)時(shí)域顯式表達(dá)式時(shí)，所考慮的持時(shí)為T(mén)= 1 ms，時(shí)間步長(zhǎng)取為 ?t= 0.002 ms.

分別采用多尺度時(shí)域顯式法和單一尺度時(shí)域顯式法，考察圖1所示單胞右側(cè)孔邊中點(diǎn)A的位移和應(yīng)力.在圖2所示荷載樣本作用下，點(diǎn)A的位移時(shí)程uxA,uyA和應(yīng)力時(shí)程 σyA分別如圖3～5所示.此外，在式(42)所示的非平穩(wěn)隨機(jī)激勵(lì)作用下，點(diǎn)A的位移方差時(shí)程D(uxA),D(uyA)和應(yīng)力方差時(shí)程D(σyA)分別如圖6～8所示.由圖3～8 可見(jiàn)，兩種方法的計(jì)算結(jié)果相當(dāng)吻合，驗(yàn)證了多尺度時(shí)域顯式法的正確性.

圖3 點(diǎn)A 水平位移時(shí)程Fig.3 Time histories of the horizontal displacement at point A

圖4 點(diǎn)A 豎向位移時(shí)程Fig.4 Time histories of the vertical displacement at point A

圖5 點(diǎn)A 豎向正應(yīng)力時(shí)程Fig.5 Time histories of the vertical normal stress at point A

圖6 點(diǎn)A 水平位移方差時(shí)程Fig.6 Time histories of variance of the horizontal displacement at point A

圖7 點(diǎn)A 豎向位移方差時(shí)程Fig.7 Time histories of variance of the vertical displacement at point A

圖8 點(diǎn)A 豎向正應(yīng)力方差時(shí)程Fig.8 Time histories of variance of the vertical normal stress at point A

在計(jì)算效率方面，多尺度時(shí)域顯式法和單一尺度時(shí)域顯式法的計(jì)算時(shí)間分別列于表2 中.計(jì)算在 Intel core i5 處理器和 8 GB 內(nèi)存的臺(tái)式計(jì)算機(jī)上完成.由表2 可見(jiàn)，兩種方法的計(jì)算時(shí)間由兩部分組成，分別用于構(gòu)建動(dòng)力響應(yīng)時(shí)域顯式表達(dá)式和統(tǒng)計(jì)矩運(yùn)算.兩種方法在統(tǒng)計(jì)矩運(yùn)算方面的計(jì)算時(shí)間是一樣的，主要區(qū)別在于動(dòng)力響應(yīng)顯式表達(dá)式系數(shù)矩陣的計(jì)算時(shí)間.由于采用了多尺度漸近均勻化法，多尺度時(shí)域顯式法用于構(gòu)建動(dòng)力響應(yīng)顯式表達(dá)式的計(jì)算時(shí)間僅為單一尺度方法的0.94%，從而顯著提高了多孔結(jié)構(gòu)非平穩(wěn)隨機(jī)振動(dòng)問(wèn)題的計(jì)算效率.

表2 兩種方法的計(jì)算時(shí)間（單位：s）Table 2 Time costs of the 2 methods (unit:s)

4 結(jié)論

為了正確反映多孔材料細(xì)觀非均質(zhì)特性及動(dòng)力荷載隨機(jī)性的影響，本文開(kāi)展了多孔結(jié)構(gòu)隨機(jī)振動(dòng)分析的漸近均勻化-時(shí)域顯式法研究，得到了以下結(jié)論：

1)所推導(dǎo)的周期性多孔結(jié)構(gòu)動(dòng)力問(wèn)題多尺度控制微分方程充分體現(xiàn)了多尺度漸近均勻化思想，具有嚴(yán)格的數(shù)學(xué)基礎(chǔ)，為多孔結(jié)構(gòu)多尺度動(dòng)力響應(yīng)分析奠定了理論基礎(chǔ).

2)采用多尺度有限元法建立的多孔結(jié)構(gòu)多尺度動(dòng)力響應(yīng)顯式表達(dá)式，完全揭示了多孔結(jié)構(gòu)的物理演化過(guò)程.在此基礎(chǔ)上，通過(guò)多孔結(jié)構(gòu)物理演化機(jī)制和概率演化機(jī)制的相對(duì)分離，實(shí)現(xiàn)了多孔結(jié)構(gòu)細(xì)觀動(dòng)力響應(yīng)演變統(tǒng)計(jì)矩的高效計(jì)算.

3)所提出的漸近均勻化-時(shí)域顯式法綜合發(fā)揮了多尺度漸近均勻化法和隨機(jī)振動(dòng)時(shí)域顯式法的計(jì)算優(yōu)勢(shì)，大幅提升了多孔結(jié)構(gòu)非平穩(wěn)隨機(jī)振動(dòng)分析的計(jì)算精度和計(jì)算效率，為多孔結(jié)構(gòu)隨機(jī)振動(dòng)問(wèn)題提供了一條可行的途徑.

值得指出的是，本文的方法并不局限于周期性多孔結(jié)構(gòu)的分析，亦可用于一般周期性非均質(zhì)材料的隨機(jī)振動(dòng)問(wèn)題.在本文工作基礎(chǔ)上，可以進(jìn)一步研究非均質(zhì)材料微結(jié)構(gòu)隨機(jī)振動(dòng)靈敏度及隨機(jī)拓?fù)鋬?yōu)化問(wèn)題.