⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 be2die.for

📁 两维弹性边界元程序
💻 FOR
📖 第 1 页 / 共 5 页
字号:
🔍
12345 
	R1=DSQRT((XFE(1,2)-XFE(1,1))**2+(XFE(2,2)-XFE(2,1))**2)
	R2=R1*2
	Do 110 I=1,NCmp
	Do 110 J=1,NCmp
	CLogGnU(I,J)=(GnUP2(I,J)-GnUP1(I,J))/(DLOG(R2)-DLOG(R1))
	CPlnGnU(I,J)=GnUP2(I,J)-CLogGnU(I,J)*DLOG(R2)
	CPln1=GnUP1(I,J)-CLogGnU(I,J)*DLOG(R1)
	CPln2=GnUP2(I,J)-CLogGnU(I,J)*DLOG(R2)
	CR_1GnP(I,J)=GnPP2(I,J)*R2
	Do 110 K=1,NCmp
	CR_1GnUD(I,J,K)=GnUDP2(I,J,K)*R2
110	CR_2GnPD(I,J,K)=GnPDP2(I,J,K)*R2*R2
*==============================================================
	R=DSQRT((XFE(1,2)-XFE(1,1))**2+(XFE(2,2)-XFE(2,1))**2)
	R1S=DSQRT((XSP(1)-XFE(1,1))**2+(XSP(2)-XFE(2,1))**2)
	R2S=DSQRT((XSP(1)-XFE(1,2))**2+(XSP(2)-XFE(2,2))**2)
	If(R.gt.R1S.and.R.gt.R2S) then
	 S1S=-R1S/R
	 S2S=R2S/R
	 TSP=-S1S
	elseif(R2S.gt.R.and.R2S.gt.R1S) then
	 S1S=R1S/R
	 S2S=R2S/R
	 TSP=-S1S
	elseif(R1S.gt.R.and.R1S.gt.R2S) then
	 S1S=-R1S/R
	 S2S=-R2S/R
	 TSP=-S1S
	endif
!------------------------------------------------------------
	ShpSFE(0)=ShpTFE(0)+ShpTFE(1)*TSP+ShpTFE(2)*TSP*TSP+
	1	ShpTFE(3)*TSP*TSP*TSP
	ShpSFE(1)=ShpTFE(1)+2.d0*ShpTFE(2)*TSP+
	1	3.d0*ShpTFE(3)*TSP*TSP
	ShpSFE(2)=ShpTFE(2)+3.d0*ShpTFE(3)*TSP
	ShpSFE(3)=ShpTFE(3)
!============================================================
	GIntPln=R*(ShpSFE(0)*(S2S-S1S)+
	1	ShpSFE(1)*(S2S**2/2.d0-S1S**2/2.d0)+
	2	ShpSFE(2)*(S2S**3/3.d0-S1S**3/3.d0)+
	3	ShpSFE(3)*(S2S**4/4.d0-S1S**4/4.d0))
!---------------------------------------------------------------
	DLOGS1S=DLOG(DABS(S1S))
	DLOGS2S=DLOG(DABS(S2S))
	GIntLog=R*(
	1	ShpSFE(0)*((S2S*DLOGS2S-S2S)-
     1				 (S1S*DLOGS1S-S1S))+
	1	ShpSFE(1)*((S2S**2/2.d0*DLOGS2S-S2S**2/4.d0)-
	1				 (S1S**2/2.d0*DLOGS1S-S1S**2/4.d0))+
	2	ShpSFE(2)*((S2S**3/3.d0*DLOGS2S-S2S**3/9.d0)-
	2				 (S1S**3/3.d0*DLOGS1S-S1S**3/9.d0))+
	3	ShpSFE(3)*((S2S**4/4.d0*DLOGS2S-S2S**4/16.d0)-
	3				 (S1S**4/4.d0*DLOGS1S-S1S**4/16.d0)))+
     4	GIntPln*DLOG(R)
!---------------------------------------------------------------
	GIntR_1=ShpSFE(0)*DLOG(DABS(S2S/S1S))+
	1	ShpSFE(1)*(S2S-S1S)+
	2	ShpSFE(2)*(S2S**2/2.d0-S1S**2/2.d0)+
	3	ShpSFE(3)*(S2S**3/3.d0-S1S**3/3.d0)
!---------------------------------------------------------------
	GIntR_2=1.d0/R*(
	1	ShpSFE(0)*(-1.d0/S2S+1.d0/S1S)+
	1	ShpSFE(1)*DLOG(DABS(S2S/S1S))+
	2	ShpSFE(2)*(S2S-S1S)+
	3	ShpSFE(3)*(S2S**2/2.d0-S1S**2/2.d0))
*==========================================================
*==========================================================
*==========================================================
!CfU
	Do 30 I=1,NCmp
	Do 30 J=1,NCmp
