固體火箭發動機對接裝配密封圈應力在線預測方法①

張家川,徐志剛,王軍義*,楊 嘯,王元鈺,侯賡舜,董祺成

(1.沈陽航空航天大學 機電工程學院,沈陽 110136;2.中國科學院 沈陽自動化研究所 機器人學國家重點實驗室,沈陽 110016;3.中國科學院 機器人與智能制造創新研究院,沈陽 110169)

0 引言

固體火箭發動機的密封圈位于燃燒室與噴管接口的密封槽及止口處,通過燃燒室與噴管的對接安裝,使密封圈產生擠壓變形,達到密封目的。密封圈作為發動機貯存過程中的薄弱環節和極少數的可更換部件,密封結構的壽命是確定發動機首翻期的主要因素之一[1-2]。在安裝“O”形橡膠密封圈時,橡膠圈的接觸應力不足而導致密封失效,但等效應力過大將造成表面裂紋,引發深層次的破裂[3],從而導致發動機在貯存過程中密封提前失效。而大型薄壁燃燒室殼體的結構決定其剛性較差、抵抗變形能力弱、屈強比大,裝配過程中受重力影響,殼體及接口會產生明顯變形[4],并且因存在加工誤差,同樣會使燃燒室與噴管的多層臺階接口產生復雜不規則變形。另外,由于在對接過程中密封圈狀態不可見、不可測,裝配位姿可能時刻變化調整。在這些復雜的裝配條件下,如何檢測裝配質量成為一大難題。

為解決此類問題,一般采用對裝配過程進行有限元仿真的方法,但有限元仿真速度慢,計算量大,仿真周期長,無法達到快速獲取裝配質量的目的。而代理模型則可以解決此類問題,前期進行一定數量的仿真計算獲取樣本數據集,通過樣本數據集建立裝配狀態與密封圈應力-應變之間關系的代理模型,后期裝配過程中即可利用代理模型快速預測密封圈的應力-應變,為裝配質量是否達標提供參考依據。對于產品的在線預測如健康、狀態與壽命預測,可分為基于機理模型、基于知識和基于數據驅動三類[5]：基于知識的預測方法主要是憑借檢測數據與歷史數據進行對比,或通過歷史數據等建立專家知識庫來實現預測[6];基于機理模型的預測方法主要依靠對失效機理建立數學模型,分析失效形式,并依靠數學模型推導出產品的失效概率密度函數或可靠度-時間函數等[7-9];基于數據驅動的預測方法則是利用機器學習方法或數據統計方法對設計、檢測、試驗等數據進行分析,進而建立出預測模型[10]。對于代理模型的選取,一般選用神經網絡、多項式響應面、徑向基函數等方法[11-12],各種預測方法面對不同情況時各有優劣,在設計范圍不可縮減的情形下,多項式響應面法不適用于強非線性和高維模型,神經網絡預測方法可以彌補多項式響應面法的缺點,但需耗費大量的運算時間[13]。Kriging模型在樣本數據空間較小的情況下能光滑擬合出原目標函數,可以保證預測模型的精度,并且運算時間相對較短[14]。

本文采用Kriging模型作為密封圈應力的預測模型進行研究。使用有限元軟件仿真計算出不同裝配條件參數下密封圈的應力-應變,并由此構建預測模型的樣本集,對樣本集數據進行擴充及抽樣處理,將處理后的數據構建Kriging預測模型,進而利用代理模型對固體發動機的裝配質量給出判斷。

1 預測模型樣本集的建立

1.1 確定裝配條件參數

不同的裝配條件會造成密封圈不同的應力分布。確定裝配條件參數即確定預測模型的輸入項,即確定影響密封圈應力-應變大小的因素。將影響密封圈應力大小的參數歸納為兩類：一是裝配過程中的裝配位姿;二是燃燒室或噴管本身因制造偏差和受裝配力的影響而造成的接口處不規則形變。

