Study on the Dynamic Ice Load of Offshore Wind Turbines with Installed Ice-Breaking Cones in Cold Regions

: The dynamic ice load of conical offshore wind turbines (OWTs) in cold regions is unclear. The ice force period is the key parameter used to establish an ice force model for conical structures. To obtain ice load data, a ﬁeld monitoring system was installed on an OWT in a cold region in China. Based on the monitoring data, a new formula for calculating the ice force period of conical structures was established. By comparing the period calculated with this formula and the measured ice force period, it was found that the calculated data generally agreed with the measured data. Then, a random dynamic ice force model for conical OWTs can be established. Based on this ice force model, the ice-induced vibration of an OWT was analyzed with the ANSYS ﬁnite element software. The results are in good agreement with the measured data obtained from the OWT in the time and frequency domains. Therefore, the random dynamic ice force model established in this paper can be used to evaluate the ice resistance performance of conical OWTs in cold regions.


Introduction
Marine wind energy resources have greatly expanded in recent years.Excellent progress has also been made in the design, construction and installation of large and superlarge offshore wind turbine (OWT) structures.Complex marine environmental loads are the main factors that affect the structural safety of OWTs, especially for wind farms built in ice-covered sea areas.Wang et al. analyzed the influence of ice-induced vibrations (IIVs) on OWTs in the Bohai Sea through a numerical method, and found that the acceleration response at the top of the OWT was significantly caused by IIV [1].Lou analyzed and compared the fatigue damage caused by IIVs and the fatigue damage caused by wind load on OWTs in the Bohai Sea, and found that they were equivalent [2].Based on the Matlock model and the Määttänen model, Ye et al. simulated the IIVs on OWTs, and found that IIVs had a great impact on OWTs, especially for the tower [3].Therefore, ice loads are one of the controlling loads for OWTs.Hendrikse et al. studied ice interactions between ice and vertical OWTs when the frequency lock-in phenomenon occurred [4].Based on the existing theory, Seidel et al. established an evaluation method for the frequency lock-in phenomenon of OWTs [5].At present, the ice load influences on OWTs and the ice-resistant design of OWTs have become key issues for wind farms in cold regions.When sea ice interacts with different types of structures, the breaking failure characteristics of sea ice can differ.According to the regulations of Q/HSn 3000-2002 [6], the design compression strength of sea ice in China's sea areas is greater than the design bending strength, and most of the design compression strength is about three times greater than the design bending strength.Therefore, to reduce the influence of sea ice on vertical structures, a sloped plane or cone was installed on the vertical structure [7].Installing a cone at the water level of an OWT may provide ice-resistant effects.However, the dynamic ice load is unclear for cones installed on OWT structures.Further research on the dynamic ice loads of OWTs is needed.
Studies of ice loads on conical structures began in the mid-1960s [8].Based on theoretical analysis, a computational method for ice loads on cone structures was established by Crosdale and Ralston et al. [9,10].Then, based on laboratory tests and field monitoring data, the ice load computational method for cone structures was established by Kato, Hirayama and Obara, and Yu et al. [11][12][13].In laboratory tests, Sodhi and Frederking both obtained ice force time history data for cone structures [14,15].Based on the ice force time history data obtained from field measurements of the JZ20-2MUQ platform in the Bohai Sea [16,17], a deterministic ice force function with three parameters was established by Yue et al. [18].The relationships among the ice force periods of conical structures, ice velocities and ice thicknesses were obtained.The statistical characteristics of ice force amplitudes and ice force periods based on the measured data of ice forces obtained on ice-loaded panels were proposed by Qu et al. [19].Based on the data obtained from field monitoring, Yue et al. [20] developed a random ice force spectrum of conical structures.The development of a discriminant method for narrow conical structures and the effect of ice-breaking cones (IBCs) on reducing ice loads and ice-induced vibrations was proposed by Xu [17].At present, research on the dynamic ice loads of conical structures has been mainly based on jacket oil and gas platforms or lighthouse structures.The natural frequency of jacket structures is approximately 1~2 Hz, the height above the water level is generally 20 m and the diameter of ice-breaking cones is generally 4~6 m.However, the natural frequency of OWTs is generally 0.2~0.3Hz, heights above the water level are generally above 60 m and the diameters of IBCs installed on OWTs are approximately 9 m.With the application of higher-power OWTs, tower diameters are increasing.As a result, IBC diameters are also increasing.Hence, interactions between sea ice and conical OWTs need to be studied, and the applicability of the existing ice load models of conical structures should be assessed.
In this paper, data on ice load and structural response were obtained from the field monitoring of conical OWTs in cold regions.The ice breaking features caused by different conical structures are analyzed.A new function for the periods of dynamic ice loads for conical structures is proposed, and a predictive ice load model suitable for conical OWTs is developed.The predicted ice load function is verified by field measurement data from OWTs.Finally, the ice resistance performance of conical OWTs is analyzed based on an ice force model, which can provide guidance for the development and design of wind farms in cold regions.

