Fortran Coder

查看: 239|回复: 2

[子程序] 利用fortran进行初始与边界条件的赋予

[复制链接]

2

帖子

1

主题

0

精华

新人

F 币
13 元
贡献
5 点
发表于 2024-5-23 14:52:46 | 显示全部楼层 |阅读模式
我用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

819

帖子

2

主题

0

精华

大宗师

F 币
3862 元
贡献
2301 点
发表于 2024-5-24 10:11:11 | 显示全部楼层
得调试,需要所有代码及输入文件。

2

帖子

1

主题

0

精华

新人

F 币
13 元
贡献
5 点
 楼主| 发表于 2024-5-24 14:15:59 | 显示全部楼层
li913 发表于 2024-5-24 10:11
得调试,需要所有代码及输入文件。

因为是编写有限元软件中的子程序,这就是所有的代码了,输入文件是txt文本
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-6-25 01:42

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表