ckunte.net

# VIV

While looking up my 20-year old spreadsheet, I realised to my horror that it contained a systematic error in a conditional statement that was failing silently! Here’s the condition I was checking — based on DNV-GL’s RP-C205:

 Excitation Vortex shedding resonance (lock-in) Inline 1.0 ≤ Vr ≤ 4.5, and Ks ≤ 1.8 Crossflow 3.0 ≤ Vr ≤ 16.0

The formula I had in Excel was like this:

=IF((1.0 <= Vr <= 4.5), "VIV occurs", "OK")


One actually needs to use the AND operator to get Excel to recognise the second conditional statement:

=IF(AND(Vr >= 1.0, Vr <= 4.5), "VIV occurs", "OK")


Whereas python can handle this:

if(1.0 <= Vr <= 4.5):
print "VIV occurs"
else:
print "OK"


as well as this:

if(Vr >= 1.0 and Vr <= 4.5):
print "VIV occurs"
else:
print "OK"


Anyway, I digress. It just so happened that I spotted this while doing some VIV checks by hand and I thought I’d rewrite this whole check in python, along with a range of pipe diameters — for computing and graphing in one go.

While the upper limit for cross-flow VIV lock-in is 16, I’ve set the plot below to the max. value of 6.0 — to keep the curvatures of plots more readable, since in any case, the safe area is only the non-shaded part.

The way to read the above graph is pretty simple:

1. Unshaded area (or area in white background) is the safe area.
2. The intersection of the curve (for each pipe size of D × t) with the upper boundary of the non-shaded area (i.e., Vr = 1.0) indicates the max. possible pipe length (span) that is unlikely to experience VIV from (ocean) current. For example, 762OD × 22WT pipe can span up to 40m, while a 406OD × 22WT only up to 22m before VIV occurs.
3. The shaded overlap is a zone in which both inline as well as crossflow VIV excitations occur. (See range of occurrences in the Table above.)

Here’s the code:

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
VIV check based on DNV-RP-C205
2019 ckunte

Jul 19: Initial commit
"""
import numpy as np
import matplotlib.pyplot as plt

def curr_viv(x, V, tm, Ym, f, c):
# Hydrodynamic properties:
if tm > 0:
cm = 1.2 # inertia coeff. for rough cylinder
else:
cm = 1.6 # inertia coeff. for smooth cylinder
for i in D:
# Pipe cross sectional area (m^2):
A = (np.pi / 4.0) * (i**2 - (i - 2 * t)**2)
# Pipe cross sectional moment of inertia (m^4):
I = (np.pi / 64.0) * (i**4 - (i - 2 * t)**4)
# Mass of pipe (kg/m):
Ms = A * Ys
Ma = cm * rho * np.pi * (i + 2 * tm)**2 / 4.0
# Entrained mass (kg/m):
Mi = f * rho * (np.pi / 4) * (i - 2 * t)**2
# Marine growth mass (kg/m):
Mg = np.pi * (i + tm) * tm * Ym
# Total mass (kg/m):
Mtot = Ms + Ma + Mi + Mg
# Pipe natural frequency (Roark):
fn = (0.5 * c / np.pi) * (E * I / (Mtot * x**4))**0.5
# Stability parameter (Ks):
Ks = 2 * Mtot * (2 * np.pi * beta) / (rho * (i + 2 * tm)**2)
# Reduced velocity:
Vr = V / (fn * i)
# Plotting results
lbl = '%0.0f$\\times$ %0.0f (D/t=%.1f), Ks=%.1f' %((i * 1E3), (t * 1E3), (i / t), Ks)
plt.plot(x, Vr, label=lbl, linewidth=2)
plt.title('Vortex induced vibration check for %.1fm/s current' %V)
# In-line VIV occurrence limits ( 1.0 =< Vr =< 4.5 )
plt.axhspan(1.0, 4.5, facecolor='r', alpha=0.18)
# Cross flow VIV occurrence limits (3.0 =< Vr =< 16.0 )
plt.axhspan(3.0, 6.0, facecolor='orange', alpha=0.18) # set to 6.0 for clarity
plt.grid(True)
plt.legend(loc=0, prop={'size': 12})
plt.xlabel('Length of pipe, L (m)')
plt.ylabel('Reduced velocity, Vr')
plt.savefig('vivc.png')
pass

# -- Constants --
E = 2.05E+11 # steel modulus of elasticity (N/m^2)
Ys = 7850.0  # steel density (kg/m^3)
beta = 0.05  # structural damping ratio for slender elements in water (pure steel pipes)
rho = 1025.0 # seawater density (kg/m^3)

# -- User inputs --
# Pipe length:
x = np.arange(10, 50.0, 0.05) # (m)
# Pipe sizes (D and t) to test for VIV
D = [0.4064, 0.4572, 0.508, 0.5588, 0.6096, 0.6604, 0.762] # (16"--30")
t = 0.022
# Data: V, tm, Ym, f, c
d = [
0.7,   # 1-year current velocity, V (m/s)
0.133, # Marine growth thickness, tm (m)
575.0,   # Marine growth density, Ym (kg/m^3)
1,     # Pipe flooded, f (1 -- flooded; 0 -- buoyant)
15.4    # End c => fixed: 22.2; clamped: 15.4; simply-supported: 9.87; cantilevered: 3.52
]
# -- End of user inputs --

if __name__ == '__main__':
curr_viv(x, d[0], d[1], d[2], d[3], d[4])
pass


When I’m lazy to write the various pipe diameters down1 in the user input section of the code above, I use this following to generate D array:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np

# Pipe outer diameter range (in inches)
D_ini = 16           # Initial dia. (only input required)
D_fin = D_ini * 2    # End of range
D = []               # Empty array to populate

# Gen a list between D_ini and D_fin incremented by 2"
for i in np.arange(D_ini, D_fin, 2):
D_metric = round(i * 25.4 * 1E-3, 4)
D.append(D_metric)

print "D =", D


This generates:

D = [0.4064, 0.4572, 0.508, 0.5588, 0.6096, 0.6604, 0.7112, 0.762]


I’ve come to appreciate that using numpy syntax is consistent between python and python3, which use xrange and range respectively, whereas in numpy, it’s numpy.arange in both, with the only difference being that the print function in python3 is print("D =", D) in the above example. Even simpler option is with seq <initval> <incr> <endval>:

seq 0.4064 0.0508 0.7620


This generates:

0.4064
0.4572
0.508
0.5588
0.6096
0.6604
0.7112
0.762


1. Imperial units, esp. round figures are often convenient for riser diameters.