Merge pull request #37 from WISDEM/hydro_2nd_order

Second-order wave loads in RAFT

examples/example-WAMITQTF.py

@@ -0,0 +1,37 @@
+# example script for running RAFT with QTFs precomputed with a frequency-domain hydrodynamics code (in WAMIT format)
+
+import numpy as np
+import matplotlib.pyplot as plt
+import yaml
+import raft
+import os.path as path
+
+# open the design YAML file and parse it into a dictionary for passing to raft
+flNm = 'OC3spar-WAMITQTF'
+with open('./examples/' + flNm + '.yaml') as file:
+    design = yaml.load(file, Loader=yaml.FullLoader)
+
+# Create the RAFT model (will set up all model objects based on the design dict)
+model = raft.Model(design)  
+
+# Evaluate the system properties and equilibrium position before loads are applied
+model.analyzeUnloaded()
+
+# Compute natural frequencie
+model.solveEigen()
+
+# Simule the different load cases
+model.analyzeCases(display=1)
+
+# Plot the power spectral densities from the load cases
+model.plotResponses()
+
+# Visualize the system in its most recently evaluated mean offset position
+model.plot(hideGrid=True)
+
+model.saveResponses(path.join(model.fowtList[0].outFolderQTF, flNm))
+
+plt.show()
+
+# 0.02
+# 12.37

examples/example-slenderBodyQTF.py

@@ -0,0 +1,39 @@
+# example script for running RAFT with second-order loads computed internally with the slender-body approximation based on Rainey's equation
+
+import numpy as np
+import matplotlib.pyplot as plt
+import yaml
+import raft
+import os.path as path
+
+# open the design YAML file and parse it into a dictionary for passing to raft
+flNm = 'OC3spar-SlenderBodyQTF'
+with open('./examples/' + flNm + '.yaml') as file:
+    design = yaml.load(file, Loader=yaml.FullLoader)
+
+# Create the RAFT model (will set up all model objects based on the design dict)
+model = raft.Model(design)  
+
+# Evaluate the system properties and equilibrium position before loads are applied
+model.analyzeUnloaded()
+
+# Compute natural frequencie
+model.solveEigen()
+
+# Simule the different load cases
+model.analyzeCases(display=1)
+
+# Plot the power spectral densities from the load cases
+model.plotResponses()
+
+# Visualize the system in its most recently evaluated mean offset position
+model.plot(hideGrid=True)
+
+# Save the response to a given output folder
+outFolder = './examples/'
+model.saveResponses(path.join(outFolder, flNm))
+
+plt.show()
+
+# 0.02
+# 12.37

raft/helpers.py

@@ -121,7 +121,7 @@ def getWaveKin(zeta0, beta, w, k, h, r, nw, rho=1025.0, g=9.81):
         z = r[2]
 
         # only process wave kinematics if this node is submerged
-        if z < 0:
+        if z <= 0:
 
             # Calculate SINH( k*( z + h ) )/SINH( k*h ) and COSH( k*( z + h ) )/SINH( k*h ) and COSH( k*( z + h ) )/COSH( k*h )
             # given the wave number, k, water depth, h, and elevation z, as inputs.
@@ -153,6 +153,147 @@ def getWaveKin(zeta0, beta, w, k, h, r, nw, rho=1025.0, g=9.81):
 
     return u, ud, pDyn
 
