NUMERICAL SIMULATION OF ICE MILLING LOADS ON PROPELLER BLADE WITH COHESIVE ELEMENT METHOD

In ice-infested waters, propellers of a polar ship are likely to be exposed to ice loads in different scenarios. Propeller milling with ice is one of the most dangerous cases for icepropeller interaction. In this study, we try to simulate dynamic milling process of icepropeller and reproduce resulting physical phenomena. Cohesive element method is used to model ice in the simulation. To simulate material properties of ice, an elastoplastic softening constitutive law is developed. Both crushing and fracture failures are included in the icepropeller milling process. The ice loads in 6 Dofs acting on blades of a propeller are calculated in time domain. The average and standard deviations of simulated dominant ice loads are compared with those from model test. A good agreement is achieved. By varying propeller rotation speed, advance velocity and cutting depth on ice block, the sensitivity study has been carried out. The results show that dominant ice loads are affected much by the three parameters. It is shown that decreasing rotation speed, or increasing advance velocity and cutting depth may lead to higher ice loads. Care should be taken to avoid over-loading on propeller when operating in ice for polar ship.


Introduction
As global warming and ice shrinking in the polar waters, increasing voyages are expected to come for polar ships.As an essential equipment, propeller installed on polar ships may be subjected to ice loads in different ways since sea ice exists widely.It is not enough to check only the hydrodynamic loads on the propeller as studied by Ghassabzadeh et al. [1], but it is also necessary to estimate the ice loads on the propeller.The impact load could result in damage to structure of the propeller and endanger shipping safety.Therefore, it is worth to investigate the ice-propeller interaction.
Feng Wang, Zao-Jian Zou, Li Zhou Numerical simulation of ice milling loads on propeller Yang Wang, Hao Yu, Haihua Zhang blade with cohesive element method 110 Generally speaking, ice-propeller interaction could be divided into three scenarios.The first one is called ice blockage, where ice block keeps a certain distance from propeller.The existing of ice has an influence on stream field close to propeller.The propulsion efficiency would decrease accordingly.The second one is direct collision of ice pieces with propeller.In this case, the ice pieces are often small enough and taken as rigid bodies so that they are not able to break during collision.The last one is that propeller mills with large ice block.The resulting ice forces are usually an order of magnitude larger than the forces occurred in the first two cases.Consequently, it is meaningful to do research on propeller-ice milling process.
Since the 1990s, full scale and model scale tests for propeller-ice interaction have been performed and used for prediction of ice loads.As a part of Joint Research Project Arrangement between Canada and Finland, many researchers carried out full-scale measurements by using simplified blades and sea ice [2,3,4,5,6].Considering the cost problem, scale-model tests were also carried out with refrigerated model ice [7,8,9,10,11].As for numerical method, only a few researchers have studied ice-propeller milling in complex conditions.Wang et al. [12,13] introduced a numerical prediction method to calculate the ice milling action.Ye et al. [14] and Wang et al. [15] investigated propeller-ice contact using peridynamics method, and focused on analyzing the effect of factors influencing the contact load.Khan et al. [16] developed a numerical tool to assess the effects of propellerice interaction on the loads acting on the propellers.The tool was then calibrated based on the results from model tests.
In this paper, a numerical model is built up by combining cohesive element method with an elastoplastic softening constitutive law to simulate dynamic propeller milling with ice process.In order to validate the accuracy of the present model, a verification study is accomplished by comparing simulated result with the measurement from model test carried out at the ice tank of Institute for Ocean Technology (IOT) by Wang et al. [11].The main factors that may affect the milling process, including rotation speed of screw, advance velocity and cutting depth, are considered in the sensitivity analysis.Finally, some conclusions are drawn based on simulation results.

Numerical models
When the propeller blade contacts with the ice, the ice failure includes crushing, fracture and crack.Wang et al. [11] have found that ice cracking arises in the interaction between ice and propeller in ice milling tests, and Khan et al. [16] observed crushing happened in propeller milling ice tests when the blade edge cut the ice block.Such crushing and fracture failures are also found in the interaction between ship and level ice [17,18,19,20].These dis-continuous problems are difficult to be solved by simply using traditional FEM method that is based on the continuous assumption.Thus, cohesive element method (CEM) is introduced to solve the occurrence and propagation of ice fracture and crack during the icestructure interaction.
Cohesive element method (CEM) is an extension of cohesive zone model (CZM) in FEM simulations.CZM was firstly proposed by Hillerborg et al. [21] to model the local crack propagation of an unreinforced concrete beam.Due to the good performance of CEM in modeling the initiation and propagation of crack, many researchers have introduced CEM into ice mechanics to simulate the fracture of ice material [22,23,24,25].
In the framework of CEM, the ice elements are divided into two kinds of elements: bulk elements and cohesive elements.As shown in Fig. 1, the cohesive elements are inserted into every internal faces between the neighboring finite bulk elements.The bulk elements and the cohesive elements share nodes on the contact surfaces, thus deformation and stress could transmit directly between them.The cohesive elements respond to the external tensile force, and cohesive surface will separate based on a traction-separation law.After reaching failure criteria, the cohesive surfaces are completely separated and the cohesive elements will be damaged.Thus, the crack is generated and propagates along the path of cohesive surfaces.On the other hand, the bulk elements also respond to the external force based on its own constitutive model.

Fig. 1. Illustration of bulk elements and cohesive elements
The most commonly used traction-separation law in CEM is linear softening model, which is firstly proposed by Mi et al. [26].The function of linear softening model is listed as follows, where T and S are the cohesive forces in normal direction (mode I fracture) and tangential direction (mode II fracture), respectively; δn and δτ are crack separations in normal and tangential direction; T0 and S0 are the maximum stress in normal and tangential direction which are called fracture stress, while δn 1 and δτ 1 are the corresponding crack separations.δn f and δτ f are the maximum crack separations in normal and tangential direction.Fig. 2 shows the relationship between normal separation and cohesive force for mode I and II fracture.Taking mode I fracture as an example, T at the crack tip zone firstly linearly increases with the growth of δn under the external loads.When T reaches maximum stress T0, the crack begins to arise from this time.With the continuous development of δn, the material stiffness presents linear softening with the gradually expansion of the crack.When T decreases to zero, δn reaches maximum crack separations δn f , which means the cohesive element is damaged and  Besides the fracture and crack, local crushing should also be considered in the simulations of ice-propeller interaction.However, the mesh applied in the FEM simulation is always too coarse to present the feature of local crushing, i.e., such microscopic failure is difficult to be explicitly modeled in the simulation.Hilding et al. [27] proposed a numerical approach, which took local crushing as a product of multiple micro-crack branching in the contact zone.An elastoplastic linear softening constitutive law, taken as a homogenization approach, was applied and have been proved to be effective by Wang et al. [25].This constitutive law is applied to the ice bulk elements to model local crushing in the simulations, and its hardening curve is shown in Fig. 3. σY is compressive strength which is the initial point of crushing.After experiencing a linear softening phase, crushed strain εc is obtained and henceforward the ice material is seen as totally crushed.When the strain comes to failure strain εf, the ice bulk element is completely damaged and deleted.

Numerical setup for propeller ice milling
Numerical model was established according to ice milling model tests performed by Wang et al. [11] in the ice basin of Institute for Ocean Technology, National Research Council of Canada (IOT).The configuration of the model tests is shown in Fig. 4. Podded propulsor was fixed to a carriage and towed forward at a constant speed to mill intact ice block.The propeller was designed with reference to Stone Marine Meridian series [28].The blades of the propeller were thicken for the safety of shipping in ice-covered waters.The  There are four blades in total.The diameter of propeller is 300mm.The mean pitch/diameter ratio is 0.76 and the expanded blade area ratio is 0.669.The model of propeller after meshing is shown in Fig. 5.The propeller is taken as rigid body so that deformation of propeller under ice action is neglected.The material properties of the propeller are given in Table 1.The model of ice block is built up as shown in Fig. 6.The size is 300mm in length, 70mm in width and height, which was determined by a prior verification study to eliminate the influence of boundary effect and meantime guarantee the calculation efficiency.The ice model was made of hexahedral ice bulk elements and cohesive elements.Both ice local crushing failure and fracture failure are considered.The elastoplastic constitutive law is used to simulate local crushing failure, while cohesive element is used to simulate fracture failure.The ice-propeller milling model is set up in LS-DYNA for numerical simulation.The propeller rotates at a constant angular speed of 5 rps clockwise and moves along minus X axis at 0.5 m/s.The blades mill ice with back side.The cutting depth is 35 mm, which is the same as model test setup.Moreover, in order to simulate propeller milling with infinitely wide ice block, the boundary of both sides of ice block are fixed in 6 degrees of freedom as shown in Fig. 7.The most commonly used contact algorithm built-in LS-DYNA is penalty function method, which is defined by the keyword of CONTACT.Command of CONTACT-ERODING-NODES-TO-SURFACE is suitable for modeling the contact between two bodies with considerably different stiffness, which is used to detect ice-propeller contact and calculate the interaction force.The friction coefficient between ice and propeller is set to 0.2.The ice-ice contact is detected by using another command named CONTACT-ERODING-SINGLE-SURFACE, which is suitable for modeling the contact between two bodies with similar stiffness.The corresponding friction coefficient is set to 0.1.The explicit time As the size of ice decreases, the ice loads tend to increase.However, this tendency becomes weak and the ice loads tend to be convergent when the size decreases to a certain value.It is clear from Fig. 8 that the simulated ice loads using S2 mesh are close to those using the finest mesh, S1 mesh.Table 3 shows the average value and the standard deviation of calculated ice loads.The relative errors compared to the corresponding ice loads using S1 mesh are also presented in brackets.It should be noted that the relative errors for S2 mesh are less than 5%.Therefore, balancing the computation efficiency and accuracy, S2 mesh is used in the following simulation cases.In the model tests carried out by Wang et al. [11], a load cell was installed on one of the four blades, which is also named as key blade.The resultant ice load and moment acting on the key blade from the simulations are used for comparison with model test results.The time series of simulated and measured ice loads within 1 second are presented for comparison in Fig. 9.It can be seen that there are five typical peaks for both resultant forces and moments during the time interval.The peaks and trends of ice load signals from simulation agree well with model test results.Table 4 and Table 5 show the peak values and the moment when they occur for both simulated and measured ice loads.The relative errors between simulation and model test results are also included in the bracket.It is found that the peak values and time from simulation coincide well with model test results except the second peak.The relative errors between averaged ice forces and moments are as low as 0.4% and 9.1% respectively, which demonstrates that the simulated results are valid and effective.The trends of the simulated and measured ice loads for one blade cutting process differ from each other, which reflects the inherent random nature of ice-structure interaction process.It should be noted that the second peak of measured moment is obviously smaller than the other peaks.The reason may be that the material property of ice block changes during milling process.Moreover, during the time when the key blade is not interacting with the ice, it is exposed to hydrodynamic load at low level.Since the effect of seawater is not considered in the numerical model, the resulting load is taken as zero in the simulation.In order to give a deeper comparison of the simulated and measured ice loads for one blade cutting process, Fig. 10 shows the simulated and measured ice loads between 0.48s~0.60sas an example.The fluctuated behaviour of simulated ice loads is caused by the local crushing and fracture of ice material, which reflects well the actual ice failure mode during ice-propeller interaction.The curves of measured loads are much more smooth than the simulated ones.This indicates that the measured curves are the filtered results, however, their sampling and filtered criterion are unknown.In order to facilitate the comparison, the simulated results are filtered by a low frequency filter and the filtered results are also shown in Fig. 10.It can be seen that the shape and trend of main peak in the smoothed load curves obtained from simulation approach those measured in the experiment.The low peaks arisen in the experimental loads are caused by the hydrodynamic loads, which are not considered in the simulations.Thus, it is reasonable that the smoothed simulation load levels are slightly lower than the measured ones due to neglect of hydrodynamic loads in the simulations.The stress nephograms at different time instants are given in Fig. 11, from which the physical process of propeller milling with ice can be observed clearly.At the initial stage of the interaction, the blade starts to intrude ice and crushes against the ice block to cut the ice.As the intrusion process continues, the edge of the blade keeps crushing with ice block and paves a way to pass through.Simultaneously, the contact area between ice block and blade surface increases as the propeller disc advances into ice block.The blade surface may collide with ice block and ice pieces chopped from milling process to induce additional contact force as well.The deformation of ice fragments during the ice crushing process can be clearly observed in Fig. 11.As the edge of blade leaves from ice block, a piece of ice is cut off from the original ice block.It is clear to see from Fig. 11(d) that the cutting edge of ice block after milling is irregular which shows that the fracture and cracking happen during milling process.The reason is that the relative speed between ice block and blade is high at the contact point when the screw rotates fast.This leads to a high strain rate that may exceed the limit (10 -3 ) of ductile-brittle transition area according to Carney et al. [30].Thus, the ice block behaves brittle failure mode.The ice pieces chopped from ice block would continue to collide with rotating blade and get shattered into smaller pieces.In addition, it can be seen from Fig. 11 Numerical simulation of ice milling loads on propeller Yang Wang, Hao Yu, Haihua Zhang blade with cohesive element method 120 that the highest effective stress appears at the contact area between the blade and ice block.
With the continuous intrusion of blade into the ice block, more medium effective stress zones discretely distribute in the whole area of the ice block.This result is reasonable and reflects the propagation of stress between ice elements within the ice block.Fig. 12 shows several time series of ice forces and moments in 6 Dofs during time interval of 0.15s ~ 0.2s.It can be seen clearly that the contact force has a trend of increasing first and decreasing after a peak value.This is mainly dependent on contact area between ice block and blade.The contact force augments to maximum when the contact area reaches the largest.Moreover, increasing rate of ice force is higher than decreasing rate clearly.This may be affected by shape of the blade.On the other side, ice pieces cut off may still interact with propeller until they leave away from the propeller disk completely.The resulting impact force, as additional action, prevents the total ice load from dropping quickly.The longitudinal force FX and moment MY are relatively larger than the other loads, which may threaten the structural safety of propeller.Therefore, these two loads are selected as the objects for comparison in the following study to analyze the influence of different key factors on the simulated loads.
0.15 0.16 0.17 0.18 0.19 Fig. 12.Time series of simulated ice force and moment in 6 DoFs exposed to propeller

