module TypeDef !定义混凝土材性模块
type :: typ_Concrete
real*8 fc, ft,E0, ENU;
C 抗压强度+，抗拉强度+，初始弹性模量，初始泊松比
real*8 Es, Et, ENUs !割线模量，切线模量，割线泊松比
real*8 SIG(6), EPS(6),dEPS(6); !开始时应力，应变，应变增量
real*8 Stress(6), Strain(6) !结束时应力，应变
real*8 Dt(6,6), Dela(6,6), Ds(6,6)
C 切线刚度矩阵，弹性刚度矩阵，割线刚度矩阵
end type typ_Concrete
end module TypeDef
subroutine HYPELA(D,G,E,DE,S,TEMP0,
1 DTEMP,NGENS,N,NN,KC,MATS,NDI,NSHEAR)
! D(6x6) 迭代本构矩阵(out)
! G(6) 由于状态改变引起的应力变化，不用(out)
! E(6) 开始时刻的应变(in)
! DE(6) 应变增量(in)
! S(6) 开始时刻的应变(in & out)
! Temp0 温度(in)
! DTEMP 温度变化(in)
! NGENS 应变维数(in)
! N(2) 单元编号(in)
! NN 积分点编号(in)
! KC 层号(in)
! MATS 材料编号(in)
! NDI 正应力维数(in)
! NSHEAR 剪应力维数(in)
use typedef
implicit none
integer :: ngens,nn,kc,mats,ndi,nshear
real*8 :: e(ngens),de(ngens),temp0(1),dtemp(1),g(ngens)
1 ,d(ngens,ngens),s(ngens)
integer :: n(2)
type(typ_concrete) :: C
real*8 Beta1
real*8 s_m,J2,J3,r,sita,TempA, TempB, TempC;
c 赋值
C%fc=30.; C%ft=3.; C%E0=30e3; C%ENU=0.18;
C%SIG=s; C%EPS=e; C%dEPS=de;
D=0.
c-----------------------------------------
c 计算弹性矩阵
call Get_D_Matrix(C%Dela,C%E0,C%ENU)
C 计算主应力和割线刚度
C%Stress=C%SIG
s_m=(C%Stress(1)+C%Stress(2)+C%Stress(3))/3. !计算平均应力
s(1:3)=C%Stress(1:3)-s_m
s(4:6)=C%Stress(4:6) !计算应力偏量
J2=-s(1)*s(2)-s(2)*s(3)-s(3)*s(1)+s(4)**2+s(5)*2+s(6)**2 !计算J2
J3=s(1)*s(2)*s(3)+2.*s(4)*s(5)*s(6) ! 计算J3
1 -s(1)*s(5)**2-s(2)*s(6)**2-s(3)*s(4)**2
r=sqrt(4.*J2/3.)
if(r.ne.0.) then
sita=acos(4.*J3/r**3)/3.
if(4.*J3/r**3>1.d0) sita=acos(1.d0)/3.
if(4.*J3/r**3<-1.d0) sita=acos(-1.d0)/3.
else
sita=0.
end if
!根据江见鲸模型求解Beta
if(J2>0.) then
TempA=1.2856/C%fc**2;
TempB=(1.4268+10.2551*cos(sita))/C%fc;
TempC=3.2128*s_m*3./C%fc-1.;
Beta1=-TempB+sqrt(TempB**2-4.*TempA*TempC)
Beta1=Beta1/2./TempA
Beta1=sqrt(J2)/Beta1
else
Beta1=0.
end if
if(Beta1>1.) Beta1=1.
C%Es=C%E0/2.*(1.+sqrt(1.-Beta1)) ! 计算割线模量
C%ENUs=C%ENU ! 计算割线泊松比
call Get_D_Matrix(C%Ds,C%Es,C%ENUs) ! 计算割线刚度矩阵
C%Et=sqrt(1.-Beta1)*C%E0 ! 计算切线刚度
call Get_D_Matrix(C%Dt,C%Et,C%ENU) ! 计算切线刚度矩阵
! 计算当前应力
if(Beta1<1.d0) then
C%Strain=C%EPS+C%dEPS;
C%Stress=matmul(C%Ds,C%Strain);
else
C%Dt=1.d-10*C%Dela
end if
C-----------------------------------------
D=C%Dt ! 设置迭代刚度矩阵
s=C%Stress ! 输出应力
return
end subroutine
subroutine Get_D_Matrix(D,E,MU) !各向同性刚度矩阵
implicit none
real*8 :: D(6,6), E,MU
integer :: K1, K2
D=0.
do K1=1, 3
do K2=1, 3
D(K2,K1)=MU
end do
D(K1,K1)=1.-MU
END do
do K1=4,6
D(K1,K1)=(1.-2.*MU)*0.5
end do
D=D*E/(1.+MU)/(1.-2.*MU)
return
end subroutine