!060115
!通過節點反力循環求和方法以及FSUM方法2種方式計算軸力
!先計算過盈配合。
!然后抽出試樣,分析摩擦。
!不考慮熱膨脹
FINI
/CLE
/FILE,COLLET
/TITLE,COLLET ANALYSIS
/CONFIG,NRES,10000		!設置最大子步數限制
*AFUN,DEG			!使用角度為單位
*SET,THETA,3			!定義錐度
*SET,LONG,70			!工作錐長度
*SET,DIN,9.5/2			!末端內半徑 與試樣相同
*SET,DOUT,22/2			!末端外半徑11 
*SET,LL,60			!最大夾持長度(深度)
*SET,LM,25			!套筒長度
*SET,DD,5			!套筒離末端初始長度
*SET,DV,30/2			!套筒末端外半徑
*SET,UU,0.3			!套筒向右位移量
*SET,XX,30			!試樣夾持長度
*SET,LS,80			!試樣右端到原點距離
*SET,NN,9			!設置錐的縫隙角度(一般為3,5,9度)
*SET,FKNN,1			!設置剛度因子
*SET,STEP1,20			!設置第一載荷步子步數
*SET,STEP2,50			!設置第二載荷步子步數
!前處理
/PREP7
ET,1,MESH200			!平面模型單元
KEYOPT,1,1,6			!四邊形單元
ET,2,SOLID185			!實體模型單元
ET,3,CONTA173,1			!3D四節點接觸單元
KEYOPT,3,1,0			!設定自由度
KEYOPT,3,12,0			!3號單元接觸面行為
NROPT,UNSYM			!指定非對稱接觸剛度
ET,4,TARGE170			!目標單元
!材料1為試樣
MP,EX,1,1.9E5			!材料1彈性模量
MP,NUXY,1,0.3			!材料1泊松比
MP,ALPX,1,1.7E-5		!材料1熱膨脹系數
MP,MU,1,0.2			!材料1摩擦系數
!材料2為錐體
MP,EX,2,2.05E5			!材料2彈性模量		
MP,NUXY,2,0.3			!材料2泊松比
MP,ALPX,2,1.8E-5		!材料2熱膨脹系數
MP,MU,2,0.2			!材料2摩擦系數
!材料3為套筒
MP,EX,3,1.95E5			!材料3彈性模量					
MP,NUXY,3,0.3			!材料3泊松比
MP,ALPX,3,1.0E-5		!材料3熱膨脹系數
MP,MU,3,0.2			!材料3摩擦系數
R,10,,,FKNN,0.1			!接觸對10,接觸剛度FKNN,罰系數0.1
R,11,,,FKNN,0.1			!接觸對11,接觸剛度FKNN,罰系數0.1
/PNUM,KP,1			!打開關鍵點編號開關
/PNUM,LINE,1			!打開線編號開關
/PNUM,AREA,1			!打開面編號開關
/VIEW,1,1,1,1			!正等軸測顯示
!創建模型
!創建關鍵點
K,1				!模型左端界限
K,2,LS				!模型右端界限
K,3,,DIN			!左端錐體圓柱控制點
K,4,LONG,DIN			!錐體右邊界與內半徑
K,5,LONG-LL,DIN			!錐體空心部分左邊界控制點
K,6,,DOUT-LONG*TAN(THETA)	!錐體左端外半徑
K,7,LONG,DOUT			!錐體右端外半徑
K,8,LONG-(DD-UU),DOUT-DD*TAN(THETA)		!套筒右端
K,9,LONG-LL,DOUT-LL*TAN(THETA)			
K,10,LONG-(DD-UU)-LM,DOUT-(LM+DD)*TAN(THETA)	!套筒左端
K,11,LONG-(DD-UU),DV		!套筒外半徑
K,12,LONG-(DD-UU)-LM,DV		
K,50,LONG-XX,DIN		!試樣左邊界與半徑
K,51,LONG-XX			!試樣軸心
!建立平面模型
A,3,5,9,6			
A,5,4,7,9
A,10,8,11,12
WPROT,,,-90			!繞Y軸旋轉-90度
CYL4,,,,,DIN			!建立實心圓面
WPOFFS,,,XX-LONG		!移動工作平面到錐內試樣底部
CYL4,,,,,DIN			!建立實心圓面
WPCSYS,				!工作平面回到笛卡兒坐標系
ASEL,S,AREA,,4,5,1		!選擇面
ASBW,ALL,,DELETE		!用工作平面切分面
WPROT,,-90			!繞X軸旋轉-90度
ASEL,S,AREA,,6,9,1		!選擇面
ASBW,ALL,,DELETE		!用工作平面切分面
ALLS				!選擇所有的元素
!劃分網格
!劃分面網格(錐體)
TYPE,1				
MAT,2
LSEL,S,LINE,,12,15,1		!選擇控制線
LESIZE,ALL,,,90/NN		!設置劃分數90/NN
LSEL,S,LINE,,22,25,1		!選擇控制線
LESIZE,ALL,,,90/NN		!設置劃分數90/NN
ASEL,S,AREA,,4,5,1		!選擇面用于分網
ASEL,A,AREA,,10,11,1		!增選面用于分網
AMESH,ALL			!劃分網格
LSEL,S,LINE,,2,6,2
LESIZE,ALL,,,3			!設置劃分數3
LSEL,S,LINE,,1,3,2
LESIZE,ALL,,,2			!設置圓劃分數2
LSEL,S,LINE,,5,7,2
LESIZE,ALL,,,10			!設置劃分數10
ASEL,S,AREA,,1,2,1
AMESH,ALL
!實體模型(錐體)
TYPE,2				!設置單元類型			
MAT,2				!設置材料屬性
EXTOPT,ACLEAR,1			!旋轉后清除面
EXTOPT,ESIZE,90/NN-1		!旋轉后單元份數
ASEL,S,AREA,,2,2,1
VROTAT,ALL,,,,,,1,2,90-NN,1	!由面旋轉生成單元和實體
LOCAL,11,1,,,,,,90		!定義局部坐標系用于復制實體
CSYS,11
VGEN,4,ALL,,,,90		!復制時角度增量90
CSYS,0				!回到地球坐標系
EXTOPT,ACLEAR,1			!旋轉后清除面
EXTOPT,ESIZE,90/NN		!旋轉后單元份數
ASEL,S,AREA,,1,1,1
VROTAT,ALL,,,,,,1,2,,4		!由面旋轉生成單元和實體
EXTOPT,ACLEAR,1			!清除面
EXTOPT,ESIZE,2			!延伸后單元份數
ASEL,S,AREA,,4,5,1
ASEL,A,AREA,,10,11,1
VEXT,ALL,,,LONG-LL		!由面延伸生成單元和實體
NUMMRG,ALL			!合并元素
NUMCMP,ALL			!壓縮編號
ESEL,S,MAT,,2
CM,ZHUITI,ELEM			!定義錐體單元ZHUITI組件
VSEL,S,MAT,,2
CM,ZT,VOLU			!定義錐體ZHUITI組件
!劃分面網格(套筒)
TYPE,1				
MAT,3
LSEL,S,LINE,,8,10,2
LESIZE,ALL,,,5			!設置長度劃分數
LSEL,S,LINE,,9,11,2
LESIZE,ALL,,,3			!設置厚度劃分數
ASEL,S,AREA,,3,3,1
AMESH,ALL
!實體模型(套筒)
TYPE,2				
MAT,3
EXTOPT,ACLEAR,1			!清除面
EXTOPT,ESIZE,10			!旋轉方向每個實體單元份數
ASEL,S,AREA,,3,3,1
VROTAT,ALL,,,,,,1,2,,4
ESEL,S,MAT,,3
CM,TAOTONG,ELEM			!定義套筒單元TAOTONG組件
VSEL,S,MAT,,3
CM,TAOT,VOLU
!劃分面網格(試樣)
TYPE,1	
MAT,1				
LSEL,S,LINE,,16,19,1
LSEL,A,LINE,,26,29,1
LESIZE,ALL,,,6
ASEL,S,AREA,,12,15,1
AMESH,ALL
!實體模型(試樣)
TYPE,2				
MAT,1
EXTOPT,ACLEAR,1			!清除面
EXTOPT,ESIZE,5			!延伸后單元份數
ASEL,S,AREA,,12,15,1
VEXT,ALL,,,LS-(LONG-XX)
ESEL,S,MAT,,1
CM,SHIYANG,ELEM			!定義試樣單元SHIYANG組件
VSEL,S,MAT,,1
CM,SY,VOLU
!接觸單元設置
!注意事項:
!1、當凸面與平面或凹面接觸時，應指定平面或凹面為目標面
!2、應指定細網格所在面作為接觸面，粗網格所在面為目標面
!3、應指定較硬的面作為目標面，較軟的面為接觸面
!4、應指定高階的基礎單元所在面作為接觸面
!錐體外表面接觸單元
ALLS
CMSEL,S,ZT			!選擇組件ZT
CMPLOT				!顯示組件
ASEL,S,AREA,,8,8,1
ASEL,A,AREA,,19,31,6
CM,OUT_Z,AREA			!定義組元OUT_Z,錐體外表面
TYPE,3				!接觸單元
MAT,2
REAL,11
NSLA,S,1
NSEL,R,LOC,X,LONG-(DD-UU)-LM-10,LS
ESLN,S,0
ESURF,ALL			!在已存在單元自由表面覆蓋產生單元
!套筒內表面目標單元	
CMSEL,S,TAOT			!選擇組件TAOT
CMPLOT				!顯示組件
ASEL,S,AREA,,62,77,5
CM,TONG,AREA			!定義組元TONG,套筒內表面
TYPE,4				!目標單元
MAT,3
REAL,11
NSLA,S,1
NSEL,R,LOC,X,LONG-(DD-UU)-LM-10,LS
ESLN,S,0
ESURF,ALL                     
!錐體內表面目標單元
CMSEL,S,ZT			!選擇組件ZT
CMPLOT				!顯示組件
ASEL,S,AREA,,6,6,1
ASEL,A,AREA,,17,29,6
CM,IN_Z,AREA			!定義組元IN_Z,錐體內表面
TYPE,4				!目標單元
MAT,2
REAL,10
NSLA,S,1
NSEL,R,LOC,X,LONG-(DD-UU)-LM-10,LS
ESLN,S,0
ESURF,ALL  
!試樣外表面接觸單元			
CMSEL,S,SY			!選擇組件SY
CMPLOT				!顯示組件
ASEL,S,AREA,,82,86,4
ASEL,A,AREA,,89,92,3
CM,SHIY,AREA			!定義組元SHIY,試樣外表面
TYPE,3 				!接觸單元
MAT,1
REAL,10
NSLA,S,1
NSEL,R,LOC,X,LONG-(DD-UU)-LM-10,LS
ESLN,S,0
ESURF,ALL			!在已存在單元自由表面覆蓋產生單元
SAVE				!保存模型
FINI
!求解
/SOLU
ANTYPE,STATIC
NSEL,S,LOC,X,0
D,ALL,ALL			!左端夾頭末端固定約束
NSEL,S,LOC,X,LS
D,ALL,ALL			!右端試樣末端固定約束
CSYS,11
NSEL,S,LOC,Z,LONG-(DD-UU)-LM,LONG-(DD-UU)
NSEL,R,LOC,X,DV
D,ALL,ALL,0			!套筒外表面約束
CSYS
!過盈配合分析
TIME,1				
AUTOTS,ON			!打開自動時間步
NSUB,STEP1			!指定此荷載步的子步數
NROPT,FULL
OUTRES,ALL,ALL			!規定寫入數據庫的求解信息
NLGEOM,ON			!打開幾何非線性
PRED,ON                 
LNSRCH,ON   
ALLS
SOLVE
!摩擦分析
TIME,2				
AUTOTS,ON			!打開自動時間步
NSUB,STEP2			!指定此荷載步的子步數
NROPT,FULL
OUTRES,ALL,ALL			!規定寫入數據庫的求解信息
NLGEOM,ON			!打開幾何非線性
PRED,ON                 
LNSRCH,ON
ALLS
NSEL,S,LOC,X,LS
D,ALL,UX,LS+XX+5		!將柱體抽出
ALLS
SOLVE
SAVE
FINI
!夾頭構件的應力應變云圖
/POST1
SET,1,LAST			!過盈配合最后的子步
ALLS
PLNSOL, EPTO,EQV, 0,1.0    	!顯示等效應變
PLNSOL, S,EQV, 0,1.0    	!顯示等效應力(節點解)
PLESOL, S,EQV, 0,1.0    	!顯示等效應力(單元解)
CMSEL,S,ZHUITI  		!錐體單元
CMPLOT
PLNSOL, S,EQV, 0,1.0		!顯示等效應力
CMSEL,S,TAOTONG			!套筒單元
CMPLOT
PLNSOL, S,EQV, 0,1.0		!顯示等效應力
CMSEL,S,SHIYANG			!試樣單元
CMPLOT
PLNSOL, S,EQV, 0,1.0		!顯示等效應力
ESEL,S,ENAME,,173		!選擇接觸單元
PLNSOL, S,EQV, 0,1.0		!顯示等效應力
PLNSOL, CONT,PRES, 0,1.0	!顯示接觸壓力
FINI
!獲取夾頭提供的摩擦合力
!方法一
/POST1
AXIAL=0
TEMP=0
TEMP1=0
SET,1,LAST
*GET,TEMP,ACTIVE,0,SOLU,NCMSS		!獲取累積的載荷子步
SET,2,LAST
*GET,TEMP1,ACTIVE,0,SOLU,NCMSS		!獲取累積的載荷子步
TEMP=TEMP+TEMP1
*DIM,FORCE1,ARRAY,TEMP
*DO,I,1,TEMP
	*IF,I,EQ,1,THEN
	SET,1,1
	NSEL,S,LOC,X,LS			!選擇節點
	*GET,NMAX,NODE,,NUM,MAX
	*GET,NMIN,NODE,,NUM,MIN
	*DO,KK,NMIN,NMAX
		*IF,NSEL(KK),EQ,1,THEN
		*GET,AXIAL,NODE,KK,RF,FX
		FORCE1(I)=FORCE1(I)+AXIAL
		*ENDIF
	*ENDDO
	*ELSEIF,I,NE,1
	SET,NEXT
	NSEL,S,LOC,X,LS			!選擇節點
	*GET,NMAX,NODE,,NUM,MAX
	*GET,NMIN,NODE,,NUM,MIN
	*DO,KK,NMIN,NMAX
		*IF,NSEL(KK),EQ,1,THEN
		*GET,AXIAL,NODE,KK,RF,FX
		FORCE1(I)=FORCE1(I)+AXIAL
		*ENDIF
	*ENDDO
	*ENDIF
