One dimensional effervescence modeling of an extraterrestrial submarine in the Saturn Titan Seas

Planetary and Space Science xxx (xxxx) xxx

Contents lists available at ScienceDirect

Planetary and Space Science journal homepage: www.elsevier.com/locate/pss

One dimensional effervescence modeling of an extraterrestrial submarine in the Saturn Titan Seas Jason Hartwig a, *, Peter Meyerhofer b, R. Balasubramaniam b, Mariano Mercado c, Ralph Lorenz d, Justin Walsh e, Steve Oleson a a

NASA Glenn Research Center, Cleveland, OH, 44135, USA Case Western Reserve University, Mechanical and Aerospace Engineering, Cleveland, OH, 44106, USA University of Texas El Paso, Mechanical Engineering, El Paso, TX, 79968, USA d Johns Hopkins Applied Physics Laboratory, Laurel, MD, 20723, USA e Penn State Applied Research Lab, State College, PA, 16801, USA b c

A R T I C L E I N F O

A B S T R A C T

Keywords: Titan Submarine Effervescence Nitrogen/methane/ethane Bubble incipience Bubble growth

An extraterrestrial submarine has been proposed to explore Saturn's moon Titan. A major design concern is the effect of effervescence on submarine operation. Because nitrogen gas has a relatively high solubility in the liquid ethane and methane seas, waste heat generated by the submarine power system may cause dissolved nitrogen gas to come out of solution. In a non-moving case, bubbles that form may interfere with science measurements. In a moving case, bubbles that form along the body may coalesce at the aft end of the craft and cause cavitation in the propellers. This paper presents a concise effervescence model as a function of the waste heat, location of submarine within the seas (liquid temperature and pressure), and contact angle between bubbles and submarine skin. Based on results, effervescence is a strong function of heat flux into the liquid as well as contact angle, possibly requiring limits on the power system or a redesign of the thermal management system. The current submarine operating waste heat flux is 370 W/m2, which when compared to recent bubble incipience data suggests a small margin in safety factor, however.

1.1. Titan and the submarine

1. Introduction The purpose of this paper is to present practical models needed to determine the effect of effervescence on a submersible body within the seas of Saturn's moon Titan, in both a non-moving and moving case, as a function of relevant parameters, such as heat flux into the liquid, mole fraction of ethane and methane, and temperature and pressure of the sea. The outline of the paper is as follows: First, a general introduction is given regarding Titan, its seas, missions to Titan, and the potential problem of effervescence. Next, fundamental models for solubility, bubble incipience, growth, and area and volume coverage are presented. Then, numerical results are presented, and effervescence is quantified as a function of the key driving parameters. Finally, implications for future extraterrestrial submersible designs are discussed. Details of the skin temperature calculations are presented in Appendix A.

A submarine has been proposed to explore the cryogenic liquid hydrocarbon seas of Saturn's moon Titan (Oleson et al., 2014) (Hartwig et al., 2016). The design evolved in a 1-year Phase I NASA Innovative Advanced Concepts (NIAC) study is depicted in Fig. 1. The concept developed was for a one year mission on Titan of a unique submersible, fully autonomous design of a submarine capable of operating within the Titan seas to carry out detailed scientific investigations. The purpose of the mission is to study the history and evolution of hydrocarbons in the solar system, to determine if life is possible on Titan, and to provide a pathfinder for later design of submersibles in the seas hidden beneath the ice crust of other outer planetary moons (e.g. Europa, Ganymede, and Enceladus). Titan is unique within the Solar System because: 1. It is the only known body other than Earth that has stable, accessible seas. While it is anticipated that there are internal 'oceans' elsewhere in the Solar System such as Europa, none are readily accessible without first drilling through many kilometers of ice crust. The seas of

* Corresponding author. E-mail address: [email protected] (J. Hartwig). https://doi.org/10.1016/j.pss.2019.04.002 Received 29 September 2018; Received in revised form 28 March 2019; Accepted 2 April 2019 Available online xxxx 0032-0633/© 2019 Published by Elsevier Ltd.

Please cite this article as: Hartwig, J. et al., One dimensional effervescence modeling of an extraterrestrial submarine in the Saturn Titan Seas, Planetary and Space Science, https://doi.org/10.1016/j.pss.2019.04.002

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

Nomenclature and units

Greek Bubble growth constant β dimensionless Surface tension γ N/m Similarity variable η dimensionless Contact angle θ degrees or radians Density ρ kg/m3 Density ratio ρ* dimensionless Solidity ratio σ . dimensionless Thrust τ N Kinematic viscosity ξ m2/s

Area A m2 Normalized mass concentration C dimensionless Drag coefficient CD dimensionless Diffusion coefficient D m2/s Force F Newtons Gravitational acceleration g m/s2 Height H meters Length L meters Nucleation site density Nn 1/m2 Pressure P MegaPascals Excess above vapor pressure P* MegaPascals Heat transfer Q_ Watts Bubble radius R meters Reynolds number Re dimensionless Supersaturation S dimensionless Temperature T Kelvin Time t seconds Width W meters Volume V m3 Volume flow rate V_ m3/s Velocity v m/s Propeller velocity change Δv m/s Lengthwise coordinate x meters Mole fraction xðÞ dimensionless Vertical coordinate y meters

Subscripts b c cr dr front g l nc p skin sub wh

buoyancy force critical bubble radius maximum growth time in quiescent case drag force submarine front (projected) nitrogen vapor in bubble Titan seas induced velocity due to natural convection low-pressure region behind the propellers submarine skin submarine dimensions or speed submarine waste heat

Superscripts ‘ representative lengths in area/volume calculations

due to the thermal management system, which is based on an energy balance between the heat generated by the radioisotope power system, and the heat loss to the surroundings through the submarine exterior (Hartwig et al., 2016). Waste heat (Lorenz, 2016) is distributed from the power system at the aft end of the sub to the forward end; and along with the appropriately sized insulation thickness, the science equipment is maintained at ambient. The rest of the heat is rejected into the surrounding sea. Fig. 2 illustrates the anticipated temperature gradients between the submarine and the Titan seas. As shown, there is a significant temperature difference between power system and sea. This temperature difference is not likely to cause local boiling of the liquid, due to the fact that the seas are so close to the freezing point, but it may be enough to cause nitrogen gas that is dissolved in the liquid to come out of solution. On Titan, the submarine may see pressures in excess of 1.03 MPa (150 psia). For terrestrial submarines, the solubility of air in water at that pressure is negligible (<0.1% mol from (“Air Solubility in Water, 2016)). On Titan however, the atmospheric pressure is 1.5 times higher. Solubility of gaseous nitrogen in liquid ethane or methane is thus expected to be a concern for all submersible designs on Titan. Specifically as it relates to the submarine operation, bubbles that come out of solution due to

Fig. 1. Illustration of submarine during a dive on Titan (Oleson et al., 2014).

Titan are primarily composed of liquid methane, liquid ethane, and nitrogen at a surface temperature between 90 and 96K. Gravity on Titan is roughly 1/6 that of Earth. 2. It is the only known moon with a significant atmosphere. The same hydrological cycles that occur on Earth also occur on Titan, except the working fluid is methane. The atmospheric pressure near the seas is 1.5 times that of atmospheric pressure on Earth, with the primary constituents being gaseous nitrogen (GN2) (95%) and gaseous methane (5%). The cryogenic conditions of Titan represent a unique and unprecedented design challenge for submersibles capable of exploring the seas.

1.2. Problem statement The submarine will be subject to cryogenic temperatures, therefore all external equipment needs to endure, operate, and cycle at temperatures below 96K. Meanwhile, all internal systems are maintained at 290K

