土豆--冥想空间

      PROGRAM point
      
      integer m,mq,lc
      Parameter (m=500,mq=150,LC=10,PI=3.1415926)
 DIMENSION kp_c(0:(m-1),0:mq),CTA_MESH(0:m,0:mQ),
     1R_MESH(0:m,0:mQ),hc(1:(m+1)*(mq+1))
 dimension kp_a(0:(m-1),0:mq),kp_b(0:(m-1),0:mq),kp_d(0:(m-1),0:mq)
 1,kp_e(0:(m-1),0:mq)
c dimension kkp_a(0:m,0:mq),kkp_b(0:m,0:mq),kkp_d(0:m,0:mq)
c 1,kkp_e(0:(m-1),0:mq)
 DIMENSION CYLINDER_A(0:m,0:mQ)
 DIMENSION CYLINDER(0:m,0:mQ)
 real HIN(0:M),HOUT(0:M)

! difine the initial value
!need change some parameters;
 temperature=300
  Pa=101400
c      Pa=10.0E+5
 AMU=0.000001458*temperature**1.5/(temperature+110.4)
c      AMU=18.e-6
 SPEED_RPM=-1000000
  SPEED_HZ=SPEED_RPM/60
 OMEGA=SPEED_HZ*2*3.1416

 ALPHA=30.0 !sprial groove
 GAMMA=1.0 !RATIO OF ridge width AND groove width

c R2=75.0E-3
      R2=1.4E-3
 ROUT=R2
 ALAMBDA=0.4    !ALAMBDA=r1/r2
 R1=ALAMBDA*R2
c     R1=52.5E-3
 RIN=R1
 DELTA=0.1       !DELTA=H2/H0
 H0=3.0E-6      !H0:groove depth
 H2=H0*DELTA     !film depth
 HG=H2+H0        !depth from the botton OF groove to shaft botton
 HL=H2           !HL=H2--film height

!compressibility number CALCULATION
      comp=6*AMU*abs(OMEGA)*R2**2/(Pa*H0**2)

