Study on the Dynamic Responses of a Large Caisson during Wet-Towing Transportation

Large caissons are extensively applied as deep-water foundations in marine engineering. In fact, caissons are generally prefabricated and transported to project site by wet towing. Motion responses of large caissons and those occurring during the towing process were investigated, and CO2 emissions under various conditions were calculated. These are all considered to ensure towing safety and environmental protection. The caisson resistance coefficient was simulated via Ansys Fluent software. The effects of towrope length, towing speed, and drift depth on the motion responses of caissons under the combined action of wind and wave were evaluated via Ansys AQWA software. Maximum heave value was dominantly affected by rope length and draft depth, and its fluctuation was highly influenced by towing speed and draft depth. However, all of the above mentioned factors had insignificant influences on pitch response. When towing existed, rope tension was rapidly increased from zero to a constant value that depended on towing speed and drift depth. However, the speed of achieving this stable phase depended on the length of the towrope.


Introduction
To response to the requirements of today's globalized world, oceanographic engineering has rapidly progressed [1]. In addition, caissons have been extensively applied in offshore structures such as wind turbines, cross-sea bridges, and gas and oil platforms. Today, high quality onshore caissons are first constructed, wet-towed with a small tugboat to project site, and installed on the seabed via a water-displacement self-sinking method. Towing operations have been performed for several years as an inexpensive transportation method for a great quantity of commodities [2]. Kim et al. compared the costs of wet and dry towing methods for a semisubmersible platform. The obtained results showed that wet towing was the cheapest option for their project [3,4]. However, in some environments, such as port, open sea, coastal waters, and internal waters, water towing operations may be the only choice. Each of the abovementioned environments have a variety of potential hazards [5]. In addition, regardless of planning feasibility, accidents can take place in towing operations, even in mild weather, inflicting enormous casualties and economical losses [6][7][8]. Hence, accurate studies have to be conducted to address the requirement of wet towing operations and ensure the practicability of the developed integrated transportation method.
Towing operations, which include the transportation of self-floating objects using one or more towing tugs, are very common in ocean engineering projects [9]. Towing dynamics are very similar to those of moored objects in currents [5]. In recent years, several authors, some of whom have applied traditional theoretical methods, have investigated this topic. Bernitsas and Kekridis [10] proposed a 3-Degrees of Freedom (3-DOF) fully analytical model for the investigation of the towing stabilities of some vessels towed by

Towing System
In the current study, the towing system consisted of a tug boat, towing rope, and caisson and was simulated using Ansys software ( Figure 2). Both global and local coordinate systems were applied. The tug boat and caisson were modeled in a local coordinate system, and the towing operation was simulated by Ansys AQWA. The towing rope was made of a mix of steel wire with a large diameter, produced by JIAN FENG SLING CO., LTD (Guangzhou, China) [19], and its cross-section is shown in Figure 3. It should be noted that, because only one towrope was used in the towing process, at high towing speeds, multiple tugs might need to be tied to the end of the towrope to provide enough power. However, when studying the dynamic responses of the caisson, we simplified the model and assumed that it was towed by one tugboat.

Towing System
In the current study, the towing system consisted of a tug boat, towing rope, and caisson and was simulated using Ansys software ( Figure 2). Both global and local coordinate systems were applied. The tug boat and caisson were modeled in a local coordinate system, and the towing operation was simulated by Ansys AQWA. The towing rope was made of a mix of steel wire with a large diameter, produced by JIAN FENG SLING CO., LTD (Guangzhou, China) [19], and its cross-section is shown in Figure 3. It should be noted that, because only one towrope was used in the towing process, at high towing speeds, multiple tugs might need to be tied to the end of the towrope to provide enough power. However, when studying the dynamic responses of the caisson, we simplified the model and assumed that it was towed by one tugboat.

Intact Stability Analysis
Metacentric height (GM) is an essential key parameter for the measurement of floating structure initial static stability, and this was calculated as the distance from the gravity center to the metacenter. Higher metacentric values indicated favorable initial stabilities against overturning. Sometimes the height of the metacentric also reflects the rolling model frequency of the structure [20]. In the current work, metacentric height was obtained by: where BM is the metacentric radius obtained as distance from buoyancy center to metacenter, ∇ is displacement volume, I is inertia moment from floating structure waterplane to axis, B Z is buoyancy center, and G Z is gravity center.

