马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？注册
x
请多指教！
该程序根据网格点和速度矢量来刻画流场，最后输出的画图文件用Tecplot打开。
!----------------------------------------------------------------------------
!* FileName : StreamLine.f90
!* Written by : Dr.Tinghe Lee
!* Purpose : Calculate position of stream lines
!* Originally complete date : 08/07/2002
!* Usage :
!* History : complete - 2002.8.7 - Dr.Tinghe Lee
!----------------------------------------------------------------------------
! Common data are saved in module for convenience
module mdlsm
implicit none
end module mdlsm
!----------------------------------------------------------------------------
program main
use mdlsm
implicit none
integer im,jm,me,ms,i,j
real dt,ei0,ej0
real,dimension( ,allocatable:: ex,ey,ez,es
real,dimension(:,,allocatable:: x,y,z,u,v,w
me = 1000
dt = 0.1
open(3,file='streamline.dat')
read(3,*) im,jm
allocate(x(im,jm),y(im,jm),z(im,jm))
allocate(u(im,jm),v(im,jm),w(im,jm))
do i = 1, im
do j = 1, jm
read(3,*) x(i,j),y(i,j),z(i,j),u(i,j),v(i,j),w(i,j)
enddo
enddo
close(3)
allocate(ex(me),ey(me),ez(me),es(me+1))
!open(1,file='tecress.dat')
!write(1,*) 'variables="x","y","z","u","v","w"'
!write(1,*) 'zone i=',im,' j=',jm,' f=point'
!do j = 1, jm
!do i = 1, im
! write(1,*) x(i,j),y(i,j),z(i,j),u(i,j),v(i,j),w(i,j)
!enddo
!enddo
!close(1)
open(2,file='tecline.dat')
write(2,*) 'variables="x","y","u","v"'
write(2,*) 'zone i=',im,' j=',jm,' f=point'
do j = 1, jm
do i = 1, im
write(2,*) x(i,j),y(i,j),u(i,j),v(i,j)
enddo
enddo
do j = 1, 20
ei0 = 1
ej0 = j
call streamline(x,y,z,u,v,w,im,jm,ei0,ej0,dt,me,ex,ey,ez,es,ms)
! output streamline for tecplot
write(2,*) 'zone i=',ms,'j=',1,' f=point'
do i = 1, ms
write(2,*) ex(i),ey(i),1,1
enddo
enddo
close(2)
deallocate(x,y,z,u,v,w,ex,ey,ez,es)
end program main
!----------------------------------------------------------------------------
! Bi-interpolation in logical field
! f00,f10,f11,f01 - function values on four corners [in]
! u,v - coordinates of logical cell [in]
! fuv - value at point (u,v) [out]
subroutine Binter(f00,f10,f11,f01,u,v,fuv)
implicit none
real f00,f10,f11,f01,u,v,fuv
fuv = (1-u)*(1-v)*f00 + u*(1-v)*f10 + u*v*f11 + (1-u)*v*f01
return
end subroutine Binter
!----------------------------------------------------------------------------
! Calculate stream line on a given surface
! Parameters:
! x,y,z - surface coordinates in phisical space [in]
! u,v,w - velocity on surface in phisical space [in]
! im,jm - dimension of surface grid [in]
! ei0,ej0 - start position of stream line [in]
! em - number of point on stream line[in]
! ei,ej - surface coordinates in logical space [out]
! es - local length along stream line[out]
! eh - streamline metric [out]
subroutine streamline(x,y,z,u,v,w,im,jm,ei0,ej0,dt,me,ex,ey,ez,es,ms)
implicit none
integer im,jm,me,i,j,m,ms,blOutBnd
real,dimension(im,jm):: x,y,z,u,v,w,ui,uj
real,dimension(im,jm):: Remi,Remj,Rems
real,dimension(me):: ex,ey,ez
real,dimension(me+1):: es
real:: ei0,ej0,ei,ej,fui,fuj,dt,ai,aj
call RightValue(x,y,z,im,jm,u,v,w,Remi,Remj,Rems)
ei = ei0
ej = ej0
ms = 0
blOutBnd = 0
do m = 1, me
if(ei<1)then
ei = 1
blOutBnd = 1
endif
if(ej<1)then
ej = 1
blOutBnd = 1
endif
if(ei>im)then
ei = im
blOutBnd = 1
endif
if(ej>jm)then
ej = jm
blOutBnd = 1
endif
i = int(ei)
j = int(ej)
ai = ei-i
aj = ej-j
ms = ms + 1
call Binter(Remi(i,j),Remi(i+1,j),Remi(i+1,j+1),Remi(i,j+1),ai,aj,fui)
call Binter(Remj(i,j),Remj(i+1,j),Remj(i+1,j+1),Remj(i,j+1),ai,aj,fuj)
ei = ei+dt*fui
ej = ej+dt*fuj
call Binter(x(i,j),x(i+1,j),x(i+1,j+1),x(i,j+1),ai,aj,ex(m))
call Binter(y(i,j),y(i+1,j),y(i+1,j+1),y(i,j+1),ai,aj,ey(m))
call Binter(z(i,j),z(i+1,j),z(i+1,j+1),z(i,j+1),ai,aj,ez(m))
es(m+1) = es(m)+dt*sqrt(fui*fui+fuj*fuj)
if(blOutBnd==1) goto 10
enddo
10 continue
return
end subroutine streamline
!----------------------------------------------------------------------------
! Function for calculating the value of vector
function FmSqrt(dx,dy,dz)
implicit none
real FmSqrt,dx,dy,dz
FmSqrt = sqrt(dx*dx+dy*dy+dz*dz)
return
end function FmSqrt
!----------------------------------------------------------------------------
! Calculate right value of difference equation
subroutine RightValue(x,y,z,im,jm,u,v,w,Remi,Remj,Rems)
implicit none
integer im,jm,i,j,in,jn
real,dimension(im,jm):: x,y,z,u,v,w,Remi,Remj,Rems
real,dimension(im,jm):: dxdi,dxdj,dxdk,dydi,dydj,dydk,dzdi,dzdj,dzdk
real:: Jack,didx,didy,didz,djdx,djdy,djdz,dkdx,dkdy,dkdz,Vel,dhdj
in = im-1
jn = jm-1
do j = 1, jm
dxdi(1,j) = x(2,j)-x(1,j)
dydi(1,j) = y(2,j)-y(1,j)
dzdi(1,j) = z(2,j)-z(1,j)
do i = 2, in
dxdi(i,j) = 0.5*(x(i+1,j)-x(i-1,j))
dydi(i,j) = 0.5*(y(i+1,j)-y(i-1,j))
dzdi(i,j) = 0.5*(z(i+1,j)-z(i-1,j))
enddo
dxdi(im,j) = x(im,j)-x(in,j)
dydi(im,j) = y(im,j)-y(in,j)
dzdi(im,j) = z(im,j)-z(in,j)
enddo
do i = 1, im
dxdj(i,1) = x(i,2)-x(i,1)
dydj(i,1) = y(i,2)-y(i,1)
dzdj(i,1) = z(i,2)-z(i,1)
do j = 2, jn
dxdj(i,j) = 0.5*(x(i,j+1)-x(i,j-1))
dydj(i,j) = 0.5*(y(i,j+1)-y(i,j-1))
dzdj(i,j) = 0.5*(z(i,j+1)-z(i,j-1))
enddo
dxdj(i,jm) = x(i,jm)-x(i,jn)
dydj(i,jm) = y(i,jm)-y(i,jn)
dzdj(i,jm) = z(i,jm)-z(i,jn)
enddo
do i = 1, im
do j = 1, jm
dxdk(i,j) = dydi(i,j)*dzdj(i,j)-dydj(i,j)*dzdi(i,j)
dydk(i,j) = dxdj(i,j)*dzdi(i,j)-dxdi(i,j)*dzdj(i,j)
dzdk(i,j) = dxdi(i,j)*dydj(i,j)-dxdj(i,j)*dydi(i,j)
enddo
enddo
do i = 1, im
do j = 1, jm
Jack = dxdi(i,j)*(dydk(i,j)*dzdj(i,j)-dydj(i,j)*dzdk(i,j))- &
dxdk(i,j)*(dydi(i,j)*dzdj(i,j)-dydj(i,j)*dzdi(i,j))+ &
dxdj(i,j)*(dydi(i,j)*dzdk(i,j)-dydk(i,j)*dzdi(i,j))
Jack = 1.0/Jack
didx = Jack*(dydk(i,j)*dzdj(i,j)-dydj(i,j)*dzdk(i,j))
didy = Jack*(dxdj(i,j)*dzdk(i,j)-dxdk(i,j)*dzdj(i,j))
didz = Jack*(dxdk(i,j)*dydj(i,j)-dxdj(i,j)*dydk(i,j))
djdx = Jack*(dydi(i,j)*dzdk(i,j)-dydk(i,j)*dzdi(i,j))
djdy = Jack*(dxdk(i,j)*dzdi(i,j)-dxdi(i,j)*dzdk(i,j))
djdz = Jack*(dxdi(i,j)*dydk(i,j)-dxdk(i,j)*dydi(i,j))
dkdx = Jack*(dydj(i,j)*dzdi(i,j)-dydi(i,j)*dzdj(i,j))
dkdy = Jack*(dxdi(i,j)*dzdj(i,j)-dxdj(i,j)*dzdi(i,j))
dkdz = Jack*(dxdj(i,j)*dydi(i,j)-dxdi(i,j)*dydj(i,j))
dhdj = sqrt(dxdj(i,j)*dxdj(i,j)+dydj(i,j)*dydj(i,j)+dzdj(i,j)*dzdj(i,j))
Vel = sqrt(u(i,j)*u(i,j)+v(i,j)*v(i,j)+w(i,j)*w(i,j))
Remi(i,j) = dhdj*(u(i,j)*didx+v(i,j)*didy+w(i,j)*didz)
Remj(i,j) = dhdj*(u(i,j)*djdx+v(i,j)*djdy+w(i,j)*djdz)
Rems(i,j) = dhdj*Vel
enddo
enddo
return
end subroutine RightValue
!----------------------------------------------------------------------------
|
IP卡
狗仔卡
发表于 2003-9-30 13:37:03
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
楼主
