From 45d2827863b4067c2f0195c9dcf992fcb44a78d0 Mon Sep 17 00:00:00 2001 From: serag_abdelmoteleb Date: Wed, 3 Jan 2024 16:22:23 +0100 Subject: [PATCH 01/14] Add: Line dynamcis + new CompositeLine class --- moorpy/__init__.py | 2 + moorpy/compositeLine.py | 111 +++++++ moorpy/dynamic_tension_functions.py | 451 ++++++++++++++++++++++++++++ moorpy/line.py | 25 +- moorpy/lineType.py | 6 +- moorpy/system.py | 10 +- 6 files changed, 601 insertions(+), 4 deletions(-) create mode 100644 moorpy/compositeLine.py create mode 100644 moorpy/dynamic_tension_functions.py diff --git a/moorpy/__init__.py b/moorpy/__init__.py index ef6e014..0f440ac 100644 --- a/moorpy/__init__.py +++ b/moorpy/__init__.py @@ -9,6 +9,8 @@ from moorpy.lineType import LineType from moorpy.system import System from moorpy.subsystem import Subsystem +#from moorpy.compositeLine import CompositeLine +#from moorpy.dynamic_tension_functions import * from moorpy.helpers import * from moorpy.Catenary import catenary diff --git a/moorpy/compositeLine.py b/moorpy/compositeLine.py new file mode 100644 index 0000000..bb45e31 --- /dev/null +++ b/moorpy/compositeLine.py @@ -0,0 +1,111 @@ +import numpy as np +from moorpy.dynamic_tension_functions import get_dynamic_tension, get_modes + +class CompositeLine(): + def __init__(self, sys, point_id, rho=1025.): + + self.rho = rho + point = sys.pointList[point_id-1] # Starting point id + + # check starting point to make sure it is either coupled or fixed and that it is not connected to more than 1 line. + if point.type == 0: + raise ValueError(f'Starting point ({point.number}) must not be free (change point type to -1 or 1).') + elif len(point.attached)>1: + raise ValueError(f'Starting point ({point.number}) cannot be connected to more than 1 line') + + self.pointList = [] + self.lineList = [] + + self.pointA = point + + # Move through the points along the composite + while True: + # make sure that the point is free + if len(point.attached)>2: + raise ValueError(f'Point {point.number} is attached to more than two lines.') + + # get the line id and line object + line_id = point.attached[-1] # get next line's id + line = sys.lineList[line_id - 1] # get line object + + # append the starting point and the line object + self.pointList.append(point) + self.lineList.append(line) + + # get the next point + attached_points = line.attached.copy() # get the IDs of the points attached to the line + pointA_id = point.number # get first point ID + attached_points.remove(pointA_id) # remove first point from attached point list + pointB_id = attached_points[0] # get second point id + point = sys.pointList[pointB_id-1] # get second point object + + # break from the loop when a point with a single attachment is reached + if len(point.attached) == 1: + self.pointList.append(point) + break + # make sure that the last point is not a free point + if point.type == 0: + raise ValueError(f'Last point ({point.number}) is a free point.') # make sure that the last point is not free + + self.pointB = point + self.nNodes = (np.sum([line.nNodes for line in self.lineList]) - (len(self.lineList)-1)) + self.sys = sys + + def getDynamicMatrices(self, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=1e-4): + n_nodes = self.nNodes # number of nodes + EA_segs = np.zeros(n_nodes-1) # extensional stiffness of the segments + n_dofs = 3*n_nodes # number of dofs + M = np.zeros([n_dofs,n_dofs], dtype='float') + A = np.zeros([n_dofs,n_dofs], dtype='float') + B = np.zeros([n_dofs,n_dofs], dtype='float') + K = np.zeros([n_dofs,n_dofs], dtype='float') + n = 0 + r_mean = np.zeros([n_nodes,3], dtype='float') + r_dynamic = np.ones((len(omegas),n_nodes,3),dtype='float')*r_dynamic + v_dynamic = 1j*omegas[:,None,None]*r_dynamic + + for line in self.lineList: + line_type = line.type + n1 = int(n/3) + n2 = n1 + line.nNodes + + # Filling matrices for line (dof n to dof 3xline_nodes+n) + M_n,A_n,B_n,K_n,r_n,EA_segs_n = line.getDynamicMatrices(omegas, S_zeta,r_dynamic[:,n1:n2,:],depth,kbot,cbot,seabed_tol=seabed_tol) + M[n:3*line.nNodes+n,n:3*line.nNodes+n] += M_n + A[n:3*line.nNodes+n,n:3*line.nNodes+n] += A_n + B[n:3*line.nNodes+n,n:3*line.nNodes+n] += B_n + K[n:3*line.nNodes+n,n:3*line.nNodes+n] += K_n + + # Attachment point properties + attachment = self.pointList[line.attached[-1]-1] # attachment point + attachment_idx = n2 - 1 # last node index + sigma_vp = np.sqrt(np.trapz(np.abs(v_dynamic[:,attachment_idx,:])**2*S_zeta[:,None],omegas,axis=0)) # standard deviations of the global components of the attachment point's velocity + + M[3*line.nNodes - 3:3*line.nNodes, 3*line.nNodes - 3:3*line.nNodes] += attachment.m*np.eye(3) + A[3*line.nNodes - 3:3*line.nNodes, 3*line.nNodes - 3:3*line.nNodes] += attachment.Ca* self.rho * attachment.v * np.eye(3) + B[3*line.nNodes - 3:3*line.nNodes, 3*line.nNodes - 3:3*line.nNodes] += 0.5* self.rho * attachment.CdA * np.pi*(3*attachment.v/4/np.pi)**(2/3) \ + * np.eye(3) * np.sqrt(8/np.pi) * np.diag(sigma_vp) + + # Static line properties + r_mean[n1:n2,:] = r_n + EA_segs[n1:n2-1] = EA_segs_n + + n += 3*line.nNodes - 3 # next line starting node add the number of dofs of the current line minus 3 to get the last shared node + + return M,A,B,K,r_mean,EA_segs + + def dynamicSolve(self,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_tol=1e-4,tol = 0.01,iters=100, w = 0.8): + T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X = get_dynamic_tension(self,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot, + seabed_tol=seabed_tol,tol = tol,iters=iters, w = w) + return T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X + + def getModes(self,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view = False,kbot=3E+06,cbot=3E+05,seabed_tol=1E-04): + + if plot_modes: + freqs,mode_shapes,r_nodes,M,A,K,fig,ax = get_modes(self,fix_A=fix_A,fix_B=fix_B,plot_modes=plot_modes,amp_factor=amp_factor,adj_view = adj_view, + kbot=kbot,cbot=cbot,seabed_tol=seabed_tol) + return freqs,mode_shapes,r_nodes,M,A,K,fig,ax + else: + freqs,mode_shapes,r_nodes,M,A,K = get_modes(self,fix_A=fix_A,fix_B=fix_B,plot_modes=plot_modes,amp_factor=amp_factor,adj_view = adj_view, + kbot=kbot,cbot=cbot,seabed_tol=seabed_tol) + return freqs,mode_shapes,r_nodes,M,A,K \ No newline at end of file diff --git a/moorpy/dynamic_tension_functions.py b/moorpy/dynamic_tension_functions.py new file mode 100644 index 0000000..22acdd5 --- /dev/null +++ b/moorpy/dynamic_tension_functions.py @@ -0,0 +1,451 @@ +import numpy as np +from numpy.linalg import solve,norm +import scipy.linalg as la +from datetime import datetime +import matplotlib.pyplot as plt +from collections.abc import Iterable + +def get_horizontal_oop_vec(p1,p2): + """ + Evaluates the horizontal out of plane vector given the coordinates of two points. + """ + hor_vec = p2 - np.array([p1[0],p1[1],p2[2]]) + ver_vec = p1 - np.array([p1[0],p1[1],p2[2]]) + + if np.isclose(np.linalg.norm(hor_vec),0): # vertical line + n_op = np.array([1,0,0]) + elif np.isclose(np.linalg.norm(ver_vec),0): # horizontal line + oop_vec = np.cross(hor_vec,np.array([0,0,1])) + n_op = oop_vec/np.linalg.norm(oop_vec) + else: + oop_vec = np.cross(hor_vec,ver_vec) + n_op = oop_vec/np.linalg.norm(oop_vec) + return n_op + +def get_dynamic_matrices(Line, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=1e-4): + """ + Evaluates dynamic matrices for a Line object. + + Parameters + ---------- + Line : Line + An instance of MoorPy's Line class + omegas : ndarray + Array of frequencies in rad/s. + S_zeta : ndarray + Wave spectrum array in m^2/(rad/s), must be of the same length as omegas. + r_dynamic : ndarray + A 3d array of the frequency dependent complex amplitudes of line nodes. The array has a shape of (m,n,3) where m is the number of frequencies, + n is the number of nodes, and the last dimension of 3 correspond to the x,y,z components. + depth : float + Water depth. + kbot : float + Vertical stiffness for points lying on the seabed. + cbot : float + Vertical damping for points lying on the seabed. + seabed_tol : float, optional + Distance from seabed within which a node is considered to be lying on the seabed, by default 1e-4 m. + + Returns + ------- + M: ndarray + Mass matrix. + A: ndarray + Added mass matrix. + B: ndarray + Damping matrix. + K: ndarray + Stiffness matrix. + M: ndarray + Mass matrix. + r_mean: ndarray + Mean static positions of the nodes given as a (m,3) array where m is the number of nodes. + EA_segs: ndarray + Extensional stiffness of segments. + """ + # extract line properties + N = Line.nNodes + mden = Line.type['m'] # line mass density function + deq = Line.type['d_vol'] # line volume equivalent diameter + EA = Line.type['EA'] # extensional stiffness + Ca = Line.type['Ca'] # normal added mass coeff + CaAx = Line.type['CaAx'] # tangential added mass coeff + Cd = Line.type['Cd'] # normal drag coeff + CdAx = Line.type['CdAx'] # tangential drag coeff + + # extract mean node coordinates + X_mean,Y_mean,Z_mean,T_mean = Line.getLineCoords(0.0,n=N) # coordinates of line nodes and tension values + r_mean = np.vstack((X_mean,Y_mean,Z_mean)).T # coordinates of line nodes + + # evaluate node velocities + r_dynamic = np.ones((len(omegas),N,3),dtype='float')*r_dynamic + v_dynamic = 1j*omegas[:,None,None]*r_dynamic + + # define out of plane normal + h_op = get_horizontal_oop_vec(r_mean[0],r_mean[-1]) # horizontal out-of-plane vector + hh_op = np.outer(h_op,h_op) + + # intialize line matrices + M = np.zeros([3*N, 3*N], dtype='float') # mass matrix + A = np.zeros([3*N, 3*N], dtype='float') # added mass matrix + B = np.zeros([3*N, 3*N], dtype='float') # linearized viscous damping matrix + K = np.zeros([3*N, 3*N], dtype='float') # stiffness matrix + + # Node 0 + dr_e1 = r_mean[1] - r_mean[0] + L_e1 = la.norm(dr_e1) # element 1 length + t_e1 = (dr_e1)/L_e1 # tangential unit vector + p_e1 = np.cross(t_e1,h_op) # in plane normal unit vector + + + ut_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],t_e1) # tangential velocity + uh_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],h_op) # normal horizontal out of plane velocity + up_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],p_e1) # normal in plane velocity + + sigma_ut_e1 = np.sqrt(np.trapz(np.abs(ut_e1)**2*S_zeta,omegas)) + sigma_uh_e1 = np.sqrt(np.trapz(np.abs(uh_e1)**2*S_zeta,omegas)) + sigma_up_e1 = np.sqrt(np.trapz(np.abs(up_e1)**2*S_zeta,omegas)) + + tt_e1 = np.outer(t_e1,t_e1) # local tangential to global components transformation matrix + pp_e1 = np.outer(p_e1,p_e1) # local normal inplane to global components transformation matrix + + M[0:3,0:3] += mden*L_e1/2*np.eye(3) # element 1 mass contribution + + A_e1 = 1025*np.pi/4*deq**2*L_e1/2*(Ca*(hh_op+pp_e1) + CaAx*tt_e1) # element 1 added mass contribution + + B_e1 = 0.5*1025*deq*L_e1/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_e1*hh_op + sigma_up_e1*pp_e1) + + CdAx*sigma_ut_e1*tt_e1) # element 1 damping contribution + + K_e1 = EA/L_e1*tt_e1 + (T_mean[0]/L_e1)*(hh_op+pp_e1) # element 1 stiffness (axial + geometric) + + ## assembling element 1 contributions (rows corresponding to node 0) + A[0:3,0:3] += A_e1 + B[0:3,0:3] += B_e1 + K[0:3,0:3] += K_e1 + K[0:3,3:6] += -K_e1 + + ## add seabed contribution to node 0 + if np.isclose(r_mean[0,2],-depth,seabed_tol): + K[2,2] += kbot + B[2,2] += cbot + + # Internal nodes loop (each internal node has contributions from two elements n-1/2 and n+1/2) + for n in range(1, N-1): + + ## backward element (n-1/2) contributions + dr_bw = r_mean[n-1] - r_mean[n] + L_bw = la.norm(dr_bw) # element 1 length + t_bw = (dr_bw)/L_bw # tangential unit vector + p_bw = np.cross(t_bw,h_op) # in plane normal unit vector + + ut_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],t_bw) # tangential velocity + uh_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],h_op) # normal horizontal out of plane velocity + up_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],p_bw) # normal in plane velocity + + sigma_ut_bw = np.sqrt(np.trapz(np.abs(ut_bw)**2*S_zeta,omegas)) + sigma_uh_bw = np.sqrt(np.trapz(np.abs(uh_bw)**2*S_zeta,omegas)) + sigma_up_bw = np.sqrt(np.trapz(np.abs(up_bw)**2*S_zeta,omegas)) + + tt_bw = np.outer(t_bw,t_bw) # local tangential to global components transformation matrix + pp_bw = np.outer(p_bw,p_bw) # local normal inplane to global components transformation matrix + + M[3*n:3*n+3,3*n:3*n+3] += mden*L_bw/2*np.eye(3) # mass contribution from adjacent elements + + A_bw = 1025*np.pi/4*deq**2*L_bw/2*(Ca*(hh_op+pp_bw) + CaAx*tt_bw) # backward element added mass contribution + + B_bw = 0.5*1025*deq*L_bw/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_bw*hh_op + sigma_up_bw*pp_bw) + + CdAx*sigma_ut_bw*tt_bw) # backward element damping contribution + + K_bw = EA/L_bw*tt_bw + (T_mean[n]/L_bw)*(hh_op+pp_bw) # backward element stiffness (axial + geometric) + + ## forward element (n+1/2) contributions + dr_fw = r_mean[n+1] - r_mean[n] + L_fw = la.norm(dr_fw) # element 1 length + t_fw = (dr_fw)/L_fw # tangential unit vector + p_fw = np.cross(t_fw,h_op) # in plane normal unit vector + + + ut_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],t_fw) # tangential velocity + uh_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],h_op) # normal horizontal out of plane velocity + up_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],p_fw) # normal in plane velocity + + sigma_ut_fw = np.sqrt(np.trapz(np.abs(ut_fw)**2*S_zeta,omegas)) + sigma_uh_fw = np.sqrt(np.trapz(np.abs(uh_fw)**2*S_zeta,omegas)) + sigma_up_fw = np.sqrt(np.trapz(np.abs(up_fw)**2*S_zeta,omegas)) + + tt_fw = np.outer(t_fw,t_fw) # local tangential to global components transformation matrix + pp_fw = np.outer(p_fw,p_fw) # local normal inplane to global components transformation matrix + + M[3*n:3*n+3,3*n:3*n+3] += mden*L_fw/2*np.eye(3) # mass contribution from adjacent elements + + A_fw = 1025*np.pi/4*deq**2*L_fw/2*(Ca*(hh_op+pp_fw) + CaAx*tt_fw) # forward element added mass contribution + + B_fw = 0.5*1025*deq*L_fw/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_fw*hh_op + sigma_up_fw*pp_fw) + + CdAx*sigma_ut_fw*tt_fw) # forward element damping contribution + + K_fw = EA/L_fw*tt_fw + (T_mean[n]/L_fw)*(hh_op+pp_fw) # forward element stiffness (axial + geometric) + + ## assembling bwd and fwd elements contributions (rows corresponding to node n) + A[3*n:3*n+3,3*n:3*n+3] += A_bw + A_fw + B[3*n:3*n+3,3*n:3*n+3] += B_bw + B_fw + K[3*n:3*n+3,3*n:3*n+3] += K_bw + K_fw + K[3*n:3*n+3,3*n-3:3*n] += -K_bw + K[3*n:3*n+3,3*n+3:3*n+6] += -K_fw + + ## add seabed contribution to node n + if np.isclose(r_mean[n,2],-depth,seabed_tol): + K[3*n+2,3*n+2] += kbot + B[3*n+2,3*n+2] += cbot + + # Node N + dr_eN = r_mean[N-1] - r_mean[N-2] + L_eN = la.norm(dr_eN) # element N length + t_eN = (dr_eN)/L_eN # tangential unit vector + p_eN = np.cross(t_eN,h_op) # in plane normal unit vector + + ut_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],t_eN) # tangential velocity + uh_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],h_op) # normal horizontal out of plane velocity + up_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],p_eN) # normal in plane velocity + + sigma_ut_eN = np.sqrt(np.trapz(np.abs(ut_eN)**2*S_zeta,omegas)) + sigma_uh_eN = np.sqrt(np.trapz(np.abs(uh_eN)**2*S_zeta,omegas)) + sigma_up_eN = np.sqrt(np.trapz(np.abs(up_eN)**2*S_zeta,omegas)) + + tt_eN = np.outer(t_eN,t_eN) # local tangential to global components transformation matrix + pp_eN = np.outer(p_eN,p_eN) # local normal inplane to global components transformation matrix + + M[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += mden*L_eN/2*np.eye(3) # element N mass contribution + + A_eN = 1025*np.pi/4*deq**2*L_eN/2*(Ca*(hh_op+pp_eN) + CaAx*tt_eN) # element N added mass contribution + + B_eN = 0.5*1025*deq*L_eN/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_eN*hh_op + sigma_up_eN*pp_eN) + + CdAx*sigma_ut_eN*tt_eN) # element N damping contribution + + K_eN = EA/L_eN*tt_eN + (T_mean[N-1]/L_eN)*(hh_op+pp_eN) # element N stiffness (axial + geometric) + + ## assembling element N contributions (rows corresponding to node N) + A[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += A_eN + B[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += B_eN + K[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += K_eN + K[3*(N-1):3*(N-1)+3,3*(N-1)-3:3*(N-1)] += -K_eN + + ## add seabed contribution to node N + if np.isclose(r_mean[N-1,2],-depth,seabed_tol): + K[3*(N-1)+2,3*(N-1)+2] += kbot + B[3*(N-1)+2,3*(N-1)+2] += cbot + + EA_segs = Line.type['EA']*np.ones(Line.nNodes - 1) + + return M,A,B,K,r_mean,EA_segs + +def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_tol=1e-4,tol = 0.01,iters=100, w = 0.8, conv_time=True): + """ + Evaluates dynamic tension at all the nodes for an instance of MoorPy's Line or CompositeLine classes. + + Parameters + ---------- + Line : Line/CompositeLine + An instance of MoorPy's Line or CompositeLine classes. + omegas : ndarray + Array of frequencies in rad/s. + S_zeta : ndarray + Wave spectrum array in m^2/(rad/s), must be of the same length as omegas. + RAO_A : ndarray + Translational RAOs for end node A given as a (m,3) array where m is the number of frequencies (must be equal to the length of omegas) . + RAO_B : ndarray + Translational RAOs for end node B given as a (m,3) array where m is the number of frequencies (must be equal to the length of omegas) . + depth : float + Water depth. + kbot : float + Vertical stiffness for points lying on the seabed. + cbot : float + Vertical damping for points lying on the seabed. + seabed_tol : float, optional + Distance from seabed within which a node is considered to be lying on the seabed, by default 1e-4 m. + tol : float, optional + Relative tolerance for iteration, by default 0.01 + iters : int, optional + Maximum number of iterations, by default 100 + w : float, optional + Succesive relaxation coefficient, by default 0.8 + + Returns + ------- + T_nodes_psd: ndarray + Tension spectra at nodes given as (m,n) array where m is the number of frequencies and n is the number of nodes. + T_nodes_std: ndarray + Tension standard deviation at nodes. + s: ndarray: + Node locations along the line. + r_static: ndarray + Nodes mean static position given as (n,3) array where n the number of nodes. + r_dynamic: ndarray + Nodes complex dynamic amplitudes given as (m,n,3) array where m the number of frequencies, n is the number of nodes + r_total: ndarray + Combined static and dynamic positions. + X: ndarray + Solution of the dynamic problem. + """ + N = Line.nNodes + n_dofs = 3*N + + if np.all(RAO_A == 0): + RAO_A = np.zeros_like(RAO_B) + if np.all(RAO_B == 0): + RAO_B = np.zeros_like(RAO_A) + + # intialize iteration matrices + r_dynamic_init = np.ones((len(omegas),N,3)) + M,A,B,K,r_static,EA_segs = Line.getDynamicMatrices(omegas,S_zeta,r_dynamic_init,depth,kbot,cbot,seabed_tol=seabed_tol) # TODO: return EA_segs + X = np.zeros((len(omegas),n_dofs),dtype = 'complex') + r_dynamic = np.zeros(((len(omegas),int(n_dofs/3),3)),dtype = 'complex') + S_Xd = np.zeros((len(omegas),n_dofs),dtype = 'float') + sigma_Xd = np.zeros(n_dofs,dtype = 'float') + sigma_Xd0 = np.zeros(n_dofs,dtype = 'float') + X[:, :3] = RAO_A + X[:,-3:] = RAO_B + + # solving dynamics + start = datetime.now() + for ni in range(iters): + H = - omegas[:,None,None]**2*(M+A)[None,:,:]\ + + 1j*omegas[:,None,None]*B[None,:,:]\ + + K[None,:,:]\ + + F_A = np.einsum('nij,njk->ni',-H[:,3:-3, :3],RAO_A[:,:,None]) + F_B = np.einsum('nij,njk->ni',-H[:,3:-3,-3:],RAO_B[:,:,None]) + F = F_A + F_B + + X[:,3:-3] = solve(H[:,3:-3,3:-3],F) + S_Xd[:] = np.abs(1j*omegas[:,None]*X)**2*S_zeta[:,None] + sigma_Xd[:] = np.sqrt(np.trapz(S_Xd,omegas,axis=0)) + r_dynamic[:] = X.reshape(X.shape[0],int(X.shape[1]/3),3) + if (np.abs(sigma_Xd-sigma_Xd0) <= tol*np.abs(sigma_Xd0)).all(): + break + else: + sigma_Xd0[:] = w * sigma_Xd + (1.-w) * sigma_Xd0 + _,_,B[:],_,_,_ = Line.getDynamicMatrices(omegas,S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=seabed_tol) + if conv_time: + print(f'Finished {ni} dynamic tension iterations in {datetime.now()-start} seconds (w = {w}).') + + # evaluate tension + dw = np.diff(omegas, + prepend= omegas[0] - (omegas[1]-omegas[0]), + append= omegas[-1] + (omegas[-1]-omegas[-2])) + dw = (dw[1:]+dw[:-1])/2 + wave_amps = np.sqrt(S_zeta*dw) #evaluate wave amplitudes of harmonic components from wave spectrum + + r_dynamic *= wave_amps[:,None,None] + r_total = r_static[None,:,:] + r_dynamic + dr_static = r_static[:-1] - r_static[1:] + dr_dynamic = r_dynamic[:,:-1,:] - r_dynamic[:,1:,:] + tangents = dr_static/la.norm(r_static[:-1] - r_static[1:], axis=-1)[:,None] + L_static = la.norm(dr_static, axis=-1) + dL_dynamic = np.einsum('mni,ni->mn', dr_dynamic, tangents) + eps_segs = np.abs(dL_dynamic)/L_static + + T_segs = EA_segs * eps_segs + T_nodes = np.zeros((len(omegas),N)) + T_nodes[:,0] = T_segs[:,0] + T_nodes[:,1:-1] = (T_segs[:,:-1] + T_segs[:,1:])/2 + T_nodes[:,-1] = T_segs[:,-1] + + # S_T = np.zeros((len(omegas),N)) + # S_T[:,1:] = T_e**2/dw[:,None] + # S_T[:,0] = S_T[:,1] + + T_nodes_psd = T_nodes**2/dw[:,None] + T_nodes_std = np.sqrt(np.trapz(T_nodes_psd,omegas,axis=0)) + + + dr = np.diff(r_static,axis=0) + ds = la.norm(dr,axis=1) + s = np.zeros_like(T_nodes_std) + s = np.cumsum(ds) + + return T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X + +def get_modes(line,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view = False,kbot=3E+06,cbot=3E+05,seabed_tol=1E-04): + M,A,_,K,r_nodes,_ = line.getDynamicMatrices(np.ones(1), np.ones(1),0.,line.sys.depth,kbot,cbot,seabed_tol=seabed_tol) + + if fix_A: + n1 = 1 + else: + n1 = 0 + + if fix_B: + n2 = -1 + else: + n2 = r_nodes.shape[0] + + eigvals,eigvecs = la.eig(K[3*n1:3*n2,3*n1:3*n2],M[3*n1:3*n2,3*n1:3*n2]+A[3*n1:3*n2,3*n1:3*n2]) + stable_eigvals = eigvals[np.real(eigvals)>0] + stable_eigvecs = eigvecs[:,np.real(eigvals)>0] + + idx = stable_eigvals.argsort()[::-1] + stable_eigvals = stable_eigvals[idx] + stable_eigvecs = stable_eigvecs[:,idx] + + freqs = np.sqrt(np.real(stable_eigvals))/2/np.pi + mode_shapes = np.zeros(stable_eigvecs.shape,dtype='float') + + for i in range(stable_eigvecs.shape[1]): + mode_shapes[:,i] = r_nodes[n1:n2].flatten('C') + stable_eigvecs[:,i] + + freqs = np.flip(freqs) + mode_shapes = np.flip(mode_shapes,axis=1) + + if plot_modes: + cols = 4 + rows = plot_modes//cols + bool(plot_modes%cols) + fig,ax = plt.subplots(rows,cols,subplot_kw={"projection": "3d"},figsize=(5*cols,5*rows)) + + i = 0 + for axes in ax: + if not isinstance(axes,Iterable): + axes = [axes] + + for axis in axes: + if i >= plot_modes: + break + + r = r_nodes.copy() + r[n1:n2] = mode_shapes[:,i].reshape([int(len(mode_shapes[:,i])/3),3]) + r = (r-r_nodes)*amp_factor + x = r[:,0] + y = r[:,1] + z = r[:,2] + + x_0 = r_nodes[:,0] + y_0 = r_nodes[:,1] + z_0 = r_nodes[:,2] + + axis.plot(x_0,y_0,z_0,'-ko',label='initial') + axis.plot(x+x_0,y+y_0,z+z_0,'--ro',label='mode shape') + + # R_0 = np.sqrt(x_0**2 + y_0**2) + if adj_view: + # h_min = np.min((x_0,y_0)) + # h_max = np.max((x_0,y_0)) + # axis.set_xlim(h_min,h_max) + # axis.set_ylim(h_min,h_max) + # axis.set_zlim(z_0.min(),z_0.max()) + sigma_x = x.std() + sigma_y = y.std() + sigma_z = z.std() + azim = np.arctan2(sigma_x,sigma_y)*180/np.pi + elev = np.arctan2(np.hypot(sigma_x,sigma_y),sigma_z)*180/np.pi + axis.view_init(elev=elev,azim=azim) + + # axis.set_box_aspect([np.ptp(coord) for coord in [x, y, z]]) + axis.set_xlabel('X (m)') + axis.set_ylabel('Y (m)') + axis.set_zlabel('Z (m)') + axis.set_title(f'f = {freqs[i]:.3f} Hz, T = {1/freqs[i]:.3f} s') + + i+=1 + + fig.tight_layout() + return freqs,mode_shapes,r_nodes,M,A,K,fig,ax + else: + return freqs,mode_shapes,r_nodes,M,A,K diff --git a/moorpy/line.py b/moorpy/line.py index ed97576..379603c 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -7,6 +7,7 @@ from moorpy.helpers import (unitVector, LineError, CatenaryError, rotationMatrix, makeTower, read_mooring_file, quiver_data_to_segments, printVec, printMat) +from moorpy.dynamic_tension_functions import get_dynamic_matrices,get_dynamic_tension,get_modes from os import path @@ -41,6 +42,8 @@ def __init__(self, mooringSys, num, L, lineType, nSegs=100, cb=0, isRod=0, attac self.sys = mooringSys # store a reference to the overall mooring system (instance of System class) + self.attached = attachments # ID numbers of the Points at the Line ends [a,b] >>> NOTE: not fully supported <<<< + self.number = num self.isRod = isRod @@ -75,7 +78,7 @@ def __init__(self, mooringSys, num, L, lineType, nSegs=100, cb=0, isRod=0, attac self.lw=0.5 self.fCurrent = np.zeros(3) # total current force vector on the line [N] - + def loadData(self, dirname, rootname, sep='.MD.', id=0): '''Loads line-specific time series data from a MoorDyn output file''' @@ -1001,6 +1004,26 @@ def revertToStaticStiffness(self): # revert to original line length self.L = self.L0 + def getDynamicMatrices(self, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=1e-4): + M,A,B,K,r_mean,EA_segs = get_dynamic_matrices(self, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=seabed_tol) + return M,A,B,K,r_mean,EA_segs + + def dynamicSolve(self,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_tol=1e-4,tol = 0.01,iters=100, w = 0.8): + T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X = get_dynamic_tension(self,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot, + seabed_tol=seabed_tol,tol = tol,iters=iters, w = w) + return T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X + + def getModes(self,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view = False,kbot=3E+06,cbot=3E+05,seabed_tol=1E-04): + + if plot_modes: + freqs,mode_shapes,r_nodes,M,A,K,fig,ax = get_modes(self,fix_A=fix_A,fix_B=fix_B,plot_modes=plot_modes,amp_factor=amp_factor,adj_view = adj_view, + kbot=kbot,cbot=cbot,seabed_tol=seabed_tol) + return freqs,mode_shapes,r_nodes,M,A,K,fig,ax + else: + freqs,mode_shapes,r_nodes,M,A,K = get_modes(self,fix_A=fix_A,fix_B=fix_B,plot_modes=plot_modes,amp_factor=amp_factor,adj_view = adj_view, + kbot=kbot,cbot=cbot,seabed_tol=seabed_tol) + return freqs,mode_shapes,r_nodes,M,A,K + def from2Dto3Drotated(K2D, F, L, R): '''Initialize a line end's analytic stiffness matrix in the diff --git a/moorpy/lineType.py b/moorpy/lineType.py index 03b5c7e..186def4 100644 --- a/moorpy/lineType.py +++ b/moorpy/lineType.py @@ -7,7 +7,7 @@ class LineType(): '''A class to hold the various properties of a mooring line type''' - def __init__(self, name, d, massden, EA, MBL=0.0, cost=0.0, notes="", input_d=0.0, input_type=''): + def __init__(self, name, d, massden, EA, MBL=0.0, cost=0.0, notes="", input_d=0.0, input_type='',Cd=0.,Ca=0.,CdAx=0.,CaAx=0.): '''Initialize LineType attributes Parameters @@ -45,6 +45,10 @@ def __init__(self, name, d, massden, EA, MBL=0.0, cost=0.0, notes="", input_d=0. self.notes = notes # optional notes/description string self.input_d = input_d # the non-volume-equivalent, input diameter [m] self.input_type = input_type # line type string (e.g. chain, polyester) + self.Cd = Cd + self.Ca = Ca + self.CdAx = CdAx + self.CaAx = CaAx diff --git a/moorpy/system.py b/moorpy/system.py index d555854..4bf01a6 100644 --- a/moorpy/system.py +++ b/moorpy/system.py @@ -18,6 +18,7 @@ from moorpy.point import Point from moorpy.line import Line from moorpy.lineType import LineType +from moorpy.compositeLine import CompositeLine import matplotlib as mpl #import moorpy.MoorSolve as msolve from moorpy.helpers import (rotationMatrix, rotatePosition, getH, printVec, @@ -271,7 +272,7 @@ def addLine(self, lUnstr, lineType, nSegs=40, pointA=0, pointB=0, cb=0): else: raise ValueError(f"The specified lineType name ({lineType}) does not correspond with any lineType stored in this MoorPy System") - self.lineList.append( Line(self, len(self.lineList)+1, lUnstr, lineType, nSegs=nSegs, cb=cb) ) + self.lineList.append( Line(self, len(self.lineList)+1, lUnstr, lineType, nSegs=nSegs, cb=cb, attachments = [pointA,pointB]) ) if pointA > 0: if pointA <= len(self.pointList): @@ -377,6 +378,11 @@ def addLineType(self, type_string, d, mass, EA, name=""): # <<< the "name" keyword in this method is confusing in that it isn't the index key. Does it have a purpose? <<< + def captureCompositeLine(self,point_id): + if not hasattr(self,'compositeLineList'): + self.compositeLineList = [] + + self.compositeLineList.append(CompositeLine(self,point_id)) def setLineType(self, dnommm, material, source=None, name="", **kwargs): '''Add or update a System lineType using the new dictionary-based method. @@ -767,7 +773,7 @@ def load(self, filename, clear=True): nSegs = np.int_(entries[5]) #lineList.append( Line(dirName, num, lUnstr, dia, nSegs) ) - self.lineList.append( Line(self, num, lUnstr, lineType, nSegs=nSegs)) #attachments = [int(entries[4]), int(entries[5])]) ) + self.lineList.append( Line(self, num, lUnstr, lineType, nSegs=nSegs, attachments = [int(entries[2]), int(entries[3])]) ) # attach end A numA = int("".join(filter(str.isdigit, entries[2]))) # get number from the attachA string From 14fc51d95557336969dc344e8bee0c88c209ddfb Mon Sep 17 00:00:00 2001 From: Matt Hall <5151457+mattEhall@users.noreply.github.com> Date: Tue, 6 Feb 2024 23:54:11 -0700 Subject: [PATCH 02/14] New dymamic functions from Serag now working with Subsystem: - Edited new Line dynamic methods to streamline passing of arguments and to only load the new functions on demand. - Added Subsystem.getDynamicMatrices method - copied closely from the CompositeLine method. This allows the new methods to work with Subsystem (all other methods are already inherited from Line). - Added Subsystem.makeFromSystem method - adapted from CompositeLine initialization to allow similar functionality for Subsystem. Not tested. - Made Examples/dynamic.py script to test some things out simply. --- examples/dynamic.py | 129 ++++++++++++++++++++++++++ moorpy/dynamic_tension_functions.py | 6 +- moorpy/line.py | 38 ++++---- moorpy/subsystem.py | 137 +++++++++++++++++++++++++++- 4 files changed, 290 insertions(+), 20 deletions(-) create mode 100644 examples/dynamic.py diff --git a/examples/dynamic.py b/examples/dynamic.py new file mode 100644 index 0000000..b78a99e --- /dev/null +++ b/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() diff --git a/moorpy/dynamic_tension_functions.py b/moorpy/dynamic_tension_functions.py index 22acdd5..720743c 100644 --- a/moorpy/dynamic_tension_functions.py +++ b/moorpy/dynamic_tension_functions.py @@ -5,6 +5,7 @@ import matplotlib.pyplot as plt from collections.abc import Iterable + def get_horizontal_oop_vec(p1,p2): """ Evaluates the horizontal out of plane vector given the coordinates of two points. @@ -22,6 +23,7 @@ def get_horizontal_oop_vec(p1,p2): n_op = oop_vec/np.linalg.norm(oop_vec) return n_op + def get_dynamic_matrices(Line, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=1e-4): """ Evaluates dynamic matrices for a Line object. @@ -238,6 +240,7 @@ def get_dynamic_matrices(Line, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_t return M,A,B,K,r_mean,EA_segs + def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_tol=1e-4,tol = 0.01,iters=100, w = 0.8, conv_time=True): """ Evaluates dynamic tension at all the nodes for an instance of MoorPy's Line or CompositeLine classes. @@ -365,6 +368,7 @@ def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_to return T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X + def get_modes(line,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view = False,kbot=3E+06,cbot=3E+05,seabed_tol=1E-04): M,A,_,K,r_nodes,_ = line.getDynamicMatrices(np.ones(1), np.ones(1),0.,line.sys.depth,kbot,cbot,seabed_tol=seabed_tol) @@ -377,7 +381,7 @@ def get_modes(line,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view n2 = -1 else: n2 = r_nodes.shape[0] - + eigvals,eigvecs = la.eig(K[3*n1:3*n2,3*n1:3*n2],M[3*n1:3*n2,3*n1:3*n2]+A[3*n1:3*n2,3*n1:3*n2]) stable_eigvals = eigvals[np.real(eigvals)>0] stable_eigvecs = eigvecs[:,np.real(eigvals)>0] diff --git a/moorpy/line.py b/moorpy/line.py index 379603c..2e2457b 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -7,7 +7,7 @@ from moorpy.helpers import (unitVector, LineError, CatenaryError, rotationMatrix, makeTower, read_mooring_file, quiver_data_to_segments, printVec, printMat) -from moorpy.dynamic_tension_functions import get_dynamic_matrices,get_dynamic_tension,get_modes + from os import path @@ -1004,27 +1004,29 @@ def revertToStaticStiffness(self): # revert to original line length self.L = self.L0 - def getDynamicMatrices(self, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=1e-4): - M,A,B,K,r_mean,EA_segs = get_dynamic_matrices(self, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=seabed_tol) - return M,A,B,K,r_mean,EA_segs - - def dynamicSolve(self,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_tol=1e-4,tol = 0.01,iters=100, w = 0.8): - T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X = get_dynamic_tension(self,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot, - seabed_tol=seabed_tol,tol = tol,iters=iters, w = w) - return T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X - def getModes(self,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view = False,kbot=3E+06,cbot=3E+05,seabed_tol=1E-04): + def getDynamicMatrices(self, *args, **kwargs): + '''Compute M,A,B,K matrices for Line. See get_dynamic_matrices().''' + from moorpy.dynamic_tension_functions import get_dynamic_matrices - if plot_modes: - freqs,mode_shapes,r_nodes,M,A,K,fig,ax = get_modes(self,fix_A=fix_A,fix_B=fix_B,plot_modes=plot_modes,amp_factor=amp_factor,adj_view = adj_view, - kbot=kbot,cbot=cbot,seabed_tol=seabed_tol) - return freqs,mode_shapes,r_nodes,M,A,K,fig,ax - else: - freqs,mode_shapes,r_nodes,M,A,K = get_modes(self,fix_A=fix_A,fix_B=fix_B,plot_modes=plot_modes,amp_factor=amp_factor,adj_view = adj_view, - kbot=kbot,cbot=cbot,seabed_tol=seabed_tol) - return freqs,mode_shapes,r_nodes,M,A,K + return get_dynamic_matrices(self, *args, **kwargs) + def dynamicSolve(self, *args, **kwargs): + '''Compute complex amplitudes of line nodes. See get_dynamic_tension().''' + from moorpy.dynamic_tension_functions import get_dynamic_tension + + return get_dynamic_tension(self, *args, **kwargs) + + + def getModes(self,*args, **kwargs): + '''Compute (and optionally plot) the line's mode shapes. + See get_modes().''' + from moorpy.dynamic_tension_functions import get_modes + + return get_modes(self, *args, **kwargs) + + def from2Dto3Drotated(K2D, F, L, R): '''Initialize a line end's analytic stiffness matrix in the plane of the catenary then rotate the matrix to be about the diff --git a/moorpy/subsystem.py b/moorpy/subsystem.py index bd0f1fa..f708f44 100644 --- a/moorpy/subsystem.py +++ b/moorpy/subsystem.py @@ -160,7 +160,7 @@ def makeGeneric(self, lengths, types, points=[], shared=False): if shared: self.addPoint(-1, rA, DOFs=[2]) # add shared line point, free only to move in z else: - self.addPoint(-1, rA, DOFs=[0,2]) # add anchor point + self.addPoint(-1, rA, DOFs=[0,2]) # add anchor point # Go through each line segment and add its upper point, add the line, and connect the line to the points for i in range(self.nLines): @@ -194,6 +194,88 @@ def makeGeneric(self, lengths, types, points=[], shared=False): # if a points list is provided, apply any mass or other properties it contains? + self.nLines = len(self.lineList) + self.nNodes = np.sum([line.nNodes for line in self.lineList]) - self.nLines + 1 + + + def makeFromSystem(self, sys, point_id): + '''Populate a Subsystem based on a series of mooring lines in an + existing system (sys), starting at an end point (point_id). + Taken from Serag's CompositeLine class. + In this version, the subsystem still needs to be added to a system + afterward. + ''' + + if len(self.pointList)+len(self.lineList) > 0: + raise Exception('makeFromSystem can only be called for an empty Subsystem.') + + point = sys.pointList[point_id-1] # Starting point id + + # check starting point to make sure it is either coupled or fixed and that it is not connected to more than 1 line. + if point.type == 0: + raise Exception(f'Starting point ({point.number}) must not be free (change point type to -1 or 1).') + elif len(point.attached)>1: + raise Exception(f'Starting point ({point.number}) cannot be connected to more than 1 line') + + # Save point A info + self.rA = np.array(point.r) + self.addPoint(1, self.rA) # add anchor point + + # Move through the points along the composite + while True: + # make sure that the point is free + if len(point.attached) > 2: + raise Exception(f'Point {point.number} is attached to more than two lines.') + + # get the line id and line object + line_id = point.attached[-1] # get next line's id + line = sys.lineList[line_id - 1] # get line object + self.lineTypes[line.type['name']] = dict(line.type) # copy the lineType dict + self.addLine(line.L0, self.lineTypes[line.type['name']]) # add the line + + + # get the next point + attached_points = line.attached.copy() # get the IDs of the points attached to the line + pointA_id = point.number # get first point ID + attached_points.remove(pointA_id) # remove first point from attached point list + pointB_id = attached_points[0] # get second point id + point = sys.pointList[pointB_id-1] # get second point object + + if line(point.attached) == 1: # must be the endpoint + self.addPoint(-1, point.r) + else: # intermediate point along line + self.addPoint(0, point.r, DOFs=[0,2]) # may need to ensure there's no y component + + # Attach the line ends to the points + self.pointList[-2].attachLine(len(self.lineList), 0) + self.pointList[-1].attachLine(len(self.lineList), 1) + + # break from the loop when a point with a single attachment is reached + if len(point.attached) == 1: + break + + # make sure that the last point is not a free point + if point.type == 0: + raise Exception(f'Last point ({point.number}) is a free point.') + + self.rB = np.array(point.r) + self.nLines = len(self.lineList) + self.nNodes = np.sum([line.nNodes for line in self.lineList]) - self.nLines + 1 + + + def initialize(self): + '''Initializes the subsystem objects to their initial positions, and + counts number of nodes.''' + + self.nDOF, self.nCpldDOF, _ = self.getDOFs() + + self.nNodes = np.sum([line.nNodes for line in self.lineList]) - self.nLines + 1 + + for point in self.pointList: + point.setPosition(point.r) + + self.staticSolve() + def setEndPosition(self, r, endB, sink=False): '''Sets either end position of the subsystem in the global/system @@ -395,6 +477,59 @@ def revertToStaticStiffness(self): System.revertToStaticStiffness(self) + # ----- Function for dynamic frequency-domain tensions ----- + + def getDynamicMatrices(self, omegas, S_zeta, r_dynamic, depth, kbot, cbot, + seabed_tol=1e-4): + '''Compute M,A,B,K matrices for the Subsystem. This calls + get_dynamic_matrices() for each Line in the Subsystem then combines + the results. Note that this method overrides the Line method. Other + Line methods used for dynamics can be used directly in Subsystem. + ''' + self.nNodes = np.sum([line.nNodes for line in self.lineList]) - self.nLines + 1 + + EA_segs = np.zeros(self.nNodes-1) # extensional stiffness of the segments + n_dofs = 3*self.nNodes # number of dofs + M = np.zeros([n_dofs,n_dofs], dtype='float') + A = np.zeros([n_dofs,n_dofs], dtype='float') + B = np.zeros([n_dofs,n_dofs], dtype='float') + K = np.zeros([n_dofs,n_dofs], dtype='float') + r_mean = np.zeros([self.nNodes,3], dtype='float') + r_dynamic = np.ones((len(omegas),self.nNodes,3),dtype='float')*r_dynamic + v_dynamic = 1j*omegas[:,None,None]*r_dynamic + + n = 0 # starting index of the next line's entries in the matrices + + for line in self.lineList: + n1 = int(n/3) + n2 = n1 + line.nNodes + + # Filling matrices for line (dof n to dof 3xline_nodes+n) + M_n,A_n,B_n,K_n,r_n,EA_segs_n = line.getDynamicMatrices(omegas, S_zeta,r_dynamic[:,n1:n2,:],depth,kbot,cbot,seabed_tol=seabed_tol) + M[n:3*line.nNodes+n,n:3*line.nNodes+n] += M_n + A[n:3*line.nNodes+n,n:3*line.nNodes+n] += A_n + B[n:3*line.nNodes+n,n:3*line.nNodes+n] += B_n + K[n:3*line.nNodes+n,n:3*line.nNodes+n] += K_n + + # Attachment point properties + attachment = self.pointList[line.attached[-1]-1] # attachment point + attachment_idx = n2 - 1 # last node index + sigma_vp = np.sqrt(np.trapz(np.abs(v_dynamic[:,attachment_idx,:])**2*S_zeta[:,None],omegas,axis=0)) # standard deviations of the global components of the attachment point's velocity + + M[3*line.nNodes - 3:3*line.nNodes, 3*line.nNodes - 3:3*line.nNodes] += attachment.m*np.eye(3) + A[3*line.nNodes - 3:3*line.nNodes, 3*line.nNodes - 3:3*line.nNodes] += attachment.Ca* self.rho * attachment.v * np.eye(3) + B[3*line.nNodes - 3:3*line.nNodes, 3*line.nNodes - 3:3*line.nNodes] += 0.5* self.rho * attachment.CdA * np.pi*(3*attachment.v/4/np.pi)**(2/3) \ + * np.eye(3) * np.sqrt(8/np.pi) * np.diag(sigma_vp) + + # Static line properties + r_mean[n1:n2,:] = r_n + EA_segs[n1:n2-1] = EA_segs_n + + n += 3*line.nNodes - 3 # next line starting node add the number of dofs of the current line minus 3 to get the last shared node + + return M,A,B,K,r_mean,EA_segs + + # ---- Extra convenience functions (subsystem should be in equilibrium) ----- From 663b66d7ac217904fa4fd5d79e43428e6f61263e Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Wed, 22 May 2024 13:33:48 -0600 Subject: [PATCH 03/14] Moved dynamic tension and compositeLine functions to existing files - functions from dynamic_tension_functions.py are now on helpers - Because the dynamic calculations depend on the number of nodes in the subsystem/line, I slightly modified functions lines2subsystem and makeGeneric to specify the number of segments to be used in addLine. - I have not used subsystem.makeFromSystem or getNodes yet --- moorpy/compositeLine.py | 111 ------- moorpy/dynamic_tension_functions.py | 455 --------------------------- moorpy/helpers.py | 457 +++++++++++++++++++++++++++- moorpy/line.py | 9 +- moorpy/subsystem.py | 16 +- moorpy/system.py | 1 - 6 files changed, 469 insertions(+), 580 deletions(-) delete mode 100644 moorpy/compositeLine.py delete mode 100644 moorpy/dynamic_tension_functions.py diff --git a/moorpy/compositeLine.py b/moorpy/compositeLine.py deleted file mode 100644 index bb45e31..0000000 --- a/moorpy/compositeLine.py +++ /dev/null @@ -1,111 +0,0 @@ -import numpy as np -from moorpy.dynamic_tension_functions import get_dynamic_tension, get_modes - -class CompositeLine(): - def __init__(self, sys, point_id, rho=1025.): - - self.rho = rho - point = sys.pointList[point_id-1] # Starting point id - - # check starting point to make sure it is either coupled or fixed and that it is not connected to more than 1 line. - if point.type == 0: - raise ValueError(f'Starting point ({point.number}) must not be free (change point type to -1 or 1).') - elif len(point.attached)>1: - raise ValueError(f'Starting point ({point.number}) cannot be connected to more than 1 line') - - self.pointList = [] - self.lineList = [] - - self.pointA = point - - # Move through the points along the composite - while True: - # make sure that the point is free - if len(point.attached)>2: - raise ValueError(f'Point {point.number} is attached to more than two lines.') - - # get the line id and line object - line_id = point.attached[-1] # get next line's id - line = sys.lineList[line_id - 1] # get line object - - # append the starting point and the line object - self.pointList.append(point) - self.lineList.append(line) - - # get the next point - attached_points = line.attached.copy() # get the IDs of the points attached to the line - pointA_id = point.number # get first point ID - attached_points.remove(pointA_id) # remove first point from attached point list - pointB_id = attached_points[0] # get second point id - point = sys.pointList[pointB_id-1] # get second point object - - # break from the loop when a point with a single attachment is reached - if len(point.attached) == 1: - self.pointList.append(point) - break - # make sure that the last point is not a free point - if point.type == 0: - raise ValueError(f'Last point ({point.number}) is a free point.') # make sure that the last point is not free - - self.pointB = point - self.nNodes = (np.sum([line.nNodes for line in self.lineList]) - (len(self.lineList)-1)) - self.sys = sys - - def getDynamicMatrices(self, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=1e-4): - n_nodes = self.nNodes # number of nodes - EA_segs = np.zeros(n_nodes-1) # extensional stiffness of the segments - n_dofs = 3*n_nodes # number of dofs - M = np.zeros([n_dofs,n_dofs], dtype='float') - A = np.zeros([n_dofs,n_dofs], dtype='float') - B = np.zeros([n_dofs,n_dofs], dtype='float') - K = np.zeros([n_dofs,n_dofs], dtype='float') - n = 0 - r_mean = np.zeros([n_nodes,3], dtype='float') - r_dynamic = np.ones((len(omegas),n_nodes,3),dtype='float')*r_dynamic - v_dynamic = 1j*omegas[:,None,None]*r_dynamic - - for line in self.lineList: - line_type = line.type - n1 = int(n/3) - n2 = n1 + line.nNodes - - # Filling matrices for line (dof n to dof 3xline_nodes+n) - M_n,A_n,B_n,K_n,r_n,EA_segs_n = line.getDynamicMatrices(omegas, S_zeta,r_dynamic[:,n1:n2,:],depth,kbot,cbot,seabed_tol=seabed_tol) - M[n:3*line.nNodes+n,n:3*line.nNodes+n] += M_n - A[n:3*line.nNodes+n,n:3*line.nNodes+n] += A_n - B[n:3*line.nNodes+n,n:3*line.nNodes+n] += B_n - K[n:3*line.nNodes+n,n:3*line.nNodes+n] += K_n - - # Attachment point properties - attachment = self.pointList[line.attached[-1]-1] # attachment point - attachment_idx = n2 - 1 # last node index - sigma_vp = np.sqrt(np.trapz(np.abs(v_dynamic[:,attachment_idx,:])**2*S_zeta[:,None],omegas,axis=0)) # standard deviations of the global components of the attachment point's velocity - - M[3*line.nNodes - 3:3*line.nNodes, 3*line.nNodes - 3:3*line.nNodes] += attachment.m*np.eye(3) - A[3*line.nNodes - 3:3*line.nNodes, 3*line.nNodes - 3:3*line.nNodes] += attachment.Ca* self.rho * attachment.v * np.eye(3) - B[3*line.nNodes - 3:3*line.nNodes, 3*line.nNodes - 3:3*line.nNodes] += 0.5* self.rho * attachment.CdA * np.pi*(3*attachment.v/4/np.pi)**(2/3) \ - * np.eye(3) * np.sqrt(8/np.pi) * np.diag(sigma_vp) - - # Static line properties - r_mean[n1:n2,:] = r_n - EA_segs[n1:n2-1] = EA_segs_n - - n += 3*line.nNodes - 3 # next line starting node add the number of dofs of the current line minus 3 to get the last shared node - - return M,A,B,K,r_mean,EA_segs - - def dynamicSolve(self,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_tol=1e-4,tol = 0.01,iters=100, w = 0.8): - T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X = get_dynamic_tension(self,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot, - seabed_tol=seabed_tol,tol = tol,iters=iters, w = w) - return T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X - - def getModes(self,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view = False,kbot=3E+06,cbot=3E+05,seabed_tol=1E-04): - - if plot_modes: - freqs,mode_shapes,r_nodes,M,A,K,fig,ax = get_modes(self,fix_A=fix_A,fix_B=fix_B,plot_modes=plot_modes,amp_factor=amp_factor,adj_view = adj_view, - kbot=kbot,cbot=cbot,seabed_tol=seabed_tol) - return freqs,mode_shapes,r_nodes,M,A,K,fig,ax - else: - freqs,mode_shapes,r_nodes,M,A,K = get_modes(self,fix_A=fix_A,fix_B=fix_B,plot_modes=plot_modes,amp_factor=amp_factor,adj_view = adj_view, - kbot=kbot,cbot=cbot,seabed_tol=seabed_tol) - return freqs,mode_shapes,r_nodes,M,A,K \ No newline at end of file diff --git a/moorpy/dynamic_tension_functions.py b/moorpy/dynamic_tension_functions.py deleted file mode 100644 index 720743c..0000000 --- a/moorpy/dynamic_tension_functions.py +++ /dev/null @@ -1,455 +0,0 @@ -import numpy as np -from numpy.linalg import solve,norm -import scipy.linalg as la -from datetime import datetime -import matplotlib.pyplot as plt -from collections.abc import Iterable - - -def get_horizontal_oop_vec(p1,p2): - """ - Evaluates the horizontal out of plane vector given the coordinates of two points. - """ - hor_vec = p2 - np.array([p1[0],p1[1],p2[2]]) - ver_vec = p1 - np.array([p1[0],p1[1],p2[2]]) - - if np.isclose(np.linalg.norm(hor_vec),0): # vertical line - n_op = np.array([1,0,0]) - elif np.isclose(np.linalg.norm(ver_vec),0): # horizontal line - oop_vec = np.cross(hor_vec,np.array([0,0,1])) - n_op = oop_vec/np.linalg.norm(oop_vec) - else: - oop_vec = np.cross(hor_vec,ver_vec) - n_op = oop_vec/np.linalg.norm(oop_vec) - return n_op - - -def get_dynamic_matrices(Line, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=1e-4): - """ - Evaluates dynamic matrices for a Line object. - - Parameters - ---------- - Line : Line - An instance of MoorPy's Line class - omegas : ndarray - Array of frequencies in rad/s. - S_zeta : ndarray - Wave spectrum array in m^2/(rad/s), must be of the same length as omegas. - r_dynamic : ndarray - A 3d array of the frequency dependent complex amplitudes of line nodes. The array has a shape of (m,n,3) where m is the number of frequencies, - n is the number of nodes, and the last dimension of 3 correspond to the x,y,z components. - depth : float - Water depth. - kbot : float - Vertical stiffness for points lying on the seabed. - cbot : float - Vertical damping for points lying on the seabed. - seabed_tol : float, optional - Distance from seabed within which a node is considered to be lying on the seabed, by default 1e-4 m. - - Returns - ------- - M: ndarray - Mass matrix. - A: ndarray - Added mass matrix. - B: ndarray - Damping matrix. - K: ndarray - Stiffness matrix. - M: ndarray - Mass matrix. - r_mean: ndarray - Mean static positions of the nodes given as a (m,3) array where m is the number of nodes. - EA_segs: ndarray - Extensional stiffness of segments. - """ - # extract line properties - N = Line.nNodes - mden = Line.type['m'] # line mass density function - deq = Line.type['d_vol'] # line volume equivalent diameter - EA = Line.type['EA'] # extensional stiffness - Ca = Line.type['Ca'] # normal added mass coeff - CaAx = Line.type['CaAx'] # tangential added mass coeff - Cd = Line.type['Cd'] # normal drag coeff - CdAx = Line.type['CdAx'] # tangential drag coeff - - # extract mean node coordinates - X_mean,Y_mean,Z_mean,T_mean = Line.getLineCoords(0.0,n=N) # coordinates of line nodes and tension values - r_mean = np.vstack((X_mean,Y_mean,Z_mean)).T # coordinates of line nodes - - # evaluate node velocities - r_dynamic = np.ones((len(omegas),N,3),dtype='float')*r_dynamic - v_dynamic = 1j*omegas[:,None,None]*r_dynamic - - # define out of plane normal - h_op = get_horizontal_oop_vec(r_mean[0],r_mean[-1]) # horizontal out-of-plane vector - hh_op = np.outer(h_op,h_op) - - # intialize line matrices - M = np.zeros([3*N, 3*N], dtype='float') # mass matrix - A = np.zeros([3*N, 3*N], dtype='float') # added mass matrix - B = np.zeros([3*N, 3*N], dtype='float') # linearized viscous damping matrix - K = np.zeros([3*N, 3*N], dtype='float') # stiffness matrix - - # Node 0 - dr_e1 = r_mean[1] - r_mean[0] - L_e1 = la.norm(dr_e1) # element 1 length - t_e1 = (dr_e1)/L_e1 # tangential unit vector - p_e1 = np.cross(t_e1,h_op) # in plane normal unit vector - - - ut_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],t_e1) # tangential velocity - uh_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],h_op) # normal horizontal out of plane velocity - up_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],p_e1) # normal in plane velocity - - sigma_ut_e1 = np.sqrt(np.trapz(np.abs(ut_e1)**2*S_zeta,omegas)) - sigma_uh_e1 = np.sqrt(np.trapz(np.abs(uh_e1)**2*S_zeta,omegas)) - sigma_up_e1 = np.sqrt(np.trapz(np.abs(up_e1)**2*S_zeta,omegas)) - - tt_e1 = np.outer(t_e1,t_e1) # local tangential to global components transformation matrix - pp_e1 = np.outer(p_e1,p_e1) # local normal inplane to global components transformation matrix - - M[0:3,0:3] += mden*L_e1/2*np.eye(3) # element 1 mass contribution - - A_e1 = 1025*np.pi/4*deq**2*L_e1/2*(Ca*(hh_op+pp_e1) + CaAx*tt_e1) # element 1 added mass contribution - - B_e1 = 0.5*1025*deq*L_e1/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_e1*hh_op + sigma_up_e1*pp_e1) + - CdAx*sigma_ut_e1*tt_e1) # element 1 damping contribution - - K_e1 = EA/L_e1*tt_e1 + (T_mean[0]/L_e1)*(hh_op+pp_e1) # element 1 stiffness (axial + geometric) - - ## assembling element 1 contributions (rows corresponding to node 0) - A[0:3,0:3] += A_e1 - B[0:3,0:3] += B_e1 - K[0:3,0:3] += K_e1 - K[0:3,3:6] += -K_e1 - - ## add seabed contribution to node 0 - if np.isclose(r_mean[0,2],-depth,seabed_tol): - K[2,2] += kbot - B[2,2] += cbot - - # Internal nodes loop (each internal node has contributions from two elements n-1/2 and n+1/2) - for n in range(1, N-1): - - ## backward element (n-1/2) contributions - dr_bw = r_mean[n-1] - r_mean[n] - L_bw = la.norm(dr_bw) # element 1 length - t_bw = (dr_bw)/L_bw # tangential unit vector - p_bw = np.cross(t_bw,h_op) # in plane normal unit vector - - ut_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],t_bw) # tangential velocity - uh_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],h_op) # normal horizontal out of plane velocity - up_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],p_bw) # normal in plane velocity - - sigma_ut_bw = np.sqrt(np.trapz(np.abs(ut_bw)**2*S_zeta,omegas)) - sigma_uh_bw = np.sqrt(np.trapz(np.abs(uh_bw)**2*S_zeta,omegas)) - sigma_up_bw = np.sqrt(np.trapz(np.abs(up_bw)**2*S_zeta,omegas)) - - tt_bw = np.outer(t_bw,t_bw) # local tangential to global components transformation matrix - pp_bw = np.outer(p_bw,p_bw) # local normal inplane to global components transformation matrix - - M[3*n:3*n+3,3*n:3*n+3] += mden*L_bw/2*np.eye(3) # mass contribution from adjacent elements - - A_bw = 1025*np.pi/4*deq**2*L_bw/2*(Ca*(hh_op+pp_bw) + CaAx*tt_bw) # backward element added mass contribution - - B_bw = 0.5*1025*deq*L_bw/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_bw*hh_op + sigma_up_bw*pp_bw) + - CdAx*sigma_ut_bw*tt_bw) # backward element damping contribution - - K_bw = EA/L_bw*tt_bw + (T_mean[n]/L_bw)*(hh_op+pp_bw) # backward element stiffness (axial + geometric) - - ## forward element (n+1/2) contributions - dr_fw = r_mean[n+1] - r_mean[n] - L_fw = la.norm(dr_fw) # element 1 length - t_fw = (dr_fw)/L_fw # tangential unit vector - p_fw = np.cross(t_fw,h_op) # in plane normal unit vector - - - ut_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],t_fw) # tangential velocity - uh_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],h_op) # normal horizontal out of plane velocity - up_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],p_fw) # normal in plane velocity - - sigma_ut_fw = np.sqrt(np.trapz(np.abs(ut_fw)**2*S_zeta,omegas)) - sigma_uh_fw = np.sqrt(np.trapz(np.abs(uh_fw)**2*S_zeta,omegas)) - sigma_up_fw = np.sqrt(np.trapz(np.abs(up_fw)**2*S_zeta,omegas)) - - tt_fw = np.outer(t_fw,t_fw) # local tangential to global components transformation matrix - pp_fw = np.outer(p_fw,p_fw) # local normal inplane to global components transformation matrix - - M[3*n:3*n+3,3*n:3*n+3] += mden*L_fw/2*np.eye(3) # mass contribution from adjacent elements - - A_fw = 1025*np.pi/4*deq**2*L_fw/2*(Ca*(hh_op+pp_fw) + CaAx*tt_fw) # forward element added mass contribution - - B_fw = 0.5*1025*deq*L_fw/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_fw*hh_op + sigma_up_fw*pp_fw) + - CdAx*sigma_ut_fw*tt_fw) # forward element damping contribution - - K_fw = EA/L_fw*tt_fw + (T_mean[n]/L_fw)*(hh_op+pp_fw) # forward element stiffness (axial + geometric) - - ## assembling bwd and fwd elements contributions (rows corresponding to node n) - A[3*n:3*n+3,3*n:3*n+3] += A_bw + A_fw - B[3*n:3*n+3,3*n:3*n+3] += B_bw + B_fw - K[3*n:3*n+3,3*n:3*n+3] += K_bw + K_fw - K[3*n:3*n+3,3*n-3:3*n] += -K_bw - K[3*n:3*n+3,3*n+3:3*n+6] += -K_fw - - ## add seabed contribution to node n - if np.isclose(r_mean[n,2],-depth,seabed_tol): - K[3*n+2,3*n+2] += kbot - B[3*n+2,3*n+2] += cbot - - # Node N - dr_eN = r_mean[N-1] - r_mean[N-2] - L_eN = la.norm(dr_eN) # element N length - t_eN = (dr_eN)/L_eN # tangential unit vector - p_eN = np.cross(t_eN,h_op) # in plane normal unit vector - - ut_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],t_eN) # tangential velocity - uh_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],h_op) # normal horizontal out of plane velocity - up_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],p_eN) # normal in plane velocity - - sigma_ut_eN = np.sqrt(np.trapz(np.abs(ut_eN)**2*S_zeta,omegas)) - sigma_uh_eN = np.sqrt(np.trapz(np.abs(uh_eN)**2*S_zeta,omegas)) - sigma_up_eN = np.sqrt(np.trapz(np.abs(up_eN)**2*S_zeta,omegas)) - - tt_eN = np.outer(t_eN,t_eN) # local tangential to global components transformation matrix - pp_eN = np.outer(p_eN,p_eN) # local normal inplane to global components transformation matrix - - M[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += mden*L_eN/2*np.eye(3) # element N mass contribution - - A_eN = 1025*np.pi/4*deq**2*L_eN/2*(Ca*(hh_op+pp_eN) + CaAx*tt_eN) # element N added mass contribution - - B_eN = 0.5*1025*deq*L_eN/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_eN*hh_op + sigma_up_eN*pp_eN) + - CdAx*sigma_ut_eN*tt_eN) # element N damping contribution - - K_eN = EA/L_eN*tt_eN + (T_mean[N-1]/L_eN)*(hh_op+pp_eN) # element N stiffness (axial + geometric) - - ## assembling element N contributions (rows corresponding to node N) - A[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += A_eN - B[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += B_eN - K[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += K_eN - K[3*(N-1):3*(N-1)+3,3*(N-1)-3:3*(N-1)] += -K_eN - - ## add seabed contribution to node N - if np.isclose(r_mean[N-1,2],-depth,seabed_tol): - K[3*(N-1)+2,3*(N-1)+2] += kbot - B[3*(N-1)+2,3*(N-1)+2] += cbot - - EA_segs = Line.type['EA']*np.ones(Line.nNodes - 1) - - return M,A,B,K,r_mean,EA_segs - - -def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_tol=1e-4,tol = 0.01,iters=100, w = 0.8, conv_time=True): - """ - Evaluates dynamic tension at all the nodes for an instance of MoorPy's Line or CompositeLine classes. - - Parameters - ---------- - Line : Line/CompositeLine - An instance of MoorPy's Line or CompositeLine classes. - omegas : ndarray - Array of frequencies in rad/s. - S_zeta : ndarray - Wave spectrum array in m^2/(rad/s), must be of the same length as omegas. - RAO_A : ndarray - Translational RAOs for end node A given as a (m,3) array where m is the number of frequencies (must be equal to the length of omegas) . - RAO_B : ndarray - Translational RAOs for end node B given as a (m,3) array where m is the number of frequencies (must be equal to the length of omegas) . - depth : float - Water depth. - kbot : float - Vertical stiffness for points lying on the seabed. - cbot : float - Vertical damping for points lying on the seabed. - seabed_tol : float, optional - Distance from seabed within which a node is considered to be lying on the seabed, by default 1e-4 m. - tol : float, optional - Relative tolerance for iteration, by default 0.01 - iters : int, optional - Maximum number of iterations, by default 100 - w : float, optional - Succesive relaxation coefficient, by default 0.8 - - Returns - ------- - T_nodes_psd: ndarray - Tension spectra at nodes given as (m,n) array where m is the number of frequencies and n is the number of nodes. - T_nodes_std: ndarray - Tension standard deviation at nodes. - s: ndarray: - Node locations along the line. - r_static: ndarray - Nodes mean static position given as (n,3) array where n the number of nodes. - r_dynamic: ndarray - Nodes complex dynamic amplitudes given as (m,n,3) array where m the number of frequencies, n is the number of nodes - r_total: ndarray - Combined static and dynamic positions. - X: ndarray - Solution of the dynamic problem. - """ - N = Line.nNodes - n_dofs = 3*N - - if np.all(RAO_A == 0): - RAO_A = np.zeros_like(RAO_B) - if np.all(RAO_B == 0): - RAO_B = np.zeros_like(RAO_A) - - # intialize iteration matrices - r_dynamic_init = np.ones((len(omegas),N,3)) - M,A,B,K,r_static,EA_segs = Line.getDynamicMatrices(omegas,S_zeta,r_dynamic_init,depth,kbot,cbot,seabed_tol=seabed_tol) # TODO: return EA_segs - X = np.zeros((len(omegas),n_dofs),dtype = 'complex') - r_dynamic = np.zeros(((len(omegas),int(n_dofs/3),3)),dtype = 'complex') - S_Xd = np.zeros((len(omegas),n_dofs),dtype = 'float') - sigma_Xd = np.zeros(n_dofs,dtype = 'float') - sigma_Xd0 = np.zeros(n_dofs,dtype = 'float') - X[:, :3] = RAO_A - X[:,-3:] = RAO_B - - # solving dynamics - start = datetime.now() - for ni in range(iters): - H = - omegas[:,None,None]**2*(M+A)[None,:,:]\ - + 1j*omegas[:,None,None]*B[None,:,:]\ - + K[None,:,:]\ - - F_A = np.einsum('nij,njk->ni',-H[:,3:-3, :3],RAO_A[:,:,None]) - F_B = np.einsum('nij,njk->ni',-H[:,3:-3,-3:],RAO_B[:,:,None]) - F = F_A + F_B - - X[:,3:-3] = solve(H[:,3:-3,3:-3],F) - S_Xd[:] = np.abs(1j*omegas[:,None]*X)**2*S_zeta[:,None] - sigma_Xd[:] = np.sqrt(np.trapz(S_Xd,omegas,axis=0)) - r_dynamic[:] = X.reshape(X.shape[0],int(X.shape[1]/3),3) - if (np.abs(sigma_Xd-sigma_Xd0) <= tol*np.abs(sigma_Xd0)).all(): - break - else: - sigma_Xd0[:] = w * sigma_Xd + (1.-w) * sigma_Xd0 - _,_,B[:],_,_,_ = Line.getDynamicMatrices(omegas,S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=seabed_tol) - if conv_time: - print(f'Finished {ni} dynamic tension iterations in {datetime.now()-start} seconds (w = {w}).') - - # evaluate tension - dw = np.diff(omegas, - prepend= omegas[0] - (omegas[1]-omegas[0]), - append= omegas[-1] + (omegas[-1]-omegas[-2])) - dw = (dw[1:]+dw[:-1])/2 - wave_amps = np.sqrt(S_zeta*dw) #evaluate wave amplitudes of harmonic components from wave spectrum - - r_dynamic *= wave_amps[:,None,None] - r_total = r_static[None,:,:] + r_dynamic - dr_static = r_static[:-1] - r_static[1:] - dr_dynamic = r_dynamic[:,:-1,:] - r_dynamic[:,1:,:] - tangents = dr_static/la.norm(r_static[:-1] - r_static[1:], axis=-1)[:,None] - L_static = la.norm(dr_static, axis=-1) - dL_dynamic = np.einsum('mni,ni->mn', dr_dynamic, tangents) - eps_segs = np.abs(dL_dynamic)/L_static - - T_segs = EA_segs * eps_segs - T_nodes = np.zeros((len(omegas),N)) - T_nodes[:,0] = T_segs[:,0] - T_nodes[:,1:-1] = (T_segs[:,:-1] + T_segs[:,1:])/2 - T_nodes[:,-1] = T_segs[:,-1] - - # S_T = np.zeros((len(omegas),N)) - # S_T[:,1:] = T_e**2/dw[:,None] - # S_T[:,0] = S_T[:,1] - - T_nodes_psd = T_nodes**2/dw[:,None] - T_nodes_std = np.sqrt(np.trapz(T_nodes_psd,omegas,axis=0)) - - - dr = np.diff(r_static,axis=0) - ds = la.norm(dr,axis=1) - s = np.zeros_like(T_nodes_std) - s = np.cumsum(ds) - - return T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X - - -def get_modes(line,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view = False,kbot=3E+06,cbot=3E+05,seabed_tol=1E-04): - M,A,_,K,r_nodes,_ = line.getDynamicMatrices(np.ones(1), np.ones(1),0.,line.sys.depth,kbot,cbot,seabed_tol=seabed_tol) - - if fix_A: - n1 = 1 - else: - n1 = 0 - - if fix_B: - n2 = -1 - else: - n2 = r_nodes.shape[0] - - eigvals,eigvecs = la.eig(K[3*n1:3*n2,3*n1:3*n2],M[3*n1:3*n2,3*n1:3*n2]+A[3*n1:3*n2,3*n1:3*n2]) - stable_eigvals = eigvals[np.real(eigvals)>0] - stable_eigvecs = eigvecs[:,np.real(eigvals)>0] - - idx = stable_eigvals.argsort()[::-1] - stable_eigvals = stable_eigvals[idx] - stable_eigvecs = stable_eigvecs[:,idx] - - freqs = np.sqrt(np.real(stable_eigvals))/2/np.pi - mode_shapes = np.zeros(stable_eigvecs.shape,dtype='float') - - for i in range(stable_eigvecs.shape[1]): - mode_shapes[:,i] = r_nodes[n1:n2].flatten('C') + stable_eigvecs[:,i] - - freqs = np.flip(freqs) - mode_shapes = np.flip(mode_shapes,axis=1) - - if plot_modes: - cols = 4 - rows = plot_modes//cols + bool(plot_modes%cols) - fig,ax = plt.subplots(rows,cols,subplot_kw={"projection": "3d"},figsize=(5*cols,5*rows)) - - i = 0 - for axes in ax: - if not isinstance(axes,Iterable): - axes = [axes] - - for axis in axes: - if i >= plot_modes: - break - - r = r_nodes.copy() - r[n1:n2] = mode_shapes[:,i].reshape([int(len(mode_shapes[:,i])/3),3]) - r = (r-r_nodes)*amp_factor - x = r[:,0] - y = r[:,1] - z = r[:,2] - - x_0 = r_nodes[:,0] - y_0 = r_nodes[:,1] - z_0 = r_nodes[:,2] - - axis.plot(x_0,y_0,z_0,'-ko',label='initial') - axis.plot(x+x_0,y+y_0,z+z_0,'--ro',label='mode shape') - - # R_0 = np.sqrt(x_0**2 + y_0**2) - if adj_view: - # h_min = np.min((x_0,y_0)) - # h_max = np.max((x_0,y_0)) - # axis.set_xlim(h_min,h_max) - # axis.set_ylim(h_min,h_max) - # axis.set_zlim(z_0.min(),z_0.max()) - sigma_x = x.std() - sigma_y = y.std() - sigma_z = z.std() - azim = np.arctan2(sigma_x,sigma_y)*180/np.pi - elev = np.arctan2(np.hypot(sigma_x,sigma_y),sigma_z)*180/np.pi - axis.view_init(elev=elev,azim=azim) - - # axis.set_box_aspect([np.ptp(coord) for coord in [x, y, z]]) - axis.set_xlabel('X (m)') - axis.set_ylabel('Y (m)') - axis.set_zlabel('Z (m)') - axis.set_title(f'f = {freqs[i]:.3f} Hz, T = {1/freqs[i]:.3f} s') - - i+=1 - - fig.tight_layout() - return freqs,mode_shapes,r_nodes,M,A,K,fig,ax - else: - return freqs,mode_shapes,r_nodes,M,A,K diff --git a/moorpy/helpers.py b/moorpy/helpers.py index b6d6ad9..03b789e 100644 --- a/moorpy/helpers.py +++ b/moorpy/helpers.py @@ -1011,7 +1011,8 @@ def lines2subsystem(lines,ms,span=None,case=0): ss.lineTypes[types[-1]] = ms.lineTypes[types[-1]] # use makeGeneric to build the subsystem line - ss.makeGeneric(lengths,types,suspended=case) + nSegs = [ms.lineList[i].nNodes-1 for i in lines] + ss.makeGeneric(lengths,types,suspended=case, nSegs=nSegs) ss.setEndPosition(ms.lineList[lines[0]].rA,endB=0) ss.setEndPosition(ms.lineList[lines[-1]].rB,endB=1) @@ -1467,4 +1468,456 @@ def read_output_file(dirName,fileName, skiplines=-1, hasunits=1, chanlim=999, di unitDict[channels[i]] = units[i] return dataDict, unitDict else: - return data3, ch, channels, units \ No newline at end of file + return data3, ch, channels, units + + +def get_horizontal_oop_vec(p1,p2): + """ + Evaluates the horizontal out of plane vector given the coordinates of two points. + """ + hor_vec = p2 - np.array([p1[0],p1[1],p2[2]]) + ver_vec = p1 - np.array([p1[0],p1[1],p2[2]]) + + if np.isclose(np.linalg.norm(hor_vec),0): # vertical line + n_op = np.array([1,0,0]) + elif np.isclose(np.linalg.norm(ver_vec),0): # horizontal line + oop_vec = np.cross(hor_vec,np.array([0,0,1])) + n_op = oop_vec/np.linalg.norm(oop_vec) + else: + oop_vec = np.cross(hor_vec,ver_vec) + n_op = oop_vec/np.linalg.norm(oop_vec) + return n_op + + +def get_dynamic_matrices(Line, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=1e-4): + """ + Evaluates dynamic matrices for a Line object. + + Parameters + ---------- + Line : Line + An instance of MoorPy's Line class + omegas : ndarray + Array of frequencies in rad/s. + S_zeta : ndarray + Wave spectrum array in m^2/(rad/s), must be of the same length as omegas. + r_dynamic : ndarray + A 3d array of the frequency dependent complex amplitudes of line nodes. The array has a shape of (m,n,3) where m is the number of frequencies, + n is the number of nodes, and the last dimension of 3 correspond to the x,y,z components. + depth : float + Water depth. + kbot : float + Vertical stiffness for points lying on the seabed. + cbot : float + Vertical damping for points lying on the seabed. + seabed_tol : float, optional + Distance from seabed within which a node is considered to be lying on the seabed, by default 1e-4 m. + + Returns + ------- + M: ndarray + Mass matrix. + A: ndarray + Added mass matrix. + B: ndarray + Damping matrix. + K: ndarray + Stiffness matrix. + M: ndarray + Mass matrix. + r_mean: ndarray + Mean static positions of the nodes given as a (m,3) array where m is the number of nodes. + EA_segs: ndarray + Extensional stiffness of segments. + """ + # extract line properties + N = Line.nNodes + mden = Line.type['m'] # line mass density function + deq = Line.type['d_vol'] # line volume equivalent diameter + EA = Line.type['EA'] # extensional stiffness + Ca = Line.type['Ca'] # normal added mass coeff + CaAx = Line.type['CaAx'] # tangential added mass coeff + Cd = Line.type['Cd'] # normal drag coeff + CdAx = Line.type['CdAx'] # tangential drag coeff + + # extract mean node coordinates + X_mean,Y_mean,Z_mean,T_mean = Line.getLineCoords(0.0,n=N) # coordinates of line nodes and tension values + r_mean = np.vstack((X_mean,Y_mean,Z_mean)).T # coordinates of line nodes + + # evaluate node velocities + r_dynamic = np.ones((len(omegas),N,3),dtype='float')*r_dynamic + v_dynamic = 1j*omegas[:,None,None]*r_dynamic + + # define out of plane normal + h_op = get_horizontal_oop_vec(r_mean[0],r_mean[-1]) # horizontal out-of-plane vector + hh_op = np.outer(h_op,h_op) + + # intialize line matrices + M = np.zeros([3*N, 3*N], dtype='float') # mass matrix + A = np.zeros([3*N, 3*N], dtype='float') # added mass matrix + B = np.zeros([3*N, 3*N], dtype='float') # linearized viscous damping matrix + K = np.zeros([3*N, 3*N], dtype='float') # stiffness matrix + + # Node 0 + dr_e1 = r_mean[1] - r_mean[0] + L_e1 = np.linalg.norm(dr_e1) # element 1 length + t_e1 = (dr_e1)/L_e1 # tangential unit vector + p_e1 = np.cross(t_e1,h_op) # in plane normal unit vector + + + ut_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],t_e1) # tangential velocity + uh_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],h_op) # normal horizontal out of plane velocity + up_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],p_e1) # normal in plane velocity + + sigma_ut_e1 = np.sqrt(np.trapz(np.abs(ut_e1)**2*S_zeta,omegas)) + sigma_uh_e1 = np.sqrt(np.trapz(np.abs(uh_e1)**2*S_zeta,omegas)) + sigma_up_e1 = np.sqrt(np.trapz(np.abs(up_e1)**2*S_zeta,omegas)) + + tt_e1 = np.outer(t_e1,t_e1) # local tangential to global components transformation matrix + pp_e1 = np.outer(p_e1,p_e1) # local normal inplane to global components transformation matrix + + M[0:3,0:3] += mden*L_e1/2*np.eye(3) # element 1 mass contribution + + A_e1 = 1025*np.pi/4*deq**2*L_e1/2*(Ca*(hh_op+pp_e1) + CaAx*tt_e1) # element 1 added mass contribution + + B_e1 = 0.5*1025*deq*L_e1/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_e1*hh_op + sigma_up_e1*pp_e1) + + CdAx*sigma_ut_e1*tt_e1) # element 1 damping contribution + + K_e1 = EA/L_e1*tt_e1 + (T_mean[0]/L_e1)*(hh_op+pp_e1) # element 1 stiffness (axial + geometric) + + ## assembling element 1 contributions (rows corresponding to node 0) + A[0:3,0:3] += A_e1 + B[0:3,0:3] += B_e1 + K[0:3,0:3] += K_e1 + K[0:3,3:6] += -K_e1 + + ## add seabed contribution to node 0 + if np.isclose(r_mean[0,2],-depth,seabed_tol): + K[2,2] += kbot + B[2,2] += cbot + + # Internal nodes loop (each internal node has contributions from two elements n-1/2 and n+1/2) + for n in range(1, N-1): + + ## backward element (n-1/2) contributions + dr_bw = r_mean[n-1] - r_mean[n] + L_bw = np.linalg.norm(dr_bw) # element 1 length + t_bw = (dr_bw)/L_bw # tangential unit vector + p_bw = np.cross(t_bw,h_op) # in plane normal unit vector + + ut_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],t_bw) # tangential velocity + uh_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],h_op) # normal horizontal out of plane velocity + up_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],p_bw) # normal in plane velocity + + sigma_ut_bw = np.sqrt(np.trapz(np.abs(ut_bw)**2*S_zeta,omegas)) + sigma_uh_bw = np.sqrt(np.trapz(np.abs(uh_bw)**2*S_zeta,omegas)) + sigma_up_bw = np.sqrt(np.trapz(np.abs(up_bw)**2*S_zeta,omegas)) + + tt_bw = np.outer(t_bw,t_bw) # local tangential to global components transformation matrix + pp_bw = np.outer(p_bw,p_bw) # local normal inplane to global components transformation matrix + + M[3*n:3*n+3,3*n:3*n+3] += mden*L_bw/2*np.eye(3) # mass contribution from adjacent elements + + A_bw = 1025*np.pi/4*deq**2*L_bw/2*(Ca*(hh_op+pp_bw) + CaAx*tt_bw) # backward element added mass contribution + + B_bw = 0.5*1025*deq*L_bw/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_bw*hh_op + sigma_up_bw*pp_bw) + + CdAx*sigma_ut_bw*tt_bw) # backward element damping contribution + + K_bw = EA/L_bw*tt_bw + (T_mean[n]/L_bw)*(hh_op+pp_bw) # backward element stiffness (axial + geometric) + + ## forward element (n+1/2) contributions + dr_fw = r_mean[n+1] - r_mean[n] + L_fw = np.linalg.norm(dr_fw) # element 1 length + t_fw = (dr_fw)/L_fw # tangential unit vector + p_fw = np.cross(t_fw,h_op) # in plane normal unit vector + + + ut_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],t_fw) # tangential velocity + uh_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],h_op) # normal horizontal out of plane velocity + up_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],p_fw) # normal in plane velocity + + sigma_ut_fw = np.sqrt(np.trapz(np.abs(ut_fw)**2*S_zeta,omegas)) + sigma_uh_fw = np.sqrt(np.trapz(np.abs(uh_fw)**2*S_zeta,omegas)) + sigma_up_fw = np.sqrt(np.trapz(np.abs(up_fw)**2*S_zeta,omegas)) + + tt_fw = np.outer(t_fw,t_fw) # local tangential to global components transformation matrix + pp_fw = np.outer(p_fw,p_fw) # local normal inplane to global components transformation matrix + + M[3*n:3*n+3,3*n:3*n+3] += mden*L_fw/2*np.eye(3) # mass contribution from adjacent elements + + A_fw = 1025*np.pi/4*deq**2*L_fw/2*(Ca*(hh_op+pp_fw) + CaAx*tt_fw) # forward element added mass contribution + + B_fw = 0.5*1025*deq*L_fw/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_fw*hh_op + sigma_up_fw*pp_fw) + + CdAx*sigma_ut_fw*tt_fw) # forward element damping contribution + + K_fw = EA/L_fw*tt_fw + (T_mean[n]/L_fw)*(hh_op+pp_fw) # forward element stiffness (axial + geometric) + + ## assembling bwd and fwd elements contributions (rows corresponding to node n) + A[3*n:3*n+3,3*n:3*n+3] += A_bw + A_fw + B[3*n:3*n+3,3*n:3*n+3] += B_bw + B_fw + K[3*n:3*n+3,3*n:3*n+3] += K_bw + K_fw + K[3*n:3*n+3,3*n-3:3*n] += -K_bw + K[3*n:3*n+3,3*n+3:3*n+6] += -K_fw + + ## add seabed contribution to node n + if np.isclose(r_mean[n,2],-depth,seabed_tol): + K[3*n+2,3*n+2] += kbot + B[3*n+2,3*n+2] += cbot + + # Node N + dr_eN = r_mean[N-1] - r_mean[N-2] + L_eN = np.linalg.norm(dr_eN) # element N length + t_eN = (dr_eN)/L_eN # tangential unit vector + p_eN = np.cross(t_eN,h_op) # in plane normal unit vector + + ut_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],t_eN) # tangential velocity + uh_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],h_op) # normal horizontal out of plane velocity + up_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],p_eN) # normal in plane velocity + + sigma_ut_eN = np.sqrt(np.trapz(np.abs(ut_eN)**2*S_zeta,omegas)) + sigma_uh_eN = np.sqrt(np.trapz(np.abs(uh_eN)**2*S_zeta,omegas)) + sigma_up_eN = np.sqrt(np.trapz(np.abs(up_eN)**2*S_zeta,omegas)) + + tt_eN = np.outer(t_eN,t_eN) # local tangential to global components transformation matrix + pp_eN = np.outer(p_eN,p_eN) # local normal inplane to global components transformation matrix + + M[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += mden*L_eN/2*np.eye(3) # element N mass contribution + + A_eN = 1025*np.pi/4*deq**2*L_eN/2*(Ca*(hh_op+pp_eN) + CaAx*tt_eN) # element N added mass contribution + + B_eN = 0.5*1025*deq*L_eN/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_eN*hh_op + sigma_up_eN*pp_eN) + + CdAx*sigma_ut_eN*tt_eN) # element N damping contribution + + K_eN = EA/L_eN*tt_eN + (T_mean[N-1]/L_eN)*(hh_op+pp_eN) # element N stiffness (axial + geometric) + + ## assembling element N contributions (rows corresponding to node N) + A[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += A_eN + B[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += B_eN + K[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += K_eN + K[3*(N-1):3*(N-1)+3,3*(N-1)-3:3*(N-1)] += -K_eN + + ## add seabed contribution to node N + if np.isclose(r_mean[N-1,2],-depth,seabed_tol): + K[3*(N-1)+2,3*(N-1)+2] += kbot + B[3*(N-1)+2,3*(N-1)+2] += cbot + + EA_segs = Line.type['EA']*np.ones(Line.nNodes - 1) + + return M,A,B,K,r_mean,EA_segs + + +def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_tol=1e-4,tol = 0.01,iters=100, w = 0.8, conv_time=True): + """ + Evaluates dynamic tension at all the nodes for an instance of MoorPy's Line or CompositeLine classes. + + Parameters + ---------- + Line : Line/CompositeLine + An instance of MoorPy's Line or CompositeLine classes. + omegas : ndarray + Array of frequencies in rad/s. + S_zeta : ndarray + Wave spectrum array in m^2/(rad/s), must be of the same length as omegas. + RAO_A : ndarray + Translational RAOs for end node A given as a (m,3) array where m is the number of frequencies (must be equal to the length of omegas) . + RAO_B : ndarray + Translational RAOs for end node B given as a (m,3) array where m is the number of frequencies (must be equal to the length of omegas) . + depth : float + Water depth. + kbot : float + Vertical stiffness for points lying on the seabed. + cbot : float + Vertical damping for points lying on the seabed. + seabed_tol : float, optional + Distance from seabed within which a node is considered to be lying on the seabed, by default 1e-4 m. + tol : float, optional + Relative tolerance for iteration, by default 0.01 + iters : int, optional + Maximum number of iterations, by default 100 + w : float, optional + Succesive relaxation coefficient, by default 0.8 + + Returns + ------- + T_nodes_psd: ndarray + Tension spectra at nodes given as (m,n) array where m is the number of frequencies and n is the number of nodes. + T_nodes_std: ndarray + Tension standard deviation at nodes. + s: ndarray: + Node locations along the line. + r_static: ndarray + Nodes mean static position given as (n,3) array where n the number of nodes. + r_dynamic: ndarray + Nodes complex dynamic amplitudes given as (m,n,3) array where m the number of frequencies, n is the number of nodes + r_total: ndarray + Combined static and dynamic positions. + X: ndarray + Solution of the dynamic problem. + """ + N = Line.nNodes + n_dofs = 3*N + + if np.all(RAO_A == 0): + RAO_A = np.zeros_like(RAO_B) + if np.all(RAO_B == 0): + RAO_B = np.zeros_like(RAO_A) + + # intialize iteration matrices + r_dynamic_init = np.ones((len(omegas),N,3)) + M,A,B,K,r_static,EA_segs = Line.getDynamicMatrices(omegas,S_zeta,r_dynamic_init,depth,kbot,cbot,seabed_tol=seabed_tol) # TODO: return EA_segs + X = np.zeros((len(omegas),n_dofs),dtype = 'complex') + r_dynamic = np.zeros(((len(omegas),int(n_dofs/3),3)),dtype = 'complex') + S_Xd = np.zeros((len(omegas),n_dofs),dtype = 'float') + sigma_Xd = np.zeros(n_dofs,dtype = 'float') + sigma_Xd0 = np.zeros(n_dofs,dtype = 'float') + X[:, :3] = RAO_A + X[:,-3:] = RAO_B + + # solving dynamics + start = time.time() + for ni in range(iters): + H = - omegas[:,None,None]**2*(M+A)[None,:,:]\ + + 1j*omegas[:,None,None]*B[None,:,:]\ + + K[None,:,:]\ + + F_A = np.einsum('nij,njk->ni',-H[:,3:-3, :3],RAO_A[:,:,None]) + F_B = np.einsum('nij,njk->ni',-H[:,3:-3,-3:],RAO_B[:,:,None]) + F = F_A + F_B + + X[:,3:-3] = np.linalg.solve(H[:,3:-3,3:-3],F) + S_Xd[:] = np.abs(1j*omegas[:,None]*X)**2*S_zeta[:,None] + sigma_Xd[:] = np.sqrt(np.trapz(S_Xd,omegas,axis=0)) + r_dynamic[:] = X.reshape(X.shape[0],int(X.shape[1]/3),3) + if (np.abs(sigma_Xd-sigma_Xd0) <= tol*np.abs(sigma_Xd0)).all(): + break + else: + sigma_Xd0[:] = w * sigma_Xd + (1.-w) * sigma_Xd0 + _,_,B[:],_,_,_ = Line.getDynamicMatrices(omegas,S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=seabed_tol) + if conv_time: + print(f'Finished {ni} dynamic tension iterations in {time.time()-start} seconds (w = {w}).') + + # evaluate tension + dw = np.diff(omegas, + prepend= omegas[0] - (omegas[1]-omegas[0]), + append= omegas[-1] + (omegas[-1]-omegas[-2])) + dw = (dw[1:]+dw[:-1])/2 + wave_amps = np.sqrt(S_zeta*dw) #evaluate wave amplitudes of harmonic components from wave spectrum + + r_dynamic *= wave_amps[:,None,None] + r_total = r_static[None,:,:] + r_dynamic + dr_static = r_static[:-1] - r_static[1:] + dr_dynamic = r_dynamic[:,:-1,:] - r_dynamic[:,1:,:] + tangents = dr_static/np.linalg.norm(r_static[:-1] - r_static[1:], axis=-1)[:,None] + L_static = np.linalg.norm(dr_static, axis=-1) + dL_dynamic = np.einsum('mni,ni->mn', dr_dynamic, tangents) + eps_segs = np.abs(dL_dynamic)/L_static + + T_segs = EA_segs * eps_segs + T_nodes = np.zeros((len(omegas),N)) + T_nodes[:,0] = T_segs[:,0] + T_nodes[:,1:-1] = (T_segs[:,:-1] + T_segs[:,1:])/2 + T_nodes[:,-1] = T_segs[:,-1] + + # S_T = np.zeros((len(omegas),N)) + # S_T[:,1:] = T_e**2/dw[:,None] + # S_T[:,0] = S_T[:,1] + + T_nodes_psd = T_nodes**2/dw[:,None] + T_nodes_std = np.sqrt(np.trapz(T_nodes_psd,omegas,axis=0)) + + + dr = np.diff(r_static,axis=0) + ds = np.linalg.norm(dr,axis=1) + s = np.zeros_like(T_nodes_std) + s = np.cumsum(ds) + + return T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X + + +def get_modes(line,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view = False,kbot=3E+06,cbot=3E+05,seabed_tol=1E-04): + import matplotlib.pyplot as plt + from collections.abc import Iterable + + M,A,_,K,r_nodes,_ = line.getDynamicMatrices(np.ones(1), np.ones(1),0.,line.sys.depth,kbot,cbot,seabed_tol=seabed_tol) + + if fix_A: + n1 = 1 + else: + n1 = 0 + + if fix_B: + n2 = -1 + else: + n2 = r_nodes.shape[0] + + eigvals,eigvecs = la.eig(K[3*n1:3*n2,3*n1:3*n2],M[3*n1:3*n2,3*n1:3*n2]+A[3*n1:3*n2,3*n1:3*n2]) + stable_eigvals = eigvals[np.real(eigvals)>0] + stable_eigvecs = eigvecs[:,np.real(eigvals)>0] + + idx = stable_eigvals.argsort()[::-1] + stable_eigvals = stable_eigvals[idx] + stable_eigvecs = stable_eigvecs[:,idx] + + freqs = np.sqrt(np.real(stable_eigvals))/2/np.pi + mode_shapes = np.zeros(stable_eigvecs.shape,dtype='float') + + for i in range(stable_eigvecs.shape[1]): + mode_shapes[:,i] = r_nodes[n1:n2].flatten('C') + stable_eigvecs[:,i] + + freqs = np.flip(freqs) + mode_shapes = np.flip(mode_shapes,axis=1) + + if plot_modes: + cols = 4 + rows = plot_modes//cols + bool(plot_modes%cols) + fig,ax = plt.subplots(rows,cols,subplot_kw={"projection": "3d"},figsize=(5*cols,5*rows)) + + i = 0 + for axes in ax: + if not isinstance(axes,Iterable): + axes = [axes] + + for axis in axes: + if i >= plot_modes: + break + + r = r_nodes.copy() + r[n1:n2] = mode_shapes[:,i].reshape([int(len(mode_shapes[:,i])/3),3]) + r = (r-r_nodes)*amp_factor + x = r[:,0] + y = r[:,1] + z = r[:,2] + + x_0 = r_nodes[:,0] + y_0 = r_nodes[:,1] + z_0 = r_nodes[:,2] + + axis.plot(x_0,y_0,z_0,'-ko',label='initial') + axis.plot(x+x_0,y+y_0,z+z_0,'--ro',label='mode shape') + + # R_0 = np.sqrt(x_0**2 + y_0**2) + if adj_view: + # h_min = np.min((x_0,y_0)) + # h_max = np.max((x_0,y_0)) + # axis.set_xlim(h_min,h_max) + # axis.set_ylim(h_min,h_max) + # axis.set_zlim(z_0.min(),z_0.max()) + sigma_x = x.std() + sigma_y = y.std() + sigma_z = z.std() + azim = np.arctan2(sigma_x,sigma_y)*180/np.pi + elev = np.arctan2(np.hypot(sigma_x,sigma_y),sigma_z)*180/np.pi + axis.view_init(elev=elev,azim=azim) + + # axis.set_box_aspect([np.ptp(coord) for coord in [x, y, z]]) + axis.set_xlabel('X (m)') + axis.set_ylabel('Y (m)') + axis.set_zlabel('Z (m)') + axis.set_title(f'f = {freqs[i]:.3f} Hz, T = {1/freqs[i]:.3f} s') + + i+=1 + + fig.tight_layout() + return freqs,mode_shapes,r_nodes,M,A,K,fig,ax + else: + return freqs,mode_shapes,r_nodes,M,A,K diff --git a/moorpy/line.py b/moorpy/line.py index 721b5bf..ffd917c 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -6,7 +6,8 @@ from moorpy.nonlinear import nonlinear from moorpy.helpers import (unitVector, LineError, CatenaryError, rotationMatrix, makeTower, read_mooring_file, - quiver_data_to_segments, printVec, printMat) + quiver_data_to_segments, printVec, printMat, + get_dynamic_matrices, get_dynamic_tension, get_modes) from os import path @@ -1013,23 +1014,17 @@ def revertToStaticStiffness(self): def getDynamicMatrices(self, *args, **kwargs): '''Compute M,A,B,K matrices for Line. See get_dynamic_matrices().''' - from moorpy.dynamic_tension_functions import get_dynamic_matrices - return get_dynamic_matrices(self, *args, **kwargs) def dynamicSolve(self, *args, **kwargs): '''Compute complex amplitudes of line nodes. See get_dynamic_tension().''' - from moorpy.dynamic_tension_functions import get_dynamic_tension - return get_dynamic_tension(self, *args, **kwargs) def getModes(self,*args, **kwargs): '''Compute (and optionally plot) the line's mode shapes. See get_modes().''' - from moorpy.dynamic_tension_functions import get_modes - return get_modes(self, *args, **kwargs) diff --git a/moorpy/subsystem.py b/moorpy/subsystem.py index f5874b1..0ffe509 100644 --- a/moorpy/subsystem.py +++ b/moorpy/subsystem.py @@ -119,7 +119,7 @@ def __init__(self, mooringSys=None, num=0, depth=0, rho=1025, g=9.81, self.qs = 1 # flag that it's a MoorPy analysis, so System methods don't complain. One day should replace this <<< - def makeGeneric(self, lengths, types, suspended=0): + def makeGeneric(self, lengths, types, suspended=0, nSegs=40): '''Creates a cable of n components going between an anchor point and a floating body (or a bridle point). If shared, it goes between two floating bodies. @@ -184,8 +184,16 @@ def makeGeneric(self, lengths, types, suspended=0): else: raise Exception(f"Can't find lineType '{types[i]}' in the SubSystem.") + # add the line segment using the reference to its lineType dict # add the line segment using the reference to its lineType dict - self.addLine(lengths[i], self.lineTypes[i]) + if nSegs is None: + self.addLine(lengths[i], self.lineTypes[i]) + elif isinstance(nSegs, (int, float)): + self.addLine(lengths[i], self.lineTypes[i], nSegs=nSegs) + elif isinstance(nSegs, list): + self.addLine(lengths[i], self.lineTypes[i], nSegs=nSegs[i]) + else: + raise ValueError("Invalid type for nSegs. Expected None, a number, or a list.") # add the upper end point of the segment if i==self.nLines-1: # if this is the upper-most line @@ -239,7 +247,7 @@ def makeFromSystem(self, sys, point_id): line_id = point.attached[-1] # get next line's id line = sys.lineList[line_id - 1] # get line object self.lineTypes[line.type['name']] = dict(line.type) # copy the lineType dict - self.addLine(line.L0, self.lineTypes[line.type['name']]) # add the line + self.addLine(line.L0, self.lineTypes[line.type['name']], nSegs=line.nNodes-1) # add the line # get the next point @@ -249,7 +257,7 @@ def makeFromSystem(self, sys, point_id): pointB_id = attached_points[0] # get second point id point = sys.pointList[pointB_id-1] # get second point object - if line(point.attached) == 1: # must be the endpoint + if len(point.attached) == 1: # must be the endpoint self.addPoint(-1, point.r) else: # intermediate point along line self.addPoint(0, point.r, DOFs=[0,2]) # may need to ensure there's no y component diff --git a/moorpy/system.py b/moorpy/system.py index 31763eb..7563101 100644 --- a/moorpy/system.py +++ b/moorpy/system.py @@ -18,7 +18,6 @@ from moorpy.point import Point from moorpy.line import Line from moorpy.lineType import LineType -from moorpy.compositeLine import CompositeLine import matplotlib as mpl #import moorpy.MoorSolve as msolve from moorpy.helpers import (rotationMatrix, rotatePosition, getH, printVec, From c2b2454fc1d303370a8f7ae614d08f3bba65b4ae Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Thu, 23 May 2024 14:42:52 -0600 Subject: [PATCH 04/14] Functions to get system dynamic matrices Applied the same idea of getSystemStiffnessA to lump the dynamic matrices (inertia, added mass, damping, and stiffness) in 3x3 matrices at the fairlead --- moorpy/point.py | 47 ++++++++++ moorpy/system.py | 225 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 272 insertions(+) diff --git a/moorpy/point.py b/moorpy/point.py index 83a8750..637ca10 100644 --- a/moorpy/point.py +++ b/moorpy/point.py @@ -340,6 +340,53 @@ def getStiffnessA(self, lines_only=False, xyz=False): return K else: # otherwise only return rows/columns of active DOFs return K[:,self.DOFs][self.DOFs,:] + + def getDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot, xyz=False): + '''Gets inertia, added mass, damping, and stiffness matrices of Point due only to mooring lines (with other objects fixed) + using a lumped mass model. + + Returns + ------- + K : matrix + 3x3 analytic stiffness matrix. + + ''' + + #print("Getting Point "+str(self.number)+" analytic stiffness matrix...") + + def lump_matrix(matrix, endB): + if endB == 1: + matrix_woAnc = matrix[3:, 3:] + matrix_inv = np.linalg.inv(matrix_woAnc) + matrix_inv_coupled = matrix_inv[-3:, -3:] + else: + matrix_woAnc = matrix[:-3, :-3] + matrix_inv = np.linalg.inv(matrix_woAnc) + matrix_inv_coupled = matrix_inv[0:3, 0:3] + return np.linalg.inv(matrix_inv_coupled) + + M = np.zeros([3,3]) + A = np.zeros([3,3]) + B = np.zeros([3,3]) + K = np.zeros([3,3]) + + # append the stiffness matrix of each line attached to the point + for lineID,endB in zip(self.attached,self.attachedEndB): + line = self.sys.lineList[lineID-1] + if r_dynamic_init == None: + r_init = np.ones((len(omegas),line.nNodes,3)) + M_all, A_all, B_all, K_all, _, _ = line.getDynamicMatrices(omegas,S_zeta,r_init,depth,kbot,cbot) + + M += lump_matrix(M_all, endB) + A += lump_matrix(A_all, endB) + B += lump_matrix(B_all, endB) + K += lump_matrix(K_all, endB) + + if sum(np.isnan(K).ravel()) > 0: breakpoint() + if xyz: # if asked to output all DOFs, do it + return M, A, B, K + else: # otherwise only return rows/columns of active DOFs + return M[:,self.DOFs][self.DOFs,:], A[:,self.DOFs][self.DOFs,:], B[:,self.DOFs][self.DOFs,:], K[:,self.DOFs][self.DOFs,:] def getCost(self): diff --git a/moorpy/system.py b/moorpy/system.py index 7563101..90570f5 100644 --- a/moorpy/system.py +++ b/moorpy/system.py @@ -2877,6 +2877,231 @@ def getSystemStiffnessA(self, DOFtype="free", lines_only=False, rho=1025, g=9.81 return K + def getCoupledDynamicMatrices(self, omegas, S_zeta, depth, kbot, cbot, r_dynamic_init=None, lines_only=False): + '''Write something here later + ''' + self.nDOF, self.nCpldDOF, DOFtypes = self.getDOFs() + + n = self.nDOF + self.nCpldDOF + + if self.display > 2: + print("Getting mooring system stiffness matrix...") + + # get full system stiffness matrix + M_all, A_all, B_all, K_all = self.getSystemDynamicMatrices(omegas, S_zeta, r_dynamic_init, depth, kbot, cbot, DOFtype="both", lines_only=lines_only) + + # invert matrix + M_inv_all = np.linalg.inv(M_all) + A_inv_all = np.linalg.inv(A_all) + B_inv_all = np.linalg.inv(B_all) + K_inv_all = np.linalg.inv(K_all) + + # remove free DOFs (this corresponds to saying that the sum of forces on these DOFs will remain zero) + #indices = list(range(n)) # list of DOF indices that will remain active for this step + mask = [True]*n # this is a mask to be applied to the array K indices + + for i in range(n-1, -1, -1): # go through DOFs and flag free ones for exclusion + if DOFtypes[i] == 0: + mask[i] = False + #del indices[i] + + M_inv_coupled = M_inv_all[mask,:][:,mask] + A_inv_coupled = A_inv_all[mask,:][:,mask] + B_inv_coupled = B_inv_all[mask,:][:,mask] + K_inv_coupled = K_inv_all[mask,:][:,mask] + + # invert reduced matrix to get coupled stiffness matrix (with free DOFs assumed to equilibrate linearly) + M_coupled = np.linalg.inv(M_inv_coupled) + A_coupled = np.linalg.inv(A_inv_coupled) + B_coupled = np.linalg.inv(B_inv_coupled) + K_coupled = np.linalg.inv(K_inv_coupled) + + #if tensions: + # return K_coupled, J + #else: + return M_coupled, A_coupled, B_coupled, K_coupled + + + def getSystemDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot, DOFtype="free", lines_only=False, rho=1025, g=9.81): + '''Write something here later + ''' + + # note: This is missing some pieces, and needs to check more. + # So far this seems to not capture yaw stiffness for non-bridle configs... + # it would require proper use of chain rule for the derivatives + + + # find the total number of free and coupled DOFs in case any object types changed + self.nDOF, self.nCpldDOF, _ = self.getDOFs() + + #self.solveEquilibrium() # should we make sure the system is in equilibrium? + + # allocate stiffness matrix according to the DOFtype specified + if DOFtype=="free": + M = np.zeros([self.nDOF, self.nDOF]) + A = np.zeros([self.nDOF, self.nDOF]) + B = np.zeros([self.nDOF, self.nDOF]) + K = np.zeros([self.nDOF, self.nDOF]) + d = [0] + elif DOFtype=="coupled": + M = np.zeros([self.nCpldDOF, self.nCpldDOF]) + A = np.zeros([self.nCpldDOF, self.nCpldDOF]) + B = np.zeros([self.nCpldDOF, self.nCpldDOF]) + K = np.zeros([self.nCpldDOF, self.nCpldDOF]) + d = [-1] + elif DOFtype=="both": + M = np.zeros([self.nDOF+self.nCpldDOF, self.nDOF+self.nCpldDOF]) + A = np.zeros([self.nDOF+self.nCpldDOF, self.nDOF+self.nCpldDOF]) + B = np.zeros([self.nDOF+self.nCpldDOF, self.nDOF+self.nCpldDOF]) + K = np.zeros([self.nDOF+self.nCpldDOF, self.nDOF+self.nCpldDOF]) + d = [0,-1] + else: + raise ValueError("getSystemDynamicMatrices called with invalid DOFtype input. Must be free, coupled, or both") + + + # The following will go through and get the lower-triangular stiffness terms, + # calculated as the force/moment on Body/Point 2 from translation/rotation of Body/Point 1. + + # go through DOFs, looking for lines that couple to anchors or other DOFs + + i = 0 # start counting number of DOFs at zero + + # go through each movable body in the system + for body1 in self.bodyList: + if body1.type in d: # >>>> when DOFtype==both, this approach gives different indexing than what is in setPositions/getForces and getSystemStiffness <<<<< + + #i = (body1.number-1)*6 # start counting index for body DOFs based on body number to keep indexing consistent + + # get body's self-stiffness matrix (now only cross-coupling terms will be handled on a line-by-line basis) + K6 = body1.getStiffnessA(lines_only=lines_only) + K[i:i+body1.nDOF, i:i+body1.nDOF] += K6 + + + # go through each attached point + for pointID1,rPointRel1 in zip(body1.attachedP,body1.rPointRel): + point1 = self.pointList[pointID1-1] + + r1 = rotatePosition(rPointRel1, body1.r6[3:]) # relative position of Point about body ref point in unrotated reference frame + H1 = getH(r1) # produce alternator matrix of current point's relative position to current body + + + for lineID in point1.attached: # go through each attached line to the Point, looking for when its other end is attached to something that moves + + endFound = 0 # simple flag to indicate when the other end's attachment has been found + j = i + body1.nDOF # first index of the DOFs this line is attached to. Start it off at the next spot after body1's DOFs + + # get cross-coupling stiffness of line: force on end attached to body1 due to motion of other end + if point1.attachedEndB == 1: + KB = self.lineList[lineID-1].KBA + else: + KB = self.lineList[lineID-1].KBA.T + + # look through Bodies further on in the list (coupling with earlier Bodies will already have been taken care of) + for body2 in self.bodyList[self.bodyList.index(body1)+1: ]: + if body2.type in d: + + # go through each attached Point + for pointID2,rPointRel2 in zip(body2.attachedP,body2.rPointRel): + point2 = self.pointList[pointID2-1] + + if lineID in point2.attached: # if the line is also attached to this Point2 in Body2 + + # following are analagous to what's in functions getH and translateMatrix3to6 except for cross coupling between two bodies + r2 = rotatePosition(rPointRel2, body2.r6[3:]) # relative position of Point about body ref point in unrotated reference frame + H2 = getH(r2) + + # loads on body1 due to motions of body2 + K66 = np.block([[ KB , np.matmul(KB, H2)], + [np.matmul(H1.T, KB), np.matmul(np.matmul(H2, KB), H1.T)]]) + + # Trim for only enabled DOFs of the two bodies + K66 = K66[body1.DOFs,:][:,body2.DOFs] + + K[i:i+body1.nDOF, j:j+body2.nDOF] += K66 # f on B1 due to x of B2 + K[j:j+body2.nDOF, i:i+body1.nDOF] += K66.T # mirror + + # note: the additional rotational stiffness due to change in moment arm does not apply to this cross-coupling case + endFound = 1 # signal that the line has been handled so we can move on to the next thing + break + + j += body2.nDOF # if this body has DOFs we're considering, then count them + + + # look through free Points + if endFound==0: # if the end of this line hasn't already been found attached to a body + for point2 in self.pointList: + if point2.type in d: # if it's a free point and + if lineID in point2.attached: # the line is also attached to it + + # only add up one off-diagonal sub-matrix for now, then we'll mirror at the end + K63 = np.vstack([KB, np.matmul(H1.T, KB)]) + + # Trim for only enabled DOFs of the point and body + K63 = K63[body1.DOFs,:][:,point2.DOFs] + + K[i:i+ body1.nDOF, j:j+point2.nDOF] += K63 # f on B1 due to x of P2 + K[j:j+point2.nDOF, i:i+ body1.nDOF] += K63.T # mirror + + break + + j += point2.nDOF # if this point has DOFs we're considering, then count them + + # note: No cross-coupling with fixed points. The body's own stiffness matrix is now calculated at the start. + + i += body1.nDOF # moving along to the next body... + + + # go through each movable point in the system + for point in self.pointList: + if point.type in d: + + n = point.nDOF + + # >>> TODO: handle case of free end point resting on seabed <<< + + # get point's self-stiffness matrix + M1, A1, B1, K1 = point.getDynamicMatrices(omegas, S_zeta, r_dynamic_init, depth, kbot, cbot) + M[i:i+n,i:i+n] += M1 + A[i:i+n,i:i+n] += A1 + B[i:i+n,i:i+n] += B1 + K[i:i+n,i:i+n] += K1 + + # # go through attached lines and add cross-coupling terms + for lineID in point.attached: + + j = i + n + + # go through movable points to see if one is attached + for point2 in self.pointList[self.pointList.index(point)+1: ]: + if point2.type in d: + if lineID in point2.attached: # if this point is at the other end of the line + + # get cross-coupling stiffness of line: force on end attached to point1 due to motion of other end + MB, AB, BB, KB = point2.getDynamicMatrices(omegas, S_zeta, r_dynamic_init, depth, kbot, cbot) + if point.attachedEndB != 1: + MB = MB.T + AB = AB.T + BB = BB.T + KB = KB.T + + # Trim stiffness matrix to only use the enabled DOFs of each point + MB = MB[point.DOFs,:][:,point2.DOFs] + AB = AB[point.DOFs,:][:,point2.DOFs] + BB = BB[point.DOFs,:][:,point2.DOFs] + KB = KB[point.DOFs,:][:,point2.DOFs] + + M[i:i+n , j:j+point2.nDOF] += MB # force on P1 due to movement of P2 + M[j:j+point2.nDOF, i:i+n ] += MB.T # mirror (f on P2 due to x of P1) + A[i:i+n , j:j+point2.nDOF] += AB + A[j:j+point2.nDOF, i:i+n ] += AB.T + B[i:i+n , j:j+point2.nDOF] += BB + B[j:j+point2.nDOF, i:i+n ] += BB.T + K[i:i+n , j:j+point2.nDOF] += KB + K[j:j+point2.nDOF, i:i+n ] += KB.T + j += point2.nDOF # if this point has DOFs we're considering, then count them + i += n + + return M, A, B, K def getAnchorLoads(self, sfx, sfy, sfz, N): ''' Calculates anchor loads From dbacf6d07fb5609269b73737635a50443e5da3f7 Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Fri, 24 May 2024 17:42:23 -0600 Subject: [PATCH 05/14] Fix double counting of an effect + deal with matrices w/ null rows/cols - In system.getSystemDynamicMatrices, the effect of other attached points (the loop on point2) is already lumped in point.getDynamicMatrices - Replaced np.linalg.inv with np.linalg.pinv (which computes the pseudo-inverse of the matrix) because we can have some cases where rows and columns of the matrix are zero (e.g., when adopting Ca=0 in a given direction) - The procedure to lump the dynamic matrices in a point involves inverting the matrices twice. If the input matrix is not exactly symmetrical, this can greatly amplify the spurious asymmetry. For this reason, we make the matrices computed with the lumped mass symmetrical. - For some reason that I don't know yet, the procedure does not work if we use the whole matrix. We need to remove the rows and columns associated with the dofs of the node opposite to the node where we want to lump the matrices. --- moorpy/point.py | 22 +++++++++++----------- moorpy/system.py | 35 +---------------------------------- 2 files changed, 12 insertions(+), 45 deletions(-) diff --git a/moorpy/point.py b/moorpy/point.py index 637ca10..358439a 100644 --- a/moorpy/point.py +++ b/moorpy/point.py @@ -341,7 +341,7 @@ def getStiffnessA(self, lines_only=False, xyz=False): else: # otherwise only return rows/columns of active DOFs return K[:,self.DOFs][self.DOFs,:] - def getDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot, xyz=False): + def getDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot): '''Gets inertia, added mass, damping, and stiffness matrices of Point due only to mooring lines (with other objects fixed) using a lumped mass model. @@ -355,15 +355,19 @@ def getDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot, #print("Getting Point "+str(self.number)+" analytic stiffness matrix...") def lump_matrix(matrix, endB): - if endB == 1: - matrix_woAnc = matrix[3:, 3:] - matrix_inv = np.linalg.inv(matrix_woAnc) + # The matrices should be symmetrical, but they can be slightly off due to numerical errors. + # Because we are going to invert them twice, we force them to be symmetrical to avoid amplifying the errors. + matrix = (matrix + matrix.T)/2 + if endB == 1: + matrix_woAnc = matrix[3:, 3:] # For some reason that I don't know yet, we need to remove the elements in the opposite side of the matrix. Otherwise it doesn't work + matrix_inv = np.linalg.pinv(matrix_woAnc) matrix_inv_coupled = matrix_inv[-3:, -3:] else: matrix_woAnc = matrix[:-3, :-3] - matrix_inv = np.linalg.inv(matrix_woAnc) + matrix_inv = np.linalg.pinv(matrix_woAnc) matrix_inv_coupled = matrix_inv[0:3, 0:3] - return np.linalg.inv(matrix_inv_coupled) + matrix_out = np.linalg.inv(matrix_inv_coupled) + return matrix_out M = np.zeros([3,3]) A = np.zeros([3,3]) @@ -382,11 +386,7 @@ def lump_matrix(matrix, endB): B += lump_matrix(B_all, endB) K += lump_matrix(K_all, endB) - if sum(np.isnan(K).ravel()) > 0: breakpoint() - if xyz: # if asked to output all DOFs, do it - return M, A, B, K - else: # otherwise only return rows/columns of active DOFs - return M[:,self.DOFs][self.DOFs,:], A[:,self.DOFs][self.DOFs,:], B[:,self.DOFs][self.DOFs,:], K[:,self.DOFs][self.DOFs,:] + return M, A, B, K def getCost(self): diff --git a/moorpy/system.py b/moorpy/system.py index 90570f5..e831c02 100644 --- a/moorpy/system.py +++ b/moorpy/system.py @@ -3065,40 +3065,7 @@ def getSystemDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, A[i:i+n,i:i+n] += A1 B[i:i+n,i:i+n] += B1 K[i:i+n,i:i+n] += K1 - - # # go through attached lines and add cross-coupling terms - for lineID in point.attached: - - j = i + n - - # go through movable points to see if one is attached - for point2 in self.pointList[self.pointList.index(point)+1: ]: - if point2.type in d: - if lineID in point2.attached: # if this point is at the other end of the line - - # get cross-coupling stiffness of line: force on end attached to point1 due to motion of other end - MB, AB, BB, KB = point2.getDynamicMatrices(omegas, S_zeta, r_dynamic_init, depth, kbot, cbot) - if point.attachedEndB != 1: - MB = MB.T - AB = AB.T - BB = BB.T - KB = KB.T - - # Trim stiffness matrix to only use the enabled DOFs of each point - MB = MB[point.DOFs,:][:,point2.DOFs] - AB = AB[point.DOFs,:][:,point2.DOFs] - BB = BB[point.DOFs,:][:,point2.DOFs] - KB = KB[point.DOFs,:][:,point2.DOFs] - - M[i:i+n , j:j+point2.nDOF] += MB # force on P1 due to movement of P2 - M[j:j+point2.nDOF, i:i+n ] += MB.T # mirror (f on P2 due to x of P1) - A[i:i+n , j:j+point2.nDOF] += AB - A[j:j+point2.nDOF, i:i+n ] += AB.T - B[i:i+n , j:j+point2.nDOF] += BB - B[j:j+point2.nDOF, i:i+n ] += BB.T - K[i:i+n , j:j+point2.nDOF] += KB - K[j:j+point2.nDOF, i:i+n ] += KB.T - j += point2.nDOF # if this point has DOFs we're considering, then count them + i += n return M, A, B, K From 4bddb152cbaf607ca24f0303f93ed10b6ce197de Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Sat, 1 Jun 2024 18:05:30 -0600 Subject: [PATCH 06/14] Functions to lump matrices at line ends and body + code reorganization - Implemented getDynamicMatricesLumped method in the line class that lumps the added mass, inertia, damping, and stiffness matrices of a line at its extremities based on the matrices for each node (applied for a line modeled using a lumped mass model) - Method getDynamicMatrices in the body class similar to getStiffnessA, but for the dynamic matrices - Reorganized the code in point.getDynamicMatrices to use the new line.getDynamicMatricesLumped method - Changed some np.liang.inv for np.linalg.pinv. We need that for cases where the matrices are singular. For instance, the added mass matrix of the line is singular if the axial added mass coefficient is zero. WIP: There seems to be something wrong in line.getDynamicMatricesLumped if part of the line rests on the seabed, as the resulting stiffness matrix is not consistent with getStiffnessA. The agreement is much better when we totally remove those nodes, but they are still off. Need to figure out what is the correct way of dealing with that. --- moorpy/body.py | 76 ++++++++++++++++++++++++++++++++++- moorpy/line.py | 57 ++++++++++++++++++++++++++ moorpy/point.py | 37 ++++------------- moorpy/system.py | 102 +++++++++++++++++++++++++++++++++++------------ 4 files changed, 217 insertions(+), 55 deletions(-) diff --git a/moorpy/body.py b/moorpy/body.py index 84bcccc..8c76314 100644 --- a/moorpy/body.py +++ b/moorpy/body.py @@ -370,8 +370,82 @@ def getStiffnessA(self, lines_only=False, all_DOFs=False): return K else: # only return rows/columns of active DOFs return K[:,self.DOFs][self.DOFs,:] - + + + def getDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dynamic_init=None, lines_only=False, all_DOFs=False): + '''Gets the dynamic matrices of the Body with other objects fixed + using a lumped mass approach. + + Parameters + ---------- + lines_only : boolean, optional + If true, the Body's contribution to its stiffness is ignored. + all_DOFs : boolean, optional + True: return all forces/moments; False: only those in DOFs list. + + Returns + ------- + K : matrix + 6x6 analytic stiffness matrix. + ''' + + M, A, B, K = (np.zeros([6,6]) for _ in range(4)) + + # Contributions from attached points (and any of their attached lines) + for PointID,rPointRel in zip(self.attachedP,self.rPointRel): + r = rotatePosition(rPointRel, self.r6[3:]) # relative position of Point about body ref point in unrotated reference frame + f3 = self.sys.pointList[PointID-1].getForces() # total force on point (for additional rotational stiffness term due to change in moment arm) + + M3, A3, B3, K3 = self.sys.pointList[PointID-1].getDynamicMatrices(omegas, S_zeta, depth, kbot, cbot, r_dynamic_init=r_dynamic_init) # local 3D dynamic matrices of the point + + # following are from functions translateMatrix3to6 + H = getH(r) + K[:3,:3] += K3 + K[:3,3:] += np.matmul(K3, H) # only add up one off-diagonal sub-matrix for now, then we'll mirror at the end + K[3:,3:] += -np.matmul(getH(f3), H) - np.matmul(H, np.matmul(K3,H)) # updated 2023-05-02 + + # Inertia, added mass, and damping matrices are similar. Except for the lower right submatrix, which does not have the first term because + # the moment arm does not explictly change with the velocity or acceleration of the body. + M[:3,:3] += M3 + M[:3,3:] += np.matmul(M3, H) + M[3:,3:] += - np.matmul(H, np.matmul(M3,H)) + + A[:3,:3] += A3 + A[:3,3:] += np.matmul(A3, H) + A[3:,3:] += - np.matmul(H, np.matmul(A3,H)) + + B[:3,:3] += B3 + B[:3,3:] += np.matmul(B3, H) + B[3:,3:] += - np.matmul(H, np.matmul(B3,H)) + + K[3:,:3] = K[:3,3:].T # copy over other off-diagonal sub-matrix + M[3:,:3] = M[:3,3:].T + A[3:,:3] = A[:3,3:].T + B[3:,:3] = B[:3,3:].T + # body's own stiffness components + if lines_only == False: + + # rotational stiffness effect of weight + rCG_rotated = rotatePosition(self.rCG, self.r6[3:]) # relative position of CG about body ref point in unrotated reference frame + Kw = -np.matmul( getH([0,0, -self.m*self.sys.g]) , getH(rCG_rotated) ) + + # rotational stiffness effect of buoyancy at metacenter + rM_rotated = rotatePosition(self.rM, self.r6[3:]) # relative position of metacenter about body ref point in unrotated reference frame + Kb = -np.matmul( getH([0,0, self.sys.rho*self.sys.g*self.v]) , getH(rM_rotated) ) + + # hydrostatic heave stiffness (if AWP is nonzero) + Kwp = self.sys.rho*self.sys.g*self.AWP + + K[3:,3:] += Kw + Kb + K[2 ,2 ] += Kwp + + # Return stiffness matrix + if all_DOFs: + return M, A, B, K + else: # only return rows/columns of active DOFs + return M[:,self.DOFs][self.DOFs,:], A[:,self.DOFs][self.DOFs,:], B[:,self.DOFs][self.DOFs,:], K[:,self.DOFs][self.DOFs,:] + def draw(self, ax): '''Draws the reference axis of the body diff --git a/moorpy/line.py b/moorpy/line.py index ffd917c..e77c6c2 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -1016,6 +1016,63 @@ def getDynamicMatrices(self, *args, **kwargs): '''Compute M,A,B,K matrices for Line. See get_dynamic_matrices().''' return get_dynamic_matrices(self, *args, **kwargs) + def getDynamicMatricesLumped(self, *args, **kwargs): + '''Lump M,A,B,K matrices for Line at its extremities, returning 6x6 matrices''' + + def lump_matrix(matrix, nodes2remove=None): + # The matrices should be symmetrical, but they can be slightly off due to numerical errors. + # Because we are going to invert them twice, we force them to be symmetrical to avoid amplifying the errors. + matrix = (matrix + matrix.T)/2 + zeros = np.zeros((3,3)) + + # Remove rows and columns corresponding to the dofs of nodes2remove. Each node has 3 dofs + if nodes2remove.size==0: + nodes2remove = None + + if nodes2remove is not None: + # Convert nodes2remove to dofs + dofs2remove = np.array([(node*3, node*3+1, node*3+2) for node in nodes2remove]).flatten() + + # Create a mask for the dofs to keep + mask = np.ones(matrix.shape[0], dtype=bool) + mask[dofs2remove] = False + # Remove the rows and columns + matrix = matrix[mask][:, mask] + + matrix_inv = np.linalg.pinv(matrix) + + top_left = matrix_inv[:3, :3] + top_right = matrix_inv[:3, -3:] + bottom_left = matrix_inv[-3:, :3] + bottom_right = matrix_inv[-3:, -3:] + + # If we are not removing the extremities, we fill the whole 6x6 matrix_inv_coupled + if nodes2remove is None or (nodes2remove[0] !=0 and nodes2remove[-1] != self.nNodes-1): + matrix_inv_coupled = np.block([[top_left, top_right], [bottom_left, bottom_right]]) + + # if we are removing the first node, we fill the bottom right 3x3 matrix + if nodes2remove is not None and nodes2remove[0] == 0: + matrix_inv_coupled = np.block([[np.zeros((3,3)), np.zeros((3,3))], [np.zeros((3,3)), bottom_right]]) + + # if we are removing the last node, we fill the top left 3x3 matrix + if nodes2remove is not None and nodes2remove[-1] == self.nNodes-1: + matrix_inv_coupled = np.block([[top_left, np.zeros((3,3))], [np.zeros((3,3)), np.zeros((3,3))]]) + + return np.linalg.pinv(matrix_inv_coupled) + + # Remove the nodes that are lying on the seabed + X_mean,Y_mean,Z_mean,T_mean = self.getLineCoords(0.0,n=self.nNodes) # coordinates of line nodes and tension values + idx2remove = np.where(Z_mean <= -self.sys.depth+1e-06)[0] + + M, A, B, K, _, _ = self.getDynamicMatrices(*args, **kwargs) + Ml = lump_matrix(M, nodes2remove=idx2remove) + Al = lump_matrix(A, nodes2remove=idx2remove) + Bl = lump_matrix(B, nodes2remove=idx2remove) + Kl = lump_matrix(K, nodes2remove=idx2remove) + return Ml, Al, Bl, Kl + + + def dynamicSolve(self, *args, **kwargs): '''Compute complex amplitudes of line nodes. See get_dynamic_tension().''' diff --git a/moorpy/point.py b/moorpy/point.py index 358439a..031e11c 100644 --- a/moorpy/point.py +++ b/moorpy/point.py @@ -341,54 +341,33 @@ def getStiffnessA(self, lines_only=False, xyz=False): else: # otherwise only return rows/columns of active DOFs return K[:,self.DOFs][self.DOFs,:] - def getDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot): + def getDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dynamic_init=None): '''Gets inertia, added mass, damping, and stiffness matrices of Point due only to mooring lines (with other objects fixed) using a lumped mass model. Returns ------- - K : matrix - 3x3 analytic stiffness matrix. - ''' - - #print("Getting Point "+str(self.number)+" analytic stiffness matrix...") - - def lump_matrix(matrix, endB): - # The matrices should be symmetrical, but they can be slightly off due to numerical errors. - # Because we are going to invert them twice, we force them to be symmetrical to avoid amplifying the errors. - matrix = (matrix + matrix.T)/2 - if endB == 1: - matrix_woAnc = matrix[3:, 3:] # For some reason that I don't know yet, we need to remove the elements in the opposite side of the matrix. Otherwise it doesn't work - matrix_inv = np.linalg.pinv(matrix_woAnc) - matrix_inv_coupled = matrix_inv[-3:, -3:] - else: - matrix_woAnc = matrix[:-3, :-3] - matrix_inv = np.linalg.pinv(matrix_woAnc) - matrix_inv_coupled = matrix_inv[0:3, 0:3] - matrix_out = np.linalg.inv(matrix_inv_coupled) - return matrix_out - M = np.zeros([3,3]) A = np.zeros([3,3]) B = np.zeros([3,3]) K = np.zeros([3,3]) - + # append the stiffness matrix of each line attached to the point for lineID,endB in zip(self.attached,self.attachedEndB): line = self.sys.lineList[lineID-1] if r_dynamic_init == None: r_init = np.ones((len(omegas),line.nNodes,3)) - M_all, A_all, B_all, K_all, _, _ = line.getDynamicMatrices(omegas,S_zeta,r_init,depth,kbot,cbot) + + M_all, A_all, B_all, K_all = line.getDynamicMatricesLumped(omegas,S_zeta,r_init,depth,kbot,cbot) - M += lump_matrix(M_all, endB) - A += lump_matrix(A_all, endB) - B += lump_matrix(B_all, endB) - K += lump_matrix(K_all, endB) + M += M_all[-3:,-3:] if endB == 1 else M_all[:3, :3] + A += A_all[-3:,-3:] if endB == 1 else A_all[:3, :3] + B += B_all[-3:,-3:] if endB == 1 else B_all[:3, :3] + K += K_all[-3:,-3:] if endB == 1 else K_all[:3, :3] return M, A, B, K - def getCost(self): '''Fill in and returns a cost dictionary for this Point object. So far it only applies for if the point is an anchor.''' diff --git a/moorpy/system.py b/moorpy/system.py index e831c02..1443bb1 100644 --- a/moorpy/system.py +++ b/moorpy/system.py @@ -2877,9 +2877,15 @@ def getSystemStiffnessA(self, DOFtype="free", lines_only=False, rho=1025, g=9.81 return K - def getCoupledDynamicMatrices(self, omegas, S_zeta, depth, kbot, cbot, r_dynamic_init=None, lines_only=False): + def getCoupledDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dynamic_init=None, lines_only=False): '''Write something here later - ''' + ''' + # The methods used to get the coupled dynamic matrices do not work well if we treat each line section separately - I don't know why yet + # Hence, we group the lines into subsystems + + + + self.nDOF, self.nCpldDOF, DOFtypes = self.getDOFs() n = self.nDOF + self.nCpldDOF @@ -2888,13 +2894,13 @@ def getCoupledDynamicMatrices(self, omegas, S_zeta, depth, kbot, cbot, r_dynamic print("Getting mooring system stiffness matrix...") # get full system stiffness matrix - M_all, A_all, B_all, K_all = self.getSystemDynamicMatrices(omegas, S_zeta, r_dynamic_init, depth, kbot, cbot, DOFtype="both", lines_only=lines_only) + M_all, A_all, B_all, K_all = self.getSystemDynamicMatrices(omegas, S_zeta, depth, kbot, cbot, DOFtype="both", lines_only=lines_only, r_dynamic_init=r_dynamic_init) # invert matrix - M_inv_all = np.linalg.inv(M_all) - A_inv_all = np.linalg.inv(A_all) - B_inv_all = np.linalg.inv(B_all) - K_inv_all = np.linalg.inv(K_all) + M_inv_all = np.linalg.pinv(M_all) + A_inv_all = np.linalg.pinv(A_all) + B_inv_all = np.linalg.pinv(B_all) + K_inv_all = np.linalg.pinv(K_all) # remove free DOFs (this corresponds to saying that the sum of forces on these DOFs will remain zero) #indices = list(range(n)) # list of DOF indices that will remain active for this step @@ -2911,18 +2917,15 @@ def getCoupledDynamicMatrices(self, omegas, S_zeta, depth, kbot, cbot, r_dynamic K_inv_coupled = K_inv_all[mask,:][:,mask] # invert reduced matrix to get coupled stiffness matrix (with free DOFs assumed to equilibrate linearly) - M_coupled = np.linalg.inv(M_inv_coupled) - A_coupled = np.linalg.inv(A_inv_coupled) - B_coupled = np.linalg.inv(B_inv_coupled) - K_coupled = np.linalg.inv(K_inv_coupled) - - #if tensions: - # return K_coupled, J - #else: + M_coupled = np.linalg.pinv(M_inv_coupled) + A_coupled = np.linalg.pinv(A_inv_coupled) + B_coupled = np.linalg.pinv(B_inv_coupled) + K_coupled = np.linalg.pinv(K_inv_coupled) + return M_coupled, A_coupled, B_coupled, K_coupled - def getSystemDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot, DOFtype="free", lines_only=False, rho=1025, g=9.81): + def getSystemDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, DOFtype="free", lines_only=False, r_dynamic_init=None, rho=1025, g=9.81): '''Write something here later ''' @@ -2973,9 +2976,11 @@ def getSystemDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, #i = (body1.number-1)*6 # start counting index for body DOFs based on body number to keep indexing consistent # get body's self-stiffness matrix (now only cross-coupling terms will be handled on a line-by-line basis) - K6 = body1.getStiffnessA(lines_only=lines_only) - K[i:i+body1.nDOF, i:i+body1.nDOF] += K6 - + M6, A6, B6, K6 = body1.getDynamicMatrices(omegas, S_zeta, depth, kbot=kbot, cbot=cbot, r_dynamic_init=r_dynamic_init, lines_only=lines_only) + M[i:i+body1.nDOF, i:i+body1.nDOF] += M6 + A[i:i+body1.nDOF, i:i+body1.nDOF] += A6 + B[i:i+body1.nDOF, i:i+body1.nDOF] += B6 + K[i:i+body1.nDOF, i:i+body1.nDOF] += K6 # go through each attached point for pointID1,rPointRel1 in zip(body1.attachedP,body1.rPointRel): @@ -2984,17 +2989,20 @@ def getSystemDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, r1 = rotatePosition(rPointRel1, body1.r6[3:]) # relative position of Point about body ref point in unrotated reference frame H1 = getH(r1) # produce alternator matrix of current point's relative position to current body - for lineID in point1.attached: # go through each attached line to the Point, looking for when its other end is attached to something that moves endFound = 0 # simple flag to indicate when the other end's attachment has been found j = i + body1.nDOF # first index of the DOFs this line is attached to. Start it off at the next spot after body1's DOFs # get cross-coupling stiffness of line: force on end attached to body1 due to motion of other end - if point1.attachedEndB == 1: - KB = self.lineList[lineID-1].KBA - else: - KB = self.lineList[lineID-1].KBA.T + if r_dynamic_init == None: + r_init = np.ones((len(omegas),self.lineList[lineID-1].nNodes,3)) + M_all, A_all, B_all, K_all = self.lineList[lineID-1].getDynamicMatricesLumped(omegas,S_zeta,r_init,depth,kbot,cbot) + + MB = M_all[-3:,:3] if point1.attachedEndB == 1 else M_all[:3,-3:] + AB = A_all[-3:,:3] if point1.attachedEndB == 1 else A_all[:3,-3:] + BB = B_all[-3:,:3] if point1.attachedEndB == 1 else B_all[:3,-3:] + KB = K_all[-3:,:3] if point1.attachedEndB == 1 else K_all[:3,-3:] # look through Bodies further on in the list (coupling with earlier Bodies will already have been taken care of) for body2 in self.bodyList[self.bodyList.index(body1)+1: ]: @@ -3060,12 +3068,56 @@ def getSystemDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, # >>> TODO: handle case of free end point resting on seabed <<< # get point's self-stiffness matrix - M1, A1, B1, K1 = point.getDynamicMatrices(omegas, S_zeta, r_dynamic_init, depth, kbot, cbot) + M1, A1, B1, K1 = point.getDynamicMatrices(omegas, S_zeta, depth, kbot, cbot, r_dynamic_init=r_dynamic_init) M[i:i+n,i:i+n] += M1 A[i:i+n,i:i+n] += A1 B[i:i+n,i:i+n] += B1 K[i:i+n,i:i+n] += K1 + + for lineID in point.attached: + + j = i + n + + # go through movable points to see if one is attached + for point2 in self.pointList[self.pointList.index(point)+1: ]: + if point2.type in d: + if lineID in point2.attached: # if this point is at the other end of the line + + # get cross-coupling stiffness of line: force on end attached to point1 due to motion of other end + # The names are somewhat confusing, but AB: Added mass matrix (A) at point B due to motions of the other extremity + if r_dynamic_init == None: + r_init = np.ones((len(omegas), self.lineList[lineID-1].nNodes, 3)) + M_all, A_all, B_all, K_all = self.lineList[lineID-1].getDynamicMatricesLumped(omegas,S_zeta,r_init,depth,kbot,cbot) + MB = M_all[-3:,:3] if point.attachedEndB == 1 else M_all[:3,-3:] + AB = A_all[-3:,:3] if point.attachedEndB == 1 else A_all[:3,-3:] + BB = B_all[-3:,:3] if point.attachedEndB == 1 else B_all[:3,-3:] + KB = K_all[-3:,:3] if point.attachedEndB == 1 else K_all[:3,-3:] + + # Trim stiffness matrix to only use the enabled DOFs of each point + MB = MB[point.DOFs,:][:,point2.DOFs] + AB = AB[point.DOFs,:][:,point2.DOFs] + BB = BB[point.DOFs,:][:,point2.DOFs] + KB = KB[point.DOFs,:][:,point2.DOFs] + + # Include in total inertia matrix + M[i:i+n , j:j+point2.nDOF] += MB # force on P1 due to movement of P2 + M[j:j+point2.nDOF, i:i+n ] += MB.T # mirror (f on P2 due to x of P1) + + # Include in total added mass matrix + A[i:i+n , j:j+point2.nDOF] += AB # force on P1 due to movement of P2 + A[j:j+point2.nDOF, i:i+n ] += AB.T # mirror (f on P2 due to x of P1) + + # Include in total damping matrix + B[i:i+n , j:j+point2.nDOF] += BB # force on P1 due to movement of P2 + B[j:j+point2.nDOF, i:i+n ] += BB.T # mirror (f on P2 due to x of P1) + + # Include in total stiffness matrix + K[i:i+n , j:j+point2.nDOF] += KB # force on P1 due to movement of P2 + K[j:j+point2.nDOF, i:i+n ] += KB.T # mirror (f on P2 due to x of P1) + + j += point2.nDOF # if this point has DOFs we're considering, then count them + i += n return M, A, B, K From 016b11955f0b7104a5f6e5580cb8d2dcf4b1efd6 Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Wed, 26 Jun 2024 16:04:20 -0600 Subject: [PATCH 07/14] RAFT needs tension amplitude at nodes + removed nodes on the seabed There is still something weird with the lumped mass approach regarding the nodes on the seabed. I'm still working on it, but for now we get the best results if we neglect those nodes when lumping the matrices at the extremities of the line. --- moorpy/helpers.py | 12 ++++++------ moorpy/line.py | 30 ++++++++++++++++++++++++++++++ moorpy/system.py | 16 ++-------------- 3 files changed, 38 insertions(+), 20 deletions(-) diff --git a/moorpy/helpers.py b/moorpy/helpers.py index 03b789e..10d4fb3 100644 --- a/moorpy/helpers.py +++ b/moorpy/helpers.py @@ -1813,16 +1813,16 @@ def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_to eps_segs = np.abs(dL_dynamic)/L_static T_segs = EA_segs * eps_segs - T_nodes = np.zeros((len(omegas),N)) - T_nodes[:,0] = T_segs[:,0] - T_nodes[:,1:-1] = (T_segs[:,:-1] + T_segs[:,1:])/2 - T_nodes[:,-1] = T_segs[:,-1] + T_nodes_amp = np.zeros((len(omegas),N)) + T_nodes_amp[:,0] = T_segs[:,0] + T_nodes_amp[:,1:-1] = (T_segs[:,:-1] + T_segs[:,1:])/2 + T_nodes_amp[:,-1] = T_segs[:,-1] # S_T = np.zeros((len(omegas),N)) # S_T[:,1:] = T_e**2/dw[:,None] # S_T[:,0] = S_T[:,1] - T_nodes_psd = T_nodes**2/dw[:,None] + T_nodes_psd = T_nodes_amp**2/dw[:,None] T_nodes_std = np.sqrt(np.trapz(T_nodes_psd,omegas,axis=0)) @@ -1831,7 +1831,7 @@ def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_to s = np.zeros_like(T_nodes_std) s = np.cumsum(ds) - return T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X + return T_nodes_amp, T_nodes_psd,T_nodes_std, s,r_static,r_dynamic,r_total,X def get_modes(line,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view = False,kbot=3E+06,cbot=3E+05,seabed_tol=1E-04): diff --git a/moorpy/line.py b/moorpy/line.py index e77c6c2..4ed44be 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -1033,6 +1033,9 @@ def lump_matrix(matrix, nodes2remove=None): # Convert nodes2remove to dofs dofs2remove = np.array([(node*3, node*3+1, node*3+2) for node in nodes2remove]).flatten() + # Create a mask for the dofs to keep + # dofs2remove = np.array([node*3 + 2 for node in nodes2remove]) + # Create a mask for the dofs to keep mask = np.ones(matrix.shape[0], dtype=bool) mask[dofs2remove] = False @@ -1060,10 +1063,37 @@ def lump_matrix(matrix, nodes2remove=None): return np.linalg.pinv(matrix_inv_coupled) + + # # Get list of points + # pointList = self.pointList if hasattr(self, 'pointList') else self.sys.pointList # Difference between line and subsystem + # rList = [p.r for p in pointList] + + # # Check if end A is fixed + # idxA = next((i for i, x in enumerate(rList) if np.allclose(x, self.rA, atol=1e-8)), -1) # Find the index of the point A in the point list + # removeA = False + # if pointList[idxA].type == 1: + # removeA = True + + # # Check if end B is fixed + # idxB = next((i for i, x in enumerate(rList) if np.allclose(x, self.rB, atol=1e-8)), -1) + # removeB = False + # if pointList[idxB].type == 1: + # removeB = True + + # # Cannot remove both ends. Return matrices with zeros if so. + # if removeA and removeB: + # return M1, A1, B1, K1 + # Remove the nodes that are lying on the seabed X_mean,Y_mean,Z_mean,T_mean = self.getLineCoords(0.0,n=self.nNodes) # coordinates of line nodes and tension values idx2remove = np.where(Z_mean <= -self.sys.depth+1e-06)[0] + # Keep one of the nodes to remove, which is the one in the interface with the nodes to keep + # is_seabed = Z_mean <= -self.sys.depth+1e-06 + # transition_indices = np.where(np.abs(np.diff(is_seabed.astype(int))) == 1)[0] # Find all indices where is_seabed changes from True to False + # idx2remove_all = np.where(is_seabed)[0] # Get all indices where is_seabed is True + # idx2remove = np.setdiff1d(idx2remove_all, transition_indices) # Remove transition_indices from idx2remove_all + M, A, B, K, _, _ = self.getDynamicMatrices(*args, **kwargs) Ml = lump_matrix(M, nodes2remove=idx2remove) Al = lump_matrix(A, nodes2remove=idx2remove) diff --git a/moorpy/system.py b/moorpy/system.py index 1443bb1..9125faf 100644 --- a/moorpy/system.py +++ b/moorpy/system.py @@ -2879,13 +2879,7 @@ def getSystemStiffnessA(self, DOFtype="free", lines_only=False, rho=1025, g=9.81 def getCoupledDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dynamic_init=None, lines_only=False): '''Write something here later - ''' - # The methods used to get the coupled dynamic matrices do not work well if we treat each line section separately - I don't know why yet - # Hence, we group the lines into subsystems - - - - + ''' self.nDOF, self.nCpldDOF, DOFtypes = self.getDOFs() n = self.nDOF + self.nCpldDOF @@ -2927,13 +2921,7 @@ def getCoupledDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dyn def getSystemDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, DOFtype="free", lines_only=False, r_dynamic_init=None, rho=1025, g=9.81): '''Write something here later - ''' - - # note: This is missing some pieces, and needs to check more. - # So far this seems to not capture yaw stiffness for non-bridle configs... - # it would require proper use of chain rule for the derivatives - - + ''' # find the total number of free and coupled DOFs in case any object types changed self.nDOF, self.nCpldDOF, _ = self.getDOFs() From c9ed16b88bc3c4aa0a758c4b4ca768ac4a87d229 Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Thu, 11 Jul 2024 15:10:29 -0600 Subject: [PATCH 08/14] Bug fix to remove nodes on the sea bed The previous code to remove nodes on the sea bed wasn't working properly for subsystem --- moorpy/line.py | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/moorpy/line.py b/moorpy/line.py index 2d12316..f9b8bac 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -1026,8 +1026,8 @@ def lump_matrix(matrix, nodes2remove=None): matrix = (matrix + matrix.T)/2 zeros = np.zeros((3,3)) - # Remove rows and columns corresponding to the dofs of nodes2remove. Each node has 3 dofs - if nodes2remove.size==0: + # Empty nodes2remove is the same thing as not removing any nodes + if nodes2remove is not None and nodes2remove.size==0: nodes2remove = None if nodes2remove is not None: @@ -1086,8 +1086,32 @@ def lump_matrix(matrix, nodes2remove=None): # return M1, A1, B1, K1 # Remove the nodes that are lying on the seabed - X_mean,Y_mean,Z_mean,T_mean = self.getLineCoords(0.0,n=self.nNodes) # coordinates of line nodes and tension values - idx2remove = np.where(Z_mean <= -self.sys.depth+1e-06)[0] + # X_mean,Y_mean,Z_mean,T_mean = self.getLineCoords(0.0,n=self.nNodes) # coordinates of line nodes and tension values + # idx2remove = np.where(Z_mean <= -self.sys.depth+1e-06)[0] + + if not hasattr(self, 'lineList'): + X_mean,Y_mean,Z_mean,T_mean = self.getLineCoords(0.0,n=self.nNodes) # coordinates of line nodes and tension values + depth = self.sys.depth + else: + X_mean, Y_mean, Z_mean = [np.zeros(self.nNodes) for _ in range(3)] + depth = self.depth + idxNodeA = 0 + for iline, line in enumerate(self.lineList): + x, y, z, _ = line.getLineCoords(0.0) + n = len(x) + if iline == 0: # For the first line, include all nodes + X_mean[idxNodeA:idxNodeA+n] = x + Y_mean[idxNodeA:idxNodeA+n] = y + Z_mean[idxNodeA:idxNodeA+n] = z + idxNodeA += n + else: # For subsequent lines, exclude the first node to avoid repetition + X_mean[idxNodeA:idxNodeA+n-1] = x[1:] + Y_mean[idxNodeA:idxNodeA+n-1] = y[1:] + Z_mean[idxNodeA:idxNodeA+n-1] = z[1:] + idxNodeA += (n - 1) + idx2remove = np.where(Z_mean <= -depth+1e-06)[0] + # idx2remove=None + # Keep one of the nodes to remove, which is the one in the interface with the nodes to keep # is_seabed = Z_mean <= -self.sys.depth+1e-06 From b286b7f889c66d9a0930cfe3896f5e6c1f43a721 Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Thu, 11 Jul 2024 15:14:50 -0600 Subject: [PATCH 09/14] Added nSegs in makeGeneric + initialize subsystem in system.initalize - Added number of segments (nSegs) in subsystem.makeGeneric to make sure that we are using the correct number of nodes (i.e., the one specified by the user) in the lumped mass approach - Initialize any subsystem that is within a system in system.initialize --- moorpy/subsystem.py | 12 +++++++++--- moorpy/system.py | 2 ++ 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/moorpy/subsystem.py b/moorpy/subsystem.py index f77f515..859e845 100644 --- a/moorpy/subsystem.py +++ b/moorpy/subsystem.py @@ -150,7 +150,7 @@ def initialize(self, daf_dict={}): self.damage = np.zeros(self.nNodes) - def makeGeneric(self, lengths, types, connectors=[], suspended=0): + def makeGeneric(self, lengths, types, connectors=[], suspended=0, nSegs=40): '''Creates a cable of n components going between an anchor point and a floating body (or a bridle point). If shared, it goes between two @@ -223,9 +223,15 @@ def makeGeneric(self, lengths, types, connectors=[], suspended=0): else: raise Exception(f"Can't find lineType '{types[i]}' in the SubSystem.") - # add the line segment using the reference to its lineType dict # add the line segment using the reference to its lineType dict - self.addLine(lengths[i],self.lineTypes[i]) + if nSegs is None: + self.addLine(lengths[i], self.lineTypes[i]) + elif isinstance(nSegs, (int, float)): + self.addLine(lengths[i], self.lineTypes[i], nSegs=nSegs) + elif isinstance(nSegs, list): + self.addLine(lengths[i], self.lineTypes[i], nSegs=nSegs[i]) + else: + raise ValueError("Invalid type for nSegs. Expected None, a number, or a list.") # add the upper end point of the segment if i==self.nLines-1: # if this is the upper-most line diff --git a/moorpy/system.py b/moorpy/system.py index 8634a4e..5b90c04 100644 --- a/moorpy/system.py +++ b/moorpy/system.py @@ -1510,6 +1510,8 @@ def initialize(self, plots=0): for line in self.lineList: line.staticSolve(profiles=1) # flag to enable additional line outputs used for plotting, tension results, etc. + if hasattr(line, 'lineList'): + line.initialize() for point in self.pointList: point.getForces() From 7ee3fe4c9cf48d95735c3852d51f1596830b958e Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Tue, 16 Jul 2024 09:12:28 -0600 Subject: [PATCH 10/14] Cleaned parts of the code + removed *args and **kwargs Removed *args and **kwargs to make the inputs more explicit --- moorpy/body.py | 3 +-- moorpy/line.py | 37 ++++++------------------------------- moorpy/subsystem.py | 5 +---- 3 files changed, 8 insertions(+), 37 deletions(-) diff --git a/moorpy/body.py b/moorpy/body.py index b4e8dc6..c112ed7 100644 --- a/moorpy/body.py +++ b/moorpy/body.py @@ -397,8 +397,7 @@ def getDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dynamic_in Returns ------- - K : matrix - 6x6 analytic stiffness matrix. + M, A, B, K: 6x6 inertia, added mass, damping, and stiffness matrices. ''' M, A, B, K = (np.zeros([6,6]) for _ in range(4)) diff --git a/moorpy/line.py b/moorpy/line.py index e7a132f..711ee73 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -1015,11 +1015,11 @@ def revertToStaticStiffness(self): self.L = self.L0 - def getDynamicMatrices(self, *args, **kwargs): + def getDynamicMatrices(self, omegas, S_zeta, r_dynamic, depth, kbot=0, cbot=0, seabed_tol=1e-4): '''Compute M,A,B,K matrices for Line. See get_dynamic_matrices().''' - return get_dynamic_matrices(self, *args, **kwargs) + return get_dynamic_matrices(self, omegas, S_zeta, r_dynamic, depth, kbot, cbot, seabed_tol=1e-4) - def getDynamicMatricesLumped(self, *args, **kwargs): + def getDynamicMatricesLumped(self, omegas, S_zeta, r_dynamic, depth, kbot=0, cbot=0, seabed_tol=1e-4): '''Lump M,A,B,K matrices for Line at its extremities, returning 6x6 matrices''' def lump_matrix(matrix, nodes2remove=None): @@ -1065,32 +1065,8 @@ def lump_matrix(matrix, nodes2remove=None): matrix_inv_coupled = np.block([[top_left, np.zeros((3,3))], [np.zeros((3,3)), np.zeros((3,3))]]) return np.linalg.pinv(matrix_inv_coupled) - - - # # Get list of points - # pointList = self.pointList if hasattr(self, 'pointList') else self.sys.pointList # Difference between line and subsystem - # rList = [p.r for p in pointList] - - # # Check if end A is fixed - # idxA = next((i for i, x in enumerate(rList) if np.allclose(x, self.rA, atol=1e-8)), -1) # Find the index of the point A in the point list - # removeA = False - # if pointList[idxA].type == 1: - # removeA = True - - # # Check if end B is fixed - # idxB = next((i for i, x in enumerate(rList) if np.allclose(x, self.rB, atol=1e-8)), -1) - # removeB = False - # if pointList[idxB].type == 1: - # removeB = True - - # # Cannot remove both ends. Return matrices with zeros if so. - # if removeA and removeB: - # return M1, A1, B1, K1 - + # Remove the nodes that are lying on the seabed - # X_mean,Y_mean,Z_mean,T_mean = self.getLineCoords(0.0,n=self.nNodes) # coordinates of line nodes and tension values - # idx2remove = np.where(Z_mean <= -self.sys.depth+1e-06)[0] - if not hasattr(self, 'lineList'): X_mean,Y_mean,Z_mean,T_mean = self.getLineCoords(0.0,n=self.nNodes) # coordinates of line nodes and tension values depth = self.sys.depth @@ -1112,8 +1088,7 @@ def lump_matrix(matrix, nodes2remove=None): Z_mean[idxNodeA:idxNodeA+n-1] = z[1:] idxNodeA += (n - 1) idx2remove = np.where(Z_mean <= -depth+1e-06)[0] - # idx2remove=None - + # idx2remove = None # Keep one of the nodes to remove, which is the one in the interface with the nodes to keep # is_seabed = Z_mean <= -self.sys.depth+1e-06 @@ -1121,7 +1096,7 @@ def lump_matrix(matrix, nodes2remove=None): # idx2remove_all = np.where(is_seabed)[0] # Get all indices where is_seabed is True # idx2remove = np.setdiff1d(idx2remove_all, transition_indices) # Remove transition_indices from idx2remove_all - M, A, B, K, _, _ = self.getDynamicMatrices(*args, **kwargs) + M, A, B, K, _, _ = self.getDynamicMatrices(omegas, S_zeta, r_dynamic, depth, kbot, cbot, seabed_tol=seabed_tol) Ml = lump_matrix(M, nodes2remove=idx2remove) Al = lump_matrix(A, nodes2remove=idx2remove) Bl = lump_matrix(B, nodes2remove=idx2remove) diff --git a/moorpy/subsystem.py b/moorpy/subsystem.py index b15193c..5d9d0f1 100644 --- a/moorpy/subsystem.py +++ b/moorpy/subsystem.py @@ -503,10 +503,7 @@ def revertToStaticStiffness(self): System.revertToStaticStiffness(self) - # ----- Function for dynamic frequency-domain tensions ----- - - def getDynamicMatrices(self, omegas, S_zeta, r_dynamic, depth, kbot, cbot, - seabed_tol=1e-4): + def getDynamicMatrices(self, omegas, S_zeta, r_dynamic, depth, kbot, cbot, seabed_tol=1e-4): '''Compute M,A,B,K matrices for the Subsystem. This calls get_dynamic_matrices() for each Line in the Subsystem then combines the results. Note that this method overrides the Line method. Other From acae46136a8a774f92d3cc83b0abb8a9d7932811 Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Wed, 17 Jul 2024 15:16:47 -0600 Subject: [PATCH 11/14] Reformulation to compute dynamic matrices of lines separately Previously, each call of system.getSystemDynamicMatrices() (and other equivalent functions in the body and point classes) would recompute the dynamic matrices of the line (inertia, added mass, damping, and stiffness). That required the sea state and the motions of the line nodes at each call. Now, we need to precompute the dynamic matrices of each line using line.updateLumpedMass(). The advantages are: 1. Fewer computations 2. It is simpler to integrate the lumped mass approach in RAFT because now it's easier to account for motions of the fairlead Also, the only matrix that requires the sea state and motions of the line is the damping matrix. So this commit reformulates the system, body, and line usage of the lumped mass approach to be able to provide inertia, added mass, and stiffness matrices without knowledge of the sea state or motions of the line. I will add some examples, but you just need to call system.updateSystemDynamicMatrices() and then system.getSystemDynamicMatrices() --- moorpy/body.py | 4 +- moorpy/helpers.py | 30 ++++++++++----- moorpy/line.py | 98 ++++++++++++++++++++++++++++++++++++----------- moorpy/point.py | 8 ++-- moorpy/system.py | 40 +++++++++++++------ 5 files changed, 130 insertions(+), 50 deletions(-) diff --git a/moorpy/body.py b/moorpy/body.py index c112ed7..4044598 100644 --- a/moorpy/body.py +++ b/moorpy/body.py @@ -384,7 +384,7 @@ def getStiffnessA(self, lines_only=False, all_DOFs=False): return K[:,self.DOFs][self.DOFs,:] - def getDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dynamic_init=None, lines_only=False, all_DOFs=False): + def getDynamicMatrices(self, lines_only=False, all_DOFs=False): '''Gets the dynamic matrices of the Body with other objects fixed using a lumped mass approach. @@ -407,7 +407,7 @@ def getDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dynamic_in r = rotatePosition(rPointRel, self.r6[3:]) # relative position of Point about body ref point in unrotated reference frame f3 = self.sys.pointList[PointID-1].getForces() # total force on point (for additional rotational stiffness term due to change in moment arm) - M3, A3, B3, K3 = self.sys.pointList[PointID-1].getDynamicMatrices(omegas, S_zeta, depth, kbot, cbot, r_dynamic_init=r_dynamic_init) # local 3D dynamic matrices of the point + M3, A3, B3, K3 = self.sys.pointList[PointID-1].getDynamicMatrices() # local 3D dynamic matrices of the point # following are from functions translateMatrix3to6 H = getH(r) diff --git a/moorpy/helpers.py b/moorpy/helpers.py index 02a49e0..1f97f28 100644 --- a/moorpy/helpers.py +++ b/moorpy/helpers.py @@ -1879,14 +1879,14 @@ def get_dynamic_matrices(Line, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_t return M,A,B,K,r_mean,EA_segs -def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_tol=1e-4,tol = 0.01,iters=100, w = 0.8, conv_time=True): +def get_dynamic_tension(line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_tol=1e-4,tol = 0.01,iters=100, w = 0.8, conv_time=False, returnMatrices=False): """ - Evaluates dynamic tension at all the nodes for an instance of MoorPy's Line or CompositeLine classes. + Evaluates dynamic tension at all the nodes for an instance of MoorPy's line or CompositeLine classes. Parameters ---------- - Line : Line/CompositeLine - An instance of MoorPy's Line or CompositeLine classes. + line : line/CompositeLine + An instance of MoorPy's line or CompositeLine classes. omegas : ndarray Array of frequencies in rad/s. S_zeta : ndarray @@ -1910,8 +1910,10 @@ def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_to w : float, optional Succesive relaxation coefficient, by default 0.8 - Returns + Returns (if returnMatrices = False) ------- + T_nodes_psd: ndarray + Tension amplitude at nodes given as (m,n) array where m is the number of frequencies and n is the number of nodes. T_nodes_psd: ndarray Tension spectra at nodes given as (m,n) array where m is the number of frequencies and n is the number of nodes. T_nodes_std: ndarray @@ -1926,8 +1928,15 @@ def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_to Combined static and dynamic positions. X: ndarray Solution of the dynamic problem. + + Returns (if returnMatrices = True) + ------- + M: (n,n) inertia matrix, where n is the number of nodes + A: (n,n) added mass matrix + B: (n,n) damping matrix + K: (n,n) stiffness matrix """ - N = Line.nNodes + N = line.nNodes n_dofs = 3*N if np.all(RAO_A == 0): @@ -1937,7 +1946,7 @@ def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_to # intialize iteration matrices r_dynamic_init = np.ones((len(omegas),N,3)) - M,A,B,K,r_static,EA_segs = Line.getDynamicMatrices(omegas,S_zeta,r_dynamic_init,depth,kbot,cbot,seabed_tol=seabed_tol) # TODO: return EA_segs + M,A,B,K,r_static,EA_segs = line.getDynamicMatrices(omegas,S_zeta,r_dynamic_init,depth,kbot,cbot,seabed_tol=seabed_tol) # TODO: return EA_segs X = np.zeros((len(omegas),n_dofs),dtype = 'complex') r_dynamic = np.zeros(((len(omegas),int(n_dofs/3),3)),dtype = 'complex') S_Xd = np.zeros((len(omegas),n_dofs),dtype = 'float') @@ -1965,7 +1974,7 @@ def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_to break else: sigma_Xd0[:] = w * sigma_Xd + (1.-w) * sigma_Xd0 - _,_,B[:],_,_,_ = Line.getDynamicMatrices(omegas,S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=seabed_tol) + _,_,B[:],_,_,_ = line.getDynamicMatrices(omegas,S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=seabed_tol) if conv_time: print(f'Finished {ni} dynamic tension iterations in {time.time()-start} seconds (w = {w}).') @@ -2004,7 +2013,10 @@ def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_to s = np.zeros_like(T_nodes_std) s = np.cumsum(ds) - return T_nodes_amp, T_nodes_psd,T_nodes_std, s,r_static,r_dynamic,r_total,X + if returnMatrices: + return M, A, B, K + else: + return T_nodes_amp, T_nodes_psd, T_nodes_std, s, r_static, r_dynamic, r_total, X def get_modes(line,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view = False,kbot=3E+06,cbot=3E+05,seabed_tol=1E-04): diff --git a/moorpy/line.py b/moorpy/line.py index 711ee73..c619bcf 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -80,7 +80,6 @@ def __init__(self, mooringSys, num, L, lineType, nSegs=100, cb=0, isRod=0, attac self.fCurrent = np.zeros(3) # total current force vector on the line [N] - def loadData(self, dirname, rootname, sep='.MD.', id=0): '''Loads line-specific time series data from a MoorDyn output file''' @@ -1013,14 +1012,40 @@ def revertToStaticStiffness(self): # revert to original line length self.L = self.L0 - - - def getDynamicMatrices(self, omegas, S_zeta, r_dynamic, depth, kbot=0, cbot=0, seabed_tol=1e-4): - '''Compute M,A,B,K matrices for Line. See get_dynamic_matrices().''' - return get_dynamic_matrices(self, omegas, S_zeta, r_dynamic, depth, kbot, cbot, seabed_tol=1e-4) + + + + def updateLumpedMass(self, omegas, S_zeta, depth, kbot=0, cbot=0, seabed_tol=1e-4, RAO_A=None, RAO_B=None): + '''Updates the inertia, added mass, damping, and stiffness matrices for the line using a lumped mass approach. + + The updated values are stored in: + - self.Ml_allNodes, self.Al_allNodes, self.Bl_allNodes, and self.Kl_allNodes: 3N x 3N matrices, where N is the number of nodes (3 dofs for each line node) + - self.Ml, self.Al, self.Bl, and self.Kl: 6 x 6 matrices that lump the matrices _allNodes at the extremities of the line. For example, + self.Kl[:3, :3] is analogous to self.KA, self.Kl[3:, 3:] is analogous to self.KB, and self.Kl[:3, 3:] is analogous to self.KBA. + + Only the damping matrix depends on sea state and RAOs. Hence, this function can be called without the RAOs to compute the inertia, added mass, and + stiffness matrices. But do NOT trust the values of the damping matrix if the RAOs are not provided!!! + + Inputs + ---------- + omegas : ndarray + Array of frequencies in rad/s. + S_zeta : ndarray + Wave spectrum array in m^2/(rad/s), must be of the same length as omegas. + depth : float + Water depth. + kbot : float + Vertical stiffness for points lying on the seabed. + cbot : float + Vertical damping for points lying on the seabed. + seabed_tol : float, optional + Distance from seabed within which a node is considered to be lying on the seabed, by default 1e-4 m. + RAO_A : ndarray + Translational RAOs for end node A given as a (m,3) array where m is the number of frequencies (must be equal to the length of omegas) . + RAO_B : ndarray + Translational RAOs for end node B given as a (m,3) array where m is the number of frequencies (must be equal to the length of omegas) . - def getDynamicMatricesLumped(self, omegas, S_zeta, r_dynamic, depth, kbot=0, cbot=0, seabed_tol=1e-4): - '''Lump M,A,B,K matrices for Line at its extremities, returning 6x6 matrices''' + ''' def lump_matrix(matrix, nodes2remove=None): # The matrices should be symmetrical, but they can be slightly off due to numerical errors. @@ -1090,25 +1115,54 @@ def lump_matrix(matrix, nodes2remove=None): idx2remove = np.where(Z_mean <= -depth+1e-06)[0] # idx2remove = None - # Keep one of the nodes to remove, which is the one in the interface with the nodes to keep - # is_seabed = Z_mean <= -self.sys.depth+1e-06 - # transition_indices = np.where(np.abs(np.diff(is_seabed.astype(int))) == 1)[0] # Find all indices where is_seabed changes from True to False - # idx2remove_all = np.where(is_seabed)[0] # Get all indices where is_seabed is True - # idx2remove = np.setdiff1d(idx2remove_all, transition_indices) # Remove transition_indices from idx2remove_all - - M, A, B, K, _, _ = self.getDynamicMatrices(omegas, S_zeta, r_dynamic, depth, kbot, cbot, seabed_tol=seabed_tol) - Ml = lump_matrix(M, nodes2remove=idx2remove) - Al = lump_matrix(A, nodes2remove=idx2remove) - Bl = lump_matrix(B, nodes2remove=idx2remove) - Kl = lump_matrix(K, nodes2remove=idx2remove) - return Ml, Al, Bl, Kl + # # Keep two of the nodes to remove, which are the ones in the interface with the nodes to keep + # is_seabed = Z_mean <= -depth + 1e-06 + # transition_indices = np.where(np.abs(np.diff(is_seabed.astype(int))) == 1)[0] + + # # # Adjust to keep an additional node at each transition + # # # For True to False transitions, keep the transition index and the one before it + # # # For False to True transitions, keep the transition index and the one after it + # additional_indices = [] + # for idx in transition_indices: + # if is_seabed[idx]: # True to False transition + # additional_indices.append(idx - 1 if idx > 0 else idx) # Ensure idx-1 is not negative + # else: # False to True transition + # additional_indices.append(idx + 1 if idx + 1 < len(is_seabed) else idx) # Ensure idx+1 is within bounds + # transition_and_additional_indices = np.unique(transition_indices.tolist() + additional_indices) # Combine transition indices with their additional indices + # idx2remove_all = np.where(is_seabed)[0] # Get all indices where is_seabed is True + # idx2remove = np.setdiff1d(idx2remove_all, transition_and_additional_indices) # Remove transition and additional indices from idx2remove_all + + # If we have the motions of the extremities of the lines, we can iterate the quadratic drag to compute the damping matrix. + # Otherwise, we use the default values for the motion amplitude (unitary displacement for all nodes) to linearize the quadratic drag. This is clearly wrong, but affects only drag. + if RAO_A is None and RAO_B is None: + r_dynamic = np.ones((len(omegas), self.nNodes, 3)) + self.Ml_allNodes, self.Al_allNodes, self.Bl_allNodes, self.Kl_allNodes, _, _ = self.getDynamicMatrices(omegas, S_zeta, r_dynamic, depth, kbot=kbot, cbot=cbot, seabed_tol=seabed_tol) + else: + if RAO_A is None: + RAO_A = np.zeros([len(omegas), 3]) + if RAO_B is None: + RAO_B = np.zeros([len(omegas), 3]) + self.Ml_allNodes, self.Al_allNodes, self.Bl_allNodes, self.Kl_allNodes = self.dynamicSolve(omegas, S_zeta, RAO_A, RAO_B, depth, kbot, cbot, seabed_tol=seabed_tol, returnMatrices=True) # Dynamic solve calls get_dynamic_matrices internally + self.Ml = lump_matrix(self.Ml_allNodes, nodes2remove=idx2remove) + self.Al = lump_matrix(self.Al_allNodes, nodes2remove=idx2remove) + self.Bl = lump_matrix(self.Bl_allNodes, nodes2remove=idx2remove) + self.Kl = lump_matrix(self.Kl_allNodes, nodes2remove=idx2remove) + + def getDynamicMatricesLumped(self): + '''Return the lumped M,A,B,K matrices for the Line object.''' + if not hasattr(self, 'Ml'): + raise LineError("Call updateLumpedMass first") + return self.Ml, self.Al, self.Bl, self.Kl + def getDynamicMatrices(self, omegas, S_zeta, r_dynamic, depth, kbot=0, cbot=0, seabed_tol=1e-4): + '''Compute M,A,B,K matrices for Line. See get_dynamic_matrices().''' + return get_dynamic_matrices(self, omegas, S_zeta, r_dynamic, depth, kbot, cbot, seabed_tol=seabed_tol) - def dynamicSolve(self, *args, **kwargs): + def dynamicSolve(self, omegas, S_zeta, RAO_A, RAO_B, depth, kbot, cbot, seabed_tol=1e-4, tol = 0.01, iters=100, w = 0.8, conv_time=False, returnMatrices=False): '''Compute complex amplitudes of line nodes. See get_dynamic_tension().''' - return get_dynamic_tension(self, *args, **kwargs) + return get_dynamic_tension(self, omegas, S_zeta, RAO_A, RAO_B, depth, kbot, cbot, seabed_tol=seabed_tol, tol=tol, iters=iters, w=w, conv_time=conv_time, returnMatrices=returnMatrices) def getModes(self,*args, **kwargs): diff --git a/moorpy/point.py b/moorpy/point.py index 3d74383..6c0ef15 100644 --- a/moorpy/point.py +++ b/moorpy/point.py @@ -337,7 +337,7 @@ def getStiffnessA(self, lines_only=False, xyz=False): else: # otherwise only return rows/columns of active DOFs return K[:,self.DOFs][self.DOFs,:] - def getDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dynamic_init=None): + def getDynamicMatrices(self): '''Gets inertia, added mass, damping, and stiffness matrices of Point due only to mooring lines (with other objects fixed) using a lumped mass model. @@ -351,11 +351,9 @@ def getDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dynamic_in # append the stiffness matrix of each line attached to the point for lineID,endB in zip(self.attached,self.attachedEndB): - line = self.sys.lineList[lineID-1] - if r_dynamic_init == None: - r_init = np.ones((len(omegas),line.nNodes,3)) + line = self.sys.lineList[lineID-1] - M_all, A_all, B_all, K_all = line.getDynamicMatricesLumped(omegas,S_zeta,r_init,depth,kbot,cbot) + M_all, A_all, B_all, K_all = line.getDynamicMatricesLumped() M += M_all[-3:,-3:] if endB == 1 else M_all[:3, :3] A += A_all[-3:,-3:] if endB == 1 else A_all[:3, :3] diff --git a/moorpy/system.py b/moorpy/system.py index 0576037..d87f6bd 100644 --- a/moorpy/system.py +++ b/moorpy/system.py @@ -2932,7 +2932,7 @@ def getSystemStiffnessA(self, DOFtype="free", lines_only=False, rho=1025, g=9.81 return K - def getCoupledDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dynamic_init=None, lines_only=False): + def getCoupledDynamicMatrices(self, lines_only=False): '''Write something here later ''' self.nDOF, self.nCpldDOF, DOFtypes = self.getDOFs() @@ -2943,7 +2943,7 @@ def getCoupledDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dyn print("Getting mooring system stiffness matrix...") # get full system stiffness matrix - M_all, A_all, B_all, K_all = self.getSystemDynamicMatrices(omegas, S_zeta, depth, kbot, cbot, DOFtype="both", lines_only=lines_only, r_dynamic_init=r_dynamic_init) + M_all, A_all, B_all, K_all = self.getSystemDynamicMatrices(DOFtype="both", lines_only=lines_only) # invert matrix M_inv_all = np.linalg.pinv(M_all) @@ -2974,7 +2974,7 @@ def getCoupledDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, r_dyn return M_coupled, A_coupled, B_coupled, K_coupled - def getSystemDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, DOFtype="free", lines_only=False, r_dynamic_init=None, rho=1025, g=9.81): + def getSystemDynamicMatrices(self, DOFtype="free", lines_only=False): '''Write something here later ''' # find the total number of free and coupled DOFs in case any object types changed @@ -3019,7 +3019,7 @@ def getSystemDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, DOFtyp #i = (body1.number-1)*6 # start counting index for body DOFs based on body number to keep indexing consistent # get body's self-stiffness matrix (now only cross-coupling terms will be handled on a line-by-line basis) - M6, A6, B6, K6 = body1.getDynamicMatrices(omegas, S_zeta, depth, kbot=kbot, cbot=cbot, r_dynamic_init=r_dynamic_init, lines_only=lines_only) + M6, A6, B6, K6 = body1.getDynamicMatrices(lines_only=lines_only) M[i:i+body1.nDOF, i:i+body1.nDOF] += M6 A[i:i+body1.nDOF, i:i+body1.nDOF] += A6 B[i:i+body1.nDOF, i:i+body1.nDOF] += B6 @@ -3038,9 +3038,7 @@ def getSystemDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, DOFtyp j = i + body1.nDOF # first index of the DOFs this line is attached to. Start it off at the next spot after body1's DOFs # get cross-coupling stiffness of line: force on end attached to body1 due to motion of other end - if r_dynamic_init == None: - r_init = np.ones((len(omegas),self.lineList[lineID-1].nNodes,3)) - M_all, A_all, B_all, K_all = self.lineList[lineID-1].getDynamicMatricesLumped(omegas,S_zeta,r_init,depth,kbot,cbot) + M_all, A_all, B_all, K_all = self.lineList[lineID-1].getDynamicMatricesLumped() MB = M_all[-3:,:3] if point1.attachedEndB == 1 else M_all[:3,-3:] AB = A_all[-3:,:3] if point1.attachedEndB == 1 else A_all[:3,-3:] @@ -3111,7 +3109,7 @@ def getSystemDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, DOFtyp # >>> TODO: handle case of free end point resting on seabed <<< # get point's self-stiffness matrix - M1, A1, B1, K1 = point.getDynamicMatrices(omegas, S_zeta, depth, kbot, cbot, r_dynamic_init=r_dynamic_init) + M1, A1, B1, K1 = point.getDynamicMatrices() M[i:i+n,i:i+n] += M1 A[i:i+n,i:i+n] += A1 B[i:i+n,i:i+n] += B1 @@ -3128,9 +3126,7 @@ def getSystemDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, DOFtyp # get cross-coupling stiffness of line: force on end attached to point1 due to motion of other end # The names are somewhat confusing, but AB: Added mass matrix (A) at point B due to motions of the other extremity - if r_dynamic_init == None: - r_init = np.ones((len(omegas), self.lineList[lineID-1].nNodes, 3)) - M_all, A_all, B_all, K_all = self.lineList[lineID-1].getDynamicMatricesLumped(omegas,S_zeta,r_init,depth,kbot,cbot) + M_all, A_all, B_all, K_all = self.lineList[lineID-1].getDynamicMatricesLumped() MB = M_all[-3:,:3] if point.attachedEndB == 1 else M_all[:3,-3:] AB = A_all[-3:,:3] if point.attachedEndB == 1 else A_all[:3,-3:] @@ -3164,7 +3160,27 @@ def getSystemDynamicMatrices(self, omegas, S_zeta, depth, kbot=0, cbot=0, DOFtyp i += n return M, A, B, K - + + + def updateSystemDynamicMatrices(self, omegas=np.array([0.0]), S_zeta=np.array([0.0]), seabed_tol=1e-4): + '''Updates the dynamic matrices of all the lines in the system. Works only with lumped mass approach (MoorDyn input file) + This function can only properly update/compute the inertia, added mass, and stiffness matrices of each line in the system. + + This function also updates the damping matrix, but this is done considering unitary amplitude motions of the nodes, which + is not correct. + If `omegas` and `S_zeta` are provided, the damping matrix will account for wave kinematics as well, but that would still + be wrong + + To properly compute the damping matrix, you should call line.updateLumpedMass externally providing the motion RAOs of the + line ends. + ''' + kbot=0 # Should read those from the input file + cbot=0 + + for line in self.lineList: + line.updateLumpedMass(omegas, S_zeta, self.depth, kbot=kbot, cbot=cbot, seabed_tol=seabed_tol) + + def getAnchorLoads(self, sfx, sfy, sfz, N): ''' Calculates anchor loads Parameters From e4e3386b17c739a0b5756445456904e6380fdd3a Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Fri, 2 Aug 2024 08:55:38 -0600 Subject: [PATCH 12/14] Changed amplitudes to be complex and have a factor of two - The wave amplitude is now 2*Sw*dw (it was missing the factor 2) - The tension amplitude is now a complex value which is important to get the correct phase for shared moorings --- moorpy/helpers.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/moorpy/helpers.py b/moorpy/helpers.py index 1f97f28..aa9a94c 100644 --- a/moorpy/helpers.py +++ b/moorpy/helpers.py @@ -1983,7 +1983,7 @@ def get_dynamic_tension(line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_to prepend= omegas[0] - (omegas[1]-omegas[0]), append= omegas[-1] + (omegas[-1]-omegas[-2])) dw = (dw[1:]+dw[:-1])/2 - wave_amps = np.sqrt(S_zeta*dw) #evaluate wave amplitudes of harmonic components from wave spectrum + wave_amps = np.sqrt(2*S_zeta*dw) #evaluate wave amplitudes of harmonic components from wave spectrum r_dynamic *= wave_amps[:,None,None] r_total = r_static[None,:,:] + r_dynamic @@ -1992,10 +1992,10 @@ def get_dynamic_tension(line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_to tangents = dr_static/np.linalg.norm(r_static[:-1] - r_static[1:], axis=-1)[:,None] L_static = np.linalg.norm(dr_static, axis=-1) dL_dynamic = np.einsum('mni,ni->mn', dr_dynamic, tangents) - eps_segs = np.abs(dL_dynamic)/L_static + eps_segs = dL_dynamic/L_static T_segs = EA_segs * eps_segs - T_nodes_amp = np.zeros((len(omegas),N)) + T_nodes_amp = np.zeros((len(omegas),N), dtype='complex') T_nodes_amp[:,0] = T_segs[:,0] T_nodes_amp[:,1:-1] = (T_segs[:,:-1] + T_segs[:,1:])/2 T_nodes_amp[:,-1] = T_segs[:,-1] @@ -2004,7 +2004,7 @@ def get_dynamic_tension(line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_to # S_T[:,1:] = T_e**2/dw[:,None] # S_T[:,0] = S_T[:,1] - T_nodes_psd = T_nodes_amp**2/dw[:,None] + T_nodes_psd = np.abs(T_nodes_amp)**2/(2*dw[:,None]) T_nodes_std = np.sqrt(np.trapz(T_nodes_psd,omegas,axis=0)) From 75e6d39213272a36338c4d00e6013aaebac27f5b Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Fri, 2 Aug 2024 15:45:32 -0600 Subject: [PATCH 13/14] Default values for kbot and cbot --- moorpy/line.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/moorpy/line.py b/moorpy/line.py index c619bcf..303e1c6 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -1160,7 +1160,7 @@ def getDynamicMatrices(self, omegas, S_zeta, r_dynamic, depth, kbot=0, cbot=0, s '''Compute M,A,B,K matrices for Line. See get_dynamic_matrices().''' return get_dynamic_matrices(self, omegas, S_zeta, r_dynamic, depth, kbot, cbot, seabed_tol=seabed_tol) - def dynamicSolve(self, omegas, S_zeta, RAO_A, RAO_B, depth, kbot, cbot, seabed_tol=1e-4, tol = 0.01, iters=100, w = 0.8, conv_time=False, returnMatrices=False): + def dynamicSolve(self, omegas, S_zeta, RAO_A, RAO_B, depth, kbot=0, cbot=0, seabed_tol=1e-4, tol = 0.01, iters=100, w = 0.8, conv_time=False, returnMatrices=False): '''Compute complex amplitudes of line nodes. See get_dynamic_tension().''' return get_dynamic_tension(self, omegas, S_zeta, RAO_A, RAO_B, depth, kbot, cbot, seabed_tol=seabed_tol, tol=tol, iters=iters, w=w, conv_time=conv_time, returnMatrices=returnMatrices) From 27d86593e17f012b6fa0152b38dcfa2f07e054db Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Fri, 2 Aug 2024 15:46:13 -0600 Subject: [PATCH 14/14] Example to show how to use the lumped mass approach This example is a modified/simplified version of one of the tests we made. As this example is meant to show how to use this new functionality only, I didn't check the results of the example case thoroughly. --- examples/RAO_fl.pkl | Bin 0 -> 3352 bytes examples/lumped_mass.py | 74 +++++++++++ examples/volturn_chain.dat | 248 +++++++++++++++++++++++++++++++++++++ 3 files changed, 322 insertions(+) create mode 100644 examples/RAO_fl.pkl create mode 100644 examples/lumped_mass.py create mode 100644 examples/volturn_chain.dat diff --git a/examples/RAO_fl.pkl b/examples/RAO_fl.pkl new file mode 100644 index 0000000000000000000000000000000000000000..0002ae9e54404be63f2f3501b764e1d9136ea6e0 GIT binary patch literal 3352 zcmXxnc|4Tc9|v%T5wcUH)I^tQ6KNsL_n>4;C?Z0kR0t8nElVLxrJ8PbDhg55ja&W9 z38m~>C_;9!WT|8czixBx_mAhF=XuU~pU<50eLWM(e@#@B=XbFddHH*KGgVzq_`0cj z`g`~tJ?QIukjZ9=IQhD{obWp7=j-p{$G*yDQGY+kW-ah?88AX7MWuu32A&)9w{h9(z>5!r{)7 zWfr5L@w>`k=t1@B;(KARYEfm{6EQ{Lp?K6K$&ZouD}a7rS%$=1K3*5Ioks#m0ehLQ zUSIGeaMP#<1@VUOLe}@;r<=@6&Pq)X@zZM4(H-#2ITMYi2sZypN%Ho5Y{rsc#p1od(3>%BKvkp8w_?Ol*H z)Qnb`$Xb*T-xv6x%1X=%@`4xqS6IxNSd4G%VAjtU5#K1p{n#IuS$GILzES?Z(mVlQ zd^}q6HuMo0e-xhBn5^F0>4WKN>iI|6X;}VJ*)g}`zlnbg6qZbL{_YFF)DW#WMZ?>8 zg65(d{h@%&KMu{g+jkcxFtMoj%qg*lxi~W1DQ?B(T(Z8ekgk!d;W(rP+ven6=G)!E zW}(@q%{Sa8zHeYCqF-|GLI8x8?w5}~9D_qx>LHHqH;C^WG%_DjhECdJ<@x6(Gkg4T zTzw@{#Lge#5uwL#;xGFU-wfo|s+>lL41rOVGG5?u0tQsg z&3*We65mfyyv3Oj;E~4m70SFuL1AFuJJQ|o#hLhif<-cOrp3<EqiMSOFh?6b18NX8gv?{YR?8gmX7ug=cU z`Nxp>=Har^c7cU!OhHLq+$AsJ47?6?=;$!eBR(GF-z0QtY#$9*{3&@MJn&C=OI`B% zq`U?hp9giH?N10wQv{=A|I7lbbMPv{%DuIcPW(KmhSnphC(y;H4bJ+W6BQ2oB0taZ z7|W6Qc~Q7$QcG3Z3P{*gtuk>s9M+$2n%bEpMb^iQv|~F>Vr>D$9v` zbAktR`yU0FPZhmu@3F>1vv7SZ zC#k`qhr#V%6x7yzPTkdQe5;D z0CS=*>%|y@>*GfT7i@N_%gY)i)3Sh2}A8B6pT=eLCJ7?nas7c3$ z0MNOg>QkySz~TC+=tqD>TJ!@2Tt1VirWbPr-ZYp^OL69jkBSV#qGy_(baEtHZiWRI z9fZ4YKiYUif8%iT6F@;4#f3Yg_`vC#&A#4Jdq_=UiSr2$bGSYMBnjusZf)2C7Zo2V zi2T?C6WR1;U;jChp8%@Y759+oY-0>i*Rh9nwm}j5n(_N&n)o?qTRe7uDW!K!z#Uk zN^m#8fn`pEOAdDupCGChxO4NoohqDl45xZG>*3POvDvk@lf);8ymvgiYyBjaA(*k% zDRIskQ<8PH6JF;NpAZr&h%8)!P#9Kn zZr&5HmK>cbY%@lDq9|F2zD_b#2B!whxT<7rz)ufo8f7ebLH3_0Dr9KH)+c`D8iJhPYjLc%Nwq+lmM4f*BXz$_3&7!=R#a#KZCpf#L)J+sE90QaeOmh z;e)rMI!F|*oN|$fBR(w5X=-qe?a4X9iqMIaW-3;z;vzxaaD<9h@x( zx>{%-b-`7B@|CB-Es~!&k`>yS-QOz(J052NmfFtXo(~DMIr*K?)Ddx< zWA?1;ix$8s2hWM({ARL$Bv6fPsX?b|E5~v~^r-xR7R0Fu49*FP3> z_*-ACPq)*7+?Mi~jPO2^j|8e`&L^=N#ekM(xB2wU7O1G06FPKmoUHHH{EX>$CMKht zjS*R$67SW3Ua-S(_vbv~lSFbkEzfH+8#q-=ji+bsD#J(Xn;y~)|M?`*QvtKEpK(;2 zk~EzAcM%OrK1&@?Kif-ul1N+I;whtC6nC@6CkGk&V0iAc{XrH zzbz}A2_NKWcuy&#`lyg9kfHIh))VxCs7tCTYchC%~c;RxlPA>i|I0#PyBbjq|lcn zdRX1ARu2Eyi)sDOmt*C&3|XhEtt3Bb)O5hrR6?ndqe$x=#)C8*FJR=^n^#Zvk2Gq8 zTIsuOZH%=|es(GE7GZaBo0<4(7V$}=uq$6{nQPiPZEdGjYL7|4pCw9xvx~10pEQcN z)I+mKo@6NGL>tRT3PIe`E0cK|i6lR1ROcS7zFGJYN8hNHxyMifLTY_tSIFj*{AAF$ g&&6dT9;1w){0xU=Qw2D_vF&3mr= 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() \ No newline at end of file diff --git a/examples/volturn_chain.dat b/examples/volturn_chain.dat new file mode 100644 index 0000000..ff9d413 --- /dev/null +++ b/examples/volturn_chain.dat @@ -0,0 +1,248 @@ +--------------------- MoorDyn Input File ------------------------------------ +MoorDyn input file +----------------------- LINE TYPES ------------------------------------------ +TypeName Diam Mass/m EA BA/-zeta EI Cd Ca CdAx CaAx +(name) (m) (kg/m) (N) (N-s/-) (N-m^2) (-) (-) (-) (-) +chain 3.33000e-01 6.85000e+02 3.27000e+09 -1.00000e+00 0 1.11000e+00 8.20000e-01 2.00000e-01 2.70000e-01 +---------------------- POINT PROPERTIES -------------------------------- +ID Type X Y Z Mass Volume CdA Ca +(#) (-) (m) (m) (m) (kg) (m^3) (m^2) (-) +1 Vessel -5.80000e+01 0.00000e+00 -1.40000e+01 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 +2 Fixed -8.37600e+02 0.00000e+00 -2.00000e+02 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 +---------------------- LINES ---------------------------------------- +ID LineType AttachA AttachB UnstrLen NumSegs LineOutputs +(#) (name) (#) (#) (m) (-) (-) +1 chain 2 1 8.50000e+02 50 - +---------------------- OPTIONS ----------------------------------------- +2 writeLog Write a log file +0.001 dtM time step to use in mooring integration (s) +3.0e6 kBot bottom stiffness (Pa/m) +3.0e5 cBot bottom damping (Pa-s/m) +1025.0 WtrDnsty water density (kg/m^3) +200.0 WtrDpth water depth (m) +1.0 dtIC time interval for analyzing convergence during IC gen (s) +60.0 TmaxIC max time for ic gen (s) +4.0 CdScaleIC factor by which to scale drag coefficients during dynamic relaxation (-) +0.001 threshIC threshold for IC convergence (-) +------------------------ OUTPUTS -------------------------------------------- +FAIRTEN1 +POINT1PX +POINT1PY +POINT1PZ +POINT1FX +POINT1FY +POINT1FZ +POINT1T +POINT2PX +POINT2PY +POINT2PZ +POINT2FX +POINT2FY +POINT2FZ +POINT2T +LINE1N0PX +LINE1N0PY +LINE1N0PZ +LINE1N0T +LINE1N1PX +LINE1N1PY +LINE1N1PZ +LINE1N1T +LINE1N2PX +LINE1N2PY +LINE1N2PZ +LINE1N2T +LINE1N3PX +LINE1N3PY +LINE1N3PZ +LINE1N3T +LINE1N4PX +LINE1N4PY +LINE1N4PZ +LINE1N4T +LINE1N5PX +LINE1N5PY +LINE1N5PZ +LINE1N5T +LINE1N6PX +LINE1N6PY +LINE1N6PZ +LINE1N6T +LINE1N7PX +LINE1N7PY +LINE1N7PZ +LINE1N7T +LINE1N8PX +LINE1N8PY +LINE1N8PZ +LINE1N8T +LINE1N9PX +LINE1N9PY +LINE1N9PZ +LINE1N9T +LINE1N10PX +LINE1N10PY +LINE1N10PZ +LINE1N10T +LINE1N11PX +LINE1N11PY +LINE1N11PZ +LINE1N11T +LINE1N12PX +LINE1N12PY +LINE1N12PZ +LINE1N12T +LINE1N13PX +LINE1N13PY +LINE1N13PZ +LINE1N13T +LINE1N14PX +LINE1N14PY +LINE1N14PZ +LINE1N14T +LINE1N15PX +LINE1N15PY +LINE1N15PZ +LINE1N15T +LINE1N16PX +LINE1N16PY +LINE1N16PZ +LINE1N16T +LINE1N17PX +LINE1N17PY +LINE1N17PZ +LINE1N17T +LINE1N18PX +LINE1N18PY +LINE1N18PZ +LINE1N18T +LINE1N19PX +LINE1N19PY +LINE1N19PZ +LINE1N19T +LINE1N20PX +LINE1N20PY +LINE1N20PZ +LINE1N20T +LINE1N21PX +LINE1N21PY +LINE1N21PZ +LINE1N21T +LINE1N22PX +LINE1N22PY +LINE1N22PZ +LINE1N22T +LINE1N23PX +LINE1N23PY +LINE1N23PZ +LINE1N23T +LINE1N24PX +LINE1N24PY +LINE1N24PZ +LINE1N24T +LINE1N25PX +LINE1N25PY +LINE1N25PZ +LINE1N25T +LINE1N26PX +LINE1N26PY +LINE1N26PZ +LINE1N26T +LINE1N27PX +LINE1N27PY +LINE1N27PZ +LINE1N27T +LINE1N28PX +LINE1N28PY +LINE1N28PZ +LINE1N28T +LINE1N29PX +LINE1N29PY +LINE1N29PZ +LINE1N29T +LINE1N30PX +LINE1N30PY +LINE1N30PZ +LINE1N30T +LINE1N31PX +LINE1N31PY +LINE1N31PZ +LINE1N31T +LINE1N32PX +LINE1N32PY +LINE1N32PZ +LINE1N32T +LINE1N33PX +LINE1N33PY +LINE1N33PZ +LINE1N33T +LINE1N34PX +LINE1N34PY +LINE1N34PZ +LINE1N34T +LINE1N35PX +LINE1N35PY +LINE1N35PZ +LINE1N35T +LINE1N36PX +LINE1N36PY +LINE1N36PZ +LINE1N36T +LINE1N37PX +LINE1N37PY +LINE1N37PZ +LINE1N37T +LINE1N38PX +LINE1N38PY +LINE1N38PZ +LINE1N38T +LINE1N39PX +LINE1N39PY +LINE1N39PZ +LINE1N39T +LINE1N40PX +LINE1N40PY +LINE1N40PZ +LINE1N40T +LINE1N41PX +LINE1N41PY +LINE1N41PZ +LINE1N41T +LINE1N42PX +LINE1N42PY +LINE1N42PZ +LINE1N42T +LINE1N43PX +LINE1N43PY +LINE1N43PZ +LINE1N43T +LINE1N44PX +LINE1N44PY +LINE1N44PZ +LINE1N44T +LINE1N45PX +LINE1N45PY +LINE1N45PZ +LINE1N45T +LINE1N46PX +LINE1N46PY +LINE1N46PZ +LINE1N46T +LINE1N47PX +LINE1N47PY +LINE1N47PZ +LINE1N47T +LINE1N48PX +LINE1N48PY +LINE1N48PZ +LINE1N48T +LINE1N49PX +LINE1N49PY +LINE1N49PZ +LINE1N49T +LINE1N50PX +LINE1N50PY +LINE1N50PZ +LINE1N50T +END +------------------------- need this line -------------------------------------- \ No newline at end of file