Changed output QTF names + mean forces in KAY + 2nd ord in lineariz loop

- The QTFs computed by RAFT now include the case number and the turbine number. Similar to other parts of the code, the case number starts at 1. The turbine number starts at 0, but perhaps we should changed that across the code for consistency.
- The correction for second-order wave diffraction (based on the paper by Kim and Yue, 1990) now also includes the mean drift. It is a simple addition, we just needed to do a limit for a term that goes to 0/0 when w1-w2=0
- We changed the way that the second-order load motions are included in the linearization loop. Now, we perform the linearization twice. First, we perform the linearization loop WITHOUT the second-order loads/motions. After the loop converges, we compute the QTFs, the second-order loads, and then the motions. Since these motions may be very relevant, we redo the linearization loop including the second-order motions. This is only done for the first wave of the case, like before.

raft/helpers.py

@@ -677,15 +677,15 @@ def getRAO(Xi, zeta):
     if Xi.shape[-1] != len(zeta):
         raise Exception("The last dimension of Xi must be the same length as zeta")
 
-    # Is there an eps value to use instead?
-    idx = np.where(np.abs(zeta)>1e-20)    
+    # # Is there an eps value to use instead?
+    idx = np.where(np.abs(zeta)>1e-6)    
 
     # Reshape zeta to be able to broadcast along the frequency dimension
     RAO = np.zeros_like(Xi, dtype=complex)
     RAO[..., idx] = Xi[..., idx] / zeta[idx]
-
     return RAO
 
+
 def printMat(mat):
     '''Print a matrix'''
     for i in range(mat.shape[0]):

raft/raft_fowt.py

@@ -1357,7 +1357,7 @@ def calcCurrentLoads(self, case):
         return D_hydro
 
 
