Merge branch 'rc1.3.2' into dev

This included updating true values of tests that were affected by bug fixes in rc1.3.2 but which were not included in that branch (e.g. tests with lumped mass moorings using the VolturnUS semi)

Author: lucas-carmo

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

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

@@ -1583,14 +1583,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
@@ -1601,6 +1600,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,:])

raft/raft_model.py

@@ -570,7 +570,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 self.staticsMod == 0:
                 K_hydrostatic.append(fowt.C_struc + fowt.C_hydro)
@@ -586,8 +588,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)
 
@@ -685,6 +686,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                        
@@ -847,7 +849,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

raft/raft_rotor.py

@@ -391,14 +391,12 @@ def setPosition(self, r6=np.zeros(6), R=None):
         # Store platform heading for use with nacelle yaw
         self.platform_heading = r6[5]
 
-        # Apply nacelle yaw depending on the yaw mode 
-        self.setYaw()
-        
         # Update RNA point locations [m] w.r.t. PRP in global orientations
         self.r_RRP_rel = np.matmul(self.R_ptfm, self.r_rel) # RNA ref point
-        self.r_CG_rel = self.r_RRP_rel + self.q*self.xCG_RNA # RNA CG location
-        self.r_hub_rel = self.r_RRP_rel + self.q*self.overhang # rotor hub location
-        
+
+        # Apply nacelle yaw depending on the yaw mode 
+        self.setYaw()
+                
         '''
         self.r_RRP = ? # RNA reference point
         self.r_CG  = ? # RNA CG location
@@ -456,6 +454,10 @@ def setYaw(self, yaw=None):
         # Compute shaft axis unit vector in FOWT and global frames 
         self.q_rel = np.matmul(R_q_rel, np.array([1,0,0]) )
         self.q = np.matmul(self.R_ptfm, self.q_rel) # Write in the global frame 
+
+        # Update RNA point locations [m] w.r.t. PRP in global orientations
+        self.r_CG_rel = self.r_RRP_rel + self.q*self.xCG_RNA # RNA CG location
+        self.r_hub_rel = self.r_RRP_rel + self.q*self.overhang # rotor hub location
 
         return self.yaw
 
@@ -636,7 +638,7 @@ def calcHydroConstants(self, dgamma=0, rho=1025, g=9.81):
         return A_hydro, I_hydro
 
 
-    def calcCavitation(self, case, azimuth=0, clearance_margin=1.0, Patm=101325, Pvap=2500, error_on_cavitation=False):
+    def calcCavitation(self, case, azimuth=0, clearance_margin=1.0, Patm=101325, Pvap=2300, error_on_cavitation=False):
         ''' Method to calculate the cavitation number of the rotor
         (wind speed (m/s), rotor speed (RPM), pitch angle (deg), azimuth (deg))
 
@@ -838,7 +840,7 @@ def calcAero(self, case, current=False, display=0):
         # Set up vectors in axis frame. Assuming CCBlade forces (but not 
         # moments) are relative to the rotor axis
         forces_axis = np.array([loads["T"][0], loads["Y"][0], loads["Z" ][0]])
-        moments_axis = np.array([loads["My"][0], loads["Q"][0], loads["Mz"][0]])        
+        moments_axis = np.array([loads["Q"][0], loads["My"][0], loads["Mz"][0]])        
 
         # Rotate forces and moments to be relative to global orientation (but still wrt hub)
         self.f0[:3] = np.matmul(self.R_q, forces_axis)

tests/test_data/IEA15MW.yaml

@@ -14,10 +14,10 @@ site:
     shearExp    : 0.12       #          shear exponent
 
 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]
     Fthrust       :       1500.0E3    #  [N]        temporary thrust force to use