非线性弹性全量本构模型程序(THURC3A)

module TypeDef !定义混凝土材性模块
type :: typ_Concrete
real*8 fc, ft,E0, ENU, EPS_Crush !抗压强度+，抗拉强度+，初始弹性模量，初始泊松比，压碎应变-
real*8 Es, ENUs !割线模量，割线泊松比
real*8 T(6,6) !坐标转换矩阵
integer NCrack (3), Pre_NCrack(3), Pre_inc, Pre_incsub; !开裂记录，前次迭代开裂记录，前次增量步，前次增量子步
real*8 SIG(6), EPS(6),dEPS(6); !开始时应力，应变，应变增量
real*8 StressP(6), StrainP(6); !主应力，主应变
real*8 Stress(6), Strain(6) !结束时应力，应变
real*8 Beta,Pre_Beta !非线性指标，前次迭代非线性指标
real*8 D(6,6), Dela(6,6), Ds(6,6) !刚度矩阵，弹性刚度矩阵，割线刚度矩阵
end type typ_Concrete
end module TypeDef
module My_MOD !开辟公共变量空间
use TypeDef
type(typ_Concrete) :: My_Con(1000,8) !定义混凝土数组
end module


subroutine Get_DS(D,G,E,DE,S,TEMP0,
1 DTEMP,NGENS,N,NN,KC,MATS,NDI,NSHEAR,inc,incsub,ncycle)
! 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)
! inc 当前增量步(in)
! incsub 当前增量子步(in)
! ncycle 当前循环数(in)

use IMSL !引用IMSL函数库
use typedef
use My_Mod
implicit none
integer :: ngens,nn,kc,mats,ndi,nshear,inc,incsub,ncycle
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,strain_m
real*8 s_m,J2,J3,r,sita,TempA, TempB, TempC;
integer NSubStep !子步积分步数
integer I,J,K1,K2


C=My_Con(n(1),nn) !得到内存中保留的数据

c ----------------------------------------------------
c 初次计算，清零并赋值
if(inc==0.and.incsub==0.and.ncycle==0) then
open(77,file='debug.txt')
write(77,*)
close(77)

C%fc=30.; C%ft=3.; C%E0=30e3; C%ENU=0.18;
C%EPS_Crush=-0.0033;
C%T=0.; C%NCrack=0; C%Pre_NCrack=0;
C%Beta=0; C%Pre_Beta=0;
C%Pre_inc=0;
C%Pre_incsub=0;
C%SIG=0.; C%EPS=0.; C%Stress=0.; C%Strain=0.;
end if
c ----------------------------------------------------
c 如果新的增量步开始，则更新相应变量
if(inc>C%Pre_inc .or. incsub>C%Pre_incsub) then
C%Pre_inc=inc; C%Pre_incsub=incsub
C%NCrack=C%Pre_NCrack; !修正裂缝状态
C%Beta=C%Pre_Beta; !修正非线性指标状态
! 判断是否压坏
strain_m=(C%EPS(1)+C%EPS(2)+C%EPS(3))/3.
if(Strain_m>0.) Strain_m=0.
if(minval(C%EPS(1:3)-Strain_m)<C%EPS_Crush) then
C%NCrack=100 !彻底破坏
end if
end if

c 数据赋值
open(77,file='debug.txt',position='append')

C%SIG=s; C%EPS=e; C%dEPS=de;
C%Pre_NCrack=C%NCrack
C%Pre_Beta=C%Beta
NSubStep=4
c ------------------------------------
c 计算弹性矩阵
C%Dela=0.
do K1=1, 3
do K2=1, 3
C%Dela(K2,K1)=C%ENU
end do
C%Dela(K1,K1)=1.-C%ENU
END do
do K1=4,6
C%Dela(K1,K1)=(1.-2.*C%ENU)*0.5
end do
C%Dela=C%Dela*C%E0/(1.+C%ENU)/(1.-2.*C%ENU)
c --------------------------------------
c 如果已经压碎，应力清零，刚度为很小值，结束计算
if(maxval(C%NCrack)==100) then
C%D=0.0001*C%Dela
C%SIG=0.;
C%Stress=0.;
s=0.
return
end if

