libearingfilm.for

边界远程序,加入油膜压力的接触问题求解的边界元程序

!==============================================================================
! 3D non-firction elastic contact (4节点) BEM program (Fortran)
! Maximum nodes number:1500
! Maximum element number:1500
! Maximum contact node and element number in assumed contact area: 500
!------------------------------------------------------------------------------
! main symbol instrution
!
! nodbs(i)--node total number
! numbs(i)--element total number
! nodps(i)-----预接触区上假定接触区以外的节点总数
! nodcs(i)-----预接触区上假定接触区的单元总数	
! nedps(i)-----预接触区上假定接触区以外的单元数
! nedcs(i)-----预接触区上假定接触区的节点数 
! mors(i,j)----预接触区上假定接触区的节点对应的节点号
! morp(i,j)----预接触区上假定接触区以外的节点对应的节点号
! morc(i,J)----第J个接触节点编号为 M
! merp(i,j)----编号为 K的单元是第NEDPS(i)个预接触单元
! merc(i,j)----编号为 K的单元是第NEDPS(i)个接触单元	
! MORD(i,k)----预接触区以外的节点数
! MORT(i,k)----编号为 K的节点是第NODCS(i)个接触节点 
! KOD(i,k)-----节点边界条件代码
! KODT(i,k)----单元类型码 
! KODE(i,k,l)--K单元内L节点的边界条件代码
! NODE(i,k,l)--K号接触节点是M单元的L个节点
! NORD(i,k,l)--K单元第L个节点的节点编号
! UP(numbe)----接触点对间的法向间隙
! XN(n)--------第N个节点的X坐标
! YN(n)--------第N个节点的Y坐标
! ZN(n)--------第N个节点的Z坐标
! BTX(m,l)-----单元上第L个节点X方向的面力分量(M组已知面力)
! BTY(m,l)-----单元上第L个节点Y方向的面力分量(M组已知面力)
! BTZ(m,l)-----单元上第L个节点Z方向的面力分量(M组已知面力)
! RUX(i,m,l)---单元上第L个节点X方向的位移分量(M组已知位移)
! RUY(i,m,l)---单元上第L个节点Y方向的位移分量(M组已知位移)
! RUZ(i,m,l)---单元上第L个节点Z方向的位移分量(M组已知位移)
! RTX(i,m,l)
! RTY(i,m,l)
! RTZ(i,m,l)
! NTTYP(i,j)---单元面力的组号
! NUTYP(i,j)---单元位移的组号                                    
!-------------------------------------------------------------
	Program main
	PARAMETER(M=8,N=18)
	COMMON/S1/PI,PR(2),PR1(2),PR2(2),PR3(2),PR4(2),CON(2),DSMIN2(2),
     1         DSMAX2(2),E(2),G(2)
	COMMON/S2/AXX(8),AXY(8),AXZ(8),AYX(8),AYY(8),AYZ(8),
     1            AZX(8),AZY(8),AZZ(8),BXX(8),BXY(8),BXZ(8),
     2            BYX(8),BYY(8),BYZ(8),BZX(8),BZY(8),BZZ(8)
    	COMMON/S4/HJACOB,COSAX,COSAY,COSAZ
 	COMMON/S6/FJACOB,COSBX,COSBY,COSBZ
	COMMON/S7/NODBS(2),NUMBS(2),NODPS(2),NODCS(2),NEDPS(2),NEDCS(2)
	COMMON/S8/MS(2),MD(6),NS(2),ND(2)
      common/s30/BTX(2000,8),BTY(2000,8),BTZ(2000,8),RTX(2,2000,8),
	1RTY(2,2000,8),RTZ(2,2000,8),RUX(2,2000,8),RUY(2,2000,8),
     #RUZ(2,2000,8)
	DIMENSION P0(0:M+1,0:N+1),H0(1:M+1,0:N),film(9,19),TFORCE(9,19),
 	1          PC(0:M+1,0:N+1),DM(M,N),p1(0:m,0:n),h1(0:m,0:n),
     #Atx(m,n),Aty(m,n),Atz(m,n)
     #,dtx(m,n),dty(m,n),dtz(m,n)

      DOUBLE PRECISION P0,H0,PC,DM,p1,h1,film,Atx,Aty,
     $Atz,DTX,DTY,DTZ,TFORCE
	common pii,xdjx,oilnd/c2/rjx,sp
      REAL L,KW,hmin,zs,v,jd,jd6,jdx
	dimension XN(2,3000),YN(2,3000),ZN(2,3000),kod(2,3000) 
	dimension COSBEX(2,2000,3),COSBEY(2,2000,3),COSBEZ(2,2000,3) 
	DIMENSION TITLE(20)
 	DIMENSION NTTYP(2,3000),NUTYP(2,3000),NODAL(8),U0(2000),mtyp(2)
	dimension Nord(2,3000,8),Kodt(2,3000),Kode(2,3000,8)
	dimension Mors(2,3000),Morp(2,3000),Morc(2,3000),Merp(2,3000),
     #          Merc(2,3000)
      dimension Mord(2,3000),Mort(2,3000),Node(2,3000,8)
	dimension MNB(2,3000)
	INTEGER(2)  tmphour1, tmpminute1, tmpsecond1, tmphund1
	INTEGER(2)  tmphour2, tmpminute2, tmpsecond2, tmphund2
	INTEGER(2)  tmphour3, tmpminute3, tmpsecond3, tmphund3
		  
     	open(1,file='para.in')
	OPEN(2,FILE='out.dat')
	OPEN(3,FILE='SUJU.DAT')
	open(4,file='suju1.dat')
	open(5,file='sujuzg.dat')
	OPEN(6,FILE="WWW.DAT",STATUS="OLD")
