Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Ice penetration by a bluff-body melting probe

Ice penetration by a bluff-body melting probe Kerstin Weinberg Department of Mechanical Engineering Universitat ¨ Siegen, 57068 Siegen, Germany Email: kerstin.weinberg@uni-siegen.de Michael Ortiz Division of Engineering and Applied Science California Institute of Technology, Pasadena, CA 91125, USA Email: ortiz@caltech.edu ABSTRACT We analyze the operation of melting probes as a Stefan problem for the liquid/solid interface sur- rounding the probe. We assume that the liquid layer is thin and, therefore, amenable to analysis by lubrication theory. The resulting Stefan problem is solvable in closed form. The solution determines the dependence of the penetration speed on the temperature differential between the probe and the surrounding ice, the size, shape and weight of the probe, the viscosity of liquid water and the thermal properties of solid ice. Keywords: melting probe, Stefan problem, liquid layer, ice penetration 1 Introduction Melting probes are hot probes that melt the surrounding ice and descend through it under the action of gravity. The probe carries its own heat source and can operate automatically in the absence of controls. The idea to use melting probes for exploring the subsurface layers of planetary ice sheets goes back, at least, to the pioneering work of Karl Philberth [1], who constructed such probes and performed field tests in Greenland. Since then, a number of programs and studies have developed the concept further, cf. [2] for a literature review. The use of melting probes for penetrating thick icy layers is particulary appealing in places where conventional drilling is infeasible, e. g., in extraterrestrial icy habitats on Mars and on the moons Europa, Enceladus and Titan, cf., e. g., [3, 4, 5, 6, 7], or in isolated terrestrial cryoenvironments, such as permafrost cryopegs and subglacial ecosystems, cf., e. g., [8]. The conventional analysis of melting probes is based on thermal analysis and balance of energy. Typically, three processes are considered in the analysis [9]: i) heating of the ice surrounding the probe to its melting temperature; ii) melting the solid ice to liquid water; and iii) conductive and radiative losses. The power required for the first two processes can be estimated simply in terms of the heat capacity of ice and the latent heat of melting. The determination of conductive losses requires the solution of a transient or steady-state heat conduction problem for the exterior domain of the probe. The speed of descent of the probe is then estimated through a balance of energy between heat supply and the latent heat required to melt the ice. However, in a more comprehensive description of the process the speed of descent results from a competition between gravity, which pushes the probe down, viscous dissipation in the liquid layer that surrounds the probe, latent heat and radiative losses. The conventional thermal analysis thus neglects potentially important parameters such as the weight of the probe and the viscosity of liquid water. In the present work, we analyze the penetration of melting probes as a Stefan problem for the liquid/solid interface surrounding the probe. We start with an analysis the temperature field of a moving probe. The liquid layer around the probe is assumed to be thin and, therefore, amenable to analysis by lubrication theory. The corresponding author arXiv:2002.02178v1 [physics.app-ph] 6 Feb 2020 resulting Stefan problem is solvable in closed form within the lubrication theory approximation. The solution determines the dependence of the penetration speed on the temperature differential between the probe and the surrounding ice, the size, shape and weight of the probe, the viscosity of liquid water and the thermal properties of solid ice. 2 Supporting experiments We begin by demonstrating the melting probe concept by means of a simple supporting experiment. The setup consists of a water ice block in a container of size 0.5 0.5 1 m, Fig. 1 (left). The walls of the container are made of transparent PMMA plates, which allows the evaluation of the ice quality, and are affixed to a steel frame. The box is filled with water and cooled down slowly to 18 C in order to obtain an ice block of good quality, i. e., without layers or large bubbles. Once frozen, the block is insolated with polystyrene panels. Fig. 1: Experimental setup: Rack with an ice block (left); lower half of the spherical penetrator (middle); detail of the heating cartridges (right). The probe is spherical of outer diameter 80 mm, is made of an aluminum (AlCuMgPb) alloy with high thermal conductivity and has a weight W = 7.24 N. The probe surface is treated by milling, drilling, scrubbing and a final surface polishing to achieve a smooth appearance. The probe consists of two identical halves, Fig. 1 (middle and right), and carries five heat cartridges. The cartridges are cylindrical (of size 10 mm  50 mm) and supply a maximum nominal power of 180 W each. They are connected to an external electric circuit that controls the power supply. In the experiment the hot probe is positioned on the top of the ice block. Driven by its own weight, it then moves through the ice, Fig. 2. In order to avoid friction with the electrical cord, the water is evacuated as is melts. We estimate the heat supply requirements by a simple energy balance Q = (c (T T ) + L) Vol (1) m ¥ ice where T = 18 C is the ice temperature, T = 0 C is the melting temperature, c = $c is the heat capacity ¥ m ice p per unit volume of ice with mass density $, L is the latent heat of melting per unit volume and Vol is the total volume of ice to melt. Reference values of these and related properties are collected in Table 1. Inserting the values into (1) and taking the dimensions of the probe and containing box into account gives an estimate of Q = 934 J. If we further assume an energy loss of /3 and a duration of the experiment of about 20 30 min, the c k c k L h ice ice water water 3 3 3 2 [J/Kcm ] [W/Km] [J/Kcm ] [W/Km] [J/cm ] [mN s/m ] 2.04 2.22 4.18 0.608 335 0.889 Table 1: Reference material constants for ice and water [10, 11]: c heat capacity per unit volume; k thermal conductivity; L latent heat of fusion per unit volume, h viscosity Fig. 2: Snapshot of the experiment as the heated probe starts to melt the ice. required power supply is estimated as Q = 900W. We thus verify that the heating cartridges located in the probe indeed have ample capacity to supply the requisite power over the duration of the experiment. ˙ ˙ Experiments were performed at two different power levels, Q = 250 W and Q = 1, 000 W. Fig. 3 shows the penetration depth as a function of time for the two cases. As may be seen from the figure, following an initial transient the penetration depth is a linear function of time, which is in turn indicative of a steady-state speed of ˙ ˙ descent. The measured speeds are: i) V = 6.9 mm/min for Q = 250 W; and ii) V = 23.9 mm/min for Q = 1, 000 W. These data correspond to the averaged values of [2]; they suggest increasing power levels are required to sustain increasing speeds of descent. Fig. 3: Measured depth vs. time for two experiments. The simple supporting experiment demonstrates the operation of melting probes and suggests a steady state of constant speed of descent of the probe as well as a dependence of the power supply on the speed. A number of other relations, such as the dependence of the power supply and speed of descent on the size and shape of the probe and on material properties such heat capacity, latent heat, thermal conductivity and viscosity remain to be elucidated. The analysis that follows is aimed at addressing these questions. 3 Thermal analysis We begin by determining the dependence of the radiated heat rate on the sinking speed. We assume that the probe is surrounded by a layer of liquid water of thickness much smaller than the dimensions of the probe. In addition, the interface between the liquid layer and the surrounding solid ice is at the melting temperature Fig. 4: Melting probe system for thermal analysis. T . Under these assumptions, the temperature field T(x,t) in the solid ice is determined by the problem ¶T c = kDT, in W(t), (2a) ¶t T = T , on G(t), (2b) T = T , r! +¥, (2c) where t denotes time, D is the Laplacian operator, W(t) is the domain occupied by the solid ice at time t, G(t) is the boundary of W(t) and the interface between the probe and solid ice, c is the heat capacity per unit volume of ice, k is the thermal conductivity of ice, T is the melting temperature of ice and T is the background ice m ¥ temperature. For simplicity, we take all material properties to be constant. We further assume a steady state wherein the probe descends through the solid ice at constant speed V. To facilitate analysis of this steady state, we introduce a reference frame moving with the probe. In this moving frame, the governing equations become ¶T cV + kDT = 0, in W, (3a) ¶x T = T , on G, (3b) T = T , r! +¥, (3c) where W is the exterior domain of the probe in an infinite body of ice, G is its boundary and x is the direction of descent. 3.1 Spherical probe The case of a spherical probe, such as used in the supporting experiments reported in Section 2, can con- veniently be solved analytically in closed form. In this case, W =fr  ag and G =fr = ag, where r is the radial distance to the origin and a is the radius of the probe. Introducing spherical coordinates x = r sinq cos j, x = r sinq sin j, x = r cosq, (4) 1 2 3 and the change of variables lx lr cosq ˜ 3 ˜ T = T + Te = T + Te , (5) ¥ ¥ eq. (3a) becomes ˜ ˜ 1 ¶ ¶T 1 ¶ ¶T 2 2 l T = r + sinq , (6) 2 2 r ¶r ¶r r sinq ¶q ¶q where cV l = (7) 2k is a normalized speed of descent. Introducing the representation T(r,q) = R(r)S(q), (8) eq. (6) further separates into the ordinary differential equations 1 d dR 2 2 2 r l r = n(n + 1), (9a) R dr dr 1 1 d dS sinq = n(n + 1), (9b) S sinq dq dq with n 0. The solutions of these equations are R (r) = A i (lr) + B k (lr), (10a) n n n n n S (q) = C P (cosq) + D Q (cosq), (10b) n n n n n where i (x) and k (x) are the modified spherical Bessel functions of the first and second kind and P (x), Q (x) n n n n are the Legendre polynomials of the first and second kind, reqectively. Boundedness at r! +¥ requires A = 0, whereas boundedness at q =p/2 requires D = 0. Under these conditions, the general solution takes the form T(r,q) = T k (lr)P (cosq). (11) n n n n=0 At the boundary r = a, we must have la cosq T(a,q) = T k (la)P (cosq) = (T T )e , (12) n n n m ¥ n=0 By the orthogonality of the Legendre polynomials, we have 2n + 1 lax T = (T T ) e P (x) dx. (13) n m ¥ n 2k (la) n 1 An appeal to Rodrigues representation, 1 d 2 n P (x) = (x 1) , (14) n n 2 n! dx further gives 2n + 1 1 d lax 2 n T = (T T ) e (x 1) dx. (15) n m ¥ n n 2k (la) 2 n! dx n 1 Integrating by parts n times, we obtain 2n + 1 (la) lax 2 n T = (T T ) e (x 1) dx, (16) n m ¥ 2k (la) 2 n! which evaluates to 2 2 1 p la 3 l a T = n + F + n, (T T ) (17) n m ¥ 0 1 2 k (la) 2 2 4 where F (a,z) is the regularized hypergeometric function F (a,z)/G(a) and G(z) is the Euler gamma function. 0 1 0 1 The corresponding heat flux into the solid ice at r = a is ¶T ¶T la cosq q(q) =k =k l cosq T e , (18) ¶r ¶r which coincides with the Dirichlet-to-Neumann map for the exterior of the spherical domain specialized to a constant surface temperature T . We expect the temperature far-field to tend asymptotically to that of a point heat source moving through the solid ice at the speed of the probe. We recall that the steady-state solution for a moving point heat source is [12] lr(1+cosq) T = T + e , (19) 4pr with l as in (7) and C > 0 quantifying the strength of the source. A straightforward calculation shows that (19) indeed satisfies (3a). We determine C by matching (19) asymptotically to the full solution (5) and (11). To leading order as r! +¥, we have lr cosq T  T + e A k (lr)P (cosq), (20) ¥ 0 0 0 with lr A = (la) (T T )(1 + coth(la)), k (lr) = , P (x) = 1, (21) 0 m ¥ 0 0 lr whereupon (20) evaluates to lr(1+cosq) T  T + (la) (T T )(1 + coth(la)) . (22) ¥ m ¥ lr Matching constants between (19) and (22), we find C = 4pla (T T )(1 + coth(la)), (23) m ¥ which fully characterizes the temperature far-field. We may now compute the total rate of heat radiated from the probe as ¶T Q = lim cVT + k 2pr sinq dq = Ck, (24) r!¥ ¶r 0 Fig. 5: Normalized radiated heat rate vs. speed of descent. which, in view of (23), evaluates to Q = 4pkla (T T )(1 + coth(la)). (25) ¥ m ¥ This identity supplies the sought relation between the speed of descent V and the radiated heat rate Q from the probe, Fig. 5. We note that the radiated heat rate increases with increasing speed of descent V, the dependence being approximately linear for large V, Q = 4pa cV(T T ). (26) ¥ m ¥ We also see from this expression that, in the large V regime the radiated heat rate scales with the square of the probe radius. For the supporting experiment, with T T = 18 K and a spherical probe of radius a = 4 cm and a velocity m ¥ of V = 7 mm/min, the radiated evaluates to Q = 88 W, or about /3 of the actual power supply. This gap suggests the operation of dissipation mechanisms other than radiative losses. We investigate one possible such mechanism, viscous dissipation in the liquid layer, next. 4 Lubrication analysis of the Stefan problem Fig. 6: Schematic of probe and liquid layer. We analyze the effect of the layer of liquid water surrounding the probe by recourse to a lubrication ap- proximation, Fig. 6. To this end, we introduce a system of cylindrical coordinates fr,q,zg aligned with the axis of descent. We further assume cylindrical symmetry and describe the boundary G of the probe by means parametric equations r(s) and z(s) in terms of the arc length s measured from the axis. The unit tangent and outward normals are 0 0 T 0 0 T m(s) = (r (s),z (s)) , n(s) = (z (s),r (s)) , (27) where the prime denotes the derivative with respect to the arc-length s. In this parametrization, we have the identity 02 02 r (s) + z (s) = 1. (28) We assume throughout steady state conditions. We note that the geometry of the liquid/solid interface, or, equivalently, the thickness of the liquid layer, is not known a priori and is an outcome of the analysis, in the spirit of Stefan problems. The main objective of the analysis is to determine the extent of viscous dissipation in the liquid layer as a function of the speed of descent. 4.1 Mass conservation Conservation of mass in a control volume dV = rh dq ds coincident with the liquid layer gives 0 0 0 (Ve ) n(s)(r(s) dq ds) = (v(s) + v (s) ds))(h(s) + h (s))(r(s) + r (s)ds) dq v(s)h(s)r(s) dq, (29) where v(s) is the mean velocity within the liquid layer and h(s) is its thickness. Simplifying, we obtain the relation Vr(s)r (s) = (r(s)h(s)v(s)). (30) ds With s measured from the axis, integration with respect to s gives r(s) = h(s)v(s), (31) whence the mean velocity of the liquid water follows as V r(s) v(s) = . (32) 2 h(s) 4.2 Weight-drag equilibrium The descent of the melting probe is driven by gravity. We assume static equilibrium between the weight W of the probe, its buoyancy $gW and the total drag force F resulting from the viscosity of the fluid layer, i. e., F = W $gW 0, (33) where $g is the mass density of water, W is the henceforth volume displaced by the probe and we assume that the probe is heavier than water. In requiring this equilibrium, we implicitly regard point-load reactions, such as would occur at a point-contact between the bottom of the probe and the ice, as unphysical. The total drag force on the probe is F = t(s)m (s)2pr(s) ds, (34) G where t(s) is the shear stress on the probe and the integral is effected over a probe meridian. Further assuming a linear velocity profile within the fluid layer, corresponding to local Couette flow, and using (27), eq. (34) becomes 2v(s) F = h z (s)2pr(s) ds. (35) h(s) Finally, inserting (32) we obtain 2 0 r (s)z (s) F = 2phV ds, (36) h (s) which scales linearly with the speed of descent for fixed fluid layer thickness. 4.3 Viscous dissipation Assuming again local Couette flow, the total dissipation evaluates to 2v(s) D = h h(s)2pr(s) ds. (37) h(s) Alternatively, inserting relation (32) we obtain r (s) D = 2phV ds, (38) h (s) which, as expected, scales with the square of the speed of descent for fixed fluid layer thickness. 4.4 Optimal fluid layer profile It remains to determine the fluid layer profile h(s). We note that h(s) defines the configuration of the fluid layer and, therefore, it may be expected to be the result of configurational equilibrium. Specifically, we assume that h(s) accommodates itself so as to minimize the total viscous dissipation (38) subject to the equilibrium constraint (34). Enforcing this constraint by means of a Lagrange multiplier m, we obtain the Lagrangian, Z Z 3 2 0 r (s) r (s)z (s) S = 2phV ds m 2phV ds . (39) 3 2 h (s) h (s) G G Stationary of S with respect to h(s) gives 3V r(s) h(s) = , (40) 2m z (s) which, inserted into (36), gives 8ph` F = m , (41) 9V with ` = z (s) ds (42) G defining an effective probe size. The equilibrium condition (33) then determines the Lagrange multiplier as 9V(W $gW) m = , (43) 8ph` whence (40) and (38) become 2ph`V r(s) h(s) = , (44) W $gW z (s) and V(W $gW) D(V) = (45) 2ph` respectively. For a spherical probe of radius a, we have r(s) = a sin(s/a), z (s) = sin(s/a), (46) with s measured from the bottom of the probe, and pG(5/4) ` = a 1.74804 a, (47) G(7/4) which, inserted into (48) and (45), fully determines the fluid layer thickness and the viscous dissipation, respec- tively. We note, in particular that the fluid layer thickness evaluates to 2ph`V h(s) = a = constant, (48) W $gW i. e., the fluid layer around a spherical probe is of uniform thickness to within the accuracy of the analysis. Fig. 7: Thickness of the liquid layer surrounding a solid spherical probe as a function of speed of descent for two different weights and under the conditions of the supporting experiment. The dependence of the liquid layer thickness on the speed of descent for the conditions of the supporting experiment are shown in Fig. 7. In all cases, the thickness of the layer is exceedingly small, which justifies a posteriori the lubrication approximation adopted in the analysis. As may be seen from the figure, the thickness increases with speed and decreases with the weight of the probe. We note that, in (48) and (45), the velocity of descent V is itself a function of the reduced weight W $gW and other factors. This dependence is determined subsequently by recourse to energy balance. 5 Speed of descent We assume that the speed of descent of the probe is the result of a balance between heat supply, radiative heat losses, viscous dissipation, latent heat and gravity in a control volume coincident with the liquid layer. From the previous analyses, this balance of energy takes the expression ˙ ˙ Q + (W $gW)V = D(V) + Q (V) + LAV (49) 0 ¥ where Q is the power supply, W is the weight of the probe, $ is the mass density of ice, g is the acceleration of gravity, W is the volume of the probe, V is the speed of descent, D is the viscous dissipation, Q (V) is the radiated heat rate, L is the latent heat rate per unit volume required to melt the solid ice and A is the cross- sectional area of the probe. Equation (49) supplies an equation that can be solved for the velocity of descent V. Equation (49) simply requires that the power supplied by the probe and the gravitational power be exactly spent on viscous dissipation, radiation losses and latent heat of melting. We note that, although the liquid water refreezes behind the probe, we assume that the latent heat of melting is not recovered by the liquid layer. This assumption is justified in view of the larger thermal conductivity of ice compared to that of water, cf. Table 1. Under these conditions, most of the recovered latent heat behind the probe leaks into the ice and does not contribute to the balance of energy of the liquid layer. Fig. 8: Speed of descent V as a function of power supply Q for the conditions of the supporting experiment. For a spherical probe, using (7), (45) and (25, eq. (49) specializes to V(W $gW) cVa 2 2 (50) Q + W $gW V = + 2pcVa (T T )(1 + coth( )) + pa LV, 0 m ¥ 2ph` 2k with ` as in (47) and W = 4pa /3, to be solved for V. For water ice, the viscosity is relatively low and the radiative losses dominate in general. For example, for the values in Table 1 and the conditions of the support- 2 5 5 ing experiment, we have: i) 4pa c(T T )  7.4 10 N, and ii) AL  16.8 10 N, which overwhelm the m ¥ remaining terms. The dependence of the resulting speed of descent on power supply is shown in Fig. 8. Also shown in the figure are two corresponding measured data points. Despite the paucity of experimental data, the agreement between theory and experiment is remarkable. 6 Concluding remarks In summary, an energy-balance argument determines the speed of descent to be the result of: i) The power supplied by the probe; ii) the buoyancy and weight of the probe; iii) viscous dissipation in the liquid layer surrounding the probe; iv) radiative heat losses, and v) latent heat of solid ice. The present analysis differs from the conventional model based on latent heat in that it additionally accounts for radiation losses, the effect of the viscous layer of liquid water surrounding the probe and the weight of the probe. In particular, the thickness of the liquid layer is unknown at the outset and its determination is part of the analysis, in the spirit of Stefan problems. Remarkably, due to the thinness of the liquid layer the Stefan problem can be solved analytically by recourse to a lubrication approximation. For water ice, the viscosity correction is small and radiative losses and latent heat dominate. The agreement between the velocity of descent predicted by the theory as a function of power supply and experimental measurements is quite remarkable and provides a modicum of validation of the theory. 7 Acknowledgment We thank the PEP-120217 student group of the University of Siegen, Germany, under the guidance of R. Notzel, ¨ for their participation in the experimental part of the study. References [1] Philberth, K., 1962. “Une methode pour mesurer les temperatures a l’interieur d’un nlandsis (a method to measure temperatures inside an ice sheet)”. Comptes Rendus des Seances de l’ Academie der Sciences, 254, p. 22. [2] Komle, N. I., Tiefenbacher, P., Weiss, P., and Bendiukova, A., 2018. “Melting probes revisited ice penetration experiments under mars surface pressure conditions”. Icarus, 308, p. 117127. [3] Chyba, C. F., and Phillips, C. B., 2001. “Possible ecosystems and the search for life on europa”. Proceedings of the National Academy of Sciences of the United States of America, 98(3), pp. 801–804. [4] Marion, G. M., Fritsen, C. H., Eicken, H., and Payne, M. C., 2003. “The search for life on europa: Limiting environmental factors, potential habitats, and earth analogues”. Astrobiology, 3(4), pp. 785–811. [5] Powell, J., Powell, J., Maise, G., and Paniagua, J., 2005. “Nemo: A mission to search for and return to earth possible life forms on europa”. Acta Astronautica, 57(2-8), pp. 579–593. [6] Hand, K. P., 2009. “Unmasking europa: The search for life on jupiter ’s ocean moon”. Nature, 457(7228), pp. 384–385. [7] Bulat, S. A., Alekhina, I. A., Marie, D., Martins, J., and Petit, J. R., 2011. “Searching for life in extreme envi- ronments relevant to jovian’s europa: Lessons from subglacial ice studies at lake vostok (east antarctica)”. Advances in Space Research, 48(4), pp. 697–701. [8] Murray, A. E., Kenig, F., Fritsen, C. H., McKay, C. P., Cawley, K. M., Edwards, R., Kuhn, E., McKnight, D. M., Ostrom, N. E., Peng, V., Ponce, A., Priscu, J. C., Samarkin, V., Townsend, A. T., Wagh, P., Young, S. A., Yung, P. T., and Doran, P. T., 2012. “Microbial life at -13 degrees c in the brine of an ice-sealed antarctic lake”. Proceedings of the National Academy of Sciences of the United States of America, 109(50), pp. 20626–20631. [9] Ulamec, S., Biele, J., Funke, O., and Engelhardt, M., 2007. “Access to glacial and subglacial environments in the solar system by meltingprobe technology”. Reviews in Environmental Science and Biotechnology, 6(1), pp. 71–94. [10] Bauer, W., Benenson, W., and Westfall, G., 2013. cliXX Physik, 2nd ed. Harri Deutsch. [11] Horvath, A. L., 1975. Physical properties of inorganic compounds: SI units. Arnold, Crane, Russak London: New York. [12] Carslaw, H. S., and Jaeger, J. C., 1986. Conduction of heat in solids, 2nd ed. Clarendon Press; Oxford University Press, Oxford Oxfordshire New York. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Physics arXiv (Cornell University)

Ice penetration by a bluff-body melting probe

Weinberg, Kerstin; Ortiz, Michael
Physics , Volume 2020 (2002) – Feb 6, 2020
Download PDF
12 pages

Loading next page...
 
/lp/arxiv-cornell-university/ice-penetration-by-a-bluff-body-melting-probe-wOrgfcyWF1

References (12)

ISSN
0021-8936
eISSN
ARCH-3341
DOI
10.1115/1.4046633
Publisher site
See Article on Publisher Site

Abstract

Kerstin Weinberg Department of Mechanical Engineering Universitat ¨ Siegen, 57068 Siegen, Germany Email: kerstin.weinberg@uni-siegen.de Michael Ortiz Division of Engineering and Applied Science California Institute of Technology, Pasadena, CA 91125, USA Email: ortiz@caltech.edu ABSTRACT We analyze the operation of melting probes as a Stefan problem for the liquid/solid interface sur- rounding the probe. We assume that the liquid layer is thin and, therefore, amenable to analysis by lubrication theory. The resulting Stefan problem is solvable in closed form. The solution determines the dependence of the penetration speed on the temperature differential between the probe and the surrounding ice, the size, shape and weight of the probe, the viscosity of liquid water and the thermal properties of solid ice. Keywords: melting probe, Stefan problem, liquid layer, ice penetration 1 Introduction Melting probes are hot probes that melt the surrounding ice and descend through it under the action of gravity. The probe carries its own heat source and can operate automatically in the absence of controls. The idea to use melting probes for exploring the subsurface layers of planetary ice sheets goes back, at least, to the pioneering work of Karl Philberth [1], who constructed such probes and performed field tests in Greenland. Since then, a number of programs and studies have developed the concept further, cf. [2] for a literature review. The use of melting probes for penetrating thick icy layers is particulary appealing in places where conventional drilling is infeasible, e. g., in extraterrestrial icy habitats on Mars and on the moons Europa, Enceladus and Titan, cf., e. g., [3, 4, 5, 6, 7], or in isolated terrestrial cryoenvironments, such as permafrost cryopegs and subglacial ecosystems, cf., e. g., [8]. The conventional analysis of melting probes is based on thermal analysis and balance of energy. Typically, three processes are considered in the analysis [9]: i) heating of the ice surrounding the probe to its melting temperature; ii) melting the solid ice to liquid water; and iii) conductive and radiative losses. The power required for the first two processes can be estimated simply in terms of the heat capacity of ice and the latent heat of melting. The determination of conductive losses requires the solution of a transient or steady-state heat conduction problem for the exterior domain of the probe. The speed of descent of the probe is then estimated through a balance of energy between heat supply and the latent heat required to melt the ice. However, in a more comprehensive description of the process the speed of descent results from a competition between gravity, which pushes the probe down, viscous dissipation in the liquid layer that surrounds the probe, latent heat and radiative losses. The conventional thermal analysis thus neglects potentially important parameters such as the weight of the probe and the viscosity of liquid water. In the present work, we analyze the penetration of melting probes as a Stefan problem for the liquid/solid interface surrounding the probe. We start with an analysis the temperature field of a moving probe. The liquid layer around the probe is assumed to be thin and, therefore, amenable to analysis by lubrication theory. The corresponding author arXiv:2002.02178v1 [physics.app-ph] 6 Feb 2020 resulting Stefan problem is solvable in closed form within the lubrication theory approximation. The solution determines the dependence of the penetration speed on the temperature differential between the probe and the surrounding ice, the size, shape and weight of the probe, the viscosity of liquid water and the thermal properties of solid ice. 2 Supporting experiments We begin by demonstrating the melting probe concept by means of a simple supporting experiment. The setup consists of a water ice block in a container of size 0.5 0.5 1 m, Fig. 1 (left). The walls of the container are made of transparent PMMA plates, which allows the evaluation of the ice quality, and are affixed to a steel frame. The box is filled with water and cooled down slowly to 18 C in order to obtain an ice block of good quality, i. e., without layers or large bubbles. Once frozen, the block is insolated with polystyrene panels. Fig. 1: Experimental setup: Rack with an ice block (left); lower half of the spherical penetrator (middle); detail of the heating cartridges (right). The probe is spherical of outer diameter 80 mm, is made of an aluminum (AlCuMgPb) alloy with high thermal conductivity and has a weight W = 7.24 N. The probe surface is treated by milling, drilling, scrubbing and a final surface polishing to achieve a smooth appearance. The probe consists of two identical halves, Fig. 1 (middle and right), and carries five heat cartridges. The cartridges are cylindrical (of size 10 mm  50 mm) and supply a maximum nominal power of 180 W each. They are connected to an external electric circuit that controls the power supply. In the experiment the hot probe is positioned on the top of the ice block. Driven by its own weight, it then moves through the ice, Fig. 2. In order to avoid friction with the electrical cord, the water is evacuated as is melts. We estimate the heat supply requirements by a simple energy balance Q = (c (T T ) + L) Vol (1) m ¥ ice where T = 18 C is the ice temperature, T = 0 C is the melting temperature, c = $c is the heat capacity ¥ m ice p per unit volume of ice with mass density $, L is the latent heat of melting per unit volume and Vol is the total volume of ice to melt. Reference values of these and related properties are collected in Table 1. Inserting the values into (1) and taking the dimensions of the probe and containing box into account gives an estimate of Q = 934 J. If we further assume an energy loss of /3 and a duration of the experiment of about 20 30 min, the c k c k L h ice ice water water 3 3 3 2 [J/Kcm ] [W/Km] [J/Kcm ] [W/Km] [J/cm ] [mN s/m ] 2.04 2.22 4.18 0.608 335 0.889 Table 1: Reference material constants for ice and water [10, 11]: c heat capacity per unit volume; k thermal conductivity; L latent heat of fusion per unit volume, h viscosity Fig. 2: Snapshot of the experiment as the heated probe starts to melt the ice. required power supply is estimated as Q = 900W. We thus verify that the heating cartridges located in the probe indeed have ample capacity to supply the requisite power over the duration of the experiment. ˙ ˙ Experiments were performed at two different power levels, Q = 250 W and Q = 1, 000 W. Fig. 3 shows the penetration depth as a function of time for the two cases. As may be seen from the figure, following an initial transient the penetration depth is a linear function of time, which is in turn indicative of a steady-state speed of ˙ ˙ descent. The measured speeds are: i) V = 6.9 mm/min for Q = 250 W; and ii) V = 23.9 mm/min for Q = 1, 000 W. These data correspond to the averaged values of [2]; they suggest increasing power levels are required to sustain increasing speeds of descent. Fig. 3: Measured depth vs. time for two experiments. The simple supporting experiment demonstrates the operation of melting probes and suggests a steady state of constant speed of descent of the probe as well as a dependence of the power supply on the speed. A number of other relations, such as the dependence of the power supply and speed of descent on the size and shape of the probe and on material properties such heat capacity, latent heat, thermal conductivity and viscosity remain to be elucidated. The analysis that follows is aimed at addressing these questions. 3 Thermal analysis We begin by determining the dependence of the radiated heat rate on the sinking speed. We assume that the probe is surrounded by a layer of liquid water of thickness much smaller than the dimensions of the probe. In addition, the interface between the liquid layer and the surrounding solid ice is at the melting temperature Fig. 4: Melting probe system for thermal analysis. T . Under these assumptions, the temperature field T(x,t) in the solid ice is determined by the problem ¶T c = kDT, in W(t), (2a) ¶t T = T , on G(t), (2b) T = T , r! +¥, (2c) where t denotes time, D is the Laplacian operator, W(t) is the domain occupied by the solid ice at time t, G(t) is the boundary of W(t) and the interface between the probe and solid ice, c is the heat capacity per unit volume of ice, k is the thermal conductivity of ice, T is the melting temperature of ice and T is the background ice m ¥ temperature. For simplicity, we take all material properties to be constant. We further assume a steady state wherein the probe descends through the solid ice at constant speed V. To facilitate analysis of this steady state, we introduce a reference frame moving with the probe. In this moving frame, the governing equations become ¶T cV + kDT = 0, in W, (3a) ¶x T = T , on G, (3b) T = T , r! +¥, (3c) where W is the exterior domain of the probe in an infinite body of ice, G is its boundary and x is the direction of descent. 3.1 Spherical probe The case of a spherical probe, such as used in the supporting experiments reported in Section 2, can con- veniently be solved analytically in closed form. In this case, W =fr  ag and G =fr = ag, where r is the radial distance to the origin and a is the radius of the probe. Introducing spherical coordinates x = r sinq cos j, x = r sinq sin j, x = r cosq, (4) 1 2 3 and the change of variables lx lr cosq ˜ 3 ˜ T = T + Te = T + Te , (5) ¥ ¥ eq. (3a) becomes ˜ ˜ 1 ¶ ¶T 1 ¶ ¶T 2 2 l T = r + sinq , (6) 2 2 r ¶r ¶r r sinq ¶q ¶q where cV l = (7) 2k is a normalized speed of descent. Introducing the representation T(r,q) = R(r)S(q), (8) eq. (6) further separates into the ordinary differential equations 1 d dR 2 2 2 r l r = n(n + 1), (9a) R dr dr 1 1 d dS sinq = n(n + 1), (9b) S sinq dq dq with n 0. The solutions of these equations are R (r) = A i (lr) + B k (lr), (10a) n n n n n S (q) = C P (cosq) + D Q (cosq), (10b) n n n n n where i (x) and k (x) are the modified spherical Bessel functions of the first and second kind and P (x), Q (x) n n n n are the Legendre polynomials of the first and second kind, reqectively. Boundedness at r! +¥ requires A = 0, whereas boundedness at q =p/2 requires D = 0. Under these conditions, the general solution takes the form T(r,q) = T k (lr)P (cosq). (11) n n n n=0 At the boundary r = a, we must have la cosq T(a,q) = T k (la)P (cosq) = (T T )e , (12) n n n m ¥ n=0 By the orthogonality of the Legendre polynomials, we have 2n + 1 lax T = (T T ) e P (x) dx. (13) n m ¥ n 2k (la) n 1 An appeal to Rodrigues representation, 1 d 2 n P (x) = (x 1) , (14) n n 2 n! dx further gives 2n + 1 1 d lax 2 n T = (T T ) e (x 1) dx. (15) n m ¥ n n 2k (la) 2 n! dx n 1 Integrating by parts n times, we obtain 2n + 1 (la) lax 2 n T = (T T ) e (x 1) dx, (16) n m ¥ 2k (la) 2 n! which evaluates to 2 2 1 p la 3 l a T = n + F + n, (T T ) (17) n m ¥ 0 1 2 k (la) 2 2 4 where F (a,z) is the regularized hypergeometric function F (a,z)/G(a) and G(z) is the Euler gamma function. 0 1 0 1 The corresponding heat flux into the solid ice at r = a is ¶T ¶T la cosq q(q) =k =k l cosq T e , (18) ¶r ¶r which coincides with the Dirichlet-to-Neumann map for the exterior of the spherical domain specialized to a constant surface temperature T . We expect the temperature far-field to tend asymptotically to that of a point heat source moving through the solid ice at the speed of the probe. We recall that the steady-state solution for a moving point heat source is [12] lr(1+cosq) T = T + e , (19) 4pr with l as in (7) and C > 0 quantifying the strength of the source. A straightforward calculation shows that (19) indeed satisfies (3a). We determine C by matching (19) asymptotically to the full solution (5) and (11). To leading order as r! +¥, we have lr cosq T  T + e A k (lr)P (cosq), (20) ¥ 0 0 0 with lr A = (la) (T T )(1 + coth(la)), k (lr) = , P (x) = 1, (21) 0 m ¥ 0 0 lr whereupon (20) evaluates to lr(1+cosq) T  T + (la) (T T )(1 + coth(la)) . (22) ¥ m ¥ lr Matching constants between (19) and (22), we find C = 4pla (T T )(1 + coth(la)), (23) m ¥ which fully characterizes the temperature far-field. We may now compute the total rate of heat radiated from the probe as ¶T Q = lim cVT + k 2pr sinq dq = Ck, (24) r!¥ ¶r 0 Fig. 5: Normalized radiated heat rate vs. speed of descent. which, in view of (23), evaluates to Q = 4pkla (T T )(1 + coth(la)). (25) ¥ m ¥ This identity supplies the sought relation between the speed of descent V and the radiated heat rate Q from the probe, Fig. 5. We note that the radiated heat rate increases with increasing speed of descent V, the dependence being approximately linear for large V, Q = 4pa cV(T T ). (26) ¥ m ¥ We also see from this expression that, in the large V regime the radiated heat rate scales with the square of the probe radius. For the supporting experiment, with T T = 18 K and a spherical probe of radius a = 4 cm and a velocity m ¥ of V = 7 mm/min, the radiated evaluates to Q = 88 W, or about /3 of the actual power supply. This gap suggests the operation of dissipation mechanisms other than radiative losses. We investigate one possible such mechanism, viscous dissipation in the liquid layer, next. 4 Lubrication analysis of the Stefan problem Fig. 6: Schematic of probe and liquid layer. We analyze the effect of the layer of liquid water surrounding the probe by recourse to a lubrication ap- proximation, Fig. 6. To this end, we introduce a system of cylindrical coordinates fr,q,zg aligned with the axis of descent. We further assume cylindrical symmetry and describe the boundary G of the probe by means parametric equations r(s) and z(s) in terms of the arc length s measured from the axis. The unit tangent and outward normals are 0 0 T 0 0 T m(s) = (r (s),z (s)) , n(s) = (z (s),r (s)) , (27) where the prime denotes the derivative with respect to the arc-length s. In this parametrization, we have the identity 02 02 r (s) + z (s) = 1. (28) We assume throughout steady state conditions. We note that the geometry of the liquid/solid interface, or, equivalently, the thickness of the liquid layer, is not known a priori and is an outcome of the analysis, in the spirit of Stefan problems. The main objective of the analysis is to determine the extent of viscous dissipation in the liquid layer as a function of the speed of descent. 4.1 Mass conservation Conservation of mass in a control volume dV = rh dq ds coincident with the liquid layer gives 0 0 0 (Ve ) n(s)(r(s) dq ds) = (v(s) + v (s) ds))(h(s) + h (s))(r(s) + r (s)ds) dq v(s)h(s)r(s) dq, (29) where v(s) is the mean velocity within the liquid layer and h(s) is its thickness. Simplifying, we obtain the relation Vr(s)r (s) = (r(s)h(s)v(s)). (30) ds With s measured from the axis, integration with respect to s gives r(s) = h(s)v(s), (31) whence the mean velocity of the liquid water follows as V r(s) v(s) = . (32) 2 h(s) 4.2 Weight-drag equilibrium The descent of the melting probe is driven by gravity. We assume static equilibrium between the weight W of the probe, its buoyancy $gW and the total drag force F resulting from the viscosity of the fluid layer, i. e., F = W $gW 0, (33) where $g is the mass density of water, W is the henceforth volume displaced by the probe and we assume that the probe is heavier than water. In requiring this equilibrium, we implicitly regard point-load reactions, such as would occur at a point-contact between the bottom of the probe and the ice, as unphysical. The total drag force on the probe is F = t(s)m (s)2pr(s) ds, (34) G where t(s) is the shear stress on the probe and the integral is effected over a probe meridian. Further assuming a linear velocity profile within the fluid layer, corresponding to local Couette flow, and using (27), eq. (34) becomes 2v(s) F = h z (s)2pr(s) ds. (35) h(s) Finally, inserting (32) we obtain 2 0 r (s)z (s) F = 2phV ds, (36) h (s) which scales linearly with the speed of descent for fixed fluid layer thickness. 4.3 Viscous dissipation Assuming again local Couette flow, the total dissipation evaluates to 2v(s) D = h h(s)2pr(s) ds. (37) h(s) Alternatively, inserting relation (32) we obtain r (s) D = 2phV ds, (38) h (s) which, as expected, scales with the square of the speed of descent for fixed fluid layer thickness. 4.4 Optimal fluid layer profile It remains to determine the fluid layer profile h(s). We note that h(s) defines the configuration of the fluid layer and, therefore, it may be expected to be the result of configurational equilibrium. Specifically, we assume that h(s) accommodates itself so as to minimize the total viscous dissipation (38) subject to the equilibrium constraint (34). Enforcing this constraint by means of a Lagrange multiplier m, we obtain the Lagrangian, Z Z 3 2 0 r (s) r (s)z (s) S = 2phV ds m 2phV ds . (39) 3 2 h (s) h (s) G G Stationary of S with respect to h(s) gives 3V r(s) h(s) = , (40) 2m z (s) which, inserted into (36), gives 8ph` F = m , (41) 9V with ` = z (s) ds (42) G defining an effective probe size. The equilibrium condition (33) then determines the Lagrange multiplier as 9V(W $gW) m = , (43) 8ph` whence (40) and (38) become 2ph`V r(s) h(s) = , (44) W $gW z (s) and V(W $gW) D(V) = (45) 2ph` respectively. For a spherical probe of radius a, we have r(s) = a sin(s/a), z (s) = sin(s/a), (46) with s measured from the bottom of the probe, and pG(5/4) ` = a 1.74804 a, (47) G(7/4) which, inserted into (48) and (45), fully determines the fluid layer thickness and the viscous dissipation, respec- tively. We note, in particular that the fluid layer thickness evaluates to 2ph`V h(s) = a = constant, (48) W $gW i. e., the fluid layer around a spherical probe is of uniform thickness to within the accuracy of the analysis. Fig. 7: Thickness of the liquid layer surrounding a solid spherical probe as a function of speed of descent for two different weights and under the conditions of the supporting experiment. The dependence of the liquid layer thickness on the speed of descent for the conditions of the supporting experiment are shown in Fig. 7. In all cases, the thickness of the layer is exceedingly small, which justifies a posteriori the lubrication approximation adopted in the analysis. As may be seen from the figure, the thickness increases with speed and decreases with the weight of the probe. We note that, in (48) and (45), the velocity of descent V is itself a function of the reduced weight W $gW and other factors. This dependence is determined subsequently by recourse to energy balance. 5 Speed of descent We assume that the speed of descent of the probe is the result of a balance between heat supply, radiative heat losses, viscous dissipation, latent heat and gravity in a control volume coincident with the liquid layer. From the previous analyses, this balance of energy takes the expression ˙ ˙ Q + (W $gW)V = D(V) + Q (V) + LAV (49) 0 ¥ where Q is the power supply, W is the weight of the probe, $ is the mass density of ice, g is the acceleration of gravity, W is the volume of the probe, V is the speed of descent, D is the viscous dissipation, Q (V) is the radiated heat rate, L is the latent heat rate per unit volume required to melt the solid ice and A is the cross- sectional area of the probe. Equation (49) supplies an equation that can be solved for the velocity of descent V. Equation (49) simply requires that the power supplied by the probe and the gravitational power be exactly spent on viscous dissipation, radiation losses and latent heat of melting. We note that, although the liquid water refreezes behind the probe, we assume that the latent heat of melting is not recovered by the liquid layer. This assumption is justified in view of the larger thermal conductivity of ice compared to that of water, cf. Table 1. Under these conditions, most of the recovered latent heat behind the probe leaks into the ice and does not contribute to the balance of energy of the liquid layer. Fig. 8: Speed of descent V as a function of power supply Q for the conditions of the supporting experiment. For a spherical probe, using (7), (45) and (25, eq. (49) specializes to V(W $gW) cVa 2 2 (50) Q + W $gW V = + 2pcVa (T T )(1 + coth( )) + pa LV, 0 m ¥ 2ph` 2k with ` as in (47) and W = 4pa /3, to be solved for V. For water ice, the viscosity is relatively low and the radiative losses dominate in general. For example, for the values in Table 1 and the conditions of the support- 2 5 5 ing experiment, we have: i) 4pa c(T T )  7.4 10 N, and ii) AL  16.8 10 N, which overwhelm the m ¥ remaining terms. The dependence of the resulting speed of descent on power supply is shown in Fig. 8. Also shown in the figure are two corresponding measured data points. Despite the paucity of experimental data, the agreement between theory and experiment is remarkable. 6 Concluding remarks In summary, an energy-balance argument determines the speed of descent to be the result of: i) The power supplied by the probe; ii) the buoyancy and weight of the probe; iii) viscous dissipation in the liquid layer surrounding the probe; iv) radiative heat losses, and v) latent heat of solid ice. The present analysis differs from the conventional model based on latent heat in that it additionally accounts for radiation losses, the effect of the viscous layer of liquid water surrounding the probe and the weight of the probe. In particular, the thickness of the liquid layer is unknown at the outset and its determination is part of the analysis, in the spirit of Stefan problems. Remarkably, due to the thinness of the liquid layer the Stefan problem can be solved analytically by recourse to a lubrication approximation. For water ice, the viscosity correction is small and radiative losses and latent heat dominate. The agreement between the velocity of descent predicted by the theory as a function of power supply and experimental measurements is quite remarkable and provides a modicum of validation of the theory. 7 Acknowledgment We thank the PEP-120217 student group of the University of Siegen, Germany, under the guidance of R. Notzel, ¨ for their participation in the experimental part of the study. References [1] Philberth, K., 1962. “Une methode pour mesurer les temperatures a l’interieur d’un nlandsis (a method to measure temperatures inside an ice sheet)”. Comptes Rendus des Seances de l’ Academie der Sciences, 254, p. 22. [2] Komle, N. I., Tiefenbacher, P., Weiss, P., and Bendiukova, A., 2018. “Melting probes revisited ice penetration experiments under mars surface pressure conditions”. Icarus, 308, p. 117127. [3] Chyba, C. F., and Phillips, C. B., 2001. “Possible ecosystems and the search for life on europa”. Proceedings of the National Academy of Sciences of the United States of America, 98(3), pp. 801–804. [4] Marion, G. M., Fritsen, C. H., Eicken, H., and Payne, M. C., 2003. “The search for life on europa: Limiting environmental factors, potential habitats, and earth analogues”. Astrobiology, 3(4), pp. 785–811. [5] Powell, J., Powell, J., Maise, G., and Paniagua, J., 2005. “Nemo: A mission to search for and return to earth possible life forms on europa”. Acta Astronautica, 57(2-8), pp. 579–593. [6] Hand, K. P., 2009. “Unmasking europa: The search for life on jupiter ’s ocean moon”. Nature, 457(7228), pp. 384–385. [7] Bulat, S. A., Alekhina, I. A., Marie, D., Martins, J., and Petit, J. R., 2011. “Searching for life in extreme envi- ronments relevant to jovian’s europa: Lessons from subglacial ice studies at lake vostok (east antarctica)”. Advances in Space Research, 48(4), pp. 697–701. [8] Murray, A. E., Kenig, F., Fritsen, C. H., McKay, C. P., Cawley, K. M., Edwards, R., Kuhn, E., McKnight, D. M., Ostrom, N. E., Peng, V., Ponce, A., Priscu, J. C., Samarkin, V., Townsend, A. T., Wagh, P., Young, S. A., Yung, P. T., and Doran, P. T., 2012. “Microbial life at -13 degrees c in the brine of an ice-sealed antarctic lake”. Proceedings of the National Academy of Sciences of the United States of America, 109(50), pp. 20626–20631. [9] Ulamec, S., Biele, J., Funke, O., and Engelhardt, M., 2007. “Access to glacial and subglacial environments in the solar system by meltingprobe technology”. Reviews in Environmental Science and Biotechnology, 6(1), pp. 71–94. [10] Bauer, W., Benenson, W., and Westfall, G., 2013. cliXX Physik, 2nd ed. Harri Deutsch. [11] Horvath, A. L., 1975. Physical properties of inorganic compounds: SI units. Arnold, Crane, Russak London: New York. [12] Carslaw, H. S., and Jaeger, J. C., 1986. Conduction of heat in solids, 2nd ed. Clarendon Press; Oxford University Press, Oxford Oxfordshire New York.

Journal

Physics – arXiv (Cornell University)

Published: Feb 6, 2020

There are no references for this article.