基于絕對節點坐標法的平面梁有限變形下變形重構1)

吳懋琦 譚述君 , 高飛雄

* (大連理工大學工業裝備結構分析國家重點實驗室,遼寧大連 116024)

? (大連理工大學遼寧省空天飛行器前沿技術重點實驗室,遼寧大連 116024)

引言

隨著智能材料和傳感器技術的發展,對各類柔性結構的健康監測和主動控制得到了越來越多的應用.在航天領域,對于大型柔性空間結構,尤其是大型天線及反射器,其結構變形的掌握程度直接關乎其性能和精度[1],在結構完整性監測和故障診斷上也很重要;在醫療器械領域,先進柔性穿刺針和其他微創手術工具對于精確地識別其在人體內的載荷狀況和變形具有很高的要求[2];在機器人領域,可靠高效的變形感知技術是相對傳統剛性機器人具有結構簡單、成本低廉、強適應力等許多優勢的連續介質機器人實用化必須的技術[3].對于以上問題,利用分布式光纖網絡[4]進行結構健康監測[5]和主動控制[6]方面的信息獲取是一個非常具有前景的解決方案.其中,如何利用離散的應變信息進行變形重構是相關問題的核心,但已有的利用應變信息對柔體結構進行姿態監測的研究較少,且往往是利用單純基于曲率的幾何關系而非將力學中更完備的應變-位移關系也納入考慮來進行曲面或曲線的變形重構,這限制了相關方法在獲取更完整精確的變形信息方面的能力,也不利于后續將變形重構得出的信息直接應用于載荷識別和健康監測等方面.

早期,Davis 等[7]進行了利用光纖光柵傳感器進行形態感知的初步嘗試,利用梁在特定載荷下的變形模式對測量的應變進行擬合,進而實現對變形的預測;Foss 和Haugse[8]利用位移振型和應變振型的對應關系,通過模態轉換給出了位移-應變轉換矩陣,將擬合用的變形模式范圍由特定載荷情況擴展到了動力學的模態范疇;Tessler 等[9-10]和Gherlone等[11]基于有限元方法中應變—位移的變分關系,用最小二乘列式的形式進行求解,發展出了一種逆有限元方法,進一步拓寬了相關問題的求解空間,相較之前的方法在準確性和理論完備性上得到了顯著的提升.Tessler 及其團隊在后續對iFEM 方法在單元類型[12-13]和應用方面[14-15]進行了大量的延伸,逆有限元方法因其不受材料性質及載荷特征影響的特性和其與有限元方法高度相關帶來的可擴展性引起了許多其他研究者關注[16-17].

以上提到的方法均可以視作變形重構問題在線性的試探函數法框架下求解方法的演進,也意味著它們均只適用于應變與位移映射關系可線性近似的小變形情況,對于在航空航天、精密醫療器械等領域得到越來越廣泛運用的柔性結構無法取得足夠好的效果,而后者相較于前者對于變形感知也具有更高的需求.為了應對這些需求,Ko 等[18]基于經典梁理論,對梁進行簡單的分段和線性近似,并采用類似于浮動坐標系的方法連接各段,從而在低成本線性計算的基礎上完成非線性的近似,對于作為其研究對象的機翼的較大變形取得了不錯的效果,但這種較為粗糙的近似無法滿足較高的精度要求,且實際未考慮軸向變形的影響,而對于那些更大變形的情況,此方法則變得不再適用;Yi 等[19]將應變向曲率轉化,再將曲率在被測梁長度方向上進行連續化處理來獲取梁上各點的空間坐標,這對于大變形情況具有極佳的適應性,但在曲率的連續性方面處理有一定困難,文中采取的數值積分方法在計算效率方面也有局限性.基于曲率的方法在變形重構問題及相關應用的研究中實際上得到了大量的應用,如Roesthuis 等[20]將FBG 傳感器陣列采集的應變信息利用基于曲率的重構方法轉化為形態信息,以實現柔性針頭刺入軟組織后的精確控制;Xu 等[21]使用正交網狀的FBG 傳感器陣列結構建立了一種柔性板三維形狀測量系統,同樣也是通過應變信息到曲率信息再到形態信息的方式實現的.

