江熒，彭炎榮

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

0 前言

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

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

圖1 三種復(fù)雜有心扇形滑移線場(chǎng)案例

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

1 數(shù)值積分法

1.1 基本原理

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

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

基于滑移線方程來(lái)構(gòu)建復(fù)雜滑移線場(chǎng)的數(shù)值積分法基本原理簡(jiǎn)述如下:

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

改寫(xiě)為差分方程的形式:

式(2)中所有下標(biāo)均表示圖2中的對(duì)應(yīng)節(jié)點(diǎn)。

由漢基第一定理知:

由式(2)可解得點(diǎn)(m，n)的坐標(biāo)值:

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

1.2 討論

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

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

2 解析求解法

2.1 基本原理

圖3為解析計(jì)算法求解有心扇形場(chǎng)。

圖3 解析計(jì)算法求解有心扇形場(chǎng)

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

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

從而有:

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

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

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

R(0，0)也就是節(jié)點(diǎn)(0，0)處 α1基線圓弧的半徑R0(考慮到α1弧線方向?yàn)轫槙r(shí)針走向，曲率半徑規(guī)定為負(fù)值，故R(0，0)=-R0，而β1基線圓弧的半徑則相反，其S(0，0)=+S0。

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

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

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

為計(jì)算簡(jiǎn)便，在以下推導(dǎo)中，令式(10)的R0=1(單位長(zhǎng)度)，于是 S0=c。

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

其中:

將式(10)代入式(11)進(jìn)行分部積分，最后得到滑移線場(chǎng)中任意節(jié)點(diǎn) P(θ，ψ)處的坐標(biāo)，)的表達(dá)式如下:

其中:

2.2 求解方法

在具體計(jì)算有心扇形場(chǎng)中所有節(jié)點(diǎn)的坐標(biāo)時(shí)，一般先將α1和β1基線圓弧按Δ微小角度等分，于是有:

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

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

2.3 討論

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

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

3 矩陣算子法

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

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

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

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

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

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

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

3.2 求AP滑移線上各點(diǎn)(以P點(diǎn)為例來(lái)說(shuō)明)的坐標(biāo)(x'，y')

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

其中，tn由AP線列向量XAP的元素Sn來(lái)確定，其計(jì)算公式為:

從而求得:

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

圖4 矩陣算子法求解有心扇形場(chǎng)

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

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

AP線上任意一點(diǎn)P的坐標(biāo)(xp，yp)為:

3.3 求BP滑移線上各點(diǎn)(仍以P點(diǎn)為例)的坐標(biāo)

求解BP滑移線上各點(diǎn)坐標(biāo)的方法同上，但須注意以下幾點(diǎn):

1)應(yīng)以 B點(diǎn)為原點(diǎn)建立 BP線的(ξ，η)和(x'，y')專(zhuān)用坐標(biāo)系。

2)因BP線屬于順時(shí)針旋轉(zhuǎn)走向，曲率半徑應(yīng)為負(fù)值，故對(duì)式(20)的符號(hào)應(yīng)作如下修改:

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

而式(22)則應(yīng)改為:

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

而B(niǎo)P線上任意一點(diǎn)P的坐標(biāo)為:

當(dāng)然，由式(28)求得的結(jié)果，應(yīng)與式(24)的結(jié)果一致，因?yàn)槎叨际峭稽c(diǎn)P的坐標(biāo)值。

3.4 節(jié)點(diǎn)(m，n)的坐標(biāo)值

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

3.5 討論

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

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

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

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

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

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

實(shí)際應(yīng)用中，θ值一般不超過(guò)90°，因此三種理論的求解方法精度都相當(dāng)高，完全滿(mǎn)足工程應(yīng)用和精細(xì)研究的要求。

表1 有心扇形場(chǎng)對(duì)稱(chēng)線上節(jié)點(diǎn)坐標(biāo)x用三種方法求解結(jié)果的比較

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

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

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

［4］羅中華，羅文波，彭炎榮.一種高精度有心扇形滑移線場(chǎng)的近似解析解［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］王祖唐，關(guān)廷棟，肖景容，等.金屬塑性成形理論［M］.北京:機(jī)械工業(yè)出版社，1989:232-240.