利用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
得调试，需要所有代码及输入文件。
因为是编写有限元软件中的子程序，这就是所有的代码了，输入文件是txt文本
页:
[1]