Field Monitoring of an OWT in Ice Regions
To obtain data on sea ice conditions, ice load, hydrology and structural response, a field monitoring system was installed on an OWT in the Zhuanghe sea area in the North Yellow Sea of China in the winter of 2017-2018 (39 • 32.5 N, 123 • 19 E).The structure is a monopile foundation OWT, and a forward and backward ice-breaking cone was installed at the waterline.Its basic properties are shown in Table 1.The field monitoring system (FMS) consists of structural-response monitoring system and a sea-ice-conditions monitoring system.The integrated monitoring system of OWT foundations is shown in Figure 1a.The structural-response monitoring system consists of four sets of acceleration tilt sensors (Figure 1b), and the sea-ice-conditions monitoring system consists of two sets of ice-image cameras (Figure 1c,d, model number: DS-2CD2T86FWDV2-I5S).The cameras conform to the main current direction at 180-degree angles from each other.Each system includes two cameras: the first set for monitoring the ice speed and ice conditions, and the second set for monitoring the ice thickness and ice breaking.The ice thickness and ice velocity data were obtained from the cameras based on the broken ice pieces that had turned over.The exposed section can be captured by a camera, as shown in Figure 2. The yellow part is the ice-broken cone of the OWT.The structural-response monitoring system consists of four sets of acceleration tilt sensors (Figure 1b), and the sea-ice-conditions monitoring system consists of two sets of ice-image cameras (Figure 1c,d, model number: DS-2CD2T86FWDV2-I5S).The cameras conform to the main current direction at 180-degree angles from each other.Each system includes two cameras: the first set for monitoring the ice speed and ice conditions, and the second set for monitoring the ice thickness and ice breaking.
The ice thickness and ice velocity data were obtained from the cameras based on the broken ice pieces that had turned over.The exposed section can be captured by a camera, as shown in Figure 2. The yellow part is the ice-broken cone of the OWT.By calibrating the camera with fixed focal lengths with an object of a known size ratio n is obtained, which relates the length of the object being measured, S, to the length on the image, s, as follows:  By calibrating the camera with fixed focal lengths with an object of a known size, the ratio n is obtained, which relates the length of the object being measured, S, to the pixel length on the image, s, as follows: Then, the expression of the measured ice thickness h and ice speed v is obtained as follows: where L is the ice breaking length, h and l are the pixel length of ice thickness and ice breaking length on the image, and t is the time that ice takes to pass the breaking length.
During the monitoring period, the level ice lasted for 15 days.It began on 25 January and ended on 8 February of 2018.Only the images with level ice were selected for analysis.The sea ice thickness and sea ice velocity were obtained by collecting multiple sets of data on a group of images and taking the mean value.The maximum level ice thickness measured was 14.9 cm, and the maximum level ice velocity measured was 74.8 cm/s (as shown in Figures 3 and 4).By calibrating the camera with fixed focal lengths with an object of a known size, the ratio n is obtained, which relates the length of the object being measured, S, to the pixel length on the image, s, as follows: Then, the expression of the measured ice thickness h and ice speed v is obtained as follows: L n l v t t where L is the ice breaking length, ' h and l are the pixel length of ice thickness and ice breaking length on the image, and t is the time that ice takes to pass the breaking length.
During the monitoring period, the level ice lasted for 15 days.It began on 25 January and ended on 8 February of 2018.Only the images with level ice were selected for analysis.The sea ice thickness and sea ice velocity were obtained by collecting multiple sets of data on a group of images and taking the mean value.The maximum level ice thickness measured was 14.9 cm, and the maximum level ice velocity measured was 74.8 cm/s (as shown in Figures 3 and 4).

Deterministic Ice Force Function
Based on long-term monitoring of the Bohai conical jacket platform, Yue et al. established an ice load model for a cone-shaped structure, namely, a deterministic ice force function [16].The ice force duration is simplified into the form of a right-angled trigonometric function of equal-period amplitudes (Figure 5).This load model assumes that sea ice has the same physical and mechanical properties in the loading section and simultaneously pulses the ice force in each cycle.The peak is used as the starting point of the cycle.In a very short loading time, the sea ice reaches its breaking strength, and the ice force then begins to gradually unload to zero.According to the conclusions obtained from the analysis of measured data of the Bohai oil and gas platform by Yue and Bi [16], as well

Deterministic Ice Force Function
Based on long-term monitoring of the Bohai conical jacket platform, Yue et al. established an ice load model for a cone-shaped structure, namely, a deterministic ice force function [16].The ice force duration is simplified into the form of a right-angled trigonometric function of equal-period amplitudes (Figure 5).This load model assumes that sea ice has the same physical and mechanical properties in the loading section and simultaneously Energies 2022, 15, 3357 5 of 20 pulses the ice force in each cycle.The peak is used as the starting point of the cycle.In a very short loading time, the sea ice reaches its breaking strength, and the ice force then begins to gradually unload to zero.According to the conclusions obtained from the analysis of measured data of the Bohai oil and gas platform by Yue and Bi [16], as well as Ji and Yue [21], the unloading period of the ice force is 1/3 of the full ice force cycle.

Deterministic Ice Force Function
Based on long-term monitoring of the Bohai conical jacket platform, Yue et al. established an ice load model for a cone-shaped structure, namely, a deterministic ice force function [16].The ice force duration is simplified into the form of a right-angled trigonometric function of equal-period amplitudes (Figure 5).This load model assumes that sea ice has the same physical and mechanical properties in the loading section and simultaneously pulses the ice force in each cycle.The peak is used as the starting point of the cycle.In a very short loading time, the sea ice reaches its breaking strength, and the ice force then begins to gradually unload to zero.According to the conclusions obtained from the analysis of measured data of the Bohai oil and gas platform by Yue and Bi [16], as well as Ji and Yue [21], the unloading period of the ice force is 1/3 of the full ice force cycle.According to experimental analysis, Li and Yue found that the bending fracture length of sea ice is mainly controlled by the sea ice thickness, and the two can be expressed as a linear relationship [22].
where L is the length of the sea ice break, h is the ice thickness and k is the ratio of the icebreak length to the ice thickness, which is referred to as the long-thickness ratio coefficient.Through an analysis of the ice-image data obtained from this field monitoring system, the ice thickness h, length of the sea ice break L and the ratio k were calculated according to the method above.The main data are shown in Table 2.According to experimental analysis, Li and Yue found that the bending fracture length of sea ice is mainly controlled by the sea ice thickness, and the two can be expressed as a linear relationship [22].
where L is the length of the sea ice break, h is the ice thickness and k is the ratio of the icebreak length to the ice thickness, which is referred to as the long-thickness ratio coefficient.Through an analysis of the ice-image data obtained from this field monitoring system, the ice thickness h, length of the sea ice break L and the ratio k were calculated according to the method above.The main data are shown in Table 2.

Sea Ice Breaking Characteristics of Conical Structure
Because the physical and mechanical properties of sea ice exhibit a certain degree of randomness, the fragmentation of sea ice is non-simultaneous on a wider cone structure.
Additionally, the greater the width of a structure, and the thinner the sea ice, the more irregular broken ice is.As shown in Figure 6, when the ratio of the diameter of the cone structure to the ice thickness (the ratio of diameter to thickness) increases, the size of the sea ice fracture decreases, and sea ice breaking occurs more non-simultaneously (photos of a and b were taken from JZ20-2MUQ in the Bohai Sea [17], and photos of c and d were taken from an OWT in the Zhuanghe sea area).Therefore, to more accurately simulate ice forces on conical OWTs, it is necessary to analyze the randomness of the ice load.From on-site image data processing of a conical OWT, sea ice breaking undergoes three stages during the bending failure process: sea ice from the contact cone to the stage of bending failure (ice loading stage), and the broken ice climbs along the cone due to pushing by the subsequent ice plates.The thrust of the subsequent ice plates on the broken ice becomes increasingly smaller (ice force unloading stage).At the same time, the surface on the cone is relatively smooth, and the friction is small.It can be considered that the ice force unloads to 0; the broken ice slides away from both sides of the cone, and the subsequent ice plate has not yet been in contact with the surface of the cone; the ice force is 0 during this period.The fracture characteristics of the sea ice in front of OWTs were the same as those of the sea ice in front of jacket structures in Bohai.
Based on the ice force measurement data of the pressure box of the Bohai conical jacket structure, Qu et al. proposed a model of the random ice force function which simplified the cone ice load into an isosceles triangle with random peaks and periodic distributions [19].In this function, the time from loading to unloading of the ice forces is 1/3 of the entire period (as shown in Figure 7).From on-site image data processing of a conical OWT, sea ice breaking undergoes three stages during the bending failure process: sea ice from the contact cone to the stage of bending failure (ice loading stage), and the broken ice climbs along the cone due to pushing by the subsequent ice plates.The thrust of the subsequent ice plates on the broken ice becomes increasingly smaller (ice force unloading stage).At the same time, the surface on the cone is relatively smooth, and the friction is small.It can be considered that the ice force unloads to 0; the broken ice slides away from both sides of the cone, and the subsequent ice plate has not yet been in contact with the surface of the cone; the ice force is 0 during this period.The fracture characteristics of the sea ice in front of OWTs were the same as those of the sea ice in front of jacket structures in Bohai.
Based on the ice force measurement data of the pressure box of the Bohai conical jacket structure, Qu et al. proposed a model of the random ice force function which simplified the cone ice load into an isosceles triangle with random peaks and periodic distributions [19].In this function, the time from loading to unloading of the ice forces is 1/3 of the entire period (as shown in Figure 7).jacket structure, Qu et al. proposed a model of the random ice force f plified the cone ice load into an isosceles triangle with random peaks butions [19].In this function, the time from loading to unloading of th the entire period (as shown in Figure 7).The ice force function throughout one cycle is defined by: The ice force function throughout one cycle is defined by: where T i is the i-th ice force cycle size and F 0i is the ice force amplitude in the cycle.On the basis of Equation ( 5), a random ice force function of any length of time can be extended: where N represents the number of ice cycles that need to be established.When i = 0, t 1 0 = 0 and when i > 1, When the ice force period is expanded to N, the ice force amplitude is F 0i for each period, and the ice force period T i constitutes two random series.
To obtain a random ice force function that can be applied, it is necessary to calculate the ice force period and ice force amplitude of the real ice force and determine the random distribution of the ice force period and amplitude sequence used in Equation (1).
Based on the measured data of the Bohai platform, the statistical characteristics of the ice force amplitude and ice force period of the cone ice force were analyzed.It was found that both the ice force period and the ice force amplitude obeyed normal distributions, as shown in Figures 8 and 9.
where Ti is the i-th ice force cycle size and F0i is the ice force amplitude in the cycle.On the basis of Equation ( 5), a random ice force function of any length of time can be extended: where N represents the number of ice cycles that need to be established.When i = 0 When the ice force period is expanded to N, the ice force amplitude is F0i for each period, and the ice force period Ti constitutes two random series.To obtain a random ice force function that can be applied, it is necessary to calculate the ice force period and ice force amplitude of the real ice force and determine the random distribution of the ice force period and amplitude sequence used in Equation (1).
Based on the measured data of the Bohai platform, the statistical characteristics of the ice force amplitude and ice force period of the cone ice force were analyzed.It was found that both the ice force period and the ice force amplitude obeyed normal distributions, as shown in Figures 8 and 9.
Based on the study by Qu et al. [19], the standard deviation of ice forces is u to 70% of the average value during the ice cycle.An average of 0.5 can be tak value of the standard deviation in the application. .T = At the same time, the coefficients of the variations in ice force amplitude domly distributed between 0.2 and 0.6, and the average value is approximately application, the standard deviation of the ice force amplitude was considered times the average value of the ice force, as shown in Equation (10):

Average Period of Ice Force
To predict the ice load model suitable for conical OWT, the ice breakin caused by different conical structures were analyzed.The ice force period, T, i of ice breaking length L to the ice speed v, which is expressed as follows: For the ice cycle sequence T i ∼ N(T, σ 2 T ), the probability density function can be expressed as: For the ice force amplitude F i ∼ N(F, σ 2 F ), the probability density function can be expressed as: Based on the study by Qu et al. [19], the standard deviation of ice forces is usually 40 to 70% of the average value during the ice cycle.An average of 0.5 can be taken as the value of the standard deviation in the application.
At the same time, the coefficients of the variations in ice force amplitudes are randomly distributed between 0.2 and 0.6, and the average value is approximately 0.4.In this application, the standard deviation of the ice force amplitude was considered to be 0.4 times the average value of the ice force, as shown in Equation (10):

