Development of a Numerical Ice Tank Based on DEM and Physical Model Testing: Methods, Validations and Applications

Abstract

The determination of ice loads on polar vessels and offshore structures is important for ice-resistant design, safe operation, and management of structural integrity in ice-infested waters. Physical model testing carried out in an ice tank/basin is usually an important technical approach for evaluating the ice loads. However, the high cost and time consumption make it difficult to perform multiple repetitions or numerous trials. Recently, the rapid development of high-performance computation techniques provides a usable alternative where the numerical methods represented by the discrete element method (DEM) have made remarkable contributions to the ice load predictions. Based on DEM simulations validated by physical model tests, numerical ice tanks can be developed as an effective complement to their counterparts. In this paper, a numerical ice tank based on 3D spherical DEM was established with respect to the small ice model basin of China Ship Scientific Research Center (CSSRC-SIMB). Based on spherical DEM with parallel bond model, the model tests of typical structures (vertical cylinder and inclined plate) in level ice sheets were established in the numerical ice tank, and the ice–structure interaction process under the same initial conditions was simulated. The accuracy of the simulations is verified by comparing the simulated ice loads with the measured ice loads from the model tests in the CSSRC-SIMB. Furthermore, the application of the numerical ice tank was extended to simulate the navigation of a Wass bow in level ice and broken ice conditions. The value of the break resistance of the Wass bow in level ice was evaluated, and the numerical ice tank produced results that were found to be consistent with those obtained from Lindqvist’s formula. The statistical properties of the bow load for different broken ice fields with the same initial physical conditions are analyzed by performing a repeatability test on the broken ice fields.

1. Introduction

In recent years, the advent of climate change has brought about a number of different changes in the Arctic, including accessibility to the opening of new Arctic shipping lanes and development of oil and gas resources in polar regions. Among them, various countries have vigorously developed polar marine engineering equipment and actively participated in Arctic oil and gas development and transportation projects. For different polar marine engineering equipment, the determination of ice load on polar ships and offshore structures is important for ice-resistant design, safe operation, and structural integrity management in ice-infested waters.

The ice load of polar ships and offshore structures in cold regions is influenced by sea environment (temperature, salinity, wind field, ocean current, etc.), ice type (floe ice, broken ice, level ice, ridge, ice rubble, etc.), structural motion characteristics (rigid and elastic bodies, fixed and floating), and other factors [1,2]. The most reliable method to determine the ice load on marine engineering structures is full-scale measurements. For example, ice load measurement projects for various marine structures (jacket platforms in the Bohai Sea [3], lighthouses in the Baltic Sea [4], Confederation Bridge in Canada [5], caisson platforms [6], polar icebreakers [7], etc.) were carried out during the last century. However, expensive installation costs and time make it difficult to implement. Correspondingly, the implementation of physical model tests in an ice tank/basin is an attractive alternative to determine the ice load on marine structures. The design of the physical model testing in an ice tank/basin usually adopts the Froude and Cauchy scaling laws to ensure the geometry, kinematics, and dynamics of the test object (structure) in the field ice environment, which are far more complicated than those in open water. The document of General Guidance and Introduction to Ice Model Testing in the International Towing Tank Conference (ITTC) presents information on the nine main overviews of grain structure and chemical composition of the model ice in operating facilities [8]. Due to the different production methods, model ice can be divided into fine-grained ice and columnar ice. The fine-grained ice is represented by Aalto University (Finland) and Aker Arctic (Finland), while columnar ice is represented by the Hamburg Ship Model Basin (Germany) and National Research Council Canada, Ocean, Coastal and Engineering (Canada). In general, physical model tests have been widely used to investigate global ice action and structural ice loads. For example, studying the ice–structure interaction mode (ice failure mode, ice pile-up, ice accumulation, ice blockage, etc.) and providing an icebreaker ship design scheme (maneuverability, stability, rapidity, design load, etc.) [9,10,11].

Due to the high cost and time of each test run, it is not common to perform multiple repetitions of a physical model test under the same conditions. However, when analyzing test results from ice tank tests, it is very important that the test results be reproducible or repeatable [12,13]. The repeatability of test results is closely related to the random uncertainty of the test, with the higher test repeatability corresponding to lower random uncertainty. Although the Procedure for assessing the experimental uncertainty in ship resistance testing in ice document [14], provided by the ITTC, has described how to independently evaluate the random uncertainty of the results by splitting up the time histories of measured results from a single test, the differences of test conditions in each segment will still affect the assessment of the random uncertainty of the final results. It would be much more convenient for us in the interpretation of the results obtained from the test of the ice tank if we could repeat it.

