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

黃 沖，劉長(zhǎng)青，趙智偉，劉明芳，郭立杰

（1.南京航空航天大學(xué)直升機(jī)傳動(dòng)技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室，江蘇南京 210016；2.上海航天設(shè)備制造總廠有限公司，上海 201109）

0 引言

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

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

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

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

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

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

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

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

故針對(duì)優(yōu)化目標(biāo)：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表2 應(yīng)力場(chǎng)重構(gòu)結(jié)果與施加的初始?xì)堄鄳?yīng)力對(duì)比Tab.2 Comparison of the stress field reconstruction result and the applied initial residual stress

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

4 結(jié)束語(yǔ)

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

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