-    def calcQTF_slenderBody(self, waveHeadInd, Xi0=None, verbose=False, iCase=None):
+    def calcQTF_slenderBody(self, waveHeadInd, Xi0=None, verbose=False, iCase=None, iWT=None):
         '''This computes the second-order strip-theory-hydrodynamics terms
            Only long crested seas are currently considered here.
            Options:
@@ -1369,7 +1369,8 @@ def calcQTF_slenderBody(self, waveHeadInd, Xi0=None, verbose=False, iCase=None):
         # TODO: I think a lot of the calculations below should move to the member class
 
         if Xi0 is None:
-            Xi0 = np.zeros([self.nDOF, len(self.w)], dtype=complex)        
+            Xi0 = np.zeros([self.nDOF, len(self.w)], dtype=complex)
+        # Xi0 = np.zeros([self.nDOF, len(self.w)], dtype=complex)
         rho = self.rho_water
         g   = self.g
 
@@ -1386,8 +1387,8 @@ def calcQTF_slenderBody(self, waveHeadInd, Xi0=None, verbose=False, iCase=None):
         # Print Xi in the same format as a WAMIT .4 file for comparison
         if self.outFolderQTF is not None and verbose:
             # Check if the file exists. If so, append a number to the end of the file name
-            if isinstance(iCase, int):
-                outPath = os.path.join(self.outFolderQTF, f"raos-slender_body_Head{whead}_Case{iCase}.4")
+            if isinstance(iCase, int) and isinstance(iWT, int):
+                outPath = os.path.join(self.outFolderQTF, f"raos-slender_body_Head{whead}_Case{iCase+1}_WT{iWT}.4")
             else:
                 outPath = os.path.join(self.outFolderQTF, f"raos-slender_body_Head{whead}.4")
 
@@ -1416,7 +1417,7 @@ def calcQTF_slenderBody(self, waveHeadInd, Xi0=None, verbose=False, iCase=None):
         self.qtf = np.zeros([len(self.w1_2nd), len(self.w2_2nd), 1, self.nDOF], dtype=complex)  # Need this fourth dimension for conformity with the case where the QTFs are read from a file
 
         if verbose:
-            print(f"Computing QTF for heading {beta:.2f}")
+            print(f" Computing QTF for heading {beta:.2f}")
 
         # Start with the force due to rotation of first-order wave forces
         # This component depends on the forces on the whole body, hence it is outside of the member loop
@@ -1500,10 +1501,6 @@ def calcQTF_slenderBody(self, waveHeadInd, Xi0=None, verbose=False, iCase=None):
                         Ca_p1  = np.interp( mem.ls[il], mem.stations, mem.Ca_p1 )
                         Ca_p2  = np.interp( mem.ls[il], mem.stations, mem.Ca_p2 )
                         Ca_End = np.interp( mem.ls[il], mem.stations, mem.Ca_End)
-                        # Ca_p1=1
-                        # Ca_p2=1
-                        # Ca_End=1
-
 
                         # ----- compute side effects ---------------------------------------------------------
                         if circ:
@@ -1610,18 +1607,17 @@ def calcQTF_slenderBody(self, waveHeadInd, Xi0=None, verbose=False, iCase=None):
                     qtf_nabla[i1,i2,waveHeadInd,:]  += F_nabla
                     qtf_rslb[i1,i2,waveHeadInd,:]   += F_rslb
 
-                    # Add Kim and Yue correction for the horizontal force
-                    # We add only the real part because the imaginary part of the force is already taken care by the slender-body approximation                                   
-                    self.qtf[i1,i2,waveHeadInd,:] += np.real(mem.correction_KAY(self.depth, w1, w2, beta, rho=rho, g=g, k1=k1, k2=k2, Nm=10))
+                    # Add Kim and Yue correction                                                     
+                    self.qtf[i1,i2,waveHeadInd,:] += mem.correction_KAY(self.depth, w1, w2, beta, rho=rho, g=g, k1=k1, k2=k2, Nm=10)
 
         # Need a complete matrix due to interpolations that are used to computed the forces
         for i in range(self.nDOF):
             self.qtf[:,:, waveHeadInd, i] = self.qtf[:,:, waveHeadInd,i] + np.conj(self.qtf[:,:, waveHeadInd,i]).T - np.diag(np.diag(np.conj(self.qtf[:,:, waveHeadInd,i])))
 
         if self.outFolderQTF is not None and verbose:
             # Check if the file exists. If so, append a number to the end of the file name
-            if isinstance(iCase, int):
-                outPath = os.path.join(self.outFolderQTF, f"qtf-slender_body-total_Head{whead}_Case{iCase}.12d")
+            if isinstance(iCase, int) and isinstance(iWT, int):
+                outPath = os.path.join(self.outFolderQTF, f"qtf-slender_body-total_Head{whead}_Case{iCase+1}_WT{iWT}.12d")
             else:
                 outPath = os.path.join(self.outFolderQTF, f"qtf-slender_body-total_Head{whead}.12d")
             self.writeQTF(self.qtf, outPath)
@@ -1685,7 +1681,7 @@ def writeQTF(self, qtfIn, outPath, w=None):
                             f.write(f"{2*np.pi/w1[i1]: 8.4e} {2*np.pi/w2[i2]: 8.4e} {self.heads_2nd[ih]: 8.4e} {self.heads_2nd[ih]: 8.4e} {iDoF+1} {np.abs(F): 8.4e} {np.angle(F): 8.4e} {F.real: 8.4e} {F.imag: 8.4e}\n")
 
 
-    def calcHydroForce_2ndOrd(self, beta, S0):
+    def calcHydroForce_2ndOrd(self, beta, S0, iCase=None, iWT=None):
         ''' Compute force due to 2nd order hydrodynamic loads
         See Pinkster (1980), Section IV.3
         We are losing the phases when computing forces from the spectrum.
@@ -1759,6 +1755,16 @@ def calcHydroForce_2ndOrd(self, beta, S0):
         # Displace f by one frequency so that it aligns with the frequency vector
         f[:, 0:-1] = f[:, 1:]
         f[:, -1] = 0
+
+        ident = 'WAMIT'
+        if self.potSecOrder == 1:
+            ident = 'slenderBody'
+
+        if self.outFolderQTF is not None:             
+            with open(os.path.join(self.outFolderQTF, f'f_2nd-{ ident }_Case{ iCase }_WT{ iWT }.txt'), 'w') as file:
+                for w, frow in zip(self.w, f.T):
+                    file.write(f'{w:.5f} {frow[0]:.5f} {frow[1]:.5f} {frow[2]:.5f} {frow[3]:.5f} {frow[4]:.5f} {frow[5]:.5f}\n')
+
         return f_mean, f
 
 

raft/raft_member.py

@@ -97,7 +97,7 @@ def __init__(self, mi, nw, BEM=[], heading=0):
         # We disable the MCF correction if the element is not circular
         if self.MCF:
             if not self.shape=='circular':
-                print(f'MacCamy-Fuchs correction not applicable to member {self.name}. Member needs to be circular. Disabling MCF correction.')
+                print(f'MacCamy-Fuchs correction not applicable to member {self.name}. Member needs to be circular. Disabling MCF.')
                 self.MCF = False
 
 
@@ -1135,10 +1135,6 @@ def omega(k1R, k2R, n):
         F = np.zeros(6, dtype=complex)
         if not self.MCF:
             return F
