A Few Frontier Issues in Ocean Engineering Mechanics

In this review, we primarily address the present state of the arts and latest progresses in a few frontier issues mostly relevant to free surface/interface in ocean engineering. They include TC (tropical cyclone) induced extreme surface wave, sloshing of LNG (liquefied natural gas), cavitation/bubble dynamics and VIM (vortex-induced motion) and VIV (vortex-induced vibration). In addition to general description, we mainly focus on the recent advances and challenging aspects of above-mentioned topics. Inspired by the achievements in the previous 70 years, mankind starts a new round of ocean exploration activities. Then, we can find obvious trends: the realm of ocean engineering is expanding from sea surface to deep sea, from low and middle latitude to polar region and from fossil to renewable energy in near future.


Introduction
The azure ocean is alluring, the roaring ocean is awe-inspiring, while the Arctic Ocean is daunting. Claiming marine resource and space has never been an easy task. Although regional maritime commerce activities on the Aegean Sea became active during Ancient Greece Age, it was not until 15th century for Zheng He's fleet to reach Southeast Asia and East Africa, and for Magellan to accomplish circumnavigation. The modern marine science only germinated at the turn of 19−20th centuries, and then the physical oceanography grew up in 1950s. Ocean engineering featured by exploration and production of marine petroleum and gas can be dated back to the launching of the first drilling platform in the Mexico Gulf out of sight of land in 1947. China began its marine exploration over the Yinggehai in1956; however, the first fixed platform in the Bohai Sea was erected in 1966.
With the great efforts of more than fifty years, platforms have already towered aloft in great numbers in the Gulf of Mexico, North Sea and West Africa. Moreover, the marching to the deep sea is quite inspiring (see Fig. 1). By the end of 2014, the records of water depth were 2934 m, 2500 m and 3174 m for subsea tree, exploration ship and production equipment, respectively. In China, platforms have expanded from shallow water to the Beibu Gulf and Dongsha islands. Deep-sea Semi drilling platforms have been designed and manufactured to reach the regions of South China Sea over the depth of 3000 m. Proud achievements such as manned deep-water exploration, stable extraction of combustible ice, polar navigation and scientific investigation in this regard are inspiring.
During previous decades, ocean engineering mechanics has made marvelous advances in many aspects. Marine environment parameters of wind, wave and current, external inertial and viscous loads on the structure and corresponding dynamic response of platforms such as springing, ringing or slow drift can be calculated and measured. As a result, safe and efficient operation of drilling and exploitation can be guaranteed. On the other hand, there are still a few long-term hard nuts and new areas waiting for exploration. Hence, we would address a few frontiers mostly concerned with free surface/interface issues such as TC induced extreme surface wave, sloshing of LNG, cavitation/ bubble dynamics, current induced VIV and VIM in this review.
Based on the comprehensive investigation on the recent exploitation in this field, we have seen the trends from near shore to far sea, from shallow water to deep sea, from low/middle to high latitude region, from fossil to renewable energy. Then, deep sea engineering, arctic ocean engineering and hydrate exploitation become most urgent tasks of ocean engineers.

Extreme TC wind and sea surface wave
Tropical cyclone (TC, known as typhoon over Pacific Ocean and hurricane over Atlantic Ocean) turns out an enormous atmospheric vortex of a few hundred kilometers, which consists of a lower pressure eye in the center and rainbelt around. Usually, they are generated in the sea area with sea surface temperature (SST) higher than 26°C beyond latitude of 5° where intense unstable convection and heat transfer as high as 250 W/m 2 occur. As soon as TC is formed, this air mass revolving counterclockwise/clockwise on the northern/southern hemisphere would translate along a certain trajectory. Following steering flow in atmospheric circulation in westerlies, a normal trajectory always moves west/northwest at the average speed of 18−21 km/h (Kossin, 2018). Usually, typhoon is specified as tropical depression, tropical storm, strong tropical storm, typhoon, strong typhoon and superstrong typhoon corresponding to its wind speed larger than 10.8, 17.2, 24.5, 32.7, 41.5 and 51.0 m/s, respectively. In contrast, hurricane is divided into five classes of Saffir-Simpson scale for wind speed larger than 119, 154, 178, 210 and 249 km/h, respectively. Hence, strong TC indeed brings severe threat to ocean engineering structures due to the giant wind speed and huge ocean wave induced.

