基于子程序二次開發(fā)的雙向曲率板應(yīng)變分布計算方法

蔡向東，趙耀，魏振帥，常利春

華中科技大學(xué) 船舶與海洋工程學(xué)院，湖北 武漢 430074

0 引 言

船舶雙向曲率板成型的本質(zhì)是提供目標(biāo)形狀需要的應(yīng)變成分。在實際成型的過程中需要考慮：由平板到目標(biāo)形狀的應(yīng)變分布；什么樣的成型工具和加載方案可以提供目標(biāo)形狀需要的應(yīng)變。應(yīng)變分布是連接變形形狀和加載方案之間的橋梁，計算機(jī)可以根據(jù)應(yīng)變分布的數(shù)據(jù)生成機(jī)器可以執(zhí)行的加載方案，從而實現(xiàn)雙向曲率板的自動化成型?？梢钥闯鰬?yīng)變分布的計算是船舶雙向曲率板自動化成型過程中的重要內(nèi)容，因此提高應(yīng)變分布的計算精度具有重要意義。

大量科研工作者的研究[1-3]表明，線成型加工產(chǎn)生的應(yīng)變分布與加工路徑之間具有明顯的關(guān)系：面內(nèi)應(yīng)變及面外應(yīng)變的主應(yīng)變方向與加工路徑垂直。據(jù)此提出一種自動化成型方案[1-2,4-6]：首先計算初始平板到目標(biāo)形狀的應(yīng)變分布，然后根據(jù)應(yīng)變分布確定加工路徑，最終根據(jù)設(shè)計路徑上的應(yīng)變確定加工參數(shù)。

在上述自動化成型方案中，首先需要根據(jù)目標(biāo)形狀計算應(yīng)變分布，因此應(yīng)變分布的計算精度直接影響到后續(xù)加工路徑和工藝參數(shù)的準(zhǔn)確性。根據(jù)目標(biāo)形狀計算應(yīng)變分布主要有力學(xué)計算方法和微分幾何方法。Cheng 等[1]指出利用微分幾何方法通常只能得到面內(nèi)應(yīng)變分布，很難得到面外應(yīng)變分布，而基于力學(xué)的計算方法可以同時準(zhǔn)確計算面內(nèi)應(yīng)變和面外應(yīng)變的分布；Liu 等[2]使用接觸算法，將目標(biāo)形狀的板置于2 塊剛性平板之間，設(shè)置無摩擦接觸，調(diào)整2 塊剛性平板之間的距離等于板厚，將目標(biāo)形狀壓平，將應(yīng)變數(shù)據(jù)符號反號，可獲得初始平板到目標(biāo)形狀的應(yīng)變分布，但是接觸分析計算效率低，而且容易不收斂；Ueda 等[4-6]使用彈性大變形有限元方法計算初始平板到目標(biāo)形狀的應(yīng)變分布，但是具體的計算過程沒有提及；胡昌成等[7]通過直接施加位移場的方式將位移載荷施加到平板上計算位移場，從而獲得目標(biāo)形狀的應(yīng)變分布，該方法只在目標(biāo)形狀曲率較小時適用，在目標(biāo)曲率較大時直接施加位移場計算的應(yīng)變分布不準(zhǔn)確；Yu 等[8]提出了一種基于非線性優(yōu)化的雙曲面展開算法，通過求解一個有約束的非線性規(guī)劃問題得到最小面內(nèi)應(yīng)變的分布，但是該算法無法計算面外應(yīng)變的分布；Shin 等[9]提出了板材變形的運動學(xué)分析方法：面內(nèi)應(yīng)變由初始板和目標(biāo)曲面的映射展開關(guān)系決定，面外應(yīng)變直接來源于目標(biāo)曲面的主曲率。由于復(fù)雜曲面的展開及其主曲率計算困難，因而該方法實際運用過程中計算量大，計算效率較低。

本文將討論直接施加位移場和模具接觸分析這2 種方法計算應(yīng)變分布的局限性，將從以下幾個方面對應(yīng)變分布的計算方法進(jìn)行研究：使用ABAQUS 子程序DISP 進(jìn)行二次開發(fā)，改進(jìn)了位移場的加載方式，提高了位移場計算精度，從而能準(zhǔn)確獲取應(yīng)變分布；考慮到實際加工過程是彈塑性的變形過程，因而討論了材料非線性對應(yīng)變分布的影響；考慮到計算效率和計算精度之間的矛盾關(guān)系，對幾種計算應(yīng)變分布方法的適用性進(jìn)行討論，根據(jù)目標(biāo)形狀的曲率半徑選取合適的計算方法來提高計算效率，同時滿足計算精度的要求；最后通過實例驗證子程序方法對于計算船舶典型雙向曲率板應(yīng)變的準(zhǔn)確性。

