Analysis of the Descent Process and Multi-Objective Optimization Design of a Benthic Lander

Abstract

The growing need for deep-sea biological research and environmental monitoring has expanded the demand for benthic landers. Compared with remotely operated vehicles (ROVs) and autonomous underwater vehicles (AUVs), benthic landers can reduce overall operation cost and also possess longer endurance. Configuring a suitable descent velocity is important for benthic lander designs, helping them avoid retrieval failure and improve sea trial efficiencies. In this study, an effective scheme for the configuration and optimization of a self-developed benthic lander was outlined. First, the structural characteristics of the benthic lander were analyzed, and then a dynamic model was established. Second, the hydrodynamic coefficients of the benthic lander during its descent process were calculated using computational fluid dynamics (CFD) methods. Third, the MATLAB Simulink simulation environment was used to solve the dynamic model, and then the multi-objective optimization algorithm was introduced for the optimization design. Finally, the model was validated based on sea trial data, which demonstrated that the designed configuration and optimization scheme were correct and efficient. Collectively, this work provides a useful reference for the rational configuration and practical application of benthic landers.

1. Introduction

The deep sea is rich in energy, minerals and biogenetic resources; thus, it has significant economic value [1,2]. The impact of human activities on the marine environment is continually increasing [3,4], especially in offshore areas, but the response of the deep-sea to human activities remains unclear. In recent years, a large number of underwater vehicles have been used in marine climate change exploration, deep-sea observation, bioprospecting, mineral resource development and other fields. However, the mobility of conventional ROVs is severely limited due to the connection of the umbilical cable, while the cost of ROV sea trials is relatively high [5,6]. Although AUVs compensate for the abovementioned ROV shortcomings, they have poor endurance [7]. In contrast, benthic landers are playing an increasingly important role due to their relatively low cost, longer endurance and reusability.

Benthic landers have been widely utilized in in situ studies of deep-sea ecology, resources and environments. In 1975, Smith et al. designed a free vehicle respirometer (FVR) that could in situ measure the oxygen consumption of benthic communities relative to abyssal depths [8]. The other two types of benthic landers—the Autonomous Lander for Biological Experiments (ALBEX) [9,10] and the deep-sea benthic environmental observation system (Benvir) [11]—were developed for studies on deep-sea sediment community oxygen consumption. With the development of marine science and engineering technology, benthic landers can be equipped with various types of sampling modules for the collection of deep-sea seawater, sediment, microfauna and macrofauna. In 2009, the sediments incubated in the benthic lander named FLUFO [12] contained seep-associated fauna, and this study reveals the relationship between methane emissions from sediments with the seep-associated fauna. In the same year, a project called the Hadal Environment and Educational Program (HADEEP) was launched; it aimed to examine hadal biology and extend ecology research [13]. In this project, benthic landers were equipped with biological traps and seawater samplers to collect hadal samples. Thereafter, two function-similar benthic landers—a 7000-m-depth-rated lander (HaiJiao) [14] and a full-ocean-depth-rated lander—were separately produced by the Shenyang Institute of Automation of Chinese Academy of Sciences and the University of California San Diego (UCSD) [15]. Recently, Wei et al. successfully mounted their newly designed multiple in situ nucleic acid collections (MISNAC) on a benthic lander (Phoenix) to collect deep-sea microbial nucleic acid samples under in situ conditions [16]. In our previous study, a novel genus-level deep-sea bacterial species was isolated from an in situ enriched consortium collected by a benthic lander [17], and the lander used there is the prototype of the lander intended for optimization in this study. Benthic landers also have important applications in deep-sea acoustics [18], geohazards [19] and ocean-bottom seismometers (OBS) [20]. Collectively, in order to have an intuitive understanding of the background, onboard sensors and sampling modules of related landers, some typical cases of deep-sea environmental and biological studies using benthic landers are summarized in

Table 1. Summary of typical benthic landers used in in situ studies of deep-sea ecology, resources and environments.

	FVR	ALBEX	FLUFO	HADEEP	Benvir	HaiJiao	UCSD-Lander
Background
Country	U.S.A.	Netherland	Germany	Japan and UK	China	China	U.S.A.
Time	1975	1998	2009	2009	2009	2015	2019
Sensor
CTD		✓	✓	✓	✓	✓
Dissolved Oxygen (DO)	✓	✓			✓	✓
PH					✓
Sampling modules
Camera	✓	✓	✓	✓		✓	✓
Seawater	✓			✓		✓	✓
Sediments		✓	✓			✓
(Microbe Enrichment)						✓
Biological trap				✓		✓	✓

