Numerical Simulation of the Force Acting on the Riser by Two Internal Solitary Waves

Abstract

An internal wave is a typical dynamic process. As an internal wave, an internal solitary wave usually occurs between two layers of fluids with different densities. Compared with general internal waves, internal solitary waves have large amplitudes, fast propagation speeds, short-wave periods, and often have tremendous energy. The propagation causes strong convergence and divergence of seawater and generates a sudden strong current. Due to its various characteristics, the propagation of internal solitary waves can cause serious harm to offshore engineering structures. Therefore, studying the effects of internal solitary waves on risers is vital in preventing environmental pollution caused by riser damage. Although the research on internal solitary waves has achieved very fruitful results, the research on structures is mostly focused on a single condition, and the occurrence of internal solitary wave, as a complex ocean phenomenon, is often accompanied by many situations. Therefore, this paper constructs a numerical simulation of the interaction between two columns of internal solitary waves and risers. This study explores the force and flow field changes of the riser under the condition of multiple internal solitary waves using the Star-CCM+ software in the simulation. The improved K-epsilon turbulence model was adopted to close the three-dimensional incompressible Navier–Stokes equation, and the solitary wave solution of the eKdV equation was used as the initial and boundary conditions. The interaction between single and double internal solitary waves and a riser was calculated, compared, and analyzed using numerical analysis. The experiment results indicate that the conditions of two internal solitary waves differ from those of a single internal solitary wave. After colliding at the riser, the waves gradually merge into a single wave, and the flow field reaches its minimum velocity. Under the two-wave condition, the horizontal force on the riser as a whole is less than the single-wave condition. As the amplitude difference between the two internal solitary waves gradually decreases, the horizontal opposing force received by the riser first increases and then decreases, while the horizontal positive force gradually decreases.

1. Introduction

Internal waves are a relatively common dynamic process. Internal solitary waves, as a form of internal waves, generally occur between two layers of fluids with different densities. Due to its short period and other propagation characteristics, its wave has huge energy, it will cause the water flow to produce a fast bunching phenomenon during the propagation process and will produce a fast, sudden, strong speed. This will cause serious harm to marine engineering structures. Therefore, it is of great significance to study the effect of internal solitary waves on the riser to prevent damage to the riser from causing huge environmental pollution [1,2,3]. In 2014, a typical internal wave generated by strong tidal interaction was observed in the southeast of the Belcher islands. The research shows that the internal wave may be generated by the high tides of Hudson Bay passing through the narrow channel with hill-like features and a dense cline. At the same time, internal solitary wave movement was observed in this area [4]. To study the relationship between internal waves and heat, nutrient flux, and environment, some scholars use field survey and numerical simulation methods to analyze the motion of internal waves in different watersheds, capture key flow characteristics such as shock and flow, and focus on the internal wave breaking and watershed temperature change [5,6]. To research the difference between static and non-static effects on internal waves, some scholars use the nonlinear and high-precision numerical model with non-static pressure term to simulate the generation and propagation of short-period internal waves in Gibraltar [7]. The results show that there are obvious differences between static and non-static results, and the rupture of the waveform is more obvious. To better discuss and study the interaction among the atmosphere, internal waves, and surface waves, an ocean–atmosphere coupled prediction system is used to simulate the semidiurnal internal waves generated by the Andaman Sea archipelago. The model shows that the semidiurnal internal waves interact with the circulation and surface waves, which changes the flow and stratification in the ocean, and increases the temperature of the thermocline [8]. Some scholars use the CTD (conductivity, temperature, depth) measurement system and Doppler velocimeter to measure the vertical displacement of water quality points near some typical areas and observe the propagation of solitary waves to shallow water areas in a short period through satellite images [9,10].

The internal wave observation method in the nearshore area has been constantly innovated [11]. Some scholars have confirmed that the pressure gradient in the nearshore basin is very important to the internal wave propagation through the influence of the single point model and the idealized cross-section model on the generation of internal waves and the front on the edge of the continental shelf [12,13]. After that, some scholars studied the changes in wave velocity and energy conversion on the edge of the continental shelf through acoustic and microscopic data and infrared heat conduction maps. Through the combination of field observation and numerical simulation, some scholars have studied the dynamic characteristics of internal waves in tropical reservoirs with two basin structures in different flood seasons [14]. The results show that the internal waves in the basin are generated by the natural oscillation of wind, and the amplitude of internal waves in the two basins is different due to the different oscillation frequencies.