Extreme TC wind prediction
Routine operational forecast of TC has become fairly mature with ever-growing accuracy (Cangialosi and Franklin, 2016) relying on 4-dimension assimilation technique by integrating numerical modeling on supercomputer and monitoring by meteorological satellites and ground-based radars available. However, it is still very hard to exactly predict ab-rupt variation in TC's strength and path. In particular, the number and maximum wind speed of strong TC are evergrowing though, we have accumulated little knowledge about extreme wind speed of a few decades of return period under the background of climate change. This circumstance would no doubt exert considerable effects on the accurate design in ocean engineering and bring about inevitable losses. For instance, the very maximum wind speed of Hurricane Katrina is that of 50-year-return period rather than 200 years, namely, lower wind speed was used in the structure design. It is not strange that 116 fixed platforms were destroyed directly and 150−200 ones removed in the following years in the Gulf of Mexico.
The simulation of climate system can be implemented to foresee future TC activities. Recently, more aspects such as atmosphere, land surface, ocean, sea ice and biogeochemical cycles have been integrated into the latest climate models. In 2016 around 100 distinct climate models are selected for CMIP6 (Coupled Model Inter-comparison Project Phase 6) to eliminate errors in individual model (Eyring et al., 2016). In GCMs (Global Climate Models), climate variables such as air flow, temperature, moisture and precipitation are solved in each three-dimensional grid of cells with the parameterization of source terms and interface fluxes for mesoscale physical, chemical and biological processes. Validated by the hind-casting tests, GCMs are driven by future forcing scenarios. Say, four RCP (Representative Concentration Pathway) scenarios are investigated in the fifth AR of IPCC focusing on the future warming, in which RCP6.0 and RCP8.5 mean little or no mitigation actions are taken, and RCP2.6 of aggressive mitigation generally limits carbon emission resulting in global warming no more than 2°C. To provide more details for a particular location, RCMs (Regional Climate Models) are embedded. Moreover, dynamic downscaling approach was implemented to fill the gap between results of GCMs and RCMs (Tapiador et al., 2020).
On the other hand, purely historical data-based method with sufficient samples in TC prone areas is initiated to avoid using time-consuming climate modeling. Specifically, the extreme wind speed can be evaluated by statistical analysis of TC wind speed data observed in history regarding TC as random events. Considering the effects of climate change, nonstationary stochastic process is assumed when estimating the distribution functions of their annual occurrence and maximal wind speed. The nonstationary extreme wind speed can be expressed as: where V PR is the extreme wind speed of return period of P R year. (t) is the annual average number of TC, (t), (t) and (t) are the shape, location and scale parameters in Weibull distribution of wind speed, respectively. All of them are re- garded as statistical parameters varying with time. Relied on the nonstationary extreme theory, the extreme wind speeds with 20−50 years return periods are 4.1%-4.4% higher than that predicted by stationary ones in the South China Sea (Wang and Li, 2016).
Many parametric TC models are proposed to reconstruct the pressure and wind fields ten meters above the sea surface usually. Those models present the relationship between the distribution of TC wind/pressure and the main TC parameters such as the central and ambient pressures, maximal wind speed, translation velocity, radius of maximum wind speed and a few coefficients. Holland model is a popular one, in which the pressure p and wind U at distance r from the center can be expressed as (Holland, 1980): ρ a where, p c and p n are the central and ambient pressures; , the air density; R, the radius of maximal wind speed. B is a coefficient reflecting the shape of the wind field vortex. Recently, variation in the wind equation exponent is further introduced into the original Holland model to make it more stable and less sensitive to the input parameters (Holland et al., 2010).

TC induced sea surface wave models
The features of TC wave field are unique and different from conventional wind waves. Firstly, the ranges of parameters of TC wave are much wider. For instance, the significant wave height can be as high as tens of meters in open sea. And, it is hard to reach a fully-developed state since the TC wind field is moving. Moreover, the distribution of wave field is complicated in this kind of moving asymmetric wind field, in which translational speed is one of the most important influential factors.
As shown in Fig. 2, waves generated in the region where wind vectors aligned with the translation direction tend to move forward with TC resulting in the extended fetch effect, i.e., waves can remain in the high wind regions for ε ν longer periods to receive energy from wind. If the translational speed closes to the group velocity of wave, the trapped fetch would result (Bowyer and MacAfee, 2005). Conversely, waves generated on the other side will propagate out of the intense wind areas rapidly. That implies the degree of asymmetry of the wave field would be more than the wind field. Meanwhile, swell with larger propagating speed will be radiated out from the intense wind region. Based on the in-situ data and statistical analysis, the spectra of TC waves are typical unimodal, which is consistent with the relation for fetch-limited growth as expressed in Eq. (3) (Young, 2017), where and are non-dimensional total wave energy and spectral peak frequency, respectively. ε = 6.356 × 10 −6 ν −3.3 .
(3) A half empirical and half theoretical method for predicting wave elements by wind elements was proposed during World War II to meet the need of D-Day invasion. Based on the concept of wave spectrum transmission equation, the first wave model was established in 1960s with defects in overestimated wind energy input and underestimated energy transmission. In spite of numerous modifications, the second generation of wave models developed in 1970s still needs a given spectral shape in advance to simplify nonlinear interaction computation and is unable to apply for modeling waves generated in the rapidly varying wind field such as storms, cyclones or fronts. To avoid these shortcomings, wide international cooperation under SWAMP (Sea Wave Modeling Project) gave birth to the third generation (3G) wave model WAM -the first global forecast model in 1980s. On this basis, the WAVEWATCH III model by NOAA and the SWAN model by Delft university of technology in the Netherlands appeared in succession with optional governing equation, numerical scheme, parametrization schemes provided.

