Priority: Normal
   Topic: Biflap mfile
    Sent: 8/18/97
      To: Dalrymple, Robert A
    From: Campbell, Craig L ( 8/18/97 )


                                                   Time:
                       Biflap mfile
                       Subject:                            2:20 PM
                                                   Date:   8/18/97
Hello again Tony,

I believe that I have finished working with you MatLab code.  I have 
fully implemented it into a different environment, and the results 
match.  The only thing left from a verification standpoint is to 
compare some inertia and damping numbers to Hudspeth's work.

I have found a couple more things, though.  We needed to get results 
at lower waveperiods, so I spent some time massaging the evanescent 
wavenumber subprogram.  I found that if the initial value is set to a 
value slightly more than pi/2, that it will converge for very low 
periods.  In fact, the error in the trig functions becomes the 
limiting factor, not the wavenumber convergence at these low periods. 
 So, I have added the following code:

maxomega = sqrt(-G*(pi/2.0+1e-6)*tan(pi/2.0+1e-6)/h)
minperiod = 2.0*pi/maxomega
if(om>maxomega)
   fprintf('Caution: Wave periods less than %f can result in 
inaccurate evanescent wavenumbers. \n',minperiod)
end

kpp=(n*pi - (pi/2.0 - 1e-6)) ....

I believe I may have caught another typo error in wvnumev:

If k<0
  kh=abs(kh)
end

should be
if kh<0
  kh=abs(kh)
end

I also added some checks to make sure that kh was between pi/2 and 
pi.  These messages have not come up though, so I believe that the 
iteration is stable enough to prevent exceeding these bounds.

The last thing that I ran into is in nlwave in the if / else if 
sections that define the parameters for the different modes of 
operation.  The kinematic (Dalzell's) mode uses a variable called 
Omega.  I had confused this variable with your variable om.  As I had 
in my notes, om is 2 pi/period, but I did some editing and ended up 
using omega = 2 pi / period.  This changed the frequency of the 
desired waves and resulted in incorrect values.  Therefore, I changed 
the Omega variable in the if / else if sections to alpha to avoid any 
future confusion.

As I mentioned, we will also be doing a comparison of these equations 
versus Hudspeth's equations.  We will send you some of the results as 
well as the final mfiles with my modifications.

Please send me an email if the k<0 to kh<0 change is not valid.

Regards,

Craig L Campbell





Priority: Normal
   Topic: RE>Biflap mfile
    Sent: 8/18/97
      To: craig.l.campbell@mts.com
    From: Robert A. Dalrymple ( 8/18/97 )


Mail*Link SMTP          RE>Biflap mfile
Craig:  

Your  kh<0 fix to wvnumev.m is correct.   

Regarding maxomega:  I'm not sure I understand the problem.
Maxomega is 1.483e+03 for all cases, corresponding to a wave period (T) 
of 
0.0042 seconds.  That's too small for any type of wave tank motions. 

I apparently have finished the Stokes correction addition to nlwave.m  
(That is 
to say I think it is correctly worked out and programmed).  This means 
we have 
the complete second order wave theory for the dual flap wavemaker, 
allowing us 
to generate clean waves to second order.  It is attached.  It would be 
nice to 
test it in a laboratory.

(I will be out of the office until Friday)

Cheers, Tony




Priority: Normal
   Topic: MATLAB Waves
    Sent: 8/21/97
      To: Dalrymple, Robert A
    From: Campbell, Craig L ( 8/21/97 )


                                                   Time:
                       MATLAB Waves
                       Subject:                            9:44 AM
                                                   Date:   8/21/97
Hi Tony,

I was in the process of eliminating some confusion of variables, 
namely k and kf, when I discovered a couple of potential bugs.

k is the desired progressive wavenumber of height = Ht
kf is the resulting free second harmonic wavenumber of height H2

I changed the calculations of kh, denom1, g1, g2, and bet to use kf 
instead of k in the free second harmonic section.  The use of k here 
was OK, since k at this point in the code was really kf.  Then I 
added a check on the recalculated free 2nd harmonic wave height, and 
found that for the linear regression case, this recalculated wave 
height was incorrect.  After reviewing the code, I believe I may have 
discovered why.

The free second harmonic section in mode=6 (linear regression) is 
using the variable "a" to calculate b1, b2, t1, and t2.  The only 
place that a is defined gives it the value Ht/2 after the user input 
for the desired wave period.  This seems to mean that the value being 
calculated for the free 2nd harmonic waveheight is wrong, since it is 
based on the desired progressive wave.

I compared the free 2nd harmonic waveheights for the seven different 
modes, and found that the linear regression value was up to eight 
times higher than the other modes.  So, I changed the calculation of 
b1, b2, t1, and t2 to use H2/2 instead of a.  This gave values that 
were comparable with the other modes.

So, in summary, I believe the mode=6 section for the free second 
harmonic calculation should have the following (substitute H2/2 for 
a):  Is this change correct??
b1 = (H2/2)*g1/(kf*kf*sinh(kf*h))
b2 = (H2/2)*g2/(kf*kf*sinh(kf*h))

temp = bet*(g1*t1 + g2*t2)
t1 = t1*H2/temp
t2 = t2*H2/temp

Upon further review, I also found another equation for the free 
second harmonic section that uses k (desired progressive wavenumber). 
 

%progressive component due to in-phase motion
B = 2*Ap*k*(t1*ai+t2*as)/(2*kf*h + sinh(2*kf*h))

I hope that this definition is correct.  B is used to calculate 
neta2c, which is used to calculate abseta2 (free second harmonic 
waveheight), which is used to calculate Ap for the free second 
harmonic.

So the question is, should B be a function of Ap(desired prog) and 
k(desired prog) in the numerator??  Upon review of your paper based 
on my limited interpretation, it appears that the equation for B is 
correct, but I wanted to verify it.

Please send me an email with your response, positive or negative.

Regards,
Craig L Campbell
MTS Systems