需要注意的是,以上的非線性變形重構都是單純基于幾何關系的變形重構,且都是以曲率-轉角的微分關系為基礎建立相關算法.此類方法并非基于系統完整的應變-位移場力學關系,實際上忽視了被測體長度方向上變化的影響,也不能方便地與力學相關方面的理論與應用相結合,這都限制了上述方法的適用范圍.

絕對節點坐標法(absolute nodal coordinate formulation,ANCF)是Shabana 等[22-24]為了處理強非線性耦合問題而發展出的一種方法,通過取各單元節點在整體坐標系下完整的位置和角度作為廣義坐標實現了較為精確的非線性建模.絕對節點坐標法建立的質量陣為常數陣,廣義力表達簡潔,無科氏力和離心力項,在大變形下的動力學問題中得到了廣泛地應用.對于本文涉及的有限變形下平面梁變形重構問題而言,絕對節點坐標法提供的應變-位移關系模型簡潔,能夠精確地將應變與節點位移間的關系以標準的二次型形式表達,進而幫助建立起可靠易行的求解算法.

本文部分繼承了變形重構逆有限元方法的思想,將梁的應變-位移重構問題視為一種最優化問題,旨在對被測對象進行有限元離散后,通過引入ANCF 形函數對大變形下的非線性應變-節點位移關系進行描述,構建目標函數并以位移作為決策變量.繼而采用最優化方法進行位移的求解,建立一種能夠具有可擴展性和較高精確性的有限變形下應變-位移重構方法.相對于已有的曲率連續化方法而言,利用有限元方法,將應變-位移場進行較為系統的重構,不僅在軸向大變形條件下具有較好的適應性,對在此方面具有需求的相關工作,例如文獻[25]中涉及的振動控制和文獻[26]中涉及的載荷識別也具有相當的意義.

1 有限變形下的平面梁變形重構問題

如圖1 所示,一個典型的基于應變的平面梁變形重構模型由對稱布置在被測梁上下表面各點的應變傳感器組成.當平面梁受載荷影響發生變形時,應變傳感器獲取變形導致的應變,利用相關應變信息理論上可以準確地重構出被測梁的變形狀態.

圖1 基于應變信息的平面梁變形重構模型Fig.1 Shape reconstruction model of plane beam based on strain information

各種變形重構方法中對分布式的應變測量數據進行處理的方式各不相同,如上文提到的基于曲率連續性的方法中將應變轉化為曲率,或是模態疊加法中將應變轉化為模態坐標.但這些方法無疑都是利用連續化的應變-位移關系描述來處理離散的應變測量值與位移之間的關系.而基于此,考慮到用一種普適性的方式對平面梁重構問題進行描述有助于更完整地展現問題并在同一種背景下對各種算法的核心部分進行比較,可以將梁的應變-位移重構問題采用連續形式以下式進行表達

其中,l為梁上某點的坐標,ε0(l) 代表連續化的測量應變,ε (r,l) 是通過某種應變-位移關系建立的、與位移相關的向量形式的應變函數.err[r(l)] 是以位移函數為宗量的泛函形式的變形重構誤差函數.對于不同變形重構方法,式(1)中的 ε (r,l) 可替換為各方法中所采用的連續化的應變-位移關系函數.可以看出,問題的核心是 ε (r,l) 所體現的應變-位移關系.

在現有的平面梁變形重構研究中,被最廣泛采用的思路是利用線性化的平面梁應變、轉角與撓度的微分關系進行逐點的遞推計算,即依賴于Euler-Bernoulli 梁方程進行討論,如Ko 法等等.

此類方法,包括那些基于相似的線性微分關系,通過待定參數法[27]或逆有限元方法[28]等構建變形重構算法的方法,其共有的局限性在于由于忽略了梁上各點軸向的位移[28-29],或將軸向位移進行簡單的線性化[27],進而忽視了軸向位移對曲率的影響,或者說軸向位移與橫向位移的耦合效應.圖2 展示了這兩種情況,其中紅色代表實際的被測對象變形,黑色代表在對應的處理方式下重構出的變形狀態.由于這兩種處理方式中對變形模式的簡化,導致了變形重構的誤差.這在小變形條件下誤差可以忽略,且計算簡便,但隨變形增大誤差也顯著擴大,顯然不能適用于柔性梁結構的變形重構.針對柔性結構的變形重構經典方法是借助單純基于曲率的幾何關系,通過Frenet 標架將由應變結算的曲率、撓率和曲線切向、法向和副法向向量納入到統一的微分方程中,再由數值積分方式進行重構[30].朱曉錦等[31]采用同樣是基于曲率的方法,不同之處在于是借助運動坐標系和密切平面將曲率與形狀聯系起來.這類算法在軸向變形較大的情況下會引起嚴重誤差,在應用于空間結構時,由于其對梁變形模式的過度簡化,局部扭轉和橫向剪切也都會導致此類方法失去適用性.