-
-        # Check if member is circular
-        if self.shape != 'circular':
-            return F
 
         # Compute k1 and k2 if not provided
         if k1 is None:
@@ -1170,11 +1166,12 @@ def omega(k1R, k2R, n):
             # Compute force lumped at the intersection with the mean waterline
             k1R, k2R = k1*R, k2*R
             Fwl = 0+0j
-            if k1R >= np.pi/10 or k2R >= np.pi/10:
-                for nn in range(Nm + 1):
-                    Fwl += -rho*g*R*2j/np.pi/(k1R*k2R) * omega(k1R, k2R, nn)
+            # if k1R >= np.pi/10 or k2R >= np.pi/10:
+            for nn in range(Nm + 1):
+                Fwl += -rho*g*R*2j/np.pi/(k1R*k2R) * omega(k1R, k2R, nn)
 
-            Fwl *= np.exp(-1j*np.dot(k1_k2, rwl))
+            Fwl = np.real(Fwl) # Get only the part related to diffraction effects to avoid double counting with Rainey's equation
+            Fwl *= np.exp(-1j*np.dot(k1_k2, rwl)) # Solution considers cylinder at (0,0). Displace it to actual location
             F += translateForce3to6DOF(Fwl*pforce, rwl)
 
 
@@ -1210,30 +1207,28 @@ def omega(k1R, k2R, n):
                 # We do not want to use the diffraction correction for very long waves
                 # because it is already well captured by the slender-body approximation
                 # and because it could render the results worse.
-                if k1R < np.pi/10 and k2R < np.pi/10:
-                    continue
+                # if k1R < np.pi/10 and k2R < np.pi/10:
+                #     continue
 
                 # For mean loads
                 dF = 0+0j  # Force per unit length. Assumed to be aligned with the wave propagation direction
-                if w1 == w2:                
-                    for nn in range(Nm + 1):
-                        dF += 0 # TODO: Implement the mean load
+                if w1 == w2:
+                    Im = 0.5 * (np.sinh((k1 + k2)*(z2+h)) / (k1h + k2h) - (z2+h)/h - np.sinh((k1 + k2)*(z1+h)) / (k1h + k2h) + (z1+h)/h)
+                    Ip = 0.5 * (np.sinh((k1 + k2)*(z2+h)) / (k1h + k2h) + (z2+h)/h - np.sinh((k1 + k2)*(z1+h)) / (k1h + k2h) - (z1+h)/h)                    
 
-                # For difference-frequency loads
                 else:
-                    coshk1h, coshk2h = np.cosh(k1h), np.cosh(k2h)                    
                     Im = 0.5 * (np.sinh((k1 + k2)*(z2+h)) / (k1h + k2h) - np.sinh((k1 - k2)*(z2+h)) / (k1h - k2h) - np.sinh((k1 + k2)*(z1+h)) / (k1h + k2h) + np.sinh((k1 - k2)*(z1+h)) / (k1h - k2h))
-                    Ip = 0.5 * (np.sinh((k1 + k2)*(z2+h)) / (k1h + k2h) + np.sinh((k1 - k2)*(z2+h)) / (k1h - k2h) - np.sinh((k1 + k2)*(z1+h)) / (k1h + k2h) - np.sinh((k1 - k2)*(z1+h)) / (k1h - k2h))
-                    for nn in range(Nm + 1):
-                        coshk1h = np.cosh(k1h)
-                        coshk2h = np.cosh(k2h)
+                    Ip = 0.5 * (np.sinh((k1 + k2)*(z2+h)) / (k1h + k2h) + np.sinh((k1 - k2)*(z2+h)) / (k1h - k2h) - np.sinh((k1 + k2)*(z1+h)) / (k1h + k2h) - np.sinh((k1 - k2)*(z1+h)) / (k1h - k2h))                    
 
-                        dF += rho*g*R*2j/np.pi/(k1R*k2R) * omega(k1R, k2R, nn) * (k1h*k2h/np.sqrt(k1h*np.tanh(k1h))/np.sqrt(k2h*np.tanh(k2h)) * (Im + Ip*nn*(nn+1)/k1R/k2R)/coshk1h/coshk2h)
+                coshk1h, coshk2h = np.cosh(k1h), np.cosh(k2h)                    
+                for nn in range(Nm + 1):
+                    dF += rho*g*R*2j/np.pi/(k1R*k2R) * omega(k1R, k2R, nn) * (k1h*k2h/np.sqrt(k1h*np.tanh(k1h))/np.sqrt(k2h*np.tanh(k2h)) * (Im + Ip*nn*(nn+1)/k1R/k2R)/coshk1h/coshk2h)
 
 
                 # The force calculated above considers a cylinder at (x,y)=(0,0), so we need to account for the wave phase
                 r = 0.5*(r1 + r2)