Fig. 2. Internal temperatures and parts of the submarine, with insulation lines. Blocks are shown for subsystems (Purple ¼ science, blue ¼ communication and data handling, red ¼ power). 2

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

excessive waste heat will cause problems in two ways:

where A, B and C are constants. The fitted values of these constants are shown in Table 2 with units and sources. Those results only apply for the pure fluids methane and ethane, hence for the binary solubility cases. A decent approximation for ternary solubility is simple linear interpolation by methane mole fraction. Above the critical temperature of methane (190 K), however, such interpolation has no methane endpoint and approximation is more difficult. The vapor pressure of nitrogen is also included as a reference; it is not needed to calculate solubility directly, but it provides an upper bound for the absolute pressure within which the model is fitted (below the critical point, 126 K).

1. In a quiescent or hovering case, bubbles that form may interfere with instrumentation such as dielectric constant measurements, turbidity, acoustic transmission (Malaska et al., 2017), sonar (Arvelo and Lorenz, 2013), depth, and any visualization. This would jeopardize critical mission objectives. 2. In a moving case, bubbles that form at the forward end and along the craft may coalesce at the aft end and cause cavitation in the propellers, potentially hindering control and navigation of the sub. In order to faithfully quantify effervescence, sub-models are needed for skin temperature, solubility, bubble nucleation, bubble growth, and bubble development as a function of location with the Titan sea system.

2.2. General bubble nucleation Given that solubility limits can determine how much gas is available in the liquid, conditions for when a bubble forms, or nucleates, are next required. Jones et al. (1999) presented a thorough review of bubble nucleation in liquids supersaturated with a gas. Four types of nucleation are classified:

2. Modeling 2.1. Solubility Supersaturation is defined as the ratio of the difference in the dissolved solute concentration in the liquid and the bubble surface, to that at the bubble surface: S¼

xb  xi xi

1. Classical homogenous nucleation where bubbles are formed in the bulk liquid. 2. Classical heterogeneous nucleation where bubbles form on the solid surface in contact with the liquid. 3. Nucleation at pre-existing gas cavities on solid surfaces with a nucleation energy barrier. 4. Nucleation at pre-existing gas cavities on solid surfaces without an energy barrier.

(1)

where x is the mole fraction of the solute (dissolved gas) in the liquid, and subscripts b and i denote bulk liquid (Titan sea) and interface (at the skin surface of the submarine) values, respectively. xb and xi are functions of the temperature of the bulk sea and the submarine temperature, respectively, as well as the pressure. Calculations for skin temperature as a function of waste heat are in Appendix A. First, solubility models are needed to determine the potential amount of dissolved nitrogen gas in the bulk Titan seas xb and near the surface of the submarine xi , respectively. Recently, a set of universal solubility correlations (Hartwig et al., 2017) was developed based on 2300 data points for predicting the solubility of nitrogen gas in binary and ternary mixtures of liquid ethane and methane. For ternary ethane/methane/nitrogen mixtures representative of Titan seas, the mole fraction solubility of nitrogen is predicted through the expression from:  a6
 xN2 xC1 ¼ a0 exp a1 T þ a2 T 2 þ a3 P* þ a4 TP* þ a5 T 2 P* P* xC1 þ xC2

For bubbles to grow, the initial nucleus size must be larger than some critical radius. Smaller bubbles will dissolve back into solution. Wilt et al. (Wilt, 1986) calculated nucleation rates for Type 1 and 2 given above for water/carbon dioxide solutions and found the nucleation rates to be quite small and unlikely to lead to experimentally observable bubbles unless a) the supersaturation ratio is very high (1100–1700) for homogeneous nucleation, or b) the liquid is very non-wetting for a supersaturation ratio of 5. The supersaturation ratio for nitrogen gas in the liquid methane/ethane sea is less than 1 for all reasonable submarine conditions. Further, the liquid methane/ethane mixture is assumed to be quite wetting (near-zero contact angle), as is the case for all cryogenic liquids in contact with metallic surfaces (Hartwig and Mann, 2014). Under these conditions, Type 1 or Type 2 nucleation is highly unlikely to occur for the Titan submarine; the most probable type of nucleation is at pre-existing gas cavities on the surface without an energy barrier (Type 4). The active nucleation site density on surfaces in boiling systems is well researched, see for example (Hibiki and Tshii, 2003). Very few quantitative studies appear to exist for gas evolution on surfaces in contact with supersaturated solutions (see for example (Qi and Klausner, 2006) (Liger-Belair et al., 1999)). For analysis in the current work, the model of (Hibiki and Tshii, 2003) will be employed, modified for

(2)

where xC1 and xC2 are the mole fraction of LCH4 and LC2H6, respectively. The values a0 to a6 are coefficients listed in Table 1. The pressure P* ¼ ðP  P0 Þ is the excess above vapor pressure, in MPa. The vapor pressure P0 used in Equation (2) is that of the binary methane/ethane mixture. This work implements the value given by the software REFPROP, which is the pressure corresponding to thermodynamic quality at a given temperature. A convenient alternative, without special software, is to use a correlation of the Antoine equation (Methane, 2016) (Ethane, 2016): log10 P0 ¼ A 

B T þC

Table 2 Constant values for the Antoine equation for vapor pressure, evaluated in units of K and MPa. These values are fitted by NIST for nitrogen (National Institute of Standards and Technology, 2016), methane (Methane, 2016) and ethane (Ethane, 2016) from the sources cited, and adapted from pressure units of bar.

(3)

Table 1 Fitted coefficients for nitrogen solubility in methane/ethane solvent, 90–100 K.

a0 ðMPa1 Þ a1 ðK 1 Þ a2 ðK 2 Þ a3 ðMPa1 Þ a4 ðMPa1 K 1 Þ a5 ðMPa1 K 2 Þ a6

Methane/ethane sea

Ethane sea

46.5 0.0439 0.00001 3.07 0.000044 0.00001 0.933

3.509 -.0307 0 0.556 0.00369 0 0

3

Fluid

T Range (K)

A

B

C

Nitrogen (Edejer and Thodos, 1967) Methane (Prydz and Goodwin, 1972) Ethane (Carruth and Kobayashi, 1973) Ethane (Loomis and Walters, 1926)

63 to 126

2.7362

264.651

6.788

90 to 190

2.9895

443.028

0.49

91 to 140

3.507

791.3

6.422

140 to 200

2.9385

659.739

16.719

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

studied by (Cable and Frade, 1987). Thus, Scriven's theory is used to assess the rate of growth of nitrogen bubbles once they are nucleated on the surface of the submarine. Two assumptions are needed:

effervescence in place of boiling. 2.3. Bubble incipience/nucleation site density Hibiki and Ishii (Hibiki and Tshii, 2003) developed a mechanistic model for the active nucleation site density on boiling surfaces utilizing the size and cone angle distributions of cavities on the surface. They improved upon the model of (Yang and Kim, 1988) by adopting physically sound approximations for cavity radius and cone angle distributions. They validated their model by comparing predictions with the boiling data from numerous investigators obtained from surfaces of different materials and experimental fluids, covering a wide range of contact angles, pressure flow conditions, and number of bubbles produced. The model equations used by (Hibiki and Tshii, 2003) are:      
 λ θ2 Nn ¼ N n 1  exp  2 exp f ρþ 1 Rc 8μ

(4)


2
3
 f ρþ ¼ 0:01064 þ 0:48246ρþ  0:22712 ρþ þ 0:05468 ρþ

(5)




1. The effect of the submarine surface on the bubble growth is insignificant, or the bubble grows in an infinite liquid space. 2. The effect of neighboring bubbles on the growth rate of a bubble is insignificant. Relaxing these two assumptions will negate the applicability of Scriven's theory of spherically symmetric bubble growth, and the analysis is non-trivial. The model for the growth of a bubble from zero initial size can be described as follows. Let C ¼ ρρbρρ be the normalized mass concentration b