Intact Stability Analysis
Metacentric height (GM) is an essential key parameter for the measurement of floating structure initial static stability, and this was calculated as the distance from the gravity center to the metacenter. Higher metacentric values indicated favorable initial stabilities against overturning. Sometimes the height of the metacentric also reflects the rolling model frequency of the structure [20]. In the current work, metacentric height was obtained by: where BM is the metacentric radius obtained as distance from buoyancy center to metacenter, ∇ is displacement volume, I is inertia moment from floating structure waterplane to axis, B Z is buoyancy center, and G Z is gravity center.

Intact Stability Analysis
Metacentric height (GM) is an essential key parameter for the measurement of floating structure initial static stability, and this was calculated as the distance from the gravity center to the metacenter. Higher metacentric values indicated favorable initial stabilities against overturning. Sometimes the height of the metacentric also reflects the rolling model frequency of the structure [20]. In the current work, metacentric height was obtained by: where BM is the metacentric radius obtained as distance from buoyancy center to metacenter, ∇ is displacement volume, I is inertia moment from floating structure waterplane to axis, Z B is buoyancy center, and Z G is gravity center.

Hydrodynamics Analysis
Wind, currents, and waves are the three key factors affecting the stability of the floating structure during wet towing. Hence, the responses of the structures in different environments is a great problem that has been solved in this work. The main goal of this study was to provide proper insight into practical applications such as offshore structure design, installation plans, and transportation formulation. The time histories of caisson motion with six degrees of freedom (6-DOF) were stated in the reference system by adopting the water plane as the origin. Based on 3D potential theory, hydrodynamic moments and forces on the caisson were obtained.
The dynamic equation for a floating body under various complex loadings is stated as: ..
where M is floating body mass matrix, ∆M is floating body added mass matrix, B vis is viscous damping matrix, B rad is radiation damping matrix, K stillwater is hydrostatic stiffness matrix, K mooring is mooring system stiffness matrix, F 1 is first order frequency load matrix, F 2Low is second order low frequency load matrix, F 2High is second order high frequency load matrix, F wind is wind load matrix, F current is current load matrix, and F others is the remaining load matrix. Floating body mass matrix is stated as: where (x G , y G , z G ) is gravity center position and I ij is inertia mass. Hydrostatic stiffness matrix is stated as: where (X B , Y B , Z B ) is buoyancy center position, S is water plane area, and S i and S ij are the first and second order moments of the water plane area. In the current work, ∆M, B rad , F 1 , F 2Low , and F 2High are obtained by AQWA-Line software and B vis is determined by Morison theory. Wind and current drag were obtained in a similar way. Environment load coefficients are defined as: where F is drag force, A is the area of body incident to flow, and v is velocity relative to wind or current positions. Therefore, force is determined as: In the current work, load coefficients have been determined using Ansys Fluent software. Figure 4 shows the computational domain obtained based on the findings of Lee et al. [21] and Liu et al. [22]. Fore body was 1.5 L away from the velocity inlet and aft body was 2.5 L away from the pressure outlet. In addition, lateral distance between wall and body was 2 L. First-order temporal and second-order convection schemes were applied for temporal and momentum equations, respectively, to be applied to perform spatial discretization. Jin et al. [23] and He et al. [24] applied the k − ω Shear Stress Transfer (SST) turbulence model to analyze the Reynolds averaged Navier-Stokes (RANS) equation and obtained satisfactory results. In the current work, the k − ω SST turbulence model was also applied.
In the current work, load coefficients have been determined using Ansys Fluent software. Figure 4 shows the computational domain obtained based on the findings of Lee et al. [21] and Liu et al. [22]. Fore body was 1.5 L away from the velocity inlet and aft body was 2.5 L away from the pressure outlet. In addition, lateral distance between wall and body was 2 L. First-order temporal and second-order convection schemes were applied for temporal and momentum equations, respectively, to be applied to perform spatial discretization. Jin et al. [23] and He et al. [24] applied the − Shear Stress Transfer (SST) turbulence model to analyze the Reynolds averaged Navier-Stokes (RANS) equation and obtained satisfactory results. In the current work, the − SST turbulence model was also applied. For the verification of the accuracy of the developed method, square cylinder load coefficients applied in the Tang's [25] experimental research were obtained, and the obtained results are summarized in Table 2. The result errors of the two methods were lower than 5%, which meant that the developed numerical method could effectively calculate load coefficients. Environment load coefficients of the investigated caisson are shown in Figure 5.  For the verification of the accuracy of the developed method, square cylinder load coefficients applied in the Tang's [25] experimental research were obtained, and the obtained results are summarized in Table 2. The result errors of the two methods were lower than 5%, which meant that the developed numerical method could effectively calculate load coefficients. Environment load coefficients of the investigated caisson are shown in Figure 5.
In the current work, load coefficients have been determined using Ansys Fluent software. Figure 4 shows the computational domain obtained based on the findings of Lee et al. [21] and Liu et al. [22]. Fore body was 1.5 L away from the velocity inlet and aft body was 2.5 L away from the pressure outlet. In addition, lateral distance between wall and body was 2 L. First-order temporal and second-order convection schemes were applied for temporal and momentum equations, respectively, to be applied to perform spatial discretization. Jin et al. [23] and He et al. [24] applied the − Shear Stress Transfer (SST) turbulence model to analyze the Reynolds averaged Navier-Stokes (RANS) equation and obtained satisfactory results. In the current work, the − SST turbulence model was also applied. For the verification of the accuracy of the developed method, square cylinder load coefficients applied in the Tang's [25] experimental research were obtained, and the obtained results are summarized in Table 2. The result errors of the two methods were lower than 5%, which meant that the developed numerical method could effectively calculate load coefficients. Environment load coefficients of the investigated caisson are shown in Figure 5. Angle of attack (deg)   The Joint North Sea Wave Project (JONSWAP) spectrum was applied as the wave spectrum. Empirical parameters γ and α, as well as peak frequency, were also applied. Spectral ordinate at any frequency is stated as: where ω P is peak frequency, γ is peak enhancement factor, and α is a constant value that depends on wave spectrum peak frequency and wind speed, which can be determined as: where ω > ω P Starting and finishing frequencies are also expressed as: where F(γ) is a weighting function, and weighting function values against γ ∈ [1.0, 20.0] are available in the AQWA Theory Manual.