1 基本理論

1.1 變形與應(yīng)變的關(guān)系

船舶雙向曲率板成型的本質(zhì)是將塑性變形施加到初始平板上來獲得目標(biāo)形狀，在由初始板變?yōu)槟繕?biāo)板的過程中，應(yīng)變是常用控制要素，可以精確地表征物體的相對變形，因此獲得應(yīng)變與變形之間的關(guān)系非常重要；船舶雙向曲率板成型屬于薄板大撓度問題，基于薄板的大變形理論，其變形與應(yīng)變的關(guān)系為

式中：x，y，z分別為板的長度方向、寬度方向、厚度方向；u，v，w分 別為x，y，z方 向的位移； ε和 γ及其下標(biāo)分別表示某一個方向的軸向應(yīng)變和剪應(yīng)變。

1.2 面內(nèi)應(yīng)變與面外應(yīng)變

通過有限元計算可以獲取初始平板到目標(biāo)形狀的應(yīng)變分布，計算得到的應(yīng)變可以分解為面內(nèi)應(yīng)變和面外應(yīng)變。面內(nèi)應(yīng)變是沿著板厚方向均勻分布的應(yīng)變，面外應(yīng)變是沿著厚度方向線性變化的應(yīng)變[10]，通過式（2）和式（3）可以計算各個方向的面內(nèi)應(yīng)變成分和面外應(yīng)變成分。

1.3 主應(yīng)變的計算

目前主要是根據(jù)面內(nèi)主應(yīng)變和面外主應(yīng)變的大小和方向來確定加工路徑和工藝參數(shù)，因此需要對主應(yīng)變進(jìn)行計算?？梢愿鶕?jù)式（2）和式（3）計算得到的面內(nèi)應(yīng)變成分和面外應(yīng)變成分，按照式（4）計算得到面內(nèi)應(yīng)變或面外應(yīng)變的第一主應(yīng)變和第二主應(yīng)變。

由式（1）可知，應(yīng)變是幾何性質(zhì)的表現(xiàn)，僅僅涉及變形的大小而與材料的力學(xué)參數(shù)無關(guān)，在計算曲面板的應(yīng)變分布時可通過式（1）計算獲得，然而實際成型過程中，板的曲面表達(dá)為離散的點，基于數(shù)學(xué)公式的推導(dǎo)過程復(fù)雜，因此本文采用大變形有限元對計算應(yīng)變分布的方法進(jìn)行研究。

2 應(yīng)變分布計算方法及比較分析

2.1 局限性分析

根據(jù)目標(biāo)形狀計算應(yīng)變分布有力學(xué)計算方法和微分幾何方法。微分幾何方法通常只能計算得到面內(nèi)應(yīng)變，而力學(xué)計算方法能同時得到面內(nèi)應(yīng)變和面外應(yīng)變的分布，在計算應(yīng)變分布方面具有一定的優(yōu)勢，本文只針對力學(xué)計算應(yīng)變分布的方法進(jìn)行研究。

目前力學(xué)獲取應(yīng)變的方法有直接施加位移場進(jìn)行力學(xué)計算和使用剛體模具進(jìn)行接觸計算這2 種方法。從式(1)可以看出：應(yīng)變分布與形狀具有一一對應(yīng)的關(guān)系，應(yīng)變分布計算的準(zhǔn)確性取決于位移場計算結(jié)果的準(zhǔn)確性；為了方便對計算方法進(jìn)行描述，將上述2 種方法分別命名為直接施加位移場和模具施加位移場的計算方法。圖1 為直接施加位移場計算應(yīng)變分布的示意圖，RP為參考點；圖2 為模具施加位移場計算應(yīng)變分布的示意圖。直接施加位移場是根據(jù)目標(biāo)曲面散點的坐標(biāo)擬合目標(biāo)曲面，根據(jù)目標(biāo)曲面計算每個節(jié)點z方向的位移，將z方向的位移作為邊界條件進(jìn)行計算；模具施加位移場是根據(jù)目標(biāo)形狀建立剛體模具，將平板置于2 個剛體模具中間，調(diào)整模具間的距離等于板厚，使用接觸算法進(jìn)行計算。

圖1 直接施加位移場方法Fig. 1 Direct displacement field application

圖2 模具施加位移場方法Fig. 2 Mould displacement field application

