基于拓撲優化的結構彈性成像方法1)

付君健 *, 李帥虎 李 好 ** 高 亮 ** 周祥曼 *, 田啟華 *,

* (水電機械設備設計與維護湖北省重點實驗室,湖北宜昌 443002)

? (三峽大學機械與動力學院,湖北宜昌 443002)

** (華中科技大學數字制造裝備與技術國家重點實驗室,武漢 430074)

引言

彈性成像源自醫學影像領域,是一種通過彈性模量表征生物體組織物理特性的成像方法[1-2].彈性成像能夠獲得常規成像技術難以獲取的彈性模量信息,輔助臨床診斷生物體器官病理變化的位置和形狀.對于工程中的機械裝備結構,結構損傷會造成特定部位彈性模量發生改變,繼而影響結構的整體物理特性.將彈性模量與損傷表征關聯,發展結構缺陷識別的彈性成像技術,對于機械工程裝備的服役安全評價具有重要意義.

與超聲成像[3]、電阻抗成像[4-5]、熱成像[6-7]、太赫茲成像[8-9]等先進的成像方法類似,彈性成像是電子技術、計算機技術和反演數學相結合的產物,它利用計算機對穿越物體的物理信號進行處理,通過算法重建物體內部信息.彈性成像的基本原理是對生物體組織施加一個外部靜態或動態激勵.根據胡克定律,生物體組織在外力作用下將產生位移、應變等物理響應,結合超聲波測距法和數字信號處理技術,可以預測生物體內部組織的彈性模量信息.盡管超聲成像、X 射線等無損檢測技術可實現結構的缺陷和損傷識別,但多數檢測技術屬于局部識別技術,對環境要求高,需要先驗信息,難以檢測整體結構性能缺陷.對于大型機械裝備結構,損傷一般發生在結構內部或難以觀測位置,發展全局性的損傷識別技術具有重要的工程意義.作為一種全局性的損傷識別方法,彈性成像是成像參數、外界激勵和反演模型三大要素的綜合表達.彈性成像的參數是材料的彈性模量,一般情況下彈性階段的力學響應即可實現成像,如結構靜態位移[10]和動力響應[11],不必拓展至塑性分析階段.通過對比仿真響應和實測響應,可以建立彈性成像的反演模型.彈性成像屬于約束優化問題,一般作為反問題求解,其關鍵技術是彈性成像的反演算法.

拓撲優化是一種在給定設計空間和約束下,尋求滿足性能最優的材料分布形式的結構優化方法[12-14].拓撲優化與彈性成像均屬于約束優化問題,其本質都是為了尋找一種滿足目標函數的彈性模量分布.因此,拓撲優化也可用于彈性成像,逆向反演結構的內部損傷所對應的彈性模量.在進行結構損傷識別時,由于拓撲優化幾何描述能力強,不需要預設損傷形狀,可自動反演復雜的損傷形狀.例如對于帶有孔洞缺陷的結構,采用水平集法[15-17]或移動變形組件法[18]及類似方法[19],可有效描述缺陷幾何形狀.雖然在靜力響應和動力響應下,拓撲優化方法都表現出了較好的損傷識別效果[10,20],但仍存在兩種問題需要解決.其一,水平集法具有清晰幾何邊界,一般用于帶有孔洞缺陷的結構損傷識別.當結構存在多種不同彈性模量的材料時,或結構存在缺損、夾雜、裂紋等多種損傷形式,水平集法則難以對結構缺陷進行有效識別,需要在拓撲優化框架中同時引入幾何形狀和彈性模量作為反演變量[21-22].其二,目前采用水平集法或移動變形組件法進行結構損傷識別時,雖然無需預設缺陷形狀,但需預估缺陷數量,初始孔洞或初始組件數量對反演結果有重要影響[23].因此,研究適用于多損傷并存、初始缺陷不敏感的彈性成像方法,提高彈性成像的局部表征和全局識別能力,具有更廣泛的理論和工程意義.

由于彈性成像的數學本質是優化驅動的反問題,與結構拓撲優化方法數學原理相同,本文提出了一種基于拓撲優化的彈性成像方法.通過構建損傷表征、結構模型與物理響應的映射關系,建立并求解結構彈性成像的優化模型,并探討了在不同結構維度、損傷類型、缺陷數量下結構彈性成像的效果.

1 彈性屬性與缺陷表征

拓撲優化本質是尋求滿足設計需求的最優材料分布,基于人工材料插值模型將單元密度與彈性模量耦合,調控單元彈性模量對結構性能的貢獻[24-26],這種對彈性模量的全局調控機制與結構彈性成像技術不謀而合.彈性模量可以衡量物體抵抗彈性變形的能力,結構缺陷會造成特定部位材料彈性模量發生改變,繼而影響結構的整體物理特性.若將彈性模量與損傷表征關聯,則有望通過彈性模量的空間分布識別結構多種損傷,反映損傷對結構物理特性的綜合影響.

