【分享】简单的混凝土的本构程序

molen

module MConcrete !混凝土模板
  implicit none
  
  type::typ_Concrete
   !混凝土抗拉强度,抗压强度,初始弹性模量,初始泊松比
   !最终泊松比,初始剪切模量
   real*8 Ft,Fc,E0,MU0,MUU,G0
   !抗拉下降段参数,裂面剪力折减系数
   real*8 A1,A2
   !Crack=1,一条裂缝,=2两条裂缝,AddLoad=1加载,=0,卸载
   integer*4 Crack, AddLoad
   !裂缝角度
   real*8 ANGLE
   !t时的应力,主应力,应力增量,t+dt时的应力
   real*8 SIG(3),SIGP(3),dSIG(3),Stress(3)
   real*8 EPS(3),EPSP(3),dEPS(3),Strain(3)
   !非线性指标, 破坏面,最大非线性指标
   real*8 Beta,J2f,BetaMax
   !弹性本构矩阵,割线本构矩阵,本构矩阵
   real*8 De(3,3), Ds(3,3), D(3,3)
   !坐标转换矩阵
   real*8 N(3,3)
   integer(4) INC, NCycle
  end type typ_Concrete
  contains
  
  subroutine Con_Initial(C) !初始化混凝土参数
   type(typ_Concrete) :: C
   C%Fc=30d6; C%Ft=3d6; C%E0=30d9; C%MU0=0.2d0; C%MUU=0.2d0
   C%A1=3000; C%A2=0.5;
   C%G0=C%E0/(2.d0*(1.d0+C%MU0))
   C%Crack=0; C%Angle=0.d0; C%AddLoad=1
   call Con_Get_De(C)
   return
  end subroutine Con_Initial
  subroutine Con_Get_D(C)
   type(typ_Concrete) :: C
   call Con_Get_De(C)
   if(C%Crack<1) then
    call MAXMIN(C%SIG,C%SIGP,C%ANGLE)
   end if
!     RCM程序
!   if(C%Crack>1) then
!    call MAXMIN(C%EPS,C%EPSP,C%ANGLE)
!   end if
   call Con_Get_N(C)
   Call Con_Add_Load(C) !判断是否为加载
  
   if(C%AddLoad==0) then !如果是卸载
    call Con_UnLoad(C)
    !return
   end if
   if(C%AddLoad==1) then
    call Con_Get_Beta(C)
    if(C%Beta<=C%BetaMax) then
     call Con_UnLoad(C)
    else
     call Con_Get_Ds(C)
     C%BetaMax=C%Beta
    end if
    if(C%Crack<1) then
     call MAXMIN(C%Stress,C%SIGP,C%ANGLE)
    end if