.

In recent years, the frame structure of the benthic lander has also been greatly developed. The early lander mainly comprised a tripod frame [8,9,13] was constructed around a centrally located acoustic release, similarly to the FVR lander; on top of the three legs is an aluminum plate. The BIGO [21,22] lander also has a tripod frame, but a cuboid frame is added on the top of the floatation sphere array of the lander to carry more scientific devices. The three legs make up the simplest stable structure, and the plate provides a larger contact area to avoid retrieval failure caused by excessive subsidence depth with the seafloor. Then, the three legs gradually developed into four legs [23,24,25] for higher stability, which greatly increased the ability of the lander to carry scientific devices, similarly to the RAP2 [26] lander. The four-legged Gothenburg [27] lander adopted a two-frame solution, and this was significant because it allowed for the flexibility to deploy only the inner-frame from small vessels. However, landers used for ocean-bottom seismic research usually have flat boards [28] or annular structures [20,29].

The control of descent velocity is a key concern in the design and application of benthic landers. After the benthic lander separated from the ship, its descent velocity was mainly determined by the system’s configuration [30]. When the descent velocity is too low, the efficiency of the sea trial will be reduced, and the lander may also experience uncontrolled lateral drift driven by ocean currents. Simultaneously, excessive descent velocities will also lead to a series of problems, such as lander damage or sinking too deep into the seabed sediment, which may all lead to a failure in lander retrieval. To solve these problems, Mortensen et al. [31] designed crossed wings to control the descent velocity of the lander. However, the velocity was too low, which resulted in lower efficiency and higher costs in sea trials. Jun [32] and Gang et al. [33] controlled the descent velocity by integrating a hydrofoil into the lander, but the effectiveness needs to be further verified. Therefore, achieving a more direct control scheme over the descent process is a problem that needs to be considered when designing a benthic lander.

Based on the above analysis, an effective scheme for the configuration and optimization of a self-developed benthic lander was outlined in this study. This scheme started with dynamic modelling and hydrodynamic coefficients calculation. Second, the CFD methods and multi-objective optimization algorithms were used to comprehensively study the influence mechanism of several variables on the lander’s descent velocity, including the number of floatation spheres, the mass of the weight stack and the contact area of the weight stack. Finally, the mutual restriction relationship between the lander’s descent time and subsidence depth was revealed. This study provides a useful reference for the rational configuration and practical application of a benthic lander.

2. Model Establishment of the Benthic Lander

2.1. Frame Structure and Working Principle of the Benthic Lander

The lander mainly comprises three parts: an upper floatation sphere array that provides positive buoyancy, a frame that supports the scientific device and a lower weight stack that provides descending gravity (see Figure 1a). The floatation sphere array is connected in series with the fiber rope, and each floatation sphere (model: NMS-FS-6700-17RO, Nautilus Marine Service GmbH) provides 26 kg of positive buoyancy; moreover, the number of floatation spheres is configurable depending on the weight of the scientific device. The frame is made of titanium alloys. Loaded within the frame is a syntactic foam that provides positive buoyancy and structural support to the frame, and on top of the frame is a lifting bail that facilitates the process of deployment and retrieval. At the center of the frame is an acoustic release (model: OCEANO 2500S Light, iXblue, Paris, France) that is responsible for releasing the disposable weight stack. The scientific device includes the microbe enrichment device [34] array, which is the main usage of the lander, and a dual-axis inclinometer sensor (model: LCT100A, Guigangshi Lecheng Information Technology Co., Ltd., Guigang, China) with an angle resolution of 0.0001° and a 1–100 Hz configurable output frequency. The weight in water of the frame (including foam, release and science devices without the weight stack) is −85 kg, which was obtained from the water tank experiment (see Figure A1). Based on the above explanation of the lander’s frame structure, the working principle of the benthic lander is described as follows: At the end of the experiment, an acoustic remote-control unit (model: TT801, iXblue, Paris, France) is used to release the hook of the acoustic release; see Figure 1b. Then, the weight stack and bearing bracket (see Figure A2) are automatically separated from the frame, and the positive buoyancy provided by the floatation sphere array and syntactic foam allows the frame to ascend from the seafloor.

2.2. Benthic Lander Dynamic Modeling

The dynamic model of an underwater vehicle in six degrees of freedom can be described as a nonlinear equation of motion in matrix form [35]:

$$M\dot{\upsilon }+C\left(\upsilon \right)\upsilon +D\left(\upsilon \right)\upsilon +g\left(\eta \right)=\tau$$

