PHOTON USE
p;;;;;;
vi 1 2 3
gr x m
gr y m
gr z m
pause
vi z
gr z m
dump gr
pause
gr off
red
set prop off
con prps z 1 fi;.001
dump prps
pause
con off
red
vec z 1 sh
dump vecs
pause
vec off
red
con u1 z 1 fi;.001
dump conu1
pause
con off
red
con p1 z 1 fi;.001
dump conp1
ENDUSE
TEXT(Parameterised inclined-plated flow :B103
Display
Air flows steadily between eight slightly-inclined steel plates,
which do not quite touch at their their ends.
Body-fitted coordinates are used; and the turbulence model is
LVEL.
Some problem-defining parameters can be set interactively, via
readvdu commands
Enddis
Groups 2 to 6. Geometry and grid
Declarations
REAL(XLEN,XCO,XI,X0,XDX,XDY,XDZ,FNX,LPERCENT,DUMMY)
REAL(YLEN,YCO,YI,Y0,YDX,YDY,YDZ,FNY,YDISP)
REAL(ZLEN,ZCO,ZI,Z0,ZDX,ZDY,ZDZ,FNZ)
REAL(INVEL)
INTEGER(IXP,IYP,IZP,IX1,IX2,IY1,IY2,IFINE)
BOOLEAN(READGRID)
grid-creation setting
READGRID=T
inlet velocity setting
INVEL=1.0
geometry settings
XLEN=1.0; YLEN=0.1; ZLEN=1.0 ; YDISP=YLEN
LPERCENT=100
grid fineness setting
IFINE=1
Interactive setting of parameters
if(iqalib.ne.0) then
+ ans=y
endif
label ask
mesgm(Steady flow of air between slightly-inclined steel plates
----------- messages to displau current options
mesg(length in flow direction is :xlen: m
mesg(thickness of air space plus plates is :ylen: m
mesg(displacement of plate ends normal to flow is :ydisp: m
mesg(length per cent is :lpercent:
mesg(inlet velocity of air is :invel:
if(ifine.eq.1) then
mesg(grid fineness is coarse
endif
if(ifine.eq.2) then
mesg(grid fineness is medium
endif
if(ifine.eq.3) then
mesg(grid fineness is finest
endif
mesg(readgrid option is :readgrid:
---------------------------- end of messages
----------------- start of interactive input
mesgm(data OK ?(Y/n)
readvdu(ans,char,y)
if(:ans:.eq.y) then
goto continue
endif ! all options will now be presented, in turn
mesg(insert new flow length?
readvdu(xlen,real,xlen)
mesg(insert new thickness?
readvdu(ylen,real,ylen)
mesg(insert new displacement?
readvdu(ydisp,real,ydisp)
mesg(insert new length percentage?
readvdu(lpercent,real,lpercent)
mesg(insert new inlet velocity?
readvdu(invel,real,invel)
mesg(change grid fineness? 1 coarse, 2 medium, 3 finest
readvdu(ifine,int,ifine)
mesg(change grid setting method? (y/N)
readvdu(ans,char,n)
if(:ans:.eq.y) then
readgrid=.not.:readgrid:
endif
goto ask
----------------- end of interactive input
label continue
----------------- start of execution
if(ifine.eq.1) then
NX=40; NY=20 ; NZ=1 ! coarse grid
endif
if(ifine.eq.2) then
NX=80;NY=40;NZ=1 ! medium grid
endif
if(ifine.eq.3) then
NX=160;NY=80;NZ=1 ! finest grid
endif
FNX=NX;FNY=NY;FNZ=NZ
BFC=T
if(readgrid) then
if(ifine.eq.1) then
readco(*/../d_earth/d_opt/d_bfc/inplib/103_1)
endif
if(ifine.eq.2) then
readco(*/../d_earth/d_opt/d_bfc/inplib/103_2)
endif
if(ifine.eq.3) then
readco(*/../d_earth/d_opt/d_bfc/inplib/103_3)
endif
else
mesg(creating bfc grid; Please wait
x0 = 0.; y0 = 0.; z0 = 0.
xdx = xlen/fnx; ydy = ylen/fny; zdz = zlen/fnz
ydx = ydisp/fnx
goto next
do ixx =0,nx
xi=ixx
if(ixx.lt.(nx/2+1)) then ! upstream half
if(ixx.lt.(nx/4+1)) then ! ist quarter
y0 = y0 + ydx
else ! 2nd quarter
y0= y0 - ydx
endif
else
if(ixx.lt.(3*nx/4+1)) then ! downstream half
y0 = y0 + ydx ! third quarter
else
y0= y0 - ydx ! fourth quarter
endif
endif
do izz =0,nz
zi=izz
do iyy =0,ny
yi=iyy
xco = x0 + xdx*xi ! corner coordinates
yco = y0 + ydy*yi
zco = z0 + zdz*zi
ixp=ixx+1
iyp=iyy+1
izp=izz+1
xc(ixp,iyp,izp)=xco
yc(ixp,iyp,izp)=yco
zc(ixp,iyp,izp)=zco
setpt(ixp,iyp,izp,xco,yco,zco) ! equivalent to the above 3 settings
enddo
enddo
enddo
endif
label next
Group 7
store(prps,enut,wdis)
#solvel
turmod(lvel)
gcv=t
Group 11
fiinit(u1)=2.0*invel
fiinit(prps)=air20
dummy=1.0 - lpercent*0.01 ! fraction cut off
dummy
dummy=dummy*fnx*0.25 ! proportion of cells cut off
dummy
ixp=dummy ! number of cells cut off
ix1=1+ixp; ix2 = nx/4 -ixp
iy1=1; iy2=ny/4
patch(slab1,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab1,prps,0.0,steel)
iy1=3*ny/4+1; iy2=ny
patch(slab2,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab2,prps,0.0,steel)
ix1=nx/4+1+ixp; ix2 = nx/2 -ixp
iy1=1; iy2=ny/4
patch(slab3,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab3,prps,0.0,steel)
iy1=3*ny/4+1; iy2=ny
patch(slab4,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab4,prps,0.0,steel)
ix1=nx/2+1+ixp; ix2 = 3*nx/4-ixp
iy1=1; iy2=ny/4
patch(slab5,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab5,prps,0.0,steel)
iy1=3*ny/4+1; iy2=ny
patch(slab6,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab6,prps,0.0,steel)
ix1=3*nx/4+1+ixp; ix2 = nx -ixp
iy1=1; iy2=ny/4
patch(slab7,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab7,prps,0.0,steel)
iy1=3*ny/4+1; iy2=ny
patch(slab8,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab8,prps,0.0,steel)
Group 13
patch(westin,west,1,1,1,ny,1,nz,1,1)
coval(westin,p1,fixflu,invel*1.189)
coval(westin,u1,onlyms,invel)
patch(eastout,east,nx,nx,1,ny,1,nz,1,1)
coval(eastout,p1,fixp,0.0)
if(lpercent.ge.100) then
patch(westin,west,1,1,1,ny/2,1,nz,1,1)
patch(eastout,east,nx,nx,1,ny/2,1,nz,1,1)
endif
#conprom
varmin(u1)=-1.e11;varmax(u1)=0.1*invel
varmin(v1)=-1.e11;varmax(v1)=0.1*invel
#maxmin
#endpause
tstswp=-1
lsweep=200
nxprin=1;nyprin=1;ixprl=10;iyprl=10
----------------- end of execution
stop