圖2 基于線性應變-位移關系的變形重構Fig.2 Deformation reconstruction based on linear strain displacement relationship

因此,在有限變形的條件下,一種能夠更加精確、更加完整系統地描述測試對象變形狀態的應變-位移關系需要被引入變形重構方法之中.

2 變形重構問題的ANCF 描述

在前述的平面梁變形重構問題的形式下,通過有限元方法對被測物進行離散可以將相關問題改造成與分布式的應變測量條件相適應的離散形式.而由于已有的相關研究因其使用的線性的應變-位移關系在大變形方面存在局限性,本文通過引入ANCF單元,獲得梁在有限變形條件下系統性的應變-位移關系,繼而構建出基于ANCF 的變形重構求解問題來對已有方法做出改進.相較于已有的方法,使用ANCF 可以實現對有限變形條件下被測體應變-位移關系的精確描述.

2.1 變形重構問題的有限元離散化

為了適應分布式的應變測量條件,引入有限元法中的形函數概念,對被測梁進行分段離散并使用節點位移進行表示,其中q為整體的節點位移,qe表示單個單元的節點位移

考慮到應變測量值實際上不是連續的,實際上的各應變傳感器測量的空間范圍通常也無法簡化為點,且單個單元內應變變化應較小.因此參考A.Tessler 發展的逆有限元法[32]中的處理方式,將測量值視為對應單元內的平均應變,并設置 εe(ξ,qei) 由軸向應變與彎曲應變組成,則平面梁的應變-位移重構問題可改寫為下式

2.2 逆梯度縮減ANCF 平面索梁單元

接下來需要找到一種合適的形函數及應變表達,即確定 εe(ξ,qe) 的具體函數表達.為了適應有限變形條件下的實際需要,再考慮到計算復雜度和泛用性方面的因素,本文采用歐拉梁假設,使用縮減ANCF 平面索梁單元中應用的形函數形式來進行推導.

在縮減ANCF 平面梁單元中,單元節點坐標可表示為

單元內部的坐標插值關系可表示為

需要注意的是,如圖3 所示,梯度縮減ANCF 索梁單元因軸向應變及彎曲應變耦合會導致在此描述中單元內的物質點(圖3 中用紅點表示)在變形時向兩端密集,導致此單元中部在此描述下相對實際情況有軸向伸長的趨勢,而兩端有軸向壓縮的趨勢,即引起虛假軸向應變.而對應變進行單元內的平均可以在很大程度上避免此類耦合效應引起的單元內部畸變問題[33].這也是采用平均應變而非由特定點利用對應的形函數值進行擬合的另一個原因.

圖3 梯度縮減ANCF 梁單元內部的虛假軸向應變Fig.3 False axial strain in gradient reduction ANCF beam element

使用Green-Lagrangian 應變張量作為應變測度,并根據相關推導[34],得到如下形式的軸向變形導致的線應變和彎曲變形導致的線應變

其中f(ξ) 表示伸長率.進一步進行單元內的平均

式中的S1和S2分別為形函數S的第一行和第二行,全式是運用矩陣運算的相關定義推導整理而來.式(9)中的為伸長率的單元內均值.如果伸長率f采用精確值,對于結果精確度的貢獻較小,但對計算的復雜度影響很大,且易受之前提到的軸向偽應變效應的影響.為簡化求解過程,參照ANCF 在處理動力學問題時已有的簡化方法[35]將單元內f簡化為一由測量值直接確定的常數值

同樣是為了計算效率而作的簡化.另外為了適應軸向大變形下可能出現的梁高度方向的變化,在變形前后體積不變、同一垂直法線的截面內部受軸向變形的影響面積發生成比例放大或縮小的近似假設下,使用均勻化處理的伸長率對梁的初始高度h(ξ)進行近似校正

