The Mechanical Mechanism and Inﬂuencing Factors of Ice Adhesion Strength on Ice-Phobic Coating

: Ice accretion can cause problems on polar ships, ocean platforms, and in other marine industries. It is important to understand the interface debonding behavior between ice and the surface of equipment. In this work, we created a mechanical model to analyze the interface debonding behavior between a square-based ice cuboid and an elastic coating base, using contact mechanics and fracture mechanics. Three-dimensional (3D) ﬁnite element (FE) simulation was used to simulate the interface debonding for normal and shear separation. A bilinear cohesive zone model (CZM) was used to simulate the interface between the ice cuboid and the elastic coating. We investigated the effect of the elastic modulus E of an elastic ﬁlm on the critical detachment force F c for normal and shear separation. The results showed that F c increases with an increase of the elastic modulus of the elastic ﬁlm. When E exceeds a certain level, F c achieves a constant value and then remains stable. Finally, a series of epoxy/polydimethylsiloxane (PDMS) interpenetrating polymer-network (IPN) gel coatings with different elastic moduli were prepared. The ice tensile and shear adhesion strengths ( σ ice and τ ice ) of the coatings were measured. The results were roughly consistent with the results of the numerical simulation when E < 1 MPa.


Introduction
As the climate warms, new shipping routes open in the polar regions. In these regions, sea spray icing or atmospheric icing cause economic as well as safety problems, especially on polar ships and ocean platforms [1,2]. Conventional ice-phobic coatings are achieved mainly by using lubricating materials and by surface microstructural modification (e.g., slippery liquid-infused surfaces [3] and super-hydrophobic coating [4]). Inspired by the mechanics of fracture, researchers have been interested in soft materials with a low elastic modulus for ice-phobic application in recent years [5][6][7][8][9]. For these soft-elastic-coating applications, it is important to understand the interface debonding behavior between ice and the elastic coating. Due to the huge difference in elastic modulus between the soft coating and the ice, the ice can be regarded as a rigid body. Interfacial debonding of a rigid object from an elastic film through normal separation has been studied by many researchers. Based on Griffith's energy criterion, Kendall, [10] decades ago, derived the critical detachment force F c required to separate a rigid flat cylindrical punch from an elastic film for normal separation. It is where r is the radius of the cylindrical punch; W ad is the adhesion work; and E, v, and t are elastic modulus, Poisson's ratio, and thickness of the film, respectively. This formula gives the proportionality F c ∝ W ad E t , a relationship that was also derived by other researchers for normal loading [11][12][13]. A few researchers also investigated shear separation. Chaudhury [14] demonstrated that σ s ∼ a l W ad G t in the case of a rigid cuboid in contact with a thin film under shear loading, where σ s is the critical shear stress, l is the distance from the loading point to the cuboid-film interface, G is the shear modulus of the film, and a is the size of the interface length. These studies have shown that the scaling law F c ∝ W ad E t is always obeyed, for both normal and shear loading. According to this scaling law, when t decreases or E increases without limit, F c grows toward infinity. However, in actual situations, adhesion energy W ad and interface adhesion strength σ max are limited, preventing the infinite growth of F c . Some researchers [15,16] have shown that when the thickness of the film is small enough, F c remains nearly constant. Interfacial debonding behavior is governed by t, E, v, W ad , σ max , etc. For soft ice-phobic coating, effects of the elastic modulus on interfacial debonding behavior is worthy of attention. In addition, although most existing research is concerned with normal loading, for ice-phobic coating applications, the ice is always under shear loading-such as with wing icing, windturbine-blade icing, Ref. [17] ice adhesion strength test, Ref. [18] and vertical plane icing. The ice under those conditions will be under shear loading from the airflow aerodynamic force, the centrifugal force, and gravity.
In the work reported in this paper, theoretical and numerical methods were used to analyze the interface debonding behavior between an ice cuboid and elastic coating. First, a mechanical model was created to analyze the critical detachment force F c required to separate an ice cuboid from an elastic coating under normal loading based on contact mechanics and fracture mechanics. Then, three-dimensional (3D) finite element (FE) simulation was used to analyze the interfacial debonding behavior of an ice cuboid from an elastic coating when it is bonded to a rigid substrate under normal and shear loading by a bilinear cohesive zone model (CZM). The effect of E on F c was studied. Finally, we prepared a series of epoxy/polydimethylsiloxane (PDMS) interpenetrating polymer network (IPN) gel coatings. Silicone oil can reduce the crosslink density of the IPN gel. By adjusting the additive amount of silicone oil, the elastic modulus of coatings can be reduced from hundreds of MPa to less than 1 MPa. The elastic moduli of those coatings were measured. In addition, the ice tensile and shear adhesion strength (σ ice and τ ice ) were measured using a force gauge and a centrifugal adhesion test at −10 • C. The experimental results and the simulation results were compared.

