Floating OWC example using the two rigid body approach v2

OWC/FloatingOWC_W2W/Mooring/lines.txt

@@ -0,0 +1,55 @@
+Mooring line data file for MoorDyn in libmoordyn.dll
+---------------------- LINE TYPES -----------------------------------------------------
+LineType  Diam      Mass/m        EA      BA/-zeta   EI         Cd     Ca    CdAx     CaAx    
+(-)       (m)       (kg/m)        (N)     (N-s/-)  (N-m^2)      (-)     (-)    (-)     (-) 
+Chain    0.032      34.82       1.5e6       -0.8     1e4        1.6     1.0    0.5     0.5
+---------------------------- BODIES -----------------------------------------------------
+ID   Attachment  X0     Y0    Z0     r0     p0     y0      Mass  CG*       I*      Volume     CdA*   Ca
+(#)     (-)      (m)    (m)   (m)   (deg)  (deg)  (deg)    (kg)  (m)    (kg-m^2)    (m^3)    (m^2)  (-)
+1     Coupled     0      0     -2.58      0      0      0    2.914e6 -31.96    1.53e9  2.9013e3    0     0
+---------------------- POINTS  -----------------------------------------------------
+ID      Attachment    X      Y       Z            Mass   Volume     CdA     CA
+(#)     (word/ID)    (m)    (m)     (m)           (kg)   (mˆ3)      (m^2)   (-)
+1         Body1     -9.28	 0	    -2.58        0         0         0      0
+2         Fixed     -211.2	 0	    -80.0        0         0         0      0
+3         Body1     4.64     8.04   -2.58        0         0         0      0
+4         Fixed     105.6    182.9   -80.0       0         0         0      0
+5         Body1     4.64    -8.04    -2.58       0         0         0      0
+6         Fixed     105.6   -182.9   -80.00      0         0         0      0
+7         Free     -53	     0	     -25.00     4030.46     0.0305    0      1
+8         Free      26.5     45.9    -25.00     4030.46     0.0305    0      1
+9         Free      27.5    -45.9    -25.00     4030.46     0.0305    0      1
+10        Free      -85	     0	     -15.00     36139.83    0.2241    0      1
+11        Free      42.5     73.61	 -15.00     36139.83    0.2241    0      1
+12        Free      42.5    -73.61   -15.00     36139.83    0.2241    0      1
+---------------------- LINES -----------------------------------------------------
+ID     LineType  AttachA  AttachB  UnstrLen     NumSegs   LineOutputs
+(#)      (-)       (#)      (#)       (m)         (-)       (-)
+1        Chain     2         10        143.28      15        tp
+2        Chain     4         11        143.28      15        tp
+3        Chain     6         12        143.28      15        tp
+4        Chain     10        7         37.01        5        tp
+5        Chain     11        8         37.01        5        tp
+6        Chain     12        9         37.01        5        tp
+7        Chain     7         1         50.4         8        tp
+8        Chain     8         3         50.4         8        tp
+9        Chain     9         5         50.4         8        tp
+---------------------- SOLVER OPTIONS-----------------------------------------
+0.0005   dtM          - time step to use in mooring integration
+0        WaveKin      - wave kinematics flag (0=neglect, the only option currently supported)
+3.0e6    kBot         - bottom stiffness
+3.0e5    cBot         - bottom damping
+80       WtrDpth      - water depth
+5.0      CdScaleIC    - factor by which to scale drag coefficients during dynamic relaxation IC gen
+0.001    threshIC     - threshold for IC con
+300.0    TmaxIC       - max time for ic gen (s)
+2        writeLog     - Write a log file
+0        WriteUnits   - 0: do not write the units header on the output files
+-------------------------- OUTPUTS --------------------------------
+FairTen1
+FairTen2
+FairTen3
+FairTen3
+FairTen5
+FairTen6
+--------------------- need this line ------------------

OWC/FloatingOWC_W2W/README.md

@@ -0,0 +1,31 @@
+# Oscillating Water Column
+
+**Authors:** Mohamed Shabara, Jeff Grasberger and Jorge Leon-Quiroga
+
+**Version:** WEC-Sim v6.1.x
+
+**Geometry:** Floating OWC Rigid Body Approach
+
+This model simulates a Floating Oscillating Water Column (OWC) device using a rigid-body approach. 
+It incorporates performance curves for both a Wells Turbine and a generator, providing a realistic 
+representation of their dynamic behavior. Additionally, the model accounts for the presence of a
+ mooring system to ensure the stability of the floating body.
+
+An optimal control strategy is implemented to maximize turbine efficiency by dynamically adjusting system parameters based on operating speed.
+
+This model is based on the research detailed in:
+Shabara, Mohamed A., et al. "Optimal Control of Floating Oscillating Water Column Wave Energy Converters." 2025 American Control Conference (ACC), IEEE, 2025.
+
+Citation:
+If this model is used in your work, please include the following citation:
+@inproceedings{shabara2025OWC,
+  title={Optimal Control of Floating Oscillating Water Column Wave Energy Converters},
+  author={Shabara, Mohamed A et al.},
+  booktitle={2025 American Control Conference (ACC)},
+  year={2025},
+  organization={IEEE}
+}
+
+**Acknowledgment:** This material is based upon work supported by the U.S.
+Department of Energy’s Office of Energy Efficiency and Renewable Energy 
+under the Water Power Technologies Office Award Number DE-EE0008895

OWC/FloatingOWC_W2W/calcTurbPower.m

@@ -0,0 +1,22 @@
+function pGen = calcTurbPower(omega, controlTorque, PgenRated, TgenMax)
+
+pGen = zeros(1);
+
+% Absolute values were added because the generator only spins in one direction
+omega = abs(omega);
+T_ctrl = abs(controlTorque);
+
+term1 = 6.280e-10*omega^5 - 1.003e-6*omega^4 + 4.416e-4*omega^3 - 0.07028*omega^2 + 2.106*omega;
+term2 = -0.06941*T_ctrl^2 - 6.014*T_ctrl - 4.664e-8*omega^4*T_ctrl + 1.79e-8*omega^3*T_ctrl^2;
+term3 = + 2.089e-5*omega^3*T_ctrl - 9.816e-6*omega^2*T_ctrl^2 - 0.003354*omega^2*T_ctrl;
+term4 = + 0.001292*omega*T_ctrl^2 + 1.189*omega*T_ctrl - 139.3;
+
+tmp = term1 + term2 + term3 + term4;
+
+% Power can not be more than the maximum allowable
+tmp = min([tmp, TgenMax * omega, ones(size(tmp))*PgenRated]);
+
+% We can't have negative power
+pGen = max(tmp , 0);
+
+end

OWC/FloatingOWC_W2W/fitCurves.m

@@ -0,0 +1,18 @@
+function [phiInt, piInt, etaInt] = fitCurves(psiInstantaneous, psiVec, etaTurbVec, piVarVec, phiVec)
+    % psiInstantaneous = clip(psiInstantaneous, min(psi), max(psi));
+    if psiInstantaneous > max(psiVec)
+        psiInstantaneous = max(psiVec);
+    elseif psiInstantaneous < min(psiVec)
+        psiInstantaneous = min(psiVec);
+    end
+
+    % Interpolate for the phi Value
+    phiInt = interp1(psiVec, phiVec, psiInstantaneous);
+
+    % Interpolate for the Pi Value
+    piInt = interp1(psiVec, piVarVec, psiInstantaneous);
+
+    % Interpolate for the eta Value
+    etaInt = interp1(psiVec, etaTurbVec, psiInstantaneous);
+
+end

OWC/FloatingOWC_W2W/fittingFunctions.m

@@ -0,0 +1,103 @@
+function  [psi, etaTurb, phi, piVar] = fittingFunctions()
+    % Curve splitting
+    psi1 = [0.008, 0.009];
+    mu_turb1 = [0, 0.100];
+    m1 = (mu_turb1(2) - mu_turb1(1)) / (psi1(2) - psi1(1));
+    q1 = mu_turb1(1) - m1 * psi1(1);
+
+    psi2 = [0.009, 0.0128];
+    mu_turb2 = [0.100, 0.200];
+    m2 = (mu_turb2(2) - mu_turb2(1)) / (psi2(2) - psi2(1));
+    q2 = mu_turb2(1) - m2 * psi2(1);
+
+    psi3 = [0.0128, 0.0199];
+    mu_turb3 = [0.200, 0.300];
+    m3 = (mu_turb3(2) - mu_turb3(1)) / (psi3(2) - psi3(1));
+    q3 = mu_turb3(1) - m3 * psi3(1);
+
+    psi4 = [0.0199, 0.0404];
+    mu_turb4 = [0.300, 0.500];
+    m4 = (mu_turb4(2) - mu_turb4(1)) / (psi4(2) - psi4(1));
+    q4 = mu_turb4(1) - m4 * psi4(1);
+
+    psi5 = [0.0404, 0.0634, 0.0878];
+    mu_turb5 = [0.500, 0.595, 0.500];
+    pp = polyfit(psi5, mu_turb5, 2);
+
+    psi6 = [0.0878, 0.1141];
+    mu_turb6 = [0.500, 0.300];
+    m6 = (mu_turb6(2) - mu_turb6(1)) / (psi6(2) - psi6(1));
+    q6 = mu_turb6(1) - m6 * psi6(1);
+
+    psi7 = [0.1141, 0.1192];
+    mu_turb7 = [0.300, 0.273];
+    m7 = (mu_turb7(2) - mu_turb7(1)) / (psi7(2) - psi7(1));
+    q7 = mu_turb7(1) - m7 * psi7(1);
+
+    psi9 = [0.2122, 0.3];
+    mu_turb9 = [0.0016, 0.0];
+    m9 = (mu_turb9(2) - mu_turb9(1)) / (psi9(2) - psi9(1));
+    q9 = mu_turb9(1) - m9 * psi9(1);
+
+    psi10 = [0.1, 0.114];
+    Pi_turb = [0.3156, 0.3170];
+    m10 = (Pi_turb(2) - Pi_turb(1)) / (psi10(2) - psi10(1));
+    q10 = Pi_turb(1) - m10*psi10(1);
+
+    m = [m1, m2, m3, m4, 0, m6, m7, 0, m9, m10];
+    q = [q1, q2, q3, q4, 0, q6, q7, 0, q9, q10];
+
+    % Compute the dimensionless functions
+    psi = linspace(-0.3, 0.3, 5000);
+    etaTurb = zeros('like',psi);
+    phi = zeros('like',psi);
+    piVar = zeros('like',psi);
+
+    for ii = 1:length(psi)
+        PSI = abs(psi(ii));
+
+        if PSI >= .008 && PSI < 0.009
+            ETA_turb = m(1) * PSI + q(1);
+        elseif PSI  >= .009 && PSI < 0.0128
+            ETA_turb = m(2) * PSI + q(2);
+        elseif PSI  >= .0128 && PSI < 0.0199
+            ETA_turb = m(3) * PSI + q(3);
+        elseif PSI >= .0199 && PSI <= 0.0404
+            ETA_turb = m(4) * PSI + q(4);
+        elseif PSI > .0404 && PSI < 0.0878
+            ETA_turb = pp(1) * PSI^2 + pp(2) * PSI + pp(3);
+        elseif PSI >= .0878 && PSI < 0.1141
+            ETA_turb = m(6) * PSI + q(6);
+        elseif PSI >= .1141 && PSI <= 0.1192
+            ETA_turb = m(7) * PSI + q(7);
+        elseif PSI >= .1192 && PSI <= 0.2122
+            ETA_turb = 200.6*exp(PSI*-55.35);
+        elseif PSI > .2122 && PSI <= 0.3
+            ETA_turb = m(9) * PSI + q(9);
+        else
+            ETA_turb = 0;
+        end
+
+        etaTurb(ii) = ETA_turb;
+
+    end
+
+    for ii = 1:length(psi)
+        PSI = abs(psi(ii));
+
+        if PSI >= 0 && PSI <= 0.1
+            PHI = 0.7750*PSI;
+        else
+            corr = 0.00167;
+            PHI = -2.503*PSI^2 + 1.608*PSI - 0.05664 - corr;
+        end
+
+        phi(ii) = PHI;
+
+        if PSI > .1 && PSI <= 0.114
+            piVar(ii) = (m(10) * PSI + q(10)) / 100;
+        else
+            piVar(ii) = etaTurb(ii) * PSI * PHI;
+        end
+    end
+end

OWC/FloatingOWC_W2W/genPower.m

@@ -0,0 +1,12 @@
+function P_gen = genPower(omega, T_ctrl, TgenMax, PgenTated)
+
+term1 = 6.280e-10 * omega^5 - 1.003e-6 * omega^4 + 4.416e-4 * omega^3 - 0.07028 * omega^2 + 2.106 * omega;
+term2 = -0.06941 * T_ctrl^2 - 6.014 * T_ctrl - 4.664e-8 * omega^4 * T_ctrl + 1.79e-8 * omega^3 * T_ctrl^2;
+term3 = + 2.089e-5 * omega^3 * T_ctrl - 9.816e-6 * omega^2 * T_ctrl^2 - 0.003354 * omega^2 * T_ctrl;
+term4 = + 0.001292 * omega * T_ctrl^2 + 1.189 * omega * T_ctrl - 139.3;
+
+tmp = term1 + term2 + term3 + term4;
+
+P_gen = min(tmp, TgenMax * omega);
+
+end