Fortran Coder

标题: 利用fortran进行初始与边界条件的赋予 [打印本页]

作者: wangkejin 时间: 2024-5-23 14:52
标题: 利用fortran进行初始与边界条件的赋予
我用fortran编写了一段程序用于初始以及边界条件的赋予，但在实际运行过程中出现了错误，求解器不进行求解计算即壁钟时间与分析时间一直为0，而不报错。代码如下
MODULE Mandrel datas !芯轴边界条件数据
   IMPLICIT NONE
   SAVE
   real*8,dimension (320600)::InitdarhoLis,InitdaThetaCorLis,InitdazC
   $ orLis,InitdanxLis,InitdanyLis,InitdafxLis,InitdafyLis
   real*8,dimension (7,200,229)::tdarhoLis,tdazCorLis,tdaThetaCorLis,
   $ tdanxLis,tdanyLis,tdafxLis,tdafyLis
   END MODULE Mandrel datas

   subroutine usinc(f,n,ndeg,iflag)
   use Mandrel datas
   implicit real*8(a-h,o-z)
   integer n,iflag,ndeg
   real*8 f
   integer u,v
   real*8,dimension (80000,3)::cor1,cor2
   real*8,dimension (100,53)::RhoCorLis,ThetaCorLis,TempLis
   real*8,dimension (6300)::IniTempLis,IniRhoCorLis,IniThetaCorLis
   real*8,dimension (3)::cor
   integer Flag,Numrho,Numtheta,TotaNodeNum,i,j,NodeID,tdNodeID
   $  ,rFlag,tFlag
   real*8:: rho,theta

   dimension f(*)
   Numrho=53
   Numtheta=100
   TotaNodeNum=5300

   if(flag==0) then
   !!温度数据导入与处理
   open (unit=10,file='2-IniTemp.txt')
   open (unit=20,file='2-RhoCor.txt')
   open (unit=30,file='2-ThetaCor.txt')
   do i=1,5300
      read(10,*) IniTempLis(i)
      read(20,*) IniRhoCorLis(i)
      read(30,*) IniThetaCorLis(i)
   enddo
   close (unit=10)
   close (unit=20)
   close (unit=30)

   do i=1,100,1
      do j=1,53,1
      NodeID=(i-1)*53+j
         TempLis(i,j)=IniTempLis(NodeID)
         RhoCorLis(i,j)=IniRhoCorLis(NodeID)
         ThetaCorLis(i,j)=IniThetaCorLis(NodeID)
      enddo
   enddo

   !!辊芯非均匀力载荷导入与处理
   open (unit=11,file='tdarho.txt')
   open (unit=12,file='tdatheta.txt')
   open (unit=13,file='tdaz.txt')
   open (unit=14,file='tdanx.txt')
   open (unit=15,file='tdany.txt')
   open (unit=16,file='tdafx.txt')
   open (unit=17,file='tdafy.txt')

   do i=1,320600
      read(11,*) InitdarhoLis(i)
      read(12,*) InitdaThetaCorLis(i)
      read(13,*) InitdazCorLis(i)
      read(14,*) InitdanxLis(i)
      read(15,*) InitdanyLis(i)
      read(16,*) InitdafxLis(i)
      read(17,*) InitdafyLis(i)
   enddo

   close (unit=11)
   close (unit=12)
   close (unit=13)
   close (unit=14)
   close (unit=15)
   close (unit=16)
   close (unit=17)

      do i=1,229,1
      do j=1,200,1
         do k=1,7,1
      tdNodeID=(i-1)*1400+(j-1)*7+k
         tdarhoLis(k,j,i)=InitdarhoLis(tdNodeID)
         tdazCorLis(k,j,i)=InitdazCorLis(tdNodeID)
         tdaThetaCorLis(k,j,i)=InitdaThetaCorLis(tdNodeID)
         tdanxLis(k,j,i)=InitdanxLis(tdNodeID)
         tdanyLis(k,j,i)=InitdanyLis(tdNodeID)
         tdafxLis(k,j,i)=InitdafxLis(tdNodeID)
         tdafyLis(k,j,i)=InitdafyLis(tdNodeID)
      enddo
      enddo
      enddo

   Flag=1
   endif

   CALL NODVAR(0,n,cor,NQNCOMP,NQDATATYPE)

   rho= sqrt(cor(1)**2 + cor(2)**2)

   if(cor(2)>=0)then
               theta = acos(cor(1)/rho)
            elseif(cor(2)<0.and.cor(1)<0)then
               theta = acos(abs(cor(1))/rho)+3.14
            elseif(cor(2)<0.and.cor(1)>0)then
               theta = acos(-(cor(1))/rho)+3.14
   endif

   do j=1,53
      if (rho>=RhoCorLis(1,j).and.rho<=RhoCorLis(1,j+1)) then

         if (abs(rho-RhoCorLis(1,j))<=abs(rho-RhoCorLis(1,j+1)))then
            rFlag=j
         else
            rFlag=j+1
         endif
      elseif (rho<=RhoCorLis(1,1)) then
         rFlag=1
      elseif (rho>=RhoCorLis(1,Numrho)) then
         rFlag=Numrho
      endif
   enddo

   do i=1,100
      if (theta>=ThetaCorLis(i,1).and.theta<=ThetaCorLis(i+1,1)) then
         if (abs(theta-ThetaCorLis(i,1))<=abs(theta-ThetaCorLis(i+1,i)
   $ ))then
            tFlag=i
         else
            tFlag=i+1
         endif
      elseif (theta<=ThetaCorLis(1,1)) then
         tFlag=1
      elseif (theta>=ThetaCorLis(Numtheta,1)) then
         tFlag=Numtheta
      endif
   enddo

   if (cor(3)>=-2719.and.cor(3)<=-2294.or.cor(3)>=-2134.and.cor
   $ (3)<=-1674.or.cor(3)>=-1504.and.cor(3)<=-1044.or.cor(3)>
   $ =-884.and.cor(3)<=-649) then
      iflag=3
   f(1)=TempLis(tFlag,rFlag)
   else
      iflag=3
   f(1)=50
   endif

   return
   end

   subroutine forcdt(u,v,a,dp,du,time,dtime,ndeg,node,
   *                ug,xord,ncrd,iacflg,inc,ipass)

   use Mandrel datas