Parametrization scheme of source terms σ θ
Taking the SWAN model for the shallow water as an example, it solves the energy balance equation as Eq. (4), since the wave action density N( , ) is conserved during propagation in the presence of ambient current.

∂N ∂t
where, and are the longitude and latitude components of group velocity of water wave, respectively. U and V are the ambient currents.
( ) is propagation velocity of refraction (frequency shift). The energy input, dissipation and transfer due to wave-wave interactions are considered by the source term S tot as Eq. (5) (Zijema, 2010): (5) Specifically, S in is the wave energy inputted by the wind; S ds,w , S ds,b and S ds,br are wave dissipation terms caused LI Jia-chun, NIE Bing-chuan China Ocean Eng., 2021, Vol. 35, No. 1, P. 1-11 3 by white capping, bottom friction and depth-induced wave breaking, respectively; S nl,3 and S nl,4 represent the threewave and four-wave interactions, respectively. Parametrizations of the source terms, especially the energy input and dissipation under extremely high TC wind speed are far from completeness. Recently, a few progresses based on in-situ observation are reported. For instance, the drag coefficient for estimating wind stress is newly parameterized at wind speed 25−40 m/s based on the in-situ measurements taken beneath five TCs (Hsu et al., 2019). In addition, the effects of wave age, namely, wind fetch and duration should also be taken into consideration. People also find that wave breaking is so significant in both wind energy input and wave energy dissipation that further in-depth theoretical and experimental efforts are extremely needed. To consider the interactions among wave, current and surge in the near coast, the wave model can be further coupled with the circulation models (Wuxi et al., 2018). On the other hand, soft computing methods such as ANN (artificial neural network) and SVM (support vector machine) are adopted for predicting significant wave height in recent years (Mafi and Amirinia, 2017).

Sloshing of internal liquid in partially-filled containers
Sloshing of LNG has become one of the most challenging issues for internal flow loads at present, where LNG indicates a natural gas liquefied at −163°C with the reduction of volume by a factor of 600 for easy transportation and storage. Because of the advantages of clean, high quality, safe and easy use in existing burning devices, the demand of LNG is remarkably increasing in recent years. In fact, sloshing phenomena of different scales from small to very large volume can often be observed for tanks in spacecraft, road vehicles and ships with partially-filled oil inside. However, LNG is far more difficult to handle than water and oil due to the extreme low boiling temperature at the risk in global motions and structure safety. Moreover, the control system in ocean engineering is affected by the sloshing frequency and total forces as well. The specific challenges in this regard include: violent or even chaotic sloshing could take place in increasing size of LNG tanks; extreme slamming pressure would cause severe local damage; hydroelasticity effect matters for the container system under this circumstance; coupling between sloshing and six-DOF motions of the ship subjected to external loads is of significance. In addition, the thermodynamic and compressibility effects due to the air entrapment may also become nonnegligible.