圖3 所示為目標(biāo)曲面的曲率半徑R分別為1.5，3.0 和4.5 m 時，2 種計算方法得到的曲面中心截面處的位移場計算結(jié)果。可以看出：曲率小時，直接施加位移場計算得到的形狀與目標(biāo)方程很接近，應(yīng)變的計算結(jié)果滿足要求；曲率大時，計算得到的形狀與目標(biāo)形狀的方程出現(xiàn)較大的偏差。產(chǎn)生這種現(xiàn)象的原因是：z方向位移場是根據(jù)初始節(jié)點的坐標(biāo)計算得到，大變形時節(jié)點的x，y位置發(fā)生較大偏移，而計算過程中每個節(jié)點的z值不會發(fā)生改變，導(dǎo)致計算的曲面與目標(biāo)曲面發(fā)生偏離；模具施加位移場的方法在計算精度上可以滿足要求，但在建模的復(fù)雜程度和計算效率上存在明顯的不足。

圖3 位移場計算結(jié)果Fig. 3 Calculation results of displacement field

2.2 子程序施加位移場的計算方法

由于目前準(zhǔn)確計算位移場的方法在計算精度和計算效率方面還存在局限性，使用ABAQUS 子程序DISP 進(jìn)行二次開發(fā)，改進(jìn)位移場的加載方式，可以準(zhǔn)確獲取目標(biāo)形狀的位移場，從而獲取準(zhǔn)確的應(yīng)變分布。為了方便描述，將該方法命名為子程序施加位移場的計算方法。子程序施加位移場的計算方法是根據(jù)目標(biāo)形狀擬合曲面，逐步施加位移場載荷，每個載荷步計算結(jié)束后獲取節(jié)點當(dāng)前坐標(biāo)，根據(jù)當(dāng)前坐標(biāo)重新計算位移載荷，可以使計算的形狀接近目標(biāo)形狀，使計算的應(yīng)變結(jié)果滿足要求。圖4 所示為子程序施加位移場的流程，子程序的主要作用是更新當(dāng)前節(jié)點坐標(biāo)、根據(jù)當(dāng)前節(jié)點坐標(biāo)計算節(jié)點位移場。圖中：h為總步長； Δh為計算允許的最大步長，用來控制計算精度； Δh′為有限元自動生成步長；k為累計步長。

圖4 子程序計算位移場流程Fig. 4 The flow of calculating displacement field using subroutine

2.3 比較分析

為了比較這3 種計算方法的計算精度和計算效率，選取半徑為r=1 m，板厚t=0.02 m 的圓形平板施加位移場，由于目標(biāo)曲面曲率半徑R=1.5 m時位移場的結(jié)果存在明顯差異，因此選取目標(biāo)曲面曲率半徑R=1.5 m 對子程序計算方法進(jìn)行驗證。使用3 種方法進(jìn)行彈性大變形計算，計算結(jié)果如圖5 所示。從計算結(jié)果看，直接施加位移場的方法相對于目標(biāo)形狀會出現(xiàn)明顯的偏差；模具施加位移場與子程序施加位移場都能很好地接近目標(biāo)形狀；同時比較了3 種計算方法的效率，表1所示為3 種方法在相同計算條件下的計算時長，子程序施加位移場的計算效率只是略低于直接施加位移場的計算效率，遠(yuǎn)高于模具施加位移場的計算效率?？紤]建模的復(fù)雜性以及計算的效率，子程序施加位移場的方法具有明顯優(yōu)勢。

圖5 位移場計算結(jié)果Fig. 5 Calculation results of displacement field

表1 直接施加位移場、子程序施加位移場和模具施加位移場三種方法計算時長的比較Table 1 Comparison of the calculation times for the three methods: direct displacement field application, subroutine displacement field application and mould displacement field application

3 材料非線性對應(yīng)變分布的影響

船舶雙向曲率板實際成型的過程是彈塑性的過程，在討論獲取應(yīng)變計算方法的同時應(yīng)該討論材料非線性對應(yīng)變分布的影響。通過前面的計算，基于子程序方法具有計算簡單且滿足計算精度要求的優(yōu)點，因此選取該方法來討論材料非線性對應(yīng)變分布的影響。

選取模型尺寸為半徑r=1 m，板厚t=0.02 m 的圓形平板，目標(biāo)曲面為曲率半徑R=2 m 的球型殼，材料為Q235 鋼，使用各向同性硬化模型，其彈性階段以及塑性階段的參數(shù)見表2[11]以及表3[12]。在該模型上分別進(jìn)行彈性大變形計算和彈塑性大變形計算，比較兩者應(yīng)變分布場的差異。

表2 Q235 鋼材料參數(shù)Table 2 Q235 steel material parameters

