基于變形力數(shù)據(jù)的結(jié)構(gòu)件殘余應力場重構(gòu)方法

黃 沖，劉長青，趙智偉，劉明芳，郭立杰

（1.南京航空航天大學直升機傳動技術(shù)國家級重點實驗室，江蘇南京 210016；2.上海航天設備制造總廠有限公司，上海 201109）

0 引言

在毛坯的制造、零件的加工和使用過程中，殘余應力是不可避免的，由于熱、力等效應，它會影響材料的變形穩(wěn)定性、加工精度、抗疲勞性能、脆性斷裂、應力腐蝕開裂和硬度［1-5］。以飛機結(jié)構(gòu)件的加工為例，由于其零件尺寸大，加工中材料去除率高（最高可達95%以上［6-7］），在完成加工后，殘余應力的釋放會導致飛機結(jié)構(gòu)件產(chǎn)生彎曲、扭曲等各種形式加工變形［8］，嚴重者會導致零件報廢。若能了解殘余應力的分布，即零件內(nèi)部的殘余應力場，有助于消除和控制機械零件在設計、制造、評價和使用過程中的殘余應力，所以殘余應力場的重建具有很高的實用價值。

隨著殘余應力測試技術(shù)與計算機計算能力的提高，人們對殘余應力的分析能力也日愈加強。國內(nèi)外的學者做了大量的相關(guān)研究，主要可以分為3 個方面：實驗法、數(shù)學模型法、有限元法。

實驗測試作為驗證理論分析的有效方法，一直以來對于通過使用不同的測試手段來分析殘余應力的探索從未間斷。殘余應力的實驗法包括破壞性實驗法和非破壞性實驗法。破壞性實驗主要有鉆孔法［9］、逐層剖離法［10］、納米壓痕法［11］等，這類方法的機理就是通過釋放存在材料內(nèi)部的原始應力，進而通過測量該過程所產(chǎn)生的應變位移反求出殘余應力值。該方法的問題是殘余應力測量后零件的結(jié)構(gòu)也會遭到破壞，無法進一步使用。非破壞性實驗主要有X 射線衍射法［12］、超聲法［13］、中子衍射法［14］等無損檢測方法，這類方法目前只能檢測表面或近表面的殘余應力。

在使用數(shù)學模型來構(gòu)建殘余應力方面，JACOBUS 等［15］通過建立平面應變熱彈塑性模型，研究材料的切削運動以預測工件表面和亞表面的殘余應力分布。EL-AXIR［16］建立了一個考慮材料彈性強度、切削速度和進給量的模型用來分析以上工藝參數(shù)對于最大殘余應力的影響，最終確定在加工過程中最大殘余應力的位置和深度。LAZOGLU 等［17］采用松弛法建立了一個增強的解析彈塑性模型，并用X 射線衍射測量得到了驗證，有效提高了殘余應力預測效率和精度。國內(nèi)對使用數(shù)學模型構(gòu)建殘余應力的研究較少，其中華中科技大學的毛寬民等［18］提出了一種利用表面應力測量數(shù)據(jù)預測零件完整的殘余應力狀態(tài)的方法，求解了單元和整體的殘余應力平衡方程，得到了結(jié)構(gòu)的殘余應力場。但是數(shù)學模型法受其預測精度和求解難度的限制，仍需要開展更多的相關(guān)的研究工作。

為了分析零件殘余應力的狀態(tài)，許多學者探索了基于有限元的求解方法，在切屑-刀具接觸模型、切屑分離準則和網(wǎng)格重劃分技術(shù)等方面都取得了很大的進展。LIN 等［19］通過在有限元仿真環(huán)境中研究了切深、切削速度等工藝參數(shù)對殘余應力的生成的影響。YANG 等［20］采用“生死單元”技術(shù)模擬材料去除，通過有限元仿真分析了殘余應力對飛機結(jié)構(gòu)件加工變形的影響。GUO 等［21］利用有限元法分析了在車削和磨削中殘余應力的分布規(guī)律。但是有限元法局限于計算耗時、模型精度低，所以該方法依然需進一步提高。