Linear and nonlinear sloshing
Infinite sloshing modes exist for liquid inside a partly filled container. Based on the inviscid and incompressible assumptions, the fluid flow satisfies the Laplace equation with specified boundary conditions. Thus, analytical expres-sions of the corresponding natural frequencies can be provided for typical tank geometries. Taking the 2D sloshing in a rectangular tank as an example, the lowest natural frequency f 1 can be estimated by where, l and h are the tank length and water depth, respectively, as ketches in Fig. 3. When the excitation amplitude is relatively small and its frequency is remote from the natural frequencies, the linear models are applied. The mass-spring-dashpot (Bauer, 1984) and pendulum (Sayar and Baumgarten, 1981) models are probably the most popular ones, in which the part of liquid moving in unison with the tank is regarded as rigid mass, while the sloshing part is modeled as a series of spring-mass systems or pendulums. Relied on the linear equivalent models, the total force and moment of sloshing liquid exerting on the tank can be estimated.
Resonance occurs as the excitation frequency is close to the natural frequencies. Nonlinear effects arise owing to finite free surface motion. Usually, two types of phenomena can be observed during resonance: dramatic increase of the sloshing amplitude and the switch between different sloshing modes (for instance, from the planar sloshing to rotary mode). Those phenomena are closely related to the coupling and instabilities of different sloshing modes. In general, the resonance of sloshing modes with the lowest natural frequencies are of primary interests, in which rotary sloshing is concerned most due to its significant horizontal force or overturning moment.
Besides the excitation frequency, the liquid depth can also affect the sloshing motion by nonlinear steepening and viscous/wave breaking dissipation. Taking sloshing in a rectangular tank as an example (see Fig. 3), traveling hydraulic jumps may occur in shallow liquid, standing wave with largest elevations at the side walls could be observed for finite depth. They can be predicted by the shallow-fluid theory and modal theory for finite depth, respectively. As for the intermediate case, asymptotic modal approximation matching aforementioned two theories are proposed (Faltinsen and Timokha, 2002). Breaking waves are more vis-ible for shallow and intermediate depths. Increasing the water depth further, the influence of the water depth becomes insignificant; however, roof impact may occur. Moreover, there exists a critical depth (h/l = 0.3368). The sloshing response has the ''hard-spring'' behavior (response amplitude grows with increasing frequency) when h/l < 0.3368 and ''soft-spring'' behavior (response amplitude reduces with increasing frequency) when h/l > 0.3368 (Waterhouse, 1994). Near the critical depth, high order sloshing modes may become dominant and intense amplification of the liquid behavior may be observed.
A few nonlinear models and analytical techniques have been developed for predicting near resonance behavior. Typical equivalent models are the spherical pendulum for the rotary sloshing and the compound pendulum (Kana, 1989) for low-liquid-filled tanks in which liquid sloshing has the properties of both rotary and normal linear sloshing. For stronger nonlinearity, more sloshing modes could be involved so as to experience complex evolution of free surface. The energy would be nonlinearly transferred from lower to higher modes, and strong nonlinear amplification of higher modes may be observed (Faltinsen et al., 2005). Nonlinear multimodal method can be used to study wave regimes, multi-branched solutions and hydrodynamic stability properly (Faltinsen, 2017). For violent sloshing under certain specific containers, the bifurcation diagrams in terms of excitation amplitude versus frequency are provided (Saeki et al., 2001). Based on that, sloshing types such as irregular beating, asymmetric, quasi-periodic and chaotic scenarios can be identified.

Effects of sloshing on ship's global motion
The hydrostatic pressure of sloshing liquid can exert so dramatic forces and moments on the ship hull as to affect its global motion. At the same time, the liquid itself is subjected to the excitation of ship motion in turn. Although the analysis of ship motion is quite classic, the consideration of coupling effect between ship motion and sloshing is far from mature. (1) When the sloshing is assumed to be small, linear ship motion is thought to be adequate (Kim et al., 2007). Hence, frequency-domain approaches can be used.
(2) When the nonlinearity of sloshing motion becomes dominant, hybrid approaches combing linear model of ship hull motion and nonlinear model of sloshing motion are preferred. More specifically, the sloshing motion is obtained by solving the Navier−Stokes equations; the sloshing induced forces and moments are incorporated into the terms of wave excitation forces and moments; the six-DOF motions of ship hull can be obtained via the time-domain approaches such as Rankine source panel methods and IRF (Impulse-response function). (3) When the sloshing motion become more violent, the global motion responses may not be assumed linear any more, especially for roll motion due to its smaller inertial moment. A framework numerically dealing with the coupled nonlinear ship motion and nonlinear sloshing has preliminarily been established. And then a few attempts with a certain degree of approximation have been implemented (Bulian and Cercos-Pita, 2018). However, a systematic discussion on the coupling effects is still lacking to date.

Sloshing induced local impact force
During a slamming event, impact between fluid and wall may occur violently. Its extreme high pressure in short interval could damage the container system of LNG tanks and other internal structures such as pump towers. Since the peak pressure of slamming will be filtered and mitigated by the container system prior to the moment when it reaches the hull, its influence on the carrier's global motion is conventionally estimated by amplifying the quasi-static pressure by a safety factor. The loading process of slamming can be described in detail by the concept of ELPs (Elementary Loading Processes) consisting of ELP1, ELP2 and ELP3 (Lafeber et al., 2012). ELP1 represents the direct impact, in which hemispherical pressure waves centered on the contacting point between fluid and wall will be generated both in the liquid and solid. The run-up/run-down of jet liquid along the wall is referred as ELP2. On the other hand, ELP3 stands for the compression process of the entrapped or escaping gas (gas pocket) like the role of a spring. Generally, the peak pressure of ELP1 and ELP2 is comparable, and they are much larger than ELP3. Nevertheless, compressibility effect should be considered in ELP1 or ELP3. A few analytical models have been proposed to evaluate the slamming load, such as the Rankine-Hugoniot model (Eq. (7)), Wagner model (Eq. (8)) and Bagnold model (Eq. (9)) for ELP1, ELP2 and ELP3, respectively.
where p max is the pressure peak; V, the normal velocity during impact; c, the sound speed; , the liquid density; Ma 0 , the Mach number; C p , the parameter related to the angle between the free and wet surfaces; p i , the initial pressure in the gas; , the adiabatic exponent; m e , the mass of impact of liquid; and , the initial volume of the entrapped gas. However, precise prediction of the impact load is still a challenge at present, especially for ELP1 and ELP 3 since they are dependent on the processes such as phase transition, bubble production and free-surface instability (Dias and Ghidaglia, 2018).

