userj_crktip.for

两维弹性边界元程序

!For User to determine 
! whether a crack tip advances, and if advancing,
! in what direction and how long a step it advances.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!The subroutine and all other subroutines it calls apply to
! QUADRATIC ELEMENTS only. Therefore, quadratic element is demanded
! for a crack-tip element, which is internally used for new tip element 
! when crack advances.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!NETip: Nth element of tip element in global numbering
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!OUTPUT PARAMETERS:
!JAdvY1N0: =1 if the crack tip is determined to advance;
!	     =0 if not.
!OCrkAdv: orientation angle for crack advance
!ACrkAdv: step length for crack advance
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!OTHER PARAMETERS IN THE PARAMETER ZONE BELOW
!Mode1or2: =1 if opening mode by the MTS criterion;
!		 =2 if shear mode by the MSS criterion.
!SIF1/2Crit: critical SIFs for crack growth
!To learn the MTS (maximum tangential stress) and MSS (maximum shear stress)
! criteria, one may refer to following paper and original references cited therein.
! B. Yang and S. Mall, 2001, "On crack initiation mechanisms 
! in fretting fatigue," ASME Journal of Applied Mechanics 68, 76-80.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!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)
!	To calculate SIFs and T-stress at the NETip crack tip.
! SIF1, SIF2,TSs: stress intensity factors and T-stress
!SUBROUTINE IsoOCrkAdv(Mode1or2,SIF1,SIF2,OEle,SIF1Eff,SIF2Eff,OCrkAdv)
!	To compute direction of advance and corresponding effective SIFs 
!	 for a chosen mode I or II
!One may write REPLACEMENT subroutine for mixed mode crack growth instead of pure 
! opening or shear mode crack growth. When one does that, only thing to make sure 
! is not to change the interface with MAIN PROGRAM.
*********************************************************
	Subroutine UserCrkTip(
	1	NETip,JAdvY1N0,OCrkAdv,ACrkAdv,
     1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	UNod,PNod,UHomoNod,PHomoNod) !output
	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)
	Common/Example/NthExample
*================================================================
	If(NthExample.eq.1) then
	 Call Ex1_UserCrkTip(
	1	NETip,JAdvY1N0,OCrkAdv,ACrkAdv,
     1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	UNod,PNod,UHomoNod,PHomoNod)
	 return
	elseif(NthExample.eq.2) then
	 Call Ex2_UserCrkTip(
	1	NETip,JAdvY1N0,OCrkAdv,ACrkAdv,
     1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	UNod,PNod,UHomoNod,PHomoNod)
	 return
	elseif(NthExample.eq.3) then
	 return !no crack
	elseif(NthExample.eq.4) then
	 Call Ex4_UserCrkTip(
	1	NETip,JAdvY1N0,OCrkAdv,ACrkAdv,
     1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	UNod,PNod,UHomoNod,PHomoNod)
	 return
	endif
*================================================================
!PARAMETER ZONE!PARAMETER ZONE!PARAMETER ZONE!PARAMETER ZONE
!either OPENING or SHEAR crack growth mode
	Mode1or2=1
!critical SIFs in mode I and II
	SIF1Crit=1.d5
	SIF2Crit=1.d5
!crack advance step if advancing
	ACrkAdv=2.d-2
!END of PARAMETER ZONE
*================================================================
	PI=4.d0*DATAN(1.d0)
!==============================================================
!compute SIFs and T-stress
	Call IsoCrkTipSifT(
	1	NETip,SIF1,SIF2,TSs,	!output
     1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	UNod,PNod,UHomoNod,PHomoNod)
!================================================================
!compute direction of advance and corresponding effective SIFs in that direction
! for a chosen mode I or II
	Call IsoOCrkAdv(Mode1or2,SIF1,SIF2,OEle(NETip),
	1	SIF1Eff,SIF2Eff,OCrkAdv) !output
!================================================================
!determine whether crack advances
	If((Mode1or2.eq.1.and.SIF1Eff.ge.SIF1Crit).or.
     1	(Mode1or2.eq.2.and.SIF2Eff.ge.SIF2Crit)) then
	 JAdvY1N0=1
	endif
*================================================================
	write(*,1010) 
1010	format('NthGrpTip: SIF1,2&T&O',1x,I6)
	write(*,1020) NthGrpEle(NETip),SIF1,SIF2,TSs,OCrkAdv/PI*180.
	write(11,1020) NthGrpEle(NETip),XEle(1,2,NETip),XEle(2,2,NETip),
	1	SIF1,SIF2,SIF1Eff,SIF2Eff,OCrkAdv/PI*180.
1020	format(1x,I6,10(1x,f13.6))
*================================================================
	return
	end
*********************************************************
	Subroutine Ex1_UserCrkTip(
	1	NETip,JAdvY1N0,OCrkAdv,ACrkAdv,
     1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	UNod,PNod,UHomoNod,PHomoNod) !output
	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)
*================================================================
!PARAMETER ZONE!PARAMETER ZONE!PARAMETER ZONE!PARAMETER ZONE
!either OPENING or SHEAR crack growth mode
	Mode1or2=1
!critical SIFs in mode I and II
	SIF1Crit=1.d5
	SIF2Crit=1.d5
!crack advance step if advancing
	ACrkAdv=2.d-2
!END of PARAMETER ZONE
*================================================================
	PI=4.d0*DATAN(1.d0)