ρþ ¼ log ρ*

(6)

ρ  ρg ρ* ¼ l ρg

(7)



∂C ρ* R2 dR ∂C D ∂ 2 ∂C þ ¼ r ∂t 1 þ ρ* r 2 dt ∂r r 2 ∂r ∂r

2γ Pa S

 (9)

where RðtÞ denotes the time-dependent radius of the bubble. Two boundary conditions and two initial conditions are required to solve Equation (9):

where Rc is the critical bubble radius for nucleation, N n ¼ 4:72*105 sites/m2 from (Hibiki and Tshii, 2003) (used as m2 since counting sites is dimensionless), θ is the liquid-solid contact angle at the liquid-vapor interface (rad), μ ¼ 0:722 rad, λ ¼ 2:5*106 m, ρl is the saturated liquid density of the pure solvent (methane, ethane, or their mixture) and ρg is the saturated vapor density of the pure solute (nitrogen). Densities are determined by using Equation (2) at a given sea depth and then using the REFPROP mixture equation of state. Equations (4)–(7) are used together with an expression for critical radius to compute nucleation site density, in m2. Hibiki and Ishii (Hibiki and Tshii, 2003) originally used this set of equations for boiling; here the same model equations are used for effervescence, but the critical radius is changed to represent effervescence of dissolved gas: Rc ¼

i

of the dissolved gas in the liquid, where ρ is the liquid mixture density. Subscripts b and i denote bulk liquid and surface values. Denoting the diffusivity of the solute (nitrogen gas) in the solvent (liquid methane/ ethane) by D, the governing equation for solute transport, assuming spherical symmetry, is:

Rðt ¼ 0Þ ¼ 0

(10)

Cðt ¼ 0; rÞ ¼ 0

(11)

Cðt; r ¼ RðtÞÞ ¼ 1

(12)

Cðt; r → ∞Þ ¼ 0

(13)

The initial conditions state that the critical radius of the bubble is zero, and the concentration is initially zero. The boundary conditions state that the concentration at the bubble radius is 1, and that far away from the bubble surface, the concentration goes to 0. The solution for the solute concentration field is obtained using a similarity solutionCðr; tÞ ¼ FðηÞ, where the similarity variable is:

(8)

r 2 Dt

η ¼ pffiffiffiffiffi

where γ is the surface tension of the pure solvent, Pa is the ambient pressure (written in Pa, not MPa, to resolve units to meters), and S is the supersaturation (Equation (1)). Note that Qi and Klausner (2006) use a similar expression for critical radius to predict nucleation site density of dissolved gas. Curiously, they start with the same formulation as (Yang and Kim, 1988) but use different approximations for the cavity size and cone angle distributions. The constants in these probability distributions are specific to the materials comprising the nucleation surfaces, and not universal as in the model proposed by (Hibiki and Tshii, 2003).

(14)

Scriven defines the growth of bubble radius as: pffiffiffiffiffi RðtÞ ¼ 2β Dt

(15)

where β is the bubble growth constant. Rewriting Equation (9) using the similarity variable, and plugging in Equation (15), Equation (9) becomes:   1 2β3 ρ* 2 dξ 2 exp  ξ  ξ 1 þ ρ* η ξ FðηÞ ¼ Z ∞   3 * 1 2β ρ exp  ξ2  dξ ξ 1 þ ρ* ξ2 β Z

2.4. Bubble growth Scriven (1959) developed a model for the spherically symmetric growth of a bubble in a liquid of infinite extent driven by both heat and mass transfer. For heat transfer driven growth, the liquid is assumed to be superheated, in which a vapor bubble in introduced whose surface temperature is the equilibrium saturation temperature at the prevailing pressure. The bubble is assumed to grow from zero size. For mass transfer driven growth, the liquid is assumed to be a two component mixture (liquid solvent and vapor solute), with the solute diffusing into the bubble and causing it to grow. In this work, Scriven's theory is restricted to the mass transfer driven growth of a single component gas. The diffusion controlled growth of multicomponent gas bubbles was later

∞

(16)

Equation (16) will yield the concentration as a function of time and radial distance for a given constant β. To determine the value of β, a mass balance is taken at the vapor/liquid interface: 



ρg  ρi 1 

ρ*

1 þ ρ*



dR ∂C ¼ ðρb  ρi ÞD dt ∂r r¼RðtÞ

which yields the following expression for β:

4

(17)

J. Hartwig et al.

2β2

Z

1 0

Planetary and Space Science xxx (xxxx) xxx

" exp  β2

!#

 *  ρ 2 x1 1 þ ρ* ð1  xÞ 1

2

dx ¼

ρb  ρi

ρg  ρi 1  1þρ ρ* *

 (18)

Using Equations (1) and (7), the right hand side of Equation (18) can be rewritten as: ρbρρi ¼ ρρl xi S ¼ ð1 þ ρ* Þxi S. Therefore, Equation (18) g

g

becomes: 2β

2

Z

1 0

" exp  β

2

!#  *  1 ρ ð1 þ ρ* Þxi S

 2 x  1 dx ¼ * 1 þ ρ* ð1  xÞ2 1  ρi 1  ρ * ρg

Fig. 3. The bottom-view dimensions of the rectangle considered in the moving case. The projected area Wsub Lsub (see Table A1 for the submarine dimensions) is used in place of skin area.

1þρ

ð1 þ ρ* Þxi S ¼ 1  ρρi

For bubbles that nucleate at a position 0  x'  xon the surface section, the growth time before reaching position x is:

l

(19) For given values of ρb , ρi , and ρg , or equivalently for given values of ρ* , xi , and S, Equation (19) can be solved numerically for β.

tx' ¼

x' vsub

(20)

Thus, the maximum bubble diameter is estimated by the growth time sub . The sum of the volume of all bubbles that cross position x is tx' ¼ Lvsub

2.5. Bubble volume and area coverage modeling

determined by integrating the bubble volume due to each width “strip” Nn Wsub dx' from the submarine nose to position x:

Now, models are needed for how the bubbles accumulate and move along the sub to determine both the volume gas fraction in the sea and also bubble area coverage along the submarine. General assumptions used in this section are:

Vbub ðxÞ ¼

4π Nn Wsub 3

Z

x 0

ðRðtx' ÞÞ3 dx'

(21)

where the volume of each bubble is that of a sphere with radius from Equation (15) at time tx' . Evaluating this integral yields:

 Nucleation sites are uniformly distributed according to Equation (4).  The production of bubbles at each site is a continuous stream.  The growth of each bubble is spherically symmetric and not influenced by the surface it nucleated from, or by other bubbles, in line with Scriven's assumptions. Bubble growth is therefore given by Equation (15) whether the bubble is moving or stationary. This is a conservative approach and will essentially lead to overestimating the amount of gas.  The bubble area used in area coverage calculations is the frontal area of a (bubble) sphere facing a flat surface (the submarine), a circle of the bubble radius. In other words, each bubble projects an area of π R2bubble on the bottom of the submarine. This is also conservative, since some bubbles will block the full projection of other bubbles onto the sub surface.  In the moving case, bubble buoyancy is neglected and all bubbles are swept uniformly at the specified submarine velocity to the aft end. This is a conservative approach and over-predicts bubble coverage because it ignores how some bubbles will rise up the submarine side due to buoyancy.  In the nonmoving case, buoyancy will however, carry bubbles away that form on the top side of the sub; therefore only bubbles that form on the bottom side are considered in calculations.  The projected area of the bottom of the submarine is used for calculations, rather than the actual area; this again over-predicts bubble coverage.  The influence of the boundary layer is neglected; the flow over the skin is at vsub .

Vbub ðxÞ ¼