空間中位姿存在六個自由度,分別是右手笛卡爾坐標系中沿三個坐標軸方向的平移,以及繞三個坐標軸的旋轉,即翻滾、俯仰、偏航。由于對接過程中所產生的翻滾動作很小,所以只將空間位姿簡化為角度變化及偏移變化。考慮到密封接口處為完全軸對稱結構,因此對于接口處的不規則形變,將其簡化為沿圓形接口法線方向上的跳動,假設接口處形狀變化連續,呈橢圓狀,那么可用接口的橢圓度大小來表示形變對密封圈應力-應變造成的影響。圖1簡要示意了內插式結構在裝配時的主要裝配條件參數,其剖面視圖表達了插入件的徑向變形。

圖1 裝配條件參數示意圖Fig.1 Assembly condition parameters

三個裝配條件參數中,d為燃燒室軸線與噴管軸線間的偏心距,分為dy和dz;α為軸傾角,即燃燒室軸線與噴管軸線之間的夾角,分為αy和αz;e為變形后接口的相對橢圓度,其表達式如下：

e=(dzn-dyn)-(dzc-dyc)

(1)

式中dzn、dyn分別為噴管接口處z向和y向上的直徑大小;dzc、dyc分別為燃燒室接口處z向和y向上的直徑大小。

1.2 密封圈應力有限元分析

密封圈的材料為橡膠。橡膠材料是一種分子構成較復雜的高分子材料,是由多個重復單元鏈接而成的長鏈分子,屬于非線性、超彈性材料[15]。橡膠本構關系的研究較多地采用唯象理論描述其應力-應變關系,唯象理論假定橡膠在未變形狀態下是各向同性的,并認為橡膠是不可壓縮材料[16],而對于超彈性材料的本構關系模型一般有Mooney-Rivlin模型、Yeoh模型和Ogden模型等[17]。對于150%以內的小應變而言,Mooney-Rivlin模型可得到合理的近似,而Yeoh模型更適合大變形行為[18]。對于橡膠的2參數Mooney-Rivlin模型應變能密度函數如下：

(2)

(3)

I2=(λ1λ2)2+(λ2λ3)2+(λ3λ1)2

(4)

λi=1+εi(i=1,2,3)

(5)

式中W為應變能密度;C10和C01是Rivlin系數;I1和I2分別為材料的第1、2階Green應變張量不變量;D1為材料的不可壓縮率;J為體積比,對于不可壓縮材料J=1;λi為3個方向上的主伸長率;εi為材料的應變。

密封圈的橡膠材料參數如表1,其應力-應變曲線見圖2。在確定了密封圈材料參數后,對發動機對接的過程進行有限元仿真,繪制燃燒室、噴管、密封圈的三維模型,為了縮短仿真周期,將模型進行簡化：由于裝配時燃燒室由吊裝與支架完全固定,噴管由六自由度并聯平臺固定,對接過程中噴管擴張段及矢量組件等對密封圈受力的影響較小,所以對于噴管組件只保留噴管接口處以及噴管喉部的部分;對于燃燒室,由于薄壁燃燒室更易受重力影響,所以將工裝支架定位位置到燃燒室接口的部分進行保留。在設計域內均勻選取幾種典型的裝配參數,利用ANSYS WORKBENCH有限元仿真軟件對對接過程進行仿真計算。

圖2 密封圈應力-應變曲線圖Fig.2 Stress-strain curve of the sealing ring

表1 密封圈材料參數Table 1 Material parameters of the sealing ring

根據圖3(a)所示發動機模型的密封結構,仿真的接觸設置如下：

(1)燃燒室內表面與噴管的外表面構成接觸;

(2)燃燒室外端面與噴管止口端面構成接觸;

(3)密封槽處的密封圈與噴管密封槽的底面和兩個側面構成接觸;

(4)止口處的密封圈與噴管止口的軸向表面和徑向端面構成接觸;

(5)密封圈與燃燒室內表面及外端面構成接觸。