Mesh Convergence Analyses
The analysis of mesh density had to be performed when calculating mean wave drift moments and forces on the structure using both Pinkster (near-field) and Maruo-Newman (far-field) methods to obtained optimum accuracy and efficiency. Mesh size had to be smaller than both one seventh of the wave length size and the minimum dimensional feature of the caisson. A full-sized caisson model was developed with mesh sizes of 0.8, 1.0, 1.2, and 1.5 m. Mean wave drift moments and forces obtained from the two methods were compared. All parameters except element size were the same when using both models. The obtained results are shown in Figure 6.
As can be seen from Figure 6, the mean wave drift moments and forces obtained by the various methods were almost the same at time periods longer than 20 s. However, at periods of less than 20 s, the findings of the two methods were different. Finally, at a mesh size of 1 m, the results obtained from the two methods were similar. Based on the abovementioned findings, a 1 m mesh size was adopted in this analysis.  Period (second) ( ) Ⅲ Pinkster Method Maruo-Newman method Yaw (N/m)

Estimation of CO2 Emissions
Currently, an increasing number of engineering projects are required to take proper measures regarding air emissions duo to global warming. In this study, we also estimated CO2 emissions during towing. It was assumed from experience that a tugboat could provide a pulling force of 199.92 kN/1000 kW. We then calculated the number of required tugboats based on rope tension. Emissions were estimated according to the following formula [26]: where MCR is the maximum continuous of the combustion engine in use, kW; LF is the engine load factor during specific activity; A is activity time, h; and EF is emissions factor, kg kWh .

Stability Analysis Results
Intact stability parameters of the caisson were calculated using theoretical and numerical simulation methods, and the obtained results are shown in Table 3. From the data given in the table, it was found that the difference of GM obtained by the two methods was less than 5%. Additionally, the GM value was decreased by increasing drift, which agreed with the findings of Rawson [20], reported in his book. All above conclusions proved the accuracy and feasibility of the developed numerical method.
The caisson frequency domain was also analyzed. Figure 7 shows caisson's response amplitude operators (RAOs) for a wave direction of 0 degrees. First-order wave frequency was about 0.5-0.8 rad/s. As the drift increases, the peak of RAOs increases accordingly.

