一个混凝土本构的程序

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) !得到非线性指标,
   !过程参见<<钢筋混凝土结构非线性有限元分析>>P56
   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),B
   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
    R=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 (
  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:S(,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