Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Slope fix #36

Merged
merged 5 commits into from
Dec 26, 2024
Merged

Slope fix #36

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
3867ccf
Preliminary support for lines out of the water:
mattEhall Nov 20, 2024
3b1b85b
More robust seabed slope support
mattEhall Nov 20, 2024
969381d
Support for lines laying fully along sloped seabed
mattEhall Nov 22, 2024
8973cfc
Line and Subsystem shadows now plot on seabed
mattEhall Dec 23, 2024
f55b264
Subsystem bathymetry functionality
shousner Nov 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
177 changes: 119 additions & 58 deletions moorpy/Catenary.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -158,6 +158,7 @@ def catenary(XF, ZF, L, EA, W, CB=0, alpha=0, HF0=0, VF0=0, Tol=0.000001,
# ProfileType=5: Similar to above but both ends are off seabed, so it's U shaped and fully slack
# ProfileType=6: Completely vertical line without seabed contact (on the seabed is handled by 4 and 5)
# ProfileType = 7: Portion of the line is resting on the seabed, and the seabed has a slope
# ProfileType = 8: line resting on the seabed, and the seabed has a slope

EA_W = EA/W

Expand All @@ -181,68 +182,120 @@ def dV_dZ_s(z0, H): # height off seabed to evaluate at (infinite if 0), horizo
# because a large value here risks adding a bad cross coupling term in the system stiffness matrix

# ProfileType 0 case - entirely along seabed
if ZF==0.0 and hA == 0.0 and W > 0:
if hB==0.0 and hA == 0.0 and W > 0:

# >>> this case should be extended to support sloped seabeds <<<
if not alpha == 0: raise CatenaryError("Line along seabed but seabed is sloped - not yet supported")

ProfileType = 0
if alpha == 0: # flat seabed case

if CB==0 or XF <= L: # case 1: no friction, or zero tension
HF = np.max([0, (XF/L - 1.0)*EA])
HA = 1.0*HF
elif 0.5*L + EA/CB/W*(1-XF/L) <= 0: # case 2: seabed friction but tension at anchor (xB estimate < 0)
HF = (XF/L -1.0)*EA + 0.5*CB*W*L
HA = np.max([0.0, HF - CB*W*L])
else: # case 3: seabed friction and zero anchor tension
if EA*CB*W*(XF-L) < 0: # something went wrong
breakpoint()
HF = np.sqrt(2*EA*CB*W*(XF-L))
HA = 0.0
ProfileType = 0

VF = 0.0
VA = 0.0

if HF > 0: # if taut
dHF_dXF = EA/L # approximation <<< what about friction? <<<<<<<<
#dVF_dZF = W + HF/L # vertical stiffness <<< approximation a
dVF_dZF = dV_dZ_s(Tol, HF) # vertical stiffness <<< approximation b
else: # if slack
dHF_dXF = 0.0
dVF_dZF = W # vertical stiffness

info["HF"] = HF # solution to be used to start next call (these are the solved variables, may be for anchor if line is reversed)
info["VF"] = 0.0
info["stiffnessB"] = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]])
info["stiffnessA"] = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]])
info["stiffnessBA"] = np.array([[-dHF_dXF, 0.0], [0.0, 0.0]])
info["LBot"] = L
info['ProfileType'] = 0
info['Zextreme'] = 0

if plots > 0:
if CB==0 or XF <= L: # case 1: no friction, or zero tension
HF = np.max([0, (XF/L - 1.0)*EA])
HA = 1.0*HF
elif 0.5*L + EA/CB/W*(1-XF/L) <= 0: # case 2: seabed friction but tension at anchor (xB estimate < 0)
HF = (XF/L -1.0)*EA + 0.5*CB*W*L
HA = np.max([0.0, HF - CB*W*L])
else: # case 3: seabed friction and zero anchor tension
if EA*CB*W*(XF-L) < 0: # something went wrong
breakpoint()
HF = np.sqrt(2*EA*CB*W*(XF-L))
HA = 0.0

VF = 0.0
VA = 0.0

if HF > 0: # if taut
dHF_dXF = EA/L # approximation <<< what about friction? <<<<<<<<
#dVF_dZF = W + HF/L # vertical stiffness <<< approximation a
dVF_dZF = dV_dZ_s(Tol, HF) # vertical stiffness <<< approximation b
else: # if slack
dHF_dXF = 0.0
dVF_dZF = W # vertical stiffness

