Impedance Matching and WECSim Validation

riasatmorshed commented 2 weeks ago

Hello!

I was trying to assess the concept of impedance matching matching for my OSWEC so that maximum available power output can be achieved. I am using the following equation to calculate the Damping and Stiffness value of a PTO for a certain wave frequency:

Kp = radiation damping Ki= (\omega^2)*(mass+added_mass)-hydrostatic_stiffness

The values for radiation damping, added mass etc. are taken, say, for a particular frequency.

My geometry is 0.8m in z axis and has 0.1 meter freeboard. I have set my center of gravity to be (0,0,-0.35m) and the rotation point are set at (0,0,-0.57m). In the attached code, I have obatained corresponding Kp and Ki value for each frequency.

Now, when I input those value into WECSim, for PTO location to be at the rotation point I mentioned, used the obtained value for stiffness and damping, and used the corresponding wave period(or, frequency), the power output does not seem to be the maximum output. Apparently, higher damping values are yielding more power (when I took the mean of "powerInternalMechanics" array of WECSim). Even though the question has overlapping issues with WECSim, I hope you guys can prove some insight into this.I am attaching my code herewith. Thanks in advance! Impedance_verification2.zip

dtgaebe commented 2 weeks ago

Hi @riasatmorshed ,

Would you be so kind and also provide the mesh WecOptTool_height_adjusted.dat ?

riasatmorshed commented 2 weeks ago

Oh sorry! My bad that I forgot to upload the accessory files. Here is the .dat file attached. Please let me know if you need anything else. WecOptTool_height_adjusted.zip

dtgaebe commented 2 weeks ago

Hi @riasatmorshed ,

Your code looks good! There are a few things that come to my mind that you should double check

Do you use the same BEM coefficients and degrees of freedom in WEC-Sim?
- The optimal damping (for mechanical power the way you calculate it) is immediately a function the radiation damping. If you would use different BEM results this will change the power output.
- Also, just to make sure, limit your WEC-Sim simulation to pitch only
- Is there any other friction $B_f$ considered in WEC-Sim? This will also change the optimal damping
Did you double check the WEC-Sim sign conventions?

I will re-derive your controller so that other users can follow and you can check that the WEC-Sim PTO actually applies the forces in the assumed directions:

Z_iU = F_{ex} - F_p

with

\begin{align}
Z_i &= B(\omega) + B_f +  i\bigg (\omega(m+A(\omega)) -\frac{K_{HS}}{\omega} \bigg)\\
& = R_i +i X_i
\end{align}

say you want a PI controller

\begin{align}
F_p &= K_p U + K_iX \\
      &= K_p U +K_i \frac{U}{i\omega} \\
      &= C_{PI}U,\\
\textrm{with } C_{PI} &= K_p + \frac{K_i}{i\omega}
\end{align}

For the dynamics defined so that the PTO force opposes the excitation force, the optimal control force for mechanical power is (as you know)

F_p = Z_i^*U \\

or,

\begin{align}
Z_i^* &= C_{PI} \\
R_i -i X_i &= K_p + \frac{K_i}{i\omega} \\
R_i +i (- X_i) &= K_p + i(\frac{-K_i}{\omega}) \\
\end{align}

thus,

\begin{align}
K_p &= R_i = \Re{(Z_i^*)} = B(\omega) + B_f \\
K_i &= \omega X_i = -\omega \Im{(Z^*_i)} = \omega \Im{(Z_i)} = \omega \bigg( \omega(m+A(\omega)) -\frac{K_{HS}}{\omega}  \bigg)
\end{align}

@jtgrasb offered help with your WEC-Sim model if you need it. You should use the WEC-Sim github issues if you do. Feel free to reference this issue.

riasatmorshed commented 2 weeks ago

@dtgaebe,

