Macro-phenomenological high-strain-rate elastic fracture model for ice-impact simulations

Jing FAN , Qingho YUAN , Fulei JING , Gungchen BAI , Xiuli SHEN ,*

a School of Energy and Power Engineering, Beihang University, Beijing 100191, China

b Beijing Key Laboratory of Aero-Engine Structure and Strength, Beijing 100191, China

c Collaborative Innovation Center of Advanced Aero-Engine, Beijing 100191, China

d Aero-Engine Academy of China, Aero Engine (Group) Corporation of China, Beijing 101304, China

KEYWORDS Crack propagation;Elastic fracture;High strain rate;Ice impact;Impact force;Material model

Abstract The ice impact can cause a severe damage to an aircraft’s exposed structure,thus,requiring its prevention.The numerical simulation represents an effective method to overcome this challenge.The establishment of the ice material model is critical.However, ice is not a common structural material and exhibits an extremely complex material behavior.The material models of ice reported so far are not able to accurately simulate the ice behavior at high strain rates.This study proposes a novel high-precision macro-phenomenological elastic fracture model based on the brittle behavior of ice at high strain rates.The developed model has been compared with five reported models by using the smoothed particle hydrodynamics method so as to simulate the ice-impact process with respect to the impact speeds and ice shapes.The important metrics and phenomena (impact force history,deformation and fragmentation of the ice projectile and deflection of the target)were compared with the experimental data reported in the literature.The findings obtained from the developed model are observed to be most consistent with the experimental data,which demonstrates that the model represents the basic physics and phenomena governing the ice impact at high strain rates.The developed model includes a relatively fewer number of material parameters.Further,the used parameters have a clear physical meaning and can be directly obtained through experiments.Moreover,no adjustment of any material parameter is needed,and the consumption duration is also acceptable.These advantages indicate that the developed model is suitable for simulating the iceimpact process and can be applied for the anti-ice impact design in aviation.

1.Introduction

The ice impact and resultant dynamic loading have always been the prominent challenges in aviation, which can cause a significant plastic deformation of the external surfaces and exposed components of aircrafts.1–3The ice ingestion in the engines can lead to the power loss.4,5Even a brief impact during takeoff, flight and landing can pose a serious threat to the flight safety.In recent years, the threat has become more serious with the application of the composite materials in the exposed structures, as the damage in such materials is barely visible and is concealed below the surface.6Therefore, such structures should be designed with sufficient tolerance to the impact of ice.Thus, in-depth studies are needed to explore the influence of ice (mass, impact angle and velocity) on the structural damage.In order to simulate the ice-impact process,the material model of ice is the most commonly employed approach.The ice impacts in aviation represent the highstrain-rate impact issues.3,6,7However, due to the complex material behavior exhibited by ice, the appropriate material models to accurately simulate the dynamic behavior of ice during the high-strain-rate impacts are still lacking.As a result,there is a certain gap between the calculated results (e.g., the impact maximum load and structural damage) and actual findings.

Ice possesses at least 13 crystalline phases and 2 amorphous states.8Of particular interest in aviation is the Ih polycrystalline ice(hexagonal crystal structure).It is the most common ice form on Earth and forms on cooling the liquid water below 0°C at ambient pressure.3,7,9The other crystal forms of ice are thermodynamically stable only at high pressures3(thus, the term‘‘ice”in this study refers to Ih ice).The mechanical behavior of ice has been widely studied for glacier research,10,11arctic ship transportation12,13and offshore facilities.14It has been reported that the mechanical properties of ice exhibit a strong strain-rate dependency.Ice undergoes a ductile to brittle transition as the strain rate increases from a low value to 10-3s-1.9,15,16The mechanical properties of ice are also affected by the temperature,microstructure,and salinity.16–18However, a majority of the studies have been conducted at low strain rates (10-7–10-1s-1),9,16,19–21such as creep and quasi-static experiments, which are not suitable for exploring the ice-impact in aviation.

Only a small number of literature studies have focused on the mechanical behavior of ice at high strain rates.In 1960s,the British Royal Aircraft Establishment (RAE) conducted a series of hail hitting aluminum plate tests to study the structural damage caused by hail.22However,the mechanical properties of ice at high strain rates were not studied.Several studies have reported the mechanical properties of ice at high strain rates in recent years.Jones,23Combescure et al.24studied the compression properties of ice by using a high-speed universal testing machine up to 50 s-1.Fasanella et al.25conducted the dynamic crush testing of ice by employing a bungee-assisted drop tower to obtain the stress–strain behavior of ice at 100–300 s-1.Dutta et al.,26,27, Shazly et al.28and Kim and Keune18obtained the compressive strength of ice by using a split Hopkinson bar at 1–100 s-1, 90–882 s-1,400–2600 s-1, respectively.It was reported3,9,15,16,29that the response of ice at high strain rates appeared to be independent of the microstructure, and ice exhibited a brittle behavior attributed to the tensile and compressive properties.The tensile strength of ice is noted to be independent of the strain rate,whereas the compressive strength is generally enhanced with the strain rate.Meanwhile, ice demonstrates the fluid properties on enhancing the strain rate.