#ifdef _IMPLICITNONE
   implicit none
#else
   implicit logical (a-z)
#endif
c    ** Start of generated type statements **
   real*8 a, dp, dtime, du
   integer iacflg, inc, ipass, ncrd, ndeg, node
   real*8 time, u, ug, v, xord
   real*8 rho1,theta1

c    ** End of generated type statements **
   dimension u(ndeg),v(ndeg),a(ndeg),dp(ndeg),du(ndeg),ug(ndeg),
   *       xord(ncrd)
   integer i,j,k,zFlag,Numtdarho,Numtdatheta,Numtdaz,rFlag,tFlag

   include 'BCLABEL'
   Numtdarho=7
   Numtdatheta=200
   Numtdaz=229

   IF (bcname=='apply6') then

   rho1= sqrt(xord(1)**2 + xord(2)**2)

   if(xord(2)>=0)then
               theta1 = acos(xord(1)/rho1)
            elseif(xord(2)<0.and.xord(1)<0)then
               theta1 = acos(abs(xord(1))/rho1)+3.14
            elseif(xord(2)<0.and.xord(1)>0)then
               theta1 = acos(-(xord(1))/rho1)+3.14
            endif

   do i=1,229
      if (xord(3)>=tdazCorLis(1,1,i).and.xord(3)<=tdazCorLis(1,1,i+1))
   $ then
         if (abs(xord(3)-tdazCorLis(1,1,i))<=abs(xord(3)-tdazCorLis(1,
   $ 1,i+1))) then
            zFlag=i
         else
            zFlag=i+1
         endif
      elseif (xord(3)<=tdazCorLis(1,1,1)) then
         zFlag=1
      elseif (xord(3)>=tdazCorLis(1,1,Numtdaz)) then
         zFlag=Numtdaz
      endif
   enddo

      do j=1,200
      if (theta1>=tdaThetaCorLis(1,j,1).and.theta1<=tdaThetaCorLis(1,j
   $ +1,1)) then
         if (abs(theta1-tdaThetaCorLis(1,j,1))<=abs(theta1-tdaThetaCor
   $ Lis(1,j+1,1)))then
            tFlag=j
         else
            tFlag=j+1
         endif
      elseif (theta1<=tdaThetaCorLis(1,1,1)) then
         tFlag=1
      elseif (theta1>=tdaThetaCorLis(1,Numtdatheta,1)) then
         tFlag=Numtdatheta
      endif
      enddo

      do k=1,7
      if (rho1>=tdarhoLis(k,1,1).and.rho1<=tdarhoLis(k+1,1,1)) then
         if (abs(rho1-tdarhoLis(k,1,1))<=abs(rho1-tdarhoLis(k+1,1,1)))
   $ then
            rFlag=k
         else
            rFlag=k+1
         endif
      elseif (rho1<=tdarhoLis(1,1,1)) then
         rFlag=1
      elseif (rho1>=tdarhoLis(Numtdarho,1,1)) then
         rFlag=Numtdarho
      endif
   enddo

   dp(1) = tdanxLis(rFlag,tFlag,zFlag)+tdafxLis(rFlag,tFlag,zFlag)
   dp(2) = tdanyLis(rFlag,tFlag,zFlag)+tdafyLis(rFlag,tFlag,zFlag)
   endif
   return
   end

作者: li913 时间: 2024-5-24 10:11
得调试，需要所有代码及输入文件。

作者: wangkejin 时间: 2024-5-24 14:15

li913 发表于 2024-5-24 10:11
得调试，需要所有代码及输入文件。

因为是编写有限元软件中的子程序，这就是所有的代码了，输入文件是txt文本

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)