基于结合面基础特性参数的球面结合部模型研究

刘洁

(延安大学 化学与化工学院，陕西 延安 716000)

0 引言

对于结合面的各种结合条件，可以通过大量的基础实验来获得不同情况下(如不同的载荷大小和性质，结合面材料，结合面的加工方法及表面完整性，润滑油性质等)结合面的基础特性参数。关于如何运用基础特性参数解析计算实际机械结合部特性参数，即如何把结合面的基本假设，数学模型的理论运用于计算结合部特性参数的研究方面，国内外不少学者都在研究。文献[1]的作者提出假设，组成结合部的构件在结合部处的刚度比结合面的刚度高的情况下，忽略构件本身在结合部处变形的影响。这就使结合部复杂的变形计算问题转化成仅计算结合部的接触表面层的变形问题，从而使结合部的变形可由构件间的相对线位移和相对角位移来表示。这是在处理结合部变形这一复杂的非线性问题上的一个突破。文献[2]的作者，根据以上结合部的物理假设，建立了平面移动导轨结合部变形的通用数学模型，并编制了计算平面移动导轨结合部变形的通用计算软件。并以JISGNC数控车床横溜板为对象进行了实验验证，证明了上述结合部物理假设为基础所建立的数学模型是适合的，计算程序是正确的。文献[3]的作者，以满足机械系统解析对结合部的上述要求为出发点，建立了圆柱结合部变形的通用数学模型，编制了计算圆柱结合部变形通用计算程序，并通过实验验证了其理论和程序是正确的。

在解析球面结合部时，也应用了文献[5]提出的结合部物理假设。并通过借鉴圆柱面结合部的研究成果，进一步发展研究，提出了球面结合部的解析方法，编制了相应的解析程序，并通过了理论对比论证了其理论及相应的解析程序也是正确的。

1 球面结合部理论解析

1.1 球面结合部数学模型

图1所示的球面结合部的总体坐标系xyz-O，xyz-o为球面上任一点的局部坐标系。

图1 球面结合部计算坐标系

1.2 球面结合部位移条件及载荷结合条件的确定

1) 设球面上任意一点A处坐标(R,φ,θ)，其中x3表示在XYZ-O中z向坐标，圆球的半径为R，其中X3=R·cosφ，r=R·sinφ。球面结合部坐标系∑J—xyz-O，球面结合面坐标系∑S—xyz-o。

从结合部坐标系∑J到结合面坐标系∑S的变换可以依次通过以下变换得到：

一次平移变换：

(rcosθ,rsinθ,X3)；

两次旋转变换：

则球面结合面坐标系∑S相对于球面结合部坐标系∑J的齐次变换矩阵为：

(1)

其中：

(2)

(3)

式中：SRJ是球面结合面坐标系与球面结合部坐标系的旋转变换矩阵；

JPS是球面结合部坐标系与球面结合面坐标系的平移变换矩阵。

2) 设球面结合部变形位移列阵为：

{xJ}={xJ1,xJ2,xJ3,xJ4,xJ5,xJ6}T

(4)

式中：xJ1,xJ2,xJ3——分别为球面结合部x，y，z方向的线位移；

xJ4,xJ5,xJ6——分别为球面结合部x，y，z方向的角位移。

{xJS}=[SJJ]·{XJ}

(5)

式中： [SJJ]——球面结合面S与球面结合部J的坐标变换雅克比矩阵；

{XJS}——球面结合面的位移列阵。

则可由式(4)、式(5)可以求出球面结合面位移。

3) 由第一章和式(5)可以求出球面结合面上任一点处的接触变形：

λτ=(xJ1+Rcosφ·xJ5)sinθ+(Rcosφ·xJ4-xJ2)cosθ-Rsinφ·xJ6

λn=cosθsinφ·xJ1+sinθsinφ·xJ2+cosφ·xJ3λ3=cosθcosφ·

xJ1+sinθcosφ·xJ2-sinφ·xJ3-Rsinθ·xJ4+Rcosθ·xJ5

(6)

4) 球面结合部反力及反力矩：

FRK=JJS·FRS

(7)

式中：FRK—作用在球面结合部上的反作用力及反力矩(K=1,2…6)；

