Simulation of Comet-Solar Wind Interaction in the AMROC Framework with the IGLM-MHD Scheme

Abstract

In this work we present 2-dimensional, 1-fluid magnetohydrodynamic (MHD) simulations of the interaction of a comet with the solar wind. Even though our model can in principle run simulations with time-dependent physical parameters representing the comet travelling along its elliptic orbit, here we focus on “snapshots” at four heliocentric distances, namely: 1.3, 2.0, 2.7 and 3.25 AU. This strategy allows us to analyse some key features of the comet at relatively low computational cost. Our aim here is to present our implementation of such a scenario with the AMROC computational code and study the main features of the comet, which has the physical characteristics of 67P/Churyumov-Gerasimenko (67P/CG). The dependence of the physical parameters on the heliocentric distance is given by means of Parker’s model and of functions obtained from curve fitting of data collected in the literature. In addition, we model comets with three forms: circular, elliptic and bilobed (or “peanut-shaped”.) AMROC is in development and the model shown here is an valuable contribution to these efforts, as it represents a new kind of problem implemented in this code. Furthermore, our model provides the possibility of creation of new and interesting scenarios with regard to comets and other sources of plasma.

Appendices

Appendix A: On the normalized units used in this paper

We use a normalization scheme such that the original CGS values of mass m, time t and length l are converted into new dimensionless ones ($m_{\text {\tiny {N}}}$, $t_{\text {\tiny {N}}}$ and $l_{\text {\tiny {N}}}$) by

$$\begin{aligned} \begin{aligned} m_{\text {\tiny {N}}}&=\frac{m}{10^{-23+3\delta }{g}},\\ t_{\text {\tiny {N}}}&=\frac{t}{10^{(-13+2\delta )/2}{s}},\\ l_{\text {\tiny {N}}}&=\frac{l}{10^{\delta }~\text {cm}}, \end{aligned} \end{aligned}$$

(A1)

where the appropriate value of $\delta$ is defined for each case, bearing in mind the size of the domain.

With the units given by (A1), the density $\rho$, pressure p and the components $B_{j}$ of the magnetic field are normalized in the following way:

$$\begin{aligned} \begin{aligned} \rho _{\text {\tiny {N}}}&=\frac{\rho }{10^{-23}{g\, \text {cm}^{-3}}}, \\ p_{\text {\tiny {N}}}&=\frac{p}{10^{-10}{g\,\text {cm}^{-1}\,s^{-2}}}, \\ B_{j}^{\text {\tiny {N}}}&=\frac{B_{j}}{10^{-5}{G}}. \end{aligned} \end{aligned}$$

(A2)

In addition, the magnetic permeability assumes the value $\mu ^{\text {\tiny {N}}}_{0}=4\pi$.

Appendix B: Accuracy Verification of the SAMR Implementation

Deiterding [21] presents a two-dimensional test created to verify the accuracy and the order of convergence of the SAMR scheme and compare them with the ones for uniform grid simulations. This test consists in a simulation in a computational domain $\Omega =[-1,1]^2$ with periodic boundary conditions at all borders. Regarding the physical conditions, the density $\rho$ has the initial profile given by the Gaussian function

$$\begin{aligned} \rho _{0}(x,y)=1+\textrm{e}^{-(x^2+y^2)/\sigma ^2} , \end{aligned}$$

(B3)

with $\sigma =0.25$. On the other hand, the initial values of velocity and pressure are given respectively by $\textbf{v}=(1,1)$ and $p_0=1$. The simulation runs until the maximum time of $t_\textrm{max}=2$, when the initial conditions are replicated, and this provides a straightforward way to calculate the error norms on uniform grids.

Here, Deiterding uses the $L_1$-error, evaluated for $\rho$ by means of the expression

$$\begin{aligned} L_1(\rho )=L_1(\Delta x_{\,l_\textrm{max}},\Delta y_{\, l_\textrm{max}},G_{l_\textrm{max}})+\sum ^{l_\textrm{max}-1}_{l=0}L_1(\Delta x_{\,l},\Delta y_{\,l},G_{l}\setminus G_{l+1}), \end{aligned}$$

(B4)

where $l_\textrm{max}$ labels the highest level of refinement and $G_{l}$ is the domain of level l (Section 4.) In addition, in Eq. (B4)

$$\begin{aligned} L_{1}(\Delta x, \Delta y, \cdot )=\sum ^{l}_{i,j}|\rho _{ij}-\rho _{0}(x^{i}_{l},y^{j}_{l})|\Delta x_{\,l}\Delta y_{\,l}. \end{aligned}$$

(B5)

That is, $L_1(\Delta x, \Delta y, \cdot )$ is the error norm on a sub-domain of level l. Note that, in the SAMR case, (B4) is calculated on each $G_{l}$ without considering higher levels of refinement. More details can be found in [21].

Figure 13 shows the $L_1$-errors and the respective orders of convergence for the aforementioned simulation for two cases: with uniform grid and with the SAMR implementation. The last runs with two levels of refinement and $r=2$ (see Section 4). Each entry in the x-axis refers to a domain with $N^2$ cells for the uniform grid. In the SAMR case, even if the simulation starts with $N=20$, we have at most $N=80$ after two levels of refinement. Therefore, for the SAMR cases, there is no data for $N=20$ and 40. The order of convergence is defined by $\log _2(L^k_1/L^{k+1}_1)$, where the integer $k\geqslant 1$ denotes the position in the x-axis, with $k=1$ corresponding to $N=20$. In addition, it is analysed the role of the conservative flux correction (“fixup”, see details in Section 2.1.6 of [21]) in the SAMR simulations.

Fig. 13

Upper diagram: $L_{1}$-Error norm for uniform grid, SAMR with fixup and SAMR without fixup; lower panel: Order of convergence for the same three cases

As observed by Deiterding, the absolute errors for the SAMR case are slightly higher than in the uniform grid. However, in all cases the errors decrease as N increases, which indicates convergence. Without fixup, the errors (order) become slightly higher (lower) than in the situation with fixup.

Still according to Deiterding, a second-order accurate method can be inferred from the values of the order of convergence, even in the SAMR cases. Furthermore, note that the uniform grid simulation has the higher order of convergence for the highest resolution.

Appendix C: Elliptic and Bilobed Comets

1.1 C.1 Elliptic Comet

We start with the elliptic coordinates given by

$$\begin{aligned} \begin{aligned} x&= \xi \cosh \mu \cos \nu \;, \\ y&= \xi \sinh \mu \sin \nu \;, \end{aligned} \end{aligned}$$

(C6)

and their inverse

$$\begin{aligned} \begin{aligned} \mu&= {{\,\textrm{arccosh}\,}}\left[ \frac{\sqrt{(x+\xi )^2+y^2}+\sqrt{(x-\xi )^2+y^2}}{2\xi }\right] \; , \\ \nu&= \arccos \left[ \frac{\sqrt{(x+\xi )^2+y^2}-\sqrt{(x-\xi )^2+y^2}}{2\xi }\right] . \end{aligned} \end{aligned}$$

(C7)

The above equations represent a family of confocal ellipses where the foci are at $-\xi$ and $\xi$ on the x-axis. In (C6) and (C7), $\nu \in [0,2\pi )$ and each value of $\mu \geqslant 0$ corresponds to an ellipse with eccentricity $e=1/\cosh \upmu$ ($\upmu =0$ is the degenerate case). Therefore, each point (x, y) of the domain belongs to one and only one ellipse, whose properties are worked out as follows: a point (x, y) is given; $\upmu$ (and therefore e) is calculated with the first of (C7); the semi-major and semi-minor axes are then given by $a=\xi /e$ and $b=\xi \sqrt{1-e^2}/e$, respectively. Note that $\alpha$ is constant and is defined for each object being simulated.

Now, for each point (x, y) we can calculate the velocity field $\mathbf {u_n}(x,y)$ around the comet, in the form

$$\begin{aligned} \mathbf {u_n}(x,y)=|\mathbf {u_{n}}|~\frac{\nabla f(x,y)}{|\nabla f(x,y)|}, \end{aligned}$$

(C8)

where $f(x,y)=e^2x^2/\xi ^2+e^2y^2/[\xi ^2(1-e^2)]-1$. The velocity $\textbf{v}$ in (4) and (6) is treated similarly.

On the other hand, we must adapt (5) to elliptic comets, once it was originally intended for spherical sources. Here we propose a simple approach: we interpret r as the average value of the radial coordinate of the ellipse with center at the origin, that is,

$$\begin{aligned} r~:= \overline{r}(e)=\frac{a}{2\pi }\int ^{2\pi }_{0}\sqrt{\frac{(1-e^2)}{1-e^2\cos ^2(\theta )}}\ d\theta = \frac{a}{2\pi }I(e). \end{aligned}$$

(C9)

With this, for each point (x, y) we have an defined value of  $\overline{r}(e)$. Unfortunately, (C9) does not have a solution in terms of elementary functions for $e \ne 0$. However, by numerical integration, we can calculate I(e) for prescribed values of e and perform a curve fitting of the data to obtain an expression for $\overline{r}(e)$. Table 3 shows the selected values of e and their corresponding I(e).

Table 3 The integral I(e) of (C9) calculated for some values of e

The 4th degree polynomial

$$\begin{aligned} I(e)=-30.789 e^4+51.206 e^3-30.283 e^2+5.810 e+5.974 \end{aligned}$$

(C10)

is a good fit for the data of Table 3 in the interval $e \in [0,1)$, giving 0.22561 and 0.0150407 for the final sum of the squares of residuals and the reduced chi-squares, respectively¹.

Finally, we must bear in mind that, in our model, Eq. (C8) is not calculated inside the object, that is, it is only valid if the following condition is satisfied:

$$\begin{aligned} \sqrt{x^2+\frac{y^2}{1-e_\textrm{com}^2}}\;\geqslant \;a_\textrm{com}, \end{aligned}$$

(C11)

where $e_\textrm{com}$ and $a_\textrm{com}$ are, respectively, the eccentricity and the semi-major axis of the comet.

1.2 C.2 Bilobed Comet

We employ the polar equation $r(\theta )=p\rho +q\cos (2\theta )$, where the radial coordinate $r(\theta )$ is measured from the origin, p and q are constants and $\rho \in [1,\infty )$. This equation describes a family of concentric curves (one for each value of $\rho$) that intersect the x-axis and the y-axis at $(\pm \alpha ,0)$ and $(0,\pm \beta )$, where $\alpha =p\rho +q$ and $\beta =p\rho -q$. Note that the bilobed curves become more circular, that is, $\alpha \approx \beta \approx p\rho$ for $\rho \gg 1$. On the other hand, in our model the surface of the comet corresponds to $\rho =1$ by definition, such that the parameters a and b (which define the dimensions of the object) are $\alpha _\textrm{com}=p+q$ and $\beta _\textrm{com}=p-q$.

From the aforementioned polar equation, we can deduce the parametric equations

$$\begin{aligned} \begin{aligned} x&= [p\rho +q\cos (2\theta )]\cos \theta , \\ y&= [p\rho +q\cos (2\theta )]\sin \theta , \end{aligned} \end{aligned}$$

(C12)

and their inverses, in the form

$$\begin{aligned} \begin{aligned} \rho&= \frac{1}{p}\sqrt{x^2+y^2}-\frac{q}{p}\left( \frac{x^2-y^2}{x^2+y^2}\right) , \\ \theta&= \arctan \left( \frac{y}{x}\right) . \end{aligned} \end{aligned}$$

(C13)

The initial velocity field of the ejected plasma is calculated with (C8) and $f(r,\theta )=r-p\rho -q\cos (2\theta )$. With this and (C13), we obtain

$$\begin{aligned} \textbf{u}(x,y)=\frac{u_0}{|\nabla f(x,y)|}\left\{ \left[ \frac{x}{r}-\frac{2(\alpha -\beta )xy^2}{r^4}\right] \hat{\textbf{x}}+\left[ \frac{y}{r}+\frac{2(\alpha -\beta )x^2y}{r^4}\right] \hat{\textbf{y}}\right\} , \end{aligned}$$

(C14)

where

$$\begin{aligned} |\nabla f(x,y)|=\left\{ 1+\left[ \frac{2(\alpha -\beta )xy}{r^3}\right] ^2\right\} ^{1/2}, \end{aligned}$$

(C15)

and $r=(x^2+y^2)^{1/2}$. Note that Eqs. (C14) and (C15) are written in the more convenient parameters $\alpha$ and $\beta$.

As in the previous case, $\mathbf {u_n}(x,y)$ is not calculated inside the body. With this and bearing in mind that $\alpha _\textrm{com}$ and $\beta _\textrm{com}$ can be related by $\alpha _\textrm{com}=c\beta _\textrm{com}$ where c is a constant, the condition of validity of Eq. (C14) is

$$\begin{aligned} \frac{c(x^2+y^2)^{3/2}}{cx^2+y^2}\;\geqslant \;\alpha _\textrm{com}. \end{aligned}$$

(C16)

Regarding the calculation of (5), we again consider r as the average value of $r(\theta )$ over $[0,2\pi )$, that is,

$$\begin{aligned} r :=\overline{r}=\frac{1}{2\pi }\int ^{2\pi }_{0}[p\rho +q\cos (2\theta )]d\theta =\rho \frac{\alpha +\beta }{2}. \end{aligned}$$

(C17)

As an example, the left panel of Fig. 14 shows the vector field of $\mathbf {u_n}(x,y)$ for an ellipse with $e_\textrm{com}=0.9$ and $a_\textrm{com}=0.25$, in normalized units. In the right panel we have the bilobed form, with $\alpha _\textrm{com}=0.25$ and $\beta _\textrm{com}=\alpha _\textrm{com}/3$.

Fig. 14

Left: vector field of $\mathbf {u_n}(x,y)$ for the elliptic case, with $\mathbf {u_n}(x,y)$ for $e_\textrm{com}=0.9$. Right: the same for the bilobular case with $\alpha _\textrm{com}=0.25$ and $\beta _\textrm{com}=\alpha _\textrm{com}/3$