64π Nn Wsub D3=2 β3 x5=2 3=2

15vsub

(22)

The units in this equation are m3. In the moving case, the most important location for bubble volume accumulation is at the aft end, at the entrance region of the propellers. The volume occupied in which a strip of bubbles resides is taken to be a rectangular box as shown in Fig. 4. The bottom surface bubbles are assumed to occupy a cross-section equal to the height times the width of the submarine, and a length equal to the propeller length. Bubbles coming off the top surface are assumed to rise above the height of the propellers and are therefore ignored; bubbles from the bottom surface are assumed to distribute evenly over the square Wsub Lsub , reflecting the varying height of surface elements on the submarine bottom, as opposed to concentrating in the propeller. The fraction of that reference volume occupied by bubbles is what affects the operation of the propellers, and is

2.5.1. Moving case First the moving case is addressed. At each nucleation site, new bubbles form and accumulate with bubbles that formed at sites upstream. The amount of accumulation at any location is determined by the submarine velocity, as well as the nucleation and growth rates. Each bubble has had time to grow, according to the distance from its nucleation site to the position of calculation, and the total bubble coverage at any position is the sum of the coverage of all individual bubbles. Illustrated in Fig. 3 is the coordinate system and dimensions used for computations. Bubbles move in the positive x direction at speed vsub .

Fig. 4. The reference volume considered to be occupied by bubbles after they come off the bottom of the submarine. 5

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

given by: Volume fraction ¼

Vbub ðLsub Þ 64π Nn D3=2 β3 L5=2 sub ¼ 3=2 Lprop Hsub Wsub 15vsub Hsub Lprop

(23)

This is a ratio of volumes, and is therefore dimensionless. The submarine width Wsub cancels out of this ratio because the bubbles from each infinitesimally-thin strip (based on projected area) down the length Lsub occupy, at the aft end, a thin box of cross-section Hsub by Lprop ; whether more or fewer such strips are accumulated to reach Wsub is immaterial in this model. The total area covered of the bubbles at position x along the bottom of the submarine at any moment in time is computed similarly to the bubble volume. The main difference is that, instead of multiplying by the volume per spherical bubble, one multiplies by the cross-sectional area per bubble: Abub ðxÞ ¼ π Nn Wsub

Z 0

x

ðRðtx' ÞÞ2 dx'

Fig. 6. The reference area considered to be occupied by bubbles after they come off the bottom of the submarine.

2.5.2. Quiescent case The bubbles generated when the submarine is stationary in the Titan seas are driven upward both by buoyancy relative to the liquid around them and by the bulk flow of that liquid due to natural convection. Such convection is induced by the submarine waste heat. Adding the bulk velocity of the liquid due to natural convection to the relative velocity of the bubble in the liquid gives the bubble velocity relative to the submarine. For the quiescent case, two cases are examined: the distribution of bubble area coverage up either side of the submarine over the instrument section, and the volume fraction of bubbles that accumulate around the top propellers before they turn on. Fig. 7 shows the 1D bases for area and volume estimates. The speed of the flow up the submarine side due to natural convection is based on a laminar similarity solution for a vertical flat plate (Incropera and DeWitt, 2002):

(24)

Evaluating this integral with the radius from Equation (15) gives, in units of m2, Abub ðxÞ ¼

2π Nn Wsub Dβ2 x2 vsub

(25)

For illustration, the functional form of area coverage across the submarine bottom surface, from forward to aft end, is that of a parabola (Fig. 5) with the highest estimated bubble coverage at x ¼ Lsub . Therefore, any science instruments that are sensitive to high bubble concentrations (e.g. cameras, sonars) should be installed as far toward the nose of the submarine as possible. The reference area to which the bubble area is compared is shown in Fig. 6, the strip with width equal to the submarine and length equal to the propeller length. This length is chosen to be consistent with how volume fraction is estimated. Therefore, the fraction of bubble area to reference area is: Area fraction ¼

Abub ðLsub Þ 2π Nn Dβ2 L2sub ¼ Wsub Lprop vsub Lprop

1=2

vnc ¼

2ξf GrHsub Hsub

(27)

where GrHsub ¼

gβl ðTskin  Tsea ÞH 3sub ξ2

(28)

By numerical approximation, the peak values of f (based on a function arising in the solution of the differential equation) depend on Prandtl number. f ¼ maxf 'ðζÞ ¼ 0:2119Pr0:302 (Pr > 0:5) This speed estimate will be considered constant in both space and time, neglecting boundary layer growth. Buoyancy is not neglected in this case, it is added at the end after assuming vnc ¼ 0 for the solution process. Error! Reference source not found. depicts the forces on a bubble in a free stream at any time t with the reference frame chosen so the sea far from the submarine is at rest. Let y point vertically up with y ¼ 0 the point where the bubble nucleates; y ¼ 0 is not fixed relative to the submarine, because different bubbles nucleate at different points on the skin and thus grow to different sizes

(26)

This is an area ratio and is therefore dimensionless. The submarine width Wsub cancels out for the same reason it cancels in the volume computation. Any area can be inserted into the numerator in Equation (26) to determine the amount of bubble coverage at any particular location.

1

0.4

a)

A

bub

bub

sub

(L )

0.6

(x)/A

0.8

0.2 b)

0

0

0.2

0.4

x/L

0.6

0.8

Fig. 7. For the quiescent case, a) the area coverage of bubbles over the side of the submarine is based on half the vehicle width (bubbles traveling up one side do not interact with bubbles traveling up the other side) b) the volume fraction of bubbles around the propellers is based on the projected submarine height. An example panel for calculation is shown in blue in each case.

1

sub

Fig. 5. Relative bubble area concentration on the bottom surface as a function of distance from submarine nose for the moving case. 6

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

pushes liquid out of the way is:

before leaving the submarine side. The bubble buoyancy is driven by the density difference between gas and liquid over the assumed spherical volume (Fig. 8):
 4π Fb ¼ g ρl  ρg RðtÞ3 3

  2    1 4π d y 3 4π dvnc 1 F1 ¼  ρ l RðtÞ3 RðtÞ3 þ ρl vnc 2 þ ρl 2 2 2 3 3 dt dt    dy d 4π RðtÞ3  dt dt 3   2π d 2 y π RðtÞ3 dy ¼  ρl RðtÞ3 2 þ ρl vnc  dt 3 t dt

(29)

where ρl is the mixture liquid density and ρg is the density of the saturated nitrogen vapor in the bubble. The resistance is assumed to come from drag due to flow around the bubble, and the associated area is a plane circle with the radius of the bubble:  2 1 dy  vnc CD π RðtÞ2 Fdr ¼ ρl 2 dt

The resulting equation of motion is:
 d2 y 9μ þ 3ρl β2 D 1 dy 2 ρl  ρg 9μ þ 3ρl β2 D vnc
  ¼ þ gþ
2 2 dt ρl þ 2ρg 2ρl þ 4ρg β D t dt 2ρl þ 4ρg β2 D t

(30)

area π RðtÞ2 . The drag coefficient used is that for shear-free flow (of the liquid over the bubble), CD ¼ 48=Re2R (Harper, 1972), where the Reynolds number of the individual bubble is based on diameter:
dy dt

  vnc 2RðtÞ ξ

yðtÞ ¼

tcr ¼

4π d2 y RðtÞ3 ρg 2 dt 3

(36)

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi# "

  9μ þ 5ρl þ 4ρg β2 D 8 ρl  ρg β2 D
 2  2 gH
 vnc þ v2nc þ 4g ρl  ρg β D 9μ þ 5ρl þ 4ρg β D (37)

(32)