The study of internal waves on the resuspension of sediment at the bottom of the boundary, sediment transport, and boundary layer mixing is also a research focus [15]. A 2014 environmental assessment report in the waters of Gibraltar noted that the vertical mixing of ocean currents depends on the interaction of internal tides and short-period internal waves [16]. Internal wave motion and the collision process have also been described in detail; Boegman and Lamb pointed out that the mixing of background media appears at the stratification and bottom boundary of the watershed driven by turbulence and proposed that with the enhancement of the nonlinearity of an internal solitary wave, the requirements for accuracy in its simulation also increase greatly [17,18]. Steady stratified flow through hills with a Gaussian shape is also a focus in the field of environmental hydrodynamics. Some scholars use the high-resolution finite-difference numerical method to study Lee waves passing through smooth axisymmetric hills and discuss the effects of the Reynolds number and Froude number on the results [19].

Many theories and experiments have studied the causes of internal waves in the actual marine environment for a long time. Many experts and scholars have also conducted much tracking research on the occurrence of internal waves in different regions [20]. At present, many researchers have conducted extensive studies on the generation of solitary waves and their interaction with ocean structures [21]. A modified Morison equation that includes a convection term has shown that the force of a cylindrical structure could be better predicted by reducing the inertia coefficient and the resistance coefficient [22]. The large eddy simulation (LES) model can be used to study the force of internal solitary waves on a pair of cylinders arranged in series [23], which can then be used to understand internal solitary waves’ effect on different cylinder diameters, their evolution process, and influence [24]. Morison’s empirical formula and modal separation method estimate the force and torque exerted by internal solitons on cylindrical piles [25,26]. The horizontal velocity and acceleration on the vertical wave profile can then be estimated to simulate the force and torque on the riser. The waveform of internally isolated waves will change to different degrees after passing through a sloping topography with a certain angle. Keisuke Nakayama et al. first applied the Reynolds number of internally isolated waves (Re_ISW) to collapsing and plunging breakers and demonstrated that the density ratio of flow field stratification has less effect on waveform breaking [27]. Another study characterizes the evolution of internal isolated waves over slope terrain with different parameters such as slope and initial wave amplitude, and the results show that the wave amplitude increases and then decreases during the propagation process, and the maximum wave amplitude is at the shoulder of the slope [28].

For seasonal water layers, the load imposed on the riser by the inner solitary wave depends mainly on the variation in the water layer [29]. Numerical simulation has been used to investigate the interaction between internal solitary waves and transverse and extended cylindrical structures [30]. In another study, Open FOAM was used to numerically simulate the interaction between internal solitary waves in three cylinders and characterize the effect on the impact coefficient [31]. One study analyzing the stress at different sections of the riser under the action of the internal solitary wave using PIV technology found the variation trends for Cd and Cm hydrodynamic coefficients, where Cd decreases with an increasing Reynolds number, and Cm does not have a significant trend with an increasing Reynolds number [32]. The flow evolution of a cylinder with a single groove is also a focus of study. Different grooving angles have a significant effect on the fluid separation in the boundary layer, and the Strohal number (ST) reaches a maximum when the grooving angle θ = 90 degrees [33,34].

In addition, numerous studies have investigated internal solitary waves and offshore structures. Large gravity density stratified tanks have been used to analyze the interaction between an internal solitary wave and a tension leg platform and its load characteristics. Studies have found that the drag force coefficient in the Morison formula is in a power function relationship with the KC number [35,36]. One study established a time-domain numerical model to calculate the action of internal solitary waves on offshore structures and structural motion response [37]. The reaction force caused by mooring lines was analyzed using the conduit theory, while the motion equation for offshore structures was solved in time-domain by the fourth-order Runge–Kutta method. Another study constructed an internal solitary wave propagation model in a continuous density stratified ocean and simulated a numerical flume with a submerged dike [38,39]. One study experimentally investigated the effect of the splitter plate perpendicular to the waveform on the interaction between internal solitary wave and riser [40]. Semi-submersible platforms have also been used to simulate the nonlinear interaction between ISW and semi-submersible platforms. Studies have shown that the stress is mainly composed of wave pressure difference forces, viscous pressure difference forces, and frictional force [41]. Bricker conducted a numerical simulation on the interaction between an internal solitary wave and a floating tunnel. This research shows that changing the tunnel section shape can effectively reduce the vertical force and floating weight ratio [42].