St及Sb分別稱為軸向應變矩陣和彎曲應變矩陣,在此種簡化下,可以避免在求解過程中為了獲得St及Sb進行復雜的數值積分運算,二者可預先計算并存儲,極大地簡化了相關求解過程,同時使應變對節點位移的函數為標準二次型格式,利于后續優化求解過程中Jacobian 矩陣與Hesse 矩陣的求取.而與通過測量的應變直接得出,同樣能使后續計算盡量簡化.

將式(9)代入式(3)進行整理,可將平面梁的變形重構問題寫為以下非線性最小二乘形式

通過式(14)可以構建出單個縮減ANCF 平面索梁單元內從應變到位移的推導關系,本文將由此種關系建立起的單元格式稱之為逆縮減ANCF 平面索梁單元.

3 基于逆ANCF 單元的變形重構方法

顯然,通過逆縮減ANCF 平面索梁單元對變形重構問題進行求解是一個非線性優化問題的求解過程,需要借助相關數值優化方法建立一個可靠的求解算法.由于式(14)所提的變形重構問題的非線性特征與其已知量和決策變量間數目間的不匹配特性,本問題符合文獻[36]中的不適定性數值最優化問題的定義,此類問題的根源是解的確定性不足,問題的提法不夠完備.問題的不適定性在求解時表現為無法正確收斂和收斂過程的不穩定,例如出現目標函數Hesse 矩陣不正定的情況.對于一般的不適定性最優化問題可以通過增加適當的限制使之改造為適定性問題,而對于本問題,可以通過引入縮減ANCF 單元外的、由實際情況帶來的限制以使本問題實現穩定的求解.所引入的關系分別為節點自由度間自相關性條件、節點處的C2連續性條件以及邊界條件,這些條件的引入同時也能使本方法在單元密度較低時取得較好的效果.而在求解過程中,由于不同工況對求解算法提出了不同的要求,通過建立逐單元和多單元兩種求解格式,并通過遞推的方法組裝單元,提升了本方法對不同應用工況的泛用性.

3.1 對逆ANCF 單元的修正

之前提到了將單元內應變平均化會導致一些問題,最主要的就是因待求解單元自由度與測量信息數目上不匹配導致的不適定性問題.針對此問題,并結合變形重構問題的特點,下文進行了討論并給出了修正方法.

3.1.1 簡化節點自由度

首先需要注意的是,節點的4 個坐標分量中,平動坐標對自然坐標的偏導數項不是獨立的.如果直接使用節點四自由度的形式作為數值最優化求解的決策變量會導致考慮其自相關約束的求解列式較為復雜,也影響計算效率.此外,偏導數項的物理意義不夠明確,這對邊界條件的參數化和后續在載荷感知等領域的擴展也不利.因此本文使用以下縮減的節點位置向量進行替代

原始節點位置向量與縮減的節點位置向量二者間的聯系可以用下式表示

式中f(0)及f(1) 兩項分別代表著兩個節點處的伸長率,在實際求解時均使用簡化的、由應變測量值可以直接給出的fˉ 替代.在進行后續的優化求解時,由于目標函數與決策變量間的微分關系是迭代過程中的主要依據,上式表示的二者聯系可以使縮減的節點位移向量在進行迭代完全替代原來的節點位移向量.

3.1.2 邊界約束及曲率連續性約束

為了唯一地確定在某一邊界條件下的被測對象的變形,需要同時引入邊界條件和曲率連續性條件.利用外點罰函數方法可以方便地實現這一點.一個典型的受邊界約束和曲率連續性約束的逆ANCF 單元如圖4 所示.

圖4 邊界約束及曲率連續性約束示意圖Fig.4 Schematic diagram of boundary constraints and curvature continuity constraints

記邊界條件為

在分段求解時,可以如圖4 通過已經求解并確定位移的節點對未經求解的相應單元施加邊界條件.

