Skip to content

Commit

Permalink
Merge pull request #36 from NREL/slope_fix
Browse files Browse the repository at this point in the history 
Slope fix
  • Loading branch information
@shousner
shousner authored Dec 26, 2024
2 parents 91b088f + f55b264 commit 0940ecc
Show file tree
Hide file tree
Showing 3 changed files with 234 additions and 75 deletions.
177 changes: 119 additions & 58 deletions moorpy/Catenary.py
View file
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
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

0 comments on commit 0940ecc

Please sign in to comment.