*ENDDO
FINI
/POST26
ALLS
NUMVAR,20			!指定允許的變量數
ANSOL,2,4012,EPEL,X,ST_N4012	!獲取節點X方向應力
ANSOL,3,4015,EPEL,X,ST_N4015	!獲取節點X方向應力
ANSOL,4,4006,EPEL,X,ST_N4006	!獲取節點X方向應力
ANSOL,5,4009,EPEL,X,ST_N4009	!獲取節點X方向應力
ANSOL,6,4023,EPEL,X,ST_N4023	!獲取節點X方向應力
NSOL,7,4012,U,X,UX_N4012  	!取出節點的X向位移數據
/AXLAB,X,TIME[S]
/AXLAB,Y,STRESS[MPA] DISPLACE[MM]
XVAR,1
PLVAR,2,3,4,5,6,7		!顯示2～7變量隨時間變化的圖形
*DIM,TIME1,ARRAY,TEMP
VGET,TIME1(1),1
*DIM,LABEL,CHAR,1,2		!載荷-時間曲線保存
LABEL(1,1)='TIME1'
LABEL(1,2)='FORCE1'
*CREATE,ANSUITMP
*CFOPEN,TIME1_FORCE1,DAT
*VWRITE,LABEL(1,1),LABEL(1,2)
(1X,A8,A8)
*VWRITE,TIME1(1),FORCE1(1)
(1X,F18.6,F18.6)
*CFCLOSE
*END
/INPUT,ANSUITMP
!方法二
/POST1
TEMP=0
TEMP1=0
SET,1,LAST
*GET,TEMP,ACTIVE,0,SOLU,NCMSS		!獲取累積的載荷子步
SET,2,LAST
*GET,TEMP1,ACTIVE,0,SOLU,NCMSS		!獲取累積的載荷子步
TEMP=TEMP+TEMP1
*DIM,FORCE2,ARRAY,TEMP
*DO,I,1,TEMP
	*IF,I,EQ,1,THEN
	SET,1,1
	NSEL,S,LOC,X,LS			!選擇節點
	FSUM
	*GET,FORCE2(I),FSUM,0,ITEM,FX
	*ELSEIF,I,NE,1
	SET,NEXT
	NSEL,S,LOC,X,LS			!選擇節點
	FSUM
	*GET,FORCE2(I),FSUM,0,ITEM,FX
	*ENDIF