ANCF 單元在節點處只具有C1連續性,這對利用基于ANCF 的逆有限元方法重構被測對象完整連續的應變場是不利的.如圖4 所示,可以利用施加額外約束條件的方法迫使本方法中ANCF 單元節點處具有C2連續性,即保證除平動坐標和轉角外,曲率和彎曲引起的應變在相鄰單元間也是連續的,以消除ANCF 單元的這種局限性.將前一個單元的末端曲率作為下一個單元起始點曲率的目標值,同樣通過外點罰函數形式進行約束條件的施加.記前一個單元的末端曲率為,則相關限制條件可以表示為

加入曲率連續性限制條件的算法不僅解決了適定性問題,同時相較于經典的ANCF 單元,在節點處C2連續的本方法顯然更加精確且接近實際.這帶來了能夠以較少的測點與單元進行精確的變形重構的優勢.

這里需要說明的是,部分求解分段可能沒有已求解共節點分段的節點曲率信息,例如最初求解的分段.這時它們的曲率約束條件可以根據實際條件給出,例如簡支端和無彎矩作用的自由端曲率可設為0.對于通常沒有作用于內部的大載荷因而曲率分布較為均勻的固支端單元,可將由其測量應變直接轉化為的曲率作為端部曲率,這樣既保證了計算的穩定性也不會引起大的誤差.在此基礎上,還可以使固支端單元短一些,亦即使對應的測點靠近端部,能夠進一步保證這種近似的合理性從而避免誤差.

3.2 逐單元求解過程與多單元整體求解過程

綜上,利用經過上述修正的逆ANCF 單元,可以構建出平面梁變形重構的方法.鑒于整體求解時決策變量維度過大,準確收斂難度較大,考慮將被測體進行分段求解,并利用ANCF 單元的特性構造遞推組裝算法.基于此,建立了逐單元求解和多單元整體求解兩種求解算法.其中逐單元求解算法由于每次求解的決策變量維度小且每一迭代步的時間復雜度約為,計算效率通常高于多單元求解算法;而多單元整體求解方法可以利用已知的邊界條件對求解進行校正,以防止誤差在長度范圍上的積累.二者實際可以在同一測量工作中混合使用,以克服二者各自的局限性.

3.2.1 逐單元迭代求解過程

對于單個單元,將測量的應變信息 εt,εb代入式(14)中,并將已求解共節點分段的共用節點處曲率κ0和位置信息re|ξ=0通過式(18)和式(19)引入,可以得到基于單個單元的變形重構最小二乘列式

從而求出此單元未求解節點處的位移re|ξ=1,如圖5所示.

圖5 逐單元變形重構Fig.5 Shape reconstruction in a single element

式中 σBC及 σAD分別為邊界約束條件和曲率連續性條件對應的罰因子,理論上會影響方法的收斂速度和收斂性.但在數值仿真中證明,本文涉及的算例中二者取值在10-1～102間均能取得較好的效果,在保證精度的情況下二者的取值范圍相當大,即本方法的變形重構效果對罰因子的取值不敏感.本文之后的數值仿真中罰因子項均取為1.

在單個單元條件下進行如上式描述的最小二乘列式迭代求解復雜度有限,最小二乘列式的梯度向量和Hesse 矩陣形式較為簡單,因此可以使用基本的Newton 法進行穩定便捷的求解.

則殘差的Jacobian 矩陣為

進而有

最小二乘列式中其他兩項罰函數的計算過程均類似,結合式(18)及式(19),分別通過邊界條件或一端的曲率求出各項的殘差并結合各項Jacobian 矩陣進行計算,此處不贅述.

因此總體最小二乘列式的相關參數可表示為

則Newton 法列式可記為

其中k表示第k次迭代,dk表示第k迭代步.在實際數值仿真中,逐單元的變形重構過程中Newton 法在各類工況下都相當穩定,不會出現奇異等情況.

對于單個單元,通常在求解時已知的曲率和其他邊界條件能夠讓對單個單元進行的變形重構問題變為適定的.

3.2.2 多單元整體迭代求解過程

多單元求解算法的優勢在于能夠適應分散的邊界條件,在處理此類問題時能結合邊界條件進行求解.另外對于迭代求解過程中剩余量較小的情況,這種算法在某種程度上能夠減少迭代的次數,這預示著此種求解算法在對被測物變形情況進行動態跟蹤時具有一定優勢.