C 计算主应力和割线刚度
C%Stress=C%SIG
do I=1, NSubStep

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.
else
sita=0.
end if

if(maxval(abs(C%Pre_NCrack))==0) then !没有裂缝
call Get_T_Matrix(C%SIG,C%T) !计算坐标转换矩阵
end if

C%StressP=matmul(transpose(C%T),C%Stress); !计算主应力
C%StrainP=matmul(transpose(C%T),C%Strain); !计算主应变

if(C%StressP(1)<0.05*C%fc.and.
1 maxval(abs(C%Pre_NCrack))==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 !根据江见鲸模型求解Beta

if(Beta1>C%Beta) then !Beta应该始终增大(对于全量模型)
C%Pre_Beta=Beta1
else
Beta1=C%Beta
end if
! 计算割线刚度和泊松比
if(Beta1>1.) Beta1=1.
C%Es=C%E0/2.*(1.+sqrt(1.-Beta1))
if(Beta1<0.8) C%ENUs=C%ENU
if(Beta1.ge.0.8) C%ENUs=0.42-(0.42-C%ENU)*
1 sqrt(1.-((Beta1-.8)/.2)**2)
!计算割线刚度矩阵
C%Ds=0.
do K1=1, 3
do K2=1, 3
C%Ds(K2,K1)=C%ENUs
end do
C%Ds(K1,K1)=1.-C%ENUs
end do
do K1=4,6
C%Ds(K1,K1)=(1.-2.*C%ENUs)*0.5
end do
C%Ds=C%Ds*C%Es/(1.+C%ENUs)/(1.-2.*C%ENUs)
else !如果处于开裂控制区
C%Ds=C%Dela
do K1=1,3
if(C%StressP(K1)>C%ft .OR. C%Pre_NCrack(K1)>0) then !按开裂处理
C%Pre_NCrack(K1)=1;
Call Crack_Open(C,K1) !计算开裂矩阵
end if
end do

C%Ds=matmul(C%T,matmul(C%Ds,transpose(C%T))); !计算割线刚度矩阵

end if
if(Beta1<0.99999d0) then !如果没有达到极限应力
C%Strain=C%EPS+C%dEPS*real(I/NSubStep);
C%Stress=matmul(C%Ds,C%Strain);

else !达到极限应力后应力不变
C%Stress=C%SIG
C%Ds=1.d-6*C%Ds
end if
end do

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 设置迭代刚度矩阵
D=c%Ds+1.d-6*C%Dela
s=C%Stress
My_Con(N(1),nn)=C;

c
c write(77,*) 'Step: ',inc, incsub, ncycle
c write(77,*) Beta1
c write(77,*) C%Pre_NCrack
c write(77,*) C%Stress(1:3)
c write(77,*)
close(77)
return
end subroutine

C 根据开裂修正刚度矩阵
subroutine Crack_Open(C,I0)
use typeDef
implicit none
type (typ_Concrete) :: C
integer,intent(in) :: I0
real*8 FACTOR,ECRA1,ECRA2,ECRA12

C%Ds(I0,:)=0.
C%Ds(:,I0)=0.

ECRA1=C%StrainP(I0);
ECRA2=0.;
ECRA12=maxval(abs(C%StrainP(4:6)))
c 计算裂面剪力传递系数
call USHRET0 (FACTOR,ECRA1,ECRA2,ECRA12)
C%Ds(4,4)=FACTOR*C%Ds(4,4)
C%Ds(5,5)=FACTOR*C%Ds(5,5)
C%Ds(6,6)=FACTOR*C%Ds(6,6)

return
end subroutine


c 计算裂面剪力传递系数
SUBROUTINE USHRET0 (FACTOR,ECRA1,ECRA2,ECRA12)
IMPLICIT REAL *8 (A-H, O-Z)
factor=0.4;
RETURN
END

c 计算主应力及坐标转换矩阵
subroutine Get_T_Matrix(olds,T)
use IMSL
implicit none
real*8 olds(6), T(6,6)
real*8 SIG(3,3),eval_r(3), EVEC(3,3)
real*8 l(3),m(3),n(3)
real*8 SIGP(3)
integer I
SIG(1,1)=olds(1)
SIG(2,2)=olds(2)
SIG(3,3)=olds(3)
SIG(1,2)=olds(4); SIG(2,1)=olds(4);
SIG(2,3)=olds(5); SIG(3,2)=olds(5);
SIG(1,3)=olds(6); SIG(3,1)=olds(6);
call DEVCSF(3, SIG, 3, EVAL, EVEC, 3)
SIGP(1)=maxval(EVAL)
do I=1,3
if(eval_r(I)==SIGP(1)) then
l=EVec(:,I)
eval_r(I)=minval(EVAL)-10;
exit
end if
end do
SIGP(2)=maxval(EVAL)
do I=1,3
if(eval_r(I)==SIGP(2)) then
m=EVec(:,I)
eval_r(I)=minval(EVAL)-10;
exit
end if
end do
SIGP(3)=maxval(EVAL)
do I=1,3
if(eval_r(I)==SIGP(3)) then
n=EVec(:,I)
eval_r(I)=minval(EVAL)-10;
exit
end if
end do
do I=1,3
T(I,:)=(/l(I)**2,m(I)**2,n(I)**2,l(I)*m(I),
1 m(I)*n(I),n(I)*l(I)/)
end do
T(4,:)=(/2.d0*l(1)*l(2),2.d0*m(1)*m(2),2.d0*n(1)*n(2),
1 l(1)*m(2)+l(2)*m(1),m(1)*n(2)+m(2)*n(1),n(1)*l(2)+n(2)*l(1)/)
T(5,:)=(/2.d0*l(2)*l(3),2.d0*m(2)*m(3),2.d0*n(2)*n(3),
1 l(2)*m(3)+l(3)*m(2),m(2)*n(3)+m(3)*n(2),n(2)*l(3)+n(3)*l(2)/)
T(6,:)=(/2.d0*l(3)*l(1),2.d0*m(3)*m(1),2.d0*n(3)*n(1),
1 l(3)*m(1)+l(1)*m(3),m(3)*n(1)+m(1)*n(3),n(3)*l(1)+n(1)*l(3)/)
return
end subroutine


c marc 接口程序
SUBROUTINE HYPELA(D,G,E,DE,S,TEMP0,
1 DTEMP,NGENS,N,NN,KC,MATS,NDI,NSHEAR)
implicit real*8 (a-h,o-z)
INCLUDE '../common/concom' ! 通过concom模块得到当前的计算步数
integer :: ngens,nn,kc,mats,ndi,nshear
real*8 :: e(1),de(1),temp0(*),dtemp(*),g(1),d(ngens,ngens),s(1)
integer :: n(2)
if(mats==1) then !如果材料编号是1
call Get_DS(D,G,E,DE,S,TEMP0,
1 DTEMP,NGENS,N,NN,KC,MATS,NDI,NSHEAR,inc,incsub,ncycle)
end if
return
end subroutine

c 后处理子程序
subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,
* nshear,jpltcd)
c* * * * * *
c
c select a variable contour plotting (user subroutine).
c
c v variable
c s (idss) stress array
c sp stresses in preferred direction
c etot total strain (generalized)
c eplas total plastic strain
c ecreep total creep strain
c t current temperature
c m(1) user element number
c m(2) internal element number
c nn integration point number
c layer layer number
c ndi (3) number of direct stress components
c nshear (3) number of shear stress components
c
c* * * * * *
use My_Mod
implicit real*8 (a-h,o-z) dp
dimension s(*),etot(*),eplas(*),ecreep(*),sp(*),m(2)
type(typ_Concrete) :: C
C=My_Con(m(1),nn)
c 后处理变量1：输出是否开裂
if(jpltcd==1) then
if(C%NCrack(1).ne.0) then
v=1
else
v=0
end if
end if
c 后处理变量2-4：输出裂缝状态
do I=1,3
if(jpltcd==I+1) then
v=C%NCrack(I)
end if
end do
c 后处理变量5：输出非线性指标
if(jpltcd==5) then
v=C%Beta
end if

c 后处理变量8：输出裂缝数量
if(jpltcd==8) then
v=0.
do I=1,3
if(C%NCrack(I).ne.0) v=v+1
end do
end if

return
end