Anelastic convection-driven dynamo benchmarks: A finite element model

Icarus 218 (2012) 345–347

Contents lists available at SciVerse ScienceDirect

Icarus journal homepage: www.elsevier.com/locate/icarus

Note

Anelastic convection-driven dynamo benchmarks: A finite element model Xiaoya Zhan a, Gerald Schubert a,⇑, Keke Zhang b a b

Department of Earth and Space Sciences, University of California, Los Angeles, CA 90095-1567, USA Department of Mathematical Sciences, University of Exeter, Exeter EX4 4QF, UK

a r t i c l e

i n f o

Article history: Received 28 October 2011 Revised 14 December 2011 Accepted 17 December 2011 Available online 29 December 2011

a b s t r a c t We present a fully non-linear, anelastic convection and dynamo model in a rotating spherical shell that is based on a finite element method. It is shown that both the hydrodynamic steadily drifting convection and convection-driven steadily drifting dynamo are in satisfactory agreement with the anelastic convection and dynamo benchmark based on the pseudo-spectral method (Jones, C.A., Boronski, P., Brun, A.S., Glatzmaier, G.A., Gastine, T., Miesch, M.S., Wicht, J. [2011]. Icarus, 216, 120–135).

Ó 2012 Elsevier Inc. All rights reserved.

Keywords: Magnetic fields Jovian planets Interiors

1. Introduction Over the past several decades, the anelastic approximation has been widely employed in numerical models of thermal convection and convection-driven dynamos in the interiors of giant planets and stars (e.g., Gilman and Glatzmaier, 1981; Glatzmaier, 1984, 1985; Miesch et al., 2000; Brun et al., 2004; Jones and Kuzanyan, 2009). An anelastic, convection-driven dynamo benchmark for a rotating spherical shell, following a highly successful benchmark for the Boussinesq problem (Christensen et al., 2001), is now available (Jones et al., 2011). The benchmark is based on four different spherical dynamo models – the Leeds code, the MAGIC code, the ASH code and the Glatzmaier code – all of which have employed the global pseudo-spectral method using spherical harmonic expansions. There are advantages and disadvantages to using the spherical spectral approximation in the era of modern massively parallel computers. An essential feature in the pseudo-spectral method for the anelastic approximation is to use the poloidal– toroidal decomposition by writing

q u ¼ r  r  rq P þ r  rq T ;  is the basic state density, and P and T represent, where r is the position vector, q respectively, the poloidal and toroidal component of the velocity u. A computational advantage of using the poloidal–toroidal decomposition is that both the anelastic continuity equation and the magnetic flux conservation equation are automatically satisfied. Additionally, the magnetic field boundary condition at the interface between the electrically conducting fluid core and the electrically insulating exterior of a planet can be readily implemented in an analytical way. Moreover, the spherical pseudo-spectral method is not only highly accurate but also easy for coding, because spherical harmonics are directly associated with some differential operators in the governing equations. However, the spherical harmonic expansion leads to a global integration that requires intensive global communication, making it less efficient on modern massively parallel computers. It is also well-known that the Legendre transform used in evaluating the surface integral limits the effectiveness of the pseudo-spectral method. Moreover, the pseudo-spectral method cannot be easily used in dealing with non-spherical geometry, the effect of which plays an essential role in dynamo action driven by planetary precession, libration or tides.

⇑ Corresponding author. E-mail address: [email protected] (G. Schubert). 0019-1035/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. doi:10.1016/j.icarus.2011.12.016