對于多個單元同時求解的情況,借鑒有限元方法中剛度矩陣的組集方式,引入布爾映射矩陣Bi建立整體坐標至局部坐標的映射.將各單元測量的應變信息 εti,εbi代入式(14)中同時將各單元由此得到的罰函數項求和,并將求解共節點分段的共用節點處曲率 κ0和位置信息re1由式(18)和式(19)引入,即可得到多單元整體求解列式

從而求出此分段中各未求解節點的坐標rei,如圖6所示.

圖6 多單元變形重構Fig.6 Shape reconstruction in multi element

在式(28)中,整體坐標的形式為

邊界條件罰函數項PmBC及曲率連續性條件罰函數項PmAD為

式(28)中的基本最小二乘列式Fm0(q) 為

式中

將式(32)寫成矩陣形式,即可將殘差記為

殘差的Jacobian 矩陣為

邊界條件罰函數項PmBC及曲率連續性條件罰函數項PmAD的梯度向量及Hesse 陣以相同形式寫出,此處不贅述.

由于整體求解算法下決策變量維度不可避免的增大,在大剩余量條件下收斂穩定性往往較差,很多時候需要采用更加復雜、魯棒性更高的最優化求解方法才能進行順利求解,此處不作深入討論,后文中整體求解方法的實現使用由Dogleg 方法求解子問題的信賴域方法.

3.3 單元遞推組裝

進行分段求解之后,還需要一種方法將被測對象的各自獨立的分段連接起來.已經知道每個節點的3 個坐標分量(r1,r2,θ)都是線性獨立的,同一節點的3 個坐標分量互不影響,且具有線性可加性,各段內部的變形實際上決定了從一個節點指向另一個節點的位移向量.如圖7,經過方向余弦矩陣進行坐標轉換處理以后,得到的位移向量從一個節點坐標指向另一個節點坐標.因此,可以在對各段分別進行求解之后再進行連接,由于之前提到的各坐標分量間的獨立可加性,這種連接不會破壞一維二節點ANCF 單元在節點處的連續.

圖7 單元遞推組裝Fig.7 Unit recursively assemble

具體算法是將每段均視為受一端固支的邊界條件限制,根據測量應變和已進行求解的共節點相鄰段的末端曲率,進行自由端位移的求解.當然,這還依賴于事先規定的遞推求解順序.

通過這種形式的求解,能夠得到從固支端指向自由端的位置向量,記為

進行逐段拼接,即依次確定各節點的位置向量,記為

拼接時,其中的平動坐標分量需要經過方向余弦矩陣進行旋轉處理,角度可以直接相加,計算如下

由此,就獲得了逐節點的求解結果.如果需要整體連續的變形狀態,可以將形函數代入,即可得到連續的位移函數.需要注意的是,此處得到的是基于初始構型和自然坐標的位移描述.

由于實際應用中具體情況的不同,以上單元組裝方法可以根據需要而預先設計規劃并靈活使用.

3.4 算法流程

綜合以上,將整體的算法流程進行總結:

(1)將被測對象進行單元劃分.在這一過程中,要考慮工作中載荷位置及大小對應變分布的影響及應變測量設備的實際限制.單元及測點密度選取的標準是使在預計的變形范圍內單元內的應變都可以通過本文使用的ANCF 方法較好地描述,這也意味著在應變變化較為劇烈的部分應增大測點密度才能使變形重構效果理想;

(2)根據之前對逐單元求解方式和多單元求解方式的各自特點及適應范圍的討論,進行各獨立求解分段的劃分;

(4)根據測量的應變信息,各求解分段按順序應用其對應的求解迭代算法和事先設定的收斂終止規則,在依次傳遞節點曲率信息的基礎上進行獨立的求解.其中單單元分段由式(20)所示的最小二乘列式進行求解;多單元分段由式(28)所示的最小二乘列式進行求解;

(5)借由式(39)進行整體坐標的組裝.

4 數值仿真

本節進行本方法的仿真驗證.下述全部算例所需的原始數據,包括應變及位移均由ABAQUS 計算給出.即對圖8～ 圖12 所示的仿真工況,即均使用ABAQUS 進行仿真分析,從而得到各個工況下的應變及對應的位移數據,將由應變重構出的位移數據與仿真得到的位移數據相比較.下列算例的ABAQUS分析均使用B22 單元,進行幾何非線性條件下的靜態通用分析步.此算例運行于MATLAB R2020a 軟件環境下,處理器型號為i7-10700.

