Planetesimal-Erosion

Flow around a planetesimal

This project aims to investigate the flow around and erosion of a planetisimal by a viscous fluid under physical conditions that typically occur in protoplanetary discs. As a starting point, the erosion model presented in Cedenblad et al. (2021) will be numerically reproduced and systematically extended in perspective.

A first step is to investigate the stresses acting on the surface of a stationary sphere flown around by a gas. In particular, the influence of the Mach number in the range \(0<\mathcal{M}<2\) on the shear and normal stress distribution is to be investigated in order to gain a basic understanding of the flow conditions relevant for erosion.

Theoretical Framework

Required is the distribution of normal and shear stresses acting on the surface of a sphere at rest with radius \(R\) as a function of the Mach number \(\mathcal{M}\). The Mach number is defined as the ratio between the velocity \(\boldsymbol{u}\) of the flowing fluid (or an object moving through the fluid, depending on the chosen reference system) and the speed of sound \(c_\mathrm{S}\) of the considered medium:

\begin{align} \mathcal{M} = \frac{|\boldsymbol{u}|}{c_\mathrm{S}} \end{align}

The normal stress describes the component of the mechanical stress acting perpendicular to the surface and contains in particular the contribution of the static pressure. Shear stress, on the other hand, acts tangentially to the surface and is closely linked to velocity gradients in the fluid. Both quantities are components of the stress tensor \(\boldsymbol{\sigma}\) and play a central role in the description of force transmissions between fluid and body, especially with regard to possible erosion processes on the object surface.

Normal and shear stress

In fluid mechanics, the mechanical stress acting on a surface element is described by the stress tensor \(\boldsymbol{\hat\sigma}\). It specifies the internal stresses acting at each point of a continuum due to external forces. The stress tensor is a second-rank tensor whose components \(\sigma_{ij}\) represent the stresses acting on a surface with a normal in the direction of the \(j\)-th coordinate axis due to forces in the direction of the \(i\)-th coordinate axis. The stress tensor describes both the normal stress \(\sigma\) and the shear stress \(\tau\) acting on a surface element with unit normal vector \(\boldsymbol{n}\).

The mechanical stress that a fluid exerts on a certain surface within the continuum is described by the so-called traction vector \(\boldsymbol{T}\). This specifies the force per unit area that acts on a surface with a unit normal vector \(\boldsymbol{n}\). Mathematically, the traction vector results from the transpose of the application of the transposed normal vector to the stress tensor:

\begin{align} \boldsymbol{T} = \boldsymbol{n} \cdot \boldsymbol{\hat\sigma} \end{align}

The normal stress is the projection of this traction vector onto the normal coordinate axis:

\begin{align} \sigma = \boldsymbol{T} \cdot \boldsymbol{n}= (\boldsymbol{n} \cdot \boldsymbol{\hat\sigma}) \cdot \boldsymbol{n} \end{align}

In component-wise notation this results in:

\begin{align} \sigma = T_j \, n_j = \sigma_{ij} \, n_i \, n_j \end{align}

The shear stress corresponds to the magnitude of the shear stress vector \(\boldsymbol{\tau}\), which results from the subtraction of the normal stress vector \(\boldsymbol{\sigma}=\sigma\boldsymbol{n}\) from the traction vector: \begin{align} \boldsymbol{\tau} = \boldsymbol{T} - \sigma\boldsymbol{n} \end{align}

In component-wise notation this results in:

\begin{align} \tau_{\mathrm{n},j} = T_j - \sigma_{\mathrm{n}} \, n_j = \sigma_{ij} \, n_i - (\sigma_{kl} \, n_k \, n_l) \, n_j \end{align}

The shear stress is then given by :