Problem Description
We investigated the interface debonding of an ice cuboid (square-based, size: a × a × 2l) from an elastic coating which is perfectly bonded to a rigid substrate under normal and shear loading as shown in Figure 1a,b. The ice cuboid can be regarded as a rigid body because of the large difference between the elastic moduli of ice and the elastic coating. The thicknesses, elastic modulus, and Poison ratio of the coating are, respectively, t, E, and v. The adhesion energy of the interface is W ad ; this is the energy per unit area required to separate the ice cuboid and coating. The normal or shear concentrated force is applied to the centroid of the ice cuboid. Figure 1. Schematics of (a) a rigid square-based ice cuboid adhered to an elastic coating that is bonded to a rigid substrate; (b) the ice cuboid under normal and shear loading. E: elastic modulus; v: Poisson's ratio; l: the distance from the loading point to the cuboid-film interface; a: the size of the interface length.

Rigid Cuboid Adhered to an Elastic Half-Space
First, we consider a concentrated normal force P applied to an elastic half-space as shown in Figure 2a. M is an arbitrary point on the surface of the elastic half-space. Using contact mechanics, the vertical displacement u z of point M due to the point force P is [19] where ρ is the distance from P to M. Utilizing the principle of superposition, u z of point M due to the whole contact area A can be obtained by integration. First, distribution pressure at the contact area caused by F can be expressed as p(x, y). By using Equation (2), the displacement u z of point M due to p(x, y) at an area element dA = dξdη (see Figure 2b,c) is In addition, because the ice cuboid is a flat rigid body, u z of point M in the contact area A has a constant value d 0 (see Figure 2b), and one gets The equilibrium equation is p(x, y) = 0 , (x, y) / ∈ A The problem evolves into solving integral Equations (4) and (5). with p(x, y) and d 0 being the unknowns. If the contact area A is circular, p(x, y) has an analytical solution, and d 0 can also be found [19]: where r is the radius of the circular contact area and d 0 is equal to Kendall's result [10]. However, for the square contact area treated in this paper, p(x, y) does not have an analytical solution; but it can be found by a semi-inverse method [20] or orthogonal polynomials method [21]. The quantity d 0 has a similar form to that of Equation (8); it is where k depends on the method used to solve for p(x, y). The quantity k is found to be 0.4265 [20] or 0.4363 [21]. From Equation (9) the stored strain energy U E is where C 1 is a constant of integration. The potential energy U P of the external force is Using Griffith's energy equilibrium theory, the energy release rate G c is When Equation (12) is satisfied, interfacial debonding occurs. Using dA = da 2 = 2ada, one obtains from Equations (10) and (11): The critical detachment force F c is found to be Equation (14) is the solution for a rigid flat square-based cuboid adhered to an elastic half-space for t/a 1.

Rigid Cuboid Adhered to an Elastic Coating That Is Bonded to a Rigid Substrate
An elastic coating bonded to rigid substrate acts as a linear Winkler foundation. Then, p(x, y) and u z (x, y, 0) are related by [22]: According to the force equilibrium equation The strain energy is then where C 1 is a constant of integration. The potential energy is Using Griffith's energy equilibrium, Equation (12), G c is found to be Finally, F c can be obtained: Equation (21) is similar to Yang's solution for a cylindrical punch adhered to a thin film bonded to a rigid substrate under normal loading [23]. It is only necessary to replace a 2 with πr 2 , where r is the radius of the cylindrical punch. Equation (21) shows that F c is proportional to W ad 1/2 and E 1/2 and inversely proportional to t 1/2 . For the case of shear loading and a similar scaling law, F c ∝ W ad 1/2 or F c ∝ E 1/2 is also derived by Chaudhury [14] who demonstrated that the critical shear stress (σ s ) needed to detach the cuboid from the thin film is where G is the shear modulus of the film (G = 0.5E/(1 + v)) and l is the distance from the loading point to the interface. As described in the next section, the 3D FE method can be used to simulate both normal and shear separation.