Cavitation dynamics and bubble oscillation
In ocean engineering, the phenomena of cavitation and gas bubble widely occur. Cavitation may take place when LI Jia-chun, NIE Bing-chuan China Ocean Eng., 2021, Vol. 35, No. 1, P. 1-11 the local cavitation number drops to the critical threshold, which usually appears around fast-running propellers, high speed vehicles and water entry interface. In contrast, gas bubbles with external gas supply can be observed during underwater explosions or wave-breaking. Both cavitation and gas bubbles can change external load, emit pressure pulse or produce noise during collapse, and thus reduce the efficiency or even threat the safety of components. On the other hand, the bubbly flow exhibiting lower viscosity and density can find various applications in engineering such as drag reduction and mixing enhancement. The exploration of cavitation and gas bubble evolution belongs to the realm of interface dynamics.

Oscillation and violent collapse of interface
The interface of cavitation and gas bubbles tends to be a spherical shape owing to the minimization of surface tension energy. They may display both forced and free oscillations. Various spherical interface oscillation models integrating the effects of compressibility, surface tension, viscosity and even mass/heat transfer have been developed. Most well-known ones probably are the Rayleigh model (Rayleigh, 1917), Gilmore model (Gilmore, 1952) and Keller− Miksis model (Keller and Miksis, 1980). Generally, these spherical interface models are presented in an ODE (ordinary differential equation) of the radius R(t), and the oscillation behavior can be solved once specific external force and initial radius/speed of bubble are given. However, it is noteworthy that the nonlinearity of the bubble model can trigger complicated dynamic behaviors corresponding to a fixed point, limit cycle or strange attractor in the state space (Parlitz et al., 1990).
Non-spherical instability might occur when spherical interface contracts strongly during oscillation. The classical Rayleigh−Taylor instability is responsible for this phenomenon. And interface modes and corresponding critical frequencies can be estimated by linear instability theory via expanding the non-spherical bubble interface into a series of spherical harmonics.
During violent collapse, the seriously distorted nonspherical interface could further develop into a toroidal shape resulting in the failure of the linear instability analysis. Fortunately, the law that the variation rate of Kelvin impulse equals external force holds (Blake, 1988). Since the interface usually collapses in a period so short that the variation of external force integral over time is insignificant, the Kelvin impulse can be regarded as a constant. On the other hand, the Kelvin impulse is only related to the global migration and collapse. Thus, the Lagrangian equation of the bubble-fluid system, Eq. (10) can be derived (Benjamin and Elli, 1966).
where M is the added mass, and x 0 is the position of bubble centroid. J, the Kelvin impulse, is dependent on the deform-ation relative to the centroid. J increases with the degree of the skewness of interface and velocity distribution, and vanishes for the spherical interface. To compensate the sharp drop of M caused by volume reduction during collapse, the translation motion must speed up. However, it is still not enough to maintain the Kelvin impulse J due to the limited translation speed. Hence, the only possible way is to transform single connected liquid into multi-connected one, which is responsible for the formation of liquid jet and toroidal interface.

Bubble−bubble and bubble−boundary interactions
When bubbles present in droves, direct and indirect bubble−bubble interactions may occur. They are mainly governed by the normalized distance (i.e. d bb defined as d/(R 1max + R 2max )), oscillating phase difference and bubble size. When d bb ≥ 10, the bubble−bubble interaction is weak, and a bubble may migrate or be trapped in the pressure field emitted by the other bubbles; when d bb < 10, behaviors such as catapult, coalescence, jet towards each other or jet away can be observed when the normalized distance and oscillating phase difference vary (Fong et al., 2009).
If the cavitation or bubble collapse nearby a boundary such as free and solid surfaces, the interface deformation will be remarkably influenced, as shown in Fig. 4. Based on the Kelvin impulse theory of semi-infinite domain, the migration and jet directions are determined by the product of stand-off parameter and buoyancy parameter , in which and are defined as:

γδ γδ
where d is the distance from the bubble center to the boundary, and R m is the maximal radius. For bubbles above a rigid surface (gravity downwards), the bubble will be repelled-away with a jet directed away when >0.442, and attracted with a jet towards the surface when < 0.442. As γδ for a bubble below the free surface, the neutral curve dividing the directions of bubble migration and jet is = 0.442 as well but in opposite directions (Blake et al., 2015). The maximum jet volume normalized by maximal volume during rebound is proportional to the nondimensional Kelvin impulse (Supponen et al., 2016).
The free and solid surfaces will be affected by the nearby oscillation bubble in turn. Taking bubble bursting under the free surface as an example: (1) a concave free surface with capillary waves is induced when the bubble suddenly opens the water layer above, (2) the capillary wave fronts collapsing on the axis bring the singularity of curvature and (3) reversal of the curvature due to singularity triggers the jet (Gañán-Calvo, 2017). In other words, the singularity of free surface and jet is driven by the surface tension and inertia force, and stopped by the viscosity (Zeff et al., 2000).