Most studies have focused on the numerical simulation of a single internal solitary wave opposite a pipe. As a complex ocean phenomenon, the occurrence of internal solitary waves is often accompanied by many other phenomena. To address the current gaps in research, this paper numerically simulates the interaction between the symmetric and asymmetric internal solitary waves and the riser. Complex waveform changes and flow fields around the neutral and the force on the riser during relative motion are also explored and discussed.

2. Materials and Methods

2.1. Study Design and Data Collection

With the rapid development of modern computing technology, using computers combined with numerical algorithms to solve complex problems has increasingly become ubiquitous in many fields [43,44]. Hydrodynamic simulation software can discretize unsolvable nonlinear equations in space or time and solve the discretized equations to obtain the approximate solution. In this paper, the Simcenter Star-CCM+, a CFD simulation software, is used for data processing and calculations. It is a complete multi-physical field solution that can simulate the products and designs working under actual conditions. The software, which adopts the most advanced computational continuum mechanics algorithms, integrates CAD design, CFD solution, and result post-processing functions [45,46,47] It does not need multiple software links for coupling and can solve highly complex flow field models.

In this experiment, the specific parameters of the calculation model were entered into the CAD plate in Star CCM+, and the numerical flume was established. The boundary conditions, initial conditions, and appropriate physical parameters were then imported into the model to approximate the natural state. After meshing, the propagation of the internal solitary wave and the interaction between the internal solitary wave and riser were calculated using a separated and coupled finite-volume-method fluid solver. This step is called the CFD solution process. Finally, the results were exported in data form, and the vector picture with data was obtained through the post-processing function.

2.2. Geometric Model and Mesh Parameters

In internal solitary wave research, the commonly used simplification method involves simplifying seawater into two layers of fluid. The ocean, with density a thermocline, is regarded as two layers with uniform density, and a sudden density change can be observed between layers [48,49]. The two-layer density model assumes that the two-phase fluids are impermeable, non-interfering, and incompatible with each other. For this study, the numerical simulation domain of internal solitary waves was 6 m long, 0.35 m wide, and 0.57 m high. According to the two-layer density model, the calculation domain is divided into upper and lower layers. The height of the upper water body was h₁ = 0.095 m and the density was ρ₁ = 999 kg/m³; for the lower water body, h₂ = 0.475 m and ρ₂ = 1024 kg/m³. The dynamic viscosity coefficient for the two water layers was 1.003 × 10⁻⁶ pa-s, while the riser diameter was 8 cm. The center of the riser was 3 m from the left end of the calculation domain and 0.175 m from the front and back of the calculation domain. Figure 1 shows the specific structure.

To avoid the interference of grid size to the numerical simulation results, the grid independence needs to be verified first. Figure 2 shows the accuracy of locally encrypted mesh in the computational domain under three different conditions. The mesh generation uses cutting volume mesh and surface reconstruction.

Table 1. Specific parameters of grid division.

Mesh	Δx (m)	Δy (m)	Δz (m)	N
Sparse	0.0075	0.01	0.0075	1,019,752
Medium	0.005	0.01	0.005	3,171,556
Dense	0.0025	0.01	0.005	6,109,244

shows the grid size, and the overall grid number under three conditions, where Δx, Δy, and Δz represent the grid size in x, y, and z directions, respectively, and N represents the total number of flow field element grids. As shown in Figure 3, there was considerable difference in waveform change under sparse precision, while in the medium and dense precisions, the differences were small. Therefore, the medium grid precision was selected in this study.

Figure 4 shows the whole grid division. The foundation size was set to 0.05 m, and the grid size of the wave surface area used the medium parameter presented in Table 1. The Z-direction was used as the central axis at the center of the riser, and the dense cylindrical areas were set with diameters 0.04 m, 0.08 m, and 0.16 m. The grid sizes were 0.002 m, 0.004 m, and 0.008 m. In this way, the stepped encryption method reduces the number of overall grids, speeds up the calculation time, and does not distort the flow field near the riser that could result in erroneous results.

2.3. Mathematical Models

2.3.1. Governing Equation

Using the Navier–Stokes equation, the propagation of the internal solitary wave and the action of the riser were numerically simulated. The governing equations for the flow field are as follows:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho \overline{u_i}\right)}{\partial x_i}=0$$

(1)

$$\frac{\partial \rho }{\partial t}=0$$

(2)

