各向同性彈塑性本構的vumat源代碼：通過修改umat
瀏覽：3175 收藏：1
1 vumat與umat的區別
從程序實現的角度，我們重點關注以下幾點區別：
? vumat不需要輸出一致性切線剛度矩陣
? vumat中應力應變存儲順序與umat不同
? vumat中存儲的應變值為張量應變值，而umat中為工程應變
? vumat的應力和狀態變量的更新方式不同，其分為old和new兩個數組
Abaqus/Explicit在啟動計算前，會進行數據檢查，在檢查的過程中會給定一組虛假的應變來檢查程序的流程。在檢查時，傳入vumat的totalTime和stepTime都為0，根據用戶給定的本構關系，程序進行計算并得到初始的穩定時間增量。如果這個穩定時間增量太大，就會導致計算不穩定（不收斂），所以需要給出彈性的計算過程，以保證得到一個比較合適的初始穩定時間增量。
vumat與umat的更多對比見下圖
2 代碼修改
從umat的源代碼出發，作出相應修改可以得到vumat。
首先增加應力和應變分量的轉換函數
!******************************************************************************* ! transfer_strain_vumat2umat：將vumat接口中的應變變量轉化為umat接口的應變變量 ! vumat的應變變量為張量應變，且存儲順序不同 ! 變量說明 ! 輸入: ! strain_vumat : vumat接口中應變張量的存儲方式(e11,e22,e33,e12,e23,e31) ! ! 輸出： ! strain_umat ：umat接口中應變張量的存儲方式(e11,e22,e33,2*e12,2*e13,2*e23) !******************************************************************************* function transfer_strain_vumat2umat(strain_vumat) result(strain_umat) real(8),intent(in) :: strain_vumat(6) real(8) :: strain_umat(6) integer :: i integer,parameter :: order(6)=[1,2,3,4,6,5] do i=1,6 strain_umat(i) = strain_vumat(order(i)) enddo strain_umat(4:6) = 2.0d0 * strain_umat(4:6) return end function transfer_strain_vumat2umat !******************************************************************************* ! transfer_stress_vumat2umat：將vumat接口中的應力變量轉化為umat接口的應力變量 ! vumat的存儲順序不同 ! 變量說明 ! 輸入: ! stress_vumat : vumat接口中應力張量的存儲方式(s11,s22,s33,s12,s23,s31) ! ! 輸出： ! stress_umat ：umat接口中應力張量的存儲方式(s11,s22,s33,s12,s13,s23) !******************************************************************************* function transfer_stress_vumat2umat(stress_vumat) result(stress_umat) real(8),intent(in) :: stress_vumat(6) real(8) :: stress_umat(6) integer :: i integer,parameter :: order(6)=[1,2,3,4,6,5] do i=1,6 stress_umat(i) = stress_vumat(order(i)) enddo return end function transfer_stress_vumat2umat
然后增加vumat函數
!******************************************************************************* ! vumat：更新應力以及狀態變量 ! 變量說明 ! 輸入: ! strain_n0 : n時刻的應變 ! strain_inc : 應變增量 ! stress_n0 ：應力張量 ! statev_n0 : 狀態變量 ! ! 輸出： ! stress_n1 ：n+1時刻的應力 ! statev-n1 : 狀態變量 ! 注意： ! 這里調用的函數都是遵循之前umat的接口，所以一切按照umat的定義標準來 ! 最后在主程序中調用時，先將vumat的變量轉化為umat相關的變量，計算完畢之后再轉化為vumat !******************************************************************************* subroutine plastic_iso_vumat(mu,kappa, n_pt,hard_data, strain_n0,strain_inc, stress_n0,stress_n1,statev_n0,statev_n1) real(8),intent(in) :: mu,kappa integer,intent(in) :: n_pt real(8),intent(in) :: hard_data(n_pt,2) real(8),intent(in) :: strain_n0(6) real(8),intent(in) :: strain_inc(6) real(8),intent(in) :: stress_n0(6) real(8),intent(out) :: stress_n1(6) real(8),intent(in) :: statev_n0(8) real(8),intent(out) :: statev_n1(8) ! 與應變相關的變量 real(8) :: strain_n1(6) real(8) :: strain_el_n1_trial(6) real(8) :: strain_el_n1_trial_vol real(8) :: strain_el_n1_trial_dev(6) ! 與應力相關的變量 real(8) :: s_n1_trial(6) real(8) :: p_n1_trial real(8) :: q_n1_trial real(8) :: n_trial(6) real(8) :: s_n1(6) ! 狀態變量 real(8) :: alpha_n0,alpha_n1 real(8) :: strain_pl_n0(6),strain_pl_n1(6) real(8) :: Dp_0,Dp_1 ! 屈服函數 real(8) :: sig_y0 real(8) :: f_trial real(8) :: gamma_inc ! 彈性矩陣 real(8) :: E,nu real(8) :: De(6,6) ! 計算全應變 strain_n1 = strain_n0 + strain_inc ! 計算彈性矩陣 E = 9.0 * kappa * mu / (3.0 * kappa + mu) nu = ( 3.0 * kappa - 2.0 * mu ) / (3.0 * kappa + mu) / 2.0 De = plastic_iso_compute_elasD(E,nu) ! 讀取n時刻的狀態變量 alpha_n0 = statev_n0(1) strain_pl_n0 = statev_n0(2:7) Dp_0 = statev_n0(8) ! 計算試驗狀態的值 call plastic_iso_compute_eps_trial(strain_n0, strain_inc, strain_pl_n0, & strain_el_n1_trial,strain_el_n1_trial_vol,strain_el_n1_trial_dev) call plastic_iso_compute_sig_trial(mu,kappa,strain_el_n1_trial_vol,strain_el_n1_trial_dev,& p_n1_trial,s_n1_trial,q_n1_trial,n_trial) sig_y0 = plastic_iso_compute_sigma_y(n_pt,hard_data,alpha_n0) f_trial = plastic_iso_compute_f_trial(q_n1_trial,sig_y0) ! 開始return-mapping算法 if (f_trial >0) then ! 發生了塑性屈服，需要計算塑性乘子增量 gamma_inc = plastic_iso_compute_gamma_inc(mu, n_pt,hard_data,alpha_n0,q_n1_trial) ! 更新應力 s_n1 = plastic_iso_compute_sigDev_n1(mu,s_n1_trial,gamma_inc) stress_n1 = plastic_iso_compute_stress_n1(mu,p_n1_trial,s_n1_trial,gamma_inc) ! 更新狀態變量 alpha_n1 = alpha_n0 + gamma_inc strain_pl_n1 = plastic_iso_compute_strain_pl(mu,strain_n1,s_n1,strain_el_n1_trial_vol) Dp_1 = plastic_iso_compute_dissipation(Dp_0, stress_n0,stress_n1, strain_pl_n0,strain_pl_n1) statev_n1(1) = alpha_n1 statev_n1(2:7) = strain_pl_n1 statev_n1(8) = Dp_1 else ! 沒有發生塑性屈服 ! 更新應力 stress_n1 = matmul(De,strain_el_n1_trial) ! 更新狀態變量 alpha_n1 = alpha_n0 strain_pl_n1 = strain_pl_n0 Dp_1 = Dp_0 statev_n1(1) = alpha_n1 statev_n1(2:7) = strain_pl_n1 statev_n1(8) = Dp_1 endif end subroutine plastic_iso_vumat
3 算例
3.1 單單元拉伸測試
對單個單元進行單軸拉伸，邊界條件如下：
von Mises應力對比結果如下（左圖為Abaqus材料庫計算，右圖為vumat子程序計算結果）：
等效塑性應變對比結果如下（左圖為Abaqus材料庫計算，右圖為vumat子程序計算結果）：
反力曲線對比如下：
塑性耗散曲線對比如下：
3.2 圓棒拉伸測試
對一圓棒骨料進行單軸拉伸，其邊界條件如下：
von Mises應力對比結果如下（左圖為Abaqus材料庫計算，右圖為vumat子程序計算結果）：
等效塑性應變對比結果如下（左圖為Abaqus材料庫計算，右圖為vumat子程序計算結果）：
反力曲線對比如下：
算例cae模型
以下內容為付費內容，請購買后觀看
包含1個文件
付費內容包括vumat的源代碼和文中所有算例文件
release.zip
11.50MB
技術鄰APP
工程師必備
工程師必備
- 項目客服
- 培訓客服
- 平臺客服
TOP
1