將噴管與燃燒室的接觸類型選擇為摩擦,其中噴管外表面為接觸面,燃燒室內表面為目標面。由于存在預壓縮,密封圈與噴管之間的摩擦很大,所以將密封圈與噴管之間的接觸類型均選為粗糙,密封圈與燃燒室之間的接觸類型定為摩擦,其中密封圈表面為接觸面,噴管及燃燒室表面為目標面。密封圈的彈性模量E=6.54 MPa,剪切模量G=1.85 MPa,由于橡膠為不可壓縮材料,取泊松比ν=0.5。由于其他材料參數為定值,而仿真的摩擦系數因實際裝配工況復雜而難以精準確定,且仿真的摩擦系數的選取對仿真結果的準確性有一定的影響,所以采取先根據接觸材料,按經驗選擇摩擦系數;再根據仿真結果修正選取最優摩擦系數的方式,使仿真結果貼近試驗結果。最終選取噴管與燃燒室之間的接觸摩擦系數為0.1,密封圈與燃燒室之間的接觸摩擦系數為0.4。

接觸算法采用增強拉格朗日法,接觸行為采取不對稱行為。由于噴管及燃燒室的剛度遠高于密封圈,所以主要考慮密封圈的變形,將密封圈的網格大小設置為0.4 mm,共13 473 181個節點,3 616 757個單元,噴管及燃燒室網格大小為20 mm,密封圈網格使用四面體單元,單元物理環境選擇非線性力學。為簡化仿真分析設置,將噴管底面施加固定支撐,并對所有幾何體施加重力約束。將載荷步分為兩階段：第一階段對密封圈施加法向位移模擬密封圈的預壓縮過程;第二階段對燃燒室底部施加位移模擬對接裝配過程。

對結構進行有限元分析非線性問題包括材料非線性、幾何非線性和狀態邊界非線性[19]。密封圈安裝于噴管的密封槽和止口處,密封槽內的密封圈在裝配過程中受到剪切擠壓,密封圈會逐漸填充占滿密封槽底面,其邊界條件會發生變化,存在狀態邊界非線性;密封圈的橡膠材料為超彈性體,本身存在材料非線性;且超彈性體在受力時會產生大應變及大位移,存在幾何非線性。

密封圈載荷分兩階段加載,而且非線性較為強烈,計算不容易收斂。為保證分析計算可靠、節約仿真周期,在仿真設置中采取混合U-P單元防止模型體積自鎖,同時采用非線性自適應網格法,并開啟大變形開關,并利用重啟動分析,幫助計算收斂。密封圈仿真結果見圖3(b)。

(a)3D model of rabbet section plane

2 擴增數據集

因為有限元仿真計算十分耗費時間,仿真出大量樣本顯然是不現實的,而數據擴充方法可以有效解決數據量不夠的問題。對于數據擴充,或數據增強,有仿射變換、生成對抗網絡(Generative Adversarial Networks,GAN)[20]等方法。對于密封圈預測模型,選用生成對抗網絡模型來擴充應力數據集,生成對抗網絡可以較好學習到原始的數據分布,并創造出新的、與原始數據分布相近的數據。

如圖4所示,生成對抗網絡由生成G和判別器D組成,生成模型通過輸入噪聲偽造具有樣本分布的數據,判別模型用于區分真實數據和偽造數據[21]。

圖4 GAN模型結構示意圖Fig.4 Structure of GAN model

整個過程相當于極大極小博弈,生成模型和判別模型相互對抗,相互優化,并不斷提升自身的生成能力及判別能力,最終判別模型無法辨別生成模型所生成的數據真偽時,生成模型和判別模型達到了納什平衡。