圖8 位移載荷下的外伸梁示意圖Fig.8 Schematic of overhanging beam under displacement load

本節使用3 個算例來展示并驗證本算法在有限變形條件下的變形重構效果.首先是在一個給定邊界條件下,將逐單元求解過程和多單元求解過程混合使用,進行不同單元密度下的效果驗證,以初步展示本方法的重構效果;第2 個算例是在一端固支、自由端受橫向載荷的簡單載荷下,將不同的經典平面梁變形重構方法進行對比,驗證本方法在大變形條件下的優越性;最后一個算例將本方法與基于曲率連續化的方法進行對比,以驗證本方法在軸向大變形條件下具有更佳的適應性.

4.1 大變形工況下的重構效果驗證

為了驗證本方法中邊界條件相關部分的有效性,設計了一種分散邊界條件下的變形工況,如圖8所示,以一根40 cm 的細長梁為對象,在自由端施加豎直方向的位移載荷u.此被測對象一端簡支,在其中部20 cm 處進行豎直方向的鉸接約束.

為驗證本方法在不同測量條件下的性能,模擬使用不同的測點個數及對應的單元密度進行建模和求解,在各種單元密度下,都采用均勻劃分單元的方式.此處對左側兩處鉸接中的部分進行整體求解,并將邊界條件使用式(18)的形式代入求解過程中,右側自由部分則進行逐單元求解.逐單元求解部分的各單元間及其與整體求解部分的連接處使用單元組裝方法進行連接,為了直觀地看出重構效果,為重構得到的節點坐標定義如下的相對誤差表征

其中n為節點數,di為重構得到的各節點坐標分量,為仿真得到的對應位移分量,為在整個被測對象上對應項的模的均值.

本算例與后續算例以兩次迭代殘差間的差值小于設定值為迭代終止條件.

圖9 中依次給出了20 測點、10 測點、4 測點數下施加不同位移載荷時本方法重構出的變形狀態(圖9 中以散點表示的節點坐標)與仿真得到的變形狀態(圖9 中以實線表示)的對比.為模擬平面梁的大變形情況,將位移載荷分別設置為50 mm,100 mm,200 mm,其中200 mm 的位移載荷已達到被測梁長度的一半,是相當極端的變形情況.從圖9(a)中可以看出,在20 測點時,各種載荷條件下本方法由應變重構出的變形狀態均與仿真得到的變形狀態均基本重合,這說明本方法能夠較為精確地實現有限變形條件下的變形重構;在10 測點(圖9(b))下,重構效果依然良好;而在最為極端的4 測點(圖9(c))下,位移載荷在50 mm 和100 mm 情況下的重構效果仍然很好,只有位移載荷200 mm 時存在較明顯誤差.正如前面分析,誤差產生的原因是4 個測點已經不足以精確反映被測對象內部的應變分布.上述結果表明,本文方法對不同程度的大變形重構都具有較高的精度和收斂性.

圖9 大變形條件下不同測點數的重構效果Fig.9 Reconstruction performance of different number of measuring points under large deformation

進一步在表1～ 表3 給出各仿真工況下重構的相對誤差及計算用時.從表1 可以看出,在各位移載荷下,總體的重構相對誤差接近,200 mm 位移載荷工況比另兩種位移載荷工況相對誤差略微增大;在計算時間方面,各位移載荷工況計算用時也都相近.這證明了本方法在不同變形條件下的精度和收斂性.表2 除同樣體現上述規律外,可以看出10 測點數與20 測點數條件下重構的相對誤差及計算用時仍然差距不大,這說明在合適的測點數范圍內,測點數對重構效果影響有限.表3 中展示了在較為極端的4 測點條件下的仿真情況,此測點數下重構相對誤差較之前有較大增加,特別是200 mm 位移載荷條件下出現較明顯的誤差,而計算用時方面則變化很小.如之前分析,這說明4 測點已經不足以精確反映被測對象內部應變分布,而200 mm 位移載荷工況下相對誤差顯著高于此測點數下另兩種位移載荷工況,則證明此種原因導致的誤差對于變形程度大的情況更加顯著.在計算用時方面,不同測點數及位移載荷工況下用時均較為接近,這說明本方法在各工況下均具有較好的收斂性.