(1)

$$\dot{\eta }=J\left(\eta \right)\upsilon$$

(2)

where $M=M_{RB}+M_A$ is the inertia matrix for the rigid-body mass matrix and added mass matrix, respectively;

$C\left(\upsilon \right)=C_{RB}\left(\upsilon \right)+C_A\left(\upsilon \right)$ is the Coriolis and centripetal matrix for rigid-body mass and added mass, respectively;

$D\left(\upsilon \right)=D_{quad.}\left(\upsilon \right)+D_{lin.}$ is the quadratic and linear damping matrix, respectively;

$g\left(\eta \right)$ is the hydrostatic restoring force matrix resulting from the different positions of the center of buoyancy and the center of gravity;

$\tau$ is the sum of the acting forces and moments produced by the propulsion system of the vehicle, which is equal to zero in this case because the lander is a freely descending vehicle;

$\upsilon$ is the generalized velocity vector; accordingly, $\dot{\upsilon }$ is the generalized acceleration vector;

$\eta$ denotes the generalized position coordinates, and $J$ is the rotation matrix for the transformation from the body-fixed frame to the earth-fixed frame. Note: $M,C\left(\upsilon \right),D\left(\upsilon \right)\in {\mathrm{R}}^{6\mathrm{x}6}$, $g\left(\eta \right),\tau ,\upsilon ,\eta \in {\mathrm{R}}^{6\mathrm{x}1}$.

Due to the relatively low descent velocity (please refer to the data presented in “6 Sea Trials”.), the Coriolis and centripetal matrix can be ignored, and its symmetrical structure design means that the lander is relatively stable with respect to its attitude (please refer to the data presented in “3.2. Calculation of Damping Force Coefficients”.); thus, the quadratic and linear viscous coupling coefficients between the degrees of freedom are small and can be ignored. Finally, the descent process of the next analysis was only carried out for the lander’s heave degree of freedom, and the model described in Equations (1) and (2) can be simplified as follows:

$$\left({\mathrm{M}}_{\mathrm{RB}}+{\mathrm{M}}_{\mathrm{A}}\right)\dot{\upsilon }+({\mathrm{D}}_{\mathrm{quad}.}\left(\upsilon \right)+{\mathrm{D}}_{\mathrm{lin}.})\upsilon =-\mathrm{g}$$

(3)

where ${\mathrm{M}}_{\mathrm{RB}}$ is the mass of the lander; ${\mathrm{D}}_{\mathrm{quad}.}$ and ${\mathrm{D}}_{\mathrm{lin}.}$ are the quadratic and linear damping coefficients, respectively; υ is the descent velocity in the heave degree of freedom; −g is the net gravity for the free descending vehicle (lander). Moreover, gravity is composed of three parts: weight in seawater of the weight stack (positive), the buoyancy of the floatation sphere array (negative) and the buoyancy of the frame (negative).

3. Calculating Hydrodynamic Coefficients

3.1. Hydrodynamic Modeling

The benthic lander is an open frame structure with different configurations of the scientific device, which increases the difficulty of solving the hydrodynamic coefficients in Equation (3). The Favre-averaged Navier–Stokes (FANS) equation is a widely used method for hydrodynamic simulations of underwater vehicles with turbulent flows. In this study, the CFD software add-in, SOLIDWORKS Flow Simulation, is used to implement the FANS to compute the hydrodynamic coefficients. Assuming that the fluid comprises continuous media and does not exhibit energy conservation, the conservation laws for mass and angular momentum can be written in the conservation form as follows [36]:

$$\frac{\partial \mathsf{\rho }}{\partial \mathrm{t}}+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left({\mathsf{\rho }\mathrm{u}}_{\mathrm{i}}\right)=0$$

(4)

$$\frac{\partial {\mathsf{\rho }\mathrm{u}}_{\mathrm{i}}}{\partial \mathrm{t}}+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{j}}}\left({\mathsf{\rho }\mathrm{u}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{j}}\right)+\frac{\partial \mathrm{p}}{\partial {\mathrm{x}}_{\mathrm{i}}}=\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{j}}}\left({\mathsf{\tau }}_{\mathrm{ij}}+{\mathsf{\tau }}_{\mathrm{ij}}^{\mathrm{R}}\right)+{\mathrm{S}}_{\mathrm{i}} \mathrm{i}=1,2,3$$

(5)

