hi @MShabara

Can you please provide an example of the WAMIT input files that will generate QTFs? It would be great to have some means to test a few versions of these files, and especially when ready to release since searching for "QTF" doesn't come up with anything in the WAMIT user manual.

@dforbush2 ,
Hi Dom,
WAMIT supports QTF only for v6.4

1. Chapter 11 in WAMIT v6.4S manual is dedicated to the QTFs "VERSION 6.4S (Second-order module)"
   https://www.wamit.com/manual6.4/Chap11.pdf
2. WAMIT files has an example that solves the QTFs (second order forces).
3. The current WAMIT output files versions that this pull request supports are OPTN.10s/d, OPTN.qcs/d, OPTN.11s/d and OPTN.12s/d
4. Future pull requests will include: (a) BIMIO examples will include the QTFs generation for WAMIT and NEMOH solvers.
   (b) Example in the WEC-Sim Applications will be added for the QTFs.

@salhus ah i see good to know. In that case...I'll see if i can get an install here of that version. If not, it would be good to check that the file format is what you expect for multiple bodies, additional GBMs or free surface modes, etc. on the BEMIO side. That's more what I meant @MShabara. If it can't handle all of those things yet thats fine, we just need to understand what its current limitations there are for documentation.

@MShabara whats the recommended way to test this? There isn't a modified wecSimInputFile that shows how to call this. I can try to piece it together but that might be less efficient than making a test example

I've sent over an example of the output files for WAMIT to utilize. If you're currently testing out the modified functions, feel free to experiment with them using the RM3 example. Simply add the following line to the wecSimInputFile:

body(i).QTFs = 1; % 0: no QTFs, 1: Full QTFs calculations (more calculation options will be added in the future).

Additionally, you can find a rundown of the major limitations of this functionality listed in the new secondOrderExt.m file, located within the WEC-Sim\source\functions\simulink\model directory.

There is a lot here and thank you @MShabara for making it so readable and well-commented already.
Some minor comments:

1. bodyClass line 598: Fmax = 10 seems hard-coded. I don't think there will ever be a need to get higher than 10 Hz(? unit ambiguous), but it is there a good reason to hard-code this limit? Looking ahead this probably isn't intentional because the frequency vector f is calculated from this on line 609 and this should actually depend on waveAmpTime(:,2).
2. Similarly, the index_min, index_max etc on bodyClass line 626 seem a bit arbitrary (e.g., 0.8, 2 * omega_max * 1.2). If these are integral to the method ignore me, but it does seem like these thresholds would be tunable or specific to the frequency spacing of the BEM run.
3. I'm very wary of the 'symmetric' flag here on the IFFTs because it can cover up minor errors in indexing by forcing the ifft to have strictly real parts. it looks like everywhere is fine to my read but I'm interested if testing w/ symmetric removed is warranted.

And I know I've asked before and probably already have this somewhere, but if you could please provide a link to the citation of the methodology you're implementing here it would be a great reference (email is fine if unpublished). Will continue to review later.

Thanks @dforbush2, these are really good questions/comments:

1. You are right, it should be F_max = 1/(waveAmpTime(2,1) - waveAmpTime(1,1));
   Fmax = 10 was for a specific test case I was working on which corresponds to simu.dt = 0.1; and should be removed.

2. The standard practice in vibration analysis is to use a margin of 20% around the signal of interest, making this margin smaller would affect the result's accuracy, but further increasing this margin should has minimal effect on the solution but increases the computational time.

3. The code underwent testing without the symmetric flag initially. Symmetry was manually verified for all vectors operated on by the ifft() function. However, in certain instances where length(N) is not even, warnings were encountered. As a resolution, 'symmetric' was implemented. Additionally, to ensure smooth operation of ifft, an if condition was added in lines 605-608.

4. The reference used in the Full QTF implementation is
   Duarte, Tiago M., António J. Sarmento, and Jason M. Jonkman. "Effects of second-order hydrodynamic forces on floating offshore wind turbines." 32nd ASME Wind Energy Symposium. 2014.

Running RM3 example with regular waves and body(1).QTFs = 1;

Unrecognized field name "QTFs".

Error in bodyClass/QTF_excitation (line 615)
            Omega_coarse = 2 * pi ./ obj.hydroData.hydro_coeffs.excitation.QTFs.Sum(1).PER_i;

Error in bodyClass/hydroForcePre (line 379)
                obj.QTF_excitation(waveAmpTime);

