Peridynamic modeling and simulation of thermo-mechanical de-icing process with modified ice failure criterion

Ying Song .Shofn Li .Shui Zhng

a College of Shipbuilding Engineering.Harbin Engineering University.Harbin.China

b Department of Civil and Environmental Engineering.University of California.Berkeley.USA

Keywords: Crack growth De-icing Peridynamics Failure criteria Temperature effect Thermal mechanical coupling

ABSTRACT De-icing technology has become an increasingly important subject in numerous applications in recent years.However.the direct numerical modeling and simulation the physical process of thermomechanical deicing is limited.This work is focusing on developing a numerical model and tool to direct simulate the de-icing process in the framework of the coupled thermo-mechanical peridynamics theory.Here.we adopted the fully coupled thermo-mechanical bond-based peridynamics (TM-BB-PD)method for modeling and simulation of de-icing.Within the framework of TM-BB-PD.the ice constitutive model is established by considering the influence of the temperature difference between two material points.and a modified failure criteria is proposed.which takes into account temperature effect to predict the damage of quasi-brittle ice material.Moreover.thermal boundary condition is used to simulate the thermal load in the de-icing process.By comparing with the experimental results and the previous reported finite element modeling.our numerical model shows good agreement with the previous predictions.Based on the numerical results.we find that the developed method can not only predict crack initiation and propagation in the ice,but also predict the temperature distribution and heat conduction during the de-icing process.Furthermore.the influence of the temperature for the ice crack growth pattern is discussed accordingly.In conclusion.the coupled thermal-mechanical peridynamics formulation with modified failure criterion is capable of providing a modeling tool for engineering applications of de-icing technology.

1.Introduction

Polar regions has abundant natural resources and energy,which motivated a generation of scientists to explore [6].However.its severe climatic conditions making the polar navigation vessels and offshore platforms easily covered with ice.The main reason for the accumulation of ice on the hull and deck of the polar ships is the sea spray.Meanwhile.the cooling droplets in the atmosphere encountered a large amount of ice after melting on the surface of the ship equipment including cables.life-saving equipment.vents[11],which not only affects the vitality and the stability of the polar sailing vessels.but also causes severe damage to machines in cold working areas.including freezing capstan.cranes and valves,covering windows.blocking vents.Therefore.the de-icing technique has become increasingly necessary and important,and it has been attracted more and more attention in recent years.There are serval methods for deicing.including mechanical deicing.electric heating deicing.high-speed heat deicing.infrared deicing.ultrasonic guided deicing.chemical deicing.sacrificial coating deicing,and super hydrophobic coating deicing method.The electric heating deicing method is to melt the ice layer by using a resistance wire device such as a resistance heating element as a heating source; the high-speed heat flow deicing method is to melt and peel off the ice on the surface of the structure by spraying highpressure steam water.and the energy that is utilized can come from the extra energy produced from engine.boiler.or other equipments;infrared deicing is the area where the ice is irradiated with infrared rays of different wavelengths.and the material is dehydrated by absorbing heat energy[38].The common feature of these three deicing methods is that the ice layer of the structuredsurface is melted and peeled off by the input of thermal energy.In order to utilize the thermal deicing more effectively with low cost,the experimental tests have been conducted [5,43]to study this subject,which is the most direct and effective approach.However,the cost of the experimental tests is high.and it hardly cover all different scenarios.Numerical modeling and simulation.on the other hand.provides a cost-effective solution to study de-icing technology.Numerical modeling and simulation has now become a main approach to understand the thermal deicing mechanism,which can be applied for the basic research of the deicing technology [10].

Technically.simulating thermal de-icing procedure is a challenging task,because it is a multiple physical coupling problem that is involved in many factors in the ice-structure interaction.Moreover,the complex properties of ice material,including the grid size,temperature.and strain rate.make the ice-structure interaction difficult to model and to be captured.The ice material behavior from ductile deformation to the brittle failure depending on the different strain rates [24].Besides the strain rate effect.the most important aspect is thermal effect during the de-icing process.So far.the main approach is to use experimental data to characterize ice material properties and hence to determine the effectiveness of ice removal[9].As early as 1959,LaPlaca and Post measured the ice samples with a temperature range from -10°C to -180°C to determine the thermal coefficient expansion of ice[14].Recently.a significant discovery was made by Schulson [29]that the compressive fracture of ice occurred within the ductile-to-brittle zone [29].which is depending on the strain rate and the temperature.and he concluded that the ice fracture mode may be independent of the size scale by studying numerous freshwater ice and sea ice samples for experimental data.More recently.using experimental data acquired and accumulated in the literature for generations.Kellne et al.adopted a machine learning technique[28,42]to establish a predictive ice material model large,which is a great progress in providing ice material properties under different temperature for numerical modeling and simulation of ice material[12].Such practical tools link the numerous ice physical influence parameters from experimental data to the simulation input of the ice material model.Based on these studies,we can now identify the main material characteristics of ice.including strength.fracture toughness,fracture mechanism et al.under different temperatures.Several numerical simulation approaches were developed to model the ice material during the thermal de-icing procedure.