where $\mathrm{u}$ is the fluid velocity, and $\mathsf{\rho }$ is the fluid density. Moreover, ${\mathrm{S}}_1$ denotes a mass-distributed external force per unit mass due to a porous media resistance, and ${\mathrm{S}}_2$ denotes buoyancy (${\mathrm{S}}_2=-{\mathsf{\rho }\mathrm{g}}_{\mathrm{i}}$, where ${\mathrm{g}}_{\mathrm{i}}$ is the gravitational acceleration component along the i-th coordinate direction). ${\mathrm{S}}_3$ denotes the coordinate system’s rotation. Moreover, ${\mathsf{\tau }}_{\mathrm{ij}}$ denotes the viscous shear stress tensor, accordingly, ${\mathsf{\tau }}_{\mathrm{ij}}^{\mathrm{R}}$ denotes the rotation form of ${\mathsf{\tau }}_{\mathrm{ij}}$. Significantly, with the exception of ${\mathrm{S}}_{\mathrm{i}}$, other subscripts represent the corresponding coordinate direction.

The κ-ε turbulence model is used for the simulation, and the fluid volume is selected as $8{\mathrm{L}}_{\mathrm{lander}}\times 8{\mathrm{W}}_{\mathrm{lander}}\times 10{\mathrm{H}}_{\mathrm{lander}}$(length, width, height (LWH)); the dimension of the lander is 0.99 $\times$ 0.97 $\times 1.04 \mathrm{m}$. The wall conditions and wall roughness comprise an adiabatic wall and is rated at zero, respectively. The two-scale wall functions (2SWF) model [36] was used to describe boundary layers on a fine mesh; the number of cells across a boundary layer is 8, and the corresponding approach is called the ”thick boundary layer”. To select an appropriate basic mesh size and obtain high calculation accuracy with as few cells as possible, six different basic mesh sizes were used to calculate the hydrodynamic drag force during the lander’s descent process, and the incoming flow velocity is set at 0.1 m/s and 1.0 m/s. The so-called mesh irrelevance verification is shown in

Table 2. The mesh irrelevance verification.

Basic Mesh 1/(m)	Cells	Drag/(N) in 0.1 m/s	Drag/(N) in 1.0 m/s
0.17	131,459	7.18	721.98
0.15	210,683	6.02	601.49
0.13	311,779	6.98	704.71
0.11	525,830	6.97	696.15
0.09	898,422	6.83	686.12
0.07	1,974,685	6.81	681.81

¹ The basic mesh is used for the x/y/z direction.

.

As observed in Table 2, when the basic mesh is less than 0.09 m, the drag force gradually stabilizes and undergoes a small change. Considering limited computing resources, a basic mesh of 0.09 m is finally selected for the next CFD simulation.

3.2. Calculation of Damping Force Coefficients

In the CFD simulation, the descent velocity of the benthic lander is changed by adjusting the incoming flow velocity. The incoming flow velocity is set as 0.1–1.0 m/s, and the interval is 0.1 m/s. In contrast, Figure 2a, a cut plot, provides a contoured view of the velocity parameter distribution in which the range is divided by colored intervals; thus, each interval has its own color. Obviously, the symmetrical velocity distribution of the left and right sides also means that the damping force (D(υ)) of the lander is essentially symmetrical. Therefore, when only hydrodynamic forces are considered, the lander is relatively stable in its attitude. The damping force is a quadratic function of velocity only. From the simulated forces, the least mean square method is used to determine the hydrodynamic parameters based on Equation (3), and the results are shown in Figure 2b. The corresponding quadratic function is shown in Equation (6). As with the lander, the damping force function of a single floatation sphere is shown in Equation (7). For the full CFD data, we refer the readers to Supplementary Table S1: CFD data.

$$\mathrm{D}\left(\mathsf{\upsilon }\right)=686.5{\mathsf{\upsilon }}^2-0.3\mathsf{\upsilon }$$

(6)

$$\mathrm{D}\left(\mathsf{\upsilon }\right)=119.4{\mathsf{\upsilon }}^2+0.5\mathsf{\upsilon }$$

(7)

3.3. Calculation of Added Mass Coefficients

When an underwater vehicle moves in a fluid, the inertia of the fluid opposes the motion; that effect is equal to having a virtual mass added to the mass of the vehicle. In other words, when the underwater vehicle exhibits unsteady motion, it will produce added mass; we refer the reader to Equation (3). The velocities with sine functions [37] or ramp functions [38] are widely used in unsteady motions. In this study, the velocity with ramp functions was selected; we refer the reader to Equation (8), where “a” is the constant acceleration and the parameter is set as 3 groups: 0.100, 0.125 and 0.150 $\mathrm{m}/{\mathrm{s}}^2$.