With the rapid development of computer technology, some numerical methods have been able to reproduce well the phenomena and results of ice tank tests, providing a solution to the problem of the repeatability of ice tank tests. Numerical simulations have potential to overcome their repeatability limitations due to their stepwise nature and capability for extensive computation [15]. For example, Jeon and Kim (2021) [16] and Wang et al. (2018) [17] applied the finite element method (FEM) with the element erosion technique to simulate the model tests of the interaction between conical structure with level ice, which was performed by Kärnä et al. (2010) [18] in the Hamburg Ship Model Basin (HSVA). The applicability of this approach is verified by the analysis of ice load and ice failure length between numerical simulations and model tests. Long et al. (2020) [19] and Jang and Kim (2021) [20] also performed the discrete element method (DEM) to simulate the model tests of conical structure, which all verify the accuracy of the DEM in analyzing the ice load and process between ice and structure. In addition, smoothed-particle hydrodynamics (SPH) method [21], moving particle semi-implicit method (MPS) [22], extended finite element method (XFEM) [23], peridynamics (PD) [24], semi-empirical formula method [25], etc., are also used in ice–structure simulation of the model tests or full-scale tests.

This study proposes an example of a numerical ice tank based on the DEM, which is not only limited to the numerical repeatability of physical ice tank tests, but can be extended to simulation scenarios and calculational parameters over some physical environments [26]. DEM simulations proposed by Cundall and Strack (1979) [27] have been widely used in ice engineering research, such as demonstrating ice failure behaviors, structural ice loads, ice rubble pile-up, and sea ice drift. Among them, the elemental morphology of DEM has also developed from 2D disk to complex 3D shapes such as spheres, polyhedral, and extended polyhedral [28,29,30]. Compared to polyhedron elements, 3D sphere elements can easily establish high-performance numerical algorithms for a numerical ice tank based on the graphics processing units (GPU) technology to achieve rapid simulation [31]. In addition, considering the coupled effect of structural dynamic response or the wave/current and sea ice damage and movement, the coupled method of 3D spherical DEM with FEM [32] or computational fluid dynamics (CFD) [33] has also been developed accordingly.

In the following contest, Section 2 sequentially states the model of numerical ice tank corresponding to the small ice model basin of China Ship Scientific Research Center (CSSRC-SIMB), the ice model based on 3D spherical DEM, and interaction model between model ice and offshore structure. Section 3 introduces the physical model tests for typical offshore structures represented by a vertical cylinder and inclined plate and compares them with the results obtained from numerical ice tank tests. In Section 4, the numerical ice tank based on the CSSRC-SIMB is extended to simulate the navigation of the new Wass bow in both level ice and broken ice field conditions. Finally, Section 5 concludes the paper.

2. Numerical Ice Tank Description Based on DEM

2.1. Numerical Ice Tank Corresponding to the CSSRC-SIMB

The main body of the CSSRC-SIMB has a dimension of 8 m in length, 2 m in width, and 1 m in depth [34,35]. The ice tank is housed in an insulated room that can be cooled down to an air temperature of −16 °C, as shown in Figure 1a. In addition, a small cold room is close to the ice tank, with a size of 3 m (length) × 2.2 m (width) × 2.8 m (height), which is mainly used to test the physical and mechanical properties of model ice. The main carriage on the tank was designed for tractive force or thrust up to 3 KN and speeds ranging from 0.01 m/s to 1 m/s. The model ice used is naturally grown columnar ice constructed from sodium chloride solutions. Another in situ test of offshore structures was designed by Froude and Cauchy scaling laws to carry for physical model testing. According to the dimensions and ice modeling methodology, the physical and mechanical parameters of ice in both full and model scales are listed in

Table 1. Physical and mechanical parameters of ice at the full and model scales.

Parameter	One-Year Columnar Sea Ice	Model Ice
Thickness (m)	0.5~2.0	0.03~0.05
Density (g/cm3)	0.91	0.9~0.92
Young’s modulus (GPa)	2~5	0.06~0.2
Compression strength (MPa)	0.5~12	0.05~0.2
Flexure strength (MPa)	0.5~2	0.03~0.1
Tensile strength (MPa)	0.2~0.8	0.01~0.03

. Based on the Froude–Cauchy similarity, the scale ratio in the full and model scales can be controlled between 10 and 20 in the CSSRC-SIMB.

Using the physical ice tank of the CSSRC-SIMB as a template, this paper describes a numerical ice tank based on the DEM to conduct virtual offline-mode tests of various offshore structure, as shown in Figure 1b. The geometric model and dimensions of the numerical ice tank are the same as those of the physical ice tank, and the physical and mechanical properties of the model ice and offshore structure in numerical and model tests are also the same, that is, both use the Froude–Cauchy similarity in the scheme design. A 3D spherical DEM code from Dalian University of Technology was used for the model ice simulation in the numerical ice tank [31,36]. The details of the method and the interaction model for the ice–structure in the numerical ice tank are presented below.

2.2. Ice Model Based on the 3D Spherical DEM in the Numerical Ice Tank

The model ice in the numerical ice tank is simplified to a particle system of 3D spherical elements with the same size and mass, which are regularly arranged in space according to Hexagonal Close Packing (HCP) arrangement, as shown in Figure 2. Essentially, the mechanical behavior of the numerical model ice is described by the motion of the elements and the forces and moments acting between them. The motion is performed per explicit calculation by Newton’s second law, and the contact force between two the elements is calculated using an elastic-viscous contact model based on the springs and dashpots. In addition, all 3D spherical elements are bonded together using a parallel bonding model to take into account the macroscopic continuum properties of the model ice [37]. In the parallel bonding model, a parallel-bonded disk with circular cross-section is set up to transfer through forces and moments between two adjacent elements. The maximum normal and shear stresses acting on the profile of the parallel-bonded disk are determined using beam theory and can be written as follows:

$${\sigma }_{\max }=\frac{-F_{\mathrm{n}}}{A}+\frac{\left|M_{\mathrm{s}}\right|}{I}R,{\tau }_{\max }=\frac{\left|F_{\mathrm{s}}\right|}{A}+\frac{\left|M_{\mathrm{n}}\right|}{J}R$$

(1)

where $F_{\mathrm{n}}$ and $F_{\mathrm{s}}$ are the normal and shear forces, respectively; $M_{\mathrm{n}}$ and $M_{\mathrm{s}}$ are the normal and shear moment, respectively; $R$ is the radius of spherical element; $A=\pi R^2$, $I=\frac{1}{4}\pi R^4$, and $J=\frac{1}{2}\pi R^4$ are the area, moment of inertia, and polar moment of inertia of the parallel-bonded disk, respectively.

The macroscopic failure of model ice is closely related to the generation and propagation of internal cracks. In the DEM of the model ice, the failure of the bonding disk between elements simulates the generation of internal cracks in the material. During the failure process, the bond failure modes between elements can be classified into tensile failure and shear failure based on the normal and tangential stresses of parallel-bonded disk. When the maximum tensile stress (${\sigma }_{\max }$) exceeds the tensile strength (${\sigma }^{\mathrm{t}}$) between adjacent elements, ${\sigma }_{\max }>{\sigma }^{\mathrm{t}}$, the parallel-bonded disk exhibits tensile failure; when the maximum shear stress (${\tau }_{\max }$) exceeds the shear strength (${\tau }^{\mathrm{s}}$) between adjacent elements, ${\tau }_{\max }>{\tau }^{\mathrm{s}}$, the parallel-bonded disk exhibits shear failure, as shown in Figure 3. Based on the Mohr–Coulomb shear friction law, the tensile and shear strengths (${\sigma }^{\mathrm{t}}$ and ${\tau }^{\mathrm{s}}$) can be determined as follows:

$${\sigma }^{\mathrm{t}}={\sigma }_{\mathrm{b}}^{\mathrm{n}}$$

(2)

$${\tau }^{\mathrm{s}}={\sigma }_{\mathrm{b}}^{\mathrm{s}}+{\mu }_{\mathrm{b}}·\max (0,{\sigma }_{\max })$$

(3)

where ${\mu }_{\mathrm{b}}$ is the friction coefficient between the bonding elements; ${\sigma }_{\mathrm{b}}^{\mathrm{n}}$ and ${\sigma }_{\mathrm{b}}^{\mathrm{s}}$ are the normal and shear bonding strengths of the parallel-bonded disk, respectively.

2.3. Interaction Model of Ice–Structure in Numerical Ice Tank