info["HF"] = HF # solution to be used to start next call (these are the solved variables, may be for anchor if line is reversed)
info["VF"] = 0.0
info["stiffnessB"] = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]])
info["stiffnessA"] = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]])
info["stiffnessBA"] = np.array([[-dHF_dXF, 0.0], [0.0, 0.0]])
info["LBot"] = L
info['ProfileType'] = 0
info['Zextreme'] = 0

if plots > 0:

if CB > 0 and XF > L:
xB = L - HF/W/CB # location of point at which line tension reaches zero
else:
xB = 0.0

# z values remain zero in this case

if CB==0 or XF <= L: # case 1: no friction, or zero tension
Xs = XF/L*s # X values uniformly distributed
Te = Te + np.max([0, (XF/L - 1.0)*EA]) # uniform tension
elif xB <= 0: # case 2: seabed friction but tension at anchor
Xs = s*(1+CB*W/EA*(0.5*s-xB))
Te = HF + CB*W*(s-L)
else: # case 3: seabed friction and zero anchor tension
for I in range(nNodes):
if s[I] <= xB: # if this node is in the zero tension range
Xs[I] = s[I]; # x is unstretched, z and Te remain zero

else: # the tension is nonzero
Xs[I] = s[I] + CB*W/EA*(s[I] - xB)**2
Te[I] = HF - CB*W*(L-s[I])

if CB > 0 and XF > L:
xB = L - HF/W/CB # location of point at which line tension reaches zero
else:
xB = 0.0
else: # sloped seabed case with full contact

ProfileType = 8

# z values remain zero in this case
cos_alpha = np.cos(np.radians(alpha))
tan_alpha = np.tan(np.radians(alpha))
sin_alpha = np.sin(np.radians(alpha))

if CB==0 or XF <= L: # case 1: no friction, or zero tension
Xs = XF/L*s # X values uniformly distributed
Te = Te + np.max([0, (XF/L - 1.0)*EA]) # uniform tension
elif xB <= 0: # case 2: seabed friction but tension at anchor
Xs = s*(1+CB*W/EA*(0.5*s-xB))
Te = HF + CB*W*(s-L)
else: # case 3: seabed friction and zero anchor tension
for I in range(nNodes):
if s[I] <= xB: # if this node is in the zero tension range
Xs[I] = s[I]; # x is unstretched, z and Te remain zero

else: # the tension is nonzero
Xs[I] = s[I] + CB*W/EA*(s[I] - xB)**2
Te[I] = HF - CB*W*(L-s[I])
LAB = XF/cos_alpha # calling this the lenght along the seabed for this case
if LAB < 0:
breakpoint()
LAB = max(LAB, 0) # Ensure LAB is not less than zero


Ts = np.max([0, (LAB/L - 1.0)*EA]) # approximate tension based on streatch
Tg = L*W*sin_alpha # effect of weight and slope
# ignoring friction for now

HF = cos_alpha*(Ts + 0.5*Tg)
VF = sin_alpha*(Ts + 0.5*Tg)
HA = cos_alpha*(Ts - 0.5*Tg)
VA = sin_alpha*(Ts - 0.5*Tg)

if plots > 0:

for I in range(nNodes):
Xs[I] = np.min([s[I]*cos_alpha, XF])
Zs[I] = Xs[I]*tan_alpha
Te[I] = Ts - 0.5*Tg + Tg*float(I/(nNodes-1))

if HF > 0: # if taut
dHF_dXF = EA/L # CRUDE approximation
dVF_dZF = dV_dZ_s(Tol, HF) # same
else: # if slack
dHF_dXF = 0.0
dVF_dZF = W # vertical stiffness

info["HF"] = HF
info["VF"] = VF
info["stiffnessB"] = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]])
info["stiffnessA"] = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]])
info["stiffnessBA"] = np.array([[-dHF_dXF, 0.0], [0.0, 0.0]])
info["LBot"] = L
info['ProfileType'] = 8
info['Zextreme'] = 0

if plots > 0:

if CB > 0 and XF > L:
xB = L - HF/W/CB # location of point at which line tension reaches zero
else:
xB = 0.0


# ProfileType 4 case - fully slack
Expand Down Expand Up @@ -1187,7 +1240,7 @@ def eval_func_cat(X, args):
info['error'] = True
info['message'] = "ProfileType 7: VF_HF + SQRT1VF_HF2 <= 0"