從以上分析可以知道，殘余應力測量方法大多比較復雜昂貴，測量范圍和精度受限，并且數(shù)學模型與有限元法預測精度有待提高。針對這些問題，本文提出了一種通過零件加工過程中的變形力數(shù)據(jù)實現(xiàn)零件殘余應力場重構(gòu)的方法，該方法首先提取了由于材料去除引起的應力重新分布帶來的裝夾力變化數(shù)據(jù)，即變形力，殘余應力到變形力之間通過協(xié)方差矩陣自適應進化策略（CMA-ES 算法）反復迭代求解出映射關(guān)系，從而構(gòu)建出零件的殘余應力場，較好地實現(xiàn)了預測殘余應力場。

1 基于CMA-ES算法殘余應力場重構(gòu)

因為殘余應力場受到多種因素的耦合作用，很難精準地建立一個物理模型來實現(xiàn)殘余應力場的重構(gòu)，所以我們把它處理成一個黑盒問題，而CMAES 算法在這種問題上能展現(xiàn)出非常好的實用價值，因此，我們在殘余應力到變形力之間采用此算法來求解出映射關(guān)系，從而構(gòu)建出殘余應力場。協(xié)方差矩陣自適應進化策略（Covariance Matrix Adaptation Evolutionary Strategies，CMA-ES）［22-23］，該算法使用高斯分布在優(yōu)化問題的解空間中進行采樣，并且根據(jù)某種樣本選擇機制對高斯分布進行采樣和更新過程持續(xù)迭代，當搜索到滿意的解或者達到采樣最大次數(shù)等停止條件時就會停止優(yōu)化過程。

故針對優(yōu)化目標：

式中：f(x)為優(yōu)化目標函數(shù)；x為輸入樣本；S為D維的實數(shù)空間。

CMA-ES 使用高斯分布N(mt，σ2t Ct)在D維的S空間進行采樣，其中mt為第t代的高斯分布的均值，σt為第t代的搜索步長，是用于控制高斯分布的局部搜索能力，Ct為第t代的高斯分布的協(xié)方差矩陣。在CMA-ES 中，每一次迭代從N(mt，σ2t Ct)中產(chǎn)生λ個樣本（一般λ取值為logn的量級，n為解空間的維度），然后對樣本進行適應度計算，也就是計算f(x)，然后使用選擇策略從樣本中選出適應度最好的前樣本用于更新分布參數(shù)。CMA-ES算法是一個不斷迭代的優(yōu)化算法，在采樣操作和高斯分布更新兩個階段中不斷循環(huán)，當滿足停止條件時就會停止計算。

使用所選擇的樣本分別對分布參數(shù)mt、σt、Ct進行獨立的更新。大體來說，對mt、Ct的更新使用的是最大似然估計（ML-update）。分布的均值即為所選擇的μ個樣本的加權(quán)的最大似然估計：

式中：mt+1為第t+1 代均值；wi為第t代每個樣本的權(quán)重系數(shù)，這里權(quán)重系數(shù)取相同的，所以此式就是使用所選擇的樣本的最大似然估計；xi：λ為第t代λ個樣本中適應度最好的前μ個樣本。更新項可以理解為對m的自然梯度信息幾何優(yōu)化、隨機優(yōu)化與進化策略。

CMA-ES 默認使用累積式步長調(diào)整（Cumulative step size adaptation，CSA）。CSA 是當前最成功、用得最多的步長調(diào)整方式，其原理是相繼搜索的方向應該是共軛的。如果相繼搜索方向之間正相關(guān)，則表明步長太小，應該增大；如果相繼搜索方向之間負相關(guān)，則表明步長太大，應該減小。所以步長的更新公式為

式 中：σt+1為 第t+1 代搜索步長；cσ與dσ為調(diào)整步長變化幅度的控制參數(shù)，通常設置為，dσ＞1；st+1為第t+1 代共軛進化路徑，其搜索路徑可以看成是一個n維標準正態(tài)分布的隨機向量，在平穩(wěn)性條件下有st+1～N（0，I）。因此，其模長服從卡方分 布，而。

高斯分布協(xié)方差矩陣C的更新結(jié)合了rank-1 和rank-μ，更新公式為