It is thus desirable to develop an independent, anelastic convection-driven dynamo code using local non-spectral methods. This Note reports the result of an anelastic convection-driven dynamo code that is based on a finite element method using a fully three-dimensional tetrahedralization that verifies the pseudo-spectral benchmark (Jones et al., 2011). There are also advantages and disadvantages in a finite element approximation. Finite element methods are local and without the need of intensive global communication and, hence, can take full advantage of modern massively parallel computers. Important theoretical aspects of our finite element approximation, such as the existence, uniqueness and stability of the numerical scheme, have been established for the dynamo problem (Chan et al., 2006). Moreover, a spherical finite element code can be readily modified for non-spherical geometries such as spheroids and triaxial ellipsoids. A major disadvantage is the conflict between the local nature of finite element methods and the global nature of magnetic field boundary conditions and, consequently, the potential field, which can be treated analytically in the pseudo-spectral method, has to be treated numerically entailing extra computational time. 2. Model and governing equations As considered by Jones et al. (2011), we consider a perfect gas confined within a spherical shell of an inner radius ri and an outer radius ro with the radius ratio denoted by v = ri/ro. By assuming that gravity is proportional to 1/r2, a polytropic solu ¼ qc fn , T ¼ T c f, P ¼ P c fnþ1 ; tion for the basic state can be written as q f ¼ c0 þ c1 d=r. Here n is the polytropic index, qc, Pc and Tc are the reference values of density, pressure and temperature at mid-depth of the shell, respectively, and constants c0 and c1 are given by c0 ¼ ð2f0  v  1Þ=ð1  vÞ, c1 ¼ ð1 þ vÞð1  f0 Þ=  ðr i Þ=q  ðr o ÞÞ representing the numð1  vÞ2 , fo ¼ ð1 þ vÞ=ðveNq =n þ 1Þ, with N q ¼ lnðq ber of scale heights of density and fo denoting the value of f at the outer boundary. The details of the physical model are discussed in (Jones et al., 2011). The system rotates uniformly with an angular velocity X, and has constant magnetic diffusivity k, magnetic permeability l, kinematic viscosity m and entropy diffusivity j. Convection takes place when the entropy contrast DS between the inner and outer boundaries becomes sufficiently large (Jones and Kuzanyan, 2009). With length scaled by d, time scaled by d2/k, qc as the unit of density, Tc as the unit of temperature, DS as the unit of entropy S and magnetic field B scaled by pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Xqc lk, the governing equations for the magnetohydrodynamic process can be written as

346

fn



X. Zhan et al. / Icarus 218 (2012) 345–347



@u Pm n þ ðu  rÞu þ 2 f ^z  u  Pmr  Fm @t E

¼ fn rP þ

Pm2 Ra n S Pm f 2 ^r þ ðr  BÞ  B; Pr r E

r  ðfn uÞ ¼ 0; fnþ1



ð1Þ

Table 1 Comparison between solution characteristics obtained from the finite element (FE) method and the pseudo-spectral (PS) method of Jones et al. (2011). The steadily drifting convection is denoted by subscript hydro while the steadily drifting convection-driven dynamo is denoted by subscript dynamo.

ð2Þ

 @S Pm c Pr c Pr r  fnþ1 rS þ 1 fn Q m þ 1 þ ðu  rÞSÞ ¼ ðr  BÞ2 ; @t Pr PmRa RaPmE

@B ¼ r  ðu  BÞ  r  r  B; @t r  B ¼ 0;

ð3Þ

ð4Þ ð5Þ

Ekin Emag L u/ S Bh

x

PShydro

FEhydro

PSdynamo

FEdynamo

81.858 – 4.1989 0.8618 0.9330 – 17.6427

81.367(0.6%) – 4.1901(0.21%) 0.8549(0.8%) 0.9314(0.17%) – 17.31(1.9%)

419,407 320,185 11.5030 91.780 0.78646 0.03395 611.864

424,439(1.2%) 341,958(6.8%) 11.413(0.78%) 88.751(3.3%) 0.78025(0.79%) 0.03443(1.4%) 592.892(3.1%)

where

 @ui @uj 2 @uj  fn þ and @xj @xi 3 @xj     1 1 @ui @uj with eij ¼ : þ Q m ¼ 2 eij eij  ðr  uÞ2 3 2 @xj @xi Fm ¼ fn



ð6Þ

The problem is primarily characterized by four dimensionless parameters, the Rayleigh number Ra, the Ekman number E, the magnetic Prandtl number Pm, and the Prandtl number Pr, defined as

Ra ¼

GMdDS ; mjcp

E¼

m X d2

;

m

Pm ¼ ; k

Pr ¼

m : j

We assume impenetrable, stress-free and constant entropy conditions at the bounding surfaces of the spherical shell. In comparison to the benchmark dynamo model (Jones et al., 2011), the magnetic diffusivity in the exterior of the gas shell is assumed to be sufficiently large that the electrically insulating condition is approximately satisfied (Zhang et al., 2003; Chan et al., 2007). More precisely, the dynamo equation with u = 0 for the exterior is solved numerically and the exterior magnetic diffusivity is taken to be a factor 10,000 larger than that of the fluid. This represents a major difference between our finite element dynamo and the pseudo-spectral benchmark model.

3. Finite element convection and dynamo model We first construct three different tetrahedra finite element meshes for the inner sphere, the fluid shell and the outer exterior. The three-dimensional tetrahedralization of the whole system produces a uniform mesh distribution on any spherical surface without the pole problem and a small number of nodes with a nearly uniform mesh distribution in the neighborhood of the center r = 0 without the origin problem. For the anelastic simulation reported here, we have used 34,734,080 tetrahedral elements with a total of 417,554,793 unknowns. A semi-implicit time stepping scheme is adopted for the time advancement while a second-order extrapolation is applied to all the non-linear terms. A Galerkin weighted residual approach is used in the finite element formulation. This makes use of mixed finite elements in which piecewise quadratic polynomials are employed to approximate the velocity and magnetic field while piecewise linear polynomials are used for approximating the pressure (Hood and Taylor, 1974). An important requirement for any dynamo benchmark is its simplicity, reliability and reproducibility (Christensen et al., 2001), all of which have been successfully achieved by the pseudo-spectral benchmark (Jones et al., 2011). This benchmark adopts a set of parameters for which fully non-linear convection-driven dynamos are steady in a drifting coordinate frame. In other words, the main characteristics of a time-dependent solution such as the kinetic and magnetic energies are constants that can be easily validated by independent non-spectral methods. We will focus on the two steady benchmark solutions, a steadily drifting hydrodynamic convection and a steadily drifting convection-driven dynamo, because an accurate comparison can be easily made. According to Jones et al. (2011), a steady hydrodynamic convection is found by using the non-dimensional parameters Ra = 351,806, E = 103, Pr = 1, Nq = 5, n = 2, v = 0.35. In simulating non-magnetic convection, we set B = 0 and Pm = Pr in the governing equations. There are five primary quantities that can be used to characterize a steadily drifting convection: kinetic energy Ekin, luminosity L, the drifting frequency x, the values of u/ and S at a particular point where ur = 0. It is significant to note that these five quantities can be identified from any convection solution independent of the numerical method employed. Table 1 shows the five characteristic values for the steadily drifting convection obtained from our finite element model along with those obtained by Jones et al. (2011). The agreement between the finite element and pseudo-spectral models, considering two fundamentally different codes, is surprisingly satisfactory. For example, Jones et al. (2011) yields the luminosity L ¼ 0:9330 while our model gives L ¼ 0:9314; Jones et al. (2011) produce the kinetic energy Ekin = 81.86 while ours is Ekin = 81.37 with less than 1% difference. There are some differences between the profiles of the two different convection solutions that are displayed in Fig. 1. Similarly, a steadily drifting convection-driven dynamo, according to Jones et al. (2011), can be found by setting Ra = 80,000, E = 2  103, Pr = 1, Pm = 50, Nq = 5,

n = 2, v = 0.35. Table 1 also shows seven characteristic values for the steady dynamo obtained from our finite element code, along with those given by Jones et al. (2011). It can be seen that while the agreement for the hydrodynamic quantities such as the kinetic energy is satisfactory, there exists about a 7% difference between the values of magnetic energy. This discrepancy is likely attributable to the approximate magnetic boundary condition used in our finite element code which allows a weak current to enter the poorly conducting exteriors of the fluid shell. Fig. 2 illustrates the profiles of two different dynamo solutions without showing any noticeable differences between them. 4. Summary and remarks We have developed a fully non-linear, anelastic code for thermal convection and convection-driven dynamo action in rotating spherical shells based on a finite element method. There is satisfactory agreement between our finite element model and the anelastic convection and dynamo benchmark of Jones et al. (2011) that is based on the pseudo-spectral method making use of a spherical harmonic expansion with the toroidal–poloidal decomposition. We focused on two steady cases because the non-uniqueness and uncertainty of a chaotic solution make it difficult to compare fundamentally different codes. It should be pointed out that while some quantities such as the kinetic energy are well defined and easily reproducible independent of numerical codes, some solution values like the drift frequency x cannot be explicitly defined in our finite element solutions and, hence, is harder to calculate accurately. In consequence, its value cannot be accurately measured and determined from our time-dependent

Fig. 1. A comparison of steadily drifting convection between our finite element model and the pseudo-spectral model of Jones et al. (2011). Contours of ur in the equatorial plane: (a) from the PS model and (c) from the FE model. Contours of S on a spherical surface r = 1.4386 (b) from the PS model and (d) from the FE model. The models are computed for the same parameter values.

X. Zhan et al. / Icarus 218 (2012) 345–347

347

Acknowledgments We thank Professors C.A. Jones and G.A. Glatzmaier for helpful discussions. G.S. is supported by the National Science Foundation under grant NSF AST-0909206 and KZ is supported by UK NERC, STFC and Leverhulme grants. References

Fig. 2. A comparison of steadily drifting convection-driven dynamo between our finite element model and the pseudo-spectral model of Jones et al. (2011). Contours of ur in the equatorial plane: (a) from the PS model and (c) from the FE model. Contours of Br on a spherical surface r = 1.3303 (b) from the PS model and (d) from the FE model. The models are computed for the same parameter values.

non-linear solution. Finally, it is also worth noting that in contrast to a pseudo-spectral dynamo model, the existence of an electrically conducting heterogeneous exterior fluid shell within a planet, if needed, can be readily incorporated into our finite element dynamo model.

Brun, A.S., Miesch, M.S., Toomre, J., 2004. Global-scale turbulent convection and magnetic dynamo action in the solar envelope. Astrophys. J. 614, 1073–1098. doi:10.1086/423835. Chan, K.H., Zhang, K., Zou, J., 2006. Spherical interface dynamos: Mathematical theory, finite element approximation, and application. SIAM J. Numer. Anal. 44, 1877–1902. Chan, K.H., Zhang, K., Li, L., Liao, X., 2007. A new generation of convection-driven spherical dynamos using EBE finite element method. Phys. Earth Planet. Inter. 163, 251–265. doi:10.1016/j.pepi.2007.04.017. Christensen, U.R. et al., 2001. A numerical dynamo benchmark. Phys. Earth Planet. Inter. 128, 25–34. doi:10.1016/S0031-9201(01)00275-8. Gilman, P.A., Glatzmaier, G.A., 1981. Compressible convection in a rotating spherical shell. I. Anelastic equations. Astrophys. J. Suppl. Ser. 45, 335–349. doi:10.1086/ 190714. Glatzmaier, G.A., 1984. Numerical simulations of stellar convective dynamos. I. The model and method. J. Comput. Phys. 55, 461–484. doi:10.1016/0021– 9991(84)90033–0. Glatzmaier, G.A., 1985. Numerical simulations of stellar convective dynamos. II. Field propagation in the convection zone. Astrophys. J. 291, 300–307. doi:10.1086/163069. Hood, P., Taylor, C., 1974. Navier–Stokes equations using mixed-interpolation. In: Oden, J.T., Zienkiewicz, O.C., Gallagher, R.H., Taylor, C. (Eds.), Finite Element Methods in Flow Problems. UAH Press, Huntsville, AL, pp. 57–66. Jones, C., Kuzanyan, K.M., 2009. Compressible convection in the deep atmospheres of giant planets. Icarus 204, 227–238. doi:10.1016/j.icarus.2009.05.022. Jones, C.A., Boronski, P., Brun, A.S., Glatzmaier, G.A., Gastine, T., Miesch, M.S., Wicht, J., 2011. Anelastic convection-driven dynamo benchmarks. Icarus 216, 120–135. doi:10.1016/j.icarus.2011.08.014. Miesch, M.S., Elliott, J.R., Toomre, J., Clune, T.L., Glatzmaier, G.A., Gilman, P.A., 2000. Three-dimensional spherical simulations of solar convection. I. Differential rotation and pattern evolution achieved with laminar and turbulent states. Astrophys. J. 532, 593–615. doi:10.1086/308555. Zhang, K., Chan, K.H., Zou, J., Liao, X., Schubert, G., 2003. A three-dimensional spherical nonlinear interface dynamo. Astrophys. J. 596, 663–679. doi:10.1086/ 377600.