JJS—球面结合部∑J与球面结合面∑S的力雅克比坐标变换矩阵；

FRS—球面结合面上的反作用力。

当结合部受力位移后，作用于球面上的反载荷可由求出：

(8)

1.3 球面结合部刚度矩阵的确定

当结合部受外力作用发生变形位移{XJ}后，结合部所受的力{FJ}={F1,F2,F3,F4,F5,F6}T

可表示如下：

F1=-FR1=

F2=-FR2=

[(ατλnβτλ3sinα+αnλnβncosα)sinθ-ατλnβτλτcosθ]X3}ds

[(ατλnβτλ3sinα+αnλnβncosα)cosθ-ατλnβτλτsinθ]X3}ds

(9)

其中：式(9)的λn，λτ，λ3与式(6)相同。

根据刚度的定义，可以用Fi(i=1,2,…6)分别对XJj(j=1,2,…6)求偏导数而求出刚度矩阵中的每个元素kij。即：

(10)

若以λn，λτ，λ3为中间变量，则kij可表示如下：

(11)

现以k21为例推导出如下：

cosθcosα+αnβnλnβn-1cos2αsinθcosθ-ατβτλnβτ-1λτcosθcosθcosα+ατλnβτsin2αsinθcosθ)φ(λn)rds

(12)

由第一章可知：

(13)

ατλnβτ=kτ

(14)

所以：

(15)

根据上面的方法可以求出圆锥面结合部刚度矩阵[KJ]中的每个元素kij。因篇幅所限，这里省略其他刚度值的推导过程。总之，根据式(10)可以推导出球面结合部的刚度矩阵。

2 球面结合部解析程序及理论验证

对前面提出的关于球面结合部的数学模型中，含有六个待求的未知量XJi(i=1,2…6),但是由于这六个未知量都是隐含的，无法从方程组中分离出来。对于隐含的非线性方程组的求解问题，目前一直没有一种通用的计算方法来加以解决，在借鉴了文献[5]中的计算程序所用的计算方法的基础上，采用了求解变量轮换弦截法，对球面结合部的数学模型，编制FORTRAN语言程序来求解。所谓坐标轮换弦截法的基本思想就是在给出欲求量赋初值后，通过弦截法来改进欲求量的值，对求解的数学模型中的每个方程轮换地进行验算，一直到每个方程都满足给定的精度后，输出欲求量。

2.1 球面结合部解析程序编制说明

a) 程序的基本功能

1) 确定球面结合部的位移结合条件；

2) 确定球面结合部的载荷结合条件；

3) 确定球面结合部的静刚度及动刚度；

4) 确定球面结合部的阻尼特性。

b) 程序结构流程图(图2)

图2 程序流程图

c) 程序各模块功能

1) 主程序模块：主要用来输入控制参数、基础特性参数、结构参数、载荷参数等。

2) 位移解析模块DISPLACE：用于确定结合部位移结合条件。

3) 反力计算模块REACTF：用于计算结合部反力及确定结合部载荷结合条件。

4) 面压计算模块PRESSURE：用于确定各结合面的载荷结合条件。

5) 刚度计算模块STIFF：用于确定各结合面的刚度及结合部总刚度。

6) 阻尼计算模块DAMP：用于确定结合面的阻尼特性。

2.2 球面结合部实例验证

a) 实例验证(图3)

图3 球面结合部坐标系

设作用于球面及圆柱面结合部的作用力均为：Fx=1400N，其余方向的作用力及力矩全为零。通过编制的程序分别比较这两种情况的计算结果：

1) 设球面的半径R=0.82m，其中取φa=89°,φb=91°，则通过编制的球面结合部程序计算输出结果为：

STATIC STIFFNESS MATRIX OF JOINT SURFACE:

IJS=1

0.12451E+06 -0.68760E-02 -0.10690E-03 0.0000 -0.11573E-02 0.63331E-01

-0.68760E-02 0.12451E+06 0.16288E-05 0.11573E-02 0.0000 492.06

-0.10690E-03 0.16288E-05 25403.-0.63331E-01 -492.06 0.0000

0.0000 0.11573E-02 -0.63331E-01 8533.5 0.41546E-03 0.22257E-04

-0.11573E-02 0.0000 -492.06 0.41546E-03 8536.0-0.13031E-05

