近期在开展基于GMC/umat的复合材料宏细观渐近损伤分析，一些技术细节分享如下：
1.理论基础
针对连续纤维增强复合材料，可以通过离散化获得如下的模型：
（a）（b）（c）
图1 连续纤维增强复合材料细观离散化（a）代表性体积单元示意图（b）2*2子胞划分（c）Nr*Nb子胞划分
在子胞中采用一阶线性位移模式：
子胞的应变由几何方程给出：
单胞的平均应变可由体积平均给出：
本构关系如下所示：
单胞的平均应力可由体积平均给出：
需要注意的是，由于采用的是线性位移模式，则子胞内为常应变和常应力。子胞内应变为子胞平均应变，而应力为平均应力。
由子胞边界的位移连续可以得到单胞应变与子胞应变的关系：
写成矩阵形式如下：
由子胞间的应力连续可得到单胞应力与子胞应力的关系，如下所示：
为了进一步包含子胞的损伤信息，进而开展损伤后力学分析，为每个子胞分配一个状态变量 ，其表征了损伤的状态，0 为未损伤，0.9999 为损伤。从而可推导出如下考虑损伤的子胞间应力连续方程：
写成矩阵形式有：
将子胞位移连续与应力连续方程合并，可得如下完备的通用单胞计算模型：
求解可得：
将子胞应变矩阵分解为子矩阵 ：
由此可获得任意子胞的应变表达式：
可获得复合材料宏观本构关系：
2.umat实现（Fortran代码）
以下为 usermat 接口及主程序，主要包含参数传递、当前刚度（考虑损伤）计算、切线刚度矩阵赋值、更新应力，考虑篇幅，星号为部分省略区，只给出关键部分。
1 *deck,usermat USERDISTRIB parallel gal 2 subroutine usermat( 3 & matId, elemId,kDomIntPt, kLayer, kSectPt, 4 & ldstep,isubst,keycut, 5 **************************************************** 6 **************************************************** 7 c------------------------------------1.自定义变量声明------------------------------------- 8 real*8 Q(6,6),mic_sig(2,2,6),mic_dam(2,2),mic_dam2(2),sstrain(6,1) 9 INTEGER i,j,ii,jj,KK1,Kk2,Nb,Nr 10 real*8 ex_f,ey_f,ez_f,gxy_f,gxz_f,gyz_f,prxy_f,prxz_f,pryz_f 11 & ,s_f(6,6),c_f(6,6), 12 & ex_m,prxy_m,g_m,s_m(6,6),c_m(6,6) 13 real*8 vf,k1,a,b,L,H 14 real*8 cell_H(2),cell_L(2) 15 integer cforcm(2,2) 16 real*8 mac_test_s(6,6),mac_e(6,1),s_temp(2,2,6) 17 character(5)::num 18 real*8 xt,xc,s,xt_f,xc_f 19 real*8 dam_ori(2,2) 20 c-------------------------------------动态数组 21 real (kind=8), allocatable:: ag(:,:), agg(:,:) 22 allocate(ag(294,294),agg(294,294)) 23 **************************************************** 24 **************************************************** 25 c------------------------------------2.目前刚度计算----------------------------------------- 26 do i=1,6 27 sstrain(i,1)=strain(i) 28 end do 29 call mic_sig_cal(sstrain,c_f,c_m,cell_H,cell_L, 30 & cforcm,mic_dam,mic_sig) 31 call mic_dam_cal(mic_sig,cforcm,xt,xc,s,xt_f, 32 & xc_f,cell_L,cell_H,mic_dam,dam_ori) 33 mac_test_s=0. 34 call mac_con_cal(c_f,c_m,cell_H,cell_L, 35 & cforcm,mic_dam,mac_test_s) 36 Q=mac_test_sc------------------------------------3.一致切线算子矩阵---------------------------------- 37 dsdePl=Q 38 c------------------------------------4.更新应力---------------------------------- 39 DO Kk1=1, ncomp 40 DO Kk2=1, ncomp 41 stress(Kk2)=stress(Kk2)+dsdePl(Kk2,Kk1)*dStrain(Kk1) 42 END DO 43 END DO 44 **************************************************** 45 **************************************************** 46 RETURN 47 END 48 c------------------------------------------------------------------------------------------------
以下为主要子程序之细观应力计算：
1 subroutine mic_sig_cal(sstrain,c_f,c_m,cell_H,cell_L, 2 & cforcm,mic_dam,mic_sig) 3 #include "impcom.inc" 4 real*8 sstrain(6,1),mic_sig(2,2,6),mic_dam(2,2) 5 integer Nb,Nr 6 real*8 c_f(6,6),c_m(6,6) 7 real*8 cell_H(2),cell_L(2) 8 integer cforcm(2,2) 9 real*8 ag(24,24),bg(24,6) 10 real*8 bgg(24,1) 11 real*8 temp(1,24) 12 real*8 mic_eps(24,1) 13 integer i,j,k,ii,jj,kk,iii,jjj,kkk,countt 14 **************************************************** 15 **************************************************** 16 c eps12 17 do j=1,Nb 18 countt=countt+1 19 do i=1,Nr 20 ii=(i-1)*Nb+j 21 ag(countt,(ii-1)*6+4)=cell_L(i) 22 end do 23 bg(countt,4)=1. 24 end do**************************************************** 25 **************************************************** 26 c sig21 27 do j=1,Nb 28 do i=1,Nr-1 29 countt=countt+1 30 temp=0. 31 ii=(i-1)*Nb+j 32 iii=(i+1-1)*Nb+j 33 if (cforcm(i,j).eq.1)then 34 temp(1,(ii-1)*6+1:(ii-1)*6+6)=c_f(4,:)*(1-mic_dam(i,j)) 35 else 36 temp(1,(ii-1)*6+1:(ii-1)*6+6)=c_m(4,:)*(1-mic_dam(i,j)) 37 end if 38 if(cforcm(i+1,1).eq.1)then 39 temp(1,(iii-1)*6+1:(iii-1)*6+6)=-1*c_f(4,:)*(1-mic_dam(i,j)) 40 else 41 temp(1,(iii-1)*6+1:(iii-1)*6+6)=-1*c_m(4,:)*(1-mic_dam(i,j)) 42 end if 43 ag(countt,:)=temp(1,:) 44 bg(countt,:)=0. 45 end do 46 end do 47 **************************************************** 48 **************************************************** 49 RETURN 50 END 51 c----------------------------------------------------------------------------------
以下为主要子程序之细观损伤计算：
1 subroutine mic_dam_cal(mic_sig,cforcm,xt,xc,s,xt_f, 2 & xc_f,cell_L,cell_H,mic_dam,dam_ori) 3 #include "impcom.inc" 4 integer i,j,k 5 real*8 mic_sig(2,2,6) ,mic_dam(2,2),dam_ori(2,2) 6 integer cforcm(2,2) 7 real*8 cri 8 real*8 cri_flag,xt,xc,s,f1,f11,f44,f12 9 real*8 xt_f,xc_freal*8 matrix_f 10 real*8 cell_H(2),cell_L(2) 11 **************************************************** 12 **************************************************** 13 if( (cforcm(i,j).eq.0).and. 14 & ( mic_sig(i,j,1).ge.(1.0*xt).or. 15 & mic_sig(i,j,1).le.(-1.0*xc) ) ) then 16 mic_dam=0.9999 17 dam_ori(i,j)=0.9999 18 end if 19 **************************************************** 20 **************************************************** 21 if( (cforcm(i,j).eq.1).and. 22 & ( (mic_sig(i,j,1).ge.xt_f*1.0).or. 23 & (mic_sig(i,j,1).le.xc_f*1.0)) ) then 24 mic_dam=0.9999 25 dam_ori(i,j)=0.9999 26 end if 27 end do 28 end do 29 RETURN 30 END 31 c---------------------------------------------------------------------
以下为主要子程序之宏观刚度计算：
1 subroutine mac_con_cal(c_f,c_m,cell_H,cell_L, 2 & cforcm,mic_dam,mac_test_s) 3 #include "impcom.inc" 4 real*8 c_f(6,6),c_m(6,6),cell_H(2),cell_L(2), 5 & mic_dam(2,2),mac_test_s(6,6) 6 integer cforcm(2,2) 7 real*8 ag(24,24),ag_inv(24,24),a(24,6),bg(24,6) 8 real*8 temp(1,24),temp_a(6,6),temp_c(6,6) 9 real*8 temp_ca(6,6) 10 integer i,j,k,ii,jj,kk,iii,jjj,kkk,countt 11 integer Nb,Nr 12 **************************************************** 13 **************************************************** 14 c sig21do j=1,Nb 15 do i=1,Nr-1 16 countt=countt+1 17 temp=0. 18 ii=(i-1)*Nb+j 19 iii=(i+1-1)*Nb+j 20 if (cforcm(i,j).eq.1)then 21 temp(1,(ii-1)*6+1:(ii-1)*6+6)=c_f(4,:)*(1-mic_dam(i,j)) 22 else 23 temp(1,(ii-1)*6+1:(ii-1)*6+6)=c_m(4,:)*(1-mic_dam(i,j)) 24 end if 25 if(cforcm(i+1,1).eq.1)then 26 temp(1,(iii-1)*6+1:(iii-1)*6+6)=-1*c_f(4,:)*(1-mic_dam(i,j)) 27 else 28 temp(1,(iii-1)*6+1:(iii-1)*6+6)=-1*c_m(4,:)*(1-mic_dam(i,j)) 29 end if 30 ag(countt,:)=temp(1,:) 31 bg(countt,:)=0. 32 end do 33 end do 34 **************************************************** 35 **************************************************** 36 temp_ca=matmul(temp_c,temp_a) 37 mac_test_s=mac_test_s+cell_L(i)*cell_H(j)*temp_ca 38 end do 39 end do 40 return 41 end
3.参考文献
[1] Aboudi J. A continuum theory for fiber-reinforced elastic-viscoplastic composites[J]. International Journal of Engineering Science. 1982, 20(5): 605-621.
[2] Aboudi J. Elastoplasticity theory for porous materials ☆[J]. Mechanics of Materials. 1984, 3(1): 81-94.
[3] Aboudi J. Closed form constitutive equations for metal matrix composites[J]. International Journal of Engineering Science. 1987, 25(9): 1229-1240.
[4] Aboudi J. Micro-Failure Prediction of the Strength of Composite Materials under Combined Loading[J]. Journal of Reinforced Plastics & Composites. 1991, 10(5): 495-503.
[5] 雷友锋. 纤维增强金属基复合材料宏-细观统一本构模型及应用研究[D]. 南京航空航天大学, 2002.
[6] 孙志刚. 复合材料高精度宏-细观统一本构模型及其应用研究[D]. 南京航空航天大学, 2005.
[7] 高希光. 陶瓷基复合材料损伤耦合的宏细观统一本构模型研究[D]. 南京航空航天大学, 2007.
[8] Tang Z, Zhang B. Prediction of biaxial failure envelopes for composite laminates based on Generalized Method of Cells[J]. Composites Part B Engineering. 2012, 43(3): 914-925.
[9] 唐占文. 考虑界面相的复合材料宏-细观渐进损伤解析模型研究[D]. 哈尔滨工业大学, 2013.
[10]Pineda E J, Bednarcyk B A, Waas A M, et al. Progressive failure of a unidirectional fiber-reinforced composite using the method of cells: Discretization objective computational results[J]. International Journal of Solids & Structures. 2013, 50(9): 1203-1216.
扫一扫
1488
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