Estimation of CO 2 Emissions
Currently, an increasing number of engineering projects are required to take proper measures regarding air emissions duo to global warming. In this study, we also estimated CO 2 emissions during towing. It was assumed from experience that a tugboat could provide a pulling force of 199.92 kN/1000 kW. We then calculated the number of required tugboats based on rope tension. Emissions were estimated according to the following formula [26]: where MCR is the maximum continuous of the combustion engine in use, kW; LF is the engine load factor during specific activity; A is activity time, h; and EF is emissions factor, kg /kWh.

Stability Analysis Results
Intact stability parameters of the caisson were calculated using theoretical and numerical simulation methods, and the obtained results are shown in Table 3. From the data given in the table, it was found that the difference of GM obtained by the two methods was less than 5%. Additionally, the GM value was decreased by increasing drift, which agreed with the findings of Rawson [20], reported in his book. All above conclusions proved the accuracy and feasibility of the developed numerical method.
The caisson frequency domain was also analyzed. Figure 7 shows caisson's response amplitude operators (RAOs) for a wave direction of 0 degrees. First-order wave frequency was about 0.5-0.8 rad/s. As the drift increases, the peak of RAOs increases accordingly. The above findings lay a solid foundation for follow-up time domain analysis.

Results of Dynamic Response
Towing float hydrodynamic responses rely on many factors such as sea load tions and drafts as well as towing system parameters such as tow speed, towrope etc., and a safe towing plan can be easily developed by analyzing the influences of v factors on the dynamic responses. Random sea conditions according to a JONSWA trum with high corresponding frequency and wave height were applied in this wo the basis of the real sea conditions of a towing project, wave height and frequency of 1 m and 1.03 rad/s, respectively, were adopted in this work. Wind was assumed steady flow with a speed of 6 m/s. where Max is maximum value in the whole process, Mean is arithmetic mean val STD is standard deviation.
As seen in Table 4, STD and normalized max. error were generally less than 1 different simulation times, and mean value barely changed. From Table 5, it is show caisson pitch for various simulation times also presented the same variation trend. F 8 and 9 demonstrate motion responses for various time durations. Figure 8 compare response in different durations; it was found that pitch amplitude varied significa the beginning of towing operation (before 600 s). After 1 h, however, pitch amplitu ied periodically. It can be clearly seen from Figure 9 that heave amplitude changed ularly at early towing stages (before 600 s) and, after 1 h, heave change tended to ular.

Results of Dynamic Response
Towing float hydrodynamic responses rely on many factors such as sea load conditions and drafts as well as towing system parameters such as tow speed, towrope length, etc., and a safe towing plan can be easily developed by analyzing the influences of various factors on the dynamic responses. Random sea conditions according to a JONSWAP spectrum with high corresponding frequency and wave height were applied in this work. On the basis of the real sea conditions of a towing project, wave height and frequency values of 1 m and 1.03 rad/s, respectively, were adopted in this work. Wind was assumed to be a steady flow with a speed of 6 m/s.
where Max is maximum value in the whole process, Mean is arithmetic mean value, and STD is standard deviation. As seen in Table 4, STD and normalized max. error were generally less than 10% for different simulation times, and mean value barely changed. From Table 5, it is shown that caisson pitch for various simulation times also presented the same variation trend. Figures 8 and 9 demonstrate motion responses for various time durations. Figure 8 compares pitch response in different durations; it was found that pitch amplitude varied significantly at the beginning of towing operation (before 600 s). After 1 h, however, pitch amplitude varied periodically. It can be clearly seen from Figure 9 that heave amplitude changed irregularly at early towing stages (before 600 s) and, after 1 h, heave change tended to be regular.

Drift Depth and Towing Speed
Because of the large size of the caisson, its wet surface and water plane could change significantly under violent heave and pitch motions. Four towing speeds and three drift conditions were studied to compare cable tension, heave, and pitch. It is noteworthy thet the above three drift depths were greater than those for self-weight balancing and, therefore, ballast weight was applied. It should also be kept in mind that ballast weight was assumed to be placed at the caisson gravity center.
Dynamic response analyses for various drift depths and towing speeds were performed for a towrope length of 100 m. Figure 10 shows the normalized maximums of heave motion and pitch angle. Figure 11 presents bar charts of the STD of caisson responses. Figure 12 illustrates towrope tension variation during towing.