$$\frac{\partial }{\partial t}(\rho \overline{u_i})+\frac{\partial }{\partial x_j}(\rho \overline{u_i}\overline{u_j})=-\frac{\partial \overline{P}}{\partial x_i}+\mu (\frac{{\partial }^2\overline{u_i}}{\partial x_j\partial x_j})+\frac{\partial (-\rho \overline{u_i^{\text{'}}{u}_j^{\text{'}}})}{\partial x_j}+\rho g_i$$

(3)

$$-\rho \overline{u_i^{\text{'}}{u}_j^{\text{'}}}={\mu }_t(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i})-\frac{2}{3}({\mu }_t\frac{\partial \overline{u_i}}{\partial x_i}){\delta }_{ij}$$

(4)

$${\mu }_t=\rho C_{\mu }\frac{k^2}{\epsilon }$$

(5)

$$\frac{dk}{dt}=\frac{\partial }{\partial x_j}\left[\left(c_k\frac{k^2}{\epsilon }+v\right)\frac{\partial k}{\partial x_j}\right]+P_k-\epsilon$$

(6)

$$\frac{d\epsilon }{dt}=\frac{\partial }{\partial x_j}\left[\left(c_{\epsilon }\frac{k^2}{\epsilon }+v\right)\frac{\partial \epsilon }{\partial x_j}\right]+c_{\epsilon 1}\frac{\epsilon }{k}{P}_k-c_{\epsilon 2}\frac{{\epsilon }^2}{k}$$

(7)

where i and j are three directional components in the Cartesian coordinate system; t is the time, p is the fluid density; µ is hydrodynamic viscosity coefficient; F_i is the volume force; P_k is the turbulence generation term; and c_µ, c_k, c_ε, c_ε1, and c_ε2 are empirical constants.

2.3.2. Internal Isolated Wave Theory

The common inner isolated wave theory solutions are KdV, mKdV, eKdV, and MCC theories, and Cui et al. have investigated the applicability of different theories based on the inner isolated wave dispersion and nonlinear parameters. Based on the established dimensions of the simulated computational domain, the eKdV theory is chosen as the theoretical equation to describe the internally isolated waves in this paper [50].

The full name of the eKdV equation is the extended Korteweg–de Vries equation. It includes nonlinear and dispersive cubic terms in the wave amplitude derived from the water–wave equations and the Lagrangian for the water–wave equations. Since there is not necessarily a balanced relationship between nonlinearity and dispersion, the KdV equation does not apply to isolated waves within a large amplitude, where a cubic nonlinear term can be introduced, i.e., a nonlinear term of significant order is added after KdV to obtain the eKdV equation:

$$\frac{\partial H}{\partial t}+(c+{\alpha }_{0}H+{\alpha }_1{H}^2)\frac{\partial H}{\partial x}+\beta \frac{{\partial }^{3}H}{\partial x^3}=0$$

(8)

Of which:

$${\alpha }_1=\frac{3c}{h_1{}^2{h}_2{}^2}\left[\frac{7}{8}{\left(\frac{{\rho }_1{h}_2{}^2-{\rho }_2{h}_1{}^2}{{\rho }_1{h}_2+{\rho }_2{h}_1}\right)}^2-\frac{{\rho }_1{h}_2{}^3-{\rho }_2{h}_1{}^3}{{\rho }_1{h}_2+{\rho }_2{h}_1}\right]$$

(9)

Equations (2)–(10) has the following inner isolated wave solution, called the eKdV theoretical solution:

$$H(x,t)=\frac{H_0}{B+(1-B){\cosh }^2\left[\frac{\left(x-c_{eKdV}t\right)}{l_{eKdV}}\right]}$$

(10)

Of which:

$$c_{eKdV}=c+\frac{H_0}{3}\left({\alpha }_0+\frac{1}{2}{\alpha }_1{H}_0\right)$$

(11)

$$l_{eKdV}{}^2=\frac{24\beta }{H_0\left(2{\alpha }_0+{\alpha }_1{H}_0\right)}$$

(12)

$$B=\frac{-H_0\beta }{2{\alpha }_0+{\alpha }_1{H}_0}$$

(13)

$c_{eKdV}$ is the phase velocity of the internal isolated wave calculated with theory. The theory corresponds to the extreme value of the internal isolated wave amplitude as:

$$H_{\max }^{eKdV}=\frac{4h_1{h}_2(h_1-h_2)}{h_1^2+h_2^2+6h_1{h}_2}$$

(14)

2.3.3. Turbulence Method and Liquid Level Capture