+# Matrix of gradient of wave velocity
+def getWaveKin_grad_u1(w, k, beta, h, r):
+    grad = np.zeros([3, 3], dtype="complex")
+    x,y,z = r[0], r[1], r[2]
+
+    # More friendly notation
+    cosBeta = np.cos(deg2rad(beta))
+    sinBeta = np.sin(deg2rad(beta))
+    khz_xy = 0
+    khz_z = 0
+
+    if z <= 0 and k > 0:
+        # When k*h is too high, which happens for deep water/short waves, sinh(k*h) and cosh(k*h) become too large and are considered "inf".
+        # Hence, we chose a threshold of 10, above which the deep water approximation is employed.
+        if k*h >= 10:
+            khz_xy = np.exp(k*z)
+            khz_z = khz_xy
+        else:
+            khz_xy = np.cosh(k * (z + h)) / np.sinh(k*h)
+            khz_z = np.sinh(k * (z + h)) / np.sinh(k*h)
+
+
+            # u[0,i] =    w[i]* zeta[i]*COSHNumOvrSIHNDen *np.cos(beta)
+            # u[1,i] =    w[i]* zeta[i]*COSHNumOvrSIHNDen *np.sin(beta)
+            # u[2,i] = 1j*w[i]* zeta[i]*SINHNumOvrSIHNDen
+
+        # Along x direction
+        aux = w * cosBeta * np.exp( -1j*(k*(np.cos(beta)*r[0] + np.sin(beta)*r[1])))
+        grad[0,0] = -1j*aux * khz_xy * k*cosBeta # du/dx
+        grad[0,1] = -1j*aux * khz_xy * k*sinBeta # du/dy
+        grad[0,2] = aux * k * khz_z              # du/dz
+
+        # Along y direction
+        aux = w * sinBeta * np.exp( -1j*(k*(np.cos(beta)*r[0] + np.sin(beta)*r[1])))
+        grad[1,0] = grad[0,1]                    # dv/dx = du/dy
+        grad[1,1] = -1j*aux * khz_xy * k*sinBeta # dv/dy
+        grad[1,2] = aux * k * khz_z              # dv/dz
+
+        # Along z direction
+        aux = 1j * w * np.exp( -1j*(k*(np.cos(beta)*r[0] + np.sin(beta)*r[1])))
+        grad[2,0] = grad[0,2]         # dw/dx = du/dz
+        grad[2,1] = grad[0,1]         # dw/dy = dv/dz
+        grad[2,2] = aux * k * khz_xy  # dw/dz
+
+    return grad
+
+# Matrix of gradient of wave acceleration
+def getWaveKin_grad_dudt(w, k, beta, h, r):
+    return 1j*w*getWaveKin_grad_u1(w, k, beta, h, r)
+
+# Gradient of first order pressure
+def getWaveKin_grad_pres1st(k, beta, h, r, rho=1025, g=9.81):
+    grad = np.zeros(3, dtype="complex")
+    z = r[2]
+
+    # More friendly notation
+    cosBeta = np.cos(deg2rad(beta))
+    sinBeta = np.sin(deg2rad(beta))
+    khz_xy = 0
+    khz_z = 0
+
+    if z <= 0 and k > 0:
+        # When k*h is too high, which happens for deep water/short waves, sinh(k*h) and cosh(k*h) become too large and are considered "inf".
+        # Hence, we chose a threshold of 10, above which the deep water approximation is employed.
+        if k*h >= 10:
+            khz_xy = np.exp(k*z)
+            khz_z = khz_xy
+        else:
+            khz_xy = np.cosh(k * (z + h)) / np.cosh(k*h)
+            khz_z = np.sinh(k * (z + h)) / np.cosh(k*h)
+
+        grad[0] = rho*g*khz_xy*np.exp(-1j*(k*(cosBeta*r[0] + sinBeta*r[1])))*(-1j*k*cosBeta)
+        grad[1] = rho*g*khz_xy*np.exp(-1j*(k*(cosBeta*r[0] + sinBeta*r[1])))*(-1j*k*sinBeta)
+        grad[2] = rho*g*khz_z *np.exp(-1j*(k*(cosBeta*r[0] + sinBeta*r[1])))*k
+    return grad
+
+# Axial-divergence acceleration
+def getWaveKin_axdivAcc(w1, w2, k1, k2, beta1, beta2, h, r, vel1, vel2, q, g=9.81):
+    acc = np.zeros(3, dtype=complex)
+
+    aux = np.matmul(getWaveKin_grad_u1(w1, k1, beta1, h, r), q)
+    dwdz1 = np.dot(np.squeeze(aux), np.squeeze(q))
+    u1, _, _ = getWaveKin(np.ones([1,1]), beta1, w1*np.ones([1,1]), k1*np.ones([1,1]), h, r, 1, rho=1025.0, g=g)
+    u1 = np.squeeze(u1)
+
+    aux = np.matmul(getWaveKin_grad_u1(w2, k2, beta2, h, r), q)
+    dwdz2 = np.dot(np.squeeze(aux), np.squeeze(q))
+    u2, _, _ = getWaveKin(np.ones([1,1]), beta2, w2*np.ones([1,1]), k2*np.ones([1,1]), h, r, 1, rho=1025.0, g=g)
+    u2 = np.squeeze(u2)
+
+    vel1 -= np.dot(vel1, q) * q
+    vel2 -= np.dot(vel2, q) * q
+    u1   -= np.dot(u1, q) * q
+    u2   -= np.dot(u2, q) * q
+
+    acc = 0.25*(dwdz1*np.conj(u2 - vel2)+np.conj(dwdz2)*(u1 - vel1))
+
+    # There can't be any axial-divergence acceleration in the axial direction
+    acc -= np.dot(acc, q) * q
+
+    return acc
+
+# Acceleration and pressure due to second-order potential
+def getWaveKin_pot2ndOrd(w1, w2, k1, k2, beta1, beta2, h, r, g=9.81, rho=1025.0):
+    acc = np.zeros(3, dtype=complex)
+    p = 0+0j
+
+    if w1 == w2: # The difference-frequency second-order potential does not contribute to the mean forces 
+        return acc, p
+
+    b1 = deg2rad(beta1)
+    cosB1 = np.cos(b1)
+    sinB1 = np.sin(b1)
+
+    b2 = deg2rad(beta2)
+    cosB2 = np.cos(b2)
+    sinB2 = np.sin(b2)
+
+    z = r[2]
+
+    if (z <= 0 and k1 > 0 and k2 > 0):
+        k1_k2 = np.array([k1 * cosB1 - k2 * cosB2, k1 * sinB1 - k2 * sinB2, 0])
+        norm_k1_k2 = np.linalg.norm(k1_k2)
+
+        gamma_12 = (-1j*g/(2*w1)) * ( (k1**2)*(1-np.tanh(k1*h)**2) - 2*k1*k2*(1+np.tanh(k1*h)*np.tanh(k2*h)) ) / ((w1-w2)**2/g - norm_k1_k2*np.tanh(norm_k1_k2*h))
+        gamma_21 = (-1j*g/(2*w2)) * ( (k2**2)*(1-np.tanh(k2*h)**2) - 2*k2*k1*(1+np.tanh(k2*h)*np.tanh(k1*h)) ) / ((w2-w1)**2/g - norm_k1_k2*np.tanh(norm_k1_k2*h))
+        aux   = 0.5*(gamma_21+np.conj(gamma_12))
+
+        khz_xy = np.cosh(norm_k1_k2 * (z + h)) / np.cosh(norm_k1_k2 * h)
+        khz_z  = np.sinh(norm_k1_k2 * (z + h)) / np.cosh(norm_k1_k2 * h)
+
+        acc[0:2] = aux * khz_xy * np.exp(-1j*np.dot(k1_k2, r))
+        acc[0]  *= (w1 - w2) * (k1 * cosB1 - k2 * cosB2)
+        acc[1]  *= (w1 - w2) * (k1 * sinB1 - k2 * sinB2)
+
+        acc[2] = aux * khz_z * np.exp(-1j*np.dot(k1_k2, r))
+        acc[2] *= 1j * (w1 - w2) * norm_k1_k2
+
+        p = aux * khz_xy * np.exp(-1j*np.dot(k1_k2, r))
+        p *= -1j * rho * (w1 - w2)
+    return acc, p
 
 
 # calculate wave number based on wave frequency in rad/s and depth