Due to the complex geometrical configuration of ships and offshore structures, discretizing the plane or surface shape of typical offshore structures into triangular panels is the most commonly used method to describe complex geometries, as shown in Figure 4. Simplifying the ice–structure interaction as a contact problem between a sphere and a triangular panel, the ice load on the structure can be obtained by setting a contact criterion between the sphere and the triangular panel.

The contact search between the sphere and the triangular panels can be determined by the position of the sphere center point $P$ and triangle vertex $A$, $B$, and $C$ in space. The contact type between the sphere and triangular panels can be subdivided into three contact modes: sphere–face, sphere–edge, and sphere–vertex, as shown in Figure 5. The projection point (contact point) of the points $P$ on the plane where the triangular panels are located is $Q$. The distance between the vector $\overrightarrow{PQ}$ and radius $R_{\mathrm{DEM}}$ of the sphere determines the contact and contact overlap amount between the sphere and the plane. If $\left|\overrightarrow{PQ}\right|<R_{\mathrm{DEM}}$, it indicates that there is an overlap between the sphere and plane, and the amount of overlap $∆L$ can be calculated as follows:

$$∆L=\left|\overrightarrow{PQ}\right|-R_{\mathrm{DEM}}$$

(4)

The normal vector ${\mathit{n}}_{\mathrm{w}}$ between the sphere and contact point of the triangular panels can be expressed as:

$${\mathit{n}}_{\mathrm{w}}=\frac{\overrightarrow{PQ}}{\left|\overrightarrow{PQ}\right|}$$

(5)

Relative displacement $∆\mathit{x}$ can be expressed as:

$$∆\mathit{x}=\left({\mathit{v}}_{\mathrm{DEM}}-{\mathit{v}}_{\mathrm{w}}\right)∆t$$

(6)

where ${\mathit{v}}_{\mathrm{DEM}}$ is the speed of the sphere, ${\mathit{v}}_{\mathrm{w}}$ is the speed of the triangular panels, and $∆t$ is the timestep of DEM. The relative displacement along normal vector ${\mathit{n}}_{\mathrm{w}}$ is decomposed into normal displacement $∆{\mathit{x}}_{\mathrm{n}}$ and tangential displacement $∆{\mathit{x}}_{\mathrm{s}}$:

$$∆{\mathit{x}}_{\mathrm{n}}=\left({\mathit{n}}_{\mathrm{w}}·∆\mathit{x}\right){\mathit{n}}_{\mathrm{w}},∆{\mathit{x}}_{\mathrm{s}}=∆\mathit{x}-∆{\mathit{x}}_{\mathrm{n}}$$

(7)

The interaction force between sphere and boundary elements is calculated using the linear contact model, in which the normal force ${}_{}{}^{\mathrm{t}}\mathit{F}{}_{\mathrm{n}}^{}$ and tangential force ${}_{}{}^{\mathrm{t}}\mathit{F}{}_{\mathrm{s}}^{}$ at time $t$ can be expressed as:

$${}_{}{}^{\mathrm{t}}\mathit{F}{}_{\mathrm{n}}^{}=k_{\mathrm{w}}^{\mathrm{n}}∆{\mathit{x}}_{\mathrm{n}}$$

(8)

$${}_{}{}^{\mathrm{t}}\mathit{F}{}_{\mathrm{s}}^{}=\min \left({}_{}{}^{\mathrm{t}-∆\mathrm{t}}\mathit{F}{}_{\mathrm{s}}^{}+k_{\mathrm{w}}^{\mathrm{s}}∆{\mathit{x}}_{\mathrm{s}},{\mu }_{\mathrm{w}}{\mathit{F}}_{\mathrm{n}}\right)$$

(9)

where $k_{\mathrm{w}}^{\mathrm{n}}$ and $k_{\mathrm{w}}^{\mathrm{s}}$ are the normal and tangential contact stiffness between the sphere and triangular panels, respectively. ${\mu }_{\mathrm{w}}$ is the friction coefficient of ice–structure. ${}_{}{}^{\mathrm{t}-∆\mathrm{t}}\mathit{F}{}_{\mathrm{s}}^{}$ is the tangential force at time step $\mathrm{t}-∆\mathrm{t}$. The total ice load on the structure ${\mathit{F}}_{\mathrm{total}}$ is expressed as:

$${\mathit{F}}_{\mathrm{total}}={\displaystyle \sum }_{i=1}^N\left({}_{}{}^{\mathrm{t}}\mathit{F}{}_{\mathrm{s}}^i+{}_{}{}^{\mathrm{t}}\mathit{F}{}_{\mathrm{n}}^i\right)$$