3D FE Model
In this study, 3D FE simulation was performed using the commercial FE software ABAQUS. The ice cuboid is modeled by a square-based rigid cuboid. It is treated as a discrete rigid body using a 3D rigid quadrilateral element (R3D4). The side length of the square is 25 mm. The height of the cuboid is 15 mm. The elastic coating is modeled by an isotropous elastic film using 3D linear full-integration hexahedral element (C3D8) with material properties v, E, and t. The values of t are set to 0.25, 1, and 5 mm (t/a = 0.01, 0.04, and 0.2); v is 0.3. Since we are mainly concerned with the effect of the elastic modulus E on F c , the range of E is chosen to be 0.2 to 10,000 MPa. The size of the coating is 50 mm × 50 mm. The coating is bonded to a rigid substrate (perfect bonding, no debonding and no slip), so that a fixed boundary condition is set for the bottom surface of the coating. Displacement control normal or shear force is applied to the center of ice cuboid. The critical detachment force (F c ) can be obtained by reading the maximum force (F max ) of the displacement-reaction force curve. A bilinear cohesive zone model (CZM) [24] was used to simulate the interfacial interaction between the ice cuboid and the elastic coating. A static implicit solver was used. Quarter-model and half-model were used for the normal separation and shear separation, respectively, to reduce the amount of calculation, taking advantage of the symmetry of geometry and loading, as shown in Figure 3a,b. Therefore, F c = 4F max and 2F max for normal and shear separation, respectively. To prevent stress singularity and element distortion under shear loading, a 90 • fillet with a radius of 0.5 mm was set at the edge of the cuboid. A total of 10 nonuniform seeds are set at the fillet (bias ratio is 5) (see inset in Figure 3b).

Cohesive Zone Model (CZM)
In this study, the CZM was performed by a simple bilinear traction-separation law using surface-based cohesive behavior in ABAQUS, as illustrated in Figure 4. The bilinear traction-separation law is described by where σ and δ are current stress and current displacement, respectively. σ max is the strength of the adhesive stress, δ 0 is the damage initiation displacement, and δ f is the separation displacement at failure. D is a damage variable, K is interface stiffness. When δ ≤ δ 0 , the traction-separation law is linear and the interface stiffness K = K 0 . When δ > δ 0 , damage occurs and the interface stiffness decreases (K = (1 − D)K 0 ). When D = 1 and K = 0, the interface fails completely.
In this study, we assume that all CZM parameters are the same in three directions. In general, the F c is insensitive to the shape of the traction-separation curve when W ad is constant and δ f is much smaller than the characteristic length scales of the contacting region (a in this study) [15,16]. In this study, δ f = 0.02 mm, δf/a = 0.0008. σmax = 0.04 MPa (40 kPa) is roughly estimated from the value of ice adhesion strength (0~100 KPa) of a typical ice-phobic coating. δf is 0.02 mm, and δ 0 = 0.5δ f = 0.01 mm. The adhesion energy (W ad ) is equal to the area under the traction-separation curve (σ max × δ0, 4 × 10 −4 N/mm or mJ/mm 2 ). The initial interface stiffness (K 0 ) is equal to the slope of the front part of the traction-separation curve (K 0 = σ max /δ 0 = 4 MPa/mm). The criterion for damage initiation is described by the quadratic stress criterion where σ n , σ s , and σ t are normal, shear, and tear stress, respectively. The Macaulay brackets 〈〉 are used to signify that a purely compressive stress state does not initiate damage (i.e., 〈σ n 〉 = 0 if σ n < 0 and 〈σ n 〉 = σ n if σ n > 0). When this criterion is satisfied, damage initiates. Damage evolution and fracture criterion in mixed mode are governed by a power law with an exponent of two: where W n , W s , and W t are the work done by stresses in the normal, first, and second shear directions, respectively.