In this paper, the improved realizable K-epsilon two-layer model (RKE—2L) is selected to close the Reynolds time average equation [51]. This method was first proposed by Rodi. The realizable K-epsilon model introduces the variable damping function f_μ (expressed as turbulence and uniform flow properties) to apply to the critical coefficient c_μ of the model.

Therefore, the turbulent eddy viscosity µ_t can be expressed as:

$${\mu }_t=\rho C_{\mu }{f}_{\mu }\frac{k^2}{\epsilon }$$

(15)

$$f_{\mu }=\frac{1}{C_{\mu }\left\{4+\sqrt{6}\cos \left[\frac{1}{3}{\cos }^{-1}\left(\sqrt{6}\frac{S^{{*}^3}}{\sqrt{S^{*}:S^{*}}}\right)\right]\frac{k}{\epsilon }\sqrt{\mathrm{S}:\mathrm{S}+\mathrm{W}:\mathrm{W}}\right\}}$$

(16)

$$S^{*}=\mathrm{S}-\frac{1}{3}\mathrm{tr}(\mathrm{S})\mathrm{I}$$

(17)

where S represents the average strain rate tensor; W represents the mean vortex tensor; I represents the identity matrix.

In this way, the numerical model and the realizable physical model can have the same mathematical constraints, so the results of this model are better than the standard K-epsilon model in many applications. In this paper, the fluid volume method (VOF) is selected as the liquid level tracking method. The fluid volume (VOF) multiphase model is suitable for simulating multiple immiscible fluid flows on the numerical grid (which can solve the interface between the phases of the mixture). The model is suitable for each phase with a large structure and relatively small contact area between each phase.

The traditional high-order schemes, such as the central difference scheme and the upwind scheme, have poor calculation results when the phase volume fraction changes greatly. In this paper, the high-resolution interface capture (HRIC) method is used to track the interface between the two phases, which can effectively reduce the error caused by the liquid level capture.

The phase distribution and position of the interface are described by the field of phase volume fraction α_i. The volume fraction of the phase is defined as follows:

$${\alpha }_i=\frac{V_i}{V}$$

(18)

where V_i is the volume of phase I in the grid cell, and V is the volume of the grid cell. The sum of the volume fractions of all phases in the grid cell must be 1:

$${\displaystyle \sum _{i=1}^N{\alpha }_i=1}$$

(19)

N is the total number of phases. According to the value of volume fraction, whether there are different phases or fluids in the grid cell can be distinguished.

• α_i = 0 represents that the grid element has no phase i at all;
• α_i = 1 represents that the grid cell is filled by phase I;
• 0 < α_i < 1 represents that the grid element has an interphase interface.

2.3.4. Numerical Methods and Boundary Condition

The initial conditions and the left and right boundaries were determined using the theoretical solution H (x, t) of the internal solitary wave of the eKdV equation [52]. The relationship between velocity and H (x, t) can be expressed as:

$$\overline{u_1}=-\frac{{\mathrm{c}}_{ekdv}H(x,t)}{h_1}$$

(20)

$$\overline{u_2}=\frac{{\mathrm{c}}_{ekdv}H(x,t)}{h_2}$$

(21)

The left boundary was set to x = 0, while the right boundary velocity condition was set similar to the left but in the opposite direction. Since the amplitude scale of internally isolated waves is much larger than that of surface waves, the study of internally isolated waves often assumes that the upper surface is an undeformable wall, called the “rigid cover assumption.” In this paper, we consider the computational domain for the force of the riser pipe running through the bottom of the sea, while doing comparison with the experimental water tank results, the bottom boundary is also set as the wall boundary. The front and rear boundaries were set as standard, the velocity derivative was zero, and the riser surface was set as a wall, speed at the wall is 0 The initial internal solitary wave positions are located at the left and right boundaries. The fluid is incompressible, the interface between the upper and lower layers uses the fluid volume method (VOF), and the turbulence model is the K-epsilon model. The characteristic wave parameters are summarized in

Table 2. Internal solitary wave parameters.

Working Condition	Left Amplitude, η0 (cm)	Right Amplitude, η0 (cm)
1	7.9	None
2	7.9	3.2
3	7.9	4.2
4	7.9	5.5
5	7.9	6.5
6	7.9	7.9

.

2.3.5. Model Validation