Now that you say it, I used NEMOH result for my WECSim run and I noticed that the hydrodynamic result from WecOptTool and NEMOH are not matching. They are off by around a multiple of 4 to 5. But my question is, even if the frequency range is different for both the simulation, for example, for NEMOH, I ran the simulation for 0.04 rad/s to 10 rad/s and for WecOptTool, I ran the simulation for 0.188 rad/s to 3.78 rad/s. Even if that, the value for a particular frequency should match, right?

Just for reference, the rotation center for both the WecOptTool and NEMOH atr -0.57m and CG used in Bemio function to convert NEMOH data to WECSim is also (0,0,-0,35). I can't see any discrepancies between them- sure, I might be wrong!

I restricted the WECSim to pitch only, meaning, my Simulink file only contains a rotational PTO for pitch.
However, for @jtgrasb's reference I am adding my WECSim files here. WECSim_Files.zip

riasatmorshed commented 6 days ago

Hello!

I would just like to follow up on the issue since I have a partial fix to the problem. Does it have to do anything with the unit? for example, in the WECSim, the damping and stiffness value should be in Nsm/rad, Nm/rad respectively. But these are the rotational damping and stiffness and am I getting the damping and stiffness from the WecOptTool as Ns/m, and N/m respectively?

I am not sure. Just wondering if that might be the case! If not, I'd really appreciate any explanations regarding that. Thanks!

dtgaebe commented 6 days ago

Hi @riasatmorshed ,

They are off by around a multiple of 4 to 5 This is odd and can't originate from WEC-Sim presenting the BEM results in dimensionless form (https://wec-sim.github.io/WEC-Sim/main/theory/theory.html#boundary-element-method-bem)

But my question is, even if the frequency range is different for both the simulation, for example, for NEMOH, I ran the simulation for 0.04 rad/s to 10 rad/s and for WecOptTool, I ran the simulation for 0.188 rad/s to 3.78 rad/s. Even if that, the value for a particular frequency should match, right?

Yes, the range doesn't matter and they should match (within reason) at a particular frequency.

Just for reference, the rotation center for both the WecOptTool and NEMOH atr -0.57m and CG used in Bemio function to convert NEMOH data to WECSim is also (0,0,-0,35)

For clarification, currently you're not yet using the core of WecOptTool (e.g. wec.solve()). You are using the WecOptTool functions that help you run the BEM solver capytaine https://github.com/capytaine/capytaine (this avoids needing additional functions like BEMIO to convert the BEM results into a form that can readily be used to define a WecOptTool wec object and you use linear frequency domain analysis for your controller gains (instead of wec.solve())

am I getting the damping and stiffness from the WecOptTool as Ns/m, and N/m respectively?

No, your Kp and Ki are in rotational units since your intrinsic impedance is in rotational units.

since I have a partial fix to the problem

What is your partial fix?

riasatmorshed commented 6 days ago

Hi @dtgaebe!

Thanks for your clarification on the WecOptTool using capytaine engine and the rotational units also. Pardon my wording on "partial fix", I actually meant "possible" fix- autocorrect!

However, I have run the hydrodynamic case using capytaine and produced .h5 using WECSim's BEMIO function. But the pto impedance concept still doesn't seem to add up. I am using "mean(output.ptos.powerInternalMechanics(:,5))" command to assess power output.

I am attaching the .nc file, .h5 file, and the corresponding WECSim module. I really appreciate your help and patience! Thanks!

capytaine_run.zip

jtgrasb commented 5 days ago

Hi @riasatmorshed. I just took a look at both your WecOptTool model and the WEC-Sim model. It is difficult to tell exactly where the mismatch is coming from without seeing the script you are using for your WEC-Sim Capytaine run, but I believe that may be the cause of your issues. Unfortunately, it is difficult to compare the WEC-Sim Capytaine data to the WecOptTool Capytaine data because they should be different for a body rotating about a hinge. The Capytaine run for WEC-Sim should be run with the center of rotation equal to the center of gravity. WEC-Sim (in Simulink) then converts the hydrodynamics to the point about which the rotation is specified (in this case, the PTO). I would recommend checking out the oswec.py script in WEC-Sim/examples/BEMIO and call_capytaine.py as an example of how to set up the Capytaine run for WEC-Sim (you should just need one body for your run). You could also try to replicate your WecOptTool Capytaine run to make sure everything matches, then only change the center of rotation to make sure everything else matches.