生成模型與判別模型的構成均由卷積神經網絡(Convolution Neural Network,CNN)構成,隨機噪聲z通過生成模型G輸出G(z),盡可能地使G(z)的分布接近于真實樣本x的分布p(x)。而判別模型則是要判斷生成的數據是否屬于真實樣本x的分布,若數據符合樣本分布,則輸出1,否則輸出0。所以生成模型需要一直生成數據使自己接近真實樣本,直到判別模型無法判斷其是否為生成的數據,其目標函數為

Ez～pz(z)[ln(1-D(G(z)))]

(6)

式中pdata(x)為真實樣本x的分布;p(z)為噪聲z的先驗分布;E(·)為期望。

然而,傳統GAN只能學習一種數據的整體,不能生成指定的數據分布,對于多個類別的數據學習能力較差,且效率低。而條件生成對抗網絡(Conditional Generative Adversarial Network,CGAN)[22]則可以生成不同類別的數據。CGAN的目標函數如下：

Ez～pz(z)[ln(1-D(G(z|y)))]

(7)

CGAN在GAN的基礎上做了一些修改,添加了標簽y作為約束信息,生成模型由噪聲z和標簽y作為輸入,判別模型判別生成數據G(z|y)和真實樣本x相似度并通過標簽y來調節生成模型所生成的數據。

CGAN模型一般用于圖像及視覺等領域,處理表格類數據時一般使用表格數據生成對抗網絡(Tabular Generative Adversarial Networks,TGAN)[23]。條件表格數據生成對抗網絡(Conditional tabular generative adversarial network,CTGAN)[24]將CGAN與TGAN結合起來,可以用條件分類的形式處理表格數據。采用CTGAN來處理仿真數據集,以處理不同裝配條件參數下的應力數據。所使用的CTGAN的生成器由兩個全連接層和兩個反卷積層構成,判別器由兩個卷積層和兩個全連接層構成。

由于在裝配位姿改變及接口變形時不同位置的密封圈單元所受應力區別較大,所以將不同位置的密封圈單元作為分類依據,按不同位置進行密封圈應力的數據擴充。

3 基于Kriging的應力預測模型

由于GAN所生成的數據精度有限,不能保證其與真實仿真計算的數據完全相同,同時防止過擬合等問題,建立Kriging預測模型之前先要對樣本在參數邊界范圍內進行采樣,常用的實驗設計方法有正交實驗法、均勻實驗法、交叉實驗法、隨機區組法、蒙特卡洛法、拉丁超立方法等。由于拉丁超立方法(Latin Hypercube Sampling,LHS)可以用少量數據表示全局較多的特征,所以使用拉丁超立方法選取樣本點。

3.1 拉丁超立方抽樣

抽樣的結果(圖5)需要保證所選樣本點能夠表達整個空間的特性,不損失模型特征;而且也不宜冗余點過多,使得預測模型過擬合,影響模型精度。拉丁超立方采樣優于隨機采樣的原因是它確保采樣值能夠覆蓋輸入隨機變量的整個分布區間[25],使抽樣點更具有代表性。

圖5 拉丁超立方抽樣示意圖Fig.5 Latin hypercube sampling

因為各個裝配參數狀態下應力數據均勻分布,且重要性等同,所以采用均勻分布的拉丁超立方抽樣。拉丁超立方抽樣的原理是,從p維空間中抽取n個點。具體做法是,將每一維等分成n個不重疊的區間,這n個區間概率相同,從每一維的區間中選取點,再將各維中所選取的點在空間中隨機組合。

3.2 Kriging預測模型的建立

Kriging模型[26]是最早由南非礦業工程師KRIGE D G提出的一種空間自協方差最優插值法,是一種利用已知點信息來預測未知點上響應的無偏估計模型。Kriging模型與響應面法有類似之處,但響應面法對高維非線性的數據擬合精度較差。而且Kriging模型相比于神經網絡、響應面法等模型而言,可在數據點較少的情況下獲得較高的精度[27]。Kriging相比于其他插值方法,只使用估計點附近的信息,而不是全部信息;而且Kriging模型同時具有局部和全局的統計特性[28]。另外,Kriging方法很容易實現局部加權插值,這樣就克服了一般距離加權插值方法插值結果的不穩定性[29]。