Note that ð0Þ ¼ vnc as the bubble is swept up in the flow. The method of turning the time tcr from Equation (37) into bubble area or volume estimates is similar to that used in the moving case. The submarine sides are discretized into n panels with representative length L' (which cancels out later) and area 1n L'W2sub for area coverage and 1 n L'ð1:3Hsub Þfor volume coverage calculations, and the bubble rise time is calculated from the center of each panel. The number of bubbles from each panel is based on the site density given by Equation (3); using tcr and the radius Rðtcr Þ from Equation (15) and, the areas and volumes of bubbles from each section are added to estimate the total coverage at any point. The area fraction of bubbles is taken relative to the width of the side scan sonar (1/8 of the submarine width); the volume fraction of bubbles is taken relative to a rectangular prism based on the propeller size. Therefore the area and volume fractions are: dy dt

The acceleration of the bubble mass (based on gas density) through the liquid is expressed by Newton's second law: Fb  Fdr þ F1 ¼


 2 ρl  ρg β2 D
 gt 2 þ vnc t 9μ þ 5ρl þ 4ρg β2 D

The time tcr for the solution to reach a height y ¼ H (taking the positive root of the quadratic equation) is

(31)

and ξ is the kinematic viscosity of the sea liquid. Therefore the drag force on the bubble is   dy  vnc Fdr ¼ 12πμRðtÞ dt

(35)

with initial conditions yð0Þ ¼ 0 and dy dt ð0Þ finite. The analytical solution that fits these criteria, by the integrating factor method, is:

where CD is the drag coefficient of a sphere based on the projected frontal

Re2R ¼

(34)

(33)

The “added mass” force (Brennan, 2005) as the bubble grows and

n P

Area fraction ¼

i¼1

Nn L'W' *π ðRðtc ÞÞ2 n L'W8sub n P

Volume fraction ¼

i¼1

¼

Nn L'H'n *43πðRðtc ÞÞ3 L'Hprop Wprop

n 8π Nn W' X ðRðtcr ÞÞ2 nWsub i¼1

¼

n 4π Nn H' X ðRðtcr ÞÞ3 3nHprop Wprop i¼1

(38)

(39)

In this computation the heights used are W' ¼ 0:5Wsub (area coverage) and H' ¼ 1:3Hsub (volume fraction) according to Fig. 7. For comparison, the detailed solution is checked against a quasisteady solution obtained by assuming no acceleration and no bubble growth Fb ¼ Fdr . When this expression is expanded and integrated from yð0Þ ¼ 0, the result is: yðtÞ ¼


 2 ρl  ρg β2 D 2 gt þ vnc t 9μ

(40)

The time tcr to reach y ¼ H (taking the positive root of the quadratic equation) is 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi3
 8 ρl  ρg β2 D 9μ 2 4  vnc þ v þ  tcr ¼
gH 5 nc 9μ 4g ρl  ρg β2 D Fig. 8. Forces on a gas bubble rising through liquid. 7

(41)

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

    1 1 Δv Δv τ ¼ ρl v2sub Afront * CD ¼ ρl Ap vsub þ 2 4 2

The deviation in area and volume fraction between Equations (37) and (41) is less than 3% within the range of submarine operation. The detailed solution, Equation (37), is used for the figures below.

(42)

where Ap ¼ 4π H 2prop is the area of the propeller disk. Since CD ¼ 0:4 and Afront ¼ 0:5 m2 for the submarine (Oleson et al., 2014), one can solve for the required velocity change:

2.6. Pressure drop in propellers The last major component of effervescence on the submarine is that due to the pressure drop through the propellers. Boiling-based cavitation is not considered because the Titan seas are near freezing. Such a pressure drop reduces nitrogen solubility further than the thermally driven case alone, and as a result allows more bubbles to come out of solution at a nucleation site. This nucleation occurs on the aft surface of the propeller blade, and any bubbles so formed travel at the forward speed of the submarine (Fig. 9). The bubbles travel at the speed vsub through half the propeller length, while not interacting with any walls around the propeller (so bubbles come from the propellers exclusively). The pressure is assumed to decrease by an amount P0  P1 between the freestream and the propeller surface. The surface area for nucleation used for this estimate is some fraction σ of the total disk area 4π H 2prop , which is determined by the propeller blade size. The solidity of the propeller, σ , accounts for the fact that bubbles will only nucleate on the blade and not in the open area portion of the propeller cavity; at fixed Hprop , σ is larger for more/ wider blades (i.e. looks more like a solid disk). It is also assumed that bubbles do not nucleate on the side walls of the cavity, behind the blade, because the pressure drop is concentrated around the hydrodynamic surface of the propeller itself. The pressure P1 is estimated according to the required thrust of the submarine using the Bernoulli method upstream of the propeller. The thrust is assumed to divide equally among the 4 propellers, while the sum of all propellers equals the vehicle drag. It is also related, by impulsemomentum, to the total velocity change Δv through the propeller:

Δv ¼ vsub

1þ

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! Afront 1 þ CD 4Ap

(43)

Finally, the pressure drop corresponding to half of such velocity change (half of the total cavity length) is approximated by the Bernoulli relation along a streamline leading to the propeller: " # 2 Δv 2  vsub vsub þ 2 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! 1 Afront ¼ ρl v2sub 4  1 þ 1 þ CD 2 4Ap sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2 3 1 Afront 5  1 þ 1 þ CD þ 4 4Ap

1 P0  P1 ¼ ρl 2

(44)

The first step in the estimate of additional volume is to obtain the nucleation site density of new bubbles on the low-pressure side of the propeller. The equilibrium nitrogen solubility xN2 is updated for the reduced pressure, at the skin temperature, then used to get supersaturation by Equation (1). Substituting the pressure P1 into Equation (8), with the new supersaturation, gives the new critical radius. Finally, the new critical radius goes into Equation (4) to estimate the total nucleation site density at the reduced pressure, Nn;p . This includes, however, the nucleation sites that would already be active at the sea pressure, so the count of new bubbles behind the propeller is based on the difference between low- and high-pressure concentrations:
 π #bubbles ¼ σ H 2prop Nn;p  Nn 4

(45)

All the new bubbles, which nucleate only on the propeller blade, have a growth time equivalent to the time taken to travel half the length of the housing. In this time, the bubbles grow to a size determined by Equation (15), and integrating such growth along the streamline from propeller surface to open sea (see Fig. 10) gives the volume of this line of bubbles in m3: Z Vline ¼

0

0:5Lprop

π ðRðtx ÞÞ2 dx ¼

Z 0

0:5Lprop

 

π R

x

vsub

2 dx ¼

π β2 DL2prop 2vsub

(46)

Finally, the volume to which the bubble stream is compared, is the disk filling the aft half of the propeller casing: 0:5Lprop 4π H 2prop . Adding the total volume fraction of the bubbles behind the propeller to the volume fraction of bubbles at the propeller blade, the result in the moving case is

Fig. 9. A diagram of the model used for pressure-based effervescence, with a plot of the assumed pressure distribution.

Fig. 10. An illustration of how bubbles grow in the propeller cavity. 8

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

Volume fractionp ¼ Volume fractionjat blade þ


5=2

¼

64π Nn D3=2 β3 Lsub 3=2

15vsub Hsub Lprop

þ

#bubbles*Vline 0:5Lprop π4H 2prop

