等节距截锥涡卷螺旋弹簧

等节距截锥涡卷螺旋弹簧静力
及模态分析
姓名：
学号：
摘要
变刚度弹簧以其独特的刚度特性在工业各领域内得到广泛应用，截锥涡卷弹簧式这种非线性螺旋弹簧的典型代表。

本文建立了截锥涡卷弹簧（特殊的等节距螺旋弹簧，其中心展开线为抛物线，螺旋线在底面投影为阿基米德螺旋线）的有限元模型，并在弹簧顶端利用mass 单元等效施加一质量块，计算质量弹簧系统在压并前的自振频率以及前三阶模态。

另外，在弹簧的静力分析中，通过子步指定荷载逐渐增加，分析弹簧的变形过程及其应力变化。

建模和分析采用APDL（ANSYS PARAMETER DESIGN LANGUAGE)，所以结构尺寸和材料很容易根据特定设计要求更改，进而为截锥涡卷弹簧的设计提供了数值分析模拟的一定参考。

关键词:有限元仿真, 等节距截锥螺旋压缩弹簧, 大变形，非线性，模态分析，接触分析，变刚度
一、背景
1、弹簧
弹簧是一种利用弹性来工作的机械零件。

一般用弹簧钢制成。

用以控制机件的运动、缓和冲击或震动、贮蓄能量、测量力的大小等，广泛用于机器、仪表中。

弹簧的种类复杂多样，按形状分，主要有螺旋弹簧、涡卷弹簧、板弹簧等。

其主要功能①控制机械的运动，如内燃机中的阀门弹簧、离合器中的控制弹簧等。

②吸收振动和冲击能量，如汽车、火车车厢下的缓冲弹簧、联轴器中的吸振弹簧等。

③储存及输出能量作为动力，如钟表弹簧、枪械中的弹簧等。

④用作测力元件，如测力器、弹簧秤中的弹簧等。

弹簧的载荷与变形之比称为弹簧刚度，刚度越大，则弹簧越硬。

过去的机械部件等曾经由于没有给予弹簧足够的重视而经常损坏，造成极为严重的后果，随着科学技术的进步，人们逐渐认识到弹簧对机器的精度，工作能力和寿命都有着极其重要的作用。

现代的数值仿真方法为弹簧的设计提供了一个新的平台，可以模拟各种弹簧在各种工况下的静力状态和动力学响应。

2、等节距截锥螺旋压缩弹簧
变刚度弹簧以其独特的刚度特性在工业各领域内得到广泛应用，截锥涡卷弹簧式这种非线性螺旋弹簧的典型代表。

截锥涡卷弹簧又名宝塔弹簧，主要特点是体积小、载荷大、变刚度，广泛用于空间小、载荷大的场合和减震装置，它可分为等螺旋角、等节距、等应力三种
形式。

截锥涡卷弹簧是由加工成扁薄的弹簧钢板热卷成。

主要应用于外形尺寸受限制而又需要吸收巨大冲击能量的机械上。

它的刚性是变化的，在刚接触负载时刚性较低继而转强．这一特性正好缓和了物体碰撞瞬间所产生的强烈冲击作用，因而广泛应用于某些重型机械上因工作荷载范围大、振动冲击较强、安装空间又受限制的减震装置上。

这种弹簧在受力后特性线段式直线，当荷载增大时弹簧的大圈开始接触，有效圈数减少而刚度增大，一直到所有的弹簧圈完全压并为止，这种弹簧的刚度是变值，所以自振频率也是变值，有利于防止共振，因而多用于需要减震的场合。

对于特殊的等节距螺旋弹簧，其中心展开线为抛物线，螺旋线在底面投影为阿基米德螺旋线（见下文ANSYS建模过程）。

3、阿基米德螺旋线
阿基米德螺旋线，亦称“等速螺线”。

当一点P沿动射线OP一等速率运动的同时，这射线有以等角速度绕点O旋转，点P的轨迹称为“阿基米德螺线。

它的极坐标方程为：r = a θ。

阿基米德弹簧，顾名思义就是弹簧中心线为空间阿基米德曲线的弹簧。

本文通过理论分析并以ANSYS为技术支持平台，讨论了等节距截锥涡卷螺旋弹簧的刚度特性以及动力特性。

二、理论分析
本文中用到的主要变量及其物理意义和具体数值如下表：
表1 变量表
变量符号物理意义值t弹簧节距9.6mm
r弹簧最内层半径14mm
1
r弹簧最外层半径43mm
2
n弹簧有效圈数 5
a弹簧板厚度4mm
b弹簧板宽度28mm
F重物重力1450N
1、弹簧刚度求解
推导截锥螺旋弹簧刚度计算近似公式，并计算弹簧的刚度曲线和质量弹簧系统在压并前的自振频率，将结果与有限元结果对比
（1）截锥涡卷螺旋弹簧及其特性线
未有弹簧圈接触的变形与强度计算：
由于弹簧的螺旋角α比较小，当弹簧受到轴向荷载F作用时，可以近似看做各弹
簧圈材料截面只受到如下扭矩：
Tt ≈F R
在单位轴向载荷作用时，材料截面所受的扭矩：
T 1t ≈1⋅R
弹簧轴向变形的一般积分式
f =∫TtT 1tds /GIp =∫FR 2/GIpds
由 ds=Rd θ, d θ=2πnR 2−R 1dR 得 d s=2πnR 2−R 1RdR 代入上式可得变形计算式：
)/(2*)()(2124
14
22
1
123
R R GI R R F n dR GI R R nFR t t R R --=-⎰ππ
对于矩形截面材料，由于b/a=7，查表则相当极惯性矩近似取为It ≈0.3ba 3 代入得： f =πnF (R 42−R 41)/2G *0.3ba 3(R 2−R 1) 则弹簧的刚度计算式为：
K =2Gba 3(R 2−R 1)*0.3/πn (R 42−R 41) （1）
有弹簧圈接触后的变形与强度计算： 从大端数起到弹簧圈i 的自由高度
Hi =nt (R 2−RiR 2−R 1)+b
弹簧圈i 接触时所受的载荷
Fi =Gba 3t 6πRi 3 （2）
载荷Fi 与变形fi 的关系
fi =nπR 2−R 1[(R 2−Ri )tπ+3(Ri 4−R 142Gba 3)Fi ] （3）
120.014,0.043,0.04,0.026,0.0096,5R m R m a m b m t m n ======
错误！未找到引用源。

第一圈半径： 12211
()0.043(0.0430.014)/50.0372R R R R m n
=--⨯=--= 错误！未找到引用源。

第一圈接触底面时所施加的力：
11
3312100.0260.0040.0096/(60.0372)1266.342(10.3)
F N π⨯=⨯⨯⨯⨯=⨯+
第二圈半径： '
222120.0430.014()0.04320.03145
R R R R m n -=--⨯=-⨯=
第二圈接触底面时所施加的力：
11
3322100.0260.0040.0096/(6*0.0314)2105.672(10.3)
F N π⨯=⨯⨯⨯=⨯+
第三圈半径： '
322130.0430.014()0.04330.02565
R R R R m n -=--⨯=-⨯= 错误！未找到引用源。

第三圈接触底面时所施加的力：
11
3332100.0260.0040.0096/(60.0256)3885.622(10.3)
F N π⨯=⨯⨯⨯⨯=⨯+
第四圈半径：'
4221()4/0.043(0.0430.014)4/50.0198R R R R n m =--⨯=--⨯=
错误！未找到引用源。

第四圈接触底面时所施加的力：
11
3342100.0260.0040.0096/(60.0198)8398.172(10.3)
F N π⨯=⨯⨯⨯⨯=⨯+
第五圈半径：'
5221()5/0.043(0.0430.014)/550.014R R R R n m =--⨯=--⨯=
第五圈接触底面式所施加的力（即完全压并时所施加的力）：
113352100.0260.0040.0096/(60.014)23757.242(10.3)
F N π⨯=⨯⨯⨯⨯=⨯+
三、ANSYS 数值模拟
1、建模
弹簧的材料参数为：
E ——弹性模量，2e11 Pa λ——泊松比，0.3
ρ——材料密度，7800kg/m3
设计模型中的各物理参数：
a、b——截锥涡卷螺旋弹簧弹簧板材料厚度和宽度
R2——为弹簧最外层的半径R1——为弹簧最内层的半径
t——为弹簧节距n——为弹簧的有效圈数
nz——弹簧支撑圈数F——弹簧所受外荷载（重物重力）
等节距截锥涡卷螺旋弹簧，各参数详见表1
表1 弹簧参数
R1(mm) R2(mm) a、b(mm) n nz t(mm) F(N)
14 43 a=4,b=26 5 0.75+0.75 9.6 1450
首先利用循环语句取出空间阿基米德曲线上的关键点，然后在柱坐标下利用line语句光滑连接各个关键点形成空间曲线。

由于设计要求在最上方和最下方有3/4段弹簧，以用来连接其他构件，为避免曲线的不光滑，连接处的曲线仍为阿基米德螺旋线，而连接段的高度坐标不变。

为了减少实体模型的面积的个数，以使后面的加载工作更加方便，利用LCOMB 语句，合并各条曲线。

这样空间阿基米德曲线便建立完毕了。

在曲线的一端建立弹簧实体的截面，通过拖拽操作，将截面沿着空间阿基米德曲线拖拽形成实体模型。

然后利用布尔操作，将各个体积合成。

截面尺寸高度为28mm，宽为4mm，弹簧总高度为48mm，共5圈，加上连接圈（1.5圈）为6.5圈。

建立的阿基米德螺旋线如下图所示。

空间阿基米德螺旋线
实体模型
2、网格划分
将截面的线设置分割数为2和4，然后进行扫略操作。

利用SOLID45八节点六面体单元将实体模型剖分，获得有限元模型，如下图。

有限元模型
3、边界条件问题分析
（1）、模态分析
模态分析中，对于等节距截锥涡卷螺旋弹簧，边界条件为约束使弹簧最下面的连接面A157的X、Y、Z方向位移。

在上3/4圈连接处将某单元（E5517）转换为MASS21质量单元。

(2)、静力分析
静力分析中，等节距截锥涡卷螺旋弹簧在工作时，当荷载增大的一定程度，弹簧的大圈逐渐接触，有效工作圈数减少，刚度逐渐增大，一直到所有的弹簧圈完全压并，刚度理论上可达到与材料弹性模量相当的值。

本问题中，大变形和接触问题是导
致问题非线性的两个因素。

其中接触问
题是高度非线性的而且在求解过程中需
要占用大量的计算机资源，因此对接触
的设置应当慎重。

弹簧被“压平”是采
用接触分析。

在弹簧底部创建一平板（厚
度为6mm）然后利用DOF命令约束底面
接触面位移和平板的上下两面的个各方
向位移。

力的大小24210N等效为面力
118583293.2N/M2施加于上接触圈（有
限元中约为17/25圈）的上表面。

目标面
为底板上表面，接触面为弹簧的各下底
面。

接触分析如下图所示。

4、加载求解
模态分析中定义分析类型为模态分析，ANTYPE，MODAL。

求前3阶模态。

静力分析涉及到大变形非线性问题，所以在求解时将大变形开关打开，NLGEOM，ON。

设置求解类型为静态求解，ANTYPE，0。

设置时间为1s，加载的荷载子步为20步。

相关荷载大小等见前述边界条件问题分析。

5、仿真结果分析及讨论
模态分析求解固有频率
在对阿基米德弹簧进行模态分析过程中，计算前三阶模态。

利用PLDISP，2命令显示模态和弹簧原始形状。

并且利用ANMODE,10,0.15,,0语句形成动画显示，从而可以清晰而直观的了解弹簧的前三阶模态，如下图所示
第一阶自振频率：FREQ1 = 2.8284（z方向的振动）
第二阶自振频率：FREQ2 =52.396 （y方向的振动）
第三阶自振频率：FREQ3 =4.798
利用*GE命令求出弹簧前三阶模态对应的固有频率，弹簧的前三阶固有频率经过计算如下：
FREQ1 2.8284 SCALAR
FREQ2 52.396 SCALAR
FREQ3 54.798 SCALAR 前三阶模态描述：通过动画显示，可以看到，第一阶模态为沿着X方向弹簧的横向震动，第二阶模态为大概沿着Z方向的前后震动。

第三阶模态为沿着XOZ平面二四象限方向的横向震动。

经过对系统进行模态分析，有限元模拟得到，系统固有频率，即一阶主频为2.72Hz。

与理论解2.80Hz相比，相对误差P=2.8%，误差是比较小的，说明之前计算固有频率的方法也是可行的，而且非常简单，同时也验证了数值模拟的准确性。

静态分析求弹簧特性曲线
从弹簧刚度位移曲线可以看出，再受力最初阶段，弹簧处于线性阶段，之后随着加载刚度逐渐增大。

上图基本上能说明等节距截锥涡卷螺旋弹簧的受力特性。

由变形图知，第一荷载步1210N作用下，此时刚度为
K=1210/0.024912=48571N/mm
与理论值的相对误差为5.9%，误差较小。

误差主要来源于第一荷载步1210N作用下可能有部分弹簧底部圈接触致使刚度变大，即若第一荷载步作用力若较1210小些，结果会更精准。

与理论值的相对误差为3.8%,误差较小。

说明之前的理论解法是可行的，而且也进一步验证了数值模拟的可信度。

5、结论
本文利用ANSYS对等节距截锥涡卷螺旋弹簧进行了模态分析和静力分析，求得了弹簧的前三阶模态和固有频率和静力分析。

解决了以下几个问题：
1. 绘制弹簧的特性曲线（荷载-变形曲线），并给出弹簧的刚度曲线图（刚度-变形曲线）；
2. 在弹簧顶端施加一质量块，计算质量弹簧系统在压并前的自振频率。

（认为质量块与弹簧连接完好）；
3. 手算弹簧刚度及质量弹簧系统在压并前的自振频率，将结果与有限元结果对比；
4. 在实际应用过程中存在质量块偏心压弹簧的情况，如何防止偏心。

参考文献
1. 《ANSYS 1
2.0从入门到实战》何波机械工业出版社,张洪才等
2.ANSYS7.0分析实例与工程应用.机械工业出版社. 邢静忠,王永刚,陈晓霞等
3.弹簧设计手册。

上海科学技术文献出版社.汪曾祥,魏先英,刘祥至
4.ANSYS参数化编程与命令手册.机械工业出版社.龚曙光,谢桂兰,黄云清
5. 变刚度螺旋弹簧的设计方法和精确建模初探.中国机械工程第13卷.第13其.祁宏钟,雷雨
成,冯晋祥
6. 大螺距圆锥形弹簧受横向力时的挠度刘东旭，刘志刚现代机械1995年6月第3期
7. 等螺旋升角截锥螺旋弹簧朱炎，朱琪，郑方芳2010年全国弹簧学术会论文集
8.机车车辆截锥螺旋压缩弹簧减震设计与研究凌凯夫，刘建秀，毛蔚，王军湖南工业大学
学报，VOL.23 NO.4
9.基于ANSYS的变刚度螺旋弹簧优化设计石丽娜，曾红，刘淑芬辽宁工业大学学报
VOL29,NO.1,Feb，2009
10. 在ANSYS中实现圆柱螺旋弹簧的参数化建模刘晓波，沈予洪煤矿机械2005年第4
命令流
模态分析：
FINI
/CLE
/FILNAME, Mode Analysis of Archimedes Spring
/PREP7
/TITLE, Mode Analysis of Archimedes Spring
ET,1,PLANE182
ET,2,SOLID45
MP,EX,1,2E11
MP,NUXY,1,0.3
MP,DENS,1,7800
*do,i,0,120
*set,r,0.043-0.0060416*i/25
*set,x,r*cos(i*3.14159/12)
*set,y,0.01*i/25
*set,z,r*sin(i*3.14159/12)
k,i+1,x,y,z
*enddo
CSYS,5
*do,i,1,120
L,i,i+1
*enddo
CSYS,0
*do,i,-1,-18,-1
*set,r,0.043-0.0060416*i/25
*set,x,r*cos(i*3.14159/12)
*set,y,0
*set,z,r*sin(i*3.14159/12)
k,121-i,x,y,z
*enddo
*do,i,1,17
*set,r,0.014-0.0060416*i/25
*set,x,r*cos(120*3.14159/12+i*3.14159/12) *set,y,0.048
*set,z,r*sin(120*3.14159/12+i*3.14159/12) k,139+i,x,y,z
*enddo
CSYS,5
L,1,122
*DO,I,122,138
L,I,I+1
*ENDDO
L,121,140
*DO,I,140,155
L,I,I+1
*ENDDO
*DO,I,1,60,1
LCOMB,2*I-1,2*I,0
*ENDDO
*DO,I,61,69,1
LCOMB,2*I-1,2*I,0
*ENDDO
*DO,I,70,77,1
LCOMB,2*I-1,2*I,0
*ENDDO
*DO,I,1,30,1
LCOMB,4*I-3,4*I-1,0
*ENDDO
*DO,I,123,137,2
LCOMB,121,I,0
*ENDDO
*DO,I,141,155,2
LCOMB,139,I,0
*ENDDO
kwpave,1
rectng,-0.002,0.002,-0.014,0.014 vdrag, 1 ,,,,,, 1
vdrag, 6 ,,,,,, 5
vdrag, 11 ,,,,,, 9
vdrag, 16 ,,,,,, 13 vdrag, 21 ,,,,,, 17 vdrag, 26 ,,,,,, 21 vdrag, 31 ,,,,,, 25 vdrag, 36 ,,,,,, 29 vdrag, 41 ,,,,,, 33 vdrag, 46 ,,,,,, 37 vdrag, 51 ,,,,,, 41 vdrag, 56 ,,,,,, 45 vdrag, 61 ,,,,,, 49 vdrag, 66 ,,,,,, 53 vdrag, 71 ,,,,,, 57 vdrag, 76 ,,,,,, 61 vdrag, 81 ,,,,,, 65 vdrag, 86 ,,,,,, 69 vdrag, 91 ,,,,,, 73 vdrag, 96 ,,,,,, 77 vdrag, 101 ,,,,,, 81 vdrag, 106 ,,,,,, 85 vdrag, 111 ,,,,,, 89 vdrag, 116 ,,,,,, 93 vdrag, 121 ,,,,,, 97 vdrag, 126 ,,,,,, 101 vdrag, 131 ,,,,,, 105 vdrag, 136 ,,,,,, 109 vdrag, 141 ,,,,,, 113 vdrag, 146 ,,,,,, 117 vdrag, 151 ,,,,,, 139 vdrag, 1 ,,,,,, 121
V ADD,ALL
CSYS,0
LESIZE,290,,,2
LESIZE,288,,,4
LESIZE,285,,,2
LESIZE,292,,,4
TYPE,2
MA T,1
vsweep,33
FINISH
/PREP7
ET,3,MASS21
R,1, ,148, , , , ,
TYPE, 3
MA T, 1
REAL, 1
ESYS, 0
SECNUM,
TSHAP,LINE
E, 5517
/SOLU
DA,157,ALL
ANTYPE,MODAL
MODOPT,SUBSP,3 ！计算前三阶模态
SOLVE
/POST1
*GET,FREQ1,MODE,1,FREQ
*GET,FREQ2,MODE,2,FREQ
*GET,FREQ3,MODE,3,FREQ
*STATUS,PARAMETERS
SET,FIRST ！从结果文件中读出第一个荷载步的数据
PLDISP,2
SET,NEXT
PLDISP,2
SET,1,3
PLDISP,2
ANMODE,10,0.15,,0
FINISH
1个MASS21单元在上面3/4连接处
FREQ1 2.72078500 SCALAR FREQ2 51.1174469 SCALAR FREQ3 52.6870141 SCALAR 一个MASS21单元在最顶端的顶点处
/PREP7
ET,3,MASS21
R,1, ,148, , , , ,
TYPE, 3
MA T, 1
REAL, 1
ESYS, 0
SECNUM,
TSHAP,LINE
E, 5516
FREQ1 2.70705719 SCALAR FREQ2 51.0119281 SCALAR FREQ3 53.6286777 SCALAR I 155.000000 SCALAR R 9.891712000E-03 SCALAR X -2.560452913E-03 SCALAR Y 4.800000000E-02 SCALAR Z -9.554582522E-03 SCALAR 静力分析：
FINI
/CLE
/FILNAME, Static And Mode Analysis of Archimedes Spring
/PREP7
/TITLE, Static And Mode Analysis of Archimedes Spring
ET,1,PLANE182
ET,2,SOLID45
MP,EX,1,2E11
MP,NUXY,1,0.3
MP,DENS,1,7800
*do,i,0,120
*set,r,0.043-0.0060416*i/25
*set,x,r*cos(i*3.14159/12)
*set,y,0.01*i/25
*set,z,r*sin(i*3.14159/12)
k,i+1,x,y,z
*enddo
CSYS,5
*do,i,1,120
L,i,i+1
*enddo
CSYS,0
*do,i,-1,-18,-1
*set,r,0.043-0.0060416*i/25
*set,x,r*cos(i*3.14159/12)
*set,y,0
*set,z,r*sin(i*3.14159/12)
k,121-i,x,y,z
*enddo
*do,i,1,17
*set,r,0.014-0.0060416*i/25
*set,x,r*cos(120*3.14159/12+i*3.14159/12) *set,y,0.048
*set,z,r*sin(120*3.14159/12+i*3.14159/12) k,139+i,x,y,z
*enddo
CSYS,5
L,1,122
*DO,I,122,138
L,I,I+1
*ENDDO
L,121,140
*DO,I,140,155
L,I,I+1
*ENDDO
*DO,I,1,60,1
LCOMB,2*I-1,2*I,0
*ENDDO
*DO,I,61,69,1
LCOMB,2*I-1,2*I,0
*ENDDO
*DO,I,70,77,1
LCOMB,2*I-1,2*I,0
*ENDDO
*DO,I,1,30,1
LCOMB,4*I-3,4*I-1,0
*ENDDO
*DO,I,123,137,2
LCOMB,121,I,0
*ENDDO
*DO,I,141,155,2
LCOMB,139,I,0
*ENDDO
kwpave,1
rectng,-0.002,0.002,-0.014,0.014
vdrag, 1 ,,,,,, 1
vdrag, 6 ,,,,,, 5
vdrag, 11 ,,,,,, 9
vdrag, 16 ,,,,,, 13
vdrag, 21 ,,,,,, 17
vdrag, 26 ,,,,,, 21
vdrag, 31 ,,,,,, 25
vdrag, 36 ,,,,,, 29
vdrag, 41 ,,,,,, 33
vdrag, 46 ,,,,,, 37
vdrag, 51 ,,,,,, 41
vdrag, 56 ,,,,,, 45
vdrag, 61 ,,,,,, 49
vdrag, 66 ,,,,,, 53
vdrag, 71 ,,,,,, 57
vdrag, 76 ,,,,,, 61
vdrag, 81 ,,,,,, 65
vdrag, 86 ,,,,,, 69
vdrag, 91 ,,,,,, 73
vdrag, 96 ,,,,,, 77
vdrag, 101 ,,,,,, 81
vdrag, 106 ,,,,,, 85
vdrag, 111 ,,,,,, 89
vdrag, 116 ,,,,,, 93
vdrag, 121 ,,,,,, 97
vdrag, 126 ,,,,,, 101
vdrag, 131 ,,,,,, 105
vdrag, 136 ,,,,,, 109
vdrag, 141 ,,,,,, 113
vdrag, 146 ,,,,,, 117
vdrag, 151 ,,,,,, 139
vdrag, 1 ,,,,,, 121
V ADD,ALL
WPA VE,0,0,0
wpro,,-90,
wpoff,0,0,-0.014
RECTNG,-0.050,0.050,-0.050,0.050 wpro,,90,
VEXT,1, , ,0,0,-0.06,,,,
CSYS,0
LESIZE,290,,,2
LESIZE,288,,,4
LESIZE,285,,,2
LESIZE,292,,,4
TYPE,2
MA T,1
VSWEEP,33
TYPE,2
MA T,1
esize,10,
vMESH,1
FINISH
/PREP7
ET,3,TARGE170
ET,4,CONTA174
/COM, CONTACT PAIR CREATION - START CM,_NODECM,NODE
CM,_ELEMCM,ELEM
CM,_KPCM,KP
CM,_LINECM,LINE
CM,_AREACM,AREA
CM,_VOLUCM,VOLU
/GSA V,cwz,gsav,,temp
MP,MU,1,0.3
MA T,1
MP,EMIS,1,7.88860905221e-031
R,3
REAL,3
ET,3,170
ET,4,174
R,3,,,0.1,0.1,0,
RMORE,,,1.0E20,0.0,1.0,
RMORE,0.0,0,1.0,,1.0,0.5
RMORE,0,1.0,1.0,0.0,,1.0
KEYOPT,4,4,0
KEYOPT,4,5,0
NROPT,UNSYM KEYOPT,4,7,0 KEYOPT,4,8,0 KEYOPT,4,9,0 KEYOPT,4,10,2 KEYOPT,4,11,0 KEYOPT,4,12,0 KEYOPT,4,2,0 KEYOPT,3,5,0
! Generate the target surface
ASEL,S,,,1
ASEL,A,,,1
CM,_TARGET,AREA TYPE,3
NSLA,S,1
ESLN,S,0
ESURF
CMSEL,S,_ELEMCM
! Generate the contact surface
*DO,I,2,157,5
ASEL,S,,,1
ASEL,A,,,I
CM,_CONTACT,AREA TYPE,4
NSLA,S,1
ESLN,S,0
ESURF
*ENDDO
ALLSEL
ESEL,ALL
ESEL,S,TYPE,,3
ESEL,A,TYPE,,4
ESEL,R,REAL,,3
/PSYMB,ESYS,1
/PNUM,TYPE,1
/NUM,1
EPLOT
ESEL,ALL
ESEL,S,TYPE,,3
ESEL,A,TYPE,,4
ESEL,R,REAL,,3
CMSEL,A,_NODECM
CMDEL,_NODECM
CMSEL,A,_ELEMCM
CMDEL,_ELEMCM
CMSEL,S,_KPCM
CMDEL,_KPCM
CMSEL,S,_LINECM
CMDEL,_LINECM
CMSEL,S,_AREACM
CMDEL,_AREACM
CMSEL,S,_VOLUCM
CMDEL,_VOLUCM
/GRES,cwz,gsav
CMDEL,_TARGET
CMDEL,_CONTACT
/COM, CONTACT PAIR CREATION - END FINISH
/SOLU
DA,157,ALL
DA,1,ALL
DA,6,ALL
SFA,154,1,PRES,118583293.2
/SOLU
ALLSEL,ALL
TIME,1
ANTYPE,0
NSUB,10
KBC,0 ！
OUTRES,ERASE
OUTRES,ALL,ALL
SOLVE
FINISH。