In summary,ice possesses the characteristics of strain-ratesensitivity, ductile-to-brittle transition, solid-to-fluid transition, distinct tension and compression behaviors and correspondingly separate failure modes.As can be observed, the material behavior of ice is extremely complex and exhibits a variation in the properties as a function of the strain-rate stages.Such a behavior is difficult to be accurately characterized by an ideal model.In earlier studies,the efforts were made to establish a phenomenological elastic–plastic model to simulate the ice-impact behavior in aviation.Kim and Kedward1first established a simple J2 elastic–plastic model of ice based on the plastic strain and tensile pressure failure criteria.Anghileri et al.30proposed a J2 elastic–plastic hydrodynamic model of ice with the tensile pressure failure.However, these models did not consider the strain-rate-sensitivity, along with the tension and compression behaviors.In actual applications, a few key parameters of these models are needed to be adjusted as a function of the impact speeds to match the corresponding experimental data.This significantly limits their ability to predict the actual impact behavior at the unknown impact speeds.Tippmann et al.6improved the Kim’s model and proposed a rate-sensitive perfectly elastic–plastic model with the tensile pressure failure.However,it did not consider the different tension and compression behaviors.The calculated impact force was noted to be too high at certain velocities, and the simulation of the large deformation, crushing and spreading of the broken debris was not ideal.Carney et al.7proposed a ratedependent elastic–plastic model with the tensile and compressive pressure failures.Sain and Narasimhan31further added a damage evolution law to the Carney’s model.Based on the Drucker-Prager plasticity criteria,Pernas-Sanchez et al.3developed a rate-sensitivity elastic–plastic model to allow for the different tension and compression behaviors of ice.Even though these models consider relatively comprehensive characteristics of ice, however, the calculation results at various impact speeds of interest in aviation are noted to be inferior to the results obtained from the previous simple models.

In case the speed of the ice projectile is above 25 m/s (the terminal velocity of the hail ice with a diameter of 25.4 mm when its weight is balanced by the air resistance), the corresponding strain rate exceeds 10 s-1,18which is much higher than the transition strain rate of 10-3s-1, and the relative impact velocity of an aircraft during flight is significantly higher than this value.Ice mainly exhibits the brittle behavior during hitting an aircraft.Therefore,it is believed that the elastic fracture model can effectively simulate the ice-impact challenge in aviation,and the plastic stage of ice can be ignored.At present, a few research studies have proposed the elastic fracture models, however, further improvements are needed to be implemented.Chuzel32and Ortiz et al.33proposed a ratedependent elastic damage model based on the Maras concrete failure model to simulate the ice impact in aviation.It included a damage evolution model,an equivalent strain to simulate the material failure and the delayed damage effect to avoid the challenges of the infinite strain rate and zero fracture energy caused by the mesh size.However, the model only considers the tensile failure, and the complex damage evolution model and introduction of the delayed effect significantly increases the complexity of the model and introduces many parameters whose values are difficult to be determined,thus,requiring the manual adjusting based on experience.Wang et al.34proposed an elastic-brittle model considering the energy release rate failure to simulate the fragmentation of the ice cover caused by the underwater explosion.In this process, the speed of the shock wave reaches 1500 m/s, which corresponds to a high strain rate.However,the model did not consider the strain rate effect, which required a manual parameter adjustment.

Owing to these reasons, this study aims to establish a novel high-precision macro-phenomenological elastic fracture model of ice to simulate the ice-impact process at high strain rates without adjusting any parameters, for aeronautical engineering applications.The model comprehensively considers the characteristics of ice at high strain rates, e.g., brittle fracture, strong strain-rate-sensitivity, solid-to-fluid transition,tension and compression behaviors, corresponding separate failure modes, etc.Further, the parameters used in the model have a clear physical meaning and can be obtained directly through experiments without any adjustment.Moreover, as the plasticity is not considered, the model possesses a small number of parameters.On the premise of the simplicity of the model, the accuracy of the numerical simulation can also be guaranteed.These advantages make it convenient for use in the aeronautical applications.The new model has also been compared with four typical elastic–plastic models as well as an elastic model.The obtained results have also been verified by the corresponding experimental data from the literature so as to verify its performance and wide applicability.In addition, the current ice models lack the horizontal quantitative comparison under unified conditions.Thus, the comparison can effectively explore the advantages and disadvantages of the model itself,thus,helping to propose and establish a superior material model for ice.

The layout of this study is as follows.In Section 2, the elastic–plastic models as well as the elastic model are briefly introduced.Afterwards, the novel elastic fracture model has been proposed in detail.In Section 3, the new model has been compared with the elastic–plastic models as well as the elastic model, respectively.The results are verified by the corresponding experimental data obtained from the literature.In Section 4, the key challenges and new discoveries when establishing an ice material model at high strain rates have been discussed.In Section 5, the conclusions are summarized.

2.Ice material model

The material model generally includes a constitutive relation and a failure criterion.The constitutive relation essentially reflects the relation between the stress and strain.For the linear elastic materials, the stress and strain obey the generalized Hooke’s law.For the elastic–plastic materials, the plastic stress–strain relations are needed to be considered, such as the perfectly plastic model, linear plastic strain hardening model and power law plastic strain hardening model.Currently, two types of constitutive relations, linear elastic model and elastic–plastic model,are used to simulate the deformation of ice at high strain rates.

The failure criterion describes the failure behavior of a material.As is well known,during an impact,a total stress tensor σ can be partitioned into a deviatoric stress component and a hydrostatic pressure component.

where sijis the deviatoric stress component, and pδijis the hydrostatic pressure component.The subscripts i,j and k indicate the coordinate directions.It should be noted that p is negative during tension, while it is positive during compression.Currently, three types of failure criteria, including the failure pressure, failure strain and failure energy release rate, are employed to simulate the failure of ice at high strain rates.Based on the different constitutive relations and failure criteria, a number of ice material models have been proposed.

In the earlier studies, the efforts were made to establish an elastic–plastic fracture model to characterize the properties of ice.However, the material properties of ice vary as a function of the strain rate stages.Fig.1 presents the transition of ice from the ductile fracture to the brittle fracture as the strain rate increases, with 10-3s-1observed as the critical strain rate at which the transition occurs.During ice impact in aviation,the strain rate generated by the collision is much higher than 10-3s-1,and the ice mainly exhibits the brittle fracture.Therefore, the elastic fracture model can effectively reflect the dynamic behavior of the ice hitting aircrafts.In view of the fact that the elastic fracture model needs further improvement,this study proposes a novel elastic fracture model.In order to quantitatively compare the difference between the elastic fracture and elastic–plastic fracture models in the aviation field as well as to identify the advantages of the developed model compared to the existing elastic fracture model, five typical literature reported material models (including four elastic–plastic models and an elastic fracture model) have been compared with respect to their advantages and limitations.The typical material models of ice have been briefly summarized first, followed by the description of the new material model.

2.1.The elastic–plastic failure models without considering the rate dependency

Kim and Kedward1first proposed a simple J2 elastic–plastic failure model to simulate the hailstone impact in aviation.The linear plastic strain strengthening and fluid effect of ice after failure were considered in the model, and the failure criteria included the plastic failure strain and tensile failure pressure.A linear Equation Of State (EOS) was used to calculate the ice pressure.Once one of the two criteria is achieved, the unit is considered to have failed.Meanwhile, the shear stress components are relaxed to zero, and the material can sustain compressive stresses only, thus, simulating the fluid properties of ice.The corresponding material parameters are shown in Table 11.

The J2 elastic–plastic hydrodynamic model was first employed by Anghileri et al.30to simulate the ice impact.It differs from the Kim’s model as it uses a nonlinear polynomial EOS and a tensile pressure failure criterion.The values of the material parameters are shown in Table 2.30

Table 1 Parameters of J2 elastic–plastic failure model.1

Table 2 Parameters of J2 elastic–plastic hydrodynamic model.30

The limitations of the two models are:

(1) The yield stress is not a function of the strain rate.

(2) The plastic hardening modulus and failure metrics of the material are tuned arbitrarilyto match the ballistic testdata.

(3) The tension and compression behaviors are symmetrical.

(4) The failure strain in the Kim’s model under tension and compression is the same and is not a function of the strain rate.

(5) The compressive failure is not considered in the Anghileri’s model.

2.2.The rate sensitive perfectly plastic model with failure

Tippmann et al.6improved the Kim’s model by taking the strain rate effects into account and proposed a strain ratesensitive perfectly plastic model.Based on the experimental statistics, the relation between the strain rate and compressive strength was fitted6.Like Anghileri’s model, this model only employs a tensile pressure failure criterion.The main material parameters are listed in Table 3.6The limitations of this model are:

Table 3 Parameters of rate sensitive perfectly plastic model.6

1.The tension and compression behaviors are symmetrical.

2.The compressive failure is not considered.

3.The value of the tensile failure pressure is chosen based on the parametric study with a certain degree of uncertainty.

2.3.The rate sensitive elastic–plastic model with separate failure modes

As ice exhibits different behaviors under tensile and compressive loads, Carney et al.7proposed a model by simultaneously considering the strain rates, tensile and compressive constitutive relations and failure criteria.The failure criterion was based on the failure pressure.The pressure was evaluated by employing a tabulated EOS with compaction.The values of the material parameters are shown in Table 4.7The relation between the strain rate and compressive strength as well as the EOS parameter values can be found in Ref.7.

The limitation of this model is that the values of the initial flow stress and failure pressure in the model are chosen based on the parametric analysis.

Table 4 Parameters of the rate sensitive elastic–plastic model for ice with separate failure modes.7

2.4.The rate sensitive elastic damage model

Different from the elastic–plastic model defined earlier,Chuzel et al.32proposed a rate sensitive elastic model with an equivalent strain failure based on the brittle behavior of ice at high strain rates.The model is based on the Maras concrete failure model, which has the following hypotheses:

? The elastic behavior is coupled with the isotropic damage;

? The damage is only caused by the tensile strength;

? The damage in tension and compression evolves differently.

The damage evolution in traction and compression can be calculated as:

The damage evolution is described by the equivalent strainε～, initial strain εD0and four parameters At, Bt, Acand Bc.The calculated traction and compression damage is updated using the delay effect in order to reduce the damage rate and mesh dependence.

where τ is the characteristic time; lcaris the characteristic length, which corresponds to the minimum possible size of a completely damaged zone; and c is the sound speed in the ice material,Dncis the uncorrected damage with delay effect being taken into account.Table 532presents the ice model parameter values for lcar= 16 mm.

The strain-rate effect as well as the tension and compression behaviors are considered simultaneously in this model.However, the model does not consider the compressive failure.The damage evolution is too complex,and the delayed damage effect further enhances the model complexity.Thus,it includes many parameters (At, Bt, Ac, Bc, lcar) whose values have a significant influence on the simulation results and are difficult to be determined, thus, requiring manual adjustments to fit the experimental data.

Table 5 Parameters of rate sensitive elastic damage model for lcar = 16 mm.32

2.5.The new rate sensitive elastic fracture model

This study aims to establish a new elastic fracture model based on the brittle behavior of ice at high strain rates.The reason is that there is almost no plastic stage (see Fig.1) in the stress–strain curve of the brittle materials, which can be considered as the fracture occurring in the elastic range.Therefore, this model considers that ice will fail in the elastic stage in case it reaches its tensile or compressive strength.The tensile strength does not depend on the strain rate, which is a constant.The specific value is obtained from the experimental data in Ref.9.The compressive strength depends on the strain rate.Therefore, the strain-rate effect on the compressive strength is needed to be considered.The literature studies3,6–7,19have pointed out that the compressive strength and strain rate of ice present a power exponential relation.Thus, the linear-log curves can fit the compressive strength vs.strain rate data.Based on the experimental data reported by Jones23, Kim18and Shazly et al.,28the relationship between the compressive strength and strain rate of ice in the high strain rate stage(>10-1s-1) has been fitted.It is worth noting that the Kim’s model18represents the highest ice strain rate test data in the existing literature.Fig.2 and Table 6 show the relation between the compressive strength amplification factor and strain rate.In summary,the constitutive relations of the model are,

where σ is the stress tensor;ε is the strain tensor;C is the elastic stiffness tensor;the superscripts c and t represent compression and tension,respectively;σcy0is the initial compressive strength of ice; σcyis the dynamic compressive strength; ˙ε is the strain rate; and A is the compressive strength ratio.

Fig.2 Fitting of compressive strength versus strain rate data.

Table 6 Strain rate dependence of compressive strength.

The failure criterion in the developed model is based on the hydrostatic pressure evaluated by using an EOS, which describes the relation between pressure, material density and temperature.This study adopts a polynomial EOS owing to its simple form.The linear EOS generally leads to the unrealistic pressure oscillations in the damaged ice.7In order to avoid this, a high-order polynomial EOS has been employed herein:

where ρ is the current density(the ratio of mass to current volume);ρ0is the reference density at the non-deformed state;E is the internal energy;and C0–C6are the polynomial coefficients.Given that the duration of an impact experiment is in milliseconds,the ice temperature is assumed to be constant during the simulation.Therefore,the temperature influence is not considered in this study,and C4–C6are set to zero.The values of the other parameters are listed in Table 7.35

In case the hydrostatic pressure exceeds the tensile or compressive strength, the material is considered to be invalid.As the tensile and compression behaviors of ice are different, this study defines the failure criteria in tension and compression separately.Ice exhibits the fluid properties on increasing the strain rate, and the failure of ice at high strain rates is mainly manifested as the brittle fracture, which still has a certain degree of residual compressive strength after failure.

In this model, in case the compressive hydrostatic pressure exceeds the compressive strength,the deviatoric stress tensor is set to zero.However, the compressive hydrostatic pressure is not limited to the compressive strength, thus, allowing to maintain the load on the impact target.Thus, it can simulate the solid-to-fluid transition as well as the ability of ice, which can still sustain a certain degree of compression after a fracture.As the compressive strength depends on the strain rate,the compressive failure pressure should be also a function of the strain rate.

Table 7 Values of coefficients in Eq.(11).35

In case the tensile hydrostatic pressure exceeds the tensile strength, the deviatoric stress tensor and tensile hydrostatic pressure are set to zero.As the tensile strength of ice is not sensitive to the strain rate, the tensile failure pressure is noted to be a constant.In summary, the failure criterion in this model can be written as Eqs.(12) and (13).

Finally, the parameters used in the developed elastic fracture model are shown in Table 8.It is worth mentioning that the data in Table 8 can be replaced with respect to the ice sample data at a different temperature,salinity or impurity.As the anti-ice impact safety design needs to consider the most dangerous state,the values here represent the pure water ice without any impurities and bubbles.Generally, the formation temperature of hail in nature is in the range -20 °C to-10 °C.A majority of the mechanical performance analysis of ice under high strain rates is carried out in this temperature range.The basic material data of ice in Table 8 is the published experimental data.

Table 9 presents a summary of the characteristics and shortcomings of the various material models.As compared with the previous models, the developed model presents the following notable features:

(1).Based on the brittle fracture characteristics of ice at high strain rates,the absence of plasticity leads to a relatively few number of parameters.

(2).All parameters have a clear physic meaning and can be measured experimentally.

Table 8 Parameters of new elastic fracture model at -20 °C to -10 °C.

Table 9 Summary of characteristics and shortcomings of the material models.

(3).It simultaneously considers the strain-rate-sensitivity,solid-to-fluid transition, tension and compression characteristics and corresponding separate failure modes.

(4).The failure criterion is based on the hydrostatic pressure,which follows the failure mode of the brittle materials.

(5).The strength and failure hydrostatic pressure under compression are functions of the strain rate.

(6).The strain-rate-related constitutive equation, dynamic failure criterion and failure pressure equal to the failure strength do not require the model to adjust any parameter for simulating the ice projectiles of different sizes and speeds.

3.Numerical simulations and comparison with the experiments

Large deformation and fragmentation will occur during ice impact,and the Finite Element Method (FEM) and Arbitrary Lagrangian Eulerian (ALE) method suffer from severe limitations in dealing with the high-speed impacts, such as the mesh distortion resulting in the remeshing and remapping of the state variables in case the domain shape deformations are large.29Meanwhile, crack propagation needs to be simulated during the simulation of fragmentation.However, FEM and ALE depend on the grid.Thus,these are not suitable for dealing with the discontinuity problem inconsistent with the original grid line.The Smoothed Particle Hydrodynamics (SPH)method was originally used to overcome the astrophysical and cosmological challenges,which was subsequently extended to the Computational Fluid Dynamics (CFD) and continuum mechanics.The SPH method avoids the dependence on grid,where trial functions are constructed in the neighborhood of the discrete points,therefore,no mesh generation and remeshing issues are observed.The basic idea is to approximate a function u(x) on domain Ω by a convolution as.where uh(x) is the approximation of u(x), and w(x-x-) is a compactly-supported function, usually called the kernel function.Three commonly used kernel functions are the cubic spline, Gaussian spline and quartic spline.Further, x-is the spatial coordinate of each point in the neighborhood of the point x.

The SPH method does not need a mesh,and it is essentially a Lagrangian method.Therefore, SPH can accurately capture the material interface without the issue of grid distortion,so as to solve the difficulties of crack growth and fracture simulation of materials.As a mature meshless method, the SPH method has been widely applied to the large deformation and failure caused by impacts, and it has been proven to simulate the ice impacts well.29,30In this study, the realization of the new material model and simulations based on the SPH method has been achieved by using the commercial software LSDYNA.To validate the proposed material model, the experimental data of ice impact under high strain rates was needed.The comparison with the verification tests with respect to the previously proposed elastic–plastic models and elastic models has been subsequently made in this study.

3.1.Comparison with elastic–plastic failure model

A majority of the classic elastic–plastic models have been verified based on the test device shown in Fig.3.During the test,the ice projectiles are launched against a force measurement‘‘rigid”bar using a nitrogen gas gun.1,6The purpose of such experiments is to obtain the impact force history and fracture images of ice.Further, the target has almost no deformation during the test.The bubble-free freshwater hail ice is generated using a cavity mold forming two hemispherical halves.Tippmann et al.6conducted a series of impact tests (19.3–193.7 m/s)using the polycrystalline hail ice at an average temperature of-17.5°C based on this device.For this reason,the test data reported by Tippmann et al.6at the typical speeds(37.3 m/s, 81.2 m/s, 189.2 m/s) was selected as the verification data for analyzing the real speed range of the hail ice impacting an aircraft.The four elastic–plastic models described in Section 2 as well as the material model proposed in this study have been used to simulate the impact processes.It is worth mentioning that the simulation results of Tippmann et al.were directly obtained from the literature study,6whereas the other elastic–plastic models were reproduced in the LS-DYNA software.The established geometric discrete model is shown in Fig.4.The diameter of the spherical ice was 61 mm, whereas the dimensions of the target were 200 mm × 200 mm × 25.4 mm.The projectile was discretized by the SPH particles, and the flat plate was discretized by using the Lagrangian 8-node hexahedral elements with a central encryption.After the sensitivity analysis of the SPH particle number(as mentioned in the next paragraph), the SPH particle number was 268096.The Lagrangian mesh number was 63700.An automatic pointsurface contact algorithm based on the penalty function method was adopted between the ice and plate.The friction between the ice and metal plate was ignored.36The linear artificial viscosity was 1.2,6and the secondary artificial viscosity was 037.The total simulation time was 0.33 ms, and an automatic time step was adopted.The target was made of the 6061-T6 aluminum alloy.An elastic–plastic material model was used to simulate its behavior.The stress–strain curve of the target is shown in Fig.5.38Here, the Poisson’s ratio and density were 0.33 and 2700 kg/m3, respectively.

Fig.3 Schematic of nitrogen gas cannon setup.

Fig.4 Geometric discrete model.

Fig.5 Stress–strain curve of the 6061-T6 aluminum alloy.38

Fig.6 The impact peak force and consumption time of different SPH particle numbers.

Fig.7 Impact force history of different SPH particle numbers of the new model.

The effect of the SPH particle number on the impact force was investigated to eliminate the influence of the geometric discretization on the calculations.Fig.6 presents the peak force and consumption time (8 Inter Core i7-4790 3.6 GHz CPUs)of four models for different SPH particle numbers at an impact velocity of 60 m/s.The consumption time of Kim’s model takes too much time(44–2233 min)and is not shown in the figure.Fig.7 presents the impact force history of the new model for different SPH particle numbers.As can be observed,in case the number of the discrete particles is small, the impact force curve exhibits the unrealistic spikes.These unrealistic spikes tend to disappear and the peak force curves tend to converge gradually as the number of the SPH particles increases.On the other hand,an increase in the number of the SPH particles reduces the computational efficiency.Thus, considering the convergence and computational efficiency, an ice model with 268,096 particles was finally adopted.

Jackson et al.39concluded that the dynamic load generated during the ice impact is a key metric in evaluating the damage to the target.Therefore,the impact force history from the simulation results was compared with the experimental results reported by Tippmann et al.6,40As observed from the impact force histories (Figs.8–10), the findings from the developed model at three velocities are most consistent with the experimental data among the five models.It is specifically manifested in the rising and falling trends, magnitude of the peak force,time at the peak force, degree of attenuation, etc.The impact force is noted to swiftly rise to the maximum(around 0.05 ms),followed by a gradual decay.The relative error in the peak force is noted in the range 6.0%–15.3% (Fig.11), whereas the relative error in the duration for achieving the peak force lies in the range 0–14.8% (Fig.12).The calculated durations are in the range 18–22 min (Table 10), which represent an acceptable level.No material parameters were adjusted in the simulations at three speeds.

Fig.8 Impact force history curves at 37.3 m/s.

Fig.9 Impact force history curves at 81.2 m/s.

In contrast, the impact force histories of the four elastic–plastic models are not ideal.The impact peak force in the case of the Kim’s model is too large with a maximum relative error of 35.6% (Fig.11), whereas the peak force appears too late at relatively low velocities (37.3 m/s, 81.2 m/s) (Fig.8, Fig.9)with a maximum relative error of 196.3% (Fig.12), although it is observed to improve to 189.2 m/s (Fig.10)with a relative error of 36.0% (Fig.12).In addition, a long calculation time(569–1847 min) (Table 10) is needed, which is not acceptable for the engineering applications.This may be due to the reason that the model uses two types of failure criteria, which significantly increases the computational complexity.

Fig.10 Impact force history curves at 189.2 m/s.

Fig.11 Relative error in the peak force for different models.

Fig.12 Relative error in the time to attain the peak force for different models.

Table 10 Time consumed by different models for calculation.

Fig.13 Impact sequence in the elastic fracture model and experiments6 at 61.8 m/s.

The impact peak force in the Anghileri’s model is noted to be closer to the test values at relatively low impact velocities(37.3 m/s, 81.2 m/s) with relative errors of 14.4% and 2.7%,respectively (Fig.11).However, at 189.2 m/s, the peak impact force is observed to be extremely large (70.9% relative error),and the impact force history shows obvious oscillations(Fig.10).In addition, the peak force is noted to appear too early (Figs.8–10).

The impact peak force in the case of the Tippmann’s model is too large (relative errors of 36.4% and 43.0%) at relatively low impact velocities (37.3 m/s, 81.2 m/s), and the calculation results at 189.2 m/s are improved (relative error of 18.0%)(Fig.11).In terms of the time to attain the peak force(Fig.12), the relative error increases gradually (14.5%–48.0%) with the velocity.It should be noted that the results in this case were directly obtained from the Ref.6,40, which exhibited the impact force history of 0.15 ms at 81.2 m/s and 189.2 m/s.At the same time, no information about the calculation duration was revealed in the studies.6,40

The Carney’s model exhibits a lower impact force at low impact velocities (37.3 m/s, 81.2 m/s) with relative errors of 60.6% and 33.7%, respectively.It is observed to be improved further to 189.2 m/s with a relative error of 9.6% (Fig.11).In terms of the time to attain the peak force,it is only close to the experimental data at 81.2 m/s with a relative error of 14.8%,and the relative errors for other velocities are too large(Fig.12).

The reason for the gap in the impact force history for the four published models and experimental data is that the mechanical properties of ice at high strain rates are not fully considered in these models.On the other hand,it may be possible that a few key material parameters need to be adjusted in the models, thus, making them fit only at certain velocities.

The failure criterion essentially describes the fracture effect of ice under ultimate load.The dynamic fragmentation of ice can be used to verify whether the failure criterion meets the physical nature of the material.Fig.13 shows the fragmentation of ice at 61.8 m/s as a function of time, using a highspeed camera.As illustrated from the camera observation in Fig.13, the micro cracks first grew away from the initial impact location.As the ice surface had no external constraints,the micro cracks in the contact area gradually converged along the ice surface and formed main cracks longitudinally oriented in the direction of travel.At the same time, the impact force increased rapidly and reached the maximum at around 50 μs.With the transmission of the stress wave, the main cracks further expanded to the rear end.Finally, the instability cracking and collapse were observed to occur.The fracture effect of ice in this study is demonstrated by the variation of material density.It can be observed from the equation of state of ice that the volume of ice will change dynamically during the collision,which is caused by the extrusion and fracture failure.However,the total mass of hail is unchanged, so the density will change dynamically.There is a neighborhood around each particle of SPH, so the density variation within each particle neighborhood can be accurately captured.Subsequently, the process of crack growth and failure can be accurately demonstrated.As observed from the simulation in Fig.13, at t = 0 μs, the initial value of density was 900 kg/m3.After contacting with the metal plate,the ice was squeezed,thus,enhancing the overall density of ice(the density of ice increased from 900 kg/m3to 937 kg/m3at t = 40 μs).The internal hydrostatic pressure of ice increased at the same time.As the pressure reached the given failure criteria, the cracks were observed to occur.Due to the generation of cracks, the density was reduced, thus,enabling the crack extension.The decrease of local density from initial value 900–0 kg/m3reflected the process of crack initiation to complete fracture.For example, the local density of ice decreased from 900 kg/m3to 728 kg/m3at t = 40 μs,which reflected the process of crack initiation.At t = 320 μs, the local density decreased to 7.2 kg/m3, which was close to 0 kg/m3, i.e., it could be considered as the complete fracture.The simulation results are mostly consistent with the tests, which demonstrates that the elastic fracture model with pressure failure proposed in this study reflects the physical mechanism of ice fracture to a certain extent.The total 760 μs calculation cost 40 min.

3.2.Comparison with elastic damage model

The elastic damage model proposed by Chuzel et al.32was verified by using two experiments in the original study.One was carried out by the Office National d’e′tudes Et de Recherches Ae′rospatiales(ONERA)to demonstrate the damage evolution during impact.The experiment was carried out using the mono-crystalline cylindrical ice (height 30 mm, diameter 25 mm) impacting an aluminum AU4G panel of dimensions 200 mm×80 mm×1.2 mm at 62.5 m/s.The plate was simply supported (180 mm distance between the supports), as shown in Fig.14.Different impact angles were tested in order to demonstrate the failure modes: fragmentation at 0°and cleavage at an impact angle of 20°.In the original study, Chuzel et al.32used the SPH method to simulate the cylindrical ice impact process by using the proposed model, followed by the comparison of the simulation results with the experimental data.In order to illustrate the superiority of the model proposed in the current study, the SPH method was used to simulate the impact process, followed by the comparison of the simulation results with the simulation results reported by Chuzel et al.as well as the experimental results.The simulation results reported by Chuzel et al.and experimental results were obtained from the original literature.32The calculation model was established as per the experiment.The ice cylinder and target plate were discretized by the SPH particles and shell elements, respectively.After the SPH particle convergence analysis,an ice model with 169680 SPH particles was adopted.The shell element number was 16000.The target employed an elastic–plastic material model.The stress–strain curve of the target is shown in Fig.15.41Here,the Poisson’s ratio and density were 0.31 and 2800 kg/m3, respectively.The friction coefficient, contact algorithm, time step, etc.are noted to be the same as the simulation in Section 3.1.The comparison between the developed model, Chuzel’s model and experiments is shown in Fig.16.

Fig.14 Geometric discrete model.

Fig.15 Stress–strain curve of aluminum AU4G.41

It can be observed that under normal impact, the cracks gradually initiate and propagates from the head to the tail of the ice cylinder along the direction of motion.Eventually,due to the geometric limitation of the tail, the cracks can not expand further, and the stress wave is reversed forming the straight cracks perpendicular to the direction of motion around the tail.During the impact, the ice front gradually reveals the fluid diffusion.The calculation results obtained using the developed model are significantly consistent with the experimental phenomena.Under the oblique impact, the crack pattern appears to be different from the normal impact,and the crack extends from the beginning till the end.The crack propagation angle is an acute angle with the impact direction and is also consistent with the experiments.

In contrast, in the findings obtained from the Chuzel’s model, the crack initiation and propagation are not obvious in both normal and oblique impact.During the normal impact,the front of the ice cylinder does not show the fluid diffusion form.The observed phenomenon may result from the equivalent strain failure criterion and tensile failure criterion used in the Chuzel’s model.The separate pressure failure criteria under tension and compression adopted in this study seem to be more in line with the physical nature of the ice fracture.

