# This file is part of OpenDrift.
#
# OpenDrift is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 2
#
# OpenDrift is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with OpenDrift. If not, see <https://www.gnu.org/licenses/>.
#
# Copyright 2015, Knut-Frode Dagestad, MET Norway
# This ShipDrift model is based on
# Soergaard, E and T Vada (1998);
# Observations and modelling of drifting ships. Report DnV 96-2011
# Reprogrammed in Python for OpenDrift by Knut-Frode Dagestad Dec 2016
import os
import numpy as np
import scipy
import logging; logger = logging.getLogger(__name__)
from opendrift.models.basemodel import OpenDriftSimulation
from opendrift.elements import LagrangianArray
from opendrift.config import CONFIG_LEVEL_ESSENTIAL, CONFIG_LEVEL_BASIC, CONFIG_LEVEL_ADVANCED
[docs]
class ShipObject(LagrangianArray):
"""Extending LagrangianArray with variables relevant for leeway objects.
"""
variables = LagrangianArray.add_variables([
('orientation', {'dtype': np.uint8,
'units': '1',
'default': 1}),
('length', {'dtype': np.float32,
'units': 'm',
'min': 1,
'max': 500,
'description': 'Length of ship',
'level': CONFIG_LEVEL_ESSENTIAL,
'default': 80}),
('height', {'dtype': np.float32,
'units': 'm',
'min': 1,
'max': 100,
'description': 'Total height of ship (above and below waterline)',
'level': CONFIG_LEVEL_ESSENTIAL,
'default': 8}),
('draft', {'dtype': np.float32, # wet part of ship [m]
'units': 'm',
'min': 1,
'max': 30,
'description': 'Draft of ship (depth below water)',
'level': CONFIG_LEVEL_ESSENTIAL,
'default': 4.0}),
('beam', {'dtype': np.float32, # width of ship
'min': 1,
'max': 70,
'units': 'm',
'description': 'Beam (width) of ship',
'level': CONFIG_LEVEL_ESSENTIAL,
'default': 10}),
('wind_drag_coeff', {'dtype': np.float32, # Cf
'units': '1',
'default': 1}),
('water_drag_coeff', {'dtype': np.float32, # Cd
'units': '1',
'default': 1}),
('jibeProbability', {'dtype': np.float32,
'units': '1/h',
'default': 0.04}),
])
[docs]
class ShipDrift(OpenDriftSimulation):
"""Demonstration trajectory model based on OpenDrift framework.
Simply advects a particle (passive tracer with
no properties except for position) with the ambient wind.
"""
ElementType = ShipObject
required_variables = {
'x_wind': {'fallback': None},
'y_wind': {'fallback': None},
'land_binary_mask': {'fallback': None},
'x_sea_water_velocity': {'fallback': None},
'y_sea_water_velocity': {'fallback': None},
'sea_surface_wave_stokes_drift_x_velocity': {'fallback': 0},
'sea_surface_wave_stokes_drift_y_velocity': {'fallback': 0},
'sea_surface_wave_significant_height': {'fallback': 0},
'sea_surface_wave_mean_period_from_variance_spectral_density_second_frequency_moment': {'fallback': 0}
}
winwav_angle = 20 # Angular offset in degrees
def __init__(self, *args, **kwargs):
# Read ship properties
d = os.path.dirname(os.path.realpath(__file__))
w = open(d + '/wforce.dat', 'r')
self.wforce = {}
w.readline()
nbeam = int(w.readline().split()[0])
self.wforce['nbeam'] = nbeam
self.wforce['BL'] = np.array(w.readline().split()[0:nbeam],
dtype=float)
ndraft = int(w.readline().split()[0])
self.wforce['ndraft'] = ndraft
self.wforce['DL'] = np.array(w.readline().split()[0:ndraft],
dtype=float)
nomega = int(w.readline().split()[0])
self.wforce['nomega'] = nomega
self.wforce['omega'] = np.zeros((nomega))
self.wforce['F'] = np.zeros((nomega, nbeam, ndraft))
self.wforce['D'] = np.zeros((nomega, nbeam, ndraft))
for o in range(nomega):
self.wforce['omega'][o] = float(w.readline().split()[0])
for i in range(ndraft):
l = w.readline()
self.wforce['F'][o, i, :] = l.split()[0:nbeam]
for i in range(ndraft):
l = w.readline()
self.wforce['D'][o, i, :] = l.split()[0:nbeam]
w.close()
wi_omega, wi_BL, wi_DL = \
np.meshgrid(self.wforce['omega'],
self.wforce['BL'], self.wforce['DL'], indexing='ij')
self.wforce_interpolator_F = scipy.interpolate.LinearNDInterpolator(
(wi_omega.ravel(), wi_BL.ravel(), wi_DL.ravel()),
self.wforce['F'].ravel())
self.wforce_interpolator_D = scipy.interpolate.LinearNDInterpolator(
(wi_omega.ravel(), wi_BL.ravel(), wi_DL.ravel()),
self.wforce['D'].ravel())
super(ShipDrift, self).__init__(*args, **kwargs)
self._add_config({'seed:orientation': {'type': 'enum',
'enum':['left', 'right', 'random'], 'default': 'random',
'level': CONFIG_LEVEL_ESSENTIAL,
'description': 'If ships are oriented to the left or right of the downwind direction,'
'or whether this is unknown. Left/right means that wind will hit ship from backboard/steerboard'}})
# Since the ShipDrift model is deterministic for given ship size,
# (in contrast to the Leeway model), we use a default diffusivity
# to yield some variability.
self._set_config_default('drift:horizontal_diffusivity', 100)
self._set_config_default('drift:max_speed', 2)
[docs]
def seed_elements(self, *args, **kwargs):
if 'number' in kwargs:
num = kwargs['number']
else:
num = self.get_config('seed:number')
for var in ['length', 'height', 'draft', 'beam']:
if var not in kwargs:
kwargs[var] = self.get_config('seed:' + var)
kwargs[var] = np.atleast_1d(kwargs[var])
if len(kwargs[var] == 1):
kwargs[var] = kwargs[var]*np.ones(num)
# Check that beam and height vs length are within expected range
dl = kwargs['draft']/kwargs['length']
if dl.min() < 0.025 or dl.max() > 0.07:
logger.warning('Ratio of draft to length should be in range '
'0.025 to 0.07, given range is %s-%s. '
'Using border value.' %
(dl.min(), dl.max()))
dl = np.clip(dl, 0.025, 0.07)
bl = kwargs['beam']/kwargs['length']
if bl.min() < 0.12 or bl.max() > 0.18:
logger.warning('Ratio of beam to length should be in range '
'0.12 to 0.18, given range is %s-%s. '
'Using border value.' %
(bl.min(), bl.max()))
# Calculate drag coefficients based on given ship dimensions
# Wind drag coefficient
exposed = kwargs['height'] - kwargs['draft']
Cf = np.zeros(num)
Cf[exposed>37.2] = 1.4
Cf[exposed<=37.2] = 1.045 + 0.016*(exposed[exposed<=37.2] - 15.)
Cf[exposed<=15] = 0.700 + 0.023*exposed[exposed<=15]
kwargs['wind_drag_coeff'] = Cf
# Water drag coefficient
beta = 2.0*dl
Cd = np.zeros(num)
Cd[beta>.12] = 1.27
Cd[beta<=.12] = 1.32 + (1.27-1.32)/0.02 * (beta[beta<=.12]-0.10)
Cd[beta<=.10] = 1.38 + (1.32-1.38)/0.02 * (beta[beta<=.10]-0.08)
Cd[beta<=.08] = 1.44 + (1.38-1.44)/0.02 * (beta[beta<=.08]-0.06)
Cd[beta<=.06] = 1.50 + (1.44-1.50)/0.01 * (beta[beta<=.06]-0.05)
kwargs['water_drag_coeff'] = Cd
if 'orientation' not in kwargs:
oc = self.get_config('seed:orientation')
if oc == 'left':
kwargs['orientation'] = np.ones(num)*0
elif oc == 'right':
kwargs['orientation'] = np.ones(num)*1
else:
kwargs['orientation'] = np.r_[:num] % 2 # Random 0 or 1
# Calling general constructor with calculated values
super(ShipDrift, self).seed_elements(*args, **kwargs)
[docs]
def update(self):
Tm = self.wave_period()
Hs = self.significant_wave_height()
dl = self.elements.draft/self.elements.length
bl = self.elements.beam/self.elements.length
# Clip to allowed range
bl = np.clip(bl, 0.12, 0.18)
dl = np.clip(dl, 0.025, 0.07)
# Additional clipping to avoid NaN from interpolator
bl = np.clip(bl, 0.121, 0.179)
dl = np.clip(dl, 0.0251, 0.069)
# Simply move particles with ambient current
self.update_positions(self.environment.x_sea_water_velocity,
self.environment.y_sea_water_velocity)
# Wind force
rho_air = 1.25 # to be checked
area_dry = self.elements.length*(self.elements.height -
self.elements.draft)
area_wet = self.elements.length*self.elements.draft
# Wind force
F_wind = (0.5*rho_air*self.elements.wind_drag_coeff*
area_dry*np.power(self.wind_speed(), 2))
# Decompose wind
F_wind_x = F_wind*self.environment.x_wind/self.wind_speed()
F_wind_y = F_wind*self.environment.y_wind/self.wind_speed()
F_wind_x[self.wind_speed()==0] = 0
F_wind_y[self.wind_speed()==0] = 0
# Wave force
rho_water = 1025
NSPEC = 100
ommin2 = 2.25
ommin3 = 7.0
ommax = 12.0
dom = (ommax-ommin2)/(NSPEC-1)
F_wave = np.zeros(self.num_elements_active()) # For each ship
beta2 = 0
scale1 = np.sqrt(9.81/self.elements.length)
tmp = np.power(2.0*np.pi/Tm, 4)
d = tmp*Hs*Hs/(4*np.pi)
b = tmp/np.pi
s = np.zeros((NSPEC, self.num_elements_active()))
for i in range(NSPEC):
omi = ommin2 + i*dom
omi = omi*scale1
s[i, :] = d * np.exp(-b/np.power(omi, 4)) / np.power(omi, 5) # m2s
f2 = 0.0
d2 = 0.0
for i in range(NSPEC):
omi = ommin2 + i*dom
f1 = f2
d1 = d2
#print omi - ommin3, 'should be positive'
#print omi
if (omi < ommin3):
f2 = self.wforce_interpolator_F(omi, bl, dl)
d2 = self.wforce_interpolator_D(omi, bl, dl)
else:
# Interval 3
f2 = 0.5
d2 = 4.0*omi*f2
F_wave = F_wave + 0.5*(f1+f2)*dom*scale1*np.power(s[i,:], 2)
beta2 = beta2 + 0.5*(d1+d2)*dom*scale1*np.power(s[i,:], 2)
F_wave = F_wave*rho_water*9.81*self.elements.length
beta2 = beta2*rho_water*np.sqrt(9.81*self.elements.length)
# Add calculated wave and wind drift
longperiod = Tm > 8.55
F_wave_b = F_wave.copy()
F_wave[longperiod] = F_wave[longperiod]*.66
beta2[longperiod] = beta2[longperiod]*.60
medperiod = ((Tm >= 5.7) & (Tm <= 8.55))
F_wave[medperiod] = F_wave[medperiod]*(1.0 - 0.34*(Tm[medperiod]-5.7)/2.85)
beta2[medperiod] = beta2[medperiod]*(1.0 - 0.4*(Tm[medperiod]-5.7)/2.85)
# Form drag (water resistance)
beta1 = 0.5*rho_water*self.elements.water_drag_coeff*area_wet
# Wave direction is taken as wind direction plus offset +/- 20 degrees
# Using minus, since angles are here defined counter-clockwise from x-axis
offset = -self.winwav_angle*2*(self.elements.orientation - 0.5)
if (self.environment.sea_surface_wave_stokes_drift_x_velocity.max() == 0 and
self.environment.sea_surface_wave_stokes_drift_y_velocity.max() == 0):
logger.info('Using wind direction as wave direction')
wave_dir = np.radians(offset) + np.arctan2(self.environment.y_wind,
self.environment.x_wind)
else:
logger.info('Using Stokes drift direction as wave direction')
wave_dir = np.radians(offset) + np.arctan2(
self.environment.sea_surface_wave_stokes_drift_y_velocity,
self.environment.sea_surface_wave_stokes_drift_x_velocity)
F_wave_x = F_wave*np.cos(wave_dir)
F_wave_y = F_wave*np.sin(wave_dir)
F_total = np.sqrt(np.power(F_wind_x + F_wave_x, 2) +
np.power(F_wind_y + F_wave_y, 2))
# Iterate 4 times in order to estimate the effect of
# wave damping and form drag
uw_tot = 0
uw_dir = 0
for i in range(4):
# Calc wave damping force in x and y (Ref (1), eq. 6.6.4)
f2x = beta2*uw_tot*np.cos(wave_dir)
f2y = beta2*uw_tot*np.sin(wave_dir)
# Calc new drift direction, including effect of
# wave damping (Ref (1), eq. 6.6.1)
uw_dir = np.arctan2(F_wind_y+F_wave_y-f2y, F_wind_x+F_wave_x-f2x)
# Ref (1), eq. 6.6.3
bet2c = beta2*np.cos((wave_dir-uw_dir))
# Ref (1), eq. 6.4.3
uw_tot = -bet2c/(2.*beta1) + np.sqrt(bet2c*bet2c +
4.*beta1*F_total)/(2.*beta1)
# Finally advect according to wind-wave forces
velocity_u = uw_tot*np.cos(uw_dir)
velocity_v = uw_tot*np.sin(uw_dir)
self.update_positions(velocity_u, velocity_v)
# Stranding
self.deactivate_elements(self.environment.land_binary_mask == 1,
reason='ship stranded')