Numerous studies have dealt with and analyzed the problems and limitations in model validation. In this paper, the same parameters used by Wang Fei (Table 2—working condition 1) were selected, and the calculations and experimental results were compared to their findings [53]. Figure 5 compares the numerical simulation waveform recorded at the detection point at x = 1 m and the theoretical waveform. In the figure, the y-axis is the numerical change in the two-phase junction point in the Z-direction, and the x-axis is time. To compare with the literature results, the riser was divided into six sections, each at 9.5 cm high and numbered 1–6 from top to bottom (see Figure 6). The horizontal stress for each riser was analyzed, and the results are presented in Figure 7. In the previous study, data for the No. 3 riser were missing due to equipment failure.

The comparative results show that the simulated horizontal force on riser section 1 was much smaller than the experimental value. In the actual experiment, there was air interference in the upper layer of the liquid level. However, in the numerical simulation, the top boundary of the calculation domain was assumed to be the wall boundary based on the rigid cover, and there was no surface disturbance. Therefore, the horizontal force for riser 1 obtained by numerical simulation was less than that obtained by experiment.

The results show that this paper’s numerical method’s solitary internal waveform is comparable with the theoretical solution waveform. The simulated solution peak value is slightly smaller than the theoretical solution peak value, but the characteristic width is consistent with the theoretical value within a reasonable error range. The overall trend of horizontal stress of each riser is comparable with the experimental and simulation results in the literature. In conclusion, the numerical simulation from this study produced small errors and more accurate results.

2.3.6. Numerical Methods and Boundary Condition

Figure 8 shows the waveform outline for the three working conditions. The left image shows the waveform transmitted to the riser, while the right image shows the waveform transmitted through the riser. As shown in the figure, the propagation of the internal solitary wave for conditions 1 and condition 3 is affected by the riser, resulting in the unsmooth waveform. Under condition 6, the waveform is hardly affected by the riser, and the waveform remains smooth and symmetrical. In the wave collision scenario, the waves propagate to each other before t = 30 s, while the inner sides of the waves continue to fuse. The wave amplitude increases and gradually fuses into a single stable wave, and the single wave amplitude is slightly less than the sum of the two internal solitary wave amplitudes before fusion. After a short fusion, the amplitude of the single wave gradually decreases while the wavelength steadily increases. Finally, it separates into two internal solitary waves, which move to each other. The separated internal solitary wave amplitude is almost the same as before the collision.

3. Results

3.1. Flow Field Analysis

Figure 9 shows the variation in the wave surface near the riser over time for the three working conditions. For conditions 1 and 3, when the wavefront passes through the riser, the flow surface on the riser’s right side appears as folds on both wings. With the propagation of the internal solitary wave, the back end propagates to the riser, and folds appear on both sides of the riser’s upstream surface. The wave surface change in condition 1 is more evident than in condition 3. For condition 6, when the front end of the two waves passes through the riser, the riser’s wave surface is smoother than that of the single wave condition. The wave surface height around the riser decreases continuously and peaks at t = 30 s. At this time, the wave height of the upstream surface on the left and right sides of the riser is almost the same. With the propagation of the internal solitary wave, the height of the flow field near the riser increases continuously and returns to the horizontal state.

Given that the flow field under the combined action of two internal solitary waves is complex, the section velocity for the three working conditions is further analyzed. The calculation domain cuts along the center of the riser in the x-z direction, as shown in Figure 10.

From Figure 10, whenever the internal solitary wave is layered up and down along the wave surface, the velocity direction is opposite. For the single-wave scenario (condition 1), the front end of the solitary inner wave propagates to the riser at t = 28 s, and then a vortex appears on the right side of the riser. After the trough propagates through the riser in the t = 30 s, the vortex on the right side of the riser disappears, and the vortex appears on the left side. As the wave trough propagates away from the riser, the vortex phenomenon increases on the left side.

For condition 3, the collision between the two waves is located on the right side of the riser. At t = 28 s, eddy currents appear at the bottom right side of the collision between the two waves; at t = 30 s, they also appear at the collision point. The left wave still dominates the force. The two waves separate at 32 s. The left wave moves away from the riser faster than the right wave, and the vortex becomes smaller and moves closer to the interface.

For condition 6, the collision of two or three waves is located at the riser, and a slight vortex phenomenon appears in the flow field near the bottom of the riser. At t = 30 s, as the trough of the solitary inner wave propagates to the riser, the velocity is smaller than at t = 28 s. At the same time, the eddy current in the flow field near the bottom of the riser is more serious. At t = 32 s, the two waves separate and move away from the riser, and the upper and lower velocities still follow. The wave surface is stratified, the velocity is almost unchanged when the collision occurs (t = 28 s), and the vortex phenomenon at the bottom of the riser disappears.