30	Cf(I,J,0)=CPlnGnU(I,J)*GIntPln+CLogGnU(I,J)*GIntLog
!=========================================================
!CfP
	Do 60 I=1,NCmp
	Do 60 J=1,NCmp
60	Cf(I,J+NCmp,0)=CR_1GnP(I,J)*GIntR_1	
*==========================================================
	If(JDeriv.eq.0) goto 9999
*========================================================
!CfUD
	Do 130 J=1,NCmp
	Do 130 K=1,NCmp
	Cf(J,K,1)=0.d0
	Do 130 I=1,NCmp
130	Cf(J,K,1)=Cf(J,K,1)+CR_1GnUD(J,I,K)*GIntR_1*VNormSP(I)
!=========================================================
!CfPD
	Do 160 J=1,NCmp
	Do 160 K=1,NCmp
	Cf(J,K+NCmp,1)=0.d0
	Do 160 I=1,NCmp
160	Cf(J,K+NCmp,1)=Cf(J,K+NCmp,1)+CR_2GnPD(J,I,K)*GIntR_2*VNormSP(I)
*==============================================================
9999	DeAllocate(ShpSFE,VNormSP,VNormFE)
	DeAllocate(GnUP1,GnPP1,GnUDP1,GnPDP1)
	DeAllocate(GnUP2,GnPP2,GnUDP2,GnPDP2)
	DeAllocate(CPlnGnU,CLogGnU,CR_1GnP,CR_1GnUD,CR_2GnPD)
	return
	end
*************************************************************
	Subroutine TransCftoLocal(NCmp,CfOld,VNormS,VNormF,CfNew)
      Implicit Real*8 (A-G,O-Y)
	Implicit Real*4 (H)
	Implicit Integer (I-N)
	Implicit Complex*16 (Z)
	Dimension CfOld(NCmp,NCmp),CfNew(NCmp,NCmp)
	Dimension VNormS(NCmp),VNormF(NCmp)
	Allocatable TransfmS(:,:),TransfmF(:,:)
	Allocate(TransfmS(NCmp,NCmp),TransfmF(NCmp,NCmp))
!transform to local coordinates
 	TransfmS(1,1)=VNormS(1)
	TransfmS(1,2)=VNormS(2)
	TransfmS(2,1)=-VNormS(2)
	TransfmS(2,2)=VNormS(1)
	TransfmF(1,1)=VNormF(1)
	TransfmF(1,2)=VNormF(2)
	TransfmF(2,1)=-VNormF(2)
	TransfmF(2,2)=VNormF(1)
	If(NCmp.eq.3) then
	 TransfmS(1,3)=0.d0
	 TransfmS(2,3)=0.d0
	 TransfmS(3,1)=0.d0
	 TransfmS(3,2)=0.d0
	 TransfmS(3,3)=1.d0
	 TransfmF(1,3)=0.d0
	 TransfmF(2,3)=0.d0
	 TransfmF(3,1)=0.d0
	 TransfmF(3,2)=0.d0
	 TransfmF(3,3)=1.d0
	endif
	Do 520 I=1,NCmp
	Do 520 J=1,NCmp
	CfNew(I,J)=0.d0
	Do 520 K=1,NCmp
	Do 520 L=1,NCmp
520	CfNew(I,J)=CfNew(I,J)+TransfmS(I,K)*CfOld(K,L)*TransfmF(J,L)
	DeAllocate(TransfmS,TransfmF)
	return
	end