elif HF*1000000 < VF:
elif HF*1000000000 < VF:
info['error'] = True
info['message'] = "ProfileType 7: HF << VF, line is slack, not supported yet"
breakpoint()
Expand Down Expand Up @@ -1467,7 +1520,15 @@ def step_func_cat(X, args, Y, info, Ytarget, err, tols, iter, maxIter):
#(fAH1, fAV1, fBH1, fBV1, info1) = catenary(0.16525209064463553, 1000.1008645356192, 976.86359774357, 861890783.955385, -5.663702119364548, 0.0, -0.0, 0.0, 17325264.3383085, 5.000000000000001e-05, 101, 1, 1000)
#(fAH1, fAV1, fBH1, fBV1, info1) = catenary(0.0, 997.8909672113768, 1006.0699998926483, 361256000.00000006, 2.4952615384615333, 0, 0.0, 0, 0, 0.0001, 101, 100, 1, 1000)
#(fAH1, fAV1, fBH1, fBV1, info1) = catenary(148.60942473375343, 1315.624259105325, 1279.7911844766932, 17497567581.802254, -81.8284437736437, CB=0.0, alpha=-0.0, HF0=795840.9590915733, VF0=6399974.444658389, Tol=0.0005, MaxIter=100, depth=1300.0, plots=1)
(fAH1, fAV1, fBH1, fBV1, info1) = catenary(2887.113193120885, -838.1577017360617, 2979.7592390255295, 41074520109.87981, 1104.878355564403, CB=0.0, alpha=81.41243787249468, HF0=0.0, VF0=3564574.0622291695, Tol=0.0005, MaxIter=100, depth=3000.0, plots=1)
#(fAH1, fAV1, fBH1, fBV1, info1) = catenary(2887.113193120885, -838.1577017360617, 2979.7592390255295, 41074520109.87981, 1104.878355564403, CB=0.0, alpha=81.41243787249468, HF0=0.0, VF0=3564574.0622291695, Tol=0.0005, MaxIter=100, depth=3000.0, plots=1)


#(fAH1, fAV1, fBH1, fBV1, info1) = catenary(20.386176076554875, -83.33919733899663, 85.0, 469622033.24936056, -166.5252156894419,
# CB=-897.8894086817654, alpha=0, HF0=1045436.0319202082, VF0=4280849.602990574, Tol=0.0001, MaxIter=100, depth=800, plots=1)


(fAH1, fAV1, fBH1, fBV1, info1) = catenary(242.43063215088046, 16.438515754164882, 242.94180148966598, 1977564385.9456, 3941.831324315529, CB=0.0, alpha=3.879122219600495, HF0=400446.3427158061, VF0=26452.686481212237, Tol=0.0001, MaxIter=100, depth=820.2155229016184, plots=1)
#(fAH1, fAV1, fBH1, fBV1, info1) = catenary(113.98763202381406, 135.64003587014417, 222.9088982113214, 1799620189.0374997, 3587.1279256290677, CB=0.0, alpha=-1.6131140432335964, HF0=737382.1987996304, VF0=1014826.8890259801, Tol=0.0001, MaxIter=100, depth=914.6515127916663, plots=1)

#(fAH1, fAV1, fBH1, fBV1, info1) = catenary(274.9, 15.6, 328.5, 528887323., -121.,
# CB=-48., alpha=0, HF0=17939., VF0=20596., Tol=2e-05, MaxIter=100, plots=1)
Expand Down
53 changes: 40 additions & 13 deletions moorpy/line.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -564,7 +564,16 @@ def drawLine(self, Time, ax, color="k", endpoints=False, shadow=True, colortensi
linebit.append(ax.plot(Xs, Ys, Zs, color=color, lw=lw, zorder=100))

if shadow:
ax.plot(Xs, Ys, np.zeros_like(Xs)-self.sys.depth, color=[0.5, 0.5, 0.5, 0.2], lw=lw, zorder = 1.5) # draw shadow
if self.sys.seabedMod == 0:
Zs = np.zeros_like(Xs)-self.sys.depth
elif self.seabedMod == 1:
Zs = self.sys.depth - self.sys.xSlope*Xs - self.sys.ySlope*Ys
elif self.seabedMod == 2:
Zs = np.zeros(len(Xs))
for i in range(len(Xs)):
Zs[i] = self.sys.getDepthAtLocation(Xs[i], Ys[i])

ax.plot(Xs, Ys, Zs, color=[0.5, 0.5, 0.5, 0.2], lw=lw, zorder = 1.5) # draw shadow

if endpoints == True:
linebit.append(ax.scatter([Xs[0], Xs[-1]], [Ys[0], Ys[-1]], [Zs[0], Zs[-1]], color = color))
Expand Down Expand Up @@ -743,6 +752,21 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0):
else:
self.cb = -depthA - self.rA[2] # when cb < 0, -cb is defined as height of end A off seabed (in catenary)

