吳懋琦 譚述君 , 高飛雄

* (大連理工大學工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室,遼寧大連 116024)

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

引言

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

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

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

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

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

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

1 有限變形下的平面梁變形重構(gòu)問題

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

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

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

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

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

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

圖2 基于線性應變-位移關(guān)系的變形重構(gòu)Fig.2 Deformation reconstruction based on linear strain displacement relationship

因此,在有限變形的條件下,一種能夠更加精確、更加完整系統(tǒng)地描述測試對象變形狀態(tài)的應變-位移關(guān)系需要被引入變形重構(gòu)方法之中.

2 變形重構(gòu)問題的ANCF 描述

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

2.1 變形重構(gòu)問題的有限元離散化

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

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

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

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

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

單元內(nèi)部的坐標插值關(guān)系可表示為

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

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

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

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

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

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

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

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

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

3 基于逆ANCF 單元的變形重構(gòu)方法

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

3.1 對逆ANCF 單元的修正

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

3.1.1 簡化節(jié)點自由度

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

原始節(jié)點位置向量與縮減的節(jié)點位置向量二者間的聯(lián)系可以用下式表示

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

3.1.2 邊界約束及曲率連續(xù)性約束

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

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

記邊界條件為

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

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

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

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

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

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

3.2.1 逐單元迭代求解過程

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

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

圖5 逐單元變形重構(gòu)Fig.5 Shape reconstruction in a single element

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

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

則殘差的Jacobian 矩陣為

進而有

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

因此總體最小二乘列式的相關(guān)參數(shù)可表示為

則Newton 法列式可記為

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

對于單個單元,通常在求解時已知的曲率和其他邊界條件能夠讓對單個單元進行的變形重構(gòu)問題變?yōu)檫m定的.

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

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

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

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

圖6 多單元變形重構(gòu)Fig.6 Shape reconstruction in multi element

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

邊界條件罰函數(shù)項PmBC及曲率連續(xù)性條件罰函數(shù)項PmAD為

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

式中

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

殘差的Jacobian 矩陣為

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

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

3.3 單元遞推組裝

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

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

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

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

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

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

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

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

3.4 算法流程

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

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

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

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

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

4 數(shù)值仿真

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

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

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

4.1 大變形工況下的重構(gòu)效果驗證

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

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

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

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

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

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

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

表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

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

4.2 與已有方法的對比

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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