the submarine has waste heat of 3800 W and skin area Askin ¼ 10:34 m2 , for an average heat flux of 370 W/m2. At this value of heat flux, the skin temperature is expected to be less than 96 K (see Error! Reference source not found.). Note, however, that this is an average value – typically on real hardware there are sites (e.g. instrument penetrations through insulation) where local heat fluxes may be a factor of several higher than this. First, Fig. 11 presents submarine-specific values of skin temperature, as a function of waste heat flux. Fig. 11 shows that, in comparing moving to quiescent cases, forced convection clearly cools the submarine skin more effectively than natural convection. The skin is also noticeably warmer in ethane, because in the lower viscosity of a methane sea Reynolds and Rayleigh numbers are substantially higher, which leads to higher heat transfer coefficients. The solubility limit in Equation (2) is shown graphically in Fig. 12 for two bounds on the methane to ethane mole ratio: 5 to 95 and 85 to 15 (the latter is based on (Hayes, 2016) to represent Ligeia Mare, the former is a generic low bound to represent Kraken Mare). The trends are that solubility rises with increasing pressure, decreasing temperature, and increasing methane mole fraction. If the bulk liquid is already saturated, heating the submarine surface will reduce solubility and cause supersaturation. Fig. 13 plots supersaturation as a function of submarine skin temperature and liquid pressure for two different sea concentrations. Clearly, a hotter surface generates a greater excess of gas over the solubility limit. The reason for the substantially greater supersaturation at the bottom of Fig. 13b is that, at lower pressures, the solubility values are lower. Therefore, the difference in solubility is greater relative to the value at the skin temperature, which raises the fractional supersaturation. The supersaturation of the liquid mixtures is relatively small, but the surface could generate a high concentration of bubbles anyway depending on the contact angle between the liquid and the submarine surface Fig. 14. Next, nucleation site density as a function of skin temperature and contact angle is plotted in Fig. 14. Contact angle is highly dependent on surface properties of the submarine and results show that the number of nucleation sites is highly sensitive to contact angle. While the contact angle between cryogenic liquids and metallic surfaces is expected to be small (Hartwig and Mann, 2014), it is still necessary to quantify the effect of non-zero contact angle. Results show that the number of sites changes by several orders of magnitude for a contact angle of 15 . Comparing plots, the number of sites is higher in the methane rich sea over the ethane rich case. Bubble growth rate is plotted in Fig. 15 (changed from meters to mm to display fewer zeros). Bubbles grow more quickly at warmer skin temperatures and higher methane concentrations because both conditions raise the supersaturation of sea liquid against the submarine skin. This in turn raises the value of the growth constantβ. Figs. 16 and 17 plot the volume fraction (at the propeller inlet) and area coverage, respectively. Table 3 presents results of simulations at



πσ Nn;p  Nn Dβ2 Lprop

(47)

vsub

In the quiescent case, the only flow through the propellers is entrained flow, which is more specific than this model covers; this work does not give volume fraction behind the propeller in such cases. If the entrained flow speed is known, this speed can be used in place of vsub with the volume fraction from Section 2.5.2. 2.7. Method of solution The method of solution is as follows: 2.7.1. Thermally driven effervescence 1. Choose whether the submarine is quiescent or moving. This specifies the Nusselt number correlation to use for skin temperature. For the moving case, also specify the submarine velocity. 2. Choose a sea (e.g. Ligeia Mare, Kraken Mare) and a location within the seas. This determines P, Tbulk , xCH4 , and xC2H6 ; these four variables uniquely specify ρ* and D. 3. Use Appendix A to compute the submarine skin temperature as a function of waste heat flux. 4. Determine the critical radius from Equation (8), then use Equation (4) to determine the number of nucleation sites per square meter. 5. Use Equation (2), with the coefficients in Table 1, to determine xi and xb , and thus S from Equation (1). 6. Solve Equation (19) for the bubble growth constant β. Finally estimate bubble area and volume fractions due to thermal effects. In the moving case, these are given by Equations (26) and (23), respectively. In the quiescent case, follow the numerical solution of Section 2.5.2, leading to Equation (38) for area and Equation (39) for volume. 2.7.2. Pressure driven effervescence In the moving case, include the volume fraction effect of the pressure drop through the propellers according to Equation (47). 3. Numerical results and discussion Simulations were run for two different sea concentrations representative of Ligeia and Kraken Mare. The diffusion coefficient value D used in the calculations is quoted in Table A1. Ten panels are used in the quiescent case at a time step of t ¼ 0:01 seconds which was shown to achieve a grid and time step independent solution. The present design of 104

Quiescent Moving, 0.5 m/s

Quiescent Moving, 0.5 m/s

99

[K] skin

98

97

T

T

[ K]

100

skin

b) 85% Methane

101

a) 5% Methane

102

96

95

94 93 92

0

200

400

600

Q /A wh

800

1000

1200

1400

0

1600

2

skin

[W/m ]

200

400

600

Q /A wh

800 skin

1000

1200

1400

1600

[W/m2]

Fig. 11. Submarine skin temperature as a function of waste heat flux for a) ethane-rich sea and b) methane-rich sea, with sea temperature 93 K and given geometry. 9

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

Fig. 12. Equilibrium solubility of nitrogen in methane and ethane for a) ethane-rich sea and b) methane-rich sea. Color represents amount of dissolved nitrogen gas. Both figures are plotted on same color scale for comparison.

5

10

4

10

3

2

a) 5% Methane 10

1

10

8

10

7

10

6

10

5

10

4

10

3

10

2

n

10

2

10

N [1/m ]

n

2

N [1/m ]

Fig. 13. Supersaturation at submarine skin for a) ethane-rich sea and b) methane-rich sea, with sea temperature 93 K.

1 Degree 15 Degrees

93

94

95

96

T

97 skin

98

99

b) 85% Methane 1 Degree 15 Degrees

10 1

100

93

94

95

[K]

96

T

97 skin

98

99

100

[K]

Fig. 14. Nucleation site density as a function of skin temperature and contact angle for a) ethane-rich sea and b) methane-rich sea, with sea temperature 93 K.

7

0.6

a) 5% Methane 0.5

0.3 0.2

4 3 2

0.1 0

T - 94 K skin T - 96 K skin T - 98 K skin

5

R(t) [mm]

R(t) [mm]

0.4

b) 85% Methane

6

T - 94 K skin T - 96 K skin T - 98 K skin

1 0

10

20

30

40

50

0

60

Time [s]

0

10

20

30

40

50

60

Time [s]

Fig. 15. Bubble radius (in mm) at several skin temperatures for a) ethane-rich sea and b) methane-rich sea, with sea conditions 93 K and 0.4 MPa. 10

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

Fig. 16. Volume fraction of bubbles around the propellers (log scale, at the propeller inlet) for a) ethane-rich quiescent case, b) ethane-rich moving case, c) methanerich quiescent case and d) methane-rich moving case. For all plots, sea temperature is 93 K and θ ¼ 150 .

Fig. 17. Highest computed area fraction of bubbles (log scale) on the submarine skin for a) ethane-rich quiescent case, b) ethane-rich moving case, c) methane-rich quiescent case and d) methane-rich moving case, with sea temperature 93 K and θ ¼ 15 degrees.

the propeller Lsub Wprop , and for the quiescent case, an area from the bottom centerline of the submarine to the end of the ballast tank, using any representative length, L'W8 (from the side scan sonar). Results in Table 3 show that effervescence has a much greater propensity in methane-rich seas than it is in ethane-rich seas (skin

various vehicle velocities, contact angles, and sea temperature and pressure (which simulates the effect of depth). For both the moving and non-moving case, the volume fraction is computed for a rectangular cube of liquid the size of a propeller Lprop Hprop Wprop . The area fraction is taken at the point where gas coverage is maximum. For the moving case, this occurs at the aft end of the submarine using an area equal to the area of 11

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

Table 3 Results of the effervescence model for several cases of interest to the Titan submarine. vsub (m/s)

Methane fraction

Qwh (W/m2) Askin

P (MPa)

θ (deg)

Tskin (K)

Nn (1/m2)

Area fraction

Volume fraction (before prop)

Volume fraction (after prop)

0 0 0 0 0.5 1 0.5 1 0.5 1 0.5 1 0 0 0 0 0.5 1 0.5 1 0.5 1 0.5 1