!	open(7,file='data.dat',status='old')
      open(11,file='result23.dat',status='old')
      open(101,file='chock.scr',status='old')
	open(100,file='roll.scr',status='old')
	open(102,file='FILM.DAT',status='old')
	open(103,file='TFORCE.DAT',status='old')
!	open(88,file='Course.dat',status='old')
!	open(99,file='Check.dat',status='old')
	write(11,*) ' Program name: Bearing_v2.0'
!-----------Calling system time---------------------------------------
	
	PIi=3.1415926
 	xdjx=0.0012
	oilnd0=0.137
	oilnd=oilnd0/1.5
	X=0.1674
	Y=0.0719
	D=215
	L=278
	sp=1.4
	zs=sp*60
	v=1.41
 	rjx=xdjx*d/2.0
	PXL=0.85
	WP=1.75
	SM0=0.0
	SM1=0.0

115   format(1x,'轴承的直径D=',f8.2,'mm',4x,'轴承宽度L=',f6.2,'mm',/1x,
     1    '转速n=',f4.0,'rpm',4x,'线速度v=',f5.2,'米/秒',/1x,
     1    '润滑油粘度=',f6.4,'kgf.s/m2',4x,'半径间隙c=',f6.4,'mm',/1x,
     1     '相对间隙=',f6.4)
 130  format(2x,'油膜轴承计算结果')


       DO 20 I=1,M
	   DO 10 J=1,N
	   PC(I,J)=0.0
10     CONTINUE
20     CONTINUE
21     PWJ=1.2494-1.0453*PXL
       write(*,*)pxl
       CALL CSMH(M,N,X,PXL,pwj,H0)
22     CALL REYN(M,N,X,Y,D,L,Pwj,WP,H0,PC,P0,Kl)
!	   call Hkl(H0,P0,Hk,X,Y,M,N,D)

    	DO 40 I=1,M
	    DO 30 J=1,N
		  DM(I,J)=abs(PC(I,J)-P0(I,J))
		  SM0=SM0+DM(I,J)
		  SM1=SM1+PC(I,J)
30     CONTINUE
40     CONTINUE
       DP=SM0/SM1
		IF(ABS(DP).GT.1E-3)THEN
		DO 70 I=1,M
		DO 60 J=1,N
		PC(I,J)=P0(I,J)
