弹簧单元combin39如何实现单向受拉或者单向受压

        弹簧单元一直是ANSYS中使用频率较高的单元，尤其是非线性弹簧单元combin39, 39单元功能较多，单元选项设置复杂，在很多方面都有其独特的运用，今天水哥介绍下39单元比较特殊的一面，也即单元的单向行为。

        所谓单向行为，也即单元的单向受压或者单向受拉，在土木地下结构计算领域，单向弹簧的使用频率最高，常常用于其模拟弹性地基，等效围岩等。

        要利用该单元实现单向弹簧，首先要读懂该单元各个单元关键项的意思，该单元有很多关键项，不同的设置会有不同的单元表现。该单元一共有八种单元表现，罗列如下：

         从上述单元表现可见，第B种和第e种情况可实现单向弹簧的功能，这两者的主要区别在于一个是卸载路径与原加载路径相同，一种是卸载路径与加载路劲的原点段平行。

        细心的同学可以发现，这儿combin39所谓的单向是指受拉单向，也即是该单元只提供单向受拉的功能，如果要实现我们口中所谓的单向受压，则需要一定的建模技巧。

        为验证该单元的单向功能，下面我们做一个小实验。

        命令流如下：

finish

/clear

/prep7

et,1,combin39

!Z方向的单向弹簧

keyopt,1,4,0

keyopt,1,3,3

keyopt,1,1,0

keyopt,1,2,1

n,1

n,2,0,0,1.0

!弹簧的初始弹性模量为100

r,1,0.1,100*0.1

e,1,2

d,1,all,0

allsel,all

!弹簧受到拉力100，可以预见最大位移为1

f,2,fz,100

/solu

allsel,all

nsubst,10,,1

autots,on

solve

        结果位移云图如下：

在上述命令流下，我们将拉力改为压力，可以预见，由于弹簧只是受拉弹簧，因此结果应该为0.

        修改的命令流如下：

!==============

f,2,fz,-100

!=============

        修改之后求解，虽然软件显示的是solution is done，然而并没有位移结果，这也表明，该弹簧不能模拟受压。

        从这个位移云图可见，弹簧被压缩了1，而这与我们口中所说的受压情况很符合。

        综合这个小实验，我们可以得到一些有用的结论：

        1、Combin39只能模拟单向受拉.

        2、Combin39所谓的拉应变是指沿着单元坐标系X的正方向，如果我们需要模拟单向受压，只需调整该单元的单元坐标系。

        为体现该单元使用上的一些技能，设有某矩形大板，长80m，宽40m，在中间3m范围内作用有均布荷载，地基刚度为20MPa,采用ANSYS模拟该过程。

个可以预见的结果便是 由于该板相对于加载面积来讲，其长边太长，会有一种中间凹两边凸的直观效果，因此如果采用弹簧模拟弹性地基，两边的弹簧会处于受拉状态，然而事实是我们只希望地基弹簧只受压不受拉，受拉弹簧退出工作。

        我们采用combin39来模拟该过程。命令流如下：

!===========

finish

/clear

/prep7

A=80$B=40$H=0.6

A1=3$B1=2

Q=150e3

ESF=2.0e7            !弹性地基刚度

ET,1,solid186

et,2,surf154         !辅助表面效应单元

et,3,combin39

keyopt,3,4,0

keyopt,3,3,3

keyopt,3,1,0

keyopt,3,2,1

!建模

!===========

!===================

allsel,all

*get,tmmax,node,,num,max  !整个模型的最大节点号

esel,s,type,,2

nsle,s,corner

*get,ntol,node,,count !节点总数

*get,etol,elem,,count !单元总数

*get,nmax,node,,num,max !最大节点号

*dim,nodno,,ntol !存储节点号

*dim,nodk,,nmax  !存储各节点面积、弹簧刚度

*dim,eleno,,etol !存储单元号

!获得各个单元号

*get,e1,elem,,num,min

eleno(1)=e1

*do,i,2,etol

e1=elnext(e1)

eleno(i)=e1

*enddo

!获得各节点号

*get,n1,node,,num,min

nodno(1)=n1

*do,i,2,ntol

n1=ndnext(n1)

nodno(i)=n1

*enddo

!求每个节点对应的面积

*do,i,1,etol

ei=eleno(i)

n1=nelem(ei,1)

n2=nelem(ei,2)

n3=nelem(ei,3)

n4=nelem(ei,4)

*get,ai,elem,ei,area

*if,n3,ne,n4,then

ai=ai/4.0

nodk(n1)=nodk(n1)+ai

nodk(n2)=nodk(n2)+ai

nodk(n3)=nodk(n3)+ai

nodk(n4)=nodk(n4)+ai

*else

ai=ai/3.0

nodk(n1)=nodk(n1)+ai

nodk(n2)=nodk(n2)+ai

nodk(n3)=nodk(n3)+ai

*endif

*enddo

!求得弹簧刚度

*do,i,1,nmax

nodk(i)=nodk(i)*ESF

*enddo

*do,i,1,ntol

ni=nodno(i)

r,i+1,0.01,nodk(ni)*0.01

*enddo

!创建弹簧单元

type,3

*do,i,1,ntol

ni=nodno(i)

n,tmmax+i,nx(ni),ny(ni),nz(ni)-0.5 !创建节点(重合)

real,i+1

!注意此处生成弹簧单元节点顺序，改变了顺序也即是改变了单元坐标系X的方向。

e,ni,tmmax+i

*enddo

nsel,s,,,tmmax+1,tmmax+ntol

d,all,all

finish

/solu

allsel,all

nsubst,100,,1

autots,on

solve

        最后的竖向位移云图如下：

可见其位移趋势符合常识。弹簧单元的轴力云图如下：

        从图中可见，只有中部弹簧受力，两边受拉弹簧轴力为0，说明该过程有效的实现了单向弹簧作用。

        祝好

        ANSYS结构院

        2018.01.18