📄 libearingfilm.for
字号:
DO 44 MM=1,8
XQP=XQP+XK(MM)*FN(MM)
YQP=YQP+YK(MM)*FN(MM)
44 ZQP=ZQP+ZK(MM)*FN(MM)
XQP=XQP-XP
YQP=YQP-YP
ZQP=ZQP-ZP
RQ2=XQP*XQP+YQP*YQP+ZQP*ZQP
RQ1=SQRT(RQ2)
RXX2=XQP*XQP/RQ2
RYY2=YQP*YQP/RQ2
RZZ2=ZQP*ZQP/RQ2
RXY2=XQP*YQP/RQ2
RYZ2=YQP*ZQP/RQ2
RXZ2=XQP*ZQP/RQ2
CALL JACOB(NODAL,XW,YW,XK,YK,ZK,0)
FL1=CON(ID)*WX*WY*FJACOB*FJ(K)/RQ1
FL2=PR1(ID)*FL1/RQ2
FL3=FL2*(XQP*COSBX+YQP*COSBY+ZQP*COSBZ)
FL4=-3.*FL3/PR1(ID)
AXXT=-FL3+FL4*RXX2
AYYT=-FL3+FL4*RYY2
AZZT=-FL3+FL4*RZZ2
AXYT=FL4*RXY2+FL2*(XQP*COSBY-YQP*COSBX)
AXZT=FL4*RXZ2+FL2*(XQP*COSBZ-ZQP*COSBX)
AYXT=FL4*RXY2+FL2*(YQP*COSBX-XQP*COSBY)
AYZT=FL4*RYZ2+FL2*(YQP*COSBZ-ZQP*COSBY)
AZXT=FL4*RXZ2+FL2*(ZQP*COSBX-XQP*COSBZ)
AZYT=FL4*RYZ2+FL2*(ZQP*COSBY-YQP*COSBZ)
BXXT=FL1*(PR4(ID)+RXX2)
BYYT=FL1*(PR4(ID)+RYY2)
BZZT=FL1*(PR4(ID)+RZZ2)
BXYT=FL1*RXY2
BXZT=FL1*RXZ2
BYZT=FL1*RYZ2
BYXT=BXYT
BZXT=BXZT
BZYT=BYZT
C
DO 90 J=1,8
IF(J.EQ.I) GOTO 80
AXX(J)=AXX(J)+FN(J)*AXXT
AXY(J)=AXY(J)+FN(J)*AXYT
AXZ(J)=AXZ(J)+FN(J)*AXZT
AYX(J)=AYX(J)+FN(J)*AYXT
AYY(J)=AYY(J)+FN(J)*AYYT
AYZ(J)=AYZ(J)+FN(J)*AYZT
AZX(J)=AZX(J)+FN(J)*AZXT
AZY(J)=AZY(J)+FN(J)*AZYT
AZZ(J)=AZZ(J)+FN(J)*AZZT
80 BXX(J)=BXX(J)+FN(J)*BXXT
BXY(J)=BXY(J)+FN(J)*BXYT
BXZ(J)=BXZ(J)+FN(J)*BXZT
BYX(J)=BYX(J)+FN(J)*BYXT
BYY(J)=BYY(J)+FN(J)*BYYT
BYZ(J)=BYZ(J)+FN(J)*BYZT
BZX(J)=BZX(J)+FN(J)*BZXT
BZY(J)=BZY(J)+FN(J)*BZYT
BZZ(J)=BZZ(J)+FN(J)*BZZT
90 CONTINUE
170 RETURN
END
SUBROUTINE LOCALE(NUMBE,xn,yn,zn,msi,cosbex,cosbey,cosbez,
# nord,mors,kodt)
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),
1 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)
dimension COSBEX(2,numbe,3),COSBEY(2,numbe,3),COSBEZ(2,numbe,3)
dimension XN(2,msi/3),YN(2,msi/3),ZN(2,msi/3)
dimension nord(2,msi/3,8),mors(2,numbe),kodt(2,msi/3)
DIMENSION NODAL(8)
CALL INITL(1)
DO 90 I=1,2
DO 90 N=1,NUMBE
XM=0.
YM=0.
ZM=0.
XH=0.
YH=0.
ZH=0.
DO 30 K=1,NUMBS(I)
L=0
DO 10 J=1,8
10 IF(NORD(I,K,J).EQ.MORS(I,N)) L=J
IF(L.EQ.0.OR.KODT(I,K).EQ.1) GOTO 30
DO 20 J=1,8
NODAL(J)=NORD(I,K,J)
IF(NORD(I,K,J).EQ.0) GOTO 15
BXX(J)=XN(I,NORD(I,K,J))
BYY(J)=YN(I,NORD(I,K,J))
BZZ(J)=ZN(I,NORD(I,K,J))
GOTO 20
15 BXX(J)=0.
BYY(J)=0.
BZZ(J)=0.
20 CONTINUE
CALL JACOBC(NODAL,AXX(L),AYY(L),BXX,BYY,BZZ)
35 XM=XM+COSBX
YM=YM+COSBY
ZM=ZM+COSBZ
XH=XH+COSAX
YH=YH+COSAY
ZH=ZH+COSAZ
30 CONTINUE
FJACOB=SQRT(XM*XM+YM*YM+ZM*ZM)
COSBEX(I,N,3)=XM/FJACOB
COSBEY(I,N,3)=YM/FJACOB
COSBEZ(I,N,3)=ZM/FJACOB
HJACOB=SQRT(XH*XH+YH*YH+ZH*ZH)
COSBEX(I,N,1)=XH/HJACOB
COSBEY(I,N,1)=YH/HJACOB
COSBEZ(I,N,1)=ZH/HJACOB
XM=COSBEY(I,N,3)*COSBEZ(I,N,1)-COSBEZ(I,N,3)*COSBEY(I,N,1)
YM=COSBEZ(I,N,3)*COSBEX(I,N,1)-COSBEX(I,N,3)*COSBEZ(I,N,1)
ZM=COSBEX(I,N,3)*COSBEY(I,N,1)-COSBEY(I,N,3)*COSBEX(I,N,1)
FJACOB=SQRT(XM*XM+YM*YM+ZM*ZM)
COSBEX(I,N,2)=XM/FJACOB
COSBEY(I,N,2)=YM/FJACOB
COSBEZ(I,N,2)=ZM/FJACOB
90 CONTINUE
RETURN
END
SUBROUTINE JACOBC(NODAL,XW,YW,XK,YK,ZK)
COMMON/S4/HJACOB,COSAX,COSAY,COSAZ
COMMON/S6/FJACOB,COSBX,COSBY,COSBZ
DIMENSION NODAL(8),DFXW(8),DFYW(8),XK(8),YK(8),ZK(8)
C
DO 10 J=5,8
DFXW(J)=0.
10 DFYW(J)=0.
IF(NODAL(5).NE.0) DFXW(5)=-XW*(1.-YW)
IF(NODAL(6).NE.0) DFYW(6)=-YW*(1.+XW)
IF(NODAL(7).NE.0) DFXW(7)=-XW*(1.+YW)
IF(NODAL(8).NE.0) DFYW(8)=-YW*(1.-XW)
IF(NODAL(5).NE.0) DFYW(5)=-0.5*(1.-XW*XW)
IF(NODAL(6).NE.0) DFXW(6)= 0.5*(1.-YW*YW)
IF(NODAL(7).NE.0) DFYW(7)= 0.5*(1.-XW*XW)
IF(NODAL(8).NE.0) DFXW(8)=-0.5*(1.-YW*YW)
DFXW(1)=-0.25*(1.-YW)-0.5*(DFXW(5)+DFXW(8))
DFXW(2)= 0.25*(1.-YW)-0.5*(DFXW(5)+DFXW(6))
DFXW(3)= 0.25*(1.+YW)-0.5*(DFXW(6)+DFXW(7))
DFXW(4)=-0.25*(1.+YW)-0.5*(DFXW(7)+DFXW(8))
DFYW(1)=-0.25*(1.-XW)-0.5*(DFYW(5)+DFYW(8))
DFYW(2)=-0.25*(1.+XW)-0.5*(DFYW(5)+DFYW(6))
DFYW(3)= 0.25*(1.+XW)-0.5*(DFYW(6)+DFYW(7))
DFYW(4)= 0.25*(1.-XW)-0.5*(DFYW(7)+DFYW(8))
XD1=DFXW(1)*XK(1)+DFXW(2)*XK(2)+DFXW(3)*XK(3)+DFXW(4)*XK(4)
1 +DFXW(5)*XK(5)+DFXW(6)*XK(6)+DFXW(7)*XK(7)+DFXW(8)*XK(8)
YD1=DFXW(1)*YK(1)+DFXW(2)*YK(2)+DFXW(3)*YK(3)+DFXW(4)*YK(4)
1 +DFXW(5)*YK(5)+DFXW(6)*YK(6)+DFXW(7)*YK(7)+DFXW(8)*YK(8)
ZD1=DFXW(1)*ZK(1)+DFXW(2)*ZK(2)+DFXW(3)*ZK(3)+DFXW(4)*ZK(4)
1 +DFXW(5)*ZK(5)+DFXW(6)*ZK(6)+DFXW(7)*ZK(7)+DFXW(8)*ZK(8)
XD2=DFYW(1)*XK(1)+DFYW(2)*XK(2)+DFYW(3)*XK(3)+DFYW(4)*XK(4)
1 +DFYW(5)*XK(5)+DFYW(6)*XK(6)+DFYW(7)*XK(7)+DFYW(8)*XK(8)
YD2=DFYW(1)*YK(1)+DFYW(2)*YK(2)+DFYW(3)*YK(3)+DFYW(4)*YK(4)
1 +DFYW(5)*YK(5)+DFYW(6)*YK(6)+DFYW(7)*YK(7)+DFYW(8)*YK(8)
ZD2=DFYW(1)*ZK(1)+DFYW(2)*ZK(2)+DFYW(3)*ZK(3)+DFYW(4)*ZK(4)
1 +DFYW(5)*ZK(5)+DFYW(6)*ZK(6)+DFYW(7)*ZK(7)+DFYW(8)*ZK(8)
G1=YD1*ZD2-ZD1*YD2
G2=ZD1*XD2-XD1*ZD2
G3=XD1*YD2-YD1*XD2
FJACOB=SQRT(G1*G1+G2*G2+G3*G3)
COSBX=G1/FJACOB
COSBY=G2/FJACOB
COSBZ=G3/FJACOB
H1=YD2*G3-ZD2*G2
H2=ZD2*G1-XD2*G3
H3=XD2*G2-YD2*G1
HJACOB=SQRT(H1*H1+H2*H2+H3*H3)
COSAX=H1/HJACOB
COSAY=H2/HJACOB
COSAZ=H3/HJACOB
RETURN
END
SUBROUTINE CONTACT(ITERL,MULTI,numbe,U0,nttyp,nutyp,
# A,P,A1,P1,AC,PC,msi,nsi,ksi,xn,yn,zn,cosbex,cosbey,cosbez,
# kod,nord,kodt,kode,mors,morp,morc,merp,merc,mord,mort,node,
# MNB,FILM,TFORCE)
COMMON/S1/PI,PR(2),PR1(2),PR2(2),PR3(2),PR4(2),
1 CON(2),DSMIN2(2),DSMAX2(2),E(2),G(2)
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)
dimension COSBEX(2,numbe,3),COSBEY(2,numbe,3),COSBEZ(2,numbe,3)
dimension A(msi,nsi),P(msi),AC(2,3*numbe,ksi),PC(2,3*numbe)
dimension A1(msi,nsi),p1(msi)
COMMON/S30/BTX(2000,8),BTY(2000,8),BTZ(2000,8),
1 RTX(2,2000,8),RTY(2,2000,8),RTZ(2,2000,8),
1 RUX(2,2000,8),RUY(2,2000,8),RUZ(2,2000,8)
dimension Y(msi),TN(msi/3),UP(numbe),UN(2,numbe)
dimension C(3*numbe,ksi),B(3*numbe),X(8*numbe)
dimension NTTYP(2,msi/3),NUTYP(2,msi/3),FILM(9,19),TFORCE(9,19)
DIMENSION U0(numbe),psign(2,numbe),xt(3000)
DIMENSION XX(3000),pb(6*numbe),AF2(6*numbe,4*numbe)
DIMENSION XX1(3000),Y1(msi),TN1(msi/3),xt1(3000)
dimension XN(2,msi/3),YN(2,msi/3),ZN(2,msi/3)
dimension Nord(2,msi/3,8),KODT(2,msi/3),Kode(2,msi/3,8)
dimension kod(2,msi/3)
dimension Mors(2,numbe),Morp(2,numbe),Morc(2,numbe),Merp(2,numbe),
# Merc(2,numbe)
dimension Mord(2,msi/3),Mort(2,msi/3),Node(2,msi/3,8)
DIMENSION MNB(2,3000)
DOUBLE PRECISION film,TFORCE
!-----------------------------------------------------------------
! WRITE(99,*) 'THE NUMBER OF CONTACT NODES IS ',NUMBE
C
C DEFINE THE LOCAL COORDINATE SYSTEMS AND COMPUTE NORMAL GAP U0
C BETWEEN THE TWO CONTACT BOUNDARIES.
C
CALL LOCAL3(NUMBE,xn,yn,zn,msi,cosbex,cosbey,cosbez,
# nord,kodt,mors)
WRITE(*,29)
WRITE(11,29)
do 300 LM=1,numbe
up(LM)=U0(LM)
WRITE(*,30) LM,MORS(1,LM),MORS(2,LM),UP(LM)
WRITE(11,30) LM,MORS(1,LM),MORS(2,LM),UP(LM)
300 CONTINUE
ITER=0
540 ITER=ITER+1
77 WRITE(*,31) ITER,NEDPS(2),NEDCS(2)
C IF(ITER.GT.ITERL) GOTO 640
MD(1)=3*NUMBE
MD(2)=3*NODPS(2)
MD(3)=3*NODCS(2)
MD(4)=2*NODCS(2)
MD(5)=3*NUMBE+NODCS(2)
MD(6)=3*NUMBE-NODCS(2)
REWIND 10
DO 560 I=1,2
550 CALL ASSEM(I,NUMBE,MULTI,AC,PC,msi,ksi,C,B,kod,mors,mort,node,
# nord,kode,kodt,UP)
c WRITE(*,*)' THE SUBROUTION OF ASSEM HAS BEEN FINISHED'
! WRITE(88,*)' THE SUBROUTION OF ASSEM HAS BEEN FINISHED'
CALL SOLVE1(MD(1),MD(5),MD(6),ksi,C,B,numbe)
DO 560 J=1,MD(1)
K=(I-1)*MD(1)+J
PB(K)=B(J)
DO 560 N=1,MD(5)
AF2(K,N)=C(J,N)
560 CONTINUE
c WRITE(*,*) 'THE FIRST SAVING STAGE HAS BEEN FINISHED'
! WRITE(88,*) 'THE FIRST SAVING STAGE HAS BEEN FINISHED'
CALL SOLVE(PB,numbe,AF2,ksi,C,B,X)
!
! PRINT CONTACT COUPLES,BOUNDARY ELEMENT DISPLACEMENTS AND PRESSURE
!
600 WRITE(*,32) NEDCS(2),ITER
DO 630 I=1,2
WRITE(*,33) I
DO 630 K=1,NUMBE
609 L=MORT(I,MORS(I,K))
IF(KOD(I,MORS(I,K)).EQ.10) GOTO 610
KZ=(I-1)*MD(6)+3*L
KY=KZ-1
KX=KZ-2
UN(i,k)=X(KZ)
PSIGN(I,K)=0.
GOTO 620
610 KZ=2*MD(6)+L
KY=(I-1)*MD(6)+MD(2)+2*L
KX=KY-1
LZ=6*NUMBE-NODCS(2)+L
UN(i,k)=REAL(I-1)*UP(K)+REAL(3-2*I)*X(KZ)
PSIGN(I,K)=2.*G(1)*X(LZ)
620 continue
IF(I.EQ.1)THEN
IUZ=3*K
IUY=IUZ-1
IUX=IUZ-2
XX1(IUZ)=UN(i,k)
XX1(IUY)=X(KY)
XX1(IUX)=X(KX)
Xt1(K)=X(LZ)
KK=K+ND(1)/3
KUX=3*KK-2
KUY=3*KK-1
KUZ=3*KK
Y1(KUX)=XX1(IUX)
Y1(KUY)=XX1(IUY)
Y1(KUZ)=XX1(IUZ)
TN1(KK)=2.0*G(1)*X(LZ)
ENDIF
IF(I.EQ.2)THEN
IUZ=3*K
IUY=IUZ-1
IUX=IUZ-2
XX(IUZ)=UN(i,k)
XX(IUY)=X(KY)
XX(IUX)=X(KX)
Xt(K)=X(LZ)
KK=K+ND(2)/3
KUX=3*KK-2
KUY=3*KK-1
KUZ=3*KK
Y(KUX)=XX(IUX)
Y(KUY)=XX(IUY)
Y(KUZ)=XX(IUZ)
TN(KK)=2.0*G(2)*X(LZ)
END IF
C write(112,34)MORS(I,K),X(KX),X(KY),UN(i,k),PSIGN(I,K)
630 WRITE(*,34) MORS(I,K),X(KX),X(KY),UN(i,k),PSIGN(I,K)! CALL GAP(NUMBE,U0,PSIGN,msi,kod,mors,UP)
call valid(numbe,X,msi,nord,kod,merc,merp,mors,morp,
# kode,kodt,mort,morc,mord,up)
C if(nedcs(2).ne.nedcs(1)) then
C goto 540
C else
WRITE(11,31) ITER,NEDPS(2),NEDCS(2)
WRITE(11,32) NEDCS(2),ITER
WRITE(11,33) I
do i=1,2
do k=1,numbe
KZ=(I-1)*MD(6)+3*MORT(I,MORS(I,K))
ky=kz-1
kx=kz-2
WRITE(11,34) MORS(I,K),X(KX),X(KY),UN(i,k),PSIGN(I,K)
end do
end do
C end if
C K1=1
C write(5,971)((j-1)*10,j=1,18)
C DO 221 I=1,20
C write(5,972)k1,(UN(1,(K-1)*21+I),K=1,17)
C K1=K1+10
C221 CONTINUE
C K1=1
C write(6,971)((j-1)*10,j=1,18)
C DO 222 I=1,20
C write(6,972)k1,(PSIGN(1,(K-1)*21+I),K=1,17)
C K1=K1+10
C222 CONTINUE
971 FORMAT(1X,21(I4,' '))972 FORMAT(1X,I4,' ',21(E12.6,' ')) CALL ASSEMM(1,NUMBE,MULTI,XX1,xt1,A1,P1,msi,nsi,node,nord,mors
#,kodt)
CALL SOLVE5(G,A1,P1,msi,nsi,Y1,kod,MNB)
CALL BSTRESS1(numbe,nttyp,nutyp,msi,Y1,TN1,xn,yn,zn,
# cosbex,cosbey,cosbez,nord,mord,kode,kodt)
CALL ASSEMM(2,NUMBE,MULTI,XX,xt,A,P,msi,nsi,node,nord,mors,kodt)
CALL SOLVE4(G,A,P,msi,nsi,Y,kod,MNB)
CALL BSTRESS(numbe,nttyp,nutyp,msi,Y,TN,xn,yn,zn,
# cosbex,cosbey,cosbez,nord,mord,kode,kodt,FILM,TFORCE)
!---------------------------------------------------------------------
29 FORMAT(//,47H ASSUMED PRE_CONTACT OR CONTACT ELEMENT INITIAL,
1 6H DATA.,//,46H CONTACT COUPLE ELEMENT(BODY_1) ELEM,
2 28HENT(BODY_2) INITIAL GAP_UD,/)
30 FORMAT(3I20,F15.4)
31 FORMAT(1X,8H ITER = ,I5,/,9H NEDPS = ,I5,/,9H NEDCS = ,I5,/)
32 FORMAT(1X,33H SOLUTION OF THE CONTACT PROBLEM.,/,11H NUMBERS OF,
1 19H CONTACT ELEMENTS =,I5,/,23H NUMBERS OF ITERATION =,I5,/)
33 FORMAT(/," DISPLACEMENT AND STRESS IN THE CONTACT AREA_BODY",
1 1H(,I2,1H),//," NODAL US1 US2 ",
2 "UN PSIGN",/)
34 FORMAT(I10,4E14.5)
640 RETURN
END
SUBROUTINE VALID(NUMBE,X,msi,nord,kod,merc,merp,mors,morp,
# kode,kodt,mort,morc,mord,UP)
COMMON/S7/NODBS(2),NUMBS(2),NODPS(2),NODCS(2),NEDPS(2),NEDCS(2)
dimension Nord(2,msi/3,8),KODT(2,msi/3),Kode(2,msi/3,8)
dimension kod(2,msi/3)
dimension Mors(2,numbe),Morp(2,numbe),Morc(2,numbe),Merp(2,numbe),
# Merc(2,numbe),Mort(2,msi/3),Mord(2,msi/3)
DIMENSION NODAL(8),XK(8),YK(8),ZK(8)
dimension X(8*numbe),UP(numbe)
C
NEDPS(1)=NEDPS(2)
NEDCS(1)=NEDCS(2)
NEDCS(2)=0
DO 45 N=1,NEDCS(1)
DO 15 J=1,8
NODAL(J)=NORD(1,MERC(1,N),J)
IF(NODAL(J).EQ.0) GOTO 10
KZ=6*NUMBE-NODCS(2)+MORT(1,NODAL(J))
XK(J)=X(KZ)
GOTO 15
10 XK(J)=0.
15 CONTINUE
YK(1)=XK(1)+XK(2)+XK(3)+XK(4)+XK(5)+XK(6)+XK(7)+XK(8)
if(yk(1).LT.0.0) go to 30
NEDPS(2)=NEDPS(2)+1
DO 20 I=1,2
20 MERP(I,NEDPS(2))=MERC(I,N)
GOTO 45
30 NEDCS(2)=NEDCS(2)+1
DO 40 I=1,2
40 MERC(I,NEDCS(2))=MERC(I,N)
45 CONTINUE
IF(NEDCS(2).NE.NEDCS(1)) GOTO 90
IF(NEDPS(2).EQ.0) GOTO 90
NEDPS(2)=0
DO 80 N=1,NEDPS(1)
DO 60 I=1,2
DO 55 J=1,8
NODAL(J)=NORD(I,MERP(I,N),J)⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -