📄 reconfea2005.f90
字号:
do j=1,NFelem
EFJtno=G_EFJtno(:,j)
EFStif=0.0
tmp1=abs(FYB(1,j))
tmp2=pi*feltprop(1,FType(j))*feltprop(2,FType(j))
if (FState(j)==1) then
if (curvetype==1) then
EsL=9.81e2-11.5e4*tmp1+2.51e6* tmp1**2
else
Esl=(5.3e2-5.04e4*tmp1+17.6e5*tmp1**2-21.9e6*tmp1**3)/sqrt(Fcu*0.67/40.7)
end if
KH=Esl*tmp2
Kv=feltprop(4,FType(j))*feltprop(2,FType(j)) !kv＝Ec×Lm
else
KH=10.0
kv=feltprop(4,FType(j))*feltprop(2,FType(j)) !kv＝Ec×Lm
end if
c=cos(feltprop(3,FType(j))*pi/180)
s=sin(feltprop(3,FType(j))*pi/180)
call shape_flet_stif(c,s,kh,kv,EFStif)
EFdof=G_EFdof(:,j)
call shape_Astif(kcol,EFStif,EFdof,4)
end do
!集装总刚结束
do j=1,ntotv
if (kcol(j)%row(j)<1e-4) then
go to 100
endif
end do
G_load=Ub_G_load
Ub_G_load=0.0
call Varbandsolv(Kcol,G_load,Ntotv)
G_load(0)=0.0
All_dis=All_dis+G_load
write(11,*) "+++++++=="
write(11,*) all_dis(124)
do j=1,ntotv
if (abs(All_dis(j))>1000) then
goto 100
end if
end do
do j=1,NSelem !钢筋单元的分析
ESJtno=G_ESJtno(:,j)
ESdof=G_ESdof(:,j)
ESdis(:)=G_load(ESdof)
ESCoord=G_Pcoord(:,ESJtno)
if (SState(j)==1) then !弹性阶段
Es=steprop(1,SType(j))
else
Es=10 !水平段了
end if
c=ESCoord(1,2)-ESCoord(1,1)
s=ESCoord(2,2)-ESCoord(2,1)
L=sqrt(c*c+s*s)
c=c/L
s=s/l
SYL(j)=SYL(j)+Es*(-c*ESdis(1)-s*ESdis(2)+c*ESdis(3)+s*ESdis(4))/L !存第二点轴应力,拉为正
tmp1=10
if (abs(SYL(j))<steprop(3,SType(j))) then
SState(j)=1
else
if (sstate(j)==1) then
destroynum=destroynum+1
end if
SState(j)=2
if (SYL(j)>0) then
tmp1=(SYL(j)-steprop(3,SType(j)))*steprop(2,SType(j))
SYL(j)=steprop(3,SType(j))
else
tmp1=(SYL(j)+steprop(3,SType(j)))*steprop(2,SType(j))
SYL(j)=-steprop(3,SType(j))
end if
ESF(1)=-tmp1*c
ESF(2)=-tmp1*s
ESF(3)=tmp1*c
ESF(4)=tmp1*s
Ub_G_load(ESdof)=Ub_G_load(ESdof)+ESF
continue
end if
end do
do j=1,NFelem !粘结单元
EFJtno=G_EFJtno(:,j)
EFdof=G_EFdof(:,j)
EFdis(:)=G_load(EFdof)
tmp1=abs(FYB(1,j))
tmp2=pi*feltprop(1,FType(j))*feltprop(2,FType(j))
if (FState(j)==1) then
if (curvetype==1) then
EsL=9.81e2-11.5e4*tmp1+2.51e6* tmp1**2
else
Esl=(5.3e2-5.04e4*tmp1+17.6e5*tmp1**2-21.9e6*tmp1**3)/sqrt(Fcu*0.67/40.7)
end if
KH=Esl*tmp2
Kv=feltprop(4,FType(j))*feltprop(2,FType(j)) !kv＝Ec×Lm
else
KH=10.0
kv=feltprop(4,FType(j))*feltprop(2,FType(j)) !kv＝Ec×Lm
end if
c=cos(feltprop(3,FType(j))*pi/180)
s=sin(feltprop(3,FType(j))*pi/180)
EFYB(1)=-c*EFdis(1)-s*EFdis(2)+c*EFdis(3)+s*EFdis(4) !增量
EFYB(2)=s*EFdis(1)-c*EFdis(2)-s*EFdis(3)+c*EFdis(4)
FYB(1,j)=FYB(1,j)+EFYB(1) !粘结单元的应变H或称位移 总的
FYB(2,j)=FYB(2,j)+EFYB(2) !粘结单元的应变v
!此fyb是某个荷载步的应变，不是总应变,这样才是总应变
if (Fstate(j)==1) then
FForec(1,j)=FForec(1,j)+KH*EFYB(1) !!粘结单元的力H
FForec(2,j)=FForec(2,j)+Kv*EFYB(2) !!粘结单元的力V
FYL(j)=FForec(1,j)/tmp2 !H方向的应力
end if
!以下是粘结单元的分析
if (curvetype==1) then
s0=0.011
tmax=4.96
else
s0=0.03
tmax=4.64*sqrt(Fcu*0.67/40.7)
end if
tmp1=abs(FYB(1,j)) ! 刷新应变
if (tmp1<s0) then
FState(j)=1
if(curvetype==1) then
R_FYL(j)=9.81e2*tmp1-5.74e4*tmp1**2+0.837e6*tmp1**3
else
R_FYL(j)=(5.3e2*tmp1-2.52e4*tmp1**2+5.87e5*tmp1**3-5.47e6*tmp1**4)*sqrt(Fcu*0.67/40.7)
end if
if (FYB(1,j)<0) R_FYL(j)=-R_FYL(j)
tmp3=FForec(1,j)-R_FYL(j)*tmp2
FForec(1,j)=R_FYL(j)*tmp2
FYL(j)=R_FYL(j)
else
if (fstate(j)==1) then
destroynum=destroynum+1
end if
FState(j)=2
if (FYB(1,j)<0) Tmax=-tmax
tmp3=FForec(1,j)-tmax*tmp2
FForec(1,j)=tmax*tmp2
FYL(j)=tmax
end if
EFF(1)=-tmp3*c
EFF(2)=-tmp3*s
EFF(3)=tmp3*c
EFF(4)=tmp3*s
Ub_G_load(EFdof)=Ub_G_load(EFdof)+EFF
end do
do j=1,NCelem !混凝土单元单元分析
if (j==9) then
continue
end if
ECJtno=G_ECJtno(:,j)
Ecdof=G_Ecdof(:,j)
ECdis(:)=G_load(Ecdof)
Eccoord=G_Pcoord(:,ECJtno)
call shape_T3_B(Eccoord,Bmatx,area)
ECYB=matmul(Bmatx,ECdis) !单元当前步的应变
CYB(:,1,j)=CYB(:,1,j)+ECYB !所有单元当前步的应变
tmp1=0.5*sqrt((CYB(1,1,j)-CYB(2,1,j))**2+CYB(3,1,j)**2)
C_main_YB(1,1,j)=(CYB(1,1,j)+CYB(2,1,j))*0.5+tmp1
C_main_YB(2,1,j)=(CYB(1,1,j)+CYB(2,1,j))*0.5-tmp1 !主应变
thera=CYL(6,1,j)
call D_matrix_2d_plane(type2,Dmatx,e(1,1,j),e(2,1,j),pu(1,j))
call shape_YL_translate_matrix(thera,T)
Dmatx=matmul(matmul(T,dmatx),transpose(T))
ECYL=matmul(Dmatx,ECYB)
CYL(1:3,1,j)=CYL(1:3,1,j)+ECYL !累加的应力
tmp1=CYL(1,1,j)+CYL(2,1,j)
tmp2=sqrt((CYL(1,1,j)-CYL(2,1,j))**2+4*CYL(3,1,j)**2)
CYL(4,1,j)=(tmp1+tmp2)*0.5
CYL(5,1,j)=(tmp1-tmp2)*0.5
Tx=CState(1,j)/10
Ty=CState(1,j)-10*Tx
oldstate=cstate(1,j)
if(Tx>=2 .or. Ty>=2) then
angle=cyl(6,1,j)
goto 2000 !非弹性状态11
end if
tmp3=CYL(2,1,j)-CYL(5,1,j)
if (abs(tmp3)<1e-4) then
angle=0.0
else
angle=90-57.296*atan(CYL(3,1,j)/tmp3)
end if !求角度
2000 tmp1= ECYL(1)+ECYL(2)
tmp2=sqrt((ECYL(1)-ECYL(2))**2+4*ECYL(3)**2)
tmp3=(tmp1+tmp2)*0.5
tmp4=(tmp1-tmp2)*0.5
if ( abs(cyl(6,1,j)-angle)<45 .or. abs(cyl(6,1,j)-angle+180)<45 .or. abs(cyl(6,1,j)-angle-180)<45) then
siu(1,1,j)=siu(1,1,j)+tmp3/E(1,1,j)
siu(2,1,j)=siu(2,1,j)+tmp4/E(2,1,j) !等效单向应变的累加
else
tmp5=siu(1,1,j)
siu(1,1,j)=siu(2,1,j)+tmp3/e(1,1,j)
siu(2,1,j)=tmp5+tmp4/e(2,1,j)
end if
cyl(6,1,j)=angle
if (CState(1,j)==11) then
call kupf(Fcu,CYL(:,:,j),Fz(:,:,j)) !算破坏应力
end if
fc=-Fcu*0.667
if (CYL(5,1,j)<0.0 .and. CYL(4,1,j)<0.0) then !双压
if (abs(CYL(5,1,j))<abs(fz(2,1,j))) then
call saenz(1,Fz(:,:,j),Fc,sic(:,:,J),siu(:,:,J),E(:,:,j),R_CYL(:,:,J),Eo,cyl(:,1,j)) !算E和真实应力
CState(1,j)=11
pu(1,j)=0.2
sl(1)=R_CYL(1,1,j)
sl(2)=R_CYL(2,1,j)
else
do k=1,2
if (abs(siu(k,1,j))<abs(4*sic(k,1,j)) .and. abs(siu(k,1,j))>=abs(sic(k,1,j))) then
sl(k)=fz(k,1,j)-(0.8*fz(k,1,j)*(siu(k,1,j)-sic(k,1,j)))/(3.0*sic(k,1,j))
CState(1,j)=22
else if (abs(siu(k,1,j))>abs(4*sic(k,1,j))) then
sl(k)=0.0
CState(1,j)=33
end if
enddo
E(:,1,j)=10
pu(1,j)=0.2
end if
elseif (CYL(4,1,j)>=0.0 .and. CYL(5,1,j)<0.0) then !拉压
if (abs(CYL(4,1,j))<abs(fZ(1,1,j))) then
call saenz(2,Fz(:,:,j),Fc,sic(:,:,J),siu(:,:,J),E(:,:,j),R_CYL(:,:,J),Eo,cyl(:,1,j))
R_CYL(1,1,j)=CYL(4,1,j)
sl(1)=R_CYL(1,1,j)
sl(2)=R_CYL(2,1,j)
CState(1,j)=11
pu(1,j)=0.2
else
Eiu=abs(siu(1,1,j))
Etu=0.00015
if (eiu<etu .and. eiu>=abs(Fz(1,1,j)/Eo)) then
sl(1)=fz(1,1,j)
tx=2
else if(eiu>=etu) then
sl(1)=0.0
tx=3
end if
E(1,1,j)=10
pu(1,j)=0.2
!以下是受压
Eiu=abs(siu(2,1,j))
Ecu=0.0033
Ec0=0.002
if (eiu<Ec0) then
Es=abs(Fc/Ec0)
q=eiu/Ec0
s=1+(Eo/es-2.0)*q+q**2
E(2,1,j)=eo*(1.0-q**2)/s**2
sl(2)=-eo*eiu/s
ty=1
elseif (eiu>=ec0 .and. eiu<ecu)then
e(2,1,j)=10
sl(2)=fc-0.8*fc*(eiu-ec0)/(ecu-ec0)
ty=2
else
sl(2)=0.0
e(2,1,j)=10
ty=3
end if
end if
pu(1,j)=0.2
CState(1,j)=10*tx+ty
elseif (CYL(4,1,j)>=0.0 .and. CYL(5,1,j)>=0.0) then !双拉
if (abs(CYL(4,1,j))<abs(fz(1,1,j))) then
sl(1)=CYL(4,1,j)
sl(2)=CYL(5,1,j)
CState(1,j)=11
else
eiu=abs(siu(1,1,j))
etu=0.00015
if (eiu>=abs(fz(1,1,j)/Eo) .and. eiu<etu) then
sl(1)=fz(1,1,j)
tx=2
elseif (eiu>=etu) then
sl(1)=0.0
tx=3
end if
e(1,1,j)=0.0
pu(1,j)=0.2
!小主应力方向
eiu=abs(siu(2,1,j))
if (abs(CYL(5,1,j))<abs(fz(2,1,j)) .or. eiu<abs(fz(2,1,j)/Eo)) then
sl(2)=CYL(5,1,j)
ty=1
e(2,1,j)=eo !这句不写我觉得更好
pu(1,j)=0.2
elseif (eiu>=abs(fz(2,1,j)/eo) .and. eiu<etu) then
sL(2)=fz(2,1,j)
ty=2
E(2,1,j)=10
pu(1,j)=0.2
else
sl(2)=0
ty=3
E(2,1,j)=10
pu(1,j)=0.2
end if
CState(1,j)=10*tx+ty
end if
end if
if (oldstate/=cstate(1,j)) then
destroynum=destroynum+1
end if
thera=cyl(6,1,j)
call shape_YL_translate_matrix(thera,T)
yl4(1)=CYL(4,1,j)-sl(1)
yl4(2)=CYL(5,1,j)-sl(2)
yl4(3)=0.0
yl1=matmul(T,yl4)
call shape_T3_B(Eccoord,Bmatx,area)
yl3=area*Con_t*matmul(transpose(bmatx),yl1)
if (max(abs(maxval(yl3)),abs(minval(yl3)))>20) then
continue
end if
Ub_G_load(Ecdof)=Ub_G_load(Ecdof)+yl3
yl4(1)=sl(1)
yl4(2)=sl(2)
yl4(3)=0.0
yl1=matmul(T,yl4)
CYL(1:3,1,j)=yl1
CYL(4,1,j)=sl(1)
CYL(5,1,j)=sl(2) !这里为什么不进一步算角度 我求
IF (CState(1,j)==11)THEN
tmp3=CYL(2,1,j)-CYL(5,1,j)
if (abs(tmp3)<1e-4) then
angle=0.0
else
angle=90-57.296*atan(CYL(3,1,j)/tmp3)
end if
CYL(6,1,J)=ANGLE
END IF
end do
goto 200
end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
100 continue
end do
end program RcConFea2005
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -