# Solar Convection Simulations

using a B-spline method

## 1 Motivation and objectives

Convection plays a key role in energy transport and global circulation in the outer layers of the Sun, in generation of solar magnetic fields and in many phenomena associated with solar activity and variability. Observations of the solar surface reveal structures that have been classified as granules, mesogranules, and supergranules. The nomenclature reflects organization at three spatial scales ranging from about 1 Mm to 30 Mm.

Numerical simulations of the near surface region of the Sun Stein & Nordlund (2000) capture structures on the granular scale, but have not been able to detect organization at large scales. The physical mechanism of supergranulation is presently unknown. It has been suggested that supergranulation corresponds to a large convective cells which develop due to enhanced convective instability in the HeII ionization layer. This layer lies deeper than Stein & Nordlund have investigated so far.

We know that as we go deeper from the surface of the sun, the turbulence structures become large, and the Mach number decreases. It is then advantageous to be able to change the spatial resolution in all three coordinates as a function of depth. In addition, it is numerically advantageous to use the anelastic approximation Ogura & Phillips (1962); Gough (1969); Gilman & Glatzmaier (1981) to the compressible Navier-Stokes equations for deep domains. Using B-splines in one coordinate direction and Fourier methods in the other two coordinate directions, Kravchenko et al. (1996) and Loulou et al. (1997) have developed B-spline–spectral numerical schemes that enable changing the resolution as a function of height. Kravchenko et al. (1996) applied the scheme to simulate the fully developed turbulent channel flow by resolving the near-wall region and relaxing the resolution in the core region. Loulou et al. (1997) simulated a fully developed turbulent pipe flow. They applied the scheme to remove the centerline singularity in cylindrical coordinates.

## 2 Numerical method

Under this effort, a new code to simulate a rectilinear section of the solar convection zone is being implemented. As a first step, the numerical method is implemented for a simple Boussinesq fluid Oberbeck (1879); Boussinesq (1903) and is presented in this paper for this case only. Later, once the new code is fully functional and tested, equations based on the anelastic approximation will be used.

### 2.1 Basic equations

Using the temperature difference across the fluid layer, , the layer thickness and the thermal diffusion time as units of temperature, length and time, the dimensionless Boussinesq equations read:

(1) |

(2) |

(3) |

, and denote the velocity of the fluid and the deviations of the temperature and pressure from their static profiles, i.e. the temperature and pressure profiles that would be observed in the case of no motion. The equations contain two dimensionless parameters: the Prandtl number and the Rayleigh number , which are defined as

(4) |

with , , and being the coefficients of thermal conductivity, kinematic viscosity and thermal expansion, and the accelleration due to gravity, which points in negative direction. The Boussinesq approximation assumes all these properties to be constant throughout the layer.

### 2.2 Velocity decomposition

A poloidal/toroidal decomposition is used for the velocity to automatically fulfill the continuity equation . Even though the divergence of the velocity does not vanish in the anelatic approximation, the equations that we are ultimately going to use, such a representation can still be used. Just not for the velocity itself, but for the quantity , where denotes the reference density.

In the Boussinesq case, the representation is in the form of

(5) |

In this way the poloidal and toroidal scalars and are periodic in and , and represents the horizontal average of the velocity, i.e

(6) |

The component of must be zero to satisfy the continuity equation. is thus a toroidal flow. However, the corresponding toroidal scalar varies linearly in and and is unbounded. Obviously, we require to be periodic and to remain finite in the numerical representation. Therefore needs to be included in the decomposition of Schmitt & von Wahl (1992).

### 2.3 Spectral method for the horizontal directions

Periodic boundary conditions are used in the horizontal directions to reduce the influence of the horizontal boundaries on the flow. A natural choice is the use of Fourier modes in those directions, i.e.

(7) |

with being the horizontal wave vector. Multiplying the equations of motion with weight functions and integrating over the entire horizontal plane then yields equations for the and dependent Fourier coefficients:

(9) | |||||

(11) | |||||

(12) |

All non-linear terms are contained in the right-hand sides . The equation for is very similar:

(13) |

If fact, we can discuss the equation for , and together, since they all can be written in the form

(14) |

with some parameters and .

### 2.4 B-spline method for the vertical direction

yy1 \psfragyy2 \psfragyy3 \psfragyy4 \psfragyy5 \psfragyy6 \psfraga \psfragb \psfragc \psfragd \psfrage \psfragx

y

Spatial discretization in the vertical direction is done by a B-spline method Kravchenko & Moin (1998), i.e. the unknowns, and , are expanded in terms of -order piecewise polynomials called basis splines, or B-splines. With a given set of knot points that divide the domain into subintervals, one can construct of these spline functions using the recursive expression

(15) |

Here, denotes the -th B-spline of order . The -order splines are given by

(16) |

As it turns out, additional points at each side of the domain ( and ) are needed in the construction of B-splines near the boundaries. These virtual points are chosen to be equal to the boundary points and , i.e. and .

Second-order B-spline functions for two different sets of knot points are shown in figure 2.4 to illustrate how these functions look like. As can be seen from the figure, B-splines have compact support and are non-negative:

(17) |

and can be tuned to the specific needs at hand by the choice of knot points. In areas where small structures need to be resolve, e.g. in boundary layers, one can choose closely spaced knot points, while in areas where such high resolution is not required, the grid can be coarse to save computational costs.

Incorporating a B-spline expansion of the unknowns, i.e.

(18) |

into the governing equations (9) and (14), multiplying with weight functions and integrating over the domain yields equations for the expansion coefficients :

(19) |

(20) |

Three matrices, , and , occur in these equations, which contain integrals of the product of two B-splines and the product of a B-spline with a derivative of a B-spline:

(21) |

It follows from the property (17) that all of these matrices are of band-diagonal form. Specifically, they have only subdiagonals and superdiagonals which contain non-zero entries.

### 2.5 Time advancement

A mixed implicit/explicit method using a Crank-Nicolson scheme for the diffusive terms and a second-order Adams-Bashforth scheme for the advection and buoyancy terms is used for the time advancement.

In the end, one has to solve a matrix equation for each Fourier mode of each field (poloidal, toroidal, mean flow and temperature) in the form

(22) |

with representing the vector of B-spline coefficients at the new time step , i.e. ( again stands for either , , or ), and the matrix is composed of the previously introduced matrices , and . is therefore band-diagonal and a fast method for solving band-diagonal systems can be used.

### 2.6 Imposing boundary conditions

So far, no boundary conditions have been imposed in the vertical direction.
Since only one B-spline, and the derivative of only two B-splines are non-zero at the boundary,
it is rather easy to implement various boundary conditions.
If for example one needs to prescribe the value of at the boundary , e.g.
, one simply has to set the first B-spline coefficient to that value, i.e.
.
The original equation for has to be removed from the linear system (22),
and ’s contribution has to be deducted from the right hand, i.e.
equation (22) is modified in the following way:

\begin{picture}(13.0,7.0)\par\put(0.0,5.2){$\left(\begin{array}[]{cccccc}\star%
&\star&&&&\\
\star&\star&\star&&&\\
&\star&\star&\star&&\\
&&\ddots&\ddots&\ddots&\\
&&&\star&\star&\star\\
&&&&\star&\star\\
\end{array}\right)\left(\begin{array}[]{c}\alpha_{f,1}\\
\alpha_{f,2}\\
\\
\vdots\\
\\
\alpha_{f,n+k}\\
\end{array}\right)=\left(\begin{array}[]{c}\beta_{f,1}\\
\beta_{f,2}\\
\\
\vdots\\
\\
\beta_{f,n+k}\\
\end{array}\right)$}
\par\put(2.17,1.3){$\left(\begin{array}[]{ccccc}\star&\star&&&\\
\star&\star&\star&&\\
&\ddots&\ddots&\ddots&\\
&&\star&\star&\star\\
&&&\star&\star\\
\end{array}\right)\left(\begin{array}[]{c}\alpha_{f,2}\\
\\
\vdots\\
\\
\alpha_{f,n+k}\\
\end{array}\right)=\left(\begin{array}[]{c}\beta_{f,2}\\
\\
\vdots\\
\\
\beta_{f,n+k}\\
\end{array}\right)-\alpha_{f,1}\left(\begin{array}[]{c}\star\\
\\
\\
\\
\\
\end{array}\right)$}
\par\par\put(0.38,6.25){\line(0,-1){2.35}}
\put(0.38,6.25){\line(1,0){0.39}}
\put(0.77,6.25){\line(0,-1){2.35}}
\put(0.38,3.9){\line(1,0){0.39}}
\par\put(0.88,6.25){\line(0,-1){2.35}}
\put(0.88,6.25){\line(1,0){3.22}}
\put(4.1,6.25){\line(0,-1){2.35}}
\put(0.88,3.9){\line(1,0){3.22}}
\par\put(4.96,6.25){\line(0,-1){2.35}}
\put(4.96,6.25){\line(1,0){1.1}}
\put(6.06,6.25){\line(0,-1){2.35}}
\put(4.96,3.9){\line(1,0){1.1}}
\par\put(7.4,6.25){\line(0,-1){2.35}}
\put(7.4,6.25){\line(1,0){1.1}}
\put(8.5,6.25){\line(0,-1){2.35}}
\put(7.4,3.9){\line(1,0){1.1}}
\par\par\put(6.62,2.55){\line(0,-1){2.35}}
\put(6.62,2.55){\line(1,0){1.1}}
\put(7.72,2.55){\line(0,-1){2.35}}
\put(6.62,0.2){\line(1,0){1.1}}
\par\put(9.08,2.55){\line(0,-1){2.35}}
\put(9.08,2.55){\line(1,0){1.1}}
\put(10.18,2.55){\line(0,-1){2.35}}
\put(9.08,0.2){\line(1,0){1.1}}
\par\put(2.55,2.55){\line(0,-1){2.35}}
\put(2.55,2.55){\line(1,0){3.22}}
\put(5.77,2.55){\line(0,-1){2.35}}
\put(2.55,0.2){\line(1,0){3.22}}
\par\put(12.08,2.55){\line(0,-1){2.35}}
\put(12.08,2.55){\line(1,0){0.39}}
\put(12.47,2.55){\line(0,-1){2.35}}
\put(12.08,0.2){\line(1,0){0.39}}
\par\par\par\put(2.69,3.75){\line(2,-1){0.9}}
\put(3.79,3.2){\line(2,-1){0.9}}
\put(4.79,2.7){\vector(2,-1){0.0}}
\par\put(5.6,3.75){\line(2,-1){0.9}}
\put(6.7,3.2){\line(2,-1){0.9}}
\put(7.7,2.7){\vector(2,-1){0.0}}
\par\put(7.77,3.75){\line(2,-1){0.9}}
\put(8.87,3.2){\line(2,-1){0.9}}
\put(9.87,2.7){\vector(2,-1){0.0}}
\par\put(0.58,3.75){\line(0,-1){0.5}}
\put(0.58,3.25){\line(1,0){11.7}}
\put(12.28,3.25){\line(0,-1){0.48}}
\put(12.28,2.7){\vector(0,-1){0.0}}
\end{picture}