Merge pull request #68 from WISDEM/rc1.3.2

RAFT v1.3.2

Author: lucas-carmo

.github/workflows/CI_RAFT.yml

@@ -18,6 +18,7 @@ jobs:
                 os: ["ubuntu-latest", "macOS-latest", "windows-latest"]
                 python-version: ["3.10", "3.11"]
 
+
         steps:
             - name: checkout repository
               uses: actions/checkout@v4

examples/VolturnUS-S_example.yaml

@@ -27,10 +27,10 @@ cases:
 
 turbine:
 
-    mRNA          :     991000        #  [kg]       RNA mass 
+    mRNA          :     1001505.907   #  [kg]       RNA mass: m_blades = 3*68076.969 = 204230.907 kg (from ED.sum), m_yawbearing = 100000 kg, m_hub = 190000 kg, m_nacelle = 507275 kg
     IxRNA         :          0        #  [kg-m2]    RNA moment of inertia about local x axis (assumed to be identical to rotor axis for now, as approx) [kg-m^2]
     IrRNA         :          0        #  [kg-m2]    RNA moment of inertia about local y or z axes [kg-m^2]
-    xCG_RNA       :          0        #  [m]        x location of RNA center of mass [m] (Actual is ~= -0.27 m)
+    xCG_RNA       :         -7.27     #  [m]        x location of RNA center of mass [m]: xg_blades = -13.0 m, xg_hub = -11.0 m, xg_nac = -4.99 m, xg_yaw = 0.0 m
     hHub          :     150.0       #  [m]        hub height above water line [m]
     #rRNA          :  [0,0, 148.742] #  [m] Can use the position of the RNA reference point (which the RNA yaws about) with respect to the FOWT reference point. If also providing hHub, the z coord be overwritten.
     Fthrust       :       1500.0E3    #  [N]        temporary thrust force to use
@@ -43,7 +43,7 @@ turbine:
     precone     : 4.0        # [deg]
     shaft_tilt  : 6.0        # [deg]
     overhang    : -12.0313    # [m]
-    aeroMod     : 1          # 0 aerodynamics off; 1 aerodynamics on
+    aeroServoMod     : 1          # 0 aerodynamics off; 1 aerodynamics on
 
 
     blade: 

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "OpenRAFT"
-version = "1.3.1"
+version = "1.3.2"
 description = "RAFT: Response Amplitudes of Floating Turbines"
 readme = "README.md"
 requires-python = ">=3.9"

raft/omdao_raft.py

@@ -600,7 +600,7 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
                 ring_spacing = inputs[m_name+'ring_spacing']
                 n_stiff = 0 if ring_spacing == 0.0 else int(np.floor(s_height / ring_spacing))
                 s_ring = (np.arange(1, n_stiff + 0.1) - 0.5) * (ring_spacing / s_height)
-                if s_ring:
+                if np.any(s_ring):
                     if not m_shape == 'rect':
                         d_ring = np.interp(s_ring, s_grid, design['platform']['members'][i]['d'])
                     else:
@@ -623,7 +623,7 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
                     t_cap = t_cap[:-1]
                     di_cap = di_cap[:-1]
                 # Combine with ring stiffeners
-                if s_ring:
+                if np.any(s_ring):
                     s_cap = np.r_[s_ring, s_cap]
                     t_cap = np.r_[inputs[m_name+'ring_t']*np.ones(n_stiff), t_cap]
                     di_cap = np.r_[d_ring-2*inputs[m_name+'ring_h'], di_cap]

raft/raft_fowt.py

@@ -3,7 +3,7 @@
 import os
 import matplotlib.pyplot as plt
 import numpy as np
-from scipy.interpolate import interp1d, interp2d, griddata
+from scipy.interpolate import interp1d, RegularGridInterpolator, griddata
 
 import raft.member2pnl as pnl
 from raft.helpers import *
@@ -1575,14 +1575,13 @@ def calcQTF_slenderBody(self, waveHeadInd, Xi0=None, verbose=False, iCase=None,
                         f_rslb  += - rho*v_i * aux
 
 
-                        # ----- axial/end effects  ------
+                        # ----- end effects  ------                               
                         # note : v_i and a_i work out to zero for non-tapered sections or non-end sections
+                        # TODO: Should probably skip this part if along the length of the element, i.e. if a_i and v_i are different from zero due to tapering
                         if circ:
                             v_i = np.pi/12.0 * abs((mem.ds[il]+mem.drs[il])**3 - (mem.ds[il]-mem.drs[il])**3)  # volume assigned to this end surface
-                            a_i = np.pi*mem.ds[il] * mem.drs[il]   # signed end area (positive facing down) = mean diameter of strip * radius change of strip
                         else:
                             v_i = np.pi/12.0 * ((np.mean(mem.ds[il]+mem.drs[il]))**3 - (np.mean(mem.ds[il]-mem.drs[il]))**3)    # so far just using sphere eqn and taking mean of side lengths as d
-                            a_i = (mem.ds[il,0]+mem.drs[il,0])*(mem.ds[il,1]+mem.drs[il,1]) - (mem.ds[il,0]-mem.drs[il,0])*(mem.ds[il,1]-mem.drs[il,1])
 
                         f_2ndPot += mem.a_i[il]*p_2nd*mem.q # 2nd order pressure
                         f_2ndPot += rho*v_i*Ca_End*np.matmul(mem.qMat, acc_2ndPot) # 2nd order axial acceleration
@@ -1593,6 +1592,11 @@ def calcQTF_slenderBody(self, waveHeadInd, Xi0=None, verbose=False, iCase=None,
                         p_drop    = -2*0.25*0.5*rho*np.dot(np.matmul(mem.p1Mat + mem.p2Mat, u[:,i1, il]-nodeV[:, i1, il]), np.conj(np.matmul(Ca_p1*mem.p1Mat + Ca_p2*mem.p2Mat, u[:,i2, il]-nodeV[:, i2, il])))
                         f_conv   += mem.a_i[il]*p_drop*mem.q
 
+                        u1_aux    = np.matmul(Ca_p1*mem.p1Mat + Ca_p2*mem.p2Mat, u1_aux)
+                        u2_aux    = np.matmul(Ca_p1*mem.p1Mat + Ca_p2*mem.p2Mat, u2_aux)
+                        f_transv  = 0.25*mem.a_i[il]*rho*(np.conj(u1_aux)*nodeV_axial_rel[i2, il] + u2_aux*np.conj(nodeV_axial_rel[i1, il]))
+                        f_conv   += f_transv
+
                         F_2ndPot += translateForce3to6DOF(f_2ndPot, mem.r[il,:])
                         F_conv   += translateForce3to6DOF(f_conv, mem.r[il,:])                   
                         F_axdv   += translateForce3to6DOF(f_axdv, mem.r[il,:])
@@ -1788,10 +1792,15 @@ def calcHydroForce_2ndOrd(self, beta, S0, iCase=None, iWT=None, interpMode='qtf'
         else:
             f = np.zeros([self.nDOF, self.nw]) # Force amplitude
             for idof in range(0,self.nDOF):
-                # Interpolate the QTF matrix to the same frequencies as the wave spectrum. Need to interpolate real and imaginary part separately.
-                qtf_interp_Re = interp2d(self.w1_2nd, self.w1_2nd, qtf_interpBeta[:,:, idof].real, bounds_error=False, fill_value=(0))(self.w, self.w)
-                qtf_interp_Im = interp2d(self.w1_2nd, self.w1_2nd, qtf_interpBeta[:,:, idof].imag, bounds_error=False, fill_value=(0))(self.w, self.w)
-                qtf_interp = qtf_interp_Re + 1j*qtf_interp_Im
+                qtf_interp_Re_interpolator = RegularGridInterpolator((self.w1_2nd, self.w1_2nd), qtf_interpBeta[:, :, idof].real, bounds_error=False, fill_value=0)
+                qtf_interp_Im_interpolator = RegularGridInterpolator((self.w1_2nd, self.w1_2nd), qtf_interpBeta[:, :, idof].imag, bounds_error=False, fill_value=0)
+
+                w_mesh = np.meshgrid(self.w, self.w, indexing='ij')
+                points = np.array([w_mesh[0].ravel(), w_mesh[1].ravel()]).T
+
+                qtf_interp_Re = qtf_interp_Re_interpolator(points).reshape(len(self.w), len(self.w))
+                qtf_interp_Im = qtf_interp_Im_interpolator(points).reshape(len(self.w), len(self.w))
+                qtf_interp = qtf_interp_Re + 1j * qtf_interp_Im
 
                 for imu in range(1, self.nw): # Loop the difference frequencies
                     Saux = np.zeros(self.nw)

raft/raft_member.py

@@ -150,7 +150,11 @@ def __init__(self, mi, nw, BEM=[], heading=0):
             self.cap_stations = []
         else:
             self.cap_t        = getFromDict(mi, 'cap_t'   , shape=cap_stations.shape[0])   # thicknesses [m]
-            self.cap_d_in     = getFromDict(mi, 'cap_d_in', shape=cap_stations.shape[0])   # inner diameter (if it isn't a solid plate) [m]
+            if self.shape == 'circular':
+                self.cap_d_in = getFromDict(mi, 'cap_d_in', shape=cap_stations.shape[0])   # inner diameter (if it isn't a solid plate) [m]
+            elif self.shape == 'rectangular':
+                self.cap_d_in = getFromDict(mi, 'cap_d_in', shape=[cap_stations.shape[0],2])   # inner diameter (if it isn't a solid plate) [m]
+
             self.cap_stations = (cap_stations - st[0])/(st[-1] - st[0])*self.l             # calculate station positions along the member axis from 0 to l [m]
 
 
@@ -449,7 +453,7 @@ def RectangularFrustumMOI(La, Wa, Lb, Wb, H, p):
                     v_shell = V_outer-V_inner               # volume of hollow frustum with shell thickness [m^3]
                     m_shell = v_shell*rho_shell             # mass of hollow frustum [kg]
 
-                    hc_shell = ((hco*V_outer)-(hci*V_inner))/(V_outer-V_inner)  # center of volume of hollow frustum with shell thickness [m]
+                    hc_shell = ((hco*V_outer)-(hci*V_inner))/(V_outer-V_inner) if V_outer-V_inner!=0 else 0  # center of volume of hollow frustum with shell thickness [m]
 
                     dBi_fill = (dBi-dAi)*(l_fill/l) + dAi   # interpolated inner diameter of frustum that ballast is filled to [m] 
                     v_fill, hc_fill = FrustumVCV(dAi, dBi_fill, l_fill)         # volume and center of volume of solid inner frustum that ballast occupies [m^3] [m]
@@ -459,7 +463,7 @@ def RectangularFrustumMOI(La, Wa, Lb, Wb, H, p):
                     # then the ballast sits on top of the end cap. Depending on the thickness of the end cap, this can affect m_fill, hc_fill, and MoI_fill >>>>>
 
                     mass = m_shell + m_fill                 # total mass of the submember [kg]
-                    hc = ((hc_fill*m_fill) + (hc_shell*m_shell))/mass       # total center of mass of the submember from the submember's rA location [m]
+                    hc = ((hc_fill*m_fill) + (hc_shell*m_shell))/mass if mass!=0 else 0      # total center of mass of the submember from the submember's rA location [m]
 
 
                     # MOMENT OF INERTIA
@@ -492,14 +496,14 @@ def RectangularFrustumMOI(La, Wa, Lb, Wb, H, p):
                     v_shell = V_outer-V_inner                   # volume of hollow frustum with shell thickness [m^3]
                     m_shell = v_shell*rho_shell                 # mass of hollow frustum [kg]
 
-                    hc_shell = ((hco*V_outer)-(hci*V_inner))/(V_outer-V_inner)  # center of volume of the hollow frustum with shell thickness [m]
+                    hc_shell = ((hco*V_outer)-(hci*V_inner))/(V_outer-V_inner) if V_outer-V_inner!=0 else 0  # center of volume of the hollow frustum with shell thickness [m]
 
                     slBi_fill = (slBi-slAi)*(l_fill/l) + slAi   # interpolated side lengths of frustum that ballast is filled to [m]
                     v_fill, hc_fill = FrustumVCV(slAi, slBi_fill, l_fill)   # volume and center of volume of inner frustum that ballast occupies [m^3]
                     m_fill = v_fill*rho_fill                    # mass of ballast in the submember [kg]
 
                     mass = m_shell + m_fill                     # total mass of the submember [kg]
-                    hc = ((hc_fill*m_fill) + (hc_shell*m_shell))/mass       # total center of mass of the submember from the submember's rA location [m]
+                    hc = ((hc_fill*m_fill) + (hc_shell*m_shell))/mass if mass !=0 else 0     # total center of mass of the submember from the submember's rA location [m]
 
 
                     # MOMENT OF INERTIA
@@ -602,7 +606,7 @@ def RectangularFrustumMOI(La, Wa, Lb, Wb, H, p):
                 V_inner, hci = FrustumVCV(dAi, dBi, h)
                 v_cap = V_outer-V_inner
                 m_cap = v_cap*rho_cap    # assume it's made out of the same material as the shell for now (can add in cap density input later if needed)
-                hc_cap = ((hco*V_outer)-(hci*V_inner))/(V_outer-V_inner)
+                hc_cap = ((hco*V_outer)-(hci*V_inner))/(V_outer-V_inner) if V_outer-V_inner!=0 else 0
 
                 I_rad_end_outer, I_ax_outer = FrustumMOI(dA, dB, h, rho_cap)
                 I_rad_end_inner, I_ax_inner = FrustumMOI(dAi, dBi, h, rho_cap)
@@ -622,11 +626,12 @@ def RectangularFrustumMOI(La, Wa, Lb, Wb, H, p):
 
                 if L==self.stations[0]:  # if the cap is on the bottom end of the member
                     slA = sl[0,:]
-                    slB = np.interp(L+h, self.stations, sl)
+                    slB = np.zeros(slA.shape)
+                    slB = np.array([np.interp(L+h, self.stations, sl[:,0]), np.interp(L+h, self.stations, sl[:,1])])
                     slAi = sl_hole
                     slBi = slB*(slAi/slA)       # keep the same proportion in d_hole from bottom to top
                 elif L==self.stations[-1]:    # if the cap is on the top end of the member
-                    slA = np.interp(L-h, self.stations, sl)
+                    slA = np.array([np.interp(L-h, self.stations, sl[:,0]), np.interp(L-h, self.stations, sl[:,1])])
                     slB = sl[-1,:]
                     slAi = slA*(slBi/slB)
                     slBi = sl_hole
@@ -659,10 +664,10 @@ def RectangularFrustumMOI(La, Wa, Lb, Wb, H, p):
                 V_inner, hci = FrustumVCV(slAi, slBi, h)
                 v_cap = V_outer-V_inner
                 m_cap = v_cap*rho_cap    # assume it's made out of the same material as the shell for now (can add in cap density input later if needed)
-                hc_cap = ((hco*V_outer)-(hci*V_inner))/(V_outer-V_inner)
+                hc_cap = ((hco*V_outer)-(hci*V_inner))/(V_outer-V_inner) if V_outer-V_inner!=0 else 0
 
-                Ixx_end_outer, Iyy_end_outer, Izz_end_outer = RectangularFrustumMOI(slA, slB, h, rho_cap)
-                Ixx_end_inner, Iyy_end_inner, Izz_end_inner = RectangularFrustumMOI(slAi, slBi, h, rho_cap)
+                Ixx_end_outer, Iyy_end_outer, Izz_end_outer = RectangularFrustumMOI(slA[0], slA[1], slB[0], slB[1], h, rho_cap)
+                Ixx_end_inner, Iyy_end_inner, Izz_end_inner = RectangularFrustumMOI(slAi[0], slAi[1], slBi[0], slBi[1], h, rho_cap)
                 Ixx_end = Ixx_end_outer-Ixx_end_inner
                 Iyy_end = Iyy_end_outer-Iyy_end_inner
                 Izz_end = Izz_end_outer-Izz_end_inner
@@ -699,7 +704,8 @@ def RectangularFrustumMOI(La, Wa, Lb, Wb, H, p):
             # translate this submember's local inertia matrix to the PRP and add it to the total member's M_struc matrix
             self.M_struc += translateMatrix6to6DOF(Mmat, center_cap) # mass matrix of the member about the PRP
 
-
+        self.mshell = mshell
+        self.mfill  = mfill
         mass = self.M_struc[0,0]        # total mass of the entire member [kg]
         center = mass_center/mass       # total center of mass of the entire member from the PRP [m]
 

raft/raft_model.py

@@ -525,7 +525,9 @@ def solveStatics(self, case, display=0):
 
             X_initial[6*i:6*i+6] = np.array([fowt.x_ref, fowt.y_ref,0,0,0,0])
             fowt.setPosition(X_initial[6*i:6*i+6])      # zero platform offsets
-            fowt.calcStatics()
+            if case:
+                fowt.calcTurbineConstants(case, ptfm_pitch=0)  # for turbine forces >>> still need to update to use current fowt pose <<<
+            fowt.calcStatics() # Recompute statics because turbine heading may have changed due to yaw control
 
             if statics_mod == 0:
                 K_hydrostatic.append(fowt.C_struc + fowt.C_hydro)
@@ -541,8 +543,7 @@ def solveStatics(self, case, display=0):
                     if display > 1: 
                         print('Fowt ' + str(i))
                         print(case)
-                
-                fowt.calcTurbineConstants(case, ptfm_pitch=0)  # for turbine forces >>> still need to update to use current fowt pose <<<
+
                 fowt.calcHydroConstants()
                 F_env_constant[6*i:6*i+6] = np.sum(fowt.f_aero0, axis=1) + fowt.calcCurrentLoads(case)
 
@@ -639,6 +640,7 @@ def eval_func_equil(X, args):
                             case['wind_speed'] = caseorig['wind_speed'][i]
 
                         fowt.calcTurbineConstants(case, ptfm_pitch=r6[4])  # for turbine forces >>> still need to update to use current fowt pose <<<
+                        fowt.calcStatics() # Recompute statics because turbine heading may have changed due to yaw control
                         fowt.calcHydroConstants()  # prep for drag force and mean drift
 
                         Fnet[6*i:6*i+6] += np.sum(fowt.f_aero0, axis=1)  # sum mean turbine force across turbines                        
@@ -789,7 +791,7 @@ def step_func_equil(X, args, Y, oths, Ytarget, err, tol_, iter, maxIter):
 
         for i, fowt in enumerate(self.fowtList):
             print(f"Found mean offets of FOWT {i+1} with surge = {fowt.Xi0[0]: .2f} m,  sway  = {fowt.Xi0[1]: .2f},  and heave = {fowt.Xi0[2]: .2f} m")
-            print(f"                                 roll  = {fowt.Xi0[3]*180/np.pi: .2f} deg, pitch = {fowt.Xi0[4]*180/np.pi: .2f}, and yaw   = {fowt.Xi0[5]*180/np.pi: .2f} deg")
+            print(f"                                 roll  = {fowt.Xi0[3]*180/np.pi: .2f} deg, pitch = {fowt.Xi0[4]*180/np.pi: .2f} deg, and yaw   = {fowt.Xi0[5]*180/np.pi: .2f} deg")
 
 
         #dsolvePlot(info) # plot solver convergence trajectories