PHOTON USE
p
parphi
1 10 1
vi -x
con den1 x 1 fi;1
msg density
ENDUSE
DISPLAY
Oscillating water column, invented by Dennis Carey, 1970
Two side-by-side water columns oscillate in a structure on the
ocean floor influenced by sinusoidal pressure variations caused
by surface waves.
Water spills from the left (donor) column into the right
(receptor) column and flows thence through a turbine back to the
sea.
This input file represents its idealised behaviour by
use of In-Form, so as to:
1. calculate heights of water in each column
2. calculate corresponding density distributions
3. estimate maximum possible power generation,
4. (if findfreq = T) to determine the natural frequency
of the donor column
/////////
!.......!
!.air...!
!.......!
!..!....!
donor !..!....! receptor
column !..!....! column
!..! ! ^
!..! !
^ ! ! ! Hr
! ! !
Hd ! ! !
--- in/outflow--___!_____----outflow through turbine
^^^^^^^^^^^^^^ ocean floor ^^^^^^^^^^^^^^^^^^^
In the following Q1 file, economies of calculation and display
are achieved by treating the receptor column as standing on top
of the donor, so that the system can be regerded as one-
dimensional but time-dependent.
Explanations are provided as comments in the Q1 file.
ENDDIS
TEXT(CAREY Oscillating water column
LIBREF=737
--------------------------------------------------------------
geometric and other data
--------------------------------------------------------------
declarations
REAL(HSD,HSR,HNEUT,HWAL,WAVPE,G,OMEG,WAVAMP,WAVPER)
REAL(ADON,BDON) ! donor column coefs. for north area porosity
! npor=adon+bdon*y
REAL(AREC,BREC) ! recpt column coefs. for north area porosity
! npor=arec+brec*y
! note that non-zero bdon and brec may no longer be satisfactory
REAL(OUTCOR) ! outflo coefficient of recpt column
REAL(HDON,HREC) ! height of left (donor) and right (receptor) cols.
REAL(PCDC,PCRC,PBOTTOM,MASSAREA)
REAL(XBD,XBR) ! x-dimension of base of donor and receptor
INTEGER(NXD,NXR) ! no. of x-cells in donor and receptor
INTEGER(NYDC,NYRC) ! no. of y-cells in donor and receptor (equal)
INTEGER(NUMC,NST) ! number of cycles, and time steps per cycle
INTEGER(NWAIT) ! number of cycles to wait before opening receptor
! outlet
REAL(DTFAL)
BOOLEAN(FINDFREQ,NONRET)
settings
ZWLAST=12.5 ! z-wise width of haslar wave tank
ZWLAST=1.0 ! z-wise width of haslar wave tank
XULAST=1.0 ! nominal width
XBD=0.35 ! x-wise width of donor column at base
XBR=0.7 ! x-wise width of donor column at base
nx=1 ! number of cells in x-direction
nxd=1 ! nx in donor, = nx
nxr=1 ! nx in receptor, = nx
ny=400 ! total number of y-direction intervals
ny=200
ny=20
ny=8
nydc=ny/2 ! number of y intervals in donor
nyrc=ny/2 ! number of y intervals in receptor
HWAL=1.66 ! height of wall between donor and recpt cols
HWAL=1.5 ! height of wall between donor and recpt cols
HWAL=1.25 ! height of wall between donor and recpt cols
HNEUT = 1.0 ! height of liquid with no surface waves
HDON=2.0; HREC=2.0 ! estimated heights of both columns
G=9.81 ! gravitational acceleration
G=10.0 ! value to use when 'round numbers' needed to check
NUMC=10 ! number of cycles
NUMC=1 ! number of cycles
nst=250 ! number of time steps per cycle
lstep=nst*numc
lsweep=10 ! number of sweeps per time step
STEADY=F
PCDC=0.; PCRC=0. ! d porosity/dy for donor and receptor
ADON=xbd ! north porosity at base of donor
AREC=xbr ! north porosity at base of receptor
BDON=-PCDC/HDON
BREC=PCRC/HREC
mesgm(:title:
mesga(Two alternatives are provided for:
mesg(1. There are no surface waves to drive the
mesg( oscillation; instead the water in the donor
mesg( is raised above its equilibrium position,
mesg( but not so high as the spillover wall; it
mesg( falls to the equilibrium level, by way of
mesg( damped oscillations, revealing the natural
mesg( frequency of the system.
mesga(2. Surface waves exist, and cause the donor
mesg( column to oscillate at their frequency.
mesg( If the amplitude is sufficient, water spills
mesg( into the receptor column.
mesgm(Graphical results are best displayed by PHOTON
mesgm(file INFOROUT contains special print-out
mesgm(Find natural frequency (y/N)
readvdu(ans,char,n)
if(:ans:.eq.Y) then
findfreq=t
else
findfreq=F
mesgm(Water columns are driven by surface waves
endif
IF(FINDFREQ) THEN ! if the task is to find the natural frequency set
WAVAMP=1.e-15 ! wavamp=0, hwal=hdon, outcor=0.0
wavper=2.5 ! wavamp and wavper set to nominal values
HWAL=HDON ! oscillation takes place in donor alone
OUTCOR=0.0
HNEUT=1.0
HSD=HNEUT+0.1 ! Starting water level in donor column
! findings so far (07.06.05)
! hneut frequency
! 1.0 6/12.5
! 1.5 5/12.5
! 1.8 4.5/12.2
ELSE ! oscillation forced by surface waves
WAVAMP=0.25 ! wave amplitude (peak to trough))
WAVPER=2.5 ! wave period in seconds
HNEUT=1.0 ! neutral height of water i.e. that when wavamp=0.0
HSD=HNEUT
OUTCOR=0.5
ENDIF
NONRET=T
NWAIT=0
HSR=HNEUT ! Starting water level in recptor column
label start
integer(index)
messages to the vdu, inviting data changes
mesgm(units employed are meters and seconds
if(.not.findfreq) then
mesg(surface-wave properties
mesg(1. number of waves to simulate is :numc:
mesg(2. surface wave amplitude is :wavamp:
mesg(3. wave period is :WAVPER: seconds
endif
mesg(geometry of OWC device
mesg(4. total height of space is :HDON:
mesg(5. height of dividing wall is :hwal:
mesg(6. x-wise width of donor column is :xbd:
mesg(7. x-wise width of receptor column is :xbr:
mesg(8. z-wise width is :zwlast:
mesg(initial conditions
mesg(9. equilibrium height of donor column is :HNEUT:
mesg(10. starting height of donor column is :HSD:
mesg(11. starting height of receptor column is :HSR:
if(.not.findfreq) then
mesg(receptor outlet conditions)
mesg(12. non-return is :nonret:
mesg(13. outlet coefficient is :outcor:
mesg(14. no. cycles before opening outlet is :NWAIT:
endif
mesg(numerical parameters
mesg(15. no. of cells in y-direction (height) is :ny:
mesg(16. no. of time steps per wave cycle is :nst:
mesg(17. no. of cells in x-direction in donor is :nxd:
mesg(18. no. of cells in x-direction in recep is :nxr:
mesg(19. no. of sweeps per timestep is :lsweep:
mesg(If data ok, press return)
mesg(else indicate which you wish to change)
readvdu(index,int,0)
if(index.eq.0) then
goto end
endif
mesg(insert required value
if(index.eq.1) then
readvdu(numc,int,numc)
endif
if(index.eq.2) then
readvdu(wavamp,real,wavamp)
endif
if(index.eq.3) then
readvdu(wavper,real,wavper)
endif
if(index.eq.4) then
readvdu(hdon,real,hdon)
endif
if(index.eq.5) then
readvdu(hwal,real,hwal)
endif
if(index.eq.6) then
readvdu(xbd,real,xbd)
endif
if(index.eq.7) then
readvdu(xbr,real,xbr)
endif
if(index.eq.8) then
readvdu(zwlast,real,zwlast)
endif
if(index.eq.9) then
readvdu(HNEUT,real,HNEUT)
endif
if(index.eq.10) then
readvdu(hsd,real,hsd)
endif
if(index.eq.11) then
readvdu(hsr,real,hsr)
endif
if(index.eq.12) then
mesg(Remove non-return valve? (Y/n)
readvdu(ans,char,Y)
if(:ans:.eq.y) then
nonret=f
endif
endif
if(index.eq.13) then
readvdu(outcor,real,outcor)
endif
if(index.eq.14) then
readvdu(nwait,int,nwait)
endif
if(index.eq.15) then
readvdu(ny,int,ny)
endif
if(index.eq.16) then
readvdu(nst,int,nst)
endif
if(index.eq.17) then
readvdu(nxd,int,nxd)
endif
if(index.eq.18) then
readvdu(nxr,int,nxr)
endif
if(index.eq.19) then
readvdu(lsweep,int,lsweep)
endif
goto start
label end
Derived quantities
TLAST=WAVPER*NUMC
YVLAST=HDON+HREC
OMEG=2.*3.1416/WAVPER radians per second
messages to the vdu
mesgm(units employed are meters and seconds
mesg(total height of donor space is :HDON:
mesg(starting height of donor column is :HSD:
mesg(total height of receptor space is :HREC:
mesg(starting height of receptor column is :HSR:
mesg(height of dividing wall is :hwal:
mesg(x-wise width of donor column is :xulast:
mesg(z-wise width is :zwlast:
mesg(time for :numc: cycles is :TLAST:)
mesg(angular velocity is :omeg: radians/s
mesg(cycle frequency is :1/wavper: sec**-1
mesg(wave amplitude is :wavamp:
mesg(d porosity/dy's are :pcdc: and :pcrc:
mesg(
g
omeg
tlast
wavamp
structure of OWC example
///////
!.....!
!.....!
!.....!
!..!..!
donor !..!..! recpt
column !..!..! column
!..! ! ^
!..! !
^ ! ! ! Hr
! ! !
Hl ! ! !
!--!--!
ocean floor
lsweep=10
resfac=0.
resref(1)=0.0
endit(1) =0.0
resref(5)=0.0
endit(5) =0.0
isg21=lsweep
#unigrid ! use uniform grid
solve(p1,v1)
store(den1)
store(npor,vpor)
output(npor,p,p,n,p,p,p)
output(vpor,p,p,n,p,p,p)
output(p1,p,p,n,p,p,p)
output(v1,p,p,n,p,p,p)
fiinit(vpor)=0. ! Initial volume porosity
fiinit(npor)=0. ! Initial porosity values
! npor=0 between donor and recpt column
char(form)
patch(npordon,cell,1,nx,1,NYDC-1,1,nz,1,1)
form=adon+(bdon)*yv
(initial of NPOR at npordon is :form:) ! North porosity in donor column
patch(vpordon,cell,1,nx,1,NYDC,1,nz,1,1)
(initial of VPOR at vpordon is :form:) ! Volume porosity in donor column
patch(nporrec,cell,1,nx,NYDC+1,ny,1,nz,1,1)
form=arec+brec*(yv-HDON)
(initial of NPOR at nporrec is :form:) ! North porosity in recpt column
(initial of VPOR at nporrec is :form:) ! Volume porosity in recpt column
rho1=1.000e3
rho1=1.000e4
fiinit(den1)=rho1
gala=t
dtfal=5/wavper ! cycle-time component
dtfal=dtfal + 10/((hdon/g)**0.5) ! plus gravitational component
dtfal=1/dtfal
relax(p1,linrlx,0.5)
relax(p1,linrlx,1.0)
relax(v1,falsdt,dtfal)
varmax(v1)=0.1;varmin(v1)=-1.e11
Pressure at donor Top
patch(topd,cell,1,nx,NYDC,NYDC,1,nz,1,lstep)
coval(topd,p1,1.e7,0.0) ! Pressure at top of donor column set to zero
patch(topr,cell,1,nx,ny,ny,1,nz,1,lstep)
coval(topr,p1,fixval,0.0) ! Pressure at top of receptor column likewise
pbottom=(rho1-1.0)*g*hneut ! bottom pressure
Pressure at donor bottom
patch(donbot,cell,1,nx,1,1,1,nz,1,lstep)
form=pbottom+(rho1-1.0)*g*0.5*wavamp*sin(omeg*tim)
(source of p1 at donbot is coval(1.E7,:form:))
Pressure at recpt bottom
INTEGER(IWAIT)
IWAIT=1+NWAIT*NST
IF(NONRET) THEN
PATCH(RECBOT,OUTFLO,1,NX,NYDC+1,NYDC+1,1,NZ,IWAIT,LSTEP) ! One-way valve
ELSE
PATCH(RECBOT,CELL,1,NX,NYDC+1,NYDC+1,1,NZ,IWAIT,LSTEP) ! No one-way valve
ENDIF
(source of p1 at recbot is coval(outcor,:form:))
The next statement ensures that the file nochkflo exists, so as
to prevent satellite from making unwanted velocity patches
intrpt(w,nochkflo)
Initial pressures
patch(doninit,inival,1,1,1,nydc/2,1,1,1,1)
(initial of p1 at doninit is pbottom*$
(1 - (yg-yg[,1])/hsd))
patch(recinit,inival,1,1,nydc+1,nydc+1+nyrc/2,1,1,1,1)
(initial of p1 at recinit is pbottom*$
(1 - (yg-yg[,:nydc+1:])/hsr))
Gravity term for both donor and receptor
PATCH(GRAVITY,VOLUME,1,NX,1,NY,1,NZ,1,LSTEP)
form=-g*vpor*(DEN1-1.0) ! subtract 1.0, which is air density
(source of v1 at GRAVITY is coval(fixflu,:form:))
PATCH(DONCOL,CELL,1,nx,1,1,1,nz,1,lstep) ! location of the column
FORM=V1[1,1]*NPOR/(adon+(bdon)*HOD)
(STORE1 of VDC at DONCOL is :FORM:)
PATCH(RECCOL,CELL,1,nx,NYDC+1,NYDC+1,1,1,1,lstep) ! location of the column
FORM=V1[1,:NYRC:+1]*NPOR/(arec+brec*HOR)
!!!! inform statements which represent water-column !!!!
HND, HNR are new heights; HOD, HOR are old heights,
VDC, VRC are velocities of donor and receptor columns
The MAKE statements declare the variables and iitialise them to 0
(MAKE HND is 0) ! Water level of donor column; new
(MAKE HNR is 0) ! Water level of recpt column; new
(MAKE HOD is 0) ! Water level of donor column; previous
(MAKE HOR is 0) ! Water level of recpt column; previous
(MAKE VDC is 0) ! Water velocity of donor column at base
(MAKE VRC is 0) ! Water velocity of recpt column at base
(MAKE DHWL is 0) ! HND-HWAL
(MAKE HDMX is 0) ! maximum HND
(MAKE HRMX is 0) ! maximum HNR
(MAKE PBOT is 0) ! pbottom
(STORE1 of VRC at RECCOL is :FORM:)
(STORE1 of HOD is HSD with TSTSTR!IF(ISTEP.EQ.1)) ! 1st donor ht.
(STORE1 of HOD is HND with TSTSTR!IF(ISTEP.GT.1)) ! old donor ht.
(STORE1 of HOR is HSR with TSTSTR!IF(ISTEP.EQ.1)) ! 1st recep ht.
(STORE1 of HOR is HNR with TSTSTR!IF(ISTEP.GT.1)) ! old donor ht.
(STORE1 of HND is HOD+VDC*DT) ! New donor ht.
(STORE1 of HNR is HOR+VRC*DT) ! New recep ht.
(STORE1 of pbot is pbottom) ! New recep ht.
The above statements suffice if there is no flow between donor
and receptor; but flow can occur in three ways, namely:
1. from donor to receptor over the dividing wall, when the
receptor column height is below that of the wall; then the
donor height remains at wall height.
2. from donor to receptor when the heights are equal, greater
then the wall height and increasing.
3. from receptor to donor when the heights are equal, greater
then the wall height and decreasing.
Case 3 is notyet dealt with
Over-spilling
(STORE1 of HND is min(hdon,HOD+VDC*DT) with tstfin) ! New donor ht.
(STORE1 of HNR is HOR+VRC*DT with tstfin) ! New recep ht.
DHWL = ht. of donor column above higher of wall and receptor
(STORE1 of DHWL is MAX(0,HND-max(:hwal:,HNR)) with TSTFIN)
(STORE1 of HND is HND-DHWL with TSTFIN)
FORM=HNR+DHWL*0.5*(2*adon+(bdon)*(HND+:hwal:))/$
(arec+brec*HNR)
(STORE1 of HNR is :FORM: with TSTFIN)
maximum heights
(STORE1 of HRMX is max(HNR,HRMX) with TSTFIN)
(STORE1 of HDMX is max(HND,HDMX) with TSTFIN)
(make SPIL is 0) ! water spilled from donor to receptor
(make RATE is 0) ! rate of spill from donor to receptor
massarea=xbd*zwlast*rho1
massarea
(store1 of SPIL is SPIL+DHWL*:massarea: with TSTFIN)
(store1 of RATE is SPIL/TIM with TSTFIN)
Set density
VARMIN(DEN1)= 1.0
PATCH(COLUMND,VOLUME,1,NX,1,NYDC,1,1,1,LSTEP)
(property den1 at columnd is 1.0 with if(yg.gt.HND))
PATCH(COLUMNR,VOLUME,1,NX,NYRC+1,NY,1,1,1,LSTEP)
(property den1 at columnr is 1.0 with if(yg.gt.hnr+:HDON:))
Print-out in inforout file
(PRINT of rate is RATE)
(PRINT of hdmx is hdmx)
(PRINT of hrmx is hrmx)
idispa=1 ! for dumping to parphi
#maxabs
ixmon=1;iymon=3;izmon=1
tstswp=-1
ntprin=lstep/10
Profiles
patch(P,profil,1,1,1,1,1,1,1,lstep)
patch(P2,profil,1,1,NY/2+1,NY/2+1,1,1,1,lstep)
(print of hnd at p is hnd)
(print of hnr at p2 is hnr)
orsiz=0.2
patch(Q,profil,1,1,1,1,1,1,1,lstep)
patch(Q2,profil,1,1,NY/2+1,NY/2+1,1,1,1,lstep)
(print of dhwl at q is dhwl)
(print of spil at q2 is spil)
patch(donbottm,profil,1,1,1,1,1,1,1,lstep)
coval(donbottm,p1,0.,0.)
coval(donbottm,v1,0.,0.)
patch(recbottm,profil,1,1,NYRC+1,NYRC+1,1,1,1,lstep)
coval(recbottm,p1,0.,0.)
inifld=t
nyprin=ny/10
stop