Cavitation evolution and adverse effects
Hydrodynamics of cavitation occurring at the suction side of surfaces is the most concerned one in practical experiments and applications. Relied on platforms such as cavitation tunnel, new techniques of measurements such as X-ray (Ganesh et al., 2016) along with numerical simulations, the appearance, evolution and impacts of cavitation can be revealed. In general, incipient cavitation, sheet cavitation, cloud cavitation and supercavitation present in succession with the decreasing of cavitation number. Incipient cavitation is the beginning of cavitation with nucleation boom and a few small individual bubbles. When the population and size of bubbles increase, they could merge together to form sheet cavitation. Whereas, cloud cavitation is always characterized of unclear vapor-liquid interface. However, a big stable bubble encompassing the object can be observed for supercavitation. It is noteworthy that the transition between them is also dependent on the Reynold number in addition to cavitation number.
Cavitation can experience periodic formation, shedding and collapse. They will cause the low frequency pressure pulsation in the flow, which is highly related to adverse effects such as structure vibration, noise radiation and structure damage. At present, two mechanisms have been identified for the quasi-periodic shedding. One is the re-entrant jet mechanism (Wang et al., 2012). It works together with the main flow like a pair of scissors cutting the sheet cavitation. The other one is the bubbly shock propagation induced shedding (Ganesh et al., 2016). Namely, the collapse of cavitation cloud downstream can generate a shock upstream to induce the detachment of the next cloud cavitation. Furthermore, collapse of cloud cavitation can produce jet and pressure wave to cause the material damage and noise. In this process, pressure waves with peak as high as tens MPa could be detected. And the nearby material might experience cavitation erosion consisting of plastic deformation hardening stage (known as pits) and mass loss stage including acceleration (or accumulation), deceleration (or attenuation) and steady-state (or terminal) periods.

Structural damage by a large underwater explosion bubble
In an underwater explosion, about half of the explosion energy is converted to a spherical shock wave propagating away. When the shock strikes on a deformable structural surface or free surface, cavitation and subsequent reloading effect may occur. On the other hand, the chemical reaction will turn the charge into explosion product to form explosion bubble of scale of meters. This hot explosion bubble will expand due to the high pressure inside, and over-expand because of inertia, then contract and collapse. Although, this expanding-contracting cycle may repeat several times, majority of the total energy is consumed in the first period. The time scales of the shock and pulsating bubble differ a lot, the former one is of order of millisecond, whereas the bubble can last seconds. In general, these two processes are investigated separately.
Probably, the most important parameters in an underwater explosion are the peak pressure of the shock, the maximal radius and oscillating period of the explosion bubble. They can be obtained using the scaling laws Eqs. (12)−(14) by Cole (1965) with surprisingly good precision.
where, P m is pressure peak in megapascal, a m and T b are maximal radius (m) and oscillating period (s) of underwater explosion bubble in free field, respectively. C is the charge weight (kg), S, the explosion distance (m), and H, the water depth (m) of explosion.
It is not surprising that if the bubble frequency is close to the natural frequency of slender structures, whipping response will be triggered. A dominant dimensionless parameter E n expressed as Eq. (15), which incorporating the influence of bubble load, structural resistance ability and fluid-structure interaction can be used to evaluating the severity of global structural response.
ρ where E, I, and A are the nondimensional elastic modulus, inertial moment, material density, and area of the cross-section, respectively. M a , the nondimensional added mass, and R/L, the slenderness ratio. On this basis, the concepts of near-, middle-and far-fields can be defined, that is: the ma-terial starts to yield in the near field, thus leading to structure breakdown immediately; remarkable elastic response can be found in the middle field in addition to substantial movement; in contrast, attentions should be paid to overall structural motion as a rigid body in the far-field. Further, the thresholds dividing the near-, middle-and far-fields can be determined for a given structural material (Nie et al., 2015).

External current induced structural VIV and VIM
Current induced structure response is quite common in ocean engineering for marine vehicles, platform foundations, risers, umbilical tubes and mooring lines in surface/ internal wave induced flow and tidal current. Then, flow separation occurs time and again under the adverse pressure gradient in the boundary layer. Vortices will be shed from the body regularly or irregularly. Transverse lift and longitudinal in-line drag would lead to dramatic structure movement, oscillation or even damage.

