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