Drift Depth and Towing Speed
Because of the large size of the caisson, its wet surface and water plane could change significantly under violent heave and pitch motions. Four towing speeds and three drift conditions were studied to compare cable tension, heave, and pitch. It is noteworthy thet the above three drift depths were greater than those for self-weight balancing and, therefore, ballast weight was applied. It should also be kept in mind that ballast weight was assumed to be placed at the caisson gravity center.
Dynamic response analyses for various drift depths and towing speeds were performed for a towrope length of 100 m. Figure 10 shows the normalized maximums of heave motion and pitch angle. Figure 11 presents bar charts of the STD of caisson responses. Figure 12 illustrates towrope tension variation during towing.

Drift Depth and Towing Speed
Because of the large size of the caisson, its wet surface and water plane could change significantly under violent heave and pitch motions. Four towing speeds and three drift conditions were studied to compare cable tension, heave, and pitch. It is noteworthy thet the above three drift depths were greater than those for self-weight balancing and, therefore, ballast weight was applied. It should also be kept in mind that ballast weight was assumed to be placed at the caisson gravity center.
Dynamic response analyses for various drift depths and towing speeds were performed for a towrope length of 100 m. Figure 10 shows the normalized maximums of heave motion and pitch angle. Figure 11 presents bar charts of the STD of caisson responses. Figure 12 illustrates towrope tension variation during towing.    As can be seen from Figure 10a, at a constant towing speed, the normalized maximums of pitch angle barely changed (around 2.6-3.2) when increasing drift depth, especially at a towing speed of 6 kn. Additionally, at constant drift depth, normalized maximums of pitch angle slightly fluctuated when increasing towing speed. Figure 10b reveals that, at drift depths of less than 13.5 m, normalized maximums of heave motion barely changed when increasing towing speed. However, at a drift depth of 13.5 m, normalized maximum of heave was less than the former one and increased with towing speed.
As can be seen in Figure 11a, towing speed and drift depth had insignificant effects on pitch angle STD because the caisson water plane was constant no matter the drift depth changes. So, large wet surface areas could have small effects on pitch motion [18]. As can be seen in the bars presented in Figure 11b, heave motion fluctuations by increasing towing speed were significant. Additionally, at a drift depth of 12.5 m and a constant towing speed, heave motion STD was smaller than those of the other two conditions. Towrope tension is another key index for the measurement of towing safety. Towrope tension variations at various towing speeds and drift depths were evaluated, and the obtained results are shown in Figure 12. Towrope tension tended to remain constant after towing for about 1 h. By increasing towing speed, the time required for stabilization was increased. According to the towrope tension [28] calculation equation, flow surface areas and relative velocity greatly influenced the towrope. Therefore, an increase of towing speed increased the towrope tension.  As can be seen from Figure 10a, at a constant towing speed, the normalized maximums of pitch angle barely changed (around 2.6-3.2) when increasing drift depth, especially at a towing speed of 6 kn. Additionally, at constant drift depth, normalized maximums of pitch angle slightly fluctuated when increasing towing speed. Figure 10b reveals that, at drift depths of less than 13.5 m, normalized maximums of heave motion barely changed when increasing towing speed. However, at a drift depth of 13.5 m, normalized maximum of heave was less than the former one and increased with towing speed.
As can be seen in Figure 11a, towing speed and drift depth had insignificant effects on pitch angle STD because the caisson water plane was constant no matter the drift depth changes. So, large wet surface areas could have small effects on pitch motion [18]. As can be seen in the bars presented in Figure 11b, heave motion fluctuations by increasing towing speed were significant. Additionally, at a drift depth of 12.5 m and a constant towing speed, heave motion STD was smaller than those of the other two conditions. Towrope tension is another key index for the measurement of towing safety. Towrope tension variations at various towing speeds and drift depths were evaluated, and the obtained results are shown in Figure 12. Towrope tension tended to remain constant after towing for about 1 h. By increasing towing speed, the time required for stabilization was increased. According to the towrope tension [28] calculation equation, flow surface areas and relative velocity greatly influenced the towrope. Therefore, an increase of towing speed increased the towrope tension. As can be seen from Figure 10a, at a constant towing speed, the normalized maximums of pitch angle barely changed (around 2.6-3.2) when increasing drift depth, especially at a towing speed of 6 kn. Additionally, at constant drift depth, normalized maximums of pitch angle slightly fluctuated when increasing towing speed. Figure 10b reveals that, at drift depths of less than 13.5 m, normalized maximums of heave motion barely changed when increasing towing speed. However, at a drift depth of 13.5 m, normalized maximum of heave was less than the former one and increased with towing speed.
As can be seen in Figure 11a, towing speed and drift depth had insignificant effects on pitch angle STD because the caisson water plane was constant no matter the drift depth changes. So, large wet surface areas could have small effects on pitch motion [18]. As can be seen in the bars presented in Figure 11b, heave motion fluctuations by increasing towing speed were significant. Additionally, at a drift depth of 12.5 m and a constant towing speed, heave motion STD was smaller than those of the other two conditions. Towrope tension is another key index for the measurement of towing safety. Towrope tension variations at various towing speeds and drift depths were evaluated, and the obtained results are shown in Figure 12. Towrope tension tended to remain constant after towing for about 1 h. By increasing towing speed, the time required for stabilization was increased. According to the towrope tension [28] calculation equation, flow surface areas and relative velocity greatly influenced the towrope. Therefore, an increase of towing speed increased the towrope tension.

