0 引言
对于圆盘类零件的轴对称塑性力学问题，LEE-WU[1]基于薄板理论并在平面应力及简单加载等假设条件下，给出了直接积分求解的方法，从而改进了需不断迭代才能得到解析解的方法[2]。刘幺和等[3]、高国华[4]将这种解法用于求解轴对称冲压成形中法兰区应力应变。秦泗吉等[5-6]在此基础上，将直接积分解法从求解平面内的轴对称问题推广至一般轴对称曲面问题。寻求理论解有助于进一步分析板材成形过程中的破裂及起皱等基础问题[7-8]。
对板材轴对称成形问题，采用平面应力假设比采用平面应变假设条件更符合实际[5,9]，但解法中所采用的比例加载条件在大变形情况下与实际相差较大，有时不能得到令人满意的结果。
文献[1,3-5]选择了现时构形为参考构形，在比例加载等假设条件下采用直接积分解法进行求解，这种方法可使所描述的问题直观化，且由于应变的计算是一次加载得到的，不需要跟踪变形质点，因此方程式和求解过程也相对简单。而当采用增量理论进行分析计算时，需要跟踪变形质点，以方便累加计算变形质点的应变，显然，这种情况下，应选择初始构形作为参考构形。
本文基于增量理论，在平面应力和薄板理论等假设条件下，选择初始构形为参考构形，通过引入参变量建立应力和应变之间的关系，分析得到了轴对称成形问题的应力参数方程、等效应变增量参数方程、协调方程，以及微分平衡方程等。在此基础上，化简消元得到了等效应变增量和参变量的二元方程。以圆筒形件和圆锥形件的成形为例，采用直接积分数值解法求出了各变形区的应变分布，分别对增量理论和比例加载假设条件下的两种方法得到的解析结果进行了分析，并与实验结果进行了比较。
1 基于增量理论的基本方程
1.1 应力参数方程和应变增量参数方程
假设板材面内同性、厚向异性，当厚度方向的应力为0时，根据等效应力的定义[9]， 有
(1)
式中，σρ、σθ、σ分别为径(经)向应力、周向应力和等效应力；R为板厚方向性系数。
根据增量理论，应变增量与应力之间的关系如下：
(2)
式中,
ερ、εθ和ε分别为同一变形质点经历了时间T后的径(经)向应变增量、周向应变增量和等效应变增量。
令
引进参数ω，式(1)可以用如下的参数方程表示：
(3)
将式(2)代入式(3)，可以得到应变增量的参数方程：
(4)
可以验证，式(4)满足等效应变增量的定义：
(5)
1.2 基于初始构形的变形协调方程
如图1所示，从圆盘或圆环的平板毛坯变形为轴对称壳体，采用初始构形为参考构形，设变形前某一微圆环初始内外径分别为ρ和ρ+dρ，若对应半径为ρ的质点变形后的径向位移为u，则变形后微曲面形圆环的外缘径向尺寸为ρ+u+dρ+du。α、α+dα分别为上下纬端面处的母线切线与对称轴的夹角。
图1 初始构形下的轴对称成形问题变形分析图
Fig.1 Deformation analysis diagram of axisymmetri c forming problem under initial configuration
由应变的定义可以得到径向应变ερ和周向应变εθ分别为
(6)
(7)
消去u后，得
(8)
将α=
/2代入式(8)可得平面轴对称问题的协调方程：
(9)
式(8)可用增量应变形式表示如下：
(10)
式中，ερ0和εθ0分别为增量加载前的径向和周向应变，它们是质点初始位置的函数。
1.3 微分平衡方程
在轴对称壳体零件成形中，每一个变形质点的主轴方向为经向、纬向及法向，对应的三个方向的应力分别表示为σρ、σθ和σz。如图2所示，参照文献[5-6]的分析方法，在二个相邻的纬锥面上截取一微锥壳体，然后沿轴对称线剖开。图中，ds为微锥壳体的经向弧长；ρ+u、ρ+u+dρ+du分别为微锥壳体的上下端纬面至对称轴的距离； σρ、σρ+dσρ分别为上下纬端面上作用的经向应力；σθ为作用于微锥壳体上的纬向应力。
图2 半微锥形圆环的受力分析图
Fig.2 Force analysis diagram of semi-micr o conical ring
分别以t、t+dt表示上下纬端面的厚度。设作用于壳体内表面的单位压力为p，以半微锥环为研究对象，分别在轴线方向和剖面的法线方向列平衡条件，得
(11)
消去p可得
d[σρt(ρ+u)]=σθtsinαds
(12)
由于ds =dρ/sinα，因此
d[σρt(ρ+u)]=σθtd(ρ+u)
(13)
展开式(13)，利用式(6)和式(7)消去u，并根据体积不变条件，得
(14)
采用增量求解方法时，式(14)可进一步表示为
ρdσρ/dρ-ρσρd(δερ+δεθ+ερ0+εθ0)/dρ =
(σθ-σρ)sinαexp(δερ-δεθ+ερ0-εθ0)
(15)
显然，分别以初始构形和现时构形为参考构形得到的微分平衡方程是完全不同的。
1.4 材料应力应变关系
设材料应力应变关系符合Hollomon硬化法则
σ=Bεn
(16)
式中，B为材料的强度系数；n为硬化指数。
式(16)中的等效应变ε可由下式得到：
(17)
(18)
ερ、εθ可用式(4)代入，因此，ερ和εθ可表示成等效应变增量和参变量的函数。
若用ε0表示增量加载前的等效应变，由于一般情况下应变不成比例增加，因此，ε≠ε0+ε。即需要先按式(18)计算出各应变分量，才能由式(17)计算得到等效应变。这说明，同一质点的各应变分量可以累加计算，而等效应变的计算则不能累加。这使得采用增量理论方法求解更加复杂。
式(3)、式(4)、式(10) 、式(15)和式(17)给出了包含σ、ε、σρ、σθ、ερ、εθ、ρ以及ω共8个变量7个方程，当边界条件和参数ω给定时，方程可解。
因变形质点初始位置ρ与变形瞬间的坐标位置有一一对应关系，因此求得的结果可容易地转换成各变量与坐标位置的关系。
2 直接积分解法
2.1 求解方程
由式(3)、式(16)分别可得
(19)
(20)
由式(10)和式(15)消去ρ，并将式(3)、式(19)和式(20)代入，化简后可得
(21)
另由式(4)可得
(22)
由式(17)可知，等效应变及其一阶导数都可以表示成等效应变增量和参变量的函数，设
可表示为
(23)
将式(22)和式(23)代入式(21)，得
(24)
a0=dεθ0/dω a1=sin(ω+β) a2=cos(ω+β)
fa=exp(ερ0-εθ0-2sinωcosβδε)/sinα-1
b0=-tan(ω+β)-d(ερ0+εθ0)/dω+nf0
b1=-2sinβcosω+nf1 b2=2sinβsinω+nf2
fb=2sinωsinβsinαexp(ερ0-εθ0-2sinωcosβδε)/a2
式(24)可进一步表示为
(25)
令
是包含了变量ε、dε/ω及ω的函数，则式(25)可写成
(26)
式(26)是含有ε和ω的二元方程，由于方程复杂，一般需根据边界条件采用数值方法才能求解。
2.2 直接积分解法及收敛性分析
对平面内的轴对称问题，采用比例加载假设条件时，因应变的一阶导数可以表示成应变和参变量的显式函数，因此可方便地采用直接积分解法进行求解[1,3-5]。同样地，采用增量理论时，对轴对称平面内的成形，如拉深过程中的法兰区，α为定值，式(26)等式右端不含dε/ω项，这样可以由ε及ω直接求出dε/ω，因此也能采用积分解法直接求解。
对一般轴对称曲面零件的成形，当成形制件的形状一定时，α可以表示成质点位置ρ的函数。而变形质点ρ又是ε、dε/ω及ω的函数，因而α也是ε、dε/ω及ω的函数。这样，式(26)的左右端都包含应变的一阶导数，一般来说，已知ε和ω，需要反复迭代才能求出dε/ω，这给求解过程带来不便。
参照文献[5]在比例加载假设条件下的分析方法，采用增量理论对一般轴对称曲面零件成形问题进行求解，其求解过程和收敛性分析如下：
根据泰勒级数展开式，若应变增量的一阶导数dε/dω存在，则有
(27)
式中，
ω为参变量ω的增量；o(ω)为ω的高次项。
将式(26)代入式(27)，得
(28)
设ωi+1是ωi的邻近点，由式(28)可知
(29)
ωi+1=ωi+1-ωi
对于平面内的轴对称成形问题，应变增量的一阶导数dε/dω可直接表示成ε、ω的函数，即式(29)中的函数f仅是应变增量和参变量的函数，因而可由εi和
ωi+1直接求出，进而可得到εi+1的近似值。
一般情况下，因方程的左右端都含等效应变增量的一阶导数，由式(26)不能直接得到dε/dω与ε、ω的关系，需重新考查积分求解过程。
由泰勒展开式可知，若ε的二阶导数存在，即函数f(dε/dω,ε,ω)对dε/dω的一阶偏导数存在，则
(30)
将式(30)代入式(29)，得
(31)
假设ωi+1=ωi+1-ωi=ω，式(31)可进一步简化为
(32)
这样，εi+1可近似表示成εi、ωi和dεi-1/dωi-1的函数：
(33)
对一般轴对称成形问题，当ε0、ω0和dε0/dω0初始值给定时，可根据直接积分法逐步计算出ε1、ε2直至εn。对于任意给定的ω，等效应变增量ε的计算步骤如下：
(1)将区间(ω0, ω)等分成N段，则ω=(ω-ω0)/N，ωi=ω0+iω(i=1，2，…,N)。εi是对应ωi的等效应变增量。
(2)ε1可近似用下式计算：
(34)
(3)由式(33)可知，ε2可近似用下式计算:
(35)
(4)一般地，i(i≥2)为任意值时，εi都可近似用式(33)计算。考虑到
和
则对于任意ω(ω>ω1)对应的应变增量ε近似为
(36)
式(36)中省略掉的是比ω高级的微分项，N越大，计算精度越高，因此，式(36)是收敛的，可根据需要，得到足够精确的求解结果。
3 应用举例
为了进一步说明和验证上文的理论分析过程，首先以圆筒形件的拉深成形为对象，采用理论分析方法计算法兰区和凹模圆角区的应变分布并与实验结果进行对照。
理论分析和实验所用的材料为ST16。性能参数为：强度系数B=511.4 MPa，硬化指数n=0.26，板厚方向性系数R=2.243。初始板坯尺寸：直径为110 mm，厚度为0.87 mm。成形凸模外径为50 mm，凸模圆角为5 mm，凹模圆角为9.1 mm。压边力为10 kN。
设压边力为F，摩擦因数为μ，若摩擦力全部作用于法兰的外缘，即当位置半径为Rw时，径向应力σρ=F/(Rwtw)，其中，tw为法兰外缘对应的板坯厚度。边界条件为：等效应变ε0=ln(R0/Rw)，ω0=2-arccosγ-β。其中，γ=
当不考虑摩擦，径向应力初值为0时，ω0=3/2-β。考虑实验中采用了薄膜润滑条件，与文献[5]一致，取μ=0.06。
由于在法兰区α=/2，根据给出的初始应变和参变量边界条件，可由式(26)直接求出等效应变增量的一阶导数，进而逐步求出法兰区的应变增量和应变。当计算进行至凹模圆角区时，因其法兰区邻近点的等效应变增量和一阶导数都是已知的，故仍可根据前面的求解方法得到等效应变增量的一阶导数，从而完成后续点的应变求解过程。在每一个加载步，N值全部取1 000(计算表明，N值继续增大时，计算精度没有显著变化)。
在材料模型参数、模具尺寸、板坯尺寸以及压边力等成形条件与实验相同的情况下，分别采用增量理论和比例加载假设条件的理论方法进行分析，将计算结果与实验结果进行比较。实验采用文献[5]给出的方法和结果，限于篇幅，具体实验和测量方法本文不再赘述。
采用增量理论进行计算时，法兰外缘相对位置半径r/R0加载至0.858时，增量加载步数取5步，相对位置半径至0.807时，增量加载步数再增加2步。每一个增量步按几何平均值选取。
图3为采用增量理论逐步加载计算得到的径(经)向应变和周向应变沿坐标位置的分布曲线。曲线1～7分别对应1～7次的加载步。横坐标表示径向相对坐标值r/R0。计算结果显示，变形程度较小时，凹模圆角区应变绝对值小于邻近的法兰区应变绝对值，但随着变形程度的不断增大，凹模圆角区的应变绝对值逐渐增大直至最大，这与拉深过程中的实际情况是一致的。
图3 应变分布图(增量理论计算值)
Fig.3 Strain distribution curves
(by incremental theory)
可将两向应变以应变状态图的形式表示，如图4所示,图中曲线1～7分别对应1～7次的加载步。借助于应变状态图，更容易判断两向应变间的关系、加载路径和变形方式，以及便于进行成形极限分析等。图4表明，在加载过程中，法兰区的应变增加较平缓，而凹模圆角区应变增加较剧烈。图4中还给出了初始相对位置ρ/R0在0.607～1之间的5个变形质点在不同加载步下的应变值。ρ/R0=1表示法兰外缘的质点，ρ/R0=0.607表示接近凹模口的质点,ρ/R0介于0.864～1之间为法兰区，ρ/R0介于0.864～0.607之间为凹模圆角区。
图4 应变状态图(增量理论计算值)
Fig.4 Diagram of strain state
(by incremental theory)
为了更清楚地表示同一质点在加载过程中两向应变的关系，将各质点在逐次加载中的应变值另表示在图5中。分析表明，随着变形过程的进行，在法兰区的变形质点应变接近成比例增大，而在凹模圆角区，变形质点应变不完全符合按比例增大的条件。
图5 部分变形质点两向应变关系
Fig.5 Relationship between two-direction strain o f partially deformed particles
图6所示为拉深件法兰外缘相对半径Rw/ R0分别为0.858和0.807时，分别采用比例加载假设条件和增量理论两种方法计算得到的应变值以及实验值沿径向坐标位置r/R0的分布情况。可以看出，两种理论方法计算得到的结果在法兰区相差很小，这进一步验证了图5的分析结果。周向应变在整个区域都相差很小(两条曲线基本重合)，径向应变在法兰区相差较小，在凹模圆角区，采用增量理论求解得到的计算结果稍小。图6中的离散点为实验结果，即采用增量理论计算得到的结果更接近实验值。
图7还给出了圆锥形件拉深成形各变形区应变分布的理论值和实验值，法兰外缘相对半径Rw/R0分别为0.878和0.838。理论计算分别采用了比例加载和增量理论的方法，其中增量理论总加载次数为7。实验值采用了文献[6]的实验结果，理论计算的边界条件也与文献[6]一致。图中光滑曲线和离散点分别表示理论值和实验值。
图7表明，在法兰区、凹模圆角区及悬空侧壁区的应变分布理论计算值和实验结果分布趋势吻合，采用增量理论得到的计算值更接近于实验结果，其中法兰区和凹模圆角区应变分布规律与圆筒形件的结果一致。
(a)Rw/R0=0.807
(b)Rw/R0=0.858
图6 圆筒形件应变分布理论值和实验结果的比较
Fig.6 Comparison of theoretical and experimenta l values of strain distributions of cylindrical part
(a)Rw/R0=0.878
(b)Rw/R0=0.838
图7 圆锥形件应变分布理论值和实验结果的比较
Fig.7 Comparison of theoretical and experimenta l values of strain distributions of conical part
图6和图7还表明，在凹模圆角区的径向应变计算值与实验值相差稍大，这主要是因为板坯在凹模圆角区发生了弯曲变形。而对于圆锥形件而言，由于板坯离开凹模圆角后又产生了反向弯曲，使得悬空侧壁区的计算值又稍接近于实验值。对于圆筒形件的成形，文献[5]考虑了凹模圆角弯曲的影响，分析得到的理论结果与实验值更吻合。限于篇幅，这里不讨论板坯内外层的区别以及板坯经过凹模圆角因产生弯曲和反弯曲而引起应变的变化。
前文分析过程中所给出的基本方程适用于轴对称拉深、胀形和翻边等成形方式，因而该方法在添加适当的边界条件后，也适用于这些变形方式对应问题的求解。
4 结论
(1)对于一般轴对称成形问题，在平面应力和薄板理论等假设条件下，以初始构形为参考构形，采用参数分析方法，根据平衡方程、变形协调方程、增量理论以及材料的等效应力应变关系等方法分析得到了等效应变增量的微分方程。
(2)对于轴对称曲面零件的成形问题，给出了当等效应变增量一阶导数不能表示成显式函数时的直接积分解法，并对收敛性进行了分析。
(3)以圆筒形件和圆锥形件的拉深成形为例，将总变形分成7个增量加载步，采用增量理论解法，求解得到了法兰区和凹模圆角区的应变分布。结果表明，变形过程中法兰区的变形质点接近比例加载条件，而凹模区的变形质点采用增量理论求解更接近实验值。
[1] LEE-WU M H. A Simple Method of Determining Plastic Stresses and Strains in Rotating Disk with Nonuniform Metal Properties[J]. Journal of Applied Mechanics, 1952,19(4):489-495.
[2] STOWELL E Z. Stress and Strain Concentration at a Circular Hole in an Ininite Plate[J]. National Advisory Committee Aeronautics, 1950,19(9):1-14.
[3] 刘幺和,胡世光. 板料的n值R值对应力应变分布的影响[J]. 锻压技术, 1982(3):11-18.
LIU Yaohe, HU Shiguang. Influence of Values of n and R on Stress and Strain Distribution[J]. Forging, 1982(3): 11-18.
[4] 高国华. 轴对称平面应力问题的塑性解及其在板料冲压中的应用[J]. 应用力学学报, 1987, 4(4): 109-132.
GAO Guohua. An Analysis for Axisymmetric Plane Stress Problem and Its Application to Sheet Metal Forming Process[J]. Journal of Applied Mechanics, 1987, 4(4): 109-132.
[5] QIN Siji, GAI Binbin, KONG Xiaohua, et al. Analytical Solutions of Strain of Axisymmetric Curved Part in Sheet Metal Forming Process Using Direct Integral Method[J]. International Journal of Mechanical Sciences, 2015, 101/102:49-58.
[6] 秦泗吉, 邓超, 杨莉, 等. 圆锥形件拉深成形应力应变的直接积分解法[J]. 中国机械工程, 2016, 27(2):251-257.
QIN Siji, DENG Chao, YANG Li, et al. A Direct Integral Method to Solve Stress and Strain of Conical Parts in Deep Drawing Process[J]. China Mechanical Engineering, 2016, 27(2): 251-257.
[7] HOUSSEM B, CARL L, KHEMAIS S. Ductile Damage Prediction in Sheet and Bulk Metal Forming[J]. Comptes Rendus Mecanique, 2016, 344(4/5): 296-318.
[8] MAZIAR K, ABDOLHAMID G, MOHAMMAD B, et al. Investigation of Wrinkling in Hydrodynamic Deep Drawing Assisted by Radial Pressure with Inward Flowing Liquid[J]. Procedia Engineering, 2017, 183:65-70.
[9] 梁炳文,胡世光. 板料成形塑性理论[M]. 北京:机械工业出版社, 1987.
LIANG Bingwen, HU Shiguang. Sheet Metal Forming Theory of Plasticity[M].Beijing: China Machine Press, 1987.
