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{\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_{\mathrm{n}}\) and the shear stress \(\tau_{\mathrm{n}}\) 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}^{(\boldsymbol{n})}\). 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})} = (\boldsymbol{n}^\mathsf{T} \cdot \boldsymbol{\sigma})^\mathsf{T} \end{align}

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

\begin{align} \sigma_{\mathrm{n}} = \boldsymbol{T}^{(\boldsymbol{n})} \cdot \boldsymbol{n}= (\boldsymbol{n}^\mathsf{T} \cdot \boldsymbol{\sigma})^\mathsf{T} \cdot \boldsymbol{n} \end{align}

In component-wise notation this results in:

\begin{align} \sigma_{\mathrm{n}} = \big(\boldsymbol{T}^{(\boldsymbol{n})}\big)_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}_{\mathrm{n}}\), which results from the subtraction of the normal stress vector \(\boldsymbol{\sigma}_{\mathrm{n}}=\sigma_{\mathrm{n}}\boldsymbol{n}\) from the traction vector: \begin{align} \boldsymbol{\tau}_{\mathrm{n}} = \boldsymbol{T}^{(\boldsymbol{n})} - \sigma_{\mathrm{n}}\boldsymbol{n} \end{align}

In component-wise notation this results in:

\begin{align} \tau_{\mathrm{n},j} = \big(\boldsymbol{T}^{(\boldsymbol{n})}\big)_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_{\mathrm{n}}^2 = |\boldsymbol{\tau}_{\mathrm{n}}|^2 = \tau_{\mathrm{n},j} \, \tau_{\mathrm{n},j} \\[6pt] \Leftrightarrow\;& \tau_{\mathrm{n}}^2 = (\sigma_{ij} \, n_i - \sigma_{\mathrm{n}} \, n_j) (\sigma_{mj} \, n_m - \sigma_{\mathrm{n}} \, n_j) \\[6pt] \Leftrightarrow\;& \tau_{\mathrm{n}}^2 = \sigma_{ij} \, \sigma_{mj} \, n_i \, n_m - \sigma_{ij} \, \sigma_{\mathrm{n}} \, n_i \, n_j - \sigma_{mj} \, \sigma_{\mathrm{n}} \, n_m \, n_j + \sigma_{\mathrm{n}}^2 \\[6pt] \Leftrightarrow\;& \tau_{\mathrm{n}}^2 = (\sigma_{ij} \, n_i) (\sigma_{mj} \, n_m) - (\sigma_{ij} \, n_i \, n_j) \, \sigma_{\mathrm{n}} - (\sigma_{mj} \, n_m \, n_j) \, \sigma_{\mathrm{n}} + \sigma_{\mathrm{n}}^2 \\[6pt] \Leftrightarrow\;& \tau_{\mathrm{n}}^2 = \big(\boldsymbol{T}^{(\boldsymbol{n})}\big)_j \, \big(\boldsymbol{T}^{(\boldsymbol{n})}\big)_j - \sigma_{\mathrm{n}}^2 \\[6pt] \Leftrightarrow\;& \tau_{\mathrm{n}}^2 = \big(\boldsymbol{T}^{(\boldsymbol{n})}\big)^2 - \boldsymbol{\sigma}_{\mathrm{n}}^2 \\[6pt] \Rightarrow\;& \tau_{\mathrm{n}} = \sqrt{ \big(\boldsymbol{T}^{(\boldsymbol{n})}\big)^2 - \boldsymbol{\sigma}_{\mathrm{n}}^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)} = (\boldsymbol{e}_r^\mathsf{T} \cdot \boldsymbol{\sigma})^\mathsf{T} = \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_{\mathrm{n}} = \boldsymbol{T}^{(\boldsymbol{e}_r)} \cdot \boldsymbol{e}_r = \sigma_{rr} \\[6pt] &\tau_{\mathrm{n}} = \sqrt{ \big(\boldsymbol{T}^{(\boldsymbol{e}_r)}\big)^2 - \boldsymbol{\sigma}_{\mathrm{n}}^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}\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}\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}\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 + \bigg(\zeta + \frac{4}{3}\mu\bigg) \frac{\partial u_r}{\partial r}\Big|_{r=R} \end{align}

To be continued ...