基于VB的ANSYS二次開發之Duncan-Chang本構模型算法介紹
所謂的本構模型是指材料的應力應變關系，比如遵循虎克定律的線彈性應力應變關系就是一種常用的本構模型。ANSYS為用戶提供了許多本構模型，但在某些特殊領域，現有的本構模型卻很少，完全不能滿足分析要求。為了解決這個問題，ANSYS為用戶提供了usermat等UPFs用戶子程序，這些用戶子程序擁有強大的二次開發功能，可以實現各種復雜本構模型的開發。但是，對于一些簡單的本構模型，用戶也可以利用APDL語言進行開發，比如Duncan-Chang本構模型。
Duncan-Chang本構模型介紹
對于應力應變關系復雜的材料而言，變形主要是由顆粒間位置的變化引起的，不同應力水平下相同的應力增量引起的應變增量并不相同，從而表現出復雜的非線性特征。為了反映材料特性，人們提出了許多本構模型，以期更好的反映材料的應力應變規律。
Duncan-Chang本構模型主要優點是可以利用常規三軸固結排水剪試驗確定模型參數，因此在工程實踐中得到了廣泛應用。Duncan-Chang本構模型是非線性彈性模型，彈性矩陣中的彈性系數不再是常量，而是隨應力狀態而改變。由于不考慮塑性變形，所以一般只適用于載荷不大的情況（即不太接近破壞的情況）。Duncan-Chang模型有E-V和E-B模型兩類。
當然，該模型也存在一些缺陷，如無法反應不同應力路徑的影響、加卸載判斷不明確等，不可避免造成了計算分析誤差，長期以來許多學者試圖來對其進行修正。有限元軟件通常采用傳統塑性理論的應力符號，以拉應力為正，下面對Duncan-Chang模型采用有限元軟件的應力形式進行說明。
Duncan-Chang本構模型算法
該模型是Duncan和Chang根據Konder關于巖土材料的三軸試驗的偏應力與軸向應變近似呈雙曲線的假定而提出的。根據Konder的建對應變硬化的模型，其試驗的應力應變可用雙曲線關系來近似地描述，如圖3-1所示，即在不變時：
式中a和b是試驗常數。
上式也可以寫成
其斜率為b,截距為a。
Duncan和Chang利用上述關系推導出了彈性模量公式：
可見增量虎克定律中所用的彈性模量實際是常規試驗曲線的切線斜率。這樣的模量叫切線彈性模量，可用來表示，
詳細介紹見附件文檔?。。。?/strong>
生成并調用宏文件
在編程實現本構模型的過程中，需要重復執行某一部分，用戶可以將該部分獨立編寫后放入指定位置，由APDL命令來調用，也可以通過*CREATE命令創建宏文件，并用*END命令結束宏的創建。利用*USE命令調用宏文件，并向宏文件傳遞參數：
*USE,Name,ARG1,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,AR10,AR11,AR12,AR13,AR14,AR15,AR16,AR17,AR18
其中，Name是宏文件名，ARGI到AR18是宏文件用到的參數值。我們將用*CREATE命令創建名為Duncan-Chang的宏文件，其中包含9個參數，使用*USE命令對模型內的每個單元反復調用Duncan-Chang宏文件，不斷計算得到新的切線模量。
APDL實現過程
Duncan-Chang E-v模型是一種建立在增量廣義虎克定律基礎上的非線性變彈性模型，是通過不斷改變其切線彈性模量來實現非線性的，完全可以通過ANSYS APDL進行編程分析。計算過程中主要通過如下方式來實現：取初始材料參數，施加第一步載荷，計算并讀取單元應力，根據單元的當前應力調用Duncan-Chang模型宏命令計算新的材料參數（主要是材料的彈性模量和泊松比），代替初始材料參數；施加第二步載荷，計算并讀取應力增量，根據單元的當前應力調用Duncan-Chang模型宏命令計算新的材料參數，以此類推。下面給出了 Duncan-Chang E-v模型創建宏文件，并用來模擬壓縮試驗的APDL詳細代碼(\chp6\Duncan-E-v.inp):
?。?）定義單元類型：
/PREP7
ET,1,185
?。?）設置材料屬性
MPTEMP,,,,,,,,
MPTEMP,1,0
MPDATA,EX,1,,2e11
MPDATA,PRXY,1,,0.3
MPTEMP,,,,,,,,
MPTEMP,1,0
MPDATA,DENS,1,,7800
!3.建立幾何模型
！（1）顯示工作平面
WPSYYLE,,,,,,,,1
!(2)建立模型
BLOCK,0,10,0,20,0,2,
/REPLO
wpoff,3,5,0
/REPLO
CYLIND,1, ,0,10,0,360,
/REPLO
wpoff,0,10,0
CYLIND,1, ,0,10,0,360,
wpoff,0,-5,0
CYLIND,1, ,0,10,0,360,
wpoff,4,0,0
CYLIND,1, ,0,10,0,360,
wpoff,0,5,0
CYLIND,1, ,0,10,0,360,
wpoff,0,-10,0
CYLIND,1, ,0,10,0,360,
?。?）體布爾操作
FLST,3,6,6,ORDE,2
FITEM,3,2
FITEM,3,-7
VSBV, 1,P51X
！4.生成有限元模型
/REPLO
FINISH
/AUX12
FINISH
/PREP7
MSHAPE,1,3D
MSHKEY,0
CM,_Y,VOLU
VSEL, , , , 8
CM,_Y1,VOLU
CHKMSH,'VOLU'
CMSEL,S,_Y
VMESH,_Y1
CMDELE,_Y
CMDELE,_Y1
CMDELE,_Y2
！5.加載以及求解
?。?）加載
//REPLO
VPLOT
FINISH
/SOL
FLST,2,1,5,ORDE,1
FITEM,2,3
/GO
DA,P51X,ALL,0
/REPLO
FLST,2,1,4,ORDE,1
FITEM,2,7
/GO
FLST,2,1,5,ORDE,1
FITEM,2,4
/GO
!*
SFA,P51X,1,PRES,1000
!*
/STATUS,SOLU
SOLVE
(2) 創建Duncan-Chang計算
！
！***********
*dim,Str_max,array,120
*dim,S_max,ARRAY,120
*create,Duncan-Chang.mac !創建Duncan-Chang計算的宏命令
*afun，deg ！角度單位用度，不是弧度
*set，Pa，le5
*set，sigm1，-ArrS3(i) ！ArrS3為第三主應力，拉負壓正
*set，sigm3，-ArrS3(i) ！ArrS3為第一主應力，拉負壓正
*if，sigm3，LT,0.1*Pa，then
sigm3=0.1*Pa
*endif
str=2*(c*cos（Fai）+sigm3*sin（Fai））/(1-sin(Fai)) !破壞時的主應力差
S=(sigml-sigm3)/Slr !應力水平
*if,S.GT,0.95,then
S=0.95 !應力水平最大取值
*endif
!判斷加卸綾荷，如果(sigml-sigm3)、S均小于歷史最大值時視為卸載一再加載過程
*if,Str_max(i),gl,sigml-sigm3,and,S_max(i),gt,S,thcn Et=Kur*Pa*(sigm3/Pa)**n !卸載模置
*elseif，Str_max(i),gt,sigm1-sigm3，and，S_max(i),gt,S,then Ei=k*Pa*(sigm3/Pa)**n
Et=Ei*(l-Rf*S)**2 !加載情況的切線模量
S_max(i)=S ！保存歷史最大應力水平
*elseif,Str_max(i),le,sigml-sigm3,and,S_max(i),gt,S,then Ei=k*Pa*(sigm3/Pa)**n
Et=Ei*(1-Rf*S)**2 !加載情況的切線模量
Str_max(i)=sigm1-sigm3 ！保存歷史最大應力
*elseif，Str_max(i),gt,sigm1-sigm3，and，S_max(i),gt,S,then Ei=k*Pa*(sigm3/Pa)**n
Et=Ei*(l-Rf*S)**2 !加載情況的切線模量
Str_max(i)=sigm1-sigm3 !保存歷史最大主應力差
S_max(i)=S ！保存歷史最大應力水平
*endif
Ei=k*Pa*(sigm3/Pa)**n
A=(sigm 1 -sigm3)*D/(Ei*( 1 -Rf*S))
Mu=(G-F*loglO(sigm3/Pa)/(l-A)**2 ！計算泊松比
*if,Mu,GE,0.49,then
Mu=0.49
*endif
Mp,ex,I,Et !修改單元i的Et
Mp,nuxy,I,Mu
Mpchg,i,i !修改i單元的材料屬性，mpchg,mat,elem
!
*vwrite,Et !格式化輸出Et，用于觀察中間計算結果
（E13.5）
*end !結束創建宏文件
！
！*****************************************
*dim ArrSl,array,120
*dim ArrS3,array,120
*do,j,2,9
!取出計算結果，修改彈性模*
/post1
etable,etabs 1,s,1 !取各單元第一主應力
etable,etabs 3,s,3 !取各單元第一主應力
*do,num,1,120 !num為單元編號
*get，Arrs1(num),elem,num,etab,etabs1 !將單元結果存入數組
*get，Arrs3(num),elem,num,etab,etabs3
*enddo
/prep7
*do,i,1,120,1 !各單元依次計算
c=130E3
Fai=31.2
Rf=0.9
k=538
n=0.35
Kur=645
G-0.33
F=0.033
D=3.2
!c一粘聚力，Fai—內摩擦角，Rf—破壞比，k—初始切線模最數，
!n—指數，Mu-泊松比， Kur—考慮卸載再加載時的模量數,G、F、D一參數
*use,Duncan-Chang.mac,c,Fai,Rf,k,n,Kur,G,F,D!調用Duncan-Chang宏文件
*enddo
！*******************************************
/solu
Time,j !定義第j載荷步
Nsubst,20,100,5 !指定子步數
Sfl,3,pres,3e5+le5*j !施加載荷，增量100kPa，最終上表面壓力為1.2MPa
Outres,all,all
Solve
*enddo
參考文獻
[1] 師訪編.ANSYS二次開發及應用實例講解[M].中國水利水電出版社.2012.1
[2] 隋麗娜，遲劍，郭立峰編. Visual Basic范例開發大全[M].清華大學出版社.
[3] 廖孟柯編. 基于VB的ANSYS二次開發與應用[J].
[4] 張艷.Visual Basic程序設計教程[M].徐州:中國礦業大學出版社,2001.
[5] 張晉西.用VB增強ANSYS前處理能力[J].計算機應用,2002,22(3).
[6] 龔曙光,謝桂蘭. ANSYS操作命令與參數化編程[M].北京:機械工業出版社,2004.
工程師必備
- 項目客服
- 培訓客服
- 平臺客服
TOP
