Modelling Ocean Thermal Energy Conversion pipe’s deployment from a barge

The below model is a part of Ocean Intella software — open-source software for modelling offshore pipelaying.

According to Wikipedia, Ocean thermal energy conversion (OTEC) is a renewable energy technology that harnesses the temperature difference between the warm surface waters of the ocean and the cold depths to run a heat engine to produce electricity. It is a unique form of clean energy generation that has the potential to provide a consistent and sustainable source of power. OTEC uses the ocean thermal gradient between cooler deep and warmer shallow or surface seawaters to run a heat engine and produce useful work, usually in the form of electricity.

Attempts to develop and refine OTEC technology started in the 1880s. The 1970s saw an uptick in OTEC research and development during the post 1973 Arab Israeli War, which caused oil prices to triple. The U.S. federal government poured $260 million into OTEC research after President Carter signed a law that committed the US to a production goal of 10,000 MW of electricity from OTEC systems by 1999.

1981 also saw a major development in OTEC technology when Russian engineer, Dr. Alexander Kalina, used a mixture of ammonia and water to produce electricity.

Nowadays, Japan is a major contributor to the development of OTEC technology. Beginning in 1970 the Tokyo Electric Power Company successfully built and deployed a 100-kW closed-cycle OTEC plant on the island of Nauru. The plant became operational on 14 October 1981, producing about 120 kW of electricity; 90 kW was used to power the plant, and the remaining electricity was used to power a school and other places. This set a world record for power output from an OTEC system where the power was sent to a real (as opposed to an experimental) power grid.

In this work, I modelled a vertical pipe of about 1000m in length and 10 to 20 m in diameter that is required in typical ocean thermal energy conversion units. Such a cold-water pipe (CWP) is used to raise the cooler water from the ocean depths to the warmer surface water where the resulting fluid temperature difference of 10 to 20 deg C is sufficient to produce net power through heat exchange. Motion analyses of CWP system under excitation of ocean currents, waves, and the deployment barge to which the upper end is attached is presented in [1]. The analysis, based on the work of Wilson et al, addresses the problem of dynamic stability for a uniform, continuous, vertical CWP which allows for an arbitrary elastic rotational restraint at the barge end. The barge heave and sway motions are included as the end excitation parameters, and Morison’s equation is employed to account for the wave drag and inertial forces along the pipe length. The wave environment is simulated by discretising the Pierson-Moskowitz wave spectrum.

Figure 1. Mathematical model of the OTEC pipeline-barge system, [1]

In this article, I present Python code for dynamic analysis of the OTEC pipe system.

The transverse displacement for the pipeline v(x,t) can be approximated by Bernoulli-Euler equation:

Bernoulli-Euler equation, (1)

where P(x,t) — the longitudinal tension, q(x,t) — the transverse wave loading per unit length, EI — bending stiffness, m - pipe mass per unit length. Complete list of parameters required to solve the problem can be found in [1].

Motion is restricted to the plane in which q is maximum, which is in directions of wave.

The Pierson-Moskowitz wave is presented by:

def S_eta(omega):
    return A*omega**(-5)*np.exp(-B/omega**4)

The sea surface wave amplitude is given by:

def eta(omega_n, del_omega):
    return (2*S_eta(omega_n)*del_omega)**0.5

The heave motion of the barge is :

def H(t, del_omega):
    return np.sum([R(i*del_omega)*eta(i*del_omega,del_omega)*cos(i*del_omega*t+uniform(0, 2*pi)) 
                   for i in range(1,N+1)])

Where R is given by:

def R(omega_n):
    if 0<omega_n<=omega_b :
        return 1
    elif omega_b<omega_n<=omega_N:
        return 1-(omega_n-omega_b)/(omega_N-omega_b)
    else:
        return 0

Sway motion of the barge can be presented as:

def S(t, del_omega):
    return (eta(omega_p, del_omega)*np.sin(e1*omega_p*t+uniform(0, 2*np.pi))+
            np.sum([e2*eta(i*del_omega,del_omega)*np.sin(e1*i*del_omega*t+uniform(0, 2*np.pi))
                    for i in range(int(omega_p/del_omega+1),N+1)]))

The longitudinal tension equals:

def P(x,t, del_omega):
    x_=Symbol('x_')
    t_=Symbol('t_')
    def P_(x_,t_, del_omega):
        return M1*g + (L-x_)*m1_hat*g + (L-x_)*m_hat*diff(H(t_, del_omega),t_,2) + M0*diff(H(t_, del_omega),t_,2)
    return lambdify((x_, t_), P_(x_,t_,del_omega), 'numpy')(x,t)

Whereas, the transverse wave loading is :

def q_hat(t, x, del_omega, ddv):
    z = -x 

    def du(z,t,del_omega):
        return np.sum( [(i*del_omega)**2*eta(i*del_omega, del_omega)*np.exp((i*del_omega)**2/g*(z+d))*np.sin(i*del_omega*t+uniform(0, 2*np.pi))/np.sinh((i*del_omega)**2/g*d)
                   for i in range(1,N+1)])

    return 1.1*0.25*np.pi*CM*rho*D**2*du(z,t,del_omega)-0.25*np.pi*CA*rho*D**2*ddv

I have utilised a finite difference method to solve equation (1).

Some time integration parameters were as follows:

T = 600  # Maximum time bound
del_t =  L/np.sqrt(E/rho_steel)
m = 16  # No. of intervals in space
n = int(T/del_t)  # No. of intervals in time
h = L / m  # Step size in space
k = T / n  # Step size in time

With the initial conditions given as :

# Initialization of solution
v = np.zeros((m + 1, n + 1))

# Initial condition
ic1 = lambda t: S(t, omega_N/N)  

# Time discretization
for j in range(1, n + 2):
    v[0, j - 1] = ic1((j - 1) * k)

Finite difference scheme for the equation (1) is :

def ddddv_ddddx(y0,y1,y2,y3,y4,h):
    return (y0-4*y1+6*y2-4*y3+y4)/h**4

def ddv_ddx(y1,y2,y3, h):
    return (y3-2*y2+y1)/h**2

def dv_dx(y1,y3,h):
    return (y3-y1)/(2*h)

def dP_dx(x1,x3,h,t,del_omega):
    return (P(x3,t,del_omega)-P(x1,t,del_omega))/(2*h)

def gov_eq(y0, y1, y2, y3, y4, t3, t5, h, k, x,t):
    d_omega=omega_N/N
    return (EI*ddddv_ddddx(y0,y1,y2,y3,y4,h)
            - dP_dx(x-h,x+h,h,t,d_omega)*dv_dx(y2,y4,h)
            - P(x,t,d_omega)*ddv_ddx(y2,y3,y4, h)
            + m_hat*ddv_ddx(t3,y4,t5, k)
            - q_hat(t, x, d_omega, ddv_ddx(t3,y4,t5, k)))

Solution is obtained by:

y4=Symbol('y4')

for j in tqdm(range(1, n)):
    for i in range(-1, m-1):
        v[i+2][j] = solve(gov_eq(v[i - 2, j], 
                                 v[i - 1, j], 
                                 v[i, j], 
                                 v[i + 1, j],
                                 y4, 
                                 v[i , j - 1], 
                                 v[i , j + 1], 
                                 h, 
                                 k, 
                                 (i+2)*h,
                                 j*k), y4)[0]

Results of the modelling for this particular setting produced a bounded response of the pipe and were validated against [1].

Figure 2. Results of the modelling

Codebase for the model can be found here: https://github.com/gigatskhondia/ocean_intella

To keep up with the project please visit https://gigala.io/

[1] Dynamics of offshore structures, J. F. Wilson (Chapter 10)