基于子程序二次開發的雙向曲率板應變分布計算方法

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

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

0 引 言

船舶雙向曲率板成型的本質是提供目標形狀需要的應變成分。在實際成型的過程中需要考慮：由平板到目標形狀的應變分布；什么樣的成型工具和加載方案可以提供目標形狀需要的應變。應變分布是連接變形形狀和加載方案之間的橋梁，計算機可以根據應變分布的數據生成機器可以執行的加載方案，從而實現雙向曲率板的自動化成型?？梢钥闯鰬兎植嫉挠嬎闶谴半p向曲率板自動化成型過程中的重要內容，因此提高應變分布的計算精度具有重要意義。

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

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

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

1 基本理論

1.1 變形與應變的關系

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

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

1.2 面內應變與面外應變

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

1.3 主應變的計算

目前主要是根據面內主應變和面外主應變的大小和方向來確定加工路徑和工藝參數，因此需要對主應變進行計算?？梢愿鶕剑?）和式（3）計算得到的面內應變成分和面外應變成分，按照式（4）計算得到面內應變或面外應變的第一主應變和第二主應變。

由式（1）可知，應變是幾何性質的表現，僅僅涉及變形的大小而與材料的力學參數無關，在計算曲面板的應變分布時可通過式（1）計算獲得，然而實際成型過程中，板的曲面表達為離散的點，基于數學公式的推導過程復雜，因此本文采用大變形有限元對計算應變分布的方法進行研究。

2 應變分布計算方法及比較分析

2.1 局限性分析

根據目標形狀計算應變分布有力學計算方法和微分幾何方法。微分幾何方法通常只能計算得到面內應變，而力學計算方法能同時得到面內應變和面外應變的分布，在計算應變分布方面具有一定的優勢，本文只針對力學計算應變分布的方法進行研究。

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

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

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

圖3 所示為目標曲面的曲率半徑R分別為1.5，3.0 和4.5 m 時，2 種計算方法得到的曲面中心截面處的位移場計算結果?？梢钥闯觯呵市r，直接施加位移場計算得到的形狀與目標方程很接近，應變的計算結果滿足要求；曲率大時，計算得到的形狀與目標形狀的方程出現較大的偏差。產生這種現象的原因是：z方向位移場是根據初始節點的坐標計算得到，大變形時節點的x，y位置發生較大偏移，而計算過程中每個節點的z值不會發生改變，導致計算的曲面與目標曲面發生偏離；模具施加位移場的方法在計算精度上可以滿足要求，但在建模的復雜程度和計算效率上存在明顯的不足。

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

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

由于目前準確計算位移場的方法在計算精度和計算效率方面還存在局限性，使用ABAQUS 子程序DISP 進行二次開發，改進位移場的加載方式，可以準確獲取目標形狀的位移場，從而獲取準確的應變分布。為了方便描述，將該方法命名為子程序施加位移場的計算方法。子程序施加位移場的計算方法是根據目標形狀擬合曲面，逐步施加位移場載荷，每個載荷步計算結束后獲取節點當前坐標，根據當前坐標重新計算位移載荷，可以使計算的形狀接近目標形狀，使計算的應變結果滿足要求。圖4 所示為子程序施加位移場的流程，子程序的主要作用是更新當前節點坐標、根據當前節點坐標計算節點位移場。圖中：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 的圓形平板施加位移場，由于目標曲面曲率半徑R=1.5 m時位移場的結果存在明顯差異，因此選取目標曲面曲率半徑R=1.5 m 對子程序計算方法進行驗證。使用3 種方法進行彈性大變形計算，計算結果如圖5 所示。從計算結果看，直接施加位移場的方法相對于目標形狀會出現明顯的偏差；模具施加位移場與子程序施加位移場都能很好地接近目標形狀；同時比較了3 種計算方法的效率，表1所示為3 種方法在相同計算條件下的計算時長，子程序施加位移場的計算效率只是略低于直接施加位移場的計算效率，遠高于模具施加位移場的計算效率?？紤]建模的復雜性以及計算的效率，子程序施加位移場的方法具有明顯優勢。

圖5 位移場計算結果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 材料非線性對應變分布的影響

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

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

表2 Q235 鋼材料參數Table 2 Q235 steel material parameters

表3 Q235 鋼塑性應力-塑性應變數據Table 3 Data of Q235 steel plastic stress-plastic strain

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

圖6 彈性大變形計算面內應變分布Fig. 6 In-plane strain distribution of elastic calculation

圖7 彈塑性大變形計算面內應變分布Fig. 7 In-plane strain distribution of elastoplastic calculation

圖8 面內應變分布曲線Fig. 8 Distribution of in-plane strain

圖9 面外應變分布曲線Fig. 9 Distribution of out-plane strain

圖10 面內位移分布曲線Fig. 10 Distribution of in-plane displacement

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

4 適用性分析

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

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

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

圖11 曲率半徑與面內應變Fig. 11 Radius of curvature and in-plane strain

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

5 實例驗證

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

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

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

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

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

6 結 論

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

1） 目標形狀的曲率較大時，使用直接施加位移場的方法會導致位移場的計算結果不準確，使用模具施加位移場的方法建模復雜，計算時間長，子程序施加位移場的方法能很好解決建模復雜和計算時間長的問題。

2） 材料非線性會導致彈性大變形與彈塑性大變形計算的位移場結果出現細微的差別，導致應變分布出現差異，但應變分布規律基本一致，考慮到計算的效率，建議使用彈性大變形的計算方法來計算應變。

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