Patterns of vortex shedding in VIM and VIV
Vortex-induced motion (VIM) and vortex-induced vibration (VIV) are two kinds of the most concerned responses caused by vortex shedding. For instance, VIM of platforms in current can cause the fatigue of mooring system (Fujarra et al., 2012), while VIV is able to damage risers (Bearman, 1984). As a matter of fact, the fundamental mechanisms of VIM and VIV are essentially the same and share the similar self-excited and self-limited behaviors. However, differences between them are obvious as well. VIM corresponds to the motion of large bluff rigid bodies with smaller structural aspect ratio and larger Reynolds number. Whereas VIV corresponds to the dynamic response of flexible slender structures with smaller Reynolds number. Besides, the oscillation period of VIM more than 10 s is much longer than that of VIV.
The vortex shedding in the wake behind a bluff body under different Reynold numbers have been examined extensively. Perhaps, the most well-known one could be the Kármán vortex street, which takes place widely in both nature and engineering. When a flow is passing through a structure, either free to move or flexible to vibrate, the patterns of vortices are much more complicated due to fluidstructure interaction. The typical patterns include 2S, 2P, 2C and 2T modes (Flemming and Williamson, 2005), as sketched in Fig. 5. Specifically, 2S mode stands for two single vortices shed per cycle, i.e., the classic Kármán vor-tex street mode. However, 2P mode represents two pairs of vortices shed per cycle with opposite sign in each pair, while two co-rotating vortex pairs are formed per cycle for 2C modes. As for the 2T mode, two triplets of vortices are shed per cycle. Besides, hybrid modes of previous basic modes can also be found. For example, when a pivoted cylinder subjected to non-uniform current, 2P mode nearby the tip and 2S mode around the pivot can be observed. Moreover, the transitions between different modes such as 2P to 2S or 2S to 2P may happen as well.

