Merge pull request #29 from NREL/dynamics

Dynamic tensions on mooring lines

examples/dynamic.py

@@ -0,0 +1,129 @@
+# Hasty example with a Subsystem and the new dynamic features
+
+import numpy as np
+import matplotlib.pyplot as plt
+import moorpy as mp
+
+
+# ===== Create a MoorPy system to test with =====
+
+# ----- choose some system geometry parameters -----
+
+depth     = 200                             # water depth [m]
+angles    = np.radians([60, 300])      # line headings list [rad]
+rAnchor   = 600                            # anchor radius/spacing [m]
+zFair     = -21                             # fairlead z elevation [m]
+rFair     = 20                              # fairlead radius [m]
+lineLength= 650                            # line unstretched length [m]
+typeName  = "chain1"                        # identifier string for the line type
+
+
+# ----- Now a Subsystem in a larger System with a floating body -----
+
+# Create new MoorPy System and set its depth
+ms = mp.System(depth=depth)
+
+# add a line type
+ms.setLineType(dnommm=120, material='chain', name=typeName, Cd=1.4, Ca=1, CdAx=0.4, CaAx=0)  # this would be 120 mm chain
+
+# Add a free, body at [0,0,0] to the system (including some properties to make it hydrostatically stiff)
+ms.addBody(0, np.zeros(6), m=1e6, v=1e3, rM=100, AWP=1e3)
+
+# For each line heading, set the anchor point, the fairlead point, and the line itself
+for i, angle in enumerate(angles):
+
+    # create end Points for the line
+    ms.addPoint(1, [rAnchor*np.cos(angle), rAnchor*np.sin(angle), -depth])   # create anchor point (type 0, fixed)
+    ms.addPoint(1, [  rFair*np.cos(angle),   rFair*np.sin(angle),  zFair])   # create fairlead point (type 0, fixed)
+
+    # attach the fairlead Point to the Body (so it's fixed to the Body rather than the ground)
+    ms.bodyList[0].attachPoint(2*i+2, [rFair*np.cos(angle), rFair*np.sin(angle), zFair]) 
+
+    # add a Line going between the anchor and fairlead Points
+    ms.addLine(lineLength, typeName, pointA=2*i+1, pointB=2*i+2)
+
+# ----- Now add a SubSystem line! -----
+ss = mp.Subsystem(mooringSys=ms, depth=depth, spacing=rAnchor, rBfair=[10,0,-20])
+
+# set up the line types
+ms.setLineType(180, 'chain', name='one', Cd=1.4, Ca=1, CdAx=0.4, CaAx=0)
+ms.setLineType( 50, 'chain', name='two', Cd=1.4, Ca=1, CdAx=0.4, CaAx=0)
+
+
+# set up the lines and points and stuff
+ls = [350, 300]
+ts = ['one', 'two']
+ss.makeGeneric(lengths=ls, types=ts)
+ss.lineList[1].nNodes = 10  # reduce nodes on rope line for easier viewing
+ss.initialize()  # update things after changing node number
+
+# add points that the subSystem will attach to...
+ms.addPoint(1, [-rAnchor, 100, -depth])   # Point 5 - create anchor point (type 0, fixed)
+ms.addPoint(1, [ -rFair ,   0,  zFair])   # Point 6 - create fairlead point (type 0, fixed)
+ms.bodyList[0].attachPoint(6, [-rFair, 0, zFair])  # attach the fairlead Point to the Body 
+
+# string the Subsystem between the two points!
+ms.lineList.append(ss)  # add the SubSystem to the System's lineList
+ss.number = 3
+ms.pointList[4].attachLine(3, 0)  # attach it to the respective points
+ms.pointList[5].attachLine(3, 1)  # attach it to the respective points
+
+
+# ----- run the model to check that the Subsystem is working -----
+
+ms.initialize()                                             # make sure everything's connected
+
+ms.solveEquilibrium()                                       # equilibrate
+fig, ax = ms.plot()                                         # plot the system in original configuration
+
+ms.bodyList[0].f6Ext = np.array([3e6, 0, 0, 0, 0, 0])       # apply an external force on the body 
+ms.solveEquilibrium()                                      # equilibrate
+fig, ax = ms.plot(ax=ax, color='red')                       # plot the system in displaced configuration (on the same plot, in red)
+
+print(f"Body offset position is {ms.bodyList[0].r6}")
+
+
+
+# ===== Now try out Serag's dynamic tension functions with it =====
+
+# dynamic solve of some lines
+kbot = 3E+06
+cbot = 3E+05
+
+moorpy_freqs = []
+moorpy_fairten_psds = []
+moorpy_ten_stds = []
+
+omegas = np.array([ 0.02,  0.04,  0.06,  0.08 ])
+S_zeta = np.array([ 10.0 ,  10.0 ,  10.0 ,  10.0  ])
+RAO_fl = np.array([[ 2.0 ,  0.0 ,  0.0 ],
+                   [ 2.0 ,  0.0 ,  0.0 ],
+                   [ 2.0 ,  0.0 ,  0.0 ],
+                   [ 2.0 ,  0.0 ,  0.0 ]])
+
+T_nodes_psd_fd,T_nodes_std_fd,s,r_static,r_dynamic,r_total,X = ms.lineList[1].dynamicSolve(
+                                omegas,S_zeta,RAO_A=0,RAO_B=RAO_fl,depth=-ms.depth,
+                                kbot=kbot,cbot=cbot, seabed_tol=1e-4, tol=0.01, iters=100, w=0.8)
+
+fig2,ax2 = plt.subplots(1,1)
+
+plt.plot(T_nodes_std_fd,'-r',label = 'MoorPy (FD)')
+plt.xlabel('Node number (anchor to fairlead)')
+plt.ylabel('Fairlead tension std. dev [$N$]')
+
+# now try a Subsystem
+T_nodes_psd_fd,T_nodes_std_fd,s,r_static,r_dynamic,r_total,X = ms.lineList[2].dynamicSolve(
+                                omegas,S_zeta,RAO_A=0,RAO_B=RAO_fl,depth=-ms.depth,
+                                kbot=kbot,cbot=cbot, seabed_tol=1e-4, tol=0.01, iters=100, w=0.8)
+
+fig2,ax2 = plt.subplots(1,1)
+plt.plot(T_nodes_std_fd,'-r',label = 'MoorPy (FD)')
+plt.xlabel('Node number (anchor to fairlead)')
+plt.ylabel('Fairlead tension std. dev [$N$]')
+
+
+# Some mode shape plots
+ms.lineList[1].getModes(plot_modes=7, amp_factor=1000)  # with a Line
+ms.lineList[2].getModes(plot_modes=7, amp_factor=1000)  # with a Subsystem
+
+plt.show()