\begin{align} \phantom{\Rightarrow}\;& \tau^2 = |\boldsymbol{\tau}|^2 = \tau_j \, \tau_j \\[6pt] \Leftrightarrow\;& \tau^2 = (\sigma_{ij} \, n_i - \sigma \, n_j) (\sigma_{mj} \, n_m - \sigma \, n_j) \\[6pt] \Leftrightarrow\;& \tau^2 = \sigma_{ij} \, \sigma_{mj} \, n_i \, n_m - \sigma_{ij} \, \sigma \, n_i \, n_j - \sigma_{mj} \, \sigma \, n_m \, n_j + \sigma^2 \\[6pt] \Leftrightarrow\;& \tau^2 = (\sigma_{ij} \, n_i) (\sigma_{mj} \, n_m) - (\sigma_{ij} \, n_i \, n_j) \, \sigma - (\sigma_{mj} \, n_m \, n_j) \, \sigma + \sigma^2 \\[6pt] \Leftrightarrow\;& \tau^2 = T_j \, T_j - \sigma^2 \\[6pt] \Leftrightarrow\;& \tau^2 = \boldsymbol{T}^2 - \boldsymbol{\sigma}^2 \\[6pt] \Rightarrow\;& \tau = \sqrt{ \boldsymbol{T}^2 - \boldsymbol{\sigma}^2 } \end{align}

Normal and shear stresses on a sphere in a subsonic to slightly supersonic flow

The aim is now to determine the normal and shear stresses on the surface of the sphere. Since it is a sphere with centre at the origin, the unit normal vector at each point of the surface is identical to the unit vector in radial direction \(\boldsymbol{e}_r\):

\begin{align} \boldsymbol{n} = \boldsymbol{e}_r \end{align}

The traction vector on the surface of the sphere is therefore as follows:

\begin{align} \boldsymbol{T} = \boldsymbol{e}_r \cdot \boldsymbol{\hat\sigma} = \sigma_{rr} \boldsymbol{e}_r + \sigma_{r\theta} \boldsymbol{e}_\theta + \sigma_{r\varphi} \boldsymbol{e}_\varphi \end{align}

Die Normal- bzw.~ Scherspannung sind dann gegeben durch:

\begin{align} &\sigma = \boldsymbol{T} \cdot \boldsymbol{e}_r = \sigma_{rr} \\[6pt] &\tau_{\mathrm{n}} = \sqrt{ \boldsymbol{T}^2 - \boldsymbol{\sigma}^2 } = \sqrt{ \sigma_{r\theta}^2 + \sigma_{r\varphi}^2 } \end{align}

To explicitly determine the normal and shear stresses acting on the spherical surface, the stress tensor of the surrounding fluid is now required. This is generally given in a fluid by:

\begin{align} \boldsymbol{\sigma} = - \bigg[p - \zeta(\boldsymbol{\nabla} \cdot \boldsymbol{u})\bigg] \mathbf{I} + \mu \bigg[\boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^\mathsf{T} - \frac{2}{3} (\boldsymbol{\nabla} \cdot \boldsymbol{u})\mathbf{I}\bigg] \end{align}

Where \(p\) is the pressure field, \(\zeta\) is the volume viscosity, \(\boldsymbol{u}\) is the velocity field of the fluid and \(\mathbf{I}\) is the unit tensor.

The \(rr\)-component of the stress tensor is then given by: \begin{align} \phantom{\Rightarrow}\;& \sigma_{rr} = - \bigg[p - \zeta (\boldsymbol{\nabla} \cdot \boldsymbol{u})\bigg] + \mu \bigg[2 (\boldsymbol{\nabla} \boldsymbol{u})_{rr} - \frac{2}{3} (\boldsymbol{\nabla} \cdot \boldsymbol{u})\bigg] \\[6pt] \Leftrightarrow\;& \sigma_{rr} = - p + 2 \mu (\boldsymbol{\nabla} \boldsymbol{u})_{rr} + \bigg(\zeta - \frac{2}{3}\mu\bigg) (\boldsymbol{\nabla} \cdot \boldsymbol{u}) \end{align}

\((\boldsymbol{\nabla} \boldsymbol{u})_{rr}\) is here the \(rr\)-component of the vector gradient \(\boldsymbol{\nabla} \boldsymbol{u}\), which is given by:

\begin{align} (\boldsymbol{\nabla} \boldsymbol{u})_{rr}=\frac{\partial u_r}{\partial r} \end{align}

The divergence of the velocity field generally has the following form in spherical coordinates:

\begin{align} \boldsymbol{\nabla} \cdot \boldsymbol{u} = \frac{1}{r^2} \frac{\partial r^2 u_r}{\partial r} +\frac{1}{r \sin(\theta)} \frac{\partial \sin(\theta) u_{\theta}}{\partial \theta} +\frac{1}{r \sin(\theta)} \frac{\partial u_{\varphi}}{\partial \varphi} \end{align}

Due to the fact that the fluid flows in homogeneously and flows around a sphere, the velocity field around the sphere has cylindrical symmetry. Therefore, the azimuthal derivative is zero when considering the velocity field around the sphere. In addition, it follows from the given initial condition that the fluid flows in along the \(z\)-axis, together with the spherical symmetry of the sphere, that the velocity field and all other relevant vector fields describing this problem cannot have an azimuthal component. This implies for the divergence of the velocity field:

\begin{align} \phantom{\Rightarrow}\;& \boldsymbol{\nabla} \cdot \boldsymbol{u} = \frac{1}{r^2} \frac{\partial r^2 u_r}{\partial r} + \frac{1}{r \sin(\theta)} \frac{\partial \sin(\theta) u_{\theta}}{\partial \theta} \\[6pt] \Leftrightarrow\;& \boldsymbol{\nabla} \cdot \boldsymbol{u} = \frac{1}{r^2} \bigg(2 r u_r + r^2 \frac{\partial u_r}{\partial r}\bigg) + \frac{1}{r \sin(\theta)} \bigg( \cos(\theta) u_\theta + \sin(\theta) \frac{\partial u_{\theta}}{\partial \theta} \bigg) \end{align}

Inserting yields for the \(rr\)-component of the stress tensor:

\begin{align} \sigma_{rr} = - p + 2 \mu \frac{\partial u_r}{\partial r} + \bigg(\zeta - \frac{2}{3}\mu\bigg) \bigg(\frac{1}{r^2} \bigg(2 r u_r + r^2 \frac{\partial u_r}{\partial r}\bigg) + \frac{1}{r \sin(\theta)} \bigg( \cos(\theta) u_\theta + \sin(\theta) \frac{\partial u_{\theta}}{\partial \theta} \bigg)\bigg) \end{align}

Since the wanted stresses act on the surface of the sphere, the stress tensor must be evaluated at exactly this position. For the \(rr\)-component of the stress tensor on the surface of the sphere, this results in:

\begin{align} \sigma_{rr}\Big|_{r=R}= - p(r=R) + 2 \mu \frac{\partial u_r}{\partial r}\Big|_{r=R} + \bigg(\zeta - \frac{2}{3}\mu\bigg) \bigg(\frac{1}{R^2} \bigg(2 R u_r\Big|_{r=R} + R^2 \frac{\partial u_r}{\partial r}\Big|_{r=R}\bigg) + \frac{1}{R \sin(\theta)} \bigg( \cos(\theta) u_\theta\Big|_{r=R} + \sin(\theta) \frac{\partial u_{\theta}}{\partial \theta}\Big|_{r=R} \bigg)\bigg) \end{align}

It follows from the viscosity of the fluid under consideration that the relative velocity between the fluid and the object under consideration must be zero at the boundary surfaces. This is a result of the internal friction in viscous fluids, as the shear forces do not allow any relative movement at the boundary surfaces. Due to the fact that the sphere considered here is at rest, the following must therefore apply at the surface of the sphere:

\begin{align} \boldsymbol{u}(r=R,\theta)=\boldsymbol{0}\quad \forall\;\theta\in[0,\pi] \end{align}

This results in the following for the components of the velocity field:

\begin{align} u_r\Big|_{r=R}=u_{\theta}\Big|_{r=R}=0 \end{align}

Furthermore, the disappearance of the velocity vectors on the spherical surface and its spherical symmetry implies that the polar derivatives on the surface also disappear:

\begin{align} \frac{\partial u_r}{\partial\theta}\Big|_{r=R}= \frac{\partial u_{\theta}}{\partial\theta}\Big|_{r=R}=0 \end{align}

This means that the \(rr\)-component of the stress tensor, evaluated at the surface of the sphere, is given as follows:

\begin{align} \sigma_{rr}\Big|_{r=R} = - p\Big|_{r=R} + \bigg(\zeta + \frac{4}{3}\mu\bigg) \frac{\partial u_r}{\partial r}\Big|_{r=R} \end{align}

The \(r \theta\)-component of the stress tensor is given by:

\begin{align} \phantom{\Rightarrow}\;& \sigma_{r \theta} = \mu\bigg[(\boldsymbol{\nabla}\boldsymbol{u})_{r \theta} + (\boldsymbol{\nabla}\boldsymbol{u})_{\theta r}\bigg] \\[6pt] \Leftrightarrow\;& \sigma_{r \theta} = \mu\bigg[\frac{1}{r}\frac{\partial u_r}{\partial\theta} - \frac{u_\theta}{r} + \frac{\partial u_\theta}{\partial r}\bigg] \end{align}

On the surface of the sphere, the corresponding component of the stress tensor is obtained:

\begin{align} \sigma_{r \theta}\Big|_{r=R} = \mu\frac{\partial u_\theta}{\partial r}\Big|_{r=R} \end{align}

The \(r \varphi\)-component of the stress tensor is given by:

\begin{align} \phantom{\Rightarrow}\;& \sigma_{r \varphi} = \mu\bigg[(\boldsymbol{\nabla}\boldsymbol{u})_{r \varphi} + (\boldsymbol{\nabla}\boldsymbol{u})_{\varphi r}\bigg] \\[6pt] \Leftrightarrow\;& \sigma_{r \varphi} = \mu\bigg[\frac{1}{r\sin(\theta)}\frac{\partial u_r}{\partial\varphi} - \frac{u_\varphi}{r} + \frac{\partial u_\varphi}{\partial r}\bigg] \end{align}

Taking into account the symmetry properties of the system, the corresponding component of the stress tensor is calculated on the surface of the sphere:

\begin{align} \sigma_{r \varphi}\Big|_{r=R} = 0 \end{align}

Therefore, the normal and shear stresses acting on the surface of the sphere are given by:

\begin{align} &\sigma\Big|_{r=R} = \sigma_{rr}\Big|_{r=R} = - p\Big|_{r=R} + \bigg(\zeta + \frac{4}{3}\mu\bigg) \frac{\partial u_r}{\partial r}\Big|_{r=R} \\[6pt] &\tau\Big|_{r=R} = \sqrt{\sigma_{r\theta}\Big|_{r=R}^2 + \sigma_{r\varphi}\Big|_{r=R}^2} = \mu\bigg|\frac{\partial u_\theta}{\partial r}\Big|_{r=R}\bigg| \end{align}

Normal and shear stresses on a sphere in a protoplanetary disc

Due to the fact that the gases in protoplanetary discs consist to \(99 \%\) of hydrogen and helium (see Owen et al. 2020) and their volume viscosities are negligibly small or zero (see Landau & Lifshitz, Fluid Mechanics, 2nd ed., 1987), it is assumed in the following that the volume viscosity in protoplanetary discs can be considered to be zero:

\begin{align} \zeta = 0 \end{align}

Under these assumptions, the normal and shear stresses acting on a spherical surface in a protoplanetary disc are given by:

\begin{align} &\sigma\Big|_{r=R} = - p\Big|_{r=R} + \frac{4}{3}\mu \frac{\partial u_r}{\partial r}\Big|_{r=R} \\[6pt] &\tau\Big|_{r=R} = \mu\bigg|\frac{\partial u_\theta}{\partial r}\Big|_{r=R}\bigg| \end{align}

Natural units

For the given problem, it makes sense to measure the relevant quantities in units of the object radius \(R\), the velocity at infinity \(u_{\infty}\) and the density at infinity \(\rho_{\infty}\). The corresponding base units, i.e. the units of mass, length and time, are then given by:

\begin{align} \rho_{\infty}R^3 \, , \; R \, , \; \frac{R}{u_{\infty}} \end{align}

It follows then, that the dynamic viscosity \(\mu_{\infty}^*\) at infinity, in the chosen natural units, is given only by the Reynolds number \(\text{Re}_{\infty}\) at infinity:

\begin{align} \phantom{\Rightarrow}\;& \mu_{\infty}^* = \frac{\mu_{\infty}}{\rho_{\infty} u_{\infty} R} \; , \quad \mu_{\infty} = \frac{\rho_{\infty} u_{\infty} R}{\text{Re}_{\infty}} \\[6pt] \Leftrightarrow\;& \mu_{\infty}^* = \frac{1}{\text{Re}_{\infty}} \end{align}

Since the dynamic viscosity is set constantly to the value \(\mu_{\infty}\) in the entire simulation area, the following relations result for the normal and shear stresses in the chosen natural units:

\begin{align} &\sigma^*\Big|_{r^*=1} = - p^*\Big|_{r^*=1} + \frac{4}{3\text{Re}_{\infty}} \frac{\partial u_r^*}{\partial r^*}\Big|_{r^*=1} \\[6pt] &\tau^*\Big|_{r^*=1} = \frac{1}{\text{Re}_{\infty}}\bigg|\frac{\partial u_\theta^*}{\partial r^*}\Big|_{r^*=1}\bigg| \end{align}

In addition, the speed of sound \(c_{S,\infty}^*\) at infinity, measured in the chosen natural units, is given by:

\begin{align} c_{S,\infty}^* = \frac{c_{S,\infty}}{u_{\infty}} = \frac{1}{\mathcal{M}_\infty} \end{align}

Numerical methods

The simulations carried out as part of this project are performed using the PLUTO simulation software (Mignone et al., 2007). PLUTO is an open source software developed by Mignone et al. (2007). It is used for the numerical simulation of real fluids, i.e. viscous, compressible and thermally conductive fluids. The simulations are carried out by numerical integration of the differential equations underlying the hydrodynamic processes. In this project, a third-order Runge-Kutta method is used for time integration and a third-order finite volume method for spatial discretisation.

Computational domain

As explained a sphere with the radius \(R\) around which a viscous fluid flows with an inflow velocity \(\boldsymbol{u}=-u_{\infty}\boldsymbol{e}_z\) is investigated. Due to the spherical symmetry of the problem, the computational domain is modelled in spherical coordinates, with the centre of the sphere at the origin of the coordinate system. However, due to the axial symmetry, it is sufficient to set up a two-dimensional simulation domain, whereby only the radial component \(r\) and the polar component \(\theta\) are considered.

The grid extends from an inner radius \(R_{\text{I}}\), which corresponds to the spherical radius \(R\), to an outer radius \(R_{\text{A}}\). In the polar direction, the grid starts at an angle of \(\theta_\text{S}=0\) and extends to \(\theta_{\text{E}}=\pi\). The number of cells \(N_r\) in the radial direction and \(N_\theta\) in the polar direction is selected so that the cells have an approximately square shape. This configuration is crucial for the numerical precision of the simulations, as highly distorted cells are more susceptible to numerical errors.

To maximise the efficiency of the simulations, a linearly increasing cell size is used in the radial direction. This ensures a high resolution in the area of particular interest around the sphere, while areas further away from the sphere are displayed with a lower resolution. It follows from the PLUTO User's Guide that the length \(\Delta r_i\) of the ith grid cell in the radial direction is given by the following equation:

\begin{align} \Delta r_i=r_{i-\frac{1}{2}} \big(10^{\Delta\xi}-1\big) \end{align}