Drift Depth and Towrope Length
The towrope should have adequate length to guarantee that the caisson would not collide with the tugboat during the whole towing process. Taking account of the length of the caisson and the tugboat, as well as the towing speed, the towrope length adopted in this research is 100 m, 150 m, and 200 m. Dynamic responses were analyzed for three towrope lengths and two drift depths at a constant towing speed of 3 knots.
Normalized maximums for various depths and towrope lengths in pitch and heave are shown in Figure 13. It could be easily concluded that variations in maximum pitch angle (around 3.2-3.6) were not greatly affected by towrope length by comparing the normalized maximum of the pitch angle with different towrope lengths at a given drift depth. However, the heave motion peak was influenced by towrope length. As seen from Figure 14a, pitch angle STD was barely changed with towrope length, indicating that the pitch angle fluctuation was hardly influenced by towrope length. Figure 14b reveals that towrope length had unobvious effects on heave response fluctuation. Therefore, it can be easily concluded that towrope length had a limited effect on the normalized maximum values and STDs of heave and pitch.

Drift Depth and Towrope Length
The towrope should have adequate length to guarantee that the caisson would not collide with the tugboat during the whole towing process. Taking account of the length of the caisson and the tugboat, as well as the towing speed, the towrope length adopted in this research is 100 m, 150 m, and 200 m. Dynamic responses were analyzed for three towrope lengths and two drift depths at a constant towing speed of 3 knots.
Normalized maximums for various depths and towrope lengths in pitch and heave are shown in Figure 13. It could be easily concluded that variations in maximum pitch angle (around 3.2-3.6) were not greatly affected by towrope length by comparing the normalized maximum of the pitch angle with different towrope lengths at a given drift depth. However, the heave motion peak was influenced by towrope length. As seen from Figure  14a, pitch angle STD was barely changed with towrope length, indicating that the pitch angle fluctuation was hardly influenced by towrope length. Figure 14b reveals that towrope length had unobvious effects on heave response fluctuation. Therefore, it can be easily concluded that towrope length had a limited effect on the normalized maximum values and STDs of heave and pitch.   Figure 15 illustrates the towrope tension time domain curves for various towrope lengths and drift depths. Towrope tension was increased from 0 to a maximum value, then tended to remain constant. This constant value was similar at given drift depths and towing speeds, regardless of towrope length, which totally complied with the proposed

Drift Depth and Towrope Length
The towrope should have adequate length to guarantee that the caisson would not collide with the tugboat during the whole towing process. Taking account of the length of the caisson and the tugboat, as well as the towing speed, the towrope length adopted in this research is 100 m, 150 m, and 200 m. Dynamic responses were analyzed for three towrope lengths and two drift depths at a constant towing speed of 3 knots.
Normalized maximums for various depths and towrope lengths in pitch and heave are shown in Figure 13. It could be easily concluded that variations in maximum pitch angle (around 3.2-3.6) were not greatly affected by towrope length by comparing the normalized maximum of the pitch angle with different towrope lengths at a given drift depth. However, the heave motion peak was influenced by towrope length. As seen from Figure  14a, pitch angle STD was barely changed with towrope length, indicating that the pitch angle fluctuation was hardly influenced by towrope length. Figure 14b reveals that towrope length had unobvious effects on heave response fluctuation. Therefore, it can be easily concluded that towrope length had a limited effect on the normalized maximum values and STDs of heave and pitch.   Figure 15 illustrates the towrope tension time domain curves for various towrope lengths and drift depths. Towrope tension was increased from 0 to a maximum value, then tended to remain constant. This constant value was similar at given drift depths and towing speeds, regardless of towrope length, which totally complied with the proposed Figure 14. STD of dynamic responses at different drift depths and towrope lengths. Figure 15 illustrates the towrope tension time domain curves for various towrope lengths and drift depths. Towrope tension was increased from 0 to a maximum value, then tended to remain constant. This constant value was similar at given drift depths and towing speeds, regardless of towrope length, which totally complied with the proposed calculation equation (CCS, 2011). However, by increasing the length of the towrope, tension was increased faster at the beginning of towing (less than 1 h).