!Output important information
100   FORMAT (1x,'m= ',i4/,'mq= ',i4/,'槽数= ',i4//,'速度=
     1',e10.2/,'初始温度= ',f10.2/,'初始压力= ',f10.3/,'粘度AMU= ',
     1f11.9/,'螺旋角ALPHA= ',f6.3//,'台槽比GAMMA= ',f6.3/,
     1'半径比ALAMBDA=r1/r2= ',f6.3/,'DELTA=H2/H0= ',f6.3//,
 1'内径= ',f8.5/,'外径= ',f8.5//,
 1'H0:groove depth= ',f10.7/,'HL=H2:film height= ',f10.7/,
 1'HG=H2+H0:总深度= ',f10.7//,'轴承数= ',f16.10/)
      open (1,file='important parameters.dat')
      write (1,100) m,mq,LC,abs(SPEED_RPM),temperature,Pa,AMU,
 1ALPHA,GAMMA,DELTA,r1,r2,H0,HL,HG,comp
      close(1)

c   nondimensional parameter
 R1=R1/RIN
 R2=R2/RIN
 HG=HG/H0
 HL=HL/H0

 DCTA=pi*2./m
 DR=(R2-R1)/mq                
!计算每个单元值
 DO I=0,m
   DO J=0,mQ
     R_mesh(I,J)=r1+J*dr
          CTA_MESH(I,J)=I*DCTA
   END DO
 END DO
!转化成为直角坐标
      DO I=0,m
   DO J=0,mQ
     CYLINDER_A(I,J)=R_mesh(I,J)*cos(CTA_mesh(I,J))    
     CYLINDER(I,J)=R_mesh(I,J)*sin(CTA_mesh(I,J))       
   ENDDO
      END DO

 DCTA_G=2.*PI/(LC*(1+GAMMA))
  DCTA_L=2.*PI/(LC*(1+1./GAMMA))

      CALL POSITION(M,mQ,LC,ALPHA,HG,HL,DCTA_G,DCTA_L,KP_C,
 1DCTA,DR,R1,hin,hout,r2,kp_a,kp_b,kp_d,kp_e)

 CALL hhc(m,mq,KP_C,hin,hout,hc,HG,HL)

101   FORMAT(1X, F11.7)
      open (1,file='h.dat')
      do i=0,m
        do j=0,mQ
    t=(m+1)*j+i+1 
        if (j<mQ) then
          write (1,advance="no",fmt=101) hc(t) !不换行
     else
             write (1,101) hc(t)
     end if
     end do
      end do
      close(1)
 open (1,file='x.dat')
      do i=0,m
        do j=0,mQ
        if (j<mQ) then
          write (1,advance="no",fmt=101) CYLINDER_A(I,J) !不换行
     else
             write (1,101) R2*cos(CTA_mesh(I,mQ))
     end if
     end do
      end do
      close(1)
 open (1,file='y.dat')
      do i=0,m
        do j=0,mQ
        if (j<mQ) then
          write (1,advance="no",fmt=101) CYLINDER(I,J) !不换行
     else
              write (1,101) R2*sin(CTA_mesh(I,mQ))       
     end if
     end do
      end do 

      write(*,*) "*****************************************"
       write(*,*) ""
     
c do i=1,(m+1)*(mq+1)
cc write(*,*)'h= ',i,hc(i)
c      end do

 write(*,*) ""
      write(*,*) "*****************************************"

      stop
 end
 
!psolver
 
 SUBROUTINE POSITION(M,mQ,LC,ALPHA,HG,HL,DCTA_G,DCTA_L,KP_C,
 1DCTA,DR,R1,hin,hout,r2,kp_a,kp_b,kp_d,kp_e)
C***********************************************************************
 PARAMETER (PI=3.1415926)
      dimension kp_a(0:(m-1),0:mq),kp_b(0:(m-1),0:mq),kp_d(0:(m-1),0:mq)
 1,kp_e(0:(m-1),0:mq)
 DIMENSION CTA_start(0:2*LC)
 dimension rgl_c(0:2*LC),kp_c(0:(m-1),0:mQ)
 dimension RGL_a(0:2*LC),RGL_b(0:2*LC),RGL_d(0:2*LC),RGL_e(0:2*LC)
 dimension HIN(0:M),HOUT(0:M)
 
!所有的点都确定出具体的是在台上还是在槽上，所以定义的kP_需要是到mq，具体的系数是和kp_搭配的，
!所以也是需要到mq，但是注意的是：在后面的使用过程中我们实际上使用的不是那么多，只是一部分。
C --- Calculating start point of spiral curve

 ALPHA=ALPHA*PI/180.

 KK=0

 CTA_start(0)=0.
 CTA_start(1)=CTA_start(0)+DCTA_L  !第一个是台
 CTA_start(2*LC)=2.*PI
 
      DO KK=1,LC-1
 IK=2*KK
 CTA_start(IK)=CTA_start(IK-1)+DCTA_G
 CTA_start(IK+1)=CTA_start(IK)+DCTA_L
 ENDDO

*******************************************************************************
*                  CTA,R --- C POINT                                      *
*开始定义c点，定义的c点直接关系到计算，因为在内外边界上的值我能够确定的，所以 *
*在定义c点位置的时候，就只是使用了内部的点，把内部和外部的点排除在外。        *
*******************************************************************************
 DO I=0,m-1
 DO J=0,mq

 CTA=I*DCTA
 R=R1+J*DR
************************************************************
*下面开始找包围在c点周围的四个点                           *
c  CTA_a,R_a--Position of Point A when Point C is CTA,R*
************************************************************
 IF(I.EQ.0)THEN
 CTA_a=2*PI-0.5*DCTA
    CTA_b=2*PI-0.5*DCTA
 ELSE
 CTA_a=CTA-0.5*DCTA
    CTA_b=CTA-0.5*DCTA
 ENDIF
     
  IF(I.EQ.M)THEN
  CTA_d=0.+0.5*DCTA
      CTA_e=0.+0.5*DCTA
 ELSE
  CTA_d=CTA+0.5*DCTA
      CTA_e=CTA+0.5*DCTA
 ENDIF

      R_a=R-0.5*DR
  R_b=R+0.5*DR
   R_d=R+0.5*DR
    R_e=R-0.5*DR
********************************************************************************
*上面的程序是先找出周边各个点的极径和极角,求出的是所有的节点和周围四点的位置****
********************************************************************************

      do KJ=0,2*LC
 IF(CTA.GE.CTA_start(KJ))THEN
      RGL_C(KJ)=EXP((CTA-CTA_start(KJ))*TAN(ALPHA))
 ELSE
      RGL_C(KJ)=EXP((CTA+2.*PI-CTA_start(KJ))*TAN(ALPHA))
 ENDIF          !得到所有c点的无量纲r
 end do

 KRC=0
C    ANGLE(JJA/JJB=KKJ;JJD/JJE=KKJ-1); RADIUS(RRGL_A/RRGL_E=KKJ-1; RRGL_B/RRGL_D=KKJ)
 
      DO 121 KKJ=2*LC,1,-1
 IF(R.GT.RGL_C(KKJ).AND.R.LT.RGL_C(KKJ-1))THEN
 KGL_C=KKJ
C     KGL_a=EVEN NUMBER
        IF(MOD(KGL_C,2).EQ.0)THEN
C KP_a=1----  A is in groove area ;RRGL_a----groove length
  KP_C(I,J)=1
 RRGL_C=RGL_C(KKJ-1)
 JJC=KKJ
   ELSE
 KP_C(I,J)=0
   ENDIF
 ENDIF
 IF(R.GT.RGL_C(KKJ))KRC=1
121  CONTINUE
 
      IF(KRC.EQ.0)THEN
 DO 131 KKKJ=1,2*LC
   IF(CTA.GT.CTA_start(KKKJ-1).AND.CTA.LT.CTA_start(KKKJ))THEN
 KKGL_C=KKKJ
C     KKGL_a=EVEN NUMBER
        IF(MOD(KKGL_C,2).EQ.0)THEN
C KP_a=1----  A is in groove area
     KP_C(I,J)=1
      RRGL_C=EXP((CTA-CTA_start(KKKJ-1))*TAN(ALPHA))
   JJC=KKKJ
   ELSE
     KP_C(I,J)=0
   ENDIF

 ENDIF
131   CONTINUE
    ENDIF

C       RGL_a: spiral curve radius from different start point*
**************************************************************
      DO 11 KJ=0,2*LC

 IF(CTA_a.GE.CTA_start(KJ))THEN
      RGL_a(KJ)=EXP((CTA_a-CTA_start(KJ))*TAN(ALPHA))
 ELSE
      RGL_a(KJ)=EXP((CTA_a+2.*PI-CTA_start(KJ))*TAN(ALPHA))
 ENDIF

 RGL_b(KJ)=RGL_a(KJ)

  IF(CTA_d.GE.CTA_start(KJ))THEN
      RGL_d(KJ)=EXP((CTA_d-CTA_start(KJ))*TAN(ALPHA))
 ELSE
      RGL_d(KJ)=EXP((CTA_d+2.*PI-CTA_start(KJ))*TAN(ALPHA))
 ENDIF

    RGL_e(KJ)=RGL_d(KJ)
 11 CONTINUE
**************************************************************************************
C=================== POSITION OF A ==================== 确定各个点的位置**************
**************************************************************************************
   KRA=0
C    ANGLE(JJA/JJB=KKJ;JJD/JJE=KKJ-1); RADIUS(RRGL_A/RRGL_E=KKJ-1; RRGL_B/RRGL_D=KKJ)
 
      DO 12 KKJ=2*LC,1,-1
 IF(R_a.GT.RGL_a(KKJ).AND.R_a.LT.RGL_a(KKJ-1))THEN
 KGL_a=KKJ
C     KGL_a=EVEN NUMBER
        IF(MOD(KGL_a,2).EQ.0)THEN
C KP_a=1----  A is in groove area ;RRGL_a----groove length
  KP_a(I,J)=1
 RRGL_a=RGL_a(KKJ-1)
 JJA=KKJ
   ELSE
 KP_a(I,J)=0
   ENDIF
 ENDIF
 IF(R_a.GT.RGL_a(KKJ))KRA=1
 12 CONTINUE
 
      IF(KRA.EQ.0)THEN
 DO 13 KKKJ=1,2*LC
   IF(CTA_a.GT.CTA_start(KKKJ-1).AND.CTA_a.LT.CTA_start(KKKJ))THEN
 KKGL_a=KKKJ
C     KKGL_a=EVEN NUMBER
        IF(MOD(KKGL_a,2).EQ.0)THEN
C KP_a=1----  A is in groove area
     KP_a(I,J)=1
      RRGL_a=EXP((CTA_a-CTA_start(KKKJ-1))*TAN(ALPHA))
   JJA=KKKJ
   ELSE
     KP_a(I,J)=0
   ENDIF

 ENDIF
  13  CONTINUE
    ENDIF
C=================== POSITION OF B =====================   
  KRB=0
C    ANGLE(JJA/JJB=KKJ;JJD/JJE=KKJ-1); RADIUS(RRGL_A/RRGL_E=KKJ-1; RRGL_B/RRGL_D=KKJ)
      DO 14 KKJ=2*LC,1,-1
 IF(R_b.GT.RGL_b(KKJ).AND.R_b.LT.RGL_b(KKJ-1))THEN
 KGL_b=KKJ
C     KGL_b=EVEN NUMBER
        IF(MOD(KGL_b,2).EQ.0)THEN
C KP_b=1----  B is in groove area
  KP_b(I,J)=1
 RRGL_b=RGL_b(KKJ)
 JJB=KKJ
   ELSE
 KP_b(I,J)=0
   ENDIF
 ENDIF
 IF(R_b.GT.RGL_b(KKJ))KRB=1
 14 CONTINUE
 
      IF(KRB.EQ.0)THEN
 DO 15 KKKJ=1,2*LC
   IF(CTA_b.GT.CTA_start(KKKJ-1).AND.CTA_b.LT.CTA_start(KKKJ))THEN
 KKGL_b=KKKJ
C     KKGL_b=EVEN NUMBER
        IF(MOD(KKGL_b,2).EQ.0)THEN
C KP_b=1----  B is in groove area
     KP_b(I,J)=1
      RRGL_b=EXP((CTA_b-CTA_start(KKKJ))*TAN(ALPHA))
   JJB=KKKJ
   ELSE
     KP_b(I,J)=0
   ENDIF
 ENDIF
  15  CONTINUE
 ENDIF

C=================== POSITION OF E ==================== 
      KRE=0
C    ANGLE(JJA/JJB=KKJ;JJD/JJE=KKJ-1); RADIUS(RRGL_A/RRGL_E=KKJ-1; RRGL_B/RRGL_D=KKJ)
      DO 16 KKJ=2*LC,1,-1
 IF(R_e.GT.RGL_e(KKJ).AND.R_e.LT.RGL_e(KKJ-1))THEN
 KGL_e=KKJ
C     KGL_e=EVEN NUMBER
        IF(MOD(KGL_e,2).EQ.0)THEN
C KP_e=0----  B is in groove area
  KP_e(I,J)=1
 RRGL_e=RGL_e(KKJ-1)
  JJE=KKJ-1

   ELSE
 KP_e(I,J)=0
   ENDIF
 ENDIF
 IF(R_e.GT.RGL_e(KKJ))KRE=1
 16 CONTINUE
 
      IF(KRE.EQ.0)THEN
 DO 17 KKKJ=1,2*LC
   IF(CTA_e.GT.CTA_start(KKKJ-1).AND.CTA_e.LT.CTA_start(KKKJ))THEN
 KKGL_e=KKKJ
C     KKGL_e=EVEN NUMBER
        IF(MOD(KKGL_e,2).EQ.0)THEN
C KP_e=0----  E is in groove area
     KP_e(I,J)=1
      RRGL_e=EXP((CTA_e-CTA_start(KKKJ-1))*TAN(ALPHA))
   JJE=KKKJ-1
   ELSE
     KP_e(I,J)=0
   ENDIF
 ENDIF
  17  CONTINUE
 ENDIF

C=================== POSITION OF D ==================== 
      KRD=0
C    ANGLE(JJA/JJB=KKJ;JJD/JJE=KKJ-1); RADIUS(RRGL_A/RRGL_E=KKJ-1; RRGL_B/RRGL_D=KKJ)
      DO 18 KKJ=2*LC,1,-1
 IF(R_d.GT.RGL_d(KKJ).AND.R_d.LT.RGL_d(KKJ-1))THEN
 KGL_d=KKJ
C     KGL_d=EVEN NUMBER
        IF(MOD(KGL_d,2).EQ.0)THEN
C KP_d=0----  D is in groove area
  KP_d(I,J)=1
 RRGL_d=RGL_d(KKJ)
  JJD=KKJ-1
   ELSE
 KP_d(I,J)=0
   ENDIF
 ENDIF
 IF(R_d.GT.RGL_d(KKJ))KRD=1
 18 CONTINUE
 
      IF(KRD.EQ.0)THEN
 DO 19 KKKJ=1,2*LC
   IF(CTA_d.GT.CTA_start(KKKJ-1).AND.CTA_d.LT.CTA_start(KKKJ))THEN
 KKGL_d=KKKJ
C     KKGL_d=EVEN NUMBER
        IF(MOD(KKGL_d,2).EQ.0)THEN
C KP_d=0----  D is in groove area
     KP_d(I,J)=1
      RRGL_d=EXP((CTA_d-CTA_start(KKKJ))*TAN(ALPHA))
   JJD=KKKJ-1
   ELSE
     KP_d(I,J)=0
   ENDIF
   ENDIF
  19  CONTINUE
 ENDIF
      end do
 end do

c========================================================
!hin calculating
 DO 35 I=0,M-1
 CTA=I*DCTA

 DO JJ=0,2*LC-1
 IF(CTA.GE.CTA_start(JJ).AND.CTA.LE.CTA_start(JJ+1))THEN
 IF(MOD(JJ,2).EQ.0)THEN
     HIN(I)=HL
 ELSE
 HIN(I)=HG
 ENDIF
 ENDIF
 ENDDO

35 CONTINUE
  HIN(M)=HIN(0)

C=======================================
!HOUT CALCULATING HOUT
 DO I=0,M-1
 CTA=I*DCTA  

 DO JJ=0,2*LC-1
  CTA_OUT_JJ=LOG(R2)/TAN(ALPHA)+CTA_start(JJ)
      CTA_OUT_JJP1=LOG(R2)/TAN(ALPHA)+CTA_start(JJ+1)

 IF(CTA_OUT_JJP1.GT.2*PI.AND.CTA.LT.CTA_OUT_JJ)THEN
 CTA_OUT_JJP1=CTA_OUT_JJP1-2*PI
 CTA_OUT_JJ=CTA_OUT_JJ-2*PI
 ENDIF

 IF(CTA.GE.CTA_OUT_JJ.AND.CTA.LE.CTA_OUT_JJP1)THEN
 IF(MOD(JJ,2).EQ.0)THEN
     HOUT(I)=HL
 ELSE
 HOUT(I)=HG
 ENDIF
 ENDIF
 ENDDO

      end do
  HOUT(M)=HOUT(0)
 RETURN
 END

!hc CALCULATION the height of ALL mesh nodes.

 subroutine hhc(m,mq,KP_C,hin,hout,hc,HG,HL)
 integer m,mq
 real HIN(0:M),HOUT(0:M)
      dimension hc(1:(m+1)*(mq+1)),KP_C(0:(m-1),0:mQ)

 do i=0,m
 hc(i+1)=hin(i)
 end do
 
 j=mq
 do i=0,m
      t=(m+1)*j+i+1
 hc(t)=hout(i)
 end do   

      do i=0,m-1
   do j=1,mq-1
         t=(m+1)*j+i+1
     if(KP_C(i,j).eq.1)then
           hc(t)=HG
     else
      hc(t)=HL
     end if
   end do
 end do

 do j=1,mq-1
 t=(m+1)*j+1
      hc(t+m)=hc(t)
 end do

 return
 end