0.85 0.85 0.05 0.05 0.85 0.85 0.85 0.85 0.05 0.05 0.05 0.05 0.85 0.85 0.05 0.05 0.85 0.85 0.85 0.85 0.05 0.05 0.05 0.05

370 370 370 370 370 370 370 370 370 370 370 370 370 370 370 370 370 370 370 370 370 370 370 370

0.15 0.3 0.15 0.3 0.15 0.15 0.3 0.3 0.15 0.15 0.3 0.3 0.15 0.3 0.15 0.3 0.15 0.15 0.3 0.3 0.15 0.15 0.3 0.3

1 1 1 1 1 1 1 1 1 1 1 1 15 15 15 15 15 15 15 15 15 15 15 15

95.5 95.5 96.3 96.2 93.7 93.4 93.6 93.4 94.1 93.6 94.0 93.6 95.5 95.5 96.3 96.2 93.7 93.4 93.6 93.4 94.1 93.6 94.0 93.6

34.0 83.7 21.3 43.7 6.21 3.38 11.27 6.00 5.47 3.01 9.70 5.21 7.60E03 1.87E04 4.76E03 9.75E03 1.39E03 756 2.52E03 1.34E03 1.22E03 671 2.16E03 1.16E03

6.94E-05 0.000244 4.24E-06 1.53E-05 2.56E-05 3.33E-06 9.65E-05 1.16E-05 1.66–06 2.47E-07 5.13E-06 7.38E-07 0.0133 0.0545 9.47E-04 0.0034 0.0057 7.44E-04 0.0215 0.0026 3.69E-04 5.51E-05 0.0011 1.64E-04

6.65E-07 3.36E-06 2.08E-09 9.71E-08 1.13E-08 1.04E-09 4.27E-08 3.64E-09 7.32E-10 7.73E-11 2.27E-09 2.31E-10 1.49E-04 7.50E-04 4.64E-06 2.17E-05 2.53E-06 2.33E-07 9.53E-06 8.13E-07 1.63E-07 1.73E-08 5.07E-07 5.14E-08

NA NA NA NA 1.14E-08 1.08E-09 4.28E-08 3.76E-09 7.35E-10 8.02E-11 2.28E-09 2.37E-10 NA NA NA NA 2.54E-06 2.41E-07 9.56E-06 8.39E-07 1.64E-07 1.79E-08 5.08E-07 5.28E-08

residence time to grow. Eventually as velocity is further increased, pressure drop through the propellers increases, causing more nucleation sites, and thus more bubbles to come out of solution. Therefore, operationally for the submarine, there is an optimal velocity range which minimizes both the low and high speed limits. Finally, this model can be compared to recent data on the heat flux values that trigger effervescence in similar mixtures (Richardson et al., 2018) in Fig. 20. The submarine operating point, about 370 W/m2, is an order of magnitude less than the data, suggesting a safety factor of at least 5 for the present design. As noted in (Lorenz, 2016), fluxes of the order of 90 W/m2, as encountered with the Huygens probe, are enough to cause detectable environmental perturbations on Titan, but considerably larger fluxes (~50,000 W/m2) are needed to cause actual boiling of liquid methane.

temperature being equal). This is attributed to the fact that methane rich seas have a higher supersaturation than ethane rich seas, which amplifies both the number of nucleation sites and the bubble growth rate. Comparing quiescent to moving cases, the highest accumulations are in the quiescent cases. The highest area coverage for the cases in Table 3 is 5.5%, and the highest volume fraction is 0.0758%. These extremes occur in the quiescent, methane-rich sea, high pressure, θ ¼ 15 degree case; other cases (especially where θ ¼ 1 degree) are orders of magnitude smaller. The gas volume fraction after the propellers is determined using baseline propeller parameters. Based on pictures of the submarine (Oleson et al., 2014), the solidity of the propellers is about 0.5. The results of this assumption are shown in Fig. 18. Compared to Fig. 16b and (d), the difference is small (on the order of 1% of the value before the propellers). The conclusion that methane produces more bubbles than ethane still holds. Fig. 19 summarizes operational limits for the submarine in Kraken Mare and Ligeia Mare, respectively by plotting gas void faction at the aft end of the propellers as a function of submarine velocity. Comapring the two plots, for the same velocity, void fraction is higher in the methane rich sea over the ethane rich sea, is as expected. At zero speed, volume fraction is at the thermally driven limit where only natural convection sweeps away bubbles that form. As velocity increases, volume fraction decreases in two ways: bubbles that nucleate do not have enough

4. Conclusion This paper presented a combined model by which the incipience and growth of bubbles around a submersible body may be approximated. It uses a nucleation site model adapted from boiling studies for effervescence bubbles. The growth of these bubbles is based on undisturbed spherical symmetry and bubbles accumulate around the submarine. The results show that, depending on heat flux, Titan sea composition and pressure, bubble fractions may occupy substantial area and volume

Fig. 18. Volume fraction after the propellers, moving at 0.5 m/s, for a) ethane-rich sea and b) methane-rich sea. 12

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

10

-4

10

-5

1 Degree Contact Angle 15 Degree Contact Angle

-6

10

-7

10

-8

10

-9

10

-10

10

-11

Void Fraction

Void Fraction

10

Kraken Mare

0

5

10

15

20

10

-3

10

-4

Ligeia Mare

10

-5

1 Degree Contact Angle 15 Degree Contact Angle

10

-6

10

-7

10

-8

10

-9

10

-10

10

-11

0

5

Speed [m/s]

10

15

20

Speed [m/s]

10

5

10

4

2

Waste Heat Flux [W/m ]

Fig. 19. Volume fraction after the propellers as a function of velocity for Kraken Mare (ethane rich sea) at 400m depth and Ligeia Mare (methane rich sea) at 150m.

Data, Richardson et. al. Submarine Operating Point

1000

100

0

0.05

0.1

0.15

xN

0.2

0.25

0.3

2

Fig. 20. Comparison of submarine operating point to experimental heat flux values that triggered effervescence.

small change between Figs. 18b and 16d suggests that the expected pressure drop, given the propeller size, does not substantially impact performance concerns at the lower anticipated operation speeds of the submarine.

fractions. These bubble effects are large enough to potentially impact submarine design. The highest area coverage for the quiescent case was approximately 5% at depth in Ligeia Mare. However, the larger area fractions that are calculated for the moving case are a conservative estimate based on neglecting buoyancy-induced sweeping of bubbles off the sides of the submarine. Even so, a methane-rich sea might still cause difficulty (much more so than an ethane-rich sea). If bubble accumulation is an issue, excess heat may be concentrated at a point on the skin away from propellers and instruments. With respect to propellers, the

Acknowledgement This work was funded by the NASA Innovative Advanced Concepts Office under the Space Technology Mission Directorate at NASA.

Appendix A. Submarine Skin Temperature Model The parameters used in calculating the submarine skin temperature are shown in Table A1, and adopted from (Oleson et al., 2014). The diffusion coefficient of nitrogen in methane is the only fluid property taken as a constant, because few studies of it exist in the literature. All other properties are computed from the software REFPROP.

13

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

Table A.1 Parameters for computing skin temperature and general submarine dimensions. Quantity

Units

Value

Propeller length Lprop Propeller width Wprop Propeller height Hprop Length Lsub Width Wsub Height of main body (excluding communication array) Hsub Skin area Askin Skin temperature Tskin Sea temperature Tbulk Travel speed vsub Waste heat Qwh Titan gravity g Convection coefficient h Diffusion coefficient of nitrogen in liquid methane D

m m m m m m m2 K K m/s W m/s2 W/m2-K m2/s