Experiment
Previous studies [10][11][12][13][14], as well as theoretical solutions and simulation results in this study, all suggest that the critical detachment force F c required to separate an ice cuboid from an elastic coating can be reduced by decreasing the elastic modulus (E) of the coating and the adhesion energy (W ad ). For ice-phobic coatings, W ad can be decreased by introducing a silicone polymer such as PDMS because of its low surface energy. For pure epoxy, it is difficult to adjust the elastic modulus. Similarly, for pure PDMS, the upper limit of its elastic modulus is no more than 10 MPa. In this study, we use epoxy/PDMS IPN gel as an elastic coating. The elastic modulus of epoxy/PDMS IPN gel can be tuned by adjusting the addition of silicone oil.

Ice Tensile Adhesion Strength
The ice tensile adhesion strength (σ ice ) test apparatus is shown in Figure 5a-d. First, the coated aluminum alloy substrate was fixed on a base with 4 screws in a refrigerator at −10 • C). A silicone mold with two open ends (cavity size: 25 × 25 × 20 mm) was put on the coating. A probe was connected to a digital force gauge, which was fixed on a hand-driven ball screw (see Figure 5c). The distance between the probe and coating surface was adjusted to~4 mm by the ball screw. Deionized water was dripped onto the silicone mold to a depth of about 1-2 mm by a dropper. The water surface tension prevents water from getting out of the gap between the mold and coating surface. After the water of 1-2 mm depth was frozen, another 13-14 mm depth of water was dripped. Two hours later, the probe was frozen into an ice block (25 × 25 ×~15 mm). Because the force is applied normally, the height of the ice block has no effect on σ ice . The normal force was applied by slowly turning the hand wheel of the ball screw (see Figure 5c). Figure 5d shows the probe with the ice block before and after the test. The maximum force (F max ) at separation can be read on the display of the digital force gauge, and σ ice can be calculated using σ ice = F max /A (A is adhesion area, 25 × 25 mm 2 ). Five equally coated substrates were tested for each coating to obtain the average and standard deviation of σ ice .

Ice Shear Adhesion Strength
Ice shear adhesion strength (τ ice ) was measured by a centrifugal adhesion test, as described in our previous work [25]. Photographs of the refrigeration centrifuge and centrifuge test setup are shown in Figure 6a,b. The ice cuboid (25 × 25 × 15 mm 3 ) was frozen onto a coated aluminum alloy substrate by using a silicone mold at −10 • C. The substrate with ice was bolted to the rotor (radius 15 cm). With an increase of rotation speed, the centrifugal force overcomes τ ice . Then, the ice detaches from the coating and impacts the protecting mesh. The rotation speed in rpm was recorded. F max can then be calculated using the centrifugal formula (F max = mω 2 r, where m is the mass of ice and r is the centrifugal radius) using SI units. Then, τ ice = F max /A. Five equally coated substrates were tested for each coating to obtain the average and standard deviation of τ ice .

Contact Angle (CA)and Elastic Modulus Measurement
The CA was measured by a contact angle measurement (YIKE-360A, Chengde Yike Instrument Factory, Chengde, China). Images of deionized water drop (~10 µL) were photographed by a camera sensor and analyzed with ImageJ plug-in Dropsnake software. The contour of the drop was determined using the piecewise polynomial fit method, and then the CA was calculated. The elastic modulus of the coatings was measured using a compression test by a universal testing machine (INSTRON 4505, USA). The test samples were cylindrical (height: 50 mm, diameter: 30 mm) and made by the casting method. The cured cylindrical samples were mounted on the universal testing machine and tested with an indenter at a speed of 0.2 mm/min. The elastic modulus data were obtained from analyzing the displacement-load curve. 3 of the same samples were prepared. Figure 7 shows the 3D FE simulation results for F c versus E for (a) normal and (b) shear separation with t/a = 0.01, 0.04, and 0.2 (a = 25 mm, t = 0.25, 1, and 5 mm). For both normal and shear separation, when E is small (<100 MPa), F c increases monotonically with E. When E reaches 100 MPa, F c is close to a 2 σ max (25 mm 2 × 0.04 MPa = 25 N), as indicated by the gray horizontal dotted line. For a film with the same E, a higher F c is required to detach the ice cuboid from the thin film (small t/a) than from the thick film (large t/a) for both normal and shear separation. For normal separation, Equation (14) (case with half-space k = 0.4265) and Equation (21) (case with film) are also plotted in Figure 7a. It can be seen that Equation (21) agrees well with the simulation results for low E and thick films (t/a = 0.04 and 0.2). For t/a = 0.01 and a large E (>~10 MPa), Equation (21) cannot match the simulation results. The F c results of Equation (21) increase rapidly with E and extend to infinity.