3.2. Riser Force Analysis

The stress of the riser under the conditions of the two kinds of internal solitary waves was then analyzed.

$${\mathrm{f}}_f^{shear}=-T_f·a_f$$

(22)

$a_f$ represents the area vector of surface grid, $T_f$ is the stress tensor on the surface.

The amplitude of the internal solitary wave on the left side of the calculation domain is 7.9 cm, and the amplitude of solitary waves on the right side increases from zero to 7.9 cm. The comparative analysis of the riser’s horizontal direction (horizontal) for the different conditions is shown in Figure 11. In the image, the left wave propagation is presented in the positive direction, while the right wave propagation is at the opposite direction. The horizontal force of the riser under the combined action of the two waves is less than that of the standpipe under the single-wave condition—opposing force. The riser receives an opposing X force first, then a positive X force, and finally a tiny opposing X force. In Figure 12, when the single wave condition and the amplitude difference on both sides are small, the change law of the riser torque curve is like that of the horizontal force curve. The torque direction is the first negative, then positive, and finally negative, with two deflections in the direction. After working condition 3, the riser is only affected by the negative torque, and the peak value occurs at about t = 30 s, that is, before the two trains of waves merge into a stable single wave. With the continuous reduction in the amplitude difference between the two trains of waves, the occurrence time of the peak value of the torque goes backwards, and the extreme value of the torque also decreases.

The extreme absolute force values (Fxmax) were compared for the different working conditions. In Figure 13, the opposing force in the first section is line 1, the positive force in the second section is line 2, and the opposing force in the third section is line 3. Under the combined action of two internal solitary waves, the increasing amplitude on the right side causes the opposing force on the riser in the horizontal to first increase and then decrease, while the positive force decreases continuously. The amplitude of the first negative force on the riser is greater than that of the second positive force on the combined action of two waves. For the section 3, the opposing force amplitude is small. Here, the force on the riser is dominated by the right wave since the propagation of the left wave is faster than that of the right wave.

After analyzing the magnitude of the extreme value, the occurrence time of the extreme value is compared, to analyze the impact of the collision time on the force of the opposite tube. The amplitude difference of incident waves on both sides decreases continuously, resulting in the forward movement of the fusion time of the two waves. The time of the first peak gradually moves backwards with the increasing amplitude on the right side. When the amplitude on the right side is small, the time goes faster. After the amplitude on the right side reaches 4.2 cm, the time of the peak tends to be gentle. Peak 2 and peak 3 move forward with the increase in wave amplitude on the right, and the linear relationship between the passage time of peak 3 and the change in wave amplitude is obvious. With the decrease in amplitude difference between the two sides, the time of the riser subjected to negative horizontal force increases and the time of the riser subjected to positive horizontal force decreases. The following figure shows the stress comparison of riser sections at different times under conditions 2, 3, and 4. The ox axis represents the stress of the riser, in which the positive direction is the propagation direction of the incident wave on the left, the negative direction is the propagation direction of the incident wave on the right, and the oy axis is the distribution of riser sections along with the water depth.

In Figure 14, Figure 15 and Figure 16, We can see from the figure above that the stress distribution of risers 2, 3, and 4 under the 28 s condition is the same as the horizontal velocity distribution trend of a single internal solitary wave. The riser at the upper layer of the interface is subject to X positive force, and the riser at the lower layer of the interface is subject to X negative force. Due to the fast phase velocity of the incident wave on the left, the force on the opposite pipe plays a leading role. The greater the difference in wave amplitude, the greater the force on each section of the riser. At t = 31 s, the two trains of waves will converge. Currently, the horizontal force direction of each section of the riser is opposite to that at t = 28 s, and the horizontal force amplitude of section 2 of the riser exceeds that of section 1 of the riser. The force amplitude increases with the decrease in the amplitude difference between the two trains of waves. At t = 34 s, the single wave gradually separated, and the horizontal force direction of each section of the riser did not change compared with t = 31 s. The horizontal force on the first section of the riser increased, and the horizontal force on the other sections of the riser decreased. It can be further learned from the above analysis that under the joint action of two trains of waves, the horizontal force of each section of the riser becomes smaller as a whole compared with the single wave working condition, but the stress time is prolonged, the joint of section 2 and section 3 of the riser is repeatedly subjected to shear action, and the section 1 of the riser is subjected to horizontal force with changed direction and large amplitude.