(10)

where $N$ is the numerical of interaction force between sphere and boundary elements.

3. Verification of Physical Model Tests in the Numerical Ice Tank

This section describes model tests of two typical offshore structures (vertical cylinder and inclined plate) carried out on physical and numerical ice tanks. By comparing the failure modes of the model ice and the structural load characteristics of two, the correctness of the construction of the numerical ice tank based on DEM is verified.

3.1. Description of Physical Model Tests

Model testing of typical offshore structures represented by vertical cylinder and inclined plate was performed at the CSSRC-SIMB as part of the exploration of the model testing methods. The similarity between model ice and sea ice is emphasized in the design of relevant tests to ensure that the physical and mechanical parameters of sea ice are within a certain range, while the size and carriage speed of structures do not need to strictly meet the similarity law of model tests. Considering the test capacity of the CSSRC-SIMB, the model tests with vertical cylinder and inclined plate structures are selected as test models, as shown in Figure 6. The relevant test information for each condition, including sea ice strength, ice density, speed, structural geometry, etc., is listed in

Table 2. Test matrix and ice conditions of interaction between the vertical cylinder and level ice.

Test No.	Compressive Strength (kPa)	Ice Density (kg/m3)	Ice Thickness (mm)	Diameter of Cylinder (mm)	Speed (mm/s)
#101	51.6	901	36	100	50
#102					100
#103					150
#201	55.3	902	37	150	50
#202					100
#203					150
#301	57.3	915	37	200	50
#302					100
#303					150

and

Table 3. Test matrix and ice conditions of interaction between the inclined plate and level ice.

Test No.	Flexural Strength (kPa)	Ice Density (kg/m3)	Ice Thickness (mm)	Angle of Inclined Plate (°)	Speed (mm/s)
#401	33.4	924	37	120	50
#402					100
#403					150

. The former is used to analyze the brittle crushing of the ice and the latter is used to analyze the bending failure of the ice.

3.2. Comparison of the Results

In order to verify the validity of the numerical ice tanks presented in this paper, a visual and numerical comparison of the simulation results of the above model tests was performed. In addition to the numerical comparison, the comparison of visual results during the ice–structure interaction is the main purpose of the construction of the numerical ice tank, which is a virtual and repeated model test scenario. The computational parameters of the numerical simulations in the numerical ice tank are listed in

Table 4. Computational parameters of the numerical simulations in the numerical ice tank.

Definition	Symbol	Value
Young’s modulus of elements	$E_{\mathrm{p}}$	7.5 MPa
Friction coefficient of elements	${\mu }_{\mathrm{b}}$	0.2
Normal and shear bonding strengths	${\sigma }_{\mathrm{b}}^{\mathrm{n}}$$,{\sigma }_{\mathrm{b}}^{\mathrm{s}}$	62 kPa
Ice–structure friction coefficient	${\mu }_{\mathrm{w}}$	0.15
Number of elements	N	56 W

.

A visual comparison of the results for test cases #202 and #402 is presented in Figure 7. Figure 7 shows the interaction process, including crushing failure, bending failure, and channel ice accumulation, that occurs in both physical and numerical ice tanks. During the interaction process between the ice sheet and vertical cylinder, the ice sheet is subjected to crushing brittle failure, and the broken ice fragments are extruded from the upper and lower surfaces of the ice sheet to form a temporary accumulation body. During the interaction between the ice sheet and inclined plane, the ice sheets have periodic bending failure and the length of each failure of the ice sheets is basically the same. In addition, the phenomenon of sea ice climbing and accumulation during the action is very obvious.

Time history of the ice loads for vertical cylindrical structures with different diameters is presented in Figure 8. In both the physical and numerical ice tank results, the ice load of the structure presents continuous brittle extrusion and crushing, and the load presents loading and unloading form but without periodicity. Meanwhile, comparison of time history of ice loads for inclined plate structures with different speeds is shown in Figure 9. The dynamic load on the inclined plate is attributed to the continuous bending failure that occurs around the front of the structure. During each loading and uploading period, the ice load gradually increases over time and then drops rapidly once it reaches a certain level.

The mean load and maximum loads are comparable in magnitude. The mean and maximum values of the ice load for the vertical cylindrical structures in the simulations and tests are compared and analyzed. Here, the mean value and maximum value are the time histories of the ice load during the stable period, as shown by the dotted line in Figure 9. The mean peak load and maximum load of the plate structure are also analyzed. The peak load represents the maximum value of the structure within a single loading and unloading cycle.