表1 20 個測點下的相對誤差及計算用時Table 1 Reconstruction error and calculation time under 20 measuring points condition

表2 10 個測點下的相對誤差及計算用時Table 2 Reconstruction error and calculation time under 10 measuring points condition

表3 4 個測點下的相對誤差及計算用時Table 3 Reconstruction error and calculation time under 4 measuring points condition

單純從單次重構所需的時間可以看出,本方法計算需時實際較短.考慮到以相近的變形狀態作為初值能夠進一步縮減迭代步數,以及儲存的軸向應變矩陣及彎曲應變矩陣在同一被測對象的不同變形狀態下無需重復計算,本方法在變形狀態動態追蹤方面應該也具有潛力.

4.2 與已有方法的對比

之前提到了線性逆有限元法、Ko 法以及基于曲率連續化的方法的特點及其局限性,在此處以較為簡單的懸臂梁受自由端較大數值方向位移載荷時的工況,給出以上方法與逆ANCF 單元方法的實際效果的比較.使用圖10 所示的仿真設置進行不同變形重構方法的對比.

圖10 位移載荷下的懸臂梁示意圖Fig.10 Schematic of cantilever beam under displacement load

本算例下各種仿真條件下的不同方法均采用20 測點,單元以均勻方式布置.具體效果如圖11.

圖11 不同方法的懸臂梁變形重構仿真結果Fig.11 Reconstruction performance of cantilever beam in different methods

圖11 中給出了施加不同位移載荷時本方法、曲率連續化方法、Ko 方法和逆有限元重構效果的對比.從圖中可以看出,Ko 法與逆有限元法由于其線性化特性及忽略了軸向位移影響,隨變形程度增大顯著提升.逆有限元方法相較于用分段去對橫向位移進行近似非線性計算的Ko 法而言誤差更大.這充分說明了在有限變形條件下基于線性假設的變形重構方法的局限性.曲率連續化方法和本方法在此仿真算例的各種位移載荷設置下重構得到的變形狀態均能與仿真給出的變形狀態高度重合.

在此種軸向變形不顯著的情況下,單純基于曲率的方法與逆ANCF 單元法均可以準確地重構出被測體的變形,效果差異差距不顯著.因此設計如圖12的仿真驗證軸向大變形情況下ANCF 方法的表現.

圖12 軸向大變形條件下的簡支梁示意圖Fig.12 Schematic of simply supported beam under large axial deformation

此算例依舊使用20 測點,其中軸向和橫向的位移載荷設為相等,均為u.結果見圖13,隨橫向和軸向位移載荷的增大,被測對象的軸向變形愈加顯著,此時,曲率連續化方法的重構誤差明顯增加,而本方法仍然能保證較高的精度.可以看出,本方法對軸向大變形情況的適應性明顯高于單純基于曲率的方法.

圖13 軸向大變形條件下曲率連續化方法和iANCF 方法的對比Fig.13 Comparison of curvature-based method and iANCF under large axial deformation

5 結論

本文繼承了逆有限元法的思想,利用非線性的縮減ANCF 平面梁單元對梁的應變-位移映射關系進行描述,進而采用數值優化方法的思想進行位移的求解.將求解格式分為逐單元與整體兩種,前者求解穩定且計算量小,后者能夠實現對邊界條件的靈活適應.而通過單元的連接算法,二者實際上可以在實際應用中相互融合.

通過多個數值算例證明了本算法的有效性,在大變形條件下能夠用較少的測點實現精確的變形重構.此外,相較于已有的柔性結構變形重構算法而言,本方法引入梯度縮減ANCF 單元來刻畫被測對象的應變場和位移場間的關系,進而成功地將能夠實現精確變形重構的變形模式范圍擴展到了能用ANCF 方法描述的任意變形模式.結合已有的不同類型ANCF 單元,本文提出的利用ANCF 單元構造優化問題進而建立變形重構算法的思想能夠推廣到不同類型的被測對象及場景.而由于本方法與有限元方法的高度融合,也能夠自然地將應用范圍擴展到與有限元相關的動靜載荷識別、振動控制及健康監測方面.