VIV
While looking up my 20year 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 DNVGL’s RPC205:
Excitation  Vortex shedding resonance (lockin) 
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 crossflow VIV lockin 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 nonshaded part.
The way to read the above graph is pretty simple:
 Unshaded area (or area in white background) is the safe area.
 The intersection of the curve (for each pipe size of D × t) with the upper boundary of the nonshaded 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.
 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: utf8 *
"""
VIV check based on DNVRPC205
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
# Added mass (kg/m):
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)
# Inline 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, # 1year 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; simplysupported: 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 down^{1} in the user input section of the code above, I use this following to generate D
array:
#!/usr/bin/env python
# * coding: utf8 *
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 * 1E3, 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

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