為了解決結構變形、裂紋、缺損等多損傷的識別問題,不同損傷類型的統一表征是關鍵.受拓撲優化理論啟發,采用結構離散單元的相對密度(彈性模量系數)作為彈性成像參數來表征損傷程度,建立結構缺陷、彈性屬性和力學模型之間的耦合關系,成像參數與損傷表征的數學定義見式(1).如圖1 所示,塊狀分布的中間彈性模量值可表征夾雜,帶狀分布的中/低彈性模量值可表征裂紋,塊狀分布的低彈性模量值可表征孔洞.通過缺陷的形態和成像參數的數值,即可實現對結構缺陷的評估.

圖1 損傷表征Fig.1 Damage characterization

基于拓撲優化理論的固體材料各向同性懲罰模型(solid isotropic material with penalization,SIMP)[27],構建成像參數x與彈性模量的非線性插值模型如下

式中,x為成像參數,取值范圍為 0 .001 ≤x≤1 ,x最小值取0.001 為避免剛度矩陣奇異.p為懲罰系數,其作用是為了加速優化模型收斂.E0為實體材料的彈性模量,E為被懲罰后的單元彈性模量.

進一步構建彈性屬性和力學模型之間的關系,將融合了結構缺陷信息的彈性模量引入結構有限元平衡方程,以結構靜力學問題為例,有限元平衡方程如下

式中,K為關于成像參數x的結構全局剛度矩陣,U為結構的全局位移響應,F為激勵載荷向量,xe為單元e的成像參數,ke為單元的剛度矩陣,B為應變位移矩陣,D為實體材料彈性矩陣,Ω 代表結構域,Ωe代表單元e的結構域.通過求解式(3),可得到受到缺陷影響的結構位移響應,繼而構建了損傷表征、結構模型與物理響應的映射關系.

2 彈性成像優化模型

為體現缺陷對結構物理響應的影響,以損傷結構響應和無損結構響應的最小二乘為目標函數.再定義與結構自由度維度相同的測量位置向量 ψa,通過在最小二乘函數前乘以測量位置向量,可實現結構特定位置響應誤差的最小化

式中,J(x) 為關于成像參數的目標函數,ψa為結構響應的測量位置向量(在自由度a處元素為1,其余元素為0),a代表結構響應測量位置的自由度編號,上標 T 代表矩陣轉置,U0為無損結構響應,U為含缺陷信息的損傷結構響應.

考慮到結構自由度上的響應值的數量級可能不同,以結構響應變化率的平方為目標函數,定義如下

在研究靜力學結構彈性成像理論時,僅考慮結構位移響應,以成像模型和真實模型位移響應變化率的平方為目標函數,以靜力學平衡方程為約束條件,以成像參數上下限為約束條件,構建式(8)所示的結構彈性成像模型

式中,x為元素取值在[0.001,1]之間的成像參數向量,N代表成像參數的數量,J(x) 表示目標函數,U0為初始結構(無損結構)位移響應,U為含缺陷信息的結構(損傷結構)位移響應,K為結構全局剛度矩陣,F為全局載荷向量.

從優化模型(8)可知,傳統拓撲優化方法屬于正向設計,即給定設計目標計算新型結構形式[28-31];而彈性成像屬于逆向設計,通過已知結構響應反演結構的拓撲構型(彈性模量分布),從而實現結構缺陷形態和位置、損傷程度的識別.

3 優化模型求解

3.1 敏度分析

為了求解彈性成像優化模型,本文基于伴隨法推導目標函數關于成像參數的導數[32],構造拉格朗日函數如下

式中,λ 為靜力學平衡方程的拉格朗日乘子列向量.

求解拉格朗日乘子列向量得

則拉格朗日函數關于成像參數的導數為

3.2 數值實施

考慮成像參數的有界性約束,成像參數的迭代更新策略如下

為了增加收斂曲線的數值穩定性,將第n步敏度和n+1步敏度加權平均,以修正靈敏度信息

經過敏度修正后的成像參數更新策略如下

結構彈性成像優化模型的收斂準則定義如下

式中,ε 代表一個極小正數,nmax代表最大迭代數.

4 數值算例

4.1 二維均質材料結構彈性成像

為驗證彈性成像的有效性,首先對二維均質材料結構進行彈性成像研究.如圖2 所示的二維懸臂梁結構,結構域長寬為 2 m×1 m,懸臂梁左端固定,右端中心處施加有F=-1 N 的集中載荷,將結構域劃分為 4 0×20 個雙線性四邊形單元,材料彈性模量E=100 GPa ,泊松比v=0.3 (如無特殊說明,本文僅考慮各向同性材料) .本算例的位移響應測量點設置在懸臂梁上、下邊緣中點處.