!==============================================================
!compute SIFs and T-stress
	Call IsoCrkTipSifT(
	1	NETip,SIF1,SIF2,TSs,	!output
     1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	UNod,PNod,UHomoNod,PHomoNod)
!================================================================
!compute direction of advance and corresponding effective SIFs in that direction
! for a chosen mode I or II
	Call IsoOCrkAdv(Mode1or2,SIF1,SIF2,OEle(NETip),
	1	SIF1Eff,SIF2Eff,OCrkAdv) !output
!================================================================
!determine whether crack advances
	If((Mode1or2.eq.1.and.SIF1Eff.ge.SIF1Crit).or.
     1	(Mode1or2.eq.2.and.SIF2Eff.ge.SIF2Crit)) then
	 JAdvY1N0=1
	endif
*================================================================
	write(*,1010) 
1010	format('NthGrpTip: SIF1,2&T&O',1x,I6)
	write(*,1020) NthGrpEle(NETip),SIF1,SIF2,TSs,OCrkAdv/PI*180.
	write(11,1020) NthGrpEle(NETip),XEle(1,2,NETip),XEle(2,2,NETip),
	1	SIF1,SIF2,SIF1Eff,SIF2Eff,TSs,OCrkAdv/PI*180.
1020	format(1x,I6,10(1x,f13.6))
*================================================================
	return
	end
*********************************************************
	Subroutine Ex2_UserCrkTip(
	1	NETip,JAdvY1N0,OCrkAdv,ACrkAdv,
     1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	UNod,PNod,UHomoNod,PHomoNod) !output
	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)
*================================================================
!PARAMETER ZONE!PARAMETER ZONE!PARAMETER ZONE!PARAMETER ZONE
!either OPENING or SHEAR crack growth mode
	Mode1or2=1
!critical SIFs in mode I and II
	SIF1Crit=1.d-5
	SIF2Crit=1.d-5
!crack advance step if advancing
	ACrkAdv=3.d-2
!END of PARAMETER ZONE
*================================================================
	PI=4.d0*DATAN(1.d0)
!==============================================================
!compute SIFs and T-stress
	Call IsoCrkTipSifT(
	1	NETip,SIF1,SIF2,TSs,	!output
     1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	UNod,PNod,UHomoNod,PHomoNod)
!================================================================
!compute direction of advance and corresponding effective SIFs in that direction
! for a chosen mode I or II
	Call IsoOCrkAdv(Mode1or2,SIF1,SIF2,OEle(NETip),
	1	SIF1Eff,SIF2Eff,OCrkAdv) !output
!================================================================
!determine whether crack advances
	If((Mode1or2.eq.1.and.SIF1Eff.ge.SIF1Crit).or.
     1	(Mode1or2.eq.2.and.SIF2Eff.ge.SIF2Crit)) then
	 JAdvY1N0=1
	endif
*================================================================
	write(*,1010) 
1010	format('NthGrpTip: SIF1,2&T&O',1x,I6)
	write(*,1020) NthGrpEle(NETip),SIF1,SIF2,TSs,OCrkAdv/PI*180.
	write(11,1020) NthGrpEle(NETip),XEle(1,2,NETip),XEle(2,2,NETip),
	1	SIF1,SIF2,SIF1Eff,SIF2Eff,TSs,OCrkAdv/PI*180.
1020	format(1x,I6,10(1x,f13.8))
*================================================================
	return
	end
*********************************************************
	Subroutine Ex4_UserCrkTip(
	1	NETip,JAdvY1N0,OCrkAdv,ACrkAdv,
     1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	UNod,PNod,UHomoNod,PHomoNod) !output
	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)
*================================================================
!PARAMETER ZONE!PARAMETER ZONE!PARAMETER ZONE!PARAMETER ZONE
!either OPENING or SHEAR crack growth mode
	Mode1or2=1
!critical SIFs in mode I and II
	If(NthGrpEle(NETip).eq.101) then
	 SIF1Crit=1.d-5
	 SIF2Crit=1.d-5
	else
	 SIF1Crit=1.d5
	 SIF2Crit=1.d5
	endif
!crack advance step if advancing
	ACrkAdv=3.d-2
!END of PARAMETER ZONE
*================================================================
	PI=4.d0*DATAN(1.d0)
!==============================================================
!compute SIFs and T-stress
	Call IsoCrkTipSifT(
	1	NETip,SIF1,SIF2,TSs,	!output
     1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	UNod,PNod,UHomoNod,PHomoNod)
!================================================================
!compute direction of advance and corresponding effective SIFs in that direction
! for a chosen mode I or II
	Call IsoOCrkAdv(Mode1or2,SIF1,SIF2,OEle(NETip),
	1	SIF1Eff,SIF2Eff,OCrkAdv) !output
!================================================================
!determine whether crack advances
	If((Mode1or2.eq.1.and.SIF1Eff.ge.SIF1Crit).or.
     1	(Mode1or2.eq.2.and.SIF2Eff.ge.SIF2Crit)) then
	 JAdvY1N0=1
	endif
*================================================================
	write(*,1010) 
1010	format('NthGrpTip: SIF1,2&T&O',1x,I6)
	write(*,1020) NthGrpEle(NETip),SIF1,SIF2,TSs,OCrkAdv/PI*180.
	write(11,1020) NthGrpEle(NETip),XEle(1,2,NETip),XEle(2,2,NETip),
	1	SIF1,SIF2,SIF1Eff,SIF2Eff,TSs,OCrkAdv/PI*180.
1020	format(1x,I6,10(1x,f13.6))
*================================================================
	return
	end