!     RCM程序
!    if(C%Crack>1) then
!     call MAXMIN(C%EPS,C%EPSP,C%ANGLE)
!    end if
    call Con_Get_N(C)
    call Con_Crack(C)
   end if
   return
  end subroutine Con_Get_D
  subroutine Con_Crack(C) !处理裂缝
   type(typ_Concrete) :: C
   real*8 EPSC,EPST
   integer(4) :: CrackState(3)
   real*8 :: E1,E2,E12,G
  
   EPSC=-C%Fc/C%E0*2.d0 !峰值压应变
   EPST=C%Ft/C%E0       !峰值拉应变
   CrackState=0
  
   C%SIGP=matmul(transpose(C%N),C%Stress)
   if(C%SIGP(1)>C%Ft.and.C%Crack<1) then
    C%Crack=1
   end if
   if(C%SIGP(2)>C%Ft.and.C%Crack<2) then
    C%Crack=2
   end if
   if(C%Crack>0) then
    C%EPSP=matmul(transpose(C%N),(C%EPS+C%dEPS))
    if(C%EPSP(1)<=0.d0) then
     if(abs(C%EPSP(1))<abs(EPSC)) then
      C%SIGP(1)=2.d0*(C%EPSP(1)/EPSC)
  1       -(C%EPSP(1)/EPSC)**2
      C%SIGP(1)=-C%SIGP(1)*C%Fc
     else
      C%SIGP(1)=-C%Fc
     end if
    else
     if(C%EPSP(1)<=EPSt) then
      C%SIGP(1)=C%EPSP(1)*C%E0
     else
      C%SIGP(1)=C%Ft*exp(-C%A1*(C%EPSP(1)-EPSt))
      CrackState(1)=1
     end if
    end if
    if(C%EPSP(2)<=0.d0) then
     if(abs(C%EPSP(2))<abs(EPSC)) then
      C%SIGP(2)=2.d0*(C%EPSP(2)/EPSC)
  1       -(C%EPSP(2)/EPSC)**2
      C%SIGP(2)=-C%SIGP(2)*C%Fc
     else
      C%SIGP(2)=-C%Fc
     end if
    else
     if(C%EPSP(2)<=EPSt) then
      C%SIGP(2)=C%EPSP(2)*C%E0
     else
      C%SIGP(2)=C%Ft*exp(-C%A1*(C%EPSP(2)-EPSt))
      CrackState(2)=1
     end if
    end if
    C%SIGP(3)=C%G0*C%EPSP(3)*C%A2
  
    C%Stress=matmul(matinv(transpose(C%N)),C%SIGP)
    C%Strain=C%EPS+C%dEPS
    if(CrackState(1)==1) then
     E1=-0.01*C%E0
    else
     E1=C%E0
    end if
  
    if(CrackState(2)==1) then
     E2=-0.01*C%E0
    else
     E2=C%E0
    end if
  
    E12=0;
    G=C%G0*C%A2
  
    C%D(1,=(/E1,E12,0.d0/)
    C%D(2,:)=(/E12,E2,0.d0/)
    C%D(3,:)=(/0.0d0,0.0d0,G/)
    C%D=matmul(C%N,matmul(C%D,transpose(C%N)))
   end if
   return
  end subroutine Con_Crack
  subroutine Con_Get_Ds(C) !得到割线模量
   type(typ_Concrete) :: C
   real*8 Es, MUs
   if(C%Beta<=1.d0)  then
    Es=C%E0*(1.d0+sqrt(1.d0-C%Beta))/2.d0
    MUs=C%MU0
    if(C%Beta>0.8d0) then
     MUs=C%MUU-(C%MUU-C%MU0)*
  1     sqrt(1.d0-((C%Beta-0.8d0)/0.2d0)**2)
    end if
    C%Ds(1,:)=(/1.d0,MUs,0.0d0/)
    C%Ds(2,:)=(/MUs,1.d0,0.0d0/)
    C%Ds(3,:)=(/0.d0,0.d0,(1.d0-MUs)/2.d0/)
    C%Ds=C%Ds*Es/(1.d0-MUs**2)
    C%Stress=matmul(C%Ds,(C%EPS+C%dEPS))
    C%Strain=C%EPS+C%dEPS
    C%D=C%De
   else
    C%D=0.d0
    C%Stress=C%SIG
    C%Strain=C%EPS+C%dEPS
   end if
   return
  end subroutine Con_Get_Ds
  subroutine Con_Get_Beta(C) !得到非线性指标,
   !过程参见<<钢筋混凝土结构非线性有限元分析>&gt56
   type(Typ_Concrete) :: C
   real*8 SIGMA(6),S(6)
   real*8 I1,J2,J3,r,sita
   real*8 S_P(3)
   real*8 PI
   real*8 A,B,C1
   PI=atan(1.d0)*4.d0
   SIGMA=0.d0
   SIGMA(1:2)=C%SIG(1:2)+C%dSIG(1:2)/2.d0
   SIGMA(4)=C%SIG(3)+C%dSIG(3)/2.d0
   I1=SIGMA(1)+SIGMA(2)+SIGMA(3)
   S=SIGMA
   S(1)=S(1)-I1/3.d0
   S(2)=S(2)-I1/3.d0
   S(3)=S(3)-I1/3.d0
   J2=-S(1)*S(2)-S(2)*S(3)-S(3)*S(1)+S(4)**2+S(5)**2+S(6)**2
   J3=S(1)*S(2)*S(3)+2.d0*S(4)*S(5)*S(6)-S(1)*S(5)**2-S(2)
  1   *S(6)**2-S(3)*S(4)**2
   r=sqrt(4.d0*J2/3.d0)
   if(r.ne.0.d0) then
    sita=acos(4.d0*J3/r**3)/3.d0
   else
    sita=0.d0
   end if
   S_P(1)=2.d0*sqrt(J2)/sqrt(3.d0)*cos(sita)+I1/3.d0
   S_P(2)=2.d0*sqrt(J2)/sqrt(3.d0)*cos(sita-2.0d0*PI/3.d0)
  1   +I1/3.d0
   S_P(3)=2.d0*sqrt(J2)/sqrt(3.d0)*cos(sita+2.0d0*PI/3.d0)
  1   +I1/3.d0
   A=1.8148d0/C%Fc**2
   B=(1.180d0+13.2566d0*Cos(sita))/C%Fc
   C1=4.1145d0*I1/C%Fc-1.d0
   C%J2f=((-B+sqrt(B**2-4.d0*A*C1))/(2.d0*A))**2
   C%Beta=sqrt(J2)/sqrt(C%J2f)
   return
  end subroutine Con_Get_Beta
  
  subroutine Con_UnLoad(C) !卸载
   type(typ_Concrete) :: C
   C%D=C%De
   C%Stress=C%SIG+matmul(C%De,C%dEPS)
   C%Strain=C%EPS+C%dEPS
   return
  end subroutine Con_UnLoad
  subroutine Con_Add_Load(C) !判断加卸载
   type(typ_Concrete) :: C
   real*8 X(3),XP(3),J0,J1
   C%dSIG=matmul(C%De,C%dEPS)
   C%SIGP=matmul(transpose(C%N),C%SIG)
   X=C%SIG+C%dSIG
   XP=matmul(transpose(C%N),X)
   J0=(C%SIGP(1)-C%SIGP(2))**2+C%SIGP(2)**2+C%SIGP(1)**2
   J1=(XP(1)-XP(2))**2+XP(1)**2+XP(2)**2
   if(J0<=J1) then
    C%AddLoad=1
   else
    C%AddLoad=0
   end if
   return
  end subroutine Con_Add_Load
  subroutine Con_Get_N(C) !得到坐标转换矩阵
   type(typ_Concrete) :: C
   real*8 :: SinA,COSA
   COSA=cos(C%Angle); SINA=sin(C%Angle)
   C%N(1,:)=(/COSA**2,SINA**2,SINA*COSA/);
   C%N(2,:)=(/SINA**2,COSA**2,-SINA*COSA/);
   C%N(3,:)=(/-2d0*COSA*SINA,2.0d0*SINA*COSA,
  1  COSA**2-SINA**2/);
   return
  end subroutine Con_Get_N
  subroutine Con_Get_De(C) !得到弹性本构矩阵
   type(typ_Concrete) :: C
   C%G0=C%E0/(2.d0*(1.d0+C%MU0))
   C%De(1,:)=(/1.d0,C%MU0,0.d0/)
   C%De(2,:)=(/C%MU0,1.d0,0.d0/)
   C%De(3,:)=(/0.d0,0.d0,(1.d0-C%MU0)/2.d0/)
   C%De=C%De*C%E0/(1.d0-2.d0*C%MU0**2)
   return
  end subroutine Con_Get_De
  SUBROUTINE MAXMIN (STRESS,P,AG) !得到主应力(应变方向)
   implicit real*8 (A-H,O-Z)
   real*8 STRESS(3),P(3) !标量, 主方向
   real*8 T(3,3)  !转换矩阵
   PI=atan(1.0d0)*4.0d0 !得到PI
   CC = (STRESS(1)+STRESS(2)) * 0.5
   BB = (STRESS(1)-STRESS(2)) * 0.5
   CR =  SQRT(BB**2 + STRESS(3)**2)
   AG=PI/4.d0
   IF(BB.NE.0.0d0) Then
    AG = 0.5d0* ATAN2(-STRESS(3),BB)
   end if
   SINA=SIN(AG); COSA=COS(AG)
   T(1,:)=(/COSA**2,SINA**2,SINA*COSA/);
   T(2,:)=(/SINA**2,COSA**2,-SINA*COSA/);
   T(3,:)=(/-2d0*COSA*SINA,2.0d0*SINA*COSA,COSA**2-SINA**2/);
   P=matmul(transpose(T),STRESS)
   if(P(1)<P(2)) then
    CR=P1; P1=P2; P2=CR;
    AG=PI/2+AG;
   end if
   if(P(1)==0.0.and.P(2)==0.0) then
    AG=0;
   end if
       RETURN
  end subroutine MAXMIN
  function matinv(A) result (B)
  real(8) ,intent (in)::A(:,:)
  !real(8) , allocatable::B(:,:)
  real(8) , pointer::B(:,:)
  integer(4):: N,I,J,K
  real(8):,T
  real(8), allocatable::IS(:),JS(:)
  N=size(A,dim=2)
  allocate(B(N,N))
  allocate(IS(N));allocate(JS(N))
  B=A
  do  K=1,N
   D=0.0D0
   do I=K,N
    do J=K,N
     if(abs(B(I,J))>D) then
      D=abs(B(I,J))
      IS(K)=I
      JS(K)=J
     end if
    end do
   end do
   do J=1,N
    T=B(K,J)
    B(K,J)=B(int(IS(K)),J)
    B(int(IS(K)),J)=T
   end do
   do I=1,N
    T=B(I,K)
    B(I,K)=B(I,int(JS(K)))
    B(I,JS(K))=T
   end do
   B(K,K)=1/B(K,K)
   do J=1,N
    if(J.NE.K) then
     B(K,J)=B(K,J)*B(K,K)
    end if
   end do
   do I=1,N
    if(I.NE.K) then
     do J=1,N
      if(J.NE.K) then
       B(I,J)=B(I,J)-B(I,K)*B(K,J)
      end if
     end do
    end if
   end do
   do I=1,N
    if(I.NE.K) then
     B(I,K)=-B(I,K)*B(K,K)
    end if
   end do
  end do
  do K=N,1,-1
   do J=1,N
    T=B(K,J)
    B(K,J)=B(int(JS(K)),J)
    B(int(JS(K)),J)=T
   end do
   do I=1,N
    T=B(I,K)
    B(I,K)=B(I,int(IS(K)))
    B(I,int(IS(K)))=T
   end do
  end do
  return
  end function matinv
  end module MConcrete

zyx

可惜积分不够，看不见！
哪位好心人Email给我：
puppy@zju.edu.cn
不胜感激！

volunteer

为什么要设置积分，这个很没意思，象我这样的初学者就有很多看不到的，很麻烦

nieqi

这是别人写的吧

magiconly_2010

能不能传个pdf 的清楚额版本啊，谢谢

lipup

学习哈      

wanghui546

这个是什么啊，好像没有程序啊

catking2000

学习下，还得靠自己摸索啊

zh123168

可能值得看下，但是太乱啦！我认为如果刚贴到本论坛，就那么多乱码的话，绝非是原创！