圖2 二維均質懸臂梁結構Fig.2 Two-dimensional homogeneous cantilever beam structure

對于均質各向同性材料結構,因為只存在唯一成像參數,需對敏度信息進行均勻化處理

為了得到響應點的位移響應,圖3 所示,本文僅通過有限元法進行結構變形的數值模擬,得到兩個測量點X和Y方向的位移響應作為模型輸入.優化模型的初始彈性模量設置為 2 00 GPa,成像參數初始值為1,迭代步長為 ξ =0.1 .由圖4 可知,經過280 步迭代后,目標函數收斂至6.594 × 10-5,成像參數收斂至0.5.將成像參數與初始彈性模量相乘,得到此時的結構彈性模量約為 1 00 GPa,反演得到的結構彈性模量與設置的彈性模量一致.對于均質材料結構,結構響應與成像參數之間呈線型關系,即使以少量測量點的位移響應為參照,也能準確反演整體結構的彈性模量,此方法可在結構材料未知的情況下對彈性模量進行成像.

圖3 二維均質懸臂梁彈性成像結構Fig.3 Two-dimensional homogeneous cantilever beam elastography structure

圖4 二維均質懸臂梁彈性成像迭代圖Fig.4 Iterative graph of two-dimensional homogeneous cantilever beam elastography

4.2 二維非均質材料結構彈性成像

為驗證本方法對非均質材料結構成像的有效性,如圖5 所示,采用算例4.1 中的結構參數和邊界條件,材料彈性模量設置為E=200 GPa,泊松比均v=0.3 ,將結構域劃分為4 0×20 個雙線性四邊形單元.彈性成像的最小分辨率為1 個單元尺寸,在懸臂梁長度方向1/4 和3/4 的位置及寬度方向1/2 位置設置尺寸為 0 .2 m×0.2 m 的孔洞缺陷,孔洞缺陷區域材料的彈性模量設置為 0 .001E,迭代步長為 ξ =0.1,分別設置單孔洞缺陷和雙孔洞缺陷形式,進行非均質結構的彈性成像.本算例的位移響應測量點設置在懸臂梁上、下邊緣處.

圖5 內置缺陷的二維懸臂梁結構Fig.5 Two-dimensional cantilever beam structure with built-in damaged

圖6(a)所示為存在A缺陷結構的成像結果,優化模型迭代259 步后目標函數收斂,結構內部出現較低的彈性模量區域,結構彈性成像結果較清晰,最終彈性成像呈現的缺陷位置和大小與設定缺陷基本一致.圖6(b)所示為存在B缺陷結構的成像結果,優化模型迭代93 步后目標函數收斂,結構內部缺陷位置的彈性模量與設定缺陷有一定誤差.雖然B缺陷彈性成像結果較模糊,但已初步表明缺陷存在的大致位置.圖7 所示為同時考慮A和B兩處缺陷的結構成像過程.在優化模型迭代至100 步時,左邊A缺陷的彈性模量先識別出來,右邊B缺陷的彈性模量識別稍慢.隨著迭代的進行,A缺陷的彈性模量成像逐漸清晰,可正確識別A缺陷的位置和大小,B缺陷彈性模量成像的清晰度一般,但仍可以表征此處存在一定缺陷,成像結果與A、B缺陷單獨成像的結果基本一致.在前述研究中[10,20-21,23],若結構存在多個缺陷,在算法識別時需要提前設定缺陷數量.而本文在單元上定義了成像參數,成像參數的數量與單元數量相同,每一處的單元都可用以表征缺陷.因此,本文方法無需提前人為指定缺陷數量,算法可自動識別缺陷處的彈性模量,通過彈性模量的變化判斷是否存在缺陷.

圖6 二維懸臂梁結構單缺陷彈性成像結果Fig.6 Single damage elastography results of the two-dimensional cantilever beam structure

圖7 二維懸臂梁結構多缺陷彈性成像過程Fig.7 Multi-damage elastography process of the two-dimensional cantilever beam structure

由圖6～ 圖8 可知,不同位置缺陷的彈性成像質量有一定差異,主要存在兩方面的原因.其一,成像質量與缺陷對結構性能的影響程度有關.以靜力學為例,不同位置的材料對測量點處結構變形的貢獻度不同,若缺陷對測量點處的結構變形影響較小(如缺陷B),則彈性成像效果會偏弱,若缺陷對結構測量點處的變形影響較大,則成像質量效果會偏強(如缺陷A).在本例中若將約束置于右邊界,載荷置于左邊界中點處,則B缺陷的識別效果會更好.其二,彈性成像質量與測量點數量有關.不同于均質材料結構的彈性成像,非均勻材料結構的結構響應與成像參數之間呈現高度非線性,采用較少的測量點無法準確對含有缺陷的結構進行彈性成像.本算例僅采用懸臂梁上、下邊界來測量結構響應,沒有用到結構內部的測量點.在實際工程中,測量點的數量有限、位置特殊,有限的測量點會一定程度影響彈性成像的結果.

圖8 二維懸臂梁結構彈性成像迭代圖Fig.8 Iterative diagram of elastography of two-dimensional cantilever beam structure

4.3 邊界條件對成像結果的影響

為驗證彈性成像方法的通用性,研究邊界條件對彈性成像的影響.定義圖9 所示的二維米歇爾(Michell)梁結構,結構域長寬為 2 m×1 m,米歇爾梁左下角固定約束,右下角支撐約束.在下端中心處施加F=-1 N 的集中載荷,將結構域劃分為 4 0×20 個雙線性四邊形單元.材料彈性模量E=200 GPa,泊松比v=0.3 .在米歇爾梁長度方向1/4 和3/4 的位置及寬度方向1/2 位置設置尺寸為 0 .2 m×0.2 m 的孔洞缺陷,孔洞缺陷區域材料的彈性模量設置為 0 .001E,分別設置單孔洞缺陷和雙孔洞缺陷形式,進行不同邊界條件的非均質結構的彈性成像.

圖9 內置缺陷的二維Michell 梁結構Fig.9 Two-dimensional Michell-like beam structure with built-in damaged

對比圖6 和圖10 可知,同樣的二維長方形結構,在米歇爾梁的約束條件下,依然可以得到缺陷A和缺陷B的彈性模量圖,說明本方法在不同的邊界條件下本方法仍具有良好的成像特性.在不同的邊界條件下,彈性成像結果會存在微小的差異,因為同樣的缺陷在不同邊界條件下,對結構性能的影響不同,進而會影響彈性成像優化模型的反演.

圖10 二維Michell 梁結構彈性成像結果Fig.10 Two-dimensional Michell-like beam structure elastography results

4.4 三維結構彈性成像

為進一步驗證本方法對非均質材料結構成像的通用性,對三維非均質材料結構進行彈性成像研究.如圖11 所示的三維懸臂梁結構,結構域長、寬、高為L×H×W=2 m×1 m×1 m,懸臂梁左端固定,右端面中心處沿Y軸方向施加F=-1 N 的線載荷,將結構域劃分為 4 0×20×20 個六面體單元,材料彈性模量E=200 GPa ,泊松比v=0.3 .在懸臂梁中心設置尺

圖11 內置缺陷的三維懸臂梁結構Fig.11 Three-dimensional cantilever beam structure with built-in damaged

寸為 0 .2 m×0.2 m×0.2 m 的缺陷,考慮孔洞和夾雜兩種缺陷形式,孔洞缺陷彈性模量設置為 0 .001E,夾雜缺陷彈性模量設置為 0 .5E.在懸臂梁上、下邊緣處設置位移響應測量點,進行三維非均質結構的彈性成像.

以懸臂梁上、下端面節點的X和Y自由度位移響應為參照.圖12(a)所示為夾雜缺陷的三維彈性成像結果,缺陷處的成像參數接近0.5.圖12(b)所示為孔洞缺陷的三維彈性成像結果,缺陷處的成像參數接近0.由于缺陷位于結構內部,圖12 僅顯示半邊結構空間的成像結果以方便觀察.圖13 所示優化迭代曲線光滑,表明算法具有良好的穩定性.由圖12 可知,本方法可對不同損傷程度的缺陷形式進行彈性成像,不僅能識別缺陷的位置和相對大小,還能識別缺陷的近似彈性模量值.由二維彈性成像擴展至三維彈性成像時,僅需要定義新的三維結構有限元模型.三維彈性成像計算成本高于二維問題,彈性成像的變量定義、優化模型和求解算法均與二維彈性成像一致.

圖12 三維懸臂梁結構彈性成像結果Fig.12 Three-dimensional cantilever beam structure elasticity elastography

圖13 三維懸臂梁結構彈性成像迭代圖Fig.13 Three-dimensional cantilever beam structure elastography iteration diagram

5 結論

本文提出一種基于拓撲優化的結構彈性成像方法,構建了損傷表征、結構模型與物理響應的映射關系,建立并求解了結構彈性成像的優化模型.研究結果表明,彈性成像方法在二維、三維問題上具有通用性,可有效實現均質/非均質、單/多缺陷結構的彈性成像,且彈性成像結果不依賴于特定邊界條件.本文雖然僅提出了靜態響應下的彈性成像理論及求解方法,但該方法可進一步拓展至復雜動態響應的結構彈性成像.