📄 be2die.for
字号:
else
WD=DSQRT((-UNod(1,NNSlv,NESlv)-UNod(1,NN,NE))**2
1 +(-UNod(2,NNSlv,NESlv)-UNod(2,NN,NE))**2
1 +(UNod(3,NNSlv,NESlv)-UNod(3,NN,NE))**2)
endif
endif
100 If(WD.gt.WDmgNod(NN,NE)) WDmgNod(NN,NE)=WD
*=========================================================
return
end
*************************************************************
*************************************************************
*************************************************************
Subroutine AdvCrkTip(NLoad,NE,
1 NCmp,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
1 JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
1 NtotCfKCap,NtotCfK,NthCfk,CfUP,
1 WDmgNod,UNod,PNod,UHomoNod,PHomoNod,
1 OCrkAdv,ACrkAdv)
Implicit Real*8 (A-H,O-Z)
Implicit Integer (I-N)
Dimension XEle(2,2,NECap),NodEle(NECap),OEle(NECap)
Dimension NthBdEle(NECap),NthGrpEle(NECap)
Dimension JCrkN0C1Tip2Ele(NECap),NthSlvMstEle(NECap)
Dimension XNod(2,3,NECap),ShpTNod(0:3,3,NECap)
Dimension NthCfK(NECap,3,NECap,3,0:1)
Dimension CfUP(NCmp,NCmp*2,0:NtotCfKCap)
Dimension WDmgNod(3,NECap),UNod(NCmp,3,NECap),PNod(NCmp,3,NECap)
Dimension UHomoNod(NCmp,3,NECap),PHomoNod(NCmp,3,NECap)
*=========================================================
!new tip element
X11=XEle(1,2,NE)
Y11=XEle(2,2,NE)
X12=X11+ACrkAdv*DCOS(OCrkAdv)
Y12=Y11+ACrkAdv*DSIN(OCrkAdv)
Do 277 ME=1,NtotEle
If(NE.eq.ME) goto 277
If(NthBdEle(ME).ne.NthBdEle(NE)) goto 277
X21=XEle(1,1,ME)
Y21=XEle(2,1,ME)
X22=XEle(1,2,ME)
Y22=XEle(2,2,ME)
call Intersect
1 (X11,Y11,X12,Y12,X21,Y21,X22,Y22,INTERSEC,IENDOF2)
If(INTERSEC.ne.0) then
IEINTERS=ME
goto 279
endif
277 continue
279 X1Stt=XEle(1,2,NE)
X2Stt=XEle(2,2,NE)
If(INTERSEC.eq.0) then
X1End=X1Stt+ACrkAdv*DCOS(OCrkAdv)
X2End=X2Stt+ACrkAdv*DSIN(OCrkAdv)
JBnd0Crk1Tip2=2
else
X1End=XEle(1,IENDOF2,IEINTERS)
X2End=XEle(2,IENDOF2,IEINTERS)
JBnd0Crk1Tip2=1
endif
NthBd=NthBdEle(NE)
NDiv=1
NodE=3
Gradient=1.d0
NthGrp=NthGrpEle(NE)
Call SegElemMst(NtotEle,NthBd,
1 JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE,
1 Gradient,NthGrp,
2 NECap,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
2 JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod)
Do I=1,3
WDmgNod(I,NtotEle)=
1 UserInitWDmg(I,NtotEle,NthGrpEle(NE),XNod(1,I,NtotEle))
enddo
Call HomoFldsEle(NLoad,
1 NCmp,NodEle(NtotEle),OEle(NtotEle),NthBdEle(NtotEle),
1 XNod(1,1,NtotEle),
2 UHomoNod(1,1,NtotEle),PHomoNod(1,1,NtotEle)) !output
!revision of old tip element
X1Stt=XEle(1,1,NE)
X2Stt=XEle(2,1,NE)
X1End=XEle(1,2,NE)
X2End=XEle(2,2,NE)
JBnd0Crk1Tip2=1
NthBd=NthBdEle(NE)
NDiv=1
NodE=3
Gradient=1.d0
NthGrp=NthGrpEle(NE)
NE_1=NE-1
Call SegElemMst(NE_1,NthBd,
1 JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE,
1 Gradient,NthGrp,
2 NECap,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
2 JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod)
!update influence coefficient
Do 310 NES=1,NtotEle
Do 310 NNS=1,NodEle(NES)
Do 310 NEF=1,NtotEle
Do 310 NNF=1,NodEle(NEF)
If(NthBdEle(NES).ne.NthBdEle(NEF)) goto 310
If(NES.ne.NE.and.NES.ne.NtotEle.
1 and.NEF.ne.NE.and.NEF.ne.NtotEle) goto 310
JDeriv=0
If(JCrkN0C1Tip2Ele(NES).gt.0) JDeriv=1
If(NthCfK(NES,NNS,NEF,NNF,0).ne.0) then
NthCfKRec=NthCfK(NES,NNS,NEF,NNF,0)
else
NthCfKRec=NtotCfK+1
endif
Call CfIntInf(
1 NCmp,NthBdEle(NES),XNod(1,NNS,NES),OEle(NES),
3 XEle(1,1,NEF),OEle(NEF),ShpTNod(0,NNF,NEF),
5 JDeriv,CfUP(1,1,NthCfKRec))
If(NthCfK(NES,NNS,NEF,NNF,0).ne.0) then
else
NtotCfK=NtotCfK+1
NthCfK(NES,NNS,NEF,NNF,0)=NtotCfK
If(JDeriv.eq.1) NtotCfK=NtotCfK+1
If(JDeriv.eq.1) NthCfK(NES,NNS,NEF,NNF,1)=NtotCfK
endif
310 continue
*==============================================================
return
end
*******************SUBROUTINE INTERSECT*********************
SUBROUTINE INTERSECT(X11,Y11,X12,Y12,X21,Y21,X22,Y22,
1 INTERSEC,IENDOF2)
IMPLICIT DOUBLE PRECISION (A-G,O-Z)
R11=(X21-X11)*(Y12-Y11)-(X12-X11)*(Y21-Y11)
R12=(X22-X11)*(Y12-Y11)-(X12-X11)*(Y22-Y11)
R21=(X11-X21)*(Y22-Y21)-(X22-X21)*(Y11-Y21)
R22=(X12-X21)*(Y22-Y21)-(X22-X21)*(Y12-Y21)
IF(R11*R12.LT.0.d0.AND.R21*R22.LT.0.d0) THEN
INTERSEC=1
R121=DSQRT((X21-X12)**2+(Y21-Y12)**2)
R122=DSQRT((X22-X12)**2+(Y22-Y12)**2)
IF(R122.GE.R121) THEN
IENDOF2=2
ELSE
IENDOF2=1
ENDIF
RETURN
ELSE
INTERSEC=0
IENDOF2=100
ENDIF
*=========================================================
return
end
******************************************************************
Subroutine IsoOCrkAdv(
1 Mode1or2,SIF1,SIF2,OEle,SIF1Eff,SIF2Eff,OCrkAdv)
Implicit Real*8 (A-H,O-Z)
Implicit Integer (I-N)
!---------maximum principal stress criterion---------------
If(Mode1or2.eq.1) then
If(DABS(SIF2).lt.1.d-10) then
OCrkAdv=OEle
else
OCrkAdv=OEle
1 +2.d0*DATAN((DSQRT(SIF1**2+8.d0*SIF2**2)-SIF1)/SIF2/4.d0)
endif
SIF1Eff=SIF1*DCOS((OCrkAdv-OEle)/2.d0)**3+3.d0*SIF2
1 *DCOS((OCrkAdv-OEle)/2.d0)**2*DSIN((OCrkAdv-OEle)/2.d0)
SIF2Eff=0.d0
!---------maximum shear stress criterion---------------
elseif(Mode1or2.eq.2) then
Call MSS_KField(SIF1,SIF2,SIF1Eff,SIF2Eff,AngleK2,LSuccess)
OCrkAdv=OEle+AngleK2
endif
*==============================================================
return
end
C******************SUBROUTINE MSS_KFIELD*****************
SUBROUTINE MSS_KFIELD(SIF1,SIF2,SIF1EFF,SIF2EFF,ANGLEK2,LSUCCESS)
Implicit Real*8 (A-H,O-Y)
Implicit Complex*16 (Z)
Implicit Integer (I-N)
DIMENSION A(4),ZX(3),XR(3),XI(3)
PI=4.d0*DATAN(1.d0)
TINY=1.0d-10
!==========================================================
IF(SIF1.LE.TINY.AND.DABS(SIF2).LE.TINY) THEN
SIF2EFF=SIF2
SIF1EFF=SIF1
ANGLEK2=0.d0
RETURN
ELSEIF(SIF1.LE.TINY) THEN
SIF2EFF=SIF2
SIF1EFF=SIF1
ANGLEK2=0.d0
RETURN
ELSEIF(DABS(SIF2).LE.TINY) THEN
AHALF=DATAN(DSQRT(2.d0)/2.d0)
IF(SIF2.LT.0.d0) AHALF=-DATAN(DSQRT(2.d0)/2.d0)
SIF2EFF=SIF1*(DSIN(AHALF)+DSIN(3.d0*AHALF))/4.d0+
1 SIF2*(DCOS(AHALF)+3.d0*DCOS(3.d0*AHALF))/4.d0
SIF1EFF=SIF1/4.d0*(3.d0*DCOS(AHALF)+DCOS(3.d0*AHALF))
1 -SIF2*0.75d0*(DSIN(AHALF)+DSIN(3.d0*AHALF))
ANGLEK2=AHALF*2.d0
RETURN
ENDIF
!===========================================================
A(1)=SIF2*2.d0
A(2)=-SIF1*2.d0
A(3)=-SIF2*7.d0
A(4)=SIF1
CALL DZPLRC(3,A,ZX) !SRRT(A,XR,XI,N,M,L,B)
XR=ZX
XI=ZX*(0.d0,1.d0)
!===========================================================
LSUCCESS=0
DO K=1,3
IF(XI(K)+1.0.EQ.1.0) THEN
AHALF=DATAN(XR(K))
DERIV2=-SIF1/16.d0*(DSIN(AHALF)+9.d0*DSIN(3.d0*AHALF))
1 -SIF2/16.d0*(DCOS(AHALF)+27.d0*DCOS(3.d0*AHALF))
IF(DERIV2.LT.0.d0) THEN
SIF2EFF=SIF1*(DSIN(AHALF)+DSIN(3.d0*AHALF))/4.d0+
1 SIF2*(DCOS(AHALF)+3.0*DCOS(3.d0*AHALF))/4.d0
SIF1EFF=SIF1/4.d0*(3.d0*DCOS(AHALF)+DCOS(3.d0*AHALF))
1 -SIF2*0.75d0*(DSIN(AHALF)+DSIN(3.d0*AHALF))
ANGLEK2=AHALF*2.d0
LSUCCESS=LSUCCESS+1
ENDIF
ENDIF
ENDDO
*=========================================================
RETURN
END
*********************************************************
Subroutine IsoCrkTipSifT(
1 NETip,SIF1,SIF2,TSs,
1 NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
1 JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
1 UNod,PNod,UHomoNod,PHomoNod)
Implicit Real*8 (A-H,O-Z)
Implicit Integer (I-N)
Dimension XEle(2,2,NECap),NodEle(NECap),OEle(NECap)
Dimension NthBdEle(NECap),NthGrpEle(NECap)
Dimension JCrkN0C1Tip2Ele(NECap),NthSlvMstEle(NECap)
Dimension XNod(2,3,NECap),ShpTNod(0:3,3,NECap)
Dimension UNod(2,3,NECap),PNod(2,3,NECap)
Dimension UHomoNod(2,3,NECap),PHomoNod(2,3,NECap)
!local
Dimension UPnt(2),SsPnt(2,2)
*================================================================
PI=4.d0*DATAN(1.d0)
*================================================================
NthBdTip=NthBdEle(NETip)
Call UserElasConst(ES0,PR0,NthBdTip)
*================================================================
!SIFs
DSE=DSQRT((XEle(1,2,NETip)-XEle(1,1,NETip))**2
1 +(XEle(2,2,NETip)-XEle(2,1,NETip))**2)
Do 30 IC=1,2
U1=UNod(IC,3,NETip)
U2=UNod(IC,2,NETip)
U3=UNod(IC,1,NETip)
CALL WilliamsExp(DSE,U1,U2,U3,AE0,AE1,AE2,AESIF)
If(IC.eq.1) SIF1=DSQRT(2.d0*PI)*ES0/(1.d0-PR0)/4.d0*AESIF
30 If(IC.eq.2) SIF2=DSQRT(2.d0*PI)*ES0/(1.d0-PR0)/4.d0*AESIF
*==============================================================
!T-stress
Ratio=1.001d0
X1Pnt=XEle(1,1,NETip)+Ratio*(XEle(1,2,NETip)-XEle(1,1,NETip))
X2Pnt=XEle(2,1,NETip)+Ratio*(XEle(2,2,NETip)-XEle(2,1,NETip))
JAcqSsY1N0=1
Call InterDisplStress(
1 X1Pnt,X2Pnt,NthBdTip,UPnt,JAcqSsY1N0,SsPnt,
1 NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
1 JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
1 UNod,PNod,UHomoNod,PHomoNod)
TSs=SsPnt(1,1)-SsPnt(2,2)
*==============================================================
return
end
*********************WILLIAMSEXP********************
SUBROUTINE WilliamsExp(RR,U1,U2,U3,AE0,AE1,AE2,AESIF)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION A(3,3),X(3),B(3)
DIMENSION R(3)
R(1)=RR/6.d0
R(2)=RR/2.d0
R(3)=RR/6.d0*5.d0
DO I=1,3
A(I,1)=R(I)**0.5d0
A(I,2)=R(I)**1.5d0
A(I,3)=R(I)**2.5d0
ENDDO
B(1)=U1
B(2)=U2
B(3)=U3
N=3
CALL DLSLRG(N,A,N,B,1,X)
AE0=X(1)
AE1=X(2)
AE2=X(3)
AESIF=X(1)
RETURN
END
*********************************************************
*********************************************************
*********************************************************
*********************************************************
Subroutine EPostProcess(NLoad,
1 NCmp,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
1 JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
1 NtotCfKCap,NthCfk,CfUP,
1 JBCU1P2,UNod,PNod,UHomoNod,PHomoNod)
Implicit Real*8 (A-H,O-Z)
Implicit Integer (I-N)
!master
Dimension XEle(2,2,NECap),NodEle(NECap),OEle(NECap)
Dimension NthBdEle(NECap),NthGrpEle(NECap)
Dimension JCrkN0C1Tip2Ele(NECap),NthSlvMstEle(NECap)
Dimension XNod(2,3,NECap),ShpTNod(0:3,3,NECap)
Dimension NthCfK(NECap,3,NECap,3,0:1)
Dimension CfUP(NCmp,NCmp*2,0:NtotCfKCap)
Dimension JBCU1P2(NCmp,3,NECap)
Dimension UNod(NCmp,3,NECap),PNod(NCmp,3,NECap)
Dimension UHomoNod(NCmp,3,NECap),PHomoNod(NCmp,3,NECap)
!local
Allocatable UNd(:,:,:),UCrkPair1(:,:,:),UCrkPair2(:,:,:)
Allocate(UNd(NCmp,3,NECap))
Allocate(UCrkPair1(NCmp,3,NECap),UCrkPair2(NCmp,3,NECap))
*********************************************************
PI=4.d0*DATAN(1.d0)
*********************************************************
!crack faces
Do 6400 NE=1,NtotEle
NthBdNE=NthBdEle(NE)
Do 6400 NN=1,NodEle(NE)
If(JCrkN0C1Tip2Ele(NE).gt.0) then
Do 6300 IC=1,NCmp
UNd(IC,NN,NE)=UHomoNod(IC,NN,NE)
Do 6200 ME=1,NtotEle
If(NthBdEle(ME).ne.NthBdNE) goto 6200
Do 6100 MN=1,NodEle(ME)
NthK=NthCfK(NE,NN,ME,MN,0)
Do 6100 JC=1,NCmp
If(JCrkN0C1Tip2Ele(ME).eq.0) then
CfUSF=CfUP(IC,JC,NthK)
CfPSF=CfUP(IC,JC+NCmp,NthK)
UNd(IC,NN,NE)=UNd(IC,NN,NE)
1 +CfUSF*(PNod(JC,MN,ME)-PHomoNod(JC,MN,ME))
1 -CfPSF*(UNod(JC,MN,ME)-UHomoNod(JC,MN,ME))
else
CfPSF=CfUP(IC,JC+NCmp,NthK)
UNd(IC,NN,NE)=UNd(IC,NN,NE)+CfPSF*UNod(JC,MN,ME)
endif
6100 continue
UCrkPair1(IC,NN,NE)=(UNd(IC,NN,NE)*2+UNod(IC,NN,NE))/2
UCrkPair2(IC,NN,NE)=(UNd(IC,NN,NE)*2-UNod(IC,NN,NE))/2
6200 continue
6300 continue
endif
6400 continue
*********************************************************
!Save!Save!Save!Save!Save!Save!Save!Save!Save!Save!Save!Save!Save
*================================================================
!deformed
Open(51,file='ZDefmd.txt',status='unknown')
Do 110 NE=1,NtotEle
VNorm1=DCOS(OEle(NE)-Pi/2.d0)
VNorm2=DSIN(OEle(NE)-Pi/2.d0)
Do 110 NN=⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -