📄 be2die.for
字号:
!============================================================
!crack!crack!crack!crack!crack!crack!crack!crack!crack!crack
elseif(JCrkN0C1Tip2Ele(NE).gt.0) then
!-----------------------------------------------------------
!resid!resid!resid!resid!resid!resid!resid!resid!resid!resid!resid
Do 1301 IC=1,NCmp
JCrkN0C1Resid=1
1301 Call ResidCrkN0C1(
1 JCrkN0C1Resid,IC,NN,NE,Resid(IC),ResidMax,
1 NCmp,NECap,NtotEle,NodEle,NthBdEle,JCrkN0C1Tip2Ele,
1 NtotCfKCap,NthCfK,CfUP,
1 UNod,PNod,UHomoNod,PHomoNod)
!--------------------------------------------------------------
!solving!solving!solving!solving!solving!solving!solving
If(NCmp.eq.2) then
WDmg=DSQRT(UNod(1,NN,NE)**2+UNod(2,NN,NE)**2+1.d-14)
else
WDmg=DSQRT(UNod(1,NN,NE)**2+UNod(2,NN,NE)**2+1.d-14
1 +UNod(3,NN,NE)**2)
endif
If(WDmg.lt.WDmgNod(NN,NE)) WDmg=WDmgNod(NN,NE)
NthGrp=NthGrpEle(NE)
EK=UserEKLineSpring(NthGrp,WDmg)
Call UserFrictParam(NthGrp,FC,WTaoD)
NCallNN=-1
1321 NCallNN=NCallNN+1
If(PNod(1,NN,NE).lt.0.d0) then !contact
IC=1
UNod(IC,NN,NE)=0.d0
CfPDiag=CfUP(IC,IC+NCmp,NthCf1Diag)
PNod(IC,NN,NE)=Resid(IC)+PNod1(IC,NN,NE)-CfPDiag*UNod(IC,NN,NE)
If(PNod(1,NN,NE).gt.0.d0.and.NCallNN.le.3) goto 1321
If(NCmp.eq.2) then
WTaoDTp=DABS(UNod(2,NN,NE)-UNod0(2,NN,NE))
else
WTaoDTp=DSQRT((UNod(2,NN,NE)-UNod0(2,NN,NE))**2
1 +(UNod(3,NN,NE)-UNod0(3,NN,NE))**2+1.d-14)
endif
IF(WTaoDTp.lt.WTaoD) WTaoDTp=WTaoD
Do IC=2,NCmp
CfPDiag=CfUP(IC,IC+NCmp,NthCf1Diag)
UNod(IC,NN,NE)=UNod(IC,NN,NE)+
1 AlphaC*Resid(IC)/(EK-FC*PNod(1,NN,NE)/WTaoDTp-CfPDiag)
enddo
If(NCmp.eq.2) then
WTaoDTp=DABS(UNod(2,NN,NE)-UNod0(2,NN,NE))
else
WTaoDTp=DSQRT((UNod(2,NN,NE)-UNod0(2,NN,NE))**2
1 +(UNod(3,NN,NE)-UNod0(3,NN,NE))**2+1.d-14)
endif
IF(WTaoDTp.lt.WTaoD) WTaoDTp=WTaoD
Do IC=2,NCmp
PNod(IC,NN,NE)=EK*UNod(IC,NN,NE)-
1 FC*PNod(1,NN,NE)*(UNod(IC,NN,NE)-UNod0(IC,NN,NE))/WTaoDTp
enddo
else !opening
Do IC=1,NCmp
CfPDiag=CfUP(IC,IC+NCmp,NthCf1Diag)
UNod(IC,NN,NE)=UNod1(IC,NN,NE)+AlphaC*Resid(IC)/(EK-CfPDiag)
PNod(IC,NN,NE)=EK*UNod(IC,NN,NE)
enddo
If(PNod(1,NN,NE).lt.0.d0.and.NCallNN.le.3) goto 1321
endif
If(NCallNN.gt.NCallCntctMax) NCallCntctMax=NCallNN
!============================================================
!non crack!non crack!non crack!non crack!non crack!non crack!non crack
!MstSlv!MstSlv!MstSlv!MstSlv!MstSlv!MstSlv!MstSlv!MstSlv!MstSlv!MstSlv
elseif(NthSlvMstEle(NE).gt.0) then
NESlv=NthSlvMstEle(NE)
NNSlv=NodEle(NE)+1-NN
NthCf0DiagSlv=NthCfK(NESlv,NNSlv,NESlv,NNSlv,0)
!-----------------------------------------------------------
!resid!resid!resid!resid!resid!resid!resid!resid!resid!resid!resid
Do 1401 IC=1,NCmp
JCrkN0C1Resid=0
Call ResidCrkN0C1(
1 JCrkN0C1Resid,IC,NN,NE,Resid(IC),ResidMax,
1 NCmp,NECap,NtotEle,NodEle,NthBdEle,JCrkN0C1Tip2Ele,
1 NtotCfKCap,NthCfK,CfUP,
1 UNod,PNod,UHomoNod,PHomoNod)
1401 Call ResidCrkN0C1(
1 JCrkN0C1Resid,IC,NNSlv,NESlv,ResidSlv(IC),ResidMax,
1 NCmp,NECap,NtotEle,NodEle,NthBdEle,JCrkN0C1Tip2Ele,
1 NtotCfKCap,NthCfK,CfUP,
1 UNod,PNod,UHomoNod,PHomoNod)
!----------------------------------------------------------------
!solving!solving!solving!solving!solving!solving!solving!solving
Theta=OEle(NESlv)-OEle(NE)-PI
If(NCmp.eq.2) then
WDmg=DSQRT(
1 (-UNod(1,NNSlv,NESlv)-UNod(1,NN,NE)+DNormGapMstNod(NN,NE))**2
1 +(-UNod(2,NNSlv,NESlv)-UNod(2,NN,NE))**2 +1.d-14)
else
WDmg=DSQRT(
1 (-UNod(1,NNSlv,NESlv)-UNod(1,NN,NE)+DNormGapMstNod(NN,NE))**2
1 +(-UNod(2,NNSlv,NESlv)-UNod(2,NN,NE))**2
1 +(UNod(3,NNSlv,NESlv)-UNod(3,NN,NE))**2+1.d-14)
endif
If(WDmg.lt.WDmgNod(NN,NE)) WDmg=WDmgNod(NN,NE)
NthGrp=NthGrpEle(NE)
EK=UserEKLineSpring(NthGrp,WDmg)
Call UserFrictParam(NthGrp,FC,WTaoD)
NCallNN=-1
1411 NCallNN=NCallNN+1
If(PNod(1,NN,NE).lt.0.d0) then !contact
IC=1
CfUDiag=CfUP(IC,IC,NthCf0Diag)
CfUDiagSlv=CfUP(IC,IC,NthCf0DiagSlv)
PNod(IC,NN,NE)=PNod1(IC,NN,NE)-Alpha*
1 (Resid(IC)+ResidSlv(IC))/(CfUDiag+CfUDiagSlv)
If(PNod(1,NN,NE).ge.0.d0.and.NCallNN.le.3) goto 1411
U1NodDif=-DNormGapMstNod(NN,NE) !U1Dif=-(UnSlv*COS(qSlv-qMst-pi)-UtSlv*SIN)-UnMst
CfPDiag=CfUP(IC,IC+NCmp,NthCf0Diag)
CfPDiagSlv=CfUP(IC,IC+NCmp,NthCf0DiagSlv)
U1NodSum=-UNod1(IC,NNSlv,NESlv)+UNod1(IC,NN,NE)+ !USum=-USlave+UMaster
1 Alpha*(-ResidSlv(IC)+Resid(IC))/(1.d0+CfPDiag+CfPDiagSlv)
!
If(NCmp.eq.2) then
U2Frct0=-UNod0(2,NN,NE)-(UNod0(2,NNSlv,NESlv)*DCOS(Theta)+
1 UNod0(1,NNSlv,NESlv)*DSIN(Theta))
U2Frct1=-UNod1(2,NN,NE)-(UNod1(2,NNSlv,NESlv)*DCOS(Theta)+
1 UNod1(1,NNSlv,NESlv)*DSIN(Theta))
WTaoDTp=DABS(U2Frct1-U2Frct0)
else
U2Frct0=-UNod0(2,NN,NE)-(UNod0(2,NNSlv,NESlv)*DCOS(Theta)+
1 UNod0(1,NNSlv,NESlv)*DSIN(Theta))
U2Frct1=-UNod1(2,NN,NE)-(UNod1(2,NNSlv,NESlv)*DCOS(Theta)+
1 UNod1(1,NNSlv,NESlv)*DSIN(Theta))
U3Frct0=UNod0(3,NNSlv,NESlv)-UNod0(3,NN,NE)
U3Frct1=UNod1(3,NNSlv,NESlv)-UNod1(3,NN,NE)
WTaoDTp=DSQRT((U2Frct1-U2Frct0)**2+(U3Frct1-U3Frct0)**2)
endif
If(WTaoDTp.lt.WTaoD) WTaoDTp=WTaoD
!
IC=2
CfUDiag=CfUP(IC,IC,NthCf0Diag)
CfPDiag=CfUP(IC,IC+NCmp,NthCf0Diag)
CfUDiagSlv=CfUP(IC,IC,NthCf0DiagSlv)
CfPDiagSlv=CfUP(IC,IC+NCmp,NthCf0DiagSlv)
U2NodDif=U2Frct1+ !U2Dif=-(UnSlv*SIN(qSlv-qMst-pi)+UtSlv*COS)-U2Mst
1 Alpha*(-(ResidSlv(1)*DSIN(Theta)+ResidSlv(2)*DCOS(Theta))
1 -Resid(IC))/(1.d0+CfPDiag+CfPDiagSlv+
1 (EK-FC*PNod(1,NN,NE)/WTaoDTp)*(CfUDiag+CfUDiagSlv))
U2NodSum=-UNod1(IC,NNSlv,NESlv)+UNod1(IC,NN,NE)+ !USum=-USlave+UMaster
1 Alpha*(-ResidSlv(IC)+Resid(IC))/(1.d0+CfPDiagSlv+CfPDiag)
UNod(1,NNSlv,NESlv)=((U1NodSum+U1NodDif)*(1.d0+DCOS(Theta))+
1 (U2NodSum+U2NodDif)*DSIN(Theta))/(-2.d0-2.d0*DCOS(Theta))
UNod(2,NNSlv,NESlv)=(U2NodSum+U2NodDif+
1 UNod(1,NNSlv,NESlv)*DSIN(Theta))/(-1.d0-DCOS(Theta))
UNod(1,NN,NE)=U1NodSum+UNod(1,NNSlv,NESlv)
UNod(2,NN,NE)=U2NodSum+UNod(2,NNSlv,NESlv)
If(NCmp.eq.3) then
IC=3
CfUDiag=CfUP(IC,IC,NthCf0Diag)
CfPDiag=CfUP(IC,IC+NCmp,NthCf0Diag)
CfUDiagSlv=CfUP(IC,IC,NthCf0DiagSlv)
CfPDiagSlv=CfUP(IC,IC+NCmp,NthCf0DiagSlv)
UNodDif=UNod1(IC,NNSlv,NESlv)-UNod1(IC,NN,NE)+ !UDif=USlave-UMaster
1 Alpha*(ResidSlv(IC)-Resid(IC))/(1.d0+CfPDiag+CfPDiagSlv+
1 (EK-FC*PNod(1,NN,NE)/WTaoDTp)*(CfUDiag+CfUDiagSlv))
UNodSum=UNod1(IC,NNSlv,NESlv)+UNod1(IC,NN,NE)+ !USum=USlave+UMaster
1 Alpha*(ResidSlv(IC)+Resid(IC))/(1.d0+CfPDiagSlv+CfPDiag)
UNod(IC,NN,NE)=(UNodSum-UNodDif)/2.d0
UNod(IC,NNSlv,NESlv)=(UNodSum+UNodDif)/2.d0
endif
If(NCmp.eq.2) then
U2Frct=-UNod(2,NN,NE)-(UNod(2,NNSlv,NESlv)*DCOS(Theta)+
1 UNod(1,NNSlv,NESlv)*DSIN(Theta))
else
U2Frct=-UNod(2,NN,NE)-(UNod(2,NNSlv,NESlv)*DCOS(Theta)+
1 UNod(1,NNSlv,NESlv)*DSIN(Theta))
U3Frct=UNod(3,NNSlv,NESlv)-UNod(3,NN,NE)
endif
IC=2
PNod(2,NN,NE)=EK*U2Frct
1 -FC*PNod(1,NN,NE)*((U2Frct-U2Frct0)/WTaoDTp)
PNod(1,NNSlv,NESlv)=PNod(1,NN,NE)*DCOS(Theta)+
1 PNod(2,NN,NE)*DSIN(Theta)
PNod(2,NNSlv,NESlv)=-PNod(1,NN,NE)*DSIN(Theta)+
1 PNod(2,NN,NE)*DCOS(Theta)
If(NCmp.eq.3) then
PNod(3,NN,NE)=EK*U3Frct
1 -FC*PNod(1,NN,NE)*((U3Frct-U3Frct0)/WTaoDTp)
PNod(3,NNSlv,NESlv)=-PNod(3,NN,NE)
endif
else !opening
Do 1415 IC=1,2
CfUDiag=CfUP(IC,IC,NthCf0Diag)
CfPDiag=CfUP(IC,IC+NCmp,NthCf0Diag)
CfUDiagSlv=CfUP(IC,IC,NthCf0DiagSlv)
CfPDiagSlv=CfUP(IC,IC+NCmp,NthCf0DiagSlv)
UNodDif=-UNod1(IC,NNSlv,NESlv)-UNod1(IC,NN,NE)+ !UDif=-USlave-UMaster
1 Alpha*(-ResidSlv(IC)-Resid(IC))/
1 (1.d0+CfPDiagSlv+CfPDiag+EK*(CfUDiagSlv+CfUDiag))
UNodSum=-UNod1(IC,NNSlv,NESlv)+UNod1(IC,NN,NE)+ !USum=-USlave+UMaster
1 Alpha*(-ResidSlv(IC)+Resid(IC))/(1.d0+CfPDiagSlv+CfPDiag)
UNod(IC,NN,NE)=(UNodSum-UNodDif)/2.d0
1415 UNod(IC,NNSlv,NESlv)=(-UNodSum-UNodDif)/2.d0
U1NodDif=-UNod(1,NN,NE) !U1Dif=-(UnSlv*COS(qSlv-qMst-pi)-UtSlv*SIN)-UnMst
1 -(UNod(1,NNSlv,NESlv)*DCOS(Theta)-
1 UNod(2,NNSlv,NESlv)*DSIN(Theta))
U2NodDif=-UNod(2,NN,NE) !U2Dif=-(UnSlv*SIN(qSlv-qMst-pi)+UtSlv*COS)-U2Mst
1 -(UNod(1,NNSlv,NESlv)*DSIN(Theta)+
1 UNod(2,NNSlv,NESlv)*DCOS(Theta))
PNod(1,NN,NE)=EK*(U1NodDif+DNormGapMstNod(NN,NE))
PNod(2,NN,NE)=EK*U2NodDif
PNod(1,NNSlv,NESlv)=PNod(1,NN,NE)*DCOS(Theta)+
1 PNod(2,NN,NE)*DSIN(Theta)
PNod(2,NNSlv,NESlv)=-PNod(1,NN,NE)*DSIN(Theta)+
1 PNod(2,NN,NE)*DCOS(Theta)
If(PNod(1,NN,NE).lt.0.d0.and.NCallNN.le.3) goto 1411
If(NCmp.eq.3) then
IC=3
CfUDiag=CfUP(IC,IC,NthCf0Diag)
CfPDiag=CfUP(IC,IC+NCmp,NthCf0Diag)
CfUDiagSlv=CfUP(IC,IC,NthCf0DiagSlv)
CfPDiagSlv=CfUP(IC,IC+NCmp,NthCf0DiagSlv)
UNodDif=UNod1(IC,NNSlv,NESlv)-UNod1(IC,NN,NE)+ !UDif=USlave-UMaster
1 Alpha*(ResidSlv(IC)-Resid(IC))/
1 (1.d0+CfPDiagSlv+CfPDiag+EK*(CfUDiagSlv+CfUDiag))
UNodSum=UNod1(IC,NNSlv,NESlv)+UNod1(IC,NN,NE)+ !USum=USlave+UMaster
1 Alpha*(ResidSlv(IC)+Resid(IC))/(1.d0+CfPDiagSlv+CfPDiag)
UNod(IC,NN,NE)=(UNodSum-UNodDif)/2.d0
UNod(IC,NNSlv,NESlv)=(UNodSum+UNodDif)/2.d0
PNod(IC,NN,NE)=EK*(UNod(IC,NNSlv,NESlv)-UNod(IC,NN,NE))
PNod(IC,NNSlv,NESlv)=-PNod(IC,NN,NE)
endif
endif
If(NCallNN.gt.NCallCntctMax) NCallCntctMax=NCallNN
!=================================================================
endif
!=================================================================
1980 continue
*===================================================================
!error control!error control!error control!error control!error control
UErMax=0.d0
PErMax=0.d0
Do 2010 NE=1,NtotEle
Do 2010 NN=1,NodEle(NE)
Do 2010 I=1,NCmp
UEr=DABS(UNod(I,NN,NE)-UNod1(I,NN,NE))
If(UEr.gt.UErMax) UErMax=UEr
PEr=DABS(PNod(I,NN,NE)-PNod1(I,NN,NE))
2010 If(PEr.gt.PErMax) PErMax=PEr
If((Iterat-1)/IteratShow*IteratShow.eq.Iterat-1) then
write(*,3110) Iterat,UErMax,PErMax,NCallCntctMax
endif
3110 format(1x,I4,2(1x,F13.9),1x,I3)
If(UErMax.gt.UErrAllow.or.PErMax.gt.PErrAllow) goto 1101
!
Write(*,3120) Iterat,UErMax,PErMax,NCallCntctMax
3120 format(1x,'Finish: Iterat=',1x,I4,2(1x,F13.9),1x,I3)
write(*,*)
*********************************************************
Call UpdateDamage(
1 NCmp,NECap,NtotEle,NodEle,NthSlvMstEle,UNod,WDmgNod)
*********************************************************
Call 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)
*********************************************************
Open(51,file='DUPSolut.txt',status='unknown')
Do NE=1,NtotEle
Do NN=1,NodEle(NE)
write(51,1110)
1 (UNod(IC,NN,NE),IC=1,NCmp),(PNod(IC,NN,NE),IC=1,NCmp)
enddo
enddo
close(51)
1110 format(20(1x,f22.9))
*********************************************************
!crack tip fields
NtotEleTp=NtotEle
Do NE=1,NtotEleTp
If(JCrkN0C1Tip2Ele(NE).eq.2) then
JAdvY1N0=0
Call UserCrkTip(
1 NE,JAdvY1N0,OCrkAdv,ACrkAdv,
1 NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
1 JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
1 UNod,PNod,UHomoNod,PHomoNod)
If(JAdvY1N0.eq.1.and.NtotEle.lt.NECap)
1 Call 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)
endif
enddo
!=====================================================
If(NtotEle.gt.NtotEleTp) goto 1001
*********************************************************
return
end
*********************************************************
*********************************************************
*********************************************************
Subroutine ResidCrkN0C1(
1 JCrkN0C1,IC,NN,NE,Resid,ResidMax,
1 NCmp,NECap,NtotEle,NodEle,NthBdEle,JCrkN0C1Tip2Ele,
1 NtotCfKCap,NthCfK,CfUP,
1 UNod,PNod,UHomoNod,PHomoNod)
Implicit Real*8 (A-H,O-Z)
Implicit Integer (I-N)
Dimension NodEle(NECap),NthBdEle(NECap)
Dimension JCrkN0C1Tip2Ele(NECap)
Dimension NthCfK(NECap,3,NECap,3,0:1)
Dimension CfUP(NCmp,NCmp*2,0:NtotCfKCap)
Dimension UNod(NCmp,3,NECap),PNod(NCmp,3,NECap)
Dimension UHomoNod(NCmp,3,NECap),PHomoNod(NCmp,3,NECap)
*========================================================
NthBdNE=NthBdEle(NE)
If(JCrkN0C1.eq.0) Resid=-(UNod(IC,NN,NE)-UHomoNod(IC,NN,NE))/2.d0
If(JCrkN0C1.eq.1) Resid=-(PNod(IC,NN,NE)-PHomoNod(IC,NN,NE))
Do 1211 ME=1,NtotEle
If(NthBdEle(ME).ne.NthBdNE) goto 1211
Do 1210 MN=1,NodEle(ME)
NthK=NthCfK(NE,NN,ME,MN,JCrkN0C1)
Do 1210 JC=1,NCmp
If(JCrkN0C1Tip2Ele(ME).eq.0) then
CfUSF=CfUP(IC,JC,NthK)
CfPSF=CfUP(IC,JC+NCmp,NthK)
Resid=Resid+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)
Resid=Resid+CfPSF*UNod(JC,MN,ME)
endif
1210 continue
1211 continue
If(DABS(Resid).gt.ResidMax) ResidMax=DABS(Resid)
*=========================================================
return
end
*************************************************************
*************************************************************
*************************************************************
Subroutine UpdateDamage(
1 NCmp,NECap,NtotEle,NodEle,NthSlvMstEle,UNod,WDmgNod)
Implicit Real*8 (A-H,O-Z)
Implicit Integer (I-N)
Dimension NodEle(NECap),NthSlvMstEle(NECap)
Dimension UNod(NCmp,3,NECap),WDmgNod(3,NECap)
*=========================================================
Do 100 NE=1,NtotEle
Do 100 NN=1,NodEle(NE)
If(NthSlvMstEle(NE).eq.0) then
If(NCmp.eq.2) then
WD=DSQRT(UNod(1,NN,NE)**2+UNod(2,NN,NE)**2)
else
WD=DSQRT(UNod(1,NN,NE)**2+UNod(2,NN,NE)**2
1 +UNod(3,NN,NE)**2)
endif
elseif(NthSlvMstEle(NE).gt.0) then
NESlv=NthSlvMstEle(NE)
NNSlv=NodEle(NE)+1-NN
If(NCmp.eq.2) then
WD=DSQRT((-UNod(1,NNSlv,NESlv)-UNod(1,NN,NE))**2
1 +(-UNod(2,NNSlv,NESlv)-UNod(2,NN,NE))**2)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -