江熒，彭炎榮

(1.南京航空航天大學能源與動力學院，江蘇 南京 210016;2.蕪湖職業(yè)技術(shù)學院機械工程系，安徽 蕪湖 241006;3.湘潭大學機械工程學院，湖南 湘潭 411105)

0 前言

滑移線場理論是利用理想彈塑性材料塑性變形過程中最大剪應力跡線的性質(zhì)，求解塑性力學邊值問題的一種嚴密方法。它滿足塑性力學平面問題的所有控制方程和邊界條件，作為經(jīng)典塑性理論的重要組成部分，廣泛地應用于金屬塑性成型、結(jié)構(gòu)塑性極限分析和土力學等工程領域［1］，文獻［2］和文獻［3］描述了滑移線場理論在金屬成型領域中的具體應用實例。利用滑移線方法模擬金屬成形過程，關鍵在于建立準確的滑移線場。若滑移線場是由均勻場、中心場、漸開線場、對數(shù)螺線場和擺線場等構(gòu)成，其建立并不困難，因為這些場中的滑移線均可用精確的解析式來描述。但對于復雜滑移線場，如，有心扇形場，其滑移線就沒有準確的方程來表達，這就給求解塑性成型問題帶來很大困難。而這方面的研究尚未有詳細、系統(tǒng)的分析資料。本文將系統(tǒng)介紹構(gòu)建復雜滑移線場的三種理論求解方法，并給出具體計算方法和解題步驟。

復雜滑移線場中的有心扇形場在求解實際問題時應用極為普遍，而且形式各異，板坯在粗糙平砧間鐓粗，或者在光滑楔形模中擠壓時，有心扇形場兩基線圓弧半徑r1和r2相等，如圖1(a)和圖1(b)所示，而板坯在粗糙楔形模中擠壓時，兩基線圓弧半徑r1和r2卻不相等，如圖1(c)所示;板坯在粗糙平砧間鐓粗時，兩圓弧對應的中心張角θ和ψ相等，如圖1(a)所示，在光滑或者粗糙楔形模中擠壓時，兩圓弧對應的中心張角θ和ψ又不相等，如圖1(b)、圖1(c)所示。

圖1 三種復雜有心扇形滑移線場案例

根據(jù)兩基線圓弧半徑r1、r2和對應的中心張角θ和ψ建立有心扇形場，其數(shù)學理論屬于求解雙曲型偏微分方程組的第一類邊值問題(黎曼問題)，可以通過理論計算求解或作圖方法求解?，F(xiàn)以有心扇形場為例，系統(tǒng)介紹三種理論計算求解方法?？紤]到圖1(c)所示的非對稱型(r1≠r2，θ≠ψ)扇形場的應用更為普遍，而文獻中論述又較少，故本文選擇非對稱有心扇形場為例，其求解的基本原理和方法也可以推廣應用于其他復雜形式的滑移線場。

1 數(shù)值積分法

1.1 基本原理

圖2為滑移線單元網(wǎng)格。

圖2 滑移線單元網(wǎng)格

基于滑移線方程來構(gòu)建復雜滑移線場的數(shù)值積分法基本原理簡述如下:

先將圖2所示滑移線方程:

改寫為差分方程的形式:

式(2)中所有下標均表示圖2中的對應節(jié)點。

由漢基第一定理知:

由式(2)可解得點(m，n)的坐標值:

如果已知圖2中曲邊矩形表示的單元網(wǎng)絡中3個節(jié)點的ω值和坐標位置，利用式(3)和式(4)就可以求得第4個未知節(jié)點的ωm，n和坐標(x，y)值。

1.2 討論

數(shù)值積分的特點是所用的差分方程原理簡單易懂，容易掌握。且分度角△越小則網(wǎng)絡越細，求解結(jié)果越準確［4］，但計算工作量隨之增加。

該方法的缺點是必須知道單元網(wǎng)格的3個節(jié)點，才能求得第4個節(jié)點。因此，計算只能從節(jié)點(0，0)開始，先沿基圓進行，再依次向有心扇形場內(nèi)部逐點“擴散”，不能直接計算場中某指定節(jié)點(m，n)的坐標位置，也不能直接確定塑性區(qū)的邊界滑移線，這在實際應用中不便。

2 解析求解法

2.1 基本原理

圖3為解析計算法求解有心扇形場。

圖3 解析計算法求解有心扇形場

解析求解法在一般書刊文獻中都沒有介紹。其基本原理是先通過解偏微分方程求得滑移線曲率半徑的近似解析解，再利用曲率半徑與節(jié)點位置的關系求得各節(jié)點的坐標。求解過程簡述如下:

由漢基第二定理知，圖3中α和β兩族滑移線的曲率半徑R和S之間有下述關系:

從而有:

將式(6)帶入式(5)，即得到一組雙曲型方程:

式(7)又稱為電報方程，可用純粹的解析方法求解。如對式(7)中的第一式求解，即得到α線在節(jié)點 P(θ，ψ)處的曲率半徑 R(θ，ψ)為:

式(8)中J0為零階貝塞爾函數(shù)，角度α'和β'如圖3。

R(0，0)也就是節(jié)點(0，0)處 α1基線圓弧的半徑R0(考慮到α1弧線方向為順時針走向，曲率半徑規(guī)定為負值，故R(0，0)=-R0，而β1基線圓弧的半徑則相反，其S(0，0)=+S0。

由漢基第二定理知，沿β1圓弧上各點的α線的曲率半徑R為:

將式(9)代入式(8)積分可得:

其中，I0和I1為第一類變型貝塞爾函數(shù)。

為計算簡便，在以下推導中，令式(10)的R0=1(單位長度)，于是 S0=c。

現(xiàn)選定(x'，y')直角坐標系，如圖3。由微分幾何原理知，圖3中與張角(θ，ψ)對應的節(jié)點P的坐標可通過該點的曲率半徑 R 來確定，其關系式為:

其中:

將式(10)代入式(11)進行分部積分，最后得到滑移線場中任意節(jié)點 P(θ，ψ)處的坐標，)的表達式如下:

其中:

2.2 求解方法

在具體計算有心扇形場中所有節(jié)點的坐標時，一般先將α1和β1基線圓弧按Δ微小角度等分，于是有:

將式(15)代入(14)式得:

為了便于和數(shù)值積分法比較，應將圖3中的(x'，y')坐標進行轉(zhuǎn)換，坐標轉(zhuǎn)換公式如下:

2.3 討論

因式(13)中的級數(shù)收斂性很好，故這種方法求解過程簡單，獲解迅速，求解精度很高。

這種方法突出的優(yōu)點是不需依次逐點進行數(shù)值計算。由式(13)就可以直接求出滑移線場中與任意一組(θ，ψ)或(m，n)值對應的某節(jié)點P的坐標值(xmn，ymn)，也可以直接求出塑性區(qū)的邊界滑移線，而不必計算出場中整個網(wǎng)絡的全部節(jié)點，這在求解實際問題時是非常適用的。

3 矩陣算子法

矩陣算子法是將α和β兩族滑移線的曲率半徑R和S用均勻收斂的雙冪級數(shù)表示，級數(shù)的系數(shù)則用列向量X表示，應用矩陣算子法和疊加原理即可求解滑移線場。此法的基本原理在文獻［5］和文獻［6］中已有詳細的論述，本文只針對有心扇形場的建立，說明其求解過程。

3.1 先求各滑移線的系數(shù)列向量X

圖4中，已知基線圓弧OA的半徑R0=1和圓弧OB的半徑S0=c，故其系數(shù)列向量XOA和XOB分別為:

與任意一組張角(θ，ψ)對應的滑移線AP和BP的列向量分別為:

XBP和XAP與基線圓弧的列向量XOA和XOB有以下關系:

其中，P*和Q*均是矩陣算子，定義如下:

以上矩陣中的元素ψn和θn分別為:

3.2 求AP滑移線上各點(以P點為例來說明)的坐標(x'，y')

先對AP線建立專用的所謂Mikhlim坐標系。這是一種轉(zhuǎn)動坐標，原點為A點，而ξ和η兩坐標軸分別與AP線在P點處的切線和法線平行(圖4)。顯然，隨著P點在AP線上位置的變動，ξ軸和η軸也隨之轉(zhuǎn)動。P點的坐標(ξ，η)值按下式計算:

其中，tn由AP線列向量XAP的元素Sn來確定，其計算公式為:

從而求得:

將這一系列t值代入式(20)，即可求得P點的坐標(ξ，η)值。

圖4 矩陣算子法求解有心扇形場

再將坐標(ξ，η)換成固定坐標(x'，y')，固定坐標(x'，y')仍以A為原點，但以AP線在A處的切向和法向作為x'軸和y'軸的方向(圖4)。將式(20)求得的(ξ，η)值代入下述坐標變換公式，即得 P 點的(x'，y')值。

最后將專用的固定坐標(x'，y')換成圖4所示的統(tǒng)一坐標(x，y)。如，A點坐標(xA，yA)為:

AP線上任意一點P的坐標(xp，yp)為:

3.3 求BP滑移線上各點(仍以P點為例)的坐標

求解BP滑移線上各點坐標的方法同上，但須注意以下幾點:

1)應以 B點為原點建立 BP線的(ξ，η)和(x'，y')專用坐標系。

2)因BP線屬于順時針旋轉(zhuǎn)走向，曲率半徑應為負值，故對式(20)的符號應作如下修改:

式(25)中的系數(shù)t則應由BP線列向量XBP的元素rn來確定，即:

而式(22)則應改為:

3)由式(26)求得專用固定坐標(x'，y')后，應按式(27)變換成統(tǒng)一的坐標(x，y)。如，B點坐標為:

而BP線上任意一點P的坐標為:

當然，由式(28)求得的結(jié)果，應與式(24)的結(jié)果一致，因為二者都是同一點P的坐標值。

3.4 節(jié)點(m，n)的坐標值

若將兩基線圓弧按△微小角度等分，則只要將θ=mΔ、ψ=nΔ代入上述計算過程，由式(28)和式(24)求得的(xp，yp)值也就是節(jié)點(m，n)的坐標(xmn，ymn)值。

3.5 討論

所有列向量X和矩陣算子P*、Q*，雖屬無窮級數(shù)，但實際計算中，取6×6階矩陣已夠精確。

矩陣算子法不僅可以求解滑移線場，還可以計算滑移線場中力的分布和速度分布，特別是求解沒有已知滑移線，需根據(jù)速端圖或滑移線幾何性質(zhì)來推導的所謂間接型問題。在用滑移線方法模擬塑性成形過程時，矩陣算子法將起重要作用。

矩陣算子法求解過程看似繁瑣，但若能熟悉矩陣運算方法，特別是可以將種類繁多的矩陣算子預先編好子程序以待隨時調(diào)用，則方法仍然很便捷。

4 三種理論解法求解結(jié)果的比較

為了便于比較，以書刊文獻中常用的對稱型有心扇形場為例，并令，滑移線網(wǎng)絡按Δ=5°等分，張角θ=ψ=135°。筆者用三種理論解法對全部節(jié)點逐點進行了計算，現(xiàn)僅將對稱軸線(x軸)上節(jié)點坐標x用三種方法計算的結(jié)果列出，如表1所示。

由表1可見，三種方法所得結(jié)果相差極小。如果以數(shù)值積分法求解的結(jié)果為標準，則θ≤45°時，解析法求解的結(jié)果與數(shù)值積分法求解結(jié)果相比相差0.04%，矩陣算子法的求解結(jié)果與數(shù)值積分法的求解結(jié)果相比相差0.06%。θ≤90°時，解析求解法求解結(jié)果與數(shù)值積分法求解結(jié)果相比相差0.18%，矩陣算子法求解結(jié)果與數(shù)值積分法求解結(jié)果相比相差0.53%。即使θ=135°時，解析求解法的求解結(jié)果與數(shù)值積分法的求解結(jié)果相比也只差0.34%，矩陣算子法的求解結(jié)果與數(shù)值積分法的求解結(jié)果相比相差1.93%。

實際應用中，θ值一般不超過90°，因此三種理論的求解方法精度都相當高，完全滿足工程應用和精細研究的要求。

表1 有心扇形場對稱線上節(jié)點坐標x用三種方法求解結(jié)果的比較

［1］曾錢幫，劉彤，馬平.廣義Hoek-Brown破壞準則平面應變問題的滑移線場理論［J］.土木建筑與環(huán)境工程，2010，32(1):4-11.

［3］秦小瓊，劉德學.杯形件反擠壓變形力的滑移線場數(shù)值分析解［J］.山東大學學報:工學版，2011，41(6):66-69.

［3］嚴志遠，雷步芳，李永堂，等.基于滑移線理論計算板料彎曲變形力［J］.機械工程與自動化，2011，168(5):176-178.

［4］羅中華，羅文波，彭炎榮.一種高精度有心扇形滑移線場的近似解析解［J］.鍛壓技術(shù)，1999，24(1):43-45.

［5］Johnson W.Plane-strain Slip-line Field for Metal-deformation Processes［M］.Oxford:Pergamon Press，1982.

［6］王祖唐，關廷棟，肖景容，等.金屬塑性成形理論［M］.北京:機械工業(yè)出版社，1989:232-240.