式中：Ct+1為第t+1 代的高斯分布協(xié)方差矩陣；pt+1為 第t+1 代的平滑指 數(shù)；c1與cμ為學習率；yi：λ為 第t代λ個樣本搜索方向中適應度最好的前μ個。C的更新原理是增大沿成功搜索方向的方差，即增大沿這些方向采樣的概率。rank-1 使用平滑指數(shù)p對均值的偏差進行記憶，隨著優(yōu)化過程歷史信息逐代衰減；rank-μ為對C的自然梯度信息的幾何優(yōu)化的策略。這里的每一個yi∈RD是一個搜索方向，一般情況下y可以通過協(xié)方差矩陣C的特征分解或者Cholesky 分解得到。學習率c1與cμ的設計原理是，即學習率與所調(diào)整的變量維度（參數(shù)個數(shù)）成反比。

根據(jù)變量維度從CMA-EA 算法的解空間賦予材料各層的殘余應力值，在各層材料去除的過程中提取由于殘余應力重分布而產(chǎn)生的變形力。將變形力數(shù)據(jù)與原始變形力差的平方和作為此算法的優(yōu)化目標，在優(yōu)化目標最小化的過程中，從解空間中得到的應力值就會無限逼近初始施加的應力值，從而實現(xiàn)殘余應力的重構(gòu)。

2 基于Abaqus 二次開發(fā)仿真環(huán)境的建立與變形力數(shù)據(jù)獲取

在零件制造過程中，變形的產(chǎn)生大部分為彈性變形，本文在不考慮零件塑性變形的情況下，根據(jù)線性彈性問題的解是唯一的，即基于基爾霍夫唯一性定理［24］可以知道應力場的解是唯一的，也就是說，當一個毛坯的內(nèi)部應力場確定，按照同樣的方式去除材料在同一點得到的變形力是唯一的。所以本文可以通過在二維仿真環(huán)境中獲取材料去除過程中的變形力，從而構(gòu)建出零件的殘余應力場。以下分為仿真環(huán)境的搭建與二次開發(fā)獲取變形力數(shù)據(jù)來論述。

2.1 仿真環(huán)境的搭建

通過在Python 腳本在Abaqus 軟件在平臺上建立一個70 mm×10 mm 二維平面，本文以廣泛應用于飛機結(jié)構(gòu)件中的7050-T7451 鋁合金為例，材料的力學屬性見表1。在仿真環(huán)境中為所建平面賦予該鋁合金的材料屬性。邊界條件則是限定平面中間區(qū)域2 mm×0.1 mm 的6 個自由度，在加工仿真中將其作為固定區(qū)域。網(wǎng)格劃分則是以大小為1 mm2四邊形在平面上進行劃分，以保證材料去除的完整性，如圖1 所示。

根據(jù)黃曉明［25］提出7050-T7451 鋁合金坯料的初始殘余應力場函數(shù)，見式（1），在平面上沿Y軸方向上逐層施加沿X方向的主應力σx，如圖1 所示。

表1 7050-T7451 鋁合金力學屬性Tab.1 Mechanical properties of 7050-T7451 aluminum alloy

式中：Z為毛坯上某一位置到毛坯中性層的相應距離。

圖1 基于Abaqus 平臺仿真環(huán)境的建立Fig.1 Establishment of the simulation environment based on Abaqus platform

施加方法是首先將初始殘余應力場離散化，然后將離散后的應力值逐層施加到網(wǎng)格單元中，將初始殘余應力與平面的寬度（毛坯厚度）之間的函數(shù)關(guān)系通過Python 腳本施加到鋁合金預拉伸毛坯板材的二維平面中，這樣保證了初始殘余應力的連續(xù)性，提高了計算精度，同時實現(xiàn)了初始殘余應力施加的自動化。

2.2 變形力數(shù)據(jù)的獲取

本文通過模擬實際的加工過程，提取每去除一層材料所產(chǎn)生的變形力。實際加工時材料的去除是一個動態(tài)的過程，刀具走過的位置材料就會被去除，在Abaqus 仿真平臺采用“生死單元”技術(shù)來實現(xiàn)材料的去除過程。通過Python 腳本在仿真平臺上將所去除的單元的剛度矩陣乘以一個約為零的縮減因數(shù)，來達到殺死單元的目的，具體的表現(xiàn)是將與單元相關(guān)的單元屬性，比如質(zhì)量、阻尼、比熱容等全部置為零。在結(jié)構(gòu)件有限元變形計算中，去除材料單元一旦被殺死其應變即為零，從而避免了從材料變形的微觀角度來討論被切除單元的分離標準。