表3 Q235 鋼塑性應(yīng)力-塑性應(yīng)變數(shù)據(jù)Table 3 Data of Q235 steel plastic stress-plastic strain

根據(jù)式（3）將應(yīng)變數(shù)據(jù)分解為面內(nèi)應(yīng)變與面外應(yīng)變，按照式（4）計算主應(yīng)變的大小和方向，繪制其主應(yīng)變的矢量圖，由于模型的對稱性，只繪制了1/4 模型的結(jié)果。圖6 和圖7 所示分別為彈性和彈塑性計算的面內(nèi)應(yīng)變分布，從圖中可以看出二者應(yīng)變分布的形式與分布規(guī)律相似。由于模型的對稱性，只提取模型對稱線位置上的面內(nèi)主應(yīng)變的值進(jìn)行對比，如圖8 與圖9 所示。彈性和彈塑性計算的周向面內(nèi)應(yīng)變分布在數(shù)值上接近，彈性和彈塑性計算的徑向面內(nèi)應(yīng)變的分布有一定的差異。由于材料的非線性，導(dǎo)致2 種計算方法的面內(nèi)位移分布不一樣。圖10 所示為面內(nèi)位移分布，由式（1）可以得知，徑向面內(nèi)位移分布不一樣會導(dǎo)致應(yīng)變分布的不一樣。如圖9 所示，面外應(yīng)變的值均位于0.004 9～0.005 1，接近面外應(yīng)變理論值t/2R=0.005。

圖6 彈性大變形計算面內(nèi)應(yīng)變分布Fig. 6 In-plane strain distribution of elastic calculation

圖7 彈塑性大變形計算面內(nèi)應(yīng)變分布Fig. 7 In-plane strain distribution of elastoplastic calculation

圖8 面內(nèi)應(yīng)變分布曲線Fig. 8 Distribution of in-plane strain

圖9 面外應(yīng)變分布曲線Fig. 9 Distribution of out-plane strain

圖10 面內(nèi)位移分布曲線Fig. 10 Distribution of in-plane displacement

彈性大變形計算與彈塑性大變形計算的面內(nèi)應(yīng)變分布結(jié)果有一定的差異，根本原因是這2 種計算方法得到的變形位移場存在細(xì)微的區(qū)別，而這2 種計算方法得到的位移場均接近目標(biāo)形狀，滿足船舶制造精度要求，因此它們得到的應(yīng)變分布圖均能作為指導(dǎo)加工的依據(jù)。而實際加工很難做到一次成型，會涉及多次應(yīng)變的計算，考慮到計算的效率，建議使用彈性大變形的計算方法來計算應(yīng)變。

4 適用性分析

船舶雙向曲率板的實際成型過程中，很難做到一次成型，實際成型過程中會多次計算應(yīng)變的分布，計算的效率會直接影響實際加工的總時長，因此對于應(yīng)變計算方法的選取要處理計算精度和計算效率之間的矛盾，在保證計算精度的同時，盡量選取計算效率高的方法。

前文討論了3 種方法的計算精度，可以看出模具施加位移場和子程序施加位移場的計算精度相當(dāng)，直接施加位移場的方法在曲率較大時會出現(xiàn)較大的誤差。從建模復(fù)雜程度以及計算效率的角度來看，直接施加位移場的方法建模最簡單、計算效率最高；子程序施加位移場的方法建模簡單，計算效率高；模具施加位移場的方法建模最復(fù)雜，計算效率遠(yuǎn)低于前2 種計算方法。實際加工過程中可能會涉及多次應(yīng)變分布的計算，因此計算效率是一個重要的因素，會影響實際加工的總時長，因此不建議使用模具施加位移場的方法。

直接施加位移場和子程序施加位移場在目標(biāo)曲面曲率較小時計算的位移場結(jié)果接近，為了提高計算效率，可以確定一個曲率半徑的臨界值，曲率半徑小于臨界值時使用子程序的計算方法，曲率半徑大于臨界值時使用直接施加位移場的計算方法。選取半徑r=1 m，板厚t=0.02 m 的圓形平板為研究對象，對于不同曲率半徑的目標(biāo)曲面進(jìn)行計算，設(shè)置目標(biāo)曲面的曲率半徑分別為2，4，6，8 和10 m，該系列曲率半徑幾乎涵蓋90%以上的船舶雙向曲率板。選取面內(nèi)應(yīng)變與面外應(yīng)變的最大值作為目標(biāo)，如圖11 和圖12 所示，當(dāng)曲率半徑大于6 m 時，面內(nèi)應(yīng)變的最大值與面外應(yīng)變的最大值誤差在5%以內(nèi)，這意味著曲率半徑大于6 m時，可以使用直接施加位移場方法，從而提高計算效率。

