diff --git a/raft/helpers.py b/raft/helpers.py index b6964a5..ed30b8e 100644 --- a/raft/helpers.py +++ b/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]): diff --git a/raft/raft_fowt.py b/raft/raft_fowt.py index a001e75..52aa6ac 100644 --- a/raft/raft_fowt.py +++ b/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,9 +1607,8 @@ 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): @@ -1620,8 +1616,8 @@ def calcQTF_slenderBody(self, waveHeadInd, Xi0=None, verbose=False, iCase=None): 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 diff --git a/raft/raft_member.py b/raft/raft_member.py index df00bfa..a01d6fe 100644 --- a/raft/raft_member.py +++ b/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) diff --git a/raft/raft_model.py b/raft/raft_model.py index 8bd969c..c0dbfca 100644 --- a/raft/raft_model.py +++ b/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,7 +1046,27 @@ 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) @@ -1060,7 +1074,8 @@ def solveDynamics(self, case, iCase, tol=0.01, conv_plot=0, RAO_plot=0): 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):