In the past decade,many numerical methods,e.g.Ref.[8,21,41],have been developed to simulate ice-structure interactions including the standard Finite Element Method (FEM) [25].the extended Finite Element Method (XFEM).Moving Particle Semiimplicit (MPS) method [25].the Smoothed Particle Hydrodynamics (SPH) method [4,31,32,45].and Molecular Dynamics (MD)method[1].For the numerical studies of the thermal de-icing,Xie Y et al.established a finite-element mathematical ice model to study the thermal de-icing process [39].by implementing the ice detachment model into the balance equation to acquire the thermo-mechanical behavior of ice.Bu X et al.adopted a tightcoupling method [3]to simulate the thermal de-icing process on the aircraft system,building a heat and mass transfer ice model to study the distribution of the temperature on the surface of the wing.Although many efforts have been made to improve the accuracy of simulating the thermal de-icing process.most of the studies were mainly focusing on the dynamic response of the ice[30]; there exists no appropriate numerical model to study the actual physical process of the deicing.that is the initiation and propagation of ice cracking under various thermal de-icing conditions.

Due to the limitation of the traditional numerical method [44],Peridynamics(PD)method[7]has obvious advantages to simulate the failure of the brittle material based on its non-local thermal mechanical formulation [17,18,40].which can directly treat discontinuous crack growth and propagation problems.In the past two decades,peridynamics has been developing rapidly,and it has been extensively applied to many fields.including multi-physics coupling problems.Moreover.PD can solve the dynamic thermomechanical problems by coupling the mechanical governing equation to thermal diffusion in single non-local peridynamic formulation [13].Alpay and Madenci introduced a basic theory of the fully coupled thermo-mechanical peridynamics [19,20].inwhich the thermo-mechanical constitutive model is implemented to include the linear momentum equation and the energy equation to solve coupled mechanical deformation and heat conduction problem simultaneously.Wang et al.proposed an improved coupled thermo-mechanical model within the framework of bondbased peridynamics to model the brittle fracture under the condition of thermal shock loadings [33,35].Moreover.the thermal mechanical failure of the brittle materials.such as cracking initiation and energy propagation under the condition of power cycles or heating.has been simulated based on the coupled thermomechanical peridynamics [34,37].To treated the wave reflection and analyze the crack pattern distribution.Ren at al [26,27].have developed an inhomogeneous particle discretization approach to improve the stability of PD algorithm to accurately capture crack initiation and propagation.Thus.peridynamics can provide a suitable tool for direct modeling and simulation thermal de-icing problem.which not only provides an accurate brittle ice material model but also is capable of capturing the ice cracking process.

In this work.we are focusing on applying the developed thermal-mechanical peridynamics method to model and simulate the thermal de-icing process.Based on the theory of the coupled thermo-mechanical bond-based peridynamics,we have developed a PD de-icing model that is able to perform direct simulation of thermal de-icing process.

This paper is organized into five parts.In Section2.we first reviewed the coupled bond-based peridynamics method,and then,we introduced the material properties and physical characteristics of ice material.Subsequently.we proposed a modified failure criterion for ice fracture taking into the temperature effect.At last,we introduce a set of thermo-mechanical boundary conditions that can be efficiently used in the de-icing simulation.In Section 3.we present the detailed numerical solution procedures.We first introduced the numerical discretization procedure.and then determined the numerical stability conditions accordingly.Finally,we discussed the surface effect on the de-icing process.In Section 4,several numerical examples under different conditions are presented,and some of them are compared with experimental results or the results of finite element modeling and simulation.In Section 5.we concluded the study with a few conclusions.

2.Peridynamic model for de-icing system

To build an appropriate thermo-mechanical model for ice materials under the condition of thermal de-icing.The theory of the coupled thermo-mechanical bond-based peridynamics is adopted and introduced.

2.1.Review of the coupled bond-based peridynamic thermomechanics

The basic assumption of the non-local peridynamics is the material points in the reference configuration Ω0only interact with the material points within a certain radius,which is called the horizon.The interaction between the material points only affective within the horizon.As shown in Fig.1(a),material points xiexert its bond force on the material points xjthrough the bond ξij.which can be denoted as

When the deformation occurs.the bond connected the two material points stretches accordingly.and the relative displacement vector ηijcan be expressed as follows,

Fig.1.Thermal bonds in thermo-mechanical bond-based peridynamics.

In the bond-based peridynamics.the bond force among different particles.f(ξij，ηij，θij，t).is assumed to be pairwise and equal at the current time t[35],and it can be written in general as follows,

The non-local equation of motion can be then expressed as,

where ρ0is the density of the material points xi.u(xi，t) is the displacement vector.b(xi，t) is the external body force density.

The material points within the reference body region is connected each other through the thermal conductors.which can be called as “thermal bonds”.or “t-bond” (see Fig.1) in short [2].The heat flux between each material points xiand xj.which is depending on the relative position ξij.and the temperature difference θijbetween the material points xiand xj,the heat flux can be presented as (see Ref.[33]),

in which κ(xi，xj) is the microscale thermal conductivity of the material,is the temperature difference between the two material points xiand xj.andis the unit vector along with the thermal bond.According to Bobaru"s studies.the heat flux transfer along the direction e between the material points xiand the family material points xjwithin a certain radius horizon δ.The “t-bond” is assumed to be insulated.which means the heat transfer is not existed within the thermal bonds.Based on the energy law.the heat transfer in the peridynamic formulation in the direction e can be defined as (see Ref.[16]),

where ρ(xi) is the density constant.cv(xi) is the specific heat capacity constant.is the average temperature on the thermal bond,which connects the material points xiand xj.Dividing Eq.(7)by ξij·e and we can get the integration as follows:

where κ(xi，xj) is the micro thermal conductivity.

Assuming that the temperature Θi(xi，t) at material points xiand the average temperaturefollows the relationship given below:

We can then determine the micro thermal conductivity by equating its equivalent average with the macros thermal conductivity (see Ref.[33]),

In 3D cases.the micro-conductivity can be linked to the macro thermal conductivity as:

where k is the macro thermal conductivity of the ice material.

Thus.we can obtain the heat conduction equation at the material point xiin the bond-based peridynamic formulation as follows (see Ref.[34]):

where VHxiis the volume of the material points xiwith its neighboring region.The heat flux under the temperature change can be defined as a pairwise thermal bond force density as follows [37],

Then.the bond-based peridynamic equation with thermal conductivity can be written as:

By taking into account the heat source,the thermal-mechanical bond-based peridynamic formulation becomes:

where hs(xi，t) = ρ0sbis the heat generation rate.

2.2.Ice constitutive model in peridynamics

2.2.1.Ice accumulation and de-icing

The accumulation of ice on the hull and deck of the polar ships is mainly due to the sea spray [23].As shown in Fig.2.the cooling droplets in the atmosphere encountered a large amount of ice after meting the surface of the polar equipment.including cables.lifesaving equipment.vents under the condition of the freezing temperature (below -2°C).The cold atmosphere brings the droplets into contact with the structure formed and built-up ice and theliquid membrane.with the accumulation of the ice.the sea salt sediment,and the ice deposits into pure ice and saline.Fig.2 shows a schematic illustration of the phenomenon of the ice growth process on the surface of the structure.

Fig.2.The heat transfer of icing on the surface of the structure [24].

There are many techniques to detect ice accumulation,and such work is ongoing.The phenomenon of the ice accumulation is complex,which mainly caused by the multiple physical properties of ice,including the thermal properties.The physical parameters of ice are affected by temperature.icing strength.water droplet diameter.water composition.and perimeter factors such as ambient air velocity.As the temperature of the ice layer increases,the adhesion strength of the ice will increase.As it increases.the adhesion strength of ice decreases.The shear strength of ice is lower than the tensile strength.and the adhesion strength of impact ice is stronger than that of static ice.The freezing strength of the frozen ice is low.The elastic modulus of ice is also a constantly changing quantity.It is necessary to choose the appropriate physical properties of ice for the de-icing simulation.The thermal properties of ice in this paper are shown in Table 1 according to the previous reference studies.

What"s more,since the previous experimental results show that the adhesion strength between the ice layer and the aluminum plate depending on the icing conditions[15],the shape of the icing form under different temperature conditions has been studied.Fig.3 shows several types of ice under typical conditions.The icing environment in Fig.3(a) is: the initial freezing temperature of the aluminum plate is-20°C,the diameter of the water droplets is 48 μm.and the ice shape is frost; Fig.3(b) shows the initial freezing temperature of the aluminum plate at -15°C.the diameter of the water droplets is 58 μm.and the ice shape is lumpy and opaque.When the temperature is higher than -10°C.ice is formed as a gravel-like shape.the ice layer firmly adheres to the aluminum plate.which is delicate to completely clean under the shock of electrical pulses.and the electric pulse deicing power needs to be increased to remove the ice layer.Thick ice layers are more accessible to be removed than thin ice layers within the range of 6mm.Because thick ice layers tend to be more brittle to break.

2.2.2.Constitutive model of ice

Ice is a quasi-brittle materials.and the commonly used peridynamics ice model is based on the bond-based peridynamics,which is a PMB(Prototype Micro-elastic Brittle)material model.For a PMB peridynamics model.the pairwise force density function can be defined as followings:

in which g(ξij，ηij，θij，t) is a scalar-valued function of the micro modulus c.thermal expansion coefficient α.and the temperature difference θijbetween the material points.It can be written as:

in which s is the bond stretch,and c is the micro-modulus,and they can be determined or expressed as follows:

where E is Young"s modulus,δ is the size of the horizon.The bond stretch represents the deformation state or strain state between two material points.and in peridynamics when the bond stretch reaches to a critical value,the mechanical bond between two points will break and carried no force.as is shown in Fig.1(b).

In Peridynamics formulation.there is a history-dependent function μ(ξij，ηij，θij，t) that can represent the damage state of the bond,

In this work.to build an appropriate constitutive model for quasi-brittle material ice.an elastic-brittle damage model is adopted based on the quasi-brittle model proposed by Wang et al.[36].As shown in Fig.4.the history-dependent function in the elastic damage model can be expressed as:

where s0is the critical bond critical stretch.γ∈(0，1)is the “strainsoftening” coefficient [36].Then by using Eq.(22).we can replace the function density function f(ξij，θij，t) as follows (in one dimensional):

The comparison between f and fhis displayed in Fig.4.

Similarly,we can also replace the thermal force density(Eq.(14)as:

It may be noted that the scalar functions μ and μhare damage functions for mechanical field and thermal field.As shown in Eqs.(21)and(22),they are not the same,but their cutoffs are both set at s = s0.When a mechanical bond between two ice particles is stretching,the material in between may become a mushy state,but it is still capable of heat conduction until it is broken.Thus the damage process is not an abrupt even but a gradually depletion process.which we use a linear “softening” process to represent it.This implies that we still have heat conduction between two ice particles during this stage.When the mechanical bond betweentwo ice particles is broken.the material space in between may be occupied by air.or snow.or water.At this stage.even if the mechanical bond between the two has been broken,there is still a heat conduction in this complex state of a mixture medium.However,one may find from Table 1 that the thermal conductivities for water,snow.and air are almost one order of magnitude less than that of ice,and a good approximation is completely ignore them after the mechanical bond is broken.That is the reason or justification that we can set the critical stretch for both mechanical bond and thermal bond at s=s0(see Fig.4).

Table 1 Thermal properties of air.fresh water.and ice at 0°C[14,16,22].

Fig.3.Ice formations under different temperature conditions (From Ref.[15]).

Fig.4.Damage models for mechanical bond and thermal bond.

2.3.Modified ice failure criteria with temperature difference

In peridynamics.the critical bond stretch s0is used to predict the damage of the bond.The traditional method to calculate s0is using the energy release rate G0.and can be expressed as followings in generally:

where K is the bulk modulus.G0the energy release rate of the ice material.and δ is radius of the horizon.

To model the elastic-brittle failure of ice.G0can reflect the fracture resistance of the ice material,the primary failure mode for ice is usually under tensile loading.because ice exhibits strong compressive strength than the tensile strength.thus the energy release.G0.can be obtained by the fracture toughness KIvia G0=/E .and subsequently the critical bond stretch can be written as:

The fracture toughness KIcan be obtained by the experimental tests as is shown in Table 1.

Considering the thermal effect,we proposed a modified critical stretchto replace s0as a new failure criterion to predict the thermo-mechanical damage of the ice material.and the modified bond critical stretch is defined as follows:

in which.Θ0is the initial temperature condition.Θcis a critical temperature difference.and Θ is the current temperaturedifference at the material particle under the consideration.As Θ?Θc.the temperature-dependent critical stretchwill decreases rapidly indicating the melting ice is much easier to fracture than deeply frozen ice.In the computation.we choose Θ0= - 25°C and Θc= 0.We introduce a local damage index at each material particle.which is dependent on the number of the broken bonds between the material point with other material points within the horizon.and it can be expressed as:

and

Finally.the fully coupled bond-based thermo-mechanical peridynamic equations can be written as:

and

2.4.Boundary conditions in peridynamic de-icing model

2.4.1.Displacement and velocity

A virtual region ?mwith depth δ is introduced as the boundary region.The initial displacement vector U0and the initial velocity vector V0are exerted on material points xiin the virtual region ?mas follows respectively,

and

2.4.2.Body force density

For a distributed pressure P(xi，t)or a concentrated force P(t)on the boundary region ?b,the vector of the body force density can be calculated as:

where Δ is the thickness of the nonlocal boundary.

2.4.3.Temperature (nonlocal dirichlet boundary)

The nonlocal Dirichlet boundary condition can be imposed in a virtual region ?tas follows:

where y*is the position of the material point in the actual region?,y is the position of the material point in the virtual region ?t.x is the position of the material point on the boundary surface ?t.y*is the location of the mirror image of y.which can be obtained by solving the equation,y*=y+2|x -y|n,n =(x -y)/|x -y|.In the case of Θ*(x，tTH) =0,the temperature in the virtual region?tturns to be the negative value of the temperature in the actual region ?.

2.4.4.Heat flux (nonlocal neumann boundary)

Constant heat flux q is imposed in terms of the heat source hsalong the surface ?fof the boundary region ?f.which can be written as:

where Δis the distance of the boundary region usually taken as the horizon size δ.

3.Numerical solutions

3.1.Numerical procedure for de-icing system

The fully coupled bond-based thermos-mechanical peridynamic equation of motion can be discretized numerically as follows:

and

where i is the index of the material point;j is the index of the family material points within the horizon.and V(j)is the volume of the material points x(i).

Coupled with the each other.the displacement field and temperature field can be updated as:

and

The stretch between the material points x(i)and x(j)at the time step n can be discretized as:

and the temperature difference in the reference configuration at the time step n can be discretized as:

Thus.we can obtain the velocity.the displacement and the temperature at the material points x(i)at the time step n+ 1 expressed as the following algorithm:

where Δt is time step increment.ΔtTHis the thermal time step increment.The fully coupled thermo-mechanical interaction of the material points is illustrated in Fig.5.The discretized equation of the linear scalar-valued function and the thermal response function can be prescribed as:

3.2.Numerical stability conditions

For thermo-mechanical numerical simulations,it is necessary to determine the size of the thermal time step ΔtTHto prevent the diffusion of the convergence from the numerical solutions.The condition of stability is defined according to the von Neumann method,and the temperature at each time step ΔtTHis assumed to obey the following equation,

where Γ is a positive real number,which represents the number of the transferred wave,and ζ is a complex number obey the equation of |ζ|≤ 1 for every Γ.By omitting the heat source.the discretized peridynamic equation with heat diffusion can be rewritten as,

Substituting Eq.(49) to Eq.(48).we can obtain the thermal diffusion equation as follows:

For simplification.we define MΓ as,

The complex number ζ is valid for every transferred wave under the condition of |ζ|≤ 1.and can be cast as followings:

Therefore,we can obtain the time step size ΔtTHrestricted to the following equation:

where MΓis defined in Eq.(52).and it must satisfy the following inequality,

Fig.5.Coupled thermo-mechanical interaction between the discretized material points.

By substituting Eq.(55) to Eq.(54).we can finally obtain the time stability condition yields the following criterion:

where ΔtTHis the stable time step for the heat transfer procedure.In addition.a stability condition for the mechanical stability condition yields the following criterion:

where ΔtMEis the stable time step for the time integration of mechanical field.

The general numerical procedure or computational flowchart is illustrated in Fig.6.

3.3.Treatment of surface effect on ice-aluminum interface and volume correction

The value of the thermal response function fhis related to the micro-conductivity κ,which is calculated by the thermal potential z within the entire domain of the material points x(i).Thus.the parameter of the micro-conductivity κ is determined by the horizon of the numerical integration domain.If the material points x(i)are located near to the surface or the interface between the two materials,the correction of the parameter of the micro-conductivity is required,as is shown in Fig.7.

The material points x(i)interact with the material point x(j)and material point x(k).x(j)is the material points within the region of Ω1.which presents the ice material,and x(k)is the material points within the region of Ω2.which presents the aluminum material.Their micro-conductivity κ(i)(j).can be presented as an equivalent thermal conductivity k(i)(k)for the material points x(i)and x(k),

where lIceand lAlare expressed the distance between the material points x(i)and x(j)within the region of ice material.and the distance between the material points x(i)and x(k)within the region of aluminum material.respectively.As is illustrated in Fig.7.kIceand kAlare the thermal conductivity of ice material and aluminum material.respectively.As is shown in Table 2.

The thermal micro-potential z(i)(j)between the two material points x(i)and x(j)can be corrected as,

The surface and the volume correction factor vcis denoted as:

where r is the line segment cut by the boundary along the radial direction of the horizon.Meanwhile,the thermal response function fnT(i)(j) between the material points x(i)and x(j).fnT(i)(k)between the material points x(i)and x(k).can be modified accordingly as,

and

Fig.6.TM-BB-PD numerical procedure for de-icing simulation.

Fig.7.Surface effect and volume correction on the ice-aluminum interface.

Table 2 The physical properties of the numerical model.

4.Numerical simulation for de-icing system

In order to demonstrate the accuracy and feasibility of the coupled TM-BB-PD method for modeling and simulation of thethermo-mechanical ice behavior during the de-icing procedure.We conducted four numerical examples in this section.

Firstly.in order to validate and demonstrate the numerical deicing model.we investigate the experimental results from the open literature,and compared it with our numerical simulation in subsection 4.1.Secondly.we investigate the convergence of the numerical model in subsection 4.2.including the δ-convergence and the m-convergence.the effect of the factor m and the size of δ are discussed.Thirdly,we implemented the de-icing model under several different loading conditions to predict the crack propagation and the temperature distribution during the coupled thermomechanical de-icing process in subsection 4.3.Finally.the temperature effect is discussed in subsection 4.4.

4.1.Validation of the numerical model for de-icing system

In this section.we focused on validating of Peridynamics deicing model and studying of accuracy and practicability for the de-icing model.

The de-icing technology is wildly used and many de-icing methods have been developed.One of the most popular de-icing technique is the electro-impulse de-icing (EIDI).Many researches including both experimental and the analytical studies have been conducted for modeling and simulation of EIDI systems.The basic principle of the EIDI system is that the electric pulse coil placed inside aluminum plate will generate a strong pulse current in a short time (on the order of milliseconds).The impulse will cause the aluminum skin to vibrate with small amplitude and high acceleration,and the ice layer covered on the aluminum layer will be broken and loosen by the airflow.In the de-icing thermal-mechanical formulation,this is equivalently to impose a concentrated load and a fixed temperature jump at the specific location of ice layer where the electro impulse is applied.

In this work.we employed the coupled thermo-mechanical peridynamics model to study a EIDI system.that is the same EIDI system studied in Ref.[15].We used the same dimension and equivalent load conditions for the EIDI model as that in Ref.[15].To determine the loading case for numerical simulations to represent the electro-impulse de-icing process,the tensile strength of tensile adhesion strength of ice is needed.Chu et al.shows that when the ice layer is very thick,it takes a very strong force to peel off the ice from the aluminum plate.because the adhesive strength between the ice and the aluminum plate is much greater than the bond strength of ice.Kandagal et al.demonstrated that the tensile adhesion strength was set to be 10MPa as a criterion to judge ice stripping process.

Fig.8.Geometrical model for de-icing system.

The peridynamics de-icing experiment set-up is illustrated in Fig.8.The de-icing model contains two interacted layers: the ice layer is on the top of an aluminum layer.The dimension of the numerical model is set as 0.90m× 0.30m× 0.006m.the thickness of the ice layer is chosen to be 0.004mm on the top of the aluminum layer.and the smallest distance between two particles or the grid size is Δx= Δy= Δz=2mm with uniform particle distribution.The total particle number is 202500.and the element number is 133802.The length of the horizon is δ=3.015Δx The material properties of the numerical model are shown in Table 2.The mechanical time step is chosen to be Δt =1.0× 10-8.The four edges of the plate are fixed,and the bottom center region of the plate,i.e.and r =0.03m is subjected to a sinusoidal pressure wave loading.The pressure applied on the loading zone is chosen as P(t)=2381sin(2433πt)combined with the temperature difference of 0.5°C to represent the electro-impulse load studied in Ref.[15].The total time step is 2500.The simulation is conducted with a number of threads 64(N)in a CPU model of Intel(R)Core(TM)i7-6500U CPU of main frequency 2.5 GHz.The fortran90 code is compiled by using the OpenMP-based Intel Parallel Studio XE 2017.The total simulation time is around 2.5 h.

The initial conditions of the problem are:

where Θ0= 0°C.

The mechanical boundary conditions are:

where ΔΘ = 0.5°C and P(t) = 2381sin(2433πt).The pressure profile used in this study is shown in Fig.9.

Fig.10 presents the comparison results between (a) the peridynamics modeling simulation.(b) the experimental observation,and (c) the finite element modeling and simulation.From Fig.10,one may find that the numerical prediction show good agreement with the experimental observation and FEM simulation [15].

In order to analyze the de-icing efficiency.the definition of deicing rate is introduced.the ratio of the cover area is expressed as a formulation

Fig.9.Loading amplitude plot.

where d is the de-icing rate,Aris the ice removal area,and A0is the total area covered by ice.

In the numerical analysis.it is not difficult to find the area covered by the remaining ice particles in the model and the area covered by the original ice layer.Based on these data.we can calculate the de-icing rate,and we find that the de-icing rate in the present simulation example is able to remove 38.6%of the ice layer from the aluminum plate.

Fig.11 presents the simulation results of EIDI de-icing process by using TM-BB-PD de-icing model.We can clearly obtained the damage pattern of the ice crack initiation and propagation during the de-icing process.which can not be easily catched by FEM simulations.

Fig.10.Comparison of the numerical result to the previous study.

Fig.11.TM-BB-PD simulation of EIDI de-icing process.

4.2.Convergence studies for numerical stability

In this section.we conducted a numerical study on the convergence of the peridynamics simulation.The de-icing model in the numerical simulation is chosen to be a square shaped plate with different sizes or dimensions and different loading conditions.The de-icing model again contains two interacted layers.the ice layer on top and the aluminum layer at bottom,as shown in Fig.12.The physical properties of the material including thermal properties and mechanical properties are shown in Table 2.

In the convergence study.four numerical tests have been conducted to investigate convergence.stability.and accuracy of the TM-BB-PD de-icing numerical model.The numerical specimen has a dimension of 0.02m× 0.02m× 0.006m.The detailed parameters in the convergence examples are shown in Table 3.including the particle distribution and the grid size.For δ-convergence study(examples 1,2).the nonlocal ratio m is fixed as m = 3,the horizon size δ is chosen as 1.5mm and 3.0mm.respectively.For m-convergence study (examples 3,4).the size of the horizon δ is fixed as 3mm,and ratio m is chosen as 12 and 15,respectively.The thermal loading condition in the numerical tests is prescribed temperature difference in each thermal bond as Δθij= 0.5°C.In the numerical simulations.we first assumed the boundary sets of the plate to be free,and the effect of the grid size of ice and aluminum surface on the verticle displacement uzis investigated to test the accuracy and stability of the simulation.The vertical displacement uzalong the interface of the ice and the aluminum with different grid size are shown in Fig.13.As is shown in Fig.13.the constant temperature difference loading causes a bending along the interface between the ice material and the aluminum material.The grid size has some effects on the peak displacement value at the center of the plate,which is in agreement between the FEM simulation conducted by using Abaqus[15].Based on the convergence study,we recommend Example 3 (δ = 3，m = 12) as the optimal grid size.

In the convergence study.we set the displacement and the velocity boundary conditions along the four edges of the plate to be zero.and we assumed that the initial temperature of the plate is(0°C).The distribution of the temperature and the damage contourson the ice-aluminum plate are shown in Figs.14-15.In the numerical calculations.the mechanical time step is chosen to be ΔtME= 1.0× 10-9s.and the thermal time step is chosen to be ΔtTH=1.0×10-12s,respectively.The total time step is 10000.The simulation is conducted with a number of threads 64(N) in a CPU model of Intel(R) Core(TM) i7-6500U CPU of main frequency 2.5 GHz.The fortran90 code is compiled by using the Intel Parallel Studio XE 2017.The total simulation time is around 1 h for Example 1; 3 h for Example 2; 11 h for Example 3.and 14 h for Example 4.

Fig.12.Illustration of the numerical specimen for de-icing experiment.

Table 3 Discretized details of the numerical convergence studies.

The comparison of the temperature distribution at time tTH=1.0×10-9s are shown in Fig.14.As seen from Fig.14.the temperature distribution is sensitive to δ and nonlocal ratio m.it can be obtained under the condition of m=15(δ =3)most clarify.

The damage pattern of ice layer for different grid size at time t =1.0×10-5s is plotted and compared in Fig.15.From Fig.15,we can find that with the increasing number of particle.the crack pattern of the ice block can be represent more clearly.

Fig.13.Comparison of the verticle displacement with different grid size.

4.3.Ice crack propagation in the de-icing procedure

In this section,a de-icing TM-BB-PD model is built to verify the accuracy and practicability for the PD simulation of thermomechanical ice removal process under the combined mechanical and thermal conditions.In this example.the dimension of the numerical model is set as:0.02m×0.02m×0.006m,the thickness of the ice is 0.004m on the top of the aluminum plate.and the thickness of the bottom aluminum plate is 0.002m.The grid size is Δx=Δy= Δz=0.2mm with uniform particle number,as shown in Fig.12.The horizon size is δ=3.015Δx.The properties of the numerical simulation is shown in Table 2.In the numerical test.we assume the steady-state heat conduction.Thus,we assume that the temperature difference in each bond θijis 0.5°C outward.The mechanical time step is chosen to be Δt = 1.0× 10-9s.and the thermal time step is chosen to be ΔtTH=1.0×10-12s,respectively.The model is fixed along the four edge of the plate,and it also is subject to a sinusoidal wave pressure at the bottom center region of the plate.i.e.and r = 0.003m.

The initial mechanical combined with thermal conditions are:

where Θ0= 0°C.

The mechanical boundary conditions are:

Fig.14.Comparison of the results of temperature distribution at time.tTH = 1.0× 10-9s

Fig.15.Convergence studies of the damage pattern (top view) at time.t = 1.0× 10-5s

where P(t) =12sin(1000πt),T(t) =1× 106t.The pressure used in the present study is shown in Fig.9.

In this section,we investigate the de-icing process by predicting its temperature distribution and the damage pattern.The temperature distribution and the damage distribution are plotted at different time steps are shown in Figs.16 and 17.From Figs.16 and 17.we can obtain the pattern of the crack wave initiation and propagation.As can be seen in Fig.16,the value of the temperature increasing with the time.The temperature contours and the damage patterns have the same distribution which can be seen obviously in Figs.16(f) and Fig.17(f).The phenomenon reflect the thermal energy to some extent affects the total energy of the material,which causes the damage and crack.In the view of the energy wave propagation.we can obtain the results from the cracking pattern.As is shown in Figs.16 and 17.due to the zero constraint forced by the boundary edge of the plate,the crack initiated around the edge of the model.the boundary edge is under the tensile stress.and the center of the plate is imposed to the compressive stress,because the specialty of the ice material which have stronger compressive strength than the tensile strength.therefore.the initial crack occurred on the boundary of the ice plate.What"s more,with the time increase,the energy wave transfer from the boundary to the center of the loading zone.the crack length increasing correspondingly around the center of the plate.The maximum temperature occurred on the central area of the plate.Which represent the detail of the damage pattern under the certain fully coupled mechanical and thermal loading pressure conditions.

Fig.16.Temperature contours.

The wave speed in the ice can be calculated as:

in which we can obtain the wave speed by using the material properties listed in Table 2.Based on Fig.17,the crack propagation speed is around 0.007/6.5×10-6=1076 m/s.which is less than the sound speed 2476 m/s in the ice obtained by Eq.(77).

4.4.Effect of the temperature change

To investigate the temperature effect on the thermo-mechanical de-icing process,we conducted five numerical tests under different temperature conditions.The loading applied on the specimen is a uniform temperature change.chosen to be 10°C.20°C.30°C.40°C,and 50°C.respectively.The boundary conditions at four edges are set to be traction free.To make the numerical results more prominent.we set the dimension of the specimen as 0.03m× 0.004m×0.001m.as shown in Fig.18.The ice material at the upper part of the plate (top view) with the width of 0.003m.and the aluminum material at the lower of the plate (top view) with the width of 0.001m.The properties of the material are shown in Table 2.The ice layer and the aluminum layer have the same length and depth,but the two materials have different thermal properties and differentwidth.The properties of the interface between the two materials is the same as the ice material.The total number of the material points is 12000,the horizon size is δ =3.015Δ.We use the adaptive dynamic relaxation algorithm in the computation with Δt = 1.0s,and the total number of time steps is 20000.

Fig.17.Ice crack initiation and propagation in de-icing procedure.

The results of the displacement field of the plate under the different thermal conditions are shown in Fig.19.From Fig.19(b),we can find a bending occurred on the interface of the ice and aluminum material in x directions.which because of the thermal expansion coefficients of the ice material and the aluminummaterial are different.With the increasing of the temperature,this phenomenon exhibit more obviously.which to some extent demonstrates the sensitivity of the temperature factor on the deformation of the ice-aluminum plate.

Fig.18.Schematic illustration of Ice-aluminum plate subjected to temperature change (top view).

Fig.19.Effect of the temperature change on the displacement distribution.

To further investigate the temperature effect on the ice crack initiation and propagation.For simplicity.we chose a cantilever beam model,the displacement and the velocity are set to be zero on the left side and a sine wave load is imposed on the right side of the beam,as is shown in Fig.20.The dimension of the beam is 0.12m×0.03m× 0.006m,the grid size is Δx =0.01m,and the horizon size is δ =3.015Δx.We conducted two constant temperature change of 8°C and 10°C,respectively.The mechanical time step is chosen to be ΔtME= 1.0× 10-9s.

The results of the damage pattern under different temperature conditions at time t =3.5×10-4are shown in Fig.21.As seen from Fig.21,the damage pattern is sensitive to the temperature factor θij.Under the prescribed coupled pressure load f(t) and temperature change θij,the ice crack initiate and propagating along the straight line of the beam.with the increasing of the temperature.the damage area expanded rapidly.

5.Conclusions

In this work.a three dimensional fully coupled thermal mechanical bond-based peridynamics (TM-BB-PD) method was used to model and simulate the thermo-mechanical de-icing process.The ice constitutive model is established by considering the influence of the temperature effect in the framework of TM-BB-PD,and a novel modified failure criteria including the temperature factor is proposed to predict the damage of quasi-brittle ice material.A nonlocal thermal boundary condition is developed to simulate the thermal loading of the thermal de-icing process.From the numerical results.we can draw the following conclusions:

· The thermo-mechanical peridynamics model can efficiently predict the de-icing process.and the Peridynamic results show good agreement with the experimental results;

· The peridynamics δ-convergence and m-convergence do have an effect on the thermo-mechanical fracture of the ice material,not only in terms of temperature distribution but also in terms of the damage or crack pattern.The ratio size m = 12.δ=3 is recommended in de-icing modeling;

· The developed Peridynamics approach can capture the initiation and propagation of the ice crack as well as the distribution of the temperature in the ice.The temperature distribution or contour and the damage distribution or crack patterns have a similar pattern of the distribution as time increases;

· The temperature change have a significant effect on the mechanical behavior of the ice plate.Both the displacements in the x-direction and the deformation in the y-direction of the ice plate increase with the temperature change.The damage pattern is sensitive to the temperature change.

In conclusion.the modified fully coupled thermo-mechanical bond-based peridynamics (TM-BB-PD) method can provide an efficient tool to model de-icing process on frozen structures.

Fig.20.Schematic of the cantilever beam model.

Fig.21.Temperature effect on the damage pattern.

Declaration of competing interests

The authors declare that they have no conflict of interests.

Acknowledgements

This work is performed at the University of California at Berkeley.Ms.Y.Song gratefully acknowledges the financial support from the Chinese Scholar Council (CSC Grant No.201706680094).