Table 5. Comparison of the ice loads between the simulations and tests with different conditions.

Structure	Test No.	Mean Load (Mean Peak Load)/N			Maximum Load/N
		Simulation	Test	Relative Error	Simulation	Test	Relative Error
vertical cylinder	#101	90.6	92.3	1.9%	233.2	246.8	5.5%
	#102	105.6	110.4	4.4%	247.4	253.3	2.3%
	#103	145.6	116.7	24.7%	299.8	220.9	35.7%
	#201	105.4	119.7	12.0%	260.6	260.5	0.1%
	#202	117.6	128.8	8.7%	283.0	308.8	8.4%
	#203	157.2	153.2	2.6%	326.0	302.0	7.9%
	#301	110.8	151.1	26.7%	286.9	334.0	14.1%
	#302	165.4	178.3	7.2%	297.6	311.0	4.3%
	#303	183.0	204.8	10.6%	342.9	357.8	4.2%
	Mean error	-	-	11.0%	-	-	9.2%
inclined plate	#401	115.2	152.7	24.6%	172.2	187.1	8.0%
	#402	194.0	176.1	10.2%	225.3	206.1	9.3%
	#403	204.2	187.5	8.9%	286.8	259.7	10.4%
	Mean error	-	-	14.5%	-	-	9.2%