*********************************************************************
*********************************************************************
*********************************************************************
	Subroutine CBCond(NLoad,NLoadMax,
	1	NCmp,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	JBCU1P2,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 JBCU1P2(NCmp,3,NECap)
	Dimension UNod(NCmp,3,NECap),PNod(NCmp,3,NECap)
	Dimension UHomoNod(NCmp,3,NECap),PHomoNod(NCmp,3,NECap)
*================================================
	Call UserBCond(NLoad,NLoadMax,
	1	NECap,NtotEle,NodEle,OEle,NthBdEle,NthGrpEle,XNod,
	1	JBCU1P2,UNod,PNod)	!output
*===================================================
	Do NE=1,NtotEle
	 Call HomoFldsEle(NLoad,
     1	NCmp,NodEle(NE),OEle(NE),NthBdEle(NE),XNod(1,1,NE),
     2		UHomoNod(1,1,NE),PHomoNod(1,1,NE))	!output
	enddo
*================================================================
	return
	end
***************************************************
	Subroutine HomoFldsEle(NLoad,
     1	NCmp,NodEle,OEle,NthBdEle,XNod,UHomoNod,PHomoNod)
	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Dimension XNod(2,3)
	Dimension UHomoNod(NCmp,3),PHomoNod(NCmp,3)
!
	Allocatable SsHomo(:,:),SnHomo(:,:)
	Allocate(SsHomo(NCmp,NCmp),SnHomo(NCmp,NCmp))
*===================================================
	Pi=4.d0*DATAN(1.d0)
*====================================================
	NthBd=NthBdEle
	Call UserElasConst(ES0,PR0,NthBd)
	SnHomo=0.d0
	SsHomo=0.d0
!SsHomo due to EigenSn
	Const=ES0*(1.d0+PR0)/(0.5d0-PR0)*UserIsoEigenSn(NLoad,NthBd)
	Do 10 I=1,NCmp
	Do 10 J=1,NCmp
10	If(I.eq.J) SsHomo(I,J)=SsHomo(I,J)-Const
!SsHomo and SnHomo due to far field loading
	Do 20 I=1,NCmp
	Do 20 J=1,NCmp
20	SsHomo(I,J)=SsHomo(I,J)+UserFarSs(NLoad,I,J,NthBd)
	SsKK=(UserFarSs(NLoad,1,1,NthBd)
	1	+UserFarSs(NLoad,2,2,NthBd))*(1.d0+PR0)
	Do 30 I=1,NCmp
	Do 30 J=1,NCmp
	SnHomo(I,J)=UserFarSs(NLoad,I,J,NthBd)/2.d0/ES0
30	If(I.eq.J) SnHomo(I,J)=SnHomo(I,J)-PR0/(1.d0+PR0)*SsKK/2.d0/ES0
*===================================================
	VNorm1=DCOS(OEle-Pi/2.d0)
	VNorm2=DSIN(OEle-Pi/2.d0)
	Do 60 NN=1,NodEle
	UHomo1=SnHomo(1,1)*XNod(1,NN)+SnHomo(1,2)*XNod(2,NN)/2
	UHomo2=SnHomo(2,1)*XNod(1,NN)/2+SnHomo(2,2)*XNod(2,NN)
	UHomoNod(1,NN)=UHomo1*VNorm1+UHomo2*VNorm2
	UHomoNod(2,NN)=-UHomo1*VNorm2+UHomo2*VNorm1
	If(NCmp.eq.3) UHomoNod(3,NN)=
	1		SnHomo(3,1)*XNod(1,NN)/2+SnHomo(3,2)*XNod(2,NN)/2
	PHomo1=SsHomo(1,1)*VNorm1+SsHomo(1,2)*VNorm2
	PHomo2=SsHomo(2,1)*VNorm1+SsHomo(2,2)*VNorm2
	PHomoNod(1,NN)=PHomo1*VNorm1+PHomo2*VNorm2
	PHomoNod(2,NN)=-PHomo1*VNorm2+PHomo2*VNorm1
60	If(NCmp.eq.3) PHomoNod(3,NN)=
	1		SsHomo(3,1)*VNorm1+SsHomo(3,2)*VNorm2
100	continue
*================================================================
	DeAllocate(SsHomo,SnHomo)
	return
	end
***************************************************
	Subroutine HomoFldsPnt(NLoad,NCmp,NthBd,XNod,UHomoPnt,SsHomoPnt)
	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Dimension XNod(2)
	Dimension UHomoPnt(NCmp),SsHomoPnt(NCmp,NCmp)
!
	Allocatable SsHomo(:,:),SnHomo(:,:)
	Allocate(SsHomo(NCmp,NCmp),SnHomo(NCmp,NCmp))
*===================================================
	Pi=4.d0*DATAN(1.d0)
*====================================================
	Call UserElasConst(ES0,PR0,NthBd)
	Do 5 I=1,NCmp
	Do 5 J=1,NCmp
	SnHomo(I,J)=0.d0
5	SsHomo(I,J)=0.d0
!SsHomo due to EigenSn
	Const=ES0*(1.d0+PR0)/(0.5d0-PR0)*UserIsoEigenSn(NLoad,NthBd)
	Do 10 I=1,NCmp
	Do 10 J=1,NCmp
10	If(I.eq.J) SsHomo(I,J)=SsHomo(I,J)-Const
!SsHomo and SnHomo due to far field loading
	Do 20 I=1,NCmp
	Do 20 J=1,NCmp
20	SsHomo(I,J)=SsHomo(I,J)+UserFarSs(NLoad,I,J,NthBd)
	SsKK=(UserFarSs(NLoad,1,1,NthBd)
	1	+UserFarSs(NLoad,2,2,NthBd))*(1.d0+PR0)
	Do 30 I=1,NCmp
	Do 30 J=1,NCmp
	SnHomo(I,J)=UserFarSs(NLoad,I,J,NthBd)/2.d0/ES0
30	If(I.eq.J) SnHomo(I,J)=SnHomo(I,J)-PR0/(1.d0+PR0)*SsKK/2.d0/ES0
*===================================================
	UHomoPnt(1)=SnHomo(1,1)*XNod(1)+SnHomo(1,2)*XNod(2)/2
	UHomoPnt(2)=SnHomo(2,1)*XNod(1)/2+SnHomo(2,2)*XNod(2)
	If(NCmp.eq.3) 
	1UHomoPnt(3)=SnHomo(3,1)*XNod(1)/2+SnHomo(3,2)*XNod(2)/2
	Do 60 I=1,NCmp
	Do 60 J=1,NCmp
60	SsHomoPnt(I,J)=SsHomo(I,J)
*================================================================
	DeAllocate(SsHomo,SnHomo)
	return
	end
*********************************************************
*********************************************************
*********************************************************
	Subroutine DSolver(NLoad,
	1	NCmp,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,DNormGapMstNod,
	1	NtotCfKCap,NtotCfK,NthCfk,CfUP,
	1	JBCU1P2,UNod,PNod,UHomoNod,PHomoNod)
	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Allocatable WDmgNod(:,:),UNod0(:,:,:),PNod0(:,:,:)
	Allocatable UNod1(:,:,:),PNod1(:,:,:)
	Allocatable Resid(:),ResidSlv(:)
*********************************************************************
	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 DNormGapMstNod(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)
*********************************************************************
	Allocate(WDmgNod(3,NECap))
	Allocate(UNod0(NCmp,3,NECap),PNod0(NCmp,3,NECap))
	Allocate(UNod1(NCmp,3,NECap),PNod1(NCmp,3,NECap))
	Allocate(Resid(NCmp),ResidSlv(NCmp))
*********************************************************************
	PI=4.d0*DATAN(1.d0)
*********************************************************************
	Call UserSolvParam(Alpha,AlphaC,UErrAllow,PErrAllow,IteratShow)
c	Alpha=0.9
c	AlphaC=1.2d0
c	UErrAllow=1.d-5
c	PErrAllow=1.d-5
c	IteratShow=50
*********************************************************************
!initial state
	If(NLoad.eq.1) then
	 Do NE=1,NtotEle
	 Do NN=1,NodEle(NE)
	  WDmgNod(NN,NE)=UserInitWDmg(NN,NE,NthGrpEle(NE),XNod(1,NN,NE))
	 enddo
	 enddo
	 UNod0=0.d0
	 PNod0=0.d0
	else
	 Open(51,file='DUPSolut.txt',status='old')
	 Do NE=1,NtotEle
	 Do NN=1,NodEle(NE)
	 read(1,*) (UNod0(IC,NN,NE),IC=1,NCmp),(PNod0(IC,NN,NE),IC=1,NCmp)
	 enddo
	 enddo
	 close(51)
	endif
*********************************************************************
*********************************************************************
!iterative solver!iterative solver!iterative solver!iterative solver
	NCrkMoved=-1
1001	NCrkMoved=NCrkMoved+1
	write(*,*)
	write(*,3310) NLoad,NCrkMoved,NtotEle
3310	format('NLd=',I5,', NtimesCrkAdv=',I3,', NtotE=',I4)
*================================================================
	Iterat=0
	NCallCntctMax=0
*================================================================
1101	Iterat=Iterat+1
	UNod1=UNod
	PNod1=PNod
!======================================================
	ResidMax=0.d0
	Do 1980 NE=1,NtotEle
	Do 1980 NN=1,NodEle(NE)
	NthCf0Diag=NthCfK(NE,NN,NE,NN,0)
	NthCf1Diag=NthCfK(NE,NN,NE,NN,1)
!=========================================================
!non crack!non crack!non crack!non crack!non crack!non crack!non crack
!no slave!no slave!no slave!no slave!no slave!no slave!no slave!no slave
	If(JCrkN0C1Tip2Ele(NE).eq.0.and.NthSlvMstEle(NE).eq.0) then
!-----------------------------------------------------------
!resid!resid!resid!resid!resid!resid!resid!resid!resid!resid!resid
	 Do 1201 IC=1,NCmp
	 JCrkN0C1Resid=0 
1201	 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!solving---
	 Do 1210 IC=1,NCmp
	 If(JBCU1P2(IC,NN,NE).eq.1) then
	  CfUDiag=CfUP(IC,IC,NthCf0Diag)
	  PNod(IC,NN,NE)=PNod(IC,NN,NE)-Alpha*Resid(IC)/CfUDiag
	 elseif(JBCU1P2(IC,NN,NE).eq.2) then
	  CfPDiag=CfUP(IC,IC+NCmp,NthCf0Diag)
	  UNod(IC,NN,NE)=UNod(IC,NN,NE)+Alpha*Resid(IC)/(0.5d0+CfPDiag)
	 endif
1210	 continue
12345

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -