© J. Markkanen and J. Agarwal 2020

1 Introduction

As comets approach the Sun, they start to eject gas and dust from their nuclei. The ejection of dust is mainly driven by the sublimation of volatile ices. Thermal modelling suggests that the sublimation of volatiles such as CO and CO₂ ices can lift large particles containing water ice off the surface, whereas sublimation of water ice can only lift smaller, ice-free particles (Gundlach et al. 2020). In addition, other mechanisms, such as the crystallization of amorphous ices and impacts, may be capable of ejecting smaller icy particles. The ice-containing particles ejected from the surface may behave like mini-comets whose trajectories are affected by the recoil force due to anisotropic outgassing of water vapour. Indeed, particles whose trajectories cannot solely be explained by gravity, gas drag, and radiation pressure have been detected for the comets 103P/Hartley 2 (Kelley et al. 2013, 2015) and 67P/Churyumov-Gerasimenko (Agarwal et al. 2016). Furthermore, the sublimation of ices can build up pressure in the interior of a particle, which may lead to their disintegration while they are travelling from the inner to outer coma.

Understanding the outgassing of dust particles requires rigorous multiphysical modelling. Typical thermophysical models applied to the cometary nuclei, however, introduce various simplifications which make them inapplicable to small dust particles. More specifically, the models are often one dimensional, they assume wavelength-independent surface absorption and emission of electromagnetic radiation, the radiative heat transfer is treated as a one-dimensional transport by defining an effective heat transfer coefficient, and the gas flow is treated as the Knudsen flow (Huebner et al. 2006; Prialnik 2004; Gundlach et al. 2020). For small particles, these approximations may introduce significant errors. The thermal evolution of small cometary dust particles is often modelled by assuming isothermal spherical particles and using the analytical Lorenz-Mie solution for the absorbed and emitted radiation and assuming that the sublimation only occurs on the surface of the particle (Lichtenegger & Kömle 1991; Gicquel et al. 2012, 2016). This is a valid assumption for small compact particles, but inside porous particles larger than tens of micrometres, temperature and pressure gradients may occur in a typical coma environment. Thus, a more realistic numerical model is needed to bridge the gap between small and very large particles.

We recently introduced a numerical method for analysing the thermal properties of ice-free particles by treating the radiative heat transfer with the radiative transfer with the reciprocal transactions framework (Markkanen & Agarwal 2019). Here, we extend the method for ice-containing particles by introducing an additional source term to the energy equation that accounts for the latent heat of sublimation and condensation. Because this term depends on the temperature and pressure, we solve the energy balance equation simultaneously with the mass balance equation. We apply the finite-element method (FEM) for the energy balance equation and the direct simulation Monte Carlo (DSMC) method with the dusty gas model for the mass transport.

The article is organized as follows. In Sect. 2, we present the governing equations describing the physics of the problem. The numerical methods we used are explained in Sect. 3. Section 4 introduces the particle model. In Sect. 5, we present the validation of the developed method against the continuous gas flow model and apply the method to study thermal properties of icy cometary dust particles at 1.35 AU. Finally, we present our conclusions in Sect. 6.

2 Governing equations

Recently, we introduced a thermal model for the aggregated ice-free particles composed of submicrometre-sized dust grains in Markkanen & Agarwal (2019). Here, we will extend the method for dust particles containing ice inclusions. We assume that the dust and icegrains are intimately mixed in the microscopic scale. We also assume local equilibrium, in the case of which the solid and gas phases have the same temperature locally and the energy exchange between the solid and gas phases is neglected.

2.1 Energy balance

The energy balance equation for the solid phase is written as $$c_{\text{p}}{\rho }_{\text{s}}\frac{\partial T}{\partial t}-\nabla \cdot \kappa \nabla T=Q_{\text{r}}+Q_{\text{l}},$$(1)

where c_p is the specific heat capacity, ρ_s is the density, T is the temperature, and κ is the conductive heat transfer coefficient. The right-hand side includes the volumetric source terms, the absorbed solar and thermally emitted and reabsorbed radiation Q_r, and the latent heat of sublimation and condensation of ices Q_l. The above values are macroscopic, that is, they are averaged over a small volume element larger than the microstructure of dust and ice grains. We also assume that the density and specific heat capacity of gas are many orders of magnitude lower than those of the solid phase and ignore the heat transport in the gas phase.

On the particle surface, we impose the Neumann boundary condition for the conductive flux as n ⋅ κ∇T = 0, where n is the outer unit normal vector. Thus, only the radiative heat flow can cross the boundary as the convective flow is neglected.

Energy related to the sublimation and condensation processes can be written as $$Q_{\text{l}}=qH,$$(2)

where q is the sublimation rate and H is the latent heat. By defining the saturation vapour pressure as $$P_{\text{sat}}=A{\text{e}}^{-B/T}$$(3)

and using the experimentally obtained coefficients, A and B (Fanale & Salvail 1984), the heat of sublimation is obtained via the Clausius-Clapeyron relation and given by $$H=B\frac{R_{\text{g}}}{\mu },$$(4)

in which R_g is the universal gas constant and μ is the molar mass of the gas.

The production rate is obtained via the Hertz–Knudsen formula as $$q=S(P_{\text{sat}}-P)\sqrt{\frac{\mu }{2\pi R_{\text{g}}T}},$$(5)

where S is the surface-to-volume ratio in the microscale, and P is the pressure in the micropores.

2.2 Massbalance

The mass balance equation for the gas and solid phases is written as $$\begin{array}{l}\frac{\partial {\rho }_{\text{g}}}{\partial t}+\nabla \cdot J=q,\\ \frac{\partial {\rho }_{\text{s}}}{\partial t}=-q,\end{array}$$

where ρ_g and ρ_s are the gas and solid phase densities, respectively, and J is the gas flux density. The macroscopic thermodynamic properties of gas such as the density, velocity, pressure, and flux emerge from the stochastic microscopic description of the state of gas given by the Boltzmann equation: $$\frac{\partial f}{\partial t}+v\cdot {\nabla }_{r}f=C(f,f)+f_{\text{q}},$$(8)

where f(t, r, v) is the time (t) position (r) and velocity (v) dependent distribution density function, f_q is the source function, and C(f, f) is the collision operator. Formally, the collision operator reads $$C(f,f)={\displaystyle \int {\displaystyle \int |}}v-v_{*}|(f_{{*}^{\prime }}{f}^{\prime }-f_{*}f)\sigma (\Omega )\text{d}\Omega \text{d}{v}_{*},$$(9)

where σ is the kernel that describes scattering, the primed functions are the post collision distribution densities, and the subscript * denotes the other particle in the collision pair.

The macroscopic number density n can be computed by integrating the distribution density function over the velocity space as $$n(t,r)={\displaystyle \int f}(t,r,v)\text{d}v,$$(10)

and the gas flux as $$J(t,r)={\displaystyle \int v}f(t,r,v)\text{d}v.$$(11)

We apply the DSMC method to solve Eq. (8), as described in Sect. 3.

3 Numerical solution

We discretize the domain of interest Ω with tetrahedral elements. Then we employ the FEM for the spatial dimension of the energy balance Eq. (1) with the nodal testing wⁿ and basis functions u^m. For the temporal dimension we will use the central finite difference formula $$\frac{\partial T}{\partial t}{|}_{t+\frac{1}{2}}\approx \frac{T_{t+1}-T_t}{\tau }=L_{t+\frac{1}{2}}$$(12)

and interpolate the other terms at $t+\frac{1}{2}$ as $$L_{t+\frac{1}{2}}\approx \frac{L_{t+1}+L_t}{2}.$$(13)

This results into the Crank–Nickolson scheme for the unknown coefficient vector x_t+1 given by $$x_{t+1}={\left(M+\frac{\tau }{2}S\right)}^{-1}\left(M{x}_t+\frac{\tau }{2}F{x}_t-\frac{\tau }{2}S{x}_t+F{x}_{t+1}\right),$$(14)

where the mass and the stiffness matrices are defined as $$M={\rho }_{\text{s}}{c}_{\text{p}}{\displaystyle {\int }_{\Omega }{w}^n}{u}^m\text{d}V,$$(15) $$S={\displaystyle {\int }_{\Omega }\nabla }{w}^n\cdot \kappa \nabla u^m\text{d}V,$$(16)

and the force vector as $$F={\displaystyle {\int }_{\Omega }{w}^n}(Q_{\text{r}}+Q_{\text{l}})\text{d}V.$$(17)

Equation (14) is strongly non-linear because of the force term. Hence, we use an iterative method with under-relaxation to solve the unknown coefficient vector, x_t+1.

We compute the absorbed solar and thermally emitted radiation, Q_r, by using the radiative transfer with reciprocal transactions framework, as described by Markkanen & Agarwal (2019). Energy related to sublimation and condensation requires evaluation of the gas production rate, q, which depends on the temperature, T, and pressure, P. Thus, we need to solve the mass balance equation simultaneously with the energy balance equation. We apply the DSMC method, and its solution thus satisfies the Boltzmann equation, to solve the gas flow by assuming that the flow is a free molecular flow, that is, the mean free path of molecule-molecule collisions is much larger than the mean free path of molecule-grain collisions. The grains are assumed to be much heavier than the gas molecules and, thus, they are stationary, giving us the so-called dusty gas model. We consider the molecule-grain collisions stochastically, as in Ahmadian et al. (2019).

To compute the source function f_q in the DSMC, we introduce new (or remove the existing) DSMC molecules for each tetrahedral element with the index, n, based on the production rate, qⁿ, at any given time step, Δt_DSMC, with weight, qⁿΔt_DSMC. The velocity of the new DSMC molecule is drawn from the local solid-phase temperature-dependent Maxwell-Boltzmann distribution. Then the trajectories of the DSMC molecules are traced within a time step. To account for the collision operator, C, the molecule-grain scattering distance is drawn from the exponential distribution with the mean free path, l_DSMC, and a new velocity for the gas molecule is drawn the Maxwell-Boltzmann distribution defined by the local temperature. Finally, the macroscopic thermodynamical properties, for instance, density, pressure, and flux are sampled by summing up the DSMC molecules in each tetrahedron and averaging them over the time frame. The process is repeated until the steady state or the energy equation time step is reached.

The characteristic time scales for energy and mass transport can be very different. This means that we need to use different time steps (τ, Δt_DSMC) for the energy and mass transport to get a sufficiently fast numerical method. The DSMC time step, Δt_DSMC, must be small enough such that the DSMC molecules cannot travel across a finite cell within a single time step, whereas τ is limited by the non-linearity and the convergence of the iterative method used to solve the energy balance equation. In practice, for cometary dust applications, the DSMC time step must be many orders of magnitudes smaller than the energy balance equation time step. Often, the gas flow reaches the steady state within the energy balance equation time step. In such a case, we stop the DSMC simulation and extrapolate the total gas production rate and the remaining ice mass at t + τ. The extrapolation does not conserve ice and gas mass exactly when a tetrahedron runs out of ice within the time step τ as the ice mass cannot be negative. To avoid the problem, we run the DSMC again by setting the ice mass to zero for the tetrahedron that has run out of ice and calculate the production rates using the total mass loss at the given tetrahedron as $$q_{t+1}=\frac{m_t}{\tau },$$(18)

where m_t is the ice mass at time t.

4 Particle model

We used a particle microstructure model derived by fitting the phase function of the coma of 67P/Churyumov-Gerasimenko in Markkanen et al. (2018). The heat conduction coefficient was derived by fitting the superheating phase function of the same comet in Markkanen & Agarwal (2019).

The microstructure consists of submicrometre-sized organic monomers and micrometre-sized silicate monomers randomly deposited inside a particle with the porosity, Φ = 0.6. The effective heat capacity, c_p = 750 J kg⁻¹ K⁻¹, and the effective density, ρ_s = 1000 kg m⁻³, are estimated from the porosity Φ and the silicate-to-organic volumetric ratio of 1/5. The size distributions of monomers follow a differential power law with the index of − 3, and the minimum and maximum cutoff limits are a_min = 65 nm and a_max = 125 nm for the organic and a_min = 0.6 μm and a_max = 1.3 μm for the silicate monomers. The refractive index for the organic monomers is approximated to be the same as amorphous carbon and is taken from Jäger et al. (1998). For the silicate monomers, the refractive index is taken from Dorschner et al. (1995). Here, we assume that ice is uniformly deposited in the micropores and the amount of ice is so low that it does not have an effect on the scattering properties. The microstructure model gives rise to the incoherent scattering properties of the volume elements which are used to evaluate the radiative part of the energy equation as described by Markkanen & Agarwal (2019). In addition, the gas mean free path, l_DSMC, and the surface to volume ratio, S, are derived from the microstructure model. The parameters used in the model are presented in Table 1.

The macrostructure is described by the tetrahedral mesh in which each tetrahedron can be assigned to different microstructure, effective material properties c_p, ρ_s, κ, and ice content. To generate the macrostructure model, we applied the hierarchical Voronoi partitioning algorithm with the parameters, N₁ = 100, N₂ = 1000, $P_1^{\text{rm}}=1$, and $P_2^{\text{rm}}=0.4,$ from Markkanen et al. (2015). This creates particles with the macroporosity of 0.6 and the mean unit Voronoi cell size of one tenth of the particle size. Figure 1 shows an example computational mesh generated by the two-level hierarchical Voronoi algorithm. In the remainder of the paper, the size of the particle refers to the radius, r, of the smallest sphere containing the particle if not defined otherwise.

5 Results

5.1 Validation

To validate the DSMC implementation, we compared the results with the continuous gas flow model in porous media. The Fick’s law states that the gas flux J is proportional to the gas density ρ_g gradient as $$J=-D\nabla {\rho }_{\text{g}},$$(19)

where D is the diffusion coefficient. Here, the gas density ρ_g is defined asthe gas mass per total volume including dust grains. Using the continuity equation, we have $$\frac{\partial {\rho }_{\text{g}}}{\partial t}-\nabla \cdot D\nabla {\rho }_{\text{g}}=S(P_{\text{sat}}-P)\sqrt{\frac{\mu }{2\pi R_{\text{g}}T}}.$$(20)

Writing the gas density in terms of pressure inside the micropores and using the ideal gas assumption, $$P=\frac{{\rho }_{\text{g}}{R}_{\text{g}}T}{\mu \Phi }$$(21)

leads to the following equation for the gas density $$\frac{\partial {\rho }_{\text{g}}}{\partial t}-\nabla \cdot D\nabla {\rho }_{\text{g}}=S\left(P_{\text{sat}}-\frac{{\rho }_{\text{g}}{R}_{\text{g}}T}{\mu \Phi }\right)\sqrt{\frac{\mu }{2\pi R_{\text{g}}T}}.$$(22)

The diffusion coefficient for the random walk in three dimensions with the Maxwell–Boltzmann mean velocity v_th is given by $$D=\frac{1}{6}{l}_{\text{DSMC}}{v}_{\text{th}}=\frac{1}{6}{l}_{\text{DSMC}}\sqrt{\frac{8R_{\text{g}}T}{\pi \mu }}.$$(23)

Furthermore, by assuming the half Maxwell-Boltzmann distribution above the particle interface, giving the normal mean velocity of v_n = v_th∕4 (Huebner & Markiewicz 2000), the Neumann boundary condition is given as $$n\cdot J={\rho }_{\text{g}}\sqrt{\frac{R_{\text{g}}T}{2\pi \mu }}.$$(24)

Equation (22), supplemented with the boundary condition (24) can be solved with the FEM analogously to the FEM solution of the energy balance Eq. (14) by changing the coefficients and the right-hand side accordingly and introducing the boundary integral for the Neumann boundary condition (which is zero for the energy equation).

First, we considered the simplest possible case; an isothermal T = 200 K spherical particle with the radius r = 0.1 mm and S = 10 000 m⁻¹ in a steady state. No other energy sources were included. We solved the problem using the DSMC and the FEM. Two different DSMC solutions were computed, namely, DSMC 1 and DSMC 2, the latter having twice the number of the DSMC molecules than the former. Figure 2 shows the pressure sampled at the barycentre of each tetrahedron as a function of the distance from the centre of the spherical particle. The solutions are in good agreement with each other but the DSMC solutions have more noise than the FEM solution. The DSMC noise decreases with the increasing number of the DSMC molecules, as expected.

It is interesting to note that when S is constant, Eq. (22) is reduced to the Helmholtz equation in a steady state with the wavenumber squared given by $$k^2=-\sqrt{\frac{R_{\text{g}}T}{2\pi \mu }}S{\Phi }^{-1}{D}^{-1}.$$(25)

The fundamental solution, that is, the Green’s function for the Helmholtz equation is written as $$G(r,r^{\prime })=\frac{e^{ik|r-r^{\prime }|}}{4\pi |r-r^{\prime }|},$$(26)

where r and r′ are the source and observation points, respectively, and i is the imaginary unit. Thus, with the imaginary k the solution is an exponentially decaying evanescent wave, and the pressure drop near the surface can be very steep when SD⁻¹ is large. This also means that in order to solve Eq. (22) with the standard continuous Galerkin FEM, an extremely fine mesh is needed at the sublimation front, which can make computational time prohibitively long.

Next, we simulated an irregularly-shaped dust particle under solar radiation at 1.35 AU. The initial temperature of the particle was 160 K and the volumetric ice content was 0.05. The surface-to-volume ratio parameter was set to S = 1 m⁻¹ to allow for an efficient FEM solution. The volume equivalent radius of the particle was 0.5 mm and it was spinning 1 rpm with the spin axis pointing perpendicular to the direction of the Sun. The solar illumination was switched on at t = 0 s. The time evolution of the average temperature, gas production rate, maximum pressure averaged over a tetrahedron inside the particle and the ice volume fraction are plotted in Fig. 3 computed by using the DSMC and FEM for the gas transport. We observed an excellent agreement between the two different methods. We note, however, that the boundary condition used in the FEM does not account for the gas flow back to the particle, whereas it is taken into account in the DSMC method. The gas production rate and maximum pressure peak at t ≈ 200 s and then they decrease. This happens because sublimation creates an insulating dry dust layer on the particle’s surface that dampens the energy transport into the interior of the particle.

Fig. 1 Crosscut of an example computational mesh generated by the two-level Voronoi partitioning. The yellow tetrahedra correspond to the particle with a specified microstructure and the blue ones correspond to free space.

Fig. 2 Pressure as a function of distance from the centre of the sphere calculated by the DSMC and FEM. The DSMC 2 solution has twice the number of the DSMC molecules than the DSMC 1 solution to demonstrate convergence.

Fig. 3 Comparisons between thermal modelling results obtained by using the DSMC and FEM solutions for gas transport. The test particle has an irregular shape, with the volume equivalent radius of 0.5 mm. The surface-to-volume ratio was set to S = 1 m−1.

5.2 Application to cometary dust

Next, we studied thermal properties of dust in a cometary coma at 1.35 AU. We assumed that the ejected dust particles have temperature T₀ = 160 K and contain5 vol. of water ice uniformly distributed into the micropores at t = 0. Thus, the dust particle was assumed to originate below the hot surface layer in order to contain ice. Here, we do not account for the ejection mechanism as our goal is to study the thermal evolution of icy dust particles once they have been ejected from the nucleus and exposed to the direct sunlight. In the simulation, we assumed that the particles are rotating around the axis perpendicular to the solar direction with the angular speed of 1 rpm. The results were not averaged over an ensemble of particles due to computational time restrictions.

Figure 4 plots the ice volume fraction as a function of time for 0.1 mm-, 1 mm-, and 10 mm-sized particles. We observed that the particleswith r = 0.1 mm run out of ice in approximately a few tens of seconds after the ejection, r = 1 mm in ten minutes, and r = 10 mm particles can be extrapolated to stop subliming after a few hours. We expect that either they run out of ice or create an insulation surface layer, completely dumping sublimation. The time integration was stopped after 20 days of computing using 24 cores for the particle with r = 10 mm. The smaller particles took less than a week each with 24 CPUs to complete the simulation. The biggest bottleneck in the simulations of large particles is the non-linearity of Eq. (14). A very small time step τ is required to solve the non-linear equation prohibiting long time evolution simulations. Thus, using a more sophisticated non-linear iterative solver or finding an efficient preconditioner would be an interesting topic for future research. Also, the DSMC simulation becomes slower with the increasing particle size as the diffusion time is proportional to r². Thus, for large particles, a continuous gas transport model would be preferred, assuming that the stability problem appearing for high S due to the exponentially decaying pressure profiles can be solved.

The sublimation of ices increases pressure inside the particles as demonstrated in Fig. 5. If the pressure is high enough, the particle may disintegrate. An experimentally verified model (Skorov & Blum 2012; Blum et al. 2014; Brisset et al. 2016) for the effective tensile strength of a particle made of aggregates is given by $T_{\text{eff}}\approx {(r_{\text{a}}/1\text{mm})}^{-2/3}$ Pa where r_a is the radius of aggregates. In our particle model, these aggregates can be considered to correspond the Voronoi cells with r_a = r∕10. The tensile strengths for the particles made of aggregates with r_a = 0.1 mm and 1 mm are also presented in Fig. 5 as dashed lines. The maximum pressure for particles with r < 1 mm is less than 3 Pa which is smaller than the corresponding tensile strength of the experimental particles. Thus, it is unlikely that a particle with r < 1 mm will disintegrate into smaller pieces because of pressure. For larger particles, the pressure can be higher. This happens as the sublimation creates a dry hot layer on the particle’s surface which helps to build up pressure at the sublimation front. As shown in Fig. 5, for the particle with r = 10 mm, the maximum pressure can reach the tensile strength of a particle made of aggregates with r_a = 0.1 mm. This indicates that the particleswith r > 10 mm may start fragmenting into smaller pieces but we expect the disintegration to stop at mm-sized fragments. However, a detailed investigation of the fragmentation of particles is beyond the scope of this study and would require a detailed structural mechanical analysis. Albeit, we note that if such fragmentation would occur, it would have a significant effect on the sublimation rate as the insulating layer would be removed from time to time, allowing heat to directly access to the icy part.

Since the particles hold temperature gradients, outgassing of water vapour is anisotropic. This introduces a rocket force that accelerates the particles. The acceleration due to the rocket force is plotted in Fig. 6. The acceleration decreases with increasing particle size. This is clear as the sublimation rate is proportional to the illuminated area and the mass is proportionalto the volume, the acceleration roughly scales linearly with the inverse of size. For small particles, the acceleration is high but it only occurs for a few seconds, whereas for large particles, the acceleration is small but it acts for a considerabletime, which may have an effect on the particle’s trajectory. The direction of the acceleration depends on the rotationstate of the particle. For a slowly rotating particle, the force points towards the antisolar direction whereas for a fast-rotating particle the direction is shifted towards the direction perpendicular to the rotation axis, e_rot, and to the axis towards the Sun, e_sun, as presentedin Fig. 7.

Finally, we studied how the surface-to-volume ratio, S, in Eq. (5) affects the results. Since sublimation requires that P_sat − P > 0 and it depends on the surface-to-volume ratio of the ice surface, S_ice, whereas condensation requires that P_sat − P < 0 and it depends on the total surface-to-volume ratio S_tot = S_ice + S_dust, we define $$S=\{\begin{array}{c}{S}_{\text{ice}}\text{if}{P}_{\text{sat}}-P>0\\ S_{\text{tot}}\text{if}{P}_{\text{sat}}-P<0.\end{array}$$(27)

If ice fully covers the dust grains, it is clear that S_ice = S_tot. If the particle is an aggregate of ice and dust grains, then S_ice < S_tot. For equisized icegrains: $$S_{\text{ice}}=3v_{\text{ice}}/r_{\text{ice}}$$(28)

where v_ice is the volumetric filling factor of ice grains of radius r_ice. Thus, S depends on how ice is distributed inside the particle. We assumed that v_ice = 0.05 and compared four different ice distributions. We also assumed that the ice distribution does not change the scattering, absorption and emission properties. In the first case, ice covered all the dust grains giving S_ice = 1 × 10⁷ m⁻¹. In the second, the ice grains of radius r_ice = 0.2 μm were evenly distributed inside the particle giving S_ice = 7.5 × 10⁵ m⁻¹. In this case, the ice grains were approximately the same size as the dust grains. In the third, the ice grains were bigger, r_ice = 2 μm, than the dust grains and S_ice = 7.5 × 10⁴ m⁻¹. In the last case, the ice grains were much bigger than dust grains and r_ice = 20 μm.

Comparisons for the total ice volume fraction and the maximum pressure as a function of time are presented in Fig. 8. We observed that when S_ice is large enough, the system is saturated and increasing S_ice does not affect the results. This is because gas diffusion gives the upper limit to the total sublimation rate. When S_ice is small enough, the total sublimation rate is smaller and the system is no longer saturated and limited by gas diffusion but is limited by the ice surface area. In such a case, the received energy is used more to increase the temperature of the particle which in turn allows higher local pressure. In fact, as seen from Eq. (22), the thickness of the sublimation front depends on the length scale, h, and the diffusion coefficient, D, and is proportional to the factor, hSD⁻¹. When the factor is large, the sublimation front is thin compared to the particle’s size and the total sublimation rate depends on the macroscopic icy surface area rather than the microscopic surface-to-volume ratio, S. It is, thus, difficult to retrieve information on how the ice grains are mixed with the dust grains in the microscopic level from the total gas production rate.

Fig. 4 Ice volume fraction as a function of time for different particle sizes.

Fig. 5 Maximum pressure as a function of time for different particle sizes.

Fig. 6 Acceleration due to outgassing as a function of time for different particle sizes.

Fig. 7 Acceleration vector components for a particle with r = 1 cm rotating 1 rpm about the axis, erot, perpendicular to the solar direction, esun, due to outgassing as a function of time.

Fig. 8 Ice volume fraction and the maximum pressure as a function of time computed for varying ice surface-to-volume ratios, Sice, corresponding to fully coated dust grains, and ice grains of sizes 0.2 μm, 2 μm, and 20 μm uniformly mixed with dust grains.

6 Conclusions

We presented a novel thermophysical model and its numerical solution for icy cometary dust particles. The model employs the radiative transfer with reciprocal transactions for the radiative heat transport and the Boltzmann equation together with the Hertz-Knudsen formula for gas transport solved by the DSMC methods. Energy changes related to absorption, thermal emission, and phase changes were incorporated as source terms into the energy balance equation with the Fourier’s heat conduction and solved by the FEM.

The developed method allows for thermal analysis of up to cm-sized particles composed of submicrometre-sized grains mixed with ice with a moderate computing power. The method can find applications in understanding the thermal physics of cometary dust particles and explaining experimental laboratory measurements, and help us to develop more accurate approximate thermal modelling methods.

We also showed that water ice sublimation inside large dust particles may generate enough pressure to reach the tensile strength and, thus, to possibly disintegrate the particles when comets are close to the Sun. Also, outgassing of water vapour can play a crucial role in the dynamics of icy dust particles after having been lifted off the comet’s nucleus. Both of these processes should leave observable effects on the remote observables via changes in the size and spatial distributions.

Finally, we showed that the thermal evolution of large dust particles is quite insensitive to the microscopic mixing scale of dust and ice, which effectively makes modelling easier. On the other hand, getting information on the mixing scale and ratio of dust and icefrom the total gas production rate is not trivial.

Acknowledgements

This work has been funded by the ERC Starting Grant No. 757390 Comet and Asteroid Re-Shaping through Activity (CAstRA). Computational resources have been provided by Gesellschaft für Wissenschaftliche Datenverarbeitung mbH Göttingen (GWDG).

References

All Tables

All Figures

	Fig. 1 Crosscut of an example computational mesh generated by the two-level Voronoi partitioning. The yellow tetrahedra correspond to the particle with a specified microstructure and the blue ones correspond to free space.
In the text

	Fig. 2 Pressure as a function of distance from the centre of the sphere calculated by the DSMC and FEM. The DSMC 2 solution has twice the number of the DSMC molecules than the DSMC 1 solution to demonstrate convergence.
In the text

	Fig. 3 Comparisons between thermal modelling results obtained by using the DSMC and FEM solutions for gas transport. The test particle has an irregular shape, with the volume equivalent radius of 0.5 mm. The surface-to-volume ratio was set to S = 1 m−1.
In the text

	Fig. 4 Ice volume fraction as a function of time for different particle sizes.
In the text

	Fig. 5 Maximum pressure as a function of time for different particle sizes.
In the text

	Fig. 6 Acceleration due to outgassing as a function of time for different particle sizes.
In the text

	Fig. 7 Acceleration vector components for a particle with r = 1 cm rotating 1 rpm about the axis, erot, perpendicular to the solar direction, esun, due to outgassing as a function of time.
In the text

	Fig. 8 Ice volume fraction and the maximum pressure as a function of time computed for varying ice surface-to-volume ratios, Sice, corresponding to fully coated dust grains, and ice grains of sizes 0.2 μm, 2 μm, and 20 μm uniformly mixed with dust grains.
In the text