基于GMC/umat的复合材料宏细观渐近损伤分析(一)

 

近期在开展基于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.

 

转载于:https://www.cnblogs.com/xym16805/p/11261842.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值