圖11 曲率半徑與面內(nèi)應(yīng)變Fig. 11 Radius of curvature and in-plane strain

圖12 曲率半徑與面外應(yīng)變Fig. 12 Radius of curvature and out-plane strain

5 實例驗證

上文的主要研究對象為圓形平板，為了證明計算方法的有效性，選取矩形平板為研究對象，對理想形狀的馬鞍型和風(fēng)帆型進(jìn)行了計算，矩形板尺寸為2 m×1 m×0.02 m，楊氏模量為2.1×1011Pa，泊松比為0.3。馬鞍型的目標(biāo)形狀方程為z=-0.1x2+0.5y2，風(fēng)帆型的目標(biāo)形狀方程為z=0.1x2+0.5y2， 2 個 目 標(biāo) 形 狀y方 向 和x方 向 的 曲率半徑分別為1 和5 m，該曲率半徑滿足船舶上大曲率板的要求。選用模具施加位移場和子程序施加位移場2 種方法進(jìn)行彈性大變形計算，以模具施加位移場的計算結(jié)果為標(biāo)準(zhǔn)，驗證了子程序施加位移場的準(zhǔn)確性。圖13 所示為使用模具施加位移場方法計算得到的馬鞍型板的面內(nèi)和面外主應(yīng)變分布，圖14 所示為使用子程序方法計算得到的面內(nèi)和面外主應(yīng)變分布，二者的分布規(guī)律一致。圖13 和圖14 中標(biāo)注了2 種計算方法的主應(yīng)變最大值 |max|，經(jīng)過對比可知，2 種計算方法的面內(nèi)主應(yīng)變相對誤差在4%以內(nèi)，面外主應(yīng)變的相對誤差在1%以內(nèi)。圖15 所示為使用模具施加位移場方法計算得到的風(fēng)帆型板材的面內(nèi)和面外主應(yīng)變分布，圖16所示為使用子程序方法計算得到的面內(nèi)和面外主應(yīng)變，二者的分布規(guī)律一致。同樣，根據(jù)圖15 和圖16 所標(biāo)注的主應(yīng)變最大值可知，2 種計算方法的面內(nèi)主應(yīng)變的相對誤差在3%以內(nèi)，面外主應(yīng)變的相對誤差在2%以內(nèi)。這2 種典型形狀的計算結(jié)果進(jìn)一步說明了子程序施加位移場計算方法的準(zhǔn)確性。其中最大誤差的位置在板材的中間區(qū)域。產(chǎn)生誤差的主要原因在于：子程序計算方法是通過逐步施加位移場逼近精確的位移場，計算結(jié)果的精度與選取的計算步長有關(guān)系。

圖13 模具施加位移場方法計算的主應(yīng)變分布Fig. 13 Principal strain distribution calculated by mould displacement field application

圖14 子程序施加位移場方法計算的主應(yīng)變分布Fig. 14 Principal strain distribution calculated by subroutine displacement field application

圖15 模具施加位移場方法計算的主應(yīng)變分布Fig. 15 Principal strain distribution calculated by mould displacement field application

圖16 子程序施加位移場方法計算的主應(yīng)變分布Fig. 16 Principal strain distribution calculated by subroutine displacement field application

6 結(jié) 論

本文指出了目前力學(xué)獲取應(yīng)變方法的局限性，使用ABAQUS 子程序進(jìn)行二次開發(fā)，改進(jìn)了位移場施加的方法，提高了應(yīng)變分布的計算精度，并討論了材料非線性對應(yīng)變分布的影響，對于各種計算方法的適用性進(jìn)行了分析。通過本文研究得出了以下結(jié)論：

1） 目標(biāo)形狀的曲率較大時，使用直接施加位移場的方法會導(dǎo)致位移場的計算結(jié)果不準(zhǔn)確，使用模具施加位移場的方法建模復(fù)雜，計算時間長，子程序施加位移場的方法能很好解決建模復(fù)雜和計算時間長的問題。

2） 材料非線性會導(dǎo)致彈性大變形與彈塑性大變形計算的位移場結(jié)果出現(xiàn)細(xì)微的差別，導(dǎo)致應(yīng)變分布出現(xiàn)差異，但應(yīng)變分布規(guī)律基本一致，考慮到計算的效率，建議使用彈性大變形的計算方法來計算應(yīng)變。

3） 目標(biāo)形狀的曲率較小時，直接施加位移場與子程序施加位移場計算的位移場結(jié)果差異很小，可以使用直接施加位移場的方法來獲取應(yīng)變，提高計算效率。