利用Kriging模型對固體發動機對接裝配時的密封圈應力進行預測,由于應力預測與應變預測的方法流程一致,所以下文均用應力代替。選用經過拉丁超立方抽樣出的樣本集作為代理模型輸入的初始樣本,包含偏心距、軸傾角、橢圓度和應力(應變)值。將初始樣本代入到Kriging模型中進行擬合,其本質流程是在估計方差最小和最小估計變異數的約束條件下,計算樣本點的權重系數,將樣本點和權重進行線性組合,求得任意點的內插估計值[30]。

Kriging的理論模型是由多項式回歸模型和隨機部分共同組成：

(8)

其中R(…)為空間相關方程,其決定了預測模型的精度,一般使用高斯相關方程等[31]。高斯相關方程具體形式為

(9)

式中n為樣本點數;θj為相關性參數。

對于樣本點x的樣本集X=[x1,x2,…,xm]T的預測輸出值Y=[y1,y2,…,ym]T,樣本點的預測輸出值即為

(10)

其中,c=c(x)∈Rm為預測方差系數。預測輸出的誤差為

(11)

式中Z=[z1,…,zm]T為設計域的誤差集;F為樣本點處多項式f(x)的估計值向量。

為保證無偏估計,需要滿足：

FTc(x)=f(x)

(12)

即FTc-f(x)=0,此時預測誤差的最后一項為零。在此情況下,預測均方差(Mean Squared Error,MSE)為

(13)

x處的響應值可表示為

(14)

其中,Rγ*=Y-Fβ*;R為樣本設計域與隨機誤差函數z之間的隨機過程相關方程矩陣;r為由隨機過程相關方程組成的樣本點x間的相關性矩陣;β*為無偏條件下的回歸系數,由關于R的廣義最小二乘法得到。

滿足無偏條件的Kriging模型預測值的MSE值為

(15)

式中u=FTR-1r-f。

對應方差的最大似然估計量為

(16)

選擇利用MATLAB來實現Kriging模型的建立,并將經過LHS抽樣的應力樣本數據代入到建立好的Kriging模型之中,生成應力預測模型。

3.3 通過增加樣本點優化Kriging模型

當Kriging模型精度沒有滿足要求時,會導致預測模型判斷的誤差[32]。利用LHS抽樣出的樣本經過Kriging模型擬合出的模型可能出現精度不夠的問題,主要原因是LHS抽樣會分布于整個設計域空間,但由于隨著裝配參數的改變,密封圈的應力變化分布并不一定是平均的,所以模型可能會損失局部特征。為此,通過局部自動加點的方式提高預測模型的精度,根據MSE為依據判斷精度,在模型精度最差的位置重新從CTGAN生成的數據中挑選LHS抽樣以外的新數據,并自動加入新數據點,加入數據后重新建模重新評估模型精度,并找到新模型精度最差位置再次添加數據,以此循環迭代,直到模型精度達到一定指標,以此達到模型主動學習的目的,Kriging模型的具體建立流程如圖6所示。

圖6 Kriging預測模型流程圖Fig.6 Kriging prediction process

對LHS抽樣的數據按單元位置逐點輸入到通過MATLAB建模的Kriging模型中,并加入預測模型精度判斷,預測模型判斷標準為MSE小于0.000 5即迭代終止。圖7為單個單元不同裝配參數下的應變模型。

圖7 單個單元的Kriging預測模型Fig.7 Kriging model of single element

為方便觀測,以下分別以軸傾角α、偏心距d、橢圓度e兩個參數為一組與應力值對應的關系做圖,如圖8所示。

(a)Inclination angle-Eccentric distance (b)Eccentric distance-Ellipticity (c)Inclination angle-Ellipticity圖8 各裝配參數對應預測關系Fig.8 Corresponding prediction correlation between each parameter

4 試驗驗證