# Handle case of line above the water
if self.rA[2] > 0 or self.rB[2] > 0: # end B out of water << or do this in catenary?
ww = (self.type['m'] - np.pi/4*self.type['d_vol']**2 *self.sys.rho) * self.sys.g
wa = (self.type['m']) * self.sys.g # weight in air

z_top = max(self.rA[2], self.rB[2])
z_bot = min(self.rA[2], self.rB[2])

# overwrite the weight based on the wet-dry weight ratio based on end z vals
self.w = ww + (wa - ww) * ( -min(0, z_bot) )/( z_top - z_bot )
else:
# reset to the normal weight if the line hasn't left the water
self.w = (self.type['m'] - np.pi/4*self.type['d_vol']**2 *self.sys.rho) * self.sys.g

# >>> TODO: should replace all this with a more versatile catenary solve in future <<<

# ----- Perform rotation/transformation to 2D plane of catenary -----

Expand All @@ -752,9 +776,9 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0):
if np.sum(np.abs(self.fCurrent)) > 0:

# total line exernal force per unit length vector (weight plus current drag)
w_vec = self.fCurrent/self.L + np.array([0, 0, -self.type["w"]])
w = np.linalg.norm(w_vec)
w_hat = w_vec/w
w_vec = self.fCurrent/self.L + np.array([0, 0, -self.w])
w_total = np.linalg.norm(w_vec)
w_hat = w_vec/w_total

# >>> may need to adjust to handle case of buoyant vertical lines <<<

Expand All @@ -773,7 +797,12 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0):
# if no current force, things are simple
else:
R_curr = np.eye(3,3)
w = self.type["w"]
w_total = self.w


# horizontal and vertical dimensions of line profile (end A to B)
LH = np.linalg.norm(dr[:2])
LV = dr[2]


# apply a rotation about Z' to align the line profile with the X'-Z' plane
Expand All @@ -787,8 +816,10 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0):
if self.rA[2] <= -depthA or self.rB[2] <= -depthB:
nvecA_prime = np.matmul(R, nvecA)

dz_dx = -nvecA_prime[0]*(1.0/nvecA_prime[2]) # seabed slope components
dz_dy = -nvecA_prime[1]*(1.0/nvecA_prime[2]) # seabed slope components
#dz_dx = -nvecA_prime[0]*(1.0/nvecA_prime[2]) # seabed slope components
#dz_dy = -nvecA_prime[1]*(1.0/nvecA_prime[2]) # seabed slope components
# Seabed slope along line direction (based on end A/B depth)
dz_dx = (-depthB + depthA)/LH
# we only care about dz_dx since the line is in the X-Z plane in this rotated situation
alpha = np.degrees(np.arctan(dz_dx))
cb = self.cb
Expand All @@ -800,25 +831,21 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0):
alpha = 0
cb = self.cb

# horizontal and vertical dimensions of line profile (end A to B)
LH = np.linalg.norm(dr[:2])
LV = dr[2]


# ----- call catenary function or alternative and save results -----

#If EA is found in the line properties we will run the original catenary function
if 'EA' in self.type:
try:
(fAH, fAV, fBH, fBV, info) = catenary(LH, LV, self.L, self.EA,
w, CB=cb, alpha=alpha, HF0=self.HF, VF0=self.VF, Tol=tol,
w_total, CB=cb, alpha=alpha, HF0=self.HF, VF0=self.VF, Tol=tol,
nNodes=self.nNodes, plots=profiles, depth=self.sys.depth)

except CatenaryError as error:
raise LineError(self.number, error.message)
#If EA isnt found then we will use the ten-str relationship defined in the input file
else:
(fAH, fAV, fBH, fBV, info) = nonlinear(LH, LV, self.L, self.type['Str'], self.type['Ten'],np.linalg.norm(w))
(fAH, fAV, fBH, fBV, info) = nonlinear(LH, LV, self.L, self.type['Str'], self.type['Ten'],np.linalg.norm(w_total))


# save line profile coordinates in global frame (involves inverse rotation)
Expand Down
Loading
Loading