$$\upsilon =\mathrm{a}\cdot t$$

(8)

In the SOLIDWORKS Flow Simulation, the simulation project selects the time-dependent physical feature. The fluid’s unsteady motion requires a certain time period (for example, in group 3, the time approx. 13 s) to reach relatively stable acceleration velocities, and the actual acceleration is different from the set value. The actual velocity development is shown in Figure 3a, where the black dashed line represents the actual acceleration (or slope) when the steady state is attained. As shown in Equation (3), the total drag force comprises two parts: the damping force calculated in Equation (6) and the difference force generated by the added mass. The results of the explained difference force are shown in Figure 3b (in group 3, the data from groups 1 and 2 are shown in Figure A3 and Figure A4, respectively). The dashed line is the boundary at which the fluid reaches a relatively stable acceleration; obviously, the difference force is relatively constant at this boundary. Finally, bringing simulation data into Equation (3), we can obtain the added mass; we refer the reader to

Table 3. Summary of added mass calculation.

Group	$\mathbf{Setting} \mathbf{Acceleration} (\mathbf{m}/{\mathbf{s}}^2)$	$\mathbf{Actual} \mathbf{Acceleration} (\mathbf{m}/{\mathbf{s}}^2)$	Boundary(s)	${\mathbf{M}}_{\mathbf{A}}$(kg) 1
1	0.100	0.143	14–18	4008.4
2	0.125	0.173	12–16	4138.2
3	0.150	0.212	12–14	4163.9
The mean value of ${\mathrm{M}}_{\mathrm{A}}$(kg) in 3 groups				4103.5

¹ Note: the ${\mathrm{M}}_{\mathrm{RB}}$ is 588.6 kg.

.

4. MATLAB Simulink

To study the influence of the lander’s configuration parameters on the variables of the descent process and based on the established dynamic model and calculated hydrodynamic coefficients, the MATLAB Simulink simulation environment was used to design the block diagram model (as we can see from Figure 4), and the corresponding file is shown in Supplementary Model S1: Simulink dynamic model. This simulation model has three input variables—the number of floatation spheres (N), the mass of the weight stack (${\mathrm{m}}_{\mathrm{stack}}$) and the bottom area of the weight stack (${\mathrm{A}}_{\mathrm{stack}}$) —and three output variables—the descent velocity (${\mathrm{V}}_{\mathrm{des}}$), descent time (${\mathrm{T}}_{\mathrm{des}}$) and subsidence depth (${\mathrm{D}}_{\mathrm{sub}}$) with the seafloor. The variable explanation is as follows, and the nomenclature is shown in

Table 4. Nomenclature.

Symbol	Description	Unit	Symbol	Description	Unit
N	number of floatation spheres	-	${\mathrm{T}}_{\mathrm{des}}$	descent time	H
${\mathrm{m}}_{\mathrm{stack}}$	mass of the weight stack	kg	${\mathrm{V}}_{\mathrm{des}}$	descent velocity	m/s
${\mathrm{A}}_{\mathrm{stack}}$	bottom area of the weight stack	${\mathrm{m}}^2$	${\mathrm{D}}_{\mathrm{sub}}$	subsidence depth	mm

.

• ${\mathrm{A}}_{\mathrm{stack}}$ influences the hydrodynamic coefficients (as we can see from Figure A5): The corresponding function relationship is calculated by using SOLIDWORKS Flow Simulation. The results are shown in Equations (9) and (10);
• ${\mathrm{T}}_{\mathrm{des}}$ takes the water depth of 5133 m during the sea trial as the background: The corresponding function is shown in Equation (11);
• ${\mathrm{D}}_{\mathrm{sub}}$ was tested by a previous study [39,40] for the same benthic lander: The corresponding function is shown in Equation (12).

  $${\mathrm{D}}_{\mathrm{quad}.}=14.29{\mathrm{A}}_{\mathrm{stack}}{}^2-23.17{\mathrm{A}}_{\mathrm{stack}}+9.92$$

  (9)

  $${\mathrm{D}}_{\mathrm{lin}.}=14.29{\mathrm{A}}_{\mathrm{stack}}{}^2-23.17{\mathrm{A}}_{\mathrm{stack}}+9.92$$

  (10)

  $${\mathrm{T}}_{\mathrm{des}}=\frac{5133}{3600\mathsf{\upsilon }}$$

  (11)

  $${\mathrm{D}}_{\mathrm{sub}}=30.7\mathsf{\upsilon }+2.65$$

  (12)