0.117 0.14 0.14 6.54* 1.23* 0.78* 10.34 Variable 93 Variable Variable 1.35 Variable 2*109 Mesli et al. (2011)

*Resized from the original design to allow for operation in methane-rich seas.

The equilibrium of the waste heat flux with convection heat flux in the Titan seas, determines the skin temperature of the vehicle. The latter heat flux is expressed as Q_ conv ¼ hðTskin  Tbulk Þ Askin

(A.1)

for the convection coefficient h. In the case where the submarine is quiescent, natural convection occurs, and the correlation used for the convection coefficient h is (Incropera and DeWitt, 2002): NuHsub ¼

" hHsub ¼ 0:6 þ h k

0:387Ra1=6 Hsub
0:5599=16 i8=27 1 þ Pr

#2 RaHsub < 1012

(A.2)

In the case where the submarine is in moving through the Titan Sea, turbulent forced convection occurs and the correlation used is (Incropera and DeWitt, 2002): NuLsub ¼

hLsub 4=5 ¼ 0:037ReLsub Pr1=3 ; 0:6 < Pr < 60 k

(A.3)

In these correlations, the Reynolds number is based on submarine length: ReLsub ¼

vsub Lsub ξ

(A.4)

while the Rayleigh number is based on submarine height: RaHsub ¼

gβl ðTskin  Tbulk ÞH 3sub Pr ξ2

(A.5)

The required liquid properties in these computations – thermal conductivity k, thermal expansion coefficient βl , kinematic viscosity ξ and Prandtl number Pr- are computed at the film temperature around the vehicle: 1 Tfilm ¼ ðTskin þ Tbulk Þ 2

(A.6)

Finally, the 1st law power balance between the submarine waste heat and Titan sea convection is: Q_ wh Q_ conv ¼ Askin Askin

(A.7)

where expressing the balance in terms of heat flux allows results for any geometric configuration. To determine the resultant skin temperature, Equation (A.7) is solved using an iterative, guess-and-check method: 1. Initiate a value for Tskin , around 100 K. 2. Solve the equations for all heat transfer terms. 3. Check the energy balance (Equation (A.7)) and iterate until equality is met. If Q_ wh < Q_ conv , then the skin temperature is too high. If Q_ wh > Q_ conv , then the skin temperature is too low.

14

J. Hartwig et al.

Planetary and Space Science xxx (xxxx) xxx

References

Liger-Belair, G., Marchal, R., Robillard, B., Vignes-Adler, M., Maujean, A., Jeandet, P., January 1999. Study of effervescence in a glass of champagne: frequencies of bubble formation, growth rates, and velocities of rising bubbles. Am. J. Enol. Vitic. 50, 317–323. Loomis, A.G., Walters, J.E., 1926. The vapor pressure of ethane near the normal boiling point. J. Am. Chem. Soc. 48 (8), 2051–2055. Lorenz, R.D., 2016. Heat rejection in the titan surface environment: potential impact on science investigations. J. Thermophys. Heat Transf. 30 (2), 257–265. Malaska, M., Hodyss, R., Lunine, J.I., Hayes, A.G., Hofgartner, J.D., Hollyday, G., Lorenz, R.D., June 2017. Laboratory measurements of nitrogen dissolution in titan lake fluids. Icarus 289, 94–105. Mesli, F., Mahboub, R., Mahboub, M., 2011. Molecular dynamics comparative study of methane-nitrogen and methane-nitrogen-ethane systems. Arabian Journal of Chemistry 4, 211–222. Methane, 2016. National Institute of Standards and Technology [Online]. Available: http ://webbook.nist.gov/cgi/cbook.cgi?ID¼C74828&Mask¼4&Type¼ANTOINE&Pl ot¼on. (Accessed 30 July 2016). National Institute of Standards and Technology, 2016. Nitrogen. National Institute of Standards and Technology [Online]. Available: http://webbook.nist.gov/cgi/ cbook.cgi?ID¼C7727379&Mask¼4&Type¼ANTOINE&Plot¼on. (Accessed 14 September 2016). Oleson, S., et al., 2014. COMPASS Final Report: Titan Submarine. NASA CD-2014-114. Prydz, R., Goodwin, R.D., January 1972. Experimental melting and vapor pressures of methane. J. Chem. Thermodyn. 4 (1), 127–133. Qi, Y., Klausner, J.G., 2006. Comparison of nucleation site density for pool boiling and gas nucleation. J. Heat Transf. 128, 13–20. Richardson, I.A., Hartwig, J.W., Leachman, J.W., 2018. Experimental PρT-x measurements of liquid methane-ethane-nitrogen mixtures to support an extraterrestrial submarine for titan. Fluid Phase Equilib. 462, 38–43. Scriven, L.E., 1959. On the dynamics of Phase growth. Chem. Eng. Sci. 10, 1–18. Wilt, P.M., 1986. Nucleation rates and bubble stability in water-carbon dioxide solutions. J. Colloid Interface Sci. 112, 530–538. Yang, S.R., Kim, R.H., 1988. A mathematical model of the pool boiling nucleation site density in terms of the surface characteristics. Int. J. Heat Mass Transf. 31, 1127–1135.

Air Solubility in Water [Online]. Available: http://www.engineeringtoolbox.com/air-sol ubility-water-d_639.html. (Accessed 13 September 2016). Arvelo, J.I., Lorenz, R.D., 2013. Plumbing the depths of Ligeia: considerations for depth sounding in Titan's hydrocarbon seas. J. Acoust. Soc. Am. 134, 4335–4350. Brennan, C.E., 2005. Fundamentals of Multiphase Flows. Cambridge University Press. Cable, M., Frade, J.R., 1987. Diffusion-controlled growth of multi-component gas bubbles. J. Mater. Sci. 22, 919–924. Carruth, G.F., Kobayashi, R., 1973. Vapor pressure of normal paraffins ethane through ndecane from their triple points to about 10 Mm Hg. J. Chem. Eng. Data 18 (2), 115–126. Edejer, M.P., Thodos, G., 1967. Vapor pressures of liquid nitrogen between the triple and critical points. J. Chem. Eng. Data 12 (2), 206–209. Ethane, 2016. National Institute of Standards and Technology [Online]. Available: http ://webbook.nist.gov/cgi/cbook.cgi?ID¼C74840&Mask¼4&Type¼ANTOINE&Pl ot¼on. . (Accessed 30 July 2016). Harper, J.F., 1972. The motion of bubbles and drops through liquids. Adv. Appl. Mech. 12, 59–129. Hartwig, J.W., Mann, J.A., February 2014. Bubble point pressures of binary methanol/ water in fine-mesh screens. American Institute of Chemical Engineers Journal 60 (2), 730–739. Hartwig, J.W., Colozza, A., Lorenz, R.D., Oleson, S., Landis, G., Schmitz, P., Paul, M., Walsh, J., 2016. Exploring the depths of Kraken Mare – power, thermal analysis, and ballast control for the Saturn Titan submarine. Cryogenics 74, 31–46. Hartwig, J.W., Meyerhofer, P., Lorenz, R., Lemmon, E., 2017. An analytical solubility model for nitrogen/ethane/methane ternary mixtures. Icarus 299, 175–186. Hayes, A.G., 2016. The lakes and seas of titan. Annu. Rev. Earth Planet Sci. 44, 57–83. Hibiki, T., Tshii, M., 2003. Active nucleation site density in boiling systems. Int. J. Heat Mass Transf. 46, 2587–2601. Incropera, F.P., DeWitt, D.P., 2002. Fundamental of Heat and Mass Transfer, fifth ed. John Wiley and Sons, New York. Jones, S.F., Evans, G.M., Galvin, K.P., 1999. Bubble nucleation from gas cavities - a review. Adv. Colloid Interface Sci. 80, 27–50.

15