-                dF *= np.exp(-1j*np.dot(k1_k2, r))
+                dF = np.real(dF) # Get only the part related to diffraction effects to avoid double counting with Rainey's equation
+                dF *= np.exp(-1j*np.dot(k1_k2, rwl)) # Solution considers cylinder at (0,0). Displace it to actual location
 
                 # Project it along the the directino computed above and add it to the total force
                 F += translateForce3to6DOF(dF*pforce, r)

raft/raft_model.py

@@ -976,25 +976,31 @@ def solveDynamics(self, case, iCase, tol=0.01, conv_plot=0, RAO_plot=0):
                 B_turb = np.zeros([6,6,self.nw])
 
             # >>>> NOTE: Turbulent wind excitation is currently disabled pending formulation checks/fixes <<<<
-            print('Solving for system response to wave excitation in primary wave direction') 
+            print('Solving for system response to wave excitation in primary wave direction')
 
             # We can compute second-order hydrodynamic forces here if they are calculated using external QTF file
-            # In some cases, they may be very relevant to the motion RMS values (e.g. pitch motion of spar platforms), so should be included in the drag linearization process
+            # In some cases, they may be very relevant to the motion RMS values (e.g. pitch motion of spar platforms), so should be included in the drag linearization process            
             fowt.Fhydro_2nd = np.zeros([fowt.nWaves, fowt.nDOF, fowt.nw], dtype=complex) 
             fowt.Fhydro_2nd_mean = np.zeros([fowt.nWaves, fowt.nDOF])
             if fowt.potSecOrder==2:
-                fowt.Fhydro_2nd_mean[0, :], fowt.Fhydro_2nd[0, :, :] = fowt.calcHydroForce_2ndOrd(fowt.beta[0], fowt.S[0,:])
+                fowt.Fhydro_2nd_mean[0, :], fowt.Fhydro_2nd[0, :, :] = fowt.calcHydroForce_2ndOrd(fowt.beta[0], fowt.S[0,:], iCase=iCase, iWT=i)
+
+            # We use this flag when fowt.potSecOrder==1, so that we compute the QTFs including first-order motions only
+            flagComputedQTF = False
+
 
             # sum up all linear (non-varying) matrices up front, including potential summation across multiple rotors
             M_lin.append( M_turb + fowt.M_struc[:,:,None] + fowt.A_BEM + fowt.A_hydro_morison[:,:,None]        ) # mass
             B_lin.append( B_turb + fowt.B_struc[:,:,None] + fowt.B_BEM + np.sum(fowt.B_gyro, axis=2)[:,:,None] ) # damping
             C_lin.append(          fowt.C_struc   + fowt.C_moor        + fowt.C_hydro                          ) # stiffness
             F_lin.append( fowt.F_BEM[0,:,:] + fowt.F_hydro_iner[0,:,:] + fowt.Fhydro_2nd[0, :, :]) # consider only excitation from the primary sea state in the load case for now
+
 
             #F_lin = np.sum(fowt.F_aero, axis=2) + fowt.F_BEM + np.sum(fowt.F_hydro_iner, axis=0) # excitation
 
             # start fixed point iteration loop for dynamics of the individual FOWT
-            for iiter in range(nIter):
+            iiter = 0
+            while iiter < nIter:
 
                 # initialize/zero total system coefficient arrays
                 M_tot = np.zeros([fowt.nDOF,fowt.nDOF,self.nw])       # total mass and added mass matrix [kg, kg-m, kg-m^2]
@@ -1029,18 +1035,6 @@ def solveDynamics(self, case, iCase, tol=0.01, conv_plot=0, RAO_plot=0):
                     # solve response (complex amplitude)
                     Xi[:,ii] = np.linalg.solve(Z[:,:,ii], F_tot[:,ii])
 
-                # If we are computing the QTFs internally, we need to consider the motions induced by first-order hydrodynamic forces.
-                # Rigorously, the QTFs should be computed at each iteration, but that would be very expensive. 
-                # We will assume that the QTFs won't change much with the first-order motions within this linearization loop, so they are computed only once.
-                # In any case, they are reevaluated after this loop considering the final first-order motions.
-                if fowt.potSecOrder == 1 and iiter == 0:
-                    print(' Computing QTF for drag linearization process.') 
-
-                    Xi0 = getRAO(Xi[i1:i2, :], fowt.zeta[0,:])
-                    fowt.calcQTF_slenderBody(waveHeadInd=0, Xi0=Xi0, verbose=False)
-                    _, fowt.Fhydro_2nd[0, :, :] = fowt.calcHydroForce_2ndOrd(fowt.beta[0], fowt.S[0,:])
-                    F_lin[i1:i2] += fowt.Fhydro_2nd[0, :, :]
-
                 if conv_plot:
                     # Convergence Plotting
                     # plots of surge response at each iteration for observing convergence