0.63331E-01 492.06 0.0000 0.22257E-04 -0.13031E-05 17066.

GENERAL DYNA.STIFFNESS MATRIX OF JOINT:

KJ1[N/μm]= 0.103E+06

KJ2[N/μm]= 0.103E+06

KJ3[N/μm]= 0.215E+05

KJ4[N·m/(μm/m)]= 0.724E+04

KJ5[N·m/(μm/m)]= 0.724E+04

KJ6[N·m/(μm/m)]= 0.145E+05

GENERAL DYNA.DAMPING MATRIX OF JOINT:

CJ1[N·s/m/mm2]=138.

CJ2[N·s/m/mm2]=138.

CJ3[N·s/m/mm2]=18.7

CJ4[N·m·s/rad]=6.29

CJ5[N·m·s/rad]=6.29

CJ6[N·m·s/rad]=12.6

2) 设圆柱面的半径为R'，则球台面等效为圆柱面的半径近似为：

R'=R·sinφa=0.819875m , 而且Za=R·cosφa=0.014311m，

Zb=R·cosφb=-0.014311m。通过圆柱结合部程序计算输出结果为：

STATIC STIFFNESS MATRIX OF JOINT SURFACE:

IJS=1

0.12449E+06 -0.12315E-02 0.00000.25746E-07

-0.41531E-05 0.61525E-01

-0.12315E-02 0.12449E+06 0.0000-0.19629E-05 -0.25746E-07 491.78

0.0000 0.0000 25378.-0.61525E-01 -491.78 0.0000

0.25746E-07 -0.19629E-05 -0.61525E-01 8534.5 0.53026E-03 0.37521E-06

-0.41531E-05 -0.25746E-07 -491.78 0.53026E-03

8537.0 -0.36427E-06

0.61525E-01 491.78 0.0000 0.37521E-06 -0.36427E-06 17059.

GENERAL DYNA.STIFFNESS MATRIX OF JOINT:

KJ1[N/μm]= 0.103E+06

KJ2[N/μm]= 0.103E+06

KJ3[N/μm]= 0.215E+05

KJ4[N·m/(μm/m)]= 0.724E+04

KJ5[N·m/(μm/m)]= 0.724E+04

KJ6[N·m/(μm/m)]= 0.145E+05

GENERAL DYNA.DAMPING MATRIX OF JOINT:

CJ1[N·s/m/mm2]=138.

CJ2[N·s/m/mm2]=138.

CJ3[N·s/m/mm2]=18.7

CJ4[N·m·s/rad]=6.30

CJ5[N·m·s/rad]=6.30

CJ6[N·m·s/rad]=12.6

b) 实例验证结论

由此可以看出，当球面半径为R时，其中φ=89°-91°，可以把球台面近似为半径R'=R·sinφ的圆柱面求解，通过对圆柱面结合部计算的对比就可以论证本文中球面结合部理论解析及解析程序是正确的。

3 总结

通过建立球面结合部数学模型，然后推导出球面结合部位移及载荷结合条件，并给出了球面结合部刚度矩阵的求法，为以后机械系统性能分析提供通用性参数打下基础。最后介绍了编制基于FORTRAN语言的球面结合部解析程序说明，并通过实例对比验证了球面结合部解析方法及解析程序是正确的。

[1] 黄玉美,张广鹏,高峰.虚拟样机整机结构特性边界元仿真[M].北京：机械工业出版社,2004:152-173.

[2] 董献国.导结合部的变形解析、加工误差综合及其结构优化设计[D].西安：陕西机械学院硕士论文.1989.

[3] 王晓春,孔祥安.接触力学及其计算方法[J].西南交通大学报,1996,31(3):230-233.

[4] 温卫东,高德平.接触问题数值分析方法的研究现状与发展[J].南京航空航天大学学报, 1994，26(5):664-667.

[5] Parson B,Wilson E A.A method for determning the surface contact stresses resulting from interference fits[J].ASME J Engng for Industry,1970;208-218.

[6] Ohtakake K,Oden J T,Kikuchi N.Analysis of certain unilateral problem in von karman plate theory by a Penalty method-Part I:a variational pririple with penalty[J].ComPut & Meth Appl Mech and Engng,1980;24:187-213.