The dynamic model in Equation (3) is a nonlinear, first-order differential equation that may be integrated (in the Simulink diagram model, the Integrator block is used; we refer the readers to Figure A5) numerically to yield vehicle translational velocities when provided with the suitable initial conditions. The fixed-step size is 0.02 s by using the ode4 (Runge–Kutta) solver to solve the nonlinear differential equation.

Keeping the number of floatation spheres (N) constant at five in the descent process, the steady state motion was calculated by changing the mass of the weight stack (${\mathrm{m}}_{\mathrm{stack}}$) and the bottom area of the weight stack (${\mathrm{A}}_{\mathrm{stack}}$). The developmental rules of the subsidence depth (${\mathrm{D}}_{\mathrm{sub}}$) and descent time (${\mathrm{T}}_{\mathrm{des}})$ are shown in Figure 5. Obviously, the subsidence depth increases with ${\mathrm{m}}_{\mathrm{stack}}$; the smaller the bottom area, the greater ${\mathrm{m}}_{\mathrm{stack}}$ affects the subsidence depth. This shows that when the bottom area is small, it is easier to change the subsidence depth of the lander by adjusting ${\mathrm{m}}_{\mathrm{stack}}$. In contrast, the descent time decreases with ${\mathrm{m}}_{\mathrm{stack}}$; therefore, ${\mathrm{m}}_{\mathrm{stack}}$ shall be appropriately increased in the sea trial in pursuit of operational efficiency.

Keeping the bottom area of the weight stack (${\mathrm{A}}_{\mathrm{stack}}$) constant at 0.96 ${\mathrm{m}}^2$ in the descent process, the steady state motion was calculated by changing the mass of the weight stack (${\mathrm{m}}_{\mathrm{stack}}$) and the number of floatation spheres (N). The results of the descent time (${\mathrm{T}}_{\mathrm{des}}$) and subsidence depth (${\mathrm{D}}_{\mathrm{sub}}$) are shown in Figure 6. As the number of floatation spheres increases, ${\mathrm{T}}_{\mathrm{des}}$ increases but ${\mathrm{D}}_{\mathrm{sub}}$ decreases. When ${\mathrm{m}}_{\mathrm{stack}}$ is small, ${\mathrm{T}}_{\mathrm{des}}$ and ${\mathrm{D}}_{\mathrm{sub}}$ are more sensitive to changes in the number of floatation spheres. For example, when N is six, ${\mathrm{T}}_{\mathrm{des}}$ increased dramatically and lost its practical significance; thus, this situation was deleted from this Figure. Therefore, ${\mathrm{m}}_{\mathrm{stack}}$ should not be too small in value in order to reduce the influence of ${\mathrm{T}}_{\mathrm{des}}$ and ${\mathrm{D}}_{\mathrm{sub}}$.

5. Multi-Objective Optimization

Based on the analysis of MATLAB Simulink, we know the relationship between specific configuration parameters and descent time (${\mathrm{T}}_{\mathrm{des}}$) and subsidence depth (${\mathrm{D}}_{\mathrm{sub}}$) of the lander, but making a choice is still difficult because there are various configuration schemes. Second, in practical applications, the descent time should be as small as possible to increase the operation efficiency in sea trials, and the subsidence depth should be as small as possible to avoid retrieval failure. As we can see from Figure 5, the descent time and subsidence depth are two conflicting objectives; thus, it is a multi-objective optimization problem (MOP). In this case, the objective of optimization is to minimize the descent time and subsidence depth simultaneously, which can be formulated [41] as follows:

$$\left\{\begin{array}{c}\underset{\mathrm{X}}{\min }\mathrm{F}\left(\mathrm{X}\right)={\left({\mathrm{f}}_1\left(\mathrm{X}\right),{\mathrm{f}}_2\left(\mathrm{X}\right)\right)}^{\mathrm{T}}\\ \mathrm{s}.\mathrm{t}. {\mathrm{g}}_{\mathrm{i}}\left(\mathrm{X}\right)\le 0,\mathrm{i}=1,2,3\end{array}\right.$$

(13)

where $\mathrm{X}={\left({\mathrm{x}}_1,{\mathrm{x}}_2,{\mathrm{x}}_3\right)}^{\mathrm{T}}$ is the decision vector; decision variable ${\mathrm{x}}_1\in \left[4,5,6\right]$ is the number of floatation spheres (N), which is configured by previous sea trial experience; ${\mathrm{x}}_2\in \left(370-580\right)$ is the mass of the weight stack (${\mathrm{m}}_{\mathrm{stack}}$), which is limited by the boundary of the bottom area of the weight stack (${\mathrm{A}}_{\mathrm{stack}}$); and ${\mathrm{x}}_3\in \left(0.8-1.2\right)$ is the bottom area of the weight stack (${\mathrm{A}}_{\mathrm{stack}}$), which is limited by the number of scientific device carried by the lander. $\mathrm{F}\left(\mathrm{X}\right)$ is the objective vector, objective value ${\mathrm{f}}_1\left(\mathrm{X}\right)$ is the descent time (${\mathrm{T}}_{\mathrm{des}}$) defined in Equation (11), and ${\mathrm{f}}_2\left(\mathrm{X}\right)$ is the subsidence depth (${\mathrm{D}}_{\mathrm{sub}}$) defined in Equation (12). ${\mathrm{g}}_{\mathrm{i}}\left(\mathrm{X}\right)$ denotes inequality constraints, where ${\mathrm{g}}_1\left(\mathrm{X}\right)$ ensures that the total net gravity is positive, ${\mathrm{g}}_2\left(\mathrm{X}\right)$ ensures that the descent time is less than 3 h, and ${\mathrm{g}}_3\left(\mathrm{X}\right)$ ensures that the subsidence depth is less than 45 mm. The details are described as follows, and the corresponding MATLAB code is shown in Algorithms B1.

$${\mathrm{g}}_1\left(\mathrm{X}\right)=26\mathrm{N}-\frac{1869.7}{2869.7}{\mathrm{x}}_2+85\le 0$$

(14)

$${\mathrm{g}}_2\left(\mathrm{X}\right)={\mathrm{f}}_1\left(\mathrm{X}\right)-3\le 0$$

(15)

$${\mathrm{g}}_3\left(\mathrm{X}\right)={\mathrm{f}}_2\left(\mathrm{X}\right)-45\le 0$$

(16)

Based on the above MOP, the platform for evolutionary multi-objective optimization (PlatEMO) [42] was used to solve this problem, and the corresponding multi-objective evolutionary algorithms (MOEAs) comprise NSGA-II; [43]. The population size is set to 100, and the evaluation number is 10,000. The Pareto-optimal front (objective value) is obtained as shown in Figure 7a, and the corresponding Pareto-optimal set (decision variable) is shown in Figure 7b. The black, red and blue points distinguish the different numbers of floatation spheres. Finally, parts of the configuration from the Pareto-optimal front and set are shown in

Table 5. Summary of optimization results.

Decision Variables			Objective Values		Decision Variables			Objective Values
N	$m_{stack}$ (kg)	$A_{stack}$$\left({\mathrm{m}}^2\right)$	$T_{des}$ (H)	$D_{sub}$ (mm)	N	$m_{stack}$ (kg)	$A_{stack}$$\left(m^2\right)$	$T_{des}$ (H)	$D_{sub}$ (mm)
4	557.7	0.89	1.03	45.0	4	511.8	0.92	1.18	39.8
4	510.0	0.83	1.10	42.3	4	508.1	0.87	1.13	41.3
5	392.0	0.96	1.76	27.5	5	441.9	1.20	2.63	19.3
5	468.3	1.19	2.35	21.3	5	444.5	1.10	2.24	22.2
6	410.4	0.89	3.00	17.2	6	553.6	0.92	1.44	33.1
6	436.5	0.88	2.32	21.5	6	414.0	0.89	2.88	17.9

(Supplementary Table S2: multi-objective optimization), which provides an accurate configuration reference depending on the specific ${\mathrm{T}}_{\mathrm{des}}$ and ${\mathrm{D}}_{\mathrm{sub}}$.

The goal of MOEAs is to provide the last population with good convergence and diversity when the evaluation ends; convergence and diversity are also two important performance indicators for evaluating the optimization results. Diversity is usually measured by the diversity metric (DM) [44], and larger values are better. Convergence is usually measured by the inverted generational distance (IGD) [45], and smaller values are better. Finally, the results of the DM and IGD indicators in this study are shown in Figure 8. Obviously, the IGD tends to be stable when the number of evaluations is greater than 500; the DM has some oscillation, but the average value is approximately 0.7, both of which meet the requirements of engineering practice.

6. Sea Trials

The benthic lander sea trial was conducted in the South China Sea in July 2018 [39,40], and the newest sea trial was conducted in the western Pacific Ocean (location: 134.84°E, 16.96°N) in January 2021 (Figure 9). The corresponding video is shown in Supplementary Video S1: sea trial. The newest sea trial will be discussed in this study, and the basic configuration is N = 5, ${\mathrm{m}}_{\mathrm{stack}}=392 \mathrm{kg}$, and ${\mathrm{A}}_{\mathrm{stack}}=0.96{ \mathrm{m}}^2$, one of the optimization results from Table 5. After the lifting device is released from the benthic lander, gravity provided by the weight stack drives the lander to freely land on the bottom of the seafloor. The original document that shows the range function of the acoustic release that obtains the relationship between the descent depth and descent time (

Table 6. The depth information comes from acoustic release.

Descent Time(min)	Descent Depth(m)	Corresponding Velocity(m/s)	Descent Time(min)	Descent Depth(m)	Corresponding Velocity(m/s)
0	0	-	28	1536	0.89
2	88	0.73	41	2083	0.70
4	343	2.13	58	2908	0.81
6	434	0.76	70	3396	0.68
11	694	0.87	95	5071	1.12
17	951	0.71	105	5133	0.10

and Figure 10a) can be found in Supplementary Document S1: Lander deployment log sheet, and the following details are noteworthy:

• The depth information reflected by the range data of the acoustic release is relatively accurate and reliable. First, the total mass of the lander system (603 kg) is large, and the sea current has little influence on it. Second, the successful retrieval of the lander at the original location in the later sea trials proved that the horizontal drifts during the descent process can be ignored;
• Removal of instability during the initial launch and final landing, and the steady descent (6–70 min) velocity fluctuates between 0.68 and 0.89 m/s, which meets the expected results.

As observed in Figure 10b, the Simulink simulation velocity is basically consistent with the mean descent velocity in the sea trial; the correctness of the Simulink simulation results and CFD hydrodynamic calculation was verified; the mean descent velocity was calculated in the steady descent (6–70 min) process. As the added mass is considered in Section 2.1, the acceleration process of the lander is relatively long; at approx. 50 s, its velocity reaches the maximum and remains relatively stable. The descent time (${\mathrm{T}}_{\mathrm{des}}$) calculated based on the average descent velocity of the sea trial is 1.85 h. Table 5 shows that, under the sea trial configuration, the descent time obtained by using multi-objective optimizations is 1.76 h, and the relative error between the two is 5.1%. Thus, the correctness of multi-objective optimizations was verified.

7. Discussion

In this study, based on the analysis of MATLAB Simulink and CFD methods, the multi-objective optimization algorithm was used to achieve a rational configuration scheme for a benthic lander that realizes reasonable control of descent velocity. This scheme makes up the shortcomings of previous studies, for example, using crossed wings (Mortensen et al. [31]) or integrating a hydrofoil (Jun [32] and Gang et al. [33]) into the benthic lander. On the other hand, the main reasons for the higher Simulink simulation velocity (see Figure 10b) can be concluded as follows. First, the hydrodynamic drag coefficients of the lander and the floatation spheres are considered, but the coupling effects between the two are neglected. Second, the CFD damping force calculation is based on the static steady descent velocity, but the actual descent velocity has some oscillations (Figure 10a). Finally, as described in Section 6, the horizontal drifts of the lander during the descent process can be ignored, but the drifts exist, which leads to the actual descent depth being smaller than the range data from the acoustic release. However, we can conclude that the overall trend is basically correct.

In specific applications, we need to decide the descent time (${\mathrm{T}}_{\mathrm{des}}$) tolerated by the sea trials and the subsidence depth (${\mathrm{D}}_{\mathrm{sub}}$) tolerated by the lander according to Figure 7a, and the two values were used to choose a rational configuration from Table 5 (the full table see Supplementary Table S2: multi-objective optimization). Significantly, the proper number of floatation spheres (N) should be carefully selected because N has a large influence on ${\mathrm{T}}_{\mathrm{des}}$ and ${\mathrm{D}}_{\mathrm{sub}}$ (Figure 6). The population size and number of evaluations should be appropriately increased when using multi-objective optimization, which can improve the choice of configuration scheme and enhance the convergence/diversity of the population, respectively.

The results were calculated by a multi-objective optimization algorithm, and only configurations with five floatation spheres were verified in the sea trial. In future work, the number of floatation spheres at four and six still require further verification, which can comprehensively verify the correctness of the proposed scheme. The applicability of CFD simulation also needs to be tested more directly, especially since the calculation of the added mass coefficients usually has a relatively large error. In future studies, a spherical example should be introduced to test the relevant settings of SOLIDWORKS Flow Simulation in order to reduce simulation errors.