為驗證所建立的預測模型精準度,搭建試驗樣機進行對接裝配試驗,如圖9所示。

圖9 試驗樣機Fig.9 Experimental prototype

試驗用發動機的密封圈尺寸為φ676.5 mm×φ5.5 mm和φ687.5 mm×φ5.5 mm,壓縮率20%,測量采用內置薄膜柔性傳感器的方式,傳感器置于密封圈與密封槽、密封圈與止口之間,其放置方式如圖10所示,以90°為相位均勻選取四點作為采樣點,測量在裝配對接過程中密封圈的受力情況。

圖10 傳感器安裝位置示意圖Fig.10 Installation locations of sensors

將對接過程從起始到終止分為2個時間點,即對接進程的的50%和100%,并分別進行采樣;進行2次試驗。由于試驗樣機數量有限,采用一個噴管并旋轉不同角度以模擬橢圓度的變化進行試驗,試驗所采用的參數及結果如表2所示。

表2 試驗參數與結果Table 2 Experimental parameters and results

試驗結果與預測結果的對比情況如圖11所示,可以看出,密封圈試驗應力分布與預測模型所預測的應力分布趨勢基本一致,并隨裝配條件參數的變化規律一致;4個采樣點預測值與試驗值的誤差ep均處于8%以下,證明本方法較為準確和合理。

圖11 試驗結果與預測結果對比Fig.11 Comparison between the experimental and predicted results

由試驗數據可看出,試驗得到的應力值稍大于預測模型所預測的應力值,誤差的來源可能主要有幾個方面：

(1)對橢圓度的簡化造成了誤差,實際接口變形為不規則形狀;

(2)傳感器具有一定的厚度,造成試驗誤差,且傳感器精度有限;

(3)仿真將部分接觸設置成了粗糙,而實際在一些極端位姿下,密封圈可能出現了滑移;

(4)試驗溫度與仿真溫度不同造成了微小誤差。

5 預測流程

在裝配對接過程中,由六自由度并聯平臺進行噴管的裝配,并聯平臺帶有力位傳感器,可實時傳輸裝配位姿信號;由三維激光掃描儀獲取接口的點云數據,以經過處理的形貌數據獲取接口的橢圓度。

預測流程如圖12所示,以裝配位姿及橢圓度作為輸入,數據實時傳輸到預測模型中,預測模型可在數秒內求得預測值,并生成該裝配位姿下的云圖,并顯示最大應力位置與數值。

圖12 預測流程Fig.12 Workflow of the prediction model

由于經過試驗測試,雖然試驗值與真實值存在偏差,但預測值仍可能會略小于真實值,為提高預測的魯棒性與可靠性,根據式(17)計算考慮預測誤差ep所需的預警值σwarning,將結果向下取整,選擇當最大應力到達設計許用應力值90%時提出警告,以此達到指導裝配的目的。

(17)

6 結論

本文利用有限元仿真和機器學習相結合的方式證明了此類方法適用于裝配應力等的在線預測,得出的預測模型表明三種裝配條件參數與密封圈應力分布均呈現一定相關性。有別于傳統的神經網絡模型,使用了更適合小樣本的Kriging模型,并使用GAN輔助加點的方式優化預測模型。

本方法取代了原本需耗時2 d以上的人工預裝配,并將預測時長縮短至秒級。經試驗驗證對比,該預測方法精度達到90%以上,分析了預測誤差產生的來源,并針對誤差修正所需預警值,證明該預測模型可實現發動機裝配過程中密封圈失效情況的在線預測。

由于使用簡化模型,且仿真本身的設置也影響仿真結果的準確性,使得目前仿真結果與實際仍存在偏差。本文在如下方面還有待深入研究：

(1)進一步細化模型,減小仿真誤差,并考慮裝配時螺栓擰緊等細節方面對密封圈的影響;

(2)在線預測的同時,以減小裝配應力為目標,優化并聯平臺裝配路徑。