60      CONTINUE
70      CONTINUE
		GOTO 22
  		ENDIF
	
		call cznl(m,n,x,y,pwj,p0,w1,kw,o)
!     write(*,*)kw
!	write(*,*)o
		DF=O/KW
!	write(*,*)df
		IF(DF.lt.2E-3)THEN
		PWJ=PWJ-DF
		GOTO  22
		ENDIF

        w=oilnd*v*(l/1000)*w1/xdjx**2.0*1e-4
		IF(PXL.LT.0.85)THEN
		PXL=PXL+0.015
 		else
		PXL=PXL+0.01
		endif
!	if(pxl.gt.0.960.and.pxl.lt.0.995)then
!	write(*,*) pxl
		do 85 j=1,n
		p1(0,j)=j
		h0(0,j)=j
		h1(0,j)=j

 85		continue
		do 86 i=1,M
		p1(i,0)=i
		h0(i,0)=i
		h1(i,0)=i
 86		continue
!	write(3,200)((p1(i,j),j=0,n),i=0,M)
 		WRITE(4,200)((h0(i,J),J=0,N),i=0,m)
!	endif 
		Call Get_time(tmphour1,tmpminute1,tmpsecond1,tmphund1) 
!---------------------------------------------------------------------
	   	open(7,file='data.dat',status='old')
		READ(7,1) (TITLE(I),I=1,20)
    	    READ(7,*) MODEX,ITERL,MULTI,node1,ment1,node2,ment2
	    msi=3*max(node1,ment1,node2,ment2)
  		IF(MULTI.LE.0.OR.MULTI.GT.6) MULTI=4
 		call input(numbe,nodal,nttyp,nutyp,U0, mtyp,xn,yn,zn,msi,
     #	       cosbex,cosbey,cosbez,kod,nord,kodt,kode,mors,morp,
     #           morc,merp,merc,mord,mort,node,MNB)
! 	call input(modex,iterl,multi,numbe,nodal,nttyp,nutyp,U0,mtyp,
!    #	xn,yn,zn)
 	IF(MODEX.EQ.1) STOP	'Checing data'
      CLOSE(7)
	do i=1,2
		MS(I)=3*NODBS(I)
		NS(I)=3*NODBS(I)+MULTI*NUMBE
		ND(I)=3*(NODBS(I)-NUMBE)
	end do
!	msi=max(ms(1),ms(2))
		nsi=max(ns(1),ns(2)) 
		ksi=(3 + multi) * numbe 
!	write(*,800) msi,nsi
800   format(1x,'Coeficent matrix maximum size: ',i5,' x',i5)
	 
!---------------------------------------------------------------------
    	call matrix(iterl,multi,numbe,nodal,nttyp,nutyp,U0,mtyp,msi,nsi,
     #	        ksi,xn,yn,zn,cosbex,cosbey,cosbez,kod,nord,kodt,kode,
     #          mors,morp,morc,merp,merc,mord,mort,node,MNB,FILM,TFORCE)
!-------- Call system time ---------------------------------------
!	  do 119 i=0,8
!	  do 120 j=0,18
!	  read(4,360)(h0(i,j),j=0,18)
!	  read(102,360)(FILM(i,j),j=0,18)
!120  continue
!119   continue		   
      
	do 111 i=1,8
	do 122 j=1,18
	h0(i,j)=h0(i,j)+film(i,j)/rjx
122   continue
111   continue
	 
      DO 95 I=1,M
	DO 90 J=1,N
      P1(I,J)=2.0*SP*OILND*P0(I,J)/(((1000*XDJX)**2.0))
!	P1(I,J)=2.0*V*OILND*P0(I,J)*(D/2/1000)**2/(((1000*XDJX)**2.0))
	h1(i,j)=rjx*h0(i,j)
90    CONTINUE
95    CONTINUE

	DO 915 I=1,M
	DO 910 J=1,N
      DTX(I,J)=0
	DTY(I,J)=0
	DTZ(I,J)=0
910    CONTINUE
915    CONTINUE


	do 105 i=1,m
	jd6=6.283185/12
	i1=i+1 
	do 106 j=1,n
	k=k+1
 	k1=j+1
	if(j.le.9)then	  
	Atx(i,j)=p1(i,j)*cos(jd6)/12
	else
	Atx(i,j)=-p1(i,j)*cos(jd6)/12
	endif
	Aty(i,j)=-p1(i,j)*sin(jd6)/12 
	AtZ(i,j)=0
	jd6=jd6+jdx
      write(5,201)k,Atx(i1,j),Aty(i1,j),Atz(i1,j),
	#Atx(i,j),Aty(i,j),Atz(i,j),Atx(i,k1),Aty(i,k1),Atz(i,k1),
	#Atx(i1,k1),Aty(i1,k1),Atz(i1,k1),dtx(i1,j),dty(i1,j),dtz(i1,j),
	#dtx(i,j),dty(i,j),dtz(i,j),dtx(i,k1),dty(i,k1),dtz(i,k1),
	#dtx(i1,k1),dty(i1,k1),dtz(i1,k1)
106   continue
105	continue

	   	hmin=h1(1,1)
		do 819 i=1,m
		do 818 j=1,n
		if(h1(i,j).lt.hmin)then
		hmin=h1(i,j)
		endif
818      continue
819      continue
	    hmin=1000*hmin
!      WRITE(2,121)PXL	 
121    FORMAT(1X,F4.2)
	!	WRITE(4,200)((H1(i,J),J=0,N),i=0,m)
	!120   FORMAT(1X,'偏心率=',F4.2)
	!	WRITE(2,*)'刚性间隙H0='
	!	WRITE(2,100)(H0(1,J),J=1,N)
	! 	WRITE(3,*)'初始压力P0='
	! 	WRITE(3,100)((P0(I,J),J=1,N),I=1,M)
	!	WRITE(2,110)W,kw,kl,hmin
		WRITE(2,*)pxl,W
		WRITE(2,110)W,kw,kl,hmin
110	    FORMAT(1X,'轧制压力w=',F10.2,'吨',2x,'承载能力kw=',f6.2,
     1       /1X,'迭代次数K=',I6.0,2x,'最小油膜厚度Hmin=',f8.2,'微米')
 !	WRITE(2,150)
	    
		write(6,360) ((h1(i,j),J=0,18),i=0,8)
		write(3,200)((p1(i,j),j=0,n),i=0,M)
c	  不同偏心率下的油膜计算和接触计算耦合迭代
		if(pxl.lt.0.899)then
		goto 21
		endif

      write(6,360) ((h1(i,j),J=0,18),i=0,8)

150   FORMAT(1X)
100   FORMAT(1X,48F10.6)
200   FORMAT(1X,19F12.6)
201	FORMAT(1X,i4,6(1x,f10.4)/5x,6(1x,f10.4))
202	FORMAT(1X,3(1x,f10.6))
360   format(1X,19F12.6)


	close(99)
	CALL GET_TIME (tmphour2, tmpminute2, tmpsecond2, tmphund2)
	IF(tmphund2.GE.tmphund1) THEN
		tmphund3=tmphund2-tmphund1
	ELSE
		tmpsecond2=tmpsecond2 - 1
		tmphund3=100+tmphund2 - tmphund1
	END IF			
	IF(tmpsecond2.GT.tmpsecond1) THEN
		tmpsecond3=tmpsecond2 - tmpsecond1
	ELSE
		tmpminute2=tmpminute2 - 1
		tmpsecond3=60+tmpsecond2-tmpsecond1
	END IF
	IF(tmpminute2.GT.tmpminute1) THEN
		tmpminute3=tmpminute2 - tmpminute1
	ELSE
		tmphour2=tmphour2 - 1
		tmpminute3=60+tmpminute2 - tmpminute1
	END IF

	if((tmphour2-tmphour1).le.0.and.tmphour1.ne.0) then
		tmphour2=24 + tmphour2
	end if
		tmphour3=tmphour2 - tmphour1
	write(11,*) 
	WRITE (*, 901) tmphour3,tmpminute3,tmpsecond3,tmphund3
	WRITE (11, 901) tmphour3,tmpminute3,tmpsecond3,tmphund3
901   FORMAT(2X,'Total time consumed during executive program is',
     #       I2, ':', I2.2, ':', I2.2, ':', I2.2)
1     FORMAT(20A4)
!--------------FORMAT STATEMENTS----------------------------------------
		close(4)
		close(3)
		close(2)
		CLOSE(7)
		CLOSE(9)
		close(11)
		close(88)	
		Close(100)	
		STOP
		END

!     计算初始膜厚

	SUBROUTINE CSMH(M,N,X,PXL,pwj,H0)
	DIMENSION  H0(1:M+1,0:N),FI(M,0:N)
	DOUBLE PRECISION H0,FI
      common 	pii,xdjx,oilnd
      DO 2 I=1,M
      DO 1 J=0,N
	FI(I,J)=2*PIi/3.0+X*(J-1)-PWJ
1     CONTINUE
2     CONTINUE
      DO 5 I=1,M
      DO 5 J=0,N
	H0(I,J)=1.0+PXL*COS(FI(I,J))
5     CONTINUE
      DO 7 J=1,N
	H0(M+1,J)=H0(M,J)
7     CONTINUE
	RETURN
	END
    
!    计算刚流雷诺方程
    
	SUBROUTINE REYN(M,N,X,Y,D,L,Pwj,WP,H0,P,PS,K)
	DIMENSION P(0:M+1,0:N+1),H0(1:M+1,0:N),H3(M,N),H4(M,N),H5(M,N),
     1	      H6(M,N),H7(M,N),PT(M,N),DIFF(M,N),K3(M,N),B2(M,N),
 	1		  D2(M,N),C2(M,N),K2(M,N),F3(M,N),FI(M,0:N),
     1          PS(0:M+1,0:N+1) 
	DOUBLE PRECISION P,H0,H3,H4,H5,H6,H7,PT,DIFF,F3,K2,D2,C2,B2,FI,K3,
	1                 PS
!     PS——返回值
	REAL L
	INTEGER K
	common pii,xdjx,oilnd/c2/rjx,sp											 

1000  S0=0
	  S1=0
	
	DO 2 I=1,M
      DO 1 J=0,N
	 FI(I,J)=2*PIi/3.0+X*(J-1)-PWJ
1     CONTINUE
2     CONTINUE

      DO 20 I=1,M
      DO 10 J=1,N
      B2(I,J)=3*H0(I,J)**2.0*(H0(I,J)-H0(I,J-1))/(X*X)
      C2(I,J)=H0(I,J)**3.0/(X*X)
      D2(I,J)=(D/L)**2.0*3.0*H0(I,J)**2.0*(H0(I+1,J)-H0(I,J))/(Y*Y)
		K2(I,J)=(D/L)**2.0*H0(I,J)**3.0/(Y*Y)
		F3(I,J)=3*(H0(I,J)-H0(I,J-1))/X
      K3(I,J)=B2(I,J)+2*C2(I,J)+2*K2(I,J)+D2(I,J)
      H3(I,J)=(B2(I,J)+C2(I,J))/K3(I,J)
      H4(I,J)=(D2(I,J)+K2(I,J))/K3(I,J)
      H5(I,J)=C2(I,J)/K3(I,J)
      H6(I,J)=K2(I,J)/K3(I,J)
      H7(I,J)=F3(I,J)/K3(I,J)
10    CONTINUE
20    CONTINUE

      DO 25 I=1,M
		P(I,0)=0.0
		P(I,N+1)=0.0
25    CONTINUE
		DO 30 J=1,N
		P(0,J)=0.0					   
		P(M+1,J)=0