Extmated of CO2 Emissions
Currently, the power of general tugboats is in the range of 4000-5000 kW. In the current work, we assumed that the tugboat worked with a constant power of 4500 kW and towing distance was 30 nautical miles. Therefore, the actual number of tugboats required for various working conditions was calculated based on towrope tension. Additionally, CO2 emissions calculated by the equation are shown in Table 6. It can be seen from the data provided in Table 6 that the number of tugboats required for the mentioned working conditions was 1 to 3. Due to relatively small differences in draft changes for several different working conditions, draft did not affect the number of required tugboats. Comparing CO2 emissions for the same draft at different speeds showed that, when speed was 3 kn, only one tugboat could afford the project. Although towing time was 10 h, CO2 emissions were only 2766.4 kg. At the speed of 4 kn, however, one less tug was used than at a speed of 6kn, although it took 2.5 h longer to tow. Surprisingly, CO2 emissions of both cases were equal to 34,149.6 kg. At the speed of 5 kn, CO2 emissions were the largest and equal to 40,979.52 kg. In addition, from an energy point of view, by the increase of CO2 emissions, more energy was actually required for the above conditions. That is to say, at the speed of 5 kn, the highest amount of energy was consumed.

Extmated of CO 2 Emissions
Currently, the power of general tugboats is in the range of 4000-5000 kW. In the current work, we assumed that the tugboat worked with a constant power of 4500 kW and towing distance was 30 nautical miles. Therefore, the actual number of tugboats required for various working conditions was calculated based on towrope tension. Additionally, CO 2 emissions calculated by the equation are shown in Table 6. It can be seen from the data provided in Table 6 that the number of tugboats required for the mentioned working conditions was 1 to 3. Due to relatively small differences in draft changes for several different working conditions, draft did not affect the number of required tugboats. Comparing CO 2 emissions for the same draft at different speeds showed that, when speed was 3 kn, only one tugboat could afford the project. Although towing time was 10 h, CO 2 emissions were only 2766.4 kg. At the speed of 4 kn, however, one less tug was used than at a speed of 6kn, although it took 2.5 h longer to tow. Surprisingly, CO 2 emissions of both cases were equal to 34,149.6 kg. At the speed of 5 kn, CO 2 emissions were the largest and equal to 40,979.52 kg. In addition, from an energy point of view, by the increase of CO 2 emissions, more energy was actually required for the above conditions. That is to say, at the speed of 5 kn, the highest amount of energy was consumed.

Conclusions
The towing of caissons with bulkheads inside was studied in the current work. The towing system consisted of the caisson, a towrope, and tugboats. Wet towing was designed as an economical and wieldy transportation method. A numerical evaluation was performed, which provided a new method to assess the feasibility of the developed method, and the following conclusions were drawn.
Drifts and towing speeds had little effect on pitch response regardless of normalized maximums and their fluctuations, while deeper drifts result in higher heave maximums when the towing speed is given. Drifts also significantly affected heave fluctuations. At drift depth of 12.5 m, STD was smaller than those of the other two cases. At a certain drift, towing speed more significantly affected pitch fluctuations compared to heave fluctuations. Once a stable state was obtained for the towing process, towrope tensions were determined by towing speeds and drifts. Towrope length was found to have an insignificant effect on heave and pitch responses. However, it could affect the speed of achieving a stable state for towrope tension. Meanwhile, CO 2 emissions were mainly affected by towing speed. Above all, the optimum towing plan for the project was a drift depth of 12.5 m, a towing speed of 3 kn, and a towrope length of 200 m. Data Availability Statement: Data available on request due to restrictions e.g., privacy or ethical. The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the Chinese law.

Conflicts of Interest:
We declare no conflict of interest.