在加工仿真中，將材料分為9 層去除，每層去除厚度為1 mm 的兩個槽，隨著模擬加工過程中單元的去除，殘余應力逐層釋放，每去掉一層單元，在有限元中分為一個分析步計算一次，一共有9 個分析步。通過Python 腳本在平面底邊兩個頂點上施加彈簧用來提取每層材料去除產(chǎn)生的變形力，如圖2所示。

圖2 變形力數(shù)據(jù)提取Fig.2 Acquisition of deformation force data

使用Python 腳本通過Abaqus 接口訪問Abaqus對象中數(shù)據(jù)。在對象創(chuàng)建后，可使用該對象提供的方法來處理對象中的數(shù)據(jù)成員。ODB 對象是結(jié)果數(shù)據(jù)庫對象，包含了模型數(shù)據(jù)和結(jié)果數(shù)據(jù)，本文通過訪問ODB 下面的steps 對象，遍歷多個分析步和幀的位移值（U），部件實例的節(jié)點集合NodeSets，讀取變形力提取點1 和2 的節(jié)點在每一層材料去除后的位移值，再與彈簧彈性系數(shù)相乘得到變形力數(shù)據(jù)，寫入CSV 文件以作CMA-ES 算法的數(shù)據(jù)輸入。

3 重構(gòu)結(jié)果對比與分析

根據(jù)初始應力場施加的層數(shù)確定CMA-EA 算法解空間的維度，本文案例的維度為10。根據(jù)變量維度從算法的解空間賦予材料各層的殘余應力值，本文為解空間賦予的初始值為0 MPa。在算法的解空間不斷搜索的過程中使得優(yōu)化目標最小化，從而在優(yōu)化目標最小化時得到的應力值就會逼近初始施加的應力值。每層X方向主應力與訓練次數(shù)關(guān)系如圖3 所示。從圖中可以看出，在仿真次數(shù)為1 400次左右以后，基于CMA-ES 算法的殘余應力場解空間趨于收斂，各層應力值只會小幅度波動。因此就有理由認為在這1 810 次仿真中，我們優(yōu)化的目標值即每層變形力差值的平方和最小的狀態(tài)，就是各層應力最逼近原始應力的狀態(tài)。

本文方法的優(yōu)化目標與仿真次數(shù)的關(guān)系如圖4所示。從圖中可以看出，仿真次數(shù)在1 400 次左右以后優(yōu)化目標趨于收斂，變化趨勢趨于平緩，不會再有大幅度波動。因此，有理由認為基于CMA-ES算法的殘余應力場解空間中使得優(yōu)化目標最小的值包含在這1 810 次仿真中。同時，仿真次數(shù)在356次的時取到了優(yōu)化目標的最小值0.012。

本文根據(jù)式（1）施加的X方向主應力的原始值與優(yōu)化目標最小時所施加的各層應力值的對比見表2。從表中可以看出，這10 層的應力重構(gòu)值與與初始施加的原始值變化趨勢是相同的，其中最大誤差為0.8 MPa，最小誤差為0.02 MPa，平均誤差為0.38 MPa，如圖5 所示。

圖3 各層應力重構(gòu)值與仿真次數(shù)關(guān)系Fig.3 Relationship between the reconstruction value of each layer stress and the simulation time

表2 應力場重構(gòu)結(jié)果與施加的初始殘余應力對比Tab.2 Comparison of the stress field reconstruction result and the applied initial residual stress

圖4 優(yōu)化目標與仿真次數(shù)關(guān)系Fig.4 Relationship between the optimization objective and the simulation time

4 結(jié)束語

本文研究了一種通過零件加工過程中的變形力數(shù)據(jù)實現(xiàn)零件殘余應力場重構(gòu)的方法。通過少量點的變形力數(shù)據(jù)就能重構(gòu)出零件的殘余應力場，結(jié)果表明該方法能較好地預測殘余應力場，重構(gòu)出的零件的初始殘余應力場與原始施加的應力場平均誤差值約0.38 MPa。但本文只考慮二維平面的應力場，后續(xù)將研究推廣到三維實體，并將通過提取更多點的變形力數(shù)據(jù)，從而提高應力場重構(gòu)的精度。

圖5 應力值與重構(gòu)誤差值對比Fig.5 Comparison of the stress value and the reconstruction error