多孔結構多尺度隨機振動分析的漸近均勻化-時域顯式法*

蘇 成，羅俊哲，許 秩

（1.華南理工大學 土木與交通學院，廣州 510640；2.華南理工大學 亞熱帶建筑科學國家重點實驗室，廣州 510640；3.人工智能與數字經濟廣東省實驗室，廣州 510330）

引言

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

由于多孔結構具有非均質特性，因此采用單一尺度有限元法進行結構分析時，需要采用非常精細的有限元網格，導致計算量十分巨大.多尺度分析方法[7-8]通過建立不同尺度的模型，把單一尺度下的復雜問題分解為不同尺度下的簡單問題，可在保證計算精度的前提下大幅提高整體計算效率，是解決多孔結構計算問題的一類重要方法.在多尺度分析方法中，漸近均勻化法[9-11]是一種具有嚴格數學基礎的方法，當引入的特征參數ε趨近于0 時，該方法的計算結果收斂于精確解答.

在隨機動力荷載作用下多孔結構隨機振動問題的研究尚不多見.文獻[12]采用功率譜法研究了理想白噪聲平穩隨機荷載作用下多孔彈性板的隨機振動問題.然而，對于非平穩隨機振動問題，為了獲得結構響應的演化功率譜，功率譜法需要在大量的離散頻率點上進行時域積分運算，難以應用于實際工程中的多孔結構問題.近年來提出的一類非平穩隨機振動時域顯式法[13-16]，通過構建結構動力響應的時域顯式表達式，能夠在時域內直接建立非平穩響應統計矩的顯式列式，實現任意自由度上的降維計算，已成功應用于非均質結構非平穩隨機響應問題[17].

本文綜合多尺度漸近均勻化法和隨機振動時域顯式法的計算優勢，實現了非平穩隨機激勵下多孔結構隨機振動問題的高效求解.首先采用多尺度漸近均勻化法推導并求解了多孔結構動力問題的多尺度控制微分方程，建立了多孔結構宏觀和細觀動力響應的時域顯式表達式；然后結合結構隨機振動時域顯式法，獲得了非平穩隨機激勵下多孔結構動力響應統計矩的演化規律.數值算例表明，本文所提出的方法具有理想的計算精度與計算效率.

1 多尺度動力分析的漸近均勻化法

1.1 多尺度動力學模型

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

圖1 周期性多孔結構和宏觀等效結構Fig.1 The periodic porous structure and the macroscopic equivalent structure

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

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

相應地，Φ(x,y)的偏導數可以用Φε(x,y,ξ,η)的偏導數表示為

在單一尺度坐標下，多孔結構的控制微分方程可表示為

式中

按式(2)和式(3)的方式擴充坐標后，多尺度坐標下的多孔結構控制微分方程可表示為

式中

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

對多尺度坐標下的位移向量進行關于特征參數ε的多尺度漸近展開，并保留前三項得[18]

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

將式(8)代入式(6)，整理后得到多尺度坐標下多孔結構控制微分方程關于特征參數ε的展開形式為

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

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

1.2 單胞控制微分方程

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

式中X=X(ξ,η)，稱為單胞特征位移矩陣，其只與單胞內細觀坐標(ξ,η)有關，表示為

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

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

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

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

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

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

根據單胞特征位移矩陣和本構方程進一步定義單胞特征應力矩陣為

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

1.3 宏觀控制微分方程

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

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

式中V為單胞的面積.

文獻[21]證明了以下等式：

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

式中

顯然，式(23)可以理解為關于宏觀位移向量u(0)的宏觀運動微分方程，式中 ρH，cH，DH和fH分別為宏觀質量密度矩陣、宏觀阻尼系數矩陣、宏觀彈性矩陣和宏觀體力向量.該方程的求解區域為圖1所示的宏觀等效結構區域，其邊界條件與實際問題一致.式(24)～(27)稱為均勻化方程，用于求解各均勻化宏觀參數，其中宏觀彈 性矩陣DH依賴于1.2 小節獲得的單胞特征應力矩陣 Υ.

1.4 回代分析

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

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

單胞細觀應力向量 σε可根據幾何方程和本構方程由uε求得，表達式為

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

式中單胞特征應力矩陣 Υ見式(19).

2 多尺度隨機振動分析的時域顯式法

2.1 宏觀動力響應時域顯式表達式

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

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

式中M，C和K分別為宏觀質量、阻尼和剛度矩陣；分別為宏觀節點位移、速度和加速度向量；F為非平穩隨機激勵向量，L為其定位矩陣.