Critical Detachment Force F c
For shear separation, Chaudhury's result [14] (Equation (22)) is compared with simulation results. Both sides of Equation (22) can be multiplied by a 2 , giving where C is a constant of proportionality. The simulation data for t/a = 0.4 (E ranging from 0.2 to 10 MPa) were fitted using Equation (26). The fitting gives C = 0.52. Curves based on Equation (26), with three t/a values, are also plotted in Figure 7b. When E is large, Equation (26) cannot match the simulation results. Our simulation results are similar to Peng's, [15] which showed that when the film is thick-that is, large t/r (r being the radius of the rigid punch in their work)-F c increases monotonically with decreasing t/r. When the film is thin enough (t/r 1), F c will converge to Aσ max (A = πr 2 , the area of contact). For this reason, F c is determined mainly by the interface property (σ max of CZM) rather than by the material properties of the elastic film when E is large enough (as in this study) or t/r is small (as in Peng's work [15]). High E and thin film limit the deformability of the elastic film. This means that interfacial debonding will not easily occur in the form of interface crack propagation, which thus leads to higher F c . Equations (21) and (26) [14] are derived from energy equilibrium theory, and can therefore match simulation results well when E is low or t/a is not small, but cannot do so when E is high or t/a is small.
It can also be seen that F c of shear separation is smaller than F c of normal separation for the same E when E is small. This phenomenon was also found in Sun's work [16], in which the interfacial debonding of a rigid cylindrical punch from an elastic substrate using the 3D FE method was studied. This is because the shear force provides a moment (M = Fl) for the ice cuboid. This moment will increase normal stress at the left edge of the interface and make cracking easier, thus reducing F c . Obviously, increasing l will also reduce F c . When E is high enough, because the coating deformation ability is poor, separation depends on shear stress rather than normal stress at the left edge of the interface. We assume that all CZM parameters are the same in three directions, so F c of shear separation is equal to F c of normal separation (a 2 σ max , 25 N) when E is large.

Interfacial Debonding Modes
The interfacial debonding mode differs for low and high E. Force-displacement curves during detachment for normal separation, for E = 0.5 and 1000 MPa, are shown in Figure 8a,b. For E = 0.5 MPa, when the reaction force reached near F c , cracks began to appear at the edge of the interface and extended to the interior until the ice cuboid separated completely (see Figure 6b, E = 0.5 MPa). This interfacial debonding behavior can be regarded as a crack propagation mode. For E = 1000 MPa, there is no crack propagation during the whole detachment process. The damage variable distribution in the whole interface is then almost uniform (see Figure 8b, E = 1000 MPa). This interfacial debonding behavior can be regarded as a uniform separation mode. These two kinds of interfacial debonding behavior were also observed in Peng's work [15] for large and small values of t. For shear separation, the interfacial debonding behavior is similar to that for normal separation. As shown in Figure 9, when E = 0.5 MPa, cracks began to appear at the left edge of the interface and extend to the right edge until the interface completely separates. This interfacial debonding behavior can also be regarded as a crack propagation mode. When E = 1000 MPa, the evolution of the damage variable is gradual (see Figure 8b, E = 1000 MPa). Equivalently, this interfacial debonding behavior can also be regarded as a uniform separation mode.
These two interfacial debonding modes for normal and shear separation are summarized schematically in Figure 10a,b. For the crack propagation mode (Figure 10a), the crack propagates from the edge to the center (normal separation) or from the left to right edge (shear separation) due to edge stress concentration, and F c is determined mainly by the elastic modulus E of the film. For the uniform separation mode (Figure 10b), there is no obvious crack propagation, and F c is determined mainly by the interface property (σ max of CZM).   Figure 11 shows the 3D FE simulation results for F c versus v for (a) normal separation and (b) shear separation for E = 0.5, 5, and 100 MPa (t/a = 0.04). For E = 0.5 MPa in normal separation and E = 0.5 MPa and 5 MPa in shear separation, F c increases with v. This is because large v reduces the deformability and compliance of the film, resulting in higher F c . For other E values, the F c values are close to a 2 σ max (25 N). This is because the E and t/a values had already led interfacial debonding to a tendency to uniform separation mode (F c ≈ a 2 σ max ), obscuring the effect of v.  Figure 12 shows the influence of the weight ratio of silicone oil/PDMS on CA and E. CA increases monotonically with the weight ratio of silicone oil/PDMS, with values in the range of 104 • -107 • . As the fraction of silicone oil increases, E obviously decreases. When the weight ratio of silicone oil to PDMS exceeds 0.5, E is reduced by two orders of magnitude (from 118.94 MPa to less than 2.13 MPa). When the weight ratio of silicone oil to PDMS reaches 2, E falls to 0.18 MPa. This is because the silicone oil greatly decreases the crosslinking density of the matrix during the formation of epoxy/PDMS IPN gel. In general, polymer gel has a viscoelasticity behavior. In this study, we neglect the effect of viscidity (loss modulus) of the coating. In Figure 13, σ ice and τ ice are plotted vs. E on a log-log scale for both experimental results and FE simulation results.  As shown in this figure, the experimental results for σ ice and τ ice increase monotonically with E for both normal and shear separation. When E is 0.18 MPa, σ ice and τ ice are, respectively, 18.5 kPa and 8.5 kPa, which shows the good ice-phobic performance of low modulus epoxy/PDMS IPN gel coating. These values are much smaller than those of the coating with E = 119.0 MPa (σ ice = 246.6 kPa, τ ice = 150.7 kPa). In this section, the following two sets of CZM parameters were used for FE simulation results: parameter #1: W ad = 0.001 mJ/mm 2 , σ max = 0.2 MPa; parameter #2: W ad = 0.0001 mJ/mm 2 , σ max = 0.025 MPa. As shown in Figure 13a,b, the ascent phase of the FE results curve with parameter #2 matched the experimental results well when E < 1MPa for both normal and shear separation. This indicates that the interfacial debonding mode is the crack propagation mode. When E > 1 MPa, the FE results with parameter #1 only roughly match the trend of experimental results. In practice, the silicone oil not only reduces the E of epoxy/PDMS IPN gel but also, unavoidably, reduces W ad , σ max , and v. Therefore, FE simulation with a single CZM parameter cannot perfectly match the experiment.

Conclusions
In this work, we established a mechanical model to analyze the critical detachment force needed to separate an ice cuboid from an elastic coating. The critical detachment force (F c ) was analyzed by 3D FE simulation using a bilinear CZM for both normal and shear separation. The effect of the elastic modulus € of the film on F c under normal and shear loading was investigated. There are two kinds of interfacial debonding modes for different values of E, for both normal and shear separation. For a small E, the interfacial debonding mode is crack propagation, and F c increases monotonically with E. For a large E, the interfacial debonding mode is a uniform separation mode. The force (F c ) is determined by the interface property and will converge to a constant value (Aσ max ). Finally, we prepared a series of epoxy/PDMS IPN gel coatings with different elastic moduli by the addition of silicone oil. Results show that ice adhesion strengths σ ice and τ ice were 18.5 kPa and 8.5 kPa, respectively, when the elastic modulus is 0.18 MPa, showing the good ice-phobic performance of low modulus epoxy/PDMS IPN gel coating. Both simulation and experiment results can provide valuable references for a soft ice-phobic coating design.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.