|
|
发表于 2004-10-17 18:08:14
|
显示全部楼层
[求助]FVM(有限体积法)一、二维各种格式的原代码
program kneedeep
c----67-------------------------------------------------------flowtec-|
c solves the riemann problem for st.venant equations (flat bottom)
c analytically by considering the following zones:
c (e.g. for a rarefaction-shock configuration)
c
c (L) (1) (C) (2) (R)
c ------
c \\\\\\ ----------
c \\\\\ |
c \\\\ |
c \\\ | ^ state q
c \\ | |
c \ | |
c --------------------- -----> (x/t)
c----67---------------------------------------------------------------|
c 12.02.99
c BUGS: -?-
c----67---------------------------------------------------------------|
implicit double precision (a-h,o-z)
parameter(npoimx=500)
dimension x(npoimx),u(npoimx),fh(npoimx)
logical lvacuum,l1simple,l2simple,lldry,lrdry
parameter (f23=2.d0/3.d0,f13=1.d0/3.d0,f19=1.d0/9.d0)
character*20 outfile
parameter (iunit=12,outfile='riemann.dat')
c----67---------------------------------------------------------------|
parameter (g=2.,time=1.d0)
parameter (nnwtmx=20,small=1.d-6,width=1.2d0)
c----67---------------------------------------------------------------|
c /* set the two initial states L/R of diaphragm */
write(*,10)g
write(*,*)'hL= ?'
read(*,*)fhL
write(*,*)'uL= ?'
read(*,*)uL
write(*,*)'hR= ?'
read(*,*)fhR
write(*,*)'uR= ?'
read(*,*)uR
cL=dsqrt(g*fhL)
cR=dsqrt(g*fhR)
if(fhl.ne.0.d0)then
if(fhR.eq.0.d0)then
write(*,*)'dry bed on the right'
lrdry=.TRUE.
endif
B=fhR/fhL
C=(uR-uL)/cL
else
write(*,*)'dry bed on the left'
lldry=.TRUE.
B=0.d0
C=0.d0
endif
c /* check for vacuum (i.e. dry bed) */
lvacuum=.FALSE.
if(2.d0*(1.d0+dsqrt(B)).le.C)lvacuum=.TRUE.
write(*,*)'appearance of vacuum:',lvacuum
if(lvacuum)write(*,*)'will set velocity arbitrarily to zero'
c /* check of which type are the 1- and 2-families */
l1simple=.FALSE.
l2simple=.FALSE.
if(B.gt.0.d0)
$ check1=f1(dlog(B))*dsqrt(B)
if(check1.lt.C.or.lrdry)l1simple=.TRUE.
if(B.gt.0.d0)
$ check2=f1(-dlog(b))
if(check2.lt.C.or.lldry)l2simple=.TRUE.
write(*,*)'1-family:'
if(l1simple)write(*,*)'simple wave'
if(.not.l1simple.and..not.lldry)write(*,*)'shock'
write(*,*)'2-family:'
if(l2simple)write(*,*)'simple wave'
if(.not.l2simple.and..not.lrdry)write(*,*)'shock'
c /* solve the 1-family */
if(l1simple)then
if(lvacuum.or.lrdry)then
fhC=0.d0
uc=0.d0
else
x1=-2.d0*dlog((1.d0+dsqrt(B))/2.-C/4.d0)
fhC=fhL*dexp(-x1)
CC=dsqrt(g*fhc)
uC=uL+cL*2.d0*(1.d0-dexp(-x1/2.d0))
write(*,*)'1-rarefaction: hC=',fhC,' uC=',uC
endif
else
if(l2simple)then
if(lvacuum.or.lldry)then
fhC=0.d0
uc=0.d0
else
x2=-2.d0*dlog(-C/(dsqrt(B)*4.d0)+1.d0/(2.d0*dsqrt(B))+
$ .5d0)
fhC=fhR*dexp(-x2)
CC=dsqrt(g*fhc)
uC=uR-cC*2.d0*(dexp(x2/2.d0)-1.d0)
write(*,*)'2-rarefaction: hC=',fhC,' uC=',uC
endif
else
c /* we have in fact two shocks: need to iterate */
z1=1.1d0
do n=1,nnwtmx
zold=z1
z1=zold-funcz1(zold,B,C)/dfuncz1(zold,B)
if(dabs(z1-zold).le.small)goto 111
enddo
write(*,*)'newton iteration not converged',z1,zold
stop
111 continue
fhC=z1*fhL
cc=dsqrt(fhc*g)
uC=uL-cL*(z1-1.d0)*dsqrt((z1+1.d0)/(2.d0*z1))
endif
endif
c /* we now know about the center state (c) */
c /* calculate limits of zone (1) and (2) */
c /* note: "positions" x** are actually given in "xi/t",i.e. velocities */
if(l1simple)then
xL1=uL-cL
if(.not.lvacuum.and..not.lrdry)then
x1C=uC-cC
else
c /* limit zone (1) such that vacuum is just reached through */
c /* the expansion */
x1C=uL+2.d0*cL
endif
else
c /* find shock position: calc shock speed from jump cond. */
c /* note: in this case, the zone has zero width */
if(.not.lldry)then
z=fhC/fhL
write(*,*)'1-shock: hL/hC=',z
sigma=(uC*z-uL)/(z-1.d0)
xl1=sigma
x1C=sigma
else
c /* leave this bound undefined */
endif
endif
if(l2simple)then
x2R=uR+cR
if(.not.lvacuum.and..not.lldry)then
xC2=uC+cC
else
xC2=uR-2.d0*cR
endif
else
if(.not.lrdry)then
z=fhR/fhC
write(*,*)'2-shock: hR/hC=',z
sigma=(uR*z-uC)/(z-1.d0)
x2R=sigma
xC2=sigma
else
c /* since nothing happens in the dry zone, extend it a bit from (1)(C) */
x2R=x1C+dabs(.1d0*x1c)
xC2=x2R
endif
endif
c /* get back to the case of a dry bed on the left */
if(lldry)then
x1C=xC2-dabs(.1d0*xC2)
xl1=x1c
endif
c /* set a discrete grid, a bit wider than interesting zone */
flx=(x2R-xL1)
foverlap=(width-1.d0)/2.d0
do i=1,npoimx
x(i)=(xL1-foverlap*flx)+
$ dfloat(i-1)*flx*width/dfloat(npoimx-1)
enddo
c /* set the variables at grid points */
do i=1,npoimx
if(x(i).le.xL1)then
u(i)=uL
fh(i)=fhl
elseif(x(i).le.x1C)then
u(i)=f23*x(i)+f13*(uL+2.d0*cL)
fh(i)=f19/g*(uL+2.d0*cL-x(i))**2
elseif(x(i).le.xC2)then
u(i)=uC
fh(i)=fhC
elseif(x(i).le.x2R)then
u(i)=f23*x(i)+f13*(uR-2.d0*cR)
fh(i)=f19/g*(x(i)-uR+2.d0*cR)**2
else
u(i)=uR
fh(i)=fhR
endif
enddo
c /* output */
open(iunit,file=outfile)
write(iunit,112)
do i=1,npoimx
if(fh(i).ne.0.d0)then
froude=u(i)/dsqrt(g*fh(i))
else
froude=0.d0
endif
write(iunit,*)x(i),x(i)*time,fh(i),u(i),fh(i)*u(i),
$ dsqrt(fh(i)*g),
$ u(i)+2.d0*dsqrt(fh(i)*g),u(i)-2.d0*dsqrt(fh(i)*g),
$ froude
enddo
close(iunit)
c /* format statements */
10 format('Riemann problem: (L) | (R) ; gravitational ',
$ 'accel. = ',e8.3)
112 format('#1:x/t 2:x 3:h 4:u 5:h*u 6:c 7:u+2c 8:u-2c 9:Froude')
c /* finalize */
stop
end
c----67---------------------------------------------------------------|
double precision function f1(x)
implicit double precision (a-h,o-z)
c
if(x.lt.0.d0)then
f1=-(dexp(-x)-1.d0)*dsqrt((dexp(-x)+1)/(2.d0*dexp(-x)))
elseif(x.gt.0.d0)then
f1=2.d0*(1.d0-dexp(-x/2.d0))
else
f1=0.d0
endif
end
c----67---------------------------------------------------------------|
double precision function f2(x)
implicit double precision (a-h,o-z)
c
if(x.lt.0.d0)then
f2=(dexp(x)-1.d0)*dsqrt((dexp(x)+1.d0)/(2.d0*dexp(x)))
elseif(x.gt.0.d0)then
f2=2.d0*(dexp(x/2.d0)-1.d0)
else
f2=0.d0
endif
end
c----67---------------------------------------------------------------|
double precision function funcz1(z,B,C)
implicit double precision (a-h,o-z)
c
funcz1= -(z-1.D0)*dsqrt(2.D0)*dsqrt((z+1.D0)/z)/2.D0-
$ (z/B-1.D0)*dsqrt(2.D0)*dsqrt((z/B+1.D0)/z*B)*
$ dsqrt(B)/2.D0-C
c
end
c----67---------------------------------------------------------------|
double precision function dfuncz1(z,B)
implicit double precision (a-h,o-z)
c
dfuncz1= -dsqrt(2.D0)*dsqrt((z+1.D0)/z)/2.D0-
$ (z-1.D0)*dsqrt(2.D0)/dsqrt((z+1.D0)/z)*
$ (1.D0/z-(z+1.D0)/z**2)/4.D0-1.D0/dsqrt(B)*
$ dsqrt(2.D0)*dsqrt((z/B+1.D0)/z*B)/2.D0-
$ (z/B-1.D0)*dsqrt(2.D0)/dsqrt((z/B+1.D0)/z*B)*
$ dsqrt(B)*(1.D0/z-(z/B+1.D0)/z**2*B)/4.D0
c
end
c----67--------------------------------------------------------- |
|
IP卡
狗仔卡
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