riasatmorshed commented 5 days ago

Thanks a lot for taking the trouble of going through my files, @jtgrasb! I will surely compare my WecOptTool and Capytaine result. #

The Capytaine run for WEC-Sim should be run with the center of rotation equal to the center of gravity.

I think, that's where the problem is. Because i set my center of rotation to be (0,0,-0.57) in Capytaine (which is surely not the CG) and used the same coordinates for the PTO rotation point in WECsim.

It is difficult to tell exactly where the mismatch is coming from without seeing the script you are using for your WEC-Sim Capytaine run, but I believe that may be the cause of your issues.

Just letting you know that, my capytaine script is in the "capytaine_run.zip" file that is uploaded in this thread. I basically modified the script found on WECSim's "bemio" folder for my capytaine run. But meanwhile I will also cross-check everything and let you know if I found any discrepancy. Thanks for the support again!

jtgrasb commented 5 days ago

Ahh, thanks for pointing out the 'capytaine_run.zip' file. I missed that one. Hopefully changing the rotation center back to the CG for WEC-Sim works.

riasatmorshed commented 4 days ago

@jtgrasb,

Thanks for pointing out the major problem. It "sort of" fixes the issue. Meaning, I used Capytaine throughout the whole process- from generating hydrodynamic data to calculate impedance value. And surprisingly the damping and stiffness value of the PTO seem to make sense. However, it is not exactly accurate. By which I mean, for wave period 1.967s, I determined the damping and stiffness value should be 1976.662 Nsm/rad, and 11345.447Nm/rad respectively. By keeping the stiffness value that, I was playing around with the damping value to see what is the power output. I noticed that beyond the damping value of 1976.662 Nsm/rad, say for 2500 Nsm/rad, the power increases and the power output seem to decrease at around 3500 Nsm/rad. I mean, the damping and stiffness value obtained are way different than the WecOptTool output and for this current run at least I see a decreasing trend of power output beyond a certain value. That's why I said "sort of" fixes the issue. Although, if I run a parameter search using WECSimMCR, that might be of help to some extent in this regard I guess, however, the theoretical PTO damping and stiffness value should check out, right? What do you think?

Just a follow-up question on that. Isn't WecOptTool using Capytaine engine to solve hydrodynamic coefficients, right? I ran both of them for the same frequency range and all, I noticed that as the frequency increases, the discrepancy between, say, Radiation Damping Coefficients value, increases. Is it supposed to happen? If so, then why?

Thanks a lot again! I really appreciate it.

Here's everything- Capytaine script to WECSim project file in this .zip file. fullCapytaineRunForImpendance.zip

jtgrasb commented 3 days ago

Hi @riasatmorshed, I just opened your files. I apologize, I don't think I explained it well enough, but you do actually need different Capytaine runs for WEC-Sim and the impedance calculation (the only difference being the rotation center). The impedance model should have the rotation center specified so that it calculates the hydrodynamic coefficients with respect to the rotation center. The BEM run for WEC-Sim, on the other hand, should have rotation center = center of gravity as you implemented in 'oswec_impedance_matching.ipynb', because WEC-Sim converts the hydrodynamic forces to be around the hinge point afterwards. Hopefully this helps resolve the issues.

And yes, WecOptTool is using Capytaine to solve the hydrodynamic coefficients. If you call run_bem from WecOptTool for the same frequencies as calling Capytaine directly, the results should be the exact same (assuming every other factor is the same). A couple things to check that may cause slight differences are the density and water depth.

sandialabs / WecOptTool

Impedance Matching and WECSim Validation #355