The performance of the ice model was further explored by another experiment which was a hailstone impacting on an aluminum plate carried out by RAE.The experiment consisted of an impact by the spherical ice sample of 25.40 mm diameter onto a 0.91 mm thick aluminum alloy 2014-T4 panel of 305 mm side-edge(Fig.17)fixed with the blind rivets,presenting a free surface of 200 mm×200 mm.In this study,the hailstone and target were discretized by 268096 SPH particles and 40000 shell elements, respectively.The target used an elastic–plastic material model,and the stress–strain curve of the target is shown in Fig.1838.The Poisson’s ratio and density were 0.33 and 2800 kg/m3, respectively.The metric of this experiment was the deflection of the central section A-A of the aluminum plate.The experiential results are presented in Fig.19 and Table 11,along with the simulation findings from the Chuzel’s model and developed model at 192 m/s.As can be seen, the model proposed in this study demonstrates a smaller error as compared with Chuzel’s model, and the deformation of the aluminum plate is more consistent with the experiment.

Fig.16 Comparison of the impact sequence for the experiments,model developed in this study and Chuzel’s model32 at 62.5 m/s.

Fig.17 Experimental test panel.

Fig.18 Stress–strain curve of the aluminum alloy 2014-T438.

Fig.19 Comparison of the residual shapes,RAE test at 192 m/s.

It can be observed from the comparison that the new elastic fracture model proposed in this study can obtain the relatively ideal results irrespective of a spherical ice sample hitting a thick‘rigid’plate or a cylindrical ice sample hitting a thin plate to produce larger deformation.It can be inferred that the material model conforms to the physical nature of the ice deformation and fracture to a certain extent.The model also possesses the characteristics of high calculation efficiency,few model parameters, clear physical meaning of the model parameters, parameters directly obtainable from the experiments, no adjustment of any parameters, etc., which significantly facilitate its application in aeronautical applications.

4.Discussion

The recurrence of the important experimental phenomena and experimental data is attributed to the following:

1.Based on the brittle fracture characteristics of ice at high strain rates,no consideration of plasticity makes the model simple, thus, requiring a few parameters.Specifically, only density, elastic modulus, Poisson’s ratio and tension and compression strength are required, which can be obtained through experiments.

2.The strain-rate-related constitutive equation, dynamic failure criterion and failure pressure equal to the failure strength do not present any need for adjusting any parameters in the model.

3.The asymmetrical behaviors and corresponding distinct hydrostatic pressure failure criteria in tension and compression obey the physical nature of the ice brittle fracture.

4.Once a compressive failure has occurred, the ability of ice to carry the deviatoric stress is eliminated.However, the ability of the ice to carry the hydrostatic stress is not altered, thus, allowing it to maintain load on the impact target.As a result, the ice flows like a fluid.

5.Once the tensile failure has occurred, the ability of ice to carry the deviatoric and hydrostatic stresses is eliminated.

6.Generally, the strain-based failure criterion may hardly work well with the SPH method.The failure criterion on the pressure seems to be more applicable for the SPH method to simulate the constitutional model of ice.In this study, the use of the hydrostatic pressure failure characterizes the ice failure phenomenon while ensuring the high calculation efficiency.

The material model proposed in this study can represent the elastic fracture behavior of ice at high strain rates.However,it may not be applicable to static or quasi-static simulations of elastic–plastic failure at low strain rates,which is not the main purpose of this study.

5.Conclusions

This study proposes the establishment of a high-precision macro-phenomenological ice material model that can accurately characterize the mechanical behavior of ice impact at high strain rates without adjusting any parameter with acceptable calculation time needed for the aeronautical engineering applications.The proposed elastic fracture model does not include the plastic deformation stage of ice based on the brittle behavior of ice under high strain rates.It comprehensively considers the characteristics of ice at high strain rates, including the brittle fracture, strong strain-rate-sensitivity, solid-tofluid transition, distinct tension and compression behaviors and corresponding separate failure modes.The material parameters in the model have a clear physical meaning and can be measured experimentally.Further, the SPH method has been used to simulate the ice-impact processes and quantitatively compare the performance with the previously reported elastic–plastic models as well as an elastic model with respect to the different impact speeds and ice shapes, with the results verified by the corresponding experiments.The phenomenological elastic fracture model proposed in this study is observed to capture the important data and phenomena such as impactforce history and crack propagation during ice impact.Further,the model exhibits a high calculation efficiency and a significant calculation stability.Moreover, this model does not require any manual parameter adjustment.These abilities are observed to be largely missing in the previously reported material models.Thus, the developed model overcomes the challenges faced by the previously reported models, including accurate determination of the impact load and ice fracture phenomenon as well as the need to manually adjust the material model parameters.In addition, a comprehensive comparison and analysis of the advantages and limitations of the several existing models and simulation results have been carried out under the same impact conditions.Due to the limitation of the experiments, no experimental data is available at higher strain rates, therefore, it is impossible to further verify the performance of the model under supersonic speed impact conditions.However, based on the impact speed range simulated in this study and simulation results,the developed model represents an effective tool for developing the ice-impact simulation and anti-ice impact design in aviation.Future work will be carried out by using this model for establishing the anti-ice impact performance of the aero-engine fan blades and other components.

Table 11 Comparison of the maximum plate deflection.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This study was supported by the National Science and Technology Major Project, China (No.J2019-I-0013-0013).