summarizes the statistical characteristics and relative errors of the ice loads for the 12 working conditions. Although the relative error between the test and the simulation is relatively large in some cases (#103 and #301), the mean load (mean peak load) and the maximum load error of the numerical ice tanks are still within the acceptable range. These sources of errors are not only unreasonable choices of numerical simulation parameters, but also random uncertsainties of some experiments.

4. Application of Numerical Ice Tank

4.1. Breaking Resistance Evaluation of New Wass Bow

The Wass bow takes its name from Dr. Heinrich Wass of Germany [38]. The main characteristics of the simplified Thyssen Wass bow include a small stem angle, a large waterline angle, and a large outward drift angle, which can increase the proportion of bending failure components of the ice sheet and greatly reduce the ice breaking resistance when the ship is sailing through the ice region. In this section, a new Wass bow test model is designed with the simplified Wass bow as the mother ship, as show in Figure 10a. The numerical ice tank model above is used to evaluate the ice breaking resistance of the bow under different ice thickness conditions, as shown in Figure 10b. At full-scale test, the ice breaking conditions for the new Wass bow are as follows: ice flexural strength is 600 kPa, ice thickness is 0.6~1.5 m, and ship speed is 2 kn. The scale of the model test set was set to 30, based on the size of the ship’s bow model and CSSRC-SIMB. The choice of the discretization element parameters in the numerical ice tank refers to the parameters in Table 4, and the other parameters for the relevant model tests are listed in

Table 6. Model test parameters for the new Wass bow in the numerical ice tank.

Definition		Symbol	Value
Wass bow model	Length	$L_{\mathrm{ship}}$	0.96 m
	Front/Back breadth	$B_{\mathrm{ship}}$	0.40/0.46 m
	Depth	$D_{\mathrm{ship}}$	0.30 m
	Draft	$T_{\mathrm{ship}}$	0.18 m
	Stem angle	$\phi$	24 deg
	Waterline angle	$\alpha$	84 deg
	Drift angle	$\beta$	52 deg
	Ice–structure friction coefficient	${\mu }_{\mathrm{w}}$	0.1
	Speed	$v_{\mathrm{ship}}$	0.188 m/s
Model ice	Thickness	$h_{\mathrm{ice}}$	0.02~0.05 m
	Flexural strength	${\sigma }_{\mathrm{ice}}$	20.0 kPa
	Density	${\rho }_{\mathrm{ice}}$	920 kg/m3
	Poisson’s ratio	$\nu$	0.3
	Young’s modulus of elements	$E_{\mathrm{p}}$	4.5 MPa
	Normal and shear bonding strengths	${\sigma }_{\mathrm{b}}^{\mathrm{n}}$$,{\sigma }_{\mathrm{b}}^{\mathrm{s}}$	37 kPa
Scale		$\lambda$	30

.

A series of model tests of the Wass bow breaking the level ice were carried out in the numerical ice tank. Figure 11 shows a simulation snapshot of the interaction between ice sheet and Wass bow at an ice thickness of 0.03 m. The ship bow interacts with the ice sheet at a constant carriage speed, resulting in crushing and bending failure of the ice sheet, as shown in Figure 11a. The bending failure of ice plate is the most obvious, and the development of circumferential and radial cracks can be observed during the failure process. In addition, the mode of failure of the ice sheet is dominated by bending failure caused by circumferential and radial cracks, as shown in Fig. 11b. The same failure phenomenon has been observed in other simulations of the Wass bow-ice interactions [24]. In Zhang’s paper, the simplified Wass bow is also used for the icebreaker, which is the same design as the ship bow selected in this paper. Figure 11c also shows the shape of the ice channel that was created behind the ship bow.

The time histories of ice force on the ship bow in the three directions are shown in Figure 12. With the bow continuously entering the ice area, its ice force also gradually increases. About 5 s later, the bow completely enters the ice area and the variation trend of its ice force tends to be stable. Since then, it can be seen that the ice force of the ship bow has periodic load characteristics due to the continuous bending failure of ice sheet. The ice force in the x direction represents the breaking resistance of the bow, which consists of two parts: bending resistance and crushing resistance. The mean value of breaking resistance in this state is analyzed and compared, which is the red dotted line in Figure 12a.

The validity of the numerical ice tank results can be verified by empirical formulas, thus improving the ability of the numerical ice tank to comprehensively describe the ice load characteristics. In this paper, Lindqvist (1989) [39] ice resistance formula is selected to check the rationality of the numerical ice tank results. The ice breaking resistance of the ship bow in Lindqvist’s formula can be divided into bending resistance and crushing resistance, which can be expressed as:

$$R_{\mathrm{ice}}=\left(R_{\mathrm{c}}+R_{\mathrm{b}}\right)\left(1+\frac{1.4v_{\mathrm{ship}}}{\sqrt{g{h}_{\mathrm{ice}}}}\right)$$

(11)

$$R_{\mathrm{c}}=0.5{\sigma }_{\mathrm{ice}}{h}_{\mathrm{ice}}^2\frac{\tan \phi +\mu \cos \phi /\cos \psi }{1-\mu \sin \phi /\cos \psi }$$

(12)

$$R_{\mathrm{b}}=\frac{27}{64}{\sigma }_{\mathrm{ice}}{B}_{\mathrm{ship}}\frac{h_{\mathrm{ice}}{}^{1.5}}{\sqrt{\frac{E}{12\left(1-{\upsilon }^2\right)g{\rho }_{\mathrm{w}}}}}\frac{\tan \psi +\mu \cos \phi }{\cos \psi \sin \alpha }\left(1+\frac{1}{\cos \psi }\right)$$

(13)

where $R_{\mathrm{ice}}$, $R_{\mathrm{c}}$, and $R_{\mathrm{b}}$ is the breaking resistance, crushing resistance, and bending resistance, respectively. $\psi$ is normal angle, $\psi =\arctan \left(\tan \phi /\sin \alpha \right)$. $E$ and $\upsilon$ is Young’s module and Poisson’s ratio of ice. ${\rho }_{\mathrm{w}}$ is the density of water, ${\rho }_{\mathrm{w}}=1035$ kg/m³. $g$ is the gravitational acceleration.

The mean breaking resistance results of the numerical ice tank on the Wass ship bow under different ice thickness are shown in Figure 13. It can be observed that the trend of the breaking resistance is the same as that obtained with Lindqvist’s formula. Overall, the numerical ice tank results are slightly lower than those calculated using the formula. The results for the ice bath are slightly lower than those calculated by the formula but the errors are within an acceptable limit. In summary, the numerical ice tank model presented in this paper has certain advantages for the evaluation of virtual ice-breaking resistance tests.

4.2. Repeatability of Broken Ice Tests of New Wass Bow

The repeatability of ice tank tests with broken ice field are important for assessing the random uncertainty [40]. This is not often performed in physical ice tanks but can be conveniently implemented in numerical ice tanks. The broken ice field in the ice tank was created by cutting the ice sheet into square floes. The square floes with sides ranging from 0.2 to 0.5 were chosen to create a broken ice field with 70% concentration in the CSSRC-SIMB, as shown in Figure 14a. In the numerical ice tank, the broken ice fields were created by randomly filling the square floes with the same size distribution and shape in the CSSRC-SIMB, and the randomly filling method is similar to the method used in Metrikin’s paper [41], as shown in Figure 14b. Each square floe is composed of 3D spherical elements with the same physical and mechanical parameters as the level ice in Section 4.1. By randomly changing the position and angle of each floe, in Figure 15, three numerical broken ice fields are constructed under the same ice condition (floe size and concentration) to carry out repeatability tests of new Wass bow in broken ice fields.

The repeatability of the numerical ice tank test was assessed by performing three numerical simulation runs with the same initial conditions, including ice thickness of 0.03 m, carriage speed of 0.188 m/s, floe size of 0.2 to 0.5, flexural strength of 20.0 kPa, and ice concentration of 70%. The only difference between each simulation is the change in position and orientation in each broken ice field. Figure 16 shows a simulation snapshot of the interaction process between the bow and the floe at t = 20 s. During the towing process, the ship’s bow inadvertently collided with some broken ice. The impact caused the floes to be displaced and subsequently crushed, bent, and split. As a result, the floes exhibited failure. Additionally, it was noticed that there was a clear channel at the stern of the bow. Figure 17 also shows the shape of the ice channel that was created behind the ship bow in three numerical broken ice fields.

Comparisons of the ice force results of three numerical broken ice fields indicate that the ice resistance of case 1 and case 3 increases obviously near the end point of model test, as shown in Figure 18. In addition, in case 2, ice resistance is relatively stable throughout. The event that affects the ice resistance is jamming of ice. Jamming behavior occurred during the near end points of case 1 and case 3, as shown in Figure 17. It causes floes to accumulate near the bow of the ship and, as the ship continues to move forward, ice resistance becomes significantly greater.

In the initial phase, the ice force in the first 5 s of the bow is not taken into account, as the bow gradually enters the region of broken ice, which will cause the ice load to gradually increase and become unstable. Data from the next 30 s were selected to calculate the mean value and standard deviation of the ice force in each operating condition. The standard deviation is 1.39, 1.28, and 1.38 in case 1, case 2, and case 3, respectively. The average ice force is 5.11 N, 4.79 N, and 5.14 N in case 1, case 2, and case 3, respectively. The above results show that there are significant differences in the time history characteristics of ice force under the same initial field conditions. The accumulation of floes causes an increase in the ice load on the structure, but its generation is random. In case 2, there is no accumulation when the ship moves forward, so there was no increase process. The mean value and standard deviation of the force are smaller than in the other two cases. Uncontrolled conditions of tests in broken ice (the initial field of floes) can cause significant differences to the time history characteristics of the structure’s force. Therefore, the repeatability of the numerical ice tank tests is important for evaluation of loads from broken ice, and numerical ice tank models provide a solution to this problem.

5. Conclusions

In this study, a numerical ice tank model corresponding to the CSSRC-SIMB is proposed to investigate the interaction process between offshore structures and ice using the 3D spherical DEM. A series of ice tank test cases including three structure types (vertical cylinder, inclined plate, and new Wass bow) were analyzed in the numerical ice tank. Numerical reproduction of the model tests in the numerical ice tank provides convenient information on offshore structure–ice interaction processes, ice failure models, and time histories of ice force on the structures. The results of this study support the following conclusions:

(1)

The numerical ice tank model is a reasonable method to perform repetitions of the model test in a physical ice tank/basin, which can reproduce ice failure performance, including curding, bending, accumulation, and dynamic ice load on offshore structures.

(2)

The visual and numerical comparison of the simulation results of 12 ice tank tests demonstrates the rationality of the numerical ice tank. The mean load (mean peak load) and the maximum load error of the numerical ice tanks are still within the acceptable range.

(3)

The numerical ice tank model was used to evaluate the breaking resistance of the new Wass bow. In tests with various ice thicknesses, the numerical ice tank produced results that were found to be consistent with those obtained from Lindqvist’s formula.

(4)

The repeatability of Wass bow tests with broken ice fields is important for evaluating the time history characteristics, and uncontrolled conditions of tests in broken ice (the initial field of floes) can cause the ice load to gradually increase and become unstable.

In summary, the work in this paper presents a systematic introduction to the method, validation, and application of the numerical ice tanks. In numerical ice tanks, simplified hydrodynamic treatment of structures and model ice may reduce the accuracy of structural ice load prediction, which is also a limitation of this paper. Computational hydrodynamics and discrete-element coupling methods are the focus of our subsequent work on accurate modeling of numerical ice tanks.