Average Period of Ice Force
To predict the ice load model suitable for conical OWT, the ice breaking features caused by different conical structures were analyzed.The ice force period, T, is the ratio of ice breaking length L to the ice speed v, which is expressed as follows: Based on field-measured and laboratory test data (as shown in Table 3), it can be seen that the aspect ratios of structures with different scales were different.Moreover, these ratios generally increased with increased cone diameters (Figure 10).This is because as the diameter of the cone increased, the linear width of the ice plate thickness increased.At the same time, the length of sea ice breaks and the dispersion of the length of the fracture increased (ice thicknesses were from 5 cm to 40 cm).At the same time, the length of sea ice breaks and the dispersion of the leng ture increased (ice thicknesses were from 5 cm to 40 cm).Since the cone width influences the length-to-thickness ratio coefficien fitting, Xu introduced a dimensionless parameter, the width-to-thickness the analysis of ice cones (the ratio of diameter D of the cone at the sea ice to of sea ice h) [17]; the coefficient q = D/h; therefore, this coefficient was used variation law of the aspect ratio.Through video data analysis, it was found ice acted on a wind-resistant anti-ice cone, when the ice thickness was less sea ice breakage length was more dispersed; multiple fractures occurred an into small ice cubes, or splitting damage occurred.Therefore, the measurem < 160 were selected for analysis.
Through an analysis of the data, it was found that coefficient k grew with the aspect ratio when the width-to-thickness ratio q was smaller, an thickness ratio grew more slowly when the aspect ratio q was larger ( Since the cone width influences the length-to-thickness ratio coefficient obtained by fitting, Xu introduced a dimensionless parameter, the width-to-thickness ratio (D/h), in the analysis of ice cones (the ratio of diameter D of the cone at the sea ice to the thickness of sea ice h) [17]; the coefficient q = D/h; therefore, this coefficient was used to analyze the variation law of the aspect ratio.Through video data analysis, it was found that when sea ice acted on a wind-resistant anti-ice cone, when the ice thickness was less than 5 cm, the sea ice breakage length was more dispersed; multiple fractures occurred and were broken into small ice cubes, or splitting damage occurred.Therefore, the measurement data for q < 160 were selected for analysis.
Through an analysis of the data, it was found that coefficient k grew more rapidly with the aspect ratio when the width-to-thickness ratio q was smaller, and the width-thickness ratio grew more slowly when the aspect ratio q was larger (Figure 11).Therefore, the logarithmic form of the equation was used for analysis: Energies From the analysis of the measured data, it was found that the length-thickness ratio coefficient, k, was greater than 4, and the length-thickness ratio coefficient, k, on the cone, which is recommended in ISO19906 [23], ranged from 4~10.Therefore, we used a = 4 in the formula and obtained the relationship between the length-thickness ratio coefficient, k, and the width-thickness ratio coefficient, q, from the measured data; thus, the relationship between the average period of sea ice breaking and ice thickness and ice speed can be obtained: Based on the ice thickness and ice speed data obtained from field measurements, the ice force cycle was calculated by using the above-described average ice force period relationship, and the measured ice force cycle was compared and analyzed (Figure 12).The overall calculated average period was slightly longer than the measured ice force period.Sea ice inhomogeneity and secondary fractures cause the extracted ice force cycle to be shorter.The measured ice force cycle was based on pressure box data analysis, and the ice thickness and ice velocity were obtained from ice-image cameras.During the climbing process of sea ice interacting on the conical structure, secondary fractures may occur.The period of secondary fracture is less than the primary failure period of sea ice.The period measured by the ice load measured panel includes the period of secondary fracture.Then, compared with the structural response caused by the primary fracture of sea ice, the structural response caused by the secondary fracture of sea ice climbing is small, and the fatigue damage to the structure is also small.The above expression was used to calculate the average ice force cycle to meet engineering needs.From the analysis of the measured data, it was found that the length-thickness ratio coefficient, k, was greater than 4, and the length-thickness ratio coefficient, k, on the cone, which is recommended in ISO19906 [23], ranged from 4~10.Therefore, we used a = 4 in the formula and obtained the relationship between the length-thickness ratio coefficient, k, and the width-thickness ratio coefficient, q, from the measured data; thus, the relationship between the average period of sea ice breaking and ice thickness and ice speed can be obtained: Based on the ice thickness and ice speed data obtained from field measurements, the ice force cycle was calculated by using the above-described average ice force period relationship, and the measured ice force cycle was compared and analyzed (Figure 12).The overall calculated average period was slightly longer than the measured ice force period.Sea ice inhomogeneity and secondary fractures cause the extracted ice force cycle to be shorter.The measured ice force cycle was based on pressure box data analysis, and the ice thickness and ice velocity were obtained from ice-image cameras.During the climbing process of sea ice interacting on the conical structure, secondary fractures may occur.The period of secondary fracture is less than the primary failure period of sea ice.The period measured by the ice load measured panel includes the period of secondary fracture.Then, compared with the structural response caused by the primary fracture of sea ice, the structural response caused by the secondary fracture of sea ice climbing is small, and the fatigue damage to the structure is also small.The above expression was used to calculate the average ice force cycle to meet engineering needs.From the analysis of the measured data, it was found that the length-thickness ratio coefficient, k, was greater than 4, and the length-thickness ratio coefficient, k, on the cone, which is recommended in ISO19906 [23], ranged from 4~10.Therefore, we used a = 4 in the formula and obtained the relationship between the length-thickness ratio coefficient, k, and the width-thickness ratio coefficient, q, from the measured data; thus, the relationship between the average period of sea ice breaking and ice thickness and ice speed can be obtained: Based on the ice thickness and ice speed data obtained from field measurements, the ice force cycle was calculated by using the above-described average ice force period relationship, and the measured ice force cycle was compared and analyzed (Figure 12).The overall calculated average period was slightly longer than the measured ice force period.Sea ice inhomogeneity and secondary fractures cause the extracted ice force cycle to be shorter.The measured ice force cycle was based on pressure box data analysis, and the ice thickness and ice velocity were obtained from ice-image cameras.During the climbing process of sea ice interacting on the conical structure, secondary fractures may occur.The period of secondary fracture is less than the primary failure period of sea ice.The period measured by the ice load measured panel includes the period of secondary fracture.Then, compared with the structural response caused by the primary fracture of sea ice, the structural response caused by the secondary fracture of sea ice climbing is small, and the fatigue damage to the structure is also small.The above expression was used to calculate the average ice force cycle to meet engineering needs.

Ice Force Amplitude
The average ice force, F, of sea ice is related to the extreme ice force.Qu et al. obtained F = F max /1.8 from a statistical analysis of measured data.F max is the extreme static force on the ice breaking cone [19].Xu [17], Li and Yue [22], as well as Yu [24] and other scholars, studied the extreme static ice load model on conical structures and analyzed the applicability of the existing static ice load model on the conical platform in Bohai.By contrast, the empirical formula for the total horizontal component of ice breaking derived by Hirayama and Obara (1986) was applied to the narrower cone structure extreme ice load analysis [12]: where σ f is the sea ice bending strength, the average bending strength is 750 kPa, h is the sea ice thickness, D is the diameter of the cone structure at the level of ice action and B is an empirical coefficient that takes 3.7 in this paper.L c is the characteristic length of broken ice, ; E is the elastic modulus, and its recommended value is 1.5 GPa; ρ is the seawater density, 1.025 kg•m −3 ; and g is the gravitational acceleration, 9.8 m•s −2 .

Dynamic Ice Load Analysis of Conical OWT
To further analyze the rationality of the ice load model, the OWT measured in the Zhuanghe sea area was simulated numerically.The structure is a monopile foundation OWT, and a forward and backward ice-breaking cone was installed at the waterline.Its basic properties are shown in Table 1.The results of the measurement and simulation were compared to clarify the applicability of the ice load model to OWTs.

Characteristics of the OWT
According to the properties of the OWT structure in Table 1, a simplified OWT model was established based on the ANSYS finite element software.The tower and foundation were simulated by shell181 elements, and the turbine was simulated by an MASS21 element.By adjusting the model grid, the mesh size was 0.2 m, which can meet the requirements of calculation accuracy and efficiency.The equivalent embedded constraint was used to simulate the constraint of the foundation below mud surface of the OWT.The equivalent pile length of the FEM model is five times that of the pile diameter.The FEM model and modal analysis results of the OWT are shown in Figure 13.
force on the ice breaking cone [19].Xu [17], Li and Yue [22], as well as Yu [24] scholars, studied the extreme static ice load model on conical structures and ana applicability of the existing static ice load model on the conical platform in Boha trast, the empirical formula for the total horizontal component of ice breaking d Hirayama and Obara (1986) was applied to the narrower cone structure extrem analysis [12]: where  is the sea ice bending strength, the average bending strength is 750 kP sea ice thickness, D is the diameter of the cone structure at the level of ice actio an empirical coefficient that takes 3.7 in this paper.Lc is the characteristic length ice, is the seawater density, 1.025 kg•m −3 ; and g is the gravitational acceleration, 9.8

Dynamic Ice Load Analysis of Conical OWT
To further analyze the rationality of the ice load model, the OWT measu Zhuanghe sea area was simulated numerically.The structure is a monopile fo OWT, and a forward and backward ice-breaking cone was installed at the wat basic properties are shown in Table 1.The results of the measurement and s were compared to clarify the applicability of the ice load model to OWTs.

Characteristics of the OWT
According to the properties of the OWT structure in Table 1, a simplified OW was established based on the ANSYS finite element software.The tower and fo were simulated by shell181 elements, and the turbine was simulated by an MA ment.By adjusting the model grid, the mesh size was 0.2 m, which can meet th ments of calculation accuracy and efficiency.The equivalent embedded const used to simulate the constraint of the foundation below mud surface of the O equivalent pile length of the FEM model is five times that of the pile diameter.model and modal analysis results of the OWT are shown in Figure 13.The first-order mode of the OWT swings in the directions of the X-axis and Y-axis, and the natural frequency is approximately 0.265 Hz.The calculated first frequency was close to the actual frequency (0.26 Hz), and the second frequency was approximately 1.02 Hz.
The dynamic characteristics of a structure are determined by mass, stiffness and damping.However, it is difficult to determine the real damping characteristics of a complex structure.The equivalent viscous damping ratio ζ, with the same attenuation rate under free vibration, is typically used to represent structural damping.According to the measured attenuation curve of OWT free vibrations, the damping ratio of the structure calculated by the logarithmic reduction method was approximately 0.02.The Rayleigh damping model is often used in vibration analysis in ANSYS.It was assumed that the damping matrix of the structure was a linear combination of the mass matrix and stiffness matrix (C = αM + βK).In this expression, α is the mass damping coefficient, and β is the stiffness damping coefficient.The i-th modal damping ratio ζ i can be expressed as: Because the first-order vibration characteristics of OWTs are mainly considered under the action of sea ice, the first two natural frequencies were mainly analyzed.According to the analysis of the measured data, the damping ratio of the structure was 0.02.Then, the parameters α and β were calculated.

Verification of the Random Ice Load Model of the OWT
To measure the coupling effect between sea ice and other environmental loads, the structural-response monitoring system consisted of four sets of acceleration tilt sensors which were evenly distributed inside the tower from the top.The X-axes of the sensors were parallel to the main ice direction.For the monitoring of sea ice conditions and ice breaking, image monitoring methods using two sets of ice-imaging cameras were adopted.Through the analyses of the monitoring data, we extracted data from a day when the wind speed was low (lower than moderate breeze) and ice conditions were serious.At 13:10 on that day, the average ice thickness was 0.11 m, and the ice velocity was 0.51 m•s −1 (Figure 14).
damping.However, it is difficult to determine the real damping characterist plex structure.The equivalent viscous damping ratio ζ, with the same atte under free vibration, is typically used to represent structural damping.Acco measured attenuation curve of OWT free vibrations, the damping ratio of calculated by the logarithmic reduction method was approximately 0.02.T damping model is often used in vibration analysis in ANSYS.It was assum damping matrix of the structure was a linear combination of the mass matrix matrix ( C M K α β

= +
).In this expression, α is the mass damping coeffici the stiffness damping coefficient.The i-th modal damping ratio ζi can be exp Because the first-order vibration characteristics of OWTs are mainly co der the action of sea ice, the first two natural frequencies were mainly analyze to the analysis of the measured data, the damping ratio of the structure wa the parameters α and β were calculated.

Verification of the Random Ice Load Model of the OWT
To measure the coupling effect between sea ice and other environment structural-response monitoring system consisted of four sets of acceleration which were evenly distributed inside the tower from the top.The X-axes o were parallel to the main ice direction.For the monitoring of sea ice condit breaking, image monitoring methods using two sets of ice-imaging ca adopted.Through the analyses of the monitoring data, we extracted data from the wind speed was low (lower than moderate breeze) and ice conditions w At 13:10 on that day, the average ice thickness was 0.11 m, and the ice veloc m•s −1 (Figure 14).From the above analysis of ice periods and ice amplitudes, the random conical OWT can be established.Based on the Monte Carlo method, a time h lation of random ice forces on a conical OWT using the MATLAB software was conducted; the ice force time history data were applied to the ANSYS OWT to analyze the structural response.An acceleration sensor was located From the above analysis of ice periods and ice amplitudes, the random ice force of a conical OWT can be established.Based on the Monte Carlo method, a time history simulation of random ice forces on a conical OWT using the MATLAB software (Figure 15b) was conducted; the ice force time history data were applied to the ANSYS FEM of the OWT to analyze the structural response.An acceleration sensor was located in the tower and 14 m from the top of the foundation (Figure 15a).The acceleration response data of the structure corresponding to the sensor placement point in the FEM were extracted.
was 26.6%.By comparing the acceleration data, the maximum acceleration difference between the simulated value and the measured value was small, and the standard deviation had some errors.The average period of the simulated values was 1.52 s, and the average period of the measured values was 1.88 s.The error between them was 19.1%.The standard deviation error between the simulated value and the measured value was 2.8%.Because the influence of wind loads was not considered in the numerical simulation, the measured values were influenced by the wind load.Therefore, there were some errors between the amplitudes and the average periods in the measured and simulated values.From the comparison between the measured acceleration and the simulated acceleration time history curves in Figure 15c, it was found that the measured values of the acceleration amplitudes were slightly larger than the simulated values.As shown in Table 4, the maximum measured acceleration was 11.3 mm•s −2 , and that from the simulation was 10.98 mm•s −2 .The error between the two was 2.8%.In comparing the acceleration standard deviation of the measurement and the simulation, it is found that the error between them was 26.6%.By comparing the acceleration data, the maximum acceleration difference between the simulated value and the measured value was small, and the standard deviation had some errors.The average period of the simulated values was 1.52 s, and the average period of the measured values was 1.88 s.The error between them was 19.1%.The standard deviation error between the simulated value and the measured value was 2.8%.Because the influence of wind loads was not considered in the numerical simulation, the measured values were influenced by the wind load.Therefore, there were some errors between the amplitudes and the average periods in the measured and simulated values.
To further analyze the correlation between the measured and the simulated values and clarify the corresponding relationship between them in the frequency domain, the acceleration data were transformed with a Fourier transform.Fourier transforms can establish the conversion relationship between time-domain signals and frequency-domain signals and can reflect the proportion of signals at different frequencies.The fast Fourier transform (FFT) is a fast algorithm of the Fourier transform.The relationship between the amplitude and the frequency components of the acceleration data was obtained by FFT processing of the simulated and measured data.The comparative analysis results are shown in Figure 16.To further analyze the correlation between the measured and the simulate and clarify the corresponding relationship between them in the frequency dom acceleration data were transformed with a Fourier transform.Fourier transform tablish the conversion relationship between time-domain signals and frequency signals and can reflect the proportion of signals at different frequencies.The fast transform (FFT) is a fast algorithm of the Fourier transform.The relationship betw amplitude and the frequency components of the acceleration data was obtained processing of the simulated and measured data.The comparative analysis res shown in Figure 16. Figure 16 shows that the simulated values and the measured data had simila tude-frequency relationship curves and that the amplitude responses were conc at approximately 0.25 Hz and 1.0 Hz.By comparing the frequencies correspondin first three highest peaks (as shown in Table 5), the first peak frequency of simula 0.251 Hz, and that of measurement was 0.262 Hz.The error between the measu the simulated value was 4.2%.In comparing the next two peak frequencies of th lated value and the measured value, it was found that the errors between them w and 8.8%.The errors were relatively small.The amplitude responses of the num simulated accelerations were similar at both frequencies.The amplitude respon measured accelerations was mainly concentrated at 0.25 Hz, and there was a corr ing peak near 1.0 Hz.Different from the measured data, which were affected environmental conditions such as winds and currents, the numerical simulation o sidered the influence of sea ice loads; there were thus some differences in the am responses between the two, but the overall response trends were the same.
Based on the above analysis, the formulas for the average periodic ice forces random ice forces are suitable for conical OWT structures.Figure 16 shows that the simulated values and the measured data had similar amplitudefrequency relationship curves and that the amplitude responses were concentrated at approximately 0.25 Hz and 1.0 Hz.By comparing the frequencies corresponding to the first three highest peaks (as shown in Table 5), the first peak frequency of simulation was 0.251 Hz, and that of measurement was 0.262 Hz.The error between the measured and the simulated value was 4.2%.In comparing the next two peak frequencies of the simulated value and the measured value, it was found that the errors between them were 8.2% and 8.8%.The errors were relatively small.The amplitude responses of the numerically simulated accelerations were similar at both frequencies.The amplitude response of the measured accelerations was mainly concentrated at 0.25 Hz, and there was a corresponding peak near 1.0 Hz.Different from the measured data, which were affected by other environmental conditions such as winds and currents, the numerical simulation only considered the influence of sea ice loads; there were thus some differences in the amplitude responses between the two, but the overall response trends were the same.Based on the above analysis, the formulas for the average periodic ice forces and for random ice forces are suitable for conical OWT structures.

Ice-Induced Vibration of the OWT
The ultimate ice thickness in the wind farm was 0.32 m (50-year return period).Due to the lack of long-term distribution data for ice velocities, the ice speed data measured in that year were used as the analysis data for ice velocities.The measured maximum ice velocity was 0.75 m•s −1 .The ice load may induce resonance when the ice breaking period is close to the natural period of the structure where the natural period is 3.77 s.The ice velocity was 0.63 m•s −1 , which was calculated using Equations ( 13) and ( 14) in Section 4, when the ice thickness was 0.32 m.The numerical simulation of the structural response was based on the two ice conditions (as shown in Table 6).The numerical simulation curves in Figures 17 and 18 show that there was no resonance in the structural response curve at the hub when the average ice force period was equal to the first natural period of the structure, and the acceleration curve characteristics and amplitude were similar to those under ultimate ice conditions.On the one hand, because the ice velocities of the two were close and the ice velocity was relatively high, the ice force periods were relatively close, and the extreme ice forces of the two were the same; thus, the vibration responses of the structure were relatively close.On the other hand, when the ice force period was close to the structural period, the resonance amplification phenomenon of the structure response was not induced, and the vibration of the structure still presented randomness.Due to the nonuniformity of natural sea ice, the bending failure of the ice exhibited randomness.The ultimate ice thickness in the wind farm was 0.32 m (50-year return perio to the lack of long-term distribution data for ice velocities, the ice speed data mea that year were used as the analysis data for ice velocities.The measured maxim velocity was 0.75 m•s −1 .The ice load may induce resonance when the ice breaking is close to the natural period of the structure where the natural period is 3.77 s. velocity was 0.63 m•s −1 , which was calculated using Equations ( 13) and ( 14) in S when the ice thickness was 0.32 m.The numerical simulation of the structural r was based on the two ice conditions (as shown in Table 6).The numerical simulation curves in Figures 17 and 18 show that there was nance in the structural response curve at the hub when the average ice force per equal to the first natural period of the structure, and the acceleration curve charac and amplitude were similar to those under ultimate ice conditions.On the one h cause the ice velocities of the two were close and the ice velocity was relatively h ice force periods were relatively close, and the extreme ice forces of the two were th thus, the vibration responses of the structure were relatively close.On the othe when the ice force period was close to the structural period, the resonance ampl phenomenon of the structure response was not induced, and the vibration of the s still presented randomness.Due to the nonuniformity of natural sea ice, the bend ure of the ice exhibited randomness.The numerical simulation curves in Figures 19 and 20 show that the acceleration response was random at the top of the foundation (ToF) and that the acceleration curve characteristics and amplitude were similar to those under ultimate ice conditions.The reason was the same as described above.The numerical simulation curves in Figures 19 and 20 show that the acceleratio sponse was random at the top of the foundation (ToF) and that the acceleration c characteristics and amplitude were similar to those under ultimate ice conditions.reason was the same as described above.By comparing Figures 17-20, it can be seen that the response at the ToF and at the hub for the two ice conditions exhibited the same trend.The difference was that the acceleration response frequency at the ToF was higher than that at the hub, and the acceleration response at the ToF was similar to that at the hub.However, the displacement response at the ToF was far less than that at the hub (as shown in Table 7).The vibration acceleration at the ToF should be considered for OWTs in cold regions.After transient analysis, the responses with large von Mises stress and the largest amplitude changes were defined as hot spot stress.It can be seen from the hot spot stress positions (Figures 21 and 22) and stress curves (Figures 23 and 24) that the hot spot stress position was the same under resonant ice conditions (RIC) and ultimate ice conditions (UIC); the maximum stress was 116.4 MPa under RIC, and the maximum stress was 116.1 MPa under UIC.The stress results were calculated using the rain flow counting method.These results showed that the maximum stress amplitude was 37.7 MPa and the corresponding average stress was 97.5 MPa under RIC; the maximum stress amplitude was 35.6 MPa and the corresponding average stress was 97.5 MPa under UIC (shown in Table 8).
(S a /S −1 ) + (S m /S u ) = 1 (17) Energies 2022, 15, x FOR PEER REVIEW 18 of 21 response at the ToF was similar to that at the hub.However, the displacement response at the ToF was far less than that at the hub (as shown in Table 7).The vibration acceleration at the ToF should be considered for OWTs in cold regions.After transient analysis, the responses with large von Mises stress and the largest amplitude changes were defined as hot spot stress.It can be seen from the hot spot stress positions (Figures 21 and 22) and stress curves (Figures 23 and 24) that the hot spot stress position was the same under resonant ice conditions (RIC) and ultimate ice conditions (UIC); the maximum stress was 116.4 MPa under RIC, and the maximum stress was 116.1 MPa under UIC.The stress results were calculated using the rain flow counting method.These results showed that the maximum stress amplitude was 37.7 MPa and the corresponding average stress was 97.5 MPa under RIC; the maximum stress amplitude was 35.6 MPa and the corresponding average stress was 97.5 MPa under UIC (shown in Table 8).Energies 2022, 15, x FOR PEER REVIEW 18 of 21 response at the ToF was similar to that at the hub.However, the displacement response at the ToF was far less than that at the hub (as shown in Table 7).The vibration acceleration at the ToF should be considered for OWTs in cold regions.After transient analysis, the responses with large von Mises stress and the largest amplitude changes were defined as hot spot stress.It can be seen from the hot spot stress positions (Figures 21 and 22) and stress curves (Figures 23 and 24) that the hot spot stress position was the same under resonant ice conditions (RIC) and ultimate ice conditions (UIC); the maximum stress was 116.4 MPa under RIC, and the maximum stress was 116.1 MPa under UIC.The stress results were calculated using the rain flow counting method.These results showed that the maximum stress amplitude was 37.7 MPa and the corresponding average stress was 97.5 MPa under RIC; the maximum stress amplitude was 35.6 MPa and the corresponding average stress was 97.5 MPa under UIC (shown in Table 8).

Figure 1 .
Figure 1.FMS implementation on an OWT foundation: (a) FMS implementation on an OWT foundation; (b) acceleration tilt sensors; (c) first set of ice-image cameras; (d) second set of ice-image cameras.

Figure 1 .
Figure 1.FMS implementation on an OWT foundation: (a) FMS implementation on an OWT foundation; (b) acceleration tilt sensors; (c) first set of ice-image cameras; (d) second set of ice-image cameras.

Energies 2022 ,Figure 2 .
Figure 2. Image monitoring methods of sea ice factors.

Figure 2 .
Figure 2. Image monitoring methods of sea ice factors.

Figure 2 .
Figure 2. Image monitoring methods of sea ice factors.

Figure 3 .
Figure 3.The maximum level ice thickness measured per day in January and February 2018.

Figure 3 .
Figure 3.The maximum level ice thickness measured per day in January and February 2018.Energies 2022, 15, x FOR PEER REVIEW 5 of 21

Figure 4 .
Figure 4.The maximum level ice velocity measured per day in January and February 2018.

Figure 4 .
Figure 4.The maximum level ice velocity measured per day in January and February 2018.

Figure 5 .
Figure 5. Simplified deterministic ice force model of a conical structure.

Figure 5 .
Figure 5. Simplified deterministic ice force model of a conical structure.

Figure 7 .
Figure 7. Simplified random ice force model on the cone structure.

Figure 7 .
Figure 7. Simplified random ice force model on the cone structure.
σ , the probability density funct expressed as:

Figure 9 .
Figure 9. Normal distribution characteristics of flat-ice force amplitudes.

Figure 10 .
Figure 10.Length-thickness ratios, k, of sea ice bending failures at different scales.

Figure 10 .
Figure 10.Length-thickness ratios, k, of sea ice bending failures at different scales.

Figure 12 .
Figure 12.Comparison between the calculated ice force periods and the measured ice force periods.

;
E is the elastic modulus, and its recommended value is 1

Figure 13 .
Figure 13.FEM model and modal analysis results for the OWT.

Figure 14 .
Figure 14.Bending failure phenomenon of ice on a conical OWT.

Figure 14 .
Figure 14.Bending failure phenomenon of ice on a conical OWT.

Figure 15 .
Figure 15.(a) Location point of the acceleration sensor; (b) time history curve of random ice forces under typical ice conditions; (c) simulated acceleration versus measured acceleration.

Figure 15 .
Figure 15.(a) Location point of the acceleration sensor; (b) time history curve of random ice forces under typical ice conditions; (c) simulated acceleration versus measured acceleration.

Figure 16 .
Figure 16.The measured data versus analysis result of acceleration after FFT transformati

Figure 16 .
Figure 16.The measured data versus analysis result of acceleration after FFT transformation.

Figure 17 .
Figure 17.Response accelerations of the hub under resonant ice conditions (ice thickness m and ice velocity was 0.63 m•s −1 ).

Figure 17 .
Figure 17.Response accelerations of the hub under resonant ice conditions (ice thickness was 0.32 m and ice velocity was 0.63 m•s −1 ).

Figure 18 .
Figure 18.Response accelerations of the hub under ultimate ice conditions (ice thickness was m and ice velocity was 0.75 m•s −1 ).

Figure 19 .
Figure 19.Response accelerations at the ToF under resonant ice conditions (ice thickness was m and ice velocity was 0.63 m•s −1 ).

Figure 20 .
Figure 20.Response accelerations at the ToF under ultimate ice conditions (ice thickness was 0. and ice velocity was 0.75 m•s −1 ).

Figure 21 .
Figure 21.Hot spot under RIC (ice thickness was 0.32 m and ice velocity was 0.63 m•s −1 ).

Figure 22 .
Figure 22.Hot spot under UIC (ice thickness was 0.32 m and ice velocity was 0.75 m•s −1 ).

Figure 21 .
Figure 21.Hot spot under RIC (ice thickness was 0.32 m and ice velocity was 0.63 m•s −1 ).

Figure 21 .
Figure 21.Hot spot under RIC (ice thickness was 0.32 m and ice velocity was 0.63 m•s −1 ).

Figure 22 .
Figure 22.Hot spot under UIC (ice thickness was 0.32 m and ice velocity was 0.75 m•s −1 ).

Figure 22 .
Figure 22.Hot spot under UIC (ice thickness was 0.32 m and ice velocity was 0.75 m•s −1 ).

Table 1 .
Basic properties of the OWT structure.

Table 2 .
The sea ice information measured from conical OWT.

Table 3 .
Structural characteristics of different anti-ice cones and their length-thickness ratios, k.

Table 3 .
Structural characteristics of different anti-ice cones and their length-thickn FF-first-order frequency; WS-waterline stiffness; and MCD-maximum cone length-thickness ratio.

Table 4 .
Acceleration and period comparison between simulation and measurement.

Table 4 .
Acceleration and period comparison between simulation and measurement.std -standard deviation of acceleration; a max -maximum acceleration; T ave -average period; T std -standard deviation of period. a

Table 5 .
Comparison of peak frequency between simulated value and measured value.

Table 5 .
Comparison of peak frequency between simulated value and measured value.

Table 6 .
The design ice load for the simulation.

Table 6 .
The design ice load for the simulation.

Table 7 .
The responses of the OWT under UIC and RIC.-ice conditions; RIC-resonant ice conditions; UIC-ultimate ice conditions; MD-maximum displacement; and MA-maximum acceleration. IC

Table 7 .
The responses of the OWT under UIC and RIC.-ice conditions; RIC-resonant ice conditions; UIC-ultimate ice conditions; MD-maximum displacement; and MA-maximum acceleration. IC

Table 7 .
The responses of the OWT under UIC and RIC.-ice conditions; RIC-resonant ice conditions; UIC-ultimate ice conditions; MD-maximum displacement; and MA-maximum acceleration. IC