@@ -1052,15 +1046,36 @@ def solveDynamics(self, case, iCase, tol=0.01, conv_plot=0, RAO_plot=0):
                 tolCheck = np.abs(Xi - XiLast) / ((np.abs(Xi)+tol))
                 if (tolCheck < tol).all():
                     print(f" Iteration {iiter}, converged (largest change is {np.max(tolCheck):.5f} < {tol})")
-                    break
+
+                    if fowt.potSecOrder != 1 or flagComputedQTF:
+                        break
+                    else:
+                        iiter = 0
+
+                        # If we are computing the QTFs internally, we need to consider the motions induced by first-order hydrodynamic forces.
+                        # So, we compute the QTFs considering the first-order motions computed with the linearized drag and then we loop again
+                        # to make sure the second-order motions are considered in the linearization procedure. 
+                        # This is important for cases where the second-order motions are large compared to the first-order motions.
+                        print(f"Resolving for system response in primary wave direction, now with second-order wave loads.")
+                        Xi0 = getRAO(Xi[i1:i2, :], fowt.zeta[0,:])
+
+                        tic = time.perf_counter()                        
+                        fowt.calcQTF_slenderBody(waveHeadInd=0, Xi0=Xi0, verbose=True, iCase=iCase, iWT=i)
+                        toc = time.perf_counter()
+                        print(f"\n Time to compute QTFs for fowt {i}: {toc - tic:0.4f} seconds")
+
+                        fowt.Fhydro_2nd_mean[0, :], fowt.Fhydro_2nd[0, :, :] = fowt.calcHydroForce_2ndOrd(fowt.beta[0], fowt.S[0,:], iCase=iCase, iWT=i)
+                        F_lin[i1:i2] += fowt.Fhydro_2nd[0, :, :]
+                        flagComputedQTF = True                         
                 else:
                     XiLast = 0.2*XiLast + 0.8*Xi    # use a mix of the old and new response amplitudes to use for the next iteration
                                                     # (uses hard-coded successive under relaxation for now)
                     print(f" Iteration {iiter}, unconverged (largest change is {np.max(tolCheck):.5f} >= {tol})")
 
                 if iiter == nIter-1:
                     print("WARNING - solveDynamics iteration did not converge to the tolerance.")
-
+
+                iiter += 1
 
             if conv_plot:
                 # labels for convergence plots
@@ -1140,15 +1155,13 @@ def solveDynamics(self, case, iCase, tol=0.01, conv_plot=0, RAO_plot=0):
             for i, fowt in enumerate(self.fowtList):
                 i1, i2 = i*6, i*6+6
                 if fowt.potSecOrder == 1:
-                    Xi0 = getRAO(self.Xi[ih,i1:i2, :], fowt.zeta[ih,:])
-
-                    tic = time.perf_counter()
-                    fowt.calcQTF_slenderBody(waveHeadInd=ih, Xi0=Xi0, verbose=True, iCase=iCase)
-                    toc = time.perf_counter()
-                    print(f"\n Time to compute QTFs for fowt {i} and wave {ih}: {toc - tic:0.4f} seconds") 
-
-                    fowt.Fhydro_2nd_mean[ih, :], fowt.Fhydro_2nd[ih, :, :] = fowt.calcHydroForce_2ndOrd(fowt.beta[ih], fowt.S[ih,:])
-
+                    # Don't recompute the QTFs for the first wave because it was already done above.
+                    # Also, we would end up including second-order motions if we computed it again.
+                    if ih > 0: 
+                        Xi0 = getRAO(self.Xi[ih,i1:i2, :], fowt.zeta[ih,:])                        
+                        fowt.calcQTF_slenderBody(waveHeadInd=ih, Xi0=Xi0, verbose=True, iCase=iCase, iWT=i)                        
+                        fowt.Fhydro_2nd_mean[ih, :], fowt.Fhydro_2nd[ih, :, :] = fowt.calcHydroForce_2ndOrd(fowt.beta[ih], fowt.S[ih,:])
+
                     # Recompute the wave excitation forces and consequent motions to include second-order hydrodynamic forces
                     F_wave[i1:i2] = fowt.F_BEM[ih,:,:] + fowt.F_hydro_iner[ih,:,:] + F_linearized + fowt.Fhydro_2nd[ih, :, :]
                     for iw in range(self.nw):