Self-excited/self-limited oscillation and lock-in
When a uniform flow passing by a cylinder, the frequency f v of vortex shedding can be evaluated by where, U is the flow speed and D is the structural diameter.
St, Strouhal number, usually keeps invariant around 0.2 for a large range of Reynolds number. If the cylinder is free to move or flexible to vibrate, structural self-excited oscillation could be triggered when the vortex-shedding frequency f v is close to its natural frequency f s . Usually, resonance occurs in a considerable range of vortex shedding frequency neighboring f s . Then, the frequency of vortex shedding does not follow the Strouhal relationship any longer and would be locked-in to a certain structural natural frequency. During lock-in, the response is self-limited to moderate amplitude in the lock-in regime. In fact, structural response can affect vortex shedding, and vortex shedding can affect the amplitude of structural response in turn via the hydrodynamic force. Specifically, the shedding modes of this fluid-structure system can be optimized by controlling self-excitation and added mass from the very beginning, so that the shedding frequency can be maintained around the natural frequency. However, the self-adaptation fails when the amplitude is large enough. Consequently, the deviation between vortex shedding frequency and natural frequency increases resulting in the amplitude reduction and system detuning. In addition to Reynold number, the most influential factors include reduced speed V r , mass ratio m*, damping ratio and aspect ratio L/D as well. V r is written as U/(f s D). m* is defined as the ratio of structural mass and the fluid mass displaced by the structure. is the ratio of the structural damping to critical damping in fluid. Two primary issues people concerned most are when lock-in occurs and how large is the maximal amplitude during lock-in in the parameter space.
(1) The range of lock-in is mainly dependent on the relation between the vortex shedding frequency f v and actual structural oscillation frequency f o . Generally, the ratio of f v to f s can be used to judge whether lock-in occurs assuming the structure oscillates at its natural frequency (f o ~ f s ). Using the Strouhal relation for f v , f v /f s can be written as StU/(f s D), i.e., lock-in occurs when V r ~ 1/St. In practice, lock-in can be observed when 3 < V r < 10. However, the assumption f o ~ f s may fail when m* is small, because f o is not a constant any more for small m*; instead, it will be dramatically influenced by the added mass. Specifically, f o increases with V r , since the added mass decreases with V r . Consequently, the system can modulate more easily to let f o be close to the vortex shedding frequency f v , which will enlarge the range of lock-in. Say, lock-in of structure of low m* can even be observed when V r is as high as 12. To exclude the influence of m*, the so-called true reduce velocity U/(f o D) is also proposed to judge the occurrence of lock-in. It seems that lock-in presents when 5 < U/(f o D) < 8. (2) The maximal amplitude of response A * max during lock-in can be concisely represented dependent on (m*+C A ) , or even m* . Where C A is the ratio of the added mass to the displaced mass, which is usually around 0.95−1.05. Many combined parameters consisting of mass and damping ratios are proposed to evaluate the maximal response amplitude during lock-in, say, the Griffin curve versus Skop-Griffin parameter. And discussions on the valid condition of combined parameter collapsing peak amplitude data were carried out, i.e., (m*+C A ) >0.006 (Khalak et al., 1999). Incidentally, it is also suggested by many researchers that mass and damping ratios should not be combined since they can affect the system independently (Sarpkaya, 2004). However, the study on the maximal response amplitude for low mass and damping seems to be still unavailable yet thus far. 5.3 Two-/three-dimensional multi-branch oscillation. ζ ζ ζ During lock-in, multi branches of response exist with varying V r . Two distinct types of response exist depending on whether m * is high or low. Only two branches (initial excitation branch and lower branch) with one hysteretic transition between them can be observed for high and moderate m * . The so-called hysteretic transition means the curve of response amplitude is dependent on the direction of approaching the transition regime, i.e., from a lower or a higher reduced speed. The initial excitation branch and lower branch are related to the 2S mode and 2P mode of vortex shedding, respectively. As for low m * , an additional upper branch is reported (Govardhan and Williamson, 2000), and then two transitions followed. The transition between the initial excitation branch and the upper branch is hysteretic corresponding to the jump between 2S and 2P mode, whereas the transition between the upper and lower branches is intermittent corresponding to the jump of phase between total force and displacement. If m* is sufficiently small (less than a critical value 0.54), the lower branch will never be reached. As for the near-wall cylinders such as suspending submarine pipelines, VIV will be further affected by the gap-to-diameter ratio owing to the interaction between the boundary layers of the seabed and cylinders. The overlap range of reduced velocity during hysteresis transition will increases with the decrease of gap-to-diameter ratio (Liu and Gao, 2019).
As a matter of fact, the response could be two-dimensional, or even three-dimensional. In general, the transverse amplitude is larger than in-line amplitude, while transverse oscillation frequency could be nearly half of the in-line frequency. Coupling between the transverse and in-line responses brings about the various trajectories of cross-section displacement such as eight and crescent shapes. This coupling can also influence the peak amplitude, say, it was reported that the transverse amplitudes of the 8-shape trajectory can be 19% higher than that of 1DOF (Sarpkaya, 1995). Moreover, the trajectories along the long flexible cylinder can be entirely different suggesting the occurrence of multi-mode and multi-frequency VIV.
To predict VIV, force-based models and flow-based models are used. Oscillating lift and drag are derived by integration with coefficients based on empirical data in the force-based model. While, flow-based models solve the flow around to derive the forces on the structure. Because of the similarity between the self-limitation and self-excitation of nonlinear oscillators and vortex-shedding process, several wake-oscillator models are proposed. Under this circumstance, the lift force coefficient is described by an oscillator integrating a linear spring with a nonlinear damping. The van der Pol equation or Rayleigh equation are widely used. They are quite similar but with different damping terms. One typical example is expressed as (Ogink et al., 2010): ω s ε where C y is the nondimensional lift force; , the vortex shedding frequency; C yA and are the amplitude and tuning parameter of the oscillator. The hydrodynamic force f(Y, t) due to structural motion reflects the coupling of structural motion and flow. The expressions of different orders of nonlinearity can be found in literature (Qu and Metrikine, 2020). By coupling the wake-oscillator with Eq. (18), the VIV response is ready to be yielded.
where, Y is the transverse amplitude; F VY , the total instantaneous vortex force; m a , the added mass; , the structural natural frequency, and , the damping. Since all parameters in wake-oscillator models need to be determined by experiments, the booming CFD-based estimation seems to be more promising to replace the wake-oscillator models.

Perspective
In summary, ocean engineering mechanics has made remarkable progresses to meet human's growing demand for marine petroleum and natural gas. Thus far people can provide important information in this regard such as environment parameters, external loads, and structure response etc. However, extreme marine environment prediction, free surface and interface issues, fluid -solid interaction, multiphase fluid transport and measurement etc. are still challenging. However, we prefer to merely focus on a few frontier issues mostly relevant to free surface/interface including TC induced extreme surface wave, sloshing of LNG, cavitation and bubble dynamics, vortex-induced motions and vibrations, in this review.
In the post-Paris protocol era, we are confronting various challenges due to climate change. For instance, TC intensity strengthening, polar ice extent and thickness reduction, fossil energy replacement by renewable ones etc. Under these circumstances, recent trends of ocean engineering from shallow water to deep sea, from low/middle to high latitude region, from fossil to renewable energy become more obvious. Therefore, the following research topics would be put on the agenda by researchers in this field: (1) the design of underwater structure and study on the stability of long pipe over uneven seabed for integrated deep-sea workstation; (2) partly icy sea surface wave prediction in polar area; constitutive relationship of sea ice, collision of floating ice and structure for polar exploration; (3) the measurement of hydrate mechanic and thermodynamic properties, stability of hydrate sea bed and sustainable exploitation approach of hydrate, marine wind power resource investigation, design of wind power turbines etc. Based on the previous fundamental research and engineering practice, ocean engineering exhibits a bright future.