Parametric study
Based on the numerical model, the effects of rotation speed, advance velocity and cutting depth on ice load in longitudinal and transverse directions are analyzed.The parameter setups of the simulation are kept the same as those used in Section 3 if not stated else.The calculated results from all cases are compared with the basic case.The average value and standard deviation of ice load components FX and MY are of interest herein.The load components from other cases are shown in percentage compared to the results of the basic case whose components are taken as 100%.

Rotation speed
In this section, five different rotation speeds are used for investigating how much the ice load components would be affected.They are 3, 4, 5, 6 and 7rps, where 5rps was applied to the basic case.The average values and standard deviations of ice load components FX and MY under different rotation speeds are given in Fig. 13.It is found that both load components FX and MY show similar trend as rotation speed increases.The average and standard deviation descend successively as the rotation speed increase.Especially for the average, it decreases in Feng Wang, Zao-Jian Zou, Li Zhou Numerical simulation of ice loads on propeller Yang Wang, Hao Yu, Haihua Zhang blade with cohesive element method 122 a linear way approximately.Fig. 14 presents the time series of FX and MY with rotation speeds of 3, 5 and 7 rps within the same range of phase angle from simulations.It is seen from Fig. 14 that the effect of rotation speed is minor at the load increasing stage, but major in the load dropping stage.The more slowly the propeller rotates, the larger the peak value of load time series tends to become.This is mainly attributed to that it takes longer for full contact of icepropeller under the condition of slow propeller rotation.Moreover, the possibility of chopped ice pieces colliding with blades increases, which may result in slow declining of ice loads when the blade leaves off ice block.Therefore, to speed up propeller rotation to some extent could be helpful to break ice block efficiently and also reduce ice loads to improve the safety of propeller operation in ice-covered waters.Five different advance velocities are used for investigating how much the ice load components would be influenced.The screw is advancing straightly to the ice block at 0.2, 0.35, 0.5, 0.65 and 0.8 m/s respectively, where 0.5 m/s is used in the basic case.Fig. 15 shows the curves of average values and standard deviations of ice load components FX and MY as function of advance velocity.It is clear that both load components FX and MY increase significantly as the advance velocity increases, especially for average ice load components.
Fig. 16 presents time series of ice load FX and MY at advance velocities of 0.2, 0.5 and 0.8 m/s.We can see that the simulated ice loads are affected by relative velocity between ice and propeller strongly.During the ice milling process, the ice loads increase significantly with the increasing velocity.The peak of ice loads at the highest velocity 0.8 m/s is almost twice the peak loads from the basic case with 0.5 m/s.This may pose a great challenge to the safe use of propeller in ice.Therefore, it is suggested that the advance velocity should be well controlled to avoid damage of propeller.Five cutting depths are used to study the influence of cutting depth on ice loads.The cutting depths of 21, 24.5, 28, 31.5 and 35 mm are set up in the simulations, where 28 mm is used in the basic case.Fig. 17 shows variations of average values and standard deviations of ice loads as a function of cutting depth.It is clear that both load components FX and MY increase significantly as the cutting depth increases, especially for average ice load components.
Fig. 18 presents time series of ice load FX and MY at cutting depths of 21, 28 and 35 mm.It is seen that as the velocity increases, the ice load level increases at both descending and ascending stages.This trend is caused by the long interaction process and large contact area for ice and blade with big cutting depth.Moreover, the edge of blade enters ice block ahead of expected time when cutting depth become large.Simultaneously, it is found that ice loads decrease to zero at the same time for different cases.The reason is that at the exiting stage the edge of the blade leaves away from the ice block before the surface.From Fig. 11 (d), it can observed that the contact force is small and can be neglected though the surface still interacts with chopped ice pieces.Thus, the time when the load turns to zero depends on when the edge of blade exits from ice block.Fig. 19 shows the stress nephograms with different cutting depths at the same moment before the blade exits from ice block.From the figure, it can be found that the exiting time of the blade edge for different cutting depths is almost the same and thus the loads drop to zero simultaneously.This conclusion coincides with that obtained by Wang et al. [15].

Conclusion
In this paper, we develop a numerical method to simulate continuous ice-propeller milling process based on an elastoplastic softening constitutive law with cohesive element method, which is implemented in commercial finite element software LS-DYNA.The present numerical method models ice cracking and crushing behaviors happened in the ice-propeller interaction, and the numerical results are validated by comparing with the available experimental results.Then, the effects of rotation speed of screw, advance velocity and cutting depth on dominant ice load components are analyzed and discussed.The main conclusions can be drawn as follows: (1) There exists a large fluctuation of simulated ice loads due to local failures of ice and small pieces of ice generated.The simulated time series of ice loads coincide well with model test results with respect to the peak values and the time when the peak values occur, indicating that the linear softening constitutive law used in cohesive element method is capable of modeling propeller-ice milling process.(2) Ice loads are affected by setup of ice mesh size.The loads will increase with decreasing ice mesh size.The ice loads tend to be steady when the size increases to certain level.(3) According to the sensitivity study, the ice loads are affected by rotation speed of screw, advance velocity and cutting depth to varying degrees in terms of average and standard deviation values.The dominant ice loads decrease as the rotation speed increases and increase as the advance velocity or cutting depth increases.The ice loads are influenced by the advance velocity most compared to the other two factors.It is shown that the present numerical method has a good prospect in modeling icepropeller interactions.However, it should be noted that all conclusions are drawn based on the present simulation.Besides, the sensitivity of ice properties to the ice loads is not investigated herein and should be studied in the future.

Fig. 2 .
Fig. 2. Relationship between cohesive force and crack separation (Left: model I; Right: model II)

Fig. 3 .
Fig. 3. Elastoplastic constitutive law[27] Numerical simulation of ice milling loads on propeller Feng Wang, Zao-Jian Zou, Li Zhou blade with cohesive element method Yang Wang, Hao Yu, Haihua Zhang 113 EG/AD/S ice was used in the experiments.The EG/AD/S ice is generated by the cooling of a diluted aqueous solution of ethylene glycol (EG), aliphatic detergent (AD), and sugar (S) under approximately -20°C, which can provide similar material properties with the first-year columnar sea ice.

Fig. 6 .
Fig. 6.The model of ice block after meshing

Fig. 9 .
Fig. 9. Simulated and measured ice loads on the key blade

Fig. 10 .
Fig. 10.Simulated and measured ice loads for one blade cutting process

Fig. 16 .
Fig. 16.Time series of simulated ice load FX and MY with different advance velocities

17 .Fig. 18 .Fig. 19 .
Fig. 18.Time series of simulated ice load FX and MY with different cutting depths

Table 1
Material properties of propeller model Feng Wang, Zao-Jian Zou, Li Zhou Numerical simulation of ice milling loads on propeller Yang Wang, Hao Yu, Haihua Zhang blade with cohesive element method 114 Table 2 presents the properties of ice bulk elements and cohesive elements applied in the simulations.The compressive strength and shear strength (S0) are chosen based on the measured values in the model test provided by Wang et al. [11].The other parameters are determined with reference to the work by Timco et al. [29] on the engineering properties of first-year sea ice.

Table 3
Comparison of ice loads for cases with different ice mesh sizes (b) Resultant moment Fig.8.Calculated time series of ice loads with different ice mesh sizes

Table 4
Comparison of peak values for ice force

Table 5
Comparison of peak values for ice moment Statistics of ice load FX and MY based on simulations with different advance velocities