*ENDDO

/POST26
NUMVAR,20			!指定允許的變量數
ANSOL,8,3889,EPEL,EQV,ST_N3889	!獲取節點等效應力
ANSOL,9,4066,EPEL,EQV,ST_N4066	!獲取節點等效應力
ANSOL,10,4067,EPEL,EQV,ST_N4067	!獲取節點等效應力
ANSOL,11,3887,EPEL,EQV,ST_N3887	!獲取節點等效應力
ANSOL,12,4058,EPEL,EQV,ST_N4058	!獲取節點等效應力
ANSOL,13,4059,EPEL,EQV,ST_N4059	!取出節點的X向位移數據
/AXLAB,X,TIME[S]
/AXLAB,Y,STRESS[MPA]
XVAR,1
PLVAR,8,9,10,11,12,13		!顯示8～13變量隨時間變化的圖形
*DIM,TIME2,ARRAY,TEMP
VGET,TIME2(1),1
*DIM,LABEL,CHAR,1,2		!載荷-時間曲線保存
LABEL(1,1)='TIME2'
LABEL(1,2)='FORCE2'
*CREATE,ANSUITMP
*CFOPEN,TIME2_FORCE2,DAT
*VWRITE,LABEL(1,1),LABEL(1,2)
(1X,A8,A8)
*VWRITE,TIME2(1),FORCE2(1)
(1X,F18.6,F18.6)
*CFCLOSE
*END
/INPUT,ANSUITMP