The results suggest that the riser’s overall stress is minimal under the combined action of two internal solitary waves. Further analysis was then carried out by comparing the shear stress of the opposite riser. Figure 17 shows the comparison diagram of the wall shear stress of the riser at t = 30 s in the horizontal and vertical directions for the three working conditions. Under condition 1, at t = 30 s, the riser is stratified by sheer force as a whole. The upstream surface of the upper layer is subjected to the shear force in the positive direction of the x-axis. The single wave propagation and the lower layer are subjected to shear force in the negative direction of the x-axis. This means that the opposite direction of the single wave propagation and the absolute value of the shear force on the upper layer is greater than the absolute value of the shear force on the lower layer. The wall in the vertical negatively shears the upstream surface on the wave surface of the riser, and the wall in the vertical positively shears the upstream surface on the wave surface. Condition 3 has a similar shear force trend with condition 1 but has smaller amount of stress and more extensive stress area on the upper layer.

Under condition 3, due to the influence of the relative motion of two waves with the same amplitude, the wall shear stress in the horizontal and vertical is small, almost negligible. This is consistent with the horizontal stress result of the riser.

4. Discussion

Based on the N–S equation, this study explores the flow field and wave surface of the opposite pipe for different internal solitary wave conditions. The analysis shows that the interaction of multiple internal solitary waves with the riser can be accurately estimated using numerical simulation, laying the foundation for further research. The main simulation results are as follows:

• Unlike a single internal solitary wave, after the two internal solitary waves collide at the riser, they gradually fuse. When they finally fuse into a single wave, their amplitude reaches the maximum and is slightly less than the sum of the amplitudes of the two waves. At that point, the velocity of the flow field around the riser reaches a minimum. When the two waves act together, the change in the wave surface tends to be gentle. When the amplitude difference between the two sides is less than a certain degree, the wave surface will not deform;
• The force on the riser in the horizontal and vertical directions is more significant than that in the z direction, and the sectional force on the riser in the upper fluid area is greater. Overall, with the continuous reduction in the difference in amplitude, the stress direction, extreme value, and occurrence time of the extreme value in the horizontal and vertical directions will change. The overall force of the riser under the joint action of two waves is less than that under the condition of a single wave, and when the amplitude difference of waves on both sides is large, the riser will be subjected to an additional section of X negative force with small amplitude after being subjected to X positive force;
• As the amplitude difference between the two sides decreases, the amplitude and time of the positive X force on the riser decreases; the X negative force first increases and then decreases, and the duration of the negative force increases. The horizontal force of the riser section changes before and after the collision and separation of two waves. At the same time, the force amplitude of the riser section is affected by the amplitude difference. When a stable single wave is generated in the flow field, the amplitude difference on both sides has a negative correlation with the horizontal force of the riser section; Before and after single wave fusion and separation, the amplitude difference between the two sides is positively correlated with the horizontal force of the riser section.

5. Conclusions

Although the mesh dependence of the waveform is verified in this paper, a finer mesh is still needed to capture the vortex phenomenon of two-wave collision more accurately. The Courant number of the flow field stays within 0.001 in most of the calculated area, but the Courant number increases locally near the wall of the riser. Different turbulence calculation methods may affect the results to some extent, and the results obtained by different methods should be compared and analyzed in future analysis to further enhance their credibility.

The results suggest that the proposed numerical simulation of the interaction between multiple internal solitary waves and the riser is reliable and provides accurate estimates and that the influence of right-side waves on the flow field is significant. Unlike the traditional concept of wave collision canceling each other, the force on the riser is not apparent when the amplitude of the two internal solitary waves is the same but in opposite directions. When the amplitude difference of the two internal solitary waves is significant, it has a positive impact on the opposite tube. While the interface change in the two-wave collision fusion is small, it does not mean that the change in the whole flow field is not complex. The vortex phenomenon generated by the collision of the two waves in the lower fluid area significantly impacts the action of the riser. As the amplitude on the right continuous to increase, the collision of the two waves gradually approaches the riser. In this process, the force on the riser increases first and then decreases.

6. Limitations

There were several limitations in this study, given that the numerical simulation was implemented based on the condition that the riser structure is rigid. In practical engineering problems, the material properties of oil risers are usually polymer composites, which can deform after being impacted by waves in the marine environment. Subsequent experiments and simulations should explore the fluid distribution and force deformation of flexible risers interacting with multiple internal solitary waves.