examples/lumped_mass.py

@@ -0,0 +1,74 @@
+import moorpy as mp
+import os
+import numpy as np
+from moorpy.helpers import lines2subsystem, lines2ss
+try:
+    import pickle5 as pickle
+except:
+    import pickle
+import matplotlib.pyplot as plt
+
+def JONSWAP(ws, Hs, Tp, Gamma=None):
+    '''JONSWAP function copied from RAFT
+    '''
+    # If peak shape parameter gamma is not specified, use the recommendation 
+    # from IEC 61400-3 as a function of Hs and Tp. For PM spectrum, use 1.
+    if not Gamma:
+        TpOvrSqrtHs = Tp/np.sqrt(Hs)
+        if TpOvrSqrtHs <= 3.6:
+            Gamma = 5.0
+        elif TpOvrSqrtHs >= 5.0:
+            Gamma = 1.0
+        else:
+            Gamma = np.exp( 5.75 - 1.15*TpOvrSqrtHs )
+
+    # handle both scalar and array inputs
+    if isinstance(ws, (list, tuple, np.ndarray)):
+        ws = np.array(ws)
+    else:
+        ws = np.array([ws])
+
+    # initialize output array
+    S = np.zeros(len(ws))
+
+
+    # the calculations
+    f        = 0.5/np.pi * ws                         # wave frequencies in Hz
+    fpOvrf4  = pow((Tp*f), -4.0)                      # a common term, (fp/f)^4 = (Tp*f)^(-4)
+    C        = 1.0 - ( 0.287*np.log(Gamma) )          # normalizing factor
+    Sigma = 0.07*(f <= 1.0/Tp) + 0.09*(f > 1.0/Tp)    # scaling factor
+
+    Alpha = np.exp( -0.5*((f*Tp - 1.0)/Sigma)**2 )
+
+    return  0.5/np.pi *C* 0.3125*Hs*Hs*fpOvrf4/f *np.exp( -1.25*fpOvrf4 )* Gamma**Alpha
+
+
+
+current_dir = os.path.dirname(os.path.abspath(__file__))
+ms = mp.System(os.path.join(current_dir, 'volturn_chain.dat'))
+ms.initialize()
+ms.solveEquilibrium()
+# ms = lines2ss(ms) # For cases with multisegment lines, need to convert each of them to a subsystem
+
+# Updates the dynamic matrices of all the lines in the system.
+# This function can only properly update the inertia, added mass, and stiffness matrices of each line.        
+# Though it updates the damping matrix, this is done considering unitary amplitude motions of the nodes and
+# no fluid kinematics, so it is not correct.
+ms.updateSystemDynamicMatrices() 
+M, A, B, K_dyn = ms.getCoupledDynamicMatrices() # Get the dynamic matrices of the system
+K_qsA  = ms.getCoupledStiffnessA(lines_only=True) # We can also compute the stiffness matrix without the lumped mass model. K_dyn should be similar to K_qsa
+
+
+# Get the dynamic tension along the line
+line = ms.lineList[0] # Get the line object
+RAO_data = pickle.load(open(os.path.join(current_dir, 'RAO_fl.pkl'), 'rb')) # Read the nFreq x 3 RAO matrix (nFreq x 4 complex numpy array, first column are the frequencies in rad/s)
+RAO_fl = RAO_data[:, 1:] # Motion RAOs of the fairlead
+w = RAO_data[:, 0] # Frequencies of the RAO data
+Sw = JONSWAP(ws = w, Hs = 6, Tp = 8) # Arbitrary wave spectrum
+T_nodes_amp, T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X = line.dynamicSolve(w, Sw, RAO_A=0,RAO_B=RAO_fl, depth=np.abs(line.rA[2]))
+
+fig, ax = plt.subplots(1, 1)
+ax.plot(w, T_nodes_psd[:,-1], '-k')
+ax.set_xlabel('Frequency (rad/s)')
+ax.set_ylabel('PSD fairlead tension (N^2.s/rad)')
+plt.show()