記時間步長為 ?t，積分步數為n，把隨機激勵向量F離散為時刻0,t1,t2,···,tn處的荷載向量F0,F1,F2,···,Fn.假定結構響應具有零初值，采用Newmark-β數值積分格式求解式(31)，可以推導得到宏觀狀態向量的顯式表達式為

式中

其中

γ=0.5,β=0.25.

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

表1 各時刻顯式表達式的系數矩陣Table 1 Coefficient matrices for explicit formulation at different instants

2.2 細觀動力響應時域顯式表達式

為了進一步建立單胞細觀動力響應時域顯式表達式，需要采用式(28)和式(30)進行回代分析.由該兩式可以得到單胞細觀單元節點位移向量及應力向量與宏觀節點狀態向量之間的關系：

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

式中Ai,j的表達式見式(33).

在多孔結構的隨機振動分析中，通常只需關注結構的某些關鍵響應.假設r為所關注的關鍵響應，如位移、應力等，則由式(36)可直接得到ri=r(ti)的顯式表達式為

2.3 動力響應統計矩

前述工作利用多尺度漸近均勻化法建立了多孔結構動力響應的顯式表達式，揭示了多孔結構的物理演化規律.在此基礎上，可以進一步研究動力響應統計矩的演化規律.根據一階矩和二階矩的運算規則，由式(38)可得各時刻關鍵響應ri的均值和方差如下：

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

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

3 數值算例

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

該周期性多孔有機玻璃方形板在上部受到零均值非平穩隨機激勵f(t)的作用，f(t)取為均勻調制非平穩隨機過程，即

式中g(t)為調制函數，用于反映隨機過程的非平穩特性，取為

q(t)為零均值平穩隨機過程，其相關函數取為

式中τ為時間間隔，λ=144 MPa2.于是，f(t)的相關函數可以表達為

根據式(45)所示的相關函數，可以采用Cholesky 分解法生成隨機激勵f(t)的樣本，其中一個樣本如圖2所示.

圖2 隨機激勵f(t)的一個樣本Fig.2 A sample of random excitation f(t)

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

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

圖3 點A 水平位移時程Fig.3 Time histories of the horizontal displacement at point A

圖4 點A 豎向位移時程Fig.4 Time histories of the vertical displacement at point A

圖5 點A 豎向正應力時程Fig.5 Time histories of the vertical normal stress at point A

圖6 點A 水平位移方差時程Fig.6 Time histories of variance of the horizontal displacement at point A

圖7 點A 豎向位移方差時程Fig.7 Time histories of variance of the vertical displacement at point A

圖8 點A 豎向正應力方差時程Fig.8 Time histories of variance of the vertical normal stress at point A

在計算效率方面，多尺度時域顯式法和單一尺度時域顯式法的計算時間分別列于表2 中.計算在 Intel core i5 處理器和 8 GB 內存的臺式計算機上完成.由表2 可見，兩種方法的計算時間由兩部分組成，分別用于構建動力響應時域顯式表達式和統計矩運算.兩種方法在統計矩運算方面的計算時間是一樣的，主要區別在于動力響應顯式表達式系數矩陣的計算時間.由于采用了多尺度漸近均勻化法，多尺度時域顯式法用于構建動力響應顯式表達式的計算時間僅為單一尺度方法的0.94%，從而顯著提高了多孔結構非平穩隨機振動問題的計算效率.

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

4 結論

為了正確反映多孔材料細觀非均質特性及動力荷載隨機性的影響，本文開展了多孔結構隨機振動分析的漸近均勻化-時域顯式法研究，得到了以下結論：

1)所推導的周期性多孔結構動力問題多尺度控制微分方程充分體現了多尺度漸近均勻化思想，具有嚴格的數學基礎，為多孔結構多尺度動力響應分析奠定了理論基礎.

2)采用多尺度有限元法建立的多孔結構多尺度動力響應顯式表達式，完全揭示了多孔結構的物理演化過程.在此基礎上，通過多孔結構物理演化機制和概率演化機制的相對分離，實現了多孔結構細觀動力響應演變統計矩的高效計算.

3)所提出的漸近均勻化-時域顯式法綜合發揮了多尺度漸近均勻化法和隨機振動時域顯式法的計算優勢，大幅提升了多孔結構非平穩隨機振動分析的計算精度和計算效率，為多孔結構隨機振動問題提供了一條可行的途徑.

值得指出的是，本文的方法并不局限于周期性多孔結構的分析，亦可用于一般周期性非均質材料的隨機振動問題.在本文工作基礎上，可以進一步研究非均質材料微結構隨機振動靈敏度及隨機拓撲優化問題.