Error in initializeWecSim (line 284)
        body(it).hydroForcePre(waves(1).omega,waves(1).direction,simu.cicTime,waves(1).bem.count,simu.dt,...

Error in run (line 99)
evalin('caller', strcat(script, ';'));

Error in wecSim (line 40)
run('initializeWecSim');

This occurred because I didn't run bemio with the 12d and 12s files in my hydrodata repo. I added them, then re-ran bemio. Also tried deleting h5 file and re-running w/ these files in there, and saw the same error. It doesn't appear the QTFs field is being properly created in the hydroData sub-structure. It does not appear in the *h5 file.

The error seems to stem from readWAMIT

It seems to only be filling these fields on indices greater than 1, (screen shot of hydro.QTFs after completion of readWAMIT. Of the populated structures, all the sub fields that should be there seem to be!)

The problem is that the provided files were _2.##[s,d], and I think corresponded to body(2), but I mod'd the wecSimInputFile as instructed above with i=1. So this is my own mistake but something to clarify in documentation. Another point that may be worthy of discussion: how is this physically different from the mean-drift force implementation we already have the option to account for?

I recommend adding a warning flag in the try-catch in readWAMIT beginning in line 369. Maybe something to the effect of "warning: qtf file parser is looking for a file ending in extension ##s or ##d. Please remove or rename any extraneous files with matching extensions or this may cause issues in simulation"

Otherwise, I ran this for regular and irregular waves and it seems to function as intended. Still need to check the underlying mathematics in the QTF term calculation. Right now we see a mean drift in surge of ~300 m for in regular waves and ~700 m in irregular waves...I believe this is because A > 0 and X > 0 by definition in equations 14 and 15 and in regular waves all (but one?) of the H_mu elements will be zero.

Hi @dforbush2

1. According the QTF equations implemented, as well as Newman's paper, the hydrodynamic force and moment acting upon a marine vehicle or structure in regular waves will generally include a second-order nonlinear component proportional to the square of the wave amplitude. Thus, we should a non-zero force for regular waves.
   It worth noting that HydroDyn displays an error message and interrupts the solution if the QTFs were enabled in a regular waves simulation. Which I do not think is needed for WEC-Sim.

2. I will make sure about adding this to the documentation, we can also add an error message to help users identify the problem if they encounter it. While coding it I thought its intuitive that if you didn't include the QTFs for a body in the BEMIO and tried to solve for it in WEC-Sim you will get errors.

3. When two linear wave periods are the same, the quadratic forces for the difference frequency are the same as the mean drift forces. The mean drift force has two contributions:
   a) 2nd-order hydrodynamic pressure integration due to the first-order wave
   b) Interaction between the first-order motion and the first-order wave
   Currently, WEC-Sim only reads mean drift coefficients representing the first contribution. According to WAMIT documentation, the quadratic forces may be evaluated more accurately than using the conventional pressure integration method.

4. Thanks for the warring warning flag suggestion, I will add this to the code.

Please don't forget to check the run from Simulink feature as it needs a revision.

@dforbush2 I've edited my last comment to correct the first point (1) regarding the application of the QTFs in regular waves.

Thanks for these changes @MShabara. Things appear to be working intuitively, though surge displacements in regular waves are still really high (~180 m for 100s run time with 60s wave ramp). We had discussed the validity of this approach in regular waves, and it seems like a lot of drift in surge.

Looking at irregular waves, there is some unusual behavior when QTF is paired with quadratic drag. The direction of average drift is reversed (though the magnitude of this displacement is reduced as expected) This is certainly non-physical and I don't have a good explanation of how it arises: in isolation Drag seems to work fine. In REGULAR waves, there is no reversal of drift direction in CD + QTF cases (but again, not sure if its appropriate to consider regular waves in this context)

The one point of potential confusion that I recommend is clarified in the documentation is how QTF relates to mean drift (i.e., the user shouldn't use body(i).meanDrift = 1 AND body(i).QTFs = 1.

I cannot get run from Simulink to work here however, i cannot get it to work anywhere on 2023b or 2024a... @akeeste is this something you have insight on?

@dforbush2 @MShabara Yes I can look into Run From Simulink. I'm late on that but should have time next week

@dforbush2 Thanks for your feedback.
You can get large displacements given that the QTF WAMIT files are for TLP Wind Turbine, but we are testing them on an RM3 model.

I just tested the QTFs and the Quadratic drag and I noticed this behavior too, I think that is because of the current QTF calculations does not correct for the body frame motion. I will make this on the top of my to-do list after I finishing the documentation and the example.

I also noticed that the wave waveform change when enabling the quadDrag, this affects the value of the QTFs.
body(2).QTFs = 1;
body(2).quadDrag.cd = [1.25, 1.25, 0.175, 1.25, 1.25 , 0.1];

body(2).QTFs = 1; or body(2).QTFs = 0;

@akeeste Thanks for helping with the run from WEC-Sim issue, looking forward to your comments.

@MShabara please merge https://github.com/MShabara/WEC-Sim/pull/2 and/or grant me write access to your fork. Thanks!