TALK=T;RUN(1,1) 737.htm

  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